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 @ 9433

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

add two outputs: mean snow temperature and ice surface temperature per category

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