Skip to content

Netcdf parallel output extensions and bug repairs #1476

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 15 commits into from
Oct 1, 2025
Merged

Netcdf parallel output extensions and bug repairs #1476

Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
c9983b0
updated wav_history_mod.F90 with netcdf PIO bug fixes
kestonsmith-noaa Jul 22, 2025
5a6963b
Removed unnecessary comments from wav_history_mod.F90
kestonsmith-noaa Jul 23, 2025
95e5adb
Merge remote-tracking branch 'upstream/dev/ufs-weather-model' into Ne…
kestonsmith-noaa Aug 11, 2025
d30fd89
Changes to write_var3d_transpose to address bug
kestonsmith-noaa Aug 20, 2025
1911949
replace write_var3d with DW's extended write_var3d including global flag
kestonsmith-noaa Aug 20, 2025
71cfb59
Modified new write_var3d to handle directional NaN correctly in
kestonsmith-noaa Aug 21, 2025
732253e
Removed obsolete subroutines wav_history_mod
kestonsmith-noaa Aug 21, 2025
1c736cd
Reorder subroutines to match origonal order
kestonsmith-noaa Aug 21, 2025
d37e368
Merge remote-tracking branch 'upstream/dev/ufs-weather-model' into Ne…
kestonsmith-noaa Aug 21, 2025
18edc47
moved frequency from double precision to real
kestonsmith-noaa Aug 21, 2025
e80cd31
Revert "Reorder subroutines to match origonal order"
kestonsmith-noaa Aug 22, 2025
f64b92f
moved write_var3d before writevar2d
kestonsmith-noaa Aug 22, 2025
b7bc13f
Merge remote-tracking branch 'upstream/dev/ufs-weather-model' into Ne…
kestonsmith-noaa Aug 25, 2025
da8bba5
removed tabs
kestonsmith-noaa Aug 25, 2025
89710b7
fixed spacing in wav_history_mod.F90
kestonsmith-noaa Aug 25, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
126 changes: 105 additions & 21 deletions model/src/wav_history_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ module wav_history_mod
real, allocatable, target :: var3dm(:,:)
real, allocatable, target :: var3dp(:,:)
real, allocatable, target :: var3dk(:,:)
real, allocatable, target :: var3dnk(:,:) ! WN : node vs frequency dimensions reversed

! output variable for (nx,ny,nz) fields
real, pointer :: var3d(:,:)
Expand All @@ -46,6 +47,7 @@ module wav_history_mod
type(io_desc_t) :: iodesc3dm !m-axis variables
type(io_desc_t) :: iodesc3dp !p-axis variables
type(io_desc_t) :: iodesc3dk !k-axis variables
type(io_desc_t) :: iodesc3dnk !nk-axis variables for wavenumber (nk.ne.k-axis)

! variable attributes
type :: varatts
Expand Down Expand Up @@ -96,6 +98,13 @@ subroutine write_history ( timen )
use w3timemd , only : set_user_timestring
use w3odatmd , only : time_origin, calendar_name, elapsed_secs
use w3odatmd , only : user_histfname

use w3adatmd , only : wn

!needed for frequency calculation
use constants , only : TPIINV
USE W3GDATMD , ONLY : SIG

!TODO: use unstr_mesh from wav_shr_mod; currently fails due to CI
!use wav_shr_mod , only : unstr_mesh

Expand All @@ -111,13 +120,17 @@ subroutine write_history ( timen )
! indicator logfile
character(len=256) :: log_fname ! log file name
integer :: log_unit = 28888 ! unit number for log file


integer :: n, xtid, ytid, xeid, ztid, stid, mtid, ptid, ktid, timid, nmode, nktid,k

integer :: n, xtid, ytid, xeid, ztid, stid, mtid, ptid, ktid, timid, nmode
integer :: len_s, len_m, len_p, len_k
logical :: s_axis = .false., m_axis = .false., p_axis = .false., k_axis = .false.
integer :: len_s, len_m, len_p, len_k, len_nk
logical :: s_axis = .false., m_axis = .false., p_axis = .false., k_axis = .false., nk_axis=.false.

integer :: lmap(nseal_cpl)

real, allocatable :: freq_wn(:), freq_ef(:)

! -------------------------------------------------------------
! create the netcdf file
! -------------------------------------------------------------
Expand All @@ -132,7 +145,7 @@ subroutine write_history ( timen )
fname = trim(FNMGRD)//trim(user_histfname)//trim(user_timestring)//'.nc'
log_fname = trim(FNMGRD)//'log.'//trim(user_histfname)//trim(user_timestring)//'.nc.txt'
end if

pioid%fh = -1
nmode = pio_clobber
! only applies to classic NETCDF files.
Expand All @@ -147,6 +160,8 @@ subroutine write_history ( timen )
len_m = p2msf(3)-p2msf(2) + 1 ! ?
len_p = usspf(2) ! partitions
len_k = e3df(3,1) - e3df(2,1) + 1 ! frequencies
len_nk = nk ! frequencies for wavenumber


! define the dimensions required for the requested gridded fields
do n = 1,size(outvars)
Expand All @@ -155,14 +170,19 @@ subroutine write_history ( timen )
if(trim(outvars(n)%dims) == 'm')m_axis = .true.
if(trim(outvars(n)%dims) == 'p')p_axis = .true.
if(trim(outvars(n)%dims) == 'k')k_axis = .true.
end if
if(trim(outvars(n)%dims) == 'nk')nk_axis = .true.
end if
end do

! allocate arrays if needed
if (s_axis) allocate(var3ds(1:nseal_cpl,len_s))
if (m_axis) allocate(var3dm(1:nseal_cpl,len_m))
if (p_axis) allocate(var3dp(1:nseal_cpl,len_p))
if (k_axis) allocate(var3dk(1:nseal_cpl,len_k))
if (nk_axis) allocate(var3dnk(1:nseal_cpl,len_nk))

if ( k_axis ) allocate(freq_ef(1:len_k))
if ( nk_axis ) allocate(freq_wn(1:len_nk))

ierr = pio_def_dim(pioid, 'nx', nx, xtid)
ierr = pio_def_dim(pioid, 'ny', ny, ytid)
Expand All @@ -171,7 +191,10 @@ subroutine write_history ( timen )
if (s_axis) ierr = pio_def_dim(pioid, 'noswll', len_s, stid)
if (m_axis) ierr = pio_def_dim(pioid, 'nm' , len_m, mtid)
if (p_axis) ierr = pio_def_dim(pioid, 'np' , len_p, ptid)
if (k_axis) ierr = pio_def_dim(pioid, 'freq' , len_k, ktid)

if (k_axis) ierr = pio_def_dim(pioid, 'nf_ef' , len_k, ktid)
if (nk_axis) ierr = pio_def_dim(pioid, 'nf_wn' , len_nk, nktid)

if (gtype .eq. ungtype) then
ierr = pio_def_dim(pioid, 'ne' , ntri, xeid)
ierr = pio_def_dim(pioid, 'nn' , 3, ztid)
Expand Down Expand Up @@ -207,6 +230,18 @@ subroutine write_history ( timen )
ierr = pio_put_att(pioid, varid, 'long_name', 'node connectivity')
end if

! define the frequency axis variables for wavenumber(wn) and spectra(ef)
if (k_axis) then
ierr = pio_def_var(pioid, 'freq_ef', PIO_REAL, (/ktid/), varid)
call handle_err(ierr,'def_freq')
ierr = pio_put_att(pioid, varid, 'units', 's-1')
end if
if (nk_axis) then
ierr = pio_def_var(pioid, 'freq_wn', PIO_REAL, (/nktid/), varid)
call handle_err(ierr,'def_freq')
ierr = pio_put_att(pioid, varid, 'units', 's-1')
end if

Copy link
Contributor

@DeniseWorthen DeniseWorthen Aug 18, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please fix the indenting on these lines (2-space indents for if-statements).

Also, is there a reason they are PIO_DOUBLE? They're just real(kind=4) values, correct?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@kestonsmith-noaa Why are these axis values required to be double precision?

! define the variables
dimid3(1:2) = (/xtid, ytid/)
dimid4(1:2) = (/xtid, ytid/)
Expand All @@ -223,6 +258,9 @@ subroutine write_history ( timen )
else if (trim(outvars(n)%dims) == 'k') then
dimid4(3:4) = (/ktid, timid/)
dimid => dimid4
else if (trim(outvars(n)%dims) == 'nk') then ! wavenumber
dimid4(3:4) = (/nktid, timid/)
dimid => dimid4
else
dimid3(3) = timid
dimid => dimid3
Expand All @@ -244,6 +282,7 @@ subroutine write_history ( timen )
if (m_axis)call wav_pio_initdecomp(len_m, iodesc3dm)
if (p_axis)call wav_pio_initdecomp(len_p, iodesc3dp)
if (k_axis)call wav_pio_initdecomp(len_k, iodesc3dk)
if (nk_axis)call wav_pio_initdecomp(len_nk, iodesc3dnk)

