Changeset 4669 for LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd/lscp_tools_mod.F90

Timestamp:

Sep 4, 2023, 10:17:16 AM (10 months ago)

Author:

Laurent Fairhead

Message:

Merged with trunk revision 4586 corresponding to june 2023 testing

Location:

LMDZ6/branches/LMDZ_cdrag_LSCE

Files:

• . (modified) (1 prop)
• libf/phylmd/lscp_tools_mod.F90 (modified) (10 diffs)

LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd/lscp_tools_mod.F90

@@ -16 +16 @@
     ! 3212–3234. https://doi.org/10.1029/2019MS001642
     
-    
+    use lscp_ini_mod, only: iflag_vice, ffallv_con, ffallv_lsc
+    use lscp_ini_mod, only: cice_velo, dice_velo 
+
     IMPLICIT NONE
-
-    INCLUDE "nuage.h"
-    INCLUDE "fisrtilp.h"
 
     INTEGER, INTENT(IN) :: klon
@@ -33 +32 @@
 
     INTEGER i
-    REAL logvm,iwcg,tempc,phpa,cvel,dvel,fallv_tun
+    REAL logvm,iwcg,tempc,phpa,fallv_tun
     REAL m2ice, m2snow, vmice, vmsnow
     REAL aice, bice, asnow, bsnow
@@ -62 +61 @@
         
         velo(i)=fallv_tun*velo(i)/100.0 ! from cm/s to m/s
-        dvel=0.2
-        cvel=fallv_tun*65.0*(rho(i)**0.2)*(150./phpa)**0.15
 
     ELSE IF (iflag_vice .EQ. 2) THEN
@@ -94 +91 @@
         vmsnow=vmsnow*((1000./phpa)**0.35)
         velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s
-        dvel=0.2
-        cvel=velo(i)/((iwcg/1000.*rho(i))**dvel)
 
     ELSE
         ! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990
-        velo(i) = fallv_tun*3.29/2.0*((iwcg/1000.)**0.16)
-        dvel=0.16
-        cvel=fallv_tun*3.29/2.0*(rho(i)**0.16)
+        velo(i) = fallv_tun*cice_velo*((iwcg/1000.)**dice_velo)
     ENDIF
     ENDDO
@@ -109 +102 @@
 
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-SUBROUTINE ICEFRAC_LSCP(klon, temp, sig, icefrac, dicefracdT)
+SUBROUTINE ICEFRAC_LSCP(klon, temp, iflag_ice_thermo, distcltop, icefrac, dicefracdT)
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   
@@ -120 +113 @@
 
     USE print_control_mod, ONLY: lunout, prt_level
+    USE lscp_ini_mod, ONLY: t_glace_min, t_glace_max, exposant_glace, iflag_t_glace
+    USE lscp_ini_mod, ONLY : RTT, dist_liq
 
     IMPLICIT NONE
 
-
-    INCLUDE "YOMCST.h"
-    INCLUDE "nuage.h"
-    INCLUDE "clesphys.h"
-
-
-  ! nuage.h contains:
-  ! t_glace_min: if T < Tmin, the cloud is only made of water ice
-  ! t_glace_max: if T > Tmax, the cloud is only made of liquid water
-  ! exposant_glace: controls the sharpness of the transition
 
     INTEGER, INTENT(IN)                 :: klon       ! number of horizontal grid points
     REAL, INTENT(IN), DIMENSION(klon)   :: temp       ! temperature
-    REAL, INTENT(IN), DIMENSION(klon)   :: sig
+    REAL, INTENT(IN), DIMENSION(klon)   :: distcltop       ! distance to cloud top
+    INTEGER, INTENT(IN)                 :: iflag_ice_thermo
     REAL, INTENT(OUT), DIMENSION(klon)  :: icefrac
     REAL, INTENT(OUT), DIMENSION(klon)  :: dicefracdT
@@ -142 +128 @@
 
     INTEGER i
-    REAL    sig0,www,tmin_tmp,liqfrac_tmp
+    REAL    liqfrac_tmp, dicefrac_tmp
     REAL    Dv, denomdep,beta,qsi,dqsidt
-    INTEGER exposant_glace_old
-    REAL t_glace_min_old
     LOGICAL ice_thermo
 
-    sig0=0.8
-    t_glace_min_old = RTT - 15.0
-    ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
-    IF (ice_thermo) THEN
-      exposant_glace_old = 2
-    ELSE
-      exposant_glace_old = 6
+    CHARACTER (len = 20) :: modname = 'lscp_tools'
+    CHARACTER (len = 80) :: abort_message
+
+    IF ((iflag_t_glace.LT.2) .OR. (iflag_t_glace.GT.6)) THEN
+       abort_message = 'lscp cannot be used if iflag_t_glace<2 or >6'
+       CALL abort_physic(modname,abort_message,1)
     ENDIF
 
-
-! calculation of icefrac and dicefrac/dT
+    IF (.NOT.((iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3))) THEN
+       abort_message = 'lscp cannot be used without ice thermodynamics'
+       CALL abort_physic(modname,abort_message,1)
+    ENDIF
+
 
     DO i=1,klon
-
-    IF (iflag_t_glace.EQ.1) THEN
-            ! Transition to ice close to surface for T<Tmax
-            ! w=1 at the surface and 0 for sig < sig0
-            www=(max(sig(i)-sig0,0.))/(1.-sig0)
-    ELSEIF (iflag_t_glace.GE.2) THEN
-            ! No convertion to ice close to surface
-            www = 0.
-    ENDIF
-     
-    tmin_tmp=www*t_glace_max+(1.-www)*t_glace_min
-    liqfrac_tmp=  (temp(i)-tmin_tmp) / (t_glace_max-tmin_tmp)
-    liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
-       
-    IF (iflag_t_glace.GE.3) THEN
-            icefrac(i) = 1.0-liqfrac_tmp**exposant_glace
-            IF ((icefrac(i) .GT.0.) .AND. (liqfrac_tmp .GT. 0)) THEN
-                 dicefracdT(i)= exposant_glace * ((liqfrac_tmp)**(exposant_glace-1.)) &
-                           / (t_glace_min - t_glace_max)
-            ELSE
-
-                 dicefracdT(i)=0.
-            ENDIF
-
-    ELSE
+ 
+        ! old function with sole dependence upon temperature
+        IF (iflag_t_glace .EQ. 2) THEN
+            liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
+            liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
             icefrac(i) = (1.0-liqfrac_tmp)**exposant_glace
             IF (icefrac(i) .GT.0.) THEN
@@ -196 +162 @@
             ENDIF
 
-    ENDIF
-    
-    ENDDO
-
-
-    RETURN
-
+        ENDIF
+
+        ! function of temperature used in CMIP6 physics
+        IF (iflag_t_glace .EQ. 3) THEN
+            liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
+            liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
+            icefrac(i) = 1.0-liqfrac_tmp**exposant_glace
+            IF ((icefrac(i) .GT.0.) .AND. (liqfrac_tmp .GT. 0.)) THEN
+                dicefracdT(i)= exposant_glace * ((liqfrac_tmp)**(exposant_glace-1.)) &
+                          / (t_glace_min - t_glace_max)
+            ELSE
+                dicefracdT(i)=0.
+            ENDIF
+        ENDIF
+
+        ! for iflag_t_glace .GE. 4, the liquid fraction depends upon temperature at cloud top
+        ! and then decreases with decreasing height 
+
+        !with linear function of temperature at cloud top
+        IF (iflag_t_glace .EQ. 4) THEN  
+                liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
+                liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
+                icefrac(i)    =  MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
+                dicefrac_tmp = - temp(i)/(t_glace_max-t_glace_min)
+                dicefracdT(i) = dicefrac_tmp*exp(-distcltop(i)/dist_liq)
+                IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
+                        dicefracdT(i) = 0.
+                ENDIF
+        ENDIF
+
+        ! with CMIP6 function of temperature at cloud top
+        IF (iflag_t_glace .EQ. 5) THEN 
+                liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
+                liqfrac_tmp =  MIN(MAX(liqfrac_tmp,0.0),1.0)
+                liqfrac_tmp = liqfrac_tmp**exposant_glace
+                icefrac(i)  =  MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
+                IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
+                        dicefracdT(i) = 0.
+                ELSE
+                        dicefracdT(i) = exposant_glace*((liqfrac_tmp)**(exposant_glace-1.))/(t_glace_min- t_glace_max) &
+                                        *exp(-distcltop(i)/dist_liq)
+                ENDIF
+        ENDIF
+
+        ! with modified function of temperature at cloud top 
+        ! to get largere values around 260 K, works well with t_glace_min = 241K
+        IF (iflag_t_glace .EQ. 6) THEN 
+                IF (temp(i) .GT. t_glace_max) THEN
+                        liqfrac_tmp = 1.
+                ELSE
+                        liqfrac_tmp = -((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))**2+1.
+                ENDIF
+                liqfrac_tmp  = MIN(MAX(liqfrac_tmp,0.0),1.0)
+                icefrac(i)   = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)        
+                IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
+                        dicefracdT(i) = 0.
+                ELSE
+                        dicefracdT(i) = 2*((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))/(t_glace_max-t_glace_min) &
+                                  *exp(-distcltop(i)/dist_liq)
+                ENDIF
+        ENDIF
+
+
+
+     ENDDO ! klon
+ 
+     RETURN
+ 
 END SUBROUTINE ICEFRAC_LSCP
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -276 +303 @@
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
+    use lscp_ini_mod, only: iflag_gammasat, t_glace_min, RTT
 
     IMPLICIT NONE
 
