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

create_tarball.sh

 #!/bin/bash
 
 #define version number
-version=10.3beta5
+
+githash=$(git rev-parse --short --verify HEAD)
+
+
+version=10.3beta5_$githash
 
 # define tarball name
 targetdir=../flexpart_distribution/
 tarball_tmp=${targetdir}flexpart_v$version 
+#tarball=${targetdir}flexpart_v$version.tar 
+tarball=${tarball_tmp}.tar 
 
 # clean old package
 if [ -d $tarball_tmp ]; then
-  echo $tarball_tmp exists: move to $tarball_tmp.bk   
-  mv $tarball_tmp ${tarball_tmp}.bk 
+  echo $tarball_tmp exists: move to $tarball_tmp.bk and exit  
+  mkdir $tarball_tmp.bk 
+  mv $tarball_tmp ${tarball_tmp}.bk/ 
+  mv $tarball ${tarball_tmp}.bk/ 
   exit 
 fi
 
@@ -79,6 +87,15 @@ mkdir $tarball_tmp/preprocess/flex_extract
 #cp -r flex_ecmwf_src/* $tarball_tmp/preprocess/flex_extract/
 ## cp -r flex_extract/work/EA* $tarball_tmp/preprocess/flex_extract/work   
 
+echo include flex_extract v7.0.4 b7c1c04a204c91e53759ef590504bf52dfaece64
+flex_extract=../flex_extract_v7.0.4/
+cp $flex_extract/README.md $tarball_tmp/preprocess/flex_extract
+cp -r $flex_extract/docs $tarball_tmp/preprocess/flex_extract
+cp -r $flex_extract/grib_templates $tarball_tmp/preprocess/flex_extract
+cp -r $flex_extract/python $tarball_tmp/preprocess/flex_extract
+cp -r $flex_extract/src $tarball_tmp/preprocess/flex_extract
+
+
 
 
 
@@ -174,7 +191,10 @@ mkdir $tarball_tmp/tests/ctbto
 # cp -r tests/NILU/test_1 $tarball_tmp/tests/
 # cp -r tests/default_cases $tarball_tmp/tests/
 
+tar cvf  $tarball  $tarball_tmp  
 
+echo  $tarball complete
+echo exported untarred files in $tarball_tmp 
 exit
 #return
 ###############################################################

src/FLEXPART.f90

@@ -136,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)
@@ -451,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 
@@ -466,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

@@ -145,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)

src/com_mod.f90

@@ -123,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
@@ -175,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)
@@ -360,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)
@@ -382,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

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)
+        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)
 
 
   ! Reinitialization of grid

src/concoutput_surf.f90

@@ -71,6 +71,9 @@ subroutine concoutput_surf(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,
@@ -100,6 +103,7 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
   logical :: sp_zer
   character :: adate*8,atime*6
   character(len=3) :: anspec
+  logical :: lexist
 
 
   if (verbosity.eq.1) then
@@ -179,6 +183,9 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
         jjy=max(min(nint(yl),nymin1),0)
         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
   end do
@@ -189,8 +196,14 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
     iix=max(min(nint(xl),nxmin1),0)
     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
   do kz=1,numzgrid
@@ -595,6 +608,49 @@ subroutine concoutput_surf(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,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/ &
        wetgridtotal
@@ -621,7 +677,23 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
     end do