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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icevar.F90 in NEMO/branches/UKMO/dev_r10037_GPU/src/ICE – NEMO

source: NEMO/branches/UKMO/dev_r10037_GPU/src/ICE/icevar.F90 @ 11467

Last change on this file since 11467 was 11467, checked in by andmirek, 5 years ago

Ticket #2197 allocate arrays at the beggining of the run

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