source: trunk/libf/dyn3d/advz.F @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago

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