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

Context Navigation

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 8373

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

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

Download in other formats: