chem/module_gocart_seasalt.F

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