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

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

source: branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 5965

Last change on this file since 5965 was 5965, checked in by timgraham, 8 years ago

Upgraded branch to r5518 of trunk (v3.6 stable revision)

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