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/releases/release-4.0/src/ICE – NEMO

source: NEMO/releases/release-4.0/src/ICE/icevar.F90 @ 10944

Last change on this file since 10944 was 10944, checked in by clem, 4 years ago

finish repairing the code when only ice dynamics is activated

  • Property svn:keywords set to Id
File size: 45.9 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( pdt, 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      REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step
546      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
547      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
548      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
549      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_i      ! salt content
550      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
551      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
552      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
553      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
554      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
555      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
556      !
557      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
558      REAL(wp) ::   z1_dt
559      !!-------------------------------------------------------------------
560      !
561      z1_dt = 1._wp / pdt
562      !
563      DO jl = 1, jpl       !==  loop over the categories  ==!
564         !
565         ! make sure a_i=0 where v_i<=0
566         WHERE( pv_i(:,:,:) <= 0._wp )   pa_i(:,:,:) = 0._wp
567
568         !----------------------------------------
569         ! zap ice energy and send it to the ocean
570         !----------------------------------------
571         DO jk = 1, nlay_i
572            DO jj = 1 , jpj
573               DO ji = 1 , jpi
574                  IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
575                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
576                     pe_i(ji,jj,jk,jl) = 0._wp
577                  ENDIF
578               END DO
579            END DO
580         END DO
581         !
582         DO jk = 1, nlay_s
583            DO jj = 1 , jpj
584               DO ji = 1 , jpi
585                  IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
586                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
587                     pe_s(ji,jj,jk,jl) = 0._wp
588                  ENDIF
589               END DO
590            END DO
591         END DO
592         !
593         !-----------------------------------------------------
594         ! zap ice and snow volume, add water and salt to ocean
595         !-----------------------------------------------------
596         DO jj = 1 , jpj
597            DO ji = 1 , jpi
598               IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
599                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
600                  pv_i   (ji,jj,jl) = 0._wp
601               ENDIF
602               IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
603                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
604                  pv_s   (ji,jj,jl) = 0._wp
605               ENDIF
606               IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
607                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
608                  psv_i  (ji,jj,jl) = 0._wp
609               ENDIF
610            END DO
611         END DO
612         !
613      END DO 
614      !
615      WHERE( pato_i(:,:)   < 0._wp )   pato_i(:,:)   = 0._wp
616      WHERE( poa_i (:,:,:) < 0._wp )   poa_i (:,:,:) = 0._wp
617      WHERE( pa_i  (:,:,:) < 0._wp )   pa_i  (:,:,:) = 0._wp
618      WHERE( pa_ip (:,:,:) < 0._wp )   pa_ip (:,:,:) = 0._wp
619      WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
620      !                                                        but it does not change conservation, so keep it this way is ok
621      !
622   END SUBROUTINE ice_var_zapneg
623
624   
625   SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
626      !!-------------------------------------------------------------------
627      !!                ***  ROUTINE ice_var_itd   ***
628      !!
629      !! ** Purpose :  converting 1-cat ice to multiple ice categories
630      !!
631      !!                  ice thickness distribution follows a gaussian law
632      !!               around the concentration of the most likely ice thickness
633      !!                           (similar as iceistate.F90)
634      !!
635      !! ** Method:   Iterative procedure
636      !!               
637      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
638      !!
639      !!               2) Check whether the distribution conserves area and volume, positivity and
640      !!                  category boundaries
641      !!             
642      !!               3) If not (input ice is too thin), the last category is empty and
643      !!                  the number of categories is reduced (jpl-1)
644      !!
645      !!               4) Iterate until ok (SUM(itest(:) = 4)
646      !!
647      !! ** Arguments : zhti: 1-cat ice thickness
648      !!                zhts: 1-cat snow depth
649      !!                zati: 1-cat ice concentration
650      !!
651      !! ** Output    : jpl-cat
652      !!
653      !!  (Example of application: BDY forcings when input are cell averaged) 
654      !!-------------------------------------------------------------------
655      INTEGER  :: ji, jk, jl             ! dummy loop indices
656      INTEGER  :: idim, i_fill, jl0 
657      REAL(wp) :: zarg, zV, zconv, zdh, zdv
658      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
659      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i ! output ice/snow variables
660      INTEGER , DIMENSION(4)                  ::   itest
661      !!-------------------------------------------------------------------
662      !
663      ! ----------------------------------------
664      ! distribution over the jpl ice categories
665      ! ----------------------------------------
666      ! a gaussian distribution for ice concentration is used
667      ! then we check whether the distribution fullfills
668      ! volume and area conservation, positivity and ice categories bounds
669      idim = SIZE( zhti , 1 )
670      zh_i(1:idim,1:jpl) = 0._wp
671      zh_s(1:idim,1:jpl) = 0._wp
672      za_i(1:idim,1:jpl) = 0._wp
673      !
674      DO ji = 1, idim
675         !
676         IF( zhti(ji) > 0._wp ) THEN
677            !
678            ! find which category (jl0) the input ice thickness falls into
679            jl0 = jpl
680            DO jl = 1, jpl
681               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
682                  jl0 = jl
683                  CYCLE
684               ENDIF
685            END DO
686            !
687            itest(:) = 0
688            i_fill   = jpl + 1                                            !------------------------------------
689            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
690               !                                                          !------------------------------------
691               i_fill = i_fill - 1
692               !
693               zh_i(ji,1:jpl) = 0._wp
694               za_i(ji,1:jpl) = 0._wp
695               itest(:)       = 0     
696               !
697               IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
698                  zh_i(ji,1) = zhti(ji)
699                  za_i (ji,1) = zati (ji)
700               ELSE                         !-- case ice is thicker: fill categories >1
701                  ! thickness
702                  DO jl = 1, i_fill - 1
703                     zh_i(ji,jl) = hi_mean(jl)
704                  END DO
705                  !
706                  ! concentration
707                  za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
708                  DO jl = 1, i_fill - 1
709                     IF ( jl /= jl0 ) THEN
710                        zarg        = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
711                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
712                     ENDIF
713                  END DO
714                  !
715                  ! last category
716                  za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
717                  zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
718                  zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
719                  !
720                  ! correction if concentration of upper cat is greater than lower cat
721                  !    (it should be a gaussian around jl0 but sometimes it is not)
722                  IF ( jl0 /= jpl ) THEN
723                     DO jl = jpl, jl0+1, -1
724                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
725                           zdv = zh_i(ji,jl) * za_i(ji,jl)
726                           zh_i(ji,jl    ) = 0._wp
727                           za_i (ji,jl    ) = 0._wp
728                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
729                        END IF
730                     END DO
731                  ENDIF
732                  !
733               ENDIF
734               !
735               ! Compatibility tests
736               zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) ) 
737               IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation
738               !
739               zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )
740               IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation
741               !
742               IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ?
743               !
744               itest(4) = 1
745               DO jl = 1, i_fill
746                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0                                ! Test 4: positivity of ice concentrations
747               END DO
748               !                                         !----------------------------
749            END DO                                       ! end iteration on categories
750            !                                            !----------------------------
751         ENDIF
752      END DO
753
754      ! Add Snow in each category where za_i is not 0
755      DO jl = 1, jpl
756         DO ji = 1, idim
757            IF( za_i(ji,jl) > 0._wp ) THEN
758               zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
759               ! In case snow load is in excess that would lead to transformation from snow to ice
760               ! Then, transfer the snow excess into the ice (different from icethd_dh)
761               zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 ) 
762               ! recompute h_i, h_s avoiding out of bounds values
763               zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
764               zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos )
765            ENDIF
766         END DO
767      END DO
768      !
769   END SUBROUTINE ice_var_itd
770
771
772   SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )
773      !!-------------------------------------------------------------------
774      !!                ***  ROUTINE ice_var_itd2   ***
775      !!
776      !! ** Purpose :  converting N-cat ice to jpl ice categories
777      !!
778      !!                  ice thickness distribution follows a gaussian law
779      !!               around the concentration of the most likely ice thickness
780      !!                           (similar as iceistate.F90)
781      !!
782      !! ** Method:   Iterative procedure
783      !!               
784      !!               1) Fill ice cat that correspond to input thicknesses
785      !!                  Find the lowest(jlmin) and highest(jlmax) cat that are filled
786      !!
787      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1
788      !!                   by removing 25% ice area from jlmin and jlmax (resp.)
789      !!             
790      !!               3) Expand the filling to the empty cat between jlmin and jlmax
791      !!                   by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
792      !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
793      !!
794      !! ** Arguments : zhti: N-cat ice thickness
795      !!                zhts: N-cat snow depth
796      !!                zati: N-cat ice concentration
797      !!
798      !! ** Output    : jpl-cat
799      !!
800      !!  (Example of application: BDY forcings when inputs have N-cat /= jpl) 
801      !!-------------------------------------------------------------------
802      INTEGER  ::   ji, jl, jl1, jl2             ! dummy loop indices
803      INTEGER  ::   idim, icat 
804      REAL(wp), PARAMETER ::   ztrans = 0.25_wp
805      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
806      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
807      INTEGER , DIMENSION(:,:), ALLOCATABLE   ::   jlfil, jlfil2
808      INTEGER , DIMENSION(:)  , ALLOCATABLE   ::   jlmax, jlmin
809      !!-------------------------------------------------------------------
810      !
811      idim = SIZE( zhti, 1 )
812      icat = SIZE( zhti, 2 )
813      !
814      ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays
815      ALLOCATE( jlmin(idim), jlmax(idim) )
816
817      ! --- initialize output fields to 0 --- !
818      zh_i(1:idim,1:jpl) = 0._wp
819      zh_s(1:idim,1:jpl) = 0._wp
820      za_i(1:idim,1:jpl) = 0._wp
821      !
822      ! --- fill the categories --- !
823      !     find where cat-input = cat-output and fill cat-output fields 
824      jlmax(:) = 0
825      jlmin(:) = 999
826      jlfil(:,:) = 0
827      DO jl1 = 1, jpl
828         DO jl2 = 1, icat
829            DO ji = 1, idim
830               IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
831                  ! fill the right category
832                  zh_i(ji,jl1) = zhti(ji,jl2)
833                  zh_s(ji,jl1) = zhts(ji,jl2)
834                  za_i(ji,jl1) = zati(ji,jl2)
835                  ! record categories that are filled
836                  jlmax(ji) = MAX( jlmax(ji), jl1 )
837                  jlmin(ji) = MIN( jlmin(ji), jl1 )
838                  jlfil(ji,jl1) = jl1
839               ENDIF
840            END DO
841         END DO
842      END DO
843      !
844      ! --- fill the gaps between categories --- ! 
845      !     transfer from categories filled at the previous step to the empty ones in between
846      DO ji = 1, idim
847         jl1 = jlmin(ji)
848         jl2 = jlmax(ji)
849         IF( jl1 > 1 ) THEN
850            ! fill the lower cat (jl1-1)
851            za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
852            zh_i(ji,jl1-1) = hi_mean(jl1-1)
853            ! remove from cat jl1
854            za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
855         ENDIF
856         IF( jl2 < jpl ) THEN
857            ! fill the upper cat (jl2+1)
858            za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
859            zh_i(ji,jl2+1) = hi_mean(jl2+1)
860            ! remove from cat jl2
861            za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
862         ENDIF
863      END DO
864      !
865      jlfil2(:,:) = jlfil(:,:) 
866      ! fill categories from low to high
867      DO jl = 2, jpl-1
868         DO ji = 1, idim
869            IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
870               ! fill high
871               za_i(ji,jl) = ztrans * za_i(ji,jl-1)
872               zh_i(ji,jl) = hi_mean(jl)
873               jlfil(ji,jl) = jl
874               ! remove low
875               za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
876            ENDIF
877         END DO
878      END DO
879      !
880      ! fill categories from high to low
881      DO jl = jpl-1, 2, -1
882         DO ji = 1, idim
883            IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
884               ! fill low
885               za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
886               zh_i(ji,jl) = hi_mean(jl) 
887               jlfil2(ji,jl) = jl
888               ! remove high
889               za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
890            ENDIF
891         END DO
892      END DO
893      !
894      DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
895      DEALLOCATE( jlmin, jlmax )
896      !
897   END SUBROUTINE ice_var_itd2
898
899
900   SUBROUTINE ice_var_bv
901      !!-------------------------------------------------------------------
902      !!                ***  ROUTINE ice_var_bv ***
903      !!
904      !! ** Purpose :   computes mean brine volume (%) in sea ice
905      !!
906      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
907      !!
908      !! References : Vancoppenolle et al., JGR, 2007
909      !!-------------------------------------------------------------------
910      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
911      !!-------------------------------------------------------------------
912      !
913!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
914!!   instead of setting everything to zero as just below
915      bv_i (:,:,:) = 0._wp
916      DO jl = 1, jpl
917         DO jk = 1, nlay_i
918            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
919               bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
920            END WHERE
921         END DO
922      END DO
923      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
924      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
925      END WHERE
926      !
927   END SUBROUTINE ice_var_bv
928
929
930   SUBROUTINE ice_var_enthalpy
931      !!-------------------------------------------------------------------
932      !!                   ***  ROUTINE ice_var_enthalpy ***
933      !!                 
934      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
935      !!
936      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
937      !!-------------------------------------------------------------------
938      INTEGER  ::   ji, jk   ! dummy loop indices
939      REAL(wp) ::   ztmelts  ! local scalar
940      !!-------------------------------------------------------------------
941      !
942      DO jk = 1, nlay_i             ! Sea ice energy of melting
943         DO ji = 1, npti
944            ztmelts      = - rTmlt  * sz_i_1d(ji,jk)
945            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
946                                                                !   (sometimes zdf scheme produces abnormally high temperatures)   
947            e_i_1d(ji,jk) = rhoi * ( rcpi  * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
948               &                   + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
949               &                   - rcp   * ztmelts )
950         END DO
951      END DO
952      DO jk = 1, nlay_s             ! Snow energy of melting
953         DO ji = 1, npti
954            e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
955         END DO
956      END DO
957      !
958   END SUBROUTINE ice_var_enthalpy
959
960   FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
961      !!---------------------------------------------------------------------
962      !!                   ***  ROUTINE ice_var_sshdyn  ***
963      !!                     
964      !! ** Purpose :  compute the equivalent ssh in lead when sea ice is embedded
965      !!
966      !! ** Method  :  ssh_lead = ssh + (Mice + Msnow) / rau0
967      !!
968      !! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
969      !!                Sea ice-ocean coupling using a rescaled vertical coordinate z*,
970      !!                Ocean Modelling, Volume 24, Issues 1-2, 2008
971      !!----------------------------------------------------------------------
972      !
973      ! input
974      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh            !: ssh [m]
975      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass    !: mass of snow and ice at current  ice time step [Kg/m2]
976      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b  !: mass of snow and ice at previous ice time step [Kg/m2]
977      !
978      ! output
979      REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn  ! equivalent ssh in lead [m]
980      !
981      ! temporary
982      REAL(wp) :: zintn, zintb                     ! time interpolation weights []
983      REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload  ! snow and ice load [m]
984      !
985      ! compute ice load used to define the equivalent ssh in lead
986      IF( ln_ice_embd ) THEN
987         !                                           
988         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
989         !                                               = (1/nn_fsbc)^2 * {SUM[n]        , n=0,nn_fsbc-1}
990         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
991         !
992         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   *    {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
993         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
994         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
995         !
996         zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0
997         !
998      ELSE
999         zsnwiceload(:,:) = 0.0_wp
1000      ENDIF
1001      ! compute equivalent ssh in lead
1002      ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:)
1003      !
1004   END FUNCTION ice_var_sshdyn
1005
1006
1007#else
1008   !!----------------------------------------------------------------------
1009   !!   Default option         Dummy module           NO SI3 sea-ice model
1010   !!----------------------------------------------------------------------
1011#endif
1012
1013   !!======================================================================
1014END MODULE icevar
Note: See TracBrowser for help on using the repository browser.