Context Navigation

limrhg_2.F90 @ 7646

Last change on this file since 7646 was 7646, checked in by timgraham, 7 years ago

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

Property svn:keywords set to Id

File size: 35.7 KB

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90 @ 7646

Download in other formats: