source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 5070

Last change on this file since 5070 was 5070, checked in by clem, 6 years ago

LIM3 cosmetic changes

  • Property svn:keywords set to Id
File size: 36.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   !!            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, elastic time scale 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      !!-------------------------------------------------------------------
111      INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation
112      INTEGER, INTENT(in) ::   k_jpj   ! northern j-index for ice computation
113      !!
114      INTEGER ::   ji, jj   ! dummy loop indices
115      INTEGER ::   jter     ! local integers
116      CHARACTER (len=50) ::   charout
117      REAL(wp) ::   zt11, zt12, zt21, zt22, ztagnx, ztagny, delta                         !
118      REAL(wp) ::   za, zstms          ! local scalars
119      REAL(wp) ::   zc1, zc2, zc3      ! ice mass
120
121      REAL(wp) ::   dtevp , z1_dtevp              ! time step for subcycling
122      REAL(wp) ::   dtotel, z1_dtotel, ecc2, ecci ! square of yield ellipse eccenticity
123      REAL(wp) ::   z0, zr, zcca, zccb            ! temporary scalars
124      REAL(wp) ::   zu_ice2, zv_ice1              !
125      REAL(wp) ::   zddc, zdtc                    ! delta on corners and on centre
126      REAL(wp) ::   zdst                          ! shear at the center of the grid point
127      REAL(wp) ::   zdsshx, zdsshy                ! term for the gradient of ocean surface
128      REAL(wp) ::   sigma1, sigma2                ! internal ice stress
129
130      REAL(wp) ::   zresm         ! Maximal error on ice velocity
131      REAL(wp) ::   zintb, zintn  ! dummy argument
132
133      REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh           ! temporary array for ice strength
134      REAL(wp), POINTER, DIMENSION(:,:) ::   zpreshc          ! Ice strength on grid cell corners (zpreshc)
135      REAL(wp), POINTER, DIMENSION(:,:) ::   zfrld1, zfrld2   ! lead fraction on U/V points
136      REAL(wp), POINTER, DIMENSION(:,:) ::   zmass1, zmass2   ! ice/snow mass on U/V points
137      REAL(wp), POINTER, DIMENSION(:,:) ::   zcorl1, zcorl2   ! coriolis parameter on U/V points
138      REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct , za2ct    ! temporary arrays
139      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oce1           ! ocean u/v component on U points                           
140      REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2           ! ocean u/v component on V points
141      REAL(wp), POINTER, DIMENSION(:,:) ::   u_ice2, v_ice1   ! ice u/v component on V/U point
142      REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2      ! arrays for internal stresses
143      REAL(wp), POINTER, DIMENSION(:,:) ::   zmask            ! mask ocean grid points
144     
145      REAL(wp), POINTER, DIMENSION(:,:) ::   zdt              ! tension at centre of grid cells
146      REAL(wp), POINTER, DIMENSION(:,:) ::   zds              ! Shear on northeast corner of grid cells
147      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1   , zs2      ! Diagonal stress tensor components zs1 and zs2
148      REAL(wp), POINTER, DIMENSION(:,:) ::   zs12             ! Non-diagonal stress tensor component zs12
149      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr   ! Local error on velocity
150      REAL(wp), POINTER, DIMENSION(:,:) ::   zpice            ! array used for the calculation of ice surface slope:
151                                                              !   ocean surface (ssh_m) if ice is not embedded
152                                                              !   ice top surface if ice is embedded   
153
154      REAL(wp), PARAMETER               ::   zepsi = 1.0e-20_wp ! tolerance parameter
155      REAL(wp), PARAMETER               ::   zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity
156      !!-------------------------------------------------------------------
157
158      CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
159      CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               )
160      CALL wrk_alloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  )
161      CALL wrk_alloc( jpi,jpj, zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )
162
163#if  defined key_lim2 && ! defined key_lim2_vp
164# if defined key_agrif
165      USE ice_2, vt_s => hsnm
166      USE ice_2, vt_i => hicm
167# else
168      vt_s => hsnm
169      vt_i => hicm
170# endif
171      at_i(:,:) = 1. - frld(:,:)
172#endif
173#if defined key_agrif && defined key_lim2 
174      CALL agrif_rhg_lim2_load      ! First interpolation of coarse values
175#endif
176      !
177      !------------------------------------------------------------------------------!
178      ! 1) Ice strength (zpresh)                                !
179      !------------------------------------------------------------------------------!
180      !
181      ! Put every vector to 0
182      delta_i(:,:) = 0._wp   ;
183      zpresh (:,:) = 0._wp   ; 
184      zpreshc(:,:) = 0._wp
185      u_ice2 (:,:) = 0._wp   ;   v_ice1(:,:) = 0._wp
186      divu_i (:,:) = 0._wp   ;   zdt   (:,:) = 0._wp   ;   zds(:,:) = 0._wp
187      shear_i(:,:) = 0._wp
188
189#if defined key_lim3
190      CALL lim_itd_me_icestrength( nn_icestr )      ! LIM-3: Ice strength on T-points
191#endif
192
193      DO jj = k_j1 , k_jpj       ! Ice mass and temp variables
194         DO ji = 1 , jpi
195#if defined key_lim3
196            zpresh(ji,jj) = tmask(ji,jj,1) *  strength(ji,jj)
197#endif
198#if defined key_lim2
199            zpresh(ji,jj) = tmask(ji,jj,1) *  pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) )
200#endif
201            ! zmask = 1 where there is ice or on land
202            zmask(ji,jj)    = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tmask(ji,jj,1)
203         END DO
204      END DO
205
206      ! Ice strength on grid cell corners (zpreshc)
207      ! needed for calculation of shear stress
208      DO jj = k_j1+1, k_jpj-1
209         DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1)
210            zstms          =  tmask(ji+1,jj+1,1) * wght(ji+1,jj+1,2,2) + tmask(ji,jj+1,1) * wght(ji+1,jj+1,1,2) +   &
211               &              tmask(ji+1,jj,1)   * wght(ji+1,jj+1,2,1) + tmask(ji,jj,1)   * wght(ji+1,jj+1,1,1)
212            zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) +   &
213               &               zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + zpresh(ji,jj)   * wght(ji+1,jj+1,1,1)     &
214               &             ) / MAX( zstms, zepsi )
215         END DO
216      END DO
217      CALL lbc_lnk( zpreshc(:,:), 'F', 1. )
218      !
219      !------------------------------------------------------------------------------!
220      ! 2) Wind / ocean stress, mass terms, coriolis terms
221      !------------------------------------------------------------------------------!
222      !
223      !  Wind stress, coriolis and mass terms on the sides of the squares       
224      !  zfrld1: lead fraction on U-points                                     
225      !  zfrld2: lead fraction on V-points                                     
226      !  zmass1: ice/snow mass on U-points                                   
227      !  zmass2: ice/snow mass on V-points                                   
228      !  zcorl1: Coriolis parameter on U-points                             
229      !  zcorl2: Coriolis parameter on V-points                           
230      !  (ztagnx,ztagny): wind stress on U/V points                       
231      !  v_oce1: ocean v component on u points                         
232      !  u_oce2: ocean u component on v points                         
233
234      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==!
235         !                                           
236         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
237         !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
238         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp     
239         !
240         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
241         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
242         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
243         !
244         zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0
245         !
246      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
247         zpice(:,:) = ssh_m(:,:)
248      ENDIF
249
250      DO jj = k_j1+1, k_jpj-1
251         DO ji = fs_2, fs_jpim1
252
253            zc1 = tmask(ji  ,jj  ,1) * ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
254            zc2 = tmask(ji+1,jj  ,1) * ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
255            zc3 = tmask(ji  ,jj+1,1) * ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
256
257            zt11 = tmask(ji  ,jj,1) * e1t(ji  ,jj)
258            zt12 = tmask(ji+1,jj,1) * e1t(ji+1,jj)
259            zt21 = tmask(ji,jj  ,1) * e2t(ji,jj  )
260            zt22 = tmask(ji,jj+1,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 + zepsi )
264            zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + zepsi )
265
266            ! Mass, coriolis coeff. and currents
267            zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi )
268            zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi )
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) + zepsi )
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) + zepsi )
273            !
274            ! Ocean has no slip boundary condition
275            v_oce1(ji,jj)  = 0.5 * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji,jj)      &
276               &                   + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) )  &
277               &                   / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
278
279            u_oce2(ji,jj)  = 0.5 * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj)      &
280               &                   + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) )  &
281               &                   / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
282
283            ! Wind stress at U,V-point
284            ztagnx = ( 1. - zfrld1(ji,jj) ) * utau_ice(ji,jj)
285            ztagny = ( 1. - zfrld2(ji,jj) ) * vtau_ice(ji,jj)
286
287            ! Computation of the velocity field taking into account the ice internal interaction.
288            ! Terms that are independent of the velocity field.
289
290            ! SB On utilise maintenant le gradient de la pente de l'ocean
291            ! include it later
292
293            zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
294            zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
295
296            za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx
297            za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy
298
299         END DO
300      END DO
301
302      !
303      !------------------------------------------------------------------------------!
304      ! 3) Solution of the momentum equation, iterative procedure
305      !------------------------------------------------------------------------------!
306      !
307      ! Time step for subcycling
308      dtevp  = rdt_ice / nn_nevp
309#if defined key_lim3
310      dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice )
311#else
312      dtotel = dtevp / ( 2._wp * telast )
313#endif
314      z1_dtotel = 1._wp / ( 1._wp + dtotel )
315      z1_dtevp  = 1._wp / dtevp
316      !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)
317      ecc2 = rn_ecc * rn_ecc
318      ecci = 1. / ecc2
319
320      !-Initialise stress tensor
321      zs1 (:,:) = stress1_i (:,:) 
322      zs2 (:,:) = stress2_i (:,:)
323      zs12(:,:) = stress12_i(:,:)
324
325      !                                               !----------------------!
326      DO jter = 1 , nn_nevp                           !    loop over jter    !
327         !                                            !----------------------!       
328         DO jj = k_j1, k_jpj-1
329            zu_ice(:,jj) = u_ice(:,jj)    ! velocity at previous time step
330            zv_ice(:,jj) = v_ice(:,jj)
331         END DO
332
333         DO jj = k_j1+1, k_jpj-1
334            DO ji = fs_2, fs_jpim1   !RB bug no vect opt due to zmask
335
336               
337               !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002)
338               !- divu_i(:,:), zdt(:,:): divergence and tension at centre of grid cells
339               !- zds(:,:): shear on northeast corner of grid cells
340               !
341               !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded,
342               !                      there are many repeated calculations.
343               !                      Speed could be improved by regrouping terms. For
344               !                      the moment, however, the stress is on clarity of coding to avoid
345               !                      bugs (Martin, for Miguel).
346               !
347               !- ALSO: arrays zdt, zds and delta could
348               !  be removed in the future to minimise memory demand.
349               !
350               !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of
351               !              grid cells, exactly as in the B grid case. For simplicity, the indexation on
352               !              the corners is the same as in the B grid.
353               !
354               !
355               divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
356                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
357                  &            ) * r1_e12t(ji,jj)
358
359               zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
360                  &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
361                  &         ) * r1_e12t(ji,jj)
362
363               !
364               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
365                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
366                  &         ) * r1_e12f(ji,jj) * ( 2._wp - fmask(ji,jj,1) )   &
367                  &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1)
368
369
370               v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
371                  &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   &
372                  &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
373
374               u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
375                  &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   &
376                  &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
377            END DO
378         END DO
379         CALL lbc_lnk( v_ice1  , 'U', -1. )   ;   CALL lbc_lnk( u_ice2  , 'V', -1. )      ! lateral boundary cond.
380
381         DO jj = k_j1+1, k_jpj-1
382            DO ji = fs_2, fs_jpim1
383
384               !- Calculate Delta at centre of grid cells
385               zdst          = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )   &
386                  &            + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1)   &
387                  &            ) * r1_e12t(ji,jj)
388
389               delta          = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 
390               delta_i(ji,jj) = delta + rn_creepl
391
392               !- Calculate Delta on corners
393               zddc  = (  ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
394                  &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
395                  &    ) * r1_e12f(ji,jj)
396
397               zdtc  = (- ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
398                  &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
399                  &    ) * r1_e12f(ji,jj)
400
401               zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + rn_creepl
402
403               !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide).
404               zs1(ji,jj)  = ( zs1 (ji,jj) + dtotel * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj)   * zpresh(ji,jj)   &
405                  &          ) * z1_dtotel
406               zs2(ji,jj)  = ( zs2 (ji,jj) + dtotel *         ecci * zdt(ji,jj) / delta_i(ji,jj)   * zpresh(ji,jj)   &
407                  &          ) * z1_dtotel
408               !-Calculate stress tensor component zs12 at corners
409               zs12(ji,jj) = ( zs12(ji,jj) + dtotel *         ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj)  &
410                  &          ) * z1_dtotel 
411
412            END DO
413         END DO
414         CALL lbc_lnk( zs1 , 'T', 1. )   ;   CALL lbc_lnk( zs2, 'T', 1. )
415         CALL lbc_lnk( zs12, 'F', 1. )
416
417         ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002)
418         DO jj = k_j1+1, k_jpj-1
419            DO ji = fs_2, fs_jpim1
420               !- contribution of zs1, zs2 and zs12 to zf1
421               zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)  &
422                  &             + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) * r1_e2u(ji,jj)          &
423                  &             + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) * r1_e1u(ji,jj)  &
424                  &                ) * r1_e12u(ji,jj)
425               ! contribution of zs1, zs2 and zs12 to zf2
426               zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)  &
427                  &             - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) * r1_e1v(ji,jj)          &
428                  &             + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) * r1_e2v(ji,jj)  &
429                  &               )  * r1_e12v(ji,jj)
430            END DO
431         END DO
432         !
433         ! Computation of ice velocity
434         !
435         ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly.
436         !
437         IF (MOD(jter,2).eq.0) THEN
438
439            DO jj = k_j1+1, k_jpj-1
440               DO ji = fs_2, fs_jpim1
441                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1)
442                  z0           = zmass1(ji,jj) * z1_dtevp
443
444                  ! SB modif because ocean has no slip boundary condition
445                  zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji  ,jj)     &
446                     &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   &
447                     &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
448                  za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  &
449                     &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) )
450                  zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj)
451                  zcca         = z0 + za
452                  zccb         = zcorl1(ji,jj)
453                  u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 
454               END DO
455            END DO
456
457            CALL lbc_lnk( u_ice(:,:), 'U', -1. )
458#if defined key_agrif && defined key_lim2
459            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
460#endif
461#if defined key_bdy
462         CALL bdy_ice_lim_dyn( 'U' )
463#endif         
464
465            DO jj = k_j1+1, k_jpj-1
466               DO ji = fs_2, fs_jpim1
467
468                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1)
469                  z0           = zmass2(ji,jj) * z1_dtevp
470                  ! SB modif because ocean has no slip boundary condition
471                  zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       &
472                     &                 + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   &
473                     &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
474                  za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  & 
475                     &                         ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) )
476                  zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj)
477                  zcca         = z0 + za
478                  zccb         = zcorl2(ji,jj)
479                  v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch
480               END DO
481            END DO
482
483            CALL lbc_lnk( v_ice(:,:), 'V', -1. )
484#if defined key_agrif && defined key_lim2
485            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
486#endif
487#if defined key_bdy
488         CALL bdy_ice_lim_dyn( 'V' )
489#endif         
490
491         ELSE
492            DO jj = k_j1+1, k_jpj-1
493               DO ji = fs_2, fs_jpim1
494                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1)
495                  z0           = zmass2(ji,jj) * z1_dtevp
496                  ! SB modif because ocean has no slip boundary condition
497                  zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       &
498                     &                  +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   &
499                     &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)   
500
501                  za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  &
502                     &                         ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) )
503                  zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj)
504                  zcca         = z0 + za
505                  zccb         = zcorl2(ji,jj)
506                  v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch
507               END DO
508            END DO
509
510            CALL lbc_lnk( v_ice(:,:), 'V', -1. )
511#if defined key_agrif && defined key_lim2
512            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
513#endif
514#if defined key_bdy
515         CALL bdy_ice_lim_dyn( 'V' )
516#endif         
517
518            DO jj = k_j1+1, k_jpj-1
519               DO ji = fs_2, fs_jpim1
520                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1)
521                  z0           = zmass1(ji,jj) * z1_dtevp
522                  zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji,jj)       &
523                     &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   &
524                     &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
525
526                  za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  &
527                     &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) )
528                  zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj)
529                  zcca         = z0 + za
530                  zccb         = zcorl1(ji,jj)
531                  u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 
532               END DO
533            END DO
534
535            CALL lbc_lnk( u_ice(:,:), 'U', -1. )
536#if defined key_agrif && defined key_lim2
537            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
538#endif
539#if defined key_bdy
540         CALL bdy_ice_lim_dyn( 'U' )
541#endif         
542
543         ENDIF
544         
545         IF(ln_ctl) THEN
546            !---  Convergence test.
547            DO jj = k_j1+1 , k_jpj-1
548               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
549            END DO
550            zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) )
551            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
552         ENDIF
553
554         !                                                ! ==================== !
555      END DO                                              !  end loop over jter  !
556      !                                                   ! ==================== !
557      !
558      !------------------------------------------------------------------------------!
559      ! 4) Prevent ice velocities when the ice is thin
560      !------------------------------------------------------------------------------!
561      ! If the ice volume is below zvmin then ice velocity should equal the
562      ! ocean velocity. This prevents high velocity when ice is thin
563      DO jj = k_j1+1, k_jpj-1
564         DO ji = fs_2, fs_jpim1
565            IF ( vt_i(ji,jj) <= zvmin ) THEN
566               u_ice(ji,jj) = u_oce(ji,jj)
567               v_ice(ji,jj) = v_oce(ji,jj)
568            ENDIF
569         END DO
570      END DO
571
572      CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
573      CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
574#if defined key_agrif && defined key_lim2
575      CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' )
576      CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' )
577#endif
578#if defined key_bdy
579      CALL bdy_ice_lim_dyn( 'U' )
580      CALL bdy_ice_lim_dyn( 'V' )
581#endif         
582
583      DO jj = k_j1+1, k_jpj-1 
584         DO ji = fs_2, fs_jpim1
585            IF ( vt_i(ji,jj) <= zvmin ) THEN
586               v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji,  jj-1) ) * e1t(ji+1,jj)     &
587                  &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   &
588                  &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
589
590               u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
591                  &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   &
592                  &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
593            ENDIF
594         END DO
595      END DO
596
597      CALL lbc_lnk( u_ice2(:,:), 'V', -1. ) 
598      CALL lbc_lnk( v_ice1(:,:), 'U', -1. )
599
600      ! Recompute delta, shear and div, inputs for mechanical redistribution
601      DO jj = k_j1+1, k_jpj-1
602         DO ji = fs_2, jpim1   !RB bug no vect opt due to zmask
603            !- divu_i(:,:), zdt(:,:): divergence and tension at centre
604            !- zds(:,:): shear on northeast corner of grid cells
605            IF ( vt_i(ji,jj) <= zvmin ) THEN
606
607               divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj  ) * u_ice(ji-1,jj  )   &
608                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji  ,jj-1) * v_ice(ji  ,jj-1)   &
609                  &            ) * r1_e12t(ji,jj)
610
611               zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)  &
612                  &          -( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)  &
613                  &         ) * r1_e12t(ji,jj)
614               !
615               ! SB modif because ocean has no slip boundary condition
616               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
617                  &          +( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
618                  &         ) * r1_e12f(ji,jj) * ( 2.0 - fmask(ji,jj,1) )                                     &
619                  &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1)
620
621               zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )    &
622                  &   + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1) ) * r1_e12t(ji,jj)
623
624               delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 
625               delta_i(ji,jj) = delta + rn_creepl
626           
627            ENDIF
628         END DO
629      END DO
630      !
631      !------------------------------------------------------------------------------!
632      ! 5) Store stress tensor and its invariants
633      !------------------------------------------------------------------------------!
634      ! * Invariants of the stress tensor are required for limitd_me
635      !   (accelerates convergence and improves stability)
636      DO jj = k_j1+1, k_jpj-1
637         DO ji = fs_2, fs_jpim1
638            zdst           = (  e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)  &   
639               &              + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e12t(ji,jj) 
640            shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst )
641         END DO
642      END DO
643
644      ! Lateral boundary condition
645      CALL lbc_lnk( divu_i (:,:), 'T', 1. )
646      CALL lbc_lnk( delta_i(:,:), 'T', 1. )
647      ! CALL lbc_lnk( shear_i(:,:), 'F', 1. )
648      CALL lbc_lnk( shear_i(:,:), 'T', 1. )
649
650      ! * Store the stress tensor for the next time step
651      stress1_i (:,:) = zs1 (:,:)
652      stress2_i (:,:) = zs2 (:,:)
653      stress12_i(:,:) = zs12(:,:)
654
655      !
656      !------------------------------------------------------------------------------!
657      ! 6) Control prints of residual and charge ellipse
658      !------------------------------------------------------------------------------!
659      !
660      ! print the residual for convergence
661      IF(ln_ctl) THEN
662         WRITE(charout,FMT="('lim_rhg  : res =',D23.16, ' iter =',I4)") zresm, jter
663         CALL prt_ctl_info(charout)
664         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
665      ENDIF
666
667      ! print charge ellipse
668      ! This can be desactivated once the user is sure that the stress state
669      ! lie on the charge ellipse. See Bouillon et al. 08 for more details
670      IF(ln_ctl) THEN
671         CALL prt_ctl_info('lim_rhg  : numit  :',ivar1=numit)
672         CALL prt_ctl_info('lim_rhg  : nwrite :',ivar1=nwrite)
673         CALL prt_ctl_info('lim_rhg  : MOD    :',ivar1=MOD(numit,nwrite))
674         IF( MOD(numit,nwrite) .EQ. 0 ) THEN
675            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. '
676            CALL prt_ctl_info(charout)
677            DO jj = k_j1+1, k_jpj-1
678               DO ji = 2, jpim1
679                  IF (zpresh(ji,jj) > 1.0) THEN
680                     sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
681                     sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )
682                     WRITE(charout,FMT="('lim_rhg  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
683                     CALL prt_ctl_info(charout)
684                  ENDIF
685               END DO
686            END DO
687            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 2000, numit, 0, 0, ' ch. ell. '
688            CALL prt_ctl_info(charout)
689         ENDIF
690      ENDIF
691      !
692      CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
693      CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               )
694      CALL wrk_dealloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  )
695      CALL wrk_dealloc( jpi,jpj, zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )
696
697   END SUBROUTINE lim_rhg
698
699#else
700   !!----------------------------------------------------------------------
701   !!   Default option          Dummy module           NO LIM sea-ice model
702   !!----------------------------------------------------------------------
703CONTAINS
704   SUBROUTINE lim_rhg( k1 , k2 )         ! Dummy routine
705      WRITE(*,*) 'lim_rhg: You should not have seen this print! error?', k1, k2
706   END SUBROUTINE lim_rhg
707#endif
708
709   !!==============================================================================
710END MODULE limrhg
Note: See TracBrowser for help on using the repository browser.