subroutine initstatsdump
use modmpi, only : my_real,mpierr,comm3d,mpi_logical,mpi_integer,mpi_character,cmyid
use modglobal,only : imax,jmax,kmax,cexpnr,ifnamopt,fname_options,kb,ke,ladaptive,btime,&
nsv,lslicedump,lxytdump
use modstat_nc,only: open_nc, define_nc,ncinfo,writestat_dims_nc
use modfields,only : ncstaty,ncstatyt,ncstattke,ncstatxy,ncstatslice,ncstatxyt,ncstatt
implicit none
integer :: ierr
namelist/NAMSTATSDUMP/ &
lydump,tsample,klow,khigh,tstatsdump,lytdump,ltkedump,lxydump,lxytdump,ltdump
allocate(ncstaty(nstaty,4))
allocate(ncstatyt(nstatyt,4))
allocate(ncstattke(nstattke,4))
allocate(ncstatxy(nstatxy,4))
allocate(ncstatslice(nstatslice,4))
allocate(ncstatxyt(nstatxyt,4))
allocate(ncstatt(nstatt,4))
klow=kb
khigh=ke
if(myid==0)then
open(ifnamopt,file=fname_options,status='old',iostat=ierr)
read (ifnamopt,NAMSTATSDUMP,iostat=ierr)
if (ierr > 0) then
write(0, *) 'ERROR: Problem in namoptions NAMSTATSDUMP'
write(0, *) 'iostat error: ', ierr
stop 1
endif
write(6 ,NAMSTATSDUMP)
close(ifnamopt)
end if
call MPI_BCAST(klow ,1,MPI_INTEGER,0,comm3d,ierr) !have to do this? just want nc for first CPU
call MPI_BCAST(khigh ,1,MPI_INTEGER,0,comm3d,ierr)
call MPI_BCAST(nstatt ,1,MPI_INTEGER,0,comm3d,ierr)
! call MPI_BCAST(nstaty ,1,MPI_INTEGER,0,comm3d,ierr)
call MPI_BCAST(ncstatyt ,80,MPI_CHARACTER,0,comm3d,mpierr)
call MPI_BCAST(ncstaty ,80,MPI_CHARACTER,0,comm3d,mpierr)
call MPI_BCAST(ncstattke ,80,MPI_CHARACTER,0,comm3d,mpierr)
call MPI_BCAST(ncstatxy ,80,MPI_CHARACTER,0,comm3d,mpierr)
call MPI_BCAST(ncstatxyt ,80,MPI_CHARACTER,0,comm3d,mpierr)
call MPI_BCAST(ncstatt ,80,MPI_CHARACTER,0,comm3d,mpierr)
call MPI_BCAST(ltdump ,1,MPI_LOGICAL,0,comm3d,ierr)
!> Generate y-averaged NetCDF: ydump.xxx.nc
if(lydump) then
yname(7:9) = cexpnr
call ncinfo(tncstaty(1,:),'time' ,'Time' ,'s' ,'time')
call ncinfo(ncstaty( 1,:),'uy' ,'Streamwise velocity' ,'m/s' ,'m0tt')
call ncinfo(ncstaty( 2,:),'vy' ,'Spanwise velocity' ,'m/s' ,'t0tt')
call ncinfo(ncstaty( 3,:),'wy' ,'Vertical velocity' ,'m/s' ,'t0mt')
call ncinfo(ncstaty( 4,:),'thly' ,'Temperature' ,'K' ,'t0tt')
call ncinfo(ncstaty( 5,:),'qty' ,'Moisture' ,'kg/kg' ,'t0tt')
call ncinfo(ncstaty( 6,:),'sca1y' ,'Scalar field 1' ,'kg/m^3' ,'t0tt')
call ncinfo(ncstaty( 7,:),'sca2y' ,'Scalar field 2' ,'kg/m^3' ,'t0tt')
call ncinfo(ncstaty( 8,:),'sca3y' ,'Scalar field 3' ,'kg/m^3' ,'t0tt')
call ncinfo(ncstaty( 9,:),'upwpy' ,'Turbulent mom. flux' ,'m^2/s^2','m0mt')
call ncinfo(ncstaty(10,:),'wpthlpy' ,'Turbulent heat flux' ,'K m/s' ,'t0mt')
call ncinfo(ncstaty(11,:),'usgsy' ,'SGS mom. flux' ,'m^2/s^2','m0mt')
call ncinfo(ncstaty(12,:),'thlsgsy' ,'SGS heat flux' ,'K m/s' ,'t0mt')
call ncinfo(ncstaty(13,:),'uwyik' ,'Advective mom. flux' ,'m^2/s^2','m0mt')
call ncinfo(ncstaty(14,:),'wthlyk' ,'Advective heat flux' ,'K m/s' ,'t0mt')
if (myid==0) then
call open_nc(yname, ncidy, nrecy, n1=imax, n3=khigh-klow+1)
if (nrecy==0) then
call define_nc( ncidy, 1, tncstaty)
call writestat_dims_nc(ncidy)
end if
call define_nc( ncidy, nstaty, ncstaty)
endif !myid==0
endif
!> Generate time and y averaged NetCDF: ytdump.xxx.nc
if (lytdump) then
ytname(8:10) = cexpnr
call ncinfo(tncstatyt(1,:),'time' ,'Sampling time' ,'s' ,'time')
call ncinfo(ncstatyt( 1,:),'uyt' ,'Streamwise velocity' ,'m/s' ,'m0tt')
call ncinfo(ncstatyt( 2,:),'vyt' ,'Spanwise velocity' ,'m/s' ,'t0tt')
call ncinfo(ncstatyt( 3,:),'wyt' ,'Vertical velocity' ,'m/s' ,'t0mt')
call ncinfo(ncstatyt( 4,:),'thlyt' ,'Temperature' ,'K' ,'t0tt')
call ncinfo(ncstatyt( 5,:),'qtyt' ,'Moisture' ,'kg/kg' ,'t0tt')
call ncinfo(ncstatyt( 6,:),'sca1yt' ,'Scalar field 1' ,'kg/m^3' ,'t0tt')
call ncinfo(ncstatyt( 7,:),'sca2yt' ,'Scalar field 2' ,'kg/m^3' ,'t0tt')
call ncinfo(ncstatyt( 8,:),'sca3yt' ,'Scalar field 3' ,'kg/m^3' ,'t0tt')
call ncinfo(ncstatyt( 9,:),'upwpyt' ,'Turbulent mom. flux' ,'m^2/s^2' ,'m0mt')
call ncinfo(ncstatyt( 10,:),'wpthlpyt' ,'Turbulent heat flux' ,'K m/s' ,'t0mt')
call ncinfo(ncstatyt( 11,:),'wpqtpyt' ,'Turbulent moisture flux' ,'kg/kg m/s','t0mt')
call ncinfo(ncstatyt( 12,:),'wpsca1tpyt','Turbulent scalar flux' ,'M m/s' ,'t0mt')
call ncinfo(ncstatyt( 13,:),'wpsca2tpyt','Turbulent scalar flux' ,'M m/s' ,'t0mt')
call ncinfo(ncstatyt( 14,:),'wpsca3tpyt','Turbulent scalar flux' ,'M m/s' ,'t0mt')
call ncinfo(ncstatyt( 15,:),'uwyt' ,'Kinematic mom. flux' ,'m^2/s^2' ,'m0mt')
call ncinfo(ncstatyt( 16,:),'wthlyt' ,'Kinematic heat flux' ,'K m/s' ,'t0mt')
call ncinfo(ncstatyt( 17,:),'wqtyt' ,'Kinematic moisture flux' ,'K m/s' ,'t0mt')
call ncinfo(ncstatyt( 18,:),'wsca1yt' ,'Kinematic scalar flux' ,'K m/s' ,'t0mt')
call ncinfo(ncstatyt( 19,:),'wsca2yt' ,'Kinematic scalar flux' ,'K m/s' ,'t0mt')
call ncinfo(ncstatyt( 20,:),'wsca3yt' ,'Kinematic scalar flux' ,'K m/s' ,'t0mt')
call ncinfo(ncstatyt( 21,:),'upupyt' ,'mom. variance' ,'m^2/s^2' ,'m0tt')
call ncinfo(ncstatyt( 22,:),'wpwpyt' ,'mom. variance' ,'m^2/s^2' ,'t0mt')
call ncinfo(ncstatyt( 23,:),'thlpthlpyt' ,'temp. variance' ,'K^2' ,'t0tt')
call ncinfo(ncstatyt( 24,:),'qtpqtpyt' ,'moisture. variance' ,'kg^2/kg^2','t0tt')
call ncinfo(ncstatyt( 25,:),'sca1tpsca1pyt','scalar. variance' ,'M^2' ,'t0tt')
call ncinfo(ncstatyt( 26,:),'sca2tpsca2pyt','scalar. variance' ,'M^2' ,'t0tt')
call ncinfo(ncstatyt( 27,:),'sca3tpsca3pyt','scalar. variance' ,'M^2' ,'t0tt')
call ncinfo(ncstatyt( 28,:),'usgsyt' ,'SGS mom. flux' ,'m^2/s^2' ,'m0mt')
call ncinfo(ncstatyt( 29,:),'wsgsyt' ,'SGS mom. flux' ,'m^2/s^2' ,'t0mt')
call ncinfo(ncstatyt( 30,:),'thlsgsyt' ,'SGS heat flux' ,'K m/s' ,'t0mt')
call ncinfo(ncstatyt( 31,:),'qtsgsyt' ,'SGS moisture flux' ,'kg/kg m/s','t0mt')
call ncinfo(ncstatyt( 32,:),'sca1sgsyt' ,'SGS scalar flux' ,'M m/s' ,'t0mt')
call ncinfo(ncstatyt( 33,:),'sca2sgsyt' ,'SGS scalar flux' ,'M m/s' ,'t0mt')
call ncinfo(ncstatyt( 34,:),'sca3sgsyt' ,'SGS scalar flux' ,'M m/s' ,'t0mt')
if (myid==0) then
call open_nc(ytname, ncidyt, nrecyt, n1=imax, n3=khigh-klow+1)
if (nrecyt==0) then
call define_nc( ncidyt, 1, tncstatyt)
call writestat_dims_nc(ncidyt)
end if
call define_nc( ncidyt, nstatyt, ncstatyt)
endif !myid==0
endif
!> Generate y and x averaged NetCDF: xydump.xxx.nc
if (lxydump) then
xyname(8:10) = cexpnr
call ncinfo(tncstatxy(1,:),'time' ,'Time' ,'s' ,'time')
call ncinfo(ncstatxy( 1,:),'uxy' ,'Streamwise velocity' ,'m/s' ,'tt' )
call ncinfo(ncstatxy( 2,:),'vxy' ,'Spanwise velocity' ,'m/s' ,'tt' )
call ncinfo(ncstatxy( 3,:),'wxy' ,'Vertical velocity' ,'m/s' ,'mt' )
call ncinfo(ncstatxy( 4,:),'thlxy' ,'Temperature' ,'K' ,'tt' )
call ncinfo(ncstatxy( 5,:),'qtxy' ,'Moisture' ,'kg/kg' ,'tt' )
call ncinfo(ncstatxy( 6,:),'pxy' ,'Pressure' ,'kgm/s^2','tt' )
call ncinfo(ncstatxy( 7,:),'upwpxy' ,'Mom. flux' ,'m^2/s^2','mt' )
call ncinfo(ncstatxy( 8,:),'wpthlpxy','Heat flux' ,'Km/s' ,'mt' )
call ncinfo(ncstatxy( 9,:),'vpwpxy' ,'Mom. flux' ,'Km/s' ,'mt' )
call ncinfo(ncstatxy(10,:),'usgsxy' ,'SGS mom. flux' ,'m^2/s^2','mt' )
call ncinfo(ncstatxy(11,:),'thlsgsxy','SGS heat flux' ,'Km/s' ,'mt' )
call ncinfo(ncstatxy(12,:),'vsgsxy' ,'SGS mom. flux' ,'m^2/s^2','mt' )
call ncinfo(ncstatxy(13,:),'uwxyik' ,'Advective mom. flux' ,'m^2/s^2','mt' )
call ncinfo(ncstatxy(14,:),'wthlxy' ,'Advective heat flux' ,'K m/s' ,'mt' )
call ncinfo(ncstatxy(15,:),'vwxy' ,'Advective mom. flux' ,'m^2/s^2','mt' )
if (myid==0) then
call open_nc(xyname, ncidxy, nrecxy, n3=khigh-klow+1)
if (nrecxy==0) then
call define_nc( ncidxy, 1, tncstatxy)
call writestat_dims_nc(ncidxy)
end if
call define_nc( ncidxy, nstatxy, ncstatxy)
end if
end if
!> Generate time, y and x averaged NetCDF: xytdump.xxx.nc
if (lxytdump) then
xytname(9:11) = cexpnr
call ncinfo(tncstatxyt(1,:),'time' ,'Time' ,'s' ,'time')
call ncinfo(ncstatxyt( 1,:),'uxyt' ,'Streamwise velocity' ,'m/s' ,'tt' )
call ncinfo(ncstatxyt( 2,:),'vxyt' ,'Spanwise velocity' ,'m/s' ,'tt' )
call ncinfo(ncstatxyt( 3,:),'wxyt' ,'Vertical velocity' ,'m/s' ,'mt' )
call ncinfo(ncstatxyt( 4,:),'thlxyt' ,'Temperature' ,'K' ,'tt' )
call ncinfo(ncstatxyt( 5,:),'qtxyt' ,'Moisture' ,'kg/kg' ,'tt' )
call ncinfo(ncstatxyt( 6,:),'pxyt' ,'Pressure' ,'kgm/s^2','tt' )
call ncinfo(ncstatxyt( 7,:),'upwpxyt' ,'Turbulent mom. flux' ,'m^2/s^2','mt' )
call ncinfo(ncstatxyt( 8,:),'wpthlpxyt' ,'Turbulent heat flux' ,'K m/s' ,'mt' )
call ncinfo(ncstatxyt( 9,:),'vpwpxyt' ,'Turbulent mom. flux' ,'m^2/s^2','mt' )
call ncinfo(ncstatxyt( 10,:),'upvpxyt' ,'Turbulent mom. flux' ,'m^2/s^2','mt' )
call ncinfo(ncstatxyt( 11,:),'uwxyt' ,'Kinematic mom. flux' ,'m^2/s^2','mt' )
call ncinfo(ncstatxyt( 12,:),'wthlxyt' ,'Kinematic heat flux' ,'K m/s' ,'mt' )
call ncinfo(ncstatxyt( 13,:),'uvxyt' ,'Kinematic mom. flux' ,'m^2/s^2','mt' )
call ncinfo(ncstatxyt( 14,:),'vwxyt' ,'Kinematic mom. flux' ,'m^2/s^2','mt' )
call ncinfo(ncstatxyt( 15,:),'wwxyt' ,'Kinematic mom. flux' ,'m^2/s^2','mt' )
call ncinfo(ncstatxyt( 16,:),'usgsxyt' ,'SGS mom. flux' ,'m^2/s^2','mt' )
call ncinfo(ncstatxyt( 17,:),'thlsgsxyt' ,'SGS heat flux' ,'K m/s' ,'mt' )
call ncinfo(ncstatxyt( 18,:),'vsgsxyt' ,'SGS mom. flux' ,'K m/s' ,'mt' )
call ncinfo(ncstatxyt( 19,:),'thlpthlptxy','Temp. variance' ,'K^2' ,'tt' )
call ncinfo(ncstatxyt( 20,:),'upuptxyc' ,'u variance' ,'m^2/s^2','tt' )
call ncinfo(ncstatxyt( 21,:),'vpvptxyc' ,'v variance' ,'m^2/s^2','tt' )
call ncinfo(ncstatxyt( 22,:),'wpwptxyc' ,'w variance' ,'m^2/s^2','tt' )
call ncinfo(ncstatxyt( 23,:),'tketxyc' ,'tke' ,'m^2/s^2','tt' )
if (myid==0) then
call open_nc(xytname, ncidxyt, nrecxyt, n3=khigh-klow+1)
if (nrecxyt==0) then
call define_nc( ncidxyt, 1, tncstatxyt)
call writestat_dims_nc(ncidxyt)
end if
call define_nc( ncidxyt, nstatxyt, ncstatxyt)
end if
end if
!> Generate time averaged NetCDF: tdump.xxx.nc
if (ltdump) then
tname(7:9) = cmyid
tname(11:13) = cexpnr
call ncinfo(tncstatt(1,:),'time' ,'Time' ,'s' ,'time')
call ncinfo(ncstatt( 1,:),'ut' ,'Streamwise velocity' ,'m/s' ,'mttt' )
call ncinfo(ncstatt( 2,:),'vt' ,'Spanwise velocity' ,'m/s' ,'tmtt' )
call ncinfo(ncstatt( 3,:),'wt' ,'Vertical velocity' ,'m/s' ,'ttmt' )
call ncinfo(ncstatt( 4,:),'thlt' ,'Temperature' ,'K' ,'tttt' )
call ncinfo(ncstatt( 5,:),'qtt' ,'Moisture' ,'kg/kg' ,'tttt' )
call ncinfo(ncstatt( 6,:),'pt' ,'Pressure' ,'kgm/s^2','tttt' )
call ncinfo(ncstatt( 7,:),'sca1t' ,'Concentration field 1' ,'g/m^3' ,'tttt' )
call ncinfo(ncstatt( 8,:),'sca2t' ,'Concentration field 2' ,'g/m^3' ,'tttt' )
call ncinfo(ncstatt( 9,:),'sca3t' ,'Concentration field 3' ,'g/m^3' ,'tttt' )
call ncinfo(ncstatt(10,:),'sca4t' ,'Concentration field 4' ,'g/m^3' ,'tttt' )
call ncinfo(ncstatt(11,:),'PSS' ,'PSS defect' ,'gm/s' ,'tttt' )
call ncinfo(ncstatt(12,:),'upwpt' ,'Turbulent momentum flux' ,'m^2/s^2','mtmt' )
call ncinfo(ncstatt(13,:),'vpwpt' ,'Turbulent momentum flux' ,'m^2/s^2','tmmt' )
call ncinfo(ncstatt(14,:),'upvpt' ,'Turbulent momentum flux' ,'m^2/s^2','mmtt' )
call ncinfo(ncstatt(15,:),'wpthlpt' ,'Turbulent heat flux' ,'K m/s' ,'ttmt' )
call ncinfo(ncstatt(16,:),'wpsca1pt' ,'Turbulent flux 1' ,'gm/s' ,'ttmt' )
call ncinfo(ncstatt(17,:),'wpsca2pt' ,'Turbulent flux 2' ,'gm/s' ,'ttmt' )
call ncinfo(ncstatt(18,:),'wpsca3pt' ,'Turbulent flux 3' ,'gm/s' ,'ttmt' )
call ncinfo(ncstatt(19,:),'wpsca4pt' ,'Turbulent flux 4' ,'gm/s' ,'ttmt' )
call ncinfo(ncstatt(20,:),'thlpthlpt','Temperature variance' ,'K^2' ,'tttt' )
call ncinfo(ncstatt(21,:),'upuptc' ,'u variance' ,'m^2/s^2','tttt' )
call ncinfo(ncstatt(22,:),'vpvptc' ,'v variance' ,'m^2/s^2','tttt' )
call ncinfo(ncstatt(23,:),'wpwptc' ,'w variance' ,'m^2/s^2','tttt' )
call ncinfo(ncstatt(24,:),'tketc' ,'TKE' ,'m^2/s^2','tttt' )
call ncinfo(ncstatt(25,:),'sca1psca1pt','Concentration variance 1' ,'g^2/m^6','tttt' )
call ncinfo(ncstatt(26,:),'sca2psca2pt','Concentration variance 2' ,'g^2/m^6','tttt' )
call ncinfo(ncstatt(27,:),'sca3psca3pt','Concentration variance 3' ,'g^2/m^6','tttt' )
call ncinfo(ncstatt(28,:),'sca4psca4pt','Concentration variance 4' ,'g^2/m^6','tttt' )
call ncinfo(ncstatt(29,:),'sv1sgs' ,'SGS flux 1' ,'gm/s' ,'ttmt' )
call ncinfo(ncstatt(30,:),'sv2sgs' ,'SGS flux 2' ,'gm/s' ,'ttmt' )
call ncinfo(ncstatt(31,:),'sv3sgs' ,'SGS flux 3' ,'gm/s' ,'ttmt' )
call ncinfo(ncstatt(32,:),'sv4sgs' ,'SGS flux 4' ,'gm/s' ,'ttmt' )
! call ncinfo(ncstatt(33,:),'sca1t_max','Max concentration field 1' ,'g/m^3' ,'tttt' )
! call ncinfo(ncstatt(34,:),'sca2t_max','Max concentration field 2' ,'g/m^3' ,'tttt' )
! call ncinfo(ncstatt(35,:),'sca3t_max','Max concentration field 3' ,'g/m^3' ,'tttt' )
! call ncinfo(ncstatt(36,:),'sca4t_max','Max concentration field 4' ,'g/m^3' ,'tttt' )
! if (myid==0) then
call open_nc(tname, ncidt, nrect, n1=imax, n2=jmax, n3=khigh-klow+1)
if (nrect==0) then
call define_nc( ncidt, 1, tncstatt)
call writestat_dims_nc(ncidt)
end if
call define_nc( ncidt, nstatt, ncstatt)
! end if
end if
!> Generate time, y and x averaged NetCDF for tke budget: tkedump.xxx.nc
if (ltkedump) then
tkename(9:11) = cexpnr
call ncinfo(tncstattke(1,:),'time' ,'Time' ,'s' ,'time')
call ncinfo(ncstattke( 1,:),'p_b' ,'p_bant production or consumption term', 'm^2/s^3','tt' )
call ncinfo(ncstattke( 2,:),'t_p' ,'total viscous transport (?)' , 'm^2/s^3','tt' )
call ncinfo(ncstattke( 3,:),'adv' ,'Advection by mean wind' , 'm^2/s^3','tt' )
call ncinfo(ncstattke( 4,:),'t_t' ,'Total turb???' , 'm^2/s^3','tt' )
call ncinfo(ncstattke( 5,:),'t_sgs','total SGS term' , 'm^2/s^3','tt' )
call ncinfo(ncstattke( 6,:),'p_t' ,'Shear production term' , 'm^2/s^3','tt' )
call ncinfo(ncstattke( 7,:),'t_v' ,'Resolved viscous dissipation term' , 'm^2/s^3','tt' )
call ncinfo(ncstattke( 8,:),'d_sgs','SGS dissipation term' , 'm^2/s^3','tt' )
if (myid==0) then
call open_nc(tkename, ncidtke, nrectke, n3=khigh-klow+1)
if (nrectke==0) then
call define_nc( ncidtke, 1, tncstattke)
call writestat_dims_nc(ncidtke)
end if
call define_nc( ncidtke, nstattke, ncstattke)
endif !myid==0
endif
!> Generate sliced NetCDF: slicedump.xxx.xxx.nc
if (lslicedump) then
slicename(11:13) = cmyid
slicename(15:17) = cexpnr
call ncinfo(tncstatslice(1,:),'time' ,'Time' ,'s' ,'time')
call ncinfo(ncstatslice( 1,:),'sca_kb1' ,'Scalar field at kb', '-', 'tt0t')
call ncinfo(ncstatslice( 2,:),'sca_ave1' ,'Averaged scalar field over canyon', '-', 'tt0t')
call ncinfo(ncstatslice( 3,:),'sca_kb2' ,'Scalar field at kb+1', '-', 'tt0t')
call ncinfo(ncstatslice( 4,:),'sca_ave2' ,'Averaged scalar field over canyon', '-', 'tt0t')
call ncinfo(ncstatslice( 5,:),'sca_kb3' ,'Scalar field at kb+1', '-', 'tt0t')
call ncinfo(ncstatslice( 6,:),'sca_ave3' ,'Averaged scalar field over canyon', '-', 'tt0t')
call ncinfo(ncstatslice( 7,:),'u_kb' ,'Streamwise velocity at kb', '-', 'mt0t')
call ncinfo(ncstatslice( 8,:),'v_kb' ,'Spanwise velocity at kb', '-', 'tm0t')
call open_nc(slicename, ncidslice, nrecslice, n1=imax, n2=jmax)
if (nrecslice==0) then
call define_nc( ncidslice, 1, tncstatslice)
call writestat_dims_nc(ncidslice)
end if
call define_nc( ncidslice, nstatslice, ncstatslice)
end if
!> Set times to zero so works for warm starts... could have issues with warmstarts here...
tsamplep = 0.
tstatsdumpp = 0.
end subroutine initstatsdump