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

Context Navigation

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 5059

Last change on this file since 5059 was 5053, checked in by clem, 9 years ago
LIM3 cleaning continues. No change in the physics besides the introduction of the monocategory sea ice
Property svn:keywords set to `Id`
File size: 37.7 KB

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