!>\file seas_mod.F90 !! This file contains the sea salt emission module. module seas_mod use machine , only : kind_phys ! use chem_rc_mod, only : chem_rc_test ! use chem_tracers_mod, only : p_seas_1, p_seas_2, p_seas_3, p_seas_4, p_seas_5, & ! p_eseas1, p_eseas2, p_eseas3, p_eseas4, p_eseas5, & ! config => chem_config use seas_data_mod use seas_ngac_mod implicit none integer, parameter :: CHEM_OPT_GOCART = 300 integer, parameter :: chem_opt = 300 ! -- NGAC parameters integer, parameter :: emission_scheme = 3 ! GEOSS 2012 private public :: gocart_seasalt_driver CONTAINS subroutine gocart_seasalt_driver(dt,alt,t_phy,u_phy, & v_phy,chem,rho_phy,dz8w,u10,v10,ustar,p8w,tsk, & xland,xlat,xlong,area,g,emis_seas,pi, & seashelp,num_emis_seas,num_chem,seas_opt, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) INTEGER, INTENT(IN ) :: num_emis_seas,num_chem, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte,seas_opt ! REAL(kind=kind_phys), DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), & ! INTENT(IN ) :: moist REAL(kind=kind_phys), DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & INTENT(INOUT ) :: chem REAL(kind=kind_phys), DIMENSION( ims:ime, 1, jms:jme,num_emis_seas), & INTENT(OUT ) :: & emis_seas REAL(kind=kind_phys), DIMENSION( ims:ime , jms:jme ) , & INTENT(IN ) :: & u10, & v10, & ustar,tsk, & xland, & xlat, & xlong,area REAL(kind=kind_phys), DIMENSION( ims:ime , jms:jme ), & INTENT(OUT ) :: seashelp REAL(kind=kind_phys), DIMENSION( ims:ime , kms:kme , jms:jme ), & INTENT(IN ) :: & alt, & t_phy, & dz8w,p8w, & u_phy,v_phy,rho_phy REAL(kind=kind_phys), INTENT(IN ) :: dt,g,pi ! integer, parameter :: p_seas_1=15 integer, parameter :: p_seas_2=16 integer, parameter :: p_seas_3=17 integer, parameter :: p_seas_4=18 integer, parameter :: p_seas_5=19 integer, parameter :: p_eseas1=1 integer, parameter :: p_eseas2=2 integer, parameter :: p_eseas3=3 integer, parameter :: p_eseas4=4 integer, parameter :: p_eseas5=5 ! ! local variables ! integer :: ipr,i,j,imx,jmx,lmx,n,rc,chem_config integer,dimension (1,1) :: ilwi real(kind=kind_phys) :: fsstemis, memissions, nemissions, tskin_c, ws10m real(kind=kind_phys) :: delp real(kind=kind_phys), DIMENSION (number_ss_bins) :: tc,bems real(kind=kind_phys), dimension (1,1) ::w10m,airmas,tskin real(kind=kind_phys), dimension (1) :: dxy real(kind=kind_phys), dimension(1,1,1) :: airmas1 real(kind=kind_phys), dimension(1,1,1,number_ss_bins) :: tc1 real(kind=kind_phys), dimension(1,1,number_ss_bins) :: bems1 ! ! local parameters ! real(kind=kind_phys), parameter :: conver = 1.e-9_kind_phys real(kind=kind_phys), parameter :: converi = 1.e+9_kind_phys ! ! number of dust bins ! imx=1 jmx=1 lmx=1 chem_config=CHEM_OPT_GOCART emis_seas = 0. ! select case (config % chem_opt) select case (chem_opt) case (304, 316, 317) seashelp(:,:)=0. do j=jts,jte do i=its,ite ! ! don't do dust over water!!! ! if(xland(i,j).lt.0.5)then ilwi(1,1)=0 tc(1)=chem(i,kts,j,p_seas_1)*conver tc(2)=1.e-30_kind_phys tc(3)=chem(i,kts,j,p_seas_2)*conver tc(4)=1.e-30_kind_phys w10m(1,1)=sqrt(u10(i,j)*u10(i,j)+v10(i,j)*v10(i,j)) tskin(1,1)=tsk(i,j) delp = p8w(i,kts,j)-p8w(i,kts+1,j) airmas(1,1)=area(i,j) * delp / g ! ! we don't trust the u10,v10 values, is model layers are very thin near surface ! if(dz8w(i,kts,j).lt.12.)w10m=sqrt(u_phy(i,kts,j)*u_phy(i,kts,j)+v_phy(i,kts,j)*v_phy(i,kts,j)) ! dxy(1)=area(i,j) ipr=0 airmas1(1,1,1) = airmas(1,1) tc1(1,1,1,:) = tc bems1(1,1,:) = bems call source_ss( imx, jmx, lmx, number_ss_bins, dt, tc1, pi, ilwi, dxy, w10m, airmas1, bems1,ipr) tc = tc1(1,1,1,:) chem(i,kts,j,p_seas_1)=(tc(1)+.75*tc(2))*converi chem(i,kts,j,p_seas_2)=(tc(3)+.25*tc(2))*converi seashelp(i,j)=tc(2)*converi endif enddo enddo case default select case (seas_opt) case (1) ! -- original GOCART sea salt scheme do j = jts, jte do i = its, ite ! -- only use sea salt scheme over water if (xland(i,j) < 0.5) then ! -- compute auxiliary variables delp = p8w(i,kts,j)-p8w(i,kts+1,j) if (dz8w(i,kts,j) < 12.) then w10m = sqrt(u_phy(i,kts,j)*u_phy(i,kts,j)+v_phy(i,kts,j)*v_phy(i,kts,j)) else w10m = sqrt(u10(i,j)*u10(i,j)+v10(i,j)*v10(i,j)) end if ilwi(1,1)=0 tc = 0. tskin(1,1)=tsk(i,j) airmas(1,1)=area(i,j) * delp / g dxy(1)=area(i,j) ipr=0 airmas1(1,1,1) = airmas(1,1) tc1(1,1,1,:) = tc bems1(1,1,:) = bems call source_ss( imx,jmx,lmx,number_ss_bins, dt, tc1, pi, ilwi, dxy, w10m, airmas1, bems1,ipr) tc = tc1(1,1,1,:) bems = bems1(1,1,:) ! -- add sea salt emission increments to existing airborne concentrations chem(i,kts,j,p_seas_1) = chem(i,kts,j,p_seas_1) + tc(1)*converi chem(i,kts,j,p_seas_2) = chem(i,kts,j,p_seas_2) + tc(2)*converi chem(i,kts,j,p_seas_3) = chem(i,kts,j,p_seas_3) + tc(3)*converi chem(i,kts,j,p_seas_4) = chem(i,kts,j,p_seas_4) + tc(4)*converi chem(i,kts,j,p_seas_5) = chem(i,kts,j,p_seas_5) + tc(5)*converi !print*,'hli tc(2),chem(i,kts,j,p_seas_2)',tc(2),chem(i,kts,j,p_seas_2) ! for output diagnostics emis_seas(i,1,j,p_eseas1) = bems(1) emis_seas(i,1,j,p_eseas2) = bems(2) emis_seas(i,1,j,p_eseas3) = bems(3) emis_seas(i,1,j,p_eseas4) = bems(4) emis_seas(i,1,j,p_eseas5) = bems(5) end if end do end do case (2) ! -- NGAC sea salt scheme do j = jts, jte do i = its, ite ! -- only use sea salt scheme over water if (xland(i,j) < 0.5) then ! -- compute auxiliary variables delp = p8w(i,kts,j)-p8w(i,kts+1,j) if (dz8w(i,kts,j) < 12.) then ws10m = sqrt(u_phy(i,kts,j)*u_phy(i,kts,j)+v_phy(i,kts,j)*v_phy(i,kts,j)) else ws10m = sqrt(u10(i,j)*u10(i,j)+v10(i,j)*v10(i,j)) end if ! -- compute NGAC SST correction tskin_c = tsk(i,j) - 273.15 tskin_c = min(max(tskin_c, -0.1), 36.0) ! temperature range (0, 36) C fsstemis = -1.107211 & - tskin_c*(0.010681+0.002276*tskin_c) & + 60.288927/(40.0 - tskin_c) fsstemis = min(max(fsstemis, 0.0), 7.0) do n = 1, number_ss_bins memissions = 0. nemissions = 0. call SeasaltEmission( ra(n), rb(n), emission_scheme, & ws10m, ustar(i,j), pi, memissions, nemissions, rc ) ! if (chem_rc_test((rc /= 0), msg="Error in NGAC sea salt scheme", & ! file=__FILE__, line=__LINE__)) return bems(n) = emission_scale(n) * fsstemis * memissions tc(n) = bems(n) * dt * g / delp end do ! -- add sea salt emission increments to existing airborne concentrations chem(i,kts,j,p_seas_1) = chem(i,kts,j,p_seas_1) + tc(1)*converi chem(i,kts,j,p_seas_2) = chem(i,kts,j,p_seas_2) + tc(2)*converi chem(i,kts,j,p_seas_3) = chem(i,kts,j,p_seas_3) + tc(3)*converi chem(i,kts,j,p_seas_4) = chem(i,kts,j,p_seas_4) + tc(4)*converi chem(i,kts,j,p_seas_5) = chem(i,kts,j,p_seas_5) + tc(5)*converi ! for output diagnostics kg/m2/s emis_seas(i,1,j,p_eseas1) = bems(1) emis_seas(i,1,j,p_eseas2) = bems(2) emis_seas(i,1,j,p_eseas3) = bems(3) emis_seas(i,1,j,p_eseas4) = bems(4) emis_seas(i,1,j,p_eseas5) = bems(5) end if end do end do case default ! -- no sea salt scheme end select end select end subroutine gocart_seasalt_driver SUBROUTINE source_ss(imx,jmx,lmx,nmx, dt1, tc, pi, & ilwi, dxy, w10m, airmas, & bems,ipr) ! **************************************************************************** ! * Evaluate the source of each seasalt particles size classes (kg/m3) ! * by soil emission. ! * Input: ! * SSALTDEN Sea salt density (kg/m3) ! * DXY Surface of each grid cell (m2) ! * NDT1 Time step (s) ! * W10m Velocity at the anemometer level (10meters) (m/s) ! * ! * Output: ! * DSRC Source of each sea salt bins (kg/timestep/cell) ! * ! * ! * Number flux density: Original formula by Monahan et al. (1986) adapted ! * by Sunling Gong (JGR 1997 (old) and GBC 2003 (new)). The new version is ! * to better represent emission of sub-micron sea salt particles. ! ! * dFn/dr = c1*u10**c2/(r**A) * (1+c3*r**c4)*10**(c5*exp(-B**2)) ! * where B = (b1 -log(r))/b2 ! * see c_old, c_new, b_old, b_new below for the constants. ! * number fluxes are at 80% RH. ! * ! * To calculate the flux: ! * 1) Calculate dFn based on Monahan et al. (1986) and Gong (2003) ! * 2) Assume that wet radius r at 80% RH = dry radius r_d *frh ! * 3) Convert particles flux to mass flux : ! * dFM/dr_d = 4/3*pi*rho_d*r_d^3 *(dr/dr_d) * dFn/dr ! * = 4/3*pi*rho_d*r_d^3 * frh * dFn/dr ! * where rho_p is particle density [kg/m3] ! * The factor 1.e-18 is to convert in micro-meter r_d^3 ! **************************************************************************** IMPLICIT NONE INTEGER, INTENT(IN) :: nmx,imx,jmx,lmx,ipr INTEGER, INTENT(IN) :: ilwi(imx,jmx) REAL(kind=kind_phys), INTENT(IN) :: dxy(jmx), w10m(imx,jmx), pi REAL(kind=kind_phys), INTENT(IN) :: airmas(imx,jmx,lmx) REAL(kind=kind_phys), INTENT(INOUT) :: tc(imx,jmx,lmx,nmx) REAL(kind=kind_phys), INTENT(OUT) :: bems(imx,jmx,nmx) REAL(kind=kind_phys) :: c0(5), b0(2) ! REAL(kind=kind_phys), PARAMETER :: c_old(5)=(/1.373, 3.41, 0.057, 1.05, 1.190/) ! REAL(kind=kind_phys), PARAMETER :: c_new(5)=(/1.373, 3.41, 0.057, 3.45, 1.607/) ! Change suggested by MC REAL(kind=kind_phys), PARAMETER :: c_old(5)=(/1.373, 3.2, 0.057, 1.05, 1.190/) REAL(kind=kind_phys), PARAMETER :: c_new(5)=(/1.373, 3.2, 0.057, 3.45, 1.607/) REAL(kind=kind_phys), PARAMETER :: b_old(2)=(/0.380, 0.650/) REAL(kind=kind_phys), PARAMETER :: b_new(2)=(/0.433, 0.433/) REAL(kind=kind_phys), PARAMETER :: dr=5.0D-2 ! um REAL(kind=kind_phys), PARAMETER :: theta=30.0 ! Swelling coefficient frh (d rwet / d rd) !!! REAL(kind=kind_phys), PARAMETER :: frh = 1.65 REAL(kind=kind_phys), PARAMETER :: frh = 2.d0 LOGICAL, PARAMETER :: old=.TRUE., new=.FALSE. REAL(kind=kind_phys) :: rho_d, r0, r1, r, r_w, a, b, dfn, r_d, dfm, src INTEGER :: i, j, n, nr, ir REAL(kind=kind_phys) :: dt1,fudge_fac REAL(kind=kind_phys) :: tcmw(nmx), ar(nmx), tcvv(nmx) REAL(kind=kind_phys) :: ar_wetdep(nmx), kc(nmx) CHARACTER(LEN=20) :: tcname(nmx), tcunits(nmx) LOGICAL :: aerosol(nmx) REAL(kind=kind_phys) :: tc1(imx,jmx,lmx,nmx) REAL(kind=kind_phys), TARGET :: tcms(imx,jmx,lmx,nmx) ! tracer mass (kg; kgS for sulfur case) REAL(kind=kind_phys), TARGET :: tcgm(imx,jmx,lmx,nmx) ! g/m3 !----------------------------------------------------------------------- ! sea salt specific !----------------------------------------------------------------------- ! REAL(kind=kind_phys), DIMENSION(nmx) :: ra, rb ! REAL(kind=kind_phys) :: ch_ss(nmx,12) !----------------------------------------------------------------------- ! emissions (input) !----------------------------------------------------------------------- REAL(kind=kind_phys) :: e_an(imx,jmx,2,nmx), e_bb(imx,jmx,nmx), & e_ac(imx,jmx,lmx,nmx) !----------------------------------------------------------------------- ! diagnostics (budget) !----------------------------------------------------------------------- ! ! tendencies per time step and process ! REAL(kind=kind_phys), TARGET :: bems(imx,jmx,nmx), bdry(imx,jmx,nmx), bstl(imx,jmx,nmx) ! REAL(kind=kind_phys), TARGET :: bwet(imx,jmx,nmx), bcnv(imx,jmx,nmx)! ! ! integrated tendencies per process ! REAL(kind=kind_phys), TARGET :: tems(imx,jmx,nmx), tstl(imx,jmx,nmx) ! REAL(kind=kind_phys), TARGET :: tdry(imx,jmx,nmx), twet(imx,jmx,nmx), tcnv(imx,jmx,nmx) ! global mass balance per time step REAL(kind=kind_phys) :: tmas0(nmx), tmas1(nmx) REAL(kind=kind_phys) :: dtems(nmx), dttrp(nmx), dtdif(nmx), dtcnv(nmx) REAL(kind=kind_phys) :: dtwet(nmx), dtdry(nmx), dtstl(nmx) REAL(kind=kind_phys) :: dtems2(nmx), dttrp2(nmx), dtdif2(nmx), dtcnv2(nmx) REAL(kind=kind_phys) :: dtwet2(nmx), dtdry2(nmx), dtstl2(nmx) ! detailed integrated budgets for individual emissions REAL(kind=kind_phys), TARGET :: ems_an(imx,jmx,nmx), ems_bb(imx,jmx,nmx), ems_tp(imx,jmx) REAL(kind=kind_phys), TARGET :: ems_ac(imx,jmx,lmx,nmx) REAL(kind=kind_phys), TARGET :: ems_co(imx,jmx,nmx) ! executable statements ! decrease seasalt emissions (Colarco et al. 2010) ! !fudge_fac= 1. !.5 !fudge_fac= .5 !lzhang fudge_fac= .25 !lzhang ! DO n = 1,nmx bems(:,:,n) = 0.0 rho_d = den_seas(n) r0 = ra(n)*frh r1 = rb(n)*frh r = r0 nr = INT((r1-r0)/dr+.001) DO ir = 1,nr r_w = r + dr*0.5 r = r + dr IF (new) THEN a = 4.7*(1.0 + theta*r_w)**(-0.017*r_w**(-1.44)) c0 = c_new b0 = b_new ELSE a = 3.0 c0 = c_old b0 = b_old END IF ! b = (b0(1) - LOG10(r_w))/b0(2) dfn = (c0(1)/r_w**a)*(1.0 + c0(3)*r_w**c0(4))* & 10**(c0(5)*EXP(-(b**2))) r_d = r_w/frh*1.0D-6 ! um -> m dfm = 4.0/3.0*pi*r_d**3*rho_d*frh*dfn*dr*dt1 ! 3600 !dt1 DO i = 1,imx DO j = 1,jmx ! IF (water(i,j) > 0.0) THEN IF (ilwi(i,j) == 0) THEN ! src = dfm*dxy(j)*water(i,j)*w10m(i,j)**c0(2) src = dfm*dxy(j)*w10m(i,j)**c0(2) ! src = ch_ss(n,dt(1)%mn)*dfm*dxy(j)*w10m(i,j)**c0(2) tc(i,j,1,n) = tc(i,j,1,n) + fudge_fac*src/airmas(i,j,1) ELSE src = 0.0 END IF bems(i,j,n) = bems(i,j,n) + src*fudge_fac/(dxy(j)*dt1) !kg/m2/s END DO ! i END DO ! j END DO ! ir END DO ! n END SUBROUTINE source_ss end module seas_mod