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/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 4269

Last change on this file since 4269 was 4269, checked in by clem, 10 years ago
update lim3 for bdy purpose
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
[2528]	9	!! 3.3 ! 2009-05 (G.Garric) addition of the lim2_evp cas
[3680]	10	!! 3.4 ! 2011-01 (A. Porter) dynamical allocation
	11	!! 3.5 ! 2012-08 (R. Benshila) AGRIF
[1244]	12	!!----------------------------------------------------------------------
[2528]	13	#if defined key_lim3 \|\| ( defined key_lim2 && ! defined key_lim2_vp )
[825]	14	!!----------------------------------------------------------------------
[2528]	15	!! 'key_lim3' OR LIM-3 sea-ice model
[2717]	16	!! 'key_lim2' AND NOT 'key_lim2_vp' EVP LIM-2 sea-ice model
[825]	17	!!----------------------------------------------------------------------
[3625]	18	!! lim_rhg : computes ice velocities
[825]	19	!!----------------------------------------------------------------------
[3625]	20	USE phycst ! Physical constant
	21	USE oce , ONLY : snwice_mass, snwice_mass_b
	22	USE par_oce ! Ocean parameters
	23	USE dom_oce ! Ocean domain
	24	USE sbc_oce ! Surface boundary condition: ocean fields
	25	USE sbc_ice ! Surface boundary condition: ice fields
[2528]	26	#if defined key_lim3
[3625]	27	USE ice ! LIM-3: ice variables
	28	USE dom_ice ! LIM-3: ice domain
	29	USE limitd_me ! LIM-3:
[2528]	30	#else
[3625]	31	USE ice_2 ! LIM-2: ice variables
	32	USE dom_ice_2 ! LIM-2: ice domain
[2528]	33	#endif
[3625]	34	USE lbclnk ! Lateral Boundary Condition / MPP link
	35	USE lib_mpp ! MPP library
	36	USE wrk_nemo ! work arrays
	37	USE in_out_manager ! I/O manager
	38	USE prtctl ! Print control
	39	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[3680]	40	#if defined key_agrif && defined key_lim2
	41	USE agrif_lim2_interp
	42	#endif
[4161]	43	#if defined key_bdy
	44	USE bdyice_lim
	45	#endif
[825]	46
	47	IMPLICIT NONE
	48	PRIVATE
	49
[2715]	50	PUBLIC lim_rhg ! routine called by lim_dyn (or lim_dyn_2)
[825]	51
[2528]	52	REAL(wp) :: rzero = 0._wp ! constant values
	53	REAL(wp) :: rone = 1._wp ! constant values
	54
[868]	55	!! * Substitutions
	56	# include "vectopt_loop_substitute.h90"
[825]	57	!!----------------------------------------------------------------------
[4161]	58	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]	59	!! $Id$
[2528]	60	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	61	!!----------------------------------------------------------------------
	62	CONTAINS
	63
	64	SUBROUTINE lim_rhg( k_j1, k_jpj )
	65	!!-------------------------------------------------------------------
[834]	66	!! * SUBROUTINE lim_rhg *
	67	!! EVP-C-grid
[825]	68	!!
[834]	69	!! ** purpose : determines sea ice drift from wind stress, ice-ocean
[825]	70	!! stress and sea-surface slope. Ice-ice interaction is described by
[834]	71	!! a non-linear elasto-viscous-plastic (EVP) law including shear
	72	!! strength and a bulk rheology (Hunke and Dukowicz, 2002).
[825]	73	!!
[834]	74	!! The points in the C-grid look like this, dear reader
[825]	75	!!
[834]	76	!! (ji,jj)
	77	!! \|
	78	!! \|
	79	!! (ji-1,jj) \| (ji,jj)
	80	!! ---------
	81	!! \| \|
	82	!! \| (ji,jj) \|------(ji,jj)
	83	!! \| \|
	84	!! ---------
	85	!! (ji-1,jj-1) (ji,jj-1)
[825]	86	!!
[834]	87	!! ** Inputs : - wind forcing (stress), oceanic currents
	88	!! ice total volume (vt_i) per unit area
	89	!! snow total volume (vt_s) per unit area
[825]	90	!!
[834]	91	!! ** Action : - compute u_ice, v_ice : the components of the
	92	!! sea-ice velocity vector
	93	!! - compute delta_i, shear_i, divu_i, which are inputs
	94	!! of the ice thickness distribution
[825]	95	!!
[834]	96	!! ** Steps : 1) Compute ice snow mass, ice strength
	97	!! 2) Compute wind, oceanic stresses, mass terms and
	98	!! coriolis terms of the momentum equation
	99	!! 3) Solve the momentum equation (iterative procedure)
	100	!! 4) Prevent high velocities if the ice is thin
	101	!! 5) Recompute invariants of the strain rate tensor
	102	!! which are inputs of the ITD, store stress
	103	!! for the next time step
	104	!! 6) Control prints of residual (convergence)
	105	!! and charge ellipse.
	106	!! The user should make sure that the parameters
	107	!! nevp, telast and creepl maintain stress state
	108	!! on the charge ellipse for plastic flow
	109	!! e.g. in the Canadian Archipelago
	110	!!
[2528]	111	!! References : Hunke and Dukowicz, JPO97
	112	!! Bouillon et al., Ocean Modelling 2009
	113	!! Vancoppenolle et al., Ocean Modelling 2008
	114	!!-------------------------------------------------------------------
	115	INTEGER, INTENT(in) :: k_j1 ! southern j-index for ice computation
	116	INTEGER, INTENT(in) :: k_jpj ! northern j-index for ice computation
[834]	117	!!
[2528]	118	INTEGER :: ji, jj ! dummy loop indices
	119	INTEGER :: jter ! local integers
