Changeset 2828

Timestamp:

Mar 20, 2017, 6:05:17 PM (8 years ago)

Author:

fhourdin

Message:

Modifications de yamada4 issues de la these d'Etienne Vignon (sous

l'otpion new_yamada4=y) :
1/ schema numerique inspire de celui de MAR
2/ viscosite moleculaire ajoute pour les cas ou la diffusion

turbulente s'effondre.

3/ Desactivation de la modification Holtzlag Boville pour les

couches limites tres stables.

4/ longueur de melange nulle par defaut.

File:

• LMDZ5/trunk/libf/phylmd/yamada4.F90 (modified) (10 diffs)

LMDZ5/trunk/libf/phylmd/yamada4.F90

@@ -120 +120 @@
   DATA first, ipas/.FALSE., 0/
   !$OMP THREADPRIVATE( first,ipas)
-  LOGICAL hboville
+  LOGICAL,SAVE :: hboville=.TRUE.
+  REAL,SAVE :: viscom,viscoh
+  !$OMP THREADPRIVATE( hboville,viscom,viscoh)
   INTEGER ig, k
   REAL ri, zrif, zalpha, zsm, zsn
@@ -130 +132 @@
   REAL zz(klon, klev+1)
   INTEGER iter
-  REAL dissip(klon,klev), tkeprov, shear(klon,klev), buoy(klon,klev)
+  REAL dissip(klon,klev), tkeprov,tkeexp, shear(klon,klev), buoy(klon,klev)
   REAL,SAVE :: ric0,ric,rifc, b1, kap
   !$OMP THREADPRIVATE(ric0,ric,rifc,b1,kap)
@@ -141 +143 @@
   INTEGER, SAVE :: yamada4_num
   !$OMP THREADPRIVATE(new_yamada4,yamada4_num)
-  REAL, SAVE :: yun,ydeux,b1leff
+  REAL, SAVE :: yun,ydeux,disseff
   !$OMP THREADPRIVATE(yun,ydeux)
   REAL frif, falpha, fsm
@@ -163 +165 @@
     new_yamada4=.false.
     CALL getin_p('new_yamada4',new_yamada4)
-    yamada4_num=0
-    CALL getin_p('yamada4_num',yamada4_num)
-! Pour garantir la continuite dans la mise au point de CMIP6.
-! Eliminer l'option new_yamada4 des le printemps 2016
+
     IF (new_yamada4) THEN
+! Corrections et reglages issus du travail de these d'Etienne Vignon.
        ric=0.143 ! qui donne des valeurs proches des seuils proposes
                  ! dans YAMADA 1983 : sm=0.0845 (0.085 dans Y83)
@@ -179 +179 @@
        yun=1.
        ydeux=2.
-!      yun=2.
-!      ydeux=1.
+       hboville=.FALSE.
+       viscom=1.46E-5
+       viscoh=2.06E-5
+       lmixmin=0.
+       yamada4_num=5
     ELSE
-  !!  DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.4/
        ric=0.195
        rifc=0.191
@@ -189 +191 @@
        yun=2.
        ydeux=1.
+       hboville=.TRUE.
+       viscom=0.
+       viscoh=0.
+       lmixmin=1.
+       yamada4_num=0
     ENDIF
+
     PRINT*,'YAMADA4 RIc, RIfc, Sm_min, Alpha_min',ric,rifc,seuilsm,seuilalpha
     firstcall = .FALSE.
-    lmixmin=1.
     CALL getin_p('lmixmin',lmixmin)
+    CALL getin_p('yamada4_hboville',hboville)
+    CALL getin_p('yamada4_num',yamada4_num)
   END IF
 
@@ -203 +212 @@
 
 ! On utilise ou non la routine de Holstalg Boville pour les cas tres stables
-   hboville=.TRUE.
 
 
@@ -420 +428 @@
          DO ig=1,ngrid
          tkeprov=q2(ig,k)/ydeux
-         q2(ig,k)= tkeprov*                           &
+         tkeprov= tkeprov*                           &
            &  (tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k))))/ &
            &  (tkeprov+dt*((-1.)*min(0.,buoy(ig,k))+dissip(ig,k)))
+         q2(ig,k)=tkeprov*ydeux
         ENDDO
        ENDDO
