source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icevar.F90 @ 12720

Last change on this file since 12720 was 12720, checked in by clem, 8 months ago

implementation of ice pond lids (before debugging)

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