Index: LMDZ6/trunk/libf/phylmd/StratAer/aer_sedimnt.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/StratAer/aer_sedimnt.F90 (revision 3524) +++ LMDZ6/trunk/libf/phylmd/StratAer/aer_sedimnt.F90 (revision 3526) @@ -1,2 +1,5 @@ +! +! $Id$ +! SUBROUTINE AER_SEDIMNT(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer) Index: LMDZ6/trunk/libf/phylmd/StratAer/calcaerosolstrato_rrtm.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/StratAer/calcaerosolstrato_rrtm.F90 (revision 3524) +++ LMDZ6/trunk/libf/phylmd/StratAer/calcaerosolstrato_rrtm.F90 (revision 3526) @@ -1,2 +1,5 @@ +! +! $Id$ +! SUBROUTINE calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut) Index: LMDZ6/trunk/libf/phylmd/StratAer/coagulate.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/StratAer/coagulate.F90 (revision 3524) +++ LMDZ6/trunk/libf/phylmd/StratAer/coagulate.F90 (revision 3526) @@ -1,2 +1,5 @@ +! +! $Id$ +! SUBROUTINE COAGULATE(pdtcoag,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato) ! ----------------------------------------------------------------------- Index: LMDZ6/trunk/libf/phylmd/StratAer/cond_evap_tstep_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/StratAer/cond_evap_tstep_mod.F90 (revision 3524) +++ LMDZ6/trunk/libf/phylmd/StratAer/cond_evap_tstep_mod.F90 (revision 3526) @@ -1,2 +1,5 @@ +! +! $Id$ +! MODULE cond_evap_tstep_mod Index: LMDZ6/trunk/libf/phylmd/StratAer/interp_sulf_input.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/StratAer/interp_sulf_input.F90 (revision 3524) +++ LMDZ6/trunk/libf/phylmd/StratAer/interp_sulf_input.F90 (revision 3526) @@ -1,2 +1,5 @@ +! +! $Id$ +! SUBROUTINE interp_sulf_input(debutphy,pdtphys,paprs,tr_seri) Index: LMDZ6/trunk/libf/phylmd/StratAer/micphy_tstep.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/StratAer/micphy_tstep.F90 (revision 3524) +++ LMDZ6/trunk/libf/phylmd/StratAer/micphy_tstep.F90 (revision 3526) @@ -1,4 +1,8 @@ +! +! $Id$ +! SUBROUTINE micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato) + USE geometry_mod, ONLY : latitude_deg !NL- latitude corr. to local domain USE dimphy, ONLY : klon,klev USE aerophys @@ -9,5 +13,7 @@ USE sulfate_aer_mod, ONLY : STRAACT USE YOMCST, ONLY : RPI, RD, RG - + USE print_control_mod, ONLY: lunout + USE strataer_mod + IMPLICIT NONE @@ -89,5 +95,12 @@ ! compute nucleation rate in kg(H2SO4)/kgA/s CALL nucleation_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev),rh(ilon,ilev), & - & a_xm,b_xm,c_xm,nucl_rate,ntot,x) + & a_xm,b_xm,c_xm,nucl_rate,ntot,x) + !NL - add nucleation box (if flag on) + IF (flag_nuc_rate_box) THEN + IF (latitude_deg(ilon).LE.(nuclat_min) .OR. latitude_deg(ilon).GE.(nuclat_max) & + .OR. pplay(ilon,ilev).GE.nucpres_max .AND. pplay(ilon,ilev) .LE. nucpres_min ) THEN + nucl_rate=0.0 + ENDIF + ENDIF ! compute cond/evap rate in kg(H2SO4)/kgA/s CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), & @@ -160,5 +173,5 @@ 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 + WRITE(lunout,*) 'micphy_tstep: negative concentration', tr_seri(ilon,ilev,it), ilon, ilev, it ENDIF ENDDO Index: LMDZ6/trunk/libf/phylmd/StratAer/miecalc_aer.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/StratAer/miecalc_aer.F90 (revision 3524) +++ LMDZ6/trunk/libf/phylmd/StratAer/miecalc_aer.F90 (revision 3526) @@ -1,2 +1,5 @@ +! +! $Id$ +! SUBROUTINE MIECALC_AER(tau_strat, piz_strat, cg_strat, tau_strat_wave, tau_lw_abs_rrtm, paprs, debut) Index: LMDZ6/trunk/libf/phylmd/StratAer/minmaxsimple.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/StratAer/minmaxsimple.F90 (revision 3524) +++ LMDZ6/trunk/libf/phylmd/StratAer/minmaxsimple.F90 (revision 3526) @@ -1,4 +1,4 @@ ! -! $Id: minmaxsimple.F90 1910 2013-11-29 08:40:25Z fairhead $ +! $Id$ ! SUBROUTINE minmaxsimple(zq,qmin,qmax,comment) Index: LMDZ6/trunk/libf/phylmd/StratAer/nucleation_tstep_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/StratAer/nucleation_tstep_mod.F90 (revision 3524) +++ LMDZ6/trunk/libf/phylmd/StratAer/nucleation_tstep_mod.F90 (revision 3526) @@ -1,2 +1,5 @@ +! +! $Id$ +! MODULE nucleation_tstep_mod Index: LMDZ6/trunk/libf/phylmd/StratAer/ocs_to_so2.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/StratAer/ocs_to_so2.F90 (revision 3524) +++ LMDZ6/trunk/libf/phylmd/StratAer/ocs_to_so2.F90 (revision 3526) @@ -1,2 +1,5 @@ +! +! $Id$ +! SUBROUTINE ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato) Index: LMDZ6/trunk/libf/phylmd/StratAer/so2_to_h2so4.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/StratAer/so2_to_h2so4.F90 (revision 3524) +++ LMDZ6/trunk/libf/phylmd/StratAer/so2_to_h2so4.F90 (revision 3526) @@ -1,2 +1,5 @@ +! +! $Id$ +! SUBROUTINE SO2_TO_H2SO4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato) Index: LMDZ6/trunk/libf/phylmd/StratAer/strataer_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/StratAer/strataer_mod.F90 (revision 3526) +++ LMDZ6/trunk/libf/phylmd/StratAer/strataer_mod.F90 (revision 3526) @@ -0,0 +1,233 @@ +! $Id$ +MODULE strataer_mod +! This module contains information about strato microphysic model parameters + + IMPLICIT NONE + + ! flag to constraint nucleation rate in a lat/pres box + LOGICAL,SAVE :: flag_nuc_rate_box ! Nucleation rate limit or not to a lat/pres + !$OMP THREADPRIVATE(flag_nuc_rate_box) + REAL,SAVE :: nuclat_min ! min lat to activate nuc rate + REAL,SAVE :: nuclat_max ! max lat to activate nuc rate + REAL,SAVE :: nucpres_min ! min pres to activate nuc rate + REAL,SAVE :: nucpres_max ! max pres to activate nuc rate + !$OMP THREADPRIVATE(nuclat_min, nuclat_max, nucpres_min, nucpres_max) + + ! flag for sulfur emission scenario: (0) background aer ; (1) volcanic eruption ; (2) strato aer injections (SAI) + INTEGER,SAVE :: flag_sulf_emit + !$OMP THREADPRIVATE(flag_sulf_emit) + + ! flag for sulfur emission altitude distribution: (0) gaussian; (1) uniform + INTEGER,SAVE :: flag_sulf_emit_distrib + !$OMP THREADPRIVATE(flag_sulf_emit_distrib) + + !--flag_sulf_emit=1 -- Volcanic eruption(s) + INTEGER,SAVE :: nErupt ! number of eruptions specs + REAL,SAVE :: injdur ! volcanic injection duration + !$OMP THREADPRIVATE(nErupt, injdur) + INTEGER,ALLOCATABLE,SAVE :: year_emit_vol(:) ! year of emission date + INTEGER,ALLOCATABLE,SAVE :: mth_emit_vol(:) ! month of emission date + INTEGER,ALLOCATABLE,SAVE :: day_emit_vol(:) ! day of emission date + !$OMP THREADPRIVATE(year_emit_vol, mth_emit_vol, day_emit_vol) + REAL,ALLOCATABLE,SAVE :: m_aer_emiss_vol(:) ! emitted sulfur mass in kgS, e.g. 7Tg(S)=14Tg(SO2) + REAL,ALLOCATABLE,SAVE :: altemiss_vol(:) ! emission altitude in m + REAL,ALLOCATABLE,SAVE :: sigma_alt_vol(:) ! standard deviation of emission altitude in m + !$OMP THREADPRIVATE(m_aer_emiss_vol, altemiss_vol, sigma_alt_vol) + INTEGER,ALLOCATABLE,SAVE :: ponde_lonlat_vol(:) ! lon/lat ponderation factor + REAL,ALLOCATABLE,SAVE :: xlat_min_vol(:) ! min latitude of volcano in degree + REAL,ALLOCATABLE,SAVE :: xlat_max_vol(:) ! max latitude of volcano in degree + REAL,ALLOCATABLE,SAVE :: xlon_min_vol(:) ! min longitude of volcano in degree + REAL,ALLOCATABLE,SAVE :: xlon_max_vol(:) ! max longitude of volcano in degree + !$OMP THREADPRIVATE(ponde_lonlat_vol, xlat_min_vol, xlat_max_vol, xlon_min_vol, xlon_max_vol) + + !--flag_sulf_emit=2 --SAI + REAL,SAVE :: m_aer_emiss_sai ! emitted sulfur mass in kgS, eg 1e9=1TgS, 1e10=10TgS + REAL,SAVE :: altemiss_sai ! emission altitude in m + REAL,SAVE :: sigma_alt_sai ! standard deviation of emission altitude in m + !$OMP THREADPRIVATE(m_aer_emiss_sai, altemiss_sai, sigma_alt_sai) + REAL,SAVE :: xlat_sai ! latitude of SAI in degree + REAL,SAVE :: xlon_sai ! longitude of SAI in degree + REAL,SAVE :: dlat, dlon ! delta latitude and d longitude of grid in degree + !$OMP THREADPRIVATE(xlat_sai, xlon_sai, dlat, dlon) + + !--flag_sulf_emit=3 -- SAI + REAL,SAVE :: xlat_max_sai ! maximum latitude of SAI in degrees + REAL,SAVE :: xlat_min_sai ! minimum latitude of SAI in degrees + !$OMP THREADPRIVATE(xlat_min_sai,xlat_max_sai) + +CONTAINS + + SUBROUTINE strataer_init() + USE ioipsl_getin_p_mod, ONLY : getin_p + USE print_control_mod, ONLY : lunout + USE mod_phys_lmdz_para, ONLY : is_master + + ! Local var + INTEGER :: ieru + + INTEGER :: i + + WRITE(lunout,*) 'IN STRATAER INIT WELCOME!' + + !Config Key = flag_sulf_emit + !Config Desc = aerosol emission mode + ! - 0 = background aerosol + ! - 1 = volcanic eruption + ! - 2 = geo-ingeneering design + ! - 3 = geo-engineering between two latitudes + !Config Def = 0 + !Config Help = Used in physiq.F + ! + flag_sulf_emit = 0 + nErupt = 0 ! eruption number + injdur = 0 ! init injection duration + CALL getin_p('flag_sulf_emit',flag_sulf_emit) + + IF (flag_sulf_emit==1) THEN ! Volcano + CALL getin_p('nErupt',nErupt) + CALL getin_p('injdur',injdur) + ELSEIF (flag_sulf_emit == 2) THEN ! SAI + CALL getin_p('m_aer_emiss_sai',m_aer_emiss_sai) + CALL getin_p('altemiss_sai',altemiss_sai) + CALL getin_p('sigma_alt_sai',sigma_alt_sai) + CALL getin_p('xlat_sai',xlat_sai) + CALL getin_p('xlon_sai',xlon_sai) + CALL getin_p('flag_sulf_emit_distrib',flag_sulf_emit_distrib) + ELSEIF (flag_sulf_emit == 3) THEN ! SAI between latitudes + CALL getin_p('m_aer_emiss_sai',m_aer_emiss_sai) + CALL getin_p('altemiss_sai',altemiss_sai) + CALL getin_p('sigma_alt_sai',sigma_alt_sai) + CALL getin_p('xlon_sai',xlon_sai) + CALL getin_p('xlat_max_sai',xlat_max_sai) + CALL getin_p('xlat_min_sai',xlat_min_sai) + CALL getin_p('flag_sulf_emit_distrib',flag_sulf_emit_distrib) + ENDIF + + ALLOCATE(year_emit_vol(nErupt),mth_emit_vol(nErupt),day_emit_vol(nErupt)) + ALLOCATE(m_aer_emiss_vol(nErupt),altemiss_vol(nErupt),sigma_alt_vol(nErupt)) + ALLOCATE(xlat_min_vol(nErupt),xlon_min_vol(nErupt)) + ALLOCATE(xlat_max_vol(nErupt),xlon_max_vol(nErupt)) + + year_emit_vol=0 ; mth_emit_vol=0 ; day_emit_vol=0 + m_aer_emiss_vol=0. ; altemiss_vol=0. ; sigma_alt_vol=0. + xlon_min_vol=0. ; xlon_max_vol=0. + xlat_min_vol=0. ; xlat_max_vol=0. + + CALL getin_p('year_emit_vol',year_emit_vol) + CALL getin_p('mth_emit_vol',mth_emit_vol) + CALL getin_p('day_emit_vol',day_emit_vol) + CALL getin_p('m_aer_emiss_vol',m_aer_emiss_vol) + CALL getin_p('altemiss_vol',altemiss_vol) + CALL getin_p('sigma_alt_vol',sigma_alt_vol) + CALL getin_p('xlon_min_vol',xlon_min_vol) + CALL getin_p('xlon_max_vol',xlon_max_vol) + CALL getin_p('xlat_min_vol',xlat_min_vol) + CALL getin_p('xlat_max_vol',xlat_max_vol) + !Config Key = flag_nuc_rate_box + !Config Desc = define or not a box for nucleation rate + ! - F = global nucleation + ! - T = 2D-box for nucleation need nuclat_min, nuclat_max, nucpres_min and + ! nucpres_max + ! to define its bounds. + !Config Def = F + !Config Help = Used in physiq.F + ! + flag_nuc_rate_box = .FALSE. + CALL getin_p('flag_nuc_rate_box',flag_nuc_rate_box) + CALL getin_p('nuclat_min',nuclat_min) + CALL getin_p('nuclat_max',nuclat_max) + CALL getin_p('nucpres_min',nucpres_min) + CALL getin_p('nucpres_max',nucpres_max) + + WRITE(lunout,*) 'IN STRATAER INIT2 year_emit_vol = ',year_emit_vol + WRITE(lunout,*) 'IN STRATAER INIT2 mth_emit_vol = ',mth_emit_vol + WRITE(lunout,*) 'IN STRATAER INIT2 day_emit_vol = ',day_emit_vol + + !IF (is_master) THEN + WRITE(lunout,*) 'IN STRATAER INIT2 year_emit_vol = ',year_emit_vol + WRITE(lunout,*) 'IN STRATAER INIT2 mth_emit_vol=',mth_emit_vol + WRITE(lunout,*) 'IN STRATAER INIT2 day_emit_vol=',day_emit_vol + WRITE(lunout,*) 'IN STRATAER INIT2 =m_aer_emiss_vol',m_aer_emiss_vol + WRITE(lunout,*) 'IN STRATAER INIT2 =altemiss_vol',altemiss_vol + WRITE(lunout,*) 'IN STRATAER INIT2 =sigma_alt_vol',sigma_alt_vol + WRITE(lunout,*) 'IN STRATAER INIT2 xlon_min_vol=',xlon_min_vol + WRITE(lunout,*) 'IN STRATAER INIT2 xlon_max_vol=',xlon_max_vol + WRITE(lunout,*) 'IN STRATAER INIT2 xlat_min_vol=',xlat_min_vol + WRITE(lunout,*) 'IN STRATAER INIT2 xlat_max_vol=',xlat_max_vol + WRITE(lunout,*) 'flag_nuc_rate_box = ',flag_nuc_rate_box + WRITE(lunout,*) 'nuclat_min = ',nuclat_min + WRITE(lunout,*) 'nuclat_max = ',nuclat_max + WRITE(lunout,*) 'nucpres_min = ',nucpres_min + WRITE(lunout,*) 'nucpres_max = ',nucpres_max + WRITE(lunout,*) 'flag_sulf_emit = ',flag_sulf_emit + WRITE(lunout,*) 'injdur = ',injdur + WRITE(lunout,*) 'flag_sulf_emit_distrib = ',flag_sulf_emit_distrib + WRITE(lunout,*) 'nErupt = ',nErupt + WRITE(lunout,*) 'year_emit_vol = ',year_emit_vol + WRITE(lunout,*) 'mth_emit_vol = ',mth_emit_vol + WRITE(lunout,*) 'day_emit_vol = ',day_emit_vol + WRITE(lunout,*) 'm_aer_emiss_vol = ',m_aer_emiss_vol + WRITE(lunout,*) 'altemiss_vol = ',altemiss_vol + WRITE(lunout,*) 'sigma_alt_vol = ',sigma_alt_vol + WRITE(lunout,*) 'xlat_min_vol = ',xlat_min_vol + WRITE(lunout,*) 'xlat_max_vol = ',xlat_max_vol + WRITE(lunout,*) 'xlon_min_vol = ',xlon_min_vol + WRITE(lunout,*) 'xlon_max_vol = ',xlon_max_vol + WRITE(lunout,*) 'm_aer_emiss_sai = ',m_aer_emiss_sai + WRITE(lunout,*) 'altemiss_sai = ',altemiss_sai + WRITE(lunout,*) 'sigma_alt_sai = ',sigma_alt_sai + WRITE(lunout,*) 'xlat_sai = ',xlat_sai + WRITE(lunout,*) 'xlon_sai = ',xlon_sai + WRITE(lunout,*) 'xlat_min_sai = ',xlat_min_sai + WRITE(lunout,*) 'xlat_max_sai = ',xlat_max_sai + !ENDIF + + CALL strataer_ponde_init + WRITE(lunout,*) 'IN STRATAER INT2 END' + + END SUBROUTINE strataer_init + + ! Compute the ponderation to applicate in each grid point for all eruptions and init + ! dlat & dlon variables + SUBROUTINE strataer_ponde_init() + + USE regular_lonlat_mod, ONLY: lon_reg, lat_reg + USE dimphy, ONLY: klon + USE mod_grid_phy_lmdz, ONLY: nbp_lat, nbp_lon + USE print_control_mod, ONLY : lunout + USE YOMCST, ONLY : RPI + + ! local var + REAL :: pi,lat_reg_deg,lon_reg_deg! latitude and d longitude of grid in degree + INTEGER :: ieru, i, j + + ALLOCATE(ponde_lonlat_vol(nErupt)) + + !Compute lon/lat ponderation for injection + 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 + DO ieru=1, nErupt + ponde_lonlat_vol(ieru) = 0 + DO i=1,nbp_lon + lon_reg_deg = lon_reg(i)*180./RPI + DO j=1,nbp_lat + lat_reg_deg = lat_reg(j)*180./RPI + IF ( lat_reg_deg.GE.xlat_min_vol(ieru)-dlat .AND. lat_reg_deg.LT.xlat_max_vol(ieru)+dlat .AND. & + lon_reg_deg.GE.xlon_min_vol(ieru)-dlon .AND. lon_reg_deg.LT.xlon_max_vol(ieru)+dlon ) THEN + ponde_lonlat_vol(ieru) = ponde_lonlat_vol(ieru) + 1 + ENDIF + ENDDO + ENDDO + IF(ponde_lonlat_vol(ieru) == 0) THEN + WRITE(lunout,*) 'STRATAER_INIT ERROR: no grid point found for eruption ieru=',ieru + ENDIF + ENDDO !ieru + WRITE(lunout,*) 'IN STRATAER_INIT ponde_lonlat: ',ponde_lonlat_vol + + END SUBROUTINE strataer_ponde_init + +END MODULE strataer_mod Index: LMDZ6/trunk/libf/phylmd/StratAer/sulfate_aer_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/StratAer/sulfate_aer_mod.F90 (revision 3524) +++ LMDZ6/trunk/libf/phylmd/StratAer/sulfate_aer_mod.F90 (revision 3526) @@ -1,2 +1,5 @@ +! +! $Id$ +! MODULE sulfate_aer_mod Index: LMDZ6/trunk/libf/phylmd/StratAer/traccoag_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/StratAer/traccoag_mod.F90 (revision 3524) +++ LMDZ6/trunk/libf/phylmd/StratAer/traccoag_mod.F90 (revision 3526) @@ -1,2 +1,5 @@ +! +! $Id$ +! MODULE traccoag_mod ! @@ -16,5 +19,5 @@ USE infotrac USE aerophys - USE geometry_mod, ONLY : cell_area + USE geometry_mod, ONLY : cell_area, boundslat USE mod_grid_phy_lmdz USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root @@ -24,4 +27,7 @@ USE phys_local_var_mod, ONLY: stratomask USE YOMCST + USE print_control_mod, ONLY: lunout + USE strataer_mod + USE phys_cal_mod, ONLY : year_len IMPLICIT NONE @@ -52,26 +58,7 @@ ! Local variables !---------------- -! flag for sulfur emission scenario: (0) background aerosol ; (1) volcanic eruption ; (2) stratospheric aerosol injections (SAI) - INTEGER,PARAMETER :: flag_sulf_emit=2 -! -!--flag_sulf_emit=1 --example Pinatubo - INTEGER,PARAMETER :: year_emit_vol=1991 ! year of emission date - INTEGER,PARAMETER :: mth_emit_vol=6 ! month of emission date - INTEGER,PARAMETER :: day_emit_vol=15 ! day of emission date - REAL,PARAMETER :: m_aer_emiss_vol=7.e9 ! emitted sulfur mass in kgS, e.g. 7Tg(S)=14Tg(SO2) - REAL,PARAMETER :: altemiss_vol=17.e3 ! emission altitude in m - REAL,PARAMETER :: sigma_alt_vol=1.e3 ! standard deviation of emission altitude in m - REAL,PARAMETER :: xlat_vol=15.14 ! latitude of volcano in degree - REAL,PARAMETER :: xlon_vol=120.35 ! longitude of volcano in degree - -!--flag_sulf_emit=2 --SAI - REAL,PARAMETER :: m_aer_emiss_sai=1.e10 ! emitted sulfur mass in kgS, eg 1e9=1TgS, 1e10=10TgS - REAL,PARAMETER :: altemiss_sai=17.e3 ! emission altitude in m - REAL,PARAMETER :: sigma_alt_sai=1.e3 ! standard deviation of emission altitude in m - REAL,PARAMETER :: xlat_sai=0.01 ! latitude of SAI in degree - REAL,PARAMETER :: xlon_sai=120.35 ! longitude of SAI in degree - -!--other local variables - INTEGER :: it, k, i, ilon, ilev, itime, i_int + REAL :: m_aer_emiss_vol_daily ! daily injection mass emission + REAL :: sum_emi_so2 ! Test sum of all LON for budg_emi_so2 + INTEGER :: it, k, i, ilon, ilev, itime, i_int, ieru LOGICAL,DIMENSION(klon,klev) :: is_strato ! true = above tropopause, false = below REAL,DIMENSION(klon,klev) :: m_air_gridbox ! mass of air in every grid box [kg] @@ -90,14 +77,46 @@ REAL,DIMENSION(klev) :: zdm ! mass of atm. model layer in kg REAL,DIMENSION(klon,klev) :: dens_aer ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction - REAL :: dlat, dlon ! d latitude and d longitude of grid in degree REAL :: emission ! emission + REAL :: theta_min, theta_max ! for SAI computation between two latitudes + REAL :: dlat_loc IF (is_mpi_root) THEN - PRINT *,'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour + WRITE(lunout,*) 'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour + WRITE(lunout,*) 'IN traccoag flag_sulf_emit: ',flag_sulf_emit + IF (flag_sulf_emit == 1) THEN + WRITE(lunout,*) 'IN traccoag nErupt: ',nErupt + WRITE(lunout,*) 'IN traccoag injdur: ',injdur + WRITE(lunout,*) 'IN traccoag : year_emit_vol',year_emit_vol + WRITE(lunout,*) 'IN traccoag : mth_emit_vol',mth_emit_vol + WRITE(lunout,*) 'IN traccoag : day_emit_vol',day_emit_vol + WRITE(lunout,*) 'IN traccoag : m_aer_emiss_vol',m_aer_emiss_vol + WRITE(lunout,*) 'IN traccoag : altemiss_vol',altemiss_vol + WRITE(lunout,*) 'IN traccoag : sigma_alt_vol',sigma_alt_vol + WRITE(lunout,*) 'IN traccoag : ponde_lonlat_vol',ponde_lonlat_vol + WRITE(lunout,*) 'IN traccoag : xlat_min_vol',xlat_min_vol + WRITE(lunout,*) 'IN traccoag : xlat_max_vol',xlat_max_vol + WRITE(lunout,*) 'IN traccoag : xlon_min_vol',xlon_min_vol + WRITE(lunout,*) 'IN traccoag : xlon_max_vol',xlon_max_vol + ELSEIF (flag_sulf_emit == 2) THEN + WRITE(lunout,*) 'IN traccoag : m_aer_emiss_sai',m_aer_emiss_sai + WRITE(lunout,*) 'IN traccoag : altemiss_sai',altemiss_sai + WRITE(lunout,*) 'IN traccoag : sigma_alt_sai',sigma_alt_sai + WRITE(lunout,*) 'IN traccoag : xlat_sai',xlat_sai + WRITE(lunout,*) 'IN traccoag : xlon_sai',xlon_sai + ELSEIF (flag_sulf_emit == 3) THEN + WRITE(lunout,*) 'IN traccoag : m_aer_emiss_sai',m_aer_emiss_sai + WRITE(lunout,*) 'IN traccoag : altemiss_sai',altemiss_sai + WRITE(lunout,*) 'IN traccoag : sigma_alt_sai',sigma_alt_sai + WRITE(lunout,*) 'IN traccoag : xlat_min_sai',xlat_min_sai + WRITE(lunout,*) 'IN traccoag : xlat_max_sai',xlat_max_sai + WRITE(lunout,*) 'IN traccoag : xlon_sai',xlon_sai + ENDIF + WRITE(lunout,*) 'IN traccoag : flag_nuc_rate_box = ',flag_nuc_rate_box + IF (flag_nuc_rate_box) THEN + WRITE(lunout,*) 'IN traccoag : nuclat_min = ',nuclat_min,', nuclat_max = ',nuclat_max + WRITE(lunout,*) 'IN traccoag : nucpres_min = ',nucpres_min,', nucpres_max = ',nucpres_max + ENDIF ENDIF - - dlat=180./2./FLOAT(nbp_lat) ! d latitude in degree - dlon=360./2./FLOAT(nbp_lon) ! d longitude in degree - + DO it=1, nbtr_bin r_bin(it)=mdw(it)/2. @@ -120,5 +139,5 @@ IF (debutphy .and. is_mpi_root) THEN DO it=1, nbtr_bin - PRINT *,'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')' + WRITE(lunout,*) 'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')' ENDDO ENDIF @@ -170,46 +189,82 @@ !--only emit on day of eruption ! stretch emission over one day of Pinatubo eruption - IF (year_cur==year_emit_vol.AND.mth_cur==mth_emit_vol.AND.day_cur==day_emit_vol) THEN -! - DO i=1,klon - !Pinatubo eruption at 15.14N, 120.35E - IF ( xlat(i).GE.xlat_vol-dlat .AND. xlat(i).LT.xlat_vol+dlat .AND. & - xlon(i).GE.xlon_vol-dlon .AND. xlon(i).LT.xlon_vol+dlon ) THEN -! - PRINT *,'coordinates of volcanic injection point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur -! compute altLMDz - altLMDz(:)=0.0 - DO k=1, klev - zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 - zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg - zdz=zdm(k)/zrho !thickness of layer in m - altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface - ENDDO - !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude) - f_lay_sum=0.0 - DO k=1, klev - f_lay_emiss(k)=0.0 - DO i_int=1, n_int_alt - alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) - f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_vol)* & - & exp(-0.5*((alt-altemiss_vol)/sigma_alt_vol)**2.)* & - & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) - ENDDO - f_lay_sum=f_lay_sum+f_lay_emiss(k) - ENDDO - !correct for step integration error - f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum - !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) - !vertically distributed emission - DO k=1, klev - ! 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 - ENDDO - ENDIF ! emission grid cell - ENDDO ! klon loop - ENDIF ! emission period - + DO ieru=1, nErupt + IF (is_mpi_root) THEN + sum_emi_so2 = 0.0 ! Init sum + ENDIF + IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.& + day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru)+injdur)) THEN + ! + ! daily injection mass emission - NL + m_aer_emiss_vol_daily = m_aer_emiss_vol(ieru)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru))) + WRITE(lunout,*) 'IN traccoag DD m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru), & + 'ponde_lonlat_vol(ieru)=',ponde_lonlat_vol(ieru),'(injdur*ponde_lonlat_vol(ieru))', & + (injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily,'ieru=',ieru + WRITE(lunout,*) 'IN traccoag, dlon=',dlon + DO i=1,klon + !Pinatubo eruption at 15.14N, 120.35E + dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes + WRITE(lunout,*) 'IN traccoag, dlat=',dlat_loc + IF ( xlat(i).GE.xlat_min_vol(ieru)-dlat_loc .AND. xlat(i).LT.xlat_max_vol(ieru)+dlat_loc .AND. & + xlon(i).GE.xlon_min_vol(ieru)-dlon .AND. xlon(i).LT.xlon_max_vol(ieru)+dlon ) THEN + ! + WRITE(lunout,*) 'coordinates of volcanic injection point=',xlat(i),xlon(i),day_cur,mth_cur,year_cur + WRITE(lunout,*) 'DD m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily + ! compute altLMDz + altLMDz(:)=0.0 + DO k=1, klev + zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 + zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg + zdz=zdm(k)/zrho !thickness of layer in m + altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface + ENDDO + + SELECT CASE(flag_sulf_emit_distrib) + + CASE(0) ! Gaussian distribution + !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude) + f_lay_sum=0.0 + DO k=1, klev + f_lay_emiss(k)=0.0 + DO i_int=1, n_int_alt + alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) + f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_vol(ieru))* & + & exp(-0.5*((alt-altemiss_vol(ieru))/sigma_alt_vol(ieru))**2.)* & + & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) + ENDDO + f_lay_sum=f_lay_sum+f_lay_emiss(k) + ENDDO + + CASE(1) ! Uniform distribution + ! In this case, parameter sigma_alt_vol(ieru) is considered to be half the + ! height of the injection, centered around altemiss_vol(ieru) + DO k=1, klev + f_lay_emiss(k)=max(min(altemiss_vol(ieru)+sigma_alt_vol(ieru),altLMDz(k+1))- & + & max(altemiss_vol(ieru)-sigma_alt_vol(ieru),altLMDz(k)),0.)/(2.*sigma_alt_vol(ieru)) + f_lay_sum=f_lay_sum+f_lay_emiss(k) + ENDDO + + END SELECT ! End CASE over flag_sulf_emit_distrib) + + WRITE(lunout,*) "IN traccoag m_aer_emiss_vol=",m_aer_emiss_vol(ieru) + WRITE(lunout,*) "IN traccoag f_lay_emiss=",f_lay_emiss + !correct for step integration error + f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum + !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) + !vertically distributed emission + DO k=1, klev + ! stretch emission over one day of Pinatubo eruption + emission=m_aer_emiss_vol_daily*(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 + ENDDO + sum_emi_so2 = sum_emi_so2 + budg_emi_so2(i) ! Sum all LON + ENDIF ! emission grid cell + ENDDO ! klon loop + WRITE(lunout,*) "IN traccoag (ieru=",ieru,") global sum_emi_so2=",sum_emi_so2 + WRITE(lunout,*) "IN traccoag (ieru=",ieru,") m_aer_emiss_vol_daily=",m_aer_emiss_vol_daily + ENDIF ! emission period + ENDDO ! eruption number + CASE(2) ! stratospheric aerosol injections (SAI) ! @@ -217,10 +272,12 @@ ! SAI standard scenario with continuous emission from 1 grid point at the equator ! SAI emission on single month -! IF ((mth_cur==4 .AND. & ! SAI continuous emission o - IF ( xlat(i).GE.xlat_sai-dlat .AND. xlat(i).LT.xlat_sai+dlat .AND. & + dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes + WRITE(lunout,*) "IN traccoag, dlon=",dlon + WRITE(lunout,*) "IN traccoag, dlat=",dlat_loc + IF ( xlat(i).GE.xlat_sai-dlat_loc .AND. xlat(i).LT.xlat_sai+dlat_loc .AND. & & xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN ! - PRINT *,'coordinates of SAI point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur + WRITE(lunout,*) 'coordinates of SAI point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur ! compute altLMDz altLMDz(:)=0.0 @@ -231,16 +288,33 @@ altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface ENDDO + + SELECT CASE(flag_sulf_emit_distrib) + + CASE(0) ! Gaussian distribution !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude) f_lay_sum=0.0 - DO k=1, klev - f_lay_emiss(k)=0.0 - DO i_int=1, n_int_alt - alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) - f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* & - & exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)* & - & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) - ENDDO - f_lay_sum=f_lay_sum+f_lay_emiss(k) - ENDDO + DO k=1, klev + f_lay_emiss(k)=0.0 + DO i_int=1, n_int_alt + alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) + f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* & + & exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)* & + & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) + ENDDO + f_lay_sum=f_lay_sum+f_lay_emiss(k) + ENDDO + + CASE(1) ! Uniform distribution + f_lay_sum=0.0 + ! In this case, parameter sigma_alt_vol(ieru) is considered to be half + ! the height of the injection, centered around altemiss_sai + DO k=1, klev + f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- & + & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai) + f_lay_sum=f_lay_sum+f_lay_emiss(k) + ENDDO + + END SELECT ! Gaussian or uniform distribution + !correct for step integration error f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum @@ -249,13 +323,89 @@ DO k=1, klev ! stretch emission over whole year (360d) - emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/360./86400. + emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400. 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 ENDDO + ! !emission as monodisperse particles with 0.1um dry radius (BIN21) ! !vertically distributed emission ! DO k=1, klev ! ! stretch emission over whole year (360d) -! emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/360./86400 +! emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400 +! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys +! budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol +! ENDDO + ENDIF ! emission grid cell + ENDDO ! klon loop + + CASE(3) ! --- SAI injection over a single band of longitude and between + ! lat_min and lat_max + + WRITE(lunout,*) 'IN traccoag, dlon=',dlon + DO i=1,klon +! SAI scenario with continuous emission + dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes + WRITE(lunout,*) 'IN traccoag, dlat = ',dlat_loc + theta_min = max(xlat(i)-dlat_loc,xlat_min_sai) + theta_max = min(xlat(i)+dlat_loc,xlat_max_sai) + IF ( xlat(i).GE.xlat_min_sai-dlat_loc .AND. xlat(i).LT.xlat_max_sai+dlat_loc .AND. & + & xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN +! +! compute altLMDz + altLMDz(:)=0.0 + DO k=1, klev + zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 + zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg + zdz=zdm(k)/zrho !thickness of layer in m + altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface + ENDDO + + SELECT CASE(flag_sulf_emit_distrib) + + CASE(0) ! Gaussian distribution + !compute distribution of emission to vertical model layers (based on + !Gaussian peak in altitude) + f_lay_sum=0.0 + DO k=1, klev + f_lay_emiss(k)=0.0 + DO i_int=1, n_int_alt + alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) + f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* & + & exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)* & + & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) + ENDDO + f_lay_sum=f_lay_sum+f_lay_emiss(k) + ENDDO + + CASE(1) ! Uniform distribution + f_lay_sum=0.0 + ! In this case, parameter sigma_alt_vol(ieru) is considered to be half + ! the height of the injection, centered around altemiss_sai + DO k=1, klev + f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- & + & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai) + f_lay_sum=f_lay_sum+f_lay_emiss(k) + ENDDO + + END SELECT ! Gaussian or uniform distribution + + !correct for step integration error + f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum + !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) + !vertically distributed emission + DO k=1, klev + ! stretch emission over whole year (360d) + emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/ & + & year_len/86400.*(sin(theta_max/180.*RPI)-sin(theta_min/180.*RPI))/ & + & (sin(xlat_max_sai/180.*RPI)-sin(xlat_min_sai/180.*RPI)) + 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 + ENDDO + +! !emission as monodisperse particles with 0.1um dry radius (BIN21) +! !vertically distributed emission +! DO k=1, klev +! ! stretch emission over whole year (360d) +! emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400 ! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys ! budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol @@ -291,5 +441,6 @@ IF (mdw(it) .LT. 2.5e-6) THEN !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) & - !assume that particles consist of ammonium sulfate at the surface (132g/mol) and are dry at T = 20 deg. C and 50 perc. humidity + !assume that particles consist of ammonium sulfate at the surface (132g/mol) + !and are dry at T = 20 deg. C and 50 perc. humidity surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) & & *132./98.*dens_aer_dry*4./3.*RPI*(mdw(it)/2.)**3 &