Skip to content

Commit

Permalink
This version works and has good ocean state
Browse files Browse the repository at this point in the history 
This version has good ocean state, but some issues with the sea ice.

Recent changes include the improved copying from IOB to the seaice type
used in mom6.
  • Loading branch information
Theresa Morrison authored and Theresa Morrison committed Jan 17, 2025
1 parent b8f2fd6 commit b1fa533
Show file tree
Hide file tree
Showing 5 changed files with 206 additions and 36 deletions.
123 changes: 115 additions & 8 deletions config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1172,13 +1172,118 @@ subroutine convert_merged_ice(seaice_in, seaice_out)
seaice_out = seaice_in
end subroutine

subroutine extract_merged_ice_from_IOB(iob, seaice)
subroutine extract_merged_ice_from_IOB(IOB, seaice, index_bounds, Time, G, US, CS, tau_halo)
type(ice_ocean_boundary_type), &
intent(in) :: iob !< An ice-ocean boundary type with fluxes to drive the
target, intent(in) :: IOB !< An ice-ocean boundary type with fluxes to drive
!! the ocean in a coupled model
type(SIS_dyn_state_2d), &
intent(inout) :: seaice
intent(inout) :: seaice
integer, optional, intent(in) :: tau_halo !< The halo size of wind stresses to set, 0 by default.
integer, dimension(4), intent(in) :: index_bounds !< The i- and j- size of the arrays in IOB.
type(time_type), intent(in) :: Time !< The time of the fluxes, used for interpolating the
!! salinity to the right time, when it is being restored.
type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(surface_forcing_CS),pointer :: CS !< A pointer to the control structure returned by a
!! previous call to surface_forcing_init.
! Local variables
real, dimension(SZIB_(G),SZJ_(G)) :: u_in_C ! Zonal velocity
!real, dimension(SZIB_(G),SZJ_(G)) :: uh_step !
real, dimension(SZIB_(G),SZJ_(G)) :: WindStr_x !
real, dimension(SZIB_(G),SZJ_(G)) :: WindStr_ocn_x !

real, dimension(SZI_(G),SZJB_(G)) :: v_in_C ! Meridional velocity [R Z L T-2 ~> Pa] at v points
!real, dimension(SZI_(G),SZJB_(G)) :: vh_step
real, dimension(SZI_(G),SZJB_(G)) :: WindStr_y
real, dimension(SZI_(G),SZJB_(G)) :: WindStr_ocn_y

real, dimension(SZI_(G),SZJ_(G)) :: &
mi_sum, &
ice_free, &
ice_cover, &
FIA_ice_free, &
FIA_ice_cover

integer :: i, j, is, ie, js, je, ish, ieh, jsh, jeh, Isqh, Ieqh, Jsqh, Jeqh, i0, j0, halo

halo = tau_halo
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
ish = G%isc-halo ; ieh = G%iec+halo ; jsh = G%jsc-halo ; jeh = G%jec+halo
Isqh = G%IscB-halo ; Ieqh = G%IecB+halo ; Jsqh = G%JscB-halo ; Jeqh = G%JecB+halo
i0 = is - index_bounds(1) ; j0 = js - index_bounds(3)

!is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
!Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
!isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
!IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
!isr = is-isd+1 ; ier = ie-isd+1 ; jsr = js-jsd+1 ; jer = je-jsd+1

!i0 = is - isc_bnd ; j0 = js - jsc_bnd
!do j=js,je ; do i=is,ie
! fluxes%lprec(i,j) = kg_m2_s_conversion * IOB%lprec(i-i0,j-j0) * G%mask2dT(i,j)
!enddo; enddo

! Set surface momentum stress related fields as a function of staggering.
u_in_C(:,:) = 0.0 ; v_in_C(:,:) = 0.0
ice_cover(:,:) = 0.0 ; mi_sum(:,:) = 0.0
WindStr_x(:,:) = 0.0 ; WindStr_ocn_x(:,:) = 0.0
WindStr_y(:,:) = 0.0 ; WindStr_ocn_y(:,:) = 0.0
FIA_ice_cover(:,:) = 0.0 ; FIA_ice_free(:,:) = 0.0
do j=js,je ; do i=is,ie
mi_sum(i,J) = IOB%IceDS2d%mi_sum(i-i0,j-j0)
ice_cover(i,J) = IOB%IceDS2d%ice_cover(i-i0,j-j0)
FIA_ice_free(i,J) = IOB%IceDS2d%FIA_2d%ice_free(i-i0,j-j0)
FIA_ice_cover(i,J) = IOB%IceDS2d%FIA_2d%ice_cover(i-i0,j-j0)

u_in_C(I,j) = IOB%IceDS2d%u_ice_C(i-i0,j-j0)
v_in_C(i,J) = IOB%IceDS2d%v_ice_C(i-i0,j-j0)
!uh_step(I,j,:) = IOB%IceDS2d%uh_step(i-i0,j-j0,:)
!vh_step(i,J,:) = IOB%IceDS2d%vh_step(i-i0,j-j0,:)

WindStr_x(I,j) = IOB%IceDS2d%FIA_2d%WindStr_x(i-i0,j-j0)
WindStr_y(i,J) = IOB%IceDS2d%FIA_2d%WindStr_y(i-i0,j-j0)
WindStr_ocn_x(I,j) = IOB%IceDS2d%FIA_2d%WindStr_ocn_x(i-i0,j-j0)
WindStr_ocn_y(i,J) = IOB%IceDS2d%FIA_2d%WindStr_ocn_y(i-i0,j-j0)
enddo ; enddo

if (G%symmetric) call fill_symmetric_edges(u_in_C, v_in_C, G%Domain)
!if (G%symmetric) call fill_symmetric_edges(uh_step, vh_step, G%Domain)
if (G%symmetric) call fill_symmetric_edges(WindStr_x, WindStr_y, G%Domain)
if (G%symmetric) call fill_symmetric_edges(WindStr_ocn_x, WindStr_ocn_y, G%Domain)
call pass_vector(u_in_C, v_in_C, G%Domain, halo=max(1,halo))
!call pass_vector(uh_step, vh_step, G%Domain, halo=max(1,halo))
call pass_vector(WindStr_x, WindStr_y, G%Domain, halo=max(1,halo))
call pass_vector(WindStr_ocn_x, WindStr_ocn_y, G%Domain, halo=max(1,halo))
call pass_var(mi_sum, G%Domain, halo=max(1,halo))
!call pass_var(ice_free, G%Domain, halo=max(1,halo))
call pass_var(ice_cover, G%Domain, halo=max(1,halo))
call pass_var(FIA_ice_free, G%Domain, halo=max(1,halo))
call pass_var(FIA_ice_cover, G%Domain, halo=max(1,halo))

do j=jsh,jeh ; do i=ish,ieh
seaice%ice_cover(i,j) = G%mask2dT(i,j)*ice_cover(i,j)
seaice%mi_sum(i,j) = G%mask2dT(i,j)*mi_sum(i,j)
seaice%FIA_2d%ice_free(i,j) = G%mask2dT(i,j)*FIA_ice_free(i,j)
seaice%FIA_2d%ice_cover(i,j) = G%mask2dT(i,j)*FIA_ice_cover(i,j)
enddo ; enddo
do j=jsh,jeh ; do I=Isqh,Ieqh
seaice%u_ice_C(I,j) = G%mask2dCu(I,j)*u_in_C(I,j)
!seaice%uh_step(I,j,:) = G%mask2dCu(I,j)*uh_step(I,j,:)
seaice%FIA_2d%WindStr_x(I,j) = G%mask2dCu(I,j)*WindStr_x(I,j)
seaice%FIA_2d%WindStr_ocn_x(I,j) = G%mask2dCu(I,j)*WindStr_ocn_x(I,j)
enddo ; enddo
do J=Jsqh,Jeqh ; do i=ish,ieh
seaice%v_ice_C(i,J) = G%mask2dCv(i,J)*v_in_C(i,J)
!seaice%vh_step(i,J,:) = G%mask2dCv(i,J)*vh_step(i,J,:)
seaice%FIA_2d%WindStr_y(i,J) = G%mask2dCv(i,J)*WindStr_y(i,J)
seaice%FIA_2d%WindStr_ocn_y(i,J) = G%mask2dCv(i,J)*WindStr_ocn_y(i,J)
enddo ; enddo

seaice = iob%IceDS2d
seaice%SIS_C_dyn_CSp = IOB%IceDS2d%SIS_C_dyn_CSp
seaice%dynmer_trans_CSp = IOB%IceDS2d%dynmer_trans_CSp
seaice%sG = IOB%IceDS2d%sG
seaice%IG = IOB%IceDS2d%IG
seaice%US = IOB%IceDS2d%US

! seaice%max_nts = iob%IceDS2d
! seaice%nts = iob%
Expand All @@ -1198,7 +1303,7 @@ subroutine extract_merged_ice_from_IOB(iob, seaice)
! seaice%FIA_2d = iob%
! seaice%SIS_C_dyn_CSp = iob%
! seaice%dynmer_trans_CSp = iob%
! seaice%sG = iob%
! seaice%sG = iob%y
! seaice%IG = iob%
! seaice%US = iob%

Expand Down Expand Up @@ -1342,9 +1447,11 @@ subroutine seaice_init(G, DS2d)

! Local variables
integer :: i, j, isd, ied, jsd, jed
integer :: isc, iec, jsc, jec
integer :: isdB, iedB, jsdB, jedB

isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec
isdB = G%isdB ; iedB = G%iedB ; jsdB = G%jsdB ; jedB = G%jedB

if (.not.associated(DS2d)) allocate(DS2d)
Expand All @@ -1363,8 +1470,8 @@ subroutine seaice_init(G, DS2d)
call safe_alloc(DS2d%u_ice_C, G%IsdB, G%IedB, G%jsd, G%jed)
call safe_alloc(DS2d%v_ice_C, G%isd, G%ied, G%JsdB, G%JedB)

!call safe_alloc(DS2d%uh_step, G%IsdB, G%IedB, G%jsd, G%jed)
!call safe_alloc(DS2d%vh_step, G%isd, G%ied, G%JsdB, G%JedB)
call safe_alloc(DS2d%uh_step, G%IsdB, G%IedB, G%jsd, G%jed, 3)
call safe_alloc(DS2d%vh_step, G%isd, G%ied, G%JsdB, G%JedB, 3)
!call safe_alloc(DS2d%mca_step, G%isd, G%ied, G%jsd, G%jed)
! only needed for embedded/semiembedded coupling
if (.not.associated(DS2d%FIA_2d)) allocate(DS2d%FIA_2d)
Expand All @@ -1384,7 +1491,7 @@ subroutine seaice_init(G, DS2d)
if (.not.associated(DS2d%dynmer_trans_CSp%cover_trans_CSp)) allocate(DS2d%dynmer_trans_CSp%cover_trans_CSp)
!DS2d%dynmer_trans_CSp%cover_trans_CSp = CS%cover_trans_CSp
!CS%DS2d%SIS_C_dyn_CSp => CS%SIS_C_dyn_CSp
!! if (.not.associated(DS2d%SIS_C_dyn_CSp)) allocate(DS2d%SIS_C_dyn_CSp)
if (.not.associated(DS2d%SIS_C_dyn_CSp)) allocate(DS2d%SIS_C_dyn_CSp)
!! DS2d%SIS_C_dyn_CSp = CS%SIS_C_dyn_CSp
if (.not.associated(DS2d%sG)) allocate(DS2d%sG)
!DS2d%sG = G
Expand Down
3 changes: 2 additions & 1 deletion config_src/drivers/FMS_cap/ocean_model_MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -553,7 +553,8 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
call iceberg_forces(OS%grid, OS%forces, OS%use_ice_shelf, &
OS%sfc_state, dt_coupling, OS%marine_ice_CSp)
if (OS%use_dynmer_ice) &
call extract_merged_ice_from_IOB(Ice_ocean_boundary, OS%seaice) ! where will this be put?
call extract_merged_ice_from_IOB(Ice_ocean_boundary, OS%seaice, index_bnds, OS%Time_dyn, &
OS%grid, OS%US, OS%forcing_CSp, tau_halo=1)
endif

if (do_thermo) then
Expand Down
17 changes: 17 additions & 0 deletions src/core/MOM_barotropic.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -317,6 +317,7 @@ module MOM_barotropic
!>@{ Diagnostic IDs
integer :: id_PFu_bt = -1, id_PFv_bt = -1, id_Coru_bt = -1, id_Corv_bt = -1
integer :: id_ubtforce = -1, id_vbtforce = -1, id_uaccel = -1, id_vaccel = -1
integer :: id_ubtforcetaux = -1, id_vbtforcetauy = -1
integer :: id_visc_rem_u = -1, id_visc_rem_v = -1, id_eta_cor = -1
integer :: id_ubt = -1, id_vbt = -1, id_eta_bt = -1, id_ubtav = -1, id_vbtav = -1
integer :: id_ubt_st = -1, id_vbt_st = -1, id_eta_st = -1
Expand Down Expand Up @@ -567,6 +568,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
! bt_rem_u is between 0 and 1.
BT_force_u, & ! The vertical average of all of the u-accelerations that are
! not explicitly included in the barotropic equation [L T-2 ~> m s-2].
BT_force_taux, & !
!
u_accel_bt, & ! The difference between the zonal acceleration from the
! barotropic calculation and BT_force_u [L T-2 ~> m s-2].
uhbt, & ! The zonal barotropic thickness fluxes [H L2 T-1 ~> m3 s-1 or kg s-1].
Expand Down Expand Up @@ -598,6 +601,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
! drag [nondim]. bt_rem_v is between 0 and 1.
BT_force_v, & ! The vertical average of all of the v-accelerations that are
! not explicitly included in the barotropic equation [L T-2 ~> m s-2].
BT_force_tauy, & !
!
v_accel_bt, & ! The difference between the meridional acceleration from the
! barotropic calculation and BT_force_v [L T-2 ~> m s-2].
vhbt, & ! The meridional barotropic thickness fluxes [H L2 T-1 ~> m3 s-1 or kg s-1].
Expand Down Expand Up @@ -1339,11 +1344,14 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,

if (present(taux_surf)) then
BT_force_u(I,j) = taux_surf(I,j) * GV%RZ_to_H * CS%IDatu(I,j)*visc_rem_u(I,j,1)
BT_force_taux(I,j) = taux_surf(I,j) * GV%RZ_to_H * CS%IDatu(I,j)*visc_rem_u(I,j,1)
else
BT_force_u(I,j) = forces%taux(I,j) * GV%RZ_to_H * CS%IDatu(I,j)*visc_rem_u(I,j,1)
BT_force_taux(I,j) = forces%taux(I,j) * GV%RZ_to_H * CS%IDatu(I,j)*visc_rem_u(I,j,1)
endif
else
BT_force_u(I,j) = 0.0
BT_force_taux(I,j) = 0.0
endif ; enddo ; enddo
!$OMP parallel do default(shared)
do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.0) then
Expand All @@ -1369,11 +1377,14 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,

if (present(tauy_surf)) then
BT_force_v(i,J) = tauy_surf(i,J) * GV%RZ_to_H * CS%IDatv(i,J)*visc_rem_v(i,J,1)
BT_force_tauy(i,J) = tauy_surf(i,J) * GV%RZ_to_H * CS%IDatv(i,J)*visc_rem_v(i,J,1)
else
BT_force_v(i,J) = forces%tauy(i,J) * GV%RZ_to_H * CS%IDatv(i,J)*visc_rem_v(i,J,1)
BT_force_tauy(i,J) = forces%tauy(i,J) * GV%RZ_to_H * CS%IDatv(i,J)*visc_rem_v(i,J,1)
endif
else
BT_force_v(i,J) = 0.0
BT_force_tauy(i,J) = 0.0
endif ; enddo ; enddo
if (associated(taux_bot) .and. associated(tauy_bot)) then
!$OMP parallel do default(shared)
Expand Down Expand Up @@ -2293,6 +2304,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,

if (CS%id_ubtforce > 0) call post_data(CS%id_ubtforce, BT_force_u(IsdB:IedB,jsd:jed), CS%diag)
if (CS%id_vbtforce > 0) call post_data(CS%id_vbtforce, BT_force_v(isd:ied,JsdB:JedB), CS%diag)
if (CS%id_ubtforcetaux > 0) call post_data(CS%id_ubtforcetaux, BT_force_taux(IsdB:IedB,jsd:jed), CS%diag)
if (CS%id_vbtforcetauy > 0) call post_data(CS%id_vbtforcetauy, BT_force_tauy(isd:ied,JsdB:JedB), CS%diag)
if (CS%id_uaccel > 0) call post_data(CS%id_uaccel, u_accel_bt(IsdB:IedB,jsd:jed), CS%diag)
if (CS%id_vaccel > 0) call post_data(CS%id_vaccel, v_accel_bt(isd:ied,JsdB:JedB), CS%diag)

Expand Down Expand Up @@ -5854,6 +5867,10 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
'Barotropic zonal acceleration from baroclinic terms', 'm s-2', conversion=US%L_T2_to_m_s2)
CS%id_vbtforce = register_diag_field('ocean_model', 'vbtforce', diag%axesCv1, Time, &
'Barotropic meridional acceleration from baroclinic terms', 'm s-2', conversion=US%L_T2_to_m_s2)
CS%id_ubtforcetaux = register_diag_field('ocean_model', 'ubt_taux', diag%axesCu1, Time, &
'Acceleration from surface stress, x', 'm s-2', conversion=US%L_T2_to_m_s2)
CS%id_vbtforcetauy = register_diag_field('ocean_model', 'vbt_tauy', diag%axesCv1, Time, &
'Acceleration from surface stress, y', 'm s-2', conversion=US%L_T2_to_m_s2)
CS%id_ubtdt = register_diag_field('ocean_model', 'ubt_dt', diag%axesCu1, Time, &
'Barotropic zonal acceleration', 'm s-2', conversion=US%L_T2_to_m_s2)
CS%id_vbtdt = register_diag_field('ocean_model', 'vbt_dt', diag%axesCv1, Time, &
Expand Down
Loading

0 comments on commit b1fa533

Please sign in to comment.