MODULE icedyn_rdgrft !!====================================================================== !! *** MODULE icedyn_rdgrft *** !! sea-ice : 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' ESIM sea-ice model !!---------------------------------------------------------------------- !! ice_dyn_rdgrft : ridging/rafting of sea ice !! ice_dyn_rdgrft_init : initialization of ridging/rafting of sea ice !! ice_strength : ice strength calculation !!---------------------------------------------------------------------- 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 icetab ! sea-ice: 1D <==> 2D transformation USE icevar ! sea-ice: operations USE icectl ! sea-ice: control prints ! USE in_out_manager ! I/O manager USE iom ! I/O manager library USE lib_mpp ! MPP library USE lib_fortran ! fortran utilities (glob_sum + no signed zero) USE lbclnk ! lateral boundary conditions (or mpp links) USE timing ! Timing IMPLICIT NONE PRIVATE PUBLIC ice_dyn_rdgrft ! called by icestp PUBLIC ice_dyn_rdgrft_init ! called by icedyn PUBLIC ice_strength ! called by icedyn_rhg_evp ! Variables shared among ridging subroutines REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: closing_net ! net rate at which area is removed (1/s) ! ! (ridging ice area - area of new ridges) / dt REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: opning ! rate of opening due to divergence/shear REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: closing_gross ! rate at which area removed, not counting area of new ridges REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: apartf ! participation function; fraction of ridging/closing associated w/ category n REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrmin ! minimum ridge thickness REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrmax ! maximum ridge thickness REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hraft ! thickness of rafted ice REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hi_hrdg ! thickness of ridging ice / mean ridge thickness REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: aridge ! participating ice ridging REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: araft ! participating ice rafting REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ze_i_2d REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ze_s_2d ! REAL(wp), PARAMETER :: hrdg_hi_min = 1.1_wp ! min ridge thickness multiplier: min(hrdg/hi) REAL(wp), PARAMETER :: hi_hrft = 0.5_wp ! rafting multipliyer: (hi/hraft) ! ! ** namelist (namdyn_rdgrft) ** LOGICAL :: ln_str_H79 ! ice strength parameterization (Hibler79) REAL(wp) :: rn_pstar ! determines ice strength, Hibler JPO79 REAL(wp) :: rn_crhg ! determines changes in ice strength LOGICAL :: ln_str_R75 ! ice strength parameterization (Rothrock75) REAL(wp) :: rn_perdg ! ridging work divided by pot. energy change in ridging REAL(wp) :: rn_csrdg ! 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: icedyn_rdgrft.F90 8378 2017-07-26 13:55:59Z clem $ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS INTEGER FUNCTION ice_dyn_rdgrft_alloc() !!------------------------------------------------------------------- !! *** ROUTINE ice_dyn_rdgrft_alloc *** !!------------------------------------------------------------------- ALLOCATE( closing_net(jpij), opning(jpij) , closing_gross(jpij), & & apartf(jpij,0:jpl), hrmin(jpij,jpl), hraft(jpij,jpl) , aridge(jpij,jpl), & & hrmax(jpij,jpl), hi_hrdg(jpij,jpl) , araft (jpij,jpl), & & ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc ) IF( lk_mpp ) CALL mpp_sum ( ice_dyn_rdgrft_alloc ) IF( ice_dyn_rdgrft_alloc /= 0 ) CALL ctl_warn( 'ice_dyn_rdgrft_alloc: failed to allocate arrays' ) ! END FUNCTION ice_dyn_rdgrft_alloc SUBROUTINE ice_dyn_rdgrft( kt ) !!------------------------------------------------------------------- !! *** ROUTINE ice_dyn_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, iterate_ridging ! local integer INTEGER :: nidx2 ! local integer REAL(wp) :: zfac ! local scalar INTEGER , DIMENSION(jpij) :: idxice2 ! compute ridge/raft or not REAL(wp), DIMENSION(jpij) :: zdivu_adv ! divu as implied by transport scheme (1/s) REAL(wp), DIMENSION(jpij) :: zdivu, zdelt ! 1D divu_i & delta_i ! INTEGER, PARAMETER :: nitermax = 20 !!------------------------------------------------------------------- ! controls IF( nn_timing == 1 ) CALL timing_start('icedyn_rdgrft') ! timing IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*)'ice_dyn_rdgrft: ice ridging and rafting' IF(lwp) WRITE(numout,*)'~~~~~~~~~~~~~~' ENDIF CALL ice_var_zapsmall ! Zero out categories with very small areas !-------------------------------- ! 0) Identify grid cells with ice !-------------------------------- nidx = 0 ; idxice(:) = 0 DO jj = 1, jpj DO ji = 1, jpi IF ( SUM(a_i(ji,jj,:)) > 0._wp ) THEN nidx = nidx + 1 idxice( nidx ) = (jj - 1) * jpi + ji ENDIF END DO END DO IF( nidx > 0 ) THEN ! just needed here CALL tab_2d_1d( nidx, idxice(1:nidx), zdivu(1:nidx), divu_i(:,:) ) CALL tab_2d_1d( nidx, idxice(1:nidx), zdelt(1:nidx), delta_i(:,:) ) ! needed here and in the iteration loop CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i(:,:,:) ) CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i(:,:,:) ) CALL tab_2d_1d( nidx, idxice(1:nidx), ato_i_1d(1:nidx) , ato_i(:,:) ) DO ji = 1, nidx !-----------------------------------------------------------------------------! ! 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) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 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. zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice IF( zdivu_adv(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) ) ! 2.3 opning !------------ ! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging. opning(ji) = closing_net(ji) + zdivu_adv(ji) END DO ! !------------------------------------------------ ! 2.1) Identify grid cells with nonzero ridging !------------------------------------------------ CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d ) nidx2 = 0 ; idxice2(:) = 0 DO ji = 1, nidx IF( SUM( ABS(apartf(ji,1:jpl)) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN nidx2 = nidx2 + 1 idxice2(nidx2) = idxice(ji) a_i_2d(nidx2,:) = a_i_2d(ji,:) ! adjust to new indices v_i_2d(nidx2,:) = v_i_2d(ji,:) ! adjust to new indices ato_i_1d(nidx2) = ato_i_1d(ji) ! adjust to new indices ENDIF END DO idxice(:) = idxice2(:) nidx = nidx2 ENDIF !-------------- ! Start ridging !-------------- IF( nidx > 0 ) THEN !----------- ! 2D <==> 1D !----------- ! fields used but not modified CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d(1:nidx), sss_m(:,:) ) CALL tab_2d_1d( nidx, idxice(1:nidx), sst_1d(1:nidx), sst_m(:,:) ) ! the following fields are modified in this routine !!CALL tab_2d_1d( nidx, idxice(1:nidx), ato_i_1d(1:nidx), ato_i(:,:) ) !!CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d(1:nidx,1:jpl), a_i(:,:,:) ) !!CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i (:,:,:) ) CALL tab_3d_2d( nidx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s (:,:,:) ) CALL tab_3d_2d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i(:,:,:) ) CALL tab_3d_2d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i (:,:,:) ) IF ( nn_pnd_scheme > 0 ) THEN CALL tab_3d_2d( nidx, idxice(1:nidx), a_ip_2d(1:nidx,1:jpl), a_ip(:,:,:) ) CALL tab_3d_2d( nidx, idxice(1:nidx), v_ip_2d(1:nidx,1:jpl), v_ip(:,:,:) ) ENDIF DO jl = 1, jpl DO jk = 1, nlay_s CALL tab_2d_1d( nidx, idxice(1:nidx), ze_s_2d(1:nidx,jk,jl), e_s(:,:,jk,jl) ) END DO DO jk = 1, nlay_i CALL tab_2d_1d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) ) END DO END DO CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_dyn_1d (1:nidx), sfx_dyn (:,:) ) CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_bri_1d (1:nidx), sfx_bri (:,:) ) CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_dyn_1d (1:nidx), wfx_dyn (:,:) ) CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_dyn_1d (1:nidx), hfx_dyn (:,:) ) CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_snw_dyn_1d(1:nidx), wfx_snw_dyn(:,:) ) IF ( nn_pnd_scheme > 0 ) THEN CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_pnd_1d(1:nidx), wfx_pnd(:,:) ) ENDIF !-----------------------------------------------------------------------------! ! 3) Ridging iteration !-----------------------------------------------------------------------------! niter = 1 iterate_ridging = 1 DO WHILE( iterate_ridging > 0 .AND. niter < nitermax ) CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d ) ! 3.2 Redistribute area, volume, and energy. !-----------------------------------------------------------------------------! CALL rdgrft_shift ! 3.4 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 ji = 1, nidx zfac = 1._wp - ( ato_i_1d(ji) + SUM( a_i_2d(ji,:) ) ) IF( ABS( zfac ) < epsi10 ) THEN closing_net(ji) = 0._wp opning (ji) = 0._wp ato_i_1d (ji) = MAX( 0._wp, 1._wp - SUM( a_i_2d(ji,:) ) ) ELSE iterate_ridging = 1 zdivu_adv (ji) = zfac * r1_rdtice closing_net(ji) = MAX( 0._wp, -zdivu_adv(ji) ) opning (ji) = MAX( 0._wp, zdivu_adv(ji) ) ENDIF END DO ! niter = niter + 1 IF( niter > nitermax ) CALL ctl_warn( 'icedyn_rdgrft: non-converging ridging scheme' ) ! END DO !----------- ! 1D <==> 2D !----------- CALL tab_1d_2d( nidx, idxice(1:nidx), ato_i_1d(1:nidx), ato_i(:,:) ) CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d(1:nidx,1:jpl), a_i(:,:,:) ) CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i (:,:,:) ) CALL tab_2d_3d( nidx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s (:,:,:) ) CALL tab_2d_3d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i(:,:,:) ) CALL tab_2d_3d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i (:,:,:) ) IF ( nn_pnd_scheme > 0 ) THEN CALL tab_2d_3d( nidx, idxice(1:nidx), a_ip_2d(1:nidx,1:jpl), a_ip(:,:,:) ) CALL tab_2d_3d( nidx, idxice(1:nidx), v_ip_2d(1:nidx,1:jpl), v_ip(:,:,:) ) ENDIF DO jl = 1, jpl DO jk = 1, nlay_s CALL tab_1d_2d( nidx, idxice(1:nidx), ze_s_2d(1:nidx,jk,jl), e_s(:,:,jk,jl) ) END DO DO jk = 1, nlay_i CALL tab_1d_2d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) ) END DO END DO CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_dyn_1d (1:nidx), sfx_dyn (:,:) ) CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_bri_1d (1:nidx), sfx_bri (:,:) ) CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_dyn_1d (1:nidx), wfx_dyn (:,:) ) CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_dyn_1d (1:nidx), hfx_dyn (:,:) ) CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_snw_dyn_1d(1:nidx), wfx_snw_dyn(:,:) ) IF ( nn_pnd_scheme > 0 ) THEN CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_pnd_1d(1:nidx), wfx_pnd(:,:) ) ENDIF ENDIF ! nidx > 0 CALL ice_var_agg( 1 ) ! controls IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation IF( ln_ctl ) CALL ice_prt3D ('icedyn_rdgrft') ! prints IF( nn_timing == 1 ) CALL timing_stop ('icedyn_rdgrft') ! timing ! END SUBROUTINE ice_dyn_rdgrft SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i ) !!------------------------------------------------------------------- !! *** ROUTINE rdgrft_prep *** !! !! ** Purpose : preparation for ridging calculations !! !! ** Method : Compute the thickness distribution of the ice and open water !! participating in ridging and of the resulting ridges. !!------------------------------------------------------------------- REAL(wp), DIMENSION(:) , INTENT(in) :: pato_i REAL(wp), DIMENSION(:,:), INTENT(in) :: pa_i, pv_i !! INTEGER :: ji, jl ! dummy loop indices REAL(wp) :: z1_gstar, z1_astar, zhmean, zfac ! local scalar REAL(wp), DIMENSION(jpij) :: zasum, z1_asum, zaksum ! sum of a_i+ato_i and reverse REAL(wp), DIMENSION(jpij,jpl) :: zhi ! ice thickness REAL(wp), DIMENSION(jpij,-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 ! ! Ice thickness needed for rafting WHERE( pa_i(1:nidx,:) > 0._wp ) ; zhi(1:nidx,:) = pv_i(1:nidx,:) / pa_i(1:nidx,:) ELSEWHERE ; zhi(1:nidx,:) = 0._wp END WHERE !----------------------------------------------------------------- ! 1) Participation function: a(h) = b(h).g(h) (apartf) !----------------------------------------------------------------- ! Compute the participation function apartf; 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 ! apartf(n) = total area lost due to ridging/closing ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar). ! ! The expressions for apartf are found by integrating b(h)g(h) between ! the category boundaries. ! apartf is always >= 0 and SUM(apartf(0:jpl))=1 !----------------------------------------------------------------- ! ! Compute total area of ice plus open water. ! This is in general not equal to one because of divergence during transport zasum(1:nidx) = pato_i(1:nidx) + SUM( pa_i(1:nidx,:), dim=2 ) ! WHERE( zasum(1:nidx) > 0._wp ) ; z1_asum(1:nidx) = 1._wp / zasum(1:nidx) ELSEWHERE ; z1_asum(1:nidx) = 0._wp END WHERE ! 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:nidx,-1) = 0._wp zGsum(1:nidx,0 ) = pato_i(1:nidx) * z1_asum(1:nidx) DO jl = 1, jpl zGsum(1:nidx,jl) = ( pato_i(1:nidx) + SUM( pa_i(1:nidx,1:jl), dim=2 ) ) * z1_asum(1:nidx) ! sum(1:jl) is ok (and not jpl) END DO ! IF( ln_partf_lin ) THEN !--- Linear formulation (Thorndike et al., 1975) DO jl = 0, jpl DO ji = 1, nidx IF ( zGsum(ji,jl) < rn_gstar ) THEN apartf(ji,jl) = z1_gstar * ( zGsum(ji,jl) - zGsum(ji,jl-1) ) * & & ( 2._wp - ( zGsum(ji,jl-1) + zGsum(ji,jl) ) * z1_gstar ) ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN apartf(ji,jl) = z1_gstar * ( rn_gstar - zGsum(ji,jl-1) ) * & & ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar ) * z1_gstar ) ELSE apartf(ji,jl) = 0._wp ENDIF END DO END DO ! ELSEIF( ln_partf_exp ) THEN !--- Exponential, more stable formulation (Lipscomb et al, 2007) ! zfac = 1._wp / ( 1._wp - EXP(-z1_astar) ) DO jl = -1, jpl DO ji = 1, nidx zGsum(ji,jl) = EXP( -zGsum(ji,jl) * z1_astar ) * zfac END DO END DO DO jl = 0, jpl DO ji = 1, nidx apartf(ji,jl) = zGsum(ji,jl-1) - zGsum(ji,jl) END DO END DO ! ENDIF ! !--- Ridging and rafting participation concentrations IF( ln_rafting .AND. ln_ridging ) THEN !- ridging & rafting DO jl = 1, jpl DO ji = 1, nidx aridge(ji,jl) = ( 1._wp + TANH ( rn_craft * ( zhi(ji,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jl) araft (ji,jl) = apartf(ji,jl) - aridge(ji,jl) END DO END DO ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN !- ridging alone DO jl = 1, jpl DO ji = 1, nidx aridge(ji,jl) = apartf(ji,jl) araft (ji,jl) = 0._wp END DO END DO ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN !- rafting alone DO jl = 1, jpl DO ji = 1, nidx aridge(ji,jl) = 0._wp araft (ji,jl) = apartf(ji,jl) END DO END DO ELSE !- no ridging & no rafting DO jl = 1, jpl DO ji = 1, nidx aridge(ji,jl) = 0._wp araft (ji,jl) = 0._wp END DO END DO 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 >= hrdg_hi_min*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. ! ! zaksum = 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 !----------------------------------------------------------------- zaksum(:) = apartf(:,0) zfac = 1._wp / hi_hrft ! Transfer function DO jl = 1, jpl !all categories have a specific transfer function DO ji = 1, nidx IF ( apartf(ji,jl) > 0._wp ) THEN zhmean = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min ) hrmin (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) ) hrmax (ji,jl) = 2._wp * zhmean - hrmin(ji,jl) hraft (ji,jl) = zhi(ji,jl) * zfac hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 ) ! ! Normalization factor : zaksum, ensures mass conservation zaksum(ji) = zaksum(ji) + aridge(ji,jl) * ( 1._wp - hi_hrdg(ji,jl) ) & & + araft (ji,jl) * ( 1._wp - hi_hrft ) ELSE hrmin (ji,jl) = 0._wp hrmax (ji,jl) = 0._wp hraft (ji,jl) = 0._wp hi_hrdg(ji,jl) = 1._wp ENDIF END DO END DO ! ! 3.1 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 WHERE( zaksum(1:nidx) > 0._wp ) ; closing_gross(1:nidx) = closing_net(1:nidx) / zaksum(1:nidx) ELSEWHERE ; closing_gross(1:nidx) = 0._wp END WHERE DO ji = 1, nidx ! 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. zfac = ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice IF ( zfac < 0._wp .AND. zfac > - pato_i(ji) ) THEN ! would lead to negative ato_i opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_rdtice ELSEIF( zfac > 0._wp .AND. zfac > ( zasum(ji) - pato_i(ji) ) ) THEN ! would lead to ato_i > asum opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_rdtice ENDIF 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 ji = 1, nidx zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice !! IF( zfac > pa_i(ji,jl) .AND. zfac > 0._wp ) THEN IF( zfac > pa_i(ji,jl) ) THEN closing_gross(ji) = closing_gross(ji) * pa_i(ji,jl) / zfac ENDIF END DO END DO ! END SUBROUTINE rdgrft_prep SUBROUTINE rdgrft_shift !!------------------------------------------------------------------- !! *** ROUTINE rdgrft_shift *** !! !! ** 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. !!------------------------------------------------------------------- ! INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices REAL(wp) :: hL, hR, farea ! left and right limits of integration and new area going to jl2 REAL(wp) :: vsw ! vol of water trapped into ridges REAL(wp) :: afrdg, afrft ! fraction of category area ridged/rafted REAL(wp) :: airdg1, oirdg1, aprdg1, virdg1, sirdg1 REAL(wp) :: airft1, oirft1, aprft1 REAL(wp), DIMENSION(jpij) :: airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, esrdg ! area etc of new ridges REAL(wp), DIMENSION(jpij) :: airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, esrft ! area etc of rafted ice ! REAL(wp), DIMENSION(jpij) :: ersw ! enth of water trapped into ridges REAL(wp), DIMENSION(jpij) :: zswitch, fvol ! new ridge volume going to jl2 ! REAL(wp), DIMENSION(jpij,jpl) :: z1_ai ! 1 / a ! REAL(wp), DIMENSION(jpij,nlay_i) :: eirft ! ice energy of rafting ice REAL(wp), DIMENSION(jpij,nlay_i) :: eirdg ! enth*volume of new ridges !!------------------------------------------------------------------- !------------------------------------------------------------------------------- ! 1) Compute change in open water area due to closing and opening. !------------------------------------------------------------------------------- DO ji = 1, nidx ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice ) END DO !----------------------------------------------------------------- ! 2) compute categories participating to ridging/rafting (jl1) !----------------------------------------------------------------- DO jl1 = 1, jpl CALL tab_2d_1d( nidx, idxice(1:nidx), sm_i_1d(1:nidx), sm_i(:,:,jl1) ) DO ji = 1, nidx !------------------------------------------------ ! 2.1) Identify grid cells with nonzero ridging !------------------------------------------------ IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN z1_ai(ji,jl1) = 1._wp / a_i_2d(ji,jl1) !-------------------------------------------------------------------- ! 2.2) Compute area of ridging ice (airdg1) and of new ridge (airdg2) !-------------------------------------------------------------------- airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice airft1 = araft (ji,jl1) * closing_gross(ji) * rdt_ice airdg2(ji) = airdg1 * hi_hrdg(ji,jl1) airft2(ji) = airft1 * hi_hrft !--------------------------------------------------------------- ! 2.3) Compute ridging /rafting fractions, make sure afrdg <=1 !--------------------------------------------------------------- afrdg = airdg1 * z1_ai(ji,jl1) afrft = airft1 * z1_ai(ji,jl1) ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg ersw(ji) = -rhoic * vsw * rcp * sst_1d(ji) ! clem: if sst>0, then ersw <0 (is that possible?) !--------------------------------------------------------------- ! 2.4) Compute ridging ice and new ridges (vi, vs, sm, oi, es, ei) !--------------------------------------------------------------- virdg1 = v_i_2d (ji,jl1) * afrdg virdg2(ji) = v_i_2d (ji,jl1) * afrdg * ( 1. + rn_porordg ) vsrdg(ji) = v_s_2d (ji,jl1) * afrdg sirdg1 = smv_i_2d(ji,jl1) * afrdg sirdg2(ji) = smv_i_2d(ji,jl1) * afrdg + vsw * sss_1d(ji) oirdg1 = oa_i_2d (ji,jl1) * afrdg oirdg2(ji) = oa_i_2d (ji,jl1) * afrdg * hi_hrdg(ji,jl1) esrdg(ji) = ze_s_2d (ji,1,jl1) * afrdg virft(ji) = v_i_2d (ji,jl1) * afrft vsrft(ji) = v_s_2d (ji,jl1) * afrft sirft(ji) = smv_i_2d(ji,jl1) * afrft oirft1 = oa_i_2d (ji,jl1) * afrft oirft2(ji) = oa_i_2d (ji,jl1) * afrft * hi_hrft esrft(ji) = ze_s_2d (ji,1,jl1) * afrft !MV MP 2016 IF ( nn_pnd_scheme > 0 ) THEN aprdg1 = a_ip_2d(ji,jl1) * afrdg aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1) vprdg (ji) = v_ip_2d(ji,jl1) * afrdg aprft1 = a_ip_2d(ji,jl1) * afrft aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft vprft (ji) = v_ip_2d(ji,jl1) * afrft ENDIF ! END MV MP 2016 !----------------------------------------------------------------- ! 2.5) Ice-ocean exchanges !----------------------------------------------------------------- wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoic * r1_rdtice ! increase in ice volume due to seawater frozen in voids sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoic * r1_rdtice hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_rdtice ! > 0 [W.m-2] ! Put 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. wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhosn * vsrdg(ji) * ( 1._wp - rn_fsnwrdg ) & ! fresh water source for ocean & + rhosn * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji) * ( 1._wp - rn_fsnwrdg ) & ! heat sink for ocean (<0, W.m-2) & - esrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice ! Put the melt pond water into the ocean !------------------------------------------ IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) THEN wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhofw * vprdg(ji) * ( 1._wp - rn_fpndrdg ) & ! fresh water source for ocean & + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice ENDIF ! virtual salt flux to keep salinity constant !------------------------------------------ IF( nn_icesal /= 2 ) THEN sirdg2(ji) = sirdg2(ji) - vsw * ( sss_1d(ji) - sm_i_1d(ji) ) ! ridge salinity = sm_i sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoic * r1_rdtice & ! put back sss_m into the ocean & - sm_i_1d(ji) * vsw * rhoic * r1_rdtice ! and get sm_i from the ocean ENDIF !------------------------------------------ ! 2.6) update jl1 (removing ridged/rafted area) !------------------------------------------ a_i_2d (ji,jl1) = a_i_2d (ji,jl1) - airdg1 - airft1 v_i_2d (ji,jl1) = v_i_2d (ji,jl1) - virdg1 - virft(ji) v_s_2d (ji,jl1) = v_s_2d (ji,jl1) - vsrdg(ji) - vsrft(ji) smv_i_2d(ji,jl1) = smv_i_2d(ji,jl1) - sirdg1 - sirft(ji) oa_i_2d (ji,jl1) = oa_i_2d (ji,jl1) - oirdg1 - oirft1 ! MV MP 2016 IF ( nn_pnd_scheme > 0 ) THEN a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1 - aprft1 v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji) ENDIF ! END MV MP 2016 ze_s_2d(ji,1,jl1) = ze_s_2d(ji,1,jl1) - esrdg (ji) - esrft (ji) ENDIF END DO ! ji ! special loop for e_i because of layers jk DO jk = 1, nlay_i DO ji = 1, nidx IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN ! Compute ridging /rafting fractions afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji,jl1) afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji,jl1) ! Compute ridging ice and new ridges for ei eirdg(ji,jk) = ze_i_2d (ji,jk,jl1) * afrdg + ersw(ji) * r1_nlay_i eirft(ji,jk) = ze_i_2d (ji,jk,jl1) * afrft ! Update jl1 ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) ENDIF END DO END DO !------------------------------------------------------------------------------- ! 3) Add area, volume, and energy of new ridge to each category jl2 !------------------------------------------------------------------------------- DO jl2 = 1, jpl ! DO ji = 1, nidx IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN ! Compute the fraction of ridged ice area and volume going to thickness category jl2. IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN !----------------------------------------------------------------- ! 2.8 Compute quantities used to apportion ice among categories ! in the n2 loop below !----------------------------------------------------------------- hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) ) hR = MIN( hrmax(ji,jl1), hi_max(jl2) ) farea = ( hR - hL ) / ( hrmax(ji,jl1) - hrmin(ji,jl1) ) fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) ) ELSE farea = 0._wp fvol(ji) = 0._wp ENDIF ! Compute the fraction of rafted ice area and volume going to thickness category jl2 !!gm see above IF( hraft(ji) <= hi_max(jl2) .AND. hraft(ji) > hi_max(jl2-1) ) THEN IF( hi_max(jl2-1) < hraft(ji,jl1) .AND. hraft(ji,jl1) <= hi_max(jl2) ) THEN ; zswitch(ji) = 1._wp ELSE ; zswitch(ji) = 0._wp ENDIF ! a_i_2d (ji,jl2) = a_i_2d (ji,jl2) + ( airdg2(ji) * farea + airft2(ji) * zswitch(ji) ) oa_i_2d (ji,jl2) = oa_i_2d (ji,jl2) + ( oirdg2(ji) * farea + oirft2(ji) * zswitch(ji) ) v_i_2d (ji,jl2) = v_i_2d (ji,jl2) + ( virdg2(ji) * fvol(ji) + virft (ji) * zswitch(ji) ) smv_i_2d(ji,jl2) = smv_i_2d(ji,jl2) + ( sirdg2(ji) * fvol(ji) + sirft (ji) * zswitch(ji) ) v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji) + & & vsrft (ji) * rn_fsnwrft * zswitch(ji) ) ! MV MP 2016 IF ( nn_pnd_scheme > 0 ) THEN v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + ( vprdg (ji) * rn_fpndrdg * fvol (ji) & & + vprft (ji) * rn_fpndrft * zswitch(ji) ) a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + ( aprdg2(ji) * rn_fpndrdg * farea & & + aprft2(ji) * rn_fpndrft * zswitch(ji) ) ENDIF ! END MV MP 2016 ze_s_2d(ji,1,jl2) = ze_s_2d(ji,1,jl2) + ( esrdg (ji) * rn_fsnwrdg * fvol(ji) + & & esrft (ji) * rn_fsnwrft * zswitch(ji) ) ENDIF END DO DO jk = 1, nlay_i DO ji = 1, nidx IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) & & ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji) END DO END DO ! END DO ! jl2 ! END DO ! jl1 ! ! In case ridging/rafting lead to very small negative values (sometimes it happens) WHERE( a_i_2d(1:nidx,:) < 0._wp ) a_i_2d(1:nidx,:) = 0._wp WHERE( v_i_2d(1:nidx,:) < 0._wp ) v_i_2d(1:nidx,:) = 0._wp ! END SUBROUTINE rdgrft_shift SUBROUTINE ice_strength !!---------------------------------------------------------------------- !! *** ROUTINE ice_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 !!---------------------------------------------------------------------- ! !--------------------------------------------------! IF( ln_str_R75 ) THEN ! Ice strength => Rothrock (1975) method ! ! !--------------------------------------------------! ! !--------------------------------------------------! ELSEIF( ln_str_H79 ) THEN ! Ice strength => Hibler (1979) method ! ! !--------------------------------------------------! strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) ! ismooth = 1 ! ENDIF ! !--------------------------------------------------! SELECT CASE( ismooth ) ! Smoothing ice strength ! ! !--------------------------------------------------! CASE( 1 ) !--- Spatial smoothing DO jj = 2, jpjm1 DO ji = 2, jpim1 IF ( SUM( a_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 ( SUM( a_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_strength SUBROUTINE ice_dyn_rdgrft_init !!------------------------------------------------------------------- !! *** ROUTINE ice_dyn_rdgrft_init *** !! !! ** Purpose : Physical constants and parameters linked !! to the mechanical ice redistribution !! !! ** Method : Read the namdyn_rdgrft namelist !! and check the parameters values !! called at the first timestep (nit000) !! !! ** input : Namelist namdyn_rdgrft !!------------------------------------------------------------------- INTEGER :: ios ! Local integer output status for namelist read !! NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, & & ln_str_R75, rn_perdg, & & rn_csrdg , & & 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, namdyn_rdgrft, IOSTAT = ios, ERR = 901) 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist', lwp ) ! REWIND( numnam_ice_cfg ) ! Namelist namdyn_rdgrft in configuration namelist : Ice mechanical ice redistribution READ ( numnam_ice_cfg, namdyn_rdgrft, IOSTAT = ios, ERR = 902 ) 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_rdgrft in configuration namelist', lwp ) IF(lwm) WRITE ( numoni, namdyn_rdgrft ) ! IF (lwp) THEN ! control print WRITE(numout,*) WRITE(numout,*) 'ice_dyn_rdgrft_init: ice parameters for ridging/rafting ' WRITE(numout,*) '~~~~~~~~~~~~~~~~~~' WRITE(numout,*) ' Namelist namdyn_rdgrft:' WRITE(numout,*) ' ice strength parameterization Hibler (1979) ln_str_H79 = ', ln_str_H79 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_R75 = ', ln_str_R75 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_csrdg = ', rn_csrdg 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_H79 .AND. ln_str_R75 ) .OR. ( .NOT.ln_str_H79 .AND. .NOT.ln_str_R75 ) ) THEN CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one formulation for ice strength (ln_str_H79 or ln_str_R75)' ) ENDIF ! IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' ) ENDIF ! ! allocate tke arrays IF( ice_dyn_rdgrft_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' ) ! END SUBROUTINE ice_dyn_rdgrft_init #else !!---------------------------------------------------------------------- !! Default option Empty module NO ESIM sea-ice model !!---------------------------------------------------------------------- #endif !!====================================================================== END MODULE icedyn_rdgrft