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

source: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 4506

Last change on this file since 4506 was 4506, checked in by vancop, 10 years ago

[Heat conservation in LIM3, part HC1 (LIM_SRC_3_HC17)]

  • Property svn:keywords set to Id
File size: 72.8 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
[4045]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   !!----------------------------------------------------------------------
[4045]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
[4072]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
[4332]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   !-----------------------------------------------------------------------
[4332]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
[4332]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
[4045]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      !!-----------------------------------------------------------------------------
[4072]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
[4045]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
[4045]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
[4045]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.
[4332]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
[4332]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
[4332]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
[4045]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
[4332]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
[4506]365            ! there is a unit problem for e_snow_mlt that we should fix
[3625]366            fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice     ! heat sink for ocean
[825]367
368         END DO
369      END DO
370
371      ! Check if there is a ridging error
372      DO jj = 1, jpj
373         DO ji = 1, jpi
[4332]374            IF( ABS( asum(ji,jj) - kamax)  >  epsi10 ) THEN   ! there is a bug
[825]375               WRITE(numout,*) ' '
376               WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj)
377               WRITE(numout,*) ' limitd_me '
378               WRITE(numout,*) ' POINT : ', ji, jj
379               WRITE(numout,*) ' jpl, a_i, athorn '
380               WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0)
381               DO jl = 1, jpl
382                  WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl)
383               END DO
384            ENDIF  ! asum
385
386         END DO !ji
387      END DO !jj
388
389      ! Conservation check
[834]390      IF ( con_i ) THEN
391         CALL lim_column_sum (jpl,   v_i, vt_i_final)
392         fieldid = ' v_i : limitd_me '
393         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
394      ENDIF
[825]395
[921]396      !-----------------------------------------------------------------------------!
[4332]397      ! 6) Updating state variables and trend terms (done in limupdate)
[921]398      !-----------------------------------------------------------------------------!
[834]399      CALL lim_var_glo2eqv
[825]400      CALL lim_itd_me_zapsmall
401
[4045]402      !--------------------------------
403      ! Update mass/salt fluxes (clem)
404      !--------------------------------
405      DO jl = 1, jpl
406         DO jj = 1, jpj 
407            DO ji = 1, jpi
408               diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice
409               rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic 
410               rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn 
411               sfx_mec(ji,jj) = sfx_mec(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic * r1_rdtice 
412            END DO
413         END DO
414      END DO
415
[863]416      IF(ln_ctl) THEN     ! Control print
[867]417         CALL prt_ctl_info(' ')
418         CALL prt_ctl_info(' - Cell values : ')
419         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
[863]420         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_me  : cell area :')
421         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me  : at_i      :')
422         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me  : vt_i      :')
423         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_me  : vt_s      :')
424         DO jl = 1, jpl
[867]425            CALL prt_ctl_info(' ')
[863]426            CALL prt_ctl_info(' - Category : ', ivar1=jl)
427            CALL prt_ctl_info('   ~~~~~~~~~~')
428            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_itd_me  : a_i      : ')
429            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_itd_me  : ht_i     : ')
430            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_itd_me  : ht_s     : ')
431            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_itd_me  : v_i      : ')
432            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_itd_me  : v_s      : ')
433            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_itd_me  : e_s      : ')
434            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_itd_me  : t_su     : ')
435            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_itd_me  : t_snow   : ')
436            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_itd_me  : sm_i     : ')
437            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_itd_me  : smv_i    : ')
438            DO jk = 1, nlay_i
[867]439               CALL prt_ctl_info(' ')
[863]440               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
441               CALL prt_ctl_info('   ~~~~~~~')
442               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_me  : t_i      : ')
443               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_me  : e_i      : ')
444            END DO
445         END DO
446      ENDIF
[825]447
[4045]448      ! -------------------------------
449      !- check conservation (C Rousset)
450      IF (ln_limdiahsb) THEN
451         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
452         zchk_fw  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b
453 
454         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice
455         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic )
456
457         zchk_vmin = glob_min(v_i)
458         zchk_amax = glob_max(SUM(a_i,dim=3))
459         zchk_amin = glob_min(a_i)
460       
461         IF(lwp) THEN
462            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limitd_me) = ',(zchk_v_i * rday)
463            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * rday)
464            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_me) = ',(zchk_vmin * 1.e-3)
465            IF ( zchk_amax >  kamax+epsi10  ) WRITE(numout,*) 'violation a_i>amax            (limitd_me) = ',zchk_amax
466            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limitd_me) = ',zchk_amin
467         ENDIF
468      ENDIF
469      !- check conservation (C Rousset)
470      ! -------------------------------
471
472      ENDIF  ! ln_limdyn=.true.
[3625]473      !
[3294]474      CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
[2715]475      !
[4045]476      CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zsmvold )   ! clem
477      !
[4072]478      IF( nn_timing == 1 )  CALL timing_stop('limitd_me')
[825]479   END SUBROUTINE lim_itd_me
480
481
[2715]482   SUBROUTINE lim_itd_me_icestrength( kstrngth )
[921]483      !!----------------------------------------------------------------------
484      !!                ***  ROUTINE lim_itd_me_icestrength ***
485      !!
[2715]486      !! ** Purpose :   computes ice strength used in dynamics routines of ice thickness
[921]487      !!
[2715]488      !! ** Method  :   Compute the strength of the ice pack, defined as the energy (J m-2)
489      !!              dissipated per unit area removed from the ice pack under compression,
490      !!              and assumed proportional to the change in potential energy caused
491      !!              by ridging. Note that only Hibler's formulation is stable and that
492      !!              ice strength has to be smoothed
493      !!
[921]494      !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
495      !!----------------------------------------------------------------------
[2715]496      INTEGER, INTENT(in) ::   kstrngth    ! = 1 for Rothrock formulation, 0 for Hibler (1979)
[921]497
[2715]498      INTEGER ::   ji,jj, jl   ! dummy loop indices
499      INTEGER ::   ksmooth     ! smoothing the resistance to deformation
500      INTEGER ::   numts_rm    ! number of time steps for the P smoothing
501      REAL(wp) ::   hi, zw1, zp, zdummy, zzc, z1_3   ! local scalars
[3294]502      REAL(wp), POINTER, DIMENSION(:,:) ::   zworka   ! temporary array used here
[2715]503      !!----------------------------------------------------------------------
[825]504
[3294]505      CALL wrk_alloc( jpi, jpj, zworka )
[825]506
[921]507      !------------------------------------------------------------------------------!
508      ! 1) Initialize
509      !------------------------------------------------------------------------------!
[2715]510      strength(:,:) = 0._wp
[825]511
[921]512      !------------------------------------------------------------------------------!
513      ! 2) Compute thickness distribution of ridging and ridged ice
514      !------------------------------------------------------------------------------!
[825]515      CALL lim_itd_me_ridgeprep
516
[921]517      !------------------------------------------------------------------------------!
518      ! 3) Rothrock(1975)'s method
519      !------------------------------------------------------------------------------!
[2715]520      IF( kstrngth == 1 ) THEN
521         z1_3 = 1._wp / 3._wp
[825]522         DO jl = 1, jpl
523            DO jj= 1, jpj
524               DO ji = 1, jpi
[2715]525                  !
[4332]526                  IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN
[825]527                     hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
528                     !----------------------------
529                     ! PE loss from deforming ice
530                     !----------------------------
[2715]531                     strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * hi * hi
[825]532
533                     !--------------------------
534                     ! PE gain from rafting ice
535                     !--------------------------
[2715]536                     strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * hi * hi
[825]537
538                     !----------------------------
539                     ! PE gain from ridging ice
540                     !----------------------------
[2715]541                     strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl)/krdg(ji,jj,jl)     &
542                        * z1_3 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) / ( hrmax(ji,jj,jl)-hrmin(ji,jj,jl) )   
543!!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...                   
[4332]544                  ENDIF            ! aicen > epsi10
[2715]545                  !
[825]546               END DO ! ji
547            END DO !jj
548         END DO !jl
549
[2715]550         zzc = Cf * Cp     ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and Cf accounts for frictional dissipation
551         strength(:,:) = zzc * strength(:,:) / aksum(:,:)
[921]552
[825]553         ksmooth = 1
554
[921]555         !------------------------------------------------------------------------------!
556         ! 4) Hibler (1979)' method
557         !------------------------------------------------------------------------------!
[825]558      ELSE                      ! kstrngth ne 1:  Hibler (1979) form
[2715]559         !
560         strength(:,:) = Pstar * vt_i(:,:) * EXP( - C_rhg * ( 1._wp - at_i(:,:) )  )
561         !
[825]562         ksmooth = 1
[2715]563         !
[825]564      ENDIF                     ! kstrngth
565
[921]566      !
567      !------------------------------------------------------------------------------!
568      ! 5) Impact of brine volume
569      !------------------------------------------------------------------------------!
570      ! CAN BE REMOVED
571      !
[2715]572      IF ( brinstren_swi == 1 ) THEN
[825]573
574         DO jj = 1, jpj
575            DO ji = 1, jpi
576               IF ( bv_i(ji,jj) .GT. 0.0 ) THEN
577                  zdummy = MIN ( bv_i(ji,jj), 0.10 ) * MIN( bv_i(ji,jj), 0.10 )
578               ELSE
579                  zdummy = 0.0
580               ENDIF
581               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0)))
582            END DO              ! j
583         END DO                 ! i
584
585      ENDIF
586
[921]587      !
588      !------------------------------------------------------------------------------!
589      ! 6) Smoothing ice strength
590      !------------------------------------------------------------------------------!
591      !
[825]592      !-------------------
593      ! Spatial smoothing
594      !-------------------
[2715]595      IF ( ksmooth == 1 ) THEN
[825]596
597         CALL lbc_lnk( strength, 'T', 1. )
598
[869]599         DO jj = 2, jpj - 1
600            DO ji = 2, jpi - 1
[4332]601               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN ! ice is
[921]602                  ! present
[825]603                  zworka(ji,jj) = 4.0 * strength(ji,jj)              &
[3625]604                     &          + strength(ji-1,jj) * tms(ji-1,jj) & 
605                     &          + strength(ji+1,jj) * tms(ji+1,jj) & 
606                     &          + strength(ji,jj-1) * tms(ji,jj-1) & 
607                     &          + strength(ji,jj+1) * tms(ji,jj+1)   
[825]608
[2715]609                  zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj) + tms(ji,jj-1) + tms(ji,jj+1)
[825]610                  zworka(ji,jj) = zworka(ji,jj) / zw1
611               ELSE
[3625]612                  zworka(ji,jj) = 0._wp
[825]613               ENDIF
614            END DO
615         END DO
616
[869]617         DO jj = 2, jpj - 1
618            DO ji = 2, jpi - 1
[825]619               strength(ji,jj) = zworka(ji,jj)
620            END DO
621         END DO
[869]622         CALL lbc_lnk( strength, 'T', 1. )
[825]623
624      ENDIF ! ksmooth
625
626      !--------------------
627      ! Temporal smoothing
628      !--------------------
[2715]629      IF ( numit == nit000 + nn_fsbc - 1 ) THEN
[825]630         strp1(:,:) = 0.0           
631         strp2(:,:) = 0.0           
632      ENDIF
633
[2715]634      IF ( ksmooth == 2 ) THEN
[921]635
636
[825]637         CALL lbc_lnk( strength, 'T', 1. )
[921]638
[825]639         DO jj = 1, jpj - 1
640            DO ji = 1, jpi - 1
[4332]641               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN       ! ice is present
[825]642                  numts_rm = 1 ! number of time steps for the running mean
643                  IF ( strp1(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1
644                  IF ( strp2(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1
[2715]645                  zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm
[825]646                  strp2(ji,jj) = strp1(ji,jj)
647                  strp1(ji,jj) = strength(ji,jj)
648                  strength(ji,jj) = zp
649
650               ENDIF
651            END DO
652         END DO
653
654      ENDIF ! ksmooth
[921]655
[2715]656      CALL lbc_lnk( strength, 'T', 1. )      ! Boundary conditions
[825]657
[3294]658      CALL wrk_dealloc( jpi, jpj, zworka )
[2715]659      !
[921]660   END SUBROUTINE lim_itd_me_icestrength
[825]661
662
[2715]663   SUBROUTINE lim_itd_me_ridgeprep
[921]664      !!---------------------------------------------------------------------!
665      !!                ***  ROUTINE lim_itd_me_ridgeprep ***
666      !!
[2715]667      !! ** Purpose :   preparation for ridging and strength calculations
[921]668      !!
[2715]669      !! ** Method  :   Compute the thickness distribution of the ice and open water
670      !!              participating in ridging and of the resulting ridges.
[921]671      !!---------------------------------------------------------------------!
[2715]672      INTEGER ::   ji,jj, jl    ! dummy loop indices
673      INTEGER ::   krdg_index   !
674      REAL(wp) ::   Gstari, astari, hi, hrmean, zdummy   ! local scalar
[3294]675      REAL(wp), POINTER, DIMENSION(:,:)   ::   zworka    ! temporary array used here
676      REAL(wp), POINTER, DIMENSION(:,:,:) ::   Gsum      ! Gsum(n) = sum of areas in categories 0 to n
677      !------------------------------------------------------------------------------!
[825]678
[3294]679      CALL wrk_alloc( jpi,jpj, zworka )
680      CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
[825]681
682      Gstari     = 1.0/Gstar   
683      astari     = 1.0/astar   
684      aksum(:,:)    = 0.0
685      athorn(:,:,:) = 0.0
686      aridge(:,:,:) = 0.0
687      araft (:,:,:) = 0.0
688      hrmin(:,:,:)  = 0.0 
689      hrmax(:,:,:)  = 0.0 
690      hraft(:,:,:)  = 0.0 
691      krdg (:,:,:)  = 1.0
692
[921]693      !     ! Zero out categories with very small areas
[825]694      CALL lim_itd_me_zapsmall
695
[921]696      !------------------------------------------------------------------------------!
697      ! 1) Participation function
698      !------------------------------------------------------------------------------!
[825]699
700      ! Compute total area of ice plus open water.
701      CALL lim_itd_me_asumr
702      ! This is in general not equal to one
703      ! because of divergence during transport
704
705      ! Compute cumulative thickness distribution function
706      ! Compute the cumulative thickness distribution function Gsum,
707      ! where Gsum(n) is the fractional area in categories 0 to n.
708      ! initial value (in h = 0) equals open water area
709
[2715]710      Gsum(:,:,-1) = 0._wp
[825]711
712      DO jj = 1, jpj
713         DO ji = 1, jpi
[4332]714            IF( ato_i(ji,jj) > epsi10 ) THEN   ;   Gsum(ji,jj,0) = ato_i(ji,jj)
[2715]715            ELSE                               ;   Gsum(ji,jj,0) = 0._wp
[825]716            ENDIF
717         END DO
718      END DO
719
720      ! for each value of h, you have to add ice concentration then
721      DO jl = 1, jpl
722         DO jj = 1, jpj 
723            DO ji = 1, jpi
[4332]724               IF( a_i(ji,jj,jl) .GT. epsi10 ) THEN   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl)
[2715]725               ELSE                                   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1)
[825]726               ENDIF
727            END DO
728         END DO
729      END DO
730
731      ! Normalize the cumulative distribution to 1
[2715]732      zworka(:,:) = 1._wp / Gsum(:,:,jpl)
[825]733      DO jl = 0, jpl
[2715]734         Gsum(:,:,jl) = Gsum(:,:,jl) * zworka(:,:)
[825]735      END DO
736
[921]737      ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
738      !--------------------------------------------------------------------------------------------------
739      ! Compute the participation function athorn; this is analogous to
740      ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
741      ! area lost from category n due to ridging/closing
742      ! athorn(n)   = total area lost due to ridging/closing
743      ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
744      !
745      ! The expressions for athorn are found by integrating b(h)g(h) between
746      ! the category boundaries.
747      !-----------------------------------------------------------------
[825]748
749      krdg_index = 1
750
[2715]751      IF( krdg_index == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975)
752         DO jl = 0, ice_cat_bounds(1,2)       ! only undeformed ice participates
[921]753            DO jj = 1, jpj 
754               DO ji = 1, jpi
[2715]755                  IF( Gsum(ji,jj,jl) < Gstar) THEN
[921]756                     athorn(ji,jj,jl) = Gstari * (Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * &
757                        (2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari)
758                  ELSEIF (Gsum(ji,jj,jl-1) < Gstar) THEN
759                     athorn(ji,jj,jl) = Gstari * (Gstar-Gsum(ji,jj,jl-1)) *  &
760                        (2.0 - (Gsum(ji,jj,jl-1)+Gstar)*Gstari)
761                  ELSE
762                     athorn(ji,jj,jl) = 0.0
763                  ENDIF
764               END DO ! ji
765            END DO ! jj
766         END DO ! jl
[825]767
[2715]768      ELSE                             !--- Exponential, more stable formulation (Lipscomb et al, 2007)
769         !                       
770         zdummy = 1._wp / ( 1._wp - EXP(-astari) )        ! precompute exponential terms using Gsum as a work array
[825]771
[921]772         DO jl = -1, jpl
[2715]773            Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
[921]774         END DO !jl
775         DO jl = 0, ice_cat_bounds(1,2)
[2715]776             athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
777         END DO
778         !
[825]779      ENDIF ! krdg_index
780
[2715]781      IF( raftswi == 1 ) THEN      ! Ridging and rafting ice participation functions
782         !
[825]783         DO jl = 1, jpl
784            DO jj = 1, jpj 
785               DO ji = 1, jpi
[2715]786                  IF ( athorn(ji,jj,jl) .GT. 0._wp ) THEN
787!!gm  TANH( -X ) = - TANH( X )  so can be computed only 1 time....
788                     aridge(ji,jj,jl) = ( TANH (  Craft * ( ht_i(ji,jj,jl) - hparmeter ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
789                     araft (ji,jj,jl) = ( TANH ( -Craft * ( ht_i(ji,jj,jl) - hparmeter ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
790                     IF ( araft(ji,jj,jl) < epsi06 )   araft(ji,jj,jl)  = 0._wp
791                     aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 )
[825]792                  ENDIF ! athorn
793               END DO ! ji
794            END DO ! jj
795         END DO ! jl
796
797      ELSE  ! raftswi = 0
[2715]798         !
[825]799         DO jl = 1, jpl
[2715]800            aridge(:,:,jl) = athorn(:,:,jl)
[825]801         END DO
[2715]802         !
[825]803      ENDIF
804
[2715]805      IF ( raftswi == 1 ) THEN
[825]806
[4332]807         IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi10 ) THEN
[921]808            DO jl = 1, jpl
809               DO jj = 1, jpj
810                  DO ji = 1, jpi
[4332]811                     IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. epsi10 ) THEN
[921]812                        WRITE(numout,*) ' ALERTE 96 : wrong participation function ... '
813                        WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl
814                        WRITE(numout,*) ' lat, lon   : ', gphit(ji,jj), glamt(ji,jj)
815                        WRITE(numout,*) ' aridge     : ', aridge(ji,jj,1:jpl)
816                        WRITE(numout,*) ' araft      : ', araft(ji,jj,1:jpl)
817                        WRITE(numout,*) ' athorn     : ', athorn(ji,jj,1:jpl)
818                     ENDIF
819                  END DO
[868]820               END DO
[825]821            END DO
[921]822         ENDIF
[825]823
824      ENDIF
825
[921]826      !-----------------------------------------------------------------
827      ! 2) Transfer function
828      !-----------------------------------------------------------------
829      ! Compute max and min ridged ice thickness for each ridging category.
830      ! Assume ridged ice is uniformly distributed between hrmin and hrmax.
831      !
832      ! This parameterization is a modified version of Hibler (1980).
833      ! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
834      !  and for very thick ridging ice must be >= krdgmin*hi
835      !
836      ! The minimum ridging thickness, hrmin, is equal to 2*hi
837      !  (i.e., rafting) and for very thick ridging ice is
838      !  constrained by hrmin <= (hrmean + hi)/2.
839      !
840      ! The maximum ridging thickness, hrmax, is determined by
841      !  hrmean and hrmin.
842      !
843      ! These modifications have the effect of reducing the ice strength
844      ! (relative to the Hibler formulation) when very thick ice is
845      ! ridging.
846      !
847      ! aksum = net area removed/ total area removed
848      ! where total area removed = area of ice that ridges
849      !         net area removed = total area removed - area of new ridges
850      !-----------------------------------------------------------------
[825]851
852      ! Transfer function
853      DO jl = 1, jpl !all categories have a specific transfer function
854         DO jj = 1, jpj
855            DO ji = 1, jpi
856
[4332]857               IF (a_i(ji,jj,jl) .GT. epsi10 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN
[825]858                  hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
859                  hrmean          = MAX(SQRT(Hstar*hi), hi*krdgmin)
860                  hrmin(ji,jj,jl) = MIN(2.0*hi, 0.5*(hrmean + hi))
861                  hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl)
862                  hraft(ji,jj,jl) = kraft*hi
863                  krdg(ji,jj,jl)  = hrmean / hi
864               ELSE
865                  hraft(ji,jj,jl) = 0.0
866                  hrmin(ji,jj,jl) = 0.0 
867                  hrmax(ji,jj,jl) = 0.0 
868                  krdg (ji,jj,jl) = 1.0
869               ENDIF
870
871            END DO ! ji
872         END DO ! jj
873      END DO ! jl
874
875      ! Normalization factor : aksum, ensures mass conservation
[2777]876      aksum(:,:) = athorn(:,:,0)
[825]877      DO jl = 1, jpl 
[2715]878         aksum(:,:)    = aksum(:,:) + aridge(:,:,jl) * ( 1._wp - 1._wp / krdg(:,:,jl) )    &
879            &                       + araft (:,:,jl) * ( 1._wp - 1._wp / kraft        )
[825]880      END DO
[2715]881      !
[3294]882      CALL wrk_dealloc( jpi,jpj, zworka )
883      CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
884      !
[921]885   END SUBROUTINE lim_itd_me_ridgeprep
[825]886
887
[2715]888   SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
889      !!----------------------------------------------------------------------
[921]890      !!                ***  ROUTINE lim_itd_me_icestrength ***
891      !!
[2715]892      !! ** Purpose :   shift ridging ice among thickness categories of ice thickness
[921]893      !!
[2715]894      !! ** Method  :   Remove area, volume, and energy from each ridging category
895      !!              and add to thicker ice categories.
896      !!----------------------------------------------------------------------
897      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   opning         ! rate of opening due to divergence/shear
898      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area removed, excluding area of new ridges
899      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   msnow_mlt      ! mass of snow added to ocean (kg m-2)
900      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   esnow_mlt      ! energy needed to melt snow in ocean (J m-2)
901      !
902      CHARACTER (len=80) ::   fieldid   ! field identifier
903      LOGICAL, PARAMETER ::   l_conservation_check = .true.  ! if true, check conservation (useful for debugging)
904      !
905      LOGICAL ::   neg_ato_i      ! flag for ato_i(i,j) < -puny
906      LOGICAL ::   large_afrac    ! flag for afrac > 1
907      LOGICAL ::   large_afrft    ! flag for afrac > 1
908      INTEGER ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
909      INTEGER ::   ij                ! horizontal index, combines i and j loops
910      INTEGER ::   icells            ! number of cells with aicen > puny
[4506]911      REAL(wp) ::   zindb, zsrdg2    ! local scalar
[2715]912      REAL(wp) ::   hL, hR, farea, zdummy, zdummy0, ztmelts    ! left and right limits of integration
[4506]913      REAL(wp) ::   zsstK            ! SST in Kelvin
[825]914
[3294]915      INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices
[825]916
[3294]917      REAL(wp), POINTER, DIMENSION(:,:) ::   vice_init, vice_final   ! ice volume summed over categories
918      REAL(wp), POINTER, DIMENSION(:,:) ::   eice_init, eice_final   ! ice energy summed over layers
[825]919
[3294]920      REAL(wp), POINTER, DIMENSION(:,:,:) ::   aicen_init, vicen_init   ! ice  area    & volume before ridging
921      REAL(wp), POINTER, DIMENSION(:,:,:) ::   vsnon_init, esnon_init   ! snow volume  & energy before ridging
922      REAL(wp), POINTER, DIMENSION(:,:,:) ::   smv_i_init, oa_i_init    ! ice salinity & age    before ridging
[825]923
[3294]924      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   eicen_init        ! ice energy before ridging
[825]925
[3294]926      REAL(wp), POINTER, DIMENSION(:,:) ::   afrac , fvol     ! fraction of category area ridged & new ridge volume going to n2
927      REAL(wp), POINTER, DIMENSION(:,:) ::   ardg1 , ardg2    ! area of ice ridged & new ridges
928      REAL(wp), POINTER, DIMENSION(:,:) ::   vsrdg , esrdg    ! snow volume & energy of ridging ice
929      REAL(wp), POINTER, DIMENSION(:,:) ::   oirdg1, oirdg2   ! areal age content of ridged & rifging ice
930      REAL(wp), POINTER, DIMENSION(:,:) ::   dhr   , dhr2     ! hrmax - hrmin  &  hrmax^2 - hrmin^2
[825]931
[3294]932      REAL(wp), POINTER, DIMENSION(:,:) ::   vrdg1   ! volume of ice ridged
933      REAL(wp), POINTER, DIMENSION(:,:) ::   vrdg2   ! volume of new ridges
934      REAL(wp), POINTER, DIMENSION(:,:) ::   vsw     ! volume of seawater trapped into ridges
935      REAL(wp), POINTER, DIMENSION(:,:) ::   srdg1   ! sal*volume of ice ridged
936      REAL(wp), POINTER, DIMENSION(:,:) ::   srdg2   ! sal*volume of new ridges
937      REAL(wp), POINTER, DIMENSION(:,:) ::   smsw    ! sal*volume of water trapped into ridges
[825]938
[3294]939      REAL(wp), POINTER, DIMENSION(:,:) ::   afrft            ! fraction of category area rafted
940      REAL(wp), POINTER, DIMENSION(:,:) ::   arft1 , arft2    ! area of ice rafted and new rafted zone
941      REAL(wp), POINTER, DIMENSION(:,:) ::   virft , vsrft    ! ice & snow volume of rafting ice
942      REAL(wp), POINTER, DIMENSION(:,:) ::   esrft , smrft    ! snow energy & salinity of rafting ice
943      REAL(wp), POINTER, DIMENSION(:,:) ::   oirft1, oirft2   ! areal age content of rafted ice & rafting ice
[825]944
[3294]945      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eirft      ! ice energy of rafting ice
946      REAL(wp), POINTER, DIMENSION(:,:,:) ::   erdg1      ! enth*volume of ice ridged
947      REAL(wp), POINTER, DIMENSION(:,:,:) ::   erdg2      ! enth*volume of new ridges
948      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ersw       ! enth of water trapped into ridges
949      !!----------------------------------------------------------------------
[825]950
[3294]951      CALL wrk_alloc( (jpi+1)*(jpj+1),      indxi, indxj )
952      CALL wrk_alloc( jpi, jpj,             vice_init, vice_final, eice_init, eice_final )
953      CALL wrk_alloc( jpi, jpj,             afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )
954      CALL wrk_alloc( jpi, jpj,             vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw )
955      CALL wrk_alloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
956      CALL wrk_alloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnon_init, esnon_init, smv_i_init, oa_i_init )
957      CALL wrk_alloc( jpi, jpj, jkmax,      eirft, erdg1, erdg2, ersw )
958      CALL wrk_alloc( jpi, jpj, jkmax, jpl, eicen_init )
959
[825]960      ! Conservation check
[2715]961      eice_init(:,:) = 0._wp
[825]962
[2715]963      IF( con_i ) THEN
[4332]964         CALL lim_column_sum        (jpl,    v_i,       vice_init )
[825]965         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_init )
[4332]966         DO ji = mi0(jiindx), mi1(jiindx)
967            DO jj = mj0(jjindx), mj1(jjindx)
968               WRITE(numout,*) ' vice_init  : ', vice_init(ji,jj)
969               WRITE(numout,*) ' eice_init  : ', eice_init(ji,jj)
970            END DO
971         END DO
[825]972      ENDIF
973
[921]974      !-------------------------------------------------------------------------------
975      ! 1) Compute change in open water area due to closing and opening.
976      !-------------------------------------------------------------------------------
[825]977
978      neg_ato_i = .false.
979
980      DO jj = 1, jpj
981         DO ji = 1, jpi
[2715]982            ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice        &
983               &                        + opning(ji,jj)                          * rdt_ice
[4332]984            IF( ato_i(ji,jj) < -epsi10 ) THEN
[2715]985               neg_ato_i = .TRUE.
986            ELSEIF( ato_i(ji,jj) < 0._wp ) THEN    ! roundoff error
987               ato_i(ji,jj) = 0._wp
[825]988            ENDIF
989         END DO !jj
990      END DO !ji
991
992      ! if negative open water area alert it
[2715]993      IF( neg_ato_i ) THEN       ! there is a bug
[825]994         DO jj = 1, jpj 
995            DO ji = 1, jpi
[4332]996               IF( ato_i(ji,jj) < -epsi10 ) THEN
[825]997                  WRITE(numout,*) '' 
998                  WRITE(numout,*) 'Ridging error: ato_i < 0'
999                  WRITE(numout,*) 'ato_i : ', ato_i(ji,jj)
[4332]1000               ENDIF               ! ato_i < -epsi10
[2715]1001            END DO
1002         END DO
1003      ENDIF
[825]1004
[921]1005      !-----------------------------------------------------------------
1006      ! 2) Save initial state variables
1007      !-----------------------------------------------------------------
[825]1008
1009      DO jl = 1, jpl
[2715]1010         aicen_init(:,:,jl) = a_i(:,:,jl)
1011         vicen_init(:,:,jl) = v_i(:,:,jl)
1012         vsnon_init(:,:,jl) = v_s(:,:,jl)
1013         !
1014         smv_i_init(:,:,jl) = smv_i(:,:,jl)
1015         oa_i_init (:,:,jl) = oa_i (:,:,jl)
[825]1016      END DO !jl
[868]1017
1018      esnon_init(:,:,:) = e_s(:,:,1,:)
[921]1019
[825]1020      DO jl = 1, jpl 
1021         DO jk = 1, nlay_i
[2715]1022            eicen_init(:,:,jk,jl) = e_i(:,:,jk,jl)
1023         END DO
1024      END DO
[825]1025
[921]1026      !
1027      !-----------------------------------------------------------------
1028      ! 3) Pump everything from ice which is being ridged / rafted
1029      !-----------------------------------------------------------------
1030      ! Compute the area, volume, and energy of ice ridging in each
1031      ! category, along with the area of the resulting ridge.
[825]1032
1033      DO jl1 = 1, jpl !jl1 describes the ridging category
1034
[921]1035         !------------------------------------------------
1036         ! 3.1) Identify grid cells with nonzero ridging
1037         !------------------------------------------------
[825]1038
1039         icells = 0
1040         DO jj = 1, jpj
1041            DO ji = 1, jpi
[4332]1042               IF( aicen_init(ji,jj,jl1) > epsi10 .AND. athorn(ji,jj,jl1) > 0._wp  &
1043                  &   .AND. closing_gross(ji,jj) > 0._wp ) THEN
[825]1044                  icells = icells + 1
1045                  indxi(icells) = ji
1046                  indxj(icells) = jj
1047               ENDIF ! test on a_icen_init
1048            END DO ! ji
1049         END DO ! jj
1050
1051         large_afrac = .false.
1052         large_afrft = .false.
1053
[868]1054!CDIR NODEP
[825]1055         DO ij = 1, icells
1056            ji = indxi(ij)
1057            jj = indxj(ij)
1058
[921]1059            !--------------------------------------------------------------------
1060            ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
1061            !--------------------------------------------------------------------
[825]1062
1063            ardg1(ji,jj) = aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice
1064            arft1(ji,jj) = araft (ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice
1065            ardg2(ji,jj) = ardg1(ji,jj) / krdg(ji,jj,jl1)
1066            arft2(ji,jj) = arft1(ji,jj) / kraft
1067
1068            oirdg1(ji,jj)= aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice
1069            oirft1(ji,jj)= araft (ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice
1070            oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1)
1071            oirft2(ji,jj)= oirft1(ji,jj) / kraft
1072
[921]1073            !---------------------------------------------------------------
1074            ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1
1075            !---------------------------------------------------------------
[825]1076
1077            afrac(ji,jj) = ardg1(ji,jj) / aicen_init(ji,jj,jl1) !ridging
1078            afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting
1079
[4332]1080            IF (afrac(ji,jj) > kamax + epsi10) THEN  !riging
[825]1081               large_afrac = .true.
[4045]1082            ELSEIF (afrac(ji,jj) > kamax) THEN  ! roundoff error
1083               afrac(ji,jj) = kamax
[825]1084            ENDIF
[4332]1085            IF (afrft(ji,jj) > kamax + epsi10) THEN !rafting
[825]1086               large_afrft = .true.
[4045]1087            ELSEIF (afrft(ji,jj) > kamax) THEN  ! roundoff error
1088               afrft(ji,jj) = kamax
[825]1089            ENDIF
1090
[921]1091            !--------------------------------------------------------------------------
1092            ! 3.4) Subtract area, volume, and energy from ridging
1093            !     / rafting category n1.
1094            !--------------------------------------------------------------------------
[2715]1095            vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por )
[825]1096            vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + ridge_por )
1097            vsw  (ji,jj) = vrdg1(ji,jj) * ridge_por
1098
1099            vsrdg(ji,jj) = vsnon_init(ji,jj,jl1) * afrac(ji,jj)
1100            esrdg(ji,jj) = esnon_init(ji,jj,jl1) * afrac(ji,jj)
[2715]1101            srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por )
[4506]1102            srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) !! MV HC 2014 this line seems useless
[825]1103
1104            ! rafting volumes, heat contents ...
1105            virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj)
1106            vsrft(ji,jj) = vsnon_init(ji,jj,jl1) * afrft(ji,jj)
1107            esrft(ji,jj) = esnon_init(ji,jj,jl1) * afrft(ji,jj)
1108            smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj) 
1109
1110            ! substract everything
1111            a_i(ji,jj,jl1)   = a_i(ji,jj,jl1)   - ardg1(ji,jj)  - arft1(ji,jj)
1112            v_i(ji,jj,jl1)   = v_i(ji,jj,jl1)   - vrdg1(ji,jj)  - virft(ji,jj)
1113            v_s(ji,jj,jl1)   = v_s(ji,jj,jl1)   - vsrdg(ji,jj)  - vsrft(ji,jj)
1114            e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg(ji,jj)  - esrft(ji,jj)
1115            oa_i(ji,jj,jl1)  = oa_i(ji,jj,jl1)  - oirdg1(ji,jj) - oirft1(ji,jj)
1116            smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1(ji,jj)  - smrft(ji,jj)
1117
[921]1118            !-----------------------------------------------------------------
1119            ! 3.5) Compute properties of new ridges
1120            !-----------------------------------------------------------------
[825]1121            !-------------
1122            ! Salinity
1123            !-------------
[4506]1124            smsw(ji,jj)  = vsw(ji,jj) * sss_m(ji,jj)                      ! salt content of seawater frozen in voids !! MV HC2014
[1571]1125            zsrdg2       = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge
[1572]1126
[1571]1127            srdg2(ji,jj) = MIN( s_i_max * vrdg2(ji,jj) , zsrdg2 )         ! impose a maximum salinity
[1572]1128           
[1571]1129            !                                                             ! excess of salt is flushed into the ocean
[4045]1130            !sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice
[4506]1131            ! MV HC 2014 this previous line seems ok too... Maybe units are not that good
[1572]1132
[4045]1133            !rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic    ! gurvan: increase in ice volume du to seawater frozen in voids             
[4506]1134            ! MV HC 2014 this previous line seems ok, i'm not sure at this moment of the sign convention
[4045]1135
[921]1136            !------------------------------------           
1137            ! 3.6 Increment ridging diagnostics
1138            !------------------------------------           
[825]1139
[921]1140            !        jl1 looping 1-jpl
1141            !           ij looping 1-icells
[825]1142
[2715]1143            dardg1dt   (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj)
1144            dardg2dt   (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj)
[3625]1145            opening    (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice
[825]1146
[2715]1147            IF( con_i )   vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj)
[825]1148
[921]1149            !------------------------------------------           
1150            ! 3.7 Put the snow somewhere in the ocean
1151            !------------------------------------------           
1152            !  Place part of the snow lost by ridging into the ocean.
1153            !  Note that esnow_mlt < 0; the ocean must cool to melt snow.
1154            !  If the ocean temp = Tf already, new ice must grow.
1155            !  During the next time step, thermo_rates will determine whether
1156            !  the ocean cools or new ice grows.
1157            !        jl1 looping 1-jpl
1158            !           ij looping 1-icells
1159
[2715]1160            msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosn*vsrdg(ji,jj)*(1.0-fsnowrdg)   &   ! rafting included
1161               &                                + rhosn*vsrft(ji,jj)*(1.0-fsnowrft)
[825]1162
[4506]1163            ! MV HC 2014 - I think the units are wrong here... shit there is a 10-9 in the snow...
1164            ! we should therefore do something
[2715]1165            esnow_mlt(ji,jj) = esnow_mlt(ji,jj) + esrdg(ji,jj)*(1.0-fsnowrdg)         &   !rafting included
1166               &                                + esrft(ji,jj)*(1.0-fsnowrft)         
[825]1167
[921]1168            !-----------------------------------------------------------------
1169            ! 3.8 Compute quantities used to apportion ice among categories
1170            ! in the n2 loop below
1171            !-----------------------------------------------------------------
[825]1172
[921]1173            !        jl1 looping 1-jpl
1174            !           ij looping 1-icells
[825]1175
[3625]1176            dhr (ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1)
[2715]1177            dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1)
[825]1178
1179         END DO                 ! ij
1180
[921]1181         !--------------------------------------------------------------------
1182         ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
1183         !      compute ridged ice enthalpy
1184         !--------------------------------------------------------------------
[825]1185         DO jk = 1, nlay_i
[868]1186!CDIR NODEP
[825]1187            DO ij = 1, icells
[921]1188               ji = indxi(ij)
1189               jj = indxj(ij)
1190               ! heat content of ridged ice
[2715]1191               erdg1(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por ) 
[921]1192               eirft(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj)
[2715]1193               e_i  (ji,jj,jk,jl1)  = e_i(ji,jj,jk,jl1) - erdg1(ji,jj,jk) - eirft(ji,jj,jk)
[4506]1194               
1195               
1196               ! enthalpy of the trapped seawater (J/m2, >0)
1197               zsstK  = sst_m(ji,jj) + rt0
1198               ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 )
[825]1199
[4506]1200               ! heat flux to the ocean
1201               rdq_ice(ji,jj) = rdq_ice(ji,jj) + ersw(ji,jj,jk) !! MV HC 2014
[825]1202
[921]1203               ! Correct dimensions to avoid big values
[2715]1204               ersw(ji,jj,jk)   = ersw(ji,jj,jk) * 1.e-09
[825]1205
[921]1206               ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
[4506]1207               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean
1208               !! MV HC 2014
[4045]1209               ersw (ji,jj,jk)  = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / REAL( nlay_i )
[825]1210
[921]1211               erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk)
[4506]1212
[825]1213            END DO ! ij
1214         END DO !jk
1215
1216
[2715]1217         IF( con_i ) THEN
[825]1218            DO jk = 1, nlay_i
[868]1219!CDIR NODEP
[825]1220               DO ij = 1, icells
1221                  ji = indxi(ij)
1222                  jj = indxj(ij)
[2715]1223                  eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - erdg1(ji,jj,jk)
[825]1224               END DO ! ij
1225            END DO !jk
1226         ENDIF
1227
[2715]1228         IF( large_afrac ) THEN   ! there is a bug
[868]1229!CDIR NODEP
[825]1230            DO ij = 1, icells
1231               ji = indxi(ij)
1232               jj = indxj(ij)
[4332]1233               IF( afrac(ji,jj) > kamax + epsi10 ) THEN
[825]1234                  WRITE(numout,*) ''
1235                  WRITE(numout,*) ' ardg > a_i'
[2715]1236                  WRITE(numout,*) ' ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1)
1237               ENDIF
1238            END DO
1239         ENDIF
1240         IF( large_afrft ) THEN  ! there is a bug
[868]1241!CDIR NODEP
[825]1242            DO ij = 1, icells
1243               ji = indxi(ij)
1244               jj = indxj(ij)
[4332]1245               IF( afrft(ji,jj) > kamax + epsi10 ) THEN
[825]1246                  WRITE(numout,*) ''
1247                  WRITE(numout,*) ' arft > a_i'
[2715]1248                  WRITE(numout,*) ' arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1)
1249               ENDIF
1250            END DO
1251         ENDIF
[825]1252
[921]1253         !-------------------------------------------------------------------------------
1254         ! 4) Add area, volume, and energy of new ridge to each category jl2
1255         !-------------------------------------------------------------------------------
1256         !        jl1 looping 1-jpl
[825]1257         DO jl2  = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 
[921]1258            ! over categories to which ridged ice is transferred
[868]1259!CDIR NODEP
[825]1260            DO ij = 1, icells
1261               ji = indxi(ij)
1262               jj = indxj(ij)
1263
1264               ! Compute the fraction of ridged ice area and volume going to
1265               ! thickness category jl2.
1266               ! Transfer area, volume, and energy accordingly.
1267
[3625]1268               IF( hrmin(ji,jj,jl1) >= hi_max(jl2)  .OR.        &
1269                   hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN
1270                  hL = 0._wp
1271                  hR = 0._wp
[825]1272               ELSE
[3625]1273                  hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
1274                  hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2)   )
[825]1275               ENDIF
1276
1277               ! fraction of ridged ice area and volume going to n2
[3625]1278               farea = ( hR - hL ) / dhr(ji,jj) 
1279               fvol(ji,jj) = ( hR*hR - hL*hL ) / dhr2(ji,jj)
[825]1280
[3625]1281               a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ardg2 (ji,jj) * farea
1282               v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + vrdg2 (ji,jj) * fvol(ji,jj)
1283               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * fsnowrdg
[2715]1284               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * fsnowrdg
[3625]1285               smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + srdg2 (ji,jj) * fvol(ji,jj)
1286               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirdg2(ji,jj) * farea
[825]1287
1288            END DO ! ij
1289
1290            ! Transfer ice energy to category jl2 by ridging
1291            DO jk = 1, nlay_i
[868]1292!CDIR NODEP
[825]1293               DO ij = 1, icells
1294                  ji = indxi(ij)
1295                  jj = indxj(ij)
[2715]1296                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj)*erdg2(ji,jj,jk)
1297               END DO
1298            END DO
1299            !
[825]1300         END DO                 ! jl2 (new ridges)           
1301
[2715]1302         DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 
[825]1303
[868]1304!CDIR NODEP
[825]1305            DO ij = 1, icells
1306               ji = indxi(ij)
1307               jj = indxj(ij)
1308               ! Compute the fraction of rafted ice area and volume going to
1309               ! thickness category jl2, transfer area, volume, and energy accordingly.
[3625]1310               !
1311               IF( hraft(ji,jj,jl1) <= hi_max(jl2)  .AND.        &
1312                   hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN
1313                  a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + arft2 (ji,jj)
1314                  v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + virft (ji,jj)
1315                  v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrft (ji,jj) * fsnowrft
1316                  e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrft (ji,jj) * fsnowrft
1317                  smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + smrft (ji,jj)   
1318                  oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirft2(ji,jj)   
[825]1319               ENDIF ! hraft
[3625]1320               !
[825]1321            END DO ! ij
1322
1323            ! Transfer rafted ice energy to category jl2
1324            DO jk = 1, nlay_i
[868]1325!CDIR NODEP
[825]1326               DO ij = 1, icells
1327                  ji = indxi(ij)
1328                  jj = indxj(ij)
[3625]1329                  IF(  hraft(ji,jj,jl1)  <=  hi_max(jl2)   .AND.        &
1330                       hraft(ji,jj,jl1)  >   hi_max(jl2-1)  ) THEN
[2715]1331                     e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk)
[825]1332                  ENDIF
1333               END DO           ! ij
1334            END DO !jk
1335
1336         END DO ! jl2
1337
1338      END DO ! jl1 (deforming categories)
1339
1340      ! Conservation check
1341      IF ( con_i ) THEN
1342         CALL lim_column_sum (jpl,   v_i, vice_final)
1343         fieldid = ' v_i : limitd_me '
1344         CALL lim_cons_check (vice_init, vice_final, 1.0e-6, fieldid) 
1345
1346         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_final )
1347         fieldid = ' e_i : limitd_me '
1348         CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid) 
[4332]1349
1350         DO ji = mi0(jiindx), mi1(jiindx)
1351            DO jj = mj0(jjindx), mj1(jjindx)
1352               WRITE(numout,*) ' vice_init  : ', vice_init (ji,jj)
1353               WRITE(numout,*) ' vice_final : ', vice_final(ji,jj)
1354               WRITE(numout,*) ' eice_init  : ', eice_init (ji,jj)
1355               WRITE(numout,*) ' eice_final : ', eice_final(ji,jj)
1356            END DO
1357         END DO
[825]1358      ENDIF
[2715]1359      !
[3294]1360      CALL wrk_dealloc( (jpi+1)*(jpj+1),      indxi, indxj )
1361      CALL wrk_dealloc( jpi, jpj,             vice_init, vice_final, eice_init, eice_final )
1362      CALL wrk_dealloc( jpi, jpj,             afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )
1363      CALL wrk_dealloc( jpi, jpj,             vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw )
1364      CALL wrk_dealloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
1365      CALL wrk_dealloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnon_init, esnon_init, smv_i_init, oa_i_init )
1366      CALL wrk_dealloc( jpi, jpj, jkmax,      eirft, erdg1, erdg2, ersw )
1367      CALL wrk_dealloc( jpi, jpj, jkmax, jpl, eicen_init )
1368      !
[825]1369   END SUBROUTINE lim_itd_me_ridgeshift
1370
1371
[2715]1372   SUBROUTINE lim_itd_me_asumr
[921]1373      !!-----------------------------------------------------------------------------
1374      !!                ***  ROUTINE lim_itd_me_asumr ***
1375      !!
[2715]1376      !! ** Purpose :   finds total fractional area
[921]1377      !!
[2715]1378      !! ** Method  :   Find the total area of ice plus open water in each grid cell.
1379      !!              This is similar to the aggregate_area subroutine except that the
1380      !!              total area can be greater than 1, so the open water area is
1381      !!              included in the sum instead of being computed as a residual.
1382      !!-----------------------------------------------------------------------------
1383      INTEGER ::   jl   ! dummy loop index
1384      !!-----------------------------------------------------------------------------
1385      !
1386      asum(:,:) = ato_i(:,:)                    ! open water
1387      DO jl = 1, jpl                            ! ice categories
1388         asum(:,:) = asum(:,:) + a_i(:,:,jl)
[825]1389      END DO
[2715]1390      !
[825]1391   END SUBROUTINE lim_itd_me_asumr
1392
1393
[2715]1394   SUBROUTINE lim_itd_me_init
[825]1395      !!-------------------------------------------------------------------
1396      !!                   ***  ROUTINE lim_itd_me_init ***
1397      !!
1398      !! ** Purpose :   Physical constants and parameters linked
1399      !!                to the mechanical ice redistribution
1400      !!
1401      !! ** Method  :   Read the namiceitdme namelist
1402      !!                and check the parameters values
1403      !!                called at the first timestep (nit000)
1404      !!
1405      !! ** input   :   Namelist namiceitdme
1406      !!-------------------------------------------------------------------
1407      NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft,& 
[921]1408         Gstar, astar,                                &
1409         Hstar, raftswi, hparmeter, Craft, ridge_por, &
1410         sal_max_ridge,  partfun_swi, transfun_swi,   &
1411         brinstren_swi
[825]1412      !!-------------------------------------------------------------------
[2715]1413      !
1414      REWIND( numnam_ice )                   ! read namiceitdme namelist
[825]1415      READ  ( numnam_ice , namiceitdme)
[2715]1416      !
1417      IF (lwp) THEN                          ! control print
[825]1418         WRITE(numout,*)
1419         WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution '
1420         WRITE(numout,*)' ~~~~~~~~~~~~~~~'
1421         WRITE(numout,*)'   Switch choosing the ice redistribution scheme           ridge_scheme_swi', ridge_scheme_swi 
1422         WRITE(numout,*)'   Fraction of shear energy contributing to ridging        Cs              ', Cs 
1423         WRITE(numout,*)'   Ratio of ridging work to PotEner change in ridging      Cf              ', Cf 
1424         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        fsnowrdg        ', fsnowrdg 
1425         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        fsnowrft        ', fsnowrft 
1426         WRITE(numout,*)'   Fraction of total ice coverage contributing to ridging  Gstar           ', Gstar
1427         WRITE(numout,*)'   Equivalent to G* for an exponential part function       astar           ', astar
1428         WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     Hstar           ', Hstar
1429         WRITE(numout,*)'   Rafting of ice sheets or not                            raftswi         ', raftswi
1430         WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       hparmeter       ', hparmeter
1431         WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  Craft           ', Craft 
1432         WRITE(numout,*)'   Initial porosity of ridges                              ridge_por       ', ridge_por
1433         WRITE(numout,*)'   Maximum salinity of ridging ice                         sal_max_ridge   ', sal_max_ridge
1434         WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    partfun_swi     ', partfun_swi
1435         WRITE(numout,*)'   Switch for tran. function (0) linear (1) exponential    transfun_swi    ', transfun_swi
1436         WRITE(numout,*)'   Switch for including brine volume in ice strength comp. brinstren_swi   ', brinstren_swi
1437      ENDIF
[2715]1438      !
[825]1439   END SUBROUTINE lim_itd_me_init
1440
1441
1442   SUBROUTINE lim_itd_me_zapsmall
1443      !!-------------------------------------------------------------------
1444      !!                   ***  ROUTINE lim_itd_me_zapsmall ***
1445      !!
1446      !! ** Purpose :   Remove too small sea ice areas and correct salt fluxes
1447      !!
1448      !! history :
1449      !! author: William H. Lipscomb, LANL
1450      !! Nov 2003:  Modified by Julie Schramm to conserve volume and energy
1451      !! Sept 2004: Modified by William Lipscomb; replaced normalize_state with
1452      !!            additions to local freshwater, salt, and heat fluxes
1453      !!  9.0, LIM3.0 - 02-2006 (M. Vancoppenolle) original code
1454      !!-------------------------------------------------------------------
[2715]1455      INTEGER ::   ji, jj, jl, jk   ! dummy loop indices
1456      INTEGER ::   icells           ! number of cells with ice to zap
[825]1457
[3294]1458      REAL(wp), POINTER, DIMENSION(:,:) ::   zmask   ! 2D workspace
[4332]1459      REAL(wp)                          ::   zmask_glo
[2715]1460!!gm      REAL(wp) ::   xtmp      ! temporary variable
1461      !!-------------------------------------------------------------------
[825]1462
[3294]1463      CALL wrk_alloc( jpi, jpj, zmask )
1464
[825]1465      DO jl = 1, jpl
1466
[921]1467         !-----------------------------------------------------------------
1468         ! Count categories to be zapped.
1469         ! Abort model in case of negative area.
1470         !-----------------------------------------------------------------
1471         icells = 0
[4332]1472         zmask(:,:)  = 0._wp
[921]1473         DO jj = 1, jpj
1474            DO ji = 1, jpi
[4332]1475               IF(  ( a_i(ji,jj,jl) >= -epsi10 .AND. a_i(ji,jj,jl) <  0._wp   ) .OR.   &
1476                  & ( a_i(ji,jj,jl) >  0._wp   .AND. a_i(ji,jj,jl) <= epsi10  ) .OR.   &
1477                  & ( v_i(ji,jj,jl) == 0._wp   .AND. a_i(ji,jj,jl) >  0._wp   ) .OR.   &
1478                  & ( v_i(ji,jj,jl) >  0._wp   .AND. v_i(ji,jj,jl) <= epsi10  ) )   zmask(ji,jj) = 1._wp
[921]1479            END DO
[825]1480         END DO
[4345]1481         !zmask_glo = glob_sum(zmask)
[4332]1482         !IF( ln_nicep .AND. lwp ) WRITE(numout,*) zmask_glo, ' cells of ice zapped in the ocean '
[825]1483
[921]1484         !-----------------------------------------------------------------
1485         ! Zap ice energy and use ocean heat to melt ice
1486         !-----------------------------------------------------------------
[825]1487
1488         DO jk = 1, nlay_i
[868]1489            DO jj = 1 , jpj
1490               DO ji = 1 , jpi
[3625]1491!!gm                  xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) * r1_rdtice
[2715]1492!!gm                  xtmp = xtmp * unit_fac
1493                  ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp
[4332]1494                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) )
[2715]1495               END DO
1496            END DO
1497         END DO
[825]1498
[868]1499         DO jj = 1 , jpj
1500            DO ji = 1 , jpi
[825]1501
[921]1502               !-----------------------------------------------------------------
1503               ! Zap snow energy and use ocean heat to melt snow
1504               !-----------------------------------------------------------------
1505               !           xtmp = esnon(i,j,n) / dt ! < 0
1506               !           fhnet(i,j)      = fhnet(i,j)      + xtmp
1507               !           fhnet_hist(i,j) = fhnet_hist(i,j) + xtmp
1508               ! xtmp is greater than 0
1509               ! fluxes are positive to the ocean
1510               ! here the flux has to be negative for the ocean
[3625]1511!!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice
[921]1512               !           fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp
[825]1513
[3625]1514!!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice !RB   ???????
[825]1515
[4332]1516               t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) )
[868]1517
[921]1518               !-----------------------------------------------------------------
1519               ! zap ice and snow volume, add water and salt to ocean
1520               !-----------------------------------------------------------------
[825]1521
[921]1522               !           xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) / dt
[3625]1523               !           sfx_res(ji,jj)  = sfx_res(ji,jj) + ( sss_m(ji,jj)                  )   &
1524               !                                            * rhosn * v_s(ji,jj,jl) * r1_rdtice
1525               !           sfx_res(ji,jj)  = sfx_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) )   &
1526               !                                            * rhoic * v_i(ji,jj,jl) * r1_rdtice
1527               !           sfx (i,j)      = sfx (i,j)      + xtmp
[825]1528
[2715]1529               ato_i(ji,jj)    = a_i  (ji,jj,jl) *       zmask(ji,jj)   + ato_i(ji,jj)
[4332]1530               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
1531               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
1532               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
1533               t_su (ji,jj,jl) = t_su (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj)
1534               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
1535               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
[2715]1536               !
1537            END DO
1538         END DO
1539         !
[825]1540      END DO                 ! jl
[2715]1541      !
[3294]1542      CALL wrk_dealloc( jpi, jpj, zmask )
1543      !
[825]1544   END SUBROUTINE lim_itd_me_zapsmall
1545
1546#else
[2715]1547   !!----------------------------------------------------------------------
1548   !!   Default option         Empty module          NO LIM-3 sea-ice model
1549   !!----------------------------------------------------------------------
[825]1550CONTAINS
1551   SUBROUTINE lim_itd_me           ! Empty routines
1552   END SUBROUTINE lim_itd_me
1553   SUBROUTINE lim_itd_me_icestrength
1554   END SUBROUTINE lim_itd_me_icestrength
1555   SUBROUTINE lim_itd_me_sort
1556   END SUBROUTINE lim_itd_me_sort
1557   SUBROUTINE lim_itd_me_init
1558   END SUBROUTINE lim_itd_me_init
1559   SUBROUTINE lim_itd_me_zapsmall
1560   END SUBROUTINE lim_itd_me_zapsmall
1561#endif
[2715]1562   !!======================================================================
[825]1563END MODULE limitd_me
Note: See TracBrowser for help on using the repository browser.