chem/module_cam_mam_dust_sediment.F

   1 ! module_cam_mam_dust_sediment.F
   2 ! adapted from cam3 dust_sediment_mod.F90 by r.c.easter, november 2010
   3 !
   4 ! 2010-xx-xx rce notes
   5 !
   6 !     > *** THIS IS NOT YET IMPLEMENTED ***
   7 !
   8 !     > all calls to history routines (outfld) deactivated
   9 !--------------------------------------------------------------
  10
  11 module dust_sediment_mod
  12
  13 !---------------------------------------------------------------------------------
  14 ! Purpose:
  15 !
  16 ! Contains routines to compute tendencies from sedimentation of dust
  17 !
  18 ! Author: Phil Rasch
  19 !
  20 !---------------------------------------------------------------------------------
  21
  22   use shr_kind_mod,  only: r8=>shr_kind_r8
  23 ! use ppgrid,        only: pcols, pver, pverp
  24 ! use physconst,     only: gravit, rair
  25 ! use cam_logfile,   only: iulog
  26
  27   private
  28   public :: dust_sediment_vel, dust_sediment_tend
  29
  30
  31   real (r8), parameter :: vland  = 2.8_r8            ! dust fall velocity over land  (cm/s)
  32   real (r8), parameter :: vocean = 1.5_r8            ! dust fall velocity over ocean (cm/s)
  33   real (r8), parameter :: mxsedfac   = 0.99_r8       ! maximum sedimentation flux factor
  34
  35 contains
  36
  37 !===============================================================================
  38   subroutine dust_sediment_vel ( ncol )
  39   implicit none
  40   integer :: ncol
  41   return
  42 end subroutine dust_sediment_vel !Balwinder.Singh@pnnl.gov - Added "end subroutine" as Intel compiler was complainig missing "end subroutine"
  43
  44
  45 !===============================================================================
  46   subroutine dust_sediment_tend ( ncol )
  47   implicit none
  48   integer :: ncol
  49   return
  50 end subroutine dust_sediment_tend!Balwinder.Singh@pnnl.gov - Added "end subroutine" as Intel compiler was complainig missing "end subroutine"
  51
  52
  53 #undef        ABC123_THIS_CODE_ACTIVE
  54 #if ( defined ABC123_THIS_CODE_ACTIVE )
  55
  56 !===============================================================================
  57   subroutine dust_sediment_vel (ncol,                               &
  58        icefrac , landfrac, ocnfrac , pmid    , pdel    , t       , &
  59        dustmr  , pvdust   )
  60
  61 !----------------------------------------------------------------------
  62
  63 ! Compute gravitational sedimentation velocities for dust
  64
  65     implicit none
  66
  67 ! Arguments
  68     integer, intent(in) :: ncol                     ! number of colums to process
  69
  70     real(r8), intent(in)  :: icefrac (pcols)        ! sea ice fraction (fraction)
  71     real(r8), intent(in)  :: landfrac(pcols)        ! land fraction (fraction)
  72     real(r8), intent(in)  :: ocnfrac (pcols)        ! ocean fraction (fraction)
  73     real(r8), intent(in)  :: pmid  (pcols,pver)     ! pressure of midpoint levels (Pa)
  74     real(r8), intent(in)  :: pdel  (pcols,pver)     ! pressure diff across layer (Pa)
  75     real(r8), intent(in)  :: t     (pcols,pver)     ! temperature (K)
  76     real(r8), intent(in)  :: dustmr(pcols,pver)     ! dust (kg/kg)
  77
  78     real(r8), intent(out) :: pvdust (pcols,pverp)    ! vertical velocity of dust (Pa/s)
  79 ! -> note that pvel is at the interfaces (loss from cell is based on pvel(k+1))
  80
  81 ! Local variables
  82     real (r8) :: rho(pcols,pver)                    ! air density in kg/m3
  83     real (r8) :: vfall(pcols)                       ! settling velocity of dust particles (m/s)
  84
  85     integer i,k
  86
  87     real (r8) :: lbound, ac, bc, cc
  88
  89 !-----------------------------------------------------------------------
  90 !--------------------- dust fall velocity ----------------------------
  91 !-----------------------------------------------------------------------
  92
  93     do k = 1,pver
  94        do i = 1,ncol
  95
  96           ! merge the dust fall velocities for land and ocean (cm/s)
  97           ! SHOULD ALSO ACCOUNT FOR ICEFRAC
  98           vfall(i) = vland*landfrac(i) + vocean*(1._r8-landfrac(i))
  99           !!         vfall(i) = vland*landfrac(i) + vocean*ocnfrac(i) + vseaice*icefrac(i)
 100
 101           ! fall velocity (assume positive downward)
 102           pvdust(i,k+1) = vfall(i)
 103        end do
 104     end do
 105
 106     return
 107   end subroutine dust_sediment_vel
 108
 109
 110 !===============================================================================
 111   subroutine dust_sediment_tend ( &
 112        ncol,   dtime,  pint,     pmid,    pdel,  t,   &
 113        dustmr ,pvdust, dusttend, sfdust )
 114
 115 !----------------------------------------------------------------------
 116 !     Apply Particle Gravitational Sedimentation
 117 !----------------------------------------------------------------------
 118
 119     implicit none
 120
 121 ! Arguments
 122     integer,  intent(in)  :: ncol                      ! number of colums to process
 123
 124     real(r8), intent(in)  :: dtime                     ! time step
 125     real(r8), intent(in)  :: pint  (pcols,pverp)       ! interfaces pressure (Pa)
 126     real(r8), intent(in)  :: pmid  (pcols,pver)        ! midpoint pressures (Pa)
 127     real(r8), intent(in)  :: pdel  (pcols,pver)        ! pressure diff across layer (Pa)
 128     real(r8), intent(in)  :: t     (pcols,pver)        ! temperature (K)
 129     real(r8), intent(in)  :: dustmr(pcols,pver)        ! dust (kg/kg)
 130     real(r8), intent(in)  :: pvdust (pcols,pverp)      ! vertical velocity of dust drops  (Pa/s)
 131 ! -> note that pvel is at the interfaces (loss from cell is based on pvel(k+1))
 132
 133     real(r8), intent(out) :: dusttend(pcols,pver)      ! dust tend
 134     real(r8), intent(out) :: sfdust  (pcols)           ! surface flux of dust (rain, kg/m/s)
 135
 136 ! Local variables
 137     real(r8) :: fxdust(pcols,pverp)                     ! fluxes at the interfaces, dust (positive = down)
 138
 139     integer :: i,k
 140 !----------------------------------------------------------------------
 141
 142 ! initialize variables
 143     fxdust  (:ncol,:) = 0._r8 ! flux at interfaces (dust)
 144     dusttend(:ncol,:) = 0._r8 ! tend (dust)
 145     sfdust(:ncol)     = 0._r8 ! sedimentation flux out bot of column (dust)
 146
 147 ! fluxes at interior points
 148     call getflx(ncol, pint, dustmr, pvdust, dtime, fxdust)
 149
 150 ! calculate fluxes at boundaries
 151     do i = 1,ncol
 152        fxdust(i,1) = 0
 153 ! surface flux by upstream scheme
 154        fxdust(i,pverp) = dustmr(i,pver) * pvdust(i,pverp) * dtime
 155     end do
 156
 157 ! filter out any negative fluxes from the getflx routine
 158     do k = 2,pver
 159        fxdust(:ncol,k) = max(0._r8, fxdust(:ncol,k))
 160     end do
 161
 162 ! Limit the flux out of the bottom of each cell to the water content in each phase.
 163 ! Apply mxsedfac to prevent generating very small negative cloud water/ice
 164 ! NOTE, REMOVED CLOUD FACTOR FROM AVAILABLE WATER. ALL CLOUD WATER IS IN CLOUDS.
 165 ! ***Should we include the flux in the top, to allow for thin surface layers?
 166 ! ***Requires simple treatment of cloud overlap, already included below.
 167     do k = 1,pver
 168        do i = 1,ncol
 169           fxdust(i,k+1) = min( fxdust(i,k+1), mxsedfac * dustmr(i,k) * pdel(i,k) )
 170 !!$        fxdust(i,k+1) = min( fxdust(i,k+1), dustmr(i,k) * pdel(i,k) + fxdust(i,k))
 171        end do
 172     end do
 173
 174 ! Now calculate the tendencies
 175     do k = 1,pver
 176        do i = 1,ncol
 177 ! net flux into cloud changes cloud dust/ice (all flux is out of cloud)
 178           dusttend(i,k)  = (fxdust(i,k) - fxdust(i,k+1)) / (dtime * pdel(i,k))
 179        end do
 180     end do
 181
 182 ! convert flux out the bottom to mass units Pa -> kg/m2/s
 183     sfdust(:ncol) = fxdust(:ncol,pverp) / (dtime*gravit)
 184
 185     return
 186   end subroutine dust_sediment_tend
 187
 188 !===============================================================================
 189   subroutine getflx(ncol, xw, phi, vel, deltat, flux)
 190
 191 !.....xw1.......xw2.......xw3.......xw4.......xw5.......xw6
 192 !....psiw1.....psiw2.....psiw3.....psiw4.....psiw5.....psiw6
 193 !....velw1.....velw2.....velw3.....velw4.....velw5.....velw6
 194 !.........phi1......phi2.......phi3.....phi4.......phi5.......
 195
 196
 197     implicit none
 198
 199     integer ncol                      ! number of colums to process
 200
 201     integer i
 202     integer k
 203
 204     real (r8) vel(pcols,pverp)
 205     real (r8) flux(pcols,pverp)
 206     real (r8) xw(pcols,pverp)
 207     real (r8) psi(pcols,pverp)
 208     real (r8) phi(pcols,pverp-1)
 209     real (r8) fdot(pcols,pverp)
 210     real (r8) xx(pcols)
 211     real (r8) fxdot(pcols)
 212     real (r8) fxdd(pcols)
 213
 214     real (r8) psistar(pcols)
 215     real (r8) deltat
 216
 217     real (r8) xxk(pcols,pver)
 218
 219     do i = 1,ncol
 220 !        integral of phi
 221        psi(i,1) = 0._r8
 222 !        fluxes at boundaries
 223        flux(i,1) = 0
 224        flux(i,pverp) = 0._r8
 225     end do
 226
 227 !     integral function
 228     do k = 2,pverp
 229        do i = 1,ncol
 230           psi(i,k) = phi(i,k-1)*(xw(i,k)-xw(i,k-1)) + psi(i,k-1)
 231        end do
 232     end do
 233
 234
 235 !     calculate the derivatives for the interpolating polynomial
 236     call cfdotmc_pro (ncol, xw, psi, fdot)
 237
 238 !  NEW WAY
 239 !     calculate fluxes at interior pts
 240     do k = 2,pver
 241        do i = 1,ncol
 242           xxk(i,k) = xw(i,k)-vel(i,k)*deltat
 243        end do
 244     end do
 245     do k = 2,pver
 246        call cfint2(ncol, xw, psi, fdot, xxk(1,k), fxdot, fxdd, psistar)
 247        do i = 1,ncol
 248           flux(i,k) = (psi(i,k)-psistar(i))
 249        end do
 250     end do
 251
 252
 253     return
 254   end subroutine getflx
 255
 256
 257
 258 !##############################################################################
 259
 260   subroutine cfint2 (ncol, x, f, fdot, xin, fxdot, fxdd, psistar)
 261
 262
 263     implicit none
 264
 265 ! input
 266     integer ncol                      ! number of colums to process
 267
 268     real (r8) x(pcols, pverp)
 269     real (r8) f(pcols, pverp)
 270     real (r8) fdot(pcols, pverp)
 271     real (r8) xin(pcols)
 272
 273 ! output
 274     real (r8) fxdot(pcols)
 275     real (r8) fxdd(pcols)
 276     real (r8) psistar(pcols)
 277
 278     integer i
 279     integer k
 280     integer intz(pcols)
 281     real (r8) dx
 282     real (r8) s
 283     real (r8) c2
 284     real (r8) c3
 285     real (r8) xx
 286     real (r8) xinf
 287     real (r8) psi1, psi2, psi3, psim
 288     real (r8) cfint
 289     real (r8) cfnew
 290     real (r8) xins(pcols)
 291
 292 !     the minmod function
 293     real (r8) a, b, c
 294     real (r8) minmod
 295     real (r8) medan
 296     minmod(a,b) = 0.5_r8*(sign(1._r8,a) + sign(1._r8,b))*min(abs(a),abs(b))
 297     medan(a,b,c) = a + minmod(b-a,c-a)
 298
 299     do i = 1,ncol
 300        xins(i) = medan(x(i,1), xin(i), x(i,pverp))
 301        intz(i) = 0
 302     end do
 303
 304 ! first find the interval
 305     do k =  1,pverp-1
 306        do i = 1,ncol
 307           if ((xins(i)-x(i,k))*(x(i,k+1)-xins(i)).ge.0) then
 308              intz(i) = k
 309           endif
 310        end do
 311     end do
 312
 313     do i = 1,ncol
 314        if (intz(i).eq.0._r8) then
 315           write(iulog,*) ' interval was not found for col i ', i
 316           stop
 317        endif
 318     end do
 319
 320 ! now interpolate
 321     do i = 1,ncol
 322        k = intz(i)
 323        dx = (x(i,k+1)-x(i,k))
 324        s = (f(i,k+1)-f(i,k))/dx
 325        c2 = (3*s-2*fdot(i,k)-fdot(i,k+1))/dx
 326        c3 = (fdot(i,k)+fdot(i,k+1)-2*s)/dx**2
 327        xx = (xins(i)-x(i,k))
 328        fxdot(i) =  (3*c3*xx + 2*c2)*xx + fdot(i,k)
 329        fxdd(i) = 6*c3*xx + 2*c2
 330        cfint = ((c3*xx + c2)*xx + fdot(i,k))*xx + f(i,k)
 331
 332 !        limit the interpolant
 333        psi1 = f(i,k)+(f(i,k+1)-f(i,k))*xx/dx
 334        if (k.eq.1) then
 335           psi2 = f(i,1)
 336        else
 337           psi2 = f(i,k) + (f(i,k)-f(i,k-1))*xx/(x(i,k)-x(i,k-1))
 338        endif
 339        if (k+1.eq.pverp) then
 340           psi3 = f(i,pverp)
 341        else
 342           psi3 = f(i,k+1) - (f(i,k+2)-f(i,k+1))*(dx-xx)/(x(i,k+2)-x(i,k+1))
 343        endif
 344        psim = medan(psi1, psi2, psi3)
 345        cfnew = medan(cfint, psi1, psim)
 346        if (abs(cfnew-cfint)/(abs(cfnew)+abs(cfint)+1.e-36_r8)  .gt..03_r8) then
 347 !     CHANGE THIS BACK LATER!!!
 348 !     $        .gt..1) then
 349
 350
 351 !     UNCOMMENT THIS LATER!!!
 352 !            write(iulog,*) ' cfint2 limiting important ', cfint, cfnew
 353
 354
 355        endif
 356        psistar(i) = cfnew
 357     end do
 358
 359     return
 360   end subroutine cfint2
 361
 362
 363
 364 !##############################################################################
 365
 366   subroutine cfdotmc_pro (ncol, x, f, fdot)
 367
 368 !     prototype version; eventually replace with final SPITFIRE scheme
 369
 370 !     calculate the derivative for the interpolating polynomial
 371 !     multi column version
 372
 373
 374     implicit none
 375
 376 ! input
 377     integer ncol                      ! number of colums to process
 378
 379     real (r8) x(pcols, pverp)
 380     real (r8) f(pcols, pverp)
 381 ! output
 382     real (r8) fdot(pcols, pverp)          ! derivative at nodes
 383
 384 ! assumed variable distribution
 385 !     x1.......x2.......x3.......x4.......x5.......x6     1,pverp points
 386 !     f1.......f2.......f3.......f4.......f5.......f6     1,pverp points
 387 !     ...sh1.......sh2......sh3......sh4......sh5....     1,pver points
 388 !     .........d2.......d3.......d4.......d5.........     2,pver points
 389 !     .........s2.......s3.......s4.......s5.........     2,pver points
 390 !     .............dh2......dh3......dh4.............     2,pver-1 points
 391 !     .............eh2......eh3......eh4.............     2,pver-1 points
 392 !     ..................e3.......e4..................     3,pver-1 points
 393 !     .................ppl3......ppl4................     3,pver-1 points
 394 !     .................ppr3......ppr4................     3,pver-1 points
 395 !     .................t3........t4..................     3,pver-1 points
 396 !     ................fdot3.....fdot4................     3,pver-1 points
 397
 398
 399 ! work variables
 400
 401
 402     integer i
 403     integer k
 404
 405     real (r8) a                    ! work var
 406     real (r8) b                    ! work var
 407     real (r8) c                    ! work var
 408     real (r8) s(pcols,pverp)             ! first divided differences at nodes
 409     real (r8) sh(pcols,pverp)            ! first divided differences between nodes
 410     real (r8) d(pcols,pverp)             ! second divided differences at nodes
 411     real (r8) dh(pcols,pverp)            ! second divided differences between nodes
 412     real (r8) e(pcols,pverp)             ! third divided differences at nodes
 413     real (r8) eh(pcols,pverp)            ! third divided differences between nodes
 414     real (r8) pp                   ! p prime
 415     real (r8) ppl(pcols,pverp)           ! p prime on left
 416     real (r8) ppr(pcols,pverp)           ! p prime on right
 417     real (r8) qpl
 418     real (r8) qpr
 419     real (r8) ttt
 420     real (r8) t
 421     real (r8) tmin
 422     real (r8) tmax
 423     real (r8) delxh(pcols,pverp)
 424
 425
 426 !     the minmod function
 427     real (r8) minmod
 428     real (r8) medan
 429     minmod(a,b) = 0.5_r8*(sign(1._r8,a) + sign(1._r8,b))*min(abs(a),abs(b))
 430     medan(a,b,c) = a + minmod(b-a,c-a)
 431
 432     do k = 1,pver
 433
 434
 435 !        first divided differences between nodes
 436        do i = 1, ncol
 437           delxh(i,k) = (x(i,k+1)-x(i,k))
 438           sh(i,k) = (f(i,k+1)-f(i,k))/delxh(i,k)
 439        end do
 440
 441 !        first and second divided differences at nodes
 442        if (k.ge.2) then
 443           do i = 1,ncol
 444              d(i,k) = (sh(i,k)-sh(i,k-1))/(x(i,k+1)-x(i,k-1))
 445              s(i,k) = minmod(sh(i,k),sh(i,k-1))
 446           end do
 447        endif
 448     end do
 449
 450 !     second and third divided diffs between nodes
 451     do k = 2,pver-1
 452        do i = 1, ncol
 453           eh(i,k) = (d(i,k+1)-d(i,k))/(x(i,k+2)-x(i,k-1))
 454           dh(i,k) = minmod(d(i,k),d(i,k+1))
 455        end do
 456     end do
 457
 458 !     treat the boundaries
 459     do i = 1,ncol
 460        e(i,2) = eh(i,2)
 461        e(i,pver) = eh(i,pver-1)
 462 !        outside level
 463        fdot(i,1) = sh(i,1) - d(i,2)*delxh(i,1)  &
 464             - eh(i,2)*delxh(i,1)*(x(i,1)-x(i,3))
 465        fdot(i,1) = minmod(fdot(i,1),3*sh(i,1))
 466        fdot(i,pverp) = sh(i,pver) + d(i,pver)*delxh(i,pver)  &
 467             + eh(i,pver-1)*delxh(i,pver)*(x(i,pverp)-x(i,pver-1))
 468        fdot(i,pverp) = minmod(fdot(i,pverp),3*sh(i,pver))
 469 !        one in from boundary
 470        fdot(i,2) = sh(i,1) + d(i,2)*delxh(i,1) - eh(i,2)*delxh(i,1)*delxh(i,2)
 471        fdot(i,2) = minmod(fdot(i,2),3*s(i,2))
 472        fdot(i,pver) = sh(i,pver) - d(i,pver)*delxh(i,pver)   &
 473             - eh(i,pver-1)*delxh(i,pver)*delxh(i,pver-1)
 474        fdot(i,pver) = minmod(fdot(i,pver),3*s(i,pver))
 475     end do
 476
 477
 478     do k = 3,pver-1
 479        do i = 1,ncol
 480           e(i,k) = minmod(eh(i,k),eh(i,k-1))
 481        end do
 482     end do
 483
 484
 485
 486     do k = 3,pver-1
 487
 488        do i = 1,ncol
 489
 490 !           p prime at k-0.5
 491           ppl(i,k)=sh(i,k-1) + dh(i,k-1)*delxh(i,k-1)
 492 !           p prime at k+0.5
 493           ppr(i,k)=sh(i,k)   - dh(i,k)  *delxh(i,k)
 494
 495           t = minmod(ppl(i,k),ppr(i,k))
 496
 497 !           derivate from parabola thru f(i,k-1), f(i,k), and f(i,k+1)
 498           pp = sh(i,k-1) + d(i,k)*delxh(i,k-1)
 499
 500 !           quartic estimate of fdot
 501           fdot(i,k) = pp                            &
 502                - delxh(i,k-1)*delxh(i,k)            &
 503                *(  eh(i,k-1)*(x(i,k+2)-x(i,k  ))    &
 504                + eh(i,k  )*(x(i,k  )-x(i,k-2))      &
 505                )/(x(i,k+2)-x(i,k-2))
 506
 507 !           now limit it
 508           qpl = sh(i,k-1)       &
 509                + delxh(i,k-1)*minmod(d(i,k-1)+e(i,k-1)*(x(i,k)-x(i,k-2)), &
 510                d(i,k)  -e(i,k)*delxh(i,k))
 511           qpr = sh(i,k)         &
 512                + delxh(i,k  )*minmod(d(i,k)  +e(i,k)*delxh(i,k-1),        &
 513                d(i,k+1)+e(i,k+1)*(x(i,k)-x(i,k+2)))
 514
 515           fdot(i,k) = medan(fdot(i,k), qpl, qpr)
 516
 517           ttt = minmod(qpl, qpr)
 518           tmin = min(0._r8,3*s(i,k),1.5_r8*t,ttt)
 519           tmax = max(0._r8,3*s(i,k),1.5_r8*t,ttt)
 520
 521           fdot(i,k) = fdot(i,k) + minmod(tmin-fdot(i,k), tmax-fdot(i,k))
 522
 523        end do
 524
 525     end do
 526
 527     return
 528   end subroutine cfdotmc_pro
 529 #endif
 530 end module dust_sediment_mod