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

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

prepare code to read temperature and salinity for bdy

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