MODULE limrhg !!====================================================================== !! *** MODULE limrhg *** !! Ice rheology : sea ice rheology !!====================================================================== #if defined key_lim3 !!---------------------------------------------------------------------- !! 'key_lim3' LIM3 sea-ice model !!---------------------------------------------------------------------- !! lim_rhg : computes ice velocities !!---------------------------------------------------------------------- !! * Modules used USE phycst USE par_oce USE ice_oce ! ice variables USE dom_oce USE dom_ice USE ice USE iceini USE lbclnk USE lib_mpp USE in_out_manager ! I/O manager USE limitd_me USE prtctl ! Print control IMPLICIT NONE PRIVATE !! * Routine accessibility PUBLIC lim_rhg ! routine called by lim_dyn !! * Module variables REAL(wp) :: & ! constant values rzero = 0.e0 , & rone = 1.e0 !!---------------------------------------------------------------------- !! LIM 3.0, UCL-LOCEAN-IPSL (2008) !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limrhg.F90,v 1.5 2005/03/27 18:34:42 opalod Exp $ !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt !!---------------------------------------------------------------------- CONTAINS SUBROUTINE lim_rhg( k_j1, k_jpj ) !!------------------------------------------------------------------- !! *** SUBROUTINE lim_rhg *** !! EVP-C-grid !! !! ** purpose : determines sea ice drift from wind stress, ice-ocean !! stress and sea-surface slope. Ice-ice interaction is described by !! a non-linear elasto-viscous-plastic (EVP) law including shear !! strength and a bulk rheology (Hunke and Dukowicz, 2002). !! !! The points in the C-grid look like this, dear reader !! !! (ji,jj) !! | !! | !! (ji-1,jj) | (ji,jj) !! --------- !! | | !! | (ji,jj) |------(ji,jj) !! | | !! --------- !! (ji-1,jj-1) (ji,jj-1) !! !! ** Inputs : - wind forcing (stress), oceanic currents !! ice total volume (vt_i) per unit area !! snow total volume (vt_s) per unit area !! !! ** Action : - compute u_ice, v_ice : the components of the !! sea-ice velocity vector !! - compute delta_i, shear_i, divu_i, which are inputs !! of the ice thickness distribution !! !! ** Steps : 1) Compute ice snow mass, ice strength !! 2) Compute wind, oceanic stresses, mass terms and !! coriolis terms of the momentum equation !! 3) Solve the momentum equation (iterative procedure) !! 4) Prevent high velocities if the ice is thin !! 5) Recompute invariants of the strain rate tensor !! which are inputs of the ITD, store stress !! for the next time step !! 6) Control prints of residual (convergence) !! and charge ellipse. !! The user should make sure that the parameters !! nevp, telast and creepl maintain stress state !! on the charge ellipse for plastic flow !! e.g. in the Canadian Archipelago !! !! ** References : Hunke and Dukowicz, JPO97 !! Bouillon et al., 08, in prep (update this when !! published) !! Vancoppenolle et al., OM08 !! !! History : !! 1.0 ! 07-03 (M.A. Morales Maqueda, S. Bouillon) !! 2.0 ! 08-03 M. Vancoppenolle : LIM3 !! !!------------------------------------------------------------------- ! * Arguments ! INTEGER, INTENT(in) :: & k_j1 , & !: southern j-index for ice computation k_jpj !: northern j-index for ice computation ! * Local variables INTEGER :: ji, jj !: dummy loop indices INTEGER :: & iter, jter !: temporary integers CHARACTER (len=50) :: charout REAL(wp) :: & zt11, zt12, zt21, zt22, & !: temporary scalars ztagnx, ztagny, & !: wind stress on U/V points delta ! REAL(wp) :: & za, & !: zstms, & !: temporary scalar for ice strength zsang, & !: temporary scalar for coriolis term zmask !: mask for the computation of ice mass REAL(wp),DIMENSION(jpi,jpj) :: & zpresh , & !: temporary array for ice strength zpreshc , & !: Ice strength on grid cell corners (zpreshc) zfrld1, zfrld2, & !: lead fraction on U/V points zmass1, zmass2, & !: ice/snow mass on U/V points zcorl1, zcorl2, & !: coriolis parameter on U/V points za1ct, za2ct , & !: temporary arrays zc1 , & !: ice mass zusw , & !: temporary weight for the computation !: of ice strength u_oce1, v_oce1, & !: ocean u/v component on U points u_oce2, v_oce2, & !: ocean u/v component on V points u_ice2, & !: ice u component on V point v_ice1 !: ice v component on U point REAL(wp) :: & dtevp, & !: time step for subcycling dtotel, & !: ecc2, & !: square of yield ellipse eccenticity z0, & !: temporary scalar zr, & !: temporary scalar zcca, zccb, & !: temporary scalars zu_ice2, & !: zv_ice1, & !: zddc, zdtc, & !: temporary array for delta on corners zdst, & !: temporary array for delta on centre zdsshx, zdsshy, & !: term for the gradient of ocean surface sigma1, sigma2 !: internal ice stress REAL(wp),DIMENSION(jpi,jpj) :: & zf1, zf2 !: arrays for internal stresses REAL(wp),DIMENSION(jpi,jpj) :: & zdd, zdt, & !: Divergence and tension at centre of grid cells zds, & !: Shear on northeast corner of grid cells deltat, & !: Delta at centre of grid cells deltac, & !: Delta on corners zs1, zs2, & !: Diagonal stress tensor components zs1 and zs2 zs12 !: Non-diagonal stress tensor component zs12 REAL(wp) :: & zresm , & !: Maximal error on ice velocity zindb , & !: ice (1) or not (0) zdummy !: dummy argument REAL(wp),DIMENSION(jpi,jpj) :: & zu_ice , & !: Ice velocity on previous time step zv_ice , & zresr !: Local error on velocity INTEGER :: & stress_mean_swi = 1 ! !------------------------------------------------------------------------------! ! 1) Ice-Snow mass (zc1), ice strength (zpresh) ! !------------------------------------------------------------------------------! ! ! Put every vector to 0 zpresh(:,:) = 0.0 ; zpreshc(:,:) = 0.0 ; zfrld1(:,:) = 0.0 zmass1(:,:) = 0.0 ; zcorl1(:,:) = 0.0 ; zcorl2(:,:) = 0.0 za1ct(:,:) = 0.0 ; za2ct(:,:) = 0.0 zc1(:,:) = 0.0 ; zusw(:,:) = 0.0 u_oce1(:,:) = 0.0 ; v_oce1(:,:) = 0.0 ; u_oce2(:,:) = 0.0 v_oce2(:,:) = 0.0 ; u_ice2(:,:) = 0.0 ; v_ice1(:,:) = 0.0 zf1(:,:) = 0.0 ; zf2(:,:) = 0.0 zdd(:,:) = 0.0 ; zdt(:,:) = 0.0 ; zds(:,:) = 0.0 deltat(:,:) = 0.0 ; deltac(:,:) = 0.0 zpresh(:,:) = 0.0 ! Ice strength on T-points CALL lim_itd_me_icestrength(ridge_scheme_swi) ! Ice mass and temp variables DO jj = k_j1 , k_jpj-1 DO ji = 1 , jpi zc1(ji,jj) = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) ) zpresh(ji,jj) = tms(ji,jj) * strength(ji,jj) / 2. ! tmi = 1 where ther is ice or on land ! Claude sees a bug, next line : tmi(ji,jj) tmi = 1.0 - ( 1.0 - MAX( 0.0 , SIGN ( 1.0 , vt_i(ji,jj) - & epsd ) ) ) * tms(ji,jj) END DO END DO CALL lbc_lnk( zc1(:,:), 'T', 1. ) CALL lbc_lnk( zpresh(:,:), 'T', 1. ) ! Claude sees a bug ! CALL lbc_lnk( tmi(:,:), 'T', 1. ) ! Ice strength on grid cell corners (zpreshc) ! needed for calculation of shear stress DO jj = k_j1+1, k_jpj-1 DO ji = 2, jpim1 zstms = tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & & tms(ji,jj+1) * wght(ji+1,jj+1,1,2) + & & tms(ji+1,jj) * wght(ji+1,jj+1,2,1) + & & tms(ji,jj) * wght(ji+1,jj+1,1,1) zusw(ji,jj) = 1.0 / MAX( zstms, epsd ) zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & & zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + & & zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + & & zpresh(ji,jj) * wght(ji+1,jj+1,1,1) & & ) * zusw(ji,jj) END DO END DO CALL lbc_lnk( zpreshc(:,:), 'F', 1. ) ! !------------------------------------------------------------------------------! ! 2) Wind / ocean stress, mass terms, coriolis terms !------------------------------------------------------------------------------! ! ! Wind stress, coriolis and mass terms on the sides of the squares ! zfrld1: lead fraction on U-points ! zfrld2: lead fraction on V-points ! zmass1: ice/snow mass on U-points ! zmass2: ice/snow mass on V-points ! zcorl1: Coriolis parameter on U-points ! zcorl2: Coriolis parameter on V-points ! (ztagnx,ztagny): wind stress on U/V points ! u_oce1: ocean u component on u points ! v_oce1: ocean v component on u points ! u_oce2: ocean u component on v points ! v_oce2: ocean v component on v points DO jj = k_j1+1, k_jpj-1 DO ji = 2, jpim1 zt11 = tms(ji,jj)*e1t(ji,jj) zt12 = tms(ji+1,jj)*e1t(ji+1,jj) zt21 = tms(ji,jj)*e2t(ji,jj) zt22 = tms(ji,jj+1)*e2t(ji,jj+1) ! Leads area. zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + & & zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd ) zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + & & zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd ) ! Mass, coriolis coeff. and currents zmass1(ji,jj) = ( zt12*zc1(ji,jj) + zt11*zc1(ji+1,jj) ) / (zt11+zt12+epsd) zmass2(ji,jj) = ( zt22*zc1(ji,jj) + zt21*zc1(ji,jj+1) ) / (zt21+zt22+epsd) zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + & e1t(ji,jj)*fcor(ji+1,jj) ) & / (e1t(ji,jj) + e1t(ji+1,jj) + epsd ) zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)*fcor(ji,jj) + & e2t(ji,jj)*fcor(ji,jj+1) ) & / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd ) ! u_oce1(ji,jj) = u_io(ji,jj) v_oce2(ji,jj) = v_io(ji,jj) ! Ocean has no slip boundary condition v_oce1(ji,jj) = 0.5*( (v_io(ji,jj)+v_io(ji,jj-1))*e1t(ji+1,jj) & & +(v_io(ji+1,jj)+v_io(ji+1,jj-1))*e1t(ji,jj)) & & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) u_oce2(ji,jj) = 0.5*((u_io(ji,jj)+u_io(ji-1,jj))*e2t(ji,jj+1) & & +(u_io(ji,jj+1)+u_io(ji-1,jj+1))*e2t(ji,jj)) & & / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) ! Wind stress. ztagnx = ( 1. - zfrld1(ji,jj) ) * gtaux(ji,jj) ztagny = ( 1. - zfrld2(ji,jj) ) * gtauy(ji,jj) ! Computation of the velocity field taking into account the ice internal interaction. ! Terms that are independent of the velocity field. ! SB On utilise maintenant le gradient de la pente de l'ocean ! include it later ! zdsshx = (ssh_io(ji+1,jj) - ssh_io(ji,jj))/e1u(ji,jj) ! zdsshy = (ssh_io(ji,jj+1) - ssh_io(ji,jj))/e2v(ji,jj) zdsshx = 0.0 zdsshy = 0.0 za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy END DO END DO ! !------------------------------------------------------------------------------! ! 3) Solution of the momentum equation, iterative procedure !------------------------------------------------------------------------------! ! ! Time step for subcycling dtevp = rdt_ice / nevp dtotel = dtevp / ( 2.0 * telast ) !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter) ecc2 = ecc*ecc !-Initialise stress tensor IF ( stress_mean_swi .NE. 1 ) THEN zs1(:,:) = 0.0 zs2(:,:) = 0.0 zs12(:,:) = 0.0 ELSE zs1(:,:) = stress1_i(:,:) zs2(:,:) = stress2_i(:,:) zs12(:,:) = stress12_i(:,:) ENDIF v_ice1(:,:) = 0.0 u_ice2(:,:) = 0.0 zdd(:,:) = 0.0 zdt(:,:) = 0.0 zds(:,:) = 0.0 deltat(:,:) = 0.0 !----------------------! DO iter = 1 , nevp ! loop over iter ! !----------------------! DO jj = k_j1, k_jpj-1 zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step zv_ice(:,jj) = v_ice(:,jj) END DO DO jj = k_j1+1, k_jpj-1 DO ji = 2, jpim1 ! !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002) !- zdd(:,:), zdt(:,:): divergence and tension at centre of grid cells !- zds(:,:): shear on northeast corner of grid cells ! !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded, ! there are many repeated calculations. ! Speed could be improved by regrouping terms. For ! the moment, however, the stress is on clarity of coding to avoid ! bugs (Martin, for Miguel). ! !- ALSO: arrays zdd, zdt, zds and delta could ! be removed in the future to minimise memory demand. ! !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of ! grid cells, exactly as in the B grid case. For simplicity, the indexation on ! the corners is the same as in the B grid. ! ! zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & & -e2u(ji-1,jj)*u_ice(ji-1,jj) & & +e1v(ji,jj)*v_ice(ji,jj) & & -e1v(ji,jj-1)*v_ice(ji,jj-1) & & ) & & / area(ji,jj) zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) & & -u_ice(ji-1,jj)/e2u(ji-1,jj) & & )*e2t(ji,jj)*e2t(ji,jj) & & -( v_ice(ji,jj)/e1v(ji,jj) & & -v_ice(ji,jj-1)/e1v(ji,jj-1) & & )*e1t(ji,jj)*e1t(ji,jj) & & ) & & / area(ji,jj) ! zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1) & & -u_ice(ji,jj)/e1u(ji,jj) & & )*e1f(ji,jj)*e1f(ji,jj) & & +( v_ice(ji+1,jj)/e2v(ji+1,jj) & & -v_ice(ji,jj)/e2v(ji,jj) & & )*e2f(ji,jj)*e2f(ji,jj) & & ) & & / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) & & * tmi(ji,jj) * tmi(ji,jj+1) & & * tmi(ji+1,jj) * tmi(ji+1,jj+1) v_ice1(ji,jj) = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) u_ice2(ji,jj) = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1) & & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) END DO END DO CALL lbc_lnk( zdd(:,:), 'T', 1. ) CALL lbc_lnk( zdt(:,:), 'T', 1. ) CALL lbc_lnk( zds(:,:), 'F', 1. ) DO jj = k_j1+1, k_jpj-1 DO ji = 2, jpim1 !- Calculate Delta at centre of grid cells zdst = ( e2u( ji , jj ) * v_ice1(ji,jj) & & - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & & + e1v( ji , jj ) * u_ice2(ji,jj) & & - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) & & ) & & / area(ji,jj) delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + & & ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + & (zdt(ji,jj)**2 + zdst**2)*usecc2), creepl ) !-Calculate stress tensor components zs1 and zs2 !-at centre of grid cells (see section 3.5 of CICE user's guide). zs1(ji,jj) = ( zs1(ji,jj) & & - dtotel*( ( 1.0 - alphaevp) * zs1(ji,jj) + & & ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) & * zpresh(ji,jj) ) ) & & / ( 1.0 + alphaevp * dtotel ) zs2(ji,jj) = ( zs2(ji,jj) & & - dtotel*((1.0-alphaevp)*ecc2*zs2(ji,jj) - & zdt(ji,jj)/deltat(ji,jj)*zpresh(ji,jj)) ) & & / ( 1.0 + alphaevp*ecc2*dtotel ) END DO END DO CALL lbc_lnk( zs1(:,:), 'T', 1. ) CALL lbc_lnk( zs2(:,:), 'T', 1. ) DO jj = k_j1+1, k_jpj-1 DO ji = 2, jpim1 !- Calculate Delta on corners zddc = ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1) & & -v_ice1(ji,jj)/e1u(ji,jj) & & )*e1f(ji,jj)*e1f(ji,jj) & & +( u_ice2(ji+1,jj)/e2v(ji+1,jj) & & -u_ice2(ji,jj)/e2v(ji,jj) & & )*e2f(ji,jj)*e2f(ji,jj) & & ) & & / ( e1f(ji,jj) * e2f(ji,jj) ) zdtc = (-( v_ice1(ji,jj+1)/e1u(ji,jj+1) & & -v_ice1(ji,jj)/e1u(ji,jj) & & )*e1f(ji,jj)*e1f(ji,jj) & & +( u_ice2(ji+1,jj)/e2v(ji+1,jj) & & -u_ice2(ji,jj)/e2v(ji,jj) & & )*e2f(ji,jj)*e2f(ji,jj) & & ) & & / ( e1f(ji,jj) * e2f(ji,jj) ) deltac(ji,jj) = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide). zs12(ji,jj) = ( zs12(ji,jj) & & - dtotel*( (1.0-alphaevp)*ecc2*zs12(ji,jj) - zds(ji,jj) / & & ( 2.0*deltac(ji,jj) ) * zpreshc(ji,jj))) & & / ( 1.0 + alphaevp*ecc2*dtotel ) END DO ! ji END DO ! jj CALL lbc_lnk( zs12(:,:), 'F', 1. ) ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) DO jj = k_j1+1, k_jpj-1 DO ji = 2, jpim1 !- contribution of zs1, zs2 and zs12 to zf1 zf1(ji,jj) = 0.5*( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj) & & +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj) & & +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj) & & ) / ( e1u(ji,jj)*e2u(ji,jj) ) ! contribution of zs1, zs2 and zs12 to zf2 zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) & & -(zs2(ji,jj+1)*e1t(ji,jj+1)**2 - zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) & & + 2.0*(zs12(ji,jj)*e2f(ji,jj)**2 - & zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) & & ) / ( e1v(ji,jj)*e2v(ji,jj) ) END DO END DO ! ! Computation of ice velocity ! ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly. ! IF (MOD(iter,2).eq.0) THEN DO jj = k_j1+1, k_jpj-1 DO ji = 2, jpim1 zmask = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj) zsang = SIGN ( 1.0 , fcor(ji,jj) ) * sangvg z0 = zmass1(ji,jj)/dtevp ! SB modif because ocean has no slip boundary condition zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) za = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj)) zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj)) zcca = z0+za*cangvg zccb = zcorl1(ji,jj)+za*zsang u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask END DO END DO CALL lbc_lnk( u_ice(:,:), 'U', -1. ) DO jj = k_j1+1, k_jpj-1 DO ji = 2, jpim1 zmask = (1.0-max(rzero,sign(rone,-zmass2(ji,jj))))*tmv(ji,jj) zsang = sign(1.0,fcor(ji,jj))*sangvg z0 = zmass2(ji,jj)/dtevp ! SB modif because ocean has no slip boundary condition zu_ice2 = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1) & & + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) za = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj)) zcca = z0+za*cangvg zccb = zcorl2(ji,jj)+za*zsang v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask END DO END DO CALL lbc_lnk( v_ice(:,:), 'V', -1. ) ELSE DO jj = k_j1+1, k_jpj-1 DO ji = 2, jpim1 zmask = (1.0-max(rzero,sign(rone,-zmass2(ji,jj))))*tmv(ji,jj) zsang = sign(1.0,fcor(ji,jj))*sangvg z0 = zmass2(ji,jj)/dtevp ! SB modif because ocean has no slip boundary condition zu_ice2 = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1) & & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) za = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj)) zcca = z0+za*cangvg zccb = zcorl2(ji,jj)+za*zsang v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask END DO END DO CALL lbc_lnk( v_ice(:,:), 'V', -1. ) DO jj = k_j1+1, k_jpj-1 DO ji = 2, jpim1 zmask = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj) zsang = SIGN(1.0,fcor(ji,jj))*sangvg z0 = zmass1(ji,jj)/dtevp ! SB modif because ocean has no slip boundary condition zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) za = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj)) zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj)) zcca = z0+za*cangvg zccb = zcorl1(ji,jj)+za*zsang u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask END DO ! ji END DO ! jj CALL lbc_lnk( u_ice(:,:), 'U', -1. ) ENDIF !--- Convergence test. DO jj = k_j1+1 , k_jpj-1 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , & ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) END DO zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) ! ! ==================== ! END DO ! end loop over iter ! ! ! ==================== ! ! !------------------------------------------------------------------------------! ! 4) Prevent ice velocities when the ice is thin !------------------------------------------------------------------------------! ! ! If the ice thickness is below 1cm then ice velocity should equal the ! ocean velocity, ! This prevents high velocity when ice is thin DO jj = k_j1+1, k_jpj-1 DO ji = 2, jpim1 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) IF ( zdummy .LE. 5.0e-2 ) THEN u_ice(ji,jj) = u_io(ji,jj) v_ice(ji,jj) = v_io(ji,jj) ENDIF ! zdummy END DO END DO ! !------------------------------------------------------------------------------! ! 5) Compute stress rain invariants and store strain tensor !------------------------------------------------------------------------------! ! ! Compute delta, shear and div, inputs for mechanical redistribution DO jj = k_j1+1, k_jpj-1 DO ji = 2, jpim1 !- zdd(:,:), zdt(:,:): divergence and tension at centre !- zds(:,:): shear on northeast corner of grid cells zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) IF ( zdummy .LE. 5.0e-2 ) THEN zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & & -e2u(ji-1,jj)*u_ice(ji-1,jj) & & +e1v(ji,jj)*v_ice(ji,jj) & & -e1v(ji,jj-1)*v_ice(ji,jj-1) & & ) & & / area(ji,jj) zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) & & -u_ice(ji-1,jj)/e2u(ji-1,jj) & & )*e2t(ji,jj)*e2t(ji,jj) & & -( v_ice(ji,jj)/e1v(ji,jj) & & -v_ice(ji,jj-1)/e1v(ji,jj-1) & & )*e1t(ji,jj)*e1t(ji,jj) & & ) & & / area(ji,jj) ! ! SB modif because ocean has no slip boundary condition zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) & & - u_ice(ji,jj) / e1u(ji,jj) ) & & * e1f(ji,jj) * e1f(ji,jj) & & + ( v_ice(ji+1,jj) / e2v(ji+1,jj) & & - v_ice(ji,jj) / e2v(ji,jj) ) & & * e2f(ji,jj) * e2f(ji,jj) ) & & / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) & & * tmi(ji,jj) * tmi(ji,jj+1) & & * tmi(ji+1,jj) * tmi(ji+1,jj+1) !- Calculate Delta at centre of grid cells v_ice1(ji,jj) = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) u_ice2(ji,jj) = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1) & & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) zdst = ( e2u( ji , jj ) * v_ice1(ji,jj) & & - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & & + e1v( ji , jj ) * u_ice2(ji,jj) & & - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) & & ) & & / area(ji,jj) deltat(ji,jj) = SQRT( zdd(ji,jj)*zdd(ji,jj) + & & ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 & & ) + creepl divu_i(ji,jj) = zdd(ji,jj) delta_i(ji,jj) = deltat(ji,jj) shear_i(ji,jj) = zds(ji,jj) ENDIF ! zdummy END DO !jj END DO !ji ! MV This should be removed as it is the same as previous lines ! ! dynamical invariants ! IF ( stress_mean_swi .EQ. 0 ) THEN ! the following should not be necessary DO jj = k_j1+1, k_jpj-1 DO ji = 2, jpim1 divu_i(ji,jj) = zdd(ji,jj) delta_i(ji,jj) = deltat (ji,jj) shear_i(ji,jj) = zds (ji,jj) END DO END DO ! ENDIF CALL lbc_lnk( divu_i(:,:) , 'T', 1. ) CALL lbc_lnk( delta_i(:,:), 'T', 1. ) CALL lbc_lnk( shear_i(:,:), 'F', 1. ) ! Store the stress tensor for next time step ! this accelerates convergence and improves stability IF ( stress_mean_swi .EQ. 1 ) THEN DO jj = k_j1+1, k_jpj-1 DO ji = 2, jpim1 stress1_i(ji,jj) = zs1(ji,jj) stress2_i(ji,jj) = zs2(ji,jj) stress12_i(ji,jj) = zs12(ji,jj) END DO END DO ENDIF ! Lateral boundary condition CALL lbc_lnk( stress1_i(:,:), 'T', 1. ) CALL lbc_lnk( stress2_i(:,:), 'T', 1. ) CALL lbc_lnk( stress12_i(:,:), 'F', 1. ) ! !------------------------------------------------------------------------------! ! 6) Control prints of residual and charge ellipse !------------------------------------------------------------------------------! ! ! print the residual for convergence IF(ln_ctl) THEN WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, iter CALL prt_ctl_info(charout) CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :') ENDIF ! print charge ellipse ! This can be desactivated once the user is sure that the stress state ! lie on the charge ellipse. See Bouillon et al. 08 for more details IF(ln_ctl) THEN CALL prt_ctl_info('lim_rhg : numit :',ivar1=numit) CALL prt_ctl_info('lim_rhg : nwrite :',ivar1=nwrite) CALL prt_ctl_info('lim_rhg : MOD :',ivar1=MOD(numit,nwrite)) IF( MOD(numit,nwrite) .EQ. 0 ) THEN WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. ' CALL prt_ctl_info(charout) DO jj = k_j1+1, k_jpj-1 DO ji = 2, jpim1 IF (zpresh(ji,jj) .GT. 1.0) THEN sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) WRITE(charout,FMT="('lim_rhg :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)") CALL prt_ctl_info(charout) ENDIF END DO END DO WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 2000, numit, 0, 0, ' ch. ell. ' CALL prt_ctl_info(charout) ENDIF ENDIF END SUBROUTINE lim_rhg #else !!---------------------------------------------------------------------- !! Default option Dummy module NO LIM sea-ice model !!---------------------------------------------------------------------- CONTAINS SUBROUTINE lim_rhg( k1 , k2 ) ! Dummy routine WRITE(*,*) 'lim_rhg: You should not have seen this print! error?', k1, k2 END SUBROUTINE lim_rhg #endif !!============================================================================== END MODULE limrhg