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

Context Navigation

source: branches/2011/dev_r2802_LOCEAN10_agrif_lim/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 2804

Last change on this file since 2804 was 2804, checked in by rblod, 13 years ago
dev_r2802_LOCEAN10_agrif_lim: first implementation see ticket #848
Property svn:keywords set to `Id`
File size: 38.9 KB

Rev	Line
[825]	1	MODULE limrhg
	2	!!======================================================================
	3	!! * MODULE limrhg *
[834]	4	!! Ice rheology : sea ice rheology
[825]	5	!!======================================================================
[1244]	6	!! History : - ! 2007-03 (M.A. Morales Maqueda, S. Bouillon) Original code
	7	!! 3.0 ! 2008-03 (M. Vancoppenolle) LIM3
	8	!! - ! 2008-11 (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy
[2528]	9	!! 3.3 ! 2009-05 (G.Garric) addition of the lim2_evp cas
[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
[2804]	36	#if defined key_agrif && defined key_lim2
	37	USE agrif_lim2_interp
	38	#endif
[825]	39
	40	IMPLICIT NONE
	41	PRIVATE
	42
[2715]	43	PUBLIC lim_rhg ! routine called by lim_dyn (or lim_dyn_2)
	44	PUBLIC lim_rhg_alloc ! routine called by nemo_alloc in nemogcm.F90
[825]	45
[2528]	46	REAL(wp) :: rzero = 0._wp ! constant values
	47	REAL(wp) :: rone = 1._wp ! constant values
	48
[2715]	49	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zpresh ! temporary array for ice strength
	50	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zpreshc ! Ice strength on grid cell corners (zpreshc)
	51	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfrld1, zfrld2 ! lead fraction on U/V points
	52	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zmass1, zmass2 ! ice/snow mass on U/V points
	53	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zcorl1, zcorl2 ! coriolis parameter on U/V points
	54	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: za1ct , za2ct ! temporary arrays
	55	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zc1 ! ice mass
	56	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zusw ! temporary weight for ice strength computation
	57	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_oce1, v_oce1 ! ocean u/v component on U points
	58	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_oce2, v_oce2 ! ocean u/v component on V points
	59	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_ice2, v_ice1 ! ice u/v component on V/U point
	60	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zf1 , zf2 ! arrays for internal stresses
	61
	62	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zdd , zdt ! Divergence and tension at centre of grid cells
	63	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zds ! Shear on northeast corner of grid cells
	64	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: deltat, deltac ! Delta at centre and corners of grid cells
	65	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zs1 , zs2 ! Diagonal stress tensor components zs1 and zs2
	66	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zs12 ! Non-diagonal stress tensor component zs12
	67	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zu_ice, zv_ice, zresr ! Local error on velocity
	68
[868]	69	!! * Substitutions
	70	# include "vectopt_loop_substitute.h90"
[825]	71	!!----------------------------------------------------------------------
[2715]	72	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]	73	!! $Id$
[2528]	74	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	75	!!----------------------------------------------------------------------
	76	CONTAINS
	77
[2715]	78	FUNCTION lim_rhg_alloc()
	79	!!-------------------------------------------------------------------
	80	!! * FUNCTION lim_rhg_alloc *
	81	!!-------------------------------------------------------------------
	82	INTEGER :: lim_rhg_alloc ! return value
	83	INTEGER :: ierr(2) ! local integer
	84	!!-------------------------------------------------------------------
	85	!
	86	ierr(:) = 0
	87	!
	88	ALLOCATE( zpresh (jpi,jpj) , zfrld1(jpi,jpj), zmass1(jpi,jpj), zcorl1(jpi,jpj), za1ct(jpi,jpj) , &
	89	& zpreshc(jpi,jpj) , zfrld2(jpi,jpj), zmass2(jpi,jpj), zcorl2(jpi,jpj), za2ct(jpi,jpj) , &
	90	& zc1 (jpi,jpj) , u_oce1(jpi,jpj), u_oce2(jpi,jpj), u_ice2(jpi,jpj), &
	91	& zusw (jpi,jpj) , v_oce1(jpi,jpj), v_oce2(jpi,jpj), v_ice1(jpi,jpj) , STAT=ierr(1) )
	92	!
	93	ALLOCATE( zf1(jpi,jpj) , deltat(jpi,jpj) , zu_ice(jpi,jpj) , &
	94	& zf2(jpi,jpj) , deltac(jpi,jpj) , zv_ice(jpi,jpj) , &
	95	& zdd(jpi,jpj) , zdt (jpi,jpj) , zds (jpi,jpj) , &
	96	& zs1(jpi,jpj) , zs2 (jpi,jpj) , zs12 (jpi,jpj) , zresr(jpi,jpj), STAT=ierr(2) )
[2804]	97	!
[2715]	98	lim_rhg_alloc = MAXVAL(ierr)
	99	!
	100	END FUNCTION lim_rhg_alloc
	101
	102
[825]	103	SUBROUTINE lim_rhg( k_j1, k_jpj )
	104	!!-------------------------------------------------------------------
[834]	105	!! * SUBROUTINE lim_rhg *
	106	!! EVP-C-grid
[825]	107	!!
[834]	108	!! ** purpose : determines sea ice drift from wind stress, ice-ocean
[825]	109	!! stress and sea-surface slope. Ice-ice interaction is described by
[834]	110	!! a non-linear elasto-viscous-plastic (EVP) law including shear
	111	!! strength and a bulk rheology (Hunke and Dukowicz, 2002).
[825]	112	!!
[834]	113	!! The points in the C-grid look like this, dear reader
[825]	114	!!
[834]	115	!! (ji,jj)
	116	!! \|
	117	!! \|
	118	!! (ji-1,jj) \| (ji,jj)
	119	!! ---------
	120	!! \| \|
	121	!! \| (ji,jj) \|------(ji,jj)
	122	!! \| \|
	123	!! ---------
	124	!! (ji-1,jj-1) (ji,jj-1)
[825]	125	!!
[834]	126	!! ** Inputs : - wind forcing (stress), oceanic currents
	127	!! ice total volume (vt_i) per unit area
	128	!! snow total volume (vt_s) per unit area
[825]	129	!!
[834]	130	!! ** Action : - compute u_ice, v_ice : the components of the
	131	!! sea-ice velocity vector
	132	!! - compute delta_i, shear_i, divu_i, which are inputs
	133	!! of the ice thickness distribution
[825]	134	!!
[834]	135	!! ** Steps : 1) Compute ice snow mass, ice strength
	136	!! 2) Compute wind, oceanic stresses, mass terms and
	137	!! coriolis terms of the momentum equation
	138	!! 3) Solve the momentum equation (iterative procedure)
	139	!! 4) Prevent high velocities if the ice is thin
	140	!! 5) Recompute invariants of the strain rate tensor
	141	!! which are inputs of the ITD, store stress
	142	!! for the next time step
	143	!! 6) Control prints of residual (convergence)
	144	!! and charge ellipse.
	145	!! The user should make sure that the parameters
	146	!! nevp, telast and creepl maintain stress state
	147	!! on the charge ellipse for plastic flow
	148	!! e.g. in the Canadian Archipelago
	149	!!
[2528]	150	!! References : Hunke and Dukowicz, JPO97
	151	!! Bouillon et al., Ocean Modelling 2009
	152	!! Vancoppenolle et al., Ocean Modelling 2008
	153	!!-------------------------------------------------------------------
	154	INTEGER, INTENT(in) :: k_j1 ! southern j-index for ice computation
	155	INTEGER, INTENT(in) :: k_jpj ! northern j-index for ice computation
[834]	156	!!
[2528]	157	INTEGER :: ji, jj ! dummy loop indices
	158	INTEGER :: jter ! local integers
[825]	159	CHARACTER (len=50) :: charout
[2528]	160	REAL(wp) :: zt11, zt12, zt21, zt22, ztagnx, ztagny, delta !
	161	REAL(wp) :: za, zstms, zsang, zmask ! local scalars
[825]	162
[2715]	163	REAL(wp) :: dtevp ! time step for subcycling
	164	REAL(wp) :: dtotel, ecc2 ! square of yield ellipse eccenticity
	165	REAL(wp) :: z0, zr, zcca, zccb ! temporary scalars
	166	REAL(wp) :: zu_ice2, zv_ice1 !
	167	REAL(wp) :: zddc, zdtc, zdst ! delta on corners and on centre
	168	REAL(wp) :: zdsshx, zdsshy ! term for the gradient of ocean surface
	169	REAL(wp) :: sigma1, sigma2 ! internal ice stress
[825]	170
[2715]	171	REAL(wp) :: zresm ! Maximal error on ice velocity
	172	REAL(wp) :: zindb ! ice (1) or not (0)
	173	REAL(wp) :: zdummy ! dummy argument
[2528]	174	!!-------------------------------------------------------------------
	175	#if defined key_lim2 && ! defined key_lim2_vp
	176	# if defined key_agrif
	177	USE ice_2, vt_s => hsnm
	178	USE ice_2, vt_i => hicm
	179	# else
	180	vt_s => hsnm
	181	vt_i => hicm
	182	# endif
	183	at_i(:,:) = 1. - frld(:,:)
	184	#endif
[921]	185	!
	186	!------------------------------------------------------------------------------!
	187	! 1) Ice-Snow mass (zc1), ice strength (zpresh) !
	188	!------------------------------------------------------------------------------!
	189	!
[825]	190	! Put every vector to 0
[2528]	191	zpresh (:,:) = 0._wp ; zc1 (:,:) = 0._wp
	192	zpreshc(:,:) = 0._wp
	193	u_ice2 (:,:) = 0._wp ; v_ice1(:,:) = 0._wp
	194	zdd (:,:) = 0._wp ; zdt (:,:) = 0._wp ; zds(:,:) = 0._wp
[825]	195
[2528]	196	#if defined key_lim3
	197	CALL lim_itd_me_icestrength( ridge_scheme_swi ) ! LIM-3: Ice strength on T-points
	198	#endif
[825]	199
[868]	200	!CDIR NOVERRCHK
[2528]	201	DO jj = k_j1 , k_jpj ! Ice mass and temp variables
[868]	202	!CDIR NOVERRCHK
[825]	203	DO ji = 1 , jpi
	204	zc1(ji,jj) = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) )
[2528]	205	#if defined key_lim3
	206	zpresh(ji,jj) = tms(ji,jj) * strength(ji,jj) * 0.5_wp
	207	#endif
[2580]	208	#if defined key_lim2
	209	zpresh(ji,jj) = tms(ji,jj) * pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) )
	210	#endif
[866]	211	! tmi = 1 where there is ice or on land
[2528]	212	tmi(ji,jj) = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - epsd ) ) ) * tms(ji,jj)
[825]	213	END DO
	214	END DO
	215
[834]	216	! Ice strength on grid cell corners (zpreshc)
	217	! needed for calculation of shear stress
[868]	218	!CDIR NOVERRCHK
[825]	219	DO jj = k_j1+1, k_jpj-1
[868]	220	!CDIR NOVERRCHK
	221	DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1)
[921]	222	zstms = tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
	223	& tms(ji,jj+1) * wght(ji+1,jj+1,1,2) + &
	224	& tms(ji+1,jj) * wght(ji+1,jj+1,2,1) + &
	225	& tms(ji,jj) * wght(ji+1,jj+1,1,1)
	226	zusw(ji,jj) = 1.0 / MAX( zstms, epsd )
	227	zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
	228	& zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + &
	229	& zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + &
	230	& zpresh(ji,jj) * wght(ji+1,jj+1,1,1) &
	231	& ) * zusw(ji,jj)
[825]	232	END DO
	233	END DO
	234
	235	CALL lbc_lnk( zpreshc(:,:), 'F', 1. )
[921]	236	!
	237	!------------------------------------------------------------------------------!
	238	! 2) Wind / ocean stress, mass terms, coriolis terms
	239	!------------------------------------------------------------------------------!
	240	!
[825]	241	! Wind stress, coriolis and mass terms on the sides of the squares
	242	! zfrld1: lead fraction on U-points
	243	! zfrld2: lead fraction on V-points
	244	! zmass1: ice/snow mass on U-points
	245	! zmass2: ice/snow mass on V-points
	246	! zcorl1: Coriolis parameter on U-points
	247	! zcorl2: Coriolis parameter on V-points
	248	! (ztagnx,ztagny): wind stress on U/V points
	249	! u_oce1: ocean u component on u points
	250	! v_oce1: ocean v component on u points
	251	! u_oce2: ocean u component on v points
	252	! v_oce2: ocean v component on v points
[921]	253
[825]	254	DO jj = k_j1+1, k_jpj-1
[868]	255	DO ji = fs_2, fs_jpim1
[825]	256
[2528]	257	zt11 = tms(ji ,jj) * e1t(ji ,jj)
	258	zt12 = tms(ji+1,jj) * e1t(ji+1,jj)
	259	zt21 = tms(ji,jj ) * e2t(ji,jj )
	260	zt22 = tms(ji,jj+1) * e2t(ji,jj+1)
[825]	261
	262	! Leads area.
[2528]	263	zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd )
	264	zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd )
[825]	265
	266	! Mass, coriolis coeff. and currents
[834]	267	zmass1(ji,jj) = ( zt12zc1(ji,jj) + zt11zc1(ji+1,jj) ) / (zt11+zt12+epsd)
	268	zmass2(ji,jj) = ( zt22zc1(ji,jj) + zt21zc1(ji,jj+1) ) / (zt21+zt22+epsd)
[2528]	269	zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)fcor(ji,jj) + e1t(ji,jj)fcor(ji+1,jj) ) &
	270	& / ( e1t(ji,jj) + e1t(ji+1,jj) + epsd )
	271	zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)fcor(ji,jj) + e2t(ji,jj)fcor(ji,jj+1) ) &
	272	& / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd )
[825]	273	!
[888]	274	u_oce1(ji,jj) = u_oce(ji,jj)
	275	v_oce2(ji,jj) = v_oce(ji,jj)
[825]	276
[834]	277	! Ocean has no slip boundary condition
[888]	278	v_oce1(ji,jj) = 0.5( (v_oce(ji,jj)+v_oce(ji,jj-1))e1t(ji,jj) &
[921]	279	& +(v_oce(ji+1,jj)+v_oce(ji+1,jj-1))*e1t(ji+1,jj)) &
	280	& /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
[825]	281
[888]	282	u_oce2(ji,jj) = 0.5((u_oce(ji,jj)+u_oce(ji-1,jj))e2t(ji,jj) &
[921]	283	& +(u_oce(ji,jj+1)+u_oce(ji-1,jj+1))*e2t(ji,jj+1)) &
	284	& / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
[825]	285
[1469]	286	! Wind stress at U,V-point
	287	ztagnx = ( 1. - zfrld1(ji,jj) ) * utau_ice(ji,jj)
	288	ztagny = ( 1. - zfrld2(ji,jj) ) * vtau_ice(ji,jj)
[825]	289
[834]	290	! Computation of the velocity field taking into account the ice internal interaction.
[825]	291	! Terms that are independent of the velocity field.
	292
	293	! SB On utilise maintenant le gradient de la pente de l'ocean
	294	! include it later
[834]	295
[2528]	296	zdsshx = ( ssh_m(ji+1,jj) - ssh_m(ji,jj) ) / e1u(ji,jj)
	297	zdsshy = ( ssh_m(ji,jj+1) - ssh_m(ji,jj) ) / e2v(ji,jj)
[825]	298
	299	za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx
	300	za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy
	301
	302	END DO
	303	END DO
	304
[921]	305	!
	306	!------------------------------------------------------------------------------!
	307	! 3) Solution of the momentum equation, iterative procedure
	308	!------------------------------------------------------------------------------!
	309	!
[825]	310	! Time step for subcycling
	311	dtevp = rdt_ice / nevp
[2528]	312	dtotel = dtevp / ( 2._wp * telast )
[825]	313
	314	!-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)
[2528]	315	ecc2 = ecc * ecc
[825]	316
	317	!-Initialise stress tensor
[2528]	318	zs1 (:,:) = stress1_i (:,:)
	319	zs2 (:,:) = stress2_i (:,:)
