MODULE icedyn_rhg_evp !!====================================================================== !! *** MODULE icedyn_rhg_evp *** !! Sea-Ice dynamics : rheology Elasto-Viscous-Plastic !!====================================================================== !! 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 case !! 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) !! 3.7 ! 2017 (C. Rousset) add aEVP (Kimmritz 2016-2017) !! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube] !!---------------------------------------------------------------------- #if defined key_si3 !!---------------------------------------------------------------------- !! 'key_si3' SI3 sea-ice model !!---------------------------------------------------------------------- !! ice_dyn_rhg_evp : computes ice velocities from EVP rheology !! rhg_evp_rst : read/write EVP fields in ice restart !!---------------------------------------------------------------------- USE phycst ! Physical constant 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 icedyn_rdgrft ! sea-ice: ice strength USE bdy_oce , ONLY : ln_bdy USE bdyice #if defined key_agrif USE agrif_ice_interp #endif ! USE in_out_manager ! I/O manager USE iom ! I/O manager library USE lib_mpp ! MPP library USE lib_fortran ! fortran utilities (glob_sum + no signed zero) USE lbclnk ! lateral boundary conditions (or mpp links) USE prtctl ! Print control IMPLICIT NONE PRIVATE PUBLIC ice_dyn_rhg_evp ! called by icedyn_rhg.F90 PUBLIC rhg_evp_rst ! called by icedyn_rhg.F90 !! * Substitutions # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/ICE 4.0 , NEMO Consortium (2018) !! $Id: icedyn_rhg_evp.F90 8378 2017-07-26 13:55:59Z clem $ !! Software governed by the CeCILL licence (./LICENSE) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i ) !!------------------------------------------------------------------- !! *** SUBROUTINE ice_dyn_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 : 0) compute mask at F point !! 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) Recompute delta, shear and divergence !! (which are inputs of the ITD) & store stress !! for the next time step !! 5) Diagnostics including charge ellipse !! !! ** Notes : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017) !! by setting up ln_aEVP=T (i.e. changing alpha and beta parameters). !! This is an upgraded version of mEVP from Bouillon et al. 2013 !! (i.e. more stable and better convergence) !! !! References : Hunke and Dukowicz, JPO97 !! Bouillon et al., Ocean Modelling 2009 !! Bouillon et al., Ocean Modelling 2013 !! Kimmritz et al., Ocean Modelling 2016 & 2017 !!------------------------------------------------------------------- INTEGER , INTENT(in ) :: kt ! time step REAL(wp), DIMENSION(:,:), INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i ! REAL(wp), DIMENSION(:,:), INTENT( out) :: 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) :: zalph1, z1_alph1, zalph2, z1_alph2 ! alpha coef from Bouillon 2009 or Kimmritz 2017 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) :: zbeta ! beta coef from Kimmritz 2017 ! REAL(wp), DIMENSION(jpi,jpj) :: zdt_m ! (dt / ice-snow_mass) on 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_dyn_rhg_evp: EVP sea-ice rheology' ! !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... !------------------------------------------------------------------------------! ! 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) IF( .NOT. ln_aEVP ) THEN zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp zalph2 = zalph1 * z1_ecc2 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) z1_alph2 = 1._wp / ( zalph2 + 1._wp ) ENDIF ! Initialise stress tensor zs1 (:,:) = pstress1_i (:,:) zs2 (:,:) = pstress2_i (:,:) zs12(:,:) = pstress12_i(:,:) ! Ice strength CALL ice_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) ! dt/m at T points (for alpha and beta coefficients) zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin ) ! 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_multi( zmf, 'T', 1., zdt_m, '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) = u_ice(:,jj) ! velocity at previous time step zv_ice(:,jj) = v_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) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_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, jpj ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 DO ji = 2, jpi ! 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) * 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) & & ) * r1_e1e2t(ji,jj) zdiv2 = zdiv * zdiv ! tension at T points zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_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 ) ! alpha & beta for aEVP ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m ! alpha = beta = sqrt(4*gamma) IF( ln_aEVP ) THEN zalph1 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) z1_alph1 = 1._wp / ( zalph1 + 1._wp ) zalph2 = zalph1 z1_alph2 = z1_alph1 ENDIF ! 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 ! alpha & beta for aEVP IF( ln_aEVP ) THEN zalph2 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) z1_alph2 = 1._wp / ( zalph2 + 1._wp ) zbeta(ji,jj) = zalph2 ENDIF ! 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 ! --- 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 * ( ( 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 ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) ! ! !--- v_ice at U point v_iceU(ji,jj) = 0.5_wp * ( ( 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) ) * 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 vary as in Kimmritz 2016 & 2017 ! 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( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_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) - v_ice(ji,jj) ) ! ! !--- tau_bottom/v_ice zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_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 ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_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 ) ) ) ! IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity & + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast & + ( 1._wp - rswitch ) * v_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) ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity & + zTauE + zTauO * v_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 ) * v_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) ENDIF END DO END DO CALL lbc_lnk( v_ice, 'V', -1. ) ! #if defined key_agrif !! CALL agrif_interp_ice( 'V', jter, nn_nevp ) CALL agrif_interp_ice( '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( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_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) - u_ice(ji,jj) ) ! ! !--- tau_bottom/u_ice zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_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) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_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 ) ) ) ! IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast & + ( 1._wp - rswitch ) * u_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) ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity & + zTauE + zTauO * u_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 ) * u_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) ENDIF END DO END DO CALL lbc_lnk( u_ice, 'U', -1. ) ! #if defined key_agrif !! CALL agrif_interp_ice( 'U', jter, nn_nevp ) CALL agrif_interp_ice( '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( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_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) - u_ice(ji,jj) ) ! ! !--- tau_bottom/u_ice zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_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) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_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 ) ) ) ! IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast & + ( 1._wp - rswitch ) * u_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) ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity & + zTauE + zTauO * u_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 ) * u_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) ENDIF END DO END DO CALL lbc_lnk( u_ice, 'U', -1. ) ! #if defined key_agrif !! CALL agrif_interp_ice( 'U', jter, nn_nevp ) CALL agrif_interp_ice( '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( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_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) - v_ice(ji,jj) ) ! ! !--- tau_bottom/v_ice zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_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 ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_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 ) ) ) ! IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity & + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast & + ( 1._wp - rswitch ) * v_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) ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity & + zTauE + zTauO * v_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 ) * v_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) ENDIF END DO END DO CALL lbc_lnk( v_ice, 'V', -1. ) ! #if defined key_agrif !! CALL agrif_interp_ice( 'V', jter, nn_nevp ) CALL agrif_interp_ice( '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( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_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) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) END DO END DO DO jj = 2, jpjm1 DO ji = 2, jpim1 ! no vector loop ! tension**2 at T points zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_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) * 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) & & ) * 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 --- ! CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 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('icediv') ) CALL iom_put( "icediv" , pdivu_i (:,:) * zswi(:,:) ) ! divergence IF( iom_use('iceshe') ) CALL iom_put( "iceshe" , 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 * u_ice(ji,jj) * e2u(ji,jj) * rswitch zfac_y = 0.5 * v_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., & & 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_dyn_rhg_evp SUBROUTINE rhg_evp_rst( cdrw, kt ) !!--------------------------------------------------------------------- !! *** ROUTINE rhg_evp_rst *** !! !! ** Purpose : Read or write RHG file in restart file !! !! ** Method : use of IOM library !!---------------------------------------------------------------------- CHARACTER(len=*) , INTENT(in) :: cdrw ! "READ"/"WRITE" flag INTEGER, OPTIONAL, INTENT(in) :: kt ! ice time-step ! INTEGER :: iter ! local integer INTEGER :: id1, id2, id3 ! local integers !!---------------------------------------------------------------------- ! IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialize ! ! --------------- IF( ln_rstart ) THEN !* Read the restart file ! id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. ) id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. ) id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. ) ! IF( MIN( id1, id2, id3 ) > 0 ) THEN ! fields exist CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i ) CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i ) CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i ) ELSE ! start rheology from rest IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' ==>>> previous run without rheology, set stresses to 0' stress1_i (:,:) = 0._wp stress2_i (:,:) = 0._wp stress12_i(:,:) = 0._wp ENDIF ELSE !* Start from rest IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set stresses to 0' stress1_i (:,:) = 0._wp stress2_i (:,:) = 0._wp stress12_i(:,:) = 0._wp ENDIF ! ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file ! ! ------------------- IF(lwp) WRITE(numout,*) '---- rhg-rst ----' iter = kt + nn_fsbc - 1 ! ice restarts are written at kt == nitrst - nn_fsbc + 1 ! CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i ) CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i ) CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i ) ! ENDIF ! END SUBROUTINE rhg_evp_rst #else !!---------------------------------------------------------------------- !! Default option Empty module NO SI3 sea-ice model !!---------------------------------------------------------------------- #endif !!============================================================================== END MODULE icedyn_rhg_evp