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 @ 8592

Last change on this file since 8592 was 8592, checked in by clem, 7 years ago

rewrite ice albedo routine
File size: 32.4 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      ! MV MP 2016
100      IF ( ln_pnd ) THEN                     ! Melt pond
101         at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 )
102         vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
103      ENDIF
104      ! END MP 2016
105
106      ato_i(:,:) = 1._wp - at_i(:,:)    ! open water fraction 
107
108      IF( kn > 1 ) THEN
109         !
110         ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) )
111         WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:)
112         ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp
113         END WHERE
114         WHERE( vt_i(:,:) > epsi20 )   ;   z1_vt_i(:,:) = 1._wp / vt_i(:,:)
115         ELSEWHERE                     ;   z1_vt_i(:,:) = 0._wp
116         END WHERE
117         !
118         !                          ! mean ice/snow thickness
119         hm_i(:,:) = vt_i(:,:) * z1_at_i(:,:)
120         hm_s(:,:) = vt_s(:,:) * z1_at_i(:,:)
121         !         
122         !                          ! mean temperature (K), salinity and age
123         tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
124         tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
125         om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(:,:)
126         !
127         tm_i(:,:) = 0._wp
128         sm_i(:,:) = 0._wp
129         DO jl = 1, jpl
130            DO jk = 1, nlay_i
131               tm_i(:,:) = tm_i(:,:) + r1_nlay_i * t_i (:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
132               sm_i(:,:) = sm_i(:,:) + r1_nlay_i * sz_i(:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
133            END DO
134         END DO
135         !
136         !                           ! put rt0 where there is no ice
137         WHERE( at_i(:,:)<=epsi20 )
138            tm_su(:,:) = rt0
139            tm_si(:,:) = rt0
140            tm_i (:,:) = rt0
141         END WHERE
142
143         DEALLOCATE( z1_at_i , z1_vt_i )
144      ENDIF
145      !
146   END SUBROUTINE ice_var_agg
147
148
149   SUBROUTINE ice_var_glo2eqv
150      !!-------------------------------------------------------------------
151      !!                ***  ROUTINE ice_var_glo2eqv ***
152      !!
153      !! ** Purpose :   computes equivalent variables as function of 
154      !!              global variables, i.e. it turns VGLO into VEQV
155      !!-------------------------------------------------------------------
156      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
157      REAL(wp) ::   ze_i             ! local scalars
158      REAL(wp) ::   ze_s, ztmelts, zbbb, zccc       !   -      -
159      REAL(wp) ::   zhmax, z1_zhmax                 !   -      -
160      REAL(wp) ::   zlay_i, zlay_s                  !   -      -
161      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i
162      !!-------------------------------------------------------------------
163
164!!gm Question 2:  It is possible to define existence of sea-ice in a common way between
165!!                ice area and ice volume ?
166!!                the idea is to be able to define one for all at the begining of this routine
167!!                a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
168
169      !---------------------------------------------------------------
170      ! Ice thickness, snow thickness, ice salinity, ice age and ponds
171      !---------------------------------------------------------------
172      !                                            !--- inverse of the ice area
173      WHERE( a_i(:,:,:) > epsi20 )   ;   z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
174      ELSEWHERE                      ;   z1_a_i(:,:,:) = 0._wp
175      END WHERE
176      !
177      WHERE( v_i(:,:,:) > epsi20 )   ;   z1_v_i(:,:,:) = 1._wp / v_i(:,:,:)
178      ELSEWHERE                      ;   z1_v_i(:,:,:) = 0._wp
179      END WHERE
180      !                                           !--- ice thickness
181      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
182
183      zhmax    =          hi_max(jpl)
184      z1_zhmax =  1._wp / hi_max(jpl)               
185      WHERE( h_i(:,:,jpl) > zhmax )   ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area
186         h_i  (:,:,jpl) = zhmax
187         a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 
188         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)
189      END WHERE
190      !                                           !--- snow thickness
191      h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:)
192      !                                           !--- ice age     
193      o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:)
194      !                                           !--- pond fraction and thickness     
195      a_ip_frac(:,:,:) = a_ip(:,:,:) * z1_a_i(:,:,:)
196      WHERE( a_ip_frac(:,:,:) > epsi20 )   ;   h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:)
197      ELSEWHERE                            ;   h_ip(:,:,:) = 0._wp
198      END WHERE
199      !
200      !                                           !---  salinity (with a minimum value imposed everywhere)     
201      IF( nn_icesal == 2 ) THEN
202         WHERE( v_i(:,:,:) > epsi20 )   ;   s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) * z1_v_i(:,:,:) ) )
203         ELSEWHERE                      ;   s_i(:,:,:) = rn_simin
204         END WHERE
205      ENDIF
206      CALL ice_var_salprof   ! salinity profile
207
208      !-------------------
209      ! Ice temperature   [K]   (with a minimum value (rt0 - 100.))
210      !-------------------
211      zlay_i   = REAL( nlay_i , wp )    ! number of layers
212      DO jl = 1, jpl
213         DO jk = 1, nlay_i
214            DO jj = 1, jpj
215               DO ji = 1, jpi
216                  IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area
217                     !
218                     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]
219                     ztmelts          = - sz_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C]
220                     ! Conversion q(S,T) -> T (second order equation)
221                     zbbb             = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus
222                     zccc             = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) )
223                     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
224                     !
225                  ELSE                                !--- no ice
226                     t_i(ji,jj,jk,jl) = rt0
227                  ENDIF
228               END DO
229            END DO
230         END DO
231      END DO
232
233      !--------------------
234      ! Snow temperature   [K]   (with a minimum value (rt0 - 100.))
235      !--------------------
236      zlay_s = REAL( nlay_s , wp )
237      DO jk = 1, nlay_s
238         WHERE( v_s(:,:,:) > epsi20 )        !--- icy area
239            t_s(:,:,jk,:) = MAX( -100._wp , MIN( r1_cpic * ( -r1_rhosn * (e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s) + lfus ) , 0._wp ) ) + rt0
240         ELSEWHERE                           !--- no ice
241            t_s(:,:,jk,:) = rt0
242         END WHERE
243      END DO
244
245      ! integrated values
246      vt_i (:,:) = SUM( v_i, dim=3 )
247      vt_s (:,:) = SUM( v_s, dim=3 )
248      at_i (:,:) = SUM( a_i, dim=3 )
249      !
250   END SUBROUTINE ice_var_glo2eqv
251
252
253   SUBROUTINE ice_var_eqv2glo
254      !!-------------------------------------------------------------------
255      !!                ***  ROUTINE ice_var_eqv2glo ***
256      !!
257      !! ** Purpose :   computes global variables as function of
258      !!              equivalent variables,  i.e. it turns VEQV into VGLO
259      !!-------------------------------------------------------------------
260      !
261      v_i (:,:,:) = h_i(:,:,:) * a_i(:,:,:)
262      v_s (:,:,:) = h_s(:,:,:) * a_i(:,:,:)
263      sv_i(:,:,:) = s_i(:,:,:) * v_i(:,:,:)
264      !
265   END SUBROUTINE ice_var_eqv2glo
266
267
268   SUBROUTINE ice_var_salprof
269      !!-------------------------------------------------------------------
270      !!                ***  ROUTINE ice_var_salprof ***
271      !!
272      !! ** Purpose :   computes salinity profile in function of bulk salinity     
273      !!
274      !! ** Method  : If bulk salinity greater than zsi1,
275      !!              the profile is assumed to be constant (S_inf)
276      !!              If bulk salinity lower than zsi0,
277      !!              the profile is linear with 0 at the surface (S_zero)
278      !!              If it is between zsi0 and zsi1, it is a
279      !!              alpha-weighted linear combination of s_inf and s_zero
280      !!
281      !! ** References : Vancoppenolle et al., 2007
282      !!-------------------------------------------------------------------
283      INTEGER  ::   ji, jj, jk, jl   ! dummy loop index
284      REAL(wp) ::   zsal, z1_dS
285      REAL(wp) ::   zargtemp , zs0, zs
286      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z_slope_s, zalpha    ! case 2 only
287      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
288      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
289      !!-------------------------------------------------------------------
290
291!!gm Question: Remove the option 3 ?  How many years since it last use ?
292
293      SELECT CASE ( nn_icesal )
294      !
295      !               !---------------------------------------!
296      CASE( 1 )       !  constant salinity in time and space  !
297         !            !---------------------------------------!
298         sz_i(:,:,:,:) = rn_icesal
299         s_i(:,:,:)   = rn_icesal
300         !
301         !            !---------------------------------------------!
302      CASE( 2 )       !  time varying salinity with linear profile  !
303         !            !---------------------------------------------!
304         !
305         ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
306         !
307         DO jl = 1, jpl
308            DO jk = 1, nlay_i
309               sz_i(:,:,jk,jl)  = s_i(:,:,jl)
310            END DO
311         END DO
312         !                                      ! Slope of the linear profile
313         WHERE( h_i(:,:,:) > epsi20 )   ;   z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:)
314         ELSEWHERE                       ;   z_slope_s(:,:,:) = 0._wp
315         END WHERE
316         !
317         z1_dS = 1._wp / ( zsi1 - zsi0 )
318         DO jl = 1, jpl
319            DO jj = 1, jpj
320               DO ji = 1, jpi
321                  zalpha(ji,jj,jl) = MAX(  0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp )  )
322                  !                             ! force a constant profile when SSS too low (Baltic Sea)
323                  IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) )   zalpha(ji,jj,jl) = 0._wp 
324               END DO
325            END DO
326         END DO
327
328         ! Computation of the profile
329         DO jl = 1, jpl
330            DO jk = 1, nlay_i
331               DO jj = 1, jpj
332                  DO ji = 1, jpi
333                     !                          ! linear profile with 0 surface value
334                     zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
335                     zs  = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl)     ! weighting the profile
336                     sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
337                  END DO
338               END DO
339            END DO
340         END DO
341         !
342         DEALLOCATE( z_slope_s , zalpha )
343         !
344         !            !-------------------------------------------!
345      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
346         !            !-------------------------------------------!                                   (mean = 2.30)
347         !
348         s_i(:,:,:) = 2.30_wp
349!!gm Remark: if we keep the case 3, then compute an store one for all time-step
350!!           a array  S_prof(1:nlay_i) containing the calculation and just do:
351!         DO jk = 1, nlay_i
352!            sz_i(:,:,jk,:) = S_prof(jk)
353!         END DO
354!!gm end
355         !
356         DO jl = 1, jpl
357            DO jk = 1, nlay_i
358               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
359               sz_i(:,:,jk,jl) =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
360            END DO
361         END DO
362         !
363      END SELECT
364      !
365   END SUBROUTINE ice_var_salprof
366
367   SUBROUTINE ice_var_salprof1d
368      !!-------------------------------------------------------------------
369      !!                  ***  ROUTINE ice_var_salprof1d  ***
370      !!
371      !! ** Purpose :   1d computation of the sea ice salinity profile
372      !!                Works with 1d vectors and is used by thermodynamic modules
373      !!-------------------------------------------------------------------
374      INTEGER  ::   ji, jk    ! dummy loop indices
375      REAL(wp) ::   zargtemp, zsal, z1_dS   ! local scalars
376      REAL(wp) ::   zs, zs0              !   -      -
377      !
378      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z_slope_s, zalpha   !
379      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
380      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
381      !!-------------------------------------------------------------------
382      !
383      SELECT CASE ( nn_icesal )
384      !
385      !               !---------------------------------------!
386      CASE( 1 )       !  constant salinity in time and space  !
387         !            !---------------------------------------!
388         sz_i_1d(1:npti,:) = rn_icesal
389         !
390         !            !---------------------------------------------!
391      CASE( 2 )       !  time varying salinity with linear profile  !
392         !            !---------------------------------------------!
393         !
394         ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
395         !
396         !                                      ! Slope of the linear profile
397         WHERE( h_i_1d(1:npti) > epsi20 )   ;   z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
398         ELSEWHERE                           ;   z_slope_s(1:npti) = 0._wp
399         END WHERE
400         
401         z1_dS = 1._wp / ( zsi1 - zsi0 )
402         DO ji = 1, npti
403            zalpha(ji) = MAX(  0._wp , MIN(  ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp  )  )
404            !                             ! force a constant profile when SSS too low (Baltic Sea)
405            IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) )   zalpha(ji) = 0._wp
406         END DO
407         !
408         ! Computation of the profile
409         DO jk = 1, nlay_i
410            DO ji = 1, npti
411               !                          ! linear profile with 0 surface value
412               zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
413               zs  = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
414               sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
415            END DO
416         END DO
417         !
418         DEALLOCATE( z_slope_s, zalpha )
419
420         !            !-------------------------------------------!
421      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
422         !            !-------------------------------------------!                                   (mean = 2.30)
423         !
424         s_i_1d(1:npti) = 2.30_wp
425         !
426!!gm cf remark in ice_var_salprof routine, CASE( 3 )
427         DO jk = 1, nlay_i
428            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
429            zsal =  1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
430            DO ji = 1, npti
431               sz_i_1d(ji,jk) = zsal
432            END DO
433         END DO
434         !
435      END SELECT
436      !
437   END SUBROUTINE ice_var_salprof1d
438
439
440   SUBROUTINE ice_var_zapsmall
441      !!-------------------------------------------------------------------
442      !!                   ***  ROUTINE ice_var_zapsmall ***
443      !!
444      !! ** Purpose :   Remove too small sea ice areas and correct fluxes
445      !!-------------------------------------------------------------------
446      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
447      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch
448      !!-------------------------------------------------------------------
449      !
450      DO jl = 1, jpl       !==  loop over the categories  ==!
451         !
452         !-----------------------------------------------------------------
453         ! Zap ice energy and use ocean heat to melt ice
454         !-----------------------------------------------------------------
455         WHERE( a_i(:,:,jl) > epsi20 )   ;   h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
456         ELSEWHERE                       ;   h_i(:,:,jl) = 0._wp
457         END WHERE
458         !
459         WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 )   ;   zswitch(:,:) = 0._wp
460         ELSEWHERE                                                                           ;   zswitch(:,:) = 1._wp
461         END WHERE
462         
463         DO jk = 1, nlay_i
464            DO jj = 1 , jpj
465               DO ji = 1 , jpi
466                  ! update exchanges with ocean
467                  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
468                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
469                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
470               END DO
471            END DO
472         END DO
473
474         DO jj = 1 , jpj
475            DO ji = 1 , jpi
476               ! update exchanges with ocean
477               sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoic * r1_rdtice
478               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoic * r1_rdtice
479               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhosn * r1_rdtice
480               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
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            END DO
502         END DO
503
504         IF( ln_pnd ) THEN
505            DO jj = 1 , jpj
506               DO ji = 1 , jpi
507                  IF( ln_pnd_fw )   &
508                     &   wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice
509                  a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
510                  v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
511               END DO
512            END DO
513         ENDIF
514         
515      END DO 
516
517      ! to be sure that at_i is the sum of a_i(jl)
518      at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
519      vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
520
521      ! open water = 1 if at_i=0
522      WHERE( at_i(:,:) == 0._wp )   ato_i(:,:) = 1._wp
523      !
524   END SUBROUTINE ice_var_zapsmall
525
526
527   SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
528      !!-------------------------------------------------------------------
529      !!                ***  ROUTINE ice_var_itd   ***
530      !!
531      !! ** Purpose :  converting 1-cat ice to multiple ice categories
532      !!
533      !!                  ice thickness distribution follows a gaussian law
534      !!               around the concentration of the most likely ice thickness
535      !!                           (similar as iceistate.F90)
536      !!
537      !! ** Method:   Iterative procedure
538      !!               
539      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
540      !!
541      !!               2) Check whether the distribution conserves area and volume, positivity and
542      !!                  category boundaries
543      !!             
544      !!               3) If not (input ice is too thin), the last category is empty and
545      !!                  the number of categories is reduced (jpl-1)
546      !!
547      !!               4) Iterate until ok (SUM(itest(:) = 4)
548      !!
549      !! ** Arguments : zhti: 1-cat ice thickness
550      !!                zhts: 1-cat snow depth
551      !!                zati: 1-cat ice concentration
552      !!
553      !! ** Output    : jpl-cat
554      !!
555      !!  (Example of application: BDY forcings when input are cell averaged) 
556      !!-------------------------------------------------------------------
557      INTEGER  :: ji, jk, jl             ! dummy loop indices
558      INTEGER  :: ijpij, i_fill, jl0 
559      REAL(wp) :: zarg, zV, zconv, zdh, zdv
560      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
561      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i ! output ice/snow variables
562      INTEGER , DIMENSION(4)                  ::   itest
563      !!-------------------------------------------------------------------
564      !
565      ! ----------------------------------------
566      ! distribution over the jpl ice categories
567      ! ----------------------------------------
568      ! a gaussian distribution for ice concentration is used
569      ! then we check whether the distribution fullfills
570      ! volume and area conservation, positivity and ice categories bounds
571      ijpij = SIZE( zhti , 1 )
572      zh_i(1:ijpij,1:jpl) = 0._wp
573      zh_s(1:ijpij,1:jpl) = 0._wp
574      za_i (1:ijpij,1:jpl) = 0._wp
575
576      DO ji = 1, ijpij
577         
578         IF( zhti(ji) > 0._wp ) THEN
579
580            ! find which category (jl0) the input ice thickness falls into
581            jl0 = jpl
582            DO jl = 1, jpl
583               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
584                  jl0 = jl
585                  CYCLE
586               ENDIF
587            END DO
588
589            itest(:) = 0
590            i_fill   = jpl + 1                                            !------------------------------------
591            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
592               !                                                          !------------------------------------
593               i_fill = i_fill - 1
594               !
595               zh_i(ji,1:jpl) = 0._wp
596               za_i (ji,1:jpl) = 0._wp
597               itest(:)        = 0     
598               
599               IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
600                  zh_i(ji,1) = zhti(ji)
601                  za_i (ji,1) = zati (ji)
602               ELSE                         !-- case ice is thicker: fill categories >1
603                  ! thickness
604                  DO jl = 1, i_fill - 1
605                     zh_i(ji,jl) = hi_mean(jl)
606                  END DO
607                 
608                  ! concentration
609                  za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
610                  DO jl = 1, i_fill - 1
611                     IF ( jl /= jl0 ) THEN
612                        zarg        = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
613                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
614                     ENDIF
615                  END DO
616                 
617                  ! last category
618                  za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
619                  zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
620                  zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
621                 
622                  ! clem: correction if concentration of upper cat is greater than lower cat
623                  !       (it should be a gaussian around jl0 but sometimes it is not)
624                  IF ( jl0 /= jpl ) THEN
625                     DO jl = jpl, jl0+1, -1
626                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
627                           zdv = zh_i(ji,jl) * za_i(ji,jl)
628                           zh_i(ji,jl    ) = 0._wp
629                           za_i (ji,jl    ) = 0._wp
630                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
631                        END IF
632                     ENDDO
633                  ENDIF
634               
635               ENDIF
636           
637               ! Compatibility tests
638               zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) ) 
639               IF ( zconv < epsi06 ) itest(1) = 1                                        ! Test 1: area conservation
640           
641               zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )
642               IF ( zconv < epsi06 ) itest(2) = 1                                        ! Test 2: volume conservation
643               
644               IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ?
645               
646               itest(4) = 1
647               DO jl = 1, i_fill
648                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0                                ! Test 4: positivity of ice concentrations
649               END DO
650               !                                         !----------------------------
651            END DO                                       ! end iteration on categories
652               !                                         !----------------------------
653         ENDIF
654      END DO
655
656      ! Add Snow in each category where za_i is not 0
657      DO jl = 1, jpl
658         DO ji = 1, ijpij
659            IF( za_i(ji,jl) > 0._wp ) THEN
660               zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
661               ! In case snow load is in excess that would lead to transformation from snow to ice
662               ! Then, transfer the snow excess into the ice (different from icethd_dh)
663               zdh = MAX( 0._wp, ( rhosn * zh_s(ji,jl) + ( rhoic - rau0 ) * zh_i(ji,jl) ) * r1_rau0 ) 
664               ! recompute h_i, h_s avoiding out of bounds values
665               zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
666               zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoic * r1_rhosn )
667            ENDIF
668         END DO
669      END DO
670      !
671    END SUBROUTINE ice_var_itd
672
673
674    SUBROUTINE ice_var_bv
675      !!-------------------------------------------------------------------
676      !!                ***  ROUTINE ice_var_bv ***
677      !!
678      !! ** Purpose :   computes mean brine volume (%) in sea ice
679      !!
680      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
681      !!
682      !! References : Vancoppenolle et al., JGR, 2007
683      !!-------------------------------------------------------------------
684      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
685      !!-------------------------------------------------------------------
686      !
687!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
688!!   instead of setting everything to zero as just below
689      bv_i (:,:,:) = 0._wp
690      DO jl = 1, jpl
691         DO jk = 1, nlay_i
692            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
693               bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
694            END WHERE
695         END DO
696      END DO
697      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
698      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
699      END WHERE
700      !
701   END SUBROUTINE ice_var_bv
702
703
704#else
705   !!----------------------------------------------------------------------
706   !!   Default option         Dummy module           NO ESIM sea-ice model
707   !!----------------------------------------------------------------------
708#endif
709
710   !!======================================================================
711END MODULE icevar
Note: See TracBrowser for help on using the repository browser.