source: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icevar.F90 @ 11413

Last change on this file since 11413 was 11413, checked in by gsamson, 21 months ago

dev_r11265_ABL : see #2131

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