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

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

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

Last change on this file since 3584 was 3558, checked in by rblod, 12 years ago

Fix issues when using key_nosignedzeo, see ticket #996

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