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

Context Navigation

source: branches/dev_1784_EVP/NEMO/LIM_SRC_3/limrhg.F90 @ 2362

Last change on this file since 2362 was 2154, checked in by cbricaud, 14 years ago
improvment for NEC
Property svn:keywords set to `Id`
File size: 36.4 KB

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

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

Download in other formats: