source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 3808

Last change on this file since 3808 was 3808, checked in by gm, 9 years ago

dev_MERGE_2012: #997 and #898: in LIM3, no use of qla_ice in coupled mode and its capping in forced mode

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