Changeset 4482 for LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ice_sursat_mod.F90

Timestamp:

Mar 29, 2023, 3:14:27 PM (19 months ago)

Author:

lguez

Message:

Sync latest trunk changes to branch LMDZ_ECRad

Location:

LMDZ6/branches/LMDZ_ECRad

Files:

• . (modified) (1 prop)
• libf/phylmd/ice_sursat_mod.F90 (modified) (5 diffs)

LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ice_sursat_mod.F90

@@ -3 +3 @@
 IMPLICIT NONE
 
+!--flight inventories
+!
 REAL, SAVE, ALLOCATABLE :: flight_m(:,:)    !--flown distance m s-1 per cell
 !$OMP THREADPRIVATE(flight_m)
 REAL, SAVE, ALLOCATABLE :: flight_h2o(:,:)  !--emitted kg H2O s-1 per cell
 !$OMP THREADPRIVATE(flight_h2o)
-
-!--parameters for gamma function
-!--Karcher and Lohmann (2002)
-!--gamma = 2.583 - t / 207.83
-! REAL, SAVE, PARAMETER :: gamma0=2.583, Tgamma=207.83
-!--Ren and MacKenzie (2005) reused by Kärcher
-!--gamma = 2.349 - t / 259.0
-REAL, PARAMETER :: gamma0=2.349, Tgamma=259.0
-!
-!--number of clouds in cell (needs to be parametrized at some point)
-REAL, PARAMETER :: N_cld = 100.
+!
+!--Fixed Parameters
 !
 !--safety parameters for ERF function
@@ -43 +36 @@
 !  alpha_cld ). Dans le rapport il est appelé beta. Il varie entre 0 et 5
 !  (tun_ratqs = 0 => pas de modification de ratqs).
-
+!
+!--gamma0 and Tgamma: define RHcrit limit above which heterogeneous freezing occurs as a function of T
+!--Karcher and Lohmann (2002)
+!--gamma = 2.583 - t / 207.83
+!--Ren and MacKenzie (2005) reused by Kärcher
+!--gamma = 2.349 - t / 259.0
+!
+!--N_cld: number of clouds in cell (needs to be parametrized at some point)
+!
+!--contrail cross section: typical value found in Freudenthaler et al, GRL, 22, 3501-3504, 1995
+!--in m2, 1000x200 = 200 000 m2 after 15 min
+!
 REAL, SAVE :: chi=1.1, l_turb=50.0, tun_N=1.3, tun_ratqs=3.0
-!$OMP THREADPRIVATE(chi,l_turb,tun_N,tun_ratqs)
-!
-!--contrail cross section, typical value found in Freudenthaler et al, GRL, 22, 3501-3504, 1995
-!--in m2, 1000x200 = 200 000 m2 after 15 min
-REAL, PARAMETER :: contrail_cross_section=200000.0  
+REAL, SAVE :: gamma0=2.349, Tgamma=259.0, N_cld=100, contrail_cross_section=200000.0
+!$OMP THREADPRIVATE(chi,l_turb,tun_N,tun_ratqs,gamma0,Tgamma,N_cld,contrail_cross_section)
 
 CONTAINS
@@ -66 +67 @@
   CALL getin_p('flag_tun_N',tun_N)
   CALL getin_p('flag_tun_ratqs',tun_ratqs)
-
+  CALL getin_p('gamma0',gamma0)
+  CALL getin_p('Tgamma',Tgamma)
+  CALL getin_p('N_cld',N_cld)
+  CALL getin_p('contrail_cross_section',contrail_cross_section)
+
+  WRITE(lunout,*) 'Parameters for ice_sursat param'
   WRITE(lunout,*) 'flag_chi = ', chi
   WRITE(lunout,*) 'flag_l_turb = ', l_turb
   WRITE(lunout,*) 'flag_tun_N = ', tun_N
   WRITE(lunout,*) 'flag_tun_ratqs = ', tun_ratqs
+  WRITE(lunout,*) 'gamma0 = ', gamma0
+  WRITE(lunout,*) 'Tgamma = ', Tgamma
+  WRITE(lunout,*) 'N_cld = ', N_cld
+  WRITE(lunout,*) 'contrail_cross_section = ', contrail_cross_section
 
 END SUBROUTINE ice_sursat_init
@@ -384 +394 @@
      !
      !--Estimating gamma 
-     gamma_ss = gamma0 - t_actuel/Tgamma
-     !gamma_prec = gamma0 - t_ancien(i,k)/Tgamma      !--formulation initiale d Audran
-     gamma_prec = gamma0 - t/Tgamma                   !--autre formulation possible basée sur le T du pas de temps
+     gamma_ss = MAX(1.0, gamma0 - t_actuel/Tgamma)
+     !gamma_prec = MAX(1.0, gamma0 - t_ancien(i,k)/Tgamma)      !--formulation initiale d Audran
+     gamma_prec = MAX(1.0, gamma0 - t/Tgamma)                   !--autre formulation possible basée sur le T du pas de temps
      !
      ! Initialisation de qvc : q_sat du pas de temps precedent
@@ -421 +431 @@
         ! On approxime en serie entiere erf-1(x)
         qvc = 2.*rneb-1.
-        qvc = qvc + PI/12.*qvc**3 + 7.*PI**2/480.*qvc**5 + 127.*PI**3/40320.*qvc**7 + 4369.*PI**4/5806080.*qvc**9 + 34807.*PI**5/182476800.*qvc**11
+        qvc = qvc + PI/12.*qvc**3 + 7.*PI**2/480.*qvc**5 &
+             + 127.*PI**3/40320.*qvc**7 + 4369.*PI**4/5806080.*qvc**9 &
+             + 34807.*PI**5/182476800.*qvc**11
         qvc = sqrt(PI)/2.*qvc
         qvc = (qvc-pdf_b)*pdf_k*sqrt(2.)