diff --git a/Registry/registry.chem b/Registry/registry.chem index 6cd996156b..cac0a0db9e 100644 --- a/Registry/registry.chem +++ b/Registry/registry.chem @@ -334,8 +334,8 @@ state real eseasc i{dust}jf emis_seas2 1 Z - # volcanic ash emissions # state real - ikjf emis_vol - - - - "Volcanic Emissions" "" -state real e_vash1 ikjf emis_vol 1 Z i{13}hr "E_VASH1" "Volcanic Emissions, bin1" "ug/m2/s" -state real e_vash2 ikjf emis_vol 1 Z i{13}hr "E_VASH2" "Volcanic Emissions, bin2" "ug/m2/s" +state real e_vash1 ikjf emis_vol 1 Z i{13}r "E_VASH1" "Volcanic Emissions, bin1" "ug/m2/s" +state real e_vash2 ikjf emis_vol 1 Z i{13}r "E_VASH2" "Volcanic Emissions, bin2" "ug/m2/s" state real e_vash3 ikjf emis_vol 1 Z i{13}r "E_VASH3" "Volcanic Emissions, bin3" "ug/m2/s" state real e_vash4 ikjf emis_vol 1 Z i{13}r "E_VASH4" "Volcanic Emissions, bin4" "ug/m2/s" state real e_vash5 ikjf emis_vol 1 Z i{13}r "E_VASH5" "Volcanic Emissions, bin5" "ug/m2/s" @@ -345,6 +345,9 @@ state real e_vash8 ikjf emis_vol 1 Z i{13}r state real e_vash9 ikjf emis_vol 1 Z i{13}r "E_VASH9" "Volcanic Emissions, bin9" "ug/m2/s" state real e_vash10 ikjf emis_vol 1 Z i{13}r "E_VASH10" "Volcanic Emissions, bin10" "ug/m2/s" state real e_vso2 ikjf emis_vol 1 Z i{13}r "E_VSO2" "Volcanic Emissions, SO2" "mol/m2/h" +state real e_vsulf ikjf emis_vol 1 Z i{13}r "E_VSULF" "Volcanic Emissions, SULF" "mol/m2/h" +state real e_qv ikjf emis_vol 1 Z i{13}r "E_QV" "Volcanic Emissions, water vapor" "kg/m2/s" +state real e_start ikjf emis_vol 1 Z i{13}r "E_START" "Volcanic Emissions, start" "" # # biomassburning emission arrays (used by fire plume rise) state real - ikjf ebu - - - - "Biomass burnung Emissions" "" @@ -705,7 +708,6 @@ state real br_rto ikj misc 1 - r "BRN # # following used for volcanic eruptions only # -state real ash_fall ij misc 1 - r "ASH_FALL" "VOLCANIC ASH FALL" "kg/m2" state real erup_beg ij misc 1 - i{13} "ERUP_BEG" "START TIME OF ERUPTION" "?" state real erup_end ij misc 1 - i{13} "ERUP_END" "END TIME OF ERUPTION" "?" # @@ -1100,8 +1102,21 @@ state integer num_vert_mix - misc - - r "n # # Wet deposition # -state real wd_so4_sc ij misc 1 - rdu "wd_so4_sc" "SO4 surface wet deposition, accumulated (Sc)" "mmol/m2" -state real wd_no3_sc ij misc 1 - rdu "wd_no3_sc" "NO3 surface wet deposition, accumulated (Sc)" "mmol/m2" +state real wd_so4_sc ij misc 1 - rdu "wd_so4_sc" "SO4 surface wet deposition, accumulated (Sc)" "mmol/m2" +state real wd_no3_sc ij misc 1 - rdu "wd_no3_sc" "NO3 surface wet deposition, accumulated (Sc)" "mmol/m2" + + +state real - i%jf volc_diags - - - - "Volcanic ash, so2 and sulf diagstics" "" +state real wd_ash_sc i%jf volc_diags 1 Z rdu "wd_ash_sc" "ash surface large scale wet deposition, accumulated" "g/m2" +state real wd_ash_cu i%jf volc_diags 1 Z rdu "wd_ash_cu" "ash surface wet deposition, accumulated (Cumul)" "ug/m2" +state real ash_drydep i%jf volc_diags 1 Z rdu "ash_drydep" "Accumulated volcanic ash dry deposition" "kg/m2" +state real sulf_drydep i%jf volc_diags 1 Z rdu "sulf_drydep" "Accumulated volcanic sulfate dry deposition" "mol/m2" +state real so2_drydep i%jf volc_diags 1 Z rdu "so2_drydep" "Accumulated volcanic so2 dry deposition" "mol/m2" +state real sulf_grav_setl i%jf volc_diags 1 Z rdu "sulf_grav_setl" "Accumulated volcanic sulfate grav. settling" "kg/m2" +state real so2_oh_loss i%jf volc_diags 1 Z rdu "so2_oh_loss" "Accumulated volcanic so2 loss due to OH" "kgS" +state real so2_h2o2_loss i%jf volc_diags 1 Z rdu "so2_h2o2_loss" "Accumulated volcanic so2 loss due to H2O2" "kgS" +state real ash_fall i%jf volc_diags 1 Z rdu "ash_fall" "volcanic asg fall" "kg/m2" + # added wet deposition totals for NH4 and OA for MOZART coupled version state real wd_nh4_sc ij misc 1 - rdu "wd_nh4_sc" "NH4 surface wet deposition, accumulated (Sc)" "mmol/m2" state real wd_oa_sc ij misc 1 - rdu "wd_oa_sc" "Organics surface wet deposition, accumulated (Sc)" "mmol/m2" @@ -4076,6 +4091,7 @@ package ecrimechtno emiss_opt==20 - emis_ant: # package vash emiss_opt_vol==1 - emis_vol:e_vash1,e_vash2,e_vash3,e_vash4,e_vash5,e_vash6,e_vash7,e_vash8,e_vash9,e_vash10 package vashso2 emiss_opt_vol==2 - emis_vol:e_vash1,e_vash2,e_vash3,e_vash4,e_vash5,e_vash6,e_vash7,e_vash8,e_vash9,e_vash10,e_vso2 +package vashso2qv emiss_opt_vol==3 - emis_vol:e_vash1,e_vash2,e_vash3,e_vash4,e_vash5,e_vash6,e_vash7,e_vash8,e_vash9,e_vash10,e_vso2,e_qv,e_vsulf,e_start;volc_diags:wd_ash_sc,wd_ash_cu,ash_drydep,sulf_drydep,so2_drydep,sulf_grav_setl,so2_oh_loss,so2_h2o2_loss,ash_fall # diff --git a/chem/CMakeLists.txt b/chem/CMakeLists.txt index 8c329a1397..9dfaead80d 100644 --- a/chem/CMakeLists.txt +++ b/chem/CMakeLists.txt @@ -152,6 +152,7 @@ target_sources( module_cbmz_rodas_prep.F module_ctrans_grell.F module_gocart_chem.F + module_volc_chem.F module_input_tracer.F module_lightning_nox_driver.F module_lightning_nox_ott.F diff --git a/chem/Makefile b/chem/Makefile index 9bf3acd38f..6ddc893c02 100755 --- a/chem/Makefile +++ b/chem/Makefile @@ -80,6 +80,7 @@ MODULES = \ module_gocart_dust.o \ module_gocart_dust_afwa.o \ module_gocart_seasalt.o \ + module_volc_chem.o \ module_uoc_dust.o \ module_qf03.o \ module_soilpsd.o \ diff --git a/chem/aerosol_driver.F b/chem/aerosol_driver.F index 92b424d7fd..8422cfe43a 100755 --- a/chem/aerosol_driver.F +++ b/chem/aerosol_driver.F @@ -5,6 +5,7 @@ ! ! 10/12/2011 - Ravan Ahmadov (NOAA) updated to include the RACM_SOA_VBS option ! 10/08/2014 - Kai Wang and Yang Zhang (NCSU) updated to include the CB05_MADE/SORGAM and CB05_MADE/VBS options +! 01/13/2025 - Alexander Ukhov (KAUST) added calculation of PM2.5 and PM10 for CHEM_VASH,CHEM_VOLC,CHEM_VOLC_4BIN ! SUBROUTINE aerosols_driver (id,curr_secs,ktau,dtstep,ktauc, & config_flags,dtstepc,dx, & @@ -426,6 +427,7 @@ SUBROUTINE sum_pm_driver ( config_flags, & USE module_gocart_aerosols, only: sum_pm_gocart USE module_aerosols_soa_vbs, only: sum_pm_soa_vbs USE module_aerosols_sorgam_vbs, only: sum_pm_sorgam_vbs + USE module_volc_chem, only: sum_pm_volc IMPLICIT NONE @@ -499,6 +501,15 @@ SUBROUTINE sum_pm_driver ( config_flags, & ! sum_pm_select: SELECT CASE(config_flags%chem_opt) + ! 01/13/2025 - A. Ukhov + case (chem_vash,chem_volc,chem_volc_4bin) + CALL wrf_debug(15,'sum_pm_driver: calling sum_pm_volc') + CALL sum_pm_volc ( & + alt, chem,pm2_5_dry, pm10, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte ) + CASE (GOCART_SIMPLE,GOCARTRACM_KPP,GOCARTRADM2,MOZCART_KPP,T1_MOZCART_KPP) CALL wrf_debug(15,'sum_pm_driver: calling sum_pm_gocart') CALL sum_pm_gocart ( & diff --git a/chem/chem_driver.F b/chem/chem_driver.F index 485574ffe6..9beeff0958 100755 --- a/chem/chem_driver.F +++ b/chem/chem_driver.F @@ -4,7 +4,8 @@ ! William Gustafson (PNNL), Marc Salzmann (GFDL), and Georg Grell ! 10/12/2011 - Ravan Ahmadov (NOAA) updated to include the RACM_SOA_VBS option ! 10/08/2014 - Kai Wang and Yang Zhang (NCSU) updated to include the CB05_MADE/SORGAM and CB05_MADE/VBS options -! +! 28/10/2024 - A. Ukhov (KAUST) added diagnostic for large and convective scale washout of volcanic ash +! A. Ukhov (KAUST) added calculation of optical properties of ash and sulf aerosol. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine chem_driver ( grid , config_flags & @@ -42,6 +43,7 @@ subroutine chem_driver ( grid , config_flags & USE module_wetscav_driver, only: wetscav_driver USE module_wetdep_ls, only:wetdep_ls USE module_uoc_dustwd ! Claudia, 3 April 2014 [mklose 03082015] + USE module_vash_settling, only: wetdep_ls_volc ! A. Ukhov (KAUST) 24 Oct 2024 USE module_input_chem_data, only: last_chem_time, & mozcart_lbc_set, & bdy_chem_value_top_pv, & @@ -815,7 +817,7 @@ end SUBROUTINE sum_pm_driver .or. config_flags%tracer_opt > 0 .or. config_flags%bio_emiss_opt == MEGAN2_CLM ) then call wrf_debug(15,'calling emissions driver') - call emissions_driver(grid%id,ktau,grid%dt,grid%DX, & + call emissions_driver(grid%id,ktau,grid%dt,grid%DX,grid%DY, & adapt_step_flag, curr_secs, & grid%plumerisefire_frq,grid%stepfirepl, & grid%bioemdt,grid%stepbioe, & @@ -830,7 +832,7 @@ end SUBROUTINE sum_pm_driver grid%firesize_agef,grid%firesize_agsv,grid%firesize_aggr, & grid%u10,grid%v10,grid%ivgtyp,grid%isltyp,grid%gsw,grid%vegfra,grid%rmol, & grid%ust,grid%znt,grid%dms_0,grid%erup_beg,grid%erup_end, & - grid%xland,grid%xlat,grid%xlong, & + grid%xland,grid%xlat,grid%xlong,grid%msftx,grid%msfty, & z_at_w,zmid,grid%smois,dustin,seasin, & grid%sebio_iso,grid%sebio_oli,grid%sebio_api,grid%sebio_lim, & grid%sebio_xyl,grid%sebio_hc3,grid%sebio_ete,grid%sebio_olt, & @@ -915,9 +917,9 @@ end SUBROUTINE sum_pm_driver ! if( do_photstep .and. & config_flags%chem_opt /= CHEM_TRACER .and. & - config_flags%chem_opt /= CHEM_VASH .and. & - config_flags%chem_opt /= CHEM_VOLC .and. & - config_flags%chem_opt /= CHEM_VOLC_4BIN .and. & + !config_flags%chem_opt /= CHEM_VASH .and. & ! A. Ukhov. 28 Oct 2024 + !config_flags%chem_opt /= CHEM_VOLC .and. & + !config_flags%chem_opt /= CHEM_VOLC_4BIN .and. & config_flags%chem_opt /= DUST .and. & config_flags%chem_opt /= CHEM_TRACE2 .and. & config_flags%chem_opt /= CO2_TRACER .and. & @@ -1053,7 +1055,8 @@ end SUBROUTINE sum_pm_driver grid%ccn1, grid%ccn2, grid%ccn3, grid%ccn4, grid%ccn5, grid%ccn6, & grid%qndropsource,grid%ivgtyp,grid%tsk,grid%gsw,grid%vegfra,pbl_h, & grid%rmol,grid%ust,grid%znt,grid%xlat,grid%xlong, & - zmid,z_at_w,grid%xland,grid%ash_fall, & + zmid,z_at_w,grid%xland, & + grid%volc_diags, & ! A. Ukhov 18 Oct 2024. Volcanic diagnostics grid%h2oaj,grid%h2oai,grid%nu3,grid%ac3,grid%cor3,grid%asulf, & grid%ahno3,grid%anh3,grid%cvaro1,grid%cvaro2,grid%cvalk1,grid%cvole1,& grid%cvapi1,grid%cvapi2,grid%cvlim1,grid%cvlim2,grid%dep_vel_o3, & @@ -1157,8 +1160,9 @@ end SUBROUTINE sum_pm_driver U_phy,V_phy,t_phy,moist,dz8w, & p_phy,XLV,CP,G,r_v, & z_at_w,grid%cu_co_ten, & - grid%wd_no3_cu,grid%wd_so4_cu, & + grid%wd_no3_cu,grid%wd_so4_cu, & grid%wd_nh4_cu,grid%wd_oa_cu, & + volc_diags, & !A. Ukhov 18 Oct 2024. Volcanic diagnostics grid%wd_so2_cu, grid%wd_sulf_cu, grid%wd_hno3_cu, grid%wd_nh3_cu, & grid%wd_cvasoa_cu, grid%wd_cvbsoa_cu, grid%wd_asoa_cu, grid%wd_bsoa_cu, & grid%k22_shallow,grid%kbcon_shallow,grid%ktop_shallow,grid%xmb_shallow, & @@ -1186,8 +1190,9 @@ end SUBROUTINE sum_pm_driver U_phy,V_phy,t_phy,moist,dz8w, & p_phy,XLV,CP,G,r_v, & z_at_w, grid%cu_co_ten, & - grid%wd_no3_cu,grid%wd_so4_cu, & + grid%wd_no3_cu,grid%wd_so4_cu, & grid%wd_nh4_cu,grid%wd_oa_cu, & + volc_diags, & ! A. Ukhov 18 Oct 2024. Volcanic diagnostics grid%wd_so2_cu, grid%wd_sulf_cu, grid%wd_hno3_cu, grid%wd_nh3_cu, & grid%wd_cvasoa_cu, grid%wd_cvbsoa_cu, grid%wd_asoa_cu, grid%wd_bsoa_cu, & grid%k22_shallow,grid%kbcon_shallow,grid%ktop_shallow,grid%xmb_shallow, & @@ -1219,7 +1224,7 @@ end SUBROUTINE sum_pm_driver if( do_chemstep .and. & config_flags%chem_opt /= CHEM_TRACER .and. & config_flags%chem_opt /= CHEM_VASH .and. & - config_flags%chem_opt /= CHEM_VOLC .and. & + !config_flags%chem_opt /= CHEM_VOLC .and. & ! A. Ukhov. 28 Oct 2024 config_flags%chem_opt /= CHEM_VOLC_4BIN .and. & config_flags%chem_opt /= DUST .and. & config_flags%chem_opt /= CHEM_TRACE2 .and. & @@ -1241,7 +1246,7 @@ end SUBROUTINE sum_pm_driver call mechanism_driver(grid%id,curr_secs,ktau,grid%dt,grid%ktauc,dtstepc,config_flags, & grid%gmt,ijulian,t_phy,moist,p8w,t8w,grid%gd_cldfr, & - p_phy,chem,rho,dz8w,grid%dx,g, & + p_phy,chem,rho,dz8w,grid%dx,grid%dy,grid%msftx,grid%msfty,g, & zmid,z_at_w,grid%xlat,grid%xlong, & vdrog3,vcsulf_old,vcso2_old,vch2o2_old,grid%ttday,grid%tcosz, & grid%ph_macr,grid%ph_o31d,grid%ph_o33p,grid%ph_no2, & @@ -1257,6 +1262,7 @@ end SUBROUTINE sum_pm_driver grid%xylp,grid%apip,grid%isop,grid%hc3p,grid%ethp,grid%o3p,grid%tco3, & grid%mo2,grid%o1d,grid%olnn,grid%rpho,grid%xo2, & grid%ketp,grid%olnd, & + grid%volc_diags, & ! A. Ukhov 18 Oct 2024. Volcanic diagnostics ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite,jts,jte,kts,kte ) @@ -1672,25 +1678,29 @@ end SUBROUTINE sum_pm_driver ! since we do not have h2o2 as a variable, pass in p_h2o2 as zero ! will have to use backgrund value ! - if(config_flags%chem_opt == CHEM_VOLC)then - CALL wrf_debug(15,'gocart so2-so4 conversion') - CALL so2so4(0,chem,p_so2,p_sulf,p_h2o2,p_QC,T_PHY,MOIST, & - grid%qc_cu, grid%gd_cldfr, config_flags%cu_diag, & - NUM_CHEM,NUM_MOIST, & - ids,ide, jds,jde, kds,kde, & - ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte ) - endif ! ! now do wet removal; first LS if there is no explicit aqeous phase -! - if(config_flags%wetscav_onoff<0)then - call wrf_debug(15,'calculate LS wet deposition') - call wetdep_ls(grid%dt,chem,grid%rainncv,moist,rho,num_moist, & + if(config_flags%wetscav_onoff<0)then + call wrf_debug(15,'calculate LS wet deposition') + if(config_flags%chem_opt == CHEM_VOLC .or. & + config_flags%chem_opt == chem_volc_4bin .or. & + config_flags%chem_opt == chem_vash) then + + !A. Ukhov (KAUST). large scale scaveging for Ash, volcanic SO2 and Sulf + call wetdep_ls_volc(grid%dt,chem,grid%rainncv,moist,rho,num_moist, & + num_chem,numgas,dz8w,vvel,grid%chem_opt, & + grid%wd_so2_sc,grid%wd_sulf_sc, & + grid%volc_diags, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte ) + else + call wetdep_ls(grid%dt,chem,grid%rainncv,moist,rho,num_moist, & num_chem,numgas,dz8w,vvel,grid%chem_opt, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) + endif endif ! ! Sum up the aerosol mass for radiation and diagnostic purposes. Unlike diff --git a/chem/chemics_init.F b/chem/chemics_init.F index 9856ba9dc5..a4dcda31f2 100755 --- a/chem/chemics_init.F +++ b/chem/chemics_init.F @@ -1790,6 +1790,53 @@ subroutine chem_init (id,chem,emis_ant,scalar,dt,bioemdt,photdt,chemdt,stepbioe, ! if(i.eq.19.and.j.eq.19)write(0,*)TCOSZ(i,j),ttday(i,j),julday, gmtp, sza, cosszax,xlonn,rlat enddo enddo + + ! 17 Oct. 2024. A. Ukhov. + CASE (CHEM_VOLC) + CALL wrf_debug(15,'call CHEM_VOLC initialization') + if(config_flags%phot_opt .NE. 0 )then + call wrf_error_fatal("CHEM_VOLC initialization, phot_opt MUST BE ZERO") + endif + CALL wrf_debug(15,'call volcanic aerosols initialization') + if(config_flags%chem_in_opt == 0 )then + if( .NOT. config_flags%restart )then + do j=jts,jte + do k=kts,kte + do i=its,ite + do n=1,num_chem + chem(i,k,j,n)=0. + enddo + enddo + enddo + enddo + endif + endif + + ! scale background oh in dependence on average zenith angle + ! this is done since background values are only available as average/month. + ! + ndystep=86400/ifix(dt) + do j=jts,jte + do i=its,ite + tcosz(i,j)=0. + ttday(i,j)=0. + rlat=xlat(i,j)*3.1415926535590/180. + xlonn=xlong(i,j) + do n=1,ndystep + xtime=n*dt/60. + ixhour=ifix(gmt+.01)+ifix(xtime/60.) + xhour=float(ixhour) + xmin=60.*gmt+(xtime-xhour*60.) + gmtp=mod(xhour,24.) + gmtp=gmtp+xmin/60. + CALL szangle(1, 1, julday, gmtp, sza, cosszax,xlonn,rlat) + TCOSZ(i,j)=TCOSZ(I,J)+cosszax(1,1) + if(cosszax(1,1).gt.0.)ttday(i,j)=ttday(i,j)+dt + enddo +! if(i.eq.19.and.j.eq.19)write(0,*)TCOSZ(i,j),ttday(i,j),julday, gmtp, sza, cosszax,xlonn,rlat + enddo + enddo + CASE (MOZCART_KPP,T1_MOZCART_KPP) CALL wrf_debug(15,'MOZCART dust initialization') ch_dust(:,:) = 0.8D-9 diff --git a/chem/depend.chem b/chem/depend.chem index ffe9704bfc..86fb038e81 100644 --- a/chem/depend.chem +++ b/chem/depend.chem @@ -65,7 +65,7 @@ module_gocart_chem.o: module_data_gocartchem.o module_phot_mad.o module_gocart_so2so4.o: -module_vash_settling.o: +module_vash_settling.o: module_data_sorgam.o module_gocart_settling.o: ../phys/module_data_gocart_dust.o module_data_gocart_seas.o @@ -73,6 +73,8 @@ module_gocart_aerosols.o: module_data_gocartchem.o module_gocart_drydep.o: +module_volc_chem.o : module_gocart_chem.o + module_wetdep_ls.o: module_gocart_dmsemis.o: module_data_gocartchem.o @@ -340,7 +342,7 @@ module_upper_bc_driver.o: module_tropopause.o chem_driver.o: module_radm.o ../dyn_em/module_convtrans_prep.o module_chem_utilities.o module_data_radm2.o module_dep_simple.o module_bioemi_simple.o module_vertmx_wrf.o module_phot_mad.o module_aerosols_sorgam.o module_aerosols_soa_vbs.o module_aerosols_soa_vbs_het.o module_aerosols_sorgam_vbs.o module_data_cbmz.o module_cbmz.o module_wetscav_driver.o dry_dep_driver.o emissions_driver.o module_input_tracer.o module_input_tracer_data.o module_tropopause.o module_upper_bc_driver.o module_ctrans_grell.o module_data_soa_vbs.o module_data_soa_vbs_het.o module_aer_opt_out.o module_data_sorgam.o module_gocart_so2so4.o ../phys/module_cu_camzm_driver.o module_cam_mam_gas_wetdep_driver.o module_dust_load.o module_chem_cup.o ../share/module_trajectory.o ../share/module_chem_share.o -aerosol_driver.o: module_data_sorgam.o module_aerosols_sorgam.o module_data_soa_vbs.o module_aerosols_soa_vbs.o module_data_soa_vbs_het.o module_aerosols_soa_vbs_het.o module_aerosols_sorgam_vbs.o module_mosaic_driver.o +aerosol_driver.o: module_data_sorgam.o module_aerosols_sorgam.o module_data_soa_vbs.o module_aerosols_soa_vbs.o module_data_soa_vbs_het.o module_aerosols_soa_vbs_het.o module_aerosols_sorgam_vbs.o module_mosaic_driver.o module_volc_chem.o module_sorgam_aqchem.o: ../share/module_ctrans_aqchem.o module_data_sorgam.o @@ -350,13 +352,13 @@ cloudchem_driver.o: module_mosaic_cloudchem.o module_sorgam_cloudchem.o module_s photolysis_driver.o: module_phot_mad.o module_phot_fastj.o module_ftuv_driver.o module_phot_tuv.o -mechanism_driver.o: module_data_radm2.o module_radm.o module_aerosols_sorgam.o module_aerosols_soa_vbs.o module_aerosols_soa_vbs_het.o module_data_cbmz.o module_cbmz.o +mechanism_driver.o: module_data_radm2.o module_radm.o module_aerosols_sorgam.o module_aerosols_soa_vbs.o module_aerosols_soa_vbs_het.o module_data_cbmz.o module_cbmz.o module_volc_chem.o optical_driver.o: module_optical_averaging.o module_peg_util.o module_data_mosaic_therm.o emissions_driver.o: module_add_emiss_burn.o module_data_radm2.o module_radm.o module_bioemi_simple.o module_bioemi_beis314.o module_bioemi_megan2.o module_emissions_anthropogenics.o module_cbmz_addemiss.o module_cb05_addemiss.o module_mosaic_addemiss.o module_aerosols_sorgam.o module_aerosols_soa_vbs.o module_aerosols_soa_vbs_het.o module_aerosols_sorgam_vbs.o module_plumerise1.o module_gocart_dust.o module_gocart_dust_afwa.o module_uoc_dust.o module_gocart_seasalt.o module_ghg_fluxes.o module_lightning_nox_driver.o module_cam_mam_addemiss.o -dry_dep_driver.o: module_data_radm2.o module_aer_drydep.o module_dep_simple.o module_aerosols_sorgam.o module_aerosols_soa_vbs.o module_aerosols_soa_vbs.o module_aerosols_sorgam_vbs.o module_mosaic_drydep.o module_mixactivate_wrappers.o ../phys/module_mixactivate.o module_cam_mam_drydep.o ../phys/module_data_cam_mam_asect.o ../phys/module_data_cam_mam_aero.o ../phys/module_cam_support.o +dry_dep_driver.o: module_data_radm2.o module_aer_drydep.o module_dep_simple.o module_aerosols_sorgam.o module_aerosols_soa_vbs.o module_aerosols_soa_vbs.o module_aerosols_sorgam_vbs.o module_mosaic_drydep.o module_mixactivate_wrappers.o ../phys/module_mixactivate.o module_cam_mam_drydep.o ../phys/module_data_cam_mam_asect.o ../phys/module_data_cam_mam_aero.o ../phys/module_cam_support.o module_vash_settling.o module_wetscav_driver.o: module_mosaic_wetscav.o module_aerosols_sorgam.o module_aerosols_sorgam_vbs.o module_mozcart_wetscav.o ../phys/module_data_cam_mam_aero.o module_cam_mam_wetscav.o module_aerosols_soa_vbs.o module_prep_wetscav_sorgam.o diff --git a/chem/dry_dep_driver.F b/chem/dry_dep_driver.F index b5c2f66501..41f8b43c7c 100755 --- a/chem/dry_dep_driver.F +++ b/chem/dry_dep_driver.F @@ -3,6 +3,7 @@ ! VERTMX was originally developed by Mariusz Pagowski and modified by ! Richard C. Easter (PNNL) ! 10/11/2011 - Ravan Ahmadov (NOAA) updated to include the RACM_SOA_VBS option +! 10/30/2024 - A. Ukhov (KAUST) added dry deposition of volcanic ash, so2, sulf ! !WRF:MODEL_LAYER:CHEMICS ! @@ -13,22 +14,22 @@ MODULE module_dry_dep_driver subroutine dry_dep_driver(id,curr_secs,ktau,dtstep,config_flags, & gmt,julday,t_phy,moist,scalar,p8w,t8w,w,alt, & - p_phy,chem,tracer,rho_phy,dz8w,rh,exch_h,hfx,dx, & + p_phy,chem,tracer,rho_phy,dz8w,rh,exch_h,hfx,dx, & cldfra, cldfra_old,raincv,seasin,dustin, & ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource, & ivgtyp,tsk,gsw,vegfra,pbl,rmol,ust,znt,xlat,xlong,z,z_at_w,& - xland,ash_fall,h2oaj,h2oai,nu3,ac3,cor3,asulf,ahno3, & + xland,volc_diags,h2oaj,h2oai,nu3,ac3,cor3,asulf,ahno3, & anh3,cvaro1,cvaro2, & cvalk1,cvole1,cvapi1,cvapi2,cvlim1,cvlim2,dep_vel_o3, & ddlen,ddflx, & emis_ant,ebu_in, & - sf_urban_physics,numgas,current_month,dvel,snowh, & + sf_urban_physics,numgas,current_month,dvel,snowh, & dustdrydep_1,dustdrydep_2,dustdrydep_3, & dustdrydep_4,dustdrydep_5, & - depvelocity, & + depvelocity, & dustgraset_1,dustgraset_2,dustgraset_3, & dustgraset_4,dustgraset_5, & - setvel_1,setvel_2,setvel_3,setvel_4,setvel_5, imod, & + setvel_1,setvel_2,setvel_3,setvel_4,setvel_5, imod, & is_CAMMGMP_used, & dep_vel,num_vert_mix, & ids,ide, jds,jde, kds,kde, & @@ -44,7 +45,7 @@ subroutine dry_dep_driver(id,curr_secs,ktau,dtstep,config_flags, & USE module_data_sorgam USE module_aerosols_sorgam USE module_gocart_settling - USE module_vash_settling + USE module_vash_settling, only: volc_ash_sulf_so2_drydep_driver,sulf_settling_driver,vash_settling_driver ! A.Ukhov USE module_gocart_drydep USE module_mosaic_drydep, only: mosaic_drydep_driver USE module_mixactivate_wrappers, only: mosaic_mixactivate, sorgam_mixactivate, & @@ -127,7 +128,12 @@ subroutine dry_dep_driver(id,curr_secs,ktau,dtstep,config_flags, & xlat, & xlong, & snowh, & - xland,znt,raincv,ash_fall + xland,znt,raincv + + ! A. Ukhov Volcanic diagnostics + REAL, DIMENSION( ims:ime, jms:jme,1:9), INTENT(INOUT ) :: volc_diags + + REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT(INOUT ) :: & cldfra, & ! cloud fraction current timestep @@ -265,6 +271,7 @@ subroutine dry_dep_driver(id,curr_secs,ktau,dtstep,config_flags, & config_flags%chem_opt /= GHG_TRACER .and. & config_flags%chem_opt /= CHEM_VASH .and. & config_flags%chem_opt /= CHEM_VOLC_4BIN .and. & + config_flags%chem_opt /= CHEM_VOLC .and. & ! A. Ukhov config_flags%chem_opt /= DUST .and. & config_flags%chem_opt /= GOCART_SIMPLE .and. & config_flags%chem_opt /= GOCARTRACM_KPP )THEN @@ -439,6 +446,28 @@ subroutine dry_dep_driver(id,curr_secs,ktau,dtstep,config_flags, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) + !A. UKhov Dry deposition of volcanic SO2, sulf, and ash + ELSEIF ( config_flags%chem_opt == CHEM_VOLC) then + call wesely_driver(id,ktau,dtstep, & + config_flags,current_month, & + gmt,julday,t_phy,moist,p8w,t8w,raincv, & + p_phy,chem,rho_phy,dz8w,ddvel,aer_res_def,aer_res_zcen, & + ivgtyp,tsk,gsw,vegfra,pbl,rmol,ust,znt,xlat,xlong,z,z_at_w,& + snowh, numgas, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte ) + + call volc_ash_sulf_so2_drydep_driver(dtstep, & + config_flags,numgas, & + t_phy,moist,p8w,t8w,rmol,aer_res_def, & + p_phy,chem,rho_phy,dz8w,ddvel,xland,hfx, & + tsk,pbl,ust,znt, & + volc_diags, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte ) + ELSE !Set dry deposition velocity to zero when using the !chemistry tracer mode. @@ -1395,23 +1424,39 @@ subroutine dry_dep_driver(id,curr_secs,ktau,dtstep,config_flags, & dustin,seasin,dx,g, & dustgraset_1,dustgraset_2,dustgraset_3, & dustgraset_4,dustgraset_5, & - setvel_1,setvel_2,setvel_3,setvel_4,setvel_5, imod, & - ids,ide, jds,jde, kds,kde, & - ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte ) - CASE (CHEM_VASH, CHEM_VOLC, CHEM_VOLC_4BIN) - CALL wrf_debug(15,'call vash settling routine') - call vash_settling_driver(dtstep,config_flags,t_phy,moist, & - chem,rho_phy,dz8w,p8w,p_phy, & - ash_fall,dx,g, & + setvel_1,setvel_2,setvel_3,setvel_4,setvel_5, imod, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) - CASE DEFAULT - CALL wrf_debug(15,'no settling routine') END SELECT settling_select ENDIF + +! A. Ukhov, 12 Oct 2024. added volc_diags +vash_settling_select: SELECT CASE(config_flags%chem_opt) +CASE (CHEM_VASH, CHEM_VOLC, CHEM_VOLC_4BIN) + CALL wrf_debug(15,'call vash grav. settling routine') + call vash_settling_driver(dtstep,config_flags,t_phy,moist, & + chem,rho_phy,dz8w,p8w,p_phy,volc_diags,dx,g, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte ) + + !A. Ukhov, 12 Oct 2024 + if (config_flags%chem_opt.eq.CHEM_VOLC) then + CALL wrf_debug(15,'call sulf. aerosol grav. settling routine') + call sulf_settling_driver(dtstep,config_flags,t_phy,moist, & + chem,rho_phy,dz8w,p8w,p_phy,dx,g,z_at_w,volc_diags, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte) + endif + +CASE DEFAULT + CALL wrf_debug(15,'no settling routine') +END SELECT vash_settling_select + + CALL wrf_debug(15,'end of dry_dep_driver') END SUBROUTINE dry_dep_driver diff --git a/chem/emissions_driver.F b/chem/emissions_driver.F index 47b1e6f4b2..2394518d1a 100644 --- a/chem/emissions_driver.F +++ b/chem/emissions_driver.F @@ -5,14 +5,19 @@ ! Saulo Freitas (CPTEC), and Georg Grell ! ! -! A. Ukhov, 11 March 2021, remove unused parameters in gocart_dust_driver(), +! A. Ukhov (KAUST), 11 March 2021, remove unused parameters in gocart_dust_driver(), ! gocart_dust_afwa_driver(), and uoc_dust_driver() subroutines. +! A. Ukhov (KAUST), 3 Oct 2024, added mapfac_mx and mapfac_my for +! precise calculation of cell area. This is needed for volcanic eruption. +! Corrected fractions for calculation of contribution of vash into PM2.5 and PM10. +! Added new emiss_opt_vol=3 for emissions of volcanic sulfate and water vapor. + MODULE module_emissions_driver IMPLICIT NONE CONTAINS - subroutine emissions_driver(id,ktau,dtstep,DX, & + subroutine emissions_driver(id,ktau,dtstep,DX,DY, & adapt_step_flag,curr_secs, & plumerisefire_frq,stepfirepl, & bioemdt,stepbioe, & @@ -26,7 +31,7 @@ subroutine emissions_driver(id,ktau,dtstep,DX, & firesize_agsv,firesize_aggr, & u10,v10,ivgtyp,isltyp,gsw,vegfra,rmol,ust,znt,dms_0, & erup_beg,erup_end, & - xland,xlat,xlong,z_at_w,z,smois,dustin,seasin, & + xland,xlat,xlong,mapfac_mx,mapfac_my,z_at_w,z,smois,dustin,seasin,& sebio_iso,sebio_oli,sebio_api,sebio_lim,sebio_xyl, & sebio_hc3,sebio_ete,sebio_olt,sebio_ket,sebio_ald, & sebio_hcho,sebio_eth,sebio_ora2,sebio_co,sebio_nr, & @@ -124,7 +129,7 @@ subroutine emissions_driver(id,ktau,dtstep,DX, & INTEGER, INTENT(IN ) :: & ktau,stepbioe,stepfirepl REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), & - INTENT(IN ) :: moist + INTENT(INOUT ) :: moist REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & INTENT(INOUT ) :: chem REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_tracer ), & @@ -204,7 +209,8 @@ subroutine emissions_driver(id,ktau,dtstep,DX, & rainnc, & potevp, & sfcevp, & - lu_index + lu_index, & + mapfac_mx,mapfac_my REAL, DIMENSION( ims:ime , jms:jme ) , & OPTIONAL, & @@ -312,7 +318,7 @@ subroutine emissions_driver(id,ktau,dtstep,DX, & real :: area,x1,percen_mass_umbrel,base_umbrel,ashz_above_vent REAL, INTENT(IN ) :: & - bioemdt, dtstep, dx, gmt, g + bioemdt, dtstep, dx, dy, gmt, g INTEGER, INTENT(IN ) :: & plumerisefire_frq @@ -396,7 +402,7 @@ subroutine emissions_driver(id,ktau,dtstep,DX, & REAL :: convert2(its:ite,jts:jte) #endif CHARACTER (LEN=80) :: message - LOGICAL :: do_bioemiss, do_plumerisefire,do_ex_volcanoe + LOGICAL :: do_bioemiss, do_plumerisefire INTEGER :: imod ! dust scheme option from namelist @@ -416,7 +422,11 @@ subroutine emissions_driver(id,ktau,dtstep,DX, & base_umbrel=.25 ! fraction ivolcano=0 - area=dx*dx + + !area=dx*dx + ! A. Ukhov 12 Oct 2024 + area=(dx/mapfac_mx(i,j))*(dy/mapfac_my(i,j)) + dust_emiss_active = 0 seasalt_emiss_active = 0 if(config_flags%dust_opt >= 2 )dust_emiss_active = 1 @@ -433,8 +443,8 @@ subroutine emissions_driver(id,ktau,dtstep,DX, & gmtm=mod(gmtp,60.) ! write(0,*) 'gmtp,gmtm,curr_secs = ',gmtp,gmtm,curr_secs ! - if(config_flags%emiss_opt_vol == 1 .or. config_flags%emiss_opt_vol == 2)then - do_ex_volcanoe = .false. + if(config_flags%emiss_opt_vol == 1 .or. config_flags%emiss_opt_vol == 2 .or. & + config_flags%emiss_opt_vol == 3 )then ! emiss_ash_height = config_flags%emiss_ash_hgt if(emiss_ash_height.gt. 1.)then @@ -448,7 +458,6 @@ subroutine emissions_driver(id,ktau,dtstep,DX, & so2_mass=1.5e4*3600.*1.e9/64./area eh=2600.*(emiss_ash_height*.0005)**4.1494 emiss_ash_mass=eh*1.e9/area - ! volcanic emissions ! ashz_above_vent=emiss_ash_height - z_at_w(i,kts,j) @@ -535,11 +544,12 @@ subroutine emissions_driver(id,ktau,dtstep,DX, & do j=jts,jte do i=its,ite ivolcano = 0 - if(erup_end(i,j).le.0)cycle + if(erup_end(i,j).le.0 .and. config_flags%emiss_opt_vol .ne. 3)cycle ! if(emis_vol(i,kts,j,p_e_vash1).le.0.)cycle ! ! erup_end is continuation in minutes ! + if(config_flags%emiss_opt_vol .ne. 3) then begday=int(erup_beg(i,j)/1000.)-1 beghr=int(erup_beg(i,j))-(begday+1)*1000 begmin=00. @@ -587,6 +597,13 @@ subroutine emissions_driver(id,ktau,dtstep,DX, & endif endif ! julday.ge.begday .and. julday.le.endday + ! A. Ukhov 12 Oct 2024. + elseif(config_flags%emiss_opt_vol .eq. 3) then + if (emis_vol(i,kts,j,p_e_start) .le.0) cycle + ivolcano=1 + !write(*,*)"volcano at ",ivolcano,i,j,emis_vol(i,kts,j,p_e_start),emis_vol(i,22,j,p_e_vso2),emis_vol(i,22,j,p_e_vash9) + endif !if(config_flags%emiss_opt_vol .ne. 3) + volc_select: SELECT CASE(config_flags%chem_opt) CASE (GOCART_SIMPLE,MOZCART_KPP,T1_MOZCART_KPP,GOCARTRADM2,GOCARTRACM_KPP) CALL wrf_debug(15,'Adding volcanic emissions') @@ -597,11 +614,11 @@ subroutine emissions_driver(id,ktau,dtstep,DX, & enddo do k=kts,kte conv=float(ivolcano)*alt(i,k,j)*dtstep/dz8w(i,k,j) - chem(i,k,j,p_p25)=chem(i,k,j,p_p25)+.5*emis_vol(i,k,j,p_e_vash10)*conv + chem(i,k,j,p_p25)=chem(i,k,j,p_p25)+ 0.672 * emis_vol(i,k,j,p_e_vash10)*conv !vash10:[0.0-3.906] diam. (um) => log(2.5-0)/log(3.906-0)=0.672 chem(i,k,j,p_p10)=chem(i,k,j,p_p10) & - +.5*emis_vol(i,k,j,p_e_vash10)*conv & - +emis_vol(i,k,j,p_e_vash9)*conv & - +.5*emis_vol(i,k,j,p_e_vash8)*conv + + emis_vol(i,k,j,p_e_vash10)*conv & + + emis_vol(i,k,j,p_e_vash9)*conv & + + 0.356 * emis_vol(i,k,j,p_e_vash8)*conv ! vash8:[7.812-15.625] diam. (um) => (log(10)-log(7.812))/(log(15.625)-log(7.812))=0.356 enddo CASE (RADM2SORG,RADM2SORG_AQ,RADM2SORG_KPP,RACMSORG_KPP,RACMSORG_AQ,RACM_ESRLSORG_KPP, & RACMSORG_AQCHEM_KPP,RACM_ESRLSORG_AQCHEM_KPP) @@ -630,9 +647,28 @@ subroutine emissions_driver(id,ktau,dtstep,DX, & CASE (CHEM_VOLC) CALL wrf_debug(15,'Adding volcanic emissions to case chem_volc') do k=kts,kte - conv = float(ivolcano)*4.828e-4/rho_phy(i,k,j)*dtstep/(dz8w(i,k,j)*60.) - chem(i,k,j,p_so2) = chem(i,k,j,p_so2) & - +emis_vol(i,k,j,p_e_vso2)*conv + ! p_e_vso2 in "ug/m2/min". p_so2 in ppmv + if (config_flags%emiss_opt_vol == 2) then + conv = float(ivolcano) *4.828e-4 * alt(i,k,j) * dtstep/(dz8w(i,k,j)*60.) + chem(i,k,j,p_so2) = chem(i,k,j,p_so2) + emis_vol(i,k,j,p_e_vso2) * conv + endif + + ! A. Ukhov. volcanic emissions + if (config_flags%emiss_opt_vol == 3) then + ! p_e_vso2 in "ug/m2/min". p_so2 in ppmv + conv = float(ivolcano) * 4.52656e-4 * alt(i,k,j) * dtstep/(dz8w(i,k,j)*60.) + chem(i,k,j,p_so2) = chem(i,k,j,p_so2) + emis_vol(i,k,j,p_e_vso2) * conv + + ! volcanic sulfate emission. + ! p_e_vsulf in "ug/m2/min". p_sulf in ppmv + conv = float(ivolcano) * 3.01771e-4 * alt(i,k,j) * dtstep/(dz8w(i,k,j)*60.) + chem(i,k,j,p_sulf) = chem(i,k,j,p_sulf) + emis_vol(i,k,j,p_e_vsulf) * conv + + ! volcanic water vapor emission. p_qv in [kg/kg] + ! p_e_qv in [kg/m2/sec] + conv = float(ivolcano) * alt(i,k,j) * dtstep / dz8w(i,k,j) + moist(i,k,j,p_qv) = moist(i,k,j,p_qv) + emis_vol(i,k,j,p_e_qv) * conv + endif enddo do k=kts,kte conv=float(ivolcano)*alt(i,k,j)*dtstep/dz8w(i,k,j) @@ -657,7 +693,7 @@ subroutine emissions_driver(id,ktau,dtstep,DX, & chem(i,k,j,p_vash_10)=chem(i,k,j,p_vash_10)+emis_vol(i,k,j,p_e_vash10)*conv enddo CASE (CHEM_VASH) - CALL wrf_debug(15,'Adding volcanic emissions to case chem_volc') + CALL wrf_debug(15,'Adding volcanic emissions to case chem_vash') do k=kts,kte conv=float(ivolcano)*alt(i,k,j)*dtstep/dz8w(i,k,j) chem(i,k,j,p_vash_1)=chem(i,k,j,p_vash_1)+emis_vol(i,k,j,p_e_vash1)*conv @@ -676,7 +712,7 @@ subroutine emissions_driver(id,ktau,dtstep,DX, & !!!!!! enddo enddo - ENDIF! config_flags%emiss_opt_vol == 1 .or. config_flags%emiss_opt_vol == 2 + ENDIF! config_flags%emiss_opt_vol == 1, 2 or 3 !-------------------------------------------------------------------------------------- do_plumerisefire = .false. IF ( config_flags%biomass_burn_opt == BIOMASSB_MOZC .OR. & diff --git a/chem/mechanism_driver.F b/chem/mechanism_driver.F index 9a68bfae28..7f2a7708e1 100755 --- a/chem/mechanism_driver.F +++ b/chem/mechanism_driver.F @@ -2,7 +2,7 @@ ! subroutine mechanism_driver(id,curr_secs,ktau,dtstep,ktauc,dtstepc,& config_flags,gmt,julday,t_phy,moist,p8w,t8w,gd_cldfr, & - p_phy,chem,rho_phy,dz8w,dx,g,z,z_at_w,xlat,xlong,vdrog3,& + p_phy,chem,rho_phy,dz8w,dx,dy,msftx,msfty,g,z,z_at_w,xlat,xlong,vdrog3,& vcsulf_old,vcso2_old,vch2o2_old,ttday,tcosz, & ph_macr,ph_o31d,ph_o33p,ph_no2, & ph_cl2,ph_hocl,ph_clno2,ph_fmcl, & @@ -13,7 +13,7 @@ subroutine mechanism_driver(id,curr_secs,ktau,dtstep,ktauc,dtstepc,& ph_n2o5,ph_o2,backg_oh,backg_h2o2,backg_no3, & addt,addx,addc,etep,oltp,olip,cslp,limp,hc5p,hc8p,tolp, & xylp,apip,isop,hc3p,ethp,o3p,tco3,mo2,o1d,olnn,rpho,xo2,& - ketp,olnd, & + ketp,olnd, volc_diags, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) @@ -24,6 +24,7 @@ subroutine mechanism_driver(id,curr_secs,ktau,dtstep,ktauc,dtstepc,& USE module_data_sorgam USE module_radm USE module_gocart_chem + USE module_volc_chem, only: gocart_volc_chem_driver ! A. Ukhov 24 Oct 2024 USE module_aerosols_sorgam USE module_cbmz, only: cbmz_driver IMPLICIT NONE @@ -34,7 +35,7 @@ subroutine mechanism_driver(id,curr_secs,ktau,dtstep,ktauc,dtstepc,& its,ite, jts,jte, kts,kte INTEGER, INTENT(IN ) :: ktau,ktauc REAL(KIND=8), INTENT(IN ) :: curr_secs - REAL, INTENT(IN ) :: dtstep,dtstepc,gmt,dx,g + REAL, INTENT(IN ) :: dtstep,dtstepc,gmt,dx,g,dy ! ! advected moisture variables ! @@ -47,7 +48,7 @@ subroutine mechanism_driver(id,curr_secs,ktau,dtstep,ktauc,dtstepc,& INTENT(INOUT ) :: chem REAL, DIMENSION( ims:ime, jms:jme ), & INTENT(IN ) :: & - xlat,xlong,ttday,tcosz + xlat,xlong,ttday,tcosz,msftx,msfty ! ! arrays that hold the photolysis rates @@ -70,7 +71,8 @@ subroutine mechanism_driver(id,curr_secs,ktau,dtstep,ktauc,dtstepc,& xylp,apip,isop,hc3p,ethp,o3p,tco3,mo2,o1d,olnn,rpho,xo2,& ketp,olnd - + ! A. Ukhov 24 Oct 2024, Volcanic diagnostics + REAL, DIMENSION( ims:ime , jms:jme ,1:9), INTENT(INOUT) :: volc_diags ! ! on input from meteorological model ! @@ -108,6 +110,19 @@ subroutine mechanism_driver(id,curr_secs,ktau,dtstep,ktauc,dtstepc,& ! select chemical mechanism ! chem_select: SELECT CASE(config_flags%chem_opt) + CASE (CHEM_VOLC) + !A. Ukhov 24 Oct 2024 + CALL wrf_debug(15,'calling volc. chem from mechanism_driver') + call gocart_volc_chem_driver(curr_secs,dtstepc,config_flags, & + gmt,julday,t_phy,moist, & + chem,rho_phy,dz8w,p8w,backg_oh,backg_h2o2, & + gd_cldfr,dx,dy,msftx,msfty,g, & + xlat,xlong,ttday,tcosz, & + volc_diags, & ! A. Ukhov Volc. diagnostics + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte ) + CASE (GOCART_SIMPLE,GOCARTRACM_KPP,GOCARTRADM2) CALL wrf_debug(15,'calling gocart chem from mechanism_driver') call gocart_chem_driver(curr_secs,dtstepc,config_flags, & diff --git a/chem/module_ctrans_grell.F b/chem/module_ctrans_grell.F index 1e4aa79055..1b82dd54e6 100755 --- a/chem/module_ctrans_grell.F +++ b/chem/module_ctrans_grell.F @@ -15,6 +15,7 @@ MODULE module_ctrans_grell USE module_state_description, ONLY: mozart_mosaic_4bin_kpp, & mozart_mosaic_4bin_aq_kpp, & mozcart_kpp, t1_mozcart_kpp, & + chem_vash,chem_volc,chem_volc_4bin, & ! A. Ukhov. p_hcho, p_c3h6ooh, p_onit, p_mvk, p_macr, & p_etooh, p_prooh, p_acetp, p_mgly, p_mvkooh, & p_onitr, p_isooh, p_ch3oh, p_c2h5oh, & @@ -56,7 +57,9 @@ MODULE module_ctrans_grell p_bsoaX_a04,p_bsoa1_a04,p_bsoa2_a04,p_bsoa3_a04,p_bsoa4_a04,& p_biog1_c_a04,p_biog1_o_a04,& p_cl_a04,p_co3_a04,p_nh4_a04,p_na_a04,& - p_ca_a04,p_oin_a04,p_oc_a04,p_bc_a04 + p_ca_a04,p_oin_a04,p_oc_a04,p_bc_a04,& + p_vash_1,p_vash_2,p_vash_3,p_vash_4,p_vash_5,p_vash_6,p_vash_7,& + p_vash_8,p_vash_9,p_vash_10 IMPLICIT NONE @@ -82,6 +85,7 @@ SUBROUTINE GRELLDRVCT(DT,itimestep,DX, & XLV,CP,G,r_v,z,cu_co_ten, & wd_no3,wd_so4, & wd_nh4,wd_oa, & + volc_diags, & ! A. Ukhov Volcanic diagnostics wd_so2, wd_sulf, wd_hno3, wd_nh3, & wd_cvasoa, wd_cvbsoa, wd_asoa, wd_bsoa, & k22_shallow,kbcon_shallow,ktop_shallow,xmb_shallow, & @@ -91,7 +95,7 @@ SUBROUTINE GRELLDRVCT(DT,itimestep,DX, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) ! USE module_configure -! USE module_state_description + USE module_state_description, only: p_wd_ash_cu !------------------------------------------------------------- IMPLICIT NONE !------------------------------------------------------------- @@ -140,6 +144,7 @@ SUBROUTINE GRELLDRVCT(DT,itimestep,DX, & ! Accumulated wet deposition ! REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: wd_no3,wd_so4 + REAL, DIMENSION( ims:ime , jms:jme ,1:9), INTENT(INOUT) :: volc_diags ! A. Ukhov REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: wd_nh4,wd_oa, & wd_so2, wd_sulf, wd_hno3, wd_nh3 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: & @@ -543,6 +548,32 @@ SUBROUTINE GRELLDRVCT(DT,itimestep,DX, & if (p_nh3 .gt.1) wdi_nh3(its:ite,jts:jte) = wdi_nh3(its:ite,jts:jte) + wetdep_2d(its:ite,jts:jte,p_nh3)*dt ! mmol/m2 + + ! A. Ukhov (KAUST) 18 Oct 2024 + ! Volcanic diagnostics for wet deposited ash + if (chemopt == chem_vash .OR. chemopt == chem_volc) then + ! wetdep_2d [ug/m2/s] for aerosol + volc_diags(its:ite,jts:jte,p_wd_ash_cu) = volc_diags(its:ite,jts:jte,p_wd_ash_cu) + & + (wetdep_2d(its:ite,jts:jte,p_vash_1)+ & + wetdep_2d(its:ite,jts:jte,p_vash_2) + & + wetdep_2d(its:ite,jts:jte,p_vash_3) + & + wetdep_2d(its:ite,jts:jte,p_vash_4) + & + wetdep_2d(its:ite,jts:jte,p_vash_5) + & + wetdep_2d(its:ite,jts:jte,p_vash_6) + & + wetdep_2d(its:ite,jts:jte,p_vash_7) + & + wetdep_2d(its:ite,jts:jte,p_vash_8) + & + wetdep_2d(its:ite,jts:jte,p_vash_9) + & + wetdep_2d(its:ite,jts:jte,p_vash_10))*dt ! ug/m2 + endif + + if(chemopt == chem_volc_4bin) then + volc_diags(its:ite,jts:jte,p_wd_ash_cu) = volc_diags(its:ite,jts:jte,p_wd_ash_cu) + & + (wetdep_2d(its:ite,jts:jte,p_vash_7) + & + wetdep_2d(its:ite,jts:jte,p_vash_8) + & + wetdep_2d(its:ite,jts:jte,p_vash_9) + & + wetdep_2d(its:ite,jts:jte,p_vash_10)) * dt ! ug/m2 + endif + ! Update the accumulated wet deposition: wd_no3(its:ite,jts:jte) = wd_no3(its:ite,jts:jte) + wdi_no3(its:ite,jts:jte) ! mmol/m2 diff --git a/chem/module_optical_averaging.F b/chem/module_optical_averaging.F index 278f1a65f5..72099272d3 100644 --- a/chem/module_optical_averaging.F +++ b/chem/module_optical_averaging.F @@ -301,6 +301,17 @@ subroutine optical_averaging(id,curr_secs,dtstep,config_flags, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) + + ! 28/10/2024 A. Ukhov. Volcanic case. + CASE (CHEM_VOLC) !todo: CHEM_VOLC_4BIN,CHEM_VASH + call optical_prep_volc(nbin_o, chem, alt,relhum, & + radius_core,radius_wet, number_bin, & + swrefindx,swrefindx_core, swrefindx_shell, & + lwrefindx,lwrefindx_core, lwrefindx_shell, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte ) + !gocart is now wavelength dependent --SAM 4/25/11 !and for shortwave and longwave - similar to modal --SAM 4/25/11 ! call optical_prep_gocart(nbin_o, chem, alt,relhum, & @@ -4185,6 +4196,374 @@ subroutine optical_prep_gocart(nbin_o, chem, alt,relhum, & end subroutine optical_prep_gocart + ! 28 Oct 2024. A. Ukhov (KAUST). calculation of volume-averaged refractive index + ! of volcanic ash and sulfate aerosol. Based on optical_prep_gocart + subroutine optical_prep_volc(nbin_o, chem, alt,relhum, & + radius_core,radius_wet, number_bin, & + swrefindx,swrefindx_core, swrefindx_shell, & + lwrefindx,lwrefindx_core, lwrefindx_shell, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte ) + ! + USE module_configure + USE module_model_constants + USE module_data_sorgam + + INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte, nbin_o + INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme + INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde + + REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & + INTENT(IN ) :: chem + REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & + INTENT(IN ) :: alt,relhum + REAL, DIMENSION( its:ite, kts:kte, jts:jte, 1:nbin_o), & + INTENT(OUT ) :: & + radius_wet, number_bin, radius_core + COMPLEX, DIMENSION( its:ite, kts:kte, jts:jte,1:nbin_o,nswbands), & + INTENT(OUT ) :: swrefindx, swrefindx_core, swrefindx_shell + COMPLEX, DIMENSION( its:ite, kts:kte, jts:jte,1:nbin_o,nlwbands), & + INTENT(OUT ) :: lwrefindx, lwrefindx_core, lwrefindx_shell + ! + ! local variables + ! + integer i, j, k, m, n, isize + complex ri_dum, ri_ave_a + + COMPLEX, DIMENSION(nswbands) :: & ! now only 5 aerosols have wave-dependent refr + swref_index_ash , swref_index_sulf, swref_index_h2o + + COMPLEX, DIMENSION(nlwbands) :: & ! now only 5 aerosols have wave-dependent refr + lwref_index_ash , lwref_index_sulf, lwref_index_h2o + + COMPLEX ref_index_bc,ref_index_oin + + real dens_so4 , dens_ash, dens_h2o + real mass_so4 , mass_ash, mass_so4i, mass_so4j, mass_h2o + + real vol_so4 , vol_ash, vol_h2o + real conv1a, conv1sulf + real mass_dry_a, mass_wet_a, vol_dry_a , vol_wet_a , vol_shell, & + dp_dry_a, dp_wet_a, num_a, dp_bc_a + real ifac, jfac + integer ns + real drydens,duma,dlo_um,dhi_um,sixpi,dtemp + integer iflag + real, save:: xnum_secti(8),xnum_sectj(8) + real, save:: xmas_secti(8),xmas_sectj(8) + real, dimension(1:nbin_o) :: xdia_um, xdia_cm + !REAL, PARAMETER :: FRAC2Aitken=0.25 ! Fraction of modal mass in Aitken mode + + !ukhova 23 March 2025. hardcoded value + REAL, PARAMETER :: sulf_FRAC2Aitken=0.05 ! Fraction of modal mass in Aitken mode - applied globally to each species + + !Ash bins. Diameters in microns + real*8, DIMENSION (10), PARAMETER :: da_ash(10)=(/1000d0,500d0, 250d0,125d0,62.5d0,31.25d0,15.625d0,7.8125d0,3.9065d0,0.039d0/) + real*8, DIMENSION (10), PARAMETER :: db_ash(10)=(/2000d0,1000d0,500d0,250d0,125d0, 62.5d0, 31.25d0,15.625d0,7.8125d0,3.900d0/) + real, parameter :: hygro_ash_aer = 0.1 + + + ! PARAMETERS OF TWO LOG-NORMAL DISTRIBUTIONS + ! mean diameter for volcanic sulf nuclei mode [ m ] + REAL dginin_sulf + PARAMETER (dginin_sulf=0.18E-6) + ! sigma-G for sulf nuclei mode + REAL sginin_sulf + PARAMETER (sginin_sulf=1.4) + + ! mean diameter for volcanic sulf accumulation mode [ m ] + REAL dginia_sulf + PARAMETER (dginia_sulf=0.64E-6) + ! sigma-G for sulf accumulation mode + REAL sginia_sulf + PARAMETER (sginia_sulf=1.6) + + !ASH RI + real,dimension(1:nswbands),save :: refrsw_ash,refisw_ash + real,dimension(1:nlwbands),save :: refrlw_ash,refilw_ash + !data refrsw_dust /nswbands*1.530/ + data refrsw_ash /nswbands*1.550/ + ! data refisw_ash /0.024,0.0135,0.0063,0.004/ + data refisw_ash /nswbands*0.001/ ! RI = 1.550 + i0.001 for ash. As in JGR, 2021 paper "How Does a Pinatubo-Size Volc..." + data refrlw_ash /2.340,2.904,1.748,1.508,1.911,1.822,2.917,1.557, & + 1.242,1.447,1.432,1.473,1.495,1.5,1.5,1.51/ + data refilw_ash /0.7,0.857,0.462,0.263,0.319,0.26,0.65,0.373,0.093, & + 0.105,0.061,0.0245,0.011,0.008,0.0068,0.018/ + + real dhi, dlo, xlo, xhi, dxbin, relh_frc + real dlo_sectm(nbin_o), dhi_sectm(nbin_o) + integer, parameter :: nbin_omoz=8 + real, save :: ashfrc_vol10bin(10,nbin_omoz) ! volcanic ash size distibution - mass fracs in MOSAIC 8-bins + real*8 dlogoc, dhigoc + integer, save :: kcall + data kcall / 0 / + + ! real sginin,sginia,sginic from module_data_sorgam.F + ! + ! Mass from modal distribution is divided into individual sections before + ! being passed back into the Mie routine. + ! use the same size bins as 8 default MOSAIC size bins + ! dlo_um and dhi_um define the lower and upper bounds of individual sections + ! used to compute optical properties + + sixpi=6.0/3.14159265359 + dlo_um=0.0390625 + dhi_um=10.0 + drydens=1.8 + iflag=2 + duma=1.0 + dtemp=dlo_um + + do isize=1,nbin_o + xdia_um(isize)=(dtemp+dtemp*2.0)/2.0 + dtemp=dtemp*2.0 + enddo + + if (kcall .eq. 0) then + ! calculate sectional contributions from ash + dlo = dlo_um*1.0e-6 + dhi = dhi_um*1.0e-6 + xlo = log( dlo ) + xhi = log( dhi ) + dxbin = (xhi - xlo)/nbin_o + + do n = 1, nbin_o + dlo_sectm(n) = exp( xlo + dxbin*(n-1) ) + dhi_sectm(n) = exp( xlo + dxbin*n ) + end do + + ashfrc_vol10bin=0. + do m = 1, 10 ! loop over ash size bins + dlogoc = da_ash(m)*1.E-6 ! low diameter limit (m) + dhigoc = db_ash(m)*1.E-6 ! hi diameter limit (m) + do n = 1, nbin_o + ashfrc_vol10bin(m,n)=max(DBLE(0.),min(DBLE(log(dhi_sectm(n))),log(dhigoc))- & + max(log(dlogoc),DBLE(log(dlo_sectm(n)))) )/(log(dhigoc)-log(dlogoc)) + !write(*,*)dlogoc,dhigoc,dlo_sectm(n),dhi_sectm(n),ashfrc_vol10bin(m,n) + end do + end do + kcall=kcall+1 + + !----- + WRITE(*,*)nbin_o + WRITE(*,*)'Ash redistribution:' + do m =1, 10 + WRITE(*,*)m,ashfrc_vol10bin(m,:) + end do + + ! Now divide mass into sections which is done by sect02: + ! * sect02 expects input in um + ! * pass in generic mass of 1.0 just to get a percentage distribution of mass among bins + + !aiken mode + call sect02(dginin_sulf*1.E6,sginin_sulf,drydens,iflag,duma,nbin_o,dlo_um,dhi_um, xnum_secti,xmas_secti) + + !accumulation mode + call sect02(dginia_sulf*1.E6,sginia_sulf,drydens,iflag,duma,nbin_o,dlo_um,dhi_um, xnum_sectj,xmas_sectj) + + write(*,*)'Sulfate Aitken mode: ' + write(*,*)xmas_secti + write(*,*)'Sulfate accumulation mode: ' + write(*,*)xmas_sectj + !----- + endif !kcall + + ! Define refractive indicies + do ns = 1, nswbands + swref_index_sulf(ns) = cmplx(refrsw_sulf(ns),refisw_sulf(ns)) + swref_index_ash(ns) = cmplx(refrsw_ash(ns),refisw_ash(ns)) + swref_index_h2o(ns) = cmplx(refrwsw(ns),refiwsw(ns)) + enddo + + do ns = 1, nlwbands + lwref_index_sulf(ns) = cmplx(refrlw_sulf(ns),refilw_sulf(ns)) + lwref_index_ash(ns) = cmplx(refrlw_ash(ns),refilw_ash(ns)) + lwref_index_h2o(ns) = cmplx(refrwlw(ns),refiwlw(ns)) + enddo + + ref_index_bc = cmplx(1.85,0.71) + ref_index_oin = cmplx(1.55,0.006) ! JCB, Feb. 20, 2008: "other inorganics" + + ! densities in g/cc + dens_so4 = 1.8 ! used + dens_ash = 2.5 ! used + dens_h2o = 1.0 + ! + swrefindx=0.0 + lwrefindx=0.0 + radius_wet=0.0 + number_bin=0.0 + radius_core=0.0 + swrefindx_core=0.0 + swrefindx_shell=0.0 + lwrefindx_core=0.0 + lwrefindx_shell=0.0 + ! + ! units: + ! * mass - g/cc(air) + ! * number - #/cc(air) + ! * volume - cc(air)/cc(air) + ! * diameter - cm + ! + do j = jts, jte + do k = kts, kte + do i = its, ite + mass_so4i = 0.0 + mass_so4j = 0.0 + mass_ash = 0.0 + + ! convert ug / kg dry air to g / cc air + conv1a = (1.0/alt(i,k,j)) * 1.0e-12 + ! convert ppmv sulfate to g / cc air + conv1sulf = (1.0/alt(i,k,j)) * 1.0e-9 * 96./28.97 + + ! Accumulation mode... + ! fraction of sulfate mass into modal accumulation mode + mass_so4j= (1.-sulf_FRAC2Aitken)*chem(i,k,j,p_sulf)*conv1sulf + + ! Aitken mode... + ! fraction of sulfate mass into modal Aitken mode + mass_so4i= sulf_FRAC2Aitken*chem(i,k,j,p_sulf)*conv1sulf + + !loop over 8 MOSAIC bins + do isize = 1, nbin_o + xdia_cm(isize)=xdia_um(isize)*1.0e-04 + mass_so4 = mass_so4i * xmas_secti(isize) + mass_so4j * xmas_sectj(isize) + + n = 0 + mass_ash = 0.0 + do m =p_vash_1, p_vash_10 ! loop over ash size bins less than 10 um diam + n = n+1 + mass_ash = mass_ash + ashfrc_vol10bin(n,isize) * chem(i,k,j,m) + end do + + mass_ash= mass_ash * conv1a + vol_so4 = mass_so4 / dens_so4 + vol_ash = mass_ash / dens_ash + + ! 7/23/09 SAM calculate vol_h2o from kappas in Petters and Kreidenweis ACP, 2007, vol. 7, 1961-1971. + ! Their kappas are the hygroscopicities used in Abdul-Razzak and Ghan, 2004, JGR, V105, p. 6837-6844. + ! These kappas are defined in module_data_sorgam and module_data_mosaic_asect. + ! Note that hygroscopicities are at 298K and specific surface tension - further refinement could + ! include temperature dependence in Petters and Kreidenweis + + relh_frc=amin1(.9,relhum(i,k,j)) ! Put in fractional relative humidity, max of .9, here + vol_h2o = vol_so4 * hygro_so4_aer + vol_ash * hygro_ash_aer + + vol_h2o = relh_frc*vol_h2o/(1.-relh_frc) + mass_h2o = vol_h2o*dens_h2o + mass_dry_a = mass_so4 + mass_ash + mass_wet_a = mass_dry_a + mass_h2o + + vol_dry_a = vol_so4 + vol_ash + vol_wet_a = vol_dry_a + vol_h2o + vol_shell = vol_wet_a + num_a = vol_wet_a / (0.52359877*xdia_cm(isize)*xdia_cm(isize)*xdia_cm(isize)) + + !SW + do ns=1,nswbands + ri_dum = (0.0,0.0) + ri_dum = (swref_index_sulf(ns) * mass_so4 / dens_so4) + & + (swref_index_ash(ns) * mass_ash / dens_ash) + & + (swref_index_h2o(ns) * mass_h2o / dens_h2o) + + ! if zero aerosols so add a check here to avoid divide by zero + IF(num_a .lt. 1.0e-20 .or. vol_wet_a .lt. 1.0e-20) then + dp_dry_a = xdia_cm(isize) + dp_wet_a = xdia_cm(isize) + dp_bc_a = 0.0 + ri_ave_a = 0.0 + ri_dum = 0.0 + else + dp_dry_a = (1.90985*vol_dry_a/num_a)**0.3333333 + dp_wet_a = (1.90985*vol_wet_a/num_a)**0.3333333 + dp_bc_a = 0.0 + ri_ave_a = ri_dum/vol_wet_a + ri_dum = (swref_index_sulf(ns) * mass_so4 / dens_so4) + & + (swref_index_ash(ns) * mass_ash / dens_ash) + & + (swref_index_h2o(ns) * mass_h2o / dens_h2o) + endif + + if(dp_wet_a/2.0 .lt. dlo_um*1.0e-4/2.0) then + swrefindx(i,k,j,isize,ns) = (1.5,0.0) + radius_wet(i,k,j,isize) = dlo_um*1.0e-4/2.0 + number_bin(i,k,j,isize) = num_a + radius_core(i,k,j,isize)= 0.0 + swrefindx_core(i,k,j,isize,ns) = ref_index_bc + swrefindx_shell(i,k,j,isize,ns) = ref_index_oin + elseif(num_a .lt. 1.e-20 .or. vol_shell .lt. 1.0e-20) then + swrefindx(i,k,j,isize,ns) = (1.5,0.0) + radius_wet(i,k,j,isize) =dlo_um*1.0e-4/2.0 + number_bin(i,k,j,isize) =num_a + radius_core(i,k,j,isize) =0.0 + swrefindx_core(i,k,j,isize,ns) = ref_index_bc + swrefindx_shell(i,k,j,isize,ns) = ref_index_oin + else + swrefindx(i,k,j,isize,ns) =ri_ave_a + radius_wet(i,k,j,isize) =dp_wet_a/2.0 + number_bin(i,k,j,isize) =num_a + radius_core(i,k,j,isize) =0.0 + swrefindx_core(i,k,j,isize,ns) =ref_index_bc + swrefindx_shell(i,k,j,isize,ns) =ri_dum/vol_shell + endif + enddo ! ns shortwave + + !LW + do ns=1,nlwbands + ri_dum = (0.0,0.0) + ri_dum = (lwref_index_sulf(ns) * mass_so4 / dens_so4) + & + (lwref_index_ash(ns) * mass_ash / dens_ash) + & + (lwref_index_h2o(ns) * mass_h2o / dens_h2o) + + ! if zero aerosols so add a check here to avoid divide by zero + IF(num_a .lt. 1.0e-20 .or. vol_wet_a .lt. 1.0e-20) then + dp_dry_a = xdia_cm(isize) + dp_wet_a = xdia_cm(isize) + ri_ave_a = 0.0 + ri_dum = 0.0 + else + dp_dry_a = (1.90985*vol_dry_a/num_a)**0.3333333 + dp_wet_a = (1.90985*vol_wet_a/num_a)**0.3333333 + ri_ave_a = ri_dum/vol_wet_a + ri_dum = (lwref_index_sulf(ns) * mass_so4 / dens_so4) + & + (lwref_index_ash(ns) * mass_ash / dens_ash) + & + (lwref_index_h2o(ns) * mass_h2o / dens_h2o) + endif + + if(dp_wet_a/2.0 .lt. dlo_um*1.0e-4/2.0) then + lwrefindx(i,k,j,isize,ns) = (1.5,0.0) + radius_wet(i,k,j,isize) = dlo_um*1.0e-4/2.0 + number_bin(i,k,j,isize) = num_a + radius_core(i,k,j,isize)= 0.0 + lwrefindx_core(i,k,j,isize,ns) = ref_index_bc + lwrefindx_shell(i,k,j,isize,ns) = ref_index_oin + elseif(num_a .lt. 1.e-20 .or. vol_shell .lt. 1.0e-20) then + lwrefindx(i,k,j,isize,ns) = (1.5,0.0) + radius_wet(i,k,j,isize) =dlo_um*1.0e-4/2.0 + number_bin(i,k,j,isize) =num_a + radius_core(i,k,j,isize) =0.0 + lwrefindx_core(i,k,j,isize,ns) = ref_index_bc + lwrefindx_shell(i,k,j,isize,ns) = ref_index_oin + else + lwrefindx(i,k,j,isize,ns) =ri_ave_a + radius_wet(i,k,j,isize) =dp_wet_a/2.0 + number_bin(i,k,j,isize) =num_a + radius_core(i,k,j,isize) =0.0 + lwrefindx_core(i,k,j,isize,ns) =ref_index_bc + lwrefindx_shell(i,k,j,isize,ns) =ri_dum/vol_shell + endif + enddo ! ns longwave + enddo !isize + + enddo !i + enddo !j + enddo !k + + return + + end subroutine optical_prep_volc + !below is the detail calculation for MIE code !czhao diff --git a/chem/module_vash_settling.F b/chem/module_vash_settling.F index 67205c68c6..65655f6100 100755 --- a/chem/module_vash_settling.F +++ b/chem/module_vash_settling.F @@ -1,17 +1,37 @@ +! 07/10/24 - A. Ukhov (KAUST), bug fix: ash mass balance was +! violated in the "settling" routine. I.e. prev. discrtetisation scheme +! did not conserve the mass of ash. Code clean-up. + +! 30/10/2024. A. UKhov (KAUST) Added dry deposition for volcanic ash, so2, sulf +! 30/10/2024. A. UKhov (KAUST) Added large scale scaveging for ash, so2, sulf + MODULE MODULE_VASH_SETTLING + USE module_state_description, only: p_vash_1,p_vash_2, & + p_vash_3,p_vash_4,p_vash_5,p_vash_6,p_vash_7, & + p_vash_8,p_vash_9,p_vash_10,p_sulf,p_so2, & + chem_vash,chem_volc,chem_volc_4bin,p_qc ! A. Ukhov. + +!Vertical grid indeces +! _ +!i+3 |_| l2-3 +!i+2 |_| l2-2 +!i+1 |_| l2-1 +! i |_| l2 +!i-1 |_| l2+1 +!i-2 |_| l2+2 +!i-3 |_| l2+3 +!//////////// CONTAINS SUBROUTINE vash_settling_driver(dt,config_flags,t_phy,moist, & - chem,rho_phy,dz8w,p8w,p_phy, & - ash_fall,dx,g, & + chem,rho_phy,dz8w,p8w,p_phy,volc_diags,dx,g, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) + USE module_configure USE module_state_description -! USE module_data_gocart_dust -! USE module_data_gocart_seas USE module_model_constants, ONLY: mwdry IMPLICIT NONE TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags @@ -26,12 +46,13 @@ SUBROUTINE vash_settling_driver(dt,config_flags,t_phy,moist, & INTENT(INOUT ) :: chem REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & INTENT(IN ) :: t_phy,p_phy,dz8w,p8w,rho_phy - REAL, DIMENSION( ims:ime , jms:jme ), & - INTENT(INOUT ) :: ash_fall + + ! A. Ukhov, Volcanic diagnostics + REAL, DIMENSION( ims:ime, jms:jme,1:9), INTENT(INOUT ) :: volc_diags REAL, INTENT(IN ) :: dt,dx,g - integer :: nmx,i,j,k,kk,lmx,iseas,idust - real*8, DIMENSION (1,1,kte-kts+1) :: tmp,airden,airmas,p_mid,delz,rh + integer :: nmx,i,j,k,kk,lmx + real*8, DIMENSION (1,1,kte-kts+1) :: tmp,airden,p_mid,delz,rh real*8, DIMENSION (1,1,kte-kts+1,4) :: sea_salt !srf real*8, DIMENSION (1,1,kte-kts+1,10) :: ash @@ -39,96 +60,94 @@ SUBROUTINE vash_settling_driver(dt,config_flags,t_phy,moist, & 2500.,2500.,2500.,2500.,2500. /) real*8, DIMENSION (10), PARAMETER :: reff_ash(10)=(/0.5000D-3,&! 1.00 mm diameter 0.3750D-3,&! 0.75 mm - 0.1875D-3,&! - 93.750D-6,&! - 46.875D-6,&! - 23.437D-6,&! - 11.719D-6,&! - 05.859D-6,&! - 02.930D-6,&! - 00.975D-6 /)! 3.9 um + 0.1875D-3,&! + 93.750D-6,&! + 46.875D-6,&! + 23.437D-6,&! + 11.719D-6,&! + 05.859D-6,&! + 02.930D-6,&! + 00.975D-6 /)! 3.9 um real*8, DIMENSION (10) :: bstl_ash - integer iash -!srf - -! -! bstl is for budgets -! + !real*8, DIMENSION (10,55) :: ash_speed - real*8 conver,converi - conver=1.e-9 - converi=1.e9 lmx=kte-kts+1 + do j=jts,jte do i=its,ite - kk=0 - bstl_ash(:)=0. - ash(:,:,:,:)=0. + kk=0 + + bstl_ash(:)=0. + ash(:,:,:,:)=0. + do k=kts,kte - kk=kk+1 - p_mid(1,1,kk)=.01*p_phy(i,kte-k+kts,j) - delz(1,1,kk)=dz8w(i,kte-k+kts,j) - airmas(1,1,kk)=-(p8w(i,k+1,j)-p8w(i,k,j))/g - airden(1,1,kk)=rho_phy(i,k,j) - tmp(1,1,kk)=t_phy(i,k,j) - rh(1,1,kk) = .95 - rh(1,1,kk) = MIN( .95, moist(i,k,j,p_qv) / & - (3.80*exp(17.27*(t_phy(i,k,j)-273.)/ & - (t_phy(i,k,j)-36.))/(.01*p_phy(i,k,j)))) - rh(1,1,kk)=max(1.0D-1,rh(1,1,kk)) + kk=kk+1 + p_mid(1,1,kk)=.01*p_phy(i,kte-k+kts,j) + delz(1,1,kk)=dz8w(i,kte-k+kts,j) + airden(1,1,kk)=rho_phy(i,k,j) + tmp(1,1,kk)=t_phy(i,k,j) + rh(1,1,kk) = .95 + rh(1,1,kk) = MIN( .95, moist(i,k,j,p_qv) / & + (3.80*exp(17.27*(t_phy(i,k,j)-273.)/ & + (t_phy(i,k,j)-36.))/(.01*p_phy(i,k,j)))) + rh(1,1,kk)=max(1.0D-1,rh(1,1,kk)) enddo -!ash settling - - iseas=0 - idust=0 - iash =1 - kk=0 -! write(0,*)'1',chem(i,1,j,p_dust_4) do k=kts,kte - kk=kk+1 - ash(1,1,kk,7)=chem(i,k,j,p_vash_7)*conver - ash(1,1,kk,8)=chem(i,k,j,p_vash_8)*conver - ash(1,1,kk,9)=chem(i,k,j,p_vash_9)*conver - ash(1,1,kk,10)=chem(i,k,j,p_vash_10)*conver + kk=kk+1 + ash(1,1,kk,7)=chem(i,k,j,p_vash_7) + ash(1,1,kk,8)=chem(i,k,j,p_vash_8) + ash(1,1,kk,9)=chem(i,k,j,p_vash_9) + ash(1,1,kk,10)=chem(i,k,j,p_vash_10) enddo - if(config_flags%chem_opt == 400 .or. config_flags%chem_opt == 402 ) then + + if(config_flags%chem_opt == CHEM_VASH .or. config_flags%chem_opt == CHEM_VOLC ) then kk=0 - do k=kts,kte - kk=kk+1 - ash(1,1,kk,1)=chem(i,k,j,p_vash_1)*conver - ash(1,1,kk,2)=chem(i,k,j,p_vash_2)*conver - ash(1,1,kk,3)=chem(i,k,j,p_vash_3)*conver - ash(1,1,kk,4)=chem(i,k,j,p_vash_4)*conver - ash(1,1,kk,5)=chem(i,k,j,p_vash_5)*conver - ash(1,1,kk,6)=chem(i,k,j,p_vash_6)*conver - enddo + do k=kts,kte + kk=kk+1 + ash(1,1,kk,1)=chem(i,k,j,p_vash_1) + ash(1,1,kk,2)=chem(i,k,j,p_vash_2) + ash(1,1,kk,3)=chem(i,k,j,p_vash_3) + ash(1,1,kk,4)=chem(i,k,j,p_vash_4) + ash(1,1,kk,5)=chem(i,k,j,p_vash_5) + ash(1,1,kk,6)=chem(i,k,j,p_vash_6) + enddo endif - call vsettling(1, 1, lmx, 10, g,& - ash, tmp, p_mid, delz, airmas, & - den_ash, reff_ash, dt, bstl_ash, rh, idust, iseas,iash) - kk=0 - ash_fall(i,j)=ash_fall(i,j)+sum(bstl_ash(1:10)) - do k=kts,kte - kk=kk+1 - chem(i,k,j,p_vash_7)=ash(1,1,kk,7)*converi - chem(i,k,j,p_vash_8)=ash(1,1,kk,8)*converi - chem(i,k,j,p_vash_9)=ash(1,1,kk,9)*converi - chem(i,k,j,p_vash_10)=ash(1,1,kk,10)*converi - enddo - if(config_flags%chem_opt == 400 .or. config_flags%chem_opt == 402 ) then + call vsettling(1, 1, lmx, 10, g, ash, tmp, p_mid, delz, & + den_ash, reff_ash, dt, bstl_ash, rh, airden)!,ash_speed) + + !if(i==20.and.j==20) then + ! write(*,*)'1',ash_speed(1,:) + ! write(*,*)'2',ash_speed(2,:) + ! ... + ! write(*,*)'10',ash_speed(10,:) + !endif + + !bstl_ash [ug/m2]. Ash fall in [kg/m2] + volc_diags(i,j,p_ash_fall) = volc_diags(i,j,p_ash_fall) + sum(1.e-9 * bstl_ash(1:10)) + kk=0 do k=kts,kte - kk=kk+1 - chem(i,k,j,p_vash_1)=ash(1,1,kk,1)*converi - chem(i,k,j,p_vash_2)=ash(1,1,kk,2)*converi - chem(i,k,j,p_vash_3)=ash(1,1,kk,3)*converi - chem(i,k,j,p_vash_4)=ash(1,1,kk,4)*converi - chem(i,k,j,p_vash_5)=ash(1,1,kk,5)*converi - chem(i,k,j,p_vash_6)=ash(1,1,kk,6)*converi + kk=kk+1 + chem(i,k,j,p_vash_7)=ash(1,1,kk,7) + chem(i,k,j,p_vash_8)=ash(1,1,kk,8) + chem(i,k,j,p_vash_9)=ash(1,1,kk,9) + chem(i,k,j,p_vash_10)=ash(1,1,kk,10) enddo + + if(config_flags%chem_opt == CHEM_VASH .or. config_flags%chem_opt == CHEM_VOLC ) then + kk=0 + do k=kts,kte + kk=kk+1 + chem(i,k,j,p_vash_1)=ash(1,1,kk,1) + chem(i,k,j,p_vash_2)=ash(1,1,kk,2) + chem(i,k,j,p_vash_3)=ash(1,1,kk,3) + chem(i,k,j,p_vash_4)=ash(1,1,kk,4) + chem(i,k,j,p_vash_5)=ash(1,1,kk,5) + chem(i,k,j,p_vash_6)=ash(1,1,kk,6) + enddo endif !ash settling end @@ -138,74 +157,51 @@ SUBROUTINE vash_settling_driver(dt,config_flags,t_phy,moist, & END SUBROUTINE vash_settling_driver - subroutine vsettling(imx,jmx, lmx, nmx,g0, & - tc, tmp, p_mid, delz, airmas, & - den, reff, dt, bstl, rh, idust, iseas,iash) -! **************************************************************************** -! * * -! * Calculate the loss by settling, using an implicit method * -! * * -! * Input variables: * -! * SIGE(k) - sigma coordinate of the vertical edges * -! * PS(i,j) - Surface pressure (mb) * -! * TMP(i,j,k) - Air temperature (K) * -! * CT(i,j) - Surface exchange coeff for moisture -! * * -! **************************************************************************** - +subroutine vsettling(imx,jmx, lmx, nmx, g0, & + tc, tmp, p_mid, delz, den, reff, dt, bstl, rh,airden)!,ash_speed) IMPLICIT NONE - INTEGER, INTENT(IN) :: imx, jmx, lmx, nmx,iseas,idust,iash + INTEGER, INTENT(IN) :: imx, jmx, lmx, nmx INTEGER :: ntdt - REAL, INTENT(IN) :: dt,g0 ! ,dyn_visc - REAL*8, INTENT(IN) :: tmp(imx,jmx,lmx), delz(imx,jmx,lmx), & - airmas(imx,jmx,lmx), rh(imx,jmx,lmx), & - den(nmx), reff(nmx), p_mid(imx,jmx,lmx) + REAL, INTENT(IN) :: dt,g0 + REAL*8, INTENT(IN) :: tmp(imx,jmx,lmx), delz(imx,jmx,lmx), & + rh(imx,jmx,lmx), den(nmx), reff(nmx), & + p_mid(imx,jmx,lmx),airden(imx,jmx,lmx) REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx) REAL*8, INTENT(OUT) :: bstl(imx,jmx,nmx) + !real*8, INTENT(INOUT) :: ash_speed(10,55) - REAL*8 :: tc1(imx,jmx,lmx,nmx), dt_settl(nmx), rcm(nmx), rho(nmx) + REAL*8 :: dt_settl(nmx), rho(nmx) INTEGER :: ndt_settl(nmx) - REAL*8 :: dzmin, vsettl, dtmax, pres, rhb, rwet(nmx), ratio_r(nmx) - REAL*8 :: addmass,c_stokes, free_path, c_cun, viscosity, vd_cor, growth_fac + REAL*8 :: vd_cor(lmx) + REAL*8 :: dzmin, vsettl, dtmax, rwet(nmx) + REAL*8 :: addmass,c_stokes, free_path, c_cun, viscosity,R_tilde REAL, PARAMETER :: dyn_visc = 1.5E-5 INTEGER :: k, n, i, j, l, l2 - ! for sea-salt: - REAL*8, PARAMETER :: c1=0.7674, c2=3.079, c3=2.573E-11, c4=-1.424 ! for OMP: REAL*8 :: rwet_priv(nmx), rho_priv(nmx) ! executable statements -! IF (type) /= 'dust' .AND. TRIM(aero_type) /= 'sea_salt') RETURN - if(idust.ne.1.and.iseas.ne.1.and.iash.ne.1)return - WHERE (tc(:,:,:,:) < 0.0) tc(:,:,:,:) = 1.0d-32 dzmin = MINVAL(delz(:,:,:)) - IF (idust == 1) growth_fac = 1.0 - IF (iseas == 1) growth_fac = 3.0 - IF (iash == 1) growth_fac = 1.0 - - DO k = 1,nmx + + DO k = 1,nmx ! number of bins ! Settling velocity (m/s) for each tracer (Stokes Law) ! DEN density (kg/m3) ! REFF effective radius (m) ! dyn_visc dynamic viscosity (kg/m/s) ! g0 gravity (m/s2) - ! 3.0 corresponds to a growth of a factor 3 of radius with 100% RH - ! 0.5 upper limit with temp correction - tc1(:,:,:,k) = tc(:,:,:,k) - vsettl = 2.0/9.0 * g0 * den(k) * (growth_fac*reff(k))**2 / & - (0.5*dyn_visc) + vsettl = 4.0/9.0 * g0 * den(k) * reff(k)**2 / dyn_visc ! Determine the maximum time-step satisying the CFL condition: ! dt <= (dz)_min / v_settl - ntdt=INT(dt) + ntdt = INT(dt) dtmax = dzmin / vsettl ndt_settl(k) = MAX( 1, INT( ntdt /dtmax) ) ! limit maximum number of iterations @@ -213,99 +209,748 @@ subroutine vsettling(imx,jmx, lmx, nmx,g0, & dt_settl(k) = REAL(ntdt) / REAL(ndt_settl(k)) ! Particles radius in centimeters - IF (iseas.eq.1)rcm(k) = reff(k)*100.0 -!srf IF (idust.eq.1)then - IF (idust.eq.1 .or. iash==1)then - rwet(k) = reff(k) - ratio_r(k) = 1.0 - rho(k) = den(k) - endif + rwet(k) = reff(k) + rho(k) = den(k) END DO - ! Solve the bidiagonal matrix (l,l) - !$OMP PARALLEL DO & !$OMP DEFAULT( SHARED ) & -!$OMP PRIVATE( i, j, l, l2, n, k, rhb, rwet_priv, ratio_r, c_stokes)& +!$OMP PRIVATE( i, j, l, l2, n, k, rwet_priv, c_stokes)& !$OMP PRIVATE( free_path, c_cun, viscosity, rho_priv, vd_cor ) - ! Loop over latitudes - DO j = 1,jmx - - DO k = 1,nmx - IF (idust.eq.1 .or. iash==1) THEN - rwet_priv(k) = rwet(k) - rho_priv(k) = rho(k) - END IF + + DO j = 1,jmx ! lat loop + DO i = 1,imx ! lon loop + DO k = 1,nmx ! bin loop + rwet_priv(k) = rwet(k) + rho_priv(k) = rho(k) - DO n = 1,ndt_settl(k) + bstl(i,j,k)=0. - ! Solve each vertical layer successively (layer l) - - DO l = lmx,1,-1 + DO n = 1,ndt_settl(k) + DO l = lmx,1,-1 ! height loop, from top l2 = lmx - l + 1 -! DO j = 1,jmx - DO i = 1,imx - - ! Dynamic viscosity - c_stokes = 1.458E-6 * tmp(i,j,l)**1.5/(tmp(i,j,l) + 110.4) - - ! Mean free path as a function of pressure (mb) and - ! temperature (K) - ! order of p_mid is top->sfc - free_path = 1.1E-3/p_mid(i,j,l2)/SQRT(tmp(i,j,l)) -!!! free_path = 1.1E-3/p_edge(i,j,l2)/SQRT(tmp(i,j,l)) - - ! Slip Correction Factor - c_cun = 1.0+ free_path/rwet_priv(k)* & - (1.257 + 0.4*EXP(-1.1*rwet_priv(k)/free_path)) - - ! Corrected dynamic viscosity (kg/m/s) - viscosity = c_stokes / c_cun - - ! Settling velocity -! IF (iseas.eq.1) THEN -! rho_priv(k) = ratio_r(k)*den(k) + (1.0 - ratio_r(k))*1000.0 -! END IF - - vd_cor = 2.0/9.0*g0*rho_priv(k)*rwet_priv(k)**2/viscosity - - ! Update mixing ratio - ! Order of delz is top->sfc - IF (l == lmx) THEN - tc(i,j,l,k) = tc(i,j,l,k) / & - (1.0 + dt_settl(k)*vd_cor/delz(i,j,l2)) - ELSE - tc(i,j,l,k) = 1.0/(1.0+dt_settl(k)*vd_cor/delz(i,j,l2))& - *(tc(i,j,l,k) + dt_settl(k)*vd_cor /delz(i,j,l2-1) & - * tc(i,j,l+1,k)) - END IF - END DO !i -! END DO !j - END DO !l - - END DO !n - END DO !k - + ! Dynamic viscosity + c_stokes = 1.458E-6 * tmp(i,j,l)**1.5/(tmp(i,j,l) + 110.4) + + ! Mean free path as a function of pressure (mb) and + ! temperature (K) + ! order of p_mid is top->sfc + free_path = 1.1E-3/p_mid(i,j,l2)/SQRT(tmp(i,j,l)) + + ! Slip Correction Factor + c_cun = 1.0+ free_path/rwet_priv(k)*(1.257 + 0.4*EXP(-1.1*rwet_priv(k)/free_path)) + + ! Corrected dynamic viscosity (kg/m/s) + viscosity = c_stokes / c_cun + + ! Settling velocity, depends on temp + vd_cor(l) = 2.0/9.0*g0*rho_priv(k)*rwet_priv(k)**2/viscosity + + ! Correction of ash settling speed for partciles with radii>20 micron + ! details in Mailler, Sylvain, et al. 2023 + if (k < 7) then + R_tilde = airden(i,j,l) * (2 * rwet_priv(k)) * vd_cor(l) / (2.*viscosity) + if(R_tilde>0.0116)then + vd_cor(l) = vd_cor(l) * (1 - (1 + (R_tilde /2.440)**(-0.4335))**(-1.905)) + endif + END IF + + ! Update mixing ratio. Order of delz is top->sfc + IF (l == lmx) THEN + tc(i,j,l,k) = tc(i,j,l,k) / (1.0 + dt_settl(k) * vd_cor(l)/delz(i,j,l2)) + ELSE + if (l==1) then + bstl(i,j,k) = bstl(i,j,k) + (tc(i,j,l,k) * dt_settl(k) * vd_cor(l)/delz(i,j,l2)) * & + airden(i,j,l) * delz(i,j,l2) !ug/m2 + endif + tc(i,j,l,k) = 1.0/(1.0 + dt_settl(k) * vd_cor(l)/delz(i,j,l2)) * & + (tc(i,j,l,k) + tc(i,j,l+1,k) * dt_settl(k) * vd_cor(l+1) /delz(i,j,l2-1) * & + ((airden(i,j,l+1)*delz(i,j,l2-1))/(airden(i,j,l)*delz(i,j,l2)))) + END IF + + !ash_speed(k,l)=vd_cor(l) + !write(*,*)'SA ',i,j,l2,k,vd_cor(l)!,delz(i,j,l2) + + END DO !l, height + END DO !n, time + END DO !k, bin + + END DO !i END DO !j !$OMP END PARALLEL DO - DO n = 1,nmx - DO i = 1,imx - DO j = 1,jmx - bstl(i,j,n) = 0.0 - addmass=0. - DO l = 1,lmx - addmass=addmass+(tc(i,j,l,n) - tc1(i,j,l,n)) * airmas(i,j,l) - IF (tc(i,j,l,n) < 0.0) tc(i,j,l,n) = 1.0D-32 - END DO - if(addmass.gt.0.)addmass=0 - bstl(i,j,n) = bstl(i,j,n) - addmass - END DO - END DO - END DO - END SUBROUTINE vsettling +! A. Ukhov 30 October. +! Volcanic ash, so2, sulf dry deposition. Based on gocart_drydep_driver() +subroutine volc_ash_sulf_so2_drydep_driver(dtstep, & + config_flags,numgas, & + t_phy,moist,p8w,t8w,rmol,aer_res, & + p_phy,chem,rho_phy,dz8w,ddvel,xland,hfx, & + tsk,pbl,ust,znt,volc_diags, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte ) + +USE module_model_constants +USE module_configure +USE module_state_description + +IMPLICIT NONE + +INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte,numgas +!INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN) :: ivgtyp + +REAL, INTENT(IN ) :: dtstep +REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), & +INTENT(IN ) :: moist +REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & +INTENT(INOUT ) :: chem +REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , & +INTENT(IN ) :: t_phy, p_phy,dz8w,t8w,p8w,rho_phy +REAL, DIMENSION( its:ite, jts:jte, num_chem ), & +INTENT(INOUT) :: ddvel +REAL, DIMENSION( ims:ime , jms:jme ) , & +INTENT(IN) :: tsk,pbl,ust,rmol,xland,znt,hfx +REAL, DIMENSION( its:ite, jts:jte ), INTENT(IN) :: aer_res +REAL, DIMENSION( ims:ime, jms:jme,1:9), INTENT(INOUT ) :: volc_diags + +!! .. Local Scalars .. +INTEGER :: n, nr, ipr, imx,jmx,lmx + +integer :: ii,jj,kk,i,j,k,nv +integer,dimension (1,1) :: ilwi,ireg + +REAL :: rad, ta, z1 +real*8, dimension (1,1) :: z0,airden,delz_sfc,hflux,ts,pblz,ustar,ps +REAL*8 :: dvel(1,1) +LOGICAL :: highnh3, rainflag, vegflag, wetflag + +TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags + + do nv=numgas+1,num_chem + do j=jts,jte + do i=its,ite + ddvel(i,j,nv) = 0. + enddo + enddo + enddo + + imx = 1 + jmx = 1 + lmx = 1 + do j=jts,jte + do i=its,ite + ipr = 0 + dvel(1,1) = 0._8 + + if( xland(i,j) > 1.5 ) then + ilwi(1,1) = 1 + else + ilwi(1,1) = 0 + endif + ii = 1 + ireg(1,1) = 1 + airden(1,1) = real( rho_phy(i,kts,j),kind=8 ) + delz_sfc(1,1) = real( dz8w(i,kts,j),kind=8 ) + ustar(1,1) = max( 1.e-1_8,real(ust(i,j),kind=8) ) + hflux(1,1) = real( hfx(i,j),kind=8 ) + pblz(1,1) = real( pbl(i,j),kind=8 ) + ps(1,1) = real(p8w(i,kts,j),kind=8)*.01_8 + z0(1,1) = real( znt(i,j),kind=8 ) + ts(1,1) = real( tsk(i,j),kind=8 ) + + call dry_depvel_ash_sulf_so2(config_flags,ipr,ii,imx,jmx,lmx,& + airden, delz_sfc, pblz, ts, ustar, hflux, ilwi, & + ps, z0, dvel,g,rmol(i,j),aer_res(i,j)) + + if(config_flags%chem_opt == chem_vash .or. config_flags%chem_opt == chem_volc ) then + do nv = p_vash_1,p_vash_10 + ddvel(i,j,nv) = real( dvel(1,1),kind=4 ) + enddo + endif + + if(config_flags%chem_opt == chem_volc_4bin) then + do nv = p_vash_7,p_vash_10 + ddvel(i,j,nv) = real( dvel(1,1),kind=4 ) + enddo + endif + + if(config_flags%chem_opt == chem_volc ) then + ddvel(i,j,p_sulf) = real( dvel(1,1),kind=4 ) + !SO2 + volc_diags(i,j,p_so2_drydep)= volc_diags(i,j,p_so2_drydep) + chem(i,1,j,p_so2) * airden(1,1) * ddvel(i,j,p_so2) * dtstep * 1.E-6 / 28.97 ![mol m^-2] + !sulf + volc_diags(i,j,p_sulf_drydep)= volc_diags(i,j,p_sulf_drydep) + chem(i,1,j,p_sulf) * airden(1,1) * ddvel(i,j,p_sulf) * dtstep * 1.E-6 / 28.97 ![mol m^-2] + endif + + !drydep [kg/m2] = drydep [kg/m2]+1.e-9*dt[s]*dvel [m/s] * airden [kg/m3] * chem [ug/kg] + !ashdrydep(i,j) = ashdrydep(i,j) + 1.e-9 * dtstep * dvel(1,1) * airden(1,1) * chem(i,1,j,p_vash_10) + + volc_diags(i,j,p_ash_drydep) = volc_diags(i,j,p_ash_drydep) + 1.e-9 * dtstep * dvel(1,1) * airden(1,1) * & + (chem(i,1,j,p_vash_7) + chem(i,1,j,p_vash_8) + chem(i,1,j,p_vash_9) + chem(i,1,j,p_vash_10)) + + if(config_flags%chem_opt == chem_vash .or. config_flags%chem_opt == chem_volc ) then + volc_diags(i,j,p_ash_drydep) = volc_diags(i,j,p_ash_drydep) + 1.e-9 * dtstep * dvel(1,1) * airden(1,1) * & + (chem(i,1,j,p_vash_1) + chem(i,1,j,p_vash_2) + chem(i,1,j,p_vash_3) + chem(i,1,j,p_vash_4)+ & + chem(i,1,j,p_vash_5) + chem(i,1,j,p_vash_6)) + + endif + + enddo + enddo + +end subroutine volc_ash_sulf_so2_drydep_driver + +SUBROUTINE dry_depvel_ash_sulf_so2( config_flags,ipr,ii,imx,jmx,lmx, & + airden, delz_sfc, pblz, ts, ustar, hflux, ilwi, & + ps, z0, dvel,g0,rmol,aer_res) + + ! **************************************************************************** + ! * * + ! * Calculate dry deposition velocity. * + ! * * + ! * Input variables: * + ! * AEROSOL(k) - Logical, T = aerosol species, F = gas species * + ! * IREG(i,j) - # of landtypes in grid square * + ! * ILAND(i,j,ldt) - Land type ID for element ldt =1,IREG(i,j) * + ! * IUSE(i,j,ldt) - Fraction of gridbox area occupied by land type * + ! * element ldt * + ! * USTAR(i,j) - Friction velocity (m s-1) * + ! * DELZ_SFC(i,j) - Thickness of layer above surface * + ! * PBLZ(i,j) - Mixing depth (m) * + ! * Z0(i,j) - Roughness height (m) * + ! * * + ! * Determined in this subroutine (local): * + ! * OBK - Monin-Obukhov length (m): set to 1.E5 m under * + ! * neutral conditions * + ! * Rs(ldt) - Bulk surface resistance(s m-1) for species k to * + ! * surface ldt * + ! * Ra - Aerodynamic resistance. * + ! * Rb - Sublayer resistance. * + ! * Rs - Surface resistance. * + ! * Rttl - Total deposition resistance (s m-1) for species k * + ! * Rttl(k) = Ra + Rb + Rs. * + ! * * + ! * Returned: * + ! * DVEL(i,j,k) - Deposition velocity (m s-1) of species k * + ! * DRYDf(i,j,k) - Deposition frequency (s-1) of species k, * + ! * = DVEL / DELZ_SFC * + ! * * + ! **************************************************************************** + + USE module_configure + + IMPLICIT NONE + + REAL*8, INTENT(IN) :: airden(imx,jmx), delz_sfc(imx,jmx) + REAL*8, INTENT(IN) :: hflux(imx,jmx), ts(imx,jmx) + REAL*8, INTENT(IN) :: ustar(imx,jmx), pblz(imx,jmx) + REAL*8, INTENT(IN) :: ps(imx,jmx) + INTEGER, INTENT(IN) :: ilwi(imx,jmx) + INTEGER, INTENT(IN) :: imx,jmx,lmx + REAL*8, INTENT(IN) :: z0(imx,jmx) + REAL, INTENT(IN) :: g0,rmol,aer_res + REAL*8, INTENT(OUT) :: dvel(imx,jmx)!, drydf(imx,jmx) + + TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags + + REAL*8 :: obk, vds, czh, rttl, frac, logmfrac, psi_h, cz, eps + REAL*8 :: vd, ra, rb, rs + INTEGER :: ipr,i, j, k, ldt, iolson, ii + CHARACTER(LEN=50) :: msg + REAL*8 :: prss, tempk, tempc, xnu, ckustr, reyno, aird, diam, xm, z + REAL*8 :: frpath, speed, dg, dw, rt + REAL*8 :: rad0, rix, gfact, gfaci, rdc, rixx, rluxx, rgsx, rclx + REAL*8 :: dtmp1, dtmp2, dtmp3, dtmp4 + REAL*8 :: biofit,vk + + ! executable statements + j_loop: DO j = 1,jmx + i_loop: DO i = 1,imx + vk = .4_8 + vd = 0._8 + ra = 0._8 + rb = 0._8 ! only required for gases (SO2) + rs = 0.0_8 + + ! **************************************************************************** + ! * Compute the the Monin-Obhukov length. * + ! * The direct computation of the Monin-Obhukov length is: * + ! * * + ! * - Air density * Cp * T(surface air) * Ustar^3 * + ! * OBK = ---------------------------------------------- * + ! * vK * g * Sensible Heat flux * + ! * * + ! * Cp = 1000 J/kg/K = specific heat at constant pressure * + ! * vK = 0.4 = von Karman's constant * + ! **************************************************************************** + + IF (abs(hflux(i,j)) <= 1.e-5_8) THEN + obk = 1.0E5_8 + ELSE + ! MINVAL(hflux), MINVAL(airden), MINVAL(ustar) =?? + obk = -airden(i,j) * 1000.0_8 * ts(i,j) * (ustar(i,j))**3 & + / (vk * real(g0,kind=8) * hflux(i,j)) + ! -- debug: + IF ( obk == 0.0_8 ) WRITE(*,211) obk, i, j + 211 FORMAT(1X,'OBK=', E11.2, 1X,' i,j = ', 2I4) + END IF + ! write(0,*)1./obk,rmol + + if(rmol.ne.0.)then + obk=1._8/real(rmol,kind=8) + else + obk=1.e5_8 + endif + + ! cz = delz_sfc(i,j) / 2.0_8 ! center of the grid box above surface + cz = 2._8 + + ! **************************************************************************** + ! * (1) Aerosodynamic resistance Ra and sublayer resistance Rb. * + ! * * + ! * The Reynolds number REYNO diagnoses whether a surface is * + ! * aerodynamically rough (REYNO > 10) or smooth. Surface is * + ! * rough in all cases except over water with low wind speeds. * + ! * * + ! * For gas species over land and ice (REYNO >= 10) and for aerosol * + ! * species for all surfaces: * + ! * * + ! * Ra = 1./VT (VT from GEOS Kzz at L=1, m/s). * + ! * * + ! * The following equations are from Walcek et al, 1986: * + ! * * + ! * For gas species when REYNO < 10 (smooth), Ra and Rb are combined * + ! * as Ra: * + ! * * + ! * Ra = { ln(ku* z1/Dg) - Sh } / ku* eq.(13) * + ! * * + ! * where z1 is the altitude at the center of the lowest model layer * + ! * (CZ); * + ! * Sh is a stability correction function; * + ! * k is the von Karman constant (0.4, vK); * + ! * u* is the friction velocity (USTAR). * + ! * * + ! * Sh is computed as a function of z1 and L eq ( 4) and (5)): * + ! * * + ! * 0 < z1/L <= 1: Sh = -5 * z1/L * + ! * z1/L < 0: Sh = exp{ 0.598 + 0.39*ln(E) - 0.09(ln(E))^2 } * + ! * where E = min(1,-z1/L) (Balkanski, thesis). * + ! * * + ! * For gas species when REYNO >= 10, * + ! * * + ! * Rb = 2/ku* (Dair/Dg)**(2/3) eq.(12) * + ! * where Dg is the gas diffusivity, and * + ! * Dair is the air diffusivity. * + ! * * + ! * For aerosol species, Rb is combined with surface resistance as Rs. * + ! * * + ! **************************************************************************** + + frac = cz / obk + IF (frac > 1.0_8) frac = 1.0_8 + IF (frac > 0.0_8 .AND. frac <= 1.0_8) THEN + psi_h = -5.0_8*frac + ELSE IF (frac < 0.0_8) THEN + eps = MIN(1.0_8, -frac) + logmfrac = LOG(eps) + psi_h = EXP( 0.598_8 + 0.39_8 * logmfrac - 0.09_8 * (logmfrac)**2 ) + END IF + !-------------------------------------------------------------- + ! Aerosol species, Rs here is the combination of Rb and Rs. + + ra = (LOG(cz/z0(i,j)) - psi_h) / (vk*ustar(i,j)) + + vds = 0.002_8*ustar(i,j) + IF (obk < 0.0_8) vds = vds * (1.0_8+(-300.0_8/obk)**0.6667_8) + + czh = pblz(i,j)/obk + IF (czh < -30.0_8) vds = 0.0009_8*ustar(i,j)*(-czh)**0.6667_8 + + IF( config_flags%chem_opt /= CHEM_VASH .and. & + config_flags%chem_opt /= chem_volc_4bin )THEN + ra = real(aer_res,kind=8) + ENDIF + + ! --Set Vds to be less than VDSMAX (entry in input file divided -- + ! by 1.E4). VDSMAX is taken from Table 2 of Walcek et al. [1986]. + ! Invert to get corresponding R + ! if(ii.eq.1) then + ! rs=1.0_8/MIN(vds,2.0e-2_8) + ! else + rs=1.0_8/MIN(vds,2.0e-3_8) + ! endif + + ! ------ Set max and min values for bulk surface resistances ------ + + rs= MAX(1.0_8, MIN(rs, 9.9990e+3_8)) + + ! **************************************************************************** + ! * * + ! * Compute dry deposition velocity. * + ! * * + ! * IUSE is the fraction of the grid square occupied by surface ldt in * + ! * units of per mil (IUSE=500 -> 50% of the grid square). Add the * + ! * contribution of surface type ldt to the deposition velocity; this is * + ! * a loop over all surface types in the gridbox. * + ! * * + ! * Total resistance = Ra + Rb + Rs. + ! * * + ! **************************************************************************** + + rttl = ra + rb + rs + vd = vd + 1._8/rttl + + ! ------ Load array DVEL ------ + !if(ipr.eq.1) write(0,*)rs,ra,rb,vd + + dvel(i,j) = vd !* 1.2 + + ! -- Set a minimum value for DVEL + ! MIN(VdSO2) = 2.0e-3 m/s over ice + ! = 3.0e-3 m/s over land + ! MIN(vd_aerosol) = 1.0e-4 m/s + + IF (dvel(i,j) < 1.0E-4_8) dvel(i,j) = 1.0E-4_8 + + !drydf(i,j) = dvel(i,j) / delz_sfc(i,j) + + END DO i_loop + END DO j_loop + +END SUBROUTINE dry_depvel_ash_sulf_so2 + + +! A. Ukhov 6 Nov 2024 +! Gravitational settling for volcanic sulfate aerosol (accumulation mode only) +! Based on gocart_settling_driver and on +! Ukhov, A., et al, Improving dust simulations in WRF-Chem v4.1.3 coupled..., +! Geosci. Model Dev., 14, 473–493, https://doi.org/10.5194/gmd-14-473-2021, 2021. +SUBROUTINE sulf_settling_driver(dt,config_flags,t_phy,moist, & + chem,rho_phy,dz8w,p8w,p_phy,dx,g,z_at_w,volc_diags, & + ids,ide,jds,jde,kds,kde, & + ims,ime,jms,jme,kms,kme, & + its,ite,jts,jte,kts,kte) + +USE module_configure +USE module_state_description +USE module_data_sorgam, only: rhoso4 +IMPLICIT NONE + +TYPE(grid_config_rec_type), INTENT(IN) :: config_flags + +INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte + +REAL, DIMENSION(ims:ime,kms:kme,jms:jme,num_moist), INTENT(IN) :: moist +REAL, DIMENSION(ims:ime,kms:kme,jms:jme,num_chem), INTENT(INOUT) :: chem +REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN) :: t_phy,p_phy,dz8w,p8w,rho_phy,z_at_w +REAL, DIMENSION(ims:ime,jms:jme,1:9), INTENT(INOUT) :: volc_diags + +REAL, INTENT(IN) :: dt,dx,g +INTEGER :: i,j,k,kk,lmx +REAL*8, DIMENSION (1,1,kte-kts+1) :: tmp,airden,p_mid,delz,rh +REAL*8, DIMENSION (1,1,kte-kts+1) :: dsulf +real*8 :: bstl_sulf +REAL reff_sulf +REAL, PARAMETER :: dyn_visc = 1.5E-5 +REAL*8 convert_to_ppmv, convert_to_mixing_ratio + +!real*8, DIMENSION (1,55) :: sulf_speed + +! mean diameter for volcanic sulf accumulation mode [ m ] +! same definition in module_optical_averaging.F::optical_prep_volc() +REAL dginia_sulf +PARAMETER (dginia_sulf=0.64E-6) +! sigma-G for sulf accumulation mode +REAL sginia_sulf +PARAMETER (sginia_sulf=1.6) + + +convert_to_mixing_ratio = 1.e-6 * 96.066/28.97 +convert_to_ppmv = 1.e6 * 28.97/96.066 + +! sulfate aerosol volume median radius = mu+3 * sigma^2 +!reff_sulf=(dginia_sulf/2.0) + exp(3.*log(sginia_sulf)**2) +!reff_sulf=exp(log(dginia_sulf/2.0) + 3.*log(sginia_sulf)**2) + + +!ATTENTION here 3.5*log(...) instead of 3*log(...) +reff_sulf=exp(log(dginia_sulf/2.0) + 3.5*log(sginia_sulf)**2) + + +lmx=kte-kts+1 +do j=jts,jte +do i=its,ite + + kk=0 + DO k=kts,kte + kk=kk+1 + ! sulf in [ppmv], converting it into [kg/kg] + dsulf(1,1,kk)=convert_to_mixing_ratio * chem(i,k,j,p_sulf) + + p_mid(1,1,kk) = .01 * p_phy(i,kte-k+kts,j) + delz(1,1,kk) = dz8w(i,kte-k+kts,j) ! delz(1) = dz8w(kte), delz(lmx)=dz8w(kts) + airden(1,1,kk)= rho_phy(i,k,j) + tmp(1,1,kk) = t_phy(i,k,j) + rh(1,1,kk) = .95 + rh(1,1,kk) = MIN( .95, moist(i,k,j,p_qv) / & + (3.80*exp(17.27*(t_phy(i,k,j)-273.)/ & + (t_phy(i,k,j)-36.))/(.01*p_phy(i,k,j)))) + rh(1,1,kk)=max(1.0D-1,rh(1,1,kk)) + ENDDO + + bstl_sulf=0. + CALL sulf_settling(1,1,lmx,g,dyn_visc,dsulf,tmp,p_mid,delz,rhoso4,reff_sulf,dt,rh,airden,bstl_sulf)!,sulf_speed) + + kk = 0 + DO k = kts,kte + kk = kk+1 + chem(i,k,j,p_sulf) = convert_to_ppmv * dsulf(1,1,kk) !converting sufl [kg/kg] into [ppmv] + ENDDO + + !if(i==20.and.j==20) write(*,*)'1',sulf_speed(1,:) + + !!!!bstl_sulf [ug/m2]. total sulf grav. set [kg/m2] + !!!!sulf_fall(i,j) = sulf_fall(i,j) + sum(1.e-9 * bstl_sulf) + + !bstl_sulf in [kg/kg][m/s] + volc_diags(i,j,p_sulf_grav_setl) = volc_diags(i,j,p_sulf_grav_setl) + bstl_sulf * airden(1,1,1) * dt ![kg/m2] + +enddo ! i +enddo ! j + +END SUBROUTINE sulf_settling_driver + + +subroutine sulf_settling(imx,jmx,lmx,g0,dyn_visc,tc,tmp, & + p_mid,delz, den,reff,dt,rh,airden,bstl)!,sulf_speed) + +IMPLICIT NONE + +INTEGER, INTENT(IN) :: imx, jmx, lmx +INTEGER :: ntdt +REAL, INTENT(IN) :: dt,g0,dyn_visc +REAL*8, INTENT(IN) :: tmp(imx,jmx,lmx), delz(imx,jmx,lmx), & + rh(imx,jmx,lmx), p_mid(imx,jmx,lmx),airden(imx,jmx,lmx) + +REAL, INTENT(IN) :: den, reff + +REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx) +REAL*8, INTENT(OUT) :: bstl + +!real*8, INTENT(INOUT) :: sulf_speed(1,55) +!REAL*8, INTENT(OUT) :: grasetvel(imx,jmx,lmx) + +REAL*8 :: dt_settl +INTEGER :: ndt_settl +REAL*8 :: dzmin, vsettl, dtmax, rwet +REAL*8 :: c_stokes, free_path, c_cun, viscosity +REAL*8 :: vd_cor(lmx), vd_wk1 + +REAL*8 :: transfer_to_below_level,cell_tc +INTEGER :: k, n, i, j, l, l2 + +! for OMP: +REAL*8 :: rwet_priv, rho_priv + +REAL(4),PARAMETER :: onethird = 1.0/3.0 +REAL(4),PARAMETER :: kappa = 1.19 + +! Determine the maximum time-step satisying CFL: dt <= (dz)_min / v_settl +dzmin = MINVAL(delz(:,:,:)) +vsettl = 4.0/9.0 * g0 * den * reff**2 / dyn_visc + +ntdt = INT(dt) +dtmax = dzmin / vsettl +ndt_settl = MAX( 1,INT(ntdt/dtmax) ) + +! Limit maximum number of iterations +IF (ndt_settl > 12) ndt_settl = 12 +dt_settl = REAL(ntdt) / REAL(ndt_settl) + +! assign effective radius and density +rwet_priv = reff +rho_priv = den + +!alnsg2=log(sginia_sulf)**2 + +!$OMP PARALLEL DO & +!$OMP DEFAULT( SHARED ) & +!$OMP PRIVATE( i, j, l, l2, n, k, rwet_priv, c_stokes)& +!$OMP PRIVATE( free_path, c_cun, viscosity, rho_priv, vd_cor ) + +DO j = 1,jmx ! lat loop +DO i = 1,imx ! lon loop + +DO n = 1,ndt_settl ! time loop + + transfer_to_below_level=0 + + DO l = lmx,1,-1 ! height loop, from top + l2 = lmx - l + 1 + + rwet_priv=reff*((1+rh(i,j,l)*(kappa-1))/(1-rh(i,j,l)))**onethird ! Aerosol growth with relative humidity + + c_stokes = 1.458E-6*tmp(i,j,l)**1.5/(tmp(i,j,l) + 110.4) ! Dynamic viscosity + free_path = 1.1E-3/p_mid(i,j,l2)/SQRT(tmp(i,j,l)) ! Free path as func of pres(mb) and temp(K); order of p_mid: top->sfc + + c_cun = 1.0+free_path/rwet_priv*(1.257 + 0.4*EXP(-1.1*rwet_priv/free_path)) ! Slip correction; This formula is from Davis (1945) + viscosity = c_stokes / c_cun ! Corrected dynamic viscosity (kg/m/s) + + vd_cor(l) = 2./9.*g0*rho_priv*rwet_priv**2/viscosity ! Settling velocity, depends on temp [m/s] + + ! Update mixing ratio; order of delz: top->sfc + cell_tc=tc(i,j,l) ! cell_tc - for temporal storage [kg/kg] + vd_wk1 = dt_settl*vd_cor(l)/delz(i,j,l2) ! fraction to leave level + + tc(i,j,l) = cell_tc * (1.- vd_wk1) + transfer_to_below_level ! [kg/kg] + + IF (l==1) THEN + bstl = bstl + vd_cor(l) * cell_tc/ndt_settl ! [kg/kg][m/s] + !grasetvel(i,j,k)=vd_cor(l) ! [m/s] + ELSE + transfer_to_below_level = (cell_tc*vd_wk1)*((delz(i,j,l2)*airden(i,j,l))/(delz(i,j,l2+1)*airden(i,j,l-1))) ! [kg/kg] + ENDIF + + !sulf_speed(1,l) = vd_cor(l) + + ENDDO !l, height +ENDDO !n, time +ENDDO !i +END DO !j +!$OMP END PARALLEL DO + +END SUBROUTINE sulf_settling + + + +! A. Ukhov 30/10/2024. Large Scale scaveging Ash, SO2 and Sulf aerosols +! Based on wetdep_ls() for GOCART aerosols +subroutine wetdep_ls_volc(dt,var,rain,moist,rho,num_moist, & + num_chem,numgas,dz8w,vvel,chem_opt, & + wd_so2_sc,wd_sulf_sc,volc_diags, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte ) + + USE module_state_description, only: p_wd_ash_sc + USE module_model_constants, ONLY: mwdry +IMPLICIT NONE + +INTEGER, INTENT(IN ) :: num_chem,numgas,num_moist, & + chem_opt, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte +real, INTENT(IN ) :: dt +REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), & + INTENT(IN ) :: moist +REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & + INTENT(IN ) :: rho,dz8w,vvel +REAL, DIMENSION( ims:ime , kms:kme , jms:jme ,1:num_chem), & + INTENT(INOUT) :: var +REAL, DIMENSION( ims:ime , jms:jme ), & + INTENT(IN ) :: rain + +!Accumulated large scale scaveged so2 and sulf +REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: wd_so2_sc,wd_sulf_sc +REAL, DIMENSION( ims:ime , jms:jme ,1:9), INTENT(INOUT) :: volc_diags + +REAL, DIMENSION( its:ite , jts:jte ) :: var_sum,var_rmv +REAL, DIMENSION( its:ite , jts:jte ) :: frc,var_sum_clw,rain_clw +real :: dvar,factor,clsum,alpha!,rho_water +integer :: nv,i,j,k,km + +!write(0,*) 'in wetdepls, numgas,num_chem = ',numgas,num_chem,chem_opt +!in wetdepls, numgas,num_chem = 3 13 402 + do nv=1,num_chem + + !if(nv.le. numgas .and. nv.ne.p_sulf)cycle + alpha = .5 ! scavenging factor + + if(nv.eq.p_sulf .or. nv.eq.p_so2)alpha=1. + + do i=its,ite + do j=jts,jte + var_sum_clw(i,j)=0. + var_sum(i,j)=0. + !var_rmvl(i,:,j)=0. + frc(i,j)=0. + rain_clw(i,j)=0. + if(rain(i,j).gt.1.e-10)then + ! convert rain back to rate + rain_clw(i,j)=rain(i,j)/dt + ! total cloud water + ! + do k=1,kte-1 + dvar=max(0.,moist(i,k,j,p_qc)*rho(i,k,j)*vvel(i,k,j)*dz8w(i,k,j)) + var_sum_clw(i,j)=var_sum_clw(i,j)+dvar + var_sum(i,j)=var_sum(i,j)+var(i,k,j,nv)*rho(i,k,j) + enddo + if(var_sum(i,j).gt.1.e-10 .and. var_sum_clw(i,j).gt.1.e-10 ) then + ! assuming that frc is onstant, it is my conversion factor + ! (just like in convec. parameterization + frc(i,j)=rain_clw(i,j)/var_sum_clw(i,j) + ! write(0,*)'frc ', frc(i,j),var_sum_clw(i,j),var_sum(i,j) + frc(i,j)=max(1.e-6,min(frc(i,j),.005)) + endif + endif + enddo !j + enddo !i + ! + ! get rid of it + ! + do i=its,ite + do j=jts,jte + if(rain(i,j).gt.1.e-10 .and. var_sum(i,j).gt.1.e-10 .and. var_sum_clw(i,j).gt.1.e-10)then + do k=kts,kte-2 + if(var(i,k,j,nv).gt.1.e-16 .and. moist(i,k,j,p_qc).gt.0.)then + factor = max(0.,frc(i,j)*rho(i,k,j)*dz8w(i,k,j)*vvel(i,k,j)) + ! write(0,*)'var before ',k,km,var(i,k,j,nv),factor + dvar=alpha*factor/(1+factor)*var(i,k,j,nv) + var(i,k,j,nv)=max(1.e-16,var(i,k,j,nv)-dvar) + ! write(0,*)'var after ',km,var(i,k,j,nv),dvar + + ! accumulated scaveged volcanic SO2 and Sulf + if (chem_opt .eq. chem_volc) then + if(nv.eq.p_so2) then + wd_so2_sc(i,j) = wd_so2_sc(i,j) + dvar * rho(i,k,j) * dz8w(i,k,j)/mwdry ! mmol/m2 + endif + + if(nv.eq.p_sulf) then + wd_sulf_sc(i,j) = wd_sulf_sc(i,j) + dvar * rho(i,k,j) * dz8w(i,k,j)/mwdry ! mmol/m2 + endif + endif + + ! accumulated scaveged volcanic ash + if (chem_opt .eq. chem_vash .OR. chem_opt .eq. chem_volc) then + if(nv.eq.p_vash_1 .or. nv.eq.p_vash_2 .or. nv.eq.p_vash_3 .or. nv.eq.p_vash_4 .or. & + nv.eq.p_vash_5 .or. nv.eq.p_vash_6 .or. nv.eq.p_vash_7 .or. nv.eq.p_vash_8 .or. & + nv.eq.p_vash_9 .or. nv.eq.p_vash_10) then + volc_diags(i,j,p_wd_ash_sc) = volc_diags(i,j,p_wd_ash_sc) + 1.E-6 * dvar * rho(i,k,j) * dz8w(i,k,j) !g/m2 + endif + endif + + if(chem_opt == chem_volc_4bin) then + if(nv.eq.p_vash_7 .or. nv.eq.p_vash_8 .or. nv.eq.p_vash_9 .or. nv.eq.p_vash_10) then + volc_diags(i,j,p_wd_ash_sc) = volc_diags(i,j,p_wd_ash_sc) + 1.E-6 * dvar * rho(i,k,j) * dz8w(i,k,j) !g/m2 + endif + endif + endif + enddo !k + ! var_rmv(i,j)=var_rmv(i,j)+var_rmvl(i,j) + endif !rain... + enddo !j + enddo !i + enddo !nv +END SUBROUTINE WETDEP_LS_VOLC + + END MODULE MODULE_VASH_SETTLING diff --git a/chem/module_volc_chem.F b/chem/module_volc_chem.F new file mode 100644 index 0000000000..05e6cfd9e0 --- /dev/null +++ b/chem/module_volc_chem.F @@ -0,0 +1,337 @@ +! 17 Oct. 2024. A. Ukhov (KAUST). Oxidation of SO2 into Sulf using OH and H2O2 +! based on module_gocart_chem. +! 13 Dec. 2024. A. Ukhov (KAUST). Contribution of volcanic Sulf and ash into PM2.5 and PM10 + +! More details in: +! A. Ukhov, G. Stenchikov, J. Schnell, R. Ahmadov, U. Rizza, G. Grell, and I. Hoteit. +! Enhancing volcanic eruption simulations with the wrf-chem v4.8. +! Geoscientific Model Development, 18(23):9805–9825, 2025. doi: 10.5194/ +! gmd-18-9805-2025. URL https://gmd.copernicus.org/articles/18/9805/2025/. + +MODULE MODULE_VOLC_CHEM + + CONTAINS + + subroutine gocart_volc_chem_driver(curr_secs,dt,config_flags, & + gmt,julday,t_phy,moist, & + chem,rho_phy,dz8w,p8w,backg_oh,backg_h2o2, & + gd_cldf,dx,dy,mapfac_mx,mapfac_my,g,xlat,xlong,ttday,tcosz, & + volc_diags , & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte ) + USE module_configure + USE module_state_description + USE module_gocart_chem, only: szangle + IMPLICIT NONE + TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags + + INTEGER, INTENT(IN ) :: julday, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte + REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), & + INTENT(IN ) :: moist + REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & + INTENT(INOUT ) :: chem + REAL, DIMENSION( ims:ime , jms:jme ), & + INTENT(IN ) :: & + xlat,xlong,ttday,tcosz,mapfac_mx,mapfac_my + + REAL, DIMENSION( ims:ime , jms:jme ,1:9), INTENT(INOUT) :: volc_diags + + REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & + OPTIONAL, & + INTENT(IN ) :: gd_cldf + REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & + INTENT(IN ) :: t_phy,backg_oh,backg_h2o2,dz8w,p8w, & + rho_phy + REAL(KIND=8), INTENT(IN) :: curr_secs + REAL, INTENT(IN ) :: dt,dx,dy,g,gmt + integer :: nmx,i,j,k,imx,jmx,lmx + real*8, DIMENSION (1,1,1) :: tmp,airden,airmas,oh,h2o2, & + chlso2_oh,chlso2_aq, cldf!chpso4 + real*8, DIMENSION (1,1) :: cossza + real, DIMENSION (1,1) :: sza,cosszax + real*8, DIMENSION (1,1,1,2) :: tc + real(kind=8) :: xtime,xhour + real:: rlat,xlonn,area + real :: xmin,gmtp + integer(kind=8) :: ixhour + + imx=1 + jmx=1 + lmx=1 + nmx=4 + + xtime=curr_secs/60._8 + ixhour=int(gmt+.01,8)+int(xtime/60._8,8) + xhour=real(ixhour,8) + xmin=60.*gmt+real(xtime-xhour*60._8,8) + gmtp=mod(xhour,24._8) + gmtp=gmtp+xmin/60. + + do j=jts,jte + do i=its,ite + rlat=xlat(i,j)*3.1415926535590/180. + xlonn=xlong(i,j) + CALL szangle(1, 1, julday, gmtp, sza, cosszax,xlonn,rlat) + cossza(1,1)=cosszax(1,1) + ! + do k=kts,kte-1 + !chpso4=0. + chlso2_oh=0. + chlso2_aq=0. + if (present(gd_cldf) ) then + cldf(1,1,1)=gd_cldf(i,k,j) + else + cldf(1,1,1)=0. + endif + + if(p_qc.gt.1 .and. p_qi.gt.1)then + if(moist(i,k,j,p_qc).gt.0.or.moist(i,k,j,p_qi).gt.0.)cldf(1,1,1)=1. + elseif(p_qc.gt.1 .and. p_qi.le.1)then + if(moist(i,k,j,p_qc).gt.0.)cldf(1,1,1)=1. + endif + + area=(dx/mapfac_mx(i,j))*(dy/mapfac_my(i,j)) + airmas(1,1,1)=-(p8w(i,k+1,j)-p8w(i,k,j))*area/g ! air mass inside cell (kg) + !airmas(1,1,1)=-(p8w(i,k+1,j)-p8w(i,k,j))*dx*dx/g + + airden(1,1,1)=rho_phy(i,k,j) + tmp(1,1,1)=t_phy(i,k,j) + !oh(1,1,1)=86400./dt*cossza(1,1)*backg_oh(i,k,j)/tcosz(i,j) !backg_oh in ppv + oh(1,1,1)=86400./dt*cossza(1,1)*backg_oh(i,k,j)/tcosz(i,j) !backg_oh in molecules/cm3 + + h2o2(1,1,1)=backg_h2o2(i,k,j) !ppv + tc(1,1,1,1)=chem(i,k,j,p_so2) *1.d-6 !ppv + tc(1,1,1,2)=chem(i,k,j,p_sulf)*1.d-6 !ppv + + call chmdrv_volc_su( imx,jmx,lmx, nmx, dt, tmp, airden, airmas, & + oh, h2o2, cldf, tc, chlso2_oh, chlso2_aq)!chpso4 + + chem(i,k,j,p_so2)= tc(1,1,1,1)*1.e6 !ppmv + chem(i,k,j,p_sulf)=tc(1,1,1,2)*1.e6 !ppmv + + volc_diags(i,j,p_so2_oh_loss) = volc_diags(i,j,p_so2_oh_loss) + chlso2_oh(1,1,1) !loss (kg/dt) + volc_diags(i,j,p_so2_h2o2_loss) = volc_diags(i,j,p_so2_h2o2_loss) + chlso2_aq(1,1,1) !loss (kg/dt) + enddo + enddo + enddo + end subroutine gocart_volc_chem_driver + + + SUBROUTINE chmdrv_volc_su( imx,jmx,lmx, nmx, dt1, tmp, airden, airmas, & + oh, h2o2, cldf, tc, chlso2_oh, chlso2_aq)!chpso4 + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nmx,imx,jmx,lmx + integer :: ndt1 + real, intent(in) :: dt1 + REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: tmp, airden, airmas + REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: oh, cldf + REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: h2o2 + REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx) + REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: chlso2_oh, chlso2_aq!chpso4 + + ndt1=ifix(dt1) + if(ndt1.le.0)stop + + CALL chem_volc_so2_so4(imx,jmx,lmx, nmx, ndt1, tmp, airden, airmas, & + oh, h2o2, cldf, tc, chlso2_oh, chlso2_aq)!, chpso4,pso4_so2) + + END SUBROUTINE chmdrv_volc_su + + + SUBROUTINE chem_volc_so2_so4(imx,jmx,lmx, nmx, ndt1, tmp, airden, airmas, & + oh, h2o2, cldf, tc, chlso2_oh, chlso2_aq)!,chpso4, pso4_so2) + + ! **************************************************************************** + ! * * + ! * This is SO2 chemistry subroutine. * + ! * * + ! * SO2 loss: * + ! * SO2 + OH -> SO4 * + ! * SO2 + H2O2 -> SO4 * + ! * * + ! * If there is cloud in the gridbox (fraction = fc), then the aqueous * + ! * phase chemistry also takes place in cloud. The amount of SO2 oxidized * + ! * by H2O2 in cloud is limited by the available H2O2; * + ! * * + ! **************************************************************************** + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nmx, ndt1,imx,jmx,lmx + REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: tmp, airden,airmas + REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: cldf, oh + REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: h2o2 + REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx) + + REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: chlso2_oh, chlso2_aq !chpso4 + ! REAL*8, INTENT(OUT) :: pso4_so2(imx,jmx,lmx) + REAL*8, DIMENSION(imx,jmx,lmx) :: pso4_so2 + + + REAL*8 :: k0, kinf, kk, m, l1, l2!, ld + INTEGER :: NSO2=1 + INTEGER :: NSO4=2 + REAL, PARAMETER :: airmw = 28.97 + REAL, PARAMETER :: smw = 32.00 + ! Factor to convert AIRDEN from kgair/m3 to molecules/cm3: + REAL*8, PARAMETER :: f = 1000. / airmw * 6.022D23 * 1.0D-6 + !REAL*8, PARAMETER :: ki = 1.5D-12 + INTEGER :: i, j, l + REAL*8 :: so20,so40, tk, f1, rk1, rkt, so2_cd, fc, so2,so4 + + ! executable statements + + DO l = 1,lmx + DO j = 1,jmx + DO i = 1,imx + + so20 = tc(i,j,l,NSO2) !mol/mol + so40 = tc(i,j,l,NSO4) !mol/mol + + ! RK: SO2 + OH(g), in s-1 + tk = tmp(i,j,l) + k0 = 2.9D-31 * (298/tk)**4.1 !cm6 molecule−2 s−1 + kinf = 1.7D-12 * (298/tk)**(-0.2) !cm3 molecule−1 s−1 + m = airden(i,j,l) * f + kk = k0 * m / kinf + f1 = ( 1.0+ ( LOG10(kk) )**2 )**(-1) + + ! if OH in ppv convert it to molecules/cm3 + !rk1 = ( k0 * m / (1.0 + kk) ) * 0.6**f1 * oh(i,j,l) * airden(i,j,l)*f + + ! if OH in [molecules/cm3] do nothing + rk1 = ( k0 * m / (1.0 + kk) ) * 0.6**f1 * oh(i,j,l) + !rk2 = 0.0 + + !rk = (rk1 + rk2) + rkt = rk1 * REAL(ndt1) + + ! **************************************************************************** + ! * Update SO2 concentration after gas phase chemistry + ! **************************************************************************** + + IF (rkt > 0.0) THEN + so2_cd = so20 * EXP(-rkt) !mol/mol + l1 = (so20 - so2_cd) !* rk1/rk ! l1 loss SO2 [mol/mol] + ELSE + so2_cd = so20 + l1 = 0.0 + END IF + + ! **************************************************************************** + ! * Update SO2 concentration after cloud chemistry. * + ! * SO2 chemical loss rate = SO4 production rate (MixingRatio/timestep). * + ! **************************************************************************** + + ! Cloud chemistry (above 258K): + fc = cldf(i,j,l) + IF (fc > 0.0 .AND. so2_cd > 0.0 .AND. tk > 258.0) THEN + + IF (so2_cd > h2o2(i,j,l)) THEN !mol/mol + fc = fc * (h2o2(i,j,l)/so2_cd) + !h2o2(i,j,l) = h2o2(i,j,l) * (1.0 - cldf(i,j,l)) + !ELSE + !h2o2(i,j,l) = h2o2(i,j,l) * & + ! (1.0 - cldf(i,j,l)*so2_cd/h2o2(i,j,l)) + END IF + so2 = so2_cd * (1.0 - fc) + ! Aqueous phase SO2 loss rate (MixingRatio/timestep): + l2 = so2_cd * fc ! l2 loss SO2 !mol/mol + ELSE + so2 = so2_cd + l2 = 0.0 + END IF + + so2 = MAX(so2, 1.0D-32) + tc(i,j,l,NSO2) = so2 + + ! **************************************************************************** + ! * SO2 chemical loss rate = SO4 production rate (MixingRatio/timestep). * + ! **************************************************************************** + + pso4_so2(i,j,l) = max(0.0D0,l1 + l2) !mol/mol + so4 = so40 + pso4_so2(i,j,l) + + if(pso4_so2(i,j,l).lt.0.)then + write(0,*)'so4 routine, pso4 = ',pso4_so2(i,j,l),so4,so40 + endif + + so4 = MAX(so4, 1.0D-32) + tc(i,j,l,NSO4) = so4 + + ! --------------------------------------------------------------- + ! DIAGNOSTICS: SO2 gas-phase loss (kgS/timestep) + ! SO2 aqueous-phase loss (kgS/timestep) + ! !SO4 production (kgS/timestep) + ! --------------------------------------------------------------- + chlso2_oh(i,j,l) = chlso2_oh(i,j,l) + l1 * airmas(i,j,l) / airmw * smw + chlso2_aq(i,j,l) = chlso2_aq(i,j,l) + l2 * airmas(i,j,l) / airmw * smw + + !chpso4(i,j,l) = chpso4(i,j,l) + pso4_so2(i,j,l) * airmas(i,j,l) / airmw * smw + + END DO + END DO + END DO + + END SUBROUTINE chem_volc_so2_so4 + + ! Compute pm2_5 and pm10 from volcanic ash and sulfate + subroutine sum_pm_volc ( & + alt, chem,pm2_5_dry, pm10, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte ) + USE module_configure + USE module_state_description + USE module_data_gocartchem, only: nh4_mfac + USE module_model_constants, only: mwdry + IMPLICIT NONE + + REAL, PARAMETER :: mwso4 = 96.0576 + INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte + + REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT ) :: pm2_5_dry, pm10 + REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: alt + REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), INTENT(IN ) :: chem + + real sulfate + integer i,j,k,ii,jj + + pm2_5_dry(its:ite, kts:kte, jts:jte) = 0. + pm10(its:ite, kts:kte, jts:jte) = 0. + + do j=jts,jte + jj=min(jde-1,j) + do k=kts,kte + do i=its,ite + ii=min(ide-1,i) + sulfate = chem(ii,k,jj,p_sulf) * (mwso4/mwdry) * 1.e3 !ppmv -> (ug/kg) + + !vash10:[0.0-3.906] diam. (um) => log(2.5-0)/log(3.906-0)=0.672 + pm2_5_dry(i,k,j) = pm2_5_dry(i,k,j) + sulfate * nh4_mfac & + + 0.672 * chem(ii,k,jj,p_vash_10) + + ! vash8:[7.812-15.625] diam. (um) => (log(10)-log(7.812))/(log(15.625)-log(7.812))=0.356 + pm10(i,k,j) = pm10(i,k,j) + sulfate * nh4_mfac & + + chem(ii,k,jj,p_vash_10) + chem(ii,k,jj,p_vash_9) & + + 0.356 * chem(ii,k,jj,p_vash_8) + + !Convert from mixing ratio (ug/kg) to concentration (ug m^-3) + pm2_5_dry(i,k,j) = pm2_5_dry(i,k,j) / alt(ii,k,jj) + pm10(i,k,j) = pm10(i,k,j)/ alt(ii,k,jj) + enddo + enddo + enddo + + end subroutine sum_pm_volc + + END MODULE MODULE_volc_CHEM + diff --git a/chem/module_wetdep_ls.F b/chem/module_wetdep_ls.F index df3c4e8737..d6225a6c8d 100644 --- a/chem/module_wetdep_ls.F +++ b/chem/module_wetdep_ls.F @@ -21,7 +21,7 @@ subroutine wetdep_ls(dt,var,rain,moist,rho,num_moist, & INTENT(IN ) :: moist REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & INTENT(IN ) :: rho,dz8w,vvel - REAL, DIMENSION( ims:ime , kms:kme , jms:jme ,1:num_chem), & + REAL, DIMENSION( ims:ime , kms:kme , jms:jme ,1:num_chem), & INTENT(INOUT) :: var REAL, DIMENSION( ims:ime , jms:jme ), & INTENT(IN ) :: rain diff --git a/chem/optical_driver.F b/chem/optical_driver.F index fff5ad8644..097f715eef 100755 --- a/chem/optical_driver.F +++ b/chem/optical_driver.F @@ -114,6 +114,7 @@ SUBROUTINE optical_driver(id,curr_secs,dtstep,config_flags,haveaer,& case ( RADM2SORG, RADM2SORG_KPP, RADM2SORG_AQ, RADM2SORG_AQCHEM, & GOCART_SIMPLE, RACMSORG_KPP, RACMSORG_AQ, RACMSORG_AQCHEM_KPP, & RACM_ESRLSORG_AQCHEM_KPP, RACM_SOA_VBS_KPP, RACM_SOA_VBS_AQCHEM_KPP, & + CHEM_VOLC,CHEM_VOLC_4BIN,CHEM_VASH, & ! A. Ukhov 28 Oct 2024 RACM_SOA_VBS_HET_KPP, & GOCARTRACM_KPP, GOCARTRADM2, & RACM_ESRLSORG_KPP, MOZCART_KPP, T1_MOZCART_KPP, & @@ -146,6 +147,7 @@ SUBROUTINE optical_driver(id,curr_secs,dtstep,config_flags,haveaer,& ! A. Ukhov 09/17/24 !GOCARTRACM_KPP, GOCARTRADM2, GOCART_SIMPLE, MOZCART_KPP, & RACMSORG_KPP, RACMSORG_AQ, RACMSORG_AQCHEM_KPP, & + CHEM_VOLC, CHEM_VASH,CHEM_VOLC_4BIN, & RACM_ESRLSORG_AQCHEM_KPP, RACM_SOA_VBS_KPP, RACM_SOA_VBS_AQCHEM_KPP, & RACM_SOA_VBS_HET_KPP, CBMZSORG, CBMZSORG_AQ, & !T1_MOZCART_KPP, & CBMZ_CAM_MAM3_NOAQ, CBMZ_CAM_MAM7_NOAQ, CBMZ_CAM_MAM3_AQ, CBMZ_CAM_MAM7_AQ, & diff --git a/share/mediation_integrate.F b/share/mediation_integrate.F index cf0981e7a5..56f07d2f52 100644 --- a/share/mediation_integrate.F +++ b/share/mediation_integrate.F @@ -209,8 +209,15 @@ SUBROUTINE med_before_solve_io ( grid , config_flags ) ELSE IF( ialarm .EQ. AUXINPUT13_ALARM .AND. config_flags%chem_opt > 0 ) THEN IF( config_flags%emiss_opt_vol /= 0 ) THEN IF( WRFU_AlarmIsRinging( grid%alarms( ialarm ), rc=rc ) ) THEN - call wrf_debug(15,' CALL med_read_wrf_volc_emiss ') - CALL med_read_wrf_volc_emiss ( grid , config_flags ) + + IF( config_flags%emiss_opt_vol == 3 ) THEN + call wrf_debug(15,' CALL med_read_wrf_volc_emiss_opt_vol_eq_3 ') + CALL med_read_wrf_volc_emiss_opt_vol_eq_3 ( grid , config_flags ) + ELSE + call wrf_debug(15,' CALL med_read_wrf_volc_emiss ') + CALL med_read_wrf_volc_emiss ( grid , config_flags ) + ENDIF + CALL WRFU_AlarmRingerOff( grid%alarms( ialarm ), rc=rc ) call wrf_debug(15,' Back from CALL med_read_wrf_volc_emiss ') ENDIF @@ -2388,6 +2395,109 @@ SUBROUTINE med_read_wrf_volc_emiss ( grid , config_flags ) END SUBROUTINE med_read_wrf_volc_emiss +SUBROUTINE med_read_wrf_volc_emiss_opt_vol_eq_3 ( grid , config_flags ) + ! Driver layer + USE module_domain , ONLY : domain , domain_clock_get + USE module_io_domain + USE module_timing + USE module_configure , ONLY : grid_config_rec_type + ! Model layer + USE module_bc_time_utilities +#ifdef DM_PARALLEL + USE module_dm +#endif + USE module_date_time + USE module_utility + + IMPLICIT NONE + + ! Arguments + TYPE(domain) :: grid + + TYPE (grid_config_rec_type) , INTENT(IN) :: config_flags + + ! Local data + !LOGICAL, EXTERNAL :: wrf_dm_on_monitor + INTEGER :: ierr + CHARACTER (LEN=256) :: message, current_date_char + CHARACTER (LEN=256) :: inpname + + INTEGER :: day_of_year, julyr + REAL :: gmt + + character(len=50) :: trimmed_auxinput13_inname + INTEGER :: trimmed_day_of_year, trimmed_julyr + REAL :: trimmed_gmt + + logical :: first_call + save first_call + data first_call / .true. / + +#include "wrf_io_flags.h" + + CALL domain_clock_get( grid, current_timestr=current_date_char ) + + IF( grid%auxinput13_oid == 0 ) then + CALL get_julgmt ( current_date_char , julyr , day_of_year , gmt ) + + !print *, 'current_date_char: ', current_date_char + !print *, 'julyr, day_of_year, gmt : ', julyr, day_of_year, gmt + trimmed_auxinput13_inname = grid%auxinput13_inname(14:) + + CALL get_julgmt ( trimmed_auxinput13_inname , trimmed_julyr , trimmed_day_of_year , trimmed_gmt ) + !print *, 'trimmed_current_date_char: ', trimmed_auxinput13_inname + !print *, 'trimmed_julyr, trimmed_day_of_year, trimmed_gmt : ', trimmed_julyr, trimmed_day_of_year, trimmed_gmt + + IF (trimmed_julyr .eq. julyr .and. trimmed_day_of_year .eq. day_of_year .and. gmt.ge.trimmed_gmt .and. first_call) THEN + inpname=TRIM(grid%auxinput13_inname) + WRITE(message,*)'mediation_integrate: med_read_wrf_volc_emiss_opt_vol_eq_3: Open file ',TRIM(inpname) + CALL wrf_message( TRIM(message) ) + + CALL open_r_dataset ( grid%auxinput13_oid, TRIM(inpname) , grid , config_flags, "DATASET=AUXINPUT13", ierr ) + IF ( ierr .NE. 0 ) THEN + WRITE( message , * ) 'med_read_wrf_volc_emiss_opt_vol_eq_3: error opening ', TRIM( inpname ) + CALL wrf_error_fatal( TRIM( message ) ) + ENDIF + + WRITE(message,*)'mediation_integrate: med_read_wrf_volc_emiss_opt_vol_eq_3: Read volcanic ash, so2, sulf and qv emissions ',& + TRIM(current_date_char),' frame ',grid%emissframes + CALL wrf_message( TRIM(message) ) + ! + ! Read the emissions data. + ! + CALL wrf_debug (100 , 'mediation_integrate: calling input_auxinput13' ) + CALL input_auxinput13 ( grid%auxinput13_oid, grid , config_flags , ierr ) + ! + grid%emissframes = grid%emissframes + 1 + ENDIF + ELSE + WRITE(message,*)'mediation_integrate: med_read_wrf_volc_emiss_opt_vol_eq_3: Read volcanic ash, so2, sulf and qv emissions ',& + TRIM(current_date_char),' frame ',grid%emissframes + CALL wrf_message( TRIM(message) ) +! +! Read the emissions data. +! + CALL wrf_debug (100 , 'mediation_integrate: calling input_auxinput13' ) + CALL input_auxinput13 ( grid%auxinput13_oid, grid , config_flags , ierr ) +! +! If reached the indicated number of frames in the emissions file, close it. +! + grid%emissframes = grid%emissframes + 1 + + IF ( grid%emissframes >= config_flags%frames_per_auxinput13 ) THEN + CALL close_dataset ( grid%auxinput13_oid , config_flags , "DATASET=AUXINPUT13" ) + grid%auxinput13_oid = 0 + WRITE(message,*)'mediation_integrate: med_read_wrf_volc_emiss_opt_vol_eq_3: close emissions file ',TRIM(inpname),& + TRIM(current_date_char),' frame ',grid%emissframes + CALL wrf_message( TRIM(message) ) + + first_call = .false. + ENDIF + ENDIF !( grid%auxinput13_oid == 0 ) + + CALL wrf_debug (100 , 'mediation_integrate: med_read_wrf_volc_emiss_opt_vol_eq_3: exit' ) + + END SUBROUTINE med_read_wrf_volc_emiss_opt_vol_eq_3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE med_read_wrf_chem_emissopt3 ( grid , config_flags ) ! Driver layer