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/trunk/src/ICE – NEMO

source: NEMO/trunk/src/ICE/icedyn_rdgrft.F90 @ 9922

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

sea-ice: remove negative values mainly created by the advection schemes (UMx and Prather) + adapt BDY to ensure that all the fields are coherent with the rest of the code (max concentration etc)

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