source: LMDZ5/branches/testing/libf/dyn3dmem/advz.F @ 1695

Last change on this file since 1695 was 1669, checked in by Laurent Fairhead, 12 years ago

Version testing basée sur la r1668

http://lmdz.lmd.jussieu.fr/utilisateurs/distribution-du-modele/versions-intermediaires

Testing release based on r1668
File size: 8.9 KB
RevLine 
[1632]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
107               sqi = sqi + S0(i,j,l,ntra)
108            ENDDO
109         ENDDO
110      ENDDO
111      PRINT*,'-------- DIAG DANS ADVZ - ENTREE ---------'
112      PRINT*,'sqi=',sqi
113
114C-----------------------------------------------------------------
115C  Interface : adaptation nouveau modele
116C  -------------------------------------
117C
118C  Conversion du flux de masse en kg.s-1
119
120      DO 500 l = 1,llm
121         DO 500 j = 1,jjp1
122            DO 500 i = 1,iip1 
123c            wgri (i,j,llm+1-l) =  w (i,j,l) / g
124               wgri (i,j,llm+1-l) =  w (i,j,l)
125c             wgri (i,j,0) = 0.                ! a detruire ult.
126c             wgri (i,j,l) = 0.1               !    w (i,j,l)
127c             wgri (i,j,llm) = 0.              ! a detruire ult.
128  500 CONTINUE
129         DO  j = 1,jjp1
130            DO i = 1,iip1 
131               wgri(i,j,0)=0.
132            enddo
133         enddo
134
135C-----------------------------------------------------------------
136 
137C  start here         
138C  boucle sur les latitudes
139C
140      DO 1 K=1,LAT
141C
142C  place limits on appropriate moments before transport
143C      (if flux-limiting is to be applied)
144C
145      IF(.NOT.LIMIT) GO TO 101
146C
147      DO 10 JV=1,NTRA
148      DO 10 L=1,NIV
149         DO 100 I=1,LON
150            sz(I,K,L,JV)=SIGN(AMIN1(AMAX1(S0(I,K,L,JV),0.),
151     +                              ABS(sz(I,K,L,JV))),sz(I,K,L,JV))
152 100     CONTINUE
153 10   CONTINUE
154C
155 101  CONTINUE
156C
157C  boucle sur les niveaux intercouches de 1 a NIV-1
158C   (flux nul au sommet L=0 et a la base L=NIV)
159C
160C  calculate flux and moments between adjacent boxes
161C     (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)
162C  1- create temporary moments/masses for partial boxes in transit
163C  2- reajusts moments remaining in the box
164C
165      DO 11 L=1,NIV-1
166      LP=L+1
167C
168      DO 110 I=1,LON
169C
170         IF(WGRI(I,K,L).LT.0.) THEN
171           FM(I,L)=-WGRI(I,K,L)*DTZ
172           ALF(I)=FM(I,L)/SM(I,K,LP)
173           SM(I,K,LP)=SM(I,K,LP)-FM(I,L)
174         ELSE
175           FM(I,L)=WGRI(I,K,L)*DTZ
176           ALF(I)=FM(I,L)/SM(I,K,L)
177           SM(I,K,L)=SM(I,K,L)-FM(I,L)
178         ENDIF
179C
180         ALFQ (I)=ALF(I)*ALF(I)
181         ALF1 (I)=1.-ALF(I)
182         ALF1Q(I)=ALF1(I)*ALF1(I)
183C
184 110  CONTINUE
185C
186      DO 111 JV=1,NTRA
187      DO 1110 I=1,LON
188C
189         IF(WGRI(I,K,L).LT.0.) THEN
190C
191           F0(I,L,JV)=ALF (I)*( S0(I,K,LP,JV)-ALF1(I)*sz(I,K,LP,JV) )
192           FZ(I,L,JV)=ALFQ(I)*sz(I,K,LP,JV)
193           FX(I,L,JV)=ALF (I)*sx(I,K,LP,JV)
194           FY(I,L,JV)=ALF (I)*sy(I,K,LP,JV)
195C
196           S0(I,K,LP,JV)=S0(I,K,LP,JV)-F0(I,L,JV)
197           sz(I,K,LP,JV)=ALF1Q(I)*sz(I,K,LP,JV)
198           sx(I,K,LP,JV)=sx(I,K,LP,JV)-FX(I,L,JV)
199           sy(I,K,LP,JV)=sy(I,K,LP,JV)-FY(I,L,JV)
200C
201         ELSE
202C
203           F0(I,L,JV)=ALF (I)*(S0(I,K,L,JV)+ALF1(I)*sz(I,K,L,JV) )
204           FZ(I,L,JV)=ALFQ(I)*sz(I,K,L,JV)
205           FX(I,L,JV)=ALF (I)*sx(I,K,L,JV)
206           FY(I,L,JV)=ALF (I)*sy(I,K,L,JV)
207C
208           S0(I,K,L,JV)=S0(I,K,L,JV)-F0(I,L,JV)
209           sz(I,K,L,JV)=ALF1Q(I)*sz(I,K,L,JV)
210           sx(I,K,L,JV)=sx(I,K,L,JV)-FX(I,L,JV)
211           sy(I,K,L,JV)=sy(I,K,L,JV)-FY(I,L,JV)
212C
213         ENDIF
214C
215 1110 CONTINUE
216 111  CONTINUE
217C
218 11   CONTINUE
219C
220C  puts the temporary moments Fi into appropriate neighboring boxes
221C
222      DO 12 L=1,NIV-1
223      LP=L+1
224C
225      DO 120 I=1,LON
226C
227         IF(WGRI(I,K,L).LT.0.) THEN
228           SM(I,K,L)=SM(I,K,L)+FM(I,L)
229           ALF(I)=FM(I,L)/SM(I,K,L)
230         ELSE
231           SM(I,K,LP)=SM(I,K,LP)+FM(I,L)
232           ALF(I)=FM(I,L)/SM(I,K,LP)
233         ENDIF
234C
235         ALF1(I)=1.-ALF(I)
236         ALFQ(I)=ALF(I)*ALF(I)
237         ALF1Q(I)=ALF1(I)*ALF1(I)
238C
239 120  CONTINUE
240C
241      DO 121 JV=1,NTRA
242      DO 1210 I=1,LON
243C
244         IF(WGRI(I,K,L).LT.0.) THEN
245C
246           TEMPTM=-ALF(I)*S0(I,K,L,JV)+ALF1(I)*F0(I,L,JV)
247           S0(I,K,L,JV)=S0(I,K,L,JV)+F0(I,L,JV)
248           sz(I,K,L,JV)=ALF(I)*FZ(I,L,JV)+ALF1(I)*sz(I,K,L,JV)+3.*TEMPTM
249           sx(I,K,L,JV)=sx(I,K,L,JV)+FX(I,L,JV)
250           sy(I,K,L,JV)=sy(I,K,L,JV)+FY(I,L,JV)
251C
252         ELSE
253C
254           TEMPTM=ALF(I)*S0(I,K,LP,JV)-ALF1(I)*F0(I,L,JV)
255           S0(I,K,LP,JV)=S0(I,K,LP,JV)+F0(I,L,JV)
256           sz(I,K,LP,JV)=ALF(I)*FZ(I,L,JV)+ALF1(I)*sz(I,K,LP,JV)
257     +                  +3.*TEMPTM
258           sx(I,K,LP,JV)=sx(I,K,LP,JV)+FX(I,L,JV)
259           sy(I,K,LP,JV)=sy(I,K,LP,JV)+FY(I,L,JV)
260C
261         ENDIF
262C
263 1210 CONTINUE
264 121  CONTINUE
265C
266 12   CONTINUE
267C
268C  fin de la boucle principale sur les latitudes
269C
270 1    CONTINUE
271C
272C-------------------------------------------------------------
273C
274C ----------- AA Test en fin de ADVX ------ Controle des S*
275
276c     DO 9999 l = 1, llm
277c     DO 9999 j = 1, jjp1
278c     DO 9999 i = 1, iip1
279c        IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN
280c           PRINT*, '-------------------'
281c           PRINT*, 'En fin de ADVZ'
282c           PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
283c           print*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)
284c           print*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)
285c           print*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)
286c           WRITE (*,*) 'On arrete !! - pbl en fin de ADVZ1'
287c            STOP
288c        ENDIF
289 9999 CONTINUE
290
291C *** ------------------- bouclage cyclique  en X ------------
292     
293c      DO l = 1,llm
294c         DO j = 1,jjp1
295c            SM(iip1,j,l) = SM(1,j,l)
296c            S0(iip1,j,l,ntra) = S0(1,j,l,ntra)
297C            sx(iip1,j,l,ntra) = sx(1,j,l,ntra)
298c            sy(iip1,j,l,ntra) = sy(1,j,l,ntra)
299c            sz(iip1,j,l,ntra) = sz(1,j,l,ntra)
300c         ENDDO
301c      ENDDO
302           
303C-------------------------------------------------------------
304C *** Test : diag de la qqtite totale de traceur
305C            dans l'atmosphere avant l'advection en z
306      DO l = 1,llm
307         DO j = 1,jjp1
308            DO i = 1,iim
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.