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

source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icevar.F90 @ 10345

Last change on this file since 10345 was 10345, checked in by smasson, 5 years ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10344, see #2133

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