From 3c1da523f53497adfac8379a95ba9e93216ad2c3 Mon Sep 17 00:00:00 2001
From: Ignacio Pisso <Ignacio.Pisso@nilu.no>
Date: Wed, 15 Jul 2015 11:23:56 +0200
Subject: [PATCH] from svn branch petra r37 | pesei | 2015-04-22 17:42:35 +0200
(Wed, 22 Apr 2015) | 1 line Changed paths: M
/branches/petra/src/FLEXPART.f90 M /branches/petra/src/com_mod.f90 M
/branches/petra/src/getfields.f90 M /branches/petra/src/gridcheck.f90 M
/branches/petra/src/interpol_rain.f90 M
/branches/petra/src/interpol_rain_nests.f90 M
/branches/petra/src/par_mod.f90 M /branches/petra/src/readcommand.f90 M
/branches/petra/src/readspecies.f90 M /branches/petra/src/readwind.f90
M /branches/petra/src/readwind_nests.f90 M
/branches/petra/src/timemanager.f90 A /branches/petra/src/verboutput.f90
M /branches/petra/src/verttransform.f90 M
/branches/petra/src/verttransform_nests.f90 M
/branches/petra/src/wetdepo.f90 M /branches/petra/src/wetdepokernel.f90
M /branches/petra/src/wetdepokernel_nest.f90
Wet dep quick fix and other small changes. Wet depo quick fix not final yet.
---
src/FLEXPART.f90 | 2 +-
src/com_mod.f90 | 15 +--
src/getfields.f90 | 24 ++++
src/gridcheck.f90 | 189 ++++++++++++++++--------------
src/interpol_rain.f90 | 81 +++++++++++--
src/interpol_rain_nests.f90 | 92 ++++++++++++---
src/par_mod.f90 | 21 +++-
src/readcommand.f90 | 27 ++---
src/readspecies.f90 | 218 +++++++++++++++++-----------------
src/readwind.f90 | 207 ++++++++++++++++++---------------
src/readwind_nests.f90 | 226 +++++++++++++++++++-----------------
src/timemanager.f90 | 10 +-
src/verboutput.f90 | 55 +++++++++
src/verttransform.f90 | 116 +++++++++++-------
src/verttransform_nests.f90 | 135 +++++++++++++--------
src/wetdepo.f90 | 174 ++++++++++++++-------------
src/wetdepokernel.f90 | 79 ++++++++-----
src/wetdepokernel_nest.f90 | 89 ++++++++------
18 files changed, 1063 insertions(+), 697 deletions(-)
create mode 100644 src/verboutput.f90
diff --git a/src/FLEXPART.f90 b/src/FLEXPART.f90
index 40c5722d..c2d9aabb 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 cd61d7a1..bc6ad539 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 fad4e969..38a4ed85 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 32b11fe0..72c9e4fc 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 e9d71c5a..c6b13856 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 53e64b88..7d2e9437 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 df5d9882..edc7bcf8 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 8941a792..e76d2a8c 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 820aea54..a4cfc8e4 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 b3c08a1f..039b9765 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 ada53329..9865d588 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 1d70ce48..04b50191 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 00000000..74e9d84b
--- /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 650af3a5..95ebb4bd 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 07fcd35d..d1cdc9b4 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 7f99daa1..0084826c 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 e90fca67..9075b357 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 f3aea688..a187720f 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
--
GitLab