Changeset 2143 for trunk/LMDZ.GENERIC/libf/phystd/thermcell_closure.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_closure.F90 (modified) (9 diffs)

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

@@ -3 +3 @@
 !
 SUBROUTINE thermcell_closure(ngrid,nlay,ptimestep,rho,zlev,                   &
-                             lmax,alim_star,zmax,wmax,f)
+                             lmax,alim_star,zmin,zmax,wmax,f)
+      
       
 !===============================================================================
@@ -13 +14 @@
 ! 2. On ne garde qu'une version des couples wmax,zmax et wmax_sec,zmax_sec
 ! l'idee etant que le choix se fasse a l'appel de thermcell_closure
-! 3. Vectorisation en mettant les boucles en l l'exterieur avec des if
+! 3. Vectorisation en mettant les boucles en l a l'exterieur avec des if
 !===============================================================================
       
@@ -35 +36 @@
       REAL zlev(ngrid,nlay)
       REAL alim_star(ngrid,nlay)
+      REAL zmin(ngrid)
       REAL zmax(ngrid)
       REAL wmax(ngrid)
@@ -51 +53 @@
       REAL alim_star_tot(ngrid)
       REAL alim_star2(ngrid)
-      REAL zdenom(ngrid)
+      REAL plume_height(ngrid)
       
 !===============================================================================
@@ -64 +66 @@
       llmax = 1
       
+      DO ig=1,ngrid
+         plume_height(ig) = zmax(ig) - zmin(ig)
+      ENDDO
+      
 !===============================================================================
 ! Closure
@@ -73 +79 @@
       
       DO ig=1,ngrid
-         IF (lmax(ig).GT.llmax) THEN
+         IF (lmax(ig) > llmax) THEN
             llmax = lmax(ig)
          ENDIF
@@ -86 +92 @@
             IF (l < lmax(ig)) THEN
                alim_star2(ig) = alim_star2(ig) + alim_star(ig,l)**2           &
-               &              / (rho(ig,l) * (zlev(ig,l+1) - zlev(ig,l)))
+               &              / (rho(ig,l) * (zlev(ig,l+1) - zlev(ig,l)))        ! => intergration is ok because alim_star = a* dz
                alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig,l)
             ENDIF
@@ -93 +99 @@
       
       DO ig=1,ngrid
-         IF (alim_star2(ig).GT.0.) THEN
-            f(ig) = wmax(ig) * alim_star_tot(ig)                              &
-            &     / (max(1.,zmax(ig)) * r_aspect_thermals * alim_star2(ig))
+         IF ((alim_star2(ig) > 0.).and.(plume_height(ig) > 0.)) THEN
+            f(ig) = wmax(ig) * alim_star_tot(ig)                              &  ! => normalization is ok
+            &     / (plume_height(ig) * r_aspect_thermals * alim_star2(ig))
          ELSE
             f(ig) = 0.
@@ -101 +107 @@
       ENDDO
       
+      print *, 'A*2 ', alim_star2(1)
+      print *, 'A*  ', alim_star_tot(1)
+      print *, 'H   ', plume_height(1)
+      print *, 'wmax', wmax(1)
       
 RETURN