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

Context Navigation

source: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 6853

Last change on this file since 6853 was 6853, checked in by clem, 8 years ago
correct bugs on LIM3 rheology and outputs. Make sure the ice thickness distribution (used in initialization and bdy boundary conditions) is a gaussian.
Property svn:keywords set to `Id`
File size: 39.4 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
[6746]	11	!! 3.5 ! 2012-08 (R. Benshila) AGRIF
[6796]	12	!! 3.6 ! 2016-06 (C. Rousset) Rewriting + landfast ice + possibility to use mEVP (Bouillon 2013)
[1244]	13	!!----------------------------------------------------------------------
[6746]	14	#if defined key_lim3
[825]	15	!!----------------------------------------------------------------------
[6746]	16	!! 'key_lim3' LIM-3 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
[6746]	26	USE ice ! ice variables
	27	USE limitd_me ! ice strength
[3625]	28	USE lbclnk ! Lateral Boundary Condition / MPP link
	29	USE lib_mpp ! MPP library
	30	USE wrk_nemo ! work arrays
	31	USE in_out_manager ! I/O manager
	32	USE prtctl ! Print control
	33	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[6746]	34	#if defined key_agrif
[6584]	35	USE agrif_lim3_interp
	36	#endif
[4161]	37	#if defined key_bdy
	38	USE bdyice_lim
	39	#endif
[6853]	40	!!!clem
	41	!! USE iom
[825]	42
	43	IMPLICIT NONE
	44	PRIVATE
	45
[6746]	46	PUBLIC lim_rhg ! routine called by lim_dyn
[825]	47
[868]	48	!! * Substitutions
	49	# include "vectopt_loop_substitute.h90"
[825]	50	!!----------------------------------------------------------------------
[4161]	51	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]	52	!! $Id$
[2528]	53	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	54	!!----------------------------------------------------------------------
	55	CONTAINS
	56
[6746]	57	SUBROUTINE lim_rhg
[825]	58	!!-------------------------------------------------------------------
[834]	59	!! * SUBROUTINE lim_rhg *
	60	!! EVP-C-grid
[825]	61	!!
[834]	62	!! ** purpose : determines sea ice drift from wind stress, ice-ocean
[825]	63	!! stress and sea-surface slope. Ice-ice interaction is described by
[834]	64	!! a non-linear elasto-viscous-plastic (EVP) law including shear
	65	!! strength and a bulk rheology (Hunke and Dukowicz, 2002).
[825]	66	!!
[834]	67	!! The points in the C-grid look like this, dear reader
[825]	68	!!
[834]	69	!! (ji,jj)
	70	!! \|
	71	!! \|
	72	!! (ji-1,jj) \| (ji,jj)
	73	!! ---------
	74	!! \| \|
	75	!! \| (ji,jj) \|------(ji,jj)
	76	!! \| \|
	77	!! ---------
	78	!! (ji-1,jj-1) (ji,jj-1)
[825]	79	!!
[834]	80	!! ** Inputs : - wind forcing (stress), oceanic currents
	81	!! ice total volume (vt_i) per unit area
	82	!! snow total volume (vt_s) per unit area
[825]	83	!!
[834]	84	!! ** Action : - compute u_ice, v_ice : the components of the
	85	!! sea-ice velocity vector
	86	!! - compute delta_i, shear_i, divu_i, which are inputs
	87	!! of the ice thickness distribution
[825]	88	!!
[834]	89	!! ** Steps : 1) Compute ice snow mass, ice strength
	90	!! 2) Compute wind, oceanic stresses, mass terms and
	91	!! coriolis terms of the momentum equation
	92	!! 3) Solve the momentum equation (iterative procedure)
	93	!! 4) Prevent high velocities if the ice is thin
	94	!! 5) Recompute invariants of the strain rate tensor
	95	!! which are inputs of the ITD, store stress
	96	!! for the next time step
	97	!! 6) Control prints of residual (convergence)
	98	!! and charge ellipse.
	99	!! The user should make sure that the parameters
[5123]	100	!! nn_nevp, elastic time scale and rn_creepl maintain stress state
[834]	101	!! on the charge ellipse for plastic flow
	102	!! e.g. in the Canadian Archipelago
	103	!!
[6746]	104	!! ** Notes : There is the possibility to use mEVP from Bouillon 2013
	105	!! (by uncommenting some lines in part 3 and changing alpha and beta parameters)
	106	!! but this solution appears very unstable (see Kimmritz et al 2016)
	107	!!
[2528]	108	!! References : Hunke and Dukowicz, JPO97
	109	!! Bouillon et al., Ocean Modelling 2009
[6746]	110	!! Bouillon et al., Ocean Modelling 2013
[2528]	111	!!-------------------------------------------------------------------
[6746]	112	INTEGER :: ji, jj ! dummy loop indices
	113	INTEGER :: jter ! local integers
[825]	114	CHARACTER (len=50) :: charout
	115
[6796]	116	REAL(wp) :: zdtevp, z1_dtevp ! time step for subcycling
[6789]	117	REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity
	118	REAL(wp) :: zbeta, zalph1, z1_alph1, zalph2, z1_alph2 ! alpha and beta from Bouillon 2009 and 2013
[6796]	119	REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV ! ice/snow mass
[6789]	120	REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars
[6796]	121	REAL(wp) :: zTauO, zTauB, zTauE, zCor, zvel ! temporary scalars
[825]	122
[6789]	123	REAL(wp) :: zsig1, zsig2 ! internal ice stress
	124	REAL(wp) :: zresm ! Maximal error on ice velocity
	125	REAL(wp) :: zintb, zintn ! dummy argument
	126
	127	REAL(wp), POINTER, DIMENSION(:,:) :: z1_e1t0, z1_e2t0 ! scale factors
	128	REAL(wp), POINTER, DIMENSION(:,:) :: zp_delt ! P/delta at T points
	129	!
[6796]	130	REAL(wp), POINTER, DIMENSION(:,:) :: zaU , zaV ! ice fraction on U/V points
	131	REAL(wp), POINTER, DIMENSION(:,:) :: zmU_t, zmV_t ! ice/snow mass/dt on U/V points
	132	REAL(wp), POINTER, DIMENSION(:,:) :: zmf ! coriolis parameter at T points
	133	REAL(wp), POINTER, DIMENSION(:,:) :: zTauU_ia , ztauV_ia ! ice-atm. stress at U-V points
	134	REAL(wp), POINTER, DIMENSION(:,:) :: zspgU , zspgV ! surface pressure gradient at U/V points
	135	REAL(wp), POINTER, DIMENSION(:,:) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points
	136	REAL(wp), POINTER, DIMENSION(:,:) :: zfU , zfV ! internal stresses
[3294]	137
[6789]	138	REAL(wp), POINTER, DIMENSION(:,:) :: zds ! shear
	139	REAL(wp), POINTER, DIMENSION(:,:) :: zs1, zs2, zs12 ! stress tensor components
	140	REAL(wp), POINTER, DIMENSION(:,:) :: zu_ice, zv_ice, zresr ! check convergence
	141	REAL(wp), POINTER, DIMENSION(:,:) :: zpice ! array used for the calculation of ice surface slope:
	142	! ocean surface (ssh_m) if ice is not embedded
	143	! ice top surface if ice is embedded
[6853]	144	REAL(wp), POINTER, DIMENSION(:,:) :: zswitchU, zswitchV ! dummy arrays
	145	REAL(wp), POINTER, DIMENSION(:,:) :: zmaskU, zmaskV ! mask for ice presence
[6789]	146	REAL(wp), POINTER, DIMENSION(:,:) :: zfmask, zwf ! mask at F points for the ice
[5123]	147
[6789]	148	REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter
[6853]	149	REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity equals ocean velocity
[6789]	150	REAL(wp), PARAMETER :: zshlat = 2._wp ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip)
[2528]	151	!!-------------------------------------------------------------------
[3294]	152
[6796]	153	CALL wrk_alloc( jpi,jpj, z1_e1t0, z1_e2t0, zp_delt )
	154	CALL wrk_alloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
	155	CALL wrk_alloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
[6763]	156	CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
[6853]	157	CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
[3294]	158
[6746]	159	#if defined key_agrif
[6584]	160	CALL agrif_interp_lim3('U') ! First interpolation of coarse values
[6746]	161	CALL agrif_interp_lim3('V')
[6584]	162	#endif
[921]	163	!
	164	!------------------------------------------------------------------------------!
[6789]	165	! 0) mask at F points for the ice
[921]	166	!------------------------------------------------------------------------------!
[6789]	167	! ocean/land mask
	168	DO jj = 1, jpjm1
	169	DO ji = 1, jpim1 ! NO vector opt.
	170	zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
	171	END DO
	172	END DO
	173	CALL lbc_lnk( zfmask, 'F', 1._wp )
	174
	175	! Lateral boundary conditions on velocity (modify zfmask)
	176	zwf(:,:) = zfmask(:,:)
	177	DO jj = 2, jpjm1
	178	DO ji = fs_2, fs_jpim1 ! vector opt.
	179	IF( zfmask(ji,jj) == 0._wp ) THEN
	180	zfmask(ji,jj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
	181	ENDIF
	182	END DO
	183	END DO
	184	DO jj = 2, jpjm1
	185	IF( zfmask(1,jj) == 0._wp ) THEN
	186	zfmask(1 ,jj) = zshlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
	187	ENDIF
	188	IF( zfmask(jpi,jj) == 0._wp ) THEN
	189	zfmask(jpi,jj) = zshlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
	190	ENDIF
	191	END DO
	192	DO ji = 2, jpim1
	193	IF( zfmask(ji,1) == 0._wp ) THEN
	194	zfmask(ji,1 ) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
	195	ENDIF
	196	IF( zfmask(ji,jpj) == 0._wp ) THEN
	197	zfmask(ji,jpj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
	198	ENDIF
	199	END DO
	200	CALL lbc_lnk( zfmask, 'F', 1._wp )
	201
	202	!------------------------------------------------------------------------------!
	203	! 1) define some variables and initialize arrays
	204	!------------------------------------------------------------------------------!
[6746]	205	! ecc2: square of yield ellipse eccenticrity
	206	ecc2 = rn_ecc * rn_ecc
	207	z1_ecc2 = 1._wp / ecc2
[825]	208
[6746]	209	! Time step for subcycling
[6796]	210	zdtevp = rdt_ice / REAL( nn_nevp )
	211	z1_dtevp = 1._wp / zdtevp
[825]	212
[6746]	213	! alpha parameters (Bouillon 2009)
	214	zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
	215	zalph2 = zalph1 * z1_ecc2
[825]	216
[6746]	217	! alpha and beta parameters (Bouillon 2013)
	218	!!zalph1 = 40.
	219	!!zalph2 = 40.
	220	!!zbeta = 3000.
	221	!!zbeta = REAL( nn_nevp ) ! close to classical EVP of Hunke (2001)
	222
	223	z1_alph1 = 1._wp / ( zalph1 + 1._wp )
	224	z1_alph2 = 1._wp / ( zalph2 + 1._wp )
	225
	226	! Initialise stress tensor
	227	zs1 (:,:) = stress1_i (:,:)
	228	zs2 (:,:) = stress2_i (:,:)
	229	zs12(:,:) = stress12_i(:,:)
	230
	231	! Ice strength
	232	CALL lim_itd_me_icestrength( nn_icestr )
	233
	234	! scale factors
	235	DO jj = 2, jpjm1
	236	DO ji = fs_2, fs_jpim1
	237	z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj ) + e1t(ji,jj ) )
	238	z1_e2t0(ji,jj) = 1._wp / ( e2t(ji ,jj+1) + e2t(ji,jj ) )
[825]	239	END DO
	240	END DO
[6746]	241
[921]	242	!
	243	!------------------------------------------------------------------------------!
	244	! 2) Wind / ocean stress, mass terms, coriolis terms
	245	!------------------------------------------------------------------------------!
	246
[3625]	247	IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==!
[5123]	248	!
	249	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
	250	! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
[3625]	251	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
[5123]	252	!
	253	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
	254	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
[3625]	255	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
[5123]	256	!
[6746]	257	zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
[5123]	258	!
[3625]	259	ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==!
	260	zpice(:,:) = ssh_m(:,:)
	261	ENDIF
	262
[6746]	263	DO jj = 2, jpjm1
[868]	264	DO ji = fs_2, fs_jpim1
[825]	265
[6746]	266	! ice fraction at U-V points
[6801]	267	zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e12t(ji,jj) + at_i(ji+1,jj) * e12t(ji+1,jj) ) * r1_e12u(ji,jj) * umask(ji,jj,1)
	268	zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e12t(ji,jj) + at_i(ji,jj+1) * e12t(ji,jj+1) ) * r1_e12v(ji,jj) * vmask(ji,jj,1)
[4990]	269
[6746]	270	! Ice/snow mass at U-V points
	271	zm1 = ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) )
	272	zm2 = ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) )
	273	zm3 = ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) )
[6801]	274	zmassU = 0.5_wp * ( zm1 * e12t(ji,jj) + zm2 * e12t(ji+1,jj) ) * r1_e12u(ji,jj) * umask(ji,jj,1)
	275	zmassV = 0.5_wp * ( zm1 * e12t(ji,jj) + zm3 * e12t(ji,jj+1) ) * r1_e12v(ji,jj) * vmask(ji,jj,1)
[825]	276
[6746]	277	! Ocean currents at U-V points
[6796]	278	v_oceU(ji,jj) = 0.5_wp * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji+1,jj) &
	279	& + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
[6746]	280
[6796]	281	u_oceV(ji,jj) = 0.5_wp * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj+1) &
	282	& + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
[825]	283
[6796]	284	! Coriolis at T points (m*f)
	285	zmf(ji,jj) = zm1 * ff_t(ji,jj)
[825]	286
[6796]	287	! m/dt
	288	zmU_t(ji,jj) = zmassU * z1_dtevp
	289	zmV_t(ji,jj) = zmassV * z1_dtevp
	290
	291	! Drag ice-atm.
	292	zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
	293	zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
	294
[6746]	295	! Surface pressure gradient (- mgGRAD(ssh)) at U-V points
[6796]	296	zspgU(ji,jj) = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
	297	zspgV(ji,jj) = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
[825]	298
[6853]	299	! masks
	300	zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice
	301	zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice
	302
[6746]	303	! switches
[6853]	304	zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
	305	zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
[825]	306
	307	END DO
	308	END DO
[6796]	309	CALL lbc_lnk( zmf, 'T', 1. )
[921]	310	!
	311	!------------------------------------------------------------------------------!
	312	! 3) Solution of the momentum equation, iterative procedure
	313	!------------------------------------------------------------------------------!
	314	!
[2528]	315	! !----------------------!
[5123]	316	DO jter = 1 , nn_nevp ! loop over jter !
[2528]	317	! !----------------------!
[6763]	318	IF(ln_ctl) THEN ! Convergence test
[6746]	319	DO jj = 1, jpjm1
	320	zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
	321	zv_ice(:,jj) = v_ice(:,jj)
	322	END DO
	323	ENDIF
[825]	324
[6746]	325	! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
[6763]	326	DO jj = 1, jpjm1 ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
[6796]	327	DO ji = 1, jpim1
[825]	328
[6746]	329	! shear at F points
[5123]	330	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) &
	331	& + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
[6789]	332	& ) * r1_e12f(ji,jj) * zfmask(ji,jj)
[825]	333
[921]	334	END DO
	335	END DO
[6746]	336	CALL lbc_lnk( zds, 'F', 1. )
[921]	337
[6746]	338	DO jj = 2, jpjm1
	339	DO ji = 2, jpim1 ! no vector loop
[825]	340
[6746]	341	! shear**2 at T points (doc eq. A16)
	342	zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e12f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e12f(ji-1,jj ) &
	343	& + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1) &
	344	& ) * 0.25_wp * r1_e12t(ji,jj)
	345
[6763]	346	! divergence at T points
	347	zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
	348	& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
	349	& ) * r1_e12t(ji,jj)
	350	zdiv2 = zdiv * zdiv
	351
	352	! tension at T points
	353	zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) &
	354	& - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
	355	& ) * r1_e12t(ji,jj)
	356	zdt2 = zdt * zdt
	357
[6746]	358	! delta at T points
[6763]	359	zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * usecc2 )
[4205]	360
[6746]	361	! P/delta at T points
[6796]	362	zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
[6763]	363
	364	! stress at T points
	365	zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
	366	zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
	367
[6746]	368	END DO
	369	END DO
	370	CALL lbc_lnk( zp_delt, 'T', 1. )
[825]	371
[6763]	372	DO jj = 1, jpjm1
	373	DO ji = 1, jpim1
	374
[6746]	375	! P/delta at F points
	376	zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )
	377
[6763]	378	! stress at F points
	379	zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
	380
[5123]	381	END DO
	382	END DO
[6746]	383	CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
	384
[6763]	385
[6746]	386	! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
	387	DO jj = 2, jpjm1
	388	DO ji = fs_2, fs_jpim1
[6763]	389
[6746]	390	! U points
[6796]	391	zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) &
[6746]	392	& + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) &
	393	& ) * r1_e2u(ji,jj) &
	394	& + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) &
	395	& ) * 2._wp * r1_e1u(ji,jj) &
	396	& ) * r1_e12u(ji,jj)
[825]	397
[6746]	398	! V points
[6796]	399	zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) &
[6746]	400	& - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) &
	401	& ) * r1_e1v(ji,jj) &
	402	& + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) &
	403	& ) * 2._wp * r1_e2v(ji,jj) &
	404	& ) * r1_e12v(ji,jj)
[6763]	405
	406	! u_ice at V point
[6796]	407	u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) &
[6763]	408	& + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
	409
	410	! v_ice at U point
[6796]	411	v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) &
[6763]	412	& + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
	413
[921]	414	END DO
	415	END DO
[825]	416	!
[6746]	417	! --- Computation of ice velocity --- !
[6796]	418	! Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp
[6746]	419	! Bouillon et al. 2009 (eq 34-35) => stable
	420	IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations
	421
	422	DO jj = 2, jpjm1
[921]	423	DO ji = fs_2, fs_jpim1
[825]	424
[6796]	425	! tau_io/(v_oce - v_ice)
	426	zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) &
	427	& + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
	428
	429	! tau_bottom/v_ice
	430	zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
[6801]	431	zTauB = - tau_icebfr(ji,jj) / zvel
[6796]	432
	433	! Coriolis at V-points (energy conserving formulation)
	434	zCor = - 0.25_wp * r1_e2v(ji,jj) * &
	435	& ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) &
	436	& + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
	437
	438	! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
	439	zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
	440
	441	! landfast switch => 0 = static friction ; 1 = sliding friction
[6801]	442	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
[6763]	443
[6796]	444	! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
[6853]	445	v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity
	446	& + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	447	& ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	448	& + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
	449	& ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
	450	& ) * zmaskV(ji,jj)
[6746]	451	! Bouillon 2013
[6796]	452	!!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) ) &
	453	!! & + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj) &
[6853]	454	!! & ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
[6746]	455
[921]	456	END DO
	457	END DO
[6746]	458	CALL lbc_lnk( v_ice, 'V', -1. )
	459
	460	#if defined key_agrif
	461	CALL agrif_interp_lim3( 'V' )
[3680]	462	#endif
[4333]	463	#if defined key_bdy
[6746]	464	CALL bdy_ice_lim_dyn( 'V' )
[4333]	465	#endif
[825]	466
[6746]	467	DO jj = 2, jpjm1
[921]	468	DO ji = fs_2, fs_jpim1
[6796]	469
	470	! tau_io/(u_oce - u_ice)
	471	zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) &
	472	& + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
[834]	473
[6796]	474	! tau_bottom/u_ice
	475	zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
[6801]	476	zTauB = - tau_icebfr(ji,jj) / zvel
[6796]	477
	478	! Coriolis at U-points (energy conserving formulation)
	479	zCor = 0.25_wp * r1_e1u(ji,jj) * &
	480	& ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) &
	481	& + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
[6746]	482
[6796]	483	! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
	484	zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
	485
	486	! landfast switch => 0 = static friction ; 1 = sliding friction
[6801]	487	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
[6796]	488
	489	! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
[6853]	490	u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity
	491	& + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	492	& ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	493	& + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
	494	& ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
	495	& ) * zmaskU(ji,jj)
[6746]	496	! Bouillon 2013
[6796]	497	!!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) ) &
	498	!! & + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj) &
[6853]	499	!! & ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
[921]	500	END DO
	501	END DO
[6746]	502	CALL lbc_lnk( u_ice, 'U', -1. )
	503
	504	#if defined key_agrif
	505	CALL agrif_interp_lim3( 'U' )
[3680]	506	#endif
[4333]	507	#if defined key_bdy
[6746]	508	CALL bdy_ice_lim_dyn( 'U' )
[4333]	509	#endif
[825]	510
[6746]	511	ELSE ! odd iterations
	512
	513	DO jj = 2, jpjm1
