Major changes

Add type that contains EVP state, make changes to "smuggle" this
information into btstep via the mech_forcing type.

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

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, 

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