Merge branch 'dev' into release-10

options/COMMAND

@@ -18,7 +18,7 @@
  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   
+ IPOUT=                 0, ! Particle position output: 0]no 1]every output 2]only at end 3]time averaged
  LSUBGRID=              0, ! Increase of ABL heights due to sub-grid scale orographic variations;[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  

src/FLEXPART.f90

@@ -67,13 +67,7 @@ program flexpart
   integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
   integer :: detectformat
 
-
-
-  ! Initialize arrays in com_mod
-  !*****************************
-  call com_mod_allocate_part(maxpart)
   
-
   ! Generate a large number of random numbers
   !******************************************
 
@@ -142,7 +136,7 @@ program flexpart
     print*,'nxshift=',nxshift 
     write(*,*) 'call readpaths'
   endif 
-  call readpaths(pathfile)
+  call readpaths
  
   if (verbosity.gt.1) then !show clock info 
      !print*,'length(4)',length(4)
@@ -171,6 +165,11 @@ program flexpart
     endif     
   endif
 
+  ! Initialize arrays in com_mod
+  !*****************************
+  call com_mod_allocate_part(maxpart)
+
+
   ! Read the age classes to be used
   !********************************
   if (verbosity.gt.0) then
@@ -452,7 +451,9 @@ program flexpart
      print*,'call timemanager'
   endif
 
+  if (verbosity.gt.0) write (*,*) 'timemanager> call wetdepo'
   call timemanager(metdata_format)
+ 
 
   if (verbosity.gt.0) then
 ! NIK 16.02.2005 
@@ -467,7 +468,6 @@ program flexpart
          write(*,*) '**********************************************'
       endif
     end do
-    write (*,*) 'timemanager> call wetdepo'
   endif
   
   write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&

src/FLEXPART_MPI.f90

@@ -76,12 +76,7 @@ program flexpart
 
   if (mp_measure_time) call mpif_mtime('flexpart',0)
 
-  ! Initialize arrays in com_mod 
-  !*****************************
-
-  if(.not.(lmpreader.and.lmp_use_reader)) call com_mod_allocate_part(maxpart_mpi)
-
-
+  
   ! Generate a large number of random numbers
   !******************************************
 
@@ -150,7 +145,7 @@ program flexpart
   if (verbosity.gt.0) then
     write(*,*) 'call readpaths'
   endif 
-  call readpaths(pathfile)
+  call readpaths
  
   if (verbosity.gt.1) then !show clock info 
      !print*,'length(4)',length(4)
@@ -179,6 +174,11 @@ program flexpart
     endif
   endif
 
+  ! Initialize arrays in com_mod 
+  !*****************************
+
+  if(.not.(lmpreader.and.lmp_use_reader)) call com_mod_allocate_part(maxpart_mpi)
+
 
 ! Read the age classes to be used
 !********************************
@@ -412,7 +412,7 @@ program flexpart
     if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
   end if ! (mpif_pid == 0) 
 
-  if (mp_measure_time) call mpif_mtime('iotime',0)
+  if (mp_measure_time) call mpif_mtime('iotime',1)
 
   if (verbosity.gt.0 .and. lroot) then
     print*,'call openreceptors'

src/com_mod.f90

@@ -18,6 +18,8 @@ module com_mod
 
   implicit none
 
+
+
   !****************************************************************
   ! Variables defining where FLEXPART input/output files are stored
   !****************************************************************
@@ -68,7 +70,7 @@ module com_mod
   ! outstep = real(abs(loutstep))
 
   real :: ctl,fine
-  integer :: ifine,iout,ipout,ipin,iflux,mdomainfill
+  integer :: ifine,iout,ipout,ipin,iflux,mdomainfill,ipoutfac
   integer :: mquasilag,nested_output,ind_source,ind_receptor
   integer :: ind_rel,ind_samp,ioutputforeachrelease,linit_cond,surf_only
   logical :: turbswitch
@@ -81,6 +83,7 @@ module com_mod
   ! iflux    flux calculation options: 1 calculation of fluxes, 2 no fluxes
   ! iout     output options: 1 conc. output (ng/m3), 2 mixing ratio (pptv), 3 both
   ! ipout    particle dump options: 0 no, 1 every output interval, 2 only at end
+  ! ipoutfac increase particle dump interval by factor (default 1)
   ! ipin     read in particle positions from dumped file from a previous run
   ! fine     real(ifine)
   ! mdomainfill 0: normal run
@@ -120,6 +123,9 @@ module com_mod
   integer :: lnetcdfout
   ! lnetcdfout   1 for netcdf grid output, 0 if not. Set in COMMAND (namelist input)
 
+  integer :: linversionout
+  ! linversionout 1 for one grid_time output file for each release containing all timesteps
+
   integer :: nageclass,lage(maxageclass)
 
   ! nageclass               number of ageclasses for the age spectra calculation
@@ -127,7 +133,6 @@ module com_mod
 
 
   logical :: gdomainfill
-
   ! gdomainfill             .T., if domain-filling is global, .F. if not
 
 !ZHG SEP 2015 wheather or not to read clouds from GRIB
@@ -173,7 +178,7 @@ 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)
+  real :: rm(maxspec),dryvel(maxspec),kao(maxspec)
   real :: ohcconst(maxspec),ohdconst(maxspec),ohnconst(maxspec)
 
   real :: area_hour(maxspec,24),point_hour(maxspec,24)
