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

Context Navigation

source: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 6584

Last change on this file since 6584 was 6584, checked in by clem, 8 years ago
LIM3 and Agrif compatibility
Property svn:keywords set to `Id`
File size: 37.1 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
[6584]	43	#if defined key_agrif && defined key_lim3
	44	USE agrif_lim3_interp
	45	#endif
[4161]	46	#if defined key_bdy
	47	USE bdyice_lim
	48	#endif
[825]	49
	50	IMPLICIT NONE
	51	PRIVATE
	52
[2715]	53	PUBLIC lim_rhg ! routine called by lim_dyn (or lim_dyn_2)
[825]	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
[5123]	107	!! nn_nevp, elastic time scale and rn_creepl maintain stress state
[834]	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	!!-------------------------------------------------------------------
	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 !
[5123]	121	REAL(wp) :: za, zstms ! local scalars
	122	REAL(wp) :: zc1, zc2, zc3 ! ice mass
[825]	123
[5123]	124	REAL(wp) :: dtevp , z1_dtevp ! time step for subcycling
	125	REAL(wp) :: dtotel, z1_dtotel, ecc2, ecci ! square of yield ellipse eccenticity
	126	REAL(wp) :: z0, zr, zcca, zccb ! temporary scalars
	127	REAL(wp) :: zu_ice2, zv_ice1 !
	128	REAL(wp) :: zddc, zdtc ! delta on corners and on centre
	129	REAL(wp) :: zdst ! shear at the center of the grid point
	130	REAL(wp) :: zdsshx, zdsshy ! term for the gradient of ocean surface
	131	REAL(wp) :: sigma1, sigma2 ! internal ice stress
[825]	132
[2715]	133	REAL(wp) :: zresm ! Maximal error on ice velocity
[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
[5123]	142	REAL(wp), POINTER, DIMENSION(:,:) :: v_oce1 ! ocean u/v component on U points
	143	REAL(wp), POINTER, DIMENSION(:,:) :: u_oce2 ! ocean u/v component on V points
[3294]	144	REAL(wp), POINTER, DIMENSION(:,:) :: u_ice2, v_ice1 ! ice u/v component on V/U point
	145	REAL(wp), POINTER, DIMENSION(:,:) :: zf1 , zf2 ! arrays for internal stresses
[5123]	146	REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! mask ocean grid points
[3294]	147
[4990]	148	REAL(wp), POINTER, DIMENSION(:,:) :: zdt ! tension at centre of grid cells
[3294]	149	REAL(wp), POINTER, DIMENSION(:,:) :: zds ! Shear on northeast corner of grid cells
	150	REAL(wp), POINTER, DIMENSION(:,:) :: zs1 , zs2 ! Diagonal stress tensor components zs1 and zs2
	151	REAL(wp), POINTER, DIMENSION(:,:) :: zs12 ! Non-diagonal stress tensor component zs12
	152	REAL(wp), POINTER, DIMENSION(:,:) :: zu_ice, zv_ice, zresr ! Local error on velocity
[3625]	153	REAL(wp), POINTER, DIMENSION(:,:) :: zpice ! array used for the calculation of ice surface slope:
	154	! ocean surface (ssh_m) if ice is not embedded
[3791]	155	! ice top surface if ice is embedded
[5123]	156
	157	REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter
	158	REAL(wp), PARAMETER :: zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity
[2528]	159	!!-------------------------------------------------------------------
[3294]	160
	161	CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
[5123]	162	CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask )
[4990]	163	CALL wrk_alloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds )
[5888]	164	CALL wrk_alloc( jpi,jpj, zs1 , zs2 , zs12 , zresr , zpice )
[3294]	165
[2528]	166	#if defined key_lim2 && ! defined key_lim2_vp
	167	# if defined key_agrif
[5123]	168	USE ice_2, vt_s => hsnm
	169	USE ice_2, vt_i => hicm
[2528]	170	# else
[5123]	171	vt_s => hsnm
	172	vt_i => hicm
[2528]	173	# endif
[5123]	174	at_i(:,:) = 1. - frld(:,:)
[2528]	175	#endif
[3680]	176	#if defined key_agrif && defined key_lim2
[5123]	177	CALL agrif_rhg_lim2_load ! First interpolation of coarse values
[3680]	178	#endif
[6584]	179	#if defined key_agrif && defined key_lim3
	180	CALL agrif_interp_lim3('U') ! First interpolation of coarse values
	181	CALL agrif_interp_lim3('V') ! First interpolation of coarse values
	182	#endif
[921]	183	!
	184	!------------------------------------------------------------------------------!
[4990]	185	! 1) Ice strength (zpresh) !
[921]	186	!------------------------------------------------------------------------------!
	187	!
[825]	188	! Put every vector to 0
[4990]	189	delta_i(:,:) = 0._wp ;
	190	zpresh (:,:) = 0._wp ;
[2528]	191	zpreshc(:,:) = 0._wp
	192	u_ice2 (:,:) = 0._wp ; v_ice1(:,:) = 0._wp
[4990]	193	divu_i (:,:) = 0._wp ; zdt (:,:) = 0._wp ; zds(:,:) = 0._wp
	194	shear_i(:,:) = 0._wp
[825]	195
[2528]	196	#if defined key_lim3
[5123]	197	CALL lim_itd_me_icestrength( nn_icestr ) ! LIM-3: Ice strength on T-points
[2528]	198	#endif
[825]	199
[2528]	200	DO jj = k_j1 , k_jpj ! Ice mass and temp variables
[825]	201	DO ji = 1 , jpi
[2528]	202	#if defined key_lim3
[5123]	203	zpresh(ji,jj) = tmask(ji,jj,1) * strength(ji,jj)
[2528]	204	#endif
[2580]	205	#if defined key_lim2
[5123]	206	zpresh(ji,jj) = tmask(ji,jj,1) * pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) )
[2580]	207	#endif
[5123]	208	! zmask = 1 where there is ice or on land
	209	zmask(ji,jj) = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tmask(ji,jj,1)
[825]	210	END DO
	211	END DO
	212
[834]	213	! Ice strength on grid cell corners (zpreshc)
	214	! needed for calculation of shear stress
[825]	215	DO jj = k_j1+1, k_jpj-1
[868]	216	DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1)
[5123]	217	zstms = tmask(ji+1,jj+1,1) * wght(ji+1,jj+1,2,2) + tmask(ji,jj+1,1) * wght(ji+1,jj+1,1,2) + &
	218	& tmask(ji+1,jj,1) * wght(ji+1,jj+1,2,1) + tmask(ji,jj,1) * wght(ji+1,jj+1,1,1)
	219	zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + &
	220	& zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + zpresh(ji,jj) * wght(ji+1,jj+1,1,1) &
	221	& ) / MAX( zstms, zepsi )
[825]	222	END DO
	223	END DO
	224	CALL lbc_lnk( zpreshc(:,:), 'F', 1. )
[921]	225	!
	226	!------------------------------------------------------------------------------!
	227	! 2) Wind / ocean stress, mass terms, coriolis terms
	228	!------------------------------------------------------------------------------!
	229	!
[825]	230	! Wind stress, coriolis and mass terms on the sides of the squares
	231	! zfrld1: lead fraction on U-points
	232	! zfrld2: lead fraction on V-points
	233	! zmass1: ice/snow mass on U-points
	234	! zmass2: ice/snow mass on V-points
	235	! zcorl1: Coriolis parameter on U-points
	236	! zcorl2: Coriolis parameter on V-points
	237	! (ztagnx,ztagny): wind stress on U/V points
	238	! v_oce1: ocean v component on u points
	239	! u_oce2: ocean u component on v points
[921]	240
[3625]	241	IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==!
[5123]	242	!
	243	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
	244	! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
[3625]	245	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
[5123]	246	!
	247	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
	248	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
[3625]	249	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
[5123]	250	!
[3625]	251	zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
[5123]	252	!
[3625]	253	ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==!
	254	zpice(:,:) = ssh_m(:,:)
	255	ENDIF
	256
[825]	257	DO jj = k_j1+1, k_jpj-1
[868]	258	DO ji = fs_2, fs_jpim1
[825]	259
[5123]	260	zc1 = tmask(ji ,jj ,1) * ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) )
	261	zc2 = tmask(ji+1,jj ,1) * ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) )
	262	zc3 = tmask(ji ,jj+1,1) * ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) )
