Changes in trunk/src/wetdepo.f90 [24:20]

File:

• trunk/src/wetdepo.f90 (modified) (9 diffs)

trunk/src/wetdepo.f90

@@ -21 +21 @@
 
 subroutine wetdepo(itime,ltsample,loutnext)
-  !                  i      i        i
+  !                     i      i        i
   !*****************************************************************************
   !                                                                            *
@@ -39 +39 @@
   ! Code may not be correct for decay of deposition!                           *
   !                                                                            *
+  ! Modification by Sabine Eckhart to introduce a new in-/below-cloud
+  ! scheme, not dated
+  ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
   !*****************************************************************************
   !                                                                            *
@@ -71 +74 @@
 
   integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
-  integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
-  integer :: ks, kp
+  integer :: ngrid,itage,nage,hz,il,interp_time, n , clouds_v !NIK scheme
+  integer :: kz !PS scheme
+  integer :: ks, kp, n1,n2, icbot,ictop, indcloud
+  integer :: scheme_number ! NIK==1, PS ==2
   real :: S_i, act_temp, cl, cle ! in cloud scavenging
   real :: clouds_h ! cloud height for the specific grid point
-  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav
+  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav, precsub,f
   real :: wetdeposit(maxspec),restmass
   real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
@@ -98 +103 @@
 
   do jpart=1,numpart
-
     if (itra1(jpart).eq.-999999999) goto 20
     if(ldirect.eq.1)then
@@ -110 +114 @@
       if (itage.lt.lage(nage)) goto 33
     end do
-33  continue
+33   continue
 
 
@@ -119 +123 @@
     do j=numbnests,1,-1
       if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. &
-      (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
+           (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
         ngrid=j
         goto 23
       endif
     end do
-23  continue
+23   continue
 
 
@@ -147 +151 @@
     interp_time=nint(itime-0.5*ltsample)
 
-    if (ngrid.eq.0) then
-      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
-      1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
-      memtime(1),memtime(2),interp_time,lsp,convp,cc)
-    else
-      call interpol_rain_nests(lsprecn,convprecn,tccn, &
-      nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
-      memtime(1),memtime(2),interp_time,lsp,convp,cc)
-    endif
-
-    if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
+!    PS nest case still needs to be implemented!!    
+!    if (ngrid.eq.0) then
+!      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
+!           1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
+!           memtime(1),memtime(2),interp_time,lsp,convp,cc)
+           call interpol_rain(lsprec,convprec,tcc,     &
+            icloudbot,icloudthck,nxmax,nymax,1,nx,ny,         &
+            memind,sngl(xtra1(jpart)),sngl(ytra1(jpart)),1,memtime(1), &
+            memtime(2),interp_time,lsp,convp,cc,icbot,ictop,icmv)  
+!    else
+!      call interpol_rain_nests(lsprecn,convprecn,tccn, &
+!           nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
+!           memtime(1),memtime(2),interp_time,lsp,convp,cc)
+!    endif
+
+
+ !   if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
+ !PS 2012: subtract a small value, eg 0.01 mm/h, to remove spurious precip
+        prec = lsp+convp
+        precsub = 0.01
+        if (prec .lt. precsub) then
+          goto 20
+        else
+          f = (prec-precsub)/prec
+          lsp = f*lsp
+          convp = f*convp
+        endif
+
+
   ! get the level were the actual particle is in
-    do il=2,nz
-      if (height(il).gt.ztra1(jpart)) then
-        hz=il-1
-        goto 26
-      endif
-    end do
-26  continue
-
-    n=memind(2)
-    if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
-    n=memind(1)
+      do il=2,nz
+        if (height(il).gt.ztra1(jpart)) then
+          !hz=il-1
+          kz=il-1
+          goto 26
+        endif
+      end do
+26     continue
+
+  n=memind(2)
+  if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
+       n=memind(1)
 
   ! if there is no precipitation or the particle is above the clouds no
   ! scavenging is done
 
-    if (ngrid.eq.0) then
-      clouds_v=clouds(ix,jy,hz,n)
-      clouds_h=cloudsh(ix,jy,n)
-    else
-      clouds_v=cloudsn(ix,jy,hz,n,ngrid)
-      clouds_h=cloudsnh(ix,jy,n,ngrid)
-    endif
-  !write(*,*) 'there is
-  !    + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
-    if (clouds_v.le.1) goto 20
-  !write (*,*) 'there is scavenging'
+!old scheme
+!  if (ngrid.eq.0) then
+!     clouds_v=clouds(ix,jy,hz,n)
+!     clouds_h=cloudsh(ix,jy,n)
+!  else
+!     clouds_v=cloudsn(ix,jy,hz,n,ngrid)
+!     clouds_h=cloudsnh(ix,jy,n,ngrid)
+!  endif
+!  !write(*,*) 'there is
+!  !    + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
+!  if (clouds_v.le.1) goto 20
+!  !write (*,*) 'there is scavenging'
+
+  ! PS: part of 2011/2012 fix 
+
+        if (ztra1(jpart) .le. float(ictop)) then
+          if (ztra1(jpart) .gt. float(icbot)) then
+            indcloud = 2 ! in-cloud
+          else
+            indcloud = 1 ! below-cloud
+          endif
+        elseif (ictop .eq. icmv) then
+          indcloud = 0 ! no cloud found, use old scheme
+        else
+          goto 20 ! above cloud
+        endif
+
 
   ! 1) Parameterization of the the area fraction of the grid cell where the
@@ -228 +267 @@
   !**********************************************************
 
-    do ks=1,nspec      ! loop over species
-      wetdeposit(ks)=0. 
-      wetscav=0.   
-
-  ! NIK 09.12.13: allowed to turn off either below cloud or in-cloud by including TWO separate blocks of if tests
-
-
-  !****i*******************
-  ! BELOW CLOUD SCAVENGING
-  !***********************  
-      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
- !       for aerosols and not highliy soluble substances weta=5E-6
-        wetscav=weta(ks)*prec**wetb(ks)  ! scavenging coefficient 
-      endif !end below-cloud
-
-
-  !********************
-  ! IN CLOUD SCAVENGING
-  !********************
-      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
-
-        if (ngrid.gt.0) then
-          act_temp=ttn(ix,jy,hz,n,ngrid)
-        else
-          act_temp=tt(ix,jy,hz,n)
-        endif
+    do ks=1,nspec            ! loop over species
+      wetdeposit(ks)=0.
+
+    
+    !conflicting changes to the same routine: 1=NIK 2 =PS
+    scheme_number=2   
+    if (scheme_number.eq.1) then !NIK 
+
+      if (weta(ks).gt.0.) then
+        if (clouds_v.ge.4) then
+  !          BELOW CLOUD SCAVENGING
+  !          for aerosols and not highliy soluble substances weta=5E-6
+          wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
+  !        write(*,*) 'bel. wetscav: ',wetscav
+
+        else ! below_cloud clouds_v is lt 4 and gt 1 -> in cloud scavenging
+  !        IN CLOUD SCAVENGING
+  ! BUGFIX tt for nested fields should be ttn
+  ! sec may 2008
+          if (ngrid.gt.0) then
+             act_temp=ttn(ix,jy,hz,n,ngrid)
+          else
+             act_temp=tt(ix,jy,hz,n)
+          endif
 
 ! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening 
@@ -260 +297 @@
 ! wetc_in=0.9 (default)
 ! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0-no scaling)
-        cl=weta_in(ks)*prec**wetb_in(ks)
-        if (dquer(ks).gt.0) then ! is particle
-          S_i=wetc_in(ks)/cl
-        else ! is gas
-          cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
-          S_i=1/cle
-        endif
-        wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
-  !     write(*,*) 'in. wetscav:'
-  !     +  ,wetscav,cle,cl,act_temp,prec,clouds_h
-       
-      endif ! end in-cloud
-
-
-      
-  !**************************************************
-  ! CALCULATE DEPOSITION 
-  !**************************************************
-  !     if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
-  !     +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
-
-      if (wetscav.gt.0.) then
+          cl=weta_in(ks)*prec**wetb_in(ks)
+          if (dquer(ks).gt.0) then ! is particle
+            S_i=wetc_in(ks)/cl
+           else ! is gas
+            cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
+            S_i=1/cle
+           endif
+           wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
+  !         write(*,*) 'in. wetscav:'
+  !    +          ,wetscav,cle,cl,act_temp,prec,clouds_h
+        endif
+
+
+  !      if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
+  !    +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
         wetdeposit(ks)=xmass1(jpart,ks)* &
