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_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE – NEMO

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

Last change on this file since 11334 was 11334, checked in by flemarie, 5 years ago

Follow-on to the ABL implementation (see ticket #2131)

  • Modification of diawri.F90 to output ABL variables with IOIPSL (i.e. without key_iomput)
  • Modification of ice_var_agg in icevar.F90 to allow the use of tm_su in ABL and bulk
  • Error handling in case ln_abl = .true. and ABL sources were not compiled
  • ABL now works with SI3 considering an average over ice categories
  • Update reference namelist to avoid runtime warnings due to nam_tide
  • Sanity checks with ORCA2 for the following configurations :

1 - ABL src + ln_blk = .true.
2 - ABL src + ln_abl = .true.
3 - no ABL src + ln_blk = .true.

All configurations are MPP reproducible and configurations 1 and 3 results are identical

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