source: trunk/LMDZ.COMMON/libf/dyn3d_common/advz.F

Last change on this file was 1422, checked in by milmd, 10 years ago

In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).

File size: 8.9 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE advz(limit,dtz,w,sm,s0,sx,sy,sz)
5      IMPLICIT NONE
6
7CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
8C                                                                C
9C  first-order moments (FOM) advection of tracer in Z direction  C
10C                                                                C
11C  Source : Pascal Simon (Meteo,CNRM)                            C
12C  Adaptation : A.Armengaud (LGGE) juin 94                       C
13C                                                                C
14C                                                                C
15C  sont des arguments d'entree pour le s-pg...                   C
16C                                                                C
17C  dq est l'argument de sortie pour le s-pg                      C
18C                                                                C
19CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
20C
21C  parametres principaux du modele
22C
23#include "dimensions.h"
24#include "paramet.h"
25
26C    #include "traceur.h"
27
28C  Arguments :
29C  -----------
30C  dtz : frequence fictive d'appel du transport
31C  w : flux de masse en z en Pa.m2.s-1
32
33      INTEGER ntra
34      PARAMETER (ntra = 1)
35
36      REAL dtz
37      REAL w ( iip1,jjp1,llm )
38   
39C  moments: SM  total mass in each grid box
40C           S0  mass of tracer in each grid box
41C           Si  1rst order moment in i direction
42C
43      REAL SM(iip1,jjp1,llm)
44     +    ,S0(iip1,jjp1,llm,ntra)
45      REAL sx(iip1,jjp1,llm,ntra)
46     +    ,sy(iip1,jjp1,llm,ntra)
47     +    ,sz(iip1,jjp1,llm,ntra)
48
49
50C  Local :
51C  -------
52
53C  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
54C  mass fluxes in kg
55C  declaration :
56
57      REAL WGRI(iip1,jjp1,0:llm)
58
59C
60C  the moments F are used as temporary  storage for
61C  portions of grid boxes in transit at the current latitude
62C
63      REAL FM(iim,llm)
64      REAL F0(iim,llm,ntra),FX(iim,llm,ntra)
65      REAL FY(iim,llm,ntra),FZ(iim,llm,ntra)
66C
67C  work arrays
68C
69      REAL ALF(iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)
70      REAL TEMPTM            ! Just temporal variable
71      REAL sqi,sqf
72C
73      LOGICAL LIMIT
74      INTEGER lon,lat,niv
75      INTEGER i,j,jv,k,l,lp
76
77      lon = iim
78      lat = jjp1
79      niv = llm
80
81C *** Test de passage d'arguments ******
82 
83c     DO 399 l = 1, llm
84c     DO 399 j = 1, jjp1
85c     DO 399 i = 1, iip1
86c        IF (S0(i,j,l,ntra) .lt. 0. ) THEN
87c           PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
88c           print*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)
89c           print*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)
90c           print*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)
91c           PRINT*, 'AIE !! debut ADVZ - pbl arg. passage dans ADVZ'
92c            STOP
93c        ENDIF
94  399 CONTINUE
95
96C-----------------------------------------------------------------
97C *** Test : diag de la qqtite totale de traceur
98C            dans l'atmosphere avant l'advection en z
99      sqi = 0.
100      sqf = 0.
101
102      DO l = 1,llm
103         DO j = 1,jjp1
104            DO i = 1,iim
105cIM 240305            sqi = sqi + S0(i,j,l,9)
106               sqi = sqi + S0(i,j,l,ntra)
107            ENDDO
108         ENDDO
109      ENDDO
110      PRINT*,'-------- DIAG DANS ADVZ - ENTREE ---------'
111      PRINT*,'sqi=',sqi
112
113C-----------------------------------------------------------------
114C  Interface : adaptation nouveau modele
115C  -------------------------------------
116C
117C  Conversion du flux de masse en kg.s-1
118
119      DO 500 l = 1,llm
120         DO 500 j = 1,jjp1
121            DO 500 i = 1,iip1 
122c            wgri (i,j,llm+1-l) =  w (i,j,l) / g
123               wgri (i,j,llm+1-l) =  w (i,j,l)
124c             wgri (i,j,0) = 0.                ! a detruire ult.
125c             wgri (i,j,l) = 0.1               !    w (i,j,l)
126c             wgri (i,j,llm) = 0.              ! a detruire ult.
127  500 CONTINUE
128         DO  j = 1,jjp1
129            DO i = 1,iip1 
130               wgri(i,j,0)=0.
131            enddo
132         enddo
133
134C-----------------------------------------------------------------
135 
136C  start here         
137C  boucle sur les latitudes
138C
139      DO 1 K=1,LAT
140C
141C  place limits on appropriate moments before transport
142C      (if flux-limiting is to be applied)
143C
144      IF(.NOT.LIMIT) GO TO 101
145C
146      DO 10 JV=1,NTRA
147      DO 10 L=1,NIV
148         DO 100 I=1,LON
149            sz(I,K,L,JV)=SIGN(AMIN1(AMAX1(S0(I,K,L,JV),0.),
150     +                              ABS(sz(I,K,L,JV))),sz(I,K,L,JV))
151 100     CONTINUE
152 10   CONTINUE
153C
154 101  CONTINUE
155C
156C  boucle sur les niveaux intercouches de 1 a NIV-1
157C   (flux nul au sommet L=0 et a la base L=NIV)
158C
159C  calculate flux and moments between adjacent boxes
160C     (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)
161C  1- create temporary moments/masses for partial boxes in transit
162C  2- reajusts moments remaining in the box
163C
164      DO 11 L=1,NIV-1
165      LP=L+1
166C
167      DO 110 I=1,LON
168C
169         IF(WGRI(I,K,L).LT.0.) THEN
170           FM(I,L)=-WGRI(I,K,L)*DTZ
171           ALF(I)=FM(I,L)/SM(I,K,LP)
172           SM(I,K,LP)=SM(I,K,LP)-FM(I,L)
173         ELSE
174           FM(I,L)=WGRI(I,K,L)*DTZ
175           ALF(I)=FM(I,L)/SM(I,K,L)
176           SM(I,K,L)=SM(I,K,L)-FM(I,L)
177         ENDIF
178C
179         ALFQ (I)=ALF(I)*ALF(I)
180         ALF1 (I)=1.-ALF(I)
181         ALF1Q(I)=ALF1(I)*ALF1(I)
182C
183 110  CONTINUE
184C
185      DO 111 JV=1,NTRA
186      DO 1110 I=1,LON
187C
188         IF(WGRI(I,K,L).LT.0.) THEN
189C
190           F0(I,L,JV)=ALF (I)*( S0(I,K,LP,JV)-ALF1(I)*sz(I,K,LP,JV) )
191           FZ(I,L,JV)=ALFQ(I)*sz(I,K,LP,JV)
192           FX(I,L,JV)=ALF (I)*sx(I,K,LP,JV)
193           FY(I,L,JV)=ALF (I)*sy(I,K,LP,JV)
194C
195           S0(I,K,LP,JV)=S0(I,K,LP,JV)-F0(I,L,JV)
196           sz(I,K,LP,JV)=ALF1Q(I)*sz(I,K,LP,JV)
197           sx(I,K,LP,JV)=sx(I,K,LP,JV)-FX(I,L,JV)
198           sy(I,K,LP,JV)=sy(I,K,LP,JV)-FY(I,L,JV)
199C
200         ELSE
201C
202           F0(I,L,JV)=ALF (I)*(S0(I,K,L,JV)+ALF1(I)*sz(I,K,L,JV) )
203           FZ(I,L,JV)=ALFQ(I)*sz(I,K,L,JV)
204           FX(I,L,JV)=ALF (I)*sx(I,K,L,JV)
205           FY(I,L,JV)=ALF (I)*sy(I,K,L,JV)
206C
207           S0(I,K,L,JV)=S0(I,K,L,JV)-F0(I,L,JV)
208           sz(I,K,L,JV)=ALF1Q(I)*sz(I,K,L,JV)
209           sx(I,K,L,JV)=sx(I,K,L,JV)-FX(I,L,JV)
210           sy(I,K,L,JV)=sy(I,K,L,JV)-FY(I,L,JV)
211C
212         ENDIF
213C
214 1110 CONTINUE
215 111  CONTINUE
216C
217 11   CONTINUE
218C
219C  puts the temporary moments Fi into appropriate neighboring boxes
220C
221      DO 12 L=1,NIV-1
222      LP=L+1
223C
224      DO 120 I=1,LON
225C
226         IF(WGRI(I,K,L).LT.0.) THEN
227           SM(I,K,L)=SM(I,K,L)+FM(I,L)
228           ALF(I)=FM(I,L)/SM(I,K,L)
229         ELSE
230           SM(I,K,LP)=SM(I,K,LP)+FM(I,L)
231           ALF(I)=FM(I,L)/SM(I,K,LP)
232         ENDIF
233C
234         ALF1(I)=1.-ALF(I)
235         ALFQ(I)=ALF(I)*ALF(I)
236         ALF1Q(I)=ALF1(I)*ALF1(I)
237C
238 120  CONTINUE
239C
240      DO 121 JV=1,NTRA
241      DO 1210 I=1,LON
242C
243         IF(WGRI(I,K,L).LT.0.) THEN
244C
245           TEMPTM=-ALF(I)*S0(I,K,L,JV)+ALF1(I)*F0(I,L,JV)
246           S0(I,K,L,JV)=S0(I,K,L,JV)+F0(I,L,JV)
247           sz(I,K,L,JV)=ALF(I)*FZ(I,L,JV)+ALF1(I)*sz(I,K,L,JV)+3.*TEMPTM
248           sx(I,K,L,JV)=sx(I,K,L,JV)+FX(I,L,JV)
249           sy(I,K,L,JV)=sy(I,K,L,JV)+FY(I,L,JV)
250C
251         ELSE
252C
253           TEMPTM=ALF(I)*S0(I,K,LP,JV)-ALF1(I)*F0(I,L,JV)
254           S0(I,K,LP,JV)=S0(I,K,LP,JV)+F0(I,L,JV)
255           sz(I,K,LP,JV)=ALF(I)*FZ(I,L,JV)+ALF1(I)*sz(I,K,LP,JV)
256     +                  +3.*TEMPTM
257           sx(I,K,LP,JV)=sx(I,K,LP,JV)+FX(I,L,JV)
258           sy(I,K,LP,JV)=sy(I,K,LP,JV)+FY(I,L,JV)
259C
260         ENDIF
261C
262 1210 CONTINUE
263 121  CONTINUE
264C
265 12   CONTINUE
266C
267C  fin de la boucle principale sur les latitudes
268C
269 1    CONTINUE
270C
271C-------------------------------------------------------------
272C
273C ----------- AA Test en fin de ADVX ------ Controle des S*
274
275c     DO 9999 l = 1, llm
276c     DO 9999 j = 1, jjp1
277c     DO 9999 i = 1, iip1
278c        IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN
279c           PRINT*, '-------------------'
280c           PRINT*, 'En fin de ADVZ'
281c           PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
282c           print*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)
283c           print*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)
284c           print*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)
285c           WRITE (*,*) 'On arrete !! - pbl en fin de ADVZ1'
286c            STOP
287c        ENDIF
288 9999 CONTINUE
289
290C *** ------------------- bouclage cyclique  en X ------------
291     
292c      DO l = 1,llm
293c         DO j = 1,jjp1
294c            SM(iip1,j,l) = SM(1,j,l)
295c            S0(iip1,j,l,ntra) = S0(1,j,l,ntra)
296C            sx(iip1,j,l,ntra) = sx(1,j,l,ntra)
297c            sy(iip1,j,l,ntra) = sy(1,j,l,ntra)
298c            sz(iip1,j,l,ntra) = sz(1,j,l,ntra)
299c         ENDDO
300c      ENDDO
301           
302C-------------------------------------------------------------
303C *** Test : diag de la qqtite totale de traceur
304C            dans l'atmosphere avant l'advection en z
305      DO l = 1,llm
306         DO j = 1,jjp1
307            DO i = 1,iim
308cIM 240305            sqf = sqf + S0(i,j,l,9)
309               sqf = sqf + S0(i,j,l,ntra)
310            ENDDO
311         ENDDO
312      ENDDO
313      PRINT*,'-------- DIAG DANS ADVZ - SORTIE ---------'
314      PRINT*,'sqf=', sqf
315
316C-------------------------------------------------------------
317      RETURN
318      END
319C_______________________________________________________________
320C_______________________________________________________________
Note: See TracBrowser for help on using the repository browser.