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.
limthd_dh.F90 in trunk/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMO/LIM_SRC_3/limthd_dh.F90 @ 1571

Last change on this file since 1571 was 1571, checked in by ctlod, 15 years ago

Correct 3 bugs in LIM-3, see ticket: #505
  • Property svn:keywords set to Id
File size: 35.0 KB
Line 
1MODULE limthd_dh
2#if defined key_lim3
3   !!----------------------------------------------------------------------
4   !!   'key_lim3'                                      LIM3 sea-ice model
5   !!----------------------------------------------------------------------
6   !!======================================================================
7   !!                       ***  MODULE limthd_dh ***
8   !!                thermodynamic growth and decay of the ice
9   !!======================================================================
10
11   !!----------------------------------------------------------------------
12   !!   lim_thd_dh  : vertical accr./abl. and lateral ablation of sea ice
13   !! * Modules used
14
15   USE par_oce          ! ocean parameters
16   USE phycst           ! physical constants (OCE directory)
17   USE sbc_oce          ! Surface boundary condition: ocean fields
18   USE thd_ice
19   USE iceini
20   USE limistate
21   USE in_out_manager
22   USE ice
23   USE par_ice
24   USE lib_mpp
25
26   IMPLICIT NONE
27   PRIVATE
28
29   !! * Routine accessibility
30   PUBLIC lim_thd_dh        ! called by lim_thd
31
32   !! * Module variables
33   REAL(wp)  ::           &  ! constant values
34      epsi20 = 1e-20   ,  &
35      epsi13 = 1e-13   ,  &
36      epsi16 = 1e-16   ,  &
37      zzero  = 0.e0    ,  &
38      zone   = 1.e0
39
40   !!----------------------------------------------------------------------
41   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2005)
42   !! $Id$
43   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45
46CONTAINS
47
48   SUBROUTINE lim_thd_dh(kideb,kiut,jl)
49      !!------------------------------------------------------------------
50      !!                ***  ROUTINE lim_thd_dh  ***
51      !!------------------------------------------------------------------
52      !! ** Purpose :
53      !!           This routine determines variations of ice and snow thicknesses.
54      !! ** Method  :
55      !!           Ice/Snow surface melting arises from imbalance in surface fluxes
56      !!           Bottom accretion/ablation arises from flux budget
57      !!           Snow thickness can increase by precipitation and decrease by
58      !!              sublimation
59      !!           If snow load excesses Archmiede limit, snow-ice is formed by
60      !!              the flooding of sea-water in the snow
61      !! ** Steps 
62      !!           1) Compute available flux of heat for surface ablation
63      !!           2) Compute snow and sea ice enthalpies
64      !!           3) Surface ablation and sublimation
65      !!           4) Bottom accretion/ablation
66      !!           5) Case of Total ablation
67      !!           6) Snow ice formation
68      !!
69      !! ** Arguments
70      !!
71      !! ** Inputs / Outputs
72      !!
73      !! ** External
74      !!
75      !! ** References : Bitz and Lipscomb, JGR 99
76      !!                 Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
77      !!                 Vancoppenolle, Fichefet and Bitz, GRL 2005
78      !!                 Vancoppenolle et al., OM08
79      !!
80      !! ** History  :
81      !!   original code    01-04 (LIM)
82      !!   original routine
83      !!               (05-2003) M. Vancoppenolle, Louvain-La-Neuve, Belgium
84      !!               (06/07-2005) 3D version
85      !!               (03-2008)    Clean code
86      !!
87      !!------------------------------------------------------------------
88      !! * Arguments
89      INTEGER , INTENT (IN) ::  &
90         kideb     ,         &  !: Start point on which the  the computation is applied
91         kiut      ,         &  !: End point on which the  the computation is applied
92         jl                     !: Thickness cateogry number
93
94      !! * Local variables
95      INTEGER ::             & 
96         ji        ,         &  !: space index
97         jk        ,         &  !: ice layer index
98         isnow     ,         &  !: switch for presence (1) or absence (0) of snow
99         zji       ,         &  !: 2D corresponding indices to ji
100         zjj       ,         &
101         isnowic   ,         &  !: snow ice formation not
102         i_ice_switch   ,    &  !: ice thickness above a certain treshold or not
103         iter
104
105      REAL(wp) ::            &
106         zzfmass_i ,         &
107         zzfmass_s ,         &
108         zhsnew    ,         &  !: new snow thickness
109         zihgnew   ,         &  !: switch for total ablation
110         ztmelts   ,         &  !: melting point
111         zhn       ,         &
112         zdhcf     ,         &
113         zdhbf     ,         &
114         zhni      ,         &
115         zhnfi     ,         &
116         zihg      ,         &
117         zdhnm     ,         &
118         zhnnew    ,         &
119         zeps = 1.0e-13,     &
120         zhisn     ,         &
121         zfracs    ,         &  !: fractionation coefficient for bottom salt
122                                !: entrapment
123         zds       ,         &  !: increment of bottom ice salinity
124         zcoeff    ,         &  !: dummy argument for snowfall partitioning
125                                !: over ice and leads
126         zsm_snowice,        &  !: snow-ice salinity
127         zswi1     ,         &  !: switch for computation of bottom salinity
128         zswi12    ,         &  !: switch for computation of bottom salinity
129         zswi2     ,         &  !: switch for computation of bottom salinity
130         zgrr      ,         &  !: bottom growth rate
131         zihic     ,         &  !:
132         ztform                 !: bottom formation temperature
133
134      REAL(wp) , DIMENSION(jpij) ::  &
135         zh_i      ,         &  ! ice layer thickness
136         zh_s      ,         &  ! snow layer thickness
137         ztfs      ,         &  ! melting point
138         zhsold    ,         &  ! old snow thickness
139         zqprec    ,         &  !: energy of fallen snow
140         zqfont_su ,         &  ! incoming, remaining surface energy
141         zqfont_bo              ! incoming, bottom energy
142
143      REAL(wp) , DIMENSION(jpij) :: &
144         z_f_surf,           &  ! surface heat for ablation
145         zhgnew  ,           &  ! new ice thickness
146         zfmass_i
147
148      REAL(wp), DIMENSION(jpij) :: &
149         zdh_s_mel  ,        &  ! snow melt
150         zdh_s_pre  ,        &  ! snow precipitation
151         zdh_s_sub  ,        &  ! snow sublimation
152         zfsalt_melt            ! salt flux due to ice melt
153
154      REAL(wp) , DIMENSION(jpij,jkmax) :: &
155         zdeltah
156
157      ! Pathological cases
158      REAL(wp), DIMENSION(jpij) :: &
159         zfdt_init  ,        &  !: total incoming heat for ice melt
160         zfdt_final ,        &  !: total remaing heat for ice melt
161         zqt_i      ,        &  !: total ice heat content
162         zqt_s      ,        &  !: total snow heat content
163         zqt_dummy              !: dummy heat content
164
165      REAL(wp), DIMENSION(jpij,jkmax) :: &
166         zqt_i_lay              !: total ice heat content
167
168      ! Heat conservation
169      REAL(wp), DIMENSION(jpij) :: &
170         zfbase,             &
171         zdq_i
172
173      INTEGER, DIMENSION(jpij) ::  &
174         innermelt
175
176      REAL(wp) :: &
177         meance_dh
178
179      INTEGER ::                   &
180         num_iter_max,       &
181         numce_dh
182
183      zfsalt_melt(:)  = 0.0
184      ftotal_fin(:)   = 0.0
185      zfdt_init(:)    = 0.0
186      zfdt_final(:)   = 0.0
187
188      DO ji = kideb, kiut
189         old_ht_i_b(ji) = ht_i_b(ji)
190         old_ht_s_b(ji) = ht_s_b(ji)
191      END DO
192      !
193      !------------------------------------------------------------------------------!
194      !  1) Calculate available heat for surface ablation                            !
195      !------------------------------------------------------------------------------!
196      !
197      DO ji = kideb,kiut
198         isnow         = INT( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ) )
199         ztfs(ji)      = isnow * rtt + ( 1.0 - isnow ) * rtt
200         z_f_surf(ji)  = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) *                 & 
201            qsr_ice_1d(ji) - fc_su(ji)
202         z_f_surf(ji)  = MAX( zzero , z_f_surf(ji) ) *                        &
203            MAX( zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) )
204         zfdt_init(ji) = ( z_f_surf(ji) + &
205            MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) )  &
206            * rdt_ice
207      END DO ! ji
208
209      zqfont_su(:) = 0.0
210      zqfont_bo(:) = 0.0
211      dsm_i_se_1d(:) = 0.0     
212      dsm_i_si_1d(:) = 0.0     
213      !
214      !------------------------------------------------------------------------------!
215      !  2) Computing layer thicknesses and  snow and sea-ice enthalpies.            !
216      !------------------------------------------------------------------------------!
217      !
218      ! Layer thickness
219      DO ji = kideb,kiut
220         zh_i(ji) = ht_i_b(ji) / nlay_i
221         zh_s(ji) = ht_s_b(ji) / nlay_s
222      END DO
223
224      ! Total enthalpy of the snow
225      zqt_s(:) = 0.0
226      DO jk = 1, nlay_s
227         DO ji = kideb,kiut
228            zqt_s(ji)    =  zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s
229         END DO
230      END DO
231
232      ! Total enthalpy of the ice
233      zqt_i(:) = 0.0
234      DO jk = 1, nlay_i
235         DO ji = kideb,kiut
236            zqt_i(ji)        =  zqt_i(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
237            zqt_i_lay(ji,jk) =  q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
238         END DO
239      END DO
240      !
241      !------------------------------------------------------------------------------|
242      !  3) Surface ablation and sublimation                                         |
243      !------------------------------------------------------------------------------|
244      !
245      !-------------------------
246      ! 3.1 Snow precips / melt
247      !-------------------------
248      ! Snow accumulation in one thermodynamic time step
249      ! snowfall is partitionned between leads and ice
250      ! if snow fall was uniform, a fraction (1-at_i) would fall into leads
251      ! but because of the winds, more snow falls on leads than on sea ice
252      ! and a greater fraction (1-at_i)^beta of the total mass of snow
253      ! (beta < 1) falls in leads.
254      ! In reality, beta depends on wind speed,
255      ! and should decrease with increasing wind speed but here, it is
256      ! considered as a constant. an average value is 0.66
257      ! Martin Vancoppenolle, December 2006
258
259      ! Snow fall
260      DO ji = kideb, kiut
261         zcoeff = ( 1.0 - ( 1.0 - at_i_b(ji) )**betas ) / at_i_b(ji) 
262         zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn
263      END DO
264      zdh_s_mel(:) =  0.0
265
266      ! Melt of fallen snow
267      DO ji = kideb, kiut
268         ! tatm_ice is now in K
269         zqprec(ji)     =  rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus ) 
270         zqfont_su(ji)  =  z_f_surf(ji) * rdt_ice
271         zdeltah(ji,1)  =  MIN( 0.0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) )
272         zqfont_su(ji)  =  MAX( 0.0 , - zdh_s_pre(ji) - zdeltah(ji,1) )      * &
273            zqprec(ji)
274         zdeltah(ji,1)  =  MAX( - zdh_s_pre(ji) , zdeltah(ji,1) )
275         zdh_s_mel(ji)  =  zdh_s_mel(ji) + zdeltah(ji,1)
276         ! heat conservation
277         qt_s_in(ji,jl) =  qt_s_in(ji,jl) + zqprec(ji) * zdh_s_pre(ji)
278         zqt_s(ji)      =  zqt_s(ji) + zqprec(ji) * zdh_s_pre(ji)
279         zqt_s(ji)      =  MAX ( zqt_s(ji) - zqfont_su(ji) , 0.0 ) 
280      END DO
281
282
283      ! Snow melt due to surface heat imbalance
284      DO jk = 1, nlay_s
285         DO ji = kideb, kiut
286            zdeltah(ji,jk) = - zqfont_su(ji) / q_s_b(ji,jk)
287            zqfont_su(ji)  =  MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * &
288               q_s_b(ji,jk) 
289            zdeltah(ji,jk) =  MAX( zdeltah(ji,jk) , - zh_s(ji) )
290            zdh_s_mel(ji)  =  zdh_s_mel(ji) + zdeltah(ji,jk) !resulting melt of snow   
291         END DO
292      END DO
293
294      ! Apply snow melt to snow depth
295      DO ji = kideb, kiut
296         dh_s_tot(ji)   =  zdh_s_mel(ji) + zdh_s_pre(ji)
297         ! Old and new snow depths
298         zhsold(ji)     =  ht_s_b(ji)
299         zhsnew         =  ht_s_b(ji) + dh_s_tot(ji)
300         ! If snow is still present zhn = 1, else zhn = 0
301         zhn            =  1.0 - MAX( zzero , SIGN( zone , - zhsnew ) )
302         ht_s_b(ji)     =  MAX( zzero , zhsnew )
303         ! Volume and mass variations of snow
304         dvsbq_1d(ji)   =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji)              &
305            - zdh_s_mel(ji) )
306         dvsbq_1d(ji)   =  MIN( zzero, dvsbq_1d(ji) )
307         rdmsnif_1d(ji) =  rdmsnif_1d(ji) + rhosn * dvsbq_1d(ji)
308      END DO ! ji
309
310      !--------------------------
311      ! 3.2 Surface ice ablation
312      !--------------------------
313      DO ji = kideb, kiut 
314         dh_i_surf(ji) =  0.0
315         ! For heat conservation test
316         z_f_surf(ji)  =  zqfont_su(ji) / rdt_ice ! heat conservation test
317         zdq_i(ji)     =  0.0
318      END DO ! ji
319
320      DO jk = 1, nlay_i
321         DO ji = kideb, kiut 
322            ! melt of layer jk
323            zdeltah(ji,jk)      = - zqfont_su(ji) / q_i_b(ji,jk)
324            ! recompute heat available
325            zqfont_su(ji)       = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) *  &
326               q_i_b(ji,jk) 
327            ! melt of layer jk cannot be higher than its thickness
328            zdeltah(ji,jk)      = MAX( zdeltah(ji,jk) , - zh_i(ji) )
329            ! update surface melt
330            dh_i_surf(ji)       = dh_i_surf(ji) + zdeltah(ji,jk) 
331            ! for energy conservation
332            zdq_i(ji)           = zdq_i(ji) + zdeltah(ji,jk) *                &
333               q_i_b(ji,jk) / rdt_ice
334            ! contribution to ice-ocean salt flux
335            zji                 = MOD( npb(ji) - 1, jpi ) + 1
336            zjj                 = ( npb(ji) - 1 ) / jpi + 1
337            zfsalt_melt(ji)     = zfsalt_melt(ji) +                           &
338               ( sss_m(zji,zjj) - sm_i_b(ji)   ) *         &
339               a_i_b(ji) *                                 &
340               MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 
341         END DO ! ji
342      END DO ! jk
343
344      !-------------------
345      ! Conservation test
346      !-------------------
347      IF ( con_i ) THEN
348         numce_dh = 0
349         meance_dh = 0.0
350         DO ji = kideb, kiut
351
352            IF ( ( z_f_surf(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN
353               numce_dh  = numce_dh + 1
354               meance_dh = meance_dh + z_f_surf(ji) + zdq_i(ji)
355            ENDIF
356
357            IF ( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN!
358               WRITE(numout,*) ' ALERTE heat loss for surface melt '
359               WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl
360               WRITE(numout,*) ' ht_i_b  : ', ht_i_b(ji)
361               WRITE(numout,*) ' z_f_surf  : ', z_f_surf(ji)
362               WRITE(numout,*) ' zdq_i   : ', zdq_i(ji)
363               WRITE(numout,*) ' ht_i_b  : ', ht_i_b(ji)
364               WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji)
365               WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji)
366               WRITE(numout,*) ' qlbbq_1d: ', qlbbq_1d(ji)
367               WRITE(numout,*) ' s_i_new : ', s_i_new(ji)
368               WRITE(numout,*) ' sss_m   : ', sss_m(zji,zjj)
369            ENDIF
370
371         END DO ! ji
372
373         IF ( numce_dh .GT. 0 ) meance_dh = meance_dh / numce_dh
374         WRITE(numout,*) ' Error report - Category : ', jl
375         WRITE(numout,*) ' ~~~~~~~~~~~~ '
376         WRITE(numout,*) ' Number of points where there is sur. me. error : ', numce_dh
377         WRITE(numout,*) ' Mean basal growth error on error points : ', meance_dh
378
379      ENDIF ! con_i
380
381      !----------------------
382      ! 3.3 Snow sublimation
383      !----------------------
384
385      DO ji = kideb, kiut
386         ! if qla is positive (upwards), heat goes to the atmosphere, therefore
387         ! snow sublimates, if qla is negative (downwards), snow condensates
388         zdh_s_sub(ji)   =  - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice
389         dh_s_tot(ji)    =  dh_s_tot(ji) + zdh_s_sub(ji)
390         zdhcf           =  ht_s_b(ji) + zdh_s_sub(ji) 
391         ht_s_b(ji)      =  MAX( zzero , zdhcf )
392         ! we recompute dh_s_tot
393         dh_s_tot(ji)    =  ht_s_b(ji) - zhsold(ji)
394         qt_s_in(ji,jl) =  qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1)
395      END DO  !ji
396
397      zqt_dummy(:) = 0.0
398      DO jk = 1, nlay_s
399         DO ji = kideb,kiut
400            q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
401            ! heat conservation
402            zqt_dummy(ji)    =  zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s
403         END DO
404      END DO
405
406      DO jk = 1, nlay_s !n
407         DO ji = kideb, kiut !n
408            ! In case of disparition of the snow, we have to update the snow
409            ! temperatures
410            zhisn  =  MAX( zzero , SIGN( zone, - ht_s_b(ji) ) )
411            t_s_b(ji,jk) = ( 1.0 - zhisn ) * t_s_b(ji,jk) + zhisn * rtt
412            q_s_b(ji,jk) = ( 1.0 - zhisn ) * q_s_b(ji,jk)
413         END DO
414      END DO
415
416      !
417      !------------------------------------------------------------------------------!
418      ! 4) Basal growth / melt                                                       !
419      !------------------------------------------------------------------------------!
420      !
421      ! Ice basal growth / melt is given by the ratio of heat budget over basal
422      ! ice heat content.  Basal heat budget is given by the difference between
423      ! the inner conductive flux  (fc_bo_i), from the open water heat flux
424      ! (qlbbqb) and the turbulent ocean flux (fbif).
425      ! fc_bo_i is positive downwards. fbif and qlbbq are positive to the ice
426
427      !-----------------------------------------------------
428      ! 4.1 Basal growth - (a) salinity not varying in time
429      !-----------------------------------------------------
430      IF ( ( num_sal .NE. 2 ) .AND. ( num_sal .NE. 4 ) ) THEN
431         DO ji = kideb, kiut
432            IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN
433               s_i_new(ji)         =  sm_i_b(ji)
434               ! Melting point in K
435               ztmelts             =  - tmut * s_i_new(ji) + rtt 
436               ! New ice heat content (Bitz and Lipscomb, 1999)
437               ztform              =  t_i_b(ji,nlay_i)  ! t_bo_b crashes in the
438               ! Baltic
439               q_i_b(ji,nlay_i+1)  =  rhoic * &
440                  ( cpic * ( ztmelts - ztform     )       &
441                  + lfus * ( 1.0 - ( ztmelts - rtt ) /    &
442                  ( ztform     - rtt ) )       &
443                  - rcp * ( ztmelts-rtt ) )
444               ! Basal growth rate = - F*dt / q
445               dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + &
446                  qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
447            ENDIF ! heat budget
448         END DO ! ji
449      ENDIF ! num_sal
450
451      !-------------------------------------------------
452      ! 4.1 Basal growth - (b) salinity varying in time
453      !-------------------------------------------------
454      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN
455         ! the growth rate (dh_i_bott) is function of the new ice
456         ! heat content (q_i_b(nlay_i+1)). q_i_b depends on the new ice
457         ! salinity (snewice). snewice depends on dh_i_bott
458         ! it converges quickly, so, no problem
459         ! See Vancoppenolle et al., OM08 for more info on this
460
461         ! Initial value (tested 1D, can be anything between 1 and 20)
462         num_iter_max = 4
463         s_i_new(:) = 4.0
464
465         ! Iterative procedure
466         DO iter = 1, num_iter_max
467            DO ji = kideb, kiut
468               IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN
469                  zji = MOD( npb(ji) - 1, jpi ) + 1
470                  zjj = ( npb(ji) - 1 ) / jpi + 1
471                  ! Melting point in K
472                  ztmelts             =   - tmut * s_i_new(ji) + rtt 
473                  ! New ice heat content (Bitz and Lipscomb, 1999)
474                  q_i_b(ji,nlay_i+1)  =  rhoic * &
475                     ( cpic * ( ztmelts - t_bo_b(ji) )       &
476                     + lfus * ( 1.0 - ( ztmelts - rtt ) /    &
477                     ( t_bo_b(ji) - rtt ) )       &
478                     - rcp * ( ztmelts-rtt ) )
479                  ! Bottom growth rate = - F*dt / q
480                  dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) &
481                     + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)
482                  ! New ice salinity ( Cox and Weeks, JGR, 1988 )
483                  ! zswi2  (1) if dh_i_bott/rdt .GT. 3.6e-7
484                  ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8
485                  ! zswi1  (1) if dh_i_bott/rdt .LT. 2.0e-8
486                  zgrr   = MIN( 1.0e-3, MAX ( dh_i_bott(ji) / rdt_ice , zeps ) )
487                  zswi2  = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) ) 
488                  zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
489                  zswi1  = 1. - zswi2 * zswi12 
490                  zfracs = zswi1  * 0.12 +  &
491                     zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) + &
492                     zswi2  * 0.26 /  &
493                     ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 
494                  zds         = zfracs*sss_m(zji,zjj) - s_i_new(ji)
495                  s_i_new(ji) = zfracs * sss_m(zji,zjj)
496               ENDIF ! fc_bo_i
497            END DO ! ji
498         END DO ! iter
499
500         ! Final values
501         DO ji = kideb, kiut
502            IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN
503               ! New ice salinity must not exceed 15 psu
504               s_i_new(ji) = MIN( s_i_new(ji), s_i_max )
505               ! Metling point in K
506               ztmelts     =   - tmut * s_i_new(ji) + rtt 
507               ! New ice heat content (Bitz and Lipscomb, 1999)
508               q_i_b(ji,nlay_i+1)  =  rhoic *                              &
509                  ( cpic * ( ztmelts - t_bo_b(ji) )       &
510                  + lfus * ( 1.0 - ( ztmelts - rtt ) /    &
511                  ( t_bo_b(ji) - rtt ) )       &
512                  - rcp * ( ztmelts-rtt ) )
513               ! Basal growth rate = - F*dt / q
514               dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + &
515                  qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
516               ! Salinity update
517               ! entrapment during bottom growth
518               dsm_i_se_1d(ji) =  ( s_i_new(ji)*dh_i_bott(ji) +              & 
519                  sm_i_b(ji)*ht_i_b(ji) ) /                 & 
520                  MAX( ht_i_b(ji) + dh_i_bott(ji) ,zeps )   &
521                  - sm_i_b(ji)
522            ENDIF ! heat budget
523         END DO ! ji
524      ENDIF ! num_sal
525
526      !----------------
527      ! 4.2 Basal melt
528      !----------------
529      meance_dh = 0.0
530      numce_dh = 0
531      innermelt(:) = 0
532
533      DO ji = kideb, kiut
534         ! heat convergence at the surface > 0
535         IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN
536
537            s_i_new(ji)   =  s_i_b(ji,nlay_i)
538            zqfont_bo(ji) =  rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) )
539
540            zfbase(ji)    =  zqfont_bo(ji) / rdt_ice ! heat conservation test
541            zdq_i(ji)     =  0.0
542
543            dh_i_bott(ji) =  0.0
544         ENDIF
545      END DO
546
547      DO jk = nlay_i, 1, -1
548         DO ji = kideb, kiut
549            IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN
550               ztmelts             =   - tmut * s_i_b(ji,jk) + rtt 
551               IF ( t_i_b(ji,jk) .GE. ztmelts ) THEN
552                  zdeltah(ji,jk)  = - zh_i(ji)
553                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)
554                  innermelt(ji)   = 1
555               ELSE  ! normal ablation
556                  zdeltah(ji,jk)  = - zqfont_bo(ji) / q_i_b(ji,jk)
557                  zqfont_bo(ji)   = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * &
558                     q_i_b(ji,jk)
559                  zdeltah(ji,jk)  = MAX(zdeltah(ji,jk), - zh_i(ji) )
560                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)
561                  zdq_i(ji)       = zdq_i(ji) + zdeltah(ji,jk) * &
562                     q_i_b(ji,jk) / rdt_ice
563                  ! contribution to salt flux
564                  zji             = MOD( npb(ji) - 1, jpi ) + 1
565                  zjj             = ( npb(ji) - 1 ) / jpi + 1
566                  zfsalt_melt(ji) = zfsalt_melt(ji) +                         &
567                     ( sss_m(zji,zjj) - sm_i_b(ji)   ) *        &
568                     a_i_b(ji) * &
569                     MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 
570               ENDIF
571            ENDIF
572         END DO ! ji
573      END DO ! jk
574
575      !-------------------
576      ! Conservation test
577      !-------------------
578      IF ( con_i ) THEN
579         DO ji = kideb, kiut
580            IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN
581               IF ( ( zfbase(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN
582                  numce_dh = numce_dh + 1
583                  meance_dh = meance_dh + zfbase(ji) + zdq_i(ji)
584               ENDIF
585               IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN
586                  WRITE(numout,*) ' ALERTE heat loss for basal  melt '
587                  WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl
588                  WRITE(numout,*) ' ht_i_b  : ', ht_i_b(ji)
589                  WRITE(numout,*) ' zfbase  : ', zfbase(ji)
590                  WRITE(numout,*) ' zdq_i   : ', zdq_i(ji)
591                  WRITE(numout,*) ' ht_i_b  : ', ht_i_b(ji)
592                  WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji)
593                  WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji)
594                  WRITE(numout,*) ' qlbbq_1d: ', qlbbq_1d(ji)
595                  WRITE(numout,*) ' s_i_new : ', s_i_new(ji)
596                  WRITE(numout,*) ' sss_m   : ', sss_m(zji,zjj)
597                  WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
598                  WRITE(numout,*) ' innermelt : ', innermelt(ji)
599               ENDIF
600            ENDIF ! heat convergence at the surface
601         END DO ! ji
602
603         IF ( numce_dh .GT. 0 ) meance_dh = meance_dh / numce_dh
604         WRITE(numout,*) ' Number of points where there is bas. me. error : ', numce_dh
605         WRITE(numout,*) ' Mean basal melt error on error points : ', meance_dh
606         WRITE(numout,*) ' Remaining bottom heat : ', zqfont_bo(jiindex_1d)
607
608      ENDIF ! con_i
609
610      !
611      !------------------------------------------------------------------------------!
612      !  5) Pathological cases                                                       !
613      !------------------------------------------------------------------------------!
614      !
615      !----------------------------------------------
616      ! 5.1 Excessive ablation in a 1-category model
617      !----------------------------------------------
618
619      DO ji = kideb, kiut
620         ! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15)
621         zdhbf = dh_i_bott(ji) 
622         IF (jpl.EQ.1) zdhbf = MAX( hmelt , dh_i_bott(ji) )
623         ! excessive energy is sent to lateral ablation
624         fsup(ji)            =  rhoic*lfus * at_i_b(ji) / MAX( ( 1.0 - at_i_b(ji) ),epsi13) &
625            *  ( zdhbf - dh_i_bott(ji) ) / rdt_ice
626
627         dh_i_bott(ji)  = zdhbf
628         !since ice volume is only used for outputs, we keep it global for all categories
629         dvbbq_1d(ji)   = a_i_b(ji)*dh_i_bott(ji)
630         !new ice thickness
631         zhgnew(ji)     = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji)
632
633         ! diagnostic ( bottom ice growth )
634         zji                 = MOD( npb(ji) - 1, jpi ) + 1
635         zjj                 = ( npb(ji) - 1 ) / jpi + 1
636         diag_bot_gr(zji,zjj) = diag_bot_gr(zji,zjj) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) &
637            / rdt_ice
638         diag_sur_me(zji,zjj) = diag_sur_me(zji,zjj) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) &
639            / rdt_ice
640         diag_bot_me(zji,zjj) = diag_bot_me(zji,zjj) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) &
641            / rdt_ice
642      END DO
643
644      !-----------------------------------
645      ! 5.2 More than available ice melts
646      !-----------------------------------
647      ! then heat applied minus heat content at previous time step
648      ! should equal heat remaining
649      !
650      DO ji = kideb, kiut
651         ! Adapt the remaining energy if too much ice melts
652         !--------------------------------------------------
653         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice
654         ! 0 if no more ice
655         zhgnew(ji) =  zihgnew * zhgnew(ji) ! ice thickness is put to 0
656         ! remaining heat
657         zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) +  zqfont_bo(ji) ) 
658
659         ! If snow remains, energy is used to melt snow
660         zhni       =  ht_s_b(ji) ! snow depth at previous time step
661         zihg       =  MAX( zzero , SIGN ( zone , - ht_s_b(ji) ) ) ! 0 if snow
662
663         ! energy of melting of remaining snow
664         zqt_s(ji)  =  ( 1. - zihg) * zqt_s(ji) / MAX( zhni, zeps )
665         zdhnm      =  - ( 1. - zihg ) * ( 1. - zihgnew ) * ( zfdt_final(ji) /  &
666            MAX( zqt_s(ji) , zeps ) )
667         zhnfi          =  zhni + zdhnm
668         zfdt_final(ji) =  MAX ( zfdt_final(ji) + zqt_s(ji) * zdhnm , 0.0 )
669         ht_s_b(ji)     =  MAX( zzero , zhnfi )
670         zqt_s(ji)      =  zqt_s(ji) * ht_s_b(ji)
671
672         ! Mass variations of ice and snow
673         !---------------------------------
674         !
675         zzfmass_s = - a_i_b(ji) * ( zhni       - ht_s_b(ji) ) * rhosn   ! snow
676         zzfmass_i =   a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic   ! ice 
677         !
678         zfmass_i(ji) = zzfmass_i                       ! ice variation saved to compute salt flux (see below)
679         !
680         !                                              ! mass variation cumulated over category
681         rdmsnif_1d(ji) = rdmsnif_1d(ji) + zzfmass_s                     ! snow
682         rdmicif_1d(ji) = rdmicif_1d(ji) + zzfmass_i                     ! ice
683
684         ! Remaining heat to the ocean
685         !---------------------------------
686         ! focea is in W.m-2 * dt
687         focea(ji)  = - zfdt_final(ji) / rdt_ice
688
689      END DO
690
691      ftotal_fin (:) = zfdt_final(:)  / rdt_ice
692
693      !---------------------------
694      ! Salt flux and heat fluxes                   
695      !---------------------------
696      DO ji = kideb, kiut
697         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice
698
699         ! Salt flux
700         zji           = MOD( npb(ji) - 1, jpi ) + 1
701         zjj           = ( npb(ji) - 1 ) / jpi + 1
702         IF ( num_sal .NE. 4 ) &
703            fseqv_1d(ji)  = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) +           &
704            (1.0 - zihgnew) * zfmass_i(ji) *                  &
705            ( sss_m(zji,zjj) - sm_i_b(ji) ) / rdt_ice
706         ! new lines
707         IF ( num_sal .EQ. 4 ) &
708            fseqv_1d(ji)  = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) +           &
709            (1.0 - zihgnew) * zfmass_i(ji) *                  &
710            ( sss_m(zji,zjj) - bulk_sal ) / rdt_ice
711         ! Heat flux
712         ! excessive bottom ablation energy (fsup) - 0 except if jpl = 1
713         ! excessive total ablation energy (focea) sent to the ocean
714         qfvbq_1d(ji)  = qfvbq_1d(ji) + &
715            fsup(ji) + ( 1.0 - zihgnew ) *        & 
716            focea(ji) * a_i_b(ji) * rdt_ice
717
718         zihic   = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) )
719         ! equals 0 if ht_i = 0, 1 if ht_i gt 0
720         fscbq_1d(ji) =  a_i_b(ji) * fstbif_1d(ji)
721         qldif_1d(ji)  = qldif_1d(ji)                                         &
722            + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) & 
723            * rdt_ice                               &
724            + ( 1.0 - zihic ) * fscbq_1d(ji) * rdt_ice
725      END DO  ! ji
726
727      !-------------------------------------------
728      ! Correct temperature, energy and thickness
729      !-------------------------------------------
730      DO ji = kideb, kiut
731         zihgnew        =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 
732         t_su_b(ji)     =  zihgnew * t_su_b(ji) + ( 1.0 - zihgnew ) * rtt
733      END DO  ! ji
734
735      DO jk = 1, nlay_i
736         DO ji = kideb, kiut
737            zihgnew       =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 
738            t_i_b(ji,jk)  =  zihgnew * t_i_b(ji,jk) + ( 1.0 - zihgnew ) * rtt
739            q_i_b(ji,jk)  =  zihgnew * q_i_b(ji,jk)
740         END DO
741      END DO  ! ji
742
743      DO ji = kideb, kiut
744         ht_i_b(ji) = zhgnew(ji)
745      END DO  ! ji
746      !
747      !------------------------------------------------------------------------------|
748      !  6) Snow-Ice formation                                                       |
749      !------------------------------------------------------------------------------|
750      !
751      ! When snow load excesses Archimede's limit, snow-ice interface goes down
752      ! under sea-level, flooding of seawater transforms snow into ice
753      ! dh_snowice is positive for the ice
754      DO ji = kideb, kiut
755
756         dh_snowice(ji) = MAX(zzero,(rhosn*ht_s_b(ji)+(rhoic-rau0) &
757            * ht_i_b(ji))/(rhosn+rau0-rhoic))
758         zhgnew(ji)     = MAX(zhgnew(ji),zhgnew(ji)+dh_snowice(ji))
759         zhnnew            = MIN(ht_s_b(ji),ht_s_b(ji)-dh_snowice(ji))
760
761         !  Changes in ice volume and ice mass.
762         dvnbq_1d(ji)      = a_i_b(ji) * (zhgnew(ji)-ht_i_b(ji))
763         dmgwi_1d(ji)      = dmgwi_1d(ji) + a_i_b(ji) &
764            *(ht_s_b(ji)-zhnnew)*rhosn
765
766         rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) & 
767            * ( zhgnew(ji) - ht_i_b(ji) )*rhoic 
768         rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) & 
769            * ( zhnnew       - ht_s_b(ji) )*rhosn
770
771         !        Equivalent salt flux (1) Snow-ice formation component
772         !        -----------------------------------------------------
773         zji                 = MOD( npb(ji) - 1, jpi ) + 1
774         zjj                 = ( npb(ji) - 1 ) / jpi + 1
775
776         zsm_snowice  = ( rhoic - rhosn ) / rhoic *            &
777            sss_m(zji,zjj) 
778
779         IF ( num_sal .NE. 2 ) zsm_snowice = sm_i_b(ji)
780
781         IF ( num_sal .NE. 4 ) &
782            fseqv_1d(ji)   = fseqv_1d(ji)   + &
783            ( sss_m(zji,zjj) - zsm_snowice ) * &
784            a_i_b(ji)   * &
785            ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice
786         ! new lines
787         IF ( num_sal .EQ. 4 ) &
788            fseqv_1d(ji)   = fseqv_1d(ji)   + &
789            ( sss_m(zji,zjj) - bulk_sal    ) * &
790            a_i_b(ji)   * &
791            ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice
792
793         ! entrapment during snow ice formation
794         i_ice_switch = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i_b(ji) + 1.0e-6 ) )
795         isnowic      = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - dh_snowice(ji) ) ) * &
796            i_ice_switch
797         IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
798            dsm_i_si_1d(ji)  = ( zsm_snowice*dh_snowice(ji) &
799            + sm_i_b(ji) * ht_i_b(ji)                          & 
800            / MAX( ht_i_b(ji) + dh_snowice(ji), zeps)          &
801            - sm_i_b(ji) ) * isnowic     
802
803         !  Actualize new snow and ice thickness.
804         ht_s_b(ji)  = zhnnew
805         ht_i_b(ji)  = zhgnew(ji)
806
807         ! Total ablation ! new lines added to debug
808         IF( ht_i_b(ji).LE.0.0 ) a_i_b(ji) = 0.0
809
810         ! diagnostic ( snow ice growth )
811         zji                 = MOD( npb(ji) - 1, jpi ) + 1
812         zjj                 = ( npb(ji) - 1 ) / jpi + 1
813         diag_sni_gr(zji,zjj)  = diag_sni_gr(zji,zjj) + dh_snowice(ji)*a_i_b(ji) / &
814            rdt_ice
815
816      END DO !ji
817
818   END SUBROUTINE lim_thd_dh
819#else
820   !!======================================================================
821   !!                       ***  MODULE limthd_dh   ***
822   !!                              no sea ice model 
823   !!======================================================================
824CONTAINS
825   SUBROUTINE lim_thd_dh          ! Empty routine
826   END SUBROUTINE lim_thd_dh
827#endif
828END MODULE limthd_dh
Note: See TracBrowser for help on using the repository browser.