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/2021/ticket2680_C1D_PAPA/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2021/ticket2680_C1D_PAPA/src/ICE/icedyn_rhg_evp.F90 @ 15020

Last change on this file since 15020 was 15020, checked in by gsamson, 3 years ago
merge trunk into branch (#2680)
Property svn:keywords set to `Id`
File size: 60.9 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
[14072]	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
[14072]	10	!! 3.4 ! 2011-01 (A. Porter) dynamical allocation
[9604]	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
[14072]	30	USE bdy_oce , ONLY : ln_bdy
	31	USE bdyice
[8813]	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
[13472]	43	USE netcdf ! NetCDF library for convergence test
[8586]	44	IMPLICIT NONE
	45	PRIVATE
	46
	47	PUBLIC ice_dyn_rhg_evp ! called by icedyn_rhg.F90
	48	PUBLIC rhg_evp_rst ! called by icedyn_rhg.F90
	49
[13472]	50	!! for convergence tests
	51	INTEGER :: ncvgid ! netcdf file id
	52	INTEGER :: nvarid ! netcdf variable id
[15020]	53	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fimask ! mask at F points for the ice
	54
	55	!! * Substitutions
	56	# include "do_loop_substitute.h90"
[8586]	57	!!----------------------------------------------------------------------
[9598]	58	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	59	!! $Id$
[10068]	60	!! Software governed by the CeCILL license (see ./LICENSE)
[8586]	61	!!----------------------------------------------------------------------
	62	CONTAINS
	63
[12377]	64	SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )
[8586]	65	!!-------------------------------------------------------------------
	66	!! * SUBROUTINE ice_dyn_rhg_evp *
[8813]	67	!! EVP-C-grid
[8586]	68	!!
	69	!! ** purpose : determines sea ice drift from wind stress, ice-ocean
[14072]	70	!! stress and sea-surface slope. Ice-ice interaction is described by
	71	!! a non-linear elasto-viscous-plastic (EVP) law including shear
	72	!! strength and a bulk rheology (Hunke and Dukowicz, 2002).
[8586]	73	!!
	74	!! The points in the C-grid look like this, dear reader
	75	!!
	76	!! (ji,jj)
	77	!! \|
	78	!! \|
	79	!! (ji-1,jj) \| (ji,jj)
[14072]	80	!! ---------
[8586]	81	!! \| \|
	82	!! \| (ji,jj) \|------(ji,jj)
	83	!! \| \|
[14072]	84	!! ---------
[8586]	85	!! (ji-1,jj-1) (ji,jj-1)
	86	!!
	87	!! ** Inputs : - wind forcing (stress), oceanic currents
	88	!! ice total volume (vt_i) per unit area
	89	!! snow total volume (vt_s) per unit area
	90	!!
[14072]	91	!! ** Action : - compute u_ice, v_ice : the components of the
[8586]	92	!! sea-ice velocity vector
	93	!! - compute delta_i, shear_i, divu_i, which are inputs
	94	!! of the ice thickness distribution
	95	!!
	96	!! ** Steps : 0) compute mask at F point
[14072]	97	!! 1) Compute ice snow mass, ice strength
[8586]	98	!! 2) Compute wind, oceanic stresses, mass terms and
	99	!! coriolis terms of the momentum equation
	100	!! 3) Solve the momentum equation (iterative procedure)
	101	!! 4) Recompute delta, shear and divergence
	102	!! (which are inputs of the ITD) & store stress
	103	!! for the next time step
	104	!! 5) Diagnostics including charge ellipse
	105	!!
[8813]	106	!! ** Notes : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017)
	107	!! by setting up ln_aEVP=T (i.e. changing alpha and beta parameters).
	108	!! This is an upgraded version of mEVP from Bouillon et al. 2013
	109	!! (i.e. more stable and better convergence)
[8586]	110	!!
	111	!! References : Hunke and Dukowicz, JPO97
	112	!! Bouillon et al., Ocean Modelling 2009
	113	!! Bouillon et al., Ocean Modelling 2013
[8813]	114	!! Kimmritz et al., Ocean Modelling 2016 & 2017
[8586]	115	!!-------------------------------------------------------------------
[8813]	116	INTEGER , INTENT(in ) :: kt ! time step
[12377]	117	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
[8813]	118	REAL(wp), DIMENSION(:,:), INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i !
	119	REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear_i , pdivu_i , pdelta_i !
[8586]	120	!!
	121	INTEGER :: ji, jj ! dummy loop indices
	122	INTEGER :: jter ! local integers
[8813]	123	!
[12489]	124	REAL(wp) :: zrhoco ! rho0 * rn_cio
[9049]	125	REAL(wp) :: zdtevp, z1_dtevp ! time step for subcycling
	126	REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity
	127	REAL(wp) :: zalph1, z1_alph1, zalph2, z1_alph2 ! alpha coef from Bouillon 2009 or Kimmritz 2017
[13472]	128	REAl(wp) :: zbetau, zbetav
[10413]	129	REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV, zvU, zvV ! ice/snow mass and volume
[13550]	130	REAL(wp) :: zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars
[11536]	131	REAL(wp) :: zTauO, zTauB, zRHS, zvel ! temporary scalars
[10413]	132	REAL(wp) :: zkt ! isotropic tensile strength for landfast ice
	133	REAL(wp) :: zvCr ! critical ice volume above which ice is landfast
