Skip to content

GitLab

Commit 505a96e3 authored by Ignacio Pisso's avatar Ignacio Pisso
Browse files

Pull: Merge branch 'dev' of https://git.nilu.no/flexpart/flexpart into dev

parents 3d048456 20963b15
Hide whitespace changes
Inline Side-by-side
src/FLEXPART.f90
View file @ 505a96e3
......@@ -51,10 +51,14 @@ program flexpart
use par_mod
use com_mod
use conv_mod
use netcdf_output_mod, only: writeheader_netcdf
use random_mod, only: gasdev1
use class_gribfile
#ifdef USE_NCF
use netcdf_output_mod, only: writeheader_netcdf
#endif
implicit none
integer :: i,j,ix,jy,inest
......@@ -351,6 +355,7 @@ program flexpart
! and open files that are to be kept open throughout the simulation
!******************************************************************
#ifdef USE_NCF
if (lnetcdfout.eq.1) then
call writeheader_netcdf(lnest=.false.)
else
......@@ -364,6 +369,7 @@ program flexpart
call writeheader_nest
endif
endif
#endif
if (verbosity.gt.0) then
print*,'call writeheader'
......
src/FLEXPART_MPI.f90
View file @ 505a96e3
......@@ -52,10 +52,13 @@ program flexpart
use com_mod
use conv_mod
use mpi_mod
use netcdf_output_mod, only: writeheader_netcdf
use random_mod, only: gasdev1
use class_gribfile
#ifdef USE_NCF
use netcdf_output_mod, only: writeheader_netcdf
#endif
implicit none
integer :: i,j,ix,jy,inest
......@@ -63,7 +66,8 @@ program flexpart
character(len=256) :: inline_options !pathfile, flexversion, arg2
integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
integer :: detectformat
integer(selected_int_kind(16)), dimension(maxspec) :: tot_b=0, &
& tot_i=0
! Initialize mpi
......@@ -202,11 +206,11 @@ program flexpart
metdata_format = detectformat()
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
print *,'ECMWF metdata detected'
if (lroot) print *,'ECMWF metdata detected'
elseif (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
print *,'NCEP metdata detected'
if (lroot) print *,'NCEP metdata detected'
else
print *,'Unknown metdata format'
if (lroot) print *,'Unknown metdata format'
stop
endif
......@@ -377,23 +381,24 @@ program flexpart
!******************************************************************
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
call writeheader
end if
if (nested_output.eq.1) then
if (lnetcdfout.eq.1) then
call writeheader_netcdf(lnest=.true.)
else
call writeheader_nest
if (lroot) then ! MPI: this part root process only
#ifdef USE_NCF
if (lnetcdfout.eq.1) then
call writeheader_netcdf(lnest=.false.)
else
call writeheader
end if
if (nested_output.eq.1) then
if (lnetcdfout.eq.1) then
call writeheader_netcdf(lnest=.true.)
else
call writeheader_nest
endif
endif
endif
#endif
!
if (verbosity.gt.0) then
print*,'call writeheader'
endif
......@@ -401,7 +406,7 @@ program flexpart
call writeheader
! FLEXPART 9.2 ticket ?? write header in ASCII format
call writeheader_txt
!if (nested_output.eq.1) call writeheader_nest
if (nested_output.eq.1.and.surf_only.ne.1) call writeheader_nest
if (nested_output.eq.1.and.surf_only.eq.1) call writeheader_nest_surf
if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
......@@ -409,8 +414,6 @@ program flexpart
if (mp_measure_time) call mpif_mtime('iotime',0)
!open(unitdates,file=path(2)(1:length(2))//'dates')
if (verbosity.gt.0 .and. lroot) then
print*,'call openreceptors'
endif
......@@ -480,28 +483,24 @@ program flexpart
! NIK 16.02.2005
if (lroot) then
call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
if (mp_partgroup_pid.ge.0) then ! Skip for readwind process
call MPI_Reduce(tot_blc_count, tot_b, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
call MPI_Reduce(tot_inc_count, tot_i, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
else
if (mp_partgroup_pid.ge.0) then ! Skip for readwind process
call MPI_Reduce(tot_blc_count, 0, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
call MPI_Reduce(tot_inc_count, 0, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
end if
end if
if (lroot) then
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)
& tot_b(i)
! & tot_blc_count(i)
write(*,*) 'Total number of occurences of in-cloud scavenging', &
& tot_inc_count(i)
& tot_i(i)
! & tot_inc_count(i)
write(*,*) '**********************************************'
end do
......
src/concoutput.f90
View file @ 505a96e3
......@@ -139,6 +139,8 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
write(unitdates,'(a)') adate//atime
close(unitdates)
! Overwrite existing dates file on first call, later append to it
! This fixes a bug where the dates file kept growing across multiple runs
IF (init) THEN
file_stat='OLD'
init=.false.
......@@ -312,7 +314,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
drygridsigma(ix,jy)= &
drygridsigma(ix,jy)* &
sqrt(real(nclassunc))
125 drygridsigmatotal=drygridsigmatotal+ &
drygridsigmatotal=drygridsigmatotal+ &
drygridsigma(ix,jy)
endif
......
src/concoutput_mpi.f90
View file @ 505a96e3
......@@ -103,7 +103,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
real,parameter :: weightair=28.97
logical :: sp_zer
LOGICAL,save :: init=.true.
logical,save :: init=.true.
character :: adate*8,atime*6
character(len=3) :: anspec
integer :: mind
......@@ -112,7 +112,6 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
logical :: ldates_file
integer :: ierr
character(LEN=100) :: dates_char
! character :: dates_char
! Measure execution time
if (mp_measure_time) call mpif_mtime('rootonly',0)
......@@ -128,7 +127,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! Overwrite existing dates file on first call, later append to it
! This fixes a bug where the dates file kept growing across multiple runs
! If 'dates' file exists, make a backup
! If 'dates' file exists in output directory, make a backup
inquire(file=path(2)(1:length(2))//'dates', exist=ldates_file)
if (ldates_file.and.init) then
open(unit=unitdates, file=path(2)(1:length(2))//'dates',form='formatted', &
......@@ -257,23 +256,34 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
do ks=1,nspec
write(anspec,'(i3.3)') ks
if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
if (ldirect.eq.1) then
open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// &
atime//'_'//anspec,form='unformatted')
else
open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
atime//'_'//anspec,form='unformatted')
endif
write(unitoutgrid) itime
endif
if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
if (DRYBKDEP.or.WETBKDEP) then !scavdep output
if (DRYBKDEP) &
open(unitoutgrid,file=path(2)(1:length(2))//'grid_drydep_'//adate// &
atime//'_'//anspec,form='unformatted')
if (WETBKDEP) &
open(unitoutgrid,file=path(2)(1:length(2))//'grid_wetdep_'//adate// &
atime//'_'//anspec,form='unformatted')
write(unitoutgrid) itime
else
if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
if (ldirect.eq.1) then
open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// &
atime//'_'//anspec,form='unformatted')
else
open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
atime//'_'//anspec,form='unformatted')
endif
write(unitoutgrid) itime
endif
write(unitoutgridppt) itime
endif
if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
atime//'_'//anspec,form='unformatted')
write(unitoutgridppt) itime
endif
endif ! if deposition output
do kp=1,maxpointspec_act
do nage=1,nageclass
......@@ -353,7 +363,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! Concentration output
!*********************
if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5).or.(iout.eq.6)) then
! Wet deposition
sp_count_i=0
......@@ -448,10 +458,18 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
sp_fact=sp_fact*(-1.)
endif
sp_count_r=sp_count_r+1
if (lparticlecountoutput) then
sparse_dump_r(sp_count_r)= &
sp_fact* &
grid(ix,jy,kz)
else
sparse_dump_r(sp_count_r)= &
sp_fact* &
grid(ix,jy,kz)* &
factor3d(ix,jy,kz)/tot_mu(ks,kp)
end if
! if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
! + write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
! sparse_dump_u(sp_count_r)=
......@@ -637,24 +655,8 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! Reinitialization of grid
!*************************
do ks=1,nspec
do kp=1,maxpointspec_act
do i=1,numreceptor
creceptor(i,ks)=0.
end do
do jy=0,numygrid-1
do ix=0,numxgrid-1
do l=1,nclassunc
do nage=1,nageclass
do kz=1,numzgrid
gridunc(ix,jy,kz,ks,kp,l,nage)=0.
end do
end do
end do
end do
end do
end do
end do
creceptor(:,:)=0.
gridunc(:,:,:,:,:,:,:)=0.
if (mp_measure_time) call mpif_mtime('rootonly',1)
......
src/getfields_mpi.f90
View file @ 505a96e3
......@@ -232,6 +232,10 @@ subroutine getfields(itime,nstop,memstat,metdata_format)
end do
40 indmin=indj
if (WETBKDEP.and.(lmpreader.or.(.not.lmp_use_reader.and.lroot))) then
call writeprecip(itime,memind(1))
endif
else
! No wind fields, which can be used, are currently in memory
......@@ -286,7 +290,10 @@ subroutine getfields(itime,nstop,memstat,metdata_format)
end do
60 indmin=indj
mind3=memstat
! if (WETBKDEP.and.lroot) then
if (WETBKDEP.and.(lmpreader.or.(.not.lmp_use_reader.and.lroot))) then
call writeprecip(itime,memind(1))
endif
endif
......
src/makefile
View file @ 505a96e3
......@@ -30,6 +30,10 @@ SHELL = /bin/bash
# Compile for debugging parallel FLEXPART
# make [-j] mpi-dbg
#
# NETCDF OUTPUT
# To add support for output in netCDF format, append `ncf=yes` to the
# `make` command
#
################################################################################
## PROGRAMS
......@@ -56,8 +60,7 @@ ifeq ($(gcc), 4.9)
INCPATH1 = ${ROOT_DIR}/gcc-4.9.1/include
INCPATH2 = ${ROOT_DIR}/include
LIBPATH1 = ${ROOT_DIR}/lib
else #ifeq ($(gcc), 5.4)
else
# Compiled libraries under user ~flexpart, gfortran v5.4
ROOT_DIR = /homevip/flexpart/
......@@ -67,18 +70,18 @@ else #ifeq ($(gcc), 5.4)
INCPATH1 = ${ROOT_DIR}/gcc-5.4.0/include
INCPATH2 = /usr/include
LIBPATH1 = ${ROOT_DIR}/gcc-5.4.0/lib
endif
#else
# Default: System libraries at NILU, gfortran v4.6
# F90 = /usr/bin/gfortran
# MPIF90 = /usr/bin/mpif90.openmpi
# INCPATH1 = /xnilu_wrk/projects/FLEXPART/flex_wrk/bin64/grib_api/include
# INCPATH2 = /usr/include
# LIBPATH1 = /xnilu_wrk/projects/FLEXPART/flex_wrk/bin64/grib_api/lib
### Enable netCDF output?
ifeq ($(ncf), yes)
NCOPT = -DUSE_NCF -lnetcdff
else
NCOPT = -UUSE_NCF
endif
# path to gributils used to detect meteodata format
VPATH = gributils/
......@@ -88,11 +91,12 @@ O_LEV = 2 # [0,1,2,3,g,s,fast]
O_LEV_DBG = g # [0,g]
## LIBRARIES
LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper -lnetcdff # -fopenmp
#LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper -lnetcdff
LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper $(NCOPT)
FFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER) #-Warray-bounds -fcheck=all # -march=native
FFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(NCOPT) $(FUSER) #-Warray-bounds -fcheck=all # -march=native
DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace -Wall -fdump-core $(FUSER) # -ffpe-trap=invalid,overflow,denormal,underflow,zero -Warray-bounds -fcheck=all
DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) $(NCOPT) -fbacktrace -Wall -fdump-core $(FUSER) # -ffpe-trap=invalid,overflow,denormal,underflow,zero -Warray-bounds -fcheck=all
LDFLAGS = $(FFLAGS) -L$(LIBPATH1) -Wl,-rpath,$(LIBPATH1) $(LIBS) #-L$(LIBPATH2)
LDDEBUG = $(DBGFLAGS) -L$(LIBPATH1) $(LIBS) #-L$(LIBPATH2)
......@@ -139,6 +143,8 @@ OBJECTS_MPI = releaseparticles_mpi.o partoutput_mpi.o \
getfields_mpi.o \
readwind_ecmwf_mpi.o
OBJECTS_NCF = netcdf_output_mod.o
OBJECTS = \
advance.o initialize.o \
writeheader.o writeheader_txt.o \
......@@ -195,7 +201,11 @@ ohreaction.o getvdep_nests.o \
initial_cond_calc.o initial_cond_output.o \
dynamic_viscosity.o get_settling.o \
initialize_cbl_vel.o re_initialize_particle.o \
cbl.o netcdf_output_mod.o
cbl.o
ifeq ($(ncf), yes)
OBJECTS := $(OBJECTS) $(OBJECTS_NCF)
endif
%.o: %.mod
......
src/readOHfield.f90
View file @ 505a96e3
......@@ -21,25 +21,27 @@
subroutine readOHfield
!*****************************************************************************
! *
! Reads the OH field into memory *
! *
! AUTHOR: R.L. Thompson, Nov 2014 *
! *
!*****************************************************************************
! *
! Variables: *
! *
! path(numpath) contains the path names *
! lonOH(nxOH) longitude of OH fields *
! latOH(nyOH) latitude of OH fields *
! altOH(nzOH) altitude of OH fields *
! etaOH(nzOH) eta-levels of OH fields *
! OH_field(nxOH,nyOH,nzOH,m) OH concentration (molecules/cm3) *
! *
! *
!*****************************************************************************
!*****************************************************************************
! *
! Reads the OH field into memory *
! *
! AUTHOR: R.L. Thompson, Nov 2014 *
! *
! UPDATES: *
! 03/2018 SEC: Converted original netCDF files to binary format *
!*****************************************************************************
! *
! Variables: *
! *
! path(numpath) contains the path names *
! lonOH(nxOH) longitude of OH fields *
! latOH(nyOH) latitude of OH fields *
! altOH(nzOH) altitude of OH fields *
! etaOH(nzOH) eta-levels of OH fields *
! OH_field(nxOH,nyOH,nzOH,m) OH concentration (molecules/cm3) *
! *
! *
!*****************************************************************************
use oh_mod
use par_mod
......@@ -47,11 +49,7 @@ subroutine readOHfield
implicit none
include 'netcdf.inc'
character(len=150) :: thefile
character(len=2) :: mm
integer :: nid,ierr,xid,yid,zid,vid,m
integer :: i,j,k,l,ierr
real, dimension(:), allocatable :: etaOH
! real, parameter :: gasct=8.314 ! gas constant
......@@ -60,144 +58,39 @@ subroutine readOHfield
! real, parameter :: lrate=0.0065 ! K m-1
real, parameter :: scalehgt=7000. ! scale height in metres
! Read OH fields and level heights
!********************************
do m=1,12
! open netcdf file
write(mm,fmt='(i2.2)') m
! 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 (*,*) 'The OH field at: '//thefile// ' could not be opened'
write (*,*) 'please copy the OH fields there, or change the path in the'
write (*,*) 'COMMAND namelist!'
write(*,*) nf_strerror(ierr)
stop
endif
! inquire about variables
ierr=nf_inq_dimid(nid,'Lon-000',xid)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_dimid(nid,'Lat-000',yid)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_dimid(nid,'Alt-000',zid)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
if(m.eq.1) then
! read dimension sizes
ierr=nf_inq_dimlen(nid,xid,nxOH)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_dimlen(nid,yid,nyOH)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_dimlen(nid,zid,nzOH)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
! allocate variables
allocate(lonOH(nxOH))
allocate(latOH(nyOH))
allocate(etaOH(nzOH))
allocate(altOH(nzOH))
allocate(OH_field(nxOH,nyOH,nzOH,12))
allocate(OH_hourly(nxOH,nyOH,nzOH,2))
! read dimension variables
ierr=nf_inq_varid(nid,'LON',xid)
ierr=nf_get_var_real(nid,xid,lonOH)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_varid(nid,'LAT',yid)
ierr=nf_get_var_real(nid,yid,latOH)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_varid(nid,'ETAC',zid)
ierr=nf_get_var_real(nid,zid,etaOH)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
! convert eta-level to altitude (assume surface pressure of 1010 hPa)
altOH=log(1010./(etaOH*1010.))*scalehgt
endif ! m.eq.1
! read OH_field
ierr=nf_inq_varid(nid,'CHEM-L_S__OH',vid)
ierr=nf_get_var_real(nid,vid,OH_field(:,:,:,m))
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_close(nid)
end do
deallocate(etaOH)
! Read J(O1D) photolysis rates
!********************************
open(unitOH,file=trim(ohfields_path) &
//'OH_FIELDS/OH_variables.bin',status='old', &
form='UNFORMATTED', iostat=ierr, convert='little_endian')
! open netcdf file
! 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)
write(*,*) 'Cannot read binary OH fields in ',trim(ohfields_path)//'OH_FIELDS/OH_variables.bin'
stop
endif
! read dimension variables
ierr=nf_inq_varid(nid,'longitude',xid)
ierr=nf_get_var_real(nid,xid,lonjr)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_varid(nid,'latitude',yid)
ierr=nf_get_var_real(nid,yid,latjr)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
! read jrate_average
ierr=nf_inq_varid(nid,'jrate',vid)
ierr=nf_get_var_real(nid,vid,jrate_average)
if(ierr.ne.0) then