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

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 4444

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

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

Download in other formats: