source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advx.f90 @ 5209

Last change on this file since 5209 was 5159, checked in by abarral, 3 months ago

Put dimensions.h and paramet.h into modules

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.8 KB
Line 
1
2! $Header$
3
4SUBROUTINE  advx(limit,dtx,pbaru,sm,s0, &
5        sx,sy,sz,lati,latf)
6  USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
7  USE lmdz_paramet
8  IMPLICIT NONE
9
10  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
11  !                                                                C
12  !  first-order moments (FOM) advection of tracer in X direction  C
13  !                                                                C
14  !  Source : Pascal Simon (Meteo,CNRM)                            C
15  !  Adaptation : A.Armengaud (LGGE) juin 94                       C
16  !                                                                C
17  !  limit,dtx,pbaru,pbarv,sm,s0,sx,sy,sz                       C
18  !  sont des arguments d'entree pour le s-pg...                   C
19  !                                                                C
20  !  sm,s0,sx,sy,sz                                                C
21  !  sont les arguments de sortie pour le s-pg                     C
22  !                                                                C
23  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
24
25  !  parametres principaux du modele
26  !
27
28
29
30  !  Arguments :
31  !  -----------
32  !  dtx : frequence fictive d'appel du transport
33  !  pbaru, pbarv : flux de masse en x et y en Pa.m2.s-1
34
35   INTEGER :: ntra
36   PARAMETER (ntra = 1)
37
38  ! ATTENTION partout ou on trouve ntra, insertion de boucle
39        ! possible dans l'avenir.
40
41  REAL :: dtx
42  REAL :: pbaru ( iip1,jjp1,llm )
43
44  !  moments: SM  total mass in each grid box
45        ! S0  mass of tracer in each grid box
46        ! Si  1rst order moment in i direction
47
48  REAL :: SM(iip1,jjp1,llm),S0(iip1,jjp1,llm,ntra)
49  REAL :: sx(iip1,jjp1,llm,ntra) &
50        ,sy(iip1,jjp1,llm,ntra)
51  REAL :: sz(iip1,jjp1,llm,ntra)
52
53  !  Local :
54  !  -------
55
56  !  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
57  !  mass fluxes in kg
58  !  declaration :
59
60  REAL :: UGRI(iip1,jjp1,llm)
61
62  !  Rem : VGRI et WGRI ne sont pas utilises dans
63  !  cette SUBROUTINE ( advection en x uniquement )
64
65  !  Ti are the moments for the current latitude and level
66
67  REAL :: TM(iim)
68  REAL :: T0(iim,ntra),TX(iim,ntra)
69  REAL :: TY(iim,ntra),TZ(iim,ntra)
70  REAL :: TEMPTM                ! just a temporary variable
71
72  !  the moments F are similarly defined and used as temporary
73  !  storage for portions of the grid boxes in transit
74
75  REAL :: FM(iim)
76  REAL :: F0(iim,ntra),FX(iim,ntra)
77  REAL :: FY(iim,ntra),FZ(iim,ntra)
78
79  !  work arrays
80
81  REAL :: ALF(iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)
82
83  REAL :: SMNEW(iim),UEXT(iim)
84
85  REAL :: sqi,sqf
86
87  LOGICAL :: LIMIT
88  INTEGER :: NUM(jjp1),LONK,NUMK
89  INTEGER :: lon,lati,latf,niv
90  INTEGER :: i,i2,i3,j,jv,l,k,itrac
91
92  lon = iim
93  niv = llm
94
95  ! *** Test de passage d'arguments ******
96
97
98  !  -------------------------------------
99  DO j = 1,jjp1
100     NUM(j) = 1
101  END DO
102  sqi = 0.
103  sqf = 0.
104
105  DO l = 1,llm
106     DO j = 1,jjp1
107        DO i = 1,iim
108  !IM 240305            sqi = sqi + S0(i,j,l,9)
109           sqi = sqi + S0(i,j,l,ntra)
110        ENDDO
111     ENDDO
112  ENDDO
113  PRINT*,'-------- DIAG DANS ADVX - ENTREE ---------'
114  PRINT*,'sqi=',sqi
115
116
117  !  Interface : adaptation nouveau modele
118  !  -------------------------------------
119
120  !  ---------------------------------------------------------
121  !  Conversion des flux de masses en kg/s
122  !  pbaru est en N/s d'ou :
123  !  ugri est en kg/s
124
125  DO l = 1,llm
126     DO j = 1,jjm+1
127        DO i = 1,iip1
128         ! ugri (i,j,llm+1-l) = pbaru (i,j,l) * ( dsig(l) / g )
129         ugri (i,j,llm+1-l) = pbaru (i,j,l)
130  END DO
131  END DO
132  END DO
133
134
135  !  ---------------------------------------------------------
136  !  ---------------------------------------------------------
137  !  ---------------------------------------------------------
138
139  !  start here
140
141  !  boucle principale sur les niveaux et les latitudes
142
143  DO L=1,NIV
144  DO K=lati,latf
145
146  !  initialisation
147
148  !  program assumes periodic boundaries in X
149
150  DO I=2,LON
151     SMNEW(I)=SM(I,K,L)+(UGRI(I-1,K,L)-UGRI(I,K,L))*DTX
152  END DO
153  SMNEW(1)=SM(1,K,L)+(UGRI(LON,K,L)-UGRI(1,K,L))*DTX
154
155  !  modifications for extended polar zones
156
157  NUMK=NUM(K)
158  LONK=LON/NUMK
159
160  IF(NUMK>1) THEN
161
162  DO I=1,LON
163     TM(I)=0.
164  END DO
165  DO JV=1,NTRA
166  DO I=1,LON
167     T0(I,JV)=0.
168     TX(I,JV)=0.
169     TY(I,JV)=0.
170     TZ(I,JV)=0.
171  END DO
172  END DO
173
174  DO I2=1,NUMK
175
176     DO I=1,LONK
177        I3=(I-1)*NUMK+I2
178        TM(I)=TM(I)+SM(I3,K,L)
179        ALF(I)=SM(I3,K,L)/TM(I)
180        ALF1(I)=1.-ALF(I)
181  END DO
182
183     DO  JV=1,NTRA
184     DO  I=1,LONK
185        I3=(I-1)*NUMK+I2
186        TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I) &
187              *S0(I3,K,L,JV)
188        T0(I,JV)=T0(I,JV)+S0(I3,K,L,JV)
189        TX(I,JV)=ALF(I)  *sx(I3,K,L,JV)+ &
190              ALF1(I)*TX(I,JV) +3.*TEMPTM
191        TY(I,JV)=TY(I,JV)+sy(I3,K,L,JV)
192        TZ(I,JV)=TZ(I,JV)+sz(I3,K,L,JV)
193     ENDDO
194     ENDDO
195
196  END DO
197
198  ELSE
199
200  DO I=1,LON
201     TM(I)=SM(I,K,L)
202  END DO
203  DO JV=1,NTRA
204  DO I=1,LON
205     T0(I,JV)=S0(I,K,L,JV)
206     TX(I,JV)=sx(I,K,L,JV)
207     TY(I,JV)=sy(I,K,L,JV)
208     TZ(I,JV)=sz(I,K,L,JV)
209  END DO
210  END DO
211
212  ENDIF
213
214  DO I=1,LONK
215     UEXT(I)=UGRI(I*NUMK,K,L)
216  END DO
217
218  !  place limits on appropriate moments before transport
219  !  (if flux-limiting is to be applied)
220
221  IF(.NOT.LIMIT) GO TO 13
222
223  DO JV=1,NTRA
224  DO I=1,LONK
225    TX(I,JV)=SIGN(AMIN1(AMAX1(T0(I,JV),0.),ABS(TX(I,JV))),TX(I,JV))
226  END DO
227  END DO
228
229 13   CONTINUE
230
231  !  calculate flux and moments between adjacent boxes
232  !  1- create temporary moments/masses for partial boxes in transit
233  !  2- reajusts moments remaining in the box
234
235  !  flux from IP to I if U(I).lt.0
236
237  DO I=1,LONK-1
238     IF(UEXT(I)<0.) THEN
239       FM(I)=-UEXT(I)*DTX
240       ALF(I)=FM(I)/TM(I+1)
241       TM(I+1)=TM(I+1)-FM(I)
242     ENDIF
243  END DO
244
245  I=LONK
246  IF(UEXT(I)<0.) THEN
247    FM(I)=-UEXT(I)*DTX
248    ALF(I)=FM(I)/TM(1)
249    TM(1)=TM(1)-FM(I)
250  ENDIF
251
252  !  flux from I to IP if U(I).gt.0
253
254  DO I=1,LONK
255     IF(UEXT(I)>=0.) THEN
256       FM(I)=UEXT(I)*DTX
257       ALF(I)=FM(I)/TM(I)
258       TM(I)=TM(I)-FM(I)
259     ENDIF
260  END DO
261
262  DO I=1,LONK
263     ALFQ(I)=ALF(I)*ALF(I)
264     ALF1(I)=1.-ALF(I)
265     ALF1Q(I)=ALF1(I)*ALF1(I)
266  END DO
267
268  DO JV=1,NTRA
269  DO I=1,LONK-1
270
271     IF(UEXT(I)<0.) THEN
272
273       F0(I,JV)=ALF (I)* ( T0(I+1,JV)-ALF1(I)*TX(I+1,JV) )
274       FX(I,JV)=ALFQ(I)*TX(I+1,JV)
275       FY(I,JV)=ALF (I)*TY(I+1,JV)
276       FZ(I,JV)=ALF (I)*TZ(I+1,JV)
277
278       T0(I+1,JV)=T0(I+1,JV)-F0(I,JV)
279       TX(I+1,JV)=ALF1Q(I)*TX(I+1,JV)
280       TY(I+1,JV)=TY(I+1,JV)-FY(I,JV)
281       TZ(I+1,JV)=TZ(I+1,JV)-FZ(I,JV)
282
283     ENDIF
284
285  END DO
286  END DO
287
288  I=LONK
289  IF(UEXT(I)<0.) THEN
290
291    DO JV=1,NTRA
292
293       F0 (I,JV)=ALF (I)* ( T0(1,JV)-ALF1(I)*TX(1,JV) )
294       FX (I,JV)=ALFQ(I)*TX(1,JV)
295       FY (I,JV)=ALF (I)*TY(1,JV)
296       FZ (I,JV)=ALF (I)*TZ(1,JV)
297
298       T0(1,JV)=T0(1,JV)-F0(I,JV)
299       TX(1,JV)=ALF1Q(I)*TX(1,JV)
300       TY(1,JV)=TY(1,JV)-FY(I,JV)
301       TZ(1,JV)=TZ(1,JV)-FZ(I,JV)
302
303  END DO
304
305  ENDIF
306
307  DO JV=1,NTRA
308  DO I=1,LONK
309
310     IF(UEXT(I)>=0.) THEN
311
312       F0(I,JV)=ALF (I)* ( T0(I,JV)+ALF1(I)*TX(I,JV) )
313       FX(I,JV)=ALFQ(I)*TX(I,JV)
314       FY(I,JV)=ALF (I)*TY(I,JV)
315       FZ(I,JV)=ALF (I)*TZ(I,JV)
316
317       T0(I,JV)=T0(I,JV)-F0(I,JV)
318       TX(I,JV)=ALF1Q(I)*TX(I,JV)
319       TY(I,JV)=TY(I,JV)-FY(I,JV)
320       TZ(I,JV)=TZ(I,JV)-FZ(I,JV)
321
322     ENDIF
323
324  END DO
325  END DO
326
327  !  puts the temporary moments Fi into appropriate neighboring boxes
328
329  DO I=1,LONK
330     IF(UEXT(I)<0.) THEN
331       TM(I)=TM(I)+FM(I)
332       ALF(I)=FM(I)/TM(I)
333     ENDIF
334  END DO
335
336  DO I=1,LONK-1
337     IF(UEXT(I)>=0.) THEN
338       TM(I+1)=TM(I+1)+FM(I)
339       ALF(I)=FM(I)/TM(I+1)
340     ENDIF
341  END DO
342
343  I=LONK
344  IF(UEXT(I)>=0.) THEN
345    TM(1)=TM(1)+FM(I)
346    ALF(I)=FM(I)/TM(1)
347  ENDIF
348
349  DO I=1,LONK
350     ALF1(I)=1.-ALF(I)
351  END DO
352
353  DO JV=1,NTRA
354  DO I=1,LONK
355
356     IF(UEXT(I)<0.) THEN
357
358       TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)*F0(I,JV)
359       T0(I,JV)=T0(I,JV)+F0(I,JV)
360       TX(I,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(I,JV)+3.*TEMPTM
361       TY(I,JV)=TY(I,JV)+FY(I,JV)
362       TZ(I,JV)=TZ(I,JV)+FZ(I,JV)
363
364     ENDIF
365
366  END DO
367  END DO
368
369  DO JV=1,NTRA
370  DO I=1,LONK-1
371
372     IF(UEXT(I)>=0.) THEN
373
374       TEMPTM=ALF(I)*T0(I+1,JV)-ALF1(I)*F0(I,JV)
375       T0(I+1,JV)=T0(I+1,JV)+F0(I,JV)
376       TX(I+1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(I+1,JV)+3.*TEMPTM
377       TY(I+1,JV)=TY(I+1,JV)+FY(I,JV)
378       TZ(I+1,JV)=TZ(I+1,JV)+FZ(I,JV)
379
380     ENDIF
381
382  END DO
383  END DO
384
385  I=LONK
386  IF(UEXT(I)>=0.) THEN
387    DO JV=1,NTRA
388       TEMPTM=ALF(I)*T0(1,JV)-ALF1(I)*F0(I,JV)
389       T0(1,JV)=T0(1,JV)+F0(I,JV)
390       TX(1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(1,JV)+3.*TEMPTM
391       TY(1,JV)=TY(1,JV)+FY(I,JV)
392       TZ(1,JV)=TZ(1,JV)+FZ(I,JV)
393  END DO
394  ENDIF
395
396  !  retour aux mailles d'origine (passage des Tij aux Sij)
397
398  IF(NUMK>1) THEN
399
400  DO I2=1,NUMK
401
402     DO I=1,LONK
403
404        I3=I2+(I-1)*NUMK
405        SM(I3,K,L)=SMNEW(I3)
406        ALF(I)=SMNEW(I3)/TM(I)
407        TM(I)=TM(I)-SMNEW(I3)
408
409        ALFQ(I)=ALF(I)*ALF(I)
410        ALF1(I)=1.-ALF(I)
411        ALF1Q(I)=ALF1(I)*ALF1(I)
412
413  END DO
414  END DO
415
416     DO  JV=1,NTRA
417     DO  I=1,LONK
418
419        I3=I2+(I-1)*NUMK
420        S0(I3,K,L,JV)=ALF (I) &
421              * (T0(I,JV)-ALF1(I)*TX(I,JV))
422        sx(I3,K,L,JV)=ALFQ(I)*TX(I,JV)
423        sy(I3,K,L,JV)=ALF (I)*TY(I,JV)
424        sz(I3,K,L,JV)=ALF (I)*TZ(I,JV)
425
426  !   reajusts moments remaining in the box
427
428        T0(I,JV)=T0(I,JV)-S0(I3,K,L,JV)
429        TX(I,JV)=ALF1Q(I)*TX(I,JV)
430        TY(I,JV)=TY(I,JV)-sy(I3,K,L,JV)
431        TZ(I,JV)=TZ(I,JV)-sz(I3,K,L,JV)
432      ENDDO
433      ENDDO
434
435
436  ELSE
437
438  DO I=1,LON
439     SM(I,K,L)=TM(I)
440  END DO
441  DO JV=1,NTRA
442  DO I=1,LON
443     S0(I,K,L,JV)=T0(I,JV)
444     sx(I,K,L,JV)=TX(I,JV)
445     sy(I,K,L,JV)=TY(I,JV)
446     sz(I,K,L,JV)=TZ(I,JV)
447  END DO
448  END DO
449
450  ENDIF
451
452  END DO
453  END DO
454
455  ! ----------- AA Test en fin de ADVX ------ Controle des S*
456  ! OK
457  !  DO 9998 l = 1, llm
458  !  DO 9998 j = 1, jjp1
459  !  DO 9998 i = 1, iip1
460  !     IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN
461  !        PRINT*, '-------------------'
462  !        PRINT*, 'En fin de ADVX'
463  !        PRINT*,'SM(',i,j,l,')=',SM(i,j,l)
464  !        PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
465  !        PRINT*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)
466  !        PRINT*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)
467  !        PRINT*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)
468  !        WRITE (*,*) 'On arrete !! - pbl en fin de ADVX1'
469  !c            STOP
470  !     ENDIF
471  ! 9998 CONTINUE
472
473  ! ---------- bouclage cyclique
474  DO itrac=1,ntra
475  DO l = 1,llm
476    DO j = lati,latf
477       SM(iip1,j,l) = SM(1,j,l)
478       S0(iip1,j,l,itrac) = S0(1,j,l,itrac)
479       sx(iip1,j,l,itrac) = sx(1,j,l,itrac)
480       sy(iip1,j,l,itrac) = sy(1,j,l,itrac)
481       sz(iip1,j,l,itrac) = sz(1,j,l,itrac)
482    END DO
483  END DO
484  ENDDO
485
486  ! ----------- qqtite totale de traceur dans tte l'atmosphere
487  DO l = 1, llm
488    DO j = 1, jjp1
489      DO i = 1, iim
490  !IM 240405          sqf = sqf + S0(i,j,l,9)
491         sqf = sqf + S0(i,j,l,ntra)
492      END DO
493    END DO
494  END DO
495
496  PRINT*,'------ DIAG DANS ADVX - SORTIE -----'
497  PRINT*,'sqf=',sqf
498  !-------------
499
500  RETURN
501END SUBROUTINE advx
502!_________________________________________________________________
503!_________________________________________________________________
Note: See TracBrowser for help on using the repository browser.