source: branches/UKMO/r6232_tracer_advection/NEMOGCM/TOOLS/GRIDGEN/src/cfg_tools.f90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

File size: 22.3 KB
Line 
1MODULE 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   !
162print*, 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 utilisŽs 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 rŽpŽtŽes dans toute la grille.
269   !Il n'est donc nŽcessaire de calculer qu'une fois les 4 coefficients
270   !qui seront utilisŽs 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!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
633END MODULE cfg_tools
Note: See TracBrowser for help on using the repository browser.