source: NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icevar.F90 @ 10411

Last change on this file since 10411 was 10411, checked in by clem, 23 months ago

make all the sette tests successful

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