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

Last change on this file since 11223 was 11223, checked in by smasson, 2 years ago

dev_r10984_HPC-13 : cleaning of rewriting of bdydta, see #2285

  • Property svn:keywords set to Id
File size: 51.7 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      IF( jpl == 1 ) THEN
856         CALL ice_var_itd_1c1c( zhti, zhts, zati, zh_i(:,1), zh_s(:,1), za_i(:,1) )
857         RETURN
858      ENDIF
859      !
860      DO ji = 1, idim
861         !
862         IF( zhti(ji) > 0._wp ) THEN
863            !
864            ! find which category (jl0) the input ice thickness falls into
865            jl0 = jpl
866            DO jl = 1, jpl
867               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
868                  jl0 = jl
869                  CYCLE
870               ENDIF
871            END DO
872            !
873            itest(:) = 0
874            i_fill   = jpl + 1                                            !------------------------------------
875            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
876               !                                                          !------------------------------------
877               i_fill = i_fill - 1
878               !
879               zh_i(ji,1:jpl) = 0._wp
880               za_i(ji,1:jpl) = 0._wp
881               itest(:)       = 0     
882               !
883               IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
884                  zh_i(ji,1) = zhti(ji)
885                  za_i (ji,1) = zati (ji)
886               ELSE                         !-- case ice is thicker: fill categories >1
887                  ! thickness
888                  DO jl = 1, i_fill - 1
889                     zh_i(ji,jl) = hi_mean(jl)
890                  END DO
891                  !
892                  ! concentration
893                  za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
894                  DO jl = 1, i_fill - 1
895                     IF ( jl /= jl0 ) THEN
896                        zarg        = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
897                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
898                     ENDIF
899                  END DO
900                  !
901                  ! last category
902                  za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
903                  zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
904                  zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
905                  !
906                  ! correction if concentration of upper cat is greater than lower cat
907                  !    (it should be a gaussian around jl0 but sometimes it is not)
908                  IF ( jl0 /= jpl ) THEN
909                     DO jl = jpl, jl0+1, -1
910                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
911                           zdv = zh_i(ji,jl) * za_i(ji,jl)
912                           zh_i(ji,jl    ) = 0._wp
913                           za_i (ji,jl    ) = 0._wp
914                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
915                        END IF
916                     END DO
917                  ENDIF
918                  !
919               ENDIF
920               !
921               ! Compatibility tests
922               zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) ) 
923               IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation
924               !
925               zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )
926               IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation
927               !
928               IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ?
929               !
930               itest(4) = 1
931               DO jl = 1, i_fill
932                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0                                ! Test 4: positivity of ice concentrations
933               END DO
934               !                                         !----------------------------
935            END DO                                       ! end iteration on categories
936            !                                            !----------------------------
937         ENDIF
938      END DO
939
940      ! Add Snow in each category where za_i is not 0
941      DO jl = 1, jpl
942         DO ji = 1, idim
943            IF( za_i(ji,jl) > 0._wp ) THEN
944               zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
945               ! In case snow load is in excess that would lead to transformation from snow to ice
946               ! Then, transfer the snow excess into the ice (different from icethd_dh)
947               zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 ) 
948               ! recompute h_i, h_s avoiding out of bounds values
949               zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
950               zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos )
951            ENDIF
952         END DO
953      END DO
954      !
955   END SUBROUTINE ice_var_itd_1cMc
956
957   SUBROUTINE ice_var_itd_NcMc( zhti, zhts, zati, zh_i, zh_s, za_i )
958      !!-------------------------------------------------------------------
959      !!
960      !! ** Purpose :  converting N-cat ice to jpl ice categories
961      !!
962      !!                  ice thickness distribution follows a gaussian law
963      !!               around the concentration of the most likely ice thickness
964      !!                           (similar as iceistate.F90)
965      !!
966      !! ** Method:   Iterative procedure
967      !!               
968      !!               1) Fill ice cat that correspond to input thicknesses
969      !!                  Find the lowest(jlmin) and highest(jlmax) cat that are filled
970      !!
971      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1
972       !!                   by removing 25% ice area from jlmin and jlmax (resp.)
973      !!             
974      !!               3) Expand the filling to the empty cat between jlmin and jlmax
975      !!                   by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
976      !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
977      !!
978      !! ** Arguments : zhti: N-cat ice thickness
979      !!                zhts: N-cat snow depth
980      !!                zati: N-cat ice concentration
981      !!
982      !! ** Output    : jpl-cat
983      !!
984      !!  (Example of application: BDY forcings when inputs have N-cat /= jpl) 
985      !!-------------------------------------------------------------------
986      INTEGER  ::   ji, jl, jl1, jl2             ! dummy loop indices
987      INTEGER  ::   idim, icat 
988      REAL(wp), PARAMETER ::   ztrans = 0.25_wp
989      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
990      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
991      INTEGER , DIMENSION(:,:), ALLOCATABLE   ::   jlfil, jlfil2
992      INTEGER , DIMENSION(:)  , ALLOCATABLE   ::   jlmax, jlmin
993      !!-------------------------------------------------------------------
994      !
995      idim = SIZE( zhti, 1 )
996      icat = SIZE( zhti, 2 )
997      !                                 ! ---------------------- !
998      IF( icat == jpl ) THEN            ! input cat = output cat !
999         !                              ! ---------------------- !
1000         zh_i(:,:) = zhti(:,:)
1001         zh_s(:,:) = zhts(:,:)
1002         za_i(:,:) = zati(:,:)
1003         !                              ! ---------------------- !
1004      ELSEIF( icat == 1 ) THEN          ! specific case if N = 1 !
1005         !                              ! ---------------------- !
1006         !
1007         CALL ice_var_itd_1cMc( zhti(:,1), zhts(:,1), zati(:,1), zh_i, zh_s, za_i )
1008         !
1009         !                              ! ---------------------- !
1010      ELSEIF( jpl == 1 ) THEN           ! specific case if M = 1 !
1011         !                              ! ---------------------- !
1012         !
1013         CALL ice_var_itd_Nc1c( zhti, zhts, zati, zh_i(:,1), zh_s(:,1), za_i(:,1) )
1014         !
1015         !                              ! ----------------------- !
1016      ELSE                              ! input cat /= output cat !
1017         !                              ! ----------------------- !
1018         
1019         ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays
1020         ALLOCATE( jlmin(idim), jlmax(idim) )
1021
1022         ! --- initialize output fields to 0 --- !
1023         zh_i(1:idim,1:jpl) = 0._wp
1024         zh_s(1:idim,1:jpl) = 0._wp
1025         za_i(1:idim,1:jpl) = 0._wp
1026         !
1027         ! --- fill the categories --- !
1028         !     find where cat-input = cat-output and fill cat-output fields 
1029         jlmax(:) = 0
1030         jlmin(:) = 999
1031         jlfil(:,:) = 0
1032         DO jl1 = 1, jpl
1033            DO jl2 = 1, icat
1034               DO ji = 1, idim
1035                  IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
1036                     ! fill the right category
1037                     zh_i(ji,jl1) = zhti(ji,jl2)
1038                     zh_s(ji,jl1) = zhts(ji,jl2)
1039                     za_i(ji,jl1) = zati(ji,jl2)
1040                     ! record categories that are filled
1041                     jlmax(ji) = MAX( jlmax(ji), jl1 )
1042                     jlmin(ji) = MIN( jlmin(ji), jl1 )
1043                     jlfil(ji,jl1) = jl1
1044                  ENDIF
1045               END DO
1046            END DO
1047         END DO
1048         !
1049         ! --- fill the gaps between categories --- ! 
1050         !     transfer from categories filled at the previous step to the empty ones in between
1051         DO ji = 1, idim
1052            jl1 = jlmin(ji)
1053            jl2 = jlmax(ji)
1054            IF( jl1 > 1 ) THEN
1055               ! fill the lower cat (jl1-1)
1056               za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
1057               zh_i(ji,jl1-1) = hi_mean(jl1-1)
1058               ! remove from cat jl1
1059               za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
1060            ENDIF
1061            IF( jl2 < jpl ) THEN
1062               ! fill the upper cat (jl2+1)
1063               za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
1064               zh_i(ji,jl2+1) = hi_mean(jl2+1)
1065               ! remove from cat jl2
1066               za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
1067            ENDIF
1068         END DO
1069         !
1070         jlfil2(:,:) = jlfil(:,:) 
1071         ! fill categories from low to high
1072         DO jl = 2, jpl-1
1073            DO ji = 1, idim
1074               IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
1075                  ! fill high
1076                  za_i(ji,jl) = ztrans * za_i(ji,jl-1)
1077                  zh_i(ji,jl) = hi_mean(jl)
1078                  jlfil(ji,jl) = jl
1079                  ! remove low
1080                  za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
1081               ENDIF
1082            END DO
1083         END DO
1084         !
1085         ! fill categories from high to low
1086         DO jl = jpl-1, 2, -1
1087            DO ji = 1, idim
1088               IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
1089                  ! fill low
1090                  za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
1091                  zh_i(ji,jl) = hi_mean(jl) 
1092                  jlfil2(ji,jl) = jl
1093                  ! remove high
1094                  za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
1095               ENDIF
1096            END DO
1097         END DO
1098         !
1099         DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
1100         DEALLOCATE( jlmin, jlmax )
1101         !
1102      ENDIF
1103      !
1104   END SUBROUTINE ice_var_itd_NcMc
1105
1106#else
1107   !!----------------------------------------------------------------------
1108   !!   Default option         Dummy module           NO SI3 sea-ice model
1109   !!----------------------------------------------------------------------
1110#endif
1111
1112   !!======================================================================
1113END MODULE icevar
Note: See TracBrowser for help on using the repository browser.