[921]	514	DO ji = fs_2, fs_jpim1
[825]	515
[6796]	516	! tau_io/(u_oce - u_ice)
	517	zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) &
	518	& + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
[6746]	519
[6796]	520	! tau_bottom/u_ice
	521	zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
[6801]	522	zTauB = - tau_icebfr(ji,jj) / zvel
[6746]	523
[6796]	524	! Coriolis at U-points (energy conserving formulation)
	525	zCor = 0.25_wp * r1_e1u(ji,jj) * &
	526	& ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) &
	527	& + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
	528
	529	! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
	530	zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
	531
	532	! landfast switch => 0 = static friction ; 1 = sliding friction
[6801]	533	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
[6796]	534
	535	! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
[6853]	536	u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity
	537	& + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	538	& ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	539	& + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
	540	& ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
	541	& ) * zmaskU(ji,jj)
[6746]	542	! Bouillon 2013
[6796]	543	!!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) ) &
	544	!! & + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj) &
[6853]	545	!! & ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
[921]	546	END DO
	547	END DO
[6746]	548	CALL lbc_lnk( u_ice, 'U', -1. )
	549
	550	#if defined key_agrif
	551	CALL agrif_interp_lim3( 'U' )
[3680]	552	#endif
[4333]	553	#if defined key_bdy
[6746]	554	CALL bdy_ice_lim_dyn( 'U' )
[4333]	555	#endif
[825]	556
[6746]	557	DO jj = 2, jpjm1
[921]	558	DO ji = fs_2, fs_jpim1
[6796]	559
	560	! tau_io/(v_oce - v_ice)
	561	zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) &
	562	& + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
[825]	563
[6796]	564	! tau_bottom/v_ice
	565	zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
	566	ztauB = - tau_icebfr(ji,jj) / zvel
	567
	568	! Coriolis at V-points (energy conserving formulation)
	569	zCor = - 0.25_wp * r1_e2v(ji,jj) * &
	570	& ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) &
	571	& + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
[6746]	572
[6796]	573	! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
	574	zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
	575
[6853]	576	! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction
	577	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
[6746]	578
[6796]	579	! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
[6853]	580	v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity
	581	& + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	582	& ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
	583	& + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
	584	& ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
	585	& ) * zmaskV(ji,jj)
[6746]	586	! Bouillon 2013
[6796]	587	!!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) ) &
	588	!! & + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj) &
[6853]	589	!! & ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
[5123]	590	END DO
	591	END DO
[6746]	592	CALL lbc_lnk( v_ice, 'V', -1. )
	593
	594	#if defined key_agrif
	595	CALL agrif_interp_lim3( 'V' )
[3680]	596	#endif
[4333]	597	#if defined key_bdy
[6746]	598	CALL bdy_ice_lim_dyn( 'V' )
[4333]	599	#endif
[825]	600
[921]	601	ENDIF
[4161]	602
[6763]	603	IF(ln_ctl) THEN ! Convergence test
[6746]	604	DO jj = 2 , jpjm1
[5123]	605	zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
[921]	606	END DO
[6746]	607	zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
[921]	608	IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain
	609	ENDIF
[6746]	610	!
[4161]	611	! ! ==================== !
[868]	612	END DO ! end loop over jter !
[825]	613	! ! ==================== !
[921]	614	!
	615	!------------------------------------------------------------------------------!
[6853]	616	! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
[921]	617	!------------------------------------------------------------------------------!
[6763]	618	DO jj = 1, jpjm1
[6796]	619	DO ji = 1, jpim1
[6746]	620
	621	! shear at F points
	622	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) &
	623	& + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
[6789]	624	& ) * r1_e12f(ji,jj) * zfmask(ji,jj)
[6763]	625
[868]	626	END DO
[6763]	627	END DO
	628	CALL lbc_lnk( zds, 'F', 1. )
[6746]	629
	630	DO jj = 2, jpjm1
	631	DO ji = 2, jpim1 ! no vector loop
[4161]	632
[6763]	633	! tension**2 at T points
	634	zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) &
	635	& - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
	636	& ) * r1_e12t(ji,jj)
	637	zdt2 = zdt * zdt
	638
[6746]	639	! shear**2 at T points (doc eq. A16)
	640	zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e12f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e12f(ji-1,jj ) &
	641	& + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1) &
	642	& ) * 0.25_wp * r1_e12t(ji,jj)
	643
[6763]	644	! shear at T points
	645	shear_i(ji,jj) = SQRT( zdt2 + zds2 )
	646
	647	! divergence at T points
	648	divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
	649	& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
	650	& ) * r1_e12t(ji,jj)
	651
[6746]	652	! delta at T points
[6763]	653	zdelta = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * usecc2 )
[6853]	654	rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
	655	delta_i(ji,jj) = zdelta + rn_creepl * rswitch
[6763]	656
[5123]	657	END DO
	658	END DO
[6763]	659	CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. )
[6746]	660
	661	! --- Store the stress tensor for the next time step --- !
[868]	662	stress1_i (:,:) = zs1 (:,:)
	663	stress2_i (:,:) = zs2 (:,:)
	664	stress12_i(:,:) = zs12(:,:)
[6746]	665	!
[868]	666
[921]	667	!------------------------------------------------------------------------------!
[6853]	668	! 5) Control prints of residual and charge ellipse
[921]	669	!------------------------------------------------------------------------------!
	670	!
[834]	671	! print the residual for convergence
	672	IF(ln_ctl) THEN
[868]	673	WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, jter
[834]	674	CALL prt_ctl_info(charout)
	675	CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
	676	ENDIF
[825]	677
[834]	678	! print charge ellipse
	679	! This can be desactivated once the user is sure that the stress state
	680	! lie on the charge ellipse. See Bouillon et al. 08 for more details
[825]	681	IF(ln_ctl) THEN
	682	CALL prt_ctl_info('lim_rhg : numit :',ivar1=numit)
	683	CALL prt_ctl_info('lim_rhg : nwrite :',ivar1=nwrite)
	684	CALL prt_ctl_info('lim_rhg : MOD :',ivar1=MOD(numit,nwrite))
	685	IF( MOD(numit,nwrite) .EQ. 0 ) THEN
	686	WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. '
	687	CALL prt_ctl_info(charout)
[6746]	688	DO jj = 2, jpjm1
[825]	689	DO ji = 2, jpim1
[6796]	690	IF (strength(ji,jj) > 1.0) THEN
	691	zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)*2 + 4zs12(ji,jj)*2 ) ) / ( 2strength(ji,jj) )
	692	zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)*2 + 4zs12(ji,jj)*2 ) ) / ( 2strength(ji,jj) )
[825]	693	WRITE(charout,FMT="('lim_rhg :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
	694	CALL prt_ctl_info(charout)
	695	ENDIF
	696	END DO
	697	END DO
	698	WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 2000, numit, 0, 0, ' ch. ell. '
	699	CALL prt_ctl_info(charout)
	700	ENDIF
	701	ENDIF
[2715]	702	!
[6853]	703
	704	!! ! velocity
	705	!! IF ( iom_use( "uice_ipa" ) .OR. iom_use( "vice_ipa" ) .OR. iom_use( "icevel" ) ) THEN
	706	!! CALL iom_put( "uice_ipa" , u_ice ) ! ice velocity u component
	707	!! CALL iom_put( "vice_ipa" , v_ice ) ! ice velocity v component
	708	!! ENDIF
	709	!! IF ( iom_use( "iceconc_cat" ) ) CALL iom_put( "iceconc_cat" , a_i ) ! area for categories
	710	!! IF ( iom_use( "icethic_cat" ) ) CALL iom_put( "icethic_cat" , ht_i ) ! thickness for categories
	711
	712
[6796]	713	CALL wrk_dealloc( jpi,jpj, z1_e1t0, z1_e2t0, zp_delt )
	714	CALL wrk_dealloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
	715	CALL wrk_dealloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
[6763]	716	CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
[6853]	717	CALL wrk_dealloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
[3294]	718
[825]	719	END SUBROUTINE lim_rhg
	720
	721	#else
	722	!!----------------------------------------------------------------------
	723	!! Default option Dummy module NO LIM sea-ice model
	724	!!----------------------------------------------------------------------
	725	CONTAINS
[6746]	726	SUBROUTINE lim_rhg ! Dummy routine
	727	WRITE(,) 'lim_rhg: You should not have seen this print! error?'
[825]	728	END SUBROUTINE lim_rhg
	729	#endif
	730
	731	!!==============================================================================
	732	END MODULE limrhg

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

Download in other formats: