source: CPL/oasis3/trunk/src/lib/scrip/src/gradient_bicubic.f @ 1677

Last change on this file since 1677 was 1677, checked in by aclsce, 12 years ago

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