Skip to content
Snippets Groups Projects
Commit c5409402 authored by Lucie Bakels's avatar Lucie Bakels
Browse files

added Pressure to netcdf output, added option for restarting from netcdf partoutput file

parent 98156287
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/netcdf_output_mod.f90
+ 262
− 20
View file @ c5409402
......@@ -66,17 +66,18 @@ module netcdf_output_mod
public :: writeheader_netcdf,concoutput_surf_nest_netcdf,concoutput_netcdf,&
&concoutput_nest_netcdf,concoutput_surf_netcdf,writeheader_partoutput,partoutput_netcdf,&
open_partoutput_file,close_partoutput_file
open_partoutput_file,close_partoutput_file,create_particles_initialoutput,&
write_particles_initialoutput,partinitpointer,readpartpositions_netcdf
! include 'netcdf.inc'
! parameter for data compression (1-9, 9 = most aggressive)
integer, parameter :: deflate_level = 9
integer, parameter :: deflate_level = 5
logical, parameter :: min_size = .false. ! if set true, redundant fields (topography) are not written to minimize file size
character(len=255), parameter :: institution = 'NILU'
integer :: tpointer=0,tpointer_part=0,ppointer_part=0
character(len=255) :: ncfname, ncfnamen, ncfname_part!(maxpoint)
integer :: tpointer=0,tpointer_part=0,ppointer_part=0,partinitpointer=0
character(len=255) :: ncfname, ncfnamen, ncfname_part, ncfname_partinit!(maxpoint)
! netcdf dimension and variable IDs for main and nested output grid
integer, dimension(maxspec) :: specID,specIDppt, wdspecID,ddspecID
......@@ -86,9 +87,10 @@ module netcdf_output_mod
integer, dimension(5) :: depdimids, depdimidsn
!IDs for partoutput
integer :: partID
integer :: itramemID,topoID,pvID,qvID,rhoID
integer :: partID,filenumber=1
integer :: itramemID,topoID,pvID,qvID,rhoID,prID
integer :: hmixID,trID,ttID,lonID,latID,levID,massID(maxspec)
integer :: partIDi,tIDi,lonIDi,latIDi,levIDi ! For initial particle outputs
real,parameter :: eps=nxmax/3.e5
......@@ -1429,6 +1431,134 @@ subroutine concoutput_surf_nest_netcdf(itime,outnum)
print*,'Netcdf output for surface only not yet implemented'
end subroutine concoutput_surf_nest_netcdf
subroutine create_particles_initialoutput(itime,idate,itime_start,idate_start)
implicit none
integer, intent(in) :: itime,idate,itime_start,idate_start
! integer, intent(in) :: irelease
integer :: cache_size,ncid,j,totpart
integer :: partDimID
character(len=11) :: fprefix
character(len=3) :: anspec,arelease
character :: adate*8,atime*6,adate_start*8,atime_start*6,timeunit*32
character(len=255) :: fname_partoutput
real :: fillval
write(adate,'(i8.8)') idate
write(atime,'(i6.6)') itime
write(adate_start,'(i8.8)') idate_start
write(atime_start,'(i6.6)') itime_start
! write(arelease, '(i3.3)') irelease
fprefix = 'partinit_'!rel'//arelease//'_'
fname_partoutput = path(2)(1:length(2))//trim(fprefix)//adate//atime//'.nc'
!ncfname_part(irelease) = fname_partoutput
ncfname_partinit = fname_partoutput
call nf90_err(nf90_create(trim(fname_partoutput), cmode = nf90_hdf5, ncid = ncid))!, &
! cache_size = cache_size))
! create dimensions:
!*************************
! particle
partinitpointer=0
call nf90_err(nf90_def_dim(ncid, 'particle', nf90_unlimited, partDimID))
! create variables
!*************************
! particles
call nf90_err(nf90_def_var(ncid, 'particle', nf90_int, (/ partDimID/), partIDi))
call nf90_err(nf90_put_att(ncid, partIDi, 'long_name', 'particle index'))
fillval = -1.
! time
timeunit = 'seconds since '//adate_start(1:4)//'-'//adate_start(5:6)// &
'-'//adate_start(7:8)//' '//atime_start(1:2)//':'//atime_start(3:4)
call nf90_err(nf90_def_var(ncid, 'time', nf90_int, (/ partDimID /), tIDi))
call nf90_err(nf90_def_var_deflate(ncid,tIDi,shuffle=0,deflate=1,deflate_level=1))
call nf90_err(nf90_put_att(ncid, tIDi, 'long_name', 'time of release'))
call nf90_err(nf90_put_att(ncid, tIDi, '_FillValue', -1))
call nf90_err(nf90_put_att(ncid, tIDi, 'axis', 't'))
call nf90_err(nf90_put_att(ncid, tIDi, 'units', timeunit))
call nf90_err(nf90_put_att(ncid, tIDi, 'calendar', 'proleptic_gregorian'))
call nf90_err(nf90_put_att(ncid, tIDi, 'standard_name', 'time'))
call nf90_err(nf90_put_att(ncid, tIDi, 'description', 'time of release'))
! lon
call nf90_err(nf90_def_var(ncid, 'longitude', nf90_float, (/ partDimID /), lonIDi))
call nf90_err(nf90_def_var_deflate(ncid,lonIDi,shuffle=0,deflate=1,deflate_level=1))
call nf90_err(nf90_put_att(ncid, lonIDi, 'long_name', 'longitude in degree east'))
call nf90_err(nf90_put_att(ncid, lonIDi, '_FillValue', fillval))
call nf90_err(nf90_put_att(ncid, lonIDi, 'axis', 'Lon'))
call nf90_err(nf90_put_att(ncid, lonIDi, 'units', 'degrees_east'))
call nf90_err(nf90_put_att(ncid, lonIDi, 'standard_name', 'longitude'))
call nf90_err(nf90_put_att(ncid, lonIDi, 'description', 'longitude of particles'))
! lat
call nf90_err(nf90_def_var(ncid, 'latitude', nf90_float, (/ partDimID /), latIDi))
call nf90_err(nf90_def_var_deflate(ncid,latIDi,shuffle=0,deflate=1,deflate_level=1))
call nf90_err(nf90_put_att(ncid, latIDi, 'long_name', 'latitude in degree north'))
call nf90_err(nf90_put_att(ncid, latIDi, 'axis', 'Lat'))
call nf90_err(nf90_put_att(ncid, latIDi, '_FillValue', fillval))
call nf90_err(nf90_put_att(ncid, latIDi, 'units', 'degrees_north'))
call nf90_err(nf90_put_att(ncid, latIDi, 'standard_name', 'latitude'))
call nf90_err(nf90_put_att(ncid, latIDi, 'description', 'latitude of particles'))
! height
call nf90_err(nf90_def_var(ncid, 'height', nf90_float, (/ partDimID /), levIDi))
call nf90_err(nf90_def_var_deflate(ncid,levIDi,shuffle=0,deflate=1,deflate_level=1))
call nf90_err(nf90_put_att(ncid, levIDi, 'units', 'meters'))
call nf90_err(nf90_put_att(ncid, levIDi, '_FillValue', fillval))
call nf90_err(nf90_put_att(ncid, levIDi, 'positive', 'up'))
call nf90_err(nf90_put_att(ncid, levIDi, 'standard_name', 'height'))
call nf90_err(nf90_put_att(ncid, levIDi, 'long_name', 'height above ground'))
! moves the file from define to data mode
call nf90_err(nf90_enddef(ncid))
call nf90_err(nf90_close(ncid))
end subroutine create_particles_initialoutput
subroutine write_particles_initialoutput(itime,istart,iend)
use com_mod
implicit none
integer, intent(in) :: &
itime, & ! time of particle release
istart, & ! index of first newly released particle
iend ! index of last newly released partile
integer, allocatable :: partindices(:),releasetimes(:)
integer :: newpart,ncid,j
newpart = iend-istart
call nf90_err(nf90_open(trim(ncfname_partinit), nf90_write, ncid))
allocate ( partindices(newpart) )
do j=1,newpart
partindices(j)=j+partinitpointer
end do
call nf90_err(nf90_put_var(ncid,partIDi,partindices,(/ partinitpointer+1 /),(/ newpart /)))
deallocate (partindices)
allocate ( releasetimes(newpart) )
releasetimes=itime
call nf90_err(nf90_put_var(ncid,tIDi,releasetimes,(/ partinitpointer+1 /),(/ newpart /)))
deallocate (releasetimes)
call nf90_err(nf90_put_var(ncid,lonIDi,xlon0+xtra1(partinitpointer+1:iend)*dx, (/ partinitpointer+1 /),(/ newpart /)))
call nf90_err(nf90_put_var(ncid,latIDi,ylat0+ytra1(partinitpointer+1:iend)*dy, (/ partinitpointer+1 /),(/ newpart /)))
call nf90_err(nf90_put_var(ncid,levIDi,ztra1(partinitpointer+1:iend), (/ partinitpointer+1 /),(/ newpart /)))
call nf90_err(nf90_close(ncid))
partinitpointer = partinitpointer+newpart
end subroutine write_particles_initialoutput
subroutine writeheader_partoutput(itime,idate,itime_start,idate_start)!,irelease)
implicit none
......@@ -1437,11 +1567,10 @@ subroutine writeheader_partoutput(itime,idate,itime_start,idate_start)!,irelease
! integer, intent(in) :: irelease
integer,allocatable :: partindices(:)
integer :: cache_size,ncid,j,totpart
!integer, dimension(:), allocatable :: partindices
integer :: timeDimID,partDimID,tID,memDimID
character(len=11) :: fprefix
character(len=3) :: anspec,arelease
character :: adate*8,atime*6,adate_start*8,atime_start*6,timeunit*32
character :: adate*8,atime*6,adate_start*8,atime_start*6,timeunit*32,afilenumber*8
character(len=255) :: fname_partoutput
real :: fillval
......@@ -1453,6 +1582,7 @@ subroutine writeheader_partoutput(itime,idate,itime_start,idate_start)!,irelease
write(atime,'(i6.6)') itime
write(adate_start,'(i8.8)') idate_start
write(atime_start,'(i6.6)') itime_start
write(afilenumber,'(i8)') filenumber
! write(arelease, '(i3.3)') irelease
fprefix = 'partoutput_'!rel'//arelease//'_'
......@@ -1542,6 +1672,15 @@ subroutine writeheader_partoutput(itime,idate,itime_start,idate_start)!,irelease
call nf90_err(nf90_put_att(ncid, topoID, 'standard_name', 'topography'))
call nf90_err(nf90_put_att(ncid, topoID, 'long_name', 'topography above sealevel'))
! Pressure
call nf90_err(nf90_def_var(ncid, 'pr', nf90_float, (/ timeDimID,partDimID /), prID))
call nf90_err(nf90_def_var_chunking(ncid,prID,NF90_CHUNKED,chunksizes=(/ 1,totpart /)))
call nf90_err(nf90_def_var_deflate(ncid,prID,shuffle=0,deflate=1,deflate_level=1))
call nf90_err(nf90_put_att(ncid, prID, 'units', 'Pa'))
call nf90_err(nf90_put_att(ncid, prID, '_FillValue', fillval))
call nf90_err(nf90_put_att(ncid, prID, 'standard_name', 'pressure'))
call nf90_err(nf90_put_att(ncid, prID, 'long_name', 'pressure'))
! Potential vorticity
call nf90_err(nf90_def_var(ncid, 'pv', nf90_float, (/ timeDimID,partDimID /), pvID))
call nf90_err(nf90_def_var_chunking(ncid,pvID,NF90_CHUNKED,chunksizes=(/ 1,totpart /)))
......@@ -1622,6 +1761,9 @@ subroutine writeheader_partoutput(itime,idate,itime_start,idate_start)!,irelease
call nf90_err(nf90_close(ncid))
! Update filenumber for next file
!********************************
filenumber = filenumber + 1
return
110 write(*,FMT='(80("#"))')
write(*,*) 'ERROR: output directory ', trim(path(2)(1:length(2))), ' does not exist&
......@@ -1694,33 +1836,133 @@ subroutine partoutput_netcdf(itime,field,fieldname,imass,ncid)
ppointer_part = numpart
endif
case('LO')
case('LO') ! Longitude
call nf90_err(nf90_put_var(ncid,lonID,field, (/ tpointer_part,1 /),(/ 1,numpart /)))
case('LA')
case('LA') ! Latitude
call nf90_err(nf90_put_var(ncid,latID,field, (/ tpointer_part,1 /),(/ 1,numpart /)))
case('ZZ')
case('ZZ') ! Height
call nf90_err(nf90_put_var(ncid,levID,field, (/ tpointer_part,1 /),(/ 1,numpart /)))
case('IT')
case('IT') ! Itramem (not in use atm)
call nf90_err(nf90_put_var(ncid,itramemID,field, (/ tpointer_part,1 /),(/ 1,numpart /)))
case('TO')
case('TO') ! Topography
call nf90_err(nf90_put_var(ncid,topoID,field, (/ tpointer_part,1 /),(/ 1,numpart /)))
case('PV')
case('PV') ! Potential vorticity
call nf90_err(nf90_put_var(ncid,pvID,field, (/ tpointer_part,1 /),(/ 1,numpart /)))
case('QV')
case('PR') ! Pressure
call nf90_err(nf90_put_var(ncid,prID,field, (/ tpointer_part,1 /),(/ 1,numpart /)))
case('QV') ! Specific humidity
call nf90_err(nf90_put_var(ncid,qvID,field, (/ tpointer_part,1 /),(/ 1,numpart /)))
case('RH')
case('RH') ! Air density
call nf90_err(nf90_put_var(ncid,rhoID,field, (/ tpointer_part,1 /),(/ 1,numpart /)))
case('HM')
case('HM') ! Mixing height
call nf90_err(nf90_put_var(ncid,hmixID,field, (/ tpointer_part,1 /),(/ 1,numpart /)))
case('TR')
case('TR') ! Tropopause
call nf90_err(nf90_put_var(ncid,trID,field, (/ tpointer_part,1 /),(/ 1,numpart /)))
case('TT')
case('TT') ! Temperature
call nf90_err(nf90_put_var(ncid,ttID,field, (/ tpointer_part,1 /),(/ 1,numpart /)))
case('MA')
case('MA') ! Mass
call nf90_err(nf90_put_var(ncid,massID(imass),field, (/ tpointer_part,1 /),(/ 1,numpart /)))
end select
! call nf90_err(nf90_close(ncid))
end subroutine partoutput_netcdf
subroutine readpartpositions_netcdf(ibtime_l,ibdate_l)
use par_mod
use com_mod
use random_mod
implicit none
integer, intent(in) :: ibtime_l,ibdate_l
integer :: ncidend,tIDend,pIDend,tempIDend,j
integer :: tlen,plen,tend,i
integer :: idate_start,itime_start
character :: adate*8,atime*6,timeunit*32,adate_start*8,atime_start*6
real(kind=dp) :: julin,julcommand,julpartin,juldate
character(len=3) :: anspec
integer :: idummy = -8
write(adate,'(i8.8)') ibdate_l
write(atime,'(i6.6)') ibtime_l
if (mquasilag.ne.0) then
write(*,*) 'Combination of ipin, netcdf partoutput, and mquasilag!=0 does not work yet'
stop
endif
! Open partoutput_end.nc file
call nf90_err(nf90_open(trim('partoutput_end.nc'), mode=NF90_NOWRITE,ncid=ncidend))
! Take the positions of the particles at the last timestep in the file
! It needs to be the same as given in the COMMAND file, this is arbitrary
! and should be removed in the future for easier use
! First get the time dimension
call nf90_err(nf90_inq_dimid(ncid=ncidend,name='time',dimid=tIDend))
call nf90_err(nf90_inquire_dimension(ncid=ncidend,dimid=tIDend,len=tlen))
! Check if the time corresponds to the one given in the COMMAND file
call nf90_err(nf90_inq_varid(ncid=ncidend,name='time',varid=tIDend))
call nf90_err(nf90_get_att(ncid=ncidend,varid=tIDend,name='units',values=timeunit))
call nf90_err(nf90_get_var(ncid=ncidend,varid=tIDend,values=tend,start=(/ tlen /)))!,count=(/ 1 /)))
adate_start(1:4) = timeunit(15:18)
adate_start(5:6) = timeunit(20:21)
adate_start(7:8) = timeunit(23:24)
atime_start = '000000'
atime_start(1:2) = timeunit(26:27)
atime_start(3:4) = timeunit(29:30)
read(adate_start,*) idate_start
read(atime_start,*) itime_start
julin = juldate(idate_start,itime_start)+real(tend,kind=dp)/86400._dp
julcommand = juldate(ibdate_l,ibtime_l)
if (abs(julin-julcommand).gt.1.e-5) then
write(*,*) 'ERROR: The given starting time and date do not correspond to'
write(*,*) 'the last timestep of partoutput_end.nc:'
write(*,*) julin,julcommand,tend
stop
endif
! Then the particle dimension
call nf90_err(nf90_inq_dimid(ncid=ncidend,name='particle',dimid=pIDend))
call nf90_err(nf90_inquire_dimension(ncid=ncidend,dimid=pIDend,len=plen))
! Now spawn the correct number of particles
write(*,*) 'Npart:',plen
numpart=plen
! And give them the correct positions
! Longitude
call nf90_err(nf90_inq_varid(ncid=ncidend,name='longitude',varid=tempIDend))
call nf90_err(nf90_get_var(ncid=ncidend,varid=tempIDend,values=xtra1(1:numpart), &
start=(/ tlen, 1 /),count=(/ 1, plen /)))
xtra1(1:numpart)=(xtra1(1:numpart)-xlon0)/dx
! Latitude
call nf90_err(nf90_inq_varid(ncid=ncidend,name='latitude',varid=tempIDend))
call nf90_err(nf90_get_var(ncid=ncidend,varid=tempIDend,values=ytra1(1:numpart), &
start=(/ tlen, 1 /),count=(/ 1, plen /)))
ytra1(1:numpart)=(ytra1(1:numpart)-ylat0)/dy
! Height
call nf90_err(nf90_inq_varid(ncid=ncidend,name='height',varid=tempIDend))
call nf90_err(nf90_get_var(ncid=ncidend,varid=tempIDend,values=ztra1(1:numpart), &
start=(/ tlen, 1 /),count=(/ 1, plen /)))
! Mass
do j=1,nspec
write(anspec, '(i3.3)') j
call nf90_err(nf90_inq_varid(ncid=ncidend,name='mass'//anspec,varid=tempIDend))
call nf90_err(nf90_get_var(ncid=ncidend,varid=tempIDend,values=xmass1(1:numpart,j), &
start=(/ tlen, 1 /),count=(/ 1, plen /)))
end do
do i=1,numpart
nclass(i)=min(int(ran1(idummy)*real(nclassunc))+1, &
nclassunc)
idt(i)=mintime
npoint(i)=i
itramem(i)=0
end do
write(*,*) 'TEST: ',xtra1(1),ytra1(1),ztra1(1),xmass1(1,1),nclass(1),idt(1),npoint(1)
call nf90_err(nf90_close(ncidend))
end subroutine readpartpositions_netcdf
end module netcdf_output_mod
src/partoutput.f90
+ 19
− 3
View file @ c5409402
......@@ -31,13 +31,13 @@ subroutine partoutput(itime)
integer :: itime,i,j,jjjjmmdd,ihmmss
integer :: ix,jy,ixp,jyp,indexh,m,il,ind,indz,indzp
real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
real :: hm(2),pv1(2),pvprof(2),qv1(2),qvprof(2)
real :: hm(2),pv1(2),pvprof(2),qv1(2),qvprof(2),pr1(2),prprof(2)
real :: tt1(2),ttprof(2),rho1(2),rhoprof(2)
real :: tr(2)
character :: adate*8,atime*6
real :: xlon(numpart),ylat(numpart)
real :: tti(numpart),rhoi(numpart),pvi(numpart),qvi(numpart)
real :: tti(numpart),rhoi(numpart),pvi(numpart),qvi(numpart),pri(numpart)
real :: topo(numpart),hmixi(numpart),tri(numpart)
#ifdef USE_NCF
......@@ -65,6 +65,8 @@ subroutine partoutput(itime)
topo(i)=-1.
hmixi(i)=-1.
tri(i)=-1.
pri(i)=-1.
if (itra1(i).eq.itime) then
xlon(i)=xlon0+xtra1(i)*dx
ylat(i)=ylat0+ytra1(i)*dy
......@@ -141,16 +143,23 @@ subroutine partoutput(itime)
+p2*rho(ixp,jy ,ind,indexh) &
+p3*rho(ix ,jyp,ind,indexh) &
+p4*rho(ixp,jyp,ind,indexh)
! Pressure
pr1(m)=p1*prs(ix ,jy ,ind,indexh) &
+p2*prs(ixp,jy ,ind,indexh) &
+p3*prs(ix ,jyp,ind,indexh) &
+p4*prs(ixp,jyp,ind,indexh)
end do
pvprof(ind-indz+1)=(pv1(1)*dt2+pv1(2)*dt1)*dtt
qvprof(ind-indz+1)=(qv1(1)*dt2+qv1(2)*dt1)*dtt
ttprof(ind-indz+1)=(tt1(1)*dt2+tt1(2)*dt1)*dtt
rhoprof(ind-indz+1)=(rho1(1)*dt2+rho1(2)*dt1)*dtt
prprof(ind-indz+1)=(pr1(1)*dt2+pr1(2)*dt1)*dtt
end do
pvi(i)=(dz1*pvprof(2)+dz2*pvprof(1))*dz
qvi(i)=(dz1*qvprof(2)+dz2*qvprof(1))*dz
tti(i)=(dz1*ttprof(2)+dz2*ttprof(1))*dz
rhoi(i)=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz
pri(i)=(dz1*prprof(2)+dz2*prprof(1))*dz
! Tropopause and PBL height
!**************************
......@@ -178,7 +187,13 @@ subroutine partoutput(itime)
end do
! Determine current calendar date, needed for the file name
!**********************************************************
if (numpart.gt.0) then
write(*,*) 'topo: ', topo(1), 'z:', ztra1(1)
write(*,*) 'xtra: ', xlon(1)
write(*,*) 'ytra: ', ylat(1)
write(*,*) 'mass: ', xmass1(1,1)
write(*,*) pvi(1),qvi(1),tti(1),rhoi(1)
endif
jul=bdate+real(itime,kind=dp)/86400._dp
call caldate(jul,jjjjmmdd,ihmmss)
write(adate,'(i8.8)') jjjjmmdd
......@@ -201,6 +216,7 @@ subroutine partoutput(itime)
call partoutput_netcdf(itime,hmixi,'HM',j,ncid)
call partoutput_netcdf(itime,tri,'TR',j,ncid)
call partoutput_netcdf(itime,tti,'TT',j,ncid)
call partoutput_netcdf(itime,pri,'PR',j,ncid)
do j=1,nspec
call partoutput_netcdf(itime,xmass1(:,j),'MA',j,ncid)
end do
......
......
src/readpartpositions.f90
+ 20
− 0
View file @ c5409402
......@@ -22,6 +22,7 @@ subroutine readpartpositions
use par_mod
use com_mod
use random_mod
use netcdf_output_mod
implicit none
......@@ -37,6 +38,25 @@ subroutine readpartpositions
! Open header file of dumped particle data
!*****************************************
! Open header file of dumped particle data
!*****************************************
if (lnetcdfout.eq.1) then
#ifdef USE_NCF
call readpartpositions_netcdf(ibtime,ibdate)
write(*,*) 'TEST2: ',xtra1(1),ytra1(1),ztra1(1),xmass1(1,1),nclass(1),idt(1),npoint(1)
do i=1,numpart
nclass(i)=min(int(ran1(idummy)*real(nclassunc))+1, &
nclassunc)
idt(i)=mintime
itra1(i)=0
itrasplit(i)=ldirect*itsplit
end do
numparticlecount=numpart
return
#endif
endif
open(unitpartin,file=path(2)(1:length(2))//'header', &
form='unformatted',err=998)
......
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment