source: LMDZ6/branches/Amaury_dev/libf/dyn3d/vlsplt.F90 @ 5137

Last change on this file since 5137 was 5136, checked in by abarral, 11 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: 25.0 KB
Line 
1!
2! $Id: vlsplt.F90 5136 2024-07-28 14:17:54Z abarral $
3!
4
5SUBROUTINE vlsplt(q, pente_max, masse, w, pbaru, pbarv, pdt, iq)
6  USE infotrac, ONLY: nqtot, tracers
7  USE lmdz_ssum_scopy, ONLY: scopy
8  !
9  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
10  !
11  !    ********************************************************************
12  ! Shema  d'advection " pseudo amont " .
13  !    ********************************************************************
14  ! q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
15  !
16  !   pente_max facteur de limitation des pentes: 2 en general
17  !                                           0 pour un schema amont
18  !   pbaru,pbarv,w flux de masse en u ,v ,w
19  !   pdt pas de temps
20  !
21  !   --------------------------------------------------------------------
22  IMPLICIT NONE
23  !
24  INCLUDE "dimensions.h"
25  INCLUDE "paramet.h"
26
27  !
28  !   Arguments:
29  !   ----------
30  REAL :: masse(ip1jmp1, llm), pente_max
31  REAL :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
32  REAL :: q(ip1jmp1, llm, nqtot)
33  REAL :: w(ip1jmp1, llm), pdt
34  INTEGER :: iq ! CRisi
35  !
36  !  Local
37  !   ---------
38  !
39  INTEGER :: ij, l
40  !
41  REAL :: zm(ip1jmp1, llm, nqtot)
42  REAL :: mu(ip1jmp1, llm)
43  REAL :: mv(ip1jm, llm)
44  REAL :: mw(ip1jmp1, llm + 1)
45  REAL :: zq(ip1jmp1, llm, nqtot)
46  REAL :: zzpbar, zzw
47  INTEGER :: ifils, iq2 ! CRisi
48
49  REAL :: qmin, qmax
50  DATA qmin, qmax/0., 1.e33/
51
52  zzpbar = 0.5 * pdt
53  zzw = pdt
54  DO l = 1, llm
55    DO ij = iip2, ip1jm
56      mu(ij, l) = pbaru(ij, l) * zzpbar
57    ENDDO
58    DO ij = 1, ip1jm
59      mv(ij, l) = pbarv(ij, l) * zzpbar
60    ENDDO
61    DO ij = 1, ip1jmp1
62      mw(ij, l) = w(ij, l) * zzw
63    ENDDO
64  ENDDO
65
66  DO ij = 1, ip1jmp1
67    mw(ij, llm + 1) = 0.
68  ENDDO
69
70  CALL SCOPY(ijp1llm, q(1, 1, iq), 1, zq(1, 1, iq), 1)
71  CALL SCOPY(ijp1llm, masse, 1, zm(1, 1, iq), 1)
72
73  do ifils = 1, tracers(iq)%nqDescen
74    iq2 = tracers(iq)%iqDescen(ifils)
75    CALL SCOPY(ijp1llm, q(1, 1, iq2), 1, zq(1, 1, iq2), 1)
76  enddo
77
78  CALL vlx(zq, pente_max, zm, mu, iq)
79  CALL vly(zq, pente_max, zm, mv, iq)
80  CALL vlz(zq, pente_max, zm, mw, iq)
81  CALL vly(zq, pente_max, zm, mv, iq)
82  CALL vlx(zq, pente_max, zm, mu, iq)
83
84  DO l = 1, llm
85    DO ij = 1, ip1jmp1
86      q(ij, l, iq) = zq(ij, l, iq)
87    ENDDO
88    DO ij = 1, ip1jm + 1, iip1
89      q(ij + iim, l, iq) = q(ij, l, iq)
90    ENDDO
91  ENDDO
92  ! CRisi: aussi pour les fils
93  do ifils = 1, tracers(iq)%nqDescen
94    iq2 = tracers(iq)%iqDescen(ifils)
95    DO l = 1, llm
96      DO ij = 1, ip1jmp1
97        q(ij, l, iq2) = zq(ij, l, iq2)
98      ENDDO
99      DO ij = 1, ip1jm + 1, iip1
100        q(ij + iim, l, iq2) = q(ij, l, iq2)
101      ENDDO
102    ENDDO
103  enddo
104
105END SUBROUTINE vlsplt
106RECURSIVE SUBROUTINE vlx(q, pente_max, masse, u_m, iq)
107  USE infotrac, ONLY: nqtot, tracers, & ! CRisi
108          min_qParent, min_qMass, min_ratio ! MVals et CRisi
109  USE lmdz_iniprint, ONLY: lunout, prt_level
110
111  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
112  !
113  !    ********************************************************************
114  ! Shema  d'advection " pseudo amont " .
115  !    ********************************************************************
116  ! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
117  !
118  !
119  !   --------------------------------------------------------------------
120  IMPLICIT NONE
121  !
122  INCLUDE "dimensions.h"
123  INCLUDE "paramet.h"
124  !
125  !
126  !   Arguments:
127  !   ----------
128  REAL :: masse(ip1jmp1, llm, nqtot), pente_max
129  REAL :: u_m(ip1jmp1, llm)
130  REAL :: q(ip1jmp1, llm, nqtot)
131  INTEGER :: iq ! CRisi
132  !
133  !  Local
134  !   ---------
135  !
136  INTEGER :: ij, l, j, i, iju, ijq, indu(ip1jmp1), niju
137  INTEGER :: n0, iadvplus(ip1jmp1, llm), nl(llm)
138  !
139  REAL :: new_m, zu_m, zdum(ip1jmp1, llm)
140  REAL :: dxq(ip1jmp1, llm), dxqu(ip1jmp1)
141  REAL :: zz(ip1jmp1)
142  REAL :: adxqu(ip1jmp1), dxqmax(ip1jmp1, llm)
143  REAL :: u_mq(ip1jmp1, llm)
144
145  ! CRisi
146  REAL :: masseq(ip1jmp1, llm, nqtot), Ratio(ip1jmp1, llm, nqtot)
147  INTEGER :: ifils, iq2 ! CRisi
148
149  LOGICAL, SAVE :: first
150  DATA first/.TRUE./
151
152  !   calcul de la pente a droite et a gauche de la maille
153
154  IF (pente_max>-1.e-5) THEN
155    ! IF (pente_max.gt.10) THEN
156
157    !   calcul des pentes avec limitation, Van Leer scheme I:
158    !   -----------------------------------------------------
159
160    !   calcul de la pente aux points u
161    DO l = 1, llm
162      DO ij = iip2, ip1jm - 1
163        dxqu(ij) = q(ij + 1, l, iq) - q(ij, l, iq)
164      ENDDO
165      DO ij = iip1 + iip1, ip1jm, iip1
166        dxqu(ij) = dxqu(ij - iim)
167        ! sigu(ij)=sigu(ij-iim)
168      ENDDO
169
170      DO ij = iip2, ip1jm
171        adxqu(ij) = abs(dxqu(ij))
172      ENDDO
173
174      !   calcul de la pente maximum dans la maille en valeur absolue
175
176      DO ij = iip2 + 1, ip1jm
177        dxqmax(ij, l) = pente_max * &
178                min(adxqu(ij - 1), adxqu(ij))
179        ! limitation subtile
180        !    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
181
182      ENDDO
183
184      DO ij = iip1 + iip1, ip1jm, iip1
185        dxqmax(ij - iim, l) = dxqmax(ij, l)
186      ENDDO
187
188      DO ij = iip2 + 1, ip1jm
189        IF(dxqu(ij - 1) * dxqu(ij)>0) THEN
190          dxq(ij, l) = dxqu(ij - 1) + dxqu(ij)
191        ELSE
192          !   extremum local
193          dxq(ij, l) = 0.
194        ENDIF
195        dxq(ij, l) = 0.5 * dxq(ij, l)
196        dxq(ij, l) = &
197                sign(min(abs(dxq(ij, l)), dxqmax(ij, l)), dxq(ij, l))
198      ENDDO
199
200    ENDDO ! l=1,llm
201    !print*,'Ok calcul des pentes'
202
203  ELSE ! (pente_max.lt.-1.e-5)
204
205    !   Pentes produits:
206    !   ----------------
207
208    DO l = 1, llm
209      DO ij = iip2, ip1jm - 1
210        dxqu(ij) = q(ij + 1, l, iq) - q(ij, l, iq)
211      ENDDO
212      DO ij = iip1 + iip1, ip1jm, iip1
213        dxqu(ij) = dxqu(ij - iim)
214      ENDDO
215
216      DO ij = iip2 + 1, ip1jm
217        zz(ij) = dxqu(ij - 1) * dxqu(ij)
218        zz(ij) = zz(ij) + zz(ij)
219        IF(zz(ij)>0) THEN
220          dxq(ij, l) = zz(ij) / (dxqu(ij - 1) + dxqu(ij))
221        ELSE
222          !   extremum local
223          dxq(ij, l) = 0.
224        ENDIF
225      ENDDO
226
227    ENDDO
228
229  ENDIF ! (pente_max.lt.-1.e-5)
230
231  !   bouclage de la pente en iip1:
232  !   -----------------------------
233
234  DO l = 1, llm
235    DO ij = iip1 + iip1, ip1jm, iip1
236      dxq(ij - iim, l) = dxq(ij, l)
237    ENDDO
238    DO ij = 1, ip1jmp1
239      iadvplus(ij, l) = 0
240    ENDDO
241
242  ENDDO
243  !   calcul des flux a gauche et a droite
244
245  !   on cumule le flux correspondant a toutes les mailles dont la masse
246  !   au travers de la paroi pENDant le pas de temps.
247  DO l = 1, llm
248    DO ij = iip2, ip1jm - 1
249      IF (u_m(ij, l)>0.) THEN
250        zdum(ij, l) = 1. - u_m(ij, l) / masse(ij, l, iq)
251        u_mq(ij, l) = u_m(ij, l) * (q(ij, l, iq) + 0.5 * zdum(ij, l) * dxq(ij, l))
252      ELSE
253        zdum(ij, l) = 1. + u_m(ij, l) / masse(ij + 1, l, iq)
254        u_mq(ij, l) = u_m(ij, l) * (q(ij + 1, l, iq) &
255                - 0.5 * zdum(ij, l) * dxq(ij + 1, l))
256      ENDIF
257    ENDDO
258  ENDDO
259
260  !   detection des points ou on advecte plus que la masse de la
261  !   maille
262  DO l = 1, llm
263    DO ij = iip2, ip1jm - 1
264      IF(zdum(ij, l)<0) THEN
265        iadvplus(ij, l) = 1
266        u_mq(ij, l) = 0.
267      ENDIF
268    ENDDO
269  ENDDO
270  DO l = 1, llm
271    DO ij = iip1 + iip1, ip1jm, iip1
272      iadvplus(ij, l) = iadvplus(ij - iim, l)
273    ENDDO
274  ENDDO
275
276
277  !   traitement special pour le cas ou on advecte en longitude plus que le
278  !   contenu de la maille.
279  !   cette partie est mal vectorisee.
280
281  !  calcul du nombre de maille sur lequel on advecte plus que la maille.
282
283  n0 = 0
284  DO l = 1, llm
285    nl(l) = 0
286    DO ij = iip2, ip1jm
287      nl(l) = nl(l) + iadvplus(ij, l)
288    ENDDO
289    n0 = n0 + nl(l)
290  ENDDO
291
292  IF(n0>0) THEN
293    IF (prt_level > 2) PRINT *, &
294            'Nombre de points pour lesquels on advect plus que le' &
295            , 'contenu de la maille : ', n0
296
297    DO l = 1, llm
298      IF(nl(l)>0) THEN
299        iju = 0
300        !   indicage des mailles concernees par le traitement special
301        DO ij = iip2, ip1jm
302          IF(iadvplus(ij, l)==1.AND.mod(ij, iip1)/=0) THEN
303            iju = iju + 1
304            indu(iju) = ij
305          ENDIF
306        ENDDO
307        niju = iju
308
309        !  traitement des mailles
310        DO iju = 1, niju
311          ij = indu(iju)
312          j = (ij - 1) / iip1 + 1
313          zu_m = u_m(ij, l)
314          u_mq(ij, l) = 0.
315          IF(zu_m>0.) THEN
316            ijq = ij
317            i = ijq - (j - 1) * iip1
318            !   accumulation pour les mailles completements advectees
319            do while(zu_m>masse(ijq, l, iq))
320              u_mq(ij, l) = u_mq(ij, l) + q(ijq, l, iq) &
321                      * masse(ijq, l, iq)
322              zu_m = zu_m - masse(ijq, l, iq)
323              i = mod(i - 2 + iim, iim) + 1
324              ijq = (j - 1) * iip1 + i
325            ENDDO
326            !   ajout de la maille non completement advectee
327            u_mq(ij, l) = u_mq(ij, l) + zu_m * &
328                    (q(ijq, l, iq) + 0.5 * (1. - zu_m / masse(ijq, l, iq)) &
329                            * dxq(ijq, l))
330          ELSE
331            ijq = ij + 1
332            i = ijq - (j - 1) * iip1
333            !   accumulation pour les mailles completements advectees
334            do while(-zu_m>masse(ijq, l, iq))
335              u_mq(ij, l) = u_mq(ij, l) - q(ijq, l, iq) &
336                      * masse(ijq, l, iq)
337              zu_m = zu_m + masse(ijq, l, iq)
338              i = mod(i, iim) + 1
339              ijq = (j - 1) * iip1 + i
340            ENDDO
341            !   ajout de la maille non completement advectee
342            u_mq(ij, l) = u_mq(ij, l) + zu_m * (q(ijq, l, iq) - &
343                    0.5 * (1. + zu_m / masse(ijq, l, iq)) * dxq(ijq, l))
344          ENDIF
345        ENDDO
346      ENDIF
347    ENDDO
348  ENDIF  ! n0.gt.0
349
350
351  !   bouclage en latitude
352  !print*,'cvant bouclage en latitude'
353  DO l = 1, llm
354    DO ij = iip1 + iip1, ip1jm, iip1
355      u_mq(ij, l) = u_mq(ij - iim, l)
356    ENDDO
357  ENDDO
358
359  ! CRisi: appel récursif de l'advection sur les fils.
360  ! Il faut faire ça avant d'avoir mis à jour q et masse
361  !WRITE(*,*) 'vlsplt 326: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen
362
363  do ifils = 1, tracers(iq)%nqDescen
364    iq2 = tracers(iq)%iqDescen(ifils)
365    DO l = 1, llm
366      DO ij = iip2, ip1jm
367        ! On a besoin de q et masse seulement entre iip2 et ip1jm
368        !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
369        !           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
370        !Mvals: veiller a ce qu'on n'ait pas de denominateur nul
371        masseq(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass)
372        IF (q(ij, l, iq)>min_qParent) THEN
373          Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq)
374        else
375          Ratio(ij, l, iq2) = min_ratio
376        endif
377      enddo
378    enddo
379  enddo
380  do ifils = 1, tracers(iq)%nqChildren
381    iq2 = tracers(iq)%iqDescen(ifils)
382    CALL vlx(Ratio, pente_max, masseq, u_mq, iq2)
383  enddo
384  ! end CRisi
385
386
387  !   calcul des tENDances
388
389  DO l = 1, llm
390    DO ij = iip2 + 1, ip1jm
391      !MVals: veiller a ce qu'on ait pas de denominateur nul
392      new_m = max(masse(ij, l, iq) + u_m(ij - 1, l) - u_m(ij, l), min_qMass)
393      q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) + &
394              u_mq(ij - 1, l) - u_mq(ij, l)) &
395              / new_m
396      masse(ij, l, iq) = new_m
397    ENDDO
398    DO ij = iip1 + iip1, ip1jm, iip1
399      q(ij - iim, l, iq) = q(ij, l, iq)
400      masse(ij - iim, l, iq) = masse(ij, l, iq)
401    ENDDO
402  ENDDO
403
404  ! retablir les fils en rapport de melange par rapport a l'air:
405  ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio
406  ! puis on boucle en longitude
407  do ifils = 1, tracers(iq)%nqDescen
408    iq2 = tracers(iq)%iqDescen(ifils)
409    DO l = 1, llm
410      DO ij = iip2 + 1, ip1jm
411        q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2)
412      enddo
413      DO ij = iip1 + iip1, ip1jm, iip1
414        q(ij - iim, l, iq2) = q(ij, l, iq2)
415      enddo ! DO ij=ijb+iip1-1,ije,iip1
416    enddo !DO l=1,llm
417  enddo
418
419END SUBROUTINE vlx
420RECURSIVE SUBROUTINE vly(q, pente_max, masse, masse_adv_v, iq)
421  USE infotrac, ONLY: nqtot, tracers, & ! CRisi
422          min_qParent, min_qMass, min_ratio ! MVals et CRisi
423  USE lmdz_ssum_scopy, ONLY: ssum
424  USE lmdz_comgeom
425
426 !
427  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
428  !
429  !    ********************************************************************
430  ! Shema  d'advection " pseudo amont " .
431  !    ********************************************************************
432  ! q,masse_adv_v,w sont des arguments d'entree  pour le s-pg ....
433  ! dq         sont des arguments de sortie pour le s-pg ....
434  !
435  !
436  !   --------------------------------------------------------------------
437  USE comconst_mod, ONLY: pi
438  IMPLICIT NONE
439  !
440  INCLUDE "dimensions.h"
441  INCLUDE "paramet.h"
442  !
443  !
444  !   Arguments:
445  !   ----------
446  REAL :: masse(ip1jmp1, llm, nqtot), pente_max
447  REAL :: masse_adv_v(ip1jm, llm)
448  REAL :: q(ip1jmp1, llm, nqtot)
449  INTEGER :: iq ! CRisi
450  !
451  !  Local
452  !   ---------
453  !
454  INTEGER :: i, ij, l
455  !
456  REAL :: airej2, airejjm, airescb(iim), airesch(iim)
457  REAL :: dyq(ip1jmp1, llm), dyqv(ip1jm)
458  REAL :: adyqv(ip1jm), dyqmax(ip1jmp1)
459  REAL :: qbyv(ip1jm, llm)
460
461  REAL :: qpns, qpsn, dyn1, dys1, dyn2, dys2, newmasse, fn, fs
462  LOGICAL, SAVE :: first
463
464  REAL :: convpn, convps, convmpn, convmps
465  REAL :: massepn, masseps, qpn, qps
466  REAL :: sinlon(iip1), sinlondlon(iip1)
467  REAL :: coslon(iip1), coslondlon(iip1)
468  SAVE sinlon, coslon, sinlondlon, coslondlon
469  SAVE airej2, airejjm
470
471  REAL :: masseq(ip1jmp1, llm, nqtot), Ratio(ip1jmp1, llm, nqtot) ! CRisi
472  INTEGER :: ifils, iq2 ! CRisi
473
474  DATA first/.TRUE./
475
476  ! !WRITE(*,*) 'vly 578: entree, iq=',iq
477
478  IF(first) THEN
479    PRINT*, 'Shema  Amont nouveau  appele dans  Vanleer   '
480    first = .FALSE.
481    do i = 2, iip1
482      coslon(i) = cos(rlonv(i))
483      sinlon(i) = sin(rlonv(i))
484      coslondlon(i) = coslon(i) * (rlonu(i) - rlonu(i - 1)) / pi
485      sinlondlon(i) = sinlon(i) * (rlonu(i) - rlonu(i - 1)) / pi
486    ENDDO
487    coslon(1) = coslon(iip1)
488    coslondlon(1) = coslondlon(iip1)
489    sinlon(1) = sinlon(iip1)
490    sinlondlon(1) = sinlondlon(iip1)
491    airej2 = SSUM(iim, aire(iip2), 1)
492    airejjm = SSUM(iim, aire(ip1jm - iim), 1)
493  ENDIF
494
495  !
496  !PRINT*,'CALCUL EN LATITUDE'
497
498  DO l = 1, llm
499    !
500    !   --------------------------------
501    !  CALCUL EN LATITUDE
502    !   --------------------------------
503
504    !   On commence par calculer la valeur du traceur moyenne sur le premier cercle
505    !   de latitude autour du pole (qpns pour le pole nord et qpsn pour
506    !    le pole nord) qui sera utilisee pour evaluer les pentes au pole.
507
508    DO i = 1, iim
509      airescb(i) = aire(i + iip1) * q(i + iip1, l, iq)
510      airesch(i) = aire(i + ip1jm - iip1) * q(i + ip1jm - iip1, l, iq)
511    ENDDO
512    qpns = SSUM(iim, airescb, 1) / airej2
513    qpsn = SSUM(iim, airesch, 1) / airejjm
514
515    !   calcul des pentes aux points v
516
517    DO ij = 1, ip1jm
518      dyqv(ij) = q(ij, l, iq) - q(ij + iip1, l, iq)
519      adyqv(ij) = abs(dyqv(ij))
520    ENDDO
521
522    !   calcul des pentes aux points scalaires
523
524    DO ij = iip2, ip1jm
525      dyq(ij, l) = .5 * (dyqv(ij - iip1) + dyqv(ij))
526      dyqmax(ij) = min(adyqv(ij - iip1), adyqv(ij))
527      dyqmax(ij) = pente_max * dyqmax(ij)
528    ENDDO
529
530    !   calcul des pentes aux poles
531
532    DO ij = 1, iip1
533      dyq(ij, l) = qpns - q(ij + iip1, l, iq)
534      dyq(ip1jm + ij, l) = q(ip1jm + ij - iip1, l, iq) - qpsn
535    ENDDO
536
537    !   filtrage de la derivee
538    dyn1 = 0.
539    dys1 = 0.
540    dyn2 = 0.
541    dys2 = 0.
542    DO ij = 1, iim
543      dyn1 = dyn1 + sinlondlon(ij) * dyq(ij, l)
544      dys1 = dys1 + sinlondlon(ij) * dyq(ip1jm + ij, l)
545      dyn2 = dyn2 + coslondlon(ij) * dyq(ij, l)
546      dys2 = dys2 + coslondlon(ij) * dyq(ip1jm + ij, l)
547    ENDDO
548    DO ij = 1, iip1
549      dyq(ij, l) = dyn1 * sinlon(ij) + dyn2 * coslon(ij)
550      dyq(ip1jm + ij, l) = dys1 * sinlon(ij) + dys2 * coslon(ij)
551    ENDDO
552
553    !   calcul des pentes limites aux poles
554
555    goto 8888
556    fn = 1.
557    fs = 1.
558    DO ij = 1, iim
559      IF(pente_max * adyqv(ij)<abs(dyq(ij, l))) THEN
560        fn = min(pente_max * adyqv(ij) / abs(dyq(ij, l)), fn)
561      ENDIF
562      IF(pente_max * adyqv(ij + ip1jm - iip1)<abs(dyq(ij + ip1jm, l))) THEN
563        fs = min(pente_max * adyqv(ij + ip1jm - iip1) / abs(dyq(ij + ip1jm, l)), fs)
564      ENDIF
565    ENDDO
566    DO ij = 1, iip1
567      dyq(ij, l) = fn * dyq(ij, l)
568      dyq(ip1jm + ij, l) = fs * dyq(ip1jm + ij, l)
569    ENDDO
570    8888   continue
571    DO ij = 1, iip1
572      dyq(ij, l) = 0.
573      dyq(ip1jm + ij, l) = 0.
574    ENDDO
575
576    !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
577    !  En memoire de dIFferents tests sur la
578    !  limitation des pentes aux poles.
579    !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
580    ! PRINT*,dyq(1)
581    ! PRINT*,dyqv(iip1+1)
582    ! appn=abs(dyq(1)/dyqv(iip1+1))
583    ! PRINT*,dyq(ip1jm+1)
584    ! PRINT*,dyqv(ip1jm-iip1+1)
585    ! apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
586    ! DO ij=2,iim
587    !    appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
588    !    apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
589    ! ENDDO
590    ! appn=min(pente_max/appn,1.)
591    ! apps=min(pente_max/apps,1.)
592    !
593    !
594    !   cas ou on a un extremum au pole
595    !
596    ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
597    !    &   appn=0.
598    ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
599    !    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
600    !    &   apps=0.
601    !
602    !   limitation des pentes aux poles
603    ! DO ij=1,iip1
604    !    dyq(ij)=appn*dyq(ij)
605    !    dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
606    ! ENDDO
607    !
608    !   test
609    !  DO ij=1,iip1
610    !     dyq(iip1+ij)=0.
611    !     dyq(ip1jm+ij-iip1)=0.
612    !  ENDDO
613    !  DO ij=1,ip1jmp1
614    !     dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
615    !  ENDDO
616    !
617    ! changement 10 07 96
618    ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
619    !    &   THEN
620    !    DO ij=1,iip1
621    !       dyqmax(ij)=0.
622    !    ENDDO
623    ! ELSE
624    !    DO ij=1,iip1
625    !       dyqmax(ij)=pente_max*abs(dyqv(ij))
626    !    ENDDO
627    ! ENDIF
628    !
629    ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
630    !    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
631    !    &THEN
632    !    DO ij=ip1jm+1,ip1jmp1
633    !       dyqmax(ij)=0.
634    !    ENDDO
635    ! ELSE
636    !    DO ij=ip1jm+1,ip1jmp1
637    !       dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
638    !    ENDDO
639    ! ENDIF
640    !   fin changement 10 07 96
641    !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
642
643    !   calcul des pentes limitees
644
645    DO ij = iip2, ip1jm
646      IF(dyqv(ij) * dyqv(ij - iip1)>0.) THEN
647        dyq(ij, l) = sign(min(abs(dyq(ij, l)), dyqmax(ij)), dyq(ij, l))
648      ELSE
649        dyq(ij, l) = 0.
650      ENDIF
651    ENDDO
652
653  ENDDO
654
655  ! !WRITE(*,*) 'vly 756'
656  DO l = 1, llm
657    DO ij = 1, ip1jm
658      IF(masse_adv_v(ij, l)>0) THEN
659        qbyv(ij, l) = q(ij + iip1, l, iq) + dyq(ij + iip1, l) * &
660                0.5 * (1. - masse_adv_v(ij, l) &
661                / masse(ij + iip1, l, iq))
662      ELSE
663        qbyv(ij, l) = q(ij, l, iq) - dyq(ij, l) * &
664                0.5 * (1. + masse_adv_v(ij, l) &
665                / masse(ij, l, iq))
666      ENDIF
667      qbyv(ij, l) = masse_adv_v(ij, l) * qbyv(ij, l)
668    ENDDO
669  ENDDO
670
671  ! CRisi: appel récursif de l'advection sur les fils.
672  ! Il faut faire ça avant d'avoir mis à jour q et masse
673  ! WRITE(*,*) 'vly 689: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen
674
675  do ifils = 1, tracers(iq)%nqDescen
676    iq2 = tracers(iq)%iqDescen(ifils)
677    DO l = 1, llm
678      DO ij = 1, ip1jmp1
679        ! ! attention, chaque fils doit avoir son masseq, sinon, le 1er
680        ! ! fils ecrase le masseq de ses freres.
681        ! !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
682        !           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
683        ! !MVals: veiller a ce qu'on n'ait pas de denominateur nul
684        masseq(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass)
685        IF (q(ij, l, iq)>min_qParent) THEN
686          Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq)
687        else
688          Ratio(ij, l, iq2) = min_ratio
689        endif
690      enddo
691    enddo
692  enddo
693
694  do ifils = 1, tracers(iq)%nqDescen
695    iq2 = tracers(iq)%iqDescen(ifils)
696    CALL vly(Ratio, pente_max, masseq, qbyv, iq2)
697  enddo
698
699  DO l = 1, llm
700    DO ij = iip2, ip1jm
701      newmasse = masse(ij, l, iq) &
702              + masse_adv_v(ij, l) - masse_adv_v(ij - iip1, l)
703      q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) + qbyv(ij, l) &
704              - qbyv(ij - iip1, l)) / newmasse
705      masse(ij, l, iq) = newmasse
706    ENDDO
707    convpn = SSUM(iim, qbyv(1, l), 1)
708    convmpn = ssum(iim, masse_adv_v(1, l), 1)
709    massepn = ssum(iim, masse(1, l, iq), 1)
710    qpn = 0.
711    do ij = 1, iim
712      qpn = qpn + masse(ij, l, iq) * q(ij, l, iq)
713    enddo
714    qpn = (qpn + convpn) / (massepn + convmpn)
715    do ij = 1, iip1
716      q(ij, l, iq) = qpn
717    enddo
718    convps = -SSUM(iim, qbyv(ip1jm - iim, l), 1)
719    convmps = -ssum(iim, masse_adv_v(ip1jm - iim, l), 1)
720    masseps = ssum(iim, masse(ip1jm + 1, l, iq), 1)
721    qps = 0.
722    do ij = ip1jm + 1, ip1jmp1 - 1
723      qps = qps + masse(ij, l, iq) * q(ij, l, iq)
724    enddo
725    qps = (qps + convps) / (masseps + convmps)
726    do ij = ip1jm + 1, ip1jmp1
727      q(ij, l, iq) = qps
728    enddo
729  ENDDO
730
731  ! retablir les fils en rapport de melange par rapport a l'air:
732  do ifils = 1, tracers(iq)%nqDescen
733    iq2 = tracers(iq)%iqDescen(ifils)
734    DO l = 1, llm
735      DO ij = 1, ip1jmp1
736        q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2)
737      enddo
738    enddo
739  enddo
740
741  ! !WRITE(*,*) 'vly 853: sortie'
742
743END SUBROUTINE vly
744RECURSIVE SUBROUTINE vlz(q, pente_max, masse, w, iq)
745  USE infotrac, ONLY: nqtot, tracers, & ! CRisi
746          min_qParent, min_qMass, min_ratio ! MVals et CRisi
747  !
748  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
749  !
750  !    ********************************************************************
751  ! Shema  d'advection " pseudo amont " .
752  !    ********************************************************************
753  !    q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
754  ! dq         sont des arguments de sortie pour le s-pg ....
755  !   --------------------------------------------------------------------
756  IMPLICIT NONE
757  !
758  INCLUDE "dimensions.h"
759  INCLUDE "paramet.h"
760  !
761  !
762  !   Arguments:
763  !   ----------
764  REAL :: masse(ip1jmp1, llm, nqtot), pente_max
765  REAL :: q(ip1jmp1, llm, nqtot)
766  REAL :: w(ip1jmp1, llm + 1)
767  INTEGER :: iq
768  !
769  !  Local
770  !   ---------
771  !
772  INTEGER :: ij, l
773  !
774  REAL :: wq(ip1jmp1, llm + 1), newmasse
775
776  REAL :: dzq(ip1jmp1, llm), dzqw(ip1jmp1, llm), adzqw(ip1jmp1, llm), dzqmax
777  REAL :: sigw
778
779  REAL :: masseq(ip1jmp1, llm, nqtot), Ratio(ip1jmp1, llm, nqtot) ! CRisi
780  INTEGER :: ifils, iq2 ! CRisi
781
782#ifdef BIDON
783  REAL :: temps0,temps1,second
784  SAVE temps0,temps1
785
786  DATA temps0,temps1/0.,0./
787#endif
788
789  !    On oriente tout dans le sens de la pression c'est a dire dans le
790  !    sens de W
791  DO l = 2, llm
792    DO ij = 1, ip1jmp1
793      dzqw(ij, l) = q(ij, l - 1, iq) - q(ij, l, iq)
794      adzqw(ij, l) = abs(dzqw(ij, l))
795    ENDDO
796  ENDDO
797
798  DO l = 2, llm - 1
799    DO ij = 1, ip1jmp1
800      IF(dzqw(ij, l) * dzqw(ij, l + 1)>0.) THEN
801        dzq(ij, l) = 0.5 * (dzqw(ij, l) + dzqw(ij, l + 1))
802      ELSE
803        dzq(ij, l) = 0.
804      ENDIF
805      dzqmax = pente_max * min(adzqw(ij, l), adzqw(ij, l + 1))
806      dzq(ij, l) = sign(min(abs(dzq(ij, l)), dzqmax), dzq(ij, l))
807    ENDDO
808  ENDDO
809
810  ! !WRITE(*,*) 'vlz 954'
811  DO ij = 1, ip1jmp1
812    dzq(ij, 1) = 0.
813    dzq(ij, llm) = 0.
814  ENDDO
815
816  ! ---------------------------------------------------------------
817  !   .... calcul des termes d'advection verticale  .......
818  ! ---------------------------------------------------------------
819
820  ! calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
821
822  DO l = 1, llm - 1
823    do  ij = 1, ip1jmp1
824      IF(w(ij, l + 1)>0.) THEN
825        sigw = w(ij, l + 1) / masse(ij, l + 1, iq)
826        wq(ij, l + 1) = w(ij, l + 1) * (q(ij, l + 1, iq) &
827                + 0.5 * (1. - sigw) * dzq(ij, l + 1))
828      ELSE
829        sigw = w(ij, l + 1) / masse(ij, l, iq)
830        wq(ij, l + 1) = w(ij, l + 1) * (q(ij, l, iq) - 0.5 * (1. + sigw) * dzq(ij, l))
831      ENDIF
832    ENDDO
833  ENDDO
834
835  DO ij = 1, ip1jmp1
836    wq(ij, llm + 1) = 0.
837    wq(ij, 1) = 0.
838  ENDDO
839
840  ! CRisi: appel récursif de l'advection sur les fils.
841  ! Il faut faire ça avant d'avoir mis à jour q et masse
842  ! !WRITE(*,*) 'vlsplt 942: iq,nqChildren(iq)=',iq,nqChildren(iq)
843  do ifils = 1, tracers(iq)%nqDescen
844    iq2 = tracers(iq)%iqDescen(ifils)
845    DO l = 1, llm
846      DO ij = 1, ip1jmp1
847        ! !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
848        !           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
849        ! !MVals: veiller a ce qu'on n'ait pas de denominateur nul
850        masseq(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass)
851        IF (q(ij, l, iq)>min_qParent) THEN
852          Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq)
853        else
854          Ratio(ij, l, iq2) = min_ratio
855        endif
856      enddo
857    enddo
858  enddo
859
860  do ifils = 1, tracers(iq)%nqChildren
861    iq2 = tracers(iq)%iqDescen(ifils)
862    CALL vlz(Ratio, pente_max, masseq, wq, iq2)
863  enddo
864  ! end CRisi
865
866  DO l = 1, llm
867    DO ij = 1, ip1jmp1
868      newmasse = masse(ij, l, iq) + w(ij, l + 1) - w(ij, l)
869      q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) + wq(ij, l + 1) - wq(ij, l)) &
870              / newmasse
871      masse(ij, l, iq) = newmasse
872    ENDDO
873  ENDDO
874
875  ! retablir les fils en rapport de melange par rapport a l'air:
876  do ifils = 1, tracers(iq)%nqDescen
877    iq2 = tracers(iq)%iqDescen(ifils)
878    DO l = 1, llm
879      DO ij = 1, ip1jmp1
880        q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2)
881      enddo
882    enddo
883  enddo
884
885END SUBROUTINE vlz
886
887SUBROUTINE minmaxq(zq, qmin, qmax, comment)
888
889  INCLUDE "dimensions.h"
890  INCLUDE "paramet.h"
891
892  CHARACTER(LEN = 20) :: comment
893  REAL :: qmin, qmax
894  REAL :: zq(ip1jmp1, llm)
895  REAL :: zzq(iip1, jjp1, llm)
896
897END SUBROUTINE  minmaxq
898
899
900
Note: See TracBrowser for help on using the repository browser.