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/UKMO/dev_r5518_area_calc/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/UKMO/dev_r5518_area_calc/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90 @ 6740

Last change on this file since 6740 was 6740, checked in by jrioual, 8 years ago

Remove svn keywords

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