MODULE icerhg_evp !!====================================================================== !! *** MODULE icerhg_evp *** !! Ice rheology : sea ice rheology !!====================================================================== !! History : - ! 2007-03 (M.A. Morales Maqueda, S. Bouillon) Original code !! 3.0 ! 2008-03 (M. Vancoppenolle) LIM3 !! - ! 2008-11 (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy !! 3.3 ! 2009-05 (G.Garric) addition of the evp cas !! 3.4 ! 2011-01 (A. Porter) dynamical allocation !! 3.5 ! 2012-08 (R. Benshila) AGRIF !! 3.6 ! 2016-06 (C. Rousset) Rewriting + landfast ice + possibility to use mEVP (Bouillon 2013) !!---------------------------------------------------------------------- #if defined key_lim3 !!---------------------------------------------------------------------- !! 'key_lim3' LIM-3 sea-ice model !!---------------------------------------------------------------------- !! ice_rhg_evp : computes ice velocities !!---------------------------------------------------------------------- USE phycst ! Physical constant USE par_oce ! Ocean parameters USE dom_oce ! Ocean domain USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b USE ice ! sea-ice: ice variables USE icerdgrft ! sea-ice: ice strength ! USE lbclnk ! Lateral Boundary Condition / MPP link USE lib_mpp ! MPP library USE in_out_manager ! I/O manager USE prtctl ! Print control USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) USE iom ! #if defined key_agrif USE agrif_lim3_interp #endif USE bdy_oce , ONLY: ln_bdy USE bdyice IMPLICIT NONE PRIVATE PUBLIC ice_rhg_evp ! routine called by icerhg.F90 !! * Substitutions # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/ICE 4.0 , NEMO Consortium (2017) !! $Id: icerhg_evp.F90 8378 2017-07-26 13:55:59Z clem $ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE ice_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i ) !!------------------------------------------------------------------- !! *** SUBROUTINE ice_rhg_evp *** !! 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 !! nn_nevp, elastic time scale and rn_creepl maintain stress state !! on the charge ellipse for plastic flow !! e.g. in the Canadian Archipelago !! !! ** Notes : There is the possibility to use mEVP from Bouillon 2013 !! (by uncommenting some lines in part 3 and changing alpha and beta parameters) !! but this solution appears very unstable (see Kimmritz et al 2016) !! !! References : Hunke and Dukowicz, JPO97 !! Bouillon et al., Ocean Modelling 2009 !! Bouillon et al., Ocean Modelling 2013 !!------------------------------------------------------------------- INTEGER, INTENT(in) :: kt ! time step REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i !! INTEGER :: ji, jj ! dummy loop indices INTEGER :: jter ! local integers REAL(wp) :: zrhoco ! rau0 * rn_cio REAL(wp) :: zdtevp, z1_dtevp ! time step for subcycling REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity REAL(wp) :: zbeta, zalph1, z1_alph1, zalph2, z1_alph2 ! alpha and beta from Bouillon 2009 and 2013 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV ! ice/snow mass REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars REAL(wp) :: zTauO, zTauB, zTauE, zvel ! temporary scalars REAL(wp) :: zresm ! Maximal error on ice velocity REAL(wp) :: zintb, zintn ! dummy argument REAL(wp) :: zfac_x, zfac_y REAL(wp) :: zshear, zdum1, zdum2 REAL(wp), DIMENSION(jpi,jpj) :: z1_e1t0, z1_e2t0 ! scale factors REAL(wp), DIMENSION(jpi,jpj) :: zp_delt ! P/delta at T points ! REAL(wp), DIMENSION(jpi,jpj) :: zaU , zaV ! ice fraction on U/V points REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! ice/snow mass/dt on U/V points REAL(wp), DIMENSION(jpi,jpj) :: zmf ! coriolis parameter at T points REAL(wp), DIMENSION(jpi,jpj) :: zTauU_ia , ztauV_ia ! ice-atm. stress at U-V points REAL(wp), DIMENSION(jpi,jpj) :: zspgU , zspgV ! surface pressure gradient at U/V points REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points REAL(wp), DIMENSION(jpi,jpj) :: zfU , zfV ! internal stresses REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice, zresr ! check convergence REAL(wp), DIMENSION(jpi,jpj) :: zpice ! array used for the calculation of ice surface slope: ! ocean surface (ssh_m) if ice is not embedded ! ice top surface if ice is embedded REAL(wp), DIMENSION(jpi,jpj) :: zCorx, zCory ! Coriolis stress array REAL(wp), DIMENSION(jpi,jpj) :: ztaux_oi, ztauy_oi ! Ocean-to-ice stress array REAL(wp), DIMENSION(jpi,jpj) :: zswitchU, zswitchV ! dummy arrays REAL(wp), DIMENSION(jpi,jpj) :: zmaskU, zmaskV ! mask for ice presence REAL(wp), DIMENSION(jpi,jpj) :: zfmask, zwf ! mask at F points for the ice REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity equals ocean velocity !! --- diags REAL(wp), DIMENSION(jpi,jpj) :: zswi REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig1, zsig2, zsig3 !! --- SIMIP diags REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_sig1 ! Average normal stress in sea ice REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_sig2 ! Maximum shear stress in sea ice REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_dssh_dx ! X-direction sea-surface tilt term (N/m2) REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_dssh_dy ! X-direction sea-surface tilt term (N/m2) REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_corstrx ! X-direction coriolis stress (N/m2) REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_corstry ! Y-direction coriolis stress (N/m2) REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_intstrx ! X-direction internal stress (N/m2) REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_intstry ! Y-direction internal stress (N/m2) REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_utau_oi ! X-direction ocean-ice stress REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_vtau_oi ! Y-direction ocean-ice stress REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s) REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s) REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s) REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xatrp ! X-component of area transport (m2/s) REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_yatrp ! Y-component of area transport (m2/s) !!------------------------------------------------------------------- IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_rhg_evp: EVP sea-ice rheology' #if defined key_agrif CALL agrif_interp_lim3( 'U', 0, nn_nevp ) ! First interpolation of coarse values CALL agrif_interp_lim3( 'V', 0, nn_nevp ) #endif ! !------------------------------------------------------------------------------! ! 0) mask at F points for the ice !------------------------------------------------------------------------------! ! ocean/land mask DO jj = 1, jpjm1 DO ji = 1, jpim1 ! NO vector opt. zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) END DO END DO CALL lbc_lnk( zfmask, 'F', 1._wp ) ! Lateral boundary conditions on velocity (modify zfmask) zwf(:,:) = zfmask(:,:) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. IF( zfmask(ji,jj) == 0._wp ) THEN zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) ) ENDIF END DO END DO DO jj = 2, jpjm1 IF( zfmask(1,jj) == 0._wp ) THEN zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) ) ENDIF IF( zfmask(jpi,jj) == 0._wp ) THEN zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) ) ENDIF END DO DO ji = 2, jpim1 IF( zfmask(ji,1) == 0._wp ) THEN zfmask(ji,1 ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) ) ENDIF IF( zfmask(ji,jpj) == 0._wp ) THEN zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) ) ENDIF END DO CALL lbc_lnk( zfmask, 'F', 1._wp ) !------------------------------------------------------------------------------! ! 1) define some variables and initialize arrays !------------------------------------------------------------------------------! zrhoco = rau0 * rn_cio ! ecc2: square of yield ellipse eccenticrity ecc2 = rn_ecc * rn_ecc z1_ecc2 = 1._wp / ecc2 ! Time step for subcycling zdtevp = rdt_ice / REAL( nn_nevp ) z1_dtevp = 1._wp / zdtevp ! alpha parameters (Bouillon 2009) zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp zalph2 = zalph1 * z1_ecc2 ! alpha and beta parameters (Bouillon 2013) !!zalph1 = 40. !!zalph2 = 40. !!zbeta = 3000. !!zbeta = REAL( nn_nevp ) ! close to classical EVP of Hunke (2001) z1_alph1 = 1._wp / ( zalph1 + 1._wp ) z1_alph2 = 1._wp / ( zalph2 + 1._wp ) ! Initialise stress tensor zs1 (:,:) = pstress1_i (:,:) zs2 (:,:) = pstress2_i (:,:) zs12(:,:) = pstress12_i(:,:) ! Ice strength CALL ice_rdgrft_strength ! scale factors DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj ) + e1t(ji,jj ) ) z1_e2t0(ji,jj) = 1._wp / ( e2t(ji ,jj+1) + e2t(ji,jj ) ) END DO END DO ! !------------------------------------------------------------------------------! ! 2) Wind / ocean stress, mass terms, coriolis terms !------------------------------------------------------------------------------! IF( ln_ice_embd ) THEN !== embedded sea ice: compute representative ice top surface ==! ! ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} ! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1} zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp ! ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} ! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp ! zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0 ! ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==! zpice(:,:) = ssh_m(:,:) ENDIF DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! ice fraction at U-V points zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) ! Ice/snow mass at U-V points zm1 = ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) ) zm2 = ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) ) zm3 = ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) ) zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) ! Ocean currents at U-V points v_oceU(ji,jj) = 0.5_wp * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji+1,jj) & & + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) u_oceV(ji,jj) = 0.5_wp * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj+1) & & + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) ! Coriolis at T points (m*f) zmf(ji,jj) = zm1 * ff_t(ji,jj) ! m/dt zmU_t(ji,jj) = zmassU * z1_dtevp zmV_t(ji,jj) = zmassV * z1_dtevp ! Drag ice-atm. zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points zspgU(ji,jj) = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) zspgV(ji,jj) = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) ! masks zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice ! switches zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin END DO END DO CALL lbc_lnk( zmf, 'T', 1. ) ! !------------------------------------------------------------------------------! ! 3) Solution of the momentum equation, iterative procedure !------------------------------------------------------------------------------! ! ! !----------------------! DO jter = 1 , nn_nevp ! loop over jter ! ! !----------------------! IF(ln_ctl) THEN ! Convergence test DO jj = 1, jpjm1 zu_ice(:,jj) = pu_ice(:,jj) ! velocity at previous time step zv_ice(:,jj) = pv_ice(:,jj) END DO ENDIF ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! DO jj = 1, jpjm1 ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points DO ji = 1, jpim1 ! shear at F points zds(ji,jj) = ( ( pu_ice(ji,jj+1) * r1_e1u(ji,jj+1) - pu_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & & + ( pv_ice(ji+1,jj) * r1_e2v(ji+1,jj) - pv_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) END DO END DO CALL lbc_lnk( zds, 'F', 1. ) DO jj = 2, jpjm1 DO ji = 2, jpim1 ! no vector loop ! shear**2 at T points (doc eq. A16) zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & & ) * 0.25_wp * r1_e1e2t(ji,jj) ! divergence at T points zdiv = ( e2u(ji,jj) * pu_ice(ji,jj) - e2u(ji-1,jj) * pu_ice(ji-1,jj) & & + e1v(ji,jj) * pv_ice(ji,jj) - e1v(ji,jj-1) * pv_ice(ji,jj-1) & & ) * r1_e1e2t(ji,jj) zdiv2 = zdiv * zdiv ! tension at T points zdt = ( ( pu_ice(ji,jj) * r1_e2u(ji,jj) - pu_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & & - ( pv_ice(ji,jj) * r1_e1v(ji,jj) - pv_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & & ) * r1_e1e2t(ji,jj) zdt2 = zdt * zdt ! delta at T points zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) ! P/delta at T points zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) ! stress at T points zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1 zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2 END DO END DO CALL lbc_lnk( zp_delt, 'T', 1. ) DO jj = 1, jpjm1 DO ji = 1, jpim1 ! P/delta at F points zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) ! stress at F points zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2 END DO END DO CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! !--- U points zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) & & ) * r1_e2u(ji,jj) & & + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) & & ) * 2._wp * r1_e1u(ji,jj) & & ) * r1_e1e2u(ji,jj) ! ! !--- V points zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) & & ) * r1_e1v(ji,jj) & & + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) & & ) * 2._wp * r1_e2v(ji,jj) & & ) * r1_e1e2v(ji,jj) ! ! !--- u_ice at V point u_iceV(ji,jj) = 0.5_wp * ( ( pu_ice(ji,jj ) + pu_ice(ji-1,jj ) ) * e2t(ji,jj+1) & & + ( pu_ice(ji,jj+1) + pu_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) ! ! !--- v_ice at U point v_iceU(ji,jj) = 0.5_wp * ( ( pv_ice(ji ,jj) + pv_ice(ji ,jj-1) ) * e1t(ji+1,jj) & & + ( pv_ice(ji+1,jj) + pv_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) ! END DO END DO ! ! --- Computation of ice velocity --- ! ! Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp ! Bouillon et al. 2009 (eq 34-35) => stable IF( MOD(jter,2) == 0 ) THEN ! even iterations ! DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! !--- tau_io/(v_oce - v_ice) zTauO = zaV(ji,jj) * zrhoco * SQRT( ( pv_ice(ji,jj) - v_oce (ji,jj) ) * ( pv_ice(ji,jj) - v_oce (ji,jj) ) & & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) ! !--- Ocean-to-Ice stress ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - pv_ice(ji,jj) ) ! ! !--- tau_bottom/v_ice zvel = MAX( zepsi, SQRT( pv_ice(ji,jj) * pv_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) zTauB = - tau_icebfr(ji,jj) / zvel ! ! !--- Coriolis at V-points (energy conserving formulation) zCory(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & & ( zmf(ji,jj ) * ( e2u(ji,jj ) * pu_ice(ji,jj ) + e2u(ji-1,jj ) * pu_ice(ji-1,jj ) ) & & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * pu_ice(ji,jj+1) + e2u(ji-1,jj+1) * pu_ice(ji-1,jj+1) ) ) ! ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) ! ! !--- landfast switch => 0 = static friction ; 1 = sliding friction rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) ! ! !--- ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) pv_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * pv_ice(ji,jj) & ! previous velocity & + zTauE + zTauO * pv_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast & + ( 1._wp - rswitch ) * pv_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin & ) * zmaskV(ji,jj) ! ! Bouillon 2013 !!pv_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * pv_ice(ji,jj) + pv_ice_b(ji,jj) ) & !! & + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj) & !! & ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) ! END DO END DO CALL lbc_lnk( pv_ice, 'V', -1. ) ! #if defined key_agrif !! CALL agrif_interp_lim3( 'V', jter, nn_nevp ) CALL agrif_interp_lim3( 'V' ) #endif IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) ! DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! tau_io/(u_oce - u_ice) zTauO = zaU(ji,jj) * zrhoco * SQRT( ( pu_ice(ji,jj) - u_oce (ji,jj) ) * ( pu_ice(ji,jj) - u_oce (ji,jj) ) & & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) ! Ocean-to-Ice stress ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - pu_ice(ji,jj) ) ! tau_bottom/u_ice zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + pu_ice(ji,jj) * pu_ice(ji,jj) ) ) zTauB = - tau_icebfr(ji,jj) / zvel ! Coriolis at U-points (energy conserving formulation) zCorx(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & & ( zmf(ji ,jj) * ( e1v(ji ,jj) * pv_ice(ji ,jj) + e1v(ji ,jj-1) * pv_ice(ji ,jj-1) ) & & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * pv_ice(ji+1,jj) + e1v(ji+1,jj-1) * pv_ice(ji+1,jj-1) ) ) ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) ! landfast switch => 0 = static friction ; 1 = sliding friction rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) pu_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * pu_ice(ji,jj) & ! previous velocity & + zTauE + zTauO * pu_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast & + ( 1._wp - rswitch ) * pu_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin & ) * zmaskU(ji,jj) ! Bouillon 2013 !!pu_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * pu_ice(ji,jj) + pu_ice_b(ji,jj) ) & !! & + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj) & !! & ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) END DO END DO CALL lbc_lnk( pu_ice, 'U', -1. ) ! #if defined key_agrif !! CALL agrif_interp_lim3( 'U', jter, nn_nevp ) CALL agrif_interp_lim3( 'U' ) #endif IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) ! ELSE ! odd iterations ! DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! tau_io/(u_oce - u_ice) zTauO = zaU(ji,jj) * zrhoco * SQRT( ( pu_ice(ji,jj) - u_oce (ji,jj) ) * ( pu_ice(ji,jj) - u_oce (ji,jj) ) & & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) ! Ocean-to-Ice stress ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - pu_ice(ji,jj) ) ! tau_bottom/u_ice zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + pu_ice(ji,jj) * pu_ice(ji,jj) ) ) zTauB = - tau_icebfr(ji,jj) / zvel ! Coriolis at U-points (energy conserving formulation) zCorx(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & & ( zmf(ji ,jj) * ( e1v(ji ,jj) * pv_ice(ji ,jj) + e1v(ji ,jj-1) * pv_ice(ji ,jj-1) ) & & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * pv_ice(ji+1,jj) + e1v(ji+1,jj-1) * pv_ice(ji+1,jj-1) ) ) ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) ! landfast switch => 0 = static friction ; 1 = sliding friction rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) pu_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * pu_ice(ji,jj) & ! previous velocity & + zTauE + zTauO * pu_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast & + ( 1._wp - rswitch ) * pu_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin & ) * zmaskU(ji,jj) ! Bouillon 2013 !!pu_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * pu_ice(ji,jj) + pu_ice_b(ji,jj) ) & !! & + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj) & !! & ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) END DO END DO CALL lbc_lnk( pu_ice, 'U', -1. ) ! #if defined key_agrif !! CALL agrif_interp_lim3( 'U', jter, nn_nevp ) CALL agrif_interp_lim3( 'U' ) #endif IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) ! DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! tau_io/(v_oce - v_ice) zTauO = zaV(ji,jj) * zrhoco * SQRT( ( pv_ice(ji,jj) - v_oce (ji,jj) ) * ( pv_ice(ji,jj) - v_oce (ji,jj) ) & & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) ! Ocean-to-Ice stress ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - pv_ice(ji,jj) ) ! tau_bottom/v_ice zvel = MAX( zepsi, SQRT( pv_ice(ji,jj) * pv_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) ztauB = - tau_icebfr(ji,jj) / zvel ! Coriolis at V-points (energy conserving formulation) zCory(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & & ( zmf(ji,jj ) * ( e2u(ji,jj ) * pu_ice(ji,jj ) + e2u(ji-1,jj ) * pu_ice(ji-1,jj ) ) & & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * pu_ice(ji,jj+1) + e2u(ji-1,jj+1) * pu_ice(ji-1,jj+1) ) ) ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) pv_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * pv_ice(ji,jj) & ! previous velocity & + zTauE + zTauO * pv_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast & + ( 1._wp - rswitch ) * pv_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin & ) * zmaskV(ji,jj) ! Bouillon 2013 !!pv_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * pv_ice(ji,jj) + pv_ice_b(ji,jj) ) & !! & + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj) & !! & ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) END DO END DO CALL lbc_lnk( pv_ice, 'V', -1. ) ! #if defined key_agrif !! CALL agrif_interp_lim3( 'V', jter, nn_nevp ) CALL agrif_interp_lim3( 'V' ) #endif IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) ! ENDIF IF(ln_ctl) THEN ! Convergence test DO jj = 2 , jpjm1 zresr(:,jj) = MAX( ABS( pu_ice(:,jj) - zu_ice(:,jj) ), ABS( pv_ice(:,jj) - zv_ice(:,jj) ) ) END DO zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain ENDIF ! ! ! ==================== ! END DO ! end loop over jter ! ! ! ==================== ! ! !------------------------------------------------------------------------------! ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) !------------------------------------------------------------------------------! DO jj = 1, jpjm1 DO ji = 1, jpim1 ! shear at F points zds(ji,jj) = ( ( pu_ice(ji,jj+1) * r1_e1u(ji,jj+1) - pu_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & & + ( pv_ice(ji+1,jj) * r1_e2v(ji+1,jj) - pv_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) END DO END DO CALL lbc_lnk( zds, 'F', 1. ) DO jj = 2, jpjm1 DO ji = 2, jpim1 ! no vector loop ! tension**2 at T points zdt = ( ( pu_ice(ji,jj) * r1_e2u(ji,jj) - pu_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & & - ( pv_ice(ji,jj) * r1_e1v(ji,jj) - pv_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & & ) * r1_e1e2t(ji,jj) zdt2 = zdt * zdt ! shear**2 at T points (doc eq. A16) zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & & ) * 0.25_wp * r1_e1e2t(ji,jj) ! shear at T points pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) ! divergence at T points pdivu_i(ji,jj) = ( e2u(ji,jj) * pu_ice(ji,jj) - e2u(ji-1,jj) * pu_ice(ji-1,jj) & & + e1v(ji,jj) * pv_ice(ji,jj) - e1v(ji,jj-1) * pv_ice(ji,jj-1) & & ) * r1_e1e2t(ji,jj) ! delta at T points zdelta = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch END DO END DO CALL lbc_lnk_multi( pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) ! --- Store the stress tensor for the next time step --- ! pstress1_i (:,:) = zs1 (:,:) pstress2_i (:,:) = zs2 (:,:) pstress12_i(:,:) = zs12(:,:) ! !------------------------------------------------------------------------------! ! 5) diagnostics !------------------------------------------------------------------------------! DO jj = 1, jpj DO ji = 1, jpi zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice END DO END DO ! --- divergence, shear and strength --- ! IF( iom_use('idive') ) CALL iom_put( "idive" , pdivu_i(:,:) * zswi(:,:) ) ! divergence IF( iom_use('ishear') ) CALL iom_put( "ishear" , pshear_i(:,:) * zswi(:,:) ) ! shear IF( iom_use('icestr') ) CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) ) ! Ice strength ! --- charge ellipse --- ! IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN ! ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) ! DO jj = 2, jpjm1 DO ji = 2, jpim1 zdum1 = ( zswi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point & zswi(ji ,jj) * pstress12_i(ji ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & & / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) ) zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) ) !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) !! zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) !! zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) END DO END DO CALL lbc_lnk_multi( zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) ! IF( iom_use('isig1') ) CALL iom_put( "isig1" , zsig1 ) IF( iom_use('isig2') ) CALL iom_put( "isig2" , zsig2 ) IF( iom_use('isig3') ) CALL iom_put( "isig3" , zsig3 ) ! DEALLOCATE( zsig1 , zsig2 , zsig3 ) ENDIF ! --- SIMIP --- ! IF ( iom_use( 'normstr' ) .OR. iom_use( 'sheastr' ) .OR. iom_use( 'dssh_dx' ) .OR. iom_use( 'dssh_dy' ) .OR. & & iom_use( 'corstrx' ) .OR. iom_use( 'corstry' ) .OR. iom_use( 'intstrx' ) .OR. iom_use( 'intstry' ) .OR. & & iom_use( 'utau_oi' ) .OR. iom_use( 'vtau_oi' ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. & & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp' ) .OR. iom_use( 'yatrp' ) ) THEN ALLOCATE( zdiag_sig1 (jpi,jpj) , zdiag_sig2 (jpi,jpj) , zdiag_dssh_dx (jpi,jpj) , zdiag_dssh_dy (jpi,jpj) , & & zdiag_corstrx (jpi,jpj) , zdiag_corstry (jpi,jpj) , zdiag_intstrx (jpi,jpj) , zdiag_intstry (jpi,jpj) , & & zdiag_utau_oi (jpi,jpj) , zdiag_vtau_oi (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , & & zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp (jpi,jpj) , zdiag_yatrp (jpi,jpj) ) DO jj = 2, jpjm1 DO ji = 2, jpim1 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice ! Stress tensor invariants (normal and shear stress N/m) zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch ! normal stress zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch ! shear stress ! Stress terms of the momentum equation (N/m2) zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch ! sea surface slope stress term zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch ! Coriolis stress term zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch zdiag_intstrx(ji,jj) = zfU(ji,jj) * rswitch ! internal stress term zdiag_intstry(ji,jj) = zfV(ji,jj) * rswitch zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch ! oceanic stress zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch ! 2D ice mass, snow mass, area transport arrays (X, Y) zfac_x = 0.5 * pu_ice(ji,jj) * e2u(ji,jj) * rswitch zfac_y = 0.5 * pv_ice(ji,jj) * e1v(ji,jj) * rswitch zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component zdiag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) ! '' Y- '' zdiag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component zdiag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) ! '' Y- '' zdiag_xatrp(ji,jj) = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) ) ! area transport, X-component zdiag_yatrp(ji,jj) = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) ) ! '' Y- '' END DO END DO CALL lbc_lnk_multi( zdiag_sig1 , 'T', 1., zdiag_sig2 , 'T', 1., & & zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1., & & zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1., & & zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1. ) CALL lbc_lnk_multi( zdiag_utau_oi, 'U', -1., zdiag_vtau_oi, 'V', -1. ) CALL lbc_lnk_multi( zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1., & & zdiag_xatrp , 'U', -1., zdiag_ymtrp_ice, 'V', -1., & & zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp , 'V', -1. ) IF ( iom_use( 'normstr' ) ) CALL iom_put( 'normstr' , zdiag_sig1(:,:) ) ! Normal stress IF ( iom_use( 'sheastr' ) ) CALL iom_put( 'sheastr' , zdiag_sig2(:,:) ) ! Shear stress IF ( iom_use( 'dssh_dx' ) ) CALL iom_put( 'dssh_dx' , zdiag_dssh_dx(:,:) ) ! Sea-surface tilt term in force balance (x) IF ( iom_use( 'dssh_dy' ) ) CALL iom_put( 'dssh_dy' , zdiag_dssh_dy(:,:) ) ! Sea-surface tilt term in force balance (y) IF ( iom_use( 'corstrx' ) ) CALL iom_put( 'corstrx' , zdiag_corstrx(:,:) ) ! Coriolis force term in force balance (x) IF ( iom_use( 'corstry' ) ) CALL iom_put( 'corstry' , zdiag_corstry(:,:) ) ! Coriolis force term in force balance (y) IF ( iom_use( 'intstrx' ) ) CALL iom_put( 'intstrx' , zdiag_intstrx(:,:) ) ! Internal force term in force balance (x) IF ( iom_use( 'intstry' ) ) CALL iom_put( 'intstry' , zdiag_intstry(:,:) ) ! Internal force term in force balance (y) IF ( iom_use( 'utau_oi' ) ) CALL iom_put( 'utau_oi' , zdiag_utau_oi(:,:) ) ! Ocean stress term in force balance (x) IF ( iom_use( 'vtau_oi' ) ) CALL iom_put( 'vtau_oi' , zdiag_vtau_oi(:,:) ) ! Ocean stress term in force balance (y) IF ( iom_use( 'xmtrpice' ) ) CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice(:,:) ) ! X-component of sea-ice mass transport (kg/s) IF ( iom_use( 'ymtrpice' ) ) CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice(:,:) ) ! Y-component of sea-ice mass transport IF ( iom_use( 'xmtrpsnw' ) ) CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw(:,:) ) ! X-component of snow mass transport (kg/s) IF ( iom_use( 'ymtrpsnw' ) ) CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw(:,:) ) ! Y-component of snow mass transport IF ( iom_use( 'xatrp' ) ) CALL iom_put( 'xatrp' , zdiag_xatrp(:,:) ) ! X-component of ice area transport IF ( iom_use( 'yatrp' ) ) CALL iom_put( 'yatrp' , zdiag_yatrp(:,:) ) ! Y-component of ice area transport DEALLOCATE( zdiag_sig1 , zdiag_sig2 , zdiag_dssh_dx , zdiag_dssh_dy , & & zdiag_corstrx , zdiag_corstry , zdiag_intstrx , zdiag_intstry , & & zdiag_utau_oi , zdiag_vtau_oi , zdiag_xmtrp_ice , zdiag_ymtrp_ice , & & zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp ) ENDIF ! END SUBROUTINE ice_rhg_evp #else !!---------------------------------------------------------------------- !! Default option Empty module NO LIM-3 sea-ice model !!---------------------------------------------------------------------- #endif !!============================================================================== END MODULE icerhg_evp