diff --git a/model/src/CMakeLists.txt b/model/src/CMakeLists.txt index 9242409e7..9eaa67441 100644 --- a/model/src/CMakeLists.txt +++ b/model/src/CMakeLists.txt @@ -1,4 +1,3 @@ - # Open switch file file(STRINGS ${CMAKE_BINARY_DIR}/switch switch_strings) separate_arguments(switches UNIX_COMMAND ${switch_strings}) @@ -85,7 +84,7 @@ elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$") set(compile_flags -g -fno-second-underscore -ffree-line-length-none) set(compile_flags_release -O3) set(compile_flags_debug -Wall -fcheck=all -ffpe-trap=invalid,zero,overflow -frecursive -fbacktrace) - + if(${CMAKE_Fortran_COMPILER_VERSION} VERSION_GREATER_EQUAL 10) target_compile_options(ww3_lib PUBLIC -fallow-argument-mismatch) endif() @@ -94,7 +93,7 @@ elseif(CMAKE_Fortran_COMPILER_ID MATCHES "PGI") set(compile_flags -g -i4 -r4 -Kieee) set(compile_flags_release -O3) set(compile_flags_debug -O0 -Mbounds -Mchkfpstk -Mchkstk -Mdalign -Mdclchk -Mdepchk -Miomutex -Ktrap=fp -Mrecursive -traceback) - + if(${CMAKE_Fortran_COMPILER_VERSION} VERSION_GREATER_EQUAL 10) target_compile_options(ww3_lib PUBLIC -fallow-argument-mismatch -fallow-invalid-boz) endif() @@ -152,7 +151,7 @@ if(UFS_CAP) elseif(UFS_CAP STREQUAL "NUOPC_MESH") set(cap_src ${nuopc_mesh_cap_src}) endif() - + target_sources(ww3_lib PRIVATE ${cap_src}) target_link_libraries(ww3_lib PUBLIC esmf) # Don't build executables when building WW3 ESMF library @@ -231,7 +230,7 @@ install( install(FILES ${CMAKE_BINARY_DIR}/switch DESTINATION ${CMAKE_INSTALL_PREFIX}) install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/mod DESTINATION ${CMAKE_INSTALL_PREFIX}) - + export(EXPORT WW3Exports NAMESPACE WW3:: diff --git a/model/src/cmake/src_list.cmake b/model/src/cmake/src_list.cmake index dcab88a09..825f5d637 100644 --- a/model/src/cmake/src_list.cmake +++ b/model/src/cmake/src_list.cmake @@ -55,8 +55,7 @@ set(ftn_src wmupdtmd.F90 wmwavemd.F90 w3tidemd.F90 - wav_grdout.F90 - w3iogoncdmd.F90 + wav_history_mod.F90 wav_shr_flags.F90 ) @@ -67,6 +66,8 @@ set(nuopc_mesh_cap_src wav_comp_nuopc.F90 wav_import_export.F90 wav_wrapper_mod.F90 + wav_pio_mod.F90 + wav_restart_mod.F90 ) set(esmf_multi_cap_src diff --git a/model/src/w3gridmd.F90 b/model/src/w3gridmd.F90 index 052b1f2c9..aa703d704 100644 --- a/model/src/w3gridmd.F90 +++ b/model/src/w3gridmd.F90 @@ -836,7 +836,7 @@ MODULE W3GRIDMD ! #ifdef W3_ST4 INTEGER :: SWELLFPAR, SDSISO, SDSBRFDF - REAL :: SDSBCHOICE + REAL :: SDSBCHOICE REAL :: ZWND, ALPHA0, Z0MAX, BETAMAX, SINTHP,& ZALP, Z0RAT, TAUWSHELTER, SWELLF, & SWELLF2,SWELLF3,SWELLF4, SWELLF5, & @@ -3280,7 +3280,7 @@ SUBROUTINE W3GRID() JGS_TERMINATE_DIFFERENCE, & JGS_TERMINATE_NORM, & JGS_LIMITER, & - JGS_LIMITER_FUNC, & + JGS_LIMITER_FUNC, & JGS_USE_JACOBI, & JGS_BLOCK_GAUSS_SEIDEL, & JGS_MAXITER, & @@ -3617,7 +3617,7 @@ SUBROUTINE W3GRID() END SELECT IF (FSTOTALIMP .or. FSTOTALEXP) THEN - LPDLIB = .TRUE. + LPDLIB = .TRUE. ENDIF ! IF (SUM(UNSTSCHEMES).GT.1) WRITE(NDSO,1035) diff --git a/model/src/w3initmd.F90 b/model/src/w3initmd.F90 index 9f7c62de8..0f1362277 100644 --- a/model/src/w3initmd.F90 +++ b/model/src/w3initmd.F90 @@ -445,6 +445,10 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD, FLGR2, FLGD, #ifdef W3_UOST USE W3UOSTMD, ONLY: UOST_SETGRID #endif + use w3timemd, only : set_user_timestring + use w3odatmd, only : runtype, restart_from_binary, use_restartnc, user_restfname + use w3odatmd, only : logfile_is_assigned + use wav_restart_mod, only : read_restart !/ #ifdef W3_MPI INCLUDE "mpif.h" @@ -512,7 +516,10 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD, FLGR2, FLGD, #ifdef W3_PDLIB INTEGER :: IScal(1), IPROC #endif + logical :: exists integer :: memunit + character(len=16) :: user_timestring !YYYY-MM-DD-SSSSS + character(len=1024) :: fname !/ !/ ------------------------------------------------------------------- / ! @@ -639,53 +646,52 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD, FLGR2, FLGD, IF (FSTOTALIMP .and. .NOT. LPDLIB) THEN WRITE(NDSE,*) 'IMPTOTAL is selected' WRITE(NDSE,*) 'But PDLIB is not' - CALL FLUSH(NDSE) - STOP + CALL FLUSH(NDSE) + STOP ELSE IF (FSTOTALEXP .and. .NOT. LPDLIB) THEN WRITE(NDSE,*) 'EXPTOTAL is selected' WRITE(NDSE,*) 'But PDLIB is not' - CALL FLUSH(NDSE) - STOP + CALL FLUSH(NDSE) + STOP END IF #ifdef W3_PDLIB IF (B_JGS_BLOCK_GAUSS_SEIDEL .AND. .NOT. B_JGS_USE_JACOBI) THEN WRITE(NDSE,*) 'B_JGS_BLOCK_GAUSS_SEIDEL is used but the Jacobi solver is not choosen' WRITE(NDSE,*) 'Please set JGS_USE_JACOBI .eqv. .true.' - CALL FLUSH(NDSE) - STOP + CALL FLUSH(NDSE) + STOP ENDIF #endif - ! ! 1.c Open files without unpacking MDS ,,, ! - IE = LEN_TRIM(FEXT) - LFILE = 'log.' // FEXT(:IE) - IFL = LEN_TRIM(LFILE) + if (.not. logfile_is_assigned) then + IE = LEN_TRIM(FEXT) + LFILE = 'log.' // FEXT(:IE) + IFL = LEN_TRIM(LFILE) #ifdef W3_SHRD - TFILE = 'test.' // FEXT(:IE) + TFILE = 'test.' // FEXT(:IE) #endif #ifdef W3_DIST - IW = 1 + INT ( LOG10 ( REAL(NAPROC) + 0.5 ) ) - IW = MAX ( 3 , MIN ( 9 , IW ) ) - WRITE (FORMAT,'(A5,I1.1,A1,I1.1,A4)') & - '(A4,I', IW, '.', IW, ',2A)' - WRITE (TFILE,FORMAT) 'test', & - OUTPTS(IMOD)%IAPROC, '.', FEXT(:IE) -#endif - IFT = LEN_TRIM(TFILE) - J = LEN_TRIM(FNMPRE) - ! -#ifndef W3_CESMCOUPLED - IF ( OUTPTS(IMOD)%IAPROC .EQ. OUTPTS(IMOD)%NAPLOG ) & - OPEN (MDS(1),FILE=FNMPRE(:J)//LFILE(:IFL),ERR=888,IOSTAT=IERR) -#endif - ! - IF ( MDS(3).NE.MDS(1) .AND. MDS(3).NE.MDS(4) .AND. TSTOUT ) THEN - INQUIRE (MDS(3),OPENED=OPENED) - IF ( .NOT. OPENED ) OPEN (MDS(3),FILE=FNMPRE(:J)//TFILE(:IFT), ERR=889, & - IOSTAT=IERR) - END IF + IW = 1 + INT ( LOG10 ( REAL(NAPROC) + 0.5 ) ) + IW = MAX ( 3 , MIN ( 9 , IW ) ) + WRITE (FORMAT,'(A5,I1.1,A1,I1.1,A4)') & + '(A4,I', IW, '.', IW, ',2A)' + WRITE (TFILE,FORMAT) 'test', & + OUTPTS(IMOD)%IAPROC, '.', FEXT(:IE) +#endif + IFT = LEN_TRIM(TFILE) + J = LEN_TRIM(FNMPRE) + ! + IF ( OUTPTS(IMOD)%IAPROC .EQ. OUTPTS(IMOD)%NAPLOG ) & + OPEN (MDS(1),FILE=FNMPRE(:J)//LFILE(:IFL),ERR=888,IOSTAT=IERR) + ! + IF ( MDS(3).NE.MDS(1) .AND. MDS(3).NE.MDS(4) .AND. TSTOUT ) THEN + INQUIRE (MDS(3),OPENED=OPENED) + IF ( .NOT. OPENED ) OPEN (MDS(3),FILE=FNMPRE(:J)//TFILE(:IFT), ERR=889, & + IOSTAT=IERR) + END IF + end if ! if (.not. logfile_is_assigned) ! ! 1.d Dataset unit numbers ! @@ -725,6 +731,7 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD, FLGR2, FLGD, ! 2.a Read model definition file ! CALL W3IOGR ( 'READ', NDS(5), IMOD, FEXT ) + IF (GTYPE .eq. UNGTYPE) THEN CALL SPATIAL_GRID CALL NVECTRI @@ -952,40 +959,64 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD, FLGR2, FLGD, ! 3.a Read restart file ! VA(:,:) = 0. + if (use_restartnc) then + if (runtype == 'continue' )then + call set_user_timestring(time,user_timestring) + if (restart_from_binary) then + fname = trim(user_restfname)//trim(user_timestring) + else + fname = trim(user_restfname)//trim(user_timestring)//'.nc' + endif + inquire(file=trim(fname), exist=exists) + if (exists) then + if (restart_from_binary) then + call w3iors('READ', nds(6), sig(nk), imod, filename=trim(fname)) + else + call read_restart(trim(fname), va=va, mapsta=mapsta, mapst2=mapst2) + end if + else + call extcde (60, msg="required restart file " // trim(fname) // " does not exist") + end if + else + call read_restart('none') + ! mapst2 is module variable defined in read of mod_def; maptst is from 2.b above + flcold = .true. + end if + else #ifdef W3_DEBUGCOH - CALL ALL_VA_INTEGRAL_PRINT(IMOD, "Before W3IORS call", 1) + CALL ALL_VA_INTEGRAL_PRINT(IMOD, "Before W3IORS call", 1) #endif #ifdef W3_TIMINGS - CALL PRINT_MY_TIME("Before W3IORS") + CALL PRINT_MY_TIME("Before W3IORS") #endif - CALL W3IORS ( 'READ', NDS(6), SIG(NK), IMOD) + CALL W3IORS ( 'READ', NDS(6), SIG(NK), IMOD) #ifdef W3_TIMINGS - CALL PRINT_MY_TIME("After W3IORS") + CALL PRINT_MY_TIME("After W3IORS") #endif - call print_memcheck(memunit, 'memcheck_____:'//' WW3_INIT SECTION 3a') + call print_memcheck(memunit, 'memcheck_____:'//' WW3_INIT SECTION 3a') #ifdef W3_DEBUGCOH - CALL ALL_VA_INTEGRAL_PRINT(IMOD, "After W3IORS call", 1) -#endif - FLCOLD = RSTYPE.LE.1 .OR. RSTYPE.EQ.4 - IF ( IAPROC .EQ. NAPLOG ) THEN - IF (RSTYPE.EQ.0) THEN - WRITE (NDSO,930) 'cold start (idealized).' - ELSE IF ( RSTYPE .EQ. 1 ) THEN - WRITE (NDSO,930) 'cold start (wind).' - ELSE IF ( RSTYPE .EQ. 4 ) THEN - WRITE (NDSO,930) 'cold start (calm).' - ELSE - WRITE (NDSO,930) 'full restart.' + CALL ALL_VA_INTEGRAL_PRINT(IMOD, "After W3IORS call", 1) +#endif + FLCOLD = RSTYPE.LE.1 .OR. RSTYPE.EQ.4 + IF ( IAPROC .EQ. NAPLOG ) THEN + IF (RSTYPE.EQ.0) THEN + WRITE (NDSO,930) 'cold start (idealized).' + ELSE IF ( RSTYPE .EQ. 1 ) THEN + WRITE (NDSO,930) 'cold start (wind).' + ELSE IF ( RSTYPE .EQ. 4 ) THEN + WRITE (NDSO,930) 'cold start (calm).' + ELSE + WRITE (NDSO,930) 'full restart.' + END IF END IF - END IF #ifdef W3_DEBUGCOH - CALL ALL_VA_INTEGRAL_PRINT(IMOD, "W3INIT, step 4.2", 1) + CALL ALL_VA_INTEGRAL_PRINT(IMOD, "W3INIT, step 4.2", 1) #endif #ifdef W3_TIMINGS - CALL PRINT_MY_TIME("After restart inits") + CALL PRINT_MY_TIME("After restart inits") #endif - + end if ! if (use_restartnc) ! ! 3.b Compare MAPSTA from grid and restart ! @@ -1263,7 +1294,6 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD, FLGR2, FLGD, ! MAPTST = MOD(MAPST2/2,2) MAPST2 = MAPST2 - 2*MAPTST - ! !Li For multi-resolution SMC grid, these 1-NX and 1-NY nested loops !Li may miss the refined cells as they are not 1-1 corresponding to @@ -1337,12 +1367,10 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD, FLGR2, FLGD, CALL SET_IOBDP_PDLIB ENDIF #endif - ! #ifdef W3_DEBUGCOH CALL ALL_VA_INTEGRAL_PRINT(IMOD, "W3INIT, step 8.2", 1) #endif - ! MAPST2 = MAPST2 + 2*MAPTST ! @@ -1396,7 +1424,6 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD, FLGR2, FLGD, ! END DO END DO - ! ! 6. Initialize arrays ---------------------------------------------- / ! Some initialized in W3IORS @@ -1413,7 +1440,7 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD, FLGR2, FLGD, ! ! 7. Write info to log file ----------------------------------------- / ! - IF ( IAPROC .EQ. NAPLOG ) THEN + IF ( IAPROC .EQ. NAPLOG) THEN ! WRITE (NDSO,970) GNAME IF ( FLLEV ) WRITE (NDSO,971) 'Prescribed' @@ -1500,7 +1527,9 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD, FLGR2, FLGD, WRITE (NDSO,990) DTME21 END IF ! - WRITE (NDSO,984) + if (.not. logfile_is_assigned) then + WRITE (NDSO,984) + end if ! END IF ! @@ -2171,6 +2200,7 @@ SUBROUTINE W3MPIO ( IMOD ) #endif USE W3GDATMD, ONLY: GTYPE, UNGTYPE USE CONSTANTS, ONLY: LPDLIB + use w3odatmd, only : restart_from_binary, use_restartnc, use_historync !/ #ifdef W3_MPI INCLUDE "mpif.h" @@ -2202,6 +2232,7 @@ SUBROUTINE W3MPIO ( IMOD ) #ifdef W3_MPIT CHARACTER(LEN=5) :: STRING #endif + logical :: do_rstsetup !/ !/ ------------------------------------------------------------------- / !/ @@ -2225,7 +2256,7 @@ SUBROUTINE W3MPIO ( IMOD ) IROOT = NAPFLD - 1 ! ! - IF ((FLOUT(1) .OR. FLOUT(7)) .and. (.not. LPDLIB)) THEN + IF ((FLOUT(1) .OR. FLOUT(7)) .and. (.not. LPDLIB) .and. (.not. use_historync)) THEN ! ! NRQMAX is the maximum number of output fields that require MPI communication, ! aimed to gather field values stored in each processor into one processor in @@ -4760,7 +4791,7 @@ SUBROUTINE W3MPIO ( IMOD ) CALL EXTCDE (11) END IF ! - END IF ! IF ((FLOUT(1) .OR. FLOUT(7)) .and. (.not. LPDLIB)) THEN + END IF ! IF ((FLOUT(1) .OR. FLOUT(7)) .and. (.not. LPDLIB) .and. (.not. use_historync)) THEN ! ! 2. Set-up for W3IORS ---------------------------------------------- / ! 2.a General preparations @@ -4769,7 +4800,17 @@ SUBROUTINE W3MPIO ( IMOD ) IH = 0 IROOT = NAPRST - 1 ! - IF ((FLOUT(4) .OR. FLOUT(8)) .and. (.not. LPDLIB)) THEN + if (use_restartnc) then + if (restart_from_binary) then + do_rstsetup = .true. + else + do_rstsetup = .false. + end if + else + do_rstsetup = .true. + end if + ! + IF ((FLOUT(4) .OR. FLOUT(8)) .and. (.not. LPDLIB) .and. do_rstsetup) THEN IF (OARST) THEN ALLOCATE ( OUTPTS(IMOD)%OUT4%IRQRS(34*NAPROC) ) ELSE @@ -5647,7 +5688,7 @@ SUBROUTINE W3MPIO ( IMOD ) ! END IF ! - END IF ! IF ((FLOUT(4) .OR. FLOUT(8)) .and. (.not. LPDLIB)) THEN + END IF ! IF ((FLOUT(4) .OR. FLOUT(8)) .and. (.not. LPDLIB) .and. do_rstsetup) THEN #endif ! ! 3. Set-up for W3IOBC ( SENDs ) ------------------------------------ / diff --git a/model/src/w3iogomd.F90 b/model/src/w3iogomd.F90 index 61495a6fe..03810ca27 100644 --- a/model/src/w3iogomd.F90 +++ b/model/src/w3iogomd.F90 @@ -2570,8 +2570,6 @@ SUBROUTINE W3IOGO ( INXOUT, NDSOG, IOTST, IMOD ) #ifdef W3_S USE W3SERVMD, ONLY: STRACE #endif - use w3timemd , only: set_user_timestring - use w3odatmd , only: use_user_histname, user_histfname ! !/ !/ ------------------------------------------------------------------- / @@ -2600,8 +2598,6 @@ SUBROUTINE W3IOGO ( INXOUT, NDSOG, IOTST, IMOD ) #endif CHARACTER(LEN=30) :: IDTST, TNAME CHARACTER(LEN=10) :: VERTST - CHARACTER(len=512) :: FNAME - character(len=16) :: user_timestring !YYYY-MM-DD-SSSSS !/ !/ ------------------------------------------------------------------- / !/ @@ -2653,25 +2649,15 @@ SUBROUTINE W3IOGO ( INXOUT, NDSOG, IOTST, IMOD ) IF ( IPASS.EQ.1 .AND. OFILES(1) .EQ. 0) THEN I = LEN_TRIM(FILEXT) J = LEN_TRIM(FNMPRE) - if (use_user_histname) then - if (len_trim(user_histfname) == 0 ) then - call extcde (60, MSG="user history filename requested"// & - " but not provided") - end if - call set_user_timestring(time,user_timestring) - fname = trim(user_histfname)//trim(user_timestring) - else - fname = 'out_grd.'//FILEXT(:I) - end if ! #ifdef W3_T - WRITE (NDST,9001) FNMPRE(:J)//trim(fname) + WRITE (NDST,9001) FNMPRE(:J)//'out_grd.'//FILEXT(:I) #endif IF ( WRITE ) THEN - OPEN (NDSOG,FILE=FNMPRE(:J)//trim(fname), & + OPEN (NDSOG,FILE=FNMPRE(:J)//'out_grd.'//FILEXT(:I), & form='UNFORMATTED', convert=file_endian,ERR=800,IOSTAT=IERR) ELSE - OPEN (NDSOG,FILE=FNMPRE(:J)//trim(fname), & + OPEN (NDSOG,FILE=FNMPRE(:J)//'out_grd.'//FILEXT(:I), & form='UNFORMATTED', convert=file_endian,ERR=800,IOSTAT=IERR,STATUS='OLD') END IF ! @@ -2730,30 +2716,19 @@ SUBROUTINE W3IOGO ( INXOUT, NDSOG, IOTST, IMOD ) CALL EXTCDE ( 2 ) END IF END IF - ! IF ( IPASS.GE.1 .AND. OFILES(1) .EQ. 1) THEN I = LEN_TRIM(FILEXT) J = LEN_TRIM(FNMPRE) - if (use_user_histname) then - if (len_trim(user_histfname) == 0 ) then - call extcde (60, MSG="user history filename requested"// & - " but not provided") - end if - call set_user_timestring(time,user_timestring) - fname = trim(user_histfname)//trim(user_timestring) - else - ! - ! Create TIMETAG for file name using YYYYMMDD.HHMMS prefix - WRITE(TIMETAG,"(i8.8,'.'i6.6)")TIME(1),TIME(2) + ! + ! Create TIMETAG for file name using YYYYMMDD.HHMMS prefix + WRITE(TIMETAG,"(i8.8,'.'i6.6)")TIME(1),TIME(2) #ifdef W3_T - WRITE (NDST,9001) FNMPRE(:J)//TIMETAG//'.out_grd.'//FILEXT(:I) + WRITE (NDST,9001) FNMPRE(:J)//TIMETAG//'.out_grd.'//FILEXT(:I) #endif - fname = TIMETAG//'.out_grd.'//FILEXT(:I) - end if IF ( WRITE ) THEN - OPEN (NDSOG,FILE=FNMPRE(:J)//trim(fname), & - form='UNFORMATTED', convert=file_endian,ERR=800,IOSTAT=IERR) + OPEN (NDSOG,FILE=FNMPRE(:J)//TIMETAG//'.out_grd.' & + //FILEXT(:I),form='UNFORMATTED', convert=file_endian,ERR=800,IOSTAT=IERR) ELSE OPEN (NDSOG,FILE=FNMPRE(:J)//'out_grd.'//FILEXT(:I), & form='UNFORMATTED', convert=file_endian,ERR=800,IOSTAT=IERR,STATUS='OLD') diff --git a/model/src/w3iogoncdmd.F90 b/model/src/w3iogoncdmd.F90 deleted file mode 100644 index 813aa28d2..000000000 --- a/model/src/w3iogoncdmd.F90 +++ /dev/null @@ -1,538 +0,0 @@ -!> @file w3iogoncmd -!! -!> @brief Write gridded model output as netCDF -!! -!> @author mvertens@ucar.edu, Denise.Worthen@noaa.gov -!> @date 01-05-2022 -#include "w3macros.h" - -module w3iogoncdmd - - use w3gdatmd , only : nk, nx, ny, mapsf, mapsta, nsea - use w3odatmd , only : noswll, undef - use w3odatmd , only : nds, iaproc, napout - use netcdf - - implicit none - - private - - public :: w3iogoncd - - ! used/reused in module - - integer :: isea, ierr, ncid, varid - integer :: len_s, len_m, len_p, len_k - character(len=1024) :: fname - - real, allocatable, target :: var3ds(:,:,:) - real, allocatable, target :: var3dm(:,:,:) - real, allocatable, target :: var3dp(:,:,:) - real, allocatable, target :: var3dk(:,:,:) - - real, pointer :: var3d(:,:,:) - - !=============================================================================== -contains - !=============================================================================== - - subroutine w3iogoncd () - - use w3odatmd , only : fnmpre - use w3gdatmd , only : filext, trigp, ntri, ungtype, gtype - use w3servmd , only : extcde - use w3wdatmd , only : w3setw, w3dimw, time, wlv, ice, icef, iceh, berg, ust, ustdir, asf, rhoair - use w3gdatmd , only : xgrd, ygrd - use w3gdatmd , only : e3df, p2msf, us3df, usspf, w3setg - use w3odatmd , only : nogrp, ngrpp, idout, ndst, ndse, noswll, w3seto - use w3adatmd , only : w3seta, w3dima, w3xeta - use w3adatmd , only : ainit, dw, ua, ud, as, cx, cy, wn, taua, tauadir - use w3adatmd , only : hs, wlm, t02, t0m1, t01, fp0, thm, ths, thp0, wbt, wnmean - use w3adatmd , only : dtdyn - use w3adatmd , only : fcut, aba, abd, uba, ubd, sxx, syy, sxy - use w3adatmd , only : phs, ptp, plp, pdir, psi, pws, pwst, pnr - use w3adatmd , only : pthp0, pqp, ppe, pgw, psw, ptm1, pt1, pt2 - use w3adatmd , only : pep, usero, tauox, tauoy, tauwix, tauwiy - use w3adatmd , only : phiaw, phioc, tusx, tusy, prms, tpms - use w3adatmd , only : ussx, ussy, mssx, mssy, mssd, mscx, mscy - use w3adatmd , only : mscd, qp, tauwnx, tauwny, charn, tws, bhd - use w3adatmd , only : phibbl, taubbl, whitecap, bedforms, cge, ef - use w3adatmd , only : cflxymax, cflthmax, cflkmax, p2sms, us3d - use w3adatmd , only : th1m, sth1m, th2m, sth2m, hsig, phice, tauice - use w3adatmd , only : stmaxe, stmaxd, hmaxe, hcmaxe, hmaxd, hcmaxd, ussp, tauocx, tauocy - use w3adatmd , only : usshx, usshy - use wav_grdout , only : varatts, outvars - use w3timemd , only : set_user_timestring - use w3odatmd , only : time_origin, calendar_name, elapsed_secs - use w3odatmd , only : use_user_histname, user_histfname - !TODO: use unstr_mesh from wav_shr_mod; currently fails due to CI - !use wav_shr_mod, only : unstr_mesh - - ! local variables - integer :: igrd - integer ,target :: dimid3(3) - integer ,target :: dimid4(4) - integer ,pointer :: dimid(:) - character(len=12) :: vname - character(len=16) :: user_timestring !YYYY-MM-DD-SSSSS - - integer :: n, xtid, ytid, xeid, ztid, stid, mtid, ptid, ktid, timid, varid - logical :: s_axis = .false., m_axis = .false., p_axis = .false., k_axis = .false. - - !------------------------------------------------------------------------------- - - igrd = 1 - call w3seto ( igrd, ndse, ndst ) - call w3setg ( igrd, ndse, ndst ) - call w3seta ( igrd, ndse, ndst ) ! sets pointers into wadats in w3adatmd - call w3xeta ( igrd, ndse, ndst ) ! sets pointers into wadats in w3adatmd - call w3setw ( igrd, ndse, ndst ) ! sets pointers into wdatas in w3wdatmd - - ! ------------------------------------------------------------- - ! create the netcdf file - ! ------------------------------------------------------------- - - if (use_user_histname) then - if (len_trim(user_histfname) == 0 ) then - call extcde (60, msg="user history filename requested but not provided") - end if - call set_user_timestring(time,user_timestring) - fname = trim(user_histfname)//trim(user_timestring)//'.nc' - else - write(fname,'(a,i8.8,a1,i6.6,a)')trim(fnmpre),time(1),'.',time(2),'.out_grd.'//trim(filext)//'.nc' - end if - - len_s = noswll + 1 ! 0:noswll - len_m = p2msf(3)-p2msf(2) + 1 ! ? - len_p = usspf(2) ! partitions - len_k = e3df(3,1) - e3df(2,1) + 1 ! frequencies - - ! define the dimensions required for the requested gridded fields - do n = 1,size(outvars) - if (outvars(n)%validout) then - if(trim(outvars(n)%dims) == 's')s_axis = .true. - 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 - end do - - ! allocate arrays if needed - if (s_axis) allocate(var3ds(1:nx,1:ny,len_s)) - if (m_axis) allocate(var3dm(1:nx,1:ny,len_m)) - if (p_axis) allocate(var3dp(1:nx,1:ny,len_p)) - if (k_axis) allocate(var3dk(1:nx,1:ny,len_k)) - - ! create the netcdf file - ierr = nf90_create(trim(fname), nf90_clobber, ncid) - call handle_err(ierr, 'nf90_create') - ierr = nf90_def_dim(ncid, 'nx', nx, xtid) - ierr = nf90_def_dim(ncid, 'ny', ny, ytid) - ierr = nf90_def_dim(ncid, 'time', nf90_unlimited, timid) - - if (s_axis) ierr = nf90_def_dim(ncid, 'noswll', len_s, stid) - if (m_axis) ierr = nf90_def_dim(ncid, 'nm' , len_m, mtid) - if (p_axis) ierr = nf90_def_dim(ncid, 'np' , len_p, ptid) - if (k_axis) ierr = nf90_def_dim(ncid, 'freq' , len_k, ktid) - if (gtype .eq. ungtype) then - ierr = nf90_def_dim(ncid, 'ne' , ntri, xeid) - ierr = nf90_def_dim(ncid, 'nn' , 3, ztid) - end if - - ! define the time variable - ierr = nf90_def_var(ncid, 'time', nf90_double, timid, varid) - call handle_err(ierr,'def_timevar') - ierr = nf90_put_att(ncid, varid, 'units', trim(time_origin)) - call handle_err(ierr,'def_time_units') - ierr = nf90_put_att(ncid, varid, 'calendar', trim(calendar_name)) - call handle_err(ierr,'def_time_calendar') - - ! define the spatial axis variables (lat,lon) - ierr = nf90_def_var(ncid, 'lon', nf90_double, (/xtid,ytid/), varid) - call handle_err(ierr,'def_lonvar') - ierr = nf90_put_att(ncid, varid, 'units', 'degrees_east') - ierr = nf90_def_var(ncid, 'lat', nf90_double, (/xtid,ytid/), varid) - call handle_err(ierr,'def_latvar') - ierr = nf90_put_att(ncid, varid, 'units', 'degrees_north') - - ! add mapsta - ierr = nf90_def_var(ncid, 'mapsta', nf90_int, (/xtid, ytid, timid/), varid) - call handle_err(ierr, 'def_mapsta') - ierr = nf90_put_att(ncid, varid, 'units', 'unitless') - ierr = nf90_put_att(ncid, varid, 'long_name', 'map status') - - if (gtype .eq. ungtype) then - ierr = nf90_def_var(ncid, 'nconn', nf90_int, (/ztid,xeid/), varid) - call handle_err(ierr,'def_nodeconnections') - ierr = nf90_put_att(ncid, varid, 'units', 'unitless') - ierr = nf90_put_att(ncid, varid, 'long_name', 'node connectivity') - end if - - ! define the variables - dimid3(1:2) = (/xtid, ytid/) - dimid4(1:2) = (/xtid, ytid/) - do n = 1,size(outvars) - if (trim(outvars(n)%dims) == 's') then - dimid4(3:4) = (/stid, timid/) - dimid => dimid4 - else if (trim(outvars(n)%dims) == 'm') then - dimid4(3:4) = (/mtid, timid/) - dimid => dimid4 - else if (trim(outvars(n)%dims) == 'p') then - dimid4(3:4) = (/ptid, timid/) - dimid => dimid4 - else if (trim(outvars(n)%dims) == 'k') then - dimid4(3:4) = (/ktid, timid/) - dimid => dimid4 - else - dimid3(3) = timid - dimid => dimid3 - end if - - ierr = nf90_def_var(ncid, trim(outvars(n)%var_name), nf90_float, dimid, varid) - call handle_err(ierr, 'define variable '//trim((outvars(n)%var_name))) - ierr = nf90_put_att(ncid, varid, 'units' , trim(outvars(n)%unit_name)) - ierr = nf90_put_att(ncid, varid, 'long_name' , trim(outvars(n)%long_name)) - ierr = nf90_put_att(ncid, varid, '_FillValue', undef) - end do - ! end variable definitions - ierr = nf90_enddef(ncid) - call handle_err(ierr, 'end variable definition') - - ! write the time and spatial axis values (lat,lon,time) - ierr = nf90_inq_varid(ncid, 'lat', varid) - call handle_err(ierr, 'inquire variable lat ') - ierr = nf90_put_var(ncid, varid, transpose(ygrd)) - call handle_err(ierr, 'put lat') - - ierr = nf90_inq_varid(ncid, 'lon', varid) - call handle_err(ierr, 'inquire variable lon ') - ierr = nf90_put_var(ncid, varid, transpose(xgrd)) - call handle_err(ierr, 'put lon') - - ierr = nf90_inq_varid(ncid, 'time', varid) - call handle_err(ierr, 'inquire variable time ') - ierr = nf90_put_var(ncid, varid, elapsed_secs) - call handle_err(ierr, 'put time') - - if (gtype .eq. ungtype) then - ierr = nf90_inq_varid(ncid, 'nconn', varid) - call handle_err(ierr, 'inquire variable nconn ') - ierr = nf90_put_var(ncid, varid, trigp) - call handle_err(ierr, 'put trigp') - end if - - !maps - ierr = nf90_inq_varid(ncid, 'mapsta', varid) - call handle_err(ierr, 'inquire variable mapsta ') - ierr = nf90_put_var(ncid, varid, transpose(mapsta)) - call handle_err(ierr, 'put mapsta') - - ! close the file - ierr = nf90_close(ncid) - - ! write the requested variables - do n = 1,size(outvars) - vname = trim(outvars(n)%var_name) - if (trim(outvars(n)%dims) == 's') then - var3d => var3ds - ! Group 4 - if(vname .eq. 'PHS') call write_var3d(vname, phs (1:nsea,0:noswll) ) - if(vname .eq. 'PTP') call write_var3d(vname, ptp (1:nsea,0:noswll) ) - if(vname .eq. 'PLP') call write_var3d(vname, plp (1:nsea,0:noswll) ) - if(vname .eq. 'PDIR') call write_var3d(vname, pdir (1:nsea,0:noswll) ) - if(vname .eq. 'PSI') call write_var3d(vname, psi (1:nsea,0:noswll) ) - if(vname .eq. 'PWS') call write_var3d(vname, pws (1:nsea,0:noswll) ) - if(vname .eq. 'PDP') call write_var3d(vname, pthp0 (1:nsea,0:noswll) ) - if(vname .eq. 'PQP') call write_var3d(vname, pqp (1:nsea,0:noswll) ) - if(vname .eq. 'PPE') call write_var3d(vname, ppe (1:nsea,0:noswll) ) - if(vname .eq. 'PGW') call write_var3d(vname, pgw (1:nsea,0:noswll) ) - if(vname .eq. 'PSW') call write_var3d(vname, psw (1:nsea,0:noswll) ) - if(vname .eq. 'PTM1') call write_var3d(vname, ptm1 (1:nsea,0:noswll) ) - if(vname .eq. 'PT1') call write_var3d(vname, pt1 (1:nsea,0:noswll) ) - if(vname .eq. 'PT2') call write_var3d(vname, pt2 (1:nsea,0:noswll) ) - if(vname .eq. 'PEP') call write_var3d(vname, pep (1:nsea,0:noswll) ) - - else if (trim(outvars(n)%dims) == 'm') then ! m axis - var3d => var3dm - ! Group 6 - if (vname .eq. 'P2SMS') call write_var3d(vname, p2sms (1:nsea,p2msf(2):p2msf(3)) ) - - else if (trim(outvars(n)%dims) == 'p') then ! partition axis - var3d => var3dp - ! Group 6 - if (vname .eq. 'USSPX') call write_var3d(vname, ussp (1:nsea, 1:usspf(2)) ) - if (vname .eq. 'USSPY') call write_var3d(vname, ussp (1:nsea,nk+1:nk+usspf(2)) ) - - else if (trim(outvars(n)%dims) == 'k') then ! freq axis - var3d => var3dk - ! Group 3 - if(vname .eq. 'EF') call write_var3d(vname, ef (1:nsea,e3df(2,1):e3df(3,1)) ) - if(vname .eq. 'TH1M') call write_var3d(vname, ef (1:nsea,e3df(2,2):e3df(3,2)) ) - if(vname .eq. 'STH1M') call write_var3d(vname, ef (1:nsea,e3df(2,3):e3df(3,3)) ) - if(vname .eq. 'TH2M') call write_var3d(vname, ef (1:nsea,e3df(2,4):e3df(3,4)) ) - if(vname .eq. 'STH2M') call write_var3d(vname, ef (1:nsea,e3df(2,5):e3df(3,5)) ) - !TODO: wn has reversed indices (1:nk, 1:nsea) - ! Group 6 - if (vname .eq. 'US3DX') call write_var3d(vname, us3d (1:nsea, us3df(2):us3df(3)) ) - if (vname .eq. 'US3DY') call write_var3d(vname, us3d (1:nsea,nk+us3df(2):nk+us3df(3)) ) - - else - ! Group 1 - if (vname .eq. 'DW') call write_var2d(vname, dw (1:nsea), init0='false') - if (vname .eq. 'CX') call write_var2d(vname, cx (1:nsea), init0='false') - if (vname .eq. 'CY') call write_var2d(vname, cy (1:nsea), init0='false') - if (vname .eq. 'UAX') call write_var2d(vname, ua (1:nsea), dir=cos(ud(1:nsea)), init0='false') - if (vname .eq. 'UAY') call write_var2d(vname, ua (1:nsea), dir=sin(ud(1:nsea)), init0='false') - if (vname .eq. 'AS') call write_var2d(vname, as (1:nsea), init0='false') - if (vname .eq. 'WLV') call write_var2d(vname, wlv (1:nsea), init0='false') - if (vname .eq. 'ICE') call write_var2d(vname, ice (1:nsea), init0='false') - if (vname .eq. 'BERG') call write_var2d(vname, berg (1:nsea), init0='false') - if (vname .eq. 'TAUX') call write_var2d(vname, taua (1:nsea), dir=cos(tauadir(1:nsea)), init0='false') - if (vname .eq. 'TAUY') call write_var2d(vname, taua (1:nsea), dir=sin(tauadir(1:nsea)), init0='false') - if (vname .eq. 'RHOAIR') call write_var2d(vname, rhoair (1:nsea), init0='false') - if (vname .eq. 'ICEH') call write_var2d(vname, iceh (1:nsea), init0='false') - if (vname .eq. 'ICEF') call write_var2d(vname, icef (1:nsea), init0='false') - - ! Group 2 - if (vname .eq. 'HS') call write_var2d(vname, hs (1:nsea) ) - if (vname .eq. 'WLM') call write_var2d(vname, wlm (1:nsea) ) - if (vname .eq. 'T02') call write_var2d(vname, t02 (1:nsea) ) - if (vname .eq. 'T0M1') call write_var2d(vname, t0m1 (1:nsea) ) - if (vname .eq. 'T01') call write_var2d(vname, t01 (1:nsea) ) - if (vname .eq. 'FP0') call write_var2d(vname, fp0 (1:nsea) ) - if (vname .eq. 'THM') call write_var2d(vname, thm (1:nsea) ) - if (vname .eq. 'THS') call write_var2d(vname, ths (1:nsea) ) - if (vname .eq. 'THP0') call write_var2d(vname, thp0 (1:nsea) ) - if (vname .eq. 'HSIG') call write_var2d(vname, hsig (1:nsea) ) - if (vname .eq. 'STMAXE') call write_var2d(vname, stmaxe (1:nsea) ) - if (vname .eq. 'STMAXD') call write_var2d(vname, stmaxd (1:nsea) ) - if (vname .eq. 'HMAXE') call write_var2d(vname, hmaxe (1:nsea) ) - if (vname .eq. 'HCMAXE') call write_var2d(vname, hcmaxe (1:nsea) ) - if (vname .eq. 'HMAXD') call write_var2d(vname, hmaxd (1:nsea) ) - if (vname .eq. 'HCMAXD') call write_var2d(vname, hcmaxd (1:nsea) ) - if (vname .eq. 'WBT') call write_var2d(vname, wbt (1:nsea) ) - if (vname .eq. 'WNMEAN') call write_var2d(vname, wnmean (1:nsea), init0='false') - - ! Group 4 - if(vname .eq. 'PWST') call write_var2d(vname, pwst (1:nsea) ) - if(vname .eq. 'PNR') call write_var2d(vname, pnr (1:nsea) ) - - ! Group 5 - if (vname .eq. 'USTX') call write_var2d(vname, ust (1:nsea)*asf(1:nsea), dir=cos(ustdir(1:nsea)), usemask='true') - if (vname .eq. 'USTY') call write_var2d(vname, ust (1:nsea)*asf(1:nsea), dir=sin(ustdir(1:nsea)), usemask='true') - if (vname .eq. 'CHARN') call write_var2d(vname, charn (1:nsea) ) - if (vname .eq. 'CGE') call write_var2d(vname, cge (1:nsea) ) - if (vname .eq. 'PHIAW') call write_var2d(vname, phiaw (1:nsea), init2='true') - if (vname .eq. 'TAUWIX') call write_var2d(vname, tauwix (1:nsea), init2='true') - if (vname .eq. 'TAUWIY') call write_var2d(vname, tauwiy (1:nsea), init2='true') - if (vname .eq. 'TAUWNX') call write_var2d(vname, tauwnx (1:nsea), init2='true') - if (vname .eq. 'TAUWNY') call write_var2d(vname, tauwny (1:nsea), init2='true') - if (vname .eq. 'WCC') call write_var2d(vname, whitecap (1:nsea,1), init2='true') - if (vname .eq. 'WCF') call write_var2d(vname, whitecap (1:nsea,2), init2='true') - if (vname .eq. 'WCH') call write_var2d(vname, whitecap (1:nsea,3), init2='true') - if (vname .eq. 'WCM') call write_var2d(vname, whitecap (1:nsea,4), init2='true') - if (vname .eq. 'TWS') call write_var2d(vname, tws (1:nsea) ) - - ! Group 6 - if (vname .eq. 'SXX') call write_var2d(vname, sxx (1:nsea) ) - if (vname .eq. 'SYY') call write_var2d(vname, syy (1:nsea) ) - if (vname .eq. 'SXY') call write_var2d(vname, sxy (1:nsea) ) - if (vname .eq. 'TAUOX') call write_var2d(vname, tauox (1:nsea), init2='true') - if (vname .eq. 'TAUOY') call write_var2d(vname, tauoy (1:nsea), init2='true') - if (vname .eq. 'BHD') call write_var2d(vname, bhd (1:nsea) ) - if (vname .eq. 'PHIOC') call write_var2d(vname, phioc (1:nsea), init2='true') - if (vname .eq. 'TUSX') call write_var2d(vname, tusx (1:nsea) ) - if (vname .eq. 'TUSY') call write_var2d(vname, tusy (1:nsea) ) - if (vname .eq. 'USSX') call write_var2d(vname, ussx (1:nsea) ) - if (vname .eq. 'USSY') call write_var2d(vname, ussy (1:nsea) ) - if (vname .eq. 'PRMS') call write_var2d(vname, prms (1:nsea) ) - if (vname .eq. 'TPMS') call write_var2d(vname, tpms (1:nsea) ) - if (vname .eq. 'TAUICEX') call write_var2d(vname, tauice (1:nsea,1) ) - if (vname .eq. 'TAUICEY') call write_var2d(vname, tauice (1:nsea,2) ) - if (vname .eq. 'PHICE') call write_var2d(vname, phice (1:nsea) ) - if (vname .eq. 'TAUOCX') call write_var2d(vname, tauocx (1:nsea) ) - if (vname .eq. 'TAUOCY') call write_var2d(vname, tauocy (1:nsea) ) - if (vname .eq. 'USSHX') call write_var2d(vname, usshx (1:nsea) ) - if (vname .eq. 'USSHY') call write_var2d(vname, usshy (1:nsea) ) - ! Group 7 - if (vname .eq. 'ABAX') call write_var2d(vname, aba (1:nsea), dir=cos(abd(1:nsea)) ) - if (vname .eq. 'ABAY') call write_var2d(vname, aba (1:nsea), dir=sin(abd(1:nsea)) ) - if (vname .eq. 'UBAX') call write_var2d(vname, uba (1:nsea), dir=cos(ubd(1:nsea)) ) - if (vname .eq. 'UBAY') call write_var2d(vname, uba (1:nsea), dir=sin(ubd(1:nsea)) ) - if (vname .eq. 'BED') call write_var2d(vname, bedforms (1:nsea,1), init2='true') - if (vname .eq. 'RIPPLEX') call write_var2d(vname, bedforms (1:nsea,2), init2='true') - if (vname .eq. 'RIPPLEY') call write_var2d(vname, bedforms (1:nsea,3), init2='true') - if (vname .eq. 'PHIBBL') call write_var2d(vname, phibbl (1:nsea), init2='true') - if (vname .eq. 'TAUBBLX') call write_var2d(vname, taubbl (1:nsea,1), init2='true') - if (vname .eq. 'TAUBBLY') call write_var2d(vname, taubbl (1:nsea,2), init2='true') - - ! Group 8 - if (vname .eq. 'MSSX') call write_var2d(vname, mssx (1:nsea) ) - if (vname .eq. 'MSSY') call write_var2d(vname, mssy (1:nsea) ) - if (vname .eq. 'MSCX') call write_var2d(vname, mscx (1:nsea) ) - if (vname .eq. 'MSCY') call write_var2d(vname, mscy (1:nsea) ) - !TODO: remaining variables have inconsistency between shel_inp listing and iogo code - - ! Group 9 - if (vname .eq. 'DTDYN') call write_var2d(vname, dtdyn (1:nsea) ) - if (vname .eq. 'FCUT') call write_var2d(vname, fcut (1:nsea) ) - if (vname .eq.'CFLXYMAX') call write_var2d(vname, cflxymax (1:nsea) ) - if (vname .eq.'CFLTHMAX') call write_var2d(vname, cflthmax (1:nsea) ) - if (vname .eq. 'CFLKMAX') call write_var2d(vname, cflkmax (1:nsea) ) - - ! Group 10 - end if - end do - - if (s_axis) deallocate(var3ds) - if (m_axis) deallocate(var3dm) - if (p_axis) deallocate(var3dp) - if (k_axis) deallocate(var3dk) - - ! Flush the buffers for write - call w3seta ( igrd, ndse, ndst ) - - end subroutine w3iogoncd - - !=============================================================================== - subroutine write_var2d(vname, var, dir, usemask, init0, init2) - ! write (nsea) array as (nx,ny) - ! if dir is present, write x or y component of (nsea) array as (nx,ny) - ! if mask is present and true, use mapsta=1 to mask values - ! if init0 is present and false, do not initialize values - ! for mapsta<0. this prevents group 1 variables being set undef over - ! ice. if init2 is present and true, apply a second initialization to - ! a subset of variables for where mapsta==2 - - character(len=*), intent(in) :: vname - real , intent(in) :: var(:) - real, optional , intent(in) :: dir(:) - character(len=*), optional, intent(in) :: usemask - character(len=*), optional, intent(in) :: init0 - character(len=*), optional, intent(in) :: init2 - - ! local variables - real, dimension(nx,ny) :: var2d - logical :: lmask, linit0, linit2 - real :: varloc - - lmask = .false. - if (present(usemask)) then - lmask = (trim(usemask) == "true") - end if - linit0 = .true. - if (present(init0)) then - linit0 = (trim(init0) == "true") - end if - linit2 = .false. - if (present(init2)) then - linit2 = (trim(init2) == "true") - end if - - ! DEBUG - ! write(nds(1),'(a)')' writing variable ' //trim(vname)//' to history file '//trim(fname) - - var2d = undef - do isea = 1,nsea - - ! initialization - varloc = var(isea) - if (linit0) then - if (mapsta(mapsf(isea,2),mapsf(isea,1)) < 0) varloc = undef - end if - if (linit2) then - if (mapsta(mapsf(isea,2),mapsf(isea,1)) == 2) varloc = undef - end if - - if (present(dir)) then - if (varloc .ne. undef) then - if (lmask) then - if (mapsta(mapsf(isea,2),mapsf(isea,1)) == 1) then - var2d(mapsf(isea,1),mapsf(isea,2)) = varloc*dir(isea) - end if - else - var2d(mapsf(isea,1),mapsf(isea,2)) = varloc*dir(isea) - end if - end if - else - var2d(mapsf(isea,1),mapsf(isea,2)) = varloc - end if - end do - - ierr = nf90_open(trim(fname), nf90_write, ncid) - call handle_err(ierr, 'open '//trim(fname)//' for writing') - ierr = nf90_inq_varid(ncid, trim(vname), varid) - call handle_err(ierr, 'inquire variable '//trim(vname)) - ierr = nf90_put_var(ncid, varid, var2d) - call handle_err(ierr, 'put variable '//trim(vname)) - ierr = nf90_close(ncid) - - end subroutine write_var2d - - !=============================================================================== - subroutine write_var3d(vname, var, init2) - ! write (nsea,:) array as (nx,ny,:) - ! if init2 is present and true, apply a second initialization to - ! a subset of variables for where mapsta==2 - - character(len=*), intent(in) :: vname - real , intent(in) :: var(:,:) - character(len=*), optional, intent(in) :: init2 - - ! local variables - real, allocatable, dimension(:) :: varloc - logical :: linit2 - integer :: lb, ub - - linit2 = .false. - if (present(init2)) then - linit2 = (trim(init2) == "true") - end if - - lb = lbound(var,2) - ub = ubound(var,2) - allocate(varloc(lb:ub)) - - ! DEBUG - ! write(nds(1),'(a,2i6)')' writing variable ' //trim(vname)//' to history file ' & - ! //trim(fname)//' with bounds ',lb,ub - - var3d = undef - do isea = 1,nsea - ! initialization - varloc(:) = var(isea,:) - 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 - var3d(mapsf(isea,1),mapsf(isea,2),:) = varloc(:) - end do - - ierr = nf90_open(trim(fname), nf90_write, ncid) - call handle_err(ierr, 'open '//trim(fname)//' for writing') - ierr = nf90_inq_varid(ncid, trim(vname), varid) - call handle_err(ierr, 'inquire variable '//trim(vname)) - ierr = nf90_put_var(ncid, varid, var3d) - call handle_err(ierr, 'put variable '//trim(vname)) - ierr = nf90_close(ncid) - - deallocate(varloc) - end subroutine write_var3d - - !=============================================================================== - subroutine handle_err(ierr,string) - use w3odatmd , only : ndse - use w3servmd , only : extcde - - ! input/output variables - integer , intent(in) :: ierr - character(len=*), intent(in) :: string - - if (ierr /= nf90_noerr) then - write(ndse,*) "*** WAVEWATCH III netcdf error: ",trim(string),':',trim(nf90_strerror(ierr)) - call extcde ( 49 ) - end if - end subroutine handle_err - -end module w3iogoncdmd diff --git a/model/src/w3iorsmd.F90 b/model/src/w3iorsmd.F90 index 24d9a280c..f983903eb 100644 --- a/model/src/w3iorsmd.F90 +++ b/model/src/w3iorsmd.F90 @@ -111,7 +111,7 @@ MODULE W3IORSMD !> !> @author H. L. Tolman @date 22-Mar-2021 !> - SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT ) + SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT , filename) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | @@ -327,7 +327,7 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT ) #ifdef W3_TIMINGS USE W3PARALL, ONLY: PRINT_MY_TIME #endif - USE w3odatmd, ONLY : RUNTYPE, INITFILE + USE w3odatmd, ONLY : RUNTYPE USE w3adatmd, ONLY : USSHX, USSHY #ifdef W3_PDLIB USE PDLIB_FIELD_VEC @@ -336,9 +336,6 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT ) USE W3SERVMD, ONLY: STRACE #endif ! - use w3timemd, only: set_user_timestring - use w3odatmd, only: use_user_restname, user_restfname, ndso - ! #ifdef W3_MPI INCLUDE "mpif.h" #endif @@ -352,6 +349,7 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT ) REAL, INTENT(INOUT) :: DUMFPI CHARACTER, INTENT(IN) :: INXOUT*(*) LOGICAL, INTENT(IN),OPTIONAL :: FLRSTRT + character(len=*), intent(in), optional :: filename !/ !/ ------------------------------------------------------------------- / !/ Local parameters @@ -382,12 +380,10 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT ) LOGICAL :: NDSROPN CHARACTER(LEN=4) :: TYPE CHARACTER(LEN=10) :: VERTST - CHARACTER(LEN=512) :: FNAME + CHARACTER(LEN=40) :: FNAME CHARACTER(LEN=26) :: IDTST CHARACTER(LEN=30) :: TNAME CHARACTER(LEN=15) :: TIMETAG - character(len=16) :: user_timestring !YYYY-MM-DD-SSSSS - logical :: exists !/ !/ ------------------------------------------------------------------- / !/ @@ -465,46 +461,9 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT ) ! ! open file ---------------------------------------------------------- * ! - if (use_user_restname) then - ierr = -99 - if (.not. write) then - if (runtype == 'initial') then - if (len_trim(initfile) == 0) then - ! no IC file, use startup option - goto 800 - else - ! IC file exists - use it - fname = trim(initfile) - end if - else - call set_user_timestring(time,user_timestring) - fname = trim(user_restfname)//trim(user_timestring) - inquire( file=trim(fname), exist=exists) - if (.not. exists) then - call extcde (60, msg="required initial/restart file " // trim(fname) // " does not exist") - end if - end if - else - call set_user_timestring(time,user_timestring) - fname = trim(user_restfname)//trim(user_timestring) - end if - ! write out filename - if (iaproc == naprst) then - IF ( WRITE ) THEN - write (ndso,'(a)') 'WW3: writing restart file '//trim(fname) - else - write (ndso,'(a)') 'WW3: reading initial/restart file '//trim(fname) - end if - end if - if ( write ) then - IF ( .NOT.IOSFLG .OR. IAPROC.EQ.NAPRST ) & - open (ndsr,file=trim(fname), form='unformatted', convert=file_endian, & - ACCESS='STREAM',ERR=800,IOSTAT=IERR) - ELSE ! READ - open (ndsr, file=trim(fname), form='unformatted', convert=file_endian, & - ACCESS='STREAM',ERR=800,IOSTAT=IERR, & - STATUS='OLD',ACTION='READ') - END IF + if (present(filename)) then ! only when restart_nc and restart_from_binary=true + open (ndsr,file=trim(filename),form='unformatted', convert=file_endian, & + access='stream',err=800,iostat=ierr, status='old',action='read') else I = LEN_TRIM(FILEXT) J = LEN_TRIM(FNMPRE) @@ -530,10 +489,9 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT ) IFILE = IFILE + 1 ! #ifdef W3_T - WRITE (NDST,9001) trim(FNAME), LRECL + WRITE (NDST,9001) FNAME, LRECL #endif ! - IF(NDST.EQ.NDSR)THEN IF ( IAPROC .EQ. NAPERR ) & WRITE(NDSE,'(A,I8)')'UNIT NUMBERS OF RESTART FILE AND '& @@ -543,14 +501,14 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT ) IF ( WRITE ) THEN IF ( .NOT.IOSFLG .OR. IAPROC.EQ.NAPRST ) & - OPEN (NDSR,FILE=FNMPRE(:J)//trim(FNAME),form='UNFORMATTED', convert=file_endian, & + OPEN (NDSR,FILE=FNMPRE(:J)//FNAME,form='UNFORMATTED', convert=file_endian, & ACCESS='STREAM',ERR=800,IOSTAT=IERR) ELSE - OPEN (NDSR,FILE=FNMPRE(:J)//trim(FNAME),form='UNFORMATTED', convert=file_endian, & + OPEN (NDSR,FILE=FNMPRE(:J)//FNAME,form='UNFORMATTED', convert=file_endian, & ACCESS='STREAM',ERR=800,IOSTAT=IERR, & STATUS='OLD',ACTION='READ') END IF - end if + end if ! if (present(filename)) ! ! test info ---------------------------------------------------------- * ! @@ -638,21 +596,11 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT ) END IF ELSE READ (NDSR,POS=RPOS,ERR=802,IOSTAT=IERR) TTIME -#ifdef W3_CESMCOUPLED - if (runtype == 'branch' .or. runtype == 'continue') then - IF (TIME(1).NE.TTIME(1) .OR. TIME(2).NE.TTIME(2)) THEN - IF ( IAPROC .EQ. NAPERR ) & - WRITE (NDSE,906) TTIME, TIME - CALL EXTCDE ( 20 ) - END IF - end if -#else IF (TIME(1).NE.TTIME(1) .OR. TIME(2).NE.TTIME(2)) THEN IF ( IAPROC .EQ. NAPERR ) & WRITE (NDSE,906) TTIME, TIME CALL EXTCDE ( 20 ) END IF -#endif END IF ! #ifdef W3_T @@ -687,7 +635,7 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT ) ! Original non-server version writing of spectra ! IF ( .NOT.IOSFLG .OR. (NAPROC.EQ.1.AND.NAPRST.EQ.1) ) THEN -#ifdef W3_MPI +#ifdef W3_MPI DO JSEA=1, NSEAL CALL INIT_GET_ISEA(ISEA, JSEA) NREC = ISEA + 2 @@ -696,7 +644,7 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT ) WRITEBUFF(1:NSPEC) = VA(1:NSPEC,JSEA) WRITE (NDSR,POS=RPOS,ERR=803,IOSTAT=IERR) WRITEBUFF END DO -#else +#else DO JSEA=1, NSEA ISEA = JSEA NREC = ISEA + 2 @@ -705,7 +653,7 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT ) WRITEBUFF(1:NSPEC) = VA(1:NSPEC,JSEA) WRITE (NDSR,POS=RPOS,ERR=803,IOSTAT=IERR) WRITEBUFF END DO -#endif +#endif ! ! I/O server version writing of spectra ( !/MPI ) ! diff --git a/model/src/w3odatmd.F90 b/model/src/w3odatmd.F90 index 5c1c7d239..1f0d804ac 100644 --- a/model/src/w3odatmd.F90 +++ b/model/src/w3odatmd.F90 @@ -558,24 +558,36 @@ MODULE W3ODATMD LOGICAL, POINTER :: FLFORM, FLCOMB, O6INIT INTEGER, POINTER :: PTMETH ! C. Bunney; Partitioning method REAL, POINTER :: PTFCUT ! C. Bunney; Part. 5 freq cut - character(len=8) :: runtype = '' !< @public the run type (startup,branch,continue) - character(len=256) :: initfile = '' !< @public name of wave initial condition file - !! if runtype is startup or branch run, then initfile is used - logical :: use_user_histname = .false. !<@public logical flag for user set history filenames - logical :: use_user_restname = .false. !<@public logical flag for user set restart filenames - character(len=512) :: user_histfname = '' !<@public user history filename prefix, timestring - !! YYYY-MM-DD-SSSSS will be appended - character(len=512) :: user_restfname = '' !<@public user restart filename prefix, timestring - !! YYYY-MM-DD-SSSSS will be appended - logical :: histwr = .false. !<@public logical to trigger history write - !! if true => write history file (snapshot) - logical :: rstwr = .false. !<@public logical to trigger restart write - !! if true => write restart - logical :: user_netcdf_grdout = .false. !<@public logical flag to use netCDF for gridded - !! field output - character(len= 36) :: time_origin = '' !< @public the time_origin used for netCDF output - character(len= 36) :: calendar_name = '' !< @public the calendar used for netCDF output - integer(kind=8) :: elapsed_secs = 0 !< @public the time in seconds from the time_origin + + character(len=8) :: runtype = '' !< @public the run type (startup,branch,continue) + character(len=256) :: initfile = '' !< @public name of wave initial condition file + !! if runtype is startup or branch run, then initfile is used + character(len=512) :: user_histfname = '' !< @public user history filename prefix, timestring + !! YYYY-MM-DD-SSSSS will be appended + character(len=512) :: user_restfname = '' !< @public user restart filename prefix, timestring + !! YYYY-MM-DD-SSSSS will be appended + logical :: histwr = .false. !< @public logical to trigger history write + !! if true => write history file (snapshot) + logical :: rstwr = .false. !< @public logical to trigger restart write + !! if true => write restart + logical :: use_historync = .false. !< @public logical flag to use netCDF for gridded + !! field output + logical :: use_restartnc = .false. !< @public logical flag to read and write netCDF restarts + logical :: restart_from_binary = .false. !< @public logical flag for restarting from binary restart + ! when use_restartnc is true + logical :: logfile_is_assigned = .false. !< @public logical flag for assignment of nds(1) to specified + !! log file in mesh cap + logical :: verboselog = .true. !< @public logical flag to enable verbose WW3 native logging + logical :: addrstflds = .false. !< @public logical flag for additional restart fields + integer :: rstfldcnt = 0 !< @public the actual number of additional restart fields + character(len=10), dimension(10) :: rstfldlist = '' !< @public a list of additional fields for the restart file, + !! currently set to a maximum of 10. Additional restart fields + !! are required only when waves are in the slow loop and ice + !! is present. Note that waves should not be in the slow loop + !! if coupling to CICE is set + character(len=36) :: time_origin = '' !< @public the time_origin used for netCDF output + character(len=36) :: calendar_name = '' !< @public the calendar used for netCDF output + integer(kind=8) :: elapsed_secs = 0 !< @public the time in seconds from the time_origin !/ CONTAINS !/ ------------------------------------------------------------------- / diff --git a/model/src/w3triamd.F90 b/model/src/w3triamd.F90 index 9fac503b6..a97bed7e4 100644 --- a/model/src/w3triamd.F90 +++ b/model/src/w3triamd.F90 @@ -578,12 +578,18 @@ SUBROUTINE GET_BOUNDARY_STATUS(STATUS) !/ ! integer*2, intent(out) :: STATUS(NX) - INTEGER :: COLLECTED(NX), NEXTVERT(NX), PREVVERT(NX) + integer, allocatable :: collected(:), nextvert(:), prevvert(:) INTEGER :: ISFINISHED, INEXT, IPREV INTEGER :: IPNEXT, IPPREV, ZNEXT, IP, I, IE #ifdef W3_S CALL STRACE (IENT, 'VA_SETUP_IOBPD') #endif + allocate(collected(nx)) + allocate(nextvert(nx)) + allocate(prevvert(nx)) + nextvert = 0 + prevvert = 0 + STATUS(:) = 0 DO IE=1,NTRI DO I=1,3 @@ -650,6 +656,11 @@ SUBROUTINE GET_BOUNDARY_STATUS(STATUS) EXIT END IF END DO + + deallocate(collected) + deallocate(nextvert) + deallocate(prevvert) + END SUBROUTINE GET_BOUNDARY_STATUS !/ -------------------------------------------------------------------/ @@ -852,9 +863,10 @@ SUBROUTINE UG_GETOPENBOUNDARY(TMPSTA,ZBIN,ZLIM) !/ Local parameters !/ INTEGER :: IBC, IX - INTEGER :: MASK(NX) + integer, allocatable :: mask(:) INTEGER*2 :: STATUS(NX) ! + allocate(mask(nx)) MASK(:)=1 CALL SET_IOBP (MASK, STATUS) ! @@ -870,6 +882,8 @@ SUBROUTINE UG_GETOPENBOUNDARY(TMPSTA,ZBIN,ZLIM) IF ( (TMPSTA(1,IX).EQ.1) .AND. (STATUS(IX).EQ.0) .AND. (ZBIN(1,IX) .LT. ZLIM)) TMPSTA(1,IX) = 2 END IF END DO + + deallocate(mask) ! END SUBROUTINE UG_GETOPENBOUNDARY !/ ------------------------------------------------------------------- / @@ -964,14 +978,14 @@ SUBROUTINE SPATIAL_GRID I2 = TRIGP(2,K) I3 = TRIGP(3,K) -!AR: todo call this only for global grid +!AR: todo call this only for global grid CALL FIX_PERIODCITY(I1,I2,I3,XGRD,YGRD,PT) ! ! cross product of edge-vector (orientated anticlockwise) ! - TRIA(K) = REAL( (PT(2,2)-PT(1,2)) & - *(PT(1,1)-PT(3,1)) & - +(PT(3,2)-PT(1,2)) & + TRIA(K) = REAL( (PT(2,2)-PT(1,2)) & + *(PT(1,1)-PT(3,1)) & + +(PT(3,2)-PT(1,2)) & *(PT(2,1)-PT(1,1)) )*0.5 ! ! test on negative triangle area, which means that the orientiation is not as assumed to be anticw. @@ -1193,8 +1207,8 @@ SUBROUTINE COUNT(TRIGPTEMP) !/ ------------------------------------------------------------------- / !/ local parameter - INTEGER :: CONN(NX) - INTEGER :: COUNTER, IP, IE, I, J, N(3) + integer, allocatable :: conn(:) + INTEGER :: COUNTER, IP, IE, I, J, N(3) #ifdef W3_S INTEGER :: IENT = 0 #endif @@ -1203,7 +1217,7 @@ SUBROUTINE COUNT(TRIGPTEMP) #ifdef W3_S CALL STRACE (IENT, 'COUNT') #endif - + allocate(conn(nx)) COUNTRI=0 COUNTOT=0 CONN(:)= 0 @@ -1234,6 +1248,7 @@ SUBROUTINE COUNT(TRIGPTEMP) ENDDO COUNTOT=J + deallocate(conn) END SUBROUTINE COUNT !/---------------------------------------------------------------------------- @@ -1395,12 +1410,11 @@ SUBROUTINE AREA_SI(IMOD) INTEGER :: COUNTER,ifound,alreadyfound INTEGER :: I, J, K, II INTEGER :: IP, IE, POS, POS_I, POS_J, POS_K, IP_I, IP_J, IP_K - INTEGER :: I1, I2, I3, IP2, CHILF(NX) - INTEGER :: TMP(NX), CELLVERTEX(NX,COUNTRI,2) + INTEGER :: I1, I2, I3, IP2 INTEGER :: COUNT_MAX DOUBLE PRECISION :: TRIA03 INTEGER, ALLOCATABLE :: PTABLE(:,:) - + integer, allocatable :: cellvertex(:,:,:), tmp(:) #ifdef W3_S INTEGER :: IENT = 0 #endif @@ -1425,18 +1439,20 @@ SUBROUTINE AREA_SI(IMOD) SI(I2) = SI(I2) + TRIA03 SI(I3) = SI(I3) + TRIA03 ENDDO + allocate(cellvertex(nx,countri,2)) + allocate(tmp(nx)) CELLVERTEX(:,:,:) = 0 ! Stores for each node the Elementnumbers of the connected Elements ! and the Position of the Node in the Element Index - CHILF = 0 + tmp = 0 DO IE = 1, NTRI DO J=1,3 I = TRIGP(J,IE)!INE(J,IE) - CHILF(I) = CHILF(I)+1 - CELLVERTEX(I,CHILF(I),1) = IE - CELLVERTEX(I,CHILF(I),2) = J + TMP(I) = TMP(I)+1 + CELLVERTEX(I,TMP(I),1) = IE + CELLVERTEX(I,TMP(I),2) = J END DO ENDDO ! @@ -1454,6 +1470,7 @@ SUBROUTINE AREA_SI(IMOD) END DO INDEX_CELL(IP+1)=J+1 END DO + deallocate(cellvertex) IF (.NOT. FSNIMP) RETURN @@ -1573,6 +1590,7 @@ SUBROUTINE AREA_SI(IMOD) END DO END DO + deallocate(tmp) DEALLOCATE(PTABLE) END SUBROUTINE AREA_SI @@ -2105,9 +2123,11 @@ SUBROUTINE UG_GRADIENTS (PARAM, DIFFX, DIFFY) REAL :: DIFFXTMP, DIFFYTMP REAL :: DEDX(3), DEDY(3) REAL :: DVDXIE, DVDYIE - REAL :: WEI(NX), WEI_LOCAL(NSEAL) + REAL :: WEI_LOCAL(NSEAL) + real, allocatable :: wei(:) REAL*8 :: RTMP(NSEAL) + allocate(wei(nx)) DIFFX = 0. DIFFY = 0. ! @@ -2166,6 +2186,7 @@ SUBROUTINE UG_GRADIENTS (PARAM, DIFFX, DIFFY) CALL PDLIB_exchange1Dreal(DIFFX(1,:)) CALL PDLIB_exchange1Dreal(DIFFY(1,:)) #endif + deallocate(wei) ! END SUBROUTINE UG_GRADIENTS !/ ------------------------------------------------------------------- / @@ -2382,14 +2403,21 @@ SUBROUTINE SET_IOBP (MASK, STATUS) INTEGER, INTENT(IN) :: MASK(NX) INTEGER*2, INTENT(OUT) :: STATUS(NX) ! - INTEGER :: COLLECTED(NX), NEXTVERT(NX), PREVVERT(NX) - INTEGER :: ISFINISHED !, INEXT, IPREV - INTEGER :: INEXT(3), IPREV(3) - INTEGER :: ZNEXT, IP, I, IE, IPNEXT, IPPREV, COUNT + integer, allocatable :: collected(:), nextvert(:), prevvert(:) + INTEGER :: ISFINISHED !, INEXT, IPREV + INTEGER :: INEXT(3), IPREV(3) + INTEGER :: ZNEXT, IP, I, IE, IPNEXT, IPPREV, COUNT integer nb0, nb1, nbM1 STATUS = -1 INEXT=(/ 2, 3, 1 /) !IPREV=1+MOD(I+1,3) IPREV=(/ 3, 1, 2 /) !INEXT=1+MOD(I,3) + + allocate(collected(nx)) + allocate(nextvert(nx)) + allocate(prevvert(nx)) + nextvert = 0 + prevvert = 0 + DO IE=1,NTRI ! If one of the points of the triangle is masked out (land) then do as if triangle does not exist... ! IF ((MASK(TRIGP(1,IE)).GT.0).AND.(MASK(TRIGP(2,IE)).GT.0).AND.(MASK(TRIGP(3,IE)).GT.0)) THEN @@ -2456,6 +2484,9 @@ SUBROUTINE SET_IOBP (MASK, STATUS) STATUS = 1 CALL GET_BOUNDARY(NX, NTRI, TRIGP, STATUS, PREVVERT, NEXTVERT) + deallocate(collected) + deallocate(nextvert) + deallocate(prevvert) !#ifdef MPI_PARALL_GRID ! CALL exchange_p2di(STATUS) !#endif @@ -2796,7 +2827,7 @@ SUBROUTINE TRIANG_INDEXES(I, INEXT, IPREV) END SUBROUTINE TRIANG_INDEXES !/ ------------------------------------------------------------------- / - + !> !> @brief Redefines the values of the boundary points and angle pointers !> based on the MAPSTA array. @@ -2903,7 +2934,7 @@ SUBROUTINE SET_UG_IOBP() REAL (KIND = 8) :: DYP1, DYP2, DYP3, eDet1, eDet2, EVX, EVY REAL(KIND=8), PARAMETER :: THR = TINY(1.) INTEGER :: I1, I2, I3 - INTEGER :: ITMP(NX), NEXTVERT(NX), PREVVERT(NX) + integer, allocatable :: itmp(:), nextvert(:), prevvert(:) CHARACTER(60) :: FNAME #ifdef W3_S INTEGER, SAVE :: IENT = 0 @@ -2916,6 +2947,11 @@ SUBROUTINE SET_UG_IOBP() #ifdef W3_S CALL STRACE (IENT, 'SETUGIOBP') #endif + allocate(itmp(nx)) + allocate(nextvert(nx)) + allocate(prevvert(nx)) + nextvert = 0 + prevvert = 0 ! !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! 2. Searches for boundary points @@ -3029,6 +3065,9 @@ SUBROUTINE SET_UG_IOBP() END IF END DO #endif + deallocate(itmp) + deallocate(nextvert) + deallocate(prevvert) ! ! Recomputes the angles used in the gradients estimation ! diff --git a/model/src/w3wavemd.F90 b/model/src/w3wavemd.F90 index b11aeb3b2..c6b5520d6 100644 --- a/model/src/w3wavemd.F90 +++ b/model/src/w3wavemd.F90 @@ -491,8 +491,11 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & #ifdef W3_TIMINGS USE W3PARALL, only : PRINT_MY_TIME #endif - use w3iogoncdmd , only : w3iogoncd - use w3odatmd , only : histwr, rstwr, user_netcdf_grdout + use wav_restart_mod , only : write_restart + use wav_history_mod , only : write_history + use w3odatmd , only : histwr, rstwr, use_historync, use_restartnc, user_restfname + use w3odatmd , only : verboselog + use w3timemd , only : set_user_timestring ! #ifdef W3_MPI INCLUDE "mpif.h" @@ -501,7 +504,7 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & !/ ------------------------------------------------------------------- / !/ Parameter list !/ - INTEGER, INTENT(IN) :: IMOD, TEND(2),ODAT(35) + INTEGER, INTENT(IN) :: IMOD, TEND(2),ODAT(40) LOGICAL, INTENT(IN), OPTIONAL :: STAMP, NO_OUT #ifdef W3_OASIS INTEGER, INTENT(IN), OPTIONAL :: ID_LCOMM @@ -511,6 +514,9 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & !/ ------------------------------------------------------------------- / !/ Local parameters : !/ +#ifdef W3_T + INTEGER :: ILEN +#endif #ifdef W3_S INTEGER, SAVE :: IENT = 0 #endif @@ -523,12 +529,15 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & INTEGER :: TTEST(2),DTTEST REAL :: ICEDAVE ! +#ifdef W3_MPI LOGICAL :: SBSED - LOGICAL :: CPLWRTFLG +#endif #ifdef W3_SEC1 INTEGER :: ISEC1 #endif +#ifdef W3_SBS INTEGER :: JJ, NDSOFLG +#endif #ifdef W3_MPI INTEGER :: IERR_MPI, NRQMAX INTEGER, ALLOCATABLE :: STATCO(:,:), STATIO(:,:) @@ -561,8 +570,8 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & ! LOGICAL :: FLACT, FLZERO, FLFRST, FLMAP, TSTAMP,& SKIP_O, FLAG_O, FLDDIR, READBC, & - FLAG0 = .FALSE., FLOUTG, FLPFLD, & - FLPART, LOCAL, FLOUTG2 + FLAG0 = .FALSE., FLOUTG = .false., FLPFLD, & + FLPART, LOCAL, FLOUTG2 = .false. ! #ifdef W3_MPI LOGICAL :: FLGMPI(0:8) @@ -583,7 +592,9 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & REAL :: VD_SPEC(NSPEC) #endif ! +#ifdef W3_SBS CHARACTER(LEN=30) :: FOUTNAME +#endif ! #ifdef W3_T REAL :: INDSORT(NSEA), DTCFL1(NSEA) @@ -594,28 +605,9 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & REAL, ALLOCATABLE :: BACSPEC(:) REAL :: BACANGL #endif - ! locally defined flags -#ifdef W3_SBS - logical, parameter :: w3_sbs_flag = .true. -#else - logical, parameter :: w3_sbs_flag = .false. -#endif -#ifdef W3_CESMCOUPLED - logical, parameter :: w3_cesmcoupled_flag = .true. -#else - logical, parameter :: w3_cesmcoupled_flag = .false. -#endif - integer :: memunit - logical :: do_gridded_output - logical :: do_point_output - logical :: do_track_output - logical :: do_restart_output - logical :: do_sf_output - logical :: do_coupler_output - logical :: do_wavefield_separation_output - logical :: do_startall - logical :: do_w3outg - + integer :: memunit + character(len=16) :: user_timestring !YYYY-MM-DD-SSSSS + character(len=256) :: fname !/ ------------------------------------------------------------------- / ! 0. Initializations ! @@ -707,12 +699,15 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & FACX = 1. END IF ! +#ifdef W3_SBS + NDSOFLG = 99 +#endif +#ifdef W3_MPI SBSED = .FALSE. - if (w3_sbs_flag) then - NDSOFLG = 99 - SBSED = .TRUE. - end if - +#endif +#ifdef W3_SBS + SBSED = .TRUE. +#endif ! TAUWX = 0. TAUWY = 0. @@ -720,7 +715,8 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & ! 0.d Test output ! #ifdef W3_T - WRITE (NDST,9000) IMOD, trim(FILEXT), TEND + ILEN = LEN_TRIM(FILEXT) + WRITE (NDST,9000) IMOD, FILEXT(:ILEN), TEND #endif ! ! 1. Check the consistency of the input ----------------------------- / @@ -2305,10 +2301,6 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & END DO IF (IT.GT.0) DTG=DTGTEMP #endif - - - - ! ! ! 3.8 Update global time step. @@ -2321,7 +2313,7 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & DTG = DTTST / REAL(NT-IT) END IF ! - IF ( FLACT .AND. IT.NE.NT .AND. IAPROC.EQ.NAPLOG ) THEN + IF ( FLACT .AND. IT.NE.NT .AND. IAPROC.EQ.NAPLOG .and. verboselog) THEN CALL STME21 ( TIME , IDTIME ) IF ( IDLAST .NE. TIME(1) ) THEN WRITE (NDSO,900) ITIME, IPASS, IDTIME(01:19), IDACT, OUTID @@ -2341,7 +2333,7 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & #endif ! ! - END DO + END DO ! DO IT = IT0, NT #ifdef W3_TIMINGS CALL PRINT_MY_TIME("W3WAVE, step 6.21.1") @@ -2362,6 +2354,26 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & ! Delay if data assimilation time. ! ! + if (dsec21(time,tend) == 0.0) then ! req'd in case waves are running in slow loop + if (use_historync) then + floutg = .false. + floutg2 = .false. + if (histwr) then + call w3cprt (imod) + call w3outg (va, flpfld, .true., .false. ) + call write_history(tend) + end if + end if + + if (use_restartnc) then + if (rstwr) then + call set_user_timestring(tend,user_timestring) + fname = trim(user_restfname)//trim(user_timestring)//'.nc' + call write_restart(trim(fname), va, mapsta+8*mapst2) + end if + end if + end if + IF ( TOFRST(1) .EQ. -1 ) THEN DTTST = 1. ELSE @@ -2389,88 +2401,82 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & ! ! 4.b Processing and MPP preparations ! - IF ( FLOUT(1) ) THEN - FLOUTG = DSEC21(TIME,TONEXT(:,1)).EQ.0. - ELSE - FLOUTG = .FALSE. - END IF - ! - IF ( FLOUT(7) ) THEN - FLOUTG2 = DSEC21(TIME,TONEXT(:,7)).EQ.0. - ELSE - FLOUTG2 = .FALSE. - END IF - ! - FLPART = .FALSE. - IF ( FLOUT(1) .AND. FLPFLD ) FLPART = FLPART .OR. DSEC21(TIME,TONEXT(:,1)).EQ.0. - IF ( FLOUT(6) ) FLPART = FLPART .OR. DSEC21(TIME,TONEXT(:,6)).EQ.0. - ! + if (.not. use_historync) then + IF ( FLOUT(1) ) THEN + FLOUTG = DSEC21(TIME,TONEXT(:,1)).EQ.0. + ELSE + FLOUTG = .FALSE. + END IF + ! + IF ( FLOUT(7) ) THEN + FLOUTG2 = DSEC21(TIME,TONEXT(:,7)).EQ.0. + ELSE + FLOUTG2 = .FALSE. + END IF + ! + FLPART = .FALSE. + IF ( FLOUT(1) .AND. FLPFLD ) FLPART = FLPART .OR. DSEC21(TIME,TONEXT(:,1)).EQ.0. + IF ( FLOUT(6) ) FLPART = FLPART .OR. DSEC21(TIME,TONEXT(:,6)).EQ.0. + ! #ifdef W3_T - WRITE (NDST,9042) LOCAL, FLPART, FLOUTG + WRITE (NDST,9042) LOCAL, FLPART, FLOUTG #endif + ! + IF ( LOCAL .AND. FLPART ) CALL W3CPRT ( IMOD ) + IF ( LOCAL .AND. (FLOUTG .OR. FLOUTG2) ) then + CALL W3OUTG ( VA, FLPFLD, FLOUTG, FLOUTG2 ) + end if + end if ! if (.not. use_historync) then ! - IF ( LOCAL .AND. FLPART ) then - CALL W3CPRT ( IMOD ) - end IF - - do_w3outg = .false. - if (w3_cesmcoupled_flag .and. histwr) then - do_w3outg = .true. - else if ( LOCAL .AND. (FLOUTG .OR. FLOUTG2) ) then - do_w3outg = .true. - end if - if (do_w3outg) then - CALL W3OUTG ( VA, FLPFLD, FLOUTG, FLOUTG2 ) - end if ! #ifdef W3_MPI FLGMPI = .FALSE. NRQMAX = 0 +#endif ! - do_startall = .false. - if (w3_cesmcoupled_flag .and. histwr) then - IF ( FLOUT(1) .OR. FLOUT(7) ) THEN - do_startall = .true. - end IF - else - CPLWRTFLG=.FALSE. - IF ( FLOUT(7) .AND. SBSED ) THEN - IF (DSEC21(TIME,TONEXT(:,7)).EQ.0.) THEN - CPLWRTFLG=.TRUE. - END IF - END IF - IF ( ( (DSEC21(TIME,TONEXT(:,1)).EQ.0.) .AND. FLOUT(1) ) .OR. & - ( CPLWRTFLG ) ) THEN - do_startall = .true. - end IF - end if - if (do_startall) then +#ifdef W3_MPI + IF ( (FLOUTG) .OR. (FLOUTG2 .AND. SBSED) ) THEN IF (.NOT. LPDLIB) THEN IF (NRQGO.NE.0 ) THEN +#endif +#ifdef W3_MPI CALL MPI_STARTALL ( NRQGO, IRQGO , IERR_MPI ) +#endif +#ifdef W3_MPI FLGMPI(0) = .TRUE. NRQMAX = MAX ( NRQMAX , NRQGO ) +#endif #ifdef W3_MPIT WRITE (NDST,9043) '1a', NRQGO, NRQMAX, NAPFLD #endif +#ifdef W3_MPI END IF +#endif ! +#ifdef W3_MPI IF (NRQGO2.NE.0 ) THEN +#endif +#ifdef W3_MPI CALL MPI_STARTALL ( NRQGO2, IRQGO2, IERR_MPI ) - +#endif +#ifdef W3_MPI FLGMPI(1) = .TRUE. NRQMAX = MAX ( NRQMAX , NRQGO2 ) +#endif #ifdef W3_MPIT WRITE (NDST,9043) '1b', NRQGO2, NRQMAX, NAPFLD #endif +#ifdef W3_MPI END IF ELSE +#endif #ifdef W3_PDLIB CALL DO_OUTPUT_EXCHANGES(IMOD) #endif +#ifdef W3_MPI END IF ! IF (.NOT. LPDLIB) THEN - END IF ! if (do_startall) + END IF #endif call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE AFTER TIME LOOP 1') ! @@ -2490,34 +2496,36 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & #endif ! #ifdef W3_MPI - IF ( FLOUT(4) .AND. NRQRS.NE.0 ) THEN - IF ( DSEC21(TIME,TONEXT(:,4)).EQ.0. ) THEN - CALL MPI_STARTALL ( NRQRS, IRQRS , IERR_MPI ) - FLGMPI(4) = .TRUE. - NRQMAX = MAX ( NRQMAX , NRQRS ) + if (.not. use_restartnc) then + IF ( FLOUT(4) .AND. NRQRS.NE.0 ) THEN + IF ( DSEC21(TIME,TONEXT(:,4)).EQ.0. ) THEN + CALL MPI_STARTALL ( NRQRS, IRQRS , IERR_MPI ) + FLGMPI(4) = .TRUE. + NRQMAX = MAX ( NRQMAX , NRQRS ) #endif #ifdef W3_MPIT - WRITE (NDST,9043) '4 ', NRQRS, NRQMAX, NAPRST + WRITE (NDST,9043) '4 ', NRQRS, NRQMAX, NAPRST #endif #ifdef W3_MPI + END IF END IF - END IF #endif - ! + ! #ifdef W3_MPI - IF ( FLOUT(8) .AND. NRQRS.NE.0 ) THEN - IF ( DSEC21(TIME,TONEXT(:,8)).EQ.0. ) THEN - CALL MPI_STARTALL ( NRQRS, IRQRS , IERR_MPI ) - FLGMPI(8) = .TRUE. - NRQMAX = MAX ( NRQMAX , NRQRS ) + IF ( FLOUT(8) .AND. NRQRS.NE.0 ) THEN + IF ( DSEC21(TIME,TONEXT(:,8)).EQ.0. ) THEN + CALL MPI_STARTALL ( NRQRS, IRQRS , IERR_MPI ) + FLGMPI(8) = .TRUE. + NRQMAX = MAX ( NRQMAX , NRQRS ) #endif #ifdef W3_MPIT - WRITE (NDST,9043) '8 ', NRQRS, NRQMAX, NAPRST + WRITE (NDST,9043) '8 ', NRQRS, NRQMAX, NAPRST #endif #ifdef W3_MPI + END IF END IF - END IF #endif + end if ! if (.not. use_restartnc) ! #ifdef W3_MPI IF ( FLOUT(5) .AND. NRQBP.NE.0 ) THEN @@ -2554,7 +2562,6 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE AFTER TIME LOOP 2') ! ! 4.c Reset next output time - ! TOFRST(1) = -1 TOFRST(2) = 0 @@ -2562,29 +2569,6 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & DO J=1, NOTYPE IF ( FLOUT(J) ) THEN - ! - ! - ! Determine output flags - ! - if (w3_sbs_flag) then - do_gridded_output = ( j .eq. 1 ) .or. ( j .eq. 7 ) - else - if (w3_cesmcoupled_flag) then - do_gridded_output = ( j .eq. 1 ) .and. histwr - else - do_gridded_output = ( j .eq. 1 ) - end if - end if - do_point_output = (j .eq. 2) - do_track_output = (j .eq. 3) - if (w3_cesmcoupled_flag) then - do_restart_output = (j .eq. 4) .and. rstwr - else - do_restart_output = (j .eq. 4) - end if - do_wavefield_separation_output = (j .eq. 5) - do_sf_output = (j .eq. 6) - do_coupler_output = (j .eq. 7) ! ! 4.d Perform output ! @@ -2595,75 +2579,82 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & DTTST = DSEC21 ( TIME, TOUT ) ! IF ( DTTST .EQ. 0. ) THEN - if (do_gridded_output) then - if (user_netcdf_grdout) then -#ifdef W3_MPI - IF ( FLGMPI(0) )CALL MPI_WAITALL( NRQGO, IRQGO, STATIO, IERR_MPI ) - FLGMPI(0) = .FALSE. + IF ( ( J .EQ. 1 ) & +#ifdef W3_SBS + .OR. ( J .EQ. 7 ) & #endif - IF ( IAPROC .EQ. NAPFLD ) THEN + .and. .not. use_historync) THEN + IF ( IAPROC .EQ. NAPFLD ) THEN #ifdef W3_MPI - IF ( FLGMPI(1) ) CALL MPI_WAITALL( NRQGO2, IRQGO2, STATIO, IERR_MPI ) - FLGMPI(1) = .FALSE. + IF ( FLGMPI(1) ) CALL MPI_WAITALL ( NRQGO2, IRQGO2, STATIO, IERR_MPI ) + FLGMPI(1) = .FALSE. #endif - CALL W3IOGONCD () - END IF - else - ! default (binary) output - IF ( IAPROC .EQ. NAPFLD ) THEN -#ifdef W3_MPI - IF ( FLGMPI(1) ) CALL MPI_WAITALL( NRQGO2, IRQGO2, STATIO, IERR_MPI ) - FLGMPI(1) = .FALSE. -#endif - if (w3_sbs_flag) then - IF ( J .EQ. 1 ) THEN - CALL W3IOGO( 'WRITE', NDS(7), ITEST, IMOD ) - ENDIF - - ! Generate output flag file for fields and SBS coupling. - CALL STME21 ( TIME, IDTIME ) - FOUTNAME = 'Field_done.' // IDTIME(1:4) & - // IDTIME(6:7) // IDTIME(9:10) & - // IDTIME(12:13) // '.' // TRIM(FILEXT) - OPEN( UNIT=NDSOFLG, FILE=FOUTNAME) - CLOSE( NDSOFLG ) - else - CALL W3IOGO( 'WRITE', NDS(7), ITEST, IMOD ) - endif - end if - end if ! user_netcdf_grdout - - ELSE IF ( do_point_output ) THEN + ! +#ifdef W3_SBS + IF ( J .EQ. 1 ) THEN +#endif + CALL W3IOGO( 'WRITE', NDS(7), ITEST, IMOD ) +#ifdef W3_SBS + ENDIF +#endif + ! +#ifdef W3_SBS + ! + ! Generate output flag file for fields and SBS coupling. + ! + JJ = LEN_TRIM ( FILEXT ) + CALL STME21 ( TIME, IDTIME ) + FOUTNAME = 'Field_done.' // IDTIME(1:4) & + // IDTIME(6:7) // IDTIME(9:10) & + // IDTIME(12:13) // '.' // FILEXT(1:JJ) +#endif + ! +#ifdef W3_SBS + OPEN( UNIT=NDSOFLG, FILE=FOUTNAME) + CLOSE( NDSOFLG ) +#endif + END IF + ! + ELSE IF ( J .EQ. 2 ) THEN + ! + ! Point output + ! IF ( IAPROC .EQ. NAPPNT ) THEN + ! + ! Gets the necessary spectral data + ! CALL W3IOPE ( VA ) CALL W3IOPO ( 'WRITE', NDS(8), ITEST, IMOD ) END IF - - ELSE IF ( do_track_output ) THEN + ! + ELSE IF ( J .EQ. 3 ) THEN + ! + ! Track output + ! CALL W3IOTR ( NDS(11), NDS(12), VA, IMOD ) - - ELSE IF ( do_restart_output ) THEN + ELSE IF ( J .EQ. 4 .and. .not. use_restartnc) THEN CALL W3IORS ('HOT', NDS(6), XXX, IMOD, FLOUT(8) ) ITEST = RSTYPE - - ELSE IF ( do_wavefield_separation_output ) THEN + ELSE IF ( J .EQ. 5 ) THEN IF ( IAPROC .EQ. NAPBPT ) THEN #ifdef W3_MPI IF (NRQBP2.NE.0) CALL MPI_WAITALL ( NRQBP2, IRQBP2,STATIO, IERR_MPI ) #endif - CALL W3IOBC ( 'WRITE', NDS(10), TIME, TIME, ITEST, IMOD ) + CALL W3IOBC ( 'WRITE', NDS(10), & + TIME, TIME, ITEST, IMOD ) END IF - ELSE IF ( do_sf_output ) THEN + ELSE IF ( J .EQ. 6 ) THEN CALL W3IOSF ( NDS(13), IMOD ) #ifdef W3_OASIS - ELSE IF ( do_coupler_output ) THEN + ELSE IF ( J .EQ. 7 ) THEN ! ! Send variables to atmospheric or ocean circulation or ice model ! IF (DTOUT(7).NE.0) THEN IF ( (MOD(ID_OASIS_TIME,NINT(DTOUT(7))) .EQ. 0 ) .AND. & (DSEC21 (TIME00, TIME) .GT. 0.0) ) THEN - IF ( (CPLT0 .AND. (DSEC21 (TIME, TIMEN) .GT. 0.0)) .OR. .NOT. CPLT0 ) THEN + IF ( (CPLT0 .AND. (DSEC21 (TIME, TIMEN) .GT. 0.0)) .OR. & + .NOT. CPLT0 ) THEN IF (CPLT0) ID_OASIS_TIME = NINT(DSEC21 ( TIME00 , TIME )) #endif @@ -2722,7 +2713,7 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & ! If there is a second stream of restart files then J=8 and FLOUT(8)=.TRUE. J=8 - IF ( FLOUT(J) ) THEN + IF ( FLOUT(J) .and. .not. use_restartnc) THEN ! ! 4.d Perform output ! @@ -2767,11 +2758,6 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & ! #ifdef W3_MPI IF ( FLGMPI(0) ) CALL MPI_WAITALL ( NRQGO, IRQGO , STATIO, IERR_MPI ) - if (user_netcdf_grdout) then - IF ( FLGMPI(1) .and. ( IAPROC .EQ. NAPFLD ) ) then - CALL MPI_WAITALL ( NRQGO2, IRQGO2 , STATIO, IERR_MPI ) - end if - end if IF ( FLGMPI(2) ) CALL MPI_WAITALL ( NRQPO, IRQPO1, STATIO, IERR_MPI ) IF ( FLGMPI(4) ) CALL MPI_WAITALL ( NRQRS, IRQRS , STATIO, IERR_MPI ) IF ( FLGMPI(8) ) CALL MPI_WAITALL ( NRQRS, IRQRS , STATIO, IERR_MPI ) @@ -2790,7 +2776,7 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & ! ! 5. Update log file ------------------------------------------------ / ! - IF ( IAPROC.EQ.NAPLOG ) THEN + IF ( IAPROC.EQ.NAPLOG .and. verboselog) THEN ! CALL STME21 ( TIME , IDTIME ) IF ( FLCUR ) THEN @@ -2843,7 +2829,7 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & WRITE (SCREEN,951) STTIME END IF - IF ( IAPROC .EQ. NAPLOG ) WRITE (NDSO,902) + IF ( IAPROC .EQ. NAPLOG .and. verboselog) WRITE (NDSO,902) ! DEALLOCATE(FIELD) DEALLOCATE(TAUWX, TAUWY) diff --git a/model/src/wav_comp_nuopc.F90 b/model/src/wav_comp_nuopc.F90 index c7bb14ef2..4280b3b14 100644 --- a/model/src/wav_comp_nuopc.F90 +++ b/model/src/wav_comp_nuopc.F90 @@ -44,8 +44,8 @@ module wav_comp_nuopc use wav_shr_mod , only : wav_coupling_to_cice, nwav_elev_spectrum use wav_shr_mod , only : merge_import, dbug_flag use w3odatmd , only : nds, iaproc, napout - use w3odatmd , only : runtype, use_user_histname, user_histfname, use_user_restname, user_restfname - use w3odatmd , only : user_netcdf_grdout + use w3odatmd , only : runtype, user_histfname, user_restfname, verboselog + use w3odatmd , only : use_historync, use_restartnc, restart_from_binary, logfile_is_assigned use w3odatmd , only : time_origin, calendar_name, elapsed_secs use wav_shr_mod , only : casename, multigrid, inst_suffix, inst_index, unstr_mesh use wav_wrapper_mod , only : ufs_settimer, ufs_logtimer, ufs_file_setlogunit, wtime @@ -92,14 +92,6 @@ module wav_comp_nuopc #endif integer, allocatable :: tend(:,:) !< the ending time of ModelAdvance when !! run with multigrid=true - logical :: user_histalarm = .false. !< logical flag for user to set history alarms - !! using ESMF. If history_option is present as config - !! option, user_histalarm will be true and will be - !! set using history_option, history_n and history_ymd - logical :: user_restalarm = .false. !< logical flag for user to set restart alarms - !! using ESMF. If restart_option is present as config - !! option, user_restalarm will be true and will be - !! set using restart_option, restart_n and restart_ymd integer :: ymd !< current year-month-day integer :: tod !< current time of day (sec) integer :: time0(2) !< start time stored as yyyymmdd,hhmmss @@ -229,6 +221,7 @@ end subroutine InitializeP0 subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) use wav_shr_flags, only : w3_pdlib_flag + ! input/output arguments type(ESMF_GridComp) :: gcomp type(ESMF_State) :: importState, exportState @@ -383,6 +376,14 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) if (runtimelog) then call ufs_file_setLogUnit('./log.ww3.timer',nu_timer,runtimelog) end if + + ! Determine verbose native WW3 logging + call NUOPC_CompAttributeGet(gcomp, name="verboselog", value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) verboselog=(trim(cvalue)=="true") + write(logmsg,*) verboselog + call ESMF_LogWrite('WW3_cap: Verbose WW3 native logging is = '//trim(logmsg), ESMF_LOGMSG_INFO) + call advertise_fields(importState, exportState, flds_scalar_name, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -415,23 +416,24 @@ end subroutine InitializeAdvertise !> @date 01-05-2022 subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) - use w3odatmd , only : w3nout, w3seto, naproc, iaproc, naperr, napout - use w3timemd , only : stme21 - use w3adatmd , only : w3naux, w3seta - use w3idatmd , only : w3seti, w3ninp - use w3gdatmd , only : nk, nseal, nsea, nx, ny, mapsf, w3nmod, w3setg - use w3gdatmd , only : rlgtype, ungtype, gtype - use w3wdatmd , only : va, time, w3ndat, w3dimw, w3setw - use w3parall , only : init_get_isea + use w3odatmd , only : w3nout, w3seto, naproc, naperr + use w3timemd , only : stme21 + use w3adatmd , only : w3naux, w3seta + use w3idatmd , only : w3seti, w3ninp + use w3gdatmd , only : nk, nseal, nsea, nx, ny, mapsf, w3nmod, w3setg + use w3gdatmd , only : rlgtype, ungtype, gtype + use w3wdatmd , only : va, time, w3ndat, w3dimw, w3setw + use w3parall , only : init_get_isea #ifndef W3_CESMCOUPLED - use wminitmd , only : wminit, wminitnml - use wmunitmd , only : wmuget, wmuset + use wminitmd , only : wminit, wminitnml + use wmunitmd , only : wmuget, wmuset #endif - use wav_shel_inp , only : set_shel_io - use wav_grdout , only : wavinit_grdout - use wav_shr_mod , only : diagnose_mesh, write_meshdecomp + use wav_shel_inp , only : set_shel_io + use wav_history_mod , only : wav_history_init + use wav_pio_mod , only : wav_pio_init + use wav_shr_mod , only : diagnose_mesh, write_meshdecomp, wav_loginit #ifdef W3_PDLIB - use yowNodepool , only : ng + use yowNodepool , only : ng #endif ! input/output variables @@ -450,6 +452,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) type(ESMF_TimeInterval) :: TimeOffset type(ESMF_TimeInterval) :: TimeStep type(ESMF_Calendar) :: calendar + type(ESMF_Info) :: info character(CL) :: cvalue integer :: shrlogunit integer :: yy,mm,dd,hh,ss @@ -471,7 +474,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) integer(i4) :: maskmin integer(i4), pointer :: meshmask(:) character(23) :: dtme21 - integer :: iam, mpi_comm + integer :: iam, mpi_comm, num_threads character(ESMF_MAXSTR) :: msgString character(ESMF_MAXSTR) :: diro character(CL) :: logfile @@ -481,6 +484,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) integer :: stdout integer :: petcount real(r8) :: toff + logical :: isPresent, isSet character(ESMF_MAXSTR) :: preamb = './' character(ESMF_MAXSTR) :: ifname = 'ww3_multi.inp' character(len=*), parameter :: subname = '(wav_comp_nuopc:InitializeRealize)' @@ -517,6 +521,12 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) call ESMF_VMGet(vm, mpiCommunicator=mpi_comm, peCount=petcount, localPet=iam, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_InfoGetFromHost(gcomp, info=info, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_InfoGet(info, key="/NUOPC/Hint/PePerPet/MaxCount", value=num_threads, default=1, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + #ifndef W3_CESMCOUPLED nmproc = petcount #else @@ -555,11 +565,17 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) call NUOPC_CompAttributeGet(gcomp, name="logfile", value=logfile, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return open (newunit=stdout, file=trim(diro)//"/"//trim(logfile)) + logfile_is_assigned = .true. else stdout = 6 endif else - stdout = 6 + if ( root_task ) then + open (newunit=stdout, file='log.ww3') + logfile_is_assigned = .true. + else + stdout = 6 + end if end if if (.not. multigrid) call set_shel_io(stdout,mds,ntrace) @@ -567,6 +583,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) if ( root_task ) then write(stdout,'(a)')' *** WAVEWATCH III Program shell *** ' write(stdout,'(a)')'===============================================' + write(stdout,'(/)') write(stdout,'(a,l)')' Wave wav_coupling_to_cice setting is ',wav_coupling_to_cice end if @@ -584,7 +601,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) runtype = "branch" end if if ( root_task ) then - write(stdout,*) 'WW3 runtype is '//trim(runtype) + write(stdout,'(a)') ' WW3 runtype is '//trim(runtype) end if call ESMF_LogWrite('WW3 runtype is '//trim(runtype), ESMF_LOGMSG_INFO) @@ -623,7 +640,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return endif ! Determine time attributes for history output - call ESMF_TimeGet( esmfTime, timeString=time_origin, calendar=calendar, rc=rc ) + call ESMF_TimeGet( startTime, timeString=time_origin, calendar=calendar, rc=rc ) if (ChkErr(rc,__LINE__,u_FILE_u)) return time_origin = 'seconds since '//time_origin(1:10)//' '//time_origin(12:19) !call ESMF_ClockGet(clock, calendar=calendar) @@ -659,13 +676,82 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) call stme21 ( time0 , dtme21 ) if ( root_task ) then write (stdout,'(a)')' Starting time : '//trim(dtme21) - write (stdout,'(a,i8,2x,i8)') 'start_ymd, stop_ymd = ',start_ymd, stop_ymd + write (stdout,'(a,i8,2x,i8)') ' start_ymd, stop_ymd = ',start_ymd, stop_ymd end if #ifndef W3_CESMCOUPLED stime = time0 etime = timen #endif + !-------------------------------------------------------------------- + ! Initialize PIO. This needs to be done prior to the call to w3init + ! in order to read the restart file. The filename strings must also + ! be available + !-------------------------------------------------------------------- + + if (cesmcoupled) then + if (len_trim(inst_suffix) > 0) then + user_restfname = trim(casename)//'.ww3'//trim(inst_suffix)//'.r.' + user_histfname = trim(casename)//'.ww3'//trim(inst_suffix)//'.hi.' + else + user_restfname = trim(casename)//'.ww3.r.' + user_histfname = trim(casename)//'.ww3.hi.' + endif + + ! netcdf is used for CESM history and restart + use_historync = .true. + use_restartnc = .true. + else + call NUOPC_CompAttributeGet(gcomp, name='use_restartnc', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + use_restartnc=(trim(cvalue)=="true") + end if + if (root_task) write(stdout,'(a,l4)') trim(subname)//': Wave use_restartnc setting is ',use_restartnc + + ! user filenaming is required with netcdf restarts or restart_from_binary. If netcdf restarts are not used, + ! only native WW3 file naming is possible + if (use_restartnc) then + user_restfname = trim(casename)//'.ww3.r.' + if (root_task) write(stdout,'(a)') trim(subname)//': Custom restart prefix is '//trim(user_restfname) + end if + + call NUOPC_CompAttributeGet(gcomp, name='use_historync', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + use_historync=(trim(cvalue)=="true") + end if + if (root_task) write(stdout,'(a,l4)') trim(subname)//': Wave use_historync setting is ',use_historync + + ! user filenaming is optional with netcdf output. If netcdf history is not used, only native WW3 + ! naming is possible + if (use_historync) then + call NUOPC_CompAttributeGet(gcomp, name='user_histname', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (trim(cvalue)=="true") then + user_histfname = trim(casename)//'.ww3.hi.' + if (root_task) write(stdout,'(a)') trim(subname)//': Custom history prefix is '//trim(user_histfname) + else + user_histfname = '' + end if + end if + end if ! if (cesmcoupled) + + ! allow startup from binary restarts as special case + if (use_restartnc) then + call NUOPC_CompAttributeGet(gcomp, name='restart_from_binary', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + restart_from_binary=(trim(cvalue)=="true") + end if + if (root_task) write(stdout,'(a,l4)') trim(subname)//': Wave restart_from_binary setting is ',restart_from_binary + end if + + if (use_restartnc .or. use_historync) then + call wav_pio_init(gcomp, mpi_comm, stdout, naproc/num_threads, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + !-------------------------------------------------------------------- ! Wave model initialization !-------------------------------------------------------------------- @@ -694,7 +780,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) end do call ESMF_LogWrite(trim(subname)//' done = wminit', ESMF_LOGMSG_INFO) else - call waveinit_ufs(gcomp, ntrace, mpi_comm, mds, rc) + call waveinit_ufs(gcomp, stdout, ntrace, mpi_comm, mds, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if #else @@ -704,11 +790,10 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) call waveinit_cesm(gcomp, ntrace, mpi_comm, mds, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return #endif - - ! call mpi_barrier ( mpi_comm, ierr ) + !call mpi_barrier ( mpi_comm, ierr ) if ( root_task ) then - inquire(unit=nds(1), name=logfile) - print *,'WW3 log written to '//trim(logfile) + inquire(unit=stdout, name=logfile) + write(*,'(a)')'WW3 log written to '//trim(logfile) end if if (wav_coupling_to_cice) then @@ -718,14 +803,6 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) end if end if - !-------------------------------------------------------------------- - ! Intialize the list of requested output variables for netCDF output - !-------------------------------------------------------------------- - - if (user_netcdf_grdout) then - call wavinit_grdout - end if - !-------------------------------------------------------------------- ! Mesh initialization !-------------------------------------------------------------------- @@ -761,7 +838,6 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) ! create distGrid from global index array of sea points with no ghost points DistGrid = ESMF_DistGridCreate(arbSeqIndexList=gindex_sea, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - deallocate(gindex_sea) else ! create a global index array for non-sea (i.e. land points) allocate(mask_global(nx*ny), mask_local(nx*ny)) @@ -807,8 +883,6 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) gindex(ncnt) = gindex_lnd(ncnt-nseal_cpl) end if end do - deallocate(gindex_sea) - deallocate(gindex_lnd) ! create distGrid from global index array DistGrid = ESMF_DistGridCreate(arbSeqIndexList=gindex, rc=rc) @@ -818,14 +892,23 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) ! get the mesh file name call NUOPC_CompAttributeGet(gcomp, name='mesh_wav', value=cvalue, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - - ! read in the mesh with the above DistGrid + ! read in the mesh with the the DistGrid EMesh = ESMF_MeshCreate(filename=trim(cvalue), fileformat=ESMF_FILEFORMAT_ESMFMESH, & elementDistgrid=Distgrid,rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (dbug_flag > 5) then - call diagnose_mesh(EMesh, size(gindex), 'EMesh', rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (unstr_mesh) then + call diagnose_mesh(EMesh, size(gindex_sea), 'EMesh', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + deallocate(gindex_sea) + else + call diagnose_mesh(EMesh, size(gindex), 'EMesh', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + deallocate(gindex) + deallocate(gindex_sea) + deallocate(gindex_lnd) + end if end if if (.not. unstr_mesh) then @@ -857,7 +940,6 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if deallocate(meshmask) - deallocate(gindex) end if if (dbug_flag > 5) then @@ -887,6 +969,24 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) enddo end if #endif + !-------------------------------------------------------------------- + ! Intialize the list of requested output variables for netCDF output. + ! This needs to occur after mod_def has been read in w3init since + ! some variables are available only if they are defined in the mod_def + !-------------------------------------------------------------------- + + if (use_historync) then + call wav_history_init(stdout) + end if + + !-------------------------------------------------------------------- + ! Write the header string for WW3 native logging + !-------------------------------------------------------------------- + + if (root_task) then + if (verboselog) call wav_loginit(stdout) + end if + if (root_task) call ufs_logtimer(nu_timer,time,start_tod,'InitializeRealize time: ',runtimelog,wtime) if (dbug_flag > 5) call ESMF_LogWrite(trim(subname)//' done', ESMF_LOGMSG_INFO) @@ -1026,7 +1126,7 @@ subroutine ModelAdvance(gcomp, rc) type(ESMF_Time) :: currTime, nextTime, startTime, stopTime integer :: yy,mm,dd,hh,ss integer :: imod - integer :: shrlogunit ! original log unit and level + !integer :: shrlogunit ! original log unit and level character(ESMF_MAXSTR) :: msgString character(len=*),parameter :: subname = '(wav_comp_nuopc:ModelAdvance) ' !------------------------------------------------------- @@ -1062,8 +1162,8 @@ subroutine ModelAdvance(gcomp, rc) ss = tod - (hh*3600) - (mm*60) time0(1) = ymd time0(2) = hh*10000 + mm*100 + ss - if ( root_task ) then - write(nds(1),'(a,3i4,i10)') 'ymd2date currTime wav_comp_nuopc hh,mm,ss,ymd', hh,mm,ss,ymd + if (dbug_flag > 5) then + if ( root_task ) write(nds(1),'(a,3i4,i10)') 'ymd2date currTime wav_comp_nuopc hh,mm,ss,ymd', hh,mm,ss,ymd end if if (root_task) call ufs_logtimer(nu_timer,time,tod,'ModelAdvance time since last step: ',runtimelog,wtime) call ufs_settimer(wtime) @@ -1108,41 +1208,29 @@ subroutine ModelAdvance(gcomp, rc) !------------ if(profile_memory) call ESMF_VMLogMemInfo("Entering WW3 Run : ") - if (user_restalarm) then - ! Determine if time to write ww3 restart files - call ESMF_ClockGetAlarm(clock, alarmname='alarm_restart', alarm=alarm, rc=rc) + ! Determine if time to write ww3 restart files + call ESMF_ClockGetAlarm(clock, alarmname='alarm_restart', alarm=alarm, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (ESMF_AlarmIsRinging(alarm, rc=rc)) then + if (ChkErr(rc,__LINE__,u_FILE_u)) return + rstwr = .true. + call ESMF_AlarmRingerOff( alarm, rc=rc ) if (ChkErr(rc,__LINE__,u_FILE_u)) return - if (ESMF_AlarmIsRinging(alarm, rc=rc)) then - if (ChkErr(rc,__LINE__,u_FILE_u)) return - rstwr = .true. - call ESMF_AlarmRingerOff( alarm, rc=rc ) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - else - rstwr = .false. - endif else rstwr = .false. - end if + endif - if (user_histalarm) then - ! Determine if time to write ww3 history files - call ESMF_ClockGetAlarm(clock, alarmname='alarm_history', alarm=alarm, rc=rc) + ! Determine if time to write ww3 history files + call ESMF_ClockGetAlarm(clock, alarmname='alarm_history', alarm=alarm, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (ESMF_AlarmIsRinging(alarm, rc=rc)) then + if (ChkErr(rc,__LINE__,u_FILE_u)) return + histwr = .true. + call ESMF_AlarmRingerOff( alarm, rc=rc ) if (ChkErr(rc,__LINE__,u_FILE_u)) return - if (ESMF_AlarmIsRinging(alarm, rc=rc)) then - if (ChkErr(rc,__LINE__,u_FILE_u)) return - histwr = .true. - call ESMF_AlarmRingerOff( alarm, rc=rc ) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - else - histwr = .false. - endif else histwr = .false. - end if - if ( root_task ) then - ! write(nds(1),*) 'wav_comp_nuopc time', time, timen - ! write(nds(1),*) 'ww3 hist flag ', histwr, hh - end if + endif ! Advance the wave model #ifndef W3_CESMCOUPLED @@ -1179,6 +1267,7 @@ end subroutine ModelAdvance !> @date 01-05-2022 subroutine ModelSetRunClock(gcomp, rc) + use wav_shel_inp , only : odat ! input/output variables type(ESMF_GridComp) :: gcomp integer, intent(out) :: rc @@ -1189,6 +1278,7 @@ subroutine ModelSetRunClock(gcomp, rc) type(ESMF_Time) :: mstoptime type(ESMF_Time) :: mstarttime type(ESMF_TimeInterval) :: mtimestep, dtimestep + character(ESMF_MAXSTR) :: msgString logical :: isPresent logical :: isSet character(len=256) :: cvalue @@ -1273,12 +1363,6 @@ subroutine ModelSetRunClock(gcomp, rc) call ESMF_AlarmSet(restart_alarm, clock=mclock, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - user_restalarm = .true. - else - ! If attribute is not present - write restarts at native WW3 freq - restart_option = 'none' - restart_n = -999 - user_restalarm = .false. end if !---------------- @@ -1331,14 +1415,22 @@ subroutine ModelSetRunClock(gcomp, rc) call ESMF_AlarmSet(history_alarm, clock=mclock, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - user_histalarm = .true. else - ! If attribute is not present - write history output at native WW3 frequency - history_option = 'none' - history_n = -999 - user_histalarm = .false. + ! If attribute is not present - write history output at stride frequency + history_option = 'nseconds' + history_n = odat(3) + history_ymd = -999 + call alarmInit(mclock, history_alarm, history_option, & + opt_n = history_n, & + opt_ymd = history_ymd, & + RefTime = mStartTime, & + alarmname = 'alarm_history', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_AlarmSet(history_alarm, clock=mclock, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + write(msgString,'(a,i10)')' History will be written at field%stride freq ',history_n + call ESMF_LogWrite(trim(subname)//trim(msgString), ESMF_LOGMSG_INFO) end if - end if !-------------------------------- @@ -1428,7 +1520,7 @@ subroutine waveinit_cesm(gcomp, ntrace, mpi_comm, mds, rc) ! local variables integer :: ierr integer :: unitn ! namelist unit number - integer :: shrlogunit + !integer :: shrlogunit logical :: isPresent, isSet real(r8) :: dtmax_in ! Maximum overall time step. real(r8) :: dtmin_in ! Minimum dynamic time step for source @@ -1545,23 +1637,6 @@ subroutine waveinit_cesm(gcomp, ntrace, mpi_comm, mds, rc) inflags2(-3) = .false. ! ice floe size end if - ! custom restart and history file names are used for CESM - use_user_histname = .true. - use_user_restname = .true. - - ! if runtype=initial, the initfile will be read in w3iorsmd - if (len_trim(inst_suffix) > 0) then - user_restfname = trim(casename)//'.ww3'//trim(inst_suffix)//'.r.' - user_histfname = trim(casename)//'.ww3'//trim(inst_suffix)//'.hi.' - else - user_restfname = trim(casename)//'.ww3.r.' - user_histfname = trim(casename)//'.ww3.hi.' - endif - - ! netcdf gridded output is used for CESM - user_netcdf_grdout = .true. - ! restart and history alarms are set for CESM by default through config - ! Read in initial/restart data and initialize the model ! ww3 read initialization occurs in w3iors (which is called by initmd in module w3initmd) ! ww3 always starts up from a 'restart' file type @@ -1592,6 +1667,7 @@ end subroutine waveinit_cesm !! ww3_shel.nml file. Calls w3init to initialize the wave model !! !! @param[in] gcomp an ESMF_GridComp object + !! @param[in] stdout the logfile unit on the root task !! @param[in] ntrace unit numbers for trace !! @param[in] mpi_comm an mpi communicator !! @param[in] mds unit numbers @@ -1599,78 +1675,62 @@ end subroutine waveinit_cesm !! !> @author mvertens@ucar.edu, Denise.Worthen@noaa.gov !> @date 01-05-2022 - subroutine waveinit_ufs( gcomp, ntrace, mpi_comm, mds, rc) + subroutine waveinit_ufs( gcomp, stdout, ntrace, mpi_comm, mds, rc) ! Initialize ww3 for ufs (called from InitializeRealize) - use w3odatmd , only : fnmpre + use w3odatmd , only : fnmpre, addrstflds, rstfldlist, rstfldcnt use w3gdatmd , only : dtcfl, dtcfli, dtmax, dtmin use w3initmd , only : w3init + use w3servmd , only : strsplit + use w3timemd , only : set_user_timestring use wav_shel_inp , only : read_shel_config use wav_shel_inp , only : npts, odat, iprt, x, y, pnames, prtfrm use wav_shel_inp , only : flgrd, flgd, flgr2, flg2 ! input/output variables type(ESMF_GridComp) :: gcomp + integer, intent(in) :: stdout integer, intent(in) :: ntrace(:) integer, intent(in) :: mpi_comm integer, intent(in) :: mds(:) integer, intent(out) :: rc ! local variables - character(len=CL) :: logmsg - logical :: isPresent, isSet - character(len=CL) :: cvalue - integer :: dt_in(4) + logical :: isPresent, isSet + character(len=CL) :: cvalue + character(len=CL) :: logmsg + character(len=CL) :: fldrst = '' + character(len=100) :: tmplist(100) = '' + integer :: dt_in(4) + integer :: i, cnt character(len=*), parameter :: subname = '(wav_comp_nuopc:wavinit_ufs)' ! ------------------------------------------------------------------- rc = ESMF_SUCCESS if (dbug_flag > 5) call ESMF_LogWrite(trim(subname)//' called', ESMF_LOGMSG_INFO) - ! restart and history alarms are optional for UFS and used via allcomp config settings - call NUOPC_CompAttributeGet(gcomp, name='user_sets_histname', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - if (isPresent .and. isSet) then - use_user_histname=(trim(cvalue)=="true") - end if - write(logmsg,'(A,l)') trim(subname)//': Custom history names in use ',use_user_histname - call ESMF_LogWrite(trim(logmsg), ESMF_LOGMSG_INFO) - - call NUOPC_CompAttributeGet(gcomp, name='user_sets_restname', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - if (isPresent .and. isSet) then - use_user_restname=(trim(cvalue)=="true") - end if - write(logmsg,'(A,l)') trim(subname)//': Custom restart names in use ',use_user_restname - call ESMF_LogWrite(trim(logmsg), ESMF_LOGMSG_INFO) - - call NUOPC_CompAttributeGet(gcomp, name='gridded_netcdfout', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - if (isPresent .and. isSet) then - user_netcdf_grdout=(trim(cvalue)=="true") - end if - write(logmsg,'(A,l)') trim(subname)//': Gridded netcdf output is requested ',user_netcdf_grdout - call ESMF_LogWrite(trim(logmsg), ESMF_LOGMSG_INFO) - - if (use_user_histname) then - user_histfname = trim(casename)//'.ww3.hi.' - end if - if (use_user_restname) then - user_restfname = trim(casename)//'.ww3.r.' - end if - fnmpre = './' + if (root_task) write(stdout,'(a)') trim(subname)//' call read_shel_config' + call read_shel_config(mpi_comm, mds, time0_overwrite=time0, timen_overwrite=timen, rstfldlist=fldrst) + + ! Define any additional restart fields + if(len_trim(fldrst) > 0) then + addrstflds = .true. + call strsplit(fldrst, tmplist) + do i = 1,size(rstfldlist) + rstfldlist(i) = trim(tmplist(i)) + if (len_trim(rstfldlist(i)) > 0) rstfldcnt = rstfldcnt + 1 + end do + end if - call ESMF_LogWrite(trim(subname)//' call read_shel_config', ESMF_LOGMSG_INFO) - call read_shel_config(mpi_comm, mds, time0_overwrite=time0, timen_overwrite=timen) - - call ESMF_LogWrite(trim(subname)//' call w3init', ESMF_LOGMSG_INFO) + if (root_task) write(stdout,'(a,/)') trim(subname)//' call w3init' call w3init ( 1, .false., 'ww3', mds, ntrace, odat, flgrd, flgr2, flgd, flg2, & npts, x, y, pnames, iprt, prtfrm, mpi_comm ) - write(logmsg,'(A,4f10.2)') trim(subname)//': mod_def timesteps file ',dtmax,dtcfl,dtcfli,dtmin - call ESMF_LogWrite(trim(logmsg), ESMF_LOGMSG_INFO) + write(cvalue,'(4f10.1)')dtmax,dtcfl,dtcfli,dtmin + write(logmsg,'(a)')trim(subname)//': WW3 timesteps from mod_def '//trim(cvalue) + call NUOPC_CompAttributeGet(gcomp, name='dt_in', isPresent=isPresent, isSet=isSet, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return if (isPresent .and. isSet) then @@ -1681,9 +1741,23 @@ subroutine waveinit_ufs( gcomp, ntrace, mpi_comm, mds, rc) dtcfl = real(dt_in(2),4) dtcfli = real(dt_in(3),4) dtmin = real(dt_in(4),4) - write(logmsg,'(A,4f10.2)') trim(subname)//': mod_def timesteps reset ',dtmax,dtcfl,dtcfli,dtmin - call ESMF_LogWrite(trim(logmsg), ESMF_LOGMSG_INFO) end if + + ! log info + if (root_task) then + write(stdout,'(a)') trim(logmsg) + write(cvalue,'(4f10.1)')dtmax,dtcfl,dtcfli,dtmin + write(stdout,'(a)') trim(subname)//': WW3 timesteps '//trim(cvalue) + + if (addrstflds) then + do i = 1,rstfldcnt + write(stdout,'(a,i3,a)') trim(subname)//': WW3 additional restart field : ',i,' '//trim(rstfldlist(i)) + end do + else + write(stdout,'(/,a)') trim(subname)//': WW3 NO additional restart fields will be written ' + end if + end if + if (dbug_flag > 5) call ESMF_LogWrite(trim(subname)//' done', ESMF_LOGMSG_INFO) end subroutine waveinit_ufs diff --git a/model/src/wav_grdout.F90 b/model/src/wav_grdout.F90 deleted file mode 100644 index 4583070d7..000000000 --- a/model/src/wav_grdout.F90 +++ /dev/null @@ -1,295 +0,0 @@ -module wav_grdout - - use w3odatmd , only: nogrp, ngrpp - - implicit none - - integer, parameter :: maxvars = 25 ! maximum number of variables/group - - private ! except - - public :: varatts - public :: outvars - public :: wavinit_grdout - - ! tag read from inp file and is used to set flogrd flags - ! var_name is the name of the variable - type :: varatts - character(len= 5) :: tag - character(len=10) :: var_name - character(len=48) :: long_name - character(len=10) :: unit_name - character(len= 2) :: dims - logical :: validout - end type varatts - - type(varatts), dimension(nogrp,maxvars) :: gridoutdefs - - type(varatts), dimension(:), allocatable :: outvars - - !=============================================================================== -contains - !=============================================================================== - - !==================================================================================== - subroutine wavinit_grdout - - use w3gdatmd , only: e3df, p2msf, us3df, usspf - use w3odatmd , only: nds, iaproc, napout - use w3iogomd , only: fldout - use w3servmd , only: strsplit - - ! local variables - character(len=100) :: inptags(100) = '' - integer :: j,k,n,nout - character(len= 12) :: ttag - - ! obtain all possible output variable tags and attributes - call initialize_gridout - - ! obtain the tags for the requested output variables - call strsplit(fldout,inptags) - - ! determine which variables are tagged for output - do k = 1,nogrp - do j = 1,maxvars - if (len_trim(gridoutdefs(k,j)%tag) > 0) then - do n = 1,len(inptags) - if (len_trim(inptags(n)) > 0) then - if (trim(inptags(n)) == trim(gridoutdefs(k,j)%tag)) gridoutdefs(k,j)%validout = .true. - end if - end do - end if - end do - end do - - ! remove requested variables which are only allocated if specific - ! options are set in mod_def (see w3adatmd, '3D arrays') - do k = 1,nogrp - do j = 1,maxvars - if (gridoutdefs(k,j)%validout) then - ttag = trim(gridoutdefs(k,j)%tag) - if (ttag == 'EF' .and. e3df(1,1) == 0) gridoutdefs(k,j)%validout = .false. - if (ttag == 'TH1M' .and. e3df(1,2) == 0) gridoutdefs(k,j)%validout = .false. - if (ttag == 'STH1M' .and. e3df(1,3) == 0) gridoutdefs(k,j)%validout = .false. - if (ttag == 'TH2M' .and. e3df(1,4) == 0) gridoutdefs(k,j)%validout = .false. - if (ttag == 'STH2M' .and. e3df(1,5) == 0) gridoutdefs(k,j)%validout = .false. - - if (ttag == 'P2L' .and. p2msf(1) == 0) gridoutdefs(k,j)%validout = .false. - if (ttag == 'USF' .and. us3df(1) == 0) gridoutdefs(k,j)%validout = .false. - if (ttag == 'USP' .and. usspf(1) == 0) gridoutdefs(k,j)%validout = .false. - end if - end do - end do - - ! determine number of output variables (not the same as the number of tags) - n = 0 - do k = 1,nogrp - do j = 1,maxvars - if (gridoutdefs(k,j)%validout) n = n+1 - end do - end do - nout = n - allocate(outvars(1:nout)) - - ! subset variables requested - n = 0 - do k = 1,nogrp - do j = 1,maxvars - if (gridoutdefs(k,j)%validout) then - n = n+1 - outvars(n) = gridoutdefs(k,j) - end if - enddo - end do - - ! check - if ( iaproc == napout ) then - write(nds(1),*) - write(nds(1),'(a)')' --------------------------------------------------' - write(nds(1),'(a)')' Requested gridded output variables : ' - write(nds(1),'(a)')' --------------------------------------------------' - write(nds(1),*) - do n = 1,nout - write(nds(1),'(i5,2a12,a50)')n,' '//trim(outvars(n)%tag), & - ' '//trim(outvars(n)%var_name), & - ' '//trim(outvars(n)%long_name) - end do - write(nds(1),*) - end if - - end subroutine wavinit_grdout - - !==================================================================================== - subroutine initialize_gridout - - gridoutdefs(:,:)%tag = "" - gridoutdefs(:,:)%var_name = "" - gridoutdefs(:,:)%long_name = "" - gridoutdefs(:,:)%unit_name = "" - gridoutdefs(:,:)%dims = "" - gridoutdefs(:,:)%validout = .false. - - ! TODO: confirm unit values - ! 1 Forcing Fields - gridoutdefs(1,1:14) = [ & - varatts( "DPT ", "DW ", "Water depth ", "m ", " ", .false.) , & - varatts( "CUR ", "CX ", "Mean current, x-component ", "m s-1 ", " ", .false.) , & - varatts( "CUR ", "CY ", "Mean current, y-component ", "m s-1 ", " ", .false.) , & - varatts( "WND ", "UAX ", "Mean wind, x-component ", "m s-1 ", " ", .false.) , & - varatts( "WND ", "UAY ", "Mean wind, y-component ", "m s-1 ", " ", .false.) , & - varatts( "AST ", "AS ", "Air-sea temperature difference ", "K ", " ", .false.) , & - varatts( "WLV ", "WLV ", "Water levels ", "m ", " ", .false.) , & - varatts( "ICE ", "ICE ", "Ice coverage ", "nd ", " ", .false.) , & - varatts( "IBG ", "BERG ", "Iceberg-induced damping ", "km-1 ", " ", .false.) , & - varatts( "TAUA ", "TAUAX ", "Atm momentum x ", "Pa ", " ", .false.) , & - varatts( "TAUA ", "TAUAY ", "Atm momentum y ", "Pa ", " ", .false.) , & - varatts( "RHO ", "RHOAIR ", "Air density ", "kg m-3 ", " ", .false.) , & - varatts( "IC1 ", "ICEH ", "Ice thickness ", "m ", " ", .false.) , & - varatts( "IC5 ", "ICEF ", "Ice floe diameter ", "m ", " ", .false.) & - ] - - ! 2 Standard mean wave Parameters - gridoutdefs(2,1:18) = [ & - varatts( "HS ", "HS ", "Significant wave height ", "m ", " ", .false.) , & - varatts( "LM ", "WLM ", "Mean wave length ", "m ", " ", .false.) , & - varatts( "T02 ", "T02 ", "Mean wave period (Tm0,2) ", "s ", " ", .false.) , & - varatts( "T0M1 ", "T0M1 ", "Mean wave period (Tm0,-1) ", "s ", " ", .false.) , & - varatts( "T01 ", "T01 ", "Mean wave period (Tm0,1) ", "s ", " ", .false.) , & - varatts( "FP ", "FP0 ", "Peak frequency ", "s-1 ", " ", .false.) , & - varatts( "DIR ", "THM ", "Mean wave direction ", "rad ", " ", .false.) , & - varatts( "SPR ", "THS ", "Mean directional spread ", "rad ", " ", .false.) , & - varatts( "DP ", "THP0 ", "Peak direction ", "rad ", " ", .false.) , & - varatts( "HIG ", "HSIG ", "Infragravity height ", "m ", " ", .false.) , & - varatts( "MXE ", "STMAXE ", "Max surface elev (STE) ", "m ", " ", .false.) , & - varatts( "MXES ", "STMAXD ", "St Dev Max surface elev (STE) ", "m ", " ", .false.) , & - varatts( "MXH ", "HMAXE ", "Max wave height (S.) ", "m ", " ", .false.) , & - varatts( "MXHC ", "HCMAXE ", "Max wave height from crest (STE) ", "m ", " ", .false.) , & - varatts( "SDMH ", "HMAXD ", "St Dev of MXC (STE) ", "m ", " ", .false.) , & - varatts( "SDMHC", "HCMAXD ", "St Dev of MXHC (STE) ", "m ", " ", .false.) , & - varatts( "WBT ", "WBT ", "Dominant wave breaking probability (b_T) ", "nd ", " ", .false.) , & - varatts( "WNM ", "WNMEAN ", "Mean wave number ", "m-1 ", " ", .false.) & - ] - - ! 3 Spectral Parameters - gridoutdefs(3,1:6) = [ & - varatts( "EF ", "EF ", "1D spectral density ", "m2 s ", "k ", .false.) , & - varatts( "TH1M ", "TH1M ", "Mean wave direction from a1,b2 ", "deg ", "k ", .false.) , & - 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.) & - ] - - ! 4 Spectral Partition Parameters - gridoutdefs(4,1:17) = [ & - varatts( "PHS ", "PHS ", "Partitioned wave heights ", "m ", "s ", .false.) , & - varatts( "PTP ", "PTP ", "Partitioned peak period ", "s ", "s ", .false.) , & - varatts( "PLP ", "PLP ", "Partitioned peak wave length ", "m ", "s ", .false.) , & - varatts( "PDIR ", "PDIR ", "Partitioned mean direction ", "deg ", "s ", .false.) , & - varatts( "PSPR ", "PSI ", "Partitioned mean directional spread ", "deg ", "s ", .false.) , & - varatts( "PWS ", "PWS ", "Partitioned wind sea fraction ", "nd ", "s ", .false.) , & - varatts( "PDP ", "PTHP0 ", "Peak wave direction of partition ", "deg ", "s ", .false.) , & - varatts( "PQP ", "PQP ", "Goda peakdedness parameter of partition ", "nd ", "s ", .false.) , & - varatts( "PPE ", "PPE ", "JONSWAP peak enhancement factor of partition ", "s-1 ", "s ", .false.) , & - varatts( "PGW ", "PGW ", "Gaussian frequency width of partition ", "nd ", "s ", .false.) , & - varatts( "PSW ", "PSW ", "Spectral width of partition ", "nd ", "s ", .false.) , & - varatts( "PTM10", "PTM1 ", "Mean wave period (m-1,0) of partition ", "s ", "s ", .false.) , & - varatts( "PT01 ", "PT1 ", "Mean wave period (m0,1) of partition ", "s ", "s ", .false.) , & - varatts( "PT02 ", "PT2 ", "Mean wave period (m0,2) of partition ", "s ", "s ", .false.) , & - varatts( "PEP ", "PEP ", "Peak spectral density of partition ", "m2 s rad-1", "s ", .false.) , & - varatts( "TWS ", "PWST ", "Total wind sea fraction ", "nd ", " ", .false.) , & - varatts( "PNR ", "PNR ", "Number of partitions ", "nd ", " ", .false.) & - ] - - ! 5 Atmosphere-waves layer - gridoutdefs(5,1:14) = [ & - 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( "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.) , & - varatts( "TWA ", "TAUWNX ", "Negative part of the wave-supported stress x ", "m2 s-2 ", " ", .false.) , & - varatts( "TWA ", "TAUWNY ", "Negative part of the wave-supported stress y ", "m2 s-2 ", " ", .false.) , & - varatts( "WCC ", "WCC ", "Whitecap coverage ", "nd ", " ", .false.) , & - varatts( "WCF ", "WCF ", "Whitecap foam thickness ", "m ", " ", .false.) , & - varatts( "WCH ", "WCH ", "Mean breaking wave heigh ", "m ", " ", .false.) , & - varatts( "WCM ", "WCM ", "Whitecap moment ", "nd ", " ", .false.) , & - varatts( "FWS ", "TWS ", "Wind sea mean period ", "s ", " ", .false.) & - ] - - ! 6 Wave-ocean layer - gridoutdefs(6,1:25) = [ & - varatts( "SXY ", "SXX ", "Radiation stresses xx ", "N m-1 ", " ", .false.) , & - varatts( "SXY ", "SYY ", "Radiation stresses yy ", "N m-1 ", " ", .false.) , & - varatts( "SXY ", "SXY ", "Radiation stresses xy ", "N m-1 ", " ", .false.) , & - varatts( "TWO ", "TAUOX ", "Wave to ocean momentum flux x ", "m2 s-2 ", " ", .false.) , & - varatts( "TWO ", "TAUOY ", "Wave to ocean momentum flux y ", "m2 s-2 ", " ", .false.) , & - varatts( "BHD ", "BHD ", "Bernoulli head (J term) ", "m2 s-2 ", " ", .false.) , & - varatts( "FOC ", "PHIOC ", "Wave to ocean energy flux ", "W m-2 ", " ", .false.) , & - varatts( "TUS ", "TUSX ", "Stokes transport x ", "m2 s-1 ", " ", .false.) , & - varatts( "TUS ", "TUSY ", "Stokes transport y ", "m2 s-1 ", " ", .false.) , & - varatts( "USS ", "USSX ", "Surface Stokes drift x ", "m s-1 ", " ", .false.) , & - varatts( "USS ", "USSY ", "Surface Stokes drift y ", "m s-1 ", " ", .false.) , & - varatts( "P2S ", "PRMS ", "Second-order sum pressure ", "m4 ", " ", .false.) , & - varatts( "P2S ", "TPMS ", "Second-order sum pressure ", "s-1 ", " ", .false.) , & - varatts( "USF ", "US3DX ", "Spectrum of surface Stokes drift x ", "m s-1 Hz-1", "k ", .false.) , & - varatts( "USF ", "US3DY ", "Spectrum of surface Stokes drift y ", "m s-1 Hz-1", "k ", .false.) , & - varatts( "P2L ", "P2SMS ", "Micro seism source term ", "Pa2 m2 s ", "m ", .false.) , & - varatts( "TWI ", "TAUICEX ", "Wave to sea ice stress x ", "m2 s-2 ", " ", .false.) , & - varatts( "TWI ", "TAUICEY ", "Wave to sea ice stress y ", "m2 s-2 ", " ", .false.) , & - varatts( "FIC ", "PHICE ", "Wave to sea ice energy flux ", "W m-2 ", " ", .false.) , & - varatts( "USP ", "USSPX ", "Partitioned surface Stokes drift x ", "m s-1 ", "p ", .false.) , & - varatts( "USP ", "USSPY ", "Partitioned surface Stokes drift y ", "m s-1 ", "p ", .false.) , & - varatts( "TWC ", "TAUOCX ", "Total wave to ocean stress x ", "Pa ", " ", .false.) , & - varatts( "TWC ", "TAUOCY ", "Total wave to ocean stress y ", "Pa ", " ", .false.) , & - varatts( "USSH ", "USSHX ", "Surface layer averaged Stokes drift x ", "m s-1 ", " ", .false.) , & - varatts( "USSH ", "USSHY ", "Surface layer averaged Stokes drift y ", "m s-1 ", " ", .false.) & - ] - - ! 7 Wave-bottom layer - gridoutdefs(7,1:10) = [ & - varatts( "ABR ", "ABAX ", "Near bottom rms wave excursion amplitudes x ", "m ", " ", .false.) , & - varatts( "ABR ", "ABAY ", "Near bottom rms wave excursion amplitudes y ", "m ", " ", .false.) , & - varatts( "UBR ", "UBAX ", "Near bottom rms wave velocities x ", "m s-1 ", " ", .false.) , & - varatts( "UBR ", "UBAY ", "Near bottom rms wave velocities y ", "m s-1 ", " ", .false.) , & - varatts( "BED ", "BED ", "Bottom roughness ", "m ", " ", .false.) , & - varatts( "BED ", "RIPPLEX ", "Sea bottom ripple wavelength x ", "m ", " ", .false.) , & - varatts( "BED ", "RIPPLEY ", "Sea bottom ripple wavelength y ", "m ", " ", .false.) , & - varatts( "FBB ", "PHIBBL ", "Energy flux due to bottom friction ", "W m-2 ", " ", .false.) , & - varatts( "TBB ", "TAUBBLX ", "Momentum flux due to bottom friction x ", "m2 s-2 ", " ", .false.) , & - varatts( "TBB ", "TAUBBLY ", "Momentum flux due to bottom friction y ", "m2 s-2 ", " ", .false.) & - ] - - ! 8 Spectrum parameters - gridoutdefs(8,1:9) = [ & - varatts( "MSS ", "MSSX ", "Surface mean square slope x ", "nd ", " ", .false.) , & - varatts( "MSS ", "MSSY ", "Surface mean square slope y ", "nd ", " ", .false.) , & - varatts( "MSC ", "MSCX ", "Spectral level at high frequency tail x ", "nd ", " ", .false.) , & - varatts( "MSC ", "MSCY ", "Spectral level at high frequency tail y ", "nd ", " ", .false.) , & - varatts( "WL02 ", "WL02X ", "East/X North/Y mean wavelength component ", "nd ", " ", .false.) , & - varatts( "WL02 ", "WL02Y ", "East/X North/Y mean wavelength component ", "nd ", " ", .false.) , & - varatts( "AXT ", "ALPXT ", "Correl sea surface gradients (x,t) ", "nd ", " ", .false.) , & - varatts( "AYT ", "ALPYT ", "Correl sea surface gradients (y,t) ", "nd ", " ", .false.) , & - varatts( "AXY ", "ALPXY ", "Correl sea surface gradients (x,y) ", "nd ", " ", .false.) & - ] - - ! 9 Numerical diagnostics - gridoutdefs(9,1:5) = [ & - varatts( "DTD ", "DTDYN ", "Average time step in integration ", "min ", " ", .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.) , & - varatts( "CFK ", "CFLKMAX ", "Max. CFL number for k-advection ", "nd ", " ", .false.) & - ] - - ! 10 User defined - gridoutdefs(10,1:2) = [ & - varatts( "U1 ", "U1 ", "User defined 1 ", "nd ", " ", .false.) , & - varatts( "U2 ", "U2 ", "User defined 2 ", "nd ", " ", .false.) & - ] - end subroutine initialize_gridout -end module wav_grdout diff --git a/model/src/wav_history_mod.F90 b/model/src/wav_history_mod.F90 new file mode 100644 index 000000000..e7a3e71a9 --- /dev/null +++ b/model/src/wav_history_mod.F90 @@ -0,0 +1,914 @@ +!> @file wav_history_mod +!! +!> @brief Manage gridded model output as netCDF using PIO +!! +!> @author mvertens@ucar.edu, Denise.Worthen@noaa.gov +!> @date 01-05-2022 +module wav_history_mod + + use constants , only : rade + use w3parall , only : init_get_isea + use w3gdatmd , only : xgrd, ygrd + use w3gdatmd , only : nk, nx, ny, mapsf, mapsta, nsea + use w3odatmd , only : undef + use w3adatmd , only : mpi_comm_wave + use wav_import_export , only : nseal_cpl + use wav_pio_mod , only : pio_iotype, pio_ioformat, wav_pio_subsystem + use wav_pio_mod , only : handle_err, wav_pio_initdecomp + use pio + use netcdf + + implicit none + + private + + public :: wav_history_init + public :: write_history + public :: varatts + public :: outvars + + ! used/reused in module + integer :: isea, jsea, ix, iy, ierr + + real, allocatable, target :: var3ds(:,:) + real, allocatable, target :: var3dm(:,:) + real, allocatable, target :: var3dp(:,:) + real, allocatable, target :: var3dk(:,:) + + ! output variable for (nx,ny,nz) fields + real, pointer :: var3d(:,:) + + type(file_desc_t) :: pioid + type(var_desc_t) :: varid + type(io_desc_t) :: iodesc2d !2d only + type(io_desc_t) :: iodesc2dint !2d only, integer + type(io_desc_t) :: iodesc3ds !s-axis variables + type(io_desc_t) :: iodesc3dm !m-axis variables + type(io_desc_t) :: iodesc3dp !p-axis variables + type(io_desc_t) :: iodesc3dk !k-axis variables + + ! variable attributes + type :: varatts + character(len= 5) :: tag + character(len=10) :: var_name + character(len=48) :: long_name + character(len=10) :: unit_name + character(len= 2) :: dims + logical :: validout + end type varatts + + type(varatts), dimension(:), allocatable :: outvars + + !=============================================================================== +contains + !=============================================================================== + !> Write the requested list of fields using parallel netCDF via PIO + !! + !! @param[in] timen the timestamp for the file + !! + !> author DeniseWorthen@noaa.gov + !> @date 08-26-2024 + subroutine write_history ( timen ) + + use w3odatmd , only : fnmpre + use w3gdatmd , only : filext, trigp, ntri, ungtype, gtype + use w3servmd , only : extcde + use w3wdatmd , only : wlv, ice, icef, iceh, berg, ust, ustdir, asf, rhoair + use w3gdatmd , only : e3df, p2msf, us3df, usspf + use w3odatmd , only : noswll + use w3odatmd , only : ndso, iaproc + use w3adatmd , only : dw, ua, ud, as, cx, cy, taua, tauadir + use w3adatmd , only : hs, wlm, t02, t0m1, t01, fp0, thm, ths, thp0, wbt, wnmean + use w3adatmd , only : dtdyn + use w3adatmd , only : fcut, aba, abd, uba, ubd, sxx, syy, sxy + use w3adatmd , only : phs, ptp, plp, pdir, psi, pws, pwst, pnr + use w3adatmd , only : pthp0, pqp, ppe, pgw, psw, ptm1, pt1, pt2 + use w3adatmd , only : pep, tauox, tauoy, tauwix, tauwiy + use w3adatmd , only : phiaw, phioc, tusx, tusy, prms, tpms + use w3adatmd , only : ussx, ussy, mssx, mssy, mscx, mscy + use w3adatmd , only : tauwnx, tauwny, charn, tws, bhd + use w3adatmd , only : phibbl, taubbl, whitecap, bedforms, cge, ef + use w3adatmd , only : cflxymax, cflthmax, cflkmax, p2sms, us3d + use w3adatmd , only : hsig, phice, tauice + use w3adatmd , only : stmaxe, stmaxd, hmaxe, hcmaxe, hmaxd, hcmaxd, ussp, tauocx, tauocy + use w3adatmd , only : usshx, usshy + + use w3timemd , only : set_user_timestring + use w3odatmd , only : time_origin, calendar_name, elapsed_secs + use w3odatmd , only : user_histfname + !TODO: use unstr_mesh from wav_shr_mod; currently fails due to CI + !use wav_shr_mod , only : unstr_mesh + + integer, intent(in) :: timen(2) + + ! local variables + integer ,target :: dimid3(3) + integer ,target :: dimid4(4) + integer ,pointer :: dimid(:) + character(len=1024) :: fname + character(len=12) :: vname + character(len=16) :: user_timestring !YYYY-MM-DD-SSSSS + + 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 :: lmap(nseal_cpl) + + ! ------------------------------------------------------------- + ! create the netcdf file + ! ------------------------------------------------------------- + + ! native WW3 file naming + if (len_trim(user_histfname) == 0) then + write(fname,'(a,i8.8,a1,i6.6,a)')trim(fnmpre),timen(1),'.',timen(2),'.out_grd.ww3.nc' + else + call set_user_timestring(timen,user_timestring) + fname = trim(user_histfname)//trim(user_timestring)//'.nc' + end if + + pioid%fh = -1 + nmode = pio_clobber + ! only applies to classic NETCDF files. + if (pio_iotype == PIO_IOTYPE_NETCDF .or. pio_iotype == PIO_IOTYPE_PNETCDF) then + nmode = ior(nmode,pio_ioformat) + endif + ierr = pio_createfile(wav_pio_subsystem, pioid, pio_iotype, trim(fname), nmode) + call handle_err(ierr, 'pio_create') + if (iaproc == 1) write(ndso,'(a)')' Writing history file '//trim(fname) + + len_s = noswll + 1 ! 0:noswll + len_m = p2msf(3)-p2msf(2) + 1 ! ? + len_p = usspf(2) ! partitions + len_k = e3df(3,1) - e3df(2,1) + 1 ! frequencies + + ! define the dimensions required for the requested gridded fields + do n = 1,size(outvars) + if (outvars(n)%validout) then + if(trim(outvars(n)%dims) == 's')s_axis = .true. + 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 + 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)) + + ierr = pio_def_dim(pioid, 'nx', nx, xtid) + ierr = pio_def_dim(pioid, 'ny', ny, ytid) + ierr = pio_def_dim(pioid, 'time', PIO_UNLIMITED, timid) + + 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 (gtype .eq. ungtype) then + ierr = pio_def_dim(pioid, 'ne' , ntri, xeid) + ierr = pio_def_dim(pioid, 'nn' , 3, ztid) + end if + + ! define the time variable + ierr = pio_def_var(pioid, 'time', PIO_DOUBLE, (/timid/), varid) + call handle_err(ierr,'def_timevar') + ierr = pio_put_att(pioid, varid, 'units', trim(time_origin)) + call handle_err(ierr,'def_time_units') + ierr = pio_put_att(pioid, varid, 'calendar', trim(calendar_name)) + call handle_err(ierr,'def_time_calendar') + + ! define the spatial axis variables (lat,lon) + ierr = pio_def_var(pioid, 'lon', PIO_DOUBLE, (/xtid,ytid/), varid) + call handle_err(ierr,'def_lonvar') + ierr = pio_put_att(pioid, varid, 'units', 'degrees_east') + ierr = pio_def_var(pioid, 'lat', PIO_DOUBLE, (/xtid,ytid/), varid) + call handle_err(ierr,'def_latvar') + ierr = pio_put_att(pioid, varid, 'units', 'degrees_north') + + ! add mapsta + ierr = pio_def_var(pioid, 'mapsta', PIO_INT, (/xtid, ytid, timid/), varid) + call handle_err(ierr, 'def_mapsta') + ierr = pio_put_att(pioid, varid, 'units', 'unitless') + ierr = pio_put_att(pioid, varid, 'long_name', 'map status') + ierr = pio_put_att(pioid, varid, '_FillValue', nf90_fill_int) + + if (gtype .eq. ungtype) then + ierr = pio_def_var(pioid, 'nconn', PIO_INT, (/ztid,xeid/), varid) + call handle_err(ierr,'def_nodeconnections') + ierr = pio_put_att(pioid, varid, 'units', 'unitless') + ierr = pio_put_att(pioid, varid, 'long_name', 'node connectivity') + end if + + ! define the variables + dimid3(1:2) = (/xtid, ytid/) + dimid4(1:2) = (/xtid, ytid/) + do n = 1,size(outvars) + if (trim(outvars(n)%dims) == 's') then + dimid4(3:4) = (/stid, timid/) + dimid => dimid4 + else if (trim(outvars(n)%dims) == 'm') then + dimid4(3:4) = (/mtid, timid/) + dimid => dimid4 + else if (trim(outvars(n)%dims) == 'p') then + dimid4(3:4) = (/ptid, timid/) + dimid => dimid4 + else if (trim(outvars(n)%dims) == 'k') then + dimid4(3:4) = (/ktid, timid/) + dimid => dimid4 + else + dimid3(3) = timid + dimid => dimid3 + end if + + ierr = pio_def_var(pioid, trim(outvars(n)%var_name), PIO_REAL, dimid, varid) + call handle_err(ierr, 'define variable '//trim((outvars(n)%var_name))) + ierr = pio_put_att(pioid, varid, 'units' , trim(outvars(n)%unit_name)) + ierr = pio_put_att(pioid, varid, 'long_name' , trim(outvars(n)%long_name)) + ierr = pio_put_att(pioid, varid, '_FillValue', undef) + end do + ! end variable definitions + ierr = pio_enddef(pioid) + call handle_err(ierr, 'end variable definition') + + call wav_pio_initdecomp(iodesc2d) + call wav_pio_initdecomp(iodesc2dint, use_int=.true.) + if (s_axis)call wav_pio_initdecomp(len_s, iodesc3ds) + 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) + + ! write the time and spatial axis values (lat,lon,time) + ierr = pio_inq_varid(pioid, 'lat', varid) + call handle_err(ierr, 'inquire variable lat ') + ierr = pio_put_var(pioid, varid, transpose(ygrd)) + call handle_err(ierr, 'put lat') + + ierr = pio_inq_varid(pioid, 'lon', varid) + call handle_err(ierr, 'inquire variable lon ') + ierr = pio_put_var(pioid, varid, transpose(xgrd)) + call handle_err(ierr, 'put lon') + + ierr = pio_inq_varid(pioid, 'time', varid) + 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 (gtype .eq. ungtype) then + ierr = pio_inq_varid(pioid, 'nconn', varid) + call handle_err(ierr, 'inquire variable nconn ') + ierr = pio_put_var(pioid, varid, trigp) + call handle_err(ierr, 'put trigp') + end if + + ! TODO: tried init decomp w/ use_int=.true. but getting garbage + ! land values....sea values OK + ! mapsta is global + lmap(:) = 0 + do jsea = 1,nseal_cpl + call init_get_isea(isea, jsea) + ix = mapsf(isea,1) + iy = mapsf(isea,2) + lmap(jsea) = mapsta(iy,ix) + end do + ierr = pio_inq_varid(pioid, 'mapsta', varid) + call handle_err(ierr, 'inquire variable mapsta ') + call pio_setframe(pioid, varid, int(1,kind=Pio_Offset_Kind)) + call pio_write_darray(pioid, varid, iodesc2dint, lmap, ierr) + call handle_err(ierr, 'put variable mapsta') + + ! write the requested variables + do n = 1,size(outvars) + 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. '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) ) + if(vname .eq. 'PPE') call write_var3d(iodesc3ds, vname, ppe (1:nseal_cpl,0:noswll) ) + if(vname .eq. 'PGW') call write_var3d(iodesc3ds, vname, pgw (1:nseal_cpl,0:noswll) ) + if(vname .eq. 'PSW') call write_var3d(iodesc3ds, vname, psw (1:nseal_cpl,0:noswll) ) + if(vname .eq. 'PTM1') call write_var3d(iodesc3ds, vname, ptm1 (1:nseal_cpl,0:noswll) ) + if(vname .eq. 'PT1') call write_var3d(iodesc3ds, vname, pt1 (1:nseal_cpl,0:noswll) ) + if(vname .eq. 'PT2') call write_var3d(iodesc3ds, vname, pt2 (1:nseal_cpl,0:noswll) ) + if(vname .eq. 'PEP') call write_var3d(iodesc3ds, vname, pep (1:nseal_cpl,0:noswll) ) + + else if (trim(outvars(n)%dims) == 'm') then ! m axis + var3d => var3dm + ! Group 6 + if (vname .eq. 'P2SMS') call write_var3d(iodesc3dm, vname, p2sms (1:nseal_cpl,p2msf(2):p2msf(3)) ) + + else if (trim(outvars(n)%dims) == 'p') then ! partition axis + var3d => var3dp + ! Group 6 + 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) == '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. '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)) ) + if(vname .eq. 'STH2M') call write_var3d(iodesc3dk, vname, ef (1:nseal_cpl,e3df(2,5):e3df(3,5)) ) + !TODO: wn has reversed indices (1:nk, 1:nseal_cpl) + ! Group 6 + if (vname .eq. 'US3DX') call write_var3d(iodesc3dk, vname, us3d (1:nseal_cpl, us3df(2):us3df(3)) ) + if (vname .eq. 'US3DY') call write_var3d(iodesc3dk, vname, us3d (1:nseal_cpl,nk+us3df(2):nk+us3df(3)) ) + + else + ! Group 1 + if (vname .eq. 'DW') call write_var2d(vname, dw (1:nsea), init0='false', global='true') + if (vname .eq. 'CX') call write_var2d(vname, cx (1:nsea), init0='false', global='true') + if (vname .eq. 'CY') call write_var2d(vname, cy (1:nsea), init0='false', global='true') + if (vname .eq. 'UAX') call write_var2d(vname, ua (1:nsea), dir=cos(ud(1:nsea)), init0='false', global='true') + if (vname .eq. 'UAY') call write_var2d(vname, ua (1:nsea), dir=sin(ud(1:nsea)), init0='false', global='true') + if (vname .eq. 'AS') call write_var2d(vname, as (1:nsea), init0='false', global='true') + if (vname .eq. 'WLV') call write_var2d(vname, wlv (1:nsea), init0='false', global='true') + if (vname .eq. 'ICE') call write_var2d(vname, ice (1:nsea), init0='false', global='true') + if (vname .eq. 'IBG') call write_var2d(vname, berg (1:nsea), init0='false', global='true') + if (vname .eq. 'TAUX') call write_var2d(vname, taua (1:nsea), dir=cos(tauadir(1:nsea)), init0='false', global='true') + if (vname .eq. 'TAUY') call write_var2d(vname, taua (1:nsea), dir=sin(tauadir(1:nsea)), init0='false', global='true') + if (vname .eq. 'RHOAIR') call write_var2d(vname, rhoair (1:nsea), init0='false', global='true') + if (vname .eq. 'ICEH') call write_var2d(vname, iceh (1:nsea), init0='false', global='true') + if (vname .eq. 'ICEF') call write_var2d(vname, icef (1:nsea), init0='false', global='true') + + ! Group 2 + if (vname .eq. 'HS') call write_var2d(vname, hs (1:nseal_cpl) ) + if (vname .eq. 'WLM') call write_var2d(vname, wlm (1:nseal_cpl) ) + if (vname .eq. 'T02') call write_var2d(vname, t02 (1:nseal_cpl) ) + if (vname .eq. 'T0M1') call write_var2d(vname, t0m1 (1:nseal_cpl) ) + 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. '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) ) + if (vname .eq. 'STMAXD') call write_var2d(vname, stmaxd (1:nseal_cpl) ) + if (vname .eq. 'HMAXE') call write_var2d(vname, hmaxe (1:nseal_cpl) ) + if (vname .eq. 'HCMAXE') call write_var2d(vname, hcmaxe (1:nseal_cpl) ) + if (vname .eq. 'HMAXD') call write_var2d(vname, hmaxd (1:nseal_cpl) ) + if (vname .eq. 'HCMAXD') call write_var2d(vname, hcmaxd (1:nseal_cpl) ) + if (vname .eq. 'WBT') call write_var2d(vname, wbt (1:nseal_cpl) ) + if (vname .eq. 'WNMEAN') call write_var2d(vname, wnmean (1:nseal_cpl), init0='false') + + ! Group 4 + if(vname .eq. 'PWST') call write_var2d(vname, pwst (1:nseal_cpl) ) + 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. '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') + if (vname .eq. 'TAUWIX') call write_var2d(vname, tauwix (1:nseal_cpl), init2='true') + if (vname .eq. 'TAUWIY') call write_var2d(vname, tauwiy (1:nseal_cpl), init2='true') + if (vname .eq. 'TAUWNX') call write_var2d(vname, tauwnx (1:nseal_cpl), init2='true') + if (vname .eq. 'TAUWNY') call write_var2d(vname, tauwny (1:nseal_cpl), init2='true') + if (vname .eq. 'WCC') call write_var2d(vname, whitecap (1:nseal_cpl,1), init2='true') + if (vname .eq. 'WCF') call write_var2d(vname, whitecap (1:nseal_cpl,2), init2='true') + if (vname .eq. 'WCH') call write_var2d(vname, whitecap (1:nseal_cpl,3), init2='true') + if (vname .eq. 'WCM') call write_var2d(vname, whitecap (1:nseal_cpl,4), init2='true') + if (vname .eq. 'TWS') call write_var2d(vname, tws (1:nseal_cpl) ) + + ! Group 6 + if (vname .eq. 'SXX') call write_var2d(vname, sxx (1:nseal_cpl) ) + if (vname .eq. 'SYY') call write_var2d(vname, syy (1:nseal_cpl) ) + if (vname .eq. 'SXY') call write_var2d(vname, sxy (1:nseal_cpl) ) + if (vname .eq. 'TAUOX') call write_var2d(vname, tauox (1:nseal_cpl), init2='true') + if (vname .eq. 'TAUOY') call write_var2d(vname, tauoy (1:nseal_cpl), init2='true') + if (vname .eq. 'BHD') call write_var2d(vname, bhd (1:nseal_cpl) ) + if (vname .eq. 'PHIOC') call write_var2d(vname, phioc (1:nseal_cpl), init2='true') + if (vname .eq. 'TUSX') call write_var2d(vname, tusx (1:nseal_cpl) ) + if (vname .eq. 'TUSY') call write_var2d(vname, tusy (1:nseal_cpl) ) + if (vname .eq. 'USSX') call write_var2d(vname, ussx (1:nseal_cpl) ) + if (vname .eq. 'USSY') call write_var2d(vname, ussy (1:nseal_cpl) ) + if (vname .eq. 'PRMS') call write_var2d(vname, prms (1:nseal_cpl) ) + if (vname .eq. 'TPMS') call write_var2d(vname, tpms (1:nseal_cpl) ) + if (vname .eq. 'TAUICEX') call write_var2d(vname, tauice (1:nseal_cpl,1) ) + if (vname .eq. 'TAUICEY') call write_var2d(vname, tauice (1:nseal_cpl,2) ) + if (vname .eq. 'PHICE') call write_var2d(vname, phice (1:nseal_cpl) ) + if (vname .eq. 'TAUOCX') call write_var2d(vname, tauocx (1:nseal_cpl) ) + if (vname .eq. 'TAUOCY') call write_var2d(vname, tauocy (1:nseal_cpl) ) + if (vname .eq. 'USSHX') call write_var2d(vname, usshx (1:nseal_cpl) ) + if (vname .eq. 'USSHY') call write_var2d(vname, usshy (1:nseal_cpl) ) + ! Group 7 + if (vname .eq. 'ABAX') call write_var2d(vname, aba (1:nseal_cpl), dir=cos(abd(1:nseal_cpl)) ) + if (vname .eq. 'ABAY') call write_var2d(vname, aba (1:nseal_cpl), dir=sin(abd(1:nseal_cpl)) ) + if (vname .eq. 'UBAX') call write_var2d(vname, uba (1:nseal_cpl), dir=cos(ubd(1:nseal_cpl)) ) + if (vname .eq. 'UBAY') call write_var2d(vname, uba (1:nseal_cpl), dir=sin(ubd(1:nseal_cpl)) ) + if (vname .eq. 'BED') call write_var2d(vname, bedforms (1:nseal_cpl,1), init2='true') + if (vname .eq. 'RIPPLEX') call write_var2d(vname, bedforms (1:nseal_cpl,2), init2='true') + if (vname .eq. 'RIPPLEY') call write_var2d(vname, bedforms (1:nseal_cpl,3), init2='true') + if (vname .eq. 'PHIBBL') call write_var2d(vname, phibbl (1:nseal_cpl), init2='true') + if (vname .eq. 'TAUBBLX') call write_var2d(vname, taubbl (1:nseal_cpl,1), init2='true') + if (vname .eq. 'TAUBBLY') call write_var2d(vname, taubbl (1:nseal_cpl,2), init2='true') + + ! Group 8 + if (vname .eq. 'MSSX') call write_var2d(vname, mssx (1:nseal_cpl) ) + if (vname .eq. 'MSSY') call write_var2d(vname, mssy (1:nseal_cpl) ) + if (vname .eq. 'MSCX') call write_var2d(vname, mscx (1:nseal_cpl) ) + if (vname .eq. 'MSCY') call write_var2d(vname, mscy (1:nseal_cpl) ) + !TODO: remaining variables have inconsistency between shel_inp listing and iogo code + + ! Group 9 + if (vname .eq. 'DTDYN') call write_var2d(vname, dtdyn (1:nseal_cpl) ) + if (vname .eq. 'FCUT') call write_var2d(vname, fcut (1:nseal_cpl) ) + if (vname .eq.'CFLXYMAX') call write_var2d(vname, cflxymax (1:nseal_cpl) ) + if (vname .eq.'CFLTHMAX') call write_var2d(vname, cflthmax (1:nseal_cpl) ) + if (vname .eq. 'CFLKMAX') call write_var2d(vname, cflkmax (1:nseal_cpl) ) + + ! Group 10 + end if + end do + + if (s_axis) deallocate(var3ds) + if (m_axis) deallocate(var3dm) + if (p_axis) deallocate(var3dp) + if (k_axis) deallocate(var3dk) + + 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) + + call pio_closefile(pioid) + + end subroutine write_history + + !=============================================================================== + !> Write an array of (nseal) points as (nx,ny) + !! + !! @details If dir is present, the written variable will represent either the X + !! or Y component of the variable. If mask is present and true, use mapsta=1 to + !! mask values. If init0 is present and false, do not initialize values for mapsta<0. + !! This prevents group 1 variables being set undef over ice. If init2 is present and + !! true, apply a second initialization where mapsta==2. If fldir is present and true + !! then the directions will be converted to degrees. If global is present and true, + !! write pe-local copy of global field + !! + !! @param[in] vname the variable name + !! @param[in] var the variable array + !! @param[in] dir the direction array, optional + !! @param[in] usemask a flag for variable masking, optional + !! @param[in] init0 a flag for variable initialization, optional + !! @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 DeniseWorthen@noaa.gov + !> @date 08-26-2024 + subroutine write_var2d(vname, var, dir, usemask, init0, init2, fldir, global) + + character(len=*), intent(in) :: vname + real , intent(in) :: var(:) + real, optional , intent(in) :: dir(:) + character(len=*), optional, intent(in) :: usemask + character(len=*), optional, intent(in) :: init0 + character(len=*), optional, intent(in) :: init2 + character(len=*), optional, intent(in) :: fldir + character(len=*), optional, intent(in) :: global + + ! local variables + real, dimension(nseal_cpl) :: varout + logical :: lmask, linit0, linit2, lfldir, lglobal + real :: varloc + + lmask = .false. + if (present(usemask)) then + lmask = (trim(usemask) == "true") + end if + linit0 = .true. + if (present(init0)) then + linit0 = (trim(init0) == "true") + end if + linit2 = .false. + if (present(init2)) then + linit2 = (trim(init2) == "true") + end if + lfldir = .false. + if (present(fldir)) then + lfldir = (trim(fldir) == "true") + end if + lglobal = .false. + if (present(global)) then + lglobal = (trim(global) == "true") + end if + + varout = undef + do jsea = 1,nseal_cpl + call init_get_isea(isea, jsea) + if (lglobal) then + varloc = var(isea) + else + varloc = var(jsea) + end if + + if (linit0) then + if (mapsta(mapsf(isea,2),mapsf(isea,1)) < 0) varloc = undef + end if + if (linit2) then + if (mapsta(mapsf(isea,2),mapsf(isea,1)) == 2) varloc = undef + end if + + if (lfldir) then + if (varloc .ne. undef) then + varloc = mod(630. - rade*varloc, 360.) + end if + end if + if (present(dir)) then + if (varloc .ne. undef) then + if (lmask) then + if (mapsta(mapsf(isea,2),mapsf(isea,1)) == 1) then + if (lglobal) then + varout(jsea) = varloc*dir(isea) + else + varout(jsea) = varloc*dir(jsea) + end if + end if + else + if (lglobal) then + varout(jsea) = varloc*dir(isea) + else + varout(jsea) = varloc*dir(jsea) + end if + end if + end if + else + varout(jsea) = varloc + end if + end do + + ierr = pio_inq_varid(pioid, trim(vname), varid) + call handle_err(ierr, 'inquire variable '//trim(vname)) + call pio_setframe(pioid, varid, int(1,kind=Pio_Offset_Kind)) + call pio_write_darray(pioid, varid, iodesc2d, varout, ierr) + call handle_err(ierr, 'put variable '//trim(vname)) + + end subroutine write_var2d + + !=============================================================================== + !> Write an array of (nseal,:) points as (nx,ny,:) + !! + !! @details If init2 is present and true, apply a second initialization to a + !! subset of variables for where mapsta==2. If fldir is present and true then + !! the directions will be converted to degrees. + !! + !! @param[in] iodesc the PIO decomposition handle + !! @param[in] vname the variable name + !! @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 + !! + !> author DeniseWorthen@noaa.gov + !> @date 08-26-2024 + subroutine write_var3d(iodesc, vname, var, init2, fldir) + + 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 + + ! local variables + real, allocatable, dimension(:) :: varloc + logical :: linit2, lfldir + integer :: lb, ub + + linit2 = .false. + if (present(init2)) then + linit2 = (trim(init2) == "true") + end if + lfldir = .false. + if (present(fldir)) then + lfldir = (trim(fldir) == "true") + end if + + lb = lbound(var,2) + ub = ubound(var,2) + allocate(varloc(lb:ub)) + + var3d = undef + do jsea = 1,nseal_cpl + call init_get_isea(isea, jsea) + ! initialization + varloc(:) = var(jsea,:) + 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.) + end if + end if + var3d(jsea,:) = varloc(:) + end do + + ierr = pio_inq_varid(pioid, trim(vname), varid) + call handle_err(ierr, 'inquire variable '//trim(vname)) + call pio_setframe(pioid, varid, int(1,kind=PIO_OFFSET_KIND)) + call pio_write_darray(pioid, varid, iodesc, var3d, ierr) + + deallocate(varloc) + end subroutine write_var3d + + !=============================================================================== + !> Scan through all possible fields to determine a list of requested variables + !! + !> @details Utilizes the list of variables specified via WW3's flout array + !! to define the variables written to the file. The elements of the flout + !! array are here referred to as "tags". Tags are not 1:1 with output fields. + !! For example, the tag "CUR" means that the ocean currents, comprising two + !! directional values are requested. They will both be requested via the single + !! CUR tag. The tag "SXY" means that three components of radiation stresses + !! are requested (XX,YY,XY). + !! + !! @param[in] stdout the logfile unit on the root_task + !! + !> @author Denise.Worthen@noaa.gov + !> @date 09-19-2022 + subroutine wav_history_init(stdout) + + use w3gdatmd, only: e3df, p2msf, us3df, usspf + use w3odatmd, only: iaproc, nogrp, ngrpp + use w3iogomd, only: fldout + use w3servmd, only: strsplit + + integer, intent(in) :: stdout + + ! local variables + integer, parameter :: maxvars = 25 ! maximum number of variables/group + type(varatts), dimension(nogrp,maxvars) :: gridoutdefs + + character(len=100) :: inptags(100) = '' + integer :: j,k,n,nout + character(len= 12) :: ttag + + ! obtain all possible output variable tags and attributes + call define_fields(gridoutdefs) + + ! obtain the tags for the requested output variables + call strsplit(fldout,inptags) + + ! determine which variables are tagged for output + do k = 1,nogrp + do j = 1,maxvars + if (len_trim(gridoutdefs(k,j)%tag) > 0) then + do n = 1,len(inptags) + if (len_trim(inptags(n)) > 0) then + if (trim(inptags(n)) == trim(gridoutdefs(k,j)%tag)) gridoutdefs(k,j)%validout = .true. + end if + end do + end if + end do + end do + + ! remove requested variables which are only allocated if specific + ! options are set in mod_def (see w3adatmd, '3D arrays') + do k = 1,nogrp + do j = 1,maxvars + if (gridoutdefs(k,j)%validout) then + ttag = trim(gridoutdefs(k,j)%tag) + if (ttag == 'EF' .and. e3df(1,1) == 0) gridoutdefs(k,j)%validout = .false. + if (ttag == 'TH1M' .and. e3df(1,2) == 0) gridoutdefs(k,j)%validout = .false. + if (ttag == 'STH1M' .and. e3df(1,3) == 0) gridoutdefs(k,j)%validout = .false. + if (ttag == 'TH2M' .and. e3df(1,4) == 0) gridoutdefs(k,j)%validout = .false. + if (ttag == 'STH2M' .and. e3df(1,5) == 0) gridoutdefs(k,j)%validout = .false. + + if (ttag == 'P2L' .and. p2msf(1) == 0) gridoutdefs(k,j)%validout = .false. + if (ttag == 'USF' .and. us3df(1) == 0) gridoutdefs(k,j)%validout = .false. + if (ttag == 'USP' .and. usspf(1) == 0) gridoutdefs(k,j)%validout = .false. + end if + end do + end do + + ! determine number of output variables (not the same as the number of tags) + n = 0 + do k = 1,nogrp + do j = 1,maxvars + if (gridoutdefs(k,j)%validout) n = n+1 + end do + end do + nout = n + allocate(outvars(1:nout)) + + ! subset variables requested + n = 0 + do k = 1,nogrp + do j = 1,maxvars + if (gridoutdefs(k,j)%validout) then + n = n+1 + outvars(n) = gridoutdefs(k,j) + end if + enddo + end do + + ! check + if ( iaproc == 1 ) then + write(stdout,*) + write(stdout,'(a)')' --------------------------------------------------' + write(stdout,'(a)')' Requested gridded output variables : ' + write(stdout,'(a)')' --------------------------------------------------' + write(stdout,*) + do n = 1,nout + write(stdout,'(i5,2a12,a50)')n,' '//trim(outvars(n)%tag), & + ' '//trim(outvars(n)%var_name), & + ' '//trim(outvars(n)%long_name) + end do + write(stdout,*) + call flush (stdout) + end if + + end subroutine wav_history_init + + !==================================================================================== + !> Define the available output fields and their attributes + !! + !> @author Denise.Worthen@noaa.gov + !> @date 09-19-2022 + subroutine define_fields (gridoutdefs) + + type(varatts), dimension(:,:), intent(inout) :: gridoutdefs + + gridoutdefs(:,:)%tag = "" + gridoutdefs(:,:)%var_name = "" + gridoutdefs(:,:)%long_name = "" + gridoutdefs(:,:)%unit_name = "" + gridoutdefs(:,:)%dims = "" + gridoutdefs(:,:)%validout = .false. + + ! 1 Forcing Fields + gridoutdefs(1,1:14) = [ & + varatts( "DPT ", "DW ", "Water depth ", "m ", " ", .false.) , & + varatts( "CUR ", "CX ", "Mean current, x-component ", "m s-1 ", " ", .false.) , & + varatts( "CUR ", "CY ", "Mean current, y-component ", "m s-1 ", " ", .false.) , & + varatts( "WND ", "UAX ", "Mean wind, x-component ", "m s-1 ", " ", .false.) , & + varatts( "WND ", "UAY ", "Mean wind, y-component ", "m s-1 ", " ", .false.) , & + varatts( "AST ", "AS ", "Air-sea temperature difference ", "K ", " ", .false.) , & + varatts( "WLV ", "WLV ", "Water levels ", "m ", " ", .false.) , & + varatts( "ICE ", "ICE ", "Ice coverage ", "nd ", " ", .false.) , & + varatts( "IBG ", "BERG ", "Iceberg-induced damping ", "km-1 ", " ", .false.) , & + varatts( "TAUA ", "TAUAX ", "Atm momentum x ", "Pa ", " ", .false.) , & + varatts( "TAUA ", "TAUAY ", "Atm momentum y ", "Pa ", " ", .false.) , & + varatts( "RHO ", "RHOAIR ", "Air density ", "kg m-3 ", " ", .false.) , & + varatts( "IC1 ", "ICEH ", "Ice thickness ", "m ", " ", .false.) , & + varatts( "IC5 ", "ICEF ", "Ice floe diameter ", "m ", " ", .false.) & + ] + + ! 2 Standard mean wave Parameters + gridoutdefs(2,1:18) = [ & + varatts( "HS ", "HS ", "Significant wave height ", "m ", " ", .false.) , & + varatts( "LM ", "WLM ", "Mean wave length ", "m ", " ", .false.) , & + varatts( "T02 ", "T02 ", "Mean wave period (Tm0,2) ", "s ", " ", .false.) , & + varatts( "T0M1 ", "T0M1 ", "Mean wave period (Tm0,-1) ", "s ", " ", .false.) , & + varatts( "T01 ", "T01 ", "Mean wave period (Tm0,1) ", "s ", " ", .false.) , & + varatts( "FP ", "FP0 ", "Peak frequency ", "s-1 ", " ", .false.) , & + varatts( "DIR ", "THM ", "Mean wave direction ", "rad ", " ", .false.) , & + varatts( "SPR ", "THS ", "Mean directional spread ", "rad ", " ", .false.) , & + varatts( "DP ", "THP0 ", "Peak direction ", "rad ", " ", .false.) , & + varatts( "HIG ", "HSIG ", "Infragravity height ", "m ", " ", .false.) , & + varatts( "MXE ", "STMAXE ", "Max surface elev (STE) ", "m ", " ", .false.) , & + varatts( "MXES ", "STMAXD ", "St Dev Max surface elev (STE) ", "m ", " ", .false.) , & + varatts( "MXH ", "HMAXE ", "Max wave height (S.) ", "m ", " ", .false.) , & + varatts( "MXHC ", "HCMAXE ", "Max wave height from crest (STE) ", "m ", " ", .false.) , & + varatts( "SDMH ", "HMAXD ", "St Dev of MXC (STE) ", "m ", " ", .false.) , & + varatts( "SDMHC", "HCMAXD ", "St Dev of MXHC (STE) ", "m ", " ", .false.) , & + varatts( "WBT ", "WBT ", "Dominant wave breaking probability (b_T) ", "nd ", " ", .false.) , & + varatts( "WNM ", "WNMEAN ", "Mean wave number ", "m-1 ", " ", .false.) & + ] + + ! 3 Spectral Parameters + gridoutdefs(3,1:6) = [ & + varatts( "EF ", "EF ", "1D spectral density ", "m2 s ", "k ", .false.) , & + varatts( "TH1M ", "TH1M ", "Mean wave direction from a1,b2 ", "deg ", "k ", .false.) , & + 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.) & + ] + + ! 4 Spectral Partition Parameters + gridoutdefs(4,1:17) = [ & + varatts( "PHS ", "PHS ", "Partitioned wave heights ", "m ", "s ", .false.) , & + varatts( "PTP ", "PTP ", "Partitioned peak period ", "s ", "s ", .false.) , & + varatts( "PLP ", "PLP ", "Partitioned peak wave length ", "m ", "s ", .false.) , & + varatts( "PDIR ", "PDIR ", "Partitioned mean direction ", "deg ", "s ", .false.) , & + varatts( "PSPR ", "PSI ", "Partitioned mean directional spread ", "deg ", "s ", .false.) , & + varatts( "PWS ", "PWS ", "Partitioned wind sea fraction ", "nd ", "s ", .false.) , & + varatts( "PDP ", "PTHP0 ", "Peak wave direction of partition ", "deg ", "s ", .false.) , & + varatts( "PQP ", "PQP ", "Goda peakdedness parameter of partition ", "nd ", "s ", .false.) , & + varatts( "PPE ", "PPE ", "JONSWAP peak enhancement factor of partition ", "s-1 ", "s ", .false.) , & + varatts( "PGW ", "PGW ", "Gaussian frequency width of partition ", "nd ", "s ", .false.) , & + varatts( "PSW ", "PSW ", "Spectral width of partition ", "nd ", "s ", .false.) , & + varatts( "PTM10", "PTM1 ", "Mean wave period (m-1,0) of partition ", "s ", "s ", .false.) , & + varatts( "PT01 ", "PT1 ", "Mean wave period (m0,1) of partition ", "s ", "s ", .false.) , & + varatts( "PT02 ", "PT2 ", "Mean wave period (m0,2) of partition ", "s ", "s ", .false.) , & + varatts( "PEP ", "PEP ", "Peak spectral density of partition ", "m2 s rad-1", "s ", .false.) , & + varatts( "TWS ", "PWST ", "Total wind sea fraction ", "nd ", " ", .false.) , & + varatts( "PNR ", "PNR ", "Number of partitions ", "nd ", " ", .false.) & + ] + + ! 5 Atmosphere-waves layer + gridoutdefs(5,1:14) = [ & + 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( "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.) , & + varatts( "TWA ", "TAUWNX ", "Negative part of the wave-supported stress x ", "m2 s-2 ", " ", .false.) , & + varatts( "TWA ", "TAUWNY ", "Negative part of the wave-supported stress y ", "m2 s-2 ", " ", .false.) , & + varatts( "WCC ", "WCC ", "Whitecap coverage ", "nd ", " ", .false.) , & + varatts( "WCF ", "WCF ", "Whitecap foam thickness ", "m ", " ", .false.) , & + varatts( "WCH ", "WCH ", "Mean breaking wave heigh ", "m ", " ", .false.) , & + varatts( "WCM ", "WCM ", "Whitecap moment ", "nd ", " ", .false.) , & + varatts( "FWS ", "TWS ", "Wind sea mean period ", "s ", " ", .false.) & + ] + + ! 6 Wave-ocean layer + gridoutdefs(6,1:25) = [ & + varatts( "SXY ", "SXX ", "Radiation stresses xx ", "N m-1 ", " ", .false.) , & + varatts( "SXY ", "SYY ", "Radiation stresses yy ", "N m-1 ", " ", .false.) , & + varatts( "SXY ", "SXY ", "Radiation stresses xy ", "N m-1 ", " ", .false.) , & + varatts( "TWO ", "TAUOX ", "Wave to ocean momentum flux x ", "m2 s-2 ", " ", .false.) , & + varatts( "TWO ", "TAUOY ", "Wave to ocean momentum flux y ", "m2 s-2 ", " ", .false.) , & + varatts( "BHD ", "BHD ", "Bernoulli head (J term) ", "m2 s-2 ", " ", .false.) , & + varatts( "FOC ", "PHIOC ", "Wave to ocean energy flux ", "W m-2 ", " ", .false.) , & + varatts( "TUS ", "TUSX ", "Stokes transport x ", "m2 s-1 ", " ", .false.) , & + varatts( "TUS ", "TUSY ", "Stokes transport y ", "m2 s-1 ", " ", .false.) , & + varatts( "USS ", "USSX ", "Surface Stokes drift x ", "m s-1 ", " ", .false.) , & + varatts( "USS ", "USSY ", "Surface Stokes drift y ", "m s-1 ", " ", .false.) , & + varatts( "P2S ", "PRMS ", "Second-order sum pressure ", "m4 ", " ", .false.) , & + varatts( "P2S ", "TPMS ", "Second-order sum pressure ", "s-1 ", " ", .false.) , & + varatts( "USF ", "US3DX ", "Spectrum of surface Stokes drift x ", "m s-1 Hz-1", "k ", .false.) , & + varatts( "USF ", "US3DY ", "Spectrum of surface Stokes drift y ", "m s-1 Hz-1", "k ", .false.) , & + varatts( "P2L ", "P2SMS ", "Micro seism source term ", "Pa2 m2 s ", "m ", .false.) , & + varatts( "TWI ", "TAUICEX ", "Wave to sea ice stress x ", "m2 s-2 ", " ", .false.) , & + varatts( "TWI ", "TAUICEY ", "Wave to sea ice stress y ", "m2 s-2 ", " ", .false.) , & + varatts( "FIC ", "PHICE ", "Wave to sea ice energy flux ", "W m-2 ", " ", .false.) , & + varatts( "USP ", "USSPX ", "Partitioned surface Stokes drift x ", "m s-1 ", "p ", .false.) , & + varatts( "USP ", "USSPY ", "Partitioned surface Stokes drift y ", "m s-1 ", "p ", .false.) , & + varatts( "TWC ", "TAUOCX ", "Total wave to ocean stress x ", "Pa ", " ", .false.) , & + varatts( "TWC ", "TAUOCY ", "Total wave to ocean stress y ", "Pa ", " ", .false.) , & + varatts( "USSH ", "USSHX ", "Surface layer averaged Stokes drift x ", "m s-1 ", " ", .false.) , & + varatts( "USSH ", "USSHY ", "Surface layer averaged Stokes drift y ", "m s-1 ", " ", .false.) & + ] + + ! 7 Wave-bottom layer + gridoutdefs(7,1:10) = [ & + varatts( "ABR ", "ABAX ", "Near bottom rms wave excursion amplitudes x ", "m ", " ", .false.) , & + varatts( "ABR ", "ABAY ", "Near bottom rms wave excursion amplitudes y ", "m ", " ", .false.) , & + varatts( "UBR ", "UBAX ", "Near bottom rms wave velocities x ", "m s-1 ", " ", .false.) , & + varatts( "UBR ", "UBAY ", "Near bottom rms wave velocities y ", "m s-1 ", " ", .false.) , & + varatts( "BED ", "BED ", "Bottom roughness ", "m ", " ", .false.) , & + varatts( "BED ", "RIPPLEX ", "Sea bottom ripple wavelength x ", "m ", " ", .false.) , & + varatts( "BED ", "RIPPLEY ", "Sea bottom ripple wavelength y ", "m ", " ", .false.) , & + varatts( "FBB ", "PHIBBL ", "Energy flux due to bottom friction ", "W m-2 ", " ", .false.) , & + varatts( "TBB ", "TAUBBLX ", "Momentum flux due to bottom friction x ", "m2 s-2 ", " ", .false.) , & + varatts( "TBB ", "TAUBBLY ", "Momentum flux due to bottom friction y ", "m2 s-2 ", " ", .false.) & + ] + + ! 8 Spectrum parameters + gridoutdefs(8,1:9) = [ & + varatts( "MSS ", "MSSX ", "Surface mean square slope x ", "nd ", " ", .false.) , & + varatts( "MSS ", "MSSY ", "Surface mean square slope y ", "nd ", " ", .false.) , & + varatts( "MSC ", "MSCX ", "Spectral level at high frequency tail x ", "nd ", " ", .false.) , & + varatts( "MSC ", "MSCY ", "Spectral level at high frequency tail y ", "nd ", " ", .false.) , & + varatts( "WL02 ", "WL02X ", "East/X North/Y mean wavelength component ", "nd ", " ", .false.) , & + varatts( "WL02 ", "WL02Y ", "East/X North/Y mean wavelength component ", "nd ", " ", .false.) , & + varatts( "AXT ", "ALPXT ", "Correl sea surface gradients (x,t) ", "nd ", " ", .false.) , & + varatts( "AYT ", "ALPYT ", "Correl sea surface gradients (y,t) ", "nd ", " ", .false.) , & + varatts( "AXY ", "ALPXY ", "Correl sea surface gradients (x,y) ", "nd ", " ", .false.) & + ] + + ! 9 Numerical diagnostics + gridoutdefs(9,1:5) = [ & + varatts( "DTD ", "DTDYN ", "Average time step in integration ", "min ", " ", .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.) , & + varatts( "CFK ", "CFLKMAX ", "Max. CFL number for k-advection ", "nd ", " ", .false.) & + ] + + ! 10 User defined + gridoutdefs(10,1:2) = [ & + varatts( "U1 ", "U1 ", "User defined 1 ", "nd ", " ", .false.) , & + varatts( "U2 ", "U2 ", "User defined 2 ", "nd ", " ", .false.) & + ] + end subroutine define_fields +end module wav_history_mod diff --git a/model/src/wav_pio_mod.F90 b/model/src/wav_pio_mod.F90 new file mode 100644 index 000000000..4f9b637a4 --- /dev/null +++ b/model/src/wav_pio_mod.F90 @@ -0,0 +1,420 @@ +!> @file wav_pio +!! +!> @brief Manage PIO for WW3 +!! +!> @author Denise.Worthen@noaa.gov +!> @date 08-02-2024 +module wav_pio_mod + + use w3gdatmd , only : nk, nx, ny, mapsf + use w3parall , only : init_get_isea + use w3gdatmd , only : nseal + use pio + use netcdf +#ifdef W3_PDLIB + use yowNodepool , only : ng +#endif + implicit none + + private + + interface wav_pio_initdecomp + module procedure wav_pio_initdecomp_2d + module procedure wav_pio_initdecomp_3d + end interface wav_pio_initdecomp + + integer :: pio_iotype + integer :: pio_ioformat + type(iosystem_desc_t), pointer :: wav_pio_subsystem + + public :: wav_pio_init + public :: pio_iotype + public :: pio_ioformat + public :: wav_pio_subsystem + public :: wav_pio_initdecomp + public :: handle_err + + !=============================================================================== +contains + !=============================================================================== + !> Configure PIO for WW3 + !! + !> @details Use either CESM shr code or configuration variables to configure PIO. + !! This configuration code is lifted from CMEPS. + !! + !! @param gcomp an ESMF_GridComp object + !! @param mpi_comm the MPI communicator + !! @param[in] stdout the logfile unit on the root_task + !! @param[in] numprocs naproc/nthrds + !! @param[out] rc a return code + !! + !> @author Denise.Worthen@noaa.gov + !> @date 08-02-2024 + subroutine wav_pio_init(gcomp, mpi_comm, stdout, numprocs, rc) + +#ifdef CESMCOUPLED + use shr_pio_mod, only : shr_pio_getiosys, shr_pio_getiotype, shr_pio_getioformat +#endif + use ESMF , only : ESMF_GridComp, ESMF_UtilStringUpperCase, ESMF_VM, ESMF_FAILURE + use ESMF , only : ESMF_SUCCESS, ESMF_LogWrite, ESMF_LOGMSG_ERROR + use NUOPC , only : NUOPC_CompAttributeGet + use wav_kind_mod , only : CL=>SHR_KIND_CL, CS=>SHR_KIND_CS + use w3odatmd , only : iaproc + use wav_shr_mod , only : chkerr + + ! input/output arguments + type(ESMF_GridComp), intent(in) :: gcomp + integer , intent(in) :: mpi_comm + integer , intent(in) :: stdout + integer , intent(in) :: numprocs + integer , intent(out) :: rc + + integer :: pio_numiotasks + integer :: pio_stride + integer :: pio_rearranger + integer :: pio_root + integer :: pio_debug_level + character(len=CS) :: cvalue + logical :: isPresent, isSet + integer :: my_task, master_task + character(len=CS) :: subname='wav_pio_init' + character(*), parameter :: u_FILE_u = & !< a character string for an ESMF log message + __FILE__ + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + +#ifdef CESMCOUPLED + ! TODO: needs testing + wav_pio_subsystem => shr_pio_getiosys(inst_name) + pio_iotype = shr_pio_getiotype(inst_name) + if ((pio_iotype==PIO_IOTYPE_NETCDF).or.(pio_iotype==PIO_IOTYPE_PNETCDF)) then + nmode0 = shr_pio_getioformat(inst_name) + else + nmode0 = 0 + endif + + call pio_seterrorhandling(wav_pio_subsystem, PIO_RETURN_ERROR) +#else + my_task = iaproc - 1 + master_task = 0 + + ! code lifted from CMEPS med_io_mod.F90 + ! query component specific PIO attributes + ! pio_netcdf_format + call NUOPC_CompAttributeGet(gcomp, name='pio_netcdf_format', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + if (isPresent .and. isSet) then + cvalue = ESMF_UtilStringUpperCase(cvalue) + if (trim(cvalue) .eq. 'CLASSIC') then + pio_ioformat = 0 + else if (trim(cvalue) .eq. '64BIT_OFFSET') then + pio_ioformat = PIO_64BIT_OFFSET + else if (trim(cvalue) .eq. '64BIT_DATA') then + pio_ioformat = PIO_64BIT_DATA + else + call ESMF_LogWrite(trim(subname)//': need to provide valid option for pio_ioformat ' & + //'(CLASSIC|64BIT_OFFSET|64BIT_DATA)', ESMF_LOGMSG_ERROR) + rc = ESMF_FAILURE + return + end if + else + cvalue = '64BIT_OFFSET' + pio_ioformat = PIO_64BIT_OFFSET + end if + if (my_task == 0) write(stdout,*) trim(subname), ' : pio_netcdf_format = ', trim(cvalue), pio_ioformat + + ! pio_typename + call NUOPC_CompAttributeGet(gcomp, name='pio_typename', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + if (isPresent .and. isSet) then + cvalue = ESMF_UtilStringUpperCase(cvalue) + if (trim(cvalue) .eq. 'NETCDF') then + pio_iotype = PIO_IOTYPE_NETCDF + else if (trim(cvalue) .eq. 'PNETCDF') then + pio_iotype = PIO_IOTYPE_PNETCDF + else if (trim(cvalue) .eq. 'NETCDF4C') then + pio_iotype = PIO_IOTYPE_NETCDF4C + else if (trim(cvalue) .eq. 'NETCDF4P') then + pio_iotype = PIO_IOTYPE_NETCDF4P + else + call ESMF_LogWrite(trim(subname)//': need to provide valid option for pio_typename ' & + //'(NETCDF|PNETCDF|NETCDF4C|NETCDF4P)', ESMF_LOGMSG_ERROR) + rc = ESMF_FAILURE + return + end if + else + cvalue = 'NETCDF' + pio_iotype = PIO_IOTYPE_NETCDF + end if + if (my_task == 0) write(stdout,*) trim(subname), ' : pio_typename = ', trim(cvalue), pio_iotype + + ! pio_root + call NUOPC_CompAttributeGet(gcomp, name='pio_root', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + if (isPresent .and. isSet) then + read(cvalue,*) pio_root + if (pio_root < 0) then + pio_root = 1 + endif + pio_root = min(pio_root, numprocs-1) + else + pio_root = 1 + end if + if (my_task == 0) write(stdout,*) trim(subname), ' : pio_root = ', pio_root + + ! pio_stride + call NUOPC_CompAttributeGet(gcomp, name='pio_stride', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + if (isPresent .and. isSet) then + read(cvalue,*) pio_stride + else + pio_stride = -99 + end if + if (my_task == 0) write(stdout,*) trim(subname), ' : pio_stride = ', pio_stride + + ! pio_numiotasks + call NUOPC_CompAttributeGet(gcomp, name='pio_numiotasks', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + if (isPresent .and. isSet) then + read(cvalue,*) pio_numiotasks + else + pio_numiotasks = -99 + end if + if (my_task == 0) write(stdout,*) trim(subname), ' : pio_numiotasks = ', pio_numiotasks + + ! check for parallel IO, it requires at least two io pes + if (numprocs > 1 .and. pio_numiotasks == 1 .and. & + (pio_iotype .eq. PIO_IOTYPE_PNETCDF .or. pio_iotype .eq. PIO_IOTYPE_NETCDF4P)) then + pio_numiotasks = 2 + pio_stride = min(pio_stride, numprocs/2) + if (my_task == 0) then + write(stdout,*) ' parallel io requires at least two io pes - following parameters are updated:' + write(stdout,*) trim(subname), ' : pio_stride = ', pio_stride + write(stdout,*) trim(subname), ' : pio_numiotasks = ', pio_numiotasks + end if + endif + + ! check/set/correct io pio parameters + if (pio_stride > 0 .and. pio_numiotasks < 0) then + pio_numiotasks = max(1, numprocs/pio_stride) + if (my_task == 0) write(stdout,*) trim(subname), ' : update pio_numiotasks = ', pio_numiotasks + else if(pio_numiotasks > 0 .and. pio_stride < 0) then + pio_stride = max(1, numprocs/pio_numiotasks) + if (my_task == 0) write(stdout,*) trim(subname), ' : update pio_stride = ', pio_stride + else if(pio_numiotasks < 0 .and. pio_stride < 0) then + pio_stride = max(1,numprocs/4) + pio_numiotasks = max(1,numprocs/pio_stride) + if (my_task == 0) write(stdout,*) trim(subname), ' : update pio_numiotasks = ', pio_numiotasks + if (my_task == 0) write(stdout,*) trim(subname), ' : update pio_stride = ', pio_stride + end if + if (pio_stride == 1) then + pio_root = 0 + endif + + if (pio_root + (pio_stride)*(pio_numiotasks-1) >= numprocs .or. & + pio_stride <= 0 .or. pio_numiotasks <= 0 .or. pio_root < 0 .or. pio_root > numprocs-1) then + if (numprocs < 100) then + pio_stride = max(1, numprocs/4) + else if(numprocs < 1000) then + pio_stride = max(1, numprocs/8) + else + pio_stride = max(1, numprocs/16) + end if + if(pio_stride > 1) then + pio_numiotasks = numprocs/pio_stride + pio_root = min(1, numprocs-1) + else + pio_numiotasks = numprocs + pio_root = 0 + end if + if (my_task == 0) then + write(stdout,*) 'pio_stride, iotasks or root out of bounds - resetting to defaults:' + write(stdout,*) trim(subname), ' : pio_root = ', pio_root + write(stdout,*) trim(subname), ' : pio_stride = ', pio_stride + write(stdout,*) trim(subname), ' : pio_numiotasks = ', pio_numiotasks + end if + end if + + ! pio_rearranger + call NUOPC_CompAttributeGet(gcomp, name='pio_rearranger', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + if (isPresent .and. isSet) then + cvalue = ESMF_UtilStringUpperCase(cvalue) + if (trim(cvalue) .eq. 'BOX') then + pio_rearranger = PIO_REARR_BOX + else if (trim(cvalue) .eq. 'SUBSET') then + pio_rearranger = PIO_REARR_SUBSET + else + call ESMF_LogWrite(trim(subname)//': need to provide valid option for pio_rearranger (BOX|SUBSET)', ESMF_LOGMSG_ERROR) + rc = ESMF_FAILURE + return + end if + else + cvalue = 'SUBSET' + pio_rearranger = PIO_REARR_SUBSET + end if + if (my_task == 0) write(stdout,*) trim(subname), ' : pio_rearranger = ', trim(cvalue), pio_rearranger + + ! init PIO + if (my_task == 0) then + write(stdout,*) trim(subname),' calling pio init' + write(stdout,*) trim(subname), ' : pio_root = ', pio_root + write(stdout,*) trim(subname), ' : pio_stride = ', pio_stride + write(stdout,*) trim(subname), ' : pio_numiotasks = ', pio_numiotasks + end if + + allocate(wav_pio_subsystem) + call pio_init(my_task, mpi_comm, pio_numiotasks, master_task, pio_stride, pio_rearranger, & + wav_pio_subsystem, base=pio_root) + + ! PIO debug related options + ! pio_debug_level + call NUOPC_CompAttributeGet(gcomp, name='pio_debug_level', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + read(cvalue,*) pio_debug_level + if (pio_debug_level < 0 .or. pio_debug_level > 6) then + call ESMF_LogWrite(trim(subname)//': need to provide valid option for pio_debug_level (0-6)', ESMF_LOGMSG_ERROR) + rc = ESMF_FAILURE + return + end if + else + pio_debug_level = 0 + end if + if (my_task == 0) write(stdout,*) trim(subname), ' : pio_debug_level = ', pio_debug_level + + ! set PIO debug level + call pio_setdebuglevel(pio_debug_level) + + call pio_seterrorhandling(wav_pio_subsystem, PIO_RETURN_ERROR) +#endif + end subroutine wav_pio_init + + !=============================================================================== + !> Define a decomposition for a 2d variable in WW3 + !! + !! @param[out] iodesc the PIO decomposition handle + !! @param[out] use_int define a decomposition for an integer array + !! + !> @author Denise.Worthen@noaa.gov + !> @date 08-02-2024 + subroutine wav_pio_initdecomp_2d(iodesc, use_int) + + type(io_desc_t), intent(out) :: iodesc + logical , optional, intent(in) :: use_int + + ! local variables + integer :: n, isea, jsea, ix, iy, nseal_cpl + logical :: luse_int + integer(kind=PIO_OFFSET_KIND) :: lnx,lny + integer(kind=PIO_OFFSET_KIND), allocatable :: dof2d(:) +#ifdef W3_PDLIB + nseal_cpl = nseal - ng +#else + nseal_cpl = nseal +#endif + luse_int = .false. + if (present(use_int)) luse_int = use_int + + allocate(dof2d(nseal_cpl)) + dof2d = 0 + lnx = int(nx,PIO_OFFSET_KIND) + lny = int(ny,PIO_OFFSET_KIND) + + n = 0 + do jsea = 1,nseal_cpl + call init_get_isea(isea, jsea) + ix = mapsf(isea,1) ! global ix + iy = mapsf(isea,2) ! global iy + n = n+1 + dof2d(n) = (iy-1)*lnx + ix ! local index : global index + end do + + if (luse_int) then + call pio_initdecomp(wav_pio_subsystem, PIO_INT, (/nx,ny/), dof2d, iodesc) + else + call pio_initdecomp(wav_pio_subsystem, PIO_REAL, (/nx,ny/), dof2d, iodesc) + end if + deallocate(dof2d) + + end subroutine wav_pio_initdecomp_2d + + !=============================================================================== + !> Define a decomposition for a 3d variable in WW3 + !! + !! @param[in] nz the non-spatial dimension + !! @param[out] iodesc the PIO decomposition handle + !! + !> @author Denise.Worthen@noaa.gov + !> @date 08-02-2024 + subroutine wav_pio_initdecomp_3d(nz, iodesc) + + integer , intent(in) :: nz + type(io_desc_t) , intent(out) :: iodesc + + ! local variables + integer :: n, k, isea, jsea, ix, iy, nseal_cpl + integer(kind=PIO_OFFSET_KIND) :: lnx,lny + integer(kind=PIO_OFFSET_KIND), allocatable :: dof3d(:) +#ifdef W3_PDLIB + nseal_cpl = nseal - ng +#else + nseal_cpl = nseal +#endif + allocate(dof3d(nz*nseal_cpl)) + + dof3d = 0 + lnx = int(nx,PIO_OFFSET_KIND) + lny = int(ny,PIO_OFFSET_KIND) + + n = 0 + do k = 1,nz + do jsea = 1,nseal_cpl + call init_get_isea(isea, jsea) + ix = mapsf(isea,1) ! global ix + iy = mapsf(isea,2) ! global iy + n = n+1 + dof3d(n) = ((iy-1)*lnx + ix) + (k-1)*lnx*lny ! local index : global index + end do + end do + + call pio_initdecomp(wav_pio_subsystem, PIO_REAL, (/nx,ny,nz/), dof3d, iodesc) + deallocate(dof3d) + + end subroutine wav_pio_initdecomp_3d + + !=============================================================================== + !> Handle errors + !! + !! @param[in] ierr the error code + !! @param[in] string the error message + !! + !> @author Denise.Worthen@noaa.gov + !> @date 08-02-2024 + subroutine handle_err(ierr,string) + + use w3odatmd , only : ndse + use w3servmd , only : extcde + + ! input/output variables + integer , intent(in) :: ierr + character(len=*), intent(in) :: string + + integer :: strerror_status + character(len=pio_max_name) :: err_msg + + if (ierr /= PIO_NOERR) then + strerror_status = pio_strerror(ierr, err_msg) + write(ndse,*) "*** WAVEWATCH III netcdf error: ",trim(string),':',trim(err_msg) + call extcde ( 49 ) + end if + end subroutine handle_err + +end module wav_pio_mod diff --git a/model/src/wav_restart_mod.F90 b/model/src/wav_restart_mod.F90 new file mode 100644 index 000000000..9d1c67780 --- /dev/null +++ b/model/src/wav_restart_mod.F90 @@ -0,0 +1,538 @@ +!> @file wav_restart_mod +!! +!> @brief Handle WW3 restart files as netCDF using PIO +!! +!> @author Denise.Worthen@noaa.gov +!> @date 08-26-2024 +module wav_restart_mod + + use w3parall , only : init_get_isea + use w3adatmd , only : nsealm + use w3gdatmd , only : nth, nk, nx, ny, mapsf, nspec, nseal, nsea + use w3odatmd , only : ndso, iaproc, addrstflds, rstfldlist, rstfldcnt + use w3wdatmd , only : ice + use wav_pio_mod , only : pio_iotype, pio_ioformat, wav_pio_subsystem + use wav_pio_mod , only : handle_err, wav_pio_initdecomp +#ifdef W3_PDLIB + use yowNodepool , only : ng +#endif + use pio + use netcdf + + implicit none + + private + + type(file_desc_t) :: pioid + type(var_desc_t) :: varid + type(io_desc_t) :: iodesc2dint + type(io_desc_t) :: iodesc2d + + integer(kind=Pio_Offset_Kind) :: frame + + public :: write_restart + public :: read_restart + + ! used/reused in module + character(len=4) :: cspec + character(len=12) :: vname + integer :: ik, ith, ix, iy, kk, isea, jsea, ierr, i + + !=============================================================================== +contains + !=============================================================================== + !> Write a WW3 restart file + !! + !! @details Called by w3wavemd to write a restart file at a given frequency or + !! time + !! + !! @param[in] fname the time-stamped file name + !! @param[in] va the va array + !! @param[in] mapsta the mapsta + 8*mapst2 array + !! + !> author DeniseWorthen@noaa.gov + !> @date 08-26-2024 + subroutine write_restart (fname, va, mapsta) + + use w3odatmd , only : time_origin, calendar_name, elapsed_secs + + real , intent(in) :: va(1:nspec,0:nsealm) + integer , intent(in) :: mapsta(ny,nx) + character(len=*), intent(in) :: fname + + ! local variables + integer :: timid, xtid, ytid + integer :: nseal_cpl, nmode + integer :: dimid(3) + real , allocatable :: lva(:,:) + integer, allocatable :: lmap(:) + !------------------------------------------------------------------------------- + +#ifdef W3_PDLIB + nseal_cpl = nseal - ng +#else + nseal_cpl = nseal +#endif + allocate(lva(1:nseal_cpl,1:nspec)) + allocate(lmap(1:nseal_cpl)) + lva(:,:) = 0.0 + lmap(:) = 0 + + ! create the netcdf file + frame = 1 + pioid%fh = -1 + nmode = pio_clobber + ! only applies to classic NETCDF files. + if (pio_iotype == PIO_IOTYPE_NETCDF .or. pio_iotype == PIO_IOTYPE_PNETCDF) then + nmode = ior(nmode,pio_ioformat) + endif + ierr = pio_createfile(wav_pio_subsystem, pioid, pio_iotype, trim(fname), nmode) + call handle_err(ierr, 'pio_create') + if (iaproc == 1) write(ndso,'(a)')' Writing restart file '//trim(fname) + + ierr = pio_def_dim(pioid, 'nx', nx, xtid) + ierr = pio_def_dim(pioid, 'ny', ny, ytid) + ierr = pio_def_dim(pioid, 'time', PIO_UNLIMITED, timid) + + ! define the time variable + ierr = pio_def_var(pioid, 'time', PIO_DOUBLE, (/timid/), varid) + call handle_err(ierr,'def_timevar') + ierr = pio_put_att(pioid, varid, 'units', trim(time_origin)) + call handle_err(ierr,'def_time_units') + ierr = pio_put_att(pioid, varid, 'calendar', trim(calendar_name)) + call handle_err(ierr,'def_time_calendar') + + ! define the nth,nk sizes + ierr = pio_def_var(pioid, 'nth', PIO_INT, varid) + call handle_err(ierr,'def_nth') + ierr = pio_put_att(pioid, varid, 'long_name', 'number of direction bins') + ierr = pio_def_var(pioid, 'nk', PIO_INT, varid) + call handle_err(ierr,'def_nk') + ierr = pio_put_att(pioid, varid, 'long_name', 'number of frequencies') + + ! write each nspec as separate variable + do kk = 1,nspec + write(cspec,'(i4.4)')kk + vname = 'va'//cspec + dimid = (/xtid, ytid, timid/) + ierr = pio_def_var(pioid, trim(vname), PIO_REAL, dimid, varid) + call handle_err(ierr, 'define variable '//trim(vname)) + ierr = pio_put_att(pioid, varid, '_FillValue', nf90_fill_float) + call handle_err(ierr, 'define _FillValue '//trim(vname)) + end do + + vname = 'mapsta' + ierr = pio_def_var(pioid, trim(vname), PIO_INT, (/xtid, ytid, timid/), varid) + call handle_err(ierr, 'define variable '//trim(vname)) + ierr = pio_put_att(pioid, varid, '_FillValue', nf90_fill_int) + call handle_err(ierr, 'define _FillValue '//trim(vname)) + + ! define any requested additional fields + if (addrstflds) then + do i = 1,rstfldcnt + vname = trim(rstfldlist(i)) + ierr = pio_def_var(pioid, trim(vname), PIO_REAL, (/xtid, ytid, timid/), varid) + call handle_err(ierr, 'define variable '//trim(vname)) + ierr = pio_put_att(pioid, varid, '_FillValue', nf90_fill_float) + call handle_err(ierr, 'define _FillValue '//trim(vname)) + end do + end if + ! end variable definitions + ierr = pio_enddef(pioid) + call handle_err(ierr, 'end variable definition') + + ! write the freq and direction sizes + ierr = pio_inq_varid(pioid, 'nth', varid) + call handle_err(ierr, 'inquire variable nth ') + ierr = pio_put_var(pioid, varid, nth) + call handle_err(ierr, 'put nth') + ierr = pio_inq_varid(pioid, 'nk', varid) + call handle_err(ierr, 'inquire variable nk ') + ierr = pio_put_var(pioid, varid, nk) + call handle_err(ierr, 'put nk') + + ! initialize the decomp + call wav_pio_initdecomp(iodesc2dint, use_int=.true.) + call wav_pio_initdecomp(iodesc2d) + + ! write the time + ierr = pio_inq_varid(pioid, 'time', varid) + call handle_err(ierr, 'inquire variable time ') + ierr = pio_put_var(pioid, varid, (/1/), real(elapsed_secs,8)) + call handle_err(ierr, 'put time') + + ! mapsta is global + do jsea = 1,nseal_cpl + call init_get_isea(isea, jsea) + ix = mapsf(isea,1) + iy = mapsf(isea,2) + lmap(jsea) = mapsta(iy,ix) + end do + + ! write PE local map + vname = 'mapsta' + ierr = pio_inq_varid(pioid, trim(vname), varid) + call handle_err(ierr, 'inquire variable '//trim(vname)) + call pio_setframe(pioid, varid, int(1,kind=Pio_Offset_Kind)) + call pio_write_darray(pioid, varid, iodesc2dint, lmap, ierr) + call handle_err(ierr, 'put variable '//trim(vname)) + + ! write va + do jsea = 1,nseal_cpl + kk = 0 + do ik = 1,nk + do ith = 1,nth + kk = kk + 1 + lva(jsea,kk) = va(kk,jsea) + end do + end do + end do + + do kk = 1,nspec + write(cspec,'(i4.4)')kk + vname = 'va'//cspec + ierr = pio_inq_varid(pioid, trim(vname), varid) + call handle_err(ierr, 'inquire variable '//trim(vname)) + call pio_setframe(pioid, varid, int(1,kind=PIO_OFFSET_KIND)) + call pio_write_darray(pioid, varid, iodesc2d, lva(:,kk), ierr) + call handle_err(ierr, 'put variable '//trim(vname)) + end do + + ! write requested additional global(nsea) fields + if (addrstflds) then + do i = 1,rstfldcnt + vname = trim(rstfldlist(i)) + if (vname == 'ice')call write_globalfield(vname, nseal_cpl, ice(1:nsea)) + end do + end if + + call pio_syncfile(pioid) + call pio_freedecomp(pioid, iodesc2d) + call pio_freedecomp(pioid, iodesc2dint) + call pio_closefile(pioid) + + end subroutine write_restart + + !=============================================================================== + !> Read a WW3 restart file + !! + !> @details Called by w3init to read a restart file which is known to exist or to + !! initialize a set of variables when the filename is "none". + !! + !! @param[in] fname the time-stamped file name + !! @param[out] va the va array, optional + !! @param[out] mapsta the mapsta array, optional + !! @param[inout] mapst2 the mapst2 array, optional + !! + !> author DeniseWorthen@noaa.gov + !> @date 08-26-2024 + subroutine read_restart (fname, va, mapsta, mapst2) + + use mpi_f08 + use w3adatmd , only : mpi_comm_wave + use w3gdatmd , only : sig + use w3idatmd , only : icei + use w3wdatmd , only : time, tlev, tice, trho, tic1, tic5, wlv, asf, fpis + + character(len=*) , intent(in) :: fname + real , optional , intent(out) :: va(1:nspec,0:nsealm) + integer, optional , intent(out) :: mapsta(ny,nx) + integer, optional , intent(inout) :: mapst2(ny,nx) + + ! local variables + type(MPI_Comm) :: wave_communicator ! needed for mpi_f08 + integer, allocatable :: global_input(:), global_output(:) + integer :: nseal_cpl + integer :: ifill + real :: rfill + real , allocatable :: lva(:,:) + integer, allocatable :: lmap(:) + integer, allocatable :: lmap2d(:,:) + integer, allocatable :: st2init(:,:) + !------------------------------------------------------------------------------- + + ! cold start, set initial values and return. + if (trim(fname) == 'none') then + tlev(1) = -1 + tlev(2) = 0 + tice(1) = -1 + tice(2) = 0 + trho(1) = -1 + trho(2) = 0 + tic1(1) = -1 + tic1(2) = 0 + tic5(1) = -1 + tic5(2) = 0 + wlv = 0. + ice = 0. + asf = 1. + fpis = sig(nk) + if (iaproc == 1) write(ndso,'(a)')' Initializing WW3 at rest ' + return + end if + + ! read a netcdf restart + wave_communicator%mpi_val = MPI_COMM_WAVE +#ifdef W3_PDLIB + nseal_cpl = nseal - ng +#else + nseal_cpl = nseal +#endif + allocate(lva(1:nseal_cpl,1:nspec)) + allocate(lmap(1:nseal_cpl)) + allocate(lmap2d(1:ny,1:nx)) + allocate(st2init(1:ny,1:nx)) + lva(:,:) = 0.0 + lmap(:) = 0 + lmap2d(:,:) = 0 + + ! save a copy of initial mapst2 from mod_def + st2init = mapst2 + + ! all times are restart times + tlev = time + tice = time + trho = time + tic1 = time + tic5 = time + frame = 1 + ierr = pio_openfile(wav_pio_subsystem, pioid, pio_iotype, trim(fname), pio_nowrite) + call handle_err(ierr, 'open file '//trim(fname)) + if (iaproc == 1) write(ndso,'(a)')' Reading restart file '//trim(fname) + + ! check the field dimensions and sizes against the current values + call checkfile() + + ! initialize the decomp + call wav_pio_initdecomp(iodesc2dint, use_int=.true.) + call wav_pio_initdecomp(iodesc2d) + + do kk = 1,nspec + write(cspec,'(i4.4)')kk + vname = 'va'//cspec + ierr = pio_inq_varid(pioid, trim(vname), varid) + call handle_err(ierr, 'inquire variable '//trim(vname)) + call pio_setframe(pioid, varid, frame) + ierr = pio_get_att(pioid, varid, "_FillValue", rfill) + call handle_err(ierr, 'get variable _FillValue'//trim(vname)) + call pio_read_darray(pioid, varid, iodesc2d, lva(:,kk), ierr) + call handle_err(ierr, 'get variable '//trim(vname)) + end do + + va = 0.0 + do jsea = 1,nseal_cpl + kk = 0 + do ik = 1,nk + do ith = 1,nth + kk = kk + 1 + if (lva(jsea,kk) .ne. rfill) then + va(kk,jsea) = lva(jsea,kk) + end if + end do + end do + end do + + vname = 'mapsta' + ierr = pio_inq_varid(pioid, trim(vname), varid) + call handle_err(ierr, 'inquire variable '//trim(vname)) + call pio_setframe(pioid, varid, frame) + call pio_read_darray(pioid, varid, iodesc2dint, lmap, ierr) + call handle_err(ierr, 'get variable '//trim(vname)) + ierr = pio_get_att(pioid, varid, "_FillValue", ifill) + call handle_err(ierr, 'get variable _FillValue'//trim(vname)) + + ! fill global array with PE local values + allocate(global_input(nsea)) + allocate(global_output(nsea)) + global_input = 0 + global_output = 0 + do jsea = 1,nseal_cpl + call init_get_isea(isea, jsea) + if (lmap(jsea) .ne. ifill) then + global_input(isea) = lmap(jsea) + end if + end do + ! reduce across all PEs to create global array + call MPI_AllReduce(global_input, global_output, nsea, MPI_INTEGER, MPI_SUM, wave_communicator, ierr) + + ! fill global array on each PE + do isea = 1,nsea + ix = mapsf(isea,1) + iy = mapsf(isea,2) + lmap2d(iy,ix) = global_output(isea) + end do + deallocate(global_input) + deallocate(global_output) + + mapsta = mod(lmap2d+2,8) - 2 + mapst2 = st2init + (lmap2d-mapsta)/8 + + ! read additional global(nsea) restart fields + if (addrstflds) then + do i = 1,rstfldcnt + vname = trim(rstfldlist(i)) + if (vname == 'ice')call read_globalfield(wave_communicator, vname, nseal_cpl, ice(1:nsea), icei) + end do + end if + + call pio_syncfile(pioid) + call pio_freedecomp(pioid, iodesc2d) + call pio_freedecomp(pioid, iodesc2dint) + call pio_closefile(pioid) + + end subroutine read_restart + + !=============================================================================== + !> Write a decomposed array of (nsea) global values + !! + !! @param[in] vname the variable name + !! @param[in] nseal_cpl the PE local dimension, disregarding halos + !! @param[in] global_input the global array + !! + !> author DeniseWorthen@noaa.gov + !> @date 09-22-2024 + subroutine write_globalfield(vname, nseal_cpl, global_input) + + character(len=*) , intent(in) :: vname + integer , intent(in) :: nseal_cpl + real , intent(in) :: global_input(:) + + ! local variable + real, allocatable :: lvar(:) + + allocate(lvar(1:nseal_cpl)) + + lvar(:) = 0.0 + do jsea = 1,nseal_cpl + call init_get_isea(isea, jsea) + lvar(jsea) = global_input(isea) + end do + + !write PE local field + ierr = pio_inq_varid(pioid, trim(vname), varid) + call handle_err(ierr, 'inquire variable '//trim(vname)) + call pio_setframe(pioid, varid, int(1,kind=Pio_Offset_Kind)) + call pio_write_darray(pioid, varid, iodesc2d, lvar, ierr) + call handle_err(ierr, 'put variable '//trim(vname)) + + end subroutine write_globalfield + + !=============================================================================== + !> Read a decomposed array of (nsea) global values and return a global field on + !! each DE + !! + !! @param[in] wave_communicator the MPI handle + !! @param[in] vname the variable name + !! @param[in] nseal_cpl the PE local dimension, disregarding halos + !! @param[out] global_output the global array, nsea points on each DE + !! @param[out] global_2d the global array, (nx,ny) points on each DE + !! + !> author DeniseWorthen@noaa.gov + !> @date 09-22-2024 + subroutine read_globalfield(wave_communicator, vname, nseal_cpl, global_output, global_2d) + + use mpi_f08 + + type(MPI_Comm) , intent(in) :: wave_communicator ! needed for mpi_f08 + character(len=*) , intent(in) :: vname + integer , intent(in) :: nseal_cpl + real , intent(out) :: global_output(:) + real , intent(out) :: global_2d(:,:) + + ! local variables + real, allocatable :: global_input(:) + real :: rfill + real, allocatable :: lvar(:) + + allocate(lvar(1:nseal_cpl)) + lvar(:) = 0.0 + + ierr = pio_inq_varid(pioid, trim(vname), varid) + call handle_err(ierr, 'inquire variable '//trim(vname)) + call pio_setframe(pioid, varid, frame) + call pio_read_darray(pioid, varid, iodesc2d, lvar, ierr) + call handle_err(ierr, 'get variable '//trim(vname)) + ierr = pio_get_att(pioid, varid, "_FillValue", rfill) + call handle_err(ierr, 'get variable _FillValue'//trim(vname)) + + ! fill global array with PE local values + allocate(global_input(nsea)) + global_input = 0.0 + global_output = 0.0 + do jsea = 1,nseal_cpl + call init_get_isea(isea, jsea) + if (lvar(jsea) .ne. rfill) then + global_input(isea) = lvar(jsea) + end if + end do + ! reduce across all PEs to create global array + call MPI_AllReduce(global_input, global_output, nsea, MPI_REAL, MPI_SUM, wave_communicator, ierr) + deallocate(global_input) + + global_2d = 0.0 + do isea = 1,nsea + ix = mapsf(isea,1) + iy = mapsf(isea,2) + global_2d(ix,iy) = global_output(isea) + end do + + end subroutine read_globalfield + + !=============================================================================== + !> Check that a restart file has the expected dimensions and sizes + !! + !> author DeniseWorthen@noaa.gov + !> @date 10-15-2024 + subroutine checkfile() + + use w3odatmd , only : ndse + use w3servmd , only : extcde + + integer :: dimid, ivar + integer(kind=PIO_OFFSET_KIND) :: dimlen + + ! check dimension nx + vname = 'nx' + ierr = pio_inq_dimid(pioid, vname, dimid) + call handle_err(ierr, 'inquire dimension '//trim(vname)) + ierr = pio_inq_dimlen(pioid, dimid, dimlen) + if (dimlen /= int(nx,PIO_OFFSET_KIND)) then + write(ndse,*) '*** WAVEWATCH III restart error: '//trim(vname)//' does not match expected value' + call extcde ( 49 ) + end if + + ! check dimension ny + vname = 'ny' + ierr = pio_inq_dimid(pioid, vname, dimid) + call handle_err(ierr, 'inquire dimension '//trim(vname)) + ierr = pio_inq_dimlen(pioid, dimid, dimlen) + if (dimlen /= int(ny,PIO_OFFSET_KIND)) then + write(ndse,*) '*** WAVEWATCH III restart error: '//trim(vname)//' does not match expected value' + call extcde ( 49 ) + end if + + ! check number of directions + vname = 'nth' + ierr = pio_inq_varid(pioid, vname, varid) + call handle_err(ierr, 'inquire variable '//trim(vname)) + ierr = pio_get_var(pioid, varid, ivar) + call handle_err(ierr, 'get variable '//trim(vname)) + if (ivar .ne. nth) then + write(ndse,*) '*** WAVEWATCH III restart error: '//trim(vname)//' does not match expected value' + call extcde ( 49 ) + end if + + ! check number of frequencies + vname = 'nk' + ierr = pio_inq_varid(pioid, vname, varid) + call handle_err(ierr, 'inquire variable '//trim(vname)) + ierr = pio_get_var(pioid, varid, ivar) + call handle_err(ierr, 'get variable '//trim(vname)) + if (ivar .ne. nk) then + write(ndse,*) '*** WAVEWATCH III restart error: '//trim(vname)//' does not match expected value' + call extcde ( 49 ) + end if + + end subroutine checkfile + +end module wav_restart_mod diff --git a/model/src/wav_shel_inp.F90 b/model/src/wav_shel_inp.F90 index 1590b3b38..276196a03 100644 --- a/model/src/wav_shel_inp.F90 +++ b/model/src/wav_shel_inp.F90 @@ -52,6 +52,9 @@ subroutine set_shel_io(stdout,mds,ntrace) integer , intent(in) :: stdout integer , intent(out) :: mds(13), ntrace(2) + ! local variables + integer :: i + ! Note that nds is set to mds in w3initmd.F90 - mds is a local array ! The following units are referenced in module w3initmd ! NDS(1) ! OUTPUT LOG: General output unit number ("log file") @@ -80,17 +83,13 @@ subroutine set_shel_io(stdout,mds,ntrace) ! By default, unit numbers between 50 and 99 are scanned to find an ! unopened unit number - call ESMF_UtilIOUnitGet(mds(5)) ; open(unit=mds(5) , status='scratch') - call ESMF_UtilIOUnitGet(mds(6)) ; open(unit=mds(6) , status='scratch') - call ESMF_UtilIOUnitGet(mds(7)) ; open(unit=mds(7) , status='scratch') - call ESMF_UtilIOUnitGet(mds(8)) ; open(unit=mds(8) , status='scratch') - call ESMF_UtilIOUnitGet(mds(9)) ; open(unit=mds(9) , status='scratch') - call ESMF_UtilIOUnitGet(mds(10)); open(unit=mds(10) , status='scratch') - call ESMF_UtilIOUnitGet(mds(11)); open(unit=mds(11) , status='scratch') - call ESMF_UtilIOUnitGet(mds(12)); open(unit=mds(12) , status='scratch') - call ESMF_UtilIOUnitGet(mds(13)); open(unit=mds(13) , status='scratch') - close(mds(5)); close(mds(6)); close(mds(7)); close(mds(8)); close(mds(9)); close(mds(10)) - close(mds(11)); close(mds(12)); close(mds(13)) + do i = 5,size(mds) + call ESMF_UtilIOUnitGet(mds(i)) + open(unit=mds(i), status='scratch') + end do + do i = 5,size(mds) + close(mds(i)) + end do ntrace(1) = mds(3) ntrace(2) = 10 @@ -101,10 +100,14 @@ end subroutine set_shel_io !> Read ww3_shel.inp Or ww3_shel.nml !! !! @param[in] mpi_comm mpi communicator + !! @param[in] mds an array of unit numbers + !! @param[in] time0_overwrite the initial time for overwriting the nml file, optional + !! @param[in] timen_overwrite the endding time for overwriting the nml file, optional + !! @param[out] rstfldlist a list of additional restart fields, optional !! !> @author mvertens@ucar.edu, Denise.Worthen@noaa.gov !> @date 01-05-2022 - subroutine read_shel_config(mpi_comm, mds, time0_overwrite, timen_overwrite) + subroutine read_shel_config(mpi_comm, mds, time0_overwrite, timen_overwrite, rstfldlist) use wav_shr_flags use w3nmlshelmd , only : nml_domain_t, nml_input_t, nml_output_type_t @@ -128,12 +131,14 @@ subroutine read_shel_config(mpi_comm, mds, time0_overwrite, timen_overwrite) #ifdef W3_NL5 use w3wdatmd , only : qi5tbeg #endif + use wav_kind_mod , only : CL => shr_kind_cl ! input/output parameters - integer, intent(in) :: mpi_comm - integer, intent(in) :: mds(:) - integer, intent(in), optional :: time0_overwrite(2) - integer, intent(in), optional :: timen_overwrite(2) + integer, intent(in) :: mpi_comm + integer, intent(in) :: mds(:) + integer, intent(in), optional :: time0_overwrite(2) + integer, intent(in), optional :: timen_overwrite(2) + character(len=CL), intent(out), optional :: rstfldlist ! local parameters integer, parameter :: nhmax = 200 @@ -204,8 +209,9 @@ subroutine read_shel_config(mpi_comm, mds, time0_overwrite, timen_overwrite) memunit = 740+IAPROC call print_logmsg(740+IAPROC, 'read_shel_config, step 1', w3_debuginit_flag) - ! ndso, ndse, ndst are set in w3initmd using mds; w3initmd is called by either - ! cesm_init or uwm_int after calling the read_shel_config routine + ! module variables ndso, ndse, ndst are set in w3initmd using mds; w3initmd is + ! called by either cesm_init or uwm_int after calling the read_shel_config routine. + ! these nd units are local variables here ndso = mds(1) ndse = mds(1) ndst = mds(1) @@ -639,6 +645,13 @@ subroutine read_shel_config(mpi_comm, mds, time0_overwrite, timen_overwrite) ! Extra fields to be written in the restart fldrst = nml_output_type%restart%extra call w3flgrdflag ( ndso, ndso, ndse, fldrst, flogr, flogrr, iaproc, napout, ierr ) + if (present(rstfldlist)) then + if (trim(fldrst) .ne. 'unset')then + rstfldlist = trim(fldrst) + else + rstfldlist = ' ' + end if + end if if ( ierr .ne. 0 ) goto 2222 ! force minimal allocation to avoid memory seg fault diff --git a/model/src/wav_shr_mod.F90 b/model/src/wav_shr_mod.F90 index 9bddd503e..e0a3ccc8f 100644 --- a/model/src/wav_shr_mod.F90 +++ b/model/src/wav_shr_mod.F90 @@ -47,6 +47,7 @@ module wav_shr_mod private :: field_getfldptr !< @private obtain a pointer to a field public :: diagnose_mesh !< @public write out info about mesh public :: write_meshdecomp !< @public write the mesh decomposition to a file + public :: wav_loginit !< @public write the verbose WW3 log header interface state_getfldptr module procedure state_getfldptr_1d @@ -1343,6 +1344,28 @@ subroutine ymd2date_long(year,month,day,date) if (year < 0) date = -date end subroutine ymd2date_long + !=============================================================================== + !> Write the verbose WW3 log header + !! + !! @param[in] stdout the logfile unit on the root task + !! + !> @author Denise.Worthen@noaa.gov + !> @date 09-14-2024 + + subroutine wav_loginit(stdout) + + integer, intent(in) :: stdout + + write(stdout,984) +984 format (// & + 37x,'| input | output |'/ & + 37x,'|-----------------------|------------------|'/ & + 2x,' step | pass | date time | b w l c t r i i1 i5 d | g p t r b f c r2 |'/ & + 2x,'--------|------|---------------------|-----------------------|------------------|'/ & + 2x,'--------+------+---------------------+---------------------------+--------------+') + + end subroutine wav_loginit + !=============================================================================== !> Return a logical true if ESMF_LogFoundError detects an error !! diff --git a/model/src/wmwavemd.F90 b/model/src/wmwavemd.F90 index 99eec5eca..6dffa68b2 100644 --- a/model/src/wmwavemd.F90 +++ b/model/src/wmwavemd.F90 @@ -267,7 +267,7 @@ SUBROUTINE WMWAVE ( TEND ) !/ INTEGER :: J, JJ, I, JO, TPRNT(2), TAUX(2), & II, JJJ, IX, IY, UPNEXT(2), UPLAST(2) - INTEGER :: DUMMY2(35)=0 + INTEGER :: DUMMY2(40)=0 #ifdef W3_T INTEGER :: ILOOP #endif