Changeset 2765 for LMDZ5

Timestamp:

Jan 19, 2017, 3:15:40 PM (8 years ago)

Author:

fhourdin

Message:

Inclusion d'un traitement spécifique dans l'advection verticale
pour les cas de nombre de courant w/masse > 1
ou w est le transfert de masse entre deux couches.

ce traitement garanti la convergence numérique quand le critère
n'est jamais violé.
Dans les autres cas, il fait proprement le transport vertical
qui est autrement hors de contrôle.
Pourrait contribuer a resoudre les plantages himalayens.

File:

• LMDZ5/trunk/libf/dyn3dmem/vlsplt_loc.F (modified) (4 diffs)

LMDZ5/trunk/libf/dyn3dmem/vlsplt_loc.F

@@ -885 +885 @@
 c
       INTEGER i,ij,l,j,ii
+
+      REAL,DIMENSION(ijb_u:ije_u,llm+1) :: wresi,morig,qorig,dzqorig
+      INTEGER,DIMENSION(ijb_u:ije_u,llm+1) :: lorig
+      INTEGER,SAVE :: countcfl
+!$OMP THREADPRIVATE(countcfl)
 c
       REAL newmasse
@@ -911 +916 @@
       ! vlz_loc si on veut qu'elles soient vues par tous les threads.  
       INTEGER ifils,iq2 ! CRisi
+
 
       IF (first) THEN
@@ -968 +974 @@
       ENDIF
 #endif
+
+!--------------------------------------------------------
+! On repere les points qui violent le CFL (|w| > masse)
+!--------------------------------------------------------
+
+      countcfl=0
+!     print*,'vlz nouveau'
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+      DO l = 2,llm
+         DO ij = ijb,ije
+          IF(  (w(ij,l,iq)>0.AND.w(ij,l,iq)>masse(ij,l,iq))
+     s    .OR. (w(ij,l,iq)<=0.AND.ABS(w(ij,l,iq))>masse(ij,l-1,iq)) )
+     s    countcfl=countcfl+1
+         ENDDO
+      ENDDO
+c$OMP END DO NOWAIT    
+
+c ---------------------------------------------------------------
+c  Identification des mailles ou on viole le CFL : w > masse
+c ---------------------------------------------------------------
+
+      IF (countcfl==0) THEN
+
 c ---------------------------------------------------------------
 c   .... calcul des termes d'advection verticale  .......
+c     Dans le cas où le  |w| < masse partout.
+c     Version d'origine
+c     Pourrait etre enleve si on voit que le code plus general
+c     est aussi rapide
 c ---------------------------------------------------------------
 
@@ -991 +1024 @@
 c$OMP END DO NOWAIT    
        !write(*,*) 'vlz 1001'    
+
+      ELSE ! countcfl>=1
+
+      PRINT*,'vlz passage dans le non local'
+c ---------------------------------------------------------------
+c  Debut du traitement du cas ou on viole le CFL : w > masse
+c ---------------------------------------------------------------
+
+c Initialisation
+
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+       DO l = 2,llm
+         DO ij = ijb,ije
+            wresi(ij,l)=w(ij,l,iq)
+            wq(ij,l,iq)=0.
+            IF(w(ij,l,iq).gt.0.) THEN
+               lorig(ij,l)=l
+               morig(ij,l)=masse(ij,l,iq)
+               qorig(ij,l)=q(ij,l,iq)
+               dzqorig(ij,l)=dzq(ij,l)
+            ELSE
+               lorig(ij,l)=l-1
+               morig(ij,l)=masse(ij,l-1,iq)
+               qorig(ij,l)=q(ij,l-1,iq)
+               dzqorig(ij,l)=dzq(ij,l-1)
+            ENDIF
+         ENDDO
+       ENDDO
+c$OMP END DO NO WAIT
+
+c Reindicage vertical en accumulant les flux sur
+c  les mailles qui viollent le CFL
+c  on itère jusqu'à ce que tous les poins satisfassent
+c  le critère
+      DO WHILE (countcfl>=1) 
+      print*,'On viole le CFL Vertical sur ',countcfl,' pts'
+      countcfl=0
+
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+      DO l = 2,llm
+         DO ij = ijb,ije
+          IF (ABS(wresi(ij,l))>morig(ij,l)) THEN
+             countcfl=countcfl+1
+! rm : les 8 lignes ci dessous pourraient sans doute s'ecrire
+! avec la fonction sign
+             IF(w(ij,l,iq)>0.) THEN
+                wresi(ij,l)=wresi(ij,l)-morig(ij,l)
+                wq(ij,l,iq)=wq(ij,l,iq)+morig(ij,l)*qorig(ij,l)
+                lorig(ij,l)=lorig(ij,l)+1
+             ELSE
+                wresi(ij,l)=wresi(ij,l)+morig(ij,l)
+                wq(ij,l,iq)=wq(ij,l,iq)-morig(ij,l)*qorig(ij,l)
+                lorig(ij,l)=lorig(ij,l)-1
+             ENDIF
+             morig(ij,l)=masse(ij,lorig(ij,l),iq)
+             qorig(ij,l)=q(ij,lorig(ij,l),iq)
+             dzqorig(ij,l)=dzq(ij,lorig(ij,l))
+          ENDIF
+         ENDDO
+      ENDDO
+c$OMP END DO NO WAIT
+
+      ENDDO ! WHILE (countcfl>=1)
+
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+       DO l = 2,llm
+         do  ij = ijb,ije
+          sigw=wresi(ij,l)/morig(ij,l)
+          IF(w(ij,l,iq).gt.0.) THEN
+             wq(ij,l,iq)=wq(ij,l,iq)+wresi(ij,l)*(qorig(ij,l)
+     :           +0.5*(1.-sigw)*dzqorig(ij,l))
+          ELSE
+             wq(ij,l,iq)=wq(ij,l,iq)+wresi(ij,l)*(qorig(ij,l)
+     :           -0.5*(1.+sigw)*dzqorig(ij,l))
+          ENDIF
+         ENDDO
+       ENDDO
+c$OMP END DO NOWAIT    
+
+
+       ENDIF ! councfl=0
+
+
 
 c$OMP MASTER