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

Last change on this file since 5112 was 5105, checked in by abarral, 4 months ago

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