source: NEMO/branches/UKMO/NEMO_4.0_GO8_package/src/ICE/icevar.F90 @ 11082

Last change on this file since 11082 was 11082, checked in by davestorkey, 18 months ago

UKMO/NEMO_4.0_GO8_package : update to be relative to 11081 of NEMO_4.0_mirror.

File size: 48.4 KB
Line 
1MODULE icevar
2   !!======================================================================
3   !!                       ***  MODULE icevar ***
4   !!   sea-ice:  series of functions to transform or compute ice variables
5   !!======================================================================
6   !! History :   -   !  2006-01  (M. Vancoppenolle) Original code
7   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
8   !!----------------------------------------------------------------------
9#if defined key_si3
10   !!----------------------------------------------------------------------
11   !!   'key_si3'                                       SI3 sea-ice model
12   !!----------------------------------------------------------------------
13   !!
14   !!                 There are three sets of variables
15   !!                 VGLO : global variables of the model
16   !!                        - v_i (jpi,jpj,jpl)
17   !!                        - v_s (jpi,jpj,jpl)
18   !!                        - a_i (jpi,jpj,jpl)
19   !!                        - t_s (jpi,jpj,jpl)
20   !!                        - e_i (jpi,jpj,nlay_i,jpl)
21   !!                        - e_s (jpi,jpj,nlay_s,jpl)
22   !!                        - sv_i(jpi,jpj,jpl)
23   !!                        - oa_i(jpi,jpj,jpl)
24   !!                 VEQV : equivalent variables sometimes used in the model
25   !!                        - h_i(jpi,jpj,jpl)
26   !!                        - h_s(jpi,jpj,jpl)
27   !!                        - t_i(jpi,jpj,nlay_i,jpl)
28   !!                        ...
29   !!                 VAGG : aggregate variables, averaged/summed over all
30   !!                        thickness categories
31   !!                        - vt_i(jpi,jpj)
32   !!                        - vt_s(jpi,jpj)
33   !!                        - at_i(jpi,jpj)
34   !!                        - et_s(jpi,jpj)  total snow heat content
35   !!                        - et_i(jpi,jpj)  total ice thermal content
36   !!                        - sm_i(jpi,jpj)  mean ice salinity
37   !!                        - tm_i(jpi,jpj)  mean ice temperature
38   !!                        - tm_s(jpi,jpj)  mean snw temperature
39   !!----------------------------------------------------------------------
40   !!   ice_var_agg       : integrate variables over layers and categories
41   !!   ice_var_glo2eqv   : transform from VGLO to VEQV
42   !!   ice_var_eqv2glo   : transform from VEQV to VGLO
43   !!   ice_var_salprof   : salinity profile in the ice
44   !!   ice_var_salprof1d : salinity profile in the ice 1D
45   !!   ice_var_zapsmall  : remove very small area and volume
46   !!   ice_var_zapneg    : remove negative ice fields
47   !!   ice_var_roundoff  : remove negative values arising from roundoff erros
48   !!   ice_var_itd       : convert 1-cat to jpl-cat
49   !!   ice_var_itd2      : convert N-cat to jpl-cat
50   !!   ice_var_bv        : brine volume
51   !!   ice_var_enthalpy  : compute ice and snow enthalpies from temperature
52   !!   ice_var_sshdyn    : compute equivalent ssh in lead
53   !!----------------------------------------------------------------------
54   USE dom_oce        ! ocean space and time domain
55   USE phycst         ! physical constants (ocean directory)
56   USE sbc_oce , ONLY : sss_m, ln_ice_embd, nn_fsbc
57   USE ice            ! sea-ice: variables
58   USE ice1D          ! sea-ice: thermodynamics variables
59   !
60   USE in_out_manager ! I/O manager
61   USE lib_mpp        ! MPP library
62   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
63
64   IMPLICIT NONE
65   PRIVATE
66
67   PUBLIC   ice_var_agg         
68   PUBLIC   ice_var_glo2eqv     
69   PUBLIC   ice_var_eqv2glo     
70   PUBLIC   ice_var_salprof     
71   PUBLIC   ice_var_salprof1d   
72   PUBLIC   ice_var_zapsmall
73   PUBLIC   ice_var_zapneg
74   PUBLIC   ice_var_roundoff
75   PUBLIC   ice_var_itd
76   PUBLIC   ice_var_itd2
77   PUBLIC   ice_var_bv           
78   PUBLIC   ice_var_enthalpy           
79   PUBLIC   ice_var_sshdyn
80
81   !!----------------------------------------------------------------------
82   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
83   !! $Id$
84   !! Software governed by the CeCILL license (see ./LICENSE)
85   !!----------------------------------------------------------------------
86CONTAINS
87
88   SUBROUTINE ice_var_agg( kn )
89      !!-------------------------------------------------------------------
90      !!                ***  ROUTINE ice_var_agg  ***
91      !!
92      !! ** Purpose :   aggregates ice-thickness-category variables to
93      !!              all-ice variables, i.e. it turns VGLO into VAGG
94      !!-------------------------------------------------------------------
95      INTEGER, INTENT( in ) ::   kn     ! =1 state variables only
96      !                                 ! >1 state variables + others
97      !
98      INTEGER ::   ji, jj, jk, jl   ! dummy loop indices
99      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z1_at_i, z1_vt_i, z1_vt_s
100      !!-------------------------------------------------------------------
101      !
102      !                                      ! integrated values
103      vt_i(:,:) =       SUM( v_i(:,:,:)           , dim=3 )
104      vt_s(:,:) =       SUM( v_s(:,:,:)           , dim=3 )
105      at_i(:,:) =       SUM( a_i(:,:,:)           , dim=3 )
106      et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 )
107      et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
108      !
109      at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds
110      vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
111      !
112      ato_i(:,:) = 1._wp - at_i(:,:)         ! open water fraction 
113
114      ! The following fields are calculated for diagnostics and outputs only
115      ! ==> Do not use them for other purposes
116      IF( kn > 1 ) THEN
117         !
118         ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) )
119         WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:)
120         ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp
121         END WHERE
122         WHERE( vt_i(:,:) > epsi20 )   ;   z1_vt_i(:,:) = 1._wp / vt_i(:,:)
123         ELSEWHERE                     ;   z1_vt_i(:,:) = 0._wp
124         END WHERE
125         WHERE( vt_s(:,:) > epsi20 )   ;   z1_vt_s(:,:) = 1._wp / vt_s(:,:)
126         ELSEWHERE                     ;   z1_vt_s(:,:) = 0._wp
127         END WHERE
128         !
129         !                          ! mean ice/snow thickness
130         hm_i(:,:) = vt_i(:,:) * z1_at_i(:,:)
131         hm_s(:,:) = vt_s(:,:) * z1_at_i(:,:)
132         !         
133         !                          ! mean temperature (K), salinity and age
134         tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
135         tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
136         om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(:,:)
137         sm_i (:,:) = SUM( sv_i(:,:,:)              , dim=3 ) * z1_vt_i(:,:)
138         !
139         tm_i(:,:) = 0._wp
140         tm_s(:,:) = 0._wp
141         DO jl = 1, jpl
142            DO jk = 1, nlay_i
143               tm_i(:,:) = tm_i(:,:) + r1_nlay_i * t_i (:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
144            END DO
145            DO jk = 1, nlay_s
146               tm_s(:,:) = tm_s(:,:) + r1_nlay_s * t_s (:,:,jk,jl) * v_s(:,:,jl) * z1_vt_s(:,:)
147            END DO
148         END DO
149         !
150         !                           ! put rt0 where there is no ice
151         WHERE( at_i(:,:)<=epsi20 )
152            tm_su(:,:) = rt0
153            tm_si(:,:) = rt0
154            tm_i (:,:) = rt0
155            tm_s (:,:) = rt0
156         END WHERE
157
158         DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s )
159      ENDIF
160      !
161   END SUBROUTINE ice_var_agg
162
163
164   SUBROUTINE ice_var_glo2eqv
165      !!-------------------------------------------------------------------
166      !!                ***  ROUTINE ice_var_glo2eqv ***
167      !!
168      !! ** Purpose :   computes equivalent variables as function of 
169      !!              global variables, i.e. it turns VGLO into VEQV
170      !!-------------------------------------------------------------------
171      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
172      REAL(wp) ::   ze_i             ! local scalars
173      REAL(wp) ::   ze_s, ztmelts, zbbb, zccc       !   -      -
174      REAL(wp) ::   zhmax, z1_zhmax                 !   -      -
175      REAL(wp) ::   zlay_i, zlay_s                  !   -      -
176      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i
177      !!-------------------------------------------------------------------
178
179!!gm Question 2:  It is possible to define existence of sea-ice in a common way between
180!!                ice area and ice volume ?
181!!                the idea is to be able to define one for all at the begining of this routine
182!!                a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
183
184      !---------------------------------------------------------------
185      ! Ice thickness, snow thickness, ice salinity, ice age and ponds
186      !---------------------------------------------------------------
187      !                                            !--- inverse of the ice area
188      WHERE( a_i(:,:,:) > epsi20 )   ;   z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
189      ELSEWHERE                      ;   z1_a_i(:,:,:) = 0._wp
190      END WHERE
191      !
192      WHERE( v_i(:,:,:) > epsi20 )   ;   z1_v_i(:,:,:) = 1._wp / v_i(:,:,:)
193      ELSEWHERE                      ;   z1_v_i(:,:,:) = 0._wp
194      END WHERE
195      !                                           !--- ice thickness
196      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
197
198      zhmax    =          hi_max(jpl)
199      z1_zhmax =  1._wp / hi_max(jpl)               
200      WHERE( h_i(:,:,jpl) > zhmax )   ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area
201         h_i  (:,:,jpl) = zhmax
202         a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 
203         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)
204      END WHERE
205      !                                           !--- snow thickness
206      h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:)
207      !                                           !--- ice age     
208      o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:)
209      !                                           !--- pond fraction and thickness     
210      a_ip_frac(:,:,:) = a_ip(:,:,:) * z1_a_i(:,:,:)
211      WHERE( a_ip_frac(:,:,:) > epsi20 )   ;   h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:)
212      ELSEWHERE                            ;   h_ip(:,:,:) = 0._wp
213      END WHERE
214      !
215      !                                           !---  salinity (with a minimum value imposed everywhere)     
216      IF( nn_icesal == 2 ) THEN
217         WHERE( v_i(:,:,:) > epsi20 )   ;   s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) * z1_v_i(:,:,:) ) )
218         ELSEWHERE                      ;   s_i(:,:,:) = rn_simin
219         END WHERE
220      ENDIF
221      CALL ice_var_salprof   ! salinity profile
222
223      !-------------------
224      ! Ice temperature   [K]   (with a minimum value (rt0 - 100.))
225      !-------------------
226      zlay_i   = REAL( nlay_i , wp )    ! number of layers
227      DO jl = 1, jpl
228         DO jk = 1, nlay_i
229            DO jj = 1, jpj
230               DO ji = 1, jpi
231                  IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area
232                     !
233                     ze_i             =   e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i             ! Energy of melting e(S,T) [J.m-3]
234                     ztmelts          = - sz_i(ji,jj,jk,jl) * rTmlt                                 ! Ice layer melt temperature [C]
235                     ! Conversion q(S,T) -> T (second order equation)
236                     zbbb             = ( rcp - rcpi ) * ztmelts + ze_i * r1_rhoi - rLfus
237                     zccc             = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts , 0._wp) )
238                     t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts
239                     !
240                  ELSE                                   !--- no ice
241                     t_i(ji,jj,jk,jl) = rt0
242                  ENDIF
243               END DO
244            END DO
245         END DO
246      END DO
247
248      !--------------------
249      ! Snow temperature   [K]   (with a minimum value (rt0 - 100.))
250      !--------------------
251      zlay_s = REAL( nlay_s , wp )
252      DO jk = 1, nlay_s
253         WHERE( v_s(:,:,:) > epsi20 )        !--- icy area
254            t_s(:,:,jk,:) = rt0 + MAX( -100._wp ,  &
255                 &                MIN( r1_rcpi * ( -r1_rhos * ( e_s(:,:,jk,:) / v_s(:,:,:) * zlay_s ) + rLfus ) , 0._wp ) )
256         ELSEWHERE                           !--- no ice
257            t_s(:,:,jk,:) = rt0
258         END WHERE
259      END DO
260      !
261      ! integrated values
262      vt_i (:,:) = SUM( v_i, dim=3 )
263      vt_s (:,:) = SUM( v_s, dim=3 )
264      at_i (:,:) = SUM( a_i, dim=3 )
265      !
266   END SUBROUTINE ice_var_glo2eqv
267
268
269   SUBROUTINE ice_var_eqv2glo
270      !!-------------------------------------------------------------------
271      !!                ***  ROUTINE ice_var_eqv2glo ***
272      !!
273      !! ** Purpose :   computes global variables as function of
274      !!              equivalent variables,  i.e. it turns VEQV into VGLO
275      !!-------------------------------------------------------------------
276      !
277      v_i (:,:,:) = h_i (:,:,:) * a_i (:,:,:)
278      v_s (:,:,:) = h_s (:,:,:) * a_i (:,:,:)
279      sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:)
280      v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
281      !
282   END SUBROUTINE ice_var_eqv2glo
283
284
285   SUBROUTINE ice_var_salprof
286      !!-------------------------------------------------------------------
287      !!                ***  ROUTINE ice_var_salprof ***
288      !!
289      !! ** Purpose :   computes salinity profile in function of bulk salinity     
290      !!
291      !! ** Method  : If bulk salinity greater than zsi1,
292      !!              the profile is assumed to be constant (S_inf)
293      !!              If bulk salinity lower than zsi0,
294      !!              the profile is linear with 0 at the surface (S_zero)
295      !!              If it is between zsi0 and zsi1, it is a
296      !!              alpha-weighted linear combination of s_inf and s_zero
297      !!
298      !! ** References : Vancoppenolle et al., 2007
299      !!-------------------------------------------------------------------
300      INTEGER  ::   ji, jj, jk, jl   ! dummy loop index
301      REAL(wp) ::   zsal, z1_dS
302      REAL(wp) ::   zargtemp , zs0, zs
303      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z_slope_s, zalpha    ! case 2 only
304      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
305      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
306      !!-------------------------------------------------------------------
307
308!!gm Question: Remove the option 3 ?  How many years since it last use ?
309
310      SELECT CASE ( nn_icesal )
311      !
312      !               !---------------------------------------!
313      CASE( 1 )       !  constant salinity in time and space  !
314         !            !---------------------------------------!
315         sz_i(:,:,:,:) = rn_icesal
316         s_i (:,:,:)   = rn_icesal
317         !
318         !            !---------------------------------------------!
319      CASE( 2 )       !  time varying salinity with linear profile  !
320         !            !---------------------------------------------!
321         !
322         ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
323         !
324         DO jl = 1, jpl
325            DO jk = 1, nlay_i
326               sz_i(:,:,jk,jl)  = s_i(:,:,jl)
327            END DO
328         END DO
329         !                                      ! Slope of the linear profile
330         WHERE( h_i(:,:,:) > epsi20 )   ;   z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:)
331         ELSEWHERE                      ;   z_slope_s(:,:,:) = 0._wp
332         END WHERE
333         !
334         z1_dS = 1._wp / ( zsi1 - zsi0 )
335         DO jl = 1, jpl
336            DO jj = 1, jpj
337               DO ji = 1, jpi
338                  zalpha(ji,jj,jl) = MAX(  0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp )  )
339                  !                             ! force a constant profile when SSS too low (Baltic Sea)
340                  IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) )   zalpha(ji,jj,jl) = 0._wp 
341               END DO
342            END DO
343         END DO
344         !
345         ! Computation of the profile
346         DO jl = 1, jpl
347            DO jk = 1, nlay_i
348               DO jj = 1, jpj
349                  DO ji = 1, jpi
350                     !                          ! linear profile with 0 surface value
351                     zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
352                     zs  = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl)     ! weighting the profile
353                     sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
354                  END DO
355               END DO
356            END DO
357         END DO
358         !
359         DEALLOCATE( z_slope_s , zalpha )
360         !
361         !            !-------------------------------------------!
362      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
363         !            !-------------------------------------------!                                   (mean = 2.30)
364         !
365         s_i(:,:,:) = 2.30_wp
366!!gm Remark: if we keep the case 3, then compute an store one for all time-step
367!!           a array  S_prof(1:nlay_i) containing the calculation and just do:
368!         DO jk = 1, nlay_i
369!            sz_i(:,:,jk,:) = S_prof(jk)
370!         END DO
371!!gm end
372         !
373         DO jl = 1, jpl
374            DO jk = 1, nlay_i
375               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
376               sz_i(:,:,jk,jl) =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
377            END DO
378         END DO
379         !
380      END SELECT
381      !
382   END SUBROUTINE ice_var_salprof
383
384
385   SUBROUTINE ice_var_salprof1d
386      !!-------------------------------------------------------------------
387      !!                  ***  ROUTINE ice_var_salprof1d  ***
388      !!
389      !! ** Purpose :   1d computation of the sea ice salinity profile
390      !!                Works with 1d vectors and is used by thermodynamic modules
391      !!-------------------------------------------------------------------
392      INTEGER  ::   ji, jk    ! dummy loop indices
393      REAL(wp) ::   zargtemp, zsal, z1_dS   ! local scalars
394      REAL(wp) ::   zs, zs0              !   -      -
395      !
396      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z_slope_s, zalpha   !
397      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
398      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
399      !!-------------------------------------------------------------------
400      !
401      SELECT CASE ( nn_icesal )
402      !
403      !               !---------------------------------------!
404      CASE( 1 )       !  constant salinity in time and space  !
405         !            !---------------------------------------!
406         sz_i_1d(1:npti,:) = rn_icesal
407         !
408         !            !---------------------------------------------!
409      CASE( 2 )       !  time varying salinity with linear profile  !
410         !            !---------------------------------------------!
411         !
412         ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
413         !
414         !                                      ! Slope of the linear profile
415         WHERE( h_i_1d(1:npti) > epsi20 )   ;   z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
416         ELSEWHERE                          ;   z_slope_s(1:npti) = 0._wp
417         END WHERE
418         
419         z1_dS = 1._wp / ( zsi1 - zsi0 )
420         DO ji = 1, npti
421            zalpha(ji) = MAX(  0._wp , MIN(  ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp  )  )
422            !                             ! force a constant profile when SSS too low (Baltic Sea)
423            IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) )   zalpha(ji) = 0._wp
424         END DO
425         !
426         ! Computation of the profile
427         DO jk = 1, nlay_i
428            DO ji = 1, npti
429               !                          ! linear profile with 0 surface value
430               zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
431               zs  = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
432               sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
433            END DO
434         END DO
435         !
436         DEALLOCATE( z_slope_s, zalpha )
437
438         !            !-------------------------------------------!
439      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
440         !            !-------------------------------------------!                                   (mean = 2.30)
441         !
442         s_i_1d(1:npti) = 2.30_wp
443         !
444!!gm cf remark in ice_var_salprof routine, CASE( 3 )
445         DO jk = 1, nlay_i
446            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
447            zsal =  1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
448            DO ji = 1, npti
449               sz_i_1d(ji,jk) = zsal
450            END DO
451         END DO
452         !
453      END SELECT
454      !
455   END SUBROUTINE ice_var_salprof1d
456
457
458   SUBROUTINE ice_var_zapsmall
459      !!-------------------------------------------------------------------
460      !!                   ***  ROUTINE ice_var_zapsmall ***
461      !!
462      !! ** Purpose :   Remove too small sea ice areas and correct fluxes
463      !!-------------------------------------------------------------------
464      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
465      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch
466      !!-------------------------------------------------------------------
467      !
468      DO jl = 1, jpl       !==  loop over the categories  ==!
469         !
470         WHERE( a_i(:,:,jl) > epsi10 )   ;   h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
471         ELSEWHERE                       ;   h_i(:,:,jl) = 0._wp
472         END WHERE
473         !
474         WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 )   ;   zswitch(:,:) = 0._wp
475         ELSEWHERE                                                                           ;   zswitch(:,:) = 1._wp
476         END WHERE
477         !
478         !-----------------------------------------------------------------
479         ! Zap ice energy and use ocean heat to melt ice
480         !-----------------------------------------------------------------
481         DO jk = 1, nlay_i
482            DO jj = 1 , jpj
483               DO ji = 1 , jpi
484                  ! update exchanges with ocean
485                  hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
486                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
487                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
488               END DO
489            END DO
490         END DO
491         !
492         DO jk = 1, nlay_s
493            DO jj = 1 , jpj
494               DO ji = 1 , jpi
495                  ! update exchanges with ocean
496                  hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
497                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj)
498                  t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
499               END DO
500            END DO
501         END DO
502         !
503         !-----------------------------------------------------------------
504         ! zap ice and snow volume, add water and salt to ocean
505         !-----------------------------------------------------------------
506         DO jj = 1 , jpj
507            DO ji = 1 , jpi
508               ! update exchanges with ocean
509               sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoi * r1_rdtice
510               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoi * r1_rdtice
511               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhos * r1_rdtice
512               !
513               a_i  (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
514               v_i  (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
515               v_s  (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
516               t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
517               oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
518               sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
519               !
520               h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
521               h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
522               !
523               a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
524               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
525               !
526            END DO
527         END DO
528         !
529      END DO 
530
531      ! to be sure that at_i is the sum of a_i(jl)
532      at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
533      vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
534
535      ! open water = 1 if at_i=0
536      WHERE( at_i(:,:) == 0._wp )   ato_i(:,:) = 1._wp
537      !
538   END SUBROUTINE ice_var_zapsmall
539
540
541   SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
542      !!-------------------------------------------------------------------
543      !!                   ***  ROUTINE ice_var_zapneg ***
544      !!
545      !! ** Purpose :   Remove negative sea ice fields and correct fluxes
546      !!-------------------------------------------------------------------
547      REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step
548      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
549      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
550      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
551      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_i      ! salt content
552      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
553      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
554      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
555      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
556      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
557      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
558      !
559      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
560      REAL(wp) ::   z1_dt
561      !!-------------------------------------------------------------------
562      !
563      z1_dt = 1._wp / pdt
564      !
565      DO jl = 1, jpl       !==  loop over the categories  ==!
566         !
567         ! make sure a_i=0 where v_i<=0
568         WHERE( pv_i(:,:,:) <= 0._wp )   pa_i(:,:,:) = 0._wp
569
570         !----------------------------------------
571         ! zap ice energy and send it to the ocean
572         !----------------------------------------
573         DO jk = 1, nlay_i
574            DO jj = 1 , jpj
575               DO ji = 1 , jpi
576                  IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
577                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
578                     pe_i(ji,jj,jk,jl) = 0._wp
579                  ENDIF
580               END DO
581            END DO
582         END DO
583         !
584         DO jk = 1, nlay_s
585            DO jj = 1 , jpj
586               DO ji = 1 , jpi
587                  IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
588                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
589                     pe_s(ji,jj,jk,jl) = 0._wp
590                  ENDIF
591               END DO
592            END DO
593         END DO
594         !
595         !-----------------------------------------------------
596         ! zap ice and snow volume, add water and salt to ocean
597         !-----------------------------------------------------
598         DO jj = 1 , jpj
599            DO ji = 1 , jpi
600               IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
601                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
602                  pv_i   (ji,jj,jl) = 0._wp
603               ENDIF
604               IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
605                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
606                  pv_s   (ji,jj,jl) = 0._wp
607               ENDIF
608               IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
609                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
610                  psv_i  (ji,jj,jl) = 0._wp
611               ENDIF
612            END DO
613         END DO
614         !
615      END DO 
616      !
617      WHERE( pato_i(:,:)   < 0._wp )   pato_i(:,:)   = 0._wp
618      WHERE( poa_i (:,:,:) < 0._wp )   poa_i (:,:,:) = 0._wp
619      WHERE( pa_i  (:,:,:) < 0._wp )   pa_i  (:,:,:) = 0._wp
620      WHERE( pa_ip (:,:,:) < 0._wp )   pa_ip (:,:,:) = 0._wp
621      WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
622      !                                                        but it does not change conservation, so keep it this way is ok
623      !
624   END SUBROUTINE ice_var_zapneg
625
626
627   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i )
628      !!-------------------------------------------------------------------
629      !!                   ***  ROUTINE ice_var_roundoff ***
630      !!
631      !! ** Purpose :   Remove negative sea ice values arising from roundoff errors
632      !!-------------------------------------------------------------------
633      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
634      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_i       ! ice volume
635      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_s       ! snw volume
636      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   psv_i      ! salt content
637      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   poa_i      ! age content
638      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
639      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
640      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
641      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
642      !!-------------------------------------------------------------------
643      !
644      WHERE( pa_i (1:npti,:)   < 0._wp .AND. pa_i (1:npti,:)   > -epsi10 )   pa_i (1:npti,:)   = 0._wp   !  a_i must be >= 0
645      WHERE( pv_i (1:npti,:)   < 0._wp .AND. pv_i (1:npti,:)   > -epsi10 )   pv_i (1:npti,:)   = 0._wp   !  v_i must be >= 0
646      WHERE( pv_s (1:npti,:)   < 0._wp .AND. pv_s (1:npti,:)   > -epsi10 )   pv_s (1:npti,:)   = 0._wp   !  v_s must be >= 0
647      WHERE( psv_i(1:npti,:)   < 0._wp .AND. psv_i(1:npti,:)   > -epsi10 )   psv_i(1:npti,:)   = 0._wp   ! sv_i must be >= 0
648      WHERE( poa_i(1:npti,:)   < 0._wp .AND. poa_i(1:npti,:)   > -epsi10 )   poa_i(1:npti,:)   = 0._wp   ! oa_i must be >= 0
649      WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0
650      WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0
651      IF ( ln_pnd_H12 ) THEN
652         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0
653         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0
654      ENDIF
655      !
656   END SUBROUTINE ice_var_roundoff
657   
658   
659   SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
660      !!-------------------------------------------------------------------
661      !!                ***  ROUTINE ice_var_itd   ***
662      !!
663      !! ** Purpose :  converting 1-cat ice to multiple ice categories
664      !!
665      !!                  ice thickness distribution follows a gaussian law
666      !!               around the concentration of the most likely ice thickness
667      !!                           (similar as iceistate.F90)
668      !!
669      !! ** Method:   Iterative procedure
670      !!               
671      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
672      !!
673      !!               2) Check whether the distribution conserves area and volume, positivity and
674      !!                  category boundaries
675      !!             
676      !!               3) If not (input ice is too thin), the last category is empty and
677      !!                  the number of categories is reduced (jpl-1)
678      !!
679      !!               4) Iterate until ok (SUM(itest(:) = 4)
680      !!
681      !! ** Arguments : zhti: 1-cat ice thickness
682      !!                zhts: 1-cat snow depth
683      !!                zati: 1-cat ice concentration
684      !!
685      !! ** Output    : jpl-cat
686      !!
687      !!  (Example of application: BDY forcings when input are cell averaged) 
688      !!-------------------------------------------------------------------
689      INTEGER  :: ji, jk, jl             ! dummy loop indices
690      INTEGER  :: idim, i_fill, jl0 
691      REAL(wp) :: zarg, zV, zconv, zdh, zdv
692      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input  ice/snow variables
693      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
694      INTEGER , DIMENSION(4)                  ::   itest
695      !!-------------------------------------------------------------------
696      !
697      ! ----------------------------------------
698      ! distribution over the jpl ice categories
699      ! ----------------------------------------
700      ! a gaussian distribution for ice concentration is used
701      ! then we check whether the distribution fullfills
702      ! volume and area conservation, positivity and ice categories bounds
703      idim = SIZE( zhti , 1 )
704      zh_i(1:idim,1:jpl) = 0._wp
705      zh_s(1:idim,1:jpl) = 0._wp
706      za_i(1:idim,1:jpl) = 0._wp
707      !
708      DO ji = 1, idim
709         !
710         IF( zhti(ji) > 0._wp ) THEN
711            !
712            ! find which category (jl0) the input ice thickness falls into
713            jl0 = jpl
714            DO jl = 1, jpl
715               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
716                  jl0 = jl
717                  CYCLE
718               ENDIF
719            END DO
720            !
721            itest(:) = 0
722            i_fill   = jpl + 1                                            !------------------------------------
723            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
724               !                                                          !------------------------------------
725               i_fill = i_fill - 1
726               !
727               zh_i(ji,1:jpl) = 0._wp
728               za_i(ji,1:jpl) = 0._wp
729               itest(:)       = 0     
730               !
731               IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
732                  zh_i(ji,1) = zhti(ji)
733                  za_i (ji,1) = zati (ji)
734               ELSE                         !-- case ice is thicker: fill categories >1
735                  ! thickness
736                  DO jl = 1, i_fill - 1
737                     zh_i(ji,jl) = hi_mean(jl)
738                  END DO
739                  !
740                  ! concentration
741                  za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
742                  DO jl = 1, i_fill - 1
743                     IF ( jl /= jl0 ) THEN
744                        zarg        = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
745                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
746                     ENDIF
747                  END DO
748                  !
749                  ! last category
750                  za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
751                  zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
752                  zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
753                  !
754                  ! correction if concentration of upper cat is greater than lower cat
755                  !    (it should be a gaussian around jl0 but sometimes it is not)
756                  IF ( jl0 /= jpl ) THEN
757                     DO jl = jpl, jl0+1, -1
758                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
759                           zdv = zh_i(ji,jl) * za_i(ji,jl)
760                           zh_i(ji,jl    ) = 0._wp
761                           za_i (ji,jl    ) = 0._wp
762                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
763                        END IF
764                     END DO
765                  ENDIF
766                  !
767               ENDIF
768               !
769               ! Compatibility tests
770               zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) ) 
771               IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation
772               !
773               zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )
774               IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation
775               !
776               IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ?
777               !
778               itest(4) = 1
779               DO jl = 1, i_fill
780                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0                                ! Test 4: positivity of ice concentrations
781               END DO
782               !                                         !----------------------------
783            END DO                                       ! end iteration on categories
784            !                                            !----------------------------
785         ENDIF
786      END DO
787
788      ! Add Snow in each category where za_i is not 0
789      DO jl = 1, jpl
790         DO ji = 1, idim
791            IF( za_i(ji,jl) > 0._wp ) THEN
792               zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
793               ! In case snow load is in excess that would lead to transformation from snow to ice
794               ! Then, transfer the snow excess into the ice (different from icethd_dh)
795               zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 ) 
796               ! recompute h_i, h_s avoiding out of bounds values
797               zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
798               zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos )
799            ENDIF
800         END DO
801      END DO
802      !
803   END SUBROUTINE ice_var_itd
804
805
806   SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )
807      !!-------------------------------------------------------------------
808      !!                ***  ROUTINE ice_var_itd2   ***
809      !!
810      !! ** Purpose :  converting N-cat ice to jpl ice categories
811      !!
812      !!                  ice thickness distribution follows a gaussian law
813      !!               around the concentration of the most likely ice thickness
814      !!                           (similar as iceistate.F90)
815      !!
816      !! ** Method:   Iterative procedure
817      !!               
818      !!               1) Fill ice cat that correspond to input thicknesses
819      !!                  Find the lowest(jlmin) and highest(jlmax) cat that are filled
820      !!
821      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1
822       !!                   by removing 25% ice area from jlmin and jlmax (resp.)
823      !!             
824      !!               3) Expand the filling to the empty cat between jlmin and jlmax
825      !!                   by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
826      !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
827      !!
828      !! ** Arguments : zhti: N-cat ice thickness
829      !!                zhts: N-cat snow depth
830      !!                zati: N-cat ice concentration
831      !!
832      !! ** Output    : jpl-cat
833      !!
834      !!  (Example of application: BDY forcings when inputs have N-cat /= jpl) 
835      !!-------------------------------------------------------------------
836      INTEGER  ::   ji, jl, jl1, jl2             ! dummy loop indices
837      INTEGER  ::   idim, icat 
838      REAL(wp), PARAMETER ::   ztrans = 0.25_wp
839      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
840      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
841      INTEGER , DIMENSION(:,:), ALLOCATABLE   ::   jlfil, jlfil2
842      INTEGER , DIMENSION(:)  , ALLOCATABLE   ::   jlmax, jlmin
843      !!-------------------------------------------------------------------
844      !
845      idim = SIZE( zhti, 1 )
846      icat = SIZE( zhti, 2 )
847      !
848      ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays
849      ALLOCATE( jlmin(idim), jlmax(idim) )
850
851      ! --- initialize output fields to 0 --- !
852      zh_i(1:idim,1:jpl) = 0._wp
853      zh_s(1:idim,1:jpl) = 0._wp
854      za_i(1:idim,1:jpl) = 0._wp
855      !
856      ! --- fill the categories --- !
857      !     find where cat-input = cat-output and fill cat-output fields 
858      jlmax(:) = 0
859      jlmin(:) = 999
860      jlfil(:,:) = 0
861      DO jl1 = 1, jpl
862         DO jl2 = 1, icat
863            DO ji = 1, idim
864               IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
865                  ! fill the right category
866                  zh_i(ji,jl1) = zhti(ji,jl2)
867                  zh_s(ji,jl1) = zhts(ji,jl2)
868                  za_i(ji,jl1) = zati(ji,jl2)
869                  ! record categories that are filled
870                  jlmax(ji) = MAX( jlmax(ji), jl1 )
871                  jlmin(ji) = MIN( jlmin(ji), jl1 )
872                  jlfil(ji,jl1) = jl1
873               ENDIF
874            END DO
875         END DO
876      END DO
877      !
878      ! --- fill the gaps between categories --- ! 
879      !     transfer from categories filled at the previous step to the empty ones in between
880      DO ji = 1, idim
881         jl1 = jlmin(ji)
882         jl2 = jlmax(ji)
883         IF( jl1 > 1 ) THEN
884            ! fill the lower cat (jl1-1)
885            za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
886            zh_i(ji,jl1-1) = hi_mean(jl1-1)
887            ! remove from cat jl1
888            za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
889         ENDIF
890         IF( jl2 < jpl ) THEN
891            ! fill the upper cat (jl2+1)
892            za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
893            zh_i(ji,jl2+1) = hi_mean(jl2+1)
894            ! remove from cat jl2
895            za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
896         ENDIF
897      END DO
898      !
899      jlfil2(:,:) = jlfil(:,:) 
900      ! fill categories from low to high
901      DO jl = 2, jpl-1
902         DO ji = 1, idim
903            IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
904               ! fill high
905               za_i(ji,jl) = ztrans * za_i(ji,jl-1)
906               zh_i(ji,jl) = hi_mean(jl)
907               jlfil(ji,jl) = jl
908               ! remove low
909               za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
910            ENDIF
911         END DO
912      END DO
913      !
914      ! fill categories from high to low
915      DO jl = jpl-1, 2, -1
916         DO ji = 1, idim
917            IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
918               ! fill low
919               za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
920               zh_i(ji,jl) = hi_mean(jl) 
921               jlfil2(ji,jl) = jl
922               ! remove high
923               za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
924            ENDIF
925         END DO
926      END DO
927      !
928      DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
929      DEALLOCATE( jlmin, jlmax )
930      !
931   END SUBROUTINE ice_var_itd2
932
933
934   SUBROUTINE ice_var_bv
935      !!-------------------------------------------------------------------
936      !!                ***  ROUTINE ice_var_bv ***
937      !!
938      !! ** Purpose :   computes mean brine volume (%) in sea ice
939      !!
940      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
941      !!
942      !! References : Vancoppenolle et al., JGR, 2007
943      !!-------------------------------------------------------------------
944      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
945      !!-------------------------------------------------------------------
946      !
947!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
948!!   instead of setting everything to zero as just below
949      bv_i (:,:,:) = 0._wp
950      DO jl = 1, jpl
951         DO jk = 1, nlay_i
952            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
953               bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
954            END WHERE
955         END DO
956      END DO
957      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
958      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
959      END WHERE
960      !
961   END SUBROUTINE ice_var_bv
962
963
964   SUBROUTINE ice_var_enthalpy
965      !!-------------------------------------------------------------------
966      !!                   ***  ROUTINE ice_var_enthalpy ***
967      !!                 
968      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
969      !!
970      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
971      !!-------------------------------------------------------------------
972      INTEGER  ::   ji, jk   ! dummy loop indices
973      REAL(wp) ::   ztmelts  ! local scalar
974      !!-------------------------------------------------------------------
975      !
976      DO jk = 1, nlay_i             ! Sea ice energy of melting
977         DO ji = 1, npti
978            ztmelts      = - rTmlt  * sz_i_1d(ji,jk)
979            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point => likely conservation issue
980                                                                !   (sometimes zdf scheme produces abnormally high temperatures)   
981            e_i_1d(ji,jk) = rhoi * ( rcpi  * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
982               &                   + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
983               &                   - rcp   * ztmelts )
984         END DO
985      END DO
986      DO jk = 1, nlay_s             ! Snow energy of melting
987         DO ji = 1, npti
988            e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
989         END DO
990      END DO
991      !
992   END SUBROUTINE ice_var_enthalpy
993
994   FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
995      !!---------------------------------------------------------------------
996      !!                   ***  ROUTINE ice_var_sshdyn  ***
997      !!                     
998      !! ** Purpose :  compute the equivalent ssh in lead when sea ice is embedded
999      !!
1000      !! ** Method  :  ssh_lead = ssh + (Mice + Msnow) / rau0
1001      !!
1002      !! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
1003      !!                Sea ice-ocean coupling using a rescaled vertical coordinate z*,
1004      !!                Ocean Modelling, Volume 24, Issues 1-2, 2008
1005      !!----------------------------------------------------------------------
1006      !
1007      ! input
1008      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh            !: ssh [m]
1009      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass    !: mass of snow and ice at current  ice time step [Kg/m2]
1010      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b  !: mass of snow and ice at previous ice time step [Kg/m2]
1011      !
1012      ! output
1013      REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn  ! equivalent ssh in lead [m]
1014      !
1015      ! temporary
1016      REAL(wp) :: zintn, zintb                     ! time interpolation weights []
1017      REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload  ! snow and ice load [m]
1018      !
1019      ! compute ice load used to define the equivalent ssh in lead
1020      IF( ln_ice_embd ) THEN
1021         !                                           
1022         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
1023         !                                               = (1/nn_fsbc)^2 * {SUM[n]        , n=0,nn_fsbc-1}
1024         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
1025         !
1026         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   *    {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
1027         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
1028         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
1029         !
1030         zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0
1031         !
1032      ELSE
1033         zsnwiceload(:,:) = 0.0_wp
1034      ENDIF
1035      ! compute equivalent ssh in lead
1036      ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:)
1037      !
1038   END FUNCTION ice_var_sshdyn
1039
1040
1041#else
1042   !!----------------------------------------------------------------------
1043   !!   Default option         Dummy module           NO SI3 sea-ice model
1044   !!----------------------------------------------------------------------
1045#endif
1046
1047   !!======================================================================
1048END MODULE icevar
Note: See TracBrowser for help on using the repository browser.