diff --git a/columnphysics/icepack_itd.F90 b/columnphysics/icepack_itd.F90 index 75f60aa95..f745c9915 100644 --- a/columnphysics/icepack_itd.F90 +++ b/columnphysics/icepack_itd.F90 @@ -31,10 +31,10 @@ module icepack_itd use icepack_parameters, only: rhosi, sk_l, hs_ssl, min_salin, rsnw_fall 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, tr_iso, tr_snow, nt_smice, nt_rsnw, nt_rhos + use icepack_tracers, only: n_iso, tr_iso, tr_snow, nt_smice, nt_rsnw, nt_rhos, nt_sice 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_parameters, only: kcatbound, kitd, saltflux_option 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 @@ -1089,7 +1089,7 @@ subroutine zap_small_areas (dt, ntrcr, & n, k, it, & !counting indices blevels - real (kind=dbl_kind) :: xtmp ! temporary variables + real (kind=dbl_kind) :: xtmp, sicen ! temporary variables real (kind=dbl_kind) , dimension (1):: trcr_skl real (kind=dbl_kind) , dimension (nblyr+1):: bvol @@ -1187,7 +1187,15 @@ subroutine zap_small_areas (dt, ntrcr, & xtmp = (rhoi*vicen(n)) / dt dfresh = dfresh + xtmp - xtmp = rhoi*vicen(n)*ice_ref_salinity*p001 / dt + if (saltflux_option == 'prognostic') then + sicen = c0 + do k=1,nilyr + sicen = sicen + trcrn(nt_sice+k-1,n) / real(nilyr,kind=dbl_kind) + enddo + xtmp = rhoi*vicen(n)*sicen*p001 / dt + else + xtmp = rhoi*vicen(n)*ice_ref_salinity*p001 / dt + endif dfsalt = dfsalt + xtmp aice0 = aice0 + aicen(n) @@ -1301,13 +1309,22 @@ subroutine zap_small_areas (dt, ntrcr, & !----------------------------------------------------------------- xtmp = (rhoi*vicen(n) + rhos*vsnon(n)) & - * (aice-c1)/aice / dt - dfresh = dfresh + xtmp - - xtmp = rhoi*vicen(n)*ice_ref_salinity*p001 & - * (aice-c1)/aice / dt - dfsalt = dfsalt + xtmp - + * (aice-c1)/aice / dt + dfresh = dfresh + xtmp + + if (saltflux_option == 'prognostic') then + sicen = c0 + do k=1,nilyr + sicen = sicen + trcrn(nt_sice+k-1,n) / real(nilyr,kind=dbl_kind) + enddo + xtmp = rhoi*vicen(n)*sicen*p001 & + * (aice-c1)/aice / dt + else + xtmp = rhoi*vicen(n)*ice_ref_salinity*p001 & + * (aice-c1)/aice / dt + endif + dfsalt = dfsalt + xtmp + if (solve_zsal) then do k = 1,nblyr xtmp = rhosi*trcrn(nt_fbri,n)*vicen(n)*p001& diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index b6b1f197d..43bd2c21a 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/columnphysics/icepack_parameters.F90 @@ -165,6 +165,11 @@ module icepack_parameters ! 'linear_salt' = -depressT * sss ! 'mushy' conforms with ktherm=2 + character(len=char_len), public :: & + saltflux_option = 'constant'! Salt flux computation + ! 'constant' reference value of ice_ref_salinity + ! 'prognostic' prognostic salt flux + !----------------------------------------------------------------------- ! Parameters for radiation !----------------------------------------------------------------------- @@ -454,6 +459,7 @@ subroutine icepack_init_parameters( & atmbndy_in, calc_strair_in, formdrag_in, highfreq_in, natmiter_in, & atmiter_conv_in, calc_dragio_in, & tfrz_option_in, kitd_in, kcatbound_in, hs0_in, frzpnd_in, & + saltflux_option_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, & @@ -557,12 +563,17 @@ subroutine icepack_init_parameters( & 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=*), intent(in), optional :: & + tfrz_option_in ! form of ocean freezing temperature + ! 'minus1p8' = -1.8 C + ! 'linear_salt' = -depressT * sss + ! 'mushy' conforms with ktherm=2 - character(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 + character(len=*), intent(in), optional :: & + saltflux_option_in ! Salt flux computation + ! 'constant' reference value of ice_ref_salinity + ! 'prognostic' prognostic salt flux !----------------------------------------------------------------------- ! Parameters for radiation @@ -917,6 +928,7 @@ subroutine icepack_init_parameters( & 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(saltflux_option_in) ) saltflux_option = saltflux_option_in if (present(kitd_in) ) kitd = kitd_in if (present(kcatbound_in) ) kcatbound = kcatbound_in if (present(floeshape_in) ) floeshape = floeshape_in @@ -1145,6 +1157,7 @@ subroutine icepack_query_parameters( & atmbndy_out, calc_strair_out, formdrag_out, highfreq_out, natmiter_out, & atmiter_conv_out, calc_dragio_out, & tfrz_option_out, kitd_out, kcatbound_out, hs0_out, frzpnd_out, & + saltflux_option_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, & @@ -1264,6 +1277,12 @@ subroutine icepack_query_parameters( & ! 'linear_salt' = -depressT * sss ! 'mushy' conforms with ktherm=2 + character(len=*), intent(out), optional :: & + saltflux_option_out ! Salt flux computation + ! 'constant' reference value of ice_ref_salinity + ! 'prognostic' prognostic salt flux + + !----------------------------------------------------------------------- ! Parameters for radiation !----------------------------------------------------------------------- @@ -1653,6 +1672,7 @@ subroutine icepack_query_parameters( & 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(saltflux_option_out) ) saltflux_option_out = saltflux_option if (present(kitd_out) ) kitd_out = kitd if (present(kcatbound_out) ) kcatbound_out = kcatbound if (present(floeshape_out) ) floeshape_out = floeshape @@ -1856,6 +1876,7 @@ subroutine icepack_write_parameters(iounit) write(iounit,*) " natmiter = ", natmiter write(iounit,*) " atmiter_conv = ", atmiter_conv write(iounit,*) " tfrz_option = ", tfrz_option + write(iounit,*) " saltflux_option = ", saltflux_option write(iounit,*) " kitd = ", kitd write(iounit,*) " kcatbound = ", kcatbound write(iounit,*) " floeshape = ", floeshape diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index d93a7c47a..db76fcaa7 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -27,11 +27,12 @@ module icepack_therm_itd use icepack_parameters, only: rhosi, conserv_check, rhosmin use icepack_parameters, only: kitd, ktherm use icepack_parameters, only: z_tracers, solve_zsal, hfrazilmin + use icepack_parameters, only: saltflux_option use icepack_tracers, only: ntrcr, nbtrcr use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice use icepack_tracers, only: nt_apnd, nt_hpnd, nt_aero, nt_isosno, nt_isoice - use icepack_tracers, only: nt_Tsfc, nt_iage, nt_FY, nt_fsd, nt_rhos + use icepack_tracers, only: nt_Tsfc, nt_iage, nt_FY, nt_fsd, nt_rhos, nt_sice use icepack_tracers, only: nt_alvl, nt_vlvl use icepack_tracers, only: tr_pond_lvl, tr_pond_topo, tr_snow use icepack_tracers, only: tr_iage, tr_FY, tr_lvl, tr_aero, tr_iso, tr_brine, tr_fsd @@ -984,7 +985,7 @@ subroutine lateral_melt (dt, ncat, & real (kind=dbl_kind), intent(in) :: & sss real (kind=dbl_kind) :: & - Ti, Si0, qi0, & + Ti, Si0, qi0, sicen, & elapsed_t, & ! FSD subcycling subdt ! FSD timestep (s) @@ -1093,7 +1094,15 @@ subroutine lateral_melt (dt, ncat, & ! dfresh > 0, dfsalt > 0, dfpond > 0 dfresh = (rhoi*vicen(n) + rhos*vsnon(n)) * rsiden(n) / dt - dfsalt = rhoi*vicen(n)*ice_ref_salinity*p001 * rsiden(n) / dt + if (saltflux_option == 'prognostic') then + sicen = c0 + do k=1,nilyr + sicen = sicen + trcrn(nt_sice+k-1,n) / real(nilyr,kind=dbl_kind) + enddo + dfsalt = rhoi*vicen(n)*sicen*p001 * rsiden(n) / dt + else + dfsalt = rhoi*vicen(n)*ice_ref_salinity*p001 * rsiden(n) / dt + endif fresh = fresh + dfresh fsalt = fsalt + dfsalt @@ -1596,18 +1605,23 @@ subroutine add_new_ice (ncat, nilyr, & !----------------------------------------------------------------- if (update_ocn_f) then - if (ktherm <= 1) then - dfresh = -rhoi*vi0new/dt + dfresh = -rhoi*vi0new/dt + if (saltflux_option == 'prognostic') then + dfsalt = Si0new*p001*dfresh + else dfsalt = ice_ref_salinity*p001*dfresh - fresh = fresh + dfresh - fsalt = fsalt + dfsalt - ! elseif (ktherm == 2) the fluxes are added elsewhere endif + fresh = fresh + dfresh + fsalt = fsalt + dfsalt else ! update_ocn_f = false if (ktherm == 2) then ! return mushy-layer frazil to ocean (POP) vi0tmp = fnew*dt / (rhoi*Lfresh) dfresh = -rhoi*(vi0new - vi0tmp)/dt - dfsalt = ice_ref_salinity*p001*dfresh + if (saltflux_option == 'prognostic') then + dfsalt = Si0new*p001*dfresh + else + dfsalt = ice_ref_salinity*p001*dfresh + endif fresh = fresh + dfresh fsalt = fsalt + dfsalt frazil_diag = frazil - vi0tmp diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 17e0c1fd1..0dba94f73 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -27,6 +27,7 @@ module icepack_therm_vertical use icepack_parameters, only: ustar_min, fbot_xfer_type, formdrag, calc_strair use icepack_parameters, only: rfracmin, rfracmax, dpscale, frzpnd, snwgrain, snwlvlfac use icepack_parameters, only: phi_i_mushy, floeshape, floediam, use_smliq_pnd, snwredist + use icepack_parameters, only: saltflux_option use icepack_tracers, only: tr_iage, tr_FY, tr_aero, tr_pond, tr_fsd, tr_iso use icepack_tracers, only: tr_pond_lvl, tr_pond_topo @@ -239,7 +240,7 @@ subroutine thermo_vertical (nilyr, nslyr, & einter ! intermediate energy real (kind=dbl_kind) :: & - fadvocn ! advective heat flux to ocean + fadvocn, saltvol, dfsalt ! advective heat flux to ocean character(len=*),parameter :: subname='(thermo_vertical)' @@ -292,6 +293,14 @@ subroutine thermo_vertical (nilyr, nslyr, & worki = hin works = hsn + ! Save initial salt volume for prognostic flux + if (saltflux_option == 'prognostic') then + saltvol = c0 + do k=1,nilyr + saltvol = saltvol + rhoi*zSin(k)*hin*p001 / real(nilyr,kind=dbl_kind) + enddo + endif + !----------------------------------------------------------------- ! Compute new surface temperature and internal ice and snow ! temperatures. @@ -435,8 +444,16 @@ subroutine thermo_vertical (nilyr, nslyr, & dhs = hsn - works - hsn_new freshn = freshn + evapn - (rhoi*dhi + rhos*dhs) / dt - fsaltn = fsaltn - rhoi*dhi*ice_ref_salinity*p001/dt - fhocnn = fhocnn + fadvocn ! for ktherm=2 + if (saltflux_option == 'prognostic') then + dfsalt = c0 + do k=1,nilyr + dfsalt = dfsalt + rhoi*zSin(k)*hin*p001 / real(nilyr,kind=dbl_kind) + enddo + fsaltn = fsaltn - (dfsalt - saltvol) / dt + else + fsaltn = fsaltn - rhoi*dhi*ice_ref_salinity*p001/dt + endif + fhocnn = fhocnn + fadvocn ! for ktherm=2 if (hin == c0) then if (tr_pond_topo) fpond = fpond - aicen * apond * hpond @@ -2569,7 +2586,6 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & enddo endif - !----------------------------------------------------------------- ! Update the neutral drag coefficients to account for form drag ! Oceanic and atmospheric drag coefficients diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90 index 200eba29e..75c8ac5d9 100644 --- a/configuration/driver/icedrv_init.F90 +++ b/configuration/driver/icedrv_init.F90 @@ -97,7 +97,7 @@ subroutine input_data natmiter, kitd, kcatbound character (len=char_len) :: shortwave, albedo_type, conduct, fbot_xfer_type, & - tfrz_option, frzpnd, atmbndy, wave_spec_type, snwredist, snw_aging_table + tfrz_option, saltflux_option, frzpnd, atmbndy, wave_spec_type, snwredist, snw_aging_table logical (kind=log_kind) :: sw_redist, use_smliq_pnd, snwgrain real (kind=dbl_kind) :: sw_frac, sw_dtemp @@ -105,6 +105,9 @@ subroutine input_data ! Flux convergence tolerance real (kind=dbl_kind) :: atmiter_conv + ! Ice reference salinity for fluxes + real (kind=dbl_kind) :: ice_ref_salinity + logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair, calc_dragio logical (kind=log_kind) :: conserv_check @@ -168,7 +171,8 @@ subroutine input_data fbot_xfer_type, oceanmixed_ice, emissivity, & formdrag, highfreq, natmiter, & atmiter_conv, calc_dragio, & - tfrz_option, default_season, wave_spec_type, & + tfrz_option, saltflux_option, ice_ref_salinity, & + default_season, wave_spec_type, & precip_units, fyear_init, ycycle, & atm_data_type, ocn_data_type, bgc_data_type, & atm_data_file, ocn_data_file, bgc_data_file, & @@ -216,7 +220,8 @@ subroutine input_data dSdt_slow_mode_out=dSdt_slow_mode, & phi_c_slow_mode_out=phi_c_slow_mode, & phi_i_mushy_out=phi_i_mushy, conserv_check_out=conserv_check, & - tfrz_option_out=tfrz_option, kalg_out=kalg, & + tfrz_option_out=tfrz_option, saltflux_option_out=saltflux_option, & + ice_ref_salinity_out=ice_ref_salinity, kalg_out=kalg, & fbot_xfer_type_out=fbot_xfer_type, puny_out=puny, & wave_spec_type_out=wave_spec_type, & sw_redist_out=sw_redist, sw_frac_out=sw_frac, sw_dtemp_out=sw_dtemp, & @@ -573,6 +578,13 @@ subroutine input_data 'WARNING: For consistency, set tfrz_option = mushy' endif + if (ktherm == 1 .and. trim(saltflux_option) /= 'constant') then + write (nu_diag,*) & + 'WARNING: ktherm = 1 and saltflux_option = ',trim(saltflux_option) + write (nu_diag,*) & + 'WARNING: For consistency, set saltflux_option = constant' + endif + if (ktherm == 0) then write (nu_diag,*) 'WARNING: ktherm = 0 zero-layer thermodynamics' write (nu_diag,*) 'WARNING: has been deprecated' @@ -761,6 +773,11 @@ subroutine input_data write(nu_diag,1010) ' oceanmixed_ice = ', oceanmixed_ice write(nu_diag,*) ' tfrz_option = ', & trim(tfrz_option) + write(nu_diag,*) ' saltflux_option = ', & + trim(saltflux_option) + if (trim(saltflux_option) == 'constant') then + write(nu_diag,1005) ' ice_ref_salinity = ', ice_ref_salinity + endif write(nu_diag,1010) ' restore_ocn = ', restore_ocn if (restore_ocn) & write(nu_diag,1005) ' trestore = ', trestore @@ -942,7 +959,8 @@ subroutine input_data dSdt_slow_mode_in=dSdt_slow_mode, & phi_c_slow_mode_in=phi_c_slow_mode, & phi_i_mushy_in=phi_i_mushy, conserv_check_in=conserv_check, & - tfrz_option_in=tfrz_option, kalg_in=kalg, & + tfrz_option_in=tfrz_option, saltflux_option_in=saltflux_option, & + ice_ref_salinity_in=ice_ref_salinity, kalg_in=kalg, & fbot_xfer_type_in=fbot_xfer_type, & wave_spec_type_in=wave_spec_type, wave_spec_in=wave_spec, & sw_redist_in=sw_redist, sw_frac_in=sw_frac, sw_dtemp_in=sw_dtemp, & diff --git a/configuration/scripts/icepack_in b/configuration/scripts/icepack_in index ce345648f..505dad4c6 100644 --- a/configuration/scripts/icepack_in +++ b/configuration/scripts/icepack_in @@ -109,6 +109,8 @@ update_ocn_f = .false. l_mpond_fresh = .false. tfrz_option = 'mushy' + ice_ref_salinity = 4.0 + saltflux_option = 'constant' oceanmixed_ice = .true. wave_spec_type = 'none' restore_ocn = .false. diff --git a/configuration/scripts/options/set_nml.saltflux b/configuration/scripts/options/set_nml.saltflux new file mode 100644 index 000000000..8fdb3940e --- /dev/null +++ b/configuration/scripts/options/set_nml.saltflux @@ -0,0 +1,2 @@ + ktherm = 2 + saltflux_option = 'prognostic' diff --git a/configuration/scripts/tests/base_suite.ts b/configuration/scripts/tests/base_suite.ts index 93f456524..771381355 100644 --- a/configuration/scripts/tests/base_suite.ts +++ b/configuration/scripts/tests/base_suite.ts @@ -31,6 +31,7 @@ restart col 1x1 alt01 restart col 1x1 alt02 restart col 1x1 alt03 restart col 1x1 alt04 +restart col 1x1 saltflux restart col 1x1 dyn restart col 1x1 fsd12,short restart col 1x1 snwitdrdg,snwgrain diff --git a/doc/source/icepack_index.rst b/doc/source/icepack_index.rst index 80be4abf6..439e4c74e 100755 --- a/doc/source/icepack_index.rst +++ b/doc/source/icepack_index.rst @@ -403,6 +403,7 @@ either Celsius or Kelvin units). Deprecated parameters are listed at the end. "rsnw_sig", "standard deviation of snow grain radius", "250. :math:`\times` 10\ :math:`^{-6}` m" "**S**", "", "" "salinz", "ice salinity profile", "ppt" + "saltflux_option", "constant or prognostic salinity fluxes","constant" "saltmax", "max salinity, at ice base (:cite:`Bitz99`)", "3.2 ppt" "scale_factor", "scaling factor for shortwave radiation components", "" "sec", "seconds elasped into idate", "" diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 837997a0a..85e75d4ce 100644 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -358,7 +358,9 @@ forcing_nml "", "``mm_per_sec``", "(same as MKS units)", "" "", "``m_per_sec``", "", "" "``restore_ocn``", "logical", "restore sst to data", "``.false.``" - "``tfrz_option``", "``linear_salt``", "linear functino of salinity (ktherm=1)", "``mushy``" + "``saltflux_option``", "``constant``","salt flux is referenced to a constant salinity","``constant``" + "","``prognostic``","use actual sea ice bulk salinity in flux" + "``tfrz_option``", "``linear_salt``", "linear function of salinity (ktherm=1)", "``mushy``" "", "``minus1p8``", "constant ocean freezing temperature (:math:`-1.8^{\circ} C`)", "" "", "``mushy``", "matches mushy-layer thermo (ktherm=2)", "" "``trestore``", "integer", "sst restoring time scale (days)", "90"