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

Context Navigation

source: branches/UKMO/DEV_r5107_dynvor_updates/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 5256

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

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

Download in other formats: