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 branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90 @ 8623

Last change on this file since 8623 was 8623, checked in by clem, 6 years ago

rewrite melt ponds

File size: 32.2 KB
Line 
1MODULE icevar
2   !!======================================================================
3   !!                       ***  MODULE icevar ***
4   !!   sea-ice:     Different sets of ice model variables
5   !!                   how to switch from one to another
6   !!
7   !!                 There are three sets of variables
8   !!                 VGLO : global variables of the model
9   !!                        - v_i (jpi,jpj,jpl)
10   !!                        - v_s (jpi,jpj,jpl)
11   !!                        - a_i (jpi,jpj,jpl)
12   !!                        - t_s (jpi,jpj,jpl)
13   !!                        - e_i (jpi,jpj,nlay_i,jpl)
14   !!                        - sv_i(jpi,jpj,jpl)
15   !!                        - oa_i(jpi,jpj,jpl)
16   !!                 VEQV : equivalent variables sometimes used in the model
17   !!                        - h_i(jpi,jpj,jpl)
18   !!                        - h_s(jpi,jpj,jpl)
19   !!                        - t_i(jpi,jpj,nlay_i,jpl)
20   !!                        ...
21   !!                 VAGG : aggregate variables, averaged/summed over all
22   !!                        thickness categories
23   !!                        - vt_i(jpi,jpj)
24   !!                        - vt_s(jpi,jpj)
25   !!                        - at_i(jpi,jpj)
26   !!                        - et_s(jpi,jpj)  !total snow heat content
27   !!                        - et_i(jpi,jpj)  !total ice thermal content
28   !!                        - sm_i(jpi,jpj) !mean ice salinity
29   !!                        - tm_i (jpi,jpj) !mean ice temperature
30   !!======================================================================
31   !! History :   -   ! 2006-01 (M. Vancoppenolle) Original code
32   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
33   !!            3.5  ! 2012    (M. Vancoppenolle)  add ice_var_itd
34   !!            3.6  ! 2014-01 (C. Rousset) add ice_var_zapsmall, rewrite ice_var_itd
35   !!----------------------------------------------------------------------
36#if defined key_lim3
37   !!----------------------------------------------------------------------
38   !!   'key_lim3'                                       ESIM sea-ice model
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 multiple cat
47   !!   ice_var_bv        : brine volume
48   !!----------------------------------------------------------------------
49   USE dom_oce        ! ocean space and time domain
50   USE phycst         ! physical constants (ocean directory)
51   USE sbc_oce , ONLY : sss_m
52   USE ice            ! sea-ice: variables
53   USE ice1D          ! sea-ice: thermodynamics variables
54   !
55   USE in_out_manager ! I/O manager
56   USE lib_mpp        ! MPP library
57   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
58
59   IMPLICIT NONE
60   PRIVATE
61
62   PUBLIC   ice_var_agg         
63   PUBLIC   ice_var_glo2eqv     
64   PUBLIC   ice_var_eqv2glo     
65   PUBLIC   ice_var_salprof     
66   PUBLIC   ice_var_salprof1d   
67   PUBLIC   ice_var_zapsmall
68   PUBLIC   ice_var_itd
69   PUBLIC   ice_var_bv           
70
71   !!----------------------------------------------------------------------
72   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
73   !! $Id: icevar.F90 8422 2017-08-08 13:58:05Z clem $
74   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
75   !!----------------------------------------------------------------------
76CONTAINS
77
78   SUBROUTINE ice_var_agg( kn )
79      !!-------------------------------------------------------------------
80      !!                ***  ROUTINE ice_var_agg  ***
81      !!
82      !! ** Purpose :   aggregates ice-thickness-category variables to
83      !!              all-ice variables, i.e. it turns VGLO into VAGG
84      !!-------------------------------------------------------------------
85      INTEGER, INTENT( in ) ::   kn     ! =1 state variables only
86      !                                 ! >1 state variables + others
87      !
88      INTEGER ::   ji, jj, jk, jl   ! dummy loop indices
89      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z1_at_i, z1_vt_i
90      !!-------------------------------------------------------------------
91      !
92      !                                      ! integrated values
93      vt_i(:,:) =       SUM( v_i(:,:,:)           , dim=3 )
94      vt_s(:,:) =       SUM( v_s(:,:,:)           , dim=3 )
95      at_i(:,:) =       SUM( a_i(:,:,:)           , dim=3 )
96      et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 )
97      et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
98
99      at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds
100      vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
101
102      ato_i(:,:) = 1._wp - at_i(:,:)    ! open water fraction 
103
104      IF( kn > 1 ) THEN
105         !
106         ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) )
107         WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:)
108         ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp
109         END WHERE
110         WHERE( vt_i(:,:) > epsi20 )   ;   z1_vt_i(:,:) = 1._wp / vt_i(:,:)
111         ELSEWHERE                     ;   z1_vt_i(:,:) = 0._wp
112         END WHERE
113         !
114         !                          ! mean ice/snow thickness
115         hm_i(:,:) = vt_i(:,:) * z1_at_i(:,:)
116         hm_s(:,:) = vt_s(:,:) * z1_at_i(:,:)
117         !         
118         !                          ! mean temperature (K), salinity and age
119         tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
120         tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
121         om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(:,:)
122         !
123         tm_i(:,:) = 0._wp
124         sm_i(:,:) = 0._wp
125         DO jl = 1, jpl
126            DO jk = 1, nlay_i
127               tm_i(:,:) = tm_i(:,:) + r1_nlay_i * t_i (:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
128               sm_i(:,:) = sm_i(:,:) + r1_nlay_i * sz_i(:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
129            END DO
130         END DO
131         !
132         !                           ! put rt0 where there is no ice
133         WHERE( at_i(:,:)<=epsi20 )
134            tm_su(:,:) = rt0
135            tm_si(:,:) = rt0
136            tm_i (:,:) = rt0
137         END WHERE
138
139         DEALLOCATE( z1_at_i , z1_vt_i )
140      ENDIF
141      !
142   END SUBROUTINE ice_var_agg
143
144
145   SUBROUTINE ice_var_glo2eqv
146      !!-------------------------------------------------------------------
147      !!                ***  ROUTINE ice_var_glo2eqv ***
148      !!
149      !! ** Purpose :   computes equivalent variables as function of 
150      !!              global variables, i.e. it turns VGLO into VEQV
151      !!-------------------------------------------------------------------
152      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
153      REAL(wp) ::   ze_i             ! local scalars
154      REAL(wp) ::   ze_s, ztmelts, zbbb, zccc       !   -      -
155      REAL(wp) ::   zhmax, z1_zhmax                 !   -      -
156      REAL(wp) ::   zlay_i, zlay_s                  !   -      -
157      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i
158      !!-------------------------------------------------------------------
159
160!!gm Question 2:  It is possible to define existence of sea-ice in a common way between
161!!                ice area and ice volume ?
162!!                the idea is to be able to define one for all at the begining of this routine
163!!                a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
164
165      !---------------------------------------------------------------
166      ! Ice thickness, snow thickness, ice salinity, ice age and ponds
167      !---------------------------------------------------------------
168      !                                            !--- inverse of the ice area
169      WHERE( a_i(:,:,:) > epsi20 )   ;   z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
170      ELSEWHERE                      ;   z1_a_i(:,:,:) = 0._wp
171      END WHERE
172      !
173      WHERE( v_i(:,:,:) > epsi20 )   ;   z1_v_i(:,:,:) = 1._wp / v_i(:,:,:)
174      ELSEWHERE                      ;   z1_v_i(:,:,:) = 0._wp
175      END WHERE
176      !                                           !--- ice thickness
177      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
178
179      zhmax    =          hi_max(jpl)
180      z1_zhmax =  1._wp / hi_max(jpl)               
181      WHERE( h_i(:,:,jpl) > zhmax )   ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area
182         h_i  (:,:,jpl) = zhmax
183         a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 
184         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)
185      END WHERE
186      !                                           !--- snow thickness
187      h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:)
188      !                                           !--- ice age     
189      o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:)
190      !                                           !--- pond fraction and thickness     
191      a_ip_frac(:,:,:) = a_ip(:,:,:) * z1_a_i(:,:,:)
192      WHERE( a_ip_frac(:,:,:) > epsi20 )   ;   h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:)
193      ELSEWHERE                            ;   h_ip(:,:,:) = 0._wp
194      END WHERE
195      !
196      !                                           !---  salinity (with a minimum value imposed everywhere)     
197      IF( nn_icesal == 2 ) THEN
198         WHERE( v_i(:,:,:) > epsi20 )   ;   s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) * z1_v_i(:,:,:) ) )
199         ELSEWHERE                      ;   s_i(:,:,:) = rn_simin
200         END WHERE
201      ENDIF
202      CALL ice_var_salprof   ! salinity profile
203
204      !-------------------
205      ! Ice temperature   [K]   (with a minimum value (rt0 - 100.))
206      !-------------------
207      zlay_i   = REAL( nlay_i , wp )    ! number of layers
208      DO jl = 1, jpl
209         DO jk = 1, nlay_i
210            DO jj = 1, jpj
211               DO ji = 1, jpi
212                  IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area
213                     !
214                     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]
215                     ztmelts          = - sz_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C]
216                     ! Conversion q(S,T) -> T (second order equation)
217                     zbbb             = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus
218                     zccc             = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) )
219                     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
220                     !
221                  ELSE                                !--- no ice
222                     t_i(ji,jj,jk,jl) = rt0
223                  ENDIF
224               END DO
225            END DO
226         END DO
227      END DO
228
229      !--------------------
230      ! Snow temperature   [K]   (with a minimum value (rt0 - 100.))
231      !--------------------
232      zlay_s = REAL( nlay_s , wp )
233      DO jk = 1, nlay_s
234         WHERE( v_s(:,:,:) > epsi20 )        !--- icy area
235            t_s(:,:,jk,:) = MAX( -100._wp , MIN( r1_cpic * ( -r1_rhosn * (e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s) + lfus ) , 0._wp ) ) + rt0
236         ELSEWHERE                           !--- no ice
237            t_s(:,:,jk,:) = rt0
238         END WHERE
239      END DO
240
241      ! integrated values
242      vt_i (:,:) = SUM( v_i, dim=3 )
243      vt_s (:,:) = SUM( v_s, dim=3 )
244      at_i (:,:) = SUM( a_i, dim=3 )
245      !
246   END SUBROUTINE ice_var_glo2eqv
247
248
249   SUBROUTINE ice_var_eqv2glo
250      !!-------------------------------------------------------------------
251      !!                ***  ROUTINE ice_var_eqv2glo ***
252      !!
253      !! ** Purpose :   computes global variables as function of
254      !!              equivalent variables,  i.e. it turns VEQV into VGLO
255      !!-------------------------------------------------------------------
256      !
257      v_i (:,:,:) = h_i (:,:,:) * a_i (:,:,:)
258      v_s (:,:,:) = h_s (:,:,:) * a_i (:,:,:)
259      sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:)
260      v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
261      !
262   END SUBROUTINE ice_var_eqv2glo
263
264
265   SUBROUTINE ice_var_salprof
266      !!-------------------------------------------------------------------
267      !!                ***  ROUTINE ice_var_salprof ***
268      !!
269      !! ** Purpose :   computes salinity profile in function of bulk salinity     
270      !!
271      !! ** Method  : If bulk salinity greater than zsi1,
272      !!              the profile is assumed to be constant (S_inf)
273      !!              If bulk salinity lower than zsi0,
274      !!              the profile is linear with 0 at the surface (S_zero)
275      !!              If it is between zsi0 and zsi1, it is a
276      !!              alpha-weighted linear combination of s_inf and s_zero
277      !!
278      !! ** References : Vancoppenolle et al., 2007
279      !!-------------------------------------------------------------------
280      INTEGER  ::   ji, jj, jk, jl   ! dummy loop index
281      REAL(wp) ::   zsal, z1_dS
282      REAL(wp) ::   zargtemp , zs0, zs
283      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z_slope_s, zalpha    ! case 2 only
284      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
285      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
286      !!-------------------------------------------------------------------
287
288!!gm Question: Remove the option 3 ?  How many years since it last use ?
289
290      SELECT CASE ( nn_icesal )
291      !
292      !               !---------------------------------------!
293      CASE( 1 )       !  constant salinity in time and space  !
294         !            !---------------------------------------!
295         sz_i(:,:,:,:) = rn_icesal
296         s_i(:,:,:)   = rn_icesal
297         !
298         !            !---------------------------------------------!
299      CASE( 2 )       !  time varying salinity with linear profile  !
300         !            !---------------------------------------------!
301         !
302         ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
303         !
304         DO jl = 1, jpl
305            DO jk = 1, nlay_i
306               sz_i(:,:,jk,jl)  = s_i(:,:,jl)
307            END DO
308         END DO
309         !                                      ! Slope of the linear profile
310         WHERE( h_i(:,:,:) > epsi20 )   ;   z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:)
311         ELSEWHERE                       ;   z_slope_s(:,:,:) = 0._wp
312         END WHERE
313         !
314         z1_dS = 1._wp / ( zsi1 - zsi0 )
315         DO jl = 1, jpl
316            DO jj = 1, jpj
317               DO ji = 1, jpi
318                  zalpha(ji,jj,jl) = MAX(  0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp )  )
319                  !                             ! force a constant profile when SSS too low (Baltic Sea)
320                  IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) )   zalpha(ji,jj,jl) = 0._wp 
321               END DO
322            END DO
323         END DO
324
325         ! Computation of the profile
326         DO jl = 1, jpl
327            DO jk = 1, nlay_i
328               DO jj = 1, jpj
329                  DO ji = 1, jpi
330                     !                          ! linear profile with 0 surface value
331                     zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
332                     zs  = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl)     ! weighting the profile
333                     sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
334                  END DO
335               END DO
336            END DO
337         END DO
338         !
339         DEALLOCATE( z_slope_s , zalpha )
340         !
341         !            !-------------------------------------------!
342      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
343         !            !-------------------------------------------!                                   (mean = 2.30)
344         !
345         s_i(:,:,:) = 2.30_wp
346!!gm Remark: if we keep the case 3, then compute an store one for all time-step
347!!           a array  S_prof(1:nlay_i) containing the calculation and just do:
348!         DO jk = 1, nlay_i
349!            sz_i(:,:,jk,:) = S_prof(jk)
350!         END DO
351!!gm end
352         !
353         DO jl = 1, jpl
354            DO jk = 1, nlay_i
355               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
356               sz_i(:,:,jk,jl) =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
357            END DO
358         END DO
359         !
360      END SELECT
361      !
362   END SUBROUTINE ice_var_salprof
363
364   SUBROUTINE ice_var_salprof1d
365      !!-------------------------------------------------------------------
366      !!                  ***  ROUTINE ice_var_salprof1d  ***
367      !!
368      !! ** Purpose :   1d computation of the sea ice salinity profile
369      !!                Works with 1d vectors and is used by thermodynamic modules
370      !!-------------------------------------------------------------------
371      INTEGER  ::   ji, jk    ! dummy loop indices
372      REAL(wp) ::   zargtemp, zsal, z1_dS   ! local scalars
373      REAL(wp) ::   zs, zs0              !   -      -
374      !
375      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z_slope_s, zalpha   !
376      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
377      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
378      !!-------------------------------------------------------------------
379      !
380      SELECT CASE ( nn_icesal )
381      !
382      !               !---------------------------------------!
383      CASE( 1 )       !  constant salinity in time and space  !
384         !            !---------------------------------------!
385         sz_i_1d(1:npti,:) = rn_icesal
386         !
387         !            !---------------------------------------------!
388      CASE( 2 )       !  time varying salinity with linear profile  !
389         !            !---------------------------------------------!
390         !
391         ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
392         !
393         !                                      ! Slope of the linear profile
394         WHERE( h_i_1d(1:npti) > epsi20 )   ;   z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
395         ELSEWHERE                           ;   z_slope_s(1:npti) = 0._wp
396         END WHERE
397         
398         z1_dS = 1._wp / ( zsi1 - zsi0 )
399         DO ji = 1, npti
400            zalpha(ji) = MAX(  0._wp , MIN(  ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp  )  )
401            !                             ! force a constant profile when SSS too low (Baltic Sea)
402            IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) )   zalpha(ji) = 0._wp
403         END DO
404         !
405         ! Computation of the profile
406         DO jk = 1, nlay_i
407            DO ji = 1, npti
408               !                          ! linear profile with 0 surface value
409               zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
410               zs  = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
411               sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
412            END DO
413         END DO
414         !
415         DEALLOCATE( z_slope_s, zalpha )
416
417         !            !-------------------------------------------!
418      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
419         !            !-------------------------------------------!                                   (mean = 2.30)
420         !
421         s_i_1d(1:npti) = 2.30_wp
422         !
423!!gm cf remark in ice_var_salprof routine, CASE( 3 )
424         DO jk = 1, nlay_i
425            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
426            zsal =  1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
427            DO ji = 1, npti
428               sz_i_1d(ji,jk) = zsal
429            END DO
430         END DO
431         !
432      END SELECT
433      !
434   END SUBROUTINE ice_var_salprof1d
435
436
437   SUBROUTINE ice_var_zapsmall
438      !!-------------------------------------------------------------------
439      !!                   ***  ROUTINE ice_var_zapsmall ***
440      !!
441      !! ** Purpose :   Remove too small sea ice areas and correct fluxes
442      !!-------------------------------------------------------------------
443      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
444      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch
445      !!-------------------------------------------------------------------
446      !
447      DO jl = 1, jpl       !==  loop over the categories  ==!
448         !
449         !-----------------------------------------------------------------
450         ! Zap ice energy and use ocean heat to melt ice
451         !-----------------------------------------------------------------
452         WHERE( a_i(:,:,jl) > epsi20 )   ;   h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
453         ELSEWHERE                       ;   h_i(:,:,jl) = 0._wp
454         END WHERE
455         !
456         WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 )   ;   zswitch(:,:) = 0._wp
457         ELSEWHERE                                                                           ;   zswitch(:,:) = 1._wp
458         END WHERE
459         
460         DO jk = 1, nlay_i
461            DO jj = 1 , jpj
462               DO ji = 1 , jpi
463                  ! update exchanges with ocean
464                  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
465                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
466                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
467               END DO
468            END DO
469         END DO
470
471         DO jj = 1 , jpj
472            DO ji = 1 , jpi
473               ! update exchanges with ocean
474               sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoic * r1_rdtice
475               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoic * r1_rdtice
476               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhosn * r1_rdtice
477               hfx_res(ji,jj)  = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s (ji,jj,1,jl) * r1_rdtice ! W.m-2 <0
478               IF( ln_pnd_fwb ) THEN
479                  wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice
480               ENDIF
481               !-----------------------------------------------------------------
482               ! Zap snow energy
483               !-----------------------------------------------------------------
484               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
485               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zswitch(ji,jj)
486
487               !-----------------------------------------------------------------
488               ! zap ice and snow volume, add water and salt to ocean
489               !-----------------------------------------------------------------
490               ato_i(ji,jj)    = a_i (ji,jj,jl) * ( 1._wp - zswitch(ji,jj) ) + ato_i(ji,jj)
491               a_i  (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
492               v_i  (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
493               v_s  (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
494               t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
495               oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
496               sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
497
498               h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
499               h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
500
501               a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
502               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
503
504            END DO
505         END DO
506       
507      END DO 
508
509      ! to be sure that at_i is the sum of a_i(jl)
510      at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
511      vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
512
513      ! open water = 1 if at_i=0
514      WHERE( at_i(:,:) == 0._wp )   ato_i(:,:) = 1._wp
515
516      !
517   END SUBROUTINE ice_var_zapsmall
518
519
520   SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
521      !!-------------------------------------------------------------------
522      !!                ***  ROUTINE ice_var_itd   ***
523      !!
524      !! ** Purpose :  converting 1-cat ice to multiple ice categories
525      !!
526      !!                  ice thickness distribution follows a gaussian law
527      !!               around the concentration of the most likely ice thickness
528      !!                           (similar as iceistate.F90)
529      !!
530      !! ** Method:   Iterative procedure
531      !!               
532      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
533      !!
534      !!               2) Check whether the distribution conserves area and volume, positivity and
535      !!                  category boundaries
536      !!             
537      !!               3) If not (input ice is too thin), the last category is empty and
538      !!                  the number of categories is reduced (jpl-1)
539      !!
540      !!               4) Iterate until ok (SUM(itest(:) = 4)
541      !!
542      !! ** Arguments : zhti: 1-cat ice thickness
543      !!                zhts: 1-cat snow depth
544      !!                zati: 1-cat ice concentration
545      !!
546      !! ** Output    : jpl-cat
547      !!
548      !!  (Example of application: BDY forcings when input are cell averaged) 
549      !!-------------------------------------------------------------------
550      INTEGER  :: ji, jk, jl             ! dummy loop indices
551      INTEGER  :: ijpij, i_fill, jl0 
552      REAL(wp) :: zarg, zV, zconv, zdh, zdv
553      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
554      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i ! output ice/snow variables
555      INTEGER , DIMENSION(4)                  ::   itest
556      !!-------------------------------------------------------------------
557      !
558      ! ----------------------------------------
559      ! distribution over the jpl ice categories
560      ! ----------------------------------------
561      ! a gaussian distribution for ice concentration is used
562      ! then we check whether the distribution fullfills
563      ! volume and area conservation, positivity and ice categories bounds
564      ijpij = SIZE( zhti , 1 )
565      zh_i(1:ijpij,1:jpl) = 0._wp
566      zh_s(1:ijpij,1:jpl) = 0._wp
567      za_i (1:ijpij,1:jpl) = 0._wp
568
569      DO ji = 1, ijpij
570         
571         IF( zhti(ji) > 0._wp ) THEN
572
573            ! find which category (jl0) the input ice thickness falls into
574            jl0 = jpl
575            DO jl = 1, jpl
576               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
577                  jl0 = jl
578                  CYCLE
579               ENDIF
580            END DO
581
582            itest(:) = 0
583            i_fill   = jpl + 1                                            !------------------------------------
584            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
585               !                                                          !------------------------------------
586               i_fill = i_fill - 1
587               !
588               zh_i(ji,1:jpl) = 0._wp
589               za_i (ji,1:jpl) = 0._wp
590               itest(:)        = 0     
591               
592               IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
593                  zh_i(ji,1) = zhti(ji)
594                  za_i (ji,1) = zati (ji)
595               ELSE                         !-- case ice is thicker: fill categories >1
596                  ! thickness
597                  DO jl = 1, i_fill - 1
598                     zh_i(ji,jl) = hi_mean(jl)
599                  END DO
600                 
601                  ! concentration
602                  za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
603                  DO jl = 1, i_fill - 1
604                     IF ( jl /= jl0 ) THEN
605                        zarg        = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
606                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
607                     ENDIF
608                  END DO
609                 
610                  ! last category
611                  za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
612                  zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
613                  zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
614                 
615                  ! clem: correction if concentration of upper cat is greater than lower cat
616                  !       (it should be a gaussian around jl0 but sometimes it is not)
617                  IF ( jl0 /= jpl ) THEN
618                     DO jl = jpl, jl0+1, -1
619                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
620                           zdv = zh_i(ji,jl) * za_i(ji,jl)
621                           zh_i(ji,jl    ) = 0._wp
622                           za_i (ji,jl    ) = 0._wp
623                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
624                        END IF
625                     ENDDO
626                  ENDIF
627               
628               ENDIF
629           
630               ! Compatibility tests
631               zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) ) 
632               IF ( zconv < epsi06 ) itest(1) = 1                                        ! Test 1: area conservation
633           
634               zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )
635               IF ( zconv < epsi06 ) itest(2) = 1                                        ! Test 2: volume conservation
636               
637               IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ?
638               
639               itest(4) = 1
640               DO jl = 1, i_fill
641                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0                                ! Test 4: positivity of ice concentrations
642               END DO
643               !                                         !----------------------------
644            END DO                                       ! end iteration on categories
645               !                                         !----------------------------
646         ENDIF
647      END DO
648
649      ! Add Snow in each category where za_i is not 0
650      DO jl = 1, jpl
651         DO ji = 1, ijpij
652            IF( za_i(ji,jl) > 0._wp ) THEN
653               zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
654               ! In case snow load is in excess that would lead to transformation from snow to ice
655               ! Then, transfer the snow excess into the ice (different from icethd_dh)
656               zdh = MAX( 0._wp, ( rhosn * zh_s(ji,jl) + ( rhoic - rau0 ) * zh_i(ji,jl) ) * r1_rau0 ) 
657               ! recompute h_i, h_s avoiding out of bounds values
658               zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
659               zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoic * r1_rhosn )
660            ENDIF
661         END DO
662      END DO
663      !
664    END SUBROUTINE ice_var_itd
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) - tmut * 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#else
698   !!----------------------------------------------------------------------
699   !!   Default option         Dummy module           NO ESIM sea-ice model
700   !!----------------------------------------------------------------------
701#endif
702
703   !!======================================================================
704END MODULE icevar
Note: See TracBrowser for help on using the repository browser.