source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_2/limrhg.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

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