source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icevar.F90 @ 11518

Last change on this file since 11518 was 11518, checked in by clem, 14 months ago

add the final touch to the famous gaston's branch. More precisely, add the possibility to have melt ponds as input file when using bdy

  • Property svn:keywords set to Id
File size: 59.5 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               !
536            END DO
537         END DO
538         !
539      END DO 
540
541      ! to be sure that at_i is the sum of a_i(jl)
542      at_i (:,:) = SUM( a_i (:,:,:), dim=3 )
543      vt_i (:,:) = SUM( v_i (:,:,:), dim=3 )
544!!clem add?
545!      vt_s (:,:) = SUM( v_s (:,:,:), dim=3 )
546!      st_i (:,:) = SUM( sv_i(:,:,:), dim=3 )
547!      et_s(:,:)  = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 )
548!      et_i(:,:)  = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 )
549!!clem
550
551      ! open water = 1 if at_i=0
552      WHERE( at_i(:,:) == 0._wp )   ato_i(:,:) = 1._wp
553      !
554   END SUBROUTINE ice_var_zapsmall
555
556
557   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 )
558      !!-------------------------------------------------------------------
559      !!                   ***  ROUTINE ice_var_zapneg ***
560      !!
561      !! ** Purpose :   Remove negative sea ice fields and correct fluxes
562      !!-------------------------------------------------------------------
563      REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step
564      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
565      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
566      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
567      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_i      ! salt content
568      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
569      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
570      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
571      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
572      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
573      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
574      !
575      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
576      REAL(wp) ::   z1_dt
577      !!-------------------------------------------------------------------
578      !
579      z1_dt = 1._wp / pdt
580      !
581      DO jl = 1, jpl       !==  loop over the categories  ==!
582         !
583         ! make sure a_i=0 where v_i<=0
584         WHERE( pv_i(:,:,:) <= 0._wp )   pa_i(:,:,:) = 0._wp
585
586         !----------------------------------------
587         ! zap ice energy and send it to the ocean
588         !----------------------------------------
589         DO jk = 1, nlay_i
590            DO jj = 1 , jpj
591               DO ji = 1 , jpi
592                  IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
593                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
594                     pe_i(ji,jj,jk,jl) = 0._wp
595                  ENDIF
596               END DO
597            END DO
598         END DO
599         !
600         DO jk = 1, nlay_s
601            DO jj = 1 , jpj
602               DO ji = 1 , jpi
603                  IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
604                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
605                     pe_s(ji,jj,jk,jl) = 0._wp
606                  ENDIF
607               END DO
608            END DO
609         END DO
610         !
611         !-----------------------------------------------------
612         ! zap ice and snow volume, add water and salt to ocean
613         !-----------------------------------------------------
614         DO jj = 1 , jpj
615            DO ji = 1 , jpi
616               IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
617                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
618                  pv_i   (ji,jj,jl) = 0._wp
619               ENDIF
620               IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
621                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
622                  pv_s   (ji,jj,jl) = 0._wp
623               ENDIF
624               IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
625                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
626                  psv_i  (ji,jj,jl) = 0._wp
627               ENDIF
628            END DO
629         END DO
630         !
631      END DO 
632      !
633      WHERE( pato_i(:,:)   < 0._wp )   pato_i(:,:)   = 0._wp
634      WHERE( poa_i (:,:,:) < 0._wp )   poa_i (:,:,:) = 0._wp
635      WHERE( pa_i  (:,:,:) < 0._wp )   pa_i  (:,:,:) = 0._wp
636      WHERE( pa_ip (:,:,:) < 0._wp )   pa_ip (:,:,:) = 0._wp
637      WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
638      !                                                        but it does not change conservation, so keep it this way is ok
639      !
640   END SUBROUTINE ice_var_zapneg
641
642
643   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i )
644      !!-------------------------------------------------------------------
645      !!                   ***  ROUTINE ice_var_roundoff ***
646      !!
647      !! ** Purpose :   Remove negative sea ice values arising from roundoff errors
648      !!-------------------------------------------------------------------
649      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
650      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_i       ! ice volume
651      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_s       ! snw volume
652      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   psv_i      ! salt content
653      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   poa_i      ! age content
654      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
655      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
656      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
657      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
658      !!-------------------------------------------------------------------
659      !
660      WHERE( pa_i (1:npti,:)   < 0._wp .AND. pa_i (1:npti,:)   > -epsi10 )   pa_i (1:npti,:)   = 0._wp   !  a_i must be >= 0
661      WHERE( pv_i (1:npti,:)   < 0._wp .AND. pv_i (1:npti,:)   > -epsi10 )   pv_i (1:npti,:)   = 0._wp   !  v_i must be >= 0
662      WHERE( pv_s (1:npti,:)   < 0._wp .AND. pv_s (1:npti,:)   > -epsi10 )   pv_s (1:npti,:)   = 0._wp   !  v_s must be >= 0
663      WHERE( psv_i(1:npti,:)   < 0._wp .AND. psv_i(1:npti,:)   > -epsi10 )   psv_i(1:npti,:)   = 0._wp   ! sv_i must be >= 0
664      WHERE( poa_i(1:npti,:)   < 0._wp .AND. poa_i(1:npti,:)   > -epsi10 )   poa_i(1:npti,:)   = 0._wp   ! oa_i must be >= 0
665      WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0
666      WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0
667      IF( ln_pnd_H12 ) THEN
668         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0
669         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0
670      ENDIF
671      !
672   END SUBROUTINE ice_var_roundoff
673   
674
675   SUBROUTINE ice_var_bv
676      !!-------------------------------------------------------------------
677      !!                ***  ROUTINE ice_var_bv ***
678      !!
679      !! ** Purpose :   computes mean brine volume (%) in sea ice
680      !!
681      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
682      !!
683      !! References : Vancoppenolle et al., JGR, 2007
684      !!-------------------------------------------------------------------
685      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
686      !!-------------------------------------------------------------------
687      !
688!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
689!!   instead of setting everything to zero as just below
690      bv_i (:,:,:) = 0._wp
691      DO jl = 1, jpl
692         DO jk = 1, nlay_i
693            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
694               bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
695            END WHERE
696         END DO
697      END DO
698      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
699      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
700      END WHERE
701      !
702   END SUBROUTINE ice_var_bv
703
704
705   SUBROUTINE ice_var_enthalpy
706      !!-------------------------------------------------------------------
707      !!                   ***  ROUTINE ice_var_enthalpy ***
708      !!                 
709      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
710      !!
711      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
712      !!-------------------------------------------------------------------
713      INTEGER  ::   ji, jk   ! dummy loop indices
714      REAL(wp) ::   ztmelts  ! local scalar
715      !!-------------------------------------------------------------------
716      !
717      DO jk = 1, nlay_i             ! Sea ice energy of melting
718         DO ji = 1, npti
719            ztmelts      = - rTmlt  * sz_i_1d(ji,jk)
720            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
721                                                                !   (sometimes zdf scheme produces abnormally high temperatures)   
722            e_i_1d(ji,jk) = rhoi * ( rcpi  * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
723               &                   + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
724               &                   - rcp   * ztmelts )
725         END DO
726      END DO
727      DO jk = 1, nlay_s             ! Snow energy of melting
728         DO ji = 1, npti
729            e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
730         END DO
731      END DO
732      !
733   END SUBROUTINE ice_var_enthalpy
734
735   
736   FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
737      !!---------------------------------------------------------------------
738      !!                   ***  ROUTINE ice_var_sshdyn  ***
739      !!                     
740      !! ** Purpose :  compute the equivalent ssh in lead when sea ice is embedded
741      !!
742      !! ** Method  :  ssh_lead = ssh + (Mice + Msnow) / rau0
743      !!
744      !! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
745      !!                Sea ice-ocean coupling using a rescaled vertical coordinate z*,
746      !!                Ocean Modelling, Volume 24, Issues 1-2, 2008
747      !!----------------------------------------------------------------------
748      !
749      ! input
750      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh            !: ssh [m]
751      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass    !: mass of snow and ice at current  ice time step [Kg/m2]
752      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b  !: mass of snow and ice at previous ice time step [Kg/m2]
753      !
754      ! output
755      REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn  ! equivalent ssh in lead [m]
756      !
757      ! temporary
758      REAL(wp) :: zintn, zintb                     ! time interpolation weights []
759      REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload  ! snow and ice load [m]
760      !
761      ! compute ice load used to define the equivalent ssh in lead
762      IF( ln_ice_embd ) THEN
763         !                                           
764         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
765         !                                               = (1/nn_fsbc)^2 * {SUM[n]        , n=0,nn_fsbc-1}
766         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
767         !
768         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   *    {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
769         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
770         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
771         !
772         zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0
773         !
774      ELSE
775         zsnwiceload(:,:) = 0.0_wp
776      ENDIF
777      ! compute equivalent ssh in lead
778      ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:)
779      !
780   END FUNCTION ice_var_sshdyn
781
782   
783   !!-------------------------------------------------------------------
784   !!                ***  INTERFACE ice_var_itd   ***
785   !!
786   !! ** Purpose :  converting N-cat ice to jpl ice categories
787   !!-------------------------------------------------------------------
788   SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,                       ph_i, ph_s, pa_i, &
789      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
790      !!-------------------------------------------------------------------
791      !! ** Purpose :  converting 1-cat ice to 1 ice category
792      !!-------------------------------------------------------------------
793      REAL(wp), DIMENSION(:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
794      REAL(wp), DIMENSION(:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
795      REAL(wp), DIMENSION(:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds
796      REAL(wp), DIMENSION(:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds
797      !!-------------------------------------------------------------------
798      ! == thickness and concentration == !
799      ph_i(:) = phti(:)
800      ph_s(:) = phts(:)
801      pa_i(:) = pati(:)
802      !
803      ! == temperature and salinity and ponds == !
804      pt_i (:) = ptmi (:)
805      pt_s (:) = ptms (:)
806      pt_su(:) = ptmsu(:)
807      ps_i (:) = psmi (:)
808      pa_ip(:) = patip(:)
809      ph_ip(:) = phtip(:)
810     
811   END SUBROUTINE ice_var_itd_1c1c
812
813   SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,                       ph_i, ph_s, pa_i, &
814      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
815      !!-------------------------------------------------------------------
816      !! ** Purpose :  converting N-cat ice to 1 ice category
817      !!-------------------------------------------------------------------
818      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
819      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
820      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds
821      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds
822      !
823      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z1_ai, z1_vi, z1_vs
824      !
825      INTEGER ::   idim 
826      !!-------------------------------------------------------------------
827      !
828      idim = SIZE( phti, 1 )
829      !
830      ! == thickness and concentration == !
831      ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim) )
832      !
833      pa_i(:) = SUM( pati(:,:), dim=2 )
834
835      WHERE( ( pa_i(:) ) /= 0._wp )   ;   z1_ai(:) = 1._wp / pa_i(:)
836      ELSEWHERE                       ;   z1_ai(:) = 0._wp
837      END WHERE
838
839      ph_i(:) = SUM( phti(:,:) * pati(:,:), dim=2 ) * z1_ai(:)
840      ph_s(:) = SUM( phts(:,:) * pati(:,:), dim=2 ) * z1_ai(:)
841      !
842      ! == temperature and salinity == !
843      WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp )   ;   z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) )
844      ELSEWHERE                                 ;   z1_vi(:) = 0._wp
845      END WHERE
846      WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp )   ;   z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) )
847      ELSEWHERE                                 ;   z1_vs(:) = 0._wp
848      END WHERE
849      pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
850      pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
851      pt_su(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:)
852      ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
853
854      ! == ponds == !
855      pa_ip(:) = SUM( patip(:,:), dim=2 )
856      WHERE( pa_ip(:) /= 0._wp )   ;   ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
857      ELSEWHERE                    ;   ph_ip(:) = 0._wp
858      END WHERE
859      !
860      DEALLOCATE( z1_ai, z1_vi, z1_vs )
861      !
862   END SUBROUTINE ice_var_itd_Nc1c
863   
864   SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, &
865      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
866      !!-------------------------------------------------------------------
867      !!
868      !! ** Purpose :  converting 1-cat ice to jpl ice categories
869      !!
870      !!                  ice thickness distribution follows a gaussian law
871      !!               around the concentration of the most likely ice thickness
872      !!                           (similar as iceistate.F90)
873      !!
874      !! ** Method:   Iterative procedure
875      !!               
876      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
877      !!
878      !!               2) Check whether the distribution conserves area and volume, positivity and
879      !!                  category boundaries
880      !!             
881      !!               3) If not (input ice is too thin), the last category is empty and
882      !!                  the number of categories is reduced (jpl-1)
883      !!
884      !!               4) Iterate until ok (SUM(itest(:) = 4)
885      !!
886      !! ** Arguments : phti: 1-cat ice thickness
887      !!                phts: 1-cat snow depth
888      !!                pati: 1-cat ice concentration
889      !!
890      !! ** Output    : jpl-cat
891      !!
892      !!  (Example of application: BDY forcings when input are cell averaged) 
893      !!-------------------------------------------------------------------
894      REAL(wp), DIMENSION(:),   INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
895      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
896      REAL(wp), DIMENSION(:)  , INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds
897      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds
898      !
899      INTEGER , DIMENSION(4) ::   itest
900      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zfra
901      INTEGER  ::   ji, jk, jl
902      INTEGER  ::   idim, i_fill, jl0 
903      REAL(wp) ::   zarg, zV, zconv, zdh, zdv
904      !!-------------------------------------------------------------------
905      !
906      ! == thickness and concentration == !
907      ! distribution over the jpl ice categories
908      !    a gaussian distribution for ice concentration is used
909      !    then we check whether the distribution fullfills
910      !    volume and area conservation, positivity and ice categories bounds
911      idim = SIZE( phti , 1 )
912      !
913      ph_i(1:idim,1:jpl) = 0._wp
914      ph_s(1:idim,1:jpl) = 0._wp
915      pa_i(1:idim,1:jpl) = 0._wp
916      !
917      DO ji = 1, idim
918         !
919         IF( phti(ji) > 0._wp ) THEN
920            !
921            ! find which category (jl0) the input ice thickness falls into
922            jl0 = jpl
923            DO jl = 1, jpl
924               IF ( ( phti(ji) >= hi_max(jl-1) ) .AND. ( phti(ji) < hi_max(jl) ) ) THEN
925                  jl0 = jl
926                  CYCLE
927               ENDIF
928            END DO
929            !
930            itest(:) = 0
931            i_fill   = jpl + 1                                            !------------------------------------
932            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
933               !                                                          !------------------------------------
934               i_fill = i_fill - 1
935               !
936               ph_i(ji,1:jpl) = 0._wp
937               pa_i(ji,1:jpl) = 0._wp
938               itest(:)       = 0     
939               !
940               IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
941                  ph_i(ji,1) = phti(ji)
942                  pa_i(ji,1) = pati (ji)
943               ELSE                         !-- case ice is thicker: fill categories >1
944                  ! thickness
945                  DO jl = 1, i_fill - 1
946                     ph_i(ji,jl) = hi_mean(jl)
947                  END DO
948                  !
949                  ! concentration
950                  pa_i(ji,jl0) = pati(ji) / SQRT(REAL(jpl))
951                  DO jl = 1, i_fill - 1
952                     IF ( jl /= jl0 ) THEN
953                        zarg        = ( ph_i(ji,jl) - phti(ji) ) / ( phti(ji) * 0.5_wp )
954                        pa_i(ji,jl) =   pa_i (ji,jl0) * EXP(-zarg**2)
955                     ENDIF
956                  END DO
957                  !
958                  ! last category
959                  pa_i(ji,i_fill) = pati(ji) - SUM( pa_i(ji,1:i_fill-1) )
960                  zV = SUM( pa_i(ji,1:i_fill-1) * ph_i(ji,1:i_fill-1) )
961                  ph_i(ji,i_fill) = ( phti(ji) * pati(ji) - zV ) / MAX( pa_i(ji,i_fill), epsi10 ) 
962                  !
963                  ! correction if concentration of upper cat is greater than lower cat
964                  !    (it should be a gaussian around jl0 but sometimes it is not)
965                  IF ( jl0 /= jpl ) THEN
966                     DO jl = jpl, jl0+1, -1
967                        IF ( pa_i(ji,jl) > pa_i(ji,jl-1) ) THEN
968                           zdv = ph_i(ji,jl) * pa_i(ji,jl)
969                           ph_i(ji,jl    ) = 0._wp
970                           pa_i (ji,jl    ) = 0._wp
971                           pa_i (ji,1:jl-1) = pa_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * phti(ji), epsi10 )
972                        END IF
973                     END DO
974                  ENDIF
975                  !
976               ENDIF
977               !
978               ! Compatibility tests
979               zconv = ABS( pati(ji) - SUM( pa_i(ji,1:jpl) ) ) 
980               IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation
981               !
982               zconv = ABS( phti(ji)*pati(ji) - SUM( pa_i(ji,1:jpl)*ph_i(ji,1:jpl) ) )
983               IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation
984               !
985               IF ( ph_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                   ! Test 3: thickness of the last category is in-bounds ?
986               !
987               itest(4) = 1
988               DO jl = 1, i_fill
989                  IF ( pa_i(ji,jl) < 0._wp ) itest(4) = 0                                  ! Test 4: positivity of ice concentrations
990               END DO
991               !                                         !----------------------------
992            END DO                                       ! end iteration on categories
993            !                                            !----------------------------
994         ENDIF
995      END DO
996
997      ! Add Snow in each category where pa_i is not 0
998      DO jl = 1, jpl
999         DO ji = 1, idim
1000            IF( pa_i(ji,jl) > 0._wp ) THEN
1001               ph_s(ji,jl) = ph_i(ji,jl) * ( phts(ji) / phti(ji) )
1002               ! In case snow load is in excess that would lead to transformation from snow to ice
1003               ! Then, transfer the snow excess into the ice (different from icethd_dh)
1004               zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rau0 ) * ph_i(ji,jl) ) * r1_rau0 ) 
1005               ! recompute h_i, h_s avoiding out of bounds values
1006               ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh )
1007               ph_s(ji,jl) = MAX( 0._wp, ph_s(ji,jl) - zdh * rhoi * r1_rhos )
1008            ENDIF
1009         END DO
1010      END DO
1011      !
1012      ! == temperature and salinity == !
1013      DO jl = 1, jpl
1014         pt_i (:,jl) = ptmi (:)
1015         pt_s (:,jl) = ptms (:)
1016         pt_su(:,jl) = ptmsu(:)
1017         ps_i (:,jl) = psmi (:)
1018         ps_i (:,jl) = psmi (:)         
1019      END DO
1020      !
1021      ! == ponds == !
1022      ALLOCATE( zfra(idim) )
1023      ! keep the same pond fraction atip/ati for each category
1024      WHERE( pati(:) /= 0._wp )   ;   zfra(:) = patip(:) / pati(:)
1025      ELSEWHERE                   ;   zfra(:) = 0._wp
1026      END WHERE
1027      DO jl = 1, jpl
1028         pa_ip(:,jl) = zfra(:) * pa_i(:,jl)
1029      END DO
1030      ! keep the same v_ip/v_i ratio for each category
1031      WHERE( ( phti(:) * pati(:) ) /= 0._wp )   ;   zfra(:) = ( phtip(:) * patip(:) ) / ( phti(:) * pati(:) )
1032      ELSEWHERE                                 ;   zfra(:) = 0._wp
1033      END WHERE
1034      DO jl = 1, jpl
1035         WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1036         ELSEWHERE                       ;   ph_ip(:,jl) = 0._wp
1037         END WHERE
1038      END DO
1039      DEALLOCATE( zfra )
1040      !
1041   END SUBROUTINE ice_var_itd_1cMc
1042
1043   SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, &
1044      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
1045      !!-------------------------------------------------------------------
1046      !!
1047      !! ** Purpose :  converting N-cat ice to jpl ice categories
1048      !!
1049      !!                  ice thickness distribution follows a gaussian law
1050      !!               around the concentration of the most likely ice thickness
1051      !!                           (similar as iceistate.F90)
1052      !!
1053      !! ** Method:   Iterative procedure
1054      !!               
1055      !!               1) Fill ice cat that correspond to input thicknesses
1056      !!                  Find the lowest(jlmin) and highest(jlmax) cat that are filled
1057      !!
1058      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1
1059       !!                   by removing 25% ice area from jlmin and jlmax (resp.)
1060      !!             
1061      !!               3) Expand the filling to the empty cat between jlmin and jlmax
1062      !!                   by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
1063      !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
1064      !!
1065      !! ** Arguments : phti: N-cat ice thickness
1066      !!                phts: N-cat snow depth
1067      !!                pati: N-cat ice concentration
1068      !!
1069      !! ** Output    : jpl-cat
1070      !!
1071      !!  (Example of application: BDY forcings when inputs have N-cat /= jpl) 
1072      !!-------------------------------------------------------------------
1073      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
1074      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
1075      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds
1076      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds
1077      !
1078      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   jlfil, jlfil2
1079      INTEGER , ALLOCATABLE, DIMENSION(:)   ::   jlmax, jlmin
1080      REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   z1_ai, z1_vi, z1_vs, ztmp, zfra
1081      !
1082      REAL(wp), PARAMETER ::   ztrans = 0.25_wp
1083      INTEGER  ::   ji, jl, jl1, jl2
1084      INTEGER  ::   idim, icat 
1085      !!-------------------------------------------------------------------
1086      !
1087      idim = SIZE( phti, 1 )
1088      icat = SIZE( phti, 2 )
1089      !
1090      ! == thickness and concentration == !
1091      !                                 ! ---------------------- !
1092      IF( icat == jpl ) THEN            ! input cat = output cat !
1093         !                              ! ---------------------- !
1094         ph_i(:,:) = phti(:,:)
1095         ph_s(:,:) = phts(:,:)
1096         pa_i(:,:) = pati(:,:)
1097         !
1098         ! == temperature and salinity and ponds == !
1099         pt_i (:,:) = ptmi (:,:)
1100         pt_s (:,:) = ptms (:,:)
1101         pt_su(:,:) = ptmsu(:,:)
1102         ps_i (:,:) = psmi (:,:)
1103         pa_ip(:,:) = patip(:,:)
1104         ph_ip(:,:) = phtip(:,:)
1105         !                              ! ---------------------- !
1106      ELSEIF( icat == 1 ) THEN          ! input cat = 1          !
1107         !                              ! ---------------------- !
1108         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), &
1109            &                    ph_i(:,:), ph_s(:,:), pa_i (:,:), &
1110            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), &
1111            &                    pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:)  )
1112         !                              ! ---------------------- !
1113      ELSEIF( jpl == 1 ) THEN           ! output cat = 1         !
1114         !                              ! ---------------------- !
1115         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), &
1116            &                    ph_i(:,1), ph_s(:,1), pa_i (:,1), &
1117            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), &
1118            &                    pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1)  )
1119         !                              ! ----------------------- !
1120      ELSE                              ! input cat /= output cat !
1121         !                              ! ----------------------- !
1122         
1123         ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays
1124         ALLOCATE( jlmin(idim), jlmax(idim) )
1125
1126         ! --- initialize output fields to 0 --- !
1127         ph_i(1:idim,1:jpl) = 0._wp
1128         ph_s(1:idim,1:jpl) = 0._wp
1129         pa_i(1:idim,1:jpl) = 0._wp
1130         !
1131         ! --- fill the categories --- !
1132         !     find where cat-input = cat-output and fill cat-output fields 
1133         jlmax(:) = 0
1134         jlmin(:) = 999
1135         jlfil(:,:) = 0
1136         DO jl1 = 1, jpl
1137            DO jl2 = 1, icat
1138               DO ji = 1, idim
1139                  IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN
1140                     ! fill the right category
1141                     ph_i(ji,jl1) = phti(ji,jl2)
1142                     ph_s(ji,jl1) = phts(ji,jl2)
1143                     pa_i(ji,jl1) = pati(ji,jl2)
1144                     ! record categories that are filled
1145                     jlmax(ji) = MAX( jlmax(ji), jl1 )
1146                     jlmin(ji) = MIN( jlmin(ji), jl1 )
1147                     jlfil(ji,jl1) = jl1
1148                  ENDIF
1149               END DO
1150            END DO
1151         END DO
1152         !
1153         ! --- fill the gaps between categories --- ! 
1154         !     transfer from categories filled at the previous step to the empty ones in between
1155         DO ji = 1, idim
1156            jl1 = jlmin(ji)
1157            jl2 = jlmax(ji)
1158            IF( jl1 > 1 ) THEN
1159               ! fill the lower cat (jl1-1)
1160               pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1)
1161               ph_i(ji,jl1-1) = hi_mean(jl1-1)
1162               ! remove from cat jl1
1163               pa_i(ji,jl1  ) = ( 1._wp - ztrans ) * pa_i(ji,jl1)
1164            ENDIF
1165            IF( jl2 < jpl ) THEN
1166               ! fill the upper cat (jl2+1)
1167               pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2)
1168               ph_i(ji,jl2+1) = hi_mean(jl2+1)
1169               ! remove from cat jl2
1170               pa_i(ji,jl2  ) = ( 1._wp - ztrans ) * pa_i(ji,jl2)
1171            ENDIF
1172         END DO
1173         !
1174         jlfil2(:,:) = jlfil(:,:) 
1175         ! fill categories from low to high
1176         DO jl = 2, jpl-1
1177            DO ji = 1, idim
1178               IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
1179                  ! fill high
1180                  pa_i(ji,jl) = ztrans * pa_i(ji,jl-1)
1181                  ph_i(ji,jl) = hi_mean(jl)
1182                  jlfil(ji,jl) = jl
1183                  ! remove low
1184                  pa_i(ji,jl-1) = ( 1._wp - ztrans ) * pa_i(ji,jl-1)
1185               ENDIF
1186            END DO
1187         END DO
1188         !
1189         ! fill categories from high to low
1190         DO jl = jpl-1, 2, -1
1191            DO ji = 1, idim
1192               IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
1193                  ! fill low
1194                  pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1)
1195                  ph_i(ji,jl) = hi_mean(jl) 
1196                  jlfil2(ji,jl) = jl
1197                  ! remove high
1198                  pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1)
1199               ENDIF
1200            END DO
1201         END DO
1202         !
1203         DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
1204         DEALLOCATE( jlmin, jlmax )
1205         !
1206         ! == temperature and salinity == !
1207         !
1208         ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) )
1209         !
1210         WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp )               ;   z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 )
1211         ELSEWHERE                                               ;   z1_ai(:) = 0._wp
1212         END WHERE
1213         WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp )   ;   z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 )
1214         ELSEWHERE                                               ;   z1_vi(:) = 0._wp
1215         END WHERE
1216         WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp )   ;   z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 )
1217         ELSEWHERE                                               ;   z1_vs(:) = 0._wp
1218         END WHERE
1219         !
1220         ! fill all the categories with the same value
1221         ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
1222         DO jl = 1, jpl
1223            pt_i (:,jl) = ztmp(:)
1224         END DO
1225         ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
1226         DO jl = 1, jpl
1227            pt_s (:,jl) = ztmp(:)
1228         END DO
1229         ztmp(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:)
1230         DO jl = 1, jpl
1231            pt_su(:,jl) = ztmp(:)
1232         END DO
1233         ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
1234         DO jl = 1, jpl
1235            ps_i (:,jl) = ztmp(:)
1236         END DO
1237         !
1238         DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp )
1239         !
1240         ! == ponds == !
1241         ALLOCATE( zfra(idim) )
1242         ! keep the same pond fraction atip/ati for each category
1243         WHERE( SUM( pati(:,:), dim=2 ) /= 0._wp )   ;   zfra(:) = SUM( patip(:,:), dim=2 ) / SUM( pati(:,:), dim=2 )
1244         ELSEWHERE                                   ;   zfra(:) = 0._wp
1245         END WHERE
1246         DO jl = 1, jpl
1247            pa_ip(:,jl) = zfra(:) * pa_i(:,jl)
1248         END DO
1249         ! keep the same v_ip/v_i ratio for each category
1250         WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp )
1251            zfra(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 )
1252         ELSEWHERE
1253            zfra(:) = 0._wp
1254         END WHERE
1255         DO jl = 1, jpl
1256            WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1257            ELSEWHERE                       ;   ph_ip(:,jl) = 0._wp
1258            END WHERE
1259         END DO
1260         DEALLOCATE( zfra )
1261         !
1262      ENDIF
1263      !
1264   END SUBROUTINE ice_var_itd_NcMc
1265
1266#else
1267   !!----------------------------------------------------------------------
1268   !!   Default option         Dummy module           NO SI3 sea-ice model
1269   !!----------------------------------------------------------------------
1270#endif
1271
1272   !!======================================================================
1273END MODULE icevar
Note: See TracBrowser for help on using the repository browser.