Changeset 4950 for LMDZ6/trunk/libf/phylmd/StratAer/micphy_tstep.F90

Timestamp:

May 22, 2024, 3:16:36 PM (5 weeks ago)

Author:

lebasn

Message:

StratAer?: New model version (microphysic, composition routine, code cleaning, new params...)

File:

• LMDZ6/trunk/libf/phylmd/StratAer/micphy_tstep.F90 (modified) (9 diffs)

LMDZ6/trunk/libf/phylmd/StratAer/micphy_tstep.F90

@@ -8 +8 @@
   USE aerophys
   USE infotrac_phy, ONLY : nbtr_bin, nbtr_sulgas, nbtr, id_H2SO4_strat
-  USE phys_local_var_mod, ONLY: mdw, budg_3D_nucl, budg_3D_cond_evap, budg_h2so4_to_part, R2SO4, DENSO4, f_r_wet
+  USE phys_local_var_mod, ONLY: mdw, budg_3D_nucl, budg_3D_cond_evap, budg_h2so4_to_part, R2SO4, DENSO4, &
+       f_r_wet, R2SO4B, DENSO4B, f_r_wetB
   USE nucleation_tstep_mod
   USE cond_evap_tstep_mod
@@ -14 +15 @@
   USE YOMCST, ONLY : RPI, RD, RG
   USE print_control_mod, ONLY: lunout
-  USE strataer_local_var_mod
+  USE strataer_local_var_mod ! contains also RRSI and Vbin
   
   IMPLICIT NONE
@@ -35 +36 @@
   REAL                      :: ntot !total number of molecules in the critical cluster (ntot>4)
   REAL                      :: x    ! molefraction of H2SO4 in the critical cluster     
-  REAL Vbin(nbtr_bin)
   REAL a_xm, b_xm, c_xm
   REAL PDT, dt
   REAL H2SO4_init
   REAL ACTSO4(klon,klev)
-  REAL RRSI(nbtr_bin)
   REAL nucl_rate
   REAL cond_evap_rate
@@ -48 +47 @@
   REAL DNDR(nbtr_bin)
   REAL H2SO4_sat
-
-  DO it=1,nbtr_bin
-    Vbin(it)=4.0*RPI*((mdw(it)/2.)**3)/3.0
-  ENDDO
-
+  REAL R2SO4ik(nbtr_bin), DENSO4ik(nbtr_bin), f_r_wetik(nbtr_bin)
+  
   !coefficients for H2SO4 density parametrization used for nucleation if ntot<4
   a_xm = 0.7681724 + 1.*(2.1847140 + 1.*(7.1630022 + 1.*(-44.31447 + &
@@ -61 +57 @@
        & 1.*(7.990811e-4 + 1.*(-7.458060e-4 + 1.*2.58139e-4 )))))
 
-  ! STRAACT (R2SO4, t_seri -> H2SO4 activity coefficient (ACTSO4)) for cond/evap
-  CALL STRAACT(ACTSO4)
-
-  ! compute particle radius in cm RRSI from diameter in m 
-  DO it=1,nbtr_bin
-    RRSI(it)=mdw(it)/2.*100.
-  ENDDO
-
+  IF(.not.flag_new_strat_compo) THEN
+     ! STRAACT (R2SO4, t_seri -> H2SO4 activity coefficient (ACTSO4)) for cond/evap
+     CALL STRAACT(ACTSO4)
+  ENDIF
+  
   DO ilon=1, klon
 !
@@ -104 +97 @@
       ENDIF
       ! compute cond/evap rate in kg(H2SO4)/kgA/s
-      CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
-             & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
-             & RRSI,Vbin,FL,ASO4,DNDR)
+      IF(flag_new_strat_compo) THEN
+         R2SO4ik(:)   = R2SO4B(ilon,ilev,:)
+         DENSO4ik(:)  = DENSO4B(ilon,ilev,:)
+         f_r_wetik(:) = f_r_wetB(ilon,ilev,:)
+         CALL condens_evapor_rate_kelvin(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
+              & R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
+              & R2SO4ik,DENSO4ik,f_r_wetik,FL,ASO4,DNDR)
+      ELSE
+         CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
+              & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
+              & FL,ASO4,DNDR)
+      ENDIF
       ! Compute H2SO4 saturate vapor for big particules
       H2SO4_sat = DNDR(nbtr_bin)/(pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mH2SO4mol)
@@ -127 +129 @@
       tr_seri(ilon,ilev,id_H2SO4_strat)=MAX(0.,tr_seri(ilon,ilev,id_H2SO4_strat)-(nucl_rate+cond_evap_rate)*dt)
       ! apply cond to bins
-      CALL cond_evap_part(dt,FL,ASO4,f_r_wet(ilon,ilev),RRSI,Vbin,tr_seri(ilon,ilev,:))
+      CALL condens_evapor_part(dt,FL,ASO4,f_r_wet(ilon,ilev),tr_seri(ilon,ilev,:))
       ! apply nucl. to bins
-      CALL nucleation_part(nucl_rate,ntot,x,dt,Vbin,tr_seri(ilon,ilev,:))
+      CALL nucleation_part(nucl_rate,ntot,x,dt,tr_seri(ilon,ilev,:))
       ! compute fluxes as diagnostic in [kg(S)/m2/layer/s] (now - for evap and + for cond)
       budg_3D_cond_evap(ilon,ilev)=budg_3D_cond_evap(ilon,ilev)+mSatom/mH2SO4mol &
@@ -142 +144 @@
         & *pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mH2SO4mol
     ! compute cond/evap rate in kg(H2SO4)/kgA/s (now only evap for pdtphys)
-    CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
-           & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
-           & RRSI,Vbin,FL,ASO4,DNDR)
+    IF(flag_new_strat_compo) THEN
+       CALL condens_evapor_rate_kelvin(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
+            & R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
+            & R2SO4ik,DENSO4ik,f_r_wetik,FL,ASO4,DNDR)
+    ELSE
+       CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
+            & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
+            & FL,ASO4,DNDR)
+    ENDIF
     ! limit evaporation (negative FL) over one physics time step to H2SO4 content of the droplet
     DO it=1,nbtr_bin
@@ -159 +167 @@
     tr_seri(ilon,ilev,id_H2SO4_strat)=MAX(0.,tr_seri(ilon,ilev,id_H2SO4_strat)-evap_rate*pdtphys)
     ! apply evap to bins
-    CALL cond_evap_part(pdtphys,FL,ASO4,f_r_wet(ilon,ilev),RRSI,Vbin,tr_seri(ilon,ilev,:))
+    CALL condens_evapor_part(pdtphys,FL,ASO4,f_r_wet(ilon,ilev),tr_seri(ilon,ilev,:))
     ! compute fluxes as diagnostic in [kg(S)/m2/layer/s] (now - for evap and + for cond)
     budg_3D_cond_evap(ilon,ilev)=budg_3D_cond_evap(ilon,ilev)+mSatom/mH2SO4mol &