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

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

source: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 3625

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

Branch dev_NOC_2012_r3555. #1006. Step 7. Check in code now merged with dev_r3385_NOCS04_HAMF

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