source: branches/2017/dev_v3.20_2017_transport_max/SOURCES/source_3.20/ice_th_dh.f @ 39

Last change on this file since 39 was 39, checked in by vancop, 7 years ago

add tank mass and salt balance

File size: 16.1 KB
Line 
1      SUBROUTINE ice_th_dh(nlay_s,nlay_i,kideb,kiut)
2       !!------------------------------------------------------------------
3       !!                ***         ROUTINE ice_th_dh         ***
4       !! ** Purpose :
5       !!           This routine determines variations of ice and snow thicknesses.
6       !! ** Method  :
7       !!           Ice/Snow surface melting arises from imbalance in surface fluxes
8       !!           Bottom accretion/ablation arises from flux budget
9       !!           Snow thickness can increase by precipitation and decrease by
10       !!              sublimation
11       !!           If snow load excesses Archmiede limit, snow-ice is formed by
12       !!              the flooding of sea-water in the snow
13       !! ** Steps 
14       !!           1) Compute available flux of heat for surface ablation
15       !!           2) Compute snow and sea ice enthalpies
16       !!           3) Surface ablation and sublimation
17       !!           4) Bottom accretion/ablation
18       !!           5) Case of Total ablation
19       !!           6) Snow ice formation
20       !!
21       !! ** Inputs / Outputs
22       !!
23       !! ** External
24       !!
25       !! ** References : Bitz and Lipscomb, JGR 99
26       !!                 Fichefet T. and M. Maqueda 1997,
27       !!                 J. Geophys. Res., 102(C6), 12609-12646   
28       !!                 Vancoppenolle, Fichefet and Bitz, GRL 2005
29       !!                 Vancoppenolle et al., OM08
30       !!
31       !! ** History  :
32       !!   original code    01-04 (LIM)
33       !!   original routine
34       !!               (05-2003) M. Vancoppenolle, Louvain-La-Neuve, Belgium
35       !!               (05-2008) BIO-LIM
36       !!
37       !!------------------------------------------------------------------
38       !! * Arguments
39       !!------------------------------------------------------------------
40
41      USE lib_fortran
42
43      INCLUDE 'type.com'
44      INCLUDE 'para.com'
45      INCLUDE 'const.com'
46      INCLUDE 'ice.com'
47      INCLUDE 'thermo.com'
48      INCLUDE 'tank.com'
49
50      ! Local Variables
51      DIMENSION zrchu1(nbpt), zrchu2(nbpt), zqsat(nbpt), z_f_surf(nbpt) 
52      DIMENSION zdeltah(maxnlay)
53      LOGICAL l_write
54
55      zqt_s_ini = 0.0
56      zqt_s_fin = 0.0
57      zdqt_s    = 0.0
58      zqt_i_ini = 0.0
59      zqt_i_fin = 0.0
60      zdqt_i    = 0.0
61      s_i_max   = 15.0
62
63      ! Local Constants
64      zeps = 1.0e-20
65      l_write = .TRUE.
66
67      IF ( l_write ) THEN
68         WRITE(numout,*) ' ** ice_th_dh : '
69         WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
70         WRITE(numout,*) 
71      ENDIF
72!
73!------------------------------------------------------------------------------|
74!  1) Calculate available heat for surface ablation         
75!------------------------------------------------------------------------------|
76!
77      DO 20 ji = kideb, kiut
78     
79      z_f_surf(ji)  = fratsb(ji) + fleb(ji) + fcsb(ji) - fc_su(ji)
80     &              + ab(ji)*fsolgb(ji)
81      z_f_surf(ji)  = MAX(c_zero,z_f_surf(ji))
82      z_f_surf(ji)  = z_f_surf(ji)*MAX(c_zero,
83     &                SIGN(one,t_su_b(ji)-tfs(ji)))
84
85      IF ( l_write ) THEN
86         WRITE(numout,*) ' Available heat for surface ablation ... '
87         WRITE(numout,*) 
88         WRITE(numout,*) ' z_f_surf : ', z_f_surf(ji)
89         WRITE(numout,*) ' fratsb   : ', fratsb(ji)
90         WRITE(numout,*) ' fleb     : ', fleb(ji)
91         WRITE(numout,*) ' fcsb     : ', fcsb(ji)
92         WRITE(numout,*) ' fc_su    : ', fc_su(ji)
93         WRITE(numout,*) ' ab       : ', ab(ji) 
94         WRITE(numout,*) ' fsolgb   : ', fsolgb(ji)
95         WRITE(numout,*)
96         WRITE(numout,*) ' ht_i_b    : ', ht_i_b(ji)
97         WRITE(numout,*) ' ht_s_b    : ', ht_s_b(ji)
98         WRITE(numout,*) ' t_su_b   : ', t_su_b(ji)
99         WRITE(numout,*)
100      ENDIF
101
102 20   CONTINUE
103
104!
105!------------------------------------------------------------------------------|
106!  2) Snowfall and surface melt                                                |
107!------------------------------------------------------------------------------|
108!
109      DO 40 ji = kideb, kiut
110
111      ! total snow heat content for conservation
112      zqt_s_ini = q_s_b(ji,1) * ht_s_b(ji)
113
114      IF ( l_write ) THEN
115         WRITE(numout,*) ' Surface ablation and sublimation ... '
116         WRITE(numout,*) 
117         WRITE(numout,*) ' zqt_s_ini : ', zqt_s_ini / ddtb
118         WRITE(numout,*) ' ht_s_b    : ', ht_s_b(ji)
119         WRITE(numout,*) ' q_s_b(1)  : ', q_s_b(ji,1)
120         WRITE(numout,*)
121      ENDIF
122
123      !----------
124      ! Snowfall
125      !----------
126      dh_s_prec(ji)    =  hnpbqb(ji)
127      dh_s_melt(ji)    =  0.0
128      zqprec           =  rhon * ( cpg * ( tpw - tabqb(ji) ) + lfus )
129
130      ! Conservation update
131      zqt_s_ini        =  zqt_s_ini + zqprec*hnpbqb(ji)
132      fprec            =  - zqprec * hnpbqb(ji) / ddtb
133
134      IF ( l_write ) THEN
135         WRITE(numout,*) ' snow falls! '
136         WRITE(numout,*) ' dh_s_prec : ', dh_s_prec(ji)
137         WRITE(numout,*) ' flux of h : ',
138     &                   zqprec*hnpbqb(ji) / ddtb
139         WRITE(numout,*) ' zqt_s_ini : ', zqt_s_ini / ddtb
140         WRITE(numout,*)
141      ENDIF
142
143      !-----------
144      ! Snow melt
145      !-----------
146      ! Energy available for surface melt
147      zqfont_su =  ( z_f_surf(ji) + f_s_im(ji) ) * ddtb 
148      IF ( l_write ) THEN
149         WRITE(numout,*) ' snow melts! '
150         WRITE(numout,*) ' z_f_surf(ji) : ', z_f_surf(ji)
151         WRITE(numout,*) ' f_s_im       : ', f_s_im(ji)
152         WRITE(numout,*) ' zqfont_su    : ', zqfont_su
153      ENDIF
154
155      ! Melt of fallen snow
156      zdeltah(1)       =  MIN( 0.0 , - zqfont_su / 
157     &                    MAX( zqprec , zeps ) )
158      zqfont_su        =  MAX( 0.0 , - dh_s_prec(ji) - zdeltah(1) ) * 
159     &                    zqprec
160      zdeltah(1)       =  MAX( - dh_s_prec(ji), zdeltah(1) )
161      dh_s_melt(ji)    =  dh_s_melt(ji) + zdeltah(1)
162
163      ! Melt / evaporation of snow
164      zswi_evap = 0
165      WRITE(numout,*) ' ln_evap : ', ln_evap
166      WRITE(numout,*) ' tdewb   : ', tdewb(ji)
167      IF ( ln_evap .AND. ( qsfcb(ji) .GT. qabqb(ji) ) .AND.
168     &     tdewb(ji) .LT. 1. ) zswi_evap = 1
169     
170      DO layer = 1, nlay_s ! in case of melting of more than 1 layer
171         zlvap = zswi_evap * rhon * ( lvap + cpg * ( tabqb(ji) - tpw ) ) ! MV new evaporation
172
173!        zdeltah(layer) =  - zqfont_su / q_s_b(ji,layer)
174         zdeltah(layer) =  - zqfont_su / ( q_s_b(ji,layer) + zlvap )
175         zqfont_su      = MAX( c_zero, - deltaz_s_phy(layer) - 
176     &                    zdeltah(layer) ) * 
177!    &                    q_s_b(ji,layer) !MV new evaporation
178     &                    ( q_s_b(ji,layer) + zlvap )
179         zdeltah(layer) = MAX( zdeltah(layer), - deltaz_s_phy(layer) )
180         dh_s_melt(ji) =  dh_s_melt(ji) + zdeltah(layer) !resulting melt of snow   
181      END DO
182
183      dh_s_tot(ji)     =  dh_s_melt(ji) + dh_s_prec(ji)
184
185      ! old and new snow thicknesses
186      hsold            =  ht_s_b(ji)
187      hsnew            =  ht_s_b(ji) + dh_s_tot(ji)
188
189      ! if snow is still present zhn = 1, else zhn = 0
190      zhn              =  1.0 - MAX( c_zero , SIGN( one , - hsnew ) )
191      ht_s_b(ji)       =  MAX( c_zero , hsnew )
192
193      !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
194      ! Conservation test for snow
195      zqt_s_fin = q_s_b(ji,1) * ht_s_b(ji)
196      zdqt_s = zqt_s_fin - zqt_s_ini
197
198      WRITE(numout,*)
199      WRITE(numout,*) ' Conservation in snow... '
200      WRITE(numout,*) ' dh_s_melt : ', dh_s_melt(ji)
201      WRITE(numout,*) ' dh_s_prec : ', dh_s_prec(ji)
202      WRITE(numout,*) ' ht_s_b    : ', ht_s_b(ji)
203      WRITE(numout,*) ' zqt_s_ini : ', zqt_s_ini / ddtb
204      WRITE(numout,*) ' zqt_s_fin : ', zqt_s_fin / ddtb
205      WRITE(numout,*) ' zdqt_s    : ', zdqt_s / ddtb
206      WRITE(numout,*) ' z_f_surf  : ', - z_f_surf(ji)
207      IF ( zqt_s_fin.GT.0.0 ) THEN
208         cons_err = ABS(zdqt_s / ddtb  + z_f_surf(ji) )
209      ELSE
210         cons_err = ABS(zqt_s_ini / ddtb + zdqt_s / ddtb )
211      ENDIF
212      WRITE(numout,*) ' Cons error, snow : ', cons_err
213      WRITE(numout,*)
214      !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
215
216      !------------------
217      ! Ice surface melt
218      !------------------
219      IF ( l_write ) WRITE(numout,*) ' ice melts!  '
220
221      zzf_surf = zqfont_su / ddtb
222      zdqt_i            = 0.0
223      dh_i_surf(ji)     =  0.0
224      DO layer = 1, nlay_i
225         zdeltah(layer) =  - zqfont_su / q_i_b(ji,layer)
226         zqfont_su      =  MAX( c_zero , - deltaz_i_phy(layer) 
227     &                  - zdeltah(layer) ) *  q_i_b(ji,layer)
228         zdeltah(layer) =  MAX( zdeltah(layer) , - deltaz_i_phy(layer) )
229         dh_i_surf(ji)  =  dh_i_surf(ji) + zdeltah(layer) !resulting melt of ice
230         zdqt_i         = zdqt_i  + zdeltah(layer) * q_i_b(ji,layer) 
231     &                  / ddtb
232      END DO
233
234      cons_err = ABS( zzf_surf + zdqt_i )
235
236      IF ( l_write ) THEN
237         WRITE(numout,*) ' Conservation in sea ice, surface '
238         WRITE(numout,*) ' dh_i_surf: ', dh_i_surf(ji)
239         WRITE(numout,*) ' ht_i_b   : ', ht_i_b(ji)
240         WRITE(numout,*) ' zzf_surf : ', zzf_surf
241         WRITE(numout,*) ' zdqt_i   : ', zdqt_i
242         WRITE(numout,*) ' Cons error, ice : ', cons_err
243         WRITE(numout,*)
244      ENDIF
245     
246!
247!------------------------------------------------------------------------------|
248!  3) Sublimation at the surface                                               |
249!------------------------------------------------------------------------------|
250!
251      !------------------
252      ! Snow sublimation
253      !------------------
254
255      !------------------
256      ! Ice sublimation
257      !------------------
258
259      !-------------------
260      ! Snow condensation
261      !-------------------
262
263      ! THIS IS DEFINITELY WRONG SINCE LATENT HEAT IS ALREADY COUNTED TO
264      ! remove heat from the snow pack
265!     ! 4.3) Snow/ice sublimation
266!     !
267!     ! If fleb is negative, snow condensates at the surface.
268!     !
269      dh_s_subl(ji)    =  + parsub*fleb(ji)/(rhon*lsub)*ddtb
270      dh_s_tot(ji)     =  dh_s_tot(ji) + dh_s_subl(ji)
271      zdhcf            =  ht_s_b(ji) + dh_s_subl(ji) 
272      ht_s_b(ji)       =  MAX(c_zero,zdhcf)
273      dh_s_tot(ji)     =  ht_s_b(ji) - hsold
274      dh_i_subl(ji)    =  - MAX(c_zero,-zdhcf)*rhon/rhog
275
276      dh_i_surf(ji)    =  dh_i_surf(ji) + dh_i_subl(ji)
277
278      hsnew            =  ht_s_b(ji)
279
280      IF ( ht_s_b(ji) .LE. 0.0 ) THEN
281         dh_s_tot(ji) =  MAX( 0.0 , dh_s_tot(ji) )
282      ENDIF
283
284      IF ( l_write ) THEN
285         WRITE(numout,*) ' Snow sublimation ... '
286         WRITE(numout,*) ' '
287         WRITE(numout,*) ' parsub : ', parsub
288         WRITE(numout,*) ' dh_s_subl : ', dh_s_subl(ji)
289         WRITE(numout,*) ' dh_i_subl : ', dh_i_subl(ji)
290      ENDIF
291
292 40   CONTINUE
293!
294!------------------------------------------------------------------------------|
295!  4) Basal growth and melt                                                    |
296!------------------------------------------------------------------------------|
297!
298      DO 50 ji = kideb, kiut
299
300      IF ( l_write ) THEN
301         WRITE(numout,*) ' Basal growth and melt ... '
302         WRITE(numout,*) 
303         WRITE(numout,*) ' fbbqb     : ', fbbqb(ji)
304         WRITE(numout,*) ' fc_bo_i   : ', fc_bo_i(ji)
305         WRITE(numout,*)
306      ENDIF
307
308      ! formation / melting of ice at the base is determined by the balance of
309      ! the conductive heat flux in the ice (fc_bo_i), and the heat fluxes
310      ! from the ocean (fbbqb). Conductive heat flux is positive
311      ! downwards and fbbq is positive to the ice, i.e., upwards.
312      ! Melt/formation rate is modulated by the enthalpy of the bottom ice layer.
313      ! accretion of ice at the base
314
315      !--------------
316      ! Basal growth
317      !--------------
318      IF ( ( fc_bo_i(ji) + fbbqb(ji) ) .LT. 0.0 ) THEN
319 
320         s_i_new = rn_e_newice * s_w ! New ice salinity (g/kg)
321 
322         ztmelts  = - tmut * s_i_new + tpw ! Melting point in K
323 
324         ! New ice heat content (Bitz and Lipscomb, 1999)
325         ! should be zdE instead
326         q_i_b(ji,nlay_i+1) = rhog*( cpg*(tmelts-t_bo_b(ji))
327     &                 + lfus*( 1.0-(tmelts-tpw)
328     &                 / (t_bo_b(ji) - tpw) )
329     &                 - cpw*(tmelts-tpw) ) 
330
331         zE1 = - cpw * ( t_bo_b(ji) - tpw )          ! specific enthalpy of sea water  <0
332         zE2 = - q_i_b(ji,nlay_i+1) / rhog ! specific enthalpy of new sea ice  <0
333
334         WRITE(numout,*) ' zE1 : ', zE1
335         WRITE(numout,*) ' zE2 : ', zE2
336
337         zdE = zE2 - zE1 ! <0
338
339         WRITE(numout,*) ' zdE : ', zdE
340
341         dh_i_bott(ji)    =  ddtb*(fc_bo_i(ji) + fbbqb(ji) )
342     &                    / ( zdE * rhog )
343 
344         ! salt flux due to initial salt entrapment (keep ?)
345         fsbp = s_w * ( 1. - rn_e_newice ) * dh_i_bott(ji) / ddtb *
346     &          rhog / 1000.
347 
348         IF ( l_write ) WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
349 
350      ENDIF
351
352      !-----------------
353      ! Basal melt
354      !-----------------
355      IF ( ( fc_bo_i(ji) + fbbqb(ji) ) .GE. 0.0 ) THEN
356
357         IF ( l_write ) WRITE(numout,*) ' Energy available for
358     &                                    basal melt   : ',
359     &                  fc_bo_i(ji) + fbbqb(ji)
360
361         zqfont_bo   = ddtb * ( fc_bo_i(ji) + fbbqb(ji) )
362         zzf_base    = zqfont_bo / ddtb
363         zdqt_i      = 0.0
364
365         IF ( l_write ) WRITE(numout,*) ' zqfont_bo : ', zqfont_bo
366
367         dh_i_bott(ji)     =  0.0
368         DO layer = nlay_i, 1, -1
369            zdeltah(layer) =  - zqfont_bo / q_i_b(ji,layer)
370            zqfont_bo      =  MAX ( 0.0 , - deltaz_i_phy(layer) -
371     &                        zdeltah(layer) )
372     &                     *  q_i_b(ji,layer)
373            dh_i_bott(ji)  =  dh_i_bott(ji) + zdeltah(layer)
374            zdqt_i         =  zdqt_i + zdeltah(layer) * 
375     &                        q_i_b(ji,layer) / ddtb
376         END DO
377
378         IF ( l_write ) WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
379
380         cons_err = ABS( zzf_base + zdqt_i )
381
382         IF ( l_write ) THEN
383            WRITE(numout,*) ' Conservation in sea ice, base '
384            WRITE(numout,*) ' dh_i_bott: ', dh_i_bott(ji)
385            WRITE(numout,*) ' ht_i_b   : ', ht_i_b(ji)
386            WRITE(numout,*) ' zzf_base : ', zzf_base
387            WRITE(numout,*) ' zdqt_i   : ', zdqt_i
388!           WRITE(numout,*) ' Conservation error  ice surface : ',
389!    &                      cons_err
390!           WRITE(numout,*)
391         ENDIF
392
393      ENDIF
394
395      ! It can be than an internal temperature is greater than melt point
396      ! then, see lim3 for correction
397
398      ! new ice thickness
399      zhgnew         = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji)
400      old_ht_i_b(ji) = ht_i_b(ji)
401
402      ht_i_b(ji) = zhgnew
403
404 50   CONTINUE
405
406!
407!------------------------------------------------------------------------------|
408!  5) Formation of snow-ice                                                    |
409!------------------------------------------------------------------------------|
410!
411      ! When snow load excesses Archimede's limit, snow-ice interface goes down
412      ! under sea-level, flooding of seawater transforms snow into ice
413      ! dh_snowice is positive for the ice
414
415      DO 70 ji = kideb, kiut
416
417      dh_snowice(ji) = MAX( c_zero , ( rhon * ht_s_b(ji) + (rhog -rho0 )
418     &                 * ht_i_b(ji)) / ( rhon + rho0 - rhog ) )
419       
420      zhgnew         = MAX( zhgnew , zhgnew + dh_snowice(ji) )
421      zhnnew         = MIN( ht_s_b(ji) , ht_s_b(ji) - dh_snowice(ji) )
422
423      ht_s_b(ji)  = zhnnew
424      ht_i_b(ji)  = zhgnew
425
426      IF ( l_write ) THEN
427         WRITE(numout,*) ' dh_snowice : ', dh_snowice(ji)
428         WRITE(numout,*)
429         WRITE(numout,*) ' At the end of the routine ... '
430         WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji)
431         WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji)
432      ENDIF
433
434      ! give a minimum snow depth
435      zhsnmin        = 1.0e-10
436      isnow = 0
437      IF ( ht_s_b(ji) .GT. zhsnmin ) isnow = 1
438      ht_s_b(ji) = isnow * ht_s_b(ji)
439
440      !--- Remove ---
441      WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
442      WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji)
443      WRITE(numout,*) ' dh_i_snowice : ', dh_snowice(ji)
444      WRITE(numout,*) ' dh_s_melt : ', dh_s_melt(ji)
445      !--- Remove ---
446
447 70   CONTINUE
448
449      RETURN
450
451!------------------------------------------------------------------------------|
452! Fin de la subroutine ice_th_dh
453      END SUBROUTINE
Note: See TracBrowser for help on using the repository browser.