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 branches/dev_001_SBC/NEMO/LIM_SRC – NEMO

source: branches/dev_001_SBC/NEMO/LIM_SRC/limrhg.F90 @ 755

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

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