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

Last change on this file since 1156 was 1156, checked in by rblod, 16 years ago

Update Id and licence information, see ticket #210
  • Property svn:keywords set to Id
File size: 34.8 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 ice_oce          ! ice variables
18   USE sbc_oce          ! Surface boundary condition: ocean fields
19   USE thd_ice
20   USE iceini
21   USE limistate
22   USE in_out_manager
23   USE ice
24   USE par_ice
25   USE lib_mpp
26
27   IMPLICIT NONE
28   PRIVATE
29
30   !! * Routine accessibility
31   PUBLIC lim_thd_dh        ! called by lim_thd
32
33   !! * Module variables
34   REAL(wp)  ::           &  ! constant values
35      epsi20 = 1e-20   ,  &
36      epsi13 = 1e-13   ,  &
37      epsi16 = 1e-16   ,  &
38      zzero  = 0.e0    ,  &
39      zone   = 1.e0
40
41   !!----------------------------------------------------------------------
42   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2005)
43   !! $Id$
44   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46
47CONTAINS
48
49   SUBROUTINE lim_thd_dh(kideb,kiut,jl)
50      !!------------------------------------------------------------------
51      !!                ***  ROUTINE lim_thd_dh  ***
52      !!------------------------------------------------------------------
53      !! ** Purpose :
54      !!           This routine determines variations of ice and snow thicknesses.
55      !! ** Method  :
56      !!           Ice/Snow surface melting arises from imbalance in surface fluxes
57      !!           Bottom accretion/ablation arises from flux budget
58      !!           Snow thickness can increase by precipitation and decrease by
59      !!              sublimation
60      !!           If snow load excesses Archmiede limit, snow-ice is formed by
61      !!              the flooding of sea-water in the snow
62      !! ** Steps 
63      !!           1) Compute available flux of heat for surface ablation
64      !!           2) Compute snow and sea ice enthalpies
65      !!           3) Surface ablation and sublimation
66      !!           4) Bottom accretion/ablation
67      !!           5) Case of Total ablation
68      !!           6) Snow ice formation
69      !!
70      !! ** Arguments
71      !!
72      !! ** Inputs / Outputs
73      !!
74      !! ** External
75      !!
76      !! ** References : Bitz and Lipscomb, JGR 99
77      !!                 Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
78      !!                 Vancoppenolle, Fichefet and Bitz, GRL 2005
79      !!                 Vancoppenolle et al., OM08
80      !!
81      !! ** History  :
82      !!   original code    01-04 (LIM)
83      !!   original routine
84      !!               (05-2003) M. Vancoppenolle, Louvain-La-Neuve, Belgium
85      !!               (06/07-2005) 3D version
86      !!               (03-2008)    Clean code
87      !!
88      !!------------------------------------------------------------------
89      !! * Arguments
90      INTEGER , INTENT (IN) ::  &
91         kideb     ,         &  !: Start point on which the  the computation is applied
92         kiut      ,         &  !: End point on which the  the computation is applied
93         jl                     !: Thickness cateogry number
94
95      !! * Local variables
96      INTEGER ::             & 
97         ji        ,         &  !: space index
98         jk        ,         &  !: ice layer index
99         isnow     ,         &  !: switch for presence (1) or absence (0) of snow
100         zji       ,         &  !: 2D corresponding indices to ji
101         zjj       ,         &
102         isnowic   ,         &  !: snow ice formation not
103         i_ice_switch   ,    &  !: ice thickness above a certain treshold or not
104         iter
105
106      REAL(wp) ::            &
107         zhsnew    ,         &  !: new snow thickness
108         zihgnew   ,         &  !: switch for total ablation
109         ztmelts   ,         &  !: melting point
110         zhn       ,         &
111         zdhcf     ,         &
112         zdhbf     ,         &
113         zhni      ,         &
114         zhnfi     ,         &
115         zihg      ,         &
116         zdhnm     ,         &
117         zhnnew    ,         &
118         zeps = 1.0e-13,     &
119         zhisn     ,         &
120         zfracs    ,         &  !: fractionation coefficient for bottom salt
121                                !: entrapment
122         zds       ,         &  !: increment of bottom ice salinity
123         zcoeff    ,         &  !: dummy argument for snowfall partitioning
124                                !: over ice and leads
125         zsm_snowice,        &  !: snow-ice salinity
126         zswi1     ,         &  !: switch for computation of bottom salinity
127         zswi12    ,         &  !: switch for computation of bottom salinity
128         zswi2     ,         &  !: switch for computation of bottom salinity
129         zgrr      ,         &  !: bottom growth rate
130         zihic     ,         &  !:
131         ztform                 !: bottom formation temperature
132
133      REAL(wp) , DIMENSION(jpij) ::  &
134         zh_i      ,         &  ! ice layer thickness
135         zh_s      ,         &  ! snow layer thickness
136         ztfs      ,         &  ! melting point
137         zhsold    ,         &  ! old snow thickness
138         zqprec    ,         &  !: energy of fallen snow
139         zqfont_su ,         &  ! incoming, remaining surface energy
140         zqfont_bo              ! incoming, bottom energy
141
142      REAL(wp) , DIMENSION(jpij) :: &
143         z_f_surf,           &  ! surface heat for ablation
144         zhgnew                 ! new ice thickness
145
146      REAL(wp), DIMENSION(jpij) :: &
147         zdh_s_mel  ,        &  ! snow melt
148         zdh_s_pre  ,        &  ! snow precipitation
149         zdh_s_sub  ,        &  ! snow sublimation
150         zfsalt_melt            ! salt flux due to ice melt
151
152      REAL(wp) , DIMENSION(jpij,jkmax) :: &
153         zdeltah
154
155      ! Pathological cases
156      REAL(wp), DIMENSION(jpij) :: &
157         zfdt_init  ,        &  !: total incoming heat for ice melt
158         zfdt_final ,        &  !: total remaing heat for ice melt
159         zqt_i      ,        &  !: total ice heat content
160         zqt_s      ,        &  !: total snow heat content
161         zqt_dummy              !: dummy heat content
162
163      REAL(wp), DIMENSION(jpij,jkmax) :: &
164         zqt_i_lay              !: total ice heat content
165
166      ! Heat conservation
167      REAL(wp), DIMENSION(jpij) :: &
168         zfbase,             &
169         zdq_i
170
171      INTEGER, DIMENSION(jpij) ::  &
172         innermelt
173
174      REAL(wp) :: &
175         meance_dh
176
177      INTEGER ::                   &
178         num_iter_max,       &
179         numce_dh
180
181      zfsalt_melt(:)  = 0.0
182      ftotal_fin(:)   = 0.0
183      zfdt_init(:)    = 0.0
184      zfdt_final(:)   = 0.0
185
186      DO ji = kideb, kiut
187         old_ht_i_b(ji) = ht_i_b(ji)
188         old_ht_s_b(ji) = ht_s_b(ji)
189      END DO
190      !
191      !------------------------------------------------------------------------------!
192      !  1) Calculate available heat for surface ablation                            !
193      !------------------------------------------------------------------------------!
194      !
195      DO ji = kideb,kiut
196         isnow         = INT( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ) )
197         ztfs(ji)      = isnow * rtt + ( 1.0 - isnow ) * rtt
198         z_f_surf(ji)  = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) *                 & 
199            qsr_ice_1d(ji) - fc_su(ji)
200         z_f_surf(ji)  = MAX( zzero , z_f_surf(ji) ) *                        &
201            MAX( zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) )
202         zfdt_init(ji) = ( z_f_surf(ji) + &
203            MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) )  &
204            * rdt_ice
205      END DO ! ji
206
207      zqfont_su(:) = 0.0
208      zqfont_bo(:) = 0.0
209      dsm_i_se_1d(:) = 0.0     
210      dsm_i_si_1d(:) = 0.0     
211      !
212      !------------------------------------------------------------------------------!
213      !  2) Computing layer thicknesses and  snow and sea-ice enthalpies.            !
214      !------------------------------------------------------------------------------!
215      !
216      ! Layer thickness
217      DO ji = kideb,kiut
218         zh_i(ji) = ht_i_b(ji) / nlay_i
219         zh_s(ji) = ht_s_b(ji) / nlay_s
220      END DO
221
222      ! Total enthalpy of the snow
223      zqt_s(:) = 0.0
224      DO jk = 1, nlay_s
225         DO ji = kideb,kiut
226            zqt_s(ji)    =  zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s
227         END DO
228      END DO
229
230      ! Total enthalpy of the ice
231      zqt_i(:) = 0.0
232      DO jk = 1, nlay_i
233         DO ji = kideb,kiut
234            zqt_i(ji)        =  zqt_i(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
235            zqt_i_lay(ji,jk) =  q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
236         END DO
237      END DO
238      !
239      !------------------------------------------------------------------------------|
240      !  3) Surface ablation and sublimation                                         |
241      !------------------------------------------------------------------------------|
242      !
243      !-------------------------
244      ! 3.1 Snow precips / melt
245      !-------------------------
246      ! Snow accumulation in one thermodynamic time step
247      ! snowfall is partitionned between leads and ice
248      ! if snow fall was uniform, a fraction (1-at_i) would fall into leads
249      ! but because of the winds, more snow falls on leads than on sea ice
250      ! and a greater fraction (1-at_i)^beta of the total mass of snow
251      ! (beta < 1) falls in leads.
252      ! In reality, beta depends on wind speed,
253      ! and should decrease with increasing wind speed but here, it is
254      ! considered as a constant. an average value is 0.66
255      ! Martin Vancoppenolle, December 2006
256
257      ! Snow fall
258      DO ji = kideb, kiut
259         zcoeff = ( 1.0 - ( 1.0 - at_i_b(ji) )**betas ) / at_i_b(ji) 
260         zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn
261      END DO
262      zdh_s_mel(:) =  0.0
263
264      ! Melt of fallen snow
265      DO ji = kideb, kiut
266         ! tatm_ice is now in K
267         zqprec(ji)     =  rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus ) 
268         zqfont_su(ji)  =  z_f_surf(ji) * rdt_ice
269         zdeltah(ji,1)  =  MIN( 0.0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) )
270         zqfont_su(ji)  =  MAX( 0.0 , - zdh_s_pre(ji) - zdeltah(ji,1) )      * &
271            zqprec(ji)
272         zdeltah(ji,1)  =  MAX( - zdh_s_pre(ji) , zdeltah(ji,1) )
273         zdh_s_mel(ji)  =  zdh_s_mel(ji) + zdeltah(ji,1)
274         ! heat conservation
275         qt_s_in(ji,jl) =  qt_s_in(ji,jl) + zqprec(ji) * zdh_s_pre(ji)
276         zqt_s(ji)      =  zqt_s(ji) + zqprec(ji) * zdh_s_pre(ji)
277      END DO
278
279      ! Update total snow heat content
280      zqt_s(ji)         =  MAX ( zqt_s(ji) - zqfont_su(ji) , 0.0 ) 
281      IF( lk_mpp ) CALL mpp_max(zqt_s(ji), kcom = ncomm_ice ) 
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) =  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         rdmicif_1d(ji) =  rdmicif_1d(ji) + a_i_b(ji) *                        &
675            (zhgnew(ji)-ht_i_b(ji))*rhoic ! good
676
677         rdmsnif_1d(ji) =  rdmsnif_1d(ji) + a_i_b(ji) * &
678            (ht_s_b(ji)-zhni)*rhosn ! good too
679
680         ! Remaining heat to the ocean
681         !---------------------------------
682         ! focea is in W.m-2 * dt
683         focea(ji)  = - zfdt_final(ji) / rdt_ice
684
685      END DO
686
687      ftotal_fin (:) = zfdt_final(:)  / rdt_ice
688
689      !---------------------------
690      ! Salt flux and heat fluxes                   
691      !---------------------------
692      DO ji = kideb, kiut
693         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice
694
695         ! Salt flux
696         zji           = MOD( npb(ji) - 1, jpi ) + 1
697         zjj           = ( npb(ji) - 1 ) / jpi + 1
698         IF ( num_sal .NE. 4 ) &
699            fseqv_1d(ji)  = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) +           &
700            (1.0 - zihgnew) * rdmicif_1d(ji) *                  &
701            ( sss_m(zji,zjj) - sm_i_b(ji) ) / rdt_ice
702         ! new lines
703         IF ( num_sal .EQ. 4 ) &
704            fseqv_1d(ji)  = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) +           &
705            (1.0 - zihgnew) * rdmicif_1d(ji) *                  &
706            ( sss_m(zji,zjj) - bulk_sal ) / rdt_ice
707         ! Heat flux
708         ! excessive bottom ablation energy (fsup) - 0 except if jpl = 1
709         ! excessive total ablation energy (focea) sent to the ocean
710         qfvbq_1d(ji)  = qfvbq_1d(ji) + &
711            fsup(ji) + ( 1.0 - zihgnew ) *        & 
712            focea(ji) * a_i_b(ji) * rdt_ice
713
714         zihic   = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) )
715         ! equals 0 if ht_i = 0, 1 if ht_i gt 0
716         fscbq_1d(ji) =  a_i_b(ji) * fstbif_1d(ji)
717         qldif_1d(ji)  = qldif_1d(ji)                                         &
718            + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) & 
719            * rdt_ice                               &
720            + ( 1.0 - zihic ) * fscbq_1d(ji) * rdt_ice
721      END DO  ! ji
722
723      !-------------------------------------------
724      ! Correct temperature, energy and thickness
725      !-------------------------------------------
726      DO ji = kideb, kiut
727         zihgnew        =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 
728         t_su_b(ji)     =  zihgnew * t_su_b(ji) + ( 1.0 - zihgnew ) * rtt
729      END DO  ! ji
730
731      DO jk = 1, nlay_i
732         DO ji = kideb, kiut
733            zihgnew       =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 
734            t_i_b(ji,jk)  =  zihgnew * t_i_b(ji,jk) + ( 1.0 - zihgnew ) * rtt
735            q_i_b(ji,jk)  =  zihgnew * q_i_b(ji,jk)
736         END DO
737      END DO  ! ji
738
739      DO ji = kideb, kiut
740         ht_i_b(ji) = zhgnew(ji)
741      END DO  ! ji
742      !
743      !------------------------------------------------------------------------------|
744      !  6) Snow-Ice formation                                                       |
745      !------------------------------------------------------------------------------|
746      !
747      ! When snow load excesses Archimede's limit, snow-ice interface goes down
748      ! under sea-level, flooding of seawater transforms snow into ice
749      ! dh_snowice is positive for the ice
750      DO ji = kideb, kiut
751
752         dh_snowice(ji) = MAX(zzero,(rhosn*ht_s_b(ji)+(rhoic-rau0) &
753            * ht_i_b(ji))/(rhosn+rau0-rhoic))
754         zhgnew(ji)     = MAX(zhgnew(ji),zhgnew(ji)+dh_snowice(ji))
755         zhnnew            = MIN(ht_s_b(ji),ht_s_b(ji)-dh_snowice(ji))
756
757         !  Changes in ice volume and ice mass.
758         dvnbq_1d(ji)      = a_i_b(ji) * (zhgnew(ji)-ht_i_b(ji))
759         dmgwi_1d(ji)      = dmgwi_1d(ji) + a_i_b(ji) &
760            *(ht_s_b(ji)-zhnnew)*rhosn
761
762         rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) & 
763            * ( zhgnew(ji) - ht_i_b(ji) )*rhoic 
764         rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) & 
765            * ( zhnnew       - ht_s_b(ji) )*rhosn
766
767         !        Equivalent salt flux (1) Snow-ice formation component
768         !        -----------------------------------------------------
769         zji                 = MOD( npb(ji) - 1, jpi ) + 1
770         zjj                 = ( npb(ji) - 1 ) / jpi + 1
771
772         zsm_snowice  = ( rhoic - rhosn ) / rhoic *            &
773            sss_m(zji,zjj) 
774
775         IF ( num_sal .NE. 2 ) zsm_snowice = sm_i_b(ji)
776
777         IF ( num_sal .NE. 4 ) &
778            fseqv_1d(ji)   = fseqv_1d(ji)   + &
779            ( sss_m(zji,zjj) - zsm_snowice ) * &
780            a_i_b(ji)   * &
781            ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice
782         ! new lines
783         IF ( num_sal .EQ. 4 ) &
784            fseqv_1d(ji)   = fseqv_1d(ji)   + &
785            ( sss_m(zji,zjj) - bulk_sal    ) * &
786            a_i_b(ji)   * &
787            ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice
788
789         ! entrapment during snow ice formation
790         i_ice_switch = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i_b(ji) + 1.0e-6 ) )
791         isnowic      = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - dh_snowice(ji) ) ) * &
792            i_ice_switch
793         IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
794            dsm_i_si_1d(ji)  = ( zsm_snowice*dh_snowice(ji) &
795            + sm_i_b(ji) * ht_i_b(ji)                          & 
796            / MAX( ht_i_b(ji) + dh_snowice(ji), zeps)          &
797            - sm_i_b(ji) ) * isnowic     
798
799         !  Actualize new snow and ice thickness.
800         ht_s_b(ji)  = zhnnew
801         ht_i_b(ji)  = zhgnew(ji)
802
803         ! Total ablation ! new lines added to debug
804         IF( ht_i_b(ji).LE.0.0 ) a_i_b(ji) = 0.0
805
806         ! diagnostic ( snow ice growth )
807         zji                 = MOD( npb(ji) - 1, jpi ) + 1
808         zjj                 = ( npb(ji) - 1 ) / jpi + 1
809         diag_sni_gr(zji,zjj)  = diag_sni_gr(zji,zjj) + dh_snowice(ji)*a_i_b(ji) / &
810            rdt_ice
811
812      END DO !ji
813
814   END SUBROUTINE lim_thd_dh
815#else
816   !!======================================================================
817   !!                       ***  MODULE limthd_dh   ***
818   !!                              no sea ice model 
819   !!======================================================================
820CONTAINS
821   SUBROUTINE lim_thd_dh          ! Empty routine
822   END SUBROUTINE lim_thd_dh
823#endif
824END MODULE limthd_dh
Note: See TracBrowser for help on using the repository browser.