source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advzp.f90 @ 5139

Last change on this file since 5139 was 5136, checked in by abarral, 3 months ago

Put comgeom.h, comgeom2.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: 13.5 KB
Line 
1! $Header$
2
3SUBROUTINE ADVZP(LIMIT, DTZ, W, SM, S0, SSX, SY, SZ &
4        , SSXX, SSXY, SSXZ, SYY, SYZ, SZZ, ntra)
5  USE lmdz_comgeom
6
7  IMPLICIT NONE
8
9  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
10  !                                                                 C
11  !  second-order moments (SOM) advection of tracer in Z direction  C
12  !                                                                 C
13  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
14  !                                                                 C
15  !  Source : Pascal Simon ( Meteo, CNRM )                          C
16  !  Adaptation : A.A. (LGGE)                                       C
17  !  Derniere Modif : 19/11/95 LAST                                 C
18  !                                                                 C
19  !  sont les arguments d'entree pour le s-pg                       C
20  !                                                                 C
21  !  argument de sortie du s-pg                                     C
22  !                                                                 C
23  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
24  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
25  !
26  ! Rem : Probleme aux poles il faut reecrire ce cas specifique
27  !    Attention au sens de l'indexation
28  !
29
30  !
31  !  parametres principaux du modele
32  !
33  INCLUDE "dimensions.h"
34  INCLUDE "paramet.h"
35  !
36  !  Arguments :
37  !  ----------
38  !  dty : frequence fictive d'appel du transport
39  !  parbu,pbarv : flux de masse en x et y en Pa.m2.s-1
40  !
41  INTEGER :: lon, lat, niv
42  INTEGER :: i, j, jv, k, kp, l, lp
43  INTEGER :: ntra
44  ! PARAMETER (ntra = 1)
45  !
46  REAL :: dtz
47  REAL :: w (iip1, jjp1, llm)
48  !
49  !  moments: SM  total mass in each grid box
50  !       S0  mass of tracer in each grid box
51  !       Si  1rst order moment in i direction
52  !
53  REAL :: SM(iip1, jjp1, llm) &
54          , S0(iip1, jjp1, llm, ntra)
55  REAL :: SSX(iip1, jjp1, llm, ntra) &
56          , SY(iip1, jjp1, llm, ntra) &
57          , SZ(iip1, jjp1, llm, ntra) &
58          , SSXX(iip1, jjp1, llm, ntra) &
59          , SSXY(iip1, jjp1, llm, ntra) &
60          , SSXZ(iip1, jjp1, llm, ntra) &
61          , SYY(iip1, jjp1, llm, ntra) &
62          , SYZ(iip1, jjp1, llm, ntra) &
63          , SZZ(iip1, jjp1, llm, ntra)
64  !
65  !  Local :
66  !  -------
67  !
68  !  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
69  !  mass fluxes in kg
70  !  declaration :
71  !
72  REAL :: WGRI(iip1, jjp1, 0:llm)
73
74  ! Rem : UGRI et VGRI ne sont pas utilises dans
75  !  cette SUBROUTINE ( advection en z uniquement )
76  !  Rem 2 :le dimensionnement de VGRI depend de celui de pbarv
77  ! attention a celui de WGRI
78  !
79  !  the moments F are similarly defined and used as temporary
80  !  storage for portions of the grid boxes in transit
81  !
82  !  the moments Fij are used as temporary storage for
83  !  portions of the grid boxes in transit at the current level
84  !
85  !  work arrays
86  !
87  !
88  REAL :: F0(iim, llm, ntra), FM(iim, llm)
89  REAL :: FX(iim, llm, ntra), FY(iim, llm, ntra)
90  REAL :: FZ(iim, llm, ntra)
91  REAL :: FXX(iim, llm, ntra), FXY(iim, llm, ntra)
92  REAL :: FXZ(iim, llm, ntra), FYY(iim, llm, ntra)
93  REAL :: FYZ(iim, llm, ntra), FZZ(iim, llm, ntra)
94  REAL :: S00(ntra)
95  REAL :: SM0             ! Just temporal variable
96  !
97  !  work arrays
98  !
99  REAL :: ALF(iim), ALF1(iim)
100  REAL :: ALFQ(iim), ALF1Q(iim)
101  REAL :: ALF2(iim), ALF3(iim)
102  REAL :: ALF4(iim)
103  REAL :: TEMPTM          ! Just temporal variable
104  REAL :: SLPMAX, S1MAX, S1NEW, S2NEW
105  !
106  REAL :: sqi, sqf
107  LOGICAL :: LIMIT
108
109  lon = iim         ! rem : Il est possible qu'un pbl. arrive ici
110  lat = jjp1        ! a cause des dim. differentes entre les
111  niv = llm         !       tab. S et VGRI
112
113  !-----------------------------------------------------------------
114  ! *** Test : diag de la qtite totale de traceur dans
115  ! l'atmosphere avant l'advection en Y
116  !
117  sqi = 0.
118  sqf = 0.
119  !
120  DO l = 1, llm
121    DO j = 1, jjp1
122      DO i = 1, iim
123        sqi = sqi + S0(i, j, l, ntra)
124      END DO
125    END DO
126  END DO
127  PRINT*, '---------- DIAG DANS ADVZP - ENTREE --------'
128  PRINT*, 'sqi=', sqi
129
130  !-----------------------------------------------------------------
131  !  Interface : adaptation nouveau modele
132  !  -------------------------------------
133  !
134  !  Conversion des flux de masses en kg
135
136  DO l = 1, llm
137    DO j = 1, jjp1
138      DO i = 1, iip1
139        wgri (i, j, llm + 1 - l) = w (i, j, l)
140      END DO
141    END DO
142  END DO
143  do j = 1, jjp1
144    do i = 1, iip1
145      wgri(i, j, 0) = 0.
146    enddo
147  enddo
148  !
149  !AA rem : Je ne suis pas sur du signe
150  !AA       Je ne suis pas sur pour le 0:llm
151  !
152  !-----------------------------------------------------------------
153  !---------------------- START HERE -------------------------------
154  !
155  !  boucle sur les latitudes
156  !
157  DO K = 1, LAT
158    !
159    !  place limits on appropriate moments before transport
160    !  (if flux-limiting is to be applied)
161    !
162    IF(.NOT.LIMIT) GO TO 101
163    !
164    DO JV = 1, NTRA
165      DO L = 1, NIV
166        DO I = 1, LON
167          IF(S0(I, K, L, JV)>0.) THEN
168            SLPMAX = S0(I, K, L, JV)
169            S1MAX = 1.5 * SLPMAX
170            S1NEW = AMIN1(S1MAX, AMAX1(-S1MAX, SZ(I, K, L, JV)))
171            S2NEW = AMIN1(2. * SLPMAX - ABS(S1NEW) / 3., &
172                    AMAX1(ABS(S1NEW) - SLPMAX, SZZ(I, K, L, JV)))
173            SZ (I, K, L, JV) = S1NEW
174            SZZ(I, K, L, JV) = S2NEW
175            SSXZ(I, K, L, JV) = AMIN1(SLPMAX, AMAX1(-SLPMAX, SSXZ(I, K, L, JV)))
176            SYZ(I, K, L, JV) = AMIN1(SLPMAX, AMAX1(-SLPMAX, SYZ(I, K, L, JV)))
177          ELSE
178            SZ (I, K, L, JV) = 0.
179            SZZ(I, K, L, JV) = 0.
180            SSXZ(I, K, L, JV) = 0.
181            SYZ(I, K, L, JV) = 0.
182          ENDIF
183        END DO
184      END DO
185    END DO
186    !
187    101   CONTINUE
188    !
189    !  boucle sur les niveaux intercouches de 1 a NIV-1
190    !   (flux nul au sommet L=0 et a la base L=NIV)
191    !
192    !  calculate flux and moments between adjacent boxes
193    ! (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)
194    !  1- create temporary moments/masses for partial boxes in transit
195    !  2- reajusts moments remaining in the box
196    !
197    DO L = 1, NIV - 1
198      LP = L + 1
199      !
200      DO I = 1, LON
201        !
202        IF(WGRI(I, K, L)<0.) THEN
203          FM(I, L) = -WGRI(I, K, L) * DTZ
204          ALF(I) = FM(I, L) / SM(I, K, LP)
205          SM(I, K, LP) = SM(I, K, LP) - FM(I, L)
206        ELSE
207          FM(I, L) = WGRI(I, K, L) * DTZ
208          ALF(I) = FM(I, L) / SM(I, K, L)
209          SM(I, K, L) = SM(I, K, L) - FM(I, L)
210        ENDIF
211        !
212        ALFQ (I) = ALF(I) * ALF(I)
213        ALF1 (I) = 1. - ALF(I)
214        ALF1Q(I) = ALF1(I) * ALF1(I)
215        ALF2 (I) = ALF1(I) - ALF(I)
216        ALF3 (I) = ALF(I) * ALFQ(I)
217        ALF4 (I) = ALF1(I) * ALF1Q(I)
218        !
219      END DO
220      !
221      DO JV = 1, NTRA
222        DO I = 1, LON
223          !
224          IF(WGRI(I, K, L)<0.) THEN
225            !
226            F0 (I, L, JV) = ALF (I) * (S0(I, K, LP, JV) - ALF1(I) * &
227                    (SZ(I, K, LP, JV) - ALF2(I) * SZZ(I, K, LP, JV)))
228            FZ (I, L, JV) = ALFQ(I) * (SZ(I, K, LP, JV) - 3. * ALF1(I) * SZZ(I, K, LP, JV))
229            FZZ(I, L, JV) = ALF3(I) * SZZ(I, K, LP, JV)
230            FXZ(I, L, JV) = ALFQ(I) * SSXZ(I, K, LP, JV)
231            FYZ(I, L, JV) = ALFQ(I) * SYZ(I, K, LP, JV)
232            FX (I, L, JV) = ALF (I) * (SSX(I, K, LP, JV) - ALF1(I) * SSXZ(I, K, LP, JV))
233            FY (I, L, JV) = ALF (I) * (SY(I, K, LP, JV) - ALF1(I) * SYZ(I, K, LP, JV))
234            FXX(I, L, JV) = ALF (I) * SSXX(I, K, LP, JV)
235            FXY(I, L, JV) = ALF (I) * SSXY(I, K, LP, JV)
236            FYY(I, L, JV) = ALF (I) * SYY(I, K, LP, JV)
237            !
238            S0 (I, K, LP, JV) = S0 (I, K, LP, JV) - F0 (I, L, JV)
239            SZ (I, K, LP, JV) = ALF1Q(I) &
240                    * (SZ(I, K, LP, JV) + 3. * ALF(I) * SZZ(I, K, LP, JV))
241            SZZ(I, K, LP, JV) = ALF4 (I) * SZZ(I, K, LP, JV)
242            SSXZ(I, K, LP, JV) = ALF1Q(I) * SSXZ(I, K, LP, JV)
243            SYZ(I, K, LP, JV) = ALF1Q(I) * SYZ(I, K, LP, JV)
244            SSX (I, K, LP, JV) = SSX (I, K, LP, JV) - FX (I, L, JV)
245            SY (I, K, LP, JV) = SY (I, K, LP, JV) - FY (I, L, JV)
246            SSXX(I, K, LP, JV) = SSXX(I, K, LP, JV) - FXX(I, L, JV)
247            SSXY(I, K, LP, JV) = SSXY(I, K, LP, JV) - FXY(I, L, JV)
248            SYY(I, K, LP, JV) = SYY(I, K, LP, JV) - FYY(I, L, JV)
249            !
250          ELSE
251            !
252            F0 (I, L, JV) = ALF (I) * (S0(I, K, L, JV) &
253                    + ALF1(I) * (SZ(I, K, L, JV) + ALF2(I) * SZZ(I, K, L, JV)))
254            FZ (I, L, JV) = ALFQ(I) * (SZ(I, K, L, JV) + 3. * ALF1(I) * SZZ(I, K, L, JV))
255            FZZ(I, L, JV) = ALF3(I) * SZZ(I, K, L, JV)
256            FXZ(I, L, JV) = ALFQ(I) * SSXZ(I, K, L, JV)
257            FYZ(I, L, JV) = ALFQ(I) * SYZ(I, K, L, JV)
258            FX (I, L, JV) = ALF (I) * (SSX(I, K, L, JV) + ALF1(I) * SSXZ(I, K, L, JV))
259            FY (I, L, JV) = ALF (I) * (SY(I, K, L, JV) + ALF1(I) * SYZ(I, K, L, JV))
260            FXX(I, L, JV) = ALF (I) * SSXX(I, K, L, JV)
261            FXY(I, L, JV) = ALF (I) * SSXY(I, K, L, JV)
262            FYY(I, L, JV) = ALF (I) * SYY(I, K, L, JV)
263            !
264            S0 (I, K, L, JV) = S0 (I, K, L, JV) - F0(I, L, JV)
265            SZ (I, K, L, JV) = ALF1Q(I) * (SZ(I, K, L, JV) - 3. * ALF(I) * SZZ(I, K, L, JV))
266            SZZ(I, K, L, JV) = ALF4 (I) * SZZ(I, K, L, JV)
267            SSXZ(I, K, L, JV) = ALF1Q(I) * SSXZ(I, K, L, JV)
268            SYZ(I, K, L, JV) = ALF1Q(I) * SYZ(I, K, L, JV)
269            SSX (I, K, L, JV) = SSX (I, K, L, JV) - FX (I, L, JV)
270            SY (I, K, L, JV) = SY (I, K, L, JV) - FY (I, L, JV)
271            SSXX(I, K, L, JV) = SSXX(I, K, L, JV) - FXX(I, L, JV)
272            SSXY(I, K, L, JV) = SSXY(I, K, L, JV) - FXY(I, L, JV)
273            SYY(I, K, L, JV) = SYY(I, K, L, JV) - FYY(I, L, JV)
274            !
275          ENDIF
276          !
277        END DO
278      END DO
279      !
280    END DO
281    !
282    !  puts the temporary moments Fi into appropriate neighboring boxes
283    !
284    DO L = 1, NIV - 1
285      LP = L + 1
286      !
287      DO I = 1, LON
288        !
289        IF(WGRI(I, K, L)<0.) THEN
290          SM(I, K, L) = SM(I, K, L) + FM(I, L)
291          ALF(I) = FM(I, L) / SM(I, K, L)
292        ELSE
293          SM(I, K, LP) = SM(I, K, LP) + FM(I, L)
294          ALF(I) = FM(I, L) / SM(I, K, LP)
295        ENDIF
296        !
297        ALF1(I) = 1. - ALF(I)
298        ALFQ(I) = ALF(I) * ALF(I)
299        ALF1Q(I) = ALF1(I) * ALF1(I)
300        ALF2(I) = ALF(I) * ALF1(I)
301        ALF3(I) = ALF1(I) - ALF(I)
302        !
303      END DO
304      !
305      DO JV = 1, NTRA
306        DO I = 1, LON
307          !
308          IF(WGRI(I, K, L)<0.) THEN
309            !
310            TEMPTM = -ALF(I) * S0(I, K, L, JV) + ALF1(I) * F0(I, L, JV)
311            S0 (I, K, L, JV) = S0(I, K, L, JV) + F0(I, L, JV)
312            SZZ(I, K, L, JV) = ALFQ(I) * FZZ(I, L, JV) + ALF1Q(I) * SZZ(I, K, L, JV) &
313                    + 5. * (ALF2(I) * (FZ(I, L, JV) - SZ(I, K, L, JV)) + ALF3(I) * TEMPTM)
314            SZ (I, K, L, JV) = ALF (I) * FZ (I, L, JV) + ALF1 (I) * SZ (I, K, L, JV) &
315                    + 3. * TEMPTM
316            SSXZ(I, K, L, JV) = ALF (I) * FXZ(I, L, JV) + ALF1 (I) * SSXZ(I, K, L, JV) &
317                    + 3. * (ALF1(I) * FX (I, L, JV) - ALF  (I) * SSX (I, K, L, JV))
318            SYZ(I, K, L, JV) = ALF (I) * FYZ(I, L, JV) + ALF1 (I) * SYZ(I, K, L, JV) &
319                    + 3. * (ALF1(I) * FY (I, L, JV) - ALF  (I) * SY (I, K, L, JV))
320            SSX (I, K, L, JV) = SSX (I, K, L, JV) + FX (I, L, JV)
321            SY (I, K, L, JV) = SY (I, K, L, JV) + FY (I, L, JV)
322            SSXX(I, K, L, JV) = SSXX(I, K, L, JV) + FXX(I, L, JV)
323            SSXY(I, K, L, JV) = SSXY(I, K, L, JV) + FXY(I, L, JV)
324            SYY(I, K, L, JV) = SYY(I, K, L, JV) + FYY(I, L, JV)
325            !
326          ELSE
327            !
328            TEMPTM = ALF(I) * S0(I, K, LP, JV) - ALF1(I) * F0(I, L, JV)
329            S0 (I, K, LP, JV) = S0(I, K, LP, JV) + F0(I, L, JV)
330            SZZ(I, K, LP, JV) = ALFQ(I) * FZZ(I, L, JV) + ALF1Q(I) * SZZ(I, K, LP, JV) &
331                    + 5. * (ALF2(I) * (SZ(I, K, LP, JV) - FZ(I, L, JV)) - ALF3(I) * TEMPTM)
332            SZ (I, K, LP, JV) = ALF (I) * FZ(I, L, JV) + ALF1(I) * SZ(I, K, LP, JV) &
333                    + 3. * TEMPTM
334            SSXZ(I, K, LP, JV) = ALF(I) * FXZ(I, L, JV) + ALF1(I) * SSXZ(I, K, LP, JV) &
335                    + 3. * (ALF(I) * SSX(I, K, LP, JV) - ALF1(I) * FX(I, L, JV))
336            SYZ(I, K, LP, JV) = ALF(I) * FYZ(I, L, JV) + ALF1(I) * SYZ(I, K, LP, JV) &
337                    + 3. * (ALF(I) * SY(I, K, LP, JV) - ALF1(I) * FY(I, L, JV))
338            SSX (I, K, LP, JV) = SSX (I, K, LP, JV) + FX (I, L, JV)
339            SY (I, K, LP, JV) = SY (I, K, LP, JV) + FY (I, L, JV)
340            SSXX(I, K, LP, JV) = SSXX(I, K, LP, JV) + FXX(I, L, JV)
341            SSXY(I, K, LP, JV) = SSXY(I, K, LP, JV) + FXY(I, L, JV)
342            SYY(I, K, LP, JV) = SYY(I, K, LP, JV) + FYY(I, L, JV)
343            !
344          ENDIF
345          !
346        END DO
347      END DO
348      !
349    END DO
350    !
351    !  fin de la boucle principale sur les latitudes
352    !
353  END DO
354  !
355  DO l = 1, llm
356    DO j = 1, jjp1
357      SM(iip1, j, l) = SM(1, j, l)
358      S0(iip1, j, l, ntra) = S0(1, j, l, ntra)
359      SSX(iip1, j, l, ntra) = SSX(1, j, l, ntra)
360      SY(iip1, j, l, ntra) = SY(1, j, l, ntra)
361      SZ(iip1, j, l, ntra) = SZ(1, j, l, ntra)
362    ENDDO
363  ENDDO
364  ! C-------------------------------------------------------------
365  ! *** Test : diag de la qqtite totale de tarceur
366  ! dans l'atmosphere avant l'advection en z
367  DO l = 1, llm
368    DO j = 1, jjp1
369      DO i = 1, iim
370        sqf = sqf + S0(i, j, l, ntra)
371      ENDDO
372    ENDDO
373  ENDDO
374  PRINT*, '-------- DIAG DANS ADVZ - SORTIE ---------'
375  PRINT*, 'sqf=', sqf
376
377  RETURN
378END SUBROUTINE ADVZP
Note: See TracBrowser for help on using the repository browser.