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

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

almost useless commits

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