source: LMDZ6/branches/Amaury_dev/libf/dyn3dmem/vlspltqs_loc.f90 @ 5112

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