MODULE icerdgrft !!====================================================================== !! *** MODULE icerdgrft *** !! LIM-3 : Mechanical impact on ice thickness distribution !!====================================================================== !! History : LIM ! 2006-02 (M. Vancoppenolle) Original code !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_dyn !! 4.0 ! 2011-02 (G. Madec) dynamical allocation !!---------------------------------------------------------------------- #if defined key_lim3 !!---------------------------------------------------------------------- !! 'key_lim3' LIM-3 sea-ice model !!---------------------------------------------------------------------- USE par_oce ! ocean parameters USE dom_oce ! ocean domain USE phycst ! physical constants (ocean directory) USE sbc_oce , ONLY : sss_m, sst_m ! surface boundary condition: ocean fields USE ice1D ! sea-ice: thermodynamics USE ice ! sea-ice: variables USE icevar ! sea-ice: operations USE icectl ! sea-ice: control prints ! USE lbclnk ! lateral boundary condition - MPP exchanges USE lib_mpp ! MPP library USE in_out_manager ! I/O manager USE iom ! I/O manager USE lib_fortran ! glob_sum USE timing ! Timing IMPLICIT NONE PRIVATE PUBLIC ice_rdgrft ! called by ice_stp PUBLIC ice_rdgrft_strength ! called by icerhg_evp PUBLIC ice_rdgrft_init ! called by ice_stp PUBLIC ice_rdgrft_alloc ! called by ice_init ! Variables shared among ridging subroutines REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: asum ! sum of total ice and open water area REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: aksum ! ratio of area removed to area ridged REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: athorn ! participation function; fraction of ridging/closing associated w/ category n REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmin ! minimum ridge thickness REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmax ! maximum ridge thickness REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hraft ! thickness of rafted ice REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: krdg ! thickness of ridging ice / mean ridge thickness REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: aridge ! participating ice ridging REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: araft ! participating ice rafting ! REAL(wp), PARAMETER :: krdgmin = 1.1_wp ! min ridge thickness multiplier REAL(wp), PARAMETER :: kraft = 0.5_wp ! rafting multipliyer REAL(wp) :: zdrho ! ! ! ** namelist (namice_rdgrft) ** REAL(wp) :: rn_cs ! fraction of shearing energy contributing to ridging LOGICAL :: ln_partf_lin ! participation function linear (Thorndike et al. (1975)) REAL(wp) :: rn_gstar ! fractional area of young ice contributing to ridging LOGICAL :: ln_partf_exp ! participation function exponential (Lipscomb et al. (2007)) REAL(wp) :: rn_astar ! equivalent of G* for an exponential participation function LOGICAL :: ln_ridging ! ridging of ice or not REAL(wp) :: rn_hstar ! thickness that determines the maximal thickness of ridged ice REAL(wp) :: rn_porordg ! initial porosity of ridges (0.3 regular value) REAL(wp) :: rn_fsnwrdg ! fractional snow loss to the ocean during ridging REAL(wp) :: rn_fpndrdg ! fractional pond loss to the ocean during ridging LOGICAL :: ln_rafting ! rafting of ice or not REAL(wp) :: rn_hraft ! threshold thickness (m) for rafting / ridging REAL(wp) :: rn_craft ! coefficient for smoothness of the hyperbolic tangent in rafting REAL(wp) :: rn_fsnwrft ! fractional snow loss to the ocean during rafting REAL(wp) :: rn_fpndrft ! fractional pond loss to the ocean during rafting ! !!---------------------------------------------------------------------- !! NEMO/ICE 4.0 , NEMO Consortium (2017) !! $Id: icerdgrft.F90 8378 2017-07-26 13:55:59Z clem $ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS INTEGER FUNCTION ice_rdgrft_alloc() !!---------------------------------------------------------------------! !! *** ROUTINE ice_rdgrft_alloc *** !!---------------------------------------------------------------------! ALLOCATE( asum (jpi,jpj) , athorn(jpi,jpj,0:jpl) , aksum (jpi,jpj) , & & hrmin(jpi,jpj,jpl) , hraft (jpi,jpj,jpl) , aridge(jpi,jpj,jpl) , & & hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) , STAT=ice_rdgrft_alloc ) ! IF( ice_rdgrft_alloc /= 0 ) CALL ctl_warn( 'ice_rdgrft_alloc: failed to allocate arrays' ) ! END FUNCTION ice_rdgrft_alloc SUBROUTINE ice_rdgrft( kt ) !!---------------------------------------------------------------------! !! *** ROUTINE ice_rdgrft *** !! !! ** Purpose : computes the mechanical redistribution of ice thickness !! !! ** Method : Steps : !! 1) Thickness categories boundaries, ice / o.w. concentrations !! Ridge preparation !! 2) Dynamical inputs (closing rate, divu_adv, opning) !! 3) Ridging iteration !! 4) Ridging diagnostics !! 5) Heat, salt and freshwater fluxes !! 6) Compute increments of tate variables and come back to old values !! !! References : Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626. !! Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980. !! Rothrock, D. A., 1975: JGR, 80, 4514-4519. !! Thorndike et al., 1975, JGR, 80, 4501-4513. !! Bitz et al., JGR, 2001 !! Amundrud and Melling, JGR 2005 !! Babko et al., JGR 2002 !! !! This routine is based on CICE code and authors William H. Lipscomb, !! and Elizabeth C. Hunke, LANL are gratefully acknowledged !!--------------------------------------------------------------------! INTEGER, INTENT(in) :: kt ! number of iteration !! INTEGER :: ji, jj, jk, jl ! dummy loop index INTEGER :: niter ! local integer INTEGER :: iterate_ridging ! if =1, repeat the ridging REAL(wp) :: za ! local scalar REAL(wp), DIMENSION(jpi,jpj) :: closing_net ! net rate at which area is removed (1/s) ! ! (ridging ice area - area of new ridges) / dt REAL(wp), DIMENSION(jpi,jpj) :: divu_adv ! divu as implied by transport scheme (1/s) REAL(wp), DIMENSION(jpi,jpj) :: opning ! rate of opening due to divergence/shear REAL(wp), DIMENSION(jpi,jpj) :: closing_gross ! rate at which area removed, not counting area of new ridges ! INTEGER, PARAMETER :: nitermax = 20 ! REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b !!----------------------------------------------------------------------------- IF( nn_timing == 1 ) CALL timing_start('icerdgrft') IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*)'ice_rdgrft: ice ridging and rafting' IF(lwp) WRITE(numout,*)'~~~~~~~~~~' ! zdrho = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0 ! proport const for PE ! ENDIF ! ! conservation test IF( ln_limdiachk ) CALL ice_cons_hsm(0, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) !-----------------------------------------------------------------------------! ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons !-----------------------------------------------------------------------------! ! CALL ice_rdgrft_prep ! prepare ridging ! DO jj = 1, jpj ! Initialize arrays. DO ji = 1, jpi !-----------------------------------------------------------------------------! ! 2) Dynamical inputs (closing rate, divu_adv, opning) !-----------------------------------------------------------------------------! ! ! 2.1 closing_net !----------------- ! Compute the net rate of closing due to convergence ! and shear, based on Flato and Hibler (1995). ! ! The energy dissipation rate is equal to the net closing rate ! times the ice strength. ! ! NOTE: The NET closing rate is equal to the rate that open water ! area is removed, plus the rate at which ice area is removed by ! ridging, minus the rate at which area is added in new ridges. ! The GROSS closing rate is equal to the first two terms (open ! water closing and thin ice ridging) without the third term ! (thick, newly ridged ice). closing_net(ji,jj) = rn_cs * 0.5_wp * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) ! 2.2 divu_adv !-------------- ! Compute divu_adv, the divergence rate given by the transport/ ! advection scheme, which may not be equal to divu as computed ! from the velocity field. ! ! If divu_adv < 0, make sure the closing rate is large enough ! to give asum = 1.0 after ridging. divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice ! asum found in ridgeprep IF( divu_adv(ji,jj) < 0._wp ) closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) ) ! 2.3 opning !------------ ! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging. opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj) END DO END DO !-----------------------------------------------------------------------------! ! 3) Ridging iteration !-----------------------------------------------------------------------------! niter = 1 ! iteration counter iterate_ridging = 1 DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax ) ! 3.2 closing_gross !-----------------------------------------------------------------------------! ! Based on the ITD of ridging and ridged ice, convert the net ! closing rate to a gross closing rate. ! NOTE: 0 < aksum <= 1 closing_gross(:,:) = closing_net(:,:) / aksum(:,:) ! correction to closing rate and opening if closing rate is excessive !--------------------------------------------------------------------- ! Reduce the closing rate if more than 100% of the open water ! would be removed. Reduce the opening rate proportionately. DO jj = 1, jpj DO ji = 1, jpi za = ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice IF ( za < 0._wp .AND. za > - ato_i(ji,jj) ) THEN ! would lead to negative ato_i opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice ELSEIF( za > 0._wp .AND. za > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN ! would lead to ato_i > asum opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice ENDIF END DO END DO ! correction to closing rate / opening if excessive ice removal !--------------------------------------------------------------- ! Reduce the closing rate if more than 100% of any ice category ! would be removed. Reduce the opening rate proportionately. DO jl = 1, jpl DO jj = 1, jpj DO ji = 1, jpi za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice IF( za > a_i(ji,jj,jl) ) THEN closing_gross(ji,jj) = closing_gross(ji,jj) * a_i(ji,jj,jl) / za ENDIF END DO END DO END DO ! 3.3 Redistribute area, volume, and energy. !-----------------------------------------------------------------------------! CALL ice_rdgrft_ridgeshift( opning, closing_gross ) ! 3.4 Compute total area of ice plus open water after ridging. !-----------------------------------------------------------------------------! ! This is in general not equal to one because of divergence during transport asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 ) ! 3.5 Do we keep on iterating? !-----------------------------------------------------------------------------! ! Check whether asum = 1. If not (because the closing and opening ! rates were reduced above), ridge again with new rates. iterate_ridging = 0 DO jj = 1, jpj DO ji = 1, jpi IF( ABS( asum(ji,jj) - 1._wp ) < epsi10 ) THEN closing_net(ji,jj) = 0._wp opning (ji,jj) = 0._wp ato_i (ji,jj) = MAX( 0._wp, 1._wp - SUM( a_i(ji,jj,:) ) ) ELSE iterate_ridging = 1 divu_adv (ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) ) opning (ji,jj) = MAX( 0._wp, divu_adv(ji,jj) ) ENDIF END DO END DO IF( lk_mpp ) CALL mpp_max( iterate_ridging ) ! Repeat if necessary. ! NOTE: If strength smoothing is turned on, the ridging must be ! iterated globally because of the boundary update in the smoothing. niter = niter + 1 ! IF( iterate_ridging == 1 ) THEN CALL ice_rdgrft_prep IF( niter > nitermax ) THEN WRITE(numout,*) ' ALERTE : non-converging ridging scheme ' WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging ENDIF ENDIF ! END DO !! on the do while over iter CALL ice_var_agg( 1 ) !-----------------------------------------------------------------------------! ! control prints !-----------------------------------------------------------------------------! ! ! conservation test IF( ln_limdiachk ) CALL ice_cons_hsm(1, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) ! ! control prints IF( ln_ctl ) CALL ice_prt3D( 'icerdgrft' ) ! IF( nn_timing == 1 ) CALL timing_stop('icerdgrft') ! END SUBROUTINE ice_rdgrft SUBROUTINE ice_rdgrft_prep !!---------------------------------------------------------------------! !! *** ROUTINE ice_rdgrft_prep *** !! !! ** Purpose : preparation for ridging and strength calculations !! !! ** Method : Compute the thickness distribution of the ice and open water !! participating in ridging and of the resulting ridges. !!---------------------------------------------------------------------! INTEGER :: ji, jj, jl ! dummy loop indices REAL(wp) :: z1_gstar, z1_astar, zhmean, zdummy ! local scalar REAL(wp), DIMENSION(jpi,jpj,-1:jpl) :: zGsum ! zGsum(n) = sum of areas in categories 0 to n !------------------------------------------------------------------------------! z1_gstar = 1._wp / rn_gstar z1_astar = 1._wp / rn_astar CALL ice_var_zapsmall ! Zero out categories with very small areas ! ! Ice thickness needed for rafting WHERE( a_i(:,:,:) >= epsi20 ) ; ht_i(:,:,:) = v_i (:,:,:) / a_i(:,:,:) ELSEWHERE ; ht_i(:,:,:) = 0._wp END WHERE !------------------------------------------------------------------------------! ! 1) Participation function !------------------------------------------------------------------------------! ! ! Compute total area of ice plus open water. ! This is in general not equal to one because of divergence during transport asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 ) ! ! Compute cumulative thickness distribution function ! Compute the cumulative thickness distribution function zGsum, ! where zGsum(n) is the fractional area in categories 0 to n. ! initial value (in h = 0) equals open water area zGsum(:,:,-1) = 0._wp zGsum(:,:,0 ) = ato_i(:,:) / asum(:,:) DO jl = 1, jpl zGsum(:,:,jl) = ( ato_i(:,:) + SUM( a_i(:,:,1:jl), dim=3 ) ) / asum(:,:) ! sum(1:jl) is ok (and not jpl) END DO ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn) !-------------------------------------------------------------------------------------------------- ! Compute the participation function athorn; this is analogous to ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975). ! area lost from category n due to ridging/closing ! athorn(n) = total area lost due to ridging/closing ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar). ! ! The expressions for athorn are found by integrating b(h)g(h) between ! the category boundaries. ! athorn is always >= 0 and SUM(athorn(0:jpl))=1 !----------------------------------------------------------------- ! IF( ln_partf_lin ) THEN !--- Linear formulation (Thorndike et al., 1975) DO jl = 0, jpl DO jj = 1, jpj DO ji = 1, jpi IF ( zGsum(ji,jj,jl) < rn_gstar ) THEN athorn(ji,jj,jl) = z1_gstar * ( zGsum(ji,jj,jl) - zGsum(ji,jj,jl-1) ) * & & ( 2._wp - ( zGsum(ji,jj,jl-1) + zGsum(ji,jj,jl) ) * z1_gstar ) ELSEIF( zGsum(ji,jj,jl-1) < rn_gstar ) THEN athorn(ji,jj,jl) = z1_gstar * ( rn_gstar - zGsum(ji,jj,jl-1) ) * & & ( 2._wp - ( zGsum(ji,jj,jl-1) + rn_gstar ) * z1_gstar ) ELSE athorn(ji,jj,jl) = 0._wp ENDIF END DO END DO END DO ! ELSEIF( ln_partf_exp ) THEN !--- Exponential, more stable formulation (Lipscomb et al, 2007) ! zdummy = 1._wp / ( 1._wp - EXP(-z1_astar) ) ! precompute exponential terms using zGsum as a work array DO jl = -1, jpl zGsum(:,:,jl) = EXP( -zGsum(:,:,jl) * z1_astar ) * zdummy END DO DO jl = 0, jpl athorn(:,:,jl) = zGsum(:,:,jl-1) - zGsum(:,:,jl) END DO ! ENDIF ! !--- Ridging and rafting participation concentrations IF( ln_rafting .AND. ln_ridging ) THEN !- ridging & rafting DO jl = 1, jpl DO jj = 1, jpj DO ji = 1, jpi aridge(ji,jj,jl) = ( 1._wp + TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) ) * 0.5_wp * athorn(ji,jj,jl) araft (ji,jj,jl) = athorn(ji,jj,jl) - aridge(ji,jj,jl) END DO END DO END DO ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN !- ridging alone aridge(:,:,:) = athorn(:,:,1:jpl) araft (:,:,:) = 0._wp ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN !- rafting alone aridge(:,:,:) = 0._wp araft (:,:,:) = athorn(:,:,1:jpl) ELSE !- no ridging & no rafting aridge(:,:,:) = 0._wp araft (:,:,:) = 0._wp ENDIF !----------------------------------------------------------------- ! 2) Transfer function !----------------------------------------------------------------- ! Compute max and min ridged ice thickness for each ridging category. ! Assume ridged ice is uniformly distributed between hrmin and hrmax. ! ! This parameterization is a modified version of Hibler (1980). ! The mean ridging thickness, zhmean, is proportional to hi^(0.5) ! and for very thick ridging ice must be >= krdgmin*hi ! ! The minimum ridging thickness, hrmin, is equal to 2*hi ! (i.e., rafting) and for very thick ridging ice is ! constrained by hrmin <= (zhmean + hi)/2. ! ! The maximum ridging thickness, hrmax, is determined by ! zhmean and hrmin. ! ! These modifications have the effect of reducing the ice strength ! (relative to the Hibler formulation) when very thick ice is ! ridging. ! ! aksum = net area removed/ total area removed ! where total area removed = area of ice that ridges ! net area removed = total area removed - area of new ridges !----------------------------------------------------------------- aksum(:,:) = athorn(:,:,0) ! Transfer function DO jl = 1, jpl !all categories have a specific transfer function DO jj = 1, jpj DO ji = 1, jpi IF ( athorn(ji,jj,jl) > 0._wp ) THEN zhmean = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin ) hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( zhmean + ht_i(ji,jj,jl) ) ) hrmax(ji,jj,jl) = 2._wp * zhmean - hrmin(ji,jj,jl) hraft(ji,jj,jl) = ht_i(ji,jj,jl) / kraft krdg (ji,jj,jl) = ht_i(ji,jj,jl) / MAX( zhmean, epsi20 ) ! ! Normalization factor : aksum, ensures mass conservation aksum(ji,jj) = aksum(ji,jj) + aridge(ji,jj,jl) * ( 1._wp - krdg(ji,jj,jl) ) & & + araft (ji,jj,jl) * ( 1._wp - kraft ) ELSE hrmin(ji,jj,jl) = 0._wp hrmax(ji,jj,jl) = 0._wp hraft(ji,jj,jl) = 0._wp krdg (ji,jj,jl) = 1._wp ENDIF END DO END DO END DO ! END SUBROUTINE ice_rdgrft_prep SUBROUTINE ice_rdgrft_ridgeshift( opning, closing_gross ) !!---------------------------------------------------------------------- !! *** ROUTINE ice_rdgrft_strength *** !! !! ** Purpose : shift ridging ice among thickness categories of ice thickness !! !! ** Method : Remove area, volume, and energy from each ridging category !! and add to thicker ice categories. !!---------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: opning ! rate of opening due to divergence/shear REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: closing_gross ! rate at which area retreats, excluding area of new ridges ! INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices INTEGER :: ij ! horizontal index, combines i and j loops INTEGER :: icells ! number of cells with a_i > puny REAL(wp) :: hL, hR, farea ! left and right limits of integration and new area going to jl2 INTEGER , DIMENSION(jpij) :: indxi, indxj ! compressed indices REAL(wp), DIMENSION(jpij) :: zswitch, fvol ! new ridge volume going to jl2 REAL(wp), DIMENSION(jpij) :: afrac ! fraction of category area ridged REAL(wp), DIMENSION(jpij) :: ardg1 , ardg2 ! area of ice ridged & new ridges REAL(wp), DIMENSION(jpij) :: vsrdg , esrdg ! snow volume & energy of ridging ice ! MV MP 2016 REAL(wp), DIMENSION(jpij) :: vprdg ! pond volume of ridging ice REAL(wp), DIMENSION(jpij) :: aprdg1 ! pond area of ridging ice REAL(wp), DIMENSION(jpij) :: aprdg2 ! pond area of ridging ice ! END MV MP 2016 REAL(wp), DIMENSION(jpij) :: dhr , dhr2 ! hrmax - hrmin & hrmax^2 - hrmin^2 REAL(wp), DIMENSION(jpij) :: vrdg1 ! volume of ice ridged REAL(wp), DIMENSION(jpij) :: vrdg2 ! volume of new ridges REAL(wp), DIMENSION(jpij) :: vsw ! volume of seawater trapped into ridges REAL(wp), DIMENSION(jpij) :: srdg1 ! sal*volume of ice ridged REAL(wp), DIMENSION(jpij) :: srdg2 ! sal*volume of new ridges REAL(wp), DIMENSION(jpij) :: smsw ! sal*volume of water trapped into ridges REAL(wp), DIMENSION(jpij) :: oirdg1, oirdg2 ! ice age of ice ridged REAL(wp), DIMENSION(jpij) :: afrft ! fraction of category area rafted REAL(wp), DIMENSION(jpij) :: arft1 , arft2 ! area of ice rafted and new rafted zone REAL(wp), DIMENSION(jpij) :: virft , vsrft ! ice & snow volume of rafting ice ! MV MP 2016 REAL(wp), DIMENSION(jpij) :: vprft ! pond volume of rafting ice REAL(wp), DIMENSION(jpij) :: aprft1 ! pond area of rafted ice REAL(wp), DIMENSION(jpij) :: aprft2 ! pond area of new rafted ice ! END MV MP 2016 REAL(wp), DIMENSION(jpij) :: esrft , smrft ! snow energy & salinity of rafting ice REAL(wp), DIMENSION(jpij) :: oirft1, oirft2 ! ice age of ice rafted REAL(wp), DIMENSION(jpij,nlay_i) :: eirft ! ice energy of rafting ice REAL(wp), DIMENSION(jpij,nlay_i) :: erdg1 ! enth*volume of ice ridged REAL(wp), DIMENSION(jpij,nlay_i) :: erdg2 ! enth*volume of new ridges REAL(wp), DIMENSION(jpij,nlay_i) :: ersw ! enth of water trapped into ridges !!---------------------------------------------------------------------- !------------------------------------------------------------------------------- ! 1) Compute change in open water area due to closing and opening. !------------------------------------------------------------------------------- DO jj = 1, jpj DO ji = 1, jpi ato_i(ji,jj) = MAX( 0._wp, ato_i(ji,jj) + & & ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice ) END DO END DO !----------------------------------------------------------------- ! 3) Pump everything from ice which is being ridged / rafted !----------------------------------------------------------------- ! Compute the area, volume, and energy of ice ridging in each ! category, along with the area of the resulting ridge. DO jl1 = 1, jpl !jl1 describes the ridging category !------------------------------------------------ ! 3.1) Identify grid cells with nonzero ridging !------------------------------------------------ icells = 0 DO jj = 1, jpj DO ji = 1, jpi IF( athorn(ji,jj,jl1) > 0._wp .AND. closing_gross(ji,jj) > 0._wp ) THEN icells = icells + 1 indxi(icells) = ji indxj(icells) = jj ENDIF END DO END DO DO ij = 1, icells ji = indxi(ij) ; jj = indxj(ij) !-------------------------------------------------------------------- ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2) !-------------------------------------------------------------------- ardg1(ij) = aridge(ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice arft1(ij) = araft (ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice !--------------------------------------------------------------- ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1 !--------------------------------------------------------------- afrac(ij) = ardg1(ij) / a_i(ji,jj,jl1) !ridging afrft(ij) = arft1(ij) / a_i(ji,jj,jl1) !rafting ardg2(ij) = ardg1(ij) * krdg(ji,jj,jl1) arft2(ij) = arft1(ij) * kraft !-------------------------------------------------------------------------- ! 3.4) Substract area, volume, and energy from ridging ! / rafting category n1. !-------------------------------------------------------------------------- vrdg1(ij) = v_i(ji,jj,jl1) * afrac(ij) vrdg2(ij) = vrdg1(ij) * ( 1. + rn_porordg ) vsw (ij) = vrdg1(ij) * rn_porordg vsrdg (ij) = v_s (ji,jj, jl1) * afrac(ij) esrdg (ij) = e_s (ji,jj,1,jl1) * afrac(ij) !MV MP 2016 IF ( nn_pnd_scheme > 0 ) THEN aprdg1(ij) = a_ip(ji,jj, jl1) * afrac(ij) aprdg2(ij) = a_ip(ji,jj, jl1) * afrac(ij) * krdg(ji,jj,jl1) vprdg(ij) = v_ip(ji,jj, jl1) * afrac(ij) ENDIF ! END MV MP 2016 srdg1 (ij) = smv_i(ji,jj, jl1) * afrac(ij) oirdg1(ij) = oa_i (ji,jj, jl1) * afrac(ij) oirdg2(ij) = oa_i (ji,jj, jl1) * afrac(ij) * krdg(ji,jj,jl1) ! rafting volumes, heat contents ... virft (ij) = v_i (ji,jj, jl1) * afrft(ij) vsrft (ij) = v_s (ji,jj, jl1) * afrft(ij) !MV MP 2016 IF ( nn_pnd_scheme > 0 ) THEN aprft1(ij) = a_ip (ji,jj, jl1) * afrft(ij) aprft2(ij) = a_ip (ji,jj, jl1) * afrft(ij) * kraft vprft(ij) = v_ip(ji,jj,jl1) * afrft(ij) ENDIF ! END MV MP 2016 srdg1 (ij) = smv_i(ji,jj, jl1) * afrac(ij) esrft (ij) = e_s (ji,jj,1,jl1) * afrft(ij) smrft (ij) = smv_i(ji,jj, jl1) * afrft(ij) oirft1(ij) = oa_i (ji,jj, jl1) * afrft(ij) oirft2(ij) = oa_i (ji,jj, jl1) * afrft(ij) * kraft !----------------------------------------------------------------- ! 3.5) Compute properties of new ridges !----------------------------------------------------------------- smsw(ij) = vsw(ij) * sss_m(ji,jj) ! salt content of seawater frozen in voids srdg2(ij) = srdg1(ij) + smsw(ij) ! salt content of new ridge sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ij) * rhoic * r1_rdtice wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ij) * rhoic * r1_rdtice ! increase in ice volume due to seawater frozen in voids ! virtual salt flux to keep salinity constant IF( nn_icesal == 1 .OR. nn_icesal == 3 ) THEN srdg2(ij) = srdg2(ij) - vsw(ij) * ( sss_m(ji,jj) - sm_i(ji,jj,jl1) ) ! ridge salinity = sm_i sfx_bri(ji,jj) = sfx_bri(ji,jj) + sss_m(ji,jj) * vsw(ij) * rhoic * r1_rdtice & ! put back sss_m into the ocean & - sm_i(ji,jj,jl1) * vsw(ij) * rhoic * r1_rdtice ! and get sm_i from the ocean ENDIF !------------------------------------------ ! 3.7 Put the snow somewhere in the ocean !------------------------------------------ ! Place part of the snow lost by ridging into the ocean. ! Note that esrdg > 0; the ocean must cool to melt snow. ! If the ocean temp = Tf already, new ice must grow. ! During the next time step, thermo_rates will determine whether ! the ocean cools or new ice grows. wfx_snw_dyn(ji,jj) = wfx_snw_dyn(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnwrdg ) & & + rhosn * vsrft(ij) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice ! fresh water source for ocean hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnwrdg ) & & - esrft(ij) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice ! heat sink for ocean (<0, W.m-2) ! MV MP 2016 !------------------------------------------ ! 3.X Put the melt pond water in the ocean !------------------------------------------ ! Place part of the melt pond volume into the ocean. IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) THEN wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + ( rhofw * vprdg(ij) * ( 1._wp - rn_fpndrdg ) & & + rhofw * vprft(ij) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice ! fresh water source for ocean ENDIF ! END MV MP 2016 !----------------------------------------------------------------- ! 3.8 Compute quantities used to apportion ice among categories ! in the n2 loop below !----------------------------------------------------------------- dhr (ij) = 1._wp / ( hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) ) dhr2(ij) = 1._wp / ( hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) ) ! update jl1 (removing ridged/rafted area) a_i (ji,jj, jl1) = a_i (ji,jj, jl1) - ardg1 (ij) - arft1 (ij) v_i (ji,jj, jl1) = v_i (ji,jj, jl1) - vrdg1 (ij) - virft (ij) v_s (ji,jj, jl1) = v_s (ji,jj, jl1) - vsrdg (ij) - vsrft (ij) e_s (ji,jj,1,jl1) = e_s (ji,jj,1,jl1) - esrdg (ij) - esrft (ij) smv_i(ji,jj, jl1) = smv_i(ji,jj, jl1) - srdg1 (ij) - smrft (ij) oa_i (ji,jj, jl1) = oa_i (ji,jj, jl1) - oirdg1(ij) - oirft1(ij) ! MV MP 2016 IF ( nn_pnd_scheme > 0 ) THEN v_ip (ji,jj,jl1) = v_ip (ji,jj,jl1) - vprdg (ij) - vprft (ij) a_ip (ji,jj,jl1) = a_ip (ji,jj,jl1) - aprdg1(ij) - aprft1(ij) ENDIF ! END MV MP 2016 END DO !-------------------------------------------------------------------- ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and ! compute ridged ice enthalpy !-------------------------------------------------------------------- DO jk = 1, nlay_i DO ij = 1, icells ji = indxi(ij) ; jj = indxj(ij) ! heat content of ridged ice erdg1(ij,jk) = e_i(ji,jj,jk,jl1) * afrac(ij) eirft(ij,jk) = e_i(ji,jj,jk,jl1) * afrft(ij) ! enthalpy of the trapped seawater (J/m2, >0) ! clem: if sst>0, then ersw <0 (is that possible?) ersw(ij,jk) = - rhoic * vsw(ij) * rcp * sst_m(ji,jj) * r1_nlay_i ! heat flux to the ocean hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ij,jk) * r1_rdtice ! > 0 [W.m-2] ocean->ice flux ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean erdg2(ij,jk) = erdg1(ij,jk) + ersw(ij,jk) ! update jl1 e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk) END DO END DO !------------------------------------------------------------------------------- ! 4) Add area, volume, and energy of new ridge to each category jl2 !------------------------------------------------------------------------------- DO jl2 = 1, jpl ! over categories to which ridged/rafted ice is transferred DO ij = 1, icells ji = indxi(ij) ; jj = indxj(ij) ! Compute the fraction of ridged ice area and volume going to thickness category jl2. IF( hrmin(ji,jj,jl1) <= hi_max(jl2) .AND. hrmax(ji,jj,jl1) > hi_max(jl2-1) ) THEN hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) ) hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2) ) farea = ( hR - hL ) * dhr(ij) fvol(ij) = ( hR * hR - hL * hL ) * dhr2(ij) ELSE farea = 0._wp fvol(ij) = 0._wp ENDIF ! Compute the fraction of rafted ice area and volume going to thickness category jl2 !!gm see above IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN IF( hi_max(jl2-1) < hraft(ji,jj,jl1) .AND. hraft(ji,jj,jl1) <= hi_max(jl2) ) THEN ; zswitch(ij) = 1._wp ELSE ; zswitch(ij) = 0._wp ENDIF ! a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + ( ardg2 (ij) * farea + arft2 (ij) * zswitch(ij) ) oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + ( oirdg2(ij) * farea + oirft2(ij) * zswitch(ij) ) v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + ( vrdg2 (ij) * fvol(ij) + virft (ij) * zswitch(ij) ) smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + ( srdg2 (ij) * fvol(ij) + smrft (ij) * zswitch(ij) ) v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + ( vsrdg (ij) * rn_fsnwrdg * fvol(ij) + & & vsrft (ij) * rn_fsnwrft * zswitch(ij) ) e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnwrdg * fvol(ij) + & & esrft (ij) * rn_fsnwrft * zswitch(ij) ) ! MV MP 2016 IF ( nn_pnd_scheme > 0 ) THEN v_ip (ji,jj,jl2) = v_ip(ji,jj,jl2) + ( vprdg (ij) * rn_fpndrdg * fvol (ij) & & + vprft (ij) * rn_fpndrft * zswitch(ij) ) a_ip (ji,jj,jl2) = a_ip(ji,jj,jl2) + ( aprdg2(ij) * rn_fpndrdg * farea & & + aprft2(ij) * rn_fpndrft * zswitch(ji) ) ENDIF ! END MV MP 2016 END DO ! Transfer ice energy to category jl2 by ridging DO jk = 1, nlay_i DO ij = 1, icells ji = indxi(ij) ; jj = indxj(ij) e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + erdg2(ij,jk) * fvol(ij) + eirft(ij,jk) * zswitch(ij) END DO END DO ! END DO ! jl2 ! END DO ! jl1 (deforming categories) ! END SUBROUTINE ice_rdgrft_ridgeshift SUBROUTINE ice_rdgrft_strength !!---------------------------------------------------------------------- !! *** ROUTINE ice_rdgrft_strength *** !! !! ** Purpose : computes ice strength used in dynamics routines of ice thickness !! !! ** Method : Compute the strength of the ice pack, defined as the energy (J m-2) !! dissipated per unit area removed from the ice pack under compression, !! and assumed proportional to the change in potential energy caused !! by ridging. Note that only Hibler's formulation is stable and that !! ice strength has to be smoothed !!---------------------------------------------------------------------- INTEGER :: ji,jj, jl ! dummy loop indices INTEGER :: ismooth ! smoothing the resistance to deformation INTEGER :: itframe ! number of time steps for the P smoothing REAL(wp) :: zp, z1_3 ! local scalars REAL(wp), DIMENSION(jpi,jpj) :: zworka ! temporary array used here REAL(wp), DIMENSION(jpi,jpj) :: zstrp1, zstrp2 ! strength at previous time steps !!---------------------------------------------------------------------- ! !--------------------------------------------------! CALL ice_rdgrft_prep ! Thickness distribution of ridging and ridged ice ! ! !--------------------------------------------------! ! !--------------------------------------------------! IF( ln_str_Rot ) THEN ! Ice strength => Rothrock (1975) method ! ! !--------------------------------------------------! z1_3 = 1._wp / 3._wp DO jl = 1, jpl WHERE( athorn(:,:,jl) > 0._wp ) strength(:,:) = - athorn(:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl) & ! PE loss from deforming ice & + 2._wp * araft (:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl) & ! PE gain from rafting ice & + aridge(:,:,jl) * krdg(:,:,jl) * z1_3 * & ! PE gain from ridging ice & ( hrmax(:,:,jl) * hrmax(:,:,jl) + & & hrmin(:,:,jl) * hrmin(:,:,jl) + & & hrmax(:,:,jl) * hrmin(:,:,jl) ) ELSEWHERE strength(:,:) = 0._wp END WHERE END DO strength(:,:) = rn_perdg * zdrho * strength(:,:) / aksum(:,:) * tmask(:,:,1) ! where zdrho = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_perdg accounts for frictional dissipation ismooth = 1 ! !--------------------------------------------------! ELSEIF( ln_str_Hib ) THEN ! Ice strength => Hibler (1979) method ! ! !--------------------------------------------------! strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) ) * tmask(:,:,1) ! ismooth = 1 ! ENDIF ! !--------------------------------------------------! SELECT CASE( ismooth ) ! Smoothing ice strength ! ! !--------------------------------------------------! CASE( 1 ) !--- Spatial smoothing DO jj = 2, jpjm1 DO ji = 2, jpim1 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN zworka(ji,jj) = ( 4.0 * strength(ji,jj) & & + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & & + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) & & ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) ) ELSE zworka(ji,jj) = 0._wp ENDIF END DO END DO DO jj = 2, jpjm1 DO ji = 2, jpim1 strength(ji,jj) = zworka(ji,jj) END DO END DO CALL lbc_lnk( strength, 'T', 1. ) ! CASE( 2 ) !--- Temporal smoothing IF ( kt_ice == nit000 ) THEN zstrp1(:,:) = 0._wp zstrp2(:,:) = 0._wp ENDIF ! DO jj = 2, jpjm1 DO ji = 2, jpim1 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN itframe = 1 ! number of time steps for the running mean IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1 IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1 zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe zstrp2 (ji,jj) = zstrp1 (ji,jj) zstrp1 (ji,jj) = strength(ji,jj) strength(ji,jj) = zp ENDIF END DO END DO CALL lbc_lnk( strength, 'T', 1. ) ! END SELECT ! END SUBROUTINE ice_rdgrft_strength SUBROUTINE ice_rdgrft_init !!------------------------------------------------------------------- !! *** ROUTINE ice_rdgrft_init *** !! !! ** Purpose : Physical constants and parameters linked !! to the mechanical ice redistribution !! !! ** Method : Read the namice_rdgrft namelist !! and check the parameters values !! called at the first timestep (nit000) !! !! ** input : Namelist namice_rdgrft !!------------------------------------------------------------------- INTEGER :: ios ! Local integer output status for namelist read !! NAMELIST/namice_rdgrft/ ln_str_Hib, rn_pstar, rn_crhg, & & ln_str_Rot, rn_perdg, & & rn_cs , & & ln_partf_lin, rn_gstar, & & ln_partf_exp, rn_astar, & & ln_ridging, rn_hstar, rn_porordg, rn_fsnwrdg, rn_fpndrdg, & & ln_rafting, rn_hraft, rn_craft , rn_fsnwrft, rn_fpndrft !!------------------------------------------------------------------- ! REWIND( numnam_ice_ref ) ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution READ ( numnam_ice_ref, namice_rdgrft, IOSTAT = ios, ERR = 901) 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_rdgrft in reference namelist', lwp ) ! REWIND( numnam_ice_cfg ) ! Namelist namice_rdgrft in configuration namelist : Ice mechanical ice redistribution READ ( numnam_ice_cfg, namice_rdgrft, IOSTAT = ios, ERR = 902 ) 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_rdgrft in configuration namelist', lwp ) IF(lwm) WRITE ( numoni, namice_rdgrft ) ! IF (lwp) THEN ! control print WRITE(numout,*) WRITE(numout,*) 'ice_rdgrft_init : ice parameters for ridging/rafting ' WRITE(numout,*) '~~~~~~~~~~~~~~~' WRITE(numout,*) ' Namelist namice_rdgrft' WRITE(numout,*) ' ice strength parameterization Hibler (1979) ln_str_Hib = ', ln_str_Hib WRITE(numout,*) ' 1st bulk-rheology parameter rn_pstar = ', rn_pstar WRITE(numout,*) ' 2nd bulk-rhelogy parameter rn_crhg = ', rn_crhg WRITE(numout,*) ' ice strength parameterization Rothrock (1975) ln_str_Rot = ', ln_str_Rot WRITE(numout,*) ' Ratio of ridging work to PotEner change in ridging rn_perdg = ', rn_perdg WRITE(numout,*) ' Fraction of shear energy contributing to ridging rn_cs = ', rn_cs WRITE(numout,*) ' linear ridging participation function ln_partf_lin = ', ln_partf_lin WRITE(numout,*) ' Fraction of ice coverage contributing to ridging rn_gstar = ', rn_gstar WRITE(numout,*) ' Exponential ridging participation function ln_partf_exp = ', ln_partf_exp WRITE(numout,*) ' Equivalent to G* for an exponential function rn_astar = ', rn_astar WRITE(numout,*) ' Ridging of ice sheets or not ln_ridging = ', ln_ridging WRITE(numout,*) ' max ridged ice thickness rn_hstar = ', rn_hstar WRITE(numout,*) ' Initial porosity of ridges rn_porordg = ', rn_porordg WRITE(numout,*) ' Fraction of snow volume conserved during ridging rn_fsnwrdg = ', rn_fsnwrdg WRITE(numout,*) ' Fraction of pond volume conserved during ridging rn_fpndrdg = ', rn_fpndrdg WRITE(numout,*) ' Rafting of ice sheets or not ln_rafting = ', ln_rafting WRITE(numout,*) ' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft WRITE(numout,*) ' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft WRITE(numout,*) ' Fraction of snow volume conserved during rafting rn_fsnwrft = ', rn_fsnwrft WRITE(numout,*) ' Fraction of pond volume conserved during rafting rn_fpndrft = ', rn_fpndrft ENDIF ! IF ( ( ln_str_Hib .AND. ln_str_Rot ) .OR. ( .NOT.ln_str_Hib .AND. .NOT.ln_str_Rot ) ) THEN CALL ctl_stop( 'ice_rdgrft_init: choose one and only one formulation for ice strength (ln_str_Hib or ln_str_Rot)' ) ENDIF ! IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN CALL ctl_stop( 'ice_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' ) ENDIF ! END SUBROUTINE ice_rdgrft_init #else !!---------------------------------------------------------------------- !! Default option Empty module NO LIM-3 sea-ice model !!---------------------------------------------------------------------- #endif !!====================================================================== END MODULE icerdgrft