Changeset 37 for branches/petra/src/wetdepo.f90

Timestamp:

Apr 22, 2015, 3:42:35 PM (9 years ago)

Author:

pesei

Message:

Wet dep quick fix and other small changes. Wet depo quick fix not final yet.

File:

• branches/petra/src/wetdepo.f90 (modified) (11 diffs)

branches/petra/src/wetdepo.f90

@@ -1 +1 @@
 !**********************************************************************
-! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Copyright 1998-2015                                                 *
 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
@@ -39 +39 @@
   ! Code may not be correct for decay of deposition!                           *
   !                                                                            *
+  ! PS, 2/2015: implement wet depo quick fix from 2011/12 in this version
+  !  it implements interpolated and improved clouds
+  !  Also, certain deficiencies for thin clouds fixed also
+  !  Pass itage to wetdepokernel
+  !                                                                            *
   !*****************************************************************************
   !                                                                            *
@@ -71 +76 @@
 
   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,kz,il,interp_time,n
+  integer :: ks, kp, n1,n2, icbot,ictop, indcloud
   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,wetscavold,precsub,f
   real :: wetdeposit(maxspec),restmass
   real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
@@ -97 +101 @@
   !************************
 
+  particle_loop: &
   do jpart=1,numpart
 
     if (itra1(jpart).eq.-999999999) goto 20
-    if(ldirect.eq.1)then
+    if (ldirect.eq.1) then
       if (itra1(jpart).gt.itime) goto 20
     else
       if (itra1(jpart).lt.itime) goto 20
     endif
+    
   ! Determine age class of the particle
     itage=abs(itra1(jpart)-itramem(jpart))
@@ -123 +129 @@
         goto 23
       endif
-    end do
+    enddo
 23  continue
 
@@ -131 +137 @@
 
     if (ngrid.gt.0) then
-      xtn=(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
-      ytn=(ytra1(jpart)-yln(ngrid))*yresoln(ngrid)
+      xtn=real(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
+      ytn=real(ytra1(jpart)-yln(ngrid))*yresoln(ngrid)
       ix=int(xtn)
       jy=int(ytn)
@@ -148 +154 @@
 
     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,nxmax,nymax,1,nx,ny,memind,&
+      real(xtra1(jpart)),real(ytra1(jpart)), &
+      1,memtime(1),memtime(2),interp_time,lsp,convp,cc,icbot,ictop)
     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
+      nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn, &
+      1,memtime(1),memtime(2),interp_time,lsp,convp,cc,icbot,ictop)
+    endif
+
+! PS 2012/2015: subtract a small value, eg 0.01 mm/h, 
+! to remove spurious precip; replaces previous code
+        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
+        kz=il-1
         goto 26
       endif
-    end do
+    enddo
 26  continue
 
     n=memind(2)
     if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
-    n=memind(1)
+      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'
+  ! PS: part of 2011/2012/2015 fix, replaces previous code
+  
+    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 +247 @@
   !**********************************************************
 
+    species_loop: &
     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
+      if (weta(ks).gt.0.) then 
+        ! positive below-cloud coefficient from SPECIES file
+        
+        if (indcloud .eq. 1) then ! below-cloud scavenging
+
+         ! for aerosols and not highly soluble substances weta=5E-6
+           wetscav=weta(ks)*prec**wetb(ks)  ! scavenging coefficient 
+
+        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
 
 ! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening 
@@ -259 +274 @@
 ! wetb_in=0.36 (default)
 ! wetc_in=0.9 (default)
-! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0-no scaling)
+! 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
+          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
-
-
+        wetscav=S_i*prec/3.6E6/wetd_in(ks)/(ictop-icbot) ! 3.6e6 converts mm/h to m/s
+! PS - prevent extrem values of wetscav:
+        wetscavold = 2.e-5*prec**0.8      
+        wetscav = min( 0.1*wetscav, wetscavold ) ! here we sacle the above wetscav. This is now twice, also in wetd_in.
+
+      else ! PS: no cloud diagnosed, old scheme, 
+
+! PS    using with fixed a,b for simplicity, one may wish to change!!
+        wetscav = 2.e-5*prec**0.8 ! was before: 1.e-4*prec**0.62
+
+      endif ! end regime cases
       
-  !**************************************************
-  ! 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
-        wetdeposit(ks)=xmass1(jpart,ks)* &
-        (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! wet deposition
-      else ! if no scavenging
-        wetdeposit(ks)=0.
-      endif
-
-
+  ! calculate deposition 
+      wetdeposit(ks)=xmass1(jpart,ks)* &
+        (1.-exp(-wetscav*abs(ltsample)))*grfraction 
       restmass = xmass1(jpart,ks)-wetdeposit(ks)
       if (ioutputforeachrelease.eq.1) then
@@ -297 +305 @@
       if (restmass .gt. smallnum) then
         xmass1(jpart,ks)=restmass
-  !   depostatistic
-  !   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
-  !   depostatistic
+    !   depostatistic:
+   !!   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
       else
         xmass1(jpart,ks)=0.
@@ -305 +312 @@
   !   Correct deposited mass to the last time step when radioactive decay of
   !   gridded deposited mass was calculated
-      if (decay(ks).gt.0.) then
+      if (decay(ks).gt.0.) &
         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
-  ! Add the wet deposition to accumulated amount on output grid and nested output grid
-  !*****************************************************************************
+
+    else ! weta(k) <= 0
+
+      wetdeposit(ks)=0.
+    
+    endif
+
+  enddo species_loop 
+
+ ! Sabine Eckhardt, June 2008: write deposition only in forward runs
+ ! add the wet deposition from this step 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)
+      call wetdepokernel(nclass(jpart),wetdeposit, &
+        real(xtra1(jpart)),real(ytra1(jpart)),itage,nage,kp)
       if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
-        wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),nage,kp)
-    endif
-
-20  continue
-  end do ! all particles
+        wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),itage,nage,kp)
+    endif
+
+20  continue ! jump here for particles not to be treated
+  enddo particle_loop
 
 end subroutine wetdepo