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

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

source: NEMO/trunk/src/ICE/icevar.F90 @ 10415

Last change on this file since 10415 was 10415, checked in by clem, 6 years ago

correct the advection for ice ponds, make a couple of cosmetic changes and fix the conservation issue for sea ice advection (just a bugged control print)

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