[866]	320	zs12(:,:) = stress12_i(:,:)
[825]	321
[2528]	322	! !----------------------!
[868]	323	DO jter = 1 , nevp ! loop over jter !
[2528]	324	! !----------------------!
[825]	325	DO jj = k_j1, k_jpj-1
[2528]	326	zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
[825]	327	zv_ice(:,jj) = v_ice(:,jj)
[921]	328	END DO
[825]	329
[834]	330	DO jj = k_j1+1, k_jpj-1
[988]	331	DO ji = fs_2, jpim1 !RB bug no vect opt due to tmi
[825]	332
[921]	333	!
	334	!- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002)
	335	!- zdd(:,:), zdt(:,:): divergence and tension at centre of grid cells
	336	!- zds(:,:): shear on northeast corner of grid cells
	337	!
	338	!- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded,
	339	! there are many repeated calculations.
	340	! Speed could be improved by regrouping terms. For
	341	! the moment, however, the stress is on clarity of coding to avoid
	342	! bugs (Martin, for Miguel).
	343	!
	344	!- ALSO: arrays zdd, zdt, zds and delta could
	345	! be removed in the future to minimise memory demand.
	346	!
	347	!- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of
	348	! grid cells, exactly as in the B grid case. For simplicity, the indexation on
	349	! the corners is the same as in the B grid.
	350	!
	351	!
	352	zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) &
	353	& -e2u(ji-1,jj)*u_ice(ji-1,jj) &
	354	& +e1v(ji,jj)*v_ice(ji,jj) &
	355	& -e1v(ji,jj-1)*v_ice(ji,jj-1) &
	356	& ) &
	357	& / area(ji,jj)
[825]	358
[921]	359	zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) &
	360	& -u_ice(ji-1,jj)/e2u(ji-1,jj) &
	361	& )e2t(ji,jj)e2t(ji,jj) &
	362	& -( v_ice(ji,jj)/e1v(ji,jj) &
	363	& -v_ice(ji,jj-1)/e1v(ji,jj-1) &
	364	& )e1t(ji,jj)e1t(ji,jj) &
	365	& ) &
	366	& / area(ji,jj)
[825]	367
[921]	368	!
	369	zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1) &
	370	& -u_ice(ji,jj)/e1u(ji,jj) &
	371	& )e1f(ji,jj)e1f(ji,jj) &
	372	& +( v_ice(ji+1,jj)/e2v(ji+1,jj) &
	373	& -v_ice(ji,jj)/e2v(ji,jj) &
	374	& )e2f(ji,jj)e2f(ji,jj) &
	375	& ) &
	376	& / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) &
	377	& * tmi(ji,jj) * tmi(ji,jj+1) &
	378	& * tmi(ji+1,jj) * tmi(ji+1,jj+1)
[825]	379
	380
[921]	381	v_ice1(ji,jj) = 0.5( (v_ice(ji,jj)+v_ice(ji,jj-1))e1t(ji+1,jj) &
	382	& +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
	383	& /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
[825]	384
[921]	385	u_ice2(ji,jj) = 0.5( (u_ice(ji,jj)+u_ice(ji-1,jj))e2t(ji,jj+1) &
	386	& +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &
	387	& /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
[825]	388
[921]	389	END DO
	390	END DO
[2528]	391	CALL lbc_lnk( v_ice1, 'U', -1. ) ; CALL lbc_lnk( u_ice2, 'V', -1. ) ! lateral boundary cond.
[921]	392
[868]	393	!CDIR NOVERRCHK
[921]	394	DO jj = k_j1+1, k_jpj-1
[868]	395	!CDIR NOVERRCHK
[921]	396	DO ji = fs_2, fs_jpim1
[825]	397
[921]	398	!- Calculate Delta at centre of grid cells
[2528]	399	zdst = ( e2u(ji , jj) * v_ice1(ji ,jj) &
	400	& - e2u(ji-1, jj) * v_ice1(ji-1,jj) &
	401	& + e1v(ji, jj ) * u_ice2(ji,jj ) &
	402	& - e1v(ji, jj-1) * u_ice2(ji,jj-1) &
	403	& ) &
[921]	404	& / area(ji,jj)
[825]	405
[2528]	406	delta = SQRT( zdd(ji,jj)zdd(ji,jj) + ( zdt(ji,jj)zdt(ji,jj) + zdstzdst ) usecc2 )
	407	deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)2 + (zdt(ji,jj)2 + zdst*2)usecc2), creepl )
[825]	408
[921]	409	!-Calculate stress tensor components zs1 and zs2
	410	!-at centre of grid cells (see section 3.5 of CICE user's guide).
	411	zs1(ji,jj) = ( zs1(ji,jj) &
	412	& - dtotel( ( 1.0 - alphaevp) zs1(ji,jj) + &
	413	& ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) &
	414	* zpresh(ji,jj) ) ) &
	415	& / ( 1.0 + alphaevp * dtotel )
[825]	416
[921]	417	zs2(ji,jj) = ( zs2(ji,jj) &
	418	& - dtotel((1.0-alphaevp)ecc2*zs2(ji,jj) - &
	419	zdt(ji,jj)/deltat(ji,jj)*zpresh(ji,jj)) ) &
	420	& / ( 1.0 + alphaevpecc2dtotel )
[825]	421
[921]	422	END DO
	423	END DO
[825]	424
[921]	425	CALL lbc_lnk( zs1(:,:), 'T', 1. )
	426	CALL lbc_lnk( zs2(:,:), 'T', 1. )
[825]	427
[868]	428	!CDIR NOVERRCHK
[921]	429	DO jj = k_j1+1, k_jpj-1
[868]	430	!CDIR NOVERRCHK
[921]	431	DO ji = fs_2, fs_jpim1
	432	!- Calculate Delta on corners
	433	zddc = ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1) &
	434	& -v_ice1(ji,jj)/e1u(ji,jj) &
	435	& )e1f(ji,jj)e1f(ji,jj) &
	436	& +( u_ice2(ji+1,jj)/e2v(ji+1,jj) &
	437	& -u_ice2(ji,jj)/e2v(ji,jj) &
	438	& )e2f(ji,jj)e2f(ji,jj) &
	439	& ) &
	440	& / ( e1f(ji,jj) * e2f(ji,jj) )
[825]	441
[921]	442	zdtc = (-( v_ice1(ji,jj+1)/e1u(ji,jj+1) &
	443	& -v_ice1(ji,jj)/e1u(ji,jj) &
	444	& )e1f(ji,jj)e1f(ji,jj) &
	445	& +( u_ice2(ji+1,jj)/e2v(ji+1,jj) &
	446	& -u_ice2(ji,jj)/e2v(ji,jj) &
	447	& )e2f(ji,jj)e2f(ji,jj) &
	448	& ) &
	449	& / ( e1f(ji,jj) * e2f(ji,jj) )
[825]	450
[921]	451	deltac(ji,jj) = SQRT(zddc2+(zdtc2+zds(ji,jj)*2)usecc2) + creepl
[825]	452
[921]	453	!-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide).
	454	zs12(ji,jj) = ( zs12(ji,jj) &
	455	& - dtotel( (1.0-alphaevp)ecc2*zs12(ji,jj) - zds(ji,jj) / &
	456	& ( 2.0deltac(ji,jj) ) zpreshc(ji,jj))) &
	457	& / ( 1.0 + alphaevpecc2dtotel )
[825]	458
[921]	459	END DO ! ji
	460	END DO ! jj
[825]	461
[921]	462	CALL lbc_lnk( zs12(:,:), 'F', 1. )
[825]	463
[921]	464	! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002)
	465	DO jj = k_j1+1, k_jpj-1
	466	DO ji = fs_2, fs_jpim1
	467	!- contribution of zs1, zs2 and zs12 to zf1
	468	zf1(ji,jj) = 0.5( (zs1(ji+1,jj)-zs1(ji,jj))e2u(ji,jj) &
	469	& +(zs2(ji+1,jj)e2t(ji+1,jj)2-zs2(ji,jj)e2t(ji,jj)**2)/e2u(ji,jj) &
	470	& +2.0(zs12(ji,jj)e1f(ji,jj)*2-zs12(ji,jj-1)e1f(ji,jj-1)**2)/e1u(ji,jj) &
	471	& ) / ( e1u(ji,jj)*e2u(ji,jj) )
	472	! contribution of zs1, zs2 and zs12 to zf2
	473	zf2(ji,jj) = 0.5( (zs1(ji,jj+1)-zs1(ji,jj))e1v(ji,jj) &
	474	& -(zs2(ji,jj+1)e1t(ji,jj+1)2 - zs2(ji,jj)e1t(ji,jj)**2)/e1v(ji,jj) &
	475	& + 2.0(zs12(ji,jj)e2f(ji,jj)**2 - &
	476	zs12(ji-1,jj)e2f(ji-1,jj)*2)/e2v(ji,jj) &
	477	& ) / ( e1v(ji,jj)*e2v(ji,jj) )
	478	END DO
	479	END DO
[825]	480	!
	481	! Computation of ice velocity
	482	!
	483	! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly.
	484	!
[921]	485	IF (MOD(jter,2).eq.0) THEN
[825]	486
[868]	487	!CDIR NOVERRCHK
[921]	488	DO jj = k_j1+1, k_jpj-1
[868]	489	!CDIR NOVERRCHK
[921]	490	DO ji = fs_2, fs_jpim1
	491	zmask = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj)
	492	zsang = SIGN ( 1.0 , fcor(ji,jj) ) * sangvg
	493	z0 = zmass1(ji,jj)/dtevp
[825]	494
[921]	495	! SB modif because ocean has no slip boundary condition
	496	zv_ice1 = 0.5( (v_ice(ji,jj)+v_ice(ji,jj-1))e1t(ji,jj) &
	497	& +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) &
	498	& /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
	499	za = rhocoSQRT((u_ice(ji,jj)-u_oce1(ji,jj))*2 + &
	500	(zv_ice1-v_oce1(ji,jj))*2) (1.0-zfrld1(ji,jj))
	501	zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
	502	za(cangvgu_oce1(ji,jj)-zsang*v_oce1(ji,jj))
	503	zcca = z0+za*cangvg
	504	zccb = zcorl1(ji,jj)+za*zsang
	505	u_ice(ji,jj) = (zr+zccbzv_ice1)/(zcca+epsd)zmask
[825]	506
[921]	507	END DO
	508	END DO
[825]	509
[921]	510	CALL lbc_lnk( u_ice(:,:), 'U', -1. )
[2804]	511	#if defined key_agrif
	512	CALL agrif_dyn_lim(jter ,nevp , 'U')
	513	#endif
[868]	514	!CDIR NOVERRCHK
[921]	515	DO jj = k_j1+1, k_jpj-1
[868]	516	!CDIR NOVERRCHK
[921]	517	DO ji = fs_2, fs_jpim1
[834]	518
[921]	519	zmask = (1.0-MAX(rzero,SIGN(rone,-zmass2(ji,jj))))*tmv(ji,jj)
	520	zsang = SIGN(1.0,fcor(ji,jj))*sangvg
	521	z0 = zmass2(ji,jj)/dtevp
	522	! SB modif because ocean has no slip boundary condition
	523	zu_ice2 = 0.5( (u_ice(ji,jj)+u_ice(ji-1,jj))e2t(ji,jj) &
	524	& + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) &
	525	& /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
	526	za = rhocoSQRT((zu_ice2-u_oce2(ji,jj))*2 + &
	527	(v_ice(ji,jj)-v_oce2(ji,jj))*2)(1.0-zfrld2(ji,jj))
	528	zr = z0*v_ice(ji,jj) + zf2(ji,jj) + &
	529	za2ct(ji,jj) + za(cangvgv_oce2(ji,jj)+zsang*u_oce2(ji,jj))
	530	zcca = z0+za*cangvg
	531	zccb = zcorl2(ji,jj)+za*zsang
	532	v_ice(ji,jj) = (zr-zccbzu_ice2)/(zcca+epsd)zmask
[825]	533
[921]	534	END DO
	535	END DO
[825]	536
[921]	537	CALL lbc_lnk( v_ice(:,:), 'V', -1. )
[2804]	538	#if defined key_agrif
	539	CALL agrif_dyn_lim(jter ,nevp , 'U')
	540	#endif
[834]	541	ELSE
[868]	542	!CDIR NOVERRCHK
[921]	543	DO jj = k_j1+1, k_jpj-1
[868]	544	!CDIR NOVERRCHK
[921]	545	DO ji = fs_2, fs_jpim1
	546	zmask = (1.0-MAX(rzero,SIGN(rone,-zmass2(ji,jj))))*tmv(ji,jj)
	547	zsang = SIGN(1.0,fcor(ji,jj))*sangvg
	548	z0 = zmass2(ji,jj)/dtevp
	549	! SB modif because ocean has no slip boundary condition
	550	zu_ice2 = 0.5( (u_ice(ji,jj)+u_ice(ji-1,jj))e2t(ji,jj) &
	551	& +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) &
	552	& /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
[825]	553
[921]	554	za = rhocoSQRT((zu_ice2-u_oce2(ji,jj))*2 + &
	555	(v_ice(ji,jj)-v_oce2(ji,jj))*2)(1.0-zfrld2(ji,jj))
	556	zr = z0*v_ice(ji,jj) + zf2(ji,jj) + &
	557	za2ct(ji,jj) + za(cangvgv_oce2(ji,jj)+zsang*u_oce2(ji,jj))
	558	zcca = z0+za*cangvg
	559	zccb = zcorl2(ji,jj)+za*zsang
	560	v_ice(ji,jj) = (zr-zccbzu_ice2)/(zcca+epsd)zmask
[825]	561
[921]	562	END DO
	563	END DO
[825]	564
[921]	565	CALL lbc_lnk( v_ice(:,:), 'V', -1. )
[2804]	566	#if defined key_agrif
	567	CALL agrif_dyn_lim(jter ,nevp , 'V')
	568	#endif
[868]	569	!CDIR NOVERRCHK
[921]	570	DO jj = k_j1+1, k_jpj-1
[868]	571	!CDIR NOVERRCHK
[921]	572	DO ji = fs_2, fs_jpim1
	573	zmask = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj)
	574	zsang = SIGN(1.0,fcor(ji,jj))*sangvg
	575	z0 = zmass1(ji,jj)/dtevp
	576	! SB modif because ocean has no slip boundary condition
	577	! GG Bug
	578	! zv_ice1 = 0.5( (v_ice(ji,jj)+v_ice(ji,jj-1))e1t(ji+1,jj) &
	579	! & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
	580	! & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
	581	zv_ice1 = 0.5( (v_ice(ji,jj)+v_ice(ji,jj-1))e1t(ji,jj) &
	582	& +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) &
	583	& /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
[825]	584
[921]	585	za = rhocoSQRT((u_ice(ji,jj)-u_oce1(ji,jj))*2 + &
	586	(zv_ice1-v_oce1(ji,jj))*2)(1.0-zfrld1(ji,jj))
	587	zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
	588	za(cangvgu_oce1(ji,jj)-zsang*v_oce1(ji,jj))
	589	zcca = z0+za*cangvg
	590	zccb = zcorl1(ji,jj)+za*zsang
	591	u_ice(ji,jj) = (zr+zccbzv_ice1)/(zcca+epsd)zmask
	592	END DO ! ji
	593	END DO ! jj
[825]	594
[921]	595	CALL lbc_lnk( u_ice(:,:), 'U', -1. )
[2804]	596	#if defined key_agrif
	597	CALL agrif_dyn_lim(jter ,nevp , 'U')
	598	#endif
[921]	599	ENDIF
[825]	600
[921]	601	IF(ln_ctl) THEN
	602	!--- Convergence test.
	603	DO jj = k_j1+1 , k_jpj-1
	604	zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , &
	605	ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
	606	END DO
	607	zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) )
	608	IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain
	609	ENDIF
	610
	611	! ! ==================== !
[868]	612	END DO ! end loop over jter !
[825]	613	! ! ==================== !
	614
[921]	615	!
	616	!------------------------------------------------------------------------------!
	617	! 4) Prevent ice velocities when the ice is thin
	618	!------------------------------------------------------------------------------!
	619	!
[825]	620	! If the ice thickness is below 1cm then ice velocity should equal the
[834]	621	! ocean velocity,
	622	! This prevents high velocity when ice is thin
[868]	623	!CDIR NOVERRCHK
[825]	624	DO jj = k_j1+1, k_jpj-1
[868]	625	!CDIR NOVERRCHK
	626	DO ji = fs_2, fs_jpim1
[825]	627	zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) )
	628	zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
	629	IF ( zdummy .LE. 5.0e-2 ) THEN
[888]	630	u_ice(ji,jj) = u_oce(ji,jj)
	631	v_ice(ji,jj) = v_oce(ji,jj)
[825]	632	ENDIF ! zdummy
	633	END DO
	634	END DO
[866]	635
[869]	636	CALL lbc_lnk( u_ice(:,:), 'U', -1. )
	637	CALL lbc_lnk( v_ice(:,:), 'V', -1. )
[2804]	638	#if defined key_agrif_lim2
	639	CALL agrif_dyn_lim(nevp ,nevp , 'U')
	640	CALL agrif_dyn_lim(nevp ,nevp , 'V')
	641	#endif
[869]	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: