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_th.F90 in branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 8179

Last change on this file since 8179 was 8061, checked in by vancop, 7 years ago

Quick commit on empirical melt ponds

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