source: CPL/oasis3-mct/branches/OASIS3-MCT_5.0_branch/examples/spoc/spoc_regridding/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: 9.2 KB
Line 
1subroutine gradient_bicubic(NX1, NY1, src_array, sou_mask, &
2                                  src_latitudes, src_longitudes, &
3                                  id_per, cd_per, &
4                                  gradient_i, gradient_j, gradient_ij)
5!****
6!               *****************************
7!               * OASIS ROUTINE  -  LEVEL ? *
8!               * -------------     ------- *
9!               *****************************
10!
11!**** *gradient_bicubic*  - calculate gradients for bicubic remapping
12!
13!     Purpose:
14!     -------
15!     Calculation of gradients for bicubic interpolation. In contrast to
16!     the gradients of conservative remapping, these gradients are   
17!     calculated with respect to grid rows and grid lines.
18!
19!**   Interface:
20!     ---------
21!       *CALL*  *gradient_bicubic*(NX1, NY1, src_array, sou_mask,
22!                                  src_latitudes, src_longitudes,
23!                                  gradient_i, gradient_j, gradient_ij)
24!
25!     Input:
26!     -----
27!          NX1            : grid dimension in x-direction (integer)
28!          NY1            : grid dimension in y-direction (integer)
29!          src_array      : array on source grid (real 2D)
30!          sou_mask       : source grid mask (integer 2D)
31!          src_latitudes  : latitudes on source grid (real 2D)
32!          src_longitudes : longitudes on source grid (real 2D)
33!          id_per         : number of overlapping points for source grid
34!          cd_per         : grip periodicity type
35!
36!     Output:
37!     ------
38!          gradient_i     : gradient in i-direction (real 2D)
39!          gradient_j     : gradient in j-direction (real 2D)
40!          gradient_ij    : gradient in ij-direction (real 2D)
41!
42!     History:
43!     -------
44!       Version   Programmer     Date        Description
45!       -------   ----------     ----        ----------- 
46!       2.5       V. Gayler      2002/04/05  created
47!
48! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49
50IMPLICIT NONE
51     
52  INTEGER, PARAMETER :: wp = SELECTED_REAL_KIND(12,307) ! double
53!-----------------------------------------------------------------------
54!     INTENT(IN)
55!-----------------------------------------------------------------------
56INTEGER, INTENT(IN) :: &
57         NX1, NY1,  &           ! source grid dimensiones
58         id_per                ! nbr of overlapping grid points
59
60CHARACTER(len=8), INTENT(IN) ::  &
61         cd_per                ! grip periodicity type
62
63REAL (kind=wp), DIMENSION(NX1,NY1), INTENT(IN) :: &
64          src_array            ! array on source grid
65
66INTEGER, DIMENSION(NX1,NY1), INTENT(IN) :: &
67         sou_mask             ! source grid mask
68
69REAL (kind=wp), DIMENSION(NX1,NY1), INTENT(IN) :: &
70          src_latitudes,   &    ! source grid latitudes
71          src_longitudes        ! source grid longitudes
72
73!-----------------------------------------------------------------------
74!     INTENT(OUT)
75!-----------------------------------------------------------------------
76REAL (kind=wp), DIMENSION(NX1,NY1), INTENT(OUT) :: &
77          gradient_i,   &       ! gradient in i-direction (real 2D)
78          gradient_j,   &       ! gradient in j-direction (real 2D)
79          gradient_ij           ! gradient in ij-direction (real 2D)
80
81!-----------------------------------------------------------------------
82!     LOCAL VARIABLES
83!-----------------------------------------------------------------------
84INTEGER ::  &
85          i, j,    &            ! looping indicees
86          ip1, jp1, im1, jm1
87     
88REAL (kind=wp) ::  &
89          di, dj,        &         ! factor depending on grid cell distance
90          gradient_ij1,  &         ! gradient needed to calculate gradient_ij
91          gradient_ij2             ! gradient needed to calculate gradient_ij
92
93REAL (kind=wp), DIMENSION(NX1,NY1) :: &
94          src_lon,   &          ! source grid longitudes [radiants]
95          src_lat,   &          ! source grid latitudes [radiants]
96          pi180                 ! conversion factor: deg -> rad
97
98INTEGER, PARAMETER ::  il_maskval= 1 ! in our grids sea_value = 0 and land_value = 1
99
100!----------------------------------------------------------------------
101!
102!     Transformation from degree to radiant
103!     -------------------------------------
104      pi180 = 1.74532925199432957692e-2 ! =PI/180
105
106      src_lon = src_longitudes * pi180
107      src_lat = src_latitudes * pi180
108
109!     Initialization
110!     --------------
111      gradient_i  = 0.
112      gradient_j  = 0. 
113      gradient_ij = 0. 
114
115!     calculate gradients
116!     -------------------
117      DO i = 1, NX1
118         DO j = 1, NY1
119                   
120            IF (sou_mask (i,j) /= il_maskval) THEN
121
122               di = 0.5
123               dj = 0.5
124
125               ip1 = i + 1
126               im1 = i - 1
127               IF (i == NX1) THEN
128                   IF (cd_per == 'P') ip1 = 1 + id_per ! the 0-meridian
129                   IF (cd_per == 'R') ip1 = NX1
130               ENDIF
131               IF (i == 1 )  THEN
132                   IF (cd_per == 'P') im1 = NX1 - id_per
133                   IF (cd_per == 'R') im1 = 1
134               ENDIF
135               jp1 = j + 1
136               jm1 = j - 1
137               IF (j == NY1) THEN ! treatment of the last..
138                  jp1 = NY1 
139                  dj = 1.
140               ENDIF   
141               IF (j == 1 ) THEN  ! .. and the first grid-row
142                  jm1 = 1
143                  dj = 1.
144               ENDIF
145
146
147!              gradient i
148!              ----------
149               IF (sou_mask(ip1,j) /= il_maskval .OR. &
150                   sou_mask(im1,j) /= il_maskval) THEN
151                  IF (sou_mask(ip1,j) == il_maskval) THEN
152                     ip1 = i
153                     di = 1.
154                  ELSE IF (sou_mask(im1,j) == il_maskval) THEN
155                     im1 = i
156                     di = 1.
157                  ENDIF
158                  gradient_i(i,j) = di * (src_array(ip1,j) - src_array(im1,j))
159               ENDIF
160
161!              gradient j
162!              ----------
163               IF (sou_mask(i,jp1) /= il_maskval .OR. &
164                   sou_mask(i,jm1) /= il_maskval) THEN
165                  IF (sou_mask(i,jp1) == il_maskval) THEN
166                     jp1 = j
167                     dj = 1.
168                  ELSE IF (sou_mask(i,jm1) == il_maskval) THEN
169                     jm1 = j
170                     dj = 1.
171                  ENDIF
172                  gradient_j(i,j) = dj * (src_array(i,jp1) - src_array(i,jm1))
173               ENDIF
174!
175!              gradient ij
176!              -----------
177               di = 0.5
178               dj = 0.5
179               ip1 = i + 1
180               im1 = i - 1
181               IF (i == NX1) THEN
182                   IF (cd_per == 'P') ip1 = 1 + id_per ! the 0-meridian
183                   IF (cd_per == 'R') ip1 = NX1
184               ENDIF
185               IF (i == 1 )  THEN
186                   IF (cd_per == 'P')  im1 = NX1 - id_per
187                   IF (cd_per == 'R')  im1 = 1
188               ENDIF
189               jp1 = j + 1
190               jm1 = j - 1
191               IF (j == NY1) THEN ! treatment of the last..
192                  jp1 = NY1 
193                  dj = 1.
194               ENDIF   
195               IF (j == 1 ) THEN  ! .. and the first grid-row
196                  jm1 = 1
197                  dj = 1.
198               ENDIF
199
200               gradient_ij1 = 0.
201               IF (sou_mask(ip1,jp1) /= il_maskval .OR. &
202                   sou_mask(im1,jp1) /= il_maskval) THEN
203                  IF (sou_mask(ip1,jp1) == il_maskval .AND. &
204                      sou_mask(i,jp1) /= il_maskval) THEN
205                     ip1 = i
206                     di = 1.
207                  ELSE IF (sou_mask(im1,jp1) == il_maskval .AND. &
208                           sou_mask(i,jp1) /= il_maskval) THEN
209                     im1 = i
210                     di = 1.
211                  ELSE
212                     di = 0.
213                  ENDIF
214                  gradient_ij1 = di * (src_array(ip1,jp1) - src_array(im1,jp1))
215               ENDIF
216
217               di = 0.5
218               ip1 = i + 1
219               im1 = i - 1
220               IF (i == NX1) THEN
221                   IF (cd_per == 'P') ip1 = 1 + id_per ! the 0-meridian
222                   IF (cd_per == 'R') ip1 = NX1
223               ENDIF
224               IF (i == 1)  THEN
225                   IF (cd_per == 'P') im1 = NX1 - id_per
226                   IF (cd_per == 'R') im1 = 1
227               ENDIF
228               gradient_ij2 = 0.
229               IF (sou_mask(ip1,jm1) /= il_maskval .OR. &
230                   sou_mask(im1,jm1) /= il_maskval) THEN
231                  IF (sou_mask(ip1,jm1) == il_maskval .AND. &
232                      sou_mask(i,jm1) /= il_maskval) THEN
233                     ip1 = i
234                     di = 1.
235                  ELSE IF (sou_mask(im1,jm1) == il_maskval .AND. &
236                          sou_mask(i,jm1) /= il_maskval) THEN
237                     im1 = i
238                     di = 1.
239                  ELSE
240                     di = 0.
241                  ENDIF
242                  gradient_ij2 = di * (src_array(ip1,jm1) - src_array(im1,jm1))
243               ENDIF
244
245               IF (gradient_ij1 /= 0. .AND. gradient_ij2 /= 0.) THEN
246                  gradient_ij(i,j) = dj * (gradient_ij1 - gradient_ij2)
247               ENDIF
248            ENDIF
249           
250         ENDDO
251      ENDDO
252!
253END SUBROUTINE gradient_bicubic
Note: See TracBrowser for help on using the repository browser.