New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

icedyn_rhg_evp.F90 in NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_rhg_evp.F90 @ 13570

Last change on this file since 13570 was 13570, checked in by mocavero, 4 years ago
Added key_mpi3 to activate/deactivate the use of MPI3 neighborhood collectives lbc routines - ticket #2496
Property svn:keywords set to `Id`
File size: 56.1 KB

Rev	Line
[8586]	1	MODULE icedyn_rhg_evp
	2	!!======================================================================
	3	!! * MODULE icedyn_rhg_evp *
	4	!! Sea-Ice dynamics : rheology Elasto-Viscous-Plastic
	5	!!======================================================================
	6	!! History : - ! 2007-03 (M.A. Morales Maqueda, S. Bouillon) Original code
[9656]	7	!! 3.0 ! 2008-03 (M. Vancoppenolle) adaptation to new model
[8586]	8	!! - ! 2008-11 (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy
[9604]	9	!! 3.3 ! 2009-05 (G.Garric) addition of the evp case
	10	!! 3.4 ! 2011-01 (A. Porter) dynamical allocation
	11	!! 3.5 ! 2012-08 (R. Benshila) AGRIF
[9656]	12	!! 3.6 ! 2016-06 (C. Rousset) Rewriting + landfast ice + mEVP (Bouillon 2013)
[9604]	13	!! 3.7 ! 2017 (C. Rousset) add aEVP (Kimmritz 2016-2017)
	14	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
[8586]	15	!!----------------------------------------------------------------------
[9570]	16	#if defined key_si3
[8586]	17	!!----------------------------------------------------------------------
[9570]	18	!! 'key_si3' SI3 sea-ice model
[8586]	19	!!----------------------------------------------------------------------
[8813]	20	!! ice_dyn_rhg_evp : computes ice velocities from EVP rheology
	21	!! rhg_evp_rst : read/write EVP fields in ice restart
[8586]	22	!!----------------------------------------------------------------------
	23	USE phycst ! Physical constant
	24	USE dom_oce ! Ocean domain
	25	USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
	26	USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
	27	USE ice ! sea-ice: ice variables
[10332]	28	USE icevar ! ice_var_sshdyn
[8586]	29	USE icedyn_rdgrft ! sea-ice: ice strength
[8813]	30	USE bdy_oce , ONLY : ln_bdy
	31	USE bdyice
	32	#if defined key_agrif
[9596]	33	USE agrif_ice_interp
[8813]	34	#endif
[8586]	35	!
	36	USE in_out_manager ! I/O manager
	37	USE iom ! I/O manager library
	38	USE lib_mpp ! MPP library
	39	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
	40	USE lbclnk ! lateral boundary conditions (or mpp links)
	41	USE prtctl ! Print control
	42
	43	IMPLICIT NONE
	44	PRIVATE
	45
	46	PUBLIC ice_dyn_rhg_evp ! called by icedyn_rhg.F90
	47	PUBLIC rhg_evp_rst ! called by icedyn_rhg.F90
	48
	49	!! * Substitutions
[12377]	50	# include "do_loop_substitute.h90"
[13237]	51	# include "domzgr_substitute.h90"
[8586]	52	!!----------------------------------------------------------------------
[9598]	53	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	54	!! $Id$
[10068]	55	!! Software governed by the CeCILL license (see ./LICENSE)
[8586]	56	!!----------------------------------------------------------------------
	57	CONTAINS
	58
[12377]	59	SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )
[8586]	60	!!-------------------------------------------------------------------
	61	!! * SUBROUTINE ice_dyn_rhg_evp *
[8813]	62	!! EVP-C-grid
[8586]	63	!!
	64	!! ** purpose : determines sea ice drift from wind stress, ice-ocean
	65	!! stress and sea-surface slope. Ice-ice interaction is described by
	66	!! a non-linear elasto-viscous-plastic (EVP) law including shear
	67	!! strength and a bulk rheology (Hunke and Dukowicz, 2002).
	68	!!
	69	!! The points in the C-grid look like this, dear reader
	70	!!
	71	!! (ji,jj)
	72	!! \|
	73	!! \|
	74	!! (ji-1,jj) \| (ji,jj)
	75	!! ---------
	76	!! \| \|
	77	!! \| (ji,jj) \|------(ji,jj)
	78	!! \| \|
	79	!! ---------
	80	!! (ji-1,jj-1) (ji,jj-1)
	81	!!
	82	!! ** Inputs : - wind forcing (stress), oceanic currents
	83	!! ice total volume (vt_i) per unit area
	84	!! snow total volume (vt_s) per unit area
	85	!!
	86	!! ** Action : - compute u_ice, v_ice : the components of the
	87	!! sea-ice velocity vector
	88	!! - compute delta_i, shear_i, divu_i, which are inputs
	89	!! of the ice thickness distribution
	90	!!
	91	!! ** Steps : 0) compute mask at F point
	92	!! 1) Compute ice snow mass, ice strength
	93	!! 2) Compute wind, oceanic stresses, mass terms and
	94	!! coriolis terms of the momentum equation
	95	!! 3) Solve the momentum equation (iterative procedure)
	96	!! 4) Recompute delta, shear and divergence
	97	!! (which are inputs of the ITD) & store stress
	98	!! for the next time step
	99	!! 5) Diagnostics including charge ellipse
	100	!!
[8813]	101	!! ** Notes : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017)
	102	!! by setting up ln_aEVP=T (i.e. changing alpha and beta parameters).
	103	!! This is an upgraded version of mEVP from Bouillon et al. 2013
	104	!! (i.e. more stable and better convergence)
[8586]	105	!!
	106	!! References : Hunke and Dukowicz, JPO97
	107	!! Bouillon et al., Ocean Modelling 2009
	108	!! Bouillon et al., Ocean Modelling 2013
[8813]	109	!! Kimmritz et al., Ocean Modelling 2016 & 2017
[8586]	110	!!-------------------------------------------------------------------
[8813]	111	INTEGER , INTENT(in ) :: kt ! time step
[12377]	112	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
[8813]	113	REAL(wp), DIMENSION(:,:), INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i !
	114	REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear_i , pdivu_i , pdelta_i !
[8586]	115	!!
	116	INTEGER :: ji, jj ! dummy loop indices
	117	INTEGER :: jter ! local integers
[8813]	118	!
[12489]	119	REAL(wp) :: zrhoco ! rho0 * rn_cio
[9049]	120	REAL(wp) :: zdtevp, z1_dtevp ! time step for subcycling
	121	REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity
	122	REAL(wp) :: zalph1, z1_alph1, zalph2, z1_alph2 ! alpha coef from Bouillon 2009 or Kimmritz 2017
[10413]	123	REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV, zvU, zvV ! ice/snow mass and volume
[9049]	124	REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars
[11536]	125	REAL(wp) :: zTauO, zTauB, zRHS, zvel ! temporary scalars
[10413]	126	REAL(wp) :: zkt ! isotropic tensile strength for landfast ice
	127	REAL(wp) :: zvCr ! critical ice volume above which ice is landfast
[8813]	128	!
[9049]	129	REAL(wp) :: zresm ! Maximal error on ice velocity
	130	REAL(wp) :: zintb, zintn ! dummy argument
[8586]	131	REAL(wp) :: zfac_x, zfac_y
	132	REAL(wp) :: zshear, zdum1, zdum2
[8813]	133	!
[8586]	134	REAL(wp), DIMENSION(jpi,jpj) :: zp_delt ! P/delta at T points
[8813]	135	REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017
[8586]	136	!
[8813]	137	REAL(wp), DIMENSION(jpi,jpj) :: zdt_m ! (dt / ice-snow_mass) on T points
[11536]	138	REAL(wp), DIMENSION(jpi,jpj) :: zaU , zaV ! ice fraction on U/V points
[8813]	139	REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! (ice-snow_mass / dt) on U/V points
[8586]	140	REAL(wp), DIMENSION(jpi,jpj) :: zmf ! coriolis parameter at T points
	141	REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points
[8813]	142	!
[8586]	143	REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear
	144	REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components
[10425]	145	!!$ REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice, zresr ! check convergence
[10415]	146	REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope:
[9049]	147	! ! ocean surface (ssh_m) if ice is not embedded
[10415]	148	! ! ice bottom surface if ice is embedded
[11536]	149	REAL(wp), DIMENSION(jpi,jpj) :: zfU , zfV ! internal stresses
	150	REAL(wp), DIMENSION(jpi,jpj) :: zspgU, zspgV ! surface pressure gradient at U/V points
	151	REAL(wp), DIMENSION(jpi,jpj) :: zCorU, zCorV ! Coriolis stress array
	152	REAL(wp), DIMENSION(jpi,jpj) :: ztaux_ai, ztauy_ai ! ice-atm. stress at U-V points
	153	REAL(wp), DIMENSION(jpi,jpj) :: ztaux_oi, ztauy_oi ! ice-ocean stress at U-V points
	154	REAL(wp), DIMENSION(jpi,jpj) :: ztaux_bi, ztauy_bi ! ice-OceanBottom stress at U-V points (landfast)
	155	REAL(wp), DIMENSION(jpi,jpj) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast)
[8813]	156	!
[11536]	157	REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays
	158	REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence
[8586]	159	REAL(wp), DIMENSION(jpi,jpj) :: zfmask, zwf ! mask at F points for the ice
	160
	161	REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter
[10891]	162	REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small
	163	REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small
[8586]	164	!! --- diags
[11536]	165	REAL(wp), DIMENSION(jpi,jpj) :: zmsk00
[8586]	166	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig1, zsig2, zsig3
	167	!! --- SIMIP diags
	168	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
	169	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
	170	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
	171	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
	172	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xatrp ! X-component of area transport (m2/s)
	173	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_yatrp ! Y-component of area transport (m2/s)
	174	!!-------------------------------------------------------------------
	175
	176	IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology'
	177	!
[8813]	178	!!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization....
[8586]	179	!------------------------------------------------------------------------------!
	180	! 0) mask at F points for the ice
	181	!------------------------------------------------------------------------------!
	182	! ocean/land mask
[13295]	183	DO_2D( 1, 0, 1, 0 )
[12377]	184	zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
	185	END_2D
[13570]	186	#if defined key_mpi3
[13561]	187	CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', zfmask, 'F', 1._wp)
[13570]	188	#else
	189	CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp)
	190	#endif
[8586]	191
	192	! Lateral boundary conditions on velocity (modify zfmask)
	193	zwf(:,:) = zfmask(:,:)
[13295]	194	DO_2D( 0, 0, 0, 0 )
[12377]	195	IF( zfmask(ji,jj) == 0._wp ) THEN
	196	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) ) )
	197	ENDIF
	198	END_2D
[8586]	199	DO jj = 2, jpjm1
	200	IF( zfmask(1,jj) == 0._wp ) THEN
	201	zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
	202	ENDIF
	203	IF( zfmask(jpi,jj) == 0._wp ) THEN
	204	zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
	205	ENDIF
	206	END DO
	207	DO ji = 2, jpim1
	208	IF( zfmask(ji,1) == 0._wp ) THEN
	209	zfmask(ji,1 ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
	210	ENDIF
	211	IF( zfmask(ji,jpj) == 0._wp ) THEN
	212	zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
	213	ENDIF
	214	END DO
[13570]	215	#if defined key_mpi3
[13561]	216	CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', zfmask, 'F', 1._wp )
[13570]	217	#else
	218	CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp )
	219	#endif
[8586]	220
	221	!------------------------------------------------------------------------------!
	222	! 1) define some variables and initialize arrays
	223	!------------------------------------------------------------------------------!
[12489]	224	zrhoco = rho0 * rn_cio
[8586]	225
	226	! ecc2: square of yield ellipse eccenticrity
	227	ecc2 = rn_ecc * rn_ecc
	228	z1_ecc2 = 1._wp / ecc2
	229
	230	! Time step for subcycling
[12489]	231	zdtevp = rDt_ice / REAL( nn_nevp )
[8586]	232	z1_dtevp = 1._wp / zdtevp
	233
	234	! alpha parameters (Bouillon 2009)
[8813]	235	IF( .NOT. ln_aEVP ) THEN
[12489]	236	zalph1 = ( 2._wp * rn_relast * rDt_ice ) * z1_dtevp
[8813]	237	zalph2 = zalph1 * z1_ecc2
[8586]	238
[8813]	239	z1_alph1 = 1._wp / ( zalph1 + 1._wp )
	240	z1_alph2 = 1._wp / ( zalph2 + 1._wp )
	241	ENDIF
	242
[8586]	243	! Initialise stress tensor
	244	zs1 (:,:) = pstress1_i (:,:)
	245	zs2 (:,:) = pstress2_i (:,:)
	246	zs12(:,:) = pstress12_i(:,:)
	247
	248	! Ice strength
	249	CALL ice_strength
	250
[10413]	251	! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
[11536]	252	IF( ln_landfast_L16 ) THEN ; zkt = rn_tensile
	253	ELSE ; zkt = 0._wp
[10413]	254	ENDIF
[8586]	255	!
	256	!------------------------------------------------------------------------------!
	257	! 2) Wind / ocean stress, mass terms, coriolis terms
	258	!------------------------------------------------------------------------------!
[10415]	259	! sea surface height
	260	! embedded sea ice: compute representative ice top surface
	261	! non-embedded sea ice: use ocean surface for slope calculation
	262	zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
[8586]	263
[13295]	264	DO_2D( 0, 0, 0, 0 )
[8586]	265
[12377]	266	! ice fraction at U-V points
	267	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)
	268	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)
[8586]	269
[12377]	270	! Ice/snow mass at U-V points
	271	zm1 = ( rhos * vt_s(ji ,jj ) + rhoi * vt_i(ji ,jj ) )
	272	zm2 = ( rhos * vt_s(ji+1,jj ) + rhoi * vt_i(ji+1,jj ) )
	273	zm3 = ( rhos * vt_s(ji ,jj+1) + rhoi * vt_i(ji ,jj+1) )
	274	zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
	275	zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
[8586]	276
[12377]	277	! Ocean currents at U-V points
	278	v_oceU(ji,jj) = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1)
	279	u_oceV(ji,jj) = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1)
[8586]	280
[12377]	281	! Coriolis at T points (m*f)
	282	zmf(ji,jj) = zm1 * ff_t(ji,jj)
[8586]	283
[12377]	284	! dt/m at T points (for alpha and beta coefficients)
	285	zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin )
	286
	287	! m/dt
	288	zmU_t(ji,jj) = zmassU * z1_dtevp
	289	zmV_t(ji,jj) = zmassV * z1_dtevp
	290
	291	! Drag ice-atm.
	292	ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
	293	ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
[8586]	294
[12377]	295	! Surface pressure gradient (- mgGRAD(ssh)) at U-V points
	296	zspgU(ji,jj) = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
	297	zspgV(ji,jj) = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
[8586]	298
[12377]	299	! masks
	300	zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice
	301	zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice
[8586]	302
[12377]	303	! switches
	304	IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN ; zmsk01x(ji,jj) = 0._wp
	305	ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF
	306	IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN ; zmsk01y(ji,jj) = 0._wp
	307	ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF
[8586]	308
[12377]	309	END_2D
[13570]	310	#if defined key_mpi3
[13561]	311	CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', zmf, 'T', 1.0_wp, zdt_m, 'T', 1.0_wp )
[13570]	312	#else
	313	CALL lbc_lnk_multi( 'icedyn_rhg_evp', zmf, 'T', 1.0_wp, zdt_m, 'T', 1.0_wp )
	314	#endif
[8586]	315	!
[10413]	316	! !== Landfast ice parameterization ==!
	317	!
	318	IF( ln_landfast_L16 ) THEN !-- Lemieux 2016
[13295]	319	DO_2D( 0, 0, 0, 0 )
[12377]	320	! ice thickness at U-V points
	321	zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
	322	zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
	323	! ice-bottom stress at U points
	324	zvCr = zaU(ji,jj) * rn_depfra * hu(ji,jj,Kmm)
	325	ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
	326	! ice-bottom stress at V points
	327	zvCr = zaV(ji,jj) * rn_depfra * hv(ji,jj,Kmm)
	328	ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
	329	! ice_bottom stress at T points
	330	zvCr = at_i(ji,jj) * rn_depfra * ht(ji,jj)
	331	tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
	332	END_2D
[13570]	333	#if defined key_mpi3
[13561]	334	CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1.0_wp )
[13570]	335	#else
	336	CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1.0_wp )
	337	#endif
[10413]	338	!
[11536]	339	ELSE !-- no landfast
[13295]	340	DO_2D( 0, 0, 0, 0 )
[12377]	341	ztaux_base(ji,jj) = 0._wp
	342	ztauy_base(ji,jj) = 0._wp
	343	END_2D
[10413]	344	ENDIF
	345
[8586]	346	!------------------------------------------------------------------------------!
	347	! 3) Solution of the momentum equation, iterative procedure
	348	!------------------------------------------------------------------------------!
	349	!
[11536]	350	! ! ==================== !
[8586]	351	DO jter = 1 , nn_nevp ! loop over jter !
[11536]	352	! ! ==================== !
[10425]	353	l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
	354	!
[12377]	355	!!$ IF(sn_cfctl%l_prtctl) THEN ! Convergence test
[10425]	356	!!$ DO jj = 1, jpjm1
	357	!!$ zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
	358	!!$ zv_ice(:,jj) = v_ice(:,jj)
	359	!!$ END DO
	360	!!$ ENDIF
[8586]	361
	362	! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
[13295]	363	DO_2D( 1, 0, 1, 0 )
[8586]	364
[12377]	365	! shear at F points
	366	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) &
	367	& + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
	368	& ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
[8586]	369
[12377]	370	END_2D
[13226]	371	CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1.0_wp )
[8586]	372
[13295]	373	DO_2D( 0, 1, 0, 1 )
[8586]	374
[12377]	375	! shear**2 at T points (doc eq. A16)
	376	zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) &
	377	& + 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) &
	378	& ) * 0.25_wp * r1_e1e2t(ji,jj)
	379
	380	! divergence at T points
	381	zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
	382	& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
	383	& ) * r1_e1e2t(ji,jj)
	384	zdiv2 = zdiv * zdiv
	385
	386	! tension at T points
	387	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) &
	388	& - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
	389	& ) * r1_e1e2t(ji,jj)
	390	zdt2 = zdt * zdt
	391
	392	! delta at T points
	393	zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )
[8586]	394
[12377]	395	! P/delta at T points
	396	zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
[8813]	397
[12377]	398	! alpha & beta for aEVP
	399	! gamma = 0.5P/(delta+creepl) (cpi)2/Area dt/m
	400	! alpha = beta = sqrt(4*gamma)
	401	IF( ln_aEVP ) THEN
	402	zalph1 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
	403	z1_alph1 = 1._wp / ( zalph1 + 1._wp )
	404	zalph2 = zalph1
	405	z1_alph2 = z1_alph1
	406	ENDIF
	407
	408	! stress at T points (zkt/=0 if landfast)
	409	zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1
	410	zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
	411
	412	END_2D
[13570]	413	#if defined key_mpi3
[13561]	414	CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', zp_delt, 'T', 1.0_wp )
[13570]	415	#else
	416	CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1.0_wp )
	417	#endif
[13295]	418	DO_2D( 1, 0, 1, 0 )
[8586]	419
[12377]	420	! alpha & beta for aEVP
	421	IF( ln_aEVP ) THEN
	422	zalph2 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
	423	z1_alph2 = 1._wp / ( zalph2 + 1._wp )
	424	zbeta(ji,jj) = zalph2
	425	ENDIF
	426
	427	! P/delta at F points
	428	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) )
	429
	430	! stress at F points (zkt/=0 if landfast)
	431	zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
[8586]	432
[12377]	433	END_2D
[8586]	434
	435	! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
[13295]	436	DO_2D( 0, 0, 0, 0 )
[12377]	437	! !--- U points
	438	zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) &
	439	& + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) &
	440	& ) * r1_e2u(ji,jj) &
	441	& + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) &
	442	& ) * 2._wp * r1_e1u(ji,jj) &
	443	& ) * r1_e1e2u(ji,jj)
	444	!
	445	! !--- V points
	446	zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) &
	447	& - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) &
	448	& ) * r1_e1v(ji,jj) &
	449	& + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) &
	450	& ) * 2._wp * r1_e2v(ji,jj) &
	451	& ) * r1_e1e2v(ji,jj)
	452	!
	453	! !--- ice currents at U-V point
	454	v_iceU(ji,jj) = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1)
	455	u_iceV(ji,jj) = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1)
	456	!
	457	END_2D
[8586]	458	!
	459	! --- Computation of ice velocity --- !
[8813]	460	! Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
[8586]	461	! Bouillon et al. 2009 (eq 34-35) => stable
	462	IF( MOD(jter,2) == 0 ) THEN ! even iterations
	463	!
[13295]	464	DO_2D( 0, 0, 0, 0 )
[12377]	465	! !--- tau_io/(v_oce - v_ice)
	466	zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) &
	467	& + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
	468	! !--- Ocean-to-Ice stress
	469	ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
	470	!
	471	! !--- tau_bottom/v_ice
	472	zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
	473	zTauB = ztauy_base(ji,jj) / zvel
	474	! !--- OceanBottom-to-Ice stress
	475	ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
	476	!
	477	! !--- Coriolis at V-points (energy conserving formulation)
	478	zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * &
	479	& ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) &
	480	& + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
	481	!
	482	! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
	483	zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
	484	!
	485	! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS)
	486	! 1 = sliding friction : TauB < RHS
	487	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
	488	!
	489	IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
	490	v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity
	491	& + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	492	& / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	493	& + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
	494	& ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
	495	& ) * zmsk00y(ji,jj)
	496	ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
	497	v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity
	498	& + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	499	& / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	500	& + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
	501	& ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
	502	& ) * zmsk00y(ji,jj)
	503	ENDIF
	504	END_2D
[13570]	505	#if defined key_mpi3
[13561]	506	CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp )
[13570]	507	#else
	508	CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp )
	509	#endif
[8586]	510	!
	511	#if defined key_agrif
[9610]	512	!! CALL agrif_interp_ice( 'V', jter, nn_nevp )
	513	CALL agrif_interp_ice( 'V' )
[8586]	514	#endif
[11536]	515	IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
[8586]	516	!
[13295]	517	DO_2D( 0, 0, 0, 0 )
[12377]	518	! !--- tau_io/(u_oce - u_ice)
	519	zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) &
	520	& + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
	521	! !--- Ocean-to-Ice stress
	522	ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
	523	!
	524	! !--- tau_bottom/u_ice
	525	zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
	526	zTauB = ztaux_base(ji,jj) / zvel
	527	! !--- OceanBottom-to-Ice stress
	528	ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
	529	!
	530	! !--- Coriolis at U-points (energy conserving formulation)
	531	zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * &
	532	& ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) &
	533	& + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
	534	!
	535	! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
	536	zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
	537	!
	538	! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS)
	539	! 1 = sliding friction : TauB < RHS
	540	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
	541	!
	542	IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
	543	u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity
	544	& + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	545	& / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	546	& + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
	547	& ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
	548	& ) * zmsk00x(ji,jj)
	549	ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
	550	u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity
	551	& + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	552	& / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	553	& + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
	554	& ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
	555	& ) * zmsk00x(ji,jj)
	556	ENDIF
	557	END_2D
[13570]	558	#if defined key_mpi3
[13561]	559	CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp )
[13570]	560	#else
	561	CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp )
	562	#endif
[8586]	563	!
	564	#if defined key_agrif
[9610]	565	!! CALL agrif_interp_ice( 'U', jter, nn_nevp )
	566	CALL agrif_interp_ice( 'U' )
[8586]	567	#endif
[11536]	568	IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
[8586]	569	!
	570	ELSE ! odd iterations
	571	!
[13295]	572	DO_2D( 0, 0, 0, 0 )
[12377]	573	! !--- tau_io/(u_oce - u_ice)
	574	zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) &
	575	& + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
	576	! !--- Ocean-to-Ice stress
	577	ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
	578	!
	579	! !--- tau_bottom/u_ice
	580	zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
	581	zTauB = ztaux_base(ji,jj) / zvel
	582	! !--- OceanBottom-to-Ice stress
	583	ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
	584	!
	585	! !--- Coriolis at U-points (energy conserving formulation)
	586	zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * &
	587	& ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) &
	588	& + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
	589	!
	590	! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
	591	zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
	592	!
	593	! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS)
	594	! 1 = sliding friction : TauB < RHS
	595	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
	596	!
	597	IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
	598	u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity
	599	& + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	600	& / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	601	& + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
	602	& ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
	603	& ) * zmsk00x(ji,jj)
	604	ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
	605	u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity
	606	& + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	607	& / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	608	& + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
	609	& ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
	610	& ) * zmsk00x(ji,jj)
	611	ENDIF
	612	END_2D
[13570]	613	#if defined key_mpi3
[13561]	614	CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp )
[13570]	615	#else
	616	CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp )
	617	#endif
[8586]	618	!
	619	#if defined key_agrif
[9610]	620	!! CALL agrif_interp_ice( 'U', jter, nn_nevp )
	621	CALL agrif_interp_ice( 'U' )
[8586]	622	#endif
[11536]	623	IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
[8586]	624	!
[13295]	625	DO_2D( 0, 0, 0, 0 )
[12377]	626	! !--- tau_io/(v_oce - v_ice)
	627	zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) &
	628	& + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
	629	! !--- Ocean-to-Ice stress
	630	ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
	631	!
	632	! !--- tau_bottom/v_ice
	633	zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
	634	zTauB = ztauy_base(ji,jj) / zvel
	635	! !--- OceanBottom-to-Ice stress
	636	ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
	637	!
	638	! !--- Coriolis at v-points (energy conserving formulation)
	639	zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * &
	640	& ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) &
	641	& + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
	642	!
	643	! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
	644	zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
	645	!
	646	! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS)
	647	! 1 = sliding friction : TauB < RHS
	648	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
	649	!
	650	IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
	651	v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity
	652	& + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	653	& / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	654	& + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
	655	& ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
	656	& ) * zmsk00y(ji,jj)
	657	ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
	658	v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity
	659	& + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	660	& / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	661	& + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
	662	& ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
	663	& ) * zmsk00y(ji,jj)
	664	ENDIF
	665	END_2D
[13570]	666	#if defined key_mpi3
[13561]	667	CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp )
[13570]	668	#else
	669	CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp )
	670	#endif
[8586]	671	!
	672	#if defined key_agrif
[9610]	673	!! CALL agrif_interp_ice( 'V', jter, nn_nevp )
	674	CALL agrif_interp_ice( 'V' )
[8586]	675	#endif
[11536]	676	IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
[8586]	677	!
	678	ENDIF
[10425]	679
[12377]	680	!!$ IF(sn_cfctl%l_prtctl) THEN ! Convergence test
[10425]	681	!!$ DO jj = 2 , jpjm1
	682	!!$ zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
	683	!!$ END DO
	684	!!$ zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
	685	!!$ CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain
	686	!!$ ENDIF
[8586]	687	!
	688	! ! ==================== !
	689	END DO ! end loop over jter !
	690	! ! ==================== !
	691	!
	692	!------------------------------------------------------------------------------!
	693	! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
	694	!------------------------------------------------------------------------------!
[13295]	695	DO_2D( 1, 0, 1, 0 )
[8586]	696
[12377]	697	! shear at F points
	698	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) &
	699	& + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
	700	& ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
[8586]	701
[12377]	702	END_2D
[8586]	703
[13295]	704	DO_2D( 0, 0, 0, 0 )
[12377]	705
	706	! tension**2 at T points
	707	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) &
	708	& - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
	709	& ) * r1_e1e2t(ji,jj)
	710	zdt2 = zdt * zdt
	711
	712	! shear**2 at T points (doc eq. A16)
	713	zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) &
	714	& + 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) &
	715	& ) * 0.25_wp * r1_e1e2t(ji,jj)
	716
	717	! shear at T points
	718	pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
[8586]	719
[12377]	720	! divergence at T points
	721	pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
	722	& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
	723	& ) * r1_e1e2t(ji,jj)
	724
	725	! delta at T points
	726	zdelta = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )
	727	rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
	728	pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
[8586]	729
[12377]	730	END_2D
[13570]	731	#if defined key_mpi3
[13561]	732	CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1.0_wp, pdivu_i, 'T', 1.0_wp, pdelta_i, 'T', 1.0_wp )
[13570]	733	#else
	734	CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1.0_wp, pdivu_i, 'T', 1.0_wp, pdelta_i, 'T', 1.0_wp )
	735	#endif
[8586]	736
	737	! --- Store the stress tensor for the next time step --- !
[13570]	738	#if defined key_mpi3
[13561]	739	CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp, zs12, 'F', 1.0_wp )
[13570]	740	#else
	741	CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp, zs12, 'F', 1.0_wp )
	742	#endif
[8586]	743	pstress1_i (:,:) = zs1 (:,:)
	744	pstress2_i (:,:) = zs2 (:,:)
	745	pstress12_i(:,:) = zs12(:,:)
	746	!
	747
	748	!------------------------------------------------------------------------------!
	749	! 5) diagnostics
	750	!------------------------------------------------------------------------------!
[13295]	751	DO_2D( 1, 1, 1, 1 )
[12377]	752	zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
	753	END_2D
[8586]	754
[11536]	755	! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
	756	IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
	757	& iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
	758	!
[13570]	759	#if defined key_mpi3
[13561]	760	CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1.0_wp, ztauy_oi, 'V', -1.0_wp, ztaux_ai, 'U', -1.0_wp, ztauy_ai, 'V', -1.0_wp, &
[13226]	761	& ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp )
[13570]	762	#else
	763	CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1.0_wp, ztauy_oi, 'V', -1.0_wp, ztaux_ai, 'U', -1.0_wp, ztauy_ai, 'V', -1.0_wp, &
	764	& ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp )
	765	#endif
[11536]	766	!
	767	CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
	768	CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
	769	CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
	770	CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
	771	CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
	772	CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
	773	ENDIF
	774
[8586]	775	! --- divergence, shear and strength --- !
[11536]	776	IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence
	777	IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear
	778	IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength
[8586]	779
[11536]	780	! --- stress tensor --- !
	781	IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN
[8586]	782	!
	783	ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
	784	!
[13295]	785	DO_2D( 0, 0, 0, 0 )
[12377]	786	zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point
	787	& zmsk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) &
	788	& / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) )
[8586]	789
[12377]	790	zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress
[8586]	791
[12377]	792	zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) )
[8586]	793
	794	!! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
	795	!! zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
	796	!! 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
	797	!! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
[12377]	798	zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015
	799	zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress
	800	zsig3(ji,jj) = zdum2*2 ( ( pstress1_i(ji,jj) + strength(ji,jj) )*2 + ( rn_ecc zshear )**2 )
	801	END_2D
[13570]	802	#if defined key_mpi3
[13561]	803	CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', zsig1, 'T', 1.0_wp, zsig2, 'T', 1.0_wp, zsig3, 'T', 1.0_wp )
[13570]	804	#else
	805	CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1.0_wp, zsig2, 'T', 1.0_wp, zsig3, 'T', 1.0_wp )
	806	#endif
[8586]	807	!
[11536]	808	CALL iom_put( 'isig1' , zsig1 )
	809	CALL iom_put( 'isig2' , zsig2 )
	810	CALL iom_put( 'isig3' , zsig3 )
[8586]	811	!
[11536]	812	! Stress tensor invariants (normal and shear stress N/m)
	813	IF( iom_use('normstr') ) CALL iom_put( 'normstr' , ( zs1(:,:) + zs2(:,:) ) * zmsk00(:,:) ) ! Normal stress
	814	IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )*2 + 4zs12(:,:)*2 ) zmsk00(:,:) ) ! Shear stress
	815
[8586]	816	DEALLOCATE( zsig1 , zsig2 , zsig3 )
	817	ENDIF
	818
	819	! --- SIMIP --- !
[11536]	820	IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
	821	& iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
	822	!
[13570]	823	#if defined key_mpi3
[13561]	824	CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', zspgU, 'U', -1.0_wp, zspgV, 'V', -1.0_wp, &
[13226]	825	& zCorU, 'U', -1.0_wp, zCorV, 'V', -1.0_wp, zfU, 'U', -1.0_wp, zfV, 'V', -1.0_wp )
[13570]	826	#else
	827	CALL lbc_lnk_multi( 'icedyn_rhg_evp', zspgU, 'U', -1.0_wp, zspgV, 'V', -1.0_wp, &
	828	& zCorU, 'U', -1.0_wp, zCorV, 'V', -1.0_wp, zfU, 'U', -1.0_wp, zfV, 'V', -1.0_wp )
	829	#endif
[8586]	830
[11536]	831	CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x)
	832	CALL iom_put( 'dssh_dy' , zspgV * zmsk00 ) ! Sea-surface tilt term in force balance (y)
	833	CALL iom_put( 'corstrx' , zCorU * zmsk00 ) ! Coriolis force term in force balance (x)
	834	CALL iom_put( 'corstry' , zCorV * zmsk00 ) ! Coriolis force term in force balance (y)
	835	CALL iom_put( 'intstrx' , zfU * zmsk00 ) ! Internal force term in force balance (x)
	836	CALL iom_put( 'intstry' , zfV * zmsk00 ) ! Internal force term in force balance (y)
	837	ENDIF
	838
	839	IF( iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
	840	& iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
	841	!
	842	ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
	843	& zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
	844	!
[13295]	845	DO_2D( 0, 0, 0, 0 )
[12377]	846	! 2D ice mass, snow mass, area transport arrays (X, Y)
	847	zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
	848	zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)
[11536]	849
[12377]	850	zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
	851	zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) ! '' Y- ''
[11536]	852
[12377]	853	zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
	854	zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) ! '' Y- ''
[11536]	855
[12377]	856	zdiag_xatrp(ji,jj) = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) ) ! area transport, X-component
	857	zdiag_yatrp(ji,jj) = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) ) ! '' Y- ''
[11536]	858
[12377]	859	END_2D
[8586]	860
[13570]	861	#if defined key_mpi3
[13561]	862	CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1.0_wp, zdiag_ymtrp_ice, 'V', -1.0_wp, &
[13226]	863	& zdiag_xmtrp_snw, 'U', -1.0_wp, zdiag_ymtrp_snw, 'V', -1.0_wp, &
	864	& zdiag_xatrp , 'U', -1.0_wp, zdiag_yatrp , 'V', -1.0_wp )
[13570]	865	#else
	866	CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1.0_wp, zdiag_ymtrp_ice, 'V', -1.0_wp, &
	867	& zdiag_xmtrp_snw, 'U', -1.0_wp, zdiag_ymtrp_snw, 'V', -1.0_wp, &
	868	& zdiag_xatrp , 'U', -1.0_wp, zdiag_yatrp , 'V', -1.0_wp )
	869	#endif
[8586]	870
[11536]	871	CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice ) ! X-component of sea-ice mass transport (kg/s)
	872	CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice ) ! Y-component of sea-ice mass transport
	873	CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw ) ! X-component of snow mass transport (kg/s)
	874	CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw ) ! Y-component of snow mass transport
	875	CALL iom_put( 'xatrp' , zdiag_xatrp ) ! X-component of ice area transport
	876	CALL iom_put( 'yatrp' , zdiag_yatrp ) ! Y-component of ice area transport
	877
	878	DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , &
	879	& zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp )
	880
[8586]	881	ENDIF
	882	!
	883	END SUBROUTINE ice_dyn_rhg_evp
	884
[8813]	885
[8586]	886	SUBROUTINE rhg_evp_rst( cdrw, kt )
	887	!!---------------------------------------------------------------------
	888	!! * ROUTINE rhg_evp_rst *
	889	!!
	890	!! ** Purpose : Read or write RHG file in restart file
	891	!!
	892	!! ** Method : use of IOM library
	893	!!----------------------------------------------------------------------
	894	CHARACTER(len=*) , INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	895	INTEGER, OPTIONAL, INTENT(in) :: kt ! ice time-step
	896	!
	897	INTEGER :: iter ! local integer
	898	INTEGER :: id1, id2, id3 ! local integers
	899	!!----------------------------------------------------------------------
	900	!
	901	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialize
	902	! ! ---------------
	903	IF( ln_rstart ) THEN !* Read the restart file
	904	!
	905	id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
	906	id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
	907	id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
	908	!
	909	IF( MIN( id1, id2, id3 ) > 0 ) THEN ! fields exist
[13286]	910	CALL iom_get( numrir, jpdom_auto, 'stress1_i' , stress1_i , cd_type = 'T' )
	911	CALL iom_get( numrir, jpdom_auto, 'stress2_i' , stress2_i , cd_type = 'T' )
	912	CALL iom_get( numrir, jpdom_auto, 'stress12_i', stress12_i, cd_type = 'F' )
[8586]	913	ELSE ! start rheology from rest
[9169]	914	IF(lwp) WRITE(numout,*)
	915	IF(lwp) WRITE(numout,*) ' ==>>> previous run without rheology, set stresses to 0'
[8586]	916	stress1_i (:,:) = 0._wp
	917	stress2_i (:,:) = 0._wp
	918	stress12_i(:,:) = 0._wp
	919	ENDIF
	920	ELSE !* Start from rest
[9169]	921	IF(lwp) WRITE(numout,*)
	922	IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set stresses to 0'
[8586]	923	stress1_i (:,:) = 0._wp
	924	stress2_i (:,:) = 0._wp
	925	stress12_i(:,:) = 0._wp
	926	ENDIF
	927	!
	928	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
	929	! ! -------------------
	930	IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
	931	iter = kt + nn_fsbc - 1 ! ice restarts are written at kt == nitrst - nn_fsbc + 1
	932	!
	933	CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i )
	934	CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i )
	935	CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
	936	!
	937	ENDIF
	938	!
	939	END SUBROUTINE rhg_evp_rst
	940
	941	#else
	942	!!----------------------------------------------------------------------
[9570]	943	!! Default option Empty module NO SI3 sea-ice model
[8586]	944	!!----------------------------------------------------------------------
	945	#endif
	946
	947	!!==============================================================================
	948	END MODULE icedyn_rhg_evp

Note: See TracBrowser for help on using the repository browser.

Download in other formats: