Changes in / [4c52d0b:aa939a9] in flexpart.git

Files:

• .gitignore (modified) (2 diffs)
• src/calcfluxes.f90 (modified) (3 diffs)
• src/init_domainfill.f90 (modified) (2 diffs)
• src/netcdf_output_mod.f90 (modified) (3 diffs, 1 prop)
• src/readspecies.f90 (modified) (3 diffs)
• src/timemanager.f90 (modified) (2 diffs)
• src/timemanager_mpi.f90 (modified) (1 diff)

.gitignore

@@ -1 @@
-log*
-noh*
-diff*
 FP_ecmwf_gfortran*
 *.o
@@ -10 +7 @@
 FLEXPART
 FLEXPART_MPI
-FLEX*

src/calcfluxes.f90

@@ -64 +64 @@
        (jyave.le.numygrid-1)) then
     do kz=1,numzgrid                ! determine height of cell
-!      if (outheighthalf(kz).gt.zold) goto 11
-      if (outheight(kz).gt.zold) goto 11 !sec, use upper layer instead of mid layer
+      if (outheighthalf(kz).gt.zold) goto 11
     end do
 11   k1=min(numzgrid,kz)
     do kz=1,numzgrid                ! determine height of cell
-!      if (outheighthalf(kz).gt.ztra1(jpart)) goto 21
-      if (outheight(kz).gt.ztra1(jpart)) goto 21
+      if (outheighthalf(kz).gt.ztra1(jpart)) goto 21
     end do
 21   k2=min(numzgrid,kz)
@@ -98 +96 @@
 
     if (abs(xold-xtra1(jpart)).lt.real(nx)/2.) then
-      ix1=int((xold*dx+xoutshift)/dxout)              ! flux throught the western boundary (sec)
-      ix2=int((xtra1(jpart)*dx+xoutshift)/dxout)      ! flux throught the western boundary (sec)
-!      ix1=int((xold*dx+xoutshift)/dxout+0.5) 
-!      ix2=int((xtra1(jpart)*dx+xoutshift)/dxout+0.5)
+      ix1=int((xold*dx+xoutshift)/dxout+0.5)
+      ix2=int((xtra1(jpart)*dx+xoutshift)/dxout+0.5)
       do k=1,nspec
         do ix=ix1,ix2-1
@@ -149 +145 @@
   if ((kzave.le.numzgrid).and.(ixave.ge.0).and. &
        (ixave.le.numxgrid-1)) then
-    jy1=int((yold*dy+youtshift)/dyout)         ! flux throught the southern boundary (sec)
-    jy2=int((ytra1(jpart)*dy+youtshift)/dyout) ! flux throught the southern boundary (sec)
-!    jy1=int((yold*dy+youtshift)/dyout+0.5)
-!    jy2=int((ytra1(jpart)*dy+youtshift)/dyout+0.5)
+    jy1=int((yold*dy+youtshift)/dyout+0.5)
+    jy2=int((ytra1(jpart)*dy+youtshift)/dyout+0.5)
 
     do k=1,nspec

src/init_domainfill.f90

@@ -209 +209 @@
                    ztra1(numpart+jj)=height(nz)-0.5
 
+
 ! Interpolate PV to the particle position
 !****************************************
@@ -215 +216 @@
               ixp=ixm+1
               jyp=jym+1
-              if (jyp.gt.180) then
-                 write (*,*) 'init_domainfill, over: ',jyp,jym,ytra1(numpart+jj),jy,ran1(idummy),ny
-                 jyp=jym
-              endif
               ddx=xtra1(numpart+jj)-real(ixm)
               ddy=ytra1(numpart+jj)-real(jym)

src/netcdf_output_mod.f90

@@ -36 +36 @@
   use outg_mod,  only: outheight,oroout,densityoutgrid,factor3d,volume,&
                        wetgrid,wetgridsigma,drygrid,drygridsigma,grid,gridsigma,&
