Changeset 3780 for LMDZ6/trunk/libf/phylmd/yamada4.F90

Timestamp:

Oct 22, 2020, 2:50:18 PM (4 years ago)

Author:

evignon

Message:

Premiere comission Etienne: changements pour le 1D (forcage en Ts au dessus des continents) et inclusion drag arbres dans yamada4_num=6

File:

• LMDZ6/trunk/libf/phylmd/yamada4.F90 (modified) (13 diffs)

LMDZ6/trunk/libf/phylmd/yamada4.F90

@@ -6 +6 @@
   USE dimphy
   USE ioipsl_getin_p_mod, ONLY : getin_p
-
+  USE phys_local_var_mod, only: tke_dissip
+  
   IMPLICIT NONE
   include "iniprint.h"
@@ -56 +57 @@
   !             iflag_pbl=11 -> the model starts with NP from start files created by ce0l
   !                          -> the model can run with longer time-steps.
-  !             2016/11/30 (EV etienne.vignon@univ-grenoble-alpes.fr)
+  !             2016/11/30 (EV etienne.vignon@lmd.ipsl.fr)
   !               On met tke (=q2/2) en entr??e plut??t que q2
   !               On corrige l'update de la tke
-  !
+  !             2020/10/01 (EV)
+  !               On ajoute la dissipation de la TKE en diagnostique de sortie
+  !                  
   ! Inpout/Output :
   !==============
@@ -121 +124 @@
   REAL,SAVE :: viscom,viscoh
   !$OMP THREADPRIVATE( hboville,viscom,viscoh)
-  INTEGER ig, k
+  INTEGER ig, jg, k
   REAL ri, zrif, zalpha, zsm, zsn
   REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)
@@ -416 +419 @@
   ELSE IF (iflag_pbl>=10) THEN
 
+    shear(:,:)=0.
+    buoy(:,:)=0.
+    dissip(:,:)=0.
+    km(:,:)=0.
+        
     IF (yamada4_num>=1) THEN
  
@@ -424 +432 @@
       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))
+      dissip(ig,k)=min(max(((sqrt(q2(ig,k)))**3)/(b1*l(ig,k)),1.E-12),1.E4)
      ENDDO
     ENDDO
-
+    
     IF (yamada4_num==1) THEN ! Schema du MAR tel quel
        DO k = 2, klev - 1
@@ -478 +486 @@
         ENDDO
        ENDDO
+       
+      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+      !! Attention, yamada4_num=5 est inexacte car néglige les termes de flottabilité
+      !! en conditions instables
+      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     ELSE IF (yamada4_num==5) THEN ! version modifiee avec integration exacte pour la dissipation
        DO k = 2, klev - 1
@@ -507 +520 @@
        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 dissipeff/TKE est constant.
+         ! Puis on prend la solution exacte
+         ! 
+         ! With drag and dissipation from high vegetation (EV & FC, 05/10/2020)
+         winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2)
          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))
+         tkeprov=tkeprov+max(buoy(ig,k)+shear(ig,k)+drgpro(ig,k)*(winds(ig,k))**3,0.)*dt
+         disseff=dissip(ig,k)-min(0.,buoy(ig,k)+shear(ig,k)+drgpro(ig,k)*(winds(ig,k))**3) + drgpro(ig,k)*tkeprov
          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
+         q2(ig,k)=tkeprov*ydeux        
+         
         ENDDO
        ENDDO
@@ -533 +556 @@
 !     q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k))
       q2(1:ngrid, k) = min(max(q2(1:ngrid,k),1.E-10), 1.E4)
-       q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(yun*l(1:ngrid,k)*b1))
+      q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(yun*l(1:ngrid,k)*b1))
 !     q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1))
       q2(1:ngrid, k) = q2(1:ngrid, k)*q2(1:ngrid, k)
@@ -725 +748 @@
 
 !============================================================================
+! Diagnostique de la dissipation
+!============================================================================
+
+! Diagnostics
+ tke_dissip(1:ngrid,:,nsrf)=0.
+ DO k=2,klev
+    DO ig=1,ngrid
+       jg=ni(ig)
+       tke_dissip(jg,k,nsrf)=dissip(ig,k)
+    ENDDO
+ ENDDO
+ 
+!=============================================================================
 
   RETURN
@@ -1017 +1053 @@
 !=====================================================================
 
-
+  l1(1:ngrid,:)=0.
   IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
 
@@ -1069 +1105 @@
    l2(1:ngrid,:)=0.0
    l_mixmin(1:ngrid,:,nsrf)=0.
-   l_mix(1:ngrid,:,nsrf)=0.
+   l_mix(1:ngrid,:,nsrf)=1.E-5
 
    IF (nsrf .EQ. 1) THEN
@@ -1135 +1171 @@
 
 
- DO k=2,klev
+ DO k=1,klev+1
     DO ig=1,ngrid
        lmix(ig,k)=MAX(MAX(l1(ig,k), l2(ig,k)),lmixmin)
+       lmix(ig,k)=MAX(lmix(ig,k),1.E-5)
    ENDDO
  ENDDO
@@ -1143 +1180 @@
 ! Diagnostics
 
- DO k=2,klev
+ DO k=2,klev+1
     DO ig=1,ngrid
        jg=ni(ig)