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

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

changes in style - part5 - almost done

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