source: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icevar.F90 @ 12379

Last change on this file since 12379 was 12379, checked in by dancopsey, 14 months ago

Add meltpond lid thickness as a new prognostic.

File size: 50.8 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   !!                        - et_s(jpi,jpj)  total snow heat content
35   !!                        - et_i(jpi,jpj)  total ice thermal content
36   !!                        - sm_i(jpi,jpj)  mean ice salinity
37   !!                        - tm_i(jpi,jpj)  mean ice temperature
38   !!                        - tm_s(jpi,jpj)  mean snw temperature
39   !!----------------------------------------------------------------------
40   !!   ice_var_agg       : integrate variables over layers and categories
41   !!   ice_var_glo2eqv   : transform from VGLO to VEQV
42   !!   ice_var_eqv2glo   : transform from VEQV to VGLO
43   !!   ice_var_salprof   : salinity profile in the ice
44   !!   ice_var_salprof1d : salinity profile in the ice 1D
45   !!   ice_var_zapsmall  : remove very small area and volume
46   !!   ice_var_zapneg    : remove negative ice fields
47   !!   ice_var_roundoff  : remove negative values arising from roundoff erros
48   !!   ice_var_itd       : convert 1-cat to jpl-cat
49   !!   ice_var_itd2      : convert N-cat to jpl-cat
50   !!   ice_var_bv        : brine volume
51   !!   ice_var_enthalpy  : compute ice and snow enthalpies from temperature
52   !!   ice_var_sshdyn    : compute equivalent ssh in lead
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_itd
76   PUBLIC   ice_var_itd2
77   PUBLIC   ice_var_bv           
78   PUBLIC   ice_var_enthalpy           
79   PUBLIC   ice_var_sshdyn
80
81   !!----------------------------------------------------------------------
82   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
83   !! $Id$
84   !! Software governed by the CeCILL license (see ./LICENSE)
85   !!----------------------------------------------------------------------
86CONTAINS
87
88   SUBROUTINE ice_var_agg( kn )
89      !!-------------------------------------------------------------------
90      !!                ***  ROUTINE ice_var_agg  ***
91      !!
92      !! ** Purpose :   aggregates ice-thickness-category variables to
93      !!              all-ice variables, i.e. it turns VGLO into VAGG
94      !!-------------------------------------------------------------------
95      INTEGER, INTENT( in ) ::   kn     ! =1 state variables only
96      !                                 ! >1 state variables + others
97      !
98      INTEGER ::   ji, jj, jk, jl   ! dummy loop indices
99      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z1_at_i, z1_vt_i, z1_vt_s
100      !!-------------------------------------------------------------------
101      !
102      !                                      ! integrated values
103      vt_i(:,:) =       SUM( v_i(:,:,:)           , dim=3 )
104      vt_s(:,:) =       SUM( v_s(:,:,:)           , dim=3 )
105      at_i(:,:) =       SUM( a_i(:,:,:)           , dim=3 )
106      et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 )
107      et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
108      !
109      at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds
110      vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
111      !
112      ato_i(:,:) = 1._wp - at_i(:,:)         ! open water fraction 
113
114      ! The following fields are calculated for diagnostics and outputs only
115      ! ==> Do not use them for other purposes
116      IF( kn > 1 ) THEN
117         !
118         ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) )
119         WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:)
120         ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp
121         END WHERE
122         WHERE( vt_i(:,:) > epsi20 )   ;   z1_vt_i(:,:) = 1._wp / vt_i(:,:)
123         ELSEWHERE                     ;   z1_vt_i(:,:) = 0._wp
124         END WHERE
125         WHERE( vt_s(:,:) > epsi20 )   ;   z1_vt_s(:,:) = 1._wp / vt_s(:,:)
126         ELSEWHERE                     ;   z1_vt_s(:,:) = 0._wp
127         END WHERE
128         !
129         !                          ! mean ice/snow thickness
130         hm_i(:,:) = vt_i(:,:) * z1_at_i(:,:)
131         hm_s(:,:) = vt_s(:,:) * z1_at_i(:,:)
132         !         
133         !                          ! mean temperature (K), salinity and age
134         tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
135         tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
136         om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(:,:)
137         sm_i (:,:) = SUM( sv_i(:,:,:)              , dim=3 ) * z1_vt_i(:,:)
138         !
139         tm_i(:,:) = 0._wp
140         tm_s(:,:) = 0._wp
141         DO jl = 1, jpl
142            DO jk = 1, nlay_i
143               tm_i(:,:) = tm_i(:,:) + r1_nlay_i * t_i (:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
144            END DO
145            DO jk = 1, nlay_s
146               tm_s(:,:) = tm_s(:,:) + r1_nlay_s * t_s (:,:,jk,jl) * v_s(:,:,jl) * z1_vt_s(:,:)
147            END DO
148         END DO
149         !
150         !                           ! put rt0 where there is no ice
151         WHERE( at_i(:,:)<=epsi20 )
152            tm_su(:,:) = rt0
153            tm_si(:,:) = rt0
154            tm_i (:,:) = rt0
155            tm_s (:,:) = rt0
156         END WHERE
157
158         DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s )
159      ENDIF
160      !
161   END SUBROUTINE ice_var_agg
162
163
164   SUBROUTINE ice_var_glo2eqv
165      !!-------------------------------------------------------------------
166      !!                ***  ROUTINE ice_var_glo2eqv ***
167      !!
168      !! ** Purpose :   computes equivalent variables as function of 
169      !!              global variables, i.e. it turns VGLO into VEQV
170      !!-------------------------------------------------------------------
171      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
172      REAL(wp) ::   ze_i             ! local scalars
173      REAL(wp) ::   ze_s, ztmelts, zbbb, zccc       !   -      -
174      REAL(wp) ::   zhmax, z1_zhmax                 !   -      -
175      REAL(wp) ::   zlay_i, zlay_s                  !   -      -
176      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i
177      !!-------------------------------------------------------------------
178
179!!gm Question 2:  It is possible to define existence of sea-ice in a common way between
180!!                ice area and ice volume ?
181!!                the idea is to be able to define one for all at the begining of this routine
182!!                a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
183
184      !---------------------------------------------------------------
185      ! Ice thickness, snow thickness, ice salinity, ice age and ponds
186      !---------------------------------------------------------------
187      !                                            !--- inverse of the ice area
188      WHERE( a_i(:,:,:) > epsi20 )   ;   z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
189      ELSEWHERE                      ;   z1_a_i(:,:,:) = 0._wp
190      END WHERE
191      !
192      WHERE( v_i(:,:,:) > epsi20 )   ;   z1_v_i(:,:,:) = 1._wp / v_i(:,:,:)
193      ELSEWHERE                      ;   z1_v_i(:,:,:) = 0._wp
194      END WHERE
195      !                                           !--- ice thickness
196      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
197
198      zhmax    =          hi_max(jpl)
199      z1_zhmax =  1._wp / hi_max(jpl)               
200      WHERE( h_i(:,:,jpl) > zhmax )   ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area
201         h_i  (:,:,jpl) = zhmax
202         a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 
203         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)
204      END WHERE
205      !                                           !--- snow thickness
206      h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:)
207      !                                           !--- ice age     
208      o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:)
209      !                                           !--- pond fraction and thickness     
210      a_ip_frac(:,:,:) = a_ip(:,:,:) * z1_a_i(:,:,:)
211      WHERE( a_ip_frac(:,:,:) > epsi20 )   ;   h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:)
212      ELSEWHERE                            ;   h_ip(:,:,:) = 0._wp
213      END WHERE
214      !
215      !                                           !---  salinity (with a minimum value imposed everywhere)     
216      IF( nn_icesal == 2 ) THEN
217         WHERE( v_i(:,:,:) > epsi20 )   ;   s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) * z1_v_i(:,:,:) ) )
218         ELSEWHERE                      ;   s_i(:,:,:) = rn_simin
219         END WHERE
220      ENDIF
221
222      write(numout,*) 'icevar: sv_i(3,4,1) = ', sv_i(3,4,1)
223      write(numout,*) 'icevar: z1_v_i(3,4,1) = ', z1_v_i(3,4,1)
224
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
248                  If (to_print_2d(ji,jj) == 10 .AND. jk == 1 .AND. jl == 1) THEN
249                     write(numout,*)'ice_var_glo2eqv 0: t_i, e_i = ',t_i(ji,jj,jk,jl) - rt0, ' ', e_i (ji,jj,jk,jl)
250                  END IF
251
252               END DO
253            END DO
254         END DO
255      END DO
256
257      !--------------------
258      ! Snow temperature   [K]   (with a minimum value (rt0 - 100.))
259      !--------------------
260      zlay_s = REAL( nlay_s , wp )
261      DO jk = 1, nlay_s
262         WHERE( v_s(:,:,:) > epsi20 )        !--- icy area
263            t_s(:,:,jk,:) = rt0 + MAX( -100._wp ,  &
264                 &                MIN( r1_rcpi * ( -r1_rhos * ( e_s(:,:,jk,:) / v_s(:,:,:) * zlay_s ) + rLfus ) , 0._wp ) )
265         ELSEWHERE                           !--- no ice
266            t_s(:,:,jk,:) = rt0
267         END WHERE
268
269         If (to_print_2d(ji,jj) == 10 .AND. jk == 1 .AND. jl == 1) THEN
270             write(numout,*)'ice_var_glo2eqv 0: t_s, e_s = ',t_s(ji,jj,jk,jl), ' ', e_s (ji,jj,jk,jl)
271         END IF
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      REAL(wp), PARAMETER :: min_salin = 0.1_wp     ! Minimum salinity in sea ice
321      REAL(wp) ::            zn                     ! Fraction of ice depth
322      REAL(wp), PARAMETER :: msal    = 0.573_wp
323      REAL(wp), PARAMETER :: nsal    = 0.407_wp
324      REAL(wp), PARAMETER :: saltmax = 9.6_wp
325      ! REAL(wp), PARAMETER, DIMENSION(4) :: weights = (/0.06363419, 0.25545967, 0.33155539, 0.34935076/)
326      REAL(wp), PARAMETER, DIMENSION(4) :: weights = (/0.0703805, 0.25526203, 0.32860314, 0.34575433/)
327     
328      !!-------------------------------------------------------------------
329
330!!gm Question: Remove the option 3 ?  How many years since it last use ?
331
332      SELECT CASE ( nn_icesal )
333      !
334      !               !---------------------------------------!
335      CASE( 1 )       !  constant salinity in time and space  !
336         !            !---------------------------------------!
337         sz_i(:,:,:,:) = rn_icesal
338         s_i (:,:,:)   = rn_icesal
339         !
340         !            !---------------------------------------------!
341      CASE( 2 )       !  time varying salinity with linear profile  !
342         !            !---------------------------------------------!
343         !
344         ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
345         !
346         DO jl = 1, jpl
347            DO jk = 1, nlay_i
348               sz_i(:,:,jk,jl)  = s_i(:,:,jl)
349            END DO
350         END DO
351
352         !                                      ! Slope of the linear profile
353         WHERE( h_i(:,:,:) > epsi20 )   ;   z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:)
354         ELSEWHERE                      ;   z_slope_s(:,:,:) = 0._wp
355         END WHERE
356         !
357         z1_dS = 1._wp / ( zsi1 - zsi0 )
358         DO jl = 1, jpl
359            DO jj = 1, jpj
360               DO ji = 1, jpi
361                  zalpha(ji,jj,jl) = 1._wp
362                  !                             ! force a constant profile when SSS too low (Baltic Sea)
363                  IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) )   zalpha(ji,jj,jl) = 0._wp 
364               END DO
365            END DO
366         END DO
367         
368         !
369         ! Computation of the profile
370         DO jl = 1, jpl
371            DO jk = 1, nlay_i
372               DO jj = 1, jpj
373                  DO ji = 1, jpi
374
375                     !  Copied over from CICE
376                     ! zn = (real(jk,kind=dbl_kind)-0.5_wp) / real(nlay_i,kind=dbl_kind)
377                     ! sz_i(ji,jj,jk,jl)=(saltmax/2.0_wp)*(1.0_wp-cos(pi*zn**(nsal/(msal+zn))))
378                     ! sz_i(ji,jj,jk,jl) = max(sz_i(ji,jj,jk,jl), min_salin)
379
380                     ! Weighting method to match CICE
381                     zs0 = (s_i(ji,jj,jl)*4.0_wp) * weights(jk)
382                     zs  = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl)     ! weighting the profile
383                     sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
384
385
386                     ! zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
387                     ! zs  = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl)     ! weighting the profile
388                     ! sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
389
390                     IF ( jl == 1 .AND. jj == 26 .AND. ji == 42 ) THEN
391                       write(numout,*) 'icevar: jk, sz_i(ji,jj,jk,jl), s_i(ji,jj,jl), weights(jk) = ',jk, ' ',sz_i(ji,jj,jk,jl), ' ',s_i(ji,jj,jl), ' ',weights(jk)
392                     ENDIF
393
394                  END DO
395               END DO
396            END DO
397         END DO
398         !
399         DEALLOCATE( z_slope_s , zalpha )
400         !
401         !            !-------------------------------------------!
402      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
403         !            !-------------------------------------------!                                   (mean = 2.30)
404         !
405         s_i(:,:,:) = 2.30_wp
406!!gm Remark: if we keep the case 3, then compute an store one for all time-step
407!!           a array  S_prof(1:nlay_i) containing the calculation and just do:
408!         DO jk = 1, nlay_i
409!            sz_i(:,:,jk,:) = S_prof(jk)
410!         END DO
411!!gm end
412         !
413         DO jl = 1, jpl
414            DO jk = 1, nlay_i
415               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
416               sz_i(:,:,jk,jl) =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
417            END DO
418         END DO
419         !
420      END SELECT
421      !
422   END SUBROUTINE ice_var_salprof
423
424
425   SUBROUTINE ice_var_salprof1d
426      !!-------------------------------------------------------------------
427      !!                  ***  ROUTINE ice_var_salprof1d  ***
428      !!
429      !! ** Purpose :   1d computation of the sea ice salinity profile
430      !!                Works with 1d vectors and is used by thermodynamic modules
431      !!-------------------------------------------------------------------
432      INTEGER  ::   ji, jk    ! dummy loop indices
433      REAL(wp) ::   zargtemp, zsal, z1_dS   ! local scalars
434      REAL(wp) ::   zs, zs0              !   -      -
435      !
436      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z_slope_s, zalpha   !
437      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
438      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
439      !!-------------------------------------------------------------------
440      !
441      SELECT CASE ( nn_icesal )
442      !
443      !               !---------------------------------------!
444      CASE( 1 )       !  constant salinity in time and space  !
445         !            !---------------------------------------!
446         sz_i_1d(1:npti,:) = rn_icesal
447         !
448         !            !---------------------------------------------!
449      CASE( 2 )       !  time varying salinity with linear profile  !
450         !            !---------------------------------------------!
451         !
452         ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
453         !
454         !                                      ! Slope of the linear profile
455         WHERE( h_i_1d(1:npti) > epsi20 )   ;   z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
456         ELSEWHERE                          ;   z_slope_s(1:npti) = 0._wp
457         END WHERE
458         
459         z1_dS = 1._wp / ( zsi1 - zsi0 )
460         DO ji = 1, npti
461            zalpha(ji) = MAX(  0._wp , MIN(  ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp  )  )
462            !                             ! force a constant profile when SSS too low (Baltic Sea)
463            IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) )   zalpha(ji) = 0._wp
464         END DO
465         !
466         ! Computation of the profile
467         DO jk = 1, nlay_i
468            DO ji = 1, npti
469               !                          ! linear profile with 0 surface value
470               zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
471               zs  = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
472               sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
473            END DO
474         END DO
475         !
476         DEALLOCATE( z_slope_s, zalpha )
477
478         !            !-------------------------------------------!
479      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
480         !            !-------------------------------------------!                                   (mean = 2.30)
481         !
482         s_i_1d(1:npti) = 2.30_wp
483         !
484!!gm cf remark in ice_var_salprof routine, CASE( 3 )
485         DO jk = 1, nlay_i
486            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
487            zsal =  1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
488            DO ji = 1, npti
489               sz_i_1d(ji,jk) = zsal
490            END DO
491         END DO
492         !
493      END SELECT
494      !
495   END SUBROUTINE ice_var_salprof1d
496
497
498   SUBROUTINE ice_var_zapsmall
499      !!-------------------------------------------------------------------
500      !!                   ***  ROUTINE ice_var_zapsmall ***
501      !!
502      !! ** Purpose :   Remove too small sea ice areas and correct fluxes
503      !!-------------------------------------------------------------------
504      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
505      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch
506      !!-------------------------------------------------------------------
507      !
508      DO jl = 1, jpl       !==  loop over the categories  ==!
509         !
510         WHERE( a_i(:,:,jl) > epsi10 )   ;   h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
511         ELSEWHERE                       ;   h_i(:,:,jl) = 0._wp
512         END WHERE
513         !
514         WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 )   ;   zswitch(:,:) = 0._wp
515         ELSEWHERE                                                                           ;   zswitch(:,:) = 1._wp
516         END WHERE
517         !
518         !-----------------------------------------------------------------
519         ! Zap ice energy and use ocean heat to melt ice
520         !-----------------------------------------------------------------
521         DO jk = 1, nlay_i
522            DO jj = 1 , jpj
523               DO ji = 1 , jpi
524                  ! update exchanges with ocean
525                  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
526                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
527                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
528               END DO
529            END DO
530         END DO
531         !
532         DO jk = 1, nlay_s
533            DO jj = 1 , jpj
534               DO ji = 1 , jpi
535                  ! update exchanges with ocean
536                  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
537                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj)
538                  t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
539
540                  If (to_print_2d(ji,jj) == 10 .AND. jk == 1 .AND. jl == 1) THEN
541                    write(numout,*)'ice_var_glo2eqv 1: t_s, e_s = ',t_s(ji,jj,jk,jl), ' ', e_s (ji,jj,jk,jl)
542                  END IF
543               END DO
544            END DO
545         END DO
546         !
547         !-----------------------------------------------------------------
548         ! zap ice and snow volume, add water and salt to ocean
549         !-----------------------------------------------------------------
550         DO jj = 1 , jpj
551            DO ji = 1 , jpi
552               ! update exchanges with ocean
553               sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoi * r1_rdtice
554               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoi * r1_rdtice
555               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhos * r1_rdtice
556               !
557               a_i  (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
558               v_i  (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
559               v_s  (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
560               t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
561               oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
562               sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
563               !
564               h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
565               h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
566               !
567               a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
568               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
569               lh_ip (ji,jj,jl) = lh_ip (ji,jj,jl) * zswitch(ji,jj)
570               !
571            END DO
572         END DO
573         !
574      END DO 
575
576      ! to be sure that at_i is the sum of a_i(jl)
577      at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
578      vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
579
580      ! open water = 1 if at_i=0
581      WHERE( at_i(:,:) == 0._wp )   ato_i(:,:) = 1._wp
582      !
583   END SUBROUTINE ice_var_zapsmall
584
585
586   SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
587      !!-------------------------------------------------------------------
588      !!                   ***  ROUTINE ice_var_zapneg ***
589      !!
590      !! ** Purpose :   Remove negative sea ice fields and correct fluxes
591      !!-------------------------------------------------------------------
592      REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step
593      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
594      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
595      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
596      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_i      ! salt content
597      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
598      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
599      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
600      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
601      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness
602      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
603      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
604      !
605      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
606      REAL(wp) ::   z1_dt
607      !!-------------------------------------------------------------------
608      !
609      z1_dt = 1._wp / pdt
610      !
611      DO jl = 1, jpl       !==  loop over the categories  ==!
612         !
613         ! make sure a_i=0 where v_i<=0
614         WHERE( pv_i(:,:,:) <= 0._wp )   pa_i(:,:,:) = 0._wp
615
616         !----------------------------------------
617         ! zap ice energy and send it to the ocean
618         !----------------------------------------
619         DO jk = 1, nlay_i
620            DO jj = 1 , jpj
621               DO ji = 1 , jpi
622                  IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
623                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
624                     pe_i(ji,jj,jk,jl) = 0._wp
625                  ENDIF
626               END DO
627            END DO
628         END DO
629         !
630         DO jk = 1, nlay_s
631            DO jj = 1 , jpj
632               DO ji = 1 , jpi
633                  IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
634                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
635                     pe_s(ji,jj,jk,jl) = 0._wp
636                  ENDIF
637               END DO
638            END DO
639         END DO
640         !
641         !-----------------------------------------------------
642         ! zap ice and snow volume, add water and salt to ocean
643         !-----------------------------------------------------
644         DO jj = 1 , jpj
645            DO ji = 1 , jpi
646               IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
647                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
648                  pv_i   (ji,jj,jl) = 0._wp
649               ENDIF
650               IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
651                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
652                  pv_s   (ji,jj,jl) = 0._wp
653               ENDIF
654               IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
655                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
656                  psv_i  (ji,jj,jl) = 0._wp
657               ENDIF
658            END DO
659         END DO
660         !
661      END DO 
662      !
663      WHERE( pato_i(:,:)   < 0._wp )   pato_i(:,:)   = 0._wp
664      WHERE( poa_i (:,:,:) < 0._wp )   poa_i (:,:,:) = 0._wp
665      WHERE( pa_i  (:,:,:) < 0._wp )   pa_i  (:,:,:) = 0._wp
666      WHERE( pa_ip (:,:,:) < 0._wp )   pa_ip (:,:,:) = 0._wp
667      WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
668      !                                                        but it does not change conservation, so keep it this way is ok
669      WHERE( plh_ip (:,:,:) < 0._wp )   plh_ip (:,:,:) = 0._wp
670      !
671   END SUBROUTINE ice_var_zapneg
672
673
674   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
675      !!-------------------------------------------------------------------
676      !!                   ***  ROUTINE ice_var_roundoff ***
677      !!
678      !! ** Purpose :   Remove negative sea ice values arising from roundoff errors
679      !!-------------------------------------------------------------------
680      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
681      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_i       ! ice volume
682      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_s       ! snw volume
683      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   psv_i      ! salt content
684      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   poa_i      ! age content
685      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
686      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
687      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness
688      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
689      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
690      !!-------------------------------------------------------------------
691      !
692      WHERE( pa_i (1:npti,:)   < 0._wp .AND. pa_i (1:npti,:)   > -epsi10 )   pa_i (1:npti,:)   = 0._wp   !  a_i must be >= 0
693      WHERE( pv_i (1:npti,:)   < 0._wp .AND. pv_i (1:npti,:)   > -epsi10 )   pv_i (1:npti,:)   = 0._wp   !  v_i must be >= 0
694      WHERE( pv_s (1:npti,:)   < 0._wp .AND. pv_s (1:npti,:)   > -epsi10 )   pv_s (1:npti,:)   = 0._wp   !  v_s must be >= 0
695      WHERE( psv_i(1:npti,:)   < 0._wp .AND. psv_i(1:npti,:)   > -epsi10 )   psv_i(1:npti,:)   = 0._wp   ! sv_i must be >= 0
696      WHERE( poa_i(1:npti,:)   < 0._wp .AND. poa_i(1:npti,:)   > -epsi10 )   poa_i(1:npti,:)   = 0._wp   ! oa_i must be >= 0
697      WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0
698      WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0
699      IF ( ln_pnd_H12 ) THEN
700         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0
701         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0
702         WHERE( plh_ip(1:npti,:) < 0._wp .AND. plh_ip(1:npti,:) > -epsi10 )    plh_ip(1:npti,:)   = 0._wp   ! lh_ip must be >= 0
703      ENDIF
704      !
705   END SUBROUTINE ice_var_roundoff
706   
707   
708   SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
709      !!-------------------------------------------------------------------
710      !!                ***  ROUTINE ice_var_itd   ***
711      !!
712      !! ** Purpose :  converting 1-cat ice to multiple ice categories
713      !!
714      !!                  ice thickness distribution follows a gaussian law
715      !!               around the concentration of the most likely ice thickness
716      !!                           (similar as iceistate.F90)
717      !!
718      !! ** Method:   Iterative procedure
719      !!               
720      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
721      !!
722      !!               2) Check whether the distribution conserves area and volume, positivity and
723      !!                  category boundaries
724      !!             
725      !!               3) If not (input ice is too thin), the last category is empty and
726      !!                  the number of categories is reduced (jpl-1)
727      !!
728      !!               4) Iterate until ok (SUM(itest(:) = 4)
729      !!
730      !! ** Arguments : zhti: 1-cat ice thickness
731      !!                zhts: 1-cat snow depth
732      !!                zati: 1-cat ice concentration
733      !!
734      !! ** Output    : jpl-cat
735      !!
736      !!  (Example of application: BDY forcings when input are cell averaged) 
737      !!-------------------------------------------------------------------
738      INTEGER  :: ji, jk, jl             ! dummy loop indices
739      INTEGER  :: idim, i_fill, jl0 
740      REAL(wp) :: zarg, zV, zconv, zdh, zdv
741      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input  ice/snow variables
742      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
743      INTEGER , DIMENSION(4)                  ::   itest
744      !!-------------------------------------------------------------------
745      !
746      ! ----------------------------------------
747      ! distribution over the jpl ice categories
748      ! ----------------------------------------
749      ! a gaussian distribution for ice concentration is used
750      ! then we check whether the distribution fullfills
751      ! volume and area conservation, positivity and ice categories bounds
752      idim = SIZE( zhti , 1 )
753      zh_i(1:idim,1:jpl) = 0._wp
754      zh_s(1:idim,1:jpl) = 0._wp
755      za_i(1:idim,1:jpl) = 0._wp
756      !
757      DO ji = 1, idim
758         !
759         IF( zhti(ji) > 0._wp ) THEN
760            !
761            ! find which category (jl0) the input ice thickness falls into
762            jl0 = jpl
763            DO jl = 1, jpl
764               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
765                  jl0 = jl
766                  CYCLE
767               ENDIF
768            END DO
769            !
770            itest(:) = 0
771            i_fill   = jpl + 1                                            !------------------------------------
772            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
773               !                                                          !------------------------------------
774               i_fill = i_fill - 1
775               !
776               zh_i(ji,1:jpl) = 0._wp
777               za_i(ji,1:jpl) = 0._wp
778               itest(:)       = 0     
779               !
780               IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
781                  zh_i(ji,1) = zhti(ji)
782                  za_i (ji,1) = zati (ji)
783               ELSE                         !-- case ice is thicker: fill categories >1
784                  ! thickness
785                  DO jl = 1, i_fill - 1
786                     zh_i(ji,jl) = hi_mean(jl)
787                  END DO
788                  !
789                  ! concentration
790                  za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
791                  DO jl = 1, i_fill - 1
792                     IF ( jl /= jl0 ) THEN
793                        zarg        = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
794                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
795                     ENDIF
796                  END DO
797                  !
798                  ! last category
799                  za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
800                  zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
801                  zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
802                  !
803                  ! correction if concentration of upper cat is greater than lower cat
804                  !    (it should be a gaussian around jl0 but sometimes it is not)
805                  IF ( jl0 /= jpl ) THEN
806                     DO jl = jpl, jl0+1, -1
807                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
808                           zdv = zh_i(ji,jl) * za_i(ji,jl)
809                           zh_i(ji,jl    ) = 0._wp
810                           za_i (ji,jl    ) = 0._wp
811                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
812                        END IF
813                     END DO
814                  ENDIF
815                  !
816               ENDIF
817               !
818               ! Compatibility tests
819               zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) ) 
820               IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation
821               !
822               zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )
823               IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation
824               !
825               IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ?
826               !
827               itest(4) = 1
828               DO jl = 1, i_fill
829                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0                                ! Test 4: positivity of ice concentrations
830               END DO
831               !                                         !----------------------------
832            END DO                                       ! end iteration on categories
833            !                                            !----------------------------
834         ENDIF
835      END DO
836
837      ! Add Snow in each category where za_i is not 0
838      DO jl = 1, jpl
839         DO ji = 1, idim
840            IF( za_i(ji,jl) > 0._wp ) THEN
841               zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
842               ! In case snow load is in excess that would lead to transformation from snow to ice
843               ! Then, transfer the snow excess into the ice (different from icethd_dh)
844               zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 ) 
845               ! recompute h_i, h_s avoiding out of bounds values
846               zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
847               zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos )
848            ENDIF
849         END DO
850      END DO
851      !
852   END SUBROUTINE ice_var_itd
853
854
855   SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )
856      !!-------------------------------------------------------------------
857      !!                ***  ROUTINE ice_var_itd2   ***
858      !!
859      !! ** Purpose :  converting N-cat ice to jpl ice categories
860      !!
861      !!                  ice thickness distribution follows a gaussian law
862      !!               around the concentration of the most likely ice thickness
863      !!                           (similar as iceistate.F90)
864      !!
865      !! ** Method:   Iterative procedure
866      !!               
867      !!               1) Fill ice cat that correspond to input thicknesses
868      !!                  Find the lowest(jlmin) and highest(jlmax) cat that are filled
869      !!
870      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1
871       !!                   by removing 25% ice area from jlmin and jlmax (resp.)
872      !!             
873      !!               3) Expand the filling to the empty cat between jlmin and jlmax
874      !!                   by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
875      !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
876      !!
877      !! ** Arguments : zhti: N-cat ice thickness
878      !!                zhts: N-cat snow depth
879      !!                zati: N-cat ice concentration
880      !!
881      !! ** Output    : jpl-cat
882      !!
883      !!  (Example of application: BDY forcings when inputs have N-cat /= jpl) 
884      !!-------------------------------------------------------------------
885      INTEGER  ::   ji, jl, jl1, jl2             ! dummy loop indices
886      INTEGER  ::   idim, icat 
887      REAL(wp), PARAMETER ::   ztrans = 0.25_wp
888      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
889      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
890      INTEGER , DIMENSION(:,:), ALLOCATABLE   ::   jlfil, jlfil2
891      INTEGER , DIMENSION(:)  , ALLOCATABLE   ::   jlmax, jlmin
892      !!-------------------------------------------------------------------
893      !
894      idim = SIZE( zhti, 1 )
895      icat = SIZE( zhti, 2 )
896      !
897      ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays
898      ALLOCATE( jlmin(idim), jlmax(idim) )
899
900      ! --- initialize output fields to 0 --- !
901      zh_i(1:idim,1:jpl) = 0._wp
902      zh_s(1:idim,1:jpl) = 0._wp
903      za_i(1:idim,1:jpl) = 0._wp
904      !
905      ! --- fill the categories --- !
906      !     find where cat-input = cat-output and fill cat-output fields 
907      jlmax(:) = 0
908      jlmin(:) = 999
909      jlfil(:,:) = 0
910      DO jl1 = 1, jpl
911         DO jl2 = 1, icat
912            DO ji = 1, idim
913               IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
914                  ! fill the right category
915                  zh_i(ji,jl1) = zhti(ji,jl2)
916                  zh_s(ji,jl1) = zhts(ji,jl2)
917                  za_i(ji,jl1) = zati(ji,jl2)
918                  ! record categories that are filled
919                  jlmax(ji) = MAX( jlmax(ji), jl1 )
920                  jlmin(ji) = MIN( jlmin(ji), jl1 )
921                  jlfil(ji,jl1) = jl1
922               ENDIF
923            END DO
924         END DO
925      END DO
926      !
927      ! --- fill the gaps between categories --- ! 
928      !     transfer from categories filled at the previous step to the empty ones in between
929      DO ji = 1, idim
930         jl1 = jlmin(ji)
931         jl2 = jlmax(ji)
932         IF( jl1 > 1 ) THEN
933            ! fill the lower cat (jl1-1)
934            za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
935            zh_i(ji,jl1-1) = hi_mean(jl1-1)
936            ! remove from cat jl1
937            za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
938         ENDIF
939         IF( jl2 < jpl ) THEN
940            ! fill the upper cat (jl2+1)
941            za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
942            zh_i(ji,jl2+1) = hi_mean(jl2+1)
943            ! remove from cat jl2
944            za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
945         ENDIF
946      END DO
947      !
948      jlfil2(:,:) = jlfil(:,:) 
949      ! fill categories from low to high
950      DO jl = 2, jpl-1
951         DO ji = 1, idim
952            IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
953               ! fill high
954               za_i(ji,jl) = ztrans * za_i(ji,jl-1)
955               zh_i(ji,jl) = hi_mean(jl)
956               jlfil(ji,jl) = jl
957               ! remove low
958               za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
959            ENDIF
960         END DO
961      END DO
962      !
963      ! fill categories from high to low
964      DO jl = jpl-1, 2, -1
965         DO ji = 1, idim
966            IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
967               ! fill low
968               za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
969               zh_i(ji,jl) = hi_mean(jl) 
970               jlfil2(ji,jl) = jl
971               ! remove high
972               za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
973            ENDIF
974         END DO
975      END DO
976      !
977      DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
978      DEALLOCATE( jlmin, jlmax )
979      !
980   END SUBROUTINE ice_var_itd2
981
982
983   SUBROUTINE ice_var_bv
984      !!-------------------------------------------------------------------
985      !!                ***  ROUTINE ice_var_bv ***
986      !!
987      !! ** Purpose :   computes mean brine volume (%) in sea ice
988      !!
989      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
990      !!
991      !! References : Vancoppenolle et al., JGR, 2007
992      !!-------------------------------------------------------------------
993      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
994      !!-------------------------------------------------------------------
995      !
996!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
997!!   instead of setting everything to zero as just below
998      bv_i (:,:,:) = 0._wp
999      DO jl = 1, jpl
1000         DO jk = 1, nlay_i
1001            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
1002               bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
1003            END WHERE
1004         END DO
1005      END DO
1006      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
1007      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
1008      END WHERE
1009      !
1010   END SUBROUTINE ice_var_bv
1011
1012
1013   SUBROUTINE ice_var_enthalpy
1014      !!-------------------------------------------------------------------
1015      !!                   ***  ROUTINE ice_var_enthalpy ***
1016      !!                 
1017      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
1018      !!
1019      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
1020      !!-------------------------------------------------------------------
1021      INTEGER  ::   ji, jk   ! dummy loop indices
1022      REAL(wp) ::   ztmelts  ! local scalar
1023      !!-------------------------------------------------------------------
1024      !
1025      DO jk = 1, nlay_i             ! Sea ice energy of melting
1026         DO ji = 1, npti
1027            ztmelts      = - rTmlt  * sz_i_1d(ji,jk)
1028            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
1029                                                                !   (sometimes zdf scheme produces abnormally high temperatures)   
1030            e_i_1d(ji,jk) = rhoi * ( rcpi  * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
1031               &                   + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
1032               &                   - rcp   * ztmelts )
1033         END DO
1034      END DO
1035      DO jk = 1, nlay_s             ! Snow energy of melting
1036         DO ji = 1, npti
1037            e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
1038         END DO
1039      END DO
1040      !
1041   END SUBROUTINE ice_var_enthalpy
1042
1043   FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
1044      !!---------------------------------------------------------------------
1045      !!                   ***  ROUTINE ice_var_sshdyn  ***
1046      !!                     
1047      !! ** Purpose :  compute the equivalent ssh in lead when sea ice is embedded
1048      !!
1049      !! ** Method  :  ssh_lead = ssh + (Mice + Msnow) / rau0
1050      !!
1051      !! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
1052      !!                Sea ice-ocean coupling using a rescaled vertical coordinate z*,
1053      !!                Ocean Modelling, Volume 24, Issues 1-2, 2008
1054      !!----------------------------------------------------------------------
1055      !
1056      ! input
1057      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh            !: ssh [m]
1058      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass    !: mass of snow and ice at current  ice time step [Kg/m2]
1059      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b  !: mass of snow and ice at previous ice time step [Kg/m2]
1060      !
1061      ! output
1062      REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn  ! equivalent ssh in lead [m]
1063      !
1064      ! temporary
1065      REAL(wp) :: zintn, zintb                     ! time interpolation weights []
1066      REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload  ! snow and ice load [m]
1067      !
1068      ! compute ice load used to define the equivalent ssh in lead
1069      IF( ln_ice_embd ) THEN
1070         !                                           
1071         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
1072         !                                               = (1/nn_fsbc)^2 * {SUM[n]        , n=0,nn_fsbc-1}
1073         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
1074         !
1075         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   *    {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
1076         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
1077         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
1078         !
1079         zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0
1080         !
1081      ELSE
1082         zsnwiceload(:,:) = 0.0_wp
1083      ENDIF
1084      ! compute equivalent ssh in lead
1085      ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:)
1086      !
1087   END FUNCTION ice_var_sshdyn
1088
1089
1090#else
1091   !!----------------------------------------------------------------------
1092   !!   Default option         Dummy module           NO SI3 sea-ice model
1093   !!----------------------------------------------------------------------
1094#endif
1095
1096   !!======================================================================
1097END MODULE icevar
Note: See TracBrowser for help on using the repository browser.