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 branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90 @ 9118

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

debug ice test cases with thermo or dynamics only

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