Skip to content

Commit

Permalink
Merge pull request #32 from oceanmodeling/feature/schism_3d
Browse files Browse the repository at this point in the history 
Changes to enable coupling with WW3 using 3d vortex formulation
  • Loading branch information
@josephzhang8
josephzhang8 authored Jan 29, 2025
2 parents 2c9d178 + 3f93003 commit d3ab56d
Show file tree
Hide file tree
Showing 2 changed files with 328 additions and 4 deletions.
247 changes: 247 additions & 0 deletions src/schism/schism_esmf_util.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,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 @@ -2648,4 +2649,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
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

0 comments on commit d3ab56d

Please sign in to comment.