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

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

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

Download in other formats: