Changeset 2420 for LMDZ5/trunk/libf/phylmd/cv3p1_closure.F90

Timestamp:

Jan 5, 2016, 4:37:49 PM (8 years ago)

Author:

crio

Message:

Nouvelle option d'epluchage de l'ascendance adiabatique dans le schema d'Emanuel: epluchage fonction de B/w2 au lieu de w. S'active avec iflag_mix_adiab=1 (valeur par defaut iflag_mix_adiab=0). Fonctionne avec iflag_mix=0 et iflag_mix=1.
Correction de bugs dans le schema de convection pour le calcul de inb, cape et buoy (sous le meme flag pour l'instant).
New option for the erosion of the adiabatic ascent in the Emanuel scheme: erosion function of B/w2 instead of w. Activated by iflag_mix_adiab=1 (standard value iflag_mix_adiab=0). Should work with iflag_mix=0 and iflag_mix=1.
Various bug corrections in the convection scheme for the computation of inb, cape, buoy (protected by the same flag for now).

File:

• LMDZ5/trunk/libf/phylmd/cv3p1_closure.F90 (modified) (5 diffs)

LMDZ5/trunk/libf/phylmd/cv3p1_closure.F90

@@ -58 +58 @@
   INTEGER il, i, j, k, icbmax, i0(nloc), klfc(nloc)
   REAL deltap, fac, w, amu
-  REAL rhodp
+  REAL rhodp, dz
   REAL pbmxup
   REAL dtmin(nloc, nd), sigold(nloc, nd)
@@ -79 +79 @@
   REAL term1, term2, term3
   REAL alp2(nloc) ! Alp with offset
+  !CR: variables for new erosion of adiabiatic ascent
+  REAL mad(nloc, nd), me(nloc, nd), betalim(nloc, nd), beta_coef(nloc, nd)
+  REAL med(nloc, nd), md(nloc,nd)
+  REAL coef_peel
+  PARAMETER (coef_peel=0.25)
 
   REAL sigmax
@@ -110 +115 @@
     END DO
   END DO
+
+!CR: initializations for erosion of adiabatic ascent 
+  DO k = 1,nl
+    DO il = 1, ncum
+        mad(il,k)=0.
+        me(il,k)=0.
+        betalim(il,k)=1.
+        wlim(il,k)=0.
+    ENDDO
+  ENDDO
 
   ! -------------------------------------------------------
@@ -431 +446 @@
 
       IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
-
+        IF (iflag_mix_adiab.eq.1) THEN
+!CR:computation of cape from LCL: keep flag or to modify in all cases?
+        deltap = min(plcl(il), ph(il,k-1)) - min(plcl(il), ph(il,k))
+        ELSE
         deltap = min(pbase(il), ph(il,k-1)) - min(pbase(il), ph(il,k))
+        ENDIF
         cape(il) = cape(il) + rrd*buoy(il, k-1)*deltap/p(il, k-1)
         cape(il) = amax1(0.0, cape(il))
@@ -601 +620 @@
   IF (prt_level>=20) PRINT *, 'cv3p1_param apres w0_sig_M'
 
+!CR: new erosion of adiabatic ascent: modification of m
+!computation of the sum of ascending fluxes 
+  IF (iflag_mix_adiab.eq.1) THEN
+
+!Verification sum(me)=sum(m)
+  DO k = 1,nl+1
+    DO il = 1, ncum
+       md(il,k)=0.
+       med(il,k)=0.
+    ENDDO
+  ENDDO
+
+  DO k = nl,1,-1
+    DO il = 1, ncum
+           md(il,k)=md(il,k+1)+m(il,k+1)
+    ENDDO
+  ENDDO
+
+  DO k = nl,1,-1
+    DO il = 1, ncum
+        IF ((k>=(icb(il))) .AND. (k<=inb(il))) THEN
+           mad(il,k)=mad(il,k+1)+m(il,k+1)
+        ENDIF
+!        print*,"mad",il,k,mad(il,k)
+    ENDDO
+  ENDDO
+
+!CR: erosion of each adiabatic ascent during its ascent
+
+!Computation of erosion coefficient beta_coef
+  DO k = 1, nl
+    DO il = 1, ncum
+       IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (mlim(il,k).gt.0.)) THEN     
+!          print*,"beta_coef",il,k,icb(il),inb(il),buoy(il,k),tv(il,k),wlim(il,k),wlim(il,k+1)
+          beta_coef(il,k)=RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2
+       ELSE
+          beta_coef(il,k)=0.
+       ENDIF
+    ENDDO
+  ENDDO
+
+!  print*,"apres beta_coef"
+
+  DO k = 1, nl
+    DO il = 1, ncum
+
+      IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
+
+!        print*,"dz",il,k,tv(il, k-1)
+        dz = (ph(il,k-1)-ph(il,k))/(p(il, k-1)/(rrd*tv(il, k-1))*RG)
+        betalim(il,k)=betalim(il,k-1)*exp(-1.*beta_coef(il,k-1)*dz)
+!        betalim(il,k)=betalim(il,k-1)*exp(-RG*coef_peel*buoy(il,k-1)/tv(il,k-1)/5.**2*dz)
+!        print*,"me",il,k,mlim(il,k),buoy(il,k),wlim(il,k),mad(il,k)
+        dz = (ph(il,k)-ph(il,k+1))/(p(il, k)/(rrd*tv(il, k))*RG)
+!        me(il,k)=betalim(il,k)*(m(il,k)+RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz*mad(il,k))
+        me(il,k)=betalim(il,k)*(m(il,k)+beta_coef(il,k)*dz*mad(il,k))
+!        print*,"B/w2",il,k,RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz    
+      
+      END IF
+        
+!Modification of m
+      m(il,k)=me(il,k) 
+    END DO
+  END DO
+ 
+!  DO il = 1, ncum
+!     dz = (ph(il,icb(il))-ph(il,icb(il)+1))/(p(il, icb(il))/(rrd*tv(il, icb(il)))*RG)
+!     m(il,icb(il))=m(il,icb(il))+RG*coef_peel*buoy(il,icb(il))/tv(il,icb(il)) &
+!                  /((wlim(il,icb(il))+wlim(il,icb(il)+1))/2.)**2*dz*mad(il,icb(il))
+!     print*,"wlim(icb)",icb(il),wlim(il,icb(il)),m(il,icb(il))
+!  ENDDO
+
+!Verification sum(me)=sum(m)
+  DO k = nl,1,-1
+    DO il = 1, ncum
+           med(il,k)=med(il,k+1)+m(il,k+1)
+!           print*,"somme(me),somme(m)",il,k,icb(il),med(il,k),md(il,k),me(il,k),m(il,k),wlim(il,k)
+    ENDDO
+  ENDDO
+
+
+  ENDIF !(iflag_mix_adiab)
+!RC
+
+
+
   ! c 3. Compute final cloud base mass flux and set iflag to 3 if
   ! c    cloud base mass flux is exceedingly small and is decreasing (i.e. if