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

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

changes in style - part6 - one more round

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