[3526] | 1 | ! |
---|
| 2 | ! $Id: micphy_tstep.F90 4601 2023-06-30 22:07:30Z nfevrier $ |
---|
| 3 | ! |
---|
[2690] | 4 | SUBROUTINE micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato) |
---|
| 5 | |
---|
[3526] | 6 | USE geometry_mod, ONLY : latitude_deg !NL- latitude corr. to local domain |
---|
[2690] | 7 | USE dimphy, ONLY : klon,klev |
---|
| 8 | USE aerophys |
---|
[4293] | 9 | USE infotrac_phy, ONLY : nbtr_bin, nbtr_sulgas, nbtr, id_H2SO4_strat |
---|
[3094] | 10 | USE phys_local_var_mod, ONLY: mdw, budg_3D_nucl, budg_3D_cond_evap, budg_h2so4_to_part, R2SO4, DENSO4, f_r_wet |
---|
[2690] | 11 | USE nucleation_tstep_mod |
---|
| 12 | USE cond_evap_tstep_mod |
---|
| 13 | USE sulfate_aer_mod, ONLY : STRAACT |
---|
[2695] | 14 | USE YOMCST, ONLY : RPI, RD, RG |
---|
[3526] | 15 | USE print_control_mod, ONLY: lunout |
---|
[4601] | 16 | USE strataer_local_var_mod |
---|
[3526] | 17 | |
---|
[2690] | 18 | IMPLICIT NONE |
---|
| 19 | |
---|
| 20 | !-------------------------------------------------------- |
---|
| 21 | |
---|
| 22 | ! transfer variables when calling this routine |
---|
| 23 | REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde) |
---|
| 24 | REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] |
---|
| 25 | REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature |
---|
| 26 | REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) |
---|
| 27 | REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa) |
---|
| 28 | REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! humidite relative |
---|
| 29 | LOGICAL,DIMENSION(klon,klev),INTENT(IN) :: is_strato |
---|
| 30 | |
---|
| 31 | ! local variables in coagulation routine |
---|
| 32 | INTEGER, PARAMETER :: nbtstep=4 ! Max number of time steps in microphysics per time step in physics |
---|
[3098] | 33 | INTEGER :: it,ilon,ilev,count_tstep |
---|
[2690] | 34 | REAL :: rhoa !H2SO4 number density [molecules/cm3] |
---|
| 35 | REAL :: ntot !total number of molecules in the critical cluster (ntot>4) |
---|
| 36 | REAL :: x ! molefraction of H2SO4 in the critical cluster |
---|
| 37 | REAL Vbin(nbtr_bin) |
---|
| 38 | REAL a_xm, b_xm, c_xm |
---|
| 39 | REAL PDT, dt |
---|
| 40 | REAL H2SO4_init |
---|
| 41 | REAL ACTSO4(klon,klev) |
---|
| 42 | REAL RRSI(nbtr_bin) |
---|
| 43 | REAL nucl_rate |
---|
| 44 | REAL cond_evap_rate |
---|
| 45 | REAL evap_rate |
---|
| 46 | REAL FL(nbtr_bin) |
---|
| 47 | REAL ASO4(nbtr_bin) |
---|
| 48 | REAL DNDR(nbtr_bin) |
---|
[4293] | 49 | REAL H2SO4_sat |
---|
[2690] | 50 | |
---|
[3098] | 51 | DO it=1,nbtr_bin |
---|
| 52 | Vbin(it)=4.0*RPI*((mdw(it)/2.)**3)/3.0 |
---|
[2690] | 53 | ENDDO |
---|
| 54 | |
---|
| 55 | !coefficients for H2SO4 density parametrization used for nucleation if ntot<4 |
---|
| 56 | a_xm = 0.7681724 + 1.*(2.1847140 + 1.*(7.1630022 + 1.*(-44.31447 + & |
---|
| 57 | & 1.*(88.75606 + 1.*(-75.73729 + 1.*23.43228))))) |
---|
| 58 | b_xm = 1.808225e-3 + 1.*(-9.294656e-3 + 1.*(-0.03742148 + 1.*(0.2565321 + & |
---|
| 59 | & 1.*(-0.5362872 + 1.*(0.4857736 - 1.*0.1629592))))) |
---|
| 60 | c_xm = -3.478524e-6 + 1.*(1.335867e-5 + 1.*(5.195706e-5 + 1.*(-3.717636e-4 + & |
---|
| 61 | & 1.*(7.990811e-4 + 1.*(-7.458060e-4 + 1.*2.58139e-4 ))))) |
---|
| 62 | |
---|
| 63 | ! STRAACT (R2SO4, t_seri -> H2SO4 activity coefficient (ACTSO4)) for cond/evap |
---|
| 64 | CALL STRAACT(ACTSO4) |
---|
| 65 | |
---|
| 66 | ! compute particle radius in cm RRSI from diameter in m |
---|
| 67 | DO it=1,nbtr_bin |
---|
| 68 | RRSI(it)=mdw(it)/2.*100. |
---|
| 69 | ENDDO |
---|
| 70 | |
---|
| 71 | DO ilon=1, klon |
---|
[3094] | 72 | ! |
---|
| 73 | !--initialisation of diagnostic |
---|
| 74 | budg_h2so4_to_part(ilon)=0.0 |
---|
| 75 | ! |
---|
[2690] | 76 | DO ilev=1, klev |
---|
[3094] | 77 | ! |
---|
| 78 | !--initialisation of diagnostic |
---|
| 79 | budg_3D_nucl(ilon,ilev)=0.0 |
---|
| 80 | budg_3D_cond_evap(ilon,ilev)=0.0 |
---|
| 81 | ! |
---|
[2690] | 82 | ! only in the stratosphere |
---|
| 83 | IF (is_strato(ilon,ilev)) THEN |
---|
| 84 | ! initialize sulfur fluxes |
---|
| 85 | H2SO4_init=tr_seri(ilon,ilev,id_H2SO4_strat) |
---|
| 86 | ! adaptive timestep for nucleation and condensation |
---|
| 87 | PDT=pdtphys |
---|
| 88 | count_tstep=0 |
---|
[2695] | 89 | DO WHILE (PDT>0.0) |
---|
[2690] | 90 | count_tstep=count_tstep+1 |
---|
[2695] | 91 | IF (count_tstep .GT. nbtstep) EXIT |
---|
[2690] | 92 | ! convert tr_seri(GASH2SO4) (in kg/kgA) to H2SO4 number density (in molecules/cm3) |
---|
| 93 | rhoa=tr_seri(ilon,ilev,id_H2SO4_strat) & |
---|
| 94 | & *pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mH2SO4mol |
---|
| 95 | ! compute nucleation rate in kg(H2SO4)/kgA/s |
---|
| 96 | CALL nucleation_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev),rh(ilon,ilev), & |
---|
[3526] | 97 | & a_xm,b_xm,c_xm,nucl_rate,ntot,x) |
---|
| 98 | !NL - add nucleation box (if flag on) |
---|
| 99 | IF (flag_nuc_rate_box) THEN |
---|
[3527] | 100 | IF (latitude_deg(ilon).LE.nuclat_min .OR. latitude_deg(ilon).GE.nuclat_max & |
---|
| 101 | .OR. pplay(ilon,ilev).GE.nucpres_max .AND. pplay(ilon,ilev).LE.nucpres_min) THEN |
---|
[3526] | 102 | nucl_rate=0.0 |
---|
| 103 | ENDIF |
---|
| 104 | ENDIF |
---|
[2690] | 105 | ! compute cond/evap rate in kg(H2SO4)/kgA/s |
---|
| 106 | CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), & |
---|
| 107 | & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), & |
---|
| 108 | & RRSI,Vbin,FL,ASO4,DNDR) |
---|
[4293] | 109 | ! Compute H2SO4 saturate vapor for big particules |
---|
| 110 | H2SO4_sat = DNDR(nbtr_bin)/(pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mH2SO4mol) |
---|
[2690] | 111 | ! consider only condensation (positive FL) |
---|
[3098] | 112 | DO it=1,nbtr_bin |
---|
| 113 | FL(it)=MAX(FL(it),0.) |
---|
[2690] | 114 | ENDDO |
---|
| 115 | ! compute total H2SO4 cond flux for all particles |
---|
| 116 | cond_evap_rate=0.0 |
---|
[3098] | 117 | DO it=1, nbtr_bin |
---|
| 118 | cond_evap_rate=cond_evap_rate+tr_seri(ilon,ilev,it+nbtr_sulgas)*FL(it)*mH2SO4mol |
---|
[2690] | 119 | ENDDO |
---|
| 120 | ! determine appropriate time step |
---|
[4293] | 121 | dt=(H2SO4_init-H2SO4_sat)/float(nbtstep)/MAX(1.e-30, nucl_rate+cond_evap_rate) !cond_evap_rate pos. for cond. and neg. for evap. |
---|
[2690] | 122 | IF (dt.LT.0.0) THEN |
---|
| 123 | dt=PDT |
---|
| 124 | ENDIF |
---|
| 125 | dt=MIN(dt,PDT) |
---|
| 126 | ! update H2SO4 concentration |
---|
| 127 | tr_seri(ilon,ilev,id_H2SO4_strat)=MAX(0.,tr_seri(ilon,ilev,id_H2SO4_strat)-(nucl_rate+cond_evap_rate)*dt) |
---|
| 128 | ! apply cond to bins |
---|
| 129 | CALL cond_evap_part(dt,FL,ASO4,f_r_wet(ilon,ilev),RRSI,Vbin,tr_seri(ilon,ilev,:)) |
---|
| 130 | ! apply nucl. to bins |
---|
| 131 | CALL nucleation_part(nucl_rate,ntot,x,dt,Vbin,tr_seri(ilon,ilev,:)) |
---|
| 132 | ! compute fluxes as diagnostic in [kg(S)/m2/layer/s] (now - for evap and + for cond) |
---|
[2752] | 133 | budg_3D_cond_evap(ilon,ilev)=budg_3D_cond_evap(ilon,ilev)+mSatom/mH2SO4mol & |
---|
[2690] | 134 | & *cond_evap_rate*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*dt/pdtphys |
---|
[2752] | 135 | budg_3D_nucl(ilon,ilev)=budg_3D_nucl(ilon,ilev)+mSatom/mH2SO4mol & |
---|
[2690] | 136 | & *nucl_rate*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*dt/pdtphys |
---|
| 137 | ! update time step |
---|
| 138 | PDT=PDT-dt |
---|
| 139 | ENDDO |
---|
| 140 | ! convert tr_seri(GASH2SO4) (in kg/kgA) to H2SO4 number density (in molecules/cm3) |
---|
| 141 | rhoa=tr_seri(ilon,ilev,id_H2SO4_strat) & |
---|
| 142 | & *pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mH2SO4mol |
---|
| 143 | ! compute cond/evap rate in kg(H2SO4)/kgA/s (now only evap for pdtphys) |
---|
| 144 | CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), & |
---|
| 145 | & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), & |
---|
| 146 | & RRSI,Vbin,FL,ASO4,DNDR) |
---|
| 147 | ! limit evaporation (negative FL) over one physics time step to H2SO4 content of the droplet |
---|
[3098] | 148 | DO it=1,nbtr_bin |
---|
| 149 | FL(it)=MAX(FL(it)*pdtphys,0.-ASO4(it))/pdtphys |
---|
[2690] | 150 | ! consider only evap (negative FL) |
---|
[3098] | 151 | FL(it)=MIN(FL(it),0.) |
---|
[2690] | 152 | ENDDO |
---|
| 153 | ! compute total H2SO4 evap flux for all particles |
---|
| 154 | evap_rate=0.0 |
---|
[3098] | 155 | DO it=1, nbtr_bin |
---|
| 156 | evap_rate=evap_rate+tr_seri(ilon,ilev,it+nbtr_sulgas)*FL(it)*mH2SO4mol |
---|
[2690] | 157 | ENDDO |
---|
| 158 | ! update H2SO4 concentration after evap |
---|
| 159 | tr_seri(ilon,ilev,id_H2SO4_strat)=MAX(0.,tr_seri(ilon,ilev,id_H2SO4_strat)-evap_rate*pdtphys) |
---|
| 160 | ! apply evap to bins |
---|
| 161 | CALL cond_evap_part(pdtphys,FL,ASO4,f_r_wet(ilon,ilev),RRSI,Vbin,tr_seri(ilon,ilev,:)) |
---|
| 162 | ! compute fluxes as diagnostic in [kg(S)/m2/layer/s] (now - for evap and + for cond) |
---|
[2752] | 163 | budg_3D_cond_evap(ilon,ilev)=budg_3D_cond_evap(ilon,ilev)+mSatom/mH2SO4mol & |
---|
[2690] | 164 | & *evap_rate*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG |
---|
[3094] | 165 | ! compute vertically integrated flux due to the net effect of nucleation and condensation/evaporation |
---|
| 166 | budg_h2so4_to_part(ilon)=budg_h2so4_to_part(ilon)+(H2SO4_init-tr_seri(ilon,ilev,id_H2SO4_strat)) & |
---|
| 167 | & *mSatom/mH2SO4mol*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys |
---|
[2690] | 168 | ENDIF |
---|
| 169 | ENDDO |
---|
| 170 | ENDDO |
---|
| 171 | |
---|
[2695] | 172 | IF (MINVAL(tr_seri).LT.0.0) THEN |
---|
[2690] | 173 | DO ilon=1, klon |
---|
| 174 | DO ilev=1, klev |
---|
[3098] | 175 | DO it=1, nbtr |
---|
| 176 | IF (tr_seri(ilon,ilev,it).LT.0.0) THEN |
---|
[3526] | 177 | WRITE(lunout,*) 'micphy_tstep: negative concentration', tr_seri(ilon,ilev,it), ilon, ilev, it |
---|
[2690] | 178 | ENDIF |
---|
| 179 | ENDDO |
---|
| 180 | ENDDO |
---|
| 181 | ENDDO |
---|
| 182 | ENDIF |
---|
| 183 | |
---|
| 184 | END SUBROUTINE micphy_tstep |
---|