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 @ 3517

Last change on this file since 3517 was 3517, checked in by gm, 10 years ago

gm: Branch: dev_r3385_NOCS04_HAMF; #665. update sbccpl ; change LIM3 from equivalent salt flux to salt flux and mass flux

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