source: branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 3978

Last change on this file since 3978 was 3978, checked in by clem, 7 years ago

change min thickness by min volume to set up u_ice=u_oce

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