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

source: NEMO/trunk/src/ICE/icevar.F90 @ 14072

Last change on this file since 14072 was 14072, checked in by laurent, 3 years ago

Merging branch "2020/dev_r13648_ASINTER-04_laurent_bulk_ice", ticket #2369

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