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/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE – NEMO

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

Last change on this file since 12369 was 12369, checked in by dancopsey, 5 years ago

Add print statements

File size: 49.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
321!!gm Question: Remove the option 3 ?  How many years since it last use ?
322
323      SELECT CASE ( nn_icesal )
324      !
325      !               !---------------------------------------!
326      CASE( 1 )       !  constant salinity in time and space  !
327         !            !---------------------------------------!
328         sz_i(:,:,:,:) = rn_icesal
329         s_i (:,:,:)   = rn_icesal
330         !
331         !            !---------------------------------------------!
332      CASE( 2 )       !  time varying salinity with linear profile  !
333         !            !---------------------------------------------!
334         !
335         ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
336         !
337         DO jl = 1, jpl
338            DO jk = 1, nlay_i
339               sz_i(:,:,jk,jl)  = s_i(:,:,jl)
340            END DO
341         END DO
342
343         !                                      ! Slope of the linear profile
344         WHERE( h_i(:,:,:) > epsi20 )   ;   z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:)
345         ELSEWHERE                      ;   z_slope_s(:,:,:) = 0._wp
346         END WHERE
347         !
348         z1_dS = 1._wp / ( zsi1 - zsi0 )
349         DO jl = 1, jpl
350            DO jj = 1, jpj
351               DO ji = 1, jpi
352                  zalpha(ji,jj,jl) = MAX(  0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp )  )
353                  !                             ! force a constant profile when SSS too low (Baltic Sea)
354                  IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) )   zalpha(ji,jj,jl) = 0._wp 
355               END DO
356            END DO
357         END DO
358         
359         !
360         ! Computation of the profile
361         DO jl = 1, jpl
362            DO jk = 1, nlay_i
363               DO jj = 1, jpj
364                  DO ji = 1, jpi
365                     !                          ! linear profile with 0 surface value
366                     zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
367                     zs  = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl)     ! weighting the profile
368                     sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
369
370                     IF ( jl == 1 .AND. jj == 26 .AND. ji == 42 ) THEN
371                       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)
372                     ENDIF
373                  END DO
374               END DO
375            END DO
376         END DO
377         !
378         DEALLOCATE( z_slope_s , zalpha )
379         !
380         !            !-------------------------------------------!
381      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
382         !            !-------------------------------------------!                                   (mean = 2.30)
383         !
384         s_i(:,:,:) = 2.30_wp
385!!gm Remark: if we keep the case 3, then compute an store one for all time-step
386!!           a array  S_prof(1:nlay_i) containing the calculation and just do:
387!         DO jk = 1, nlay_i
388!            sz_i(:,:,jk,:) = S_prof(jk)
389!         END DO
390!!gm end
391         !
392         DO jl = 1, jpl
393            DO jk = 1, nlay_i
394               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
395               sz_i(:,:,jk,jl) =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
396            END DO
397         END DO
398         !
399      END SELECT
400      !
401   END SUBROUTINE ice_var_salprof
402
403
404   SUBROUTINE ice_var_salprof1d
405      !!-------------------------------------------------------------------
406      !!                  ***  ROUTINE ice_var_salprof1d  ***
407      !!
408      !! ** Purpose :   1d computation of the sea ice salinity profile
409      !!                Works with 1d vectors and is used by thermodynamic modules
410      !!-------------------------------------------------------------------
411      INTEGER  ::   ji, jk    ! dummy loop indices
412      REAL(wp) ::   zargtemp, zsal, z1_dS   ! local scalars
413      REAL(wp) ::   zs, zs0              !   -      -
414      !
415      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z_slope_s, zalpha   !
416      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
417      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
418      !!-------------------------------------------------------------------
419      !
420      SELECT CASE ( nn_icesal )
421      !
422      !               !---------------------------------------!
423      CASE( 1 )       !  constant salinity in time and space  !
424         !            !---------------------------------------!
425         sz_i_1d(1:npti,:) = rn_icesal
426         !
427         !            !---------------------------------------------!
428      CASE( 2 )       !  time varying salinity with linear profile  !
429         !            !---------------------------------------------!
430         !
431         ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
432         !
433         !                                      ! Slope of the linear profile
434         WHERE( h_i_1d(1:npti) > epsi20 )   ;   z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
435         ELSEWHERE                          ;   z_slope_s(1:npti) = 0._wp
436         END WHERE
437         
438         z1_dS = 1._wp / ( zsi1 - zsi0 )
439         DO ji = 1, npti
440            zalpha(ji) = MAX(  0._wp , MIN(  ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp  )  )
441            !                             ! force a constant profile when SSS too low (Baltic Sea)
442            IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) )   zalpha(ji) = 0._wp
443         END DO
444         !
445         ! Computation of the profile
446         DO jk = 1, nlay_i
447            DO ji = 1, npti
448               !                          ! linear profile with 0 surface value
449               zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
450               zs  = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
451               sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
452            END DO
453         END DO
454         !
455         DEALLOCATE( z_slope_s, zalpha )
456
457         !            !-------------------------------------------!
458      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
459         !            !-------------------------------------------!                                   (mean = 2.30)
460         !
461         s_i_1d(1:npti) = 2.30_wp
462         !
463!!gm cf remark in ice_var_salprof routine, CASE( 3 )
464         DO jk = 1, nlay_i
465            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
466            zsal =  1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
467            DO ji = 1, npti
468               sz_i_1d(ji,jk) = zsal
469            END DO
470         END DO
471         !
472      END SELECT
473      !
474   END SUBROUTINE ice_var_salprof1d
475
476
477   SUBROUTINE ice_var_zapsmall
478      !!-------------------------------------------------------------------
479      !!                   ***  ROUTINE ice_var_zapsmall ***
480      !!
481      !! ** Purpose :   Remove too small sea ice areas and correct fluxes
482      !!-------------------------------------------------------------------
483      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
484      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch
485      !!-------------------------------------------------------------------
486      !
487      DO jl = 1, jpl       !==  loop over the categories  ==!
488         !
489         WHERE( a_i(:,:,jl) > epsi10 )   ;   h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
490         ELSEWHERE                       ;   h_i(:,:,jl) = 0._wp
491         END WHERE
492         !
493         WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 )   ;   zswitch(:,:) = 0._wp
494         ELSEWHERE                                                                           ;   zswitch(:,:) = 1._wp
495         END WHERE
496         !
497         !-----------------------------------------------------------------
498         ! Zap ice energy and use ocean heat to melt ice
499         !-----------------------------------------------------------------
500         DO jk = 1, nlay_i
501            DO jj = 1 , jpj
502               DO ji = 1 , jpi
503                  ! update exchanges with ocean
504                  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
505                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
506                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
507               END DO
508            END DO
509         END DO
510         !
511         DO jk = 1, nlay_s
512            DO jj = 1 , jpj
513               DO ji = 1 , jpi
514                  ! update exchanges with ocean
515                  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
516                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj)
517                  t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
518
519                  If (to_print_2d(ji,jj) == 10 .AND. jk == 1 .AND. jl == 1) THEN
520                    write(numout,*)'ice_var_glo2eqv 1: t_s, e_s = ',t_s(ji,jj,jk,jl), ' ', e_s (ji,jj,jk,jl)
521                  END IF
522               END DO
523            END DO
524         END DO
525         !
526         !-----------------------------------------------------------------
527         ! zap ice and snow volume, add water and salt to ocean
528         !-----------------------------------------------------------------
529         DO jj = 1 , jpj
530            DO ji = 1 , jpi
531               ! update exchanges with ocean
532               sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoi * r1_rdtice
533               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoi * r1_rdtice
534               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhos * r1_rdtice
535               !
536               a_i  (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
537               v_i  (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
538               v_s  (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
539               t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
540               oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
541               sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
542               !
543               h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
544               h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
545               !
546               a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
547               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
548               !
549            END DO
550         END DO
551         !
552      END DO 
553
554      ! to be sure that at_i is the sum of a_i(jl)
555      at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
556      vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
557
558      ! open water = 1 if at_i=0
559      WHERE( at_i(:,:) == 0._wp )   ato_i(:,:) = 1._wp
560      !
561   END SUBROUTINE ice_var_zapsmall
562
563
564   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 )
565      !!-------------------------------------------------------------------
566      !!                   ***  ROUTINE ice_var_zapneg ***
567      !!
568      !! ** Purpose :   Remove negative sea ice fields and correct fluxes
569      !!-------------------------------------------------------------------
570      REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step
571      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
572      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
573      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
574      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_i      ! salt content
575      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
576      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
577      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
578      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
579      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
580      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
581      !
582      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
583      REAL(wp) ::   z1_dt
584      !!-------------------------------------------------------------------
585      !
586      z1_dt = 1._wp / pdt
587      !
588      DO jl = 1, jpl       !==  loop over the categories  ==!
589         !
590         ! make sure a_i=0 where v_i<=0
591         WHERE( pv_i(:,:,:) <= 0._wp )   pa_i(:,:,:) = 0._wp
592
593         !----------------------------------------
594         ! zap ice energy and send it to the ocean
595         !----------------------------------------
596         DO jk = 1, nlay_i
597            DO jj = 1 , jpj
598               DO ji = 1 , jpi
599                  IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
600                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
601                     pe_i(ji,jj,jk,jl) = 0._wp
602                  ENDIF
603               END DO
604            END DO
605         END DO
606         !
607         DO jk = 1, nlay_s
608            DO jj = 1 , jpj
609               DO ji = 1 , jpi
610                  IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
611                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
612                     pe_s(ji,jj,jk,jl) = 0._wp
613                  ENDIF
614               END DO
615            END DO
616         END DO
617         !
618         !-----------------------------------------------------
619         ! zap ice and snow volume, add water and salt to ocean
620         !-----------------------------------------------------
621         DO jj = 1 , jpj
622            DO ji = 1 , jpi
623               IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
624                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
625                  pv_i   (ji,jj,jl) = 0._wp
626               ENDIF
627               IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
628                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
629                  pv_s   (ji,jj,jl) = 0._wp
630               ENDIF
631               IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
632                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
633                  psv_i  (ji,jj,jl) = 0._wp
634               ENDIF
635            END DO
636         END DO
637         !
638      END DO 
639      !
640      WHERE( pato_i(:,:)   < 0._wp )   pato_i(:,:)   = 0._wp
641      WHERE( poa_i (:,:,:) < 0._wp )   poa_i (:,:,:) = 0._wp
642      WHERE( pa_i  (:,:,:) < 0._wp )   pa_i  (:,:,:) = 0._wp
643      WHERE( pa_ip (:,:,:) < 0._wp )   pa_ip (:,:,:) = 0._wp
644      WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
645      !                                                        but it does not change conservation, so keep it this way is ok
646      !
647   END SUBROUTINE ice_var_zapneg
648
649
650   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i )
651      !!-------------------------------------------------------------------
652      !!                   ***  ROUTINE ice_var_roundoff ***
653      !!
654      !! ** Purpose :   Remove negative sea ice values arising from roundoff errors
655      !!-------------------------------------------------------------------
656      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
657      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_i       ! ice volume
658      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_s       ! snw volume
659      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   psv_i      ! salt content
660      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   poa_i      ! age content
661      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
662      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
663      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
664      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
665      !!-------------------------------------------------------------------
666      !
667      WHERE( pa_i (1:npti,:)   < 0._wp .AND. pa_i (1:npti,:)   > -epsi10 )   pa_i (1:npti,:)   = 0._wp   !  a_i must be >= 0
668      WHERE( pv_i (1:npti,:)   < 0._wp .AND. pv_i (1:npti,:)   > -epsi10 )   pv_i (1:npti,:)   = 0._wp   !  v_i must be >= 0
669      WHERE( pv_s (1:npti,:)   < 0._wp .AND. pv_s (1:npti,:)   > -epsi10 )   pv_s (1:npti,:)   = 0._wp   !  v_s must be >= 0
670      WHERE( psv_i(1:npti,:)   < 0._wp .AND. psv_i(1:npti,:)   > -epsi10 )   psv_i(1:npti,:)   = 0._wp   ! sv_i must be >= 0
671      WHERE( poa_i(1:npti,:)   < 0._wp .AND. poa_i(1:npti,:)   > -epsi10 )   poa_i(1:npti,:)   = 0._wp   ! oa_i must be >= 0
672      WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0
673      WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0
674      IF ( ln_pnd_H12 ) THEN
675         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0
676         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0
677      ENDIF
678      !
679   END SUBROUTINE ice_var_roundoff
680   
681   
682   SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
683      !!-------------------------------------------------------------------
684      !!                ***  ROUTINE ice_var_itd   ***
685      !!
686      !! ** Purpose :  converting 1-cat ice to multiple ice categories
687      !!
688      !!                  ice thickness distribution follows a gaussian law
689      !!               around the concentration of the most likely ice thickness
690      !!                           (similar as iceistate.F90)
691      !!
692      !! ** Method:   Iterative procedure
693      !!               
694      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
695      !!
696      !!               2) Check whether the distribution conserves area and volume, positivity and
697      !!                  category boundaries
698      !!             
699      !!               3) If not (input ice is too thin), the last category is empty and
700      !!                  the number of categories is reduced (jpl-1)
701      !!
702      !!               4) Iterate until ok (SUM(itest(:) = 4)
703      !!
704      !! ** Arguments : zhti: 1-cat ice thickness
705      !!                zhts: 1-cat snow depth
706      !!                zati: 1-cat ice concentration
707      !!
708      !! ** Output    : jpl-cat
709      !!
710      !!  (Example of application: BDY forcings when input are cell averaged) 
711      !!-------------------------------------------------------------------
712      INTEGER  :: ji, jk, jl             ! dummy loop indices
713      INTEGER  :: idim, i_fill, jl0 
714      REAL(wp) :: zarg, zV, zconv, zdh, zdv
715      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input  ice/snow variables
716      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
717      INTEGER , DIMENSION(4)                  ::   itest
718      !!-------------------------------------------------------------------
719      !
720      ! ----------------------------------------
721      ! distribution over the jpl ice categories
722      ! ----------------------------------------
723      ! a gaussian distribution for ice concentration is used
724      ! then we check whether the distribution fullfills
725      ! volume and area conservation, positivity and ice categories bounds
726      idim = SIZE( zhti , 1 )
727      zh_i(1:idim,1:jpl) = 0._wp
728      zh_s(1:idim,1:jpl) = 0._wp
729      za_i(1:idim,1:jpl) = 0._wp
730      !
731      DO ji = 1, idim
732         !
733         IF( zhti(ji) > 0._wp ) THEN
734            !
735            ! find which category (jl0) the input ice thickness falls into
736            jl0 = jpl
737            DO jl = 1, jpl
738               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
739                  jl0 = jl
740                  CYCLE
741               ENDIF
742            END DO
743            !
744            itest(:) = 0
745            i_fill   = jpl + 1                                            !------------------------------------
746            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
747               !                                                          !------------------------------------
748               i_fill = i_fill - 1
749               !
750               zh_i(ji,1:jpl) = 0._wp
751               za_i(ji,1:jpl) = 0._wp
752               itest(:)       = 0     
753               !
754               IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
755                  zh_i(ji,1) = zhti(ji)
756                  za_i (ji,1) = zati (ji)
757               ELSE                         !-- case ice is thicker: fill categories >1
758                  ! thickness
759                  DO jl = 1, i_fill - 1
760                     zh_i(ji,jl) = hi_mean(jl)
761                  END DO
762                  !
763                  ! concentration
764                  za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
765                  DO jl = 1, i_fill - 1
766                     IF ( jl /= jl0 ) THEN
767                        zarg        = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
768                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
769                     ENDIF
770                  END DO
771                  !
772                  ! last category
773                  za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
774                  zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
775                  zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
776                  !
777                  ! correction if concentration of upper cat is greater than lower cat
778                  !    (it should be a gaussian around jl0 but sometimes it is not)
779                  IF ( jl0 /= jpl ) THEN
780                     DO jl = jpl, jl0+1, -1
781                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
782                           zdv = zh_i(ji,jl) * za_i(ji,jl)
783                           zh_i(ji,jl    ) = 0._wp
784                           za_i (ji,jl    ) = 0._wp
785                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
786                        END IF
787                     END DO
788                  ENDIF
789                  !
790               ENDIF
791               !
792               ! Compatibility tests
793               zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) ) 
794               IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation
795               !
796               zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )
797               IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation
798               !
799               IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ?
800               !
801               itest(4) = 1
802               DO jl = 1, i_fill
803                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0                                ! Test 4: positivity of ice concentrations
804               END DO
805               !                                         !----------------------------
806            END DO                                       ! end iteration on categories
807            !                                            !----------------------------
808         ENDIF
809      END DO
810
811      ! Add Snow in each category where za_i is not 0
812      DO jl = 1, jpl
813         DO ji = 1, idim
814            IF( za_i(ji,jl) > 0._wp ) THEN
815               zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
816               ! In case snow load is in excess that would lead to transformation from snow to ice
817               ! Then, transfer the snow excess into the ice (different from icethd_dh)
818               zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 ) 
819               ! recompute h_i, h_s avoiding out of bounds values
820               zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
821               zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos )
822            ENDIF
823         END DO
824      END DO
825      !
826   END SUBROUTINE ice_var_itd
827
828
829   SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )
830      !!-------------------------------------------------------------------
831      !!                ***  ROUTINE ice_var_itd2   ***
832      !!
833      !! ** Purpose :  converting N-cat ice to jpl ice categories
834      !!
835      !!                  ice thickness distribution follows a gaussian law
836      !!               around the concentration of the most likely ice thickness
837      !!                           (similar as iceistate.F90)
838      !!
839      !! ** Method:   Iterative procedure
840      !!               
841      !!               1) Fill ice cat that correspond to input thicknesses
842      !!                  Find the lowest(jlmin) and highest(jlmax) cat that are filled
843      !!
844      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1
845       !!                   by removing 25% ice area from jlmin and jlmax (resp.)
846      !!             
847      !!               3) Expand the filling to the empty cat between jlmin and jlmax
848      !!                   by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
849      !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
850      !!
851      !! ** Arguments : zhti: N-cat ice thickness
852      !!                zhts: N-cat snow depth
853      !!                zati: N-cat ice concentration
854      !!
855      !! ** Output    : jpl-cat
856      !!
857      !!  (Example of application: BDY forcings when inputs have N-cat /= jpl) 
858      !!-------------------------------------------------------------------
859      INTEGER  ::   ji, jl, jl1, jl2             ! dummy loop indices
860      INTEGER  ::   idim, icat 
861      REAL(wp), PARAMETER ::   ztrans = 0.25_wp
862      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
863      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
864      INTEGER , DIMENSION(:,:), ALLOCATABLE   ::   jlfil, jlfil2
865      INTEGER , DIMENSION(:)  , ALLOCATABLE   ::   jlmax, jlmin
866      !!-------------------------------------------------------------------
867      !
868      idim = SIZE( zhti, 1 )
869      icat = SIZE( zhti, 2 )
870      !
871      ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays
872      ALLOCATE( jlmin(idim), jlmax(idim) )
873
874      ! --- initialize output fields to 0 --- !
875      zh_i(1:idim,1:jpl) = 0._wp
876      zh_s(1:idim,1:jpl) = 0._wp
877      za_i(1:idim,1:jpl) = 0._wp
878      !
879      ! --- fill the categories --- !
880      !     find where cat-input = cat-output and fill cat-output fields 
881      jlmax(:) = 0
882      jlmin(:) = 999
883      jlfil(:,:) = 0
884      DO jl1 = 1, jpl
885         DO jl2 = 1, icat
886            DO ji = 1, idim
887               IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
888                  ! fill the right category
889                  zh_i(ji,jl1) = zhti(ji,jl2)
890                  zh_s(ji,jl1) = zhts(ji,jl2)
891                  za_i(ji,jl1) = zati(ji,jl2)
892                  ! record categories that are filled
893                  jlmax(ji) = MAX( jlmax(ji), jl1 )
894                  jlmin(ji) = MIN( jlmin(ji), jl1 )
895                  jlfil(ji,jl1) = jl1
896               ENDIF
897            END DO
898         END DO
899      END DO
900      !
901      ! --- fill the gaps between categories --- ! 
902      !     transfer from categories filled at the previous step to the empty ones in between
903      DO ji = 1, idim
904         jl1 = jlmin(ji)
905         jl2 = jlmax(ji)
906         IF( jl1 > 1 ) THEN
907            ! fill the lower cat (jl1-1)
908            za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
909            zh_i(ji,jl1-1) = hi_mean(jl1-1)
910            ! remove from cat jl1
911            za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
912         ENDIF
913         IF( jl2 < jpl ) THEN
914            ! fill the upper cat (jl2+1)
915            za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
916            zh_i(ji,jl2+1) = hi_mean(jl2+1)
917            ! remove from cat jl2
918            za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
919         ENDIF
920      END DO
921      !
922      jlfil2(:,:) = jlfil(:,:) 
923      ! fill categories from low to high
924      DO jl = 2, jpl-1
925         DO ji = 1, idim
926            IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
927               ! fill high
928               za_i(ji,jl) = ztrans * za_i(ji,jl-1)
929               zh_i(ji,jl) = hi_mean(jl)
930               jlfil(ji,jl) = jl
931               ! remove low
932               za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
933            ENDIF
934         END DO
935      END DO
936      !
937      ! fill categories from high to low
938      DO jl = jpl-1, 2, -1
939         DO ji = 1, idim
940            IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
941               ! fill low
942               za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
943               zh_i(ji,jl) = hi_mean(jl) 
944               jlfil2(ji,jl) = jl
945               ! remove high
946               za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
947            ENDIF
948         END DO
949      END DO
950      !
951      DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
952      DEALLOCATE( jlmin, jlmax )
953      !
954   END SUBROUTINE ice_var_itd2
955
956
957   SUBROUTINE ice_var_bv
958      !!-------------------------------------------------------------------
959      !!                ***  ROUTINE ice_var_bv ***
960      !!
961      !! ** Purpose :   computes mean brine volume (%) in sea ice
962      !!
963      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
964      !!
965      !! References : Vancoppenolle et al., JGR, 2007
966      !!-------------------------------------------------------------------
967      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
968      !!-------------------------------------------------------------------
969      !
970!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
971!!   instead of setting everything to zero as just below
972      bv_i (:,:,:) = 0._wp
973      DO jl = 1, jpl
974         DO jk = 1, nlay_i
975            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
976               bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
977            END WHERE
978         END DO
979      END DO
980      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
981      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
982      END WHERE
983      !
984   END SUBROUTINE ice_var_bv
985
986
987   SUBROUTINE ice_var_enthalpy
988      !!-------------------------------------------------------------------
989      !!                   ***  ROUTINE ice_var_enthalpy ***
990      !!                 
991      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
992      !!
993      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
994      !!-------------------------------------------------------------------
995      INTEGER  ::   ji, jk   ! dummy loop indices
996      REAL(wp) ::   ztmelts  ! local scalar
997      !!-------------------------------------------------------------------
998      !
999      DO jk = 1, nlay_i             ! Sea ice energy of melting
1000         DO ji = 1, npti
1001            ztmelts      = - rTmlt  * sz_i_1d(ji,jk)
1002            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
1003                                                                !   (sometimes zdf scheme produces abnormally high temperatures)   
1004            e_i_1d(ji,jk) = rhoi * ( rcpi  * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
1005               &                   + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
1006               &                   - rcp   * ztmelts )
1007         END DO
1008      END DO
1009      DO jk = 1, nlay_s             ! Snow energy of melting
1010         DO ji = 1, npti
1011            e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
1012         END DO
1013      END DO
1014      !
1015   END SUBROUTINE ice_var_enthalpy
1016
1017   FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
1018      !!---------------------------------------------------------------------
1019      !!                   ***  ROUTINE ice_var_sshdyn  ***
1020      !!                     
1021      !! ** Purpose :  compute the equivalent ssh in lead when sea ice is embedded
1022      !!
1023      !! ** Method  :  ssh_lead = ssh + (Mice + Msnow) / rau0
1024      !!
1025      !! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
1026      !!                Sea ice-ocean coupling using a rescaled vertical coordinate z*,
1027      !!                Ocean Modelling, Volume 24, Issues 1-2, 2008
1028      !!----------------------------------------------------------------------
1029      !
1030      ! input
1031      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh            !: ssh [m]
1032      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass    !: mass of snow and ice at current  ice time step [Kg/m2]
1033      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b  !: mass of snow and ice at previous ice time step [Kg/m2]
1034      !
1035      ! output
1036      REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn  ! equivalent ssh in lead [m]
1037      !
1038      ! temporary
1039      REAL(wp) :: zintn, zintb                     ! time interpolation weights []
1040      REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload  ! snow and ice load [m]
1041      !
1042      ! compute ice load used to define the equivalent ssh in lead
1043      IF( ln_ice_embd ) THEN
1044         !                                           
1045         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
1046         !                                               = (1/nn_fsbc)^2 * {SUM[n]        , n=0,nn_fsbc-1}
1047         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
1048         !
1049         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   *    {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
1050         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
1051         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
1052         !
1053         zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0
1054         !
1055      ELSE
1056         zsnwiceload(:,:) = 0.0_wp
1057      ENDIF
1058      ! compute equivalent ssh in lead
1059      ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:)
1060      !
1061   END FUNCTION ice_var_sshdyn
1062
1063
1064#else
1065   !!----------------------------------------------------------------------
1066   !!   Default option         Dummy module           NO SI3 sea-ice model
1067   !!----------------------------------------------------------------------
1068#endif
1069
1070   !!======================================================================
1071END MODULE icevar
Note: See TracBrowser for help on using the repository browser.