Commit 78e62dc0 authored by flexpart's avatar flexpart
Browse files

New OH parameter in SPECIES files (now 3 instead of 2). New path to OH binariy files.

parent 2f8635b3
......@@ -2,7 +2,7 @@
* *
* Input file for the Lagrangian particle dispersion model FLEXPART *
* Definition file of chemical species/radionuclides *
* ESO: edited to work with Rona's changes+Henrik *
* This file is on the format required by version 10 (05/2015) *
****************************************************************************
SO2 Tracer name
-999.9 Species half life
......@@ -13,12 +13,13 @@ SO2 Tracer name
-9.9 Dry deposition (gases) - D
Dry deposition (gases) - Henrys const.
Dry deposition (gases) - f0 (reactivity)
1.9E03 Dry deposition (particles) - rho
1.9E+03 Dry deposition (particles) - rho
4.0E-07 Dry deposition (particles) - dquer
3.0E-01 Dry deposition (particles) - dsig
-9.99 Alternative: dry deposition velocity
64.00 molweight
-1.35E-12 OH Reaction rate - C
1.0 OH Reaction rate - D
1.35E-14 OH Reaction rate - C [cm^3/molecule/sec]
120.0 OH Reaction rate - D [K]
2.0 OH Reaction rate - N (no unit)
-9 number of associated specias (neg. none)
-99.99 KOA - organic matter air partitioning
......@@ -68,7 +68,7 @@ program flexpart
call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))
! FLEXPART version string
flexversion='Version 10.0pre (2015-03-01)'
flexversion='Version 10.0beta (2015-05-01)'
verbosity=0
! Read the pathnames where input/output files are stored
......
......@@ -76,7 +76,7 @@ program flexpart
call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))
! FLEXPART version string
flexversion='Ver. 10.0pre MPI (2015-03-01)'
flexversion='Ver. 10 Beta MPI (2015-05-01)'
verbosity=0
! Read the pathnames where input/output files are stored
......@@ -335,6 +335,9 @@ program flexpart
! and open files that are to be kept open throughout the simulation
!******************************************************************
if (mp_measure_time) call mpif_mtime('iotime',0)
if (lroot) then ! MPI: this part root process only
if (lnetcdfout.eq.1) then
call writeheader_netcdf(lnest=.false.)
else
......@@ -349,7 +352,7 @@ program flexpart
endif
endif
if (lroot) then ! MPI: this part root process only
!
if (verbosity.gt.0) then
print*,'call writeheader'
endif
......@@ -363,9 +366,9 @@ program flexpart
if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
end if ! (mpif_pid == 0)
!open(unitdates,file=path(2)(1:length(2))//'dates')
if (mp_measure_time) call mpif_mtime('iotime',0)
!open(unitdates,file=path(2)(1:length(2))//'dates')
!open(unitdates,file=path(2)(1:length(2))//'dates')
if (verbosity.gt.0 .and. lroot) then
print*,'call openreceptors'
......
......@@ -25,12 +25,14 @@ module com_mod
character :: path(numpath+2*maxnests)*120
integer :: length(numpath+2*maxnests)
character(len=256) :: pathfile, flexversion, arg1, arg2
character(len=256) :: ohfields_path
! path path names needed for trajectory model
! length length of path names needed for trajectory model
! pathfile file where pathnames are stored
! flexversion version of flexpart
! arg input arguments from launch at command line
! ohfields_path path to binary files for OH fields
!********************************************************
! Variables defining the general model run specifications
......@@ -162,7 +164,7 @@ module com_mod
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 :: ohcconst(maxspec),ohdconst(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)
......
......@@ -42,10 +42,13 @@ subroutine getfields(itime,nstop,memstat)
! Function of nstop extended.
!
! eso 2014:
! MPI version. 3 fields instead of 2, to allow reading the newest in
! the background.
! Only one process (lmpreader=.true.) does the actual reading, while the
! rest call this routine only to update memind, memstat etc.
! MPI version.
! If running with number of processes >= mpi_mod::read_grp_min,
! only one process (mpi_mod::lmpreader=.true.) does the actual reading, while
! the rest call this routine only to update memind, memstat etc.
!
! If mpi_mod::lmp_sync=.true., uses 3 fields instead of 2, to allow reading
! the newest in the background.
!
! Return memstat, which is the sum of
!
......@@ -60,8 +63,8 @@ subroutine getfields(itime,nstop,memstat)
! lwindinterval [s] time difference between the two wind fields read in *
! indj indicates the number of the wind field to be read in *
! indmin remembers the number of wind fields already treated *
! memind(2) pointer, on which place the wind fields are stored *
! memtime(2) [s] times of the wind fields, which are kept in memory *
! memind(2[3]) pointer, on which place the wind fields are stored *
! memtime(2[3]) [s] times of the wind fields, which are kept in memory *
! itime [s] current time since start date of trajectory calcu- *
! lation *
! nstop > 0, if trajectory has to be terminated *
......
......@@ -119,7 +119,7 @@ subroutine gethourlyOH(itime)
m1=(jjjjmmdd-(jjjjmmdd/10000)*10000)/100
memOHtime(1)=0.
jul2=bdate+real(1./24.,kind=dp) ! date for next hour
jul2=bdate+ldirect*real(1./24.,kind=dp) ! date for next hour
call caldate(jul2,jjjjmmdd,hhmmss)
m2=(jjjjmmdd-(jjjjmmdd/10000)*10000)/100
memOHtime(2)=ldirect*3600.
......
......@@ -89,7 +89,7 @@ module mpi_mod
! MPI tags/requests for send/receive operation
integer :: tm1
integer, parameter :: nvar_async=27 !29 :DBG:
!integer, dimension(:), allocatable :: tags
!integer, dimension(:), allocatable :: tags
integer, dimension(:), allocatable :: reqs
......@@ -116,7 +116,7 @@ module mpi_mod
logical, parameter :: mp_dev_mode = .false.
logical, parameter :: mp_dbg_out = .false.
logical, parameter :: mp_time_barrier=.true.
logical, parameter :: mp_measure_time=.true.
logical, parameter :: mp_measure_time=.false.
! for measuring CPU/Wall time
real(sp) :: mp_comm_time_beg, mp_comm_time_end, mp_comm_time_total=0.
......@@ -138,6 +138,8 @@ module mpi_mod
real(dp) :: mp_advance_wtime_beg, mp_advance_wtime_end, mp_advance_wtime_total=0.
real(dp) :: mp_conccalc_time_beg, mp_conccalc_time_end, mp_conccalc_time_total=0.
real(dp) :: mp_total_wtime_beg, mp_total_wtime_end, mp_total_wtime_total=0.
real(dp) :: mp_vt_wtime_beg, mp_vt_wtime_end, mp_vt_wtime_total
real(sp) :: mp_vt_time_beg, mp_vt_time_end, mp_vt_time_total
! dat_lun logical unit number for i/o
integer, private :: dat_lun
......@@ -417,7 +419,6 @@ contains
end if
! redefine numpart as 'numpart per process' throughout the code
!**************************************************************
numpart = numpart_mpi
......@@ -438,10 +439,6 @@ contains
integer :: i
!***********************************************************************
! Time for MPI communications
!****************************
if (mp_measure_time) call mpif_mtime('commtime',0)
......@@ -1242,7 +1239,7 @@ contains
!*****************************************************
do dest=0,mp_np-1 ! mp_np-2 will also work if last proc reserved for reading
! TODO: use mp_partgroup_np here
! TODO: use mp_partgroup_np here
if (dest.eq.id_read) cycle
i=dest*nvar_async
call MPI_Isend(uu(:,:,:,mind),d3s1,mp_pp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
......@@ -1338,8 +1335,8 @@ contains
call MPI_Isend(ciwc(:,:,:,mind),d3s1,mp_pp,dest,tm1,&
&MPI_COMM_WORLD,reqs(i),mp_ierr)
if (mp_ierr /= 0) goto 600
! else
! i=i+2
! else
! i=i+2
end if
end do
......@@ -1388,10 +1385,10 @@ contains
!*******************************************************************************
! :TODO: don't need these
! d3s1=d3_size1
! d3s2=d3_size2
! d2s1=d2_size1
! d2s2=d2_size2
! d3s1=d3_size1
! d3s2=d3_size2
! d2s1=d2_size1
! d2s2=d2_size2
! At the time this immediate receive is posted, memstat is the state of
! windfield indices at the previous/current time. From this, the future
......@@ -1592,15 +1589,15 @@ contains
! if (readclouds) then
call MPI_Waitall(n_req,reqs,MPI_STATUSES_IGNORE,mp_ierr)
! endif
! else
! do i = 0, nvar_async*mp_np-1
! if (mod(i,27).eq.0 .or. mod(i,28).eq.0) then
! call MPI_Cancel(reqs(i),mp_ierr)
! cycle
! end if
! call MPI_Wait(reqs(i),MPI_STATUS_IGNORE,mp_ierr)
! end do
! end if
! else
! do i = 0, nvar_async*mp_np-1
! if (mod(i,27).eq.0 .or. mod(i,28).eq.0) then
! call MPI_Cancel(reqs(i),mp_ierr)
! cycle
! end if
! call MPI_Wait(reqs(i),MPI_STATUS_IGNORE,mp_ierr)
! end do
! end if
if (mp_ierr /= 0) goto 600
......@@ -1733,7 +1730,6 @@ contains
& mp_comm_used, mp_ierr)
end if
if ((WETDEP).and.(ldirect.gt.0)) then
call MPI_Reduce(wetgriduncn, wetgriduncn0, grid_size2d, mp_pp, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
......@@ -1853,6 +1849,20 @@ contains
& mp_io_time_beg)
end if
case ('verttransform')
if (imode.eq.0) then
mp_vt_wtime_beg = mpi_wtime()
call cpu_time(mp_vt_time_beg)
else
mp_vt_wtime_end = mpi_wtime()
call cpu_time(mp_vt_time_end)
mp_vt_wtime_total = mp_vt_wtime_total + (mp_vt_wtime_end - &
& mp_vt_wtime_beg)
mp_vt_time_total = mp_vt_time_total + (mp_vt_time_end - &
& mp_vt_time_beg)
end if
case ('readwind')
if (imode.eq.0) then
call cpu_time(mp_readwind_time_beg)
......@@ -1908,7 +1918,7 @@ contains
!***********************************************************************
if (mp_measure_time) then
IF (mp_measure_time) THEN
do ip=0, mp_np-1
call MPI_BARRIER(MPI_COMM_WORLD, mp_ierr)
......@@ -1955,23 +1965,38 @@ contains
& mp_wetdepo_time_total
write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR CONCCALC:',&
& mp_conccalc_time_total
! write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR VERTTRANSFORM:',&
! & mp_vt_wtime_total
! write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR VERTTRANSFORM:',&
! & mp_vt_time_total
! NB: the 'flush' function is possibly a gfortran-specific extension
call flush()
end if
end do
end if
! This call to barrier is for correctly formatting output
call MPI_BARRIER(MPI_COMM_WORLD, mp_ierr)
if (lroot) then
write(*,FMT='(72("#"))')
WRITE(*,*) "To turn off output of time measurements, set "
WRITE(*,*) " mp_measure_time=.false."
WRITE(*,*) "in file mpi_mod.f90"
write(*,FMT='(72("#"))')
end if
! j=mp_pid*nvar_async
! In the implementation with 3 fields, the processes may have posted
! MPI_Irecv requests that should be cancelled here
!! TODO:
! if (.not.lmp_sync) then
! r=mp_pid*nvar_async
! do j=r,r+nvar_async-1
! call MPI_Cancel(j,mp_ierr)
! if (mp_ierr /= 0) write(*,*) '#### mpif_finalize::MPI_Cancel> ERROR ####'
! end do
! end if
! if (.not.lmp_sync) then
! r=mp_pid*nvar_async
! do j=r,r+nvar_async-1
! call MPI_Cancel(j,mp_ierr)
! if (mp_ierr /= 0) write(*,*) '#### mpif_finalize::MPI_Cancel> ERROR ####'
! end do
! end if
call MPI_FINALIZE(mp_ierr)
if (mp_ierr /= 0) then
......@@ -1980,61 +2005,61 @@ contains
end if
end subroutine mpif_finalize
end subroutine mpif_finalize
subroutine get_lun(my_lun)
subroutine get_lun(my_lun)
!***********************************************************************
! get_lun:
! Starting from 100, get next free logical unit number
!***********************************************************************
implicit none
implicit none
integer, intent(inout) :: my_lun
integer, save :: free_lun=100
logical :: exists, iopen
integer, intent(inout) :: my_lun
integer, save :: free_lun=100
logical :: exists, iopen
!***********************************************************************
loop1: do
inquire(UNIT=free_lun, EXIST=exists, OPENED=iopen)
if (exists .and. .not.iopen) exit loop1
free_lun = free_lun+1
end do loop1
my_lun = free_lun
loop1: do
inquire(UNIT=free_lun, EXIST=exists, OPENED=iopen)
if (exists .and. .not.iopen) exit loop1
free_lun = free_lun+1
end do loop1
my_lun = free_lun
end subroutine get_lun
end subroutine get_lun
subroutine write_data_dbg(array_in, array_name, tstep, ident)
subroutine write_data_dbg(array_in, array_name, tstep, ident)
!***********************************************************************
! Write one-dimensional arrays to disk (for debugging purposes)
! Write one-dimensional arrays to file (for debugging purposes)
!***********************************************************************
implicit none
implicit none
real, intent(in), dimension(:) :: array_in
integer, intent(in) :: tstep
integer :: lios
character(LEN=*), intent(in) :: ident, array_name
real, intent(in), dimension(:) :: array_in
integer, intent(in) :: tstep
integer :: lios
character(LEN=*), intent(in) :: ident, array_name
character(LEN=8) :: c_ts
character(LEN=40) :: fn_1, fn_2
character(LEN=8) :: c_ts
character(LEN=40) :: fn_1, fn_2
!***********************************************************************
write(c_ts, FMT='(I8.8,BZ)') tstep
fn_1='-'//trim(adjustl(c_ts))//'-'//trim(ident)
write(c_ts, FMT='(I2.2,BZ)') mp_np
fn_2= trim(adjustl(array_name))//trim(adjustl(fn_1))//'-np'//trim(adjustl(c_ts))//'.dat'
write(c_ts, FMT='(I8.8,BZ)') tstep
fn_1='-'//trim(adjustl(c_ts))//'-'//trim(ident)
write(c_ts, FMT='(I2.2,BZ)') mp_np
fn_2= trim(adjustl(array_name))//trim(adjustl(fn_1))//'-np'//trim(adjustl(c_ts))//'.dat'
call get_lun(dat_lun)
open(UNIT=dat_lun, FILE=fn_2, IOSTAT=lios, ACTION='WRITE', &
FORM='UNFORMATTED', STATUS='REPLACE')
write(UNIT=dat_lun, IOSTAT=lios) array_in
close(UNIT=dat_lun)
call get_lun(dat_lun)
open(UNIT=dat_lun, FILE=fn_2, IOSTAT=lios, ACTION='WRITE', &
FORM='UNFORMATTED', STATUS='REPLACE')
write(UNIT=dat_lun, IOSTAT=lios) array_in
close(UNIT=dat_lun)
end subroutine write_data_dbg
end subroutine write_data_dbg
end module mpi_mod
......@@ -152,7 +152,7 @@ subroutine ohreaction(itime,ltsample,loutnext)
do k=1,nspec
if (ohcconst(k).gt.0.) then
ohrate=ohcconst(k)*temp**2*exp(-ohdconst(k)/temp)*oh_average
ohrate=ohcconst(k)*temp**ohnconst(k)*exp(-ohdconst(k)/temp)*oh_average
! new particle mass
restmass = xmass1(jpart,k)*exp(-1*ohrate*abs(ltsample))
if (restmass .gt. smallnum) then
......
......@@ -45,7 +45,6 @@ subroutine readOHfield
use par_mod
use com_mod
implicit none
include 'netcdf.inc'
......@@ -68,7 +67,8 @@ subroutine readOHfield
! open netcdf file
write(mm,fmt='(i2.2)') m
thefile=trim(path(1))//'OH_FIELDS/'//'geos-chem.OH.2005'//mm//'01.nc'
! thefile=trim(path(1))//'OH_FIELDS/'//'geos-chem.OH.2005'//mm//'01.nc'
thefile=trim(ohfields_path)//'OH_FIELDS/'//'geos-chem.OH.2005'//mm//'01.nc'
ierr=nf_open(trim(thefile),NF_NOWRITE,nid)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
......@@ -162,7 +162,8 @@ subroutine readOHfield
!********************************
! open netcdf file
thefile=trim(path(1))//'OH_FIELDS/jrate_average.nc'
! thefile=trim(path(1))//'OH_FIELDS/jrate_average.nc'
thefile=trim(ohfields_path)//'OH_FIELDS/jrate_average.nc'
ierr=nf_open(trim(thefile),NF_NOWRITE,nid)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
......
......@@ -110,7 +110,8 @@ subroutine readcommand
linit_cond, &
lnetcdfout, &
surf_only, &
cblflag
cblflag, &
ohfields_path
! Presetting namelist command
ldirect=0
......@@ -142,6 +143,7 @@ subroutine readcommand
lnetcdfout=0
surf_only=0
cblflag=0
ohfields_path="../../flexin/"
! Open the command file and read user options
! Namelist input first: try to read as namelist file
......@@ -231,6 +233,7 @@ subroutine readcommand
read(unitcommand,*) surf_only
if (old) call skplin(3,unitcommand) !added by mc
read(unitcommand,*) cblflag !added by mc
close(unitcommand)
endif ! input format
......
......@@ -384,7 +384,7 @@ subroutine readreleases
if (ohcconst(i).gt.0.) then
OHREA=.true.
write (*,*) 'OHreaction switched on: ',ohcconst(i),i
if (lroot) write (*,*) 'OHreaction switched on: ',ohcconst(i),i
endif
if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or.(dryvel(i).gt.0.)) then
......
......@@ -46,6 +46,7 @@ subroutine readspecies(id_spec,pos_spec)
! wetb_in Parameter for determining in-cloud scavenging *
! ohcconst OH reaction rate constant C *
! ohdconst OH reaction rate constant D *
! ohnconst OH reaction rate constant n *
! id_spec SPECIES number as referenced in RELEASE file *
! id_pos position where SPECIES data shall be stored *
! *
......@@ -65,7 +66,7 @@ subroutine readspecies(id_spec,pos_spec)
character(len=16) :: pspecies
real :: pdecay, pweta, pwetb, preldiff, phenry, pf0, pdensity, pdquer
real :: pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pspec_ass, pkao
real :: pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, pspec_ass, pkao
real :: pweta_in, pwetb_in, pwetc_in, pwetd_in
integer :: readerror
......@@ -74,7 +75,7 @@ subroutine readspecies(id_spec,pos_spec)
pspecies, pdecay, pweta, pwetb, &
pweta_in, pwetb_in, pwetc_in, pwetd_in, &
preldiff, phenry, pf0, pdensity, pdquer, &
pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pspec_ass, pkao
pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, pspec_ass, pkao
pspecies=" "
pdecay=-999.9
......@@ -93,6 +94,7 @@ subroutine readspecies(id_spec,pos_spec)
pdryvel=-9.99
pohcconst=-9.99
pohdconst=-9.9E-09
pohnconst=2.0
pspec_ass=-9
pkao=-99.99
pweightmolar=-789.0 ! read failure indicator value
......@@ -159,6 +161,8 @@ subroutine readspecies(id_spec,pos_spec)
! write(*,*) ohcconst(pos_spec)
read(unitspecies,'(f8.2)',end=22) ohdconst(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)
......@@ -182,6 +186,7 @@ subroutine readspecies(id_spec,pos_spec)
pweightmolar=weightmolar(pos_spec)
pohcconst=ohcconst(pos_spec)
pohdconst=ohdconst(pos_spec)
pohnconst=ohnconst(pos_spec)
pspec_ass=spec_ass(pos_spec)
pkao=kao(pos_spec)
......@@ -205,6 +210,7 @@ subroutine readspecies(id_spec,pos_spec)
weightmolar(pos_spec)=pweightmolar
ohcconst(pos_spec)=pohcconst
ohdconst(pos_spec)=pohdconst
ohnconst(pos_spec)=pohnconst
spec_ass(pos_spec)=pspec_ass
kao(pos_spec)=pkao
......
......@@ -656,7 +656,7 @@ subroutine timemanager
if (linit_cond.ge.1) call initial_cond_output(itime) ! dump initial cond. field
close(104)
!close(104)
! De-allocate memory and end
!***************************
......
......@@ -207,16 +207,22 @@ subroutine timemanager
write (*,*) 'timemanager> call getfields'
endif
! This time measure includes reading/MPI communication (for the reader process),
! or MPI communication time only (for other processes)
if (mp_measure_time) call mpif_mtime('getfields',0)
call getfields(itime,nstop1,memstat)
if (mp_measure_time) call mpif_mtime('getfields',1)
! Broadcast fields to all MPI processes
! Skip if all processes have called getfields or if no new fields
!*****************************************************************
if (mp_measure_time.and..not.(lmpreader.and.lmp_use_reader)) call mpif_mtime('getfields',0)
! Version 1 (lmp_sync=.true.) uses a read-ahead process where send/recv is done
! in sync at start of each new field time interval
if (lmp_sync.and.lmp_use_reader.and.memstat.gt.0) then
......@@ -258,6 +264,9 @@ subroutine timemanager
end if
if (mp_measure_time.and..not.(lmpreader.and.lmp_use_reader)) call mpif_mtime('getfields',1)
!*******************************************************************************
if (lmpreader.and.nstop1.gt.1) stop 'NO METEO FIELDS AVAILABLE'
......@@ -445,6 +454,7 @@ subroutine timemanager
!**************************************
call mpif_tm_reduce_grid
if (mp_measure_time) call mpif_mtime('iotime',0)
if (surf_only.ne.1) then
if (lroot) then
if (lnetcdfout.eq.1) then
......@@ -469,6 +479,7 @@ subroutine timemanager
gridunc(:,:,:,:,:,:,:)=0.
endif
endif
if (mp_measure_time) call mpif_mtime('iotime',1)
! :TODO: Correct calling of conc_surf above?
......@@ -480,6 +491,8 @@ subroutine timemanager
! MPI: Root process collects/sums nested grids
!*********************************************
call mpif_tm_reduce_grid_nest
if (mp_measure_time) call mpif_mtime('iotime',0)
if (lnetcdfout.eq.0) then
if (surf_only.ne.1) then
......@@ -514,11 +527,14 @@ subroutine timemanager