-        (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! wet deposition
-      else ! if no scavenging
-        wetdeposit(ks)=0.
-      endif
-
-
-      restmass = xmass1(jpart,ks)-wetdeposit(ks)
-      if (ioutputforeachrelease.eq.1) then
-        kp=npoint(jpart)
-      else
-        kp=1
-      endif
-      if (restmass .gt. smallnum) then
-        xmass1(jpart,ks)=restmass
-  !   depostatistic
-  !   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
-  !   depostatistic
-      else
-        xmass1(jpart,ks)=0.
-      endif
-  !   Correct deposited mass to the last time step when radioactive decay of
-  !   gridded deposited mass was calculated
-      if (decay(ks).gt.0.) then
-        wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
-      endif
-
-
-
-    end do !all species
-
-  ! Sabine Eckhardt, June 2008 create deposition runs only for forward runs
+             (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! wet deposition
+  !      new particle mass:
+  !      if (wetdeposit(ks).gt.0) then
+  !         write(*,*) 'wetdepo: ',wetdeposit(ks),ks
+  !      endif
+        restmass = xmass1(jpart,ks)-wetdeposit(ks)
+         if (ioutputforeachrelease.eq.1) then
+           kp=npoint(jpart)
+         else
+           kp=1
+         endif
+        if (restmass .gt. smallnum) then
+          xmass1(jpart,ks)=restmass
+  !ccccccccccccccc depostatistic
+  !       wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
+  !ccccccccccccccc depostatistic
+        else
+          xmass1(jpart,ks)=0.
+        endif
+  ! Correct deposited mass to the last time step when radioactive decay of
+  ! gridded deposited mass was calculated
+       if (decay(ks).gt.0.) then
+          wetdeposit(ks)=wetdeposit(ks) &
+               *exp(abs(ldeltat)*decay(ks))
+       endif
+      else  ! weta(k)
+         wetdeposit(ks)=0.
+      endif ! weta(k)
+
+     elseif (scheme_number.eq.2) then ! PS
+
+!PS          indcloud=0 ! Use this for FOR TESTING, 
+!PS          will skip the new in/below cloud method !!!             
+
+          if (weta(ks).gt.0.) then
+            if (indcloud .eq. 1) then ! BELOW CLOUD SCAVENGING
+!C               for aerosols and not highliy soluble substances weta=5E-6
+              wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
+!c             write(*,*) 'bel. wetscav: ',wetscav
+            elseif (indcloud .eq. 2) then !  IN CLOUD SCAVENGING
+              if (ngrid.gt.0) then
+                 act_temp=ttn(ix,jy,kz,n,ngrid)
+              else
+                 act_temp=tt(ix,jy,kz,n)
+              endif
+
+! from NIK              
+! weta_in=2.0E-07 (default)
+! wetb_in=0.36 (default)
+! wetc_in=0.9 (default)
+
+
+              cl=2E-7*prec**0.36
+              if (dquer(ks).gt.0) then ! is particle
+                S_i=0.9/cl
+              else ! is gas
+                cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
+                S_i=1/cle
+              endif
+              wetscav=S_i*prec/3.6E6/(ictop-icbot) ! 3.6e6 converts mm/h to m/s
+            else ! PS: no cloud diagnosed, old scheme, 
+!CPS          using with fixed a,b for simplicity, one may wish to change!!
+              wetscav = 1.e-4*prec**0.62
+            endif
+
+
+            wetdeposit(ks)=xmass1(jpart,ks)*    &
+            ! (1.-exp(-wetscav*abs(ltsample)))*fraction  ! wet deposition
+             (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! fraction = grfraction (PS) 
+            restmass = xmass1(jpart,ks)-wetdeposit(ks)
+            if (ioutputforeachrelease.eq.1) then
+              kp=npoint(jpart)
+            else
+              kp=1
+            endif
+            if (restmass .gt. smallnum) then
+              xmass1(jpart,ks)=restmass
+!cccccccccccccccc depostatistic
+!c            wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
+!cccccccccccccccc depostatistic
+            else
+              xmass1(jpart,ks)=0.
+            endif
+!C Correct deposited mass to the last time step when radioactive decay of 
+!C gridded deposited mass was calculated
+            if (decay(ks).gt.0.) then
+              wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
+            endif
+          else  ! weta(k)<0
+             wetdeposit(ks)=0.
+          endif 
+
+     endif !on scheme
+
+
+
+    end do
+
+  ! Sabine Eckhard, June 2008 create deposition runs only for forward runs
   ! Add the wet deposition to accumulated amount on output grid and nested output grid
   !*****************************************************************************
 
-    if (ldirect.eq.1) then
-      call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
-      real(ytra1(jpart)),nage,kp)
-      if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
-        wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),nage,kp)
-    endif
+!    if (ldirect.eq.1) then
+!    call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
+!         real(ytra1(jpart)),nage,kp)
+!    if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
+!         wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), &
+!         nage,kp)
+!    endif
+
+   !PS
+        if (ldirect.eq.1) then
+          call wetdepokernel(nclass(jpart),wetdeposit,     &
+             sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp)
+          if (nested_output.eq.1)                          &
+           call wetdepokernel_nest(nclass(jpart),wetdeposit, &
+             sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp)  
+        endif
+
 
 20  continue
-  end do ! all particles
+  end do
 
 end subroutine wetdepo