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

Last change on this file since 11560 was 11560, checked in by clem, 13 months ago

correct the conversion from 1cat to Ncat for ice data at the boundaries when using bdy

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