[825]	120	CHARACTER (len=50) :: charout
[2528]	121	REAL(wp) :: zt11, zt12, zt21, zt22, ztagnx, ztagny, delta !
	122	REAL(wp) :: za, zstms, zsang, zmask ! local scalars
[825]	123
[2715]	124	REAL(wp) :: dtevp ! time step for subcycling
[4205]	125	REAL(wp) :: dtotel, ecc2, ecci ! square of yield ellipse eccenticity
[2715]	126	REAL(wp) :: z0, zr, zcca, zccb ! temporary scalars
	127	REAL(wp) :: zu_ice2, zv_ice1 !
[3791]	128	REAL(wp) :: zddc, zdtc, zzdst ! delta on corners and on centre
[2715]	129	REAL(wp) :: zdsshx, zdsshy ! term for the gradient of ocean surface
	130	REAL(wp) :: sigma1, sigma2 ! internal ice stress
[825]	131
[2715]	132	REAL(wp) :: zresm ! Maximal error on ice velocity
	133	REAL(wp) :: zindb ! ice (1) or not (0)
	134	REAL(wp) :: zdummy ! dummy argument
[3625]	135	REAL(wp) :: zintb, zintn ! dummy argument
[3294]	136
	137	REAL(wp), POINTER, DIMENSION(:,:) :: zpresh ! temporary array for ice strength
	138	REAL(wp), POINTER, DIMENSION(:,:) :: zpreshc ! Ice strength on grid cell corners (zpreshc)
	139	REAL(wp), POINTER, DIMENSION(:,:) :: zfrld1, zfrld2 ! lead fraction on U/V points
	140	REAL(wp), POINTER, DIMENSION(:,:) :: zmass1, zmass2 ! ice/snow mass on U/V points
	141	REAL(wp), POINTER, DIMENSION(:,:) :: zcorl1, zcorl2 ! coriolis parameter on U/V points
	142	REAL(wp), POINTER, DIMENSION(:,:) :: za1ct , za2ct ! temporary arrays
	143	REAL(wp), POINTER, DIMENSION(:,:) :: zc1 ! ice mass
	144	REAL(wp), POINTER, DIMENSION(:,:) :: zusw ! temporary weight for ice strength computation
	145	REAL(wp), POINTER, DIMENSION(:,:) :: u_oce1, v_oce1 ! ocean u/v component on U points
	146	REAL(wp), POINTER, DIMENSION(:,:) :: u_oce2, v_oce2 ! ocean u/v component on V points
	147	REAL(wp), POINTER, DIMENSION(:,:) :: u_ice2, v_ice1 ! ice u/v component on V/U point
	148	REAL(wp), POINTER, DIMENSION(:,:) :: zf1 , zf2 ! arrays for internal stresses
	149
	150	REAL(wp), POINTER, DIMENSION(:,:) :: zdd , zdt ! Divergence and tension at centre of grid cells
	151	REAL(wp), POINTER, DIMENSION(:,:) :: zds ! Shear on northeast corner of grid cells
[3791]	152	REAL(wp), POINTER, DIMENSION(:,:) :: zdst ! Shear on centre of grid cells
[3294]	153	REAL(wp), POINTER, DIMENSION(:,:) :: deltat, deltac ! Delta at centre and corners of grid cells
	154	REAL(wp), POINTER, DIMENSION(:,:) :: zs1 , zs2 ! Diagonal stress tensor components zs1 and zs2
	155	REAL(wp), POINTER, DIMENSION(:,:) :: zs12 ! Non-diagonal stress tensor component zs12
	156	REAL(wp), POINTER, DIMENSION(:,:) :: zu_ice, zv_ice, zresr ! Local error on velocity
[3625]	157	REAL(wp), POINTER, DIMENSION(:,:) :: zpice ! array used for the calculation of ice surface slope:
	158	! ocean surface (ssh_m) if ice is not embedded
[3791]	159	! ice top surface if ice is embedded
[2528]	160	!!-------------------------------------------------------------------
[3294]	161
	162	CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
	163	CALL wrk_alloc( jpi,jpj, zc1 , u_oce1, u_oce2, u_ice2, zusw , v_oce1 , v_oce2, v_ice1 )
[3791]	164	CALL wrk_alloc( jpi,jpj, zf1 , deltat, zu_ice, zf2 , deltac, zv_ice , zdd , zdt , zds , zdst )
[3625]	165	CALL wrk_alloc( jpi,jpj, zdd , zdt , zds , zs1 , zs2 , zs12 , zresr , zpice )
[3294]	166
[2528]	167	#if defined key_lim2 && ! defined key_lim2_vp
	168	# if defined key_agrif
	169	USE ice_2, vt_s => hsnm
	170	USE ice_2, vt_i => hicm
	171	# else
	172	vt_s => hsnm
	173	vt_i => hicm
	174	# endif
	175	at_i(:,:) = 1. - frld(:,:)
	176	#endif
[3680]	177	#if defined key_agrif && defined key_lim2
	178	CALL agrif_rhg_lim2_load ! First interpolation of coarse values
	179	#endif
[921]	180	!
	181	!------------------------------------------------------------------------------!
	182	! 1) Ice-Snow mass (zc1), ice strength (zpresh) !
	183	!------------------------------------------------------------------------------!
	184	!
[825]	185	! Put every vector to 0
[2528]	186	zpresh (:,:) = 0._wp ; zc1 (:,:) = 0._wp
	187	zpreshc(:,:) = 0._wp
	188	u_ice2 (:,:) = 0._wp ; v_ice1(:,:) = 0._wp
	189	zdd (:,:) = 0._wp ; zdt (:,:) = 0._wp ; zds(:,:) = 0._wp
[825]	190
[2528]	191	#if defined key_lim3
	192	CALL lim_itd_me_icestrength( ridge_scheme_swi ) ! LIM-3: Ice strength on T-points
	193	#endif
[825]	194
[868]	195	!CDIR NOVERRCHK
[2528]	196	DO jj = k_j1 , k_jpj ! Ice mass and temp variables
[868]	197	!CDIR NOVERRCHK
[825]	198	DO ji = 1 , jpi
	199	zc1(ji,jj) = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) )
[2528]	200	#if defined key_lim3
[4205]	201	zpresh(ji,jj) = tms(ji,jj) * strength(ji,jj)
[2528]	202	#endif
[2580]	203	#if defined key_lim2
	204	zpresh(ji,jj) = tms(ji,jj) * pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) )
	205	#endif
[866]	206	! tmi = 1 where there is ice or on land
[2528]	207	tmi(ji,jj) = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - epsd ) ) ) * tms(ji,jj)
[825]	208	END DO
	209	END DO
	210
[834]	211	! Ice strength on grid cell corners (zpreshc)
	212	! needed for calculation of shear stress
[868]	213	!CDIR NOVERRCHK
[825]	214	DO jj = k_j1+1, k_jpj-1
[868]	215	!CDIR NOVERRCHK
	216	DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1)
[921]	217	zstms = tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
	218	& tms(ji,jj+1) * wght(ji+1,jj+1,1,2) + &
	219	& tms(ji+1,jj) * wght(ji+1,jj+1,2,1) + &
	220	& tms(ji,jj) * wght(ji+1,jj+1,1,1)
	221	zusw(ji,jj) = 1.0 / MAX( zstms, epsd )
	222	zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
	223	& zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + &
	224	& zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + &
	225	& zpresh(ji,jj) * wght(ji+1,jj+1,1,1) &
	226	& ) * zusw(ji,jj)
[825]	227	END DO
	228	END DO
	229
	230	CALL lbc_lnk( zpreshc(:,:), 'F', 1. )
[921]	231	!
	232	!------------------------------------------------------------------------------!
	233	! 2) Wind / ocean stress, mass terms, coriolis terms
	234	!------------------------------------------------------------------------------!
	235	!
[825]	236	! Wind stress, coriolis and mass terms on the sides of the squares
	237	! zfrld1: lead fraction on U-points
	238	! zfrld2: lead fraction on V-points
	239	! zmass1: ice/snow mass on U-points
	240	! zmass2: ice/snow mass on V-points
	241	! zcorl1: Coriolis parameter on U-points
	242	! zcorl2: Coriolis parameter on V-points
	243	! (ztagnx,ztagny): wind stress on U/V points
	244	! u_oce1: ocean u component on u points
	245	! v_oce1: ocean v component on u points
	246	! u_oce2: ocean u component on v points
	247	! v_oce2: ocean v component on v points
[921]	248
[3625]	249	IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==!
	250	!
	251	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
	252	! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
	253	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
	254	!
	255	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
	256	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
	257	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
	258	!
	259	zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
	260	!
	261	ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==!
	262	zpice(:,:) = ssh_m(:,:)
	263	ENDIF
	264
[825]	265	DO jj = k_j1+1, k_jpj-1
[868]	266	DO ji = fs_2, fs_jpim1
[825]	267
[2528]	268	zt11 = tms(ji ,jj) * e1t(ji ,jj)
	269	zt12 = tms(ji+1,jj) * e1t(ji+1,jj)
	270	zt21 = tms(ji,jj ) * e2t(ji,jj )
	271	zt22 = tms(ji,jj+1) * e2t(ji,jj+1)
[825]	272
	273	! Leads area.
[2528]	274	zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd )
	275	zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd )
[825]	276
	277	! Mass, coriolis coeff. and currents
[834]	278	zmass1(ji,jj) = ( zt12zc1(ji,jj) + zt11zc1(ji+1,jj) ) / (zt11+zt12+epsd)
	279	zmass2(ji,jj) = ( zt22zc1(ji,jj) + zt21zc1(ji,jj+1) ) / (zt21+zt22+epsd)
[2528]	280	zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)fcor(ji,jj) + e1t(ji,jj)fcor(ji+1,jj) ) &
	281	& / ( e1t(ji,jj) + e1t(ji+1,jj) + epsd )
	282	zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)fcor(ji,jj) + e2t(ji,jj)fcor(ji,jj+1) ) &
	283	& / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd )
[825]	284	!
[888]	285	u_oce1(ji,jj) = u_oce(ji,jj)
	286	v_oce2(ji,jj) = v_oce(ji,jj)
[825]	287
[834]	288	! Ocean has no slip boundary condition
[888]	289	v_oce1(ji,jj) = 0.5( (v_oce(ji,jj)+v_oce(ji,jj-1))e1t(ji,jj) &
[921]	290	& +(v_oce(ji+1,jj)+v_oce(ji+1,jj-1))*e1t(ji+1,jj)) &
	291	& /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
[825]	292
[888]	293	u_oce2(ji,jj) = 0.5((u_oce(ji,jj)+u_oce(ji-1,jj))e2t(ji,jj) &
[921]	294	& +(u_oce(ji,jj+1)+u_oce(ji-1,jj+1))*e2t(ji,jj+1)) &
	295	& / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
[825]	296
[1469]	297	! Wind stress at U,V-point
	298	ztagnx = ( 1. - zfrld1(ji,jj) ) * utau_ice(ji,jj)
	299	ztagny = ( 1. - zfrld2(ji,jj) ) * vtau_ice(ji,jj)
[825]	300
[834]	301	! Computation of the velocity field taking into account the ice internal interaction.
[825]	302	! Terms that are independent of the velocity field.
	303
	304	! SB On utilise maintenant le gradient de la pente de l'ocean
	305	! include it later
[834]	306
[3625]	307	zdsshx = ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj)
	308	zdsshy = ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj)
