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

Last change on this file since 888 was 888, checked in by ctlod, 16 years ago

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