[8813]	134	!
[9049]	135	REAL(wp) :: zintb, zintn ! dummy argument
[8586]	136	REAL(wp) :: zfac_x, zfac_y
	137	REAL(wp) :: zshear, zdum1, zdum2
[8813]	138	!
[13550]	139	REAL(wp), DIMENSION(jpi,jpj) :: zdelta, zp_delt ! delta and P/delta at T points
[8813]	140	REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017
[8586]	141	!
[8813]	142	REAL(wp), DIMENSION(jpi,jpj) :: zdt_m ! (dt / ice-snow_mass) on T points
[11536]	143	REAL(wp), DIMENSION(jpi,jpj) :: zaU , zaV ! ice fraction on U/V points
[8813]	144	REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! (ice-snow_mass / dt) on U/V points
[8586]	145	REAL(wp), DIMENSION(jpi,jpj) :: zmf ! coriolis parameter at T points
[13550]	146	REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points
[8813]	147	!
[8586]	148	REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear
[13612]	149	REAL(wp), DIMENSION(jpi,jpj) :: zten_i ! tension
[8586]	150	REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components
[10415]	151	REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope:
[9049]	152	! ! ocean surface (ssh_m) if ice is not embedded
[14072]	153	! ! ice bottom surface if ice is embedded
[11536]	154	REAL(wp), DIMENSION(jpi,jpj) :: zfU , zfV ! internal stresses
	155	REAL(wp), DIMENSION(jpi,jpj) :: zspgU, zspgV ! surface pressure gradient at U/V points
	156	REAL(wp), DIMENSION(jpi,jpj) :: zCorU, zCorV ! Coriolis stress array
	157	REAL(wp), DIMENSION(jpi,jpj) :: ztaux_ai, ztauy_ai ! ice-atm. stress at U-V points
	158	REAL(wp), DIMENSION(jpi,jpj) :: ztaux_oi, ztauy_oi ! ice-ocean stress at U-V points
	159	REAL(wp), DIMENSION(jpi,jpj) :: ztaux_bi, ztauy_bi ! ice-OceanBottom stress at U-V points (landfast)
	160	REAL(wp), DIMENSION(jpi,jpj) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast)
[8813]	161	!
[15020]	162	REAL(wp), DIMENSION(jpi,jpj) :: zmsk00, zmsk15
[11536]	163	REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays
	164	REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence
[8586]	165
	166	REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter
[10891]	167	REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small
	168	REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small
[13472]	169	!! --- check convergence
	170	REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice
[8586]	171	!! --- diags
[13601]	172	REAL(wp) :: zsig1, zsig2, zsig12, zfac, z1_strength
[14072]	173	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p
[8586]	174	!! --- SIMIP diags
	175	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
	176	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
	177	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
	178	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
	179	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xatrp ! X-component of area transport (m2/s)
[14072]	180	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_yatrp ! Y-component of area transport (m2/s)
[8586]	181	!!-------------------------------------------------------------------
	182
	183	IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology'
	184	!
[13472]	185	! for diagnostics and convergence tests
	186	DO_2D( 1, 1, 1, 1 )
	187	zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice
	188	END_2D
[15020]	189	IF( nn_rhg_chkcvg > 0 ) THEN
	190	DO_2D( 1, 1, 1, 1 )
	191	zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
	192	END_2D
	193	ENDIF
[13472]	194	!
[8586]	195	!------------------------------------------------------------------------------!
	196	! 0) mask at F points for the ice
	197	!------------------------------------------------------------------------------!
[15020]	198	IF( kt == nit000 ) THEN
	199	! ocean/land mask
	200	ALLOCATE( fimask(jpi,jpj) )
	201	IF( rn_ishlat == 0._wp ) THEN
	202	DO_2D( 0, 0, 0, 0 )
	203	fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
	204	END_2D
	205	ELSE
	206	DO_2D( 0, 0, 0, 0 )
	207	fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
	208	! Lateral boundary conditions on velocity (modify fimask)
	209	IF( fimask(ji,jj) == 0._wp ) THEN
	210	fimask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), &
	211	& vmask(ji,jj,1), vmask(ji+1,jj,1) ) )
	212	ENDIF
	213	END_2D
[12377]	214	ENDIF
[15020]	215	CALL lbc_lnk( 'icedyn_rhg_evp', fimask, 'F', 1._wp )
	216	ENDIF
[8586]	217	!------------------------------------------------------------------------------!
	218	! 1) define some variables and initialize arrays
	219	!------------------------------------------------------------------------------!
[14072]	220	zrhoco = rho0 * rn_cio
[8586]	221
	222	! ecc2: square of yield ellipse eccenticrity
	223	ecc2 = rn_ecc * rn_ecc
	224	z1_ecc2 = 1._wp / ecc2
	225
	226	! alpha parameters (Bouillon 2009)
[8813]	227	IF( .NOT. ln_aEVP ) THEN
[13472]	228	zdtevp = rDt_ice / REAL( nn_nevp )
	229	zalph1 = 2._wp * rn_relast * REAL( nn_nevp )
[8813]	230	zalph2 = zalph1 * z1_ecc2
[8586]	231
[8813]	232	z1_alph1 = 1._wp / ( zalph1 + 1._wp )
	233	z1_alph2 = 1._wp / ( zalph2 + 1._wp )
[13472]	234	ELSE
	235	zdtevp = rdt_ice
	236	! zalpha parameters set later on adaptatively
[8813]	237	ENDIF
[13472]	238	z1_dtevp = 1._wp / zdtevp
[14072]	239
	240	! Initialise stress tensor
	241	zs1 (:,:) = pstress1_i (:,:)
[8586]	242	zs2 (:,:) = pstress2_i (:,:)
	243	zs12(:,:) = pstress12_i(:,:)
	244
	245	! Ice strength
	246	CALL ice_strength
	247
[10413]	248	! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
[13472]	249	IF( ln_landfast_L16 ) THEN ; zkt = rn_lf_tensile
[11536]	250	ELSE ; zkt = 0._wp
[10413]	251	ENDIF
[8586]	252	!
	253	!------------------------------------------------------------------------------!
	254	! 2) Wind / ocean stress, mass terms, coriolis terms
	255	!------------------------------------------------------------------------------!
[10415]	256	! sea surface height
	257	! embedded sea ice: compute representative ice top surface
	258	! non-embedded sea ice: use ocean surface for slope calculation
	259	zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
[8586]	260
[13295]	261	DO_2D( 0, 0, 0, 0 )
[8586]	262
[12377]	263	! ice fraction at U-V points
	264	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)
	265	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]	266
[12377]	267	! Ice/snow mass at U-V points
	268	zm1 = ( rhos * vt_s(ji ,jj ) + rhoi * vt_i(ji ,jj ) )
	269	zm2 = ( rhos * vt_s(ji+1,jj ) + rhoi * vt_i(ji+1,jj ) )
	270	zm3 = ( rhos * vt_s(ji ,jj+1) + rhoi * vt_i(ji ,jj+1) )
	271	zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
	272	zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
[8586]	273
[12377]	274	! Ocean currents at U-V points
	275	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)
	276	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]	277
[12377]	278	! Coriolis at T points (m*f)
	279	zmf(ji,jj) = zm1 * ff_t(ji,jj)
[8586]	280
[12377]	281	! dt/m at T points (for alpha and beta coefficients)
	282	zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin )
[14072]	283
[12377]	284	! m/dt
	285	zmU_t(ji,jj) = zmassU * z1_dtevp
	286	zmV_t(ji,jj) = zmassV * z1_dtevp
[14072]	287
[12377]	288	! Drag ice-atm.
	289	ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
	290	ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
[8586]	291
[12377]	292	! Surface pressure gradient (- mgGRAD(ssh)) at U-V points
	293	zspgU(ji,jj) = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
	294	zspgV(ji,jj) = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
[8586]	295
[12377]	296	! masks
	297	zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice
	298	zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice
[8586]	299
[12377]	300	! switches
	301	IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN ; zmsk01x(ji,jj) = 0._wp
	302	ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF
	303	IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN ; zmsk01y(ji,jj) = 0._wp
	304	ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF
[8586]	305
[12377]	306	END_2D
[14433]	307	CALL lbc_lnk( 'icedyn_rhg_evp', zmf, 'T', 1.0_wp, zdt_m, 'T', 1.0_wp )
[8586]	308	!
[10413]	309	! !== Landfast ice parameterization ==!
	310	!
	311	IF( ln_landfast_L16 ) THEN !-- Lemieux 2016
[13295]	312	DO_2D( 0, 0, 0, 0 )
[12377]	313	! ice thickness at U-V points
	314	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)
	315	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)
	316	! ice-bottom stress at U points
[14005]	317	zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0
[13472]	318	ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
[12377]	319	! ice-bottom stress at V points
[14005]	320	zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0
[13472]	321	ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
[12377]	322	! ice_bottom stress at T points
[14005]	323	zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0
[13472]	324	tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
[12377]	325	END_2D
[13226]	326	CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1.0_wp )
[10413]	327	!
[11536]	328	ELSE !-- no landfast
[13295]	329	DO_2D( 0, 0, 0, 0 )
[12377]	330	ztaux_base(ji,jj) = 0._wp
	331	ztauy_base(ji,jj) = 0._wp
	332	END_2D
[10413]	333	ENDIF
	334
[8586]	335	!------------------------------------------------------------------------------!
	336	! 3) Solution of the momentum equation, iterative procedure
	337	!------------------------------------------------------------------------------!
	338	!
[11536]	339	! ! ==================== !
[8586]	340	DO jter = 1 , nn_nevp ! loop over jter !
[14072]	341	! ! ==================== !
[10425]	342	l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
	343	!
[13472]	344	! convergence test
	345	IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN
	346	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
	347	zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step
	348	zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1)
	349	END_2D
	350	ENDIF
[8586]	351
	352	! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
[13295]	353	DO_2D( 1, 0, 1, 0 )
[8586]	354
[12377]	355	! shear at F points
	356	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) &
	357	& + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
[15020]	358	& ) * r1_e1e2f(ji,jj) * fimask(ji,jj)
[8586]	359
[12377]	360	END_2D
[8586]	361
[13550]	362	DO_2D( 0, 0, 0, 0 )
[8586]	363
[12377]	364	! shear**2 at T points (doc eq. A16)
	365	zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) &
	366	& + 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) &
	367	& ) * 0.25_wp * r1_e1e2t(ji,jj)
[14072]	368
[12377]	369	! divergence at T points
	370	zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
	371	& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
	372	& ) * r1_e1e2t(ji,jj)
	373	zdiv2 = zdiv * zdiv
[14072]	374
[12377]	375	! tension at T points
	376	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) &
	377	& - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
	378	& ) * r1_e1e2t(ji,jj)
	379	zdt2 = zdt * zdt
[14072]	380
[12377]	381	! delta at T points
[14072]	382	zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )
[8586]	383
[13550]	384	END_2D
	385	CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1.0_wp )
[8813]	386
[13550]	387	! P/delta at T points
	388	DO_2D( 1, 1, 1, 1 )
	389	zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl )
	390	END_2D
	391
	392	DO_2D( 0, 1, 0, 1 ) ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2
	393
	394	! divergence at T points (duplication to avoid communications)
	395	zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
	396	& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
	397	& ) * r1_e1e2t(ji,jj)
[14072]	398
[13550]	399	! tension at T points (duplication to avoid communications)
	400	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) &
	401	& - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
	402	& ) * r1_e1e2t(ji,jj)
[14072]	403
[13472]	404	! alpha for aEVP
[12377]	405	! gamma = 0.5P/(delta+creepl) (cpi)2/Area dt/m
	406	! alpha = beta = sqrt(4*gamma)
	407	IF( ln_aEVP ) THEN
	408	zalph1 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
	409	z1_alph1 = 1._wp / ( zalph1 + 1._wp )
	410	zalph2 = zalph1
	411	z1_alph2 = z1_alph1
[13472]	412	! explicit:
	413	! z1_alph1 = 1._wp / zalph1
	414	! z1_alph2 = 1._wp / zalph1
	415	! zalph1 = zalph1 - 1._wp
	416	! zalph2 = zalph1
[12377]	417	ENDIF
[14072]	418
[12377]	419	! stress at T points (zkt/=0 if landfast)
[13550]	420	zs1(ji,jj) = ( zs1(ji,jj)zalph1 + zp_delt(ji,jj) ( zdiv(1._wp + zkt) - zdelta(ji,jj)(1._wp - zkt) ) ) * z1_alph1
	421	zs2(ji,jj) = ( zs2(ji,jj)zalph2 + zp_delt(ji,jj) ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
[14072]	422
[12377]	423	END_2D
[8586]	424
[13472]	425	! Save beta at T-points for further computations
	426	IF( ln_aEVP ) THEN
	427	DO_2D( 1, 1, 1, 1 )
	428	zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
	429	END_2D
	430	ENDIF
[14072]	431
[13295]	432	DO_2D( 1, 0, 1, 0 )
[8586]	433
[13472]	434	! alpha for aEVP
[12377]	435	IF( ln_aEVP ) THEN
[13472]	436	zalph2 = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) )
[12377]	437	z1_alph2 = 1._wp / ( zalph2 + 1._wp )
[13472]	438	! explicit:
	439	! z1_alph2 = 1._wp / zalph2
	440	! zalph2 = zalph2 - 1._wp
[12377]	441	ENDIF
[14072]	442
[12377]	443	! P/delta at F points
	444	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) )
[14072]	445
[12377]	446	! stress at F points (zkt/=0 if landfast)
	447	zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
[8586]	448
[12377]	449	END_2D
[8586]	450
	451	! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
[13295]	452	DO_2D( 0, 0, 0, 0 )
[12377]	453	! !--- U points
	454	zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) &
	455	& + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) &
	456	& ) * r1_e2u(ji,jj) &
	457	& + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) &
	458	& ) * 2._wp * r1_e1u(ji,jj) &
	459	& ) * r1_e1e2u(ji,jj)
	460	!
	461	! !--- V points
	462	zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) &
	463	& - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) &
	464	& ) * r1_e1v(ji,jj) &
	465	& + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) &
	466	& ) * 2._wp * r1_e2v(ji,jj) &
	467	& ) * r1_e1e2v(ji,jj)
	468	!
	469	! !--- ice currents at U-V point
	470	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)
	471	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)
	472	!
	473	END_2D
[8586]	474	!
	475	! --- Computation of ice velocity --- !
[8813]	476	! Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
[8586]	477	! Bouillon et al. 2009 (eq 34-35) => stable
	478	IF( MOD(jter,2) == 0 ) THEN ! even iterations
	479	!
[13295]	480	DO_2D( 0, 0, 0, 0 )
[12377]	481	! !--- tau_io/(v_oce - v_ice)
	482	zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) &
	483	& + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
	484	! !--- Ocean-to-Ice stress
	485	ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
	486	!
	487	! !--- tau_bottom/v_ice
	488	zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
	489	zTauB = ztauy_base(ji,jj) / zvel
	490	! !--- OceanBottom-to-Ice stress
	491	ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
	492	!
	493	! !--- Coriolis at V-points (energy conserving formulation)
	494	zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * &
	495	& ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) &
	496	& + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
	497	!
	498	! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
	499	zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
	500	!
	501	! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS)
	502	! 1 = sliding friction : TauB < RHS
	503	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
	504	!
	505	IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
[13472]	506	zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
	507	v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity
	508	& + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	509	& ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
[14072]	510	& + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) &
[13472]	511	& + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0
	512	& ) / ( zbetav + 1._wp ) &
	513	& ) * 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
[12377]	514	& ) * zmsk00y(ji,jj)
	515	ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
[13472]	516	v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity
	517	& + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	518	& ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	519	& + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0
	520	& ) * 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
	521	& ) * zmsk00y(ji,jj)
[12377]	522	ENDIF
	523	END_2D
[13226]	524	CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp )
[8586]	525	!
	526	#if defined key_agrif
[9610]	527	!! CALL agrif_interp_ice( 'V', jter, nn_nevp )
	528	CALL agrif_interp_ice( 'V' )
[8586]	529	#endif
[11536]	530	IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
[8586]	531	!
[13295]	532	DO_2D( 0, 0, 0, 0 )
[12377]	533	! !--- tau_io/(u_oce - u_ice)
	534	zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) &
	535	& + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
	536	! !--- Ocean-to-Ice stress
	537	ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
	538	!
	539	! !--- tau_bottom/u_ice
	540	zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
	541	zTauB = ztaux_base(ji,jj) / zvel
	542	! !--- OceanBottom-to-Ice stress
	543	ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
	544	!
	545	! !--- Coriolis at U-points (energy conserving formulation)
	546	zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * &
	547	& ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) &
	548	& + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
	549	!
	550	! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
	551	zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
	552	!
	553	! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS)
	554	! 1 = sliding friction : TauB < RHS
	555	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
	556	!
	557	IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
[13472]	558	zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
	559	u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity
	560	& + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	561	& ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	562	& + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) &
	563	& + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0
	564	& ) / ( zbetau + 1._wp ) &
[14072]	565	& ) * 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
[12377]	566	& ) * zmsk00x(ji,jj)
	567	ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
[13472]	568	u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity
	569	& + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	570	& ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	571	& + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0
[14072]	572	& ) * 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
[13472]	573	& ) * zmsk00x(ji,jj)
[12377]	574	ENDIF
	575	END_2D
[13226]	576	CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp )
[8586]	577	!
	578	#if defined key_agrif
[9610]	579	!! CALL agrif_interp_ice( 'U', jter, nn_nevp )
	580	CALL agrif_interp_ice( 'U' )
[8586]	581	#endif
[11536]	582	IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
[8586]	583	!
	584	ELSE ! odd iterations
	585	!
[13295]	586	DO_2D( 0, 0, 0, 0 )
[12377]	587	! !--- tau_io/(u_oce - u_ice)
	588	zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) &
	589	& + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
	590	! !--- Ocean-to-Ice stress
	591	ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
	592	!
	593	! !--- tau_bottom/u_ice
	594	zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
	595	zTauB = ztaux_base(ji,jj) / zvel
	596	! !--- OceanBottom-to-Ice stress
	597	ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
	598	!
	599	! !--- Coriolis at U-points (energy conserving formulation)
	600	zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * &
	601	& ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) &
	602	& + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
	603	!
	604	! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
	605	zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
	606	!
	607	! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS)
	608	! 1 = sliding friction : TauB < RHS
	609	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
	610	!
	611	IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
[13472]	612	zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
	613	u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity
	614	& + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	615	& ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	616	& + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) &
	617	& + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0
	618	& ) / ( zbetau + 1._wp ) &
[14072]	619	& ) * 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
[12377]	620	& ) * zmsk00x(ji,jj)
	621	ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
[13472]	622	u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity
	623	& + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	624	& ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	625	& + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0
	626	& ) * 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
	627	& ) * zmsk00x(ji,jj)
[12377]	628	ENDIF
	629	END_2D
[13226]	630	CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp )
[8586]	631	!
	632	#if defined key_agrif
[9610]	633	!! CALL agrif_interp_ice( 'U', jter, nn_nevp )
	634	CALL agrif_interp_ice( 'U' )
[8586]	635	#endif
[11536]	636	IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
[8586]	637	!
[13295]	638	DO_2D( 0, 0, 0, 0 )
[12377]	639	! !--- tau_io/(v_oce - v_ice)
	640	zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) &
	641	& + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
	642	! !--- Ocean-to-Ice stress
	643	ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
	644	!
	645	! !--- tau_bottom/v_ice
	646	zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
	647	zTauB = ztauy_base(ji,jj) / zvel
	648	! !--- OceanBottom-to-Ice stress
	649	ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
	650	!
	651	! !--- Coriolis at v-points (energy conserving formulation)
	652	zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * &
	653	& ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) &
	654	& + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
	655	!
	656	! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
	657	zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
	658	!
	659	! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS)
	660	! 1 = sliding friction : TauB < RHS
	661	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
	662	!
	663	IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