[825]	309
	310	za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx
	311	za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy
	312
	313	END DO
	314	END DO
	315
[921]	316	!
	317	!------------------------------------------------------------------------------!
	318	! 3) Solution of the momentum equation, iterative procedure
	319	!------------------------------------------------------------------------------!
	320	!
[825]	321	! Time step for subcycling
	322	dtevp = rdt_ice / nevp
[2528]	323	dtotel = dtevp / ( 2._wp * telast )
[825]	324
	325	!-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)
[2528]	326	ecc2 = ecc * ecc
[4205]	327	ecci = 1. / ecc2
[825]	328
	329	!-Initialise stress tensor
[2528]	330	zs1 (:,:) = stress1_i (:,:)
	331	zs2 (:,:) = stress2_i (:,:)
[866]	332	zs12(:,:) = stress12_i(:,:)
[825]	333
[2528]	334	! !----------------------!
[868]	335	DO jter = 1 , nevp ! loop over jter !
[2528]	336	! !----------------------!
[825]	337	DO jj = k_j1, k_jpj-1
[2528]	338	zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
[825]	339	zv_ice(:,jj) = v_ice(:,jj)
[921]	340	END DO
[825]	341
[834]	342	DO jj = k_j1+1, k_jpj-1
[988]	343	DO ji = fs_2, jpim1 !RB bug no vect opt due to tmi
[825]	344
[921]	345	!
	346	!- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002)
	347	!- zdd(:,:), zdt(:,:): divergence and tension at centre of grid cells
	348	!- zds(:,:): shear on northeast corner of grid cells
	349	!
	350	!- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded,
	351	! there are many repeated calculations.
	352	! Speed could be improved by regrouping terms. For
	353	! the moment, however, the stress is on clarity of coding to avoid
	354	! bugs (Martin, for Miguel).
	355	!
	356	!- ALSO: arrays zdd, zdt, zds and delta could
	357	! be removed in the future to minimise memory demand.
	358	!
	359	!- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of
	360	! grid cells, exactly as in the B grid case. For simplicity, the indexation on
	361	! the corners is the same as in the B grid.
	362	!
	363	!
	364	zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) &
	365	& -e2u(ji-1,jj)*u_ice(ji-1,jj) &
	366	& +e1v(ji,jj)*v_ice(ji,jj) &
	367	& -e1v(ji,jj-1)*v_ice(ji,jj-1) &
	368	& ) &
	369	& / area(ji,jj)
[825]	370
[921]	371	zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) &
	372	& -u_ice(ji-1,jj)/e2u(ji-1,jj) &
	373	& )e2t(ji,jj)e2t(ji,jj) &
	374	& -( v_ice(ji,jj)/e1v(ji,jj) &
	375	& -v_ice(ji,jj-1)/e1v(ji,jj-1) &
	376	& )e1t(ji,jj)e1t(ji,jj) &
	377	& ) &
	378	& / area(ji,jj)
[825]	379
[921]	380	!
	381	zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1) &
	382	& -u_ice(ji,jj)/e1u(ji,jj) &
	383	& )e1f(ji,jj)e1f(ji,jj) &
	384	& +( v_ice(ji+1,jj)/e2v(ji+1,jj) &
	385	& -v_ice(ji,jj)/e2v(ji,jj) &
	386	& )e2f(ji,jj)e2f(ji,jj) &
	387	& ) &
	388	& / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) &
	389	& * tmi(ji,jj) * tmi(ji,jj+1) &
	390	& * tmi(ji+1,jj) * tmi(ji+1,jj+1)
[825]	391
	392
[921]	393	v_ice1(ji,jj) = 0.5( (v_ice(ji,jj)+v_ice(ji,jj-1))e1t(ji+1,jj) &
	394	& +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
	395	& /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
[825]	396
[921]	397	u_ice2(ji,jj) = 0.5( (u_ice(ji,jj)+u_ice(ji-1,jj))e2t(ji,jj+1) &
	398	& +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &
	399	& /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
[825]	400
[921]	401	END DO
	402	END DO
[2528]	403	CALL lbc_lnk( v_ice1, 'U', -1. ) ; CALL lbc_lnk( u_ice2, 'V', -1. ) ! lateral boundary cond.
[921]	404
[868]	405	!CDIR NOVERRCHK
[921]	406	DO jj = k_j1+1, k_jpj-1
[868]	407	!CDIR NOVERRCHK
[921]	408	DO ji = fs_2, fs_jpim1
[825]	409
[921]	410	!- Calculate Delta at centre of grid cells
[3791]	411	zzdst = ( e2u(ji , jj) * v_ice1(ji ,jj) &
[2528]	412	& - e2u(ji-1, jj) * v_ice1(ji-1,jj) &
	413	& + e1v(ji, jj ) * u_ice2(ji,jj ) &
	414	& - e1v(ji, jj-1) * u_ice2(ji,jj-1) &
	415	& ) &
[921]	416	& / area(ji,jj)
[825]	417
[3791]	418	delta = SQRT( zdd(ji,jj)zdd(ji,jj) + ( zdt(ji,jj)zdt(ji,jj) + zzdstzzdst ) usecc2 )
[4161]	419	! MV rewriting
	420	! deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)2 + (zdt(ji,jj)2 + zzdst*2)usecc2), creepl )
	421	!!gm faster to replace the line above with simply:
	422	!! deltat(ji,jj) = MAX( delta, creepl )
	423	!!gm end
	424	deltat(ji,jj) = delta + creepl
	425	! END MV
