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/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 2332

Last change on this file since 2332 was 2332, checked in by rblod, 13 years ago

correct small bugs in LIM2, see ticket #746

  • Property svn:keywords set to Id
File size: 37.4 KB
Line 
1MODULE limrhg
2   !!======================================================================
3   !!                     ***  MODULE  limrhg  ***
4   !!   Ice rheology : sea ice rheology
5   !!======================================================================
6   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
7   !!            3.0  !  2008-03  (M. Vancoppenolle) LIM3
8   !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy
9   !!            3.3  !  2009-05  (G.Garric) addition of the lim2_evp cas
10   !!----------------------------------------------------------------------
11#if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp )
12   !!----------------------------------------------------------------------
13   !!   'key_lim3'                                      LIM3 sea-ice model
14   !!----------------------------------------------------------------------
15   !!   lim_rhg   : computes ice velocities
16   !!----------------------------------------------------------------------
17   !! * Modules used
18   USE phycst
19   USE par_oce
20   USE dom_oce
21   USE sbc_oce         ! Surface boundary condition: ocean fields
22   USE sbc_ice         ! Surface boundary condition: ice fields
23   USE lbclnk
24   USE lib_mpp
25   USE in_out_manager  ! I/O manager
26   USE prtctl          ! Print control
27#if defined key_lim3
28   USE ice
29   USE dom_ice
30   USE iceini
31   USE limitd_me
32#endif
33#if defined key_lim2 && ! defined key_lim2_vp
34   USE ice_2            ! LIM2: ice variables
35   USE dom_ice_2        ! LIM2: ice domain
36   USE iceini_2         ! LIM2: ice initialisation
37#endif
38
39
40   IMPLICIT NONE
41   PRIVATE
42
43   !! * Routine accessibility
44   PUBLIC lim_rhg  ! routine called by lim_dyn
45
46   !! * Substitutions
47#  include "vectopt_loop_substitute.h90"
48
49   !! * Module variables
50   REAL(wp)  ::           &  ! constant values
51      rzero   = 0.e0   ,  &
52      rone    = 1.e0
53   !!----------------------------------------------------------------------
54   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
55   !! $Id$
56   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
57   !!----------------------------------------------------------------------
58
59CONTAINS
60
61   SUBROUTINE lim_rhg( k_j1, k_jpj )
62
63      !!-------------------------------------------------------------------
64      !!                 ***  SUBROUTINE lim_rhg  ***
65      !!                          EVP-C-grid
66      !!
67      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
68      !!  stress and sea-surface slope. Ice-ice interaction is described by
69      !!  a non-linear elasto-viscous-plastic (EVP) law including shear
70      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).   
71      !!
72      !!  The points in the C-grid look like this, dear reader
73      !!
74      !!                              (ji,jj)
75      !!                                 |
76      !!                                 |
77      !!                      (ji-1,jj)  |  (ji,jj)
78      !!                             ---------   
79      !!                            |         |
80      !!                            | (ji,jj) |------(ji,jj)
81      !!                            |         |
82      !!                             ---------   
83      !!                     (ji-1,jj-1)     (ji,jj-1)
84      !!
85      !! ** Inputs  : - wind forcing (stress), oceanic currents
86      !!                ice total volume (vt_i) per unit area
87      !!                snow total volume (vt_s) per unit area
88      !!
89      !! ** Action  : - compute u_ice, v_ice : the components of the
90      !!                sea-ice velocity vector
91      !!              - compute delta_i, shear_i, divu_i, which are inputs
92      !!                of the ice thickness distribution
93      !!
94      !! ** Steps   : 1) Compute ice snow mass, ice strength
95      !!              2) Compute wind, oceanic stresses, mass terms and
96      !!                 coriolis terms of the momentum equation
97      !!              3) Solve the momentum equation (iterative procedure)
98      !!              4) Prevent high velocities if the ice is thin
99      !!              5) Recompute invariants of the strain rate tensor
100      !!                 which are inputs of the ITD, store stress
101      !!                 for the next time step
102      !!              6) Control prints of residual (convergence)
103      !!                 and charge ellipse.
104      !!                 The user should make sure that the parameters
105      !!                 nevp, telast and creepl maintain stress state
106      !!                 on the charge ellipse for plastic flow
107      !!                 e.g. in the Canadian Archipelago
108      !!
109      !! ** References : Hunke and Dukowicz, JPO97
110      !!                 Bouillon et al., 08, in prep (update this when
111      !!                 published)
112      !!                 Vancoppenolle et al., OM08
113      !!
114      !!-------------------------------------------------------------------
115      ! * Arguments
116      !
117      INTEGER, INTENT(in) :: &
118         k_j1 ,                      & !: southern j-index for ice computation
119         k_jpj                         !: northern j-index for ice computation
120
121      ! * Local variables
122      INTEGER ::   ji, jj              !: dummy loop indices
123
124      INTEGER  :: &
125         jter                          !: temporary integers
126
127      CHARACTER (len=50) ::   charout
128
129      REAL(wp) :: &
130         zt11, zt12, zt21, zt22,     & !: temporary scalars
131         ztagnx, ztagny,             & !: wind stress on U/V points                       
132         delta                         !
133
134      REAL(wp) :: &
135         za,                         & !:
136         zstms,                      & !: temporary scalar for ice strength
137         zsang,                      & !: temporary scalar for coriolis term
138         zmask                         !: mask for the computation of ice mass
139
140      REAL(wp),DIMENSION(jpi,jpj) :: &
141         zpresh        ,             & !: temporary array for ice strength
142         zpreshc       ,             & !: Ice strength on grid cell corners (zpreshc)
143         zfrld1, zfrld2,             & !: lead fraction on U/V points                                   
144         zmass1, zmass2,             & !: ice/snow mass on U/V points                                   
145         zcorl1, zcorl2,             & !: coriolis parameter on U/V points
146         za1ct, za2ct  ,             & !: temporary arrays
147         zc1           ,             & !: ice mass
148         zusw          ,             & !: temporary weight for the computation
149                                !: of ice strength
150         u_oce1, v_oce1,             & !: ocean u/v component on U points                           
151         u_oce2, v_oce2,             & !: ocean u/v component on V points
152         u_ice2,                     & !: ice u component on V point
153         v_ice1                        !: ice v component on U point
154
155      REAL(wp) :: &
156         dtevp,                      & !: time step for subcycling
157         dtotel,                     & !:
158         ecc2,                       & !: square of yield ellipse eccenticity
159         z0,                         & !: temporary scalar
160         zr,                         & !: temporary scalar
161         zcca, zccb,                 & !: temporary scalars
162         zu_ice2,                    & !:
163         zv_ice1,                    & !:
164         zddc, zdtc,                 & !: temporary array for delta on corners
165         zdst,                       & !: temporary array for delta on centre
166         zdsshx, zdsshy,             & !: term for the gradient of ocean surface
167         sigma1, sigma2                !: internal ice stress
168
169      REAL(wp),DIMENSION(jpi,jpj) :: &
170         zf1, zf2                      !: arrays for internal stresses
171
172      REAL(wp),DIMENSION(jpi,jpj) :: &
173         zdd, zdt,                   & !: Divergence and tension at centre of grid cells
174         zds,                        & !: Shear on northeast corner of grid cells
175         deltat,                     & !: Delta at centre of grid cells
176         deltac,                     & !: Delta on corners
177         zs1, zs2,                   & !: Diagonal stress tensor components zs1 and zs2
178         zs12                          !: Non-diagonal stress tensor component zs12
179
180      REAL(wp) :: &
181         zresm            ,          & !: Maximal error on ice velocity
182         zindb            ,          & !: ice (1) or not (0)     
183         zdummy                        !: dummy argument
184
185      REAL(wp),DIMENSION(jpi,jpj) :: &
186         zu_ice           ,          & !: Ice velocity on previous time step
187         zv_ice           ,          &
188         zresr                         !: Local error on velocity
189      !!-------------------------------------------------------------------
190
191#if  defined key_lim2 && ! defined key_lim2_vp
192     vt_s => hsnm
193     vt_i => hicm
194     at_i(:,:) = 1. - frld(:,:)
195#endif
196      !
197      !------------------------------------------------------------------------------!
198      ! 1) Ice-Snow mass (zc1), ice strength (zpresh)                                !
199      !------------------------------------------------------------------------------!
200      !
201      ! Put every vector to 0
202      zpresh(:,:) = 0.0 ; zc1(:,:) = 0.0
203      zpreshc(:,:) = 0.0
204      u_ice2(:,:)  = 0.0 ; v_ice1(:,:)  = 0.0
205      zdd(:,:)     = 0.0 ; zdt(:,:)     = 0.0 ; zds(:,:)     = 0.0
206
207#if defined key_lim3
208      ! Ice strength on T-points
209      CALL lim_itd_me_icestrength(ridge_scheme_swi)
210#endif
211
212      ! Ice mass and temp variables
213!CDIR NOVERRCHK
214      DO jj = k_j1 , k_jpj
215!CDIR NOVERRCHK
216         DO ji = 1 , jpi
217            zc1(ji,jj)    = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) )
218#if defined key_lim3
219            zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj) / 2.
220#endif
221            ! tmi = 1 where there is ice or on land
222            tmi(ji,jj)    = 1.0 - ( 1.0 - MAX( 0.0 , SIGN ( 1.0 , vt_i(ji,jj) - &
223               epsd ) ) ) * tms(ji,jj)
224         END DO
225      END DO
226
227      ! Ice strength on grid cell corners (zpreshc)
228      ! needed for calculation of shear stress
229!CDIR NOVERRCHK
230      DO jj = k_j1+1, k_jpj-1
231!CDIR NOVERRCHK
232         DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1)
233            zstms          =  tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
234               &              tms(ji,jj+1)   * wght(ji+1,jj+1,1,2) + &
235               &              tms(ji+1,jj)   * wght(ji+1,jj+1,2,1) + &
236               &              tms(ji,jj)     * wght(ji+1,jj+1,1,1)
237            zusw(ji,jj)    = 1.0 / MAX( zstms, epsd )
238            zpreshc(ji,jj) = (  zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
239               &                zpresh(ji,jj+1)   * wght(ji+1,jj+1,1,2) + &
240               &                zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + & 
241               &                zpresh(ji,jj)     * wght(ji+1,jj+1,1,1)   &
242               &             ) * zusw(ji,jj)
243         END DO
244      END DO
245
246      CALL lbc_lnk( zpreshc(:,:), 'F', 1. )
247      !
248      !------------------------------------------------------------------------------!
249      ! 2) Wind / ocean stress, mass terms, coriolis terms
250      !------------------------------------------------------------------------------!
251      !
252      !  Wind stress, coriolis and mass terms on the sides of the squares       
253      !  zfrld1: lead fraction on U-points                                     
254      !  zfrld2: lead fraction on V-points                                     
255      !  zmass1: ice/snow mass on U-points                                   
256      !  zmass2: ice/snow mass on V-points                                   
257      !  zcorl1: Coriolis parameter on U-points                             
258      !  zcorl2: Coriolis parameter on V-points                           
259      !  (ztagnx,ztagny): wind stress on U/V points                       
260      !  u_oce1: ocean u component on u points                           
261      !  v_oce1: ocean v component on u points                         
262      !  u_oce2: ocean u component on v points                         
263      !  v_oce2: ocean v component on v points                       
264
265      DO jj = k_j1+1, k_jpj-1
266         DO ji = fs_2, fs_jpim1
267
268            zt11 = tms(ji,jj)*e1t(ji,jj)
269            zt12 = tms(ji+1,jj)*e1t(ji+1,jj)
270            zt21 = tms(ji,jj)*e2t(ji,jj)
271            zt22 = tms(ji,jj+1)*e2t(ji,jj+1)
272
273            ! Leads area.
274            zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + &
275               &                        zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd )
276            zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + &
277               &                        zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd )
278
279            ! Mass, coriolis coeff. and currents
280            zmass1(ji,jj) = ( zt12*zc1(ji,jj) + zt11*zc1(ji+1,jj) ) / (zt11+zt12+epsd)
281            zmass2(ji,jj) = ( zt22*zc1(ji,jj) + zt21*zc1(ji,jj+1) ) / (zt21+zt22+epsd)
282            zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + &
283               e1t(ji,jj)*fcor(ji+1,jj) ) &
284               / (e1t(ji,jj) + e1t(ji+1,jj) + epsd )
285            zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)*fcor(ji,jj) + &
286               e2t(ji,jj)*fcor(ji,jj+1) ) &
287               / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd )
288            !
289            u_oce1(ji,jj)  = u_oce(ji,jj)
290            v_oce2(ji,jj)  = v_oce(ji,jj)
291
292            ! Ocean has no slip boundary condition
293            v_oce1(ji,jj)  = 0.5*( (v_oce(ji,jj)+v_oce(ji,jj-1))*e1t(ji,jj)    &
294               &                 +(v_oce(ji+1,jj)+v_oce(ji+1,jj-1))*e1t(ji+1,jj)) &
295               &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
296
297            u_oce2(ji,jj)  = 0.5*((u_oce(ji,jj)+u_oce(ji-1,jj))*e2t(ji,jj)     &
298               &                 +(u_oce(ji,jj+1)+u_oce(ji-1,jj+1))*e2t(ji,jj+1)) &
299               &                / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
300
301            ! Wind stress at U,V-point
302            ztagnx = ( 1. - zfrld1(ji,jj) ) * utau_ice(ji,jj)
303            ztagny = ( 1. - zfrld2(ji,jj) ) * vtau_ice(ji,jj)
304
305            ! Computation of the velocity field taking into account the ice internal interaction.
306            ! Terms that are independent of the velocity field.
307
308            ! SB On utilise maintenant le gradient de la pente de l'ocean
309            ! include it later
310
311            zdsshx =  (ssh_m(ji+1,jj) - ssh_m(ji,jj))/e1u(ji,jj)
312            zdsshy =  (ssh_m(ji,jj+1) - ssh_m(ji,jj))/e2v(ji,jj)
313
314            za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx
315            za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy
316
317         END DO
318      END DO
319
320      !
321      !------------------------------------------------------------------------------!
322      ! 3) Solution of the momentum equation, iterative procedure
323      !------------------------------------------------------------------------------!
324      !
325      ! Time step for subcycling
326      dtevp  = rdt_ice / nevp
327      dtotel = dtevp / ( 2.0 * telast )
328
329      !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)
330      ecc2 = ecc*ecc
331
332      !-Initialise stress tensor
333      zs1(:,:)  = stress1_i(:,:) 
334      zs2(:,:)  = stress2_i(:,:)
335      zs12(:,:) = stress12_i(:,:)
336
337      !----------------------!
338      DO jter = 1 , nevp                              !    loop over jter    !
339         !----------------------!       
340         DO jj = k_j1, k_jpj-1
341            zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
342            zv_ice(:,jj) = v_ice(:,jj)
343         END DO
344
345         DO jj = k_j1+1, k_jpj-1
346            DO ji = fs_2, jpim1   !RB bug no vect opt due to tmi
347
348               
349               !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002)
350               !- zdd(:,:), zdt(:,:): divergence and tension at centre of grid cells
351               !- zds(:,:): shear on northeast corner of grid cells
352               !
353               !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded,
354               !                      there are many repeated calculations.
355               !                      Speed could be improved by regrouping terms. For
356               !                      the moment, however, the stress is on clarity of coding to avoid
357               !                      bugs (Martin, for Miguel).
358               !
359               !- ALSO: arrays zdd, zdt, zds and delta could
360               !  be removed in the future to minimise memory demand.
361               !
362               !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of
363               !              grid cells, exactly as in the B grid case. For simplicity, the indexation on
364               !              the corners is the same as in the B grid.
365               !
366               !
367               zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      &
368                  &          -e2u(ji-1,jj)*u_ice(ji-1,jj)                  &
369                  &          +e1v(ji,jj)*v_ice(ji,jj)                      &
370                  &          -e1v(ji,jj-1)*v_ice(ji,jj-1)                  &
371                  &          )                                             &
372                  &         / area(ji,jj)
373
374               zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj)                    &
375                  &            -u_ice(ji-1,jj)/e2u(ji-1,jj)                &
376                  &           )*e2t(ji,jj)*e2t(ji,jj)                      &
377                  &          -( v_ice(ji,jj)/e1v(ji,jj)                    &
378                  &            -v_ice(ji,jj-1)/e1v(ji,jj-1)                &
379                  &           )*e1t(ji,jj)*e1t(ji,jj)                      &
380                  &         )                                              &
381                  &        / area(ji,jj)
382
383               !
384               zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1)                &
385                  &            -u_ice(ji,jj)/e1u(ji,jj)                    &
386                  &           )*e1f(ji,jj)*e1f(ji,jj)                      &
387                  &          +( v_ice(ji+1,jj)/e2v(ji+1,jj)                &
388                  &            -v_ice(ji,jj)/e2v(ji,jj)                    &
389                  &           )*e2f(ji,jj)*e2f(ji,jj)                      &
390                  &         )                                              &
391                  &        / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) &
392                  &        * tmi(ji,jj) * tmi(ji,jj+1)                     &
393                  &        * tmi(ji+1,jj) * tmi(ji+1,jj+1)
394
395
396               v_ice1(ji,jj)  = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)   &
397                  &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
398                  &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
399
400               u_ice2(ji,jj)  = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)   &
401                  &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &
402                  &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
403
404            END DO
405         END DO
406         CALL lbc_lnk( v_ice1(:,:), 'U', -1. )
407         CALL lbc_lnk( u_ice2(:,:), 'V', -1. )
408
409!CDIR NOVERRCHK
410         DO jj = k_j1+1, k_jpj-1
411!CDIR NOVERRCHK
412            DO ji = fs_2, fs_jpim1
413
414               !- Calculate Delta at centre of grid cells
415               zdst       = (  e2u( ji  , jj   ) * v_ice1(ji,jj)            &
416                  &          - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)          &
417                  &          + e1v( ji  , jj   ) * u_ice2(ji,jj)            &
418                  &          - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1)          &
419                  &          )                                              &
420                  &         / area(ji,jj)
421
422               delta = SQRT( zdd(ji,jj)*zdd(ji,jj) +                        & 
423                  &                       ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 
424               deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 +                    &
425                  (zdt(ji,jj)**2 + zdst**2)*usecc2), creepl )
426
427               !-Calculate stress tensor components zs1 and zs2
428               !-at centre of grid cells (see section 3.5 of CICE user's guide).
429               zs1(ji,jj) = ( zs1(ji,jj) &
430                  &          - dtotel*( ( 1.0 - alphaevp) * zs1(ji,jj) +    &
431                  &            ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) &
432                  * zpresh(ji,jj) ) )                          &       
433                  &        / ( 1.0 + alphaevp * dtotel )
434
435               zs2(ji,jj) = ( zs2(ji,jj)   &
436                  &          - dtotel*((1.0-alphaevp)*ecc2*zs2(ji,jj) -  &
437                  zdt(ji,jj)/deltat(ji,jj)*zpresh(ji,jj)) ) &
438                  &        / ( 1.0 + alphaevp*ecc2*dtotel )
439
440            END DO
441         END DO
442
443         CALL lbc_lnk( zs1(:,:), 'T', 1. )
444         CALL lbc_lnk( zs2(:,:), 'T', 1. )
445
446!CDIR NOVERRCHK
447         DO jj = k_j1+1, k_jpj-1
448!CDIR NOVERRCHK
449            DO ji = fs_2, fs_jpim1
450               !- Calculate Delta on corners
451               zddc  =      ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1)                &
452                  &            -v_ice1(ji,jj)/e1u(ji,jj)                    &
453                  &           )*e1f(ji,jj)*e1f(ji,jj)                       &
454                  &          +( u_ice2(ji+1,jj)/e2v(ji+1,jj)                &
455                  &            -u_ice2(ji,jj)/e2v(ji,jj)                    &
456                  &           )*e2f(ji,jj)*e2f(ji,jj)                       &
457                  &         )                                               &
458                  &        / ( e1f(ji,jj) * e2f(ji,jj) )
459
460               zdtc  =      (-( v_ice1(ji,jj+1)/e1u(ji,jj+1)                &
461                  &            -v_ice1(ji,jj)/e1u(ji,jj)                    &
462                  &           )*e1f(ji,jj)*e1f(ji,jj)                       &
463                  &          +( u_ice2(ji+1,jj)/e2v(ji+1,jj)                &
464                  &            -u_ice2(ji,jj)/e2v(ji,jj)                    &
465                  &           )*e2f(ji,jj)*e2f(ji,jj)                       &
466                  &         )                                               &
467                  &        / ( e1f(ji,jj) * e2f(ji,jj) )
468
469               deltac(ji,jj) = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl
470
471               !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide).
472               zs12(ji,jj) = ( zs12(ji,jj)      &
473                  &        - dtotel*( (1.0-alphaevp)*ecc2*zs12(ji,jj) - zds(ji,jj) / &
474                  &          ( 2.0*deltac(ji,jj) ) * zpreshc(ji,jj))) &
475                  &         / ( 1.0 + alphaevp*ecc2*dtotel ) 
476
477            END DO ! ji
478         END DO ! jj
479
480         CALL lbc_lnk( zs12(:,:), 'F', 1. )
481
482         ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002)
483         DO jj = k_j1+1, k_jpj-1
484            DO ji = fs_2, fs_jpim1
485               !- contribution of zs1, zs2 and zs12 to zf1
486               zf1(ji,jj) = 0.5*( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj) &
487                  &              +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj) &
488                  &              +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj) &
489                  &             ) / ( e1u(ji,jj)*e2u(ji,jj) )
490               ! contribution of zs1, zs2 and zs12 to zf2
491               zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) &
492                  &              -(zs2(ji,jj+1)*e1t(ji,jj+1)**2 - zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) &
493                  &              + 2.0*(zs12(ji,jj)*e2f(ji,jj)**2 -    &
494                  zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) &
495                  &             ) / ( e1v(ji,jj)*e2v(ji,jj) )
496            END DO
497         END DO
498         !
499         ! Computation of ice velocity
500         !
501         ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly.
502         !
503         IF (MOD(jter,2).eq.0) THEN 
504
505!CDIR NOVERRCHK
506            DO jj = k_j1+1, k_jpj-1
507!CDIR NOVERRCHK
508               DO ji = fs_2, fs_jpim1
509                  zmask        = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj)
510                  zsang        = SIGN ( 1.0 , fcor(ji,jj) ) * sangvg
511                  z0           = zmass1(ji,jj)/dtevp
512
513                  ! SB modif because ocean has no slip boundary condition
514                  zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)         &
515                     &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   &
516                     &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
517                  za           = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + &
518                     (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj))
519                  zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
520                     za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj))
521                  zcca         = z0+za*cangvg
522                  zccb         = zcorl1(ji,jj)+za*zsang
523                  u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 
524
525               END DO
526            END DO
527
528            CALL lbc_lnk( u_ice(:,:), 'U', -1. )
529
530!CDIR NOVERRCHK
531            DO jj = k_j1+1, k_jpj-1
532!CDIR NOVERRCHK
533               DO ji = fs_2, fs_jpim1
534
535                  zmask        = (1.0-MAX(rzero,SIGN(rone,-zmass2(ji,jj))))*tmv(ji,jj)
536                  zsang        = SIGN(1.0,fcor(ji,jj))*sangvg
537                  z0           = zmass2(ji,jj)/dtevp
538                  ! SB modif because ocean has no slip boundary condition
539                  zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj)     &
540                     &                 + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1))   &
541                     &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
542                  za           = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 
543                     (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj))
544                  zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + &
545                     za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj))
546                  zcca         = z0+za*cangvg
547                  zccb         = zcorl2(ji,jj)+za*zsang
548                  v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask
549
550               END DO
551            END DO
552
553            CALL lbc_lnk( v_ice(:,:), 'V', -1. )
554
555         ELSE 
556!CDIR NOVERRCHK
557            DO jj = k_j1+1, k_jpj-1
558!CDIR NOVERRCHK
559               DO ji = fs_2, fs_jpim1
560                  zmask        = (1.0-MAX(rzero,SIGN(rone,-zmass2(ji,jj))))*tmv(ji,jj)
561                  zsang        = SIGN(1.0,fcor(ji,jj))*sangvg
562                  z0           = zmass2(ji,jj)/dtevp
563                  ! SB modif because ocean has no slip boundary condition
564                  zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj)      &
565                     &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1))   &
566                     &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)   
567
568                  za           = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + &
569                     (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj))
570                  zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + &
571                     za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj))
572                  zcca         = z0+za*cangvg
573                  zccb         = zcorl2(ji,jj)+za*zsang
574                  v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask
575
576               END DO
577            END DO
578
579            CALL lbc_lnk( v_ice(:,:), 'V', -1. )
580
581!CDIR NOVERRCHK
582            DO jj = k_j1+1, k_jpj-1
583!CDIR NOVERRCHK
584               DO ji = fs_2, fs_jpim1
585                  zmask        = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj)
586                  zsang        = SIGN(1.0,fcor(ji,jj))*sangvg
587                  z0           = zmass1(ji,jj)/dtevp
588                  ! SB modif because ocean has no slip boundary condition
589                  ! GG Bug
590                  !                   zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)      &
591                  !                      &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj))   &
592                  !                      &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
593                  zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)      &
594                     &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   &
595                     &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
596
597                  za           = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + &
598                     (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj))
599                  zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
600                     za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj))
601                  zcca         = z0+za*cangvg
602                  zccb         = zcorl1(ji,jj)+za*zsang
603                  u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 
604               END DO ! ji
605            END DO ! jj
606
607            CALL lbc_lnk( u_ice(:,:), 'U', -1. )
608
609         ENDIF
610
611         IF(ln_ctl) THEN
612            !---  Convergence test.
613            DO jj = k_j1+1 , k_jpj-1
614               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) ,           &
615                  ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
616            END DO
617            zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) )
618            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
619         ENDIF
620
621         !                                                   ! ==================== !
622      END DO                                              !  end loop over jter  !
623      !                                                   ! ==================== !
624
625      !
626      !------------------------------------------------------------------------------!
627      ! 4) Prevent ice velocities when the ice is thin
628      !------------------------------------------------------------------------------!
629      !
630      ! If the ice thickness is below 1cm then ice velocity should equal the
631      ! ocean velocity,
632      ! This prevents high velocity when ice is thin
633!CDIR NOVERRCHK
634      DO jj = k_j1+1, k_jpj-1
635!CDIR NOVERRCHK
636         DO ji = fs_2, fs_jpim1
637            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 
638            zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
639            IF ( zdummy .LE. 5.0e-2 ) THEN
640               u_ice(ji,jj) = u_oce(ji,jj)
641               v_ice(ji,jj) = v_oce(ji,jj)
642            ENDIF ! zdummy
643         END DO
644      END DO
645
646      CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
647      CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
648
649      DO jj = k_j1+1, k_jpj-1 
650         DO ji = fs_2, fs_jpim1
651            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 
652            zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
653            IF ( zdummy .LE. 5.0e-2 ) THEN
654               v_ice1(ji,jj)  = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)   &
655                  &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
656                  &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
657
658               u_ice2(ji,jj)  = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)   &
659                  &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &
660                  &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
661            ENDIF ! zdummy
662         END DO
663      END DO
664
665      CALL lbc_lnk( u_ice2(:,:), 'V', -1. ) 
666      CALL lbc_lnk( v_ice1(:,:), 'U', -1. )
667
668      ! Recompute delta, shear and div, inputs for mechanical redistribution
669!CDIR NOVERRCHK
670      DO jj = k_j1+1, k_jpj-1
671!CDIR NOVERRCHK
672         DO ji = fs_2, jpim1   !RB bug no vect opt due to tmi
673            !- zdd(:,:), zdt(:,:): divergence and tension at centre
674            !- zds(:,:): shear on northeast corner of grid cells
675            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 
676            zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
677
678            IF ( zdummy .LE. 5.0e-2 ) THEN
679
680               zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      &
681                  &          -e2u(ji-1,jj)*u_ice(ji-1,jj)                  &
682                  &          +e1v(ji,jj)*v_ice(ji,jj)                      &
683                  &          -e1v(ji,jj-1)*v_ice(ji,jj-1)                  &
684                  &         )                                              &
685                  &         / area(ji,jj)
686
687               zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj)                    &
688                  &            -u_ice(ji-1,jj)/e2u(ji-1,jj)                &
689                  &           )*e2t(ji,jj)*e2t(ji,jj)                      &
690                  &          -( v_ice(ji,jj)/e1v(ji,jj)                    &
691                  &            -v_ice(ji,jj-1)/e1v(ji,jj-1)                &
692                  &           )*e1t(ji,jj)*e1t(ji,jj)                      &
693                  &         )                                              &
694                  &        / area(ji,jj)
695               !
696               ! SB modif because ocean has no slip boundary condition
697               zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1)              &
698                  &           - u_ice(ji,jj)   / e1u(ji,jj) )              &
699                  &           * e1f(ji,jj) * e1f(ji,jj)                    &
700                  &          + ( v_ice(ji+1,jj) / e2v(ji+1,jj)             &
701                  &            - v_ice(ji,jj)  / e2v(ji,jj) )              &
702                  &           * e2f(ji,jj) * e2f(ji,jj) )                  &
703                  &        / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) &
704                  &        * tmi(ji,jj) * tmi(ji,jj+1)                     &
705                  &        * tmi(ji+1,jj) * tmi(ji+1,jj+1)
706
707               zdst       = (  e2u( ji  , jj   ) * v_ice1(ji,jj)          &
708                  &          - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)         &
709                  &          + e1v( ji  , jj   ) * u_ice2(ji,jj)           &
710                  &          - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1)         & 
711                  &          )                                             &
712                  &         / area(ji,jj)
713
714               deltat(ji,jj) = SQRT( zdd(ji,jj)*zdd(ji,jj) + & 
715                  &                          ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 & 
716                  &                          ) + creepl
717
718            ENDIF ! zdummy
719
720         END DO !jj
721      END DO !ji
722      !
723      !------------------------------------------------------------------------------!
724      ! 5) Store stress tensor and its invariants
725      !------------------------------------------------------------------------------!
726      !
727      ! * Invariants of the stress tensor are required for limitd_me
728      ! accelerates convergence and improves stability
729      DO jj = k_j1+1, k_jpj-1
730         DO ji = fs_2, fs_jpim1
731            divu_i (ji,jj) = zdd   (ji,jj)
732            delta_i(ji,jj) = deltat(ji,jj)
733            shear_i(ji,jj) = zds   (ji,jj)
734         END DO
735      END DO
736
737      ! Lateral boundary condition
738      CALL lbc_lnk( divu_i (:,:), 'T', 1. )
739      CALL lbc_lnk( delta_i(:,:), 'T', 1. )
740      CALL lbc_lnk( shear_i(:,:), 'F', 1. )
741
742      ! * Store the stress tensor for the next time step
743      stress1_i (:,:) = zs1 (:,:)
744      stress2_i (:,:) = zs2 (:,:)
745      stress12_i(:,:) = zs12(:,:)
746
747      !
748      !------------------------------------------------------------------------------!
749      ! 6) Control prints of residual and charge ellipse
750      !------------------------------------------------------------------------------!
751      !
752      ! print the residual for convergence
753      IF(ln_ctl) THEN
754         WRITE(charout,FMT="('lim_rhg  : res =',D23.16, ' iter =',I4)") zresm, jter
755         CALL prt_ctl_info(charout)
756         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
757      ENDIF
758
759      ! print charge ellipse
760      ! This can be desactivated once the user is sure that the stress state
761      ! lie on the charge ellipse. See Bouillon et al. 08 for more details
762      IF(ln_ctl) THEN
763         CALL prt_ctl_info('lim_rhg  : numit  :',ivar1=numit)
764         CALL prt_ctl_info('lim_rhg  : nwrite :',ivar1=nwrite)
765         CALL prt_ctl_info('lim_rhg  : MOD    :',ivar1=MOD(numit,nwrite))
766         IF( MOD(numit,nwrite) .EQ. 0 ) THEN
767            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. '
768            CALL prt_ctl_info(charout)
769            DO jj = k_j1+1, k_jpj-1
770               DO ji = 2, jpim1
771                  IF (zpresh(ji,jj) .GT. 1.0) THEN
772                     sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
773                     sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )
774                     WRITE(charout,FMT="('lim_rhg  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
775                     CALL prt_ctl_info(charout)
776                  ENDIF
777               END DO
778            END DO
779            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 2000, numit, 0, 0, ' ch. ell. '
780            CALL prt_ctl_info(charout)
781         ENDIF
782      ENDIF
783
784   END SUBROUTINE lim_rhg
785
786#else
787   !!----------------------------------------------------------------------
788   !!   Default option          Dummy module           NO LIM sea-ice model
789   !!----------------------------------------------------------------------
790CONTAINS
791   SUBROUTINE lim_rhg( k1 , k2 )         ! Dummy routine
792      WRITE(*,*) 'lim_rhg: You should not have seen this print! error?', k1, k2
793   END SUBROUTINE lim_rhg
794#endif
795
796   !!==============================================================================
797END MODULE limrhg
Note: See TracBrowser for help on using the repository browser.