source: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 4506

Last change on this file since 4506 was 4506, checked in by vancop, 7 years ago

[Heat conservation in LIM3, part HC1 (LIM_SRC_3_HC17)]

  • Property svn:keywords set to Id
File size: 72.8 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            ! there is a unit problem for e_snow_mlt that we should fix
366            fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice     ! heat sink for ocean
367
368         END DO
369      END DO
370
371      ! Check if there is a ridging error
372      DO jj = 1, jpj
373         DO ji = 1, jpi
374            IF( ABS( asum(ji,jj) - kamax)  >  epsi10 ) THEN   ! there is a bug
375               WRITE(numout,*) ' '
376               WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj)
377               WRITE(numout,*) ' limitd_me '
378               WRITE(numout,*) ' POINT : ', ji, jj
379               WRITE(numout,*) ' jpl, a_i, athorn '
380               WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0)
381               DO jl = 1, jpl
382                  WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl)
383               END DO
384            ENDIF  ! asum
385
386         END DO !ji
387      END DO !jj
388
389      ! Conservation check
390      IF ( con_i ) THEN
391         CALL lim_column_sum (jpl,   v_i, vt_i_final)
392         fieldid = ' v_i : limitd_me '
393         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
394      ENDIF
395
396      !-----------------------------------------------------------------------------!
397      ! 6) Updating state variables and trend terms (done in limupdate)
398      !-----------------------------------------------------------------------------!
399      CALL lim_var_glo2eqv
400      CALL lim_itd_me_zapsmall
401
402      !--------------------------------
403      ! Update mass/salt fluxes (clem)
404      !--------------------------------
405      DO jl = 1, jpl
406         DO jj = 1, jpj 
407            DO ji = 1, jpi
408               diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice
409               rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic 
410               rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn 
411               sfx_mec(ji,jj) = sfx_mec(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic * r1_rdtice 
412            END DO
413         END DO
414      END DO
415
416      IF(ln_ctl) THEN     ! Control print
417         CALL prt_ctl_info(' ')
418         CALL prt_ctl_info(' - Cell values : ')
419         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
420         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_me  : cell area :')
421         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me  : at_i      :')
422         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me  : vt_i      :')
423         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_me  : vt_s      :')
424         DO jl = 1, jpl
425            CALL prt_ctl_info(' ')
426            CALL prt_ctl_info(' - Category : ', ivar1=jl)
427            CALL prt_ctl_info('   ~~~~~~~~~~')
428            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_itd_me  : a_i      : ')
429            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_itd_me  : ht_i     : ')
430            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_itd_me  : ht_s     : ')
431            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_itd_me  : v_i      : ')
432            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_itd_me  : v_s      : ')
433            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_itd_me  : e_s      : ')
434            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_itd_me  : t_su     : ')
435            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_itd_me  : t_snow   : ')
436            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_itd_me  : sm_i     : ')
437            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_itd_me  : smv_i    : ')
438            DO jk = 1, nlay_i
439               CALL prt_ctl_info(' ')
440               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
441               CALL prt_ctl_info('   ~~~~~~~')
442               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_me  : t_i      : ')
443               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_me  : e_i      : ')
444            END DO
445         END DO
446      ENDIF
447
448      ! -------------------------------
449      !- check conservation (C Rousset)
450      IF (ln_limdiahsb) THEN
451         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
452         zchk_fw  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b
453 
454         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice
455         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic )
456
457         zchk_vmin = glob_min(v_i)
458         zchk_amax = glob_max(SUM(a_i,dim=3))
459         zchk_amin = glob_min(a_i)
460       
461         IF(lwp) THEN
462            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limitd_me) = ',(zchk_v_i * rday)
463            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * rday)
464            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_me) = ',(zchk_vmin * 1.e-3)
465            IF ( zchk_amax >  kamax+epsi10  ) WRITE(numout,*) 'violation a_i>amax            (limitd_me) = ',zchk_amax
466            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limitd_me) = ',zchk_amin
467         ENDIF
468      ENDIF
469      !- check conservation (C Rousset)
470      ! -------------------------------
471
472      ENDIF  ! ln_limdyn=.true.
473      !
474      CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
475      !
476      CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zsmvold )   ! clem
477      !
478      IF( nn_timing == 1 )  CALL timing_stop('limitd_me')
479   END SUBROUTINE lim_itd_me
480
481
482   SUBROUTINE lim_itd_me_icestrength( kstrngth )
483      !!----------------------------------------------------------------------
484      !!                ***  ROUTINE lim_itd_me_icestrength ***
485      !!
486      !! ** Purpose :   computes ice strength used in dynamics routines of ice thickness
487      !!
488      !! ** Method  :   Compute the strength of the ice pack, defined as the energy (J m-2)
489      !!              dissipated per unit area removed from the ice pack under compression,
490      !!              and assumed proportional to the change in potential energy caused
491      !!              by ridging. Note that only Hibler's formulation is stable and that
492      !!              ice strength has to be smoothed
493      !!
494      !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
495      !!----------------------------------------------------------------------
496      INTEGER, INTENT(in) ::   kstrngth    ! = 1 for Rothrock formulation, 0 for Hibler (1979)
497
498      INTEGER ::   ji,jj, jl   ! dummy loop indices
499      INTEGER ::   ksmooth     ! smoothing the resistance to deformation
500      INTEGER ::   numts_rm    ! number of time steps for the P smoothing
501      REAL(wp) ::   hi, zw1, zp, zdummy, zzc, z1_3   ! local scalars
502      REAL(wp), POINTER, DIMENSION(:,:) ::   zworka   ! temporary array used here
503      !!----------------------------------------------------------------------
504
505      CALL wrk_alloc( jpi, jpj, zworka )
506
507      !------------------------------------------------------------------------------!
508      ! 1) Initialize
509      !------------------------------------------------------------------------------!
510      strength(:,:) = 0._wp
511
512      !------------------------------------------------------------------------------!
513      ! 2) Compute thickness distribution of ridging and ridged ice
514      !------------------------------------------------------------------------------!
515      CALL lim_itd_me_ridgeprep
516
517      !------------------------------------------------------------------------------!
518      ! 3) Rothrock(1975)'s method
519      !------------------------------------------------------------------------------!
520      IF( kstrngth == 1 ) THEN
521         z1_3 = 1._wp / 3._wp
522         DO jl = 1, jpl
523            DO jj= 1, jpj
524               DO ji = 1, jpi
525                  !
526                  IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN
527                     hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
528                     !----------------------------
529                     ! PE loss from deforming ice
530                     !----------------------------
531                     strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * hi * hi
532
533                     !--------------------------
534                     ! PE gain from rafting ice
535                     !--------------------------
536                     strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * hi * hi
537
538                     !----------------------------
539                     ! PE gain from ridging ice
540                     !----------------------------
541                     strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl)/krdg(ji,jj,jl)     &
542                        * z1_3 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) / ( hrmax(ji,jj,jl)-hrmin(ji,jj,jl) )   
543!!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...                   
544                  ENDIF            ! aicen > epsi10
545                  !
546               END DO ! ji
547            END DO !jj
548         END DO !jl
549
550         zzc = Cf * Cp     ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and Cf accounts for frictional dissipation
551         strength(:,:) = zzc * strength(:,:) / aksum(:,:)
552
553         ksmooth = 1
554
555         !------------------------------------------------------------------------------!
556         ! 4) Hibler (1979)' method
557         !------------------------------------------------------------------------------!
558      ELSE                      ! kstrngth ne 1:  Hibler (1979) form
559         !
560         strength(:,:) = Pstar * vt_i(:,:) * EXP( - C_rhg * ( 1._wp - at_i(:,:) )  )
561         !
562         ksmooth = 1
563         !
564      ENDIF                     ! kstrngth
565
566      !
567      !------------------------------------------------------------------------------!
568      ! 5) Impact of brine volume
569      !------------------------------------------------------------------------------!
570      ! CAN BE REMOVED
571      !
572      IF ( brinstren_swi == 1 ) THEN
573
574         DO jj = 1, jpj
575            DO ji = 1, jpi
576               IF ( bv_i(ji,jj) .GT. 0.0 ) THEN
577                  zdummy = MIN ( bv_i(ji,jj), 0.10 ) * MIN( bv_i(ji,jj), 0.10 )
578               ELSE
579                  zdummy = 0.0
580               ENDIF
581               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0)))
582            END DO              ! j
583         END DO                 ! i
584
585      ENDIF
586
587      !
588      !------------------------------------------------------------------------------!
589      ! 6) Smoothing ice strength
590      !------------------------------------------------------------------------------!
591      !
592      !-------------------
593      ! Spatial smoothing
594      !-------------------
595      IF ( ksmooth == 1 ) THEN
596
597         CALL lbc_lnk( strength, 'T', 1. )
598
599         DO jj = 2, jpj - 1
600            DO ji = 2, jpi - 1
601               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN ! ice is
602                  ! present
603                  zworka(ji,jj) = 4.0 * strength(ji,jj)              &
604                     &          + strength(ji-1,jj) * tms(ji-1,jj) & 
605                     &          + strength(ji+1,jj) * tms(ji+1,jj) & 
606                     &          + strength(ji,jj-1) * tms(ji,jj-1) & 
607                     &          + strength(ji,jj+1) * tms(ji,jj+1)   
608
609                  zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj) + tms(ji,jj-1) + tms(ji,jj+1)
610                  zworka(ji,jj) = zworka(ji,jj) / zw1
611               ELSE
612                  zworka(ji,jj) = 0._wp
613               ENDIF
614            END DO
615         END DO
616
617         DO jj = 2, jpj - 1
618            DO ji = 2, jpi - 1
619               strength(ji,jj) = zworka(ji,jj)
620            END DO
621         END DO
622         CALL lbc_lnk( strength, 'T', 1. )
623
624      ENDIF ! ksmooth
625
626      !--------------------
627      ! Temporal smoothing
628      !--------------------
629      IF ( numit == nit000 + nn_fsbc - 1 ) THEN
630         strp1(:,:) = 0.0           
631         strp2(:,:) = 0.0           
632      ENDIF
633
634      IF ( ksmooth == 2 ) THEN
635
636
637         CALL lbc_lnk( strength, 'T', 1. )
638
639         DO jj = 1, jpj - 1
640            DO ji = 1, jpi - 1
641               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN       ! ice is present
642                  numts_rm = 1 ! number of time steps for the running mean
643                  IF ( strp1(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1
644                  IF ( strp2(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1
645                  zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm
646                  strp2(ji,jj) = strp1(ji,jj)
647                  strp1(ji,jj) = strength(ji,jj)
648                  strength(ji,jj) = zp
649
650               ENDIF
651            END DO
652         END DO
653
654      ENDIF ! ksmooth
655
656      CALL lbc_lnk( strength, 'T', 1. )      ! Boundary conditions
657
658      CALL wrk_dealloc( jpi, jpj, zworka )
659      !
660   END SUBROUTINE lim_itd_me_icestrength
661
662
663   SUBROUTINE lim_itd_me_ridgeprep
664      !!---------------------------------------------------------------------!
665      !!                ***  ROUTINE lim_itd_me_ridgeprep ***
666      !!
667      !! ** Purpose :   preparation for ridging and strength calculations
668      !!
669      !! ** Method  :   Compute the thickness distribution of the ice and open water
670      !!              participating in ridging and of the resulting ridges.
671      !!---------------------------------------------------------------------!
672      INTEGER ::   ji,jj, jl    ! dummy loop indices
673      INTEGER ::   krdg_index   !
674      REAL(wp) ::   Gstari, astari, hi, hrmean, zdummy   ! local scalar
675      REAL(wp), POINTER, DIMENSION(:,:)   ::   zworka    ! temporary array used here
676      REAL(wp), POINTER, DIMENSION(:,:,:) ::   Gsum      ! Gsum(n) = sum of areas in categories 0 to n
677      !------------------------------------------------------------------------------!
678
679      CALL wrk_alloc( jpi,jpj, zworka )
680      CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
681
682      Gstari     = 1.0/Gstar   
683      astari     = 1.0/astar   
684      aksum(:,:)    = 0.0
685      athorn(:,:,:) = 0.0
686      aridge(:,:,:) = 0.0
687      araft (:,:,:) = 0.0
688      hrmin(:,:,:)  = 0.0 
689      hrmax(:,:,:)  = 0.0 
690      hraft(:,:,:)  = 0.0 
691      krdg (:,:,:)  = 1.0
692
693      !     ! Zero out categories with very small areas
694      CALL lim_itd_me_zapsmall
695
696      !------------------------------------------------------------------------------!
697      ! 1) Participation function
698      !------------------------------------------------------------------------------!
699
700      ! Compute total area of ice plus open water.
701      CALL lim_itd_me_asumr
702      ! This is in general not equal to one
703      ! because of divergence during transport
704
705      ! Compute cumulative thickness distribution function
706      ! Compute the cumulative thickness distribution function Gsum,
707      ! where Gsum(n) is the fractional area in categories 0 to n.
708      ! initial value (in h = 0) equals open water area
709
710      Gsum(:,:,-1) = 0._wp
711
712      DO jj = 1, jpj
713         DO ji = 1, jpi
714            IF( ato_i(ji,jj) > epsi10 ) THEN   ;   Gsum(ji,jj,0) = ato_i(ji,jj)
715            ELSE                               ;   Gsum(ji,jj,0) = 0._wp
716            ENDIF
717         END DO
718      END DO
719
720      ! for each value of h, you have to add ice concentration then
721      DO jl = 1, jpl
722         DO jj = 1, jpj 
723            DO ji = 1, jpi
724               IF( a_i(ji,jj,jl) .GT. epsi10 ) THEN   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl)
725               ELSE                                   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1)
726               ENDIF
727            END DO
728         END DO
729      END DO
730
731      ! Normalize the cumulative distribution to 1
732      zworka(:,:) = 1._wp / Gsum(:,:,jpl)
733      DO jl = 0, jpl
734         Gsum(:,:,jl) = Gsum(:,:,jl) * zworka(:,:)
735      END DO
736
737      ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
738      !--------------------------------------------------------------------------------------------------
739      ! Compute the participation function athorn; this is analogous to
740      ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
741      ! area lost from category n due to ridging/closing
742      ! athorn(n)   = total area lost due to ridging/closing
743      ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
744      !
745      ! The expressions for athorn are found by integrating b(h)g(h) between
746      ! the category boundaries.
747      !-----------------------------------------------------------------
748
749      krdg_index = 1
750
751      IF( krdg_index == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975)
752         DO jl = 0, ice_cat_bounds(1,2)       ! only undeformed ice participates
753            DO jj = 1, jpj 
754               DO ji = 1, jpi
755                  IF( Gsum(ji,jj,jl) < Gstar) THEN
756                     athorn(ji,jj,jl) = Gstari * (Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * &
757                        (2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari)
758                  ELSEIF (Gsum(ji,jj,jl-1) < Gstar) THEN
759                     athorn(ji,jj,jl) = Gstari * (Gstar-Gsum(ji,jj,jl-1)) *  &
760                        (2.0 - (Gsum(ji,jj,jl-1)+Gstar)*Gstari)
761                  ELSE
762                     athorn(ji,jj,jl) = 0.0
763                  ENDIF
764               END DO ! ji
765            END DO ! jj
766         END DO ! jl
767
768      ELSE                             !--- Exponential, more stable formulation (Lipscomb et al, 2007)
769         !                       
770         zdummy = 1._wp / ( 1._wp - EXP(-astari) )        ! precompute exponential terms using Gsum as a work array
771
772         DO jl = -1, jpl
773            Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
774         END DO !jl
775         DO jl = 0, ice_cat_bounds(1,2)
776             athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
777         END DO
778         !
779      ENDIF ! krdg_index
780
781      IF( raftswi == 1 ) THEN      ! Ridging and rafting ice participation functions
782         !
783         DO jl = 1, jpl
784            DO jj = 1, jpj 
785               DO ji = 1, jpi
786                  IF ( athorn(ji,jj,jl) .GT. 0._wp ) THEN
787!!gm  TANH( -X ) = - TANH( X )  so can be computed only 1 time....
788                     aridge(ji,jj,jl) = ( TANH (  Craft * ( ht_i(ji,jj,jl) - hparmeter ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
789                     araft (ji,jj,jl) = ( TANH ( -Craft * ( ht_i(ji,jj,jl) - hparmeter ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
790                     IF ( araft(ji,jj,jl) < epsi06 )   araft(ji,jj,jl)  = 0._wp
791                     aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 )
792                  ENDIF ! athorn
793               END DO ! ji
794            END DO ! jj
795         END DO ! jl
796
797      ELSE  ! raftswi = 0
798         !
799         DO jl = 1, jpl
800            aridge(:,:,jl) = athorn(:,:,jl)
801         END DO
802         !
803      ENDIF
804
805      IF ( raftswi == 1 ) THEN
806
807         IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi10 ) THEN
808            DO jl = 1, jpl
809               DO jj = 1, jpj
810                  DO ji = 1, jpi
811                     IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. epsi10 ) THEN
812                        WRITE(numout,*) ' ALERTE 96 : wrong participation function ... '
813                        WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl
814                        WRITE(numout,*) ' lat, lon   : ', gphit(ji,jj), glamt(ji,jj)
815                        WRITE(numout,*) ' aridge     : ', aridge(ji,jj,1:jpl)
816                        WRITE(numout,*) ' araft      : ', araft(ji,jj,1:jpl)
817                        WRITE(numout,*) ' athorn     : ', athorn(ji,jj,1:jpl)
818                     ENDIF
819                  END DO
820               END DO
821            END DO
822         ENDIF
823
824      ENDIF
825
826      !-----------------------------------------------------------------
827      ! 2) Transfer function
828      !-----------------------------------------------------------------
829      ! Compute max and min ridged ice thickness for each ridging category.
830      ! Assume ridged ice is uniformly distributed between hrmin and hrmax.
831      !
832      ! This parameterization is a modified version of Hibler (1980).
833      ! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
834      !  and for very thick ridging ice must be >= krdgmin*hi
835      !
836      ! The minimum ridging thickness, hrmin, is equal to 2*hi
837      !  (i.e., rafting) and for very thick ridging ice is
838      !  constrained by hrmin <= (hrmean + hi)/2.
839      !
840      ! The maximum ridging thickness, hrmax, is determined by
841      !  hrmean and hrmin.
842      !
843      ! These modifications have the effect of reducing the ice strength
844      ! (relative to the Hibler formulation) when very thick ice is
845      ! ridging.
846      !
847      ! aksum = net area removed/ total area removed
848      ! where total area removed = area of ice that ridges
849      !         net area removed = total area removed - area of new ridges
850      !-----------------------------------------------------------------
851
852      ! Transfer function
853      DO jl = 1, jpl !all categories have a specific transfer function
854         DO jj = 1, jpj
855            DO ji = 1, jpi
856
857               IF (a_i(ji,jj,jl) .GT. epsi10 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN
858                  hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
859                  hrmean          = MAX(SQRT(Hstar*hi), hi*krdgmin)
860                  hrmin(ji,jj,jl) = MIN(2.0*hi, 0.5*(hrmean + hi))
861                  hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl)
862                  hraft(ji,jj,jl) = kraft*hi
863                  krdg(ji,jj,jl)  = hrmean / hi
864               ELSE
865                  hraft(ji,jj,jl) = 0.0
866                  hrmin(ji,jj,jl) = 0.0 
867                  hrmax(ji,jj,jl) = 0.0 
868                  krdg (ji,jj,jl) = 1.0
869               ENDIF
870
871            END DO ! ji
872         END DO ! jj
873      END DO ! jl
874
875      ! Normalization factor : aksum, ensures mass conservation
876      aksum(:,:) = athorn(:,:,0)
877      DO jl = 1, jpl 
878         aksum(:,:)    = aksum(:,:) + aridge(:,:,jl) * ( 1._wp - 1._wp / krdg(:,:,jl) )    &
879            &                       + araft (:,:,jl) * ( 1._wp - 1._wp / kraft        )
880      END DO
881      !
882      CALL wrk_dealloc( jpi,jpj, zworka )
883      CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
884      !
885   END SUBROUTINE lim_itd_me_ridgeprep
886
887
888   SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
889      !!----------------------------------------------------------------------
890      !!                ***  ROUTINE lim_itd_me_icestrength ***
891      !!
892      !! ** Purpose :   shift ridging ice among thickness categories of ice thickness
893      !!
894      !! ** Method  :   Remove area, volume, and energy from each ridging category
895      !!              and add to thicker ice categories.
896      !!----------------------------------------------------------------------
897      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   opning         ! rate of opening due to divergence/shear
898      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area removed, excluding area of new ridges
899      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   msnow_mlt      ! mass of snow added to ocean (kg m-2)
900      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   esnow_mlt      ! energy needed to melt snow in ocean (J m-2)
901      !
902      CHARACTER (len=80) ::   fieldid   ! field identifier
903      LOGICAL, PARAMETER ::   l_conservation_check = .true.  ! if true, check conservation (useful for debugging)
904      !
905      LOGICAL ::   neg_ato_i      ! flag for ato_i(i,j) < -puny
906      LOGICAL ::   large_afrac    ! flag for afrac > 1
907      LOGICAL ::   large_afrft    ! flag for afrac > 1
908      INTEGER ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
909      INTEGER ::   ij                ! horizontal index, combines i and j loops
910      INTEGER ::   icells            ! number of cells with aicen > puny
911      REAL(wp) ::   zindb, zsrdg2    ! local scalar
912      REAL(wp) ::   hL, hR, farea, zdummy, zdummy0, ztmelts    ! left and right limits of integration
913      REAL(wp) ::   zsstK            ! SST in Kelvin
914
915      INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices
916
917      REAL(wp), POINTER, DIMENSION(:,:) ::   vice_init, vice_final   ! ice volume summed over categories
918      REAL(wp), POINTER, DIMENSION(:,:) ::   eice_init, eice_final   ! ice energy summed over layers
919
920      REAL(wp), POINTER, DIMENSION(:,:,:) ::   aicen_init, vicen_init   ! ice  area    & volume before ridging
921      REAL(wp), POINTER, DIMENSION(:,:,:) ::   vsnon_init, esnon_init   ! snow volume  & energy before ridging
922      REAL(wp), POINTER, DIMENSION(:,:,:) ::   smv_i_init, oa_i_init    ! ice salinity & age    before ridging
923
924      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   eicen_init        ! ice energy before ridging
925
926      REAL(wp), POINTER, DIMENSION(:,:) ::   afrac , fvol     ! fraction of category area ridged & new ridge volume going to n2
927      REAL(wp), POINTER, DIMENSION(:,:) ::   ardg1 , ardg2    ! area of ice ridged & new ridges
928      REAL(wp), POINTER, DIMENSION(:,:) ::   vsrdg , esrdg    ! snow volume & energy of ridging ice
929      REAL(wp), POINTER, DIMENSION(:,:) ::   oirdg1, oirdg2   ! areal age content of ridged & rifging ice
930      REAL(wp), POINTER, DIMENSION(:,:) ::   dhr   , dhr2     ! hrmax - hrmin  &  hrmax^2 - hrmin^2
931
932      REAL(wp), POINTER, DIMENSION(:,:) ::   vrdg1   ! volume of ice ridged
933      REAL(wp), POINTER, DIMENSION(:,:) ::   vrdg2   ! volume of new ridges
934      REAL(wp), POINTER, DIMENSION(:,:) ::   vsw     ! volume of seawater trapped into ridges
935      REAL(wp), POINTER, DIMENSION(:,:) ::   srdg1   ! sal*volume of ice ridged
936      REAL(wp), POINTER, DIMENSION(:,:) ::   srdg2   ! sal*volume of new ridges
937      REAL(wp), POINTER, DIMENSION(:,:) ::   smsw    ! sal*volume of water trapped into ridges
938
939      REAL(wp), POINTER, DIMENSION(:,:) ::   afrft            ! fraction of category area rafted
940      REAL(wp), POINTER, DIMENSION(:,:) ::   arft1 , arft2    ! area of ice rafted and new rafted zone
941      REAL(wp), POINTER, DIMENSION(:,:) ::   virft , vsrft    ! ice & snow volume of rafting ice
942      REAL(wp), POINTER, DIMENSION(:,:) ::   esrft , smrft    ! snow energy & salinity of rafting ice
943      REAL(wp), POINTER, DIMENSION(:,:) ::   oirft1, oirft2   ! areal age content of rafted ice & rafting ice
944
945      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eirft      ! ice energy of rafting ice
946      REAL(wp), POINTER, DIMENSION(:,:,:) ::   erdg1      ! enth*volume of ice ridged
947      REAL(wp), POINTER, DIMENSION(:,:,:) ::   erdg2      ! enth*volume of new ridges
948      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ersw       ! enth of water trapped into ridges
949      !!----------------------------------------------------------------------
950
951      CALL wrk_alloc( (jpi+1)*(jpj+1),      indxi, indxj )
952      CALL wrk_alloc( jpi, jpj,             vice_init, vice_final, eice_init, eice_final )
953      CALL wrk_alloc( jpi, jpj,             afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )
954      CALL wrk_alloc( jpi, jpj,             vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw )
955      CALL wrk_alloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
956      CALL wrk_alloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnon_init, esnon_init, smv_i_init, oa_i_init )
957      CALL wrk_alloc( jpi, jpj, jkmax,      eirft, erdg1, erdg2, ersw )
958      CALL wrk_alloc( jpi, jpj, jkmax, jpl, eicen_init )
959
960      ! Conservation check
961      eice_init(:,:) = 0._wp
962
963      IF( con_i ) THEN
964         CALL lim_column_sum        (jpl,    v_i,       vice_init )
965         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_init )
966         DO ji = mi0(jiindx), mi1(jiindx)
967            DO jj = mj0(jjindx), mj1(jjindx)
968               WRITE(numout,*) ' vice_init  : ', vice_init(ji,jj)
969               WRITE(numout,*) ' eice_init  : ', eice_init(ji,jj)
970            END DO
971         END DO
972      ENDIF
973
974      !-------------------------------------------------------------------------------
975      ! 1) Compute change in open water area due to closing and opening.
976      !-------------------------------------------------------------------------------
977
978      neg_ato_i = .false.
979
980      DO jj = 1, jpj
981         DO ji = 1, jpi
982            ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice        &
983               &                        + opning(ji,jj)                          * rdt_ice
984            IF( ato_i(ji,jj) < -epsi10 ) THEN
985               neg_ato_i = .TRUE.
986            ELSEIF( ato_i(ji,jj) < 0._wp ) THEN    ! roundoff error
987               ato_i(ji,jj) = 0._wp
988            ENDIF
989         END DO !jj
990      END DO !ji
991
992      ! if negative open water area alert it
993      IF( neg_ato_i ) THEN       ! there is a bug
994         DO jj = 1, jpj 
995            DO ji = 1, jpi
996               IF( ato_i(ji,jj) < -epsi10 ) THEN
997                  WRITE(numout,*) '' 
998                  WRITE(numout,*) 'Ridging error: ato_i < 0'
999                  WRITE(numout,*) 'ato_i : ', ato_i(ji,jj)
1000               ENDIF               ! ato_i < -epsi10
1001            END DO
1002         END DO
1003      ENDIF
1004
1005      !-----------------------------------------------------------------
1006      ! 2) Save initial state variables
1007      !-----------------------------------------------------------------
1008
1009      DO jl = 1, jpl
1010         aicen_init(:,:,jl) = a_i(:,:,jl)
1011         vicen_init(:,:,jl) = v_i(:,:,jl)
1012         vsnon_init(:,:,jl) = v_s(:,:,jl)
1013         !
1014         smv_i_init(:,:,jl) = smv_i(:,:,jl)
1015         oa_i_init (:,:,jl) = oa_i (:,:,jl)
1016      END DO !jl
1017
1018      esnon_init(:,:,:) = e_s(:,:,1,:)
1019
1020      DO jl = 1, jpl 
1021         DO jk = 1, nlay_i
1022            eicen_init(:,:,jk,jl) = e_i(:,:,jk,jl)
1023         END DO
1024      END DO
1025
1026      !
1027      !-----------------------------------------------------------------
1028      ! 3) Pump everything from ice which is being ridged / rafted
1029      !-----------------------------------------------------------------
1030      ! Compute the area, volume, and energy of ice ridging in each
1031      ! category, along with the area of the resulting ridge.
1032
1033      DO jl1 = 1, jpl !jl1 describes the ridging category
1034
1035         !------------------------------------------------
1036         ! 3.1) Identify grid cells with nonzero ridging
1037         !------------------------------------------------
1038
1039         icells = 0
1040         DO jj = 1, jpj
1041            DO ji = 1, jpi
1042               IF( aicen_init(ji,jj,jl1) > epsi10 .AND. athorn(ji,jj,jl1) > 0._wp  &
1043                  &   .AND. closing_gross(ji,jj) > 0._wp ) THEN
1044                  icells = icells + 1
1045                  indxi(icells) = ji
1046                  indxj(icells) = jj
1047               ENDIF ! test on a_icen_init
1048            END DO ! ji
1049         END DO ! jj
1050
1051         large_afrac = .false.
1052         large_afrft = .false.
1053
1054!CDIR NODEP
1055         DO ij = 1, icells
1056            ji = indxi(ij)
1057            jj = indxj(ij)
1058
1059            !--------------------------------------------------------------------
1060            ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
1061            !--------------------------------------------------------------------
1062
1063            ardg1(ji,jj) = aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice
1064            arft1(ji,jj) = araft (ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice
1065            ardg2(ji,jj) = ardg1(ji,jj) / krdg(ji,jj,jl1)
1066            arft2(ji,jj) = arft1(ji,jj) / kraft
1067
1068            oirdg1(ji,jj)= aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice
1069            oirft1(ji,jj)= araft (ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice
1070            oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1)
1071            oirft2(ji,jj)= oirft1(ji,jj) / kraft
1072
1073            !---------------------------------------------------------------
1074            ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1
1075            !---------------------------------------------------------------
1076
1077            afrac(ji,jj) = ardg1(ji,jj) / aicen_init(ji,jj,jl1) !ridging
1078            afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting
1079
1080            IF (afrac(ji,jj) > kamax + epsi10) THEN  !riging
1081               large_afrac = .true.
1082            ELSEIF (afrac(ji,jj) > kamax) THEN  ! roundoff error
1083               afrac(ji,jj) = kamax
1084            ENDIF
1085            IF (afrft(ji,jj) > kamax + epsi10) THEN !rafting
1086               large_afrft = .true.
1087            ELSEIF (afrft(ji,jj) > kamax) THEN  ! roundoff error
1088               afrft(ji,jj) = kamax
1089            ENDIF
1090
1091            !--------------------------------------------------------------------------
1092            ! 3.4) Subtract area, volume, and energy from ridging
1093            !     / rafting category n1.
1094            !--------------------------------------------------------------------------
1095            vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por )
1096            vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + ridge_por )
1097            vsw  (ji,jj) = vrdg1(ji,jj) * ridge_por
1098
1099            vsrdg(ji,jj) = vsnon_init(ji,jj,jl1) * afrac(ji,jj)
1100            esrdg(ji,jj) = esnon_init(ji,jj,jl1) * afrac(ji,jj)
1101            srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por )
1102            srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) !! MV HC 2014 this line seems useless
1103
1104            ! rafting volumes, heat contents ...
1105            virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj)
1106            vsrft(ji,jj) = vsnon_init(ji,jj,jl1) * afrft(ji,jj)
1107            esrft(ji,jj) = esnon_init(ji,jj,jl1) * afrft(ji,jj)
1108            smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj) 
1109
1110            ! substract everything
1111            a_i(ji,jj,jl1)   = a_i(ji,jj,jl1)   - ardg1(ji,jj)  - arft1(ji,jj)
1112            v_i(ji,jj,jl1)   = v_i(ji,jj,jl1)   - vrdg1(ji,jj)  - virft(ji,jj)
1113            v_s(ji,jj,jl1)   = v_s(ji,jj,jl1)   - vsrdg(ji,jj)  - vsrft(ji,jj)
1114            e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg(ji,jj)  - esrft(ji,jj)
1115            oa_i(ji,jj,jl1)  = oa_i(ji,jj,jl1)  - oirdg1(ji,jj) - oirft1(ji,jj)
1116            smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1(ji,jj)  - smrft(ji,jj)
1117
1118            !-----------------------------------------------------------------
1119            ! 3.5) Compute properties of new ridges
1120            !-----------------------------------------------------------------
1121            !-------------
1122            ! Salinity
1123            !-------------
1124            smsw(ji,jj)  = vsw(ji,jj) * sss_m(ji,jj)                      ! salt content of seawater frozen in voids !! MV HC2014
1125            zsrdg2       = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge
1126
1127            srdg2(ji,jj) = MIN( s_i_max * vrdg2(ji,jj) , zsrdg2 )         ! impose a maximum salinity
1128           
1129            !                                                             ! excess of salt is flushed into the ocean
1130            !sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice
1131            ! MV HC 2014 this previous line seems ok too... Maybe units are not that good
1132
1133            !rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic    ! gurvan: increase in ice volume du to seawater frozen in voids             
1134            ! MV HC 2014 this previous line seems ok, i'm not sure at this moment of the sign convention
1135
1136            !------------------------------------           
1137            ! 3.6 Increment ridging diagnostics
1138            !------------------------------------           
1139
1140            !        jl1 looping 1-jpl
1141            !           ij looping 1-icells
1142
1143            dardg1dt   (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj)
1144            dardg2dt   (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj)
1145            opening    (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice
1146
1147            IF( con_i )   vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj)
1148
1149            !------------------------------------------           
1150            ! 3.7 Put the snow somewhere in the ocean
1151            !------------------------------------------           
1152            !  Place part of the snow lost by ridging into the ocean.
1153            !  Note that esnow_mlt < 0; the ocean must cool to melt snow.
1154            !  If the ocean temp = Tf already, new ice must grow.
1155            !  During the next time step, thermo_rates will determine whether
1156            !  the ocean cools or new ice grows.
1157            !        jl1 looping 1-jpl
1158            !           ij looping 1-icells
1159
1160            msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosn*vsrdg(ji,jj)*(1.0-fsnowrdg)   &   ! rafting included
1161               &                                + rhosn*vsrft(ji,jj)*(1.0-fsnowrft)
1162
1163            ! MV HC 2014 - I think the units are wrong here... shit there is a 10-9 in the snow...
1164            ! we should therefore do something
1165            esnow_mlt(ji,jj) = esnow_mlt(ji,jj) + esrdg(ji,jj)*(1.0-fsnowrdg)         &   !rafting included
1166               &                                + esrft(ji,jj)*(1.0-fsnowrft)         
1167
1168            !-----------------------------------------------------------------
1169            ! 3.8 Compute quantities used to apportion ice among categories
1170            ! in the n2 loop below
1171            !-----------------------------------------------------------------
1172
1173            !        jl1 looping 1-jpl
1174            !           ij looping 1-icells
1175
1176            dhr (ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1)
1177            dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1)
1178
1179         END DO                 ! ij
1180
1181         !--------------------------------------------------------------------
1182         ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
1183         !      compute ridged ice enthalpy
1184         !--------------------------------------------------------------------
1185         DO jk = 1, nlay_i
1186!CDIR NODEP
1187            DO ij = 1, icells
1188               ji = indxi(ij)
1189               jj = indxj(ij)
1190               ! heat content of ridged ice
1191               erdg1(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por ) 
1192               eirft(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj)
1193               e_i  (ji,jj,jk,jl1)  = e_i(ji,jj,jk,jl1) - erdg1(ji,jj,jk) - eirft(ji,jj,jk)
1194               
1195               
1196               ! enthalpy of the trapped seawater (J/m2, >0)
1197               zsstK  = sst_m(ji,jj) + rt0
1198               ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 )
1199
1200               ! heat flux to the ocean
1201               rdq_ice(ji,jj) = rdq_ice(ji,jj) + ersw(ji,jj,jk) !! MV HC 2014
1202
1203               ! Correct dimensions to avoid big values
1204               ersw(ji,jj,jk)   = ersw(ji,jj,jk) * 1.e-09
1205
1206               ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
1207               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean
1208               !! MV HC 2014
1209               ersw (ji,jj,jk)  = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / REAL( nlay_i )
1210
1211               erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk)
1212
1213            END DO ! ij
1214         END DO !jk
1215
1216
1217         IF( con_i ) THEN
1218            DO jk = 1, nlay_i
1219!CDIR NODEP
1220               DO ij = 1, icells
1221                  ji = indxi(ij)
1222                  jj = indxj(ij)
1223                  eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - erdg1(ji,jj,jk)
1224               END DO ! ij
1225            END DO !jk
1226         ENDIF
1227
1228         IF( large_afrac ) THEN   ! there is a bug
1229!CDIR NODEP
1230            DO ij = 1, icells
1231               ji = indxi(ij)
1232               jj = indxj(ij)
1233               IF( afrac(ji,jj) > kamax + epsi10 ) THEN
1234                  WRITE(numout,*) ''
1235                  WRITE(numout,*) ' ardg > a_i'
1236                  WRITE(numout,*) ' ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1)
1237               ENDIF
1238            END DO
1239         ENDIF
1240         IF( large_afrft ) THEN  ! there is a bug
1241!CDIR NODEP
1242            DO ij = 1, icells
1243               ji = indxi(ij)
1244               jj = indxj(ij)
1245               IF( afrft(ji,jj) > kamax + epsi10 ) THEN
1246                  WRITE(numout,*) ''
1247                  WRITE(numout,*) ' arft > a_i'
1248                  WRITE(numout,*) ' arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1)
1249               ENDIF
1250            END DO
1251         ENDIF
1252
1253         !-------------------------------------------------------------------------------
1254         ! 4) Add area, volume, and energy of new ridge to each category jl2
1255         !-------------------------------------------------------------------------------
1256         !        jl1 looping 1-jpl
1257         DO jl2  = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 
1258            ! over categories to which ridged ice is transferred
1259!CDIR NODEP
1260            DO ij = 1, icells
1261               ji = indxi(ij)
1262               jj = indxj(ij)
1263
1264               ! Compute the fraction of ridged ice area and volume going to
1265               ! thickness category jl2.
1266               ! Transfer area, volume, and energy accordingly.
1267
1268               IF( hrmin(ji,jj,jl1) >= hi_max(jl2)  .OR.        &
1269                   hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN
1270                  hL = 0._wp
1271                  hR = 0._wp
1272               ELSE
1273                  hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
1274                  hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2)   )
1275               ENDIF
1276
1277               ! fraction of ridged ice area and volume going to n2
1278               farea = ( hR - hL ) / dhr(ji,jj) 
1279               fvol(ji,jj) = ( hR*hR - hL*hL ) / dhr2(ji,jj)
1280
1281               a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ardg2 (ji,jj) * farea
1282               v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + vrdg2 (ji,jj) * fvol(ji,jj)
1283               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * fsnowrdg
1284               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * fsnowrdg
1285               smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + srdg2 (ji,jj) * fvol(ji,jj)
1286               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirdg2(ji,jj) * farea
1287
1288            END DO ! ij
1289
1290            ! Transfer ice energy to category jl2 by ridging
1291            DO jk = 1, nlay_i
1292!CDIR NODEP
1293               DO ij = 1, icells
1294                  ji = indxi(ij)
1295                  jj = indxj(ij)
1296                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj)*erdg2(ji,jj,jk)
1297               END DO
1298            END DO
1299            !
1300         END DO                 ! jl2 (new ridges)           
1301
1302         DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 
1303
1304!CDIR NODEP
1305            DO ij = 1, icells
1306               ji = indxi(ij)
1307               jj = indxj(ij)
1308               ! Compute the fraction of rafted ice area and volume going to
1309               ! thickness category jl2, transfer area, volume, and energy accordingly.
1310               !
1311               IF( hraft(ji,jj,jl1) <= hi_max(jl2)  .AND.        &
1312                   hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN
1313                  a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + arft2 (ji,jj)
1314                  v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + virft (ji,jj)
1315                  v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrft (ji,jj) * fsnowrft
1316                  e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrft (ji,jj) * fsnowrft
1317                  smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + smrft (ji,jj)   
1318                  oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirft2(ji,jj)   
1319               ENDIF ! hraft
1320               !
1321            END DO ! ij
1322
1323            ! Transfer rafted ice energy to category jl2
1324            DO jk = 1, nlay_i
1325!CDIR NODEP
1326               DO ij = 1, icells
1327                  ji = indxi(ij)
1328                  jj = indxj(ij)
1329                  IF(  hraft(ji,jj,jl1)  <=  hi_max(jl2)   .AND.        &
1330                       hraft(ji,jj,jl1)  >   hi_max(jl2-1)  ) THEN
1331                     e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk)
1332                  ENDIF
1333               END DO           ! ij
1334            END DO !jk
1335
1336         END DO ! jl2
1337
1338      END DO ! jl1 (deforming categories)
1339
1340      ! Conservation check
1341      IF ( con_i ) THEN
1342         CALL lim_column_sum (jpl,   v_i, vice_final)
1343         fieldid = ' v_i : limitd_me '
1344         CALL lim_cons_check (vice_init, vice_final, 1.0e-6, fieldid) 
1345
1346         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_final )
1347         fieldid = ' e_i : limitd_me '
1348         CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid) 
1349
1350         DO ji = mi0(jiindx), mi1(jiindx)
1351            DO jj = mj0(jjindx), mj1(jjindx)
1352               WRITE(numout,*) ' vice_init  : ', vice_init (ji,jj)
1353               WRITE(numout,*) ' vice_final : ', vice_final(ji,jj)
1354               WRITE(numout,*) ' eice_init  : ', eice_init (ji,jj)
1355               WRITE(numout,*) ' eice_final : ', eice_final(ji,jj)
1356            END DO
1357         END DO
1358      ENDIF
1359      !
1360      CALL wrk_dealloc( (jpi+1)*(jpj+1),      indxi, indxj )
1361      CALL wrk_dealloc( jpi, jpj,             vice_init, vice_final, eice_init, eice_final )
1362      CALL wrk_dealloc( jpi, jpj,             afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )
1363      CALL wrk_dealloc( jpi, jpj,             vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw )
1364      CALL wrk_dealloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
1365      CALL wrk_dealloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnon_init, esnon_init, smv_i_init, oa_i_init )
1366      CALL wrk_dealloc( jpi, jpj, jkmax,      eirft, erdg1, erdg2, ersw )
1367      CALL wrk_dealloc( jpi, jpj, jkmax, jpl, eicen_init )
1368      !
1369   END SUBROUTINE lim_itd_me_ridgeshift
1370
1371
1372   SUBROUTINE lim_itd_me_asumr
1373      !!-----------------------------------------------------------------------------
1374      !!                ***  ROUTINE lim_itd_me_asumr ***
1375      !!
1376      !! ** Purpose :   finds total fractional area
1377      !!
1378      !! ** Method  :   Find the total area of ice plus open water in each grid cell.
1379      !!              This is similar to the aggregate_area subroutine except that the
1380      !!              total area can be greater than 1, so the open water area is
1381      !!              included in the sum instead of being computed as a residual.
1382      !!-----------------------------------------------------------------------------
1383      INTEGER ::   jl   ! dummy loop index
1384      !!-----------------------------------------------------------------------------
1385      !
1386      asum(:,:) = ato_i(:,:)                    ! open water
1387      DO jl = 1, jpl                            ! ice categories
1388         asum(:,:) = asum(:,:) + a_i(:,:,jl)
1389      END DO
1390      !
1391   END SUBROUTINE lim_itd_me_asumr
1392
1393
1394   SUBROUTINE lim_itd_me_init
1395      !!-------------------------------------------------------------------
1396      !!                   ***  ROUTINE lim_itd_me_init ***
1397      !!
1398      !! ** Purpose :   Physical constants and parameters linked
1399      !!                to the mechanical ice redistribution
1400      !!
1401      !! ** Method  :   Read the namiceitdme namelist
1402      !!                and check the parameters values
1403      !!                called at the first timestep (nit000)
1404      !!
1405      !! ** input   :   Namelist namiceitdme
1406      !!-------------------------------------------------------------------
1407      NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft,& 
1408         Gstar, astar,                                &
1409         Hstar, raftswi, hparmeter, Craft, ridge_por, &
1410         sal_max_ridge,  partfun_swi, transfun_swi,   &
1411         brinstren_swi
1412      !!-------------------------------------------------------------------
1413      !
1414      REWIND( numnam_ice )                   ! read namiceitdme namelist
1415      READ  ( numnam_ice , namiceitdme)
1416      !
1417      IF (lwp) THEN                          ! control print
1418         WRITE(numout,*)
1419         WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution '
1420         WRITE(numout,*)' ~~~~~~~~~~~~~~~'
1421         WRITE(numout,*)'   Switch choosing the ice redistribution scheme           ridge_scheme_swi', ridge_scheme_swi 
1422         WRITE(numout,*)'   Fraction of shear energy contributing to ridging        Cs              ', Cs 
1423         WRITE(numout,*)'   Ratio of ridging work to PotEner change in ridging      Cf              ', Cf 
1424         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        fsnowrdg        ', fsnowrdg 
1425         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        fsnowrft        ', fsnowrft 
1426         WRITE(numout,*)'   Fraction of total ice coverage contributing to ridging  Gstar           ', Gstar
1427         WRITE(numout,*)'   Equivalent to G* for an exponential part function       astar           ', astar
1428         WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     Hstar           ', Hstar
1429         WRITE(numout,*)'   Rafting of ice sheets or not                            raftswi         ', raftswi
1430         WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       hparmeter       ', hparmeter
1431         WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  Craft           ', Craft 
1432         WRITE(numout,*)'   Initial porosity of ridges                              ridge_por       ', ridge_por
1433         WRITE(numout,*)'   Maximum salinity of ridging ice                         sal_max_ridge   ', sal_max_ridge
1434         WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    partfun_swi     ', partfun_swi
1435         WRITE(numout,*)'   Switch for tran. function (0) linear (1) exponential    transfun_swi    ', transfun_swi
1436         WRITE(numout,*)'   Switch for including brine volume in ice strength comp. brinstren_swi   ', brinstren_swi
1437      ENDIF
1438      !
1439   END SUBROUTINE lim_itd_me_init
1440
1441
1442   SUBROUTINE lim_itd_me_zapsmall
1443      !!-------------------------------------------------------------------
1444      !!                   ***  ROUTINE lim_itd_me_zapsmall ***
1445      !!
1446      !! ** Purpose :   Remove too small sea ice areas and correct salt fluxes
1447      !!
1448      !! history :
1449      !! author: William H. Lipscomb, LANL
1450      !! Nov 2003:  Modified by Julie Schramm to conserve volume and energy
1451      !! Sept 2004: Modified by William Lipscomb; replaced normalize_state with
1452      !!            additions to local freshwater, salt, and heat fluxes
1453      !!  9.0, LIM3.0 - 02-2006 (M. Vancoppenolle) original code
1454      !!-------------------------------------------------------------------
1455      INTEGER ::   ji, jj, jl, jk   ! dummy loop indices
1456      INTEGER ::   icells           ! number of cells with ice to zap
1457
1458      REAL(wp), POINTER, DIMENSION(:,:) ::   zmask   ! 2D workspace
1459      REAL(wp)                          ::   zmask_glo
1460!!gm      REAL(wp) ::   xtmp      ! temporary variable
1461      !!-------------------------------------------------------------------
1462
1463      CALL wrk_alloc( jpi, jpj, zmask )
1464
1465      DO jl = 1, jpl
1466
1467         !-----------------------------------------------------------------
1468         ! Count categories to be zapped.
1469         ! Abort model in case of negative area.
1470         !-----------------------------------------------------------------
1471         icells = 0
1472         zmask(:,:)  = 0._wp
1473         DO jj = 1, jpj
1474            DO ji = 1, jpi
1475               IF(  ( a_i(ji,jj,jl) >= -epsi10 .AND. a_i(ji,jj,jl) <  0._wp   ) .OR.   &
1476                  & ( a_i(ji,jj,jl) >  0._wp   .AND. a_i(ji,jj,jl) <= epsi10  ) .OR.   &
1477                  & ( v_i(ji,jj,jl) == 0._wp   .AND. a_i(ji,jj,jl) >  0._wp   ) .OR.   &
1478                  & ( v_i(ji,jj,jl) >  0._wp   .AND. v_i(ji,jj,jl) <= epsi10  ) )   zmask(ji,jj) = 1._wp
1479            END DO
1480         END DO
1481         !zmask_glo = glob_sum(zmask)
1482         !IF( ln_nicep .AND. lwp ) WRITE(numout,*) zmask_glo, ' cells of ice zapped in the ocean '
1483
1484         !-----------------------------------------------------------------
1485         ! Zap ice energy and use ocean heat to melt ice
1486         !-----------------------------------------------------------------
1487
1488         DO jk = 1, nlay_i
1489            DO jj = 1 , jpj
1490               DO ji = 1 , jpi
1491!!gm                  xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) * r1_rdtice
1492!!gm                  xtmp = xtmp * unit_fac
1493                  ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp
1494                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) )
1495               END DO
1496            END DO
1497         END DO
1498
1499         DO jj = 1 , jpj
1500            DO ji = 1 , jpi
1501
1502               !-----------------------------------------------------------------
1503               ! Zap snow energy and use ocean heat to melt snow
1504               !-----------------------------------------------------------------
1505               !           xtmp = esnon(i,j,n) / dt ! < 0
1506               !           fhnet(i,j)      = fhnet(i,j)      + xtmp
1507               !           fhnet_hist(i,j) = fhnet_hist(i,j) + xtmp
1508               ! xtmp is greater than 0
1509               ! fluxes are positive to the ocean
1510               ! here the flux has to be negative for the ocean
1511!!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice
1512               !           fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp
1513
1514!!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice !RB   ???????
1515
1516               t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) )
1517
1518               !-----------------------------------------------------------------
1519               ! zap ice and snow volume, add water and salt to ocean
1520               !-----------------------------------------------------------------
1521
1522               !           xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) / dt
1523               !           sfx_res(ji,jj)  = sfx_res(ji,jj) + ( sss_m(ji,jj)                  )   &
1524               !                                            * rhosn * v_s(ji,jj,jl) * r1_rdtice
1525               !           sfx_res(ji,jj)  = sfx_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) )   &
1526               !                                            * rhoic * v_i(ji,jj,jl) * r1_rdtice
1527               !           sfx (i,j)      = sfx (i,j)      + xtmp
1528
1529               ato_i(ji,jj)    = a_i  (ji,jj,jl) *       zmask(ji,jj)   + ato_i(ji,jj)
1530               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
1531               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
1532               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
1533               t_su (ji,jj,jl) = t_su (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj)
1534               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
1535               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
1536               !
1537            END DO
1538         END DO
1539         !
1540      END DO                 ! jl
1541      !
1542      CALL wrk_dealloc( jpi, jpj, zmask )
1543      !
1544   END SUBROUTINE lim_itd_me_zapsmall
1545
1546#else
1547   !!----------------------------------------------------------------------
1548   !!   Default option         Empty module          NO LIM-3 sea-ice model
1549   !!----------------------------------------------------------------------
1550CONTAINS
1551   SUBROUTINE lim_itd_me           ! Empty routines
1552   END SUBROUTINE lim_itd_me
1553   SUBROUTINE lim_itd_me_icestrength
1554   END SUBROUTINE lim_itd_me_icestrength
1555   SUBROUTINE lim_itd_me_sort
1556   END SUBROUTINE lim_itd_me_sort
1557   SUBROUTINE lim_itd_me_init
1558   END SUBROUTINE lim_itd_me_init
1559   SUBROUTINE lim_itd_me_zapsmall
1560   END SUBROUTINE lim_itd_me_zapsmall
1561#endif
1562   !!======================================================================
1563END MODULE limitd_me
Note: See TracBrowser for help on using the repository browser.