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_v3.2.F90 on Ticket #587 – Attachment – NEMO

Ticket #587: limrhg_2_v3.2.F90

File limrhg_2_v3.2.F90, 33.2 KB (added by gm, 14 years ago)

limrhg_2.F90 for the v3.2 with the use of ssh_m

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