Changeset 5626 for LMDZ6

Timestamp:

Apr 23, 2025, 10:51:22 AM (7 weeks ago)

Author:

aborella

Message:

Corrections coupling with convective clouds

Location:

LMDZ6/branches/contrails/libf

Files:

• dynphy_lonlat/phylmd/etat0phys_netcdf.f90 (modified) (3 diffs)
• phylmd/create_etat0_unstruct_mod.f90 (modified) (1 diff)
• phylmd/dyn1d/old_lmdz1d.f90 (modified) (3 diffs)
• phylmd/dyn1d/scm.f90 (modified) (3 diffs)
• phylmd/lmdz_lscp.f90 (modified) (4 diffs)
• phylmd/lmdz_lscp_condensation.f90 (modified) (3 diffs)
• phylmd/phyetat0_mod.f90 (modified) (2 diffs)
• phylmd/phyredem.f90 (modified) (2 diffs)
• phylmd/phys_state_var_mod.F90 (modified) (3 diffs)
• phylmd/physiq_mod.F90 (modified) (4 diffs)

LMDZ6/branches/contrails/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.f90

@@ -44 +44 @@
           solsw, solswfdiff, radsol, t_ancien, wake_deltat, wake_s,  rain_fall, qsol, z0h, &
           sollw,sollwdown, rugoro, q_ancien, wake_deltaq, wake_pe, snow_fall, ratqs,w01, &
-    sig1, ftsol, cwcon, clwcon, fm_therm, wake_Cstar,  pctsrf,  entr_therm,radpas, f0,&
+    sig1, ftsol, clwcon, fm_therm, wake_Cstar,  pctsrf,  entr_therm,radpas, f0,&
     zmax0,fevap, rnebcon,falb_dir, falb_dif, wake_fip,    agesno,  detr_therm, pbl_tke,  &
     phys_state_var_init, ql_ancien, qs_ancien, prlw_ancien, prsw_ancien, &
@@ -50 +50 @@
     ale_bl, ale_bl_trig, alp_bl, &
     ale_wake, ale_bl_stat, AWAKE_S, &
-    cf_ancien, qvc_ancien, cfa_ancien, pcf_ancien, qva_ancien, qia_ancien
+    cf_ancien, qvc_ancien, qvcon, qccon, cfa_ancien, pcf_ancien, qva_ancien, qia_ancien
 
   USE comconst_mod, ONLY: pi, dtvr
@@ -243 +243 @@
   cf_ancien = 0.
   qvc_ancien = 0.
-  cwcon = 0.
+  qvcon = 0.
+  qccon = 0.
   cfa_ancien = 0.
   pcf_ancien = 0.

LMDZ6/branches/contrails/libf/phylmd/create_etat0_unstruct_mod.f90

@@ -260 +260 @@
     cf_ancien = 0.
     qvc_ancien = 0.
-    cwcon = 0.
+    qvcon = 0.
+    qccon = 0.
 
     wake_delta_pbl_TKE(:,:,:)=0

LMDZ6/branches/contrails/libf/phylmd/dyn1d/old_lmdz1d.f90

@@ -10 +10 @@
       USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar,getin
    USE phys_state_var_mod, ONLY : phys_state_var_init, phys_state_var_end, &
-       cwcon, clwcon, detr_therm, &
+       clwcon, detr_therm, &
        qsol, fevap, z0m, z0h, agesno, &
        du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
@@ -23 +23 @@
        zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, &
        ql_ancien, qs_ancien, qbs_ancien, &
-       cf_ancien, qvc_ancien, cfa_ancien, pcf_ancien, qva_ancien, qia_ancien, &
+       cf_ancien, qvc_ancien, qvcon, qccon, cfa_ancien, pcf_ancien, qva_ancien, qia_ancien, &
        prlw_ancien, prsw_ancien, prbsw_ancien, prw_ancien, &
        u10m,v10m,ale_wake,ale_bl_stat
@@ -874 +874 @@
           cf_ancien = 0.
           qvc_ancien = 0.
-          cwcon = 0.
+          qvcon = 0.
+          qccon = 0.
         ENDIF
         IF ( ok_plane_contrail ) THEN

LMDZ6/branches/contrails/libf/phylmd/dyn1d/scm.f90

@@ -6 +6 @@
    USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar,getin
    USE phys_state_var_mod, ONLY : phys_state_var_init, phys_state_var_end, &
-       cwcon, clwcon, detr_therm, &
+       clwcon, detr_therm, &
        qsol, fevap, z0m, z0h, agesno, &
        du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
