Changes in / [f203036:941db73] in flexpart.git

Files:

• options/COMMAND (modified) (1 diff)
• src/FLEXPART.f90 (modified) (5 diffs)
• src/FLEXPART_MPI.f90 (modified) (4 diffs)
• src/com_mod.f90 (modified) (11 diffs)
• src/concoutput.f90 (modified) (7 diffs)
• src/concoutput_inversion.f90 (deleted)
• src/concoutput_inversion_nest.f90 (deleted)
• src/concoutput_mpi.f90 (modified) (1 diff)
• src/concoutput_nest.f90 (modified) (5 diffs)
• src/concoutput_surf.f90 (modified) (6 diffs)
• src/concoutput_surf_nest.f90 (modified) (13 diffs)
• src/getfields.f90 (modified) (3 diffs)
• src/init_domainfill.f90 (modified) (2 diffs)
• src/init_domainfill_mpi.f90 (modified) (4 diffs)
• src/initial_cond_output_inversion.f90 (deleted)
• src/makefile (modified) (8 diffs)
• src/mpi_mod.f90 (modified) (17 diffs)
• src/netcdf_output_mod.f90 (modified) (1 diff)
• src/outg_mod.f90 (modified) (1 diff)
• src/outgrid_init.f90 (modified) (1 diff)
• src/par_mod.f90 (modified) (5 diffs)
• src/partoutput.f90 (modified) (1 diff)
• src/partoutput_average.f90 (deleted)
• src/partoutput_average_mpi.f90 (deleted)
• src/partoutput_mpi.f90 (modified) (1 diff)
• src/partpos_average.f90 (deleted)
• src/readcommand.f90 (modified) (8 diffs)
• src/readpaths.f90 (modified) (1 diff)
• src/readwind_gfs.f90 (modified) (3 diffs)
• src/releaseparticles.f90 (modified) (1 diff)
• src/timemanager.f90 (modified) (5 diffs)
• src/timemanager_mpi.f90 (modified) (6 diffs)
• src/verttransform_ecmwf.f90 (modified) (10 diffs)
• src/verttransform_gfs.f90 (modified) (1 diff)

options/COMMAND

@@ -19 +19 @@
  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 3]time averaged
+ IPOUT=                 0, ! Particle position output: 0]no 1]every output 2]only at end   
  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    

src/FLEXPART.f90

@@ -68 +68 @@
   integer :: detectformat
 
+
+
+  ! Initialize arrays in com_mod
+  !*****************************
+  call com_mod_allocate_part(maxpart)
   
+
   ! Generate a large number of random numbers
   !******************************************
@@ -137 +143 @@
     write(*,*) 'call readpaths'
   endif 
-  call readpaths
+  call readpaths(pathfile)
  
   if (verbosity.gt.1) then !show clock info 
@@ -165 +171 @@
     endif     
   endif
-
-  ! Initialize arrays in com_mod
-  !*****************************
-  call com_mod_allocate_part(maxpart)
-
 
   ! Read the age classes to be used
@@ -452 +453 @@
   endif
 
-  if (verbosity.gt.0) write (*,*) 'timemanager> call wetdepo'
   call timemanager(metdata_format)
- 
 
   if (verbosity.gt.0) then
@@ -469 +468 @@
       endif
     end do
+    write (*,*) 'timemanager> call wetdepo'
   endif
   

src/FLEXPART_MPI.f90

@@ -77 +77 @@
   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
   !******************************************
@@ -146 +151 @@
     write(*,*) 'call readpaths'
   endif 
-  call readpaths
+  call readpaths(pathfile)
  
   if (verbosity.gt.1) then !show clock info 
@@ -174 +179 @@
     endif
   endif
-
-  ! Initialize arrays in com_mod 
-  !*****************************
-
-  if(.not.(lmpreader.and.lmp_use_reader)) call com_mod_allocate_part(maxpart_mpi)
 
 
@@ -413 +413 @@
   end if ! (mpif_pid == 0) 
 
-  if (mp_measure_time) call mpif_mtime('iotime',1)
+  if (mp_measure_time) call mpif_mtime('iotime',0)
 
   if (verbosity.gt.0 .and. lroot) then

src/com_mod.f90

@@ -18 +18 @@
 
   implicit none
-
-
 
   !****************************************************************
@@ -71 +69 @@
 
   real :: ctl,fine
-  integer :: ifine,iout,ipout,ipin,iflux,mdomainfill,ipoutfac
+  integer :: ifine,iout,ipout,ipin,iflux,mdomainfill
   integer :: mquasilag,nested_output,ind_source,ind_receptor
   integer :: ind_rel,ind_samp,ioutputforeachrelease,linit_cond,surf_only
@@ -84 +82 @@
   ! 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)
@@ -124 +121 @@
   ! 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)
 
@@ -134 +128 @@
 
   logical :: gdomainfill
+
   ! gdomainfill             .T., if domain-filling is global, .F. if not
 
@@ -179 +174 @@
   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 :: rm(maxspec),dryvel(maxspec)
   real :: ohcconst(maxspec),ohdconst(maxspec),ohnconst(maxspec)
 
@@ -364 +359 @@
   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)
@@ -388 +381 @@
   ! 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
@@ -659 +651 @@
   real :: receptorarea(maxreceptor)
   real :: creceptor(maxreceptor,maxspec)
-  real, allocatable, dimension(:,:) :: creceptor0
   character(len=16) :: receptorname(maxreceptor)
   integer :: numreceptor
@@ -682 +673 @@
   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
