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

Last change on this file since 13226 was 13226, checked in by orioltp, 3 months ago

Merging dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation into the trunk

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