diff --git a/src/FLEXPART.f90 b/src/FLEXPART.f90
index cd3d1da0f2feaedaae43a59f07f73786c05c668d..85fc36c6fd028e31db2ec087c6f96b58e760abef 100644
--- a/src/FLEXPART.f90
+++ b/src/FLEXPART.f90
@@ -44,6 +44,11 @@ program flexpart
use com_mod
use conv_mod
+#ifdef NETCDF_OUTPUT
+ use netcdf_output_mod, only: writeheader_netcdf
+#endif
+
+
implicit none
integer :: i,j,ix,jy,inest
@@ -302,6 +307,24 @@ program flexpart
! and open files that are to be kept open throughout the simulation
!******************************************************************
+ if (lnetcdfout.eq.1) then
+#ifdef NETCDF_OUTPUT
+ call writeheader_netcdf(lnest = .false.)
+#endif
+ else
+ call writeheader
+ end if
+
+ if (nested_output.eq.1) then
+ if (lnetcdfout.eq.1) then
+#ifdef NETCDF_OUTPUT
+ call writeheader_netcdf(lnest = .true.)
+#endif
+ else
+ call writeheader_nest
+ endif
+ endif
+
if (verbosity.gt.0) then
print*,'call writeheader'
endif
@@ -314,7 +337,9 @@ program flexpart
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
- !open(unitdates,file=path(2)(1:length(2))//'dates')
+ if (lnetcdfout.ne.1) then
+ open(unitdates,file=path(2)(1:length(2))//'dates')
+ end if
if (verbosity.gt.0) then
print*,'call openreceptors'
@@ -358,11 +383,11 @@ program flexpart
if (verbosity.gt.0) then
if (verbosity.gt.1) then
- CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
- write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+ call system_clock(count_clock, count_rate, count_max)
+ write(*,*) 'System clock',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
endif
if (info_flag.eq.1) then
- print*, 'info only mode (stop)'
+ print*, 'Info only mode (stop)'
stop
endif
print*,'call timemanager'
@@ -372,4 +397,11 @@ program flexpart
write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLEXPART MODEL RUN!'
+ ! output wall time
+ if (verbosity .gt. 0) then
+ call system_clock(count_clock,count_rate)
+ tins=(count_clock - count_clock0)/real(count_rate)
+ print*,'Wall time ',tins,'s, ',tins/60,'min, ',tins/3600,'h.'
+ endif
+
end program flexpart
diff --git a/src/com_mod.f90 b/src/com_mod.f90
index 8dd85748d2f51057f52cd1061989a72a6e7465b6..45b9e9823044784e895028e1806cce7921ad0524 100755
--- a/src/com_mod.f90
+++ b/src/com_mod.f90
@@ -112,6 +112,8 @@ module com_mod
! lconvection 1 if convection parameterization switched on, 0 if not
! lagespectra 1 if age spectra calculation switched on, 2 if not
+ integer :: lnetcdfout
+ ! lnetcdfout 1 for netcdf grid output, 0 if not. Set in COMMAND (namelist input)
integer :: nageclass,lage(maxageclass)
@@ -686,6 +688,7 @@ module com_mod
integer :: verbosity=0
integer :: info_flag=0
integer :: count_clock, count_clock0, count_rate, count_max
+ real :: tins
logical :: nmlout=.true.
diff --git a/src/makefile b/src/makefile
index 584e0265189b26662443341727c315d6b9904c46..4ade87968dc921d9e1b67c84332cd16d5ee85e1d 100644
--- a/src/makefile
+++ b/src/makefile
@@ -9,7 +9,6 @@ FFLAGS = -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4
#FFLAGS = -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
-
MODOBJS = \
par_mod.o com_mod.o \
conv_mod.o hanna_mod.o \
diff --git a/src/makefile.netcdf b/src/makefile.netcdf
new file mode 100644
index 0000000000000000000000000000000000000000..452ede95dc22fd3bc217d91993df1473c3dbacfb
--- /dev/null
+++ b/src/makefile.netcdf
@@ -0,0 +1,95 @@
+SHELL = /bin/bash
+MAIN = FP_ecmwf_gfortran_netcdf
+
+FC = gfortran
+INCPATH = /xnilu_wrk/flex_wrk/bin64/grib_api/include
+LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
+LIBPATH2 = /usr/lib/x86_64-linux-gnu/
+NCLIBPATH = /usr/local/lib
+FFLAGS = -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
+# netcdf output is included by using the complier switch -DNETCDF_OUTPUT and the C preprocessor -cpp
+FFLAGS = -O2 -cpp -DNETCDF_OUTPUT -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
+#FFLAGS = -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
+#LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
+LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -L$(NCLIBPATH) -lnetcdf -lnetcdff -lhdf5 -lgrib_api_f90 -lgrib_api -lm -ljasper
+
+MODOBJS = \
+par_mod.o com_mod.o \
+conv_mod.o hanna_mod.o \
+interpol_mod.o cmapf_mod.o \
+unc_mod.o oh_mod.o \
+xmass_mod.o flux_mod.o \
+point_mod.o outg_mod.o \
+netcdf_output_mod.o
+
+OBJECTS = \
+writeheader.o writeheader_txt.o writeheader_surf.o assignland.o\
+calcpar.o part0.o \
+caldate.o partdep.o \
+coordtrafo.o psih.o \
+raerod.o \
+drydepokernel.o random.o \
+erf.o readavailable.o \
+ew.o readcommand.o \
+advance.o readdepo.o \
+releaseparticles.o psim.o \
+FLEXPART.o readlanduse.o \
+getfields.o init_domainfill.o\
+interpol_wind.o readoutgrid.o \
+interpol_all.o readpaths.o \
+getrb.o readreceptors.o \
+getrc.o readreleases.o \
+getvdep.o readspecies.o \
+interpol_misslev.o readwind.o \
+conccalc.o richardson.o \
+concoutput.o concoutput_surf.o scalev.o \
+pbl_profile.o readOHfield.o\
+juldate.o timemanager.o \
+interpol_vdep.o interpol_rain.o \
+verttransform.o partoutput.o \
+hanna.o wetdepokernel.o \
+mean.o wetdepo.o \
+hanna_short.o windalign.o \
+obukhov.o gridcheck.o \
+hanna1.o initialize.o \
+ gridcheck_nests.o \
+readwind_nests.o calcpar_nests.o \
+verttransform_nests.o interpol_all_nests.o \
+interpol_wind_nests.o interpol_misslev_nests.o \
+interpol_vdep_nests.o interpol_rain_nests.o \
+getvdep_nests.o \
+readageclasses.o readpartpositions.o \
+calcfluxes.o fluxoutput.o \
+qvsat.o skplin.o \
+convmix.o calcmatrix.o \
+convect43c.o redist.o \
+sort2.o distance.o \
+centerofmass.o plumetraj.o \
+openouttraj.o calcpv.o \
+calcpv_nests.o distance2.o \
+clustering.o interpol_wind_short.o \
+interpol_wind_short_nests.o shift_field_0.o \
+shift_field.o outgrid_init.o \
+openreceptors.o boundcond_domainfill.o\
+partoutput_short.o readoutgrid_nest.o \
+outgrid_init_nest.o writeheader_nest.o writeheader_nest_surf.o \
+concoutput_nest.o concoutput_surf_nest.o wetdepokernel_nest.o \
+drydepokernel_nest.o zenithangle.o \
+ohreaction.o getvdep_nests.o \
+initial_cond_calc.o initial_cond_output.o \
+dynamic_viscosity.o get_settling.o
+
+
+$(MAIN): $(MODOBJS) $(OBJECTS)
+ $(FC) *.o -o $(MAIN) $(LDFLAGS)
+
+$(OBJECTS): $(MODOBJS)
+
+%.o: %.f90
+ $(FC) -c $(FFLAGS) $<
+
+clean:
+ rm *.o *.mod
+
+cleanall:
+ rm *.o *.mod $(MAIN)
diff --git a/src/readcommand.f90 b/src/readcommand.f90
index 0fe2b2c06f22712f8266fcdfe7abac08c14a3928..3adbbdbe3bc4dc68dedf89e95f807bd87a590e6e 100644
--- a/src/readcommand.f90
+++ b/src/readcommand.f90
@@ -108,6 +108,7 @@ subroutine readcommand
mquasilag, &
nested_output, &
linit_cond, &
+ lnetcdfout, &
surf_only
! Presetting namelist command
@@ -137,6 +138,7 @@ subroutine readcommand
mquasilag=0
nested_output=0
linit_cond=0
+ lnetcdfout=0
surf_only=0
! Open the command file and read user options
@@ -343,12 +345,25 @@ subroutine readcommand
mintime=lsynctime
endif
+! check for netcdf output switch (use for non-namelist input only!)
+ if (iout.ge.8) then
+ lnetcdfout = 1
+ iout = iout - 8
+#ifndef NETCDF_OUTPUT
+ print*,'ERROR: netcdf output not activated during compile time but used in COMMAND file!'
+ print*,'Please recompile with netcdf library or use standard output format.'
+ stop
+#endif
+ endif
+
! Check whether a valid option for gridded model output has been chosen
!**********************************************************************
if ((iout.lt.1).or.(iout.gt.5)) then
write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### '
- write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4, OR 5! #### '
+ write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4, OR 5 FOR #### '
+ write(*,*) ' #### STANDARD FLEXPART OUTPUT OR 9 - 13 #### '
+ write(*,*) ' #### FOR NETCDF OUTPUT #### '
stop
endif
@@ -362,7 +377,6 @@ subroutine readcommand
endif
-
! For quasilag output for each release is forbidden
!*****************************************************************************
diff --git a/src/timemanager.f90 b/src/timemanager.f90
index 2859eac11107a94a12f6fca669bb36cfe43ecef9..1d70ce482e4e7f41a3b095db3f0e7481812baac2 100644
--- a/src/timemanager.f90
+++ b/src/timemanager.f90
@@ -42,14 +42,16 @@ subroutine timemanager
! Changes Petra Seibert, Nov 2002 *
! call convection BEFORE new fields are read in BWD mode *
! Changes Caroline Forster, Feb 2005 *
- !new interface between flexpart and convection scheme *
- !Emanuel's latest subroutine convect43c.f is used *
+ ! new interface between flexpart and convection scheme *
+ ! Emanuel's latest subroutine convect43c.f is used *
+ ! Changes Stefan Henne, Harald Sodemann, 2013-2014 *
+ ! added netcdf output code *
!*****************************************************************************
! *
! Variables: *
- ! DEP .true. if either wet or dry deposition is switched on *
+ ! dep .true. if either wet or dry deposition is switched on *
! decay(maxspec) [1/s] decay constant for radioactive decay *
- ! DRYDEP .true. if dry deposition is switched on *
+ ! drydep .true. if dry deposition is switched on *
! ideltas [s] modelling period *
! itime [s] actual temporal position of calculation *
! ldeltat [s] time since computation of radioact. decay of depositions*
@@ -69,7 +71,7 @@ subroutine timemanager
! outnum number of samples for each concentration calculation *
! prob probability of absorption at ground due to dry *
! deposition *
- ! WETDEP .true. if wet deposition is switched on *
+ ! wetdep .true. if wet deposition is switched on *
! weight weight for each concentration sample (1/2 or 1) *
! uap(maxpart),ucp(maxpart),uzp(maxpart) = random velocities due to *
! turbulence *
@@ -91,6 +93,9 @@ subroutine timemanager
use oh_mod
use par_mod
use com_mod
+#ifdef NETCDF_OUTPUT
+ use netcdf_output_mod, only: concoutput_netcdf, concoutput_nest_netcdf,concoutput_surf_netcdf, concoutput_surf_nest_netcdf
+#endif
implicit none
@@ -128,6 +133,7 @@ subroutine timemanager
!**********************************************************************
+ itime=0
!write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
write(*,46) float(itime)/3600,itime,numpart
if (verbosity.gt.0) then
@@ -140,7 +146,6 @@ subroutine timemanager
do itime=0,ideltas,lsynctime
-
! Computation of wet deposition, OH reaction and mass transfer
! between two species every lsynctime seconds
! maybe wet depo frequency can be relaxed later but better be on safe side
@@ -151,17 +156,17 @@ subroutine timemanager
! changed by Petra Seibert 9/02
!********************************************************************
- if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then
+ if (wetdep .and. itime .ne. 0 .and. numpart .gt. 0) then
if (verbosity.gt.0) then
write (*,*) 'timemanager> call wetdepo'
endif
call wetdepo(itime,lsynctime,loutnext)
endif
- if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) &
+ if (ohrea .and. itime .ne. 0 .and. numpart .gt. 0) &
call ohreaction(itime,lsynctime,loutnext)
- if (ASSSPEC .and. itime .ne. 0 .and. numpart .gt. 0) then
+ if (assspec .and. itime .ne. 0 .and. numpart .gt. 0) then
stop 'associated species not yet implemented!'
! call transferspec(itime,lsynctime,loutnext)
endif
@@ -239,7 +244,7 @@ subroutine timemanager
! deposited mass radioactively decays
!***********************************************************************
- if (DEP.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
+ if (dep.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
do ks=1,nspec
do kp=1,maxpointspec_act
if (decay(ks).gt.0.) then
@@ -349,25 +354,50 @@ subroutine timemanager
if ((itime.eq.loutend).and.(outnum.gt.0.)) then
if ((iout.le.3.).or.(iout.eq.5)) then
if (surf_only.ne.1) then
- call concoutput(itime,outnum,gridtotalunc, &
- wetgridtotalunc,drygridtotalunc)
+ if (lnetcdfout.eq.1) then
+#ifdef NETCDF_OUTPUT
+ call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+#endif
+ else
+ call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+ endif
else
- if (verbosity.eq.1) then
- print*,'call concoutput_surf '
- CALL SYSTEM_CLOCK(count_clock)
- WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
- endif
- call concoutput_surf(itime,outnum,gridtotalunc, &
- wetgridtotalunc,drygridtotalunc)
- if (verbosity.eq.1) then
- print*,'called concoutput_surf '
- CALL SYSTEM_CLOCK(count_clock)
- WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
- endif
+ if (verbosity.eq.1) then
+ print*,'call concoutput_surf '
+ call system_clock(count_clock)
+ write(*,*) 'system clock',count_clock - count_clock0
+ endif
+ if (lnetcdfout.eq.1) then
+#ifdef NETCDF_OUTPUT
+ call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+#endif
+ else
+ call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+ if (verbosity.eq.1) then
+ print*,'called concoutput_surf '
+ call system_clock(count_clock)
+ write(*,*) 'system clock',count_clock - count_clock0
+ endif
+ endif
endif
- if ((nested_output.eq.1).and.(surf_only.ne.1)) call concoutput_nest(itime,outnum)
- if ((nested_output.eq.1).and.(surf_only.eq.1)) call concoutput_surf_nest(itime,outnum)
+ if (nested_output .eq. 1) then
+ if (lnetcdfout.eq.0) then
+ if (surf_only.ne.1) then
+ call concoutput_nest(itime,outnum)
+ else
+ call concoutput_surf_nest(itime,outnum)
+ endif
+ else
+#ifdef NETCDF_OUTPUT
+ if (surf_only.ne.1) then
+ call concoutput_nest_netcdf(itime,outnum)
+ else
+ call concoutput_surf_nest_netcdf(itime,outnum)
+ endif
+#endif
+ endif
+ endif
outnum=0.
endif
if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
@@ -473,8 +503,8 @@ subroutine timemanager
! Memorize particle positions
!****************************
- xold=xtra1(j)
- yold=ytra1(j)
+ xold=real(xtra1(j))
+ yold=real(ytra1(j))
zold=ztra1(j)
! Integrate Lagevin equation for lsynctime seconds
@@ -514,7 +544,7 @@ subroutine timemanager
decfact=1.
endif
- if (DRYDEPSPEC(ks)) then ! dry deposition
+ if (drydepspec(ks)) then ! dry deposition
drydeposit(ks)=xmass1(j,ks)*prob(ks)*decfact
xmass1(j,ks)=xmass1(j,ks)*(1.-prob(ks))*decfact
if (decay(ks).gt.0.) then ! correct for decay (see wetdepo)
@@ -525,7 +555,6 @@ subroutine timemanager
xmass1(j,ks)=xmass1(j,ks)*decfact
endif
-
if (mdomainfill.eq.0) then
if (xmass(npoint(j),ks).gt.0.) &
xmassfract=max(xmassfract,real(npart(npoint(j)))* &
@@ -537,25 +566,29 @@ subroutine timemanager
if (xmassfract.lt.0.0001) then ! terminate all particles carrying less mass
itra1(j)=-999999999
+ if (verbosity.gt.0) then
+ print*,'terminated particle ',j,' for small mass'
+ endif
endif
! Sabine Eckhardt, June 2008
! don't create depofield for backward runs
- if (DRYDEP.AND.(ldirect.eq.1)) then
- call drydepokernel(nclass(j),drydeposit,real(xtra1(j)), &
- real(ytra1(j)),nage,kp)
- if (nested_output.eq.1) call drydepokernel_nest( &
- nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)), &
- nage,kp)
+ if (drydep.AND.(ldirect.eq.1)) then
+ call drydepokernel(nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)),nage,kp)
+ if (nested_output.eq.1) then
+ call drydepokernel_nest(nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)),nage,kp)
+ endif
endif
! Terminate trajectories that are older than maximum allowed age
!***************************************************************
if (abs(itra1(j)-itramem(j)).ge.lage(nageclass)) then
- if (linit_cond.ge.1) &
- call initial_cond_calc(itime+lsynctime,j)
+ if (linit_cond.ge.1) call initial_cond_calc(itime+lsynctime,j)
itra1(j)=-999999999
+ if (verbosity.gt.0) then
+ print*,'terminated particle ',j,' for age'
+ endif
endif
endif
@@ -583,26 +616,26 @@ subroutine timemanager
!***************************
if (iflux.eq.1) then
- deallocate(flux)
+ deallocate(flux)
endif
- if (OHREA.eqv..TRUE.) then
- deallocate(OH_field,OH_field_height)
+ if (ohrea.eqv..TRUE.) then
+ deallocate(OH_field,OH_field_height)
endif
if (ldirect.gt.0) then
- deallocate(drygridunc,wetgridunc)
+ deallocate(drygridunc,wetgridunc)
endif
deallocate(gridunc)
deallocate(xpoint1,xpoint2,ypoint1,ypoint2,zpoint1,zpoint2,xmass)
deallocate(ireleasestart,ireleaseend,npart,kindz)
deallocate(xmasssave)
if (nested_output.eq.1) then
- deallocate(orooutn, arean, volumen)
- if (ldirect.gt.0) then
- deallocate(griduncn,drygriduncn,wetgriduncn)
- endif
+ deallocate(orooutn, arean, volumen)
+ if (ldirect.gt.0) then
+ deallocate(griduncn,drygriduncn,wetgriduncn)
+ endif
endif
deallocate(outheight,outheighthalf)
- deallocate(oroout, area, volume)
+ deallocate(oroout,area,volume)
end subroutine timemanager