From 9835058e527910ccbcd219f22dd1346a45f01c1f Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 5 Mar 2021 10:59:16 -0500 Subject: [PATCH 1/3] Corrected openMP issues with landfast_ice PR Corrected two runtime issues with missing openMP directives by combining expressions to eliminate the recently added variables that were missing from openMP directive lists. All answers are bitwise identical in all cases that compiled. --- src/SIS_dyn_cgrid.F90 | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/src/SIS_dyn_cgrid.F90 b/src/SIS_dyn_cgrid.F90 index 22940969ff..ad4bd68827 100644 --- a/src/SIS_dyn_cgrid.F90 +++ b/src/SIS_dyn_cgrid.F90 @@ -632,7 +632,6 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, fxat, fyat, & real :: m_neglect2 ! A tiny mass per unit area squared [R2 Z2 ~> kg2 m-4]. real :: m_neglect4 ! A tiny mass per unit area to the 4th power [R4 Z4 ~> kg4 m-8]. real :: sum_area ! The sum of ocean areas around a vorticity point [L2 ~> m2]. - real :: drag_bot_u, drag_bot_v ! A drag with the bottom at u & v points [R Z T-1 ~> kg m-2 s-1]. type(time_type) :: & time_it_start, & ! The starting time of the iteratve steps. @@ -1060,9 +1059,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, fxat, fyat, & endif if (drag_max>0.) drag_u = min( drag_u, drag_max ) if (CS%lemieux_landfast) then - drag_bot_u = CS%Tb_u(I,j) / & - (sqrt(ui(I,j)**2 + v2_at_u_min ) + CS%lemieux_u0) - drag_u = drag_u + drag_bot_u + drag_u = drag_u + CS%Tb_u(I,j) / (sqrt(ui(I,j)**2 + v2_at_u_min ) + CS%lemieux_u0) endif ! This is a quasi-implicit timestep of Coriolis, followed by an explicit @@ -1154,9 +1151,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, fxat, fyat, & endif if (drag_max>0.) drag_v = min( drag_v, drag_max ) if (CS%lemieux_landfast) then - drag_bot_v = CS%Tb_v(i,J) / & - (sqrt(vi(i,J)**2 + u2_at_v_min ) + CS%lemieux_u0) - drag_v = drag_v + drag_bot_v + drag_v = drag_v + CS%Tb_v(i,J) / (sqrt(vi(i,J)**2 + u2_at_v_min ) + CS%lemieux_u0) endif ! This is a quasi-implicit timestep of Coriolis, followed by an explicit From 4d0d910e0f87eac6fcc15401081ee5f4cd021a9a Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Sat, 13 Mar 2021 09:14:37 -0500 Subject: [PATCH 2/3] Modifications to enable Icepack APIs for ridging. - These modifications enable ridging code for SIS2 C-grid configurations (not yet implemented for B-grid) using APIs provided by the Icepack module from the CICE consortium (https://github.com/CICE-Consortium/Icepack) - These APIs are reproduced in the config_src/external/Icepack directory in order to remove SIS2 dependency on Icepack by default. In order to activate Icepack, remove the config_src/external/Icepack diretory from your compile list, compile a shared Icepack library and link to SIS2. - Sea-ice ridging rate diagnostics (rdg_rate) report the fractional rate of sea-ice coverage decrease (s-1) due to the ridging scheme. - WARNING:Sea-ice enthalpy diagnostics indicate a decrease in numerical closure from roundoff (10-17) without ridging enabled to values greater than roundoff (10-7) in initial testing. Mass and salinity, however, are consistent to roundoff. The issue should be addressed before using this option to generate scientific results. - The default Icepack ridging options are enabled. Additional options can be enabled easily using the SIS_parameter file and calls to icepack_init_parameters --- .../Icepack_interfaces/icepack_itd.F90 | 136 +++ .../Icepack_interfaces/icepack_kinds.F90 | 38 + .../Icepack_interfaces/icepack_mechred.F90 | 194 ++++ .../Icepack_interfaces/icepack_tracers.F90 | 119 +++ .../Icepack_interfaces/icepack_warnings.F90 | 242 +++++ src/SIS_dyn_trans.F90 | 14 +- src/SIS_transport.F90 | 12 +- src/SIS_types.F90 | 6 +- src/ice_model.F90 | 3 + src/ice_ridge.F90 | 836 +++++++++--------- 10 files changed, 1159 insertions(+), 441 deletions(-) create mode 100644 config_src/external/Icepack_interfaces/icepack_itd.F90 create mode 100644 config_src/external/Icepack_interfaces/icepack_kinds.F90 create mode 100644 config_src/external/Icepack_interfaces/icepack_mechred.F90 create mode 100644 config_src/external/Icepack_interfaces/icepack_tracers.F90 create mode 100644 config_src/external/Icepack_interfaces/icepack_warnings.F90 diff --git a/config_src/external/Icepack_interfaces/icepack_itd.F90 b/config_src/external/Icepack_interfaces/icepack_itd.F90 new file mode 100644 index 0000000000..dbfcb2ad5a --- /dev/null +++ b/config_src/external/Icepack_interfaces/icepack_itd.F90 @@ -0,0 +1,136 @@ + module icepack_itd + + use icepack_kinds +! use icepack_parameters, only: c0, c1, c2, c3, c15, c25, c100, p1, p01, p001, p5, puny +! use icepack_parameters, only: Lfresh, rhos, ice_ref_salinity, hs_min, cp_ice, Tocnfrz, rhoi +! use icepack_parameters, only: rhosi, sk_l, hs_ssl, min_salin +! use icepack_tracers, only: nt_Tsfc, nt_qice, nt_qsno, nt_aero, nt_isosno, nt_isoice +! use icepack_tracers, only: nt_apnd, nt_hpnd, nt_fbri, tr_brine, nt_bgc_S, bio_index +! use icepack_tracers, only: n_iso +! use icepack_tracers, only: tr_iso +! use icepack_tracers, only: icepack_compute_tracers +! use icepack_parameters, only: solve_zsal, skl_bgc, z_tracers +! use icepack_parameters, only: kcatbound, kitd +! use icepack_therm_shared, only: Tmin, hi_min +! use icepack_warnings, only: warnstr, icepack_warnings_add +! use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted +! use icepack_zbgc_shared,only: zap_small_bgc + + implicit none + + private + public :: cleanup_itd ,icepack_init_itd + + + +!======================================================================= + + contains + +!======================================================================= + + subroutine icepack_init_itd(ncat, hin_max) + + integer (kind=int_kind), intent(in) :: & + ncat ! number of thickness categories + + real (kind=dbl_kind), intent(out) :: & + hin_max(0:ncat) ! category limits (m) + + + end subroutine icepack_init_itd + + subroutine cleanup_itd (dt, ntrcr, & + nilyr, nslyr, & + ncat, hin_max, & + aicen, trcrn, & + vicen, vsnon, & + aice0, aice, & + n_aero, & + nbtrcr, nblyr, & + tr_aero, & + tr_pond_topo, & + heat_capacity, & + first_ice, & + trcr_depend, trcr_base, & + n_trcr_strata,nt_strata, & + fpond, fresh, & + fsalt, fhocn, & + faero_ocn, fiso_ocn, & + fzsal, & + flux_bio, limit_aice_in) + + integer (kind=int_kind), intent(in) :: & + ncat , & ! number of thickness categories + nilyr , & ! number of ice layers + nblyr , & ! number of bio layers + nslyr , & ! number of snow layers + ntrcr , & ! number of tracers in use + nbtrcr, & ! number of bio tracers in use + n_aero ! number of aerosol tracers + + real (kind=dbl_kind), intent(in) :: & + dt ! time step + + real (kind=dbl_kind), dimension(0:ncat), intent(in) :: & + hin_max ! category boundaries (m) + + real (kind=dbl_kind), dimension (:), intent(inout) :: & + aicen , & ! concentration of ice + vicen , & ! volume per unit area of ice (m) + vsnon ! volume per unit area of snow (m) + + real (kind=dbl_kind), dimension (:,:), intent(inout) :: & + trcrn ! ice tracers + + real (kind=dbl_kind), intent(inout) :: & + aice , & ! total ice concentration + aice0 ! concentration of open water + + integer (kind=int_kind), dimension (:), intent(in) :: & + trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon + n_trcr_strata ! number of underlying tracer layers + + real (kind=dbl_kind), dimension (:,:), intent(in) :: & + trcr_base ! = 0 or 1 depending on tracer dependency + ! argument 2: (1) aice, (2) vice, (3) vsno + + integer (kind=int_kind), dimension (:,:), intent(in) :: & + nt_strata ! indices of underlying tracer layers + + logical (kind=log_kind), intent(in) :: & + tr_aero, & ! aerosol flag + tr_pond_topo, & ! topo pond flag + heat_capacity ! if false, ice and snow have zero heat capacity + + logical (kind=log_kind), dimension(ncat), intent(inout) :: & + first_ice ! For bgc and S tracers. set to true if zapping ice. + + ! ice-ocean fluxes (required for strict conservation) + + real (kind=dbl_kind), intent(inout), optional :: & + fpond , & ! fresh water flux to ponds (kg/m^2/s) + fresh , & ! fresh water flux to ocean (kg/m^2/s) + fsalt , & ! salt flux to ocean (kg/m^2/s) + fhocn , & ! net heat flux to ocean (W/m^2) + fzsal ! net salt flux to ocean from zsalinity (kg/m^2/s) + + real (kind=dbl_kind), dimension (:), intent(inout), optional :: & + flux_bio ! net tracer flux to ocean from biology (mmol/m^2/s) + + real (kind=dbl_kind), dimension (:), intent(inout), optional :: & + faero_ocn ! aerosol flux to ocean (kg/m^2/s) + + real (kind=dbl_kind), dimension (:), intent(inout) :: & + fiso_ocn ! isotope flux to ocean (kg/m^2/s) + + logical (kind=log_kind), intent(in), optional :: & + limit_aice_in ! if false, allow aice to be out of bounds + ! may want to allow this for unit tests + + + end subroutine cleanup_itd + + end module icepack_itd + +!======================================================================= diff --git a/config_src/external/Icepack_interfaces/icepack_kinds.F90 b/config_src/external/Icepack_interfaces/icepack_kinds.F90 new file mode 100644 index 0000000000..664643fa3f --- /dev/null +++ b/config_src/external/Icepack_interfaces/icepack_kinds.F90 @@ -0,0 +1,38 @@ +!======================================================================= + +! Defines variable precision for all common data types +! Code originally based on kinds_mod.F in POP +! +! author: Elizabeth C. Hunke and William H. Lipscomb, LANL +! 2006: ECH converted to free source form (F90) +! 2020: TC added support for NO_I8 and NO_R16 + + module icepack_kinds + +!======================================================================= + + implicit none + public + + integer, parameter :: char_len = 80, & + char_len_long = 256, & + log_kind = kind(.true.), & + int_kind = selected_int_kind(6), & +#if defined (NO_I8) + int8_kind = selected_int_kind(6), & +#else + int8_kind = selected_int_kind(13), & +#endif + real_kind = selected_real_kind(6), & + dbl_kind = selected_real_kind(13), & +#if defined (NO_R16) + r16_kind = selected_real_kind(13) +#else + r16_kind = selected_real_kind(33,4931) +#endif + +!======================================================================= + + end module icepack_kinds + +!======================================================================= diff --git a/config_src/external/Icepack_interfaces/icepack_mechred.F90 b/config_src/external/Icepack_interfaces/icepack_mechred.F90 new file mode 100644 index 0000000000..d95bc1725b --- /dev/null +++ b/config_src/external/Icepack_interfaces/icepack_mechred.F90 @@ -0,0 +1,194 @@ + module icepack_mechred + + use icepack_kinds + use icepack_tracers, only : n_iso, n_aero + + implicit none + + private + public :: ridge_ice + contains + + subroutine ridge_ice (dt, ndtd, & + ncat, n_aero, & + nilyr, nslyr, & + ntrcr, hin_max, & + rdg_conv, rdg_shear, & + aicen, trcrn, & + vicen, vsnon, & + aice0, & + trcr_depend, trcr_base, & + n_trcr_strata, & + nt_strata, & + krdg_partic, krdg_redist,& + mu_rdg, tr_brine, & + dardg1dt, dardg2dt, & + dvirdgdt, opening, & + fpond, & + fresh, fhocn, & + faero_ocn, fiso_ocn, & + aparticn, krdgn, & + aredistn, vredistn, & + dardg1ndt, dardg2ndt, & + dvirdgndt, & + araftn, vraftn, & + closing_flag,closing ) + + integer (kind=int_kind), intent(in) :: & + ndtd , & ! number of dynamics subcycles + ncat , & ! number of thickness categories + nilyr , & ! number of ice layers + nslyr , & ! number of snow layers + n_aero, & ! number of aerosol tracers + ntrcr ! number of tracers in use + + real (kind=dbl_kind), intent(in) :: & + mu_rdg , & ! gives e-folding scale of ridged ice (m^.5) + dt ! time step + + real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: & + hin_max ! category limits (m) + + real (kind=dbl_kind), intent(in) :: & + rdg_conv , & ! normalized energy dissipation due to convergence (1/s) + rdg_shear ! normalized energy dissipation due to shear (1/s) + + real (kind=dbl_kind), dimension (:), intent(inout) :: & + aicen , & ! concentration of ice + vicen , & ! volume per unit area of ice (m) + vsnon ! volume per unit area of snow (m) + + real (kind=dbl_kind), dimension (:,:), intent(inout) :: & + trcrn ! ice tracers + + real (kind=dbl_kind), intent(inout) :: & + aice0 ! concentration of open water + + integer (kind=int_kind), dimension (:), intent(in) :: & + trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon + n_trcr_strata ! number of underlying tracer layers + + real (kind=dbl_kind), dimension (:,:), intent(in) :: & + trcr_base ! = 0 or 1 depending on tracer dependency + ! argument 2: (1) aice, (2) vice, (3) vsno + + integer (kind=int_kind), dimension (:,:), intent(in) :: & + nt_strata ! indices of underlying tracer layers + + integer (kind=int_kind), intent(in) :: & + krdg_partic, & ! selects participation function + krdg_redist ! selects redistribution function + + logical (kind=log_kind), intent(in) :: & + closing_flag, &! flag if closing is valid + tr_brine ! if .true., brine height differs from ice thickness + + ! optional history fields + real (kind=dbl_kind), intent(inout), optional :: & + dardg1dt , & ! rate of fractional area loss by ridging ice (1/s) + dardg2dt , & ! rate of fractional area gain by new ridges (1/s) + dvirdgdt , & ! rate of ice volume ridged (m/s) + opening , & ! rate of opening due to divergence/shear (1/s) + closing , & ! rate of closing due to divergence/shear (1/s) + fpond , & ! fresh water flux to ponds (kg/m^2/s) + fresh , & ! fresh water flux to ocean (kg/m^2/s) + fhocn ! net heat flux to ocean (W/m^2) + + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & + dardg1ndt , & ! rate of fractional area loss by ridging ice (1/s) + dardg2ndt , & ! rate of fractional area gain by new ridges (1/s) + dvirdgndt , & ! rate of ice volume ridged (m/s) + aparticn , & ! participation function + krdgn , & ! mean ridge thickness/thickness of ridging ice + araftn , & ! rafting ice area + vraftn , & ! rafting ice volume + aredistn , & ! redistribution function: fraction of new ridge area + vredistn ! redistribution function: fraction of new ridge volume + + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & + faero_ocn ! aerosol flux to ocean (kg/m^2/s) + + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & + fiso_ocn ! isotope flux to ocean (kg/m^2/s) + + ! local variables + + real (kind=dbl_kind), dimension (ncat) :: & + eicen ! energy of melting for each ice layer (J/m^2) + + real (kind=dbl_kind), dimension (ncat) :: & + esnon, & ! energy of melting for each snow layer (J/m^2) + vbrin, & ! ice volume with defined by brine height (m) + sicen ! Bulk salt in h ice (ppt*m) + + real (kind=dbl_kind) :: & + asum , & ! sum of ice and open water area + aksum , & ! ratio of area removed to area ridged + msnow_mlt , & ! mass of snow added to ocean (kg m-2) + esnow_mlt , & ! energy needed to melt snow in ocean (J m-2) + mpond , & ! mass of pond added to ocean (kg m-2) + closing_net, & ! net rate at which area is removed (1/s) + ! (ridging ice area - area of new ridges) / dt + divu_adv , & ! divu as implied by transport scheme (1/s) + opning , & ! rate of opening due to divergence/shear + ! opning is a local variable; + ! opening is the history diagnostic variable + ardg1 , & ! fractional area loss by ridging ice + ardg2 , & ! fractional area gain by new ridges + virdg , & ! ice volume ridged + aopen ! area opening due to divergence/shear + + real (kind=dbl_kind), dimension (n_aero) :: & + maero ! aerosol mass added to ocean (kg m-2) + + real (kind=dbl_kind), dimension (n_iso) :: & + miso ! isotope mass added to ocean (kg m-2) + + real (kind=dbl_kind), dimension (0:ncat) :: & + apartic ! participation function; fraction of ridging + ! and closing associated w/ category n + + real (kind=dbl_kind), dimension (ncat) :: & + hrmin , & ! minimum ridge thickness + hrmax , & ! maximum ridge thickness (krdg_redist = 0) + hrexp , & ! ridge e-folding thickness (krdg_redist = 1) + krdg , & ! mean ridge thickness/thickness of ridging ice + ardg1n , & ! area of ice ridged + ardg2n , & ! area of new ridges + virdgn , & ! ridging ice volume + mraftn ! rafting ice mask + + real (kind=dbl_kind) :: & + vice_init, vice_final, & ! ice volume summed over categories + vsno_init, vsno_final, & ! snow volume summed over categories + eice_init, eice_final, & ! ice energy summed over layers + vbri_init, vbri_final, & ! ice volume in fbri*vicen summed over categories + sice_init ,sice_final, & ! ice bulk salinity summed over categories + esno_init, esno_final ! snow energy summed over layers + + integer (kind=int_kind), parameter :: & + nitermax = 20 ! max number of ridging iterations + + integer (kind=int_kind) :: & + n , & ! thickness category index + niter , & ! iteration counter + k , & ! vertical index + it ! tracer index + + real (kind=dbl_kind) :: & + dti ! 1 / dt + + logical (kind=log_kind) :: & + iterate_ridging ! if true, repeat the ridging + + character (len=char_len) :: & + fieldid ! field identifier + + + + end subroutine ridge_ice + + + end module icepack_mechred + +!======================================================================= diff --git a/config_src/external/Icepack_interfaces/icepack_tracers.F90 b/config_src/external/Icepack_interfaces/icepack_tracers.F90 new file mode 100644 index 0000000000..fec2b128cf --- /dev/null +++ b/config_src/external/Icepack_interfaces/icepack_tracers.F90 @@ -0,0 +1,119 @@ +!======================================================================= +! Indices and flags associated with the tracer infrastructure. +! Grid-dependent and max_trcr-dependent arrays are declared in ice_state.F90. +! +! author Elizabeth C. Hunke, LANL + + module icepack_tracers + + use icepack_kinds + + implicit none + + private + public :: icepack_init_tracer_indices + + real(kind=dbl_kind), parameter, public :: n_iso=0,n_aero=0 + + contains + + +!======================================================================= +!autodocument_start icepack_init_tracer_indices +! set the number of column tracer indices + + subroutine icepack_init_tracer_indices(& + nt_Tsfc_in, nt_qice_in, nt_qsno_in, nt_sice_in, & + nt_fbri_in, nt_iage_in, nt_FY_in, & + nt_alvl_in, nt_vlvl_in, nt_apnd_in, nt_hpnd_in, nt_ipnd_in, & + nt_fsd_in, nt_isosno_in, nt_isoice_in, & + nt_aero_in, nt_zaero_in, nt_bgc_C_in, & + nt_bgc_N_in, nt_bgc_chl_in, nt_bgc_DOC_in, nt_bgc_DON_in, & + nt_bgc_DIC_in, nt_bgc_Fed_in, nt_bgc_Fep_in, nt_bgc_Nit_in, nt_bgc_Am_in, & + nt_bgc_Sil_in, nt_bgc_DMSPp_in, nt_bgc_DMSPd_in, nt_bgc_DMS_in, nt_bgc_hum_in, & + nt_bgc_PON_in, nlt_zaero_in, nlt_bgc_C_in, nlt_bgc_N_in, nlt_bgc_chl_in, & + nlt_bgc_DOC_in, nlt_bgc_DON_in, nlt_bgc_DIC_in, nlt_bgc_Fed_in, & + nlt_bgc_Fep_in, nlt_bgc_Nit_in, nlt_bgc_Am_in, nlt_bgc_Sil_in, & + nlt_bgc_DMSPp_in, nlt_bgc_DMSPd_in, nlt_bgc_DMS_in, nlt_bgc_hum_in, & + nlt_bgc_PON_in, nt_zbgc_frac_in, nt_bgc_S_in, nlt_chl_sw_in, & + nlt_zaero_sw_in, & + bio_index_o_in, bio_index_in) + + integer, intent(in), optional :: & + nt_Tsfc_in, & ! ice/snow temperature + nt_qice_in, & ! volume-weighted ice enthalpy (in layers) + nt_qsno_in, & ! volume-weighted snow enthalpy (in layers) + nt_sice_in, & ! volume-weighted ice bulk salinity (CICE grid layers) + nt_fbri_in, & ! volume fraction of ice with dynamic salt (hinS/vicen*aicen) + nt_iage_in, & ! volume-weighted ice age + nt_FY_in, & ! area-weighted first-year ice area + nt_alvl_in, & ! level ice area fraction + nt_vlvl_in, & ! level ice volume fraction + nt_apnd_in, & ! melt pond area fraction + nt_hpnd_in, & ! melt pond depth + nt_ipnd_in, & ! melt pond refrozen lid thickness + nt_fsd_in, & ! floe size distribution + nt_isosno_in, & ! starting index for isotopes in snow + nt_isoice_in, & ! starting index for isotopes in ice + nt_aero_in, & ! starting index for aerosols in ice + nt_bgc_Nit_in, & ! nutrients + nt_bgc_Am_in, & ! + nt_bgc_Sil_in, & ! + nt_bgc_DMSPp_in,&! trace gases (skeletal layer) + nt_bgc_DMSPd_in,&! + nt_bgc_DMS_in, & ! + nt_bgc_hum_in, & ! + nt_bgc_PON_in, & ! zooplankton and detritus + nlt_bgc_Nit_in,& ! nutrients + nlt_bgc_Am_in, & ! + nlt_bgc_Sil_in,& ! + nlt_bgc_DMSPp_in,&! trace gases (skeletal layer) + nlt_bgc_DMSPd_in,&! + nlt_bgc_DMS_in,& ! + nlt_bgc_hum_in,& ! + nlt_bgc_PON_in,& ! zooplankton and detritus + nt_zbgc_frac_in,&! fraction of tracer in the mobile phase + nt_bgc_S_in, & ! Bulk salinity in fraction ice with dynamic salinity (Bio grid)) + nlt_chl_sw_in ! points to total chla in trcrn_sw + + integer (kind=int_kind), dimension(:), intent(in), optional :: & + bio_index_o_in, & + bio_index_in + + integer (kind=int_kind), dimension(:), intent(in), optional :: & + nt_bgc_N_in , & ! diatoms, phaeocystis, pico/small + nt_bgc_C_in , & ! diatoms, phaeocystis, pico/small + nt_bgc_chl_in, & ! diatoms, phaeocystis, pico/small + nlt_bgc_N_in , & ! diatoms, phaeocystis, pico/small + nlt_bgc_C_in , & ! diatoms, phaeocystis, pico/small + nlt_bgc_chl_in ! diatoms, phaeocystis, pico/small + + integer (kind=int_kind), dimension(:), intent(in), optional :: & + nt_bgc_DOC_in, & ! dissolved organic carbon + nlt_bgc_DOC_in ! dissolved organic carbon + + integer (kind=int_kind), dimension(:), intent(in), optional :: & + nt_bgc_DON_in, & ! dissolved organic nitrogen + nlt_bgc_DON_in ! dissolved organic nitrogen + + integer (kind=int_kind), dimension(:), intent(in), optional :: & + nt_bgc_DIC_in, & ! dissolved inorganic carbon + nlt_bgc_DIC_in ! dissolved inorganic carbon + + integer (kind=int_kind), dimension(:), intent(in), optional :: & + nt_bgc_Fed_in, & ! dissolved iron + nt_bgc_Fep_in, & ! particulate iron + nlt_bgc_Fed_in,& ! dissolved iron + nlt_bgc_Fep_in ! particulate iron + + integer (kind=int_kind), dimension(:), intent(in), optional :: & + nt_zaero_in, & ! black carbon and other aerosols + nlt_zaero_in, & ! black carbon and other aerosols + nlt_zaero_sw_in ! black carbon and dust in trcrn_sw + + end subroutine icepack_init_tracer_indices + + + end module icepack_tracers + +!======================================================================= diff --git a/config_src/external/Icepack_interfaces/icepack_warnings.F90 b/config_src/external/Icepack_interfaces/icepack_warnings.F90 new file mode 100644 index 0000000000..cf0faf82d3 --- /dev/null +++ b/config_src/external/Icepack_interfaces/icepack_warnings.F90 @@ -0,0 +1,242 @@ + +module icepack_warnings + + use icepack_kinds + + implicit none + + private + + ! warning messages + character(len=char_len_long), dimension(:), allocatable :: warnings + integer :: nWarnings = 0 + integer, parameter :: nWarningsBuffer = 10 ! incremental number of messages + + ! abort flag, accessed via icepack_warnings_setabort and icepack_warnings_aborted + logical :: warning_abort = .false. + + ! public string for all subroutines to use + character(len=char_len_long), public :: warnstr + + public :: & + icepack_warnings_clear, & + icepack_warnings_print, & + icepack_warnings_flush, & + icepack_warnings_aborted, & + icepack_warnings_add, & + icepack_warnings_setabort + + private :: & + icepack_warnings_getall, & + icepack_warnings_getone + +!======================================================================= + +contains + +!======================================================================= +!autodocument_start icepack_warnings_aborted +! turn on the abort flag in the icepack warnings package +! pass in an optional error message + + logical function icepack_warnings_aborted(instring) + + character(len=*),intent(in), optional :: instring + +!autodocument_end + + character(len=*),parameter :: subname='(icepack_warnings_aborted)' + + icepack_warnings_aborted = warning_abort + if (warning_abort .and. present(instring)) then + call icepack_warnings_add(subname//' ... '//trim(instring)) + endif + + end function icepack_warnings_aborted + +!======================================================================= + + subroutine icepack_warnings_setabort(abortflag,file,line) + + logical, intent(in) :: abortflag + character(len=*), intent(in), optional :: file + integer, intent(in), optional :: line + + character(len=*),parameter :: subname='(icepack_warnings_setabort)' + + ! try to capture just the first setabort call + + if (abortflag) then + write(warnstr,*) subname,abortflag + if (present(file)) write(warnstr,*) trim(warnstr)//' :file '//trim(file) + if (present(line)) write(warnstr,*) trim(warnstr)//' :line ',line + call icepack_warnings_add(warnstr) + endif + + warning_abort = abortflag + + end subroutine icepack_warnings_setabort + +!======================================================================= +!autodocument_start icepack_warnings_clear +! clear all warning messages from the icepack warning buffer + + subroutine icepack_warnings_clear() + +!autodocument_end + + character(len=*),parameter :: subname='(icepack_warnings_clear)' + + nWarnings = 0 + + end subroutine icepack_warnings_clear + +!======================================================================= + + subroutine icepack_warnings_getall(warningsOut) + + character(len=char_len_long), dimension(:), allocatable, intent(out) :: & + warningsOut + + integer :: iWarning + character(len=*),parameter :: subname='(icepack_warnings_getall)' + + if (allocated(warningsOut)) deallocate(warningsOut) + allocate(warningsOut(nWarnings)) + + do iWarning = 1, nWarnings + warningsOut(iWarning) = trim(icepack_warnings_getone(iWarning)) + enddo + + end subroutine icepack_warnings_getall + +!======================================================================= +!autodocument_start icepack_warnings_print +! print all warning messages from the icepack warning buffer + + subroutine icepack_warnings_print(iounit) + + integer, intent(in) :: iounit + +!autodocument_end + + integer :: iWarning + character(len=*),parameter :: subname='(icepack_warnings_print)' + +! tcraig +! this code intermittenly aborts on recursive IO errors with intel +! not sure if it's OMP or something else causing this +!$OMP MASTER + do iWarning = 1, nWarnings + write(iounit,*) trim(icepack_warnings_getone(iWarning)) + enddo +!$OMP END MASTER + + end subroutine icepack_warnings_print + +!======================================================================= +!autodocument_start icepack_warnings_flush +! print and clear all warning messages from the icepack warning buffer + + subroutine icepack_warnings_flush(iounit) + + integer, intent(in) :: iounit + +!autodocument_end + + character(len=*),parameter :: subname='(icepack_warnings_flush)' + + if (nWarnings > 0) then + call icepack_warnings_print(iounit) + endif + call icepack_warnings_clear() + + end subroutine icepack_warnings_flush + +!======================================================================= + + subroutine icepack_warnings_add(warning) + + character(len=*), intent(in) :: warning ! warning to add to array of warnings + + ! local + + character(len=char_len_long), dimension(:), allocatable :: warningsTmp + integer :: & + nWarningsArray, & ! size of warnings array at start + iWarning ! warning index + character(len=*),parameter :: subname='(icepack_warnings_add)' + +!$OMP CRITICAL (omp_warnings_add) + ! check if warnings array is not allocated + if (.not. allocated(warnings)) then + + ! allocate warning array with number of buffer elements + allocate(warnings(nWarningsBuffer)) + + ! set initial number of nWarnings + nWarnings = 0 + + ! already allocated + else + + ! find the size of the warnings array at the start + nWarningsArray = size(warnings) + + ! check to see if need more space in warnings array + if (nWarnings + 1 > nWarningsArray) then + + ! allocate the temporary warning storage + allocate(warningsTmp(nWarningsArray)) + + ! copy the warnings to temporary storage + do iWarning = 1, nWarningsArray + warningsTmp(iWarning) = trim(warnings(iWarning)) + enddo ! iWarning + + ! increase the size of the warning array by the buffer size + deallocate(warnings) + allocate(warnings(nWarningsArray + nWarningsBuffer)) + + ! copy back the temporary stored warnings + do iWarning = 1, nWarningsArray + warnings(iWarning) = trim(warningsTmp(iWarning)) + enddo ! iWarning + + ! deallocate the temporary storage + deallocate(warningsTmp) + + endif + + endif + + ! increase warning number + nWarnings = nWarnings + 1 +!$OMP END CRITICAL (omp_warnings_add) + + ! add the new warning + warnings(nWarnings) = trim(warning) + + end subroutine icepack_warnings_add + +!======================================================================= + + function icepack_warnings_getone(iWarning) result(warning) + + integer, intent(in) :: iWarning + + character(len=char_len_long) :: warning + + character(len=*),parameter :: subname='(icepack_warnings_getone)' + + if (iWarning <= nWarnings) then + warning = warnings(iWarning) + else + warning = "" + endif + + end function icepack_warnings_getone + +!======================================================================= + +end module icepack_warnings diff --git a/src/SIS_dyn_trans.F90 b/src/SIS_dyn_trans.F90 index ab7e2d36a9..fc29656110 100644 --- a/src/SIS_dyn_trans.F90 +++ b/src/SIS_dyn_trans.F90 @@ -612,11 +612,12 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, U if (DS2d%nts==0) then if (CS%do_ridging) then - call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, US, IG, CS%SIS_transport_CSp, & - rdg_rate=DS2d%avg_ridge_rate) + call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, US, IG, dt_slow, CS%SIS_transport_CSp, & + ! rdg_rate=DS2d%avg_ridge_rate) + rdg_rate=IST%rdg_rate) DS2d%ridge_rate_count = 0. ; DS2d%avg_ridge_rate(:,:) = 0.0 else - call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, US, IG, CS%SIS_transport_CSp) + call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, US, IG, dt_slow,CS%SIS_transport_CSp) endif endif call cpu_clock_end(iceClock8) @@ -749,11 +750,12 @@ subroutine complete_IST_transport(DS2d, CAS, IST, dt_adv_cycle, G, US, IG, CS) ! Convert the cell-averaged state back to the ice-state type, adjusting the ! category mass distributions, doing ridging, and updating the partition sizes. if (CS%do_ridging) then - call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, US, IG, CS%SIS_transport_CSp, & - rdg_rate=DS2d%avg_ridge_rate) + call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, US, IG, dt_adv_cycle, CS%SIS_transport_CSp, & + ! rdg_rate=DS2d%avg_ridge_rate) + rdg_rate=IST%rdg_rate) DS2d%ridge_rate_count = 0. ; DS2d%avg_ridge_rate(:,:) = 0.0 else - call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, US, IG, CS%SIS_transport_CSp) + call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, US, IG, dt_adv_cycle, CS%SIS_transport_CSp) endif DS2d%nts = 0 ! There is no outstanding transport to be done and IST is up-to-date. diff --git a/src/SIS_transport.F90 b/src/SIS_transport.F90 index 6454e01282..c8558a6964 100644 --- a/src/SIS_transport.F90 +++ b/src/SIS_transport.F90 @@ -222,11 +222,12 @@ end subroutine ice_cat_transport !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> finish_ice_transport completes the ice transport and thickness class redistribution -subroutine finish_ice_transport(CAS, IST, TrReg, G, US, IG, CS, rdg_rate) +subroutine finish_ice_transport(CAS, IST, TrReg, G, US, IG, dt, CS, rdg_rate) type(cell_average_state_type), intent(inout) :: CAS !< A structure with ocean-cell averaged masses. type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type + real, intent(in) :: dt !< The timestep used for ridging [T -> s]. type(SIS_tracer_registry_type), pointer :: TrReg !< The registry of SIS ice and snow tracers. type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_transport_CS), pointer :: CS !< A pointer to the control structure for this module @@ -259,12 +260,17 @@ subroutine finish_ice_transport(CAS, IST, TrReg, G, US, IG, CS, rdg_rate) ! Convert the ocean-cell averaged properties back into the ice_state_type. call cell_ave_state_to_ice_state(CAS, G, US, IG, CS, IST, TrReg) + 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) + else + ! Compress the ice where the fractional coverage exceeds 1, starting with the ! thinnest category, in what amounts to a minimalist version of a sea-ice ! ridging scheme. A more complete ridging scheme would also compress ! thicker ice and allow the fractional ice coverage to drop below 1. - call compress_ice(IST%part_size, IST%mH_ice, IST%mH_snow, IST%mH_pond, TrReg, G, US, IG, CS, CAS) - + call compress_ice(IST%part_size, IST%mH_ice, IST%mH_snow, IST%mH_pond, TrReg, G, US, IG, CS, CAS) + endif if (CS%bounds_check) call check_SIS_tracer_bounds(TrReg, G, IG, "After compress_ice") if (CS%readjust_categories) then diff --git a/src/SIS_types.F90 b/src/SIS_types.F90 index 1e605cff83..56823a4714 100644 --- a/src/SIS_types.F90 +++ b/src/SIS_types.F90 @@ -81,6 +81,8 @@ module SIS_types real, allocatable, dimension(:,:) :: & snow_to_ocn, & !< The mass per unit ocean area of snow that will be dumped into the !! ocean due to recent mechanical activities like ridging or drifting [R Z ~> kg m-2]. + water_to_ocn, & !< The mass per unit ocean area of pond water that will be dumped into the + !! ocean due to recent mechanical activities like ridging or drifting [R Z ~> kg m-2]. enth_snow_to_ocn !< The average enthalpy of the snow that will be dumped into the !! ocean due to recent mechanical activities like ridging or drifting [Q ~> J kg-1]. @@ -90,7 +92,7 @@ module SIS_types !! in each category and fractional thickness layer [Q ~> J kg-1]. real, allocatable, dimension(:,:,:,:) :: enth_snow !< The enthalpy of the snow !! in each category and snow thickness layer [Q ~> J kg-1]. - + 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]. @@ -452,7 +454,9 @@ subroutine alloc_IST_arrays(HI, IG, IST, omit_velocities, omit_Tsurf, do_ridging if (present(do_ridging)) then ; if (do_ridging) then allocate(IST%snow_to_ocn(isd:ied, jsd:jed)) ; IST%snow_to_ocn(:,:) = 0.0 + allocate(IST%water_to_ocn(isd:ied, jsd:jed)) ; IST%water_to_ocn(:,:) = 0.0 allocate(IST%enth_snow_to_ocn(isd:ied, jsd:jed)) ; IST%enth_snow_to_ocn(:,:) = 0.0 + allocate(IST%rdg_rate(isd:ied, jsd:jed)) ; IST%rdg_rate(:,:) = 0.0 allocate(IST%rdg_mice(isd:ied, jsd:jed, CatIce)) ; IST%rdg_mice(:,:,:) = 0.0 endif ; endif diff --git a/src/ice_model.F90 b/src/ice_model.F90 index e50f99792d..1b202b9acb 100644 --- a/src/ice_model.F90 +++ b/src/ice_model.F90 @@ -57,6 +57,7 @@ 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 @@ -2133,6 +2134,8 @@ 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, & diff --git a/src/ice_ridge.F90 b/src/ice_ridge.F90 index da75fea6aa..8b26e5299a 100644 --- a/src/ice_ridge.F90 +++ b/src/ice_ridge.F90 @@ -1,9 +1,19 @@ -!> A partially implmented ice ridging parameterizations -!! that does not yet work with an arbitrary number of vertical layers in the ice +!> A full implementation of Icepack ridging parameterizations module ice_ridging_mod ! This file is a part of SIS2. See LICENSE.md for the license. +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! +! replaced T. Martin code with wrapper for Icepack ridging function - mw 1/18 ! +! ! +! Prioritized to do list as of 6/4/19 (mw): ! +! ! +! 1) implement new snow_to_ocean diagnostic to record this flux. ! +! 2) implement ridging_rate diagnostics: ridging_shear, ridging_conv ! +! 3) implement "do_j" style optimization as in "compress_ice" or ! +! "adjust_ice_categories" (SIS_transport.F90) if deemed necessary ! +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! + use SIS_diag_mediator, only : post_SIS_data, query_SIS_averaging_enabled, SIS_diag_ctrl use SIS_diag_mediator, only : register_diag_field=>register_SIS_diag_field, time_type use MOM_domains, only : pass_var, pass_vector, BGRID_NE @@ -11,11 +21,26 @@ module ice_ridging_mod use MOM_file_parser, only : get_param, log_param, read_param, log_version, param_file_type use MOM_unit_scaling, only : unit_scale_type use SIS_hor_grid, only : SIS_hor_grid_type - +use SIS_types, only : ice_state_type, ist_chksum +use fms_io_mod, only : register_restart_field, restart_file_type +use SIS_tracer_registry, only : SIS_tracer_registry_type, SIS_tracer_type, get_SIS_tracer_pointer +use SIS2_ice_thm, only : get_SIS2_thermo_coefs +use ice_grid, only : ice_grid_type +!Icepack modules +use icepack_kinds +use icepack_itd, only: icepack_init_itd, cleanup_itd +use icepack_mechred, only: ridge_ice +use icepack_warnings, only: icepack_warnings_flush, icepack_warnings_aborted, & + icepack_warnings_setabort +use icepack_tracers, only: icepack_init_tracer_indices, icepack_init_tracer_sizes +use icepack_tracers, only: icepack_query_tracer_sizes +use icepack_parameters, only : icepack_init_parameters implicit none ; private #include + + public :: ice_ridging, ridge_rate, ice_ridging_init real, parameter :: hlim_unlim = 1.e8 !< Arbitrary huge number used in ice_ridging @@ -23,154 +48,17 @@ module ice_ridging_mod logical :: rdg_lipscomb = .true. !< If true, use the Lipscomb ridging scheme !! TODO: These parameters belong in a control structure -contains +! future namelist parameters? +integer (kind=int_kind), parameter :: & + krdg_partic = 0 , & ! 1 = new participation, 0 = Thorndike et al 75 + krdg_redist = 0 ! 1 = new redistribution, 0 = Hibler 80 -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! -! ice_ridging_init - initialize the ice ridging and set parameters. ! -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! +! e-folding scale of ridged ice, krdg_partic=1 (m^0.5) +real(kind=dbl_kind), parameter :: mu_rdg = 3.0 -!~~~~T. Martin, 2008~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! -!> ice_ridging_init makes some preparations for the ridging routine and is -!! called from within the subroutine ridging -subroutine ice_ridging_init(km, cn, hi, rho_ice, part_undef, part_undef_sum, & - hmin, hmax, efold, rdg_ratio, US) - integer, intent(in) :: km !< The number of ice thickness categories - real, dimension(0:km), intent(in) :: cn !< Fractional concentration of each thickness category, - !! including open water fraction [nondim] - real, dimension(km), intent(in) :: hi !< ice mass in each category [R Z ~> kg m-2] - real, intent(in) :: rho_ice !< Nominal ice density [R ~> kg m-3] - real, dimension(0:km), intent(out) :: part_undef !< fraction of undeformed ice or open water - !! participating in ridging [nondim] - real, intent(out) :: part_undef_sum !< The sum across categories of part_undef [nondim] - real, dimension(km), intent(out) :: hmin !< minimum ice thickness [R Z ~> kg m-2] involved in - !! Hibler's ridged ice distribution - real, dimension(km), intent(out) :: hmax !< maximum ice thickness [R Z ~> kg m-2] involved in - !! Hibler's ridged ice distribution - real, dimension(km), intent(out) :: efold !< e-folding scale lambda of Lipscomb's ridged ice - !! distribution [R Z ~> kg m-2] - real, dimension(km), intent(out) :: rdg_ratio !< ratio of ridged ice to total ice [nondim]; k_n in CICE - type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors - ! Local variables - real :: raft_max ! maximum weight [R Z ~> kg m-2] of ice that rafts - real, dimension(0:km) :: ccn ! cumulative ice thickness distribution [nondim] (or G in CICE documentation) - ! now set in namelist (see ice_dyn_param): - !Niki: The following two were commented out - real :: part_par ! participation function shape parameter in parent ice categories, - ! G* or a* in CICE [nondim] - real :: dist_par ! distribution function shape parameter in receiving categories, - ! H* or mu in CICE [R Z ~> kg m-2] - real :: dist_par_1 ! distribution function shape parameter in receiving categories, - ! H* or mu in CICE [R Z ~> kg m-2] - real :: part_pari ! [nondim] - real, dimension(0:km) :: rdgtmp ! [nondim] - integer :: k - logical :: rdg_lipscomb = .true. - - raft_max = 1.0*US%m_to_Z * rho_ice - - ! cumulative ice thickness distribution - ccn(0) = 0.0 ! helps to include calculations for open water in k-loops - do k=1,km - if (cn(k) > 1.e-10) then - ccn(k) = ccn(k-1)+cn(k) - else - ccn(k) = ccn(k-1) - endif - enddo - ! normalize ccn - do k=1,km ; ccn(k) = ccn(k) / ccn(km) ; enddo - - !---------------------------------------------------------------------------------------- - ! i) participation function for undeformed, thin ice: - ! amount of ice per thin ice category participating in ridging - ! ii) parameters defining the range of the thick ice categories - ! to which the newly ridged ice is redistributed - ! - ! A) scheme following Lipscomb et al., 2007, JGR - ! i) negative exponential participation function for level ice - ! ii) negative exponential distribution of newly ridged ice - ! - ! B) scheme following Thorndike et al., 1975, JGR - ! i) linear participation function with negative slope and root at (part_par,0.0) - ! and Hibler, 1980, Mon.Wea.Rev. - ! ii) uniform distribution of newly ridged ice in receiving categories - !---------------------------------------------------------------------------------------- - if (rdg_lipscomb) then - ! ************ - ! * Ai * - ! ************ - part_par = -0.05 ! this is -a* for practical reasons, - ! part_par(lipscomb)=part_par(thorn-hib)/3 for best comparability of schemes - !do k=1,km - !part_undef(k) = (exp(ccn(k-1)/part_par)-exp(ccn(k)/part_par)) / (1.-exp(1./part_par)) - !enddo - ! set in namelist: part_par = -20. ! this is -1/a*; CICE standard is a* = 0.05 - part_pari = 1. / (1.-exp(part_par)) - do k=0,km ; rdgtmp(k) = exp(ccn(k)*part_par) * part_pari ; enddo - do k=1,km ; part_undef(k) = rdgtmp(k-1) - rdgtmp(k) ; enddo - ! ************ - ! * Aii * - ! ************ - dist_par_1 = 4.0**2*US%m_to_Z * rho_ice - ! set in namelist: dist_par = 4.0 ! unit [m**0.5], e-folding scale of ridged ice, - ! for comparable results of Lipscomb and Thorn-Hib schemes choose - ! 3 & 25, 4 & 50, 5 & 75 or 6 & 100 - do k=2,km - if (hi(k)>0.0) then - efold(k) = sqrt(dist_par_1 * hi(k)) - hmin(k) = min(2.*hi(k), hi(k)+raft_max) - rdg_ratio(k) = (hmin(k)+efold(k)) / hi(k) - else - efold(k)=0.0 ; hmin(k)=0.0 ; rdg_ratio(k)=1.0 - endif - enddo - !---------------------------------------------------------------------------------------- - else - ! ************ - ! * Bi * - ! ************ - part_par = 0.15 - ! set in namelist: part_par = 0.15 ! CICE standard is 0.15 - do k=1,km - if (ccn(k) < part_par) then - part_undef(k) = (ccn(k) -ccn(k-1)) * (2.-( (ccn(k-1)+ccn(k)) /part_par )) / part_par - else if (ccn(k) >= part_par .and. ccn(k-1) < part_par) then - part_undef(k) = (part_par-ccn(k-1)) * (2.-( (ccn(k-1)+part_par)/part_par )) / part_par - else - part_undef(k) = 0.0 - endif - enddo - ! ************ - ! * Bii * - ! ************ - ! set in namelist: dist_par = 50.0 ! 25m suggested by CICE and is appropriate for multi-year ridges - ! (Flato and Hibler, 1995, JGR), - ! 50m gives better fit to first-year ridges (Amundrud et al., 2004, JGR) - dist_par = 50.0*US%m_to_Z * rho_ice - do k=2,km - if (hi(k)>0.0) then - ! minimum and maximum defining range in ice thickness space - ! that receives ice in the ridging process - hmin(k) = min(2.*hi(k), hi(k) + raft_max) - hmax(k) = max(2.*sqrt(dist_par*hi(k)), hmin(k) + 10.e-10*US%m_to_Z*rho_ice) - ! ratio of mean ridge thickness to thickness of parent level ice - rdg_ratio(k) = 0.5*(hmin(k) + hmax(k))/hi(k) - else - hmin(k)=0.0 ; hmax(k)=0.0 ; rdg_ratio(k)=1.0 - endif - enddo - - endif - !---------------------------------------------------------------------------------------- - - ! ratio of net ice area removed / total area participating - part_undef_sum = ccn(1) - do k=2,km - if (rdg_ratio(k)>0.0) part_undef_sum = part_undef_sum + part_undef(k) * (1. - 1./rdg_ratio(k)) - enddo +contains -end subroutine ice_ridging_init !TOM>~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! @@ -197,312 +85,398 @@ function ridge_rate(del2, div) result (rnet) ! dissipated in the ridging process end function ridge_rate -!TOM>~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! -!> ice_ridging parameterizes mechanical redistribution of thin (undeformed) ice -!! into thicker (deformed/ridged) ice categories -subroutine ice_ridging(km, cn, hi, hs, rho_ice, t1, t2, age, snow_to_ocn, enth_snow_to_ocn, rdg_rate, hi_rdg, & - dt, hlim_in, rdg_open, vlev, US) - ! Subroutine written by T. Martin, 2008 - integer, intent(in) :: km !< The number of ice thickness categories - real, dimension(0:km), intent(inout) :: cn !< Fractional concentration of each thickness category, - !! including open water fraction [nondim] - real, dimension(km), intent(inout) :: hi !< ice mass in each category [R Z ~> kg m-2] - real, dimension(km), intent(inout) :: hs !< snow mass in each category [R Z ~> kg m-2] - real, intent(in) :: rho_ice !< Nominal ice density [R ~> kg m-3] - ! CAUTION: these quantities are extensive here, - real, dimension(km), intent(inout) :: t1 !< Volume integrated upper layer temperature [degC m3]? - real, dimension(km), intent(inout) :: t2 !< Volume integrated upper layer temperature [degC m3]? - real, dimension(km), intent(inout) :: age !< Volume integrated ice age [m3 years]? - real, intent(out) :: enth_snow_to_ocn !< average of enthalpy of the snow dumped into - !! ocean due to this ridging event [Q ~> J kg-1] - real, intent(out) :: snow_to_ocn !< total snow mass dumped into ocean due to this - !! ridging event [R Z ~> kg m-2] - real, intent(in) :: rdg_rate !< Ridging rate from subroutine ridge_rate [T-1 ~> s-1] - real, dimension(km), intent(inout) :: hi_rdg !< A diagnostic of the ridged ice volume in each - !! category [R Z ~> kg m-2]. - real, intent(in) :: dt !< time step [T ~> s] - real, dimension(km), intent(in) :: hlim_in !< ice thickness category limits - real, intent(out) :: rdg_open !< Rate of change in open water area due to - !! newly formed ridges [T-1 ~> s-1] - real, intent(out) :: vlev !< mass of level ice participating in ridging [R Z T-1 ~> kg m-2 s-1] - type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors - ! Local variables - integer :: k, kd, kr, n_iterate - integer, parameter :: n_itermax = 10 ! maximum number of iterations for redistribution -! real, parameter :: frac_hs_rdg = 0.5 ! fraction of snow that remains on ridged ice; - real :: frac_hs_rdg ! fraction of snow that remains on ridged ice [nondim]; - ! (1.0-frac_hs_rdg)*hs falls into ocean - real, dimension(0:km) :: part_undef ! fraction of undeformed ice or open water participating in ridging - real :: area_undef ! fractional area of parent and ... - real, dimension(km) :: area_def ! ... newly ridged ice, respectively [nondim] - real, dimension(km) :: vol_def ! fractional volume of newly ridged ice [nondim] - real, dimension(km) :: cn_old ! concentrations at beginning of iteration loop [nondim] - real, dimension(km) :: hi_old ! thicknesses at beginning of iteration loop [R Z ~> kg m-2] - real, dimension(km) :: rdg_frac ! ratio of ridged and total ice volume [nondim] - real :: alev ! area of level ice participating in ridging [nondim] - real :: ardg, vrdg ! area and volume of newly formed rdiged (vlev=vrdg!!!) - real, dimension(km) :: hmin, hmax, efold ! [R Z ~> kg m-2] - real, dimension(km) :: rdg_ratio ! [nondim] - real, dimension(km) :: hlim ! [R Z ~> kg m-2] - real :: hl, hr - real :: snow_dump, enth_dump - real :: cn_tot, part_undef_sum - real :: div_adv, Rnet, Rdiv, Rtot ! [T-1 ~> s-1] - real :: rdg_area, rdgtmp, hlimtmp - real :: area_frac - real, dimension(km) :: area_rdg ! [nondim] - real, dimension(km) :: & - frac_hi, frac_hs, & ! Portion of the ice and snow that are being ridged [R Z ~> kg m-2] - frac_t1, frac_t2, & - frac_age - logical :: rdg_iterate - !------------------------------------------------------------------- - ! some preparations - !------------------------------------------------------------------- - hlimtmp = hlim_in(km) - hlim(km) = hlim_unlim ! ensures all ridged ice is smaller than thickest ice allowed - frac_hs_rdg = 1.0 - s2o_frac - snow_to_ocn = 0.0 ; enth_snow_to_ocn = 0.0 - alev=0.0 ; ardg=0.0 ; vlev=0.0 ; vrdg=0.0 +subroutine ice_ridging_init(G,IG,TrReg,US) + type(SIS_hor_grid_type), intent(inout) :: G !< G The ocean's grid structure. + type(ice_grid_type), intent(inout) :: IG !< The sea-ice-specific grid structure. + type(SIS_tracer_registry_type), pointer :: TrReg ! TrReg - The registry of registered SIS ice and snow tracers. + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors. + + integer (kind=int_kind) :: ntrcr, ncat, nilyr, nslyr, nblyr, nfsd, n_iso, n_aero + integer (kind=int_kind) :: nt_Tsfc, nt_sice, nt_qice, nt_vlvl, nt_qsno + + ncat=IG%CatIce ! The number of sea-ice thickness categories + nilyr=IG%NkIce ! The number of ice layers per category + nslyr=IG%NkSnow ! The number if snow layers per category + nblyr=0 ! The number of bio/brine layers per category + nfsd=0 ! The number of floe size distribution layers + n_iso=0 ! The number of isotopes in use + n_aero=0 ! The number of aerosols in use + nt_Tsfc=1 ! Tracer index for ice/snow surface temperatore + nt_qice=2 ! Starting index for ice enthalpy in layers + nt_qsno=2+nilyr ! Starting index for snow enthalpy + nt_sice=2+nilyr+nslyr ! Index for ice salinity + nt_vlvl=2+2*nilyr+nslyr ! Index for level ice volume fraction + ntrcr=2+2*nilyr+nslyr ! number of tracers in use + ! (1,2) snow/ice surface temperature + + ! (3) ice salinity*nilyr + (4) pond thickness + + call icepack_init_tracer_sizes(ntrcr_in=ntrcr, & + ncat_in=ncat, nilyr_in=nilyr, nslyr_in=nslyr, nblyr_in=nblyr, & + nfsd_in=nfsd, n_iso_in=n_iso, n_aero_in=n_aero) + + call icepack_init_tracer_indices(nt_Tsfc_in=nt_Tsfc, & + nt_sice_in=nt_sice, nt_qice_in=nt_qice, & + nt_qsno_in=nt_qsno, nt_vlvl_in=nt_vlvl ) + + call icepack_init_parameters(mu_rdg_in=mu_rdg,conserv_check_in=.true.) + +end subroutine ice_ridging_init +! +! ice_ridging is a wrapper for the icepack ridging routine ridge_ice +! +subroutine ice_ridging(IST, G, IG, mca_ice, mca_snow, mca_pond, TrReg, US, dt, rdg_rate) + type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice + type(SIS_hor_grid_type), intent(inout) :: G !< G The ocean's grid structure. + type(ice_grid_type), intent(inout) :: IG !< The sea-ice-specific grid structure. + real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), intent(inout) :: mca_ice, mca_snow, mca_pond + type(SIS_tracer_registry_type), pointer :: TrReg ! TrReg - The registry of registered SIS ice and snow tracers. + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors. + real (kind=dbl_kind), intent(in) :: dt !< The amount of time over which the ice dynamics are to be. + ! advanced in seconds. + real, dimension(SZI_(G),SZJ_(G)), intent(inout), optional :: rdg_rate !< Diagnostic of the rate of fractional + !! area loss-gain due to ridging (1/s) + +! logical, intent(in) :: dyn_Cgrid !< True if using C-gridd velocities, B-grid if False. + + + + ! these strain metrics are calculated here from the velocities used for advection + real :: sh_Dt ! sh_Dt is the horizontal tension (du/dx - dv/dy) including + ! all metric terms, in s-1. + real :: sh_Dd ! sh_Dd is the flow divergence (du/dx + dv/dy) including all + ! metric terms, in s-1. + real, dimension(SZIB_(G),SZJB_(G)) :: & + sh_Ds ! sh_Ds is the horizontal shearing strain (du/dy + dv/dx) + ! including all metric terms, in s-1. + + integer :: i, j, k ! loop vars + integer isc, iec, jsc, jec ! loop bounds + integer :: halo_sh_Ds ! The halo size that can be used in calculating sh_Ds. + + integer (kind=int_kind) :: & + ncat , & ! number of thickness categories + nilyr , & ! number of ice layers + nslyr ! number of snow layers + + + real (kind=dbl_kind), dimension(0:IG%CatIce) :: hin_max ! category limits (m) + + logical (kind=log_kind) :: & + closing_flag, &! flag if closing is valid + tr_brine ! if .true., brine height differs from ice thickness + + ! optional history fields + real (kind=dbl_kind) :: & + dardg1dt , & ! rate of fractional area loss by ridging ice (1/s) + dardg2dt , & ! rate of fractional area gain by new ridges (1/s) + dvirdgdt , & ! rate of ice volume ridged (m/s) + opening , & ! rate of opening due to divergence/shear (1/s) + closing , & ! rate of closing due to divergence/shear (1/s) + fpond , & ! fresh water flux to ponds (kg/m^2/s) + fresh , & ! fresh water flux to ocean (kg/m^2/s) + fhocn ! net heat flux to ocean (W/m^2) + + real (kind=dbl_kind), dimension(IG%CatIce) :: & + dardg1ndt , & ! rate of fractional area loss by ridging ice (1/s) + dardg2ndt , & ! rate of fractional area gain by new ridges (1/s) + dvirdgndt , & ! rate of ice volume ridged (m/s) + aparticn , & ! participation function + krdgn , & ! mean ridge thickness/thickness of ridging ice + araftn , & ! rafting ice area + vraftn , & ! rafting ice volume + aredistn , & ! redistribution function: fraction of new ridge area + vredistn ! redistribution function: fraction of new ridge volume + + real (kind=dbl_kind), dimension(IG%CatIce) :: & + faero_ocn ! aerosol flux to ocean (kg/m^2/s) + + real (kind=dbl_kind), dimension(IG%CatIce) :: & + fiso_ocn ! isotope flux to ocean (kg/m^2/s) + + integer (kind=int_kind) :: & + ndtd = 1 , & ! number of dynamics subcycles + n_aero = 0, & ! number of aerosol tracers + ntrcr = 0 ! number of tracer level + + + real(kind=dbl_kind) :: & + del_sh , & ! shear strain measure + rdg_conv = 0.0, & ! normalized energy dissipation from convergence (1/s) + rdg_shear= 0.0 ! normalized energy dissipation from shear (1/s) + + real(kind=dbl_kind), dimension(IG%CatIce) :: & + aicen, & ! concentration of ice + vicen, & ! volume per unit area of ice (m) + vsnon, & ! volume per unit area of snow (m) + tr_tmp ! for temporary storade + ! ice tracers; ntr*(NkIce+NkSnow) guaranteed to be enough for all (intensive) + real(kind=dbl_kind), dimension(2+2*IG%NkIce+IG%NkSnow,IG%CatIce) :: trcrn + + real(kind=dbl_kind) :: aice0 ! concentration of open water + + integer (kind=int_kind), dimension(2+2*IG%NkIce+IG%NkSnow) :: & + trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon (weighting to use) + n_trcr_strata ! number of underlying tracer layers + + real(kind=dbl_kind), dimension(2+2*IG%NkIce+IG%NkSnow,3) :: & + trcr_base ! = 0 or 1 depending on tracer dependency + ! argument 2: (1) aice, (2) vice, (3) vsno + + integer(kind=int_kind), dimension(2+2*IG%NkIce+IG%NkSnow,IG%CatIce) :: & + nt_strata ! indices of underlying tracer layers + + type(SIS_tracer_type), dimension(:), pointer :: Tr=>NULL() ! SIS2 tracers + real, dimension(:,:,:,:), pointer :: Tr_ice_enth_ptr=>NULL() !< A pointer to the named tracer + real, dimension(:,:,:,:), pointer :: Tr_snow_enth_ptr=>NULL() !< A pointer to the named tracer + real, dimension(:,:,:,:), pointer :: Tr_ice_salin_ptr=>NULL() !< A pointer to the named tracer + + real :: rho_ice, rho_snow, divu_adv + integer :: m, n ! loop vars for tracer; n is tracer #; m is tracer layer + integer(kind=int_kind) :: nt_tsfc_in, nt_qice_in, nt_qsno_in, nt_sice_in + integer(kind=int_kind) :: nL_ice, nL_snow ! number of tracer levels + integer(kind=int_kind) :: ncat_out, ntrcr_out, nilyr_out, nslyr_out ! array sizes returned from Icepack query + + nSlyr = IG%NkSnow + nIlyr = IG%NkIce + nCat = IG%CatIce + isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec + + call get_SIS2_thermo_coefs(IST%ITV, rho_ice=rho_ice) + call get_SIS2_thermo_coefs(IST%ITV, rho_snow=rho_snow) + + call icepack_query_tracer_sizes(ncat_out=ncat_out,ntrcr_out=ntrcr_out, nilyr_out=nilyr_out, nslyr_out=nslyr_out) + + + if (nIlyr .ne. nilyr_out .or. nSlyr .ne. nslyr_out ) call SIS_error(FATAL,'nilyr or nslyr mismatch with Icepack') + + ! copy strain calculation code from SIS_C_dynamics; might be a more elegant way ... ! - call ice_ridging_init(km, cn, hi, rho_ice, part_undef, part_undef_sum, & - hmin, hmax, efold, rdg_ratio, US) - - !------------------------------------------------------------------- - ! opening and closing rates of the ice cover - !------------------------------------------------------------------- - - ! update total area fraction as this may exceed 1 after transportation/advection - ! (cn_tot <= 1 after ridging!) - cn_tot = sum(cn(0:km)) - - ! dissipated energy in ridging from state of ice drift - ! after Flato and Hibler (1995, JGR) - ! (see subroutine ridge_rate in ice_dyn_mod), - ! equals net closing rate times ice strength - ! by passing to new, local variable rdg_rate is saved for diagnostic output - Rnet = rdg_rate - ! the divergence rate given by the advection scheme ... - div_adv = (1.-cn_tot) / dt - ! ... may exceed the rate derived from the drift state (Rnet) - if (div_adv < 0.) Rnet = max(Rnet, -div_adv) - ! non-negative opening rate that ensures cn_tot <=1 after ridging - Rdiv = Rnet + div_adv - ! total closing rate - Rtot = Rnet / part_undef_sum - - !------------------------------------------------------------------- - ! iteration of ridging redistribution - do n_iterate=1, n_itermax - !------------------------------------------------------------------- - - ! save initial state of ice concentration, total and ridged ice volume - ! at beginning of each iteration loop - do k=1,km - cn_old(k) = cn(k) - hi_old(k) = hi(k) - - rdg_frac(k) = 0.0 ; if (hi(k)>0.0) rdg_frac(k) = hi_rdg(k)/hi(k) - enddo - - ! reduce rates in case more than 100% of any category would be removed - do k=1,km ; if (cn(k)>1.0e-10 .and. part_undef(k)>0.0) then - rdg_area = part_undef(k) * Rtot * dt ! area ridged in category k - if (rdg_area > cn(k)) then - rdgtmp = cn(k)/rdg_area - Rtot = Rtot * rdgtmp - Rdiv = Rdiv * rdgtmp + halo_sh_Ds = min(isc-G%isd, jsc-G%jsd, 2) + ! if (dyn_Cgrid) then + do J=jsc-halo_sh_Ds,jec+halo_sh_Ds-1 ; do I=isc-halo_sh_Ds,iec+halo_sh_Ds-1 + ! This uses a no-slip boundary condition. + sh_Ds(I,J) = (2.0-G%mask2dBu(I,J)) * & + (G%dxBu(I,J)*G%IdyBu(I,J)*(IST%u_ice_C(I,j+1)*G%IdxCu(I,j+1) - IST%u_ice_C(I,j)*G%IdxCu(I,j)) + & + G%dyBu(I,J)*G%IdxBu(I,J)*(IST%v_ice_C(i+1,J)*G%IdyCv(i+1,J) - IST%v_ice_C(i,J)*G%IdyCv(i,J))) + enddo; enddo + + ! set category limits; Icepack has a max on the largest, unlimited, category (why?) + + hin_max(0)=0.0 + do k=1,nCat + hin_max(k) = IG%mH_cat_bound(k)/(Rho_ice*IG%kg_m2_to_H) + end do + + trcr_base = 0.0; n_trcr_strata = 0; nt_strata = 0; ! init some tracer vars + ! When would we use icepack tracer "strata"? + + ! set icepack tracer index "nt_lvl" to (last) pond tracer so it gets dumped when + ! ridging in ridge_ice (this is what happens to "level" ponds); first add up ntrcr; + ! then set nt_lvl to ntrcr+1; could move this to an initializer - mw + + call get_SIS_tracer_pointer("enth_ice",TrReg,Tr_ice_enth_ptr,nL_ice) + call get_SIS_tracer_pointer("enth_snow",TrReg,Tr_snow_enth_ptr,nL_snow) + call get_SIS_tracer_pointer("salin_ice",TrReg,Tr_ice_salin_ptr,nL_ice) + +! call IST_chksum('before ice ridging ',IST,G,US,IG) + + if (present(rdg_rate)) rdg_rate(:,:)=0.0 + do j=jsc,jec; do i=isc,iec + if ((G%mask2dT(i,j) .gt. 0.0) .and. (sum(IST%part_size(i,j,1:nCat)) .gt. 0.0)) then + ! feed locations to Icepack's ridge_ice + + ! start like we're putting ALL the snow and pond in the ocean + IST%snow_to_ocn(i,j) = IST%snow_to_ocn(i,j) + sum(mca_snow(i,j,:)) + IST%enth_snow_to_ocn(i,j) = IST%enth_snow_to_ocn(i,j) + sum(mca_snow(i,j,:)*TrReg%Tr_snow(1)%t(i,j,:,1)); + IST%water_to_ocn(i,j) = IST%water_to_ocn(i,j) + sum(mca_pond(i,j,:)); + aicen(1:nCat) = IST%part_size(i,j,1:nCat); + + if (sum(aicen) .eq. 0.0) then ! no ice -> no ridging + IST%part_size(i,j,0) = 1.0 + else + ! set up ice and snow volumes + vicen(1:nCat) = mca_ice(i,j,1:nCat) /Rho_ice ! volume per unit area of ice (m) + vsnon(1:nCat) = mca_snow(i,j,1:nCat)/Rho_snow ! volume per unit area of snow (m) + + sh_Dt = (G%dyT(i,j)*G%IdxT(i,j)*(G%IdyCu(I,j) * IST%u_ice_C(I,j) - & + G%IdyCu(I-1,j)*IST%u_ice_C(I-1,j)) - & + G%dxT(i,j)*G%IdyT(i,j)*(G%IdxCv(i,J) * IST%v_ice_C(i,J) - & + G%IdxCv(i,J-1)*IST%v_ice_C(i,J-1))) + sh_Dd = (G%IareaT(i,j)*(G%dyCu(I,j) * IST%u_ice_C(I,j) - & + G%dyCu(I-1,j)*IST%u_ice_C(I-1,j)) + & + G%IareaT(i,j)*(G%dxCv(i,J) * IST%v_ice_C(i,J) - & + G%dxCv(i,J-1)*IST%v_ice_C(i,J-1))) + + del_sh = sqrt(sh_Dd**2 + 0.25 * (sh_Dt**2 + & + (0.25 * ((sh_Ds(I-1,J-1) + sh_Ds(I,J)) + & + (sh_Ds(I-1,J) + sh_Ds(I,J-1))))**2 ) ) ! H&D eqn 9 + rdg_conv = -min(sh_Dd,0.0) ! energy dissipated by convergence ... + rdg_shear = 0.5*(del_sh-abs(sh_Dd)) ! ... and by shear + aice0 = IST%part_size(i,j,0) + ntrcr = 0 +! Tr_ptr=>NULL() + if (TrReg%ntr>0) then ! load tracer array + ntrcr=ntrcr+1 + trcrn(ntrcr,1) = Tr_ice_enth_ptr(i,j,1,1) ! surface temperature? this is taken from the thinnest ice category + trcr_depend(ntrcr) = 1; ! 1 means ice-based tracer + trcr_base(ntrcr,:) = 0.0; trcr_base(ntrcr,2) = 1.0; ! 2nd index for ice + do k=1,nL_ice + ntrcr=ntrcr+1 + trcrn(ntrcr,1:ncat)=Tr_ice_enth_ptr(i,j,1:nCat,k) + trcr_depend(ntrcr) = 1; ! 1 means ice-based tracer + trcr_base(ntrcr,:) = 0.0; trcr_base(ntrcr,2) = 1.0; ! 2nd index for ice + enddo + do k=1,nL_snow + ntrcr=ntrcr+1 + trcrn(ntrcr,1:nCat)=Tr_snow_enth_ptr(i,j,1:nCat,k) + trcr_depend(ntrcr) = 2; ! 2 means snow-based tracer + trcr_base(ntrcr,:) = 0.0; trcr_base(ntrcr,3) = 1.0; ! 3rd index for snow + enddo + do k=1,nL_ice + ntrcr=ntrcr+1 + trcrn(ntrcr,1:nCat)=Tr_ice_salin_ptr(i,j,1:nCat,k) + trcr_depend(ntrcr) = 1; ! 1 means ice-based tracer + trcr_base(ntrcr,:) = 0.0; trcr_base(ntrcr,2) = 1.0; ! 2nd index for ice + enddo + endif ! have tracers to load + + ! load pond on top of stack + ntrcr = ntrcr + 1 + trcrn(ntrcr,1:nCat) = IST%mH_pond(i,j,1:nCat) + trcr_depend(ntrcr) = 0; ! 0 means ice area-based tracer + trcr_base(ntrcr,:) = 0.0; trcr_base(ntrcr,1) = 1.0; ! 1st index for ice area + + if (ntrcr .ne. ntrcr_out ) call SIS_error(FATAL,'ntrcr mismatch with Icepack') + + tr_brine = .false. + dardg1dt=0.0 + dardg2dt=0.0 + opening=0.0 + fpond=0.0 + fresh=0.0 + fhocn=0.0 + faero_ocn=0.0 + fiso_ocn=0.0 + aparticn=0.0 + krdgn=0.0 + aredistn=0.0 + vredistn=0.0 + dardg1ndt=0.0 + dardg2ndt=0.0 + dvirdgndt=0.0 + araftn=0.0 + vraftn=0.0 + closing_flag=.false. + + ! call Icepack routine; how are ponds treated? + call ridge_ice (dt, ndtd, & + ncat, n_aero, & + nilyr, nslyr, & + ntrcr, hin_max, & + rdg_conv, rdg_shear, & + aicen, & + trcrn, & + vicen, vsnon, & + aice0, & + trcr_depend, & + trcr_base, & + n_trcr_strata, & + nt_strata, & + krdg_partic, krdg_redist, & + mu_rdg, tr_brine, & + dardg1dt=dardg1dt, dardg2dt=dardg2dt, & + dvirdgdt=dvirdgdt, opening=opening, & + fpond=fpond, & + fresh=fresh, fhocn=fhocn, & + faero_ocn=faero_ocn, fiso_ocn=fiso_ocn, & + aparticn=aparticn, & + krdgn=krdgn, & + aredistn=aredistn, & + vredistn=vredistn, & + dardg1ndt=dardg1ndt, & + dardg2ndt=dardg2ndt, & + dvirdgndt=dvirdgndt, & + araftn=araftn, & + vraftn=vraftn, & + closing_flag=closing_flag ,closing=closing) + + rdg_rate(i,j) = dardg1dt-dardg2dt + + if ( icepack_warnings_aborted() ) then + call icepack_warnings_flush(0); + call icepack_warnings_setabort(.false.) + call SIS_error(WARNING,'icepack ridge_ice error'); endif - endif ; enddo - !------------------------------------------------------------------- - ! redistribution of ice - !------------------------------------------------------------------- + ! pop pond off top of stack + tr_tmp(1:nCat)=trcrn(ntrcr,1:nCat) - ! changes in open water area - cn(0) = max(cn(0) + (Rdiv - part_undef(1)*Rtot) * dt, 0.0) + do k=1,nCat + IST%mH_pond(i,j,k) = tr_tmp(k) + mca_pond(i,j,k) = IST%mH_pond(i,j,k)*aicen(k) + enddo - if (Rtot>0.0) then + if (TrReg%ntr>0) then + ! unload tracer array reversing order of load -- stack-like fashion - ! area, volume and energy changes in each category - do kd=1,km ! donating category - area_undef = min(part_undef(kd)*Rtot*dt, cn_old(kd)) ! area that experiences ridging in category k, - ! make sure that not more than 100% are used - if (cn_old(kd) > 1.0e-10) then - area_frac = area_undef / cn_old(kd) ! fraction of level ice area involved in ridging - area_rdg(kd) = area_undef / rdg_ratio(kd) ! area of new ridges in category k - else - area_frac = 0.0 - area_rdg(kd) = 0.0 - endif - !if (rdg_ratio(kd) > 0.0) then ! distinguish between level and ridged ice in - !else ! each category: let only level ice ridge; - !endif ! sea also change of hi_rdg below - - ! reduce area, volume and energy of snow and ice in source category - frac_hi(kd) = hi(kd) * area_frac - frac_hs(kd) = hs(kd) * area_frac - frac_t1(kd) = t1(kd) * area_frac - frac_t2(kd) = t2(kd) * area_frac - frac_age(kd) = age(kd) * area_frac - - cn(kd) = cn(kd) - area_undef - hi(kd) = hi(kd) - frac_hi(kd) - hs(kd) = hs(kd) - frac_hs(kd) - t1(kd) = t1(kd) - frac_t1(kd) - t2(kd) = t2(kd) - frac_t2(kd) - age(kd) = age(kd) - frac_age(kd) - - alev = alev + area_undef ! diagnosing area of level ice participating in ridging - vlev = vlev + frac_hi(kd) ! diagnosing total ice volume moved due to ridging - ! (here donating categories) - - ! Here it is assumed that level and ridged ice - ! of a category participate in ridging in equal - ! measure; this also means that ridged ice may be ridged again - hi_rdg(kd) = hi_rdg(kd) - rdg_frac(kd)*frac_hi(kd) - hi_rdg(kd) = max(hi_rdg(kd), 0.0) ! ensure hi_rdg >= 0 - - ! dump part of the snow in ocean (here, sum volume, transformed to flux in update_ice_model_slow) - snow_dump = frac_hs(kd)*(1.0-frac_hs_rdg) - if (snow_to_ocn > 0.0) then - enth_dump = t1(kd) !### THIS IS WRONG, BUT IS A PLACEHOLDER FOR NOW. - enth_snow_to_ocn = (enth_snow_to_ocn*snow_to_ocn + enth_dump*snow_dump) / (snow_to_ocn + snow_dump) - snow_to_ocn = snow_to_ocn + snow_dump - endif + do k=1,nL_ice + tr_tmp(1:nCat)=trcrn(1+k,1:nCat) + Tr_ice_enth_ptr(i,j,1:nCat,k) = tr_tmp(1:nCat) + enddo - enddo + do k=1,nL_snow + tr_tmp(1:nCat)=trcrn(1+nl_ice+k,1:ncat) + Tr_snow_enth_ptr(i,j,1:nCat,k) = tr_tmp(1:nCat) + enddo - ! split loop in order to derive frac_... variables with initial status (before ridging redistribution) - do kd=1,km - - !---------------------------------------------------------------------------------------- - ! add area, volume and energy in receiving category : - ! A) after Lipscomb, 2007 (negative exponential distribution) - ! B) after Hibler, 1980, Mon. Weather Rev. (uniform distribution) - !---------------------------------------------------------------------------------------- - if (rdg_lipscomb) then - ! ************ - ! * A * - ! ************ - if (efold(kd)>0.0) then - do kr=1,km-1 ! receiving categories - if (hmin(kd) >= hlim(kr)) then - area_def(kr) = 0.0 - vol_def(kr) = 0.0 - else - hl = max(hmin(kd), hlim(kr-1)) - hr = hlim(kr) - area_def(kr) = exp((hmin(kd)-hl)/efold(kd)) & - - exp((hmin(kd)-hr)/efold(kd)) - vol_def(kr) = ( (hl+efold(kd))*exp((hmin(kd)-hl)/efold(kd)) & - - (hr+efold(kd))*exp((hmin(kd)-hr)/efold(kd)) ) & - / (hmin(kd)+efold(kd)) - endif - enddo ! k receiving - ! thickest categery is a special case: - hl = max(hmin(kd), hlim(km-1)) - area_def(km) = exp((hmin(kd)-hl)/efold(kd)) - vol_def(km) = ( (hl+efold(kd))*exp((hmin(kd)-hl)/efold(kd)) ) & - / (hmin(kd)+efold(kd)) - else - do kr=1,km - area_def(kr) = 0.0 - vol_def(kr) = 0.0 - enddo - endif - !---------------------------------------------------------------------------------------- - else ! not rdg_lipscomb - ! ************ - ! * B * - ! ************ - if (hmax(kd)==hmin(kd)) then - do kr=1,km ; area_def(kr) = 0.0 ; vol_def(kr) = 0.0 ; enddo - else - do kr=1,km ! receiving categories - if (hmin(kd) >= hlim(kr) .or. hmax(kd) <= hlim(kr-1)) then - hl = 0.0 - hr = 0.0 - else - hl = max(hmin(kd), hlim(kr-1)) - hr = min(hmax(kd), hlim(kr) ) - endif - area_def(kr) = (hr -hl ) / (hmax(kd) -hmin(kd) ) - !vol_def(kr) = (hr**2-hl**2) / (hmax(kd)**2-hmin(kd)**2) - vol_def(kr) = area_def(kr) * (hr+hl) / (hmax(kd)+hmin(kd)) - enddo ! k receiving - endif - - endif - !---------------------------------------------------------------------------------------- - - ! update ice/snow area, volume, energy for receiving categories - do kr=1,km ! receiving categories - cn(kr) = cn(kr) + area_def(kr) * area_rdg(kd) - rdgtmp = vol_def(kr) * frac_hi(kd) - hi(kr) = hi(kr) + rdgtmp - hs(kr) = hs(kr) + vol_def(kr) * frac_hs(kd) * frac_hs_rdg - t1(kr) = t1(kr) + vol_def(kr) * frac_t1(kd) - t2(kr) = t2(kr) + vol_def(kr) * frac_t2(kd) - age(kr) = age(kr) + vol_def(kr) * frac_age(kd) - - ardg = ardg + area_def(kr) * area_rdg(kd) ! diagnosing area of newly ridged ice - vrdg = vrdg + rdgtmp ! diagnosing total ice volume moved due to ridging - ! (here receiving categories, cross check with vlev) - - ! add newly ridged ice volume to total ridged ice in each category - hi_rdg(kr) = hi_rdg(kr) + rdgtmp - enddo + do k=1,nL_ice + tr_tmp(1:nCat)=trcrn(1+k+nl_ice+nl_snow,1:nCat) + Tr_ice_salin_ptr(i,j,1:nCat,k) = tr_tmp(1:nCat) + enddo - enddo ! kd loop over donating categories - - endif ! Rtot>0.0 - - ! update total area fraction and check if this is now <= 1 - ! and rerun the ice redistribution when necessary - cn_tot = sum(cn(0:km)) - rdg_iterate = .false. - if (abs(cn_tot-1.) > 1.0e-10) then - rdg_iterate = .true. - div_adv = (1.-cn_tot) / dt - Rnet = max(0.0, -div_adv) - Rdiv = max(0.0, div_adv) - call ice_ridging_init(km, cn, hi, rho_ice, part_undef, part_undef_sum, & - hmin, hmax, efold, rdg_ratio, US) - Rtot = Rnet / part_undef_sum - endif + endif ! have tracers to unload + + ! ! output: snow/ice masses/thicknesses + do k=1,nCat - !------------------------------------------------------------------- - if (.not. rdg_iterate) exit - enddo ! n_iterate - !------------------------------------------------------------------- + if (aicen(k) > 0.0) then + IST%part_size(i,j,k) = aicen(k) + mca_ice(i,j,k) = vicen(k)*Rho_ice + IST%mH_ice(i,j,k) = vicen(k)*Rho_ice/aicen(k) + mca_snow(i,j,k) = vsnon(k)*Rho_snow + IST%mH_snow(i,j,k) = vsnon(k)*Rho_snow/aicen(k) + else + IST%part_size(i,j,k) = 0.0 + mca_ice(i,j,k) = 0.0 + IST%mH_ice(i,j,k) = 0.0 + mca_snow(i,j,k) = 0.0 + IST%mH_snow(i,j,k) = 0.0 + endif + + enddo - ! check ridged ice volume for natural limits - do k=1,km - hi_rdg(k) = max(hi_rdg(k),0.0) ! ridged ice volume positive - hi_rdg(k) = min(hi_rdg(k),hi(k)) ! ridged ice volume not to exceed total ice volume - enddo + IST%part_size(i,j,0) = 1.0 - sum(IST%part_size(i,j,1:nCat)) - ! calculate opening rate of ridging - rdg_open = (alev - ardg) / dt + endif + ! subtract new snow/pond mass and energy on ice to sum net fluxes to ocean + IST%snow_to_ocn(i,j) = IST%snow_to_ocn(i,j) - sum(mca_snow(i,j,:)); + IST%enth_snow_to_ocn(i,j) = IST%enth_snow_to_ocn(i,j) - sum(mca_snow(i,j,:)*TrReg%Tr_snow(1)%t(i,j,:,1)); + IST%water_to_ocn(i,j) = IST%water_to_ocn(i,j) - sum(mca_pond(i,j,:)); - ! cross check ice volume transferred from level to ridged ice - if (abs(vlev - vrdg) > 1.0e-10*US%kg_m3_to_R*US%m_to_Z) then - print *,'WARNING: subroutine ice_ridging: parent ice volume does not equal ridge volume', vlev, vrdg - endif - ! turn vlev into a rate for diagnostics - vlev = vlev / dt + endif; enddo; enddo ! part_sz, j, i - ! return to true upper most ice thickness category limit - !hlim(km) = hlimtmp +! call IST_chksum('after ice ridging ',IST,G,US,IG) end subroutine ice_ridging + +! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> ice_ridging_end deallocates the memory associated with this module. subroutine ice_ridging_end() end subroutine ice_ridging_end + end module ice_ridging_mod From b528a23fc506eea7e01a0a1d0147259c4a450941 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Mon, 15 Mar 2021 17:58:08 -0400 Subject: [PATCH 3/3] Re-normalized ice coverage after call to ridge_ice - This appears to be needed to avoid propagation of negative ice free categories. Residual ice excess values of order 10-5 are being renormalized. --- src/ice_ridge.F90 | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/ice_ridge.F90 b/src/ice_ridge.F90 index 8b26e5299a..ceaafc2b92 100644 --- a/src/ice_ridge.F90 +++ b/src/ice_ridge.F90 @@ -292,8 +292,9 @@ subroutine ice_ridging(IST, G, IG, mca_ice, mca_snow, mca_pond, TrReg, US, dt, r IST%water_to_ocn(i,j) = IST%water_to_ocn(i,j) + sum(mca_pond(i,j,:)); aicen(1:nCat) = IST%part_size(i,j,1:nCat); + if (sum(aicen) .eq. 0.0) then ! no ice -> no ridging - IST%part_size(i,j,0) = 1.0 + IST%part_size(i,j,0) = 1.0; else ! set up ice and snow volumes vicen(1:nCat) = mca_ice(i,j,1:nCat) /Rho_ice ! volume per unit area of ice (m) @@ -313,7 +314,13 @@ subroutine ice_ridging(IST, G, IG, mca_ice, mca_snow, mca_pond, TrReg, US, dt, r (sh_Ds(I-1,J) + sh_Ds(I,J-1))))**2 ) ) ! H&D eqn 9 rdg_conv = -min(sh_Dd,0.0) ! energy dissipated by convergence ... rdg_shear = 0.5*(del_sh-abs(sh_Dd)) ! ... and by shear + aice0 = IST%part_size(i,j,0) + if (aice0<0.) then + call SIS_error(WARNING,'aice0<0. before call to ridge ice.') + aice0=0. + endif + ntrcr = 0 ! Tr_ptr=>NULL() if (TrReg%ntr>0) then ! load tracer array @@ -437,6 +444,11 @@ subroutine ice_ridging(IST, G, IG, mca_ice, mca_snow, mca_pond, TrReg, US, dt, r endif ! have tracers to unload + if (sum(aicen(1:nCat))>1.0) then + call SIS_error(WARNING,'Ice cover exceeds 1 after call to ice ridge. Renormalizing.') + aicen(1:ncat)=aicen(1:ncat)/sum(aicen(1:ncat)) + endif + ! ! output: snow/ice masses/thicknesses do k=1,nCat @@ -456,6 +468,7 @@ subroutine ice_ridging(IST, G, IG, mca_ice, mca_snow, mca_pond, TrReg, US, dt, r enddo + IST%part_size(i,j,0) = 1.0 - sum(IST%part_size(i,j,1:nCat)) endif