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

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

source: NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icevar.F90 @ 12937

Last change on this file since 12937 was 12937, checked in by dancopsey, 3 years ago

Merge in Clem's branch. It was originally here:

svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch

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 ) 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         IF( ln_pnd_lids ) THEN
686            WHERE( pv_il(1:npti,:) < 0._wp .AND. pv_il(1:npti,:) > -epsi10 ) pv_il(1:npti,:)   = 0._wp   ! v_il must be >= 0
687         ENDIF
688      ENDIF
689      !
690   END SUBROUTINE ice_var_roundoff
691   
692
693   SUBROUTINE ice_var_bv
694      !!-------------------------------------------------------------------
695      !!                ***  ROUTINE ice_var_bv ***
696      !!
697      !! ** Purpose :   computes mean brine volume (%) in sea ice
698      !!
699      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
700      !!
701      !! References : Vancoppenolle et al., JGR, 2007
702      !!-------------------------------------------------------------------
703      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
704      !!-------------------------------------------------------------------
705      !
706!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
707!!   instead of setting everything to zero as just below
708      bv_i (:,:,:) = 0._wp
709      DO jl = 1, jpl
710         DO jk = 1, nlay_i
711            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
712               bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
713            END WHERE
714         END DO
715      END DO
716      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
717      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
718      END WHERE
719      !
720   END SUBROUTINE ice_var_bv
721
722
723   SUBROUTINE ice_var_enthalpy
724      !!-------------------------------------------------------------------
725      !!                   ***  ROUTINE ice_var_enthalpy ***
726      !!                 
727      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
728      !!
729      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
730      !!-------------------------------------------------------------------
731      INTEGER  ::   ji, jk   ! dummy loop indices
732      REAL(wp) ::   ztmelts  ! local scalar
733      !!-------------------------------------------------------------------
734      !
735      DO jk = 1, nlay_i             ! Sea ice energy of melting
736         DO ji = 1, npti
737            ztmelts      = - rTmlt  * sz_i_1d(ji,jk)
738            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
739                                                                !   (sometimes zdf scheme produces abnormally high temperatures)   
740            e_i_1d(ji,jk) = rhoi * ( rcpi  * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
741               &                   + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
742               &                   - rcp   * ztmelts )
743         END DO
744      END DO
745      DO jk = 1, nlay_s             ! Snow energy of melting
746         DO ji = 1, npti
747            e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
748         END DO
749      END DO
750      !
751   END SUBROUTINE ice_var_enthalpy
752
753   
754   FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
755      !!---------------------------------------------------------------------
756      !!                   ***  ROUTINE ice_var_sshdyn  ***
757      !!                     
758      !! ** Purpose :  compute the equivalent ssh in lead when sea ice is embedded
759      !!
760      !! ** Method  :  ssh_lead = ssh + (Mice + Msnow) / rau0
761      !!
762      !! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
763      !!                Sea ice-ocean coupling using a rescaled vertical coordinate z*,
764      !!                Ocean Modelling, Volume 24, Issues 1-2, 2008
765      !!----------------------------------------------------------------------
766      !
767      ! input
768      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh            !: ssh [m]
769      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass    !: mass of snow and ice at current  ice time step [Kg/m2]
770      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b  !: mass of snow and ice at previous ice time step [Kg/m2]
771      !
772      ! output
773      REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn  ! equivalent ssh in lead [m]
774      !
775      ! temporary
776      REAL(wp) :: zintn, zintb                     ! time interpolation weights []
777      REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload  ! snow and ice load [m]
778      !
779      ! compute ice load used to define the equivalent ssh in lead
780      IF( ln_ice_embd ) THEN
781         !                                           
782         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
783         !                                               = (1/nn_fsbc)^2 * {SUM[n]        , n=0,nn_fsbc-1}
784         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
785         !
786         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   *    {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
787         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
788         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
789         !
790         zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0
791         !
792      ELSE
793         zsnwiceload(:,:) = 0.0_wp
794      ENDIF
795      ! compute equivalent ssh in lead
796      ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:)
797      !
798   END FUNCTION ice_var_sshdyn
799
800   
801   !!-------------------------------------------------------------------
802   !!                ***  INTERFACE ice_var_itd   ***
803   !!
804   !! ** Purpose :  converting N-cat ice to jpl ice categories
805   !!-------------------------------------------------------------------
806   SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,                             ph_i, ph_s, pa_i, &
807      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
808      !!-------------------------------------------------------------------
809      !! ** Purpose :  converting 1-cat ice to 1 ice category
810      !!-------------------------------------------------------------------
811      REAL(wp), DIMENSION(:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
812      REAL(wp), DIMENSION(:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
813      REAL(wp), DIMENSION(:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds
814      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
815      !!-------------------------------------------------------------------
816      ! == thickness and concentration == !
817      ph_i(:) = phti(:)
818      ph_s(:) = phts(:)
819      pa_i(:) = pati(:)
820      !
821      ! == temperature and salinity and ponds == !
822      pt_i (:) = ptmi (:)
823      pt_s (:) = ptms (:)
824      pt_su(:) = ptmsu(:)
825      ps_i (:) = psmi (:)
826      pa_ip(:) = patip(:)
827      ph_ip(:) = phtip(:)
828      ph_il(:) = phtil(:)
829     
830   END SUBROUTINE ice_var_itd_1c1c
831
832   SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,                             ph_i, ph_s, pa_i, &
833      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
834      !!-------------------------------------------------------------------
835      !! ** Purpose :  converting N-cat ice to 1 ice category
836      !!-------------------------------------------------------------------
837      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
838      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
839      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds
840      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
841      !
842      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z1_ai, z1_vi, z1_vs
843      !
844      INTEGER ::   idim 
845      !!-------------------------------------------------------------------
846      !
847      idim = SIZE( phti, 1 )
848      !
849      ! == thickness and concentration == !
850      ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim) )
851      !
852      pa_i(:) = SUM( pati(:,:), dim=2 )
853
854      WHERE( ( pa_i(:) ) /= 0._wp )   ;   z1_ai(:) = 1._wp / pa_i(:)
855      ELSEWHERE                       ;   z1_ai(:) = 0._wp
856      END WHERE
857
858      ph_i(:) = SUM( phti(:,:) * pati(:,:), dim=2 ) * z1_ai(:)
859      ph_s(:) = SUM( phts(:,:) * pati(:,:), dim=2 ) * z1_ai(:)
860      !
861      ! == temperature and salinity == !
862      WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp )   ;   z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) )
863      ELSEWHERE                                 ;   z1_vi(:) = 0._wp
864      END WHERE
865      WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp )   ;   z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) )
866      ELSEWHERE                                 ;   z1_vs(:) = 0._wp
867      END WHERE
868      pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
869      pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
870      pt_su(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:)
871      ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
872
873      ! == ponds == !
874      pa_ip(:) = SUM( patip(:,:), dim=2 )
875      WHERE( pa_ip(:) /= 0._wp )
876         ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
877         ph_il(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
878      ELSEWHERE
879         ph_ip(:) = 0._wp
880         ph_il(:) = 0._wp
881      END WHERE
882      !
883      DEALLOCATE( z1_ai, z1_vi, z1_vs )
884      !
885   END SUBROUTINE ice_var_itd_Nc1c
886   
887   SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,                             ph_i, ph_s, pa_i, &
888      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
889      !!-------------------------------------------------------------------
890      !!
891      !! ** Purpose :  converting 1-cat ice to jpl ice categories
892      !!
893      !!
894      !! ** Method:   ice thickness distribution follows a gamma function from Abraham et al. (2015)
895      !!              it has the property of conserving total concentration and volume
896      !!             
897      !!
898      !! ** Arguments : phti: 1-cat ice thickness
899      !!                phts: 1-cat snow depth
900      !!                pati: 1-cat ice concentration
901      !!
902      !! ** Output    : jpl-cat
903      !!
904      !!  Abraham, C., Steiner, N., Monahan, A. and Michel, C., 2015.
905      !!               Effects of subgrid‐scale snow thickness variability on radiative transfer in sea ice.
906      !!               Journal of Geophysical Research: Oceans, 120(8), pp.5597-5614
907      !!-------------------------------------------------------------------
908      REAL(wp), DIMENSION(:),   INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
909      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
910      REAL(wp), DIMENSION(:)  , INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds
911      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
912      !
913      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zfra, z1_hti
914      INTEGER  ::   ji, jk, jl
915      INTEGER  ::   idim
916      REAL(wp) ::   zv, zdh
917      !!-------------------------------------------------------------------
918      !
919      idim = SIZE( phti , 1 )
920      !
921      ph_i(1:idim,1:jpl) = 0._wp
922      ph_s(1:idim,1:jpl) = 0._wp
923      pa_i(1:idim,1:jpl) = 0._wp
924      !
925      ALLOCATE( z1_hti(idim) )
926      WHERE( phti(:) /= 0._wp )   ;   z1_hti(:) = 1._wp / phti(:)
927      ELSEWHERE                   ;   z1_hti(:) = 0._wp
928      END WHERE
929      !
930      ! == thickness and concentration == !
931      ! for categories 1:jpl-1, integrate the gamma function from hi_max(jl-1) to hi_max(jl)
932      DO jl = 1, jpl-1
933         DO ji = 1, idim
934            !
935            IF( phti(ji) > 0._wp ) THEN
936               ! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl)
937               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) ) &
938                  &                                   - ( phti(ji) + 2.*hi_max(jl  ) ) * EXP( -2.*hi_max(jl  )*z1_hti(ji) ) )
939               !
940               ! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl)
941               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) ) &
942                  &                            * EXP( -2.*hi_max(jl-1)*z1_hti(ji) ) &
943                  &                          - ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jl) + 2.*hi_max(jl)*hi_max(jl) ) &
944                  &                            * EXP(-2.*hi_max(jl)*z1_hti(ji)) )
945               ! thickness
946               IF( pa_i(ji,jl) > epsi06 ) THEN
947                  ph_i(ji,jl) = zv / pa_i(ji,jl)
948               ELSE
949                  ph_i(ji,jl) = 0.
950                  pa_i(ji,jl) = 0.
951               ENDIF
952            ENDIF
953            !
954         ENDDO
955      ENDDO
956      !
957      ! for the last category (jpl), integrate the gamma function from hi_max(jpl-1) to infinity
958      DO ji = 1, idim
959         !
960         IF( phti(ji) > 0._wp ) THEN
961            ! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jpl-1) to infinity
962            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) )
963
964            ! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jpl-1) to infinity
965            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) ) &
966               &                         * EXP( -2.*hi_max(jpl-1)*z1_hti(ji) )
967            ! thickness
968            IF( pa_i(ji,jpl) > epsi06 ) THEN
969               ph_i(ji,jpl) = zv / pa_i(ji,jpl)
970            else
971               ph_i(ji,jpl) = 0.
972               pa_i(ji,jpl) = 0.
973            ENDIF
974         ENDIF
975         !
976      ENDDO
977      !
978      ! Add Snow in each category where pa_i is not 0
979      DO jl = 1, jpl
980         DO ji = 1, idim
981            IF( pa_i(ji,jl) > 0._wp ) THEN
982               ph_s(ji,jl) = ph_i(ji,jl) * phts(ji) * z1_hti(ji)
983               ! In case snow load is in excess that would lead to transformation from snow to ice
984               ! Then, transfer the snow excess into the ice (different from icethd_dh)
985               zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rau0 ) * ph_i(ji,jl) ) * r1_rau0 ) 
986               ! recompute h_i, h_s avoiding out of bounds values
987               ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh )
988               ph_s(ji,jl) = MAX( 0._wp, ph_s(ji,jl) - zdh * rhoi * r1_rhos )
989            ENDIF
990         END DO
991      END DO
992      !
993      DEALLOCATE( z1_hti )
994      !
995      ! == temperature and salinity == !
996      DO jl = 1, jpl
997         pt_i (:,jl) = ptmi (:)
998         pt_s (:,jl) = ptms (:)
999         pt_su(:,jl) = ptmsu(:)
1000         ps_i (:,jl) = psmi (:)
1001         ps_i (:,jl) = psmi (:)         
1002      END DO
1003      !
1004      ! == ponds == !
1005      ALLOCATE( zfra(idim) )
1006      ! keep the same pond fraction atip/ati for each category
1007      WHERE( pati(:) /= 0._wp )   ;   zfra(:) = patip(:) / pati(:)
1008      ELSEWHERE                   ;   zfra(:) = 0._wp
1009      END WHERE
1010      DO jl = 1, jpl
1011         pa_ip(:,jl) = zfra(:) * pa_i(:,jl)
1012      END DO
1013      ! keep the same v_ip/v_i ratio for each category
1014      WHERE( ( phti(:) * pati(:) ) /= 0._wp )   ;   zfra(:) = ( phtip(:) * patip(:) ) / ( phti(:) * pati(:) )
1015      ELSEWHERE                                 ;   zfra(:) = 0._wp
1016      END WHERE
1017      DO jl = 1, jpl
1018         WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1019         ELSEWHERE                       ;   ph_ip(:,jl) = 0._wp
1020         END WHERE
1021      END DO
1022      ! keep the same v_il/v_i ratio for each category
1023      WHERE( ( phti(:) * pati(:) ) /= 0._wp )   ;   zfra(:) = ( phtil(:) * patip(:) ) / ( phti(:) * pati(:) )
1024      ELSEWHERE                                 ;   zfra(:) = 0._wp
1025      END WHERE
1026      DO jl = 1, jpl
1027         WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1028         ELSEWHERE                       ;   ph_il(:,jl) = 0._wp
1029         END WHERE
1030      END DO
1031      DEALLOCATE( zfra )
1032      !
1033   END SUBROUTINE ice_var_itd_1cMc
1034
1035   SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,                             ph_i, ph_s, pa_i, &
1036      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
1037      !!-------------------------------------------------------------------
1038      !!
1039      !! ** Purpose :  converting N-cat ice to jpl ice categories
1040      !!
1041      !!                  ice thickness distribution follows a gaussian law
1042      !!               around the concentration of the most likely ice thickness
1043      !!                           (similar as iceistate.F90)
1044      !!
1045      !! ** Method:   Iterative procedure
1046      !!               
1047      !!               1) Fill ice cat that correspond to input thicknesses
1048      !!                  Find the lowest(jlmin) and highest(jlmax) cat that are filled
1049      !!
1050      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1
1051       !!                   by removing 25% ice area from jlmin and jlmax (resp.)
1052      !!             
1053      !!               3) Expand the filling to the empty cat between jlmin and jlmax
1054      !!                   by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
1055      !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
1056      !!
1057      !! ** Arguments : phti: N-cat ice thickness
1058      !!                phts: N-cat snow depth
1059      !!                pati: N-cat ice concentration
1060      !!
1061      !! ** Output    : jpl-cat
1062      !!
1063      !!  (Example of application: BDY forcings when inputs have N-cat /= jpl) 
1064      !!-------------------------------------------------------------------
1065      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
1066      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
1067      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds
1068      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
1069      !
1070      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   jlfil, jlfil2
1071      INTEGER , ALLOCATABLE, DIMENSION(:)   ::   jlmax, jlmin
1072      REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   z1_ai, z1_vi, z1_vs, ztmp, zfra
1073      !
1074      REAL(wp), PARAMETER ::   ztrans = 0.25_wp
1075      INTEGER  ::   ji, jl, jl1, jl2
1076      INTEGER  ::   idim, icat 
1077      !!-------------------------------------------------------------------
1078      !
1079      idim = SIZE( phti, 1 )
1080      icat = SIZE( phti, 2 )
1081      !
1082      ! == thickness and concentration == !
1083      !                                 ! ---------------------- !
1084      IF( icat == jpl ) THEN            ! input cat = output cat !
1085         !                              ! ---------------------- !
1086         ph_i(:,:) = phti(:,:)
1087         ph_s(:,:) = phts(:,:)
1088         pa_i(:,:) = pati(:,:)
1089         !
1090         ! == temperature and salinity and ponds == !
1091         pt_i (:,:) = ptmi (:,:)
1092         pt_s (:,:) = ptms (:,:)
1093         pt_su(:,:) = ptmsu(:,:)
1094         ps_i (:,:) = psmi (:,:)
1095         pa_ip(:,:) = patip(:,:)
1096         ph_ip(:,:) = phtip(:,:)
1097         ph_il(:,:) = phtil(:,:)
1098         !                              ! ---------------------- !
1099      ELSEIF( icat == 1 ) THEN          ! input cat = 1          !
1100         !                              ! ---------------------- !
1101         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), &
1102            &                    ph_i(:,:), ph_s(:,:), pa_i (:,:), &
1103            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), phtil(:,1), &
1104            &                    pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:), ph_il(:,:)  )
1105         !                              ! ---------------------- !
1106      ELSEIF( jpl == 1 ) THEN           ! output cat = 1         !
1107         !                              ! ---------------------- !
1108         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), &
1109            &                    ph_i(:,1), ph_s(:,1), pa_i (:,1), &
1110            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), phtil(:,:), &
1111            &                    pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1), ph_il(:,1)  )
1112         !                              ! ----------------------- !
1113      ELSE                              ! input cat /= output cat !
1114         !                              ! ----------------------- !
1115         
1116         ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays
1117         ALLOCATE( jlmin(idim), jlmax(idim) )
1118
1119         ! --- initialize output fields to 0 --- !
1120         ph_i(1:idim,1:jpl) = 0._wp
1121         ph_s(1:idim,1:jpl) = 0._wp
1122         pa_i(1:idim,1:jpl) = 0._wp
1123         !
1124         ! --- fill the categories --- !
1125         !     find where cat-input = cat-output and fill cat-output fields 
1126         jlmax(:) = 0
1127         jlmin(:) = 999
1128         jlfil(:,:) = 0
1129         DO jl1 = 1, jpl
1130            DO jl2 = 1, icat
1131               DO ji = 1, idim
1132                  IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN
1133                     ! fill the right category
1134                     ph_i(ji,jl1) = phti(ji,jl2)
1135                     ph_s(ji,jl1) = phts(ji,jl2)
1136                     pa_i(ji,jl1) = pati(ji,jl2)
1137                     ! record categories that are filled
1138                     jlmax(ji) = MAX( jlmax(ji), jl1 )
1139                     jlmin(ji) = MIN( jlmin(ji), jl1 )
1140                     jlfil(ji,jl1) = jl1
1141                  ENDIF
1142               END DO
1143            END DO
1144         END DO
1145         !
1146         ! --- fill the gaps between categories --- ! 
1147         !     transfer from categories filled at the previous step to the empty ones in between
1148         DO ji = 1, idim
1149            jl1 = jlmin(ji)
1150            jl2 = jlmax(ji)
1151            IF( jl1 > 1 ) THEN
1152               ! fill the lower cat (jl1-1)
1153               pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1)
1154               ph_i(ji,jl1-1) = hi_mean(jl1-1)
1155               ! remove from cat jl1
1156               pa_i(ji,jl1  ) = ( 1._wp - ztrans ) * pa_i(ji,jl1)
1157            ENDIF
1158            IF( jl2 < jpl ) THEN
1159               ! fill the upper cat (jl2+1)
1160               pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2)
1161               ph_i(ji,jl2+1) = hi_mean(jl2+1)
1162               ! remove from cat jl2
1163               pa_i(ji,jl2  ) = ( 1._wp - ztrans ) * pa_i(ji,jl2)
1164            ENDIF
1165         END DO
1166         !
1167         jlfil2(:,:) = jlfil(:,:) 
1168         ! fill categories from low to high
1169         DO jl = 2, jpl-1
1170            DO ji = 1, idim
1171               IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
1172                  ! fill high
1173                  pa_i(ji,jl) = ztrans * pa_i(ji,jl-1)
1174                  ph_i(ji,jl) = hi_mean(jl)
1175                  jlfil(ji,jl) = jl
1176                  ! remove low
1177                  pa_i(ji,jl-1) = ( 1._wp - ztrans ) * pa_i(ji,jl-1)
1178               ENDIF
1179            END DO
1180         END DO
1181         !
1182         ! fill categories from high to low
1183         DO jl = jpl-1, 2, -1
1184            DO ji = 1, idim
1185               IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
1186                  ! fill low
1187                  pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1)
1188                  ph_i(ji,jl) = hi_mean(jl) 
1189                  jlfil2(ji,jl) = jl
1190                  ! remove high
1191                  pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1)
1192               ENDIF
1193            END DO
1194         END DO
1195         !
1196         DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
1197         DEALLOCATE( jlmin, jlmax )
1198         !
1199         ! == temperature and salinity == !
1200         !
1201         ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) )
1202         !
1203         WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp )               ;   z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 )
1204         ELSEWHERE                                               ;   z1_ai(:) = 0._wp
1205         END WHERE
1206         WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp )   ;   z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 )
1207         ELSEWHERE                                               ;   z1_vi(:) = 0._wp
1208         END WHERE
1209         WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp )   ;   z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 )
1210         ELSEWHERE                                               ;   z1_vs(:) = 0._wp
1211         END WHERE
1212         !
1213         ! fill all the categories with the same value
1214         ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
1215         DO jl = 1, jpl
1216            pt_i (:,jl) = ztmp(:)
1217         END DO
1218         ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
1219         DO jl = 1, jpl
1220            pt_s (:,jl) = ztmp(:)
1221         END DO
1222         ztmp(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:)
1223         DO jl = 1, jpl
1224            pt_su(:,jl) = ztmp(:)
1225         END DO
1226         ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
1227         DO jl = 1, jpl
1228            ps_i (:,jl) = ztmp(:)
1229         END DO
1230         !
1231         DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp )
1232         !
1233         ! == ponds == !
1234         ALLOCATE( zfra(idim) )
1235         ! keep the same pond fraction atip/ati for each category
1236         WHERE( SUM( pati(:,:), dim=2 ) /= 0._wp )   ;   zfra(:) = SUM( patip(:,:), dim=2 ) / SUM( pati(:,:), dim=2 )
1237         ELSEWHERE                                   ;   zfra(:) = 0._wp
1238         END WHERE
1239         DO jl = 1, jpl
1240            pa_ip(:,jl) = zfra(:) * pa_i(:,jl)
1241         END DO
1242         ! keep the same v_ip/v_i ratio for each category
1243         WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp )
1244            zfra(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 )
1245         ELSEWHERE
1246            zfra(:) = 0._wp
1247         END WHERE
1248         DO jl = 1, jpl
1249            WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1250            ELSEWHERE                       ;   ph_ip(:,jl) = 0._wp
1251            END WHERE
1252         END DO
1253         ! keep the same v_il/v_i ratio for each category
1254         WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp )
1255            zfra(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 )
1256         ELSEWHERE
1257            zfra(:) = 0._wp
1258         END WHERE
1259         DO jl = 1, jpl
1260            WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1261            ELSEWHERE                       ;   ph_il(:,jl) = 0._wp
1262            END WHERE
1263         END DO
1264         DEALLOCATE( zfra )
1265         !
1266      ENDIF
1267      !
1268   END SUBROUTINE ice_var_itd_NcMc
1269
1270#else
1271   !!----------------------------------------------------------------------
1272   !!   Default option         Dummy module           NO SI3 sea-ice model
1273   !!----------------------------------------------------------------------
1274#endif
1275
1276   !!======================================================================
1277END MODULE icevar
Note: See TracBrowser for help on using the repository browser.