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 branches/2017/dev_merge_2017/NEMOGCM/NEMO/ICE_SRC – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/ICE_SRC/icedyn_rdgrft.F90 @ 9570

Last change on this file since 9570 was 9570, checked in by nicolasmartin, 6 years ago

Global renaming for core routines (./NEMO)

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