@@ -19 +19 @@
        zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, &
        ql_ancien, qs_ancien, qbs_ancien, &
-       cf_ancien, qvc_ancien, cfa_ancien, pcf_ancien, qva_ancien, qia_ancien, &
+       cf_ancien, qvc_ancien, qvcon, qccon, cfa_ancien, pcf_ancien, qva_ancien, qia_ancien, &
        prlw_ancien, prsw_ancien, prbsw_ancien, prw_ancien, &
        u10m,v10m,ale_wake,ale_bl_stat, ratqs_inter_
@@ -616 +616 @@
           cf_ancien = 0.
           qvc_ancien = 0.
-          cwcon = 0.
+          qvcon = 0.
+          qccon = 0.
         ENDIF
         IF ( ok_plane_contrail ) THEN

LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90

@@ -156 +156 @@
   LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: cfcon_old       ! cloud fraction from deep convection from previous timestep [-]
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qvcon_old       ! in-cloud vapor specific humidity from deep convection from previous timestep [kg/kg]
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qccon_old       ! in-cloud condensed specific humidity from deep convection from previous timestep [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qvcon_old       ! in-cloud vapor specific humidity from deep convection from previous timestep [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qccon_old       ! in-cloud condensed specific humidity from deep convection from previous timestep [kg/kg]
   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: cfcon           ! cloud fraction from deep convection [-]
   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qvcon           ! in-cloud vapor specific humidity from deep convection [kg/kg]
@@ -683 +683 @@
                   zq_nodeep(i) = zq(i) - ( qvcon(i,k) + qccon(i,k) ) * cfcon(i,k)
 
-              IF ( ( cfcon(i,k) * qccon(i,k) ) .LT. ( cfcon_old(i,k) * qccon_old(i,k) ) ) THEN
+              IF ( cfcon(i,k) .LT. cfcon_old(i,k) ) THEN
                 !--If deep convection is weakening, we add the clouds that are not anymore
                 !--'in' deep convection to the advected clouds
-                cldfra_in(i) = cldfra_in(i) + MAX(0., cfcon_old(i,k) - cfcon(i,k))
-                qvc_in(i) = qvc_in(i) + qvcon_old(i,k) * MAX(0., cfcon_old(i,k) - cfcon(i,k))
-                qice_in(i) = qice_in(i) + ( qccon_old(i,k) * cfcon_old(i,k) &
-                                          - qccon(i,k) * cfcon(i,k) )
-              ELSEIF ( cfcon(i,k) .GT. cfcon_old(i,k) ) THEN
+                cldfra_in(i) = cldfra_in(i) + ( cfcon_old(i,k) - cfcon(i,k) )
+                qvc_in(i) = qvc_in(i) + qvcon_old(i,k) * ( cfcon_old(i,k) - cfcon(i,k) )
+                qice_in(i) = qice_in(i) + qccon_old(i,k) * ( cfcon_old(i,k) - cfcon(i,k) )
+              ELSE
                 !--Else if deep convection is strengthening, it consumes the existing cloud
                 !--fraction (which does not at this moment represent deep convection)
-                !--NB. if deep convection is strengthening while the fraction decreases,
-                !--clear sky water vapor will be transfered in priority
                 cldfra_in(i) = cldfra_in(i) * ( 1. &
                              - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) ) )
@@ -1231 +1228 @@
         qvc_seri(i,k) = qvc(i)
 
+        !--We keep convective clouds properties in memory, and account for
+        !--the sink of condensed water from precipitation
+        IF ( ptconv(i,k) ) THEN
+          qvcon_old(i,k) = qvcon(i,k)
+          qccon_old(i,k) = qccon(i,k) * zcond(i) / zoliq(i)
+        ELSE
+          qvcon_old(i,k) = 0.
+          qccon_old(i,k) = 0.
+        ENDIF
+
         !--Deep convection clouds properties are removed from radiative properties
         !--outputed from lscp (NB. rneb and radocond are only used for the radiative
@@ -1237 +1244 @@
         IF ( ptconv(i,k) .AND. pt_pron_clds(i) ) THEN
           rneb(i,k) = rneb(i,k) - cfcon(i,k)
-          radocond(i,k) = radocond(i,k) - qccon(i,k) * cfcon(i,k)
+          radocond(i,k) = radocond(i,k) - qccon_old(i,k) * cfcon(i,k)
         ENDIF
 

LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90

@@ -238 +238 @@
 ! for unadjusted clouds
 REAL :: qiceincld, qvapincld, qvapincld_new
+REAL :: qice_ratio
 !
 ! for deposition / sublimation
@@ -538 +539 @@
               !--Exact explicit integration (qvc exact, qice explicit)
               tauinv_depsub = depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub
+              qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime * tauinv_depsub )
             ELSE
               !--If the cloud is initially subsaturated
-              !--Exact explicit integration (qvc exact, qice explicit)
-              !--Same but depo_coef_cirrus = 1
+              !!--Exact explicit integration (qvc exact, qice explicit)
+              !!--Same but depo_coef_cirrus = 1
+              !tauinv_depsub = qiceincld**(1./3.) * kappa_depsub
+              !qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime * tauinv_depsub )
+              !--Exact explicit integration (qice exact, qvc explicit)
+              !--The barrier is set so that the resulting vapor in cloud
+              !--cannot be greater than qsat
+              !--qice_ratio is the ratio between the new ice content and
+              !--the old one, it is comprised between 0 and 1
               tauinv_depsub = qiceincld**(1./3.) * kappa_depsub
+              qice_ratio = tauinv_depsub * dtime / 1.5 / qiceinmix * ( qsat(i) - qvapincld )
+              !--The new vapor in the cloud is increased with the
+              !--sublimated ice
+              qvapincld_new = qvapincld + qiceincld * ( 1. - MAX(0., 1. - qice_ratio)**1.5 )
+              !--The new vapor in the cloud cannot be greater than qsat
+              qvapincld_new = MIN(qvapincld_new, qsat(i))
+              !--If all the ice is sublimated
+              IF ( qvapincld_new .GE. ( qvapincld + qiceincld ) ) qvapincld_new = 0.
             ENDIF ! qvapincld .GT. qsat
-            qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime * tauinv_depsub )
-            !--If all the ice is sublimated
-            IF ( qvapincld_new .GE. ( qvapincld + qiceincld ) ) qvapincld_new = 0.
           ELSE
             !--We keep the saturation adjustment hypothesis, and the vapor in the
@@ -836 +850 @@
           qiceinmix = ( qcld(i) - qvc(i) ) / cldfra(i) / ( 1. + clrfra_mix / cldfra_mix )
           tauinv_depsub = qiceinmix**(1./3.) * kappa_depsub
-          qvapinmix_lim = qsat(i) - qiceinmix / ( 1. - EXP( - dtime * tauinv_depsub ) )
+          !qvapinmix_lim = qsat(i) - qiceinmix / ( 1. - EXP( - dtime * tauinv_depsub ) )
+          qvapinmix_lim = qsat(i) - qiceinmix * MAX(1., 1.5 / ( dtime * tauinv_depsub ))
           qvapinclr_lim = qvapinmix_lim * ( 1. + cldfra_mix / clrfra_mix ) &
                         - qvc(i) / cldfra(i) * cldfra_mix / clrfra_mix

LMDZ6/branches/contrails/libf/phylmd/phyetat0_mod.f90

@@ -18 +18 @@
   USE surface_data,     ONLY : type_ocean, version_ocean
   USE phyetat0_get_mod, ONLY : phyetat0_get, phyetat0_srf
-  USE phys_state_var_mod, ONLY : ancien_ok, cwcon, clwcon, detr_therm, phys_tstep, &
+  USE phys_state_var_mod, ONLY : ancien_ok, clwcon, detr_therm, phys_tstep, &
        qsol, fevap, z0m, z0h, agesno, &
        du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
        falb_dir, falb_dif, prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien, &
        ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, qbs_ancien, &
-       cf_ancien, qvc_ancien, cfa_ancien, pcf_ancien, qva_ancien, qia_ancien, &
+       cf_ancien, qvc_ancien, qvcon, qccon, cfa_ancien, pcf_ancien, qva_ancien, qia_ancien, &
        radpas, radsol, rain_fall, &
        ratqs, rnebcon, rugoro, sig1, snow_fall, bs_fall, solaire_etat0, sollw, sollwdown, &
@@ -416 +416 @@
     ancien_ok=ancien_ok.AND.phyetat0_get(cf_ancien,"CFANCIEN","CFANCIEN",0.)
     ancien_ok=ancien_ok.AND.phyetat0_get(qvc_ancien,"QVCANCIEN","QVCANCIEN",0.)
-    found=phyetat0_get(cwcon,"CWCON","CWCON",0.)
+    found=phyetat0_get(qvcon,"QVCON","QVCON",0.)
+    found=phyetat0_get(qccon,"QCCON","QCCON",0.)
   ELSE
     cf_ancien(:,:)=0.

LMDZ6/branches/contrails/libf/phylmd/phyredem.f90

@@ -23 +23 @@
                                 prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien,      &
                                 ql_ancien, qs_ancien, qbs_ancien, cf_ancien, &
+                                qvcon, qccon,                                &
                                 qvc_ancien, cfa_ancien, pcf_ancien,          &
                                 qva_ancien, qia_ancien, u_ancien, v_ancien,  &
-                                cwcon, clwcon, rnebcon, ratqs, pbl_tke,      &
+                                clwcon, rnebcon, ratqs, pbl_tke,             &
                                 wake_delta_pbl_tke, zmax0, f0, sig1, w01,    &
                                 wake_deltat, wake_deltaq, wake_s, awake_s,   &
@@ -255 +256 @@
       CALL put_field(pass,"CFANCIEN", "CFANCIEN", cf_ancien)
       CALL put_field(pass,"QVCANCIEN", "QVCANCIEN", qvc_ancien)
-      CALL put_field(pass,"CWCON", "CWCON", cwcon)
+      CALL put_field(pass,"QVCON", "QVCON", qvcon)
+      CALL put_field(pass,"QCCON", "QCCON", qccon)
     ENDIF
 

LMDZ6/branches/contrails/libf/phylmd/phys_state_var_mod.F90

@@ -94 +94 @@
       REAL, ALLOCATABLE, SAVE :: cf_ancien(:,:), qvc_ancien(:,:)
 !$OMP THREADPRIVATE(cf_ancien, qvc_ancien)
+      REAL, ALLOCATABLE, SAVE :: qvcon(:,:), qccon(:,:)
+!$OMP THREADPRIVATE(qvcon, qccon)
       REAL, ALLOCATABLE, SAVE :: cfa_ancien(:,:), pcf_ancien(:,:)
 !$OMP THREADPRIVATE(cfa_ancien, pcf_ancien)
       REAL, ALLOCATABLE, SAVE :: qva_ancien(:,:), qia_ancien(:,:)
 !$OMP THREADPRIVATE(qva_ancien, qia_ancien)
-      REAL, ALLOCATABLE, SAVE :: cwcon(:,:), cwcon0(:,:)
-!$OMP THREADPRIVATE(cwcon, cwcon0)
 !!! RomP >>>
       REAL, ALLOCATABLE, SAVE :: tr_ancien(:,:,:)
@@ -596 +596 @@
       ALLOCATE(cfa_ancien(klon,klev), pcf_ancien(klon,klev))
       ALLOCATE(qva_ancien(klon,klev), qia_ancien(klon,klev))
-      ALLOCATE(cwcon(klon,klev), cwcon0(klon,klev))
+      ALLOCATE(qvcon(klon,klev), qccon(klon,klev))
 !!! Rom P >>>
       ALLOCATE(tr_ancien(klon,klev,nbtr))
@@ -826 +826 @@
       DEALLOCATE(cf_ancien, qvc_ancien, cfa_ancien, pcf_ancien)
       DEALLOCATE(qva_ancien, qia_ancien)
-      DEALLOCATE(cwcon, cwcon0)
+      DEALLOCATE(qvcon, qccon)
       DEALLOCATE(tr_ancien)                           !RomP
       DEALLOCATE(ratqs, pbl_tke,coefh,coefm)

LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90

@@ -1628 +1628 @@
        clwcon(:,:) = 0.0
        !--AB for prognostic clouds
-       cwcon0(:,:) = 0.0
-       cwcon(:,:) = 0.0
+       qvcon(:,:) = 0.0
+       qccon(:,:) = 0.0
 
        !
@@ -3950 +3950 @@
         DO i = 1, klon
           qvc_seri(i,k) = qvc_seri(i,k) * q_seri(i,k)
-          cwcon0(i,k) = zqsat(i,k)
         ENDDO
       ENDDO
@@ -3971 +3970 @@
     CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay,omega, &
          t_seri, q_seri, ql_seri_lscp, qi_seri_lscp, ratqs, sigma_qtherm, &
-         ptconv, rnebcon, cwcon, clwcon, rnebcon0, cwcon0, clwcon0, &
+         ptconv, rnebcon, qvcon, qccon, rnebcon0, zqsat, clwcon0, &
          d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, &
          pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb,  &
@@ -3991 +3990 @@
          dqrfreez, dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez, &
          dqised, dcfsed, dqvcsed)
-
-    IF (ok_ice_supersat) cwcon(:,:) = cwcon0(:,:)
 
     ELSE