From 5f132711fe592a7cececa13cbba0bb51cbfabb81 Mon Sep 17 00:00:00 2001 From: Theresa Morrison Date: Mon, 12 Feb 2024 14:51:16 -0500 Subject: [PATCH] Major changes Add type that contains EVP state, make changes to "smuggle" this information into btstep via the mech_forcing type. --- .../FMS_cap/MOM_surface_forcing_gfdl.F90 | 5 ++ src/core/MOM_barotropic.F90 | 29 +++++------ src/core/MOM_forcing_type.F90 | 50 +++++++++++++++++++ 3 files changed, 70 insertions(+), 14 deletions(-) diff --git a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 index 1e56486329..ea44efe5e5 100644 --- a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 +++ b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 @@ -40,6 +40,7 @@ module MOM_surface_forcing_gfdl use user_revise_forcing, only : user_alter_forcing, user_revise_forcing_init use user_revise_forcing, only : user_revise_forcing_CS use iso_fortran_env, only : int64 +use MOM_forcing_type, only : SIS_C_EVP_state implicit none ; private @@ -208,6 +209,7 @@ module MOM_surface_forcing_gfdl !! This flag may be set by the flux-exchange code, based on what !! the sea-ice model is providing. Otherwise, the value from !! the surface_forcing_CS is used. + type(SIS_C_EVP_state) :: EVP_type end type ice_ocean_boundary_type integer :: id_clock_forcing !< A CPU time clock @@ -716,6 +718,9 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_ kg_m2_s_conversion = US%kg_m2s_to_RZ_T + ! copy EVP info to something that will be able to make it to MOM_barotropic + forces%EVPT = IOB%EVP_type + ! allocation and initialization if this is the first time that this ! mechanical forcing type has been used. if (.not.forces%initialized) then diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 250e299d7b..19ee1568cb 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -543,6 +543,9 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! averaged between the beginning and end of the time step, ! relative to eta_PF, with SAL effects included [H ~> m or kg m-2]. + real, dimension(SZIB_(G),SZJ_(G)) :: fxoc ! ice-to-ocn stress + real, dimension(SZI_(G),SZJB_(G)) :: fyoc ! + ! These are always allocated with symmetric memory and wide halos. real :: q(SZIBW_(CS),SZJBW_(CS)) ! A pseudo potential vorticity [T-1 H-1 ~> s-1 m-1 or m2 s-1 kg-1] real, dimension(SZIBW_(CS),SZJW_(CS)) :: & @@ -1859,10 +1862,12 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! Update the surface velocity estimate do j=js,je ; do I=is-1,ie - uo(i,j) = uo_st(i,j) + (dt*(n*dtbt/dt))*(u_accel_bt(i,j)) + !uo(i,j) = uo_st(i,j) + (n*dtbt)*(u_accel_bt(i,j) + wt_u(I,j,1) * bc_accel_u(I,j,1)) + uo(i,j) = uo_st(i,j) + (n*dtbt)*(u_accel_bt(i,j) + BT_force_u(i,j)) enddo ; enddo do j=js-1,je ; do I=is,ie - vo(i,j) = vo_st(i,j) + (dt*(n*dtbt/dt))*(v_accel_bt(i,j)) + !vo(i,j) = vo_st(i,j) + (n*dtbt)*(v_accel_bt(i,j) + wt_u(i,j,1) * bc_accel_v(I,j,1)) + vo(i,j) = vo_st(i,j) + (n*dtbt)*(v_accel_bt(i,j) + BT_force_v(i,j)) enddo ; enddo if (find_PF) then @@ -1975,13 +1980,15 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! calculate each layer's accelerations for !call btstep_layer_accel(dt, u_accel_bt, v_accel_bt, pbce, gtot_E, gtot_W, gtot_N, gtot_S, & - ! e_anom, G, GV, CS, accel_layer_uo, accel_layer_vo) + ! eta_wtd, G, GV, CS, accel_layer_uo, accel_layer_vo) !do j=js,je ; do I=is-1,ie - ! uo_cont(i,j) = uo_st(i,j) + (dt*(n*dtbt/dt))*(accel_layer_uo(i,j,1) + u_accel_bt(i,j)) + ! !uo_cont(i,j) = uo_st(i,j) + (n*dtbt)*(accel_layer_uo(i,j,1) + u_accel_bt(i,j) + wt_u(i,j,1) * bc_accel_u(i,j,1)) + ! uo_cont(i,j) = uo_st(i,j) + (n*dtbt)*(accel_layer_uo(i,j,1) + u_accel_bt(i,j) + BT_force_u(i,j)) !enddo ; enddo !do j=js-1,je ; do I=is,ie - ! vo_cont(i,j) = vo_st(i,j) + (dt*(n*dtbt/dt))*(accel_layer_vo(i,j,1) + v_accel_bt(i,j)) + ! !vo_cont(i,j) = vo_st(i,j) + (n*dtbt)*(accel_layer_vo(i,j,1) + v_accel_bt(i,j) + wt_u(i,j,1) * bc_accel_v(i,j,1)) + ! vo_cont(i,j) = vo_st(i,j) + (n*dtbt)*(accel_layer_vo(i,j,1) + v_accel_bt(i,j) + BT_force_v(i,j)) !enddo ; enddo if (do_hifreq_output) then @@ -1995,8 +2002,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, if (CS%id_eta_pred_hifreq > 0) call post_data(CS%id_eta_pred_hifreq, eta_PF_BT(isd:ied,jsd:jed), CS%diag) if (CS%id_uo_hifreq > 0) call post_data(CS%id_uo_hifreq, uo(IsdB:IedB,jsd:jed), CS%diag) if (CS%id_vo_hifreq > 0) call post_data(CS%id_vo_hifreq, vo(isd:ied,JsdB:JedB), CS%diag) - !if (CS%id_uo_cont_hifreq > 0) call post_data(CS%id_uo_cont_hifreq, uo_cont(IsdB:IedB,jsd:jed), CS%diag) - !if (CS%id_vo_cont_hifreq > 0) call post_data(CS%id_vo_cont_hifreq, vo_cont(isd:ied,JsdB:JedB), CS%diag) + ! if (CS%id_uo_cont_hifreq > 0) call post_data(CS%id_uo_cont_hifreq, uo_cont(IsdB:IedB,jsd:jed), CS%diag) + ! if (CS%id_vo_cont_hifreq > 0) call post_data(CS%id_vo_cont_hifreq, vo_cont(isd:ied,JsdB:JedB), CS%diag) endif if (CS%debug_bt) then @@ -2155,13 +2162,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! TJM ! The EVP function from the SIS2 model -! call EVP_step_loop(dt_slow, ci, ui, vi, mice, uo, vo, & -! fxat, fyat, fxoc, fyoc, pres_mice, diag_val, del_sh_min_pr, & -! ui_min_trunc, ui_max_trunc, vi_min_trunc, vi_max_trunc, & -! mi_u, f2dt_u, I1_f2dt2_u, PFu, mi_v, f2dt_v, I1_f2dt2_v, PFv, & -! azon, bzon, czon, dzon, amer, bmer, cmer, dmer, & -! mi_ratio_A_q, & -! G, US, CS) + call EVP_step_loop(forces, uo, vo, PFu, PFv, fxoc, fyoc) ! In need from sea-ice: ci, ui, vi, mice, pres_mice, diag_val (?), del_sh_min_pr (?), ! ui_min_trunc, ui_max_trunc, vi_min_trunc, vi_max_trunc, ! mi_u, f2dt_u, I1_f2dt2_u, mi_v, f2dt_v, I1_f2dt2_v, diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index dcbf440292..c7353fea9c 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -64,6 +64,55 @@ module MOM_forcing_type ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units ! vary with the Boussinesq approximation, the Boussinesq variant is given first. +!> The control structure with the state that is used and updated in the EVP loop +type, public :: SIS_C_EVP_state + real, allocatable, dimension(:,:) :: ci !< Sea ice concentration [nondim] + real, allocatable, dimension(:,:) :: ui !< Zonal ice velocity [L T-1 ~> m s-1] + real, allocatable, dimension(:,:) :: vi !< Meridional ice velocity [L T-1 ~> m s-1] + real, allocatable, dimension(:,:) :: mice !< Mass per unit ocean area of sea ice [R Z ~> kg m-2] + real, allocatable, dimension(:,:) :: fxat !< Zonal air stress on ice [R Z L T-2 ~> Pa] + real, allocatable, dimension(:,:) :: fyat !< Meridional air stress on ice [R Z L T-2 ~> Pa] + + real, allocatable, dimension(:,:) :: & + pres_mice, & ! The ice internal pressure per unit column mass [L2 T-2 ~> N m kg-1]. + diag_val, & ! A temporary diagnostic array. + del_sh_min_pr ! When multiplied by pres_mice, this gives the minimum + ! value of del_sh that is used in the calculation of zeta [T-1 ~> s-1]. + ! This is set based on considerations of numerical stability, + ! and varies with the grid spacing. + + real, allocatable, dimension(:,:) :: & + ui_min_trunc, & ! The range of v-velocities beyond which the velocities + ui_max_trunc, & ! are truncated [L T-1 ~> m s-1], or 0 for land cells + mi_u, & ! The total ice and snow mass interpolated to u points [R Z ~> kg m-2]. + f2dt_u, &! The squared effective Coriolis parameter at u-points times a + ! time step [T-1 ~> s-1]. + PFu, & ! Zonal hydrostatic pressure driven acceleration [L T-2 ~> m s-2]. + I1_f2dt2_u ! 1 / ( 1 + f^2 dt^2) at u-points [nondim]. + + real, allocatable, dimension(:,:) :: & + vi_min_trunc, & ! The range of v-velocities beyond which the velocities + vi_max_trunc, & ! are truncated [L T-1 ~> m s-1], or 0 for land cells. + mi_v, & ! The total ice and snow mass interpolated to v points [R Z ~> kg m-2]. + f2dt_v, &! The squared effective Coriolis parameter at v-points times a + ! time step [T-1 ~> s-1]. + PFv, & ! hydrostatic pressure driven acceleration [L T-2 ~> m s-2]. + I1_f2dt2_v ! 1 / ( 1 + f^2 dt^2) at v-points [nondim]. + + real, allocatable, dimension(:,:) :: & + azon, bzon, & ! _zon & _mer are the values of the Coriolis force which + czon, dzon, & ! are applied to the neighboring values of vi & ui, + amer, bmer, & ! respectively to get the barotropic inertial rotation, + cmer, dmer ! in units of [T-1 ~> s-1]. azon and amer couple the same pair of + ! velocities, but with the influence going in opposite + ! directions. + + real, allocatable, dimension(:,:) :: & + mi_ratio_A_q ! A ratio of the masses interpolated to the faces around a + ! vorticity point that ranges between (4 mi_min/mi_max) and 1, + ! divided by the sum of the ocean areas around a point [L-2 ~> m-2]. +end type SIS_C_EVP_state + !> Structure that contains pointers to the boundary forcing used to drive the !! liquid ocean simulated by MOM. !! @@ -293,6 +342,7 @@ module MOM_forcing_type !! 3rd dimension - wavenumber logical :: initialized = .false. !< This indicates whether the appropriate arrays have been initialized. + type(SIS_C_EVP_state) :: EVPT !! the info for the sea-ice end type mech_forcing !> Structure that defines the id handles for the forcing type