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 tags/nemo_v3_2_2/NEMO/LIM_SRC_3 – NEMO

source: tags/nemo_v3_2_2/NEMO/LIM_SRC_3/limthd_dh.F90 @ 2720

Last change on this file since 2720 was 2477, checked in by cetlod, 13 years ago

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