-    ELSE ! version modifiee avec integration exacte pour la dissipation
+    ELSE IF (yamada4_num==2) THEN ! version modifiee avec integration exacte pour la dissipation
        DO k = 2, klev - 1
          DO ig=1,ngrid
          tkeprov=q2(ig,k)/ydeux
-         q2(ig,k)= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
+         disseff=dissip(ig,k)-min(0.,buoy(ig,k))
+         tkeprov = tkeprov/(1.+dt*disseff/(2.*tkeprov))**2
+         tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
+         q2(ig,k)=tkeprov*ydeux
          ! En cas stable, on traite la flotabilite comme la
          ! dissipation, en supposant que buoy/q2^3 est constant.
          ! Puis on prend la solution exacte
-         b1leff=1./(1./(b1*l(ig,k))-min(0.,buoy(ig,k)/sqrt(q2(ig,k)))**3)
-         q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(yun*b1leff))
+        ENDDO
+       ENDDO
+    ELSE IF (yamada4_num==3) THEN ! version modifiee avec integration exacte pour la dissipation
+       DO k = 2, klev - 1
+         DO ig=1,ngrid
+         tkeprov=q2(ig,k)/ydeux
+         disseff=dissip(ig,k)-min(0.,buoy(ig,k))
+         tkeprov=tkeprov*exp(-dt*disseff/tkeprov)
+         tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
+         q2(ig,k)=tkeprov*ydeux
+         ! En cas stable, on traite la flotabilite comme la
+         ! dissipation, en supposant que buoy/q2^3 est constant.
+         ! Puis on prend la solution exacte
+        ENDDO
+       ENDDO
+    ELSE IF (yamada4_num==4) THEN ! version modifiee avec integration exacte pour la dissipation
+       DO k = 2, klev - 1
+         DO ig=1,ngrid
+         tkeprov=q2(ig,k)/ydeux
+         tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
+         tkeprov= tkeprov*                           &
+           &  tkeprov/ &
+           &  (tkeprov+dt*((-1.)*min(0.,buoy(ig,k))+dissip(ig,k)))
+         q2(ig,k)=tkeprov*ydeux
+         ! En cas stable, on traite la flotabilite comme la
+         ! dissipation, en supposant que buoy/q2^3 est constant.
+         ! Puis on prend la solution exacte
+        ENDDO
+       ENDDO
+    ELSE IF (yamada4_num==5) THEN ! version modifiee avec integration exacte pour la dissipation
+       DO k = 2, klev - 1
+         DO ig=1,ngrid
+         tkeprov=q2(ig,k)/ydeux
+         disseff=dissip(ig,k)-min(0.,buoy(ig,k))
+         tkeexp=exp(-dt*disseff/tkeprov)
+         tkeprov= shear(ig,k)*tkeprov/disseff*(1.-tkeexp)+tkeprov*tkeexp
+         q2(ig,k)=tkeprov*ydeux
+         ! En cas stable, on traite la flotabilite comme la
+         ! dissipation, en supposant que buoy/q2^3 est constant.
+         ! Puis on prend la solution exacte
+        ENDDO
+       ENDDO
+    ELSE IF (yamada4_num==6) THEN ! version modifiee avec integration exacte pour la dissipation
+       DO k = 2, klev - 1
+         DO ig=1,ngrid
+         tkeprov=q2(ig,k)/ydeux
+         tkeprov=tkeprov+max(buoy(ig,k)+shear(ig,k),0.)*dt
+         disseff=dissip(ig,k)-min(0.,buoy(ig,k)+shear(ig,k))
+         tkeexp=exp(-dt*disseff/tkeprov)
+         tkeprov= tkeprov*tkeexp
+         q2(ig,k)=tkeprov*ydeux
+         ! En cas stable, on traite la flotabilite comme la
+         ! dissipation, en supposant que buoy/q2^3 est constant.
+         ! Puis on prend la solution exacte
         ENDDO
        ENDDO
@@ -441 +505 @@
     DO k = 2, klev - 1
       DO ig=1,ngrid
-      q2(ig,k)=q2(ig,k)*ydeux
       q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
       ENDDO
@@ -577 +640 @@
 
  END IF ! hboville
+
+! Ajout d'une viscosite moleculaire
+   km(:,:)=km(:,:)+viscom
+   kn(:,:)=kn(:,:)+viscoh
+   kq(:,:)=kq(:,:)+viscoh
 
   IF (prt_level>1) THEN