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

Last change on this file since 5226 was 5159, checked in by abarral, 16 months ago

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