diff --git a/src/FLEXPART.f90 b/src/FLEXPART.f90 index 40c5722ded2114884fdaf54090301f6bb6842289..c2d9aabb1bd27bdc16f558d00429e0751d3e9f22 100644 --- a/src/FLEXPART.f90 +++ b/src/FLEXPART.f90 @@ -65,7 +65,7 @@ program flexpart call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1)) ! FLEXPART version string - flexversion='Version 9.2_beta.r34 (2015-02-12)' + flexversion='Version 9.2_beta.r3? (2015-02-25)' verbosity=0 ! Read the pathnames where input/output files are stored diff --git a/src/com_mod.f90 b/src/com_mod.f90 index cd61d7a102f46904cd29d2f818b5b06b2df5c537..bc6ad539911424d6d941dfcb6a29fd55aa45913e 100755 --- a/src/com_mod.f90 +++ b/src/com_mod.f90 @@ -31,6 +31,7 @@ module com_mod ! * ! Modifications: 15 August 2013 IP, ! 2/2015 PS, add incremental deposition switch +! 2/2015 PS, change variables for cloud diagnostics ! * !******************************************************************************* @@ -364,10 +365,6 @@ module com_mod real :: tth(0:nxmax-1,0:nymax-1,nuvzmax,2) real :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,2) real :: pplev(0:nxmax-1,0:nymax-1,nuvzmax,2) - integer(kind=1) :: clouds(0:nxmax-1,0:nymax-1,nzmax,2) - integer :: cloudsh(0:nxmax-1,0:nymax-1,2) - - ! uu,vv,ww [m/2] wind components in x,y and z direction ! uupol,vvpol [m/s] wind components in polar stereographic projection ! tt [K] temperature data @@ -376,10 +373,6 @@ module com_mod ! rho [kg/m3] air density ! drhodz [kg/m2] vertical air density gradient ! tth,qvh tth,qvh on original eta levels - ! clouds: no cloud, no precipitation 0 - ! cloud, no precipitation 1 - ! rainout conv/lsp dominated 2/3 - ! washout conv/lsp dominated 4/5 ! pplev for the GFS version @@ -405,6 +398,7 @@ module com_mod real :: tropopause(0:nxmax-1,0:nymax-1,1,2) real :: oli(0:nxmax-1,0:nymax-1,1,2) real :: diffk(0:nxmax-1,0:nymax-1,1,2) + integer, dimension (0:nxmax-1,0:nymax-1,2) :: icloudbot, icloudthck ! ps surface pressure ! sd snow depth @@ -425,6 +419,8 @@ module com_mod ! tropopause [m] altitude of thermal tropopause ! oli [m] inverse Obukhov length (1/L) ! diffk [m2/s] diffusion coefficient at reference height + ! icloudbot (m) cloud bottom height + ! icloudthck (m) cloud thickness real :: vdep(0:nxmax-1,0:nymax-1,maxspec,2) @@ -483,8 +479,6 @@ module com_mod real :: ttn(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests) real :: qvn(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests) real :: pvn(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests) - integer(kind=1) :: cloudsn(0:nxmaxn-1,0:nymaxn-1,0:nzmax,2,maxnests) - integer :: cloudsnh(0:nxmaxn-1,0:nymaxn-1,2,maxnests) real :: rhon(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests) real :: drhodzn(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests) real :: tthn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,2,maxnests) @@ -513,6 +507,7 @@ module com_mod real :: olin(0:nxmaxn-1,0:nymaxn-1,1,2,maxnests) real :: diffkn(0:nxmaxn-1,0:nymaxn-1,1,2,maxnests) real :: vdepn(0:nxmaxn-1,0:nymaxn-1,maxspec,2,maxnests) + integer, dimension (0:nxmax-1,0:nymax-1,2,maxnests) :: icloudbotn, icloudthckn !************************************************* diff --git a/src/getfields.f90 b/src/getfields.f90 index fad4e969ce491d88f12f332df4370df5a3ce4ace..38a4ed852a7d74f5f838ade6d9705c2aaa7cb0d7 100644 --- a/src/getfields.f90 +++ b/src/getfields.f90 @@ -40,6 +40,9 @@ subroutine getfields(itime,nstop) ! Changes, Bernd C. Krueger, Feb. 2001: ! Variables tth,qvh,tthn,qvhn (on eta coordinates) in common block. ! Function of nstop extended. + ! + ! 3/2015, PS: put in verbosity output + ! !***************************************************************************** ! * ! Variables: * @@ -82,6 +85,7 @@ subroutine getfields(itime,nstop) integer :: indmin = 1 + call verboutput(verbosity,1,' getfields') ! Check, if wind fields are available for the current time step !************************************************************** @@ -124,11 +128,17 @@ subroutine getfields(itime,nstop) do indj=indmin,numbwf-1 if (ldirect*wftime(indj+1).gt.ldirect*itime) then + call verboutput(verbosity,1,' call readwind') call readwind(indj+1,memind(2),uuh,vvh,wwh) + call verboutput(verbosity,1,' call readwind_nests') call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn) + call verboutput(verbosity,1,' call calcpar') call calcpar(memind(2),uuh,vvh,pvh) + call verboutput(verbosity,1,' call calcpar_nests') call calcpar_nests(memind(2),uuhn,vvhn,pvhn) + call verboutput(verbosity,1,' call verttransform') call verttransform(memind(2),uuh,vvh,wwh,pvh) + call verboutput(verbosity,1,' call verttransform_nests') call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn) memtime(2)=wftime(indj+1) nstop = 1 @@ -147,19 +157,31 @@ subroutine getfields(itime,nstop) if ((ldirect*wftime(indj).le.ldirect*itime).and. & (ldirect*wftime(indj+1).gt.ldirect*itime)) then memind(1)=1 + call verboutput(verbosity,1,' call readwind') call readwind(indj,memind(1),uuh,vvh,wwh) + call verboutput(verbosity,1,' call readwind_nests') call readwind_nests(indj,memind(1),uuhn,vvhn,wwhn) + call verboutput(verbosity,1,' call calcpar') call calcpar(memind(1),uuh,vvh,pvh) + call verboutput(verbosity,1,' call calcpar_nests') call calcpar_nests(memind(1),uuhn,vvhn,pvhn) + call verboutput(verbosity,1,' call verttransform') call verttransform(memind(1),uuh,vvh,wwh,pvh) + call verboutput(verbosity,1,' call verttransform_nests') call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn) memtime(1)=wftime(indj) memind(2)=2 + call verboutput(verbosity,1,' call readwind') call readwind(indj+1,memind(2),uuh,vvh,wwh) + call verboutput(verbosity,1,' call readwind_nests') call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn) + call verboutput(verbosity,1,' call calcpar') call calcpar(memind(2),uuh,vvh,pvh) + call verboutput(verbosity,1,' call calcpar_nests') call calcpar_nests(memind(2),uuhn,vvhn,pvhn) + call verboutput(verbosity,1,' call verttransform') call verttransform(memind(2),uuh,vvh,wwh,pvh) + call verboutput(verbosity,1,' call verttransform_nests') call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn) memtime(2)=wftime(indj+1) nstop = 1 @@ -174,4 +196,6 @@ subroutine getfields(itime,nstop) if (lwindinterv.gt.idiffmax) nstop=3 + call verboutput(verbosity,1,' leaving getfields') + end subroutine getfields diff --git a/src/gridcheck.f90 b/src/gridcheck.f90 index 32b11fe026546a4db6129d7db84bb521773c5f16..72c9e4fc5a055333dd579a8248f6020048850f17 100755 --- a/src/gridcheck.f90 +++ b/src/gridcheck.f90 @@ -31,12 +31,17 @@ subroutine gridcheck ! DATE: 1997-08-06 * ! LAST UPDATE: 1997-10-10 * ! * + !********************************************************************** + ! * ! Update: 1999-02-08, global fields allowed, A. Stohl* ! CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with * ! ECMWF grib_api * ! CHANGE: 03/12/2008, Harald Sodemann, update to f90 with * ! ECMWF grib_api * - ! * + ! + ! 3/2015, PS: add GRIB2 coding for deta/dt dp/deta used by ZAMG + ! + ! !********************************************************************** ! * ! DESCRIPTION: * @@ -143,101 +148,105 @@ subroutine gridcheck if (gribVer.eq.1) then ! GRIB Edition 1 - !print*,'GRiB Edition 1' - !read the grib2 identifiers - call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret) - call grib_check(iret,gribFunction,gribErrorMsg) - call grib_get_int(igrib,'level',isec1(8),iret) - call grib_check(iret,gribFunction,gribErrorMsg) + !print*,'GRiB Edition 1' + !read the grib2 identifiers + call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'level',isec1(8),iret) + call grib_check(iret,gribFunction,gribErrorMsg) - !change code for etadot to code for omega - if (isec1(6).eq.77) then - isec1(6)=135 - endif + !change code for etadot to code for omega + if (isec1(6).eq.77) then + isec1(6)=135 + endif - !print*,isec1(6),isec1(8) + !print*,isec1(6),isec1(8) - else + else ! GRIB 2 - !print*,'GRiB Edition 2' - !read the grib2 identifiers - call grib_get_int(igrib,'discipline',discipl,iret) - call grib_check(iret,gribFunction,gribErrorMsg) - call grib_get_int(igrib,'parameterCategory',parCat,iret) - call grib_check(iret,gribFunction,gribErrorMsg) - call grib_get_int(igrib,'parameterNumber',parNum,iret) - call grib_check(iret,gribFunction,gribErrorMsg) - call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) - call grib_check(iret,gribFunction,gribErrorMsg) - call grib_get_int(igrib,'level',valSurf,iret) - call grib_check(iret,gribFunction,gribErrorMsg) - call grib_get_int(igrib,'paramId',parId,iret) - call grib_check(iret,gribFunction,gribErrorMsg) + !print*,'GRiB Edition 2' + ! read the grib2 identifiers - !print*,discipl,parCat,parNum,typSurf,valSurf - - !convert to grib1 identifiers - isec1(6)=-1 - isec1(7)=-1 - isec1(8)=-1 - isec1(8)=valSurf ! level - if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T - isec1(6)=130 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U - isec1(6)=131 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V - isec1(6)=132 ! indicatorOfParameter - elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q - isec1(6)=133 ! indicatorOfParameter - elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP - isec1(6)=134 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot - isec1(6)=135 ! indicatorOfParameter - elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot - isec1(6)=135 ! indicatorOfParameter - elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP - isec1(6)=151 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U - isec1(6)=165 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V - isec1(6)=166 ! indicatorOfParameter - elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T - isec1(6)=167 ! indicatorOfParameter - elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D - isec1(6)=168 ! indicatorOfParameter - elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD - isec1(6)=141 ! indicatorOfParameter - elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC - isec1(6)=164 ! indicatorOfParameter - elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP - isec1(6)=142 ! indicatorOfParameter - elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP - isec1(6)=143 ! indicatorOfParameter - elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF - isec1(6)=146 ! indicatorOfParameter - elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR - isec1(6)=176 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS - isec1(6)=180 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS - isec1(6)=181 ! indicatorOfParameter - elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO - isec1(6)=129 ! indicatorOfParameter - elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO - isec1(6)=160 ! indicatorOfParameter - elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. & - (typSurf.eq.1)) then ! LSM - isec1(6)=172 ! indicatorOfParameter - else - print*,'***ERROR: undefined GRiB2 message found!',discipl, & - parCat,parNum,typSurf - endif - if(parId .ne. isec1(6) .and. parId .ne. 77) then - write(*,*) 'parId',parId, 'isec1(6)',isec1(6) -! stop - endif + call grib_get_int(igrib,'discipline',discipl,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'parameterCategory',parCat,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'parameterNumber',parNum,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'level',valSurf,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'paramId',parId,iret) + call grib_check(iret,gribFunction,gribErrorMsg) - endif + !print*,discipl,parCat,parNum,typSurf,valSurf + + ! convert to grib1 identifiers + + isec1(6)=-1 + isec1(7)=-1 + isec1(8)=-1 + isec1(8)=valSurf ! level + if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T + isec1(6)=130 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U + isec1(6)=131 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V + isec1(6)=132 ! indicatorOfParameter + elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q + isec1(6)=133 ! indicatorOfParameter + elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP + isec1(6)=134 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually detadtdpdeta + isec1(6)=135 ! indicatorOfParameter + elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually detadtdpdeta + isec1(6)=135 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.8)) then ! Omega, actually detadtdpdeta + isec1(6)=135 ! indicatorOfParameter + elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP + isec1(6)=151 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U + isec1(6)=165 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V + isec1(6)=166 ! indicatorOfParameter + elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T + isec1(6)=167 ! indicatorOfParameter + elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D + isec1(6)=168 ! indicatorOfParameter + elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD + isec1(6)=141 ! indicatorOfParameter + elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC + isec1(6)=164 ! indicatorOfParameter + elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP + isec1(6)=142 ! indicatorOfParameter + elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP + isec1(6)=143 ! indicatorOfParameter + elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF + isec1(6)=146 ! indicatorOfParameter + elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR + isec1(6)=176 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS + isec1(6)=180 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS + isec1(6)=181 ! indicatorOfParameter + elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO + isec1(6)=129 ! indicatorOfParameter + elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO + isec1(6)=160 ! indicatorOfParameter + elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. & + (typSurf.eq.1)) then ! LSM + isec1(6)=172 ! indicatorOfParameter + else + print*,'***ERROR: undefined GRiB2 message found!',discipl, & + parCat,parNum,typSurf + endif + if(parId .ne. isec1(6) .and. parId .ne. 77) then + write(*,*) 'parId',parId, 'isec1(6)',isec1(6) + ! stop + endif + + endif ! GRIB 1 / GRIB 2 cases !get the size and data of the values array if (isec1(6).ne.-1) then diff --git a/src/interpol_rain.f90 b/src/interpol_rain.f90 index e9d71c5a7354c685ef3b120b2003afddfa588fbc..c6b138568d0d659afa58826240f0e6b92d67f082 100644 --- a/src/interpol_rain.f90 +++ b/src/interpol_rain.f90 @@ -1,5 +1,5 @@ !********************************************************************** -! 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 * ! * @@ -19,8 +19,9 @@ ! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * !********************************************************************** -subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, & - ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3) +subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx,ny, & + memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3, & + icbot,ictop) ! i i i i i i i !i i i i i i i i o o o !**************************************************************************** @@ -39,6 +40,10 @@ subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, & ! Author: A. Stohl * ! * ! 30 August 1996 * + !**************************************************************************** + ! Petra Seibert, 2011/2012-2015: + ! Add interpolation of cloud bottom and cloud thickness + ! for fix to SE's new wet scavenging scheme ! * !**************************************************************************** ! * @@ -65,19 +70,22 @@ subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, & ! zt current z coordinate * ! * !**************************************************************************** + + use com_mod, only: icloudbot, icloudthck + use par_mod, only: icmv implicit none integer :: nx,ny,nxmax,nymax,nzmax,memind(2),m,ix,jy,ixp,jyp - integer :: itime,itime1,itime2,level,indexh + integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4 + integer :: icbot,icthck,ictop,ipsum real :: yy1(0:nxmax-1,0:nymax-1,nzmax,2) real :: yy2(0:nxmax-1,0:nymax-1,nzmax,2) real :: yy3(0:nxmax-1,0:nymax-1,nzmax,2) - real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2) + real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),cbot(2),cthck(2) real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4 - ! If point at border of grid -> small displacement into grid !*********************************************************** @@ -111,7 +119,9 @@ subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, & ! Loop over 2 time steps !*********************** + memind_loop: & do m=1,2 + indexh=memind(m) y1(m)=p1*yy1(ix ,jy ,level,indexh) & @@ -126,7 +136,47 @@ subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, & + p2*yy3(ixp,jy ,level,indexh) & + p3*yy3(ix ,jyp,level,indexh) & + p4*yy3(ixp,jyp,level,indexh) - end do + + +! PS clouds: + ip1=1 + ip2=1 + ip3=1 + ip4=1 + if (icloudbot(ix ,jy ,indexh) .eq. icmv) ip1=0 + if (icloudbot(ixp,jy ,indexh) .eq. icmv) ip2=0 + if (icloudbot(ix ,jyp,indexh) .eq. icmv) ip3=0 + if (icloudbot(ixp,jyp,indexh) .eq. icmv) ip4=0 + ipsum = ip1+ip2+ip3+ip4 + if (ipsum .eq. 0) then + cbot(m)=icmv + else + cbot(m)=(ip1*p1*icloudbot(ix ,jy ,indexh) & + + ip2*p2*icloudbot(ixp,jy ,indexh) & + + ip3*p3*icloudbot(ix ,jyp,indexh) & + + ip4*p4*icloudbot(ixp,jyp,indexh)) / ipsum + endif + + ip1=1 + ip2=1 + ip3=1 + ip4=1 + if (icloudthck(ix ,jy ,indexh) .eq. icmv) ip1=0 + if (icloudthck(ixp,jy ,indexh) .eq. icmv) ip2=0 + if (icloudthck(ix ,jyp,indexh) .eq. icmv) ip3=0 + if (icloudthck(ixp,jyp,indexh) .eq. icmv) ip4=0 + ipsum = ip1+ip2+ip3+ip4 + if (ipsum .eq. 0) then + cthck(m)=icmv + else + cthck(m)=(ip1*p1*icloudthck(ix ,jy ,indexh) & + + ip2*p2*icloudthck(ixp,jy ,indexh) & + + ip3*p3*icloudthck(ix ,jyp,indexh) & + + ip4*p4*icloudthck(ixp,jyp,indexh)) / ipsum + endif +! PS end clouds + + enddo memind_loop !************************************ @@ -141,5 +191,22 @@ subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, & yint2=(y2(1)*dt2+y2(2)*dt1)/dt yint3=(y3(1)*dt2+y3(2)*dt1)/dt +! PS clouds: + icbot = nint( (cbot(1)*dt2 + cbot(2)*dt1)/dt ) + if (nint(cbot(1)) .eq. icmv) icbot=cbot(2) + if (nint(cbot(2)) .eq. icmv) icbot=cbot(1) + + icthck = nint( (cthck(1)*dt2 + cthck(2)*dt1)/dt ) + if (nint(cthck(1)) .eq. icmv) icthck=cthck(2) + if (nint(cthck(2)) .eq. icmv) icthck=cthck(1) + + if (icbot .ne. icmv .and. icthck .ne. icmv) then + ictop = icbot + icthck ! convert cloud thickness to cloud top + else + icbot=icmv + ictop=icmv + endif +! PS end clouds + end subroutine interpol_rain diff --git a/src/interpol_rain_nests.f90 b/src/interpol_rain_nests.f90 index 53e64b880b71f89b713a960f73dfd459b3958ad3..7d2e94377131ec5c26e327c167022a8e2256270a 100644 --- a/src/interpol_rain_nests.f90 +++ b/src/interpol_rain_nests.f90 @@ -1,5 +1,5 @@ !********************************************************************** -! 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 * ! * @@ -20,8 +20,9 @@ !********************************************************************** subroutine interpol_rain_nests(yy1,yy2,yy3,nxmaxn,nymaxn,nzmax, & - maxnests,ngrid,nxn,nyn,memind,xt,yt,level,itime1,itime2,itime, & - yint1,yint2,yint3) + maxnests,ngrid,nxn,nyn, & + memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3, & + icbot,ictop) ! i i i i i i ! i i i i i i i i i i i ! o o o @@ -45,6 +46,11 @@ subroutine interpol_rain_nests(yy1,yy2,yy3,nxmaxn,nymaxn,nzmax, & ! 15 March 2000 * ! * !**************************************************************************** + ! Petra Seibert, 2011/2012-2015: + ! Add interpolation of cloud bottom and cloud thickness + ! for fix to SE's new wet scavenging scheme + ! * + !**************************************************************************** ! * ! Variables: * ! * @@ -70,26 +76,27 @@ subroutine interpol_rain_nests(yy1,yy2,yy3,nxmaxn,nymaxn,nzmax, & ! * !**************************************************************************** + use com_mod, only: icloudbotn, icloudthckn + use par_mod, only: icmv + implicit none integer :: maxnests,ngrid - integer :: nxn(maxnests),nyn(maxnests),nxmaxn,nymaxn,nzmax,memind(2) - integer :: m,ix,jy,ixp,jyp,itime,itime1,itime2,level,indexh + integer :: nxn(maxnests),nyn(maxnests),nxmaxn,nymaxn,nzmax,memind(2),m + integer :: ix,jy,ixp,jyp + integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4 + integer :: icbot,ictop,icthck,ipsum real :: yy1(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests) real :: yy2(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests) real :: yy3(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests) - real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2) + real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),cbot(2),cthck(2) real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4 - - ! If point at border of grid -> small displacement into grid !*********************************************************** - if (xt.ge.(real(nxn(ngrid)-1)-0.0001)) & - xt=real(nxn(ngrid)-1)-0.0001 - if (yt.ge.(real(nyn(ngrid)-1)-0.0001)) & - yt=real(nyn(ngrid)-1)-0.0001 + if (xt.ge.(real(nxn(ngrid)-1)-0.0001)) xt=real(nxn(ngrid)-1)-0.0001 + if (yt.ge.(real(nyn(ngrid)-1)-0.0001)) yt=real(nyn(ngrid)-1)-0.0001 @@ -118,7 +125,9 @@ subroutine interpol_rain_nests(yy1,yy2,yy3,nxmaxn,nymaxn,nzmax, & ! Loop over 2 time steps !*********************** + memind_loop: & do m=1,2 + indexh=memind(m) y1(m)=p1*yy1(ix ,jy ,level,indexh,ngrid) & @@ -133,7 +142,46 @@ subroutine interpol_rain_nests(yy1,yy2,yy3,nxmaxn,nymaxn,nzmax, & + p2*yy3(ixp,jy ,level,indexh,ngrid) & + p3*yy3(ix ,jyp,level,indexh,ngrid) & + p4*yy3(ixp,jyp,level,indexh,ngrid) - end do + +! PS clouds: + ip1=1 + ip2=1 + ip3=1 + ip4=1 + if (icloudbotn(ix ,jy ,indexh,ngrid) .eq. icmv) ip1=0 + if (icloudbotn(ixp,jy ,indexh,ngrid) .eq. icmv) ip2=0 + if (icloudbotn(ix ,jyp,indexh,ngrid) .eq. icmv) ip3=0 + if (icloudbotn(ixp,jyp,indexh,ngrid) .eq. icmv) ip4=0 + ipsum = ip1+ip2+ip3+ip4 + if (ipsum .eq. 0) then + cbot(m)=icmv + else + cbot(m)=(ip1*p1*icloudbotn(ix ,jy ,indexh,ngrid) & + + ip2*p2*icloudbotn(ixp,jy ,indexh,ngrid) & + + ip3*p3*icloudbotn(ix ,jyp,indexh,ngrid) & + + ip4*p4*icloudbotn(ixp,jyp,indexh,ngrid)) / ipsum + endif + + ip1=1 + ip2=1 + ip3=1 + ip4=1 + if (icloudthckn(ix ,jy ,indexh,ngrid) .eq. icmv) ip1=0 + if (icloudthckn(ixp,jy ,indexh,ngrid) .eq. icmv) ip2=0 + if (icloudthckn(ix ,jyp,indexh,ngrid) .eq. icmv) ip3=0 + if (icloudthckn(ixp,jyp,indexh,ngrid) .eq. icmv) ip4=0 + ipsum = ip1+ip2+ip3+ip4 + if (ipsum .eq. 0) then + cthck(m)=icmv + else + cthck(m)=(ip1*p1*icloudthckn(ix ,jy ,indexh,ngrid) & + + ip2*p2*icloudthckn(ixp,jy ,indexh,ngrid) & + + ip3*p3*icloudthckn(ix ,jyp,indexh,ngrid) & + + ip4*p4*icloudthckn(ixp,jyp,indexh,ngrid)) / ipsum + endif +! PS end clouds + + enddo memind_loop !************************************ @@ -149,4 +197,22 @@ subroutine interpol_rain_nests(yy1,yy2,yy3,nxmaxn,nymaxn,nzmax, & yint3=(y3(1)*dt2+y3(2)*dt1)/dt +! PS clouds: + icbot = nint( (cbot(1)*dt2 + cbot(2)*dt1)/dt ) + if (nint(cbot(1)) .eq. icmv) icbot=cbot(2) + if (nint(cbot(2)) .eq. icmv) icbot=cbot(1) + + icthck = nint( (cthck(1)*dt2 + cthck(2)*dt1)/dt ) + if (nint(cthck(1)) .eq. icmv) icthck=cthck(2) + if (nint(cthck(2)) .eq. icmv) icthck=cthck(1) + + if (icbot .ne. icmv .and. icthck .ne. icmv) then + ictop = icbot + icthck ! convert cloud thickness to cloud top + else + icbot=icmv + ictop=icmv + endif +! PS end clouds + + end subroutine interpol_rain_nests diff --git a/src/par_mod.f90 b/src/par_mod.f90 index df5d988205ea9950ae99a1870245f812794ac20a..edc7bcf834fef8b3fea8f2c7eb9d0e77cf22369e 100755 --- a/src/par_mod.f90 +++ b/src/par_mod.f90 @@ -102,11 +102,12 @@ module par_mod !******************** integer,parameter :: idiffnorm=10800, idiffmax=2*idiffnorm, minstep=1 + integer,parameter :: itagekernmin=10800 ! idiffnorm [s] normal time interval between two wind fields ! idiffmax [s] maximum time interval between two wind fields ! minstep [s] minimum time step to be used within FLEXPART - + ! itagekernmin [s] minimum particle age [s] for using output kernel !***************************************************************** ! Parameters for polar stereographic projection close to the poles @@ -264,12 +265,22 @@ module par_mod integer,parameter :: unitboundcond=89 !****************************************************** -! integer code for missing values, used in wet scavenging (PS, 2012) +! for use in wet scavenging (PS, 2012, 2015) !****************************************************** - ! integer icmv - ! parameter(icmv=-9999) - integer,parameter :: icmv=-9999 + integer,parameter :: icmv=-9999 ! integer code for missing values + ! some cloud heights used when trying to construct reasonable cloud: + integer,parameter :: icloudtopconvmin = 6000 + integer,parameter :: icloudtopmin = 3000 + integer,parameter :: icloudbot1 = 500, icloudtop1 = 8000 + integer,parameter :: icloudbot2 = 0, icloudtop2 = 10000 + + real, parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics + real, parameter :: rhmininit = 0.90 ! standard condition for presence of clouds +! PS note that original by Sabine Eckhart was 80% +! PS however, for T<-20 C we consider saturation over ice +! PS so I think 90% should be enough + real, parameter :: rhminx = 0.30 ! extreme condition for presence of clouds ! Parameters for testing !******************************************* diff --git a/src/readcommand.f90 b/src/readcommand.f90 index 8941a79221f0887277be30f6eda6c93eb98d0459..e76d2a8c29ec660a33b81b930e1ea14217a45b8d 100644 --- a/src/readcommand.f90 +++ b/src/readcommand.f90 @@ -234,15 +234,18 @@ subroutine readcommand if (old) call skplin(3,unitcommand) read(unitcommand,*,iostat=icmdstat) linit_cond if (icmdstat .gt. 0) & - print*, 'readcommand: linit_cond not read properly',icmdstat,linit_cond + print*, 'readcommand: linit_cond not present or no read properly,& + & - check values: ',icmdstat,linit_cond if (old) call skplin(3,unitcommand) - read(unitcommand,*,iostat=icmdstat) surf_only + read(unitcommand,*,iostat=icmdstat) icmdstat,surf_only if (icmdstat .gt. 0) & - print*, 'readcommand: linit_cond not read properly',icmdstat,surf_only + print*, 'readcommand: surf_only not present or not read properly& + & - check values: ',icmdstat,surf_only if (old) call skplin(3,unitcommand) - read(unitcommand,*,iostat=icmdstat) ldep_incr + read(unitcommand,*,iostat=icmdstat) icmdstat,ldep_incr if (icmdstat .gt. 0) & - print*, 'readcommand: linit_cond not read properly',icmdstat, ldep_incr + print*, 'readcommand: ldep_incr not present or not read properly& + & - check values: ',icmdstat, ldep_incr close(unitcommand) endif ! input format @@ -415,19 +418,7 @@ subroutine readcommand endif - ! For backward runs one releasefield for all releases makes no sense, - ! For quasilag and domainfill ioutputforechrelease is forbidden - !***************************************************************************** - - if ((ldirect .lt. 0) .and. (ioutputforeachrelease .eq. 0)) then - write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####' - write(*,*) '#### FOR BACKWARD RUNS, IOUTPUTFOREACHRLEASE ####' - write(*,*) '#### MUST BE SET TO ONE! ####' - stop - endif - - - ! For backward runs one releasefield for all releases makes no sense, + ! For domainfill runs one releasefield for all releases makes no sense, ! and is "forbidden" !***************************************************************************** diff --git a/src/readspecies.f90 b/src/readspecies.f90 index 820aea543a891afed714d0748c961488541e50f0..a4cfc8e420ffc2bbf8cfda1c0bdc1f8dce62359a 100644 --- a/src/readspecies.f90 +++ b/src/readspecies.f90 @@ -19,22 +19,24 @@ ! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * !********************************************************************** -subroutine readspecies(id_spec,pos_spec) +subroutine readspecies(id_spec,id_pos) !***************************************************************************** ! * ! This routine reads names and physical constants of chemical species/ * - ! radionuclides given in the parameter pos_spec * + ! radionuclides given in the parameter id_pos * ! * ! Author: A. Stohl * ! * ! 11 July 1996 * ! ! Changes: * + ! * ! N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging * ! * ! HSO, 13 August 2013: added optional namelist input ! PS, 2/2015: access= -> position= + ! PS, 4/2015: quick fix for bug associated with end=22 ! * !***************************************************************************** ! * @@ -59,7 +61,7 @@ subroutine readspecies(id_spec,pos_spec) implicit none - integer :: i, pos_spec,j + integer :: i, id_pos,j integer :: idow,ihour,id_spec character(len=3) :: aspecnumb logical :: spec_found @@ -68,7 +70,7 @@ subroutine readspecies(id_spec,pos_spec) real :: pdecay, pweta, pwetb, preldiff, phenry, pf0, pdensity, pdquer real :: pdsigma, pdryvel, pweightmolar, pohreact, pspec_ass, pkao real :: pweta_in, pwetb_in, pwetc_in, pwetd_in - integer :: readerror + integer :: ios ! declare namelist namelist /species_params/ & @@ -99,124 +101,119 @@ subroutine readspecies(id_spec,pos_spec) ! Open the SPECIES file and read species names and properties !************************************************************ - specnum(pos_spec)=id_spec - write(aspecnumb,'(i3.3)') specnum(pos_spec) + specnum(id_pos)=id_spec + write(aspecnumb,'(i3.3)') specnum(id_pos) open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',form='formatted',err=998) - !write(*,*) 'reading SPECIES',specnum(pos_spec) + !write(*,*) 'reading SPECIES',specnum(id_pos) ASSSPEC=.FALSE. ! try namelist input - read(unitspecies,species_params,iostat=readerror) + read(unitspecies,species_params,iostat=ios) close(unitspecies) - if ((pweightmolar.eq.-789.0).or.(readerror.ne.0)) then ! no namelist found - - readerror=1 + if ((pweightmolar.eq.-789.0).or.(ios.ne.0)) then ! no namelist found - open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',err=998) + open(unitspecies,file=path(1) (1:length(1))//'SPECIES/SPECIES_'//aspecnumb, status='old', err=998) do i=1,6 read(unitspecies,*) end do - - read(unitspecies,'(a10)',end=22) species(pos_spec) - ! write(*,*) species(pos_spec) - read(unitspecies,'(f18.1)',end=22) decay(pos_spec) - ! write(*,*) decay(pos_spec) - read(unitspecies,'(e18.1)',end=22) weta(pos_spec) - ! write(*,*) weta(pos_spec) - read(unitspecies,'(f18.2)',end=22) wetb(pos_spec) - ! write(*,*) wetb(pos_spec) + read(unitspecies,'(a10)',end=20) species(id_pos) + ! write(*,*) species(id_pos) + read(unitspecies,'(f18.1)',end=20) decay(id_pos) + ! write(*,*) decay(id_pos) + read(unitspecies,'(e18.1)',end=20) weta(id_pos) + ! write(*,*) weta(id_pos) + read(unitspecies,'(f18.2)',end=20) wetb(id_pos) + ! write(*,*) wetb(id_pos) !*** NIK 31.01.2013: including in-cloud scavening parameters - read(unitspecies,'(e18.1)',end=22) weta_in(pos_spec) - ! write(*,*) weta_in(pos_spec) - read(unitspecies,'(f18.2)',end=22) wetb_in(pos_spec) - ! write(*,*) wetb_in(pos_spec) - read(unitspecies,'(f18.2)',end=22) wetc_in(pos_spec) - ! write(*,*) wetc_in(pos_spec) - read(unitspecies,'(f18.2)',end=22) wetd_in(pos_spec) - ! write(*,*) wetd_in(pos_spec) - - read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec) - ! write(*,*) reldiff(pos_spec) - read(unitspecies,'(e18.1)',end=22) henry(pos_spec) - ! write(*,*) henry(pos_spec) - read(unitspecies,'(f18.1)',end=22) f0(pos_spec) - ! write(*,*) f0(pos_spec) - read(unitspecies,'(e18.1)',end=22) density(pos_spec) - ! write(*,*) density(pos_spec) - read(unitspecies,'(e18.1)',end=22) dquer(pos_spec) - ! write(*,*) dquer(pos_spec) - read(unitspecies,'(e18.1)',end=22) dsigma(pos_spec) - ! write(*,*) dsigma(pos_spec) - read(unitspecies,'(f18.2)',end=22) dryvel(pos_spec) - ! write(*,*) dryvel(pos_spec) - read(unitspecies,'(f18.2)',end=22) weightmolar(pos_spec) - ! write(*,*) weightmolar(pos_spec) - read(unitspecies,'(e18.1)',end=22) ohreact(pos_spec) - ! write(*,*) ohreact(pos_spec) - read(unitspecies,'(i18)',end=22) spec_ass(pos_spec) - ! write(*,*) spec_ass(pos_spec) - read(unitspecies,'(f18.2)',end=22) kao(pos_spec) - ! write(*,*) kao(pos_spec) - - pspecies=species(pos_spec) - pdecay=decay(pos_spec) - pweta=weta(pos_spec) - pwetb=wetb(pos_spec) - pweta_in=weta_in(pos_spec) - pwetb_in=wetb_in(pos_spec) - pwetc_in=wetc_in(pos_spec) - pwetd_in=wetd_in(pos_spec) - preldiff=reldiff(pos_spec) - phenry=henry(pos_spec) - pf0=f0(pos_spec) - pdensity=density(pos_spec) - pdquer=dquer(pos_spec) - pdsigma=dsigma(pos_spec) - pdryvel=dryvel(pos_spec) - pweightmolar=weightmolar(pos_spec) - pohreact=ohreact(pos_spec) - pspec_ass=spec_ass(pos_spec) - pkao=kao(pos_spec) + read(unitspecies,'(e18.1)',end=20) weta_in(id_pos) + ! write(*,*) weta_in(id_pos) + read(unitspecies,'(f18.2)',end=20) wetb_in(id_pos) + ! write(*,*) wetb_in(id_pos) + read(unitspecies,'(f18.2)',end=20) wetc_in(id_pos) + ! write(*,*) wetc_in(id_pos) + read(unitspecies,'(f18.2)',end=20) wetd_in(id_pos) + ! write(*,*) wetd_in(id_pos) + + read(unitspecies,'(f18.1)',end=20) reldiff(id_pos) + ! write(*,*) reldiff(id_pos) + read(unitspecies,'(e18.1)',end=20) henry(id_pos) + ! write(*,*) henry(id_pos) + read(unitspecies,'(f18.1)',end=20) f0(id_pos) + ! write(*,*) f0(id_pos) + read(unitspecies,'(e18.1)',end=20) density(id_pos) + ! write(*,*) density(id_pos) + read(unitspecies,'(e18.1)',end=20) dquer(id_pos) + ! write(*,*) dquer(id_pos) + read(unitspecies,'(e18.1)',end=20) dsigma(id_pos) + ! write(*,*) dsigma(id_pos) + read(unitspecies,'(f18.2)',end=20) dryvel(id_pos) + ! write(*,*) dryvel(id_pos) + read(unitspecies,'(f18.2)',end=20) weightmolar(id_pos) + ! write(*,*) weightmolar(id_pos) + read(unitspecies,'(e18.1)',end=20) ohreact(id_pos) + ! write(*,*) ohreact(id_pos) + read(unitspecies,'(i18)',end=20) spec_ass(id_pos) + ! write(*,*) spec_ass(id_pos) + read(unitspecies,'(f18.2)',end=20) kao(id_pos) + ! write(*,*) kao(id_pos) + + pspecies=species(id_pos) + pdecay=decay(id_pos) + pweta=weta(id_pos) + pwetb=wetb(id_pos) + pweta_in=weta_in(id_pos) + pwetb_in=wetb_in(id_pos) + pwetc_in=wetc_in(id_pos) + pwetd_in=wetd_in(id_pos) + preldiff=reldiff(id_pos) + phenry=henry(id_pos) + pf0=f0(id_pos) + pdensity=density(id_pos) + pdquer=dquer(id_pos) + pdsigma=dsigma(id_pos) + pdryvel=dryvel(id_pos) + pweightmolar=weightmolar(id_pos) + pohreact=ohreact(id_pos) + pspec_ass=spec_ass(id_pos) + pkao=kao(id_pos) else - species(pos_spec)=pspecies - decay(pos_spec)=pdecay - weta(pos_spec)=pweta - wetb(pos_spec)=pwetb - weta_in(pos_spec)=pweta_in - wetb_in(pos_spec)=pwetb_in - wetc_in(pos_spec)=pwetc_in - wetd_in(pos_spec)=pwetd_in - reldiff(pos_spec)=preldiff - henry(pos_spec)=phenry - f0(pos_spec)=pf0 - density(pos_spec)=pdensity - dquer(pos_spec)=pdquer - dsigma(pos_spec)=pdsigma - dryvel(pos_spec)=pdryvel - weightmolar(pos_spec)=pweightmolar - ohreact(pos_spec)=pohreact - spec_ass(pos_spec)=pspec_ass - kao(pos_spec)=pkao + species(id_pos)=pspecies + decay(id_pos)=pdecay + weta(id_pos)=pweta + wetb(id_pos)=pwetb + weta_in(id_pos)=pweta_in + wetb_in(id_pos)=pwetb_in + wetc_in(id_pos)=pwetc_in + wetd_in(id_pos)=pwetd_in + reldiff(id_pos)=preldiff + henry(id_pos)=phenry + f0(id_pos)=pf0 + density(id_pos)=pdensity + dquer(id_pos)=pdquer + dsigma(id_pos)=pdsigma + dryvel(id_pos)=pdryvel + weightmolar(id_pos)=pweightmolar + ohreact(id_pos)=pohreact + spec_ass(id_pos)=pspec_ass + kao(id_pos)=pkao endif - i=pos_spec - - if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then - if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set + if ((weta(id_pos).gt.0).and.(henry(id_pos).le.0)) then + if (dquer(id_pos).le.0) goto 996 ! no particle, no henry set endif - if (spec_ass(pos_spec).gt.0) then + if (spec_ass(id_pos).gt.0) then spec_found=.FALSE. - do j=1,pos_spec-1 - if (spec_ass(pos_spec).eq.specnum(j)) then - spec_ass(pos_spec)=j + do j=1,id_pos-1 + if (spec_ass(id_pos).eq.specnum(j)) then + spec_ass(id_pos)=j spec_found=.TRUE. ASSSPEC=.TRUE. endif @@ -226,10 +223,10 @@ subroutine readspecies(id_spec,pos_spec) endif endif - if (dsigma(i).eq.1.) dsigma(i)=1.0001 ! avoid floating exception - if (dsigma(i).eq.0.) dsigma(i)=1.0001 ! avoid floating exception + if (dsigma(id_pos).eq.1.) dsigma(id_pos)=1.0001 ! avoid floating exception + if (dsigma(id_pos).eq.0.) dsigma(id_pos)=1.0001 ! avoid floating exception - if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then + if ((reldiff(id_pos).gt.0.).and.(density(id_pos).gt.0.)) then write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES" ####' write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH ####' write(*,*) '#### PARTICLE AND GAS. ####' @@ -245,23 +242,24 @@ subroutine readspecies(id_spec,pos_spec) ! default values set to 1 do j=1,24 ! initialize everything to no variation - area_hour(i,j)=1. - point_hour(i,j)=1. + area_hour(id_pos,j)=1. + point_hour(id_pos,j)=1. end do do j=1,7 - area_dow(i,j)=1. - point_dow(i,j)=1. + area_dow(id_pos,j)=1. + point_dow(id_pos,j)=1. end do - if (readerror.ne.0) then ! text format input + if (ios.ne.0) then ! text format input - read(unitspecies,*,end=22) + read(unitspecies,*,iostat=ios) + if (ios .ne. 0) goto 22 ! ifort does not like err= do j=1,24 ! 24 hours, starting with 0-1 local time - read(unitspecies,*) ihour,area_hour(i,j),point_hour(i,j) + read(unitspecies,*) ihour,area_hour(id_pos,j),point_hour(id_pos,j) end do read(unitspecies,*) do j=1,7 ! 7 days of the week, starting with Monday - read(unitspecies,*) idow,area_dow(i,j),point_dow(i,j) + read(unitspecies,*) idow,area_dow(id_pos,j),point_dow(id_pos,j) end do endif diff --git a/src/readwind.f90 b/src/readwind.f90 index b3c08a1fe9048e764a3efbb7f90baeee959aba4f..039b9765ce5811c39874b05f4c49cb7b00ee4aa6 100644 --- a/src/readwind.f90 +++ b/src/readwind.f90 @@ -39,6 +39,8 @@ subroutine readwind(indj,n,uuh,vvh,wwh) ! Changes, Bernd C. Krueger, Feb. 2001: ! Variables tth and qvh (on eta coordinates) in common block ! 2/2015, PS: add missing paramter iret in call to grib subr + ! 3/2015, PS: add some verbosity messages + ! 3/2015, PS: add GRIB2 coding for deta/dt dp/deta used by ZAMG ! !********************************************************************** ! * @@ -112,7 +114,11 @@ subroutine readwind(indj,n,uuh,vvh,wwh) ! ! OPENING OF DATA FILE (GRIB CODE) ! -5 call grib_open_file(ifile,path(3)(1:length(3)) & + +5 continue + call verboutput(verbosity,1,' call grib_open_file '//path(3)(1:length(3))& + //trim(wfname(indj))) + call grib_open_file(ifile,path(3)(1:length(3)) & //trim(wfname(indj)),'r',iret) if (iret.ne.GRIB_SUCCESS) then goto 888 ! ERROR DETECTED @@ -126,6 +132,7 @@ subroutine readwind(indj,n,uuh,vvh,wwh) ! ! GET NEXT FIELDS ! + call verboutput(verbosity,2,' call grib_new_from_file') call grib_new_from_file(ifile,igrib,iret) if (iret.eq.GRIB_END_OF_FILE) then goto 50 ! EOF DETECTED @@ -134,109 +141,118 @@ subroutine readwind(indj,n,uuh,vvh,wwh) endif !first see if we read GRIB1 or GRIB2 + call verboutput(verbosity,2,' call grib_get_int') call grib_get_int(igrib,'editionNumber',gribVer,iret) + call verboutput(verbosity,2,' call grib_check') call grib_check(iret,gribFunction,gribErrorMsg) if (gribVer.eq.1) then ! GRIB Edition 1 - !print*,'GRiB Edition 1' - !read the grib2 identifiers - call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret) - call grib_check(iret,gribFunction,gribErrorMsg) - call grib_get_int(igrib,'level',isec1(8),iret) - call grib_check(iret,gribFunction,gribErrorMsg) + !print*,'GRiB Edition 1' - !change code for etadot to code for omega - if (isec1(6).eq.77) then - isec1(6)=135 - endif + ! read the grib2 identifiers - conversion_factor=1. + call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'level',isec1(8),iret) + call grib_check(iret,gribFunction,gribErrorMsg) - else + !change code for etadot to code for omega + if (isec1(6).eq.77) then + isec1(6)=135 + endif - !print*,'GRiB Edition 2' - !read the grib2 identifiers - call grib_get_int(igrib,'discipline',discipl,iret) - call grib_check(iret,gribFunction,gribErrorMsg) - call grib_get_int(igrib,'parameterCategory',parCat,iret) - call grib_check(iret,gribFunction,gribErrorMsg) - call grib_get_int(igrib,'parameterNumber',parNum,iret) - call grib_check(iret,gribFunction,gribErrorMsg) - call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) - call grib_check(iret,gribFunction,gribErrorMsg) - call grib_get_int(igrib,'level',valSurf,iret) - call grib_check(iret,gribFunction,gribErrorMsg) - call grib_get_int(igrib,'paramId',parId,iret) - call grib_check(iret,gribFunction,gribErrorMsg) + conversion_factor=1. - !print*,discipl,parCat,parNum,typSurf,valSurf - - !convert to grib1 identifiers - isec1(6)=-1 - isec1(7)=-1 - isec1(8)=-1 - isec1(8)=valSurf ! level - conversion_factor=1. - if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T - isec1(6)=130 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U - isec1(6)=131 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V - isec1(6)=132 ! indicatorOfParameter - elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q - isec1(6)=133 ! indicatorOfParameter - elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP - isec1(6)=134 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot - isec1(6)=135 ! indicatorOfParameter - elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot - isec1(6)=135 ! indicatorOfParameter - elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP - isec1(6)=151 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U - isec1(6)=165 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V - isec1(6)=166 ! indicatorOfParameter - elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T - isec1(6)=167 ! indicatorOfParameter - elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D - isec1(6)=168 ! indicatorOfParameter - elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD - isec1(6)=141 ! indicatorOfParameter - conversion_factor=1000. - elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC - isec1(6)=164 ! indicatorOfParameter - elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP - isec1(6)=142 ! indicatorOfParameter - elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP - isec1(6)=143 ! indicatorOfParameter - conversion_factor=1000. - elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF - isec1(6)=146 ! indicatorOfParameter - elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR - isec1(6)=176 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS - isec1(6)=180 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS - isec1(6)=181 ! indicatorOfParameter - elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO - isec1(6)=129 ! indicatorOfParameter - elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO - isec1(6)=160 ! indicatorOfParameter - elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. & - (typSurf.eq.1)) then ! LSM - isec1(6)=172 ! indicatorOfParameter - else - print*,'***WARNING: undefined GRiB2 message found!',discipl, & - parCat,parNum,typSurf - endif - if(parId .ne. isec1(6) .and. parId .ne. 77) then - write(*,*) 'parId',parId, 'isec1(6)',isec1(6) -! stop - endif + else ! GRIB 2 - endif + !print*,'GRiB Edition 2' + + ! read the grib2 identifiers + + call grib_get_int(igrib,'discipline',discipl,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'parameterCategory',parCat,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'parameterNumber',parNum,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'level',valSurf,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'paramId',parId,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + + !print*,discipl,parCat,parNum,typSurf,valSurf + + ! convert to grib1 identifiers + + isec1(6)=-1 + isec1(7)=-1 + isec1(8)=-1 + isec1(8)=valSurf ! level + conversion_factor=1. + if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T + isec1(6)=130 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U + isec1(6)=131 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V + isec1(6)=132 ! indicatorOfParameter + elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q + isec1(6)=133 ! indicatorOfParameter + elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP + isec1(6)=134 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually detadtdpdeta + isec1(6)=135 ! indicatorOfParameter + elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually detadtdpdeta + isec1(6)=135 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.8)) then ! Omega, actually detadtdpdeta + isec1(6)=135 ! indicatorOfParameter + elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP + isec1(6)=151 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U + isec1(6)=165 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V + isec1(6)=166 ! indicatorOfParameter + elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T + isec1(6)=167 ! indicatorOfParameter + elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D + isec1(6)=168 ! indicatorOfParameter + elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD + isec1(6)=141 ! indicatorOfParameter + conversion_factor=1000. + elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC + isec1(6)=164 ! indicatorOfParameter + elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP + isec1(6)=142 ! indicatorOfParameter + elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP + isec1(6)=143 ! indicatorOfParameter + conversion_factor=1000. + elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF + isec1(6)=146 ! indicatorOfParameter + elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR + isec1(6)=176 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS + isec1(6)=180 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS + isec1(6)=181 ! indicatorOfParameter + elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO + isec1(6)=129 ! indicatorOfParameter + elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO + isec1(6)=160 ! indicatorOfParameter + elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. & + (typSurf.eq.1)) then ! LSM + isec1(6)=172 ! indicatorOfParameter + else + print*,'***WARNING: undefined GRiB2 message found!',discipl, & + parCat,parNum,typSurf + endif + if(parId .ne. isec1(6) .and. parId .ne. 77) then + write(*,*) 'parId',parId, 'isec1(6)',isec1(6) + ! stop + endif + + endif ! end GRIB 1 / GRIB 2 cases !HSO get the size and data of the values array if (isec1(6).ne.-1) then @@ -350,6 +366,7 @@ subroutine readwind(indj,n,uuh,vvh,wwh) end do end do + call verboutput(verbosity,2,' call grib_release') call grib_release(igrib) goto 10 !! READ NEXT LEVEL OR PARAMETER ! @@ -379,6 +396,7 @@ subroutine readwind(indj,n,uuh,vvh,wwh) !************************************************************************* if (xglobal) then + call verboutput(verbosity,2,' call shift_field_0') call shift_field_0(ewss,nxfield,ny) call shift_field_0(nsss,nxfield,ny) call shift_field_0(oro,nxfield,ny) @@ -428,6 +446,7 @@ subroutine readwind(indj,n,uuh,vvh,wwh) hlev1=fu*(plev1-ps(i,j,1,n)) ! HEIGTH OF FIRST MODEL LAYER ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2) fflev1=sqrt(uuh(i,j,3)**2+vvh(i,j,3)**2) + call verboutput(verbosity,2,' call pbl_profile') call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, & tt2(i,j,1,n),tth(i,j,3,n),ff10m,fflev1, & surfstr(i,j,1,n),sshf(i,j,1,n)) @@ -456,7 +475,9 @@ subroutine readwind(indj,n,uuh,vvh,wwh) if(iumax.ne.nuvz-1) stop 'READWIND: NUVZ NOT CONSISTENT' if(iwmax.ne.nwz) stop 'READWIND: NWZ NOT CONSISTENT' + call verboutput(verbosity,1,' leaving readwind') return + 888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD #### ' write(*,*) ' #### ',wfname(indj),' #### ' write(*,*) ' #### IS NOT GRIB FORMAT !!! #### ' diff --git a/src/readwind_nests.f90 b/src/readwind_nests.f90 index ada5332902fc125bdbfef874a848114865f28685..9865d58887e417ba2442bba123950dca4bbc7b5a 100644 --- a/src/readwind_nests.f90 +++ b/src/readwind_nests.f90 @@ -33,11 +33,15 @@ subroutine readwind_nests(indj,n,uuhn,vvhn,wwhn) ! Last update: 17 October 2000, A. Stohl * ! * !***************************************************************************** - ! Changes, Bernd C. Krueger, Feb. 2001: * + ! CHANGES: + ! + ! 2/2001, Bernc C. Krueger: * ! Variables tthn and qvhn (on eta coordinates) in common block * - ! CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with ECMWF grib_api * - ! CHANGE: 03/12/2008, Harald Sodemann, update to f90 with ECMWF grib_api * - ! 2/2015, PS: add missing paramter iret in call to grib subr + ! 11/01/2008, Harald Sodemann, GRIB1/2 input with ECMWF grib_api * + ! 03/12/2008, Harald Sodemann, update to f90 with ECMWF grib_api * + ! 2/2015, PS: add missing paramter iret in call to grib subr + ! 3/2015, PS: add verbosity output + ! 3/2015, PS: add GRIB2 coding for deta/dt dp/deta used by ZAMG ! !***************************************************************************** @@ -83,6 +87,8 @@ subroutine readwind_nests(indj,n,uuhn,vvhn,wwhn) !HSO grib api error messages character(len=24) :: gribErrorMsg = 'Error reading grib file' character(len=20) :: gribFunction = 'readwind_nests' + + call verboutput(verbosity,1,' entered readwind_nests') do l=1,numbnests hflswitch=.false. @@ -99,7 +105,10 @@ subroutine readwind_nests(indj,n,uuhn,vvhn,wwhn) ! OPENING OF DATA FILE (GRIB CODE) ! -5 call grib_open_file(ifile,path(numpath+2*(l-1)+1) & +5 continue + call verboutput(verbosity,1,' call grib_open_file '//& + path(numpath+2*(l-1)+1)(1:length(numpath+2*(l-1)+1))//trim(wfnamen(l,indj))) + call grib_open_file(ifile,path(numpath+2*(l-1)+1) & (1:length(numpath+2*(l-1)+1))//trim(wfnamen(l,indj)),'r') if (iret.ne.GRIB_SUCCESS) then goto 888 ! ERROR DETECTED @@ -126,103 +135,108 @@ subroutine readwind_nests(indj,n,uuhn,vvhn,wwhn) if (gribVer.eq.1) then ! GRIB Edition 1 - !print*,'GRiB Edition 1' - !read the grib2 identifiers - call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret) - call grib_check(iret,gribFunction,gribErrorMsg) - call grib_get_int(igrib,'level',isec1(8),iret) - call grib_check(iret,gribFunction,gribErrorMsg) - - !change code for etadot to code for omega - if (isec1(6).eq.77) then - isec1(6)=135 - endif + !print*,'GRiB Edition 1' + + ! read the grib2 identifiers + + call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'level',isec1(8),iret) + call grib_check(iret,gribFunction,gribErrorMsg) - conversion_factor=1. + !change code for etadot to code for omega + if (isec1(6).eq.77) then + isec1(6)=135 + endif + conversion_factor=1. else - !print*,'GRiB Edition 2' - !read the grib2 identifiers - call grib_get_int(igrib,'discipline',discipl,iret) - call grib_check(iret,gribFunction,gribErrorMsg) - call grib_get_int(igrib,'parameterCategory',parCat,iret) - call grib_check(iret,gribFunction,gribErrorMsg) - call grib_get_int(igrib,'parameterNumber',parNum,iret) - call grib_check(iret,gribFunction,gribErrorMsg) - call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) - call grib_check(iret,gribFunction,gribErrorMsg) - call grib_get_int(igrib,'level',valSurf,iret) - call grib_check(iret,gribFunction,gribErrorMsg) - call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consisitent with new readwind.f90 - call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consisitent with new readwind.f90 - - !print*,discipl,parCat,parNum,typSurf,valSurf - - !convert to grib1 identifiers - isec1(6)=-1 - isec1(7)=-1 - isec1(8)=-1 - isec1(8)=valSurf ! level - conversion_factor=1. - if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T - isec1(6)=130 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U - isec1(6)=131 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V - isec1(6)=132 ! indicatorOfParameter - elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q - isec1(6)=133 ! indicatorOfParameter - elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP - isec1(6)=134 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot ! - isec1(6)=135 ! indicatorOfParameter - elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added by mc to make it consisitent with new readwind.f90 - isec1(6)=135 ! indicatorOfParameter !added by mc to make it consisitent with new readwind.f90 - elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP - isec1(6)=151 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U - isec1(6)=165 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V - isec1(6)=166 ! indicatorOfParameter - elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T - isec1(6)=167 ! indicatorOfParameter - elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D - isec1(6)=168 ! indicatorOfParameter - elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD - isec1(6)=141 ! indicatorOfParameter - conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90 - elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC !added by mc to make it consisitent with new readwind.f90 - isec1(6)=164 ! indicatorOfParameter - elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP !added by mc to make it consisitent with new readwind.f90 - isec1(6)=142 ! indicatorOfParameter - elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP - isec1(6)=143 ! indicatorOfParameter - conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90 - elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF - isec1(6)=146 ! indicatorOfParameter - elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR - isec1(6)=176 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS !added by mc to make it consisitent with new readwind.f90 - isec1(6)=180 ! indicatorOfParameter - elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS !added by mc to make it consisitent with new readwind.f90 - isec1(6)=181 ! indicatorOfParameter - elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO - isec1(6)=129 ! indicatorOfParameter - elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO !added by mc to make it consisitent with new readwind.f90 - isec1(6)=160 ! indicatorOfParameter - elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. & - (typSurf.eq.1)) then ! LSM - isec1(6)=172 ! indicatorOfParameter - else - print*,'***WARNING: undefined GRiB2 message found!',discipl, & - parCat,parNum,typSurf - endif - if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consisitent with new readwind.f90 - write(*,*) 'parId',parId, 'isec1(6)',isec1(6) ! -! stop - endif + !print*,'GRiB Edition 2' + + ! read the grib2 identifiers + + call grib_get_int(igrib,'discipline',discipl,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'parameterCategory',parCat,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'parameterNumber',parNum,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'level',valSurf,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consistent with new readwind.f90 + call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consistent with new readwind.f90 + + !print*,discipl,parCat,parNum,typSurf,valSurf + + !convert to grib1 identifiers + isec1(6)=-1 + isec1(7)=-1 + isec1(8)=-1 + isec1(8)=valSurf ! level + conversion_factor=1. + if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T + isec1(6)=130 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U + isec1(6)=131 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V + isec1(6)=132 ! indicatorOfParameter + elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q + isec1(6)=133 ! indicatorOfParameter + elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP + isec1(6)=134 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually detadtdpdeta + isec1(6)=135 ! indicatorOfParameter + elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually detadtdpdeta + isec1(6)=135 ! indicatorOfParameter !added by mc to make it consistent with new readwind.f90 + elseif ((parCat.eq.2).and.(parNum.eq.8)) then ! Omega, actually detadtdpdeta + isec1(6)=135 ! indicatorOfParameter !added by mc to make it consistent with new readwind.f90 + elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP + isec1(6)=151 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U + isec1(6)=165 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V + isec1(6)=166 ! indicatorOfParameter + elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T + isec1(6)=167 ! indicatorOfParameter + elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D + isec1(6)=168 ! indicatorOfParameter + elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD + isec1(6)=141 ! indicatorOfParameter + conversion_factor=1000. !added by mc to make it consistent with new readwind.f90 + elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC !added by mc to make it consistent with new readwind.f90 + isec1(6)=164 ! indicatorOfParameter + elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP !added by mc to make it consistent with new readwind.f90 + isec1(6)=142 ! indicatorOfParameter + elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP + isec1(6)=143 ! indicatorOfParameter + conversion_factor=1000. !added by mc to make it consistent with new readwind.f90 + elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF + isec1(6)=146 ! indicatorOfParameter + elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR + isec1(6)=176 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS !added by mc to make it consistent with new readwind.f90 + isec1(6)=180 ! indicatorOfParameter + elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS !added by mc to make it consistent with new readwind.f90 + isec1(6)=181 ! indicatorOfParameter + elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO + isec1(6)=129 ! indicatorOfParameter + elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO !added by mc to make it consistent with new readwind.f90 + isec1(6)=160 ! indicatorOfParameter + elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. & + (typSurf.eq.1)) then ! LSM + isec1(6)=172 ! indicatorOfParameter + else + print*,'***WARNING: undefined GRiB2 message found!',discipl, & + parCat,parNum,typSurf + endif + if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consistent with new readwind.f90 + write(*,*) 'parId',parId, 'isec1(6)',isec1(6) ! + ! stop + endif endif @@ -245,15 +259,15 @@ subroutine readwind_nests(indj,n,uuhn,vvhn,wwhn) call grib_check(iret,gribFunction,gribErrorMsg) ! CHECK GRID SPECIFICATIONS if(isec2(2).ne.nxn(l)) stop & - 'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL' + 'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL' if(isec2(3).ne.nyn(l)) stop & - 'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL' - if(isec2(12)/2-1.ne.nlev_ec) stop 'READWIND: VERTICAL DISCRET& - IZATION NOT CONSISTENT FOR A NESTING LEVEL' + 'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL' + if(isec2(12)/2-1.ne.nlev_ec) stop & + 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT FOR A NESTING LEVEL' endif ! ifield !HSO get the second part of the grid dimensions only from GRiB1 messages - if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then ! !added by mc to make it consisitent with new readwind.f90 + if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then ! !added by mc to make it consistent with new readwind.f90 call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', & xauxin,iret) call grib_check(iret,gribFunction,gribErrorMsg) @@ -293,7 +307,7 @@ subroutine readwind_nests(indj,n,uuhn,vvhn,wwhn) if(isec1(6).eq.135) wwhn(i,j,nlev_ec-k+1,l)= &!! W VELOCITY zsec4(nxn(l)*(nyn(l)-j-1)+i+1) if(isec1(6).eq.141) sdn(i,j,1,n,l)= &!! SNOW DEPTH - zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90! + zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consistent with new readwind.f90! if(isec1(6).eq.151) msln(i,j,1,n,l)= &!! SEA LEVEL PRESS. zsec4(nxn(l)*(nyn(l)-j-1)+i+1) if(isec1(6).eq.164) tccn(i,j,1,n,l)= &!! CLOUD COVER @@ -311,7 +325,7 @@ subroutine readwind_nests(indj,n,uuhn,vvhn,wwhn) if (lsprecn(i,j,1,n,l).lt.0.) lsprecn(i,j,1,n,l)=0. endif if(isec1(6).eq.143) then !! CONVECTIVE PREC. - convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90 + convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consistent with new readwind.f90 if (convprecn(i,j,1,n,l).lt.0.) convprecn(i,j,1,n,l)=0. endif if(isec1(6).eq.146) sshfn(i,j,1,n,l)= &!! SENS. HEAT FLUX @@ -420,7 +434,10 @@ subroutine readwind_nests(indj,n,uuhn,vvhn,wwhn) end do + call verboutput(verbosity,1,' leaving readwind_nests') + return + 888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD #### ' write(*,*) ' #### ',wfnamen(l,indj),' FOR NESTING LEVEL #### ' write(*,*) ' #### ',l,' IS NOT GRIB FORMAT !!! #### ' @@ -430,5 +447,6 @@ subroutine readwind_nests(indj,n,uuhn,vvhn,wwhn) 999 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD #### ' write(*,*) ' #### ',wfnamen(l,indj),' #### ' write(*,*) ' #### CANNOT BE OPENED FOR NESTING LEVEL ',l,'####' + call verboutput(verbosity,1,' leaving readwind') end subroutine readwind_nests diff --git a/src/timemanager.f90 b/src/timemanager.f90 index 1d70ce482e4e7f41a3b095db3f0e7481812baac2..04b50191f17b0b52364861eb7db8714987d2d2e5 100644 --- a/src/timemanager.f90 +++ b/src/timemanager.f90 @@ -136,12 +136,10 @@ subroutine timemanager itime=0 !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc write(*,46) float(itime)/3600,itime,numpart - if (verbosity.gt.0) then - write (*,*) 'timemanager> starting simulation' - if (verbosity.gt.1) then - CALL SYSTEM_CLOCK(count_clock) - WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate) - endif + call verboutput(verbosity,1,'timemanager> starting simulation') + if (verbosity.gt.1) then + CALL SYSTEM_CLOCK(count_clock) + WRITE(*,*) 'timemanager> SYSTEM CLOCK', (count_clock-count_clock0) / real(count_rate) endif do itime=0,ideltas,lsynctime diff --git a/src/verboutput.f90 b/src/verboutput.f90 new file mode 100644 index 0000000000000000000000000000000000000000..74e9d84b5a10a73ac3bc7dc9f6bced579da89e4e --- /dev/null +++ b/src/verboutput.f90 @@ -0,0 +1,55 @@ +!********************************************************************** +! Copyright 1998-2015 * +! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * +! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * +! * +! This file is part of FLEXPART. * +! * +! FLEXPART is free software: you can redistribute it and/or modify * +! it under the terms of the GNU General Public License as published by* +! the Free Software Foundation, either version 3 of the License, or * +! (at your option) any later version. * +! * +! FLEXPART is distributed in the hope that it will be useful, * +! but WITHOUT ANY WARRANTY; without even the implied warranty of * +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +! GNU General Public License for more details. * +! * +! You should have received a copy of the GNU General Public License * +! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * +!********************************************************************** + + subroutine verboutput (iverbosity, iverbositymin, infostring) + + !********************************************************************** + ! * + ! SUBROUTINE VERBOUTPUT * + ! * + !********************************************************************** + ! * + ! AUTHOR: PETRA sEIBERT * + ! DATE: 2015-03-07 * + ! LAST UPDATE: 2000-10-17, Andreas Stohl * + ! * + !********************************************************************** + ! Changes: + ! + !********************************************************************** + ! * + ! DESCRIPTION: * + ! * + ! write debugging information to standard output + ! depending on the verbosity level + ! + !********************************************************************** + + implicit none + + integer :: iverbosity ! verbosity level of the current run + integer :: iverbositymin ! min. verbosity level required to produce output + character(*) :: infostring + + if (iverbosity .ge. iverbositymin) print*, '>'//infostring + + end subroutine verboutput + diff --git a/src/verttransform.f90 b/src/verttransform.f90 index 650af3a571decfbf83996ef392d7913c8b4c2b91..95ebb4bd2471fcb593f78d7d266fbd844aee5671 100644 --- a/src/verttransform.f90 +++ b/src/verttransform.f90 @@ -1,5 +1,5 @@ !********************************************************************** -! 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 * ! * @@ -49,10 +49,11 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh) ! Sabine Eckhardt, March 2007 ! added the variable cloud for use with scavenging - descr. in com_mod !***************************************************************************** + ! Petra Seibert, Feb 2015: Add quick fix from 2013 + !***************************************************************************** ! * ! Variables: * ! nx,ny,nz field dimensions in x,y and z direction * - ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition * ! uu(0:nxmax,0:nymax,nzmax,2) wind components in x-direction [m/s] * ! vv(0:nxmax,0:nymax,nzmax,2) wind components in y-direction [m/s] * ! ww(0:nxmax,0:nymax,nzmax,2) wind components in z-direction [deltaeta/s]* @@ -68,10 +69,10 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh) implicit none - integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym - integer :: rain_cloud_above,kz_inv + integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym,kz_inv + integer :: icloudtop real :: f_qvsat,pressure - real :: rh,lsp,convp + real :: rh,lsp,convp,prec,rhmin real :: rhoh(nuvzmax),pinmconv(nzmax) real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi real :: xlon,ylat,xlonr,dzdx,dzdy @@ -84,7 +85,7 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh) real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax) real,parameter :: const=r_air/ga - logical :: init = .true. + logical :: init = .true., lconvectprec, lsearch !************************************************************************* @@ -476,44 +477,79 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh) endif - !write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz - ! create a cloud and rainout/washout field, clouds occur where rh>80% - ! total cloudheight is stored at level 0 + ! write (*,*) 'diagnosing clouds, n:',n,nymin1,nxmin1,nz + jy_loop: & do jy=0,nymin1 do ix=0,nxmin1 - rain_cloud_above=0 + lsp=lsprec(ix,jy,1,n) convp=convprec(ix,jy,1,n) - cloudsh(ix,jy,n)=0 - do kz_inv=1,nz-1 - kz=nz-kz_inv+1 - pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n) - rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n)) - clouds(ix,jy,kz,n)=0 - if (rh.gt.0.8) then ! in cloud - if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation - rain_cloud_above=1 - cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ & - height(kz)-height(kz-1) - if (lsp.ge.convp) then - clouds(ix,jy,kz,n)=3 ! lsp dominated rainout - else - clouds(ix,jy,kz,n)=2 ! convp dominated rainout - endif - else ! no precipitation - clouds(ix,jy,kz,n)=1 ! cloud - endif - else ! no cloud - if (rain_cloud_above.eq.1) then ! scavenging - if (lsp.ge.convp) then - clouds(ix,jy,kz,n)=5 ! lsp dominated washout - else - clouds(ix,jy,kz,n)=4 ! convp dominated washout - endif + prec=lsp+convp + if (lsp .gt. convp) then ! prectype='lsp' + lconvectprec = .false. + else ! prectype='cp ' + lconvectprec = .true. + endif + icloudbot(ix,jy,n)=icmv + icloudtop=icmv ! this is just a local variable + rhmin=rhmininit ! just initialise in a way that cond is true + lsearch=.true. + + cloudsearch_loop: & + do while (rhmin .ge. rhminx .and. lsearch) + ! give up for < rhminx + + kz_loop: & + do kz_inv=1,nz-1 + kz=nz-kz_inv+1 + pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n) + rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n)) + if (rh .gt. rhmin) then + if (icloudbot(ix,jy,n) .eq. icmv) then + icloudtop=nint(height(kz)) ! use int to save memory endif - endif - end do - end do - end do + icloudbot(ix,jy,n)=nint(height(kz)) + endif + end do kz_loop + + ! PS try to get a cloud thicker than 50 m + ! PS if there is at least precmin mm/h + if (prec .gt. precmin .and. & + ( icloudbot(ix,jy,n) .eq. icmv .or. & + icloudtop-icloudbot(ix,jy,n) .lt. 50)) then + rhmin = rhmin - 0.05 + else + lsearch = .false. + endif + + enddo cloudsearch_loop + + ! PS implement a rough fix for badly represented convection + ! PS is based on looking at a limited set of comparison data + if (lconvectprec .and. icloudtop .lt. icloudtopconvmin .or. & + icloudbot(ix,jy,n) .lt. icloudtopmin .and. prec .gt. precmin) then + if (convp .lt. 0.1) then + icloudbot(ix,jy,n) = icloudbot1 + icloudtop = icloudtop1 + else + icloudbot(ix,jy,n) = icloudbot2 + icloudtop = icloudtop2 + endif + endif + + if (icloudtop .ne. icmv) then + icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n) + else + icloudthck(ix,jy,n) = icmv ! no cloud found whatsoever + endif + + ! PS get rid of too thin clouds + if (icloudthck(ix,jy,n) .lt. 50) then + icloudbot(ix,jy,n)=icmv + icloudthck(ix,jy,n)=icmv + endif + + enddo + enddo jy_loop end subroutine verttransform diff --git a/src/verttransform_nests.f90 b/src/verttransform_nests.f90 index 07fcd35d00c29fe24bd464c14dc1accd02de746e..d1cdc9b4a1b550e0e5a00b4a0436352551992181 100644 --- a/src/verttransform_nests.f90 +++ b/src/verttransform_nests.f90 @@ -1,5 +1,5 @@ !********************************************************************** -! 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 * ! * @@ -51,6 +51,8 @@ subroutine verttransform_nests(n,uuhn,vvhn,wwhn,pvhn) ! Sabine Eckhardt, March 2007 ! add the variable cloud for use with scavenging - descr. in com_mod !***************************************************************************** + ! Petra Seibert, Feb 2015: Add quick fix from 2013 + !***************************************************************************** ! * ! Variables: * ! nxn,nyn,nuvz,nwz field dimensions in x,y and z direction * @@ -68,24 +70,27 @@ subroutine verttransform_nests(n,uuhn,vvhn,wwhn,pvhn) implicit none - integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp - integer :: rain_cloud_above,kz_inv - real :: f_qvsat,pressure,rh,lsp,convp - real :: wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax) - real :: uvwzlev(0:nxmaxn-1,0:nymaxn-1,nzmax) - real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi,cosf + integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp,kz_inv + integer :: icloudtop + real :: f_qvsat,pressure + real :: rh,lsp,convp,prec,rhmin + real :: rhoh(nuvzmax),pinmconv(nzmax) + real :: wzlev(nwzmax),uvwzlev(0:nxmaxn-1,0:nymaxn-1,nzmax) + real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi real :: dzdx,dzdy - real :: dzdx1,dzdx2,dzdy1,dzdy2 + real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests) real,parameter :: const=r_air/ga + logical :: lconvectprec, lsearch ! Loop over all nests !******************** + numbnest_loop: & do l=1,numbnests ! Loop over the whole grid @@ -142,8 +147,8 @@ subroutine verttransform_nests(n,uuhn,vvhn,wwhn,pvhn) (aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l))) - ! Levels, where u,v,t and q are given - !************************************ + ! Levels where u,v,t and q are given + !*********************************** uun(ix,jy,1,n,l)=uuhn(ix,jy,1,l) vvn(ix,jy,1,n,l)=vvhn(ix,jy,1,l) @@ -193,8 +198,8 @@ subroutine verttransform_nests(n,uuhn,vvhn,wwhn,pvhn) end do - ! Levels, where w is given - !************************* + ! Levels where w is given + !************************ wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1) wwn(ix,jy,nz,n,l)=wwhn(ix,jy,nwz,l)*pinmconv(nz) @@ -274,50 +279,82 @@ subroutine verttransform_nests(n,uuhn,vvhn,wwhn,pvhn) wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*ui+dzdy*vi) end do - end do end do - - !write (*,*) 'initializing nested cloudsn, n:',n - ! create a cloud and rainout/washout field, cloudsn occur where rh>80% + ! write (*,*) 'diagnosing nested clouds, n:',n + jy_loop: & do jy=0,nyn(l)-1 do ix=0,nxn(l)-1 - rain_cloud_above=0 + lsp=lsprecn(ix,jy,1,n,l) convp=convprecn(ix,jy,1,n,l) - cloudsnh(ix,jy,n,l)=0 - do kz_inv=1,nz-1 - kz=nz-kz_inv+1 - pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l) - rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l)) - cloudsn(ix,jy,kz,n,l)=0 - if (rh.gt.0.8) then ! in cloud - if ((lsp.gt.0.01).or.(convp.gt.0.01)) then - rain_cloud_above=1 - cloudsnh(ix,jy,n,l)=cloudsnh(ix,jy,n,l)+ & - height(kz)-height(kz-1) - if (lsp.ge.convp) then - cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout - else - cloudsn(ix,jy,kz,n,l)=2 ! convp dominated rainout - endif - else ! no precipitation - cloudsn(ix,jy,kz,n,l)=1 ! cloud - endif - else ! no cloud - if (rain_cloud_above.eq.1) then ! scavenging - if (lsp.ge.convp) then - cloudsn(ix,jy,kz,n,l)=5 ! lsp dominated washout - else - cloudsn(ix,jy,kz,n,l)=4 ! convp dominated washout - endif + prec=lsp+convp + if (lsp .gt. convp) then ! prectype='lsp' + lconvectprec = .false. + else ! prectype='cp ' + lconvectprec = .true. + endif + icloudbotn(ix,jy,n,l)=icmv + icloudtop=icmv ! this is just a local variable + rhmin=rhmininit + lsearch=.true. + + cloudsearch_loop: & + do while (rhmin .ge. rhminx .and. lsearch) + ! give up for < rhminx + + kz_loop: & + do kz_inv=1,nz-1 + kz=nz-kz_inv+1 + pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l) + rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l)) + if (rh .gt. rhmin) then + if (icloudbotn(ix,jy,n,l) .eq. icmv) then + icloudtop=nint(height(kz)) ! use int to save memory endif - endif - end do - end do - end do - - end do + icloudbotn(ix,jy,n,l)=nint(height(kz)) + endif + end do kz_loop + + ! PS try to get a cloud thicker than 50 m + ! PS if there is at least precmin mm/h + if (prec .gt. precmin .and. & + ( icloudbotn(ix,jy,n,l) .eq. icmv .or. & + icloudtop-icloudbotn(ix,jy,n,l) .lt. 50)) then + rhmin = rhmin - 0.05 + else + lsearch = .false. + endif + + enddo cloudsearch_loop + + ! PS implement a rough fix for badly represented convection + ! PS is based on looking at a limited set of comparison data + if (lconvectprec .and. icloudtop .lt. icloudtopconvmin .or. & + icloudbotn(ix,jy,n,l) .lt. icloudtopmin .and. prec .gt. precmin) then + if (convp .lt. 0.1) then + icloudbotn(ix,jy,n,l) = icloudbot1 + icloudtop = icloudtop1 + else + icloudbotn(ix,jy,n,l) = icloudbot2 + icloudtop = icloudtop2 + endif + endif + if (icloudtop .ne. icmv) then + icloudthckn(ix,jy,n,l) = icloudtop-icloudbotn(ix,jy,n,l) + else + icloudthckn(ix,jy,n,l) = icmv + endif + ! PS get rid of too thin clouds + if (icloudthck(ix,jy,n) .lt. 50) then + icloudbotn(ix,jy,n,l)=icmv + icloudthckn(ix,jy,n,l)=icmv + endif + + enddo + enddo jy_loop + + enddo numbnest_loop end subroutine verttransform_nests diff --git a/src/wetdepo.f90 b/src/wetdepo.f90 index 7f99daa19c492cd627cd55ef7df42b84536659ed..0084826c5dcea3465be6ad0f282dcece0443d347 100644 --- a/src/wetdepo.f90 +++ b/src/wetdepo.f90 @@ -1,5 +1,5 @@ !********************************************************************** -! 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 * ! * @@ -38,6 +38,11 @@ subroutine wetdepo(itime,ltsample,loutnext) ! use centred precipitation data for integration * ! 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 + ! * !***************************************************************************** ! * ! Variables: * @@ -70,11 +75,10 @@ subroutine wetdepo(itime,ltsample,loutnext) implicit none 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 save lfr,cfr @@ -96,14 +100,16 @@ subroutine wetdepo(itime,ltsample,loutnext) ! Loop over all particles !************************ + 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)) do nage=1,nageclass @@ -122,7 +128,7 @@ subroutine wetdepo(itime,ltsample,loutnext) ngrid=j goto 23 endif - end do + enddo 23 continue @@ -130,8 +136,8 @@ subroutine wetdepo(itime,ltsample,loutnext) !********************************** 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) else @@ -147,43 +153,56 @@ subroutine wetdepo(itime,ltsample,loutnext) 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) + 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) + nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn, & + 1,memtime(1),memtime(2),interp_time,lsp,convp,cc,icbot,ictop) endif - if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20 +! 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) + ! 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 - clouds_v=cloudsn(ix,jy,hz,n,ngrid) - clouds_h=cloudsnh(ix,jy,n,ngrid) + goto 20 ! above cloud 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' + ! 1) Parameterization of the the area fraction of the grid cell where the ! precipitation occurs: the absolute limit is the total cloud cover, but @@ -227,67 +246,56 @@ subroutine wetdepo(itime,ltsample,loutnext) ! Computation of wet deposition !********************************************************** + 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 + 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 - !****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 + elseif (indcloud .eq. 2) then ! in-cloud scavenging - - !******************** - ! 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 (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 ! weta_in=2.0E-07 (default) ! 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, - - !************************************************** - ! 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 - +! 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 + wetdeposit(ks)=xmass1(jpart,ks)* & + (1.-exp(-wetscav*abs(ltsample)))*grfraction restmass = xmass1(jpart,ks)-wetdeposit(ks) if (ioutputforeachrelease.eq.1) then kp=npoint(jpart) @@ -296,34 +304,36 @@ subroutine wetdepo(itime,ltsample,loutnext) endif 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. endif ! 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 + else ! weta(k) <= 0 + wetdeposit(ks)=0. + + endif - end do !all species + enddo species_loop - ! 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 - !***************************************************************************** + ! 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) + wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),itage,nage,kp) endif -20 continue - end do ! all particles +20 continue ! jump here for particles not to be treated + enddo particle_loop end subroutine wetdepo diff --git a/src/wetdepokernel.f90 b/src/wetdepokernel.f90 index e90fca67ddae51b666a149a5d7e9105fa3fca103..9075b357d48b3dfe02547dd79f8e4a5606f5ea70 100644 --- a/src/wetdepokernel.f90 +++ b/src/wetdepokernel.f90 @@ -1,5 +1,5 @@ !********************************************************************** -! 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 * ! * @@ -19,16 +19,20 @@ ! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * !********************************************************************** -subroutine wetdepokernel(nunc,deposit,x,y,nage,kp) - ! i i i i i +subroutine wetdepokernel(nunc,deposit,x,y,itage,nage,kp) + ! i i i i i i !***************************************************************************** ! * ! Attribution of the deposition from an individual particle to the * - ! deposition fields using a uniform kernel with bandwidths dxout and dyout.* + ! deposition fields using a uniform kernel with bandwidths dxout * + ! and dyout. * ! * ! Author: A. Stohl * ! * ! 26 December 1996 * + ! + ! PS, 2/2015: do not use kernel for itage < 3 h + ! same as for concentration in conccalc.f90 ! * !***************************************************************************** ! * @@ -47,7 +51,7 @@ subroutine wetdepokernel(nunc,deposit,x,y,nage,kp) implicit none real :: x,y,deposit(maxspec),ddx,ddy,xl,yl,wx,wy,w - integer :: ix,jy,ixp,jyp,nunc,nage,ks,kp + integer :: ix,jy,ixp,jyp,nunc,nage,ks,kp,itage xl=(x*dx+xoutshift)/dxout yl=(y*dy+youtshift)/dyout @@ -73,38 +77,51 @@ subroutine wetdepokernel(nunc,deposit,x,y,nage,kp) endif - ! Determine mass fractions for four grid points - !********************************************** + if (itage .lt. itagekernmin) then + ! no kernel, direct attribution to grid cell + + do ks=1,nspec + if (ix.ge.0 .and. jy.ge.0 .and. & + ix.le.numxgrid-1 .and. jy.le.numygrid-1) & + wetgridunc(ix,jy,ks,kp,nunc,nage)= & + wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks) + enddo + + else + ! Determine mass fractions for four grid points & distribute - do ks=1,nspec + do ks=1,nspec - if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. & - (jy.le.numygrid-1)) then - w=wx*wy - wetgridunc(ix,jy,ks,kp,nunc,nage)= & - wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w - endif + if (ix.ge.0 .and. jy.ge.0 .and. & + ix.le.numxgrid-1 .and. jy.le.numygrid-1) then + w=wx*wy + wetgridunc(ix,jy,ks,kp,nunc,nage)= & + wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w + endif - if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. & - (jyp.le.numygrid-1)) then - w=(1.-wx)*(1.-wy) - wetgridunc(ixp,jyp,ks,kp,nunc,nage)= & - wetgridunc(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w - endif + if ((ixp.ge.0).and.(jyp.ge.0).and.& + (ixp.le.numxgrid-1).and. (jyp.le.numygrid-1)) then + w=(1.-wx)*(1.-wy) + wetgridunc(ixp,jyp,ks,kp,nunc,nage)= & + wetgridunc(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w + endif - if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgrid-1).and. & - (jy.le.numygrid-1)) then - w=(1.-wx)*wy - wetgridunc(ixp,jy,ks,kp,nunc,nage)= & + if ((ixp.ge.0).and.(jy.ge.0).and. & + (ixp.le.numxgrid-1).and.(jy.le.numygrid-1)) then + w=(1.-wx)*wy + wetgridunc(ixp,jy,ks,kp,nunc,nage)= & wetgridunc(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w - endif + endif + + if ((ix.ge.0).and.(jyp.ge.0).and. & + (ix.le.numxgrid-1).and.(jyp.le.numygrid-1)) then + w=wx*(1.-wy) + wetgridunc(ix,jyp,ks,kp,nunc,nage)= & + wetgridunc(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w + endif + + end do - if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgrid-1).and. & - (jyp.le.numygrid-1)) then - w=wx*(1.-wy) - wetgridunc(ix,jyp,ks,kp,nunc,nage)= & - wetgridunc(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w endif - end do end subroutine wetdepokernel diff --git a/src/wetdepokernel_nest.f90 b/src/wetdepokernel_nest.f90 index f3aea688452f87d4ff85b0c920c86fbaf97faa15..a187720f5d5e1595617068168f5e6c08be38400b 100644 --- a/src/wetdepokernel_nest.f90 +++ b/src/wetdepokernel_nest.f90 @@ -1,5 +1,5 @@ !********************************************************************** -! 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 * ! * @@ -19,9 +19,8 @@ ! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * !********************************************************************** -subroutine wetdepokernel_nest & - (nunc,deposit,x,y,nage,kp) - ! i i i i i i +subroutine wetdepokernel_nest (nunc,deposit,x,y,itage,nage,kp) + ! i i i i i i i !***************************************************************************** ! * ! Attribution of the deposition from an individual particle to the * @@ -34,6 +33,9 @@ subroutine wetdepokernel_nest & ! * ! 2 September 2004: Adaptation from wetdepokernel. * ! * + ! + ! PS, 2/2015: do not use kernel for itage < 3 h + ! same as for concentration in conccalc.f90 ! * !***************************************************************************** ! * @@ -52,7 +54,7 @@ subroutine wetdepokernel_nest & implicit none real :: x,y,deposit(maxspec),ddx,ddy,xl,yl,wx,wy,w - integer :: ix,jy,ixp,jyp,ks,kp,nunc,nage + integer :: ix,jy,ixp,jyp,ks,kp,nunc,nage,itage @@ -79,39 +81,50 @@ subroutine wetdepokernel_nest & wy=0.5+ddy endif + if (itage .lt. itagekernmin) then + ! no kernel, direct attribution to grid cell + + do ks=1,nspec + if (ix.ge.0 .and. jy.ge.0 .and. & + ix.le.numxgridn-1 .and. jy.le.numygridn-1) & + wetgriduncn(ix,jy,ks,kp,nunc,nage)= & + wetgriduncn(ix,jy,ks,kp,nunc,nage)+deposit(ks) + enddo + + else + ! Determine mass fractions for four grid points & distribute + + do ks=1,nspec + + if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. & + (jy.le.numygridn-1)) then + w=wx*wy + wetgriduncn(ix,jy,ks,kp,nunc,nage)= & + wetgriduncn(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w + endif + + if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgridn-1).and. & + (jyp.le.numygridn-1)) then + w=(1.-wx)*(1.-wy) + wetgriduncn(ixp,jyp,ks,kp,nunc,nage)= & + wetgriduncn(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w + endif + + if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgridn-1).and. & + (jy.le.numygridn-1)) then + w=(1.-wx)*wy + wetgriduncn(ixp,jy,ks,kp,nunc,nage)= & + wetgriduncn(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w + endif + + if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgridn-1).and. & + (jyp.le.numygridn-1)) then + w=wx*(1.-wy) + wetgriduncn(ix,jyp,ks,kp,nunc,nage)= & + wetgriduncn(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w + endif + + end do - ! Determine mass fractions for four grid points - !********************************************** - - do ks=1,nspec - - if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. & - (jy.le.numygridn-1)) then - w=wx*wy - wetgriduncn(ix,jy,ks,kp,nunc,nage)= & - wetgriduncn(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w - endif - - if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgridn-1).and. & - (jyp.le.numygridn-1)) then - w=(1.-wx)*(1.-wy) - wetgriduncn(ixp,jyp,ks,kp,nunc,nage)= & - wetgriduncn(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w - endif - - if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgridn-1).and. & - (jy.le.numygridn-1)) then - w=(1.-wx)*wy - wetgriduncn(ixp,jy,ks,kp,nunc,nage)= & - wetgriduncn(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w - endif - - if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgridn-1).and. & - (jyp.le.numygridn-1)) then - w=wx*(1.-wy) - wetgriduncn(ix,jyp,ks,kp,nunc,nage)= & - wetgriduncn(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w endif - - end do end subroutine wetdepokernel_nest