Changeset 2817 for LMDZ5/trunk/libf

Timestamp:

Mar 7, 2017, 10:05:20 AM (8 years ago)

Author:

fhourdin

Message:

En option, nouveaux schemas numeriques pour l'integration temporelle.
yamada4_num=0 : std

=1 : comme MAR

File:

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

LMDZ5/trunk/libf/phylmd/yamada4.F90

@@ -57 +57 @@
   !                          -> the model can run with longer time-steps.
   !             2016/11/30 (EV etienne.vignon@univ-grenoble-alpes.fr)
-  !               On met tke (=q2/2) en entrée plutôt que q2
+  !               On met tke (=q2/2) en entr??e plut??t que q2
   !               On corrige l'update de la tke
   !
@@ -130 +130 @@
   REAL zz(klon, klev+1)
   INTEGER iter
+  REAL dissip(klon,klev), tkeprov, shear(klon,klev), buoy(klon,klev)
   REAL,SAVE :: ric0,ric,rifc, b1, kap
   !$OMP THREADPRIVATE(ric0,ric,rifc,b1,kap)
@@ -138 +139 @@
   !$OMP THREADPRIVATE(lmixmin)
   LOGICAL, SAVE :: new_yamada4
-  !$OMP THREADPRIVATE(new_yamada4)
-  REAL, SAVE :: yun,ydeux
+  INTEGER, SAVE :: yamada4_num
+  !$OMP THREADPRIVATE(new_yamada4,yamada4_num)
+  REAL, SAVE :: yun,ydeux,b1leff
   !$OMP THREADPRIVATE(yun,ydeux)
   REAL frif, falpha, fsm
@@ -161 +163 @@
     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
@@ -169 +173 @@
        CALL getin_p('yamada4_ric',ric)
        ric0=0.19489      ! ric=0.195 originalement, mais produisait sm<0
-       ric=min(ric,ric0) ! Au delà de ric0, sm devient négatif
+       ric=min(ric,ric0) ! Au dela de ric0, sm devient n??gatif
        rifc=frif(ric)
        seuilsm=fsm(frif(ric))
@@ -290 +294 @@
   !=======================================================================
 
-  ! On commence par calculé q2 à partir de la tke
+  ! On commence par calculer q2 a partir de la tke
 
   IF (new_yamada4) THEN
@@ -400 +404 @@
   ELSE IF (iflag_pbl>=10) THEN
 
+    IF (yamada4_num>=1) THEN
+ 
+    DO k = 2, klev - 1
+      DO ig=1,ngrid
+      q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
+      km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
+      shear(ig,k)=km(ig, k)*m2(ig, k)
+      buoy(ig,k)=km(ig, k)*m2(ig, k)*(-1.*rif(ig,k))
+      dissip(ig,k)=((sqrt(q2(ig,k)))**3)/(b1*l(ig,k))
+     ENDDO
+    ENDDO
+
+    IF (yamada4_num==1) THEN ! Schema du MAR tel quel
+       DO k = 2, klev - 1
+         DO ig=1,ngrid
+         tkeprov=q2(ig,k)/ydeux
+         q2(ig,k)= tkeprov*                           &
+           &  (tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k))))/ &
+           &  (tkeprov+dt*((-1.)*min(0.,buoy(ig,k))+dissip(ig,k)))
+        ENDDO
+       ENDDO
+    ELSE ! 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)))
+         ! 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
+    ENDIF
+
+    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
+    ENDDO
+
+   ELSE 
+
     DO k = 2, klev - 1
       km(1:ngrid, k) = l(1:ngrid, k)*sqrt(q2(1:ngrid,k))*sm(1:ngrid, k)
@@ -410 +458 @@
     END DO
 
+  ENDIF
 
   ELSE
@@ -575 +624 @@
 
 !============================================================================
-! Mise à jour de la tke
+! Mise a jour de la tke
 !============================================================================