-    include "YOMCST.h"
-    include "YOETHF.h"
-    include "FCTTRE.h"
-    include "nuage.h" 
 
     INTEGER, INTENT(IN) :: klon                       ! number of horizontal grid points
@@ -348 +372 @@
 
 END SUBROUTINE CALC_GAMMASAT
-
-
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D)
+!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   
+   USE lscp_ini_mod, ONLY : rd,rg,tresh_cl
+
+   IMPLICIT NONE
+   
+   INTEGER, INTENT(IN) :: klon,klev                !number of horizontal and vertical grid points
+   INTEGER, INTENT(IN) :: k                        ! vertical index
+   REAL, INTENT(IN), DIMENSION(klon,klev) :: temp  ! temperature in K
+   REAL, INTENT(IN), DIMENSION(klon,klev) :: pplay ! pressure middle layer in Pa
+   REAL, INTENT(IN), DIMENSION(klon,klev+1) :: paprs ! pressure interfaces in Pa
+   REAL, INTENT(IN), DIMENSION(klon,klev) :: rneb  ! cloud fraction
+
+   REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D  ! distance from cloud top
+
+   REAL dzlay(klon,klev)
+   REAL zlay(klon,klev)
+   REAL dzinterf
+   INTEGER i,k_top, kvert
+   LOGICAL bool_cl
+
+
+   DO i=1,klon
+         ! Initialization height middle of first layer
+          dzlay(i,1) = Rd * temp(i,1) / rg * log(paprs(i,1)/paprs(i,2))
+          zlay(i,1) = dzlay(i,1)/2
+
+          DO kvert=2,klev
+                 IF (kvert.EQ.klev) THEN
+                       dzlay(i,kvert) = 2*(rd * temp(i,kvert) / rg * log(paprs(i,kvert)/pplay(i,kvert)))
+                 ELSE
+                       dzlay(i,kvert) = rd * temp(i,kvert) / rg * log(paprs(i,kvert)/paprs(i,kvert+1))
+                 ENDIF
+                       dzinterf       = rd * temp(i,kvert) / rg * log(pplay(i,kvert-1)/pplay(i,kvert))
+                       zlay(i,kvert)  = zlay(i,kvert-1) + dzinterf
+           ENDDO
+   ENDDO
+   
+   k_top = k
+   DO i=1,klon
+          IF (rneb(i,k) .LE. tresh_cl) THEN
+                 bool_cl = .FALSE.
+          ELSE
+                 bool_cl = .TRUE.
+          ENDIF
+
+          DO WHILE ((bool_cl) .AND. (k_top .LE. klev))
+          ! find cloud top
+                IF (rneb(i,k_top) .GT. tresh_cl) THEN
+                      k_top = k_top + 1
+                ELSE
+                      bool_cl = .FALSE.
+                      k_top   = k_top - 1
+                ENDIF
+          ENDDO
+          k_top=min(k_top,klev)
+
+          !dist to top is dist between current layer and layer of cloud top (from middle to middle) + dist middle to
+          !interf for layer of cloud top
+          distcltop1D(i) = zlay(i,k_top) - zlay(i,k) + dzlay(i,k_top)/2
+   ENDDO ! klon
+
+END SUBROUTINE DISTANCE_TO_CLOUD_TOP
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++