From 84e7816466f53054c5418e4af7c7d23ac8f2da29 Mon Sep 17 00:00:00 2001 From: Darin Scott Comeau Date: Thu, 26 Apr 2018 11:07:01 -0600 Subject: [PATCH 1/5] Incremental commit for capability for icebergs to induce sea ice ridging. Will error out unless config_berg_ridge_seaice = false. --- src/core_seaice/Registry.xml | 1 + src/core_seaice/column/Makefile | 9 +- src/core_seaice/column/ice_bergs_mechred.F90 | 571 ++++++++++++++++++ src/core_seaice/column/ice_mechred.F90 | 2 +- src/core_seaice/shared/mpas_seaice_bergs.F | 19 +- src/core_seaice/shared/mpas_seaice_column.F | 316 +++++++++- .../shared/mpas_seaice_time_integration.F | 2 +- .../configurations/icebergs/namelist.seaice | 1 + .../standard_physics/namelist.seaice | 1 + 9 files changed, 914 insertions(+), 8 deletions(-) create mode 100644 src/core_seaice/column/ice_bergs_mechred.F90 diff --git a/src/core_seaice/Registry.xml b/src/core_seaice/Registry.xml index ac2ec46e5a..431b9d68b4 100644 --- a/src/core_seaice/Registry.xml +++ b/src/core_seaice/Registry.xml @@ -1685,6 +1685,7 @@ + diff --git a/src/core_seaice/column/Makefile b/src/core_seaice/column/Makefile index 822ecc1305..ca9b303e2f 100644 --- a/src/core_seaice/column/Makefile +++ b/src/core_seaice/column/Makefile @@ -28,7 +28,8 @@ OBJS = ice_colpkg.o \ ice_brine.o \ ice_algae.o \ ice_zbgc.o \ - ice_zbgc_shared.o + ice_zbgc_shared.o \ + ice_bergs_mechred.o all: $(OBJS) @@ -40,7 +41,7 @@ ice_constants_colpkg.F90: cp constants/cice/ice_constants_colpkg.F90 . endif -ice_colpkg.o: ice_kinds_mod.o ice_constants_colpkg.o ice_warnings.o ice_colpkg_shared.o ice_therm_shared.o ice_orbital.o ice_atmo.o ice_age.o ice_firstyear.o ice_flux_colpkg.o ice_meltpond_cesm.o ice_meltpond_lvl.o ice_meltpond_topo.o ice_therm_vertical.o ice_itd.o ice_therm_itd.o ice_shortwave.o ice_mechred.o ice_colpkg_tracers.o ice_atmo.o ice_mushy_physics.o ice_zbgc.o ice_zbgc_shared.o ice_aerosol.o ice_algae.o ice_brine.o ice_zsalinity.o +ice_colpkg.o: ice_kinds_mod.o ice_constants_colpkg.o ice_warnings.o ice_colpkg_shared.o ice_therm_shared.o ice_orbital.o ice_atmo.o ice_age.o ice_firstyear.o ice_flux_colpkg.o ice_meltpond_cesm.o ice_meltpond_lvl.o ice_meltpond_topo.o ice_therm_vertical.o ice_itd.o ice_therm_itd.o ice_shortwave.o ice_mechred.o ice_colpkg_tracers.o ice_atmo.o ice_mushy_physics.o ice_zbgc.o ice_zbgc_shared.o ice_aerosol.o ice_algae.o ice_brine.o ice_zsalinity.o ice_bergs_mechred.o ice_kinds_mod.o: @@ -96,7 +97,9 @@ ice_brine.o: ice_kinds_mod.o ice_constants_colpkg.o ice_colpkg_shared.o ice_colp ice_zbgc.o: ice_kinds_mod.o ice_constants_colpkg.o ice_colpkg_shared.o ice_colpkg_tracers.o ice_zbgc_shared.o ice_therm_shared.o ice_itd.o -ice_zbgc_shared.o: ice_kinds_mod.o ice_constants_colpkg.o ice_colpkg_shared.o +ice_zbgc_shared.o: ice_kinds_mod.o ice_constants_colpkg.o ice_colpkg_shared.o + +ice_bergs_mechred.o: ice_kinds_mod.o ice_constants_colpkg.o ice_colpkg_shared.o ice_mechred.o ice_itd.o ice_colpkg_tracers.o .F90.o: diff --git a/src/core_seaice/column/ice_bergs_mechred.F90 b/src/core_seaice/column/ice_bergs_mechred.F90 new file mode 100644 index 0000000000..050b679229 --- /dev/null +++ b/src/core_seaice/column/ice_bergs_mechred.F90 @@ -0,0 +1,571 @@ +! SVN:$Id: $ +!======================================================================= + +! Ridging of sea ice by advancing ice shelves or icebergs +! based on ice_mechred.F90 +! +! authors: Elizabeth C. Hunke, LANL + + module ice_bergs_mechred + + use ice_kinds_mod + use ice_constants_colpkg, only: c0, c1, c2, c10, p5, & + puny, Lfresh, rhoi, rhos + use ice_itd, only: column_sum, & + column_conservation_check + use ice_mechred, only: asum_ridging, ridge_itd, ridge_shift, ridge_prep + use ice_colpkg_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice + use ice_warnings, only: add_warning + + implicit none + save + + private + public :: ridge_ice_by_bergs + + logical (kind=log_kind), parameter :: & + l_conservation_check = .false. ! if true, check conservation +! l_conservation_check = .true. ! if true, check conservation + +!======================================================================= + + contains + +!======================================================================= + +! Compute changes in the sea ice thickness distribution due to ice bergs. +! Based on ice_mechred.F90 +! +! NOTE: This subroutine operates over a single block. +! +! !REVISION HISTORY: +! +! author: William H. Lipscomb, LANL +! 2008: Elizabeth Hunke modified for icebergs. + + subroutine ridge_ice_by_bergs (dt, dice, & + ncat, n_aero, & + nilyr, nslyr, & + ntrcr, hin_max, & + aicen, btrcrn, & + vicen, vsnon, & + aice0, bice, & + trcr_depend, trcr_base , & + n_trcr_strata, & + nt_strata, l_stop, & + stop_label, & + krdg_partic, krdg_redist, & + mu_rdg, tr_brine, & + dardg1dt, dardg2dt, & + dvirdgdt, opening, & + fpond, & + fresh, fhocn, & + faero_ocn, & + aparticn, krdgn, & + aredistn, vredistn, & + dardg1ndt, dardg2ndt, & + dvirdgndt, & + araftn, vraftn) + + ! use ice_colpkg_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice + +! +! !INPUT/OUTPUT PARAMETERS: +! + real (kind=dbl_kind), intent(in) :: & + dt , & ! time step + dice , & ! fractional area to be ridged (displaced) + bice , & ! fractional berg (and shelf) area + mu_rdg ! gives e-folding scale of ridged ice (m^.5) + + integer (kind=int_kind), intent(in) :: & + ncat , & ! number of thickness categories + n_aero , & ! number of aerosol tracers + nilyr , & ! number of ice layers + nslyr , & ! number of snow layers + ! nblyr , & ! number of snow layers + ntrcr ! number of tracers in use + + real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: & + hin_max ! category limits (m) + + real (kind=dbl_kind), dimension (ncat), intent(inout) :: & + aicen , & ! concentration of ice + vicen , & ! volume per unit area of ice (m) + vsnon ! volume per unit area of snow (m) + + real (kind=dbl_kind), dimension (:,:), intent(inout) :: & + btrcrn ! temporary array for tracers + + real (kind=dbl_kind), intent(inout) :: & + aice0 ! concentration of open water + + integer (kind=int_kind), dimension (:), intent(in) :: & + trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon + n_trcr_strata ! number of underlying tracer layers + + real (kind=dbl_kind), dimension (:,:), intent(in) :: & + trcr_base ! = 0 or 1 depending on tracer dependency + ! argument 2: (1) aice, (2) vice, (3) vsno + + integer (kind=int_kind), dimension (:,:), intent(in) :: & + nt_strata ! indices of underlying tracer layers + + logical (kind=log_kind), intent(out) :: & + l_stop ! if true, abort on return + + character (char_len), intent(out) :: & + stop_label ! diagnostic information for abort + + integer (kind=int_kind), intent(in) :: & + krdg_partic , & ! selects participation function + krdg_redist ! selects redistribution function + + logical (kind=log_kind), intent(in) :: & + tr_brine ! if .true., brine height differs from ice thickness + + ! optional history fields + real (kind=dbl_kind), intent(inout), optional :: & + dardg1dt , & ! rate of fractional area loss by ridging ice (1/s) + dardg2dt , & ! rate of fractional area gain by new ridges (1/s) + dvirdgdt , & ! rate of ice volume ridged (m/s) + opening , & ! rate of opening due to divergence/shear (1/s) + fpond , & ! fresh water flux to ponds (kg/m^2/s) + fresh , & ! fresh water flux to ocean (kg/m^2/s) + fhocn ! net heat flux to ocean (W/m^2) + + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & + dardg1ndt , & ! rate of fractional area loss by ridging ice (1/s) + dardg2ndt , & ! rate of fractional area gain by new ridges (1/s) + dvirdgndt , & ! rate of ice volume ridged (m/s) + aparticn , & ! participation function + krdgn , & ! mean ridge thickness/thickness of ridging ice + araftn , & ! rafting ice area + vraftn , & ! rafting ice volume + aredistn , & ! redistribution function: fraction of new ridge area + vredistn ! redistribution function: fraction of new ridge volume + + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & + faero_ocn ! aerosol flux to ocean (kg/m^2/s) + + ! local variables + + real (kind=dbl_kind), dimension (ncat) :: & + eicen , & ! energy of melting for each ice layer (J/m^2) + esnon , & ! energy of melting for each snow layer (J/m^2) + vbrin , & ! ice volume with defined by brine height (m) + sicen ! Bulk salt in h ice (ppt*m) + + + ! variables for ridging routines + real (kind=dbl_kind) :: & + asum , & ! sum of ice and open water area + aksum , & ! ratio of area removed to area ridged + msnow_mlt , & ! mass of snow added to ocean (kg m-2) + esnow_mlt , & ! energy needed to melt snow in ocean (J m-2) + mpond , & ! mass of pond added to ocean (kg m-2) + closing_net, & ! net rate at which area is removed (1/s) + ! (ridging ice area - area of new ridges) / dt + divu_adv , & ! divu as implied by transport scheme (1/s) + opning , & ! rate of opening due to divergence/shear + ! opning is a local variable; + ! opening is the history diagnostic variable + ardg1 , & ! fractional area loss by ridging ice + ardg2 , & ! fractional area gain by new ridges + virdg , & ! ice volume ridged + aopen ! area opening due to divergence/shear + + real (kind=dbl_kind), dimension (n_aero) :: & + maero ! aerosol mass added to ocean (kg m-2) + + real (kind=dbl_kind), dimension (0:ncat) :: & + apartic ! participation function; fraction of ridging + ! and closing associated w/ category n + + real (kind=dbl_kind), dimension (ncat) :: & + hrmin , & ! minimum ridge thickness + hrmax , & ! maximum ridge thickness (krdg_redist = 0) + hrexp , & ! ridge e-folding thickness (krdg_redist = 1) + krdg , & ! mean ridge thickness/thickness of ridging ice + ardg1n , & ! area of ice ridged + ardg2n , & ! area of new ridges + virdgn , & ! ridging ice volume + mraftn ! rafting ice mask + + real (kind=dbl_kind) :: & + vice_init, vice_final, & ! ice volume summed over categories + vsno_init, vsno_final, & ! snow volume summed over categories + eice_init, eice_final, & ! ice energy summed over layers + vbri_init, vbri_final, & ! ice volume in fbri*vicen summed over categories + sice_init ,sice_final, & ! ice bulk salinity summed over categories + esno_init, esno_final ! snow energy summed over layers + + integer (kind=int_kind), parameter :: & + nitermax = 20 ! max number of ridging iterations + + integer (kind=int_kind) :: & + k , & ! layer index + n , & ! thickness category index + niter , & ! iteration counter + it ! tracer index + + real (kind=dbl_kind) :: & + dti ! 1 / dt + + logical (kind=log_kind) :: & + iterate_ridging, & ! if true, repeat the ridging + asum_error ! flag for asum .ne. 1 + + character (len=char_len) :: & + fieldid ! field identifier + + character(len=char_len_long) :: & + warning ! warning message + + !----------------------------------------------------------------- + ! Initialize + !----------------------------------------------------------------- + + l_stop = .false. + + msnow_mlt = c0 + esnow_mlt = c0 + maero(:) = c0 + ardg1 = c0 + ardg2 = c0 + virdg = c0 + ardg1n(:) = c0 + ardg2n(:) = c0 + virdgn(:) = c0 + mpond = c0 + aopen = c0 + + divu_adv = c0 + closing_net = c0 + opning = c0 + asum = c0 + + ! !----------------------------------------------------------------- + ! ! Compute area of ice plus open water before ridging. + ! !----------------------------------------------------------------- + + ! call asum_ridging (ncat, aicen, aice0, asum) + ! ! dcmod - adjust with berg area + ! asum = asum + bice + + ! !----------------------------------------------------------------- + ! ! Compute the area opening and closing. + ! !----------------------------------------------------------------- + + ! call ridge_prep (dt, & + ! ncat, hin_max, & + ! rdg_conv, rdg_shear, & + ! asum, closing_net, & + ! divu_adv, opning) + + !----------------------------------------------------------------- + ! Compute initial values of conserved quantities. + !----------------------------------------------------------------- +!dcmod +aice0 = aice0 + bice + if (l_conservation_check) then + + do n = 1, ncat + eicen(n) = c0 + esnon(n) = c0 + sicen(n) = c0 + vbrin(n) = c0 + + do k = 1, nilyr + eicen(n) = eicen(n) + btrcrn(nt_qice+k-1,n) & + * vicen(n)/real(nilyr,kind=dbl_kind) + enddo + + do k = 1, nslyr + esnon(n) = esnon(n) + btrcrn(nt_qsno+k-1,n) & + * vsnon(n)/real(nslyr,kind=dbl_kind) + enddo ! k + + vbrin(n) = vicen(n) + if (tr_brine) vbrin(n) = btrcrn(nt_fbri,n) * vicen(n) + + enddo ! n + + call column_sum (ncat, vicen(:), vice_init) + call column_sum (ncat, vsnon(:), vsno_init) + call column_sum (ncat, eicen(:), eice_init) + call column_sum (ncat, esnon(:), esno_init) + call column_sum (ncat, sicen(:), sice_init) + call column_sum (ncat, vbrin(:), vbri_init) + + endif ! conservation check + + do niter = 1, nitermax + + !----------------------------------------------------------------- + ! Compute the thickness distribution of ridging ice + ! and various quantities associated with the new ridged ice. + !----------------------------------------------------------------- +! dcmod - need to add bice to open water here? + call ridge_itd (ncat, aice0, & + aicen(:), vicen(:), & + krdg_partic, krdg_redist, & + mu_rdg, & + aksum, apartic(:), & + hrmin(:), hrmax(:), & + hrexp(:), krdg(:), & + aparticn(:), krdgn(:), & + mraftn(:)) + + !----------------------------------------------------------------- + ! Compute the area opening and closing. + !----------------------------------------------------------------- +!dcmod - sure sign is correct? + closing_net = dice/dt ! net closing rate + divu_adv = -closing_net + opning = c0 + + !----------------------------------------------------------------- + ! Redistribute area, volume, and energy. + !----------------------------------------------------------------- + + call ridge_shift (ntrcr, dt, & + ncat, hin_max(:), & + aicen(:), btrcrn(:,:), & + vicen(:), vsnon(:), & + aice0, trcr_depend, & + trcr_base, n_trcr_strata, & + nt_strata, krdg_redist, & + aksum, apartic(:), & + hrmin(:), hrmax(:), & + hrexp(:), krdg(:), & + closing_net, opning, & + ardg1, ardg2, & + virdg, aopen, & + ardg1n(:), ardg2n(:), & + virdgn(:), & + nslyr, n_aero, & + msnow_mlt, esnow_mlt, & + maero(:), mpond, & + l_stop, stop_label, & + aredistn(:), vredistn(:)) + if (l_stop) return + + !----------------------------------------------------------------- + ! Make sure the new area = 1. If not (because the closing + ! and opening rates were reduced above), prepare to ridge again + ! with new rates. + !----------------------------------------------------------------- + + call asum_ridging (ncat, aicen(:), aice0, asum) + ! dcmod - adjust with berg area + ! asum = asum + bice + if (abs(asum - c1) < puny) then + iterate_ridging = .false. + closing_net = c0 + opning = c0 + else + iterate_ridging = .true. + divu_adv = (c1 - asum) / dt + closing_net = max(c0, -divu_adv) + opning = max(c0, divu_adv) + endif + + !----------------------------------------------------------------- + ! If done, exit. If not, prepare to ridge again. + !----------------------------------------------------------------- + + if (iterate_ridging) then + write(warning,*) 'Repeat ridging (bergs), niter =', niter + write(warning,*) 'asum =', asum + write(warning,*) 'aicen =', asum - aice0 + write(warning,*) 'aice0 =', aice0 + call add_warning(warning) + else + ! exit rdg_iteration + exit + endif + + if (niter == nitermax) then + write(warning,*) ' ' + call add_warning(warning) + write(warning,*) 'Exceeded max number of ridging iterations' + call add_warning(warning) + write(warning,*) 'max =',nitermax + call add_warning(warning) + l_stop = .true. + stop_label = "ridge_ice: Exceeded max number of ridging iterations" + return + endif + + enddo ! niter + + !----------------------------------------------------------------- + ! Compute final values of conserved quantities. + ! Check for conservation (allowing for snow thrown into ocean). + !----------------------------------------------------------------- + + if (l_conservation_check) then + + do n = 1, ncat + eicen(n) = c0 + esnon(n) = c0 + sicen(n) = c0 + vbrin(n) = c0 + + do k = 1, nilyr + eicen(n) = eicen(n) + btrcrn(nt_qice+k-1,n) & + * vicen(n)/real(nilyr,kind=dbl_kind) + enddo + do k = 1, nslyr + esnon(n) = esnon(n) + btrcrn(nt_qsno+k-1,n) & + * vsnon(n)/real(nslyr,kind=dbl_kind) + enddo + + vbrin(n) = vicen(n) + if (tr_brine) vbrin(n) = btrcrn(nt_fbri,n) * vbrin(n) + + enddo ! ncat + + call column_sum (ncat, vicen(:), vice_final) + call column_sum (ncat, vsnon(:), vsno_final) + call column_sum (ncat, eicen(:), eice_final) + call column_sum (ncat, esnon(:), esno_final) + call column_sum (ncat, sicen(:), sice_final) + call column_sum (ncat, vbrin(:), vbri_final) + + vsno_final = vsno_final + msnow_mlt/rhos + esno_final = esno_final + esnow_mlt + + fieldid = 'vice, ridging' + call column_conservation_check (fieldid, & + vice_init, vice_final, & + puny, & + l_stop) + if (l_stop) return + + fieldid = 'vsno, ridging' + call column_conservation_check (fieldid, & + vsno_init, vsno_final, & + puny, & + l_stop) + if (l_stop) return + + fieldid = 'eice, ridging' + call column_conservation_check (fieldid, & + eice_init, eice_final, & + puny*Lfresh*rhoi, & + l_stop) + if (l_stop) return + + fieldid = 'esno, ridging' + call column_conservation_check (fieldid, & + esno_init, esno_final, & + puny*Lfresh*rhos, & + l_stop) + if (l_stop) return + + endif ! l_conservation_check + + !----------------------------------------------------------------- + ! Compute ridging diagnostics. + !----------------------------------------------------------------- + + dti = c1/dt + + if (present(dardg1dt)) then + dardg1dt = ardg1*dti + endif + if (present(dardg2dt)) then + dardg2dt = ardg2*dti + endif + if (present(dvirdgdt)) then + dvirdgdt = virdg*dti + endif + if (present(opening)) then + opening = aopen*dti + endif + + if (present(dardg1ndt)) then + do n = 1, ncat + dardg1ndt(n) = ardg1n(n)*dti + enddo + endif + if (present(dardg2ndt)) then + do n = 1, ncat + dardg2ndt(n) = ardg2n(n)*dti + enddo + endif + if (present(dvirdgndt)) then + do n = 1, ncat + dvirdgndt(n) = virdgn(n)*dti + enddo + endif + if (present(araftn)) then + do n = 1, ncat + araftn(n) = mraftn(n)*ardg2n(n) +! araftn(n) = mraftn(n)*ardg1n(n)*p5 + enddo + endif + if (present(vraftn)) then + do n = 1, ncat + vraftn(n) = mraftn(n)*virdgn(n) + enddo + endif + + !----------------------------------------------------------------- + ! Update fresh water and heat fluxes due to snow melt. + !----------------------------------------------------------------- + + if (present(fresh)) then + fresh = fresh + msnow_mlt*dti + endif + if (present(fhocn)) then + fhocn = fhocn + esnow_mlt*dti + endif + if (present(faero_ocn)) then + do it = 1, n_aero + faero_ocn(it) = faero_ocn(it) + maero(it)*dti + enddo + endif + if (present(fpond)) then + fpond = fpond - mpond ! units change later + endif + + !----------------------------------------------------------------- + ! Check for fractional ice area > 1. + !----------------------------------------------------------------- + + call asum_ridging (ncat, aicen(:), aice0, asum) + ! dcmod - adjust with berg area + ! asum = asum + bice + if (abs(asum - c1) > puny) then + l_stop = .true. + + write(warning,*) ' ' + call add_warning(warning) + write(warning,*) 'Ridging error (bergs): total area > 1' + call add_warning(warning) + write(warning,*) 'area:', asum + call add_warning(warning) + write(warning,*) 'n, aicen:' + call add_warning(warning) + write(warning,*) 0, aice0 + call add_warning(warning) + do n = 1, ncat + write(warning,*) n, aicen(n) + call add_warning(warning) + enddo + return + endif + +!dcmod +aice0 = aice0 - bice + + end subroutine ridge_ice_by_bergs + +!======================================================================= + + end module ice_bergs_mechred + +!======================================================================= diff --git a/src/core_seaice/column/ice_mechred.F90 b/src/core_seaice/column/ice_mechred.F90 index 668621fee7..daa59f1b02 100644 --- a/src/core_seaice/column/ice_mechred.F90 +++ b/src/core_seaice/column/ice_mechred.F90 @@ -46,7 +46,7 @@ module ice_mechred save private - public :: ridge_ice, asum_ridging, ridge_itd + public :: ridge_ice, asum_ridging, ridge_itd, ridge_shift, ridge_prep real (kind=dbl_kind), parameter :: & Cs = p25 , & ! fraction of shear energy contrbtng to ridging diff --git a/src/core_seaice/shared/mpas_seaice_bergs.F b/src/core_seaice/shared/mpas_seaice_bergs.F index 5c34822279..c4f6dce425 100644 --- a/src/core_seaice/shared/mpas_seaice_bergs.F +++ b/src/core_seaice/shared/mpas_seaice_bergs.F @@ -211,7 +211,7 @@ end subroutine seaice_run_berg_dynamics ! !----------------------------------------------------------------------- - subroutine seaice_run_berg_postdynamics(domain, clock, itimestep) + subroutine seaice_run_berg_postdynamics(domain, clock, itimestep, ierr) use seaice_berg_decay, only: & seaice_run_berg_decay @@ -225,6 +225,9 @@ subroutine seaice_run_berg_postdynamics(domain, clock, itimestep) use seaice_berg_velocity_solver, only: & seaice_berg_forcing_for_ice + use seaice_column, only: & + seaice_column_ridging_by_bergs + type(domain_type), intent(inout) :: & domain !< Input/Output: @@ -234,14 +237,19 @@ subroutine seaice_run_berg_postdynamics(domain, clock, itimestep) integer, intent(in) :: & itimestep !< Input: + integer, intent(inout) :: & + ierr !< Input/Output: + logical, pointer :: & config_use_berg_decay, & ! flag to turn on / off berg decay parameterization config_berg_update_open_water, & ! flag to turn on / off bergs affecting open water - config_berg_seaice_interaction ! flag to turn on / off bergs affecting open water + config_berg_seaice_interaction, &! flag to turn on / off bergs term in sea ice momentum equation + config_berg_ridge_seaice ! flag to turn on / off bergs causing sea ice to ridge call MPAS_pool_get_config(domain % configs, "config_use_berg_decay", config_use_berg_decay) call MPAS_pool_get_config(domain % configs, "config_berg_update_open_water", config_berg_update_open_water) call MPAS_pool_get_config(domain % configs, "config_berg_seaice_interaction", config_berg_seaice_interaction) + call MPAS_pool_get_config(domain % configs, "config_berg_ridge_seaice", config_berg_ridge_seaice) ! berg mass and area state variables update call mpas_timer_start("Berg update mass") @@ -265,7 +273,14 @@ subroutine seaice_run_berg_postdynamics(domain, clock, itimestep) call mpas_timer_stop("Berg sea ice interaction") ! update open water + call mpas_timer_start("Berg update open water") if (config_berg_update_open_water) call seaice_berg_update_open_water_area(domain) + call mpas_timer_stop("Berg update open water") + + ! berg sea ice interaction + call mpas_timer_start("Berg induced ridging") + if (config_berg_ridge_seaice) call seaice_column_ridging_by_bergs(domain, ierr) + call mpas_timer_stop("Berg induced ridging") end subroutine seaice_run_berg_postdynamics diff --git a/src/core_seaice/shared/mpas_seaice_column.F b/src/core_seaice/shared/mpas_seaice_column.F index 064ae259c2..1c31ff9794 100644 --- a/src/core_seaice/shared/mpas_seaice_column.F +++ b/src/core_seaice/shared/mpas_seaice_column.F @@ -42,7 +42,8 @@ module seaice_column seaice_column_reinitialize_diagnostics_dynamics, & seaice_column_reinitialize_diagnostics_bgc, & seaice_column_coupling_prep, & - seaice_column_finalize + seaice_column_finalize, & + seaice_column_ridging_by_bergs ! tracer object type, private :: ciceTracerObjectType @@ -13265,6 +13266,319 @@ subroutine column_combine_snow_ice_tracers(domain) end subroutine column_combine_snow_ice_tracers +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! seaice_column_ridging_by_bergs +! +!> \brief +!> \author Darin Comeau, LANL +!> \date 9th April 2018 +!> \details +!> +! +!----------------------------------------------------------------------- + + subroutine seaice_column_ridging_by_bergs(domain, ierr) + + use ice_colpkg, only: & + ! colpkg_step_ridge, & + colpkg_get_warnings, & + colpkg_clear_warnings + + use ice_bergs_mechred, only: & + ridge_ice_by_bergs + + type(domain_type), intent(inout) :: domain + + integer, intent(inout) :: ierr + + type(block_type), pointer :: block + + type(MPAS_pool_type), pointer :: & + mesh, & + icestate, & + tracers, & + tracers_aggregate, & + ponds, & + ocean_fluxes, & + ridging, & + aerosols, & + biogeochemistry, & + initial, & + velocity_solver, & + berg_state, & + berg_tracers_aggregate + + ! configs + logical, pointer :: & + config_use_column_biogeochemistry, & + config_use_brine + + real(kind=RKIND), pointer :: & + config_dt + + ! integer, pointer :: & + ! config_dynamics_subcycle_number + + character(len=strKIND), pointer :: & + config_ridging_participation_function, & + config_ridging_redistribution_function + + real(kind=RKIND), pointer :: & + config_ridiging_efolding_scale + + ! dimensions + integer, pointer :: & + nCellsSolve, & + nCategories, & + nIceLayers, & + nSnowLayers, & + nAerosols, & + nBioLayers + + ! variables + real(kind=RKIND), dimension(:), pointer :: & + pondFreshWaterFlux, & + oceanFreshWaterFlux, & + oceanSaltFlux, & + oceanHeatFlux, & + iceAreaCell, & + ridgeConvergence, & + ridgeShear, & + openWaterArea, & + areaLossRidge, & + areaGainRidge, & + iceVolumeRidged, & + openingRateRidge, & + categoryThicknessLimits, & + zSalinityFlux, & + bergDisplacedArea, & + bergAreaCell + + real(kind=RKIND), dimension(:,:), pointer :: & + oceanAerosolFlux, & + ridgeParticipationFunction, & + ratioRidgeThicknessToIce, & + fractionNewRidgeArea, & + fractionNewRidgeVolume, & + areaLossRidgeCategory, & + areaGainRidgeCategory, & + iceVolumeRidgedCategory, & + raftingIceArea, & + raftingIceVolume, & + oceanBioFluxes + + real(kind=RKIND), dimension(:,:,:), pointer :: & + iceAreaCategory, & + iceVolumeCategory, & + snowVolumeCategory + + ! integer, dimension(:,:), pointer :: & + ! newlyFormedIce + + integer, dimension(:), pointer :: & + indexToCellID + + real(kind=RKIND), pointer :: & + dynamicsTimeStep + + ! local + integer :: & + iCell, & + iCategory + + ! logical, dimension(:), allocatable :: & + ! newlyFormedIceLogical + + ! real(kind=RKIND), dimension(:,:) :: & + ! bergTracerArray + + logical :: & + abortFlag, & + setGetPhysicsTracers, & + setGetBGCTracers + + character(len=strKIND) :: & + abortMessage, & + abortLocation + + character(len=strKINDWarnings), dimension(:), allocatable :: & + warnings + + block => domain % blocklist + do while (associated(block)) + + call MPAS_pool_get_subpool(block % structs, "mesh", mesh) + call MPAS_pool_get_subpool(block % structs, "tracers", tracers) + call MPAS_pool_get_subpool(block % structs, "tracers_aggregate", tracers_aggregate) + call MPAS_pool_get_subpool(block % structs, "icestate", icestate) + call MPAS_pool_get_subpool(block % structs, "ponds", ponds) + call MPAS_pool_get_subpool(block % structs, "ocean_fluxes", ocean_fluxes) + call MPAS_pool_get_subpool(block % structs, "ridging", ridging) + call MPAS_pool_get_subpool(block % structs, "aerosols", aerosols) + call MPAS_pool_get_subpool(block % structs, "biogeochemistry", biogeochemistry) + call MPAS_pool_get_subpool(block % structs, "initial", initial) + call MPAS_pool_get_subpool(block % structs, "velocity_solver", velocity_solver) + call MPAS_pool_get_subpool(block % structs, "berg_state", berg_state) + call MPAS_pool_get_subpool(block % structs, "berg_tracers_aggregate", berg_tracers_aggregate) + + ! call MPAS_pool_get_config(block % configs, "config_dynamics_subcycle_number", config_dynamics_subcycle_number) + call MPAS_pool_get_config(block % configs, "config_use_column_biogeochemistry", config_use_column_biogeochemistry) + call MPAS_pool_get_config(block % configs, "config_use_brine", config_use_brine) + call MPAS_pool_get_config(block % configs, "config_ridging_participation_function", config_ridging_participation_function) + call MPAS_pool_get_config(block % configs, "config_ridging_redistribution_function", config_ridging_redistribution_function) + call MPAS_pool_get_config(block % configs, "config_ridiging_efolding_scale", config_ridiging_efolding_scale) + + call MPAS_pool_get_array(velocity_solver, "dynamicsTimeStep", dynamicsTimeStep) + + call MPAS_pool_get_dimension(mesh, "nCellsSolve", nCellsSolve) + call MPAS_pool_get_dimension(mesh, "nCategories", nCategories) + call MPAS_pool_get_dimension(mesh, "nIceLayers", nIceLayers) + call MPAS_pool_get_dimension(mesh, "nSnowLayers", nSnowLayers) + call MPAS_pool_get_dimension(mesh, "nAerosols", nAerosols) + ! call MPAS_pool_get_dimension(block % dimensions, "nBioLayers", nBioLayers) + + call MPAS_pool_get_array(mesh, "indexToCellID", indexToCellID) + + call MPAS_pool_get_array(tracers_aggregate, "iceAreaCell", iceAreaCell) + + call MPAS_pool_get_array(icestate, "openWaterArea", openWaterArea) + + call MPAS_pool_get_array(tracers, "iceAreaCategory", iceAreaCategory, 1) + call MPAS_pool_get_array(tracers, "iceVolumeCategory", iceVolumeCategory, 1) + call MPAS_pool_get_array(tracers, "snowVolumeCategory", snowVolumeCategory, 1) + + call MPAS_pool_get_array(ocean_fluxes, "oceanFreshWaterFlux", oceanFreshWaterFlux) + call MPAS_pool_get_array(ocean_fluxes, "oceanSaltFlux", oceanSaltFlux) + call MPAS_pool_get_array(ocean_fluxes, "oceanHeatFlux", oceanHeatFlux) + + call MPAS_pool_get_array(ridging, "ridgeConvergence", ridgeConvergence) + call MPAS_pool_get_array(ridging, "ridgeShear", ridgeShear) + call MPAS_pool_get_array(ridging, "areaLossRidge", areaLossRidge) + call MPAS_pool_get_array(ridging, "areaGainRidge", areaGainRidge) + call MPAS_pool_get_array(ridging, "iceVolumeRidged", iceVolumeRidged) + call MPAS_pool_get_array(ridging, "openingRateRidge", openingRateRidge) + call MPAS_pool_get_array(ridging, "ridgeParticipationFunction", ridgeParticipationFunction) + call MPAS_pool_get_array(ridging, "ratioRidgeThicknessToIce", ratioRidgeThicknessToIce) + call MPAS_pool_get_array(ridging, "fractionNewRidgeArea", fractionNewRidgeArea) + call MPAS_pool_get_array(ridging, "fractionNewRidgeVolume", fractionNewRidgeVolume) + call MPAS_pool_get_array(ridging, "areaLossRidgeCategory", areaLossRidgeCategory) + call MPAS_pool_get_array(ridging, "areaGainRidgeCategory", areaGainRidgeCategory) + call MPAS_pool_get_array(ridging, "iceVolumeRidgedCategory", iceVolumeRidgedCategory) + call MPAS_pool_get_array(ridging, "raftingIceArea", raftingIceArea) + call MPAS_pool_get_array(ridging, "raftingIceVolume", raftingIceVolume) + + call MPAS_pool_get_array(aerosols, "oceanAerosolFlux", oceanAerosolFlux) + + call MPAS_pool_get_array(ponds, "pondFreshWaterFlux", pondFreshWaterFlux) + + ! call MPAS_pool_get_array(biogeochemistry, "newlyFormedIce", newlyFormedIce) + ! call MPAS_pool_get_array(biogeochemistry, "oceanBioFluxes", oceanBioFluxes) + ! call MPAS_pool_get_array(biogeochemistry, "zSalinityFlux", zSalinityFlux) + + call MPAS_pool_get_array(initial, "categoryThicknessLimits", categoryThicknessLimits) + + call MPAS_pool_get_array(berg_state, "bergDisplacedArea", bergDisplacedArea) + + call MPAS_pool_get_array(berg_tracers_aggregate, "bergAreaCell", bergAreaCell) + + ! ! newly formed ice + ! allocate(newlyFormedIceLogical(nCategories)) + + setGetPhysicsTracers = .true. + setGetBGCTracers = config_use_column_biogeochemistry + + do iCell = 1, nCellsSolve +!dcmod - not needed, for bgc? + ! newly formed ice + ! do iCategory = 1, nCategories + ! newlyFormedIceLogical(iCategory) = (newlyFormedIce(iCategory,iCell) == 1) + ! enddo ! iCategory + + ! set the category tracer array + call set_cice_tracer_array_category(block, ciceTracerObject, iCell, setGetPhysicsTracers, setGetBGCTracers) + + abortFlag = .false. + abortMessage = "" +!dcmod temporary so as to not update + ! bergTracerArray = ciceTracerObject % tracerArrayCategory + + call colpkg_clear_warnings() + call ridge_ice_by_bergs(& + dynamicsTimeStep, & + bergDisplacedArea(iCell), & + nCategories, & + nAerosols, & + nIceLayers, & + nSnowLayers, & + ciceTracerObject % nTracers, & + categoryThicknessLimits, & + iceAreaCategory(1,:,iCell), & +! bergTracerArray, & ! trcrn !dcmod is there reason this needs to be temporary array? + ciceTracerObject % tracerArrayCategory, & ! trcrn !dcmod is there reason this needs to be temporary array? + iceVolumeCategory(1,:,iCell), & + snowVolumeCategory(1,:,iCell), & + openWaterArea(iCell), & + bergAreaCell(iCell), & + ciceTracerObject % parentIndex, & ! trcr_depend + ciceTracerObject % firstAncestorMask, & ! trcr_base + ciceTracerObject % ancestorNumber, & !n_trcr_strata + ciceTracerObject % ancestorIndices, & !nt_strata + abortFlag, & ! move these at bottom + abortMessage, & + config_cice_int("config_ridging_participation_function", config_ridging_participation_function), & + config_cice_int("config_ridging_redistribution_function", config_ridging_redistribution_function), & + config_ridiging_efolding_scale, & ! dcmod change + config_use_brine, & + areaLossRidge(iCell), & + areaGainRidge(iCell), & + iceVolumeRidged(iCell), & + openingRateRidge(iCell), & + pondFreshWaterFlux(iCell), & + oceanFreshWaterFlux(iCell), & + oceanHeatFlux(iCell), & + oceanAerosolFlux(:,iCell), & + ridgeParticipationFunction(:,iCell), & + ratioRidgeThicknessToIce(:,iCell), & + fractionNewRidgeArea(:,iCell), & + fractionNewRidgeVolume(:,iCell), & + areaLossRidgeCategory(:,iCell), & + areaGainRidgeCategory(:,iCell), & + iceVolumeRidgedCategory(:,iCell), & + raftingIceArea(:,iCell), & + raftingIceVolume(:,iCell)) + + call colpkg_get_warnings(warnings) + call column_write_warnings(warnings) + + ! ! update + ! do iCategory = 1, nCategories + ! newlyFormedIce(iCategory,iCell) = 0 + ! if (newlyFormedIceLogical(iCategory)) newlyFormedIce(iCategory,iCell) = 1 + ! enddo ! iCategory + + ! ! get category tracer array + ! call get_cice_tracer_array_category(block, ciceTracerObject, iCell, setGetPhysicsTracers, setGetBGCTracers) + + ! code abort + if (abortFlag) then + call mpas_log_write("column_ridging_by_bergs: "//trim(abortMessage) , messageType=MPAS_LOG_ERR) + call mpas_log_write("iCell: $i", messageType=MPAS_LOG_ERR, intArgs=(/indexToCellID(iCell)/)) + ierr = SEAICE_ERROR_COL_RIDGING + call seaice_critical_error_write_block(domain, block) + return + endif + + enddo ! iCell + + ! ! newly formed ice + ! deallocate(newlyFormedIceLogical) + + block => block % next + enddo + + end subroutine seaice_column_ridging_by_bergs + !----------------------------------------------------------------------- ! Warning messages !----------------------------------------------------------------------- diff --git a/src/core_seaice/shared/mpas_seaice_time_integration.F b/src/core_seaice/shared/mpas_seaice_time_integration.F index 8270b847c1..107fbbceea 100644 --- a/src/core_seaice/shared/mpas_seaice_time_integration.F +++ b/src/core_seaice/shared/mpas_seaice_time_integration.F @@ -187,7 +187,7 @@ subroutine seaice_timestep(& ! berg decay and calving check if (config_use_bergs) then call mpas_timer_start("Berg post-dynamics") - call seaice_run_berg_postdynamics(domain, clock, itimestep) + call seaice_run_berg_postdynamics(domain, clock, itimestep, ierr) call mpas_timer_stop("Berg post-dynamics") endif diff --git a/testing_and_setup/seaice/configurations/icebergs/namelist.seaice b/testing_and_setup/seaice/configurations/icebergs/namelist.seaice index 89c9a19deb..a9ff16a1b7 100644 --- a/testing_and_setup/seaice/configurations/icebergs/namelist.seaice +++ b/testing_and_setup/seaice/configurations/icebergs/namelist.seaice @@ -347,6 +347,7 @@ config_berg_advection_outside_velocity_subcycling = true config_berg_update_open_water = true config_berg_seaice_interaction = true + config_berg_ridge_seaice = true / &AM_highFrequencyOutput config_AM_highFrequencyOutput_enable = false diff --git a/testing_and_setup/seaice/configurations/standard_physics/namelist.seaice b/testing_and_setup/seaice/configurations/standard_physics/namelist.seaice index 47ad16f388..a7930b2c1b 100644 --- a/testing_and_setup/seaice/configurations/standard_physics/namelist.seaice +++ b/testing_and_setup/seaice/configurations/standard_physics/namelist.seaice @@ -347,6 +347,7 @@ config_berg_advection_outside_velocity_subcycling = false config_berg_update_open_water = false config_berg_seaice_interaction = false + config_berg_ridge_seaice = false / &AM_highFrequencyOutput config_AM_highFrequencyOutput_enable = false From 9cbb2095be711f679bb2dc184ec32172accdc4c5 Mon Sep 17 00:00:00 2001 From: Darin Scott Comeau Date: Tue, 8 May 2018 09:40:01 -0600 Subject: [PATCH 2/5] Cleaning up ice_bergs_mechred.F90; still not working correctly. --- src/core_seaice/column/ice_bergs_mechred.F90 | 102 ++++++++----------- 1 file changed, 41 insertions(+), 61 deletions(-) diff --git a/src/core_seaice/column/ice_bergs_mechred.F90 b/src/core_seaice/column/ice_bergs_mechred.F90 index 050b679229..e10b24e85d 100644 --- a/src/core_seaice/column/ice_bergs_mechred.F90 +++ b/src/core_seaice/column/ice_bergs_mechred.F90 @@ -15,7 +15,7 @@ module ice_bergs_mechred column_conservation_check use ice_mechred, only: asum_ridging, ridge_itd, ridge_shift, ridge_prep use ice_colpkg_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice - use ice_warnings, only: add_warning + use ice_warnings, only: add_warning implicit none save @@ -33,7 +33,7 @@ module ice_bergs_mechred !======================================================================= -! Compute changes in the sea ice thickness distribution due to ice bergs. +! Compute changes in the sea ice thickness distribution due to icebergs. ! Based on ice_mechred.F90 ! ! NOTE: This subroutine operates over a single block. @@ -53,7 +53,7 @@ subroutine ridge_ice_by_bergs (dt, dice, & trcr_depend, trcr_base , & n_trcr_strata, & nt_strata, l_stop, & - stop_label, & + stop_label, & krdg_partic, krdg_redist, & mu_rdg, tr_brine, & dardg1dt, dardg2dt, & @@ -67,8 +67,6 @@ subroutine ridge_ice_by_bergs (dt, dice, & dvirdgndt, & araftn, vraftn) - ! use ice_colpkg_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice - ! ! !INPUT/OUTPUT PARAMETERS: ! @@ -76,18 +74,18 @@ subroutine ridge_ice_by_bergs (dt, dice, & dt , & ! time step dice , & ! fractional area to be ridged (displaced) bice , & ! fractional berg (and shelf) area - mu_rdg ! gives e-folding scale of ridged ice (m^.5) + mu_rdg ! gives e-folding scale of ridged ice (m^.5) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories - n_aero , & ! number of aerosol tracers + n_aero , & ! number of aerosol tracers nilyr , & ! number of ice layers nslyr , & ! number of snow layers - ! nblyr , & ! number of snow layers - ntrcr ! number of tracers in use + ! nblyr , & ! number of bio layers + ntrcr ! number of tracers in use real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: & - hin_max ! category limits (m) + hin_max ! category limits (m) real (kind=dbl_kind), dimension (ncat), intent(inout) :: & aicen , & ! concentration of ice @@ -95,7 +93,7 @@ subroutine ridge_ice_by_bergs (dt, dice, & vsnon ! volume per unit area of snow (m) real (kind=dbl_kind), dimension (:,:), intent(inout) :: & - btrcrn ! temporary array for tracers + btrcrn ! temporary array for tracers real (kind=dbl_kind), intent(inout) :: & aice0 ! concentration of open water @@ -115,14 +113,14 @@ subroutine ridge_ice_by_bergs (dt, dice, & l_stop ! if true, abort on return character (char_len), intent(out) :: & - stop_label ! diagnostic information for abort + stop_label ! diagnostic information for abort integer (kind=int_kind), intent(in) :: & krdg_partic , & ! selects participation function krdg_redist ! selects redistribution function logical (kind=log_kind), intent(in) :: & - tr_brine ! if .true., brine height differs from ice thickness + tr_brine ! if .true., brine height differs from ice thickness ! optional history fields real (kind=dbl_kind), intent(inout), optional :: & @@ -130,7 +128,7 @@ subroutine ridge_ice_by_bergs (dt, dice, & dardg2dt , & ! rate of fractional area gain by new ridges (1/s) dvirdgdt , & ! rate of ice volume ridged (m/s) opening , & ! rate of opening due to divergence/shear (1/s) - fpond , & ! fresh water flux to ponds (kg/m^2/s) + fpond , & ! fresh water flux to ponds (kg/m^2/s) fresh , & ! fresh water flux to ocean (kg/m^2/s) fhocn ! net heat flux to ocean (W/m^2) @@ -146,16 +144,15 @@ subroutine ridge_ice_by_bergs (dt, dice, & vredistn ! redistribution function: fraction of new ridge volume real (kind=dbl_kind), dimension(:), intent(inout), optional :: & - faero_ocn ! aerosol flux to ocean (kg/m^2/s) + faero_ocn ! aerosol flux to ocean (kg/m^2/s) ! local variables real (kind=dbl_kind), dimension (ncat) :: & eicen , & ! energy of melting for each ice layer (J/m^2) - esnon , & ! energy of melting for each snow layer (J/m^2) + esnon , & ! energy of melting for each snow layer (J/m^2) vbrin , & ! ice volume with defined by brine height (m) - sicen ! Bulk salt in h ice (ppt*m) - + sicen ! Bulk salt in h ice (ppt*m) ! variables for ridging routines real (kind=dbl_kind) :: & @@ -163,7 +160,7 @@ subroutine ridge_ice_by_bergs (dt, dice, & aksum , & ! ratio of area removed to area ridged msnow_mlt , & ! mass of snow added to ocean (kg m-2) esnow_mlt , & ! energy needed to melt snow in ocean (J m-2) - mpond , & ! mass of pond added to ocean (kg m-2) + mpond , & ! mass of pond added to ocean (kg m-2) closing_net, & ! net rate at which area is removed (1/s) ! (ridging ice area - area of new ridges) / dt divu_adv , & ! divu as implied by transport scheme (1/s) @@ -173,7 +170,7 @@ subroutine ridge_ice_by_bergs (dt, dice, & ardg1 , & ! fractional area loss by ridging ice ardg2 , & ! fractional area gain by new ridges virdg , & ! ice volume ridged - aopen ! area opening due to divergence/shear + aopen ! area opening due to divergence/shear real (kind=dbl_kind), dimension (n_aero) :: & maero ! aerosol mass added to ocean (kg m-2) @@ -197,7 +194,7 @@ subroutine ridge_ice_by_bergs (dt, dice, & vsno_init, vsno_final, & ! snow volume summed over categories eice_init, eice_final, & ! ice energy summed over layers vbri_init, vbri_final, & ! ice volume in fbri*vicen summed over categories - sice_init ,sice_final, & ! ice bulk salinity summed over categories + sice_init ,sice_final, & ! ice bulk salinity summed over categories esno_init, esno_final ! snow energy summed over layers integer (kind=int_kind), parameter :: & @@ -220,7 +217,11 @@ subroutine ridge_ice_by_bergs (dt, dice, & fieldid ! field identifier character(len=char_len_long) :: & - warning ! warning message + warning ! warning message +! !dcmod +! ! temporary for open water +! real (kind=dbl_kind) :: & +! aice0tmp ! sum of ice and open water area (without bergs) !----------------------------------------------------------------- ! Initialize @@ -245,29 +246,10 @@ subroutine ridge_ice_by_bergs (dt, dice, & opning = c0 asum = c0 - ! !----------------------------------------------------------------- - ! ! Compute area of ice plus open water before ridging. - ! !----------------------------------------------------------------- - - ! call asum_ridging (ncat, aicen, aice0, asum) - ! ! dcmod - adjust with berg area - ! asum = asum + bice - - ! !----------------------------------------------------------------- - ! ! Compute the area opening and closing. - ! !----------------------------------------------------------------- - - ! call ridge_prep (dt, & - ! ncat, hin_max, & - ! rdg_conv, rdg_shear, & - ! asum, closing_net, & - ! divu_adv, opning) - !----------------------------------------------------------------- ! Compute initial values of conserved quantities. !----------------------------------------------------------------- -!dcmod -aice0 = aice0 + bice + if (l_conservation_check) then do n = 1, ncat @@ -296,9 +278,9 @@ subroutine ridge_ice_by_bergs (dt, dice, & call column_sum (ncat, eicen(:), eice_init) call column_sum (ncat, esnon(:), esno_init) call column_sum (ncat, sicen(:), sice_init) - call column_sum (ncat, vbrin(:), vbri_init) + call column_sum (ncat, vbrin(:), vbri_init) - endif ! conservation check + endif ! conservation check do niter = 1, nitermax @@ -306,7 +288,7 @@ subroutine ridge_ice_by_bergs (dt, dice, & ! Compute the thickness distribution of ridging ice ! and various quantities associated with the new ridged ice. !----------------------------------------------------------------- -! dcmod - need to add bice to open water here? + call ridge_itd (ncat, aice0, & aicen(:), vicen(:), & krdg_partic, krdg_redist, & @@ -320,7 +302,7 @@ subroutine ridge_ice_by_bergs (dt, dice, & !----------------------------------------------------------------- ! Compute the area opening and closing. !----------------------------------------------------------------- -!dcmod - sure sign is correct? + closing_net = dice/dt ! net closing rate divu_adv = -closing_net opning = c0 @@ -358,13 +340,15 @@ subroutine ridge_ice_by_bergs (dt, dice, & !----------------------------------------------------------------- call asum_ridging (ncat, aicen(:), aice0, asum) - ! dcmod - adjust with berg area + ! ! dcmod - adjust with berg area - note this will be < 1 ! asum = asum + bice - if (abs(asum - c1) < puny) then + if (abs(asum - c1) < puny) then iterate_ridging = .false. closing_net = c0 opning = c0 else + ! write(warning,*) 'ridging loop (bergs), abs(asum-1):', abs(asum - c1) + ! call add_warning(warning) iterate_ridging = .true. divu_adv = (c1 - asum) / dt closing_net = max(c0, -divu_adv) @@ -377,14 +361,13 @@ subroutine ridge_ice_by_bergs (dt, dice, & if (iterate_ridging) then write(warning,*) 'Repeat ridging (bergs), niter =', niter - write(warning,*) 'asum =', asum - write(warning,*) 'aicen =', asum - aice0 - write(warning,*) 'aice0 =', aice0 + call add_warning(warning) + write(warning,*) 'area (inc. bergs):', asum + bice call add_warning(warning) else ! exit rdg_iteration exit - endif + endif if (niter == nitermax) then write(warning,*) ' ' @@ -411,7 +394,7 @@ subroutine ridge_ice_by_bergs (dt, dice, & eicen(n) = c0 esnon(n) = c0 sicen(n) = c0 - vbrin(n) = c0 + vbrin(n) = c0 do k = 1, nilyr eicen(n) = eicen(n) + btrcrn(nt_qice+k-1,n) & @@ -464,7 +447,7 @@ subroutine ridge_ice_by_bergs (dt, dice, & puny*Lfresh*rhos, & l_stop) if (l_stop) return - + endif ! l_conservation_check !----------------------------------------------------------------- @@ -511,7 +494,7 @@ subroutine ridge_ice_by_bergs (dt, dice, & do n = 1, ncat vraftn(n) = mraftn(n)*virdgn(n) enddo - endif + endif !----------------------------------------------------------------- ! Update fresh water and heat fluxes due to snow melt. @@ -530,7 +513,7 @@ subroutine ridge_ice_by_bergs (dt, dice, & endif if (present(fpond)) then fpond = fpond - mpond ! units change later - endif + endif !----------------------------------------------------------------- ! Check for fractional ice area > 1. @@ -538,8 +521,8 @@ subroutine ridge_ice_by_bergs (dt, dice, & call asum_ridging (ncat, aicen(:), aice0, asum) ! dcmod - adjust with berg area - ! asum = asum + bice - if (abs(asum - c1) > puny) then + asum = asum + bice + if (abs(asum - c1) > puny) then l_stop = .true. write(warning,*) ' ' @@ -559,9 +542,6 @@ subroutine ridge_ice_by_bergs (dt, dice, & return endif -!dcmod -aice0 = aice0 - bice - end subroutine ridge_ice_by_bergs !======================================================================= From d60156f106fd169b9695b807cdc135491cef5fc2 Mon Sep 17 00:00:00 2001 From: Darin Scott Comeau Date: Tue, 8 May 2018 12:01:20 -0600 Subject: [PATCH 3/5] Adding fix to iceberg induced ridging; open water area is adjusted inside ridge_ice_by_bergs. --- src/core_seaice/column/ice_bergs_mechred.F90 | 27 ++++++++++++------- .../shared/mpas_seaice_berg_state.F | 2 +- 2 files changed, 18 insertions(+), 11 deletions(-) diff --git a/src/core_seaice/column/ice_bergs_mechred.F90 b/src/core_seaice/column/ice_bergs_mechred.F90 index e10b24e85d..ece890e7f3 100644 --- a/src/core_seaice/column/ice_bergs_mechred.F90 +++ b/src/core_seaice/column/ice_bergs_mechred.F90 @@ -218,10 +218,6 @@ subroutine ridge_ice_by_bergs (dt, dice, & character(len=char_len_long) :: & warning ! warning message -! !dcmod -! ! temporary for open water -! real (kind=dbl_kind) :: & -! aice0tmp ! sum of ice and open water area (without bergs) !----------------------------------------------------------------- ! Initialize @@ -339,16 +335,22 @@ subroutine ridge_ice_by_bergs (dt, dice, & ! with new rates. !----------------------------------------------------------------- + ! dcmod - adjust with berg area call asum_ridging (ncat, aicen(:), aice0, asum) - ! ! dcmod - adjust with berg area - note this will be < 1 - ! asum = asum + bice + asum = asum + bice + + ! dcmod - if less than 1 due to forced ridging, add missing area back to open water + if (asum < c1) aice0 = aice0 + (c1 - asum) + + ! check sum again + call asum_ridging (ncat, aicen(:), aice0, asum) + asum = asum + bice + if (abs(asum - c1) < puny) then iterate_ridging = .false. closing_net = c0 opning = c0 else - ! write(warning,*) 'ridging loop (bergs), abs(asum-1):', abs(asum - c1) - ! call add_warning(warning) iterate_ridging = .true. divu_adv = (c1 - asum) / dt closing_net = max(c0, -divu_adv) @@ -362,8 +364,13 @@ subroutine ridge_ice_by_bergs (dt, dice, & if (iterate_ridging) then write(warning,*) 'Repeat ridging (bergs), niter =', niter call add_warning(warning) - write(warning,*) 'area (inc. bergs):', asum + bice - call add_warning(warning) + ! ! dcmod debugging messages + ! write(warning,*) 'area (inc. bergs):', asum + ! call add_warning(warning) + ! write(warning,*) 'ice:', asum - bice - aice0 + ! call add_warning(warning) + ! write(warning,*) 'open water:', aice0 + ! call add_warning(warning) else ! exit rdg_iteration exit diff --git a/src/core_seaice/shared/mpas_seaice_berg_state.F b/src/core_seaice/shared/mpas_seaice_berg_state.F index 3c62d3ab17..98a88e9588 100644 --- a/src/core_seaice/shared/mpas_seaice_berg_state.F +++ b/src/core_seaice/shared/mpas_seaice_berg_state.F @@ -784,7 +784,7 @@ subroutine seaice_berg_init_displaced_area(domain) call MPAS_pool_get_subpool(block % structs, "berg_state", bergStatePool) call MPAS_pool_get_array(bergStatePool, "bergDisplacedArea", bergDisplacedArea) - bergDisplacedArea (:) = 0.0_RKIND + bergDisplacedArea(:) = 0.0_RKIND block => block % next enddo From 11b7ce36f575c1f844a41a048e53fbca7988fddd Mon Sep 17 00:00:00 2001 From: Darin Scott Comeau Date: Wed, 9 May 2018 08:40:28 -0600 Subject: [PATCH 4/5] Changing timeSeriesStatsMonthlyRestart input/output_interval to match standard_physics streams file. --- .../seaice/configurations/icebergs/streams.seaice | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/testing_and_setup/seaice/configurations/icebergs/streams.seaice b/testing_and_setup/seaice/configurations/icebergs/streams.seaice index 29676fa6d9..66af97e639 100644 --- a/testing_and_setup/seaice/configurations/icebergs/streams.seaice +++ b/testing_and_setup/seaice/configurations/icebergs/streams.seaice @@ -201,8 +201,8 @@ reference_time="0000-01-01_00:00:00" clobber_mode="truncate" packages="timeSeriesStatsMonthlyAMPKG" - input_interval="initial_only" - output_interval="stream:restart:output_interval" > + input_interval="none" + output_interval="none" > From c8c12f7bf57def5fc402fef0633e18945b45235b Mon Sep 17 00:00:00 2001 From: Darin Scott Comeau Date: Wed, 17 Oct 2018 14:55:50 -0600 Subject: [PATCH 5/5] Reverting changes made in #6. --- src/core_seaice/shared/mpas_seaice_berg_state.F | 2 +- .../seaice/configurations/icebergs/streams.seaice | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/core_seaice/shared/mpas_seaice_berg_state.F b/src/core_seaice/shared/mpas_seaice_berg_state.F index 98a88e9588..3c62d3ab17 100644 --- a/src/core_seaice/shared/mpas_seaice_berg_state.F +++ b/src/core_seaice/shared/mpas_seaice_berg_state.F @@ -784,7 +784,7 @@ subroutine seaice_berg_init_displaced_area(domain) call MPAS_pool_get_subpool(block % structs, "berg_state", bergStatePool) call MPAS_pool_get_array(bergStatePool, "bergDisplacedArea", bergDisplacedArea) - bergDisplacedArea(:) = 0.0_RKIND + bergDisplacedArea (:) = 0.0_RKIND block => block % next enddo diff --git a/testing_and_setup/seaice/configurations/icebergs/streams.seaice b/testing_and_setup/seaice/configurations/icebergs/streams.seaice index 66af97e639..29676fa6d9 100644 --- a/testing_and_setup/seaice/configurations/icebergs/streams.seaice +++ b/testing_and_setup/seaice/configurations/icebergs/streams.seaice @@ -201,8 +201,8 @@ reference_time="0000-01-01_00:00:00" clobber_mode="truncate" packages="timeSeriesStatsMonthlyAMPKG" - input_interval="none" - output_interval="none" > + input_interval="initial_only" + output_interval="stream:restart:output_interval" >