source: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icevar.F90 @ 12636

Last change on this file since 12636 was 12636, checked in by dancopsey, 8 months ago

Add:

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