source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/ICE/icevar.F90 @ 12154

Last change on this file since 12154 was 12154, checked in by cetlod, 10 months ago

commit

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