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

Context Navigation

source: NEMO/trunk/src/ICE/icedyn_rhg_evp.F90 @ 14558

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

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

Download in other formats: