source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

  • Property svn:keywords set to Id
File size: 33.5 KB
Line 
1MODULE limitd_th
2   !!======================================================================
3   !!                       ***  MODULE limitd_th ***
4   !!   LIM3 ice model : ice thickness distribution: Thermodynamics
5   !!======================================================================
6   !! History :   -   !          (W. H. Lipscomb and E.C. Hunke) CICE (c) original code
7   !!            3.0  ! 2005-12  (M. Vancoppenolle) adaptation to LIM-3
8   !!             -   ! 2006-06  (M. Vancoppenolle) adaptation to include salt, age
9   !!             -   ! 2007-04  (M. Vancoppenolle) Mass conservation checked
10   !!----------------------------------------------------------------------
11#if defined key_lim3
12   !!----------------------------------------------------------------------
13   !!   'key_lim3' :                                   LIM3 sea-ice model
14   !!----------------------------------------------------------------------
15   !!   lim_itd_th_rem   :
16   !!   lim_itd_th_reb   :
17   !!   lim_itd_fitline  :
18   !!   lim_itd_shiftice :
19   !!----------------------------------------------------------------------
20   USE par_oce          ! ocean parameters
21   USE dom_oce          ! ocean domain
22   USE phycst           ! physical constants (ocean directory)
23   USE thd_ice          ! LIM-3 thermodynamic variables
24   USE ice              ! LIM-3 variables
25   USE limvar           ! LIM-3 variables
26   USE prtctl           ! Print control
27   USE in_out_manager   ! I/O manager
28   USE lib_mpp          ! MPP library
29   USE lib_fortran      ! to use key_nosignedzero
30   USE limcons          ! conservation tests
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   lim_itd_th_rem
36   PUBLIC   lim_itd_th_reb
37
38   !!----------------------------------------------------------------------
39   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
40   !! $Id$
41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt )
46      !!------------------------------------------------------------------
47      !!                ***  ROUTINE lim_itd_th_rem ***
48      !!
49      !! ** Purpose :   computes the redistribution of ice thickness
50      !!              after thermodynamic growth of ice thickness
51      !!
52      !! ** Method  : Linear remapping
53      !!
54      !! References : W.H. Lipscomb, JGR 2001
55      !!------------------------------------------------------------------
56      INTEGER , INTENT (in) ::   klbnd   ! Start thickness category index point
57      INTEGER , INTENT (in) ::   kubnd   ! End point on which the  the computation is applied
58      INTEGER , INTENT (in) ::   kt      ! Ocean time step
59      !
60      INTEGER  ::   ji, jj, jl     ! dummy loop index
61      INTEGER  ::   ii, ij         ! 2D corresponding indices to ji
62      INTEGER  ::   nd             ! local integer
63      REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars
64      REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      -
65      REAL(wp) ::   zx3       
66      CHARACTER (len = 15) :: fieldid
67
68      INTEGER , DIMENSION(jpi,jpj,jpl-1) ::   zdonor   ! donor category index
69
70      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zdhice      ! ice thickness increment
71      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   g0          ! coefficients for fitting the line of the ITD
72      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   g1          ! coefficients for fitting the line of the ITD
73      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   hL          ! left boundary for the ITD for each thickness
74      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   hR          ! left boundary for the ITD for each thickness
75      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zht_i_b     ! old ice thickness
76      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   dummy_es
77      REAL(wp), DIMENSION(jpi,jpj,jpl-1) ::   zdaice, zdvice          ! local increment of ice area and volume
78      REAL(wp), DIMENSION((jpi+1)*(jpj+1))     ::   zvetamin, zvetamax      ! maximum values for etas
79      INTEGER , DIMENSION((jpi+1)*(jpj+1))     ::   nind_i, nind_j          ! compressed indices for i/j directions
80      INTEGER                             ::   nbrem                   ! number of cells with ice to transfer
81      REAL(wp)                            ::   zslope                  ! used to compute local thermodynamic "speeds"
82      REAL(wp), DIMENSION(jpi,jpj)   ::   zhb0, zhb1              ! category boundaries for thinnes categories
83      REAL(wp), DIMENSION(jpi,jpj)   ::   vt_i_init, vt_i_final   !  ice volume summed over categories
84      REAL(wp), DIMENSION(jpi,jpj)   ::   vt_s_init, vt_s_final   !  snow volume summed over categories
85      REAL(wp), DIMENSION(jpi,jpj)   ::   et_i_init, et_i_final   !  ice energy summed over categories
86      REAL(wp), DIMENSION(jpi,jpj)   ::   et_s_init, et_s_final   !  snow energy summed over categories
87      INTEGER , DIMENSION(jpi,jpj)   ::   zremap_flag      ! compute remapping or not ????
88      REAL(wp), DIMENSION(jpi,jpj,0:jpl) ::   zhbnew           ! new boundaries of ice categories
89      !!------------------------------------------------------------------
90
91
92      !!----------------------------------------------------------------------------------------------
93      !! 0) Conservation checkand changes in each ice category
94      !!----------------------------------------------------------------------------------------------
95      IF( con_i ) THEN
96         CALL lim_column_sum (jpl,   v_i, vt_i_init)
97         CALL lim_column_sum (jpl,   v_s, vt_s_init)
98         CALL lim_column_sum_energy (jpl, nlay_i,   e_i, et_i_init)
99         dummy_es(:,:,:) = e_s(:,:,1,:)
100         CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_init)
101      ENDIF
102
103      !!----------------------------------------------------------------------------------------------
104      !! 1) Compute thickness and changes in each ice category
105      !!----------------------------------------------------------------------------------------------
106      IF( kt == nit000 .AND. lwp) THEN
107         WRITE(numout,*)
108         WRITE(numout,*) 'lim_itd_th_rem  : Remapping the ice thickness distribution'
109         WRITE(numout,*) '~~~~~~~~~~~~~~~'
110         WRITE(numout,*) ' klbnd :       ', klbnd
111         WRITE(numout,*) ' kubnd :       ', kubnd
112      ENDIF
113
114      zdhice(:,:,:) = 0._wp
115      DO jl = klbnd, kubnd
116         DO jj = 1, jpj
117            DO ji = 1, jpi
118               rswitch           = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi10 ) )     !0 if no ice and 1 if yes
119               ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch
120               rswitch           = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) )
121               zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch
122               IF( a_i(ji,jj,jl) > epsi10 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) ! clem: useless IF statement?
123            END DO
124         END DO
125      END DO
126
127      !-----------------------------------------------------------------------------------------------
128      !  2) Compute fractional ice area in each grid cell
129      !-----------------------------------------------------------------------------------------------
130      at_i(:,:) = 0._wp
131      DO jl = klbnd, kubnd
132         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
133      END DO
134
135      !-----------------------------------------------------------------------------------------------
136      !  3) Identify grid cells with ice
137      !-----------------------------------------------------------------------------------------------
138      nbrem = 0
139      DO jj = 1, jpj
140         DO ji = 1, jpi
141            IF ( at_i(ji,jj) > epsi10 ) THEN
142               nbrem         = nbrem + 1
143               nind_i(nbrem) = ji
144               nind_j(nbrem) = jj
145               zremap_flag(ji,jj) = 1
146            ELSE
147               zremap_flag(ji,jj) = 0
148            ENDIF
149         END DO
150      END DO
151
152      !-----------------------------------------------------------------------------------------------
153      !  4) Compute new category boundaries
154      !-----------------------------------------------------------------------------------------------
155      !- 4.1 Compute category boundaries
156      zhbnew(:,:,:) = 0._wp
157
158      DO jl = klbnd, kubnd - 1
159         DO ji = 1, nbrem
160            ii = nind_i(ji)
161            ij = nind_j(ji)
162            !
163            zhbnew(ii,ij,jl) = hi_max(jl)
164            IF    ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN
165               !interpolate between adjacent category growth rates
166               zslope           = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) )
167               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) )
168            ELSEIF( a_i_b(ii,ij,jl) > epsi10) THEN
169               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl)
170            ELSEIF( a_i_b(ii,ij,jl+1) > epsi10) THEN
171               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1)
172            ENDIF
173         END DO
174
175         !- 4.2 Check that each zhbnew lies between adjacent values of ice thickness
176         DO ji = 1, nbrem
177            ii = nind_i(ji)
178            ij = nind_j(ji)
179
180            ! clem: we do not want ht_i to be too close to either HR or HL otherwise a division by nearly 0 is possible
181            ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
182            IF    ( a_i(ii,ij,jl  ) > epsi10 .AND. ht_i(ii,ij,jl  ) > ( zhbnew(ii,ij,jl) - epsi10 ) ) THEN
183               zremap_flag(ii,ij) = 0
184            ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) ) THEN
185               zremap_flag(ii,ij) = 0
186            ENDIF
187
188            !- 4.3 Check that each zhbnew does not exceed maximal values hi_max 
189            IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0
190            IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0
191            ! clem bug: why is not the following instead?
192            !!IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0
193            !!IF( zhbnew(ii,ij,jl) > hi_max(jl  ) ) zremap_flag(ii,ij) = 0
194 
195         END DO
196
197      END DO
198
199      !-----------------------------------------------------------------------------------------------
200      !  5) Identify cells where ITD is to be remapped
201      !-----------------------------------------------------------------------------------------------
202      nbrem = 0
203      DO jj = 1, jpj
204         DO ji = 1, jpi
205            IF( zremap_flag(ji,jj) == 1 ) THEN
206               nbrem         = nbrem + 1
207               nind_i(nbrem) = ji
208               nind_j(nbrem) = jj
209            ENDIF
210         END DO
211      END DO 
212
213      !-----------------------------------------------------------------------------------------------
214      !  6) Fill arrays with lowermost / uppermost boundaries of 'new' categories
215      !-----------------------------------------------------------------------------------------------
216      DO jj = 1, jpj
217         DO ji = 1, jpi
218            zhb0(ji,jj) = hi_max(0)
219            zhb1(ji,jj) = hi_max(1)
220
221            IF( a_i(ji,jj,kubnd) > epsi10 ) THEN
222               zhbnew(ji,jj,kubnd) = MAX( hi_max(kubnd-1), 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) )
223            ELSE
224!clem bug               zhbnew(ji,jj,kubnd) = hi_max(kubnd) 
225               zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) ! not used anyway
226            ENDIF
227
228            ! clem: we do not want ht_i_b to be too close to either HR or HL otherwise a division by nearly 0 is possible
229            ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
230            IF    ( zht_i_b(ji,jj,klbnd) < ( zhb0(ji,jj) + epsi10 ) )  THEN
231               zremap_flag(ji,jj) = 0
232            ELSEIF( zht_i_b(ji,jj,klbnd) > ( zhb1(ji,jj) - epsi10 ) )  THEN
233               zremap_flag(ji,jj) = 0
234            ENDIF
235
236         END DO
237      END DO
238
239      !-----------------------------------------------------------------------------------------------
240      !  7) Compute g(h)
241      !-----------------------------------------------------------------------------------------------
242      !- 7.1 g(h) for category 1 at start of time step
243      CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd),   &
244         &                  hR(:,:,klbnd), zremap_flag )
245
246      !- 7.2 Area lost due to melting of thin ice (first category,  klbnd)
247      DO ji = 1, nbrem
248         ii = nind_i(ji) 
249         ij = nind_j(ji) 
250
251         IF( a_i(ii,ij,klbnd) > epsi10 ) THEN
252
253            zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category
254
255            IF( zdh0 < 0.0 ) THEN      !remove area from category 1
256               zdh0 = MIN( -zdh0, hi_max(klbnd) )
257               !Integrate g(1) from 0 to dh0 to estimate area melted
258               zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd)
259
260               IF( zetamax > 0.0 ) THEN
261                  zx1    = zetamax
262                  zx2    = 0.5 * zetamax * zetamax 
263                  zda0   = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1                        ! ice area removed
264                  zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! Constrain new thickness <= ht_i               
265                  zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting
266                                                                                                !     of thin ice (zdamax > 0)
267                  ! Remove area, conserving volume
268                  ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 )
269                  a_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) - zda0
270                  v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ?
271               ENDIF
272
273            ELSE ! if ice accretion zdh0 > 0
274               ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1)
275               zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) ) 
276            ENDIF
277
278         ENDIF
279
280      END DO
281
282      !- 7.3 g(h) for each thickness category 
283      DO jl = klbnd, kubnd
284         CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl),  &
285            &                  g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag )
286      END DO
287
288      !-----------------------------------------------------------------------------------------------
289      !  8) Compute area and volume to be shifted across each boundary
290      !-----------------------------------------------------------------------------------------------
291
292      DO jl = klbnd, kubnd - 1
293         DO jj = 1, jpj
294            DO ji = 1, jpi
295               zdonor(ji,jj,jl) = 0
296               zdaice(ji,jj,jl) = 0.0
297               zdvice(ji,jj,jl) = 0.0
298            END DO
299         END DO
300
301         DO ji = 1, nbrem
302            ii = nind_i(ji)
303            ij = nind_j(ji)
304
305            IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1
306               ! left and right integration limits in eta space
307               zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl)
308               zvetamax(ji) = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl)
309               zdonor(ii,ij,jl) = jl
310
311            ELSE                                    ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl
312               ! left and right integration limits in eta space
313               zvetamin(ji) = 0.0
314               zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1)
315               zdonor(ii,ij,jl) = jl + 1
316
317            ENDIF
318
319            zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin
320            zetamin = zvetamin(ji)
321
322            zx1  = zetamax - zetamin
323            zwk1 = zetamin * zetamin
324            zwk2 = zetamax * zetamax
325            zx2  = 0.5 * ( zwk2 - zwk1 )
326            zwk1 = zwk1 * zetamin
327            zwk2 = zwk2 * zetamax
328            zx3  = 1.0 / 3.0 * ( zwk2 - zwk1 )
329            nd   = zdonor(ii,ij,jl)
330            zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1
331            zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd)
332
333         END DO
334      END DO
335
336      !!----------------------------------------------------------------------------------------------
337      !! 9) Shift ice between categories
338      !!----------------------------------------------------------------------------------------------
339      CALL lim_itd_shiftice ( klbnd, kubnd, zdonor, zdaice, zdvice )
340
341      !!----------------------------------------------------------------------------------------------
342      !! 10) Make sure ht_i >= minimum ice thickness hi_min
343      !!----------------------------------------------------------------------------------------------
344
345      DO ji = 1, nbrem
346         ii = nind_i(ji)
347         ij = nind_j(ji)
348         IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN
349            a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin 
350            ht_i(ii,ij,1) = rn_himin
351         ENDIF
352      END DO
353
354      !!----------------------------------------------------------------------------------------------
355      !! 11) Conservation check
356      !!----------------------------------------------------------------------------------------------
357      IF ( con_i ) THEN
358         CALL lim_column_sum (jpl,   v_i, vt_i_final)
359         fieldid = ' v_i : limitd_th '
360         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
361
362         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, et_i_final)
363         fieldid = ' e_i : limitd_th '
364         CALL lim_cons_check (et_i_init, et_i_final, 1.0e-3, fieldid) 
365
366         CALL lim_column_sum (jpl,   v_s, vt_s_final)
367         fieldid = ' v_s : limitd_th '
368         CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 
369
370         dummy_es(:,:,:) = e_s(:,:,1,:)
371         CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_final)
372         fieldid = ' e_s : limitd_th '
373         CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid) 
374      ENDIF
375
376
377   END SUBROUTINE lim_itd_th_rem
378
379
380   SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag )
381      !!------------------------------------------------------------------
382      !!                ***  ROUTINE lim_itd_fitline ***
383      !!
384      !! ** Purpose :   fit g(h) with a line using area, volume constraints
385      !!
386      !! ** Method  :   Fit g(h) with a line, satisfying area and volume constraints.
387      !!              To reduce roundoff errors caused by large values of g0 and g1,
388      !!              we actually compute g(eta), where eta = h - hL, and hL is the
389      !!              left boundary.
390      !!------------------------------------------------------------------
391      INTEGER                     , INTENT(in   ) ::   num_cat      ! category index
392      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   HbL, HbR     ! left and right category boundaries
393      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   hice         ! ice thickness
394      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   g0, g1       ! coefficients in linear equation for g(eta)
395      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   hL           ! min value of range over which g(h) > 0
396      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   hR           ! max value of range over which g(h) > 0
397      INTEGER , DIMENSION(jpi,jpj), INTENT(in   ) ::   zremap_flag  !
398      !
399      INTEGER  ::   ji,jj        ! horizontal indices
400      REAL(wp) ::   zh13         ! HbL + 1/3 * (HbR - HbL)
401      REAL(wp) ::   zh23         ! HbL + 2/3 * (HbR - HbL)
402      REAL(wp) ::   zdhr         ! 1 / (hR - hL)
403      REAL(wp) ::   zwk1, zwk2   ! temporary variables
404      !!------------------------------------------------------------------
405      !
406      DO jj = 1, jpj
407         DO ji = 1, jpi
408            !
409            IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10   &
410               &                        .AND. hice(ji,jj)        > 0._wp )  THEN
411
412               ! Initialize hL and hR
413               hL(ji,jj) = HbL(ji,jj)
414               hR(ji,jj) = HbR(ji,jj)
415
416               ! Change hL or hR if hice falls outside central third of range
417               zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) )
418               zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) )
419
420               IF    ( hice(ji,jj) < zh13 ) THEN   ;   hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj)
421               ELSEIF( hice(ji,jj) > zh23 ) THEN   ;   hL(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hR(ji,jj)
422               ENDIF
423
424               ! Compute coefficients of g(eta) = g0 + g1*eta
425               zdhr = 1._wp / (hR(ji,jj) - hL(ji,jj))
426               zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr
427               zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr
428               g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 )
429               g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 )
430               !
431            ELSE  ! remap_flag = .false. or a_i < epsi10
432               hL(ji,jj) = 0._wp
433               hR(ji,jj) = 0._wp
434               g0(ji,jj) = 0._wp
435               g1(ji,jj) = 0._wp
436            ENDIF
437            !
438         END DO
439      END DO
440      !
441   END SUBROUTINE lim_itd_fitline
442
443
444   SUBROUTINE lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
445      !!------------------------------------------------------------------
446      !!                ***  ROUTINE lim_itd_shiftice ***
447      !!
448      !! ** Purpose :   shift ice across category boundaries, conserving everything
449      !!              ( area, volume, energy, age*vol, and mass of salt )
450      !!
451      !! ** Method  :
452      !!------------------------------------------------------------------
453      INTEGER                           , INTENT(in   ) ::   klbnd    ! Start thickness category index point
454      INTEGER                           , INTENT(in   ) ::   kubnd    ! End point on which the  the computation is applied
455      INTEGER , DIMENSION(jpi,jpj,jpl-1), INTENT(in   ) ::   zdonor   ! donor category index
456      REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) ::   zdaice   ! ice area transferred across boundary
457      REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) ::   zdvice   ! ice volume transferred across boundary
458
459      INTEGER ::   ji, jj, jl, jl2, jl1, jk   ! dummy loop indices
460      INTEGER ::   ii, ij                     ! indices when changing from 2D-1D is done
461
462      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zaTsfn
463      REAL(wp), DIMENSION(jpi,jpj)   ::   zworka            ! temporary array used here
464
465      REAL(wp) ::   zdvsnow, zdesnow   ! snow volume and energy transferred
466      REAL(wp) ::   zdeice             ! ice energy transferred
467      REAL(wp) ::   zdsm_vice          ! ice salinity times volume transferred
468      REAL(wp) ::   zdo_aice           ! ice age times volume transferred
469      REAL(wp) ::   zdaTsf             ! aicen*Tsfcn transferred
470
471      INTEGER, DIMENSION((jpi+1)*(jpj+1)) ::   nind_i, nind_j   ! compressed indices for i/j directions
472
473      INTEGER  ::   nbrem             ! number of cells with ice to transfer
474      !!------------------------------------------------------------------
475
476
477      !----------------------------------------------------------------------------------------------
478      ! 1) Define a variable equal to a_i*T_su
479      !----------------------------------------------------------------------------------------------
480
481      DO jl = klbnd, kubnd
482         zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl)
483      END DO
484
485      !-------------------------------------------------------------------------------
486      ! 2) Transfer volume and energy between categories
487      !-------------------------------------------------------------------------------
488
489      DO jl = klbnd, kubnd - 1
490         nbrem = 0
491         DO jj = 1, jpj
492            DO ji = 1, jpi
493               IF (zdaice(ji,jj,jl) > 0.0 ) THEN ! daice(n) can be < puny
494                  nbrem = nbrem + 1
495                  nind_i(nbrem) = ji
496                  nind_j(nbrem) = jj
497               ENDIF
498            END DO
499         END DO
500
501         DO ji = 1, nbrem 
502            ii = nind_i(ji)
503            ij = nind_j(ji)
504
505            jl1 = zdonor(ii,ij,jl)
506            rswitch       = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi10 ) )
507            zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi10 ) * rswitch
508            IF( jl1 == jl) THEN   ;   jl2 = jl1+1
509            ELSE                  ;   jl2 = jl 
510            ENDIF
511
512            !--------------
513            ! Ice areas
514            !--------------
515            a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl)
516            a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl)
517
518            !--------------
519            ! Ice volumes
520            !--------------
521            v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl) 
522            v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl)
523
524            !--------------
525            ! Snow volumes
526            !--------------
527            zdvsnow        = v_s(ii,ij,jl1) * zworka(ii,ij)
528            v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow
529            v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow 
530
531            !--------------------
532            ! Snow heat content 
533            !--------------------
534            zdesnow            = e_s(ii,ij,1,jl1) * zworka(ii,ij)
535            e_s(ii,ij,1,jl1)   = e_s(ii,ij,1,jl1) - zdesnow
536            e_s(ii,ij,1,jl2)   = e_s(ii,ij,1,jl2) + zdesnow
537
538            !--------------
539            ! Ice age
540            !--------------
541            zdo_aice           = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl)
542            oa_i(ii,ij,jl1)    = oa_i(ii,ij,jl1) - zdo_aice
543            oa_i(ii,ij,jl2)    = oa_i(ii,ij,jl2) + zdo_aice
544
545            !--------------
546            ! Ice salinity
547            !--------------
548            zdsm_vice          = smv_i(ii,ij,jl1) * zworka(ii,ij)
549            smv_i(ii,ij,jl1)   = smv_i(ii,ij,jl1) - zdsm_vice
550            smv_i(ii,ij,jl2)   = smv_i(ii,ij,jl2) + zdsm_vice
551
552            !---------------------
553            ! Surface temperature
554            !---------------------
555            zdaTsf             = t_su(ii,ij,jl1) * zdaice(ii,ij,jl)
556            zaTsfn(ii,ij,jl1)  = zaTsfn(ii,ij,jl1) - zdaTsf
557            zaTsfn(ii,ij,jl2)  = zaTsfn(ii,ij,jl2) + zdaTsf 
558
559         END DO
560
561         !------------------
562         ! Ice heat content
563         !------------------
564
565         DO jk = 1, nlay_i
566            DO ji = 1, nbrem
567               ii = nind_i(ji)
568               ij = nind_j(ji)
569
570               jl1 = zdonor(ii,ij,jl)
571               IF (jl1 == jl) THEN
572                  jl2 = jl+1
573               ELSE             ! n1 = n+1
574                  jl2 = jl 
575               ENDIF
576
577               zdeice = e_i(ii,ij,jk,jl1) * zworka(ii,ij)
578               e_i(ii,ij,jk,jl1) =  e_i(ii,ij,jk,jl1) - zdeice
579               e_i(ii,ij,jk,jl2) =  e_i(ii,ij,jk,jl2) + zdeice 
580            END DO
581         END DO
582
583      END DO                   ! boundaries, 1 to ncat-1
584
585      !-----------------------------------------------------------------
586      ! Update ice thickness and temperature
587      !-----------------------------------------------------------------
588
589      DO jl = klbnd, kubnd
590         DO jj = 1, jpj
591            DO ji = 1, jpi 
592               IF ( a_i(ji,jj,jl) > epsi10 ) THEN
593                  ht_i(ji,jj,jl)  =  v_i   (ji,jj,jl) / a_i(ji,jj,jl) 
594                  t_su(ji,jj,jl)  =  zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 
595               ELSE
596                  ht_i(ji,jj,jl)  = 0._wp
597                  t_su(ji,jj,jl)  = rt0
598               ENDIF
599            END DO
600         END DO
601      END DO
602      !
603      !
604   END SUBROUTINE lim_itd_shiftice
605   
606
607   SUBROUTINE lim_itd_th_reb( klbnd, kubnd )
608      !!------------------------------------------------------------------
609      !!                ***  ROUTINE lim_itd_th_reb ***
610      !!
611      !! ** Purpose : rebin - rebins thicknesses into defined categories
612      !!
613      !! ** Method  :
614      !!------------------------------------------------------------------
615      INTEGER , INTENT (in) ::   klbnd   ! Start thickness category index point
616      INTEGER , INTENT (in) ::   kubnd   ! End point on which the  the computation is applied
617      !
618      INTEGER ::   ji,jj, jl   ! dummy loop indices
619      INTEGER ::   zshiftflag          ! = .true. if ice must be shifted
620      CHARACTER (len = 15) :: fieldid
621
622      INTEGER , DIMENSION(jpi,jpj,jpl) ::   zdonor           ! donor category index
623      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zdaice, zdvice   ! ice area and volume transferred
624
625      REAL(wp), DIMENSION(jpi,jpj) ::   vt_i_init, vt_i_final   ! ice volume summed over categories
626      REAL(wp), DIMENSION(jpi,jpj) ::   vt_s_init, vt_s_final   ! snow volume summed over categories
627      !!------------------------------------------------------------------
628     
629      !     
630      IF( con_i ) THEN                 ! conservation check
631         CALL lim_column_sum (jpl,   v_i, vt_i_init)
632         CALL lim_column_sum (jpl,   v_s, vt_s_init)
633      ENDIF
634
635      !
636      !------------------------------------------------------------------------------
637      ! 1) Compute ice thickness.
638      !------------------------------------------------------------------------------
639      DO jl = klbnd, kubnd
640         DO jj = 1, jpj
641            DO ji = 1, jpi 
642               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
643               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch
644            END DO
645         END DO
646      END DO
647
648      !------------------------------------------------------------------------------
649      ! 2) If a category thickness is not in bounds, shift the
650      ! entire area, volume, and energy to the neighboring category
651      !------------------------------------------------------------------------------
652      !-------------------------
653      ! Initialize shift arrays
654      !-------------------------
655      DO jl = klbnd, kubnd
656         zdonor(:,:,jl) = 0
657         zdaice(:,:,jl) = 0._wp
658         zdvice(:,:,jl) = 0._wp
659      END DO
660
661      !-------------------------
662      ! Move thin categories up
663      !-------------------------
664
665      DO jl = klbnd, kubnd - 1  ! loop over category boundaries
666
667         !---------------------------------------
668         ! identify thicknesses that are too big
669         !---------------------------------------
670         zshiftflag = 0
671
672         DO jj = 1, jpj 
673            DO ji = 1, jpi 
674               IF( a_i(ji,jj,jl) > epsi10 .AND. ht_i(ji,jj,jl) > hi_max(jl) ) THEN
675                  zshiftflag        = 1
676                  zdonor(ji,jj,jl)  = jl 
677                  ! begin TECLIM change
678                  !zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * 0.5_wp
679                  !zdvice(ji,jj,jl)  = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1)) * 0.5_wp
680                  ! end TECLIM change
681                  ! clem: how much of a_i you send in cat sup is somewhat arbitrary
682                  zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi20 ) / ht_i(ji,jj,jl) 
683                  zdvice(ji,jj,jl)  = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi20 )
684               ENDIF
685            END DO
686         END DO
687         IF(lk_mpp)   CALL mpp_max( zshiftflag )
688
689         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories
690            CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
691            ! Reset shift parameters
692            zdonor(:,:,jl) = 0
693            zdaice(:,:,jl) = 0._wp
694            zdvice(:,:,jl) = 0._wp
695         ENDIF
696         !
697      END DO
698
699      !----------------------------
700      ! Move thick categories down
701      !----------------------------
702
703      DO jl = kubnd - 1, 1, -1       ! loop over category boundaries
704
705         !-----------------------------------------
706         ! Identify thicknesses that are too small
707         !-----------------------------------------
708         zshiftflag = 0
709
710         DO jj = 1, jpj
711            DO ji = 1, jpi
712               IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN
713                  !
714                  zshiftflag = 1
715                  zdonor(ji,jj,jl) = jl + 1
716                  zdaice(ji,jj,jl) = a_i(ji,jj,jl+1) 
717                  zdvice(ji,jj,jl) = v_i(ji,jj,jl+1)
718               ENDIF
719            END DO
720         END DO
721
722         IF(lk_mpp)   CALL mpp_max( zshiftflag )
723         
724         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories
725            CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
726            ! Reset shift parameters
727            zdonor(:,:,jl) = 0
728            zdaice(:,:,jl) = 0._wp
729            zdvice(:,:,jl) = 0._wp
730         ENDIF
731
732      END DO
733
734      !------------------------------------------------------------------------------
735      ! 3) Conservation check
736      !------------------------------------------------------------------------------
737
738      IF( con_i ) THEN
739         CALL lim_column_sum (jpl,   v_i, vt_i_final)
740         fieldid = ' v_i : limitd_reb '
741         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
742
743         CALL lim_column_sum (jpl,   v_s, vt_s_final)
744         fieldid = ' v_s : limitd_reb '
745         CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 
746      ENDIF
747      !
748
749   END SUBROUTINE lim_itd_th_reb
750
751#else
752   !!----------------------------------------------------------------------
753   !!   Default option            Dummy module         NO LIM sea-ice model
754   !!----------------------------------------------------------------------
755CONTAINS
756   SUBROUTINE lim_itd_th_rem
757   END SUBROUTINE lim_itd_th_rem
758   SUBROUTINE lim_itd_fitline
759   END SUBROUTINE lim_itd_fitline
760   SUBROUTINE lim_itd_shiftice
761   END SUBROUTINE lim_itd_shiftice
762   SUBROUTINE lim_itd_th_reb
763   END SUBROUTINE lim_itd_th_reb
764#endif
765   !!======================================================================
766END MODULE limitd_th
Note: See TracBrowser for help on using the repository browser.