New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limthd_dh.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 11 years ago

Merge of 3.4beta into the trunk

  • Property svn:keywords set to Id
File size: 35.4 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   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
10   !!----------------------------------------------------------------------
11#if defined key_lim3
12   !!----------------------------------------------------------------------
13   !!   'key_lim3'                                      LIM3 sea-ice model
14   !!----------------------------------------------------------------------
15   !!   lim_thd_dh  : vertical accr./abl. and lateral ablation of sea ice
16   !!----------------------------------------------------------------------
17   USE par_oce          ! ocean parameters
18   USE phycst           ! physical constants (OCE directory)
19   USE sbc_oce          ! Surface boundary condition: ocean fields
20   USE ice              ! LIM variables
21   USE par_ice          ! LIM parameters
22   USE thd_ice          ! LIM thermodynamics
23   USE in_out_manager   ! I/O manager
24   USE lib_mpp          ! MPP library
25   USE wrk_nemo         ! work arrays
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   lim_thd_dh   ! called by lim_thd
31
32   REAL(wp) ::   epsi20 = 1e-20   ! constant values
33   REAL(wp) ::   epsi13 = 1e-13   !
34   REAL(wp) ::   epsi16 = 1e-16   !
35   REAL(wp) ::   zzero  = 0.e0    !
36   REAL(wp) ::   zone   = 1.e0    !
37
38   !!----------------------------------------------------------------------
39   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
40   !! $Id$
41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE lim_thd_dh( kideb, kiut, jl )
46      !!------------------------------------------------------------------
47      !!                ***  ROUTINE lim_thd_dh  ***
48      !!
49      !! ** Purpose :   determines variations of ice and snow thicknesses.
50      !!
51      !! ** Method  :   Ice/Snow surface melting arises from imbalance in surface fluxes
52      !!              Bottom accretion/ablation arises from flux budget
53      !!              Snow thickness can increase by precipitation and decrease by sublimation
54      !!              If snow load excesses Archmiede limit, snow-ice is formed by
55      !!              the flooding of sea-water in the snow
56      !!
57      !!                 1) Compute available flux of heat for surface ablation
58      !!                 2) Compute snow and sea ice enthalpies
59      !!                 3) Surface ablation and sublimation
60      !!                 4) Bottom accretion/ablation
61      !!                 5) Case of Total ablation
62      !!                 6) Snow ice formation
63      !!
64      !! References : Bitz and Lipscomb, 1999, J. Geophys. Res.
65      !!              Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
66      !!              Vancoppenolle, Fichefet and Bitz, 2005, Geophys. Res. Let.
67      !!              Vancoppenolle et al.,2009, Ocean Modelling
68      !!------------------------------------------------------------------
69      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
70      INTEGER , INTENT(in) ::   jl            ! Thickness cateogry number
71      !!
72      INTEGER  ::   ji , jk        ! dummy loop indices
73      INTEGER  ::   zji, zjj       ! 2D corresponding indices to ji
74      INTEGER  ::   isnow          ! switch for presence (1) or absence (0) of snow
75      INTEGER  ::   isnowic        ! snow ice formation not
76      INTEGER  ::   i_ice_switch   ! ice thickness above a certain treshold or not
77      INTEGER  ::   iter
78
79      REAL(wp) ::   zzfmass_i, zihgnew                     ! local scalar
80      REAL(wp) ::   zzfmass_s, zhsnew, ztmelts             ! local scalar
81      REAL(wp) ::   zhn, zdhcf, zdhbf, zhni, zhnfi, zihg   !
82      REAL(wp) ::   zdhnm, zhnnew, zhisn, zihic, zzc       !
83      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
84      REAL(wp) ::   zds          ! increment of bottom ice salinity
85      REAL(wp) ::   zcoeff       ! dummy argument for snowfall partitioning over ice and leads
86      REAL(wp) ::   zsm_snowice  ! snow-ice salinity
87      REAL(wp) ::   zswi1        ! switch for computation of bottom salinity
88      REAL(wp) ::   zswi12       ! switch for computation of bottom salinity
89      REAL(wp) ::   zswi2        ! switch for computation of bottom salinity
90      REAL(wp) ::   zgrr         ! bottom growth rate
91      REAL(wp) ::   ztform       ! bottom formation temperature
92      !
93      REAL(wp), POINTER, DIMENSION(:) ::   zh_i        ! ice layer thickness
94      REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness
95      REAL(wp), POINTER, DIMENSION(:) ::   ztfs        ! melting point
96      REAL(wp), POINTER, DIMENSION(:) ::   zhsold      ! old snow thickness
97      REAL(wp), POINTER, DIMENSION(:) ::   zqprec      ! energy of fallen snow
98      REAL(wp), POINTER, DIMENSION(:) ::   zqfont_su   ! incoming, remaining surface energy
99      REAL(wp), POINTER, DIMENSION(:) ::   zqfont_bo   ! incoming, bottom energy
100      REAL(wp), POINTER, DIMENSION(:) ::   z_f_surf    ! surface heat for ablation
101      REAL(wp), POINTER, DIMENSION(:) ::   zhgnew      ! new ice thickness
102      REAL(wp), POINTER, DIMENSION(:) ::   zfmass_i    !
103
104      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel     ! snow melt
105      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_pre     ! snow precipitation
106      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_sub     ! snow sublimation
107      REAL(wp), POINTER, DIMENSION(:) ::   zfsalt_melt   ! salt flux due to ice melt
108
109      REAL(wp), POINTER, DIMENSION(:,:) ::   zdeltah
110
111      ! Pathological cases
112      REAL(wp), POINTER, DIMENSION(:) ::   zfdt_init   ! total incoming heat for ice melt
113      REAL(wp), POINTER, DIMENSION(:) ::   zfdt_final  ! total remaing heat for ice melt
114      REAL(wp), POINTER, DIMENSION(:) ::   zqt_i       ! total ice heat content
115      REAL(wp), POINTER, DIMENSION(:) ::   zqt_s       ! total snow heat content
116      REAL(wp), POINTER, DIMENSION(:) ::   zqt_dummy   ! dummy heat content
117
118      REAL(wp), POINTER, DIMENSION(:,:) ::   zqt_i_lay   ! total ice heat content
119
120      ! Heat conservation
121      INTEGER  ::   num_iter_max, numce_dh
122      REAL(wp) ::   meance_dh
123      REAL(wp), POINTER, DIMENSION(:) ::   zinnermelt
124      REAL(wp), POINTER, DIMENSION(:) ::   zfbase, zdq_i
125      !!------------------------------------------------------------------
126
127      CALL wrk_alloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i )
128      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 )
129      CALL wrk_alloc( jpij, zinnermelt, zfbase, zdq_i )
130      CALL wrk_alloc( jpij, jkmax, zdeltah, zqt_i_lay )
131
132      zfsalt_melt(:)  = 0._wp
133      ftotal_fin(:)   = 0._wp
134      zfdt_init(:)    = 0._wp
135      zfdt_final(:)   = 0._wp
136
137      DO ji = kideb, kiut
138         old_ht_i_b(ji) = ht_i_b(ji)
139         old_ht_s_b(ji) = ht_s_b(ji)
140      END DO
141      !
142      !------------------------------------------------------------------------------!
143      !  1) Calculate available heat for surface ablation                            !
144      !------------------------------------------------------------------------------!
145      !
146      DO ji = kideb, kiut
147         isnow         = INT( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ) )
148         ztfs(ji)      = isnow * rtt + ( 1.0 - isnow ) * rtt
149         z_f_surf(ji)  = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)
150         z_f_surf(ji)  = MAX( zzero , z_f_surf(ji) ) * MAX( zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) )
151         zfdt_init(ji) = ( z_f_surf(ji) + MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) ) * rdt_ice
152      END DO ! ji
153
154      zqfont_su  (:) = 0._wp
155      zqfont_bo  (:) = 0._wp
156      dsm_i_se_1d(:) = 0._wp     
157      dsm_i_si_1d(:) = 0._wp   
158      !
159      !------------------------------------------------------------------------------!
160      !  2) Computing layer thicknesses and  snow and sea-ice enthalpies.            !
161      !------------------------------------------------------------------------------!
162      !
163      DO ji = kideb, kiut     ! Layer thickness
164         zh_i(ji) = ht_i_b(ji) / nlay_i
165         zh_s(ji) = ht_s_b(ji) / nlay_s
166      END DO
167      !
168      zqt_s(:) = 0._wp        ! Total enthalpy of the snow
169      DO jk = 1, nlay_s
170         DO ji = kideb, kiut
171            zqt_s(ji) =  zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s
172         END DO
173      END DO
174      !
175      zqt_i(:) = 0._wp        ! Total enthalpy of the ice
176      DO jk = 1, nlay_i
177         DO ji = kideb, kiut
178            zzc = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
179            zqt_i(ji)        =  zqt_i(ji) + zzc
180            zqt_i_lay(ji,jk) =              zzc
181         END DO
182      END DO
183      !
184      !------------------------------------------------------------------------------|
185      !  3) Surface ablation and sublimation                                         |
186      !------------------------------------------------------------------------------|
187      !
188      !-------------------------
189      ! 3.1 Snow precips / melt
190      !-------------------------
191      ! Snow accumulation in one thermodynamic time step
192      ! snowfall is partitionned between leads and ice
193      ! if snow fall was uniform, a fraction (1-at_i) would fall into leads
194      ! but because of the winds, more snow falls on leads than on sea ice
195      ! and a greater fraction (1-at_i)^beta of the total mass of snow
196      ! (beta < 1) falls in leads.
197      ! In reality, beta depends on wind speed,
198      ! and should decrease with increasing wind speed but here, it is
199      ! considered as a constant. an average value is 0.66
200      ! Martin Vancoppenolle, December 2006
201
202      ! Snow fall
203      DO ji = kideb, kiut
204         zcoeff = ( 1.0 - ( 1.0 - at_i_b(ji) )**betas ) / at_i_b(ji) 
205         zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn
206      END DO
207      zdh_s_mel(:) =  0._wp
208
209      ! Melt of fallen snow
210      DO ji = kideb, kiut
211         ! tatm_ice is now in K
212         zqprec   (ji)   =  rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus ) 
213         zqfont_su(ji)   =  z_f_surf(ji) * rdt_ice
214         zdeltah  (ji,1) =  MIN( 0.e0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) )
215         zqfont_su(ji)   =  MAX( 0.e0 , - zdh_s_pre(ji) - zdeltah(ji,1)              ) * zqprec(ji)
216         zdeltah  (ji,1) =  MAX( - zdh_s_pre(ji) , zdeltah(ji,1) )
217         zdh_s_mel(ji)   =  zdh_s_mel(ji) + zdeltah(ji,1)
218         ! heat conservation
219         qt_s_in(ji,jl)  =  qt_s_in(ji,jl) + zqprec(ji) * zdh_s_pre(ji)
220         zqt_s  (ji)     =  zqt_s  (ji)    + zqprec(ji) * zdh_s_pre(ji)
221         zqt_s  (ji)     =  MAX( zqt_s(ji) - zqfont_su(ji) , 0.e0 ) 
222      END DO
223
224
225      ! Snow melt due to surface heat imbalance
226      DO jk = 1, nlay_s
227         DO ji = kideb, kiut
228            zdeltah  (ji,jk) = - zqfont_su(ji) / q_s_b(ji,jk)
229            zqfont_su(ji)    =  MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * q_s_b(ji,jk) 
230            zdeltah  (ji,jk) =  MAX( zdeltah(ji,jk) , - zh_s(ji) )
231            zdh_s_mel(ji)    =  zdh_s_mel(ji) + zdeltah(ji,jk)        ! resulting melt of snow   
232         END DO
233      END DO
234
235      ! Apply snow melt to snow depth
236      DO ji = kideb, kiut
237         dh_s_tot(ji)   =  zdh_s_mel(ji) + zdh_s_pre(ji)
238         ! Old and new snow depths
239         zhsold(ji)     =  ht_s_b(ji)
240         zhsnew         =  ht_s_b(ji) + dh_s_tot(ji)
241         ! If snow is still present zhn = 1, else zhn = 0
242         zhn            =  1.0 - MAX( zzero , SIGN( zone , - zhsnew ) )
243         ht_s_b(ji)     =  MAX( zzero , zhsnew )
244         ! Volume and mass variations of snow
245         dvsbq_1d  (ji) =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_mel(ji) )
246         dvsbq_1d  (ji) =  MIN( zzero, dvsbq_1d(ji) )
247         rdmsnif_1d(ji) =  rdmsnif_1d(ji) + rhosn * dvsbq_1d(ji)
248      END DO ! ji
249
250      !--------------------------
251      ! 3.2 Surface ice ablation
252      !--------------------------
253      DO ji = kideb, kiut 
254         dh_i_surf(ji) =  0._wp
255         z_f_surf (ji) =  zqfont_su(ji) / rdt_ice ! heat conservation test
256         zdq_i    (ji) =  0._wp
257      END DO ! ji
258
259      DO jk = 1, nlay_i
260         DO ji = kideb, kiut 
261            !                                                    ! melt of layer jk
262            zdeltah  (ji,jk) = - zqfont_su(ji) / q_i_b(ji,jk)
263            !                                                    ! recompute heat available
264            zqfont_su(ji)    = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk) 
265            !                                                    ! melt of layer jk cannot be higher than its thickness
266            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) )
267            !                                                    ! update surface melt
268            dh_i_surf(ji)    = dh_i_surf(ji) + zdeltah(ji,jk) 
269            !                                                    ! for energy conservation
270            zdq_i    (ji)    = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice
271            !
272            ! contribution to ice-ocean salt flux
273            zji = MOD( npb(ji) - 1 , jpi ) + 1
274            zjj =    ( npb(ji) - 1 ) / jpi + 1
275            zfsalt_melt(ji) = zfsalt_melt(ji) + ( sss_m(zji,zjj) - sm_i_b(ji) ) * a_i_b(ji)    &
276               &                              * MIN( zdeltah(ji,jk) , 0.e0 ) * rhoic / rdt_ice 
277         END DO
278      END DO
279
280      !                     !-------------------
281      IF( con_i ) THEN      ! Conservation test
282         !                  !-------------------
283         numce_dh  = 0
284         meance_dh = 0._wp
285         DO ji = kideb, kiut
286            IF ( ( z_f_surf(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN
287               numce_dh  = numce_dh + 1
288               meance_dh = meance_dh + z_f_surf(ji) + zdq_i(ji)
289            ENDIF
290            IF( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN!
291               WRITE(numout,*) ' ALERTE heat loss for surface melt '
292               WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl
293               WRITE(numout,*) ' ht_i_b       : ', ht_i_b(ji)
294               WRITE(numout,*) ' z_f_surf     : ', z_f_surf(ji)
295               WRITE(numout,*) ' zdq_i        : ', zdq_i(ji)
296               WRITE(numout,*) ' ht_i_b       : ', ht_i_b(ji)
297               WRITE(numout,*) ' fc_bo_i      : ', fc_bo_i(ji)
298               WRITE(numout,*) ' fbif_1d      : ', fbif_1d(ji)
299               WRITE(numout,*) ' qlbbq_1d     : ', qlbbq_1d(ji)
300               WRITE(numout,*) ' s_i_new      : ', s_i_new(ji)
301               WRITE(numout,*) ' sss_m        : ', sss_m(zji,zjj)
302            ENDIF
303         END DO
304         !
305         IF( numce_dh > 0 )   meance_dh = meance_dh / numce_dh
306         WRITE(numout,*) ' Error report - Category : ', jl
307         WRITE(numout,*) ' ~~~~~~~~~~~~ '
308         WRITE(numout,*) ' Number of points where there is sur. me. error : ', numce_dh
309         WRITE(numout,*) ' Mean basal growth error on error points : ', meance_dh
310         !
311      ENDIF
312
313      !----------------------
314      ! 3.3 Snow sublimation
315      !----------------------
316
317      DO ji = kideb, kiut
318         ! if qla is positive (upwards), heat goes to the atmosphere, therefore
319         ! snow sublimates, if qla is negative (downwards), snow condensates
320         zdh_s_sub(ji)    =  - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice
321         dh_s_tot (ji)    =  dh_s_tot(ji) + zdh_s_sub(ji)
322         zdhcf            =  ht_s_b(ji) + zdh_s_sub(ji) 
323         ht_s_b   (ji)    =  MAX( zzero , zdhcf )
324         ! we recompute dh_s_tot
325         dh_s_tot (ji)    =  ht_s_b(ji) - zhsold(ji)
326         qt_s_in  (ji,jl) =  qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1)
327      END DO
328
329      zqt_dummy(:) = 0.e0
330      DO jk = 1, nlay_s
331         DO ji = kideb,kiut
332            q_s_b    (ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
333            zqt_dummy(ji)    =  zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s            ! heat conservation
334         END DO
335      END DO
336
337      DO jk = 1, nlay_s
338         DO ji = kideb, kiut
339            ! In case of disparition of the snow, we have to update the snow temperatures
340            zhisn  =  MAX( zzero , SIGN( zone, - ht_s_b(ji) ) )
341            t_s_b(ji,jk) = ( 1.0 - zhisn ) * t_s_b(ji,jk) + zhisn * rtt
342            q_s_b(ji,jk) = ( 1.0 - zhisn ) * q_s_b(ji,jk)
343         END DO
344      END DO
345
346      !
347      !------------------------------------------------------------------------------!
348      ! 4) Basal growth / melt                                                       !
349      !------------------------------------------------------------------------------!
350      !
351      ! Ice basal growth / melt is given by the ratio of heat budget over basal
352      ! ice heat content.  Basal heat budget is given by the difference between
353      ! the inner conductive flux  (fc_bo_i), from the open water heat flux
354      ! (qlbbqb) and the turbulent ocean flux (fbif).
355      ! fc_bo_i is positive downwards. fbif and qlbbq are positive to the ice
356
357      !-----------------------------------------------------
358      ! 4.1 Basal growth - (a) salinity not varying in time
359      !-----------------------------------------------------
360      IF(  num_sal /= 2  .AND.  num_sal /= 4  ) THEN
361         DO ji = kideb, kiut
362            IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0.0  ) THEN
363               s_i_new(ji)         =  sm_i_b(ji)
364               ! Melting point in K
365               ztmelts             =  - tmut * s_i_new(ji) + rtt 
366               ! New ice heat content (Bitz and Lipscomb, 1999)
367               ztform              =  t_i_b(ji,nlay_i)  ! t_bo_b crashes in the
368               ! Baltic
369               q_i_b(ji,nlay_i+1)  = rhoic * (  cpic * ( ztmelts - ztform )                                &
370                  &                           + lfus * (  1.0 - ( ztmelts - rtt ) / ( ztform - rtt )  )    &
371                  &                           - rcp  * ( ztmelts - rtt )                                 )
372               ! Basal growth rate = - F*dt / q
373               dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
374            ENDIF
375         END DO
376      ENDIF
377
378      !-------------------------------------------------
379      ! 4.1 Basal growth - (b) salinity varying in time
380      !-------------------------------------------------
381      IF(  num_sal == 2 .OR.  num_sal == 4  ) THEN
382         ! the growth rate (dh_i_bott) is function of the new ice
383         ! heat content (q_i_b(nlay_i+1)). q_i_b depends on the new ice
384         ! salinity (snewice). snewice depends on dh_i_bott
385         ! it converges quickly, so, no problem
386         ! See Vancoppenolle et al., OM08 for more info on this
387
388         ! Initial value (tested 1D, can be anything between 1 and 20)
389         num_iter_max = 4
390         s_i_new(:)   = 4.0
391
392         ! Iterative procedure
393         DO iter = 1, num_iter_max
394            DO ji = kideb, kiut
395               IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0  ) THEN
396                  zji = MOD( npb(ji) - 1, jpi ) + 1
397                  zjj = ( npb(ji) - 1 ) / jpi + 1
398                  ! Melting point in K
399                  ztmelts             =   - tmut * s_i_new(ji) + rtt 
400                  ! New ice heat content (Bitz and Lipscomb, 1999)
401                  q_i_b(ji,nlay_i+1)  =  rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             &
402                     &                            + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )   &
403                     &                            - rcp * ( ztmelts-rtt )                                     )
404                  ! Bottom growth rate = - F*dt / q
405                  dh_i_bott(ji) =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)
406                  ! New ice salinity ( Cox and Weeks, JGR, 1988 )
407                  ! zswi2  (1) if dh_i_bott/rdt .GT. 3.6e-7
408                  ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8
409                  ! zswi1  (1) if dh_i_bott/rdt .LT. 2.0e-8
410                  zgrr   = MIN( 1.0e-3, MAX ( dh_i_bott(ji) / rdt_ice , epsi13 ) )
411                  zswi2  = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) ) 
412                  zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
413                  zswi1  = 1. - zswi2 * zswi12 
414                  zfracs = zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
415                     &                   + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 
416                  zds         = zfracs * sss_m(zji,zjj) - s_i_new(ji)
417                  s_i_new(ji) = zfracs * sss_m(zji,zjj)
418               ENDIF ! fc_bo_i
419            END DO ! ji
420         END DO ! iter
421
422         ! Final values
423         DO ji = kideb, kiut
424            IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN
425               ! New ice salinity must not exceed 15 psu
426               s_i_new(ji) = MIN( s_i_new(ji), s_i_max )
427               ! Metling point in K
428               ztmelts     =   - tmut * s_i_new(ji) + rtt 
429               ! New ice heat content (Bitz and Lipscomb, 1999)
430               q_i_b(ji,nlay_i+1)  =  rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                              &
431                  &                            + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )    &
432                  &                            - rcp * ( ztmelts - rtt )                                    )
433               ! Basal growth rate = - F*dt / q
434               dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
435               ! Salinity update
436               ! entrapment during bottom growth
437               dsm_i_se_1d(ji) = ( s_i_new(ji) * dh_i_bott(ji) + sm_i_b(ji) * ht_i_b(ji) )    &
438                  &            / MAX( ht_i_b(ji) + dh_i_bott(ji) ,epsi13 ) - sm_i_b(ji)
439            ENDIF ! heat budget
440         END DO
441      ENDIF
442
443      !----------------
444      ! 4.2 Basal melt
445      !----------------
446      meance_dh = 0._wp
447      numce_dh  = 0
448      zinnermelt(:) = 0._wp
449
450      DO ji = kideb, kiut
451         ! heat convergence at the surface > 0
452         IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0._wp  ) THEN
453            s_i_new(ji)   =  s_i_b(ji,nlay_i)
454            zqfont_bo(ji) =  rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) )
455            zfbase(ji)    =  zqfont_bo(ji) / rdt_ice     ! heat conservation test
456            zdq_i(ji)     =  0._wp
457            dh_i_bott(ji) =  0._wp
458         ENDIF
459      END DO
460
461      DO jk = nlay_i, 1, -1
462         DO ji = kideb, kiut
463            IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN
464               ztmelts            =   - tmut * s_i_b(ji,jk) + rtt 
465               IF( t_i_b(ji,jk) >= ztmelts ) THEN
466                  zdeltah(ji,jk)  = - zh_i(ji)
467                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)
468                  zinnermelt(ji)   = 1._wp
469               ELSE  ! normal ablation
470                  zdeltah(ji,jk)  = - zqfont_bo(ji) / q_i_b(ji,jk)
471                  zqfont_bo(ji)   = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk)
472                  zdeltah(ji,jk)  = MAX(zdeltah(ji,jk), - zh_i(ji) )
473                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)
474                  zdq_i(ji)       = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice
475                  ! contribution to salt flux
476                  zji             = MOD( npb(ji) - 1, jpi ) + 1
477                  zjj             = ( npb(ji) - 1 ) / jpi + 1
478                  zfsalt_melt(ji) = zfsalt_melt(ji) + ( sss_m(zji,zjj) - sm_i_b(ji)   ) * a_i_b(ji)   &
479                     &                              * MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 
480               ENDIF
481            ENDIF
482         END DO ! ji
483      END DO ! jk
484
485      !                     !-------------------
486      IF( con_i ) THEN      ! Conservation test
487      !                     !-------------------
488         DO ji = kideb, kiut
489            IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0.e0  ) THEN
490               IF( ( zfbase(ji) + zdq_i(ji) ) >= 1.e-3 ) THEN
491                  numce_dh  = numce_dh + 1
492                  meance_dh = meance_dh + zfbase(ji) + zdq_i(ji)
493               ENDIF
494               IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN
495                  WRITE(numout,*) ' ALERTE heat loss for basal melt : 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 : ', INT( zinnermelt(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._wp )   a_i_b(ji) = 0._wp
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      CALL wrk_dealloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i )
699      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 )
700      CALL wrk_dealloc( jpij, zinnermelt, zfbase, zdq_i )
701      CALL wrk_dealloc( jpij, jkmax, zdeltah, zqt_i_lay )
702      !
703   END SUBROUTINE lim_thd_dh
704   
705#else
706   !!----------------------------------------------------------------------
707   !!   Default option                               NO  LIM3 sea-ice model
708   !!----------------------------------------------------------------------
709CONTAINS
710   SUBROUTINE lim_thd_dh          ! Empty routine
711   END SUBROUTINE lim_thd_dh
712#endif
713
714   !!======================================================================
715END MODULE limthd_dh
Note: See TracBrowser for help on using the repository browser.