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

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 4333

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

remove remaining bugs in LIM3, so that it can run in both regional and global config

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