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

Context Navigation

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 5059

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

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

Download in other formats: