From a2fd06e094ec2f9dcf08c47f143fdbfbb42ff665 Mon Sep 17 00:00:00 2001 From: "denise.worthen" Date: Wed, 31 Jan 2024 16:00:44 -0500 Subject: [PATCH 01/16] fix units and directions for netcdf gridded output --- model/src/w3iogoncdmd.F90 | 43 ++++++++++++++++++++++++++++++--------- model/src/wav_grdout.F90 | 6 +++--- 2 files changed, 36 insertions(+), 13 deletions(-) diff --git a/model/src/w3iogoncdmd.F90 b/model/src/w3iogoncdmd.F90 index 7520a0aa2..5f714e155 100644 --- a/model/src/w3iogoncdmd.F90 +++ b/model/src/w3iogoncdmd.F90 @@ -8,6 +8,7 @@ module w3iogoncdmd + use constants , only : rade use w3gdatmd , only : nk, nx, ny, mapsf, mapsta, nsea use w3odatmd , only : noswll, undef use w3odatmd , only : nds, iaproc, napout @@ -242,10 +243,10 @@ subroutine w3iogoncd () 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. 'PDIR') call write_var3d(vname, pdir (1:nsea,0:noswll), fldir='true' ) + if(vname .eq. 'PSI') call write_var3d(vname, psi (1:nsea,0:noswll), fldir='true' ) 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. 'PDP') call write_var3d(vname, pthp0 (1:nsea,0:noswll), fldir='true' ) 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) ) @@ -303,9 +304,9 @@ subroutine w3iogoncd () 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. 'THM') call write_var2d(vname, thm (1:nsea), fldir='true' ) + if (vname .eq. 'THS') call write_var2d(vname, ths (1:nsea), fldir='true' ) + if (vname .eq. 'THP0') call write_var2d(vname, thp0 (1:nsea), fldir='true' ) 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) ) @@ -399,7 +400,7 @@ subroutine w3iogoncd () end subroutine w3iogoncd !=============================================================================== - subroutine write_var2d(vname, var, dir, usemask, init0, init2) + subroutine write_var2d(vname, var, dir, usemask, init0, init2, fldir) ! 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 @@ -407,6 +408,7 @@ subroutine write_var2d(vname, var, dir, usemask, init0, init2) ! 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 + ! if fldir is true then the directions will be converted to degrees character(len=*), intent(in) :: vname real , intent(in) :: var(:) @@ -414,10 +416,11 @@ subroutine write_var2d(vname, var, dir, usemask, init0, init2) character(len=*), optional, intent(in) :: usemask character(len=*), optional, intent(in) :: init0 character(len=*), optional, intent(in) :: init2 + character(len=*), optional, intent(in) :: fldir ! local variables real, dimension(nx,ny) :: var2d - logical :: lmask, linit0, linit2 + logical :: lmask, linit0, linit2, lfldir real :: varloc lmask = .false. @@ -432,6 +435,10 @@ subroutine write_var2d(vname, var, dir, usemask, init0, init2) if (present(init2)) then linit2 = (trim(init2) == "true") end if + lfldir = .false. + if (present(fldir)) then + lfldir = (trim(fldir) == "true") + end if ! DEBUG ! write(nds(1),'(a)')' writing variable ' //trim(vname)//' to history file '//trim(fname) @@ -448,6 +455,11 @@ subroutine write_var2d(vname, var, dir, usemask, init0, init2) 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 @@ -474,24 +486,30 @@ subroutine write_var2d(vname, var, dir, usemask, init0, init2) end subroutine write_var2d !=============================================================================== - subroutine write_var3d(vname, var, init2) + subroutine write_var3d(vname, var, init2, fldir) ! 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 + ! if fldir is true then the directions will be converted to degrees 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 + 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) @@ -509,6 +527,11 @@ subroutine write_var3d(vname, var, init2) 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)).ge.0 ) then + varloc(:) = mod(630. - rade*varloc(:), 360.) + end if + end if var3d(mapsf(isea,1),mapsf(isea,2),:) = varloc(:) end do diff --git a/model/src/wav_grdout.F90 b/model/src/wav_grdout.F90 index 7e592e618..f93478330 100644 --- a/model/src/wav_grdout.F90 +++ b/model/src/wav_grdout.F90 @@ -157,9 +157,9 @@ subroutine initialize_gridout 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( "DIR ", "THM ", "Mean wave direction ", "deg ", " ", .false.) , & + varatts( "SPR ", "THS ", "Mean directional spread ", "deg ", " ", .false.) , & + varatts( "DP ", "THP0 ", "Peak direction ", "deg ", " ", .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.) , & From e8b3ef615eb332198846cd7848d093672504cf30 Mon Sep 17 00:00:00 2001 From: Ufuk Turuncoglu Date: Mon, 15 Jul 2024 16:23:01 -0500 Subject: [PATCH 02/16] fix nuopc cap for radiation stress components --- model/src/wav_import_export.F90 | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/model/src/wav_import_export.F90 b/model/src/wav_import_export.F90 index e6caa1a2d..4a9af31d7 100644 --- a/model/src/wav_import_export.F90 +++ b/model/src/wav_import_export.F90 @@ -145,6 +145,9 @@ subroutine advertise_fields(importState, ExportState, flds_scalar_name, rc) call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_vstokes') else call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_z0') + call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_wavsuu') + call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_wavsuv') + call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_wavsvv') end if call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_pstokes_x', ungridded_lbound=1, ungridded_ubound=3) call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_pstokes_y', ungridded_lbound=1, ungridded_ubound=3) @@ -753,14 +756,14 @@ subroutine export_fields (gcomp, rc) call CalcBotcur( va, wbcuru, wbcurv, wbcurp) end if - if ( state_fldchk(exportState, 'wavsuu') .and. & - state_fldchk(exportState, 'wavsuv') .and. & - state_fldchk(exportState, 'wavsvv')) then - call state_getfldptr(exportState, 'sxxn', sxxn, rc=rc) + if ( state_fldchk(exportState, 'Sw_wavsuu') .and. & + state_fldchk(exportState, 'Sw_wavsuv') .and. & + state_fldchk(exportState, 'Sw_wavsvv')) then + call state_getfldptr(exportState, 'Sw_wavsuu', sxxn, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getfldptr(exportState, 'sxyn', sxyn, rc=rc) + call state_getfldptr(exportState, 'Sw_wavsuv', sxyn, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getfldptr(exportState, 'syyn', syyn, rc=rc) + call state_getfldptr(exportState, 'Sw_wavsvv', syyn, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call CalcRadstr2D( va, sxxn, sxyn, syyn) end if From f13050bba3ff566d96fe1e180bda1db92eae0914 Mon Sep 17 00:00:00 2001 From: Ufuk Turuncoglu Date: Mon, 22 Jul 2024 11:41:03 -0500 Subject: [PATCH 03/16] fix type in the radiation stress components --- model/src/wav_grdout.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/model/src/wav_grdout.F90 b/model/src/wav_grdout.F90 index 4583070d7..4cde5e7e0 100644 --- a/model/src/wav_grdout.F90 +++ b/model/src/wav_grdout.F90 @@ -223,8 +223,8 @@ subroutine initialize_gridout ! 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( "SXX ", "SXX ", "Radiation stresses xx ", "N m-1 ", " ", .false.) , & + varatts( "SYY ", "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.) , & From 5debd4200d0334390d04ca38580026fb86ff4f09 Mon Sep 17 00:00:00 2001 From: Ufuk Turuncoglu Date: Mon, 29 Jul 2024 12:24:49 -0500 Subject: [PATCH 04/16] fix ww3 restart issue --- model/src/w3iorsmd.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/model/src/w3iorsmd.F90 b/model/src/w3iorsmd.F90 index 24d9a280c..6bc959baf 100644 --- a/model/src/w3iorsmd.F90 +++ b/model/src/w3iorsmd.F90 @@ -1376,6 +1376,7 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT ) TICE(1) = -1 TICE(2) = 0 TRHO(1) = -1 + TRHO(2) = 0 TIC1(1) = -1 TIC1(2) = 0 TIC5(1) = -1 From 6064dfc04b7657957b3c303e3cec757888d3f331 Mon Sep 17 00:00:00 2001 From: Denise Worthen Date: Wed, 31 Jul 2024 09:05:33 -0400 Subject: [PATCH 05/16] fix time axis in netcdf output files --- model/src/w3iogoncdmd.F90 | 10 ++++++---- model/src/w3wavemd.F90 | 2 +- model/src/wav_comp_nuopc.F90 | 2 +- 3 files changed, 8 insertions(+), 6 deletions(-) diff --git a/model/src/w3iogoncdmd.F90 b/model/src/w3iogoncdmd.F90 index 813aa28d2..7e85ca84a 100644 --- a/model/src/w3iogoncdmd.F90 +++ b/model/src/w3iogoncdmd.F90 @@ -36,7 +36,7 @@ module w3iogoncdmd contains !=============================================================================== - subroutine w3iogoncd () + subroutine w3iogoncd ( timen ) use w3odatmd , only : fnmpre use w3gdatmd , only : filext, trigp, ntri, ungtype, gtype @@ -60,7 +60,7 @@ subroutine w3iogoncd () 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 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 @@ -68,6 +68,8 @@ subroutine w3iogoncd () !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 :: igrd integer ,target :: dimid3(3) @@ -96,10 +98,10 @@ subroutine w3iogoncd () 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) + call set_user_timestring(timen,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' + write(fname,'(a,i8.8,a1,i6.6,a)')trim(fnmpre),timen(1),'.',timen(2),'.out_grd.'//trim(filext)//'.nc' end if len_s = noswll + 1 ! 0:noswll diff --git a/model/src/w3wavemd.F90 b/model/src/w3wavemd.F90 index b11aeb3b2..3a5760696 100644 --- a/model/src/w3wavemd.F90 +++ b/model/src/w3wavemd.F90 @@ -2606,7 +2606,7 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & IF ( FLGMPI(1) ) CALL MPI_WAITALL( NRQGO2, IRQGO2, STATIO, IERR_MPI ) FLGMPI(1) = .FALSE. #endif - CALL W3IOGONCD () + CALL W3IOGONCD ( tend ) END IF else ! default (binary) output diff --git a/model/src/wav_comp_nuopc.F90 b/model/src/wav_comp_nuopc.F90 index c7bb14ef2..f2d4a4f83 100644 --- a/model/src/wav_comp_nuopc.F90 +++ b/model/src/wav_comp_nuopc.F90 @@ -623,7 +623,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) From 1f52e1fd2d0637d8bdce7e76b1e1bb1d5b7e0787 Mon Sep 17 00:00:00 2001 From: "denise.worthen" Date: Mon, 5 Aug 2024 08:37:07 -0500 Subject: [PATCH 06/16] use new radstr routine to match history file values --- model/src/wav_import_export.F90 | 167 +++++++++++++++++++++++--------- 1 file changed, 120 insertions(+), 47 deletions(-) diff --git a/model/src/wav_import_export.F90 b/model/src/wav_import_export.F90 index 4a9af31d7..3821e13f2 100644 --- a/model/src/wav_import_export.F90 +++ b/model/src/wav_import_export.F90 @@ -94,8 +94,7 @@ subroutine advertise_fields(importState, ExportState, flds_scalar_name, rc) integer , intent(out) :: rc ! local variables - integer :: n, num - character(len=2) :: fvalue + integer :: n character(len=*), parameter :: subname='(wav_import_export:advertise_fields)' !------------------------------------------------------------------------------- @@ -138,16 +137,18 @@ subroutine advertise_fields(importState, ExportState, flds_scalar_name, rc) if (.not. unstr_mesh) then call fldlist_add(fldsFrWav_num, fldsFrWav, trim(flds_scalar_name)) end if + call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_ustokes') + call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_vstokes') + if (cesmcoupled) then call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_lamult' ) call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_lasl' ) - call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_ustokes') - call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_vstokes') else call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_z0') call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_wavsuu') call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_wavsuv') call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_wavsvv') + call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_hs') end if call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_pstokes_x', ungridded_lbound=1, ungridded_ubound=3) call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_pstokes_y', ungridded_lbound=1, ungridded_ubound=3) @@ -583,7 +584,7 @@ subroutine export_fields (gcomp, rc) !--------------------------------------------------------------------------- use wav_kind_mod, only : R8 => SHR_KIND_R8 - use w3adatmd , only : USSX, USSY, USSP + use w3adatmd , only : USSX, USSY, USSP, HS use w3adatmd , only : w3seta use w3idatmd , only : w3seti use w3wdatmd , only : va, w3setw @@ -592,7 +593,7 @@ subroutine export_fields (gcomp, rc) use w3iogomd , only : CALC_U3STOKES #ifdef W3_CESMCOUPLED use w3wdatmd , only : ASF, UST - use w3adatmd , only : USSHX, USSHY, UD, HS + use w3adatmd , only : USSHX, USSHY, UD use w3idatmd , only : HSL #else use wmmdatmd , only : mdse, mdst, wmsetm @@ -625,6 +626,7 @@ subroutine export_fields (gcomp, rc) real(r8), pointer :: sw_lasl(:) real(r8), pointer :: sw_ustokes(:) real(r8), pointer :: sw_vstokes(:) + real(r8), pointer :: sw_hs(:) ! d2 is location, d1 is frequency - nwav_elev_spectrum frequencies will be used real(r8), pointer :: wave_elevation_spectrum(:,:) @@ -767,6 +769,7 @@ subroutine export_fields (gcomp, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call CalcRadstr2D( va, sxxn, sxyn, syyn) end if + if (wav_coupling_to_cice) then call state_getfldptr(exportState, 'Sw_elevation_spectrum', wave_elevation_spectrum, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -798,6 +801,13 @@ subroutine export_fields (gcomp, rc) end if endif + if (state_fldchk(exportState, 'Sw_hs')) then + call state_getfldptr(exportState, 'Sw_hs', sw_hs, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + sw_hs(:) = fillvalue + call CalcHs(va, sw_hs) + end if + if (dbug_flag > 5) then call state_diagnose(exportState, 'at export ', rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -1182,15 +1192,16 @@ end subroutine CalcBotcur !=============================================================================== !> Calculate radiation stresses for export !! - !> @details TODO: + !> @details Calculates radiation stresses independently of w3iogomd to ensure + !! that export fields are updated at the coupling frequency !! !! @param[in] a input spectra !! @param sxxn a 1-D pointer to a field on a mesh !! @param sxyn a 1-D pointer to a field on a mesh !! @param syyn a 1-D pointer to a field on a mesh !! - !> @author T. J. Campbell, NRL - !> @date 09-Aug-2017 + !> @author Denise.Worthen@noaa.gov + !> @date 08-05-2024 subroutine CalcRadstr2D ( a, sxxn, sxyn, syyn ) ! Calculate radiation stresses for export @@ -1206,45 +1217,42 @@ subroutine CalcRadstr2D ( a, sxxn, sxyn, syyn ) real(ESMF_KIND_R8), pointer :: syyn(:) ! northward-component export field ! local variables - character(ESMF_MAXSTR) :: cname - character(128) :: msg - real(8), parameter :: half = 0.5 - real(8), parameter :: one = 1.0 - real(8), parameter :: two = 2.0 - integer :: isea, jsea, ik, ith - real(8) :: sxxs, sxys, syys - real(8) :: akxx, akxy, akyy, cgoc, facd, fack, facs - !---------------------------------------------------------------------- + integer :: isea, jsea, ik, ith + real :: factor, abxx, abyy, abxy - facd = dwat*grav - jsea_loop: do jsea = 1,nseal_cpl + sxxn = 0.0 + syyn = 0.0 + sxyn = 0.0 + do jsea = 1,nseal_cpl call init_get_isea(isea, jsea) - if ( dw(isea).le.zero ) cycle jsea_loop - sxxs = zero - sxys = zero - syys = zero - ik_loop: do ik = 1,nk - akxx = zero - akxy = zero - akyy = zero - cgoc = cg(ik,isea)*wn(ik,isea)/sig(ik) - cgoc = min(one,max(half,cgoc)) - ith_loop: do ith = 1,nth - akxx = akxx + (cgoc*(ec2(ith)+one)-half)*a(ith,ik,jsea) - akxy = akxy + cgoc*esc(ith)*a(ith,ik,jsea) - akyy = akyy + (cgoc*(es2(ith)+one)-half)*a(ith,ik,jsea) - enddo ith_loop - fack = dden(ik)/cg(ik,isea) - sxxs = sxxs + akxx*fack - sxys = sxys + akxy*fack - syys = syys + akyy*fack - enddo ik_loop - facs = (one+fte/cg(nk,isea))*facd - sxxn(jsea) = sxxs*facs - sxyn(jsea) = sxys*facs - syyn(jsea) = syys*facs - enddo jsea_loop + do ik = 1,nk + factor = max ( 0.5, cg(ik,isea)/sig(ik)*wn(ik,isea) ) + abxx = 0.0 + abyy = 0.0 + abxy = 0.0 + do ith = 1,nth + abxx = abxx + ((1.0 + ec2(ith))*factor-0.5) * a(ith,ik,jsea) + abyy = abyy + ((1.0 + es2(ith))*factor-0.5) * a(ith,ik,jsea) + abxy = abxy + esc(ith)* factor * a(ith,ik,jsea) + end do + + factor = dden(ik) / cg(ik,isea) + abxx = max ( 0.0, abxx ) * factor + abyy = max ( 0.0, abyy ) * factor + abxy = abxy * factor + + sxxn(jsea) = sxxn(jsea) + abxx + syyn(jsea) = syyn(jsea) + abyy + sxyn(jsea) = sxyn(jsea) + abxy + end do + sxxn(jsea) = sxxn(jsea) + fte * abxx/cg(nk,isea) + syyn(jsea) = syyn(jsea) + fte * abyy/cg(nk,isea) + sxyn(jsea) = sxyn(jsea) + fte * abxy/cg(nk,isea) + end do + sxxn = sxxn*dwat*grav + syyn = syyn*dwat*grav + sxyn = sxyn*dwat*grav end subroutine CalcRadstr2D !=============================================================================== @@ -1299,6 +1307,71 @@ subroutine CalcEF (a, wave_elevation_spectrum) end subroutine CalcEF + !=============================================================================== + !> Calculate significant wave height for export + !! + !> @details Calculates significant wave height independently of w3iogomd to ensure + !! that exported HS field is updated at the coupling frequency + !! + !! @param[in] a input spectra + !! @param[inout] hs significant wave height + !! + !> @author Denise.Worthen@noaa.gov + !> @date 8-02-2024 + subroutine CalcHs (a, hs) + + use constants, only : tpi + use w3gdatmd, only : nth, nk, nseal, mapsf, mapsta, dden, fte + use w3adatmd, only : nsealm, cg + use w3parall, only : init_get_isea + + ! input/output variables + real, intent(in) :: a(nth,nk,0:nseal) + real(r8), pointer :: hs(:) + + ! local variables + real :: ab(nseal), et(nseal) + real :: ebd, factor, eband + integer :: ik, ith, isea, jsea, ix, iy + + et = 0.0 + do ik = 1, nk + ab = 0.0 + do ith = 1, nth + do jsea = 1,nseal_cpl + ab(jsea) = ab(jsea) + a(ith,ik,jsea) + end do + end do + + do jsea = 1,nseal_cpl + call init_get_isea(isea, jsea) + factor = dden(ik) / cg(ik,isea) + ebd = ab(jsea) * factor + et(jsea) = et(jsea) + ebd + end do + end do !ik + + do jsea = 1,nseal_cpl + call init_get_isea(isea, jsea) + eband = ab(jsea)/cg(nk,isea) + et(jsea) = et(jsea) + fte*eband + ix = mapsf(isea,1) ! global ix + iy = mapsf(isea,2) ! global iy + if (mapsta(iy,ix) == 1) then ! active sea point +#ifdef W3_O9 + if ( et(jsea) .ge. 0.0 ) then + hs(jsea) = 4.0*sqrt ( et(jsea) ) + else + hs(jsea) = -4.0*sqrt ( -et(jsea) ) + end if +#else + hs(jsea) = 4.0*sqrt ( et(jsea) ) +#endif + end if + end do + + end subroutine CalcHs + !==================================================================================== !> Create a global field across all PEs !! @@ -1315,7 +1388,7 @@ end subroutine CalcEF !> @date 01-05-2022 subroutine SetGlobalInput(importState, fldname, vm, global_output, rc) - use w3gdatmd, only: nsea, nseal, nx, ny + use w3gdatmd, only: nsea, nx, ny use w3odatmd, only: naproc, iaproc ! input/output variables @@ -1529,7 +1602,7 @@ end subroutine set_importmask !> @date 01-05-2022 subroutine check_globaldata(gcomp, fldname, global_data, nvals, rc) - use w3gdatmd, only: nseal, nsea, mapsf, nx, ny + use w3gdatmd, only: nseal, mapsf, nx, ny use w3odatmd, only: naproc, iaproc ! input/output variables From e236af558bf87b763259f26d93b401da2e33bfc0 Mon Sep 17 00:00:00 2001 From: "denise.worthen" Date: Mon, 5 Aug 2024 12:32:50 -0500 Subject: [PATCH 07/16] fixes and placeholders * hs and 2dstr in export fields match history files * add placeholders for additional fields --- model/src/wav_import_export.F90 | 120 ++++++++++++++++++-------------- 1 file changed, 66 insertions(+), 54 deletions(-) diff --git a/model/src/wav_import_export.F90 b/model/src/wav_import_export.F90 index 3821e13f2..79b1dc125 100644 --- a/model/src/wav_import_export.F90 +++ b/model/src/wav_import_export.F90 @@ -145,10 +145,18 @@ subroutine advertise_fields(importState, ExportState, flds_scalar_name, rc) call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_lasl' ) else call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_z0') + ! coastal coupling call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_wavsuu') call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_wavsuv') call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_wavsvv') call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_hs') + call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_bhd') + call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_tauox') + call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_tauoy') + call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_taubblx') + call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_taubbly') + call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_ubax') + call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_ubay') end if call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_pstokes_x', ungridded_lbound=1, ungridded_ubound=3) call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_pstokes_y', ungridded_lbound=1, ungridded_ubound=3) @@ -618,15 +626,22 @@ subroutine export_fields (gcomp, rc) real(r8), pointer :: wbcuru(:) real(r8), pointer :: wbcurv(:) real(r8), pointer :: wbcurp(:) - real(r8), pointer :: sxxn(:) - real(r8), pointer :: sxyn(:) - real(r8), pointer :: syyn(:) - real(r8), pointer :: sw_lamult(:) real(r8), pointer :: sw_lasl(:) real(r8), pointer :: sw_ustokes(:) real(r8), pointer :: sw_vstokes(:) + + real(r8), pointer :: sxxn(:) + real(r8), pointer :: sxyn(:) + real(r8), pointer :: syyn(:) real(r8), pointer :: sw_hs(:) + real(r8), pointer :: sw_bhd(:) + real(r8), pointer :: sw_tauox(:) + real(r8), pointer :: sw_tauoy(:) + real(r8), pointer :: sw_taubblx(:) + real(r8), pointer :: sw_taubbly(:) + real(r8), pointer :: sw_ubax(:) + real(r8), pointer :: sw_ubay(:) ! d2 is location, d1 is frequency - nwav_elev_spectrum frequencies will be used real(r8), pointer :: wave_elevation_spectrum(:,:) @@ -808,6 +823,7 @@ subroutine export_fields (gcomp, rc) call CalcHs(va, sw_hs) end if + if (dbug_flag > 5) then call state_diagnose(exportState, 'at export ', rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -1206,7 +1222,7 @@ subroutine CalcRadstr2D ( a, sxxn, sxyn, syyn ) ! Calculate radiation stresses for export - use w3gdatmd, only : nseal, nk, nth, sig, es2, esc, ec2, fte, dden + use w3gdatmd, only : nseal, nk, nth, sig, es2, esc, ec2, fte, dden, mapsf, mapsta use w3adatmd, only : dw, cg, wn use w3odatmd, only : naproc, iaproc @@ -1217,7 +1233,7 @@ subroutine CalcRadstr2D ( a, sxxn, sxyn, syyn ) real(ESMF_KIND_R8), pointer :: syyn(:) ! northward-component export field ! local variables - integer :: isea, jsea, ik, ith + integer :: isea, jsea, ik, ith, ix, iy real :: factor, abxx, abyy, abxy sxxn = 0.0 @@ -1225,34 +1241,38 @@ subroutine CalcRadstr2D ( a, sxxn, sxyn, syyn ) sxyn = 0.0 do jsea = 1,nseal_cpl call init_get_isea(isea, jsea) - do ik = 1,nk - factor = max ( 0.5, cg(ik,isea)/sig(ik)*wn(ik,isea) ) - abxx = 0.0 - abyy = 0.0 - abxy = 0.0 - do ith = 1,nth - abxx = abxx + ((1.0 + ec2(ith))*factor-0.5) * a(ith,ik,jsea) - abyy = abyy + ((1.0 + es2(ith))*factor-0.5) * a(ith,ik,jsea) - abxy = abxy + esc(ith)* factor * a(ith,ik,jsea) - end do + ix = mapsf(isea,1) ! global ix + iy = mapsf(isea,2) ! global iy + if (mapsta(iy,ix) == 1) then ! active sea point + do ik = 1,nk + factor = max ( 0.5, cg(ik,isea)/sig(ik)*wn(ik,isea) ) + abxx = 0.0 + abyy = 0.0 + abxy = 0.0 + do ith = 1,nth + abxx = abxx + ((1.0 + ec2(ith))*factor-0.5) * a(ith,ik,jsea) + abyy = abyy + ((1.0 + es2(ith))*factor-0.5) * a(ith,ik,jsea) + abxy = abxy + esc(ith)* factor * a(ith,ik,jsea) + end do - factor = dden(ik) / cg(ik,isea) - abxx = max ( 0.0, abxx ) * factor - abyy = max ( 0.0, abyy ) * factor - abxy = abxy * factor + factor = dden(ik) / cg(ik,isea) + abxx = max ( 0.0, abxx ) * factor + abyy = max ( 0.0, abyy ) * factor + abxy = abxy * factor - sxxn(jsea) = sxxn(jsea) + abxx - syyn(jsea) = syyn(jsea) + abyy - sxyn(jsea) = sxyn(jsea) + abxy - end do - sxxn(jsea) = sxxn(jsea) + fte * abxx/cg(nk,isea) - syyn(jsea) = syyn(jsea) + fte * abyy/cg(nk,isea) - sxyn(jsea) = sxyn(jsea) + fte * abxy/cg(nk,isea) + sxxn(jsea) = sxxn(jsea) + abxx + syyn(jsea) = syyn(jsea) + abyy + sxyn(jsea) = sxyn(jsea) + abxy + end do + sxxn(jsea) = sxxn(jsea) + fte * abxx/cg(nk,isea) + syyn(jsea) = syyn(jsea) + fte * abyy/cg(nk,isea) + sxyn(jsea) = sxyn(jsea) + fte * abxy/cg(nk,isea) + end if end do - sxxn = sxxn*dwat*grav syyn = syyn*dwat*grav sxyn = sxyn*dwat*grav + end subroutine CalcRadstr2D !=============================================================================== @@ -1314,7 +1334,7 @@ end subroutine CalcEF !! that exported HS field is updated at the coupling frequency !! !! @param[in] a input spectra - !! @param[inout] hs significant wave height + !! @param[inout] hs a 1-D pointer to a field on a mesh !! !> @author Denise.Worthen@noaa.gov !> @date 8-02-2024 @@ -1330,46 +1350,38 @@ subroutine CalcHs (a, hs) real(r8), pointer :: hs(:) ! local variables - real :: ab(nseal), et(nseal) - real :: ebd, factor, eband + real :: factor, eband, ab, et integer :: ik, ith, isea, jsea, ix, iy et = 0.0 - do ik = 1, nk - ab = 0.0 - do ith = 1, nth - do jsea = 1,nseal_cpl - ab(jsea) = ab(jsea) + a(ith,ik,jsea) - end do - end do - - do jsea = 1,nseal_cpl - call init_get_isea(isea, jsea) - factor = dden(ik) / cg(ik,isea) - ebd = ab(jsea) * factor - et(jsea) = et(jsea) + ebd - end do - end do !ik - + ab = 0.0 do jsea = 1,nseal_cpl call init_get_isea(isea, jsea) - eband = ab(jsea)/cg(nk,isea) - et(jsea) = et(jsea) + fte*eband ix = mapsf(isea,1) ! global ix iy = mapsf(isea,2) ! global iy if (mapsta(iy,ix) == 1) then ! active sea point + et = 0.0 + do ik = 1,nk + factor = dden(ik) / cg(ik,isea) + ab = 0.0 + do ith = 1,nth + ab = ab + a(ith,ik,jsea) + end do + et = et + ab*factor + end do !ik + eband = ab/cg(nk,isea) + et = et + fte*eband #ifdef W3_O9 - if ( et(jsea) .ge. 0.0 ) then - hs(jsea) = 4.0*sqrt ( et(jsea) ) + if ( et .ge. 0.0 ) then + hs(jsea) = 4.0*sqrt ( et ) else - hs(jsea) = -4.0*sqrt ( -et(jsea) ) + hs(jsea) = -4.0*sqrt ( -et ) end if #else - hs(jsea) = 4.0*sqrt ( et(jsea) ) + hs(jsea) = 4.0*sqrt ( et ) #endif end if end do - end subroutine CalcHs !==================================================================================== From 213dc9ce8cf6db1e930dc66027c884dfaa35673f Mon Sep 17 00:00:00 2001 From: "denise.worthen" Date: Mon, 5 Aug 2024 15:51:14 -0500 Subject: [PATCH 08/16] additional fields, with notes --- model/src/wav_import_export.F90 | 283 +++++++++++++++++++++----------- 1 file changed, 187 insertions(+), 96 deletions(-) diff --git a/model/src/wav_import_export.F90 b/model/src/wav_import_export.F90 index 79b1dc125..ace994f96 100644 --- a/model/src/wav_import_export.F90 +++ b/model/src/wav_import_export.F90 @@ -592,7 +592,7 @@ subroutine export_fields (gcomp, rc) !--------------------------------------------------------------------------- use wav_kind_mod, only : R8 => SHR_KIND_R8 - use w3adatmd , only : USSX, USSY, USSP, HS + use w3adatmd , only : USSX, USSY, USSP, HS, tauox, tauoy use w3adatmd , only : w3seta use w3idatmd , only : w3seti use w3wdatmd , only : va, w3setw @@ -717,6 +717,7 @@ subroutine export_fields (gcomp, rc) enddo end if #endif + ! TODO: needs to be at export freq ! surface stokes drift if (state_fldchk(exportState, 'Sw_ustokes')) then call state_getfldptr(exportState, 'Sw_ustokes', sw_ustokes, rc=rc) @@ -761,17 +762,17 @@ subroutine export_fields (gcomp, rc) call CalcRoughl(z0rlen) endif - if ( state_fldchk(exportState, 'wbcuru') .and. & - state_fldchk(exportState, 'wbcurv') .and. & - state_fldchk(exportState, 'wbcurp')) then - call state_getfldptr(exportState, 'wbcuru', wbcuru, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getfldptr(exportState, 'wbcurv', wbcurv, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getfldptr(exportState, 'wbcurp', wbcurp, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call CalcBotcur( va, wbcuru, wbcurv, wbcurp) - end if + ! if ( state_fldchk(exportState, 'wbcuru') .and. & + ! state_fldchk(exportState, 'wbcurv') .and. & + ! state_fldchk(exportState, 'wbcurp')) then + ! call state_getfldptr(exportState, 'wbcuru', wbcuru, rc=rc) + ! if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! call state_getfldptr(exportState, 'wbcurv', wbcurv, rc=rc) + ! if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! call state_getfldptr(exportState, 'wbcurp', wbcurp, rc=rc) + ! if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! call CalcBotcur( va, wbcuru, wbcurv, wbcurp) + ! end if if ( state_fldchk(exportState, 'Sw_wavsuu') .and. & state_fldchk(exportState, 'Sw_wavsuv') .and. & @@ -823,6 +824,46 @@ subroutine export_fields (gcomp, rc) call CalcHs(va, sw_hs) end if + if (state_fldchk(exportState, 'Sw_bhd')) then + call state_getfldptr(exportState, 'Sw_bhd', sw_bhd, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + sw_bhd(:) = fillvalue + call CalcBHD(va, sw_bhd) + end if + + if ( state_fldchk(exportState, 'Sw_tauox') .and. & + state_fldchk(exportState, 'Sw_tauoy') )then + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call state_getfldptr(exportState, 'Sw_tauox', sw_tauox, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call state_getfldptr(exportState, 'Sw_tauoy', sw_tauoy, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + sw_tauox(:) = fillvalue + sw_tauoy(:) = fillvalue + do jsea=1, nseal_cpl + call init_get_isea(isea, jsea) + ix = mapsf(isea,1) + iy = mapsf(isea,2) + if (mapsta(iy,ix) == 1) then + sw_tauox(jsea) = tauox(jsea) + sw_tauoy(jsea) = tauoy(jsea) + else + sw_tauox(jsea) = 0. + sw_tauoy(jsea) = 0. + endif + enddo + end if + + if ( state_fldchk(exportState, 'Sw_ubax') .and. & + state_fldchk(exportState, 'Sw_ubay') )then + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call state_getfldptr(exportState, 'Sw_ubax', sw_ubax, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call state_getfldptr(exportState, 'Sw_ubay', sw_ubay, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! TODO: + !call CalcUVB(va, sw_ubax, sw_ubay) + end if if (dbug_flag > 5) then call state_diagnose(exportState, 'at export ', rc=rc) @@ -1122,88 +1163,88 @@ subroutine CalcRoughl ( wrln) end subroutine CalcRoughl - !=============================================================================== - !> Calculate wave-bottom currents for export - !! - !> @details TODO: - !! - !! @param[in] a input spectra - !! @param wbxn a 1-D pointer to a field on a mesh - !! @param wbyn a 1-D pointer to a field on a mesh - !! @param wbpn a 1-D pointer to a field on a mesh - !! - !> @author T. J. Campbell, NRL - !> @date 09-Aug-2017 - subroutine CalcBotcur ( a, wbxn, wbyn, wbpn ) - - ! Calculate wave-bottom currents for export - - use w3gdatmd, only : nseal, nk, nth, sig, dmin, ecos, esin, dden, mapsf, mapsta, nspec - use w3adatmd, only : dw, cg, wn - use w3odatmd, only : naproc, iaproc - - ! input/output variables - real, intent(in) :: a(nth,nk,0:nseal) ! Input spectra (in par list to change shape) - real(ESMF_KIND_R8), pointer :: wbxn(:) ! eastward-component export field pointer - real(ESMF_KIND_R8), pointer :: wbyn(:) ! northward-component export field pointer - real(ESMF_KIND_R8), pointer :: wbpn(:) ! period export field pointer - - ! local variables - real(8), parameter :: half = 0.5_r8 - real(8), parameter :: one = 1.0_r8 - real(8), parameter :: two = 2.0_r8 - real(8), parameter :: kdmin = 1e-7_r8 - real(8), parameter :: kdmax = 18.0_r8 - integer :: isea, jsea, ik, ith - real(8) :: depth - real(8) :: kd, fack, fkd, aka, akx, aky, abr, ubr, ubx, uby, dir - real(8), allocatable :: sig2(:) - !---------------------------------------------------------------------- - - allocate( sig2(1:nk) ) - sig2(1:nk) = sig(1:nk)**2 - - wbxn(:) = zero - wbyn(:) = zero - wbpn(:) = zero - - jsea_loop: do jsea = 1,nseal_cpl - call init_get_isea(isea, jsea) - if ( dw(isea).le.zero ) cycle jsea_loop - depth = max(dmin,dw(isea)) - abr = zero - ubr = zero - ubx = zero - uby = zero - ik_loop: do ik = 1,nk - aka = zero - akx = zero - aky = zero - ith_loop: do ith = 1,nth - aka = aka + a(ith,ik,jsea) - akx = akx + a(ith,ik,jsea)*ecos(ith) - aky = aky + a(ith,ik,jsea)*esin(ith) - enddo ith_loop - fack = dden(ik)/cg(ik,isea) - kd = max(kdmin,min(kdmax,wn(ik,isea)*depth)) - fkd = fack/sinh(kd)**2 - abr = abr + aka*fkd - ubr = ubr + aka*sig2(ik)*fkd - ubx = ubx + akx*sig2(ik)*fkd - uby = uby + aky*sig2(ik)*fkd - enddo ik_loop - if ( abr.le.zero .or. ubr.le.zero ) cycle jsea_loop - abr = sqrt(two*abr) - ubr = sqrt(two*ubr) - dir = atan2(uby,ubx) - wbxn(jsea) = ubr*cos(dir) - wbyn(jsea) = ubr*sin(dir) - wbpn(jsea) = tpi*abr/ubr - enddo jsea_loop - - deallocate( sig2 ) - - end subroutine CalcBotcur + ! !=============================================================================== + ! !> Calculate wave-bottom currents for export + ! !! + ! !> @details TODO: + ! !! + ! !! @param[in] a input spectra + ! !! @param wbxn a 1-D pointer to a field on a mesh + ! !! @param wbyn a 1-D pointer to a field on a mesh + ! !! @param wbpn a 1-D pointer to a field on a mesh + ! !! + ! !> @author T. J. Campbell, NRL + ! !> @date 09-Aug-2017 + ! subroutine CalcBotcur ( a, wbxn, wbyn, wbpn ) + + ! ! Calculate wave-bottom currents for export + + ! use w3gdatmd, only : nseal, nk, nth, sig, dmin, ecos, esin, dden, mapsf, mapsta, nspec + ! use w3adatmd, only : dw, cg, wn + ! use w3odatmd, only : naproc, iaproc + + ! ! input/output variables + ! real, intent(in) :: a(nth,nk,0:nseal) ! Input spectra (in par list to change shape) + ! real(ESMF_KIND_R8), pointer :: wbxn(:) ! eastward-component export field pointer + ! real(ESMF_KIND_R8), pointer :: wbyn(:) ! northward-component export field pointer + ! real(ESMF_KIND_R8), pointer :: wbpn(:) ! period export field pointer + + ! ! local variables + ! real(8), parameter :: half = 0.5_r8 + ! real(8), parameter :: one = 1.0_r8 + ! real(8), parameter :: two = 2.0_r8 + ! real(8), parameter :: kdmin = 1e-7_r8 + ! real(8), parameter :: kdmax = 18.0_r8 + ! integer :: isea, jsea, ik, ith + ! real(8) :: depth + ! real(8) :: kd, fack, fkd, aka, akx, aky, abr, ubr, ubx, uby, dir + ! real(8), allocatable :: sig2(:) + ! !---------------------------------------------------------------------- + + ! allocate( sig2(1:nk) ) + ! sig2(1:nk) = sig(1:nk)**2 + + ! wbxn(:) = zero + ! wbyn(:) = zero + ! wbpn(:) = zero + + ! jsea_loop: do jsea = 1,nseal_cpl + ! call init_get_isea(isea, jsea) + ! if ( dw(isea).le.zero ) cycle jsea_loop + ! depth = max(dmin,dw(isea)) + ! abr = zero + ! ubr = zero + ! ubx = zero + ! uby = zero + ! ik_loop: do ik = 1,nk + ! aka = zero + ! akx = zero + ! aky = zero + ! ith_loop: do ith = 1,nth + ! aka = aka + a(ith,ik,jsea) + ! akx = akx + a(ith,ik,jsea)*ecos(ith) + ! aky = aky + a(ith,ik,jsea)*esin(ith) + ! enddo ith_loop + ! fack = dden(ik)/cg(ik,isea) + ! kd = max(kdmin,min(kdmax,wn(ik,isea)*depth)) + ! fkd = fack/sinh(kd)**2 + ! abr = abr + aka*fkd + ! ubr = ubr + aka*sig2(ik)*fkd + ! ubx = ubx + akx*sig2(ik)*fkd + ! uby = uby + aky*sig2(ik)*fkd + ! enddo ik_loop + ! if ( abr.le.zero .or. ubr.le.zero ) cycle jsea_loop + ! abr = sqrt(two*abr) + ! ubr = sqrt(two*ubr) + ! dir = atan2(uby,ubx) + ! wbxn(jsea) = ubr*cos(dir) + ! wbyn(jsea) = ubr*sin(dir) + ! wbpn(jsea) = tpi*abr/ubr + ! enddo jsea_loop + + ! deallocate( sig2 ) + + ! end subroutine CalcBotcur !=============================================================================== !> Calculate radiation stresses for export @@ -1220,8 +1261,6 @@ end subroutine CalcBotcur !> @date 08-05-2024 subroutine CalcRadstr2D ( a, sxxn, sxyn, syyn ) - ! Calculate radiation stresses for export - use w3gdatmd, only : nseal, nk, nth, sig, es2, esc, ec2, fte, dden, mapsf, mapsta use w3adatmd, only : dw, cg, wn use w3odatmd, only : naproc, iaproc @@ -1384,6 +1423,58 @@ subroutine CalcHs (a, hs) end do end subroutine CalcHs + !=============================================================================== + !> Calculate Bernoulli head pressure for export + !! + !> @details Calculates Bernoulli head pressure independently of w3iogomd to ensure + !! that exported BHD field is updated at the coupling frequency + !! + !! @param[in] a input spectra + !! @param[inout] bhd a 1-D pointer to a field on a mesh + !! + !> @author Denise.Worthen@noaa.gov + !> @date 8-02-2024 + subroutine CalcBHD (a, bhd) + + use constants, only : tpi + use w3gdatmd, only : nth, nk, nseal, mapsf, mapsta, dden, fte + use w3adatmd, only : dw, nsealm, cg, wn + use w3parall, only : init_get_isea + + ! input/output variables + real, intent(in) :: a(nth,nk,0:nseal) + real(r8), pointer :: bhd(:) + + ! local variables + real :: factor, kd, ab, ebd + integer :: ik, ith, isea, jsea, ix, iy + + ebd = 0.0 + ab = 0.0 + bhd = 0.0 + do jsea = 1,nseal_cpl + call init_get_isea(isea, jsea) + ix = mapsf(isea,1) ! global ix + iy = mapsf(isea,2) ! global iy + if (mapsta(iy,ix) == 1) then ! active sea point + ebd = 0.0 + do ik = 1,nk + factor = dden(ik) / cg(ik,isea) + ab = 0.0 + do ith = 1,nth + ab = ab + a(ith,ik,jsea) + end do + ebd = ab*factor + kd = max ( 0.001 , wn(ik,isea) * dw(isea) ) + if (kd .lt. 6.0) then + bhd(jsea) = bhd(jsea) + grav*wn(ik,isea) * ebd / (sinh(2.*kd)) + end if + end do !ik + end if + end do + + end subroutine CalcBHD + !==================================================================================== !> Create a global field across all PEs !! From c40e6037d808e05777ccdf5fa00952f7ba4a4063 Mon Sep 17 00:00:00 2001 From: "denise.worthen" Date: Tue, 6 Aug 2024 12:02:02 -0500 Subject: [PATCH 09/16] add stokes, allow for fill value * allow for masking of export variables via fillvalue --- model/src/wav_import_export.F90 | 204 ++++++++++++++++++++------------ 1 file changed, 126 insertions(+), 78 deletions(-) diff --git a/model/src/wav_import_export.F90 b/model/src/wav_import_export.F90 index ace994f96..e672afcbd 100644 --- a/model/src/wav_import_export.F90 +++ b/model/src/wav_import_export.F90 @@ -717,37 +717,14 @@ subroutine export_fields (gcomp, rc) enddo end if #endif - ! TODO: needs to be at export freq - ! surface stokes drift - if (state_fldchk(exportState, 'Sw_ustokes')) then + if ( state_fldchk(exportState, 'Sw_ustokes') .and. & + state_fldchk(exportState, 'Sw_vstokes') )then + if (ChkErr(rc,__LINE__,u_FILE_u)) return call state_getfldptr(exportState, 'Sw_ustokes', sw_ustokes, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - sw_ustokes(:) = fillvalue - do jsea=1, nseal_cpl - call init_get_isea(isea, jsea) - ix = mapsf(isea,1) - iy = mapsf(isea,2) - if (mapsta(iy,ix) == 1) then - sw_ustokes(jsea) = USSX(jsea) - else - sw_ustokes(jsea) = 0. - endif - enddo - end if - if (state_fldchk(exportState, 'Sw_vstokes')) then call state_getfldptr(exportState, 'Sw_vstokes', sw_vstokes, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - sw_vstokes(:) = fillvalue - do jsea=1, nseal_cpl - call init_get_isea(isea, jsea) - ix = mapsf(isea,1) - iy = mapsf(isea,2) - if (mapsta(iy,ix) == 1) then - sw_vstokes(jsea) = USSY(jsea) - else - sw_vstokes(jsea) = 0. - endif - enddo + call CalcStokes(va, sw_ustokes, sw_vstokes, fillvalue) end if if (state_fldchk(exportState, 'Sw_ch')) then @@ -783,7 +760,7 @@ subroutine export_fields (gcomp, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call state_getfldptr(exportState, 'Sw_wavsvv', syyn, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call CalcRadstr2D( va, sxxn, sxyn, syyn) + call CalcRadstr2D( va, sxxn, sxyn, syyn, fillvalue) end if if (wav_coupling_to_cice) then @@ -821,14 +798,14 @@ subroutine export_fields (gcomp, rc) call state_getfldptr(exportState, 'Sw_hs', sw_hs, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return sw_hs(:) = fillvalue - call CalcHs(va, sw_hs) + call CalcHs(va, sw_hs, fillvalue) end if if (state_fldchk(exportState, 'Sw_bhd')) then call state_getfldptr(exportState, 'Sw_bhd', sw_bhd, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return sw_bhd(:) = fillvalue - call CalcBHD(va, sw_bhd) + call CalcBHD(va, sw_bhd, fillvalue) end if if ( state_fldchk(exportState, 'Sw_tauox') .and. & @@ -847,9 +824,6 @@ subroutine export_fields (gcomp, rc) if (mapsta(iy,ix) == 1) then sw_tauox(jsea) = tauox(jsea) sw_tauoy(jsea) = tauoy(jsea) - else - sw_tauox(jsea) = 0. - sw_tauoy(jsea) = 0. endif enddo end if @@ -862,7 +836,7 @@ subroutine export_fields (gcomp, rc) call state_getfldptr(exportState, 'Sw_ubay', sw_ubay, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! TODO: - !call CalcUVB(va, sw_ubax, sw_ubay) + !call CalcUVB(va, sw_ubax, sw_ubay, fillvalue) end if if (dbug_flag > 5) then @@ -1059,14 +1033,14 @@ subroutine CalcCharnk ( chkn ) #endif ! input/output variables - real(ESMF_KIND_R8), pointer :: chkn(:) ! 1D Charnock export field pointer + real(r8), pointer :: chkn(:) ! 1D Charnock export field pointer ! local variables - integer :: isea, jsea, ix, iy - real :: emean, fmean, fmean1, wnmean, amax, ustar, ustdr - real :: tauwx, tauwy, cd, z0, fmeanws, dlwmean - logical :: llws(nspec) - logical, save :: firstCall = .true. + integer :: isea, jsea, ix, iy + real :: emean, fmean, fmean1, wnmean, amax, ustar, ustdr + real :: tauwx, tauwy, cd, z0, fmeanws, dlwmean + logical :: llws(nspec) + logical, save :: firstCall = .true. !---------------------------------------------------------------------- !TODO: fix firstCall like for Roughl @@ -1256,33 +1230,34 @@ end subroutine CalcRoughl !! @param sxxn a 1-D pointer to a field on a mesh !! @param sxyn a 1-D pointer to a field on a mesh !! @param syyn a 1-D pointer to a field on a mesh - !! + !! @param[in] fval fill value !> @author Denise.Worthen@noaa.gov !> @date 08-05-2024 - subroutine CalcRadstr2D ( a, sxxn, sxyn, syyn ) + subroutine CalcRadstr2D ( a, sxxn, sxyn, syyn, fval) use w3gdatmd, only : nseal, nk, nth, sig, es2, esc, ec2, fte, dden, mapsf, mapsta use w3adatmd, only : dw, cg, wn use w3odatmd, only : naproc, iaproc ! input/output variables - real, intent(in) :: a(nth,nk,0:nseal) ! Input spectra (in par list to change shape) - real(ESMF_KIND_R8), pointer :: sxxn(:) ! eastward-component export field - real(ESMF_KIND_R8), pointer :: sxyn(:) ! eastward-northward-component export field - real(ESMF_KIND_R8), pointer :: syyn(:) ! northward-component export field + real, intent(in) :: a(nth,nk,0:nseal) ! Input spectra (in par list to change shape) + real(r8), intent(in) :: fval + real(r8), pointer, intent(inout) :: sxxn(:) ! eastward-component export field + real(r8), pointer, intent(inout) :: sxyn(:) ! eastward-northward-component export field + real(r8), pointer, intent(inout) :: syyn(:) ! northward-component export field ! local variables integer :: isea, jsea, ik, ith, ix, iy - real :: factor, abxx, abyy, abxy + real :: factor, abxx, abyy, abxy, sxx1, syy1, sxy1 - sxxn = 0.0 - syyn = 0.0 - sxyn = 0.0 do jsea = 1,nseal_cpl call init_get_isea(isea, jsea) ix = mapsf(isea,1) ! global ix iy = mapsf(isea,2) ! global iy if (mapsta(iy,ix) == 1) then ! active sea point + sxx1 = 0.0 + syy1 = 0.0 + sxy1 = 0.0 do ik = 1,nk factor = max ( 0.5, cg(ik,isea)/sig(ik)*wn(ik,isea) ) abxx = 0.0 @@ -1299,18 +1274,24 @@ subroutine CalcRadstr2D ( a, sxxn, sxyn, syyn ) abyy = max ( 0.0, abyy ) * factor abxy = abxy * factor - sxxn(jsea) = sxxn(jsea) + abxx - syyn(jsea) = syyn(jsea) + abyy - sxyn(jsea) = sxyn(jsea) + abxy - end do - sxxn(jsea) = sxxn(jsea) + fte * abxx/cg(nk,isea) - syyn(jsea) = syyn(jsea) + fte * abyy/cg(nk,isea) - sxyn(jsea) = sxyn(jsea) + fte * abxy/cg(nk,isea) + sxx1 = sxx1 + abxx + syy1 = syy1 + abyy + sxy1 = sxy1 + abxy + end do !ik + sxx1 = sxx1 + fte * abxx/cg(nk,isea) + syy1 = syy1 + fte * abyy/cg(nk,isea) + sxy1 = sxy1 + fte * abxy/cg(nk,isea) + end if + if (mapsta(iy,ix) == 1) then ! active sea point + sxxn(jsea) = sxx1*dwat*grav + syyn(jsea) = syy1*dwat*grav + sxyn(jsea) = sxy1*dwat*grav + else + sxxn(jsea) = fval + syyn(jsea) = fval + sxyn(jsea) = fval end if end do - sxxn = sxxn*dwat*grav - syyn = syyn*dwat*grav - sxyn = sxyn*dwat*grav end subroutine CalcRadstr2D @@ -1329,7 +1310,7 @@ subroutine CalcEF (a, wave_elevation_spectrum) use constants, only : tpi use w3gdatmd, only : nth, nk, nseal, mapsf, mapsta, dden, dsii - use w3adatmd, only : nsealm, cg + use w3adatmd, only : cg use w3parall, only : init_get_isea ! input/output variables @@ -1377,23 +1358,22 @@ end subroutine CalcEF !! !> @author Denise.Worthen@noaa.gov !> @date 8-02-2024 - subroutine CalcHs (a, hs) + subroutine CalcHs (a, hs, fval) use constants, only : tpi use w3gdatmd, only : nth, nk, nseal, mapsf, mapsta, dden, fte - use w3adatmd, only : nsealm, cg + use w3adatmd, only : cg use w3parall, only : init_get_isea ! input/output variables - real, intent(in) :: a(nth,nk,0:nseal) - real(r8), pointer :: hs(:) + real, intent(in) :: a(nth,nk,0:nseal) + real(r8), intent(in) :: fval + real(r8), pointer, intent(inout) :: hs(:) ! local variables real :: factor, eband, ab, et integer :: ik, ith, isea, jsea, ix, iy - et = 0.0 - ab = 0.0 do jsea = 1,nseal_cpl call init_get_isea(isea, jsea) ix = mapsf(isea,1) ! global ix @@ -1419,6 +1399,8 @@ subroutine CalcHs (a, hs) #else hs(jsea) = 4.0*sqrt ( et ) #endif + else + hs(jsea) = fval end if end do end subroutine CalcHs @@ -1430,34 +1412,33 @@ end subroutine CalcHs !! that exported BHD field is updated at the coupling frequency !! !! @param[in] a input spectra + !! @param[in] fval fillvalue !! @param[inout] bhd a 1-D pointer to a field on a mesh !! !> @author Denise.Worthen@noaa.gov !> @date 8-02-2024 - subroutine CalcBHD (a, bhd) + subroutine CalcBHD (a, bhd, fval) - use constants, only : tpi - use w3gdatmd, only : nth, nk, nseal, mapsf, mapsta, dden, fte - use w3adatmd, only : dw, nsealm, cg, wn + use w3gdatmd, only : nth, nk, nseal, mapsf, mapsta, dden + use w3adatmd, only : dw, cg, wn use w3parall, only : init_get_isea ! input/output variables - real, intent(in) :: a(nth,nk,0:nseal) - real(r8), pointer :: bhd(:) + real, intent(in) :: a(nth,nk,0:nseal) + real(r8), intent(in) :: fval + real(r8), pointer, intent(inout) :: bhd(:) ! local variables - real :: factor, kd, ab, ebd + real :: factor, kd, ab, ebd, bhd1 integer :: ik, ith, isea, jsea, ix, iy - ebd = 0.0 - ab = 0.0 - bhd = 0.0 do jsea = 1,nseal_cpl call init_get_isea(isea, jsea) ix = mapsf(isea,1) ! global ix iy = mapsf(isea,2) ! global iy if (mapsta(iy,ix) == 1) then ! active sea point ebd = 0.0 + bhd1 = 0.0 do ik = 1,nk factor = dden(ik) / cg(ik,isea) ab = 0.0 @@ -1467,14 +1448,81 @@ subroutine CalcBHD (a, bhd) ebd = ab*factor kd = max ( 0.001 , wn(ik,isea) * dw(isea) ) if (kd .lt. 6.0) then - bhd(jsea) = bhd(jsea) + grav*wn(ik,isea) * ebd / (sinh(2.*kd)) + bhd1 = bhd1 + grav*wn(ik,isea) * ebd / (sinh(2.*kd)) end if end do !ik + bhd(jsea) = bhd1 + else + bhd(jsea) = fval end if end do end subroutine CalcBHD + !==================================================================================== + !> Calculate Stokes drift for export + !! + !> @details Calculates Stokes drift independently of w3iogomd to ensure + !! that exported USSX and USSY fields are updated at the coupling frequency + !! + !! @param[in] a input spectra + !! @param[in] fval fill value + !! @param[inout] us a 1-D pointer to a field on a mesh + !! @param[inout] vs a 1-D pointer to a field on a mesh + !! + !> @author Denise.Worthen@noaa.gov + !> @date 8-02-2024 + subroutine CalcStokes(a, us, vs, fval) + + use w3gdatmd, only : nth, nk, nseal, mapsf, mapsta, dden, ecos, esin + use w3adatmd, only : dw, cg, wn + use w3gdatmd, only : sig + use w3parall, only : init_get_isea + + ! input/output variables + real, intent(in) :: a(nth,nk,0:nseal) + real(r8), intent(in) :: fval + real(r8), pointer, intent(inout) :: us(:), vs(:) + + ! local variables + real :: factor, kd, abx, aby, fkd, ussco, us1, vs1 + integer :: ik, ith, isea, jsea, ix, iy + + do jsea = 1,nseal_cpl + call init_get_isea(isea, jsea) + ix = mapsf(isea,1) ! global ix + iy = mapsf(isea,2) ! global iy + if (mapsta(iy,ix) == 1) then ! active sea point + us1 = 0.0 + vs1 = 0.0 + do ik = 1,nk + factor = dden(ik) / cg(ik,isea) + abx = 0.0 + aby = 0.0 + do ith = 1,nth + abx = abx + a(ith,ik,jsea)*ecos(ith) + aby = aby + a(ith,ik,jsea)*esin(ith) + end do + kd = max ( 0.001 , wn(ik,isea) * dw(isea) ) + if (kd .lt. 6.0) then + fkd = factor / sinh(kd)**2 + ussco = fkd*sig(ik)*wn(ik,isea)*cosh(2.0*kd) + else + ussco = factor*sig(ik)*2.0*wn(ik,isea) + end if + us1 = us1 + abx*ussco + vs1 = vs1 + aby*ussco + end do !ik + us(jsea) = us1 + vs(jsea) = vs1 + else + us(jsea) = fval + vs(jsea) = fval + end if + end do + + end subroutine CalcStokes + !==================================================================================== !> Create a global field across all PEs !! From e4be3138fd1e6af422cd0318d7169bd3e55ddf25 Mon Sep 17 00:00:00 2001 From: Denise Worthen Date: Tue, 6 Aug 2024 14:04:15 -0400 Subject: [PATCH 10/16] allow cesm to use history frequency of ussx,ussy * ensures consistency with calculation of lamult and lasl export fields --- model/src/wav_import_export.F90 | 39 +++++++++++++++++++++++++++++++-- 1 file changed, 37 insertions(+), 2 deletions(-) diff --git a/model/src/wav_import_export.F90 b/model/src/wav_import_export.F90 index e672afcbd..65975d611 100644 --- a/model/src/wav_import_export.F90 +++ b/model/src/wav_import_export.F90 @@ -668,7 +668,8 @@ subroutine export_fields (gcomp, rc) if (multigrid) then call wmsetm ( 1, mdse, mdst ) end if -#else +#endif +#ifdef W3_CESMCOUPLED if (state_fldchk(exportState, 'Sw_lamult')) then call state_getfldptr(exportState, 'Sw_lamult', sw_lamult, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -716,7 +717,40 @@ subroutine export_fields (gcomp, rc) endif enddo end if -#endif + + ! surface stokes drift at history frequency + if (state_fldchk(exportState, 'Sw_ustokes')) then + call state_getfldptr(exportState, 'Sw_ustokes', sw_ustokes, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + sw_ustokes(:) = fillvalue + do jsea=1, nseal_cpl + call init_get_isea(isea, jsea) + ix = mapsf(isea,1) + iy = mapsf(isea,2) + if (mapsta(iy,ix) == 1) then + sw_ustokes(jsea) = USSX(jsea) + else + sw_ustokes(jsea) = 0. + endif + enddo + end if + if (state_fldchk(exportState, 'Sw_vstokes')) then + call state_getfldptr(exportState, 'Sw_vstokes', sw_vstokes, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + sw_vstokes(:) = fillvalue + do jsea=1, nseal_cpl + call init_get_isea(isea, jsea) + ix = mapsf(isea,1) + iy = mapsf(isea,2) + if (mapsta(iy,ix) == 1) then + sw_vstokes(jsea) = USSY(jsea) + else + sw_vstokes(jsea) = 0. + endif + enddo + end if +#else + ! surface stokes drift at coupling frequency if ( state_fldchk(exportState, 'Sw_ustokes') .and. & state_fldchk(exportState, 'Sw_vstokes') )then if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -726,6 +760,7 @@ subroutine export_fields (gcomp, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call CalcStokes(va, sw_ustokes, sw_vstokes, fillvalue) end if +#endif if (state_fldchk(exportState, 'Sw_ch')) then call state_getfldptr(exportState, 'charno', charno, rc=rc) From 95010ab337ca9f50669ef2496d77a53073004a49 Mon Sep 17 00:00:00 2001 From: "denise.worthen" Date: Wed, 7 Aug 2024 10:32:41 -0500 Subject: [PATCH 11/16] add bottom currents and fix rad stress tag * add routine to calculate the ubrx,ubry * revert change to 'tag' for sxx,syy,sxy fields. the tag to get the stresses as output is sxy, which in w3iogonc then turns on the output of sxx,sxy,syy. this is like the tag 'cur' which is used to turn on the output for both cx and cy --- model/src/wav_grdout.F90 | 7 +- model/src/wav_import_export.F90 | 247 ++++++++++++++++---------------- 2 files changed, 126 insertions(+), 128 deletions(-) diff --git a/model/src/wav_grdout.F90 b/model/src/wav_grdout.F90 index 9d4211dc6..3c09913d9 100644 --- a/model/src/wav_grdout.F90 +++ b/model/src/wav_grdout.F90 @@ -130,7 +130,6 @@ subroutine initialize_gridout gridoutdefs(:,:)%dims = "" gridoutdefs(:,:)%validout = .false. - ! TODO: confirm unit values ! 1 Forcing Fields gridoutdefs(1,1:14) = [ & varatts( "DPT ", "DW ", "Water depth ", "m ", " ", .false.) , & @@ -178,7 +177,7 @@ subroutine initialize_gridout 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) + !TODO: has reverse indices (nk,nsea) varatts( "WN ", "WN ", "Wavenumber array ", "m-1 ", "k ", .false.) & ] @@ -223,8 +222,8 @@ subroutine initialize_gridout ! 6 Wave-ocean layer gridoutdefs(6,1:25) = [ & - varatts( "SXX ", "SXX ", "Radiation stresses xx ", "N m-1 ", " ", .false.) , & - varatts( "SYY ", "SYY ", "Radiation stresses yy ", "N m-1 ", " ", .false.) , & + 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.) , & diff --git a/model/src/wav_import_export.F90 b/model/src/wav_import_export.F90 index 65975d611..c9ee0a5e1 100644 --- a/model/src/wav_import_export.F90 +++ b/model/src/wav_import_export.F90 @@ -155,8 +155,11 @@ subroutine advertise_fields(importState, ExportState, flds_scalar_name, rc) call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_tauoy') call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_taubblx') call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_taubbly') - call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_ubax') - call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_ubay') + call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_ubrx') + call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_ubry') + call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_thm') + call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_t0m1') + call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_wnmean') end if call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_pstokes_x', ungridded_lbound=1, ungridded_ubound=3) call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_pstokes_y', ungridded_lbound=1, ungridded_ubound=3) @@ -606,6 +609,8 @@ subroutine export_fields (gcomp, rc) #else use wmmdatmd , only : mdse, mdst, wmsetm #endif + ! debug + use w3adatmd, only : thm, wnmean, t0m1 ! input/output/variables type(ESMF_GridComp) :: gcomp @@ -623,9 +628,6 @@ subroutine export_fields (gcomp, rc) real(r8), pointer :: z0rlen(:) real(r8), pointer :: charno(:) - real(r8), pointer :: wbcuru(:) - real(r8), pointer :: wbcurv(:) - real(r8), pointer :: wbcurp(:) real(r8), pointer :: sw_lamult(:) real(r8), pointer :: sw_lasl(:) real(r8), pointer :: sw_ustokes(:) @@ -640,8 +642,11 @@ subroutine export_fields (gcomp, rc) real(r8), pointer :: sw_tauoy(:) real(r8), pointer :: sw_taubblx(:) real(r8), pointer :: sw_taubbly(:) - real(r8), pointer :: sw_ubax(:) - real(r8), pointer :: sw_ubay(:) + real(r8), pointer :: sw_ubrx(:) + real(r8), pointer :: sw_ubry(:) + real(r8), pointer :: sw_thm(:) + real(r8), pointer :: sw_t0m1(:) + real(r8), pointer :: sw_wnmean(:) ! d2 is location, d1 is frequency - nwav_elev_spectrum frequencies will be used real(r8), pointer :: wave_elevation_spectrum(:,:) @@ -774,18 +779,6 @@ subroutine export_fields (gcomp, rc) call CalcRoughl(z0rlen) endif - ! if ( state_fldchk(exportState, 'wbcuru') .and. & - ! state_fldchk(exportState, 'wbcurv') .and. & - ! state_fldchk(exportState, 'wbcurp')) then - ! call state_getfldptr(exportState, 'wbcuru', wbcuru, rc=rc) - ! if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! call state_getfldptr(exportState, 'wbcurv', wbcurv, rc=rc) - ! if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! call state_getfldptr(exportState, 'wbcurp', wbcurp, rc=rc) - ! if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! call CalcBotcur( va, wbcuru, wbcurv, wbcurp) - ! end if - if ( state_fldchk(exportState, 'Sw_wavsuu') .and. & state_fldchk(exportState, 'Sw_wavsuv') .and. & state_fldchk(exportState, 'Sw_wavsvv')) then @@ -863,15 +856,32 @@ subroutine export_fields (gcomp, rc) enddo end if - if ( state_fldchk(exportState, 'Sw_ubax') .and. & - state_fldchk(exportState, 'Sw_ubay') )then + if ( state_fldchk(exportState, 'Sw_ubrx') .and. & + state_fldchk(exportState, 'Sw_ubry') )then + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call state_getfldptr(exportState, 'Sw_ubrx', sw_ubrx, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getfldptr(exportState, 'Sw_ubax', sw_ubax, rc=rc) + call state_getfldptr(exportState, 'Sw_ubry', sw_ubry, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getfldptr(exportState, 'Sw_ubay', sw_ubay, rc=rc) + call CalcUVBed(va, sw_ubrx, sw_ubry, fillvalue) + end if + + if (state_fldchk(exportState, 'Sw_thm')) then + call state_getfldptr(exportState, 'Sw_thm', sw_thm, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + sw_thm(:) = thm(:) + end if + + if (state_fldchk(exportState, 'Sw_t0m1')) then + call state_getfldptr(exportState, 'Sw_t0m1', sw_t0m1, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + sw_t0m1(:) = t0m1(:) + end if + + if (state_fldchk(exportState, 'Sw_wnmean')) then + call state_getfldptr(exportState, 'Sw_wnmean', sw_wnmean, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! TODO: - !call CalcUVB(va, sw_ubax, sw_ubay, fillvalue) + sw_wnmean(:) = wnmean(:) end if if (dbug_flag > 5) then @@ -1172,89 +1182,6 @@ subroutine CalcRoughl ( wrln) end subroutine CalcRoughl - ! !=============================================================================== - ! !> Calculate wave-bottom currents for export - ! !! - ! !> @details TODO: - ! !! - ! !! @param[in] a input spectra - ! !! @param wbxn a 1-D pointer to a field on a mesh - ! !! @param wbyn a 1-D pointer to a field on a mesh - ! !! @param wbpn a 1-D pointer to a field on a mesh - ! !! - ! !> @author T. J. Campbell, NRL - ! !> @date 09-Aug-2017 - ! subroutine CalcBotcur ( a, wbxn, wbyn, wbpn ) - - ! ! Calculate wave-bottom currents for export - - ! use w3gdatmd, only : nseal, nk, nth, sig, dmin, ecos, esin, dden, mapsf, mapsta, nspec - ! use w3adatmd, only : dw, cg, wn - ! use w3odatmd, only : naproc, iaproc - - ! ! input/output variables - ! real, intent(in) :: a(nth,nk,0:nseal) ! Input spectra (in par list to change shape) - ! real(ESMF_KIND_R8), pointer :: wbxn(:) ! eastward-component export field pointer - ! real(ESMF_KIND_R8), pointer :: wbyn(:) ! northward-component export field pointer - ! real(ESMF_KIND_R8), pointer :: wbpn(:) ! period export field pointer - - ! ! local variables - ! real(8), parameter :: half = 0.5_r8 - ! real(8), parameter :: one = 1.0_r8 - ! real(8), parameter :: two = 2.0_r8 - ! real(8), parameter :: kdmin = 1e-7_r8 - ! real(8), parameter :: kdmax = 18.0_r8 - ! integer :: isea, jsea, ik, ith - ! real(8) :: depth - ! real(8) :: kd, fack, fkd, aka, akx, aky, abr, ubr, ubx, uby, dir - ! real(8), allocatable :: sig2(:) - ! !---------------------------------------------------------------------- - - ! allocate( sig2(1:nk) ) - ! sig2(1:nk) = sig(1:nk)**2 - - ! wbxn(:) = zero - ! wbyn(:) = zero - ! wbpn(:) = zero - - ! jsea_loop: do jsea = 1,nseal_cpl - ! call init_get_isea(isea, jsea) - ! if ( dw(isea).le.zero ) cycle jsea_loop - ! depth = max(dmin,dw(isea)) - ! abr = zero - ! ubr = zero - ! ubx = zero - ! uby = zero - ! ik_loop: do ik = 1,nk - ! aka = zero - ! akx = zero - ! aky = zero - ! ith_loop: do ith = 1,nth - ! aka = aka + a(ith,ik,jsea) - ! akx = akx + a(ith,ik,jsea)*ecos(ith) - ! aky = aky + a(ith,ik,jsea)*esin(ith) - ! enddo ith_loop - ! fack = dden(ik)/cg(ik,isea) - ! kd = max(kdmin,min(kdmax,wn(ik,isea)*depth)) - ! fkd = fack/sinh(kd)**2 - ! abr = abr + aka*fkd - ! ubr = ubr + aka*sig2(ik)*fkd - ! ubx = ubx + akx*sig2(ik)*fkd - ! uby = uby + aky*sig2(ik)*fkd - ! enddo ik_loop - ! if ( abr.le.zero .or. ubr.le.zero ) cycle jsea_loop - ! abr = sqrt(two*abr) - ! ubr = sqrt(two*ubr) - ! dir = atan2(uby,ubx) - ! wbxn(jsea) = ubr*cos(dir) - ! wbyn(jsea) = ubr*sin(dir) - ! wbpn(jsea) = tpi*abr/ubr - ! enddo jsea_loop - - ! deallocate( sig2 ) - - ! end subroutine CalcBotcur - !=============================================================================== !> Calculate radiation stresses for export !! @@ -1275,11 +1202,11 @@ subroutine CalcRadstr2D ( a, sxxn, sxyn, syyn, fval) use w3odatmd, only : naproc, iaproc ! input/output variables - real, intent(in) :: a(nth,nk,0:nseal) ! Input spectra (in par list to change shape) - real(r8), intent(in) :: fval - real(r8), pointer, intent(inout) :: sxxn(:) ! eastward-component export field - real(r8), pointer, intent(inout) :: sxyn(:) ! eastward-northward-component export field - real(r8), pointer, intent(inout) :: syyn(:) ! northward-component export field + real, intent(in) :: a(nth,nk,0:nseal) ! Input spectra (in par list to change shape) + real(ESMF_KIND_R8), intent(in) :: fval + real(ESMF_KIND_R8), pointer, intent(inout) :: sxxn(:) ! eastward-component export field + real(ESMF_KIND_R8), pointer, intent(inout) :: sxyn(:) ! eastward-northward-component export field + real(ESMF_KIND_R8), pointer, intent(inout) :: syyn(:) ! northward-component export field ! local variables integer :: isea, jsea, ik, ith, ix, iy @@ -1349,8 +1276,8 @@ subroutine CalcEF (a, wave_elevation_spectrum) use w3parall, only : init_get_isea ! input/output variables - real, intent(in) :: a(nth,nk,0:nseal) - real(r8), pointer :: wave_elevation_spectrum(:,:) + real, intent(in) :: a(nth,nk,0:nseal) + real(ESMF_KIND_R8), pointer :: wave_elevation_spectrum(:,:) ! local variables real :: ab(nseal) @@ -1401,9 +1328,9 @@ subroutine CalcHs (a, hs, fval) use w3parall, only : init_get_isea ! input/output variables - real, intent(in) :: a(nth,nk,0:nseal) - real(r8), intent(in) :: fval - real(r8), pointer, intent(inout) :: hs(:) + real, intent(in) :: a(nth,nk,0:nseal) + real(ESMF_KIND_R8), intent(in) :: fval + real(ESMF_KIND_R8), pointer, intent(inout) :: hs(:) ! local variables real :: factor, eband, ab, et @@ -1459,9 +1386,9 @@ subroutine CalcBHD (a, bhd, fval) use w3parall, only : init_get_isea ! input/output variables - real, intent(in) :: a(nth,nk,0:nseal) - real(r8), intent(in) :: fval - real(r8), pointer, intent(inout) :: bhd(:) + real, intent(in) :: a(nth,nk,0:nseal) + real(ESMF_KIND_R8), intent(in) :: fval + real(ESMF_KIND_R8), pointer, intent(inout) :: bhd(:) ! local variables real :: factor, kd, ab, ebd, bhd1 @@ -1515,9 +1442,9 @@ subroutine CalcStokes(a, us, vs, fval) use w3parall, only : init_get_isea ! input/output variables - real, intent(in) :: a(nth,nk,0:nseal) - real(r8), intent(in) :: fval - real(r8), pointer, intent(inout) :: us(:), vs(:) + real, intent(in) :: a(nth,nk,0:nseal) + real(ESMF_KIND_R8), intent(in) :: fval + real(ESMF_KIND_R8), pointer, intent(inout) :: us(:), vs(:) ! local variables real :: factor, kd, abx, aby, fkd, ussco, us1, vs1 @@ -1558,6 +1485,78 @@ subroutine CalcStokes(a, us, vs, fval) end subroutine CalcStokes + !==================================================================================== + !> Calculate UVBed drift for export + !! + !> @details Calculates near bed orbital velocities independently of w3iogomd to + !! ensure that exported UBRX and UBRY fields are updated at the coupling frequency + !! + !! @param[in] a input spectra + !! @param[in] fval fill value + !! @param[inout] ubrx a 1-D pointer to a field on a mesh + !! @param[inout] vbry a 1-D pointer to a field on a mesh + !! + !> @author Denise.Worthen@noaa.gov + !> @date 8-02-2024 + subroutine CalcUVBed(a, ubrx, ubry, fval) + + use w3gdatmd, only : nth, nk, nseal, mapsf, mapsta, dden, ecos, esin + use w3adatmd, only : dw, cg, wn + use w3gdatmd, only : sig + use w3parall, only : init_get_isea + + ! input/output variables + real, intent(in) :: a(nth,nk,0:nseal) + real(ESMF_KIND_R8), intent(in) :: fval + real(ESMF_KIND_R8), pointer, intent(inout) :: ubrx(:), ubry(:) + + ! local variables + real :: factor, kd, ab, abx, aby, fkd, ussco, uba1, ubd1, ubr1 + integer :: ik, ith, isea, jsea, ix, iy + + do jsea = 1,nseal_cpl + call init_get_isea(isea, jsea) + ix = mapsf(isea,1) ! global ix + iy = mapsf(isea,2) ! global iy + if (mapsta(iy,ix) == 1) then ! active sea point + uba1 = 0.0 + ubd1 = 0.0 + ubr1 = 0.0 + do ik = 1,nk + factor = dden(ik) / cg(ik,isea) + ab = 0.0 + abx = 0.0 + aby = 0.0 + do ith = 1,nth + ab = ab + a(ith,ik,jsea) + abx = abx + a(ith,ik,jsea)*ecos(ith) + aby = aby + a(ith,ik,jsea)*esin(ith) + end do + kd = max ( 0.001 , wn(ik,isea) * dw(isea) ) + if (kd .lt. 6.0) then + fkd = factor / sinh(kd)**2 + ubr1 = ubr1 + ab*sig(ik)**2 * fkd + uba1 = uba1 + abx*sig(ik)**2 * fkd + ubd1 = ubd1 + aby*sig(ik)**2 * fkd + end if + end do !ik + ubr1 = sqrt(2.0*max(0.0,ubr1)) + if (ubr1 .ge. 1.0e-7) then + ubd1 = atan2(ubd1,uba1) + else + ubd1 = 0.0 + end if + uba1 = ubr1 + ubrx(jsea) = uba1*cos(ubd1) + ubry(jsea) = uba1*sin(ubd1) + else + ubrx(jsea) = fval + ubry(jsea) = fval + end if + end do + + end subroutine CalcUVBed + !==================================================================================== !> Create a global field across all PEs !! From b9c70b2811b9d58d9d338002bba2a72c27222912 Mon Sep 17 00:00:00 2001 From: "denise.worthen" Date: Thu, 8 Aug 2024 07:30:53 -0500 Subject: [PATCH 12/16] add t0m1 and thm --- model/src/w3iogoncdmd.F90 | 2 +- model/src/wav_import_export.F90 | 135 +++++++++++++++++++++++++++++--- 2 files changed, 127 insertions(+), 10 deletions(-) diff --git a/model/src/w3iogoncdmd.F90 b/model/src/w3iogoncdmd.F90 index 8481583c9..3cfd015ff 100644 --- a/model/src/w3iogoncdmd.F90 +++ b/model/src/w3iogoncdmd.F90 @@ -455,7 +455,7 @@ subroutine write_var2d(vname, var, dir, usemask, init0, init2, fldir) end if if (lfldir) then - if (varloc.ne.undef) then + if (varloc .ne. undef) then varloc = mod(630. - rade*varloc, 360.) end if end if diff --git a/model/src/wav_import_export.F90 b/model/src/wav_import_export.F90 index c9ee0a5e1..c91449704 100644 --- a/model/src/wav_import_export.F90 +++ b/model/src/wav_import_export.F90 @@ -595,7 +595,7 @@ subroutine export_fields (gcomp, rc) !--------------------------------------------------------------------------- use wav_kind_mod, only : R8 => SHR_KIND_R8 - use w3adatmd , only : USSX, USSY, USSP, HS, tauox, tauoy + use w3adatmd , only : USSX, USSY, USSP, HS, tauox, tauoy, wnmean use w3adatmd , only : w3seta use w3idatmd , only : w3seti use w3wdatmd , only : va, w3setw @@ -609,8 +609,6 @@ subroutine export_fields (gcomp, rc) #else use wmmdatmd , only : mdse, mdst, wmsetm #endif - ! debug - use w3adatmd, only : thm, wnmean, t0m1 ! input/output/variables type(ESMF_GridComp) :: gcomp @@ -826,7 +824,7 @@ subroutine export_fields (gcomp, rc) call state_getfldptr(exportState, 'Sw_hs', sw_hs, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return sw_hs(:) = fillvalue - call CalcHs(va, sw_hs, fillvalue) + call CalcHS(va, sw_hs, fillvalue) end if if (state_fldchk(exportState, 'Sw_bhd')) then @@ -869,13 +867,13 @@ subroutine export_fields (gcomp, rc) if (state_fldchk(exportState, 'Sw_thm')) then call state_getfldptr(exportState, 'Sw_thm', sw_thm, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - sw_thm(:) = thm(:) + call CalcTHM(va, sw_thm, fillvalue) end if if (state_fldchk(exportState, 'Sw_t0m1')) then call state_getfldptr(exportState, 'Sw_t0m1', sw_t0m1, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - sw_t0m1(:) = t0m1(:) + call CalcT0M1(va, sw_t0m1, fillvalue) end if if (state_fldchk(exportState, 'Sw_wnmean')) then @@ -1126,7 +1124,6 @@ end subroutine CalcCharnk subroutine CalcRoughl ( wrln) ! Calculate wave roughness length for export - use w3gdatmd, only : nseal, nk, nth, sig, dmin, ecos, esin, dden, mapsf, mapsta, nspec use w3adatmd, only : dw, cg, wn, charn, u10, u10d use w3wdatmd, only : va, ust @@ -1320,7 +1317,7 @@ end subroutine CalcEF !! !> @author Denise.Worthen@noaa.gov !> @date 8-02-2024 - subroutine CalcHs (a, hs, fval) + subroutine CalcHS (a, hs, fval) use constants, only : tpi use w3gdatmd, only : nth, nk, nseal, mapsf, mapsta, dden, fte @@ -1365,7 +1362,7 @@ subroutine CalcHs (a, hs, fval) hs(jsea) = fval end if end do - end subroutine CalcHs + end subroutine CalcHS !=============================================================================== !> Calculate Bernoulli head pressure for export @@ -1557,6 +1554,126 @@ subroutine CalcUVBed(a, ubrx, ubry, fval) end subroutine CalcUVBed + !=============================================================================== + !> Calculate mean wave direction for export + !! + !> @details Calculates mean wave direction independently of w3iogomd to ensure + !! that exported THM field is updated at the coupling frequency + !! + !! @param[in] a input spectra + !! @param[inout] thm a 1-D pointer to a field on a mesh + !! + !> @author Denise.Worthen@noaa.gov + !> @date 8-02-2024 + subroutine CalcTHM (a, thm, fval) + + use constants, only : rade + use w3gdatmd, only : nth, nk, nseal, mapsf, mapsta, dden, fte, ecos, esin + use w3adatmd, only : cg + use w3parall, only : init_get_isea + + ! input/output variables + real, intent(in) :: a(nth,nk,0:nseal) + real(ESMF_KIND_R8), intent(in) :: fval + real(ESMF_KIND_R8), pointer, intent(inout) :: thm(:) + + ! local variables + real :: factor, abx, aby, etx, ety + integer :: ik, ith, isea, jsea, ix, iy + + do jsea = 1,nseal_cpl + call init_get_isea(isea, jsea) + ix = mapsf(isea,1) ! global ix + iy = mapsf(isea,2) ! global iy + if (mapsta(iy,ix) == 1) then ! active sea point + etx = 0.0 + ety = 0.0 + do ik = 1,nk + factor = dden(ik) / cg(ik,isea) + abx = 0.0 + aby = 0.0 + do ith = 1,nth + abx = abx + a(ith,ik,jsea)*ecos(ith) + aby = aby + a(ith,ik,jsea)*esin(ith) + end do + etx = etx + abx*factor + ety = ety + aby*factor + end do !ik + etx = etx + fte * abx/cg(nk,isea) + ety = ety + fte * aby/cg(nk,isea) + if ( abs(etx) + abs(ety) .gt. 1.e-7 ) then + thm(jsea) = atan2(ety,etx) + else + thm(jsea) = 0.0 + end if + ! convert to degrees + thm(jsea) = mod(630.0 - rade*thm(jsea), 360.0) + else + thm(jsea) = fval + end if + end do + + end subroutine CalcTHM + + !=============================================================================== + !> Calculate mean wave direction for export + !! + !> @details Calculates mean wave period independently of w3iogomd to ensure + !! that exported T0M1 field is updated at the coupling frequency + !! + !! @param[in] a input spectra + !! @param[inout] thm a 1-D pointer to a field on a mesh + !! + !> @author Denise.Worthen@noaa.gov + !> @date 8-02-2024 + subroutine CalcT0M1 (a, t0m1, fval) + + use constants, only : tpi + use w3gdatmd, only : nth, nk, nseal, mapsf, mapsta, dden, fte, fttr, sig + use w3adatmd, only : cg + use w3parall, only : init_get_isea + + ! input/output variables + real, intent(in) :: a(nth,nk,0:nseal) + real(ESMF_KIND_R8), intent(in) :: fval + real(ESMF_KIND_R8), pointer, intent(inout) :: t0m1(:) + + ! local variables + real :: factor, eband, ab, et, ebd, etr + integer :: ik, ith, isea, jsea, ix, iy + + do jsea = 1,nseal_cpl + call init_get_isea(isea, jsea) + ix = mapsf(isea,1) ! global ix + iy = mapsf(isea,2) ! global iy + if (mapsta(iy,ix) == 1) then ! active sea point + etr = 0.0 + et = 0.0 + do ik = 1,nk + factor = dden(ik) / cg(ik,isea) + ab = 0.0 + do ith = 1,nth + ab = ab + a(ith,ik,jsea) + end do + ebd = ab*factor + et = et + ebd + etr = etr + ebd/sig(ik) + end do !ik + eband = ab/cg(nk,isea) + et = et + fte*eband + etr = etr + fttr*eband + if (et .gt. 1.0e-7) then + t0m1(jsea) = etr/et * tpi + else + t0m1(jsea) = tpi/sig(nk) + end if + else + t0m1(jsea) = fval + end if + end do + + end subroutine CalcT0M1 + !==================================================================================== !> Create a global field across all PEs !! From bd72ace6dbaab1f516e7896d133db976ebc597e2 Mon Sep 17 00:00:00 2001 From: "denise.worthen" Date: Thu, 8 Aug 2024 09:32:21 -0500 Subject: [PATCH 13/16] add bt4 switch and taubbl exports --- model/bin/switch_meshcap_pdlib_bt4 | 1 + model/src/wav_import_export.F90 | 13 ++++++++++++- 2 files changed, 13 insertions(+), 1 deletion(-) create mode 100644 model/bin/switch_meshcap_pdlib_bt4 diff --git a/model/bin/switch_meshcap_pdlib_bt4 b/model/bin/switch_meshcap_pdlib_bt4 new file mode 100644 index 000000000..136c274e5 --- /dev/null +++ b/model/bin/switch_meshcap_pdlib_bt4 @@ -0,0 +1 @@ +NCO PDLIB SCOTCH NOGRB DIST MPI PR3 UQ FLX0 SEED ST4 STAB0 NL1 BT4 DB1 MLIM FLD2 TR0 BS0 RWND WNX1 WNT1 CRX1 CRT1 O0 O1 O2 O3 O4 O5 O6 O7 O14 O15 IC0 IS0 REF0 diff --git a/model/src/wav_import_export.F90 b/model/src/wav_import_export.F90 index c91449704..fee6cef86 100644 --- a/model/src/wav_import_export.F90 +++ b/model/src/wav_import_export.F90 @@ -595,7 +595,7 @@ subroutine export_fields (gcomp, rc) !--------------------------------------------------------------------------- use wav_kind_mod, only : R8 => SHR_KIND_R8 - use w3adatmd , only : USSX, USSY, USSP, HS, tauox, tauoy, wnmean + use w3adatmd , only : USSX, USSY, USSP, HS, tauox, tauoy, wnmean, taubbl use w3adatmd , only : w3seta use w3idatmd , only : w3seti use w3wdatmd , only : va, w3setw @@ -882,6 +882,17 @@ subroutine export_fields (gcomp, rc) sw_wnmean(:) = wnmean(:) end if + if ( state_fldchk(exportState, 'Sw_taubblx') .and. & + state_fldchk(exportState, 'Sw_taubbly') )then + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call state_getfldptr(exportState, 'Sw_taubblx', sw_taubblx, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call state_getfldptr(exportState, 'Sw_taubbly', sw_taubbly, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + sw_taubblx(:) = taubbl(:,1) + sw_taubbly(:) = taubbl(:,2) + end if + if (dbug_flag > 5) then call state_diagnose(exportState, 'at export ', rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return From b8d0ff7f2dd2682dc33bd6c32d3f4239e0a1d51b Mon Sep 17 00:00:00 2001 From: "denise.worthen" Date: Thu, 8 Aug 2024 12:29:48 -0500 Subject: [PATCH 14/16] fix debug failure --- model/src/wav_import_export.F90 | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/model/src/wav_import_export.F90 b/model/src/wav_import_export.F90 index fee6cef86..0515f0320 100644 --- a/model/src/wav_import_export.F90 +++ b/model/src/wav_import_export.F90 @@ -879,7 +879,15 @@ subroutine export_fields (gcomp, rc) if (state_fldchk(exportState, 'Sw_wnmean')) then call state_getfldptr(exportState, 'Sw_wnmean', sw_wnmean, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - sw_wnmean(:) = wnmean(:) + sw_wnmean(:) = fillvalue + do jsea=1, nseal_cpl + call init_get_isea(isea, jsea) + ix = mapsf(isea,1) + iy = mapsf(isea,2) + if (mapsta(iy,ix) == 1) then + sw_wnmean(jsea) = wnmean(jsea) + endif + enddo end if if ( state_fldchk(exportState, 'Sw_taubblx') .and. & @@ -889,8 +897,17 @@ subroutine export_fields (gcomp, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call state_getfldptr(exportState, 'Sw_taubbly', sw_taubbly, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - sw_taubblx(:) = taubbl(:,1) - sw_taubbly(:) = taubbl(:,2) + sw_taubblx(:) = fillvalue + sw_taubbly(:) = fillvalue + do jsea=1, nseal_cpl + call init_get_isea(isea, jsea) + ix = mapsf(isea,1) + iy = mapsf(isea,2) + if (mapsta(iy,ix) == 1) then + sw_taubblx(jsea) = taubbl(jsea,1) + sw_taubbly(jsea) = taubbl(jsea,2) + endif + enddo end if if (dbug_flag > 5) then From 6b229be6dfde47f88e7a1290bb415a2eb579cb71 Mon Sep 17 00:00:00 2001 From: "denise.worthen" Date: Thu, 8 Aug 2024 13:04:49 -0500 Subject: [PATCH 15/16] uncomment So_h --- model/src/wav_import_export.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/model/src/wav_import_export.F90 b/model/src/wav_import_export.F90 index 0515f0320..8f11862ff 100644 --- a/model/src/wav_import_export.F90 +++ b/model/src/wav_import_export.F90 @@ -105,7 +105,7 @@ subroutine advertise_fields(importState, ExportState, flds_scalar_name, rc) ! Advertise import fields !-------------------------------- - !call fldlist_add(fldsToWav_num, fldsToWav, 'So_h' ) + call fldlist_add(fldsToWav_num, fldsToWav, 'So_h' ) call fldlist_add(fldsToWav_num, fldsToWav, 'Si_ifrac' ) call fldlist_add(fldsToWav_num, fldsToWav, 'So_u' ) call fldlist_add(fldsToWav_num, fldsToWav, 'So_v' ) From 19f7a276cb8935ec26a01a737b03d874cc05be81 Mon Sep 17 00:00:00 2001 From: "denise.worthen" Date: Thu, 8 Aug 2024 14:42:17 -0500 Subject: [PATCH 16/16] remove unused variables from wav_import_export --- model/src/wav_import_export.F90 | 40 ++++++++++++++------------------- 1 file changed, 17 insertions(+), 23 deletions(-) diff --git a/model/src/wav_import_export.F90 b/model/src/wav_import_export.F90 index 8f11862ff..6cb9ebee4 100644 --- a/model/src/wav_import_export.F90 +++ b/model/src/wav_import_export.F90 @@ -274,13 +274,13 @@ subroutine import_fields( gcomp, time0, timen, rc ) ! Obtain the wave input from the mediator !--------------------------------------------------------------------------- - use w3gdatmd , only: nsea, nseal, MAPSTA, NX, NY, w3setg + use w3gdatmd , only: nsea, NX, NY, w3setg use w3idatmd , only: CX0, CY0, CXN, CYN, DT0, DTN, ICEI, WLEV, INFLAGS1, ICEP1, ICEP5 use w3idatmd , only: TC0, TCN, TLN, TIN, TI1, TI5, TW0, TWN, WX0, WY0, WXN, WYN use w3idatmd , only: UX0, UY0, UXN, UYN, TU0, TUN use w3idatmd , only: tfn, w3seti use w3odatmd , only: w3seto - use w3wdatmd , only: time, w3setw + use w3wdatmd , only: w3setw #ifdef W3_CESMCOUPLED use w3idatmd , only: HSL #else @@ -310,7 +310,6 @@ subroutine import_fields( gcomp, time0, timen, rc ) integer :: mpi_comm_null = -1 real(r4), allocatable :: wxdata(:) ! only needed if merge_import real(r4), allocatable :: wydata(:) ! only needed if merge_import - character(len=CL) :: msgString character(len=*), parameter :: subname='(wav_import_export:import_fields)' !--------------------------------------------------------------------------- @@ -599,8 +598,8 @@ subroutine export_fields (gcomp, rc) use w3adatmd , only : w3seta use w3idatmd , only : w3seti use w3wdatmd , only : va, w3setw - use w3odatmd , only : w3seto, naproc, iaproc - use w3gdatmd , only : nseal, mapsf, MAPSTA, USSPF, NK, w3setg + use w3odatmd , only : w3seto + use w3gdatmd , only : mapsf, MAPSTA, USSPF, NK, w3setg use w3iogomd , only : CALC_U3STOKES #ifdef W3_CESMCOUPLED use w3wdatmd , only : ASF, UST @@ -622,7 +621,7 @@ subroutine export_fields (gcomp, rc) real(R8) :: fillvalue = zero ! special missing value #endif type(ESMF_State) :: exportState - integer :: n, jsea, isea, ix, iy, ib + integer :: jsea, isea, ix, iy, ib real(r8), pointer :: z0rlen(:) real(r8), pointer :: charno(:) @@ -1092,22 +1091,20 @@ subroutine CalcCharnk ( chkn ) ! Calculate Charnok for export - use w3gdatmd, only : nseal, nk, nth, sig, mapsf, mapsta, nspec + use w3gdatmd, only : nk, nspec use w3adatmd, only : cg, wn, charn, u10, u10d use w3wdatmd, only : va - use w3odatmd, only : naproc, iaproc #ifdef W3_ST3 use w3src3md, only : w3spr3 #endif #ifdef W3_ST4 use w3src4md, only : w3spr4 #endif - ! input/output variables real(r8), pointer :: chkn(:) ! 1D Charnock export field pointer ! local variables - integer :: isea, jsea, ix, iy + integer :: isea, jsea real :: emean, fmean, fmean1, wnmean, amax, ustar, ustdr real :: tauwx, tauwy, cd, z0, fmeanws, dlwmean logical :: llws(nspec) @@ -1152,10 +1149,10 @@ end subroutine CalcCharnk subroutine CalcRoughl ( wrln) ! Calculate wave roughness length for export - use w3gdatmd, only : nseal, nk, nth, sig, dmin, ecos, esin, dden, mapsf, mapsta, nspec - use w3adatmd, only : dw, cg, wn, charn, u10, u10d + use w3gdatmd, only : nk, mapsf, mapsta, nspec + use w3adatmd, only : cg, wn, charn, u10, u10d use w3wdatmd, only : va, ust - use w3odatmd, only : naproc, iaproc, runtype + use w3odatmd, only : runtype #ifdef W3_ST3 use w3src3md, only : w3spr3 #endif @@ -1223,8 +1220,7 @@ end subroutine CalcRoughl subroutine CalcRadstr2D ( a, sxxn, sxyn, syyn, fval) use w3gdatmd, only : nseal, nk, nth, sig, es2, esc, ec2, fte, dden, mapsf, mapsta - use w3adatmd, only : dw, cg, wn - use w3odatmd, only : naproc, iaproc + use w3adatmd, only : cg, wn ! input/output variables real, intent(in) :: a(nth,nk,0:nseal) ! Input spectra (in par list to change shape) @@ -1536,7 +1532,7 @@ subroutine CalcUVBed(a, ubrx, ubry, fval) real(ESMF_KIND_R8), pointer, intent(inout) :: ubrx(:), ubry(:) ! local variables - real :: factor, kd, ab, abx, aby, fkd, ussco, uba1, ubd1, ubr1 + real :: factor, kd, ab, abx, aby, fkd, uba1, ubd1, ubr1 integer :: ik, ith, isea, jsea, ix, iy do jsea = 1,nseal_cpl @@ -1718,8 +1714,7 @@ end subroutine CalcT0M1 !> @date 01-05-2022 subroutine SetGlobalInput(importState, fldname, vm, global_output, rc) - use w3gdatmd, only: nsea, nx, ny - use w3odatmd, only: naproc, iaproc + use w3gdatmd, only: nsea ! input/output variables type(ESMF_State) , intent(in) :: importState @@ -1729,7 +1724,7 @@ subroutine SetGlobalInput(importState, fldname, vm, global_output, rc) integer , intent(out) :: rc ! local variables - integer :: jsea, isea, ix, iy + integer :: jsea, isea real(r4) :: global_input(nsea) real(r8), pointer :: dataptr(:) character(len=*), parameter :: subname = '(wav_import_export:setGlobalInput)' @@ -1837,7 +1832,7 @@ end subroutine fillglobal_with_merge_import !> @date 01-05-2022 subroutine set_importmask(importState, clock, fldname, vm, rc) - use w3gdatmd, only: nsea, nseal, nx, ny + use w3gdatmd, only: nsea, nseal use w3odatmd, only: naproc, iaproc ! input/output variables @@ -1852,7 +1847,7 @@ subroutine set_importmask(importState, clock, fldname, vm, rc) type(ESMF_TimeInterval) :: timeStep logical :: firstCall, secondCall real(r4) :: fillValue = 9.99e20 - integer :: isea, jsea, ix, iy + integer :: isea, jsea real(r8), pointer :: dataptr(:) real(r4) :: mask_local(nsea) character(len=CL) :: msgString @@ -1945,7 +1940,7 @@ subroutine check_globaldata(gcomp, fldname, global_data, nvals, rc) ! local variables type(ESMF_Clock) :: clock type(ESMF_State) :: importState - type(ESMF_Time) :: currtime, nexttime + type(ESMF_Time) :: nexttime type(ESMF_Field) :: lfield type(ESMF_Field) :: newfield type(ESMF_MeshLoc) :: meshloc @@ -1953,7 +1948,6 @@ subroutine check_globaldata(gcomp, fldname, global_data, nvals, rc) character(len=CS) :: timestr character(ESMF_MAXSTR) ,pointer :: lfieldnamelist(:) integer :: fieldCount - integer :: lrank integer :: yr,mon,day,sec ! time units integer :: jsea, isea, ix, iy real(r8), pointer :: dataptr1d(:)