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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limrhg.F90 in trunk/NEMO/LIM_SRC – NEMO

source: trunk/NEMO/LIM_SRC/limrhg.F90 @ 700

Last change on this file since 700 was 700, checked in by smasson, 17 years ago

sbc(1): bugfix on wind stress computation #2
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 33.8 KB
Line 
1MODULE limrhg
2   !!======================================================================
3   !!                     ***  MODULE  limrhg  ***
4   !!   Ice rheology :  performs sea ice rheology
5   !!======================================================================
6   !! History :   0.0  !  93-12  (M.A. Morales Maqueda.)  Original code
7   !!             1.0  !  94-12  (H. Goosse)
8   !!             2.0  !  03-12  (C. Ethe, G. Madec)  F90, mpp
9   !!             3.0  !  07-06  (S. Masson, G. Madec)  ice/atmosphere stress
10   !!                             provided at I-point in forced and coupled mode
11   !!----------------------------------------------------------------------
12#if defined key_ice_lim
13   !!----------------------------------------------------------------------
14   !!   'key_ice_lim'                                     LIM sea-ice model
15   !!----------------------------------------------------------------------
16   !!   lim_rhg   : computes ice velocities
17   !!----------------------------------------------------------------------
18   !! * Modules used
19   USE phycst
20   USE par_oce
21   USE ice_oce         ! ice variables
22   USE dom_ice
23   USE ice
24   USE lbclnk
25   USE lib_mpp
26   USE in_out_manager  ! I/O manager
27   USE prtctl          ! Print control
28
29   IMPLICIT NONE
30   PRIVATE
31
32   !! * Routine accessibility
33   PUBLIC lim_rhg  ! routine called by lim_dyn
34
35   !! * Module variables
36   REAL(wp)  ::           &  ! constant values
37      rzero   = 0.e0   ,  &
38      rone    = 1.e0
39   !!----------------------------------------------------------------------
40   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)
41   !! $Header$
42   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
43   !!----------------------------------------------------------------------
44
45CONTAINS
46
47   SUBROUTINE lim_rhg( k_j1, k_jpj )
48      !!-------------------------------------------------------------------
49      !!                 ***  SUBROUTINR lim_rhg  ***
50      !!
51      !! ** purpose :   determines the velocity field of sea ice by using
52      !!  atmospheric (wind stress) and oceanic (water stress and surface
53      !!  tilt) forcings. Ice-ice interaction is described by a non-linear
54      !!  viscous-plastic law including shear strength and a bulk rheology.
55      !!
56      !! ** Action  : - compute u_ice, v_ice the sea-ice velocity
57      !!
58      !!-------------------------------------------------------------------
59      INTEGER, INTENT(in) ::   k_j1 ,  &  ! southern j-index for ice computation
60         &                     k_jpj      ! northern j-index for ice computation
61
62      INTEGER ::   ji, jj              ! dummy loop indices
63
64      INTEGER  :: &
65         iim1, ijm1, iip1 , ijp1   , & ! temporary integers
66         iter, jter                    !    "          "
67
68      CHARACTER (len=50) :: charout
69
70      REAL(wp) :: &
71         ze11  , ze12  , ze22  , ze21  ,   &  ! temporary scalars
72         zt11  , zt12  , zt21  , zt22  ,   &  !    "         "
73         zvis11, zvis21, zvis12, zvis22,   &  !    "         "
74         zgphsx, ztagnx, zusw  ,           &  !    "         "
75         zgphsy, ztagny                       !    "         "
76      REAL(wp) :: &
77         zresm, zunw, zvnw, zur, zvr, zmod, za, zac, &
78         zmpzas, zstms, zindu, zindu1, zusdtp, zmassdt, zcorlal,  &
79         ztrace2, zdeter, zdelta, zsang, zmask, zdgp, zdgi, zdiag, zic
80      REAL(wp),DIMENSION(jpi,jpj) ::   &
81         zpresh, zfrld, zmass, zcorl,     &
82         zu0, zv0, zviszeta, zviseta,     &
83         zc1u, zc1v, zc2u, zc2v, za1ct, za2ct, za1, za2, zb1, zb2,  &
84         zc1, zc2, zd1, zd2, zden, zu_ice, zv_ice, zresr
85      REAL(wp),DIMENSION(jpi,jpj,2,2) :: &
86         zs11, zs12, zs22, zs21
87      !!-------------------------------------------------------------------
88     
89      !  Store initial velocities
90      !  ------------------------
91      zu0(:,:) = u_ice(:,:)
92      zv0(:,:) = v_ice(:,:)
93
94      ! masked Ice mass, ice strength, Ice-cover fraction and Lead faction  at T-point
95      !-------------------------------------------------------------------
96
97      DO jj = k_j1 , k_jpj-1
98         DO ji = 1 , jpi
99            za1(ji,jj)    = tms(ji,jj) * ( rhosn * hsnm(ji,jj) + rhoic * hicm(ji,jj) )
100            zpresh(ji,jj) = tms(ji,jj) *  pstarh * hicm(ji,jj) * EXP( -c_rhg * frld(ji,jj) )
101
102            zb1(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) )   
103            zb2(ji,jj)    = tms(ji,jj) *         frld(ji,jj)     
104         END DO
105      END DO
106
107      !  Wind stress, coriolis, mass and Gradient of ice strenght at I-point
108      !---------------------------------------------------------------------------
109         
110      DO jj = k_j1+1, k_jpj-1
111         DO ji = 2, jpi
112            zstms = tms(ji,jj  ) * wght(ji,jj,2,2) + tms(ji-1,jj  ) * wght(ji,jj,1,2)   &
113               &  + tms(ji,jj-1) * wght(ji,jj,2,1) + tms(ji-1,jj-1) * wght(ji,jj,1,1)
114            zusw  = 1.0 / MAX( zstms, epsd )
115
116            ! Leads area at I-point
117            zfrld(ji,jj) = ( zb2(ji,jj  ) * wght(ji,jj,2,2) + zb2(ji-1,jj  ) * wght(ji,jj,1,2)   &
118               &           + zb2(ji,jj-1) * wght(ji,jj,2,1) + zb2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw
119
120            ! Ice cover area at I-point
121            zic          = ( zb1(ji,jj  ) * wght(ji,jj,2,2) + zb1(ji-1,jj  ) * wght(ji,jj,1,2)   &
122               &           + zb1(ji,jj-1) * wght(ji,jj,2,1) + zb1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw
123
124            ! Wind stress.
125            ztagnx = zic * gtaux(ji,jj)
126            ztagny = zic * gtauy(ji,jj)
127
128            ! Mass and coriolis coeff.
129            zmass(ji,jj) = ( za1(ji,jj  ) * wght(ji,jj,2,2) + za1(ji-1,jj  ) * wght(ji,jj,1,2)   &
130               &           + za1(ji,jj-1) * wght(ji,jj,2,1) + za1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw
131            zcorl(ji,jj) = zmass(ji,jj) * fcor(ji,jj)
132
133            ! Gradient of ice strength
134            zgphsx =   ( alambd(ji,jj,2,2,2,1) - alambd(ji,jj,2,1,2,1) ) * zpresh(ji  ,jj-1)   &
135               &     + ( alambd(ji,jj,2,2,2,2) - alambd(ji,jj,2,1,2,2) ) * zpresh(ji  ,jj  )   &
136               &     - ( alambd(ji,jj,2,2,1,1) + alambd(ji,jj,2,1,1,1) ) * zpresh(ji-1,jj-1)   &
137               &     - ( alambd(ji,jj,2,2,1,2) + alambd(ji,jj,2,1,1,2) ) * zpresh(ji-1,jj  )
138
139            zgphsy = - ( alambd(ji,jj,1,1,2,1) + alambd(ji,jj,1,2,2,1) ) * zpresh(ji  ,jj-1)   &
140               &     - ( alambd(ji,jj,1,1,1,1) + alambd(ji,jj,1,2,1,1) ) * zpresh(ji-1,jj-1)   &
141               &     + ( alambd(ji,jj,1,1,2,2) - alambd(ji,jj,1,2,2,2) ) * zpresh(ji  ,jj  )   &
142               &     + ( alambd(ji,jj,1,1,1,2) - alambd(ji,jj,1,2,1,2) ) * zpresh(ji-1,jj  )
143
144            ! Computation of the velocity field taking into account the ice-ice interaction.                                 
145            ! Terms that are independent of the velocity field.
146            za1ct(ji,jj) = ztagnx - zcorl(ji,jj) * v_io(ji,jj) - zgphsx
147            za2ct(ji,jj) = ztagny + zcorl(ji,jj) * u_io(ji,jj) - zgphsy
148         END DO
149      END DO
150
151!! inutile!!
152!!??    CALL lbc_lnk( za1ct, 'I', -1. )
153!!??    CALL lbc_lnk( za2ct, 'I', -1. )
154
155
156      ! SOLUTION OF THE MOMENTUM EQUATION.
157      !------------------------------------------
158      !                                                   ! ==================== !
159      DO iter = 1 , 2 * nbiter                            !    loop over iter    !
160         !                                                ! ==================== !       
161         zindu = MOD( iter , 2 )
162         zusdtp = ( zindu * 2.0 + ( 1.0 - zindu ) * 1.0 )  * REAL( nbiter ) / rdt_ice
163
164         ! Computation of free drift field for free slip boundary conditions.
165
166           DO jj = k_j1, k_jpj-1
167              DO ji = 1, jpim1
168                 !- Rate of strain tensor.
169                 zt11 =   akappa(ji,jj,1,1) * ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) - u_ice(ji,jj  ) - u_ice(ji  ,jj+1) )  &
170                    &   + akappa(ji,jj,1,2) * ( v_ice(ji+1,jj) + v_ice(ji+1,jj+1) + v_ice(ji,jj  ) + v_ice(ji  ,jj+1) )
171                 zt12 = - akappa(ji,jj,2,2) * ( u_ice(ji  ,jj) + u_ice(ji+1,jj  ) - u_ice(ji,jj+1) - u_ice(ji+1,jj+1) )  &
172                    &   - akappa(ji,jj,2,1) * ( v_ice(ji  ,jj) + v_ice(ji+1,jj  ) + v_ice(ji,jj+1) + v_ice(ji+1,jj+1) )
173                 zt22 = - akappa(ji,jj,2,2) * ( v_ice(ji  ,jj) + v_ice(ji+1,jj  ) - v_ice(ji,jj+1) - v_ice(ji+1,jj+1) )  &
174                    &   + akappa(ji,jj,2,1) * ( u_ice(ji  ,jj) + u_ice(ji+1,jj  ) + u_ice(ji,jj+1) + u_ice(ji+1,jj+1) )
175                 zt21 =   akappa(ji,jj,1,1) * ( v_ice(ji+1,jj) + v_ice(ji+1,jj+1) - v_ice(ji,jj  ) - v_ice(ji  ,jj+1) )  &
176                    &   - akappa(ji,jj,1,2) * ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) + u_ice(ji,jj  ) + u_ice(ji  ,jj+1) )
177
178                 !- Rate of strain tensor.
179                 zdgp = zt11 + zt22
180                 zdgi = zt12 + zt21
181                 ztrace2 = zdgp * zdgp 
182                 zdeter  = zt11 * zt22 - 0.25 * zdgi * zdgi
183
184                 !  Creep limit depends on the size of the grid.
185                 zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4.0 * zdeter ) * usecc2),  creepl)
186
187                 !-  Computation of viscosities.
188                 zviszeta(ji,jj) = MAX( zpresh(ji,jj) / zdelta, etamn )
189                 zviseta (ji,jj) = zviszeta(ji,jj) * usecc2
190              END DO
191           END DO
192!!??       CALL lbc_lnk( zviszeta, 'I', -1. )  ! or T point???   semble reellement inutile
193!!??       CALL lbc_lnk( zviseta , 'I', -1. )
194
195
196           !-  Determination of zc1u, zc2u, zc1v and zc2v.
197           DO jj = k_j1+1, k_jpj-1
198              DO ji = 2, jpim1
199                 ze11   =  akappa(ji-1,jj-1,1,1)
200                 ze12   = +akappa(ji-1,jj-1,2,2)
201                 ze22   =  akappa(ji-1,jj-1,2,1)
202                 ze21   = -akappa(ji-1,jj-1,1,2)
203                 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm
204                 zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
205                 zvis12 =       zviseta (ji-1,jj-1) + dm
206                 zvis21 =       zviseta (ji-1,jj-1)
207
208                 zdiag = zvis22 * ( ze11 + ze22 )
209                 zs11(ji,jj,1,1) =  zvis11 * ze11 + zdiag
210                 zs12(ji,jj,1,1) =  zvis12 * ze12 + zvis21 * ze21
211                 zs22(ji,jj,1,1) =  zvis11 * ze22 + zdiag
212                 zs21(ji,jj,1,1) =  zvis12 * ze21 + zvis21 * ze12
213
214                 ze11   = -akappa(ji,jj-1,1,1)
215                 ze12   = +akappa(ji,jj-1,2,2)
216                 ze22   =  akappa(ji,jj-1,2,1)
217                 ze21   = -akappa(ji,jj-1,1,2)
218                 zvis11 = 2.0 * zviseta (ji,jj-1) + dm
219                 zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1)
220                 zvis12 =       zviseta (ji,jj-1) + dm
221                 zvis21 =       zviseta (ji,jj-1)
222
223                 zdiag = zvis22 * ( ze11 + ze22 )
224                 zs11(ji,jj,2,1) =  zvis11 * ze11 + zdiag
225                 zs12(ji,jj,2,1) =  zvis12 * ze12 + zvis21 * ze21
226                 zs22(ji,jj,2,1) =  zvis11 * ze22 + zdiag
227                 zs21(ji,jj,2,1) =  zvis12 * ze21 + zvis21 * ze12
228
229                 ze11   =  akappa(ji-1,jj,1,1)
230                 ze12   = -akappa(ji-1,jj,2,2)
231                 ze22   =  akappa(ji-1,jj,2,1)
232                 ze21   = -akappa(ji-1,jj,1,2)
233                 zvis11 = 2.0 * zviseta (ji-1,jj) + dm
234                 zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj)
235                 zvis12 =       zviseta (ji-1,jj) + dm
236                 zvis21 =       zviseta (ji-1,jj)
237
238                 zdiag = zvis22 * ( ze11 + ze22 ) 
239                 zs11(ji,jj,1,2) =  zvis11 * ze11 + zdiag
240                 zs12(ji,jj,1,2) =  zvis12 * ze12 + zvis21 * ze21
241                 zs22(ji,jj,1,2) =  zvis11 * ze22 + zdiag
242                 zs21(ji,jj,1,2) =  zvis12 * ze21 + zvis21 * ze12
243
244                 ze11   = -akappa(ji,jj,1,1)
245                 ze12   = -akappa(ji,jj,2,2)
246                 ze22   =  akappa(ji,jj,2,1)
247                 ze21   = -akappa(ji,jj,1,2)
248                 zvis11 = 2.0 * zviseta (ji,jj) + dm
249                 zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj)
250                 zvis12 =       zviseta (ji,jj) + dm
251                 zvis21 =       zviseta (ji,jj)
252
253                 zdiag = zvis22 * ( ze11 + ze22 )
254                 zs11(ji,jj,2,2) =  zvis11 * ze11 + zdiag
255                 zs12(ji,jj,2,2) =  zvis12 * ze12 + zvis21 * ze21
256                 zs22(ji,jj,2,2) =  zvis11 * ze22 + zdiag 
257                 zs21(ji,jj,2,2) =  zvis12 * ze21 + zvis21 * ze12
258              END DO
259           END DO
260
261           DO jj = k_j1+1, k_jpj-1
262              DO ji = 2, jpim1
263                 zc1u(ji,jj) =   &
264                    + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2)   &
265                    - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2)   &
266                    - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1)   &
267                    + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2)   &
268                    + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1)   &
269                    + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2)   &
270                    - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1)   &
271                    - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2)
272                 
273                 zc2u(ji,jj) =   &
274                    + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2)   &
275                    - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2)   &
276                    - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1)   &
277                    + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2)   &
278                    - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1)   &
279                    - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2)   &
280                    + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1)   &
281                    + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2)
282             END DO
283           END DO
284
285           DO jj = k_j1+1, k_jpj-1
286              DO ji = 2, jpim1
287                 !  zc1v , zc2v.
288                 ze11   =  akappa(ji-1,jj-1,1,2)
289                 ze12   = -akappa(ji-1,jj-1,2,1)
290                 ze22   = +akappa(ji-1,jj-1,2,2)
291                 ze21   =  akappa(ji-1,jj-1,1,1)
292                 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm
293                 zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
294                 zvis12 =       zviseta (ji-1,jj-1) + dm
295                 zvis21 =       zviseta (ji-1,jj-1)
296
297                 zdiag = zvis22 * ( ze11 + ze22 )
298                 zs11(ji,jj,1,1) =  zvis11 * ze11 + zdiag
299                 zs12(ji,jj,1,1) =  zvis12 * ze12 + zvis21 * ze21
300                 zs22(ji,jj,1,1) =  zvis11 * ze22 + zdiag
301                 zs21(ji,jj,1,1) =  zvis12 * ze21 + zvis21 * ze12
302 
303                 ze11   =  akappa(ji,jj-1,1,2)
304                 ze12   = -akappa(ji,jj-1,2,1)
305                 ze22   = +akappa(ji,jj-1,2,2)
306                 ze21   = -akappa(ji,jj-1,1,1)
307                 zvis11 = 2.0 * zviseta (ji,jj-1) + dm
308                 zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1)
309                 zvis12 =       zviseta (ji,jj-1) + dm
310                 zvis21 =       zviseta (ji,jj-1)
311
312                 zdiag = zvis22 * ( ze11 + ze22 )
313                 zs11(ji,jj,2,1) =  zvis11 * ze11 + zdiag
314                 zs12(ji,jj,2,1) =  zvis12 * ze12 + zvis21 * ze21
315                 zs22(ji,jj,2,1) =  zvis11 * ze22 + zdiag
316                 zs21(ji,jj,2,1) =  zvis12 * ze21 + zvis21 * ze12
317
318                 ze11   =  akappa(ji-1,jj,1,2)
319                 ze12   = -akappa(ji-1,jj,2,1)
320                 ze22   = -akappa(ji-1,jj,2,2)
321                 ze21   =  akappa(ji-1,jj,1,1)
322                 zvis11 = 2.0 * zviseta (ji-1,jj) + dm
323                 zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj)
324                 zvis12 =       zviseta (ji-1,jj) + dm
325                 zvis21 =       zviseta (ji-1,jj)
326
327                 zdiag = zvis22 * ( ze11 + ze22 )
328                 zs11(ji,jj,1,2) =  zvis11 * ze11 + zdiag
329                 zs12(ji,jj,1,2) =  zvis12 * ze12 + zvis21 * ze21
330                 zs22(ji,jj,1,2) =  zvis11 * ze22 + zdiag
331                 zs21(ji,jj,1,2) =  zvis12 * ze21 + zvis21 * ze12
332
333                 ze11   =  akappa(ji,jj,1,2)
334                 ze12   = -akappa(ji,jj,2,1)
335                 ze22   = -akappa(ji,jj,2,2)
336                 ze21   = -akappa(ji,jj,1,1)
337                 zvis11 = 2.0 * zviseta (ji,jj) + dm
338                 zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj)
339                 zvis12 =       zviseta (ji,jj) + dm
340                 zvis21 =       zviseta (ji,jj)
341
342                 zdiag = zvis22 * ( ze11 + ze22 )
343                 zs11(ji,jj,2,2) =  zvis11 * ze11 + zdiag
344                 zs12(ji,jj,2,2) =  zvis12 * ze12 + zvis21 * ze21
345                 zs22(ji,jj,2,2) =  zvis11 * ze22 + zdiag
346                 zs21(ji,jj,2,2) =  zvis12 * ze21 + zvis21 * ze12
347
348              END DO
349           END DO
350
351           DO jj = k_j1+1, k_jpj-1
352              DO ji = 2, jpim1
353                 zc1v(ji,jj) =   &
354                    + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2)   &
355                    - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2)   &
356                    - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1)   &
357                    + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2)   &
358                    + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1)   &
359                    + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2)   &
360                    - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1)   &
361                    - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2)
362                 zc2v(ji,jj) =   &
363                    + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2)   &
364                    - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2)   &
365                    - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1)   &
366                    + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2)   &
367                    - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1)   &
368                    - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2)   &
369                    + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1)   &
370                    + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2)
371              END DO
372           END DO
373
374         ! Relaxation.
375           
376iflag:   DO jter = 1 , nbitdr
377
378            !  Store previous drift field.   
379            DO jj = k_j1, k_jpj-1
380               zu_ice(:,jj) = u_ice(:,jj)
381               zv_ice(:,jj) = v_ice(:,jj)
382            END DO
383
384            DO jj = k_j1+1, k_jpj-1
385               zsang  = SIGN( 1.e0, fcor(1,jj) ) * sangvg   ! only the sinus changes its sign with the hemisphere
386               DO ji = 2, jpim1
387                 zur     = u_ice(ji,jj) - u_io(ji,jj)
388                 zvr     = v_ice(ji,jj) - v_io(ji,jj)
389                 zmod    = SQRT( zur * zur + zvr * zvr) * ( 1.0 - zfrld(ji,jj) )
390                 za      = rhoco * zmod
391                 zac     = za * cangvg
392                  zmpzas  = alpha * zcorl(ji,jj) + za * zsang
393                  zmassdt = zusdtp * zmass(ji,jj)
394                  zcorlal = ( 1.0 - alpha ) * zcorl(ji,jj)
395
396                  za1(ji,jj) =  zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj)   &
397                     &        + za * ( cangvg * u_io(ji,jj) - zsang * v_io(ji,jj) )
398
399                  za2(ji,jj) =  zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj)   &
400                     &        + za * ( cangvg * v_io(ji,jj) + zsang * u_io(ji,jj) )
401
402                  zb1(ji,jj)  = zmassdt + zac - zc1u(ji,jj)
403                  zb2(ji,jj)  = zmpzas  - zc2u(ji,jj)
404                  zc1(ji,jj)  = zmpzas  + zc1v(ji,jj)
405                  zc2(ji,jj)  = zmassdt + zac  - zc2v(ji,jj) 
406                  zdeter      = zc1(ji,jj) * zb2(ji,jj) + zc2(ji,jj) * zb1(ji,jj)
407                  zden(ji,jj) = SIGN( rone, zdeter) / MAX( epsd , ABS( zdeter ) )
408               END DO
409            END DO
410
411            ! The computation of ice interaction term is splitted into two parts
412            !-------------------------------------------------------------------------
413
414            ! Terms that do not involve already up-dated velocities.
415         
416            DO jj = k_j1+1, k_jpj-1
417               DO ji = 2, jpim1
418                  iim1 = ji
419                  ijm1 = jj - 1
420                  iip1 = ji + 1
421                  ijp1 = jj
422                  ze11 =   akappa(iim1,ijm1,1,1) * u_ice(iip1,ijp1) + akappa(iim1,ijm1,1,2) * v_ice(iip1,ijp1)
423                  ze12 = + akappa(iim1,ijm1,2,2) * u_ice(iip1,ijp1) - akappa(iim1,ijm1,2,1) * v_ice(iip1,ijp1)
424                  ze22 = + akappa(iim1,ijm1,2,2) * v_ice(iip1,ijp1) + akappa(iim1,ijm1,2,1) * u_ice(iip1,ijp1)
425                  ze21 =   akappa(iim1,ijm1,1,1) * v_ice(iip1,ijp1) - akappa(iim1,ijm1,1,2) * u_ice(iip1,ijp1)
426                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
427                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
428                  zvis12 =       zviseta (iim1,ijm1) + dm
429                  zvis21 =       zviseta (iim1,ijm1)
430                  zdiag = zvis22 * ( ze11 + ze22 )
431                  zs11(ji,jj,2,1) =  zvis11 * ze11 + zdiag
432                  zs12(ji,jj,2,1) =  zvis12 * ze12 + zvis21 * ze21
433                  zs22(ji,jj,2,1) =  zvis11 * ze22 + zdiag
434                  zs21(ji,jj,2,1) =  zvis12 * ze21 + zvis21 * ze12
435
436
437                  iim1 = ji - 1
438                  ijm1 = jj
439                  iip1 = ji
440                  ijp1 = jj + 1                   
441                  ze11 =   akappa(iim1,ijm1,1,1) * ( u_ice(iip1,ijp1) - u_ice(iim1,ijp1) )   &
442                     &   + akappa(iim1,ijm1,1,2) * ( v_ice(iip1,ijp1) + v_ice(iim1,ijp1) )
443                  ze12 = + akappa(iim1,ijm1,2,2) * ( u_ice(iim1,ijp1) + u_ice(iip1,ijp1) )   &
444                     &   - akappa(iim1,ijm1,2,1) * ( v_ice(iim1,ijp1) + v_ice(iip1,ijp1) )
445                  ze22 = + akappa(iim1,ijm1,2,2) * ( v_ice(iim1,ijp1) + v_ice(iip1,ijp1) )   &
446                     &   + akappa(iim1,ijm1,2,1) * ( u_ice(iim1,ijp1) + u_ice(iip1,ijp1) )
447                  ze21 =   akappa(iim1,ijm1,1,1) * ( v_ice(iip1,ijp1) - v_ice(iim1,ijp1) )   &
448                     &   - akappa(iim1,ijm1,1,2) * ( u_ice(iip1,ijp1) + u_ice(iim1,ijp1) )
449                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
450                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
451                  zvis12 =       zviseta (iim1,ijm1) + dm
452                  zvis21 =       zviseta (iim1,ijm1)
453                  zdiag = zvis22 * ( ze11 + ze22 )
454                  zs11(ji,jj,1,2) =  zvis11 * ze11 + zdiag
455                  zs12(ji,jj,1,2) =  zvis12 * ze12 + zvis21 * ze21
456                  zs22(ji,jj,1,2) =  zvis11 * ze22 + zdiag
457                  zs21(ji,jj,1,2) =  zvis12 * ze21 + zvis21 * ze12
458
459                  iim1 = ji
460                  ijm1 = jj
461                  iip1 = ji + 1
462                  ijp1 = jj + 1
463                  ze11 =   akappa(iim1,ijm1,1,1) * ( u_ice(iip1,ijm1) + u_ice(iip1,ijp1) - u_ice(iim1,ijp1) )   &
464                     &   + akappa(iim1,ijm1,1,2) * ( v_ice(iip1,ijm1) + v_ice(iip1,ijp1) + v_ice(iim1,ijp1) )
465                  ze12 = - akappa(iim1,ijm1,2,2) * ( u_ice(iip1,ijm1) - u_ice(iim1,ijp1) - u_ice(iip1,ijp1) )   &
466                     &   - akappa(iim1,ijm1,2,1) * ( v_ice(iip1,ijm1) + v_ice(iim1,ijp1) + v_ice(iip1,ijp1) )
467                  ze22 = - akappa(iim1,ijm1,2,2) * ( v_ice(iip1,ijm1) - v_ice(iim1,ijp1) - v_ice(iip1,ijp1) )   &
468                     &   + akappa(iim1,ijm1,2,1) * ( u_ice(iip1,ijm1) + u_ice(iim1,ijp1) + u_ice(iip1,ijp1) )
469                  ze21 =   akappa(iim1,ijm1,1,1) * ( v_ice(iip1,ijm1) + v_ice(iip1,ijp1) - v_ice(iim1,ijp1) )   &
470                     &   - akappa(iim1,ijm1,1,2) * ( u_ice(iip1,ijm1) + u_ice(iip1,ijp1) + u_ice(iim1,ijp1) ) 
471                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
472                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
473                  zvis12 =       zviseta (iim1,ijm1) + dm
474                  zvis21 =       zviseta (iim1,ijm1)
475
476                  zdiag = zvis22 * ( ze11 + ze22 )
477                  zs11(ji,jj,2,2) =  zvis11 * ze11 + zdiag
478                  zs12(ji,jj,2,2) =  zvis12 * ze12 + zvis21 * ze21
479                  zs22(ji,jj,2,2) =  zvis11 * ze22 + zdiag
480                  zs21(ji,jj,2,2) =  zvis12 * ze21 + zvis21 * ze12
481
482               END DO
483            END DO
484
485            ! Terms involving already up-dated velocities.
486            !-Using the arrays zu_ice and zv_ice in the computation of the terms ze leads to JACOBI's method;
487            ! Using arrays u and v in the computation of the terms ze leads to GAUSS-SEIDEL method.
488             
489            DO jj = k_j1+1, k_jpj-1
490               DO ji = 2, jpim1
491                  iim1 = ji - 1
492                  ijm1 = jj - 1
493                  iip1 = ji
494                  ijp1 = jj
495                  ze11 =   akappa(iim1,ijm1,1,1) * ( zu_ice(iip1,ijm1) - zu_ice(iim1,ijm1) - zu_ice(iim1,ijp1) )   &
496                     &   + akappa(iim1,ijm1,1,2) * ( zv_ice(iip1,ijm1) + zv_ice(iim1,ijm1) + zv_ice(iim1,ijp1) )
497                  ze12 = - akappa(iim1,ijm1,2,2) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) - zu_ice(iim1,ijp1) )   &
498                     &   - akappa(iim1,ijm1,2,1) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) + zv_ice(iim1,ijp1) )
499                  ze22 = - akappa(iim1,ijm1,2,2) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) - zv_ice(iim1,ijp1) )   &
500                     &   + akappa(iim1,ijm1,2,1) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) + zu_ice(iim1,ijp1) )
501                  ze21 =   akappa(iim1,ijm1,1,1) * ( zv_ice(iip1,ijm1) - zv_ice(iim1,ijm1) - zv_ice(iim1,ijp1) )   &
502                     &  -  akappa(iim1,ijm1,1,2) * ( zu_ice(iip1,ijm1) + zu_ice(iim1,ijm1) + zu_ice(iim1,ijp1) )
503                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
504                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
505                  zvis12 =       zviseta (iim1,ijm1) + dm
506                  zvis21 =       zviseta (iim1,ijm1)
507
508                  zdiag = zvis22 * ( ze11 + ze22 )
509                  zs11(ji,jj,1,1) =  zvis11 * ze11 + zdiag
510                  zs12(ji,jj,1,1) =  zvis12 * ze12 + zvis21 * ze21
511                  zs22(ji,jj,1,1) =  zvis11 * ze22 + zdiag
512                  zs21(ji,jj,1,1) =  zvis12 * ze21 + zvis21 * ze12
513
514#if defined key_agrif
515             END DO
516          END DO
517
518          DO jj = k_j1+1, k_jpj-1
519             DO ji = 2, jpim1
520#endif
521
522                  iim1 = ji
523                  ijm1 = jj - 1
524                  iip1 = ji + 1
525                  ze11 =   akappa(iim1,ijm1,1,1) * ( zu_ice(iip1,ijm1) - zu_ice(iim1,ijm1) )   &
526                     &   + akappa(iim1,ijm1,1,2) * ( zv_ice(iip1,ijm1) + zv_ice(iim1,ijm1) )
527                  ze12 = - akappa(iim1,ijm1,2,2) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) )   &
528                     &   - akappa(iim1,ijm1,2,1) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) )
529                  ze22 = - akappa(iim1,ijm1,2,2) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) )   &
530                     &   + akappa(iim1,ijm1,2,1) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) )
531                  ze21 =   akappa(iim1,ijm1,1,1) * ( zv_ice(iip1,ijm1) - zv_ice(iim1,ijm1) )   &
532                     &   - akappa(iim1,ijm1,1,2) * ( zu_ice(iip1,ijm1) + zu_ice(iim1,ijm1) )
533                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
534                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
535                  zvis12 =       zviseta (iim1,ijm1) + dm
536                  zvis21 =       zviseta (iim1,ijm1)
537
538                  zdiag = zvis22 * ( ze11 + ze22 )
539                  zs11(ji,jj,2,1) =  zs11(ji,jj,2,1) + zvis11 * ze11 + zdiag
540                  zs12(ji,jj,2,1) =  zs12(ji,jj,2,1) + zvis12 * ze12 + zvis21 * ze21
541                  zs22(ji,jj,2,1) =  zs22(ji,jj,2,1) + zvis11 * ze22 + zdiag
542                  zs21(ji,jj,2,1) =  zs21(ji,jj,2,1) + zvis12 * ze21 + zvis21 * ze12
543
544
545                  iim1 = ji - 1
546                  ijm1 = jj 
547                  ze11 = - akappa(iim1,ijm1,1,1) * zu_ice(iim1,ijm1) + akappa(iim1,ijm1,1,2) * zv_ice(iim1,ijm1)
548                  ze12 = - akappa(iim1,ijm1,2,2) * zu_ice(iim1,ijm1) - akappa(iim1,ijm1,2,1) * zv_ice(iim1,ijm1)
549                  ze22 = - akappa(iim1,ijm1,2,2) * zv_ice(iim1,ijm1) + akappa(iim1,ijm1,2,1) * zu_ice(iim1,ijm1)
550                  ze21 = - akappa(iim1,ijm1,1,1) * zv_ice(iim1,ijm1) - akappa(iim1,ijm1,1,2) * zu_ice(iim1,ijm1)
551                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
552                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
553                  zvis12 =       zviseta (iim1,ijm1) + dm
554                  zvis21 =       zviseta (iim1,ijm1)
555
556                  zdiag = zvis22 * ( ze11 + ze22 )
557                  zs11(ji,jj,1,2) =  zs11(ji,jj,1,2) + zvis11 * ze11 + zdiag 
558                  zs12(ji,jj,1,2) =  zs12(ji,jj,1,2) + zvis12 * ze12 + zvis21 * ze21
559                  zs22(ji,jj,1,2) =  zs22(ji,jj,1,2) + zvis11 * ze22 + zdiag
560                  zs21(ji,jj,1,2) =  zs21(ji,jj,1,2) + zvis12 * ze21 + zvis21 * ze12
561
562#if defined key_agrif
563             END DO
564          END DO
565
566          DO jj = k_j1+1, k_jpj-1
567             DO ji = 2, jpim1
568#endif
569                  zd1(ji,jj) =   &
570                     + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2)  &
571                     - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2)  &
572                     - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1)  &
573                     + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2)  &
574                     + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1)  &
575                     + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2)  &
576                     - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1)  &
577                     - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2)
578                  zd2(ji,jj) =   &
579                     + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2)  &
580                     - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2)  &
581                     - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1)  &
582                     + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2)  &
583                     - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1)  &
584                     - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2)  &
585                     + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1)  &
586                     + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2)
587               END DO
588            END DO
589
590            DO jj = k_j1+1, k_jpj-1
591               DO ji = 2, jpim1
592                  zunw = (  ( za1(ji,jj) + zd1(ji,jj) ) * zc2(ji,jj)        &
593                     &    + ( za2(ji,jj) + zd2(ji,jj) ) * zc1(ji,jj) ) * zden(ji,jj)
594
595                  zvnw = (  ( za2(ji,jj) + zd2(ji,jj) ) * zb1(ji,jj)        &
596                     &    - ( za1(ji,jj) + zd1(ji,jj) ) * zb2(ji,jj) ) * zden(ji,jj)
597
598                  zmask = ( 1.0 - MAX( rzero, SIGN( rone , 1.0 - zmass(ji,jj) ) ) ) * tmu(ji,jj)
599
600                  u_ice(ji,jj) = ( u_ice(ji,jj) + om * ( zunw - u_ice(ji,jj) ) * tmu(ji,jj) ) * zmask
601                  v_ice(ji,jj) = ( v_ice(ji,jj) + om * ( zvnw - v_ice(ji,jj) ) * tmu(ji,jj) ) * zmask
602               END DO
603            END DO
604
605            CALL lbc_lnk( u_ice, 'I', -1. )
606            CALL lbc_lnk( v_ice, 'I', -1. )
607
608            !---  5.2.5.4. Convergence test.
609            DO jj = k_j1+1 , k_jpj-1
610               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
611            END DO
612            zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) )
613            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
614
615            IF ( zresm <= resl) EXIT iflag
616
617         END DO iflag
618
619         zindu1 = 1.0 - zindu
620         DO jj = k_j1 , k_jpj-1
621            zu0(:,jj) = zindu * zu0(:,jj) + zindu1 * u_ice(:,jj)
622            zv0(:,jj) = zindu * zv0(:,jj) + zindu1 * v_ice(:,jj)
623         END DO
624      !                                                   ! ==================== !
625      END DO                                              !  end loop over iter  !
626      !                                                   ! ==================== !
627
628      IF(ln_ctl) THEN
629         WRITE(charout,FMT="('lim_rhg  : res =',D23.16, ' iter =',I4)") zresm, jter
630         CALL prt_ctl_info(charout)
631         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
632      ENDIF
633
634   END SUBROUTINE lim_rhg
635
636#else
637   !!----------------------------------------------------------------------
638   !!   Default option          Dummy module           NO LIM sea-ice model
639   !!----------------------------------------------------------------------
640CONTAINS
641   SUBROUTINE lim_rhg( k1 , k2 )         ! Dummy routine
642      WRITE(*,*) 'lim_rhg: You should not have seen this print! error?', k1, k2
643   END SUBROUTINE lim_rhg
644#endif
645
646   !!==============================================================================
647END MODULE limrhg
Note: See TracBrowser for help on using the repository browser.