Changeset 2065 for trunk/LMDZ.GENERIC/libf

Timestamp:

Jan 11, 2019, 6:04:09 PM (6 years ago)

Author:

aboissinot

Message:

In thermcell_plume, replace a useless test on zalpha by a test on zw2m to avoid a possible
division by zero and remove useless variable zbuoybis.

In physiq_mod, save variable f0.

Location:

trunk/LMDZ.GENERIC/libf/phystd

Files:

• physiq_mod.F90 (modified) (1 diff)
• thermcell_plume.F90 (modified) (13 diffs)

trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90

@@ -262 +262 @@
       INTEGER zlcl(ngrid)
       
-      real f0(ngrid)                      ! Mass flux norm
+      real,save :: f0(ngrid)              ! Mass flux norm
       real fm0(ngrid, nlayer+1)           ! Mass flux
       real entr0(ngrid, nlayer)           ! Entrainment

trunk/LMDZ.GENERIC/libf/phystd/thermcell_plume.F90

@@ -113 +113 @@
       REAL zltdwn                               !
       REAL zltup                                ! useless here
-      REAL zbuoybis(ngrid,klev)
       
       LOGICAL active(ngrid)                     ! if the plume is active at ig,l (speed and incoming mass flux > 0 or l=lmin)
@@ -186 +185 @@
       
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! On pourrait n'appeler thermcell_alim que si la plume est active
+! AB : On pourrait n'appeler thermcell_alim que si la plume est active
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       CALL thermcell_alim(iflag_thermals_alim,ngrid,klev,ztv,d_temp,  &
-      &                            zlev,alim_star,lalim,lmin)
+      &                   zlev,alim_star,lalim,lmin)
       
 !==============================================================================
@@ -227 +226 @@
 ! AB : we decide here if the plume is still active or not. When the plume's
 !      first level is reached, we set active to "true". Otherwise, it is given
-!      by w2, f_star, alim_star and entr_star.
+!      by zw2, f_star, alim_star and entr_star.
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          DO ig=1,ngrid
@@ -240 +239 @@
          
          ztemp(:) = zpspsk(:,l) * ztla(:,l-1)
+         
          DO ig=1,ngrid
             CALL watersat(ztemp(ig), pplev(ig,l), zqsat(ig))
@@ -293 +293 @@
                &              + 0. * zbuoy(ig,l)
                
-               zbuoybis(ig,l) = RG * 0.5 * (                                  &
-               &                (ztva_est(ig,l) - ztv(ig,l+1)) / ztv(ig,l+1)  &
-               &              + (ztva_est(ig,l) - ztv(ig,l)) / ztv(ig,l) )
-               
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! AB : initial formulae
@@ -306 +302 @@
 !               w_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : own derivation for w_est (2010 formula, e=d=0)
+! AB : own derivation for w_est (Rio 2010 formula with e=d=0)
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                zw2fact = 2. * fact_epsilon * zdz
@@ -326 +322 @@
             IF (active(ig)) THEN
                
-               zw2m = w_est(ig,l+1)
                zdz = zlev(ig,l+1) - zlev(ig,l)
-               
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : The following test is added to avoid NaN when w_est becomes negative.
-!      zalpha value is of no importance because w_est > 0 means it will be the 
-!      last layer reached by the plume.
-!      Then entr will vanished and detr be set to compensate exactly the
-!      incoming mass flux in thermcell_flux.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-               IF (w_est(ig,l+1)>0.) THEN
-                  zalpha = f0(ig) * f_star(ig,l) / sqrt(w_est(ig,l+1)) / rhobarz(ig,l)
+               zalpha = f0(ig) * f_star(ig,l) / sqrt(w_est(ig,l+1)) / rhobarz(ig,l)
+               
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! AB : The next test is added to avoid a division by zero when w_est becomes
+!      negative.
+!      Indeed, entr and detr computed here are of no importance because w_est
+!      <= 0 means it will be the last layer reached by the plume and then they
+!      will be reset in thermcell_flux.
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+               IF (w_est(ig,l+1).eq.0.) THEN
+                  zw2m = 1.
                ELSE
-                  zalpha = 0.
+                  zw2m = w_est(ig,l+1)
                ENDIF
                
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : The following test is added to avoid a division by zero if there is no 
-!      water in the environment.
+! AB : The next test is added to avoid a division by zero if there is no water
+!      in the environment.
 !      In the case where there is no water in the env. but water in the plume
 !      (ascending from depth) we set the effect on detrainment equal to zero
@@ -351 +347 @@
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                IF (po(ig,l).lt.1.e-6) THEN
-!                  print *, 'WARNING: po=0 at level', l
+!                  print *, 'WARNING: po=0 in layer',l,'!'
 !                  print *, 'po,zqta', po(ig,l), zqta(ig,l-1)
                   zdqt(ig,l) = 0.0
@@ -424 +420 @@
          DO ig=1,ngrid
             IF (activetmp(ig)) THEN
-               ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+                         &  ! ztla is set to TP in plume (mixed)
-               &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))     &
-               &            /(f_star(ig,l+1)+detr_star(ig,l))
-               zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+                         &  ! zqta is set to qt in plume (mixed)
-               &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))       &
-               &            /(f_star(ig,l+1)+detr_star(ig,l))
+               ztla(ig,l) = (f_star(ig,l) * ztla(ig,l-1)                      &  ! ztla is set to TP in plume (mixed)
+               &          + (alim_star(ig,l) + entr_star(ig,l)) * zthl(ig,l)) &
+               &          / (f_star(ig,l+1) + detr_star(ig,l))
+               zqta(ig,l) = (f_star(ig,l) * zqta(ig,l-1) +                    &  ! zqta is set to qt in plume (mixed)
+               &          + (alim_star(ig,l) + entr_star(ig,l)) * po(ig,l))   &
+               &          / (f_star(ig,l+1) + detr_star(ig,l))
             ENDIF
          ENDDO
@@ -444 +440 @@
 ! Calcul de la vitesse verticale zw2 apres le melange
 !------------------------------------------------------------------------------
-         
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : f_star(l+1)<=0 implies activetmp=F then plume stops
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          
          DO ig=1,ngrid
@@ -462 +454 @@
                
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : initial formula (integrated)
+! AB : initial formula
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 !               zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
@@ -469 +461 @@
 !               zw2(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : own integrated derivation for zw2 (2010 formula)
+! AB : own derivation for zw2 (Rio 2010 formula)
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                zw2fact = 2. * (fact_epsilon * zdz  + entr_star(ig,l) / f_star(ig,l))
@@ -489 +481 @@
       
       DO ig=1,ngrid
-         alim_star_tot(ig)=0.
+         alim_star_tot(ig) = 0.
       ENDDO