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.

limrhg.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 8332

Last change on this file since 8332 was 8324, checked in by clem, 7 years ago
STEP3 (2): clean separation between sea-ice and ocean
Property svn:keywords set to `Id`
File size: 42.5 KB

Rev	Line
[825]	1	MODULE limrhg
	2	!!======================================================================
	3	!! * MODULE limrhg *
[834]	4	!! Ice rheology : sea ice rheology
[825]	5	!!======================================================================
[1244]	6	!! History : - ! 2007-03 (M.A. Morales Maqueda, S. Bouillon) Original code
	7	!! 3.0 ! 2008-03 (M. Vancoppenolle) LIM3
	8	!! - ! 2008-11 (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy
[8306]	9	!! 3.3 ! 2009-05 (G.Garric) addition of the evp cas
[3680]	10	!! 3.4 ! 2011-01 (A. Porter) dynamical allocation
[7646]	11	!! 3.5 ! 2012-08 (R. Benshila) AGRIF
	12	!! 3.6 ! 2016-06 (C. Rousset) Rewriting + landfast ice + possibility to use mEVP (Bouillon 2013)
[1244]	13	!!----------------------------------------------------------------------
[7646]	14	#if defined key_lim3
[825]	15	!!----------------------------------------------------------------------
[7646]	16	!! 'key_lim3' LIM-3 sea-ice model
[825]	17	!!----------------------------------------------------------------------
[3625]	18	!! lim_rhg : computes ice velocities
[825]	19	!!----------------------------------------------------------------------
[3625]	20	USE phycst ! Physical constant
	21	USE par_oce ! Ocean parameters
	22	USE dom_oce ! Ocean domain
[8316]	23	USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
[8324]	24	USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
[7646]	25	USE ice ! ice variables
	26	USE limitd_me ! ice strength
[3625]	27	USE lbclnk ! Lateral Boundary Condition / MPP link
	28	USE lib_mpp ! MPP library
	29	USE wrk_nemo ! work arrays
	30	USE in_out_manager ! I/O manager
	31	USE prtctl ! Print control
	32	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[7646]	33	#if defined key_agrif
	34	USE agrif_lim3_interp
[3680]	35	#endif
[7646]	36	USE bdy_oce , ONLY: ln_bdy
	37	USE bdyice_lim
[825]	38
	39	IMPLICIT NONE
	40	PRIVATE
	41
[7646]	42	PUBLIC lim_rhg ! routine called by lim_dyn
[825]	43
[868]	44	!! * Substitutions
	45	# include "vectopt_loop_substitute.h90"
[825]	46	!!----------------------------------------------------------------------
[4161]	47	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]	48	!! $Id$
[2528]	49	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	50	!!----------------------------------------------------------------------
	51	CONTAINS
	52
[7646]	53	SUBROUTINE lim_rhg
[825]	54	!!-------------------------------------------------------------------
[834]	55	!! * SUBROUTINE lim_rhg *
	56	!! EVP-C-grid
[825]	57	!!
[834]	58	!! ** purpose : determines sea ice drift from wind stress, ice-ocean
[825]	59	!! stress and sea-surface slope. Ice-ice interaction is described by
[834]	60	!! a non-linear elasto-viscous-plastic (EVP) law including shear
	61	!! strength and a bulk rheology (Hunke and Dukowicz, 2002).
[825]	62	!!
[834]	63	!! The points in the C-grid look like this, dear reader
[825]	64	!!
[834]	65	!! (ji,jj)
	66	!! \|
	67	!! \|
	68	!! (ji-1,jj) \| (ji,jj)
	69	!! ---------
	70	!! \| \|
	71	!! \| (ji,jj) \|------(ji,jj)
	72	!! \| \|
	73	!! ---------
	74	!! (ji-1,jj-1) (ji,jj-1)
[825]	75	!!
[834]	76	!! ** Inputs : - wind forcing (stress), oceanic currents
	77	!! ice total volume (vt_i) per unit area
	78	!! snow total volume (vt_s) per unit area
[825]	79	!!
[834]	80	!! ** Action : - compute u_ice, v_ice : the components of the
	81	!! sea-ice velocity vector
	82	!! - compute delta_i, shear_i, divu_i, which are inputs
	83	!! of the ice thickness distribution
[825]	84	!!
[834]	85	!! ** Steps : 1) Compute ice snow mass, ice strength
	86	!! 2) Compute wind, oceanic stresses, mass terms and
	87	!! coriolis terms of the momentum equation
	88	!! 3) Solve the momentum equation (iterative procedure)
	89	!! 4) Prevent high velocities if the ice is thin
	90	!! 5) Recompute invariants of the strain rate tensor
	91	!! which are inputs of the ITD, store stress
	92	!! for the next time step
	93	!! 6) Control prints of residual (convergence)
	94	!! and charge ellipse.
	95	!! The user should make sure that the parameters
[5123]	96	!! nn_nevp, elastic time scale and rn_creepl maintain stress state
[834]	97	!! on the charge ellipse for plastic flow
	98	!! e.g. in the Canadian Archipelago
	99	!!
[7646]	100	!! ** Notes : There is the possibility to use mEVP from Bouillon 2013
	101	!! (by uncommenting some lines in part 3 and changing alpha and beta parameters)
	102	!! but this solution appears very unstable (see Kimmritz et al 2016)
	103	!!
[2528]	104	!! References : Hunke and Dukowicz, JPO97
	105	!! Bouillon et al., Ocean Modelling 2009
[7646]	106	!! Bouillon et al., Ocean Modelling 2013
[2528]	107	!!-------------------------------------------------------------------
[7646]	108	INTEGER :: ji, jj ! dummy loop indices
	109	INTEGER :: jter ! local integers
[825]	110	CHARACTER (len=50) :: charout
	111
[7646]	112	REAL(wp) :: zrhoco ! rau0 * rn_cio
	113	REAL(wp) :: zdtevp, z1_dtevp ! time step for subcycling
	114	REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity
	115	REAL(wp) :: zbeta, zalph1, z1_alph1, zalph2, z1_alph2 ! alpha and beta from Bouillon 2009 and 2013
	116	REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV ! ice/snow mass
	117	REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars
[8239]	118	REAL(wp) :: zTauO, zTauB, zTauE, zvel ! temporary scalars
[825]	119
[7646]	120	REAL(wp) :: zsig1, zsig2 ! internal ice stress
	121	REAL(wp) :: zresm ! Maximal error on ice velocity
	122	REAL(wp) :: zintb, zintn ! dummy argument
[8239]	123	REAL(wp) :: zfac_x, zfac_y
[3294]	124
[7646]	125	REAL(wp), POINTER, DIMENSION(:,:) :: z1_e1t0, z1_e2t0 ! scale factors
	126	REAL(wp), POINTER, DIMENSION(:,:) :: zp_delt ! P/delta at T points
	127	!
	128	REAL(wp), POINTER, DIMENSION(:,:) :: zaU , zaV ! ice fraction on U/V points
	129	REAL(wp), POINTER, DIMENSION(:,:) :: zmU_t, zmV_t ! ice/snow mass/dt on U/V points
	130	REAL(wp), POINTER, DIMENSION(:,:) :: zmf ! coriolis parameter at T points
	131	REAL(wp), POINTER, DIMENSION(:,:) :: zTauU_ia , ztauV_ia ! ice-atm. stress at U-V points
	132	REAL(wp), POINTER, DIMENSION(:,:) :: zspgU , zspgV ! surface pressure gradient at U/V points
	133	REAL(wp), POINTER, DIMENSION(:,:) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points
	134	REAL(wp), POINTER, DIMENSION(:,:) :: zfU , zfV ! internal stresses
	135
	136	REAL(wp), POINTER, DIMENSION(:,:) :: zds ! shear
	137	REAL(wp), POINTER, DIMENSION(:,:) :: zs1, zs2, zs12 ! stress tensor components
	138	REAL(wp), POINTER, DIMENSION(:,:) :: zu_ice, zv_ice, zresr ! check convergence
	139	REAL(wp), POINTER, DIMENSION(:,:) :: zpice ! array used for the calculation of ice surface slope:
	140	! ocean surface (ssh_m) if ice is not embedded
	141	! ice top surface if ice is embedded
[8239]	142	REAL(wp), POINTER, DIMENSION(:,:) :: zCorx, zCory ! Coriolis stress array
	143	REAL(wp), POINTER, DIMENSION(:,:) :: ztaux_oi, ztauy_oi ! Ocean-to-ice stress array
	144
[7646]	145	REAL(wp), POINTER, DIMENSION(:,:) :: zswitchU, zswitchV ! dummy arrays
	146	REAL(wp), POINTER, DIMENSION(:,:) :: zmaskU, zmaskV ! mask for ice presence
	147	REAL(wp), POINTER, DIMENSION(:,:) :: zfmask, zwf ! mask at F points for the ice
[5123]	148
[7646]	149	REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter
	150	REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity equals ocean velocity
[2528]	151	!!-------------------------------------------------------------------
[3294]	152
[7646]	153	CALL wrk_alloc( jpi,jpj, z1_e1t0, z1_e2t0, zp_delt )
	154	CALL wrk_alloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
	155	CALL wrk_alloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
	156	CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
	157	CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
[8239]	158	CALL wrk_alloc( jpi,jpj, zCorx, zCory)
	159	CALL wrk_alloc( jpi,jpj, ztaux_oi, ztauy_oi)
[3294]	160
[7646]	161	#if defined key_agrif
	162	CALL agrif_interp_lim3( 'U', 0, nn_nevp ) ! First interpolation of coarse values
	163	CALL agrif_interp_lim3( 'V', 0, nn_nevp )
[2528]	164	#endif
[921]	165	!
	166	!------------------------------------------------------------------------------!
[7646]	167	! 0) mask at F points for the ice
[921]	168	!------------------------------------------------------------------------------!
[7646]	169	! ocean/land mask
	170	DO jj = 1, jpjm1
	171	DO ji = 1, jpim1 ! NO vector opt.
	172	zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
	173	END DO
	174	END DO
	175	CALL lbc_lnk( zfmask, 'F', 1._wp )
[825]	176
[7646]	177	! Lateral boundary conditions on velocity (modify zfmask)
[7753]	178	zwf(:,:) = zfmask(:,:)
[7646]	179	DO jj = 2, jpjm1
	180	DO ji = fs_2, fs_jpim1 ! vector opt.
	181	IF( zfmask(ji,jj) == 0._wp ) THEN
[8313]	182	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) ) )
[7646]	183	ENDIF
[825]	184	END DO
	185	END DO
[7646]	186	DO jj = 2, jpjm1
	187	IF( zfmask(1,jj) == 0._wp ) THEN
[8313]	188	zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
[7646]	189	ENDIF
	190	IF( zfmask(jpi,jj) == 0._wp ) THEN
[8313]	191	zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
[7646]	192	ENDIF
	193	END DO
	194	DO ji = 2, jpim1
	195	IF( zfmask(ji,1) == 0._wp ) THEN
[8313]	196	zfmask(ji,1 ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
[7646]	197	ENDIF
	198	IF( zfmask(ji,jpj) == 0._wp ) THEN
[8313]	199	zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
[7646]	200	ENDIF
	201	END DO
	202	CALL lbc_lnk( zfmask, 'F', 1._wp )
[825]	203
[7646]	204	!------------------------------------------------------------------------------!
	205	! 1) define some variables and initialize arrays
	206	!------------------------------------------------------------------------------!
	207	zrhoco = rau0 * rn_cio
	208
	209	! ecc2: square of yield ellipse eccenticrity
	210	ecc2 = rn_ecc * rn_ecc
	211	z1_ecc2 = 1._wp / ecc2
	212
	213	! Time step for subcycling
	214	zdtevp = rdt_ice / REAL( nn_nevp )
	215	z1_dtevp = 1._wp / zdtevp
	216
	217	! alpha parameters (Bouillon 2009)
	218	zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
	219	zalph2 = zalph1 * z1_ecc2
	220
	221	! alpha and beta parameters (Bouillon 2013)
	222	!!zalph1 = 40.
	223	!!zalph2 = 40.
	224	!!zbeta = 3000.
	225	!!zbeta = REAL( nn_nevp ) ! close to classical EVP of Hunke (2001)
	226
	227	z1_alph1 = 1._wp / ( zalph1 + 1._wp )
	228	z1_alph2 = 1._wp / ( zalph2 + 1._wp )
	229
	230	! Initialise stress tensor
[7753]	231	zs1 (:,:) = stress1_i (:,:)
	232	zs2 (:,:) = stress2_i (:,:)
	233	zs12(:,:) = stress12_i(:,:)
[7646]	234
	235	! Ice strength
	236	CALL lim_itd_me_icestrength( nn_icestr )
	237
	238	! scale factors
	239	DO jj = 2, jpjm1
	240	DO ji = fs_2, fs_jpim1
	241	z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj ) + e1t(ji,jj ) )
	242	z1_e2t0(ji,jj) = 1._wp / ( e2t(ji ,jj+1) + e2t(ji,jj ) )
[825]	243	END DO
	244	END DO
[7646]	245
[921]	246	!
	247	!------------------------------------------------------------------------------!
	248	! 2) Wind / ocean stress, mass terms, coriolis terms
	249	!------------------------------------------------------------------------------!
	250
[8313]	251	IF( ln_ice_embd ) THEN !== embedded sea ice: compute representative ice top surface ==!
[5123]	252	!
	253	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
	254	! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
[3625]	255	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
[5123]	256	!
	257	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
	258	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
[3625]	259	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
[5123]	260	!
[7753]	261	zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
[5123]	262	!
[3625]	263	ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==!
[7753]	264	zpice(:,:) = ssh_m(:,:)
[3625]	265	ENDIF
	266
[7646]	267	DO jj = 2, jpjm1
[868]	268	DO ji = fs_2, fs_jpim1
[825]	269
[7646]	270	! ice fraction at U-V points
	271	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)
	272	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)
[4990]	273
[7646]	274	! Ice/snow mass at U-V points
	275	zm1 = ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) )
	276	zm2 = ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) )
	277	zm3 = ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) )
	278	zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
	279	zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
[825]	280
[7646]	281	! Ocean currents at U-V points
	282	v_oceU(ji,jj) = 0.5_wp * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji+1,jj) &
	283	& + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
	284
	285	u_oceV(ji,jj) = 0.5_wp * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj+1) &
	286	& + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
[825]	287
[7646]	288	! Coriolis at T points (m*f)
	289	zmf(ji,jj) = zm1 * ff_t(ji,jj)
[825]	290
[7646]	291	! m/dt
	292	zmU_t(ji,jj) = zmassU * z1_dtevp
	293	zmV_t(ji,jj) = zmassV * z1_dtevp
[825]	294
[7646]	295	! Drag ice-atm.
	296	zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
	297	zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
[825]	298
[7646]	299	! Surface pressure gradient (- mgGRAD(ssh)) at U-V points
	300	zspgU(ji,jj) = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
	301	zspgV(ji,jj) = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
[825]	302
[7646]	303	! masks
	304	zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice
	305	zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice
[834]	306
[7646]	307	! switches
	308	zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
	309	zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
[825]	310
	311	END DO
	312	END DO
[7646]	313	CALL lbc_lnk( zmf, 'T', 1. )
[921]	314	!
	315	!------------------------------------------------------------------------------!
	316	! 3) Solution of the momentum equation, iterative procedure
	317	!------------------------------------------------------------------------------!
	318	!
[2528]	319	! !----------------------!
[5123]	320	DO jter = 1 , nn_nevp ! loop over jter !
[2528]	321	! !----------------------!
[7646]	322	IF(ln_ctl) THEN ! Convergence test
	323	DO jj = 1, jpjm1
[7753]	324	zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
	325	zv_ice(:,jj) = v_ice(:,jj)
[7646]	326	END DO
	327	ENDIF
[825]	328
[7646]	329	! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
	330	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
	331	DO ji = 1, jpim1
[825]	332
[7646]	333	! shear at F points
[5123]	334	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) &
	335	& + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
[7646]	336	& ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
[825]	337
[7646]	338	END DO
	339	END DO
	340	CALL lbc_lnk( zds, 'F', 1. )
[825]	341
[7646]	342	DO jj = 2, jpjm1
	343	DO ji = 2, jpim1 ! no vector loop
[825]	344
[7646]	345	! shear**2 at T points (doc eq. A16)
	346	zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) &
	347	& + 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) &
	348	& ) * 0.25_wp * r1_e1e2t(ji,jj)
	349
	350	! divergence at T points
	351	zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
	352	& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
	353	& ) * r1_e1e2t(ji,jj)
	354	zdiv2 = zdiv * zdiv
	355
	356	! tension at T points
	357	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) &
	358	& - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
	359	& ) * r1_e1e2t(ji,jj)
	360	zdt2 = zdt * zdt
	361
	362	! delta at T points
	363	zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )
	364
	365	! P/delta at T points
	366	zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
	367
	368	! stress at T points
	369	zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
	370	zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
	371
[921]	372	END DO
	373	END DO
[7646]	374	CALL lbc_lnk( zp_delt, 'T', 1. )
[921]	375
[7646]	376	DO jj = 1, jpjm1
	377	DO ji = 1, jpim1
[825]	378
[7646]	379	! P/delta at F points
	380	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) )
	381
	382	! stress at F points
	383	zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
[825]	384
[7646]	385	END DO
	386	END DO
	387	CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
	388
[4205]	389
[7646]	390	! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
	391	DO jj = 2, jpjm1
	392	DO ji = fs_2, fs_jpim1
[825]	393
[7646]	394	! U points
	395	zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) &
	396	& + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) &
	397	& ) * r1_e2u(ji,jj) &
	398	& + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) &
	399	& ) * 2._wp * r1_e1u(ji,jj) &
	400	& ) * r1_e1e2u(ji,jj)
[825]	401
[7646]	402	! V points
	403	zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) &
	404	& - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) &
	405	& ) * r1_e1v(ji,jj) &
	406	& + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) &
	407	& ) * 2._wp * r1_e2v(ji,jj) &
	408	& ) * r1_e1e2v(ji,jj)
[825]	409
[7646]	410	! u_ice at V point
	411	u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) &
	412	& + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
	413
	414	! v_ice at U point
	415	v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) &
	416	& + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
[4205]	417
[5123]	418	END DO
	419	END DO
[825]	420	!
[7646]	421	! --- Computation of ice velocity --- !
	422	! Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp
	423	! Bouillon et al. 2009 (eq 34-35) => stable
	424	IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations
	425
	426	DO jj = 2, jpjm1
[921]	427	DO ji = fs_2, fs_jpim1
[825]	428
[7646]	429	! tau_io/(v_oce - v_ice)
	430	zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) &
	431	& + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
	432
[8239]	433	! Ocean-to-Ice stress
	434	ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
	435
[7646]	436	! tau_bottom/v_ice
	437	zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
	438	zTauB = - tau_icebfr(ji,jj) / zvel
	439
	440	! Coriolis at V-points (energy conserving formulation)
[8239]	441	zCory(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * &
[7646]	442	& ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) &
	443	& + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
	444
	445	! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
[8239]	446	zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
[7646]	447
	448	! landfast switch => 0 = static friction ; 1 = sliding friction
	449	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
	450
	451	! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
	452	v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity
	453	& + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	454	& ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	455	& + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
	456	& ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
	457	& ) * zmaskV(ji,jj)
	458	! Bouillon 2013
	459	!!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) ) &