[921]	426	!-Calculate stress tensor components zs1 and zs2
	427	!-at centre of grid cells (see section 3.5 of CICE user's guide).
[4269]	428	!zs1(ji,jj) = ( zs1(ji,jj) - dtotel( ( 1._wp - alphaevp) zs1(ji,jj) + &
	429	! & ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) ) &
	430	! & / ( 1._wp + alphaevp * dtotel )
[825]	431
[4269]	432	!zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) - &
	433	! zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) ) &
	434	! & / ( 1._wp + alphaevp * ecc2 * dtotel )
[825]	435
[4205]	436	! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp)
[4269]	437	zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) ) &
	438	& * zpresh(ji,jj) ) ) / ( 1._wp + dtotel )
	439	zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) ) &
	440	& / ( 1._wp + dtotel )
[4205]	441
[921]	442	END DO
	443	END DO
[825]	444
[921]	445	CALL lbc_lnk( zs1(:,:), 'T', 1. )
	446	CALL lbc_lnk( zs2(:,:), 'T', 1. )
[825]	447
[868]	448	!CDIR NOVERRCHK
[921]	449	DO jj = k_j1+1, k_jpj-1
[868]	450	!CDIR NOVERRCHK
[921]	451	DO ji = fs_2, fs_jpim1
	452	!- Calculate Delta on corners
	453	zddc = ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1) &
	454	& -v_ice1(ji,jj)/e1u(ji,jj) &
	455	& )e1f(ji,jj)e1f(ji,jj) &
	456	& +( u_ice2(ji+1,jj)/e2v(ji+1,jj) &
	457	& -u_ice2(ji,jj)/e2v(ji,jj) &
	458	& )e2f(ji,jj)e2f(ji,jj) &
	459	& ) &
	460	& / ( e1f(ji,jj) * e2f(ji,jj) )
[825]	461
[921]	462	zdtc = (-( v_ice1(ji,jj+1)/e1u(ji,jj+1) &
	463	& -v_ice1(ji,jj)/e1u(ji,jj) &
	464	& )e1f(ji,jj)e1f(ji,jj) &
	465	& +( u_ice2(ji+1,jj)/e2v(ji+1,jj) &
	466	& -u_ice2(ji,jj)/e2v(ji,jj) &
	467	& )e2f(ji,jj)e2f(ji,jj) &
	468	& ) &
	469	& / ( e1f(ji,jj) * e2f(ji,jj) )
[825]	470
[921]	471	deltac(ji,jj) = SQRT(zddc2+(zdtc2+zds(ji,jj)*2)usecc2) + creepl
[825]	472
[921]	473	!-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide).
[4269]	474	!zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) / &
	475	! & ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) ) &
	476	! & / ( 1._wp + alphaevp * ecc2 * dtotel )
[825]	477
[4205]	478	! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp)
[4269]	479	zs12(ji,jj) = ( zs12(ji,jj) + dtotel * &
	480	& ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) ) &
	481	& / ( 1.0 + dtotel )
[4205]	482
[921]	483	END DO ! ji
	484	END DO ! jj
[825]	485
[921]	486	CALL lbc_lnk( zs12(:,:), 'F', 1. )
[825]	487
[4161]	488	!#if defined key_bdy
	489	! ! clem: change zs1, zs2, zs12 at the boundary for each iteration
	490	! CALL bdy_ice_lim_dyn( 2, zs1, zs2, zs12 )
	491	! CALL lbc_lnk( zs1 (:,:), 'T', 1. )
	492	! CALL lbc_lnk( zs2 (:,:), 'T', 1. )
	493	! CALL lbc_lnk( zs12(:,:), 'F', 1. )
	494	!#endif
	495
[921]	496	! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002)
	497	DO jj = k_j1+1, k_jpj-1
	498	DO ji = fs_2, fs_jpim1
	499	!- contribution of zs1, zs2 and zs12 to zf1
	500	zf1(ji,jj) = 0.5( (zs1(ji+1,jj)-zs1(ji,jj))e2u(ji,jj) &
	501	& +(zs2(ji+1,jj)e2t(ji+1,jj)2-zs2(ji,jj)e2t(ji,jj)**2)/e2u(ji,jj) &
	502	& +2.0(zs12(ji,jj)e1f(ji,jj)*2-zs12(ji,jj-1)e1f(ji,jj-1)**2)/e1u(ji,jj) &
	503	& ) / ( e1u(ji,jj)*e2u(ji,jj) )
	504	! contribution of zs1, zs2 and zs12 to zf2
	505	zf2(ji,jj) = 0.5( (zs1(ji,jj+1)-zs1(ji,jj))e1v(ji,jj) &
	506	& -(zs2(ji,jj+1)e1t(ji,jj+1)2 - zs2(ji,jj)e1t(ji,jj)**2)/e1v(ji,jj) &
	507	& + 2.0(zs12(ji,jj)e2f(ji,jj)**2 - &
	508	zs12(ji-1,jj)e2f(ji-1,jj)*2)/e2v(ji,jj) &
	509	& ) / ( e1v(ji,jj)*e2v(ji,jj) )
	510	END DO
	511	END DO
[825]	512	!
	513	! Computation of ice velocity
	514	!
	515	! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly.
	516	!
[921]	517	IF (MOD(jter,2).eq.0) THEN
[825]	518
[868]	519	!CDIR NOVERRCHK
[921]	520	DO jj = k_j1+1, k_jpj-1
[868]	521	!CDIR NOVERRCHK
[921]	522	DO ji = fs_2, fs_jpim1
	523	zmask = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj)
	524	zsang = SIGN ( 1.0 , fcor(ji,jj) ) * sangvg
	525	z0 = zmass1(ji,jj)/dtevp
[825]	526
[921]	527	! SB modif because ocean has no slip boundary condition
	528	zv_ice1 = 0.5( (v_ice(ji,jj)+v_ice(ji,jj-1))e1t(ji,jj) &
	529	& +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) &
	530	& /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
	531	za = rhocoSQRT((u_ice(ji,jj)-u_oce1(ji,jj))*2 + &
	532	(zv_ice1-v_oce1(ji,jj))*2) (1.0-zfrld1(ji,jj))
	533	zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
	534	za(cangvgu_oce1(ji,jj)-zsang*v_oce1(ji,jj))
	535	zcca = z0+za*cangvg
	536	zccb = zcorl1(ji,jj)+za*zsang
	537	u_ice(ji,jj) = (zr+zccbzv_ice1)/(zcca+epsd)zmask
[825]	538
[921]	539	END DO
	540	END DO
[825]	541
[921]	542	CALL lbc_lnk( u_ice(:,:), 'U', -1. )
[4161]	543	#if defined key_agrif && defined key_lim2
[3680]	544	CALL agrif_rhg_lim2( jter, nevp, 'U' )
	545	#endif
[825]	546
[868]	547	!CDIR NOVERRCHK
[921]	548	DO jj = k_j1+1, k_jpj-1
[868]	549	!CDIR NOVERRCHK
[921]	550	DO ji = fs_2, fs_jpim1
[834]	551
[921]	552	zmask = (1.0-MAX(rzero,SIGN(rone,-zmass2(ji,jj))))*tmv(ji,jj)
	553	zsang = SIGN(1.0,fcor(ji,jj))*sangvg
	554	z0 = zmass2(ji,jj)/dtevp
	555	! SB modif because ocean has no slip boundary condition
	556	zu_ice2 = 0.5( (u_ice(ji,jj)+u_ice(ji-1,jj))e2t(ji,jj) &
	557	& + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) &
	558	& /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
	559	za = rhocoSQRT((zu_ice2-u_oce2(ji,jj))*2 + &
	560	(v_ice(ji,jj)-v_oce2(ji,jj))*2)(1.0-zfrld2(ji,jj))
	561	zr = z0*v_ice(ji,jj) + zf2(ji,jj) + &
	562	za2ct(ji,jj) + za(cangvgv_oce2(ji,jj)+zsang*u_oce2(ji,jj))
	563	zcca = z0+za*cangvg
	564	zccb = zcorl2(ji,jj)+za*zsang
	565	v_ice(ji,jj) = (zr-zccbzu_ice2)/(zcca+epsd)zmask
[825]	566
[921]	567	END DO
	568	END DO
[825]	569
[921]	570	CALL lbc_lnk( v_ice(:,:), 'V', -1. )
[4161]	571	#if defined key_agrif && defined key_lim2
[3680]	572	CALL agrif_rhg_lim2( jter, nevp, 'V' )
	573	#endif
[825]	574
[834]	575	ELSE
[868]	576	!CDIR NOVERRCHK
[921]	577	DO jj = k_j1+1, k_jpj-1
[868]	578	!CDIR NOVERRCHK
[921]	579	DO ji = fs_2, fs_jpim1
	580	zmask = (1.0-MAX(rzero,SIGN(rone,-zmass2(ji,jj))))*tmv(ji,jj)
	581	zsang = SIGN(1.0,fcor(ji,jj))*sangvg
	582	z0 = zmass2(ji,jj)/dtevp
	583	! SB modif because ocean has no slip boundary condition
	584	zu_ice2 = 0.5( (u_ice(ji,jj)+u_ice(ji-1,jj))e2t(ji,jj) &
	585	& +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) &
	586	& /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
[825]	587
[921]	588	za = rhocoSQRT((zu_ice2-u_oce2(ji,jj))*2 + &
	589	(v_ice(ji,jj)-v_oce2(ji,jj))*2)(1.0-zfrld2(ji,jj))
	590	zr = z0*v_ice(ji,jj) + zf2(ji,jj) + &
	591	za2ct(ji,jj) + za(cangvgv_oce2(ji,jj)+zsang*u_oce2(ji,jj))
	592	zcca = z0+za*cangvg
	593	zccb = zcorl2(ji,jj)+za*zsang
	594	v_ice(ji,jj) = (zr-zccbzu_ice2)/(zcca+epsd)zmask
[825]	595
[921]	596	END DO
	597	END DO
[825]	598
[921]	599	CALL lbc_lnk( v_ice(:,:), 'V', -1. )
[4161]	600	#if defined key_agrif && defined key_lim2
[3680]	601	CALL agrif_rhg_lim2( jter, nevp , 'V' )
	602	#endif
[825]	603
[868]	604	!CDIR NOVERRCHK
[921]	605	DO jj = k_j1+1, k_jpj-1
[868]	606	!CDIR NOVERRCHK
[921]	607	DO ji = fs_2, fs_jpim1
	608	zmask = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj)
	609	zsang = SIGN(1.0,fcor(ji,jj))*sangvg
	610	z0 = zmass1(ji,jj)/dtevp
	611	! SB modif because ocean has no slip boundary condition
	612	! GG Bug
	613	! zv_ice1 = 0.5( (v_ice(ji,jj)+v_ice(ji,jj-1))e1t(ji+1,jj) &
	614	! & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
	615	! & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
	616	zv_ice1 = 0.5( (v_ice(ji,jj)+v_ice(ji,jj-1))e1t(ji,jj) &
	617	& +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) &
	618	& /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
[825]	619
[921]	620	za = rhocoSQRT((u_ice(ji,jj)-u_oce1(ji,jj))*2 + &
	621	(zv_ice1-v_oce1(ji,jj))*2)(1.0-zfrld1(ji,jj))
	622	zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
	623	za(cangvgu_oce1(ji,jj)-zsang*v_oce1(ji,jj))
	624	zcca = z0+za*cangvg
	625	zccb = zcorl1(ji,jj)+za*zsang
	626	u_ice(ji,jj) = (zr+zccbzv_ice1)/(zcca+epsd)zmask
	627	END DO ! ji
	628	END DO ! jj
[825]	629
[921]	630	CALL lbc_lnk( u_ice(:,:), 'U', -1. )
[4161]	631	#if defined key_agrif && defined key_lim2
[3680]	632	CALL agrif_rhg_lim2( jter, nevp, 'U' )
	633	#endif
[825]	634
[921]	635	ENDIF
[4161]	636
	637	!#if defined key_bdy
	638	! ! clem: change u_ice and v_ice at the boundary for each iteration
	639	! CALL bdy_ice_lim_dyn( 1 )
	640	!#endif
[825]	641
[921]	642	IF(ln_ctl) THEN
	643	!--- Convergence test.
	644	DO jj = k_j1+1 , k_jpj-1
	645	zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , &
	646	ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
	647	END DO
	648	zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) )
	649	IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain
	650	ENDIF
	651
[4161]	652	! ! ==================== !
[868]	653	END DO ! end loop over jter !
[825]	654	! ! ==================== !
[921]	655	!
	656	!------------------------------------------------------------------------------!
	657	! 4) Prevent ice velocities when the ice is thin
	658	!------------------------------------------------------------------------------!
[4161]	659	!clem : add hminrhg in the namelist
[921]	660	!
[4161]	661	! If the ice thickness is below hminrhg (5cm) then ice velocity should equal the
[834]	662	! ocean velocity,
	663	! This prevents high velocity when ice is thin
[868]	664	!CDIR NOVERRCHK
[825]	665	DO jj = k_j1+1, k_jpj-1
[868]	666	!CDIR NOVERRCHK
	667	DO ji = fs_2, fs_jpim1
[825]	668	zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) )
[4161]	669	!zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
	670	zdummy = vt_i(ji,jj)
	671	IF ( zdummy .LE. hminrhg ) THEN
[888]	672	u_ice(ji,jj) = u_oce(ji,jj)
	673	v_ice(ji,jj) = v_oce(ji,jj)
[825]	674	ENDIF ! zdummy
	675	END DO
	676	END DO
[866]	677
[869]	678	CALL lbc_lnk( u_ice(:,:), 'U', -1. )
	679	CALL lbc_lnk( v_ice(:,:), 'V', -1. )
[4161]	680	#if defined key_agrif && defined key_lim2
[3680]	681	CALL agrif_rhg_lim2( nevp , nevp, 'U' )
	682	CALL agrif_rhg_lim2( nevp , nevp, 'V' )
	683	#endif
[4161]	684	#if defined key_bdy
	685	! clem: change u_ice and v_ice at the boundary
	686	CALL bdy_ice_lim_dyn( 1 )
	687	#endif
[869]	688
[868]	689	DO jj = k_j1+1, k_jpj-1
	690	DO ji = fs_2, fs_jpim1
	691	zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) )
[4161]	692	!zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
	693	zdummy = vt_i(ji,jj)
	694	IF ( zdummy .LE. hminrhg ) THEN
[921]	695	v_ice1(ji,jj) = 0.5( (v_ice(ji,jj)+v_ice(ji,jj-1))e1t(ji+1,jj) &
	696	& +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
	697	& /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
[868]	698
[921]	699	u_ice2(ji,jj) = 0.5( (u_ice(ji,jj)+u_ice(ji-1,jj))e2t(ji,jj+1) &
	700	& +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &
	701	& /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
	702	ENDIF ! zdummy
[868]	703	END DO
	704	END DO
	705
[869]	706	CALL lbc_lnk( u_ice2(:,:), 'V', -1. )
	707	CALL lbc_lnk( v_ice1(:,:), 'U', -1. )
[868]	708
[866]	709	! Recompute delta, shear and div, inputs for mechanical redistribution
[868]	710	!CDIR NOVERRCHK
[825]	711	DO jj = k_j1+1, k_jpj-1
[868]	712	!CDIR NOVERRCHK
[988]	713	DO ji = fs_2, jpim1 !RB bug no vect opt due to tmi
[825]	714	!- zdd(:,:), zdt(:,:): divergence and tension at centre
	715	!- zds(:,:): shear on northeast corner of grid cells
	716	zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) )
[4161]	717	!zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
	718	zdummy = vt_i(ji,jj)
	719	IF ( zdummy .LE. hminrhg ) THEN
[825]	720
[921]	721	zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) &
	722	& -e2u(ji-1,jj)*u_ice(ji-1,jj) &
	723	& +e1v(ji,jj)*v_ice(ji,jj) &
	724	& -e1v(ji,jj-1)*v_ice(ji,jj-1) &
	725	& ) &
	726	& / area(ji,jj)
[825]	727
[921]	728	zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) &
	729	& -u_ice(ji-1,jj)/e2u(ji-1,jj) &
	730	& )e2t(ji,jj)e2t(ji,jj) &
	731	& -( v_ice(ji,jj)/e1v(ji,jj) &
	732	& -v_ice(ji,jj-1)/e1v(ji,jj-1) &
	733	& )e1t(ji,jj)e1t(ji,jj) &
	734	& ) &
	735	& / area(ji,jj)
	736	!
	737	! SB modif because ocean has no slip boundary condition
	738	zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) &
	739	& - u_ice(ji,jj) / e1u(ji,jj) ) &
	740	& * e1f(ji,jj) * e1f(ji,jj) &
	741	& + ( v_ice(ji+1,jj) / e2v(ji+1,jj) &
	742	& - v_ice(ji,jj) / e2v(ji,jj) ) &
	743	& * e2f(ji,jj) * e2f(ji,jj) ) &
	744	& / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) &
	745	& * tmi(ji,jj) * tmi(ji,jj+1) &
	746	& * tmi(ji+1,jj) * tmi(ji+1,jj+1)