-                       area,arean,volumen, orooutn, areaeast, areanorth, p0out, t0out
+                       area,arean,volumen, orooutn, p0out, t0out
   use par_mod,   only: dep_prec, sp, dp, maxspec, maxreceptor, nclassunc,&
                        unitoutrecept,unitoutreceptppt, nxmax,unittmp, &
                        write_p0t0
-  use com_mod,   only: path,length,ldirect,bdate,ibdate,ibtime,iedate,ietime, &
+  use com_mod,   only: path,length,ldirect,ibdate,ibtime,iedate,ietime, &
                        loutstep,loutaver,loutsample,outlon0,outlat0,&
                        numxgrid,numygrid,dxout,dyout,numzgrid, height, &
@@ -67 +67 @@
 
   public :: writeheader_netcdf,concoutput_surf_nest_netcdf,concoutput_netcdf,&
-       &concoutput_nest_netcdf,concoutput_surf_netcdf,fluxoutput_netcdf
+       &concoutput_nest_netcdf,concoutput_surf_netcdf
 
 !  include 'netcdf.inc'
@@ -1520 +1520 @@
 end subroutine concoutput_surf_nest_netcdf
 
-subroutine fluxoutput_netcdf(itime)
-     
-  !                            i
-  !*****************************************************************************
-  !                                                                            *
-  !     Output of the gridded fluxes.                                          *
-  !     Eastward, westward, northward, southward, upward and downward gross    *
-  !     fluxes are written to output file in either sparse matrix or grid dump *
-  !     format, whichever is more efficient.                                   *
-  !                                                                            *
-  !     Author: A. Stohl                                                       *
-  !                                                                            *
-  !     04 April 2000                                                          *
-  !     netcdfoutput S. Eckhardt, 2020                                         *
-  !                                                                            *
-  !*****************************************************************************
-
-   use flux_mod
-
-  implicit none
-
-  real(kind=dp) :: jul
-  integer :: itime,ix,jy,kz,ks,nage,jjjjmmdd,ihmmss,kp,i
-  
-  character(len=255) :: ncfname
-  character :: adate*8,atime*6,timeunit*32,anspec*3
-
-  integer :: ncid 
-  integer :: timeDimID, latDimID, lonDimID, levDimID
-  integer :: nspecDimID, npointDimID, nageclassDimID, ncharDimID, pointspecDimID
-  integer :: tID, lonID, latID, levID, lageID, fluxID
-  integer, dimension(6)       :: dIDs
-  integer :: cache_size
-  real, allocatable, dimension(:) :: coord
-
-
-  ! Determine current calendar date, needed for the file name
-  !**********************************************************
-
-  jul=bdate+real(itime,kind=dp)/86400._dp
-  call caldate(jul,jjjjmmdd,ihmmss)
-  write(adate,'(i8.8)') jjjjmmdd
-  write(atime,'(i6.6)') ihmmss
-
-  ncfname=path(2)(1:length(2))//'grid_flux_'//adate// &
-       atime//'.nc'
-
-  ! setting cache size in bytes. It is set to 4 times the largest data block that is written
-  !   size_type x nx x ny x nz
-  ! create file
-
-  cache_size = 16 * numxgrid * numygrid * numzgrid
-  call nf90_err(nf90_create(trim(ncfname), cmode = nf90_hdf5, ncid = ncid, &
-    cache_size = cache_size))  
-
-  ! create dimensions:
-  !*************************
-  ! time
-  call nf90_err(nf90_def_dim(ncid, 'time', nf90_unlimited, timeDimID))
-  timeunit = 'seconds since '//adate(1:4)//'-'//adate(5:6)// &
-     '-'//adate(7:8)//' '//atime(1:2)//':'//atime(3:4)
-
-  call nf90_err(nf90_def_dim(ncid, 'longitude', numxgrid, lonDimID))
-  call nf90_err(nf90_def_dim(ncid, 'latitude', numygrid, latDimID))
-  call nf90_err(nf90_def_dim(ncid, 'height', numzgrid, levDimID))
-  call nf90_err(nf90_def_dim(ncid, 'numspec', nspec, nspecDimID))
-  call nf90_err(nf90_def_dim(ncid, 'pointspec', maxpointspec_act, pointspecDimID))
-  call nf90_err(nf90_def_dim(ncid, 'nageclass', nageclass, nageclassDimID))
-  call nf90_err(nf90_def_dim(ncid, 'nchar', 45, ncharDimID))
-  call nf90_err(nf90_def_dim(ncid, 'numpoint', numpoint, npointDimID))
-
-
-  ! create variables
-  !*************************
-
-  ! time
-  call nf90_err(nf90_def_var(ncid, 'time', nf90_int, (/ timeDimID /), tID))
-  call nf90_err(nf90_put_att(ncid, tID, 'units', timeunit))
-  call nf90_err(nf90_put_att(ncid, tID, 'calendar', 'proleptic_gregorian'))
-  timeID = tID
-
-  ! lon
-  call nf90_err(nf90_def_var(ncid, 'longitude', nf90_float, (/ lonDimID /), lonID))
-  call nf90_err(nf90_put_att(ncid, lonID, 'long_name', 'longitude in degree east'))
-  call nf90_err(nf90_put_att(ncid, lonID, 'axis', 'Lon'))
-  call nf90_err(nf90_put_att(ncid, lonID, 'units', 'degrees_east'))
-  call nf90_err(nf90_put_att(ncid, lonID, 'standard_name', 'grid_longitude'))
-  call nf90_err(nf90_put_att(ncid, lonID, 'description', 'grid cell centers'))
-
-  ! lat
-  call nf90_err(nf90_def_var(ncid, 'latitude', nf90_float, (/ latDimID /), latID))
-  call nf90_err(nf90_put_att(ncid, latID, 'long_name', 'latitude in degree north'))
-  call nf90_err(nf90_put_att(ncid, latID, 'axis', 'Lat'))
-  call nf90_err(nf90_put_att(ncid, latID, 'units', 'degrees_north'))
-  call nf90_err(nf90_put_att(ncid, latID, 'standard_name', 'grid_latitude'))
-  call nf90_err(nf90_put_att(ncid, latID, 'description', 'grid cell centers'))
-
-  ! height
-  call nf90_err(nf90_def_var(ncid, 'height', nf90_float, (/ levDimID /), levID))
-  call nf90_err(nf90_put_att(ncid, levID, 'units', 'meters'))
-  call nf90_err(nf90_put_att(ncid, levID, 'positive', 'up'))
-  call nf90_err(nf90_put_att(ncid, levID, 'standard_name', 'height'))
-
-  if (.not.allocated(coord)) allocate(coord(numxgrid)) 
-  do i = 1,numxgrid
-     coord(i) = outlon0 + (i-0.5)*dxout
-  enddo
-  call nf90_err(nf90_put_var(ncid, lonID, coord(1:numxgrid)))
-  deallocate(coord)
-  if (.not.allocated(coord)) allocate(coord(numygrid))
-  do i = 1,numygrid
-     coord(i) = outlat0 + (i-0.5)*dyout
-  enddo
-  call nf90_err(nf90_put_var(ncid, latID, coord(1:numygrid)))
-  deallocate(coord)
-  call nf90_err(nf90_put_var(ncid, levID, outheight(1:numzgrid)))
-
-  ! write time, one field per time - different to the others!
-  call nf90_err(nf90_put_var( ncid, timeID, itime, (/ 1 /)))
-  
-  dIDs = (/ londimid, latdimid, levdimid, timedimid, pointspecdimid, nageclassdimid /)
-
-  do ks=1,nspec
-    do kp=1,maxpointspec_act
-      do nage=1,nageclass
- 
-         write(anspec,'(i3.3)') ks
-
- ! East Flux
-        call nf90_err(nf90_def_var(ncid,'flux_east_'//anspec, nf90_float, dIDs, &
-             fluxID))
-
-        do jy=0,numygrid-1
-          do ix=0,numxgrid-1
-            do kz=1, numzgrid
-               grid(ix,jy,kz)=flux(1,ix,jy,kz,ks,kp,nage) 
-            end do
-          end do
-        end do
-
-        call nf90_err(nf90_put_var(ncid,fluxid,1.e12*grid(0:numxgrid-1,0:numygrid-1,1:numzgrid)&
-            /areaeast(0:numxgrid-1,0:numygrid-1,1:numzgrid)/loutstep,&
-           (/ 1,1,1,1,kp,nage /), (/ numxgrid,numygrid,numzgrid,1,1,1 /) ))
-
- ! West Flux
-        call nf90_err(nf90_def_var(ncid,'flux_west_'//anspec, nf90_float, dIDs, &
-             fluxID))
-
-        do jy=0,numygrid-1
-          do ix=0,numxgrid-1
-            do kz=1, numzgrid
-               grid(ix,jy,kz)=flux(2,ix,jy,kz,ks,kp,nage) 
-            end do
-          end do
-        end do
-
-        call nf90_err(nf90_put_var(ncid,fluxid,1.e12*grid(0:numxgrid-1,0:numygrid-1,1:numzgrid)&
-            /areaeast(0:numxgrid-1,0:numygrid-1,1:numzgrid)/loutstep,&
-           (/ 1,1,1,1,kp,nage /), (/ numxgrid,numygrid,numzgrid,1,1,1 /) ))
-
- ! North Flux
-        call nf90_err(nf90_def_var(ncid,'flux_north_'//anspec, nf90_float, dIDs, &
-             fluxID))
-
-        do jy=0,numygrid-1
-          do ix=0,numxgrid-1
-            do kz=1, numzgrid
-               grid(ix,jy,kz)=flux(4,ix,jy,kz,ks,kp,nage) 
-            end do
-          end do
-        end do
-
-        call nf90_err(nf90_put_var(ncid,fluxid,1.e12*grid(0:numxgrid-1,0:numygrid-1,1:numzgrid)&
-            /areanorth(0:numxgrid-1,0:numygrid-1,1:numzgrid)/loutstep,&
-           (/ 1,1,1,1,kp,nage /), (/ numxgrid,numygrid,numzgrid,1,1,1 /) ))
-
- ! South Flux
-        call nf90_err(nf90_def_var(ncid,'flux_south_'//anspec, nf90_float, dIDs, &
-             fluxID))
-
-        do jy=0,numygrid-1
-          do ix=0,numxgrid-1
-            do kz=1, numzgrid
-               grid(ix,jy,kz)=flux(3,ix,jy,kz,ks,kp,nage) 
-            end do
-          end do
-        end do
-
-        call nf90_err(nf90_put_var(ncid,fluxid,1.e12*grid(0:numxgrid-1,0:numygrid-1,1:numzgrid)&
-            /areanorth(0:numxgrid-1,0:numygrid-1,1:numzgrid)/loutstep,&
-           (/ 1,1,1,1,kp,nage /), (/ numxgrid,numygrid,numzgrid,1,1,1 /) ))
-
- ! Up Flux
-        call nf90_err(nf90_def_var(ncid,'flux_up_'//anspec, nf90_float, dIDs, &
-             fluxID))
-
-        do jy=0,numygrid-1
-          do ix=0,numxgrid-1
-            do kz=1, numzgrid
-               grid(ix,jy,kz)=flux(5,ix,jy,kz,ks,kp,nage)/area(ix,jy)
-            end do
-          end do
-        end do
-
-        call nf90_err(nf90_put_var(ncid,fluxid,1.e12*grid(0:numxgrid-1,0:numygrid-1,1:numzgrid)&
-            /loutstep,&
-           (/ 1,1,1,1,kp,nage /), (/ numxgrid,numygrid,numzgrid,1,1,1 /) ))
-
- ! Down Flux
-        call nf90_err(nf90_def_var(ncid,'flux_down_'//anspec, nf90_float, dIDs, &
-             fluxID))
-
-        do jy=0,numygrid-1
-          do ix=0,numxgrid-1
-            do kz=1, numzgrid
-               grid(ix,jy,kz)=flux(6,ix,jy,kz,ks,kp,nage)/area(ix,jy)
-            end do
-          end do
-        end do
-
-        call nf90_err(nf90_put_var(ncid,fluxid,1.e12*grid(0:numxgrid-1,0:numygrid-1,1:numzgrid)&
-            /loutstep,&
-           (/ 1,1,1,1,kp,nage /), (/ numxgrid,numygrid,numzgrid,1,1,1 /) ))
-
-      end do
-    end do
-  end do
-
-  ! Close netCDF file
-  call nf90_err(nf90_close(ncid))
-
-  ! Reinitialization of grid
-  !*************************
-
-  do ks=1,nspec
-    do kp=1,maxpointspec_act
-      do jy=0,numygrid-1
-        do ix=0,numxgrid-1
-          do kz=1,numzgrid
-            do nage=1,nageclass
-              do i=1,6
-                flux(i,ix,jy,kz,ks,kp,nage)=0.
-              end do
-            end do
-          end do
-        end do
-      end do
-    end do
-  end do
-
-
-end subroutine fluxoutput_netcdf
-
 end module netcdf_output_mod
 

