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/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 9 years ago

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