MODULE limitd_me #if defined key_lim3 !!---------------------------------------------------------------------- !! 'key_lim3' : LIM3 sea-ice model !!---------------------------------------------------------------------- !! !!====================================================================== !! *** MODULE limitd_me *** !! Mechanical impact on ice thickness distribution !! computation of changes in g(h) !!====================================================================== !! !!---------------------------------------------------------------------- !! * Modules used USE dom_ice USE par_oce ! ocean parameters USE dom_oce USE lbclnk USE phycst ! physical constants (ocean directory) USE ice_oce ! ice variables USE sbc_oce ! Surface boundary condition: ocean fields USE thd_ice USE limistate USE in_out_manager USE ice USE par_ice USE limthd_lac USE limvar USE iceini USE limcons USE prtctl ! Print control USE lib_mpp IMPLICIT NONE PRIVATE !! * Routine accessibility PUBLIC lim_itd_me ! called by ice_stp PUBLIC lim_itd_me_icestrength PUBLIC lim_itd_me_ridgeprep PUBLIC lim_itd_me_ridgeshift PUBLIC lim_itd_me_asumr PUBLIC lim_itd_me_init PUBLIC lim_itd_me_zapsmall !! * Module variables REAL(wp) :: & ! constant values epsi20 = 1e-20 , & epsi13 = 1e-13 , & epsi11 = 1e-11 , & zzero = 0.e0 , & zone = 1.e0 !----------------------------------------------------------------------- ! Variables shared among ridging subroutines !----------------------------------------------------------------------- REAL(wp), DIMENSION (jpi,jpj) :: & asum , & ! sum of total ice and open water area aksum ! ratio of area removed to area ridged REAL(wp), DIMENSION(jpi,jpj,0:jpl) :: & athorn ! participation function; fraction of ridging/ ! closing associated w/ category n REAL(wp), DIMENSION(jpi,jpj,jpl) :: & hrmin , & ! minimum ridge thickness hrmax , & ! maximum ridge thickness hraft , & ! thickness of rafted ice krdg , & ! mean ridge thickness/thickness of ridging ice aridge , & ! participating ice ridging araft ! participating ice rafting REAL(wp), PARAMETER :: & krdgmin = 1.1, & ! min ridge thickness multiplier kraft = 2.0 ! rafting multipliyer REAL(wp) :: & Cp ! !----------------------------------------------------------------------- ! Ridging diagnostic arrays for history files !----------------------------------------------------------------------- ! REAL (wp), DIMENSION(jpi,jpj) :: & 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) !!---------------------------------------------------------------------- !! LIM-@ 3.0, UCL-ASTR (2008) !! (c) UCL-ASTR and Martin Vancoppenolle !!---------------------------------------------------------------------- CONTAINS !!-----------------------------------------------------------------------------! !!-----------------------------------------------------------------------------! SUBROUTINE lim_itd_me ! (subroutine 1/6) !!---------------------------------------------------------------------! !! *** ROUTINE lim_itd_me *** !! ** Purpose : !! This routine computes the mechanical redistribution !! of ice thickness !! !! ** Method : a very simple method :-) !! !! ** Arguments : !! kideb , kiut : Starting and ending points on which the !! the computation is applied !! !! ** Inputs / Ouputs : (global commons) !! !! ** External : !! !! ** 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 : There are a lot of references and can be difficult / !! boring to read !! !! Flato, G. M., and W. D. Hibler III, 1995: Ridging and strength !! in modeling the thickness distribution of Arctic sea ice, !! J. Geophys. Res., 100, 18,611-18,626. !! !! Hibler, W. D. III, 1980: Modeling a variable thickness sea ice !! cover, Mon. Wea. Rev., 108, 1943-1973, 1980. !! !! Rothrock, D. A., 1975: The energetics of the plastic deformation of !! pack ice by ridging, J. Geophys. Res., 80, 4514-4519. !! !! Thorndike, A. S., D. A. Rothrock, G. A. Maykut, and R. Colony, !! 1975: The thickness distribution of sea ice, J. Geophys. Res., !! 80, 4501-4513. !! !! Bitz et al., JGR 2001 !! !! Amundrud and Melling, JGR 2005 !! !! Babko et al., JGR 2002 !! !! ** History : !! This routine is based on CICE code !! and authors William H. Lipscomb, !! and Elizabeth C. Hunke, LANL !! are gratefully acknowledged !! !! (02-2006) Martin Vancoppenolle, UCL-ASTR !! !!--------------------------------------------------------------------! !! * Arguments !! * Local variables INTEGER :: ji, & ! spatial dummy loop index jj, & ! spatial dummy loop index jk, & ! vertical layering dummy loop index jl, & ! ice category dummy loop index niter, & ! iteration counter nitermax = 20 ! max number of ridging iterations REAL(wp) :: & ! constant values zeps = 1.0e-10, & epsi10 = 1.0e-10, & epsi06 = 1.0e-6 REAL(wp), DIMENSION(jpi,jpj) :: & 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 closing_gross, & ! rate at which area removed, not counting ! area of new ridges msnow_mlt , & ! mass of snow added to ocean (kg m-2) esnow_mlt ! energy needed to melt snow in ocean (J m-2) REAL(wp) :: & w1, & ! temporary variable tmpfac, & ! factor by which opening/closing rates are cut dti ! 1 / dt LOGICAL :: & asum_error ! flag for asum .ne. 1 INTEGER :: iterate_ridging ! if true, repeat the ridging REAL(wp) :: & big = 1.0e8 REAL (wp), DIMENSION(jpi,jpj) :: & ! vt_i_init, vt_i_final ! ice volume summed over categories CHARACTER (len = 15) :: fieldid !!-- End of declarations !-----------------------------------------------------------------------------! IF( numit == nstart ) CALL lim_itd_me_init ! Initialization (first time-step only) IF(ln_ctl) THEN CALL prt_ctl(tab2d_1=ato_i , clinfo1=' lim_itd_me: ato_i : ', tab2d_2=at_i , clinfo2=' at_i : ') CALL prt_ctl(tab2d_1=divu_i, clinfo1=' lim_itd_me: divu_i : ', tab2d_2=delta_i, clinfo2=' delta_i : ') ENDIF !-----------------------------------------------------------------------------! ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons !-----------------------------------------------------------------------------! ! Set hi_max(ncat) to a big value to ensure that all ridged ice ! is thinner than hi_max(ncat). hi_max(jpl) = 999.99 Cp = 0.5* grav * (rau0-rhoic)*rhoic/rau0 ! proport const for PE CALL lim_itd_me_ridgeprep ! prepare ridging ! conservation check IF ( con_i) CALL lim_column_sum (jpl, v_i, vt_i_init) ! Initialize arrays. DO jj = 1, jpj DO ji = 1, jpi msnow_mlt(ji,jj) = 0.0 esnow_mlt(ji,jj) = 0.0 dardg1dt(ji,jj) = 0.0 dardg2dt(ji,jj) = 0.0 dvirdgdt(ji,jj) = 0.0 opening (ji,jj) = 0.0 !-----------------------------------------------------------------------------! ! 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) = & Cs*0.5*(Delta_i(ji,jj)-ABS(divu_i(ji,jj))) - MIN(divu_i(ji,jj),0.0) ! 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.0-asum(ji,jj)) / rdt_ice ! asum found in ridgeprep IF (divu_adv(ji,jj) .LT. 0.0) & 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 ) DO jj = 1, jpj DO ji = 1, jpi ! 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(ji,jj) = closing_net(ji,jj) / aksum(ji,jj) ! 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. IF ( ato_i(ji,jj) .GT. epsi11 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice IF ( w1 .GT. ato_i(ji,jj)) THEN tmpfac = ato_i(ji,jj) / w1 closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac opning(ji,jj) = opning(ji,jj) * tmpfac ENDIF !w1 ENDIF !at0i and athorn END DO ! ji END DO ! jj ! 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 IF ( a_i(ji,jj,jl) .GT. epsi11 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice IF ( w1 .GT. a_i(ji,jj,jl) ) THEN tmpfac = a_i(ji,jj,jl) / w1 closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac opning(ji,jj) = opning(ji,jj) * tmpfac ENDIF ENDIF END DO !ji END DO ! jj END DO !jl ! 3.3 Redistribute area, volume, and energy. !-----------------------------------------------------------------------------! CALL lim_itd_me_ridgeshift (opning, closing_gross, & msnow_mlt, esnow_mlt) ! 3.4 Compute total area of ice plus open water after ridging. !-----------------------------------------------------------------------------! CALL lim_itd_me_asumr ! 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.0) .LT. epsi11) THEN closing_net(ji,jj) = 0.0 opning(ji,jj) = 0.0 ELSE iterate_ridging = 1 divu_adv(ji,jj) = (1.0 - asum(ji,jj)) / rdt_ice closing_net(ji,jj) = MAX(0.0, -divu_adv(ji,jj)) opning(ji,jj) = MAX(0.0, 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 IF (niter .GT. nitermax) THEN WRITE(numout,*) ' ALERTE : non-converging ridging scheme ' WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging ENDIF CALL lim_itd_me_ridgeprep ENDIF END DO !! on the do while over iter !-----------------------------------------------------------------------------! ! 4) Ridging diagnostics !-----------------------------------------------------------------------------! ! Convert ridging rate diagnostics to correct units. ! Update fresh water and heat fluxes due to snow melt. dti = 1.0/rdt_ice asum_error = .false. DO jj = 1, jpj DO ji = 1, jpi IF (ABS(asum(ji,jj) - 1.0) .GT. epsi11) asum_error = .true. dardg1dt(ji,jj) = dardg1dt(ji,jj) * dti dardg2dt(ji,jj) = dardg2dt(ji,jj) * dti dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * dti opening (ji,jj) = opening (ji,jj) * dti !-----------------------------------------------------------------------------! ! 5) Heat, salt and freshwater fluxes !-----------------------------------------------------------------------------! ! fresh water source for ocean fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj)*dti ! heat sink for ocean fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj)*dti END DO END DO ! Check if there is a ridging error DO jj = 1, jpj DO ji = 1, jpi IF (ABS(asum(ji,jj) - 1.0) .GT. epsi11) THEN ! there is a bug WRITE(numout,*) ' ' WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) WRITE(numout,*) ' limitd_me ' WRITE(numout,*) ' POINT : ', ji, jj WRITE(numout,*) ' jpl, a_i, athorn ' WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0) DO jl = 1, jpl WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl) END DO ENDIF ! asum END DO !ji END DO !jj ! Conservation check IF ( con_i ) THEN CALL lim_column_sum (jpl, v_i, vt_i_final) fieldid = ' v_i : limitd_me ' CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) ENDIF !-----------------------------------------------------------------------------! ! 6) Updating state variables and trend terms !-----------------------------------------------------------------------------! CALL lim_var_glo2eqv CALL lim_itd_me_zapsmall !----------------- ! Trend terms !----------------- d_u_ice_dyn(:,:) = u_ice(:,:) - old_u_ice(:,:) d_v_ice_dyn(:,:) = v_ice(:,:) - old_v_ice(:,:) d_a_i_trp(:,:,:) = a_i(:,:,:) - old_a_i(:,:,:) d_v_s_trp(:,:,:) = v_s(:,:,:) - old_v_s(:,:,:) d_v_i_trp(:,:,:) = v_i(:,:,:) - old_v_i(:,:,:) d_e_s_trp(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:) d_e_i_trp(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:) d_oa_i_trp(:,:,:) = oa_i(:,:,:) - old_oa_i(:,:,:) d_smv_i_trp(:,:,:) = 0.0 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) IF(ln_ctl) THEN ! Control print CALL prt_ctl_info(' ') CALL prt_ctl_info(' - Cell values : ') CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_me : cell area :') CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me : at_i :') CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me : vt_i :') CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_me : vt_s :') DO jl = 1, jpl CALL prt_ctl_info(' ') CALL prt_ctl_info(' - Category : ', ivar1=jl) CALL prt_ctl_info(' ~~~~~~~~~~') CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_me : a_i : ') CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_me : ht_i : ') CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_me : ht_s : ') CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_me : v_i : ') CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_me : v_s : ') CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_me : e_s : ') CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_me : t_su : ') CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_me : t_snow : ') CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_me : sm_i : ') CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_me : smv_i : ') DO jk = 1, nlay_i CALL prt_ctl_info(' ') CALL prt_ctl_info(' - Layer : ', ivar1=jk) CALL prt_ctl_info(' ~~~~~~~') CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : t_i : ') CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : e_i : ') END DO END DO ENDIF !-------------------------! ! Back to initial values !-------------------------! ! update of fields will be made later in lim update u_ice(:,:) = old_u_ice(:,:) v_ice(:,:) = old_v_ice(:,:) a_i(:,:,:) = old_a_i(:,:,:) v_s(:,:,:) = old_v_s(:,:,:) v_i(:,:,:) = old_v_i(:,:,:) e_s(:,:,:,:) = old_e_s(:,:,:,:) e_i(:,:,:,:) = old_e_i(:,:,:,:) oa_i(:,:,:) = old_oa_i(:,:,:) IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & smv_i(:,:,:) = old_smv_i(:,:,:) !----------------------------------------------------! ! Advection of ice in a free cell, newly ridged ice !----------------------------------------------------! ! to allow for thermodynamics to melt new ice ! we immediately advect ice in free cells ! heat content has to be corrected before ice volume DO jl = 1, jpl DO jk = 1, nlay_i DO jj = 1, jpj DO ji = 1, jpi IF ( ( old_v_i(ji,jj,jl) .LT. epsi06 ) .AND. & ( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN old_e_i(ji,jj,jk,jl) = d_e_i_trp(ji,jj,jk,jl) d_e_i_trp(ji,jj,jk,jl) = 0.0 ENDIF END DO END DO END DO END DO DO jl = 1, jpl DO jj = 1, jpj DO ji = 1, jpi IF ( ( old_v_i(ji,jj,jl) .LT. epsi06 ) .AND. & ( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN old_v_i(ji,jj,jl) = d_v_i_trp(ji,jj,jl) d_v_i_trp(ji,jj,jl) = 0.0 old_a_i(ji,jj,jl) = d_a_i_trp(ji,jj,jl) d_a_i_trp(ji,jj,jl) = 0.0 old_v_s(ji,jj,jl) = d_v_s_trp(ji,jj,jl) d_v_s_trp(ji,jj,jl) = 0.0 old_e_s(ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl) d_e_s_trp(ji,jj,1,jl) = 0.0 old_oa_i(ji,jj,jl) = d_oa_i_trp(ji,jj,jl) d_oa_i_trp(ji,jj,jl) = 0.0 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) d_smv_i_trp(ji,jj,jl) = 0.0 ENDIF END DO END DO END DO END SUBROUTINE lim_itd_me !=============================================================================== SUBROUTINE lim_itd_me_icestrength (kstrngth) ! (subroutine 2/6) !!---------------------------------------------------------------------- !! *** ROUTINE lim_itd_me_icestrength *** !! ** Purpose : !! This routine 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 !! !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using) !! !! ** External : !! !! ** References : !! !!---------------------------------------------------------------------- !! * Arguments INTEGER, INTENT(in) :: & kstrngth ! = 1 for Rothrock formulation, 0 for Hibler (1979) INTEGER :: & ji,jj, & !: horizontal indices jl, & !: thickness category index ksmooth, & !: smoothing the resistance to deformation numts_rm !: number of time steps for the P smoothing REAL(wp) :: & hi, & !: ice thickness (m) zw1, & !: temporary variable zp, & !: temporary ice strength zdummy REAL(wp), DIMENSION(jpi,jpj) :: & zworka !: temporary array used here !------------------------------------------------------------------------------! ! 1) Initialize !------------------------------------------------------------------------------! strength(:,:) = 0.0 !------------------------------------------------------------------------------! ! 2) Compute thickness distribution of ridging and ridged ice !------------------------------------------------------------------------------! CALL lim_itd_me_ridgeprep !------------------------------------------------------------------------------! ! 3) Rothrock(1975)'s method !------------------------------------------------------------------------------! IF (kstrngth == 1) then DO jl = 1, jpl DO jj= 1, jpj DO ji = 1, jpi IF( ( a_i(ji,jj,jl) .GT. epsi11 ) & .AND. ( athorn(ji,jj,jl) .GT. 0.0 ) ) THEN hi = v_i(ji,jj,jl) / a_i(ji,jj,jl) !---------------------------- ! PE loss from deforming ice !---------------------------- strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * & hi * hi !-------------------------- ! PE gain from rafting ice !-------------------------- strength(ji,jj) = strength(ji,jj) + 2.0 * araft(ji,jj,jl) & * hi * hi !---------------------------- ! PE gain from ridging ice !---------------------------- strength(ji,jj) = strength(ji,jj) & + aridge(ji,jj,jl)/krdg(ji,jj,jl) & * 1.0/3.0 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) & / (hrmax(ji,jj,jl)-hrmin(ji,jj,jl)) ENDIF ! aicen > epsi11 END DO ! ji END DO !jj END DO !jl DO jj = 1, jpj DO ji = 1, jpi strength(ji,jj) = Cf * Cp * strength(ji,jj) / aksum(ji,jj) ! Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) ! Cf accounts for frictional dissipation END DO ! j END DO ! i ksmooth = 1 !------------------------------------------------------------------------------! ! 4) Hibler (1979)' method !------------------------------------------------------------------------------! ELSE ! kstrngth ne 1: Hibler (1979) form DO jj = 1, jpj DO ji = 1, jpi strength(ji,jj) = Pstar*vt_i(ji,jj)*exp(-C_rhg*(1.0-at_i(ji,jj))) END DO ! j END DO ! i ksmooth = 1 ENDIF ! kstrngth ! !------------------------------------------------------------------------------! ! 5) Impact of brine volume !------------------------------------------------------------------------------! ! CAN BE REMOVED ! IF ( brinstren_swi .EQ. 1 ) THEN DO jj = 1, jpj DO ji = 1, jpi IF ( bv_i(ji,jj) .GT. 0.0 ) THEN zdummy = MIN ( bv_i(ji,jj), 0.10 ) * MIN( bv_i(ji,jj), 0.10 ) ELSE zdummy = 0.0 ENDIF strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0))) END DO ! j END DO ! i ENDIF ! !------------------------------------------------------------------------------! ! 6) Smoothing ice strength !------------------------------------------------------------------------------! ! !------------------- ! Spatial smoothing !------------------- IF ( ksmooth .EQ. 1 ) THEN CALL lbc_lnk( strength, 'T', 1. ) DO jj = 2, jpj - 1 DO ji = 2, jpi - 1 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi11) THEN ! ice is ! present zworka(ji,jj) = 4.0 * strength(ji,jj) & + strength(ji-1,jj) * tms(ji-1,jj) & + strength(ji+1,jj) * tms(ji+1,jj) & + strength(ji,jj-1) * tms(ji,jj-1) & + strength(ji,jj+1) * tms(ji,jj+1) zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj) & + tms(ji,jj-1) + tms(ji,jj+1) zworka(ji,jj) = zworka(ji,jj) / zw1 ELSE zworka(ji,jj) = 0.0 ENDIF END DO END DO DO jj = 2, jpj - 1 DO ji = 2, jpi - 1 strength(ji,jj) = zworka(ji,jj) END DO END DO CALL lbc_lnk( strength, 'T', 1. ) ENDIF ! ksmooth !-------------------- ! Temporal smoothing !-------------------- IF ( numit .EQ. nit000 + nn_fsbc - 1 ) THEN strp1(:,:) = 0.0 strp2(:,:) = 0.0 ENDIF IF ( ksmooth .EQ. 2 ) THEN CALL lbc_lnk( strength, 'T', 1. ) DO jj = 1, jpj - 1 DO ji = 1, jpi - 1 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi11) THEN ! ice is ! present numts_rm = 1 ! number of time steps for the running mean IF ( strp1(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1 IF ( strp2(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1 zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / & numts_rm strp2(ji,jj) = strp1(ji,jj) strp1(ji,jj) = strength(ji,jj) strength(ji,jj) = zp ENDIF END DO END DO ENDIF ! ksmooth ! Boundary conditions CALL lbc_lnk( strength, 'T', 1. ) END SUBROUTINE lim_itd_me_icestrength !=============================================================================== SUBROUTINE lim_itd_me_ridgeprep !(subroutine 3/6) !!---------------------------------------------------------------------! !! *** ROUTINE lim_itd_me_ridgeprep *** !! ** 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. !! !! ** Arguments : !! !! ** External : !! !!---------------------------------------------------------------------! !! * Arguments INTEGER :: & ji,jj, & ! horizontal indices jl, & ! thickness category index krdg_index ! which participation function using REAL(wp) :: & Gstari, & ! = 1.0/Gstar astari ! = 1.0/astar REAL(wp), DIMENSION(jpi,jpj,-1:jpl) :: & Gsum ! Gsum(n) = sum of areas in categories 0 to n REAL(wp) :: & hi, & ! ice thickness for each cat (m) hrmean ! mean ridge thickness (m) REAL(wp), DIMENSION(jpi,jpj) :: & zworka ! temporary array used here REAL(wp) :: & zdummy, & epsi06 = 1.0e-6 !------------------------------------------------------------------------------! Gstari = 1.0/Gstar astari = 1.0/astar aksum(:,:) = 0.0 athorn(:,:,:) = 0.0 aridge(:,:,:) = 0.0 araft (:,:,:) = 0.0 hrmin(:,:,:) = 0.0 hrmax(:,:,:) = 0.0 hraft(:,:,:) = 0.0 krdg (:,:,:) = 1.0 ! ! Zero out categories with very small areas CALL lim_itd_me_zapsmall !------------------------------------------------------------------------------! ! 1) Participation function !------------------------------------------------------------------------------! ! Compute total area of ice plus open water. CALL lim_itd_me_asumr ! This is in general not equal to one ! because of divergence during transport ! Compute cumulative thickness distribution function ! Compute the cumulative thickness distribution function Gsum, ! where Gsum(n) is the fractional area in categories 0 to n. ! initial value (in h = 0) equals open water area Gsum(:,:,-1) = 0.0 DO jj = 1, jpj DO ji = 1, jpi IF (ato_i(ji,jj) .GT. epsi11) THEN Gsum(ji,jj,0) = ato_i(ji,jj) ELSE Gsum(ji,jj,0) = 0.0 ENDIF END DO END DO ! for each value of h, you have to add ice concentration then DO jl = 1, jpl DO jj = 1, jpj DO ji = 1, jpi IF ( a_i(ji,jj,jl) .GT. epsi11 ) THEN Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl) ELSE Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) ENDIF END DO END DO END DO ! Normalize the cumulative distribution to 1 DO jj = 1, jpj DO ji = 1, jpi zworka(ji,jj) = 1.0 / Gsum(ji,jj,jpl) END DO END DO DO jl = 0, jpl DO jj = 1, jpj DO ji = 1, jpi Gsum(ji,jj,jl) = Gsum(ji,jj,jl) * zworka(ji,jj) END DO END DO 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. !----------------------------------------------------------------- krdg_index = 1 IF ( krdg_index .EQ. 0 ) THEN !--- Linear formulation (Thorndike et al., 1975) DO jl = 0, ice_cat_bounds(1,2) ! only undeformed ice participates DO jj = 1, jpj DO ji = 1, jpi IF (Gsum(ji,jj,jl) < Gstar) THEN athorn(ji,jj,jl) = Gstari * (Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * & (2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari) ELSEIF (Gsum(ji,jj,jl-1) < Gstar) THEN athorn(ji,jj,jl) = Gstari * (Gstar-Gsum(ji,jj,jl-1)) * & (2.0 - (Gsum(ji,jj,jl-1)+Gstar)*Gstari) ELSE athorn(ji,jj,jl) = 0.0 ENDIF END DO ! ji END DO ! jj END DO ! jl ELSE ! krdg_index = 1 !--- Exponential, more stable formulation (Lipscomb et al, 2007) ! precompute exponential terms using Gsum as a work array zdummy = 1.0 / (1.0-EXP(-astari)) DO jl = -1, jpl DO jj = 1, jpj DO ji = 1, jpi Gsum(ji,jj,jl) = EXP(-Gsum(ji,jj,jl)*astari)*zdummy END DO !ji END DO !jj END DO !jl ! compute athorn DO jl = 0, ice_cat_bounds(1,2) DO jj = 1, jpj DO ji = 1, jpi athorn(ji,jj,jl) = Gsum(ji,jj,jl-1) - Gsum(ji,jj,jl) END DO !ji END DO ! jj END DO !jl ENDIF ! krdg_index ! Ridging and rafting ice participation functions IF ( raftswi .EQ. 1 ) THEN DO jl = 1, jpl DO jj = 1, jpj DO ji = 1, jpi IF ( athorn(ji,jj,jl) .GT. 0.0 ) THEN aridge(ji,jj,jl) = ( TANH ( Craft * ( ht_i(ji,jj,jl) - & hparmeter ) ) + 1.0 ) / 2.0 * & athorn(ji,jj,jl) araft (ji,jj,jl) = ( TANH ( - Craft * ( ht_i(ji,jj,jl) - & hparmeter ) ) + 1.0 ) / 2.0 * & athorn(ji,jj,jl) IF ( araft(ji,jj,jl) .LT. epsi06 ) araft(ji,jj,jl) = 0.0 aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0) ENDIF ! athorn END DO ! ji END DO ! jj END DO ! jl ELSE ! raftswi = 0 DO jl = 1, jpl DO jj = 1, jpj DO ji = 1, jpi aridge(ji,jj,jl) = 1.0*athorn(ji,jj,jl) END DO END DO END DO ENDIF IF ( raftswi .EQ. 1 ) THEN IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi11 ) THEN DO jl = 1, jpl DO jj = 1, jpj DO ji = 1, jpi IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. & epsi11 ) THEN WRITE(numout,*) ' ALERTE 96 : wrong participation function ... ' WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl WRITE(numout,*) ' lat, lon : ', gphit(ji,jj), glamt(ji,jj) WRITE(numout,*) ' aridge : ', aridge(ji,jj,1:jpl) WRITE(numout,*) ' araft : ', araft(ji,jj,1:jpl) WRITE(numout,*) ' athorn : ', athorn(ji,jj,1:jpl) ENDIF END DO END DO END DO ENDIF 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, hrmean, 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 <= (hrmean + hi)/2. ! ! The maximum ridging thickness, hrmax, is determined by ! hrmean 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 !----------------------------------------------------------------- ! Transfer function DO jl = 1, jpl !all categories have a specific transfer function DO jj = 1, jpj DO ji = 1, jpi IF (a_i(ji,jj,jl) .GT. epsi11 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN hi = v_i(ji,jj,jl) / a_i(ji,jj,jl) hrmean = MAX(SQRT(Hstar*hi), hi*krdgmin) hrmin(ji,jj,jl) = MIN(2.0*hi, 0.5*(hrmean + hi)) hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl) hraft(ji,jj,jl) = kraft*hi krdg(ji,jj,jl) = hrmean / hi ELSE hraft(ji,jj,jl) = 0.0 hrmin(ji,jj,jl) = 0.0 hrmax(ji,jj,jl) = 0.0 krdg (ji,jj,jl) = 1.0 ENDIF END DO ! ji END DO ! jj END DO ! jl ! Normalization factor : aksum, ensures mass conservation DO jj = 1, jpj DO ji = 1, jpi aksum(ji,jj) = athorn(ji,jj,0) END DO END DO DO jl = 1, jpl DO jj = 1, jpj DO ji = 1, jpi aksum(ji,jj) = aksum(ji,jj) & + aridge(ji,jj,jl) * (1.0 - 1.0/krdg(ji,jj,jl)) & + araft (ji,jj,jl) * (1.0 - 1.0/kraft) END DO END DO END DO END SUBROUTINE lim_itd_me_ridgeprep !=============================================================================== SUBROUTINE lim_itd_me_ridgeshift(opning, closing_gross, & msnow_mlt, esnow_mlt) ! (subroutine 4/6) !!----------------------------------------------------------------------------- !! *** ROUTINE lim_itd_me_icestrength *** !! ** Purpose : !! This routine 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. !! !! ** Arguments : !! !! ** Inputs / Ouputs : !! !! ** External : !! REAL (wp), DIMENSION(jpi,jpj), INTENT(IN) :: & opning, & ! rate of opening due to divergence/shear closing_gross ! rate at which area removed, not counting ! area of new ridges REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: & msnow_mlt, & ! mass of snow added to ocean (kg m-2) esnow_mlt ! energy needed to melt snow in ocean (J m-2) INTEGER :: & ji, jj, & ! horizontal indices jl, jl1, jl2, & ! thickness category indices jk, & ! ice layer index ij, & ! horizontal index, combines i and j loops icells ! number of cells with aicen > puny INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: & indxi, indxj ! compressed indices REAL(wp), DIMENSION(jpi,jpj) :: & vice_init, vice_final, & ! ice volume summed over categories eice_init, eice_final ! ice energy summed over layers REAL(wp), DIMENSION(jpi,jpj,jpl) :: & aicen_init, & ! ice area before ridging vicen_init, & ! ice volume before ridging vsnon_init, & ! snow volume before ridging esnon_init, & ! snow energy before ridging smv_i_init, & ! ice salinity before ridging oa_i_init ! ice age before ridging REAL(wp), DIMENSION(jpi,jpj,jkmax,jpl) :: & eicen_init ! ice energy before ridging REAL(wp), DIMENSION(jpi,jpj) :: & afrac , & ! fraction of category area ridged ardg1 , & ! area of ice ridged ardg2 , & ! area of new ridges vsrdg , & ! snow volume of ridging ice esrdg , & ! snow energy of ridging ice oirdg1 , & ! areal age content of ridged ice oirdg2 , & ! areal age content of ridging ice dhr , & ! hrmax - hrmin dhr2 , & ! hrmax^2 - hrmin^2 fvol ! fraction of new ridge volume going to n2 REAL(wp), DIMENSION(jpi,jpj) :: & vrdg1 , & ! volume of ice ridged vrdg2 , & ! volume of new ridges vsw , & ! volume of seawater trapped into ridges srdg1 , & ! sal*volume of ice ridged srdg2 , & ! sal*volume of new ridges smsw ! sal*volume of water trapped into ridges REAL(wp), DIMENSION(jpi,jpj) :: & afrft , & ! fraction of category area rafted arft1 , & ! area of ice rafted arft2 , & ! area of new rafted zone virft , & ! ice volume of rafting ice vsrft , & ! snow volume of rafting ice esrft , & ! snow energy of rafting ice smrft , & ! salinity of rafting ice oirft1 , & ! areal age content of rafted ice oirft2 ! areal age content of rafting ice REAL(wp), DIMENSION(jpi,jpj,jkmax) :: & eirft , & ! ice energy of rafting ice erdg1 , & ! enth*volume of ice ridged erdg2 , & ! enth*volume of new ridges ersw ! enth of water trapped into ridges REAL(wp) :: & hL, hR , & ! left and right limits of integration farea , & ! fraction of new ridge area going to n2 zdummy , & ! dummy argument zdummy0 , & ! dummy argument ztmelts , & ! ice melting point sm_newridge ! new ridged ice salinity CHARACTER (len=80) :: & fieldid ! field identifier LOGICAL, PARAMETER :: & l_conservation_check = .true. ! if true, check conservation ! (useful for debugging) LOGICAL :: & neg_ato_i , & ! flag for ato_i(i,j) < -puny large_afrac , & ! flag for afrac > 1 large_afrft ! flag for afrac > 1 REAL(wp) :: & zeps , & epsi10 , & zindb ! switch for the presence of ridge poros or not !---------------------------------------------------------------------------- ! Conservation check eice_init(:,:) = 0.0 IF ( con_i ) THEN CALL lim_column_sum (jpl, v_i, vice_init ) WRITE(numout,*) ' vice_init : ', vice_init(jiindx,jjindx) CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_init ) WRITE(numout,*) ' eice_init : ', eice_init(jiindx,jjindx) ENDIF zeps = 1.0d-20 epsi10 = 1.0d-10 !------------------------------------------------------------------------------- ! 1) Compute change in open water area due to closing and opening. !------------------------------------------------------------------------------- neg_ato_i = .false. DO jj = 1, jpj DO ji = 1, jpi ato_i(ji,jj) = ato_i(ji,jj) & - athorn(ji,jj,0)*closing_gross(ji,jj)*rdt_ice & + opning(ji,jj)*rdt_ice IF (ato_i(ji,jj) .LT. -epsi11) THEN neg_ato_i = .true. ELSEIF (ato_i(ji,jj) .LT. 0.0) THEN ! roundoff error ato_i(ji,jj) = 0.0 ENDIF END DO !jj END DO !ji ! if negative open water area alert it IF (neg_ato_i) THEN ! there is a bug DO jj = 1, jpj DO ji = 1, jpi IF (ato_i(ji,jj) .LT. -epsi11) THEN WRITE(numout,*) '' WRITE(numout,*) 'Ridging error: ato_i < 0' WRITE(numout,*) 'ato_i : ', ato_i(ji,jj) ENDIF ! ato_i < -epsi11 END DO ! ji END DO ! jj ENDIF ! neg_ato_i !----------------------------------------------------------------- ! 2) Save initial state variables !----------------------------------------------------------------- DO jl = 1, jpl DO jj = 1, jpj DO ji = 1, jpi aicen_init(ji,jj,jl) = a_i(ji,jj,jl) vicen_init(ji,jj,jl) = v_i(ji,jj,jl) vsnon_init(ji,jj,jl) = v_s(ji,jj,jl) smv_i_init(ji,jj,jl) = smv_i(ji,jj,jl) oa_i_init (ji,jj,jl) = oa_i(ji,jj,jl) END DO !ji END DO ! jj END DO !jl esnon_init(:,:,:) = e_s(:,:,1,:) DO jl = 1, jpl DO jk = 1, nlay_i DO jj = 1, jpj DO ji = 1, jpi eicen_init(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) END DO !ji END DO !jj END DO !jk END DO !jl ! !----------------------------------------------------------------- ! 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 (aicen_init(ji,jj,jl1) .GT. epsi11 .AND. athorn(ji,jj,jl1) .GT. 0.0 & .AND. closing_gross(ji,jj) > 0.0) THEN icells = icells + 1 indxi(icells) = ji indxj(icells) = jj ENDIF ! test on a_icen_init END DO ! ji END DO ! jj large_afrac = .false. large_afrft = .false. !CDIR NODEP DO ij = 1, icells ji = indxi(ij) jj = indxj(ij) !-------------------------------------------------------------------- ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2) !-------------------------------------------------------------------- ardg1(ji,jj) = aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice arft1(ji,jj) = araft (ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice ardg2(ji,jj) = ardg1(ji,jj) / krdg(ji,jj,jl1) arft2(ji,jj) = arft1(ji,jj) / kraft oirdg1(ji,jj)= aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice oirft1(ji,jj)= araft (ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1) oirft2(ji,jj)= oirft1(ji,jj) / kraft !--------------------------------------------------------------- ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1 !--------------------------------------------------------------- afrac(ji,jj) = ardg1(ji,jj) / aicen_init(ji,jj,jl1) !ridging afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting IF (afrac(ji,jj) > 1.0 + epsi11) THEN !riging large_afrac = .true. ELSEIF (afrac(ji,jj) > 1.0) THEN ! roundoff error afrac(ji,jj) = 1.0 ENDIF IF (afrft(ji,jj) > 1.0 + epsi11) THEN !rafting large_afrft = .true. ELSEIF (afrft(ji,jj) > 1.0) THEN ! roundoff error afrft(ji,jj) = 1.0 ENDIF !-------------------------------------------------------------------------- ! 3.4) Subtract area, volume, and energy from ridging ! / rafting category n1. !-------------------------------------------------------------------------- vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) / & ( 1.0 + ridge_por ) vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + ridge_por ) vsw (ji,jj) = vrdg1(ji,jj) * ridge_por vsrdg(ji,jj) = vsnon_init(ji,jj,jl1) * afrac(ji,jj) esrdg(ji,jj) = esnon_init(ji,jj,jl1) * afrac(ji,jj) srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) / & ( 1. + ridge_por ) srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) ! rafting volumes, heat contents ... virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj) vsrft(ji,jj) = vsnon_init(ji,jj,jl1) * afrft(ji,jj) esrft(ji,jj) = esnon_init(ji,jj,jl1) * afrft(ji,jj) smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj) ! substract everything a_i(ji,jj,jl1) = a_i(ji,jj,jl1) - ardg1(ji,jj) - arft1(ji,jj) v_i(ji,jj,jl1) = v_i(ji,jj,jl1) - vrdg1(ji,jj) - virft(ji,jj) v_s(ji,jj,jl1) = v_s(ji,jj,jl1) - vsrdg(ji,jj) - vsrft(ji,jj) e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg(ji,jj) - esrft(ji,jj) oa_i(ji,jj,jl1) = oa_i(ji,jj,jl1) - oirdg1(ji,jj) - oirft1(ji,jj) smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1(ji,jj) - smrft(ji,jj) !----------------------------------------------------------------- ! 3.5) Compute properties of new ridges !----------------------------------------------------------------- !------------- ! Salinity !------------- smsw(ji,jj) = sss_m(ji,jj) * vsw(ji,jj) * ridge_por ! salinity of new ridge sm_newridge = ( srdg1(ji,jj) + smsw(ji,jj) ) / vrdg2(ji,jj) zdummy = sm_newridge * vrdg2(ji,jj) ! has to be smaller than s_i_max sm_newridge = MIN( s_i_max, sm_newridge ) ! salt flux due to ridge creation fsalt_rpo(ji,jj) = fsalt_rpo(ji,jj) + & MAX ( zdummy - srdg2(ji,jj) , 0.0 ) & * rhoic / rdt_ice ! sal times volume for new ridges srdg2(ji,jj) = sm_newridge * vrdg2(ji,jj) !------------------------------------ ! 3.6 Increment ridging diagnostics !------------------------------------ ! jl1 looping 1-jpl ! ij looping 1-icells dardg1dt(ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj) dardg2dt(ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj) diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) / rdt_ice opening(ji,jj) = opening (ji,jj) + opning(ji,jj)*rdt_ice IF (con_i) vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj) !------------------------------------------ ! 3.7 Put the snow somewhere in the ocean !------------------------------------------ ! Place part of the snow lost by ridging into the ocean. ! Note that esnow_mlt < 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. ! jl1 looping 1-jpl ! ij looping 1-icells msnow_mlt(ji,jj) = msnow_mlt(ji,jj) & + rhosn*vsrdg(ji,jj)*(1.0-fsnowrdg) & !rafting included + rhosn*vsrft(ji,jj)*(1.0-fsnowrft) esnow_mlt(ji,jj) = esnow_mlt(ji,jj) & + esrdg(ji,jj)*(1.0-fsnowrdg) & !rafting included + esrft(ji,jj)*(1.0-fsnowrft) !----------------------------------------------------------------- ! 3.8 Compute quantities used to apportion ice among categories ! in the n2 loop below !----------------------------------------------------------------- ! jl1 looping 1-jpl ! ij looping 1-icells dhr(ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) & - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) END DO ! ij !-------------------------------------------------------------------- ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and ! compute ridged ice enthalpy !-------------------------------------------------------------------- DO jk = 1, nlay_i !CDIR NODEP DO ij = 1, icells ji = indxi(ij) jj = indxj(ij) ! heat content of ridged ice erdg1(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) / & ( 1.0 + ridge_por ) eirft(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj) e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) & - erdg1(ji,jj,jk) & - eirft(ji,jj,jk) ! sea water heat content ztmelts = - tmut * sss_m(ji,jj) + rtt ! heat content per unit volume zdummy0 = - rcp * ( sst_m(ji,jj) + rt0 - rtt ) * vsw(ji,jj) ! corrected sea water salinity zindb = MAX( 0.0, SIGN( 1.0, vsw(ji,jj) - zeps ) ) zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / & MAX( ridge_por * vsw(ji,jj), zeps ) ztmelts = - tmut * zdummy + rtt ersw(ji,jj,jk) = - rcp * ( ztmelts - rtt ) * vsw(ji,jj) ! heat flux fheat_rpo(ji,jj) = fheat_rpo(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) / & rdt_ice ! Correct dimensions to avoid big values ersw(ji,jj,jk) = ersw(ji,jj,jk) / 1.0d+09 ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J ersw(ji,jj,jk) = ersw(ji,jj,jk) * & area(ji,jj) * vsw(ji,jj) / & nlay_i erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk) END DO ! ij END DO !jk IF ( con_i ) THEN DO jk = 1, nlay_i !CDIR NODEP DO ij = 1, icells ji = indxi(ij) jj = indxj(ij) eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - & erdg1(ji,jj,jk) END DO ! ij END DO !jk ENDIF IF (large_afrac) THEN ! there is a bug !CDIR NODEP DO ij = 1, icells ji = indxi(ij) jj = indxj(ij) IF ( afrac(ji,jj) > 1.0 + epsi11 ) THEN WRITE(numout,*) '' WRITE(numout,*) ' ardg > a_i' WRITE(numout,*) ' ardg, aicen_init : ', & ardg1(ji,jj), aicen_init(ji,jj,jl1) ENDIF ! afrac > 1 + puny ENDDO ! if ENDIF ! large_afrac IF (large_afrft) THEN ! there is a bug !CDIR NODEP DO ij = 1, icells ji = indxi(ij) jj = indxj(ij) IF ( afrft(ji,jj) > 1.0 + epsi11 ) THEN WRITE(numout,*) '' WRITE(numout,*) ' arft > a_i' WRITE(numout,*) ' arft, aicen_init : ', & arft1(ji,jj), aicen_init(ji,jj,jl1) ENDIF ! afrft > 1 + puny ENDDO ! if ENDIF ! large_afrft !------------------------------------------------------------------------------- ! 4) Add area, volume, and energy of new ridge to each category jl2 !------------------------------------------------------------------------------- ! jl1 looping 1-jpl DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! over categories to which ridged ice is transferred !CDIR NODEP DO ij = 1, icells ji = indxi(ij) jj = indxj(ij) ! Compute the fraction of ridged ice area and volume going to ! thickness category jl2. ! Transfer area, volume, and energy accordingly. IF (hrmin(ji,jj,jl1) .GE. hi_max(jl2) .OR. & hrmax(ji,jj,jl1) .LE. hi_max(jl2-1)) THEN hL = 0.0 hR = 0.0 ELSE hL = MAX (hrmin(ji,jj,jl1), hi_max(jl2-1)) hR = MIN (hrmax(ji,jj,jl1), hi_max(jl2)) ENDIF ! fraction of ridged ice area and volume going to n2 farea = (hR-hL) / dhr(ji,jj) fvol(ji,jj) = (hR*hR - hL*hL) / dhr2(ji,jj) a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + farea * ardg2(ji,jj) v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + fvol(ji,jj) * vrdg2(ji,jj) v_s(ji,jj,jl2) = v_s(ji,jj,jl2) & + fvol(ji,jj) * vsrdg(ji,jj) * fsnowrdg e_s(ji,jj,1,jl2) = e_s(ji,jj,1,jl2) & + fvol(ji,jj) * esrdg(ji,jj) * fsnowrdg smv_i(ji,jj,jl2) = smv_i(ji,jj,jl2) + fvol(ji,jj) * srdg2(ji,jj) oa_i(ji,jj,jl2) = oa_i(ji,jj,jl2) + farea * oirdg2(ji,jj) END DO ! ij ! Transfer ice energy to category jl2 by ridging DO jk = 1, nlay_i !CDIR NODEP DO ij = 1, icells ji = indxi(ij) jj = indxj(ij) e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) & + fvol(ji,jj)*erdg2(ji,jj,jk) END DO ! ij END DO !jk END DO ! jl2 (new ridges) DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2) !CDIR NODEP DO ij = 1, icells ji = indxi(ij) jj = indxj(ij) ! Compute the fraction of rafted ice area and volume going to ! thickness category jl2, transfer area, volume, and energy accordingly. IF (hraft(ji,jj,jl1) .LE. hi_max(jl2) .AND. & hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + arft2(ji,jj) v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + virft(ji,jj) v_s(ji,jj,jl2) = v_s(ji,jj,jl2) & + vsrft(ji,jj)*fsnowrft e_s(ji,jj,1,jl2) = e_s(ji,jj,1,jl2) & + esrft(ji,jj)*fsnowrft smv_i(ji,jj,jl2) = smv_i(ji,jj,jl2) & + smrft(ji,jj) oa_i(ji,jj,jl2) = oa_i(ji,jj,jl2) & + oirft2(ji,jj) ENDIF ! hraft END DO ! ij ! Transfer rafted ice energy to category jl2 DO jk = 1, nlay_i !CDIR NODEP DO ij = 1, icells ji = indxi(ij) jj = indxj(ij) IF (hraft(ji,jj,jl1) .LE. hi_max(jl2) .AND. & hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) & + eirft(ji,jj,jk) ENDIF END DO ! ij END DO !jk END DO ! jl2 END DO ! jl1 (deforming categories) ! Conservation check IF ( con_i ) THEN CALL lim_column_sum (jpl, v_i, vice_final) fieldid = ' v_i : limitd_me ' CALL lim_cons_check (vice_init, vice_final, 1.0e-6, fieldid) WRITE(numout,*) ' vice_init : ', vice_init(jiindx,jjindx) WRITE(numout,*) ' vice_final : ', vice_final(jiindx,jjindx) CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_final ) fieldid = ' e_i : limitd_me ' CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid) WRITE(numout,*) ' eice_init : ', eice_init(jiindx,jjindx) WRITE(numout,*) ' eice_final : ', eice_final(jiindx,jjindx) ENDIF END SUBROUTINE lim_itd_me_ridgeshift !============================================================================== SUBROUTINE lim_itd_me_asumr !(subroutine 5/6) !!----------------------------------------------------------------------------- !! *** ROUTINE lim_itd_me_asumr *** !! ** Purpose : !! This routine finds total fractional area !! !! ** Method : !! Find the total area of ice plus open water in each grid cell. !! !! This is similar to the aggregate_area subroutine except that the !! total area can be greater than 1, so the open water area is !! included in the sum instead of being computed as a residual. !! !! ** Arguments : INTEGER :: ji, jj, jl !----------------------------------------------------------------- ! open water !----------------------------------------------------------------- DO jj = 1, jpj DO ji = 1, jpi asum(ji,jj) = ato_i(ji,jj) END DO END DO !----------------------------------------------------------------- ! ice categories !----------------------------------------------------------------- DO jl = 1, jpl DO jj= 1, jpj DO ji = 1, jpi asum(ji,jj) = asum(ji,jj) + a_i(ji,jj,jl) END DO !ji END DO !jj END DO ! jl END SUBROUTINE lim_itd_me_asumr !============================================================================== SUBROUTINE lim_itd_me_init ! (subroutine 6/6) !!------------------------------------------------------------------- !! *** ROUTINE lim_itd_me_init *** !! !! ** Purpose : Physical constants and parameters linked !! to the mechanical ice redistribution !! !! ** Method : Read the namiceitdme namelist !! and check the parameters values !! called at the first timestep (nit000) !! !! ** input : Namelist namiceitdme !! !! history : !! 9.0, LIM3.0 - 02-2006 (M. Vancoppenolle) original code !!------------------------------------------------------------------- NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft,& Gstar, astar, & Hstar, raftswi, hparmeter, Craft, ridge_por, & sal_max_ridge, partfun_swi, transfun_swi, & brinstren_swi !!------------------------------------------------------------------- ! Define the initial parameters ! ------------------------- REWIND( numnam_ice ) READ ( numnam_ice , namiceitdme) IF (lwp) THEN WRITE(numout,*) WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution ' WRITE(numout,*)' ~~~~~~~~~~~~~~~' WRITE(numout,*)' Switch choosing the ice redistribution scheme ridge_scheme_swi', ridge_scheme_swi WRITE(numout,*)' Fraction of shear energy contributing to ridging Cs ', Cs WRITE(numout,*)' Ratio of ridging work to PotEner change in ridging Cf ', Cf WRITE(numout,*)' Fraction of snow volume conserved during ridging fsnowrdg ', fsnowrdg WRITE(numout,*)' Fraction of snow volume conserved during ridging fsnowrft ', fsnowrft WRITE(numout,*)' Fraction of total ice coverage contributing to ridging Gstar ', Gstar WRITE(numout,*)' Equivalent to G* for an exponential part function astar ', astar WRITE(numout,*)' Quantity playing a role in max ridged ice thickness Hstar ', Hstar WRITE(numout,*)' Rafting of ice sheets or not raftswi ', raftswi WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) hparmeter ', hparmeter WRITE(numout,*)' Rafting hyperbolic tangent coefficient Craft ', Craft WRITE(numout,*)' Initial porosity of ridges ridge_por ', ridge_por WRITE(numout,*)' Maximum salinity of ridging ice sal_max_ridge ', sal_max_ridge WRITE(numout,*)' Switch for part. function (0) linear (1) exponential partfun_swi ', partfun_swi WRITE(numout,*)' Switch for tran. function (0) linear (1) exponential transfun_swi ', transfun_swi WRITE(numout,*)' Switch for including brine volume in ice strength comp. brinstren_swi ', brinstren_swi ENDIF END SUBROUTINE lim_itd_me_init !============================================================================== SUBROUTINE lim_itd_me_zapsmall !!------------------------------------------------------------------- !! *** ROUTINE lim_itd_me_zapsmall *** !! !! ** Purpose : Remove too small sea ice areas and correct salt fluxes !! !! !! history : !! author: William H. Lipscomb, LANL !! Nov 2003: Modified by Julie Schramm to conserve volume and energy !! Sept 2004: Modified by William Lipscomb; replaced normalize_state with !! additions to local freshwater, salt, and heat fluxes !! 9.0, LIM3.0 - 02-2006 (M. Vancoppenolle) original code !!------------------------------------------------------------------- INTEGER :: & ji,jj, & ! horizontal indices jl, & ! ice category index jk, & ! ice layer index ! ij, & ! combined i/j horizontal index icells ! number of cells with ice to zap ! INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: & ! indxi, & ! compressed indices for i/j directions ! indxj INTEGER, DIMENSION(jpi,jpj) :: zmask REAL(wp) :: & xtmp ! temporary variable DO jl = 1, jpl !----------------------------------------------------------------- ! Count categories to be zapped. ! Abort model in case of negative area. !----------------------------------------------------------------- IF( MINVAL(a_i(:,:,jl)) .LT. -epsi11 ) THEN DO jj = 1, jpj DO ji = 1, jpi IF ( a_i(ji,jj,jl) .LT. -epsi11 ) THEN WRITE (numout,*) ' ALERTE 98 ' WRITE (numout,*) ' Negative ice area: ji, jj, jl: ', ji, jj,jl WRITE (numout,*) ' a_i ', a_i(ji,jj,jl) ENDIF END DO END DO ENDIF icells = 0 zmask = 0.e0 DO jj = 1, jpj DO ji = 1, jpi IF ( ( a_i(ji,jj,jl) .GE. -epsi11 .AND. a_i(ji,jj,jl) .LT. 0.0) & .OR. & ( a_i(ji,jj,jl) .GT. 0.0 .AND. a_i(ji,jj,jl) .LE. 1.0e-11 ) & .OR. & !new line ( v_i(ji,jj,jl) .EQ. 0.0 .AND. a_i(ji,jj,jl) .GT. 0.0 ) & .OR. & ( v_i(ji,jj,jl) .GT. 0.0 .AND. v_i(ji,jj,jl) .LT. 1.e-12 ) ) THEN zmask(ji,jj) = 1 ENDIF END DO END DO WRITE(numout,*) SUM(zmask), ' cells of ice zapped in the ocean ' !----------------------------------------------------------------- ! Zap ice energy and use ocean heat to melt ice !----------------------------------------------------------------- DO jk = 1, nlay_i DO jj = 1 , jpj DO ji = 1 , jpi xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) / rdt_ice xtmp = xtmp * unit_fac ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1 - zmask(ji,jj) ) END DO ! ji END DO ! jj END DO ! jk DO jj = 1 , jpj DO ji = 1 , jpi !----------------------------------------------------------------- ! Zap snow energy and use ocean heat to melt snow !----------------------------------------------------------------- ! xtmp = esnon(i,j,n) / dt ! < 0 ! fhnet(i,j) = fhnet(i,j) + xtmp ! fhnet_hist(i,j) = fhnet_hist(i,j) + xtmp ! xtmp is greater than 0 ! fluxes are positive to the ocean ! here the flux has to be negative for the ocean xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice !RB ??????? t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1 - zmask(ji,jj) ) !----------------------------------------------------------------- ! zap ice and snow volume, add water and salt to ocean !----------------------------------------------------------------- ! xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) / dt ! fresh(i,j) = fresh(i,j) + xtmp ! fresh_hist(i,j) = fresh_hist(i,j) + xtmp ! fsalt_res(ji,jj) = fsalt_res(ji,jj) + ( sss_m(ji,jj) ) * & ! rhosn * v_s(ji,jj,jl) / rdt_ice ! fsalt_res(ji,jj) = fsalt_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) ) * & ! rhoic * v_i(ji,jj,jl) / rdt_ice ! emps(i,j) = emps(i,j) + xtmp ! fsalt_hist(i,j) = fsalt_hist(i,j) + xtmp ato_i(ji,jj) = a_i(ji,jj,jl) * zmask(ji,jj) + ato_i(ji,jj) a_i(ji,jj,jl) = a_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) v_i(ji,jj,jl) = v_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) v_s(ji,jj,jl) = v_s(ji,jj,jl) * ( 1 - zmask(ji,jj) ) t_su(ji,jj,jl) = t_su(ji,jj,jl) * (1 -zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj) oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) END DO ! ji END DO ! jj END DO ! jl END SUBROUTINE lim_itd_me_zapsmall #else !!====================================================================== !! *** MODULE limitd_me *** !! no sea ice model !!====================================================================== CONTAINS SUBROUTINE lim_itd_me ! Empty routines END SUBROUTINE lim_itd_me SUBROUTINE lim_itd_me_icestrength END SUBROUTINE lim_itd_me_icestrength SUBROUTINE lim_itd_me_ridgeprep END SUBROUTINE lim_itd_me_ridgeprep SUBROUTINE lim_itd_th_me_ridgeshift END SUBROUTINE lim_itd_me_ridgeshift SUBROUTINE lim_itd_me_asumr END SUBROUTINE lim_itd_me_asumr SUBROUTINE lim_itd_me_sort END SUBROUTINE lim_itd_me_sort SUBROUTINE lim_itd_me_init END SUBROUTINE lim_itd_me_init SUBROUTINE lim_itd_me_zapsmall END SUBROUTINE lim_itd_me_zapsmall #endif END MODULE limitd_me