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 branches/dev_003_CPL/NEMO/LIM_SRC_3 – NEMO

source: branches/dev_003_CPL/NEMO/LIM_SRC_3/limthd_dh.F90 @ 990

Last change on this file since 990 was 990, checked in by smasson, 16 years ago

dev_003_CPL: to the trunk , see ticket #155

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-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      zfsalt_melt(:)  = 0.0
180      ftotal_fin(:)   = 0.0
181      zfdt_init(:)    = 0.0
182      zfdt_final(:)   = 0.0
183
184      DO ji = kideb, kiut
185         old_ht_i_b(ji) = ht_i_b(ji)
186         old_ht_s_b(ji) = ht_s_b(ji)
187      END DO
188      !
189      !------------------------------------------------------------------------------!
190      !  1) Calculate available heat for surface ablation                            !
191      !------------------------------------------------------------------------------!
192      !
193      DO ji = kideb,kiut
194         isnow         = INT( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ) )
195         ztfs(ji)      = isnow * rtt + ( 1.0 - isnow ) * rtt
196         z_f_surf(ji)  = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) *                 & 
197            qsr_ice_1d(ji) - fc_su(ji)
198         z_f_surf(ji)  = MAX( zzero , z_f_surf(ji) ) *                        &
199            MAX( zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) )
200         zfdt_init(ji) = ( z_f_surf(ji) + &
201            MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) )  &
202            * rdt_ice
203      END DO ! ji
204
205      zqfont_su(:) = 0.0
206      zqfont_bo(:) = 0.0
207      dsm_i_se_1d(:) = 0.0     
208      dsm_i_si_1d(:) = 0.0     
209      !
210      !------------------------------------------------------------------------------!
211      !  2) Computing layer thicknesses and  snow and sea-ice enthalpies.            !
212      !------------------------------------------------------------------------------!
213      !
214      ! Layer thickness
215      DO ji = kideb,kiut
216         zh_i(ji) = ht_i_b(ji) / nlay_i
217         zh_s(ji) = ht_s_b(ji) / nlay_s
218      END DO
219
220      ! Total enthalpy of the snow
221      zqt_s(:) = 0.0
222      DO jk = 1, nlay_s
223         DO ji = kideb,kiut
224            zqt_s(ji)    =  zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s
225         END DO
226      END DO
227
228      ! Total enthalpy of the ice
229      zqt_i(:) = 0.0
230      DO jk = 1, nlay_i
231         DO ji = kideb,kiut
232            zqt_i(ji)        =  zqt_i(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
233            zqt_i_lay(ji,jk) =  q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
234         END DO
235      END DO
236      !
237      !------------------------------------------------------------------------------|
238      !  3) Surface ablation and sublimation                                         |
239      !------------------------------------------------------------------------------|
240      !
241      !-------------------------
242      ! 3.1 Snow precips / melt
243      !-------------------------
244      ! Snow accumulation in one thermodynamic time step
245      ! snowfall is partitionned between leads and ice
246      ! if snow fall was uniform, a fraction (1-at_i) would fall into leads
247      ! but because of the winds, more snow falls on leads than on sea ice
248      ! and a greater fraction (1-at_i)^beta of the total mass of snow
249      ! (beta < 1) falls in leads.
250      ! In reality, beta depends on wind speed,
251      ! and should decrease with increasing wind speed but here, it is
252      ! considered as a constant. an average value is 0.66
253      ! Martin Vancoppenolle, December 2006
254
255      ! Snow fall
256      DO ji = kideb, kiut
257         zcoeff = ( 1.0 - ( 1.0 - at_i_b(ji) )**betas ) / at_i_b(ji) 
258         zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn
259      END DO
260      zdh_s_mel(:) =  0.0
261
262      ! Melt of fallen snow
263      DO ji = kideb, kiut
264         ! tatm_ice is now in K
265         zqprec(ji)     =  rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus ) 
266         zqfont_su(ji)  =  z_f_surf(ji) * rdt_ice
267         zdeltah(ji,1)  =  MIN( 0.0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) )
268         zqfont_su(ji)  =  MAX( 0.0 , - zdh_s_pre(ji) - zdeltah(ji,1) )      * &
269            zqprec(ji)
270         zdeltah(ji,1)  =  MAX( - zdh_s_pre(ji) , zdeltah(ji,1) )
271         zdh_s_mel(ji)  =  zdh_s_mel(ji) + zdeltah(ji,1)
272         ! heat conservation
273         qt_s_in(ji,jl) =  qt_s_in(ji,jl) + zqprec(ji) * zdh_s_pre(ji)
274         zqt_s(ji)      =  zqt_s(ji) + zqprec(ji) * zdh_s_pre(ji)
275      END DO
276
277      ! Update total snow heat content
278      zqt_s(ji)         =  MAX ( zqt_s(ji) - zqfont_su(ji) , 0.0 ) 
279      IF( lk_mpp ) CALL mpp_max(zqt_s(ji), kcom = ncomm_ice ) 
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.