Changeset 5926

Timestamp:

Dec 11, 2025, 4:17:28 PM (5 weeks ago)

Author:

lebasn

Message:

StratAer?: bugfix on last commit + SAI mode

Location:

LMDZ6/trunk/libf/phylmd/StratAer

Files:

• aer_sedimnt.f90 (modified) (1 diff)
• coagulate.f90 (modified) (3 diffs)
• micphy_tstep.f90 (modified) (3 diffs)
• strataer_emiss_mod.f90 (modified) (5 diffs)
• traccoag_mod.f90 (modified) (3 diffs)

LMDZ6/trunk/libf/phylmd/StratAer/aer_sedimnt.f90

@@ -89 +89 @@
       ZRHO=pplay(JL,JK)/(RD*t_seri(JL,JK))
 
-      ! stokes-velocity with cunnigham slip- flow correction
-      IF(flag_strat_compo) THEN
+      ! Use Tabazadeh et al. 1997 for H2SO4/H2O composition
+      IF(flag_strat_compo.EQ.2) THEN
          ! stokes-velocity with cunnigham slip- flow correction
          ZVAER(JL,JK,nb) = 2./9.*(DENSO4B(JL,JK,nb)*1000.-ZRHO)*RG/zvis(JL,JK)*(f_r_wetB(JL,JK,nb)*mdw(nb)/2.)**2.* &

LMDZ6/trunk/libf/phylmd/StratAer/coagulate.f90

@@ -129 +129 @@
   IF (is_strato(ilon,ilev)) THEN
   !compute actual wet particle radius & volume for every grid box
-  IF(flag_strat_compo) THEN
+  IF(flag_strat_compo.EQ.2) THEN
      DO i=1, nbtr_bin
         radiuswet(i)=f_r_wetB(ilon,ilev,i)*mdw(i)/2.
@@ -174 +174 @@
 !--pre-compute the thermal velocity of a particle thvelpar(i) from equation 20
   thvelpar=0.0
-  IF(flag_strat_compo) THEN
+  IF(flag_strat_compo.EQ.2) THEN
      DO i=1, nbtr_bin
         m_par(i)=4./3.*RPI*radiuswet(i)**3.*DENSO4B(ilon,ilev,i)*1000.
@@ -275 +275 @@
 !                     =(Vdry(k)*tr_t(ilon,ilev,k)+pdtcoag*num_dry)/( (1.+pdtcoag*denom)*Vdry(k) ) 
 !          with num_dry=...beta(i,j)*Vdry(i)*....
-!       so in old STRATAER (.not.flag_strat_compo), it was correct
+!       so in old STRATAER (flag_strat_compo.EQ.1), it was correct
   ENDIF
 

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

@@ -57 +57 @@
        & 1.*(7.990811e-4 + 1.*(-7.458060e-4 + 1.*2.58139e-4 )))))
 
-  IF(.not.flag_strat_compo) THEN
+  IF(flag_strat_compo.EQ.1) THEN
      ! STRAACT (R2SO4, t_seri -> H2SO4 activity coefficient (ACTSO4)) for cond/evap
      CALL STRAACT(ACTSO4)
@@ -97 +97 @@
       ENDIF
       ! compute cond/evap rate in kg(H2SO4)/kgA/s
-      IF(flag_strat_compo) THEN
+      IF(flag_strat_compo.EQ.2) THEN
          R2SO4ik(:)   = R2SO4B(ilon,ilev,:)
          DENSO4ik(:)  = DENSO4B(ilon,ilev,:)
@@ -144 +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)
-    IF(flag_strat_compo) THEN
+    IF(flag_strat_compo.EQ.2) 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), &

LMDZ6/trunk/libf/phylmd/StratAer/strataer_emiss_mod.f90

@@ -53 +53 @@
                 ALLOCATE(m_H2O_emiss_vol(nErupt))
                 ALLOCATE(m_H2O_emiss_vol_daily(nErupt))
