source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advy.f90 @ 5209

Last change on this file since 5209 was 5159, checked in by abarral, 7 weeks 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: 10.9 KB
Line 
1! $Header$
2
3SUBROUTINE advy(limit, dty, pbarv, sm, s0, sx, sy, sz)
4  USE lmdz_comgeom2
5
6  USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
7  USE lmdz_paramet
8  IMPLICIT NONE
9
10  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
11  !                                                                C
12  !  first-order moments (SOM) advection of tracer in Y direction  C
13  !                                                                C
14  !  Source : Pascal Simon ( Meteo, CNRM )                         C
15  !  Adaptation : A.A. (LGGE)                                      C
16  !  Derniere Modif : 15/12/94 LAST
17  ! C
18  !  sont les arguments d'entree pour le s-pg                      C
19  !                                                                C
20  !  argument de sortie du s-pg                                    C
21  !                                                                C
22  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
23  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
24
25  !  Rem : Probleme aux poles il faut reecrire ce cas specifique
26  !    Attention au sens de l'indexation
27
28  !  parametres principaux du modele
29
30  !
31
32
33
34  !  Arguments :
35  !  ----------
36  !  dty : frequence fictive d'appel du transport
37  !  parbu,pbarv : flux de masse en x et y en Pa.m2.s-1
38
39  INTEGER :: lon, lat, niv
40  INTEGER :: i, j, jv, k, kp, l
41  INTEGER :: ntra
42  PARAMETER (ntra = 1)
43
44  REAL :: dty
45  REAL :: pbarv (iip1, jjm, llm)
46
47  !  moments: SM  total mass in each grid box
48  ! S0  mass of tracer in each grid box
49  ! Si  1rst order moment in i direction
50
51  REAL :: SM(iip1, jjp1, llm) &
52          , S0(iip1, jjp1, llm, ntra)
53  REAL :: sx(iip1, jjp1, llm, ntra) &
54          , sy(iip1, jjp1, llm, ntra) &
55          , sz(iip1, jjp1, llm, ntra)
56
57
58  !  Local :
59  !  -------
60
61  !  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
62  !  mass fluxes in kg
63  !  declaration :
64
65  REAL :: VGRI(iip1, 0:jjp1, llm)
66
67  !  Rem : UGRI et WGRI ne sont pas utilises dans
68  !  cette SUBROUTINE ( advection en y uniquement )
69  !  Rem 2 :le dimensionnement de VGRI depend de celui de pbarv
70
71  !  the moments F are similarly defined and used as temporary
72  !  storage for portions of the grid boxes in transit
73
74  REAL :: F0(iim, 0:jjp1, ntra), FM(iim, 0:jjp1)
75  REAL :: FX(iim, jjm, ntra), FY(iim, jjm, ntra)
76  REAL :: FZ(iim, jjm, ntra)
77  REAL :: S00(ntra)
78  REAL :: SM0             ! Just temporal variable
79
80  !  work arrays
81
82  REAL :: ALF(iim, 0:jjp1), ALF1(iim, 0:jjp1)
83  REAL :: ALFQ(iim, 0:jjp1), ALF1Q(iim, 0:jjp1)
84  REAL :: TEMPTM          ! Just temporal variable
85
86  !  Special pour poles
87
88  REAL :: sbms, sfms, sfzs, sbmn, sfmn, sfzn
89  REAL :: sns0(ntra), snsz(ntra), snsm
90  REAL :: s1v(llm), slatv(llm)
91  REAL :: qy1(iim, llm, ntra), qylat(iim, llm, ntra)
92  REAL :: cx1(llm, ntra), cxLAT(llm, ntra)
93  REAL :: cy1(llm, ntra), cyLAT(llm, ntra)
94  REAL :: z1(iim), zcos(iim), zsin(iim)
95  REAL :: smpn, smps, s0pn, s0ps
96
97  REAL :: sqi, sqf
98  LOGICAL :: LIMIT
99
100  lon = iim         ! rem : Il est possible qu'un pbl. arrive ici
101  lat = jjp1        ! a cause des dim. differentes entre les
102  niv = llm
103
104
105  !  the moments Fi are used as temporary storage for
106  !  portions of the grid boxes in transit at the current level
107
108  !  work arrays
109  !
110
111  DO l = 1, llm
112    DO j = 1, jjm
113      DO i = 1, iip1
114        vgri (i, j, llm + 1 - l) = -1. * pbarv(i, j, l)
115      enddo
116    enddo
117    DO i = 1, iip1
118      vgri(i, 0, l) = 0.
119      vgri(i, jjp1, l) = 0.
120    enddo
121  enddo
122
123  DO L = 1, NIV
124
125    !  place limits on appropriate moments before transport
126    !  (if flux-limiting is to be applied)
127
128    IF(.NOT.LIMIT) GO TO 11
129
130    DO JV = 1, NTRA
131      DO K = 1, LAT
132        DO I = 1, LON
133          sy(I, K, L, JV) = SIGN(AMIN1(AMAX1(S0(I, K, L, JV), 0.), &
134                  ABS(sy(I, K, L, JV))), sy(I, K, L, JV))
135        END DO
136      END DO
137    END DO
138
139    11   CONTINUE
140
141    !  le flux a travers le pole Nord est traite separement
142
143    SM0 = 0.
144    DO JV = 1, NTRA
145      S00(JV) = 0.
146    END DO
147
148    DO I = 1, LON
149
150      IF(VGRI(I, 0, L)<=0.) THEN
151        FM(I, 0) = -VGRI(I, 0, L) * DTY
152        ALF(I, 0) = FM(I, 0) / SM(I, 1, L)
153        SM(I, 1, L) = SM(I, 1, L) - FM(I, 0)
154        SM0 = SM0 + FM(I, 0)
155      ENDIF
156
157      ALFQ(I, 0) = ALF(I, 0) * ALF(I, 0)
158      ALF1(I, 0) = 1. - ALF(I, 0)
159      ALF1Q(I, 0) = ALF1(I, 0) * ALF1(I, 0)
160
161    END DO
162
163    DO JV = 1, NTRA
164      DO I = 1, LON
165
166        IF(VGRI(I, 0, L)<=0.) THEN
167
168          F0(I, 0, JV) = ALF(I, 0) * &
169                  (S0(I, 1, L, JV) - ALF1(I, 0) * sy(I, 1, L, JV))
170
171          S00(JV) = S00(JV) + F0(I, 0, JV)
172          S0(I, 1, L, JV) = S0(I, 1, L, JV) - F0(I, 0, JV)
173          sy(I, 1, L, JV) = ALF1Q(I, 0) * sy(I, 1, L, JV)
174          sx(I, 1, L, JV) = ALF1 (I, 0) * sx(I, 1, L, JV)
175          sz(I, 1, L, JV) = ALF1 (I, 0) * sz(I, 1, L, JV)
176
177        ENDIF
178
179      END DO
180    END DO
181
182    DO I = 1, LON
183      IF(VGRI(I, 0, L)>0.) THEN
184        FM(I, 0) = VGRI(I, 0, L) * DTY
185        ALF(I, 0) = FM(I, 0) / SM0
186      ENDIF
187    END DO
188
189    DO JV = 1, NTRA
190      DO I = 1, LON
191        IF(VGRI(I, 0, L)>0.) THEN
192          F0(I, 0, JV) = ALF(I, 0) * S00(JV)
193        ENDIF
194      END DO
195    END DO
196
197    !  puts the temporary moments Fi into appropriate neighboring boxes
198
199    DO I = 1, LON
200
201      IF(VGRI(I, 0, L)>0.) THEN
202        SM(I, 1, L) = SM(I, 1, L) + FM(I, 0)
203        ALF(I, 0) = FM(I, 0) / SM(I, 1, L)
204      ENDIF
205
206      ALF1(I, 0) = 1. - ALF(I, 0)
207
208    END DO
209
210    DO JV = 1, NTRA
211      DO I = 1, LON
212
213        IF(VGRI(I, 0, L)>0.) THEN
214
215          TEMPTM = ALF(I, 0) * S0(I, 1, L, JV) - ALF1(I, 0) * F0(I, 0, JV)
216          S0(I, 1, L, JV) = S0(I, 1, L, JV) + F0(I, 0, JV)
217          sy(I, 1, L, JV) = ALF1(I, 0) * sy(I, 1, L, JV) + 3. * TEMPTM
218
219        ENDIF
220
221      END DO
222    END DO
223
224    !  calculate flux and moments between adjacent boxes
225    !  1- create temporary moments/masses for partial boxes in transit
226    !  2- reajusts moments remaining in the box
227
228    !  flux from KP to K if V(K).lt.0 and from K to KP if V(K).gt.0
229
230    DO K = 1, LAT - 1
231      KP = K + 1
232      DO I = 1, LON
233
234        IF(VGRI(I, K, L)<0.) THEN
235          FM(I, K) = -VGRI(I, K, L) * DTY
236          ALF(I, K) = FM(I, K) / SM(I, KP, L)
237          SM(I, KP, L) = SM(I, KP, L) - FM(I, K)
238        ELSE
239          FM(I, K) = VGRI(I, K, L) * DTY
240          ALF(I, K) = FM(I, K) / SM(I, K, L)
241          SM(I, K, L) = SM(I, K, L) - FM(I, K)
242        ENDIF
243
244        ALFQ(I, K) = ALF(I, K) * ALF(I, K)
245        ALF1(I, K) = 1. - ALF(I, K)
246        ALF1Q(I, K) = ALF1(I, K) * ALF1(I, K)
247
248      END DO
249    END DO
250
251    DO JV = 1, NTRA
252      DO K = 1, LAT - 1
253        KP = K + 1
254        DO I = 1, LON
255
256          IF(VGRI(I, K, L)<0.) THEN
257
258            F0(I, K, JV) = ALF (I, K) * &
259                    (S0(I, KP, L, JV) - ALF1(I, K) * sy(I, KP, L, JV))
260            FY(I, K, JV) = ALFQ(I, K) * sy(I, KP, L, JV)
261            FX(I, K, JV) = ALF (I, K) * sx(I, KP, L, JV)
262            FZ(I, K, JV) = ALF (I, K) * sz(I, KP, L, JV)
263
264            S0(I, KP, L, JV) = S0(I, KP, L, JV) - F0(I, K, JV)
265            sy(I, KP, L, JV) = ALF1Q(I, K) * sy(I, KP, L, JV)
266            sx(I, KP, L, JV) = sx(I, KP, L, JV) - FX(I, K, JV)
267            sz(I, KP, L, JV) = sz(I, KP, L, JV) - FZ(I, K, JV)
268
269          ELSE
270
271            F0(I, K, JV) = ALF (I, K) * &
272                    (S0(I, K, L, JV) + ALF1(I, K) * sy(I, K, L, JV))
273            FY(I, K, JV) = ALFQ(I, K) * sy(I, K, L, JV)
274            FX(I, K, JV) = ALF(I, K) * sx(I, K, L, JV)
275            FZ(I, K, JV) = ALF(I, K) * sz(I, K, L, JV)
276
277            S0(I, K, L, JV) = S0(I, K, L, JV) - F0(I, K, JV)
278            sy(I, K, L, JV) = ALF1Q(I, K) * sy(I, K, L, JV)
279            sx(I, K, L, JV) = sx(I, K, L, JV) - FX(I, K, JV)
280            sz(I, K, L, JV) = sz(I, K, L, JV) - FZ(I, K, JV)
281
282          ENDIF
283
284        END DO
285      END DO
286    END DO
287
288    !  puts the temporary moments Fi into appropriate neighboring boxes
289
290    DO K = 1, LAT - 1
291      KP = K + 1
292      DO I = 1, LON
293
294        IF(VGRI(I, K, L)<0.) THEN
295          SM(I, K, L) = SM(I, K, L) + FM(I, K)
296          ALF(I, K) = FM(I, K) / SM(I, K, L)
297        ELSE
298          SM(I, KP, L) = SM(I, KP, L) + FM(I, K)
299          ALF(I, K) = FM(I, K) / SM(I, KP, L)
300        ENDIF
301
302        ALF1(I, K) = 1. - ALF(I, K)
303
304      END DO
305    END DO
306
307    DO JV = 1, NTRA
308      DO K = 1, LAT - 1
309        KP = K + 1
310        DO I = 1, LON
311
312          IF(VGRI(I, K, L)<0.) THEN
313
314            TEMPTM = -ALF(I, K) * S0(I, K, L, JV) + ALF1(I, K) * F0(I, K, JV)
315            S0(I, K, L, JV) = S0(I, K, L, JV) + F0(I, K, JV)
316            sy(I, K, L, JV) = ALF(I, K) * FY(I, K, JV) + ALF1(I, K) * sy(I, K, L, JV) &
317                    + 3. * TEMPTM
318            sx(I, K, L, JV) = sx(I, K, L, JV) + FX(I, K, JV)
319            sz(I, K, L, JV) = sz(I, K, L, JV) + FZ(I, K, JV)
320
321          ELSE
322
323            TEMPTM = ALF(I, K) * S0(I, KP, L, JV) - ALF1(I, K) * F0(I, K, JV)
324            S0(I, KP, L, JV) = S0(I, KP, L, JV) + F0(I, K, JV)
325            sy(I, KP, L, JV) = ALF(I, K) * FY(I, K, JV) + ALF1(I, K) * sy(I, KP, L, JV) &
326                    + 3. * TEMPTM
327            sx(I, KP, L, JV) = sx(I, KP, L, JV) + FX(I, K, JV)
328            sz(I, KP, L, JV) = sz(I, KP, L, JV) + FZ(I, K, JV)
329
330          ENDIF
331
332        END DO
333      END DO
334    END DO
335
336    !  traitement special pour le pole Sud (idem pole Nord)
337
338    K = LAT
339
340    SM0 = 0.
341    DO JV = 1, NTRA
342      S00(JV) = 0.
343    END DO
344
345    DO I = 1, LON
346
347      IF(VGRI(I, K, L)>=0.) THEN
348        FM(I, K) = VGRI(I, K, L) * DTY
349        ALF(I, K) = FM(I, K) / SM(I, K, L)
350        SM(I, K, L) = SM(I, K, L) - FM(I, K)
351        SM0 = SM0 + FM(I, K)
352      ENDIF
353
354      ALFQ(I, K) = ALF(I, K) * ALF(I, K)
355      ALF1(I, K) = 1. - ALF(I, K)
356      ALF1Q(I, K) = ALF1(I, K) * ALF1(I, K)
357
358    END DO
359
360    DO JV = 1, NTRA
361      DO I = 1, LON
362
363        IF(VGRI(I, K, L)>=0.) THEN
364          F0 (I, K, JV) = ALF(I, K) * &
365                  (S0(I, K, L, JV) + ALF1(I, K) * sy(I, K, L, JV))
366          S00(JV) = S00(JV) + F0(I, K, JV)
367
368          S0(I, K, L, JV) = S0 (I, K, L, JV) - F0 (I, K, JV)
369          sy(I, K, L, JV) = ALF1Q(I, K) * sy(I, K, L, JV)
370          sx(I, K, L, JV) = ALF1(I, K) * sx(I, K, L, JV)
371          sz(I, K, L, JV) = ALF1(I, K) * sz(I, K, L, JV)
372        ENDIF
373
374      END DO
375    END DO
376
377    DO I = 1, LON
378      IF(VGRI(I, K, L)<0.) THEN
379        FM(I, K) = -VGRI(I, K, L) * DTY
380        ALF(I, K) = FM(I, K) / SM0
381      ENDIF
382    END DO
383
384    DO JV = 1, NTRA
385      DO I = 1, LON
386        IF(VGRI(I, K, L)<0.) THEN
387          F0(I, K, JV) = ALF(I, K) * S00(JV)
388        ENDIF
389      END DO
390    END DO
391
392    !  puts the temporary moments Fi into appropriate neighboring boxes
393
394    DO I = 1, LON
395
396      IF(VGRI(I, K, L)<0.) THEN
397        SM(I, K, L) = SM(I, K, L) + FM(I, K)
398        ALF(I, K) = FM(I, K) / SM(I, K, L)
399      ENDIF
400
401      ALF1(I, K) = 1. - ALF(I, K)
402
403    END DO
404
405    DO JV = 1, NTRA
406      DO I = 1, LON
407
408        IF(VGRI(I, K, L)<0.) THEN
409
410          TEMPTM = -ALF(I, K) * S0(I, K, L, JV) + ALF1(I, K) * F0(I, K, JV)
411          S0(I, K, L, JV) = S0(I, K, L, JV) + F0(I, K, JV)
412          sy(I, K, L, JV) = ALF1(I, K) * sy(I, K, L, JV) + 3. * TEMPTM
413
414        ENDIF
415
416      END DO
417    END DO
418
419  END DO
420
421  RETURN
422END SUBROUTINE advy
423
Note: See TracBrowser for help on using the repository browser.