[4990]	263
[5123]	264	zt11 = tmask(ji ,jj,1) * e1t(ji ,jj)
	265	zt12 = tmask(ji+1,jj,1) * e1t(ji+1,jj)
	266	zt21 = tmask(ji,jj ,1) * e2t(ji,jj )
	267	zt22 = tmask(ji,jj+1,1) * e2t(ji,jj+1)
[825]	268
	269	! Leads area.
[5123]	270	zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + zepsi )
	271	zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + zepsi )
[825]	272
	273	! Mass, coriolis coeff. and currents
[5123]	274	zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi )
	275	zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi )
	276	zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) ) &
	277	& / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi )
	278	zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) ) &
	279	& / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi )
[825]	280	!
[834]	281	! Ocean has no slip boundary condition
[5123]	282	v_oce1(ji,jj) = 0.5 * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji,jj) &
	283	& + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) ) &
	284	& / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
[825]	285
[5123]	286	u_oce2(ji,jj) = 0.5 * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj) &
	287	& + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) ) &
	288	& / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
[825]	289
[1469]	290	! Wind stress at U,V-point
	291	ztagnx = ( 1. - zfrld1(ji,jj) ) * utau_ice(ji,jj)
	292	ztagny = ( 1. - zfrld2(ji,jj) ) * vtau_ice(ji,jj)
[825]	293
[834]	294	! Computation of the velocity field taking into account the ice internal interaction.
[825]	295	! Terms that are independent of the velocity field.
	296
	297	! SB On utilise maintenant le gradient de la pente de l'ocean
	298	! include it later
[834]	299
[5123]	300	zdsshx = ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
	301	zdsshy = ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
[825]	302
	303	za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx
	304	za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy
	305
	306	END DO
	307	END DO
	308
[921]	309	!
	310	!------------------------------------------------------------------------------!
	311	! 3) Solution of the momentum equation, iterative procedure
	312	!------------------------------------------------------------------------------!
	313	!
[825]	314	! Time step for subcycling
[5123]	315	dtevp = rdt_ice / nn_nevp
	316	#if defined key_lim3
	317	dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice )
	318	#else
[2528]	319	dtotel = dtevp / ( 2._wp * telast )
[5123]	320	#endif
	321	z1_dtotel = 1._wp / ( 1._wp + dtotel )
	322	z1_dtevp = 1._wp / dtevp
[825]	323	!-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)
[5123]	324	ecc2 = rn_ecc * rn_ecc
[4205]	325	ecci = 1. / ecc2
[825]	326
	327	!-Initialise stress tensor
[2528]	328	zs1 (:,:) = stress1_i (:,:)
	329	zs2 (:,:) = stress2_i (:,:)
[866]	330	zs12(:,:) = stress12_i(:,:)
[825]	331
[2528]	332	! !----------------------!
[5123]	333	DO jter = 1 , nn_nevp ! loop over jter !
[2528]	334	! !----------------------!
[825]	335	DO jj = k_j1, k_jpj-1
[2528]	336	zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
[825]	337	zv_ice(:,jj) = v_ice(:,jj)
[921]	338	END DO
[825]	339
[834]	340	DO jj = k_j1+1, k_jpj-1
[5123]	341	DO ji = fs_2, fs_jpim1 !RB bug no vect opt due to zmask
[825]	342
[921]	343	!
	344	!- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002)
[4990]	345	!- divu_i(:,:), zdt(:,:): divergence and tension at centre of grid cells
[921]	346	!- zds(:,:): shear on northeast corner of grid cells
	347	!
	348	!- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded,
	349	! there are many repeated calculations.
	350	! Speed could be improved by regrouping terms. For
	351	! the moment, however, the stress is on clarity of coding to avoid
	352	! bugs (Martin, for Miguel).
	353	!
[4990]	354	!- ALSO: arrays zdt, zds and delta could
[921]	355	! be removed in the future to minimise memory demand.
	356	!
	357	!- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of
	358	! grid cells, exactly as in the B grid case. For simplicity, the indexation on
	359	! the corners is the same as in the B grid.
	360	!
	361	!
[5123]	362	divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
	363	& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
	364	& ) * r1_e12t(ji,jj)
[825]	365
[5123]	366	zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) &
	367	& - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
	368	& ) * r1_e12t(ji,jj)
[825]	369
[921]	370	!
[5123]	371	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) &
	372	& + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
	373	& ) * r1_e12f(ji,jj) * ( 2._wp - fmask(ji,jj,1) ) &
	374	& * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1)
[825]	375
	376
[5123]	377	v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) &
	378	& + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) &
	379	& / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
[825]	380
[5123]	381	u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) &
	382	& + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) &
	383	& / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
[921]	384	END DO
	385	END DO
	386
[5429]	387	CALL lbc_lnk_multi( v_ice1, 'U', -1., u_ice2, 'V', -1. ) ! lateral boundary cond.
	388
[921]	389	DO jj = k_j1+1, k_jpj-1
	390	DO ji = fs_2, fs_jpim1
[825]	391
[921]	392	!- Calculate Delta at centre of grid cells
[5123]	393	zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) &
	394	& + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji ,jj-1) * u_ice2(ji ,jj-1) &
	395	& ) * r1_e12t(ji,jj)
[825]	396
[5123]	397	delta = SQRT( divu_i(ji,jj)2 + ( zdt(ji,jj)2 + zdst*2 ) usecc2 )
	398	delta_i(ji,jj) = delta + rn_creepl
[4205]	399
[921]	400	!- Calculate Delta on corners
[5123]	401	zddc = ( ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &
	402	& + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
	403	& ) * r1_e12f(ji,jj)
[825]	404
[5123]	405	zdtc = (- ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &
	406	& + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
	407	& ) * r1_e12f(ji,jj)
[825]	408
[5123]	409	zddc = SQRT( zddc2 + ( zdtc2 + zds(ji,jj)*2 ) usecc2 ) + rn_creepl
[825]	410
[5123]	411	!-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide).
	412	zs1(ji,jj) = ( zs1 (ji,jj) + dtotel * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj) * zpresh(ji,jj) &
	413	& ) * z1_dtotel
	414	zs2(ji,jj) = ( zs2 (ji,jj) + dtotel * ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) &
	415	& ) * z1_dtotel
	416	!-Calculate stress tensor component zs12 at corners
	417	zs12(ji,jj) = ( zs12(ji,jj) + dtotel * ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) &
	418	& ) * z1_dtotel
[4205]	419
[5123]	420	END DO
	421	END DO
[825]	422
[5429]	423	CALL lbc_lnk_multi( zs1 , 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
	424
[921]	425	! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002)
	426	DO jj = k_j1+1, k_jpj-1
	427	DO ji = fs_2, fs_jpim1
	428	!- contribution of zs1, zs2 and zs12 to zf1
[5123]	429	zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) &
	430	& + ( zs2(ji+1,jj) * e2t(ji+1,jj)*2 - zs2(ji,jj) e2t(ji,jj)*2 ) r1_e2u(ji,jj) &
	431	& + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)*2 - zs12(ji,jj-1) e1f(ji,jj-1)*2 ) r1_e1u(ji,jj) &
	432	& ) * r1_e12u(ji,jj)
[921]	433	! contribution of zs1, zs2 and zs12 to zf2
[5123]	434	zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) &
	435	& - ( zs2(ji,jj+1) * e1t(ji,jj+1)*2 - zs2(ji,jj) e1t(ji,jj)*2 ) r1_e1v(ji,jj) &
	436	& + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)*2 - zs12(ji-1,jj) e2f(ji-1,jj)*2 ) r1_e2v(ji,jj) &
	437	& ) * r1_e12v(ji,jj)
[921]	438	END DO
	439	END DO
[825]	440	!
	441	! Computation of ice velocity
	442	!
	443	! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly.
	444	!
[921]	445	IF (MOD(jter,2).eq.0) THEN
[825]	446
[921]	447	DO jj = k_j1+1, k_jpj-1
	448	DO ji = fs_2, fs_jpim1
[5123]	449	rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1)
	450	z0 = zmass1(ji,jj) * z1_dtevp
[825]	451
[921]	452	! SB modif because ocean has no slip boundary condition
[5123]	453	zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji ,jj) &
	454	& + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) &
	455	& / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
	456	za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + &
	457	& ( zv_ice1 - v_oce1(ji,jj) )*2 ) ( 1.0 - zfrld1(ji,jj) )
	458	zr = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj)
	459	zcca = z0 + za
[4990]	460	zccb = zcorl1(ji,jj)
[5123]	461	u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch
[921]	462	END DO
	463	END DO
[825]	464
[921]	465	CALL lbc_lnk( u_ice(:,:), 'U', -1. )
[4161]	466	#if defined key_agrif && defined key_lim2
[5123]	467	CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
[3680]	468	#endif
[6584]	469	#if defined key_agrif && defined key_lim3
	470	CALL agrif_interp_lim3( 'U' )
	471	#endif
[4333]	472	#if defined key_bdy
	473	CALL bdy_ice_lim_dyn( 'U' )
	474	#endif
[825]	475
[921]	476	DO jj = k_j1+1, k_jpj-1
	477	DO ji = fs_2, fs_jpim1
[834]	478
[5123]	479	rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1)
	480	z0 = zmass2(ji,jj) * z1_dtevp
[921]	481	! SB modif because ocean has no slip boundary condition
[5123]	482	zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) &
	483	& + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) &
	484	& / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
	485	za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + &
	486	& ( v_ice(ji,jj) - v_oce(ji,jj))*2 ) ( 1.0 - zfrld2(ji,jj) )
	487	zr = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj)
	488	zcca = z0 + za
[4990]	489	zccb = zcorl2(ji,jj)
[5123]	490	v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch
[921]	491	END DO
	492	END DO
[825]	493
[921]	494	CALL lbc_lnk( v_ice(:,:), 'V', -1. )
[4161]	495	#if defined key_agrif && defined key_lim2
[5123]	496	CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
[3680]	497	#endif
[6584]	498	#if defined key_agrif && defined key_lim3
	499	CALL agrif_interp_lim3( 'V' )
	500	#endif
[4333]	501	#if defined key_bdy
	502	CALL bdy_ice_lim_dyn( 'V' )
	503	#endif
[825]	504
[834]	505	ELSE
[921]	506	DO jj = k_j1+1, k_jpj-1
	507	DO ji = fs_2, fs_jpim1
[5123]	508	rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1)
	509	z0 = zmass2(ji,jj) * z1_dtevp
[921]	510	! SB modif because ocean has no slip boundary condition
[5123]	511	zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) &
	512	& +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) &
	513	& / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
[825]	514
[5123]	515	za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + &
	516	& ( v_ice(ji,jj) - v_oce(ji,jj) )*2 ) ( 1.0 - zfrld2(ji,jj) )
	517	zr = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj)
	518	zcca = z0 + za
[4990]	519	zccb = zcorl2(ji,jj)
[5123]	520	v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch
[921]	521	END DO
	522	END DO
[825]	523
[921]	524	CALL lbc_lnk( v_ice(:,:), 'V', -1. )
[4161]	525	#if defined key_agrif && defined key_lim2
[5123]	526	CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
[3680]	527	#endif
[6584]	528	#if defined key_agrif && defined key_lim3
	529	CALL agrif_interp_lim3( 'V' )
	530	#endif
[4333]	531	#if defined key_bdy
	532	CALL bdy_ice_lim_dyn( 'V' )
	533	#endif
[825]	534
[921]	535	DO jj = k_j1+1, k_jpj-1
	536	DO ji = fs_2, fs_jpim1
[5123]	537	rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1)
	538	z0 = zmass1(ji,jj) * z1_dtevp
	539	zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji,jj) &
	540	& + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) &
	541	& / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
[825]	542
[5123]	543	za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + &
	544	& ( zv_ice1 - v_oce1(ji,jj) )*2 ) ( 1.0 - zfrld1(ji,jj) )
	545	zr = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj)
	546	zcca = z0 + za
[4990]	547	zccb = zcorl1(ji,jj)
[5123]	548	u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch
	549	END DO
	550	END DO
[825]	551
[921]	552	CALL lbc_lnk( u_ice(:,:), 'U', -1. )
[4161]	553	#if defined key_agrif && defined key_lim2
[5123]	554	CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
[3680]	555	#endif
[6584]	556	#if defined key_agrif && defined key_lim3
	557	CALL agrif_interp_lim3( 'U' )
	558	#endif
[4333]	559	#if defined key_bdy
	560	CALL bdy_ice_lim_dyn( 'U' )
	561	#endif
[825]	562
[921]	563	ENDIF
[4161]	564
[921]	565	IF(ln_ctl) THEN
	566	!--- Convergence test.
	567	DO jj = k_j1+1 , k_jpj-1
[5123]	568	zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
[921]	569	END DO
[5123]	570	zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) )
[921]	571	IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain
	572	ENDIF
	573
[4161]	574	! ! ==================== !
[868]	575	END DO ! end loop over jter !
[825]	576	! ! ==================== !
[921]	577	!
	578	!------------------------------------------------------------------------------!
	579	! 4) Prevent ice velocities when the ice is thin
	580	!------------------------------------------------------------------------------!
[5123]	581	! If the ice volume is below zvmin then ice velocity should equal the
	582	! ocean velocity. This prevents high velocity when ice is thin
[825]	583	DO jj = k_j1+1, k_jpj-1
[868]	584	DO ji = fs_2, fs_jpim1
[5123]	585	IF ( vt_i(ji,jj) <= zvmin ) THEN
[888]	586	u_ice(ji,jj) = u_oce(ji,jj)
	587	v_ice(ji,jj) = v_oce(ji,jj)
[5123]	588	ENDIF
[825]	589	END DO
	590	END DO
[866]	591
[5429]	592	CALL lbc_lnk_multi( u_ice(:,:), 'U', -1., v_ice(:,:), 'V', -1. )
	593
[4161]	594	#if defined key_agrif && defined key_lim2
[5123]	595	CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' )
	596	CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' )
[3680]	597	#endif
[6584]	598	#if defined key_agrif && defined key_lim3
	599	CALL agrif_interp_lim3( 'U' )
	600	CALL agrif_interp_lim3( 'V' )
	601	#endif
[4161]	602	#if defined key_bdy
[4333]	603	CALL bdy_ice_lim_dyn( 'U' )
	604	CALL bdy_ice_lim_dyn( 'V' )
[4161]	605	#endif
[869]	606
[868]	607	DO jj = k_j1+1, k_jpj-1
	608	DO ji = fs_2, fs_jpim1
[5123]	609	IF ( vt_i(ji,jj) <= zvmin ) THEN
	610	v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji, jj-1) ) * e1t(ji+1,jj) &
	611	& + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) &
	612	& / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
[868]	613
[5123]	614	u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) &
	615	& + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) &
	616	& / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
	617	ENDIF
[868]	618	END DO
	619	END DO
	620
[5429]	621	CALL lbc_lnk_multi( u_ice2(:,:), 'V', -1., v_ice1(:,:), 'U', -1. )
[868]	622
[866]	623	! Recompute delta, shear and div, inputs for mechanical redistribution
[825]	624	DO jj = k_j1+1, k_jpj-1
[5123]	625	DO ji = fs_2, jpim1 !RB bug no vect opt due to zmask
[4990]	626	!- divu_i(:,:), zdt(:,:): divergence and tension at centre
[825]	627	!- zds(:,:): shear on northeast corner of grid cells
[5123]	628	IF ( vt_i(ji,jj) <= zvmin ) THEN
[825]	629
[5123]	630	divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj ) * u_ice(ji-1,jj ) &
	631	& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji ,jj-1) * v_ice(ji ,jj-1) &
	632	& ) * r1_e12t(ji,jj)
[825]	633
[5123]	634	zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) &
	635	& -( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
	636	& ) * r1_e12t(ji,jj)
[921]	637	!
	638	! SB modif because ocean has no slip boundary condition
[5123]	639	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) &
	640	& +( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
	641	& ) * r1_e12f(ji,jj) * ( 2.0 - fmask(ji,jj,1) ) &
	642	& * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1)
[825]	643
[5123]	644	zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) &
	645	& + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji ,jj-1) * u_ice2(ji ,jj-1) ) * r1_e12t(ji,jj)