-                !                ALLOCATE(d_q_emiss(klon,klev))
                 ALLOCATE(budg_emi(klon,nAerErupt+1))
                 m_H2O_emiss_vol(:)=0.
                 m_H2O_emiss_vol_daily(:)=0.
-                !                d_q_emiss(:,:)=0.
              ELSE
                 ALLOCATE(budg_emi(klon,nAerErupt))
@@ -104 +102 @@
 
     ELSEIF (flag_emit == 2) THEN ! SAI
+       ALLOCATE(budg_emi(klon,1))
+       budg_emi(:,:)=0.
+       
        CALL getin_p('m_aer_emiss_sai',m_aer_emiss_sai)
        CALL getin_p('altemiss_sai',altemiss_sai)
@@ -117 +118 @@
 
     ELSEIF (flag_emit == 3) THEN ! SAI between latitudes
+       ALLOCATE(budg_emi(klon,1))
+       budg_emi(:,:)=0.
+       
        CALL getin_p('m_aer_emiss_sai',m_aer_emiss_sai)
        CALL getin_p('altemiss_sai',altemiss_sai)
@@ -124 +128 @@
        CALL getin_p('xlat_min_sai',xlat_min_sai)
     ENDIF
-
+    
     IF (flag_emit == 1) THEN
        DO ispec=1,nAerErupt
@@ -233 +237 @@
     dlat=180./2./FLOAT(nbp_lat)   ! d latitude in degree
     dlon=360./2./FLOAT(nbp_lon)   ! d longitude in degree
-    WRITE(lunout,*) 'IN STRATAER_INIT dlat=',dlat,'dlon=',dlon
-    WRITE(lunout,*) 'IN STRATAER_INIT nErupt=',nErupt
-    WRITE(lunout,*) 'IN STRATAER_INIT xlat_min=',xlat_min_vol,'xlat_max=',xlat_max_vol
-    WRITE(lunout,*) 'IN STRATAER_INIT xlon_min=',xlon_min_vol,'xlon_max=',xlon_max_vol
+    IF (flag_emit.EQ.1 .OR. flag_emit.EQ.4) THEN ! Volcano
+       WRITE(lunout,*) 'IN STRATAER_INIT dlat=',dlat,'dlon=',dlon
+       WRITE(lunout,*) 'IN STRATAER_INIT nErupt=',nErupt
+       WRITE(lunout,*) 'IN STRATAER_INIT xlat_min=',xlat_min_vol,'xlat_max=',xlat_max_vol
+       WRITE(lunout,*) 'IN STRATAER_INIT xlon_min=',xlon_min_vol,'xlon_max=',xlon_max_vol
+    ENDIF
 
     DO ieru=1, nErupt

LMDZ6/trunk/libf/phylmd/StratAer/traccoag_mod.f90

@@ -121 +121 @@
     WHERE (stratomask.GT.0.5) is_strato=.TRUE. 
 
-    IF(flag_strat_compo) THEN
+    IF(flag_strat_compo.EQ.2) THEN
        IF(debutphy) WRITE(lunout,*) 'traccoag: COMPO/DENSITY (Tabazadeh 97) + H2O kelvin effect', flag_strat_compo
        ! STRACOMP (H2O, P, t_seri, R -> R2SO4 + Kelvin effect) : Taba97, Socol, etc...
@@ -146 +146 @@
     
 !--initialise emission diagnostics
-    if (nErupt > 0 .and. (flag_emit == 1 .or. flag_emit == 4)) budg_emi(:,1)=0.0
+    budg_emi(:,1)=0.0
     budg_emi_ocs(:)=0.0
     budg_emi_so2(:)=0.0
@@ -383 +383 @@
                    & *pplay(i,j)/t_seri(i,j)/RD                           ! [air mass concentration in kg air /m3A]
               
-              IF(flag_strat_compo) THEN
+              IF(flag_strat_compo.EQ.2) THEN
                  !     SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3)
                  SAD_sulfate(i,j) = SAD_sulfate(i,j) + nd_mode(i,j,it) &     ! [DRY part/m3A (in bin it)]