diff --git a/src/SIS2_ice_thm.F90 b/src/SIS2_ice_thm.F90 index 5767731..e555975 100644 --- a/src/SIS2_ice_thm.F90 +++ b/src/SIS2_ice_thm.F90 @@ -1008,8 +1008,9 @@ end subroutine ice_check subroutine ice_resize_SIS2(a_ice, m_pond, 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, & - enthalpy_evap, enthalpy_melt, enthalpy_freeze) + snow_to_ice, s2i_i, s2i_s, salt_to_ice, ITV, US, CS, ablation, & + ablation_i, ablation_s, enthalpy_evap, enthalpy_melt, enthalpy_freeze, & + evap_i, evap_s) ! 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] @@ -1023,7 +1024,7 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, & intent(inout) :: Salin !< Conserved ice bulk salinity by layer [S ~> gSalt kg-1] real, intent(in ) :: snow !< new snow [R Z ~> kg m-2] real, intent(in ) :: rain !< rain for pond source [R Z ~> kg m-2] - not yet active - real, intent(in ) :: evap !< ice evaporation/sublimation [R Z ~> kg m-2] + real, intent(in ) :: evap !< ice and snow evaporation/sublimation [R Z ~> kg m-2] real, intent(in ) :: tmlt !< top melting energy [Q R Z ~> J m-2] real, intent(in ) :: bmlt !< bottom melting energy [Q R Z ~> J m-2] integer, intent(in) :: NkIce !< The number of ice layers. @@ -1035,18 +1036,24 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, & real, intent( out) :: h2o_ocn_to_ice !< liquid water flux from ocean [R Z ~> kg m-2] real, intent( out) :: evap_from_ocn !< evaporation flux from ocean [R Z ~> kg m-2] real, intent( out) :: snow_to_ice !< snow below waterline becomes ice [R Z ~> kg m-2] + real, intent( out) :: s2i_i !< snow below waterline becomes ice (ice flux) [R Z ~> kg m-2] + real, intent( out) :: s2i_s !< snow below waterline becomes ice (snow flux) [R Z ~> kg m-2] real, intent( out) :: salt_to_ice !< Net flux of salt to the ice [R Z S ~> gSalt m-2]. type(ice_thermo_type), intent(in) :: ITV !< The ice thermodynamic parameter structure. type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS2_ice_thm_CS), intent(in) :: CS !< The SIS2_ice_thm control structure. - real, intent( out) :: ablation !< The mass loss from bottom melt [R Z ~> kg m-2]. + real, intent( out) :: ablation !< The mass loss from bottom melt of ice and snow [R Z ~> kg m-2]. + real, intent( out) :: ablation_i !< The mass loss from bottom melt of ice [R Z ~> kg m-2]. + real, intent( out) :: ablation_s !< The mass loss from bottom melt of snow [R Z ~> kg m-2]. real, intent( out) :: enthalpy_evap !< The enthalpy loss due to the mass loss !! by evaporation / sublimation. [Q R Z ~> J m-2] real, intent( out) :: enthalpy_melt !< The enthalpy loss due to the mass loss !! by melting [Q R Z ~> J m-2]. real, intent( out) :: enthalpy_freeze !< The enthalpy gain due to the mass gain !! by freezing [Q R Z ~> J m-2]. + real, intent( out) :: evap_i !< ice evaporation/sublimation [R Z ~> kg m-2] + real, intent( out) :: evap_s !< snow evaporation/sublimation [R Z ~> kg m-2] real :: top_melt, bot_melt, melt_left ! Heating amounts, all in [Q R Z ~> J m-2] real :: mtot_ice ! The summed ice mass [R Z ~> kg m-2]. @@ -1097,9 +1104,11 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, & evap_from_ocn = 0.0 ! for excess evap-melt h2o_ocn_to_ice = 0.0 ; h2o_ice_to_ocn = 0.0 ; snow_to_ice = 0.0 + s2i_i = 0.0 ; s2i_s = 0.0 h2o_to_pond = 0.0 h2o_from_pond = 0.0 salt_to_ice = 0.0 + evap_s = 0.0 ; evap_i = 0.0 enthM_freezing = 0.0 ; enthM_melt = 0.0 ; enthM_evap = 0.0 ; enthM_snowfall = 0.0 ! raining on cold ice led to unphysical temperature oscillations @@ -1117,6 +1126,7 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, & if (evap < 0.0) then m_lay(0) = m_lay(0) - evap ! Treat frost formation like snow. enthM_snowfall = enthM_snowfall - evap*enthalpy(0) + evap_s = evap_s + evap endif if (top_melt < 0.0 .and. CS%do_pond) then ! mw/new: add fresh/0C ice to top layer @@ -1185,6 +1195,11 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, & evap_here = min(evap_left, m_lay(k)) evap_left = evap_left - evap_here m_lay(k) = m_lay(k) - evap_here + if (k == 0) then + evap_s = evap_s + evap_here + else + evap_i = evap_i + evap_here + endif ! Assume that evaporation does not make ice salty? if (k>0) salt_to_ice = salt_to_ice - Salin(k) * evap_here enthM_evap = enthM_evap + evap_here * enthalpy(k) @@ -1229,6 +1244,7 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, & ! apply bottom melt heat flux melt_left = bot_melt ; ablation = 0.0 + ablation_i = 0.0 ; ablation_s = 0.0 if (melt_left > 0.0) then ! melt ice and snow from below do k=NkIce,0,-1 if (melt_left < m_lay(k) * (enth_fr(k) - Enthalpy(k))) then @@ -1243,7 +1259,11 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, & h2o_ice_to_ocn = h2o_ice_to_ocn + M_melt enthM_melt = enthM_melt + M_melt*enth_fr(k) ablation = ablation + M_melt - + if (k > 0) then + ablation_i = ablation_i + M_melt + else + ablation_s = ablation_s + M_melt + endif if (melt_left <= 0.0) exit ! All melt energy has been used. enddo @@ -1287,6 +1307,7 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, & snow_to_ice = min(m_submerged - mtot_ice, m_lay(0)) ! need ice from snow m_lay(0) = m_lay(0) - snow_to_ice + s2i_s = s2i_s - snow_to_ice ! Add ice to the topmost layer and dilute its salinity. Enthalpy(1) = (m_lay(1)*Enthalpy(1) + snow_to_ice*Enthalpy(0)) / & @@ -1296,6 +1317,7 @@ subroutine ice_resize_SIS2(a_ice, m_pond, m_lay, Enthalpy, Sice_therm, Salin, & TrLay(1,tr) = TrLay(1,tr) * m_lay(1) / (m_lay(1) + snow_to_ice) enddo m_lay(1) = m_lay(1) + snow_to_ice + s2i_i = s2i_i + snow_to_ice else snow_to_ice = 0.0 endif diff --git a/src/SIS_ctrl_types.F90 b/src/SIS_ctrl_types.F90 index 041d073..cc057bd 100644 --- a/src/SIS_ctrl_types.F90 +++ b/src/SIS_ctrl_types.F90 @@ -238,6 +238,17 @@ subroutine ice_diagnostics_init(IOF, OSS, FIA, G, US, IG, diag, Time, Cgrid) 'frozen runoff sensible heat flux', 'W/m^2', conversion=US%QRZ_T_to_W_m2, missing_value=missing) FIA%id_evap = register_SIS_diag_field('ice_model', 'EVAP',diag%axesT1, Time, & 'evaporation', 'kg/(m^2*s)', conversion=US%RZ_T_to_kg_m2s, missing_value=missing) + + !CMOR evaporation diagnostics + FIA%id_evap_i = register_SIS_diag_field('ice_model', 'sidmassevapsubl',diag%axesT1, Time, & + 'Sea-Ice Mass Change Through Evaporation and Sublimation', 'kg m-2 s-1', & + conversion=US%RZ_T_to_kg_m2s, missing_value=missing, & + cmor_standard_name='water_evapotranspiration_flux') + FIA%id_evap_s = register_SIS_diag_field('ice_model', 'sisndmasssubl',diag%axesT1, Time, & + 'Snow Mass Rate of Change Through Evaporation or Sublimation', 'kg m-2 s-1', & + conversion=US%RZ_T_to_kg_m2s, missing_value=missing, & + cmor_standard_name='tendency_of_atmosphere_mass_content_of_water_vapor_due_to_sublimation_of_surface_snow_and_ice') + IOF%id_saltf = register_SIS_diag_field('ice_model', 'SALTF', diag%axesT1, Time, & 'ice to ocean salt flux', 'kg/(m^2*s)', conversion=US%S_to_ppt*US%RZ_T_to_kg_m2s, missing_value=missing) FIA%id_tmelt = register_SIS_diag_field('ice_model', 'TMELT', diag%axesT1, Time, & @@ -337,6 +348,12 @@ subroutine ice_diagnostics_init(IOF, OSS, FIA, G, US, IG, diag, Time, Cgrid) OSS%id_frazil = register_SIS_diag_field('ice_model', 'FRAZIL', diag%axesT1, Time, & 'energy flux of frazil formation', 'W/m^2', conversion=US%QRZ_T_to_W_m2, missing_value=missing) + + !CMOR frazil diagnostic + OSS%id_frazilmass = register_SIS_diag_field('ice_model', 'sidmassgrowthwat', diag%axesT1, Time, & + 'Sea-Ice Mass Change Through Growth in Supercooled Open Water (Frazil)', 'kg m-2 s-1', & + conversion=US%RZ_to_kg_m2, missing_value=missing, & + cmor_standard_name='tendency_of_sea_ice_amount_due_to_frazil_ice_accumulation_in_leads') if (coupler_type_initialized(OSS%tr_fields)) & call coupler_type_set_diags(OSS%tr_fields, 'ice_model', diag%axesT1%handles, Time) diff --git a/src/SIS_ice_diags.F90 b/src/SIS_ice_diags.F90 index c549b04..6382f50 100644 --- a/src/SIS_ice_diags.F90 +++ b/src/SIS_ice_diags.F90 @@ -222,7 +222,7 @@ subroutine post_ice_state_diagnostics(IDs, IST, OSS, IOF, dt_slow, Time, G, US, ! Write out diagnostics of the ocean surface state, as seen by the slow sea ice. ! These fields do not change over the course of the sea-ice time stepping. - call post_ocean_sfc_diagnostics(OSS, dt_slow, Time, G, diag) + call post_ocean_sfc_diagnostics(OSS, dt_slow, Time, G, diag, IST) if (IDs%id_e2m>0) then tmp2d(:,:) = 0.0 @@ -264,15 +264,18 @@ end subroutine post_ice_state_diagnostics !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> Offer diagnostics of the ocean surface field, as seen by the sea ice. -subroutine post_ocean_sfc_diagnostics(OSS, dt_slow, Time, G, diag) +subroutine post_ocean_sfc_diagnostics(OSS, dt_slow, Time, G, diag, IST) type(ocean_sfc_state_type), intent(in) :: OSS !< A structure containing the arrays that describe !! the ocean's surface state for the ice model. real, intent(in) :: dt_slow !< The time interval of these diagnostics [T ~> s] type(time_type), intent(in) :: Time !< The ending time of these diagnostics type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type type(SIS_diag_ctrl), pointer :: diag !< A structure that is used to regulate diagnostic output + type(ice_state_type), optional,intent(in) :: IST !< A type describing the state of the sea ice - real :: Idt_slow ! The inverse of the thermodynamic step [T-1 ~> s-1]. + real :: Idt_slow ! The inverse of the thermodynamic step [T-1 ~> s-1]. + real :: LatHtFus ! The latent heat of fusion of ice [Q ~> J kg-1]. + real :: ILatHtFus_s ! The inverse of latent heat of fusion of ice, per second [Q-1 ~> kg J-1 s-1]. Idt_slow = 0.0 ; if (dt_slow > 0.0) Idt_slow = 1.0/dt_slow ! Write out diagnostics of the ocean surface state, as seen by the slow sea ice. @@ -289,6 +292,11 @@ subroutine post_ocean_sfc_diagnostics(OSS, dt_slow, Time, G, diag) endif if (OSS%id_frazil>0) & call post_data(OSS%id_frazil, OSS%frazil*Idt_slow, diag) + if (OSS%id_frazilmass>0) then + call get_SIS2_thermo_coefs(IST%ITV, Latent_fusion=LatHtFus) + ILatHtFus_s = 1.0 / (LatHtFus * dt_slow) + call post_data(OSS%id_frazilmass, OSS%frazil*ILatHtFus_s, diag) + endif if (coupler_type_initialized(OSS%tr_fields)) & call coupler_type_send_data(OSS%tr_fields, Time) diff --git a/src/SIS_slow_thermo.F90 b/src/SIS_slow_thermo.F90 index a2d224e..b102622 100644 --- a/src/SIS_slow_thermo.F90 +++ b/src/SIS_slow_thermo.F90 @@ -141,6 +141,11 @@ module SIS_slow_thermo !>@{ Diagnostic IDs integer :: id_qflim=-1, id_qfres=-1, id_fwnudge=-1, id_net_melt=-1, id_CMOR_melt=-1 integer :: id_lsrc=-1, id_lsnk=-1, id_bsnk=-1, id_sn2ic=-1 + integer :: id_lsrc_i=-1, id_lsnk_i=-1, id_bsnk_i=-1 + integer :: id_lsrc_s=-1, id_lsnk_s=-1, id_bsnk_s=-1 + integer :: id_lsrc_c=-1, id_lsnk_c=-1 + integer :: id_bsnk_i_cmor=-1, id_tsnk_i=-1, id_bsrc_i=-1, id_net_i=-1, id_net_s=-1 + integer :: id_net_c=-1, id_sn2ic_i=-1, id_sn2ic_s=-1 !!@} end type slow_thermo_CS @@ -319,6 +324,9 @@ subroutine slow_thermodynamics(IST, dt_slow, CS, OSS, FIA, XSF, IOF, G, US, IG, ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: & h_ice_input ! The specified ice thickness, with specified_ice [m]. + real, dimension(SZI_(G),SZJ_(G)) :: & + tmp2d, & + h2o_change_c ! The change in ice area due to thermodynamics [s-1] real :: rho_ice ! The nominal density of sea ice [R ~> kg m-3]. real :: Idt_slow ! The inverse of the slow thermodynamic time step [T-1 ~> s-1] @@ -331,6 +339,7 @@ subroutine slow_thermodynamics(IST, dt_slow, CS, OSS, FIA, XSF, IOF, G, US, IG, real :: mass_part ! The mass per unit cell area in a thickness category [R Z ~> kg m-2] mi_old(:,:,:) = 0.0 + h2o_change_c(:,:) = 0.0 isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; ncat = IG%CatIce isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; NkIce = IG%NkIce nb = size(FIA%flux_sw_top,4) @@ -366,6 +375,9 @@ subroutine slow_thermodynamics(IST, dt_slow, CS, OSS, FIA, XSF, IOF, G, US, IG, ! ! Thermodynamics ! + do j=jsc,jec ; do k=1,ncat ; do i=isc,iec + h2o_change_c(i,j) = h2o_change_c(i,j) - IST%part_size(i,j,k) + enddo ; enddo ; enddo if (CS%specified_ice) then ! over-write changes with specifications. h_ice_input(:,:) = 0.0 call get_sea_surface(CS%Time, G%HI, SST=OSS%SST_C, ice_conc=IST%part_size(:,:,1), ice_thick=h_ice_input, US=US) @@ -469,6 +481,39 @@ subroutine slow_thermodynamics(IST, dt_slow, CS, OSS, FIA, XSF, IOF, G, US, IG, call adjust_ice_categories(IST%mH_ice, IST%mH_snow, IST%mH_pond, IST%part_size, & IST%TrReg, G, IG, CS%SIS_transport_CSp) !Niki: add ridging? + do j=jsc,jec ; do k=1,ncat ; do i=isc,iec + h2o_change_c(i,j) = h2o_change_c(i,j) + IST%part_size(i,j,k) + enddo ; enddo ; enddo + + do j=jsc,jec ; do k=1,ncat ; do i=isc,iec + h2o_change_c(i,j) = h2o_change_c(i,j) + IST%part_size(i,j,k) + enddo; enddo ; enddo + + call enable_SIS_averaging(US%T_to_s*dt_slow, CS%Time, CS%diag) + if (CS%id_lsrc_c>0) then + !$OMP parallel do default(shared) + do j=jsc,jec ; do i=isc,iec + tmp2d(i,j) = max(h2o_change_c(i,j),0.0) * Idt_slow + enddo ; enddo + call post_data(CS%id_lsrc_c, tmp2d(isc:iec,jsc:jec), CS%diag) + endif + if (CS%id_lsnk_c>0) then + !$OMP parallel do default(shared) + do j=jsc,jec ; do i=isc,iec + tmp2d(i,j) = min(h2o_change_c(i,j),0.0) * Idt_slow + enddo ; enddo + call post_data(CS%id_lsnk_c, tmp2d(isc:iec,jsc:jec), CS%diag) + endif + if (CS%id_net_c>0) then + !$OMP parallel do default(shared) + do j=jsc,jec ; do i=isc,iec + tmp2d(i,j) = h2o_change_c(i,j) * Idt_slow + enddo ; enddo + call post_data(CS%id_net_c, tmp2d(isc:iec,jsc:jec), CS%diag) + endif + call SIS_diag_send_complete() + call disable_SIS_averaging(CS%diag) + if (CS%column_check) & call write_ice_statistics(IST, CS%Time, CS%n_calls, G, US, IG, CS%sum_output_CSp, & message=" Post_thermo B ", check_column=.true.) @@ -573,14 +618,24 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, US, IG) ! but with a greater emphasis on enthalpy as the dominant state variable. real, dimension(SZI_(G),SZJ_(G),1:IG%CatIce) :: snow_to_ice ! The flux of snow to the ice [R Z ~> kg m-2] + real, dimension(SZI_(G),SZJ_(G),0:IG%CatIce) :: snow_to_ice_i ! The ice flux of snow to the ice, including + ! open ocean category [R Z ~> kg m-2] + real, dimension(SZI_(G),SZJ_(G),0:IG%CatIce) :: snow_to_ice_s ! The snow flux of snow to the ice, including + ! open ocean category [R Z ~> kg m-2] real, dimension(G%isc:G%iec,G%jsc:G%jec) :: Obs_h_ice ! Observed ice thickness for qflux calculation real, dimension(G%isc:G%iec,G%jsc:G%jec) :: Obs_cn_ice ! Observed total ice concentration [nondim] real, dimension(G%isc:G%iec,G%jsc:G%jec) :: icec ! Total ice concentration [nondim] real, dimension(SZI_(G),SZJ_(G)) :: & salt_change, & ! The change in integrated salinity [R Z S ~> gSalt m-2] - h2o_change, & ! The change in water in the ice [R Z ~> kg m-2] - bsnk, & ! The bottom melting mass sink [R Z T-1 ~> kg m-2 s-1] - tmp2d, & ! A temporary array for mass balance diagnostics [R Z yr-1 ~> kg m-2 yr-1] + h2o_change, & ! The change in water in the ice and snow [R Z ~> kg m-2] + h2o_change_i, & ! The change in water in the ice [R Z ~> kg m-2] + h2o_change_s, & ! The change in water in the snow [R Z ~> kg m-2] + bsnk, & ! The bottom melting mass sink of ice and snow [R Z T-1 ~> kg m-2 s-1] + bsnk_i, & ! The bottom melting mass sink of ice [R Z T-1 ~> kg m-2 s-1] + bsnk_s, & ! The bottom melting mass sink of snow [R Z T-1 ~> kg m-2 s-1] + simass_evap, & ! The mass of ice evaporated from the surface [R Z T-1 ~> kg m-2 s-1] + sisnmass_evap, & ! The mass of snow evaporated from the surface [R Z T-1 ~> kg m-2 s-1] + tmp2d, & ! A temporary array for mass balance diagnostics [R Z yr-1 ~> kg m-2 s-1] qflx_lim_ice, & ! Ice limiting heat flux [Q R Z T-1 ~> W m-2] qflx_res_ice, & ! Ice restoring heat flux [Q R Z T-1 ~> W m-2] cool_nudge, & ! A heat flux out of the sea ice that @@ -632,8 +687,8 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, US, IG) real :: rho_ice ! The nominal density of sea ice [R ~> kg m-3]. real :: Idt_slow ! The inverse of the thermodynamic step [T-1 ~> s-1]. - real :: yr_dtslow ! The ratio of 1 year to the thermodynamic time step times some scaling - ! factors, used to change the units of several diagnostics to rate yr-1. + real :: yr_dtslow ! The ratio of 1 year to the thermodynamic time step times + ! some scaling factors. Changes units to [yr-1] real :: heat_to_ocn ! The heat passed from the ice to the ocean [Q R Z ~> J m-2] real :: water_to_ocn ! The water passed to the ocean [R Z ~> kg m-2] real :: salt_to_ocn ! The salt passed to the ocean [R Z S ~> gSalt m-2] @@ -644,8 +699,12 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, US, IG) real :: snow_loss ! The loss of snow mass from transmutation [R Z ~> kg m-2] real :: h2o_ice_to_ocn ! The downward water flux from the ice to the ocean [R Z ~> kg m-2] real :: h2o_ocn_to_ice ! The upward water flux from the ocean to the ice [R Z ~> kg m-2] + real :: evap_i ! The loss of ice due to evaporation [R Z ~> kg m-2] + real :: evap_s ! The loss of snow due to evaporation [R Z ~> kg m-2] real :: evap_from_ocn ! The evaporation from the ocean [R Z ~> kg m-2] - real :: bablt ! The bottom ablation ice loss [R Z ~> kg m-2] + real :: bablt ! The bottom ablation ice and snow loss [R Z ~> kg m-2] + real :: bablt_i ! The bottom ablation ice loss [R Z ~> kg m-2] + real :: bablt_s ! The bottom ablation snow loss [R Z ~> kg m-2] real :: salt_to_ice ! The flux of salt from the ocean to the ice [R Z S ~> gSalt m-2]. ! This may be of either sign; in some places it is an ! average over the whole cell, while in others just a partition. @@ -675,9 +734,10 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, US, IG) integer :: i, j, k, l, m, n, b, nb, isc, iec, jsc, jec, ncat, NkIce, tr, npassive integer :: k_merge real :: LatHtFus ! The latent heat of fusion of ice [Q ~> J kg-1]. + real :: ILatHtFus ! The inverse latent heat of fusion of ice [Q-1 ~> kg J-1]. real :: LatHtVap ! The latent heat of vaporization of water at 0C [Q ~> J kg-1]. - real :: tot_heat_in, enth_here, enth_imb, norm_enth_imb, emic2, tot_heat_in2, enth_imb2 + real :: tot_heat_in, enth_here, enth_imb, norm_enth_imb, emic2, tot_heat_in2, enth_imb2, s2i_agg isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; ncat = IG%CatIce NkIce = IG%NkIce ; I_Nk = 1.0 / NkIce @@ -687,6 +747,7 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, US, IG) call get_SIS2_thermo_coefs(IST%ITV, ice_salinity=S_col, & rho_ice=rho_ice, spec_thermo_salin=spec_thermo_sal, & Latent_fusion=LatHtFus, Latent_vapor=LatHtVap) + ILatHtFus = 1.0 / LatHtFus S_col0(0) = 0.0 ; do m=1,NkIce ; S_col0(m) = S_col(m) ; enddo heat_fill_val = Enth_from_TS(0.0, 0.0, IST%ITV) @@ -830,9 +891,16 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, US, IG) call cpu_clock_begin(iceClock5) snow_to_ice(:,:,:) = 0.0 ; net_melt(:,:) = 0.0 + snow_to_ice_i(:,:,:) = 0.0 ; snow_to_ice_s(:,:,:) = 0.0 bsnk(:,:) = 0.0 + bsnk_i(:,:) = 0.0 + bsnk_s(:,:) = 0.0 + simass_evap(:,:) = 0.0 + sisnmass_evap(:,:) = 0.0 salt_change(:,:) = 0.0 h2o_change(:,:) = 0.0 + h2o_change_i(:,:) = 0.0 + h2o_change_s(:,:) = 0.0 salt_left_in_ocean(:,:) = 0.0 !$OMP parallel default(shared) private(part_ocn) if (CS%ice_rel_salin <= 0.0) then @@ -846,6 +914,8 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, US, IG) do j=jsc,jec ; do k=1,ncat ; do i=isc,iec h2o_change(i,j) = h2o_change(i,j) - IST%part_size(i,j,k) * & (IST%mH_snow(i,j,k) + IST%mH_ice(i,j,k)) + h2o_change_i(i,j) = h2o_change_i(i,j) - IST%part_size(i,j,k) * IST%mH_ice(i,j,k) + h2o_change_s(i,j) = h2o_change_s(i,j) - IST%part_size(i,j,k) * IST%mH_snow(i,j,k) enddo ; enddo ; enddo ! Start accumulating the fluxes at the ocean's surface. @@ -888,11 +958,13 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, US, IG) !$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,G,US,IST,S_col0,NkIce,S_col,dt_slow, & !$OMP snow_to_ice,heat_in,I_NK,enth_prev,enth_mass_in_col,bsnk, & -!$OMP Idt_slow,salt_change,net_melt,LatHtFus,LatHtVap,IG,CS,OSS, & +!$OMP bsnk_i,bsnk_s,sisnmass_evap,simass_evap,Idt_slow, & +!$OMP salt_change,net_melt,LatHtFus,LatHtVap,IG,CS,OSS, & !$OMP FIA,IOF,npassive,nb,salt_left_in_ocean) & !$OMP private(mass_prev,enthalpy,enthalpy_ocean,Salin, & !$OMP heat_to_ocn,h2o_ice_to_ocn,h2o_ocn_to_ice, & -!$OMP evap_from_ocn,salt_to_ice,bablt,enth_evap, & +!$OMP evap_from_ocn,salt_to_ice,bablt,bablt_i, & +!$OMP bablt_s,enth_evap,evap_i,evap_s, & !$OMP enth_ice_to_ocn,enth_ocn_to_ice,heat_input, & !$OMP heat_mass_in,mass_in,mass_here,enth_here, & !$OMP tot_heat_in,enth_imb,mass_imb,norm_enth_imb, & @@ -935,8 +1007,9 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, US, IG) 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, & - snow_to_ice(i,j,k), salt_to_ice, IST%ITV, US, CS%ice_thm_CSp, bablt, & - enth_evap, enth_ice_to_ocn, enth_ocn_to_ice) + snow_to_ice(i,j,k), snow_to_ice_i(i,j,k), snow_to_ice_s(i,j,k), salt_to_ice, & + IST%ITV, US, CS%ice_thm_CSp, bablt, bablt_i, bablt_s, enth_evap, & + enth_ice_to_ocn, enth_ocn_to_ice, evap_i, evap_s) IST%mH_snow(i,j,k) = m_lay(0) call rebalance_ice_layers(m_lay, mtot_ice, Enthalpy, Salin, NkIce, npassive, TrLay) @@ -1028,7 +1101,11 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, US, IG) (sw_tot * FIA%sw_abs_ocn(i,j,k)) net_melt(i,j) = net_melt(i,j) + IST%part_size(i,j,k) * & ((h2o_ice_to_ocn-h2o_ocn_to_ice)*Idt_slow) - bsnk(i,j) = bsnk(i,j) - IST%part_size(i,j,k)*bablt*Idt_slow ! bot. melt. ablation + bsnk(i,j) = bsnk(i,j) - IST%part_size(i,j,k)*bablt*Idt_slow ! bot. melt. ice and snow ablation + bsnk_i(i,j) = bsnk_i(i,j) - IST%part_size(i,j,k)*bablt_i*Idt_slow ! bot. melt. ice ablation + bsnk_s(i,j) = bsnk_s(i,j) - IST%part_size(i,j,k)*bablt_s*Idt_slow ! bot. melt. snow ablation + simass_evap(i,j) = simass_evap(i,j) - IST%part_size(i,j,k)*evap_i*Idt_slow + sisnmass_evap(i,j) = sisnmass_evap(i,j) - IST%part_size(i,j,k)*evap_s*Idt_slow endif ! Applying surface fluxes to each category. enddo ; enddo ; enddo @@ -1305,6 +1382,8 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, US, IG) do k=1,ncat ; do i=isc,iec h2o_change(i,j) = h2o_change(i,j) + IST%part_size(i,j,k) * & (IST%mH_snow(i,j,k) + IST%mH_ice(i,j,k)) + h2o_change_i(i,j) = h2o_change_i(i,j) + IST%part_size(i,j,k) * IST%mH_ice(i,j,k) + h2o_change_s(i,j) = h2o_change_s(i,j) + IST%part_size(i,j,k) * IST%mH_snow(i,j,k) enddo ; enddo do i=isc,iec ! Note the conversion here from g m-2 to kg m-2 s-1. @@ -1365,6 +1444,20 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, US, IG) enddo ; enddo call post_data(CS%id_lsnk, tmp2d(isc:iec,jsc:jec), CS%diag) endif + if (CS%id_lsnk_i>0) then + !$OMP parallel do default(shared) + do j=jsc,jec ; do i=isc,iec + tmp2d(i,j) = min(h2o_change_i(i,j),0.0) * Idt_slow + enddo ; enddo + call post_data(CS%id_lsnk_i, tmp2d(isc:iec,jsc:jec), CS%diag) + endif + if (CS%id_lsnk_s>0) then + !$OMP parallel do default(shared) + do j=jsc,jec ; do i=isc,iec + tmp2d(i,j) = min(h2o_change_s(i,j),0.0) * Idt_slow + enddo ; enddo + call post_data(CS%id_lsnk_s, tmp2d(isc:iec,jsc:jec), CS%diag) + endif if (CS%id_lsrc>0) then !$OMP parallel do default(shared) do j=jsc,jec ; do i=isc,iec @@ -1372,8 +1465,61 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, US, IG) enddo ; enddo call post_data(CS%id_lsrc, tmp2d(isc:iec,jsc:jec), CS%diag) endif + if (CS%id_lsrc_i>0) then + !$OMP parallel do default(shared) + do j=jsc,jec ; do i=isc,iec + tmp2d(i,j) = max(h2o_change_i(i,j),0.0) * Idt_slow + enddo ; enddo + call post_data(CS%id_lsrc_i, tmp2d(isc:iec,jsc:jec), CS%diag) + endif + if (CS%id_lsrc_s>0) then + !$OMP parallel do default(shared) + do j=jsc,jec ; do i=isc,iec + tmp2d(i,j) = max(h2o_change_s(i,j),0.0) * Idt_slow + enddo ; enddo + call post_data(CS%id_lsrc_s, tmp2d(isc:iec,jsc:jec), CS%diag) + endif + if (CS%id_net_i>0) then + !$OMP parallel do default(shared) + do j=jsc,jec ; do i=isc,iec + tmp2d(i,j) = h2o_change_i(i,j) * Idt_slow + enddo ; enddo + call post_data(CS%id_net_i, tmp2d(isc:iec,jsc:jec), CS%diag) + endif + if (CS%id_net_s>0) then + !$OMP parallel do default(shared) + do j=jsc,jec ; do i=isc,iec + s2i_agg = 0.0 + do k=0,ncat + s2i_agg = s2i_agg + IST%part_size(i,j,k) * snow_to_ice_s(i,j,k) + enddo + tmp2d(i,j) = (min(h2o_change_s(i,j),0.0)*Idt_slow - sisnmass_evap(i,j) - s2i_agg) + enddo ; enddo + call post_data(CS%id_net_s, tmp2d(isc:iec,jsc:jec), CS%diag) + endif + if (CS%id_bsrc_i>0) then + !$OMP parallel do default(shared) + do j=jsc,jec ; do i=isc,iec + s2i_agg = 0.0 + do k=0,ncat + s2i_agg = s2i_agg + IST%part_size(i,j,k) * snow_to_ice_i(i,j,k) + enddo + tmp2d(i,j) = (max(h2o_change_i(i,j),0.0)*Idt_slow - s2i_agg - OSS%frazil(i,j)*ILatHtFus) + enddo ; enddo + call post_data(CS%id_bsrc_i, tmp2d(isc:iec,jsc:jec), CS%diag) + endif + if (CS%id_tsnk_i>0) then + !$OMP parallel do default(shared) + do j=jsc,jec ; do i=isc,iec + tmp2d(i,j) = (min(h2o_change_i(i,j),0.0)*Idt_slow - bsnk_i(i,j) - simass_evap(i,j)) + enddo ; enddo + call post_data(CS%id_tsnk_i, tmp2d(isc:iec,jsc:jec), CS%diag) + endif if (IOF%id_saltf>0) call post_data(IOF%id_saltf, IOF%flux_salt, CS%diag) if (CS%id_bsnk>0) call post_data(CS%id_bsnk, bsnk, CS%diag) + if (CS%id_bsnk_i>0) call post_data(CS%id_bsnk_i, bsnk_i, CS%diag) + if (CS%id_bsnk_i_cmor>0) call post_data(CS%id_bsnk_i_cmor, bsnk_i, CS%diag) + if (CS%id_bsnk_s>0) call post_data(CS%id_bsnk_s, bsnk_s, CS%diag) if (FIA%id_tmelt>0) call post_avg(FIA%id_tmelt, FIA%tmelt, IST%part_size(:,:,1:), CS%diag, G=G, & scale=Idt_slow, wtd=.true.) if (FIA%id_bmelt>0) call post_avg(FIA%id_bmelt, FIA%bmelt, IST%part_size(:,:,1:), CS%diag, G=G, & @@ -1381,10 +1527,16 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, US, IG) if (FIA%id_bheat>0) call post_data(FIA%id_bheat, OSS%bheat, CS%diag) if (CS%id_sn2ic>0) call post_avg(CS%id_sn2ic, snow_to_ice, IST%part_size(:,:,1:), CS%diag, G=G, & scale=US%RZ_T_to_kg_m2s*Idt_slow) + if (CS%id_sn2ic_i>0) call post_avg(CS%id_sn2ic_i, snow_to_ice_i, IST%part_size, CS%diag, G=G, & + scale=US%RZ_T_to_kg_m2s*Idt_slow) + if (CS%id_sn2ic_s>0) call post_avg(CS%id_sn2ic_s, snow_to_ice_s, IST%part_size, CS%diag, G=G, & + scale=US%RZ_T_to_kg_m2s*Idt_slow) if (CS%id_qflim>0) call post_data(CS%id_qflim, qflx_lim_ice, CS%diag) if (CS%id_qfres>0) call post_data(CS%id_qfres, qflx_res_ice, CS%diag) if (CS%id_net_melt>0) call post_data(CS%id_net_melt, net_melt, CS%diag) if (CS%id_CMOR_melt>0) call post_data(CS%id_CMOR_melt, net_melt, CS%diag) + if (FIA%id_evap_i>0) call post_data(FIA%id_evap_i, simass_evap, CS%diag) + if (FIA%id_evap_s>0) call post_data(FIA%id_evap_s, sisnmass_evap, CS%diag) if (coupler_type_initialized(IOF%tr_flux_ocn_top)) & call coupler_type_send_data(IOF%tr_flux_ocn_top, CS%Time) @@ -1592,9 +1744,60 @@ subroutine SIS_slow_thermo_init(Time, G, US, IG, param_file, diag, CS, tracer_fl 'frozen water local sink', 'kg/(m^2*yr)', missing_value=missing) CS%id_bsnk = register_diag_field('ice_model','BSNK',diag%axesT1, Time, & 'frozen water local bottom sink', & - 'kg/(m^2*yr)', conversion= 864e2*365.*US%RZ_T_to_kg_m2s, & + 'kg/(m^2*yr)', conversion=864e2*365.*US%RZ_T_to_kg_m2s, & + missing_value=missing) + CS%id_lsrc_i = register_diag_field('ice_model','LSRCi', diag%axesT1, Time, & + 'frozen water local source (of ice)', 'kg/(m^2*s)', conversion=US%RZ_T_to_kg_m2s, & missing_value=missing) - CS%id_sn2ic = register_diag_field('ice_model','SN2IC' ,diag%axesT1,Time, & + CS%id_lsnk_i = register_diag_field('ice_model','LSNKi',diag%axesT1, Time, & + 'frozen water local sink (of ice)', 'kg/(m^2*s)', conversion=US%RZ_T_to_kg_m2s, & + missing_value=missing) + CS%id_bsnk_i = register_diag_field('ice_model','BSNKi',diag%axesT1, Time, & + 'frozen water local bottom sink (of ice)', 'kg/(m^2*s)', conversion=US%RZ_T_to_kg_m2s, & + missing_value=missing) + CS%id_lsrc_s = register_diag_field('ice_model','LSRCs', diag%axesT1, Time, & + 'frozen water local source (of snow)', 'kg/(m^2*s)', conversion=US%RZ_T_to_kg_m2s, & + missing_value=missing) + CS%id_lsnk_s = register_diag_field('ice_model','LSNKs',diag%axesT1, Time, & + 'frozen water local sink (of snow)', 'kg/(m^2*s)', conversion=US%RZ_T_to_kg_m2s, & + missing_value=missing) + CS%id_bsnk_s = register_diag_field('ice_model','BSNKs',diag%axesT1, Time, & + 'frozen water local bottom sink (of snow)', 'kg/(m^2*s)', conversion=US%RZ_T_to_kg_m2s, & + missing_value=missing) + CS%id_lsrc_c = register_diag_field('ice_model','LSRCc', diag%axesT1, Time, & + 'frozen water area local source', 's-1', conversion=US%RZ_T_to_kg_m2s, & + missing_value=missing) + CS%id_lsnk_c = register_diag_field('ice_model','LSNKc',diag%axesT1, Time, & + 'frozen water area local sink', 's-1', conversion=US%RZ_T_to_kg_m2s, & + missing_value=missing) + + !CMOR diagnostics for thermodynamics + CS%id_bsnk_i_cmor = register_diag_field('ice_model','sidmassmeltbot',diag%axesT1, Time, & + 'Sea-Ice Mass Change Through Bottom Melting', 'kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s, & + missing_value=missing, cmor_standard_name='tendency_of_sea_ice_amount_due_to_basal_melting') + CS%id_tsnk_i = register_diag_field('ice_model','sidmassmelttop', diag%axesT1, Time, & + 'Sea-Ice Mass Change Through Surface Melting', 'kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s, & + missing_value=missing, cmor_standard_name='tendency_of_sea_ice_amount_due_to_surface_melting') + CS%id_bsrc_i = register_diag_field('ice_model','sidmassgrowthbot',diag%axesT1, Time, & + 'Sea-Ice Mass Change Through Basal Growth', 'kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s, & + missing_value=missing, cmor_standard_name='tendency_of_sea_ice_amount_due_to_congelation_ice_accumulation') + CS%id_net_i = register_diag_field('ice_model','sidmassth',diag%axesT1, Time, & + 'Sea-Ice Mass Change from Thermodynamics', 'kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s, & + missing_value=missing, cmor_standard_name='tendency_of_sea_ice_amount_due_to_sea_ice_thermodynamics') + CS%id_net_s = register_diag_field('ice_model','sisndmassmelt', diag%axesT1, Time, & + 'Snow Mass Rate of Change Through Melt', 'kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s, & + missing_value=missing, cmor_standard_name='surface_snow_melt_flux') + CS%id_net_c = register_diag_field('ice_model','sidconcth', diag%axesT1, Time, & + 'Sea-Ice Area Fraction Tendency Due to Thermodynamics', 's-1', conversion=US%RZ_T_to_kg_m2s, & + missing_value=missing, cmor_standard_name='tendency_of_sea_ice_area_fraction_due_to_thermodynamics') + CS%id_sn2ic_i = register_diag_field('ice_model','sidmassgrowthsi', diag%axesT1,Time, & + 'Sea-Ice Mass Change Through Snow-to-Ice Conversion', 'kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s, & + missing_value=missing, cmor_standard_name='tendency_of_sea_ice_amount_due_to_conversion_of_snow_to_sea_ice') + CS%id_sn2ic_s = register_diag_field('ice_model','sisndmasssi', diag%axesT1,Time, & + 'Snow Mass Rate of Change Through Snow-to-Ice Conversion', 'kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s, & + missing_value=missing, cmor_standard_name='tendency_of_surface_snow_amount_due_to_conversion_of_snow_to_sea_ice') + + CS%id_sn2ic = register_diag_field('ice_model','SN2IC', diag%axesT1,Time, & 'rate of snow to ice conversion', 'kg/(m^2*s)', missing_value=missing) CS%id_net_melt = register_diag_field('ice_model','net_melt' ,diag%axesT1, Time, & 'net mass flux from ice & snow to ocean due to melting & freezing', & diff --git a/src/SIS_transport.F90 b/src/SIS_transport.F90 index d928b34..2886bfc 100644 --- a/src/SIS_transport.F90 +++ b/src/SIS_transport.F90 @@ -71,6 +71,8 @@ module SIS_transport !>@{ Diagnostic IDs integer :: id_ix_trans = -1, id_iy_trans = -1, id_xprt = -1, id_rdgr = -1 + integer :: id_xprt_i = -1, id_xprt_s = -1, id_xprt_c = -1 + integer :: id_xprt_i_cmor = -1, id_xprt_s_cmor = -1, id_xprt_c_cmor = -1 integer :: id_rdgh=-1 ! integer :: id_rdgo=-1, id_rdgv=-1 ! These do not exist yet !!@} @@ -97,6 +99,12 @@ module SIS_transport real :: dt_sum = 0.0 !< The accumulated time since the fields were populated from an ice state type [T ~> s]. real, allocatable, dimension(:,:) :: mass0 !< The total mass of ice, snow and melt pond water !! when the fields were populated [R Z ~> kg m-2]. + real, allocatable, dimension(:,:) :: mI0 !< The total mass of ice + !! when the fields were populated [R Z ~> kg m-2]. + real, allocatable, dimension(:,:) :: mS0 !< The total mass of snow and melt pond water + !! when the fields were populated [R Z ~> kg m-2]. + real, allocatable, dimension(:,:) :: cvr0 !< The total area of sea ice + !! when the fields were populated [nondim] real, allocatable, dimension(:,:) :: uh_sum !< The accumulated zonal mass fluxes of ice, snow !! and melt pond water, summed across categories, !! since the fields were populated [R Z L2 ~> kg]. @@ -258,7 +266,10 @@ subroutine finish_ice_transport(CAS, IST, TrReg, G, US, IG, dt, CS, OSS, rdg_rat ! rdg_open, & ! formation rate of open water due to ridging [T-1 ~> s-1] ! rdg_vosh ! rate of ice mass shifted from level to ridged ice [R Z T-1 ~> kg m-2 s-1] real :: yr_dt ! Tne number of timesteps in a year [nondim]. - real, dimension(SZI_(G),SZJ_(G)) :: trans_conv ! The convergence of frozen water transport [R Z ~> kg m-2]. + real, dimension(SZI_(G),SZJ_(G)) :: trans_conv ! The convergence of frozen water transport of ice and snow [R Z ~> kg m-2]. + real, dimension(SZI_(G),SZJ_(G)) :: trans_conv_i ! The convergence of frozen water transport of ice [R Z ~> kg m-2]. + real, dimension(SZI_(G),SZJ_(G)) :: trans_conv_s ! The convergence of frozen water transport of snow [R Z ~> kg m-2]. + real, dimension(SZI_(G),SZJ_(G)) :: trans_conv_c ! The convergence of frozen water area [nondim]. real, dimension(SZI_(G),SZJ_(G)) :: ice_cover ! The summed fractional ice concentration [nondim]. type(EFP_type) :: tot_ice, tot_snow, enth_ice, enth_snow real :: I_tot_ice, I_tot_snow @@ -386,6 +397,30 @@ subroutine finish_ice_transport(CAS, IST, TrReg, G, US, IG, dt, CS, OSS, rdg_rat enddo ; enddo call post_SIS_data(CS%id_xprt, trans_conv, CS%diag) endif + if ((CS%id_xprt_i>0) .or. (CS%id_xprt_i_cmor>0)) then + call get_ice_mass(IST, G, IG, trans_conv) + do j=jsc,jec ; do i=isc,iec + trans_conv(i,j) = (trans_conv(i,j) - CAS%mI0(i,j)) * Idt + enddo ; enddo + if (CS%id_xprt_i>0) call post_SIS_data(CS%id_xprt_i, trans_conv, CS%diag) + if (CS%id_xprt_i_cmor>0) call post_SIS_data(CS%id_xprt_i_cmor, trans_conv, CS%diag) + endif + if ((CS%id_xprt_s>0) .or. (CS%id_xprt_s_cmor>0)) then + call get_snow_mass(IST, G, IG, trans_conv) + do j=jsc,jec ; do i=isc,iec + trans_conv(i,j) = (trans_conv(i,j) - CAS%mS0(i,j)) * Idt + enddo ; enddo + if (CS%id_xprt_s>0) call post_SIS_data(CS%id_xprt_s, trans_conv, CS%diag) + if (CS%id_xprt_s_cmor>0) call post_SIS_data(CS%id_xprt_s_cmor, trans_conv, CS%diag) + endif + if ((CS%id_xprt_c>0) .or. (CS%id_xprt_c_cmor>0)) then + call get_ice_area(IST, G, IG, trans_conv) + do j=jsc,jec ; do i=isc,iec + trans_conv(i,j) = (trans_conv(i,j) - CAS%cvr0(i,j)) * Idt + enddo ; enddo + if (CS%id_xprt_c>0) call post_SIS_data(CS%id_xprt_c, trans_conv, CS%diag) + if (CS%id_xprt_c_cmor>0) call post_SIS_data(CS%id_xprt_c_cmor, trans_conv, CS%diag) + endif if (CS%id_ix_trans>0) then do j=jsc,jec ; do I=isc-1,iec ; uf(I,j) = Idt * CAS%uh_sum(I,j) ; enddo ; enddo call post_SIS_data(CS%id_ix_trans, uf, CS%diag) @@ -473,6 +508,9 @@ subroutine ice_state_to_cell_ave_state(IST, G, US, IG, CS, CAS) ! Handle diagnostics CAS%dt_sum = 0.0 if (allocated(CAS%mass0)) call get_cell_mass(IST, G, IG, CAS%mass0) + if (allocated(CAS%mI0)) call get_ice_mass(IST, G, IG, CAS%mI0) + if (allocated(CAS%mS0)) call get_snow_mass(IST, G, IG, CAS%mS0) + if (allocated(CAS%cvr0)) call get_ice_area(IST, G, IG, CAS%cvr0) if (allocated(CAS%uh_sum)) CAS%uh_sum(:,:) = 0.0 if (allocated(CAS%vh_sum)) CAS%vh_sum(:,:) = 0.0 @@ -1073,6 +1111,65 @@ subroutine get_cell_mass(IST, G, IG, cell_mass, scale) end subroutine get_cell_mass +!> get_ice_mass determines the integrated mass of ice in each cell +subroutine get_ice_mass(IST, G, IG, cell_mass, 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(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: cell_mass !< The total amount of ice [R Z ~> kg m-2]. + real, optional, intent(in) :: scale !< A scaling factor from H to the desired units. + + real :: H_to_units ! A conversion factor from H to the desired output units. + integer :: i, j, k, isc, iec, jsc, jec + isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec + + H_to_units = 1.0 ; if (present(scale)) H_to_units = 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_ice(i,j,k) + enddo ; enddo ; enddo + +end subroutine get_ice_mass + +!> get_snow_mass determines the integrated mass of snow and ponds in each cell +subroutine get_snow_mass(IST, G, IG, cell_mass, 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(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: cell_mass !< The total amount of snow [R Z ~> kg m-2]. + real, optional, intent(in) :: scale !< A scaling factor from H to the desired units. + + real :: H_to_units ! A conversion factor from H to the desired output units. + integer :: i, j, k, isc, iec, jsc, jec + isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec + + H_to_units = 1.0 ; if (present(scale)) H_to_units = 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)) + enddo ; enddo ; enddo + +end subroutine get_snow_mass + +!> get_ice_area determines the integrated area of ice in each grid cell +subroutine get_ice_area(IST, G, IG, cell_area) + 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(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: cell_area !< The fractional cover of ice [nondim]. + + integer :: i, j, k, isc, iec, jsc, jec + isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec + + cell_area(:,:) = 0.0 + do k=1,IG%CatIce ; do j=jsc,jec ; do i=isc,iec + cell_area(i,j) = cell_area(i,j) + IST%part_size(i,j,k) + enddo ; enddo ; enddo + +end subroutine get_ice_area + subroutine cell_mass_from_CAS(CAS, G, IG, mca, scale) type(cell_average_state_type), intent(in) :: CAS !< A structure with ocean-cell averaged masses by !! category and phase of water. @@ -1260,6 +1357,27 @@ subroutine SIS_transport_init(Time, G, IG, US, param_file, diag, CS, continuity_ CS%id_xprt = register_diag_field('ice_model', 'XPRT', diag%axesT1, Time, & 'frozen water transport convergence', 'kg/(m^2*yr)', conversion=US%RZ_to_kg_m2, & missing_value=missing) + CS%id_xprt_i = register_diag_field('ice_model', 'XPRTi', diag%axesT1, Time, & + 'frozen water transport convergence (of ice)', 'kg/(m^2*s)', conversion=US%RZ_to_kg_m2, & + missing_value=missing) + CS%id_xprt_s = register_diag_field('ice_model', 'XPRTs', diag%axesT1, Time, & + 'frozen water transport convergence (of snow)', 'kg/(m^2*s)', conversion=US%RZ_to_kg_m2, & + missing_value=missing) + CS%id_xprt_c = register_diag_field('ice_model', 'XPRTc', diag%axesT1, Time, & + 'frozen water area transport convergence', 's-1', conversion=US%RZ_to_kg_m2, & + missing_value=missing) + + !CMOR diagnostics for dynamics + CS%id_xprt_i_cmor = register_diag_field('ice_model', 'sidmassdyn', diag%axesT1, Time, & + 'Sea-Ice Mass Change from Dynamics', 'kg m-2 s-1', conversion=US%RZ_to_kg_m2, & + missing_value=missing, cmor_standard_name='tendency_of_sea_ice_amount_due_to_dynamics') + CS%id_xprt_s_cmor = register_diag_field('ice_model', 'sisndmassdyn', diag%axesT1, Time, & + 'Snow Mass Rate of Change Through Advection by Sea-Ice Dynamics', 'kg m-2 s-1', conversion=US%RZ_to_kg_m2, & + missing_value=missing, cmor_standard_name='tendency_of_surface_snow_amount_due_to_sea_ice_dynamics') + CS%id_xprt_c_cmor = register_diag_field('ice_model', 'sidconcdyn', diag%axesT1, Time, & + 'Sea-Ice Area Fraction Tendency Due to Dynamics', 's-1', conversion=US%RZ_to_kg_m2, & + missing_value=missing, cmor_standard_name='tendency_of_sea_ice_area_fraction_due_to_dynamics') + CS%id_rdgr = register_diag_field('ice_model', 'RDG_RATE', diag%axesT1, Time, & 'ice ridging rate', '1/sec', conversion=US%s_to_T, missing_value=missing) CS%id_rdgh = register_diag_field('ice_model', 'RDG_HEIGHT', diag%axesTc, Time, & @@ -1293,6 +1411,12 @@ subroutine alloc_cell_average_state_type(CAS, HI, IG, CS) if (present(CS)) then if (CS%id_xprt>0) & call safe_alloc(CAS%mass0, isd, ied, jsd, jed) + if ((CS%id_xprt_i>0) .or. (CS%id_xprt_i_cmor>0)) & + call safe_alloc(CAS%mI0, isd, ied, jsd, jed) + if ((CS%id_xprt_s>0) .or. (CS%id_xprt_s_cmor>0)) & + call safe_alloc(CAS%mS0, isd, ied, jsd, jed) + if ((CS%id_xprt_c>0) .or. (CS%id_xprt_c_cmor>0)) & + call safe_alloc(CAS%cvr0, isd, ied, jsd, jed) if (CS%id_ix_trans>0) & call safe_alloc(CAS%uh_sum, HI%IsdB, HI%IedB, jsd, jed) if (CS%id_iy_trans>0) & @@ -1307,6 +1431,9 @@ subroutine dealloc_cell_average_state_type(CAS) if (.not.associated(CAS)) return deallocate(CAS%m_ice, CAS%m_snow, CAS%m_pond, CAS%mH_ice) if (allocated(CAS%mass0)) deallocate(CAS%mass0) + if (allocated(CAS%mI0)) deallocate(CAS%mI0) + if (allocated(CAS%mS0)) deallocate(CAS%mS0) + if (allocated(CAS%cvr0)) deallocate(CAS%cvr0) if (allocated(CAS%uh_sum)) deallocate(CAS%uh_sum) if (allocated(CAS%vh_sum)) deallocate(CAS%vh_sum) deallocate(CAS) diff --git a/src/SIS_types.F90 b/src/SIS_types.F90 index b1cb332..0a92a85 100644 --- a/src/SIS_types.F90 +++ b/src/SIS_types.F90 @@ -144,7 +144,7 @@ module SIS_types logical :: Cgrid_dyn !< If true use a C-grid discretization of the sea-ice dynamics. !>@{ diagnostic IDs for ocean surface properties - integer :: id_sst=-1, id_sss=-1, id_ssh=-1, id_uo=-1, id_vo=-1, id_frazil=-1 + integer :: id_sst=-1, id_sss=-1, id_ssh=-1, id_uo=-1, id_vo=-1, id_frazil=-1, id_frazilmass=-1 !!@} end type ocean_sfc_state_type @@ -257,7 +257,7 @@ module SIS_types !SLOW ONLY !!@{ Diagnostic IDs integer :: id_sh=-1, id_lh=-1, id_sw=-1, id_slp=-1 - integer :: id_lw=-1, id_snofl=-1, id_rain=-1, id_evap=-1 + integer :: id_lw=-1, id_snofl=-1, id_rain=-1, id_evap=-1, id_evap_i=-1, id_evap_s=-1 integer :: id_sw_vis_dir=-1, id_sw_vis_dif=-1, id_sw_nir_dir=-1, id_sw_nir_dif=-1 integer :: id_sw_vis=-1, id_sw_dir=-1, id_sw_dif=-1, id_sw_dn=-1, id_albedo=-1 integer :: id_runoff=-1, id_calving=-1, id_runoff_hflx=-1, id_calving_hflx=-1