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

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 4788

Last change on this file since 4788 was 4788, checked in by clem, 10 years ago

bug in LIM3. Porosity during ridging was not taken into account properly. See ticket coming up from Barthelemy

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