New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icevar.F90 in NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE – NEMO

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

Last change on this file since 11507 was 11507, checked in by clem, 5 years ago

add a new functionality for bdy ice by allowing the user to have input files for ice/snw temperature and salinity instead of setting a constant value. Also, input files can have any number of categories

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