[825]	747
[3791]	748	zdst(ji,jj) = ( e2u( ji , jj ) * v_ice1(ji ,jj ) &
	749	& - e2u( ji-1, jj ) * v_ice1(ji-1,jj ) &
	750	& + e1v( ji , jj ) * u_ice2(ji ,jj ) &
	751	& - e1v( ji , jj-1 ) * u_ice2(ji ,jj-1) ) / area(ji,jj)
[825]	752
[4161]	753	! deltat(ji,jj) = SQRT( zdd(ji,jj)*zdd(ji,jj) &
	754	! & + ( zdt(ji,jj)zdt(ji,jj) + zdst(ji,jj)zdst(ji,jj) ) * usecc2 &
	755	! & ) + creepl
	756	! MV rewriting
	757	delta = SQRT( zdd(ji,jj)zdd(ji,jj) + ( zdt(ji,jj)zdt(ji,jj) + zdst(ji,jj)zdst(ji,jj) ) usecc2 )
	758	deltat(ji,jj) = delta + creepl
	759	! END MV
	760
[921]	761	ENDIF ! zdummy
[825]	762
	763	END DO !jj
	764	END DO !ji
[921]	765	!
	766	!------------------------------------------------------------------------------!
	767	! 5) Store stress tensor and its invariants
	768	!------------------------------------------------------------------------------!
	769	!
[866]	770	! * Invariants of the stress tensor are required for limitd_me
[3791]	771	! (accelerates convergence and improves stability)
[866]	772	DO jj = k_j1+1, k_jpj-1
[868]	773	DO ji = fs_2, fs_jpim1
	774	divu_i (ji,jj) = zdd (ji,jj)
	775	delta_i(ji,jj) = deltat(ji,jj)
[4161]	776	! begin TECLIM change
	777	zdst(ji,jj)= ( e2u( ji , jj ) * v_ice1(ji,jj) &
	778	& - e2u( ji-1, jj ) * v_ice1(ji-1,jj) &
	779	& + e1v( ji , jj ) * u_ice2(ji,jj) &
	780	& - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj)
[3791]	781	shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst(ji,jj) * zdst(ji,jj) )
[4161]	782	! end TECLIM change
[825]	783	END DO
[866]	784	END DO
[4161]	785
	786	! Lateral boundary condition
	787	CALL lbc_lnk( divu_i (:,:), 'T', 1. )
[825]	788	CALL lbc_lnk( delta_i(:,:), 'T', 1. )
[4161]	789	! CALL lbc_lnk( shear_i(:,:), 'F', 1. )
[3791]	790	CALL lbc_lnk( shear_i(:,:), 'T', 1. )
[866]	791
[868]	792	! * Store the stress tensor for the next time step
	793	stress1_i (:,:) = zs1 (:,:)
	794	stress2_i (:,:) = zs2 (:,:)
	795	stress12_i(:,:) = zs12(:,:)
	796
[921]	797	!
	798	!------------------------------------------------------------------------------!
	799	! 6) Control prints of residual and charge ellipse
	800	!------------------------------------------------------------------------------!
	801	!
[834]	802	! print the residual for convergence
	803	IF(ln_ctl) THEN
[868]	804	WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, jter
[834]	805	CALL prt_ctl_info(charout)
	806	CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
	807	ENDIF
[825]	808
[834]	809	! print charge ellipse
	810	! This can be desactivated once the user is sure that the stress state
	811	! lie on the charge ellipse. See Bouillon et al. 08 for more details
[825]	812	IF(ln_ctl) THEN
	813	CALL prt_ctl_info('lim_rhg : numit :',ivar1=numit)
	814	CALL prt_ctl_info('lim_rhg : nwrite :',ivar1=nwrite)
	815	CALL prt_ctl_info('lim_rhg : MOD :',ivar1=MOD(numit,nwrite))
	816	IF( MOD(numit,nwrite) .EQ. 0 ) THEN
	817	WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. '
	818	CALL prt_ctl_info(charout)
	819	DO jj = k_j1+1, k_jpj-1
	820	DO ji = 2, jpim1
	821	IF (zpresh(ji,jj) .GT. 1.0) THEN
	822	sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)*2 + 4zs12(ji,jj)2 )0.5 ) / ( 2*zpresh(ji,jj) )
	823	sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)*2 + 4zs12(ji,jj)2 )0.5 ) / ( 2*zpresh(ji,jj) )
	824	WRITE(charout,FMT="('lim_rhg :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
	825	CALL prt_ctl_info(charout)
	826	ENDIF
	827	END DO
	828	END DO
	829	WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 2000, numit, 0, 0, ' ch. ell. '
	830	CALL prt_ctl_info(charout)
	831	ENDIF
	832	ENDIF
[2715]	833	!
[3294]	834	CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
	835	CALL wrk_dealloc( jpi,jpj, zc1 , u_oce1, u_oce2, u_ice2, zusw , v_oce1 , v_oce2, v_ice1 )
[3791]	836	CALL wrk_dealloc( jpi,jpj, zf1 , deltat, zu_ice, zf2 , deltac, zv_ice , zdd , zdt , zds , zdst )
[3625]	837	CALL wrk_dealloc( jpi,jpj, zdd , zdt , zds , zs1 , zs2 , zs12 , zresr , zpice )
[3294]	838
[825]	839	END SUBROUTINE lim_rhg
	840
	841	#else
	842	!!----------------------------------------------------------------------
	843	!! Default option Dummy module NO LIM sea-ice model
	844	!!----------------------------------------------------------------------
	845	CONTAINS
	846	SUBROUTINE lim_rhg( k1 , k2 ) ! Dummy routine
	847	WRITE(,) 'lim_rhg: You should not have seen this print! error?', k1, k2
	848	END SUBROUTINE lim_rhg
	849	#endif
	850
	851	!!==============================================================================
	852	END MODULE limrhg

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

Download in other formats: