diff --git a/config_src/external/Icepack_interfaces/icepack_meltpond_sealvl.F90 b/config_src/external/Icepack_interfaces/icepack_meltpond_sealvl.F90 new file mode 100644 index 00000000..77f930af --- /dev/null +++ b/config_src/external/Icepack_interfaces/icepack_meltpond_sealvl.F90 @@ -0,0 +1,84 @@ +module icepack_meltpond_sealvl + + use icepack_kinds, only : int_kind, dbl_kind, char_len + use icepack_tracers, only : n_iso, n_aero + + implicit none + + private + public :: compute_ponds_sealvl, & + pond_hypsometry, & + pond_height + contains + +!> Interface for updating the melt ponds using Icepack. + subroutine compute_ponds_sealvl( dt, & + meltt, melts, frain, & + Tair, fsurfn, Tsfcn, & + dhs, ffrac, & + aicen, vicen, vsnon, & + qicen, sicen, & + apnd, hpnd, ipnd, & + meltsliqn, & + dpnd_freebdn, & + dpnd_dlidn, dpnd_flushn) + + real(kind=dbl_kind), intent(in) :: & + dt !< time step (s) + + real(kind=dbl_kind), intent(in) :: & + Tsfcn, & !< surface temperature (C) + meltt, & !< top melt rate (m/s) + melts, & !< snow melt rate (m/s) + frain, & !< rainfall rate (kg/m2/s) + Tair, & !< air temperature (K) + fsurfn,& !< atm-ice surface heat flux (W/m2) + aicen, & !< ice area fraction + vicen, & !< ice volume (m) + vsnon, & !< snow volume (m) + meltsliqn !< liquid contribution to meltponds in dt (kg/m^2) + + real(kind=dbl_kind), intent(inout) :: & + apnd, hpnd, ipnd, & !< pond tracers + dpnd_freebdn, & !< pond drainage rate due to freeboard constraint (m/step) + dpnd_dlidn, & !< pond loss/gain due to ice lid (m/step) + dpnd_flushn !< pond flushing rate due to ice permeability (m/s) + + real(kind=dbl_kind), dimension (:), intent(in) :: & + qicen, & !< ice layer enthalpy (J m-3) + sicen !< salinity (ppt) + + real(kind=dbl_kind), intent(in) :: & + dhs !< depth difference for snow on sea ice and pond ice + + real(kind=dbl_kind), intent(out) :: & + ffrac !< fraction of fsurfn over pond used to melt ipond + + end subroutine compute_ponds_sealvl + + subroutine pond_hypsometry(hpnd, apnd, dhpond, dvpond, hin) + + real(kind=dbl_kind), intent(inout) :: & + hpnd, & !< pond depth of ponded area tracer [m] + apnd !< pond fractional area of category tracer + + real(kind=dbl_kind), intent(in), optional :: & + dvpond, & !< incoming change in pond volume per category area + dhpond, & !< incoming change in pond depth [m] + hin !< category ice thickness [m] + + end subroutine pond_hypsometry + + subroutine pond_height(apond, hpnd, hin, hpsurf) + + real(kind=dbl_kind), intent(in) :: & + hin, & !< category mean ice thickness [m] + apond, & !< pond area fraction of the category + hpnd !< mean pond depth [m] + + real(kind=dbl_kind), intent(out) :: & + hpsurf !< height of pond surface above base of the ice [m] + + end subroutine pond_height + +end module icepack_meltpond_sealvl diff --git a/src/SIS2_ice_thm.F90 b/src/SIS2_ice_thm.F90 index 57677314..2e178401 100644 --- a/src/SIS2_ice_thm.F90 +++ b/src/SIS2_ice_thm.F90 @@ -1005,7 +1005,7 @@ end subroutine ice_check !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> ice_resize_SIS2 makes snow and ice thickness and temperature changes in ice !! columns due to thermodynamic forcing. -subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, & +subroutine ice_resize_SIS2(a_ice, m_pond, m_pond_ice, m_lay, Enthalpy, Sice_therm, Salin, & snow, rain, evap, tmlt, bmlt, NkIce, npassive, TrLay, & heat_to_ocn, h2o_ice_to_ocn, h2o_ocn_to_ice, evap_from_ocn, & snow_to_ice, salt_to_ice, ITV, US, CS, ablation, & @@ -1013,6 +1013,7 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, & ! mw/new - melt pond - added first two arguments & rain real, intent(in ) :: a_ice !< area of ice (1-open_water_frac) for pond retention [nondim] real, intent(inout) :: m_pond !< melt pond mass [R Z ~> kg m-2] + real, intent(inout) :: m_pond_ice !< melt pond ice mass [R Z ~> kg m-2] real, dimension(0:NkIce), & intent(inout) :: m_lay !< Snow and ice mass per unit area by layer [R Z ~> kg m-2]. real, dimension(0:NkIce+1), & diff --git a/src/SIS_continuity.F90 b/src/SIS_continuity.F90 index 6cfb4ad3..8797527e 100644 --- a/src/SIS_continuity.F90 +++ b/src/SIS_continuity.F90 @@ -617,7 +617,7 @@ end subroutine summed_continuity !! using input total mass fluxes with the fluxes proportionate to the relative upwind !! thicknesses. subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, US, IG, CS, & - h1, uh1, vh1, h2, uh2, vh2, h3, uh3, vh3) + h1, uh1, vh1, h2, uh2, vh2, h3, uh3, vh3, h4, uh4, vh4) 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, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_tot_in !< Initial total ice and snow mass per unit @@ -657,6 +657,15 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, US, IG, CS, real, dimension(SZI_(G),SZJB_(G),SZCAT_(IG)), & optional, intent(out) :: vh3 !< Meridional mass flux of medium 3 by category !! [R Z L2 T-1 ~> kg s-1]. + real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), & + optional, intent(inout) :: h4 !< Updated mass of medium 4 (pond ice?) by + !! category [R Z ~> kg m-2]. + real, dimension(SZIB_(G),SZJ_(G),SZCAT_(IG)), & + optional, intent(out) :: uh4 !< Zonal mass flux of medium 4 by category + !! [R Z L2 T-1 ~> kg s-1]. + real, dimension(SZI_(G),SZJB_(G),SZCAT_(IG)), & + optional, intent(out) :: vh4 !< Meridional mass flux of medium 4 by category + !! [R Z L2 T-1 ~> kg s-1]. ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: h_tot ! Total thicknesses [R Z ~> kg m-2]. @@ -730,6 +739,20 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, US, IG, CS, ((uh3(I,j,k) - uh3(I-1,j,k)) + (vh3(i,J,k) - vh3(i,J-1,k)))) enddo ; enddo ; enddo endif + if (present(h4)) then + if (do_mask) then + call zonal_proportionate_fluxes(uh_tot, I_htot, h4, uh4, G, IG, LB, masking_uh=uh1) + call merid_proportionate_fluxes(vh_tot, I_htot, h4, vh4, G, IG, LB, masking_vh=vh1) + else + call zonal_proportionate_fluxes(uh_tot, I_htot, h4, uh4, G, IG, LB) + call merid_proportionate_fluxes(vh_tot, I_htot, h4, vh4, G, IG, LB) + endif + !$OMP parallel do default(shared) + do j=js,je ; do k=1,nCat ; do i=is,ie + h4(i,j,k) = h4(i,j,k) - G%IareaT(i,j) * (dt * & + ((uh4(I,j,k) - uh4(I-1,j,k)) + (vh4(i,J,k) - vh4(i,J-1,k)))) + enddo ; enddo ; enddo + endif elseif (x_first) then ! First, advect zonally. @@ -767,6 +790,17 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, US, IG, CS, h3(i,j,k) = h3(i,j,k) - G%IareaT(i,j) * (dt * (uh3(I,j,k) - uh3(I-1,j,k))) enddo ; enddo ; enddo endif + if (present(h4)) then + if (do_mask) then + call zonal_proportionate_fluxes(uh_tot, I_htot, h4, uh4, G, IG, LB, masking_uh=uh1) + else + call zonal_proportionate_fluxes(uh_tot, I_htot, h4, uh4, G, IG, LB) + endif + !$OMP parallel do default(shared) + do j=LB%jsh,LB%jeh ; do k=1,nCat ; do i=LB%ish,LB%ieh + h4(i,j,k) = h4(i,j,k) - G%IareaT(i,j) * (dt * (uh4(I,j,k) - uh4(I-1,j,k))) + enddo ; enddo ; enddo + endif !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh @@ -811,6 +845,17 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, US, IG, CS, h3(i,j,k) = h3(i,j,k) - G%IareaT(i,j) * (dt * (vh3(i,J,k) - vh3(i,J-1,k)) ) enddo ; enddo ; enddo endif + if (present(h4)) then + if (do_mask) then + call merid_proportionate_fluxes(vh_tot, I_htot, h4, vh4, G, IG, LB, masking_vh=vh1) + else + call merid_proportionate_fluxes(vh_tot, I_htot, h4, vh4, G, IG, LB) + endif + !$OMP parallel do default(shared) + do j=js,je ; do k=1,nCat ; do i=is,ie + h4(i,j,k) = h4(i,j,k) - G%IareaT(i,j) * (dt * (vh4(i,J,k) - vh4(i,J-1,k)) ) + enddo ; enddo ; enddo + endif !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh @@ -857,6 +902,17 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, US, IG, CS, h3(i,j,k) = h3(i,j,k) - G%IareaT(i,j) * (dt * (vh3(i,J,k) - vh3(i,J-1,k)) ) enddo ; enddo ; enddo endif + if (present(h4)) then + if (do_mask) then + call merid_proportionate_fluxes(vh_tot, I_htot, h4, vh4, G, IG, LB, frac_neglect=CS%frac_neglect) + else + call merid_proportionate_fluxes(vh_tot, I_htot, h4, vh4, G, IG, LB) + endif + !$OMP parallel do default(shared) + do j=js,je ; do k=1,nCat ; do i=is,ie + h4(i,j,k) = h4(i,j,k) - G%IareaT(i,j) * (dt * (vh4(i,J,k) - vh4(i,J-1,k)) ) + enddo ; enddo ; enddo + endif !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh @@ -901,6 +957,17 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, US, IG, CS, h3(i,j,k) = h3(i,j,k) - G%IareaT(i,j) * (dt * (uh3(I,j,k) - uh3(I-1,j,k))) enddo ; enddo ; enddo endif + if (present(h4)) then + if (do_mask) then + call zonal_proportionate_fluxes(uh_tot, I_htot, h4, uh4, G, IG, LB, masking_uh=uh1) + else + call zonal_proportionate_fluxes(uh_tot, I_htot, h4, uh4, G, IG, LB) + endif + !$OMP parallel do default(shared) + do j=LB%jsh,LB%jeh ; do k=1,nCat ; do i=LB%ish,LB%ieh + h4(i,j,k) = h4(i,j,k) - G%IareaT(i,j) * (dt * (uh4(I,j,k) - uh4(I-1,j,k))) + enddo ; enddo ; enddo + endif !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh diff --git a/src/SIS_dyn_trans.F90 b/src/SIS_dyn_trans.F90 index e98a90df..347df0b8 100644 --- a/src/SIS_dyn_trans.F90 +++ b/src/SIS_dyn_trans.F90 @@ -425,7 +425,7 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, U !$OMP parallel do default(shared) do j=jsd,jed ; do k=1,ncat ; do i=isd,ied misp_sum(i,j) = misp_sum(i,j) + IST%part_size(i,j,k) * & - (IST%mH_snow(i,j,k) + IST%mH_pond(i,j,k)) + ((IST%mH_snow(i,j,k) + IST%mH_pond(i,j,k)) + IST%mH_pond_ice(i,j,k)) mi_sum(i,j) = mi_sum(i,j) + IST%mH_ice(i,j,k) * IST%part_size(i,j,k) ice_cover(i,j) = ice_cover(i,j) + IST%part_size(i,j,k) enddo ; enddo ; enddo @@ -873,7 +873,7 @@ subroutine convert_IST_to_simple_state(IST, DS2d, CAS, G, US, IG, CS) !$OMP parallel do default(shared) do j=jsd,jed ; do k=1,IG%CatIce ; do i=isd,ied DS2d%mca_step(i,j,0) = DS2d%mca_step(i,j,0) + IST%part_size(i,j,k) * & - ((IST%mH_snow(i,j,k) + IST%mH_pond(i,j,k))) + ((IST%mH_snow(i,j,k) + IST%mH_pond(i,j,k)) + IST%mH_pond_ice(i,j,k)) DS2d%mi_sum(i,j) = DS2d%mi_sum(i,j) + (IST%mH_ice(i,j,k)) * IST%part_size(i,j,k) DS2d%ice_cover(i,j) = DS2d%ice_cover(i,j) + IST%part_size(i,j,k) enddo ; enddo ; enddo diff --git a/src/SIS_slow_thermo.F90 b/src/SIS_slow_thermo.F90 index a2d224e7..dd5f4053 100644 --- a/src/SIS_slow_thermo.F90 +++ b/src/SIS_slow_thermo.F90 @@ -466,7 +466,7 @@ subroutine slow_thermodynamics(IST, dt_slow, CS, OSS, FIA, XSF, IOF, G, US, IG, if (CS%column_check) & call write_ice_statistics(IST, CS%Time, CS%n_calls, G, US, IG, CS%sum_output_CSp, & message=" Post_thermo A", check_column=.true.) - call adjust_ice_categories(IST%mH_ice, IST%mH_snow, IST%mH_pond, IST%part_size, & + call adjust_ice_categories(IST%mH_ice, IST%mH_snow, IST%mH_pond, IST%mH_pond_ice, IST%part_size, & IST%TrReg, G, IG, CS%SIS_transport_CSp) !Niki: add ridging? if (CS%column_check) & @@ -930,8 +930,8 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, US, IG) do m=1,NkIce ; m_lay(m) = IST%mH_ice(i,j,k) * I_Nk ; enddo ! mw/new - melt pond size is now adjusted here (rain ignored in resize, for now) - call ice_resize_SIS2(1-IST%part_size(i,j,0), IST%mH_pond(i,j,k), m_lay, & - enthalpy, S_col, Salin, FIA%fprec_top(i,j,k)*dt_slow, & + call ice_resize_SIS2(1-IST%part_size(i,j,0), IST%mH_pond(i,j,k), IST%mH_pond_ice(i,j,k), & + m_lay, enthalpy, S_col, Salin, FIA%fprec_top(i,j,k)*dt_slow, & FIA%lprec_top(i,j,k)*dt_slow, FIA%evap_top(i,j,k)*dt_slow, & FIA%tmelt(i,j,k), FIA%bmelt(i,j,k), NkIce, npassive, TrLay, & heat_to_ocn, h2o_ice_to_ocn, h2o_ocn_to_ice, evap_from_ocn, & @@ -995,7 +995,8 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, US, IG) + h2o_ocn_to_ice - h2o_ice_to_ocn & - (dt_slow*FIA%evap_top(i,j,k) - evap_from_ocn) - mass_here = IST%mH_snow(i,j,k) + IST%mH_pond(i,j,k) + IST%mH_ice(i,j,k) + mass_here = (IST%mH_snow(i,j,k) + IST%mH_pond(i,j,k)) + & + (IST%mH_pond_ice(i,j,k) + IST%mH_ice(i,j,k)) enth_here = IST%mH_snow(i,j,k) * IST%enth_snow(i,j,k,1) do m=1,NkIce enth_here = enth_here + (IST%mH_ice(i,j,k)*I_Nk) * IST%enth_ice(i,j,k,m) diff --git a/src/SIS_transport.F90 b/src/SIS_transport.F90 index d928b349..76c9e609 100644 --- a/src/SIS_transport.F90 +++ b/src/SIS_transport.F90 @@ -83,10 +83,12 @@ module SIS_transport type, public :: cell_average_state_type ; private real, allocatable, dimension(:,:,:) :: m_ice !< The mass of ice in each thickness category !! per unit total area in a cell [R Z ~> kg m-2]. - real, allocatable, dimension(:,:,:) :: m_snow !< The mass of ice in each thickness category + real, allocatable, dimension(:,:,:) :: m_snow !< The mass of snow in each thickness category !! per unit total area in a cell [R Z ~> kg m-2]. real, allocatable, dimension(:,:,:) :: m_pond !< The mass of melt pond water in each thickness !! category per unit total area in a cell [R Z ~> kg m-2]. + real, allocatable, dimension(:,:,:) :: m_pond_ice !< The mass of melt pond ice in each thickness + !! category per unit total area in a cell [R Z ~> kg m-2]. real, allocatable, dimension(:,:,:) :: mH_ice !< The mass of ice in each thickness category !! per unit of ice area in a cell [R Z ~> kg m-2]. The !! ratio of m_ice / mH_ice gives the fractional @@ -136,19 +138,23 @@ subroutine ice_cat_transport(CAS, TrReg, dt_slow, nsteps, G, US, IG, CS, uc, vc, ! Local variables real, dimension(SZIB_(G),SZJ_(G),SZCAT_(IG)) :: & - uh_ice, & ! Zonal fluxes of ice [R Z L2 T-1 ~> kg s-1]. - uh_snow, & ! Zonal fluxes of snow [R Z L2 T-1 ~> kg s-1]. - uh_pond ! Zonal fluxes of melt pond water [R Z L2 T-1 ~> kg s-1]. + uh_ice, & ! Zonal fluxes of ice [R Z L2 T-1 ~> kg s-1]. + uh_snow, & ! Zonal fluxes of snow [R Z L2 T-1 ~> kg s-1]. + uh_pond, & ! Zonal fluxes of melt pond water [R Z L2 T-1 ~> kg s-1]. + uh_pond_ice ! Zonal fluxes of melt pond ice [R Z L2 T-1 ~> kg s-1]. real, dimension(SZI_(G),SZJB_(G),SZCAT_(IG)) :: & - vh_ice, & ! Meridional fluxes of ice [R Z L2 T-1 ~> kg s-1]. - vh_snow, & ! Meridional fluxes of snow [R Z L2 T-1 ~> kg s-1]. - vh_pond ! Meridional fluxes of melt pond water [R Z L2 T-1 ~> kg s-1]. + vh_ice, & ! Meridional fluxes of ice [R Z L2 T-1 ~> kg s-1]. + vh_snow, & ! Meridional fluxes of snow [R Z L2 T-1 ~> kg s-1]. + vh_pond, & ! Meridional fluxes of melt pond water [R Z L2 T-1 ~> kg s-1]. + vh_pond_ice ! Meridional fluxes of melt pond ice [R Z L2 T-1 ~> kg s-1]. real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)) :: & - mca0_ice, & ! The initial mass of ice per unit ocean area in a cell [R Z ~> kg m-2]. - mca0_snow, & ! The initial mass of snow per unit ocean area in a cell [R Z ~> kg m-2]. - mca0_pond ! The initial mass of melt pond water per unit ocean area - ! in a cell [R Z ~> kg m-2]. - real :: dt_adv ! An advective timestep [T ~> s] + mca0_ice, & ! The initial mass of ice per unit ocean area in a cell [R Z ~> kg m-2]. + mca0_snow, & ! The initial mass of snow per unit ocean area in a cell [R Z ~> kg m-2]. + mca0_pond, & ! The initial mass of melt pond water per unit ocean area + ! in a cell [R Z ~> kg m-2]. + mca0_pond_ice ! The initial mass of melt pond ice per unit ocean area + ! in a cell [R Z ~> kg m-2]. + real :: dt_adv ! An advective timestep [T ~> s] logical :: merged_cont character(len=200) :: mesg integer :: i, j, k, n, isc, iec, jsc, jec, isd, ied, jsd, jed, nCat @@ -180,12 +186,14 @@ subroutine ice_cat_transport(CAS, TrReg, dt_slow, nsteps, G, US, IG, CS, uc, vc, call pass_var(CAS%m_ice, G%Domain, complete=.false.) call pass_var(CAS%m_snow, G%Domain, complete=.false.) call pass_var(CAS%m_pond, G%Domain, complete=.false.) + call pass_var(CAS%m_pond_ice, G%Domain, complete=.false.) call pass_var(CAS%mH_ice, G%Domain, complete=.true.) do k=1,nCat ; do j=jsd,jed ; do i=isd,ied mca0_ice(i,j,k) = CAS%m_ice(i,j,k) mca0_snow(i,j,k) = CAS%m_snow(i,j,k) mca0_pond(i,j,k) = CAS%m_pond(i,j,k) + mca0_pond_ice(i,j,k) = CAS%m_pond_ice(i,j,k) enddo ; enddo ; enddo if (merged_cont) then @@ -193,7 +201,8 @@ subroutine ice_cat_transport(CAS, TrReg, dt_slow, nsteps, G, US, IG, CS, uc, vc, dt_adv, G, US, IG, CS%continuity_CSp, & h1=CAS%m_ice, uh1=uh_ice, vh1=vh_ice, & h2=CAS%m_snow, uh2=uh_snow, vh2=vh_snow, & - h3=CAS%m_pond, uh3=uh_pond, vh3=vh_pond) + h3=CAS%m_pond, uh3=uh_pond, vh3=vh_pond, & + h4=CAS%m_pond_ice, uh4=uh_pond_ice, vh4=vh_pond_ice) else call continuity(uc, vc, mca0_ice, CAS%m_ice, uh_ice, vh_ice, dt_adv, & G, US, IG, CS%continuity_CSp, use_h_neg=.true.) ! this is hard-coded here to preserve previous answers @@ -201,6 +210,8 @@ subroutine ice_cat_transport(CAS, TrReg, dt_slow, nsteps, G, US, IG, CS, uc, vc, G, US, IG, CS%continuity_CSp, masking_uh=uh_ice, masking_vh=vh_ice, use_h_neg=h_neg_fix) call continuity(uc, vc, mca0_pond, CAS%m_pond, uh_pond, vh_pond, dt_adv, & G, US, IG, CS%continuity_CSp, masking_uh=uh_ice, masking_vh=vh_ice, use_h_neg=h_neg_fix) + call continuity(uc, vc, mca0_pond_ice, CAS%m_pond, uh_pond_ice, vh_pond_ice, dt_adv, & + G, US, IG, CS%continuity_CSp, masking_uh=uh_ice, masking_vh=vh_ice, use_h_neg=h_neg_fix) endif call advect_scalar(CAS%mH_ice, mca0_ice, CAS%m_ice, uh_ice, vh_ice, & @@ -213,10 +224,12 @@ subroutine ice_cat_transport(CAS, TrReg, dt_slow, nsteps, G, US, IG, CS, uc, vc, ! Accumulated diagnostics CAS%dt_sum = CAS%dt_sum + dt_adv if (allocated(CAS%uh_sum)) then ; do k=1,nCat ; do j=jsc,jec ; do I=isc-1,iec - CAS%uh_sum(I,j) = CAS%uh_sum(I,j) + dt_adv * ((uh_pond(I,j,k) + uh_snow(I,j,k)) + uh_ice(I,j,k)) + CAS%uh_sum(I,j) = CAS%uh_sum(I,j) + dt_adv * ((uh_pond(I,j,k) + uh_pond_ice(I,j,k)) + (uh_ice(I,j,k) + & + uh_snow(I,j,k))) enddo ; enddo ; enddo ; endif if (allocated(CAS%vh_sum)) then ; do k=1,nCat ; do J=jsc-1,jec ; do i=isc,iec - CAS%vh_sum(i,J) = CAS%vh_sum(i,J) + dt_adv * ((vh_pond(i,J,k) + vh_snow(i,J,k)) + vh_ice(i,J,k)) + CAS%vh_sum(i,J) = CAS%vh_sum(i,J) + dt_adv * ((vh_pond(i,J,k) + vh_pond_ice(i,J,k)) + (vh_ice(i,J,k) + & + vh_snow(i,J,k))) enddo ; enddo ; enddo ; endif if (CS%bounds_check) then @@ -287,7 +300,7 @@ subroutine finish_ice_transport(CAS, IST, TrReg, G, US, IG, dt, CS, OSS, rdg_rat if (CS%bounds_check) call check_SIS_tracer_bounds(TrReg, G, IG, "After compress_ice") if (CS%readjust_categories) then - call adjust_ice_categories(IST%mH_ice, IST%mH_snow, IST%mH_pond, IST%part_size, & + call adjust_ice_categories(IST%mH_ice, IST%mH_snow, IST%mH_pond, IST%mH_pond_ice, IST%part_size, & TrReg, G, IG, CS) if (CS%bounds_check) call check_SIS_tracer_bounds(TrReg, G, IG, "After adjust_ice_categories") endif @@ -344,6 +357,7 @@ subroutine finish_ice_transport(CAS, IST, TrReg, G, US, IG, dt, CS, OSS, rdg_rat call pass_var(IST%part_size, G%Domain) ! cannot be combined with the three updates below call pass_var(IST%mH_pond, G%Domain, complete=.false.) + call pass_var(IST%mH_pond_ice, G%Domain, complete=.false.) call pass_var(IST%mH_snow, G%Domain, complete=.false.) call pass_var(IST%mH_ice, G%Domain, complete=.true.) @@ -431,6 +445,7 @@ subroutine ice_state_to_cell_ave_state(IST, G, US, IG, CS, CAS) isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nCat = IG%CatIce CAS%m_ice(:,:,:) = 0.0 ; CAS%m_snow(:,:,:) = 0.0 ; CAS%m_pond(:,:,:) = 0.0 ; CAS%mH_ice(:,:,:) = 0.0 + CAS%m_pond_ice(:,:,:) = 0.0 ice_cover(:,:) = 0.0 ; mHi_avg(:,:) = 0.0 !$OMP parallel do default(shared) do j=jsc,jec @@ -439,6 +454,7 @@ subroutine ice_state_to_cell_ave_state(IST, G, US, IG, CS, CAS) CAS%m_ice(i,j,k) = IST%part_size(i,j,k) * IST%mH_ice(i,j,k) CAS%m_snow(i,j,k) = IST%part_size(i,j,k) * IST%mH_snow(i,j,k) CAS%m_pond(i,j,k) = IST%part_size(i,j,k) * IST%mH_pond(i,j,k) + CAS%m_pond_ice(i,j,k) = IST%part_size(i,j,k) * IST%mH_pond_ice(i,j,k) CAS%mH_ice(i,j,k) = IST%mH_ice(i,j,k) ice_cover(i,j) = ice_cover(i,j) + IST%part_size(i,j,k) mHi_avg(i,j) = mHi_avg(i,j) + CAS%m_ice(i,j,k) @@ -450,6 +466,7 @@ subroutine ice_state_to_cell_ave_state(IST, G, US, IG, CS, CAS) call SIS_error(FATAL, "Input to SIS_transport, non-zero pond mass rests atop no ice.") endif CAS%m_ice(i,j,k) = 0.0 ; CAS%m_snow(i,j,k) = 0.0 ; CAS%m_pond(i,j,k) = 0.0 + CAS%m_pond_ice(i,j,k) = 0.0 endif enddo ; enddo do i=isc,iec ; if (ice_cover(i,j) > 0.0) then @@ -533,6 +550,7 @@ subroutine cell_ave_state_to_ice_state(CAS, G, US, IG, CS, IST, TrReg) IST%part_size(i,j,k) = CAS%m_ice(i,j,k) / CAS%mH_ice(i,j,k) IST%mH_snow(i,j,k) = CAS%mH_ice(i,j,k) * (CAS%m_snow(i,j,k) / CAS%m_ice(i,j,k)) IST%mH_pond(i,j,k) = CAS%mH_ice(i,j,k) * (CAS%m_pond(i,j,k) / CAS%m_ice(i,j,k)) + IST%mH_pond_ice(i,j,k) = CAS%mH_ice(i,j,k) * (CAS%m_pond_ice(i,j,k) / CAS%m_ice(i,j,k)) IST%mH_ice(i,j,k) = CAS%mH_ice(i,j,k) ice_cover(i,j) = ice_cover(i,j) + IST%part_size(i,j,k) else @@ -543,7 +561,7 @@ subroutine cell_ave_state_to_ice_state(CAS, G, US, IG, CS, IST, TrReg) if (CAS%m_pond(i,j,k) > mass_neglect ) & call SIS_error(FATAL, & "Something needs to be done with positive CAS%m_pond values without ice.") - IST%mH_snow(i,j,k) = 0.0 ; IST%mH_pond(i,j,k) = 0.0 + IST%mH_snow(i,j,k) = 0.0 ; IST%mH_pond(i,j,k) = 0.0 ; IST%mH_pond_ice(i,j,k) = 0.0 endif enddo ; enddo ; enddo do j=jsc,jec ; do i=isc,iec @@ -554,7 +572,7 @@ end subroutine cell_ave_state_to_ice_state !> adjust_ice_categories moves mass between thickness categories if it is thinner or !! thicker than the bounding limits of each category. -subroutine adjust_ice_categories(mH_ice, mH_snow, mH_pond, part_sz, TrReg, G, IG, CS) +subroutine adjust_ice_categories(mH_ice, mH_snow, mH_pond, mH_pond_ice, part_sz, TrReg, G, IG, CS) type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type type(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), & @@ -566,6 +584,9 @@ subroutine adjust_ice_categories(mH_ice, mH_snow, mH_pond, part_sz, TrReg, G, IG real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), & intent(inout) :: mH_pond !< The mass per unit area of the pond !! on the ice in each category [R Z ~> kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), & + intent(inout) :: mH_pond_ice !< The mass per unit area of the pond ice + !! on the pond in each category [R Z ~> kg m-2]. real, dimension(SZI_(G),SZJ_(G),0:SZCAT_(IG)), & intent(inout) :: part_sz !< The fractional ice concentration !! within a cell in each thickness @@ -582,15 +603,16 @@ subroutine adjust_ice_categories(mH_ice, mH_snow, mH_pond, part_sz, TrReg, G, IG real :: part_trans ! The fractional area transferred between categories [nondim]. real :: snow_trans ! The cell-averaged snow transferred between categories [R Z ~> kg m-2]. real :: pond_trans ! The cell-averaged pond transferred between categories [R Z ~> kg m-2]. + real :: pond_ice_trans ! The cell-averaged pond ice transferred between categories [R Z ~> kg m-2]. real :: I_mH_lim1 ! The inverse of the lower thickness limit [R-1 Z-1 ~> m2 kg-1]. real, dimension(SZI_(G),SZCAT_(IG)) :: & ! The mass of snow, pond and ice per unit total area in a cell [R Z ~> kg m-2]. ! "mca" stands for "mass cell averaged" - mca_ice, mca_snow, mca_pond, & + mca_ice, mca_snow, mca_pond, mca_pond_ice, & ! Initial ice, snow and pond masses per unit cell area [R Z ~> kg m-2]. - mca0_ice, mca0_snow, mca0_pond, & + mca0_ice, mca0_snow, mca0_pond, mca0_pond_ice, & ! Cross-category transfers of ice, snow and pond mass [R Z ~> kg m-2]. - trans_ice, trans_snow, trans_pond + trans_ice, trans_snow, trans_pond, trans_pond_ice real, dimension(SZI_(G)) :: ice_cover ! The summed fractional ice coverage [nondim]. logical :: do_any, do_j(SZJ_(G)), resum_cat(SZI_(G), SZJ_(G)) integer :: i, j, k, m, is, ie, js, je, nCat @@ -619,6 +641,9 @@ subroutine adjust_ice_categories(mH_ice, mH_snow, mH_pond, part_sz, TrReg, G, IG call SIS_error(WARNING, mesg, all_print=.true.) call SIS_error(FATAL, "Input to adjust_ice_categories, non-zero pond mass rests atop no ice.") endif + if (mH_pond_ice(i,j,k) > 0.0) then + call SIS_error(FATAL, "Input to adjust_ice_categories, non-zero pond ice mass rests atop no ice.") + endif if (part_sz(i,j,k) > 0.0) resum_cat(i,j) = .true. part_sz(i,j,k) = 0.0 endif ; enddo ; enddo ; enddo @@ -639,12 +664,15 @@ subroutine adjust_ice_categories(mH_ice, mH_snow, mH_pond, part_sz, TrReg, G, IG mca_ice(i,k) = part_sz(i,j,k)*mH_ice(i,j,k) mca_snow(i,k) = part_sz(i,j,k)*mH_snow(i,j,k) mca_pond(i,k) = part_sz(i,j,k)*mH_pond(i,j,k) + mca_pond_ice(i,k) = part_sz(i,j,k)*mH_pond_ice(i,j,k) mca0_ice(i,k) = mca_ice(i,k) mca0_snow(i,k) = mca_snow(i,k) mca0_pond(i,k) = mca_pond(i,k) + mca0_pond_ice(i,k) = mca_pond_ice(i,k) enddo ; enddo trans_ice(:,:) = 0.0 ; trans_snow(:,:) = 0.0 ; trans_pond(:,:) = 0.0 + trans_pond_ice(:,:) = 0.0 do_any = .false. do k=1,nCat-1 ; do i=is,ie @@ -655,10 +683,12 @@ subroutine adjust_ice_categories(mH_ice, mH_snow, mH_pond, part_sz, TrReg, G, IG part_trans = part_sz(i,j,k) ! * (mca_trans / mca_ice) * (mH_ice / h_trans) snow_trans = mca_snow(i,k) ! * (part_trans / part_sz) = 1 pond_trans = mca_pond(i,k) ! * (part_trans / part_sz) = 1 + pond_ice_trans = mca_pond_ice(i,k) ! * (part_trans / part_sz) = 1 trans_ice(i,K) = mca_trans trans_snow(i,K) = snow_trans trans_pond(i,K) = pond_trans + trans_pond_ice(i,K) = pond_ice_trans do_any = .true. ! Use area-weighted remapped thicknesses so that the total ice area and @@ -690,6 +720,13 @@ subroutine adjust_ice_categories(mH_ice, mH_snow, mH_pond, part_sz, TrReg, G, IG mH_pond(i,j,k) = 0.0 ; mH_pond(i,j,k+1) = 0.0 if (part_sz(i,j,k)>0.0) mH_pond(i,j,k) = mca_pond(i,k) / part_sz(i,j,k) if (part_sz(i,j,k+1)>0.0) mH_pond(i,j,k+1) = mca_pond(i,k+1) / part_sz(i,j,k+1) + + mca_pond_ice(i,k+1) = mca_pond_ice(i,k+1) + pond_ice_trans + mca_pond_ice(i,k) = mca_pond_ice(i,k) - pond_ice_trans + + mH_pond_ice(i,j,k) = 0.0 ; mH_pond_ice(i,j,k+1) = 0.0 + if (part_sz(i,j,k)>0.0) mH_pond_ice(i,j,k) = mca_pond_ice(i,k) / part_sz(i,j,k) + if (part_sz(i,j,k+1)>0.0) mH_pond_ice(i,j,k+1) = mca_pond_ice(i,k+1) / part_sz(i,j,k+1) endif enddo ; enddo @@ -704,8 +741,10 @@ subroutine adjust_ice_categories(mH_ice, mH_snow, mH_pond, part_sz, TrReg, G, IG mca0_ice(i,k) = mca_ice(i,k) mca0_snow(i,k) = mca_snow(i,k) mca0_pond(i,k) = mca_pond(i,k) + mca0_pond_ice(i,k) = mca_pond_ice(i,k) enddo ; enddo trans_ice(:,:) = 0.0 ; trans_snow(:,:) = 0.0 ; trans_pond(:,:) = 0.0 + trans_pond_ice(:,:) = 0.0 do_any = .false. do k=nCat,2,-1 ; do i=is,ie @@ -716,11 +755,13 @@ subroutine adjust_ice_categories(mH_ice, mH_snow, mH_pond, part_sz, TrReg, G, IG part_trans = part_sz(i,j,k) ! * (mca_trans / mca_ice) * (mH_ice / h_trans) snow_trans = mca_snow(i,k) ! * (part_trans / part_sz) = 1 pond_trans = mca_pond(i,k) ! * (part_trans / part_sz) = 1 + pond_ice_trans = mca_pond_ice(i,k) ! * (part_trans / part_sz) = 1 do_any = .true. trans_ice(i,K-1) = -mca_trans ! Note the shifted index conventions! trans_snow(i,K-1) = -snow_trans trans_pond(i,K-1) = -pond_trans + trans_pond_ice(i,K-1) = -pond_ice_trans ! Use area-weighted remapped thicknesses so that the total ice area and ! mass are both conserved in the remapping operation. @@ -751,6 +792,13 @@ subroutine adjust_ice_categories(mH_ice, mH_snow, mH_pond, part_sz, TrReg, G, IG mH_pond(i,j,k) = 0.0 ; mH_pond(i,j,k-1) = 0.0 if (part_sz(i,j,k)>0.0) mH_pond(i,j,k) = mca_pond(i,k) / part_sz(i,j,k) if (part_sz(i,j,k-1)>0.0) mH_pond(i,j,k-1) = mca_pond(i,k-1) / part_sz(i,j,k-1) + + mca_pond_ice(i,k-1) = mca_pond_ice(i,k-1) + pond_ice_trans + mca_pond_ice(i,k) = mca_pond_ice(i,k) - pond_ice_trans + + mH_pond_ice(i,j,k) = 0.0 ; mH_pond_ice(i,j,k-1) = 0.0 + if (part_sz(i,j,k)>0.0) mH_pond_ice(i,j,k) = mca_pond_ice(i,k) / part_sz(i,j,k) + if (part_sz(i,j,k-1)>0.0) mH_pond_ice(i,j,k-1) = mca_pond_ice(i,k-1) / part_sz(i,j,k-1) endif enddo ; enddo @@ -771,6 +819,7 @@ subroutine adjust_ice_categories(mH_ice, mH_snow, mH_pond, part_sz, TrReg, G, IG part_sz(i,j,1) = part_sz(i,j,1) * (mH_ice(i,j,1) * I_mH_lim1) mH_snow(i,j,1) = mH_snow(i,j,1) * (IG%mH_cat_bound(1) / mH_ice(i,j,1)) mH_pond(i,j,1) = mH_pond(i,j,1) * (IG%mH_cat_bound(1) / mH_ice(i,j,1)) + mH_pond_ice(i,j,1) = mH_pond_ice(i,j,1) * (IG%mH_cat_bound(1) / mH_ice(i,j,1)) ! This is equivalent to mH_snow(i,j,1) = mca_snow(i,1) / part_sz(i,j,1) mH_ice(i,j,1) = IG%mH_cat_bound(1) resum_cat(i,j) = .true. @@ -1015,17 +1064,18 @@ subroutine compress_ice(part_sz, mH_ice, mH_snow, mH_pond, TrReg, G, US, IG, CS, end subroutine compress_ice !> get_total_mass determines the globally integrated mass of snow and ice -subroutine get_total_mass(IST, G, US, IG, tot_ice, tot_snow, tot_pond, scale) +subroutine get_total_mass(IST, G, US, IG, tot_ice, tot_snow, tot_pond, tot_pond_ice, scale) type(ice_state_type), intent(in) :: IST !< A type describing the state of the sea ice type(SIS_hor_grid_type), intent(in) :: G !< The horizontal grid type type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type type(EFP_type), intent(out) :: tot_ice !< The globally integrated total ice [kg]. type(EFP_type), intent(out) :: tot_snow !< The globally integrated total snow [kg]. - type(EFP_type),optional, intent(out) :: tot_pond !< The globally integrated total snow [kg]. + type(EFP_type),optional, intent(out) :: tot_pond !< The globally integrated total pond water [kg]. + type(EFP_type),optional, intent(out) :: tot_pond_ice !< The globally integrated total pond ice [kg]. real, optional, intent(in) :: scale !< A scaling factor from H to the desired units. - real, dimension(G%isc:G%iec, G%jsc:G%jec) :: sum_ice, sum_snow, sum_pond + real, dimension(G%isc:G%iec, G%jsc:G%jec) :: sum_ice, sum_snow, sum_pond, sum_pond_ice real :: H_to_units ! A conversion factor from H to the desired output units. real :: total integer :: i, j, k, m, isc, iec, jsc, jec @@ -1043,11 +1093,15 @@ subroutine get_total_mass(IST, G, US, IG, tot_ice, tot_snow, tot_pond, scale) if (present(tot_pond)) & sum_pond(i,j) = sum_pond(i,j) + G%areaT(i,j) * & (IST%part_size(i,j,k) * (H_to_units*IST%mH_pond(i,j,k))) + if (present(tot_pond_ice)) & + sum_pond_ice(i,j) = sum_pond_ice(i,j) + G%areaT(i,j) * & + (IST%part_size(i,j,k) * (H_to_units*IST%mH_pond_ice(i,j,k))) enddo ; enddo ; enddo total = reproducing_sum(sum_ice, EFP_sum=tot_ice) total = reproducing_sum(sum_snow, EFP_sum=tot_snow) if (present(tot_pond)) total = reproducing_sum(sum_pond, EFP_sum=tot_pond) + if (present(tot_pond_ice)) total = reproducing_sum(sum_pond_ice, EFP_sum=tot_pond_ice) end subroutine get_total_mass @@ -1068,7 +1122,8 @@ subroutine get_cell_mass(IST, G, IG, cell_mass, scale) cell_mass(:,:) = 0.0 do k=1,IG%CatIce ; do j=jsc,jec ; do i=isc,iec cell_mass(i,j) = cell_mass(i,j) + IST%part_size(i,j,k) * H_to_units * & - ((IST%mH_snow(i,j,k) + IST%mH_pond(i,j,k)) + IST%mH_ice(i,j,k)) + ((IST%mH_snow(i,j,k) + IST%mH_pond(i,j,k)) + & + (IST%mH_pond_ice(i,j,k) + IST%mH_ice(i,j,k))) enddo ; enddo ; enddo end subroutine get_cell_mass @@ -1090,7 +1145,8 @@ subroutine cell_mass_from_CAS(CAS, G, IG, mca, scale) do j=jsc,jec ; do i=isc,iec ; mca(i,j) = 0.0 ; enddo ; enddo do k=1,nCat ; do j=jsc,jec ; do i=isc,iec - mca(i,j) = mca(i,j) + H_to_units * (CAS%m_ice(i,j,k) + (CAS%m_snow(i,j,k) + CAS%m_pond(i,j,k))) + mca(i,j) = mca(i,j) + H_to_units * ((CAS%m_ice(i,j,k) + CAS%m_snow(i,j,k)) + \ + (CAS%m_pond(i,j,k) + CAS%m_pond_ice(i,j,k))) enddo ; enddo ; enddo end subroutine cell_mass_from_CAS @@ -1288,6 +1344,7 @@ subroutine alloc_cell_average_state_type(CAS, HI, IG, CS) call safe_alloc(CAS%m_ice, isd, ied, jsd, jed, ncat) call safe_alloc(CAS%m_snow, isd, ied, jsd, jed, ncat) call safe_alloc(CAS%m_pond, isd, ied, jsd, jed, ncat) + call safe_alloc(CAS%m_pond_ice, isd, ied, jsd, jed, ncat) call safe_alloc(CAS%mH_ice, isd, ied, jsd, jed, ncat) if (present(CS)) then @@ -1305,7 +1362,7 @@ subroutine dealloc_cell_average_state_type(CAS) type(cell_average_state_type), pointer :: CAS !< A structure with ocean-cell averaged masses !! that is being allocated here. if (.not.associated(CAS)) return - deallocate(CAS%m_ice, CAS%m_snow, CAS%m_pond, CAS%mH_ice) + deallocate(CAS%m_ice, CAS%m_snow, CAS%m_pond, CAS%m_pond_ice, CAS%mH_ice) if (allocated(CAS%mass0)) deallocate(CAS%mass0) if (allocated(CAS%uh_sum)) deallocate(CAS%uh_sum) if (allocated(CAS%vh_sum)) deallocate(CAS%vh_sum) diff --git a/src/SIS_types.F90 b/src/SIS_types.F90 index b1cb332c..60bfd31b 100644 --- a/src/SIS_types.F90 +++ b/src/SIS_types.F90 @@ -74,6 +74,7 @@ module SIS_types real, allocatable, dimension(:,:,:) :: & mH_pond, & !< The mass per unit area of the pond in each category [R Z ~> kg m-2]. + mH_pond_ice, & !< The mass per unit area of the pond ice in each category [R Z ~> kg m-2]. mH_snow, & !< The mass per unit area of the snow in each category [R Z ~> kg m-2]. mH_ice, & !< The mass per unit area of the ice in each category [R Z ~> kg m-2]. t_surf !< The surface temperature [Kelvin]. @@ -449,12 +450,13 @@ subroutine alloc_IST_arrays(HI, IG, US, IST, omit_velocities, omit_Tsurf, do_rid IST%valid_IST = .true. allocate(IST%part_size(isd:ied, jsd:jed, 0:CatIce), source=0.0) - allocate(IST%mH_pond( isd:ied, jsd:jed, CatIce), source=0.0) - allocate(IST%mH_snow( isd:ied, jsd:jed, CatIce), source=0.0) - allocate(IST%enth_snow(isd:ied, jsd:jed, CatIce, 1), source=0.0) - allocate(IST%mH_ice( isd:ied, jsd:jed, CatIce), source=0.0) - allocate(IST%enth_ice( isd:ied, jsd:jed, CatIce, NkIce), source=0.0) - allocate(IST%sal_ice( isd:ied, jsd:jed, CatIce, NkIce), source=0.0) + allocate(IST%mH_pond( isd:ied, jsd:jed, CatIce), source=0.0) + allocate(IST%mH_pond_ice(isd:ied, jsd:jed, CatIce), source=0.0) + allocate(IST%mH_snow( isd:ied, jsd:jed, CatIce), source=0.0) + allocate(IST%enth_snow( isd:ied, jsd:jed, CatIce, 1), source=0.0) + allocate(IST%mH_ice( isd:ied, jsd:jed, CatIce), source=0.0) + allocate(IST%enth_ice( isd:ied, jsd:jed, CatIce, NkIce), source=0.0) + allocate(IST%sal_ice( isd:ied, jsd:jed, CatIce, NkIce), source=0.0) if (present(do_ridging)) then ; if (do_ridging) then allocate(IST%snow_to_ocn(isd:ied, jsd:jed), source=0.0) @@ -503,6 +505,8 @@ subroutine ice_state_register_restarts(IST, G, IG, US, Ice_restart) endif call register_restart_field(Ice_restart, 'h_pond', IST%mH_pond, & mandatory=.false., units="kg m-2", conversion=US%RZ_to_kg_m2) + call register_restart_field(Ice_restart, 'h_pond_ice', IST%mH_pond_ice, & + mandatory=.false., units="kg m-2", conversion=US%RZ_to_kg_m2) call register_restart_field(Ice_restart, 'h_snow', IST%mH_snow, & mandatory=.true., units="kg m-2", conversion=US%RZ_to_kg_m2) call register_restart_field(Ice_restart, 'enth_snow', IST%enth_snow, dim_4='z_snow', & @@ -1072,6 +1076,7 @@ subroutine copy_IST_to_IST(IST_in, IST_out, HI_in, HI_out, IG) do k=1,ncat ; do j=jsc,jec ; do i=isc,iec ; i2 = i+i_off ; j2 = j+j_off IST_out%mH_pond(i2,j2,k) = IST_in%mH_pond(i,j,k) + IST_out%mH_pond_ice(i2,j2,k) = IST_in%mH_pond_ice(i,j,k) IST_out%mH_snow(i2,j2,k) = IST_in%mH_snow(i,j,k) IST_out%mH_ice(i2,j2,k) = IST_in%mH_ice(i,j,k) @@ -1111,6 +1116,8 @@ subroutine redistribute_IST_to_IST(IST_in, IST_out, domain_in, domain_out) endif call redistribute_data(domain_in, IST_in%mH_pond, domain_out, & IST_out%mH_pond, complete=.false.) + call redistribute_data(domain_in, IST_in%mH_pond_ice, domain_out, & + IST_out%mH_pond_ice, complete=.false.) call redistribute_data(domain_in, IST_in%mH_snow, domain_out, & IST_out%mH_snow, complete=.false.) call redistribute_data(domain_in, IST_in%mH_ice, domain_out, & @@ -1133,6 +1140,8 @@ subroutine redistribute_IST_to_IST(IST_in, IST_out, domain_in, domain_out) endif call redistribute_data(domain_in, null_ptr3D, domain_out, & IST_out%mH_pond, complete=.false.) + call redistribute_data(domain_in, null_ptr3D, domain_out, & + IST_out%mH_pond_ice, complete=.false.) call redistribute_data(domain_in, null_ptr3D, domain_out, & IST_out%mH_snow, complete=.false.) call redistribute_data(domain_in, null_ptr3D, domain_out, & @@ -1155,6 +1164,8 @@ subroutine redistribute_IST_to_IST(IST_in, IST_out, domain_in, domain_out) endif call redistribute_data(domain_in, IST_in%mH_pond, domain_out, & null_ptr3D, complete=.false.) + call redistribute_data(domain_in, IST_in%mH_pond_ice, domain_out, & + null_ptr3D, complete=.false.) call redistribute_data(domain_in, IST_in%mH_snow, domain_out, & null_ptr3D, complete=.false.) call redistribute_data(domain_in, IST_in%mH_ice, domain_out, & @@ -2041,7 +2052,7 @@ subroutine dealloc_IST_arrays(IST) type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice deallocate(IST%part_size, IST%mH_snow, IST%mH_ice) - deallocate(IST%mH_pond) ! mw/new + deallocate(IST%mH_pond, IST%mH_pond_ice) ! mw/new deallocate(IST%enth_snow, IST%enth_ice, IST%sal_ice) if (allocated(IST%snow_to_ocn)) deallocate(IST%snow_to_ocn) if (allocated(IST%enth_snow_to_ocn)) deallocate(IST%enth_snow_to_ocn) @@ -2405,6 +2416,7 @@ subroutine IST_chksum(mesg, IST, G, US, IG, haloshift) call hchksum(IST%mH_snow, trim(mesg)//" IST%mH_snow", G%HI, haloshift=hs, unscale=US%RZ_to_kg_m2) call hchksum(IST%enth_snow(:,:,:,1), trim(mesg)//" IST%enth_snow", G%HI, haloshift=hs, unscale=US%Q_to_J_kg) call hchksum(IST%mH_pond, trim(mesg)//" IST%mH_pond", G%HI, haloshift=hs, unscale=US%RZ_to_kg_m2) + call hchksum(IST%mH_pond_ice, trim(mesg)//" IST%mH_pond_ice", G%HI, haloshift=hs, unscale=US%RZ_to_kg_m2) if (allocated(IST%u_ice_B) .and. allocated(IST%v_ice_B)) then call Bchksum_pair(mesg//" IST%[uv]_ice_B", IST%u_ice_B, IST%v_ice_B, G, halos=hs, unscale=US%L_T_to_m_s) diff --git a/src/ice_model.F90 b/src/ice_model.F90 index 2eb1f9bf..68d38f3b 100644 --- a/src/ice_model.F90 +++ b/src/ice_model.F90 @@ -306,6 +306,7 @@ subroutine update_ice_dynamics_trans(Ice, time_step, start_cycle, end_cycle, cyc call pass_var(sIST%part_size, sG%Domain) call pass_var(sIST%mH_ice, sG%Domain, complete=.false.) call pass_var(sIST%mH_pond, sG%Domain, complete=.false.) + call pass_var(sIST%mH_pond_ice, sG%Domain, complete=.false.) call pass_var(sIST%mH_snow, sG%Domain, complete=.true.) endif @@ -618,7 +619,8 @@ subroutine set_ocean_top_fluxes(Ice, IST, IOF, FIA, OSS, G, US, IG, sCS) ! do j=jsc,jec ; do k=1,ncat ; do i=isc,iec ! i2 = i+i_off ; j2 = j+j_off! Use these to correct for indexing differences. ! Ice%mi(i2,j2) = Ice%mi(i2,j2) + IST%part_size(i,j,k) * & -! (G%US%RZ_to_kg_m2 * ((IST%mH_snow(i,j,k) + IST%mH_pond(i,j,k)) + IST%mH_ice(i,j,k))) +! (G%US%RZ_to_kg_m2 * ((IST%mH_snow(i,j,k) + IST%mH_pond(i,j,k)) + & +! (IST%mH_pond_ice(i,j,k) + IST%mH_ice(i,j,k)))) ! enddo ; enddo ; enddo ! This block of code is probably unnecessary. @@ -701,7 +703,7 @@ subroutine ice_mass_from_IST(IST, IOF, G, IG) !$OMP parallel do default(shared) do j=jsc,jec ; do k=1,ncat ; do i=isc,iec IOF%mass_ice_sn_p(i,j) = IOF%mass_ice_sn_p(i,j) + IST%part_size(i,j,k) * & - ((IST%mH_snow(i,j,k) + IST%mH_pond(i,j,k)) + IST%mH_ice(i,j,k)) + ((IST%mH_snow(i,j,k) + IST%mH_pond(i,j,k)) + (IST%mH_pond_ice(i,j,k) + IST%mH_ice(i,j,k))) enddo ; enddo ; enddo end subroutine ice_mass_from_IST @@ -1769,10 +1771,11 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, logical :: redo_fast_update ! If true, recalculate the thermal updates from the fast ! dynamics on the slowly evolving ice state, rather than ! copying over the slow ice state to the fast ice state. - logical :: do_mask_restart ! If true, apply the scaling and masks to mH_snow, mH_ice, part_size - ! mH_pond, t_surf, t_skin, sal_ice, enth_ice and enth_snow - ! after a restart. However this may cause answers to diverge - ! after a restart.Provide a switch to turn this option off. + logical :: do_mask_restart ! If true, apply the scaling and masks to mH_snow, mH_ice, + ! part_size, mH_pond, mH_pond_ice, t_surf, t_skin, sal_ice, + ! enth_ice and enth_snow after a restart. However this may + ! cause answers to diverge after a restart. Provide a switch + ! to turn this option off. logical :: recategorize_ice ! If true, adjust the distribution of the ice among thickness ! categories after initialization. logical :: do_brine_plume ! If true, keep track of how much salt is left in the ocean @@ -2487,6 +2490,7 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, sIST%enth_snow(i,j,k,1) = sIST%enth_snow(i,j,k,1) * sG%mask2dT(i,j) sIST%mH_ice(i,j,k) = sIST%mH_ice(i,j,k) * sG%mask2dT(i,j) sIST%mH_pond(i,j,k) = sIST%mH_pond(i,j,k) * sG%mask2dT(i,j) + sIST%mH_pond_ice(i,j,k) = sIST%mH_pond_ice(i,j,k) * sG%mask2dT(i,j) sIST%part_size(i,j,k) = sIST%part_size(i,j,k) * sG%mask2dT(i,j) enddo ; enddo ; enddo ! Since we masked out the part_size on land we should set @@ -2499,7 +2503,7 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, if (recategorize_ice) then ! Transfer ice to the correct thickness categories. This does not change the thicknesses or ! other properties where the ice is already in the correct categories. - call adjust_ice_categories(sIST%mH_ice, sIST%mH_snow, sIST%mH_pond, sIST%part_size, & + call adjust_ice_categories(sIST%mH_ice, sIST%mH_snow, sIST%mH_pond, sIST%mH_pond_ice, sIST%part_size, & sIST%TrReg, sG, sIG, SIS_dyn_trans_transport_CS(Ice%sCS%dyn_trans_CSp)) endif @@ -2512,6 +2516,7 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, call pass_var(sIST%mH_ice, sGD, complete=.false.) call pass_var(sIST%mH_snow, sGD, complete=.false.) call pass_var(sIST%mH_pond, sGD, complete=.false.) + call pass_var(sIST%mH_pond_ice, sGD, complete=.false.) do l=1,NkIce call pass_var(sIST%enth_ice(:,:,:,l), sGD, complete=.false.) enddo diff --git a/src/ice_pond.F90 b/src/ice_pond.F90 new file mode 100644 index 00000000..6adcc9bb --- /dev/null +++ b/src/ice_pond.F90 @@ -0,0 +1,513 @@ +!> A full implementation of Icepack ponds parameterizations (to come) +module ice_ponds_mod + +! This file is a part of SIS2. See LICENSE.md for the license. + +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! +! ! +! ! +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! + +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 +use MOM_error_handler, only : SIS_error=>MOM_error, FATAL, WARNING, SIS_mesg=>MOM_mesg +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 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_meltpond_lvl, only: compute_ponds_lvl +use icepack_meltpond_topo, only: compute_ponds_topo +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_ponds, ice_ponds_init + +!> Ice ponds control structure +type, public :: ice_ponds_CS ; private + logical :: & + level_pond = .false., & !< .true. = preferred ponds on level ice + topo_pond = .false. !< .true. = topographic ponds + real :: area_underflow = 0.0 !< a non-dimesional fractional area underflow limit for the sea-ice + !! ponding scheme. This is defaulted to zero, but a reasonable + !! value might be 10^-26 which for a km square grid cell + !! would equate to an Angstrom scale ice patch. +end type ice_ponds_CS + +contains + +!> Initialize the ice ponds +subroutine ice_ponds_init(G, IG, PF, CS, US) + type(SIS_hor_grid_type), intent(in) :: G !< G The ocean's grid structure. + type(ice_grid_type), intent(in) :: IG !< The sea-ice-specific grid structure. + type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters + type(ice_ponds_CS), pointer :: CS !< The ponds control structure. + 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_alvl, nt_vlvl, nt_qsno + character(len=40) :: mdl = "ice_ponds_init" ! This module's name. + + if (.not.associated(CS)) allocate(CS) + call get_param(PF, mdl, "MELTPOND_LEVEL", CS%level_pond, & + "Use level melt ponds", default=.false.) + call get_param(PF, mdl, "MELTPOND_TOPO", CS%topo_pond, & + "Use topographic melt ponds", default=.false.) + + 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 temperature + 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_alvl=2+2*nilyr+nslyr ! Index for level ice fraction +! nt_vlvl=3+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_alvl_in=nt_alvl, nt_vlvl_in=nt_vlvl ) + +! call icepack_init_parameters(mu_rdg_in=CS%mu_rdg, conserv_check_in=.true.) + +end subroutine ice_ponds_init + +!> ice_ponds is a wrapper for the Icepack pond routines +subroutine ice_ponds(IST, G, IG, mca_ice, mca_snow, mca_pond, TrReg, CS, US, dt, & + rdg_rate, rdg_height) + 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 !< mass of ice? + real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), intent(inout) :: mca_snow !< mass of snow? + real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), intent(inout) :: mca_pond !< mass of pond water? + type(SIS_tracer_registry_type), pointer :: TrReg !< TrReg - The registry of registered SIS ice and + !! snow tracers. + type(ice_ponds_CS), intent(in) :: CS !< The ponds control structure. + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors. + real, intent(in) :: dt !< The amount of time over which the ice dynamics are to be. + !! advanced in seconds. [T ~> s] + real, dimension(SZI_(G),SZJ_(G)), intent(out), optional :: rdg_rate !< Diagnostic of the rate of fractional + !! area loss-gain due to ridging (1/s) + real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), intent(inout), optional :: rdg_height !< A diagnostic of the ridged ice + !! height [Z ~> m] + +! logical, intent(in) :: dyn_Cgrid !< True if using C-grid velocities, B-grid if False. + + real :: dt_sec ! timestep in seconds + ! 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 :: & + ncat , & ! number of thickness categories + nilyr , & ! number of ice layers + nslyr ! number of snow layers + + real, dimension(0:IG%CatIce) :: hin_max ! category limits (m) + + logical :: & + closing_flag, &! flag if closing is valid + tr_brine ! if .true., brine height differs from ice thickness + + ! optional history fields + real :: & + 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, 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 (m) + aredistn , & ! redistribution function: fraction of new ridge area + vredistn ! redistribution function: fraction of new ridge volume (m) + + real, dimension(IG%CatIce) :: & + faero_ocn ! aerosol flux to ocean (kg/m^2/s) + + real, dimension(IG%CatIce) :: & + fiso_ocn ! isotope flux to ocean (kg/m^2/s) + + integer :: & + ndtd = 1 , & ! number of dynamics subcycles + n_aero = 0, & ! number of aerosol tracers + ntrcr = 0 ! number of tracer level + + real :: & + 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, 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 storage + ! ice tracers; ntr*(NkIce+NkSnow) guaranteed to be enough for all (intensive) + real, dimension(4+2*IG%NkIce+IG%NkSnow,IG%CatIce) :: trcrn + + real :: aice0 ! concentration of open water + + integer, dimension(4+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, dimension(4+2*IG%NkIce+IG%NkSnow,3) :: & + trcr_base ! = 0 or 1 depending on tracer dependency + ! argument 2: (1) aice, (2) vice, (3) vsno + + integer, dimension(4+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, dimension(:,:,:), pointer :: Tr_ice_alvl_ptr=>NULL() !< A pointer to the named tracer + real, dimension(:,:,:), pointer :: Tr_ice_mlvl_ptr=>NULL() !< A pointer to the named tracer + + real :: rho_ice, rho_snow ! Density of ice and snow [R ~> kg m-3] + real :: divu_adv + integer :: m, n ! loop vars for tracer; n is tracer #; m is tracer layer + integer :: nt_tsfc_in, nt_qice_in, nt_qsno_in, nt_sice_in + integer :: nL_ice, nL_snow ! number of tracer levels + integer :: ncat_out, ntrcr_out, nilyr_out, nslyr_out ! array sizes returned from Icepack query + character(len=1256) :: mesg + + 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) + dt_sec = dt*US%T_to_s + + 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, "Oops!! It looks like you are trying to use sea-ice ponds "//& + "but did not include the Icepack (https://github.com/CICE-Consortium/Icepack)"//& + "source code repository in your compilation procedure, and are instead using the default "//& + "stub routine contained in config_src/external. Adjust your compilation accordingly." ) + + ! copy strain calculation code from SIS_C_dynamics; might be a more elegant way ... + ! + halo_sh_Ds = min(isc-G%isd, jsc-G%jsd, 2) + ! if (dyn_Cgrid) then + do J=jsc-halo_sh_Ds,jec+halo_sh_Ds-1 ; do I=isc-halo_sh_Ds,iec+halo_sh_Ds-1 + ! This uses a no-slip boundary condition. + sh_Ds(I,J) = (2.0-G%mask2dBu(I,J)) * & + (G%dxBu(I,J)*G%IdyBu(I,J)*(IST%u_ice_C(I,j+1)*G%IdxCu(I,j+1) - IST%u_ice_C(I,j)*G%IdxCu(I,j)) + & + G%dyBu(I,J)*G%IdxBu(I,J)*(IST%v_ice_C(i+1,J)*G%IdyCv(i+1,J) - IST%v_ice_C(i,J)*G%IdyCv(i,J))) + enddo; enddo + + ! set category limits; Icepack has a max on the largest, unlimited, category (why?) + + hin_max(0)=0.0 + do k=1,nCat + hin_max(k) = US%Z_to_m * IG%mH_cat_bound(k) / Rho_ice + end do + + trcr_base = 0.0; n_trcr_strata = 0; nt_strata = 0 ! init some tracer vars + ! When would we use icepack tracer "strata"? + + ! set icepack tracer index "nt_lvl" to (last) pond tracer so it gets dumped when + ! ridging in ridge_ice (this is what happens to "level" ponds); first add up ntrcr; + ! then set nt_lvl to ntrcr+1; could move this to an initializer - mw + + call get_SIS_tracer_pointer("enth_ice", TrReg, Tr_ice_enth_ptr, nL_ice) + call get_SIS_tracer_pointer("enth_snow", TrReg, Tr_snow_enth_ptr, nL_snow) + call get_SIS_tracer_pointer("salin_ice", TrReg, Tr_ice_salin_ptr, nL_ice) +! call get_SIS_tracer_pointer("level_area", TrReg, Tr_ice_alvl_ptr, 1) +! call get_SIS_tracer_pointer("level_mass", TrReg, Tr_ice_mlvl_ptr, 1) + +! call IST_chksum('before ice ponds ', 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 ponds + 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 * US%Z_to_m ! volume per unit area of ice (m) + vsnon(1:nCat) = mca_snow(i,j,1:nCat)/Rho_snow * US%Z_to_m ! 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 ) )*US%s_to_T ! H&D eqn 9 + rdg_conv = -min(sh_Dd,0.0)*US%s_to_T ! energy dissipated by convergence ... + rdg_shear = 0.5*(del_sh-abs(sh_Dd))*US%s_to_T ! ... and by shear + + aice0 = IST%part_size(i,j,0) + if (aice0<0.) then + call SIS_error(WARNING, 'aice0<0. before call to pond ice.') + aice0=0. + endif + + ntrcr = 0 +! Tr_ptr=>NULL() + if (TrReg%ntr>0) then ! load tracer array + ntrcr=ntrcr+1 + do k=1,ncat + trcrn(ntrcr,k) = Tr_ice_enth_ptr(i,j,1,1) ! surface temperature taken from the ice-free category + ! copying across all categories. + enddo + trcr_depend(ntrcr) = 0 ! ice/snow surface temperature + trcr_base(ntrcr,:) = 0.0; trcr_base(ntrcr,1) = 1.0 ! 1st index for area + 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 +! ntrcr=ntrcr+1 +! trcrn(ntrcr,1:nCat) = Tr_ice_alvl_ptr(i,j,1:nCat,1) +! trcr_depend(ntrcr) = 0 ! 1 means area-based tracer +! trcr_base(ntrcr,:) = 0.0; trcr_base(ntrcr,1) = 1.0 ! 1st index for area +! ntrcr=ntrcr+1 +! trcrn(ntrcr,1:nCat) = Tr_ice_mlvl_ptr(i,j,1:nCat,1) +! trcr_depend(ntrcr) = 1 ! 1 means ice-based tracer +! trcr_base(ntrcr,:) = 0.0; trcr_base(ntrcr,2) = 1.0 ! 2nd index for ice + 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(:) = rdg_height(i,j,:)*US%Z_to_m + 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 compute_ponds_lvl (dt=dt, & +! nilyr=nilyr, & +! ktherm=ktherm, & +! hi_min=hi_min, & +! dpscale=dpscale, & +! frzpnd=frzpnd, & +! rfrac=rfrac, & +! meltt=melttn (n), & +! melts=meltsn (n), & +! frain=frain, & +! Tair=Tair, & +! fsurfn=fsurfn(n), & +! dhs=dhsn (n), & +! ffrac=ffracn (n), & +! aicen=aicen (n), & +! vicen=vicen (n), & +! vsnon=vsnon (n), & +! qicen=zqin (:,n), & +! sicen=zSin (:,n), & +! Tsfcn=Tsfc (n), & +! alvl=alvl (n), & +! apnd=apnd (n), & +! hpnd=hpnd (n), & +! ipnd=ipnd (n), & +! meltsliqn=l_meltsliqn(n)) + +! call compute_ponds_topo(dt, ncat, nilyr, & +! ktherm, & +! aice, aicen, & +! vice, vicen, & +! vsno, vsnon, & +! meltt, & +! fsurf, fpond, & +! Tsfc, Tf, & +! zqin, zSin, & +! apnd, hpnd, ipnd ) +! if (icepack_warnings_aborted(subname)) return + + if (present(rdg_rate)) rdg_rate(i,j) = (dardg1dt - dardg2dt)*US%T_to_s + if (present(rdg_height)) rdg_height(i,j,:) = krdgn(:)*US%m_to_Z + + if ( icepack_warnings_aborted() ) then + call icepack_warnings_flush(0) + call icepack_warnings_setabort(.false.) + call SIS_error(WARNING, 'icepack compute_ponds error') + endif + + ! pop pond off top of stack + tr_tmp(1:nCat)=trcrn(ntrcr,1:nCat) + + do k=1,nCat + IST%mH_pond(i,j,k) = tr_tmp(k) + mca_pond(i,j,k) = IST%mH_pond(i,j,k)*aicen(k) + enddo + if (any(vicen < 0)) then +! print *, "Negative ice volume after ponds: ", i+G%idg_offset, j+G%jdg_offset, vicen +! print *, "Before ponds: ", mca_ice(i,j,1:nCat) /Rho_ice +! print *, "Ice concentration before/after ponds: ", IST%part_size(i,j,1:nCat), aicen + do k=1,nCat + if (vicen(k) < 0.0 .and. aicen(k) > 0.0) then + write(mesg,'("Negative ice volume after ponds: ", i6, i6, 2x, 1pe12.4, 1pe12.4)') & + i+G%idg_offset, j+G%jdg_offset, aicen(k), vicen(k) + call SIS_error(WARNING, mesg, all_print=.true.) + endif + vicen(k) = max(vicen(k),0.0) + enddo +! write(mesg,'("Negative ice volume after ponds: ", 2i6, 2x, (1pe12.4))') & +! i+G%jdg_offset, j+G%jdg_offset, aicen, vicen +! call SIS_error(WARNING, mesg, all_print=.true.) + endif + + if (TrReg%ntr>0) then + ! unload tracer array reversing order of load -- stack-like fashion + +! tr_tmp(1:nCat)=trcrn(ntrcr-1,1:nCat) +! Tr_ice_mlvl_ptr(i,j,1:nCat,1) = tr_tmp(1:nCat) +! tr_tmp(1:nCat)=trcrn(ntrcr-2,1:nCat) +! Tr_ice_alvl_ptr(i,j,1:nCat,1) = tr_tmp(1:nCat) + + 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 + + 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 + + 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) < CS%area_underflow) then + aicen(k)=0.0 + vicen(k)=0.0 + endif + if (aicen(k) > 0.0) then + IST%part_size(i,j,k) = aicen(k) + mca_ice(i,j,k) = vicen(k)*Rho_ice * US%m_to_Z + IST%mH_ice(i,j,k) = vicen(k)*Rho_ice/aicen(k) * US%m_to_Z + mca_snow(i,j,k) = vsnon(k)*Rho_snow * US%m_to_Z + IST%mH_snow(i,j,k) = vsnon(k)*Rho_snow/aicen(k) * US%m_to_Z + 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 ponds ', IST, G, US, IG) + +end subroutine ice_ponds + +!> ice_ponds_end deallocates the memory associated with this module. +subroutine ice_ponds_end() + +end subroutine ice_ponds_end + +end module ice_ponds_mod diff --git a/src/ice_ridge.F90 b/src/ice_ridge.F90 index 6c710512..ecae7fa5 100644 --- a/src/ice_ridge.F90 +++ b/src/ice_ridge.F90 @@ -276,13 +276,13 @@ subroutine ice_ridging(IST, G, IG, mca_ice, mca_snow, mca_pond, TrReg, CS, US, d call get_SIS2_thermo_coefs(IST%ITV, Cp_Water=Cp_water) dt_sec = dt*US%T_to_s - call icepack_query_tracer_sizes(ncat_out=ncat_out,ntrcr_out=ntrcr_out, nilyr_out=nilyr_out, nslyr_out=nslyr_out) + 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,"Oops!! It looks like you are trying to use sea-ice ridging "//& - "but did not include the Icepack (https://github.com/CICE-Consortium/Icepack)"//& - "source code repository in your compilation procedure, and are instead using the default "//& - "stub routine contained in config_src/external. Adjust your compilation accordingly." ) + call SIS_error(FATAL, "Oops!! It looks like you are trying to use sea-ice ridging "//& + "but did not include the Icepack (https://github.com/CICE-Consortium/Icepack)"//& + "source code repository in your compilation procedure, and are instead using the default "//& + "stub routine contained in config_src/external. Adjust your compilation accordingly." ) ! copy strain calculation code from SIS_C_dynamics; might be a more elegant way ... ! @@ -305,20 +305,20 @@ subroutine ice_ridging(IST, G, IG, mca_ice, mca_snow, mca_pond, TrReg, CS, US, d ! is thinner than hin_max(ncat). hin_max(nCat) = big - trcr_base = 0.0; n_trcr_strata = 0; nt_strata = 0; ! init some tracer vars + trcr_base = 0.0; n_trcr_strata = 0; nt_strata = 0 ! init some tracer vars ! When would we use icepack tracer "strata"? ! set icepack tracer index "nt_lvl" to (last) pond tracer so it gets dumped when ! ridging in ridge_ice (this is what happens to "level" ponds); first add up ntrcr; ! then set nt_lvl to ntrcr+1; could move this to an initializer - mw - call get_SIS_tracer_pointer("enth_ice",TrReg,Tr_ice_enth_ptr,nL_ice) - call get_SIS_tracer_pointer("enth_snow",TrReg,Tr_snow_enth_ptr,nL_snow) - call get_SIS_tracer_pointer("salin_ice",TrReg,Tr_ice_salin_ptr,nL_ice) + call 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 get_SIS_tracer_pointer("level_area",TrReg,Tr_ice_alvl_ptr,1) ! call get_SIS_tracer_pointer("level_mass",TrReg,Tr_ice_mlvl_ptr,1) -! call IST_chksum('before ice ridging ',IST,G,US,IG) +! 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 @@ -565,7 +565,7 @@ subroutine ice_ridging(IST, G, IG, mca_ice, mca_snow, mca_pond, TrReg, CS, US, d endif ! part_sz enddo; enddo ! j, i -! call IST_chksum('after ice ridging ',IST,G,US,IG) +! call IST_chksum('after ice ridging ', IST, G, US, IG) end subroutine ice_ridging