diff --git a/src/schism/schism_esmf_util.F90 b/src/schism/schism_esmf_util.F90 index bb62c35..63767d0 100644 --- a/src/schism/schism_esmf_util.F90 +++ b/src/schism/schism_esmf_util.F90 @@ -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 @@ -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 diff --git a/src/schism/schism_nuopc_cap.F90 b/src/schism/schism_nuopc_cap.F90 index 49c8fd0..3246aa7 100644 --- a/src/schism/schism_nuopc_cap.F90 +++ b/src/schism/schism_nuopc_cap.F90 @@ -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) @@ -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) @@ -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) @@ -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 @@ -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 @@ -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, &