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

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 2717

Last change on this file since 2717 was 2717, checked in by flavoni, 13 years ago

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