Commit 07c3e71f authored by Ignacio Pisso's avatar Ignacio Pisso
Browse files

check dsigma for dry deposition velocity

parent 95a45d34
...@@ -61,7 +61,7 @@ program flexpart ...@@ -61,7 +61,7 @@ program flexpart
implicit none implicit none
integer :: i,j,ix,jy,inest integer :: i,j,ix,jy,inest, iopt
integer :: idummy = -320 integer :: idummy = -320
character(len=256) :: inline_options !pathfile, flexversion, arg2 character(len=256) :: inline_options !pathfile, flexversion, arg2
integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
...@@ -79,7 +79,7 @@ program flexpart ...@@ -79,7 +79,7 @@ program flexpart
! FLEXPART version string ! FLEXPART version string
flexversion_major = '10' ! Major version number, also used for species file names flexversion_major = '10' ! Major version number, also used for species file names
flexversion='Version '//trim(flexversion_major)//'.4 (2019-07-16)' flexversion='Version '//trim(flexversion_major)//'.4 (2019-07-23)'
verbosity=0 verbosity=0
! Read the pathnames where input/output files are stored ! Read the pathnames where input/output files are stored
...@@ -108,15 +108,43 @@ program flexpart ...@@ -108,15 +108,43 @@ program flexpart
print*,'Welcome to FLEXPART ', trim(flexversion) print*,'Welcome to FLEXPART ', trim(flexversion)
print*,'FLEXPART is free software released under the GNU General Public License.' print*,'FLEXPART is free software released under the GNU General Public License.'
! Ingest inline options
!*******************************************************
if (inline_options(1:1).eq.'-') then if (inline_options(1:1).eq.'-') then
if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then print*,'inline_options:',inline_options
print*, 'Verbose mode 1: display detailed information during run' !verbose mode
verbosity=1 iopt=index(inline_options,'v')
if (iopt.gt.0) then
verbosity=1
!print*, iopt, inline_options(iopt+1:iopt+1)
if (trim(inline_options(iopt+1:iopt+1)).eq.'2') then
!print*, 'verbosity=2'
print*, 'Verbose mode 2: display more detailed information during run'
verbosity=2
endif
endif endif
if (trim(inline_options).eq.'-v2') then
print*, 'Verbose mode 2: display more detailed information during run'
verbosity=2 !debug mode
iopt=index(inline_options,'d')
if (iopt.gt.0) then
debug_mode=.true.
endif
endif endif
! stop
! if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
! print*, 'Verbose mode 1: display detailed information during run'
! verbosity=1
! endif
! if (trim(inline_options).eq.'-v2') then
! print*, 'Verbose mode 2: display more detailed information during run'
! verbosity=2
! endif
if (trim(inline_options).eq.'-i') then if (trim(inline_options).eq.'-i') then
print*, 'Info mode: provide detailed run specific information and stop' print*, 'Info mode: provide detailed run specific information and stop'
verbosity=1 verbosity=1
...@@ -134,6 +162,9 @@ program flexpart ...@@ -134,6 +162,9 @@ program flexpart
print*, 'nymax=',nymax print*, 'nymax=',nymax
print*, 'nzmax=',nzmax print*, 'nzmax=',nzmax
print*,'nxshift=',nxshift print*,'nxshift=',nxshift
endif
if (verbosity.gt.0) then
write(*,*) 'call readpaths' write(*,*) 'call readpaths'
endif endif
call readpaths call readpaths
...@@ -192,6 +223,9 @@ program flexpart ...@@ -192,6 +223,9 @@ program flexpart
! Detect metdata format ! Detect metdata format
!********************** !**********************
if (verbosity.gt.0) then
write(*,*) 'call detectformat'
endif
metdata_format = detectformat() metdata_format = detectformat()
......
...@@ -761,6 +761,7 @@ module com_mod ...@@ -761,6 +761,7 @@ module com_mod
!******************** !********************
! Verbosity, testing flags, namelist I/O ! Verbosity, testing flags, namelist I/O
!******************** !********************
logical :: debug_mode=.false.
integer :: verbosity=0 integer :: verbosity=0
integer :: info_flag=0 integer :: info_flag=0
integer :: count_clock, count_clock0, count_rate, count_max integer :: count_clock, count_clock0, count_rate, count_max
......
...@@ -169,9 +169,6 @@ subroutine getvdep(n,ix,jy,ust,temp,pa,L,gr,rh,rr,snow,vdepo) ...@@ -169,9 +169,6 @@ subroutine getvdep(n,ix,jy,ust,temp,pa,L,gr,rh,rr,snow,vdepo)
if (reldiff(i).gt.0.) then if (reldiff(i).gt.0.) then
if ((ra+rb(i)+rc(i)).gt.0.) then if ((ra+rb(i)+rc(i)).gt.0.) then
vd=1./(ra+rb(i)+rc(i)) vd=1./(ra+rb(i)+rc(i))
! XXXXXXXXXXXXXXXXXXXXXXXXXX TEST
! vd=1./rc(i)
! XXXXXXXXXXXXXXXXXXXXXXXXXX TEST
else else
vd=9.999 vd=9.999
endif endif
...@@ -187,6 +184,10 @@ subroutine getvdep(n,ix,jy,ust,temp,pa,L,gr,rh,rr,snow,vdepo) ...@@ -187,6 +184,10 @@ subroutine getvdep(n,ix,jy,ust,temp,pa,L,gr,rh,rr,snow,vdepo)
call partdep(nspec,density,fract,schmi,vset,raquer,ust,nyl,vdepo) call partdep(nspec,density,fract,schmi,vset,raquer,ust,nyl,vdepo)
!if (debug_mode) then
! print*,'getvdep:188: vdepo=', vdepo
!stop
!endif
! 7. If no detailed parameterization available, take constant deposition ! 7. If no detailed parameterization available, take constant deposition
! velocity if that is available ! velocity if that is available
......
...@@ -102,19 +102,20 @@ subroutine part0(dquer,dsigma,density,fract,schmi,cun,vsh) ...@@ -102,19 +102,20 @@ subroutine part0(dquer,dsigma,density,fract,schmi,cun,vsh)
d01=dquer*dsigma**(-3.+delta*real(i)) d01=dquer*dsigma**(-3.+delta*real(i))
x01=alog(d01/dquer)/xdummy x01=alog(d01/dquer)/xdummy
x02=alog(d02/dquer)/xdummy x02=alog(d02/dquer)/xdummy
!print*,'part0:: d02=' , d02 , 'd01=', d01
! Area under Gauss-function is calculated and gives mass fraction of interval ! Area under Gauss-function is calculated and gives mass fraction of interval
!**************************************************************************** !****************************************************************************
fract(i)=0.5*(erf(x01)-erf(x02)) fract(i)=0.5*(erf(x01)-erf(x02))
!print*,'part0:: fract(',i,')', fract(i)
!print*,'part0:: fract', fract(i), x01, x02, erf(x01), erf(x02)
! Geometric mean diameter of interval in [m] ! Geometric mean diameter of interval in [m]
!******************************************* !*******************************************
dmean=1.E-6*exp(0.5*alog(d01*d02)) dmean=1.E-6*exp(0.5*alog(d01*d02))
!print*,'part0:: dmean=', dmean
! Calculation of time independent parameters of each interval ! Calculation of time independent parameters of each interval
!************************************************************ !************************************************************
...@@ -131,6 +132,10 @@ subroutine part0(dquer,dsigma,density,fract,schmi,cun,vsh) ...@@ -131,6 +132,10 @@ subroutine part0(dquer,dsigma,density,fract,schmi,cun,vsh)
schmi(i)=schmidt**(-2./3.) schmi(i)=schmidt**(-2./3.)
vsh(i)=ga*density*dmean*dmean*cun/(18.*myl) vsh(i)=ga*density*dmean*dmean*cun/(18.*myl)
!print*,'part0:: vsh(',i,')', vsh(i)
end do end do
!stop 'part0'
end subroutine part0 end subroutine part0
...@@ -71,6 +71,7 @@ subroutine partdep(nc,density,fract,schmi,vset,ra,ustar,nyl,vdep) ...@@ -71,6 +71,7 @@ subroutine partdep(nc,density,fract,schmi,vset,ra,ustar,nyl,vdep)
!***************************************************************************** !*****************************************************************************
use par_mod use par_mod
use com_mod, only: debug_mode
implicit none implicit none
...@@ -109,8 +110,19 @@ subroutine partdep(nc,density,fract,schmi,vset,ra,ustar,nyl,vdep) ...@@ -109,8 +110,19 @@ subroutine partdep(nc,density,fract,schmi,vset,ra,ustar,nyl,vdep)
!*********************************************************************** !***********************************************************************
vdep(ic)=vdep(ic)+vdepj*fract(ic,j) vdep(ic)=vdep(ic)+vdepj*fract(ic,j)
!print*, 'partdep:113: ic', ic, 'vdep', vdep
!stop
end do end do
endif endif
end do end do
!if (debug_mode) then
! print*, 'partdep:122:'
! write(*,*) (vdep(ic), ic=1,nc)
!stop
!endif
end subroutine partdep end subroutine partdep
...@@ -295,7 +295,13 @@ subroutine readspecies(id_spec,pos_spec) ...@@ -295,7 +295,13 @@ subroutine readspecies(id_spec,pos_spec)
end if end if
if (density(pos_spec) .gt. 0) then if (density(pos_spec) .gt. 0) then
write(*,'(a)') ' Dry deposition is turned : ON' write(*,'(a)') ' Dry deposition is turned : ON'
if (reldiff(pos_spec).gt.0) then
stop 'density>0 (SPECIES is a particle) implies reldiff <=0 '
endif
else else
if (reldiff(pos_spec).le.0) then
stop 'density<=0 (SPECIES is a gas) implies reldiff >0 '
endif
write(*,'(a)') ' Dry deposition is (density<0) : OFF' write(*,'(a)') ' Dry deposition is (density<0) : OFF'
end if end if
if (crain_aero(pos_spec).gt.10.0 .or. csnow_aero(pos_spec).gt.10.0 .or. & if (crain_aero(pos_spec).gt.10.0 .or. csnow_aero(pos_spec).gt.10.0 .or. &
...@@ -343,7 +349,21 @@ subroutine readspecies(id_spec,pos_spec) ...@@ -343,7 +349,21 @@ subroutine readspecies(id_spec,pos_spec)
end if end if
end if end if
if (dsigma(i).eq.0.) dsigma(i)=1.0001 ! avoid floating exception ! if (dsigma(i).eq.0.) dsigma(i)=1.0001 ! avoid floating exception
if (dsigma(i).le.1.) then !dsigma(i)=1.0001 ! avoid floating exception
!write(*,*) '#### FLEXPART MODEL ERROR! ####'
write(*,*) '#### FLEXPART MODEL WARNING ####'
write(*,*) '#### in SPECIES_',aspecnumb, ' ####'
write(*,*) '#### from v10.4 dsigma has to be larger than 1 ####'
write(*,*) '#### to adapt older SPECIES files, ####'
write(*,*) '#### if dsigma was < 1 ####'
write(*,*) '#### use the reciprocal of the old dsigma ####'
if (.not.debug_mode) then
stop
else
write(*,*) 'debug mode: continue'
endif
endif
if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES" ####' write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES" ####'
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment