source: LMDZ6/trunk/libf/dyn3d_common/advx.f90 @ 5272

Last change on this file since 5272 was 5272, checked in by abarral, 23 hours ago

Turn paramet.h into a module

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