[8239]	460	!! & + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj) &
[7646]	461	!! & ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
	462
[921]	463	END DO
	464	END DO
[7646]	465	CALL lbc_lnk( v_ice, 'V', -1. )
	466
	467	#if defined key_agrif
	468	!! CALL agrif_interp_lim3( 'V', jter, nn_nevp )
	469	CALL agrif_interp_lim3( 'V' )
[3680]	470	#endif
[7646]	471	IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'V' )
[825]	472
[7646]	473	DO jj = 2, jpjm1
[921]	474	DO ji = fs_2, fs_jpim1
[7646]	475
	476	! tau_io/(u_oce - u_ice)
	477	zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) &
	478	& + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
[834]	479
[8239]	480	! Ocean-to-Ice stress
	481	ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
	482
[7646]	483	! tau_bottom/u_ice
	484	zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
	485	zTauB = - tau_icebfr(ji,jj) / zvel
	486
	487	! Coriolis at U-points (energy conserving formulation)
[8239]	488	zCorx(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * &
[7646]	489	& ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) &
	490	& + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
	491
	492	! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
[8239]	493	zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
[7646]	494
	495	! landfast switch => 0 = static friction ; 1 = sliding friction
	496	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
	497
	498	! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
	499	u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity
	500	& + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	501	& ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	502	& + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
	503	& ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
	504	& ) * zmaskU(ji,jj)
	505	! Bouillon 2013
	506	!!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) ) &
[8239]	507	!! & + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj) &
[7646]	508	!! & ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
[921]	509	END DO
	510	END DO
[7646]	511	CALL lbc_lnk( u_ice, 'U', -1. )
	512
	513	#if defined key_agrif
	514	!! CALL agrif_interp_lim3( 'U', jter, nn_nevp )
	515	CALL agrif_interp_lim3( 'U' )
[3680]	516	#endif
[7646]	517	IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'U' )
[825]	518
[7646]	519	ELSE ! odd iterations
	520
	521	DO jj = 2, jpjm1
[921]	522	DO ji = fs_2, fs_jpim1
[825]	523
[7646]	524	! tau_io/(u_oce - u_ice)
	525	zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) &
	526	& + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
	527
[8239]	528	! Ocean-to-Ice stress
	529	ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
	530
[7646]	531	! tau_bottom/u_ice
	532	zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
	533	zTauB = - tau_icebfr(ji,jj) / zvel
	534
	535	! Coriolis at U-points (energy conserving formulation)
[8239]	536	zCorx(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * &
[7646]	537	& ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) &
	538	& + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
	539
	540	! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
[8239]	541	zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
[7646]	542
	543	! landfast switch => 0 = static friction ; 1 = sliding friction
	544	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
	545
	546	! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
	547	u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity
	548	& + zTauE + 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	& ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
	552	& ) * zmaskU(ji,jj)
	553	! Bouillon 2013
	554	!!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) ) &
[8239]	555	!! & + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj) &
[7646]	556	!! & ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
[921]	557	END DO
	558	END DO
[7646]	559	CALL lbc_lnk( u_ice, 'U', -1. )
	560
	561	#if defined key_agrif
	562	!! CALL agrif_interp_lim3( 'U', jter, nn_nevp )
	563	CALL agrif_interp_lim3( 'U' )
[3680]	564	#endif
[7646]	565	IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'U' )
[825]	566
[7646]	567	DO jj = 2, jpjm1
[921]	568	DO ji = fs_2, fs_jpim1
[7646]	569
	570	! tau_io/(v_oce - v_ice)
	571	zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) &
	572	& + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
[825]	573
[8239]	574	! Ocean-to-Ice stress
	575	ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
	576
[7646]	577	! tau_bottom/v_ice
	578	zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
	579	ztauB = - tau_icebfr(ji,jj) / zvel
	580
	581	! Coriolis at V-points (energy conserving formulation)
[8239]	582	zCory(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * &
[7646]	583	& ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) &
	584	& + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
	585
	586	! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
[8239]	587	zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
[7646]	588
	589	! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction
	590	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
	591
	592	! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
	593	v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity
	594	& + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	595	& ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	596	& + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
	597	& ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
	598	& ) * zmaskV(ji,jj)
	599	! Bouillon 2013
	600	!!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) ) &
[8239]	601	!! & + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj) &
[7646]	602	!! & ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
[5123]	603	END DO
	604	END DO
[7646]	605	CALL lbc_lnk( v_ice, 'V', -1. )
	606
	607	#if defined key_agrif
	608	!! CALL agrif_interp_lim3( 'V', jter, nn_nevp )
	609	CALL agrif_interp_lim3( 'V' )
[3680]	610	#endif
[7646]	611	IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'V' )
[825]	612
[921]	613	ENDIF
[4161]	614
[7646]	615	IF(ln_ctl) THEN ! Convergence test
	616	DO jj = 2 , jpjm1
[7753]	617	zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
[921]	618	END DO
[7646]	619	zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
[921]	620	IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain
	621	ENDIF
[7646]	622	!
[4161]	623	! ! ==================== !
[868]	624	END DO ! end loop over jter !
[825]	625	! ! ==================== !
[921]	626	!
	627	!------------------------------------------------------------------------------!
[7646]	628	! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
[921]	629	!------------------------------------------------------------------------------!
[7646]	630	DO jj = 1, jpjm1
	631	DO ji = 1, jpim1
[866]	632
[7646]	633	! shear at F points
	634	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) &
	635	& + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
	636	& ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
[5429]	637
[868]	638	END DO
[7646]	639	END DO
	640	CALL lbc_lnk( zds, 'F', 1. )
	641
	642	DO jj = 2, jpjm1
	643	DO ji = 2, jpim1 ! no vector loop
	644
	645	! tension**2 at T points
	646	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) &
	647	& - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
	648	& ) * r1_e1e2t(ji,jj)
	649	zdt2 = zdt * zdt
	650
	651	! shear**2 at T points (doc eq. A16)
	652	zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) &
	653	& + 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) &
	654	& ) * 0.25_wp * r1_e1e2t(ji,jj)
	655
	656	! shear at T points
	657	shear_i(ji,jj) = SQRT( zdt2 + zds2 )
[868]	658
[7646]	659	! divergence at T points
	660	divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
	661	& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
	662	& ) * r1_e1e2t(ji,jj)
	663
	664	! delta at T points
	665	zdelta = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )
	666	rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
	667	delta_i(ji,jj) = zdelta + rn_creepl * rswitch
[868]	668
[5123]	669	END DO
	670	END DO
[7646]	671	CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. )
	672
	673	! --- Store the stress tensor for the next time step --- !
[7753]	674	stress1_i (:,:) = zs1 (:,:)
	675	stress2_i (:,:) = zs2 (:,:)
	676	stress12_i(:,:) = zs12(:,:)
[7646]	677	!
[868]	678
[921]	679	!------------------------------------------------------------------------------!
[8239]	680	! 5) SIMIP diagnostics
[921]	681	!------------------------------------------------------------------------------!
[8239]	682
	683	DO jj = 2, jpjm1
	684	DO ji = 2, jpim1
	685	rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
	686
	687	! Stress tensor invariants (normal and shear stress N/m)
	688	diag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch ! normal stress
	689	diag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )*2 + 4zs12(ji,jj)*2 ) rswitch ! shear stress
	690
	691	! Stress terms of the momentum equation (N/m2)
	692	diag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch ! sea surface slope stress term
	693	diag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch
	694
	695	diag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch ! Coriolis stress term
	696	diag_corstry(ji,jj) = zCory(ji,jj) * rswitch
	697
	698	diag_intstrx(ji,jj) = zfU(ji,jj) * rswitch ! internal stress term
	699	diag_intstry(ji,jj) = zfV(ji,jj) * rswitch
	700
	701	diag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch ! oceanic stress
	702	diag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch
	703
	704	! 2D ice mass, snow mass, area transport arrays (X, Y)
	705	zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch
	706	zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch
	707
	708	diag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
	709	diag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) ! '' Y- ''
	710
	711	diag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
	712	diag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) ! '' Y- ''
	713
	714	diag_xatrp(ji,jj) = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) ) ! area transport, X-component
	715	diag_yatrp(ji,jj) = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) ) ! '' Y- ''
	716
	717	END DO
	718	END DO
	719
	720	CALL lbc_lnk_multi( diag_sig1 , 'T', 1., diag_sig2 , 'T', 1., &
	721	& diag_dssh_dx, 'U', -1., diag_dssh_dy, 'V', -1., &
	722	& diag_corstrx, 'U', -1., diag_corstry, 'V', -1., &
	723	& diag_intstrx, 'U', -1., diag_intstry, 'V', -1. )
	724
	725	CALL lbc_lnk_multi( diag_utau_oi, 'U', -1., diag_vtau_oi, 'V', -1. )
	726
[8289]	727	CALL lbc_lnk_multi( diag_xmtrp_ice, 'U', -1., diag_xmtrp_snw, 'U', -1., &
	728	& diag_xatrp , 'U', -1., diag_ymtrp_ice, 'V', -1., &
	729	& diag_ymtrp_snw, 'V', -1., diag_yatrp , 'V', -1. )
	730
[921]	731	!
[8239]	732	!------------------------------------------------------------------------------!
	733	! 6) Control prints of residual and charge ellipse
	734	!------------------------------------------------------------------------------!
	735	!
[834]	736	! print the residual for convergence
	737	IF(ln_ctl) THEN
[868]	738	WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, jter
[834]	739	CALL prt_ctl_info(charout)
	740	CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
	741	ENDIF
[825]	742
[834]	743	! print charge ellipse
	744	! This can be desactivated once the user is sure that the stress state
	745	! lie on the charge ellipse. See Bouillon et al. 08 for more details
[825]	746	IF(ln_ctl) THEN
[8319]	747	IF( MOD(kt_ice+nn_fsbc-1,nwrite) == 0 ) THEN
	748	WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 1000, kt_ice, 0, 0, ' ch. ell. '
[825]	749	CALL prt_ctl_info(charout)
[7646]	750	DO jj = 2, jpjm1
[825]	751	DO ji = 2, jpim1
[7646]	752	IF (strength(ji,jj) > 1.0) THEN
	753	zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)*2 + 4zs12(ji,jj)*2 ) ) / ( 2strength(ji,jj) )
	754	zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)*2 + 4zs12(ji,jj)*2 ) ) / ( 2strength(ji,jj) )
[825]	755	WRITE(charout,FMT="('lim_rhg :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
	756	CALL prt_ctl_info(charout)
	757	ENDIF
	758	END DO
	759	END DO
[8319]	760	WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 2000, kt_ice, 0, 0, ' ch. ell. '
[825]	761	CALL prt_ctl_info(charout)
	762	ENDIF
	763	ENDIF
[7646]	764	!
	765
	766	CALL wrk_dealloc( jpi,jpj, z1_e1t0, z1_e2t0, zp_delt )
	767	CALL wrk_dealloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
	768	CALL wrk_dealloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
	769	CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
	770	CALL wrk_dealloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
[8239]	771	CALL wrk_dealloc( jpi,jpj, zCorx, zCory )
	772	CALL wrk_dealloc( jpi,jpj, ztaux_oi, ztauy_oi )
[3294]	773
[825]	774	END SUBROUTINE lim_rhg
	775
	776	#else
	777	!!----------------------------------------------------------------------
	778	!! Default option Dummy module NO LIM sea-ice model
	779	!!----------------------------------------------------------------------
	780	CONTAINS
[7646]	781	SUBROUTINE lim_rhg ! Dummy routine
	782	WRITE(,) 'lim_rhg: You should not have seen this print! error?'
[825]	783	END SUBROUTINE lim_rhg
	784	#endif
	785
	786	!!==============================================================================
	787	END MODULE limrhg

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

Download in other formats: