1 | MODULE cfg_tools |
---|
2 | !!----------------------------------------------------------- |
---|
3 | !! |
---|
4 | !! to make that we use a 4th order polynomial interpolation |
---|
5 | !! |
---|
6 | !! Created by Brice Lemaire on 12/2009. |
---|
7 | !! |
---|
8 | !!----------------------------------------------------------- |
---|
9 | USE readwrite |
---|
10 | USE projection |
---|
11 | ! |
---|
12 | IMPLICIT NONE |
---|
13 | PUBLIC |
---|
14 | ! |
---|
15 | ! |
---|
16 | ! |
---|
17 | CONTAINS |
---|
18 | !******************************************************** |
---|
19 | ! SUBROUTINE interp_grid * |
---|
20 | ! * |
---|
21 | ! calculate polynomial interpolation at 4th order * |
---|
22 | ! * |
---|
23 | ! CALLED from create_coordinates.f90 * |
---|
24 | !******************************************************** |
---|
25 | SUBROUTINE interp_grid |
---|
26 | ! |
---|
27 | REAL*8, DIMENSION(:,:),ALLOCATABLE :: dlcoef !Array to store the coefficients of Lagrange |
---|
28 | INTEGER :: ji, jj, jk, jproj |
---|
29 | INTEGER :: istat, ip |
---|
30 | LOGICAL :: llnorth_pole = .FALSE. |
---|
31 | ! |
---|
32 | WRITE(*,*) '' |
---|
33 | WRITE(*,*) '### SUBROUTINE interp_grid ###' |
---|
34 | WRITE(*,*) '' |
---|
35 | ! |
---|
36 | jproj = 0 |
---|
37 | ! |
---|
38 | ! Calculate coefficients for interpolation along longitude |
---|
39 | ! |
---|
40 | ALLOCATE(dlcoef(nn_rhox-1,4)) |
---|
41 | istat = pol_coef(dlcoef,nn_rhox) |
---|
42 | IF (istat/=1) THEN |
---|
43 | WRITE(*,*) "ERROR WITH LAGRANGIAN COEFFICIENTS" |
---|
44 | STOP |
---|
45 | ENDIF |
---|
46 | ! |
---|
47 | WRITE(*,*) 'Interpolation along longitude' |
---|
48 | ! |
---|
49 | DO jj = nn_rhoy,nygmix,nn_rhoy |
---|
50 | DO ji = nn_rhox,nxgmix,nn_rhox |
---|
51 | ! |
---|
52 | DO jk = 1,nn_rhox-1 |
---|
53 | ! |
---|
54 | ! First, we check the +-180 discontinuity. |
---|
55 | ! In this case, we increase the negative values of 360. |
---|
56 | IF(ABS(smixgrd%glam(ji,jj)-smixgrd%glam(ji+1*nn_rhox,jj)).GT.180.0.AND.smixgrd%glam(ji,jj).LT.0.) THEN |
---|
57 | smixgrd%glam(ji,jj) = smixgrd%glam(ji,jj) + 360. |
---|
58 | ENDIF |
---|
59 | ! |
---|
60 | IF(ABS(smixgrd%glam(ji,jj)-smixgrd%glam(ji+1*nn_rhox,jj)).GT.180.0.AND.smixgrd%glam(ji+1*nn_rhox,jj).LT.0.) THEN |
---|
61 | smixgrd%glam(ji+1*nn_rhox,jj) = smixgrd%glam(ji+1*nn_rhox,jj) + 360. |
---|
62 | ENDIF |
---|
63 | ! |
---|
64 | IF(ABS(smixgrd%glam(ji+1*nn_rhox,jj)-smixgrd%glam(ji+2*nn_rhox,jj)).GT.180.0.AND.smixgrd%glam(ji+1*nn_rhox,jj).LT.0.) THEN |
---|
65 | smixgrd%glam(ji+1*nn_rhox,jj) = smixgrd%glam(ji+1*nn_rhox,jj) + 360. |
---|
66 | ENDIF |
---|
67 | ! |
---|
68 | IF(ABS(smixgrd%glam(ji+1*nn_rhox,jj)-smixgrd%glam(ji+2*nn_rhox,jj)).GT.180.0.AND.smixgrd%glam(ji+2*nn_rhox,jj).LT.0.) THEN |
---|
69 | smixgrd%glam(ji+2*nn_rhox,jj) = smixgrd%glam(ji+2*nn_rhox,jj) + 360. |
---|
70 | ENDIF |
---|
71 | ! |
---|
72 | IF(ABS(smixgrd%glam(ji+2*nn_rhox,jj)-smixgrd%glam(ji+3*nn_rhox,jj)).GT.180.0.AND.smixgrd%glam(ji+2*nn_rhox,jj).LT.0.) THEN |
---|
73 | smixgrd%glam(ji+2*nn_rhox,jj) = smixgrd%glam(ji+2*nn_rhox,jj) + 360. |
---|
74 | ENDIF |
---|
75 | ! |
---|
76 | IF(ABS(smixgrd%glam(ji+2*nn_rhox,jj)-smixgrd%glam(ji+3*nn_rhox,jj)).GT.180.0.AND.smixgrd%glam(ji+3*nn_rhox,jj).LT.0.) THEN |
---|
77 | smixgrd%glam(ji+3*nn_rhox,jj) = smixgrd%glam(ji+3*nn_rhox,jj) + 360. |
---|
78 | ENDIF |
---|
79 | ! |
---|
80 | ! If we are along north boundary, |
---|
81 | ! the variation of longitude looks like a heaviside fonction at the geographical north pole. |
---|
82 | ! Thus, we can't make an interpolation. |
---|
83 | IF(ABS(smixgrd%glam(ji,jj) - smixgrd%glam(ji+3*nn_rhox,jj)).EQ.180.0)THEN |
---|
84 | llnorth_pole = .TRUE. |
---|
85 | ENDIF |
---|
86 | ! |
---|
87 | ! Nearby the geographical north pole, |
---|
88 | ! the variation of the longitudes is too important. |
---|
89 | ! We need to make a polar stereographic projection before interpolation. |
---|
90 | !IF(.NOT.llnorth_pole.AND.ABS(smixgrd%glam(ji,jj)-smixgrd%glam(ji+3*nn_rhox,jj)).GE.80.0) THEN |
---|
91 | IF(.NOT.llnorth_pole.AND.smixgrd%gphi(ji,jj).GE.88.) THEN |
---|
92 | CALL stereo_projection(ji,jj,jk,llnorth_pole,1) |
---|
93 | jproj = 1 |
---|
94 | ENDIF |
---|
95 | ! |
---|
96 | smixgrd%glam(ji+nn_rhox+jk,jj) = dlcoef(jk,1) * smixgrd%glam(ji,jj) & |
---|
97 | + dlcoef(jk,2) * smixgrd%glam(ji+1*nn_rhox,jj) & |
---|
98 | + dlcoef(jk,3) * smixgrd%glam(ji+2*nn_rhox,jj) & |
---|
99 | + dlcoef(jk,4) * smixgrd%glam(ji+3*nn_rhox,jj) |
---|
100 | ! |
---|
101 | smixgrd%gphi(ji+nn_rhox+jk,jj) = dlcoef(jk,1) * smixgrd%gphi(ji,jj) & |
---|
102 | + dlcoef(jk,2) * smixgrd%gphi(ji+1*nn_rhox,jj) & |
---|
103 | + dlcoef(jk,3) * smixgrd%gphi(ji+2*nn_rhox,jj) & |
---|
104 | + dlcoef(jk,4) * smixgrd%gphi(ji+3*nn_rhox,jj) |
---|
105 | ! |
---|
106 | smixgrd%e1(ji+nn_rhox+jk,jj) = dlcoef(jk,1) * smixgrd%e1(ji,jj) & |
---|
107 | + dlcoef(jk,2) * smixgrd%e1(ji+1*nn_rhox,jj) & |
---|
108 | + dlcoef(jk,3) * smixgrd%e1(ji+2*nn_rhox,jj) & |
---|
109 | + dlcoef(jk,4) * smixgrd%e1(ji+3*nn_rhox,jj) |
---|
110 | ! |
---|
111 | smixgrd%e2(ji+nn_rhox+jk,jj) = dlcoef(jk,1) * smixgrd%e2(ji,jj) & |
---|
112 | + dlcoef(jk,2) * smixgrd%e2(ji+1*nn_rhox,jj) & |
---|
113 | + dlcoef(jk,3) * smixgrd%e2(ji+2*nn_rhox,jj) & |
---|
114 | + dlcoef(jk,4) * smixgrd%e2(ji+3*nn_rhox,jj) |
---|
115 | ! |
---|
116 | smixgrd%nav_lon(ji+nn_rhox+jk,jj) = smixgrd%glam(ji+nn_rhox+jk,jj) |
---|
117 | smixgrd%nav_lat(ji+nn_rhox+jk,jj) = smixgrd%gphi(ji+nn_rhox+jk,jj) |
---|
118 | ! |
---|
119 | ! We make the polar stereographic projection reverse if needs. |
---|
120 | IF(jproj.EQ.1)THEN |
---|
121 | CALL stereo_projection_inv(ji,jj,jk,llnorth_pole,1) |
---|
122 | ENDIF |
---|
123 | ! |
---|
124 | ! We replace the strong values along the north boundary. |
---|
125 | IF(llnorth_pole)THEN |
---|
126 | IF(smixgrd%glam(ji+1*nn_rhox,jj).EQ.smixgrd%glam(ji+2*nn_rhox,jj))THEN |
---|
127 | smixgrd%glam(ji+nn_rhox+jk,jj) = smixgrd%glam(ji+1*nn_rhox,jj) |
---|
128 | ELSEIF(ABS(smixgrd%glam(ji+1*nn_rhox,jj) - smixgrd%glam(ji+2*nn_rhox,jj)).EQ.180.0)THEN |
---|
129 | IF(smixgrd%gphi(ji+1*nn_rhox,jj).LT.smixgrd%gphi(ji+2*nn_rhox,jj))THEN |
---|
130 | smixgrd%glam(ji+nn_rhox+jk,jj) = smixgrd%glam(ji+1*nn_rhox,jj) |
---|
131 | ELSEIF(smixgrd%gphi(ji+1*nn_rhox,jj).GT.smixgrd%gphi(ji+2*nn_rhox,jj))THEN |
---|
132 | smixgrd%glam(ji+nn_rhox+jk,jj) = smixgrd%glam(ji+2*nn_rhox,jj) |
---|
133 | ENDIF |
---|
134 | ENDIF |
---|
135 | ENDIF |
---|
136 | ! |
---|
137 | ip = 0 |
---|
138 | jproj = 0 |
---|
139 | llnorth_pole = .FALSE. |
---|
140 | ! |
---|
141 | END DO |
---|
142 | END DO |
---|
143 | END DO |
---|
144 | ! |
---|
145 | WHERE(smixgrd%glam.GT.180) |
---|
146 | smixgrd%glam = smixgrd%glam - 360.0 |
---|
147 | ENDWHERE |
---|
148 | ! |
---|
149 | DEALLOCATE(dlcoef) |
---|
150 | ! |
---|
151 | ! Calculate coefficients for interpolation along latitude |
---|
152 | ! |
---|
153 | ALLOCATE(dlcoef(nn_rhoy-1,4)) |
---|
154 | istat = pol_coef(dlcoef,nn_rhoy) |
---|
155 | IF (istat/=1) THEN |
---|
156 | WRITE(*,*) "ERROR WITH LAGRANGIAN COEFFICIENTS" |
---|
157 | STOP |
---|
158 | ENDIF |
---|
159 | ! |
---|
160 | WRITE(*,*) 'Interpolation along latitude' |
---|
161 | ! |
---|
162 | print*, nequator |
---|
163 | DO ji = 1,nxgmix,1 |
---|
164 | ! |
---|
165 | DO jj = nn_rhoy,nygmix,nn_rhoy |
---|
166 | ! |
---|
167 | DO jk = 1,nn_rhoy-1 |
---|
168 | ! |
---|
169 | IF(ABS(smixgrd%glam(ji,jj)-smixgrd%glam(ji,jj+1*nn_rhoy)).GT.180.0.AND.smixgrd%glam(ji,jj).LT.0.) THEN |
---|
170 | smixgrd%glam(ji,jj) = smixgrd%glam(ji,jj) + 360. |
---|
171 | ENDIF |
---|
172 | ! |
---|
173 | IF(ABS(smixgrd%glam(ji,jj)-smixgrd%glam(ji,jj+1*nn_rhoy)).GT.180.0.AND.smixgrd%glam(ji,jj+1*nn_rhoy).LT.0.) THEN |
---|
174 | smixgrd%glam(ji,jj+1*nn_rhoy) = smixgrd%glam(ji,jj+1*nn_rhoy) + 360. |
---|
175 | ENDIF |
---|
176 | ! |
---|
177 | IF(ABS(smixgrd%glam(ji,jj+1*nn_rhoy)-smixgrd%glam(ji,jj+2*nn_rhoy)).GT.180.0.AND.smixgrd%glam(ji,jj+1*nn_rhoy).LT.0.) THEN |
---|
178 | smixgrd%glam(ji,jj+1*nn_rhoy) = smixgrd%glam(ji,jj+1*nn_rhoy) + 360. |
---|
179 | ENDIF |
---|
180 | ! |
---|
181 | IF(ABS(smixgrd%glam(ji,jj+1*nn_rhoy)-smixgrd%glam(ji,jj+2*nn_rhoy)).GT.180.0.AND.smixgrd%glam(ji,jj+2*nn_rhoy).LT.0.) THEN |
---|
182 | smixgrd%glam(ji,jj+2*nn_rhoy) = smixgrd%glam(ji,jj+2*nn_rhoy) + 360. |
---|
183 | ENDIF |
---|
184 | ! |
---|
185 | IF(ABS(smixgrd%glam(ji,jj+2*nn_rhoy)-smixgrd%glam(ji,jj+3*nn_rhoy)).GT.180.0.AND.smixgrd%glam(ji,jj+2*nn_rhoy).LT.0.) THEN |
---|
186 | smixgrd%glam(ji,jj+2*nn_rhoy) = smixgrd%glam(ji,jj+2*nn_rhoy) + 360. |
---|
187 | ENDIF |
---|
188 | ! |
---|
189 | IF(ABS(smixgrd%glam(ji,jj+2*nn_rhoy)-smixgrd%glam(ji,jj+3*nn_rhoy)).GT.180.0.AND.smixgrd%glam(ji,jj+3*nn_rhoy).LT.0.) THEN |
---|
190 | smixgrd%glam(ji,jj+3*nn_rhoy) = smixgrd%glam(ji,jj+3*nn_rhoy) + 360. |
---|
191 | ENDIF |
---|
192 | ! |
---|
193 | !IF(.NOT.llnorth_pole.AND.ABS(smixgrd%glam(ji,jj)-smixgrd%glam(ji,jj+3*nn_rhoy)).GE.60.0) THEN |
---|
194 | IF(.NOT.llnorth_pole.AND.smixgrd%gphi(ji,jj).GE.88.) THEN |
---|
195 | CALL stereo_projection(ji,jj,jk,llnorth_pole,2) |
---|
196 | jproj = 1 |
---|
197 | ENDIF |
---|
198 | ! |
---|
199 | smixgrd%glam(ji,jj+nn_rhoy+jk) = dlcoef(jk,1) * smixgrd%glam(ji,jj) & |
---|
200 | + dlcoef(jk,2) * smixgrd%glam(ji,jj+nn_rhoy) & |
---|
201 | + dlcoef(jk,3) * smixgrd%glam(ji,jj+2*nn_rhoy) & |
---|
202 | + dlcoef(jk,4) * smixgrd%glam(ji,jj+3*nn_rhoy) |
---|
203 | ! |
---|
204 | smixgrd%gphi(ji,jj+nn_rhoy+jk) = dlcoef(jk,1) * smixgrd%gphi(ji,jj) & |
---|
205 | + dlcoef(jk,2) * smixgrd%gphi(ji,jj+nn_rhoy) & |
---|
206 | + dlcoef(jk,3) * smixgrd%gphi(ji,jj+2*nn_rhoy) & |
---|
207 | + dlcoef(jk,4) * smixgrd%gphi(ji,jj+3*nn_rhoy) |
---|
208 | ! |
---|
209 | smixgrd%e1(ji,jj+nn_rhoy+jk) = dlcoef(jk,1) * smixgrd%e1(ji,jj) & |
---|
210 | + dlcoef(jk,2) * smixgrd%e1(ji,jj+nn_rhoy) & |
---|
211 | + dlcoef(jk,3) * smixgrd%e1(ji,jj+2*nn_rhoy) & |
---|
212 | + dlcoef(jk,4) * smixgrd%e1(ji,jj+3*nn_rhoy) |
---|
213 | ! |
---|
214 | smixgrd%e2(ji,jj+nn_rhoy+jk) = dlcoef(jk,1) * smixgrd%e2(ji,jj) & |
---|
215 | + dlcoef(jk,2) * smixgrd%e2(ji,jj+nn_rhoy) & |
---|
216 | + dlcoef(jk,3) * smixgrd%e2(ji,jj+2*nn_rhoy) & |
---|
217 | + dlcoef(jk,4) * smixgrd%e2(ji,jj+3*nn_rhoy) |
---|
218 | ! |
---|
219 | smixgrd%nav_lon(ji,jj+nn_rhoy+jk) = smixgrd%glam(ji,jj+nn_rhoy+jk) |
---|
220 | smixgrd%nav_lat(ji,jj+nn_rhoy+jk) = smixgrd%gphi(ji,jj+nn_rhoy+jk) |
---|
221 | ! |
---|
222 | IF(jproj.EQ.1)THEN |
---|
223 | CALL stereo_projection_inv(ji,jj,jk,llnorth_pole,2) |
---|
224 | ENDIF |
---|
225 | ! |
---|
226 | jproj = 0 |
---|
227 | llnorth_pole = .FALSE. |
---|
228 | ! |
---|
229 | END DO |
---|
230 | ! |
---|
231 | IF(jj+3*nn_rhoy.EQ.nygmix) EXIT |
---|
232 | ! |
---|
233 | END DO |
---|
234 | ! |
---|
235 | END DO |
---|
236 | ! |
---|
237 | WHERE(smixgrd%glam.GT.180.) |
---|
238 | smixgrd%glam = smixgrd%glam - 360.0 |
---|
239 | ENDWHERE |
---|
240 | ! |
---|
241 | WRITE(*,*) '' |
---|
242 | WRITE(*,*) '### END SUBROUTINE interp_grid ###' |
---|
243 | WRITE(*,*) '' |
---|
244 | ! |
---|
245 | END SUBROUTINE |
---|
246 | ! |
---|
247 | ! |
---|
248 | ! |
---|
249 | !******************************************************** |
---|
250 | ! FUNCTION pol_coef * |
---|
251 | ! * |
---|
252 | ! calculate the coefficients of Lagrange * |
---|
253 | ! for the polynomial interpolation * |
---|
254 | ! * |
---|
255 | ! CALLED from SUBROUTINE interp * |
---|
256 | !******************************************************** |
---|
257 | REAL FUNCTION pol_coef(kvect,kxy) |
---|
258 | ! |
---|
259 | REAL*8, DIMENSION(:,:),ALLOCATABLE :: kvect |
---|
260 | REAL*8, DIMENSION(3) :: dlv |
---|
261 | INTEGER :: ji, jm, jk |
---|
262 | REAL*8 :: dlx0 !position relative du point calculer |
---|
263 | INTEGER :: jx_k, jx_i !position relative des points utiliss pour l'interpolation |
---|
264 | REAL*8 :: dleps |
---|
265 | INTEGER :: kxy, irho |
---|
266 | ! |
---|
267 | !on parle de position relative puisque nous utilisons les positions |
---|
268 | !indiciaires, lesquelles sont rptes dans toute la grille. |
---|
269 | !Il n'est donc ncessaire de calculer qu'une fois les 4 coefficients |
---|
270 | !qui seront utiliss dans toute la grille en fonction de nn_rho |
---|
271 | ! |
---|
272 | WRITE(*,*) '' |
---|
273 | WRITE(*,*) '*** FUNCTION pol_coef ***' |
---|
274 | WRITE(*,*) '' |
---|
275 | ! |
---|
276 | jm=1 |
---|
277 | dleps = 1.-1e-8 |
---|
278 | ! |
---|
279 | irho = kxy |
---|
280 | ! |
---|
281 | DO jk = 1,irho-1 |
---|
282 | dlx0 = irho+1+jk |
---|
283 | ji=1 |
---|
284 | DO jx_i = 1,4+3*(irho-1),irho |
---|
285 | jm=1 |
---|
286 | DO jx_k=1,4+3*(irho-1),irho |
---|
287 | IF(jx_k == jx_i) THEN |
---|
288 | CYCLE |
---|
289 | ELSE |
---|
290 | dlv(jm) = (dlx0-jx_k) / (jx_i-jx_k) |
---|
291 | jm = jm + 1 |
---|
292 | END IF |
---|
293 | END DO |
---|
294 | kvect(jk,ji) = product(dlv) |
---|
295 | ji = ji + 1 |
---|
296 | END DO |
---|
297 | END DO |
---|
298 | ! |
---|
299 | IF(SUM(kvect).LT.dleps .OR. SUM(kvect).GT.dleps) THEN |
---|
300 | WRITE(*,*) '' |
---|
301 | WRITE(*,*) '*** CHECK LAGRANGE COEFFICIENTS: ***' |
---|
302 | WRITE(*,*) '' |
---|
303 | ! |
---|
304 | DO ji=1,irho-1 |
---|
305 | WRITE(*,*)'point #',ji |
---|
306 | WRITE(*,*) 'dlcoef(:)= ', kvect(ji,:) |
---|
307 | WRITE(*,*) 'SUM(dlcoef(:)) =', SUM(kvect(ji,:)) |
---|
308 | WRITE(*,*)'' |
---|
309 | END DO |
---|
310 | pol_coef = 1 |
---|
311 | ELSE |
---|
312 | ! |
---|
313 | pol_coef = 0 |
---|
314 | ENDIF |
---|
315 | ! |
---|
316 | END FUNCTION pol_coef |
---|
317 | ! |
---|
318 | ! |
---|
319 | ! |
---|
320 | !******************************************************** |
---|
321 | ! SUBROUTINE child_grid * |
---|
322 | ! * |
---|
323 | ! create the child grids from mixed grid * |
---|
324 | ! * |
---|
325 | ! CALLED from create_coordinates.f90 * |
---|
326 | !******************************************************** |
---|
327 | SUBROUTINE child_grid |
---|
328 | ! |
---|
329 | INTEGER :: ji, jj |
---|
330 | INTEGER :: ip, iq |
---|
331 | REAL*8, DIMENSION(2) :: dlgrdt, dlgrdu, dlgrdv, dlgrdf |
---|
332 | REAL*8 :: dleps |
---|
333 | ! |
---|
334 | WRITE(*,*) '' |
---|
335 | WRITE(*,*) '### SUBROUTINE child_grid ###' |
---|
336 | WRITE(*,*) '' |
---|
337 | ! |
---|
338 | IF(.NOT.nglobal.AND.(nn_rhox.GT.1.OR.nn_rhoy.GT.1)) THEN |
---|
339 | iq=1 |
---|
340 | DO jj = (2*nn_rhoy)+1-nequator,(nygmix-(2*nn_rhoy-1)-nequator),2 |
---|
341 | ip=1 |
---|
342 | DO ji = (2*nn_rhox)+1,(nxgmix-(2*nn_rhox-1)),2 |
---|
343 | ! |
---|
344 | sfingrd%nav_lon(ip,iq) = smixgrd%nav_lon(ji,jj) |
---|
345 | sfingrd%nav_lat(ip,iq) = smixgrd%nav_lat(ji,jj) |
---|
346 | ! |
---|
347 | sfingrd%glamt(ip,iq) = smixgrd%glam(ji,jj) |
---|
348 | sfingrd%glamu(ip,iq) = smixgrd%glam(ji+1,jj) |
---|
349 | sfingrd%glamv(ip,iq) = smixgrd%glam(ji,jj+1) |
---|
350 | sfingrd%glamf(ip,iq) = smixgrd%glam(ji+1,jj+1) |
---|
351 | ! |
---|
352 | sfingrd%gphit(ip,iq) = smixgrd%gphi(ji,jj) |
---|
353 | sfingrd%gphiu(ip,iq) = smixgrd%gphi(ji+1,jj) |
---|
354 | sfingrd%gphiv(ip,iq) = smixgrd%gphi(ji,jj+1) |
---|
355 | sfingrd%gphif(ip,iq) = smixgrd%gphi(ji+1,jj+1) |
---|
356 | ! |
---|
357 | sfingrd%e1t(ip,iq) = smixgrd%e1(ji,jj) |
---|
358 | sfingrd%e1u(ip,iq) = smixgrd%e1(ji+1,jj) |
---|
359 | sfingrd%e1v(ip,iq) = smixgrd%e1(ji,jj+1) |
---|
360 | sfingrd%e1f(ip,iq) = smixgrd%e1(ji+1,jj+1) |
---|
361 | ! |
---|
362 | sfingrd%e2t(ip,iq) = smixgrd%e2(ji,jj) |
---|
363 | sfingrd%e2u(ip,iq) = smixgrd%e2(ji+1,jj) |
---|
364 | sfingrd%e2v(ip,iq) = smixgrd%e2(ji,jj+1) |
---|
365 | sfingrd%e2f(ip,iq) = smixgrd%e2(ji+1,jj+1) |
---|
366 | ! |
---|
367 | ip=ip+1 |
---|
368 | ENDDO |
---|
369 | iq=iq+1 |
---|
370 | ENDDO |
---|
371 | ! |
---|
372 | ELSEIF(nglobal.AND.(nn_rhox.GT.1.OR.nn_rhoy.GT.1))THEN |
---|
373 | iq=1 |
---|
374 | DO jj = (2*nn_rhoy)+1-nequator,nygmix,2 |
---|
375 | ip=1 |
---|
376 | DO ji = (2*nn_rhox)-1,nxgmix,2 |
---|
377 | ! |
---|
378 | sfingrd%nav_lon(ip,iq) = smixgrd%nav_lon(ji,jj) |
---|
379 | sfingrd%nav_lat(ip,iq) = smixgrd%nav_lat(ji,jj) |
---|
380 | ! |
---|
381 | sfingrd%glamt(ip,iq) = smixgrd%glam(ji,jj) |
---|
382 | sfingrd%glamu(ip,iq) = smixgrd%glam(ji+1,jj) |
---|
383 | sfingrd%glamv(ip,iq) = smixgrd%glam(ji,jj+1) |
---|
384 | sfingrd%glamf(ip,iq) = smixgrd%glam(ji+1,jj+1) |
---|
385 | ! |
---|
386 | sfingrd%gphit(ip,iq) = smixgrd%gphi(ji,jj) |
---|
387 | sfingrd%gphiu(ip,iq) = smixgrd%gphi(ji+1,jj) |
---|
388 | sfingrd%gphiv(ip,iq) = smixgrd%gphi(ji,jj+1) |
---|
389 | sfingrd%gphif(ip,iq) = smixgrd%gphi(ji+1,jj+1) |
---|
390 | ! |
---|
391 | sfingrd%e1t(ip,iq) = smixgrd%e1(ji,jj) |
---|
392 | sfingrd%e1u(ip,iq) = smixgrd%e1(ji+1,jj) |
---|
393 | sfingrd%e1v(ip,iq) = smixgrd%e1(ji,jj+1) |
---|
394 | sfingrd%e1f(ip,iq) = smixgrd%e1(ji+1,jj+1) |
---|
395 | ! |
---|
396 | sfingrd%e2t(ip,iq) = smixgrd%e2(ji,jj) |
---|
397 | sfingrd%e2u(ip,iq) = smixgrd%e2(ji+1,jj) |
---|
398 | sfingrd%e2v(ip,iq) = smixgrd%e2(ji,jj+1) |
---|
399 | sfingrd%e2f(ip,iq) = smixgrd%e2(ji+1,jj+1) |
---|
400 | ! |
---|
401 | IF(ip.EQ.nxfine) EXIT |
---|
402 | ! |
---|
403 | ip=ip+1 |
---|
404 | ! |
---|
405 | ENDDO |
---|
406 | ! |
---|
407 | IF(iq.EQ.nyfine) EXIT |
---|
408 | ! |
---|
409 | iq=iq+1 |
---|
410 | ! |
---|
411 | ENDDO |
---|
412 | ! |
---|
413 | ELSE !No interpolation |
---|
414 | iq=1 |
---|
415 | DO jj = 1,nygmix-1,2 |
---|
416 | ip=1 |
---|
417 | DO ji = 1,nxgmix-1,2 |
---|
418 | ! |
---|
419 | sfingrd%nav_lon(ip,iq) = smixgrd%nav_lon(ji,jj) |
---|
420 | sfingrd%nav_lat(ip,iq) = smixgrd%nav_lat(ji,jj) |
---|
421 | ! |
---|
422 | sfingrd%glamt(ip,iq) = smixgrd%glam(ji,jj) |
---|
423 | sfingrd%glamu(ip,iq) = smixgrd%glam(ji+1,jj) |
---|
424 | sfingrd%glamv(ip,iq) = smixgrd%glam(ji,jj+1) |
---|
425 | sfingrd%glamf(ip,iq) = smixgrd%glam(ji+1,jj+1) |
---|
426 | ! |
---|
427 | sfingrd%gphit(ip,iq) = smixgrd%gphi(ji,jj) |
---|
428 | sfingrd%gphiu(ip,iq) = smixgrd%gphi(ji+1,jj) |
---|
429 | sfingrd%gphiv(ip,iq) = smixgrd%gphi(ji,jj+1) |
---|
430 | sfingrd%gphif(ip,iq) = smixgrd%gphi(ji+1,jj+1) |
---|
431 | ! |
---|
432 | sfingrd%e1t(ip,iq) = smixgrd%e1(ji,jj) |
---|
433 | sfingrd%e1u(ip,iq) = smixgrd%e1(ji+1,jj) |
---|
434 | sfingrd%e1v(ip,iq) = smixgrd%e1(ji,jj+1) |
---|
435 | sfingrd%e1f(ip,iq) = smixgrd%e1(ji+1,jj+1) |
---|
436 | ! |
---|
437 | sfingrd%e2t(ip,iq) = smixgrd%e2(ji,jj) |
---|
438 | sfingrd%e2u(ip,iq) = smixgrd%e2(ji+1,jj) |
---|
439 | sfingrd%e2v(ip,iq) = smixgrd%e2(ji,jj+1) |
---|
440 | sfingrd%e2f(ip,iq) = smixgrd%e2(ji+1,jj+1) |
---|
441 | ! |
---|
442 | ip=ip+1 |
---|
443 | END DO |
---|
444 | iq=iq+1 |
---|
445 | END DO |
---|
446 | ! |
---|
447 | ENDIF |
---|
448 | ! |
---|
449 | ! With a global domain, we check the overlap bands for have |
---|
450 | ! * Grid(1,:) = Grid(n-1,:) |
---|
451 | ! * Grid(2,:) = Grid(n,:) |
---|
452 | ! because after interpolation the first column have no values |
---|
453 | IF(nglobal.AND.nn_rhox.GT.1)THEN |
---|
454 | ! |
---|
455 | sfingrd%nav_lon(1,:) = sfingrd%nav_lon(nxfine-1,:) |
---|
456 | sfingrd%nav_lat(1,:) = sfingrd%nav_lat(nxfine-1,:) |
---|
457 | ! |
---|
458 | sfingrd%glamt(1,:) = sfingrd%glamt(nxfine-1,:) |
---|
459 | sfingrd%glamu(1,:) = sfingrd%glamu(nxfine-1,:) |
---|
460 | sfingrd%glamv(1,:) = sfingrd%glamv(nxfine-1,:) |
---|
461 | sfingrd%glamf(1,:) = sfingrd%glamf(nxfine-1,:) |
---|
462 | ! |
---|
463 | sfingrd%gphit(1,:) = sfingrd%gphit(nxfine-1,:) |
---|
464 | sfingrd%gphiu(1,:) = sfingrd%gphiu(nxfine-1,:) |
---|
465 | sfingrd%gphiv(1,:) = sfingrd%gphiv(nxfine-1,:) |
---|
466 | sfingrd%gphif(1,:) = sfingrd%gphif(nxfine-1,:) |
---|
467 | ! |
---|
468 | sfingrd%e1t(1,:) = sfingrd%e1t(nxfine-1,:) |
---|
469 | sfingrd%e1u(1,:) = sfingrd%e1u(nxfine-1,:) |
---|
470 | sfingrd%e1v(1,:) = sfingrd%e1v(nxfine-1,:) |
---|
471 | sfingrd%e1f(1,:) = sfingrd%e1f(nxfine-1,:) |
---|
472 | ! |
---|
473 | sfingrd%e2t(1,:) = sfingrd%e2t(nxfine-1,:) |
---|
474 | sfingrd%e2u(1,:) = sfingrd%e2u(nxfine-1,:) |
---|
475 | sfingrd%e2v(1,:) = sfingrd%e2v(nxfine-1,:) |
---|
476 | sfingrd%e2f(1,:) = sfingrd%e2f(nxfine-1,:) |
---|
477 | ! |
---|
478 | WRITE(*,*) '' |
---|
479 | WRITE(*,*) 'WE CHECK THE OVERLAP BANDS FOR EACH FINE GRID (T,U,V & F):' |
---|
480 | WRITE(*,*) ' ==> SUM{ grd(1,:) + grd(2,:) } - SUM{ grd(n-1) + grd(n,:) } = 0' |
---|
481 | ! |
---|
482 | dleps = 1e-2 |
---|
483 | ! |
---|
484 | WRITE(*,*) '* grid T:' |
---|
485 | dlgrdt(1) = SUM(sfingrd%glamt(1,:)) + SUM(sfingrd%glamt(2,:)) |
---|
486 | dlgrdt(1) = dlgrdt(1) + SUM(sfingrd%gphit(1,:)) + SUM(sfingrd%gphit(2,:)) |
---|
487 | dlgrdt(1) = dlgrdt(1) + SUM(sfingrd%e1t(1,:)) + SUM(sfingrd%e1t(2,:)) |
---|
488 | dlgrdt(1) = dlgrdt(1) + SUM(sfingrd%e2t(1,:)) + SUM(sfingrd%e2t(2,:)) |
---|
489 | ! |
---|
490 | dlgrdt(2) = SUM(sfingrd%glamt(nxfine-1,:)) + SUM(sfingrd%glamt(nxfine,:)) |
---|
491 | dlgrdt(2) = dlgrdt(2) + SUM(sfingrd%gphit(nxfine-1,:)) + SUM(sfingrd%gphit(nxfine,:)) |
---|
492 | dlgrdt(2) = dlgrdt(2) + SUM(sfingrd%e1t(nxfine-1,:)) + SUM(sfingrd%e1t(nxfine,:)) |
---|
493 | dlgrdt(2) = dlgrdt(2) + SUM(sfingrd%e2t(nxfine-1,:)) + SUM(sfingrd%e2t(nxfine,:)) |
---|
494 | ! |
---|
495 | IF((dlgrdt(1)-dlgrdt(2)).GT.dleps.OR.(dlgrdt(1)+dlgrdt(2)).LT.dleps) THEN |
---|
496 | WRITE(*,*) ' ERROR' |
---|
497 | print*,(dlgrdt(1)-dlgrdt(2)), (dlgrdt(1)+dlgrdt(2)) |
---|
498 | ELSE |
---|
499 | WRITE(*,*) 'OVERLAP BANDS OK' |
---|
500 | ENDIF |
---|
501 | ! |
---|
502 | WRITE(*,*) '* grid U:' |
---|
503 | dlgrdu(1) = SUM(sfingrd%glamu(1,:)) + SUM(sfingrd%glamu(2,:)) |
---|
504 | dlgrdu(1) = dlgrdu(1) + SUM(sfingrd%gphiu(1,:)) + SUM(sfingrd%gphiu(2,:)) |
---|
505 | dlgrdu(1) = dlgrdu(1) + SUM(sfingrd%e1u(1,:)) + SUM(sfingrd%e1u(2,:)) |
---|
506 | dlgrdu(1) = dlgrdu(1) + SUM(sfingrd%e2u(1,:)) + SUM(sfingrd%e2u(2,:)) |
---|
507 | ! |
---|
508 | dlgrdu(2) = SUM(sfingrd%glamu(nxfine-1,:)) + SUM(sfingrd%glamu(nxfine,:)) |
---|
509 | dlgrdu(2) = dlgrdu(2) + SUM(sfingrd%gphiu(nxfine-1,:)) + SUM(sfingrd%gphiu(nxfine,:)) |
---|
510 | dlgrdu(2) = dlgrdu(2) + SUM(sfingrd%e1u(nxfine-1,:)) + SUM(sfingrd%e1u(nxfine,:)) |
---|
511 | dlgrdu(2) = dlgrdu(2) + SUM(sfingrd%e2u(nxfine-1,:)) + SUM(sfingrd%e2u(nxfine,:)) |
---|
512 | ! |
---|
513 | IF((dlgrdu(1)-dlgrdu(2)).GT.dleps.OR.(dlgrdu(1)+dlgrdu(2)).LT.dleps) THEN |
---|
514 | WRITE(*,*) ' ERROR' |
---|
515 | print*,(dlgrdu(1)-dlgrdu(2)), (dlgrdu(1)+dlgrdu(2)) |
---|
516 | ELSE |
---|
517 | WRITE(*,*) 'OVERLAP BANDS OK' |
---|
518 | ENDIF |
---|
519 | ! |
---|
520 | WRITE(*,*) '* grid V:' |
---|
521 | dlgrdv(1) = SUM(sfingrd%glamv(1,:)) + SUM(sfingrd%glamv(2,:)) |
---|
522 | dlgrdv(1) = dlgrdv(1) + SUM(sfingrd%gphiv(1,:)) + SUM(sfingrd%gphiv(2,:)) |
---|
523 | dlgrdv(1) = dlgrdv(1) + SUM(sfingrd%e1v(1,:)) + SUM(sfingrd%e1v(2,:)) |
---|
524 | dlgrdv(1) = dlgrdv(1) + SUM(sfingrd%e2v(1,:)) + SUM(sfingrd%e2v(2,:)) |
---|
525 | ! |
---|
526 | dlgrdv(2) = SUM(sfingrd%glamv(nxfine-1,:)) + SUM(sfingrd%glamv(nxfine,:)) |
---|
527 | dlgrdv(2) = dlgrdv(2) + SUM(sfingrd%gphiv(nxfine-1,:)) + SUM(sfingrd%gphiv(nxfine,:)) |
---|
528 | dlgrdv(2) = dlgrdv(2) + SUM(sfingrd%e1v(nxfine-1,:)) + SUM(sfingrd%e1v(nxfine,:)) |
---|
529 | dlgrdv(2) = dlgrdv(2) + SUM(sfingrd%e2v(nxfine-1,:)) + SUM(sfingrd%e2v(nxfine,:)) |
---|
530 | ! |
---|
531 | IF((dlgrdv(1)-dlgrdv(2)).GT.dleps.OR.(dlgrdv(1)+dlgrdv(2)).LT.dleps) THEN |
---|
532 | WRITE(*,*) ' ERROR' |
---|
533 | print*,(dlgrdv(1)-dlgrdv(2)), (dlgrdv(1)+dlgrdv(2)) |
---|
534 | ELSE |
---|
535 | WRITE(*,*) 'OVERLAP BANDS OK' |
---|
536 | ENDIF |
---|
537 | ! |
---|
538 | WRITE(*,*) '* grid F:' |
---|
539 | dlgrdf(1) = SUM(sfingrd%glamf(1,:)) + SUM(sfingrd%glamf(2,:)) |
---|
540 | dlgrdf(1) = dlgrdf(1) + SUM(sfingrd%gphif(1,:)) + SUM(sfingrd%gphif(2,:)) |
---|
541 | dlgrdf(1) = dlgrdf(1) + SUM(sfingrd%e1f(1,:)) + SUM(sfingrd%e1f(2,:)) |
---|
542 | dlgrdf(1) = dlgrdf(1) + SUM(sfingrd%e2f(1,:)) + SUM(sfingrd%e2f(2,:)) |
---|
543 | ! |
---|
544 | dlgrdf(2) = SUM(sfingrd%glamf(nxfine-1,:)) + SUM(sfingrd%glamf(nxfine,:)) |
---|
545 | dlgrdf(2) = dlgrdf(2) + SUM(sfingrd%gphif(nxfine-1,:)) + SUM(sfingrd%gphif(nxfine,:)) |
---|
546 | dlgrdf(2) = dlgrdf(2) + SUM(sfingrd%e1f(nxfine-1,:)) + SUM(sfingrd%e1f(nxfine,:)) |
---|
547 | dlgrdf(2) = dlgrdf(2) + SUM(sfingrd%e2f(nxfine-1,:)) + SUM(sfingrd%e2f(nxfine,:)) |
---|
548 | ! |
---|
549 | IF((dlgrdf(1)-dlgrdf(2)).GT.dleps.OR.(dlgrdf(1)+dlgrdf(2)).LT.dleps) THEN |
---|
550 | WRITE(*,*) ' ERROR' |
---|
551 | print*, (dlgrdf(1)-dlgrdf(2)), (dlgrdf(1)+dlgrdf(2)) |
---|
552 | ELSE |
---|
553 | WRITE(*,*) 'OVERLAP BANDS OK' |
---|
554 | ENDIF |
---|
555 | ! |
---|
556 | WRITE(*,*) '' |
---|
557 | WRITE(*,*) ' grid T' |
---|
558 | WRITE(*,*) 'i= ',' 1 ',' 2 ',' 3 ' |
---|
559 | WRITE(*,*) sfingrd%glamt(1:3,1) |
---|
560 | WRITE(*,*) sfingrd%gphit(1:3,1) |
---|
561 | WRITE(*,*) 'i= ',' n-2 ',' n-1 ',' n ' |
---|
562 | WRITE(*,*) sfingrd%glamt(nxfine-2:nxfine,1) |
---|
563 | WRITE(*,*) sfingrd%gphit(nxfine-2:nxfine,1) |
---|
564 | WRITE(*,*) '' |
---|
565 | WRITE(*,*) ' grid U' |
---|
566 | WRITE(*,*) 'i= ',' 1 ',' 2 ',' 3 ' |
---|
567 | WRITE(*,*) sfingrd%glamu(1:3,1) |
---|
568 | WRITE(*,*) sfingrd%gphiu(1:3,1) |
---|
569 | WRITE(*,*) 'i= ',' n-2 ',' n-1 ',' n ' |
---|
570 | WRITE(*,*) sfingrd%glamu(nxfine-2:nxfine,1) |
---|
571 | WRITE(*,*) sfingrd%gphiu(nxfine-2:nxfine,1) |
---|
572 | WRITE(*,*) '' |
---|
573 | WRITE(*,*) ' grid V' |
---|
574 | WRITE(*,*) 'i= ',' 1 ',' 2 ',' 3 ' |
---|
575 | WRITE(*,*) sfingrd%glamv(1:3,1) |
---|
576 | WRITE(*,*) sfingrd%gphiv(1:3,1) |
---|
577 | WRITE(*,*) 'i= ',' n-2 ',' n-1 ',' n ' |
---|
578 | WRITE(*,*) sfingrd%glamv(nxfine-2:nxfine,1) |
---|
579 | WRITE(*,*) sfingrd%gphiv(nxfine-2:nxfine,1) |
---|
580 | WRITE(*,*) '' |
---|
581 | WRITE(*,*) ' grid F' |
---|
582 | WRITE(*,*) 'i= ',' 1 ',' 2 ',' 3 ' |
---|
583 | WRITE(*,*) sfingrd%glamf(1:3,1) |
---|
584 | WRITE(*,*) sfingrd%gphif(1:3,1) |
---|
585 | WRITE(*,*) 'i= ',' n-2 ',' n-1 ',' n ' |
---|
586 | WRITE(*,*) sfingrd%glamf(nxfine-2:nxfine,1) |
---|
587 | WRITE(*,*) sfingrd%gphif(nxfine-2:nxfine,1) |
---|
588 | WRITE(*,*) '' |
---|
589 | ENDIF |
---|
590 | ! |
---|
591 | IF(nglobal.AND.nequator.EQ.1) THEN |
---|
592 | ! *** GRID U |
---|
593 | jj = 1 |
---|
594 | DO ji = nxfine/2 + nn_rhox,nxfine,2 |
---|
595 | sfingrd%glamu(ji,nyfine-1) = sfingrd%glamu(nxfine/2+1 - jj,nyfine-1) |
---|
596 | sfingrd%gphiu(ji,nyfine-1) = sfingrd%gphiu(nxfine/2+1 - jj,nyfine-1) |
---|
597 | jj = jj+2 |
---|
598 | ENDDO |
---|
599 | ! |
---|
600 | jj = 0 |
---|
601 | DO ji = 2,nxfine |
---|
602 | sfingrd%glamu(ji,nyfine) = sfingrd%glamu(nxfine - jj,nyfine-2) |
---|
603 | sfingrd%gphiu(ji,nyfine) = sfingrd%gphiu(nxfine - jj,nyfine-2) |
---|
604 | jj = jj+1 |
---|
605 | ENDDO |
---|
606 | ! |
---|
607 | ! *** GRID T |
---|
608 | jj = 0 |
---|
609 | DO ji = nxfine/2 + nn_rhox,nxfine |
---|
610 | sfingrd%glamt(ji,nyfine-1) = sfingrd%glamt(nxfine/2+1 - jj,nyfine-1) |
---|
611 | sfingrd%gphit(ji,nyfine-1) = sfingrd%gphit(nxfine/2+1 - jj,nyfine-1) |
---|
612 | jj = jj+1 |
---|
613 | ENDDO |
---|
614 | ! |
---|
615 | jj = 0 |
---|
616 | DO ji = 3,nxfine |
---|
617 | sfingrd%glamt(ji,nyfine) = sfingrd%glamt(nxfine - jj,nyfine-2) |
---|
618 | sfingrd%gphit(ji,nyfine) = sfingrd%gphit(nxfine - jj,nyfine-2) |
---|
619 | jj = jj+1 |
---|
620 | ENDDO |
---|
621 | ENDIF |
---|
622 | ! |
---|
623 | WRITE(*,*) '' |
---|
624 | WRITE(*,*) '### END SUBROUTINE child_grid ###' |
---|
625 | WRITE(*,*) '' |
---|
626 | ! |
---|
627 | END SUBROUTINE child_grid |
---|
628 | ! |
---|
629 | ! |
---|
630 | ! |
---|
631 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
632 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
633 | END MODULE cfg_tools |
---|