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

Context Navigation

source: branches/devmercator2010/NEMO/LIM_SRC_3/limrhg.F90 @ 2076

Last change on this file since 2076 was 2076, checked in by cbricaud, 14 years ago
add change dev_1784_EVP
Property svn:keywords set to `Id`
File size: 37.3 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
[2076]	9	!! - ! 2009-05 (G.Garric) addition of the lim2_evp cas
[1244]	10	!!----------------------------------------------------------------------
[2076]	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 limitd_me
	27	USE prtctl ! Print control
[2076]	28	#if defined key_lim3
	29	USE dom_ice
	30	USE ice
	31	USE iceini
	32	#endif
	33	#if defined key_lim2 && ! defined key_lim2_vp
	34	USE dom_ice_2
	35	USE ice_2
	36	USE iceini_2
	37	#endif
[825]	38
	39	IMPLICIT NONE
	40	PRIVATE
	41
	42	!! * Routine accessibility
	43	PUBLIC lim_rhg ! routine called by lim_dyn
	44
[868]	45	!! * Substitutions
	46	# include "vectopt_loop_substitute.h90"
	47
[825]	48	!! * Module variables
	49	REAL(wp) :: & ! constant values
	50	rzero = 0.e0 , &
	51	rone = 1.e0
	52	!!----------------------------------------------------------------------
[834]	53	!! LIM 3.0, UCL-LOCEAN-IPSL (2008)
[1156]	54	!! $Id$
[888]	55	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[825]	56	!!----------------------------------------------------------------------
	57
	58	CONTAINS
	59
	60	SUBROUTINE lim_rhg( k_j1, k_jpj )
[834]	61
[825]	62	!!-------------------------------------------------------------------
[834]	63	!! * SUBROUTINE lim_rhg *
	64	!! EVP-C-grid
[825]	65	!!
[834]	66	!! ** purpose : determines sea ice drift from wind stress, ice-ocean
[825]	67	!! stress and sea-surface slope. Ice-ice interaction is described by
[834]	68	!! a non-linear elasto-viscous-plastic (EVP) law including shear
	69	!! strength and a bulk rheology (Hunke and Dukowicz, 2002).
[825]	70	!!
[834]	71	!! The points in the C-grid look like this, dear reader
[825]	72	!!
[834]	73	!! (ji,jj)
	74	!! \|
	75	!! \|
	76	!! (ji-1,jj) \| (ji,jj)
	77	!! ---------
	78	!! \| \|
	79	!! \| (ji,jj) \|------(ji,jj)
	80	!! \| \|
	81	!! ---------
	82	!! (ji-1,jj-1) (ji,jj-1)
[825]	83	!!
[834]	84	!! ** Inputs : - wind forcing (stress), oceanic currents
	85	!! ice total volume (vt_i) per unit area
	86	!! snow total volume (vt_s) per unit area
[825]	87	!!
[834]	88	!! ** Action : - compute u_ice, v_ice : the components of the
	89	!! sea-ice velocity vector
	90	!! - compute delta_i, shear_i, divu_i, which are inputs
	91	!! of the ice thickness distribution
[825]	92	!!
[834]	93	!! ** Steps : 1) Compute ice snow mass, ice strength
	94	!! 2) Compute wind, oceanic stresses, mass terms and
	95	!! coriolis terms of the momentum equation
	96	!! 3) Solve the momentum equation (iterative procedure)
	97	!! 4) Prevent high velocities if the ice is thin
	98	!! 5) Recompute invariants of the strain rate tensor
	99	!! which are inputs of the ITD, store stress
	100	!! for the next time step
	101	!! 6) Control prints of residual (convergence)
	102	!! and charge ellipse.
	103	!! The user should make sure that the parameters
	104	!! nevp, telast and creepl maintain stress state
	105	!! on the charge ellipse for plastic flow
	106	!! e.g. in the Canadian Archipelago
	107	!!
	108	!! ** References : Hunke and Dukowicz, JPO97
	109	!! Bouillon et al., 08, in prep (update this when
	110	!! published)
	111	!! Vancoppenolle et al., OM08
	112	!!
[825]	113	!!-------------------------------------------------------------------
	114	! * Arguments
	115	!
	116	INTEGER, INTENT(in) :: &
[834]	117	k_j1 , & !: southern j-index for ice computation
	118	k_jpj !: northern j-index for ice computation
[825]	119
	120	! * Local variables
[834]	121	INTEGER :: ji, jj !: dummy loop indices
[825]	122
	123	INTEGER :: &
[868]	124	jter !: temporary integers
[825]	125
	126	CHARACTER (len=50) :: charout
	127
	128	REAL(wp) :: &
[834]	129	zt11, zt12, zt21, zt22, & !: temporary scalars
	130	ztagnx, ztagny, & !: wind stress on U/V points
	131	delta !
[825]	132
	133	REAL(wp) :: &
[834]	134	za, & !:
	135	zstms, & !: temporary scalar for ice strength
	136	zsang, & !: temporary scalar for coriolis term
	137	zmask !: mask for the computation of ice mass
[825]	138
[834]	139	REAL(wp),DIMENSION(jpi,jpj) :: &
	140	zpresh , & !: temporary array for ice strength
	141	zpreshc , & !: Ice strength on grid cell corners (zpreshc)
	142	zfrld1, zfrld2, & !: lead fraction on U/V points
	143	zmass1, zmass2, & !: ice/snow mass on U/V points
	144	zcorl1, zcorl2, & !: coriolis parameter on U/V points
	145	za1ct, za2ct , & !: temporary arrays
	146	zc1 , & !: ice mass
	147	zusw , & !: temporary weight for the computation
[921]	148	!: of ice strength
[834]	149	u_oce1, v_oce1, & !: ocean u/v component on U points
	150	u_oce2, v_oce2, & !: ocean u/v component on V points
	151	u_ice2, & !: ice u component on V point
	152	v_ice1 !: ice v component on U point
[825]	153
	154	REAL(wp) :: &
[834]	155	dtevp, & !: time step for subcycling
	156	dtotel, & !:
	157	ecc2, & !: square of yield ellipse eccenticity
	158	z0, & !: temporary scalar
	159	zr, & !: temporary scalar
	160	zcca, zccb, & !: temporary scalars
	161	zu_ice2, & !:
	162	zv_ice1, & !:
	163	zddc, zdtc, & !: temporary array for delta on corners
	164	zdst, & !: temporary array for delta on centre
	165	zdsshx, zdsshy, & !: term for the gradient of ocean surface
	166	sigma1, sigma2 !: internal ice stress
[825]	167
[834]	168	REAL(wp),DIMENSION(jpi,jpj) :: &
	169	zf1, zf2 !: arrays for internal stresses
[825]	170
[834]	171	REAL(wp),DIMENSION(jpi,jpj) :: &
	172	zdd, zdt, & !: Divergence and tension at centre of grid cells
	173	zds, & !: Shear on northeast corner of grid cells
	174	deltat, & !: Delta at centre of grid cells
	175	deltac, & !: Delta on corners
	176	zs1, zs2, & !: Diagonal stress tensor components zs1 and zs2
	177	zs12 !: Non-diagonal stress tensor component zs12
[825]	178
[834]	179	REAL(wp) :: &
	180	zresm , & !: Maximal error on ice velocity
	181	zindb , & !: ice (1) or not (0)
	182	zdummy !: dummy argument
[825]	183
[834]	184	REAL(wp),DIMENSION(jpi,jpj) :: &
	185	zu_ice , & !: Ice velocity on previous time step
	186	zv_ice , &
	187	zresr !: Local error on velocity
[825]	188
[2076]	189	#if defined key_lim2 && ! defined key_lim2_vp
	190	vt_s => hsnm
	191	vt_i => hicm
	192	at_i(:,:) = 1. - frld(:,:)
	193	#endif
[921]	194	!
	195	!------------------------------------------------------------------------------!
	196	! 1) Ice-Snow mass (zc1), ice strength (zpresh) !
	197	!------------------------------------------------------------------------------!
	198	!
[825]	199	! Put every vector to 0
[868]	200	zpresh(:,:) = 0.0 ; zc1(:,:) = 0.0
	201	zpreshc(:,:) = 0.0
	202	u_ice2(:,:) = 0.0 ; v_ice1(:,:) = 0.0
[825]	203	zdd(:,:) = 0.0 ; zdt(:,:) = 0.0 ; zds(:,:) = 0.0
[2076]	204	#if defined key_lim3
[834]	205	! Ice strength on T-points
[825]	206	CALL lim_itd_me_icestrength(ridge_scheme_swi)
[2076]	207	#endif
[825]	208
[834]	209	! Ice mass and temp variables
[868]	210	!CDIR NOVERRCHK
	211	DO jj = k_j1 , k_jpj
	212	!CDIR NOVERRCHK
[825]	213	DO ji = 1 , jpi
	214	zc1(ji,jj) = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) )
[2076]	215	#if defined key_lim3
[825]	216	zpresh(ji,jj) = tms(ji,jj) * strength(ji,jj) / 2.
[2076]	217	#else
	218	zpresh(ji,jj) = tms(ji,jj) * 2. * pstar * hicm(ji,jj) * EXP( -c_rhg * frld(ji,jj) )
	219	#endif
[866]	220	! tmi = 1 where there is ice or on land
	221	tmi(ji,jj) = 1.0 - ( 1.0 - MAX( 0.0 , SIGN ( 1.0 , vt_i(ji,jj) - &
[921]	222	epsd ) ) ) * tms(ji,jj)
[825]	223	END DO
	224	END DO
	225
[834]	226	! Ice strength on grid cell corners (zpreshc)
	227	! needed for calculation of shear stress
[868]	228	!CDIR NOVERRCHK
[825]	229	DO jj = k_j1+1, k_jpj-1
[868]	230	!CDIR NOVERRCHK
	231	DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1)
[921]	232	zstms = tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
	233	& tms(ji,jj+1) * wght(ji+1,jj+1,1,2) + &
	234	& tms(ji+1,jj) * wght(ji+1,jj+1,2,1) + &
	235	& tms(ji,jj) * wght(ji+1,jj+1,1,1)
	236	zusw(ji,jj) = 1.0 / MAX( zstms, epsd )
	237	zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
	238	& zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + &
	239	& zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + &
	240	& zpresh(ji,jj) * wght(ji+1,jj+1,1,1) &
	241	& ) * zusw(ji,jj)
[825]	242	END DO
	243	END DO
	244
	245	CALL lbc_lnk( zpreshc(:,:), 'F', 1. )
[921]	246	!
	247	!------------------------------------------------------------------------------!
	248	! 2) Wind / ocean stress, mass terms, coriolis terms
	249	!------------------------------------------------------------------------------!
	250	!
[825]	251	! Wind stress, coriolis and mass terms on the sides of the squares
	252	! zfrld1: lead fraction on U-points
	253	! zfrld2: lead fraction on V-points
	254	! zmass1: ice/snow mass on U-points
	255	! zmass2: ice/snow mass on V-points
	256	! zcorl1: Coriolis parameter on U-points
	257	! zcorl2: Coriolis parameter on V-points
	258	! (ztagnx,ztagny): wind stress on U/V points
	259	! u_oce1: ocean u component on u points
	260	! v_oce1: ocean v component on u points
	261	! u_oce2: ocean u component on v points
	262	! v_oce2: ocean v component on v points
[921]	263
[825]	264	DO jj = k_j1+1, k_jpj-1
[868]	265	DO ji = fs_2, fs_jpim1
[825]	266
	267	zt11 = tms(ji,jj)*e1t(ji,jj)
	268	zt12 = tms(ji+1,jj)*e1t(ji+1,jj)
	269	zt21 = tms(ji,jj)*e2t(ji,jj)
	270	zt22 = tms(ji,jj+1)*e2t(ji,jj+1)
	271
	272	! Leads area.
	273	zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + &
[921]	274	& zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd )
[825]	275	zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + &
[921]	276	& zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd )
[825]	277
	278	! Mass, coriolis coeff. and currents
[834]	279	zmass1(ji,jj) = ( zt12zc1(ji,jj) + zt11zc1(ji+1,jj) ) / (zt11+zt12+epsd)
	280	zmass2(ji,jj) = ( zt22zc1(ji,jj) + zt21zc1(ji,jj+1) ) / (zt21+zt22+epsd)
	281	zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + &
[921]	282	e1t(ji,jj)*fcor(ji+1,jj) ) &
	283	/ (e1t(ji,jj) + e1t(ji+1,jj) + epsd )
[834]	284	zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)*fcor(ji,jj) + &
[921]	285	e2t(ji,jj)*fcor(ji,jj+1) ) &
	286	/ ( e2t(ji,jj+1) + e2t(ji,jj) + epsd )
[825]	287	!
[2076]	288	! Mass, coriolis coeff. and currents
[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: