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

Context Navigation

source: trunk/NEMO/LIM_SRC_3/limrhg.F90 @ 834

Last change on this file since 834 was 834, checked in by ctlod, 16 years ago
Clean comments and useless lines, see ticket:#72
File size: 37.1 KB

Rev	Line
[825]	1	MODULE limrhg
	2	!!======================================================================
	3	!! * MODULE limrhg *
[834]	4	!! Ice rheology : sea ice rheology
[825]	5	!!======================================================================
	6	#if defined key_lim3
	7	!!----------------------------------------------------------------------
[834]	8	!! 'key_lim3' LIM3 sea-ice model
[825]	9	!!----------------------------------------------------------------------
	10	!! lim_rhg : computes ice velocities
	11	!!----------------------------------------------------------------------
	12	!! * Modules used
	13	USE phycst
	14	USE par_oce
	15	USE ice_oce ! ice variables
	16	USE dom_oce
	17	USE dom_ice
	18	USE ice
	19	USE iceini
	20	USE lbclnk
	21	USE lib_mpp
	22	USE in_out_manager ! I/O manager
	23	USE limitd_me
	24	USE prtctl ! Print control
	25
	26
	27	IMPLICIT NONE
	28	PRIVATE
	29
	30	!! * Routine accessibility
	31	PUBLIC lim_rhg ! routine called by lim_dyn
	32
	33	!! * Module variables
	34	REAL(wp) :: & ! constant values
	35	rzero = 0.e0 , &
	36	rone = 1.e0
	37	!!----------------------------------------------------------------------
[834]	38	!! LIM 3.0, UCL-LOCEAN-IPSL (2008)
[825]	39	!! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limrhg.F90,v 1.5 2005/03/27 18:34:42 opalod Exp $
	40	!! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
	41	!!----------------------------------------------------------------------
	42
	43	CONTAINS
	44
	45	SUBROUTINE lim_rhg( k_j1, k_jpj )
[834]	46
[825]	47	!!-------------------------------------------------------------------
[834]	48	!! * SUBROUTINE lim_rhg *
	49	!! EVP-C-grid
[825]	50	!!
[834]	51	!! ** purpose : determines sea ice drift from wind stress, ice-ocean
[825]	52	!! stress and sea-surface slope. Ice-ice interaction is described by
[834]	53	!! a non-linear elasto-viscous-plastic (EVP) law including shear
	54	!! strength and a bulk rheology (Hunke and Dukowicz, 2002).
[825]	55	!!
[834]	56	!! The points in the C-grid look like this, dear reader
[825]	57	!!
[834]	58	!! (ji,jj)
	59	!! \|
	60	!! \|
	61	!! (ji-1,jj) \| (ji,jj)
	62	!! ---------
	63	!! \| \|
	64	!! \| (ji,jj) \|------(ji,jj)
	65	!! \| \|
	66	!! ---------
	67	!! (ji-1,jj-1) (ji,jj-1)
[825]	68	!!
[834]	69	!! ** Inputs : - wind forcing (stress), oceanic currents
	70	!! ice total volume (vt_i) per unit area
	71	!! snow total volume (vt_s) per unit area
[825]	72	!!
[834]	73	!! ** Action : - compute u_ice, v_ice : the components of the
	74	!! sea-ice velocity vector
	75	!! - compute delta_i, shear_i, divu_i, which are inputs
	76	!! of the ice thickness distribution
[825]	77	!!
[834]	78	!! ** Steps : 1) Compute ice snow mass, ice strength
	79	!! 2) Compute wind, oceanic stresses, mass terms and
	80	!! coriolis terms of the momentum equation
	81	!! 3) Solve the momentum equation (iterative procedure)
	82	!! 4) Prevent high velocities if the ice is thin
	83	!! 5) Recompute invariants of the strain rate tensor
	84	!! which are inputs of the ITD, store stress
	85	!! for the next time step
	86	!! 6) Control prints of residual (convergence)
	87	!! and charge ellipse.
	88	!! The user should make sure that the parameters
	89	!! nevp, telast and creepl maintain stress state
	90	!! on the charge ellipse for plastic flow
	91	!! e.g. in the Canadian Archipelago
	92	!!
	93	!! ** References : Hunke and Dukowicz, JPO97
	94	!! Bouillon et al., 08, in prep (update this when
	95	!! published)
	96	!! Vancoppenolle et al., OM08
	97	!!
[825]	98	!! History :
[834]	99	!! 1.0 ! 07-03 (M.A. Morales Maqueda, S. Bouillon)
	100	!! 2.0 ! 08-03 M. Vancoppenolle : LIM3
[825]	101	!!
	102	!!-------------------------------------------------------------------
	103	! * Arguments
	104	!
	105	INTEGER, INTENT(in) :: &
[834]	106	k_j1 , & !: southern j-index for ice computation
	107	k_jpj !: northern j-index for ice computation
[825]	108
	109	! * Local variables
[834]	110	INTEGER :: ji, jj !: dummy loop indices
[825]	111
	112	INTEGER :: &
[834]	113	iter, jter !: temporary integers
[825]	114
	115	CHARACTER (len=50) :: charout
	116
	117	REAL(wp) :: &
[834]	118	zt11, zt12, zt21, zt22, & !: temporary scalars
	119	ztagnx, ztagny, & !: wind stress on U/V points
	120	delta !
[825]	121
	122	REAL(wp) :: &
[834]	123	za, & !:
	124	zstms, & !: temporary scalar for ice strength
	125	zsang, & !: temporary scalar for coriolis term
	126	zmask !: mask for the computation of ice mass
[825]	127
[834]	128	REAL(wp),DIMENSION(jpi,jpj) :: &
	129	zpresh , & !: temporary array for ice strength
	130	zpreshc , & !: Ice strength on grid cell corners (zpreshc)
	131	zfrld1, zfrld2, & !: lead fraction on U/V points
	132	zmass1, zmass2, & !: ice/snow mass on U/V points
	133	zcorl1, zcorl2, & !: coriolis parameter on U/V points
	134	za1ct, za2ct , & !: temporary arrays
	135	zc1 , & !: ice mass
	136	zusw , & !: temporary weight for the computation
	137	!: of ice strength
	138	u_oce1, v_oce1, & !: ocean u/v component on U points
	139	u_oce2, v_oce2, & !: ocean u/v component on V points
	140	u_ice2, & !: ice u component on V point
	141	v_ice1 !: ice v component on U point
[825]	142
	143	REAL(wp) :: &
[834]	144	dtevp, & !: time step for subcycling
	145	dtotel, & !:
	146	ecc2, & !: square of yield ellipse eccenticity
	147	z0, & !: temporary scalar
	148	zr, & !: temporary scalar
	149	zcca, zccb, & !: temporary scalars
	150	zu_ice2, & !:
	151	zv_ice1, & !:
	152	zddc, zdtc, & !: temporary array for delta on corners
	153	zdst, & !: temporary array for delta on centre
	154	zdsshx, zdsshy, & !: term for the gradient of ocean surface
	155	sigma1, sigma2 !: internal ice stress
[825]	156
[834]	157	REAL(wp),DIMENSION(jpi,jpj) :: &
	158	zf1, zf2 !: arrays for internal stresses
[825]	159
[834]	160	REAL(wp),DIMENSION(jpi,jpj) :: &
	161	zdd, zdt, & !: Divergence and tension at centre of grid cells
	162	zds, & !: Shear on northeast corner of grid cells
	163	deltat, & !: Delta at centre of grid cells
	164	deltac, & !: Delta on corners
	165	zs1, zs2, & !: Diagonal stress tensor components zs1 and zs2
	166	zs12 !: Non-diagonal stress tensor component zs12
[825]	167
[834]	168	REAL(wp) :: &
	169	zresm , & !: Maximal error on ice velocity
	170	zindb , & !: ice (1) or not (0)
	171	zdummy !: dummy argument
[825]	172
[834]	173	REAL(wp),DIMENSION(jpi,jpj) :: &
	174	zu_ice , & !: Ice velocity on previous time step
	175	zv_ice , &
	176	zresr !: Local error on velocity
[825]	177
	178	INTEGER :: &
	179	stress_mean_swi = 1
	180	!
	181	!------------------------------------------------------------------------------!
[834]	182	! 1) Ice-Snow mass (zc1), ice strength (zpresh) !
[825]	183	!------------------------------------------------------------------------------!
	184	!
	185	! Put every vector to 0
	186	zpresh(:,:) = 0.0 ; zpreshc(:,:) = 0.0 ; zfrld1(:,:) = 0.0
	187	zmass1(:,:) = 0.0 ; zcorl1(:,:) = 0.0 ; zcorl2(:,:) = 0.0
	188	za1ct(:,:) = 0.0 ; za2ct(:,:) = 0.0
	189	zc1(:,:) = 0.0 ; zusw(:,:) = 0.0
	190	u_oce1(:,:) = 0.0 ; v_oce1(:,:) = 0.0 ; u_oce2(:,:) = 0.0
	191	v_oce2(:,:) = 0.0 ; u_ice2(:,:) = 0.0 ; v_ice1(:,:) = 0.0
[834]	192	zf1(:,:) = 0.0 ; zf2(:,:) = 0.0
[825]	193	zdd(:,:) = 0.0 ; zdt(:,:) = 0.0 ; zds(:,:) = 0.0
	194	deltat(:,:) = 0.0 ; deltac(:,:) = 0.0
[834]	195	zpresh(:,:) = 0.0
[825]	196
[834]	197	! Ice strength on T-points
[825]	198	CALL lim_itd_me_icestrength(ridge_scheme_swi)
	199
[834]	200	! Ice mass and temp variables
[825]	201	DO jj = k_j1 , k_jpj-1
	202	DO ji = 1 , jpi
	203	zc1(ji,jj) = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) )
	204	zpresh(ji,jj) = tms(ji,jj) * strength(ji,jj) / 2.
	205	! tmi = 1 where ther is ice or on land
[834]	206	! Claude sees a bug, next line : tmi(ji,jj)
[825]	207	tmi = 1.0 - ( 1.0 - MAX( 0.0 , SIGN ( 1.0 , vt_i(ji,jj) - &
	208	epsd ) ) ) * tms(ji,jj)
	209	END DO
	210	END DO
	211
	212	CALL lbc_lnk( zc1(:,:), 'T', 1. )
	213	CALL lbc_lnk( zpresh(:,:), 'T', 1. )
	214
[834]	215	! Claude sees a bug
	216	! CALL lbc_lnk( tmi(:,:), 'T', 1. )
	217
	218	! Ice strength on grid cell corners (zpreshc)
	219	! needed for calculation of shear stress
[825]	220	DO jj = k_j1+1, k_jpj-1
	221	DO ji = 2, jpim1
[834]	222	zstms = tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
	223	& tms(ji,jj+1) * wght(ji+1,jj+1,1,2) + &
	224	& tms(ji+1,jj) * wght(ji+1,jj+1,2,1) + &
	225	& tms(ji,jj) * wght(ji+1,jj+1,1,1)
[825]	226	zusw(ji,jj) = 1.0 / MAX( zstms, epsd )
[834]	227	zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
	228	& zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + &
	229	& zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + &
	230	& zpresh(ji,jj) * wght(ji+1,jj+1,1,1) &
	231	& ) * zusw(ji,jj)
[825]	232	END DO
	233	END DO
	234
	235	CALL lbc_lnk( zpreshc(:,:), 'F', 1. )
	236	!
	237	!------------------------------------------------------------------------------!
	238	! 2) Wind / ocean stress, mass terms, coriolis terms
	239	!------------------------------------------------------------------------------!
	240	!
	241	! Wind stress, coriolis and mass terms on the sides of the squares
	242	! zfrld1: lead fraction on U-points
	243	! zfrld2: lead fraction on V-points
	244	! zmass1: ice/snow mass on U-points
	245	! zmass2: ice/snow mass on V-points
	246	! zcorl1: Coriolis parameter on U-points
	247	! zcorl2: Coriolis parameter on V-points
	248	! (ztagnx,ztagny): wind stress on U/V points
	249	! u_oce1: ocean u component on u points
	250	! v_oce1: ocean v component on u points
	251	! u_oce2: ocean u component on v points
	252	! v_oce2: ocean v component on v points
	253
	254	DO jj = k_j1+1, k_jpj-1
	255	DO ji = 2, jpim1
	256
	257	zt11 = tms(ji,jj)*e1t(ji,jj)
	258	zt12 = tms(ji+1,jj)*e1t(ji+1,jj)
	259	zt21 = tms(ji,jj)*e2t(ji,jj)
	260	zt22 = tms(ji,jj+1)*e2t(ji,jj+1)
	261
	262	! Leads area.
	263	zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + &
	264	& zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd )
	265	zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + &
	266	& zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd )
	267
	268	! Mass, coriolis coeff. and currents
[834]	269	zmass1(ji,jj) = ( zt12zc1(ji,jj) + zt11zc1(ji+1,jj) ) / (zt11+zt12+epsd)
	270	zmass2(ji,jj) = ( zt22zc1(ji,jj) + zt21zc1(ji,jj+1) ) / (zt21+zt22+epsd)
	271	zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + &
	272	e1t(ji,jj)*fcor(ji+1,jj) ) &
	273	/ (e1t(ji,jj) + e1t(ji+1,jj) + epsd )
	274	zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)*fcor(ji,jj) + &
	275	e2t(ji,jj)*fcor(ji,jj+1) ) &
	276	/ ( e2t(ji,jj+1) + e2t(ji,jj) + epsd )
[825]	277	!
	278	u_oce1(ji,jj) = u_io(ji,jj)
[834]	279	v_oce2(ji,jj) = v_io(ji,jj)
[825]	280
[834]	281	! Ocean has no slip boundary condition
[825]	282	v_oce1(ji,jj) = 0.5( (v_io(ji,jj)+v_io(ji,jj-1))e1t(ji+1,jj) &
	283	& +(v_io(ji+1,jj)+v_io(ji+1,jj-1))*e1t(ji,jj)) &
	284	& /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
	285
[834]	286	u_oce2(ji,jj) = 0.5((u_io(ji,jj)+u_io(ji-1,jj))e2t(ji,jj+1) &
[825]	287	& +(u_io(ji,jj+1)+u_io(ji-1,jj+1))*e2t(ji,jj)) &
[834]	288	& / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
[825]	289
	290	! Wind stress.
	291	ztagnx = ( 1. - zfrld1(ji,jj) ) * gtaux(ji,jj)
	292	ztagny = ( 1. - zfrld2(ji,jj) ) * gtauy(ji,jj)
	293
[834]	294	! Computation of the velocity field taking into account the ice internal interaction.
[825]	295	! Terms that are independent of the velocity field.
	296
	297	! SB On utilise maintenant le gradient de la pente de l'ocean
	298	! include it later
	299	! zdsshx = (ssh_io(ji+1,jj) - ssh_io(ji,jj))/e1u(ji,jj)
	300	! zdsshy = (ssh_io(ji,jj+1) - ssh_io(ji,jj))/e2v(ji,jj)
[834]	301
[825]	302	zdsshx = 0.0
	303	zdsshy = 0.0
	304
	305	za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx
	306	za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy
	307
	308	END DO
	309	END DO
	310
	311	!
	312	!------------------------------------------------------------------------------!
[834]	313	! 3) Solution of the momentum equation, iterative procedure
[825]	314	!------------------------------------------------------------------------------!
	315	!
	316	! Time step for subcycling
	317	dtevp = rdt_ice / nevp
	318	dtotel = dtevp / ( 2.0 * telast )
	319
	320	!-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)
	321	ecc2 = ecc*ecc
	322
	323	!-Initialise stress tensor
	324	IF ( stress_mean_swi .NE. 1 ) THEN
	325	zs1(:,:) = 0.0
	326	zs2(:,:) = 0.0
	327	zs12(:,:) = 0.0
	328	ELSE
	329	zs1(:,:) = stress1_i(:,:)
	330	zs2(:,:) = stress2_i(:,:)
	331	zs12(:,:) = stress12_i(:,:)
	332	ENDIF
	333
	334	v_ice1(:,:) = 0.0
	335	u_ice2(:,:) = 0.0
	336
	337	zdd(:,:) = 0.0
	338	zdt(:,:) = 0.0
	339	zds(:,:) = 0.0
	340	deltat(:,:) = 0.0
	341	!----------------------!
	342	DO iter = 1 , nevp ! loop over iter !
	343	!----------------------!
	344	DO jj = k_j1, k_jpj-1
	345	zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
	346	zv_ice(:,jj) = v_ice(:,jj)
	347	END DO
	348
[834]	349	DO jj = k_j1+1, k_jpj-1
	350	DO ji = 2, jpim1
[825]	351
	352	!
	353	!- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002)
	354	!- zdd(:,:), zdt(:,:): divergence and tension at centre of grid cells
	355	!- zds(:,:): shear on northeast corner of grid cells
	356	!
[834]	357	!- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded,
	358	! there are many repeated calculations.
	359	! Speed could be improved by regrouping terms. For
[825]	360	! the moment, however, the stress is on clarity of coding to avoid
[834]	361	! bugs (Martin, for Miguel).
	362	!
	363	!- ALSO: arrays zdd, zdt, zds and delta could
	364	! be removed in the future to minimise memory demand.
	365	!
[825]	366	!- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of
	367	! grid cells, exactly as in the B grid case. For simplicity, the indexation on
	368	! the corners is the same as in the B grid.
	369	!
	370	!
	371	zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) &
	372	& -e2u(ji-1,jj)*u_ice(ji-1,jj) &
	373	& +e1v(ji,jj)*v_ice(ji,jj) &
	374	& -e1v(ji,jj-1)*v_ice(ji,jj-1) &
	375	& ) &
	376	& / area(ji,jj)
	377
	378	zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) &
	379	& -u_ice(ji-1,jj)/e2u(ji-1,jj) &
	380	& )e2t(ji,jj)e2t(ji,jj) &
	381	& -( v_ice(ji,jj)/e1v(ji,jj) &
	382	& -v_ice(ji,jj-1)/e1v(ji,jj-1) &
	383	& )e1t(ji,jj)e1t(ji,jj) &
	384	& ) &
	385	& / area(ji,jj)
	386
	387	!
	388	zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1) &
	389	& -u_ice(ji,jj)/e1u(ji,jj) &
	390	& )e1f(ji,jj)e1f(ji,jj) &
	391	& +( v_ice(ji+1,jj)/e2v(ji+1,jj) &
	392	& -v_ice(ji,jj)/e2v(ji,jj) &
	393	& )e2f(ji,jj)e2f(ji,jj) &
	394	& ) &
	395	& / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) &
	396	& * tmi(ji,jj) * tmi(ji,jj+1) &
	397	& * tmi(ji+1,jj) * tmi(ji+1,jj+1)
	398
	399
	400	v_ice1(ji,jj) = 0.5( (v_ice(ji,jj)+v_ice(ji,jj-1))e1t(ji+1,jj) &
	401	& +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
	402	& /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
	403
	404	u_ice2(ji,jj) = 0.5( (u_ice(ji,jj)+u_ice(ji-1,jj))e2t(ji,jj+1) &
	405	& +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &
	406	& /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
	407
	408	END DO
	409	END DO
	410
	411	CALL lbc_lnk( zdd(:,:), 'T', 1. )
	412	CALL lbc_lnk( zdt(:,:), 'T', 1. )
	413	CALL lbc_lnk( zds(:,:), 'F', 1. )
	414
	415	DO jj = k_j1+1, k_jpj-1
	416	DO ji = 2, jpim1
	417
	418	!- Calculate Delta at centre of grid cells
	419	zdst = ( e2u( ji , jj ) * v_ice1(ji,jj) &
	420	& - e2u( ji-1, jj ) * v_ice1(ji-1,jj) &
	421	& + e1v( ji , jj ) * u_ice2(ji,jj) &
	422	& - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) &
	423	& ) &
	424	& / area(ji,jj)
	425
	426	delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + &
	427	& ( zdt(ji,jj)zdt(ji,jj) + zdstzdst ) * usecc2 )
	428	deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + &
	429	(zdt(ji,jj)2 + zdst2)*usecc2), creepl )
	430
	431	!-Calculate stress tensor components zs1 and zs2
	432	!-at centre of grid cells (see section 3.5 of CICE user's guide).
	433	zs1(ji,jj) = ( zs1(ji,jj) &
	434	& - dtotel( ( 1.0 - alphaevp) zs1(ji,jj) + &
[834]	435	& ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) &
	436	* zpresh(ji,jj) ) ) &
[825]	437	& / ( 1.0 + alphaevp * dtotel )
	438
[834]	439	zs2(ji,jj) = ( zs2(ji,jj) &
	440	& - dtotel((1.0-alphaevp)ecc2*zs2(ji,jj) - &
	441	zdt(ji,jj)/deltat(ji,jj)*zpresh(ji,jj)) ) &
	442	& / ( 1.0 + alphaevpecc2dtotel )
[825]	443
	444	END DO
	445	END DO
	446
	447	CALL lbc_lnk( zs1(:,:), 'T', 1. )
	448	CALL lbc_lnk( zs2(:,:), 'T', 1. )
	449
	450	DO jj = k_j1+1, k_jpj-1
	451	DO ji = 2, jpim1
	452	!- Calculate Delta on corners
	453	zddc = ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1) &
	454	& -v_ice1(ji,jj)/e1u(ji,jj) &
	455	& )e1f(ji,jj)e1f(ji,jj) &
	456	& +( u_ice2(ji+1,jj)/e2v(ji+1,jj) &
	457	& -u_ice2(ji,jj)/e2v(ji,jj) &
	458	& )e2f(ji,jj)e2f(ji,jj) &
	459	& ) &
	460	& / ( e1f(ji,jj) * e2f(ji,jj) )
	461
	462	zdtc = (-( v_ice1(ji,jj+1)/e1u(ji,jj+1) &
	463	& -v_ice1(ji,jj)/e1u(ji,jj) &
	464	& )e1f(ji,jj)e1f(ji,jj) &
	465	& +( u_ice2(ji+1,jj)/e2v(ji+1,jj) &
	466	& -u_ice2(ji,jj)/e2v(ji,jj) &
	467	& )e2f(ji,jj)e2f(ji,jj) &
	468	& ) &
	469	& / ( e1f(ji,jj) * e2f(ji,jj) )
	470
[834]	471	deltac(ji,jj) = SQRT(zddc2+(zdtc2+zds(ji,jj)*2)usecc2) + creepl
[825]	472
	473	!-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide).
[834]	474	zs12(ji,jj) = ( zs12(ji,jj) &
	475	& - dtotel( (1.0-alphaevp)ecc2*zs12(ji,jj) - zds(ji,jj) / &
	476	& ( 2.0deltac(ji,jj) ) zpreshc(ji,jj))) &
	477	& / ( 1.0 + alphaevpecc2dtotel )
[825]	478
[834]	479	END DO ! ji
	480	END DO ! jj
[825]	481
	482	CALL lbc_lnk( zs12(:,:), 'F', 1. )
	483
	484	! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002)
	485	DO jj = k_j1+1, k_jpj-1
	486	DO ji = 2, jpim1
[834]	487	!- contribution of zs1, zs2 and zs12 to zf1
	488	zf1(ji,jj) = 0.5( (zs1(ji+1,jj)-zs1(ji,jj))e2u(ji,jj) &
	489	& +(zs2(ji+1,jj)e2t(ji+1,jj)2-zs2(ji,jj)e2t(ji,jj)**2)/e2u(ji,jj) &
	490	& +2.0(zs12(ji,jj)e1f(ji,jj)*2-zs12(ji,jj-1)e1f(ji,jj-1)**2)/e1u(ji,jj) &
	491	& ) / ( e1u(ji,jj)*e2u(ji,jj) )
	492	! contribution of zs1, zs2 and zs12 to zf2
	493	zf2(ji,jj) = 0.5( (zs1(ji,jj+1)-zs1(ji,jj))e1v(ji,jj) &
	494	& -(zs2(ji,jj+1)e1t(ji,jj+1)2 - zs2(ji,jj)e1t(ji,jj)**2)/e1v(ji,jj) &
	495	& + 2.0(zs12(ji,jj)e2f(ji,jj)**2 - &
	496	zs12(ji-1,jj)e2f(ji-1,jj)*2)/e2v(ji,jj) &
	497	& ) / ( e1v(ji,jj)*e2v(ji,jj) )
[825]	498	END DO
	499	END DO
	500	!
	501	! Computation of ice velocity
	502	!
	503	! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly.
	504	!
[834]	505	IF (MOD(iter,2).eq.0) THEN
[825]	506
	507	DO jj = k_j1+1, k_jpj-1
	508	DO ji = 2, jpim1
[834]	509	zmask = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj)
[825]	510	zsang = SIGN ( 1.0 , fcor(ji,jj) ) * sangvg
	511	z0 = zmass1(ji,jj)/dtevp
	512
	513	! SB modif because ocean has no slip boundary condition
	514	zv_ice1 = 0.5( (v_ice(ji,jj)+v_ice(ji,jj-1))e1t(ji+1,jj) &
	515	& +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
	516	& /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
[834]	517	za = rhocoSQRT((u_ice(ji,jj)-u_oce1(ji,jj))*2 + &
	518	(zv_ice1-v_oce1(ji,jj))*2) (1.0-zfrld1(ji,jj))
	519	zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
	520	za(cangvgu_oce1(ji,jj)-zsang*v_oce1(ji,jj))
[825]	521	zcca = z0+za*cangvg
	522	zccb = zcorl1(ji,jj)+za*zsang
	523	u_ice(ji,jj) = (zr+zccbzv_ice1)/(zcca+epsd)zmask
	524
	525	END DO
	526	END DO
	527
	528	CALL lbc_lnk( u_ice(:,:), 'U', -1. )
	529
	530	DO jj = k_j1+1, k_jpj-1
	531	DO ji = 2, jpim1
[834]	532
[825]	533	zmask = (1.0-max(rzero,sign(rone,-zmass2(ji,jj))))*tmv(ji,jj)
	534	zsang = sign(1.0,fcor(ji,jj))*sangvg
	535	z0 = zmass2(ji,jj)/dtevp
	536	! SB modif because ocean has no slip boundary condition
	537	zu_ice2 = 0.5( (u_ice(ji,jj)+u_ice(ji-1,jj))e2t(ji,jj+1) &
	538	& + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &
	539	& /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
[834]	540	za = rhocoSQRT((zu_ice2-u_oce2(ji,jj))*2 + &
	541	(v_ice(ji,jj)-v_oce2(ji,jj))*2)(1.0-zfrld2(ji,jj))
	542	zr = z0*v_ice(ji,jj) + zf2(ji,jj) + &
	543	za2ct(ji,jj) + za(cangvgv_oce2(ji,jj)+zsang*u_oce2(ji,jj))
[825]	544	zcca = z0+za*cangvg
	545	zccb = zcorl2(ji,jj)+za*zsang
	546	v_ice(ji,jj) = (zr-zccbzu_ice2)/(zcca+epsd)zmask
	547
	548	END DO
	549	END DO
	550
	551	CALL lbc_lnk( v_ice(:,:), 'V', -1. )
	552
[834]	553	ELSE
[825]	554	DO jj = k_j1+1, k_jpj-1
	555	DO ji = 2, jpim1
	556	zmask = (1.0-max(rzero,sign(rone,-zmass2(ji,jj))))*tmv(ji,jj)
	557	zsang = sign(1.0,fcor(ji,jj))*sangvg
	558	z0 = zmass2(ji,jj)/dtevp
	559	! SB modif because ocean has no slip boundary condition
	560	zu_ice2 = 0.5( (u_ice(ji,jj)+u_ice(ji-1,jj))e2t(ji,jj+1) &
	561	& +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &
	562	& /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
	563
[834]	564	za = rhocoSQRT((zu_ice2-u_oce2(ji,jj))*2 + &
	565	(v_ice(ji,jj)-v_oce2(ji,jj))*2)(1.0-zfrld2(ji,jj))
	566	zr = z0*v_ice(ji,jj) + zf2(ji,jj) + &
	567	za2ct(ji,jj) + za(cangvgv_oce2(ji,jj)+zsang*u_oce2(ji,jj))
[825]	568	zcca = z0+za*cangvg
	569	zccb = zcorl2(ji,jj)+za*zsang
	570	v_ice(ji,jj) = (zr-zccbzu_ice2)/(zcca+epsd)zmask
	571
	572	END DO
	573	END DO
	574
	575	CALL lbc_lnk( v_ice(:,:), 'V', -1. )
	576
	577	DO jj = k_j1+1, k_jpj-1
	578	DO ji = 2, jpim1
[834]	579	zmask = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj)
	580	zsang = SIGN(1.0,fcor(ji,jj))*sangvg
[825]	581	z0 = zmass1(ji,jj)/dtevp
	582	! SB modif because ocean has no slip boundary condition
	583	zv_ice1 = 0.5( (v_ice(ji,jj)+v_ice(ji,jj-1))e1t(ji+1,jj) &
	584	& +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
	585	& /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
	586
[834]	587	za = rhocoSQRT((u_ice(ji,jj)-u_oce1(ji,jj))*2 + &
	588	(zv_ice1-v_oce1(ji,jj))*2)(1.0-zfrld1(ji,jj))
	589	zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
	590	za(cangvgu_oce1(ji,jj)-zsang*v_oce1(ji,jj))
[825]	591	zcca = z0+za*cangvg
	592	zccb = zcorl1(ji,jj)+za*zsang
	593	u_ice(ji,jj) = (zr+zccbzv_ice1)/(zcca+epsd)zmask
[834]	594	END DO ! ji
	595	END DO ! jj
[825]	596
	597	CALL lbc_lnk( u_ice(:,:), 'U', -1. )
	598
[834]	599	ENDIF
[825]	600
[834]	601	!--- Convergence test.
	602	DO jj = k_j1+1 , k_jpj-1
	603	zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , &
	604	ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
	605	END DO
	606	zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) )
[825]	607
	608	! ! ==================== !
	609	END DO ! end loop over iter !
	610	! ! ==================== !
	611
[834]	612	!
	613	!------------------------------------------------------------------------------!
	614	! 4) Prevent ice velocities when the ice is thin
	615	!------------------------------------------------------------------------------!
	616	!
[825]	617	! If the ice thickness is below 1cm then ice velocity should equal the
[834]	618	! ocean velocity,
	619	! This prevents high velocity when ice is thin
[825]	620	DO jj = k_j1+1, k_jpj-1
	621	DO ji = 2, jpim1
	622	zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) )
	623	zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
	624	IF ( zdummy .LE. 5.0e-2 ) THEN
	625	u_ice(ji,jj) = u_io(ji,jj)
	626	v_ice(ji,jj) = v_io(ji,jj)
	627	ENDIF ! zdummy
	628	END DO
	629	END DO
[834]	630	!
	631	!------------------------------------------------------------------------------!
	632	! 5) Compute stress rain invariants and store strain tensor
	633	!------------------------------------------------------------------------------!
	634	!
	635	! Compute delta, shear and div, inputs for mechanical redistribution
[825]	636	DO jj = k_j1+1, k_jpj-1
	637	DO ji = 2, jpim1
	638	!- zdd(:,:), zdt(:,:): divergence and tension at centre
	639	!- zds(:,:): shear on northeast corner of grid cells
	640	zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) )
	641	zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
	642
	643	IF ( zdummy .LE. 5.0e-2 ) THEN
	644
	645	zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) &
	646	& -e2u(ji-1,jj)*u_ice(ji-1,jj) &
	647	& +e1v(ji,jj)*v_ice(ji,jj) &
	648	& -e1v(ji,jj-1)*v_ice(ji,jj-1) &
	649	& ) &
	650	& / area(ji,jj)
	651
	652	zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) &
	653	& -u_ice(ji-1,jj)/e2u(ji-1,jj) &
	654	& )e2t(ji,jj)e2t(ji,jj) &
	655	& -( v_ice(ji,jj)/e1v(ji,jj) &
	656	& -v_ice(ji,jj-1)/e1v(ji,jj-1) &
	657	& )e1t(ji,jj)e1t(ji,jj) &
	658	& ) &
	659	& / area(ji,jj)
	660	!
	661	! SB modif because ocean has no slip boundary condition
	662	zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) &
	663	& - u_ice(ji,jj) / e1u(ji,jj) ) &
	664	& * e1f(ji,jj) * e1f(ji,jj) &
	665	& + ( v_ice(ji+1,jj) / e2v(ji+1,jj) &
	666	& - v_ice(ji,jj) / e2v(ji,jj) ) &
	667	& * e2f(ji,jj) * e2f(ji,jj) ) &
	668	& / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) &
	669	& * tmi(ji,jj) * tmi(ji,jj+1) &
	670	& * tmi(ji+1,jj) * tmi(ji+1,jj+1)
	671
	672	!- Calculate Delta at centre of grid cells
	673	v_ice1(ji,jj) = 0.5( (v_ice(ji,jj)+v_ice(ji,jj-1))e1t(ji+1,jj) &
	674	& +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
	675	& /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
	676
	677	u_ice2(ji,jj) = 0.5( (u_ice(ji,jj)+u_ice(ji-1,jj))e2t(ji,jj+1) &
	678	& +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &
	679	& /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
	680
	681	zdst = ( e2u( ji , jj ) * v_ice1(ji,jj) &
	682	& - e2u( ji-1, jj ) * v_ice1(ji-1,jj) &
	683	& + e1v( ji , jj ) * u_ice2(ji,jj) &
	684	& - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) &
	685	& ) &
	686	& / area(ji,jj)
	687
	688	deltat(ji,jj) = SQRT( zdd(ji,jj)*zdd(ji,jj) + &
	689	& ( zdt(ji,jj)zdt(ji,jj) + zdstzdst ) * usecc2 &
	690	& ) + creepl
	691
	692	divu_i(ji,jj) = zdd(ji,jj)
	693	delta_i(ji,jj) = deltat(ji,jj)
	694	shear_i(ji,jj) = zds(ji,jj)
	695
	696	ENDIF ! zdummy
	697
	698	END DO !jj
	699	END DO !ji
	700
[834]	701	! MV This should be removed as it is the same as previous lines
	702	! ! dynamical invariants
[825]	703	! IF ( stress_mean_swi .EQ. 0 ) THEN
[834]	704	! the following should not be necessary
[825]	705	DO jj = k_j1+1, k_jpj-1
	706	DO ji = 2, jpim1
	707	divu_i(ji,jj) = zdd(ji,jj)
	708	delta_i(ji,jj) = deltat (ji,jj)
	709	shear_i(ji,jj) = zds (ji,jj)
	710	END DO
	711	END DO
	712	! ENDIF
	713
	714	CALL lbc_lnk( divu_i(:,:) , 'T', 1. )
	715	CALL lbc_lnk( delta_i(:,:), 'T', 1. )
	716	CALL lbc_lnk( shear_i(:,:), 'F', 1. )
	717
	718
[834]	719	! Store the stress tensor for next time step
	720	! this accelerates convergence and improves stability
[825]	721	IF ( stress_mean_swi .EQ. 1 ) THEN
	722	DO jj = k_j1+1, k_jpj-1
	723	DO ji = 2, jpim1
	724	stress1_i(ji,jj) = zs1(ji,jj)
	725	stress2_i(ji,jj) = zs2(ji,jj)
	726	stress12_i(ji,jj) = zs12(ji,jj)
	727	END DO
	728	END DO
	729	ENDIF
	730
[834]	731	! Lateral boundary condition
[825]	732	CALL lbc_lnk( stress1_i(:,:), 'T', 1. )
	733	CALL lbc_lnk( stress2_i(:,:), 'T', 1. )
	734	CALL lbc_lnk( stress12_i(:,:), 'F', 1. )
	735
[834]	736	!
	737	!------------------------------------------------------------------------------!
	738	! 6) Control prints of residual and charge ellipse
	739	!------------------------------------------------------------------------------!
	740	!
	741	! print the residual for convergence
	742	IF(ln_ctl) THEN
	743	WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, iter
	744	CALL prt_ctl_info(charout)
	745	CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
	746	ENDIF
[825]	747
[834]	748	! print charge ellipse
	749	! This can be desactivated once the user is sure that the stress state
	750	! lie on the charge ellipse. See Bouillon et al. 08 for more details
[825]	751	IF(ln_ctl) THEN
	752	CALL prt_ctl_info('lim_rhg : numit :',ivar1=numit)
	753	CALL prt_ctl_info('lim_rhg : nwrite :',ivar1=nwrite)
	754	CALL prt_ctl_info('lim_rhg : MOD :',ivar1=MOD(numit,nwrite))
	755	IF( MOD(numit,nwrite) .EQ. 0 ) THEN
	756	WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. '
	757	CALL prt_ctl_info(charout)
	758	DO jj = k_j1+1, k_jpj-1
	759	DO ji = 2, jpim1
	760	IF (zpresh(ji,jj) .GT. 1.0) THEN
	761	sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)*2 + 4zs12(ji,jj)2 )0.5 ) / ( 2*zpresh(ji,jj) )
	762	sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)*2 + 4zs12(ji,jj)2 )0.5 ) / ( 2*zpresh(ji,jj) )
	763	WRITE(charout,FMT="('lim_rhg :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
	764	CALL prt_ctl_info(charout)
	765	ENDIF
	766	END DO
	767	END DO
	768	WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 2000, numit, 0, 0, ' ch. ell. '
	769	CALL prt_ctl_info(charout)
	770	ENDIF
	771	ENDIF
	772
	773	END SUBROUTINE lim_rhg
	774
	775	#else
	776	!!----------------------------------------------------------------------
	777	!! Default option Dummy module NO LIM sea-ice model
	778	!!----------------------------------------------------------------------
	779	CONTAINS
	780	SUBROUTINE lim_rhg( k1 , k2 ) ! Dummy routine
	781	WRITE(,) 'lim_rhg: You should not have seen this print! error?', k1, k2
	782	END SUBROUTINE lim_rhg
	783	#endif
	784
	785	!!==============================================================================
	786	END MODULE limrhg

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

Download in other formats: