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

Context Navigation

source: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 6789

Last change on this file since 6789 was 6789, checked in by clem, 8 years ago
correct a bug at the coast for sea-ice velocity (and add the choice between no slip and free slip)
Property svn:keywords set to `Id`
File size: 34.4 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	!! 3.6 ! 2016-06 (C. Rousset) Rewriting and add the possibility to use mEVP from Bouillon 2013
13	!!----------------------------------------------------------------------
14	#if defined key_lim3
15	!!----------------------------------------------------------------------
16	!! 'key_lim3' LIM-3 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	USE ice ! ice variables
27	USE limitd_me ! ice strength
28	USE lbclnk ! Lateral Boundary Condition / MPP link
29	USE lib_mpp ! MPP library
30	USE wrk_nemo ! work arrays
31	USE in_out_manager ! I/O manager
32	USE prtctl ! Print control
33	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
34	#if defined key_agrif
35	USE agrif_lim3_interp
36	#endif
37	#if defined key_bdy
38	USE bdyice_lim
39	#endif
40
41	IMPLICIT NONE
42	PRIVATE
43
44	PUBLIC lim_rhg ! routine called by lim_dyn
45
46	!! * Substitutions
47	# include "vectopt_loop_substitute.h90"
48	!!----------------------------------------------------------------------
49	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
50	!! $Id$
51	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
52	!!----------------------------------------------------------------------
53	CONTAINS
54
55	SUBROUTINE lim_rhg
56	!!-------------------------------------------------------------------
57	!! * SUBROUTINE lim_rhg *
58	!! EVP-C-grid
59	!!
60	!! ** purpose : determines sea ice drift from wind stress, ice-ocean
61	!! stress and sea-surface slope. Ice-ice interaction is described by
62	!! a non-linear elasto-viscous-plastic (EVP) law including shear
63	!! strength and a bulk rheology (Hunke and Dukowicz, 2002).
64	!!
65	!! The points in the C-grid look like this, dear reader
66	!!
67	!! (ji,jj)
68	!! \|
69	!! \|
70	!! (ji-1,jj) \| (ji,jj)
71	!! ---------
72	!! \| \|
73	!! \| (ji,jj) \|------(ji,jj)
74	!! \| \|
75	!! ---------
76	!! (ji-1,jj-1) (ji,jj-1)
77	!!
78	!! ** Inputs : - wind forcing (stress), oceanic currents
79	!! ice total volume (vt_i) per unit area
80	!! snow total volume (vt_s) per unit area
81	!!
82	!! ** Action : - compute u_ice, v_ice : the components of the
83	!! sea-ice velocity vector
84	!! - compute delta_i, shear_i, divu_i, which are inputs
85	!! of the ice thickness distribution
86	!!
87	!! ** Steps : 1) Compute ice snow mass, ice strength
88	!! 2) Compute wind, oceanic stresses, mass terms and
89	!! coriolis terms of the momentum equation
90	!! 3) Solve the momentum equation (iterative procedure)
91	!! 4) Prevent high velocities if the ice is thin
92	!! 5) Recompute invariants of the strain rate tensor
93	!! which are inputs of the ITD, store stress
94	!! for the next time step
95	!! 6) Control prints of residual (convergence)
96	!! and charge ellipse.
97	!! The user should make sure that the parameters
98	!! nn_nevp, elastic time scale and rn_creepl maintain stress state
99	!! on the charge ellipse for plastic flow
100	!! e.g. in the Canadian Archipelago
101	!!
102	!! ** Notes : There is the possibility to use mEVP from Bouillon 2013
103	!! (by uncommenting some lines in part 3 and changing alpha and beta parameters)
104	!! but this solution appears very unstable (see Kimmritz et al 2016)
105	!!
106	!! References : Hunke and Dukowicz, JPO97
107	!! Bouillon et al., Ocean Modelling 2009
108	!! Bouillon et al., Ocean Modelling 2013
109	!!-------------------------------------------------------------------
110	INTEGER :: ji, jj ! dummy loop indices
111	INTEGER :: jter ! local integers
112	CHARACTER (len=50) :: charout
113
114	REAL(wp) :: dtevp, z1_dtevp ! time step for subcycling
115	REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity
116	REAL(wp) :: zbeta, zalph1, z1_alph1, zalph2, z1_alph2 ! alpha and beta from Bouillon 2009 and 2013
117	REAL(wp) :: zm1, zm2, zm3 ! ice/snow mass
118	REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars
119	REAL(wp) :: zTauO, zTauA, zTauB, ZCor, zmt, zu_ice2, zv_ice1, zvel ! temporary scalars
120
121	REAL(wp) :: zsig1, zsig2 ! internal ice stress
122	REAL(wp) :: zresm ! Maximal error on ice velocity
123	REAL(wp) :: zintb, zintn ! dummy argument
124
125	REAL(wp), POINTER, DIMENSION(:,:) :: zpresh ! temporary array for ice strength
126	REAL(wp), POINTER, DIMENSION(:,:) :: z1_e1t0, z1_e2t0 ! scale factors
127	REAL(wp), POINTER, DIMENSION(:,:) :: zp_delt ! P/delta at T points
128	!
129	REAL(wp), POINTER, DIMENSION(:,:) :: za1 , za2 ! ice fraction on U/V points
130	REAL(wp), POINTER, DIMENSION(:,:) :: zmass1, zmass2 ! ice/snow mass on U/V points
131	REAL(wp), POINTER, DIMENSION(:,:) :: zcor1 , zcor2 ! coriolis parameter on U/V points
132	REAL(wp), POINTER, DIMENSION(:,:) :: zspgu , zspgv ! surface pressure gradient at U/V points
133	REAL(wp), POINTER, DIMENSION(:,:) :: v_oce1, u_oce2, v_ice1, u_ice2 ! ocean/ice u/v component on V/U points
134	REAL(wp), POINTER, DIMENSION(:,:) :: zf1 , zf2 ! internal stresses
135
136	REAL(wp), POINTER, DIMENSION(:,:) :: zds ! shear
137	REAL(wp), POINTER, DIMENSION(:,:) :: zs1, zs2, zs12 ! stress tensor components
138	REAL(wp), POINTER, DIMENSION(:,:) :: zu_ice, zv_ice, zresr ! check convergence
139	REAL(wp), POINTER, DIMENSION(:,:) :: zpice ! array used for the calculation of ice surface slope:
140	! ocean surface (ssh_m) if ice is not embedded
141	! ice top surface if ice is embedded
142	REAL(wp), POINTER, DIMENSION(:,:) :: zswitch1, zswitch2 ! dummy arrays
143	REAL(wp), POINTER, DIMENSION(:,:) :: zfmask, zwf ! mask at F points for the ice
144
145	REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter
146	REAL(wp), PARAMETER :: zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity
147	REAL(wp), PARAMETER :: zshlat = 2._wp ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip)
148	!!-------------------------------------------------------------------
149
150	CALL wrk_alloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt )
151	CALL wrk_alloc( jpi,jpj, za1, za2, zmass1, zmass2, zcor1, zcor2 )
152	CALL wrk_alloc( jpi,jpj, zspgu, zspgv, v_oce1, u_oce2, v_ice1, u_ice2, zf1, zf2 )
153	CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
154	CALL wrk_alloc( jpi,jpj, zswitch1, zswitch2, zfmask, zwf )
155
156	#if defined key_agrif
157	CALL agrif_interp_lim3('U') ! First interpolation of coarse values
158	CALL agrif_interp_lim3('V')
159	#endif
160	!
161	!------------------------------------------------------------------------------!
162	! 0) mask at F points for the ice
163	!------------------------------------------------------------------------------!
164	! ocean/land mask
165	DO jj = 1, jpjm1
166	DO ji = 1, jpim1 ! NO vector opt.
167	zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
168	END DO
169	END DO
170	CALL lbc_lnk( zfmask, 'F', 1._wp )
171
172	! Lateral boundary conditions on velocity (modify zfmask)
173	zwf(:,:) = zfmask(:,:)
174	DO jj = 2, jpjm1
175	DO ji = fs_2, fs_jpim1 ! vector opt.
176	IF( zfmask(ji,jj) == 0._wp ) THEN
177	zfmask(ji,jj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
178	ENDIF
179	END DO
180	END DO
181	DO jj = 2, jpjm1
182	IF( zfmask(1,jj) == 0._wp ) THEN
183	zfmask(1 ,jj) = zshlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
184	ENDIF
185	IF( zfmask(jpi,jj) == 0._wp ) THEN
186	zfmask(jpi,jj) = zshlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
187	ENDIF
188	END DO
189	DO ji = 2, jpim1
190	IF( zfmask(ji,1) == 0._wp ) THEN
191	zfmask(ji,1 ) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
192	ENDIF
193	IF( zfmask(ji,jpj) == 0._wp ) THEN
194	zfmask(ji,jpj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
195	ENDIF
196	END DO
197	CALL lbc_lnk( zfmask, 'F', 1._wp )
198
199	!------------------------------------------------------------------------------!
200	! 1) define some variables and initialize arrays
201	!------------------------------------------------------------------------------!
202	! ecc2: square of yield ellipse eccenticrity
203	ecc2 = rn_ecc * rn_ecc
204	z1_ecc2 = 1._wp / ecc2
205
206	! Time step for subcycling
207	dtevp = rdt_ice / REAL( nn_nevp )
208	z1_dtevp = 1._wp / dtevp
209
210	! alpha parameters (Bouillon 2009)
211	zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
212	zalph2 = zalph1 * z1_ecc2
213
214	! alpha and beta parameters (Bouillon 2013)
215	!!zalph1 = 40.
216	!!zalph2 = 40.
217	!!zbeta = 3000.
218	!!zbeta = REAL( nn_nevp ) ! close to classical EVP of Hunke (2001)
219
220	z1_alph1 = 1._wp / ( zalph1 + 1._wp )
221	z1_alph2 = 1._wp / ( zalph2 + 1._wp )
222
223	! Initialise stress tensor
224	zs1 (:,:) = stress1_i (:,:)
225	zs2 (:,:) = stress2_i (:,:)
226	zs12(:,:) = stress12_i(:,:)
227
228	! Ice strength
229	CALL lim_itd_me_icestrength( nn_icestr )
230	zpresh(:,:) = tmask(:,:,1) * strength(:,:)
231
232	! scale factors
233	DO jj = 2, jpjm1
234	DO ji = fs_2, fs_jpim1
235	z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj ) + e1t(ji,jj ) )
236	z1_e2t0(ji,jj) = 1._wp / ( e2t(ji ,jj+1) + e2t(ji,jj ) )
237	END DO
238	END DO
239
240	!
241	!------------------------------------------------------------------------------!
242	! 2) Wind / ocean stress, mass terms, coriolis terms
243	!------------------------------------------------------------------------------!
244
245	IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==!
246	!
247	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
248	! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
249	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
250	!
251	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
252	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
253	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
254	!
255	zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
256	!
257	ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==!
258	zpice(:,:) = ssh_m(:,:)
259	ENDIF
260
261	DO jj = 2, jpjm1
262	DO ji = fs_2, fs_jpim1
263
264	! ice fraction at U-V points
265	za1(ji,jj) = ( at_i(ji,jj) * e1t(ji+1,jj) + at_i(ji+1,jj) * e1t(ji,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
266	za2(ji,jj) = ( at_i(ji,jj) * e2t(ji,jj+1) + at_i(ji,jj+1) * e2t(ji,jj) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
267
268	! Ice/snow mass at U-V points
269	zm1 = ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) )
270	zm2 = ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) )
271	zm3 = ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) )
272	zmass1(ji,jj) = ( zm1 * e1t(ji+1,jj) + zm2 * e1t(ji,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
273	zmass2(ji,jj) = ( zm1 * e2t(ji,jj+1) + zm3 * e2t(ji,jj) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
274
275	! Ocean currents at U-V points
276	v_oce1(ji,jj) = 0.5_wp * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji+1,jj) &
277	& + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
278
279	u_oce2(ji,jj) = 0.5_wp * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj+1) &
280	& + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
281
282	! Coriolis at U-V points (m*f)
283	zcor1(ji,jj) = zmass1(ji,jj) * 0.5_wp * ( ff(ji,jj) + ff(ji ,jj-1) )
284	zcor2(ji,jj) = - zmass2(ji,jj) * 0.5_wp * ( ff(ji,jj) + ff(ji-1,jj ) )
285
286	! Surface pressure gradient (- mgGRAD(ssh)) at U-V points
287	zspgu(ji,jj) = - zmass1(ji,jj) * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
288	zspgv(ji,jj) = - zmass2(ji,jj) * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
289
290	! switches
291	zswitch1(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) )
292	zswitch2(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) )
293
294	END DO
295	END DO
296	!
297	!------------------------------------------------------------------------------!
298	! 3) Solution of the momentum equation, iterative procedure
299	!------------------------------------------------------------------------------!
300	!
301	! !----------------------!
302	DO jter = 1 , nn_nevp ! loop over jter !
303	! !----------------------!
304	IF(ln_ctl) THEN ! Convergence test
305	DO jj = 1, jpjm1
306	zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
307	zv_ice(:,jj) = v_ice(:,jj)
308	END DO
309	ENDIF
310
311	! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
312	DO jj = 1, jpjm1 ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
313	DO ji = 1, fs_jpim1
314
315	! shear at F points
316	zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &
317	& + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
318	& ) * r1_e12f(ji,jj) * zfmask(ji,jj)
319
320	END DO
321	END DO
322	CALL lbc_lnk( zds, 'F', 1. )
323
324	DO jj = 2, jpjm1
325	DO ji = 2, jpim1 ! no vector loop
326
327	! shear**2 at T points (doc eq. A16)
328	zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e12f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e12f(ji-1,jj ) &
329	& + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1) &
330	& ) * 0.25_wp * r1_e12t(ji,jj)
331
332	! divergence at T points
333	zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
334	& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
335	& ) * r1_e12t(ji,jj)
336	zdiv2 = zdiv * zdiv
337
338	! tension at T points
339	zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) &
340	& - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
341	& ) * r1_e12t(ji,jj)
342	zdt2 = zdt * zdt
343
344	! delta at T points
345	zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * usecc2 )
346
347	! P/delta at T points
348	zp_delt(ji,jj) = zpresh(ji,jj) / ( zdelta + rn_creepl )
349
350	! stress at T points
351	zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
352	zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
353
354	END DO
355	END DO
356	CALL lbc_lnk( zp_delt, 'T', 1. )
357
358	DO jj = 1, jpjm1
359	DO ji = 1, jpim1
360
361	! P/delta at F points
362	zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )
363
364	! stress at F points
365	zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
366
367	END DO
368	END DO
369	CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
370
371
372	! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
373	DO jj = 2, jpjm1
374	DO ji = fs_2, fs_jpim1
375
376	! U points
377	zf1(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) &
378	& + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) &
379	& ) * r1_e2u(ji,jj) &
380	& + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) &
381	& ) * 2._wp * r1_e1u(ji,jj) &
382	& ) * r1_e12u(ji,jj)
383
384	! V points
385	zf2(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) &
386	& - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) &
387	& ) * r1_e1v(ji,jj) &
388	& + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) &
389	& ) * 2._wp * r1_e2v(ji,jj) &
390	& ) * r1_e12v(ji,jj)
391
392	! u_ice at V point
393	u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) &
394	& + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
395
396	! v_ice at U point
397	v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) &
398	& + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
399
400	END DO
401	END DO
402	!
403	! --- Computation of ice velocity --- !
404	! Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen smart and large nn_nevp
405	! Bouillon et al. 2009 (eq 34-35) => stable
406	IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations
407
408	DO jj = 2, jpjm1
409	DO ji = fs_2, fs_jpim1
410
411	zvel = SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_ice2(ji,jj) * u_ice2(ji,jj) )
412
413	zTauA = za2(ji,jj) * vtau_ice(ji,jj)
414	zTauO = za2(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) &
415	& + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) )
416	ztauB = - tau_icebfr(ji,jj) / MAX( zvel, zepsi )
417	zCor = zcor2(ji,jj) * u_ice2(ji,jj)
418	zmt = zmass2(ji,jj) * z1_dtevp
419
420	! Bouillon 2009
421	v_ice(ji,jj) = ( zmt * v_ice(ji,jj) + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj) &
422	& ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch2(ji,jj)
423	! Bouillon 2013
424	!!v_ice(ji,jj) = ( zmt * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) ) &
425	!! & + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj) &
426	!! & ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj)
427
428	END DO
429	END DO
430	CALL lbc_lnk( v_ice, 'V', -1. )
431
432	#if defined key_agrif
433	CALL agrif_interp_lim3( 'V' )
434	#endif
435	#if defined key_bdy
436	CALL bdy_ice_lim_dyn( 'V' )
437	#endif
438
439	DO jj = 2, jpjm1
440	DO ji = fs_2, fs_jpim1
441
442	! v_ice at U point
443	zv_ice1 = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) &
444	& + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
445
446	zvel = SQRT( v_ice1(ji,jj) * v_ice1(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
447
448	zTauA = za1(ji,jj) * utau_ice(ji,jj)
449	zTauO = za1(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) &
450	& + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) )
451	ztauB = - tau_icebfr(ji,jj) / MAX( zvel, zepsi )
452	zCor = zcor1(ji,jj) * zv_ice1
453	zmt = zmass1(ji,jj) * z1_dtevp
454
455	! Bouillon 2009
456	u_ice(ji,jj) = ( zmt * u_ice(ji,jj) + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj) &
457	& ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch1(ji,jj)
458	! Bouillon 2013
459	!!u_ice(ji,jj) = ( zmt * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) ) &
460	!! & + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj) &
461	!! & ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch1(ji,jj)
462	END DO
463	END DO
464	CALL lbc_lnk( u_ice, 'U', -1. )
465
466	#if defined key_agrif
467	CALL agrif_interp_lim3( 'U' )
468	#endif
469	#if defined key_bdy
470	CALL bdy_ice_lim_dyn( 'U' )
471	#endif
472
473	ELSE ! odd iterations
474
475	DO jj = 2, jpjm1
476	DO ji = fs_2, fs_jpim1
477
478	zvel = SQRT( v_ice1(ji,jj) * v_ice1(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
479
480	zTauA = za1(ji,jj) * utau_ice(ji,jj)
481	zTauO = za1(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) &
482	& + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) )
483	ztauB = - tau_icebfr(ji,jj) / MAX( zvel, zepsi )
484	zCor = zcor1(ji,jj) * v_ice1(ji,jj)
485	zmt = zmass1(ji,jj) * z1_dtevp
486
487	! Bouillon 2009
488	u_ice(ji,jj) = ( zmt * u_ice(ji,jj) + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj) &
489	& ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch1(ji,jj)
490	! Bouillon 2013
491	!!u_ice(ji,jj) = ( zmt * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) ) &
492	!! & + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj) &
493	!! & ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch1(ji,jj)
494	END DO
495	END DO
496	CALL lbc_lnk( u_ice, 'U', -1. )
497
498	#if defined key_agrif
499	CALL agrif_interp_lim3( 'U' )
500	#endif
501	#if defined key_bdy
502	CALL bdy_ice_lim_dyn( 'U' )
503	#endif
504
505	DO jj = 2, jpjm1
506	DO ji = fs_2, fs_jpim1
507
508	! u_ice at V point
509	zu_ice2 = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) &
510	& + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
511
512	zvel = SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_ice2(ji,jj) * u_ice2(ji,jj) )
513
514	zTauA = za2(ji,jj) * vtau_ice(ji,jj)
515	zTauO = za2(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) &
516	& + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) )
517	ztauB = - tau_icebfr(ji,jj) / MAX( zvel, zepsi )
518	zCor = zcor2(ji,jj) * zu_ice2
519	zmt = zmass2(ji,jj) * z1_dtevp
520
521	! Bouillon 2009
522	v_ice(ji,jj) = ( zmt * v_ice(ji,jj) + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj) &
523	& ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch2(ji,jj)
524	! Bouillon 2013
525	!!v_ice(ji,jj) = ( zmt * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) ) &
526	!! & + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj) &
527	!! & ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj)
528	END DO
529	END DO
530	CALL lbc_lnk( v_ice, 'V', -1. )
531
532	#if defined key_agrif
533	CALL agrif_interp_lim3( 'V' )
534	#endif
535	#if defined key_bdy
536	CALL bdy_ice_lim_dyn( 'V' )
537	#endif
538
539	ENDIF
540
541	IF(ln_ctl) THEN ! Convergence test
542	DO jj = 2 , jpjm1
543	zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
544	END DO
545	zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
546	IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain
547	ENDIF
548	!
549	! ! ==================== !
550	END DO ! end loop over jter !
551	! ! ==================== !
552	!
553	!------------------------------------------------------------------------------!
554	! 4) Limit ice velocities when ice is thin
555	!------------------------------------------------------------------------------!
556	! If the ice volume is below zvmin then ice velocity should equal the
557	! ocean velocity. This prevents high velocity when ice is thin
558	DO jj = 2, jpjm1
559	DO ji = fs_2, fs_jpim1
560	rswitch = MAX( 0._wp, SIGN( 1._wp, vt_i(ji,jj) - zvmin ) )
561	u_ice(ji,jj) = ( rswitch * u_ice(ji,jj) + (1._wp - rswitch) * u_oce(ji,jj) ) * zswitch1(ji,jj)
562	v_ice(ji,jj) = ( rswitch * v_ice(ji,jj) + (1._wp - rswitch) * v_oce(ji,jj) ) * zswitch2(ji,jj)
563	END DO
564	END DO
565
566	CALL lbc_lnk_multi( u_ice, 'U', -1., v_ice, 'V', -1. )
567
568	#if defined key_agrif
569	CALL agrif_interp_lim3( 'U' )
570	CALL agrif_interp_lim3( 'V' )
571	#endif
572	#if defined key_bdy
573	CALL bdy_ice_lim_dyn( 'U' )
574	CALL bdy_ice_lim_dyn( 'V' )
575	#endif
576
577	!------------------------------------------------------------------------------!
578	! 5) Recompute delta, shear and div (inputs for mechanical redistribution)
579	!------------------------------------------------------------------------------!
580	DO jj = 1, jpjm1
581	DO ji = 1, fs_jpim1
582
583	! shear at F points
584	zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &
585	& + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
586	& ) * r1_e12f(ji,jj) * zfmask(ji,jj)
587
588	END DO
589	END DO
590	CALL lbc_lnk( zds, 'F', 1. )
591
592	DO jj = 2, jpjm1
593	DO ji = 2, jpim1 ! no vector loop
594
595	! tension**2 at T points
596	zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) &
597	& - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
598	& ) * r1_e12t(ji,jj)
599	zdt2 = zdt * zdt
600
601	! shear**2 at T points (doc eq. A16)
602	zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e12f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e12f(ji-1,jj ) &
603	& + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1) &
604	& ) * 0.25_wp * r1_e12t(ji,jj)
605
606	! shear at T points
607	shear_i(ji,jj) = SQRT( zdt2 + zds2 )
608
609	! divergence at T points
610	divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
611	& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
612	& ) * r1_e12t(ji,jj)
613
614	! delta at T points
615	zdelta = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * usecc2 )
616	delta_i(ji,jj) = zdelta + rn_creepl
617
618	END DO
619	END DO
620	CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. )
621
622	! --- Store the stress tensor for the next time step --- !
623	stress1_i (:,:) = zs1 (:,:)
624	stress2_i (:,:) = zs2 (:,:)
625	stress12_i(:,:) = zs12(:,:)
626	!
627
628	!------------------------------------------------------------------------------!
629	! 6) Control prints of residual and charge ellipse
630	!------------------------------------------------------------------------------!
631	!
632	! print the residual for convergence
633	IF(ln_ctl) THEN
634	WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, jter
635	CALL prt_ctl_info(charout)
636	CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
637	ENDIF
638
639	! print charge ellipse
640	! This can be desactivated once the user is sure that the stress state
641	! lie on the charge ellipse. See Bouillon et al. 08 for more details
642	IF(ln_ctl) THEN
643	CALL prt_ctl_info('lim_rhg : numit :',ivar1=numit)
644	CALL prt_ctl_info('lim_rhg : nwrite :',ivar1=nwrite)
645	CALL prt_ctl_info('lim_rhg : MOD :',ivar1=MOD(numit,nwrite))
646	IF( MOD(numit,nwrite) .EQ. 0 ) THEN
647	WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. '
648	CALL prt_ctl_info(charout)
649	DO jj = 2, jpjm1
650	DO ji = 2, jpim1
651	IF (zpresh(ji,jj) > 1.0) THEN
652	zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)*2 + 4zs12(ji,jj)*2 ) ) / ( 2zpresh(ji,jj) )
653	zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)*2 + 4zs12(ji,jj)*2 ) ) / ( 2zpresh(ji,jj) )
654	WRITE(charout,FMT="('lim_rhg :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
655	CALL prt_ctl_info(charout)
656	ENDIF
657	END DO
658	END DO
659	WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 2000, numit, 0, 0, ' ch. ell. '
660	CALL prt_ctl_info(charout)
661	ENDIF
662	ENDIF
663	!
664	CALL wrk_dealloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt )
665	CALL wrk_dealloc( jpi,jpj, za1, za2, zmass1, zmass2, zcor1, zcor2 )
666	CALL wrk_dealloc( jpi,jpj, zspgu, zspgv, v_oce1, u_oce2, v_ice1, u_ice2, zf1, zf2 )
667	CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
668	CALL wrk_dealloc( jpi,jpj, zswitch1, zswitch2, zfmask, zwf )
669
670	END SUBROUTINE lim_rhg
671
672	#else
673	!!----------------------------------------------------------------------
674	!! Default option Dummy module NO LIM sea-ice model
675	!!----------------------------------------------------------------------
676	CONTAINS
677	SUBROUTINE lim_rhg ! Dummy routine
678	WRITE(,) 'lim_rhg: You should not have seen this print! error?'
679	END SUBROUTINE lim_rhg
680	#endif
681
682	!!==============================================================================
683	END MODULE limrhg

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

Download in other formats: