Merge branch 'ridging_CS' of https://github.com/ESMG/SIS2 into ridgin…

…g_CS

src/SIS_dyn_bgrid.F90

@@ -22,7 +22,6 @@ module SIS_dyn_bgrid
 use SIS_restart,       only : register_restart_field, SIS_restart_CS
 use SIS_framework,     only : safe_alloc_ptr
 use SIS_hor_grid,      only : SIS_hor_grid_type
-use ice_ridging_mod,   only : ridge_rate
 
 implicit none ; private
 
@@ -635,17 +634,6 @@ subroutine SIS_B_dynamics(ci, misp, mice, ui, vi, uo, vo, &
     if (CS%id_vi>0) call post_SIS_data(CS%id_vi, vi, CS%diag)
   endif
 
-  ! Compute the deformation rate for ridging
-  if (do_ridging) then
-    do j=jsc,jec ; do i=isc,iec ; if (ice_present(i,j)) then
-      del2 = (strn11(i,j)*strn11(i,j) + strn22(i,j)*strn22(i,j)) * (1+EC2I) +  &
-             4.0*EC2I*strn12(i,j)*strn12(i,j) + 2.0*strn11(i,j)*strn22(i,j)*(1-EC2I) ! H&D eqn 9
-      rdg_rate(i,j) = ridge_rate(del2, (strn11(i,j)+strn22(i,j)))
-    else
-      rdg_rate(i,j) = 0.0
-    endif ; enddo ; enddo
-  endif
-
 end subroutine SIS_B_dynamics
 
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!

src/SIS_dyn_trans.F90

@@ -153,7 +153,7 @@ module SIS_dyn_trans
   integer :: max_nts      !< The maximum number of transport steps that can be stored
                           !! before they are carried out.
   integer :: nts = 0      !< The number of accumulated transport steps since the last update.
-  real :: ridge_rate_count !< The number of contributions to av_ridge_rate
+  real :: ridge_rate_count !< The number of contributions to avg_ridge_rate
 
   real, allocatable, dimension(:,:) :: avg_ridge_rate !< The time average ridging rate in [T-1 ~> s-1].
 
@@ -2284,10 +2284,10 @@ subroutine SIS_dyn_trans_init(Time, G, US, IG, param_file, diag, CS, output_dir,
       call SIS_B_dyn_init(CS%Time, G, US, param_file, CS%diag, CS%SIS_B_dyn_CSp)
     endif
     if (CS%merged_cont) then
-      call SIS_transport_init(CS%Time, G, US, param_file, CS%diag, CS%SIS_transport_CSp, &
+      call SIS_transport_init(CS%Time, G, IG, US, param_file, CS%diag, CS%SIS_transport_CSp, &
                               continuity_CSp=CS%continuity_CSp, cover_trans_CSp=CS%cover_trans_CSp)
     else
-      call SIS_transport_init(CS%Time, G, US, param_file, CS%diag, CS%SIS_transport_CSp, &
+      call SIS_transport_init(CS%Time, G, IG, US, param_file, CS%diag, CS%SIS_transport_CSp, &
                               continuity_CSp=CS%continuity_CSp)
     endif
 

src/SIS_ice_diags.F90

@@ -49,9 +49,9 @@ module SIS_ice_diags
   integer :: id_mib = -1, id_mi = -1
   integer, dimension(:), allocatable :: id_t, id_sal
   integer :: id_cn = -1, id_hi = -1, id_hp = -1, id_hs = -1, id_tsn = -1, id_ext = -1
-  integer :: id_t_iceav = -1, id_s_iceav = -1, id_e2m = -1, id_rdgf = -1
+  integer :: id_t_iceav = -1, id_s_iceav = -1, id_e2m = -1, id_rdgf = -1, id_rdg_h = -1
 
-  integer :: id_simass = -1, id_sisnmass = -1, id_sivol = -1
+  integer :: id_simass = -1, id_simassn = -1, id_sisnmass = -1, id_sivol = -1
   integer :: id_siconc = -1, id_sithick = -1, id_sisnconc = -1, id_sisnthick = -1
   integer :: id_siconc_CMOR = -1, id_sisnconc_CMOR = -1, id_sivol_CMOR = -1
   integer :: id_siu = -1, id_siv = -1, id_sispeed = -1, id_sitimefrac = -1
@@ -137,6 +137,7 @@ subroutine post_ice_state_diagnostics(IDs, IST, OSS, IOF, dt_slow, Time, G, US,
   ! Thermodynamic state diagnostics
   !
   if (IDs%id_cn>0) call post_data(IDs%id_cn, IST%part_size(:,:,1:ncat), diag)
+  if (IDs%id_simassn>0) call post_data(IDs%id_simassn, IST%mH_ice, diag)
   if ((IDs%id_siconc>0) .or. (IDs%id_siconc_CMOR>0)) then
     diagVar(:,:) = 0.0
     do j=jsc,jec ; do i=isc,iec ; do k=1,ncat
@@ -257,6 +258,11 @@ subroutine post_ice_state_diagnostics(IDs, IST, OSS, IOF, dt_slow, Time, G, US,
     call post_data(IDs%id_rdgf, rdg_frac, diag)
   endif
 
+  ! Dermine the height of ridged ice rdg_height in each category.
+  if (IDs%id_rdg_h>0) then
+    call post_data(IDs%id_rdg_h, IST%rdg_height, diag)
+  endif
+
 end subroutine post_ice_state_diagnostics
 
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
@@ -373,6 +379,8 @@ subroutine register_ice_state_diagnostics(Time, IG, US, param_file, diag, IDs)
                'ice + snow mass', 'kg/m^2', conversion=US%RZ_to_kg_m2, missing_value=missing)
   IDs%id_simass = register_diag_field('ice_model', 'simass', diag%axesT1, Time, &
                'ice mass', 'kg/m^2', conversion=US%RZ_to_kg_m2, missing_value=missing)
+  IDs%id_simassn = register_diag_field('ice_model', 'simass_n', diag%axesTc, Time, &
+               'ice mass', 'kg/m^2', conversion=US%RZ_to_kg_m2, missing_value=missing)
   IDs%id_sisnmass = register_diag_field('ice_model', 'sisnmass', diag%axesT1, Time, &
                'snow mass', 'kg/m^2', conversion=US%RZ_to_kg_m2, missing_value=missing)
   IDs%id_mib  = register_diag_field('ice_model', 'MIB', diag%axesT1, Time, &
@@ -385,6 +393,8 @@ subroutine register_ice_state_diagnostics(Time, IG, US, param_file, diag, IDs)
   if (do_ridging) then
     IDs%id_rdgf = register_diag_field('ice_model', 'RDG_FRAC', diag%axesTc, Time, &
                    'ridged ice fraction', '0-1', missing_value=missing)
+    IDs%id_rdg_h = register_diag_field('ice_model', 'RDG_HEIGHT', diag%axesTc, Time, &
+                   'ridged ice fraction', '0-1', missing_value=missing)
   endif
 end subroutine register_ice_state_diagnostics
 

src/SIS_transport.F90

@@ -25,7 +25,7 @@ module SIS_transport
 use SIS_tracer_registry, only : check_SIS_tracer_bounds
 use SIS_types,         only : ice_state_type
 use ice_grid,          only : ice_grid_type
-use ice_ridging_mod,   only : ice_ridging
+use ice_ridging_mod,   only : ice_ridging_init, ice_ridging, ice_ridging_CS
 
 implicit none ; private
 
@@ -65,6 +65,8 @@ module SIS_transport
           !< The control structure for the SIS tracer advection module
   type(SIS_tracer_advect_CS), pointer :: SIS_thick_adv_CSp => NULL()
           !< The control structure for the SIS thickness advection module
+  type(ice_ridging_CS),       pointer :: ice_ridging_CSp => NULL()
+          !< Pointer to the control structure for the ice ridging
 
   !>@{ Diagnostic IDs
   integer :: id_ix_trans = -1, id_iy_trans = -1, id_xprt = -1, id_rdgr = -1
@@ -262,7 +264,8 @@ subroutine finish_ice_transport(CAS, IST, TrReg, G, US, IG, dt, CS, rdg_rate)
 
   if (CS%do_ridging) then
     ! Compress the ice using the ridging scheme taken from the CICE-Icepack module
-    call ice_ridging(IST, G, IG, CAS%m_ice, CAS%m_snow, CAS%m_pond, TrReg, US, dt, IST%rdg_rate)
+    call ice_ridging(IST, G, IG, CAS%m_ice, CAS%m_snow, CAS%m_pond, TrReg, CS%ice_ridging_CSp, US, &
+                     dt, rdg_rate=IST%rdg_rate, rdg_height=IST%rdg_height)
     ! Clean up any residuals
     call compress_ice(IST%part_size, IST%mH_ice, IST%mH_snow, IST%mH_pond, TrReg, G, US, IG, CS, CAS)
   else
@@ -1122,11 +1125,12 @@ end subroutine get_total_enthalpy
 
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
 !> SIS_transport_init initializes the ice transport and sets parameters.
-subroutine SIS_transport_init(Time, G, US, param_file, diag, CS, continuity_CSp, cover_trans_CSp)
+subroutine SIS_transport_init(Time, G, IG, US, param_file, diag, CS, continuity_CSp, cover_trans_CSp)
   type(time_type),     target, intent(in)    :: Time !< The sea-ice model's clock,
                                                      !! set with the current model time.
   type(SIS_hor_grid_type),     intent(in)    :: G    !< The horizontal grid type
-  type(unit_scale_type),       intent(in)    :: US  !< A structure with unit conversion factors
+  type(ice_grid_type),         intent(in)    :: IG   !< The sea-ice specific grid type
+  type(unit_scale_type),       intent(in)    :: US   !< A structure with unit conversion factors
   type(param_file_type),       intent(in)    :: param_file !< A structure to parse for run-time parameters
   type(SIS_diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic output
   type(SIS_transport_CS),      pointer       :: CS   !< The control structure for this module
@@ -1174,7 +1178,6 @@ subroutine SIS_transport_init(Time, G, US, param_file, diag, CS, continuity_CSp,
                  "in categories with less than this coverage to be discarded.", &
                  units="nondim", default=-1.0)
 
-
   call get_param(param_file, mdl, "CHECK_ICE_TRANSPORT_CONSERVATION", CS%check_conservation, &
                  "If true, use add multiple diagnostics of ice and snow "//&
                  "mass conservation in the sea-ice transport code.  This "//&
@@ -1217,6 +1220,8 @@ subroutine SIS_transport_init(Time, G, US, param_file, diag, CS, continuity_CSp,
   call SIS_continuity_init(Time, G, US, param_file, diag, CS%continuity_CSp, &
                            CS_cvr=cover_trans_CSp)
   call SIS_tracer_advect_init(Time, G, param_file, diag, CS%SIS_tr_adv_CSp)
+  if (CS%do_ridging) &
+      call ice_ridging_init(G, IG, param_file, CS%ice_ridging_CSp, US)
 
   if (present(continuity_CSp)) continuity_CSp => CS%continuity_CSp
 

src/SIS_types.F90

@@ -95,6 +95,8 @@ module SIS_types
   real, allocatable, dimension(:,:) :: rdg_rate !< The rate of fractional area loss by ridging [T-1 ~> s-1]
   real, allocatable, dimension(:,:,:) :: &
     rdg_mice    !< A diagnostic of the ice load that was formed by ridging [R Z ~> kg m-2].
+  real, allocatable, dimension(:,:,:) :: &
+    rdg_height    ! height of ridged ice per category [Z ~> m]
 
   logical :: Cgrid_dyn !< If true use a C-grid discretization of the sea-ice dynamics.
   logical :: valid_IST !< If true, this is currently the valid state of the ice.  Otherwise the ice
@@ -458,6 +460,7 @@ subroutine alloc_IST_arrays(HI, IG, IST, omit_velocities, omit_Tsurf, do_ridging
     allocate(IST%enth_snow_to_ocn(isd:ied, jsd:jed), source=0.0)
     allocate(IST%rdg_rate(isd:ied, jsd:jed), source=0.0)
     allocate(IST%rdg_mice(isd:ied, jsd:jed, CatIce), source=0.0)
+    allocate(IST%rdg_height(isd:ied, jsd:jed, CatIce), source=0.0)
   endif ; endif
 
   if (do_vel) then
@@ -1225,7 +1228,7 @@ subroutine copy_IST_to_IST(IST_in, IST_out, HI_in, HI_out, IG)
     IST_out%sal_ice(i2,j2,k,m) = IST_in%sal_ice(i,j,k,m)
   enddo ; enddo ; enddo ; enddo
 
-  ! The velocity components, rdg_mice, TrReg, and ITV are deliberately not being copied.
+  ! The velocity components, rdg_mice, rdg_height, TrReg, and ITV are deliberately not being copied.
 
 end subroutine copy_IST_to_IST
 
@@ -1241,7 +1244,7 @@ subroutine redistribute_IST_to_IST(IST_in, IST_out, domain_in, domain_out)
   real, pointer, dimension(:,:,:) :: null_ptr3D => NULL()
   real, pointer, dimension(:,:,:,:) :: null_ptr4D => NULL()
 
-  ! The velocity components, rdg_mice, TrReg, and ITV are deliberately not being copied.
+  ! The velocity components, rdg_mice, rdg_height, TrReg, and ITV are deliberately not being copied.
   if (associated(IST_out) .and. associated(IST_in)) then
     call redistribute_data(domain_in, IST_in%part_size, domain_out, &
                            IST_out%part_size, complete=.true.)

src/ice_model.F90

@@ -57,7 +57,6 @@ module ice_model_mod
 use ice_type_mod,       only : ice_type_slow_reg_restarts, ice_type_fast_reg_restarts
 use ice_type_mod,       only : Ice_public_type_chksum, Ice_public_type_bounds_check
 use ice_type_mod,       only : ice_model_restart, ice_stock_pe, ice_data_type_chksum
-use ice_ridging_mod, only : ice_ridging_init
 
 use SIS_ctrl_types,    only : SIS_slow_CS, SIS_fast_CS
 use SIS_ctrl_types,    only : ice_diagnostics_init, ice_diags_fast_init
@@ -2137,8 +2136,6 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow,
                                massless_val=massless_ice_salin, nonnegative=.true.)
     endif
 
-    call ice_ridging_init(sG,sIG,sIST%TrReg,US)
-
   !   Register any tracers that will be handled via tracer flow control for
   ! restarts and advection.
     call SIS_call_tracer_register(sG, sIG, param_file, Ice%sCS%SIS_tracer_flow_CSp, &