Commit 4c52d0b9 authored by Espen Sollum's avatar Espen Sollum
Browse files

Merge branch 'dev' into GFS_025

parents aa939a93 4138764d
log*
noh*
diff*
FP_ecmwf_gfortran* FP_ecmwf_gfortran*
*.o *.o
*_mod.mod *_mod.mod
...@@ -6,3 +9,4 @@ class_gribfile.mod ...@@ -6,3 +9,4 @@ class_gribfile.mod
output output
FLEXPART FLEXPART
FLEXPART_MPI FLEXPART_MPI
FLEX*
...@@ -63,11 +63,13 @@ subroutine calcfluxes(nage,jpart,xold,yold,zold) ...@@ -63,11 +63,13 @@ subroutine calcfluxes(nage,jpart,xold,yold,zold)
if ((ixave.ge.0).and.(jyave.ge.0).and.(ixave.le.numxgrid-1).and. & if ((ixave.ge.0).and.(jyave.ge.0).and.(ixave.le.numxgrid-1).and. &
(jyave.le.numygrid-1)) then (jyave.le.numygrid-1)) then
do kz=1,numzgrid ! determine height of cell do kz=1,numzgrid ! determine height of cell
if (outheighthalf(kz).gt.zold) goto 11 ! if (outheighthalf(kz).gt.zold) goto 11
if (outheight(kz).gt.zold) goto 11 !sec, use upper layer instead of mid layer
end do end do
11 k1=min(numzgrid,kz) 11 k1=min(numzgrid,kz)
do kz=1,numzgrid ! determine height of cell do kz=1,numzgrid ! determine height of cell
if (outheighthalf(kz).gt.ztra1(jpart)) goto 21 ! if (outheighthalf(kz).gt.ztra1(jpart)) goto 21
if (outheight(kz).gt.ztra1(jpart)) goto 21
end do end do
21 k2=min(numzgrid,kz) 21 k2=min(numzgrid,kz)
...@@ -95,8 +97,10 @@ subroutine calcfluxes(nage,jpart,xold,yold,zold) ...@@ -95,8 +97,10 @@ subroutine calcfluxes(nage,jpart,xold,yold,zold)
! 1) Particle does not cross domain boundary ! 1) Particle does not cross domain boundary
if (abs(xold-xtra1(jpart)).lt.real(nx)/2.) then if (abs(xold-xtra1(jpart)).lt.real(nx)/2.) then
ix1=int((xold*dx+xoutshift)/dxout+0.5) ix1=int((xold*dx+xoutshift)/dxout) ! flux throught the western boundary (sec)
ix2=int((xtra1(jpart)*dx+xoutshift)/dxout+0.5) 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)
do k=1,nspec do k=1,nspec
do ix=ix1,ix2-1 do ix=ix1,ix2-1
if ((ix.ge.0).and.(ix.le.numxgrid-1)) then if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
...@@ -144,8 +148,10 @@ subroutine calcfluxes(nage,jpart,xold,yold,zold) ...@@ -144,8 +148,10 @@ subroutine calcfluxes(nage,jpart,xold,yold,zold)
if ((kzave.le.numzgrid).and.(ixave.ge.0).and. & if ((kzave.le.numzgrid).and.(ixave.ge.0).and. &
(ixave.le.numxgrid-1)) then (ixave.le.numxgrid-1)) then
jy1=int((yold*dy+youtshift)/dyout+0.5) jy1=int((yold*dy+youtshift)/dyout) ! flux throught the southern boundary (sec)
jy2=int((ytra1(jpart)*dy+youtshift)/dyout+0.5) 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)
do k=1,nspec do k=1,nspec
do jy=jy1,jy2-1 do jy=jy1,jy2-1
......
...@@ -208,13 +208,16 @@ subroutine init_domainfill ...@@ -208,13 +208,16 @@ subroutine init_domainfill
if (ztra1(numpart+jj).gt.height(nz)-0.5) & if (ztra1(numpart+jj).gt.height(nz)-0.5) &
ztra1(numpart+jj)=height(nz)-0.5 ztra1(numpart+jj)=height(nz)-0.5
! Interpolate PV to the particle position ! Interpolate PV to the particle position
!**************************************** !****************************************
ixm=int(xtra1(numpart+jj)) ixm=int(xtra1(numpart+jj))
jym=int(ytra1(numpart+jj)) jym=int(ytra1(numpart+jj))
ixp=ixm+1 ixp=ixm+1
jyp=jym+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) ddx=xtra1(numpart+jj)-real(ixm)
ddy=ytra1(numpart+jj)-real(jym) ddy=ytra1(numpart+jj)-real(jym)
rddx=1.-ddx rddx=1.-ddx
......
...@@ -35,11 +35,11 @@ module netcdf_output_mod ...@@ -35,11 +35,11 @@ module netcdf_output_mod
xpoint1,ypoint1,xpoint2,ypoint2,zpoint1,zpoint2,npart,xmass xpoint1,ypoint1,xpoint2,ypoint2,zpoint1,zpoint2,npart,xmass
use outg_mod, only: outheight,oroout,densityoutgrid,factor3d,volume,& use outg_mod, only: outheight,oroout,densityoutgrid,factor3d,volume,&
wetgrid,wetgridsigma,drygrid,drygridsigma,grid,gridsigma,& wetgrid,wetgridsigma,drygrid,drygridsigma,grid,gridsigma,&
area,arean,volumen, orooutn, p0out, t0out area,arean,volumen, orooutn, areaeast, areanorth, p0out, t0out
use par_mod, only: dep_prec, sp, dp, maxspec, maxreceptor, nclassunc,& use par_mod, only: dep_prec, sp, dp, maxspec, maxreceptor, nclassunc,&
unitoutrecept,unitoutreceptppt, nxmax,unittmp, & unitoutrecept,unitoutreceptppt, nxmax,unittmp, &
write_p0t0 write_p0t0
use com_mod, only: path,length,ldirect,ibdate,ibtime,iedate,ietime, & use com_mod, only: path,length,ldirect,bdate,ibdate,ibtime,iedate,ietime, &
loutstep,loutaver,loutsample,outlon0,outlat0,& loutstep,loutaver,loutsample,outlon0,outlat0,&
numxgrid,numygrid,dxout,dyout,numzgrid, height, & numxgrid,numygrid,dxout,dyout,numzgrid, height, &
outlon0n,outlat0n,dxoutn,dyoutn,numxgridn,numygridn, & outlon0n,outlat0n,dxoutn,dyoutn,numxgridn,numygridn, &
...@@ -66,7 +66,7 @@ module netcdf_output_mod ...@@ -66,7 +66,7 @@ module netcdf_output_mod
private private
public :: writeheader_netcdf,concoutput_surf_nest_netcdf,concoutput_netcdf,& public :: writeheader_netcdf,concoutput_surf_nest_netcdf,concoutput_netcdf,&
&concoutput_nest_netcdf,concoutput_surf_netcdf &concoutput_nest_netcdf,concoutput_surf_netcdf,fluxoutput_netcdf
! include 'netcdf.inc' ! include 'netcdf.inc'
...@@ -1519,6 +1519,259 @@ subroutine concoutput_surf_nest_netcdf(itime,outnum) ...@@ -1519,6 +1519,259 @@ subroutine concoutput_surf_nest_netcdf(itime,outnum)
end subroutine concoutput_surf_nest_netcdf 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 end module netcdf_output_mod
...@@ -41,7 +41,7 @@ subroutine readspecies(id_spec,pos_spec) ...@@ -41,7 +41,7 @@ subroutine readspecies(id_spec,pos_spec)
implicit none implicit none
integer :: i, pos_spec,j integer :: i, pos_spec,j,icheck_dow_hour
integer :: idow,ihour,id_spec integer :: idow,ihour,id_spec
character(len=3) :: aspecnumb character(len=3) :: aspecnumb
logical :: spec_found logical :: spec_found
...@@ -228,13 +228,17 @@ subroutine readspecies(id_spec,pos_spec) ...@@ -228,13 +228,17 @@ subroutine readspecies(id_spec,pos_spec)
ohdconst(pos_spec)=pohdconst ohdconst(pos_spec)=pohdconst
ohnconst(pos_spec)=pohnconst ohnconst(pos_spec)=pohnconst
icheck_dow_hour=0
do j=1,24 ! 24 hours, starting with 0-1 local time do j=1,24 ! 24 hours, starting with 0-1 local time
area_hour(pos_spec,j)=parea_hour(j) area_hour(pos_spec,j)=parea_hour(j)
point_hour(pos_spec,j)=ppoint_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 end do
do j=1,7 ! 7 days of the week, starting with Monday do j=1,7 ! 7 days of the week, starting with Monday
area_dow(pos_spec,j)=parea_dow(j) area_dow(pos_spec,j)=parea_dow(j)
point_dow(pos_spec,j)=ppoint_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 end do
endif endif
...@@ -356,6 +360,11 @@ subroutine readspecies(id_spec,pos_spec) ...@@ -356,6 +360,11 @@ subroutine readspecies(id_spec,pos_spec)
endif endif
20 continue 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) 22 close(unitspecies)
......
...@@ -84,7 +84,7 @@ subroutine timemanager(metdata_format) ...@@ -84,7 +84,7 @@ subroutine timemanager(metdata_format)
use com_mod use com_mod
#ifdef USE_NCF #ifdef USE_NCF
use netcdf_output_mod, only: concoutput_netcdf,concoutput_nest_netcdf,& use netcdf_output_mod, only: concoutput_netcdf,concoutput_nest_netcdf,&
&concoutput_surf_netcdf,concoutput_surf_nest_netcdf &concoutput_surf_netcdf,concoutput_surf_nest_netcdf,fluxoutput_netcdf
#endif #endif
implicit none implicit none
...@@ -432,7 +432,10 @@ subroutine timemanager(metdata_format) ...@@ -432,7 +432,10 @@ subroutine timemanager(metdata_format)
outnum=0. outnum=0.
endif endif
if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime) if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
if (iflux.eq.1) call fluxoutput(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
write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
!CGZ-lifetime: output species lifetime !CGZ-lifetime: output species lifetime
......
...@@ -555,7 +555,8 @@ subroutine timemanager(metdata_format) ...@@ -555,7 +555,8 @@ subroutine timemanager(metdata_format)
outnum=0. outnum=0.
endif endif
if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime) if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
if (iflux.eq.1) call fluxoutput(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 (mp_measure_time) call mpif_mtime('iotime',1) if (mp_measure_time) call mpif_mtime('iotime',1)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment