Changeset 3798 for LMDZ6/branches/Ocean_skin/libf/phylmd/yamada4.F90

Timestamp:

Jan 11, 2021, 11:24:08 PM (3 years ago)

Author:

lguez

Message:

Sync latest trunk changes to Ocean_skin

Location:

LMDZ6/branches/Ocean_skin

Files:

• . (modified) (1 prop)
• libf/phylmd/yamada4.F90 (modified) (11 diffs)

LMDZ6/branches/Ocean_skin/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)
@@ -186 +189 @@
        viscom=1.46E-5
        viscoh=2.06E-5
+       !lmixmin=1.0E-3
        lmixmin=0.
        yamada4_num=5
@@ -416 +420 @@
   ELSE IF (iflag_pbl>=10) THEN
 
+    shear(:,:)=0.
+    buoy(:,:)=0.
+    dissip(:,:)=0.
+    km(:,:)=0.
+        
     IF (yamada4_num>=1) THEN
  
@@ -424 +433 @@
       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)
+      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
@@ -478 +488 @@
         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 +522 @@
        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
+         
         ENDDO
        ENDDO
@@ -725 +750 @@
 
 !============================================================================
+! 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 +1055 @@
 !=====================================================================
 
-
+  l1(1:ngrid,:)=0.
   IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
 
@@ -1135 +1173 @@
 
 
- 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)