src/readspecies.f90

@@ -42 +42 @@
   implicit none
 
-  integer :: i, pos_spec,j,icheck_dow_hour
+  integer :: i, pos_spec,j
   integer :: idow,ihour,id_spec
   character(len=3) :: aspecnumb
@@ -229 +229 @@
     ohnconst(pos_spec)=pohnconst
 
-
-    icheck_dow_hour=0
     do j=1,24     ! 24 hours, starting with 0-1 local time
       area_hour(pos_spec,j)=parea_hour(j)
       point_hour(pos_spec,j)=ppoint_hour(j)
-      if (parea_hour(j).ne.1 .or. ppoint_hour(j).ne.1) icheck_dow_hour=1
     end do
     do j=1,7      ! 7 days of the week, starting with Monday
       area_dow(pos_spec,j)=parea_dow(j)
       point_dow(pos_spec,j)=ppoint_dow(j)
-      if (parea_dow(j).ne.1 .or. ppoint_dow(j).ne.1) icheck_dow_hour=1
     end do
 
@@ -361 +357 @@
 20 continue
 
-  if ((icheck_dow_hour.eq.1).and.(ldirect.lt.0)) then
-    write(*,*) '#### FLEXPART MODEL WARNING                     ####'
-    write(*,*) '#### The variation for an emission release      ####'
-    write(*,*) '#### will have no effect in backward mode       ####'  
-  endif
 
 22 close(unitspecies)

src/timemanager.f90

@@ -85 +85 @@
 #ifdef USE_NCF
   use netcdf_output_mod, only: concoutput_netcdf,concoutput_nest_netcdf,&
-       &concoutput_surf_netcdf,concoutput_surf_nest_netcdf,fluxoutput_netcdf
+       &concoutput_surf_netcdf,concoutput_surf_nest_netcdf
 #endif
 
@@ -433 +433 @@
         endif
         if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
-        if ((iflux.eq.1).and.(lnetcdfout.eq.0)) call fluxoutput(itime)
-#ifdef USE_NCF
-        if ((iflux.eq.1).and.(lnetcdfout.eq.1)) call fluxoutput_netcdf(itime)
-#endif
+        if (iflux.eq.1) call fluxoutput(itime)
         write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
  

src/timemanager_mpi.f90

@@ -556 +556 @@
         endif
         if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
-        if ((iflux.eq.1).and.(lnetcdfout.eq.0)) call fluxoutput(itime)
-        if ((iflux.eq.1).and.(lnetcdfout.eq.1)) call fluxoutput_netcdf(itime)
+        if (iflux.eq.1) call fluxoutput(itime)
         if (mp_measure_time) call mpif_mtime('iotime',1)