New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limitd_me.F90 in branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 3402

Last change on this file since 3402 was 3402, checked in by acc, 12 years ago

Branch: dev_r3385_NOCS04_HAMF; #665. Stage 2 of 2012 development: suppression of emps array and introduction of sfx (salt flux) array with associated code to setup the options for embedding the seaice into the ocean

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