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 trunk/NEMO/LIM_SRC_3 – NEMO

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

Last change on this file since 863 was 863, checked in by ctlod, 16 years ago

Add control prints

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