source: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icevar.F90 @ 12371

Last change on this file since 12371 was 12371, checked in by dancopsey, 14 months ago

Add changes to match CICE settings.

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