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

Context Navigation

source: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 3634

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

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

Download in other formats: