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/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 3517

Last change on this file since 3517 was 3517, checked in by gm, 10 years ago

gm: Branch: dev_r3385_NOCS04_HAMF; #665. update sbccpl ; change LIM3 from equivalent salt flux to salt flux and mass flux

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