source: CPL/oasis3-mct/branches/OASIS3-MCT_5.0_branch/examples/regrid_environment/src/gradient_conserv.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.8 KB
Line 
1SUBROUTINE gradient_conserv(NX1, NY1, ibeg, jbeg, iloc, jloc, &
2                         cl_grd_src, id_per, cd_per, w_unit, &
3                         local_grad_lon, local_grad_lat, file_debug)
4
5!
6!**** *gradient_conserv*  - calculate gradients for conservative remapping
7!
8!     Purpose:
9!     -------
10!     Calculation of gradients in latitudinal and longitudinal direction.
11!     In a first step the gradients in direction of source-grid rows 
12!     and lines are calculated. Then they are rotated to longitudinal
13!     and latitudinal direction, using the scalar product.
14!     This routine works for logically rectangular grids, only.
15!
16!**   Interface:
17!     ---------
18!       *CALL*  *gradient_conserv*(NX1, NY1, jbeg, iloc, jloc,
19!                                  cl_grd_src, id_per, cd_per, w_unit,
20!                                  local_grad_lat, local_grad_lon, file_debug)
21!
22!     Input:
23!     -----
24!          NX1            : grid global dimension in x-direction (integer)
25!          NY1            : grid global dimension in y-direction (integer)
26!          ibeg           : start of local domain in global domain in x-direction
27!          jbeg           : start of local domain in global domain in y-direction
28!          iloc           : grid local dimension in x-direction (integer)
29!          jloc           : grid local dimension in y-direction (integer)
30!          cl_grd_src     : grid acronym
31!          id_per         : number of overlapping points for source grid
32!          cd_per         : grip periodicity type
33!          w_unit         : log file unit
34!          file_debug     : logical for activating debug outputs
35!
36!     Output:
37!     ------
38!          local_grad_lon       : gradient in longitudinal direction (real 2D)
39!          local_grad_lat       : gradient in latitudinal direction (real 2D)
40!
41! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
42!
43      USE read_all_data
44      USE function_ana
45!
46      IMPLICIT NONE
47
48      INTEGER, PARAMETER :: wp = SELECTED_REAL_KIND(12,307) ! double
49!-----------------------------------------------------------------------
50!      INTENT(IN)
51!-----------------------------------------------------------------------
52      INTEGER, INTENT(IN) :: &
53          NX1, NY1,  &         ! source grid dimensions
54         ibeg, jbeg, &         ! source grid local start
55         iloc, jloc            ! source grid local dimensions
56
57      CHARACTER(len=4), INTENT(IN) ::  &
58         cl_grd_src            ! grid acronym
59
60      INTEGER, INTENT(IN) :: &
61          id_per,      &       ! nbr of overlapping grid points
62          w_unit               ! log file
63 
64      CHARACTER*8, INTENT(IN) :: &
65          cd_per                ! grip periodicity type     
66
67      LOGICAL, INTENT(IN)      :: file_debug
68
69!-----------------------------------------------------------------------
70!     INTENT(OUT)
71!-----------------------------------------------------------------------
72      REAL (kind=wp), DIMENSION(iloc,jloc), INTENT(OUT) :: &
73           local_grad_lon, &             ! gradient in longitudinal direction
74           local_grad_lat                ! gradient in latitudinal direction
75
76!-----------------------------------------------------------------------
77!     LOCAL VARIABLES
78!-----------------------------------------------------------------------
79      INTEGER :: &
80           i, j, &                ! looping indicees
81           ip1, jp1, im1, jm1, iend, jend
82
83      REAL (kind=wp) :: &
84           distance_rad           ! distance in rad
85     
86      REAL (kind=wp) :: &
87           dVar_i, dVar_j, &      ! difference of Var in i / j direction
88           dlat_i, dlat_j, &      ! difference in lat in i / j direction
89           dlon_i, dlon_j, &      ! difference in lon in i / j direction
90           dist_i, dist_j, &      ! distance in i / j direction
91           grad_i, grad_j, &      ! gradient in i / j direction
92           ABSold, ABSnew, lat_factor
93
94      REAL (kind=wp), DIMENSION(:,:), POINTER :: &
95           src_lon, &             ! source grid longitudes [radiants]
96           src_lat, &             ! source grid latitudes [radiants]
97           src_array, &           ! analytical field
98           grad_lon, &            ! global gradient in latitudinal direction
99           grad_lat               ! global gradient in longitudinal direction
100
101      INTEGER, DIMENSION(:,:), POINTER :: &
102           sou_mask             ! source grid mask
103
104      REAL (kind=wp), PARAMETER :: &
105           pi  = 3.14159265358979323846, &  ! PI
106           pi2 = 2.0d0*pi, &                ! 2PI
107           pi180  = 1.74532925199432957692e-2 ! =PI/180
108
109      INTEGER, PARAMETER ::  il_maskval= 1 ! in our grids sea_value = 0 and land_value = 1
110
111!-----------------------------------------------------------------------
112!
113!     Read global grid and global mask
114!     --------------------------------
115      ALLOCATE(src_lon(NX1, NY1))
116      ALLOCATE(src_lat(NX1, NY1))
117      CALL read_grid(NX1, NY1, 1, 1, NX1, NY1, cl_grd_src, w_unit, src_lon, src_lat, file_debug)
118!
119      ALLOCATE(sou_mask(NX1, NY1))
120      CALL read_mask(NX1, NY1, 1, 1, NX1, NY1, cl_grd_src, w_unit, sou_mask, file_debug)
121
122!     Global field from analytical function
123!     -------------------------------------
124      ALLOCATE(src_array(NX1, NY1))
125#ifdef FANA1
126      CALL function_ana1(NX1, NY1, src_lon, src_lat, src_array)
127#elif defined FANA2
128      CALL function_ana2(NX1, NY1, src_lon, src_lat, src_array)
129#elif defined FANA3
130      CALL function_ana3(NX1, NY1, src_lon, src_lat, src_array)
131#endif
132
133!     Global gradient allocation
134!     --------------------------
135      ALLOCATE(grad_lon(NX1, NY1))
136      ALLOCATE(grad_lat(NX1, NY1))
137
138!     Initialization
139!     --------------
140      grad_lon  = 0.
141      grad_lat  = 0.
142
143!     transformation in radiants for gradient calculation
144!     ---------------------------------------------------
145      src_lon = src_lon * pi180
146      src_lat = src_lat * pi180 
147
148!     calculate gradients
149!     -------------------
150      DO i = 1, NX1
151         DO j = 1, NY1
152                   
153            IF (sou_mask(i,j) /= il_maskval) THEN
154
155               ip1 = i + 1
156               im1 = i - 1
157               IF (i == NX1) THEN
158                   IF (cd_per == 'P') ip1 = 1 + id_per ! the 0-meridian
159                   IF (cd_per == 'R') ip1 = NX1
160               ENDIF
161               IF (i == 1 )  THEN
162                   IF (cd_per == 'P') im1 = NX1 - id_per
163                   IF (cd_per == 'R') im1 = 1
164               ENDIF
165               jp1 = j + 1
166               jm1 = j - 1
167               IF (j == NY1) jp1 = NY1 ! treatment of the last..
168               IF (j == 1 )  jm1 = 1   ! .. and the first grid-row
169
170               IF (sou_mask(ip1,j) == il_maskval)  ip1 = i
171               IF (sou_mask(im1,j) == il_maskval)  im1 = i
172               IF (sou_mask(i,jp1) == il_maskval)  jp1 = j
173               IF (sou_mask(i,jm1) == il_maskval)  jm1 = j         
174
175!              difference between neighbouring datapoints
176               dVar_i = src_array(ip1,j) - src_array(im1,j)
177               dVar_j = src_array(i,jp1) - src_array(i,jm1)
178
179!              difference in latitudes
180               dlat_i = src_lat(ip1,j) - src_lat(im1,j)
181               dlat_j = src_lat(i,jp1) - src_lat(i,jm1)
182
183!              difference in longitudes
184               dlon_i = src_lon(ip1,j) - src_lon(im1,j)
185               IF (dlon_i > pi)  dlon_i = dlon_i - pi2
186               IF (dlon_i < (-pi)) dlon_i = dlon_i + pi2
187               dlon_j = src_lon(i,jp1) - src_lon(i,jm1)
188               IF (dlon_j >   pi)  dlon_j = dlon_j - pi2
189               IF (dlon_j < (-pi)) dlon_j = dlon_j + pi2
190               lat_factor = COS(src_lat(i,j))
191               dlon_i = dlon_i * lat_factor
192               dlon_j = dlon_j * lat_factor
193 
194!              distance
195               dist_i = distance_rad(src_lon(ip1,j), src_lat(ip1,j), &
196                                     src_lon(im1,j), src_lat(im1,j))
197               dist_j = distance_rad(src_lon(i,jp1), src_lat(i,jp1), &
198                                     src_lon(i,jm1), src_lat(i,jm1))
199
200!              gradients: dVar / distance (= vector lenght)
201               IF (dist_i /= 0.) THEN
202                  grad_i = dVar_i / dist_i
203               ELSE
204                  grad_i = 0
205               ENDIF
206               IF (dist_j /= 0.) THEN
207                  grad_j = dVar_j / dist_j
208               ELSE
209                  grad_j = 0
210               ENDIF
211
212!              projection by scalar product
213!              ----------------------------
214               grad_lon(i,j) = grad_i * dlon_i + grad_j * dlat_i
215               grad_lat(i,j) = grad_i * dlon_j + grad_j * dlat_j
216
217               IF (dist_i /= 0) then
218                  grad_lon(i,j) = grad_lon(i,j) / dist_i
219               ELSE
220                  grad_lon(i,j) = 0
221               ENDIF
222               IF (dist_j /= 0) then
223                  grad_lat(i,j) = grad_lat(i,j) / dist_j
224               ELSE
225                  grad_lat(i,j) = 0.
226               ENDIF
227             
228!              correct skale
229!              -------------
230               ABSold = SQRT(grad_i**2 + grad_j**2)
231               ABSnew = SQRT(grad_lon(i,j)**2 + grad_lat(i,j)**2)
232               IF (ABSnew > 1.E-10) THEN
233!                  grad_lon(i,j) = grad_lon(i,j)*ABSold/ABSnew
234                  grad_lon(i,j) = grad_lon(i,j)
235               ELSE
236                  grad_lon(i,j) = 0.0
237               ENDIF
238
239!              test orthogonality
240!              ------------------
241               IF ((dlon_i*dlon_j+dlat_j*dlat_i) > 0.1) THEN
242                  print*, 'ORTHOGONAL? ', i, j, (dlon_i*dlon_j+dlat_j*dlat_i)
243               ENDIF
244
245            ELSE
246           
247               grad_lat(i,j) = 0.
248               grad_lon(i,j) = 0. 
249           
250            ENDIF
251
252         ENDDO
253      ENDDO
254      !
255      iend = ibeg+iloc-1
256      jend = jbeg+jloc-1
257      local_grad_lat = grad_lat(ibeg:iend, jbeg:jend)
258      local_grad_lon = grad_lon(ibeg:iend, jbeg:jend)
259
260END SUBROUTINE gradient_conserv
Note: See TracBrowser for help on using the repository browser.