New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icevar.F90 in NEMO/trunk/src/ICE – NEMO

source: NEMO/trunk/src/ICE/icevar.F90 @ 9653

Last change on this file since 9653 was 9604, checked in by clem, 6 years ago

change history of the ice routines

File size: 39.2 KB
Line 
1MODULE icevar
2   !!======================================================================
3   !!                       ***  MODULE icevar ***
4   !!   sea-ice:  series of functions to transform or compute ice variables
5   !!======================================================================
6   !! History :   -   !  2006-01  (M. Vancoppenolle) Original code
7   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
8   !!----------------------------------------------------------------------
9#if defined key_si3
10   !!----------------------------------------------------------------------
11   !!   'key_si3'                                       SI3 sea-ice model
12   !!----------------------------------------------------------------------
13   !!
14   !!                 There are three sets of variables
15   !!                 VGLO : global variables of the model
16   !!                        - v_i (jpi,jpj,jpl)
17   !!                        - v_s (jpi,jpj,jpl)
18   !!                        - a_i (jpi,jpj,jpl)
19   !!                        - t_s (jpi,jpj,jpl)
20   !!                        - e_i (jpi,jpj,nlay_i,jpl)
21   !!                        - e_s (jpi,jpj,nlay_s,jpl)
22   !!                        - sv_i(jpi,jpj,jpl)
23   !!                        - oa_i(jpi,jpj,jpl)
24   !!                 VEQV : equivalent variables sometimes used in the model
25   !!                        - h_i(jpi,jpj,jpl)
26   !!                        - h_s(jpi,jpj,jpl)
27   !!                        - t_i(jpi,jpj,nlay_i,jpl)
28   !!                        ...
29   !!                 VAGG : aggregate variables, averaged/summed over all
30   !!                        thickness categories
31   !!                        - vt_i(jpi,jpj)
32   !!                        - vt_s(jpi,jpj)
33   !!                        - at_i(jpi,jpj)
34   !!                        - et_s(jpi,jpj)  total snow heat content
35   !!                        - et_i(jpi,jpj)  total ice thermal content
36   !!                        - sm_i(jpi,jpj)  mean ice salinity
37   !!                        - tm_i(jpi,jpj)  mean ice temperature
38   !!                        - tm_s(jpi,jpj)  mean snw temperature
39   !!----------------------------------------------------------------------
40   !!   ice_var_agg       : integrate variables over layers and categories
41   !!   ice_var_glo2eqv   : transform from VGLO to VEQV
42   !!   ice_var_eqv2glo   : transform from VEQV to VGLO
43   !!   ice_var_salprof   : salinity profile in the ice
44   !!   ice_var_salprof1d : salinity profile in the ice 1D
45   !!   ice_var_zapsmall  : remove very small area and volume
46   !!   ice_var_itd       : convert 1-cat to jpl-cat
47   !!   ice_var_itd2      : convert N-cat to jpl-cat
48   !!   ice_var_bv        : brine volume
49   !!   ice_var_enthalpy  : compute ice and snow enthalpies from temperature
50   !!----------------------------------------------------------------------
51   USE dom_oce        ! ocean space and time domain
52   USE phycst         ! physical constants (ocean directory)
53   USE sbc_oce , ONLY : sss_m
54   USE ice            ! sea-ice: variables
55   USE ice1D          ! sea-ice: thermodynamics variables
56   !
57   USE in_out_manager ! I/O manager
58   USE lib_mpp        ! MPP library
59   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
60
61   IMPLICIT NONE
62   PRIVATE
63
64   PUBLIC   ice_var_agg         
65   PUBLIC   ice_var_glo2eqv     
66   PUBLIC   ice_var_eqv2glo     
67   PUBLIC   ice_var_salprof     
68   PUBLIC   ice_var_salprof1d   
69   PUBLIC   ice_var_zapsmall
70   PUBLIC   ice_var_itd
71   PUBLIC   ice_var_itd2
72   PUBLIC   ice_var_bv           
73   PUBLIC   ice_var_enthalpy           
74
75   !!----------------------------------------------------------------------
76   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
77   !! $Id: icevar.F90 8422 2017-08-08 13:58:05Z clem $
78   !! Software governed by the CeCILL licence     (./LICENSE)
79   !!----------------------------------------------------------------------
80CONTAINS
81
82   SUBROUTINE ice_var_agg( kn )
83      !!-------------------------------------------------------------------
84      !!                ***  ROUTINE ice_var_agg  ***
85      !!
86      !! ** Purpose :   aggregates ice-thickness-category variables to
87      !!              all-ice variables, i.e. it turns VGLO into VAGG
88      !!-------------------------------------------------------------------
89      INTEGER, INTENT( in ) ::   kn     ! =1 state variables only
90      !                                 ! >1 state variables + others
91      !
92      INTEGER ::   ji, jj, jk, jl   ! dummy loop indices
93      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z1_at_i, z1_vt_i, z1_vt_s
94      !!-------------------------------------------------------------------
95      !
96      !                                      ! integrated values
97      vt_i(:,:) =       SUM( v_i(:,:,:)           , dim=3 )
98      vt_s(:,:) =       SUM( v_s(:,:,:)           , dim=3 )
99      at_i(:,:) =       SUM( a_i(:,:,:)           , dim=3 )
100      et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 )
101      et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
102      !
103      at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds
104      vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
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) , z1_vt_s(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         WHERE( vt_s(:,:) > epsi20 )   ;   z1_vt_s(:,:) = 1._wp / vt_s(:,:)
118         ELSEWHERE                     ;   z1_vt_s(:,:) = 0._wp
119         END WHERE
120         !
121         !                          ! mean ice/snow thickness
122         hm_i(:,:) = vt_i(:,:) * z1_at_i(:,:)
123         hm_s(:,:) = vt_s(:,:) * z1_at_i(:,:)
124         !         
125         !                          ! mean temperature (K), salinity and age
126         tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
127         tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
128         om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(:,:)
129         sm_i (:,:) = SUM( sv_i(:,:,:)              , dim=3 ) * z1_vt_i(:,:)
130         !
131         tm_i(:,:) = 0._wp
132         tm_s(:,:) = 0._wp
133         DO jl = 1, jpl
134            DO jk = 1, nlay_i
135               tm_i(:,:) = tm_i(:,:) + r1_nlay_i * t_i (:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
136            END DO
137            DO jk = 1, nlay_s
138               tm_s(:,:) = tm_s(:,:) + r1_nlay_s * t_s (:,:,jk,jl) * v_s(:,:,jl) * z1_vt_s(:,:)
139            END DO
140         END DO
141         !
142         !                           ! put rt0 where there is no ice
143         WHERE( at_i(:,:)<=epsi20 )
144            tm_su(:,:) = rt0
145            tm_si(:,:) = rt0
146            tm_i (:,:) = rt0
147            tm_s (:,:) = rt0
148         END WHERE
149
150         DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s )
151      ENDIF
152      !
153   END SUBROUTINE ice_var_agg
154
155
156   SUBROUTINE ice_var_glo2eqv
157      !!-------------------------------------------------------------------
158      !!                ***  ROUTINE ice_var_glo2eqv ***
159      !!
160      !! ** Purpose :   computes equivalent variables as function of 
161      !!              global variables, i.e. it turns VGLO into VEQV
162      !!-------------------------------------------------------------------
163      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
164      REAL(wp) ::   ze_i             ! local scalars
165      REAL(wp) ::   ze_s, ztmelts, zbbb, zccc       !   -      -
166      REAL(wp) ::   zhmax, z1_zhmax                 !   -      -
167      REAL(wp) ::   zlay_i, zlay_s                  !   -      -
168      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i
169      !!-------------------------------------------------------------------
170
171!!gm Question 2:  It is possible to define existence of sea-ice in a common way between
172!!                ice area and ice volume ?
173!!                the idea is to be able to define one for all at the begining of this routine
174!!                a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
175
176      !---------------------------------------------------------------
177      ! Ice thickness, snow thickness, ice salinity, ice age and ponds
178      !---------------------------------------------------------------
179      !                                            !--- inverse of the ice area
180      WHERE( a_i(:,:,:) > epsi20 )   ;   z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
181      ELSEWHERE                      ;   z1_a_i(:,:,:) = 0._wp
182      END WHERE
183      !
184      WHERE( v_i(:,:,:) > epsi20 )   ;   z1_v_i(:,:,:) = 1._wp / v_i(:,:,:)
185      ELSEWHERE                      ;   z1_v_i(:,:,:) = 0._wp
186      END WHERE
187      !                                           !--- ice thickness
188      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
189
190      zhmax    =          hi_max(jpl)
191      z1_zhmax =  1._wp / hi_max(jpl)               
192      WHERE( h_i(:,:,jpl) > zhmax )   ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area
193         h_i  (:,:,jpl) = zhmax
194         a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 
195         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)
196      END WHERE
197      !                                           !--- snow thickness
198      h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:)
199      !                                           !--- ice age     
200      o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:)
201      !                                           !--- pond fraction and thickness     
202      a_ip_frac(:,:,:) = a_ip(:,:,:) * z1_a_i(:,:,:)
203      WHERE( a_ip_frac(:,:,:) > epsi20 )   ;   h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:)
204      ELSEWHERE                            ;   h_ip(:,:,:) = 0._wp
205      END WHERE
206      !
207      !                                           !---  salinity (with a minimum value imposed everywhere)     
208      IF( nn_icesal == 2 ) THEN
209         WHERE( v_i(:,:,:) > epsi20 )   ;   s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) * z1_v_i(:,:,:) ) )
210         ELSEWHERE                      ;   s_i(:,:,:) = rn_simin
211         END WHERE
212      ENDIF
213      CALL ice_var_salprof   ! salinity profile
214
215      !-------------------
216      ! Ice temperature   [K]   (with a minimum value (rt0 - 100.))
217      !-------------------
218      zlay_i   = REAL( nlay_i , wp )    ! number of layers
219      DO jl = 1, jpl
220         DO jk = 1, nlay_i
221            DO jj = 1, jpj
222               DO ji = 1, jpi
223                  IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area
224                     !
225                     ze_i             =   e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3]
226                     ztmelts          = - sz_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C]
227                     ! Conversion q(S,T) -> T (second order equation)
228                     zbbb             = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus
229                     zccc             = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) )
230                     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
231                     !
232                  ELSE                                !--- no ice
233                     t_i(ji,jj,jk,jl) = rt0
234                  ENDIF
235               END DO
236            END DO
237         END DO
238      END DO
239
240      !--------------------
241      ! Snow temperature   [K]   (with a minimum value (rt0 - 100.))
242      !--------------------
243      zlay_s = REAL( nlay_s , wp )
244      DO jk = 1, nlay_s
245         WHERE( v_s(:,:,:) > epsi20 )        !--- icy area
246            t_s(:,:,jk,:) = rt0 + MAX( -100._wp ,  &
247                 &                MIN( r1_cpic * ( -r1_rhosn * ( e_s(:,:,jk,:) / v_s(:,:,:) * zlay_s ) + lfus ) , 0._wp ) )
248         ELSEWHERE                           !--- no ice
249            t_s(:,:,jk,:) = rt0
250         END WHERE
251      END DO
252      !
253      ! integrated values
254      vt_i (:,:) = SUM( v_i, dim=3 )
255      vt_s (:,:) = SUM( v_s, dim=3 )
256      at_i (:,:) = SUM( a_i, dim=3 )
257      !
258   END SUBROUTINE ice_var_glo2eqv
259
260
261   SUBROUTINE ice_var_eqv2glo
262      !!-------------------------------------------------------------------
263      !!                ***  ROUTINE ice_var_eqv2glo ***
264      !!
265      !! ** Purpose :   computes global variables as function of
266      !!              equivalent variables,  i.e. it turns VEQV into VGLO
267      !!-------------------------------------------------------------------
268      !
269      v_i (:,:,:) = h_i (:,:,:) * a_i (:,:,:)
270      v_s (:,:,:) = h_s (:,:,:) * a_i (:,:,:)
271      sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:)
272      v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
273      !
274   END SUBROUTINE ice_var_eqv2glo
275
276
277   SUBROUTINE ice_var_salprof
278      !!-------------------------------------------------------------------
279      !!                ***  ROUTINE ice_var_salprof ***
280      !!
281      !! ** Purpose :   computes salinity profile in function of bulk salinity     
282      !!
283      !! ** Method  : If bulk salinity greater than zsi1,
284      !!              the profile is assumed to be constant (S_inf)
285      !!              If bulk salinity lower than zsi0,
286      !!              the profile is linear with 0 at the surface (S_zero)
287      !!              If it is between zsi0 and zsi1, it is a
288      !!              alpha-weighted linear combination of s_inf and s_zero
289      !!
290      !! ** References : Vancoppenolle et al., 2007
291      !!-------------------------------------------------------------------
292      INTEGER  ::   ji, jj, jk, jl   ! dummy loop index
293      REAL(wp) ::   zsal, z1_dS
294      REAL(wp) ::   zargtemp , zs0, zs
295      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z_slope_s, zalpha    ! case 2 only
296      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
297      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
298      !!-------------------------------------------------------------------
299
300!!gm Question: Remove the option 3 ?  How many years since it last use ?
301
302      SELECT CASE ( nn_icesal )
303      !
304      !               !---------------------------------------!
305      CASE( 1 )       !  constant salinity in time and space  !
306         !            !---------------------------------------!
307         sz_i(:,:,:,:) = rn_icesal
308         s_i (:,:,:)   = rn_icesal
309         !
310         !            !---------------------------------------------!
311      CASE( 2 )       !  time varying salinity with linear profile  !
312         !            !---------------------------------------------!
313         !
314         ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
315         !
316         DO jl = 1, jpl
317            DO jk = 1, nlay_i
318               sz_i(:,:,jk,jl)  = s_i(:,:,jl)
319            END DO
320         END DO
321         !                                      ! Slope of the linear profile
322         WHERE( h_i(:,:,:) > epsi20 )   ;   z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:)
323         ELSEWHERE                      ;   z_slope_s(:,:,:) = 0._wp
324         END WHERE
325         !
326         z1_dS = 1._wp / ( zsi1 - zsi0 )
327         DO jl = 1, jpl
328            DO jj = 1, jpj
329               DO ji = 1, jpi
330                  zalpha(ji,jj,jl) = MAX(  0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp )  )
331                  !                             ! force a constant profile when SSS too low (Baltic Sea)
332                  IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) )   zalpha(ji,jj,jl) = 0._wp 
333               END DO
334            END DO
335         END DO
336         !
337         ! Computation of the profile
338         DO jl = 1, jpl
339            DO jk = 1, nlay_i
340               DO jj = 1, jpj
341                  DO ji = 1, jpi
342                     !                          ! linear profile with 0 surface value
343                     zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
344                     zs  = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl)     ! weighting the profile
345                     sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
346                  END DO
347               END DO
348            END DO
349         END DO
350         !
351         DEALLOCATE( z_slope_s , zalpha )
352         !
353         !            !-------------------------------------------!
354      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
355         !            !-------------------------------------------!                                   (mean = 2.30)
356         !
357         s_i(:,:,:) = 2.30_wp
358!!gm Remark: if we keep the case 3, then compute an store one for all time-step
359!!           a array  S_prof(1:nlay_i) containing the calculation and just do:
360!         DO jk = 1, nlay_i
361!            sz_i(:,:,jk,:) = S_prof(jk)
362!         END DO
363!!gm end
364         !
365         DO jl = 1, jpl
366            DO jk = 1, nlay_i
367               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
368               sz_i(:,:,jk,jl) =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
369            END DO
370         END DO
371         !
372      END SELECT
373      !
374   END SUBROUTINE ice_var_salprof
375
376
377   SUBROUTINE ice_var_salprof1d
378      !!-------------------------------------------------------------------
379      !!                  ***  ROUTINE ice_var_salprof1d  ***
380      !!
381      !! ** Purpose :   1d computation of the sea ice salinity profile
382      !!                Works with 1d vectors and is used by thermodynamic modules
383      !!-------------------------------------------------------------------
384      INTEGER  ::   ji, jk    ! dummy loop indices
385      REAL(wp) ::   zargtemp, zsal, z1_dS   ! local scalars
386      REAL(wp) ::   zs, zs0              !   -      -
387      !
388      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z_slope_s, zalpha   !
389      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
390      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
391      !!-------------------------------------------------------------------
392      !
393      SELECT CASE ( nn_icesal )
394      !
395      !               !---------------------------------------!
396      CASE( 1 )       !  constant salinity in time and space  !
397         !            !---------------------------------------!
398         sz_i_1d(1:npti,:) = rn_icesal
399         !
400         !            !---------------------------------------------!
401      CASE( 2 )       !  time varying salinity with linear profile  !
402         !            !---------------------------------------------!
403         !
404         ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
405         !
406         !                                      ! Slope of the linear profile
407         WHERE( h_i_1d(1:npti) > epsi20 )   ;   z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
408         ELSEWHERE                          ;   z_slope_s(1:npti) = 0._wp
409         END WHERE
410         
411         z1_dS = 1._wp / ( zsi1 - zsi0 )
412         DO ji = 1, npti
413            zalpha(ji) = MAX(  0._wp , MIN(  ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp  )  )
414            !                             ! force a constant profile when SSS too low (Baltic Sea)
415            IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) )   zalpha(ji) = 0._wp
416         END DO
417         !
418         ! Computation of the profile
419         DO jk = 1, nlay_i
420            DO ji = 1, npti
421               !                          ! linear profile with 0 surface value
422               zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
423               zs  = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
424               sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
425            END DO
426         END DO
427         !
428         DEALLOCATE( z_slope_s, zalpha )
429
430         !            !-------------------------------------------!
431      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
432         !            !-------------------------------------------!                                   (mean = 2.30)
433         !
434         s_i_1d(1:npti) = 2.30_wp
435         !
436!!gm cf remark in ice_var_salprof routine, CASE( 3 )
437         DO jk = 1, nlay_i
438            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
439            zsal =  1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
440            DO ji = 1, npti
441               sz_i_1d(ji,jk) = zsal
442            END DO
443         END DO
444         !
445      END SELECT
446      !
447   END SUBROUTINE ice_var_salprof1d
448
449
450   SUBROUTINE ice_var_zapsmall
451      !!-------------------------------------------------------------------
452      !!                   ***  ROUTINE ice_var_zapsmall ***
453      !!
454      !! ** Purpose :   Remove too small sea ice areas and correct fluxes
455      !!-------------------------------------------------------------------
456      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
457      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch
458      !!-------------------------------------------------------------------
459      !
460      DO jl = 1, jpl       !==  loop over the categories  ==!
461         !
462         !-----------------------------------------------------------------
463         ! Zap ice energy and use ocean heat to melt ice
464         !-----------------------------------------------------------------
465         WHERE( a_i(:,:,jl) > epsi10 )   ;   h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
466         ELSEWHERE                       ;   h_i(:,:,jl) = 0._wp
467         END WHERE
468         !
469         WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 )   ;   zswitch(:,:) = 0._wp
470         ELSEWHERE                                                                           ;   zswitch(:,:) = 1._wp
471         END WHERE
472         !
473         DO jk = 1, nlay_i
474            DO jj = 1 , jpj
475               DO ji = 1 , jpi
476                  ! update exchanges with ocean
477                  hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
478                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
479                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
480               END DO
481            END DO
482         END DO
483         !
484         DO jk = 1, nlay_s
485            DO jj = 1 , jpj
486               DO ji = 1 , jpi
487                  ! update exchanges with ocean
488                  hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
489                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj)
490                  t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
491               END DO
492            END DO
493         END DO
494         !
495         DO jj = 1 , jpj
496            DO ji = 1 , jpi
497               ! update exchanges with ocean
498               sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoic * r1_rdtice
499               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoic * r1_rdtice
500               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhosn * r1_rdtice
501               !
502               !-----------------------------------------------------------------
503               ! zap ice and snow volume, add water and salt to ocean
504               !-----------------------------------------------------------------
505               a_i  (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
506               v_i  (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
507               v_s  (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
508               t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
509               oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
510               sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
511               !
512               h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
513               h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
514               !
515               a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
516               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
517               !
518            END DO
519         END DO
520         !
521      END DO 
522
523      ! to be sure that at_i is the sum of a_i(jl)
524      at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
525      vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
526
527      ! open water = 1 if at_i=0
528      WHERE( at_i(:,:) == 0._wp )   ato_i(:,:) = 1._wp
529      !
530   END SUBROUTINE ice_var_zapsmall
531
532
533   SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
534      !!-------------------------------------------------------------------
535      !!                ***  ROUTINE ice_var_itd   ***
536      !!
537      !! ** Purpose :  converting 1-cat ice to multiple ice categories
538      !!
539      !!                  ice thickness distribution follows a gaussian law
540      !!               around the concentration of the most likely ice thickness
541      !!                           (similar as iceistate.F90)
542      !!
543      !! ** Method:   Iterative procedure
544      !!               
545      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
546      !!
547      !!               2) Check whether the distribution conserves area and volume, positivity and
548      !!                  category boundaries
549      !!             
550      !!               3) If not (input ice is too thin), the last category is empty and
551      !!                  the number of categories is reduced (jpl-1)
552      !!
553      !!               4) Iterate until ok (SUM(itest(:) = 4)
554      !!
555      !! ** Arguments : zhti: 1-cat ice thickness
556      !!                zhts: 1-cat snow depth
557      !!                zati: 1-cat ice concentration
558      !!
559      !! ** Output    : jpl-cat
560      !!
561      !!  (Example of application: BDY forcings when input are cell averaged) 
562      !!-------------------------------------------------------------------
563      INTEGER  :: ji, jk, jl             ! dummy loop indices
564      INTEGER  :: idim, i_fill, jl0 
565      REAL(wp) :: zarg, zV, zconv, zdh, zdv
566      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
567      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i ! output ice/snow variables
568      INTEGER , DIMENSION(4)                  ::   itest
569      !!-------------------------------------------------------------------
570      !
571      ! ----------------------------------------
572      ! distribution over the jpl ice categories
573      ! ----------------------------------------
574      ! a gaussian distribution for ice concentration is used
575      ! then we check whether the distribution fullfills
576      ! volume and area conservation, positivity and ice categories bounds
577      idim = SIZE( zhti , 1 )
578      zh_i(1:idim,1:jpl) = 0._wp
579      zh_s(1:idim,1:jpl) = 0._wp
580      za_i(1:idim,1:jpl) = 0._wp
581      !
582      DO ji = 1, idim
583         !
584         IF( zhti(ji) > 0._wp ) THEN
585            !
586            ! find which category (jl0) the input ice thickness falls into
587            jl0 = jpl
588            DO jl = 1, jpl
589               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
590                  jl0 = jl
591                  CYCLE
592               ENDIF
593            END DO
594            !
595            itest(:) = 0
596            i_fill   = jpl + 1                                            !------------------------------------
597            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
598               !                                                          !------------------------------------
599               i_fill = i_fill - 1
600               !
601               zh_i(ji,1:jpl) = 0._wp
602               za_i(ji,1:jpl) = 0._wp
603               itest(:)       = 0     
604               !
605               IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
606                  zh_i(ji,1) = zhti(ji)
607                  za_i (ji,1) = zati (ji)
608               ELSE                         !-- case ice is thicker: fill categories >1
609                  ! thickness
610                  DO jl = 1, i_fill - 1
611                     zh_i(ji,jl) = hi_mean(jl)
612                  END DO
613                  !
614                  ! concentration
615                  za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
616                  DO jl = 1, i_fill - 1
617                     IF ( jl /= jl0 ) THEN
618                        zarg        = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
619                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
620                     ENDIF
621                  END DO
622                  !
623                  ! last category
624                  za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
625                  zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
626                  zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
627                  !
628                  ! correction if concentration of upper cat is greater than lower cat
629                  !    (it should be a gaussian around jl0 but sometimes it is not)
630                  IF ( jl0 /= jpl ) THEN
631                     DO jl = jpl, jl0+1, -1
632                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
633                           zdv = zh_i(ji,jl) * za_i(ji,jl)
634                           zh_i(ji,jl    ) = 0._wp
635                           za_i (ji,jl    ) = 0._wp
636                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
637                        END IF
638                     END DO
639                  ENDIF
640                  !
641               ENDIF
642               !
643               ! Compatibility tests
644               zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) ) 
645               IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation
646               !
647               zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )
648               IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation
649               !
650               IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ?
651               !
652               itest(4) = 1
653               DO jl = 1, i_fill
654                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0                                ! Test 4: positivity of ice concentrations
655               END DO
656               !                                         !----------------------------
657            END DO                                       ! end iteration on categories
658            !                                            !----------------------------
659         ENDIF
660      END DO
661
662      ! Add Snow in each category where za_i is not 0
663      DO jl = 1, jpl
664         DO ji = 1, idim
665            IF( za_i(ji,jl) > 0._wp ) THEN
666               zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
667               ! In case snow load is in excess that would lead to transformation from snow to ice
668               ! Then, transfer the snow excess into the ice (different from icethd_dh)
669               zdh = MAX( 0._wp, ( rhosn * zh_s(ji,jl) + ( rhoic - rau0 ) * zh_i(ji,jl) ) * r1_rau0 ) 
670               ! recompute h_i, h_s avoiding out of bounds values
671               zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
672               zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoic * r1_rhosn )
673            ENDIF
674         END DO
675      END DO
676      !
677   END SUBROUTINE ice_var_itd
678
679
680   SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )
681      !!-------------------------------------------------------------------
682      !!                ***  ROUTINE ice_var_itd2   ***
683      !!
684      !! ** Purpose :  converting N-cat ice to jpl ice categories
685      !!
686      !!                  ice thickness distribution follows a gaussian law
687      !!               around the concentration of the most likely ice thickness
688      !!                           (similar as iceistate.F90)
689      !!
690      !! ** Method:   Iterative procedure
691      !!               
692      !!               1) Fill ice cat that correspond to input thicknesses
693      !!                  Find the lowest(jlmin) and highest(jlmax) cat that are filled
694      !!
695      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1
696      !!                   by removing 25% ice area from jlmin and jlmax (resp.)
697      !!             
698      !!               3) Expand the filling to the empty cat between jlmin and jlmax
699      !!                   by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
700      !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
701      !!
702      !! ** Arguments : zhti: N-cat ice thickness
703      !!                zhts: N-cat snow depth
704      !!                zati: N-cat ice concentration
705      !!
706      !! ** Output    : jpl-cat
707      !!
708      !!  (Example of application: BDY forcings when inputs have N-cat /= jpl) 
709      !!-------------------------------------------------------------------
710      INTEGER  ::   ji, jl, jl1, jl2             ! dummy loop indices
711      INTEGER  ::   idim, icat 
712      INTEGER, PARAMETER ::   ztrans = 0.25_wp
713      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
714      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
715      INTEGER , DIMENSION(:,:), ALLOCATABLE   ::   jlfil, jlfil2
716      INTEGER , DIMENSION(:)  , ALLOCATABLE   ::   jlmax, jlmin
717      !!-------------------------------------------------------------------
718      !
719      idim = SIZE( zhti, 1 )
720      icat = SIZE( zhti, 2 )
721      !
722      ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays
723      ALLOCATE( jlmin(idim), jlmax(idim) )
724
725      ! --- initialize output fields to 0 --- !
726      zh_i(1:idim,1:jpl) = 0._wp
727      zh_s(1:idim,1:jpl) = 0._wp
728      za_i(1:idim,1:jpl) = 0._wp
729      !
730      ! --- fill the categories --- !
731      !     find where cat-input = cat-output and fill cat-output fields 
732      jlmax(:) = 0
733      jlmin(:) = 999
734      jlfil(:,:) = 0
735      DO jl1 = 1, jpl
736         DO jl2 = 1, icat
737            DO ji = 1, idim
738               IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
739                  ! fill the right category
740                  zh_i(ji,jl1) = zhti(ji,jl2)
741                  zh_s(ji,jl1) = zhts(ji,jl2)
742                  za_i(ji,jl1) = zati(ji,jl2)
743                  ! record categories that are filled
744                  jlmax(ji) = MAX( jlmax(ji), jl1 )
745                  jlmin(ji) = MIN( jlmin(ji), jl1 )
746                  jlfil(ji,jl1) = jl1
747               ENDIF
748            END DO
749         END DO
750      END DO
751      !
752      ! --- fill the gaps between categories --- ! 
753      !     transfer from categories filled at the previous step to the empty ones in between
754      DO ji = 1, idim
755         jl1 = jlmin(ji)
756         jl2 = jlmax(ji)
757         IF( jl1 > 1 ) THEN
758            ! fill the lower cat (jl1-1)
759            za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
760            zh_i(ji,jl1-1) = hi_mean(jl1-1)
761            ! remove from cat jl1
762            za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
763         ENDIF
764         IF( jl2 < jpl ) THEN
765            ! fill the upper cat (jl2+1)
766            za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
767            zh_i(ji,jl2+1) = hi_mean(jl2+1)
768            ! remove from cat jl2
769            za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
770         ENDIF
771      END DO
772      !
773      jlfil2(:,:) = jlfil(:,:) 
774      ! fill categories from low to high
775      DO jl = 2, jpl-1
776         DO ji = 1, idim
777            IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
778               ! fill high
779               za_i(ji,jl) = ztrans * za_i(ji,jl-1)
780               zh_i(ji,jl) = hi_mean(jl)
781               jlfil(ji,jl) = jl
782               ! remove low
783               za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
784            ENDIF
785         END DO
786      END DO
787      !
788      ! fill categories from high to low
789      DO jl = jpl-1, 2, -1
790         DO ji = 1, idim
791            IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
792               ! fill low
793               za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
794               zh_i(ji,jl) = hi_mean(jl) 
795               jlfil2(ji,jl) = jl
796               ! remove high
797               za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
798            ENDIF
799         END DO
800      END DO
801      !
802      DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
803      DEALLOCATE( jlmin, jlmax )
804      !
805   END SUBROUTINE ice_var_itd2
806
807
808   SUBROUTINE ice_var_bv
809      !!-------------------------------------------------------------------
810      !!                ***  ROUTINE ice_var_bv ***
811      !!
812      !! ** Purpose :   computes mean brine volume (%) in sea ice
813      !!
814      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
815      !!
816      !! References : Vancoppenolle et al., JGR, 2007
817      !!-------------------------------------------------------------------
818      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
819      !!-------------------------------------------------------------------
820      !
821!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
822!!   instead of setting everything to zero as just below
823      bv_i (:,:,:) = 0._wp
824      DO jl = 1, jpl
825         DO jk = 1, nlay_i
826            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
827               bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
828            END WHERE
829         END DO
830      END DO
831      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
832      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
833      END WHERE
834      !
835   END SUBROUTINE ice_var_bv
836
837
838   SUBROUTINE ice_var_enthalpy
839      !!-------------------------------------------------------------------
840      !!                   ***  ROUTINE ice_var_enthalpy ***
841      !!                 
842      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
843      !!
844      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
845      !!-------------------------------------------------------------------
846      INTEGER  ::   ji, jk   ! dummy loop indices
847      REAL(wp) ::   ztmelts  ! local scalar
848      !!-------------------------------------------------------------------
849      !
850      DO jk = 1, nlay_i             ! Sea ice energy of melting
851         DO ji = 1, npti
852            ztmelts      = - tmut  * sz_i_1d(ji,jk)
853            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point
854                                                                !   (sometimes zdf scheme produces abnormally high temperatures)   
855            e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
856               &                    + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
857               &                    - rcp  *   ztmelts )
858         END DO
859      END DO
860      DO jk = 1, nlay_s             ! Snow energy of melting
861         DO ji = 1, npti
862            e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
863         END DO
864      END DO
865      !
866   END SUBROUTINE ice_var_enthalpy
867
868#else
869   !!----------------------------------------------------------------------
870   !!   Default option         Dummy module           NO SI3 sea-ice model
871   !!----------------------------------------------------------------------
872#endif
873
874   !!======================================================================
875END MODULE icevar
Note: See TracBrowser for help on using the repository browser.