source: NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icevar.F90 @ 12443

Last change on this file since 12443 was 12443, checked in by davestorkey, 9 months ago

2020/KERNEL-03_Storkey_Coward_RK3_stage2: More variable renaming:
atfp → rn_atfp (use namelist parameter everywhere)
rdtbt → rDt_e
nn_baro → nn_e
rn_scal_load → rn_load
rau0 → rho0

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