@@ -358,7 +363,9 @@ module com_mod
   real :: clwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !liquid   [kg/kg]
   real :: ciwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !ice      [kg/kg]
   real :: clw(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0  !combined [m3/m3]
-
+! RLT add pressure and dry air density
+  real :: prs(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
+  real :: rho_dry(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
   real :: pv(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
   real :: rho(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
   real :: drhodz(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
@@ -380,6 +387,7 @@ module com_mod
   ! uu,vv,ww [m/2]       wind components in x,y and z direction
   ! uupol,vvpol [m/s]    wind components in polar stereographic projection
   ! tt [K]               temperature data
+  ! prs                  air pressure
   ! qv                   specific humidity data
   ! pv (pvu)             potential vorticity
   ! rho [kg/m3]          air density
@@ -650,6 +658,7 @@ module com_mod
   real :: xreceptor(maxreceptor),yreceptor(maxreceptor)
   real :: receptorarea(maxreceptor)
   real :: creceptor(maxreceptor,maxspec)
+  real, allocatable, dimension(:,:) :: creceptor0
   character(len=16) :: receptorname(maxreceptor)
   integer :: numreceptor
 
@@ -673,6 +682,14 @@ module com_mod
   real, allocatable, dimension(:,:) :: xmass1
   real, allocatable, dimension(:,:) :: xscav_frac1
 
+! Variables used for writing out interval averages for partoutput
+!****************************************************************
+
+  integer, allocatable, dimension(:) :: npart_av
+  real, allocatable, dimension(:) :: part_av_cartx,part_av_carty,part_av_cartz,part_av_z,part_av_topo
+  real, allocatable, dimension(:) :: part_av_pv,part_av_qv,part_av_tt,part_av_rho,part_av_tro,part_av_hmix
+  real, allocatable, dimension(:) :: part_av_uu,part_av_vv,part_av_energy
+
   ! eso: Moved from timemanager
   real, allocatable, dimension(:) :: uap,ucp,uzp,us,vs,ws
   integer(kind=2), allocatable, dimension(:) :: cbt
@@ -779,13 +796,21 @@ contains
     allocate(itra1(nmpart),npoint(nmpart),nclass(nmpart),&
          & idt(nmpart),itramem(nmpart),itrasplit(nmpart),&
          & xtra1(nmpart),ytra1(nmpart),ztra1(nmpart),&
-         & xmass1(nmpart, maxspec),&
-         & checklifetime(nmpart,maxspec), species_lifetime(maxspec,2))!CGZ-lifetime
+         & xmass1(nmpart, maxspec))  ! ,&
+!         & checklifetime(nmpart,maxspec), species_lifetime(maxspec,2))!CGZ-lifetime
+
+    if (ipout.eq.3) then
+      allocate(npart_av(nmpart),part_av_cartx(nmpart),part_av_carty(nmpart),&
+           & part_av_cartz(nmpart),part_av_z(nmpart),part_av_topo(nmpart))
+      allocate(part_av_pv(nmpart),part_av_qv(nmpart),part_av_tt(nmpart),&
+           & part_av_rho(nmpart),part_av_tro(nmpart),part_av_hmix(nmpart))
+      allocate(part_av_uu(nmpart),part_av_vv(nmpart),part_av_energy(nmpart))
+    end if
 
 
     allocate(uap(nmpart),ucp(nmpart),uzp(nmpart),us(nmpart),&
          & vs(nmpart),ws(nmpart),cbt(nmpart))
-    
+
   end subroutine com_mod_allocate_part
 
 

src/concoutput.f90

@@ -71,6 +71,9 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
   integer :: sp_count_i,sp_count_r
   real :: sp_fact
   real :: outnum,densityoutrecept(maxreceptor),xl,yl
+! RLT
+  real :: densitydryrecept(maxreceptor)
+  real :: factor_dryrecept(maxreceptor)
 
 !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
 !    +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
@@ -105,6 +108,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
 ! mind        eso:added to ensure identical results between 2&3-fields versions
   character(LEN=8),save :: file_stat='REPLACE'
   logical :: ldates_file
+  logical :: lexist
   integer :: ierr
   character(LEN=100) :: dates_char
 
@@ -202,6 +206,9 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
 !      rho(iix,jjy,kzz-1,2)*dz2)/dz
         densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
              rho(iix,jjy,kzz-1,mind)*dz2)/dz
+! RLT
+        densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,mind)*dz1+ &
+             rho_dry(iix,jjy,kzz-1,mind)*dz2)/dz  
       end do
     end do
   end do
@@ -213,8 +220,14 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
     jjy=max(min(nint(yl),nymin1),0)
 !densityoutrecept(i)=rho(iix,jjy,1,2)
     densityoutrecept(i)=rho(iix,jjy,1,mind)
+! RLT
+    densitydryrecept(i)=rho_dry(iix,jjy,1,mind)
   end do
 
+! RLT
+! conversion factor for output relative to dry air
+  factor_drygrid=densityoutgrid/densitydrygrid
+  factor_dryrecept=densityoutrecept/densitydryrecept
 
 ! Output is different for forward and backward simulations
   do kz=1,numzgrid
@@ -352,7 +365,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
 
 ! Concentration output
 !*********************
-        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5).or.(iout.eq.6)) then
+        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
 
 ! Wet deposition
           sp_count_i=0
@@ -613,6 +626,49 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
 
   end do
 
+! RLT Aug 2017
+! Write out conversion factor for dry air
+  inquire(file=path(2)(1:length(2))//'factor_drygrid',exist=lexist)
+  if (lexist) then
+    ! open and append
+    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
+            status='old',action='write',access='append')
+  else
+    ! create new
+    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
+            status='new',action='write')
+  endif
+  sp_count_i=0
+  sp_count_r=0
+  sp_fact=-1.
+  sp_zer=.true.
+  do kz=1,numzgrid
+    do jy=0,numygrid-1
+      do ix=0,numxgrid-1
+        if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
+          if (sp_zer.eqv..true.) then ! first value not equal to one
+            sp_count_i=sp_count_i+1
+            sparse_dump_i(sp_count_i)= &
+                  ix+jy*numxgrid+kz*numxgrid*numygrid
+            sp_zer=.false.
+            sp_fact=sp_fact*(-1.)
+          endif
+          sp_count_r=sp_count_r+1
+          sparse_dump_r(sp_count_r)= &
+               sp_fact*factor_drygrid(ix,jy,kz)
+        else ! factor is one
+          sp_zer=.true.
+        endif
+      end do
+    end do
+  end do
+  write(unitoutfactor) sp_count_i
+  write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
+  write(unitoutfactor) sp_count_r
+  write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
+  close(unitoutfactor)
+
+
   if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
   if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
        wetgridtotal
