source: NEMO/trunk/src/ICE/icevar.F90 @ 9725

Last change on this file since 9725 was 9725, checked in by clem, 3 years ago

solve issue with uninitialized field (tm_si, purely diagnostic) in agrif configuration (and probably orca2_ice also). This issue only occurs for a certain mpp decomposition for a reason I do not know. 32 procs was failing while 16 procs was working.

File size: 39.3 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_itd       : convert 1-cat to jpl-cat
47   !!   ice_var_itd2      : convert N-cat to jpl-cat
48   !!   ice_var_bv        : brine volume
49   !!   ice_var_enthalpy  : compute ice and snow enthalpies from temperature
50   !!----------------------------------------------------------------------
51   USE dom_oce        ! ocean space and time domain
52   USE phycst         ! physical constants (ocean directory)
53   USE sbc_oce , ONLY : sss_m
54   USE ice            ! sea-ice: variables
55   USE ice1D          ! sea-ice: thermodynamics variables
56   !
57   USE in_out_manager ! I/O manager
58   USE lib_mpp        ! MPP library
59   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
60
61   IMPLICIT NONE
62   PRIVATE
63
64   PUBLIC   ice_var_agg         
65   PUBLIC   ice_var_glo2eqv     
66   PUBLIC   ice_var_eqv2glo     
67   PUBLIC   ice_var_salprof     
68   PUBLIC   ice_var_salprof1d   
69   PUBLIC   ice_var_zapsmall
70   PUBLIC   ice_var_itd
71   PUBLIC   ice_var_itd2
72   PUBLIC   ice_var_bv           
73   PUBLIC   ice_var_enthalpy           
74
75   !!----------------------------------------------------------------------
76   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
77   !! $Id: icevar.F90 8422 2017-08-08 13:58:05Z clem $
78   !! Software governed by the CeCILL licence     (./LICENSE)
79   !!----------------------------------------------------------------------
80CONTAINS
81
82   SUBROUTINE ice_var_agg( kn )
83      !!-------------------------------------------------------------------
84      !!                ***  ROUTINE ice_var_agg  ***
85      !!
86      !! ** Purpose :   aggregates ice-thickness-category variables to
87      !!              all-ice variables, i.e. it turns VGLO into VAGG
88      !!-------------------------------------------------------------------
89      INTEGER, INTENT( in ) ::   kn     ! =1 state variables only
90      !                                 ! >1 state variables + others
91      !
92      INTEGER ::   ji, jj, jk, jl   ! dummy loop indices
93      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z1_at_i, z1_vt_i, z1_vt_s
94      !!-------------------------------------------------------------------
95      !
96      !                                      ! integrated values
97      vt_i(:,:) =       SUM( v_i(:,:,:)           , dim=3 )
98      vt_s(:,:) =       SUM( v_s(:,:,:)           , dim=3 )
99      at_i(:,:) =       SUM( a_i(:,:,:)           , dim=3 )
100      et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 )
101      et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
102      !
103      at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds
104      vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
105      !
106      ato_i(:,:) = 1._wp - at_i(:,:)         ! open water fraction 
107
108      ! The following fields are calculated for diagnostics and outputs only
109      ! ==> Do not use them for other purposes
110      IF( kn > 1 ) THEN
111         !
112         ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) )
113         WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:)
114         ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp
115         END WHERE
116         WHERE( vt_i(:,:) > epsi20 )   ;   z1_vt_i(:,:) = 1._wp / vt_i(:,:)
117         ELSEWHERE                     ;   z1_vt_i(:,:) = 0._wp
118         END WHERE
119         WHERE( vt_s(:,:) > epsi20 )   ;   z1_vt_s(:,:) = 1._wp / vt_s(:,:)
120         ELSEWHERE                     ;   z1_vt_s(:,:) = 0._wp
121         END WHERE
122         !
123         !                          ! mean ice/snow thickness
124         hm_i(:,:) = vt_i(:,:) * z1_at_i(:,:)
125         hm_s(:,:) = vt_s(:,:) * z1_at_i(:,:)
126         !         
127         !                          ! mean temperature (K), salinity and age
128         tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
129         tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
130         om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(:,:)
131         sm_i (:,:) = SUM( sv_i(:,:,:)              , dim=3 ) * z1_vt_i(:,:)
132         !
133         tm_i(:,:) = 0._wp
134         tm_s(:,:) = 0._wp
135         DO jl = 1, jpl
136            DO jk = 1, nlay_i
137               tm_i(:,:) = tm_i(:,:) + r1_nlay_i * t_i (:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
138            END DO
139            DO jk = 1, nlay_s
140               tm_s(:,:) = tm_s(:,:) + r1_nlay_s * t_s (:,:,jk,jl) * v_s(:,:,jl) * z1_vt_s(:,:)
141            END DO
142         END DO
143         !
144         !                           ! put rt0 where there is no ice
145         WHERE( at_i(:,:)<=epsi20 )
146            tm_su(:,:) = rt0
147            tm_si(:,:) = rt0
148            tm_i (:,:) = rt0
149            tm_s (:,:) = rt0
150         END WHERE
151
152         DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s )
153      ENDIF
154      !
155   END SUBROUTINE ice_var_agg
156
157
158   SUBROUTINE ice_var_glo2eqv
159      !!-------------------------------------------------------------------
160      !!                ***  ROUTINE ice_var_glo2eqv ***
161      !!
162      !! ** Purpose :   computes equivalent variables as function of 
163      !!              global variables, i.e. it turns VGLO into VEQV
164      !!-------------------------------------------------------------------
165      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
166      REAL(wp) ::   ze_i             ! local scalars
167      REAL(wp) ::   ze_s, ztmelts, zbbb, zccc       !   -      -
168      REAL(wp) ::   zhmax, z1_zhmax                 !   -      -
169      REAL(wp) ::   zlay_i, zlay_s                  !   -      -
170      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i
171      !!-------------------------------------------------------------------
172
173!!gm Question 2:  It is possible to define existence of sea-ice in a common way between
174!!                ice area and ice volume ?
175!!                the idea is to be able to define one for all at the begining of this routine
176!!                a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
177
178      !---------------------------------------------------------------
179      ! Ice thickness, snow thickness, ice salinity, ice age and ponds
180      !---------------------------------------------------------------
181      !                                            !--- inverse of the ice area
182      WHERE( a_i(:,:,:) > epsi20 )   ;   z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
183      ELSEWHERE                      ;   z1_a_i(:,:,:) = 0._wp
184      END WHERE
185      !
186      WHERE( v_i(:,:,:) > epsi20 )   ;   z1_v_i(:,:,:) = 1._wp / v_i(:,:,:)
187      ELSEWHERE                      ;   z1_v_i(:,:,:) = 0._wp
188      END WHERE
189      !                                           !--- ice thickness
190      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
191
192      zhmax    =          hi_max(jpl)
193      z1_zhmax =  1._wp / hi_max(jpl)               
194      WHERE( h_i(:,:,jpl) > zhmax )   ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area
195         h_i  (:,:,jpl) = zhmax
196         a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 
197         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)
198      END WHERE
199      !                                           !--- snow thickness
200      h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:)
201      !                                           !--- ice age     
202      o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:)
203      !                                           !--- pond fraction and thickness     
204      a_ip_frac(:,:,:) = a_ip(:,:,:) * z1_a_i(:,:,:)
205      WHERE( a_ip_frac(:,:,:) > epsi20 )   ;   h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:)
206      ELSEWHERE                            ;   h_ip(:,:,:) = 0._wp
207      END WHERE
208      !
209      !                                           !---  salinity (with a minimum value imposed everywhere)     
210      IF( nn_icesal == 2 ) THEN
211         WHERE( v_i(:,:,:) > epsi20 )   ;   s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) * z1_v_i(:,:,:) ) )
212         ELSEWHERE                      ;   s_i(:,:,:) = rn_simin
213         END WHERE
214      ENDIF
215      CALL ice_var_salprof   ! salinity profile
216
217      !-------------------
218      ! Ice temperature   [K]   (with a minimum value (rt0 - 100.))
219      !-------------------
220      zlay_i   = REAL( nlay_i , wp )    ! number of layers
221      DO jl = 1, jpl
222         DO jk = 1, nlay_i
223            DO jj = 1, jpj
224               DO ji = 1, jpi
225                  IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area
226                     !
227                     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]
228                     ztmelts          = - sz_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C]
229                     ! Conversion q(S,T) -> T (second order equation)
230                     zbbb             = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus
231                     zccc             = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) )
232                     t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_cpic , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts
233                     !
234                  ELSE                                !--- no ice
235                     t_i(ji,jj,jk,jl) = rt0
236                  ENDIF
237               END DO
238            END DO
239         END DO
240      END DO
241
242      !--------------------
243      ! Snow temperature   [K]   (with a minimum value (rt0 - 100.))
244      !--------------------
245      zlay_s = REAL( nlay_s , wp )
246      DO jk = 1, nlay_s
247         WHERE( v_s(:,:,:) > epsi20 )        !--- icy area
248            t_s(:,:,jk,:) = rt0 + MAX( -100._wp ,  &
249                 &                MIN( r1_cpic * ( -r1_rhosn * ( e_s(:,:,jk,:) / v_s(:,:,:) * zlay_s ) + lfus ) , 0._wp ) )
250         ELSEWHERE                           !--- no ice
251            t_s(:,:,jk,:) = rt0
252         END WHERE
253      END DO
254      !
255      ! integrated values
256      vt_i (:,:) = SUM( v_i, dim=3 )
257      vt_s (:,:) = SUM( v_s, dim=3 )
258      at_i (:,:) = SUM( a_i, dim=3 )
259      !
260   END SUBROUTINE ice_var_glo2eqv
261
262
263   SUBROUTINE ice_var_eqv2glo
264      !!-------------------------------------------------------------------
265      !!                ***  ROUTINE ice_var_eqv2glo ***
266      !!
267      !! ** Purpose :   computes global variables as function of
268      !!              equivalent variables,  i.e. it turns VEQV into VGLO
269      !!-------------------------------------------------------------------
270      !
271      v_i (:,:,:) = h_i (:,:,:) * a_i (:,:,:)
272      v_s (:,:,:) = h_s (:,:,:) * a_i (:,:,:)
273      sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:)
274      v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
275      !
276   END SUBROUTINE ice_var_eqv2glo
277
278
279   SUBROUTINE ice_var_salprof
280      !!-------------------------------------------------------------------
281      !!                ***  ROUTINE ice_var_salprof ***
282      !!
283      !! ** Purpose :   computes salinity profile in function of bulk salinity     
284      !!
285      !! ** Method  : If bulk salinity greater than zsi1,
286      !!              the profile is assumed to be constant (S_inf)
287      !!              If bulk salinity lower than zsi0,
288      !!              the profile is linear with 0 at the surface (S_zero)
289      !!              If it is between zsi0 and zsi1, it is a
290      !!              alpha-weighted linear combination of s_inf and s_zero
291      !!
292      !! ** References : Vancoppenolle et al., 2007
293      !!-------------------------------------------------------------------
294      INTEGER  ::   ji, jj, jk, jl   ! dummy loop index
295      REAL(wp) ::   zsal, z1_dS
296      REAL(wp) ::   zargtemp , zs0, zs
297      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z_slope_s, zalpha    ! case 2 only
298      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
299      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
300      !!-------------------------------------------------------------------
301
302!!gm Question: Remove the option 3 ?  How many years since it last use ?
303
304      SELECT CASE ( nn_icesal )
305      !
306      !               !---------------------------------------!
307      CASE( 1 )       !  constant salinity in time and space  !
308         !            !---------------------------------------!
309         sz_i(:,:,:,:) = rn_icesal
310         s_i (:,:,:)   = rn_icesal
311         !
312         !            !---------------------------------------------!
313      CASE( 2 )       !  time varying salinity with linear profile  !
314         !            !---------------------------------------------!
315         !
316         ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
317         !
318         DO jl = 1, jpl
319            DO jk = 1, nlay_i
320               sz_i(:,:,jk,jl)  = s_i(:,:,jl)
321            END DO
322         END DO
323         !                                      ! Slope of the linear profile
324         WHERE( h_i(:,:,:) > epsi20 )   ;   z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:)
325         ELSEWHERE                      ;   z_slope_s(:,:,:) = 0._wp
326         END WHERE
327         !
328         z1_dS = 1._wp / ( zsi1 - zsi0 )
329         DO jl = 1, jpl
330            DO jj = 1, jpj
331               DO ji = 1, jpi
332                  zalpha(ji,jj,jl) = MAX(  0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp )  )
333                  !                             ! force a constant profile when SSS too low (Baltic Sea)
334                  IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) )   zalpha(ji,jj,jl) = 0._wp 
335               END DO
336            END DO
337         END DO
338         !
339         ! Computation of the profile
340         DO jl = 1, jpl
341            DO jk = 1, nlay_i
342               DO jj = 1, jpj
343                  DO ji = 1, jpi
344                     !                          ! linear profile with 0 surface value
345                     zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
346                     zs  = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl)     ! weighting the profile
347                     sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
348                  END DO
349               END DO
350            END DO
351         END DO
352         !
353         DEALLOCATE( z_slope_s , zalpha )
354         !
355         !            !-------------------------------------------!
356      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
357         !            !-------------------------------------------!                                   (mean = 2.30)
358         !
359         s_i(:,:,:) = 2.30_wp
360!!gm Remark: if we keep the case 3, then compute an store one for all time-step
361!!           a array  S_prof(1:nlay_i) containing the calculation and just do:
362!         DO jk = 1, nlay_i
363!            sz_i(:,:,jk,:) = S_prof(jk)
364!         END DO
365!!gm end
366         !
367         DO jl = 1, jpl
368            DO jk = 1, nlay_i
369               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
370               sz_i(:,:,jk,jl) =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
371            END DO
372         END DO
373         !
374      END SELECT
375      !
376   END SUBROUTINE ice_var_salprof
377
378
379   SUBROUTINE ice_var_salprof1d
380      !!-------------------------------------------------------------------
381      !!                  ***  ROUTINE ice_var_salprof1d  ***
382      !!
383      !! ** Purpose :   1d computation of the sea ice salinity profile
384      !!                Works with 1d vectors and is used by thermodynamic modules
385      !!-------------------------------------------------------------------
386      INTEGER  ::   ji, jk    ! dummy loop indices
387      REAL(wp) ::   zargtemp, zsal, z1_dS   ! local scalars
388      REAL(wp) ::   zs, zs0              !   -      -
389      !
390      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z_slope_s, zalpha   !
391      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
392      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
393      !!-------------------------------------------------------------------
394      !
395      SELECT CASE ( nn_icesal )
396      !
397      !               !---------------------------------------!
398      CASE( 1 )       !  constant salinity in time and space  !
399         !            !---------------------------------------!
400         sz_i_1d(1:npti,:) = rn_icesal
401         !
402         !            !---------------------------------------------!
403      CASE( 2 )       !  time varying salinity with linear profile  !
404         !            !---------------------------------------------!
405         !
406         ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
407         !
408         !                                      ! Slope of the linear profile
409         WHERE( h_i_1d(1:npti) > epsi20 )   ;   z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
410         ELSEWHERE                          ;   z_slope_s(1:npti) = 0._wp
411         END WHERE
412         
413         z1_dS = 1._wp / ( zsi1 - zsi0 )
414         DO ji = 1, npti
415            zalpha(ji) = MAX(  0._wp , MIN(  ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp  )  )
416            !                             ! force a constant profile when SSS too low (Baltic Sea)
417            IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) )   zalpha(ji) = 0._wp
418         END DO
419         !
420         ! Computation of the profile
421         DO jk = 1, nlay_i
422            DO ji = 1, npti
423               !                          ! linear profile with 0 surface value
424               zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
425               zs  = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
426               sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
427            END DO
428         END DO
429         !
430         DEALLOCATE( z_slope_s, zalpha )
431
432         !            !-------------------------------------------!
433      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
434         !            !-------------------------------------------!                                   (mean = 2.30)
435         !
436         s_i_1d(1:npti) = 2.30_wp
437         !
438!!gm cf remark in ice_var_salprof routine, CASE( 3 )
439         DO jk = 1, nlay_i
440            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
441            zsal =  1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
442            DO ji = 1, npti
443               sz_i_1d(ji,jk) = zsal
444            END DO
445         END DO
446         !
447      END SELECT
448      !
449   END SUBROUTINE ice_var_salprof1d
450
451
452   SUBROUTINE ice_var_zapsmall
453      !!-------------------------------------------------------------------
454      !!                   ***  ROUTINE ice_var_zapsmall ***
455      !!
456      !! ** Purpose :   Remove too small sea ice areas and correct fluxes
457      !!-------------------------------------------------------------------
458      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
459      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch
460      !!-------------------------------------------------------------------
461      !
462      DO jl = 1, jpl       !==  loop over the categories  ==!
463         !
464         !-----------------------------------------------------------------
465         ! Zap ice energy and use ocean heat to melt ice
466         !-----------------------------------------------------------------
467         WHERE( a_i(:,:,jl) > epsi10 )   ;   h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
468         ELSEWHERE                       ;   h_i(:,:,jl) = 0._wp
469         END WHERE
470         !
471         WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 )   ;   zswitch(:,:) = 0._wp
472         ELSEWHERE                                                                           ;   zswitch(:,:) = 1._wp
473         END WHERE
474         !
475         DO jk = 1, nlay_i
476            DO jj = 1 , jpj
477               DO ji = 1 , jpi
478                  ! update exchanges with ocean
479                  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
480                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
481                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
482               END DO
483            END DO
484         END DO
485         !
486         DO jk = 1, nlay_s
487            DO jj = 1 , jpj
488               DO ji = 1 , jpi
489                  ! update exchanges with ocean
490                  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
491                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj)
492                  t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
493               END DO
494            END DO
495         END DO
496         !
497         DO jj = 1 , jpj
498            DO ji = 1 , jpi
499               ! update exchanges with ocean
500               sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoic * r1_rdtice
501               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoic * r1_rdtice
502               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhosn * r1_rdtice
503               !
504               !-----------------------------------------------------------------
505               ! zap ice and snow volume, add water and salt to ocean
506               !-----------------------------------------------------------------
507               a_i  (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
508               v_i  (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
509               v_s  (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
510               t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
511               oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
512               sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
513               !
514               h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
515               h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
516               !
517               a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
518               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
519               !
520            END DO
521         END DO
522         !
523      END DO 
524
525      ! to be sure that at_i is the sum of a_i(jl)
526      at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
527      vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
528
529      ! open water = 1 if at_i=0
530      WHERE( at_i(:,:) == 0._wp )   ato_i(:,:) = 1._wp
531      !
532   END SUBROUTINE ice_var_zapsmall
533
534
535   SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
536      !!-------------------------------------------------------------------
537      !!                ***  ROUTINE ice_var_itd   ***
538      !!
539      !! ** Purpose :  converting 1-cat ice to multiple ice categories
540      !!
541      !!                  ice thickness distribution follows a gaussian law
542      !!               around the concentration of the most likely ice thickness
543      !!                           (similar as iceistate.F90)
544      !!
545      !! ** Method:   Iterative procedure
546      !!               
547      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
548      !!
549      !!               2) Check whether the distribution conserves area and volume, positivity and
550      !!                  category boundaries
551      !!             
552      !!               3) If not (input ice is too thin), the last category is empty and
553      !!                  the number of categories is reduced (jpl-1)
554      !!
555      !!               4) Iterate until ok (SUM(itest(:) = 4)
556      !!
557      !! ** Arguments : zhti: 1-cat ice thickness
558      !!                zhts: 1-cat snow depth
559      !!                zati: 1-cat ice concentration
560      !!
561      !! ** Output    : jpl-cat
562      !!
563      !!  (Example of application: BDY forcings when input are cell averaged) 
564      !!-------------------------------------------------------------------
565      INTEGER  :: ji, jk, jl             ! dummy loop indices
566      INTEGER  :: idim, i_fill, jl0 
567      REAL(wp) :: zarg, zV, zconv, zdh, zdv
568      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
569      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i ! output ice/snow variables
570      INTEGER , DIMENSION(4)                  ::   itest
571      !!-------------------------------------------------------------------
572      !
573      ! ----------------------------------------
574      ! distribution over the jpl ice categories
575      ! ----------------------------------------
576      ! a gaussian distribution for ice concentration is used
577      ! then we check whether the distribution fullfills
578      ! volume and area conservation, positivity and ice categories bounds
579      idim = SIZE( zhti , 1 )
580      zh_i(1:idim,1:jpl) = 0._wp
581      zh_s(1:idim,1:jpl) = 0._wp
582      za_i(1:idim,1:jpl) = 0._wp
583      !
584      DO ji = 1, idim
585         !
586         IF( zhti(ji) > 0._wp ) THEN
587            !
588            ! find which category (jl0) the input ice thickness falls into
589            jl0 = jpl
590            DO jl = 1, jpl
591               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
592                  jl0 = jl
593                  CYCLE
594               ENDIF
595            END DO
596            !
597            itest(:) = 0
598            i_fill   = jpl + 1                                            !------------------------------------
599            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
600               !                                                          !------------------------------------
601               i_fill = i_fill - 1
602               !
603               zh_i(ji,1:jpl) = 0._wp
604               za_i(ji,1:jpl) = 0._wp
605               itest(:)       = 0     
606               !
607               IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
608                  zh_i(ji,1) = zhti(ji)
609                  za_i (ji,1) = zati (ji)
610               ELSE                         !-- case ice is thicker: fill categories >1
611                  ! thickness
612                  DO jl = 1, i_fill - 1
613                     zh_i(ji,jl) = hi_mean(jl)
614                  END DO
615                  !
616                  ! concentration
617                  za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
618                  DO jl = 1, i_fill - 1
619                     IF ( jl /= jl0 ) THEN
620                        zarg        = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
621                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
622                     ENDIF
623                  END DO
624                  !
625                  ! last category
626                  za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
627                  zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
628                  zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
629                  !
630                  ! correction if concentration of upper cat is greater than lower cat
631                  !    (it should be a gaussian around jl0 but sometimes it is not)
632                  IF ( jl0 /= jpl ) THEN
633                     DO jl = jpl, jl0+1, -1
634                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
635                           zdv = zh_i(ji,jl) * za_i(ji,jl)
636                           zh_i(ji,jl    ) = 0._wp
637                           za_i (ji,jl    ) = 0._wp
638                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
639                        END IF
640                     END DO
641                  ENDIF
642                  !
643               ENDIF
644               !
645               ! Compatibility tests
646               zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) ) 
647               IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation
648               !
649               zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )
650               IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation
651               !
652               IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ?
653               !
654               itest(4) = 1
655               DO jl = 1, i_fill
656                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0                                ! Test 4: positivity of ice concentrations
657               END DO
658               !                                         !----------------------------
659            END DO                                       ! end iteration on categories
660            !                                            !----------------------------
661         ENDIF
662      END DO
663
664      ! Add Snow in each category where za_i is not 0
665      DO jl = 1, jpl
666         DO ji = 1, idim
667            IF( za_i(ji,jl) > 0._wp ) THEN
668               zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
669               ! In case snow load is in excess that would lead to transformation from snow to ice
670               ! Then, transfer the snow excess into the ice (different from icethd_dh)
671               zdh = MAX( 0._wp, ( rhosn * zh_s(ji,jl) + ( rhoic - rau0 ) * zh_i(ji,jl) ) * r1_rau0 ) 
672               ! recompute h_i, h_s avoiding out of bounds values
673               zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
674               zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoic * r1_rhosn )
675            ENDIF
676         END DO
677      END DO
678      !
679   END SUBROUTINE ice_var_itd
680
681
682   SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )
683      !!-------------------------------------------------------------------
684      !!                ***  ROUTINE ice_var_itd2   ***
685      !!
686      !! ** Purpose :  converting N-cat ice to jpl ice categories
687      !!
688      !!                  ice thickness distribution follows a gaussian law
689      !!               around the concentration of the most likely ice thickness
690      !!                           (similar as iceistate.F90)
691      !!
692      !! ** Method:   Iterative procedure
693      !!               
694      !!               1) Fill ice cat that correspond to input thicknesses
695      !!                  Find the lowest(jlmin) and highest(jlmax) cat that are filled
696      !!
697      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1
698      !!                   by removing 25% ice area from jlmin and jlmax (resp.)
699      !!             
700      !!               3) Expand the filling to the empty cat between jlmin and jlmax
701      !!                   by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
702      !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
703      !!
704      !! ** Arguments : zhti: N-cat ice thickness
705      !!                zhts: N-cat snow depth
706      !!                zati: N-cat ice concentration
707      !!
708      !! ** Output    : jpl-cat
709      !!
710      !!  (Example of application: BDY forcings when inputs have N-cat /= jpl) 
711      !!-------------------------------------------------------------------
712      INTEGER  ::   ji, jl, jl1, jl2             ! dummy loop indices
713      INTEGER  ::   idim, icat 
714      INTEGER, PARAMETER ::   ztrans = 0.25_wp
715      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
716      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
717      INTEGER , DIMENSION(:,:), ALLOCATABLE   ::   jlfil, jlfil2
718      INTEGER , DIMENSION(:)  , ALLOCATABLE   ::   jlmax, jlmin
719      !!-------------------------------------------------------------------
720      !
721      idim = SIZE( zhti, 1 )
722      icat = SIZE( zhti, 2 )
723      !
724      ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays
725      ALLOCATE( jlmin(idim), jlmax(idim) )
726
727      ! --- initialize output fields to 0 --- !
728      zh_i(1:idim,1:jpl) = 0._wp
729      zh_s(1:idim,1:jpl) = 0._wp
730      za_i(1:idim,1:jpl) = 0._wp
731      !
732      ! --- fill the categories --- !
733      !     find where cat-input = cat-output and fill cat-output fields 
734      jlmax(:) = 0
735      jlmin(:) = 999
736      jlfil(:,:) = 0
737      DO jl1 = 1, jpl
738         DO jl2 = 1, icat
739            DO ji = 1, idim
740               IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
741                  ! fill the right category
742                  zh_i(ji,jl1) = zhti(ji,jl2)
743                  zh_s(ji,jl1) = zhts(ji,jl2)
744                  za_i(ji,jl1) = zati(ji,jl2)
745                  ! record categories that are filled
746                  jlmax(ji) = MAX( jlmax(ji), jl1 )
747                  jlmin(ji) = MIN( jlmin(ji), jl1 )
748                  jlfil(ji,jl1) = jl1
749               ENDIF
750            END DO
751         END DO
752      END DO
753      !
754      ! --- fill the gaps between categories --- ! 
755      !     transfer from categories filled at the previous step to the empty ones in between
756      DO ji = 1, idim
757         jl1 = jlmin(ji)
758         jl2 = jlmax(ji)
759         IF( jl1 > 1 ) THEN
760            ! fill the lower cat (jl1-1)
761            za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
762            zh_i(ji,jl1-1) = hi_mean(jl1-1)
763            ! remove from cat jl1
764            za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
765         ENDIF
766         IF( jl2 < jpl ) THEN
767            ! fill the upper cat (jl2+1)
768            za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
769            zh_i(ji,jl2+1) = hi_mean(jl2+1)
770            ! remove from cat jl2
771            za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
772         ENDIF
773      END DO
774      !
775      jlfil2(:,:) = jlfil(:,:) 
776      ! fill categories from low to high
777      DO jl = 2, jpl-1
778         DO ji = 1, idim
779            IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
780               ! fill high
781               za_i(ji,jl) = ztrans * za_i(ji,jl-1)
782               zh_i(ji,jl) = hi_mean(jl)
783               jlfil(ji,jl) = jl
784               ! remove low
785               za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
786            ENDIF
787         END DO
788      END DO
789      !
790      ! fill categories from high to low
791      DO jl = jpl-1, 2, -1
792         DO ji = 1, idim
793            IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
794               ! fill low
795               za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
796               zh_i(ji,jl) = hi_mean(jl) 
797               jlfil2(ji,jl) = jl
798               ! remove high
799               za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
800            ENDIF
801         END DO
802      END DO
803      !
804      DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
805      DEALLOCATE( jlmin, jlmax )
806      !
807   END SUBROUTINE ice_var_itd2
808
809
810   SUBROUTINE ice_var_bv
811      !!-------------------------------------------------------------------
812      !!                ***  ROUTINE ice_var_bv ***
813      !!
814      !! ** Purpose :   computes mean brine volume (%) in sea ice
815      !!
816      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
817      !!
818      !! References : Vancoppenolle et al., JGR, 2007
819      !!-------------------------------------------------------------------
820      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
821      !!-------------------------------------------------------------------
822      !
823!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
824!!   instead of setting everything to zero as just below
825      bv_i (:,:,:) = 0._wp
826      DO jl = 1, jpl
827         DO jk = 1, nlay_i
828            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
829               bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
830            END WHERE
831         END DO
832      END DO
833      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
834      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
835      END WHERE
836      !
837   END SUBROUTINE ice_var_bv
838
839
840   SUBROUTINE ice_var_enthalpy
841      !!-------------------------------------------------------------------
842      !!                   ***  ROUTINE ice_var_enthalpy ***
843      !!                 
844      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
845      !!
846      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
847      !!-------------------------------------------------------------------
848      INTEGER  ::   ji, jk   ! dummy loop indices
849      REAL(wp) ::   ztmelts  ! local scalar
850      !!-------------------------------------------------------------------
851      !
852      DO jk = 1, nlay_i             ! Sea ice energy of melting
853         DO ji = 1, npti
854            ztmelts      = - tmut  * sz_i_1d(ji,jk)
855            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point
856                                                                !   (sometimes zdf scheme produces abnormally high temperatures)   
857            e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
858               &                    + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
859               &                    - rcp  *   ztmelts )
860         END DO
861      END DO
862      DO jk = 1, nlay_s             ! Snow energy of melting
863         DO ji = 1, npti
864            e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
865         END DO
866      END DO
867      !
868   END SUBROUTINE ice_var_enthalpy
869
870#else
871   !!----------------------------------------------------------------------
872   !!   Default option         Dummy module           NO SI3 sea-ice model
873   !!----------------------------------------------------------------------
874#endif
875
876   !!======================================================================
877END MODULE icevar
Note: See TracBrowser for help on using the repository browser.