source: CONFIG_DEVT/IPSLCM6.5_work_ENSEMBLES/modeles/NEMO/tools/GRIDGEN/src/cfg_tools.f90 @ 5501

Last change on this file since 5501 was 5501, checked in by aclsce, 4 years ago

First import of IPSLCM6.5_work_ENSEMBLES working configuration

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.