! write the time and spatial axis values (lat,lon,time)
ierr = pio_inq_varid(pioid, 'lat', varid)
Expand All @@ -260,7 +299,27 @@ subroutine write_history ( timen )
call handle_err(ierr, 'inquire variable time ')
ierr = pio_put_var(pioid, varid, (/1/), real(elapsed_secs,8))
call handle_err(ierr, 'put time')

if (k_axis) then
do k=1,len_k
freq_ef(k)=SIG( e3df(2,1) + k -1 ) * TPIINV
enddo
ierr = pio_inq_varid(pioid, 'freq_ef', varid)
call handle_err(ierr, 'inquire variable freq EF')
ierr = pio_put_var(pioid, varid, freq_ef(1:len_k) )
call handle_err(ierr, 'put freq EF')
end if

if (nk_axis) then
do k=1,len_nk
freq_wn(k)=SIG( e3df(2,1) + k -1 ) * TPIINV
enddo
ierr = pio_inq_varid(pioid, 'freq_wn', varid)
call handle_err(ierr, 'inquire variable freq WN')
ierr = pio_put_var(pioid, varid, freq_wn(1:len_nk) )
call handle_err(ierr, 'put freq WN')
end if

if (gtype .eq. ungtype) then
ierr = pio_inq_varid(pioid, 'nconn', varid)
call handle_err(ierr, 'inquire variable nconn ')
Expand All @@ -286,15 +345,15 @@ subroutine write_history ( timen )

! write the requested variables
do n = 1,size(outvars)
vname = trim(outvars(n)%var_name)
vname = trim(outvars(n)%var_name)
if (trim(outvars(n)%dims) == 's') then
var3d => var3ds
! Group 4
if(vname .eq. 'PHS') call write_var3d(iodesc3ds, vname, phs (1:nseal_cpl,0:noswll) )
if(vname .eq. 'PTP') call write_var3d(iodesc3ds, vname, ptp (1:nseal_cpl,0:noswll) )
if(vname .eq. 'PLP') call write_var3d(iodesc3ds, vname, plp (1:nseal_cpl,0:noswll) )
if(vname .eq. 'PDIR') call write_var3d(iodesc3ds, vname, pdir (1:nseal_cpl,0:noswll), fldir='true' )
if(vname .eq. 'PSI') call write_var3d(iodesc3ds, vname, psi (1:nseal_cpl,0:noswll), fldir='true' )
if(vname .eq. 'PSI') call write_var3d(iodesc3ds, vname, psi (1:nseal_cpl,0:noswll) )
if(vname .eq. 'PWS') call write_var3d(iodesc3ds, vname, pws (1:nseal_cpl,0:noswll) )
if(vname .eq. 'PDP') call write_var3d(iodesc3ds, vname, pthp0 (1:nseal_cpl,0:noswll), fldir='true' )
if(vname .eq. 'PQP') call write_var3d(iodesc3ds, vname, pqp (1:nseal_cpl,0:noswll) )
Expand All @@ -317,10 +376,15 @@ subroutine write_history ( timen )
if (vname .eq. 'USSPX') call write_var3d(iodesc3dp, vname, ussp (1:nseal_cpl, 1:usspf(2)) )
if (vname .eq. 'USSPY') call write_var3d(iodesc3dp, vname, ussp (1:nseal_cpl,nk+1:nk+usspf(2)) )

else if (trim(outvars(n)%dims) == 'nk') then ! freq + 1 axis for wavenumber
var3d => var3dnk
if(vname .eq. 'WN') call write_var3d(iodesc3dnk, vname, transpose(wn(1:nk,1:nsea)), global='true')

else if (trim(outvars(n)%dims) == 'k') then ! freq axis
var3d => var3dk
! Group 3
if(vname .eq. 'EF') call write_var3d(iodesc3dk, vname, ef (1:nseal_cpl,e3df(2,1):e3df(3,1)) )
if(vname .eq. 'EF') call write_var3d(iodesc3dk, vname, ef (1:nseal_cpl,e3df(2,1):e3df(3,1)) )

if(vname .eq. 'TH1M') call write_var3d(iodesc3dk, vname, ef (1:nseal_cpl,e3df(2,2):e3df(3,2)) )
if(vname .eq. 'STH1M') call write_var3d(iodesc3dk, vname, ef (1:nseal_cpl,e3df(2,3):e3df(3,3)) )
if(vname .eq. 'TH2M') call write_var3d(iodesc3dk, vname, ef (1:nseal_cpl,e3df(2,4):e3df(3,4)) )
Expand Down Expand Up @@ -355,7 +419,7 @@ subroutine write_history ( timen )
if (vname .eq. 'T01') call write_var2d(vname, t01 (1:nseal_cpl) )
if (vname .eq. 'FP0') call write_var2d(vname, fp0 (1:nseal_cpl) )
if (vname .eq. 'THM') call write_var2d(vname, thm (1:nseal_cpl), fldir='true' )
if (vname .eq. 'THS') call write_var2d(vname, ths (1:nseal_cpl), fldir='true' )
if (vname .eq. 'THS') call write_var2d(vname, ths (1:nseal_cpl) )
if (vname .eq. 'THP0') call write_var2d(vname, thp0 (1:nseal_cpl), fldir='true' )
if (vname .eq. 'HSIG') call write_var2d(vname, hsig (1:nseal_cpl) )
if (vname .eq. 'STMAXE') call write_var2d(vname, stmaxe (1:nseal_cpl) )
Expand All @@ -372,8 +436,9 @@ subroutine write_history ( timen )
if(vname .eq. 'PNR') call write_var2d(vname, pnr (1:nseal_cpl) )

! Group 5
if (vname .eq. 'USTX') call write_var2d(vname, ust (1:nseal_cpl)*asf(1:nseal_cpl), dir=cos(ustdir(1:nseal_cpl)), usemask='true')
if (vname .eq. 'USTY') call write_var2d(vname, ust (1:nseal_cpl)*asf(1:nseal_cpl), dir=sin(ustdir(1:nseal_cpl)), usemask='true')
if (vname .eq. 'USTX') call write_var2d(vname, ust (1:nsea)*asf(1:nsea), dir=cos(ustdir(1:nsea)), init0='false', global='true')
if (vname .eq. 'USTY') call write_var2d(vname, ust (1:nsea)*asf(1:nsea), dir=sin(ustdir(1:nsea)), init0='false', global='true')

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Have you run in debug mode and verified all the changed variables (ustx,usty, ths,psi) work as expected?

if (vname .eq. 'CHARN') call write_var2d(vname, charn (1:nseal_cpl) )
if (vname .eq. 'CGE') call write_var2d(vname, cge (1:nseal_cpl) )
if (vname .eq. 'PHIAW') call write_var2d(vname, phiaw (1:nseal_cpl), init2='true')
Expand Down Expand Up @@ -442,13 +507,18 @@ subroutine write_history ( timen )
if (m_axis) deallocate(var3dm)
if (p_axis) deallocate(var3dp)
if (k_axis) deallocate(var3dk)
if (nk_axis) deallocate(var3dnk)

if (k_axis ) deallocate(freq_ef)
if (nk_axis) deallocate(freq_wn)

call pio_freedecomp(pioid,iodesc2d)
call pio_freedecomp(pioid,iodesc2dint)
if (s_axis) call pio_freedecomp(pioid, iodesc3ds)
if (m_axis) call pio_freedecomp(pioid, iodesc3dm)
if (p_axis) call pio_freedecomp(pioid, iodesc3dp)
if (k_axis) call pio_freedecomp(pioid, iodesc3dk)
if (nk_axis) call pio_freedecomp(pioid, iodesc3dnk)

call pio_closefile(pioid)

Expand Down Expand Up @@ -586,21 +656,24 @@ end subroutine write_var2d
!! @param[in] var the variable array
!! @param[in] init2 a flag for a second initialization type, optional
!! @param[in] fldir a flag for unit conversion for direction, optional
!!
!! @param[in] global a flag for a global variable, optional

!> author [email protected]
!> @date 08-26-2024
subroutine write_var3d(iodesc, vname, var, init2, fldir)
subroutine write_var3d(iodesc, vname, var, init2, fldir, global)

type(io_desc_t), intent(inout) :: iodesc
character(len=*), intent(in) :: vname
real , intent(in) :: var(:,:)
character(len=*), optional, intent(in) :: init2
character(len=*), optional, intent(in) :: fldir
character(len=*), optional, intent(in) :: global

! local variables
real, allocatable, dimension(:) :: varloc
logical :: linit2, lfldir
logical :: linit2, lfldir, lglobal
integer :: lb, ub
integer :: k

linit2 = .false.
if (present(init2)) then
Expand All @@ -610,6 +683,10 @@ subroutine write_var3d(iodesc, vname, var, init2, fldir)
if (present(fldir)) then
lfldir = (trim(fldir) == "true")
end if
lglobal = .false.
if (present(global)) then
lglobal = (trim(global) == "true")
end if

lb = lbound(var,2)
ub = ubound(var,2)
Expand All @@ -619,14 +696,21 @@ subroutine write_var3d(iodesc, vname, var, init2, fldir)
do jsea = 1,nseal_cpl
call init_get_isea(isea, jsea)
! initialization
varloc(:) = var(jsea,:)
if (lglobal) then
varloc(:) = var(isea,:)
else
varloc(:) = var(jsea,:)
end if

if (mapsta(mapsf(isea,2),mapsf(isea,1)) < 0) varloc(:) = undef
if (linit2) then
if (mapsta(mapsf(isea,2),mapsf(isea,1)) == 2) varloc(:) = undef
end if
if (lfldir) then
if (mapsta(mapsf(isea,2),mapsf(isea,1)) > 0 ) then
varloc(:) = mod(630. - rade*varloc(:), 360.)
do k=1,size(varloc,1)
if (varloc(k).ne.undef) varloc(k)=mod( 630.-rade*varloc(k), 360.)
enddo
end if
end if
var3d(jsea,:) = varloc(:)
Expand All @@ -640,6 +724,7 @@ subroutine write_var3d(iodesc, vname, var, init2, fldir)
deallocate(varloc)
end subroutine write_var3d


!===============================================================================
!> Scan through all possible fields to determine a list of requested variables
!!
Expand Down Expand Up @@ -812,8 +897,7 @@ subroutine define_fields (gridoutdefs)
varatts( "STH1M", "STH1M ", "Directional spreading from a1,b2 ", "deg ", "k ", .false.) , &
varatts( "TH2M ", "TH2M ", "Mean wave direction from a2,b2 ", "deg ", "k ", .false.) , &
varatts( "STH2M", "STH2M ", "Directional spreading from a2,b2 ", "deg ", "k ", .false.) , &
!TODO: has reverse indices (nk,nsea)
varatts( "WN ", "WN ", "Wavenumber array ", "m-1 ", "k ", .false.) &
varatts( "WN ", "WN ", "Wavenumber array ", "m-1 ", "nk", .false.) &
]

! 4 Spectral Partition Parameters
Expand Down Expand Up @@ -842,7 +926,7 @@ subroutine define_fields (gridoutdefs)
varatts( "UST ", "USTX ", "Friction velocity x ", "m s-1 ", " ", .false.) , &
varatts( "UST ", "USTY ", "Friction velocity y ", "m s-1 ", " ", .false.) , &
varatts( "CHA ", "CHARN ", "Charnock parameter ", "nd ", " ", .false.) , &
varatts( "CGE ", "CGE ", "Energy flux ", "kW m-1 ", " ", .false.) , &
varatts( "CGE ", "CGE ", "Energy flux ", "W m-1 ", " ", .false.) , &
varatts( "FAW ", "PHIAW ", "Air-sea energy flux ", "W m-2 ", " ", .false.) , &
varatts( "TAW ", "TAUWIX ", "Net wave-supported stress x ", "m2 s-2 ", " ", .false.) , &
varatts( "TAW ", "TAUWIY ", "Net wave-supported stress y ", "m2 s-2 ", " ", .false.) , &
Expand Down Expand Up @@ -913,7 +997,7 @@ subroutine define_fields (gridoutdefs)

! 9 Numerical diagnostics
gridoutdefs(9,1:5) = [ &
varatts( "DTD ", "DTDYN ", "Average time step in integration ", "min ", " ", .false.) , &
varatts( "DTD ", "DTDYN ", "Average time step in integration ", "s ", " ", .false.) , &
varatts( "FC ", "FCUT ", "Cut-off frequency ", "s-1 ", " ", .false.) , &
varatts( "CFX ", "CFLXYMAX ", "Max. CFL number for spatial advection ", "nd ", " ", .false.) , &
varatts( "CFD ", "CFLTHMAX ", "Max. CFL number for theta-advection ", "nd ", " ", .false.) , &
Expand Down
Loading