Changeset 3356 for LMDZ6/branches/DYNAMICO-conv/libf/phylmd/StratAer

Timestamp:

Jun 29, 2018, 12:31:11 PM (6 years ago)

Author:

Laurent Fairhead

Message:

First attempt at merging with trunk

Location:

LMDZ6/branches/DYNAMICO-conv

Files:

• . (modified) (1 prop)
• libf/phylmd/StratAer/aerophys.F90 (modified) (1 diff)
• libf/phylmd/StratAer/interp_sulf_input.F90 (modified) (2 diffs)
• libf/phylmd/StratAer/micphy_tstep.F90 (modified) (8 diffs)
• libf/phylmd/StratAer/traccoag_mod.F90 (modified) (1 diff)

LMDZ6/branches/DYNAMICO-conv/libf/phylmd/StratAer/aerophys.F90

@@ -2 +2 @@
 !
 MODULE aerophys
+!
+  IMPLICIT NONE
 !
   REAL,PARAMETER                         :: ropx=1500.0              ! default aerosol particle mass density [kg/m3]

LMDZ6/branches/DYNAMICO-conv/libf/phylmd/StratAer/interp_sulf_input.F90

@@ -32 +32 @@
   INTEGER, SAVE :: n_lev   ! number of levels in the input data
   INTEGER n_mth   ! number of months in the input data
-  REAL OCS_tmp, SO2_tmp
   INTEGER, SAVE :: mth_pre
 !$OMP THREADPRIVATE(mth_pre)
@@ -238 +237 @@
     DO k=1, klev
       !
-      OCS_tmp=tr_seri(i,k,id_OCS_strat)
-      SO2_tmp=tr_seri(i,k,id_SO2_strat)
-      !--OCS and SO2 prescribed below p_bound
+      !--OCS and SO2 prescribed back to their clim values below p_bound
       IF (paprs(i,k).GT.p_bound) THEN
+        budg_3D_backgr_ocs(i,k)=OCS_clim(i,k)-tr_seri(i,k,id_OCS_strat)
+        budg_3D_backgr_so2(i,k)=SO2_clim(i,k)-tr_seri(i,k,id_SO2_strat)
         tr_seri(i,k,id_OCS_strat)=OCS_clim(i,k)
         tr_seri(i,k,id_SO2_strat)=SO2_clim(i,k)
       ENDIF
-      budg_3D_backgr_ocs(i,k)=tr_seri(i,k,id_OCS_strat)-OCS_tmp
-      budg_3D_backgr_so2(i,k)=tr_seri(i,k,id_SO2_strat)-SO2_tmp
     ENDDO
   ENDDO

LMDZ6/branches/DYNAMICO-conv/libf/phylmd/StratAer/micphy_tstep.F90

@@ -4 +4 @@
   USE aerophys
   USE infotrac
-  USE phys_local_var_mod, ONLY: mdw, budg_3D_nucl, budg_3D_cond_evap, 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
   USE nucleation_tstep_mod
   USE cond_evap_tstep_mod
@@ -25 +25 @@
   ! local variables in coagulation routine
   INTEGER, PARAMETER        :: nbtstep=4  ! Max number of time steps in microphysics per time step in physics
-  INTEGER                   :: it,ilon,ilev,IK,count_tstep
+  INTEGER                   :: it,ilon,ilev,count_tstep
   REAL                      :: rhoa !H2SO4 number density [molecules/cm3]
   REAL                      :: ntot !total number of molecules in the critical cluster (ntot>4)
@@ -43 +43 @@
   REAL H2SO4_sat(nbtr_bin)
 
-  DO IK=1,nbtr_bin
-    Vbin(IK)=4.0*RPI*((mdw(IK)/2.)**3)/3.0
+  DO it=1,nbtr_bin
+    Vbin(it)=4.0*RPI*((mdw(it)/2.)**3)/3.0
   ENDDO
 
@@ -64 +64 @@
 
   DO ilon=1, klon
+!
+!--initialisation of diagnostic
+  budg_h2so4_to_part(ilon)=0.0
+!
   DO ilev=1, klev
+!
+!--initialisation of diagnostic
+  budg_3D_nucl(ilon,ilev)=0.0
+  budg_3D_cond_evap(ilon,ilev)=0.0
+!
   ! only in the stratosphere
   IF (is_strato(ilon,ilev)) THEN
     ! initialize sulfur fluxes
-    budg_3D_nucl(ilon,ilev)=0.0
-    budg_3D_cond_evap(ilon,ilev)=0.0
     H2SO4_init=tr_seri(ilon,ilev,id_H2SO4_strat)
     ! adaptive timestep for nucleation and condensation
@@ -88 +95 @@
              & RRSI,Vbin,FL,ASO4,DNDR)
       ! consider only condensation (positive FL)
-      DO IK=1,nbtr_bin
-        FL(IK)=MAX(FL(IK),0.)
+      DO it=1,nbtr_bin
+        FL(it)=MAX(FL(it),0.)
       ENDDO
       ! compute total H2SO4 cond flux for all particles
       cond_evap_rate=0.0
-      DO IK=1, nbtr_bin
-        cond_evap_rate=cond_evap_rate+tr_seri(ilon,ilev,IK+nbtr_sulgas)*FL(IK)*mH2SO4mol
+      DO it=1, nbtr_bin
+        cond_evap_rate=cond_evap_rate+tr_seri(ilon,ilev,it+nbtr_sulgas)*FL(it)*mH2SO4mol
       ENDDO
       ! determine appropriate time step
@@ -124 +131 @@
            & RRSI,Vbin,FL,ASO4,DNDR)
     ! limit evaporation (negative FL) over one physics time step to H2SO4 content of the droplet
-    DO IK=1,nbtr_bin
-      FL(IK)=MAX(FL(IK)*pdtphys,0.-ASO4(IK))/pdtphys
+    DO it=1,nbtr_bin
+      FL(it)=MAX(FL(it)*pdtphys,0.-ASO4(it))/pdtphys
       ! consider only evap (negative FL)
-      FL(IK)=MIN(FL(IK),0.)
+      FL(it)=MIN(FL(it),0.)
     ENDDO
     ! compute total H2SO4 evap flux for all particles
     evap_rate=0.0
-    DO IK=1, nbtr_bin
-      evap_rate=evap_rate+tr_seri(ilon,ilev,IK+nbtr_sulgas)*FL(IK)*mH2SO4mol
+    DO it=1, nbtr_bin
+      evap_rate=evap_rate+tr_seri(ilon,ilev,it+nbtr_sulgas)*FL(it)*mH2SO4mol
     ENDDO
     ! update H2SO4 concentration after evap
@@ -141 +148 @@
     budg_3D_cond_evap(ilon,ilev)=budg_3D_cond_evap(ilon,ilev)+mSatom/mH2SO4mol &
              & *evap_rate*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG
+    ! compute vertically integrated flux due to the net effect of nucleation and condensation/evaporation
+    budg_h2so4_to_part(ilon)=budg_h2so4_to_part(ilon)+(H2SO4_init-tr_seri(ilon,ilev,id_H2SO4_strat)) &
+             & *mSatom/mH2SO4mol*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys
   ENDIF
   ENDDO
@@ -148 +158 @@
     DO ilon=1, klon
     DO ilev=1, klev    
-    DO IK=1, nbtr
-      IF (tr_seri(ilon,ilev,IK).LT.0.0) THEN
-        PRINT *, 'micphy_tstep: negative concentration', tr_seri(ilon,ilev,IK), ilon, ilev, IK
+    DO it=1, nbtr
+      IF (tr_seri(ilon,ilev,it).LT.0.0) THEN
+        PRINT *, 'micphy_tstep: negative concentration', tr_seri(ilon,ilev,it), ilon, ilev, it
       ENDIF
     ENDDO

LMDZ6/branches/DYNAMICO-conv/libf/phylmd/StratAer/traccoag_mod.F90

@@ -203 +203 @@
             !vertically distributed emission
             DO k=1, klev
-              ! stretch emission over one day of Pinatubo eruption
-              emission=m_aer_emiss_vol*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/1./86400. 
+              ! stretch emission over one day (minus one timestep) of Pinatubo eruption
+              emission=m_aer_emiss_vol*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/1./(86400.-pdtphys)
               tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys 
               budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol