Commit 2753a5c6 authored by Ignacio Pisso's avatar Ignacio Pisso

resolved merging conflicts with GFS branch

parents df96ea65 4ad96c54
......@@ -15,22 +15,22 @@
LOUTSAMPLE= 900, ! Interval of output sampling (s), higher stat. accuracy with shorter intervals
ITSPLIT= 99999999, ! Interval of particle splitting (s)
LSYNCTIME= 900, ! All processes are synchronized to this time interval (s)
CTL= -5.0000000, ! Times step smaller than Lagr. time scale, LSYNCTIME if CTL<0, but can be shorter
CTL= -5.0000000, ! CTL>1, ABL time step = (Lagrangian timescale (TL))/CTL, uses LSYNCTIME if CTL<0
IFINE= 4, ! Reduction for time step in vertical transport, used only if CTL>1
IOUT= 1, ! Output type: [1]mass 2]pptv 3]1&2 4]plume 5]1&4, +8 for NetCDF output
IPOUT= 0, ! Particle position output: 0]no 1]every output 2]only at end
LSUBGRID= 0, ! Increase of ABL heights due to sub-grid scale orographic variations;[0]off 1]on
LCONVECTION= 0, ! Switch for convection parameterization;[0]off 1]on
LCONVECTION= 1, ! Switch for convection parameterization;0]off [1]on
LAGESPECTRA= 0, ! Switch for calculation of age spectra (needs AGECLASSES);[0]off 1]on
IPIN= 0, ! Warm start from particle dump (needs previous partposit_end file); 0]no [1]yes
IPIN= 0, ! Warm start from particle dump (needs previous partposit_end file); [0]no 1]yes
IOUTPUTFOREACHRELEASE= 1, ! Separate output fields for each location in the RELEASE file; [0]no 1]yes
IFLUX= 0, ! Output of mass fluxes through output grid box boundaries
MDOMAINFILL= 0, ! Switch for domain-filling, if limited-area particles generated at boundary
IND_SOURCE= 1, ! Unit to be used at the ource ; [1]mass 2]mass mixing ratio
IND_RECEPTOR= 1, ! Unit to be used at the receptor; [1]mass 2]mass mixing ratio
IND_SOURCE= 1, ! Unit to be used at the source ; [1]mass 2]mass mixing ratio
IND_RECEPTOR= 1, ! Unit to be used at the receptor; [1]mass 2]mass mixing ratio 3]wet depo. 4]dry depo.
MQUASILAG= 0, ! Quasi-Lagrangian mode to track individual numbered particles
NESTED_OUTPUT= 0, ! Output also for a nested domain
LINIT_COND= 0, ! Output sensitivity to initial conditions (backwards mode only)
LINIT_COND= 0, ! Output sensitivity to initial conditions (bkw mode only) [0]off 1]conc 2]mmr
SURF_ONLY= 0, ! Output only for the lowest model layer, used w/ LINIT_COND=1 or 2
CBLFLAG= 0, ! Skewed, not Gaussian turbulence in the convective ABL, need large CTL and IFINE
OHFIELDS_PATH= "../../flexin/", ! Default path for OH file
......
*************************************************************************
* *
* *
* *
* Input file for the Lagrangian particle dispersion model FLEXPART *
* Please select your options *
* *
* *
* *
*************************************************************************
***************************************************************************************************************
* *
* *
* *
* Input file for the Lagrangian particle dispersion model FLEXPART *
* Please select your options *
* *
* *
* *
***************************************************************************************************************
&RELEASES_CTRL
NSPEC = 1,
SPECNUM_REL= 24,
NSPEC = 1, ! Total number of species
SPECNUM_REL= 24, ! Species numbers in directory SPECIES
/
&RELEASE
IDATE1 = 20120101,
ITIME1 = 090000,
IDATE2 = 20120101,
ITIME2 = 090000,
LON1 = 0.000,
LON2 = 0.000,
LAT1 = 20.000,
LAT2 = 20.000,
Z1 = 50.000,
Z2 = 50.000,
ZKIND = 1,
MASS = 1.0000E8,
PARTS = 10000,
COMMENT = "RELEASE 1",
&RELEASE ! For each release
IDATE1 = 20120101, ! Release start date, YYYYMMDD: YYYY=year, MM=month, DD=day
ITIME1 = 090000, ! Release start time in UTC HHMISS: HH hours, MI=minutes, SS=seconds
IDATE2 = 20120101, ! Release end date, same as IDATE1
ITIME2 = 090000, ! Release end time, same as ITIME1
LON1 = 0.000, ! Left longitude of release box -180 < LON1 <180
LON2 = 0.000, ! Right longitude of release box, same as LON1
LAT1 = 20.000, ! Lower latitude of release box, -90 < LAT1 < 90
LAT2 = 20.000, ! Upper latitude of release box same format as LAT1
Z1 = 50.000, ! Lower height of release box meters/hPa above reference level
Z2 = 50.000, ! Upper height of release box meters/hPa above reference level
ZKIND = 1, ! Reference level 1=above ground, 2=above sea level, 3 for pressure in hPa
MASS = 1.0000E8, ! Total mass emitted, only relevant for fwd simulations
PARTS = 10000, ! Total number of particles to be released
COMMENT = "RELEASE 1", ! Comment, written in the outputfile
/
......@@ -18,6 +18,8 @@
POHCCONST=-0.9E9, ! OH Reaction rate - C [cm^3/molecule/sec]
POHDCONST=-9.9, ! OH Reaction rate - D [K]
POHNCONST=2.0, ! OH Reaction rate - N (no unit)
PAREA_HOUR=0.578, 0.491, 0.428, 0.329, 0.384, 0.485, 0.763, 1.103, 1.084, 1.047, 1.096, 1.196, 1.298, 1.357, 1.447, 1.565, 1.636, 1.662, 1.401, 1.168, 1.031, 0.926, 0.816, 0.709
PAREA_DOW=1.060, 1.060, 1.060, 1.060, 1.060, 0.900, 0.000
PAREA_HOUR=0.578, 0.491, 0.428, 0.329, 0.384, 0.485, 0.763, 1.103, 1.084, 1.047, 1.096, 1.196, 1.298, 1.357, 1.447, 1.565, 1.636, 1.662, 1.401, 1.168, 1.031, 0.926, 0.816, 0.709,
PPOINT_HOUR=0.845, 0.806, 0.786, 0.779, 0.793, 0.832, 0.895, 0.977, 1.031, 1.071, 1.105, 1.118, 1.131, 1.136, 1.143, 1.141, 1.133, 1.118, 1.097, 1.091, 1.079, 1.036, 0.966, 0.892,
PAREA_DOW=1.060, 1.060, 1.060, 1.060, 1.060, 0.900, 0.000,
PPOINT_DOW=1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000,
/
......@@ -454,16 +454,21 @@ program flexpart
call timemanager(metdata_format)
if (verbosity.gt.0) then
! NIK 16.02.2005
do i=1,nspec
write(*,*) '**********************************************'
write(*,*) 'Scavenging statistics for species ', species(i), ':'
write(*,*) 'Total number of occurences of below-cloud scavenging', &
& tot_blc_count(i)
write(*,*) 'Total number of occurences of in-cloud scavenging', &
& tot_inc_count(i)
write(*,*) '**********************************************'
end do
do i=1,nspec
if (tot_inc_count(i).gt.0) then
write(*,*) '**********************************************'
write(*,*) 'Scavenging statistics for species ', species(i), ':'
write(*,*) 'Total number of occurences of below-cloud scavenging', &
& tot_blc_count(i)
write(*,*) 'Total number of occurences of in-cloud scavenging', &
& tot_inc_count(i)
write(*,*) '**********************************************'
endif
end do
write (*,*) 'timemanager> call wetdepo'
endif
write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
&XPART MODEL RUN!'
......
......@@ -173,10 +173,8 @@ module com_mod
real :: vset(maxspec,ni),schmi(maxspec,ni),fract(maxspec,ni)
real :: ri(5,numclass),rac(5,numclass),rcl(maxspec,5,numclass)
real :: rgs(maxspec,5,numclass),rlu(maxspec,5,numclass)
real :: rm(maxspec),dryvel(maxspec),kao(maxspec)
real :: rm(maxspec),dryvel(maxspec)
real :: ohcconst(maxspec),ohdconst(maxspec),ohnconst(maxspec)
! se it is possible to associate a species with a second one to make transfer from gas to aerosol
integer :: spec_ass(maxspec)
real :: area_hour(maxspec,24),point_hour(maxspec,24)
real :: area_dow(maxspec,7),point_dow(maxspec,7)
......
......@@ -150,12 +150,9 @@ subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc
do il=2,nz
if (height(il).gt.ztra1(jpart)) then
hz=il-1
! goto 26
exit
endif
end do
!26 continue
if (ngrid.eq.0) then
clouds_v=clouds(ix,jy,hz,n)
......@@ -202,26 +199,19 @@ subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc
endif
!ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp
!ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp - 2 and 3 not used removed by SE
! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms
! for now they are treated the same
grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
grfraction(2)=max(0.05,cc*(lfr(i)))
grfraction(3)=max(0.05,cc*(cfr(j)))
! 2) Computation of precipitation rate in sub-grid cell
!******************************************************
prec(1)=(lsp+convp)/grfraction(1)
prec(2)=(lsp)/grfraction(2)
prec(3)=(convp)/grfraction(3)
! 3) Computation of scavenging coefficients for all species
! Computation of wet deposition
!**********************************************************
if (ngrid.gt.0) then
act_temp=ttn(ix,jy,hz,n,ngrid)
else
......@@ -236,7 +226,6 @@ subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc
! For gas: if positive below-cloud parameters (A or B), and dquer<=0
!******************************************************************
if ((dquer(ks).le.0.).and.(weta_gas(ks).gt.0..or.wetb_gas(ks).gt.0.)) then
! if (weta(ks).gt.0. .or. wetb(ks).gt.0.) then
blc_count(ks)=blc_count(ks)+1
wetscav=weta_gas(ks)*prec(1)**wetb_gas(ks)
......@@ -270,8 +259,6 @@ subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc
endif
! write(*,*) 'bl-cloud, act_temp=',act_temp, ',prec=',prec(1),',wetscav=', wetscav, ', jpart=',jpart
endif ! gas or particle
! endif ! positive below-cloud scavenging parameters given in Species file
endif !end BELOW
......@@ -284,7 +271,6 @@ subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc
! given in species file, or if gas and positive Henry's constant
if ((ccn_aero(ks).gt.0. .or. in_aero(ks).gt.0.).or.(henry(ks).gt.0.and.dquer(ks).le.0)) then
inc_count(ks)=inc_count(ks)+1
! write(*,*) 'Incloud: ',inc_count
! if negative coefficients (turned off) set to zero for use in equation
if (ccn_aero(ks).lt.0.) ccn_aero(ks)=0.
if (in_aero(ks).lt.0.) in_aero(ks)=0.
......@@ -299,7 +285,7 @@ subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc
!ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
! sec test
! cl=1E6*1E-7*prec(1)**0.3 !Sec GFS new
cl=1E6*2E-7*prec(1)**0.36 !Sec ECMWF new
cl=1E6*2E-7*prec(1)**0.36 !Sec ECMWF new, is also suitable for GFS
! cl=2E-7*prec(1)**0.36 !Andreas
! cl=1.6E-6*prec(1)**0.36 !Henrik
endif
......@@ -321,27 +307,16 @@ subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc
!********
if (dquer(ks).gt.0.) then
S_i= frac_act/cl
! write(*,*) 'Si: ',S_i
! GAS
!****
else
cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
!REPLACE to switch old/ new scheme
! S_i=frac_act/cle
S_i=1/cle
endif ! gas or particle
! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
!OLD
if ((readclouds.and.ngrid.eq.0).or.(readclouds_this_nest.and.ngrid.gt.0)) then
!SEC wetscav fix, the cloud height is no longer needed, it gives wrong results
wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)
else
!SEC wetscav fix
wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)
! wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)/clouds_h
endif
endif ! positive in-cloud scavenging parameters given in Species file
endif !incloud
......
......@@ -66,8 +66,7 @@ module netcdf_output_mod
drydep,wetdep,decay,weta_gas,wetb_gas, numbnests, &
ccn_aero,in_aero, & ! wetc_in,wetd_in, &
reldiff,henry,f0,density,dquer,dsigma,dryvel,&
! weightmolar,ohreact,spec_ass,kao,vsetaver,&
weightmolar,ohcconst,ohdconst,spec_ass,kao,vsetaver,&
weightmolar,ohcconst,ohdconst,vsetaver,&
! for concoutput_netcdf and concoutput_nest_netcdf
nxmin1,nymin1,nz,oro,oron,rho,rhon,&
memind,xresoln,yresoln,xrn, xln, yrn,yln,nxn,nyn,&
......@@ -113,6 +112,8 @@ module netcdf_output_mod
logical, parameter :: write_vol = .false.
logical, parameter :: write_area = .false.
! coordinate transformation from internal to world coord
real :: xp1,yp1,xp2,yp2
contains
!****************************************************************
......@@ -511,9 +512,7 @@ subroutine writeheader_netcdf(lnest)
! call nf90_err(nf90_put_att(ncid, sID, 'ohreact', ohreact(i)))
call nf90_err(nf90_put_att(ncid, sID, 'ohcconst', ohcconst(i)))
call nf90_err(nf90_put_att(ncid, sID, 'ohdconst', ohdconst(i)))
call nf90_err(nf90_put_att(ncid, sID, 'kao', kao(i)))
call nf90_err(nf90_put_att(ncid, sID, 'vsetaver', vsetaver(i)))
call nf90_err(nf90_put_att(ncid, sID, 'spec_ass', spec_ass(i)))
if (lnest) then
specIDn(i) = sID
......@@ -534,9 +533,7 @@ subroutine writeheader_netcdf(lnest)
! call nf90_err(nf90_put_att(ncid, sID, 'ohreact', ohreact(i)))
call nf90_err(nf90_put_att(ncid, sID, 'ohcconst', ohcconst(i)))
call nf90_err(nf90_put_att(ncid, sID, 'ohdconst', ohdconst(i)))
call nf90_err(nf90_put_att(ncid, sID, 'kao', kao(i)))
call nf90_err(nf90_put_att(ncid, sID, 'vsetaver', vsetaver(i)))
call nf90_err(nf90_put_att(ncid, sID, 'spec_ass', spec_ass(i)))
if (lnest) then
specIDnppt(i) = sID
......@@ -672,10 +669,14 @@ subroutine writeheader_netcdf(lnest)
call nf90_err(nf90_put_var(ncid, relstartID, ireleasestart(i), (/i/)))
call nf90_err(nf90_put_var(ncid, relendID, ireleaseend(i), (/i/)))
call nf90_err(nf90_put_var(ncid, relkindzID, kindz(i), (/i/)))
call nf90_err(nf90_put_var(ncid, rellng1ID, xpoint1(i), (/i/)))
call nf90_err(nf90_put_var(ncid, rellng2ID, xpoint2(i), (/i/)))
call nf90_err(nf90_put_var(ncid, rellat1ID, ypoint1(i), (/i/)))
call nf90_err(nf90_put_var(ncid, rellat2ID, ypoint2(i), (/i/)))
xp1=xpoint1(i)*dx+xlon0
yp1=ypoint1(i)*dy+ylat0
xp2=xpoint2(i)*dx+xlon0
yp2=ypoint2(i)*dy+ylat0
call nf90_err(nf90_put_var(ncid, rellng1ID, xp1, (/i/)))
call nf90_err(nf90_put_var(ncid, rellng2ID, xp2, (/i/)))
call nf90_err(nf90_put_var(ncid, rellat1ID, yp1, (/i/)))
call nf90_err(nf90_put_var(ncid, rellat2ID, yp2, (/i/)))
call nf90_err(nf90_put_var(ncid, relzz1ID, zpoint1(i), (/i/)))
call nf90_err(nf90_put_var(ncid, relzz2ID, zpoint2(i), (/i/)))
call nf90_err(nf90_put_var(ncid, relpartID, npart(i), (/i/)))
......
......@@ -66,16 +66,18 @@ subroutine readspecies(id_spec,pos_spec)
character(len=16) :: pspecies
real :: pdecay, pweta_gas, pwetb_gas, preldiff, phenry, pf0, pdensity, pdquer
real :: pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, pkao
real :: pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst
real :: pcrain_aero, pcsnow_aero, pccn_aero, pin_aero
integer :: readerror, pspec_ass
real :: parea_dow(7), parea_hour(24), ppoint_dow(7), ppoint_hour(24)
integer :: readerror
! declare namelist
namelist /species_params/ &
pspecies, pdecay, pweta_gas, pwetb_gas, &
pcrain_aero, pcsnow_aero, pccn_aero, pin_aero, &
preldiff, phenry, pf0, pdensity, pdquer, &
pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, pspec_ass, pkao
pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, &
parea_dow, parea_hour, ppoint_dow, ppoint_hour
pspecies="" ! read failure indicator value
pdecay=-999.9
......@@ -95,10 +97,18 @@ subroutine readspecies(id_spec,pos_spec)
pohcconst=-9.99
pohdconst=-9.9E-09
pohnconst=2.0
pspec_ass=-9
pkao=-99.99
pweightmolar=-999.9
do j=1,24 ! initialize everything to no variation
area_hour(pos_spec,j)=1.
point_hour(pos_spec,j)=1.
end do
do j=1,7
area_dow(pos_spec,j)=1.
point_dow(pos_spec,j)=1.
end do
if (readerror.ne.0) then ! text format input
! Open the SPECIES file and read species names and properties
!************************************************************
specnum(pos_spec)=id_spec
......@@ -163,10 +173,18 @@ subroutine readspecies(id_spec,pos_spec)
! write(*,*) ohdconst(pos_spec)
read(unitspecies,'(f8.2)',end=22) ohnconst(pos_spec)
! write(*,*) ohnconst(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)
! Read in daily and day-of-week variation of emissions, if available
!*******************************************************************
read(unitspecies,*,end=22)
do j=1,24 ! 24 hours, starting with 0-1 local time
read(unitspecies,*) ihour,area_hour(pos_spec,j),point_hour(pos_spec,j)
end do
read(unitspecies,*)
do j=1,7 ! 7 days of the week, starting with Monday
read(unitspecies,*) idow,area_dow(pos_spec,j),point_dow(pos_spec,j)
end do
pspecies=species(pos_spec)
pdecay=decay(pos_spec)
......@@ -187,10 +205,18 @@ subroutine readspecies(id_spec,pos_spec)
pohcconst=ohcconst(pos_spec)
pohdconst=ohdconst(pos_spec)
pohnconst=ohnconst(pos_spec)
pspec_ass=spec_ass(pos_spec)
pkao=kao(pos_spec)
else
do j=1,24 ! 24 hours, starting with 0-1 local time
parea_hour(j)=area_hour(pos_spec,j)
ppoint_hour(j)=point_hour(pos_spec,j)
end do
do j=1,7 ! 7 days of the week, starting with Monday
parea_dow(j)=area_dow(pos_spec,j)
ppoint_dow(j)=point_dow(pos_spec,j)
end do
else ! namelist available
species(pos_spec)=pspecies
decay(pos_spec)=pdecay
......@@ -211,8 +237,15 @@ subroutine readspecies(id_spec,pos_spec)
ohcconst(pos_spec)=pohcconst
ohdconst(pos_spec)=pohdconst
ohnconst(pos_spec)=pohnconst
spec_ass(pos_spec)=pspec_ass
kao(pos_spec)=pkao
do j=1,24 ! 24 hours, starting with 0-1 local time
area_hour(pos_spec,j)=parea_hour(j)
point_hour(pos_spec,j)=ppoint_hour(j)
end do
do j=1,7 ! 7 days of the week, starting with Monday
area_dow(pos_spec,j)=parea_dow(j)
point_dow(pos_spec,j)=ppoint_dow(j)
end do
endif
......@@ -302,21 +335,6 @@ subroutine readspecies(id_spec,pos_spec)
end if
end if
if (spec_ass(pos_spec).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
spec_found=.TRUE.
ASSSPEC=.TRUE.
endif
end do
if (spec_found.eqv..false.) then
goto 997
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 ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
......@@ -329,30 +347,6 @@ subroutine readspecies(id_spec,pos_spec)
20 continue
! Read in daily and day-of-week variation of emissions, if available
!*******************************************************************
! HSO: This is not yet implemented as namelist parameters
do j=1,24 ! initialize everything to no variation
area_hour(i,j)=1.
point_hour(i,j)=1.
end do
do j=1,7
area_dow(i,j)=1.
point_dow(i,j)=1.
end do
if (readerror.ne.0) then ! text format input
read(unitspecies,*,end=22)
do j=1,24 ! 24 hours, starting with 0-1 local time
read(unitspecies,*) ihour,area_hour(i,j),point_hour(i,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)
end do
endif
22 close(unitspecies)
......
......@@ -448,7 +448,7 @@ subroutine timemanager(metdata_format)
!CGZ-lifetime: output species lifetime
!write(*,46) float(itime)/3600,itime,numpart
45 format(i13,' SECONDS SIMULATED: ',i13, ' PARTICLES: Uncertainty: ',3f7.3)
45 format(i13,' Seconds simulated: ',i13, ' Particles: Uncertainty: ',3f7.3)
46 format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
if (ipout.ge.1) call partoutput(itime) ! dump particle positions
loutnext=loutnext+loutstep
......
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