source: CPL/oasis3-mct/branches/OASIS3-MCT_5.0_branch/examples/regrid_environment/src/gradient_bicubic.F90 @ 6331

Last change on this file since 6331 was 6331, checked in by aclsce, 17 months ago

Moved oasis-mct_5.0 in oasis3-mct/branches directory.
File size: 10.3 KB
Line 
1subroutine gradient_bicubic(NX1, NY1, ibeg, jbeg, iloc, jloc, &
2                            cl_grd_src, id_per, cd_per, w_unit, &
3                            local_gradient_i, local_gradient_j, local_gradient_ij, file_debug)
4!
5!**** *gradient_bicubic*  - calculate gradients for bicubic remapping
6!
7!     Purpose:
8!     -------
9!     Calculation of gradients for bicubic interpolation. In contrast to
10!     the gradients of conservative remapping, these gradients are   
11!     calculated with respect to grid rows and grid lines.
12!
13!**   Interface:
14!     ---------
15!       *CALL*  *gradient_bicubic*(NX1, NY1, ibeg, jbeg, iloc, jloc,
16!                                  cl_grd_src, id_per, cd_per, w_unit,
17!                                  local_gradient_i, local_gradient_j, local_gradient_ij, file_debug)
18!
19!     Input:
20!     -----
21!          NX1            : grid global dimension in x-direction (integer)
22!          NY1            : grid global dimension in y-direction (integer)
23!          ibeg           : start of local domain in global domain in x-direction
24!          jbeg           : start of local domain in global domain in y-direction
25!          iloc           : grid local dimension in x-direction (integer)
26!          jloc           : grid local dimension in y-direction (integer)
27!          cl_grd_src     : grid acronym
28!          id_per         : number of overlapping points for source grid
29!          cd_per         : grip periodicity type
30!          w_unit         : log file unit
31!          file_debug     : logical for activating debug outputs
32!
33!     Output:
34!     ------
35!          local_gradient_i     : gradient in i-direction (real 2D)
36!          local_gradient_j     : gradient in j-direction (real 2D)
37!          local_gradient_ij    : gradient in ij-direction (real 2D)
38!
39! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
40!
41USE read_all_data
42USE function_ana
43!
44IMPLICIT NONE
45     
46INTEGER, PARAMETER :: wp = SELECTED_REAL_KIND(12,307) ! double
47!-----------------------------------------------------------------------
48!     INTENT(IN)
49!-----------------------------------------------------------------------
50INTEGER, INTENT(IN) :: &
51         NX1, NY1,  &          ! source grid global dimensions
52         ibeg, jbeg, &         ! source grid local start
53         iloc, jloc            ! source grid local dimensions
54
55CHARACTER(len=4), INTENT(IN) ::  &
56         cl_grd_src            ! grid acronym
57
58INTEGER, INTENT(IN) :: &
59         id_per, &             ! nbr of overlapping grid points
60         w_unit                ! log file
61
62CHARACTER(len=8), INTENT(IN) ::  &
63         cd_per                ! grip periodicity type
64
65LOGICAL, INTENT(IN)      :: file_debug
66!
67!-----------------------------------------------------------------------
68!     INTENT(OUT)
69!-----------------------------------------------------------------------
70REAL (kind=wp), DIMENSION(iloc, jloc), INTENT(OUT) :: &
71          local_gradient_i,   &       ! gradient in i-direction (real 2D)
72          local_gradient_j,   &       ! gradient in j-direction (real 2D)
73          local_gradient_ij           ! gradient in ij-direction (real 2D)
74
75!-----------------------------------------------------------------------
76!     LOCAL VARIABLES
77!-----------------------------------------------------------------------
78INTEGER ::  &
79          i, j,    &            ! looping indicees
80          ip1, jp1, im1, jm1, iend, jend
81     
82REAL (kind=wp) ::  &
83          di, dj,        &         ! factor depending on grid cell distance
84          gradient_ij1,  &         ! gradient needed to calculate gradient_ij
85          gradient_ij2             ! gradient needed to calculate gradient_ij
86
87REAL (kind=wp), DIMENSION(:,:), POINTER :: &
88          src_lon,   &          ! source grid longitudes [radiants]
89          src_lat,   &          ! source grid latitudes [radiants]
90          src_array, &          ! analytical field
91          gradient_i, &  ! global gradient in i-direction (real 2D)
92          gradient_j, &  ! global gradient in j-direction (real 2D)
93          gradient_ij    ! global gradient in ij-direction (real 2D)
94
95INTEGER, DIMENSION(:,:), POINTER :: &
96         sou_mask             ! source grid mask
97
98INTEGER, PARAMETER ::  il_maskval= 1 ! in our grids sea_value = 0 and land_value = 1
99
100!----------------------------------------------------------------------
101!
102!     Read global grid and global mask
103!     --------------------------------
104      ALLOCATE(src_lon(NX1, NY1))
105      ALLOCATE(src_lat(NX1, NY1))
106      CALL read_grid(NX1, NY1, 1, 1, NX1, NY1, cl_grd_src, w_unit, src_lon, src_lat, file_debug)
107!
108      ALLOCATE(sou_mask(NX1, NY1))
109      CALL read_mask(NX1, NY1, 1, 1, NX1, NY1, cl_grd_src, w_unit, sou_mask, file_debug)
110
111!     Global field from analytical function
112!     -------------------------------------
113      ALLOCATE(src_array(NX1, NY1))
114#ifdef FANA1
115      CALL function_ana1(NX1, NY1, src_lon, src_lat, src_array)
116#elif defined FANA2
117      CALL function_ana2(NX1, NY1, src_lon, src_lat, src_array)
118#elif defined FANA3
119      CALL function_ana3(NX1, NY1, src_lon, src_lat, src_array)
120#endif
121
122!     Global gradient allocation
123!     --------------------------
124      ALLOCATE(gradient_i(NX1, NY1))
125      ALLOCATE(gradient_j(NX1, NY1))
126      ALLOCATE(gradient_ij(NX1, NY1))
127
128!     Initialization
129!     --------------
130      gradient_i  = 0.
131      gradient_j  = 0. 
132      gradient_ij = 0. 
133
134!     calculate gradients
135!     -------------------
136      DO i = 1, NX1
137         DO j = 1, NY1
138                   
139            IF (sou_mask (i,j) /= il_maskval) THEN
140
141               di = 0.5
142               dj = 0.5
143
144               ip1 = i + 1
145               im1 = i - 1
146               IF (i == NX1) THEN
147                   IF (cd_per == 'P') ip1 = 1 + id_per ! the 0-meridian
148                   IF (cd_per == 'R') ip1 = NX1
149               ENDIF
150               IF (i == 1 )  THEN
151                   IF (cd_per == 'P') im1 = NX1 - id_per
152                   IF (cd_per == 'R') im1 = 1
153               ENDIF
154               jp1 = j + 1
155               jm1 = j - 1
156               IF (j == NY1) THEN ! treatment of the last..
157                  jp1 = NY1 
158                  dj = 1.
159               ENDIF   
160               IF (j == 1 ) THEN  ! .. and the first grid-row
161                  jm1 = 1
162                  dj = 1.
163               ENDIF
164
165
166!              gradient i
167!              ----------
168               IF (sou_mask(ip1,j) /= il_maskval .OR. &
169                   sou_mask(im1,j) /= il_maskval) THEN
170                  IF (sou_mask(ip1,j) == il_maskval) THEN
171                     ip1 = i
172                     di = 1.
173                  ELSE IF (sou_mask(im1,j) == il_maskval) THEN
174                     im1 = i
175                     di = 1.
176                  ENDIF
177                  gradient_i(i,j) = di * (src_array(ip1,j) - src_array(im1,j))
178               ENDIF
179
180!              gradient j
181!              ----------
182               IF (sou_mask(i,jp1) /= il_maskval .OR. &
183                   sou_mask(i,jm1) /= il_maskval) THEN
184                  IF (sou_mask(i,jp1) == il_maskval) THEN
185                     jp1 = j
186                     dj = 1.
187                  ELSE IF (sou_mask(i,jm1) == il_maskval) THEN
188                     jm1 = j
189                     dj = 1.
190                  ENDIF
191                  gradient_j(i,j) = dj * (src_array(i,jp1) - src_array(i,jm1))
192               ENDIF
193!
194!              gradient ij
195!              -----------
196               di = 0.5
197               dj = 0.5
198               ip1 = i + 1
199               im1 = i - 1
200               IF (i == NX1) THEN
201                   IF (cd_per == 'P') ip1 = 1 + id_per ! the 0-meridian
202                   IF (cd_per == 'R') ip1 = NX1
203               ENDIF
204               IF (i == 1 )  THEN
205                   IF (cd_per == 'P')  im1 = NX1 - id_per
206                   IF (cd_per == 'R')  im1 = 1
207               ENDIF
208               jp1 = j + 1
209               jm1 = j - 1
210               IF (j == NY1) THEN ! treatment of the last..
211                  jp1 = NY1 
212                  dj = 1.
213               ENDIF   
214               IF (j == 1 ) THEN  ! .. and the first grid-row
215                  jm1 = 1
216                  dj = 1.
217               ENDIF
218
219               gradient_ij1 = 0.
220               IF (sou_mask(ip1,jp1) /= il_maskval .OR. &
221                   sou_mask(im1,jp1) /= il_maskval) THEN
222                  IF (sou_mask(ip1,jp1) == il_maskval .AND. &
223                      sou_mask(i,jp1) /= il_maskval) THEN
224                     ip1 = i
225                     di = 1.
226                  ELSE IF (sou_mask(im1,jp1) == il_maskval .AND. &
227                           sou_mask(i,jp1) /= il_maskval) THEN
228                     im1 = i
229                     di = 1.
230                  ELSE
231                     di = 0.
232                  ENDIF
233                  gradient_ij1 = di * (src_array(ip1,jp1) - src_array(im1,jp1))
234               ENDIF
235
236               di = 0.5
237               ip1 = i + 1
238               im1 = i - 1
239               IF (i == NX1) THEN
240                   IF (cd_per == 'P') ip1 = 1 + id_per ! the 0-meridian
241                   IF (cd_per == 'R') ip1 = NX1
242               ENDIF
243               IF (i == 1)  THEN
244                   IF (cd_per == 'P') im1 = NX1 - id_per
245                   IF (cd_per == 'R') im1 = 1
246               ENDIF
247               gradient_ij2 = 0.
248               IF (sou_mask(ip1,jm1) /= il_maskval .OR. &
249                   sou_mask(im1,jm1) /= il_maskval) THEN
250                  IF (sou_mask(ip1,jm1) == il_maskval .AND. &
251                      sou_mask(i,jm1) /= il_maskval) THEN
252                     ip1 = i
253                     di = 1.
254                  ELSE IF (sou_mask(im1,jm1) == il_maskval .AND. &
255                          sou_mask(i,jm1) /= il_maskval) THEN
256                     im1 = i
257                     di = 1.
258                  ELSE
259                     di = 0.
260                  ENDIF
261                  gradient_ij2 = di * (src_array(ip1,jm1) - src_array(im1,jm1))
262               ENDIF
263
264               IF (gradient_ij1 /= 0. .AND. gradient_ij2 /= 0.) THEN
265                  gradient_ij(i,j) = dj * (gradient_ij1 - gradient_ij2)
266               ENDIF
267            ENDIF
268           
269         ENDDO
270      ENDDO
271      !
272      iend = ibeg+iloc-1
273      jend = jbeg+jloc-1
274      local_gradient_i = gradient_i(ibeg:iend, jbeg:jend)
275      local_gradient_j = gradient_j(ibeg:iend, jbeg:jend)
276      local_gradient_ij = gradient_ij(ibeg:iend, jbeg:jend)
277!
278END SUBROUTINE gradient_bicubic
Note: See TracBrowser for help on using the repository browser.