source: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 5086

Last change on this file since 5086 was 5086, checked in by timgraham, 7 years ago

Merged head of trunk into branch in preparation for putting code back onto the trunk
In working copy ran the command:
svn merge svn+sshtimgraham@…/ipsl/forge/projets/nemo/svn/trunk

Also recompiled NEMO_book.pdf with merged input files

  • 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   !!            3.4  !  2011-01  (A. Porter)  dynamical allocation
11   !!            3.5  !  2012-08  (R. Benshila)  AGRIF
12   !!----------------------------------------------------------------------
13#if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp )
14   !!----------------------------------------------------------------------
15   !!   'key_lim3'               OR                     LIM-3 sea-ice model
16   !!   'key_lim2' AND NOT 'key_lim2_vp'            EVP LIM-2 sea-ice model
17   !!----------------------------------------------------------------------
18   !!   lim_rhg       : computes ice velocities
19   !!----------------------------------------------------------------------
20   USE phycst         ! Physical constant
21   USE oce     , ONLY :  snwice_mass, snwice_mass_b
22   USE par_oce        ! Ocean parameters
23   USE dom_oce        ! Ocean domain
24   USE sbc_oce        ! Surface boundary condition: ocean fields
25   USE sbc_ice        ! Surface boundary condition: ice fields
26#if defined key_lim3
27   USE ice            ! LIM-3: ice variables
28   USE dom_ice        ! LIM-3: ice domain
29   USE limitd_me      ! LIM-3:
30#else
31   USE ice_2          ! LIM-2: ice variables
32   USE dom_ice_2      ! LIM-2: ice domain
33#endif
34   USE lbclnk         ! Lateral Boundary Condition / MPP link
35   USE lib_mpp        ! MPP library
36   USE wrk_nemo       ! work arrays
37   USE in_out_manager ! I/O manager
38   USE prtctl         ! Print control
39   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
40#if defined key_agrif && defined key_lim2
41   USE agrif_lim2_interp
42#endif
43#if defined key_bdy
44   USE bdyice_lim
45#endif
46
47   IMPLICIT NONE
48   PRIVATE
49
50   PUBLIC   lim_rhg        ! routine called by lim_dyn (or lim_dyn_2)
51
52   !! * Substitutions
53#  include "vectopt_loop_substitute.h90"
54   !!----------------------------------------------------------------------
55   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
56   !! $Id$
57   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   SUBROUTINE lim_rhg( k_j1, k_jpj )
62      !!-------------------------------------------------------------------
63      !!                 ***  SUBROUTINE lim_rhg  ***
64      !!                          EVP-C-grid
65      !!
66      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
67      !!  stress and sea-surface slope. Ice-ice interaction is described by
68      !!  a non-linear elasto-viscous-plastic (EVP) law including shear
69      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).   
70      !!
71      !!  The points in the C-grid look like this, dear reader
72      !!
73      !!                              (ji,jj)
74      !!                                 |
75      !!                                 |
76      !!                      (ji-1,jj)  |  (ji,jj)
77      !!                             ---------   
78      !!                            |         |
79      !!                            | (ji,jj) |------(ji,jj)
80      !!                            |         |
81      !!                             ---------   
82      !!                     (ji-1,jj-1)     (ji,jj-1)
83      !!
84      !! ** Inputs  : - wind forcing (stress), oceanic currents
85      !!                ice total volume (vt_i) per unit area
86      !!                snow total volume (vt_s) per unit area
87      !!
88      !! ** Action  : - compute u_ice, v_ice : the components of the
89      !!                sea-ice velocity vector
90      !!              - compute delta_i, shear_i, divu_i, which are inputs
91      !!                of the ice thickness distribution
92      !!
93      !! ** Steps   : 1) Compute ice snow mass, ice strength
94      !!              2) Compute wind, oceanic stresses, mass terms and
95      !!                 coriolis terms of the momentum equation
96      !!              3) Solve the momentum equation (iterative procedure)
97      !!              4) Prevent high velocities if the ice is thin
98      !!              5) Recompute invariants of the strain rate tensor
99      !!                 which are inputs of the ITD, store stress
100      !!                 for the next time step
101      !!              6) Control prints of residual (convergence)
102      !!                 and charge ellipse.
103      !!                 The user should make sure that the parameters
104      !!                 nevp, telast and creepl maintain stress state
105      !!                 on the charge ellipse for plastic flow
106      !!                 e.g. in the Canadian Archipelago
107      !!
108      !! References : Hunke and Dukowicz, JPO97
109      !!              Bouillon et al., Ocean Modelling 2009
110      !!              Vancoppenolle et al., Ocean Modelling 2008
111      !!-------------------------------------------------------------------
112      INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation
113      INTEGER, INTENT(in) ::   k_jpj   ! northern j-index for ice computation
114      !!
115      INTEGER ::   ji, jj   ! dummy loop indices
116      INTEGER ::   jter     ! local integers
117      CHARACTER (len=50) ::   charout
118      REAL(wp) ::   zt11, zt12, zt21, zt22, ztagnx, ztagny, delta                         !
119      REAL(wp) ::   za, zstms, zmask   ! local scalars
120      REAL(wp) ::   zc1, zc2, zc3             ! ice mass
121
122      REAL(wp) ::   dtevp              ! time step for subcycling
123      REAL(wp) ::   dtotel, ecc2, ecci ! square of yield ellipse eccenticity
124      REAL(wp) ::   z0, zr, zcca, zccb ! temporary scalars
125      REAL(wp) ::   zu_ice2, zv_ice1   !
126      REAL(wp) ::   zddc, zdtc         ! delta on corners and on centre
127      REAL(wp) ::   zdst               ! shear at the center of the grid point
128      REAL(wp) ::   zdsshx, zdsshy     ! term for the gradient of ocean surface
129      REAL(wp) ::   sigma1, sigma2     ! internal ice stress
130
131      REAL(wp) ::   zresm         ! Maximal error on ice velocity
132      REAL(wp) ::   zdummy        ! dummy argument
133      REAL(wp) ::   zintb, zintn  ! dummy argument
134
135      REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh           ! temporary array for ice strength
136      REAL(wp), POINTER, DIMENSION(:,:) ::   zpreshc          ! Ice strength on grid cell corners (zpreshc)
137      REAL(wp), POINTER, DIMENSION(:,:) ::   zfrld1, zfrld2   ! lead fraction on U/V points
138      REAL(wp), POINTER, DIMENSION(:,:) ::   zmass1, zmass2   ! ice/snow mass on U/V points
139      REAL(wp), POINTER, DIMENSION(:,:) ::   zcorl1, zcorl2   ! coriolis parameter on U/V points
140      REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct , za2ct    ! temporary arrays
141      REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce1, v_oce1   ! ocean u/v component on U points                           
142      REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2, v_oce2   ! ocean u/v component on V points
143      REAL(wp), POINTER, DIMENSION(:,:) ::   u_ice2, v_ice1   ! ice u/v component on V/U point
144      REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2      ! arrays for internal stresses
145     
146      REAL(wp), POINTER, DIMENSION(:,:) ::   zdt              ! tension at centre of grid cells
147      REAL(wp), POINTER, DIMENSION(:,:) ::   zds              ! Shear on northeast corner of grid cells
148      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1   , zs2      ! Diagonal stress tensor components zs1 and zs2
149      REAL(wp), POINTER, DIMENSION(:,:) ::   zs12             ! Non-diagonal stress tensor component zs12
150      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr   ! Local error on velocity
151      REAL(wp), POINTER, DIMENSION(:,:) ::   zpice            ! array used for the calculation of ice surface slope:
152                                                              !   ocean surface (ssh_m) if ice is not embedded
153                                                              !   ice top surface if ice is embedded   
154      !!-------------------------------------------------------------------
155
156      CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
157      CALL wrk_alloc( jpi,jpj, u_oce1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1                )
158      CALL wrk_alloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  )
159      CALL wrk_alloc( jpi,jpj, zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )
160
161#if  defined key_lim2 && ! defined key_lim2_vp
162# if defined key_agrif
163     USE ice_2, vt_s => hsnm
164     USE ice_2, vt_i => hicm
165# else
166     vt_s => hsnm
167     vt_i => hicm
168# endif
169     at_i(:,:) = 1. - frld(:,:)
170#endif
171#if defined key_agrif && defined key_lim2 
172    CALL agrif_rhg_lim2_load      ! First interpolation of coarse values
173#endif
174      !
175      !------------------------------------------------------------------------------!
176      ! 1) Ice strength (zpresh)                                !
177      !------------------------------------------------------------------------------!
178      !
179      ! Put every vector to 0
180      delta_i(:,:) = 0._wp   ;
181      zpresh (:,:) = 0._wp   ; 
182      zpreshc(:,:) = 0._wp
183      u_ice2 (:,:) = 0._wp   ;   v_ice1(:,:) = 0._wp
184      divu_i (:,:) = 0._wp   ;   zdt   (:,:) = 0._wp   ;   zds(:,:) = 0._wp
185      shear_i(:,:) = 0._wp
186
187#if defined key_lim3
188      CALL lim_itd_me_icestrength( ridge_scheme_swi )      ! LIM-3: Ice strength on T-points
189#endif
190
191!CDIR NOVERRCHK
192      DO jj = k_j1 , k_jpj       ! Ice mass and temp variables
193!CDIR NOVERRCHK
194         DO ji = 1 , jpi
195#if defined key_lim3
196            zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj)
197#endif
198#if defined key_lim2
199            zpresh(ji,jj) = tms(ji,jj) *  pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) )
200#endif
201            ! tmi = 1 where there is ice or on land
202            tmi(ji,jj)    = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - epsd ) ) ) * tms(ji,jj)
203         END DO
204      END DO
205
206      ! Ice strength on grid cell corners (zpreshc)
207      ! needed for calculation of shear stress
208!CDIR NOVERRCHK
209      DO jj = k_j1+1, k_jpj-1
210!CDIR NOVERRCHK
211         DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1)
212            zstms          =  tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
213               &              tms(ji,jj+1)   * wght(ji+1,jj+1,1,2) + &
214               &              tms(ji+1,jj)   * wght(ji+1,jj+1,2,1) + &
215               &              tms(ji,jj)     * wght(ji+1,jj+1,1,1)
216            zpreshc(ji,jj) = (  zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
217               &                zpresh(ji,jj+1)   * wght(ji+1,jj+1,1,2) + &
218               &                zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + & 
219               &                zpresh(ji,jj)     * wght(ji+1,jj+1,1,1)   &
220               &             ) / MAX( zstms, epsd )
221         END DO
222      END DO
223
224      CALL lbc_lnk( zpreshc(:,:), 'F', 1. )
225      !
226      !------------------------------------------------------------------------------!
227      ! 2) Wind / ocean stress, mass terms, coriolis terms
228      !------------------------------------------------------------------------------!
229      !
230      !  Wind stress, coriolis and mass terms on the sides of the squares       
231      !  zfrld1: lead fraction on U-points                                     
232      !  zfrld2: lead fraction on V-points                                     
233      !  zmass1: ice/snow mass on U-points                                   
234      !  zmass2: ice/snow mass on V-points                                   
235      !  zcorl1: Coriolis parameter on U-points                             
236      !  zcorl2: Coriolis parameter on V-points                           
237      !  (ztagnx,ztagny): wind stress on U/V points                       
238      !  u_oce1: ocean u component on u points                           
239      !  v_oce1: ocean v component on u points                         
240      !  u_oce2: ocean u component on v points                         
241      !  v_oce2: ocean v component on v points                       
242
243      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==!
244          !                                           
245          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
246          !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
247         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp     
248          !
249          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
250          !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
251         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
252          !
253         zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0
254          !
255      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
256         zpice(:,:) = ssh_m(:,:)
257      ENDIF
258
259      DO jj = k_j1+1, k_jpj-1
260         DO ji = fs_2, fs_jpim1
261
262            zc1 = tms(ji  ,jj  ) * ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
263            zc2 = tms(ji+1,jj  ) * ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
264            zc3 = tms(ji  ,jj+1) * ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
265
266            zt11 = tms(ji  ,jj) * e1t(ji  ,jj)
267            zt12 = tms(ji+1,jj) * e1t(ji+1,jj)
268            zt21 = tms(ji,jj  ) * e2t(ji,jj  )
269            zt22 = tms(ji,jj+1) * e2t(ji,jj+1)
270
271            ! Leads area.
272            zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd )
273            zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd )
274
275            ! Mass, coriolis coeff. and currents
276            zmass1(ji,jj) = ( zt12*zc1 + zt11*zc2 ) / (zt11+zt12+epsd)
277            zmass2(ji,jj) = ( zt22*zc1 + zt21*zc3 ) / (zt21+zt22+epsd)
278            zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + e1t(ji,jj)*fcor(ji+1,jj) )   &
279               &                          / ( e1t(ji,jj) + e1t(ji+1,jj) + epsd )
280            zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)*fcor(ji,jj) + e2t(ji,jj)*fcor(ji,jj+1) )   &
281               &                          / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd )
282            !
283            u_oce1(ji,jj)  = u_oce(ji,jj)
284            v_oce2(ji,jj)  = v_oce(ji,jj)
285
286            ! Ocean has no slip boundary condition
287            v_oce1(ji,jj)  = 0.5*( (v_oce(ji,jj)+v_oce(ji,jj-1))*e1t(ji,jj)    &
288               &                 +(v_oce(ji+1,jj)+v_oce(ji+1,jj-1))*e1t(ji+1,jj)) &
289               &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
290
291            u_oce2(ji,jj)  = 0.5*((u_oce(ji,jj)+u_oce(ji-1,jj))*e2t(ji,jj)     &
292               &                 +(u_oce(ji,jj+1)+u_oce(ji-1,jj+1))*e2t(ji,jj+1)) &
293               &                / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
294
295            ! Wind stress at U,V-point
296            ztagnx = ( 1. - zfrld1(ji,jj) ) * utau_ice(ji,jj)
297            ztagny = ( 1. - zfrld2(ji,jj) ) * vtau_ice(ji,jj)
298
299            ! Computation of the velocity field taking into account the ice internal interaction.
300            ! Terms that are independent of the velocity field.
301
302            ! SB On utilise maintenant le gradient de la pente de l'ocean
303            ! include it later
304
305            zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj)
306            zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj)
307
308            za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx
309            za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy
310
311         END DO
312      END DO
313
314      !
315      !------------------------------------------------------------------------------!
316      ! 3) Solution of the momentum equation, iterative procedure
317      !------------------------------------------------------------------------------!
318      !
319      ! Time step for subcycling
320      dtevp  = rdt_ice / nevp
321      dtotel = dtevp / ( 2._wp * telast )
322
323      !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)
324      ecc2 = ecc * ecc
325      ecci = 1. / ecc2
326
327      !-Initialise stress tensor
328      zs1 (:,:) = stress1_i (:,:) 
329      zs2 (:,:) = stress2_i (:,:)
330      zs12(:,:) = stress12_i(:,:)
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               !- divu_i(:,:), 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 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               divu_i(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( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 
417               delta_i(ji,jj) = delta + creepl
418               !-Calculate stress tensor components zs1 and zs2
419               !-at centre of grid cells (see section 3.5 of CICE user's guide).
420               zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( divu_i(ji,jj) / delta_i(ji,jj) - delta / delta_i(ji,jj) )  &
421                  &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel )
422               zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) ) )  &
423                  &         / ( 1._wp + dtotel )
424
425            END DO
426         END DO
427
428         CALL lbc_lnk( zs1(:,:), 'T', 1. )
429         CALL lbc_lnk( zs2(:,:), 'T', 1. )
430
431!CDIR NOVERRCHK
432         DO jj = k_j1+1, k_jpj-1
433!CDIR NOVERRCHK
434            DO ji = fs_2, fs_jpim1
435               !- Calculate Delta on corners
436               zddc  =      ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1)                &
437                  &            -v_ice1(ji,jj)/e1u(ji,jj)                    &
438                  &           )*e1f(ji,jj)*e1f(ji,jj)                       &
439                  &          +( u_ice2(ji+1,jj)/e2v(ji+1,jj)                &
440                  &            -u_ice2(ji,jj)/e2v(ji,jj)                    &
441                  &           )*e2f(ji,jj)*e2f(ji,jj)                       &
442                  &         )                                               &
443                  &        / ( e1f(ji,jj) * e2f(ji,jj) )
444
445               zdtc  =      (-( v_ice1(ji,jj+1)/e1u(ji,jj+1)                &
446                  &            -v_ice1(ji,jj)/e1u(ji,jj)                    &
447                  &           )*e1f(ji,jj)*e1f(ji,jj)                       &
448                  &          +( u_ice2(ji+1,jj)/e2v(ji+1,jj)                &
449                  &            -u_ice2(ji,jj)/e2v(ji,jj)                    &
450                  &           )*e2f(ji,jj)*e2f(ji,jj)                       &
451                  &         )                                               &
452                  &        / ( e1f(ji,jj) * e2f(ji,jj) )
453
454               zddc = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl
455
456               !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide).
457               zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  &
458                  &          ( ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) ) )  &
459                  &          / ( 1.0 + dtotel ) 
460
461            END DO ! ji
462         END DO ! jj
463
464         CALL lbc_lnk( zs12(:,:), 'F', 1. )
465
466         ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002)
467         DO jj = k_j1+1, k_jpj-1
468            DO ji = fs_2, fs_jpim1
469               !- contribution of zs1, zs2 and zs12 to zf1
470               zf1(ji,jj) = 0.5*( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj) &
471                  &              +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj) &
472                  &              +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj) &
473                  &             ) / ( e1u(ji,jj)*e2u(ji,jj) )
474               ! contribution of zs1, zs2 and zs12 to zf2
475               zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) &
476                  &              -(zs2(ji,jj+1)*e1t(ji,jj+1)**2 - zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) &
477                  &              + 2.0*(zs12(ji,jj)*e2f(ji,jj)**2 -    &
478                  zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) &
479                  &             ) / ( e1v(ji,jj)*e2v(ji,jj) )
480            END DO
481         END DO
482         !
483         ! Computation of ice velocity
484         !
485         ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly.
486         !
487         IF (MOD(jter,2).eq.0) THEN 
488
489!CDIR NOVERRCHK
490            DO jj = k_j1+1, k_jpj-1
491!CDIR NOVERRCHK
492               DO ji = fs_2, fs_jpim1
493                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj)
494                  z0           = zmass1(ji,jj)/dtevp
495
496                  ! SB modif because ocean has no slip boundary condition
497                  zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)         &
498                     &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   &
499                     &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
500                  za           = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + &
501                     (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj))
502                  zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
503                     za*(u_oce1(ji,jj))
504                  zcca         = z0+za
505                  zccb         = zcorl1(ji,jj)
506                  u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 
507
508               END DO
509            END DO
510
511            CALL lbc_lnk( u_ice(:,:), 'U', -1. )
512#if defined key_agrif && defined key_lim2
513            CALL agrif_rhg_lim2( jter, nevp, 'U' )
514#endif
515#if defined key_bdy
516         CALL bdy_ice_lim_dyn( 'U' )
517#endif         
518
519!CDIR NOVERRCHK
520            DO jj = k_j1+1, k_jpj-1
521!CDIR NOVERRCHK
522               DO ji = fs_2, fs_jpim1
523
524                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj)
525                  z0           = zmass2(ji,jj)/dtevp
526                  ! SB modif because ocean has no slip boundary condition
527                  zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj)     &
528                     &                 + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1))   &
529                     &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
530                  za           = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 
531                     (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj))
532                  zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + &
533                     za2ct(ji,jj) + za*(v_oce2(ji,jj))
534                  zcca         = z0+za
535                  zccb         = zcorl2(ji,jj)
536                  v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask
537
538               END DO
539            END DO
540
541            CALL lbc_lnk( v_ice(:,:), 'V', -1. )
542#if defined key_agrif && defined key_lim2
543            CALL agrif_rhg_lim2( jter, nevp, 'V' )
544#endif
545#if defined key_bdy
546         CALL bdy_ice_lim_dyn( 'V' )
547#endif         
548
549         ELSE 
550!CDIR NOVERRCHK
551            DO jj = k_j1+1, k_jpj-1
552!CDIR NOVERRCHK
553               DO ji = fs_2, fs_jpim1
554                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj)
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*(v_oce2(ji,jj))
565                  zcca         = z0+za
566                  zccb         = zcorl2(ji,jj)
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#if defined key_agrif && defined key_lim2
574            CALL agrif_rhg_lim2( jter, nevp, 'V' )
575#endif
576#if defined key_bdy
577         CALL bdy_ice_lim_dyn( 'V' )
578#endif         
579
580!CDIR NOVERRCHK
581            DO jj = k_j1+1, k_jpj-1
582!CDIR NOVERRCHK
583               DO ji = fs_2, fs_jpim1
584                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj)
585                  z0           = zmass1(ji,jj)/dtevp
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*(u_oce1(ji,jj))
594                  zcca         = z0+za
595                  zccb         = zcorl1(ji,jj)
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#if defined key_agrif && defined key_lim2
602            CALL agrif_rhg_lim2( jter, nevp, 'U' )
603#endif
604#if defined key_bdy
605         CALL bdy_ice_lim_dyn( 'U' )
606#endif         
607
608         ENDIF
609         
610         IF(ln_ctl) THEN
611            !---  Convergence test.
612            DO jj = k_j1+1 , k_jpj-1
613               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) ,           &
614                  ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
615            END DO
616            zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) )
617            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
618         ENDIF
619
620         !                                                ! ==================== !
621      END DO                                              !  end loop over jter  !
622      !                                                   ! ==================== !
623      !
624      !------------------------------------------------------------------------------!
625      ! 4) Prevent ice velocities when the ice is thin
626      !------------------------------------------------------------------------------!
627      ! If the ice thickness is below hminrhg (5cm) then ice velocity should equal the
628      ! ocean velocity,
629      ! This prevents high velocity when ice is thin
630!CDIR NOVERRCHK
631      DO jj = k_j1+1, k_jpj-1
632!CDIR NOVERRCHK
633         DO ji = fs_2, fs_jpim1
634            zdummy = vt_i(ji,jj)
635            IF ( zdummy .LE. hminrhg ) THEN
636               u_ice(ji,jj) = u_oce(ji,jj)
637               v_ice(ji,jj) = v_oce(ji,jj)
638            ENDIF ! zdummy
639         END DO
640      END DO
641
642      CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
643      CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
644#if defined key_agrif && defined key_lim2
645      CALL agrif_rhg_lim2( nevp , nevp, 'U' )
646      CALL agrif_rhg_lim2( nevp , nevp, 'V' )
647#endif
648#if defined key_bdy
649      CALL bdy_ice_lim_dyn( 'U' )
650      CALL bdy_ice_lim_dyn( 'V' )
651#endif         
652
653      DO jj = k_j1+1, k_jpj-1 
654         DO ji = fs_2, fs_jpim1
655            zdummy = vt_i(ji,jj)
656            IF ( zdummy .LE. hminrhg ) THEN
657               v_ice1(ji,jj)  = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)   &
658                  &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
659                  &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
660
661               u_ice2(ji,jj)  = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)   &
662                  &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &
663                  &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
664            ENDIF ! zdummy
665         END DO
666      END DO
667
668      CALL lbc_lnk( u_ice2(:,:), 'V', -1. ) 
669      CALL lbc_lnk( v_ice1(:,:), 'U', -1. )
670
671      ! Recompute delta, shear and div, inputs for mechanical redistribution
672!CDIR NOVERRCHK
673      DO jj = k_j1+1, k_jpj-1
674!CDIR NOVERRCHK
675         DO ji = fs_2, jpim1   !RB bug no vect opt due to tmi
676            !- divu_i(:,:), zdt(:,:): divergence and tension at centre
677            !- zds(:,:): shear on northeast corner of grid cells
678            zdummy = vt_i(ji,jj)
679            IF ( zdummy .LE. hminrhg ) THEN
680
681               divu_i(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)  ) / area(ji,jj)
712
713               delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 
714               delta_i(ji,jj) = delta + creepl
715           
716            ENDIF ! zdummy
717
718         END DO !jj
719      END DO !ji
720      !
721      !------------------------------------------------------------------------------!
722      ! 5) Store stress tensor and its invariants
723      !------------------------------------------------------------------------------!
724      !
725      ! * Invariants of the stress tensor are required for limitd_me
726      !   (accelerates convergence and improves stability)
727      DO jj = k_j1+1, k_jpj-1
728         DO ji = fs_2, fs_jpim1
729            ! begin TECLIM change
730            zdst= (  e2u( ji  , jj   ) * v_ice1(ji,jj)           &   
731               &          - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)         &   
732               &          + e1v( ji  , jj   ) * u_ice2(ji,jj)           &   
733               &          - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj) 
734            shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst )
735            ! end TECLIM change
736         END DO
737      END DO
738
739      ! Lateral boundary condition
740      CALL lbc_lnk( divu_i (:,:), 'T', 1. )
741      CALL lbc_lnk( delta_i(:,:), 'T', 1. )
742      ! CALL lbc_lnk( shear_i(:,:), 'F', 1. )
743      CALL lbc_lnk( shear_i(:,:), 'T', 1. )
744
745      ! * Store the stress tensor for the next time step
746      stress1_i (:,:) = zs1 (:,:)
747      stress2_i (:,:) = zs2 (:,:)
748      stress12_i(:,:) = zs12(:,:)
749
750      !
751      !------------------------------------------------------------------------------!
752      ! 6) Control prints of residual and charge ellipse
753      !------------------------------------------------------------------------------!
754      !
755      ! print the residual for convergence
756      IF(ln_ctl) THEN
757         WRITE(charout,FMT="('lim_rhg  : res =',D23.16, ' iter =',I4)") zresm, jter
758         CALL prt_ctl_info(charout)
759         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
760      ENDIF
761
762      ! print charge ellipse
763      ! This can be desactivated once the user is sure that the stress state
764      ! lie on the charge ellipse. See Bouillon et al. 08 for more details
765      IF(ln_ctl) THEN
766         CALL prt_ctl_info('lim_rhg  : numit  :',ivar1=numit)
767         CALL prt_ctl_info('lim_rhg  : nwrite :',ivar1=nwrite)
768         CALL prt_ctl_info('lim_rhg  : MOD    :',ivar1=MOD(numit,nwrite))
769         IF( MOD(numit,nwrite) .EQ. 0 ) THEN
770            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. '
771            CALL prt_ctl_info(charout)
772            DO jj = k_j1+1, k_jpj-1
773               DO ji = 2, jpim1
774                  IF (zpresh(ji,jj) .GT. 1.0) THEN
775                     sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
776                     sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )
777                     WRITE(charout,FMT="('lim_rhg  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
778                     CALL prt_ctl_info(charout)
779                  ENDIF
780               END DO
781            END DO
782            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 2000, numit, 0, 0, ' ch. ell. '
783            CALL prt_ctl_info(charout)
784         ENDIF
785      ENDIF
786      !
787      CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
788      CALL wrk_dealloc( jpi,jpj, u_oce1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1                )
789      CALL wrk_dealloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  )
790      CALL wrk_dealloc( jpi,jpj, zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )
791
792   END SUBROUTINE lim_rhg
793
794#else
795   !!----------------------------------------------------------------------
796   !!   Default option          Dummy module           NO LIM sea-ice model
797   !!----------------------------------------------------------------------
798CONTAINS
799   SUBROUTINE lim_rhg( k1 , k2 )         ! Dummy routine
800      WRITE(*,*) 'lim_rhg: You should not have seen this print! error?', k1, k2
801   END SUBROUTINE lim_rhg
802#endif
803
804   !!==============================================================================
805END MODULE limrhg
Note: See TracBrowser for help on using the repository browser.