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

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

fix previous revision. Make sure a_ip and a_ip_eff are correct

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