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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 2590

Last change on this file since 2590 was 2590, checked in by trackstand2, 13 years ago

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