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

Context Navigation

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 4634

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

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

Download in other formats: