New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limrhg.F90 in branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

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

Last change on this file since 5048 was 5047, checked in by clem, 9 years ago

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