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

source: branches/2012/dev_LOCEAN_2012/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 3584

Last change on this file since 3584 was 3584, checked in by cetlod, 11 years ago

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