[13472]	664	zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
	665	v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity
	666	& + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	667	& ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	668	& + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) &
	669	& + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0
[14072]	670	& ) / ( zbetav + 1._wp ) &
[13472]	671	& ) * 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
[12377]	672	& ) * zmsk00y(ji,jj)
	673	ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
[13472]	674	v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity
	675	& + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	676	& ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	677	& + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0
	678	& ) * 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
	679	& ) * zmsk00y(ji,jj)
[12377]	680	ENDIF
	681	END_2D
[13226]	682	CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp )
[8586]	683	!
	684	#if defined key_agrif
[9610]	685	!! CALL agrif_interp_ice( 'V', jter, nn_nevp )
	686	CALL agrif_interp_ice( 'V' )
[8586]	687	#endif
[11536]	688	IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
[8586]	689	!
	690	ENDIF
[10425]	691
[13472]	692	! convergence test
[15020]	693	IF( nn_rhg_chkcvg == 2 ) CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice, zmsk15 )
[8586]	694	!
	695	! ! ==================== !
	696	END DO ! end loop over jter !
	697	! ! ==================== !
[13472]	698	IF( ln_aEVP ) CALL iom_put( 'beta_evp' , zbeta )
[8586]	699	!
	700	!------------------------------------------------------------------------------!
[14072]	701	! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
[8586]	702	!------------------------------------------------------------------------------!
[13295]	703	DO_2D( 1, 0, 1, 0 )
[8586]	704
[12377]	705	! shear at F points
	706	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) &
	707	& + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
[15020]	708	& ) * r1_e1e2f(ji,jj) * fimask(ji,jj)
[8586]	709
[12377]	710	END_2D
[14072]	711
[13497]	712	DO_2D( 0, 0, 0, 0 ) ! no vector loop
[14072]	713
[12377]	714	! tension**2 at T points
	715	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) &
	716	& - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
	717	& ) * r1_e1e2t(ji,jj)
	718	zdt2 = zdt * zdt
[13601]	719
	720	zten_i(ji,jj) = zdt
[14072]	721
[12377]	722	! shear**2 at T points (doc eq. A16)
	723	zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) &
	724	& + 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) &
	725	& ) * 0.25_wp * r1_e1e2t(ji,jj)
[14072]	726
[12377]	727	! shear at T points
	728	pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
[8586]	729
[12377]	730	! divergence at T points
	731	pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
	732	& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
	733	& ) * r1_e1e2t(ji,jj)
[14072]	734
[12377]	735	! delta at T points
[14072]	736	zfac = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta
[13601]	737	rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0
	738	pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl
[8586]	739
[12377]	740	END_2D
[14433]	741	CALL lbc_lnk( 'icedyn_rhg_evp', pshear_i, 'T', 1._wp, pdivu_i, 'T', 1._wp, pdelta_i, 'T', 1._wp, zten_i, 'T', 1._wp, &
	742	& zs1 , 'T', 1._wp, zs2 , 'T', 1._wp, zs12 , 'F', 1._wp )
[14072]	743
[8586]	744	! --- Store the stress tensor for the next time step --- !
	745	pstress1_i (:,:) = zs1 (:,:)
	746	pstress2_i (:,:) = zs2 (:,:)
	747	pstress12_i(:,:) = zs12(:,:)
	748	!
	749
	750	!------------------------------------------------------------------------------!
	751	! 5) diagnostics
	752	!------------------------------------------------------------------------------!
[11536]	753	! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
	754	IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
	755	& iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
	756	!
[14433]	757	CALL lbc_lnk( 'icedyn_rhg_evp', ztaux_oi, 'U', -1.0_wp, ztauy_oi, 'V', -1.0_wp, &
	758	& ztaux_ai, 'U', -1.0_wp, ztauy_ai, 'V', -1.0_wp, &
	759	& ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp )
[11536]	760	!
	761	CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
	762	CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
	763	CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
	764	CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
	765	CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
	766	CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
	767	ENDIF
[14072]	768
[8586]	769	! --- divergence, shear and strength --- !
[11536]	770	IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence
	771	IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear
	772	IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength
[8586]	773
[13601]	774	! --- Stress tensor invariants (SIMIP diags) --- !
	775	IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
[8586]	776	!
[13601]	777	ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
[14072]	778	!
[13601]	779	DO_2D( 1, 1, 1, 1 )
[14072]	780
[13601]	781	! Ice stresses
	782	! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013)
	783	! These are NOT stress tensor components, neither stress invariants, neither stress principal components
	784	! I know, this can be confusing...
[14072]	785	zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl )
[13601]	786	zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
	787	zsig2 = zfac * z1_ecc2 * zten_i(ji,jj)
	788	zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj)
[14072]	789
[13601]	790	! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008)
[13612]	791	zsig_I (ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure
	792	zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress
[14072]	793
	794	END_2D
[13601]	795	!
	796	! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008)
	797	IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress
	798	IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
