Changeset 2397

Timestamp:

Nov 19, 2015, 12:07:15 AM (9 years ago)

Author:

jyg

Message:

Implementation of ice thermodynamics in
cv3p_mixing.F90 . Needs debugging.

Location:

LMDZ5/trunk/libf/phylmd

Files:

• cv3p_mixing.F90 (modified) (10 diffs)
• cva_driver.F90 (modified) (1 diff)

LMDZ5/trunk/libf/phylmd/cv3p_mixing.F90

@@ -1 +1 @@
 SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
-                       ph, t, rr, rs, u, v, tra, h, lv, qnk, &
+                       ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qnk, &
                        unk, vnk, hp, tv, tvp, ep, clw, sig, &
                        Ment, Qent, hent, uent, vent, nent, &
@@ -20 +20 @@
   include "cv3param.h"
   include "YOMCST2.h"
+  include "cvflag.h"
 
 !inputs:
@@ -32 +33 @@
   REAL, DIMENSION (nloc, nd, ntra), INTENT (IN)      :: tra ! input of convect3
   REAL, DIMENSION (nloc, na), INTENT (IN)            :: lv
+  REAL, DIMENSION (nloc, na), INTENT (IN)            :: lf
+  REAL, DIMENSION (nloc, na), INTENT (IN)            :: frac !ice fraction in condensate
   REAL, DIMENSION (nloc, na), INTENT (IN)            :: h  !liquid water static energy of environment
   REAL, DIMENSION (nloc, na), INTENT (IN)            :: hp !liquid water static energy of air shed from adiab. asc.
@@ -51 +54 @@
   INTEGER i, j, k, il, im, jm
   INTEGER num1, num2
-  REAL                               :: rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
+  REAL                               :: rti, bf2, anum, denom, dei, altem, cwat, stemp
   REAL                               :: alt, delp, delm
   REAL, DIMENSION (nloc)             :: Qmixmax, Rmixmax, sqmrmax
@@ -60 +63 @@
   REAL, DIMENSION (nloc)             :: Smid, Sjmin, Sjmax
   REAL, DIMENSION (nloc)             :: Sbef, sup, smin
-!jyg  REAL, DIMENSION (nloc)             :: ASij, smax, Scrit
   REAL, DIMENSION (nloc)             :: ASij, ASij_inv, smax, Scrit
   REAL, DIMENSION (nloc, nd, nd)     :: Sij
   REAL, DIMENSION (nloc, nd)         :: csum
   REAL                               :: awat
+  REAL                               :: cpm        !Mixed draught heat capacity
+  REAL                               :: Tm         !Mixed draught temperature
   LOGICAL, DIMENSION (nloc)          :: lwork
 
@@ -165 +169 @@
           rti = qnk(il) - ep(il, i)*clw(il, i)
           bf2 = 1. + lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
+!jyg(from aj)<
+          IF (cvflag_ice) THEN
+! print*,cvflag_ice,'cvflag_ice dans do 700'
+            IF (t(il,j)<=263.15) THEN
+              bf2 = 1. + (lf(il,j)+lv(il,j))*(lv(il,j)+frac(il,j)* &
+                   lf(il,j))*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
+            END IF
+          END IF
+!>jyg
           anum = h(il, j) - hp(il, i) + (cpv-cpd)*t(il, j)*(rti-rr(il,j))
           denom = h(il, i) - hp(il, i) + (cpd-cpv)*(rr(il,i)-rti)*t(il, j)
@@ -176 +189 @@
           stemp = Sij(il, i, j)
           IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
-            anum = anum - lv(il, j)*(rti-rs(il,j)-cwat*bf2)
-            denom = denom + lv(il, j)*(rr(il,i)-rti)
+!jyg(from aj)<
+            IF (cvflag_ice) THEN
+              anum = anum - (lv(il,j)+frac(il,j)*lf(il,j))*(rti-rs(il,j)-cwat*bf2)
+              denom = denom + (lv(il,j)+frac(il,j)*lf(il,j))*(rr(il,i)-rti)
+            ELSE
+              anum = anum - lv(il, j)*(rti-rs(il,j)-cwat*bf2)
+              denom = denom + lv(il, j)*(rr(il,i)-rti)
+            END IF
+!>jyg
             IF (abs(denom)<0.01) denom = 0.01
             Sij(il, i, j) = anum/denom
@@ -299 +319 @@
         lwork(il) = (nent(il,i)/=0)
         rti = qnk(il) - ep(il, i)*clw(il, i)
