source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icevar.F90 @ 11507

Last change on this file since 11507 was 11507, checked in by clem, 14 months ago

add a new functionality for bdy ice by allowing the user to have input files for ice/snw temperature and salinity instead of setting a constant value. Also, input files can have any number of categories

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