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

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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90 @ 2319

Last change on this file since 2319 was 2319, checked in by cbricaud, 13 years ago

put new EVP rheology lost during the merge

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