Skip to content
New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Changes to enable coupling with WW3 using 3d vortex formulation #32

Merged
merged 2 commits into from
Jan 29, 2025
Merged

Changes to enable coupling with WW3 using 3d vortex formulation #32

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
da59350
adjust cap to allow 3d coupling using vortex formulation
uturuncoglu Aug 14, 2024
3f93003
update field names
uturuncoglu Sep 10, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
247 changes: 247 additions & 0 deletions src/schism/schism_esmf_util.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,7 @@ module schism_esmf_util
public type_InternalState, type_InternalStateWrapper
public SCHISM_StateFieldCreateRealize,SCHISM_StateFieldCreate
public SCHISM_StateImportWaveTensor
public SCHISM_StateImportWave3dVortex
public SCHISM_MeshCreateNode
public SCHISM_MeshCreateElement
public SCHISM_InitializePtrMap
Expand Down Expand Up @@ -2633,4 +2634,250 @@ subroutine SCHISM_StateImportWaveTensor(state, isPtr, rc)

end subroutine SCHISM_StateImportWaveTensor

#undef ESMF_METHOD
#define ESMF_METHOD "SCHISM_StateImportWave3dVortex"
subroutine SCHISM_StateImportWave3dVortex(state, isPtr, rc)

use schism_glbl, only: np
implicit none

type(ESMF_State), intent(in) :: state
type(type_InternalState), pointer, intent(in) :: isPtr
integer(ESMF_KIND_I4), intent(out), optional :: rc

type(ESMF_Field) :: field
type(ESMF_StateItem_Flag) :: itemType
character(len=ESMF_MAXSTR) :: message
integer(ESMF_KIND_I4) :: localrc

real(ESMF_KIND_R8), pointer :: farrayPtr1(:) => null()
real(ESMF_KIND_R8), allocatable :: Sw_hs(:)
real(ESMF_KIND_R8), allocatable :: Sw_bhd(:)
real(ESMF_KIND_R8), allocatable :: Sw_tauox(:)
real(ESMF_KIND_R8), allocatable :: Sw_tauoy(:)
real(ESMF_KIND_R8), allocatable :: Sw_taubblx(:)
real(ESMF_KIND_R8), allocatable :: Sw_taubbly(:)
real(ESMF_KIND_R8), allocatable :: Sw_ubrx(:)
real(ESMF_KIND_R8), allocatable :: Sw_ubry(:)
real(ESMF_KIND_R8), allocatable :: Sw_thm(:)
real(ESMF_KIND_R8), allocatable :: Sw_t0m1(:)
real(ESMF_KIND_R8), allocatable :: Sw_wnmean(:)
real(ESMF_KIND_R8), allocatable :: Sw_ustokes(:)
real(ESMF_KIND_R8), allocatable :: Sw_vstokes(:)

if (present(rc)) rc=ESMF_SUCCESS
localrc = ESMF_SUCCESS

allocate(farrayPtr1(np), stat=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call ESMF_StateGet(state, itemname='sea_surface_wave_significant_height', itemType=itemType, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

if (itemType /= ESMF_STATEITEM_FIELD) then
write(message,'(A)') '--- skipped coupling with wave through vortex formulation'
call ESMF_LogWrite(trim(message), ESMF_LOGMSG_INFO)
return
endif

call ESMF_StateGet(state, itemname='sea_surface_wave_significant_height', field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_FieldPtrUpdate(field, farrayPtr1, isPtr=isPtr, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

allocate(Sw_hs(np), stat=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
Sw_hs(:) = 0.0d0

write(message,'(A,2g14.7)') 'sea_surface_wave_significant_height = ', &
minval(farrayPtr1), maxval(farrayPtr1)
call ESMF_LogWrite(trim(message), ESMF_LOGMSG_INFO)
Sw_hs(:) = farrayPtr1(1:np)

call ESMF_StateGet(state, itemname='sea_water_waves_effect_on_currents_bernoulli_head_adjustment', field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_FieldPtrUpdate(field, farrayPtr1, isPtr=isPtr, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

allocate(Sw_bhd(np), stat=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
Sw_bhd(:) = 0.0d0

write(message,'(A,2g14.7)') 'sea_water_waves_effect_on_currents_bernoulli_head_adjustment = ', &
minval(farrayPtr1), maxval(farrayPtr1)
call ESMF_LogWrite(trim(message), ESMF_LOGMSG_INFO)
Sw_bhd(:) = farrayPtr1(1:np)

call ESMF_StateGet(state, itemname='sea_surface_x_stress_due_to_waves', field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_FieldPtrUpdate(field, farrayPtr1, isPtr=isPtr, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

allocate(Sw_tauox(np), stat=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
Sw_tauox(:) = 0.0d0

write(message,'(A,2g14.7)') 'sea_surface_x_stress_due_to_waves = ', &
minval(farrayPtr1), maxval(farrayPtr1)
call ESMF_LogWrite(trim(message), ESMF_LOGMSG_INFO)
Sw_tauox(:) = farrayPtr1(1:np)

call ESMF_StateGet(state, itemname='sea_surface_y_stress_due_to_waves', field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_FieldPtrUpdate(field, farrayPtr1, isPtr=isPtr, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

allocate(Sw_tauoy(np), stat=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
Sw_tauoy(:) = 0.0d0

write(message,'(A,2g14.7)') 'sea_surface_y_stress_due_to_waves = ', &
minval(farrayPtr1), maxval(farrayPtr1)
call ESMF_LogWrite(trim(message), ESMF_LOGMSG_INFO)
Sw_tauoy(:) = farrayPtr1(1:np)

call ESMF_StateGet(state, itemname='sea_bottom_upward_x_stress_due_to_waves', field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_FieldPtrUpdate(field, farrayPtr1, isPtr=isPtr, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

allocate(Sw_taubblx(np), stat=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
Sw_taubblx(:) = 0.0d0

write(message,'(A,2g14.7)') 'sea_bottom_upward_x_stress_due_to_waves = ', &
minval(farrayPtr1), maxval(farrayPtr1)
call ESMF_LogWrite(trim(message), ESMF_LOGMSG_INFO)
Sw_taubblx(:) = farrayPtr1(1:np)

call ESMF_StateGet(state, itemname='sea_bottom_upward_y_stress_due_to_waves', field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_FieldPtrUpdate(field, farrayPtr1, isPtr=isPtr, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

allocate(Sw_taubbly(np), stat=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
Sw_taubbly(:) = 0.0d0

write(message,'(A,2g14.7)') 'sea_bottom_upward_y_stress_due_to_waves = ', &
minval(farrayPtr1), maxval(farrayPtr1)
call ESMF_LogWrite(trim(message), ESMF_LOGMSG_INFO)
Sw_taubbly(:) = farrayPtr1(1:np)

call ESMF_StateGet(state, itemname='sea_bed_orbital_x_velocity_due_to_waves', field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_FieldPtrUpdate(field, farrayPtr1, isPtr=isPtr, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

allocate(Sw_ubrx(np), stat=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
Sw_ubrx(:) = 0.0d0

write(message,'(A,2g14.7)') 'sea_bed_orbital_x_velocity_due_to_waves = ', &
minval(farrayPtr1), maxval(farrayPtr1)
call ESMF_LogWrite(trim(message), ESMF_LOGMSG_INFO)
Sw_ubrx(:) = farrayPtr1(1:np)

call ESMF_StateGet(state, itemname='sea_bed_orbital_y_velocity_due_to_waves', field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_FieldPtrUpdate(field, farrayPtr1, isPtr=isPtr, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

allocate(Sw_ubry(np), stat=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
Sw_ubry(:) = 0.0d0

write(message,'(A,2g14.7)') 'sea_bed_orbital_y_velocity_due_to_waves = ', &
minval(farrayPtr1), maxval(farrayPtr1)
call ESMF_LogWrite(trim(message), ESMF_LOGMSG_INFO)
Sw_ubry(:) = farrayPtr1(1:np)

call ESMF_StateGet(state, itemname='sea_surface_wave_mean_direction', field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_FieldPtrUpdate(field, farrayPtr1, isPtr=isPtr, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

allocate(Sw_thm(np), stat=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
Sw_thm(:) = 0.0d0

write(message,'(A,2g14.7)') 'sea_surface_wave_mean_direction = ', &
minval(farrayPtr1), maxval(farrayPtr1)
call ESMF_LogWrite(trim(message), ESMF_LOGMSG_INFO)
Sw_thm(:) = farrayPtr1(1:np)

call ESMF_StateGet(state, itemname='sea_surface_wave_mean_period', field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_FieldPtrUpdate(field, farrayPtr1, isPtr=isPtr, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

allocate(Sw_t0m1(np), stat=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
Sw_t0m1(:) = 0.0d0

write(message,'(A,2g14.7)') 'sea_surface_wave_mean_period = ', &
minval(farrayPtr1), maxval(farrayPtr1)
call ESMF_LogWrite(trim(message), ESMF_LOGMSG_INFO)
Sw_t0m1(:) = farrayPtr1(1:np)

call ESMF_StateGet(state, itemname='sea_surface_wave_mean_number', field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_FieldPtrUpdate(field, farrayPtr1, isPtr=isPtr, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

allocate(Sw_wnmean(np), stat=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
Sw_wnmean(:) = 0.0d0

write(message,'(A,2g14.7)') 'sea_surface_wave_mean_number = ', &
minval(farrayPtr1), maxval(farrayPtr1)
call ESMF_LogWrite(trim(message), ESMF_LOGMSG_INFO)
Sw_wnmean(:) = farrayPtr1(1:np)

call ESMF_StateGet(state, itemname='eastward_surface_stokes_drift_current', field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_FieldPtrUpdate(field, farrayPtr1, isPtr=isPtr, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

allocate(Sw_ustokes(np), stat=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
Sw_ustokes(:) = 0.0d0

write(message,'(A,2g14.7)') 'eastward_surface_stokes_drift_current = ', &
minval(farrayPtr1), maxval(farrayPtr1)
call ESMF_LogWrite(trim(message), ESMF_LOGMSG_INFO)
Sw_ustokes(:) = farrayPtr1(1:np)

call ESMF_StateGet(state, itemname='northward_surface_stokes_drift_current', field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_FieldPtrUpdate(field, farrayPtr1, isPtr=isPtr, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

allocate(Sw_vstokes(np), stat=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
Sw_vstokes(:) = 0.0d0

write(message,'(A,2g14.7)') 'northward_surface_stokes_drift_current = ', &
minval(farrayPtr1), maxval(farrayPtr1)
call ESMF_LogWrite(trim(message), ESMF_LOGMSG_INFO)
Sw_vstokes(:) = farrayPtr1(1:np)

call get_WW3_arrays(Sw_hs,Sw_thm,Sw_t0m1,Sw_wnmean,Sw_bhd,Sw_ustokes,Sw_vstokes,&
Sw_tauox,Sw_tauoy,Sw_taubblx,Sw_taubbly,Sw_ubrx,Sw_ubry)

end subroutine SCHISM_StateImportWave3dVortex

end module schism_esmf_util
85 changes: 81 additions & 4 deletions src/schism/schism_nuopc_cap.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -379,6 +379,8 @@ subroutine InitializeAdvertise(comp, importState, exportState, clock, rc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
call NUOPC_FieldDictionaryAddIfNeeded("ocn_current_merid", "m s-1", localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
call NUOPC_FieldDictionaryAddIfNeeded("sea_surface_height_above_sea_level", "m", localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
call NUOPC_FieldDictionaryAddIfNeeded("ocean_mask", "1", localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

Expand All @@ -388,6 +390,19 @@ subroutine InitializeAdvertise(comp, importState, exportState, clock, rc)
call NUOPC_Advertise(importState, "inst_merid_wind_height10m", rc=localrc)

! for coupling to WW3/WDAT
call NUOPC_Advertise(importState, "sea_surface_wave_significant_height", rc=localrc)
call NUOPC_Advertise(importState, "sea_water_waves_effect_on_currents_bernoulli_head_adjustment", rc=localrc)
call NUOPC_Advertise(importState, "sea_surface_x_stress_due_to_waves", rc=localrc)
call NUOPC_Advertise(importState, "sea_surface_y_stress_due_to_waves", rc=localrc)
call NUOPC_Advertise(importState, "sea_bottom_upward_x_stress_due_to_waves", rc=localrc)
call NUOPC_Advertise(importState, "sea_bottom_upward_y_stress_due_to_waves", rc=localrc)
call NUOPC_Advertise(importState, "sea_bed_orbital_x_velocity_due_to_waves", rc=localrc)
call NUOPC_Advertise(importState, "sea_bed_orbital_y_velocity_due_to_waves", rc=localrc)
call NUOPC_Advertise(importState, "sea_surface_wave_mean_direction", rc=localrc)
call NUOPC_Advertise(importState, "sea_surface_wave_mean_period", rc=localrc)
call NUOPC_Advertise(importState, "sea_surface_wave_mean_number", rc=localrc)
call NUOPC_Advertise(importState, "eastward_surface_stokes_drift_current", rc=localrc)
call NUOPC_Advertise(importState, "northward_surface_stokes_drift_current", rc=localrc)
call NUOPC_Advertise(importState, "eastward_wave_radiation_stress", rc=localrc)
call NUOPC_Advertise(importState, "eastward_northward_wave_radiation_stress", rc=localrc)
call NUOPC_Advertise(importState, "northward_wave_radiation_stress", rc=localrc)
Expand Down Expand Up @@ -429,6 +444,9 @@ subroutine InitializeAdvertise(comp, importState, exportState, clock, rc)
call NUOPC_FieldAdvertise(exportState, "depth-averaged_y-velocity", "m s-1", localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call NUOPC_FieldAdvertise(exportState, "sea_surface_height_above_sea_level", "m", localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

! required for coupling through CMEPS mediator
call NUOPC_FieldAdvertise(exportState, "ocean_mask", "1", localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
Expand Down Expand Up @@ -521,6 +539,58 @@ subroutine InitializeRealize(comp, importState, exportState, clock, rc)
name="downwelling_short_photosynthetic_radiation_at_water_surface", field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_StateFieldCreateRealize(comp, state=importState, &
name="sea_surface_wave_significant_height", field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_StateFieldCreateRealize(comp, state=importState, &
name="sea_water_waves_effect_on_currents_bernoulli_head_adjustment", field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_StateFieldCreateRealize(comp, state=importState, &
name="sea_surface_x_stress_due_to_waves", field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_StateFieldCreateRealize(comp, state=importState, &
name="sea_surface_y_stress_due_to_waves", field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_StateFieldCreateRealize(comp, state=importState, &
name="sea_bottom_upward_x_stress_due_to_waves", field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_StateFieldCreateRealize(comp, state=importState, &
name="sea_bottom_upward_y_stress_due_to_waves", field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_StateFieldCreateRealize(comp, state=importState, &
name="sea_bed_orbital_x_velocity_due_to_waves", field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_StateFieldCreateRealize(comp, state=importState, &
name="sea_bed_orbital_y_velocity_due_to_waves", field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_StateFieldCreateRealize(comp, state=importState, &
name="sea_surface_wave_mean_direction", field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_StateFieldCreateRealize(comp, state=importState, &
name="sea_surface_wave_mean_period", field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_StateFieldCreateRealize(comp, state=importState, &
name="sea_surface_wave_mean_number", field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_StateFieldCreateRealize(comp, state=importState, &
name="eastward_surface_stokes_drift_current", field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

call SCHISM_StateFieldCreateRealize(comp, state=importState, &
name="northward_surface_stokes_drift_current", field=field, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

!> @todo add more atmospheric fields (like humidity)

!> Wave parameters, for now we only have those from the WW3DATA cap in
Expand Down Expand Up @@ -967,8 +1037,9 @@ end subroutine SCHISM_Export
#define ESMF_METHOD "SCHISM_Import"
subroutine SCHISM_Import(comp, importState, clock, rc)

use schism_glbl , only: windx2, windy2, pr2
use schism_glbl , only: RADFLAG, windx2, windy2, pr2
use schism_esmf_util, only: SCHISM_StateImportWaveTensor
use schism_esmf_util, only: SCHISM_StateImportWave3dVortex
use schism_esmf_util, only: SCHISM_StateUpdate

implicit none
Expand Down Expand Up @@ -1001,9 +1072,15 @@ subroutine SCHISM_Import(comp, importState, clock, rc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)

!> Update fields on import state
!> Obtain radiation tensor from wave component and calculate the wave stress
call SCHISM_StateImportWaveTensor(importState, isDataPtr, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
if (RADFLAG == 'VOR') then
!> Obtain required variables from wave component to do coupling with vortex formulation
call SCHISM_StateImportWave3dVortex(importState, isDataPtr, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
else
!> Obtain radiation tensor from wave component and calculate the wave stress
call SCHISM_StateImportWaveTensor(importState, isDataPtr, rc=localrc)
_SCHISM_LOG_AND_FINALIZE_ON_ERROR_(rc)
end if

!> surface wind component in x direction
call SCHISM_StateUpdate(importState, 'inst_zonal_wind_height10m', windx2, &
Expand Down