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

source: branches/2011/dev_r2802_LOCEAN10_agrif_lim/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 2804

Last change on this file since 2804 was 2804, checked in by rblod, 13 years ago

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