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_parameters.F90 b/config_src/external/Icepack_interfaces/icepack_parameters.F90 new file mode 100644 index 0000000000..637934d85a --- /dev/null +++ b/config_src/external/Icepack_interfaces/icepack_parameters.F90 @@ -0,0 +1,1628 @@ +!========================================================================= +! +! flags for the column package +! +! authors: Elizabeth C. Hunke, LANL + + module icepack_parameters + + use icepack_kinds + use icepack_warnings, only: icepack_warnings_aborted + + implicit none + private + + public :: icepack_init_parameters + public :: icepack_query_parameters + public :: icepack_write_parameters + public :: icepack_recompute_constants + + !----------------------------------------------------------------- + ! parameter constants + !----------------------------------------------------------------- + + integer (kind=int_kind), parameter, public :: & + nspint = 3 ! number of solar spectral intervals + + real (kind=dbl_kind), parameter, public :: & + c0 = 0.0_dbl_kind, & + c1 = 1.0_dbl_kind, & + c1p5 = 1.5_dbl_kind, & + c2 = 2.0_dbl_kind, & + c3 = 3.0_dbl_kind, & + c4 = 4.0_dbl_kind, & + c5 = 5.0_dbl_kind, & + c6 = 6.0_dbl_kind, & + c8 = 8.0_dbl_kind, & + c10 = 10.0_dbl_kind, & + c15 = 15.0_dbl_kind, & + c16 = 16.0_dbl_kind, & + c20 = 20.0_dbl_kind, & + c25 = 25.0_dbl_kind, & + c100 = 100.0_dbl_kind, & + c180 = 180.0_dbl_kind, & + c1000= 1000.0_dbl_kind, & + p001 = 0.001_dbl_kind, & + p01 = 0.01_dbl_kind, & + p1 = 0.1_dbl_kind, & + p2 = 0.2_dbl_kind, & + p4 = 0.4_dbl_kind, & + p5 = 0.5_dbl_kind, & + p6 = 0.6_dbl_kind, & + p05 = 0.05_dbl_kind, & + p15 = 0.15_dbl_kind, & + p25 = 0.25_dbl_kind, & + p75 = 0.75_dbl_kind, & + p333 = c1/c3, & + p666 = c2/c3, & + spval_const= -1.0e36_dbl_kind + + real (kind=dbl_kind), public :: & + secday = 86400.0_dbl_kind ,&! seconds in calendar day + puny = 1.0e-11_dbl_kind, & + bignum = 1.0e+30_dbl_kind, & + pi = 3.14159265358979323846_dbl_kind + + !----------------------------------------------------------------- + ! derived physical constants + ! Lfresh = Lsub-Lvap ,&! latent heat of melting of fresh ice (J/kg) + ! cprho = cp_ocn*rhow ,&! for ocean mixed layer (J kg / K m^3) + ! Cp = 0.5_dbl_kind*gravit*(rhow-rhoi)*rhoi/rhow ,&! proport const for PE + !----------------------------------------------------------------- + + real (kind=dbl_kind), public :: & + pih = spval_const ,&! 0.5 * pi + piq = spval_const ,&! 0.25 * pi + pi2 = spval_const ,&! 2 * pi + rad_to_deg = spval_const ,&! conversion factor, radians to degrees + Lfresh = spval_const ,&! latent heat of melting of fresh ice (J/kg) + cprho = spval_const ,&! for ocean mixed layer (J kg / K m^3) + Cp = spval_const ! proport const for PE + + !----------------------------------------------------------------- + ! Densities + !----------------------------------------------------------------- + + real (kind=dbl_kind), public :: & + rhos = 330.0_dbl_kind ,&! density of snow (kg/m^3) + rhoi = 917.0_dbl_kind ,&! density of ice (kg/m^3) + rhosi = 940.0_dbl_kind ,&! average sea ice density + ! Cox and Weeks, 1982: 919-974 kg/m^2 + rhow = 1026.0_dbl_kind ,&! density of seawater (kg/m^3) + rhofresh = 1000.0_dbl_kind ! density of fresh water (kg/m^3) + +!----------------------------------------------------------------------- +! Parameters for thermodynamics +!----------------------------------------------------------------------- + + real (kind=dbl_kind), public :: & + hfrazilmin = 0.05_dbl_kind ,&! min thickness of new frazil ice (m) + cp_ice = 2106._dbl_kind ,&! specific heat of fresh ice (J/kg/K) + cp_ocn = 4218._dbl_kind ,&! specific heat of ocn (J/kg/K) + ! freshwater value needed for enthalpy + depressT = 0.054_dbl_kind ,&! Tf:brine salinity ratio (C/ppt) + viscosity_dyn = 1.79e-3_dbl_kind, & ! dynamic viscosity of brine (kg/m/s) + Tocnfrz = -1.8_dbl_kind ,&! freezing temp of seawater (C), + ! used as Tsfcn for open water + Tffresh = 273.15_dbl_kind ,&! freezing temp of fresh ice (K) + Lsub = 2.835e6_dbl_kind ,&! latent heat, sublimation freshwater (J/kg) + Lvap = 2.501e6_dbl_kind ,&! latent heat, vaporization freshwater (J/kg) + Timelt = 0.0_dbl_kind ,&! melting temperature, ice top surface (C) + Tsmelt = 0.0_dbl_kind ,&! melting temperature, snow top surface (C) + ice_ref_salinity =4._dbl_kind,&! (ppt) + ! kice is not used for mushy thermo + kice = 2.03_dbl_kind ,&! thermal conductivity of fresh ice(W/m/deg) + ! kseaice is used only for zero-layer thermo + kseaice = 2.00_dbl_kind ,&! thermal conductivity of sea ice (W/m/deg) + ksno = 0.30_dbl_kind ,&! thermal conductivity of snow (W/m/deg) + hs_min = 1.e-4_dbl_kind ,&! min snow thickness for computing zTsn (m) + snowpatch = 0.02_dbl_kind ,&! parameter for fractional snow area (m) + saltmax = 3.2_dbl_kind ,&! max salinity at ice base for BL99 (ppt) + ! phi_init, dSin0_frazil are for mushy thermo + phi_init = 0.75_dbl_kind ,&! initial liquid fraction of frazil + min_salin = p1 ,&! threshold for brine pocket treatment + salt_loss = 0.4_dbl_kind ,&! fraction of salt retained in zsalinity + dSin0_frazil = c3 ,&! bulk salinity reduction of newly formed frazil + dts_b = 50._dbl_kind ,&! zsalinity timestep + ustar_min = 0.005_dbl_kind ,&! minimum friction velocity for ocean heat flux (m/s) + ! mushy thermo + a_rapid_mode = 0.5e-3_dbl_kind,&! channel radius for rapid drainage mode (m) + Rac_rapid_mode = 10.0_dbl_kind,&! critical Rayleigh number + aspect_rapid_mode = 1.0_dbl_kind,&! aspect ratio (larger is wider) + dSdt_slow_mode = -1.5e-7_dbl_kind,&! slow mode drainage strength (m s-1 K-1) + phi_c_slow_mode = 0.05_dbl_kind,&! critical liquid fraction porosity cutoff + phi_i_mushy = 0.85_dbl_kind ! liquid fraction of congelation ice + + integer (kind=int_kind), public :: & + ktherm = 1 ! type of thermodynamics + ! 0 = 0-layer approximation + ! 1 = Bitz and Lipscomb 1999 + ! 2 = mushy layer theory + + character (char_len), public :: & + conduct = 'bubbly', & ! 'MU71' or 'bubbly' + fbot_xfer_type = 'constant' ! transfer coefficient type for ice-ocean heat flux + + logical (kind=log_kind), public :: & + heat_capacity = .true. ,&! if true, ice has nonzero heat capacity + ! if false, use zero-layer thermodynamics + calc_Tsfc = .true. ,&! if true, calculate surface temperature + ! if false, Tsfc is computed elsewhere and + ! atmos-ice fluxes are provided to CICE + update_ocn_f = .false. ,&! include fresh water and salt fluxes for frazil + solve_zsal = .false. ,&! if true, update salinity profile from solve_S_dt + modal_aero = .false. ,&! if true, use modal aerosal optical properties + ! only for use with tr_aero or tr_zaero + conserv_check = .false. ! if true, do conservations checks and abort + + character(len=char_len), public :: & + tfrz_option = 'mushy' ! form of ocean freezing temperature + ! 'minus1p8' = -1.8 C + ! 'linear_salt' = -depressT * sss + ! 'mushy' conforms with ktherm=2 + +!----------------------------------------------------------------------- +! Parameters for radiation +!----------------------------------------------------------------------- + + real (kind=dbl_kind), public :: & + ! (Briegleb JGR 97 11475-11485 July 1992) + emissivity = 0.985_dbl_kind,&! emissivity of snow and ice + albocn = 0.06_dbl_kind ,&! ocean albedo + vonkar = 0.4_dbl_kind ,&! von Karman constant + stefan_boltzmann = 567.0e-10_dbl_kind,&! W/m^2/K^4 + ! (Ebert, Schramm and Curry JGR 100 15965-15975 Aug 1995) + kappav = 1.4_dbl_kind ,&! vis extnctn coef in ice, wvlngth<700nm (1/m) + hi_ssl = 0.050_dbl_kind,&! ice surface scattering layer thickness (m) + hs_ssl = 0.040_dbl_kind,&! snow surface scattering layer thickness (m) + ! baseline albedos for ccsm3 shortwave, set in namelist + albicev = 0.78_dbl_kind ,&! visible ice albedo for h > ahmax + albicei = 0.36_dbl_kind ,&! near-ir ice albedo for h > ahmax + albsnowv = 0.98_dbl_kind ,&! cold snow albedo, visible + albsnowi = 0.70_dbl_kind ,&! cold snow albedo, near IR + ahmax = 0.3_dbl_kind ,&! thickness above which ice albedo is constant (m) + ! dEdd tuning parameters, set in namelist + R_ice = c0 ,&! sea ice tuning parameter; +1 > 1sig increase in albedo + R_pnd = c0 ,&! ponded ice tuning parameter; +1 > 1sig increase in albedo + R_snw = c1p5 ,&! snow tuning parameter; +1 > ~.01 change in broadband albedo + dT_mlt = c1p5 ,&! change in temp for non-melt to melt snow grain + ! radius change (C) + rsnw_mlt = 1500._dbl_kind,&! maximum melting snow grain radius (10^-6 m) + kalg = 0.60_dbl_kind ! algae absorption coefficient for 0.5 m thick layer + ! 0.5 m path of 75 mg Chl a / m2 + ! weights for albedos + ! 4 Jan 2007 BPB Following are appropriate for complete cloud + ! in a summer polar atmosphere with 1.5m bare sea ice surface: + ! .636/.364 vis/nir with only 0.5% direct for each band. + real (kind=dbl_kind), public :: & ! currently used only + awtvdr = 0.00318_dbl_kind, &! visible, direct ! for history and + awtidr = 0.00182_dbl_kind, &! near IR, direct ! diagnostics + awtvdf = 0.63282_dbl_kind, &! visible, diffuse + awtidf = 0.36218_dbl_kind ! near IR, diffuse + + character (len=char_len), public :: & + shortwave = 'dEdd', & ! shortwave method, 'ccsm3' or 'dEdd' + albedo_type = 'ccsm3' ! albedo parameterization, 'ccsm3' or 'constant' + ! shortwave='dEdd' overrides this parameter + +!----------------------------------------------------------------------- +! Parameters for dynamics, including ridging and strength +!----------------------------------------------------------------------- + + integer (kind=int_kind), public :: & ! defined in namelist + kstrength = 1, & ! 0 for simple Hibler (1979) formulation + ! 1 for Rothrock (1975) pressure formulation + krdg_partic = 1, & ! 0 for Thorndike et al. (1975) formulation + ! 1 for exponential participation function + krdg_redist = 1 ! 0 for Hibler (1980) formulation + ! 1 for exponential redistribution function + + real (kind=dbl_kind), public :: & + Cf = 17._dbl_kind ,&! ratio of ridging work to PE change in ridging + Pstar = 2.75e4_dbl_kind ,&! constant in Hibler strength formula + ! (kstrength = 0) + Cstar = 20._dbl_kind ,&! constant in Hibler strength formula + ! (kstrength = 0) + dragio = 0.00536_dbl_kind ,&! ice-ocn drag coefficient + gravit = 9.80616_dbl_kind ,&! gravitational acceleration (m/s^2) + mu_rdg = 3.0_dbl_kind ! e-folding scale of ridged ice, krdg_partic=1 (m^0.5) + ! (krdg_redist = 1) + +!----------------------------------------------------------------------- +! Parameters for atmosphere +!----------------------------------------------------------------------- + + real (kind=dbl_kind), public :: & + cp_air = 1005.0_dbl_kind ,&! specific heat of air (J/kg/K) + cp_wv = 1.81e3_dbl_kind ,&! specific heat of water vapor (J/kg/K) + zvir = 0.606_dbl_kind ,&! rh2o/rair - 1.0 + zref = 10._dbl_kind ,&! reference height for stability (m) + iceruf = 0.0005_dbl_kind ,&! ice surface roughness (m) + qqqice = 11637800._dbl_kind ,&! for qsat over ice + TTTice = 5897.8_dbl_kind ,&! for qsat over ice + qqqocn = 627572.4_dbl_kind ,&! for qsat over ocn + TTTocn = 5107.4_dbl_kind ! for qsat over ocn + + character (len=char_len), public :: & + atmbndy = 'default' ! atmo boundary method, 'default' ('ccsm3') or 'constant' + + logical (kind=log_kind), public :: & + calc_strair = .true. , & ! if true, calculate wind stress + formdrag = .false. , & ! if true, calculate form drag + highfreq = .false. ! if true, calculate high frequency coupling + + integer (kind=int_kind), public :: & + natmiter = 5 ! number of iterations for atm boundary layer calcs + + ! Flux convergence tolerance + real (kind=dbl_kind), public :: atmiter_conv = c0 + +!----------------------------------------------------------------------- +! Parameters for the ice thickness distribution +!----------------------------------------------------------------------- + + integer (kind=int_kind), public :: & + kitd = 1 ,&! type of itd conversions + ! 0 = delta function + ! 1 = linear remap + kcatbound = 1 ! 0 = old category boundary formula + ! 1 = new formula giving round numbers + ! 2 = WMO standard + ! 3 = asymptotic formula + +!----------------------------------------------------------------------- +! Parameters for the floe size distribution +!----------------------------------------------------------------------- + + integer (kind=int_kind), public :: & + nfreq = 25 ! number of frequencies + + real (kind=dbl_kind), public :: & + floeshape = 0.66_dbl_kind ! constant from Steele (unitless) + + real (kind=dbl_kind), public :: & + floediam = 300.0_dbl_kind ! effective floe diameter for lateral melt (m) + + logical (kind=log_kind), public :: & + wave_spec = .false. ! if true, use wave forcing + + character (len=char_len), public :: & + wave_spec_type = 'constant' ! 'none', 'constant', or 'random' + +!----------------------------------------------------------------------- +! Parameters for melt ponds +!----------------------------------------------------------------------- + + real (kind=dbl_kind), public :: & + hs0 = 0.03_dbl_kind ! snow depth for transition to bare sea ice (m) + + ! level-ice ponds + character (len=char_len), public :: & + frzpnd = 'cesm' ! pond refreezing parameterization + + real (kind=dbl_kind), public :: & + dpscale = c1, & ! alter e-folding time scale for flushing + rfracmin = 0.15_dbl_kind, & ! minimum retained fraction of meltwater + rfracmax = 0.85_dbl_kind, & ! maximum retained fraction of meltwater + pndaspect = 0.8_dbl_kind, & ! ratio of pond depth to area fraction + hs1 = 0.03_dbl_kind ! snow depth for transition to bare pond ice (m) + + ! topo ponds + real (kind=dbl_kind), public :: & + hp1 = 0.01_dbl_kind ! critical pond lid thickness for topo ponds + +!----------------------------------------------------------------------- +! Parameters for biogeochemistry +!----------------------------------------------------------------------- + + character(char_len), public :: & + ! skl biology parameters + bgc_flux_type = 'Jin2006' ! type of ocean-ice poston velocity (or 'constant') + + logical (kind=log_kind), public :: & + z_tracers = .false., & ! if .true., bgc or aerosol tracers are vertically resolved + scale_bgc = .false., & ! if .true., initialize bgc tracers proportionally with salinity + solve_zbgc = .false., & ! if .true., solve vertical biochemistry portion of code + dEdd_algae = .false., & ! if .true., algal absorption of shortwave is computed in the + skl_bgc = .false. ! if true, solve skeletal biochemistry + + real (kind=dbl_kind), public :: & + phi_snow = p5 , & ! snow porosity + grid_o = c5 , & ! for bottom flux + initbio_frac = c1 , & ! fraction of ocean trcr concentration in bio trcrs + l_sk = 7.0_dbl_kind , & ! characteristic diffusive scale (m) + grid_oS = c5 , & ! for bottom flux + l_skS = 7.0_dbl_kind , & ! characteristic skeletal layer thickness (m) (zsalinity) + algal_vel = 1.11e-8_dbl_kind, & ! 0.5 cm/d(m/s) Lavoie 2005 1.5 cm/day + R_dFe2dust = 0.035_dbl_kind , & ! g/g (3.5% content) Tagliabue 2009 + dustFe_sol = 0.005_dbl_kind , & ! solubility fraction + frazil_scav = c1 , & ! fraction or multiple of bgc concentrated in frazil ice + sk_l = 0.03_dbl_kind , & ! skeletal layer thickness (m) + min_bgc = 0.01_dbl_kind , & ! fraction of ocean bgc concentration in surface melt + T_max = c0 , & ! maximum temperature (C) + fsal = c1 , & ! Salinity limitation (1) + op_dep_min = p1 , & ! light attenuates for optical depths exceeding min + fr_graze_s = p5 , & ! fraction of grazing spilled or slopped + fr_graze_e = p5 , & ! fraction of assimilation excreted + fr_mort2min = p5 , & ! fractionation of mortality to Am + fr_dFe = 0.3_dbl_kind , & ! fraction of remineralized nitrogen + ! (in units of algal iron) + k_nitrif = c0 , & ! nitrification rate (1/day) + t_iron_conv = 3065.0_dbl_kind , & ! desorption loss pFe to dFe (day) + max_loss = 0.9_dbl_kind , & ! restrict uptake to % of remaining value + max_dfe_doc1 = 0.2_dbl_kind , & ! max ratio of dFe to saccharides in the ice + ! (nM Fe/muM C) + fr_resp = 0.05_dbl_kind , & ! fraction of algal growth lost due to respiration + fr_resp_s = 0.75_dbl_kind , & ! DMSPd fraction of respiration loss as DMSPd + y_sk_DMS = p5 , & ! fraction conversion given high yield + t_sk_conv = 3.0_dbl_kind , & ! Stefels conversion time (d) + t_sk_ox = 10.0_dbl_kind ! DMS oxidation time (d) + + +!----------------------------------------------------------------------- +! Parameters for shortwave redistribution +!----------------------------------------------------------------------- + + logical (kind=log_kind), public :: & + sw_redist = .false. + + real (kind=dbl_kind), public :: & + sw_frac = 0.9_dbl_kind , & ! Fraction of internal shortwave moved to surface + sw_dtemp = 0.02_dbl_kind ! temperature difference from melting + +!======================================================================= + + contains + +!======================================================================= + +!autodocument_start icepack_init_parameters +! subroutine to set the column package internal parameters + + subroutine icepack_init_parameters( & + puny_in, bignum_in, pi_in, secday_in, & + rhos_in, rhoi_in, rhow_in, cp_air_in, emissivity_in, & + cp_ice_in, cp_ocn_in, hfrazilmin_in, floediam_in, & + depressT_in, dragio_in, albocn_in, gravit_in, viscosity_dyn_in, & + Tocnfrz_in, rhofresh_in, zvir_in, vonkar_in, cp_wv_in, & + stefan_boltzmann_in, ice_ref_salinity_in, & + Tffresh_in, Lsub_in, Lvap_in, Timelt_in, Tsmelt_in, & + iceruf_in, Cf_in, Pstar_in, Cstar_in, kappav_in, & + kice_in, kseaice_in, ksno_in, & + zref_in, hs_min_in, snowpatch_in, rhosi_in, sk_l_in, & + saltmax_in, phi_init_in, min_salin_in, salt_loss_in, & + min_bgc_in, dSin0_frazil_in, hi_ssl_in, hs_ssl_in, & + awtvdr_in, awtidr_in, awtvdf_in, awtidf_in, & + qqqice_in, TTTice_in, qqqocn_in, TTTocn_in, & + ktherm_in, conduct_in, fbot_xfer_type_in, calc_Tsfc_in, dts_b_in, & + update_ocn_f_in, ustar_min_in, a_rapid_mode_in, & + Rac_rapid_mode_in, aspect_rapid_mode_in, & + dSdt_slow_mode_in, phi_c_slow_mode_in, & + phi_i_mushy_in, shortwave_in, albedo_type_in, albsnowi_in, & + albicev_in, albicei_in, albsnowv_in, & + ahmax_in, R_ice_in, R_pnd_in, R_snw_in, dT_mlt_in, rsnw_mlt_in, & + kalg_in, kstrength_in, krdg_partic_in, krdg_redist_in, mu_rdg_in, & + atmbndy_in, calc_strair_in, formdrag_in, highfreq_in, natmiter_in, & + atmiter_conv_in, & + tfrz_option_in, kitd_in, kcatbound_in, hs0_in, frzpnd_in, & + floeshape_in, wave_spec_in, wave_spec_type_in, nfreq_in, & + dpscale_in, rfracmin_in, rfracmax_in, pndaspect_in, hs1_in, hp1_in, & + bgc_flux_type_in, z_tracers_in, scale_bgc_in, solve_zbgc_in, & + modal_aero_in, skl_bgc_in, solve_zsal_in, grid_o_in, l_sk_in, & + initbio_frac_in, grid_oS_in, l_skS_in, dEdd_algae_in, & + phi_snow_in, heat_capacity_in, T_max_in, fsal_in, & + fr_resp_in, algal_vel_in, R_dFe2dust_in, dustFe_sol_in, & + op_dep_min_in, fr_graze_s_in, fr_graze_e_in, fr_mort2min_in, & + fr_dFe_in, k_nitrif_in, t_iron_conv_in, max_loss_in, & + max_dfe_doc1_in, fr_resp_s_in, conserv_check_in, & + y_sk_DMS_in, t_sk_conv_in, t_sk_ox_in, frazil_scav_in, & + sw_redist_in, sw_frac_in, sw_dtemp_in) + + !----------------------------------------------------------------- + ! parameter constants + !----------------------------------------------------------------- + + real (kind=dbl_kind), intent(in), optional :: & + secday_in, & ! + puny_in, & ! + bignum_in, & ! + pi_in ! + + !----------------------------------------------------------------- + ! densities + !----------------------------------------------------------------- + + real (kind=dbl_kind), intent(in), optional :: & + rhos_in, & ! density of snow (kg/m^3) + rhoi_in, & ! density of ice (kg/m^3) + rhosi_in, & ! average sea ice density (kg/m2) + rhow_in, & ! density of seawater (kg/m^3) + rhofresh_in ! density of fresh water (kg/m^3) + +!----------------------------------------------------------------------- +! Parameters for thermodynamics +!----------------------------------------------------------------------- + + real (kind=dbl_kind), intent(in), optional :: & + floediam_in, & ! effective floe diameter for lateral melt (m) + hfrazilmin_in, & ! min thickness of new frazil ice (m) + cp_ice_in, & ! specific heat of fresh ice (J/kg/K) + cp_ocn_in, & ! specific heat of ocn (J/kg/K) + depressT_in, & ! Tf:brine salinity ratio (C/ppt) + viscosity_dyn_in, & ! dynamic viscosity of brine (kg/m/s) + Tocnfrz_in, & ! freezing temp of seawater (C) + Tffresh_in, & ! freezing temp of fresh ice (K) + Lsub_in, & ! latent heat, sublimation freshwater (J/kg) + Lvap_in, & ! latent heat, vaporization freshwater (J/kg) + Timelt_in, & ! melting temperature, ice top surface (C) + Tsmelt_in, & ! melting temperature, snow top surface (C) + ice_ref_salinity_in, & ! (ppt) + kice_in, & ! thermal conductivity of fresh ice(W/m/deg) + kseaice_in, & ! thermal conductivity of sea ice (W/m/deg) + ksno_in, & ! thermal conductivity of snow (W/m/deg) + hs_min_in, & ! min snow thickness for computing zTsn (m) + snowpatch_in, & ! parameter for fractional snow area (m) + saltmax_in, & ! max salinity at ice base for BL99 (ppt) + phi_init_in, & ! initial liquid fraction of frazil + min_salin_in, & ! threshold for brine pocket treatment + salt_loss_in, & ! fraction of salt retained in zsalinity + dSin0_frazil_in ! bulk salinity reduction of newly formed frazil + + integer (kind=int_kind), intent(in), optional :: & + ktherm_in ! type of thermodynamics + ! 0 = 0-layer approximation + ! 1 = Bitz and Lipscomb 1999 + ! 2 = mushy layer theory + + character (char_len), intent(in), optional :: & + conduct_in, & ! 'MU71' or 'bubbly' + fbot_xfer_type_in ! transfer coefficient type for ice-ocean heat flux + + logical (kind=log_kind), intent(in), optional :: & + heat_capacity_in, &! if true, ice has nonzero heat capacity + ! if false, use zero-layer thermodynamics + calc_Tsfc_in , &! if true, calculate surface temperature + ! if false, Tsfc is computed elsewhere and + ! atmos-ice fluxes are provided to CICE + update_ocn_f_in ! include fresh water and salt fluxes for frazil + + real (kind=dbl_kind), intent(in), optional :: & + dts_b_in, & ! zsalinity timestep + ustar_min_in ! minimum friction velocity for ice-ocean heat flux + + ! mushy thermo + real(kind=dbl_kind), intent(in), optional :: & + a_rapid_mode_in , & ! channel radius for rapid drainage mode (m) + Rac_rapid_mode_in , & ! critical Rayleigh number for rapid drainage mode + aspect_rapid_mode_in , & ! aspect ratio for rapid drainage mode (larger=wider) + dSdt_slow_mode_in , & ! slow mode drainage strength (m s-1 K-1) + phi_c_slow_mode_in , & ! liquid fraction porosity cutoff for slow mode + phi_i_mushy_in ! liquid fraction of congelation ice + + character(len=char_len), intent(in), optional :: & + tfrz_option_in ! form of ocean freezing temperature + ! 'minus1p8' = -1.8 C + ! 'linear_salt' = -depressT * sss + ! 'mushy' conforms with ktherm=2 + +!----------------------------------------------------------------------- +! Parameters for radiation +!----------------------------------------------------------------------- + + real(kind=dbl_kind), intent(in), optional :: & + emissivity_in, & ! emissivity of snow and ice + albocn_in, & ! ocean albedo + vonkar_in, & ! von Karman constant + stefan_boltzmann_in, & ! W/m^2/K^4 + kappav_in, & ! vis extnctn coef in ice, wvlngth<700nm (1/m) + hi_ssl_in, & ! ice surface scattering layer thickness (m) + hs_ssl_in, & ! visible, direct + awtvdr_in, & ! visible, direct ! for history and + awtidr_in, & ! near IR, direct ! diagnostics + awtvdf_in, & ! visible, diffuse + awtidf_in ! near IR, diffuse + + character (len=char_len), intent(in), optional :: & + shortwave_in, & ! shortwave method, 'ccsm3' or 'dEdd' + albedo_type_in ! albedo parameterization, 'ccsm3' or 'constant' + ! shortwave='dEdd' overrides this parameter + + ! baseline albedos for ccsm3 shortwave, set in namelist + real (kind=dbl_kind), intent(in), optional :: & + albicev_in , & ! visible ice albedo for h > ahmax + albicei_in , & ! near-ir ice albedo for h > ahmax + albsnowv_in , & ! cold snow albedo, visible + albsnowi_in , & ! cold snow albedo, near IR + ahmax_in ! thickness above which ice albedo is constant (m) + + ! dEdd tuning parameters, set in namelist + real (kind=dbl_kind), intent(in), optional :: & + R_ice_in , & ! sea ice tuning parameter; +1 > 1sig increase in albedo + R_pnd_in , & ! ponded ice tuning parameter; +1 > 1sig increase in albedo + R_snw_in , & ! snow tuning parameter; +1 > ~.01 change in broadband albedo + dT_mlt_in , & ! change in temp for non-melt to melt snow grain + ! radius change (C) + rsnw_mlt_in , & ! maximum melting snow grain radius (10^-6 m) + kalg_in ! algae absorption coefficient for 0.5 m thick layer + + logical (kind=log_kind), intent(in), optional :: & + sw_redist_in ! redistribute shortwave + + real (kind=dbl_kind), intent(in), optional :: & + sw_frac_in , & ! Fraction of internal shortwave moved to surface + sw_dtemp_in ! temperature difference from melting + +!----------------------------------------------------------------------- +! Parameters for dynamics +!----------------------------------------------------------------------- + + real(kind=dbl_kind), intent(in), optional :: & + Cf_in, & ! ratio of ridging work to PE change in ridging + Pstar_in, & ! constant in Hibler strength formula + Cstar_in, & ! constant in Hibler strength formula + dragio_in, & ! ice-ocn drag coefficient + gravit_in, & ! gravitational acceleration (m/s^2) + iceruf_in ! ice surface roughness (m) + + integer (kind=int_kind), intent(in), optional :: & ! defined in namelist + kstrength_in , & ! 0 for simple Hibler (1979) formulation + ! 1 for Rothrock (1975) pressure formulation + krdg_partic_in, & ! 0 for Thorndike et al. (1975) formulation + ! 1 for exponential participation function + krdg_redist_in ! 0 for Hibler (1980) formulation + ! 1 for exponential redistribution function + + real (kind=dbl_kind), intent(in), optional :: & + mu_rdg_in ! gives e-folding scale of ridged ice (m^.5) + ! (krdg_redist = 1) + +!----------------------------------------------------------------------- +! Parameters for atmosphere +!----------------------------------------------------------------------- + + real (kind=dbl_kind), intent(in), optional :: & + cp_air_in, & ! specific heat of air (J/kg/K) + cp_wv_in, & ! specific heat of water vapor (J/kg/K) + zvir_in, & ! rh2o/rair - 1.0 + zref_in, & ! reference height for stability (m) + qqqice_in, & ! for qsat over ice + TTTice_in, & ! for qsat over ice + qqqocn_in, & ! for qsat over ocn + TTTocn_in ! for qsat over ocn + + character (len=char_len), intent(in), optional :: & + atmbndy_in ! atmo boundary method, 'default' ('ccsm3') or 'constant' + + logical (kind=log_kind), intent(in), optional :: & + calc_strair_in, & ! if true, calculate wind stress components + formdrag_in, & ! if true, calculate form drag + highfreq_in ! if true, use high frequency coupling + + integer (kind=int_kind), intent(in), optional :: & + natmiter_in ! number of iterations for boundary layer calculations + + ! Flux convergence tolerance + real (kind=dbl_kind), intent(in), optional :: atmiter_conv_in + +!----------------------------------------------------------------------- +! Parameters for the ice thickness distribution +!----------------------------------------------------------------------- + + integer (kind=int_kind), intent(in), optional :: & + kitd_in , & ! type of itd conversions + ! 0 = delta function + ! 1 = linear remap + kcatbound_in ! 0 = old category boundary formula + ! 1 = new formula giving round numbers + ! 2 = WMO standard + ! 3 = asymptotic formula + +!----------------------------------------------------------------------- +! Parameters for the floe size distribution +!----------------------------------------------------------------------- + + integer (kind=int_kind), intent(in), optional :: & + nfreq_in ! number of frequencies + + real (kind=dbl_kind), intent(in), optional :: & + floeshape_in ! constant from Steele (unitless) + + logical (kind=log_kind), intent(in), optional :: & + wave_spec_in ! if true, use wave forcing + + character (len=char_len), intent(in), optional :: & + wave_spec_type_in ! type of wave spectrum forcing + +!----------------------------------------------------------------------- +! Parameters for biogeochemistry +!----------------------------------------------------------------------- + + character(char_len), intent(in), optional :: & + bgc_flux_type_in ! type of ocean-ice piston velocity + ! 'constant', 'Jin2006' + + logical (kind=log_kind), intent(in), optional :: & + z_tracers_in, & ! if .true., bgc or aerosol tracers are vertically resolved + scale_bgc_in, & ! if .true., initialize bgc tracers proportionally with salinity + solve_zbgc_in, & ! if .true., solve vertical biochemistry portion of code + dEdd_algae_in, & ! if .true., algal absorptionof Shortwave is computed in the + modal_aero_in, & ! if .true., use modal aerosol formulation in shortwave + conserv_check_in ! if .true., run conservation checks and abort if checks fail + + logical (kind=log_kind), intent(in), optional :: & + skl_bgc_in, & ! if true, solve skeletal biochemistry + solve_zsal_in ! if true, update salinity profile from solve_S_dt + + real (kind=dbl_kind), intent(in), optional :: & + grid_o_in , & ! for bottom flux + l_sk_in , & ! characteristic diffusive scale (zsalinity) (m) + initbio_frac_in, & ! fraction of ocean tracer concentration used to initialize tracer + phi_snow_in ! snow porosity at the ice/snow interface + + real (kind=dbl_kind), intent(in), optional :: & + grid_oS_in , & ! for bottom flux (zsalinity) + l_skS_in ! 0.02 characteristic skeletal layer thickness (m) (zsalinity) + real (kind=dbl_kind), intent(in), optional :: & + fr_resp_in , & ! fraction of algal growth lost due to respiration + algal_vel_in , & ! 0.5 cm/d(m/s) Lavoie 2005 1.5 cm/day + R_dFe2dust_in , & ! g/g (3.5% content) Tagliabue 2009 + dustFe_sol_in , & ! solubility fraction + T_max_in , & ! maximum temperature (C) + fsal_in , & ! Salinity limitation (ppt) + op_dep_min_in , & ! Light attenuates for optical depths exceeding min + fr_graze_s_in , & ! fraction of grazing spilled or slopped + fr_graze_e_in , & ! fraction of assimilation excreted + fr_mort2min_in , & ! fractionation of mortality to Am + fr_dFe_in , & ! fraction of remineralized nitrogen + ! (in units of algal iron) + k_nitrif_in , & ! nitrification rate (1/day) + t_iron_conv_in , & ! desorption loss pFe to dFe (day) + max_loss_in , & ! restrict uptake to % of remaining value + max_dfe_doc1_in , & ! max ratio of dFe to saccharides in the ice + ! (nM Fe/muM C) + fr_resp_s_in , & ! DMSPd fraction of respiration loss as DMSPd + y_sk_DMS_in , & ! fraction conversion given high yield + t_sk_conv_in , & ! Stefels conversion time (d) + t_sk_ox_in , & ! DMS oxidation time (d) + frazil_scav_in ! scavenging fraction or multiple in frazil ice + + real (kind=dbl_kind), intent(in), optional :: & + sk_l_in, & ! skeletal layer thickness (m) + min_bgc_in ! fraction of ocean bgc concentration in surface melt + +!----------------------------------------------------------------------- +! Parameters for melt ponds +!----------------------------------------------------------------------- + + real (kind=dbl_kind), intent(in), optional :: & + hs0_in ! snow depth for transition to bare sea ice (m) + + ! level-ice ponds + character (len=char_len), intent(in), optional :: & + frzpnd_in ! pond refreezing parameterization + + real (kind=dbl_kind), intent(in), optional :: & + dpscale_in, & ! alter e-folding time scale for flushing + rfracmin_in, & ! minimum retained fraction of meltwater + rfracmax_in, & ! maximum retained fraction of meltwater + pndaspect_in, & ! ratio of pond depth to pond fraction + hs1_in ! tapering parameter for snow on pond ice + + ! topo ponds + real (kind=dbl_kind), intent(in), optional :: & + hp1_in ! critical parameter for pond ice thickness + +!autodocument_end + + character(len=*),parameter :: subname='(icepack_init_parameters)' + + if (present(rhos_in) ) rhos = rhos_in + if (present(rhoi_in) ) rhoi = rhoi_in + if (present(rhow_in) ) rhow = rhow_in + if (present(cp_air_in) ) cp_air = cp_air_in + if (present(emissivity_in) ) emissivity = emissivity_in + if (present(floediam_in) ) floediam = floediam_in + if (present(hfrazilmin_in) ) hfrazilmin = hfrazilmin_in + if (present(cp_ice_in) ) cp_ice = cp_ice_in + if (present(cp_ocn_in) ) cp_ocn = cp_ocn_in + if (present(depressT_in) ) depressT = depressT_in + if (present(dragio_in) ) dragio = dragio_in + if (present(albocn_in) ) albocn = albocn_in + if (present(gravit_in) ) gravit = gravit_in + if (present(viscosity_dyn_in) ) viscosity_dyn = viscosity_dyn_in + if (present(Tocnfrz_in) ) Tocnfrz = Tocnfrz_in + if (present(rhofresh_in) ) rhofresh = rhofresh_in + if (present(zvir_in) ) zvir = zvir_in + if (present(vonkar_in) ) vonkar = vonkar_in + if (present(cp_wv_in) ) cp_wv = cp_wv_in + if (present(stefan_boltzmann_in) ) stefan_boltzmann = stefan_boltzmann_in + if (present(Tffresh_in) ) Tffresh = Tffresh_in + if (present(Lsub_in) ) Lsub = Lsub_in + if (present(Lvap_in) ) Lvap = Lvap_in + if (present(Timelt_in) ) Timelt = Timelt_in + if (present(Tsmelt_in) ) Tsmelt = Tsmelt_in + if (present(ice_ref_salinity_in) ) ice_ref_salinity = ice_ref_salinity_in + if (present(iceruf_in) ) iceruf = iceruf_in + if (present(Cf_in) ) Cf = Cf_in + if (present(Pstar_in) ) Pstar = Pstar_in + if (present(Cstar_in) ) Cstar = Cstar_in + if (present(kappav_in) ) kappav = kappav_in + if (present(kice_in) ) kice = kice_in + if (present(kseaice_in) ) kseaice = kseaice_in + if (present(ksno_in) ) ksno = ksno_in + if (present(zref_in) ) zref = zref_in + if (present(hs_min_in) ) hs_min = hs_min_in + if (present(snowpatch_in) ) snowpatch = snowpatch_in + if (present(rhosi_in) ) rhosi = rhosi_in + if (present(sk_l_in) ) sk_l = sk_l_in + if (present(saltmax_in) ) saltmax = saltmax_in + if (present(phi_init_in) ) phi_init = phi_init_in + if (present(min_salin_in) ) min_salin = min_salin_in + if (present(salt_loss_in) ) salt_loss = salt_loss_in + if (present(min_bgc_in) ) min_bgc = min_bgc_in + if (present(dSin0_frazil_in) ) dSin0_frazil = dSin0_frazil_in + if (present(hi_ssl_in) ) hi_ssl = hi_ssl_in + if (present(hs_ssl_in) ) hs_ssl = hs_ssl_in + if (present(awtvdr_in) ) awtvdr = awtvdr_in + if (present(awtidr_in) ) awtidr = awtidr_in + if (present(awtvdf_in) ) awtvdf = awtvdf_in + if (present(awtidf_in) ) awtidf = awtidf_in + if (present(qqqice_in) ) qqqice = qqqice_in + if (present(TTTice_in) ) TTTice = TTTice_in + if (present(qqqocn_in) ) qqqocn = qqqocn_in + if (present(TTTocn_in) ) TTTocn = TTTocn_in + if (present(puny_in) ) puny = puny_in + if (present(bignum_in) ) bignum = bignum_in + if (present(pi_in) ) pi = pi_in + if (present(secday_in) ) secday = secday_in + if (present(ktherm_in) ) ktherm = ktherm_in + if (present(conduct_in) ) conduct = conduct_in + if (present(fbot_xfer_type_in) ) fbot_xfer_type = fbot_xfer_type_in + if (present(heat_capacity_in) ) heat_capacity = heat_capacity_in + if (present(calc_Tsfc_in) ) calc_Tsfc = calc_Tsfc_in + if (present(update_ocn_f_in) ) update_ocn_f = update_ocn_f_in + if (present(dts_b_in) ) dts_b = dts_b_in + if (present(ustar_min_in) ) ustar_min = ustar_min_in + if (present(a_rapid_mode_in) ) a_rapid_mode = a_rapid_mode_in + if (present(Rac_rapid_mode_in) ) Rac_rapid_mode = Rac_rapid_mode_in + if (present(aspect_rapid_mode_in) ) aspect_rapid_mode= aspect_rapid_mode_in + if (present(dSdt_slow_mode_in) ) dSdt_slow_mode = dSdt_slow_mode_in + if (present(phi_c_slow_mode_in) ) phi_c_slow_mode = phi_c_slow_mode_in + if (present(phi_i_mushy_in) ) phi_i_mushy = phi_i_mushy_in + if (present(shortwave_in) ) shortwave = shortwave_in + if (present(albedo_type_in) ) albedo_type = albedo_type_in + if (present(albicev_in) ) albicev = albicev_in + if (present(albicei_in) ) albicei = albicei_in + if (present(albsnowv_in) ) albsnowv = albsnowv_in + if (present(albsnowi_in) ) albsnowi = albsnowi_in + if (present(ahmax_in) ) ahmax = ahmax_in + if (present(R_ice_in) ) R_ice = R_ice_in + if (present(R_pnd_in) ) R_pnd = R_pnd_in + if (present(R_snw_in) ) R_snw = R_snw_in + if (present(dT_mlt_in) ) dT_mlt = dT_mlt_in + if (present(rsnw_mlt_in) ) rsnw_mlt = rsnw_mlt_in + if (present(kalg_in) ) kalg = kalg_in + if (present(kstrength_in) ) kstrength = kstrength_in + if (present(krdg_partic_in) ) krdg_partic = krdg_partic_in + if (present(krdg_redist_in) ) krdg_redist = krdg_redist_in + if (present(mu_rdg_in) ) mu_rdg = mu_rdg_in + if (present(atmbndy_in) ) atmbndy = atmbndy_in + if (present(calc_strair_in) ) calc_strair = calc_strair_in + if (present(formdrag_in) ) formdrag = formdrag_in + if (present(highfreq_in) ) highfreq = highfreq_in + if (present(natmiter_in) ) natmiter = natmiter_in + if (present(atmiter_conv_in) ) atmiter_conv = atmiter_conv_in + if (present(tfrz_option_in) ) tfrz_option = tfrz_option_in + if (present(kitd_in) ) kitd = kitd_in + if (present(kcatbound_in) ) kcatbound = kcatbound_in + if (present(floeshape_in) ) floeshape = floeshape_in + if (present(wave_spec_in) ) wave_spec = wave_spec_in + if (present(wave_spec_type_in) ) wave_spec_type = wave_spec_type_in + if (present(nfreq_in) ) nfreq = nfreq_in + if (present(hs0_in) ) hs0 = hs0_in + if (present(frzpnd_in) ) frzpnd = frzpnd_in + if (present(dpscale_in) ) dpscale = dpscale_in + if (present(rfracmin_in) ) rfracmin = rfracmin_in + if (present(rfracmax_in) ) rfracmax = rfracmax_in + if (present(pndaspect_in) ) pndaspect = pndaspect_in + if (present(hs1_in) ) hs1 = hs1_in + if (present(hp1_in) ) hp1 = hp1_in + if (present(bgc_flux_type_in) ) bgc_flux_type = bgc_flux_type_in + if (present(z_tracers_in) ) z_tracers = z_tracers_in + if (present(scale_bgc_in) ) scale_bgc = scale_bgc_in + if (present(solve_zbgc_in) ) solve_zbgc = solve_zbgc_in + if (present(dEdd_algae_in) ) dEdd_algae = dEdd_algae_in + if (present(modal_aero_in) ) modal_aero = modal_aero_in + if (present(conserv_check_in) ) conserv_check = conserv_check_in + if (present(skl_bgc_in) ) skl_bgc = skl_bgc_in + if (present(solve_zsal_in) ) solve_zsal = solve_zsal_in + if (present(grid_o_in) ) grid_o = grid_o_in + if (present(l_sk_in) ) l_sk = l_sk_in + if (present(initbio_frac_in) ) initbio_frac = initbio_frac_in + if (present(grid_oS_in) ) grid_oS = grid_oS_in + if (present(l_skS_in) ) l_skS = l_skS_in + if (present(phi_snow_in) ) phi_snow = phi_snow_in + if (present(fr_resp_in) ) fr_resp = fr_resp_in + if (present(algal_vel_in) ) algal_vel = algal_vel_in + if (present(R_dFe2dust_in) ) R_dFe2dust = R_dFe2dust_in + if (present(dustFe_sol_in) ) dustFe_sol = dustFe_sol_in + if (present(T_max_in) ) T_max = T_max_in + if (present(fsal_in) ) fsal = fsal_in + if (present(op_dep_min_in) ) op_dep_min = op_dep_min_in + if (present(fr_graze_s_in) ) fr_graze_s = fr_graze_s_in + if (present(fr_graze_e_in) ) fr_graze_e = fr_graze_e_in + if (present(fr_mort2min_in) ) fr_mort2min = fr_mort2min_in + if (present(fr_dFe_in) ) fr_dFe = fr_dFe_in + if (present(k_nitrif_in) ) k_nitrif = k_nitrif_in + if (present(t_iron_conv_in) ) t_iron_conv = t_iron_conv_in + if (present(max_loss_in) ) max_loss = max_loss_in + if (present(max_dfe_doc1_in) ) max_dfe_doc1 = max_dfe_doc1_in + if (present(fr_resp_s_in) ) fr_resp_s = fr_resp_s_in + if (present(y_sk_DMS_in) ) y_sk_DMS = y_sk_DMS_in + if (present(t_sk_conv_in) ) t_sk_conv = t_sk_conv_in + if (present(t_sk_ox_in) ) t_sk_ox = t_sk_ox_in + if (present(frazil_scav_in) ) frazil_scav = frazil_scav_in + if (present(sw_redist_in) ) sw_redist = sw_redist_in + if (present(sw_frac_in) ) sw_frac = sw_frac_in + if (present(sw_dtemp_in) ) sw_dtemp = sw_dtemp_in + + call icepack_recompute_constants() + if (icepack_warnings_aborted(subname)) return + + end subroutine icepack_init_parameters + +!======================================================================= + +!autodocument_start icepack_query_parameters +! subroutine to query the column package internal parameters + + subroutine icepack_query_parameters( & + puny_out, bignum_out, pi_out, rad_to_deg_out,& + secday_out, c0_out, c1_out, c1p5_out, c2_out, c3_out, c4_out, & + c5_out, c6_out, c8_out, c10_out, c15_out, c16_out, c20_out, & + c25_out, c100_out, c180_out, c1000_out, p001_out, p01_out, p1_out, & + p2_out, p4_out, p5_out, p6_out, p05_out, p15_out, p25_out, p75_out, & + p333_out, p666_out, spval_const_out, pih_out, piq_out, pi2_out, & + rhos_out, rhoi_out, rhow_out, cp_air_out, emissivity_out, & + cp_ice_out, cp_ocn_out, hfrazilmin_out, floediam_out, & + depressT_out, dragio_out, albocn_out, gravit_out, viscosity_dyn_out, & + Tocnfrz_out, rhofresh_out, zvir_out, vonkar_out, cp_wv_out, & + stefan_boltzmann_out, ice_ref_salinity_out, & + Tffresh_out, Lsub_out, Lvap_out, Timelt_out, Tsmelt_out, & + iceruf_out, Cf_out, Pstar_out, Cstar_out, kappav_out, & + kice_out, kseaice_out, ksno_out, & + zref_out, hs_min_out, snowpatch_out, rhosi_out, sk_l_out, & + saltmax_out, phi_init_out, min_salin_out, salt_loss_out, & + min_bgc_out, dSin0_frazil_out, hi_ssl_out, hs_ssl_out, & + awtvdr_out, awtidr_out, awtvdf_out, awtidf_out, & + qqqice_out, TTTice_out, qqqocn_out, TTTocn_out, update_ocn_f_out, & + Lfresh_out, cprho_out, Cp_out, ustar_min_out, a_rapid_mode_out, & + ktherm_out, conduct_out, fbot_xfer_type_out, calc_Tsfc_out, dts_b_out, & + Rac_rapid_mode_out, aspect_rapid_mode_out, dSdt_slow_mode_out, & + phi_c_slow_mode_out, phi_i_mushy_out, shortwave_out, & + albedo_type_out, albicev_out, albicei_out, albsnowv_out, & + albsnowi_out, ahmax_out, R_ice_out, R_pnd_out, R_snw_out, dT_mlt_out, & + rsnw_mlt_out, dEdd_algae_out, & + kalg_out, kstrength_out, krdg_partic_out, krdg_redist_out, mu_rdg_out, & + atmbndy_out, calc_strair_out, formdrag_out, highfreq_out, natmiter_out, & + atmiter_conv_out, & + tfrz_option_out, kitd_out, kcatbound_out, hs0_out, frzpnd_out, & + floeshape_out, wave_spec_out, wave_spec_type_out, nfreq_out, & + dpscale_out, rfracmin_out, rfracmax_out, pndaspect_out, hs1_out, hp1_out, & + bgc_flux_type_out, z_tracers_out, scale_bgc_out, solve_zbgc_out, & + modal_aero_out, skl_bgc_out, solve_zsal_out, grid_o_out, l_sk_out, & + initbio_frac_out, grid_oS_out, l_skS_out, & + phi_snow_out, heat_capacity_out, conserv_check_out, & + fr_resp_out, algal_vel_out, R_dFe2dust_out, dustFe_sol_out, & + T_max_out, fsal_out, op_dep_min_out, fr_graze_s_out, fr_graze_e_out, & + fr_mort2min_out, fr_resp_s_out, fr_dFe_out, & + k_nitrif_out, t_iron_conv_out, max_loss_out, max_dfe_doc1_out, & + y_sk_DMS_out, t_sk_conv_out, t_sk_ox_out, frazil_scav_out, & + sw_redist_out, sw_frac_out, sw_dtemp_out) + + !----------------------------------------------------------------- + ! parameter constants + !----------------------------------------------------------------- + + real (kind=dbl_kind), intent(out), optional :: & + c0_out, c1_out, c1p5_out, c2_out, c3_out, c4_out, & + c5_out, c6_out, c8_out, c10_out, c15_out, c16_out, c20_out, & + c25_out, c180_out, c100_out, c1000_out, p001_out, p01_out, p1_out, & + p2_out, p4_out, p5_out, p6_out, p05_out, p15_out, p25_out, p75_out, & + p333_out, p666_out, spval_const_out, pih_out, piq_out, pi2_out, & + secday_out, & ! number of seconds per day + puny_out, & ! a small number + bignum_out, & ! a big number + pi_out, & ! pi + rad_to_deg_out, & ! conversion factor from radians to degrees + Lfresh_out, & ! latent heat of melting of fresh ice (J/kg) + cprho_out, & ! for ocean mixed layer (J kg / K m^3) + Cp_out ! proport const for PE + + !----------------------------------------------------------------- + ! densities + !----------------------------------------------------------------- + + real (kind=dbl_kind), intent(out), optional :: & + rhos_out, & ! density of snow (kg/m^3) + rhoi_out, & ! density of ice (kg/m^3) + rhosi_out, & ! average sea ice density (kg/m2) + rhow_out, & ! density of seawater (kg/m^3) + rhofresh_out ! density of fresh water (kg/m^3) + +!----------------------------------------------------------------------- +! Parameters for thermodynamics +!----------------------------------------------------------------------- + + real (kind=dbl_kind), intent(out), optional :: & + floediam_out, & ! effective floe diameter for lateral melt (m) + hfrazilmin_out, & ! min thickness of new frazil ice (m) + cp_ice_out, & ! specific heat of fresh ice (J/kg/K) + cp_ocn_out, & ! specific heat of ocn (J/kg/K) + depressT_out, & ! Tf:brine salinity ratio (C/ppt) + viscosity_dyn_out, & ! dynamic viscosity of brine (kg/m/s) + Tocnfrz_out, & ! freezing temp of seawater (C) + Tffresh_out, & ! freezing temp of fresh ice (K) + Lsub_out, & ! latent heat, sublimation freshwater (J/kg) + Lvap_out, & ! latent heat, vaporization freshwater (J/kg) + Timelt_out, & ! melting temperature, ice top surface (C) + Tsmelt_out, & ! melting temperature, snow top surface (C) + ice_ref_salinity_out, & ! (ppt) + kice_out, & ! thermal conductivity of fresh ice(W/m/deg) + kseaice_out, & ! thermal conductivity of sea ice (W/m/deg) + ksno_out, & ! thermal conductivity of snow (W/m/deg) + hs_min_out, & ! min snow thickness for computing zTsn (m) + snowpatch_out, & ! parameter for fractional snow area (m) + saltmax_out, & ! max salinity at ice base for BL99 (ppt) + phi_init_out, & ! initial liquid fraction of frazil + min_salin_out, & ! threshold for brine pocket treatment + salt_loss_out, & ! fraction of salt retained in zsalinity + dSin0_frazil_out ! bulk salinity reduction of newly formed frazil + + integer (kind=int_kind), intent(out), optional :: & + ktherm_out ! type of thermodynamics + ! 0 = 0-layer approximation + ! 1 = Bitz and Lipscomb 1999 + ! 2 = mushy layer theory + + character (char_len), intent(out), optional :: & + conduct_out, & ! 'MU71' or 'bubbly' + fbot_xfer_type_out ! transfer coefficient type for ice-ocean heat flux + + logical (kind=log_kind), intent(out), optional :: & + heat_capacity_out,&! if true, ice has nonzero heat capacity + ! if false, use zero-layer thermodynamics + calc_Tsfc_out ,&! if true, calculate surface temperature + ! if false, Tsfc is computed elsewhere and + ! atmos-ice fluxes are provided to CICE + update_ocn_f_out ! include fresh water and salt fluxes for frazil + + real (kind=dbl_kind), intent(out), optional :: & + dts_b_out, & ! zsalinity timestep + ustar_min_out ! minimum friction velocity for ice-ocean heat flux + + ! mushy thermo + real(kind=dbl_kind), intent(out), optional :: & + a_rapid_mode_out , & ! channel radius for rapid drainage mode (m) + Rac_rapid_mode_out , & ! critical Rayleigh number for rapid drainage mode + aspect_rapid_mode_out , & ! aspect ratio for rapid drainage mode (larger=wider) + dSdt_slow_mode_out , & ! slow mode drainage strength (m s-1 K-1) + phi_c_slow_mode_out , & ! liquid fraction porosity cutoff for slow mode + phi_i_mushy_out ! liquid fraction of congelation ice + + character(len=char_len), intent(out), optional :: & + tfrz_option_out ! form of ocean freezing temperature + ! 'minus1p8' = -1.8 C + ! 'linear_salt' = -depressT * sss + ! 'mushy' conforms with ktherm=2 + +!----------------------------------------------------------------------- +! Parameters for radiation +!----------------------------------------------------------------------- + + real(kind=dbl_kind), intent(out), optional :: & + emissivity_out, & ! emissivity of snow and ice + albocn_out, & ! ocean albedo + vonkar_out, & ! von Karman constant + stefan_boltzmann_out, & ! W/m^2/K^4 + kappav_out, & ! vis extnctn coef in ice, wvlngth<700nm (1/m) + hi_ssl_out, & ! ice surface scattering layer thickness (m) + hs_ssl_out, & ! visible, direct + awtvdr_out, & ! visible, direct ! for history and + awtidr_out, & ! near IR, direct ! diagnostics + awtvdf_out, & ! visible, diffuse + awtidf_out ! near IR, diffuse + + character (len=char_len), intent(out), optional :: & + shortwave_out, & ! shortwave method, 'ccsm3' or 'dEdd' + albedo_type_out ! albedo parameterization, 'ccsm3' or 'constant' + ! shortwave='dEdd' overrides this parameter + + ! baseline albedos for ccsm3 shortwave, set in namelist + real (kind=dbl_kind), intent(out), optional :: & + albicev_out , & ! visible ice albedo for h > ahmax + albicei_out , & ! near-ir ice albedo for h > ahmax + albsnowv_out , & ! cold snow albedo, visible + albsnowi_out , & ! cold snow albedo, near IR + ahmax_out ! thickness above which ice albedo is constant (m) + + ! dEdd tuning parameters, set in namelist + real (kind=dbl_kind), intent(out), optional :: & + R_ice_out , & ! sea ice tuning parameter; +1 > 1sig increase in albedo + R_pnd_out , & ! ponded ice tuning parameter; +1 > 1sig increase in albedo + R_snw_out , & ! snow tuning parameter; +1 > ~.01 change in broadband albedo + dT_mlt_out , & ! change in temp for non-melt to melt snow grain + ! radius change (C) + rsnw_mlt_out , & ! maximum melting snow grain radius (10^-6 m) + kalg_out ! algae absorption coefficient for 0.5 m thick layer + + logical (kind=log_kind), intent(out), optional :: & + sw_redist_out ! redistribute shortwave + + real (kind=dbl_kind), intent(out), optional :: & + sw_frac_out , & ! Fraction of internal shortwave moved to surface + sw_dtemp_out ! temperature difference from melting + +!----------------------------------------------------------------------- +! Parameters for dynamics +!----------------------------------------------------------------------- + + real(kind=dbl_kind), intent(out), optional :: & + Cf_out, & ! ratio of ridging work to PE change in ridging + Pstar_out, & ! constant in Hibler strength formula + Cstar_out, & ! constant in Hibler strength formula + dragio_out, & ! ice-ocn drag coefficient + gravit_out, & ! gravitational acceleration (m/s^2) + iceruf_out ! ice surface roughness (m) + + integer (kind=int_kind), intent(out), optional :: & ! defined in namelist + kstrength_out , & ! 0 for simple Hibler (1979) formulation + ! 1 for Rothrock (1975) pressure formulation + krdg_partic_out, & ! 0 for Thorndike et al. (1975) formulation + ! 1 for exponential participation function + krdg_redist_out ! 0 for Hibler (1980) formulation + ! 1 for exponential redistribution function + + real (kind=dbl_kind), intent(out), optional :: & + mu_rdg_out ! gives e-folding scale of ridged ice (m^.5) + ! (krdg_redist = 1) + +!----------------------------------------------------------------------- +! Parameters for atmosphere +!----------------------------------------------------------------------- + + real (kind=dbl_kind), intent(out), optional :: & + cp_air_out, & ! specific heat of air (J/kg/K) + cp_wv_out, & ! specific heat of water vapor (J/kg/K) + zvir_out, & ! rh2o/rair - 1.0 + zref_out, & ! reference height for stability (m) + qqqice_out, & ! for qsat over ice + TTTice_out, & ! for qsat over ice + qqqocn_out, & ! for qsat over ocn + TTTocn_out ! for qsat over ocn + + character (len=char_len), intent(out), optional :: & + atmbndy_out ! atmo boundary method, 'default' ('ccsm3') or 'constant' + + logical (kind=log_kind), intent(out), optional :: & + calc_strair_out, & ! if true, calculate wind stress components + formdrag_out, & ! if true, calculate form drag + highfreq_out ! if true, use high frequency coupling + + integer (kind=int_kind), intent(out), optional :: & + natmiter_out ! number of iterations for boundary layer calculations + + ! Flux convergence tolerance + real (kind=dbl_kind), intent(out), optional :: atmiter_conv_out + +!----------------------------------------------------------------------- +! Parameters for the ice thickness distribution +!----------------------------------------------------------------------- + + integer (kind=int_kind), intent(out), optional :: & + kitd_out , & ! type of itd conversions + ! 0 = delta function + ! 1 = linear remap + kcatbound_out ! 0 = old category boundary formula + ! 1 = new formula giving round numbers + ! 2 = WMO standard + ! 3 = asymptotic formula + +!----------------------------------------------------------------------- +! Parameters for the floe size distribution +!----------------------------------------------------------------------- + + integer (kind=int_kind), intent(out), optional :: & + nfreq_out ! number of frequencies + + real (kind=dbl_kind), intent(out), optional :: & + floeshape_out ! constant from Steele (unitless) + + logical (kind=log_kind), intent(out), optional :: & + wave_spec_out ! if true, use wave forcing + + character (len=char_len), intent(out), optional :: & + wave_spec_type_out ! type of wave spectrum forcing + +!----------------------------------------------------------------------- +! Parameters for biogeochemistry +!----------------------------------------------------------------------- + + character(char_len), intent(out), optional :: & + bgc_flux_type_out ! type of ocean-ice piston velocity + ! 'constant', 'Jin2006' + + logical (kind=log_kind), intent(out), optional :: & + z_tracers_out, & ! if .true., bgc or aerosol tracers are vertically resolved + scale_bgc_out, & ! if .true., initialize bgc tracers proportionally with salinity + solve_zbgc_out, & ! if .true., solve vertical biochemistry portion of code + dEdd_algae_out, & ! if .true., algal absorptionof Shortwave is computed in the + modal_aero_out, & ! if .true., use modal aerosol formulation in shortwave + conserv_check_out ! if .true., run conservation checks and abort if checks fail + + logical (kind=log_kind), intent(out), optional :: & + skl_bgc_out, & ! if true, solve skeletal biochemistry + solve_zsal_out ! if true, update salinity profile from solve_S_dt + + real (kind=dbl_kind), intent(out), optional :: & + grid_o_out , & ! for bottom flux + l_sk_out , & ! characteristic diffusive scale (zsalinity) (m) + initbio_frac_out, & ! fraction of ocean tracer concentration used to initialize tracer + phi_snow_out ! snow porosity at the ice/snow interface + + real (kind=dbl_kind), intent(out), optional :: & + grid_oS_out , & ! for bottom flux (zsalinity) + l_skS_out ! 0.02 characteristic skeletal layer thickness (m) (zsalinity) + real (kind=dbl_kind), intent(out), optional :: & + fr_resp_out , & ! fraction of algal growth lost due to respiration + algal_vel_out , & ! 0.5 cm/d(m/s) Lavoie 2005 1.5 cm/day + R_dFe2dust_out , & ! g/g (3.5% content) Tagliabue 2009 + dustFe_sol_out , & ! solubility fraction + T_max_out , & ! maximum temperature (C) + fsal_out , & ! Salinity limitation (ppt) + op_dep_min_out , & ! Light attenuates for optical depths exceeding min + fr_graze_s_out , & ! fraction of grazing spilled or slopped + fr_graze_e_out , & ! fraction of assimilation excreted + fr_mort2min_out , & ! fractionation of mortality to Am + fr_dFe_out , & ! fraction of remineralized nitrogen + ! (in units of algal iron) + k_nitrif_out , & ! nitrification rate (1/day) + t_iron_conv_out , & ! desorption loss pFe to dFe (day) + max_loss_out , & ! restrict uptake to % of remaining value + max_dfe_doc1_out , & ! max ratio of dFe to saccharides in the ice + ! (nM Fe/muM C) + fr_resp_s_out , & ! DMSPd fraction of respiration loss as DMSPd + y_sk_DMS_out , & ! fraction conversion given high yield + t_sk_conv_out , & ! Stefels conversion time (d) + t_sk_ox_out , & ! DMS oxidation time (d) + frazil_scav_out ! scavenging fraction or multiple in frazil ice + + real (kind=dbl_kind), intent(out), optional :: & + sk_l_out, & ! skeletal layer thickness (m) + min_bgc_out ! fraction of ocean bgc concentration in surface melt + +!----------------------------------------------------------------------- +! Parameters for melt ponds +!----------------------------------------------------------------------- + + real (kind=dbl_kind), intent(out), optional :: & + hs0_out ! snow depth for transition to bare sea ice (m) + + ! level-ice ponds + character (len=char_len), intent(out), optional :: & + frzpnd_out ! pond refreezing parameterization + + real (kind=dbl_kind), intent(out), optional :: & + dpscale_out, & ! alter e-folding time scale for flushing + rfracmin_out, & ! minimum retained fraction of meltwater + rfracmax_out, & ! maximum retained fraction of meltwater + pndaspect_out, & ! ratio of pond depth to pond fraction + hs1_out ! tapering parameter for snow on pond ice + + ! topo ponds + real (kind=dbl_kind), intent(out), optional :: & + hp1_out ! critical parameter for pond ice thickness + +!autodocument_end + + character(len=*),parameter :: subname='(icepack_query_parameters)' + + if (present(puny_out) ) puny_out = puny + if (present(bignum_out) ) bignum_out = bignum + if (present(pi_out) ) pi_out = pi + + if (present(c0_out) ) c0_out = c0 + if (present(c1_out) ) c1_out = c1 + if (present(c1p5_out) ) c1p5_out = c1p5 + if (present(c2_out) ) c2_out = c2 + if (present(c3_out) ) c3_out = c3 + if (present(c4_out) ) c4_out = c4 + if (present(c5_out) ) c5_out = c5 + if (present(c6_out) ) c6_out = c6 + if (present(c8_out) ) c8_out = c8 + if (present(c10_out) ) c10_out = c10 + if (present(c15_out) ) c15_out = c15 + if (present(c16_out) ) c16_out = c16 + if (present(c20_out) ) c20_out = c20 + if (present(c25_out) ) c25_out = c25 + if (present(c100_out) ) c100_out = c100 + if (present(c180_out) ) c180_out = c180 + if (present(c1000_out) ) c1000_out = c1000 + if (present(p001_out) ) p001_out = p001 + if (present(p01_out) ) p01_out = p01 + if (present(p1_out) ) p1_out = p1 + if (present(p2_out) ) p2_out = p2 + if (present(p4_out) ) p4_out = p4 + if (present(p5_out) ) p5_out = p5 + if (present(p6_out) ) p6_out = p6 + if (present(p05_out) ) p05_out = p05 + if (present(p15_out) ) p15_out = p15 + if (present(p25_out) ) p25_out = p25 + if (present(p75_out) ) p75_out = p75 + if (present(p333_out) ) p333_out = p333 + if (present(p666_out) ) p666_out = p666 + if (present(spval_const_out) ) spval_const_out = spval_const + if (present(pih_out) ) pih_out = pih + if (present(piq_out) ) piq_out = piq + if (present(pi2_out) ) pi2_out = pi2 + if (present(secday_out) ) secday_out = secday + if (present(rad_to_deg_out) ) rad_to_deg_out = rad_to_deg + + if (present(rhos_out) ) rhos_out = rhos + if (present(rhoi_out) ) rhoi_out = rhoi + if (present(rhow_out) ) rhow_out = rhow + if (present(cp_air_out) ) cp_air_out = cp_air + if (present(emissivity_out) ) emissivity_out = emissivity + if (present(floediam_out) ) floediam_out = floediam + if (present(hfrazilmin_out) ) hfrazilmin_out = hfrazilmin + if (present(cp_ice_out) ) cp_ice_out = cp_ice + if (present(cp_ocn_out) ) cp_ocn_out = cp_ocn + if (present(depressT_out) ) depressT_out = depressT + if (present(dragio_out) ) dragio_out = dragio + if (present(albocn_out) ) albocn_out = albocn + if (present(gravit_out) ) gravit_out = gravit + if (present(viscosity_dyn_out) ) viscosity_dyn_out= viscosity_dyn + if (present(Tocnfrz_out) ) Tocnfrz_out = Tocnfrz + if (present(rhofresh_out) ) rhofresh_out = rhofresh + if (present(zvir_out) ) zvir_out = zvir + if (present(vonkar_out) ) vonkar_out = vonkar + if (present(cp_wv_out) ) cp_wv_out = cp_wv + if (present(stefan_boltzmann_out) ) stefan_boltzmann_out = stefan_boltzmann + if (present(Tffresh_out) ) Tffresh_out = Tffresh + if (present(Lsub_out) ) Lsub_out = Lsub + if (present(Lvap_out) ) Lvap_out = Lvap + if (present(Timelt_out) ) Timelt_out = Timelt + if (present(Tsmelt_out) ) Tsmelt_out = Tsmelt + if (present(ice_ref_salinity_out) ) ice_ref_salinity_out = ice_ref_salinity + if (present(iceruf_out) ) iceruf_out = iceruf + if (present(Cf_out) ) Cf_out = Cf + if (present(Pstar_out) ) Pstar_out = Pstar + if (present(Cstar_out) ) Cstar_out = Cstar + if (present(kappav_out) ) kappav_out = kappav + if (present(kice_out) ) kice_out = kice + if (present(kseaice_out) ) kseaice_out = kseaice + if (present(ksno_out) ) ksno_out = ksno + if (present(zref_out) ) zref_out = zref + if (present(hs_min_out) ) hs_min_out = hs_min + if (present(snowpatch_out) ) snowpatch_out = snowpatch + if (present(rhosi_out) ) rhosi_out = rhosi + if (present(sk_l_out) ) sk_l_out = sk_l + if (present(saltmax_out) ) saltmax_out = saltmax + if (present(phi_init_out) ) phi_init_out = phi_init + if (present(min_salin_out) ) min_salin_out = min_salin + if (present(salt_loss_out) ) salt_loss_out = salt_loss + if (present(min_bgc_out) ) min_bgc_out = min_bgc + if (present(dSin0_frazil_out) ) dSin0_frazil_out = dSin0_frazil + if (present(hi_ssl_out) ) hi_ssl_out = hi_ssl + if (present(hs_ssl_out) ) hs_ssl_out = hs_ssl + if (present(awtvdr_out) ) awtvdr_out = awtvdr + if (present(awtidr_out) ) awtidr_out = awtidr + if (present(awtvdf_out) ) awtvdf_out = awtvdf + if (present(awtidf_out) ) awtidf_out = awtidf + if (present(qqqice_out) ) qqqice_out = qqqice + if (present(TTTice_out) ) TTTice_out = TTTice + if (present(qqqocn_out) ) qqqocn_out = qqqocn + if (present(TTTocn_out) ) TTTocn_out = TTTocn + if (present(puny_out) ) puny_out = puny + if (present(bignum_out) ) bignum_out = bignum + if (present(pi_out) ) pi_out = pi + if (present(secday_out) ) secday_out = secday + if (present(ktherm_out) ) ktherm_out = ktherm + if (present(conduct_out) ) conduct_out = conduct + if (present(fbot_xfer_type_out) ) fbot_xfer_type_out = fbot_xfer_type + if (present(heat_capacity_out) ) heat_capacity_out= heat_capacity + if (present(calc_Tsfc_out) ) calc_Tsfc_out = calc_Tsfc + if (present(update_ocn_f_out) ) update_ocn_f_out = update_ocn_f + if (present(dts_b_out) ) dts_b_out = dts_b + if (present(ustar_min_out) ) ustar_min_out = ustar_min + if (present(a_rapid_mode_out) ) a_rapid_mode_out = a_rapid_mode + if (present(Rac_rapid_mode_out) ) Rac_rapid_mode_out = Rac_rapid_mode + if (present(aspect_rapid_mode_out) ) aspect_rapid_mode_out = aspect_rapid_mode + if (present(dSdt_slow_mode_out) ) dSdt_slow_mode_out = dSdt_slow_mode + if (present(phi_c_slow_mode_out) ) phi_c_slow_mode_out = phi_c_slow_mode + if (present(phi_i_mushy_out) ) phi_i_mushy_out = phi_i_mushy + if (present(shortwave_out) ) shortwave_out = shortwave + if (present(albedo_type_out) ) albedo_type_out = albedo_type + if (present(albicev_out) ) albicev_out = albicev + if (present(albicei_out) ) albicei_out = albicei + if (present(albsnowv_out) ) albsnowv_out = albsnowv + if (present(albsnowi_out) ) albsnowi_out = albsnowi + if (present(ahmax_out) ) ahmax_out = ahmax + if (present(R_ice_out) ) R_ice_out = R_ice + if (present(R_pnd_out) ) R_pnd_out = R_pnd + if (present(R_snw_out) ) R_snw_out = R_snw + if (present(dT_mlt_out) ) dT_mlt_out = dT_mlt + if (present(rsnw_mlt_out) ) rsnw_mlt_out = rsnw_mlt + if (present(kalg_out) ) kalg_out = kalg + if (present(kstrength_out) ) kstrength_out = kstrength + if (present(krdg_partic_out) ) krdg_partic_out = krdg_partic + if (present(krdg_redist_out) ) krdg_redist_out = krdg_redist + if (present(mu_rdg_out) ) mu_rdg_out = mu_rdg + if (present(atmbndy_out) ) atmbndy_out = atmbndy + if (present(calc_strair_out) ) calc_strair_out = calc_strair + if (present(formdrag_out) ) formdrag_out = formdrag + if (present(highfreq_out) ) highfreq_out = highfreq + if (present(natmiter_out) ) natmiter_out = natmiter + if (present(atmiter_conv_out) ) atmiter_conv_out = atmiter_conv + if (present(tfrz_option_out) ) tfrz_option_out = tfrz_option + if (present(kitd_out) ) kitd_out = kitd + if (present(kcatbound_out) ) kcatbound_out = kcatbound + if (present(floeshape_out) ) floeshape_out = floeshape + if (present(wave_spec_out) ) wave_spec_out = wave_spec + if (present(wave_spec_type_out) ) wave_spec_type_out = wave_spec_type + if (present(nfreq_out) ) nfreq_out = nfreq + if (present(hs0_out) ) hs0_out = hs0 + if (present(frzpnd_out) ) frzpnd_out = frzpnd + if (present(dpscale_out) ) dpscale_out = dpscale + if (present(rfracmin_out) ) rfracmin_out = rfracmin + if (present(rfracmax_out) ) rfracmax_out = rfracmax + if (present(pndaspect_out) ) pndaspect_out = pndaspect + if (present(hs1_out) ) hs1_out = hs1 + if (present(hp1_out) ) hp1_out = hp1 + if (present(bgc_flux_type_out) ) bgc_flux_type_out= bgc_flux_type + if (present(z_tracers_out) ) z_tracers_out = z_tracers + if (present(scale_bgc_out) ) scale_bgc_out = scale_bgc + if (present(solve_zbgc_out) ) solve_zbgc_out = solve_zbgc + if (present(dEdd_algae_out) ) dEdd_algae_out = dEdd_algae + if (present(modal_aero_out) ) modal_aero_out = modal_aero + if (present(conserv_check_out) ) conserv_check_out= conserv_check + if (present(skl_bgc_out) ) skl_bgc_out = skl_bgc + if (present(solve_zsal_out) ) solve_zsal_out = solve_zsal + if (present(grid_o_out) ) grid_o_out = grid_o + if (present(l_sk_out) ) l_sk_out = l_sk + if (present(initbio_frac_out) ) initbio_frac_out = initbio_frac + if (present(grid_oS_out) ) grid_oS_out = grid_oS + if (present(l_skS_out) ) l_skS_out = l_skS + if (present(phi_snow_out) ) phi_snow_out = phi_snow + if (present(fr_resp_out) ) fr_resp_out = fr_resp + if (present(algal_vel_out) ) algal_vel_out = algal_vel + if (present(R_dFe2dust_out) ) R_dFe2dust_out = R_dFe2dust + if (present(dustFe_sol_out) ) dustFe_sol_out = dustFe_sol + if (present(T_max_out) ) T_max_out = T_max + if (present(fsal_out) ) fsal_out = fsal + if (present(op_dep_min_out) ) op_dep_min_out = op_dep_min + if (present(fr_graze_s_out) ) fr_graze_s_out = fr_graze_s + if (present(fr_graze_e_out) ) fr_graze_e_out = fr_graze_e + if (present(fr_mort2min_out) ) fr_mort2min_out = fr_mort2min + if (present(fr_dFe_out) ) fr_dFe_out = fr_dFe + if (present(k_nitrif_out) ) k_nitrif_out = k_nitrif + if (present(t_iron_conv_out) ) t_iron_conv_out = t_iron_conv + if (present(max_loss_out) ) max_loss_out = max_loss + if (present(max_dfe_doc1_out) ) max_dfe_doc1_out = max_dfe_doc1 + if (present(fr_resp_s_out) ) fr_resp_s_out = fr_resp_s + if (present(y_sk_DMS_out) ) y_sk_DMS_out = y_sk_DMS + if (present(t_sk_conv_out) ) t_sk_conv_out = t_sk_conv + if (present(t_sk_ox_out) ) t_sk_ox_out = t_sk_ox + if (present(frazil_scav_out) ) frazil_scav_out = frazil_scav + if (present(Lfresh_out) ) Lfresh_out = Lfresh + if (present(cprho_out) ) cprho_out = cprho + if (present(Cp_out) ) Cp_out = Cp + if (present(sw_redist_out) ) sw_redist_out = sw_redist + if (present(sw_frac_out) ) sw_frac_out = sw_frac + if (present(sw_dtemp_out) ) sw_dtemp_out = sw_dtemp + + call icepack_recompute_constants() + if (icepack_warnings_aborted(subname)) return + + end subroutine icepack_query_parameters + +!======================================================================= + +!autodocument_start icepack_write_parameters +! subroutine to write the column package internal parameters + + subroutine icepack_write_parameters(iounit) + + integer (kind=int_kind), intent(in) :: & + iounit ! unit number for output + +!autodocument_end + + character(len=*),parameter :: subname='(icepack_write_parameters)' + + write(iounit,*) subname + write(iounit,*) " rhos = ",rhos + write(iounit,*) " rhoi = ",rhoi + write(iounit,*) " rhow = ",rhow + write(iounit,*) " cp_air = ",cp_air + write(iounit,*) " emissivity = ",emissivity + write(iounit,*) " floediam = ",floediam + write(iounit,*) " hfrazilmin = ",hfrazilmin + write(iounit,*) " cp_ice = ",cp_ice + write(iounit,*) " cp_ocn = ",cp_ocn + write(iounit,*) " depressT = ",depressT + write(iounit,*) " dragio = ",dragio + write(iounit,*) " albocn = ",albocn + write(iounit,*) " gravit = ",gravit + write(iounit,*) " viscosity_dyn = ",viscosity_dyn + write(iounit,*) " Tocnfrz = ",Tocnfrz + write(iounit,*) " rhofresh = ",rhofresh + write(iounit,*) " zvir = ",zvir + write(iounit,*) " vonkar = ",vonkar + write(iounit,*) " cp_wv = ",cp_wv + write(iounit,*) " stefan_boltzmann = ",stefan_boltzmann + write(iounit,*) " Tffresh = ",Tffresh + write(iounit,*) " Lsub = ",Lsub + write(iounit,*) " Lvap = ",Lvap + write(iounit,*) " Timelt = ",Timelt + write(iounit,*) " Tsmelt = ",Tsmelt + write(iounit,*) " ice_ref_salinity = ",ice_ref_salinity + write(iounit,*) " iceruf = ",iceruf + write(iounit,*) " Cf = ",Cf + write(iounit,*) " Pstar = ",Pstar + write(iounit,*) " Cstar = ",Cstar + write(iounit,*) " kappav = ",kappav + write(iounit,*) " kice = ",kice + write(iounit,*) " kseaice = ",kseaice + write(iounit,*) " ksno = ",ksno + write(iounit,*) " zref = ",zref + write(iounit,*) " hs_min = ",hs_min + write(iounit,*) " snowpatch = ",snowpatch + write(iounit,*) " rhosi = ",rhosi + write(iounit,*) " sk_l = ",sk_l + write(iounit,*) " saltmax = ",saltmax + write(iounit,*) " phi_init = ",phi_init + write(iounit,*) " min_salin = ",min_salin + write(iounit,*) " salt_loss = ",salt_loss + write(iounit,*) " min_bgc = ",min_bgc + write(iounit,*) " dSin0_frazil = ",dSin0_frazil + write(iounit,*) " hi_ssl = ",hi_ssl + write(iounit,*) " hs_ssl = ",hs_ssl + write(iounit,*) " awtvdr = ",awtvdr + write(iounit,*) " awtidr = ",awtidr + write(iounit,*) " awtvdf = ",awtvdf + write(iounit,*) " awtidf = ",awtidf + write(iounit,*) " qqqice = ",qqqice + write(iounit,*) " TTTice = ",TTTice + write(iounit,*) " qqqocn = ",qqqocn + write(iounit,*) " TTTocn = ",TTTocn + write(iounit,*) " puny = ",puny + write(iounit,*) " bignum = ",bignum + write(iounit,*) " secday = ",secday + write(iounit,*) " pi = ",pi + write(iounit,*) " pih = ",pih + write(iounit,*) " piq = ",piq + write(iounit,*) " pi2 = ",pi2 + write(iounit,*) " rad_to_deg = ",rad_to_deg + write(iounit,*) " Lfresh = ",Lfresh + write(iounit,*) " cprho = ",cprho + write(iounit,*) " Cp = ",Cp + write(iounit,*) " ktherm = ", ktherm + write(iounit,*) " conduct = ", conduct + write(iounit,*) " fbot_xfer_type = ", fbot_xfer_type + write(iounit,*) " heat_capacity = ", heat_capacity + write(iounit,*) " calc_Tsfc = ", calc_Tsfc + write(iounit,*) " update_ocn_f = ", update_ocn_f + write(iounit,*) " dts_b = ", dts_b + write(iounit,*) " ustar_min = ", ustar_min + write(iounit,*) " a_rapid_mode = ", a_rapid_mode + write(iounit,*) " Rac_rapid_mode = ", Rac_rapid_mode + write(iounit,*) " aspect_rapid_mode = ", aspect_rapid_mode + write(iounit,*) " dSdt_slow_mode = ", dSdt_slow_mode + write(iounit,*) " phi_c_slow_mode = ", phi_c_slow_mode + write(iounit,*) " phi_i_mushy = ", phi_i_mushy + write(iounit,*) " shortwave = ", shortwave + write(iounit,*) " albedo_type = ", albedo_type + write(iounit,*) " albicev = ", albicev + write(iounit,*) " albicei = ", albicei + write(iounit,*) " albsnowv = ", albsnowv + write(iounit,*) " albsnowi = ", albsnowi + write(iounit,*) " ahmax = ", ahmax + write(iounit,*) " R_ice = ", R_ice + write(iounit,*) " R_pnd = ", R_pnd + write(iounit,*) " R_snw = ", R_snw + write(iounit,*) " dT_mlt = ", dT_mlt + write(iounit,*) " rsnw_mlt = ", rsnw_mlt + write(iounit,*) " kalg = ", kalg + write(iounit,*) " kstrength = ", kstrength + write(iounit,*) " krdg_partic = ", krdg_partic + write(iounit,*) " krdg_redist = ", krdg_redist + write(iounit,*) " mu_rdg = ", mu_rdg + write(iounit,*) " atmbndy = ", atmbndy + write(iounit,*) " calc_strair = ", calc_strair + write(iounit,*) " formdrag = ", formdrag + write(iounit,*) " highfreq = ", highfreq + write(iounit,*) " natmiter = ", natmiter + write(iounit,*) " atmiter_conv = ", atmiter_conv + write(iounit,*) " tfrz_option = ", tfrz_option + write(iounit,*) " kitd = ", kitd + write(iounit,*) " kcatbound = ", kcatbound + write(iounit,*) " floeshape = ", floeshape + write(iounit,*) " wave_spec = ", wave_spec + write(iounit,*) " wave_spec_type= ", wave_spec_type + write(iounit,*) " nfreq = ", nfreq + write(iounit,*) " hs0 = ", hs0 + write(iounit,*) " frzpnd = ", frzpnd + write(iounit,*) " dpscale = ", dpscale + write(iounit,*) " rfracmin = ", rfracmin + write(iounit,*) " rfracmax = ", rfracmax + write(iounit,*) " pndaspect = ", pndaspect + write(iounit,*) " hs1 = ", hs1 + write(iounit,*) " hp1 = ", hp1 + write(iounit,*) " bgc_flux_type = ", bgc_flux_type + write(iounit,*) " z_tracers = ", z_tracers + write(iounit,*) " scale_bgc = ", scale_bgc + write(iounit,*) " solve_zbgc = ", solve_zbgc + write(iounit,*) " dEdd_algae = ", dEdd_algae + write(iounit,*) " modal_aero = ", modal_aero + write(iounit,*) " conserv_check = ", conserv_check + write(iounit,*) " skl_bgc = ", skl_bgc + write(iounit,*) " solve_zsal = ", solve_zsal + write(iounit,*) " grid_o = ", grid_o + write(iounit,*) " l_sk = ", l_sk + write(iounit,*) " initbio_frac = ", initbio_frac + write(iounit,*) " grid_oS = ", grid_oS + write(iounit,*) " l_skS = ", l_skS + write(iounit,*) " phi_snow = ", phi_snow + write(iounit,*) " fr_resp = ", fr_resp + write(iounit,*) " algal_vel = ", algal_vel + write(iounit,*) " R_dFe2dust = ", R_dFe2dust + write(iounit,*) " dustFe_sol = ", dustFe_sol + write(iounit,*) " T_max = ", T_max + write(iounit,*) " fsal = ", fsal + write(iounit,*) " op_dep_min = ", op_dep_min + write(iounit,*) " fr_graze_s = ", fr_graze_s + write(iounit,*) " fr_graze_e = ", fr_graze_e + write(iounit,*) " fr_mort2min = ", fr_mort2min + write(iounit,*) " fr_dFe = ", fr_dFe + write(iounit,*) " k_nitrif = ", k_nitrif + write(iounit,*) " t_iron_conv = ", t_iron_conv + write(iounit,*) " max_loss = ", max_loss + write(iounit,*) " max_dfe_doc1 = ", max_dfe_doc1 + write(iounit,*) " fr_resp_s = ", fr_resp_s + write(iounit,*) " y_sk_DMS = ", y_sk_DMS + write(iounit,*) " t_sk_conv = ", t_sk_conv + write(iounit,*) " t_sk_ox = ", t_sk_ox + write(iounit,*) " frazil_scav = ", frazil_scav + write(iounit,*) " sw_redist = ", sw_redist + write(iounit,*) " sw_frac = ", sw_frac + write(iounit,*) " sw_dtemp = ", sw_dtemp + + end subroutine icepack_write_parameters + +!======================================================================= + +!autodocument_start icepack_recompute_constants +! subroutine to reinitialize some derived constants + + subroutine icepack_recompute_constants() + +!autodocument_end + + character(len=*),parameter :: subname='(icepack_recompute_constants)' + + cprho = cp_ocn*rhow + Lfresh = Lsub-Lvap + Cp = 0.5_dbl_kind*gravit*(rhow-rhoi)*rhoi/rhow + pih = p5*pi + piq = p5*p5*pi + pi2 = c2*pi + rad_to_deg = c180/pi + + end subroutine icepack_recompute_constants + +!======================================================================= + + end module icepack_parameters + +!======================================================================= 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..0d062ff17c --- /dev/null +++ b/config_src/external/Icepack_interfaces/icepack_tracers.F90 @@ -0,0 +1,190 @@ +!======================================================================= +! 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 + public :: icepack_init_tracer_sizes + public :: icepack_query_tracer_sizes + + integer (kind=int_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 + + subroutine icepack_init_tracer_sizes(& + ncat_in, nilyr_in, nslyr_in, nblyr_in, nfsd_in , & + n_algae_in, n_DOC_in, n_aero_in, n_iso_in, & + n_DON_in, n_DIC_in, n_fed_in, n_fep_in, n_zaero_in, & + ntrcr_in, ntrcr_o_in, nbtrcr_in, nbtrcr_sw_in) + + integer (kind=int_kind), intent(in), optional :: & + ncat_in , & ! Categories + nfsd_in , & ! + nilyr_in , & ! Layers + nslyr_in , & ! + nblyr_in , & ! + n_algae_in, & ! Dimensions + n_DOC_in , & ! + n_DON_in , & ! + n_DIC_in , & ! + n_fed_in , & ! + n_fep_in , & ! + n_zaero_in, & ! + n_iso_in , & ! + n_aero_in , & ! + ntrcr_in , & ! number of tracers in use + ntrcr_o_in, & ! number of non-bio tracers in use + nbtrcr_in , & ! number of bio tracers in use + nbtrcr_sw_in ! number of shortwave bio tracers in use + + end subroutine icepack_init_tracer_sizes + + subroutine icepack_query_tracer_sizes(& + max_algae_out , max_dic_out , max_doc_out , & + max_don_out , max_fe_out , nmodal1_out , & + nmodal2_out , max_aero_out , max_nbtrcr_out , & + ncat_out, nilyr_out, nslyr_out, nblyr_out, nfsd_out, & + n_algae_out, n_DOC_out, n_aero_out, n_iso_out, & + n_DON_out, n_DIC_out, n_fed_out, n_fep_out, n_zaero_out, & + ntrcr_out, ntrcr_o_out, nbtrcr_out, nbtrcr_sw_out) + + integer (kind=int_kind), intent(out), optional :: & + max_algae_out , & ! maximum number of algal types + max_dic_out , & ! maximum number of dissolved inorganic carbon types + max_doc_out , & ! maximum number of dissolved organic carbon types + max_don_out , & ! maximum number of dissolved organic nitrogen types + max_fe_out , & ! maximum number of iron types + nmodal1_out , & ! dimension for modal aerosol radiation parameters + nmodal2_out , & ! dimension for modal aerosol radiation parameters + max_aero_out , & ! maximum number of aerosols + max_nbtrcr_out ! algal nitrogen and chlorophyll + + integer (kind=int_kind), intent(out), optional :: & + ncat_out , & ! Categories + nfsd_out , & ! + nilyr_out , & ! Layers + nslyr_out , & ! + nblyr_out , & ! + n_algae_out, & ! Dimensions + n_DOC_out , & ! + n_DON_out , & ! + n_DIC_out , & ! + n_fed_out , & ! + n_fep_out , & ! + n_zaero_out, & ! + n_iso_out , & ! + n_aero_out , & ! + ntrcr_out , & ! number of tracers in use + ntrcr_o_out, & ! number of non-bio tracers in use + nbtrcr_out , & ! number of bio tracers in use + nbtrcr_sw_out ! number of shortwave bio tracers in use + + end subroutine icepack_query_tracer_sizes + + 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..d8e403a60b --- /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 c582e9ba94..dee419ab25 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..55ff9aa7a0 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,19 @@ 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) + ! 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 + ! 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 511695b8da..3dd8ebe85f 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 52f167664f..7696a51958 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 @@ -2136,6 +2137,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..9790c56ff7 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 - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! -! ice_ridging_init - initialize the ice ridging and set parameters. ! -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! - -!~~~~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 +! 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 - ! 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 +! e-folding scale of ridged ice, krdg_partic=1 (m^0.5) +real(kind=dbl_kind), parameter :: mu_rdg = 3.0 - 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,406 @@ 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 - endif - endif ; enddo + 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 - !------------------------------------------------------------------- - ! redistribution of ice - !------------------------------------------------------------------- + ! set category limits; Icepack has a max on the largest, unlimited, category (why?) - ! changes in open water area - cn(0) = max(cn(0) + (Rdiv - part_undef(1)*Rtot) * dt, 0.0) + 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 - if (Rtot>0.0) then + trcr_base = 0.0; n_trcr_strata = 0; nt_strata = 0; ! init some tracer vars + ! When would we use icepack tracer "strata"? - ! 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 + ! 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 - enddo + 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) + if (aice0<0.) then + call SIS_error(WARNING,'aice0<0. before call to ridge ice.') + aice0=0. + endif - ! 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 + 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 - 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 + ! pop pond off top of stack + tr_tmp(1:nCat)=trcrn(ntrcr,1:nCat) - !------------------------------------------------------------------- - if (.not. rdg_iterate) exit - enddo ! n_iterate - !------------------------------------------------------------------- + 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 - ! 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 + if (TrReg%ntr>0) then + ! unload tracer array reversing order of load -- stack-like fashion - ! calculate opening rate of ridging - rdg_open = (alev - ardg) / dt + 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 - ! 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 + 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 - ! return to true upper most ice thickness category limit - !hlim(km) = hlimtmp + 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 + + endif ! have tracers to unload + + ! ! output: snow/ice masses/thicknesses + do k=1,nCat + + 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 + + + IST%part_size(i,j,0) = 1.0 - sum(IST%part_size(i,j,1:nCat)) + + 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,:)); + + endif; enddo; enddo ! part_sz, j, i + +! 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