Changeset 290 for trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90

Timestamp:

Sep 9, 2011, 4:24:40 PM (13 years ago)

Author:

acolaitis

Message:

--- AC 07/09/2011 ---

This is a major update for thermals and richardson layer parametrization. This update stabilizes thermals (further
studies might show that we can reduce the value of nsplit in gcm. To be tested.), and prevent downdrafts from descending into
the surface layer, which was acting as a cold feedback on the thermals. The Richardson surface layer now features
different gustiness coefficients for Richardson, heat and momentum so that u* and t* are correctly represented.

Upcoming updates will change surflayer_interpol.F90 to implement those changes in the interpolation scheme as well.

*
IMPORTANT : several parts of the vdifc code might want to use these new definitions for gustiness, u* and t*. exemple : dust devil routines
that recompute u* ? lifting routines ? TODO !

M 289 libf/phymars/thermcell_main_mars.F90
^{-----------------> Added iterations to the velocity / buoyancy / entrainment / detrainment}

computation to ensure correct convergence. Iteration number is for now set to
4, which is probably too much. This will be changed once several cases are tested.
The minimum is probably 2 iterations.

M 289 libf/phymars/vdifc.F
^{-----------------> Subgrid gustiness parametrization now uses different gustiness beta coefficients}

for heat and momentum. Comparisons with LES at Exomars landing site, matching u*
and teta* suggests values of beta=0.7 for momentum and beta=1.2 for heat, where
the formula for large scale horizontal winds in the first layer is :

U0^{2 = pu}2 + pv^{2 + (beta*wmax_th)}2

M 289 libf/phymars/vdif_cd.F
^{-----------------> LES data suggests that Richardson number distribution during daytime has a very large}

standart deviation due to highly negative Richardsons on several gridpoints. As a consequence,
the mean Richardson in daytime in the LES is much more negative than in LES. In the gcm,
parametrized subgrid turbulence prevents such cases, which can be dangerous in nearly unstable conditions.
Hence, we use beta=0.5 instead of one, so that we keep the anti-crazy-hfx function of beta and we increase the
norm of negative Richardsons in general for daytime conditions. This is set in conjonction with
beta settings for heat and momentum in vdifc.

M 289 libf/phymars/meso_inc/meso_inc_les.F
^{-----------------> HFX and USTM computations now uses the different betas for heat and momentum.}

---

File:

• trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90 (modified) (7 diffs)

trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90

@@ -114 +114 @@
       LOGICAL active(ngridmx),activetmp(ngridmx)
       REAL a1,b1,ae,be,ad,bd
+      INTEGER tic
 
 ! ==========================================
@@ -334 +335 @@
                alim_star(ig,1)=MAX((ztv(ig,1)-ztv(ig,2)),0.)  &
      &                       *sqrt(zlev(ig,2))
+!     &                       /sqrt(zlev(ig,2))
 !      &                       *zlev(ig,2)
                lalim(ig)=2
@@ -381 +383 @@
 ! dans un panache conservatif
           f_star(ig,1)=0.
+          
           f_star(ig,2)=alim_star(ig,1)
+
           zw2(ig,2)=2.*g*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
       &                     *(zlev(ig,2)-zlev(ig,1))  &
-      &                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
+      &                    *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
           w_est(ig,2)=zw2(ig,2)
 
@@ -420 +424 @@
 
 !                if(l .lt. lalim(ig)) then
-!                ztva_est(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
+!               ztva_est(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
 !     &            alim_star(ig,l)*ztv(ig,l))  &
 !     &            /(f_star(ig,l)+alim_star(ig,l))
@@ -522 +526 @@
 !---------------------------------------------------------------------------
 
+      DO tic=0,3
+
       activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
       do ig=1,ngridmx
@@ -547 +553 @@
       enddo
 
+! ================ RECOMPUTE ENTR, DETR, and F FROM NEW W2 ===================
+
+
+      do ig=1,ngridmx
+        if (activetmp(ig)) then
+
+          zw2m=zw2(ig,l+1)
+          if(zw2m .gt. 0) then
+          if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then
+          entr_star(ig,l)=f_star(ig,l)*zdz*  &
+        &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
+          else
+          entr_star(ig,l)=0.
+          endif
+
+          if(zbuoy(ig,l) .gt. 0.) then
+             if(l .lt. lalim(ig)) then
+                detr_star(ig,l)=0.
+             else
+                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
+            &  ad
+             endif
+          else
+          detr_star(ig,l)=f_star(ig,l)*zdz*                        &
+                &     bd*zbuoy(ig,l)/zw2m
+          endif
+          else
+          entr_star(ig,l)=0.
+          detr_star(ig,l)=0.
+          endif
+
+! En dessous de lalim, on prend le max de alim_star et entr_star pour
+! alim_star et 0 sinon
+
+        if (l.lt.lalim(ig)) then
+          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
+          entr_star(ig,l)=0.
+        endif
+
+! Calcul du flux montant normalise
+
+      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
+     &              -detr_star(ig,l)
+
+      endif
+      enddo
+
+      ENDDO
 !---------------------------------------------------------------------------
 !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
@@ -1121 +1175 @@
       do ig=1,ngridmx
          if(lmax(ig) .gt. 1) then
-         do k=1,lmax(ig)
+! No downdraft in the very-near surface layer, we begin at k=3
+         do k=3,lmax(ig)
             zdthladj_down(ig,k)=(1./masse(ig,k))*(fm_down(ig,k+1)*ztv(ig,k+1)- &
      & fm_down(ig,k)*ztv(ig,k)+fm_down(ig,k)*ztvd(ig,k)-fm_down(ig,k+1)*ztvd(ig,k+1))