[14072]	799
[13601]	800	DEALLOCATE ( zsig_I, zsig_II )
[14072]	801
[13601]	802	ENDIF
[8586]	803
[13601]	804	! --- Normalized stress tensor principal components --- !
	805	! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020
	806	! Recommendation 1 : we use ice strength, not replacement pressure
	807	! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities
	808	IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
[8586]	809	!
[14072]	810	ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
	811	!
[13601]	812	DO_2D( 1, 1, 1, 1 )
[14072]	813
	814	! Ice stresses computed with viscosities (delta, p/delta) at previous iterates
[13601]	815	! and deformations at current iterates
	816	! following Lemieux & Dupont (2020)
	817	zfac = zp_delt(ji,jj)
	818	zsig1 = zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) )
	819	zsig2 = zfac * z1_ecc2 * zten_i(ji,jj)
	820	zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj)
[14072]	821
[13601]	822	! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point
[13612]	823	zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure
	824	zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress
[14072]	825
[13601]	826	! Normalized principal stresses (used to display the ellipse)
	827	z1_strength = 1._wp / MAX( 1._wp, strength(ji,jj) )
	828	zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength
	829	zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
[14072]	830	END_2D
[8586]	831	!
[14072]	832	CALL iom_put( 'sig1_pnorm' , zsig1_p )
	833	CALL iom_put( 'sig2_pnorm' , zsig2_p )
[11536]	834
[13601]	835	DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II )
[14072]	836
[8586]	837	ENDIF
[13472]	838
[8586]	839	! --- SIMIP --- !
[11536]	840	IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
	841	& iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
	842	!
[14433]	843	CALL lbc_lnk( 'icedyn_rhg_evp', zspgU, 'U', -1.0_wp, zspgV, 'V', -1.0_wp, &
	844	& zCorU, 'U', -1.0_wp, zCorV, 'V', -1.0_wp, zfU, 'U', -1.0_wp, zfV, 'V', -1.0_wp )
[8586]	845
[11536]	846	CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x)
	847	CALL iom_put( 'dssh_dy' , zspgV * zmsk00 ) ! Sea-surface tilt term in force balance (y)
	848	CALL iom_put( 'corstrx' , zCorU * zmsk00 ) ! Coriolis force term in force balance (x)
	849	CALL iom_put( 'corstry' , zCorV * zmsk00 ) ! Coriolis force term in force balance (y)
	850	CALL iom_put( 'intstrx' , zfU * zmsk00 ) ! Internal force term in force balance (x)
	851	CALL iom_put( 'intstry' , zfV * zmsk00 ) ! Internal force term in force balance (y)
	852	ENDIF
	853
	854	IF( iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
	855	& iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
	856	!
	857	ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
	858	& zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
	859	!
[13295]	860	DO_2D( 0, 0, 0, 0 )
[12377]	861	! 2D ice mass, snow mass, area transport arrays (X, Y)
	862	zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
	863	zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)
[11536]	864
[12377]	865	zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
	866	zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) ! '' Y- ''
[11536]	867
[12377]	868	zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
	869	zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) ! '' Y- ''
[11536]	870
[12377]	871	zdiag_xatrp(ji,jj) = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) ) ! area transport, X-component
	872	zdiag_yatrp(ji,jj) = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) ) ! '' Y- ''
[11536]	873
[12377]	874	END_2D
[8586]	875
[14433]	876	CALL lbc_lnk( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1.0_wp, zdiag_ymtrp_ice, 'V', -1.0_wp, &
	877	& zdiag_xmtrp_snw, 'U', -1.0_wp, zdiag_ymtrp_snw, 'V', -1.0_wp, &
	878	& zdiag_xatrp , 'U', -1.0_wp, zdiag_yatrp , 'V', -1.0_wp )
[8586]	879
[11536]	880	CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice ) ! X-component of sea-ice mass transport (kg/s)
[14072]	881	CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice ) ! Y-component of sea-ice mass transport
[11536]	882	CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw ) ! X-component of snow mass transport (kg/s)
	883	CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw ) ! Y-component of snow mass transport
	884	CALL iom_put( 'xatrp' , zdiag_xatrp ) ! X-component of ice area transport
	885	CALL iom_put( 'yatrp' , zdiag_yatrp ) ! Y-component of ice area transport
	886
	887	DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , &
	888	& zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp )
	889
[8586]	890	ENDIF
	891	!
[13472]	892	! --- convergence tests --- !
	893	IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN
	894	IF( iom_use('uice_cvg') ) THEN
	895	IF( ln_aEVP ) THEN ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
	896	CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , &
	897	& ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) )
	898	ELSE ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
	899	CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , &
	900	& ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )
	901	ENDIF
	902	ENDIF
[14072]	903	ENDIF
[13472]	904	!
[8586]	905	END SUBROUTINE ice_dyn_rhg_evp
	906
[8813]	907
[15020]	908	SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb, pmsk15 )
[13472]	909	!!----------------------------------------------------------------------
	910	!! * ROUTINE rhg_cvg *
[14072]	911	!!
[13472]	912	!! ** Purpose : check convergence of oce rheology
	913	!!
	914	!! ** Method : create a file ice_cvg.nc containing the convergence of ice velocity
	915	!! during the sub timestepping of rheology so as:
	916	!! uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) )
	917	!! This routine is called every sub-iteration, so it is cpu expensive
	918	!!
[14072]	919	!! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)
[13472]	920	!!----------------------------------------------------------------------
	921	INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index
	922	REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pub, pvb ! now and before velocities
[15020]	923	REAL(wp), DIMENSION(:,:), INTENT(in) :: pmsk15
[13472]	924	!!
	925	INTEGER :: it, idtime, istatus
	926	INTEGER :: ji, jj ! dummy loop indices
[14072]	927	REAL(wp) :: zresm ! local real
[13472]	928	CHARACTER(len=20) :: clname
	929	!!----------------------------------------------------------------------
	930
	931	! create file
	932	IF( kt == nit000 .AND. kiter == 1 ) THEN
	933	!
	934	IF( lwp ) THEN
	935	WRITE(numout,*)
	936	WRITE(numout,*) 'rhg_cvg : ice rheology convergence control'
	937	WRITE(numout,*) '~~~~~~~'
	938	ENDIF
	939	!
	940	IF( lwm ) THEN
	941	clname = 'ice_cvg.nc'
	942	IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
	943	istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid )
	944	istatus = NF90_DEF_DIM( ncvgid, 'time' , NF90_UNLIMITED, idtime )
	945	istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid )
	946	istatus = NF90_ENDDEF(ncvgid)
	947	ENDIF
	948	!
	949	ENDIF
	950
	951	! time
	952	it = ( kt - 1 ) * kitermax + kiter
[14072]	953
[13472]	954	! convergence
	955	IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large)
	956	zresm = 0._wp
	957	ELSE
[15020]	958	zresm = 0._wp
	959	DO_2D( 0, 0, 0, 0 )
	960	zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
	961	& ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) )
[13472]	962	END_2D
	963	CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain
	964	ENDIF
	965
	966	IF( lwm ) THEN
	967	! write variables
	968	istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) )
	969	! close file
[13601]	970	IF( kt == nitend - nn_fsbc + 1 ) istatus = NF90_CLOSE(ncvgid)
[13472]	971	ENDIF
[14072]	972
[13472]	973	END SUBROUTINE rhg_cvg
	974
	975
[8586]	976	SUBROUTINE rhg_evp_rst( cdrw, kt )
	977	!!---------------------------------------------------------------------
	978	!! * ROUTINE rhg_evp_rst *
[14072]	979	!!
[8586]	980	!! ** Purpose : Read or write RHG file in restart file
	981	!!
	982	!! ** Method : use of IOM library
	983	!!----------------------------------------------------------------------
	984	CHARACTER(len=*) , INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	985	INTEGER, OPTIONAL, INTENT(in) :: kt ! ice time-step
	986	!
	987	INTEGER :: iter ! local integer
	988	INTEGER :: id1, id2, id3 ! local integers
	989	!!----------------------------------------------------------------------
	990	!
	991	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialize
	992	! ! ---------------
	993	IF( ln_rstart ) THEN !* Read the restart file
	994	!
	995	id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
	996	id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
	997	id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
	998	!
	999	IF( MIN( id1, id2, id3 ) > 0 ) THEN ! fields exist
[13286]	1000	CALL iom_get( numrir, jpdom_auto, 'stress1_i' , stress1_i , cd_type = 'T' )
	1001	CALL iom_get( numrir, jpdom_auto, 'stress2_i' , stress2_i , cd_type = 'T' )
	1002	CALL iom_get( numrir, jpdom_auto, 'stress12_i', stress12_i, cd_type = 'F' )
[8586]	1003	ELSE ! start rheology from rest
[9169]	1004	IF(lwp) WRITE(numout,*)
	1005	IF(lwp) WRITE(numout,*) ' ==>>> previous run without rheology, set stresses to 0'
[8586]	1006	stress1_i (:,:) = 0._wp
	1007	stress2_i (:,:) = 0._wp
	1008	stress12_i(:,:) = 0._wp
	1009	ENDIF
	1010	ELSE !* Start from rest
[9169]	1011	IF(lwp) WRITE(numout,*)
	1012	IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set stresses to 0'
[8586]	1013	stress1_i (:,:) = 0._wp
	1014	stress2_i (:,:) = 0._wp
	1015	stress12_i(:,:) = 0._wp
	1016	ENDIF
	1017	!
	1018	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
	1019	! ! -------------------
	1020	IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
	1021	iter = kt + nn_fsbc - 1 ! ice restarts are written at kt == nitrst - nn_fsbc + 1
	1022	!
[13970]	1023	CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i )
	1024	CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i )
[8586]	1025	CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
	1026	!
	1027	ENDIF
	1028	!
	1029	END SUBROUTINE rhg_evp_rst
	1030
[14072]	1031
[8586]	1032	#else
	1033	!!----------------------------------------------------------------------
[9570]	1034	!! Default option Empty module NO SI3 sea-ice model
[8586]	1035	!!----------------------------------------------------------------------
	1036	#endif
	1037
	1038	!!==============================================================================
	1039	END MODULE icedyn_rhg_evp

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

Download in other formats: