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/2018/dev_r9866_HPC_03_globcom/src/ICE – NEMO

source: NEMO/branches/2018/dev_r9866_HPC_03_globcom/src/ICE/icevar.F90 @ 10288

Last change on this file since 10288 was 10288, checked in by francesca, 5 years ago

reduce global communications, see #2010

  • Property svn:keywords set to Id
File size: 43.1 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      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      WHERE( pato_i(:,:)   < 0._wp )   pato_i(:,:)   = 0._wp
558      WHERE( poa_i (:,:,:) < 0._wp )   poa_i (:,:,:) = 0._wp
559      WHERE( pa_i  (:,:,:) < 0._wp )   pa_i  (:,:,:) = 0._wp
560      WHERE( pa_ip (:,:,:) < 0._wp )   pa_ip (:,:,:) = 0._wp
561      WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
562      !                                                        but it does not change conservation, so keep it this way is ok
563      !
564      DO jl = 1, jpl       !==  loop over the categories  ==!
565         !
566         !----------------------------------------
567         ! zap ice energy and send it to the ocean
568         !----------------------------------------
569         DO jk = 1, nlay_i
570            DO jj = 1 , jpj
571               DO ji = 1 , jpi
572                  IF( pe_i(ji,jj,jk,jl) < 0._wp ) THEN
573                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
574                     pe_i(ji,jj,jk,jl) = 0._wp
575                  ENDIF
576               END DO
577            END DO
578         END DO
579         !
580         DO jk = 1, nlay_s
581            DO jj = 1 , jpj
582               DO ji = 1 , jpi
583                  IF( pe_s(ji,jj,jk,jl) < 0._wp ) THEN
584                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
585                     pe_s(ji,jj,jk,jl) = 0._wp
586                  ENDIF
587               END DO
588            END DO
589         END DO
590         !
591         !-----------------------------------------------------
592         ! zap ice and snow volume, add water and salt to ocean
593         !-----------------------------------------------------
594         DO jj = 1 , jpj
595            DO ji = 1 , jpi
596              IF( pv_i(ji,jj,jl) < 0._wp ) THEN
597                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * r1_rdtice
598                  pv_i   (ji,jj,jl) = 0._wp
599               ENDIF
600               IF( pv_s(ji,jj,jl) < 0._wp ) THEN
601                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * r1_rdtice
602                  pv_s   (ji,jj,jl) = 0._wp
603               ENDIF
604               IF( psv_i(ji,jj,jl) < 0._wp ) THEN
605                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * r1_rdtice
606                  psv_i  (ji,jj,jl) = 0._wp
607               ENDIF
608            END DO
609         END DO
610         !
611      END DO 
612      !
613   END SUBROUTINE ice_var_zapneg
614
615   
616   SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
617      !!-------------------------------------------------------------------
618      !!                ***  ROUTINE ice_var_itd   ***
619      !!
620      !! ** Purpose :  converting 1-cat ice to multiple ice categories
621      !!
622      !!                  ice thickness distribution follows a gaussian law
623      !!               around the concentration of the most likely ice thickness
624      !!                           (similar as iceistate.F90)
625      !!
626      !! ** Method:   Iterative procedure
627      !!               
628      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
629      !!
630      !!               2) Check whether the distribution conserves area and volume, positivity and
631      !!                  category boundaries
632      !!             
633      !!               3) If not (input ice is too thin), the last category is empty and
634      !!                  the number of categories is reduced (jpl-1)
635      !!
636      !!               4) Iterate until ok (SUM(itest(:) = 4)
637      !!
638      !! ** Arguments : zhti: 1-cat ice thickness
639      !!                zhts: 1-cat snow depth
640      !!                zati: 1-cat ice concentration
641      !!
642      !! ** Output    : jpl-cat
643      !!
644      !!  (Example of application: BDY forcings when input are cell averaged) 
645      !!-------------------------------------------------------------------
646      INTEGER  :: ji, jk, jl             ! dummy loop indices
647      INTEGER  :: idim, i_fill, jl0 
648      REAL(wp) :: zarg, zV, zconv, zdh, zdv
649      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
650      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i ! output ice/snow variables
651      INTEGER , DIMENSION(4)                  ::   itest
652      !!-------------------------------------------------------------------
653      !
654      ! ----------------------------------------
655      ! distribution over the jpl ice categories
656      ! ----------------------------------------
657      ! a gaussian distribution for ice concentration is used
658      ! then we check whether the distribution fullfills
659      ! volume and area conservation, positivity and ice categories bounds
660      idim = SIZE( zhti , 1 )
661      zh_i(1:idim,1:jpl) = 0._wp
662      zh_s(1:idim,1:jpl) = 0._wp
663      za_i(1:idim,1:jpl) = 0._wp
664      !
665      DO ji = 1, idim
666         !
667         IF( zhti(ji) > 0._wp ) THEN
668            !
669            ! find which category (jl0) the input ice thickness falls into
670            jl0 = jpl
671            DO jl = 1, jpl
672               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
673                  jl0 = jl
674                  CYCLE
675               ENDIF
676            END DO
677            !
678            itest(:) = 0
679            i_fill   = jpl + 1                                            !------------------------------------
680            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
681               !                                                          !------------------------------------
682               i_fill = i_fill - 1
683               !
684               zh_i(ji,1:jpl) = 0._wp
685               za_i(ji,1:jpl) = 0._wp
686               itest(:)       = 0     
687               !
688               IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
689                  zh_i(ji,1) = zhti(ji)
690                  za_i (ji,1) = zati (ji)
691               ELSE                         !-- case ice is thicker: fill categories >1
692                  ! thickness
693                  DO jl = 1, i_fill - 1
694                     zh_i(ji,jl) = hi_mean(jl)
695                  END DO
696                  !
697                  ! concentration
698                  za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
699                  DO jl = 1, i_fill - 1
700                     IF ( jl /= jl0 ) THEN
701                        zarg        = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
702                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
703                     ENDIF
704                  END DO
705                  !
706                  ! last category
707                  za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
708                  zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
709                  zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
710                  !
711                  ! correction if concentration of upper cat is greater than lower cat
712                  !    (it should be a gaussian around jl0 but sometimes it is not)
713                  IF ( jl0 /= jpl ) THEN
714                     DO jl = jpl, jl0+1, -1
715                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
716                           zdv = zh_i(ji,jl) * za_i(ji,jl)
717                           zh_i(ji,jl    ) = 0._wp
718                           za_i (ji,jl    ) = 0._wp
719                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
720                        END IF
721                     END DO
722                  ENDIF
723                  !
724               ENDIF
725               !
726               ! Compatibility tests
727               zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) ) 
728               IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation
729               !
730               zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )
731               IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation
732               !
733               IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ?
734               !
735               itest(4) = 1
736               DO jl = 1, i_fill
737                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0                                ! Test 4: positivity of ice concentrations
738               END DO
739               !                                         !----------------------------
740            END DO                                       ! end iteration on categories
741            !                                            !----------------------------
742         ENDIF
743      END DO
744
745      ! Add Snow in each category where za_i is not 0
746      DO jl = 1, jpl
747         DO ji = 1, idim
748            IF( za_i(ji,jl) > 0._wp ) THEN
749               zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
750               ! In case snow load is in excess that would lead to transformation from snow to ice
751               ! Then, transfer the snow excess into the ice (different from icethd_dh)
752               zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 ) 
753               ! recompute h_i, h_s avoiding out of bounds values
754               zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
755               zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos )
756            ENDIF
757         END DO
758      END DO
759      !
760   END SUBROUTINE ice_var_itd
761
762
763   SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )
764      !!-------------------------------------------------------------------
765      !!                ***  ROUTINE ice_var_itd2   ***
766      !!
767      !! ** Purpose :  converting N-cat ice to jpl ice categories
768      !!
769      !!                  ice thickness distribution follows a gaussian law
770      !!               around the concentration of the most likely ice thickness
771      !!                           (similar as iceistate.F90)
772      !!
773      !! ** Method:   Iterative procedure
774      !!               
775      !!               1) Fill ice cat that correspond to input thicknesses
776      !!                  Find the lowest(jlmin) and highest(jlmax) cat that are filled
777      !!
778      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1
779      !!                   by removing 25% ice area from jlmin and jlmax (resp.)
780      !!             
781      !!               3) Expand the filling to the empty cat between jlmin and jlmax
782      !!                   by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
783      !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
784      !!
785      !! ** Arguments : zhti: N-cat ice thickness
786      !!                zhts: N-cat snow depth
787      !!                zati: N-cat ice concentration
788      !!
789      !! ** Output    : jpl-cat
790      !!
791      !!  (Example of application: BDY forcings when inputs have N-cat /= jpl) 
792      !!-------------------------------------------------------------------
793      INTEGER  ::   ji, jl, jl1, jl2             ! dummy loop indices
794      INTEGER  ::   idim, icat 
795      INTEGER, PARAMETER ::   ztrans = 0.25_wp
796      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
797      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
798      INTEGER , DIMENSION(:,:), ALLOCATABLE   ::   jlfil, jlfil2
799      INTEGER , DIMENSION(:)  , ALLOCATABLE   ::   jlmax, jlmin
800      !!-------------------------------------------------------------------
801      !
802      idim = SIZE( zhti, 1 )
803      icat = SIZE( zhti, 2 )
804      !
805      ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays
806      ALLOCATE( jlmin(idim), jlmax(idim) )
807
808      ! --- initialize output fields to 0 --- !
809      zh_i(1:idim,1:jpl) = 0._wp
810      zh_s(1:idim,1:jpl) = 0._wp
811      za_i(1:idim,1:jpl) = 0._wp
812      !
813      ! --- fill the categories --- !
814      !     find where cat-input = cat-output and fill cat-output fields 
815      jlmax(:) = 0
816      jlmin(:) = 999
817      jlfil(:,:) = 0
818      DO jl1 = 1, jpl
819         DO jl2 = 1, icat
820            DO ji = 1, idim
821               IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
822                  ! fill the right category
823                  zh_i(ji,jl1) = zhti(ji,jl2)
824                  zh_s(ji,jl1) = zhts(ji,jl2)
825                  za_i(ji,jl1) = zati(ji,jl2)
826                  ! record categories that are filled
827                  jlmax(ji) = MAX( jlmax(ji), jl1 )
828                  jlmin(ji) = MIN( jlmin(ji), jl1 )
829                  jlfil(ji,jl1) = jl1
830               ENDIF
831            END DO
832         END DO
833      END DO
834      !
835      ! --- fill the gaps between categories --- ! 
836      !     transfer from categories filled at the previous step to the empty ones in between
837      DO ji = 1, idim
838         jl1 = jlmin(ji)
839         jl2 = jlmax(ji)
840         IF( jl1 > 1 ) THEN
841            ! fill the lower cat (jl1-1)
842            za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
843            zh_i(ji,jl1-1) = hi_mean(jl1-1)
844            ! remove from cat jl1
845            za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
846         ENDIF
847         IF( jl2 < jpl ) THEN
848            ! fill the upper cat (jl2+1)
849            za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
850            zh_i(ji,jl2+1) = hi_mean(jl2+1)
851            ! remove from cat jl2
852            za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
853         ENDIF
854      END DO
855      !
856      jlfil2(:,:) = jlfil(:,:) 
857      ! fill categories from low to high
858      DO jl = 2, jpl-1
859         DO ji = 1, idim
860            IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
861               ! fill high
862               za_i(ji,jl) = ztrans * za_i(ji,jl-1)
863               zh_i(ji,jl) = hi_mean(jl)
864               jlfil(ji,jl) = jl
865               ! remove low
866               za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
867            ENDIF
868         END DO
869      END DO
870      !
871      ! fill categories from high to low
872      DO jl = jpl-1, 2, -1
873         DO ji = 1, idim
874            IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
875               ! fill low
876               za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
877               zh_i(ji,jl) = hi_mean(jl) 
878               jlfil2(ji,jl) = jl
879               ! remove high
880               za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
881            ENDIF
882         END DO
883      END DO
884      !
885      DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
886      DEALLOCATE( jlmin, jlmax )
887      !
888   END SUBROUTINE ice_var_itd2
889
890
891   SUBROUTINE ice_var_bv
892      !!-------------------------------------------------------------------
893      !!                ***  ROUTINE ice_var_bv ***
894      !!
895      !! ** Purpose :   computes mean brine volume (%) in sea ice
896      !!
897      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
898      !!
899      !! References : Vancoppenolle et al., JGR, 2007
900      !!-------------------------------------------------------------------
901      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
902      !!-------------------------------------------------------------------
903      !
904!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
905!!   instead of setting everything to zero as just below
906      bv_i (:,:,:) = 0._wp
907      DO jl = 1, jpl
908         DO jk = 1, nlay_i
909            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
910               bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
911            END WHERE
912         END DO
913      END DO
914      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
915      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
916      END WHERE
917      !
918   END SUBROUTINE ice_var_bv
919
920
921   SUBROUTINE ice_var_enthalpy
922      !!-------------------------------------------------------------------
923      !!                   ***  ROUTINE ice_var_enthalpy ***
924      !!                 
925      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
926      !!
927      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
928      !!-------------------------------------------------------------------
929      INTEGER  ::   ji, jk   ! dummy loop indices
930      REAL(wp) ::   ztmelts  ! local scalar
931      !!-------------------------------------------------------------------
932      !
933      DO jk = 1, nlay_i             ! Sea ice energy of melting
934         DO ji = 1, npti
935            ztmelts      = - rTmlt  * sz_i_1d(ji,jk)
936            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
937                                                                !   (sometimes zdf scheme produces abnormally high temperatures)   
938            e_i_1d(ji,jk) = rhoi * ( rcpi  * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
939               &                   + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
940               &                   - rcp   * ztmelts )
941         END DO
942      END DO
943      DO jk = 1, nlay_s             ! Snow energy of melting
944         DO ji = 1, npti
945            e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
946         END DO
947      END DO
948      !
949   END SUBROUTINE ice_var_enthalpy
950
951#else
952   !!----------------------------------------------------------------------
953   !!   Default option         Dummy module           NO SI3 sea-ice model
954   !!----------------------------------------------------------------------
955#endif
956
957   !!======================================================================
958END MODULE icevar
Note: See TracBrowser for help on using the repository browser.