Changeset 2143 for trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90

Timestamp:

Jun 20, 2019, 4:11:54 PM (5 years ago)

Author:

aboissinot

Message:

Update the thermal plume model:

• check formulae consistency between thermcell_plume and thercell_closure
• compute correctly thermal plume height
• fix alimentation computation in the first unstable layer

File:

• trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90 (modified) (12 diffs)

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

@@ -49 +49 @@
       
       USE thermcell_mod
-      USE tracer_h, ONLY: igcm_h2o_vap
-      USE print_control_mod, ONLY: lunout, prt_level
+      USE print_control_mod, ONLY: prt_level
       
       IMPLICIT NONE
@@ -96 +95 @@
       INTEGER lmix(ngrid)                       ! Layer in which plume vertical speed is maximal
       INTEGER lmin(ngrid)                       ! First unstable layer
+      
+      REAL zmix(ngrid)                          ! Altitude of maximal vertical speed
+      REAL zmax(ngrid)                          ! Maximal altitudes where plumes are active
+      REAL zmin(ngrid)                          ! Minimal altitudes where plumes are active
       
       REAL zlay(ngrid,nlay)                     ! Layers altitudes
@@ -124 +127 @@
       REAL zqa(ngrid,nlay,nq)                   ! q    plume (ql=0, qv=qt)
       
+      REAL f_star(ngrid,nlay+1)                 ! Normalized mass flux
+      REAL entr_star(ngrid,nlay)                ! Normalized entrainment (E* = e* dz)
+      REAL detr_star(ngrid,nlay)                ! Normalized detrainment (D* = d* dz)
+      
+      REAL fm(ngrid,nlay+1)                     ! Mass flux
+      REAL entr(ngrid,nlay)                     ! Entrainment (E = e dz)
+      REAL detr(ngrid,nlay)                     ! Detrainment (D = d dz)
+      
+      REAL f(ngrid)                             ! Mass flux norm
+      REAL lambda                               ! Time relaxation coefficent
+      REAL fraca(ngrid,nlay+1)                  ! Updraft fraction
       REAL linter(ngrid)                        ! Level (continuous) of maximal vertical speed
-      REAL zmix(ngrid)                          ! Altitude of maximal vertical speed
-      REAL zmax(ngrid)                          ! Plume height
       REAL wmax(ngrid)                          ! Maximal vertical speed
       REAL zw2(ngrid,nlay+1)                    ! Plume vertical speed
-      
-      REAL fraca(ngrid,nlay+1)                  ! Updraft fraction
-      
-      REAL f_star(ngrid,nlay+1)                 ! Normalized mass flux
-      REAL entr_star(ngrid,nlay)                ! Normalized entrainment
-      REAL detr_star(ngrid,nlay)                ! Normalized detrainment
-      
-      REAL f(ngrid)                             ! Mass flux norm
-      REAL fm(ngrid,nlay+1)                     ! Mass flux
-      REAL entr(ngrid,nlay)                     ! Entrainment
-      REAL detr(ngrid,nlay)                     ! Detrainment
-      
-      REAL lambda                               ! Time relaxation coefficent
-      
       REAL zdthladj(ngrid,nlay)                 ! Potential temperature variations
       REAL dummy(ngrid,nlay)                    ! Dummy argument for thermcell_dq()
-      
-      CHARACTER (len=20) :: modname='thermcell_main'
-      CHARACTER (len=80) :: abort_message
       
 !===============================================================================
@@ -178 +173 @@
       pdqadj(:,:,:) = 0.0
       
-! AB: Careful, hard-coded value from Earth tuned version of the thermal plume model!
       DO ig=1,ngrid
-         f0(ig) = max(f0(ig), 1.e-2)
-      ENDDO
-      
-      IF (prt_level.ge.20) then
-         DO ig=1,ngrid
-            print *, 'ig,f0', ig, f0(ig)
-         ENDDO
-      ENDIF
+! AB: Careful: Hard-coded value from Earth version!
+!         f0(ig) = max(f0(ig), 1.e-2)
+! AB: No pescribed minimal value for f0
+         f0(ig) = max(f0(ig), 0.)
+      ENDDO
       
 !===============================================================================
@@ -212 +203 @@
       
       DO l=1,nlay
-         zlay(:,l) = pphi(:,l)/RG
+         zlay(:,l) = pphi(:,l) / RG
       ENDDO
       
@@ -222 +213 @@
       
       IF (prt_level.ge.10) THEN
-         write(lunout,*) 'WARNING: thermcell_main rhobarz(:,1)=rho(:,1)'
+         print *, 'WARNING: density in the first layer is equal to density at the first level!'
+         print *, 'rhobarz(:,1)', rhobarz(:,1)
+         print *, 'rho(:,1)    ', rho(:,1)
       ENDIF
       
@@ -319 +312 @@
       &                    zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk,        &
       &                    detr_star,entr_star,f_star,                        &
-      &                    ztva,zhla,zqla,zqta,zta,                           &
-      &                    zw2,zqsa,lmix,lmin)
-            
+      &                    ztva,zhla,zqla,zqta,zta,zqsa,                      &
+      &                    zw2,lmix,lmin)
+      
 !-------------------------------------------------------------------------------
 ! Thermal plumes characteristics: zmax, zmix, wmax
@@ -327 +320 @@
       
 ! AB: Careful, zw2 became its square root in thermcell_height!
-      CALL thermcell_height(ngrid,nlay,lmin,linter,lmix,zw2,                  &
-      &                     zlev,lmax,zmax,zmix,wmax,f_star)
+      CALL thermcell_height(ngrid,nlay,lmin,linter,lmix,lmax,zw2,             &
+      &                     zlev,zmin,zmix,zmax,wmax,f_star)
       
 !===============================================================================
@@ -339 +332 @@
       
       CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev,                   &
-      &                      lmax,entr_star,zmax,wmax,f)
+      &                      lmax,entr_star,zmin,zmax,wmax,f)
       
       IF (tau_thermals>1.) THEN
@@ -352 +345 @@
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       IF (.not. (f0(1).ge.0.) ) THEN
-         abort_message = '.not. (f0(1).ge.0.)'
+         print *, 'ERROR: mass flux norm is not positive!'
          print *, 'f0 =', f0(1)
-         CALL abort_physic(modname,abort_message,1)
+         CALL abort
       ENDIF
       
@@ -393 +386 @@
       DO l=2,nlay
          DO ig=1,ngrid
-            IF (zw2(ig,l).gt.0.) THEN
+            IF (zw2(ig,l) > 0.) THEN
                fraca(ig,l) = fm(ig,l) / (rhobarz(ig,l) * zw2(ig,l))
             ELSE
@@ -433 +426 @@
       CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
       &                 zu,pduadj,zua,lmin)
-      
+       
       CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
       &                 zv,pdvadj,zva,lmin)
       
+!      CALL thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,fraca,    &
+!      &                  zmax,zmin,pu,pv,pduadj,pdvadj,zua,zva)
       
 RETURN