source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 4161

Last change on this file since 4161 was 4161, checked in by cetlod, 8 years ago

dev_LOCEAN_2013 : merge in the 3rd dev branch dev_r4028_CNRS_LIM3, see ticket #1169

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