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 @ 5122

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

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

Download in other formats: