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

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 2332

Last change on this file since 2332 was 2332, checked in by rblod, 14 years ago
correct small bugs in LIM2, see ticket #746
Property svn:keywords set to `Id`
File size: 37.4 KB

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

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

Download in other formats: