New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icevar.F90 in NEMO/branches/2019/fix_vvl_ticket1791/src/ICE – NEMO

source: NEMO/branches/2019/fix_vvl_ticket1791/src/ICE/icevar.F90 @ 11422

Last change on this file since 11422 was 11422, checked in by jchanut, 5 years ago

#1791, merge with trunk

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