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_2.F90 in branches/2011/dev_r2802_LOCEAN10_agrif_lim/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

Context Navigation

source: branches/2011/dev_r2802_LOCEAN10_agrif_lim/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90 @ 2804

Last change on this file since 2804 was 2804, checked in by rblod, 13 years ago
dev_r2802_LOCEAN10_agrif_lim: first implementation see ticket #848
Property svn:keywords set to `Id`
File size: 35.2 KB

Line
1	MODULE limrhg_2
2	!!======================================================================
3	!! * MODULE limrhg_2 *
4	!! Ice rheology : performs sea ice rheology
5	!!======================================================================
6	!! History : 0.0 ! 1993-12 (M.A. Morales Maqueda.) Original code
7	!! 1.0 ! 1994-12 (H. Goosse)
8	!! 2.0 ! 2003-12 (C. Ethe, G. Madec) F90, mpp
9	!! - ! 2006-08 (G. Madec) surface module, ice-stress at I-point
10	!! - ! 2009-09 (G. Madec) Huge verctor optimisation
11	!! 3.3 ! 2009-05 (G.Garric, C. Bricaud) addition of the lim2_evp case
12	!!----------------------------------------------------------------------
13	#if defined key_lim2 && defined key_lim2_vp
14	!!----------------------------------------------------------------------
15	!! 'key_lim2' AND LIM-2 sea-ice model
16	!! 'key_lim2_vp' VP ice rheology
17	!!----------------------------------------------------------------------
18	!! lim_rhg_2 : computes ice velocities
19	!!----------------------------------------------------------------------
20	USE par_oce ! ocean parameter
21	USE dom_oce ! ocean space and time domain
22	USE sbc_oce ! surface boundary condition: ocean variables
23	USE sbc_ice ! surface boundary condition: ice variables
24	USE dom_ice_2 ! LIM2: ice domain
25	USE phycst ! physical constant
26	USE ice_2 ! LIM2: ice variables
27	USE lbclnk ! lateral boundary condition - MPP exchanges
28	USE lib_mpp ! MPP library
29	USE in_out_manager ! I/O manager
30	USE prtctl ! Print control
31	#if defined key_agrif
32	USE agrif_lim2_interp ! nesting
33	#endif
34
35	IMPLICIT NONE
36	PRIVATE
37
38	PUBLIC lim_rhg_2 ! routine called by lim_dyn
39	PUBLIC lim_rhg_alloc_2 ! routine called by lim_dyn_alloc_2
40
41	REAL(wp) :: rzero = 0._wp ! constant value: zero
42	REAL(wp) :: rone = 1._wp ! and one
43
44	! 2D workspaces for lim_rhg_2. Can't use wrk_nemo module for them because
45	! extent in 2nd dimension is > jpj.
46	REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:) :: zu0, zv0
47	REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:) :: zu_n, zv_n
48	REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:) :: zu_a, zv_a
49	REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:) :: zviszeta, zviseta
50	REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:) :: zzfrld, zztms
51	REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:) :: zi1, zi2, zmasst, zpresh
52
53	!! * Substitutions
54	# include "vectopt_loop_substitute.h90"
55	!!----------------------------------------------------------------------
56	!! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
57	!! $Id$
58	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
59	!!----------------------------------------------------------------------
60	CONTAINS
61
62	INTEGER FUNCTION lim_rhg_alloc_2()
63	!!-------------------------------------------------------------------
64	!! * FUNCTION lim_rhg_alloc_2 *
65	!!-------------------------------------------------------------------
66	ALLOCATE( zu0(jpi,0:jpj+1), zv0(jpi,0:jpj+1), &
67	& zu_n(jpi,0:jpj+1), zv_n(jpi,0:jpj+1), &
68	& zu_a(jpi,0:jpj+1), zv_a(jpi,0:jpj+1), &
69	& zviszeta(jpi,0:jpj+1), zviseta(jpi,0:jpj+1), &
70	& zzfrld(jpi,0:jpj+1), zztms(jpi,0:jpj+1), &
71	& zi1(jpi,0:jpj+1), zi2(jpi,0:jpj+1), &
72	& zmasst(jpi,0:jpj+1), zpresh(jpi,0:jpj+1), &
73	& Stat=lim_rhg_alloc_2)
74	!
75	IF( lim_rhg_alloc_2 /= 0 ) CALL ctl_warn('lim_rhg_alloc_2 : failed to allocate arrays')
76	!
77	END FUNCTION lim_rhg_alloc_2
78
79
80	SUBROUTINE lim_rhg_2( k_j1, k_jpj )
81	!!-------------------------------------------------------------------
82	!! * SUBROUTINR lim_rhg_2 *
83	!!
84	!! ** purpose : determines the velocity field of sea ice by using
85	!! atmospheric (wind stress) and oceanic (water stress and surface
86	!! tilt) forcings. Ice-ice interaction is described by a non-linear
87	!! viscous-plastic law including shear strength and a bulk rheology.
88	!!
89	!! ** Action : - compute u_ice, v_ice the sea-ice velocity defined
90	!! at I-point
91	!!-------------------------------------------------------------------
92	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
93	USE wrk_nemo, ONLY: zfrld => wrk_2d_1, zmass => wrk_2d_2, zcorl => wrk_2d_3
94	USE wrk_nemo, ONLY: za1ct => wrk_2d_4, za2ct => wrk_2d_5, zresr => wrk_2d_6
95	USE wrk_nemo, ONLY: zc1u => wrk_2d_7, zc1v => wrk_2d_8, zc2u => wrk_2d_9
96	USE wrk_nemo, ONLY: zc2v => wrk_2d_10, zsang => wrk_2d_11
97	!!
98	INTEGER, INTENT(in) :: k_j1 ! southern j-index for ice computation
99	INTEGER, INTENT(in) :: k_jpj ! northern j-index for ice computation
100	!!
101	INTEGER :: ji, jj ! dummy loop indices
102	INTEGER :: iter, jter ! temporary integers
103	CHARACTER (len=50) :: charout
104	REAL(wp) :: ze11 , ze12 , ze22 , ze21 ! local scalars
105	REAL(wp) :: zt11 , zt12 , zt21 , zt22 ! - -
106	REAL(wp) :: zvis11, zvis21, zvis12, zvis22 ! - -
107	REAL(wp) :: zgphsx, ztagnx, zgsshx, zunw, zur, zusw ! - -
108	REAL(wp) :: zgphsy, ztagny, zgsshy, zvnw, zvr ! - -
109	REAL(wp) :: zresm, za, zac, zmod
110	REAL(wp) :: zmpzas, zstms, zindu, zusdtp, zmassdt, zcorlal
111	REAL(wp) :: ztrace2, zdeter, zdelta, zmask, zdgp, zdgi, zdiag
112	REAL(wp) :: za1, zb1, zc1, zd1
113	REAL(wp) :: za2, zb2, zc2, zd2, zden
114	REAL(wp) :: zs11_11, zs11_12, zs11_21, zs11_22
115	REAL(wp) :: zs12_11, zs12_12, zs12_21, zs12_22
116	REAL(wp) :: zs21_11, zs21_12, zs21_21, zs21_22
117	REAL(wp) :: zs22_11, zs22_12, zs22_21, zs22_22
118	!!-------------------------------------------------------------------
119
120	! Store initial velocities
121	! ----------------
122	zztms(:,0 ) = 0._wp ; zzfrld(:,0 ) = 0._wp
123	zztms(:,jpj+1) = 0._wp ; zzfrld(:,jpj+1) = 0._wp
124	zu0 (:,0 ) = 0._wp ; zv0 (:,0 ) = 0._wp
125	zu0 (:,jpj+1) = 0._wp ; zv0 (:,jpj+1) = 0._wp
126	zztms(:,1:jpj) = tms (:,:) ; zzfrld(:,1:jpj) = frld (:,:)
127	zu0 (:,1:jpj) = u_ice(:,:) ; zv0 (:,1:jpj) = v_ice(:,:)
128	zu_a (:, : ) = zu0 (:,:) ; zv_a (:, : ) = zv0 (:,:)
129	zu_n (:, : ) = zu0 (:,:) ; zv_n (:, : ) = zv0 (:,:)
130
131	!i
132	zi1 (:,:) = 0._wp
133	zi2 (:,:) = 0._wp
134	zpresh(:,:) = 0._wp
135	zmasst(:,:) = 0._wp
136	!i
137	!!gm violant
138	zfrld(:,:) =0._wp
139	zcorl(:,:) =0._wp
140	zmass(:,:) =0._wp
141	za1ct(:,:) =0._wp
142	za2ct(:,:) =0._wp
143	!!gm end
144
145	zviszeta(:,:) = 0._wp
146	zviseta (:,:) = 0._wp
147
148	!i zviszeta(:,0 ) = 0._wp ; zviseta(:,0 ) = 0._wp
149	!i zviszeta(:,jpj ) = 0._wp ; zviseta(:,jpj ) = 0._wp
150	!i zviszeta(:,jpj+1) = 0._wp ; zviseta(:,jpj+1) = 0._wp
151
152	#if defined key_agrif
153	! load the boundary value of velocity in special array zuive and zvice
154	CALL Agrif_dyn_lim2_load
155	#endif
156
157	! Ice mass, ice strength, and wind stress at the center \|
158	! of the grid squares. \|
159	!-------------------------------------------------------------------
160
161	!CDIR NOVERRCHK
162	DO jj = k_j1 , k_jpj-1
163	!CDIR NOVERRCHK
164	DO ji = 1 , jpi
165	! only the sinus changes its sign with the hemisphere
166	zsang(ji,jj) = SIGN( 1._wp, fcor(ji,jj) ) * sangvg ! only the sinus changes its sign with the hemisphere
167	!
168	zmasst(ji,jj) = tms(ji,jj) * ( rhosn * hsnm(ji,jj) + rhoic * hicm(ji,jj) )
169	zpresh(ji,jj) = tms(ji,jj) * pstarh * hicm(ji,jj) * EXP( -c_rhg * frld(ji,jj) )
170	!!gm :: stress given at I-point (F-point for the ocean) only compute the ponderation with the ice fraction (1-frld)
171	zi1(ji,jj) = tms(ji,jj) * ( 1._wp - frld(ji,jj) )
172	zi2(ji,jj) = tms(ji,jj) * ( 1._wp - frld(ji,jj) )
173	END DO
174	END DO
175
176
177	!---------------------------------------------------------------------------
178	! Wind stress, coriolis and mass terms at the corners of the grid squares \|
179	! Gradient of ice strenght. \|
180	!---------------------------------------------------------------------------
181
182	DO jj = k_j1+1, k_jpj-1
183	DO ji = 2, jpi ! NO vector opt.
184	zstms = zztms(ji,jj ) * wght(ji,jj,2,2) + zztms(ji-1,jj ) * wght(ji,jj,1,2) &
185	& + zztms(ji,jj-1) * wght(ji,jj,2,1) + zztms(ji-1,jj-1) * wght(ji,jj,1,1)
186	zusw = 1._wp / MAX( zstms, epsd )
187
188	zt11 = zztms(ji ,jj ) * zzfrld(ji ,jj )
189	zt12 = zztms(ji-1,jj ) * zzfrld(ji-1,jj )
190	zt21 = zztms(ji ,jj-1) * zzfrld(ji ,jj-1)
191	zt22 = zztms(ji-1,jj-1) * zzfrld(ji-1,jj-1)
192
193	! Leads area.
194	zfrld(ji,jj) = ( zt11 * wght(ji,jj,2,2) + zt12 * wght(ji,jj,1,2) &
195	& + zt21 * wght(ji,jj,2,1) + zt22 * wght(ji,jj,1,1) ) * zusw
196
197	! Mass and coriolis coeff. at I-point
198	zmass(ji,jj) = ( zmasst(ji,jj ) * wght(ji,jj,2,2) + zmasst(ji-1,jj ) * wght(ji,jj,1,2) &
199	& + zmasst(ji,jj-1) * wght(ji,jj,2,1) + zmasst(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw
200	zcorl(ji,jj) = zmass(ji,jj) &
201	& ( fcor(ji,jj ) wght(ji,jj,2,2) + fcor(ji-1,jj )*wght(ji,jj,1,2) &
202	& + fcor(ji,jj-1) * wght(ji,jj,2,1) + fcor(ji-1,jj-1)wght(ji,jj,1,1) ) zusw
203
204	! Wind stress.
205	! always provide stress at I-point
206	ztagnx = ( zi1(ji,jj ) * wght(ji,jj,2,2) + zi1(ji-1,jj ) * wght(ji,jj,1,2) &
207	& + zi1(ji,jj-1) * wght(ji,jj,2,1) + zi1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * utau_ice(ji,jj)
208	ztagny = ( zi2(ji,jj ) * wght(ji,jj,2,2) + zi2(ji-1,jj ) * wght(ji,jj,1,2) &
209	& + zi2(ji,jj-1) * wght(ji,jj,2,1) + zi2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * vtau_ice(ji,jj)
210
211	! Gradient of ice strength
212	zgphsx = ( alambd(ji,jj,2,2,2,1) - alambd(ji,jj,2,1,2,1) ) * zpresh(ji ,jj-1) &
213	& + ( alambd(ji,jj,2,2,2,2) - alambd(ji,jj,2,1,2,2) ) * zpresh(ji ,jj ) &
214	& - ( alambd(ji,jj,2,2,1,1) + alambd(ji,jj,2,1,1,1) ) * zpresh(ji-1,jj-1) &
215	& - ( alambd(ji,jj,2,2,1,2) + alambd(ji,jj,2,1,1,2) ) * zpresh(ji-1,jj )
216
217	zgphsy = - ( alambd(ji,jj,1,1,2,1) + alambd(ji,jj,1,2,2,1) ) * zpresh(ji ,jj-1) &
218	& - ( alambd(ji,jj,1,1,1,1) + alambd(ji,jj,1,2,1,1) ) * zpresh(ji-1,jj-1) &
219	& + ( alambd(ji,jj,1,1,2,2) - alambd(ji,jj,1,2,2,2) ) * zpresh(ji ,jj ) &
220	& + ( alambd(ji,jj,1,1,1,2) - alambd(ji,jj,1,2,1,2) ) * zpresh(ji-1,jj )
221
222	! Gradient of the sea surface height
223	zgsshx = ( (ssh_m(ji ,jj ) - ssh_m(ji-1,jj ))/e1u(ji-1,jj ) &
224	& + (ssh_m(ji ,jj-1) - ssh_m(ji-1,jj-1))/e1u(ji-1,jj-1) ) * 0.5_wp
225	zgsshy = ( (ssh_m(ji ,jj ) - ssh_m(ji ,jj-1))/e2v(ji ,jj-1) &
226	& + (ssh_m(ji-1,jj ) - ssh_m(ji-1,jj-1))/e2v(ji-1,jj-1) ) * 0.5_wp
227
228	! Computation of the velocity field taking into account the ice-ice interaction.
229	! Terms that are independent of the ice velocity field.
230	za1ct(ji,jj) = ztagnx - zmass(ji,jj) * grav * zgsshx - zgphsx
231	za2ct(ji,jj) = ztagny - zmass(ji,jj) * grav * zgsshy - zgphsy
232	END DO
233	END DO
234
235
236	! SOLUTION OF THE MOMENTUM EQUATION.
237	!------------------------------------------
238	! ! ==================== !
239	DO iter = 1 , 2 * nbiter ! loop over iter !
240	! ! ==================== !
241	zindu = MOD( iter , 2 )
242	zusdtp = ( zindu * 2._wp + ( 1._wp - zindu ) * 1._wp ) * REAL( nbiter ) / rdt_ice
243
244	! Computation of free drift field for free slip boundary conditions.
245
246	!CDIR NOVERRCHK
247	DO jj = k_j1, k_jpj-1
248	!CDIR NOVERRCHK
249	DO ji = 1, fs_jpim1
250	!- Rate of strain tensor.
251	zt11 = akappa(ji,jj,1,1) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) - zu_a(ji,jj ) - zu_a(ji ,jj+1) ) &
252	& + akappa(ji,jj,1,2) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) + zv_a(ji,jj ) + zv_a(ji ,jj+1) )
253	zt12 = - akappa(ji,jj,2,2) * ( zu_a(ji ,jj) + zu_a(ji+1,jj ) - zu_a(ji,jj+1) - zu_a(ji+1,jj+1) ) &
254	& - akappa(ji,jj,2,1) * ( zv_a(ji ,jj) + zv_a(ji+1,jj ) + zv_a(ji,jj+1) + zv_a(ji+1,jj+1) )
255	zt22 = - akappa(ji,jj,2,2) * ( zv_a(ji ,jj) + zv_a(ji+1,jj ) - zv_a(ji,jj+1) - zv_a(ji+1,jj+1) ) &
256	& + akappa(ji,jj,2,1) * ( zu_a(ji ,jj) + zu_a(ji+1,jj ) + zu_a(ji,jj+1) + zu_a(ji+1,jj+1) )
257	zt21 = akappa(ji,jj,1,1) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) - zv_a(ji,jj ) - zv_a(ji ,jj+1) ) &
258	& - akappa(ji,jj,1,2) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) + zu_a(ji,jj ) + zu_a(ji ,jj+1) )
259
260	!- Rate of strain tensor.
261	zdgp = zt11 + zt22
262	zdgi = zt12 + zt21
263	ztrace2 = zdgp * zdgp
264	zdeter = zt11 * zt22 - 0.25_wp * zdgi * zdgi
265
266	! Creep limit depends on the size of the grid.
267	zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4._wp * zdeter ) * usecc2 ), creepl)
268
269	!- Computation of viscosities.
270	zviszeta(ji,jj) = MAX( zpresh(ji,jj) / zdelta, etamn )
271	zviseta (ji,jj) = zviszeta(ji,jj) * usecc2
272	END DO
273	END DO
274
275	!- Determination of zc1u, zc2u, zc1v and zc2v.
276	DO jj = k_j1+1, k_jpj-1
277	DO ji = 2, fs_jpim1 ! NO vector opt.
278	!* zc1u , zc2v
279	zvis11 = 2._wp * zviseta (ji-1,jj-1) + dm
280	zvis12 = zviseta (ji-1,jj-1) + dm
281	zvis21 = zviseta (ji-1,jj-1)
282	zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
283	zdiag = zvis22 * ( akappa(ji-1,jj-1,1,1) + akappa(ji-1,jj-1,2,1) )
284	zs11_11 = zvis11 * akappa(ji-1,jj-1,1,1) + zdiag
285	zs12_11 = zvis12 * akappa(ji-1,jj-1,2,2) - zvis21 * akappa(ji-1,jj-1,1,2)
286	zs21_11 = -zvis12 * akappa(ji-1,jj-1,1,2) + zvis21 * akappa(ji-1,jj-1,2,2)
287	zs22_11 = zvis11 * akappa(ji-1,jj-1,2,1) + zdiag
288
289	zvis11 = 2._wp * zviseta (ji,jj-1) + dm
290	zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1)
291	zvis12 = zviseta (ji,jj-1) + dm
292	zvis21 = zviseta (ji,jj-1)
293	zdiag = zvis22 * ( -akappa(ji,jj-1,1,1) + akappa(ji,jj-1,2,1) )
294	zs11_21 = -zvis11 * akappa(ji,jj-1,1,1) + zdiag
295	zs12_21 = zvis12 * akappa(ji,jj-1,2,2) - zvis21 * akappa(ji,jj-1,1,2)
296	zs22_21 = zvis11 * akappa(ji,jj-1,2,1) + zdiag
297	zs21_21 = -zvis12 * akappa(ji,jj-1,1,2) + zvis21 * akappa(ji,jj-1,2,2)
298
299	zvis11 = 2._wp * zviseta (ji-1,jj) + dm
300	zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj)
301	zvis12 = zviseta (ji-1,jj) + dm
302	zvis21 = zviseta (ji-1,jj)
303	zdiag = zvis22 * ( akappa(ji-1,jj,1,1) + akappa(ji-1,jj,2,1) )
304	zs11_12 = zvis11 * akappa(ji-1,jj,1,1) + zdiag
305	zs12_12 = -zvis12 * akappa(ji-1,jj,2,2) - zvis21 * akappa(ji-1,jj,1,2)
306	zs22_12 = zvis11 * akappa(ji-1,jj,2,1) + zdiag
307	zs21_12 = -zvis12 * akappa(ji-1,jj,1,2) - zvis21 * akappa(ji-1,jj,2,2)
308
309	zvis11 = 2._wp * zviseta (ji,jj) + dm
310	zvis22 = zviszeta(ji,jj) - zviseta(ji,jj)
311	zvis12 = zviseta (ji,jj) + dm
312	zvis21 = zviseta (ji,jj)
313	zdiag = zvis22 * ( -akappa(ji,jj,1,1) + akappa(ji,jj,2,1) )
314	zs11_22 = -zvis11 * akappa(ji,jj,1,1) + zdiag
315	zs12_22 = -zvis12 * akappa(ji,jj,2,2) - zvis21 * akappa(ji,jj,1,2)
316	zs22_22 = zvis11 * akappa(ji,jj,2,1) + zdiag
317	zs21_22 = -zvis12 * akappa(ji,jj,1,2) - zvis21 * akappa(ji,jj,2,2)
318
319	zc1u(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22 &
320	& - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12 &
321	& - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11 &
322	& + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12 &
323	& + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21 &
324	& + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22 &
325	& - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21 &
326	& - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22
327
328	zc2u(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22 &
329	& - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12 &
330	& - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11 &
331	& + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12 &
332	& - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21 &
333	& - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22 &
334	& + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21 &
335	& + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22
336
337	!* zc1v , zc2v.
338	zvis11 = 2._wp * zviseta (ji-1,jj-1) + dm
339	zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
340	zvis12 = zviseta (ji-1,jj-1) + dm
341	zvis21 = zviseta (ji-1,jj-1)
342	zdiag = zvis22 * ( akappa(ji-1,jj-1,1,2) + akappa(ji-1,jj-1,2,2) )
343	zs11_11 = zvis11 * akappa(ji-1,jj-1,1,2) + zdiag
344	zs12_11 = -zvis12 * akappa(ji-1,jj-1,2,1) + zvis21 * akappa(ji-1,jj-1,1,1)
345	zs22_11 = zvis11 * akappa(ji-1,jj-1,2,2) + zdiag
346	zs21_11 = zvis12 * akappa(ji-1,jj-1,1,1) - zvis21 * akappa(ji-1,jj-1,2,1)
347
348	zvis11 = 2._wp * zviseta (ji,jj-1) + dm
349	zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1)
350	zvis12 = zviseta (ji,jj-1) + dm
351	zvis21 = zviseta (ji,jj-1)
352	zdiag = zvis22 * ( akappa(ji,jj-1,1,2) + akappa(ji,jj-1,2,2) )
353	zs11_21 = zvis11 * akappa(ji,jj-1,1,2) + zdiag
354	zs12_21 = -zvis12 * akappa(ji,jj-1,2,1) - zvis21 * akappa(ji,jj-1,1,1)
355	zs22_21 = zvis11 * akappa(ji,jj-1,2,2) + zdiag
356	zs21_21 = -zvis12 * akappa(ji,jj-1,1,1) - zvis21 * akappa(ji,jj-1,2,1)
357
358	zvis11 = 2._wp * zviseta (ji-1,jj) + dm
359	zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj)
360	zvis12 = zviseta (ji-1,jj) + dm
361	zvis21 = zviseta (ji-1,jj)
362	zdiag = zvis22 * ( akappa(ji-1,jj,1,2) - akappa(ji-1,jj,2,2) )
363	zs11_12 = zvis11 * akappa(ji-1,jj,1,2) + zdiag
364	zs12_12 = -zvis12 * akappa(ji-1,jj,2,1) + zvis21 * akappa(ji-1,jj,1,1)
365	zs22_12 = -zvis11 * akappa(ji-1,jj,2,2) + zdiag
366	zs21_12 = zvis12 * akappa(ji-1,jj,1,1) - zvis21 * akappa(ji-1,jj,2,1)
367
368	zvis11 = 2._wp * zviseta (ji,jj) + dm
369	zvis22 = zviszeta(ji,jj) - zviseta(ji,jj)
370	zvis12 = zviseta (ji,jj) + dm
371	zvis21 = zviseta (ji,jj)
372	zdiag = zvis22 * ( akappa(ji,jj,1,2) - akappa(ji,jj,2,2) )
373	zs11_22 = zvis11 * akappa(ji,jj,1,2) + zdiag
374	zs12_22 = -zvis12 * akappa(ji,jj,2,1) - zvis21 * akappa(ji,jj,1,1)
375	zs22_22 = -zvis11 * akappa(ji,jj,2,2) + zdiag
376	zs21_22 = -zvis12 * akappa(ji,jj,1,1) - zvis21 * akappa(ji,jj,2,1)
377
378	zc1v(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22 &
379	& - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12 &
380	& - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11 &
381	& + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12 &
382	& + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21 &
383	& + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22 &
384	& - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21 &
385	& - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22
386
387	zc2v(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22 &
388	& - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12 &
389	& - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11 &
390	& + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12 &
391	& - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21 &
392	& - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22 &
393	& + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21 &
394	& + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22
395	END DO
396	END DO
397
398	! GAUSS-SEIDEL method
399	! ! ================ !
400	iflag: DO jter = 1 , nbitdr ! Relaxation !
401	! ! ================ !
402	!CDIR NOVERRCHK
403	DO jj = k_j1+1, k_jpj-1
404	!CDIR NOVERRCHK
405	DO ji = 2, fs_jpim1 ! NO vector opt.
406	!
407	ze11 = akappa(ji,jj-1,1,1) * zu_a(ji+1,jj) + akappa(ji,jj-1,1,2) * zv_a(ji+1,jj)
408	ze12 = + akappa(ji,jj-1,2,2) * zu_a(ji+1,jj) - akappa(ji,jj-1,2,1) * zv_a(ji+1,jj)
409	ze22 = + akappa(ji,jj-1,2,2) * zv_a(ji+1,jj) + akappa(ji,jj-1,2,1) * zu_a(ji+1,jj)
410	ze21 = akappa(ji,jj-1,1,1) * zv_a(ji+1,jj) - akappa(ji,jj-1,1,2) * zu_a(ji+1,jj)
411	zvis11 = 2._wp * zviseta (ji,jj-1) + dm
412	zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1)
413	zvis12 = zviseta (ji,jj-1) + dm
414	zvis21 = zviseta (ji,jj-1)
415	zdiag = zvis22 * ( ze11 + ze22 )
416	zs11_21 = zvis11 * ze11 + zdiag
417	zs12_21 = zvis12 * ze12 + zvis21 * ze21
418	zs22_21 = zvis11 * ze22 + zdiag
419	zs21_21 = zvis12 * ze21 + zvis21 * ze12
420
421	ze11 = akappa(ji-1,jj,1,1) * ( zu_a(ji ,jj+1) - zu_a(ji-1,jj+1) ) &
422	& + akappa(ji-1,jj,1,2) * ( zv_a(ji ,jj+1) + zv_a(ji-1,jj+1) )
423	ze12 = + akappa(ji-1,jj,2,2) * ( zu_a(ji-1,jj+1) + zu_a(ji ,jj+1) ) &
424	& - akappa(ji-1,jj,2,1) * ( zv_a(ji-1,jj+1) + zv_a(ji ,jj+1) )
425	ze22 = + akappa(ji-1,jj,2,2) * ( zv_a(ji-1,jj+1) + zv_a(ji ,jj+1) ) &
426	& + akappa(ji-1,jj,2,1) * ( zu_a(ji-1,jj+1) + zu_a(ji ,jj+1) )
427	ze21 = akappa(ji-1,jj,1,1) * ( zv_a(ji ,jj+1) - zv_a(ji-1,jj+1) ) &
428	& - akappa(ji-1,jj,1,2) * ( zu_a(ji ,jj+1) + zu_a(ji-1,jj+1) )
429	zvis11 = 2._wp * zviseta (ji-1,jj) + dm
430	zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj)
431	zvis12 = zviseta (ji-1,jj) + dm
432	zvis21 = zviseta (ji-1,jj)
433	zdiag = zvis22 * ( ze11 + ze22 )
434	zs11_12 = zvis11 * ze11 + zdiag
435	zs12_12 = zvis12 * ze12 + zvis21 * ze21
436	zs22_12 = zvis11 * ze22 + zdiag
437	zs21_12 = zvis12 * ze21 + zvis21 * ze12
438
439	ze11 = akappa(ji,jj,1,1) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) - zu_a(ji ,jj+1) ) &
440	& + akappa(ji,jj,1,2) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) + zv_a(ji ,jj+1) )
441	ze12 = - akappa(ji,jj,2,2) * ( zu_a(ji+1,jj) - zu_a(ji ,jj+1) - zu_a(ji+1,jj+1) ) &
442	& - akappa(ji,jj,2,1) * ( zv_a(ji+1,jj) + zv_a(ji ,jj+1) + zv_a(ji+1,jj+1) )
443	ze22 = - akappa(ji,jj,2,2) * ( zv_a(ji+1,jj) - zv_a(ji ,jj+1) - zv_a(ji+1,jj+1) ) &
444	& + akappa(ji,jj,2,1) * ( zu_a(ji+1,jj) + zu_a(ji ,jj+1) + zu_a(ji+1,jj+1) )
445	ze21 = akappa(ji,jj,1,1) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) - zv_a(ji ,jj+1) ) &
446	& - akappa(ji,jj,1,2) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) + zu_a(ji ,jj+1) )
447	zvis11 = 2._wp * zviseta (ji,jj) + dm
448	zvis22 = zviszeta(ji,jj) - zviseta(ji,jj)
449	zvis12 = zviseta (ji,jj) + dm
450	zvis21 = zviseta (ji,jj)
451	zdiag = zvis22 * ( ze11 + ze22 )
452	zs11_22 = zvis11 * ze11 + zdiag
453	zs12_22 = zvis12 * ze12 + zvis21 * ze21
454	zs22_22 = zvis11 * ze22 + zdiag
455	zs21_22 = zvis12 * ze21 + zvis21 * ze12
456
457	! 2nd part
458	ze11 = akappa(ji-1,jj-1,1,1) * ( zu_a(ji ,jj-1) - zu_a(ji-1,jj-1) - zu_a(ji-1,jj) ) &
459	& + akappa(ji-1,jj-1,1,2) * ( zv_a(ji ,jj-1) + zv_a(ji-1,jj-1) + zv_a(ji-1,jj) )
460	ze12 = - akappa(ji-1,jj-1,2,2) * ( zu_a(ji-1,jj-1) + zu_a(ji ,jj-1) - zu_a(ji-1,jj) ) &
461	& - akappa(ji-1,jj-1,2,1) * ( zv_a(ji-1,jj-1) + zv_a(ji ,jj-1) + zv_a(ji-1,jj) )
462	ze22 = - akappa(ji-1,jj-1,2,2) * ( zv_a(ji-1,jj-1) + zv_a(ji ,jj-1) - zv_a(ji-1,jj) ) &
463	& + akappa(ji-1,jj-1,2,1) * ( zu_a(ji-1,jj-1) + zu_a(ji ,jj-1) + zu_a(ji-1,jj) )
464	ze21 = akappa(ji-1,jj-1,1,1) * ( zv_a(ji ,jj-1) - zv_a(ji-1,jj-1) - zv_a(ji-1,jj) ) &
465	& - akappa(ji-1,jj-1,1,2) * ( zu_a(ji ,jj-1) + zu_a(ji-1,jj-1) + zu_a(ji-1,jj) )
466	zvis11 = 2._wp * zviseta (ji-1,jj-1) + dm
467	zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
468	zvis12 = zviseta (ji-1,jj-1) + dm
469	zvis21 = zviseta (ji-1,jj-1)
470	zdiag = zvis22 * ( ze11 + ze22 )
471	zs11_11 = zvis11 * ze11 + zdiag
472	zs12_11 = zvis12 * ze12 + zvis21 * ze21
473	zs22_11 = zvis11 * ze22 + zdiag
474	zs21_11 = zvis12 * ze21 + zvis21 * ze12
475
476	ze11 = akappa(ji,jj-1,1,1) * ( zu_a(ji+1,jj-1) - zu_a(ji ,jj-1) ) &
477	& + akappa(ji,jj-1,1,2) * ( zv_a(ji+1,jj-1) + zv_a(ji ,jj-1) )
478	ze12 = - akappa(ji,jj-1,2,2) * ( zu_a(ji ,jj-1) + zu_a(ji+1,jj-1) ) &
479	& - akappa(ji,jj-1,2,1) * ( zv_a(ji ,jj-1) + zv_a(ji+1,jj-1) )
480	ze22 = - akappa(ji,jj-1,2,2) * ( zv_a(ji ,jj-1) + zv_a(ji+1,jj-1) ) &
481	& + akappa(ji,jj-1,2,1) * ( zu_a(ji ,jj-1) + zu_a(ji+1,jj-1) )
482	ze21 = akappa(ji,jj-1,1,1) * ( zv_a(ji+1,jj-1) - zv_a(ji ,jj-1) ) &
483	& - akappa(ji,jj-1,1,2) * ( zu_a(ji+1,jj-1) + zu_a(ji ,jj-1) )
484	zvis11 = 2._wp * zviseta (ji,jj-1) + dm
485	zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1)
486	zvis12 = zviseta (ji,jj-1) + dm
487	zvis21 = zviseta (ji,jj-1)
488	zdiag = zvis22 * ( ze11 + ze22 )
489	zs11_21 = zs11_21 + zvis11 * ze11 + zdiag
490	zs12_21 = zs12_21 + zvis12 * ze12 + zvis21 * ze21
491	zs22_21 = zs22_21 + zvis11 * ze22 + zdiag
492	zs21_21 = zs21_21 + zvis12 * ze21 + zvis21 * ze12
493
494	ze11 = - akappa(ji-1,jj,1,1) * zu_a(ji-1,jj) + akappa(ji-1,jj,1,2) * zv_a(ji-1,jj)
495	ze12 = - akappa(ji-1,jj,2,2) * zu_a(ji-1,jj) - akappa(ji-1,jj,2,1) * zv_a(ji-1,jj)
496	ze22 = - akappa(ji-1,jj,2,2) * zv_a(ji-1,jj) + akappa(ji-1,jj,2,1) * zu_a(ji-1,jj)
497	ze21 = - akappa(ji-1,jj,1,1) * zv_a(ji-1,jj) - akappa(ji-1,jj,1,2) * zu_a(ji-1,jj)
498	zvis11 = 2._wp * zviseta (ji-1,jj) + dm
499	zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj)
500	zvis12 = zviseta (ji-1,jj) + dm
501	zvis21 = zviseta (ji-1,jj)
502	zdiag = zvis22 * ( ze11 + ze22 )
503	zs11_12 = zs11_12 + zvis11 * ze11 + zdiag
504	zs12_12 = zs12_12 + zvis12 * ze12 + zvis21 * ze21
505	zs22_12 = zs22_12 + zvis11 * ze22 + zdiag
506	zs21_12 = zs21_12 + zvis12 * ze21 + zvis21 * ze12
507
508	zd1 = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22 &
509	& - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12 &
510	& - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11 &
511	& + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12 &
512	& + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21 &
513	& + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22 &
514	& - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21 &
515	& - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22
516
517	zd2 = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22 &
518	& - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12 &
519	& - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11 &
520	& + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12 &
521	& - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21 &
522	& - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22 &
523	& + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21 &
524	& + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22
525
526	zur = zu_a(ji,jj) - u_oce(ji,jj)
527	zvr = zv_a(ji,jj) - v_oce(ji,jj)
528	!!!!
529	zmod = SQRT( zurzur + zvrzvr ) * ( 1._wp - zfrld(ji,jj) )
530	za = rhoco * zmod
531	!!!!
532	!!gm chg resul za = rhoco * SQRT( zurzur + zvrzvr ) * ( 1._wp - zfrld(ji,jj) )
533	zac = za * cangvg
534	zmpzas = alpha * zcorl(ji,jj) + za * zsang(ji,jj)
535	zmassdt = zusdtp * zmass(ji,jj)
536	zcorlal = ( 1._wp - alpha ) * zcorl(ji,jj)
537
538	za1 = zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj) &
539	& + za * ( cangvg * u_oce(ji,jj) - zsang(ji,jj) * v_oce(ji,jj) )
540	za2 = zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj) &
541	& + za * ( cangvg * v_oce(ji,jj) + zsang(ji,jj) * u_oce(ji,jj) )
542	zb1 = zmassdt + zac - zc1u(ji,jj)
543	zb2 = zmpzas - zc2u(ji,jj)
544	zc1 = zmpzas + zc1v(ji,jj)
545	zc2 = zmassdt + zac - zc2v(ji,jj)
546	zdeter = zc1 * zb2 + zc2 * zb1
547	zden = SIGN( rone, zdeter) / MAX( epsd , ABS( zdeter ) )
548	zunw = ( ( za1 + zd1 ) * zc2 + ( za2 + zd2 ) * zc1 ) * zden
549	zvnw = ( ( za2 + zd2 ) * zb1 - ( za1 + zd1 ) * zb2 ) * zden
550	zmask = ( 1._wp - MAX( rzero, SIGN( rone , 1._wp - zmass(ji,jj) ) ) ) * tmu(ji,jj)
551
552	zu_n(ji,jj) = ( zu_a(ji,jj) + om * ( zunw - zu_a(ji,jj) ) * tmu(ji,jj) ) * zmask
553	zv_n(ji,jj) = ( zv_a(ji,jj) + om * ( zvnw - zv_a(ji,jj) ) * tmu(ji,jj) ) * zmask
554	END DO
555	END DO
556
557	CALL lbc_lnk( zu_n(:,1:jpj), 'I', -1. )
558	CALL lbc_lnk( zv_n(:,1:jpj), 'I', -1. )
559
560	#if defined key_agrif
561	! copy the boundary value from u_ice_nst and v_ice_nst to u_ice and v_ice
562	! before next interations
563	CALL Agrif_dyn_lim2_copy(zu_n,zv_n)
564	#endif
565
566	! Test of Convergence
567	DO jj = k_j1+1 , k_jpj-1
568	zresr(:,jj) = MAX( ABS( zu_a(:,jj) - zu_n(:,jj) ) , ABS( zv_a(:,jj) - zv_n(:,jj) ) )
569	END DO
570	zresm = MAXVAL( zresr(1:jpi,k_j1+1:k_jpj-1) )
571	!!!! this should be faster on scalar processor
572	! zresm = MAXVAL( MAX( ABS( zu_a(1:jpi,k_j1+1:k_jpj-1) - zu_n(1:jpi,k_j1+1:k_jpj-1) ), &
573	! & ABS( zv_a(1:jpi,k_j1+1:k_jpj-1) - zv_n(1:jpi,k_j1+1:k_jpj-1) ) ) )
574	!!!!
575	IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain
576
577	DO jj = k_j1, k_jpj
578	zu_a(:,jj) = zu_n(:,jj)
579	zv_a(:,jj) = zv_n(:,jj)
580	END DO
581
582	IF( zresm <= resl ) EXIT iflag
583
584	! ! ================ !
585	END DO iflag ! end Relaxation !
586	! ! ================ !
587
588	IF( zindu == 0 ) THEN ! even iteration
589	DO jj = k_j1 , k_jpj-1
590	zu0(:,jj) = zu_a(:,jj)
591	zv0(:,jj) = zv_a(:,jj)
592	END DO
593	ENDIF
594	! ! ==================== !
595	END DO ! end loop over iter !
596	! ! ==================== !
597
598	u_ice(:,:) = zu_a(:,1:jpj)
599	v_ice(:,:) = zv_a(:,1:jpj)
600
601	IF(ln_ctl) THEN
602	WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, jter
603	CALL prt_ctl_info(charout)
604	CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
605	ENDIF
606
607	END SUBROUTINE lim_rhg_2
608
609	#else
610	!!----------------------------------------------------------------------
611	!! Default option Dummy module NO VP & LIM-2 sea-ice model
612	!!----------------------------------------------------------------------
613	CONTAINS
614	SUBROUTINE lim_rhg_2( k1 , k2 ) ! Dummy routine
615	WRITE(,) 'lim_rhg_2: You should not have seen this print! error?', k1, k2
616	END SUBROUTINE lim_rhg_2
617	#endif
618
619	!!==============================================================================
620	END MODULE limrhg_2

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

Download in other formats: