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

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

now input air temperature data for CORE & CLIO bulks must be in Kelvin, see ticket: #242
  • Property svn:keywords set to Id
File size: 34.7 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         zqt_s(ji)      =  MAX ( zqt_s(ji) - zqfont_su(ji) , 0.0 ) 
278      END DO
279
280
281      ! Snow melt due to surface heat imbalance
282      DO jk = 1, nlay_s
283         DO ji = kideb, kiut
284            zdeltah(ji,jk) = - zqfont_su(ji) / q_s_b(ji,jk)
285            zqfont_su(ji)  =  MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * &
286               q_s_b(ji,jk) 
287            zdeltah(ji,jk) =  MAX( zdeltah(ji,jk) , - zh_s(ji) )
288            zdh_s_mel(ji)  =  zdh_s_mel(ji) + zdeltah(ji,jk) !resulting melt of snow   
289         END DO
290      END DO
291
292      ! Apply snow melt to snow depth
293      DO ji = kideb, kiut
294         dh_s_tot(ji)   =  zdh_s_mel(ji) + zdh_s_pre(ji)
295         ! Old and new snow depths
296         zhsold(ji)     =  ht_s_b(ji)
297         zhsnew         =  ht_s_b(ji) + dh_s_tot(ji)
298         ! If snow is still present zhn = 1, else zhn = 0
299         zhn            =  1.0 - MAX( zzero , SIGN( zone , - zhsnew ) )
300         ht_s_b(ji)     =  MAX( zzero , zhsnew )
301         ! Volume and mass variations of snow
302         dvsbq_1d(ji)   =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji)              &
303            - zdh_s_mel(ji) )
304         dvsbq_1d(ji)   =  MIN( zzero, dvsbq_1d(ji) )
305         rdmsnif_1d(ji) =  rhosn*dvsbq_1d(ji)
306      END DO ! ji
307
308      !--------------------------
309      ! 3.2 Surface ice ablation
310      !--------------------------
311      DO ji = kideb, kiut 
312         dh_i_surf(ji) =  0.0
313         ! For heat conservation test
314         z_f_surf(ji)  =  zqfont_su(ji) / rdt_ice ! heat conservation test
315         zdq_i(ji)     =  0.0
316      END DO ! ji
317
318      DO jk = 1, nlay_i
319         DO ji = kideb, kiut 
320            ! melt of layer jk
321            zdeltah(ji,jk)      = - zqfont_su(ji) / q_i_b(ji,jk)
322            ! recompute heat available
323            zqfont_su(ji)       = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) *  &
324               q_i_b(ji,jk) 
325            ! melt of layer jk cannot be higher than its thickness
326            zdeltah(ji,jk)      = MAX( zdeltah(ji,jk) , - zh_i(ji) )
327            ! update surface melt
328            dh_i_surf(ji)       = dh_i_surf(ji) + zdeltah(ji,jk) 
329            ! for energy conservation
330            zdq_i(ji)           = zdq_i(ji) + zdeltah(ji,jk) *                &
331               q_i_b(ji,jk) / rdt_ice
332            ! contribution to ice-ocean salt flux
333            zji                 = MOD( npb(ji) - 1, jpi ) + 1
334            zjj                 = ( npb(ji) - 1 ) / jpi + 1
335            zfsalt_melt(ji)     = zfsalt_melt(ji) +                           &
336               ( sss_m(zji,zjj) - sm_i_b(ji)   ) *         &
337               a_i_b(ji) *                                 &
338               MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 
339         END DO ! ji
340      END DO ! jk
341
342      !-------------------
343      ! Conservation test
344      !-------------------
345      IF ( con_i ) THEN
346         numce_dh = 0
347         meance_dh = 0.0
348         DO ji = kideb, kiut
349
350            IF ( ( z_f_surf(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN
351               numce_dh  = numce_dh + 1
352               meance_dh = meance_dh + z_f_surf(ji) + zdq_i(ji)
353            ENDIF
354
355            IF ( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN!
356               WRITE(numout,*) ' ALERTE heat loss for surface melt '
357               WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl
358               WRITE(numout,*) ' ht_i_b  : ', ht_i_b(ji)
359               WRITE(numout,*) ' z_f_surf  : ', z_f_surf(ji)
360               WRITE(numout,*) ' zdq_i   : ', zdq_i(ji)
361               WRITE(numout,*) ' ht_i_b  : ', ht_i_b(ji)
362               WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji)
363               WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji)
364               WRITE(numout,*) ' qlbbq_1d: ', qlbbq_1d(ji)
365               WRITE(numout,*) ' s_i_new : ', s_i_new(ji)
366               WRITE(numout,*) ' sss_m   : ', sss_m(zji,zjj)
367            ENDIF
368
369         END DO ! ji
370
371         IF ( numce_dh .GT. 0 ) meance_dh = meance_dh / numce_dh
372         WRITE(numout,*) ' Error report - Category : ', jl
373         WRITE(numout,*) ' ~~~~~~~~~~~~ '
374         WRITE(numout,*) ' Number of points where there is sur. me. error : ', numce_dh
375         WRITE(numout,*) ' Mean basal growth error on error points : ', meance_dh
376
377      ENDIF ! con_i
378
379      !----------------------
380      ! 3.3 Snow sublimation
381      !----------------------
382
383      DO ji = kideb, kiut
384         ! if qla is positive (upwards), heat goes to the atmosphere, therefore
385         ! snow sublimates, if qla is negative (downwards), snow condensates
386         zdh_s_sub(ji)   =  - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice
387         dh_s_tot(ji)    =  dh_s_tot(ji) + zdh_s_sub(ji)
388         zdhcf           =  ht_s_b(ji) + zdh_s_sub(ji) 
389         ht_s_b(ji)      =  MAX( zzero , zdhcf )
390         ! we recompute dh_s_tot
391         dh_s_tot(ji)    =  ht_s_b(ji) - zhsold(ji)
392         qt_s_in(ji,jl) =  qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1)
393      END DO  !ji
394
395      zqt_dummy(:) = 0.0
396      DO jk = 1, nlay_s
397         DO ji = kideb,kiut
398            q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
399            ! heat conservation
400            zqt_dummy(ji)    =  zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s
401         END DO
402      END DO
403
404      DO jk = 1, nlay_s !n
405         DO ji = kideb, kiut !n
406            ! In case of disparition of the snow, we have to update the snow
407            ! temperatures
408            zhisn  =  MAX( zzero , SIGN( zone, - ht_s_b(ji) ) )
409            t_s_b(ji,jk) = ( 1.0 - zhisn ) * t_s_b(ji,jk) + zhisn * rtt
410            q_s_b(ji,jk) = ( 1.0 - zhisn ) * q_s_b(ji,jk)
411         END DO
412      END DO
413
414      !
415      !------------------------------------------------------------------------------!
416      ! 4) Basal growth / melt                                                       !
417      !------------------------------------------------------------------------------!
418      !
419      ! Ice basal growth / melt is given by the ratio of heat budget over basal
420      ! ice heat content.  Basal heat budget is given by the difference between
421      ! the inner conductive flux  (fc_bo_i), from the open water heat flux
422      ! (qlbbqb) and the turbulent ocean flux (fbif).
423      ! fc_bo_i is positive downwards. fbif and qlbbq are positive to the ice
424
425      !-----------------------------------------------------
426      ! 4.1 Basal growth - (a) salinity not varying in time
427      !-----------------------------------------------------
428      IF ( ( num_sal .NE. 2 ) .AND. ( num_sal .NE. 4 ) ) THEN
429         DO ji = kideb, kiut
430            IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN
431               s_i_new(ji)         =  sm_i_b(ji)
432               ! Melting point in K
433               ztmelts             =  - tmut * s_i_new(ji) + rtt 
434               ! New ice heat content (Bitz and Lipscomb, 1999)
435               ztform              =  t_i_b(ji,nlay_i)  ! t_bo_b crashes in the
436               ! Baltic
437               q_i_b(ji,nlay_i+1)  =  rhoic * &
438                  ( cpic * ( ztmelts - ztform     )       &
439                  + lfus * ( 1.0 - ( ztmelts - rtt ) /    &
440                  ( ztform     - rtt ) )       &
441                  - rcp * ( ztmelts-rtt ) )
442               ! Basal growth rate = - F*dt / q
443               dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + &
444                  qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
445            ENDIF ! heat budget
446         END DO ! ji
447      ENDIF ! num_sal
448
449      !-------------------------------------------------
450      ! 4.1 Basal growth - (b) salinity varying in time
451      !-------------------------------------------------
452      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN
453         ! the growth rate (dh_i_bott) is function of the new ice
454         ! heat content (q_i_b(nlay_i+1)). q_i_b depends on the new ice
455         ! salinity (snewice). snewice depends on dh_i_bott
456         ! it converges quickly, so, no problem
457         ! See Vancoppenolle et al., OM08 for more info on this
458
459         ! Initial value (tested 1D, can be anything between 1 and 20)
460         num_iter_max = 4
461         s_i_new(:) = 4.0
462
463         ! Iterative procedure
464         DO iter = 1, num_iter_max
465            DO ji = kideb, kiut
466               IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN
467                  zji = MOD( npb(ji) - 1, jpi ) + 1
468                  zjj = ( npb(ji) - 1 ) / jpi + 1
469                  ! Melting point in K
470                  ztmelts             =   - tmut * s_i_new(ji) + rtt 
471                  ! New ice heat content (Bitz and Lipscomb, 1999)
472                  q_i_b(ji,nlay_i+1)  =  rhoic * &
473                     ( cpic * ( ztmelts - t_bo_b(ji) )       &
474                     + lfus * ( 1.0 - ( ztmelts - rtt ) /    &
475                     ( t_bo_b(ji) - rtt ) )       &
476                     - rcp * ( ztmelts-rtt ) )
477                  ! Bottom growth rate = - F*dt / q
478                  dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) &
479                     + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)
480                  ! New ice salinity ( Cox and Weeks, JGR, 1988 )
481                  ! zswi2  (1) if dh_i_bott/rdt .GT. 3.6e-7
482                  ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8
483                  ! zswi1  (1) if dh_i_bott/rdt .LT. 2.0e-8
484                  zgrr   = MIN( 1.0e-3, MAX ( dh_i_bott(ji) / rdt_ice , zeps ) )
485                  zswi2  = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) ) 
486                  zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
487                  zswi1  = 1. - zswi2 * zswi12 
488                  zfracs = zswi1  * 0.12 +  &
489                     zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) + &
490                     zswi2  * 0.26 /  &
491                     ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 
492                  zds         = zfracs*sss_m(zji,zjj) - s_i_new(ji)
493                  s_i_new(ji) = zfracs * sss_m(zji,zjj)
494               ENDIF ! fc_bo_i
495            END DO ! ji
496         END DO ! iter
497
498         ! Final values
499         DO ji = kideb, kiut
500            IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN
501               ! New ice salinity must not exceed 15 psu
502               s_i_new(ji) = MIN( s_i_new(ji), s_i_max )
503               ! Metling point in K
504               ztmelts     =   - tmut * s_i_new(ji) + rtt 
505               ! New ice heat content (Bitz and Lipscomb, 1999)
506               q_i_b(ji,nlay_i+1)  =  rhoic *                              &
507                  ( cpic * ( ztmelts - t_bo_b(ji) )       &
508                  + lfus * ( 1.0 - ( ztmelts - rtt ) /    &
509                  ( t_bo_b(ji) - rtt ) )       &
510                  - rcp * ( ztmelts-rtt ) )
511               ! Basal growth rate = - F*dt / q
512               dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + &
513                  qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
514               ! Salinity update
515               ! entrapment during bottom growth
516               dsm_i_se_1d(ji) =  ( s_i_new(ji)*dh_i_bott(ji) +              & 
517                  sm_i_b(ji)*ht_i_b(ji) ) /                 & 
518                  MAX( ht_i_b(ji) + dh_i_bott(ji) ,zeps )   &
519                  - sm_i_b(ji)
520            ENDIF ! heat budget
521         END DO ! ji
522      ENDIF ! num_sal
523
524      !----------------
525      ! 4.2 Basal melt
526      !----------------
527      meance_dh = 0.0
528      numce_dh = 0
529      innermelt(:) = 0
530
531      DO ji = kideb, kiut
532         ! heat convergence at the surface > 0
533         IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN
534
535            s_i_new(ji)   =  s_i_b(ji,nlay_i)
536            zqfont_bo(ji) =  rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) )
537
538            zfbase(ji)    =  zqfont_bo(ji) / rdt_ice ! heat conservation test
539            zdq_i(ji)     =  0.0
540
541            dh_i_bott(ji) =  0.0
542         ENDIF
543      END DO
544
545      DO jk = nlay_i, 1, -1
546         DO ji = kideb, kiut
547            IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN
548               ztmelts             =   - tmut * s_i_b(ji,jk) + rtt 
549               IF ( t_i_b(ji,jk) .GE. ztmelts ) THEN
550                  zdeltah(ji,jk)  = - zh_i(ji)
551                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)
552                  innermelt(ji)   = 1
553               ELSE  ! normal ablation
554                  zdeltah(ji,jk)  = - zqfont_bo(ji) / q_i_b(ji,jk)
555                  zqfont_bo(ji)   = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * &
556                     q_i_b(ji,jk)
557                  zdeltah(ji,jk)  = MAX(zdeltah(ji,jk), - zh_i(ji) )
558                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)
559                  zdq_i(ji)       = zdq_i(ji) + zdeltah(ji,jk) * &
560                     q_i_b(ji,jk) / rdt_ice
561                  ! contribution to salt flux
562                  zji             = MOD( npb(ji) - 1, jpi ) + 1
563                  zjj             = ( npb(ji) - 1 ) / jpi + 1
564                  zfsalt_melt(ji) = zfsalt_melt(ji) +                         &
565                     ( sss_m(zji,zjj) - sm_i_b(ji)   ) *        &
566                     a_i_b(ji) * &
567                     MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 
568               ENDIF
569            ENDIF
570         END DO ! ji
571      END DO ! jk
572
573      !-------------------
574      ! Conservation test
575      !-------------------
576      IF ( con_i ) THEN
577         DO ji = kideb, kiut
578            IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN
579               IF ( ( zfbase(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN
580                  numce_dh = numce_dh + 1
581                  meance_dh = meance_dh + zfbase(ji) + zdq_i(ji)
582               ENDIF
583               IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN
584                  WRITE(numout,*) ' ALERTE heat loss for basal  melt '
585                  WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl
586                  WRITE(numout,*) ' ht_i_b  : ', ht_i_b(ji)
587                  WRITE(numout,*) ' zfbase  : ', zfbase(ji)
588                  WRITE(numout,*) ' zdq_i   : ', zdq_i(ji)
589                  WRITE(numout,*) ' ht_i_b  : ', ht_i_b(ji)
590                  WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji)
591                  WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji)
592                  WRITE(numout,*) ' qlbbq_1d: ', qlbbq_1d(ji)
593                  WRITE(numout,*) ' s_i_new : ', s_i_new(ji)
594                  WRITE(numout,*) ' sss_m   : ', sss_m(zji,zjj)
595                  WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
596                  WRITE(numout,*) ' innermelt : ', innermelt(ji)
597               ENDIF
598            ENDIF ! heat convergence at the surface
599         END DO ! ji
600
601         IF ( numce_dh .GT. 0 ) meance_dh = meance_dh / numce_dh
602         WRITE(numout,*) ' Number of points where there is bas. me. error : ', numce_dh
603         WRITE(numout,*) ' Mean basal melt error on error points : ', meance_dh
604         WRITE(numout,*) ' Remaining bottom heat : ', zqfont_bo(jiindex_1d)
605
606      ENDIF ! con_i
607
608      !
609      !------------------------------------------------------------------------------!
610      !  5) Pathological cases                                                       !
611      !------------------------------------------------------------------------------!
612      !
613      !----------------------------------------------
614      ! 5.1 Excessive ablation in a 1-category model
615      !----------------------------------------------
616
617      DO ji = kideb, kiut
618         ! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15)
619         zdhbf = dh_i_bott(ji) 
620         IF (jpl.EQ.1) zdhbf = MAX( hmelt , dh_i_bott(ji) )
621         ! excessive energy is sent to lateral ablation
622         fsup(ji)            =  rhoic*lfus * at_i_b(ji) / MAX( ( 1.0 - at_i_b(ji) ),epsi13) &
623            *  ( zdhbf - dh_i_bott(ji) ) / rdt_ice
624
625         dh_i_bott(ji)  = zdhbf
626         !since ice volume is only used for outputs, we keep it global for all categories
627         dvbbq_1d(ji)   = a_i_b(ji)*dh_i_bott(ji)
628         !new ice thickness
629         zhgnew(ji)     = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji)
630
631         ! diagnostic ( bottom ice growth )
632         zji                 = MOD( npb(ji) - 1, jpi ) + 1
633         zjj                 = ( npb(ji) - 1 ) / jpi + 1
634         diag_bot_gr(zji,zjj) = diag_bot_gr(zji,zjj) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) &
635            / rdt_ice
636         diag_sur_me(zji,zjj) = diag_sur_me(zji,zjj) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) &
637            / rdt_ice
638         diag_bot_me(zji,zjj) = diag_bot_me(zji,zjj) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) &
639            / rdt_ice
640      END DO
641
642      !-----------------------------------
643      ! 5.2 More than available ice melts
644      !-----------------------------------
645      ! then heat applied minus heat content at previous time step
646      ! should equal heat remaining
647      !
648      DO ji = kideb, kiut
649         ! Adapt the remaining energy if too much ice melts
650         !--------------------------------------------------
651         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice
652         ! 0 if no more ice
653         zhgnew(ji) =  zihgnew * zhgnew(ji) ! ice thickness is put to 0
654         ! remaining heat
655         zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) +  zqfont_bo(ji) ) 
656
657         ! If snow remains, energy is used to melt snow
658         zhni       =  ht_s_b(ji) ! snow depth at previous time step
659         zihg       =  MAX( zzero , SIGN ( zone , - ht_s_b(ji) ) ) ! 0 if snow
660
661         ! energy of melting of remaining snow
662         zqt_s(ji)  =  ( 1. - zihg) * zqt_s(ji) / MAX( zhni, zeps )
663         zdhnm      =  - ( 1. - zihg ) * ( 1. - zihgnew ) * ( zfdt_final(ji) /  &
664            MAX( zqt_s(ji) , zeps ) )
665         zhnfi          =  zhni + zdhnm
666         zfdt_final(ji) =  MAX ( zfdt_final(ji) + zqt_s(ji) * zdhnm , 0.0 )
667         ht_s_b(ji)     =  MAX( zzero , zhnfi )
668         zqt_s(ji)      =  zqt_s(ji) * ht_s_b(ji)
669
670         ! Mass variations of ice and snow
671         !---------------------------------
672         rdmicif_1d(ji) =  rdmicif_1d(ji) + a_i_b(ji) *                        &
673            (zhgnew(ji)-ht_i_b(ji))*rhoic ! good
674
675         rdmsnif_1d(ji) =  rdmsnif_1d(ji) + a_i_b(ji) * &
676            (ht_s_b(ji)-zhni)*rhosn ! good too
677
678         ! Remaining heat to the ocean
679         !---------------------------------
680         ! focea is in W.m-2 * dt
681         focea(ji)  = - zfdt_final(ji) / rdt_ice
682
683      END DO
684
685      ftotal_fin (:) = zfdt_final(:)  / rdt_ice
686
687      !---------------------------
688      ! Salt flux and heat fluxes                   
689      !---------------------------
690      DO ji = kideb, kiut
691         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice
692
693         ! Salt flux
694         zji           = MOD( npb(ji) - 1, jpi ) + 1
695         zjj           = ( npb(ji) - 1 ) / jpi + 1
696         IF ( num_sal .NE. 4 ) &
697            fseqv_1d(ji)  = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) +           &
698            (1.0 - zihgnew) * rdmicif_1d(ji) *                  &
699            ( sss_m(zji,zjj) - sm_i_b(ji) ) / rdt_ice
700         ! new lines
701         IF ( num_sal .EQ. 4 ) &
702            fseqv_1d(ji)  = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) +           &
703            (1.0 - zihgnew) * rdmicif_1d(ji) *                  &
704            ( sss_m(zji,zjj) - bulk_sal ) / rdt_ice
705         ! Heat flux
706         ! excessive bottom ablation energy (fsup) - 0 except if jpl = 1
707         ! excessive total ablation energy (focea) sent to the ocean
708         qfvbq_1d(ji)  = qfvbq_1d(ji) + &
709            fsup(ji) + ( 1.0 - zihgnew ) *        & 
710            focea(ji) * a_i_b(ji) * rdt_ice
711
712         zihic   = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) )
713         ! equals 0 if ht_i = 0, 1 if ht_i gt 0
714         fscbq_1d(ji) =  a_i_b(ji) * fstbif_1d(ji)
715         qldif_1d(ji)  = qldif_1d(ji)                                         &
716            + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) & 
717            * rdt_ice                               &
718            + ( 1.0 - zihic ) * fscbq_1d(ji) * rdt_ice
719      END DO  ! ji
720
721      !-------------------------------------------
722      ! Correct temperature, energy and thickness
723      !-------------------------------------------
724      DO ji = kideb, kiut
725         zihgnew        =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 
726         t_su_b(ji)     =  zihgnew * t_su_b(ji) + ( 1.0 - zihgnew ) * rtt
727      END DO  ! ji
728
729      DO jk = 1, nlay_i
730         DO ji = kideb, kiut
731            zihgnew       =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 
732            t_i_b(ji,jk)  =  zihgnew * t_i_b(ji,jk) + ( 1.0 - zihgnew ) * rtt
733            q_i_b(ji,jk)  =  zihgnew * q_i_b(ji,jk)
734         END DO
735      END DO  ! ji
736
737      DO ji = kideb, kiut
738         ht_i_b(ji) = zhgnew(ji)
739      END DO  ! ji
740      !
741      !------------------------------------------------------------------------------|
742      !  6) Snow-Ice formation                                                       |
743      !------------------------------------------------------------------------------|
744      !
745      ! When snow load excesses Archimede's limit, snow-ice interface goes down
746      ! under sea-level, flooding of seawater transforms snow into ice
747      ! dh_snowice is positive for the ice
748      DO ji = kideb, kiut
749
750         dh_snowice(ji) = MAX(zzero,(rhosn*ht_s_b(ji)+(rhoic-rau0) &
751            * ht_i_b(ji))/(rhosn+rau0-rhoic))
752         zhgnew(ji)     = MAX(zhgnew(ji),zhgnew(ji)+dh_snowice(ji))
753         zhnnew            = MIN(ht_s_b(ji),ht_s_b(ji)-dh_snowice(ji))
754
755         !  Changes in ice volume and ice mass.
756         dvnbq_1d(ji)      = a_i_b(ji) * (zhgnew(ji)-ht_i_b(ji))
757         dmgwi_1d(ji)      = dmgwi_1d(ji) + a_i_b(ji) &
758            *(ht_s_b(ji)-zhnnew)*rhosn
759
760         rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) & 
761            * ( zhgnew(ji) - ht_i_b(ji) )*rhoic 
762         rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) & 
763            * ( zhnnew       - ht_s_b(ji) )*rhosn
764
765         !        Equivalent salt flux (1) Snow-ice formation component
766         !        -----------------------------------------------------
767         zji                 = MOD( npb(ji) - 1, jpi ) + 1
768         zjj                 = ( npb(ji) - 1 ) / jpi + 1
769
770         zsm_snowice  = ( rhoic - rhosn ) / rhoic *            &
771            sss_m(zji,zjj) 
772
773         IF ( num_sal .NE. 2 ) zsm_snowice = sm_i_b(ji)
774
775         IF ( num_sal .NE. 4 ) &
776            fseqv_1d(ji)   = fseqv_1d(ji)   + &
777            ( sss_m(zji,zjj) - zsm_snowice ) * &
778            a_i_b(ji)   * &
779            ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice
780         ! new lines
781         IF ( num_sal .EQ. 4 ) &
782            fseqv_1d(ji)   = fseqv_1d(ji)   + &
783            ( sss_m(zji,zjj) - bulk_sal    ) * &
784            a_i_b(ji)   * &
785            ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice
786
787         ! entrapment during snow ice formation
788         i_ice_switch = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i_b(ji) + 1.0e-6 ) )
789         isnowic      = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - dh_snowice(ji) ) ) * &
790            i_ice_switch
791         IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
792            dsm_i_si_1d(ji)  = ( zsm_snowice*dh_snowice(ji) &
793            + sm_i_b(ji) * ht_i_b(ji)                          & 
794            / MAX( ht_i_b(ji) + dh_snowice(ji), zeps)          &
795            - sm_i_b(ji) ) * isnowic     
796
797         !  Actualize new snow and ice thickness.
798         ht_s_b(ji)  = zhnnew
799         ht_i_b(ji)  = zhgnew(ji)
800
801         ! Total ablation ! new lines added to debug
802         IF( ht_i_b(ji).LE.0.0 ) a_i_b(ji) = 0.0
803
804         ! diagnostic ( snow ice growth )
805         zji                 = MOD( npb(ji) - 1, jpi ) + 1
806         zjj                 = ( npb(ji) - 1 ) / jpi + 1
807         diag_sni_gr(zji,zjj)  = diag_sni_gr(zji,zjj) + dh_snowice(ji)*a_i_b(ji) / &
808            rdt_ice
809
810      END DO !ji
811
812   END SUBROUTINE lim_thd_dh
813#else
814   !!======================================================================
815   !!                       ***  MODULE limthd_dh   ***
816   !!                              no sea ice model 
817   !!======================================================================
818CONTAINS
819   SUBROUTINE lim_thd_dh          ! Empty routine
820   END SUBROUTINE lim_thd_dh
821#endif
822END MODULE limthd_dh
Note: See TracBrowser for help on using the repository browser.