Changeset 20 for trunk/src/wetdepo.f90


Ignore:
Timestamp:
Dec 23, 2013, 6:23:38 PM (11 years ago)
Author:
igpis
Message:

move version 9.1.8 form branches to trunk. Contributions from HSO, saeck, pesei, NIK, RT, XKF, IP and others

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/wetdepo.f90

    r4 r20  
    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
    4144  !*****************************************************************************
    4245  !                                                                            *
     
    7174
    7275  integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
    73   integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
    74   integer :: ks, kp
     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
    7580  real :: S_i, act_temp, cl, cle ! in cloud scavenging
    7681  real :: clouds_h ! cloud height for the specific grid point
    77   real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav
     82  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav, precsub,f
    7883  real :: wetdeposit(maxspec),restmass
    7984  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
     
    146151    interp_time=nint(itime-0.5*ltsample)
    147152
    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
     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
    159181
    160182  ! get the level were the actual particle is in
    161183      do il=2,nz
    162184        if (height(il).gt.ztra1(jpart)) then
    163           hz=il-1
     185          !hz=il-1
     186          kz=il-1
    164187          goto 26
    165188        endif
     
    174197  ! scavenging is done
    175198
    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'
     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
    187226
    188227  ! 1) Parameterization of the the area fraction of the grid cell where the
     
    228267  !**********************************************************
    229268
    230     do ks=1,nspec                                  ! loop over species
     269    do ks=1,nspec            ! loop over species
    231270      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
    232277      if (weta(ks).gt.0.) then
    233278        if (clouds_v.ge.4) then
     
    236281          wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
    237282  !        write(*,*) 'bel. wetscav: ',wetscav
     283
    238284        else ! below_cloud clouds_v is lt 4 and gt 1 -> in cloud scavenging
    239285  !        IN CLOUD SCAVENGING
     
    245291             act_temp=tt(ix,jy,hz,n)
    246292          endif
    247           cl=2E-7*prec**0.36
     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)
    248300          if (dquer(ks).gt.0) then ! is particle
    249             S_i=0.9/cl
     301            S_i=wetc_in(ks)/cl
    250302           else ! is gas
    251303            cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
    252304            S_i=1/cle
    253305           endif
    254            wetscav=S_i*prec/3.6E6/clouds_h
     306           wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
    255307  !         write(*,*) 'in. wetscav:'
    256308  !    +          ,wetscav,cle,cl,act_temp,prec,clouds_h
     
    289341         wetdeposit(ks)=0.
    290342      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
    291411    end do
    292412
     
    295415  !*****************************************************************************
    296416
    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
     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
    304434
    30543520  continue
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG