source: trunk/NEMO/LIM_SRC_3/limitd_me.F90 @ 903

Last change on this file since 903 was 903, checked in by ctlod, 13 years ago

correct the wrong name of a dummy subroutine, see ticket: #120

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