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

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_rhg_evp.F90 @ 12236

Last change on this file since 12236 was 12236, checked in by acc, 4 years ago
Branch 2019/dev_r11943_MERGE_2019. Merge in changes from 2019/fix_sn_cfctl_ticket2328. Fully SETTE tested
Property svn:keywords set to `Id`
File size: 55.7 KB

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

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

Download in other formats: