source: CONFIG_DEVT/IPSLCM6.5_work_ENSEMBLES/oasis3-mct/examples/spoc/spoc_regridding/gradient_conserv.f90 @ 5725

Last change on this file since 5725 was 5725, checked in by aclsce, 3 years ago

Added new oasis3-MCT version to be used to handle ensembles simulations with XIOS.
File size: 9.0 KB
Line 
1SUBROUTINE gradient_conserv(NX1, NY1, src_array, sou_mask, &
2                         src_latitudes, src_longitudes, &
3                         id_per, cd_per, &
4                         grad_lat, grad_lon)
5
6!****
7!               *****************************
8!               * OASIS ROUTINE  -  LEVEL ? *
9!               * -------------     ------- *
10!               *****************************
11!
12!**** *gradient_conserv*  - calculate gradients for conservative remapping
13!
14!     Purpose:
15!     -------
16!     Calculation of gradients in latitudinal and longitudinal direction.
17!     In a first step the gradients in direction of source-grid rows 
18!     and lines are calculated. Then they are rotated to longitudinal
19!     and latitudinal direction, using the scalar product.
20!     This routine works for logically rectangular grids, only.
21!
22!**   Interface:
23!     ---------
24!       *CALL*  *gradient_conserv*(NX1, NY1, src_array, sou_mask, src_latitudes, &
25!                                  src_longitudes, grad_lat, grad_lon)
26!
27!     Input:
28!     -----
29!          NX1            : grid dimension in x-direction (integer)
30!          NY1            : grid dimension in y-direction (integer)
31!          src_array      : array on source grid (real 2D)
32!          sou_mask       : source grid mask (integer 2D)
33!          src_latitudes  : latitudes on source grid (real 2D)
34!          src_longitudes : longitudes on source grid (real 2D)
35!          id_per         : number of overlapping points for source grid
36!          cd_per         : grip periodicity type
37!
38!     Output:
39!     ------
40!          grad_lat       : gradient in latitudinal direction (real 2D)
41!          grad_lon       : gradient in longitudinal direction (real 2D)
42!
43!     History:
44!     -------
45!       Version   Programmer     Date        Description
46!       -------   ----------     ----        ----------- 
47!       2.5       V. Gayler      2001/09/20  created
48!
49! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
50
51IMPLICIT NONE
52
53  INTEGER, PARAMETER :: wp = SELECTED_REAL_KIND(12,307) ! double
54!-----------------------------------------------------------------------
55!     INTENT(IN)
56!-----------------------------------------------------------------------
57      INTEGER, INTENT(IN) :: &
58          NX1, NY1, &             ! source grid dimensions
59          id_per                  ! nbr of overlapping grid points
60
61      CHARACTER*8, INTENT(IN) :: &
62          cd_per                ! grip periodicity type     
63
64      REAL (kind=wp), DIMENSION(NX1,NY1), INTENT(IN) :: &
65          src_array             ! array on source grid
66
67      INTEGER, DIMENSION(NX1,NY1), INTENT(IN) :: &
68          sou_mask              ! source grid mask
69
70      REAL (kind=wp), DIMENSION(NX1,NY1), INTENT(IN) :: &
71          src_latitudes, &        ! source grid latitudes
72          src_longitudes          ! source grid longitudes
73
74!-----------------------------------------------------------------------
75!     INTENT(OUT)
76!-----------------------------------------------------------------------
77      REAL (kind=wp), DIMENSION(NX1,NY1), INTENT(OUT) :: &
78           grad_lat, &            ! gradient in latitudinal direction
79           grad_lon             ! gradient in longitudinal direction
80
81!-----------------------------------------------------------------------
82!     LOCAL VARIABLES
83!-----------------------------------------------------------------------
84      INTEGER :: &
85           i, j, &                ! looping indicees
86           ip1, jp1, im1, jm1
87
88      REAL (kind=wp) :: &
89           distance_rad           ! distance in rad
90     
91      REAL (kind=wp) :: &
92           dVar_i, dVar_j, &      ! difference of Var in i / j direction
93           dlat_i, dlat_j, &      ! difference in lat in i / j direction
94           dlon_i, dlon_j, &      ! difference in lon in i / j direction
95           dist_i, dist_j, &      ! distance in i / j direction
96           grad_i, grad_j, &      ! gradient in i / j direction
97           ABSold, ABSnew, lat_factor
98
99      REAL (kind=wp), DIMENSION(NX1,NY1) :: &
100           src_lon, &             ! source grid longitudes [radiants]
101           src_lat, &             ! source grid latitudes [radiants]
102           pi180                  ! conversion factor: deg -> rad
103      REAL (kind=wp), PARAMETER :: &
104           pi  = 3.14159265358979323846, &  ! PI
105           pi2 = 2.0d0*pi                   ! 2PI
106
107      INTEGER, PARAMETER ::  il_maskval= 1 ! in our grids sea_value = 0 and land_value = 1
108
109!-----------------------------------------------------------------------
110!
111!!      IF (nlogprt .GE. 2) THEN
112!!          WRITE (UNIT = nulou,FMT = *)' '
113!!          WRITE (UNIT = nulou,FMT = *)' Entering routine gradient_conserv   '
114!!          WRITE (UNIT = nulou,FMT = *)' '
115!!          CALL FLUSH(nulou)
116!!      ENDIF
117!
118!     Transformation from degree to radiant
119!     -------------------------------------
120      pi180  = 1.74532925199432957692e-2 ! =PI/180
121
122      src_lon = src_longitudes * pi180
123      src_lat = src_latitudes * pi180
124
125!-----------------------------------------------------------------------
126
127      DO i = 1, NX1
128
129         DO j = 1, NY1
130                   
131            IF (sou_mask(i,j) /= il_maskval) THEN
132
133               ip1 = i + 1
134               im1 = i - 1
135               IF (i == NX1) THEN
136                   IF (cd_per == 'P') ip1 = 1 + id_per ! the 0-meridian
137                   IF (cd_per == 'R') ip1 = NX1
138               ENDIF
139               IF (i == 1 )  THEN
140                   IF (cd_per == 'P') im1 = NX1 - id_per
141                   IF (cd_per == 'R') im1 = 1
142               ENDIF
143               jp1 = j + 1
144               jm1 = j - 1
145               IF (j == NY1) jp1 = NY1 ! treatment of the last..
146               IF (j == 1 )  jm1 = 1   ! .. and the first grid-row
147
148               IF (sou_mask(ip1,j) == il_maskval)  ip1 = i
149               IF (sou_mask(im1,j) == il_maskval)  im1 = i
150               IF (sou_mask(i,jp1) == il_maskval)  jp1 = j
151               IF (sou_mask(i,jm1) == il_maskval)  jm1 = j         
152
153!              difference between neighbouring datapoints
154               dVar_i = src_array(ip1,j) - src_array(im1,j)
155               dVar_j = src_array(i,jp1) - src_array(i,jm1)
156
157!              difference in latitudes
158               dlat_i = src_lat(ip1,j) - src_lat(im1,j)
159               dlat_j = src_lat(i,jp1) - src_lat(i,jm1)
160
161!              difference in longitudes
162               dlon_i = src_lon(ip1,j) - src_lon(im1,j)
163               IF (dlon_i > pi)  dlon_i = dlon_i - pi2
164               IF (dlon_i < (-pi)) dlon_i = dlon_i + pi2
165               dlon_j = src_lon(i,jp1) - src_lon(i,jm1)
166               IF (dlon_j >   pi)  dlon_j = dlon_j - pi2
167               IF (dlon_j < (-pi)) dlon_j = dlon_j + pi2
168               lat_factor = COS(src_lat(i,j))
169               dlon_i = dlon_i * lat_factor
170               dlon_j = dlon_j * lat_factor
171 
172!              distance
173               dist_i = distance_rad(src_lon(ip1,j), src_lat(ip1,j), &
174                                     src_lon(im1,j), src_lat(im1,j))
175               dist_j = distance_rad(src_lon(i,jp1), src_lat(i,jp1), &
176                                     src_lon(i,jm1), src_lat(i,jm1))
177
178!              gradients: dVar / distance (= vector lenght)
179               IF (dist_i /= 0.) THEN
180                  grad_i = dVar_i / dist_i
181               ELSE
182                  grad_i = 0
183               ENDIF
184               IF (dist_j /= 0.) THEN
185                  grad_j = dVar_j / dist_j
186               ELSE
187                  grad_j = 0
188               ENDIF
189
190!              projection by scalar product
191!              ----------------------------
192               grad_lon(i,j) = grad_i * dlon_i + grad_j * dlat_i
193               grad_lat(i,j) = grad_i * dlon_j + grad_j * dlat_j
194
195               IF (dist_i /= 0) then
196                  grad_lon(i,j) = grad_lon(i,j) / dist_i
197               ELSE
198                  grad_lon(i,j) = 0
199               ENDIF
200               IF (dist_j /= 0) then
201                  grad_lat(i,j) = grad_lat(i,j) / dist_j
202               ELSE
203                  grad_lat(i,j) = 0.
204               ENDIF
205             
206!              correct skale
207!              -------------
208               ABSold = SQRT(grad_i**2 + grad_j**2)
209               ABSnew = SQRT(grad_lon(i,j)**2 + grad_lat(i,j)**2)
210               IF (ABSnew > 1.E-10) THEN
211!                  grad_lon(i,j) = grad_lon(i,j)*ABSold/ABSnew
212                  grad_lon(i,j) = grad_lon(i,j)
213               ELSE
214                  grad_lon(i,j) = 0.0
215               ENDIF
216
217!              test orthogonality
218!              ------------------
219               IF ((dlon_i*dlon_j+dlat_j*dlat_i) > 0.1) THEN
220                  print*, 'ORTHOGONAL? ', i, j, (dlon_i*dlon_j+dlat_j*dlat_i)
221               ENDIF
222
223            ELSE
224           
225               grad_lat(i,j) = 0.
226               grad_lon(i,j) = 0. 
227           
228            ENDIF
229
230         ENDDO
231      ENDDO
232!!      IF (nlogprt .GE. 2) THEN
233!!          WRITE (UNIT = nulou,FMT = *)' '
234!!          WRITE (UNIT = nulou,FMT = *)' Leaving routine gradient   '
235!!          WRITE (UNIT = nulou,FMT = *)' '
236!!          CALL FLUSH(nulou)
237!!      ENDIF
238
239END SUBROUTINE gradient_conserv
Note: See TracBrowser for help on using the repository browser.