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

source: NEMO/branches/2020/temporary_r4_trunk/src/ICE/icevar.F90 @ 13466

Last change on this file since 13466 was 13466, checked in by smasson, 4 years ago

r4_trunk: merge r4 13280:13310, see #2523

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