| 1 | MODULE GOCART_SEASALT |
|---|
| 2 | |
|---|
| 3 | CONTAINS |
|---|
| 4 | subroutine gocart_seasalt_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u_phy, & |
|---|
| 5 | v_phy,chem,rho_phy,dz8w,u10,v10,p8w, & |
|---|
| 6 | xland,xlat,xlong,dx,g,emis_seas, & |
|---|
| 7 | ids,ide, jds,jde, kds,kde, & |
|---|
| 8 | ims,ime, jms,jme, kms,kme, & |
|---|
| 9 | its,ite, jts,jte, kts,kte ) |
|---|
| 10 | USE module_configure |
|---|
| 11 | USE module_state_description |
|---|
| 12 | USE module_model_constants, ONLY: mwdry |
|---|
| 13 | IMPLICIT NONE |
|---|
| 14 | TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags |
|---|
| 15 | |
|---|
| 16 | INTEGER, INTENT(IN ) :: julday, ktau, & |
|---|
| 17 | ids,ide, jds,jde, kds,kde, & |
|---|
| 18 | ims,ime, jms,jme, kms,kme, & |
|---|
| 19 | its,ite, jts,jte, kts,kte |
|---|
| 20 | REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), & |
|---|
| 21 | INTENT(IN ) :: moist |
|---|
| 22 | REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & |
|---|
| 23 | INTENT(INOUT ) :: chem |
|---|
| 24 | REAL, DIMENSION( ims:ime, 1, jms:jme,num_emis_seas),OPTIONAL,& |
|---|
| 25 | INTENT(INOUT ) :: & |
|---|
| 26 | emis_seas |
|---|
| 27 | REAL, DIMENSION( ims:ime , jms:jme ) , & |
|---|
| 28 | INTENT(IN ) :: & |
|---|
| 29 | u10, & |
|---|
| 30 | v10, & |
|---|
| 31 | xland, & |
|---|
| 32 | xlat, & |
|---|
| 33 | xlong |
|---|
| 34 | REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & |
|---|
| 35 | INTENT(IN ) :: & |
|---|
| 36 | alt, & |
|---|
| 37 | t_phy, & |
|---|
| 38 | dz8w,p8w, & |
|---|
| 39 | u_phy,v_phy,rho_phy |
|---|
| 40 | |
|---|
| 41 | REAL, INTENT(IN ) :: dt,dx,g |
|---|
| 42 | ! |
|---|
| 43 | ! local variables |
|---|
| 44 | ! |
|---|
| 45 | integer :: ipr,nmx,i,j,k,ndt,imx,jmx,lmx |
|---|
| 46 | integer,dimension (1,1) :: ilwi |
|---|
| 47 | real*8, DIMENSION (4) :: tc,bems |
|---|
| 48 | real*8, dimension (1,1) :: w10m,gwet,airden,airmas |
|---|
| 49 | real*8, dimension (1) :: dxy |
|---|
| 50 | real*8 conver,converi |
|---|
| 51 | conver=1.d-9 |
|---|
| 52 | converi=1.d9 |
|---|
| 53 | ! |
|---|
| 54 | ! number of dust bins |
|---|
| 55 | ! |
|---|
| 56 | imx=1 |
|---|
| 57 | jmx=1 |
|---|
| 58 | lmx=1 |
|---|
| 59 | nmx=4 |
|---|
| 60 | k=kts |
|---|
| 61 | do j=jts,jte |
|---|
| 62 | do i=its,ite |
|---|
| 63 | ! |
|---|
| 64 | ! donṫ do dust over water!!! |
|---|
| 65 | ! |
|---|
| 66 | if(xland(i,j).gt.1.5)then |
|---|
| 67 | ilwi(1,1)=0 |
|---|
| 68 | tc(1)=chem(i,kts,j,p_seas_1)*conver |
|---|
| 69 | tc(2)=chem(i,kts,j,p_seas_2)*conver |
|---|
| 70 | tc(3)=chem(i,kts,j,p_seas_3)*conver |
|---|
| 71 | tc(4)=chem(i,kts,j,p_seas_4)*conver |
|---|
| 72 | w10m(1,1)=sqrt(u10(i,j)*u10(i,j)+v10(i,j)*v10(i,j)) |
|---|
| 73 | airmas(1,1)=-(p8w(i,kts+1,j)-p8w(i,kts,j))*dx*dx/g |
|---|
| 74 | ! |
|---|
| 75 | ! we donṫ trust the u10,v10 values, is model layers are very thin near surface |
|---|
| 76 | ! |
|---|
| 77 | 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)) |
|---|
| 78 | ! |
|---|
| 79 | dxy(1)=dx*dx |
|---|
| 80 | ipr=0 |
|---|
| 81 | |
|---|
| 82 | call source_ss( imx,jmx,lmx,nmx, dt, tc,ilwi, dxy, w10m, airmas, bems,ipr) |
|---|
| 83 | chem(i,kts,j,p_seas_1)=tc(1)*converi |
|---|
| 84 | chem(i,kts,j,p_seas_2)=tc(2)*converi |
|---|
| 85 | chem(i,kts,j,p_seas_3)=tc(3)*converi |
|---|
| 86 | chem(i,kts,j,p_seas_4)=tc(4)*converi |
|---|
| 87 | ! for output diagnostics |
|---|
| 88 | emis_seas(i,1,j,p_edust1)=bems(1) |
|---|
| 89 | emis_seas(i,1,j,p_edust2)=bems(2) |
|---|
| 90 | emis_seas(i,1,j,p_edust3)=bems(3) |
|---|
| 91 | emis_seas(i,1,j,p_edust4)=bems(4) |
|---|
| 92 | endif |
|---|
| 93 | enddo |
|---|
| 94 | enddo |
|---|
| 95 | ! |
|---|
| 96 | |
|---|
| 97 | end subroutine gocart_seasalt_driver |
|---|
| 98 | ! |
|---|
| 99 | SUBROUTINE source_ss(imx,jmx,lmx,nmx, dt1, tc, & |
|---|
| 100 | ilwi, dxy, w10m, airmas, & |
|---|
| 101 | bems,ipr) |
|---|
| 102 | |
|---|
| 103 | ! **************************************************************************** |
|---|
| 104 | ! * Evaluate the source of each seasalt particles size classes (kg/m3) |
|---|
| 105 | ! * by soil emission. |
|---|
| 106 | ! * Input: |
|---|
| 107 | ! * SSALTDEN Sea salt density (kg/m3) |
|---|
| 108 | ! * DXY Surface of each grid cell (m2) |
|---|
| 109 | ! * NDT1 Time step (s) |
|---|
| 110 | ! * W10m Velocity at the anemometer level (10meters) (m/s) |
|---|
| 111 | ! * |
|---|
| 112 | ! * Output: |
|---|
| 113 | ! * DSRC Source of each sea salt bins (kg/timestep/cell) |
|---|
| 114 | ! * |
|---|
| 115 | ! * |
|---|
| 116 | ! * Number flux density: Original formula by Monahan et al. (1986) adapted |
|---|
| 117 | ! * by Sunling Gong (JGR 1997 (old) and GBC 2003 (new)). The new version is |
|---|
| 118 | ! * to better represent emission of sub-micron sea salt particles. |
|---|
| 119 | ! |
|---|
| 120 | ! * dFn/dr = c1*u10**c2/(r**A) * (1+c3*r**c4)*10**(c5*exp(-B**2)) |
|---|
| 121 | ! * where B = (b1 -log(r))/b2 |
|---|
| 122 | ! * see c_old, c_new, b_old, b_new below for the constants. |
|---|
| 123 | ! * number fluxes are at 80% RH. |
|---|
| 124 | ! * |
|---|
| 125 | ! * To calculate the flux: |
|---|
| 126 | ! * 1) Calculate dFn based on Monahan et al. (1986) and Gong (2003) |
|---|
| 127 | ! * 2) Assume that wet radius r at 80% RH = dry radius r_d *frh |
|---|
| 128 | ! * 3) Convert particles flux to mass flux : |
|---|
| 129 | ! * dFM/dr_d = 4/3*pi*rho_d*r_d^3 *(dr/dr_d) * dFn/dr |
|---|
| 130 | ! * = 4/3*pi*rho_d*r_d^3 * frh * dFn/dr |
|---|
| 131 | ! * where rho_p is particle density [kg/m3] |
|---|
| 132 | ! * The factor 1.e-18 is to convert in micro-meter r_d^3 |
|---|
| 133 | ! **************************************************************************** |
|---|
| 134 | |
|---|
| 135 | |
|---|
| 136 | USE module_data_gocart_seas |
|---|
| 137 | |
|---|
| 138 | IMPLICIT NONE |
|---|
| 139 | |
|---|
| 140 | INTEGER, INTENT(IN) :: nmx,imx,jmx,lmx,ipr |
|---|
| 141 | INTEGER, INTENT(IN) :: ilwi(imx,jmx) |
|---|
| 142 | REAL*8, INTENT(IN) :: dxy(jmx), w10m(imx,jmx) |
|---|
| 143 | REAL*8, INTENT(IN) :: airmas(imx,jmx,lmx) |
|---|
| 144 | REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx) |
|---|
| 145 | REAL*8, INTENT(OUT) :: bems(imx,jmx,nmx) |
|---|
| 146 | |
|---|
| 147 | REAL*8 :: c0(5), b0(2) |
|---|
| 148 | ! REAL*8, PARAMETER :: c_old(5)=(/1.373, 3.41, 0.057, 1.05, 1.190/) |
|---|
| 149 | ! REAL*8, PARAMETER :: c_new(5)=(/1.373, 3.41, 0.057, 3.45, 1.607/) |
|---|
| 150 | ! Change suggested by MC |
|---|
| 151 | REAL*8, PARAMETER :: c_old(5)=(/1.373, 3.2, 0.057, 1.05, 1.190/) |
|---|
| 152 | REAL*8, PARAMETER :: c_new(5)=(/1.373, 3.2, 0.057, 3.45, 1.607/) |
|---|
| 153 | REAL*8, PARAMETER :: b_old(2)=(/0.380, 0.650/) |
|---|
| 154 | REAL*8, PARAMETER :: b_new(2)=(/0.433, 0.433/) |
|---|
| 155 | REAL*8, PARAMETER :: dr=5.0D-2 ! um |
|---|
| 156 | REAL*8, PARAMETER :: theta=30.0 |
|---|
| 157 | ! Swelling coefficient frh (d rwet / d rd) |
|---|
| 158 | !!! REAL*8, PARAMETER :: frh = 1.65 |
|---|
| 159 | REAL*8, PARAMETER :: frh = 2.d0 |
|---|
| 160 | LOGICAL, PARAMETER :: old=.TRUE., new=.FALSE. |
|---|
| 161 | REAL*8 :: rho_d, r0, r1, r, r_w, a, b, dfn, r_d, dfm, src |
|---|
| 162 | INTEGER :: i, j, n, nr, ir |
|---|
| 163 | REAL :: dt1 |
|---|
| 164 | |
|---|
| 165 | |
|---|
| 166 | REAL*8 :: tcmw(nmx), ar(nmx), tcvv(nmx) |
|---|
| 167 | REAL*8 :: ar_wetdep(nmx), kc(nmx) |
|---|
| 168 | CHARACTER(LEN=20) :: tcname(nmx), tcunits(nmx) |
|---|
| 169 | LOGICAL :: aerosol(nmx) |
|---|
| 170 | |
|---|
| 171 | |
|---|
| 172 | REAL*8 :: tc1(imx,jmx,lmx,nmx) |
|---|
| 173 | REAL*8, TARGET :: tcms(imx,jmx,lmx,nmx) ! tracer mass (kg; kgS for sulfur case) |
|---|
| 174 | REAL*8, TARGET :: tcgm(imx,jmx,lmx,nmx) ! g/m3 |
|---|
| 175 | |
|---|
| 176 | !----------------------------------------------------------------------- |
|---|
| 177 | ! sea salt specific |
|---|
| 178 | !----------------------------------------------------------------------- |
|---|
| 179 | ! REAL*8, DIMENSION(nmx) :: ra, rb |
|---|
| 180 | ! REAL*8 :: ch_ss(nmx,12) |
|---|
| 181 | |
|---|
| 182 | !----------------------------------------------------------------------- |
|---|
| 183 | ! emissions (input) |
|---|
| 184 | !----------------------------------------------------------------------- |
|---|
| 185 | REAL*8 :: e_an(imx,jmx,2,nmx), e_bb(imx,jmx,nmx), & |
|---|
| 186 | e_ac(imx,jmx,lmx,nmx) |
|---|
| 187 | |
|---|
| 188 | !----------------------------------------------------------------------- |
|---|
| 189 | ! diagnostics (budget) |
|---|
| 190 | !----------------------------------------------------------------------- |
|---|
| 191 | ! ! tendencies per time step and process |
|---|
| 192 | ! REAL*8, TARGET :: bems(imx,jmx,nmx), bdry(imx,jmx,nmx), bstl(imx,jmx,nmx) |
|---|
| 193 | ! REAL*8, TARGET :: bwet(imx,jmx,nmx), bcnv(imx,jmx,nmx)! |
|---|
| 194 | |
|---|
| 195 | ! ! integrated tendencies per process |
|---|
| 196 | ! REAL*8, TARGET :: tems(imx,jmx,nmx), tstl(imx,jmx,nmx) |
|---|
| 197 | ! REAL*8, TARGET :: tdry(imx,jmx,nmx), twet(imx,jmx,nmx), tcnv(imx,jmx,nmx) |
|---|
| 198 | |
|---|
| 199 | ! global mass balance per time step |
|---|
| 200 | REAL*8 :: tmas0(nmx), tmas1(nmx) |
|---|
| 201 | REAL*8 :: dtems(nmx), dttrp(nmx), dtdif(nmx), dtcnv(nmx) |
|---|
| 202 | REAL*8 :: dtwet(nmx), dtdry(nmx), dtstl(nmx) |
|---|
| 203 | REAL*8 :: dtems2(nmx), dttrp2(nmx), dtdif2(nmx), dtcnv2(nmx) |
|---|
| 204 | REAL*8 :: dtwet2(nmx), dtdry2(nmx), dtstl2(nmx) |
|---|
| 205 | |
|---|
| 206 | ! detailed integrated budgets for individual emissions |
|---|
| 207 | REAL*8, TARGET :: ems_an(imx,jmx,nmx), ems_bb(imx,jmx,nmx), ems_tp(imx,jmx) |
|---|
| 208 | REAL*8, TARGET :: ems_ac(imx,jmx,lmx,nmx) |
|---|
| 209 | REAL*8, TARGET :: ems_co(imx,jmx,nmx) |
|---|
| 210 | |
|---|
| 211 | |
|---|
| 212 | ! executable statements |
|---|
| 213 | |
|---|
| 214 | DO n = 1,nmx |
|---|
| 215 | ! if(ipr.eq.1)write(0,*)'in seasalt',n,ipr,ilwi |
|---|
| 216 | bems(:,:,n) = 0.0 |
|---|
| 217 | rho_d = den_seas(n) |
|---|
| 218 | r0 = ra(n)*frh |
|---|
| 219 | r1 = rb(n)*frh |
|---|
| 220 | r = r0 |
|---|
| 221 | nr = INT((r1-r0)/dr+.001) |
|---|
| 222 | ! if(ipr.eq.1.and.n.eq.1)write(0,*)'in seasalt',nr,r1,r0,dr,rho_d |
|---|
| 223 | DO ir = 1,nr |
|---|
| 224 | r_w = r + dr*0.5 |
|---|
| 225 | r = r + dr |
|---|
| 226 | IF (new) THEN |
|---|
| 227 | a = 4.7*(1.0 + theta*r_w)**(-0.017*r_w**(-1.44)) |
|---|
| 228 | c0 = c_new |
|---|
| 229 | b0 = b_new |
|---|
| 230 | ELSE |
|---|
| 231 | a = 3.0 |
|---|
| 232 | c0 = c_old |
|---|
| 233 | b0 = b_old |
|---|
| 234 | END IF |
|---|
| 235 | ! |
|---|
| 236 | b = (b0(1) - LOG10(r_w))/b0(2) |
|---|
| 237 | dfn = (c0(1)/r_w**a)*(1.0 + c0(3)*r_w**c0(4))* & |
|---|
| 238 | 10**(c0(5)*EXP(-(b**2))) |
|---|
| 239 | |
|---|
| 240 | r_d = r_w/frh*1.0D-6 ! um -> m |
|---|
| 241 | dfm = 4.0/3.0*pi*r_d**3*rho_d*frh*dfn*dr*dt1 |
|---|
| 242 | DO i = 1,imx |
|---|
| 243 | DO j = 1,jmx |
|---|
| 244 | ! IF (water(i,j) > 0.0) THEN |
|---|
| 245 | IF (ilwi(i,j) == 0) THEN |
|---|
| 246 | ! src = dfm*dxy(j)*water(i,j)*w10m(i,j)**c0(2) |
|---|
| 247 | src = dfm*dxy(j)*w10m(i,j)**c0(2) |
|---|
| 248 | ! src = ch_ss(n,dt(1)%mn)*dfm*dxy(j)*w10m(i,j)**c0(2) |
|---|
| 249 | tc(i,j,1,n) = tc(i,j,1,n) + src/airmas(i,j,1) |
|---|
| 250 | ! if(ipr.eq.1)write(0,*)n,dfm,c0(2),dxy(j),w10m(i,j),src,airmas(i,j,1) |
|---|
| 251 | ELSE |
|---|
| 252 | src = 0.0 |
|---|
| 253 | END IF |
|---|
| 254 | bems(i,j,n) = bems(i,j,n) + src |
|---|
| 255 | END DO ! i |
|---|
| 256 | END DO ! j |
|---|
| 257 | END DO ! ir |
|---|
| 258 | END DO ! n |
|---|
| 259 | |
|---|
| 260 | END SUBROUTINE source_ss |
|---|
| 261 | END MODULE GOCART_SEASALT |
|---|
