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.
icedyn_rdgrft.F90 in NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE – NEMO

source: NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedyn_rdgrft.F90 @ 9923

Last change on this file since 9923 was 9923, checked in by gm, 6 years ago

#1911 (ENHANCE-04): step I.2: dev_r9838_ENHANCE04_MLF

File size: 50.6 KB
Line 
1MODULE icedyn_rdgrft
2   !!======================================================================
3   !!                       ***  MODULE icedyn_rdgrft ***
4   !!    sea-ice : Mechanical impact on ice thickness distribution     
5   !!======================================================================
6   !! History :       !  2006-02  (M. Vancoppenolle) Original code
7   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
8   !!----------------------------------------------------------------------
9#if defined key_si3
10   !!----------------------------------------------------------------------
11   !!   'key_si3'                                       SI3 sea-ice model
12   !!----------------------------------------------------------------------
13   !!   ice_dyn_rdgrft       : ridging/rafting of sea ice
14   !!   ice_dyn_rdgrft_init  : initialization of ridging/rafting of sea ice
15   !!   ice_strength         : ice strength calculation
16   !!----------------------------------------------------------------------
17   USE dom_oce        ! ocean domain
18   USE phycst         ! physical constants (ocean directory)
19   USE sbc_oce , ONLY : sss_m, sst_m   ! surface boundary condition: ocean fields
20   USE ice1D          ! sea-ice: thermodynamics
21   USE ice            ! sea-ice: variables
22   USE icetab         ! sea-ice: 1D <==> 2D transformation
23   USE icevar         ! sea-ice: operations
24   USE icectl         ! sea-ice: control prints
25   !
26   USE in_out_manager ! I/O manager
27   USE iom            ! I/O manager library
28   USE lib_mpp        ! MPP library
29   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
30   USE lbclnk         ! lateral boundary conditions (or mpp links)
31   USE timing         ! Timing
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   ice_dyn_rdgrft        ! called by icestp
37   PUBLIC   ice_dyn_rdgrft_init   ! called by icedyn
38   PUBLIC   ice_strength          ! called by icedyn_rhg_evp
39
40   ! Variables shared among ridging subroutines
41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   closing_net     ! net rate at which area is removed    (1/s)
42      !                                                 ! (ridging ice area - area of new ridges) / dt
43   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   opning         ! rate of opening due to divergence/shear
44   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   closing_gross  ! rate at which area removed, not counting area of new ridges
45   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   apartf   ! participation function; fraction of ridging/closing associated w/ category n
46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hrmin    ! minimum ridge thickness
47   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hrmax    ! maximum ridge thickness
48   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hraft    ! thickness of rafted ice
49   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hi_hrdg  ! thickness of ridging ice / mean ridge thickness
50   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   aridge   ! participating ice ridging
51   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   araft    ! participating ice rafting
52   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ze_i_2d
53   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ze_s_2d
54   !
55   REAL(wp), PARAMETER ::   hrdg_hi_min = 1.1_wp    ! min ridge thickness multiplier: min(hrdg/hi)
56   REAL(wp), PARAMETER ::   hi_hrft     = 0.5_wp    ! rafting multipliyer: (hi/hraft)
57   !
58   ! ** namelist (namdyn_rdgrft) **
59   LOGICAL  ::   ln_str_H79       ! ice strength parameterization (Hibler79)
60   REAL(wp) ::   rn_pstar         ! determines ice strength, Hibler JPO79
61   REAL(wp) ::   rn_crhg          ! determines changes in ice strength
62   REAL(wp) ::   rn_csrdg         ! fraction of shearing energy contributing to ridging           
63   LOGICAL  ::   ln_partf_lin     ! participation function linear (Thorndike et al. (1975))
64   REAL(wp) ::   rn_gstar         !    fractional area of young ice contributing to ridging
65   LOGICAL  ::   ln_partf_exp     ! participation function exponential (Lipscomb et al. (2007))
66   REAL(wp) ::   rn_astar         !    equivalent of G* for an exponential participation function
67   LOGICAL  ::   ln_ridging       ! ridging of ice or not                       
68   REAL(wp) ::   rn_hstar         !    thickness that determines the maximal thickness of ridged ice
69   REAL(wp) ::   rn_porordg       !    initial porosity of ridges (0.3 regular value)
70   REAL(wp) ::   rn_fsnwrdg       !    fractional snow loss to the ocean during ridging
71   REAL(wp) ::   rn_fpndrdg       !    fractional pond loss to the ocean during ridging
72   LOGICAL  ::   ln_rafting       ! rafting of ice or not                       
73   REAL(wp) ::   rn_hraft         !    threshold thickness (m) for rafting / ridging
74   REAL(wp) ::   rn_craft         !    coefficient for smoothness of the hyperbolic tangent in rafting
75   REAL(wp) ::   rn_fsnwrft       !    fractional snow loss to the ocean during rafting
76   REAL(wp) ::   rn_fpndrft       !    fractional pond loss to the ocean during rafting
77   !
78   !!----------------------------------------------------------------------
79   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
80   !! $Id: icedyn_rdgrft.F90 8378 2017-07-26 13:55:59Z clem $
81   !! Software governed by the CeCILL licence     (./LICENSE)
82   !!----------------------------------------------------------------------
83CONTAINS
84
85   INTEGER FUNCTION ice_dyn_rdgrft_alloc()
86      !!-------------------------------------------------------------------
87      !!                ***  ROUTINE ice_dyn_rdgrft_alloc ***
88      !!-------------------------------------------------------------------
89      ALLOCATE( closing_net(jpij), opning(jpij)   , closing_gross(jpij),   &
90         &      apartf(jpij,0:jpl), hrmin(jpij,jpl), hraft(jpij,jpl)    , aridge(jpij,jpl),   &
91         &      hrmax(jpij,jpl), hi_hrdg(jpij,jpl)  , araft (jpij,jpl),  &
92         &      ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc )
93
94      IF( lk_mpp                    )   CALL mpp_sum ( ice_dyn_rdgrft_alloc )
95      IF( ice_dyn_rdgrft_alloc /= 0 )   CALL ctl_warn( 'ice_dyn_rdgrft_alloc: failed to allocate arrays' )
96      !
97   END FUNCTION ice_dyn_rdgrft_alloc
98
99
100   SUBROUTINE ice_dyn_rdgrft( kt )
101      !!-------------------------------------------------------------------
102      !!                ***  ROUTINE ice_dyn_rdgrft ***
103      !!
104      !! ** Purpose :   computes the mechanical redistribution of ice thickness
105      !!
106      !! ** Method  :   Steps :
107      !!       0) Identify grid cells with ice
108      !!       1) Calculate closing rate, divergence and opening
109      !!       2) Identify grid cells with ridging
110      !!       3) Start ridging iterations
111      !!          - prep = ridged and rafted ice + closing_gross
112      !!          - shift = move ice from one category to another
113      !!
114      !! ** Details
115      !!    step1: The net rate of closing is due to convergence and shear, based on Flato and Hibler (1995).
116      !!           The energy dissipation rate is equal to the net closing rate times the ice strength.
117      !!
118      !!    step3: The gross closing rate is equal to the first two terms (open
119      !!           water closing and thin ice ridging) without the third term
120      !!           (thick, newly ridged ice).
121      !!
122      !! References :   Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626.
123      !!                Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980.
124      !!                Rothrock, D. A., 1975: JGR, 80, 4514-4519.
125      !!                Thorndike et al., 1975, JGR, 80, 4501-4513.
126      !!                Bitz et al., JGR, 2001
127      !!                Amundrud and Melling, JGR 2005
128      !!                Babko et al., JGR 2002
129      !!
130      !!     This routine is based on CICE code and authors William H. Lipscomb,
131      !!     and Elizabeth C. Hunke, LANL are gratefully acknowledged
132      !!-------------------------------------------------------------------
133      INTEGER, INTENT(in) ::   kt     ! number of iteration
134      !!
135      INTEGER  ::   ji, jj, jk, jl             ! dummy loop index
136      INTEGER  ::   iter, iterate_ridging      ! local integer
137      INTEGER  ::   ipti                       ! local integer
138      REAL(wp) ::   zfac                       ! local scalar
139      INTEGER , DIMENSION(jpij) ::   iptidx        ! compute ridge/raft or not
140      REAL(wp), DIMENSION(jpij) ::   zdivu_adv     ! divu as implied by transport scheme  (1/s)
141      REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i
142      !
143      INTEGER, PARAMETER ::   jp_itermax = 20   
144      !!-------------------------------------------------------------------
145      ! controls
146      IF( ln_timing    )   CALL timing_start('icedyn_rdgrft')                                                             ! timing
147      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
148
149      IF( kt == nit000 ) THEN
150         IF(lwp) WRITE(numout,*)
151         IF(lwp) WRITE(numout,*)'ice_dyn_rdgrft: ice ridging and rafting'
152         IF(lwp) WRITE(numout,*)'~~~~~~~~~~~~~~'
153      ENDIF     
154
155      CALL ice_var_zapsmall   ! Zero out categories with very small areas
156
157      !--------------------------------
158      ! 0) Identify grid cells with ice
159      !--------------------------------
160      npti = 0   ;   nptidx(:) = 0
161      ipti = 0   ;   iptidx(:) = 0
162      DO jj = 1, jpj
163         DO ji = 1, jpi
164            IF ( at_i(ji,jj) > 0._wp ) THEN
165               npti           = npti + 1
166               nptidx( npti ) = (jj - 1) * jpi + ji
167            ENDIF
168         END DO
169      END DO
170     
171      !--------------------------------------------------------
172      ! 1) Dynamical inputs (closing rate, divergence, opening)
173      !--------------------------------------------------------
174      IF( npti > 0 ) THEN
175       
176         ! just needed here
177         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu(1:npti), divu_i(:,:) )
178         CALL tab_2d_1d( npti, nptidx(1:npti), zdelt(1:npti), delta_i(:,:) )
179         ! needed here and in the iteration loop
180         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i(:,:,:) )
181         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i(:,:,:) )
182         CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i(:,:) )
183
184         DO ji = 1, npti
185            ! closing_net = rate at which open water area is removed + ice area removed by ridging
186            !                                                        - ice area added in new ridges
187            closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp )
188            !
189            ! divergence given by the advection scheme
190            !   (which may not be equal to divu as computed from the velocity field)
191            zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_Dt_ice
192            !
193            IF( zdivu_adv(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) )   ! make sure the closing rate is large enough
194            !                                                                                        ! to give asum = 1.0 after ridging
195            ! Opening rate (non-negative) that will give asum = 1.0 after ridging.
196            opning(ji) = closing_net(ji) + zdivu_adv(ji)
197         END DO
198         !
199         !------------------------------------
200         ! 2) Identify grid cells with ridging
201         !------------------------------------
202         CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
203
204         DO ji = 1, npti
205            IF( SUM( apartf(ji,1:jpl) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
206               ipti = ipti + 1
207               iptidx     (ipti)   = nptidx     (ji)
208               ! adjust to new indices
209               a_i_2d     (ipti,:) = a_i_2d     (ji,:)
210               v_i_2d     (ipti,:) = v_i_2d     (ji,:)
211               ato_i_1d   (ipti)   = ato_i_1d   (ji)
212               closing_net(ipti)   = closing_net(ji)
213               zdivu_adv  (ipti)   = zdivu_adv  (ji)
214               opning     (ipti)   = opning     (ji)
215            ENDIF
216         END DO
217
218      ENDIF
219
220      ! grid cells with ridging
221      nptidx(:) = iptidx(:)
222      npti      = ipti
223
224      !-----------------
225      ! 3) Start ridging
226      !-----------------
227      IF( npti > 0 ) THEN
228         
229         CALL ice_dyn_1d2d( 1 )            ! --- Move to 1D arrays --- !
230
231         iter            = 1
232         iterate_ridging = 1     
233         !                                                        !----------------------!
234         DO WHILE( iterate_ridging > 0 .AND. iter < jp_itermax )  !  ridging iterations  !
235            !                                                     !----------------------!
236            ! Calculate participation function (apartf)
237            !       and transfer      function
238            !       and closing_gross (+correction on opening)
239            CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
240
241            ! Redistribute area, volume, and energy between categories
242            CALL rdgrft_shift
243
244            ! Do we keep on iterating?
245            !-------------------------
246            ! Check whether a_i + ato_i = 0
247            ! If not, because the closing and opening rates were reduced above, ridge again with new rates
248            iterate_ridging = 0
249            DO ji = 1, npti
250               zfac = 1._wp - ( ato_i_1d(ji) + SUM( a_i_2d(ji,:) ) )
251               IF( ABS( zfac ) < epsi10 ) THEN
252                  closing_net(ji) = 0._wp
253                  opning     (ji) = 0._wp
254                  ato_i_1d   (ji) = MAX( 0._wp, 1._wp - SUM( a_i_2d(ji,:) ) )
255               ELSE
256                  iterate_ridging  = 1
257                  zdivu_adv  (ji) = zfac * r1_Dt_ice
258                  closing_net(ji) = MAX( 0._wp, -zdivu_adv(ji) )
259                  opning     (ji) = MAX( 0._wp,  zdivu_adv(ji) )
260               ENDIF
261            END DO
262            !
263            iter = iter + 1
264            IF( iter  >  jp_itermax )    CALL ctl_warn( 'icedyn_rdgrft: non-converging ridging scheme' )
265            !
266         END DO
267
268         CALL ice_dyn_1d2d( 2 )            ! --- Move to 2D arrays --- !
269
270      ENDIF
271   
272      CALL ice_var_agg( 1 ) 
273
274      ! controls
275      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
276      IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rdgrft')                                                             ! prints
277      IF( ln_timing    )   CALL timing_stop ('icedyn_rdgrft')                                                             ! timing
278      !
279   END SUBROUTINE ice_dyn_rdgrft
280
281
282   SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i, pclosing_net )
283      !!-------------------------------------------------------------------
284      !!                ***  ROUTINE rdgrft_prep ***
285      !!
286      !! ** Purpose :   preparation for ridging calculations
287      !!
288      !! ** Method  :   Compute the thickness distribution of the ice and open water
289      !!                participating in ridging and of the resulting ridges.
290      !!-------------------------------------------------------------------
291      REAL(wp), DIMENSION(:)  , INTENT(in) ::   pato_i, pclosing_net 
292      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pa_i, pv_i 
293      !!
294      INTEGER  ::   ji, jl                     ! dummy loop indices
295      REAL(wp) ::   z1_gstar, z1_astar, zhmean, zfac   ! local scalar
296      REAL(wp), DIMENSION(jpij)        ::   zasum, z1_asum, zaksum   ! sum of a_i+ato_i and reverse
297      REAL(wp), DIMENSION(jpij,jpl)    ::   zhi                      ! ice thickness
298      REAL(wp), DIMENSION(jpij,-1:jpl) ::   zGsum                    ! zGsum(n) = sum of areas in categories 0 to n
299      !--------------------------------------------------------------------
300
301      z1_gstar = 1._wp / rn_gstar
302      z1_astar = 1._wp / rn_astar
303
304      !                       ! Ice thickness needed for rafting
305      WHERE( pa_i(1:npti,:) > 0._wp )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
306      ELSEWHERE                         ;   zhi(1:npti,:) = 0._wp
307      END WHERE
308
309      ! 1) Participation function (apartf): a(h) = b(h).g(h)
310      !-----------------------------------------------------------------
311      ! Compute the participation function = total area lost due to ridging/closing
312      ! This is analogous to
313      !   a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
314      !   assuming b(h) = (2/Gstar) * (1 - G(h)/Gstar).
315      !
316      ! apartf = integrating b(h)g(h) between the category boundaries
317      ! apartf is always >= 0 and SUM(apartf(0:jpl))=1
318      !-----------------------------------------------------------------
319      !
320      ! Compute total area of ice plus open water.
321      ! This is in general not equal to one because of divergence during transport
322      zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 )
323      !
324      WHERE( zasum(1:npti) > 0._wp )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti)
325      ELSEWHERE                        ;   z1_asum(1:npti) = 0._wp
326      END WHERE
327      !
328      ! Compute cumulative thickness distribution function
329      ! Compute the cumulative thickness distribution function zGsum,
330      ! where zGsum(n) is the fractional area in categories 0 to n.
331      ! initial value (in h = 0) = open water area
332      zGsum(1:npti,-1) = 0._wp
333      zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti)
334      DO jl = 1, jpl
335         zGsum(1:npti,jl) = ( pato_i(1:npti) + SUM( pa_i(1:npti,1:jl), dim=2 ) ) * z1_asum(1:npti)  ! sum(1:jl) is ok (and not jpl)
336      END DO
337      !
338      IF( ln_partf_lin ) THEN          !--- Linear formulation (Thorndike et al., 1975)
339         DO jl = 0, jpl   
340            DO ji = 1, npti
341               IF    ( zGsum(ji,jl)   < rn_gstar ) THEN
342                  apartf(ji,jl) = z1_gstar * ( zGsum(ji,jl) - zGsum(ji,jl-1) ) * &
343                     &                       ( 2._wp - ( zGsum(ji,jl-1) + zGsum(ji,jl) ) * z1_gstar )
344               ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN
345                  apartf(ji,jl) = z1_gstar * ( rn_gstar     - zGsum(ji,jl-1) ) *  &
346                     &                       ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar        ) * z1_gstar )
347               ELSE
348                  apartf(ji,jl) = 0._wp
349               ENDIF
350            END DO
351         END DO
352         !
353      ELSEIF( ln_partf_exp ) THEN      !--- Exponential, more stable formulation (Lipscomb et al, 2007)
354         !                       
355         zfac = 1._wp / ( 1._wp - EXP(-z1_astar) )
356         DO jl = -1, jpl
357            DO ji = 1, npti
358               zGsum(ji,jl) = EXP( -zGsum(ji,jl) * z1_astar ) * zfac
359            END DO
360         END DO
361         DO jl = 0, jpl
362            DO ji = 1, npti
363               apartf(ji,jl) = zGsum(ji,jl-1) - zGsum(ji,jl)
364            END DO
365         END DO
366         !
367      ENDIF
368
369      !                                !--- Ridging and rafting participation concentrations
370      IF( ln_rafting .AND. ln_ridging ) THEN             !- ridging & rafting
371         DO jl = 1, jpl
372            DO ji = 1, npti
373               aridge(ji,jl) = ( 1._wp + TANH ( rn_craft * ( zhi(ji,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jl)
374               araft (ji,jl) = apartf(ji,jl) - aridge(ji,jl)
375            END DO
376         END DO
377      ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN   !- ridging alone
378         DO jl = 1, jpl
379            DO ji = 1, npti
380               aridge(ji,jl) = apartf(ji,jl)
381               araft (ji,jl) = 0._wp
382            END DO
383         END DO
384      ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN   !- rafting alone   
385         DO jl = 1, jpl
386            DO ji = 1, npti
387               aridge(ji,jl) = 0._wp
388               araft (ji,jl) = apartf(ji,jl)
389            END DO
390         END DO
391      ELSE                                               !- no ridging & no rafting
392         DO jl = 1, jpl
393            DO ji = 1, npti
394               aridge(ji,jl) = 0._wp
395               araft (ji,jl) = 0._wp         
396            END DO
397         END DO
398      ENDIF
399
400      ! 2) Transfer function
401      !-----------------------------------------------------------------
402      ! Compute max and min ridged ice thickness for each ridging category.
403      ! Assume ridged ice is uniformly distributed between hrmin and hrmax.
404      !
405      ! This parameterization is a modified version of Hibler (1980).
406      ! The mean ridging thickness, zhmean, is proportional to hi^(0.5)
407      !  and for very thick ridging ice must be >= hrdg_hi_min*hi
408      !
409      ! The minimum ridging thickness, hrmin, is equal to 2*hi
410      !  (i.e., rafting) and for very thick ridging ice is
411      !  constrained by hrmin <= (zhmean + hi)/2.
412      !
413      ! The maximum ridging thickness, hrmax, is determined by zhmean and hrmin.
414      !
415      ! These modifications have the effect of reducing the ice strength
416      ! (relative to the Hibler formulation) when very thick ice is ridging.
417      !
418      ! zaksum = net area removed/ total area removed
419      ! where total area removed = area of ice that ridges
420      !         net area removed = total area removed - area of new ridges
421      !-----------------------------------------------------------------
422      zfac = 1._wp / hi_hrft
423      zaksum(1:npti) = apartf(1:npti,0)
424      !
425      DO jl = 1, jpl
426         DO ji = 1, npti
427            IF ( apartf(ji,jl) > 0._wp ) THEN
428               zhmean         = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min )
429               hrmin  (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) )
430               hrmax  (ji,jl) = 2._wp * zhmean - hrmin(ji,jl)
431               hraft  (ji,jl) = zhi(ji,jl) * zfac
432               hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 )
433               !
434               ! Normalization factor : zaksum, ensures mass conservation
435               zaksum(ji) = zaksum(ji) + aridge(ji,jl) * ( 1._wp - hi_hrdg(ji,jl) )    &
436                  &                    + araft (ji,jl) * ( 1._wp - hi_hrft )
437            ELSE
438               hrmin  (ji,jl) = 0._wp 
439               hrmax  (ji,jl) = 0._wp 
440               hraft  (ji,jl) = 0._wp 
441               hi_hrdg(ji,jl) = 1._wp
442            ENDIF
443         END DO
444      END DO
445      !
446      ! 3) closing_gross
447      !-----------------
448      ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate. 
449      ! NOTE: 0 < aksum <= 1
450      WHERE( zaksum(1:npti) > 0._wp )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
451      ELSEWHERE                         ;   closing_gross(1:npti) = 0._wp
452      END WHERE
453     
454      ! correction to closing rate if excessive ice removal
455      !----------------------------------------------------
456      ! Reduce the closing rate if more than 100% of any ice category would be removed
457      ! Reduce the opening rate in proportion
458      DO jl = 1, jpl
459         DO ji = 1, npti
460            zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice
461            IF( zfac > pa_i(ji,jl) ) THEN
462               closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_Dt_ice
463            ENDIF
464         END DO
465      END DO     
466
467      ! 4) correction to opening if excessive open water removal
468      !---------------------------------------------------------
469      ! Reduce the closing rate if more than 100% of the open water would be removed
470      ! Reduce the opening rate in proportion
471      DO ji = 1, npti 
472         zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice
473         IF( zfac < 0._wp ) THEN           ! would lead to negative ato_i
474            opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_Dt_ice 
475         ELSEIF( zfac > zasum(ji) ) THEN   ! would lead to ato_i > asum
476            opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_Dt_ice 
477         ENDIF
478      END DO
479      !
480   END SUBROUTINE rdgrft_prep
481
482
483   SUBROUTINE rdgrft_shift
484      !!-------------------------------------------------------------------
485      !!                ***  ROUTINE rdgrft_shift ***
486      !!
487      !! ** Purpose :   shift ridging ice among thickness categories of ice thickness
488      !!
489      !! ** Method  :   Remove area, volume, and energy from each ridging category
490      !!                and add to thicker ice categories.
491      !!-------------------------------------------------------------------
492      !
493      INTEGER  ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
494      REAL(wp) ::   hL, hR, farea              ! left and right limits of integration and new area going to jl2
495      REAL(wp) ::   vsw                        ! vol of water trapped into ridges
496      REAL(wp) ::   afrdg, afrft               ! fraction of category area ridged/rafted
497      REAL(wp)                  ::   airdg1, oirdg1, aprdg1, virdg1, sirdg1
498      REAL(wp)                  ::   airft1, oirft1, aprft1
499      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg  ! area etc of new ridges
500      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft  ! area etc of rafted ice
501      !
502      REAL(wp), DIMENSION(jpij) ::   ersw             ! enth of water trapped into ridges
503      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol    ! new ridge volume going to jl2
504      REAL(wp), DIMENSION(jpij) ::   z1_ai            ! 1 / a
505      !
506      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrft     ! snow energy of rafting ice
507      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirft     ! ice  energy of rafting ice
508      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrdg     ! enth*volume of new ridges     
509      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirdg     ! enth*volume of new ridges
510      !
511      INTEGER , DIMENSION(jpij) ::   itest_rdg, itest_rft   ! test for conservation
512      !!-------------------------------------------------------------------
513
514      ! 1) Change in open water area due to closing and opening
515      !--------------------------------------------------------
516      DO ji = 1, npti
517         ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice )
518      END DO
519     
520      ! 2) compute categories in which ice is removed (jl1)
521      !----------------------------------------------------
522      DO jl1 = 1, jpl
523
524         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,jl1) )
525
526         DO ji = 1, npti
527
528            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging
529
530               z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
531
532               ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2)
533               airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice
534               airft1 = araft (ji,jl1) * closing_gross(ji) * rdt_ice
535
536               airdg2(ji) = airdg1 * hi_hrdg(ji,jl1)
537               airft2(ji) = airft1 * hi_hrft
538
539               ! ridging /rafting fractions
540               afrdg = airdg1 * z1_ai(ji)
541               afrft = airft1 * z1_ai(ji)
542
543               ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges
544               vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg
545               ersw(ji) = -rhoic * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?)
546
547               ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei)
548               virdg1     = v_i_2d (ji,jl1)   * afrdg
549               virdg2(ji) = v_i_2d (ji,jl1)   * afrdg * ( 1. + rn_porordg )
550               vsrdg(ji)  = v_s_2d (ji,jl1)   * afrdg
551               sirdg1     = sv_i_2d(ji,jl1)   * afrdg
552               sirdg2(ji) = sv_i_2d(ji,jl1)   * afrdg + vsw * sss_1d(ji)
553               oirdg1     = oa_i_2d(ji,jl1)   * afrdg
554               oirdg2(ji) = oa_i_2d(ji,jl1)   * afrdg * hi_hrdg(ji,jl1) 
555
556               virft(ji)  = v_i_2d (ji,jl1)   * afrft
557               vsrft(ji)  = v_s_2d (ji,jl1)   * afrft
558               sirft(ji)  = sv_i_2d(ji,jl1)   * afrft 
559               oirft1     = oa_i_2d(ji,jl1)   * afrft 
560               oirft2(ji) = oa_i_2d(ji,jl1)   * afrft * hi_hrft 
561
562               IF ( ln_pnd_H12 ) THEN
563                  aprdg1     = a_ip_2d(ji,jl1) * afrdg
564                  aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1)
565                  vprdg (ji) = v_ip_2d(ji,jl1) * afrdg
566                  aprft1     = a_ip_2d(ji,jl1) * afrft
567                  aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft
568                  vprft (ji) = v_ip_2d(ji,jl1) * afrft
569               ENDIF
570
571               ! Ice-ocean exchanges associated with ice porosity
572               wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw              * rhoic * r1_Dt_ice   ! increase in ice volume due to seawater frozen in voids
573               sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoic * r1_Dt_ice
574               hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji)                 * r1_Dt_ice   ! > 0 [W.m-2]
575
576               ! Put the snow lost by ridging into the ocean
577               !  Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow.
578               wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhosn * vsrdg(ji) * ( 1._wp - rn_fsnwrdg )   &   ! fresh water source for ocean
579                  &                                      + rhosn * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice
580
581               ! Put the melt pond water into the ocean
582               ! clem: I think the following lines must be commented since there
583               !       is no net mass flux between melt ponds and the ocean (see icethd_pnd.F90 for ex.)
584               !IF ( ln_pnd_fwb ) THEN
585               !   wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhofw * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean
586               !      &                              + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_Dt_ice
587               !ENDIF
588
589               ! virtual salt flux to keep salinity constant
590               IF( nn_icesal /= 2 )  THEN
591                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )        ! ridge salinity = s_i
592                  sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoic * r1_Dt_ice  &  ! put back sss_m into the ocean
593                     &                            - s_i_1d(ji) * vsw * rhoic * r1_Dt_ice     ! and get  s_i  from the ocean
594               ENDIF
595
596               ! Remove area, volume of new ridge to each category jl1
597               !------------------------------------------------------
598               a_i_2d (ji,jl1) = a_i_2d (ji,jl1) - airdg1    - airft1
599               v_i_2d (ji,jl1) = v_i_2d (ji,jl1) - virdg1    - virft(ji)
600               v_s_2d (ji,jl1) = v_s_2d (ji,jl1) - vsrdg(ji) - vsrft(ji)
601               sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - sirdg1    - sirft(ji)
602               oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1    - oirft1
603               IF ( ln_pnd_H12 ) THEN
604                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1
605                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji)
606               ENDIF
607            ENDIF
608
609         END DO ! ji
610
611         ! special loop for e_s because of layers jk
612         DO jk = 1, nlay_s
613            DO ji = 1, npti
614               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
615                  ! Compute ridging /rafting fractions
616                  afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
617                  afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
618                  ! Compute ridging /rafting ice and new ridges for es
619                  esrdg(ji,jk) = ze_s_2d (ji,jk,jl1) * afrdg
620                  esrft(ji,jk) = ze_s_2d (ji,jk,jl1) * afrft
621                  ! Put the snow lost by ridging into the ocean
622                  hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji,jk) * ( 1._wp - rn_fsnwrdg )   &                 ! heat sink for ocean (<0, W.m-2)
623                     &                                - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice
624                  !
625                  ! Remove energy of new ridge to each category jl1
626                  !-------------------------------------------------
627                  ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 
628               ENDIF
629            END DO
630         END DO
631                 
632         ! special loop for e_i because of layers jk
633         DO jk = 1, nlay_i
634            DO ji = 1, npti
635               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
636                  ! Compute ridging /rafting fractions
637                  afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
638                  afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
639                  ! Compute ridging ice and new ridges for ei
640                  eirdg(ji,jk) = ze_i_2d (ji,jk,jl1) * afrdg + ersw(ji) * r1_nlay_i
641                  eirft(ji,jk) = ze_i_2d (ji,jk,jl1) * afrft
642                  !
643                  ! Remove energy of new ridge to each category jl1
644                  !-------------------------------------------------
645                  ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 
646               ENDIF
647            END DO
648         END DO
649         
650         ! 3) compute categories in which ice is added (jl2)
651         !--------------------------------------------------
652         itest_rdg(1:npti) = 0
653         itest_rft(1:npti) = 0
654         DO jl2  = 1, jpl 
655            !
656            DO ji = 1, npti
657
658               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
659
660                  ! Compute the fraction of ridged ice area and volume going to thickness category jl2
661                  IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN
662                     hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
663                     hR = MIN( hrmax(ji,jl1), hi_max(jl2)   )
664                     farea    = ( hR      - hL      ) / ( hrmax(ji,jl1)                 - hrmin(ji,jl1)                 )
665                     fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) )
666                     !
667                     itest_rdg(ji) = 1   ! test for conservation
668                  ELSE
669                     farea    = 0._wp 
670                     fvol(ji) = 0._wp                 
671                  ENDIF
672
673                  ! Compute the fraction of rafted ice area and volume going to thickness category jl2
674                  IF( hraft(ji,jl1) <= hi_max(jl2) .AND. hraft(ji,jl1) >  hi_max(jl2-1) ) THEN
675                     zswitch(ji) = 1._wp
676                     !
677                     itest_rft(ji) = 1   ! test for conservation
678                  ELSE
679                     zswitch(ji) = 0._wp
680                  ENDIF
681                  !
682                  ! Patch to ensure perfect conservation if ice thickness goes mad
683                  ! Sometimes thickness is larger than hi_max(jpl) because of advection scheme (for very small areas)
684                  ! Then ice volume is removed from one category but the ridging/rafting scheme
685                  ! does not know where to move it, leading to a conservation issue. 
686                  IF( itest_rdg(ji) == 0 .AND. jl2 == jpl ) THEN   ;   farea = 1._wp   ;   fvol(ji) = 1._wp   ;   ENDIF
687                  IF( itest_rft(ji) == 0 .AND. jl2 == jpl )      zswitch(ji) = 1._wp
688                  !
689                  ! Add area, volume of new ridge to category jl2
690                  !----------------------------------------------
691                  a_i_2d (ji,jl2) = a_i_2d (ji,jl2) + ( airdg2(ji) * farea    + airft2(ji) * zswitch(ji) )
692                  oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ( oirdg2(ji) * farea    + oirft2(ji) * zswitch(ji) )
693                  v_i_2d (ji,jl2) = v_i_2d (ji,jl2) + ( virdg2(ji) * fvol(ji) + virft (ji) * zswitch(ji) )
694                  sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ( sirdg2(ji) * fvol(ji) + sirft (ji) * zswitch(ji) )
695                  v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji)  +  &
696                     &                                  vsrft (ji) * rn_fsnwrft * zswitch(ji) )
697                  IF ( ln_pnd_H12 ) THEN
698                     v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   &
699                        &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   )
700                     a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + (   aprdg2(ji) * rn_fpndrdg * farea         & 
701                        &                                   + aprft2(ji) * rn_fpndrft * zswitch(ji)   )
702                  ENDIF
703                 
704               ENDIF
705
706            END DO
707            ! Add snow energy of new ridge to category jl2
708            !---------------------------------------------
709            DO jk = 1, nlay_s
710               DO ji = 1, npti
711                  IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   &
712                     &   ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ( esrdg(ji,jk) * rn_fsnwrdg * fvol(ji)  +  &
713                     &                                               esrft(ji,jk) * rn_fsnwrft * zswitch(ji) )
714               END DO
715            END DO
716            ! Add ice energy of new ridge to category jl2
717            !--------------------------------------------
718            DO jk = 1, nlay_i
719               DO ji = 1, npti
720                  IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   &
721                     &   ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji)                 
722               END DO
723            END DO
724            !
725         END DO ! jl2
726         !
727      END DO ! jl1
728      !
729      ! In case ridging/rafting lead to very small negative values (sometimes it happens)
730      WHERE( a_i_2d(1:npti,:) < 0._wp )   a_i_2d(1:npti,:) = 0._wp
731      WHERE( v_i_2d(1:npti,:) < 0._wp )   v_i_2d(1:npti,:) = 0._wp
732      !
733   END SUBROUTINE rdgrft_shift
734
735
736   SUBROUTINE ice_strength
737      !!----------------------------------------------------------------------
738      !!                ***  ROUTINE ice_strength ***
739      !!
740      !! ** Purpose :   computes ice strength used in dynamics routines of ice thickness
741      !!
742      !! ** Method  :   Compute the strength of the ice pack, defined as the energy (J m-2)
743      !!              dissipated per unit area removed from the ice pack under compression,
744      !!              and assumed proportional to the change in potential energy caused
745      !!              by ridging. Note that only Hibler's formulation is stable and that
746      !!              ice strength has to be smoothed
747      !!----------------------------------------------------------------------
748      INTEGER             ::   ji, jj, jl  ! dummy loop indices
749      INTEGER             ::   ismooth     ! smoothing the resistance to deformation
750      INTEGER             ::   itframe     ! number of time steps for the P smoothing
751      REAL(wp)            ::   zp, z1_3    ! local scalars
752      REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here
753      REAL(wp), DIMENSION(jpi,jpj) ::   zstrp1, zstrp2   ! strength at previous time steps
754      !!----------------------------------------------------------------------
755      !                              !--------------------------------------------------!
756      IF( ln_str_H79 ) THEN          ! Ice strength => Hibler (1979) method             !
757      !                              !--------------------------------------------------!
758         strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
759         ismooth = 1
760         !                           !--------------------------------------------------!
761      ELSE                           ! Zero strength                                    !
762         !                           !--------------------------------------------------!
763         strength(:,:) = 0._wp
764         ismooth = 0
765      ENDIF
766      !                              !--------------------------------------------------!
767      SELECT CASE( ismooth )         ! Smoothing ice strength                           !
768      !                              !--------------------------------------------------!
769      CASE( 1 )               !--- Spatial smoothing
770         DO jj = 2, jpjm1
771            DO ji = 2, jpim1
772               IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
773                  zworka(ji,jj) = ( 4.0 * strength(ji,jj)              &
774                     &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 
775                     &                  + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
776                     &            ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
777               ELSE
778                  zworka(ji,jj) = 0._wp
779               ENDIF
780            END DO
781         END DO
782         
783         DO jj = 2, jpjm1
784            DO ji = 2, jpim1
785               strength(ji,jj) = zworka(ji,jj)
786            END DO
787         END DO
788         CALL lbc_lnk( strength, 'T', 1. )
789         !
790      CASE( 2 )               !--- Temporal smoothing
791         IF ( kt_ice == nit000 ) THEN
792            zstrp1(:,:) = 0._wp
793            zstrp2(:,:) = 0._wp
794         ENDIF
795         !
796         DO jj = 2, jpjm1
797            DO ji = 2, jpim1
798               IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
799                  itframe = 1 ! number of time steps for the running mean
800                  IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1
801                  IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1
802                  zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe
803                  zstrp2  (ji,jj) = zstrp1  (ji,jj)
804                  zstrp1  (ji,jj) = strength(ji,jj)
805                  strength(ji,jj) = zp
806               ENDIF
807            END DO
808         END DO
809         CALL lbc_lnk( strength, 'T', 1. )
810         !
811      END SELECT
812      !
813   END SUBROUTINE ice_strength
814
815   
816   SUBROUTINE ice_dyn_1d2d( kn )
817      !!-----------------------------------------------------------------------
818      !!                   ***  ROUTINE ice_dyn_1d2d ***
819      !!                 
820      !! ** Purpose :   move arrays from 1d to 2d and the reverse
821      !!-----------------------------------------------------------------------
822      INTEGER, INTENT(in) ::   kn   ! 1= from 2D to 1D   ;   2= from 1D to 2D
823      !
824      INTEGER ::   jl, jk   ! dummy loop indices
825      !!-----------------------------------------------------------------------
826      !
827      SELECT CASE( kn )
828      !                    !---------------------!
829      CASE( 1 )            !==  from 2D to 1D  ==!
830         !                 !---------------------!
831         ! fields used but not modified
832         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) )
833         CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) )
834         ! the following fields are modified in this routine
835         !!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
836         !!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) )
837         !!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i  (:,:,:) )
838         CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
839         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
840         CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
841         CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
842         CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
843         DO jl = 1, jpl
844            DO jk = 1, nlay_s
845               CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
846            END DO
847            DO jk = 1, nlay_i
848               CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
849            END DO
850         END DO
851         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) )
852         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) )
853         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) )
854         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
855         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
856         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
857
858         !                 !---------------------!
859      CASE( 2 )            !==  from 1D to 2D  ==!
860         !                 !---------------------!
861         CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
862         CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
863         CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
864         CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
865         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
866         CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
867         CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
868         CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
869         DO jl = 1, jpl
870            DO jk = 1, nlay_s
871               CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
872            END DO
873            DO jk = 1, nlay_i
874               CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
875            END DO
876         END DO
877         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) )
878         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) )
879         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) )
880         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
881         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
882         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
883         !
884      END SELECT
885      !
886   END SUBROUTINE ice_dyn_1d2d
887   
888
889   SUBROUTINE ice_dyn_rdgrft_init
890      !!-------------------------------------------------------------------
891      !!                  ***  ROUTINE ice_dyn_rdgrft_init ***
892      !!
893      !! ** Purpose :   Physical constants and parameters linked
894      !!                to the mechanical ice redistribution
895      !!
896      !! ** Method  :   Read the namdyn_rdgrft namelist
897      !!                and check the parameters values
898      !!                called at the first timestep (nit000)
899      !!
900      !! ** input   :   Namelist namdyn_rdgrft
901      !!-------------------------------------------------------------------
902      INTEGER :: ios                 ! Local integer output status for namelist read
903      !!
904      NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, &
905         &                    rn_csrdg  ,                    &
906         &                    ln_partf_lin, rn_gstar,        &
907         &                    ln_partf_exp, rn_astar,        & 
908         &                    ln_ridging, rn_hstar, rn_porordg, rn_fsnwrdg, rn_fpndrdg,  & 
909         &                    ln_rafting, rn_hraft, rn_craft  , rn_fsnwrft, rn_fpndrft
910      !!-------------------------------------------------------------------
911      !
912      REWIND( numnam_ice_ref )              ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
913      READ  ( numnam_ice_ref, namdyn_rdgrft, IOSTAT = ios, ERR = 901)
914901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist', lwp )
915      REWIND( numnam_ice_cfg )              ! Namelist namdyn_rdgrft in configuration namelist : Ice mechanical ice redistribution
916      READ  ( numnam_ice_cfg, namdyn_rdgrft, IOSTAT = ios, ERR = 902 )
917902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in configuration namelist', lwp )
918      IF(lwm) WRITE ( numoni, namdyn_rdgrft )
919      !
920      IF (lwp) THEN                          ! control print
921         WRITE(numout,*)
922         WRITE(numout,*) 'ice_dyn_rdgrft_init: ice parameters for ridging/rafting '
923         WRITE(numout,*) '~~~~~~~~~~~~~~~~~~'
924         WRITE(numout,*) '   Namelist namdyn_rdgrft:'
925         WRITE(numout,*) '      ice strength parameterization Hibler (1979)              ln_str_H79   = ', ln_str_H79 
926         WRITE(numout,*) '            1st bulk-rheology parameter                        rn_pstar     = ', rn_pstar
927         WRITE(numout,*) '            2nd bulk-rhelogy parameter                         rn_crhg      = ', rn_crhg
928         WRITE(numout,*) '      Fraction of shear energy contributing to ridging         rn_csrdg     = ', rn_csrdg 
929         WRITE(numout,*) '      linear ridging participation function                    ln_partf_lin = ', ln_partf_lin
930         WRITE(numout,*) '            Fraction of ice coverage contributing to ridging   rn_gstar     = ', rn_gstar
931         WRITE(numout,*) '      Exponential ridging participation function               ln_partf_exp = ', ln_partf_exp
932         WRITE(numout,*) '            Equivalent to G* for an exponential function       rn_astar     = ', rn_astar
933         WRITE(numout,*) '      Ridging of ice sheets or not                             ln_ridging   = ', ln_ridging
934         WRITE(numout,*) '            max ridged ice thickness                           rn_hstar     = ', rn_hstar
935         WRITE(numout,*) '            Initial porosity of ridges                         rn_porordg   = ', rn_porordg
936         WRITE(numout,*) '            Fraction of snow volume conserved during ridging   rn_fsnwrdg   = ', rn_fsnwrdg 
937         WRITE(numout,*) '            Fraction of pond volume conserved during ridging   rn_fpndrdg   = ', rn_fpndrdg 
938         WRITE(numout,*) '      Rafting of ice sheets or not                             ln_rafting   = ', ln_rafting
939         WRITE(numout,*) '            Parmeter thickness (threshold between ridge-raft)  rn_hraft     = ', rn_hraft
940         WRITE(numout,*) '            Rafting hyperbolic tangent coefficient             rn_craft     = ', rn_craft 
941         WRITE(numout,*) '            Fraction of snow volume conserved during rafting   rn_fsnwrft   = ', rn_fsnwrft 
942         WRITE(numout,*) '            Fraction of pond volume conserved during rafting   rn_fpndrft   = ', rn_fpndrft 
943      ENDIF
944      !
945      IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN
946         CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
947      ENDIF
948      !                              ! allocate tke arrays
949      IF( ice_dyn_rdgrft_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' )
950      !
951  END SUBROUTINE ice_dyn_rdgrft_init
952
953#else
954   !!----------------------------------------------------------------------
955   !!   Default option         Empty module           NO SI3 sea-ice model
956   !!----------------------------------------------------------------------
957#endif
958
959   !!======================================================================
960END MODULE icedyn_rdgrft
Note: See TracBrowser for help on using the repository browser.