@@ -797 +780 @@
          & idt(nmpart),itramem(nmpart),itrasplit(nmpart),&
          & xtra1(nmpart),ytra1(nmpart),ztra1(nmpart),&
-         & 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
+         & xmass1(nmpart, maxspec),&
+         & checklifetime(nmpart,maxspec), species_lifetime(maxspec,2))!CGZ-lifetime
 
 
     allocate(uap(nmpart),ucp(nmpart),uzp(nmpart),us(nmpart),&
          & vs(nmpart),ws(nmpart),cbt(nmpart))
-
+    
   end subroutine com_mod_allocate_part
 

src/concoutput.f90

@@ -72 +72 @@
   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),
@@ -109 +106 @@
   character(LEN=8),save :: file_stat='REPLACE'
   logical :: ldates_file
-  logical :: lexist
   integer :: ierr
   character(LEN=100) :: dates_char
@@ -207 +203 @@
         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
@@ -221 +214 @@
 !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
@@ -366 +353 @@
 ! 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
@@ -627 +614 @@
   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/ &
@@ -696 +640 @@
   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

@@ -364 +364 @@
 ! 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

src/concoutput_nest.f90

@@ -70 +70 @@
   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),
@@ -99 +96 @@
   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
@@ -168 +164 @@
         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
@@ -182 +175 @@
     !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
@@ -564 +551 @@
   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)
-        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)
 
 

src/concoutput_surf.f90

@@ -72 +72 @@
   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),
@@ -104 +101 @@
   character :: adate*8,atime*6
   character(len=3) :: anspec
-  logical :: lexist
 
 
@@ -184 +180 @@
         densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
              rho(iix,jjy,kzz-1,2)*dz2)/dz
-! RLT
-        densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,2)*dz1+ &
-             rho_dry(iix,jjy,kzz-1,2)*dz2)/dz
       end do
     end do
@@ -197 +190 @@
     jjy=max(min(nint(yl),nymin1),0)
     densityoutrecept(i)=rho(iix,jjy,1,2)
-! RLT
-    densitydryrecept(i)=rho_dry(iix,jjy,1,2)
   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
@@ -609 +596 @@
   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,1
-    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/ &
@@ -678 +622 @@
   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_surf_nest.f90

@@ -70 +70 @@
   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),
@@ -99 +96 @@
   character :: adate*8,atime*6
   character(len=3) :: anspec
-  logical :: lexist
+
 
   ! Determine current calendar date, needed for the file name
@@ -162 +159 @@
         densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
              rho(iix,jjy,kzz-1,2)*dz2)/dz
-! RLT
-        densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,2)*dz1+ &
-             rho_dry(iix,jjy,kzz-1,2)*dz2)/dz
       end do
     end do
@@ -175 +169 @@
       jjy=max(min(nint(yl),nymin1),0)
       densityoutrecept(i)=rho(iix,jjy,1,2)
-! RLT
-    densitydryrecept(i)=rho_dry(iix,jjy,1,2)
     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
@@ -329 +317 @@
          write(unitoutgrid) sp_count_r
          write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgrid) sp_count_r
-!         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
+         write(unitoutgrid) sp_count_r
+         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
 
   ! Dry deposition
@@ -366 +354 @@
          write(unitoutgrid) sp_count_r
          write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgrid) sp_count_r
-!         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
+         write(unitoutgrid) sp_count_r
+         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
 
 
@@ -411 +399 @@
          write(unitoutgrid) sp_count_r
          write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgrid) sp_count_r
-!         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
+         write(unitoutgrid) sp_count_r
+         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
          else
 
@@ -452 +440 @@
          write(unitoutgrid) sp_count_r
          write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgrid) sp_count_r
-!         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
+         write(unitoutgrid) sp_count_r
+         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
          endif ! surf_only
 
@@ -499 +487 @@
          write(unitoutgridppt) sp_count_r
          write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgridppt) sp_count_r
-!         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
+         write(unitoutgridppt) sp_count_r
+         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
 
 
@@ -538 +526 @@
          write(unitoutgridppt) sp_count_r
          write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgridppt) sp_count_r
-!         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
+         write(unitoutgridppt) sp_count_r
+         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
 
 
@@ -582 +570 @@
          write(unitoutgridppt) sp_count_r
          write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgridppt) sp_count_r
-!         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
+         write(unitoutgridppt) sp_count_r
+         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
          else
 
@@ -623 +611 @@
          write(unitoutgridppt) sp_count_r
          write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
-!         write(unitoutgridppt) sp_count_r
-!         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
+         write(unitoutgridppt) sp_count_r
+         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
          endif ! surf_only
 
@@ -636 +624 @@
 
   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,1
-    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)
-        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)
 
 

src/getfields.f90

@@ -87 +87 @@
   real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
   real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
-! RLT added partial pressure water vapor 
-  real :: pwater(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
-  integer :: kz, ix
-  character(len=100) :: rowfmt
 
   integer :: indmin = 1
@@ -157 +153 @@
 40  indmin=indj
 
-    if (WETBKDEP) then
-      call writeprecip(itime,memind(1))
-    endif
+   if (WETBKDEP) then
+        call writeprecip(itime,memind(1))
+   endif
 
   else
@@ -208 +204 @@
 60  indmin=indj
 
-    if (WETBKDEP) then
-      call writeprecip(itime,memind(1))
-    endif
-
-  end if
-
-  ! RLT calculate dry air density
-  pwater=qv*prs/((r_air/r_water)*(1.-qv)+qv)
-  rho_dry=(prs-pwater)/(r_air*tt)
-
-  ! test density
-!  write(rowfmt,'(A,I6,A)') '(',nymax,'(E11.4,1X))'
-!  if(itime.eq.0) then
-!    open(500,file=path(2)(1:length(2))//'rho_dry.txt',status='replace',action='write')
-!    do kz=1,nzmax
-!      do ix=1,nxmax
-!        write(500,fmt=rowfmt) rho_dry(ix,:,kz,1)
-!      end do
-!    end do
-!    close(500)
-!    open(500,file=path(2)(1:length(2))//'rho.txt',status='replace',action='write')
-!    do kz=1,nzmax
-!      do ix=1,nxmax
-!        write(500,fmt=rowfmt) rho(ix,:,kz,1)
-!      end do
-!    end do
-!    close(500)
-!  endif
+   if (WETBKDEP) then
+        call writeprecip(itime,memind(1))
+   endif
+
+  endif
 
   lwindinterv=abs(memtime(2)-memtime(1))

src/init_domainfill.f90

@@ -86 +86 @@
     endif
   endif
-
-! Exit here if resuming a run from particle dump
-!***********************************************
-  if (gdomainfill.and.ipin.ne.0) return
 
 ! Do not release particles twice (i.e., not at both in the leftmost and rightmost
@@ -418 +414 @@
 !***************************************************************************
 
-  if ((ipin.eq.1).and.(.not.gdomainfill)) then
+  if (ipin.eq.1) then
     open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', &
          form='unformatted')

src/init_domainfill_mpi.f90

@@ -110 +110 @@
   endif
 
-! Exit here if resuming a run from particle dump
-!***********************************************
-  if (gdomainfill.and.ipin.ne.0) return
-
 ! Do not release particles twice (i.e., not at both in the leftmost and rightmost
 ! grid cell) for a global domain
@@ -217 +213 @@
         colmass(ix,jy)=(pp(1)-pp(nz))/ga*gridarea(jy)
         colmasstotal=colmasstotal+colmass(ix,jy)
+
       end do
     end do
@@ -469 +466 @@
 
 ! eso TODO: only needed for root process
-    if ((ipin.eq.1).and.(.not.gdomainfill)) then
+    if (ipin.eq.1) then
       open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', &
            form='unformatted')
@@ -477 +474 @@
     endif
 
-    if (ipin.eq.0) then    
-      numpart = numpart/mp_partgroup_np
-      if (mod(numpart,mp_partgroup_np).ne.0) numpart=numpart+1
-    end if
-
-  else ! Allocate dummy arrays for receiving processes
-    if (ipin.eq.0) then    
-      allocate(itra1_tmp(nullsize),npoint_tmp(nullsize),nclass_tmp(nullsize),&
-           & idt_tmp(nullsize),itramem_tmp(nullsize),itrasplit_tmp(nullsize),&
-           & xtra1_tmp(nullsize),ytra1_tmp(nullsize),ztra1_tmp(nullsize),&
-           & xmass1_tmp(nullsize, nullsize))
-    end if
+    numpart = numpart/mp_partgroup_np
+    if (mod(numpart,mp_partgroup_np).ne.0) numpart=numpart+1
+
+  else ! Allocate dummy arrays for receiving processes 
+    allocate(itra1_tmp(nullsize),npoint_tmp(nullsize),nclass_tmp(nullsize),&
+         & idt_tmp(nullsize),itramem_tmp(nullsize),itrasplit_tmp(nullsize),&
+         & xtra1_tmp(nullsize),ytra1_tmp(nullsize),ztra1_tmp(nullsize),&
+         & xmass1_tmp(nullsize, nullsize))
     
-  end if ! end if(lroot)
-
+  end if ! end if(lroot)  
 
 
 ! Distribute particles to other processes (numpart is 'per-process', not total)
-! Only if not restarting from previous run
-  if (ipin.eq.0) then
-    call MPI_Bcast(numpart, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
-    call mpif_send_part_properties(npart(1)/mp_partgroup_np)
+  call MPI_Bcast(numpart, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
+! eso TODO: xmassperparticle: not necessary to send
+  call MPI_Bcast(xmassperparticle, 1, mp_sp, id_root, mp_comm_used, mp_ierr)
+  call mpif_send_part_properties(numpart)
 
 ! Deallocate the temporary arrays used for all particles
-    deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,&
+  deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,&
          & itrasplit_tmp,xtra1_tmp,ytra1_tmp,ztra1_tmp,xmass1_tmp)
-  end if
 
 

src/makefile

@@ -118 +118 @@
 OBJECTS_SERIAL = \
         releaseparticles.o      partoutput.o \
-        partoutput_average.o \
         conccalc.o \
         init_domainfill.o       concoutput.o  \
@@ -128 +127 @@
         redist.o                \
         concoutput_surf.o       concoutput_surf_nest.o  \
-        concoutput_inversion_nest.o     \
-        concoutput_inversion.o \
         getfields.o \
         readwind_ecmwf.o
@@ -135 +132 @@
 ## For MPI version 
 OBJECTS_MPI = releaseparticles_mpi.o partoutput_mpi.o \
-        partoutput_average_mpi.o conccalc_mpi.o \
+        conccalc_mpi.o \
         init_domainfill_mpi.o concoutput_mpi.o  \
         timemanager_mpi.o FLEXPART_MPI.o        \
@@ -152 +149 @@
 advance.o               initialize.o            \
 writeheader.o           writeheader_txt.o       \
-partpos_average.o       writeprecip.o \
+writeprecip.o \
 writeheader_surf.o      assignland.o\
 part0.o                 gethourlyOH.o\
@@ -202 +199 @@
 drydepokernel_nest.o    zenithangle.o \
 ohreaction.o            getvdep_nests.o \
-initial_cond_calc.o     initial_cond_output.o initial_cond_output_inversion.o \
+initial_cond_calc.o     initial_cond_output.o \
 dynamic_viscosity.o     get_settling.o  \
 initialize_cbl_vel.o    re_initialize_particle.o \
@@ -274 +271 @@
 conccalc_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o unc_mod.o
 concoutput.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
-concoutput_inversion.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
 concoutput_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o point_mod.o \
         unc_mod.o mean_mod.o
 concoutput_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
-concoutput_inversion_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
 concoutput_nest_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o point_mod.o \
         unc_mod.o mean_mod.o
@@ -325 +320 @@
 initial_cond_calc.o: com_mod.o outg_mod.o par_mod.o unc_mod.o
 initial_cond_output.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
-initial_cond_output_inversion.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
 initialize.o: com_mod.o hanna_mod.o interpol_mod.o par_mod.o random_mod.o
 initialize_cbl_vel.o: com_mod.o par_mod.o random_mod.o
@@ -354 +348 @@
 part0.o: par_mod.o
 partdep.o: par_mod.o
-partpos_average.o: com_mod.o par_mod.o
 partoutput.o: com_mod.o par_mod.o
-partoutput_average.o: com_mod.o par_mod.o
-partoutput_average_mpi.o: com_mod.o par_mod.o mpi_mod.o
 partoutput_mpi.o: com_mod.o mpi_mod.o par_mod.o
 partoutput_short.o: com_mod.o par_mod.o

src/mpi_mod.f90

@@ -88 +88 @@
 ! Variables for MPI processes in the 'particle' group
   integer, allocatable, dimension(:) :: mp_partgroup_rank
-  integer, allocatable, dimension(:) :: npart_per_process
   integer :: mp_partgroup_comm, mp_partgroup_pid, mp_partgroup_np
 
@@ -126 +125 @@
 ! mp_time_barrier   Measure MPI barrier time
 ! mp_exact_numpart  Use an extra MPI communication to give the exact number of particles
-!                   to standard output (this does not otherwise affect the simulation)
+!                   to standard output (this does *not* otherwise affect the simulation) 
   logical, parameter :: mp_dbg_mode = .false.
   logical, parameter :: mp_dev_mode = .false.
@@ -191 +190 @@
 !   mp_np       number of running processes, decided at run-time
 !***********************************************************************
-    use par_mod, only: maxpart, numwfmem, dep_prec, maxreceptor, maxspec
-    use com_mod, only: mpi_mode, verbosity, creceptor0
+    use par_mod, only: maxpart, numwfmem, dep_prec
+    use com_mod, only: mpi_mode, verbosity
 
     implicit none
@@ -338 +337 @@
 
 ! Set maxpart per process
-! ESO 08/2016: Increase maxpart per process, in case of unbalanced distribution
+! eso 08/2016: Increase maxpart per process, in case of unbalanced distribution
     maxpart_mpi=int(mp_maxpart_factor*real(maxpart)/real(mp_partgroup_np))
     if (mp_np == 1) maxpart_mpi = maxpart
@@ -366 +365 @@
     end if
 
-! Allocate array for number of particles per process    
-    allocate(npart_per_process(0:mp_partgroup_np-1))
-
-! Write whether MPI_IN_PLACE is used or not
-#ifdef USE_MPIINPLACE
-    if (lroot) write(*,*) 'Using MPI_IN_PLACE operations'
-#else
-    if (lroot) allocate(creceptor0(maxreceptor,maxspec))
-    if (lroot) write(*,*) 'Not using MPI_IN_PLACE operations'
-#endif
     goto 101
 
@@ -570 +559 @@
 ! invalid particles at the end of the arrays
 
-601 do i=numpart, 1, -1
+601 do i=num_part, 1, -1
       if (itra1(i).eq.-999999999) then
         numpart=numpart-1
@@ -609 +598 @@
     integer :: i,jj,nn, num_part=1,m,imin, num_trans
     logical :: first_iter
-    integer,allocatable,dimension(:) :: idx_arr
+    integer,allocatable,dimension(:) :: numparticles_mpi, idx_arr
     real,allocatable,dimension(:) :: sorted ! TODO: we don't really need this
 
@@ -618 +607 @@
 ! All processes exchange information on number of particles
 !****************************************************************************
-    allocate( idx_arr(0:mp_partgroup_np-1), sorted(0:mp_partgroup_np-1))
-
-    call MPI_Allgather(numpart, 1, MPI_INTEGER, npart_per_process, &
+    allocate(numparticles_mpi(0:mp_partgroup_np-1), &
+         &idx_arr(0:mp_partgroup_np-1), sorted(0:mp_partgroup_np-1))
+
+    call MPI_Allgather(numpart, 1, MPI_INTEGER, numparticles_mpi, &
          & 1, MPI_INTEGER, mp_comm_used, mp_ierr)
 
@@ -626 +616 @@
 ! Sort from lowest to highest
 ! Initial guess: correct order
-    sorted(:) = npart_per_process(:)
+    sorted(:) = numparticles_mpi(:)
     do i=0, mp_partgroup_np-1
       idx_arr(i) = i
     end do
-
-! Do not rebalance particles for ipout=3    
-    if (ipout.eq.3) return
 
 ! For each successive element in index array, see if a lower value exists
@@ -658 +645 @@
     m=mp_partgroup_np-1 ! index for last sorted process (most particles)
     do i=0,mp_partgroup_np/2-1
-      num_trans = npart_per_process(idx_arr(m)) - npart_per_process(idx_arr(i))
+      num_trans = numparticles_mpi(idx_arr(m)) - numparticles_mpi(idx_arr(i))
       if (mp_partid.eq.idx_arr(m).or.mp_partid.eq.idx_arr(i)) then
-        if ( npart_per_process(idx_arr(m)).gt.mp_min_redist.and.&
-             & real(num_trans)/real(npart_per_process(idx_arr(m))).gt.mp_redist_fract) then
+        if ( numparticles_mpi(idx_arr(m)).gt.mp_min_redist.and.&
+             & real(num_trans)/real(numparticles_mpi(idx_arr(m))).gt.mp_redist_fract) then
 ! DBG
-          ! write(*,*) 'mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, npart_per_process', &
-          !      &mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, npart_per_process
+          ! write(*,*) 'mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, numparticles_mpi', &
+          !      &mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, numparticles_mpi
 ! DBG
           call mpif_redist_part(itime, idx_arr(m), idx_arr(i), num_trans/2)
@@ -672 +659 @@
     end do
 
-    deallocate(idx_arr, sorted)
+    deallocate(numparticles_mpi, idx_arr, sorted)
 
   end subroutine mpif_calculate_part_redist
@@ -1974 +1961 @@
     if (readclouds) then
       j=j+1
-      call MPI_Irecv(ctwc(:,:,mind),d2s1*5,mp_sp,id_read,MPI_ANY_TAG,&
+      call MPI_Irecv(ctwc(:,:,mind),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
            &MPI_COMM_WORLD,reqs(j),mp_ierr)
       if (mp_ierr /= 0) goto 600
@@ -2339 +2326 @@
       if (readclouds) then
         j=j+1
-        call MPI_Irecv(ctwcn(:,:,mind,k),d2s1*5,mp_sp,id_read,MPI_ANY_TAG,&
+        call MPI_Irecv(ctwcn(:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
              &MPI_COMM_WORLD,reqs(j),mp_ierr)
         if (mp_ierr /= 0) goto 600
@@ -2475 +2462 @@
     end if
 
-  ! Receptor concentrations    
+#else
+
+      call MPI_Reduce(gridunc, gridunc0, grid_size3d, mp_sp, MPI_SUM, id_root, &
+           & mp_comm_used, mp_ierr)
+      if (lroot) gridunc = gridunc0
+
+#endif
+
+    if ((WETDEP).and.(ldirect.gt.0)) then
+      call MPI_Reduce(wetgridunc, wetgridunc0, grid_size2d, mp_cp, MPI_SUM, id_root, &
+           & mp_comm_used, mp_ierr)
+      if (mp_ierr /= 0) goto 600
+    end if
+
+    if ((DRYDEP).and.(ldirect.gt.0)) then
+      call MPI_Reduce(drygridunc, drygridunc0, grid_size2d, mp_cp, MPI_SUM, id_root, &
+           & mp_comm_used, mp_ierr)
+      if (mp_ierr /= 0) goto 600
+    end if
+
+! Receptor concentrations    
     if (lroot) then
       call MPI_Reduce(MPI_IN_PLACE,creceptor,rcpt_size,mp_sp,MPI_SUM,id_root, &
@@ -2484 +2491 @@
            & mp_comm_used,mp_ierr)
     end if
-
-#else
-
-      call MPI_Reduce(gridunc, gridunc0, grid_size3d, mp_sp, MPI_SUM, id_root, &
-           & mp_comm_used, mp_ierr)
-      if (mp_ierr /= 0) goto 600
-      if (lroot) gridunc = gridunc0
-
-      call MPI_Reduce(creceptor, creceptor0,rcpt_size,mp_sp,MPI_SUM,id_root, &
-           & mp_comm_used,mp_ierr)
-      if (mp_ierr /= 0) goto 600
-      if (lroot) creceptor = creceptor0
-
-#endif
-
-    if ((WETDEP).and.(ldirect.gt.0)) then
-      call MPI_Reduce(wetgridunc, wetgridunc0, grid_size2d, mp_cp, MPI_SUM, id_root, &
-           & mp_comm_used, mp_ierr)
-      if (mp_ierr /= 0) goto 600
-    end if
-
-    if ((DRYDEP).and.(ldirect.gt.0)) then
-      call MPI_Reduce(drygridunc, drygridunc0, grid_size2d, mp_cp, MPI_SUM, id_root, &
-           & mp_comm_used, mp_ierr)
-      if (mp_ierr /= 0) goto 600
-    end if
-
 
     if (mp_measure_time) call mpif_mtime('commtime',1)
@@ -2720 +2700 @@
       end if
 
-   case ('readwind')
-     if (imode.eq.0) then
-       call cpu_time(mp_readwind_time_beg)
-       mp_readwind_wtime_beg = mpi_wtime()
-     else
-       call cpu_time(mp_readwind_time_end)
-       mp_readwind_wtime_end = mpi_wtime()
-
-       mp_readwind_time_total = mp_readwind_time_total + &
-            &(mp_readwind_time_end - mp_readwind_time_beg)
-       mp_readwind_wtime_total = mp_readwind_wtime_total + &
-            &(mp_readwind_wtime_end - mp_readwind_wtime_beg)
-     end if
+    case ('readwind')
+      if (imode.eq.0) then
+        call cpu_time(mp_readwind_time_beg)
+        mp_readwind_wtime_beg = mpi_wtime()
+      else
+        call cpu_time(mp_readwind_time_end)
+        mp_readwind_wtime_end = mpi_wtime()
+
+        mp_readwind_time_total = mp_readwind_time_total + &
+             &(mp_readwind_time_end - mp_readwind_time_beg)
+        mp_readwind_wtime_total = mp_readwind_wtime_total + &
+             &(mp_readwind_wtime_end - mp_readwind_wtime_beg)
+      end if
 
     case ('commtime')
@@ -2808 +2788 @@
           write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR GETFIELDS:',&
                & mp_getfields_time_total
-!          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR READWIND:',&
-!               & mp_readwind_wtime_total
-!          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR READWIND:',&
-!               & mp_readwind_time_total
+          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR READWIND:',&
+               & mp_readwind_wtime_total
+          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR READWIND:',&
+               & mp_readwind_time_total
           write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR FILE IO:',&
                & mp_io_wtime_total

src/netcdf_output_mod.f90

@@ -93 +93 @@
   character(len=255), parameter :: institution = 'NILU'
 
-  integer            :: tpointer=0
+  integer            :: tpointer
   character(len=255) :: ncfname, ncfnamen
 

src/outg_mod.f90

@@ -37 +37 @@
   real,allocatable, dimension (:,:,:) :: areanorth
   real,allocatable, dimension (:,:,:) :: densityoutgrid
-  real,allocatable, dimension (:,:,:) :: densitydrygrid ! added RLT 
-  real,allocatable, dimension (:,:,:) :: factor_drygrid ! added RLT 
   real,allocatable, dimension (:,:,:) :: factor3d
   real,allocatable, dimension (:,:,:) :: grid

src/outgrid_init.f90

@@ -268 +268 @@
     if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
   allocate(densityoutgrid(0:max(numxgrid,numxgridn)-1, &
-       0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
-! RLT
-  allocate(densitydrygrid(0:max(numxgrid,numxgridn)-1, &
-       0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
-  allocate(factor_drygrid(0:max(numxgrid,numxgridn)-1, &
        0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
     if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'

src/par_mod.f90

@@ -30 +30 @@
 !        Update 15 August 2013 IP                                              *
 !                                                                              *
+!        ESO 2016:                                                             *
+!          GFS specific parameters moved to gfs_mod.f90                        *
+!          ECMWF specific parameters moved to ecmwf_mod.f90                    *
 !                                                                              *
 !*******************************************************************************
@@ -77 +80 @@
   real,parameter :: pi=3.14159265, r_earth=6.371e6, r_air=287.05, ga=9.81
   real,parameter :: cpa=1004.6, kappa=0.286, pi180=pi/180., vonkarman=0.4
-  ! additional constants RLT Aug-2017
-  real,parameter :: rgas=8.31447 
-  real,parameter :: r_water=461.495
 
   ! pi                      number "pi"
@@ -89 +89 @@
   ! kappa                   exponent of formula for potential temperature
   ! vonkarman               von Karman constant
-  ! rgas                    universal gas constant [J/mol/K]
-  ! r_water                 specific gas constant for water vapor [J/kg/K]
 
   real,parameter :: karman=0.40, href=15., convke=2.0
@@ -282 +280 @@
 
   integer,parameter :: unitpath=1, unitcommand=1, unitageclasses=1, unitgrid=1
-  integer,parameter :: unitavailab=1, unitreleases=88, unitpartout=93, unitpartout_average=105
+  integer,parameter :: unitavailab=1, unitreleases=88, unitpartout=93
   integer,parameter :: unitpartin=93, unitflux=98, unitouttraj=96
   integer,parameter :: unitvert=1, unitoro=1, unitpoin=1, unitreceptor=1
@@ -292 +290 @@
   integer,parameter :: unitboundcond=89
   integer,parameter :: unittmp=101
-! RLT
-  integer,parameter :: unitoutfactor=102
 
 !******************************************************

src/partoutput.f90

@@ -71 +71 @@
   !**************************************
 
-  if (ipout.eq.1.or.ipout.eq.3) then
+  if (ipout.eq.1) then
     open(unitpartout,file=path(2)(1:length(2))//'partposit_'//adate// &
          atime,form='unformatted')

src/partoutput_mpi.f90

@@ -78 +78 @@
   !**************************************
 
-  if (ipout.eq.1.or.ipout.eq.3) then
+  if (ipout.eq.1) then
     open(unitpartout,file=path(2)(1:length(2))//'partposit_'//adate// &
          atime,form='unformatted',status=file_stat,position='append')

src/readcommand.f90

@@ -50 +50 @@
   ! ipin                 1 continue simulation with dumped particle data, 0 no *
   ! ipout                0 no particle dump, 1 every output time, 3 only at end*
-  ! ipoutfac             increase particle dump interval by factor (default 1) *
   ! itsplit [s]          time constant for particle splitting                  *
   ! loutaver [s]         concentration output is an average over loutaver      *
@@ -98 +97 @@
   iout, &
   ipout, &
-  ipoutfac, &
   lsubgrid, &
   lconvection, &
@@ -114 +112 @@
   surf_only, &
   cblflag, &
-  linversionout, &
   ohfields_path
 
@@ -132 +129 @@
   iout=3
   ipout=0
-  ipoutfac=1
   lsubgrid=1
   lconvection=1
@@ -148 +144 @@
   surf_only=0 
   cblflag=0 ! if using old-style COMMAND file, set to 1 here to use mc cbl routine
-  linversionout=0
   ohfields_path="../../flexin/"
 
@@ -416 +411 @@
   !**********************************************************************
 
-  if ((iout.lt.1).or.(iout.gt.5)) then
+  if ((iout.lt.1).or.(iout.gt.6)) then
     write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND:     #### '
     write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4 OR 5 FOR        #### '
@@ -476 +471 @@
       write(*,*) '#### FOR DOMAIN FILLING RUNS OUTPUT FOR      ####'
       write(*,*) '#### EACH RELEASE IS FORBIDDEN !             ####'
-      stop
-  endif
-
-  ! Inversion output format only for backward runs
-  !*****************************************************************************
-  
-  if ((linversionout.eq.1).and.(ldirect.eq.1)) then
-      write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND:     ####'
-      write(*,*) '#### INVERSION OUTPUT FORMAT ONLY FOR        ####'
-      write(*,*) '#### BACKWARD RUNS                           ####'
       stop
   endif
@@ -522 +507 @@
   !****************************************************************
 
-  if ((ipout.ne.0).and.(ipout.ne.1).and.(ipout.ne.2).and.(ipout.ne.3)) then
+  if ((ipout.ne.0).and.(ipout.ne.1).and.(ipout.ne.2)) then
     write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND:     #### '
-    write(*,*) ' #### IPOUT MUST BE 0, 1, 2 OR 3!             #### '
+    write(*,*) ' #### IPOUT MUST BE 1, 2 OR 3!                #### '
     stop
   endif

src/readpaths.f90

@@ -20 +20 @@
 !**********************************************************************
 
-subroutine readpaths
+subroutine readpaths !(pathfile)
 
   !*****************************************************************************

src/readwind_gfs.f90

@@ -83 +83 @@
 
   ! NCEP
-  integer :: numpt,numpu,numpv,numpw,numprh,numpclwch
+  integer :: numpt,numpu,numpv,numpw,numprh
   real :: help, temp, ew
   real :: elev
@@ -134 +134 @@
   numpw=0
   numprh=0
-  numpclwch=0
   ifield=0
 10   ifield=ifield+1
@@ -558 +557 @@
       endif
 ! SEC & IP 12/2018 read GFS clouds
-      if((isec1(6).eq.153).and.(isec1(7).eq.100)) then  !! CLWCR  Cloud liquid water content [kg/kg] 
-         if((i.eq.0).and.(j.eq.0)) then
-            do ii=1,nuvz
-              if ((isec1(8)*100.0).eq.akz(ii)) numpclwch=ii
-            end do
-        endif
-        help=zsec4(nxfield*(ny-j-1)+i+1)
-        if(i.le.i180) then
-          clwch(i179+i,j,numpclwch,n)=help
-        else
-          clwch(i-i181,j,numpclwch,n)=help
-        endif
+      if(isec1(6).eq.153) then  !! CLWCR  Cloud liquid water content [kg/kg]
+        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
         readclouds=.true.
         sumclouds=.true.
-!        readclouds=.false.
-!       sumclouds=.false.
       endif
 

src/releaseparticles.f90

@@ -114 +114 @@
       average_timecorrect=0.
       do k=1,nspec
-        if(abs(xpoint2(i)-xpoint1(i)).lt.1.E-4.and.abs(ypoint2(i)-ypoint1(i)).lt.1.E-4) then
-!        if (zpoint1(i).gt.0.5) then      ! point source
+        if (zpoint1(i).gt.0.5) then      ! point source
           timecorrect(k)=point_hour(k,nhour)*point_dow(k,ndayofweek)
         else                             ! area source

src/timemanager.f90

@@ -408 +408 @@
 #endif
             else
-              if (linversionout.eq.1) then
-                call concoutput_inversion(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
-                if (verbosity.eq.1) then
-                  print*,'called concoutput_inversion'
-                  call system_clock(count_clock)
-                  write(*,*) 'system clock',count_clock - count_clock0 
-                endif
-              else
               call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
-              endif
               if (verbosity.eq.1) then
                 print*,'called concoutput_surf '
@@ -431 +422 @@
                 call concoutput_nest(itime,outnum)
               else 
-                if(linversionout.eq.1) then
-                  call concoutput_inversion_nest(itime,outnum)
-                else 
                 call concoutput_surf_nest(itime,outnum)
-              endif
               endif
             else
@@ -464 +451 @@
 45      format(i13,' Seconds simulated: ',i13, ' Particles:    Uncertainty: ',3f7.3)
 46      format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
-        if (ipout.ge.1) then
-          if (mod(itime,ipoutfac*loutstep).eq.0) call partoutput(itime) ! dump particle positions
-          if (ipout.eq.3) call partoutput_average(itime) ! dump particle positions
-        endif
+        if (ipout.ge.1) call partoutput(itime)    ! dump particle positions
         loutnext=loutnext+loutstep
         loutstart=loutnext-loutaver/2
@@ -625 +609 @@
 !        write (*,*) 'advance: ',prob(1),xmass1(j,1),ztra1(j)
 
-  ! Calculate average position for particle dump output
-  !****************************************************
-
-        if (ipout.eq.3) call partpos_average(itime,j)
-
-
   ! Calculate the gross fluxes across layer interfaces
   !***************************************************
@@ -752 +730 @@
   if (ipout.eq.2) call partoutput(itime)     ! dump particle positions
 
-  if (linit_cond.ge.1) then
-    if(linversionout.eq.1) then
-      call initial_cond_output_inversion(itime)   ! dump initial cond. field
-    else
-      call initial_cond_output(itime)   ! dump initial cond. fielf
-    endif
-  endif
+  if (linit_cond.ge.1) call initial_cond_output(itime)   ! dump initial cond. field
 
   !close(104)

src/timemanager_mpi.f90

@@ -113 +113 @@
   integer :: j,ks,kp,l,n,itime=0,nstop,nstop1,memstat=0
 ! integer :: ksp
-  integer :: ip,irec
+  integer :: ip
   integer :: loutnext,loutstart,loutend
   integer :: ix,jy,ldeltat,itage,nage,idummy
@@ -129 +129 @@
 ! Measure time spent in timemanager
   if (mp_measure_time) call mpif_mtime('timemanager',0)
-
 
 ! First output for time 0
@@ -533 +532 @@
                 end if
 
-              else
+              else  ! :TODO: check for zeroing in the netcdf module
                 call concoutput_surf_nest(itime,outnum)
               end if
@@ -594 +593 @@
 46      format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
         if (ipout.ge.1) then
-          if (mp_measure_time) call mpif_mtime('iotime',0)
-          irec=0
           do ip=0, mp_partgroup_np-1
-            if (ip.eq.mp_partid) then
-              if (mod(itime,ipoutfac*loutstep).eq.0) call partoutput(itime) ! dump particle positions
-              if (ipout.eq.3) call partoutput_average(itime,irec) ! dump particle positions
-            endif
-            if (ipout.eq.3) irec=irec+npart_per_process(ip)
+            if (ip.eq.mp_partid) call partoutput(itime) ! dump particle positions
             call mpif_mpi_barrier
           end do
-          if (mp_measure_time) call mpif_mtime('iotime',1)
         end if
 
@@ -765 +757 @@
         if (mp_measure_time) call mpif_mtime('advance',1)
 
-  ! Calculate average position for particle dump output
-  !****************************************************
-
-        if (ipout.eq.3) call partpos_average(itime,j)
-
 
 ! Calculate the gross fluxes across layer interfaces
@@ -908 +895 @@
     do ip=0, mp_partgroup_np-1
       if (ip.eq.mp_partid) then 
-        if (mp_dbg_mode) write(*,*) 'call partoutput(itime), proc, mp_partid',ip,mp_partid
+        !if (mp_dbg_mode) write(*,*) 'call partoutput(itime), proc, mp_partid',ip,mp_partid
         call partoutput(itime)    ! dump particle positions
       end if

src/verttransform_ecmwf.f90

@@ -73 +73 @@
   use com_mod
   use cmapf_mod, only: cc2gll
+!  use mpi_mod
 
   implicit none
@@ -81 +82 @@
   real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: rhoh,uvzlev,wzlev
   real,dimension(0:nxmax-1,0:nymax-1,nzmax) :: pinmconv
-  ! RLT added pressure
-  real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: prsh
   real,dimension(0:nxmax-1,0:nymax-1) ::  tvold,pold,pint,tv
   real,dimension(0:nymax-1) :: cosf
@@ -220 +219 @@
 !*************************
 
+
   do jy=0,nymin1
     do ix=0,nxmin1
@@ -230 +230 @@
   wzlev(:,:,1)=0.
   rhoh(:,:,1)=pold/(r_air*tvold)
-  ! RLT add pressure
-  prsh(:,:,1)=ps(:,:,1,n)
 
 
@@ -239 +237 @@
   do kz=2,nuvz
     pint=akz(kz)+bkz(kz)*ps(:,:,1,n)
-    ! RLT add pressure
-    prsh(:,:,kz)=pint
     tv=tth(:,:,kz,n)*(1.+0.608*qvh(:,:,kz,n))
     rhoh(:,:,kz)=pint(:,:)/(r_air*tv)
@@ -292 +288 @@
   pv(:,:,1,n)=pvh(:,:,1)
   rho(:,:,1,n)=rhoh(:,:,1)
-! RLT add pressure
-  prs(:,:,1,n)=prsh(:,:,1)
 
   uu(:,:,nz,n)=uuh(:,:,nuvz)
@@ -307 +301 @@
   pv(:,:,nz,n)=pvh(:,:,nuvz)
   rho(:,:,nz,n)=rhoh(:,:,nuvz)
-! RLT
-  prs(:,:,nz,n)=prsh(:,:,nuvz)
+
 
   kmin=2
@@ -328 +321 @@
           pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
           rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
-! RLT
-          prs(ix,jy,iz,n)=prs(ix,jy,nz,n)
         else
           innuvz: do kz=idx(ix,jy),nuvz
@@ -363 +354 @@
           pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
           rho(ix,jy,iz,n)=(rhoh(ix,jy,kz-1)*dz2+rhoh(ix,jy,kz)*dz1)/dz
-! RLT add pressure
-          prs(ix,jy,iz,n)=(prsh(ix,jy,kz-1)*dz2+prsh(ix,jy,kz)*dz1)/dz
         endif
       enddo
@@ -665 +654 @@
         if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
 
-          do kz=nz,2,-1 !go Bottom up!
+          do kz=nz,1,-1 !go Bottom up!
             if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
               cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1) 

src/verttransform_gfs.f90

@@ -548 +548 @@
         if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
 
-          do kz=nz,2,-1 !go Bottom up!
+          do kz=nz,1,-1 !go Bottom up!
             if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
               cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
               clouds(ix,jy,kz,n)=1                               ! is a cloud
               if (lsp.ge.convp) then
-                clouds(ix,jy,kz,n)=3                             ! lsp in-cloud
+                clouds(ix,jy,kz,n)=3                            ! lsp in-cloud
               else
                 clouds(ix,jy,kz,n)=2                             ! convp in-cloud