source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 4333

Last change on this file since 4333 was 4333, checked in by clem, 7 years ago

remove remaining bugs in LIM3, so that it can run in both regional and global config

  • Property svn:keywords set to Id
File size: 73.1 KB
Line 
1MODULE limitd_me
2   !!======================================================================
3   !!                       ***  MODULE limitd_me ***
4   !! LIM-3 : 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_mec
8   !!            4.0  ! 2011-02  (G. Madec) dynamical allocation
9   !!----------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3'                                      LIM-3 sea-ice model
13   !!----------------------------------------------------------------------
14   USE par_oce          ! ocean parameters
15   USE dom_oce          ! ocean domain
16   USE phycst           ! physical constants (ocean directory)
17   USE sbc_oce          ! surface boundary condition: ocean fields
18   USE thd_ice          ! LIM thermodynamics
19   USE ice              ! LIM variables
20   USE par_ice          ! LIM parameters
21   USE dom_ice          ! LIM domain
22   USE limthd_lac       ! LIM
23   USE limvar           ! LIM
24   USE limcons          ! LIM
25   USE in_out_manager   ! I/O manager
26   USE lbclnk           ! lateral boundary condition - MPP exchanges
27   USE lib_mpp          ! MPP library
28   USE wrk_nemo         ! work arrays
29   USE prtctl           ! Print control
30  ! Check budget (Rousset)
31   USE iom              ! I/O manager
32   USE lib_fortran     ! glob_sum
33   USE limdiahsb
34   USE timing          ! Timing
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC   lim_itd_me               ! called by ice_stp
40   PUBLIC   lim_itd_me_icestrength
41   PUBLIC   lim_itd_me_init
42   PUBLIC   lim_itd_me_zapsmall
43   PUBLIC   lim_itd_me_alloc        ! called by iceini.F90
44
45   REAL(wp) ::   epsi20 = 1.e-20_wp   ! constant values
46   REAL(wp) ::   epsi10 = 1.e-10_wp   ! constant values
47   REAL(wp) ::   epsi06 = 1.e-06_wp   ! constant values
48
49   !-----------------------------------------------------------------------
50   ! Variables shared among ridging subroutines
51   !-----------------------------------------------------------------------
52   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area
53   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged
54   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/
55   !                                                     ! closing associated w/ category n
56   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness
57   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness
58   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice
59   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   krdg     ! mean ridge thickness/thickness of ridging ice
60   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging
61   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting
62
63   REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier
64   REAL(wp), PARAMETER ::   kraft   = 2.0_wp    ! rafting multipliyer
65   REAL(wp), PARAMETER ::   kamax   = 1.0_wp    ! max of ice area authorized (clem: scheme is not stable if kamax <= 0.99)
66
67   REAL(wp) ::   Cp                             !
68   !
69   !-----------------------------------------------------------------------
70   ! Ridging diagnostic arrays for history files
71   !-----------------------------------------------------------------------
72   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dardg1dt   ! rate of fractional area loss by ridging ice (1/s)
73   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dardg2dt   ! rate of fractional area gain by new ridges (1/s)
74   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dvirdgdt   ! rate of ice volume ridged (m/s)
75   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   opening    ! rate of opening due to divergence/shear (1/s)
76   !
77   !!----------------------------------------------------------------------
78   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
79   !! $Id$
80   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
81   !!----------------------------------------------------------------------
82CONTAINS
83
84   INTEGER FUNCTION lim_itd_me_alloc()
85      !!---------------------------------------------------------------------!
86      !!                ***  ROUTINE lim_itd_me_alloc ***
87      !!---------------------------------------------------------------------!
88      ALLOCATE(                                                                     &
89         !* Variables shared among ridging subroutines
90         &      asum (jpi,jpj)     , athorn(jpi,jpj,0:jpl)                    ,     &
91         &      aksum(jpi,jpj)                                                ,     &
92         !
93         &      hrmin(jpi,jpj,jpl) , hraft(jpi,jpj,jpl) , aridge(jpi,jpj,jpl) ,     &
94         &      hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) ,     &
95         !
96         !* Ridging diagnostic arrays for history files
97         &      dardg1dt(jpi,jpj)  , dardg2dt(jpi,jpj)                        ,     & 
98         &      dvirdgdt(jpi,jpj)  , opening(jpi,jpj)                         , STAT=lim_itd_me_alloc )
99         !
100      IF( lim_itd_me_alloc /= 0 )   CALL ctl_warn( 'lim_itd_me_alloc: failed to allocate arrays' )
101      !
102   END FUNCTION lim_itd_me_alloc
103
104
105   SUBROUTINE lim_itd_me
106      !!---------------------------------------------------------------------!
107      !!                ***  ROUTINE lim_itd_me ***
108      !!
109      !! ** Purpose :   computes the mechanical redistribution of ice thickness
110      !!
111      !! ** Method  :   Steps :
112      !!       1) Thickness categories boundaries, ice / o.w. concentrations
113      !!          Ridge preparation
114      !!       2) Dynamical inputs (closing rate, divu_adv, opning)
115      !!       3) Ridging iteration
116      !!       4) Ridging diagnostics
117      !!       5) Heat, salt and freshwater fluxes
118      !!       6) Compute increments of tate variables and come back to old values
119      !!
120      !! References :   Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626.
121      !!                Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980.
122      !!                Rothrock, D. A., 1975: JGR, 80, 4514-4519.
123      !!                Thorndike et al., 1975, JGR, 80, 4501-4513.
124      !!                Bitz et al., JGR, 2001
125      !!                Amundrud and Melling, JGR 2005
126      !!                Babko et al., JGR 2002
127      !!
128      !!     This routine is based on CICE code and authors William H. Lipscomb,
129      !!  and Elizabeth C. Hunke, LANL are gratefully acknowledged
130      !!--------------------------------------------------------------------!
131      INTEGER ::   ji, jj, jk, jl   ! dummy loop index
132      INTEGER ::   niter, nitermax = 20   ! local integer
133      LOGICAL  ::   asum_error            ! flag for asum .ne. 1
134      INTEGER  ::   iterate_ridging       ! if true, repeat the ridging
135      REAL(wp) ::   w1, tmpfac            ! local scalar
136      CHARACTER (len = 15) ::   fieldid
137      REAL(wp), POINTER, DIMENSION(:,:) ::   closing_net     ! net rate at which area is removed    (1/s)
138                                                             ! (ridging ice area - area of new ridges) / dt
139      REAL(wp), POINTER, DIMENSION(:,:) ::   divu_adv        ! divu as implied by transport scheme  (1/s)
140      REAL(wp), POINTER, DIMENSION(:,:) ::   opning          ! rate of opening due to divergence/shear
141      REAL(wp), POINTER, DIMENSION(:,:) ::   closing_gross   ! rate at which area removed, not counting area of new ridges
142      REAL(wp), POINTER, DIMENSION(:,:) ::   msnow_mlt       ! mass of snow added to ocean (kg m-2)
143      REAL(wp), POINTER, DIMENSION(:,:) ::   esnow_mlt       ! energy needed to melt snow in ocean (J m-2)
144      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final  !  ice volume summed over categories
145      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)
146      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
147      ! mass and salt flux (clem)
148      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold, zsmvold   ! old ice volume...
149      !!-----------------------------------------------------------------------------
150      IF( nn_timing == 1 )  CALL timing_start('limitd_me')
151
152      CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
153
154      CALL wrk_alloc( jpi, jpj, jpl, zviold, zvsold, zsmvold )   ! clem
155
156      IF( numit == nstart  )   CALL lim_itd_me_init   ! Initialization (first time-step only)
157
158      IF(ln_ctl) THEN
159         CALL prt_ctl(tab2d_1=ato_i , clinfo1=' lim_itd_me: ato_i  : ', tab2d_2=at_i   , clinfo2=' at_i    : ')
160         CALL prt_ctl(tab2d_1=divu_i, clinfo1=' lim_itd_me: divu_i : ', tab2d_2=delta_i, clinfo2=' delta_i : ')
161      ENDIF
162
163      IF( ln_limdyn ) THEN          !   Start ridging and rafting   !
164      ! -------------------------------
165      !- check conservation (C Rousset)
166      IF (ln_limdiahsb) THEN
167         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
168         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
169         zchk_fw_b  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) )
170         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) )
171      ENDIF
172      !- check conservation (C Rousset)
173      ! -------------------------------
174
175      ! mass and salt flux init (clem)
176      zviold(:,:,:) = v_i(:,:,:)
177      zvsold(:,:,:) = v_s(:,:,:)
178      zsmvold(:,:,:) = smv_i(:,:,:)
179
180      !-----------------------------------------------------------------------------!
181      ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
182      !-----------------------------------------------------------------------------!
183      Cp = 0.5 * grav * (rau0-rhoic) * rhoic / rau0                ! proport const for PE
184      !
185      CALL lim_itd_me_ridgeprep                                    ! prepare ridging
186      !
187      IF( con_i)   CALL lim_column_sum( jpl, v_i, vt_i_init )      ! conservation check
188
189      DO jj = 1, jpj                                     ! Initialize arrays.
190         DO ji = 1, jpi
191            msnow_mlt(ji,jj) = 0._wp
192            esnow_mlt(ji,jj) = 0._wp
193            dardg1dt (ji,jj) = 0._wp
194            dardg2dt (ji,jj) = 0._wp
195            dvirdgdt (ji,jj) = 0._wp
196            opening  (ji,jj) = 0._wp
197
198            !-----------------------------------------------------------------------------!
199            ! 2) Dynamical inputs (closing rate, divu_adv, opning)
200            !-----------------------------------------------------------------------------!
201            !
202            ! 2.1 closing_net
203            !-----------------
204            ! Compute the net rate of closing due to convergence
205            ! and shear, based on Flato and Hibler (1995).
206            !
207            ! The energy dissipation rate is equal to the net closing rate
208            ! times the ice strength.
209            !
210            ! NOTE: The NET closing rate is equal to the rate that open water
211            !  area is removed, plus the rate at which ice area is removed by
212            !  ridging, minus the rate at which area is added in new ridges.
213            !  The GROSS closing rate is equal to the first two terms (open
214            !  water closing and thin ice ridging) without the third term
215            !  (thick, newly ridged ice).
216
217            closing_net(ji,jj) = Cs * 0.5 * ( Delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
218
219            ! 2.2 divu_adv
220            !--------------
221            ! Compute divu_adv, the divergence rate given by the transport/
222            ! advection scheme, which may not be equal to divu as computed
223            ! from the velocity field.
224            !
225            ! If divu_adv < 0, make sure the closing rate is large enough
226            ! to give asum = 1.0 after ridging.
227
228            divu_adv(ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice  ! asum found in ridgeprep
229
230            IF( divu_adv(ji,jj) < 0._wp )   closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) )
231
232            ! 2.3 opning
233            !------------
234            ! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging.
235            opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj)
236         END DO
237      END DO
238
239      !-----------------------------------------------------------------------------!
240      ! 3) Ridging iteration
241      !-----------------------------------------------------------------------------!
242      niter           = 1                 ! iteration counter
243      iterate_ridging = 1
244
245      DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax )
246
247         DO jj = 1, jpj
248            DO ji = 1, jpi
249
250               ! 3.2 closing_gross
251               !-----------------------------------------------------------------------------!
252               ! Based on the ITD of ridging and ridged ice, convert the net
253               !  closing rate to a gross closing rate. 
254               ! NOTE: 0 < aksum <= 1
255               closing_gross(ji,jj) = closing_net(ji,jj) / aksum(ji,jj)
256
257               ! correction to closing rate and opening if closing rate is excessive
258               !---------------------------------------------------------------------
259               ! Reduce the closing rate if more than 100% of the open water
260               ! would be removed.  Reduce the opening rate proportionately.
261               IF ( ato_i(ji,jj) .GT. epsi10 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN
262                  w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice
263                  IF ( w1 .GT. ato_i(ji,jj)) THEN
264                     tmpfac = ato_i(ji,jj) / w1
265                     closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac
266                     opning(ji,jj) = opning(ji,jj) * tmpfac
267                  ENDIF !w1
268               ENDIF !at0i and athorn
269
270            END DO ! ji
271         END DO ! jj
272
273         ! correction to closing rate / opening if excessive ice removal
274         !---------------------------------------------------------------
275         ! Reduce the closing rate if more than 100% of any ice category
276         ! would be removed.  Reduce the opening rate proportionately.
277
278         DO jl = 1, jpl
279            DO jj = 1, jpj
280               DO ji = 1, jpi
281                  IF ( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp )THEN
282                     w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
283                     IF ( w1  >  a_i(ji,jj,jl) ) THEN
284                        tmpfac = a_i(ji,jj,jl) / w1
285                        closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac
286                        opning       (ji,jj) = opning       (ji,jj) * tmpfac
287                     ENDIF
288                  ENDIF
289               END DO !ji
290            END DO ! jj
291         END DO !jl
292
293         ! 3.3 Redistribute area, volume, and energy.
294         !-----------------------------------------------------------------------------!
295
296         CALL lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
297
298         ! 3.4 Compute total area of ice plus open water after ridging.
299         !-----------------------------------------------------------------------------!
300
301         CALL lim_itd_me_asumr
302
303         ! 3.5 Do we keep on iterating ???
304         !-----------------------------------------------------------------------------!
305         ! Check whether asum = 1.  If not (because the closing and opening
306         ! rates were reduced above), ridge again with new rates.
307
308         iterate_ridging = 0
309
310         DO jj = 1, jpj
311            DO ji = 1, jpi
312               IF (ABS(asum(ji,jj) - kamax ) .LT. epsi10) THEN
313                  closing_net(ji,jj) = 0._wp
314                  opning     (ji,jj) = 0._wp
315               ELSE
316                  iterate_ridging    = 1
317                  divu_adv   (ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice
318                  closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) )
319                  opning     (ji,jj) = MAX( 0._wp,  divu_adv(ji,jj) )
320               ENDIF
321            END DO
322         END DO
323
324         IF( lk_mpp )   CALL mpp_max( iterate_ridging )
325
326         ! Repeat if necessary.
327         ! NOTE: If strength smoothing is turned on, the ridging must be
328         ! iterated globally because of the boundary update in the
329         ! smoothing.
330
331         niter = niter + 1
332
333         IF( iterate_ridging == 1 ) THEN
334            IF( niter  >  nitermax ) THEN
335               WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
336               WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
337            ENDIF
338            CALL lim_itd_me_ridgeprep
339         ENDIF
340
341      END DO !! on the do while over iter
342
343      !-----------------------------------------------------------------------------!
344      ! 4) Ridging diagnostics
345      !-----------------------------------------------------------------------------!
346      ! Convert ridging rate diagnostics to correct units.
347      ! Update fresh water and heat fluxes due to snow melt.
348
349      asum_error = .false. 
350
351      DO jj = 1, jpj
352         DO ji = 1, jpi
353
354            IF(ABS(asum(ji,jj) - kamax) > epsi10 ) asum_error = .true.
355
356            dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice
357            dardg2dt(ji,jj) = dardg2dt(ji,jj) * r1_rdtice
358            dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * r1_rdtice
359            opening (ji,jj) = opening (ji,jj) * r1_rdtice
360
361            !-----------------------------------------------------------------------------!
362            ! 5) Heat, salt and freshwater fluxes
363            !-----------------------------------------------------------------------------!
364            fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean
365            fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice     ! heat sink for ocean
366
367         END DO
368      END DO
369
370      ! Check if there is a ridging error
371      DO jj = 1, jpj
372         DO ji = 1, jpi
373            IF( ABS( asum(ji,jj) - kamax)  >  epsi10 ) THEN   ! there is a bug
374               WRITE(numout,*) ' '
375               WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj)
376               WRITE(numout,*) ' limitd_me '
377               WRITE(numout,*) ' POINT : ', ji, jj
378               WRITE(numout,*) ' jpl, a_i, athorn '
379               WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0)
380               DO jl = 1, jpl
381                  WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl)
382               END DO
383            ENDIF  ! asum
384
385         END DO !ji
386      END DO !jj
387
388      ! Conservation check
389      IF ( con_i ) THEN
390         CALL lim_column_sum (jpl,   v_i, vt_i_final)
391         fieldid = ' v_i : limitd_me '
392         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
393      ENDIF
394
395      !-----------------------------------------------------------------------------!
396      ! 6) Updating state variables and trend terms (done in limupdate)
397      !-----------------------------------------------------------------------------!
398      CALL lim_var_glo2eqv
399      CALL lim_itd_me_zapsmall
400
401      !--------------------------------
402      ! Update mass/salt fluxes (clem)
403      !--------------------------------
404      DO jl = 1, jpl
405         DO jj = 1, jpj 
406            DO ji = 1, jpi
407               diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice
408               rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic 
409               rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn 
410               sfx_mec(ji,jj) = sfx_mec(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic * r1_rdtice 
411            END DO
412         END DO
413      END DO
414
415      IF(ln_ctl) THEN     ! Control print
416         CALL prt_ctl_info(' ')
417         CALL prt_ctl_info(' - Cell values : ')
418         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
419         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_me  : cell area :')
420         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me  : at_i      :')
421         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me  : vt_i      :')
422         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_me  : vt_s      :')
423         DO jl = 1, jpl
424            CALL prt_ctl_info(' ')
425            CALL prt_ctl_info(' - Category : ', ivar1=jl)
426            CALL prt_ctl_info('   ~~~~~~~~~~')
427            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_itd_me  : a_i      : ')
428            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_itd_me  : ht_i     : ')
429            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_itd_me  : ht_s     : ')
430            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_itd_me  : v_i      : ')
431            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_itd_me  : v_s      : ')
432            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_itd_me  : e_s      : ')
433            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_itd_me  : t_su     : ')
434            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_itd_me  : t_snow   : ')
435            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_itd_me  : sm_i     : ')
436            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_itd_me  : smv_i    : ')
437            DO jk = 1, nlay_i
438               CALL prt_ctl_info(' ')
439               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
440               CALL prt_ctl_info('   ~~~~~~~')
441               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_me  : t_i      : ')
442               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_me  : e_i      : ')
443            END DO
444         END DO
445      ENDIF
446
447      ! -------------------------------
448      !- check conservation (C Rousset)
449      IF (ln_limdiahsb) THEN
450         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
451         zchk_fw  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b
452 
453         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice
454         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic )
455
456         zchk_vmin = glob_min(v_i)
457         zchk_amax = glob_max(SUM(a_i,dim=3))
458         zchk_amin = glob_min(a_i)
459       
460         IF(lwp) THEN
461            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limitd_me) = ',(zchk_v_i * rday)
462            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * rday)
463            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_me) = ',(zchk_vmin * 1.e-3)
464            IF ( zchk_amax >  kamax+epsi10  ) WRITE(numout,*) 'violation a_i>amax            (limitd_me) = ',zchk_amax
465            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limitd_me) = ',zchk_amin
466         ENDIF
467      ENDIF
468      !- check conservation (C Rousset)
469      ! -------------------------------
470
471      ENDIF  ! ln_limdyn=.true.
472      !
473      CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
474      !
475      CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zsmvold )   ! clem
476      !
477      IF( nn_timing == 1 )  CALL timing_stop('limitd_me')
478   END SUBROUTINE lim_itd_me
479
480
481   SUBROUTINE lim_itd_me_icestrength( kstrngth )
482      !!----------------------------------------------------------------------
483      !!                ***  ROUTINE lim_itd_me_icestrength ***
484      !!
485      !! ** Purpose :   computes ice strength used in dynamics routines of ice thickness
486      !!
487      !! ** Method  :   Compute the strength of the ice pack, defined as the energy (J m-2)
488      !!              dissipated per unit area removed from the ice pack under compression,
489      !!              and assumed proportional to the change in potential energy caused
490      !!              by ridging. Note that only Hibler's formulation is stable and that
491      !!              ice strength has to be smoothed
492      !!
493      !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
494      !!----------------------------------------------------------------------
495      INTEGER, INTENT(in) ::   kstrngth    ! = 1 for Rothrock formulation, 0 for Hibler (1979)
496
497      INTEGER ::   ji,jj, jl   ! dummy loop indices
498      INTEGER ::   ksmooth     ! smoothing the resistance to deformation
499      INTEGER ::   numts_rm    ! number of time steps for the P smoothing
500      REAL(wp) ::   hi, zw1, zp, zdummy, zzc, z1_3   ! local scalars
501      REAL(wp), POINTER, DIMENSION(:,:) ::   zworka   ! temporary array used here
502      !!----------------------------------------------------------------------
503
504      CALL wrk_alloc( jpi, jpj, zworka )
505
506      !------------------------------------------------------------------------------!
507      ! 1) Initialize
508      !------------------------------------------------------------------------------!
509      strength(:,:) = 0._wp
510
511      !------------------------------------------------------------------------------!
512      ! 2) Compute thickness distribution of ridging and ridged ice
513      !------------------------------------------------------------------------------!
514      CALL lim_itd_me_ridgeprep
515
516      !------------------------------------------------------------------------------!
517      ! 3) Rothrock(1975)'s method
518      !------------------------------------------------------------------------------!
519      IF( kstrngth == 1 ) THEN
520         z1_3 = 1._wp / 3._wp
521         DO jl = 1, jpl
522            DO jj= 1, jpj
523               DO ji = 1, jpi
524                  !
525                  IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN
526                     hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
527                     !----------------------------
528                     ! PE loss from deforming ice
529                     !----------------------------
530                     strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * hi * hi
531
532                     !--------------------------
533                     ! PE gain from rafting ice
534                     !--------------------------
535                     strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * hi * hi
536
537                     !----------------------------
538                     ! PE gain from ridging ice
539                     !----------------------------
540                     strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl)/krdg(ji,jj,jl)     &
541                        * z1_3 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) / ( hrmax(ji,jj,jl)-hrmin(ji,jj,jl) )   
542!!gm Optimization:  (a**3-b**3)/(a-b) = a*a+ab+b*b   ==> less costly operations even if a**3 is replaced by a*a*a...                   
543                  ENDIF            ! aicen > epsi10
544                  !
545               END DO ! ji
546            END DO !jj
547         END DO !jl
548
549         zzc = Cf * Cp     ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and Cf accounts for frictional dissipation
550         strength(:,:) = zzc * strength(:,:) / aksum(:,:)
551
552         ksmooth = 1
553
554         !------------------------------------------------------------------------------!
555         ! 4) Hibler (1979)' method
556         !------------------------------------------------------------------------------!
557      ELSE                      ! kstrngth ne 1:  Hibler (1979) form
558         !
559         strength(:,:) = Pstar * vt_i(:,:) * EXP( - C_rhg * ( 1._wp - at_i(:,:) )  )
560         !
561         ksmooth = 1
562         !
563      ENDIF                     ! kstrngth
564
565      !
566      !------------------------------------------------------------------------------!
567      ! 5) Impact of brine volume
568      !------------------------------------------------------------------------------!
569      ! CAN BE REMOVED
570      !
571      IF ( brinstren_swi == 1 ) THEN
572
573         DO jj = 1, jpj
574            DO ji = 1, jpi
575               IF ( bv_i(ji,jj) .GT. 0.0 ) THEN
576                  zdummy = MIN ( bv_i(ji,jj), 0.10 ) * MIN( bv_i(ji,jj), 0.10 )
577               ELSE
578                  zdummy = 0.0
579               ENDIF
580               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0)))
581            END DO              ! j
582         END DO                 ! i
583
584      ENDIF
585
586      !
587      !------------------------------------------------------------------------------!
588      ! 6) Smoothing ice strength
589      !------------------------------------------------------------------------------!
590      !
591      !-------------------
592      ! Spatial smoothing
593      !-------------------
594      IF ( ksmooth == 1 ) THEN
595
596         CALL lbc_lnk( strength, 'T', 1. )
597
598         DO jj = 2, jpj - 1
599            DO ji = 2, jpi - 1
600               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN ! ice is
601                  ! present
602                  zworka(ji,jj) = 4.0 * strength(ji,jj)              &
603                     &          + strength(ji-1,jj) * tms(ji-1,jj) & 
604                     &          + strength(ji+1,jj) * tms(ji+1,jj) & 
605                     &          + strength(ji,jj-1) * tms(ji,jj-1) & 
606                     &          + strength(ji,jj+1) * tms(ji,jj+1)   
607
608                  zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj) + tms(ji,jj-1) + tms(ji,jj+1)
609                  zworka(ji,jj) = zworka(ji,jj) / zw1
610               ELSE
611                  zworka(ji,jj) = 0._wp
612               ENDIF
613            END DO
614         END DO
615
616         DO jj = 2, jpj - 1
617            DO ji = 2, jpi - 1
618               strength(ji,jj) = zworka(ji,jj)
619            END DO
620         END DO
621         CALL lbc_lnk( strength, 'T', 1. )
622
623      ENDIF ! ksmooth
624
625      !--------------------
626      ! Temporal smoothing
627      !--------------------
628      IF ( numit == nit000 + nn_fsbc - 1 ) THEN
629         strp1(:,:) = 0.0           
630         strp2(:,:) = 0.0           
631      ENDIF
632
633      IF ( ksmooth == 2 ) THEN
634
635
636         CALL lbc_lnk( strength, 'T', 1. )
637
638         DO jj = 1, jpj - 1
639            DO ji = 1, jpi - 1
640               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN       ! ice is present
641                  numts_rm = 1 ! number of time steps for the running mean
642                  IF ( strp1(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1
643                  IF ( strp2(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1
644                  zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm
645                  strp2(ji,jj) = strp1(ji,jj)
646                  strp1(ji,jj) = strength(ji,jj)
647                  strength(ji,jj) = zp
648
649               ENDIF
650            END DO
651         END DO
652
653      ENDIF ! ksmooth
654
655      CALL lbc_lnk( strength, 'T', 1. )      ! Boundary conditions
656
657      CALL wrk_dealloc( jpi, jpj, zworka )
658      !
659   END SUBROUTINE lim_itd_me_icestrength
660
661
662   SUBROUTINE lim_itd_me_ridgeprep
663      !!---------------------------------------------------------------------!
664      !!                ***  ROUTINE lim_itd_me_ridgeprep ***
665      !!
666      !! ** Purpose :   preparation for ridging and strength calculations
667      !!
668      !! ** Method  :   Compute the thickness distribution of the ice and open water
669      !!              participating in ridging and of the resulting ridges.
670      !!---------------------------------------------------------------------!
671      INTEGER ::   ji,jj, jl    ! dummy loop indices
672      INTEGER ::   krdg_index   !
673      REAL(wp) ::   Gstari, astari, hi, hrmean, zdummy   ! local scalar
674      REAL(wp), POINTER, DIMENSION(:,:)   ::   zworka    ! temporary array used here
675      REAL(wp), POINTER, DIMENSION(:,:,:) ::   Gsum      ! Gsum(n) = sum of areas in categories 0 to n
676      !------------------------------------------------------------------------------!
677
678      CALL wrk_alloc( jpi,jpj, zworka )
679      CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
680
681      Gstari     = 1.0/Gstar   
682      astari     = 1.0/astar   
683      aksum(:,:)    = 0.0
684      athorn(:,:,:) = 0.0
685      aridge(:,:,:) = 0.0
686      araft (:,:,:) = 0.0
687      hrmin(:,:,:)  = 0.0 
688      hrmax(:,:,:)  = 0.0 
689      hraft(:,:,:)  = 0.0 
690      krdg (:,:,:)  = 1.0
691
692      !     ! Zero out categories with very small areas
693      CALL lim_itd_me_zapsmall
694
695      !------------------------------------------------------------------------------!
696      ! 1) Participation function
697      !------------------------------------------------------------------------------!
698
699      ! Compute total area of ice plus open water.
700      CALL lim_itd_me_asumr
701      ! This is in general not equal to one
702      ! because of divergence during transport
703
704      ! Compute cumulative thickness distribution function
705      ! Compute the cumulative thickness distribution function Gsum,
706      ! where Gsum(n) is the fractional area in categories 0 to n.
707      ! initial value (in h = 0) equals open water area
708
709      Gsum(:,:,-1) = 0._wp
710
711      DO jj = 1, jpj
712         DO ji = 1, jpi
713            IF( ato_i(ji,jj) > epsi10 ) THEN   ;   Gsum(ji,jj,0) = ato_i(ji,jj)
714            ELSE                               ;   Gsum(ji,jj,0) = 0._wp
715            ENDIF
716         END DO
717      END DO
718
719      ! for each value of h, you have to add ice concentration then
720      DO jl = 1, jpl
721         DO jj = 1, jpj 
722            DO ji = 1, jpi
723               IF( a_i(ji,jj,jl) .GT. epsi10 ) THEN   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl)
724               ELSE                                   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1)
725               ENDIF
726            END DO
727         END DO
728      END DO
729
730      ! Normalize the cumulative distribution to 1
731      zworka(:,:) = 1._wp / Gsum(:,:,jpl)
732      DO jl = 0, jpl
733         Gsum(:,:,jl) = Gsum(:,:,jl) * zworka(:,:)
734      END DO
735
736      ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
737      !--------------------------------------------------------------------------------------------------
738      ! Compute the participation function athorn; this is analogous to
739      ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
740      ! area lost from category n due to ridging/closing
741      ! athorn(n)   = total area lost due to ridging/closing
742      ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
743      !
744      ! The expressions for athorn are found by integrating b(h)g(h) between
745      ! the category boundaries.
746      !-----------------------------------------------------------------
747
748      krdg_index = 1
749
750      IF( krdg_index == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975)
751         DO jl = 0, ice_cat_bounds(1,2)       ! only undeformed ice participates
752            DO jj = 1, jpj 
753               DO ji = 1, jpi
754                  IF( Gsum(ji,jj,jl) < Gstar) THEN
755                     athorn(ji,jj,jl) = Gstari * (Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * &
756                        (2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari)
757                  ELSEIF (Gsum(ji,jj,jl-1) < Gstar) THEN
758                     athorn(ji,jj,jl) = Gstari * (Gstar-Gsum(ji,jj,jl-1)) *  &
759                        (2.0 - (Gsum(ji,jj,jl-1)+Gstar)*Gstari)
760                  ELSE
761                     athorn(ji,jj,jl) = 0.0
762                  ENDIF
763               END DO ! ji
764            END DO ! jj
765         END DO ! jl
766
767      ELSE                             !--- Exponential, more stable formulation (Lipscomb et al, 2007)
768         !                       
769         zdummy = 1._wp / ( 1._wp - EXP(-astari) )        ! precompute exponential terms using Gsum as a work array
770
771         DO jl = -1, jpl
772            Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
773         END DO !jl
774         DO jl = 0, ice_cat_bounds(1,2)
775             athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
776         END DO
777         !
778      ENDIF ! krdg_index
779
780      IF( raftswi == 1 ) THEN      ! Ridging and rafting ice participation functions
781         !
782         DO jl = 1, jpl
783            DO jj = 1, jpj 
784               DO ji = 1, jpi
785                  IF ( athorn(ji,jj,jl) .GT. 0._wp ) THEN
786!!gm  TANH( -X ) = - TANH( X )  so can be computed only 1 time....
787                     aridge(ji,jj,jl) = ( TANH (  Craft * ( ht_i(ji,jj,jl) - hparmeter ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
788                     araft (ji,jj,jl) = ( TANH ( -Craft * ( ht_i(ji,jj,jl) - hparmeter ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
789                     IF ( araft(ji,jj,jl) < epsi06 )   araft(ji,jj,jl)  = 0._wp
790                     aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 )
791                  ENDIF ! athorn
792               END DO ! ji
793            END DO ! jj
794         END DO ! jl
795
796      ELSE  ! raftswi = 0
797         !
798         DO jl = 1, jpl
799            aridge(:,:,jl) = athorn(:,:,jl)
800         END DO
801         !
802      ENDIF
803
804      IF ( raftswi == 1 ) THEN
805
806         IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi10 ) THEN
807            DO jl = 1, jpl
808               DO jj = 1, jpj
809                  DO ji = 1, jpi
810                     IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. epsi10 ) THEN
811                        WRITE(numout,*) ' ALERTE 96 : wrong participation function ... '
812                        WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl
813                        WRITE(numout,*) ' lat, lon   : ', gphit(ji,jj), glamt(ji,jj)
814                        WRITE(numout,*) ' aridge     : ', aridge(ji,jj,1:jpl)
815                        WRITE(numout,*) ' araft      : ', araft(ji,jj,1:jpl)
816                        WRITE(numout,*) ' athorn     : ', athorn(ji,jj,1:jpl)
817                     ENDIF
818                  END DO
819               END DO
820            END DO
821         ENDIF
822
823      ENDIF
824
825      !-----------------------------------------------------------------
826      ! 2) Transfer function
827      !-----------------------------------------------------------------
828      ! Compute max and min ridged ice thickness for each ridging category.
829      ! Assume ridged ice is uniformly distributed between hrmin and hrmax.
830      !
831      ! This parameterization is a modified version of Hibler (1980).
832      ! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
833      !  and for very thick ridging ice must be >= krdgmin*hi
834      !
835      ! The minimum ridging thickness, hrmin, is equal to 2*hi
836      !  (i.e., rafting) and for very thick ridging ice is
837      !  constrained by hrmin <= (hrmean + hi)/2.
838      !
839      ! The maximum ridging thickness, hrmax, is determined by
840      !  hrmean and hrmin.
841      !
842      ! These modifications have the effect of reducing the ice strength
843      ! (relative to the Hibler formulation) when very thick ice is
844      ! ridging.
845      !
846      ! aksum = net area removed/ total area removed
847      ! where total area removed = area of ice that ridges
848      !         net area removed = total area removed - area of new ridges
849      !-----------------------------------------------------------------
850
851      ! Transfer function
852      DO jl = 1, jpl !all categories have a specific transfer function
853         DO jj = 1, jpj
854            DO ji = 1, jpi
855
856               IF (a_i(ji,jj,jl) .GT. epsi10 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN
857                  hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
858                  hrmean          = MAX(SQRT(Hstar*hi), hi*krdgmin)
859                  hrmin(ji,jj,jl) = MIN(2.0*hi, 0.5*(hrmean + hi))
860                  hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl)
861                  hraft(ji,jj,jl) = kraft*hi
862                  krdg(ji,jj,jl)  = hrmean / hi
863               ELSE
864                  hraft(ji,jj,jl) = 0.0
865                  hrmin(ji,jj,jl) = 0.0 
866                  hrmax(ji,jj,jl) = 0.0 
867                  krdg (ji,jj,jl) = 1.0
868               ENDIF
869
870            END DO ! ji
871         END DO ! jj
872      END DO ! jl
873
874      ! Normalization factor : aksum, ensures mass conservation
875      aksum(:,:) = athorn(:,:,0)
876      DO jl = 1, jpl 
877         aksum(:,:)    = aksum(:,:) + aridge(:,:,jl) * ( 1._wp - 1._wp / krdg(:,:,jl) )    &
878            &                       + araft (:,:,jl) * ( 1._wp - 1._wp / kraft        )
879      END DO
880      !
881      CALL wrk_dealloc( jpi,jpj, zworka )
882      CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
883      !
884   END SUBROUTINE lim_itd_me_ridgeprep
885
886
887   SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
888      !!----------------------------------------------------------------------
889      !!                ***  ROUTINE lim_itd_me_icestrength ***
890      !!
891      !! ** Purpose :   shift ridging ice among thickness categories of ice thickness
892      !!
893      !! ** Method  :   Remove area, volume, and energy from each ridging category
894      !!              and add to thicker ice categories.
895      !!----------------------------------------------------------------------
896      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   opning         ! rate of opening due to divergence/shear
897      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area removed, excluding area of new ridges
898      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   msnow_mlt      ! mass of snow added to ocean (kg m-2)
899      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   esnow_mlt      ! energy needed to melt snow in ocean (J m-2)
900      !
901      CHARACTER (len=80) ::   fieldid   ! field identifier
902      LOGICAL, PARAMETER ::   l_conservation_check = .true.  ! if true, check conservation (useful for debugging)
903      !
904      LOGICAL ::   neg_ato_i      ! flag for ato_i(i,j) < -puny
905      LOGICAL ::   large_afrac    ! flag for afrac > 1
906      LOGICAL ::   large_afrft    ! flag for afrac > 1
907      INTEGER ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
908      INTEGER ::   ij                ! horizontal index, combines i and j loops
909      INTEGER ::   icells            ! number of cells with aicen > puny
910      REAL(wp) ::   zindb, zsrdg2   ! local scalar
911      REAL(wp) ::   hL, hR, farea, zdummy, zdummy0, ztmelts    ! left and right limits of integration
912
913      INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices
914
915      REAL(wp), POINTER, DIMENSION(:,:) ::   vice_init, vice_final   ! ice volume summed over categories
916      REAL(wp), POINTER, DIMENSION(:,:) ::   eice_init, eice_final   ! ice energy summed over layers
917
918      REAL(wp), POINTER, DIMENSION(:,:,:) ::   aicen_init, vicen_init   ! ice  area    & volume before ridging
919      REAL(wp), POINTER, DIMENSION(:,:,:) ::   vsnon_init, esnon_init   ! snow volume  & energy before ridging
920      REAL(wp), POINTER, DIMENSION(:,:,:) ::   smv_i_init, oa_i_init    ! ice salinity & age    before ridging
921
922      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   eicen_init        ! ice energy before ridging
923
924      REAL(wp), POINTER, DIMENSION(:,:) ::   afrac , fvol     ! fraction of category area ridged & new ridge volume going to n2
925      REAL(wp), POINTER, DIMENSION(:,:) ::   ardg1 , ardg2    ! area of ice ridged & new ridges
926      REAL(wp), POINTER, DIMENSION(:,:) ::   vsrdg , esrdg    ! snow volume & energy of ridging ice
927      REAL(wp), POINTER, DIMENSION(:,:) ::   oirdg1, oirdg2   ! areal age content of ridged & rifging ice
928      REAL(wp), POINTER, DIMENSION(:,:) ::   dhr   , dhr2     ! hrmax - hrmin  &  hrmax^2 - hrmin^2
929
930      REAL(wp), POINTER, DIMENSION(:,:) ::   vrdg1   ! volume of ice ridged
931      REAL(wp), POINTER, DIMENSION(:,:) ::   vrdg2   ! volume of new ridges
932      REAL(wp), POINTER, DIMENSION(:,:) ::   vsw     ! volume of seawater trapped into ridges
933      REAL(wp), POINTER, DIMENSION(:,:) ::   srdg1   ! sal*volume of ice ridged
934      REAL(wp), POINTER, DIMENSION(:,:) ::   srdg2   ! sal*volume of new ridges
935      REAL(wp), POINTER, DIMENSION(:,:) ::   smsw    ! sal*volume of water trapped into ridges
936
937      REAL(wp), POINTER, DIMENSION(:,:) ::   afrft            ! fraction of category area rafted
938      REAL(wp), POINTER, DIMENSION(:,:) ::   arft1 , arft2    ! area of ice rafted and new rafted zone
939      REAL(wp), POINTER, DIMENSION(:,:) ::   virft , vsrft    ! ice & snow volume of rafting ice
940      REAL(wp), POINTER, DIMENSION(:,:) ::   esrft , smrft    ! snow energy & salinity of rafting ice
941      REAL(wp), POINTER, DIMENSION(:,:) ::   oirft1, oirft2   ! areal age content of rafted ice & rafting ice
942
943      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eirft      ! ice energy of rafting ice
944      REAL(wp), POINTER, DIMENSION(:,:,:) ::   erdg1      ! enth*volume of ice ridged
945      REAL(wp), POINTER, DIMENSION(:,:,:) ::   erdg2      ! enth*volume of new ridges
946      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ersw       ! enth of water trapped into ridges
947      !!----------------------------------------------------------------------
948
949      CALL wrk_alloc( (jpi+1)*(jpj+1),      indxi, indxj )
950      CALL wrk_alloc( jpi, jpj,             vice_init, vice_final, eice_init, eice_final )
951      CALL wrk_alloc( jpi, jpj,             afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )
952      CALL wrk_alloc( jpi, jpj,             vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw )
953      CALL wrk_alloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
954      CALL wrk_alloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnon_init, esnon_init, smv_i_init, oa_i_init )
955      CALL wrk_alloc( jpi, jpj, jkmax,      eirft, erdg1, erdg2, ersw )
956      CALL wrk_alloc( jpi, jpj, jkmax, jpl, eicen_init )
957
958      ! Conservation check
959      eice_init(:,:) = 0._wp
960
961      IF( con_i ) THEN
962         CALL lim_column_sum        (jpl,    v_i,       vice_init )
963         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_init )
964         DO ji = mi0(jiindx), mi1(jiindx)
965            DO jj = mj0(jjindx), mj1(jjindx)
966               WRITE(numout,*) ' vice_init  : ', vice_init(ji,jj)
967               WRITE(numout,*) ' eice_init  : ', eice_init(ji,jj)
968            END DO
969         END DO
970      ENDIF
971
972      !-------------------------------------------------------------------------------
973      ! 1) Compute change in open water area due to closing and opening.
974      !-------------------------------------------------------------------------------
975
976      neg_ato_i = .false.
977
978      DO jj = 1, jpj
979         DO ji = 1, jpi
980            ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice        &
981               &                        + opning(ji,jj)                          * rdt_ice
982            IF( ato_i(ji,jj) < -epsi10 ) THEN
983               neg_ato_i = .TRUE.
984            ELSEIF( ato_i(ji,jj) < 0._wp ) THEN    ! roundoff error
985               ato_i(ji,jj) = 0._wp
986            ENDIF
987         END DO !jj
988      END DO !ji
989
990      ! if negative open water area alert it
991      IF( neg_ato_i ) THEN       ! there is a bug
992         DO jj = 1, jpj 
993            DO ji = 1, jpi
994               IF( ato_i(ji,jj) < -epsi10 ) THEN
995                  WRITE(numout,*) '' 
996                  WRITE(numout,*) 'Ridging error: ato_i < 0'
997                  WRITE(numout,*) 'ato_i : ', ato_i(ji,jj)
998               ENDIF               ! ato_i < -epsi10
999            END DO
1000         END DO
1001      ENDIF
1002
1003      !-----------------------------------------------------------------
1004      ! 2) Save initial state variables
1005      !-----------------------------------------------------------------
1006
1007      DO jl = 1, jpl
1008         aicen_init(:,:,jl) = a_i(:,:,jl)
1009         vicen_init(:,:,jl) = v_i(:,:,jl)
1010         vsnon_init(:,:,jl) = v_s(:,:,jl)
1011         !
1012         smv_i_init(:,:,jl) = smv_i(:,:,jl)
1013         oa_i_init (:,:,jl) = oa_i (:,:,jl)
1014      END DO !jl
1015
1016      esnon_init(:,:,:) = e_s(:,:,1,:)
1017
1018      DO jl = 1, jpl 
1019         DO jk = 1, nlay_i
1020            eicen_init(:,:,jk,jl) = e_i(:,:,jk,jl)
1021         END DO
1022      END DO
1023
1024      !
1025      !-----------------------------------------------------------------
1026      ! 3) Pump everything from ice which is being ridged / rafted
1027      !-----------------------------------------------------------------
1028      ! Compute the area, volume, and energy of ice ridging in each
1029      ! category, along with the area of the resulting ridge.
1030
1031      DO jl1 = 1, jpl !jl1 describes the ridging category
1032
1033         !------------------------------------------------
1034         ! 3.1) Identify grid cells with nonzero ridging
1035         !------------------------------------------------
1036
1037         icells = 0
1038         DO jj = 1, jpj
1039            DO ji = 1, jpi
1040               IF( aicen_init(ji,jj,jl1) > epsi10 .AND. athorn(ji,jj,jl1) > 0._wp  &
1041                  &   .AND. closing_gross(ji,jj) > 0._wp ) THEN
1042                  icells = icells + 1
1043                  indxi(icells) = ji
1044                  indxj(icells) = jj
1045               ENDIF ! test on a_icen_init
1046            END DO ! ji
1047         END DO ! jj
1048
1049         large_afrac = .false.
1050         large_afrft = .false.
1051
1052!CDIR NODEP
1053         DO ij = 1, icells
1054            ji = indxi(ij)
1055            jj = indxj(ij)
1056
1057            !--------------------------------------------------------------------
1058            ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
1059            !--------------------------------------------------------------------
1060
1061            ardg1(ji,jj) = aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice
1062            arft1(ji,jj) = araft (ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice
1063            ardg2(ji,jj) = ardg1(ji,jj) / krdg(ji,jj,jl1)
1064            arft2(ji,jj) = arft1(ji,jj) / kraft
1065
1066            oirdg1(ji,jj)= aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice
1067            oirft1(ji,jj)= araft (ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice
1068            oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1)
1069            oirft2(ji,jj)= oirft1(ji,jj) / kraft
1070
1071            !---------------------------------------------------------------
1072            ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1
1073            !---------------------------------------------------------------
1074
1075            afrac(ji,jj) = ardg1(ji,jj) / aicen_init(ji,jj,jl1) !ridging
1076            afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting
1077
1078            IF (afrac(ji,jj) > kamax + epsi10) THEN  !riging
1079               large_afrac = .true.
1080            ELSEIF (afrac(ji,jj) > kamax) THEN  ! roundoff error
1081               afrac(ji,jj) = kamax
1082            ENDIF
1083            IF (afrft(ji,jj) > kamax + epsi10) THEN !rafting
1084               large_afrft = .true.
1085            ELSEIF (afrft(ji,jj) > kamax) THEN  ! roundoff error
1086               afrft(ji,jj) = kamax
1087            ENDIF
1088
1089            !--------------------------------------------------------------------------
1090            ! 3.4) Subtract area, volume, and energy from ridging
1091            !     / rafting category n1.
1092            !--------------------------------------------------------------------------
1093            vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por )
1094            vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + ridge_por )
1095            vsw  (ji,jj) = vrdg1(ji,jj) * ridge_por
1096
1097            vsrdg(ji,jj) = vsnon_init(ji,jj,jl1) * afrac(ji,jj)
1098            esrdg(ji,jj) = esnon_init(ji,jj,jl1) * afrac(ji,jj)
1099            srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por )
1100            srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj)
1101
1102            ! rafting volumes, heat contents ...
1103            virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj)
1104            vsrft(ji,jj) = vsnon_init(ji,jj,jl1) * afrft(ji,jj)
1105            esrft(ji,jj) = esnon_init(ji,jj,jl1) * afrft(ji,jj)
1106            smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj) 
1107
1108            ! substract everything
1109            a_i(ji,jj,jl1)   = a_i(ji,jj,jl1)   - ardg1(ji,jj)  - arft1(ji,jj)
1110            v_i(ji,jj,jl1)   = v_i(ji,jj,jl1)   - vrdg1(ji,jj)  - virft(ji,jj)
1111            v_s(ji,jj,jl1)   = v_s(ji,jj,jl1)   - vsrdg(ji,jj)  - vsrft(ji,jj)
1112            e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg(ji,jj)  - esrft(ji,jj)
1113            oa_i(ji,jj,jl1)  = oa_i(ji,jj,jl1)  - oirdg1(ji,jj) - oirft1(ji,jj)
1114            smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1(ji,jj)  - smrft(ji,jj)
1115
1116            !-----------------------------------------------------------------
1117            ! 3.5) Compute properties of new ridges
1118            !-----------------------------------------------------------------
1119            !-------------
1120            ! Salinity
1121            !-------------
1122            smsw(ji,jj)  = sss_m(ji,jj) * vsw(ji,jj) * rhoic / rau0       ! salt content of seawater frozen in voids
1123
1124            zsrdg2       = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge
1125
1126            srdg2(ji,jj) = MIN( s_i_max * vrdg2(ji,jj) , zsrdg2 )         ! impose a maximum salinity
1127           
1128            !                                                             ! excess of salt is flushed into the ocean
1129            !sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice
1130
1131            !rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic    ! gurvan: increase in ice volume du to seawater frozen in voids             
1132
1133            !------------------------------------           
1134            ! 3.6 Increment ridging diagnostics
1135            !------------------------------------           
1136
1137            !        jl1 looping 1-jpl
1138            !           ij looping 1-icells
1139
1140            dardg1dt   (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj)
1141            dardg2dt   (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj)
1142            opening    (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice
1143
1144            IF( con_i )   vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj)
1145
1146            !------------------------------------------           
1147            ! 3.7 Put the snow somewhere in the ocean
1148            !------------------------------------------           
1149            !  Place part of the snow lost by ridging into the ocean.
1150            !  Note that esnow_mlt < 0; the ocean must cool to melt snow.
1151            !  If the ocean temp = Tf already, new ice must grow.
1152            !  During the next time step, thermo_rates will determine whether
1153            !  the ocean cools or new ice grows.
1154            !        jl1 looping 1-jpl
1155            !           ij looping 1-icells
1156
1157            msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosn*vsrdg(ji,jj)*(1.0-fsnowrdg)   &   ! rafting included
1158               &                                + rhosn*vsrft(ji,jj)*(1.0-fsnowrft)
1159
1160            esnow_mlt(ji,jj) = esnow_mlt(ji,jj) + esrdg(ji,jj)*(1.0-fsnowrdg)         &   !rafting included
1161               &                                + esrft(ji,jj)*(1.0-fsnowrft)         
1162
1163            !-----------------------------------------------------------------
1164            ! 3.8 Compute quantities used to apportion ice among categories
1165            ! in the n2 loop below
1166            !-----------------------------------------------------------------
1167
1168            !        jl1 looping 1-jpl
1169            !           ij looping 1-icells
1170
1171            dhr (ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1)
1172            dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1)
1173
1174         END DO                 ! ij
1175
1176         !--------------------------------------------------------------------
1177         ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
1178         !      compute ridged ice enthalpy
1179         !--------------------------------------------------------------------
1180         DO jk = 1, nlay_i
1181!CDIR NODEP
1182            DO ij = 1, icells
1183               ji = indxi(ij)
1184               jj = indxj(ij)
1185               ! heat content of ridged ice
1186               erdg1(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por ) 
1187               eirft(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj)
1188               e_i  (ji,jj,jk,jl1)  = e_i(ji,jj,jk,jl1) - erdg1(ji,jj,jk) - eirft(ji,jj,jk)
1189               ! sea water heat content
1190               ztmelts          = - tmut * sss_m(ji,jj) + rtt
1191               ! heat content per unit volume
1192               zdummy0          = - rcp * ( sst_m(ji,jj) + rt0 - rtt ) * vsw(ji,jj)
1193
1194               ! corrected sea water salinity
1195               zindb  = MAX( 0._wp , SIGN( 1._wp , vsw(ji,jj) - epsi20 ) )
1196               zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / MAX( ridge_por * vsw(ji,jj), epsi20 )
1197
1198               ztmelts          = - tmut * zdummy + rtt
1199               ersw(ji,jj,jk)   = - rcp * ( ztmelts - rtt ) * vsw(ji,jj)
1200
1201               ! heat flux
1202               fheat_mec(ji,jj) = fheat_mec(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) * r1_rdtice
1203
1204               ! Correct dimensions to avoid big values
1205               ersw(ji,jj,jk)   = ersw(ji,jj,jk) * 1.e-09
1206
1207               ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
1208               ersw (ji,jj,jk)  = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / REAL( nlay_i )
1209
1210               erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk)
1211            END DO ! ij
1212         END DO !jk
1213
1214
1215         IF( con_i ) THEN
1216            DO jk = 1, nlay_i
1217!CDIR NODEP
1218               DO ij = 1, icells
1219                  ji = indxi(ij)
1220                  jj = indxj(ij)
1221                  eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - erdg1(ji,jj,jk)
1222               END DO ! ij
1223            END DO !jk
1224         ENDIF
1225
1226         IF( large_afrac ) THEN   ! there is a bug
1227!CDIR NODEP
1228            DO ij = 1, icells
1229               ji = indxi(ij)
1230               jj = indxj(ij)
1231               IF( afrac(ji,jj) > kamax + epsi10 ) THEN
1232                  WRITE(numout,*) ''
1233                  WRITE(numout,*) ' ardg > a_i'
1234                  WRITE(numout,*) ' ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1)
1235               ENDIF
1236            END DO
1237         ENDIF
1238         IF( large_afrft ) THEN  ! there is a bug
1239!CDIR NODEP
1240            DO ij = 1, icells
1241               ji = indxi(ij)
1242               jj = indxj(ij)
1243               IF( afrft(ji,jj) > kamax + epsi10 ) THEN
1244                  WRITE(numout,*) ''
1245                  WRITE(numout,*) ' arft > a_i'
1246                  WRITE(numout,*) ' arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1)
1247               ENDIF
1248            END DO
1249         ENDIF
1250
1251         !-------------------------------------------------------------------------------
1252         ! 4) Add area, volume, and energy of new ridge to each category jl2
1253         !-------------------------------------------------------------------------------
1254         !        jl1 looping 1-jpl
1255         DO jl2  = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 
1256            ! over categories to which ridged ice is transferred
1257!CDIR NODEP
1258            DO ij = 1, icells
1259               ji = indxi(ij)
1260               jj = indxj(ij)
1261
1262               ! Compute the fraction of ridged ice area and volume going to
1263               ! thickness category jl2.
1264               ! Transfer area, volume, and energy accordingly.
1265
1266               IF( hrmin(ji,jj,jl1) >= hi_max(jl2)  .OR.        &
1267                   hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN
1268                  hL = 0._wp
1269                  hR = 0._wp
1270               ELSE
1271                  hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
1272                  hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2)   )
1273               ENDIF
1274
1275               ! fraction of ridged ice area and volume going to n2
1276               farea = ( hR - hL ) / dhr(ji,jj) 
1277               fvol(ji,jj) = ( hR*hR - hL*hL ) / dhr2(ji,jj)
1278
1279               a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ardg2 (ji,jj) * farea
1280               v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + vrdg2 (ji,jj) * fvol(ji,jj)
1281               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * fsnowrdg
1282               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * fsnowrdg
1283               smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + srdg2 (ji,jj) * fvol(ji,jj)
1284               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirdg2(ji,jj) * farea
1285
1286            END DO ! ij
1287
1288            ! Transfer ice energy to category jl2 by ridging
1289            DO jk = 1, nlay_i
1290!CDIR NODEP
1291               DO ij = 1, icells
1292                  ji = indxi(ij)
1293                  jj = indxj(ij)
1294                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj)*erdg2(ji,jj,jk)
1295               END DO
1296            END DO
1297            !
1298         END DO                 ! jl2 (new ridges)           
1299
1300         DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 
1301
1302!CDIR NODEP
1303            DO ij = 1, icells
1304               ji = indxi(ij)
1305               jj = indxj(ij)
1306               ! Compute the fraction of rafted ice area and volume going to
1307               ! thickness category jl2, transfer area, volume, and energy accordingly.
1308               !
1309               IF( hraft(ji,jj,jl1) <= hi_max(jl2)  .AND.        &
1310                   hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN
1311                  a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + arft2 (ji,jj)
1312                  v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + virft (ji,jj)
1313                  v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrft (ji,jj) * fsnowrft
1314                  e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrft (ji,jj) * fsnowrft
1315                  smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + smrft (ji,jj)   
1316                  oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirft2(ji,jj)   
1317               ENDIF ! hraft
1318               !
1319            END DO ! ij
1320
1321            ! Transfer rafted ice energy to category jl2
1322            DO jk = 1, nlay_i
1323!CDIR NODEP
1324               DO ij = 1, icells
1325                  ji = indxi(ij)
1326                  jj = indxj(ij)
1327                  IF(  hraft(ji,jj,jl1)  <=  hi_max(jl2)   .AND.        &
1328                       hraft(ji,jj,jl1)  >   hi_max(jl2-1)  ) THEN
1329                     e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk)
1330                  ENDIF
1331               END DO           ! ij
1332            END DO !jk
1333
1334         END DO ! jl2
1335
1336      END DO ! jl1 (deforming categories)
1337
1338      ! Conservation check
1339      IF ( con_i ) THEN
1340         CALL lim_column_sum (jpl,   v_i, vice_final)
1341         fieldid = ' v_i : limitd_me '
1342         CALL lim_cons_check (vice_init, vice_final, 1.0e-6, fieldid) 
1343
1344         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_final )
1345         fieldid = ' e_i : limitd_me '
1346         CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid) 
1347
1348         DO ji = mi0(jiindx), mi1(jiindx)
1349            DO jj = mj0(jjindx), mj1(jjindx)
1350               WRITE(numout,*) ' vice_init  : ', vice_init (ji,jj)
1351               WRITE(numout,*) ' vice_final : ', vice_final(ji,jj)
1352               WRITE(numout,*) ' eice_init  : ', eice_init (ji,jj)
1353               WRITE(numout,*) ' eice_final : ', eice_final(ji,jj)
1354            END DO
1355         END DO
1356      ENDIF
1357      !
1358      CALL wrk_dealloc( (jpi+1)*(jpj+1),      indxi, indxj )
1359      CALL wrk_dealloc( jpi, jpj,             vice_init, vice_final, eice_init, eice_final )
1360      CALL wrk_dealloc( jpi, jpj,             afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )
1361      CALL wrk_dealloc( jpi, jpj,             vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw )
1362      CALL wrk_dealloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
1363      CALL wrk_dealloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnon_init, esnon_init, smv_i_init, oa_i_init )
1364      CALL wrk_dealloc( jpi, jpj, jkmax,      eirft, erdg1, erdg2, ersw )
1365      CALL wrk_dealloc( jpi, jpj, jkmax, jpl, eicen_init )
1366      !
1367   END SUBROUTINE lim_itd_me_ridgeshift
1368
1369
1370   SUBROUTINE lim_itd_me_asumr
1371      !!-----------------------------------------------------------------------------
1372      !!                ***  ROUTINE lim_itd_me_asumr ***
1373      !!
1374      !! ** Purpose :   finds total fractional area
1375      !!
1376      !! ** Method  :   Find the total area of ice plus open water in each grid cell.
1377      !!              This is similar to the aggregate_area subroutine except that the
1378      !!              total area can be greater than 1, so the open water area is
1379      !!              included in the sum instead of being computed as a residual.
1380      !!-----------------------------------------------------------------------------
1381      INTEGER ::   jl   ! dummy loop index
1382      !!-----------------------------------------------------------------------------
1383      !
1384      asum(:,:) = ato_i(:,:)                    ! open water
1385      DO jl = 1, jpl                            ! ice categories
1386         asum(:,:) = asum(:,:) + a_i(:,:,jl)
1387      END DO
1388      !
1389   END SUBROUTINE lim_itd_me_asumr
1390
1391
1392   SUBROUTINE lim_itd_me_init
1393      !!-------------------------------------------------------------------
1394      !!                   ***  ROUTINE lim_itd_me_init ***
1395      !!
1396      !! ** Purpose :   Physical constants and parameters linked
1397      !!                to the mechanical ice redistribution
1398      !!
1399      !! ** Method  :   Read the namiceitdme namelist
1400      !!                and check the parameters values
1401      !!                called at the first timestep (nit000)
1402      !!
1403      !! ** input   :   Namelist namiceitdme
1404      !!-------------------------------------------------------------------
1405      INTEGER :: ios                 ! Local integer output status for namelist read
1406      NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft,& 
1407         Gstar, astar,                                &
1408         Hstar, raftswi, hparmeter, Craft, ridge_por, &
1409         sal_max_ridge,  partfun_swi, transfun_swi,   &
1410         brinstren_swi
1411      !!-------------------------------------------------------------------
1412      !
1413      REWIND( numnam_ice_ref )              ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
1414      READ  ( numnam_ice_ref, namiceitdme, IOSTAT = ios, ERR = 901)
1415901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in reference namelist', lwp )
1416
1417      REWIND( numnam_ice_cfg )              ! Namelist namiceitdme in configuration namelist : Ice mechanical ice redistribution
1418      READ  ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 )
1419902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in configuration namelist', lwp )
1420      WRITE ( numoni, namiceitdme )
1421      !
1422      IF (lwp) THEN                          ! control print
1423         WRITE(numout,*)
1424         WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution '
1425         WRITE(numout,*)' ~~~~~~~~~~~~~~~'
1426         WRITE(numout,*)'   Switch choosing the ice redistribution scheme           ridge_scheme_swi', ridge_scheme_swi 
1427         WRITE(numout,*)'   Fraction of shear energy contributing to ridging        Cs              ', Cs 
1428         WRITE(numout,*)'   Ratio of ridging work to PotEner change in ridging      Cf              ', Cf 
1429         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        fsnowrdg        ', fsnowrdg 
1430         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        fsnowrft        ', fsnowrft 
1431         WRITE(numout,*)'   Fraction of total ice coverage contributing to ridging  Gstar           ', Gstar
1432         WRITE(numout,*)'   Equivalent to G* for an exponential part function       astar           ', astar
1433         WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     Hstar           ', Hstar
1434         WRITE(numout,*)'   Rafting of ice sheets or not                            raftswi         ', raftswi
1435         WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       hparmeter       ', hparmeter
1436         WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  Craft           ', Craft 
1437         WRITE(numout,*)'   Initial porosity of ridges                              ridge_por       ', ridge_por
1438         WRITE(numout,*)'   Maximum salinity of ridging ice                         sal_max_ridge   ', sal_max_ridge
1439         WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    partfun_swi     ', partfun_swi
1440         WRITE(numout,*)'   Switch for tran. function (0) linear (1) exponential    transfun_swi    ', transfun_swi
1441         WRITE(numout,*)'   Switch for including brine volume in ice strength comp. brinstren_swi   ', brinstren_swi
1442      ENDIF
1443      !
1444   END SUBROUTINE lim_itd_me_init
1445
1446
1447   SUBROUTINE lim_itd_me_zapsmall
1448      !!-------------------------------------------------------------------
1449      !!                   ***  ROUTINE lim_itd_me_zapsmall ***
1450      !!
1451      !! ** Purpose :   Remove too small sea ice areas and correct salt fluxes
1452      !!
1453      !! history :
1454      !! author: William H. Lipscomb, LANL
1455      !! Nov 2003:  Modified by Julie Schramm to conserve volume and energy
1456      !! Sept 2004: Modified by William Lipscomb; replaced normalize_state with
1457      !!            additions to local freshwater, salt, and heat fluxes
1458      !!  9.0, LIM3.0 - 02-2006 (M. Vancoppenolle) original code
1459      !!-------------------------------------------------------------------
1460      INTEGER ::   ji, jj, jl, jk   ! dummy loop indices
1461      INTEGER ::   icells           ! number of cells with ice to zap
1462
1463      REAL(wp), POINTER, DIMENSION(:,:) ::   zmask   ! 2D workspace
1464      REAL(wp)                          ::   zmask_glo
1465!!gm      REAL(wp) ::   xtmp      ! temporary variable
1466      !!-------------------------------------------------------------------
1467
1468      CALL wrk_alloc( jpi, jpj, zmask )
1469
1470      DO jl = 1, jpl
1471
1472         !-----------------------------------------------------------------
1473         ! Count categories to be zapped.
1474         ! Abort model in case of negative area.
1475         !-----------------------------------------------------------------
1476         icells = 0
1477         zmask(:,:)  = 0._wp
1478         DO jj = 1, jpj
1479            DO ji = 1, jpi
1480               IF(  ( a_i(ji,jj,jl) >= -epsi10 .AND. a_i(ji,jj,jl) <  0._wp   ) .OR.   &
1481                  & ( a_i(ji,jj,jl) >  0._wp   .AND. a_i(ji,jj,jl) <= epsi10  ) .OR.   &
1482                  & ( v_i(ji,jj,jl) == 0._wp   .AND. a_i(ji,jj,jl) >  0._wp   ) .OR.   &
1483                  & ( v_i(ji,jj,jl) >  0._wp   .AND. v_i(ji,jj,jl) <= epsi10  ) )   zmask(ji,jj) = 1._wp
1484            END DO
1485         END DO
1486         zmask_glo = glob_sum(zmask)
1487         !IF( ln_nicep .AND. lwp ) WRITE(numout,*) zmask_glo, ' cells of ice zapped in the ocean '
1488
1489         !-----------------------------------------------------------------
1490         ! Zap ice energy and use ocean heat to melt ice
1491         !-----------------------------------------------------------------
1492
1493         DO jk = 1, nlay_i
1494            DO jj = 1 , jpj
1495               DO ji = 1 , jpi
1496!!gm                  xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) * r1_rdtice
1497!!gm                  xtmp = xtmp * unit_fac
1498                  ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp
1499                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) )
1500               END DO
1501            END DO
1502         END DO
1503
1504         DO jj = 1 , jpj
1505            DO ji = 1 , jpi
1506
1507               !-----------------------------------------------------------------
1508               ! Zap snow energy and use ocean heat to melt snow
1509               !-----------------------------------------------------------------
1510               !           xtmp = esnon(i,j,n) / dt ! < 0
1511               !           fhnet(i,j)      = fhnet(i,j)      + xtmp
1512               !           fhnet_hist(i,j) = fhnet_hist(i,j) + xtmp
1513               ! xtmp is greater than 0
1514               ! fluxes are positive to the ocean
1515               ! here the flux has to be negative for the ocean
1516!!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice
1517               !           fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp
1518
1519!!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice !RB   ???????
1520
1521               t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) )
1522
1523               !-----------------------------------------------------------------
1524               ! zap ice and snow volume, add water and salt to ocean
1525               !-----------------------------------------------------------------
1526
1527               !           xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) / dt
1528               !           sfx_res(ji,jj)  = sfx_res(ji,jj) + ( sss_m(ji,jj)                  )   &
1529               !                                            * rhosn * v_s(ji,jj,jl) * r1_rdtice
1530               !           sfx_res(ji,jj)  = sfx_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) )   &
1531               !                                            * rhoic * v_i(ji,jj,jl) * r1_rdtice
1532               !           sfx (i,j)      = sfx (i,j)      + xtmp
1533
1534               ato_i(ji,jj)    = a_i  (ji,jj,jl) *       zmask(ji,jj)   + ato_i(ji,jj)
1535               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
1536               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
1537               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
1538               t_su (ji,jj,jl) = t_su (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj)
1539               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
1540               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
1541               !
1542            END DO
1543         END DO
1544         !
1545      END DO                 ! jl
1546      !
1547      CALL wrk_dealloc( jpi, jpj, zmask )
1548      !
1549   END SUBROUTINE lim_itd_me_zapsmall
1550
1551#else
1552   !!----------------------------------------------------------------------
1553   !!   Default option         Empty module          NO LIM-3 sea-ice model
1554   !!----------------------------------------------------------------------
1555CONTAINS
1556   SUBROUTINE lim_itd_me           ! Empty routines
1557   END SUBROUTINE lim_itd_me
1558   SUBROUTINE lim_itd_me_icestrength
1559   END SUBROUTINE lim_itd_me_icestrength
1560   SUBROUTINE lim_itd_me_sort
1561   END SUBROUTINE lim_itd_me_sort
1562   SUBROUTINE lim_itd_me_init
1563   END SUBROUTINE lim_itd_me_init
1564   SUBROUTINE lim_itd_me_zapsmall
1565   END SUBROUTINE lim_itd_me_zapsmall
1566#endif
1567   !!======================================================================
1568END MODULE limitd_me
Note: See TracBrowser for help on using the repository browser.