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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 4444

Last change on this file since 4444 was 3837, checked in by trackstand2, 11 years ago

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