source: trunk/NEMO/LIM_SRC_3/limrhg.F90 @ 834

Last change on this file since 834 was 834, checked in by ctlod, 13 years ago

Clean comments and useless lines, see ticket:#72

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