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

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

changes in style - part6 - pure cosmetics

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