-        anum = h(il, i) - hp(il, i) - lv(il, i)*(rti-rs(il,i)) + &
-               (cpv-cpd)*t(il, i)*(rti-rr(il,i))
-        denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-rti) + &
-                (cpd-cpv)*t(il, i)*(rr(il,i)-rti)
+!jyg<
+        IF (cvflag_ice) THEN
+
+          anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)*lf(il,i))* &
+                       (rti-rs(il,i)) + (cpv-cpd)*t(il, i)*(rti-rr(il,i))
+          denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)*lf(il,i))* &
+                       (rr(il,i)-rti) + (cpd-cpv)*t(il, i)*(rr(il,i)-rti)
+        ELSE
+
+          anum = h(il, i) - hp(il, i) - lv(il, i)*(rti-rs(il,i)) + &
+                       (cpv-cpd)*t(il, i)*(rti-rr(il,i))
+          denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-rti) + &
+                       (cpd-cpv)*t(il, i)*(rr(il,i)-rti)
+        END IF
+!>jyg
         IF (abs(denom)<0.01) denom = 0.01
         Scrit(il) = min(anum/denom, 1.)
@@ -452 +483 @@
             hent(il, i, j) = (1.-Sigij(il,i,j))*hp(il, i) + Sigij(il, i, j)*h(il, i)
 
+!jyg<
+!            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
+!            elij(il, i, j) = elij(il, i, j) + &
+!                             ((h(il,j)-hent(il,i,j))*rs(il,j)*lv(il,j) / &
+!                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
+!            elij(il, i, j) = elij(il, i, j) / &
+!                             (1.+lv(il,j)*lv(il,j)*rs(il,j) / &
+!                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
+!
+!       Computation of condensate amount Elij, taking into account the ice fraction frac
+!       Warning : the same saturation humidity rs is used over both liquid water and ice; this
+!                 should be corrected.
+!
+!  Heat capacity of mixed draught
+    cpm = cpd+Qent(il,i,j)*(cpv-cpd)
+!
+    IF (cvflag_ice .and. frac(il,j) .gt. 0.) THEN
             elij(il, i, j) = Qent(il, i, j) - rs(il, j)
             elij(il, i, j) = elij(il, i, j) + &
-                             ((h(il,j)-hent(il,i,j))*rs(il,j)*lv(il,j) / &
-                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
+                             (h(il,j)-hent(il,i,j)+(cpv-cpd)*(Qent(il,i,j)-rr(il,j))*t(il,j))* &
+                             rs(il,j)*lv(il,j) / (cpm*rrv*t(il,j)*t(il,j))
+            elij(il, i, j) = elij(il, i, j) / &
+                             (1.+(lv(il,j)+frac(il,j)*lf(il,j))*lv(il,j)*rs(il,j) / &
+                              (cpm*rrv*t(il,j)*t(il,j)))
+    ELSE
+            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
+            elij(il, i, j) = elij(il, i, j) + &
+                             (h(il,j)-hent(il,i,j)+(cpv-cpd)*(Qent(il,i,j)-rr(il,j))*t(il,j))* &
+                             rs(il,j)*lv(il,j) / (cpm*rrv*t(il,j)*t(il,j))
             elij(il, i, j) = elij(il, i, j) / &
                              (1.+lv(il,j)*lv(il,j)*rs(il,j) / &
-                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
-
+                              (cpm*rrv*t(il,j)*t(il,j)))
+    ENDIF
+!>jyg
             elij(il, i, j) = max(elij(il,i,j), 0.)
 
@@ -474 +531 @@
 ! :         t(il,j))
 
-            hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*t(il,j))*awat
+!jyg<
+!            hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*t(il,j))*awat
+! Mixed draught temperature at level j
+    IF (cvflag_ice .and. frac(il,j) .gt. 0.) THEN
+          Tm = t(il,j) + (Qent(il,i,j)-elij(il,i,j)-rs(il,j))*rrv*t(il,j)*t(il,j)/(lv(il,j)*rs(il,j))
+          hent(il, i, j) = hent(il, i, j) + (lv(il,j)+frac(il,j)*lf(il,j)+(cpd-cpv)*Tm)*awat
+    ELSE
+          Tm = t(il,j) + (Qent(il,i,j)-elij(il,i,j)-rs(il,j))*rrv*t(il,j)*t(il,j)/(lv(il,j)*rs(il,j))
+          hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*Tm)*awat
+    ENDIF
+!>jyg
+
 !IM 301008 end
 

LMDZ5/trunk/libf/phylmd/cva_driver.F90

@@ -900 +900 @@
         CALL zilch(supmax, nloc*klev)
         CALL cv3p_mixing(nloc, ncum, nd, nd, ntra, icb, nk, inb, &           ! na->nd
-                         ph, t, q, qs, u, v, tra, h, lv, qnk, &
+                         ph, t, q, qs, u, v, tra, h, lv, lf, frac, qnk, &
                          unk, vnk, hp, tv, tvp, ep, clw, sig, &
                          ment, qent, hent, uent, vent, nent, &