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

Last change on this file since 9119 was 9119, checked in by nicolasmartin, 3 years ago

Fix longer lines so should be harmless (passed SETTE compilations)

File size: 39.1 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                                                                          ,   &
239                 &               MIN( r1_cpic * ( -r1_rhosn * ( e_s(:,:,jk,:) / v_s(:,:,:) * zlay_s ) + lfus ) ,       &
240                 &                    0._wp                                                                      )   ) &
241                 &        + rt0
242         ELSEWHERE                           !--- no ice
243            t_s(:,:,jk,:) = rt0
244         END WHERE
245      END DO
246      !
247      ! integrated values
248      vt_i (:,:) = SUM( v_i, dim=3 )
249      vt_s (:,:) = SUM( v_s, dim=3 )
250      at_i (:,:) = SUM( a_i, dim=3 )
251      !
252   END SUBROUTINE ice_var_glo2eqv
253
254
255   SUBROUTINE ice_var_eqv2glo
256      !!-------------------------------------------------------------------
257      !!                ***  ROUTINE ice_var_eqv2glo ***
258      !!
259      !! ** Purpose :   computes global variables as function of
260      !!              equivalent variables,  i.e. it turns VEQV into VGLO
261      !!-------------------------------------------------------------------
262      !
263      v_i (:,:,:) = h_i (:,:,:) * a_i (:,:,:)
264      v_s (:,:,:) = h_s (:,:,:) * a_i (:,:,:)
265      sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:)
266      v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
267      !
268   END SUBROUTINE ice_var_eqv2glo
269
270
271   SUBROUTINE ice_var_salprof
272      !!-------------------------------------------------------------------
273      !!                ***  ROUTINE ice_var_salprof ***
274      !!
275      !! ** Purpose :   computes salinity profile in function of bulk salinity     
276      !!
277      !! ** Method  : If bulk salinity greater than zsi1,
278      !!              the profile is assumed to be constant (S_inf)
279      !!              If bulk salinity lower than zsi0,
280      !!              the profile is linear with 0 at the surface (S_zero)
281      !!              If it is between zsi0 and zsi1, it is a
282      !!              alpha-weighted linear combination of s_inf and s_zero
283      !!
284      !! ** References : Vancoppenolle et al., 2007
285      !!-------------------------------------------------------------------
286      INTEGER  ::   ji, jj, jk, jl   ! dummy loop index
287      REAL(wp) ::   zsal, z1_dS
288      REAL(wp) ::   zargtemp , zs0, zs
289      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z_slope_s, zalpha    ! case 2 only
290      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
291      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
292      !!-------------------------------------------------------------------
293
294!!gm Question: Remove the option 3 ?  How many years since it last use ?
295
296      SELECT CASE ( nn_icesal )
297      !
298      !               !---------------------------------------!
299      CASE( 1 )       !  constant salinity in time and space  !
300         !            !---------------------------------------!
301         sz_i(:,:,:,:) = rn_icesal
302         s_i (:,:,:)   = rn_icesal
303         !
304         !            !---------------------------------------------!
305      CASE( 2 )       !  time varying salinity with linear profile  !
306         !            !---------------------------------------------!
307         !
308         ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
309         !
310         DO jl = 1, jpl
311            DO jk = 1, nlay_i
312               sz_i(:,:,jk,jl)  = s_i(:,:,jl)
313            END DO
314         END DO
315         !                                      ! Slope of the linear profile
316         WHERE( h_i(:,:,:) > epsi20 )   ;   z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:)
317         ELSEWHERE                      ;   z_slope_s(:,:,:) = 0._wp
318         END WHERE
319         !
320         z1_dS = 1._wp / ( zsi1 - zsi0 )
321         DO jl = 1, jpl
322            DO jj = 1, jpj
323               DO ji = 1, jpi
324                  zalpha(ji,jj,jl) = MAX(  0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp )  )
325                  !                             ! force a constant profile when SSS too low (Baltic Sea)
326                  IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) )   zalpha(ji,jj,jl) = 0._wp 
327               END DO
328            END DO
329         END DO
330         !
331         ! Computation of the profile
332         DO jl = 1, jpl
333            DO jk = 1, nlay_i
334               DO jj = 1, jpj
335                  DO ji = 1, jpi
336                     !                          ! linear profile with 0 surface value
337                     zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
338                     zs  = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl)     ! weighting the profile
339                     sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
340                  END DO
341               END DO
342            END DO
343         END DO
344         !
345         DEALLOCATE( z_slope_s , zalpha )
346         !
347         !            !-------------------------------------------!
348      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
349         !            !-------------------------------------------!                                   (mean = 2.30)
350         !
351         s_i(:,:,:) = 2.30_wp
352!!gm Remark: if we keep the case 3, then compute an store one for all time-step
353!!           a array  S_prof(1:nlay_i) containing the calculation and just do:
354!         DO jk = 1, nlay_i
355!            sz_i(:,:,jk,:) = S_prof(jk)
356!         END DO
357!!gm end
358         !
359         DO jl = 1, jpl
360            DO jk = 1, nlay_i
361               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
362               sz_i(:,:,jk,jl) =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
363            END DO
364         END DO
365         !
366      END SELECT
367      !
368   END SUBROUTINE ice_var_salprof
369
370
371   SUBROUTINE ice_var_salprof1d
372      !!-------------------------------------------------------------------
373      !!                  ***  ROUTINE ice_var_salprof1d  ***
374      !!
375      !! ** Purpose :   1d computation of the sea ice salinity profile
376      !!                Works with 1d vectors and is used by thermodynamic modules
377      !!-------------------------------------------------------------------
378      INTEGER  ::   ji, jk    ! dummy loop indices
379      REAL(wp) ::   zargtemp, zsal, z1_dS   ! local scalars
380      REAL(wp) ::   zs, zs0              !   -      -
381      !
382      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z_slope_s, zalpha   !
383      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
384      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
385      !!-------------------------------------------------------------------
386      !
387      SELECT CASE ( nn_icesal )
388      !
389      !               !---------------------------------------!
390      CASE( 1 )       !  constant salinity in time and space  !
391         !            !---------------------------------------!
392         sz_i_1d(1:npti,:) = rn_icesal
393         !
394         !            !---------------------------------------------!
395      CASE( 2 )       !  time varying salinity with linear profile  !
396         !            !---------------------------------------------!
397         !
398         ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
399         !
400         !                                      ! Slope of the linear profile
401         WHERE( h_i_1d(1:npti) > epsi20 )   ;   z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
402         ELSEWHERE                          ;   z_slope_s(1:npti) = 0._wp
403         END WHERE
404         
405         z1_dS = 1._wp / ( zsi1 - zsi0 )
406         DO ji = 1, npti
407            zalpha(ji) = MAX(  0._wp , MIN(  ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp  )  )
408            !                             ! force a constant profile when SSS too low (Baltic Sea)
409            IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) )   zalpha(ji) = 0._wp
410         END DO
411         !
412         ! Computation of the profile
413         DO jk = 1, nlay_i
414            DO ji = 1, npti
415               !                          ! linear profile with 0 surface value
416               zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
417               zs  = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
418               sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
419            END DO
420         END DO
421         !
422         DEALLOCATE( z_slope_s, zalpha )
423
424         !            !-------------------------------------------!
425      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
426         !            !-------------------------------------------!                                   (mean = 2.30)
427         !
428         s_i_1d(1:npti) = 2.30_wp
429         !
430!!gm cf remark in ice_var_salprof routine, CASE( 3 )
431         DO jk = 1, nlay_i
432            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
433            zsal =  1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
434            DO ji = 1, npti
435               sz_i_1d(ji,jk) = zsal
436            END DO
437         END DO
438         !
439      END SELECT
440      !
441   END SUBROUTINE ice_var_salprof1d
442
443
444   SUBROUTINE ice_var_zapsmall
445      !!-------------------------------------------------------------------
446      !!                   ***  ROUTINE ice_var_zapsmall ***
447      !!
448      !! ** Purpose :   Remove too small sea ice areas and correct fluxes
449      !!-------------------------------------------------------------------
450      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
451      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch
452      !!-------------------------------------------------------------------
453      !
454      DO jl = 1, jpl       !==  loop over the categories  ==!
455         !
456         !-----------------------------------------------------------------
457         ! Zap ice energy and use ocean heat to melt ice
458         !-----------------------------------------------------------------
459         WHERE( a_i(:,:,jl) > epsi20 )   ;   h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
460         ELSEWHERE                       ;   h_i(:,:,jl) = 0._wp
461         END WHERE
462         !
463         WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 )   ;   zswitch(:,:) = 0._wp
464         ELSEWHERE                                                                           ;   zswitch(:,:) = 1._wp
465         END WHERE
466         !
467         DO jk = 1, nlay_i
468            DO jj = 1 , jpj
469               DO ji = 1 , jpi
470                  ! update exchanges with ocean
471                  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
472                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
473                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
474               END DO
475            END DO
476         END DO
477         !
478         DO jj = 1 , jpj
479            DO ji = 1 , jpi
480               ! update exchanges with ocean
481               sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoic * r1_rdtice
482               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoic * r1_rdtice
483               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhosn * r1_rdtice
484               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
485               !-----------------------------------------------------------------
486               ! Zap snow energy
487               !-----------------------------------------------------------------
488               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
489               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zswitch(ji,jj)
490               !
491               !-----------------------------------------------------------------
492               ! zap ice and snow volume, add water and salt to ocean
493               !-----------------------------------------------------------------
494               ato_i(ji,jj)    = a_i (ji,jj,jl) * ( 1._wp - zswitch(ji,jj) ) + ato_i(ji,jj)
495               a_i  (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
496               v_i  (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
497               v_s  (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
498               t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
499               oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
500               sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
501               !
502               h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
503               h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
504               !
505               a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
506               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
507               !
508            END DO
509         END DO
510         !
511      END DO 
512
513      ! to be sure that at_i is the sum of a_i(jl)
514      at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
515      vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
516
517      ! open water = 1 if at_i=0
518      WHERE( at_i(:,:) == 0._wp )   ato_i(:,:) = 1._wp
519      !
520   END SUBROUTINE ice_var_zapsmall
521
522
523   SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
524      !!-------------------------------------------------------------------
525      !!                ***  ROUTINE ice_var_itd   ***
526      !!
527      !! ** Purpose :  converting 1-cat ice to multiple ice categories
528      !!
529      !!                  ice thickness distribution follows a gaussian law
530      !!               around the concentration of the most likely ice thickness
531      !!                           (similar as iceistate.F90)
532      !!
533      !! ** Method:   Iterative procedure
534      !!               
535      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
536      !!
537      !!               2) Check whether the distribution conserves area and volume, positivity and
538      !!                  category boundaries
539      !!             
540      !!               3) If not (input ice is too thin), the last category is empty and
541      !!                  the number of categories is reduced (jpl-1)
542      !!
543      !!               4) Iterate until ok (SUM(itest(:) = 4)
544      !!
545      !! ** Arguments : zhti: 1-cat ice thickness
546      !!                zhts: 1-cat snow depth
547      !!                zati: 1-cat ice concentration
548      !!
549      !! ** Output    : jpl-cat
550      !!
551      !!  (Example of application: BDY forcings when input are cell averaged) 
552      !!-------------------------------------------------------------------
553      INTEGER  :: ji, jk, jl             ! dummy loop indices
554      INTEGER  :: idim, i_fill, jl0 
555      REAL(wp) :: zarg, zV, zconv, zdh, zdv
556      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
557      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i ! output ice/snow variables
558      INTEGER , DIMENSION(4)                  ::   itest
559      !!-------------------------------------------------------------------
560      !
561      ! ----------------------------------------
562      ! distribution over the jpl ice categories
563      ! ----------------------------------------
564      ! a gaussian distribution for ice concentration is used
565      ! then we check whether the distribution fullfills
566      ! volume and area conservation, positivity and ice categories bounds
567      idim = SIZE( zhti , 1 )
568      zh_i(1:idim,1:jpl) = 0._wp
569      zh_s(1:idim,1:jpl) = 0._wp
570      za_i(1:idim,1:jpl) = 0._wp
571      !
572      DO ji = 1, idim
573         !
574         IF( zhti(ji) > 0._wp ) THEN
575            !
576            ! find which category (jl0) the input ice thickness falls into
577            jl0 = jpl
578            DO jl = 1, jpl
579               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
580                  jl0 = jl
581                  CYCLE
582               ENDIF
583            END DO
584            !
585            itest(:) = 0
586            i_fill   = jpl + 1                                            !------------------------------------
587            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
588               !                                                          !------------------------------------
589               i_fill = i_fill - 1
590               !
591               zh_i(ji,1:jpl) = 0._wp
592               za_i(ji,1:jpl) = 0._wp
593               itest(:)       = 0     
594               !
595               IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
596                  zh_i(ji,1) = zhti(ji)
597                  za_i (ji,1) = zati (ji)
598               ELSE                         !-- case ice is thicker: fill categories >1
599                  ! thickness
600                  DO jl = 1, i_fill - 1
601                     zh_i(ji,jl) = hi_mean(jl)
602                  END DO
603                  !
604                  ! concentration
605                  za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
606                  DO jl = 1, i_fill - 1
607                     IF ( jl /= jl0 ) THEN
608                        zarg        = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
609                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
610                     ENDIF
611                  END DO
612                  !
613                  ! last category
614                  za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
615                  zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
616                  zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
617                  !
618                  ! correction if concentration of upper cat is greater than lower cat
619                  !    (it should be a gaussian around jl0 but sometimes it is not)
620                  IF ( jl0 /= jpl ) THEN
621                     DO jl = jpl, jl0+1, -1
622                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
623                           zdv = zh_i(ji,jl) * za_i(ji,jl)
624                           zh_i(ji,jl    ) = 0._wp
625                           za_i (ji,jl    ) = 0._wp
626                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
627                        END IF
628                     END DO
629                  ENDIF
630                  !
631               ENDIF
632               !
633               ! Compatibility tests
634               zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) ) 
635               IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation
636               !
637               zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )
638               IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation
639               !
640               IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ?
641               !
642               itest(4) = 1
643               DO jl = 1, i_fill
644                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0                                ! Test 4: positivity of ice concentrations
645               END DO
646               !                                         !----------------------------
647            END DO                                       ! end iteration on categories
648            !                                            !----------------------------
649         ENDIF
650      END DO
651
652      ! Add Snow in each category where za_i is not 0
653      DO jl = 1, jpl
654         DO ji = 1, idim
655            IF( za_i(ji,jl) > 0._wp ) THEN
656               zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
657               ! In case snow load is in excess that would lead to transformation from snow to ice
658               ! Then, transfer the snow excess into the ice (different from icethd_dh)
659               zdh = MAX( 0._wp, ( rhosn * zh_s(ji,jl) + ( rhoic - rau0 ) * zh_i(ji,jl) ) * r1_rau0 ) 
660               ! recompute h_i, h_s avoiding out of bounds values
661               zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
662               zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoic * r1_rhosn )
663            ENDIF
664         END DO
665      END DO
666      !
667   END SUBROUTINE ice_var_itd
668
669
670   SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )
671      !!-------------------------------------------------------------------
672      !!                ***  ROUTINE ice_var_itd2   ***
673      !!
674      !! ** Purpose :  converting N-cat ice to jpl ice categories
675      !!
676      !!                  ice thickness distribution follows a gaussian law
677      !!               around the concentration of the most likely ice thickness
678      !!                           (similar as iceistate.F90)
679      !!
680      !! ** Method:   Iterative procedure
681      !!               
682      !!               1) Fill ice cat that correspond to input thicknesses
683      !!                  Find the lowest(jlmin) and highest(jlmax) cat that are filled
684      !!
685      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1
686      !!                   by removing 25% ice area from jlmin and jlmax (resp.)
687      !!             
688      !!               3) Expand the filling to the empty cat between jlmin and jlmax
689      !!                   by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
690      !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
691      !!
692      !! ** Arguments : zhti: N-cat ice thickness
693      !!                zhts: N-cat snow depth
694      !!                zati: N-cat ice concentration
695      !!
696      !! ** Output    : jpl-cat
697      !!
698      !!  (Example of application: BDY forcings when inputs have N-cat /= jpl) 
699      !!-------------------------------------------------------------------
700      INTEGER  ::   ji, jl, jl1, jl2             ! dummy loop indices
701      INTEGER  ::   idim, icat 
702      INTEGER, PARAMETER ::   ztrans = 0.25_wp
703      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
704      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
705      INTEGER , DIMENSION(:,:), ALLOCATABLE   ::   jlfil, jlfil2
706      INTEGER , DIMENSION(:)  , ALLOCATABLE   ::   jlmax, jlmin
707      !!-------------------------------------------------------------------
708      !
709      idim = SIZE( zhti, 1 )
710      icat = SIZE( zhti, 2 )
711      !
712      ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays
713      ALLOCATE( jlmin(idim), jlmax(idim) )
714
715      ! --- initialize output fields to 0 --- !
716      zh_i(1:idim,1:jpl) = 0._wp
717      zh_s(1:idim,1:jpl) = 0._wp
718      za_i(1:idim,1:jpl) = 0._wp
719      !
720      ! --- fill the categories --- !
721      !     find where cat-input = cat-output and fill cat-output fields 
722      jlmax(:) = 0
723      jlmin(:) = 999
724      jlfil(:,:) = 0
725      DO jl1 = 1, jpl
726         DO jl2 = 1, icat
727            DO ji = 1, idim
728               IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
729                  ! fill the right category
730                  zh_i(ji,jl1) = zhti(ji,jl2)
731                  zh_s(ji,jl1) = zhts(ji,jl2)
732                  za_i(ji,jl1) = zati(ji,jl2)
733                  ! record categories that are filled
734                  jlmax(ji) = MAX( jlmax(ji), jl1 )
735                  jlmin(ji) = MIN( jlmin(ji), jl1 )
736                  jlfil(ji,jl1) = jl1
737               ENDIF
738            END DO
739         END DO
740      END DO
741      !
742      ! --- fill the gaps between categories --- ! 
743      !     transfer from categories filled at the previous step to the empty ones in between
744      DO ji = 1, idim
745         jl1 = jlmin(ji)
746         jl2 = jlmax(ji)
747         IF( jl1 > 1 ) THEN
748            ! fill the lower cat (jl1-1)
749            za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
750            zh_i(ji,jl1-1) = hi_mean(jl1-1)
751            ! remove from cat jl1
752            za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
753         ENDIF
754         IF( jl2 < jpl ) THEN
755            ! fill the upper cat (jl2+1)
756            za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
757            zh_i(ji,jl2+1) = hi_mean(jl2+1)
758            ! remove from cat jl2
759            za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
760         ENDIF
761      END DO
762      !
763      jlfil2(:,:) = jlfil(:,:) 
764      ! fill categories from low to high
765      DO jl = 2, jpl-1
766         DO ji = 1, idim
767            IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
768               ! fill high
769               za_i(ji,jl) = ztrans * za_i(ji,jl-1)
770               zh_i(ji,jl) = hi_mean(jl)
771               jlfil(ji,jl) = jl
772               ! remove low
773               za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
774            ENDIF
775         END DO
776      END DO
777      !
778      ! fill categories from high to low
779      DO jl = jpl-1, 2, -1
780         DO ji = 1, idim
781            IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
782               ! fill low
783               za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
784               zh_i(ji,jl) = hi_mean(jl) 
785               jlfil2(ji,jl) = jl
786               ! remove high
787               za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
788            ENDIF
789         END DO
790      END DO
791      !
792      DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
793      DEALLOCATE( jlmin, jlmax )
794      !
795   END SUBROUTINE ice_var_itd2
796
797
798   SUBROUTINE ice_var_bv
799      !!-------------------------------------------------------------------
800      !!                ***  ROUTINE ice_var_bv ***
801      !!
802      !! ** Purpose :   computes mean brine volume (%) in sea ice
803      !!
804      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
805      !!
806      !! References : Vancoppenolle et al., JGR, 2007
807      !!-------------------------------------------------------------------
808      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
809      !!-------------------------------------------------------------------
810      !
811!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
812!!   instead of setting everything to zero as just below
813      bv_i (:,:,:) = 0._wp
814      DO jl = 1, jpl
815         DO jk = 1, nlay_i
816            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
817               bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
818            END WHERE
819         END DO
820      END DO
821      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
822      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
823      END WHERE
824      !
825   END SUBROUTINE ice_var_bv
826
827
828   SUBROUTINE ice_var_enthalpy
829      !!-------------------------------------------------------------------
830      !!                   ***  ROUTINE ice_var_enthalpy ***
831      !!                 
832      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
833      !!
834      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
835      !!-------------------------------------------------------------------
836      INTEGER  ::   ji, jk   ! dummy loop indices
837      REAL(wp) ::   ztmelts  ! local scalar
838      !!-------------------------------------------------------------------
839      !
840      DO jk = 1, nlay_i             ! Sea ice energy of melting
841         DO ji = 1, npti
842            ztmelts      = - tmut  * sz_i_1d(ji,jk)
843            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point
844                                                                !   (sometimes zdf scheme produces abnormally high temperatures)   
845            e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
846               &                    + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
847               &                    - rcp  *   ztmelts )
848         END DO
849      END DO
850      DO jk = 1, nlay_s             ! Snow energy of melting
851         DO ji = 1, npti
852            e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
853         END DO
854      END DO
855      !
856   END SUBROUTINE ice_var_enthalpy
857
858#else
859   !!----------------------------------------------------------------------
860   !!   Default option         Dummy module           NO ESIM sea-ice model
861   !!----------------------------------------------------------------------
862#endif
863
864   !!======================================================================
865END MODULE icevar
Note: See TracBrowser for help on using the repository browser.