[825]	646
[5123]	647	delta = SQRT( divu_i(ji,jj)2 + ( zdt(ji,jj)2 + zdst*2 ) usecc2 )
	648	delta_i(ji,jj) = delta + rn_creepl
[4161]	649
[5123]	650	ENDIF
	651	END DO
	652	END DO
[921]	653	!
	654	!------------------------------------------------------------------------------!
	655	! 5) Store stress tensor and its invariants
	656	!------------------------------------------------------------------------------!
[866]	657	! * Invariants of the stress tensor are required for limitd_me
[3791]	658	! (accelerates convergence and improves stability)
[866]	659	DO jj = k_j1+1, k_jpj-1
[868]	660	DO ji = fs_2, fs_jpim1
[5123]	661	zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj ) * v_ice1(ji-1,jj) &
	662	& + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e12t(ji,jj)
[4990]	663	shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst )
[825]	664	END DO
[866]	665	END DO
[4161]	666
	667	! Lateral boundary condition
[5429]	668	CALL lbc_lnk_multi( divu_i (:,:), 'T', 1., delta_i(:,:), 'T', 1., shear_i(:,:), 'T', 1. )
[866]	669
[868]	670	! * Store the stress tensor for the next time step
	671	stress1_i (:,:) = zs1 (:,:)
	672	stress2_i (:,:) = zs2 (:,:)
	673	stress12_i(:,:) = zs12(:,:)
	674
[921]	675	!
	676	!------------------------------------------------------------------------------!
	677	! 6) Control prints of residual and charge ellipse
	678	!------------------------------------------------------------------------------!
	679	!
[834]	680	! print the residual for convergence
	681	IF(ln_ctl) THEN
[868]	682	WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, jter
[834]	683	CALL prt_ctl_info(charout)
	684	CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
	685	ENDIF
[825]	686
[834]	687	! print charge ellipse
	688	! This can be desactivated once the user is sure that the stress state
	689	! lie on the charge ellipse. See Bouillon et al. 08 for more details
[825]	690	IF(ln_ctl) THEN
	691	CALL prt_ctl_info('lim_rhg : numit :',ivar1=numit)
	692	CALL prt_ctl_info('lim_rhg : nwrite :',ivar1=nwrite)
	693	CALL prt_ctl_info('lim_rhg : MOD :',ivar1=MOD(numit,nwrite))
	694	IF( MOD(numit,nwrite) .EQ. 0 ) THEN
	695	WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. '
	696	CALL prt_ctl_info(charout)
	697	DO jj = k_j1+1, k_jpj-1
	698	DO ji = 2, jpim1
[5123]	699	IF (zpresh(ji,jj) > 1.0) THEN
[825]	700	sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)*2 + 4zs12(ji,jj)2 )0.5 ) / ( 2*zpresh(ji,jj) )
	701	sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)*2 + 4zs12(ji,jj)2 )0.5 ) / ( 2*zpresh(ji,jj) )
	702	WRITE(charout,FMT="('lim_rhg :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
	703	CALL prt_ctl_info(charout)
	704	ENDIF
	705	END DO
	706	END DO
	707	WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 2000, numit, 0, 0, ' ch. ell. '
	708	CALL prt_ctl_info(charout)
	709	ENDIF
	710	ENDIF
[2715]	711	!
[3294]	712	CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
[5123]	713	CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask )
[4990]	714	CALL wrk_dealloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds )
[5888]	715	CALL wrk_dealloc( jpi,jpj, zs1 , zs2 , zs12 , zresr , zpice )
[3294]	716
[825]	717	END SUBROUTINE lim_rhg
	718
	719	#else
	720	!!----------------------------------------------------------------------
	721	!! Default option Dummy module NO LIM sea-ice model
	722	!!----------------------------------------------------------------------
	723	CONTAINS
	724	SUBROUTINE lim_rhg( k1 , k2 ) ! Dummy routine
	725	WRITE(,) 'lim_rhg: You should not have seen this print! error?', k1, k2
	726	END SUBROUTINE lim_rhg
	727	#endif
	728
	729	!!==============================================================================
	730	END MODULE limrhg

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

Download in other formats: