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

Context Navigation

source: branches/dev_002_LIM/NEMO/LIM_SRC_3/limrhg.F90 @ 830

Last change on this file since 830 was 825, checked in by ctlod, 16 years ago
dev_002_LIM : add the LIM 3.0 component, see ticketr: #71
File size: 37.5 KB

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

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

Download in other formats: