from svn

r30 | hasod | 2014-10-21 10:08:00 +0200 (Tue, 21 Oct 2014) | 5 lines
Changed paths:
   M /trunk/src/FLEXPART.f90
   M /trunk/src/com_mod.f90
   M /trunk/src/makefile
   A /trunk/src/makefile.netcdf
   M /trunk/src/readcommand.f90
   M /trunk/src/timemanager.f90

ADD: Optional (compressed) netcdf output added. Activated via COMMAND
file. During compile time switches -DNETCDF_OUTPUT -cpp need to be
invoked. Compliation and linking is shown in makefile.netcdf

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

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.
    
 

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 \

src/makefile.netcdf

+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)

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
   !*****************************************************************************
 

src/timemanager.f90

@@ -44,12 +44,14 @@ subroutine timemanager
   !  Changes Caroline Forster, Feb 2005                                        *
   !   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   
+             call system_clock(count_clock)
+             write(*,*) 'system clock',count_clock - count_clock0   
             endif
-          call concoutput_surf(itime,outnum,gridtotalunc, &
-               wetgridtotalunc,drygridtotalunc)
+            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   
+                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
 
@@ -585,7 +618,7 @@ subroutine timemanager
   if (iflux.eq.1) then
     deallocate(flux)
   endif
-  if (OHREA.eqv..TRUE.) then
+  if (ohrea.eqv..TRUE.) then
     deallocate(OH_field,OH_field_height)
   endif
   if (ldirect.gt.0) then