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_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

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

Last change on this file since 8564 was 8564, checked in by clem, 7 years ago

change variable names

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