phys/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,                  &
   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't 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't 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