source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icevar.F90 @ 11381

Last change on this file since 11381 was 11381, checked in by girrmann, 2 years ago

dev_r10984_HPC-13 : repair faulty commit [11380], see #2308

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