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

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 7738

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

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

Download in other formats: