!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS is distributed in the hope that it will be useful, !! !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! !! GNU General Public License for more details. !! !! !! !! You should have received a copy of the GNU General Public License !! !! along with FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! A SEA ICE MODEL for coupling through the GFDL exchange grid; this module ! ! manages fluxes, diagnostics, and ice timesteps; sea ice dynamics and ! ! thermodynamics are performed in ice_[dyn|thm].f90 - Mike Winton (Michael.Winton)! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! module ice_model_mod use mpp_mod, only: mpp_clock_begin, mpp_clock_end use mpp_domains_mod, only: mpp_update_domains, BGRID_NE, CGRID_NE use fms_mod, only: error_mesg use diag_manager_mod, only: send_data use time_manager_mod, only: time_type, operator(+), get_date, get_time use time_manager_mod, only: operator(-), set_date use astronomy_mod, only: universal_time, orbital_time, diurnal_solar, daily_mean_solar use coupler_types_mod,only: coupler_2d_bc_type, coupler_3d_bc_type use constants_mod, only: hlv, hlf, Tfreeze, grav, STEFAN use ocean_albedo_mod, only: compute_ocean_albedo ! ice sets ocean surface use ocean_rough_mod, only: compute_ocean_roughness ! properties over water use ice_type_mod, only: ice_data_type, ice_model_init, ice_model_end, hlim,& mom_rough_ice, heat_rough_ice, atmos_winds, kmelt, & slab_ice, spec_ice, ice_bulk_salin, id_cn, id_hi, & id_hs, id_t1, id_t2, id_ts, id_sh, id_lh, id_sw, & id_swdn, id_lw, id_lwdn, id_snofl, id_rain, & id_runoff_hflx, id_calving_hflx, & id_runoff, id_calving, id_evap, id_saltf, id_tmelt,& id_mi, id_bmelt, id_bheat, id_frazil, id_alb, & id_xprt, id_lsrc, id_lsnk, id_bsnk, id_strna, & id_sigi, id_sigii, id_stren, id_ui, id_vi, id_fax, & id_fay, id_fix, id_fiy, id_fcx, id_fcy, id_fwx, & id_fwy, id_sn2ic, id_ext, id_slp, id_sst, id_sss, & id_ssh, id_uo, id_vo, id_e2m, id_qflim, id_qfres, & id_ix_trans, id_iy_trans, & do_ice_restore, do_ice_limit, max_ice_limit, & ice_restore_timescale, do_init, h2o, heat, & conservation_check, slp2ocean, iceClock, verbose, & iceClock1, iceClock2, iceClock3, & iceClock4, iceClock5, iceClock6, & iceClock7, iceClock8, iceClock9, & iceClocka, iceClockb, iceClockc, & id_sw_vis, id_sw_dir, id_sw_dif, id_sw_vis_dir, & id_sw_vis_dif, id_sw_nir_dir, id_sw_nir_dif, & cm2_bugs, ice_stock_pe, do_icebergs, ice_model_restart, & add_diurnal_sw, id_mib, ice_data_type_chksum, & id_ustar, id_vstar, channel_viscosity, smag_ocn, & ssh_gravity, chan_cfl_limit, id_vocean, id_uocean, & id_vchan, id_uchan, id_wnd use ice_type_mod, only: do_sun_angle_for_alb, & id_alb_vis_dir, id_alb_vis_dif, & id_alb_nir_dir, id_alb_nir_dif use ice_grid_mod, only: uv_to_t, t_to_uv, vel_t_to_uv, cut_check, tripolar_grid use ice_grid_mod, only: Domain, isc, iec, jsc, jec, isd, ied, jsd, jed, im, jm, km use ice_grid_mod, only: ice_advect, ice_avg, all_avg, ice_line, slab_ice_advect use ice_grid_mod, only: geo_lon, geo_lat, cell_area, sin_rot, cos_rot, latitude use ice_spec_mod, only: get_sea_surface use ice_grid_mod, only: dte, dtn, dxv, dyv, dxt, dyt, dt_adv, wett use ice_grid_mod, only: reproduce_siena_201303 ! ! the following two modules are the work horses of the sea ice model ! use ice_thm_mod, only: ice_optics, ice_thm_param, ice3lay_temp, ice3lay_resize use ice_thm_mod, only: thm_pack, thm_unpack, DI, DS, MU_TS, TFI, e_to_melt use ice_dyn_mod, only: ice_dynamics, ice_dyn_param, strain_angle, ice_strength, sigI, sigII use ice_bergs, only: icebergs_run, icebergs_incr_mass implicit none private public :: ice_data_type, ocean_ice_boundary_type, atmos_ice_boundary_type, land_ice_boundary_type public :: ice_model_init, ice_model_end, ice_bottom_to_ice_top, update_ice_model_fast, ice_stock_pe, cell_area public :: update_ice_model_slow, update_ice_model_slow_up, update_ice_model_slow_dn ! for new coupler public :: ice_model_restart ! for intermediate restart public :: ocn_ice_bnd_type_chksum, atm_ice_bnd_type_chksum, & lnd_ice_bnd_type_chksum, ice_data_type_chksum interface update_ice_model_fast ! overload to support old interface module procedure update_ice_model_fast_new ! module procedure update_ice_model_fast_old !Balaji: no need to expose old interface end interface ! ! the following three types are for data exchange with the new coupler ! they are defined here but declared in coupler_main and allocated in flux_init ! type :: ocean_ice_boundary_type real, dimension(:,:), pointer :: u =>NULL() real, dimension(:,:), pointer :: v =>NULL() real, dimension(:,:), pointer :: t =>NULL() real, dimension(:,:), pointer :: s =>NULL() real, dimension(:,:), pointer :: frazil =>NULL() real, dimension(:,:), pointer :: sea_level =>NULL() real, dimension(:,:,:), pointer :: data =>NULL() ! collective field for "named" fields above integer :: xtype ! REGRID, REDIST or DIRECT used by coupler type(coupler_2d_bc_type) :: fields ! array of fields used for additional tracers end type type :: atmos_ice_boundary_type real, dimension(:,:,:), pointer :: u_flux =>NULL() real, dimension(:,:,:), pointer :: v_flux =>NULL() real, dimension(:,:,:), pointer :: u_star =>NULL() real, dimension(:,:,:), pointer :: t_flux =>NULL() real, dimension(:,:,:), pointer :: q_flux =>NULL() real, dimension(:,:,:), pointer :: lw_flux =>NULL() real, dimension(:,:,:), pointer :: sw_flux_vis_dir =>NULL() real, dimension(:,:,:), pointer :: sw_flux_vis_dif =>NULL() real, dimension(:,:,:), pointer :: sw_flux_nir_dir =>NULL() real, dimension(:,:,:), pointer :: sw_flux_nir_dif =>NULL() real, dimension(:,:,:), pointer :: lprec =>NULL() real, dimension(:,:,:), pointer :: fprec =>NULL() real, dimension(:,:,:), pointer :: dhdt =>NULL() real, dimension(:,:,:), pointer :: dedt =>NULL() real, dimension(:,:,:), pointer :: drdt =>NULL() real, dimension(:,:,:), pointer :: coszen =>NULL() real, dimension(:,:,:), pointer :: coszen_radwt =>NULL() real, dimension(:,:,:), pointer :: p =>NULL() real, dimension(:,:,:), pointer :: data =>NULL() integer :: xtype type(coupler_3d_bc_type) :: fluxes ! array of fluxes used for additional tracers end type type :: land_ice_boundary_type real, dimension(:,:), pointer :: runoff =>NULL() real, dimension(:,:), pointer :: calving =>NULL() real, dimension(:,:), pointer :: runoff_hflx =>NULL() real, dimension(:,:), pointer :: calving_hflx =>NULL() real, dimension(:,:,:), pointer :: data =>NULL() ! collective field for "named" fields above integer :: xtype ! REGRID, REDIST or DIRECT used by coupler end type real, parameter :: T_sw_freeze = Tfreeze-1.8 ! seawater freezing temperature (K) logical :: first_time = .true. ! first time ice_bottom_to_ice_top contains !####################################################################### ! ! new coupler interface to provide ocean surface data to atmosphere ! subroutine update_ice_model_slow_up ( Ocean_boundary, Ice ) type(ocean_ice_boundary_type), intent(inout) :: Ocean_boundary type (ice_data_type), intent(inout) :: Ice call mpp_clock_begin(iceClock) call mpp_clock_begin(iceClock1) call ice_bottom_to_ice_top (Ice, Ocean_boundary%t, Ocean_boundary%u, Ocean_boundary%v, & Ocean_boundary%frazil, Ocean_boundary, Ocean_boundary%s, Ocean_boundary%sea_level ) call mpp_clock_end(iceClock1) call mpp_clock_end(iceClock) end subroutine update_ice_model_slow_up ! ! new coupler interface to do slow ice processes: dynamics, transport, mass ! subroutine update_ice_model_slow_dn ( Atmos_boundary, Land_boundary, Ice ) type(atmos_ice_boundary_type), intent(inout) :: Atmos_boundary type(land_ice_boundary_type), intent(inout) :: Land_boundary type (ice_data_type), intent(inout) :: Ice call mpp_clock_begin(iceClock) call mpp_clock_begin(iceClock2) call update_ice_model_slow (Ice, Land_boundary%runoff, Land_boundary%calving, Land_boundary%runoff_hflx, Land_boundary%calving_hflx, Atmos_boundary%p ) call mpp_clock_end(iceClock2) call mpp_clock_end(iceClock) end subroutine update_ice_model_slow_dn subroutine update_ice_model_fast_new ( Atmos_boundary, Ice ) type(atmos_ice_boundary_type), intent(inout) :: Atmos_boundary type (ice_data_type), intent(inout) :: Ice call mpp_clock_begin(iceClock) call mpp_clock_begin(iceClock3) call update_ice_model_fast_old (Ice, Atmos_boundary%fluxes, & Atmos_boundary%u_flux, & Atmos_boundary%v_flux, & Atmos_boundary%u_star, & Atmos_boundary%sw_flux_nir_dir, & Atmos_boundary%sw_flux_nir_dif, & Atmos_boundary%sw_flux_vis_dir, & Atmos_boundary%sw_flux_vis_dif, & Atmos_boundary%lw_flux, & Atmos_boundary%t_flux, & Atmos_boundary%q_flux, & Atmos_boundary%dhdt, & Atmos_boundary%dedt, & Atmos_boundary%drdt, & Atmos_boundary%lprec, & Atmos_boundary%fprec, & Atmos_boundary%coszen, & Atmos_boundary%coszen_radwt, & Atmos_boundary%p ) call mpp_clock_end(iceClock3) call mpp_clock_end(iceClock) end subroutine update_ice_model_fast_new !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! zero_top_quantities - zero fluxes to begin summing in ice fast physics ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! subroutine zero_top_quantities ( Ice ) type (ice_data_type), intent(inout) :: Ice integer :: i, j, k, n, m Ice % avg_kount = 0 do k = 1, km do j = jsd, jed do i = isd, ied Ice % flux_u_top(i,j,k) = 0.0 Ice % flux_v_top(i,j,k) = 0.0 enddo enddo enddo do k = 1, km do j = jsc, jec do i = isc, iec Ice % flux_t_top(i,j,k) = 0.0 Ice % flux_q_top(i,j,k) = 0.0 Ice % flux_lw_top(i,j,k) = 0.0 Ice % flux_lh_top(i,j,k) = 0.0 Ice % flux_sw_nir_dir_top(i,j,k) = 0.0 Ice % flux_sw_nir_dif_top(i,j,k) = 0.0 Ice % flux_sw_vis_dir_top(i,j,k) = 0.0 Ice % flux_sw_vis_dif_top(i,j,k) = 0.0 Ice % lprec_top(i,j,k) = 0.0 Ice % fprec_top(i,j,k) = 0.0 do n = 1, Ice%ocean_fluxes_top%num_bcs !{ do m = 1, Ice%ocean_fluxes_top%bc(n)%num_fields !{ Ice%ocean_fluxes_top%bc(n)%field(m)%values(i,j,k) = 0.0 enddo !} m enddo !} n enddo enddo enddo do j = jsc, jec do i = isc, iec Ice % lwdn(i,j) = 0.0 Ice % swdn(i,j) = 0.0 enddo enddo return end subroutine zero_top_quantities !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! sum_top_quantities - sum fluxes for later use by ice/ocean slow physics ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! subroutine sum_top_quantities ( Ice, Atmos_boundary_fluxes, flux_u, flux_v, flux_t, flux_q, & flux_sw_nir_dir, flux_sw_nir_dif, flux_sw_vis_dir, flux_sw_vis_dif,& flux_lw, lprec, fprec, flux_lh ) type (ice_data_type), intent(inout) :: Ice type(coupler_3d_bc_type), intent(inout) :: Atmos_boundary_fluxes real, intent(in), dimension(isc:iec,jsc:jec,km) :: flux_u, flux_v, flux_t, flux_q real, intent(in), dimension(isc:iec,jsc:jec,km) :: flux_lw, lprec, fprec, flux_lh real, intent(in), dimension(isc:iec,jsc:jec,km) :: flux_sw_nir_dir real, intent(in), dimension(isc:iec,jsc:jec,km) :: flux_sw_nir_dif real, intent(in), dimension(isc:iec,jsc:jec,km) :: flux_sw_vis_dir real, intent(in), dimension(isc:iec,jsc:jec,km) :: flux_sw_vis_dif real,dimension(isc:iec,jsc:jec) :: tmp integer :: i, j, k, m, n if (Ice % avg_kount == 0) call zero_top_quantities (Ice) do k = 1, km do j = jsc, jec do i = isc, iec Ice % flux_u_top(i,j,k) = Ice % flux_u_top(i,j,k) + flux_u(i,j,k) Ice % flux_v_top(i,j,k) = Ice % flux_v_top(i,j,k) + flux_v(i,j,k) Ice % flux_t_top(i,j,k) = Ice % flux_t_top(i,j,k) + flux_t(i,j,k) Ice % flux_q_top(i,j,k) = Ice % flux_q_top(i,j,k) + flux_q(i,j,k) Ice % flux_sw_nir_dir_top(i,j,k) = Ice % flux_sw_nir_dir_top(i,j,k) + flux_sw_nir_dir(i,j,k) Ice % flux_sw_nir_dif_top(i,j,k) = Ice % flux_sw_nir_dif_top(i,j,k) + flux_sw_nir_dif(i,j,k) Ice % flux_sw_vis_dir_top(i,j,k) = Ice % flux_sw_vis_dir_top(i,j,k) + flux_sw_vis_dir(i,j,k) Ice % flux_sw_vis_dif_top(i,j,k) = Ice % flux_sw_vis_dif_top(i,j,k) + flux_sw_vis_dif(i,j,k) Ice % flux_lw_top(i,j,k) = Ice % flux_lw_top(i,j,k) + flux_lw(i,j,k) Ice % lprec_top(i,j,k) = Ice % lprec_top(i,j,k) + lprec(i,j,k) Ice % fprec_top(i,j,k) = Ice % fprec_top(i,j,k) + fprec(i,j,k) Ice % flux_lh_top(i,j,k) = Ice % flux_lh_top(i,j,k) + flux_lh(i,j,k) enddo enddo enddo do n = 1, Ice%ocean_fluxes_top%num_bcs !{ do m = 1, Ice%ocean_fluxes_top%bc(n)%num_fields !{ do k = 1, km do j = jsc, jec do i = isc, iec Ice%ocean_fluxes_top%bc(n)%field(m)%values(i,j,k) = Ice%ocean_fluxes_top%bc(n)%field(m)%values(i,j,k) + & Atmos_boundary_fluxes%bc(n)%field(m)%values(i,j,k) enddo !} m enddo !} n enddo enddo enddo if (id_lwdn > 0) then tmp = all_avg(flux_lw+STEFAN*Ice%t_surf**4, Ice%part_size(isc:iec,jsc:jec,:)) do j = jsc, jec do i = isc, iec if(Ice%mask(i,j) ) Ice%lwdn(i,j) = Ice%lwdn(i,j) + tmp(i,j) enddo enddo endif if (id_swdn > 0) then tmp = all_avg(flux_sw_vis_dir/(1-Ice%albedo_vis_dir)+flux_sw_vis_dif/(1-Ice%albedo_vis_dif)+ & flux_sw_nir_dir/(1-Ice%albedo_nir_dir)+flux_sw_nir_dif/(1-Ice%albedo_nir_dif), & Ice%part_size(isc:iec,jsc:jec,:) ) do j = jsc, jec do i = isc, iec if(Ice%mask(i,j) ) Ice%swdn(i,j) = Ice%swdn(i,j) + tmp(i,j) enddo enddo endif Ice % avg_kount = Ice % avg_kount + 1 end subroutine sum_top_quantities !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! avg_top_quantities - time average fluxes for ice and ocean slow physics ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! subroutine avg_top_quantities ( Ice ) type (ice_data_type), intent(inout) :: Ice real :: u, v, divid integer :: i, j, k, m, n logical :: sent ! ! compute average fluxes ! if (Ice % avg_kount == 0) call error_mesg ('avg_top_quantities', & 'no ocean model fluxes have been averaged', 3) divid = 1./float(Ice%avg_kount) do k=1,km do j = jsc, jec do i = isc, iec u = Ice % flux_u_top(i,j,k) * divid v = Ice % flux_v_top(i,j,k) * divid Ice % flux_u_top(i,j,k) = u*cos_rot(i,j)-v*sin_rot(i,j) ! rotate stress from lat/lon Ice % flux_v_top(i,j,k) = v*cos_rot(i,j)+u*sin_rot(i,j) ! to ocean coordinates end do end do end do if (atmos_winds) then ! put wind stress on u,v points and change sign to +down call mpp_update_domains(Ice % flux_u_top, Ice % flux_v_top, Domain ) do k=1,km call vel_t_to_uv( -Ice%flux_u_top(:,:,k),-Ice%flux_v_top(:,:,k), & Ice%flux_u_top_bgrid(isc:iec,jsc:jec,k), Ice%flux_v_top_bgrid(isc:iec,jsc:jec,k) ) end do else if(reproduce_siena_201303) then Ice%flux_u_top_bgrid(isc:iec,jsc:jec,:) = Ice%flux_u_top(isc:iec,jsc:jec,:) Ice%flux_v_top_bgrid(isc:iec,jsc:jec,:) = Ice%flux_v_top(isc:iec,jsc:jec,:) else call mpp_update_domains(Ice % flux_u_top, Ice % flux_v_top, Domain ) do k=1,km call vel_t_to_uv( Ice%flux_u_top(:,:,k),Ice%flux_v_top(:,:,k), & Ice%flux_u_top_bgrid(isc:iec,jsc:jec,k), Ice%flux_v_top_bgrid(isc:iec,jsc:jec,k) ) end do endif endif do k = 1, km do j = jsc, jec do i = isc, iec Ice % flux_t_top(i,j,k) = Ice % flux_t_top(i,j,k) * divid Ice % flux_q_top(i,j,k) = Ice % flux_q_top(i,j,k) * divid Ice % flux_sw_nir_dir_top(i,j,k) = Ice % flux_sw_nir_dir_top(i,j,k) * divid Ice % flux_sw_nir_dif_top(i,j,k) = Ice % flux_sw_nir_dif_top(i,j,k) * divid Ice % flux_sw_vis_dir_top(i,j,k) = Ice % flux_sw_vis_dir_top(i,j,k) * divid Ice % flux_sw_vis_dif_top(i,j,k) = Ice % flux_sw_vis_dif_top(i,j,k) * divid Ice % flux_lw_top(i,j,k) = Ice % flux_lw_top(i,j,k) * divid Ice % fprec_top(i,j,k) = Ice % fprec_top(i,j,k) * divid Ice % lprec_top(i,j,k) = Ice % lprec_top(i,j,k) * divid Ice % flux_lh_top(i,j,k) = Ice % flux_lh_top(i,j,k) * divid do n = 1, Ice%ocean_fluxes_top%num_bcs !{ do m = 1, Ice%ocean_fluxes_top%bc(n)%num_fields !{ Ice%ocean_fluxes_top%bc(n)%field(m)%values(i,j,k) = Ice%ocean_fluxes_top%bc(n)%field(m)%values(i,j,k) * divid enddo !} m enddo !} n enddo enddo enddo do j = jsc, jec do i = isc, iec Ice % lwdn(i,j) = Ice % lwdn(i,j)* divid Ice % swdn(i,j) = Ice % swdn(i,j)* divid enddo enddo ! ! Flux diagnostics ! if (id_sh>0) sent = send_data(id_sh, all_avg(Ice%flux_t_top,Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) if (id_lh>0) sent = send_data(id_lh, all_avg(Ice%flux_lh_top,Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) if (id_evap>0) sent = send_data(id_evap, all_avg(Ice%flux_q_top,Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) if (id_sw>0) sent = send_data(id_sw, all_avg(Ice%flux_sw_vis_dir_top+Ice%flux_sw_vis_dif_top+ & Ice%flux_sw_nir_dir_top+Ice%flux_sw_nir_dif_top,Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) if (id_lw>0) sent = send_data(id_lw, all_avg(Ice%flux_lw_top,Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) if (id_snofl>0) sent = send_data(id_snofl, all_avg(Ice%fprec_top,Ice%part_size(isc:iec,jsc:jec,:)),& Ice%Time, mask=Ice%mask) if (id_rain>0) sent = send_data(id_rain, all_avg(Ice%lprec_top,Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) if (id_lwdn>0) sent = send_data(id_lwdn, Ice%lwdn, Ice%Time, mask=Ice%mask) if (id_swdn>0) sent = send_data(id_swdn, Ice%swdn, Ice%Time, mask=Ice%mask) if (id_sw_vis>0) sent = send_data(id_sw_vis, all_avg(Ice%flux_sw_vis_dif_top+Ice%flux_sw_vis_dir_top, & Ice%part_size(isc:iec,jsc:jec,:)), Ice%Time, mask=Ice%mask) if (id_sw_nir_dir>0) sent = send_data(id_sw_nir_dir, all_avg(Ice%flux_sw_nir_dir_top,Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) if (id_sw_nir_dif>0) sent = send_data(id_sw_nir_dif, all_avg(Ice%flux_sw_nir_dif_top,Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) if (id_sw_vis_dir>0) sent = send_data(id_sw_vis_dir, all_avg(Ice%flux_sw_vis_dir_top,Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) if (id_sw_vis_dif>0) sent = send_data(id_sw_vis_dif, all_avg(Ice%flux_sw_vis_dif_top,Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) ! ! set count to zero and fluxes will be zeroed before the next sum ! Ice % avg_kount = 0 end subroutine avg_top_quantities subroutine ice_top_to_ice_bottom (Ice, part_size, part_size_uv) type (ice_data_type), intent(inout) :: Ice real, dimension (:,:,:), intent(in) :: part_size, part_size_uv integer :: m, n Ice % flux_u = all_avg( Ice % flux_u_top_bgrid(isc:iec,jsc:jec,:) , part_size_uv ) Ice % flux_v = all_avg( Ice % flux_v_top_bgrid(isc:iec,jsc:jec,:) , part_size_uv ) Ice % flux_t = all_avg( Ice % flux_t_top , part_size ) Ice % flux_q = all_avg( Ice % flux_q_top , part_size ) Ice % flux_sw_nir_dir = all_avg( Ice % flux_sw_nir_dir_top, part_size ) Ice % flux_sw_nir_dif = all_avg( Ice % flux_sw_nir_dif_top, part_size ) Ice % flux_sw_vis_dir = all_avg( Ice % flux_sw_vis_dir_top, part_size ) Ice % flux_sw_vis_dif = all_avg( Ice % flux_sw_vis_dif_top, part_size ) Ice % flux_lw = all_avg( Ice % flux_lw_top, part_size ) Ice % fprec = all_avg( Ice % fprec_top , part_size ) Ice % lprec = all_avg( Ice % lprec_top , part_size ) Ice % flux_lh = all_avg( Ice % flux_lh_top, part_size ) do n = 1, Ice%ocean_fluxes%num_bcs !{ do m = 1, Ice%ocean_fluxes%bc(n)%num_fields !{ Ice%ocean_fluxes%bc(n)%field(m)%values = & all_avg(Ice%ocean_fluxes_top%bc(n)%field(m)%values, part_size) enddo !} m enddo !} n end subroutine ice_top_to_ice_bottom !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! ice_bottom_to_ice_top - prepare surface state for atmosphere fast physics ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! subroutine ice_bottom_to_ice_top (Ice, t_surf_ice_bot, u_surf_ice_bot, v_surf_ice_bot, & frazil_ice_bot, Ocean_ice_boundary, s_surf_ice_bot, sea_lev_ice_bot ) type (ice_data_type), intent(inout) :: Ice real, dimension(isc:iec,jsc:jec), intent(in) :: t_surf_ice_bot, u_surf_ice_bot real, dimension(isc:iec,jsc:jec), intent(in) :: v_surf_ice_bot, frazil_ice_bot type(ocean_ice_boundary_type), intent(inout) :: Ocean_ice_boundary real, dimension(isc:iec,jsc:jec), intent(in), optional :: s_surf_ice_bot, sea_lev_ice_bot real, dimension(isc:iec,jsc:jec) :: sst, tmp real :: u, v integer :: i, j, k, m, n logical :: sent real, parameter :: LI = hlf ! ! pass ocean state through ice on first partition ! if (.not. spec_ice) & ! otherwise, already set by update_ice_model_slow Ice%t_surf(:,:,1) = t_surf_ice_bot if (do_init) then call get_sea_surface(Ice%Time, Ice%t_surf(:,:,1), Ice%part_size(isc:iec,jsc:jec,1:2), & Ice%h_ice (isc:iec,jsc:jec,2 ) ) call mpp_update_domains(Ice%part_size(:,:,1:2), Domain ) ! these two updates cannot be combined call mpp_update_domains(Ice%h_ice(:,:,2), Domain ) ! these two updates cannot be combined Ice%part_size_uv(:,:,1) = 1.0 do k=2,km call t_to_uv(Ice%part_size(:,:,k),Ice%part_size_uv(:,:,k)) enddo do k=2,km do j=jsc,jec do i=isc,iec Ice%part_size_uv (i,j,1) = Ice%part_size_uv(i,j,1)-Ice%part_size_uv (i,j,k) enddo end do end do do_init = .false. end if if (first_time .and. conservation_check) then h2o = 0.0 h2o(1) = sum(cell_area*all_avg(DI*Ice%h_ice(isc:iec,jsc:jec,:)+ & DS*Ice%h_snow(isc:iec,jsc:jec,:),Ice%part_size(isc:iec,jsc:jec,:))) heat = 0.0 do k=2,km do j=jsc,jec do i=isc,iec if ((Ice%part_size(i,j,k)>0.0.and.Ice%h_ice(i,j,k)>0.0)) then if (slab_ice) then heat(1) = heat(1) - cell_area(i,j) * Ice%part_size(i,j,k) * Ice%h_ice(i,j,2)*DI*LI else heat(1) = heat(1) - cell_area(i,j) * Ice%part_size(i,j,k) * & e_to_melt(Ice%h_snow(i,j,k), Ice%h_ice(i,j,k)/2, & Ice%t_ice1(i,j,k), Ice%h_ice(i,j,k)/2, Ice%t_ice2(i,j,k)) end if end if end do end do end do first_time = .false.; end if sst = Ice%t_surf(:,:,1) if (present(s_surf_ice_bot)) then do j = jsc, jec do i = isc, iec Ice%s_surf(i,j) = s_surf_ice_bot(i,j) enddo enddo else do j = jsc, jec do i = isc, iec Ice%s_surf(i,j) = 0.0 enddo enddo end if do j = jsc,jec do i = isc,iec Ice%frazil(i,j) = frazil_ice_bot(i,j) enddo enddo ! transfer the ocean state for extra tracer fluxes ! do n = 1, Ocean_ice_boundary%fields%num_bcs !{ do m = 1, Ocean_ice_boundary%fields%bc(n)%num_fields !{ Ice%ocean_fields%bc(n)%field(m)%values(:,:,1) = Ocean_ice_boundary%fields%bc(n)%field(m)%values enddo !} m enddo !} n if (present(sea_lev_ice_bot)) then do j = jsc, jec do i = isc, iec Ice%sea_lev(i,j) = sea_lev_ice_bot(i,j) enddo enddo else do j = jsc, jec do i = isc, iec Ice%sea_lev(i,j) = 0.0 enddo enddo end if call mpp_update_domains(Ice%sea_lev, Domain) do k = 2, km do j = jsc,jec do i = isc,iec Ice % tmelt(i,j,k) = 0.0 Ice % bmelt(i,j,k) = 0.0 enddo enddo enddo tmp = ice_avg(Ice%h_ice(isc:iec,jsc:jec,:),Ice%part_size(isc:iec,jsc:jec,:) ) do j = jsc, jec do i = isc, iec if( tmp(i,j) > 0.0) then Ice%bheat(i,j) = kmelt*(sst(i,j)-Tfreeze+MU_TS*Ice%s_surf(i,j)) else Ice%bheat(i,j) = 0.0 endif enddo enddo do j = jsc, jec do i = isc, iec Ice%ice_mask(i,j,1) = .false. enddo enddo do k = 2, km do j = jsc, jec do i = isc, iec if(Ice%h_ice(i,j,k) > 0.0) then Ice%ice_mask(i,j,k) = .true. else Ice%ice_mask(i,j,k) = .false. endif enddo enddo enddo do k=2,km do j=jsc, jec do i=isc, iec if (Ice%ice_mask(i,j,k)) then call ice_optics(Ice%albedo(i,j,k), Ice%pen(i,j,k), Ice%trn(i,j,k), & Ice%h_snow(i,j,k), Ice%h_ice(i,j,k), & Ice%t_surf(i,j,k)-Tfreeze, -MU_TS*Ice%s_surf(i,j) ) end if end do end do end do !!! THIS NEEDS TO BE DEFINED PROPERLY -- ARE THEY NEEDED ON RESTART ?? !!! FOR NOW, simply set all values to that of Ice%albedo, In general, !! the ability to read from restart file and initialize when none !! present should be incorporated above. if (cm2_bugs) then do k = 1, km do j = jsc, jec do i = isc, iec Ice%albedo_vis_dir(i,j,k) = Ice%albedo(i,j,k) Ice%albedo_nir_dir(i,j,k) = Ice%albedo(i,j,k) Ice%albedo_vis_dif(i,j,k) = Ice%albedo(i,j,k) Ice%albedo_nir_dif(i,j,k) = Ice%albedo(i,j,k) enddo enddo enddo else do k = 2, km do j = jsc, jec do i = isc, iec Ice%albedo_vis_dir(i,j,k) = Ice%albedo(i,j,k) Ice%albedo_nir_dir(i,j,k) = Ice%albedo(i,j,k) Ice%albedo_vis_dif(i,j,k) = Ice%albedo(i,j,k) Ice%albedo_nir_dif(i,j,k) = Ice%albedo(i,j,k) enddo enddo enddo end if do j = jsc, jec do i = isc, iec Ice%u_ocn(i,j) = u_surf_ice_bot(i,j) ! need under-ice current Ice%v_ocn(i,j) = v_surf_ice_bot(i,j) ! for water drag term enddo enddo if(reproduce_siena_201303) then call mpp_update_domains(Ice%u_ocn, Ice%v_ocn, Domain) else call mpp_update_domains(Ice%u_ocn, Ice%v_ocn, Domain, gridtype=BGRID_NE) endif ! put ocean and ice velocities into Ice%u_surf/v_surf on t-cells call uv_to_t(Ice%u_ocn, Ice%u_surf(:,:,1)) call uv_to_t(Ice%v_ocn, Ice%v_surf(:,:,1)) call uv_to_t(Ice%u_ice, Ice%u_surf(:,:,2)) call uv_to_t(Ice%v_ice, Ice%v_surf(:,:,2)) do k=1,2 do j = jsc, jec do i = isc, iec u = Ice % u_surf(i,j,k) v = Ice % v_surf(i,j,k) Ice % u_surf(i,j,k) = u*cos_rot(i,j)+v*sin_rot(i,j) ! rotate velocity from ocean if (cm2_bugs) then Ice % v_surf(i,j,k) = -v*cos_rot(i,j)+u*sin_rot(i,j) ! coord. to lat/lon coord. else Ice % v_surf(i,j,k) = v*cos_rot(i,j)-u*sin_rot(i,j) ! coord. to lat/lon coord. endif end do end do end do do k=3,km do j = jsc, jec do i = isc, iec Ice%u_surf(i,j,k) = Ice%u_surf(i,j,2) ! same ice flow on all ice partitions Ice%v_surf(i,j,k) = Ice%v_surf(i,j,2) ! end do end do end do ! ! Pre-timestep diagnostics ! if (id_sst>0) sent = send_data(id_sst, sst-Tfreeze, Ice%Time, mask=Ice%mask) if (id_sss>0) sent = send_data(id_sss, Ice%s_surf , Ice%Time, mask=Ice%mask) if (id_ssh>0) sent = send_data(id_ssh, Ice%sea_lev(isc:iec,jsc:jec), Ice%Time, mask=Ice%mask) if (id_uo >0) sent = send_data(id_uo , Ice%u_ocn(isc:iec,jsc:jec) , Ice%Time, mask=Ice%mask) if (id_vo >0) sent = send_data(id_vo , Ice%v_ocn(isc:iec,jsc:jec) , Ice%Time, mask=Ice%mask) if (id_bheat>0) sent = send_data(id_bheat, Ice%bheat, Ice%Time, mask=Ice%mask) end subroutine ice_bottom_to_ice_top !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! update_ice_model_fast - records fluxes (in Ice) and calculates ice temp. on ! ! (fast) atmospheric timestep (see coupler_main.f90) ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! subroutine update_ice_model_fast_old (Ice, Atmos_boundary_fluxes, flux_u, flux_v, u_star, & flux_sw_nir_dir, flux_sw_nir_dif, flux_sw_vis_dir, flux_sw_vis_dif,& flux_lw, flux_t, flux_q, dhdt, dedt, drdt, lprec, fprec, coszen, coszen_radwt, p_surf ) type (ice_data_type), intent(inout) :: Ice type(coupler_3d_bc_type), intent(inout) :: Atmos_boundary_fluxes real, dimension(isc:iec,jsc:jec,km), intent(in) :: flux_u, flux_v ! surface stress (+ up) real, dimension(isc:iec,jsc:jec,km), intent(in) :: u_star ! friction velocity (for ocean roughness) real, dimension(isc:iec,jsc:jec,km), intent(in) :: flux_lw ! net longwave radiation (+ down) real, dimension(isc:iec,jsc:jec,km), intent(in) :: flux_t ! sensible heat flux (+ up) real, dimension(isc:iec,jsc:jec,km), intent(in) :: flux_q ! specific humidity flux (+up) real, dimension(isc:iec,jsc:jec,km), intent(in) :: flux_sw_nir_dir ! net near IR direct shortwave radiation (+ down) real, dimension(isc:iec,jsc:jec,km), intent(in) :: flux_sw_nir_dif ! net near IR diffuse shortwave radiation (+ down) real, dimension(isc:iec,jsc:jec,km), intent(in) :: flux_sw_vis_dir ! net visible direct shortwave radiation (+ down) real, dimension(isc:iec,jsc:jec,km), intent(in) :: flux_sw_vis_dif ! net visible diffuse shortwave radiation (+ down) real, dimension(isc:iec,jsc:jec,km), intent(in) :: dhdt, dedt, drdt ! d(flux)/d(surf_temp) (+ up) real, dimension(isc:iec,jsc:jec,km), intent(in) :: lprec, fprec ! liquid & frozen precip. rate [kg/(m^2*sec)](+ down) real, dimension(isc:iec,jsc:jec,km), intent(in) :: coszen ! cosine of the zenith angle real, dimension(isc:iec,jsc:jec,km), intent(in) :: coszen_radwt ! cosine of the zenith angle, radiation weighted real, dimension(isc:iec,jsc:jec,km), intent(in), optional :: p_surf real, dimension(isc:iec,jsc:jec,km) :: flux_t_new, flux_q_new, flux_lh_new, flux_lw_new real, dimension(isc:iec,jsc:jec,km) :: flux_sw_nir_dir_new real, dimension(isc:iec,jsc:jec,km) :: flux_sw_nir_dif_new real, dimension(isc:iec,jsc:jec,km) :: flux_sw_vis_dir_new real, dimension(isc:iec,jsc:jec,km) :: flux_sw_vis_dif_new real, dimension(isc:iec,jsc:jec,km) :: flux_u_new, flux_v_new, lprec_new, fprec_new integer :: dy, sc, i, j, k real :: dt_fast, ts_new, dts, hf, hfd, latent logical :: sent real :: gmt, time_since_ae, cosz, rrsun, fracday, fracday_dt_ice, fracday_day real, dimension(isc:iec,jsc:jec) :: diurnal_factor real :: rad, cosz_day, cosz_dt_ice, rrsun_day, rrsun_dt_ice type (time_type) :: Dt_ice real, dimension(isc:iec,jsc:jec) :: cosz_alb if (id_alb>0) sent = send_data(id_alb, all_avg(Ice%albedo,Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) ! ! set up local copies of fluxes for modification ! do k = 1, km do j = jsc, jec do i = isc, iec flux_u_new(i,j,k) = flux_u(i,j,k) flux_v_new(i,j,k) = flux_v(i,j,k) flux_t_new(i,j,k) = flux_t(i,j,k) flux_q_new(i,j,k) = flux_q(i,j,k) flux_lh_new(i,j,k) = hlv*flux_q(i,j,k) flux_lw_new(i,j,k) = flux_lw(i,j,k) flux_sw_nir_dir_new(i,j,k) = flux_sw_nir_dir(i,j,k) flux_sw_nir_dif_new(i,j,k) = flux_sw_nir_dif(i,j,k) flux_sw_vis_dir_new(i,j,k) = flux_sw_vis_dir(i,j,k) flux_sw_vis_dif_new(i,j,k) = flux_sw_vis_dif(i,j,k) lprec_new(i,j,k) = lprec(i,j,k) fprec_new(i,j,k) = fprec(i,j,k) enddo enddo enddo if (add_diurnal_sw .or. do_sun_angle_for_alb) then !--------------------------------------------------------------------- ! extract time of day (gmt) from time_type variable time with ! function universal_time. !--------------------------------------------------------------------- gmt = universal_time(Ice%Time) !--------------------------------------------------------------------- ! extract the time of year relative to the northern hemisphere ! autumnal equinox (time_since_ae) from time_type variable ! time using the function orbital_time. !--------------------------------------------------------------------- time_since_ae = orbital_time(Ice%Time) !-------------------------------------------------------------------- ! call diurnal_solar_2d to calculate astronomy fields ! convert geo_lon and geo_lat to radians ! Per Rick Hemler: ! call daily_mean_solar to get cosz (over a day) ! call diurnal_solar with dtime=Dt_ice to get cosz over Dt_ice ! diurnal_factor = cosz_dt_ice*fracday_dt_ice*rrsun_dt_ice/cosz_day*fracday_day*rrsun_day !-------------------------------------------------------------------- rad = acos(-1.)/180. Dt_ice = Ice%Time_step_fast endif if (add_diurnal_sw) then do j = jsc, jec do i = isc, iec call diurnal_solar(geo_lat(i,j)*rad, geo_lon(i,j)*rad, Ice%time, cosz=cosz_dt_ice, & fracday=fracday_dt_ice, rrsun=rrsun_dt_ice, dt_time=Dt_ice) call daily_mean_solar (geo_lat(i,j)*rad, time_since_ae, cosz_day, fracday_day, rrsun_day) diurnal_factor(i,j) = cosz_dt_ice*fracday_dt_ice*rrsun_dt_ice / & max(1e-30, cosz_day*fracday_day*rrsun_day) enddo enddo do k = 1, km do j = jsc, jec do i = isc, iec flux_sw_nir_dir_new(i,j,k) = flux_sw_nir_dir(i,j,k) * diurnal_factor(i,j) flux_sw_nir_dif_new(i,j,k) = flux_sw_nir_dif(i,j,k) * diurnal_factor(i,j) flux_sw_vis_dir_new(i,j,k) = flux_sw_vis_dir(i,j,k) * diurnal_factor(i,j) flux_sw_vis_dif_new(i,j,k) = flux_sw_vis_dif(i,j,k) * diurnal_factor(i,j) enddo enddo enddo endif Ice%p_surf = 0.0 if (present(p_surf)) & Ice%p_surf = all_avg(p_surf,Ice%part_size(isc:iec,jsc:jec,:)) ! ! implicit update of ice surface temperature ! call get_time(Ice%Time_step_fast,sc,dy); dt_fast = 864e2*dy+sc do k = 2, km do j = jsc, jec do i=isc, iec if (Ice%ice_mask(i,j,k)) then if (Ice%h_snow(i,j,k)>0.0 .or. slab_ice) then flux_lh_new(i,j,k) = flux_lh_new(i,j,k) + hlf*flux_q(i,j,k) latent = hlv+hlf else flux_lh_new(i,j,k) = flux_lh_new(i,j,k)+hlf*flux_q(i,j,k)*(1-TFI/Ice%t_ice1(i,j,k)) latent = hlv+hlf*(1-TFI/Ice%t_ice1(i,j,k)) end if hfd = dhdt(i,j,k) + dedt(i,j,k)*latent + drdt(i,j,k) hf = flux_t(i,j,k) + flux_q(i,j,k)*latent - flux_lw(i,j,k) & - (1-Ice%pen(i,j,k))* (flux_sw_vis_dir(i,j,k)+flux_sw_vis_dif(i,j,k)+ & flux_sw_nir_dir(i,j,k)+flux_sw_nir_dif(i,j,k)) & - hfd*(Ice%t_surf(i,j,k)-Tfreeze) call ice3lay_temp(Ice%h_snow(i,j,k), Ice%h_ice(i,j,k), Ice%t_ice1(i,j,k), & Ice%t_ice2(i,j,k), ts_new, hf, hfd, Ice%pen(i,j,k)* & (1-Ice%trn(i,j,k))*(flux_sw_vis_dir(i,j,k)+ & flux_sw_vis_dif(i,j,k)+flux_sw_nir_dir(i,j,k)+ & flux_sw_nir_dif(i,j,k)), -MU_TS*Ice%s_surf(i,j), & Ice%bheat(i,j), dt_fast, Ice%tmelt(i,j,k), Ice%bmelt(i,j,k)) dts = ts_new-(Ice%t_surf(i,j,k)-Tfreeze) Ice%t_surf(i,j,k) = Ice%t_surf(i,j,k) + dts flux_t_new(i,j,k) = flux_t_new(i,j,k) + dts * dhdt(i,j,k) flux_q_new(i,j,k) = flux_q_new(i,j,k) + dts * dedt(i,j,k) flux_lh_new(i,j,k) = flux_lh_new(i,j,k) + dts * dedt(i,j,k) * latent flux_lw_new(i,j,k) = flux_lw_new(i,j,k) - dts * drdt(i,j,k) end if end do end do end do call compute_ocean_roughness (Ice%mask, u_star(:,:,1), Ice%rough_mom(:,:,1), & Ice%rough_heat(:,:,1), Ice%rough_moist(:,:,1) ) if(cm2_bugs) then call compute_ocean_albedo (Ice%mask, coszen(:,:,1), Ice%albedo(:,:,1), latitude ) Ice%albedo_vis_dir(:,:,1) = Ice%albedo(:,:,1) Ice%albedo_nir_dir(:,:,1) = Ice%albedo(:,:,1) Ice%albedo_vis_dif(:,:,1) = Ice%albedo(:,:,1) Ice%albedo_nir_dif(:,:,1) = Ice%albedo(:,:,1) elseif (do_sun_angle_for_alb) then call diurnal_solar(geo_lat*rad, geo_lon*rad, Ice%time, cosz=cosz_alb, & fracday=diurnal_factor, rrsun=rrsun_dt_ice, dt_time=Dt_ice) !diurnal_factor as dummy call compute_ocean_albedo (Ice%mask, cosz_alb(:,:), Ice%albedo_vis_dir(:,:,1),& Ice%albedo_vis_dif(:,:,1), Ice%albedo_nir_dir(:,:,1),& Ice%albedo_nir_dif(:,:,1), latitude ) else call compute_ocean_albedo (Ice%mask, coszen(:,:,1), Ice%albedo_vis_dir(:,:,1),& Ice%albedo_vis_dif(:,:,1), Ice%albedo_nir_dir(:,:,1),& Ice%albedo_nir_dif(:,:,1), latitude, coszen_radwt(:,:,1) ) endif if (id_alb_vis_dir>0) sent = send_data(id_alb_vis_dir, all_avg(Ice%albedo_vis_dir,Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) if (id_alb_vis_dif>0) sent = send_data(id_alb_vis_dif, all_avg(Ice%albedo_vis_dif,Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) if (id_alb_nir_dir>0) sent = send_data(id_alb_nir_dir, all_avg(Ice%albedo_nir_dir,Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) if (id_alb_nir_dif>0) sent = send_data(id_alb_nir_dif, all_avg(Ice%albedo_nir_dif,Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) call sum_top_quantities ( Ice, Atmos_boundary_fluxes, flux_u_new, flux_v_new, flux_t_new, & flux_q_new, flux_sw_nir_dir_new, flux_sw_nir_dif_new, flux_sw_vis_dir_new, & flux_sw_vis_dif_new, flux_lw_new, lprec_new, fprec_new, flux_lh_new ) Ice%Time = Ice%Time + Ice%Time_step_fast ! advance time end subroutine update_ice_model_fast_old !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! update_ice_model_slow - do ice dynamics, transport, and mass changes ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! subroutine update_ice_model_slow (Ice, runoff, calving, & runoff_hflx, calving_hflx, p_surf) type (ice_data_type), intent(inout) :: Ice real, dimension(isc:iec,jsc:jec), intent(in), optional :: runoff, calving real, dimension(isc:iec,jsc:jec), intent(in), optional :: runoff_hflx, calving_hflx real, dimension(:,:,:), intent(in), optional :: p_surf ! obsolete real, dimension(isc:iec,jsc:jec) :: fx_ice, fy_ice real, dimension(isc:iec,jsc:jec) :: fx_cor, fy_cor real, dimension(isc:iec,jsc:jec) :: fx_wat, fy_wat real, dimension(isc:iec,jsc:jec) :: hi_change, h2o_change, bsnk, x real, dimension(isc:iec,jsc:jec,2:km) :: snow_to_ice real, dimension(isc:iec,jsc:jec,km) :: part_save, part_save_uv real, dimension(isc:iec,jsc:jec) :: dum1, Obs_h_ice ! for qflux calculation real, dimension(isc:iec,jsc:jec,2) :: Obs_cn_ice ! partition 2 = ice concentration real, dimension(isd:ied,jsd:jed) :: tmp1, tmp2 real, dimension(2:km) :: e2m integer :: i, j, k, l, sc, dy, iyr, imon, iday, ihr, imin, isec real :: dt_slow, heat_to_ocn, h2o_to_ocn, h2o_from_ocn, sn2ic, bablt real :: heat_limit_ice, heat_res_ice real :: tot_heat, heating, tot_frazil logical :: sent real, parameter :: LI = hlf call get_time(Ice%Time_step_slow,sc,dy); dt_slow = 864e2*dy+sc ! ! Fluxes ! if (present( runoff)) then ! save liquid runoff for ocean do j = jsc, jec do i = isc, iec Ice % runoff(i,j) = runoff(i,j) enddo enddo if (id_runoff>0) & sent = send_data( id_runoff, Ice%runoff, Ice%Time, mask=Ice%mask ) else do j = jsc, jec do i = isc, iec Ice % runoff(i,j) = 0.0 enddo enddo endif if (present(calving)) then ! save frozen runoff for ocean do j = jsc, jec do i = isc, iec Ice % calving(i,j) = calving(i,j) enddo enddo if (id_calving>0) & sent = send_data(id_calving, Ice%calving, Ice%Time, mask=Ice%mask ) else do j = jsc, jec do i = isc, iec Ice % calving(i,j) = 0.0 enddo enddo endif if (present( runoff_hflx)) then ! save liquid runoff hflx for ocean do j = jsc, jec do i = isc, iec Ice % runoff_hflx(i,j) = runoff_hflx(i,j) enddo enddo if (id_runoff_hflx>0) & sent = send_data( id_runoff_hflx, Ice%runoff_hflx, Ice%Time, mask=Ice%mask ) else do j = jsc, jec do i = isc, iec Ice % runoff_hflx(i,j) = 0.0 enddo enddo endif if (present(calving_hflx)) then ! save frozen runoff hflx for ocean do j = jsc, jec do i = isc, iec Ice % calving_hflx(i,j) = calving_hflx(i,j) enddo enddo if (id_calving_hflx>0) & sent = send_data(id_calving_hflx, Ice%calving_hflx, Ice%Time, mask=Ice%mask ) else do j = jsc, jec do i = isc, iec Ice % calving_hflx(i,j) = 0.0 enddo enddo endif !TOM> assume that open water area is not up to date: call mpp_clock_end(iceClock) call mpp_clock_end(iceClock2) tmp1=1.-max(1.-sum(Ice%part_size(:,:,2:km),dim=3),0.0) tmp2=ice_avg(Ice%h_ice,Ice%part_size) ! Calve off icebergs and integrate forward iceberg trajectories if (do_icebergs) call icebergs_run( Ice%icebergs,Ice%Time, & Ice%calving, Ice%u_ocn, Ice%v_ocn, Ice%u_ice, Ice%v_ice, & Ice%flux_u, Ice%flux_v, Ice%sea_lev, Ice%t_surf(:,:,1), & Ice%calving_hflx, tmp1, tmp2) call mpp_clock_begin(iceClock2) call mpp_clock_begin(iceClock) call avg_top_quantities(Ice) ! average fluxes from update_ice_model_fast do k = 1, km do j = jsc, jec do i = isc, iec part_save(i,j,k) = Ice%part_size(i,j,k) part_save_uv(i,j,k) = Ice%part_size_uv(i,j,k) enddo enddo enddo ! ! conservation checks: top fluxes ! call mpp_clock_begin(iceClock7) if (conservation_check) then h2o(2) = h2o(2)+dt_slow*sum(cell_area*(Ice%runoff+Ice%calving & +all_avg(Ice%lprec_top+Ice%fprec_top-Ice%flux_q_top,Ice%part_size(isc:iec,jsc:jec,:))),& cell_area > 0 ) heat(2) = heat(2)+dt_slow*sum(cell_area*( & all_avg(Ice%flux_sw_vis_dir_top+Ice%flux_sw_vis_dif_top & +Ice%flux_sw_nir_dir_top+Ice%flux_sw_nir_dif_top, & Ice%part_size(isc:iec,jsc:jec,:)) & +all_avg(Ice%flux_lw_top, Ice%part_size(isc:iec,jsc:jec,:)) & -all_avg(Ice%flux_t_top , Ice%part_size(isc:iec,jsc:jec,:)) & -all_avg(Ice%flux_lh_top, Ice%part_size(isc:iec,jsc:jec,:)) & -LI*(all_avg(Ice%fprec_top,Ice%part_size(isc:iec,jsc:jec,:)) & +Ice%calving)), cell_area > 0 ) tot_frazil = sum(cell_area*Ice%frazil) end if call mpp_clock_end(iceClock7) ! Dynamics ! call mpp_clock_begin(iceClock4) tmp1 = ice_avg(Ice%h_snow,Ice%part_size) tmp2 = ice_avg(Ice%h_ice,Ice%part_size) call mpp_clock_begin(iceClocka) if(reproduce_siena_201303) then call ice_dynamics(1-Ice%part_size(:,:,1), tmp1, tmp2, Ice%u_ice, Ice%v_ice, & Ice%sig11, Ice%sig22, Ice%sig12, Ice%u_ocn, Ice%v_ocn, & ice_avg(Ice%flux_u_top_bgrid(isc:iec,jsc:jec,:),Ice%part_size(isc:iec,jsc:jec,:) ), & ice_avg(Ice%flux_v_top_bgrid(isc:iec,jsc:jec,:),Ice%part_size(isc:iec,jsc:jec,:) ), & Ice%sea_lev, fx_wat, fy_wat, fx_ice, fy_ice, fx_cor, fy_cor) else call ice_dynamics(1-Ice%part_size(:,:,1), tmp1, tmp2, Ice%u_ice, Ice%v_ice, & Ice%sig11, Ice%sig22, Ice%sig12, Ice%u_ocn, Ice%v_ocn, & ice_avg(Ice%flux_u_top_bgrid(isc:iec,jsc:jec,:),Ice%part_size_uv(isc:iec,jsc:jec,:) ), & ice_avg(Ice%flux_v_top_bgrid(isc:iec,jsc:jec,:),Ice%part_size_uv(isc:iec,jsc:jec,:) ), & Ice%sea_lev, fx_wat, fy_wat, fx_ice, fy_ice, fx_cor, fy_cor) endif call mpp_clock_end(iceClocka) call mpp_clock_begin(iceClockb) if(reproduce_siena_201303) then call mpp_update_domains(Ice%u_ice, Ice%v_ice, Domain) if(tripolar_grid) then call cut_check('u_ice', Ice%u_ice) ! these calls fix round off differences call cut_check('v_ice', Ice%v_ice) ! in northernmost velocities over the fold call mpp_update_domains(Ice%u_ice, Ice%v_ice, Domain) endif else call mpp_update_domains(Ice%u_ice, Ice%v_ice, Domain, gridtype=BGRID_NE) endif call mpp_clock_end(iceClockb) call mpp_clock_begin(iceClockc) ! ! Dynamics diagnostics ! if (id_fax>0) & sent = send_data(id_fax, all_avg(Ice%flux_u_top_bgrid(isc:iec,jsc:jec,:),Ice%part_size_uv), Ice%Time) if (id_fay>0) & sent = send_data(id_fay, all_avg(Ice%flux_v_top_bgrid(isc:iec,jsc:jec,:),Ice%part_size_uv), Ice%Time) if (id_fix>0) sent = send_data(id_fix, fx_ice, Ice%Time) if (id_fiy>0) sent = send_data(id_fiy, fy_ice, Ice%Time) if (id_fcx>0) sent = send_data(id_fcx, fx_cor, Ice%Time) if (id_fcy>0) sent = send_data(id_fcy, fy_cor, Ice%Time) if (id_fwx>0) sent = send_data(id_fwx, -fx_wat, Ice%Time) ! water force on ice if (id_fwy>0) sent = send_data(id_fwy, -fy_wat, Ice%Time) ! ...= -ice on water if (id_strna>0) & sent = send_data(id_strna, strain_angle(Ice%u_ice,Ice%v_ice), Ice%Time, mask=Ice%mask) if (id_sigi>0) sent = send_data(id_sigi, sigI(ice_avg(Ice%h_ice(isc:iec,jsc:jec,:), & Ice%part_size(isc:iec,jsc:jec,:)),1-Ice%part_size(isc:iec,jsc:jec,1), & Ice%sig11(isc:iec,jsc:jec),Ice%sig22(isc:iec,jsc:jec),Ice%sig12(isc:iec,jsc:jec)), & Ice%Time, mask=Ice%mask) if (id_sigii>0) sent = send_data(id_sigii, sigII(ice_avg(Ice%h_ice(isc:iec,jsc:jec,:), & Ice%part_size(isc:iec,jsc:jec,:)),1-Ice%part_size(isc:iec,jsc:jec,1), & Ice%sig11(isc:iec,jsc:jec),Ice%sig22(isc:iec,jsc:jec),Ice%sig12(isc:iec,jsc:jec)), & Ice%Time, mask=Ice%mask) if (id_stren>0) sent = send_data(id_stren, ice_strength(ice_avg(Ice%h_ice(isc:iec,jsc:jec,:), & Ice%part_size(isc:iec,jsc:jec,:)), 1-Ice%part_size(isc:iec,jsc:jec,1)), Ice%Time, mask=Ice%mask) if (id_ui>0) sent = send_data(id_ui, Ice%u_ice(isc:iec,jsc:jec), Ice%Time) if (id_vi>0) sent = send_data(id_vi, Ice%v_ice(isc:iec,jsc:jec), Ice%Time) do k=2,km do j = jsc, jec do i = isc, iec Ice%flux_u_top_bgrid(i,j,k) = fx_wat(i,j) ! stress of ice on ocean Ice%flux_v_top_bgrid(i,j,k) = fy_wat(i,j) ! enddo enddo end do call mpp_clock_end(iceClockc) call mpp_clock_end(iceClock4) ! ! Thermodynamics ! call mpp_clock_begin(iceClock5) if (id_frazil>0) sent = send_data(id_frazil, Ice%frazil/dt_slow, Ice%Time, mask=Ice%mask) snow_to_ice = 0 bsnk = 0 hi_change = all_avg(Ice%h_ice(isc:iec,jsc:jec,:), Ice%part_size(isc:iec,jsc:jec,:)) h2o_change = all_avg(DS*Ice%h_snow(isc:iec,jsc:jec,:)+DI*Ice%h_ice(isc:iec,jsc:jec,:),Ice%part_size(isc:iec,jsc:jec,:)) ! get observed ice thickness for ice restoring, if calculating qflux if (do_ice_restore) & call get_sea_surface(Ice%Time, dum1, Obs_cn_ice, Obs_h_ice) do k=2,km do j=jsc, jec do i=isc, iec if (cell_area(i,j) > 0 .and.Ice%h_ice(i,j,k) > 0) then ! reshape the ice based on fluxes h2o_from_ocn = 0.0; h2o_to_ocn = 0.0; heat_to_ocn = 0.0 call ice3lay_resize(Ice%h_snow(i,j,k), Ice%h_ice(i,j,k), & Ice%t_ice1(i,j,k), Ice%t_ice2(i,j,k), & Ice%fprec_top(i,j,k) *dt_slow, 0.0, & Ice%flux_q_top(i,j,k)*dt_slow, & Ice%tmelt (i,j,k), Ice%bmelt(i,j,k), & -MU_TS*Ice%s_surf(i,j), & heat_to_ocn, h2o_to_ocn, h2o_from_ocn, & snow_to_ice(i,j,k), bablt ) ! modify above-ice to under-ice fluxes for passing to ocean Ice%flux_q_top (i,j,k) = h2o_from_ocn/dt_slow ! no ice, evaporation left Ice%flux_lh_top(i,j,k) = hlv*Ice%flux_q_top(i,j,k) Ice%flux_lw_top(i,j,k) = 0.0 Ice%flux_t_top (i,j,k) = Ice%bheat(i,j)-heat_to_ocn/dt_slow Ice%flux_sw_vis_dif_top(i,j,k) = (Ice%flux_sw_vis_dir_top(i,j,k)+ & Ice%flux_sw_vis_dif_top(i,j,k)+Ice%flux_sw_nir_dir_top(i,j,k)+ & Ice%flux_sw_nir_dif_top(i,j,k))*Ice%pen(i,j,k)*Ice%trn(i,j,k) Ice%flux_sw_nir_dir_top(i,j,k) = 0.0 Ice%flux_sw_nir_dif_top(i,j,k) = 0.0 Ice%flux_sw_vis_dir_top(i,j,k) = 0.0 Ice%fprec_top (i,j,k) = 0.0 Ice%lprec_top (i,j,k) = Ice%lprec_top(i,j,k) + h2o_to_ocn/dt_slow bsnk(i,j) = bsnk(i,j) - Ice%part_size(i,j,k)*bablt ! bot. melt. ablation end if ! ! absorb frazil in thinest ice partition available ! if (Ice%frazil(i,j)>0.and.Ice%part_size(i,j,1)+Ice%part_size(i,j,k)>0) then Ice%h_snow (i,j,k) = Ice%h_snow (i,j,k) * Ice%part_size(i,j,k) Ice%h_ice (i,j,k) = Ice%h_ice (i,j,k) * Ice%part_size(i,j,k) Ice%t_surf (i,j,k) = (Ice%t_surf(i,j,k) *Ice%part_size(i,j,k) & + (Tfreeze-MU_TS*Ice%s_surf(i,j))*Ice%part_size(i,j,1)) Ice%part_size(i,j,k) = Ice%part_size(i,j,k) + Ice%part_size(i,j,1) Ice%part_size(i,j,1) = 0.0 Ice%h_snow (i,j,k) = Ice%h_snow (i,j,k) / Ice%part_size(i,j,k) Ice%h_ice (i,j,k) = Ice%h_ice (i,j,k) / Ice%part_size(i,j,k) Ice%t_surf (i,j,k) = Ice%t_surf (i,j,k) / Ice%part_size(i,j,k) call ice3lay_resize(Ice%h_snow(i,j,k), Ice%h_ice (i,j,k), & Ice%t_ice1(i,j,k), Ice%t_ice2(i,j,k), 0.0, & Ice%frazil(i,j)/Ice%part_size(i,j,k), 0.0, 0.0, 0.0, & -MU_TS*Ice%s_surf(i,j), & heat_to_ocn, h2o_to_ocn, h2o_from_ocn, sn2ic) Ice%frazil(i,j) = 0.0; ! ! spread frazil salinification over all partitions ! Ice%lprec_top (i,j,:) = Ice%lprec_top(i,j,:) + h2o_to_ocn*Ice%part_size(i,j,k)/dt_slow end if end do ! i-loop end do ! j-loop end do ! k-loop call mpp_clock_end(iceClock5) ! ! Calculate QFLUX's from (1) restoring to obs and (2) limiting total ice. ! call mpp_clock_begin(iceClock6) if (do_ice_restore .or. do_ice_limit) then do j = jsc, jec do i = isc, iec Ice % qflx_lim_ice(i,j) = 0.0 Ice % qflx_res_ice(i,j) = 0.0 enddo enddo do j=jsc, jec do i=isc, iec heat_res_ice = 0.0 heat_limit_ice = 0.0 ! ! calculate enthalpy ! if (slab_ice) then e2m(2) = Ice%h_ice(i,j,2)*DI*LI else do k=2,km if ((Ice%part_size(i,j,k)>0.0.and.Ice%h_ice(i,j,k)>0.0)) then e2m(k) = e_to_melt(Ice%h_snow(i,j,k), Ice%h_ice(i,j,k)/2, & Ice%t_ice1(i,j,k), Ice%h_ice(i,j,k)/2, & Ice%t_ice2(i,j,k)) * Ice%part_size(i,j,k) else e2m(k) = 0.0 end if end do end if ! ! calculate heat needed to constrain ice enthalpy ! if (do_ice_restore) then ! TK Mod: restore to observed enthalpy (no change for slab, but for ! sis ice, results in restoring toward thickness * concentration ! TK Mod for test 10/18/02 ! Concentration is 1.0 where there is ice for slab case. ! The input field Obs_cn_ice may have values less than 1, ! so put in if test... if (slab_ice) then heat_res_ice = -(LI*DI*Obs_h_ice(i,j)-sum(e2m)) & *dt_slow/(86400*ice_restore_timescale) else heat_res_ice = -(LI*DI*Obs_h_ice(i,j)*Obs_cn_ice(i,j,2)-sum(e2m)) & *dt_slow/(86400*ice_restore_timescale) end if ! heat_res_ice = -(LI*DI*Obs_h_ice(i,j)-sum(e2m)) & ! *dt_slow/(86400*ice_restore_timescale) end if if (do_ice_limit .and. (sum(e2m) > max_ice_limit*DI*LI)) then heat_limit_ice = sum(e2m)-LI*DI*max_ice_limit ! should we "heat_ice_res = 0.0" ? end if ! ! apply constraining heat to ice ! tot_heat = heat_res_ice+heat_limit_ice if (slab_ice) Ice%h_ice (i,j,2) = Ice%h_ice (i,j,2) - tot_heat/(DI*LI) if (.not. slab_ice .and. (tot_heat>0.0)) then ! add like ocean-ice heat do k=2,km if (e2m(k) > 0.0) then heating = tot_heat/sum(Ice%part_size(i,j,k:km)) if (heating*Ice%part_size(i,j,k) > e2m(k)) then ! cat. melts away Ice%h_ice (i,j,k) = 0.0 Ice%h_snow(i,j,k) = 0.0 tot_heat = tot_heat - e2m(k) else h2o_from_ocn = 0.0; h2o_to_ocn = 0.0; heat_to_ocn = 0.0 call ice3lay_resize(Ice%h_snow(i,j,k), Ice%h_ice(i,j,k), & Ice%t_ice1(i,j,k), Ice%t_ice2(i,j,k), & 0.0, 0.0, 0.0, 0.0, heating, & -MU_TS*Ice%s_surf(i,j), & heat_to_ocn, h2o_to_ocn, h2o_from_ocn, & snow_to_ice(i,j,k), bablt ) tot_heat = tot_heat - heating*Ice%part_size(i,j,k) end if end if end do end if tot_heat = heat_res_ice+heat_limit_ice if (.not. slab_ice .and. (tot_heat<0.0)) then ! add like frazil do k=2,km if (Ice%part_size(i,j,1)+Ice%part_size(i,j,k)>0) exit end do ! k is thinnest ice partition that can recieve frazil Ice%h_snow (i,j,k) = Ice%h_snow (i,j,k) * Ice%part_size(i,j,k) Ice%h_ice (i,j,k) = Ice%h_ice (i,j,k) * Ice%part_size(i,j,k) Ice%t_surf (i,j,k) = (Ice%t_surf(i,j,k) *Ice%part_size(i,j,k) & +(Tfreeze-MU_TS*Ice%s_surf(i,j))*Ice%part_size(i,j,1)) Ice%part_size(i,j,k) = Ice%part_size(i,j,k) + Ice%part_size(i,j,1) Ice%part_size(i,j,1) = 0.0 Ice%h_snow (i,j,k) = Ice%h_snow (i,j,k) / Ice%part_size(i,j,k) Ice%h_ice (i,j,k) = Ice%h_ice (i,j,k) / Ice%part_size(i,j,k) Ice%t_surf (i,j,k) = Ice%t_surf (i,j,k) / Ice%part_size(i,j,k) call ice3lay_resize(Ice%h_snow(i,j,k), Ice%h_ice (i,j,k), & Ice%t_ice1(i,j,k), Ice%t_ice2(i,j,k), 0.0, & -tot_heat/Ice%part_size(i,j,k), 0.0, 0.0, 0.0, & -MU_TS*Ice%s_surf(i,j), & heat_to_ocn, h2o_to_ocn, h2o_from_ocn, sn2ic) end if ! Convert constraining heat from energy (J/m^2) to flux (W/m^2) Ice%qflx_lim_ice(i,j) = heat_limit_ice / dt_slow Ice%qflx_res_ice(i,j) = heat_res_ice / dt_slow ! ! Check for energy conservation ! if (slab_ice) then e2m(2) = e2m(2) - Ice%h_ice(i,j,2)*DI*LI else do k=2,km if (Ice%part_size(i,j,k)>0.0.and.Ice%h_ice(i,j,k)>0.0) & e2m(k) = e2m(k)-e_to_melt(Ice%h_snow(i,j,k), Ice%h_ice(i,j,k)/2, & Ice%t_ice1(i,j,k), Ice%h_ice(i,j,k)/2, & Ice%t_ice2(i,j,k)) * Ice%part_size(i,j,k) end do end if ! if (abs(sum(e2m) - heat_res_ice - heat_limit_ice)>DI*LI*1e-3) & ! print *, 'QFLUX conservation error at', i, j, 'heat2ice=', & ! tot_heat, 'melted=', sum(e2m), 'h*part_size=', & ! Ice%h_ice(i,j,:)*Ice%part_size(i,j,:) end do end do end if call mpp_clock_end(iceClock6) ! ! Salt fluxes to ocean ! call mpp_clock_begin(iceClock8) hi_change = all_avg(Ice%h_ice(isc:iec,jsc:jec,:), Ice%part_size(isc:iec,jsc:jec,:))-hi_change Ice%flux_salt = ice_bulk_salin*DI*hi_change/dt_slow if (id_saltf>0) sent = send_data(id_saltf, Ice%flux_salt, Ice%Time, mask=Ice%mask) h2o_change = all_avg(DS*Ice%h_snow(isc:iec,jsc:jec,:)+DI*Ice%h_ice(isc:iec,jsc:jec,:), & Ice%part_size(isc:iec,jsc:jec,:))-h2o_change if (id_lsnk>0) then x = h2o_change*864e2*365/dt_slow do j = jsc, jec do i = isc, iec if(x(i,j)>0.0) x(i,j) = 0.0 enddo enddo sent = send_data(id_lsnk, x, Ice%Time, mask=Ice%mask) end if if (id_lsrc>0) then x = h2o_change*864e2*365/dt_slow do j = jsc, jec do i = isc, iec if(x(i,j)<0.0) x(i,j) = 0.0 enddo enddo sent = send_data(id_lsrc, x, Ice%Time, mask=Ice%mask) end if if (id_bsnk>0) sent = send_data(id_bsnk, bsnk*864e2*365/dt_slow, Ice%Time, mask=Ice%mask) if (id_tmelt>0) sent = send_data(id_tmelt, ice_avg(Ice%tmelt,Ice%part_size(isc:iec,jsc:jec,:))/dt_slow, & Ice%Time, mask=Ice%mask) if (id_bmelt>0) sent = send_data(id_bmelt, ice_avg(Ice%bmelt,Ice%part_size(isc:iec,jsc:jec,:))/dt_slow, & Ice%Time, mask=Ice%mask) if (id_sn2ic>0) sent = send_data(id_sn2ic, all_avg(snow_to_ice,Ice%part_size(isc:iec,jsc:jec,:))/dt_slow, & Ice%Time, mask=Ice%mask ) if (id_qflim>0) sent = send_data(id_qflim, Ice%qflx_lim_ice, Ice%Time, mask=Ice%mask) if (id_qfres>0) sent = send_data(id_qfres, Ice%qflx_res_ice, Ice%Time, mask=Ice%mask) ! ! Ice transport ... all ocean fluxes have been calculated by now ! h2o_change = all_avg(DS*Ice%h_snow(isc:iec,jsc:jec,:)+DI*Ice%h_ice(isc:iec,jsc:jec,:),Ice%part_size(isc:iec,jsc:jec,:)) call transport(Ice) x = all_avg(DS*Ice%h_snow(isc:iec,jsc:jec,:)+DI*Ice%h_ice(isc:iec,jsc:jec,:),Ice%part_size(isc:iec,jsc:jec,:)) if (id_mi>0) sent = send_data(id_mi, x, Ice%Time, mask = Ice%mask) do j = jsc, jec do i = isc, iec Ice%mi(i,j) = x(i,j) enddo enddo if (do_icebergs) call icebergs_incr_mass(Ice%icebergs, x) ! Add icebergs mass in kg/m^2 if (id_mib>0) sent = send_data(id_mib, x, Ice%Time, mask = Ice%mask) ! Diagnose total mass if (id_slp>0) sent = send_data(id_slp, Ice%p_surf, Ice%Time, mask=Ice%mask) do j = jsc, jec do i = isc, iec if (slp2ocean) then Ice%p_surf(i,j) = Ice%p_surf(i,j)-1e5 else Ice%p_surf(i,j) = 0.0 endif Ice%p_surf(i,j) = Ice%p_surf(i,j) + grav*x(i,j) enddo enddo h2o_change = x-h2o_change if (spec_ice) then ! over-write changes with specifications call get_sea_surface(Ice%Time, Ice%t_surf(:,:,1), Ice%part_size(isc:iec,jsc:jec,:), Ice%h_ice (isc:iec,jsc:jec,2)) call mpp_update_domains(Ice%part_size, Domain) endif do j = jsd, jed do i = isd, ied Ice%part_size (i,j,1) = 1.0 enddo enddo do j = jsc, jec do i = isc, iec Ice%part_size_uv(i,j,1) = 1.0 enddo enddo do k=2,km do j = jsd, jed do i = isd, ied Ice%part_size(i,j,1) = Ice%part_size (i,j,1) - Ice%part_size (i,j,k) enddo enddo enddo do k=2,km call t_to_uv(Ice%part_size(:,:,k),Ice%part_size_uv(:,:,k)) Ice%part_size_uv(:,:,1) = Ice%part_size_uv(:,:,1) - Ice%part_size_uv(:,:,k) end do call mpp_clock_end(iceClock8) ! ! Thermodynamics diagnostics ! call mpp_clock_begin(iceClock9) if (id_cn>0) sent = send_data(id_cn, Ice%part_size(isc:iec,jsc:jec,2:km), Ice%Time, mask=spread(Ice%mask,3,km-1) ) ! TK Mod: 10/18/02 ! if (id_obs_cn>0) sent = send_data(id_obs_cn, Obs_cn_ice(:,:,2), Ice%Time, & ! mask=Ice%mask ) if (id_ext>0) sent = send_data(id_ext, ext(Ice%part_size(isc:iec,jsc:jec,1)), Ice%Time, mask=Ice%mask) if (id_wnd>0) sent = send_data(id_wnd, Ice%wnd(isc:iec,jsc:jec,1), Ice%Time, mask=Ice%mask) if (id_hs>0) sent = send_data(id_hs, ice_avg(Ice%h_snow(isc:iec,jsc:jec,:),Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) if (id_hi>0) sent = send_data(id_hi, ice_avg(Ice%h_ice(isc:iec,jsc:jec,:) ,Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) ! TK Mod: 10/18/02: (commented out...does not compile yet... add later ! if (id_obs_hi>0) & ! sent = send_data(id_obs_hi, ice_avg(Obs_h_ice,Ice%part_size), Ice%Time, & ! mask=Ice%mask) if (id_ts>0) sent = send_data(id_ts, ice_avg(Ice%t_surf-Tfreeze,Ice%part_size(isc:iec,jsc:jec,:)), Ice%Time, mask=Ice%mask) if (id_t1>0) sent = send_data(id_t1, ice_avg(Ice%t_ice1(isc:iec,jsc:jec,:),Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) if (id_t2>0) sent = send_data(id_t2, ice_avg(Ice%t_ice2(isc:iec,jsc:jec,:),Ice%part_size(isc:iec,jsc:jec,:)), & Ice%Time, mask=Ice%mask) if (id_xprt>0) sent = send_data(id_xprt, h2o_change*864e2*365/dt_slow, Ice%Time, mask=Ice%mask) if (id_e2m>0) then x = 0.0 do k = 2, km do j = jsc, jec do i=isc, iec if (Ice%h_ice(i,j,k)>0.0) & x(i,j) = x(i,j) + Ice%part_size(i,j,k)*e_to_melt(Ice%h_snow(i,j,k), & Ice%h_ice(i,j,k)/2, Ice%t_ice1(i,j,k), Ice%h_ice(i,j,k)/2, Ice%t_ice2(i,j,k) ) end do end do end do sent = send_data(id_e2m, x, Ice%Time, mask=Ice%mask) end if call ice_top_to_ice_bottom(Ice, part_save, part_save_uv) ! ! conservation checks: bottom fluxes and final ! if (conservation_check) then h2o(3) = h2o(3) + dt_slow*sum(cell_area *(Ice%lprec+Ice%fprec-Ice%flux_q+Ice%runoff+Ice%calving),& cell_area > 0 ) heat(3) = heat(3) + tot_frazil + dt_slow*sum(cell_area*(Ice%flux_sw_vis_dir+Ice%flux_sw_vis_dif+ & Ice%flux_sw_nir_dir+Ice%flux_sw_nir_dif+Ice%flux_lw- & Ice%flux_t-Ice%flux_lh -LI*(Ice%fprec+Ice%calving)), cell_area > 0 ) h2o(4) = sum(cell_area*all_avg(DI*Ice%h_ice(isc:iec,jsc:jec,:)+ & DS*Ice%h_snow(isc:iec,jsc:jec,:),Ice%part_size(isc:iec,jsc:jec,:))) heat(4) = 0.0 do k=2,km do j=jsc, jec do i=isc, iec if ((Ice%part_size(i,j,k)>0.0.and.Ice%h_ice(i,j,k)>0.0)) then if (slab_ice) then heat(4) = heat(4) - cell_area(i,j) * Ice%part_size(i,j,k)*Ice%h_ice(i,j,2)*DI*LI else heat(4) = heat(4) - cell_area(i,j) * Ice%part_size(i,j,k)*e_to_melt(Ice%h_snow(i,j,k), & Ice%h_ice(i,j,k)/2, Ice%t_ice1(i,j,k), Ice%h_ice(i,j,k)/2, Ice%t_ice2(i,j,k)) end if end if end do end do end do end if call get_date(Ice%Time, iyr, imon, iday, ihr, imin, isec) call get_time(Ice%Time-set_date(iyr,1,1,0,0,0),isec,iday) if(verbose) call ice_line(iyr, iday+1, isec, Ice%part_size(isc:iec,jsc:jec,:), & Ice%t_surf(:,:,1)-Tfreeze) do j=jsc, jec do i=isc, iec if (Ice%mask(i,j).and.(abs(sum(Ice%part_size(i,j,:))-1.0)>1e-2)) & print *,'ICE%PART_SIZE=',Ice%part_size(i,j,:), 'DOES NOT SUM TO 1 AT', & geo_lon(i,j), geo_lat(i,j) end do end do call mpp_clock_end(iceClock9) end subroutine update_ice_model_slow !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! ice_redistribute - a simple ice redistribution scheme from Igor Polyakov ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! subroutine ice_redistribute(cn, hs, hi, t1, t2) !!$ real, intent(inout), dimension(2:km ) :: cn, hs, hi, t1, t2 real, intent(inout), dimension(2:) :: cn, hs, hi, t1, t2 real :: cw ! open water concentration integer :: k cw = 1-sum(cn) if (cw<0) cn(2) = cw + cn(2) ! open water has been eliminated by convergence do k=2,km-1 if (cn(k)<0) then cn(k+1) = cn(k+1)+cn(k); cn(k) = 0 ! pass concentration deficit up to hs(k+1) = hs(k+1)+hs(k); hs(k) = 0 ! next thicker category hi(k+1) = hi(k+1)+hi(k); hi(k) = 0 t1(k+1) = t1(k+1)+t1(k); t1(k) = 0 ! NOTE: here between the thm_pack and t2(k+1) = t2(k+1)+t2(k); t2(k) = 0 ! thm_unpack calls, all quantities are endif ! extensive, so we add instead of end do ! averaging do k=2,km-1 if (hi(k)>hlim(k)*cn(k)) then cn(k+1) = cn(k+1)+cn(k); cn(k) = 0 ! upper thickness limit exceeded hs(k+1) = hs(k+1)+hs(k); hs(k) = 0 ! move ice up to next thicker category hi(k+1) = hi(k+1)+hi(k); hi(k) = 0 t1(k+1) = t1(k+1)+t1(k); t1(k) = 0 t2(k+1) = t2(k+1)+t2(k); t2(k) = 0 endif end do do k=km,3,-1 if (hi(k) 0.0) then Ice%part_size(i,j,2) = 1.0 else Ice%part_size(i,j,2) = 0.0 endif enddo enddo return; end if call thm_pack(Ice%part_size(isc:iec,jsc:jec,2:km), Ice%h_snow(isc:iec,jsc:jec,:), Ice%h_ice(isc:iec,jsc:jec,:), & Ice%t_ice1(isc:iec,jsc:jec,:), Ice%t_ice2(isc:iec,jsc:jec,:) ) call mpp_update_domains(Ice%part_size, Domain) ! cannot be combined with updates below call mpp_update_domains(Ice%h_snow, Domain, complete=.false.) call mpp_update_domains(Ice%h_ice, Domain, complete=.false.) call mpp_update_domains(Ice%t_ice1, Domain, complete=.false.) call mpp_update_domains(Ice%t_ice2, Domain, complete=.true.) ! Move transporting flow to the C-grid (from the B-grid) ! Note: this block of code was moved here from within s/r ice_advect ! where the same calculations and mpp_updates were repeated for each ! class and variable. -AJA uc = 0.0; vc = 0.0 do j = jsc, jec do i = isc, iec uc(i,j) = 0.5 * ( Ice%u_ice(i,j-1) + Ice%u_ice(i,j) ) vc(i,j) = 0.5 * ( Ice%v_ice(i-1,j) + Ice%v_ice(i,j) ) enddo enddo if (channel_viscosity>0.) then ! This block of code is a parameterization of either (or both) ! i) the pressure driven oceanic flow in a narrow channel, or ! ii) pressure driven flow of ice itself. ! The latter is a speculative but both are missing due to ! masking of velocities to zero in a single-cell wide channel. tmp1=1.-max(1.-sum(Ice%part_size(:,:,2:km),dim=3),0.0) ustar(:,:)=0.; vstar(:,:)=0. ustaro(:,:)=0.; vstaro(:,:)=0. ustarv(:,:)=0.; vstarv(:,:)=0. do j = jsc, jec do i = isc, iec if (Ice%u_ice(i,j).ne.0..and.Ice%vmask(i,j)==0.) stop 'Ooops' ! Debug new vmask if (Ice%v_ice(i,j).ne.0..and.Ice%vmask(i,j)==0.) stop 'Ooops' ! Debug new vmask if (uc(i,j)==0. & ! this is a redundant test due to following line .and. Ice%vmask(i,j)+Ice%vmask(i,j-1)==0. & ! =0 => no vels .and. wett(i,j)*wett(i+1,j)>0.) then ! >0 => open for transport grad_eta=(Ice%sea_lev(i+1,j)-Ice%sea_lev(i,j)) & ! delta_i eta /(0.5*(dxv(i,j)+dxv(i,j-1))) ! /dx u_visc=-ssh_gravity*((dte(i,j)*dte(i,j))/(12.*channel_viscosity)) & ! -g*dy^2/(12*visc) *grad_eta ! d/dx eta u_ocn=sqrt( ssh_gravity*dte(i,j)*abs(grad_eta)/(36.*smag_ocn) ) ! Magnitude of ocean current u_ocn=sign(u_ocn, -grad_eta) ! Direct down the ssh gradient cnn=max(tmp1(i,j),tmp1(i+1,j))**2. ! Use the larger concentration uc(i,j)=cnn*u_visc+(1.-cnn)*u_ocn ! Limit flow to be stable for fully divergent flow if (uc(i,j)>0.) then uc(i,j)=min( uc(i,j), chan_cfl_limit*dxt(i,j)/dt_adv) else uc(i,j)=max( uc(i,j),(-1*chan_cfl_limit)*dxt(i+1,j)/dt_adv) endif if (id_ustar>0) ustar(i,j)=uc(i,j) if (id_uocean>0) ustaro(i,j)=u_ocn if (id_uchan>0) ustarv(i,j)=u_visc endif if (vc(i,j)==0. & ! this is a redundant test due to following line .and. Ice%vmask(i,j)+Ice%vmask(i-1,j)==0. & ! =0 => no vels .and. wett(i,j)*wett(i,j+1)>0.) then ! >0 => open for transport grad_eta=(Ice%sea_lev(i,j+1)-Ice%sea_lev(i,j)) & ! delta_i eta /(0.5*(dyv(i,j)+dyv(i,j-1))) ! /dy u_visc=-ssh_gravity*((dtn(i,j)*dtn(i,j))/(12.*channel_viscosity)) & ! -g*dy^2/(12*visc) *grad_eta ! d/dx eta u_ocn=sqrt( ssh_gravity*dtn(i,j)*abs(grad_eta)/(36.*smag_ocn) ) ! Magnitude of ocean current u_ocn=sign(u_ocn, -grad_eta) ! Direct down the ssh gradient cnn=max(tmp1(i,j),tmp1(i,j+1))**2. ! Use the larger concentration vc(i,j)=cnn*u_visc+(1.-cnn)*u_ocn ! Limit flow to be stable for fully divergent flow if (vc(i,j)>0.) then vc(i,j)=min( vc(i,j), chan_cfl_limit*dyt(i,j)/dt_adv) else vc(i,j)=max( vc(i,j),(-1*chan_cfl_limit)*dyt(i,j+1)/dt_adv) endif if (id_vstar>0) vstar(i,j)=vc(i,j) if (id_vocean>0) vstaro(i,j)=u_ocn if (id_vchan>0) vstarv(i,j)=u_visc endif enddo enddo endif if(reproduce_siena_201303) then call mpp_update_domains(uc, vc, Domain) else call mpp_update_domains(uc, vc, Domain, gridtype=CGRID_NE) endif uf = 0.0; vf = 0.0 do k=2,km call ice_advect(uc, vc, Ice%part_size(:,:,k)) call ice_advect(uc, vc, Ice%h_snow (:,:,k), uf0, vf0) uf = uf + DS*uf0; vf = vf + DS*vf0 call ice_advect(uc, vc, Ice%h_ice (:,:,k), uf0, vf0) uf = uf + DI*uf0; vf = vf + DI*vf0 call ice_advect(uc, vc, Ice%t_ice1 (:,:,k)) call ice_advect(uc, vc, Ice%t_ice2 (:,:,k)) end do sent = send_data(id_ix_trans, uf, Ice%Time) sent = send_data(id_iy_trans, vf, Ice%Time) do j=jsc, jec do i=isc, iec if (sum(Ice%h_ice(i,j,:))>0) & call ice_redistribute(Ice%part_size(i,j,2:km), & Ice%h_snow(i,j,:), Ice%h_ice (i,j,:), & Ice%t_ice1(i,j,:), Ice%t_ice2(i,j,:)) do k=2,km if (Ice%part_size(i,j,k)<1e-10) then Ice%h_ice (i,j,k) = 0 ! thm_unpack will zero other quantities Ice%t_surf(i,j,k) = Tfreeze-MU_TS*Ice%s_surf(i,j) end if end do end do end do call thm_unpack(Ice%part_size(isc:iec,jsc:jec,2:km), Ice%h_snow(isc:iec,jsc:jec,:), Ice%h_ice(isc:iec,jsc:jec,:), & Ice%t_ice1(isc:iec,jsc:jec,:), Ice%t_ice2(isc:iec,jsc:jec,:) ) call mpp_update_domains(Ice%part_size, Domain) ! cannot be combined with the two updates below call mpp_update_domains(Ice%h_snow, Domain, complete=.false.) call mpp_update_domains(Ice%h_ice, Domain, complete=.true.) if (id_ustar >0) sent = send_data(id_ustar , ustar(isc:iec,jsc:jec) , Ice%Time) if (id_vstar >0) sent = send_data(id_vstar , vstar(isc:iec,jsc:jec) , Ice%Time) if (id_vocean>0) sent = send_data(id_vocean, vstaro(isc:iec,jsc:jec) , Ice%Time) if (id_uocean>0) sent = send_data(id_uocean, ustaro(isc:iec,jsc:jec) , Ice%Time) if (id_vchan>0) sent = send_data(id_vchan, vstarv(isc:iec,jsc:jec) , Ice%Time) if (id_uchan>0) sent = send_data(id_uchan, ustarv(isc:iec,jsc:jec) , Ice%Time) return end subroutine transport !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! ext - modeled ice is one if ice cover >= 0.15 otherwise zero ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! function ext(cw) real, dimension(isc:,jsc:), intent(in) :: cw real, dimension(isc:iec,jsc:jec) :: ext integer :: i, j do j = jsc, jec do i = isc, iec if(cw(i,j) < 0.85) then ext(i,j) = 1 else ext(i,j) = 0 endif enddo enddo return end function ext subroutine ocn_ice_bnd_type_chksum(id, timestep, bnd_type) use fms_mod, only: stdout use mpp_mod, only: mpp_chksum character(len=*), intent(in) :: id integer , intent(in) :: timestep type(ocean_ice_boundary_type), intent(in) :: bnd_type integer :: n, m, outunit outunit = stdout() write(outunit,*) 'BEGIN CHECKSUM(ocean_ice_boundary_type):: ', id, timestep write(outunit,100) 'ocn_ice_bnd_type%u ',mpp_chksum(bnd_type%u ) write(outunit,100) 'ocn_ice_bnd_type%v ',mpp_chksum(bnd_type%v ) write(outunit,100) 'ocn_ice_bnd_type%t ',mpp_chksum(bnd_type%t ) write(outunit,100) 'ocn_ice_bnd_type%s ',mpp_chksum(bnd_type%s ) write(outunit,100) 'ocn_ice_bnd_type%frazil ',mpp_chksum(bnd_type%frazil ) write(outunit,100) 'ocn_ice_bnd_type%sea_level',mpp_chksum(bnd_type%sea_level) ! write(outunit,100) 'ocn_ice_bnd_type%data ',mpp_chksum(bnd_type%data ) 100 FORMAT("CHECKSUM::",A32," = ",Z20) do n = 1, bnd_type%fields%num_bcs !{ do m = 1, bnd_type%fields%bc(n)%num_fields !{ write(outunit,101) 'oibt%',trim(bnd_type%fields%bc(n)%name), & trim(bnd_type%fields%bc(n)%field(m)%name), & mpp_chksum(bnd_type%fields%bc(n)%field(m)%values) enddo !} m enddo !} n 101 FORMAT("CHECKSUM::",A16,a,'%',a," = ",Z20) end subroutine ocn_ice_bnd_type_chksum subroutine atm_ice_bnd_type_chksum(id, timestep, bnd_type) use fms_mod, only: stdout use mpp_mod, only: mpp_chksum character(len=*), intent(in) :: id integer , intent(in) :: timestep type(atmos_ice_boundary_type), intent(in) :: bnd_type integer :: n, outunit outunit = stdout() write(outunit,*) 'BEGIN CHECKSUM(atmos_ice_boundary_type):: ', id, timestep write(outunit,100) 'atm_ice_bnd_type%u_flux ',mpp_chksum(bnd_type%u_flux) write(outunit,100) 'atm_ice_bnd_type%v_flux ',mpp_chksum(bnd_type%v_flux) write(outunit,100) 'atm_ice_bnd_type%u_star ',mpp_chksum(bnd_type%u_star) write(outunit,100) 'atm_ice_bnd_type%t_flux ',mpp_chksum(bnd_type%t_flux) write(outunit,100) 'atm_ice_bnd_type%q_flux ',mpp_chksum(bnd_type%q_flux) write(outunit,100) 'atm_ice_bnd_type%lw_flux ',mpp_chksum(bnd_type%lw_flux) write(outunit,100) 'atm_ice_bnd_type%sw_flux_vis_dir ',mpp_chksum(bnd_type%sw_flux_vis_dir) write(outunit,100) 'atm_ice_bnd_type%sw_flux_vis_dif ',mpp_chksum(bnd_type%sw_flux_vis_dif) write(outunit,100) 'atm_ice_bnd_type%sw_flux_nir_dir ',mpp_chksum(bnd_type%sw_flux_nir_dir) write(outunit,100) 'atm_ice_bnd_type%sw_flux_nir_dif ',mpp_chksum(bnd_type%sw_flux_nir_dif) write(outunit,100) 'atm_ice_bnd_type%lprec ',mpp_chksum(bnd_type%lprec) write(outunit,100) 'atm_ice_bnd_type%fprec ',mpp_chksum(bnd_type%fprec) write(outunit,100) 'atm_ice_bnd_type%dhdt ',mpp_chksum(bnd_type%dhdt) write(outunit,100) 'atm_ice_bnd_type%dedt ',mpp_chksum(bnd_type%dedt) write(outunit,100) 'atm_ice_bnd_type%drdt ',mpp_chksum(bnd_type%drdt) write(outunit,100) 'atm_ice_bnd_type%coszen ',mpp_chksum(bnd_type%coszen) write(outunit,100) 'atm_ice_bnd_type%coszen_radwt ',mpp_chksum(bnd_type%coszen_radwt) write(outunit,100) 'atm_ice_bnd_type%p ',mpp_chksum(bnd_type%p) ! write(outunit,100) 'atm_ice_bnd_type%data ',mpp_chksum(bnd_type%data) 100 FORMAT("CHECKSUM::",A32," = ",Z20) end subroutine atm_ice_bnd_type_chksum subroutine lnd_ice_bnd_type_chksum(id, timestep, bnd_type) use fms_mod, only: stdout use mpp_mod, only: mpp_chksum character(len=*), intent(in) :: id integer , intent(in) :: timestep type(land_ice_boundary_type), intent(in) :: bnd_type integer :: n, outunit outunit = stdout() write(outunit,*) 'BEGIN CHECKSUM(land_ice_boundary_type):: ', id, timestep write(outunit,100) 'lnd_ice_bnd_type%runoff ',mpp_chksum(bnd_type%runoff) write(outunit,100) 'lnd_ice_bnd_type%calving ',mpp_chksum(bnd_type%calving) ! write(outunit,100) 'lnd_ice_bnd_type%data ',mpp_chksum(bnd_type%data) 100 FORMAT("CHECKSUM::",A32," = ",Z20) end subroutine lnd_ice_bnd_type_chksum end module ice_model_mod