@@ -639,7 +695,23 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
     end do
   endif
 
-
+! RLT Aug 2017
+! Write out conversion factor for dry air
+  if (numreceptor.gt.0) then
+    inquire(file=path(2)(1:length(2))//'factor_dryreceptor',exist=lexist)
+     if (lexist) then
+     ! open and append
+      open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
+              status='old',action='write',access='append')
+    else
+      ! create new
+      open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
+              status='new',action='write')
+    endif
+    write(unitoutfactor) itime
+    write(unitoutfactor) (factor_dryrecept(i),i=1,numreceptor)
+    close(unitoutfactor)
+  endif
 
 ! Reinitialization of grid
 !*************************

src/concoutput_mpi.f90

@@ -363,7 +363,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
 
 ! Concentration output
 !*********************
-        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5).or.(iout.eq.6)) then
+        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
 
 ! Wet deposition
           sp_count_i=0

src/concoutput_nest.f90

@@ -69,6 +69,9 @@ subroutine concoutput_nest(itime,outnum)
   integer :: sp_count_i,sp_count_r
   real :: sp_fact
   real :: outnum,densityoutrecept(maxreceptor),xl,yl
+! RLT
+  real :: densitydryrecept(maxreceptor)
+  real :: factor_dryrecept(maxreceptor)
 
   !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
   !    +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
@@ -95,6 +98,7 @@ subroutine concoutput_nest(itime,outnum)
   logical :: sp_zer
   character :: adate*8,atime*6
   character(len=3) :: anspec
+  logical :: lexist
   integer :: mind 
 ! mind        eso:added to ensure identical results between 2&3-fields versions
 
@@ -163,6 +167,9 @@ subroutine concoutput_nest(itime,outnum)
         !      rho(iix,jjy,kzz-1,2)*dz2)/dz
         densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
              rho(iix,jjy,kzz-1,mind)*dz2)/dz
+! RLT
+        densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,mind)*dz1+ &
+             rho_dry(iix,jjy,kzz-1,mind)*dz2)/dz
       end do
     end do
   end do
@@ -174,8 +181,14 @@ subroutine concoutput_nest(itime,outnum)
     jjy=max(min(nint(yl),nymin1),0)
     !densityoutrecept(i)=rho(iix,jjy,1,2)
     densityoutrecept(i)=rho(iix,jjy,1,mind)
+! RLT
+    densitydryrecept(i)=rho_dry(iix,jjy,1,mind)
   end do
 
+! RLT
+! conversion factor for output relative to dry air
+  factor_drygrid=densityoutgrid/densitydrygrid
+  factor_dryrecept=densityoutrecept/densitydryrecept
 
   ! Output is different for forward and backward simulations
     do kz=1,numzgrid
@@ -550,6 +563,47 @@ subroutine concoutput_nest(itime,outnum)
 
   end do
 
+! RLT Aug 2017
+! Write out conversion factor for dry air
+  inquire(file=path(2)(1:length(2))//'factor_drygrid_nest',exist=lexist)
+  if (lexist) then
+    ! open and append
+    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
+            status='old',action='write',access='append')
+  else
+    ! create new
+    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
+            status='new',action='write')
+  endif
+  sp_count_i=0
+  sp_count_r=0
+  sp_fact=-1.
+  sp_zer=.true.
+  do kz=1,numzgrid
+    do jy=0,numygridn-1
+      do ix=0,numxgridn-1
+        if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
+          if (sp_zer.eqv..true.) then ! first value not equal to one
+            sp_count_i=sp_count_i+1
+            sparse_dump_i(sp_count_i)= &
+                  ix+jy*numxgridn+kz*numxgridn*numygridn
+            sp_zer=.false.
+            sp_fact=sp_fact*(-1.)
+          endif
+          sp_count_r=sp_count_r+1
+          sparse_dump_r(sp_count_r)= &
+               sp_fact*factor_drygrid(ix,jy,kz)