Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/wetdepo.f90

    r13 r4  
    3939  ! Code may not be correct for decay of deposition!                           *
    4040  !                                                                            *
    41   ! Modification by Sabine Eckhart to introduce a new in-/below-cloud
    42   ! scheme, not dated
    43   ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
    4441  !*****************************************************************************
    4542  !                                                                            *
     
    7471
    7572  integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
    76   integer :: ngrid,itage,nage,hz,il,interp_time, n , clouds_v !NIK scheme
    77   integer :: kz !PS scheme
    78   integer :: ks, kp, n1,n2, icbot,ictop, indcloud
    79   integer :: scheme_number ! NIK==1, PS ==2
     73  integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
     74  integer :: ks, kp
    8075  real :: S_i, act_temp, cl, cle ! in cloud scavenging
    8176  real :: clouds_h ! cloud height for the specific grid point
    82   real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav, precsub,f
     77  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav
    8378  real :: wetdeposit(maxspec),restmass
    8479  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
     
    151146    interp_time=nint(itime-0.5*ltsample)
    152147
    153 !    PS nest case still needs to be implemented!!   
    154 !    if (ngrid.eq.0) then
    155 !      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
    156 !           1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
    157 !           memtime(1),memtime(2),interp_time,lsp,convp,cc)
    158            call interpol_rain(lsprec,convprec,tcc,     &
    159             icloudbot,icloudthck,nxmax,nymax,1,nx,ny,         &
    160             memind,sngl(xtra1(jpart)),sngl(ytra1(jpart)),1,memtime(1), &
    161             memtime(2),interp_time,lsp,convp,cc,icbot,ictop,icmv) 
    162 !    else
    163 !      call interpol_rain_nests(lsprecn,convprecn,tccn, &
    164 !           nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
    165 !           memtime(1),memtime(2),interp_time,lsp,convp,cc)
    166 !    endif
    167 
    168 
    169  !   if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
    170  !PS 2012: subtract a small value, eg 0.01 mm/h, to remove spurious precip
    171         prec = lsp+convp
    172         precsub = 0.01
    173         if (prec .lt. precsub) then
    174           goto 20
    175         else
    176           f = (prec-precsub)/prec
    177           lsp = f*lsp
    178           convp = f*convp
    179         endif
    180 
     148    if (ngrid.eq.0) then
     149      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
     150           1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
     151           memtime(1),memtime(2),interp_time,lsp,convp,cc)
     152    else
     153      call interpol_rain_nests(lsprecn,convprecn,tccn, &
     154           nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
     155           memtime(1),memtime(2),interp_time,lsp,convp,cc)
     156    endif
     157
     158    if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
    181159
    182160  ! get the level were the actual particle is in
    183161      do il=2,nz
    184162        if (height(il).gt.ztra1(jpart)) then
    185           !hz=il-1
    186           kz=il-1
     163          hz=il-1
    187164          goto 26
    188165        endif
     
    197174  ! scavenging is done
    198175
    199 !old scheme
    200 !  if (ngrid.eq.0) then
    201 !     clouds_v=clouds(ix,jy,hz,n)
    202 !     clouds_h=cloudsh(ix,jy,n)
    203 !  else
    204 !     clouds_v=cloudsn(ix,jy,hz,n,ngrid)
    205 !     clouds_h=cloudsnh(ix,jy,n,ngrid)
    206 !  endif
    207 !  !write(*,*) 'there is
    208 !  !    + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
    209 !  if (clouds_v.le.1) goto 20
    210 !  !write (*,*) 'there is scavenging'
    211 
    212   ! PS: part of 2011/2012 fix
    213 
    214         if (ztra1(jpart) .le. float(ictop)) then
    215           if (ztra1(jpart) .gt. float(icbot)) then
    216             indcloud = 2 ! in-cloud
    217           else
    218             indcloud = 1 ! below-cloud
    219           endif
    220         elseif (ictop .eq. icmv) then
    221           indcloud = 0 ! no cloud found, use old scheme
    222         else
    223           goto 20 ! above cloud
    224         endif
    225 
     176  if (ngrid.eq.0) then
     177     clouds_v=clouds(ix,jy,hz,n)
     178     clouds_h=cloudsh(ix,jy,n)
     179  else
     180     clouds_v=cloudsn(ix,jy,hz,n,ngrid)
     181     clouds_h=cloudsnh(ix,jy,n,ngrid)
     182  endif
     183  !write(*,*) 'there is
     184  !    + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
     185  if (clouds_v.le.1) goto 20
     186  !write (*,*) 'there is scavenging'
    226187
    227188  ! 1) Parameterization of the the area fraction of the grid cell where the
     
    267228  !**********************************************************
    268229
    269     do ks=1,nspec            ! loop over species
     230    do ks=1,nspec                                  ! loop over species
    270231      wetdeposit(ks)=0.
    271 
    272    
    273     !conflicting changes to the same routine: 1=NIK 2 =PS
    274     scheme_number=2   
    275     if (scheme_number.eq.1) then !NIK
    276 
    277232      if (weta(ks).gt.0.) then
    278233        if (clouds_v.ge.4) then
     
    281236          wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
    282237  !        write(*,*) 'bel. wetscav: ',wetscav
    283 
    284238        else ! below_cloud clouds_v is lt 4 and gt 1 -> in cloud scavenging
    285239  !        IN CLOUD SCAVENGING
     
    291245             act_temp=tt(ix,jy,hz,n)
    292246          endif
    293 
    294 ! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening
    295 ! weta_in=2.0E-07 (default)
    296 ! wetb_in=0.36 (default)
    297 ! wetc_in=0.9 (default)
    298 ! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0-no scaling)
    299           cl=weta_in(ks)*prec**wetb_in(ks)
     247          cl=2E-7*prec**0.36
    300248          if (dquer(ks).gt.0) then ! is particle
    301             S_i=wetc_in(ks)/cl
     249            S_i=0.9/cl
    302250           else ! is gas
    303251            cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
    304252            S_i=1/cle
    305253           endif
    306            wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
     254           wetscav=S_i*prec/3.6E6/clouds_h
    307255  !         write(*,*) 'in. wetscav:'
    308256  !    +          ,wetscav,cle,cl,act_temp,prec,clouds_h
     
    341289         wetdeposit(ks)=0.
    342290      endif ! weta(k)
    343 
    344      elseif (scheme_number.eq.2) then ! PS
    345 
    346 !PS          indcloud=0 ! Use this for FOR TESTING,
    347 !PS          will skip the new in/below cloud method !!!             
    348 
    349           if (weta(ks).gt.0.) then
    350             if (indcloud .eq. 1) then ! BELOW CLOUD SCAVENGING
    351 !C               for aerosols and not highliy soluble substances weta=5E-6
    352               wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
    353 !c             write(*,*) 'bel. wetscav: ',wetscav
    354             elseif (indcloud .eq. 2) then !  IN CLOUD SCAVENGING
    355               if (ngrid.gt.0) then
    356                  act_temp=ttn(ix,jy,kz,n,ngrid)
    357               else
    358                  act_temp=tt(ix,jy,kz,n)
    359               endif
    360 
    361 ! from NIK             
    362 ! weta_in=2.0E-07 (default)
    363 ! wetb_in=0.36 (default)
    364 ! wetc_in=0.9 (default)
    365 
    366 
    367               cl=2E-7*prec**0.36
    368               if (dquer(ks).gt.0) then ! is particle
    369                 S_i=0.9/cl
    370               else ! is gas
    371                 cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
    372                 S_i=1/cle
    373               endif
    374               wetscav=S_i*prec/3.6E6/(ictop-icbot) ! 3.6e6 converts mm/h to m/s
    375             else ! PS: no cloud diagnosed, old scheme,
    376 !CPS          using with fixed a,b for simplicity, one may wish to change!!
    377               wetscav = 1.e-4*prec**0.62
    378             endif
    379 
    380 
    381             wetdeposit(ks)=xmass1(jpart,ks)*    &
    382             ! (1.-exp(-wetscav*abs(ltsample)))*fraction  ! wet deposition
    383              (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! fraction = grfraction (PS)
    384             restmass = xmass1(jpart,ks)-wetdeposit(ks)
    385             if (ioutputforeachrelease.eq.1) then
    386               kp=npoint(jpart)
    387             else
    388               kp=1
    389             endif
    390             if (restmass .gt. smallnum) then
    391               xmass1(jpart,ks)=restmass
    392 !cccccccccccccccc depostatistic
    393 !c            wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
    394 !cccccccccccccccc depostatistic
    395             else
    396               xmass1(jpart,ks)=0.
    397             endif
    398 !C Correct deposited mass to the last time step when radioactive decay of
    399 !C gridded deposited mass was calculated
    400             if (decay(ks).gt.0.) then
    401               wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
    402             endif
    403           else  ! weta(k)<0
    404              wetdeposit(ks)=0.
    405           endif
    406 
    407      endif !on scheme
    408 
    409 
    410 
    411291    end do
    412292
     
    415295  !*****************************************************************************
    416296
    417 !    if (ldirect.eq.1) then
    418 !    call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
    419 !         real(ytra1(jpart)),nage,kp)
    420 !    if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
    421 !         wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), &
    422 !         nage,kp)
    423 !    endif
    424 
    425    !PS
    426         if (ldirect.eq.1) then
    427           call wetdepokernel(nclass(jpart),wetdeposit,     &
    428              sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp)
    429           if (nested_output.eq.1)                          &
    430            call wetdepokernel_nest(nclass(jpart),wetdeposit, &
    431              sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp) 
    432         endif
    433 
     297    if (ldirect.eq.1) then
     298    call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
     299         real(ytra1(jpart)),nage,kp)
     300    if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
     301         wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), &
     302         nage,kp)
     303    endif
    434304
    43530520  continue
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG