Changeset 24 for trunk/src/wetdepo.f90


Ignore:
Timestamp:
May 23, 2014, 11:48:41 AM (10 years ago)
Author:
igpis
Message:

version 9.2 beta. Changes from HH, AST, MC, NIK, IP. Changes in vert transform. New SPECIES input includes scavenging coefficients

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/wetdepo.f90

    r20 r24  
    2121
    2222subroutine wetdepo(itime,ltsample,loutnext)
    23   !                     i      i        i
     23  !                  i      i        i
    2424  !*****************************************************************************
    2525  !                                                                            *
    … …  
    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
    … …  
    10398
    10499  do jpart=1,numpart
     100
    105101    if (itra1(jpart).eq.-999999999) goto 20
    106102    if(ldirect.eq.1)then
    … …  
    114110      if (itage.lt.lage(nage)) goto 33
    115111    end do
    116 33   continue
     11233  continue
    117113
    118114
    … …  
    123119    do j=numbnests,1,-1
    124120      if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. &
    125            (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
     121      (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
    126122        ngrid=j
    127123        goto 23
    128124      endif
    129125    end do
    130 23   continue
     12623  continue
    131127
    132128
    … …  
    151147    interp_time=nint(itime-0.5*ltsample)
    152148
    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 
    181 
     149    if (ngrid.eq.0) then
     150      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
     151      1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
     152      memtime(1),memtime(2),interp_time,lsp,convp,cc)
     153    else
     154      call interpol_rain_nests(lsprecn,convprecn,tccn, &
     155      nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
     156      memtime(1),memtime(2),interp_time,lsp,convp,cc)
     157    endif
     158
     159    if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
    182160  ! get the level were the actual particle is in
    183       do il=2,nz
    184         if (height(il).gt.ztra1(jpart)) then
    185           !hz=il-1
    186           kz=il-1
    187           goto 26
    188         endif
    189       end do
    190 26     continue
    191 
    192   n=memind(2)
    193   if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
    194        n=memind(1)
     161    do il=2,nz
     162      if (height(il).gt.ztra1(jpart)) then
     163        hz=il-1
     164        goto 26
     165      endif
     166    end do
     16726  continue
     168
     169    n=memind(2)
     170    if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
     171    n=memind(1)
    195172
    196173  ! if there is no precipitation or the particle is above the clouds no
    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
    270       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 
    277       if (weta(ks).gt.0.) then
    278         if (clouds_v.ge.4) then
    279   !          BELOW CLOUD SCAVENGING
    280   !          for aerosols and not highliy soluble substances weta=5E-6
    281           wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
    282   !        write(*,*) 'bel. wetscav: ',wetscav
    283 
    284         else ! below_cloud clouds_v is lt 4 and gt 1 -> in cloud scavenging
    285   !        IN CLOUD SCAVENGING
    286   ! BUGFIX tt for nested fields should be ttn
    287   ! sec may 2008
    288           if (ngrid.gt.0) then
    289              act_temp=ttn(ix,jy,hz,n,ngrid)
    290           else
    291              act_temp=tt(ix,jy,hz,n)
    292           endif
     230    do ks=1,nspec      ! loop over species
     231      wetdeposit(ks)=0.
     232      wetscav=0.   
     233
     234  ! NIK 09.12.13: allowed to turn off either below cloud or in-cloud by including TWO separate blocks of if tests
     235
     236
     237  !****i*******************
     238  ! BELOW CLOUD SCAVENGING
     239  !*********************** 
     240      if (weta(ks).gt.0. .and. clouds_v.ge.4) then !if positive below-cloud coefficient given in SPECIES file, and if in below cloud regime
     241 !       for aerosols and not highliy soluble substances weta=5E-6
     242        wetscav=weta(ks)*prec**wetb(ks)  ! scavenging coefficient
     243      endif !end below-cloud
     244
     245
     246  !********************
     247  ! IN CLOUD SCAVENGING
     248  !********************
     249      if (weta_in(ks).gt.0. .and. clouds_v.lt.4) then !if positive in-cloud coefficient given in SPECIES file, and if in in-cloud regime
     250
     251        if (ngrid.gt.0) then
     252          act_temp=ttn(ix,jy,hz,n,ngrid)
     253        else
     254          act_temp=tt(ix,jy,hz,n)
     255        endif
    293256
    294257! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening
    … …  
    297260! wetc_in=0.9 (default)
    298261! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0-no scaling)
    299           cl=weta_in(ks)*prec**wetb_in(ks)
    300           if (dquer(ks).gt.0) then ! is particle
    301             S_i=wetc_in(ks)/cl
    302            else ! is gas
    303             cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
    304             S_i=1/cle
    305            endif
    306            wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
    307   !         write(*,*) 'in. wetscav:'
    308   !    +          ,wetscav,cle,cl,act_temp,prec,clouds_h
     262        cl=weta_in(ks)*prec**wetb_in(ks)
     263        if (dquer(ks).gt.0) then ! is particle
     264          S_i=wetc_in(ks)/cl
     265        else ! is gas
     266          cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
     267          S_i=1/cle
    309268        endif
    310 
    311 
    312   !      if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
    313   !    +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
     269        wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
     270  !     write(*,*) 'in. wetscav:'
     271  !     +  ,wetscav,cle,cl,act_temp,prec,clouds_h
     272       
     273      endif ! end in-cloud
     274
     275
     276     
     277  !**************************************************
     278  ! CALCULATE DEPOSITION
     279  !**************************************************
     280  !     if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
     281  !     +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
     282
     283      if (wetscav.gt.0.) then
    314284        wetdeposit(ks)=xmass1(jpart,ks)* &
    315              (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! wet deposition
    316   !      new particle mass:
    317   !      if (wetdeposit(ks).gt.0) then
    318   !         write(*,*) 'wetdepo: ',wetdeposit(ks),ks
    319   !      endif
    320         restmass = xmass1(jpart,ks)-wetdeposit(ks)
    321          if (ioutputforeachrelease.eq.1) then
    322            kp=npoint(jpart)
    323          else
    324            kp=1
    325          endif
    326         if (restmass .gt. smallnum) then
    327           xmass1(jpart,ks)=restmass
    328   !ccccccccccccccc depostatistic
    329   !       wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
    330   !ccccccccccccccc depostatistic
    331         else
    332           xmass1(jpart,ks)=0.
    333         endif
    334   ! Correct deposited mass to the last time step when radioactive decay of
    335   ! gridded deposited mass was calculated
    336        if (decay(ks).gt.0.) then
    337           wetdeposit(ks)=wetdeposit(ks) &
    338                *exp(abs(ldeltat)*decay(ks))
    339        endif
    340       else  ! weta(k)
    341          wetdeposit(ks)=0.
    342       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 
    411     end do
    412 
    413   ! Sabine Eckhard, June 2008 create deposition runs only for forward runs
     285        (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! wet deposition
     286      else ! if no scavenging
     287        wetdeposit(ks)=0.
     288      endif
     289
     290
     291      restmass = xmass1(jpart,ks)-wetdeposit(ks)
     292      if (ioutputforeachrelease.eq.1) then
     293        kp=npoint(jpart)
     294      else
     295        kp=1
     296      endif
     297      if (restmass .gt. smallnum) then
     298        xmass1(jpart,ks)=restmass
     299  !   depostatistic
     300  !   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
     301  !   depostatistic
     302      else
     303        xmass1(jpart,ks)=0.
     304      endif
     305  !   Correct deposited mass to the last time step when radioactive decay of
     306  !   gridded deposited mass was calculated
     307      if (decay(ks).gt.0.) then
     308        wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
     309      endif
     310
     311
     312
     313    end do !all species
     314
     315  ! Sabine Eckhardt, June 2008 create deposition runs only for forward runs
    414316  ! Add the wet deposition to accumulated amount on output grid and nested output grid
    415317  !*****************************************************************************
    416318
    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 
     319    if (ldirect.eq.1) then
     320      call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
     321      real(ytra1(jpart)),nage,kp)
     322      if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
     323        wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),nage,kp)
     324    endif
    434325
    43532620  continue
    436   end do
     327  end do ! all particles
    437328
    438329end subroutine wetdepo
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG