source: branches/2016/dev_v3.20_2016_platelet/SOURCES/source_3.20/ice_th_dh.f

Last change on this file was 34, checked in by vancop, 8 years ago

bug in ice_th_dh.f

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