source: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 4506

Last change on this file since 4506 was 4506, checked in by vancop, 8 years ago

[Heat conservation in LIM3, part HC1 (LIM_SRC_3_HC17)]

  • Property svn:keywords set to Id
File size: 38.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 rdm_snw & rdm_ice
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 = 1.e-20   ! constant values
35   REAL(wp) ::   epsi10 = 1.e-10   !
36   REAL(wp) ::   epsi13 = 1.e-13   !
37   REAL(wp) ::   zzero  = 0._wp    !
38   REAL(wp) ::   zone   = 1._wp    !
39
40   !!----------------------------------------------------------------------
41   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
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  ::   i_use          ! debugging flag
77      INTEGER  ::   isnow          ! switch for presence (1) or absence (0) of snow
78      INTEGER  ::   isnowic        ! snow ice formation not
79      INTEGER  ::   i_ice_switch   ! ice thickness above a certain treshold or not
80      INTEGER  ::   iter
81
82      REAL(wp) ::   zzfmass_i, zihgnew                     ! local scalar
83      REAL(wp) ::   zzfmass_s, zhsnew, ztmelts             ! local scalar
84      REAL(wp) ::   zhn, zdhcf, zdhbf, zhni, zhnfi, zihg   !
85      REAL(wp) ::   zdhnm, zhnnew, zhisn, zihic, zzc       !
86      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
87      REAL(wp) ::   zcoeff       ! dummy argument for snowfall partitioning over ice and leads
88      REAL(wp) ::   zs_snic  ! 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) ::   zt_i_new     ! bottom formation temperature
94
95      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
96      REAL(wp) ::   zEi          ! specific enthalpy of sea ice (J/kg)
97      REAL(wp) ::   zEw          ! specific enthalpy of exchanged water (J/kg)
98      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg)
99      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean
100      REAL(wp) ::   zsstK        ! SST in Kelvin
101
102      REAL(wp), POINTER, DIMENSION(:) ::   zh_i        ! ice layer thickness
103      REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness
104      REAL(wp), POINTER, DIMENSION(:) ::   ztfs        ! melting point
105      REAL(wp), POINTER, DIMENSION(:) ::   zhsold      ! old snow thickness
106      REAL(wp), POINTER, DIMENSION(:) ::   zqprec      ! energy of fallen snow
107      REAL(wp), POINTER, DIMENSION(:) ::   zqfont_su   ! incoming, remaining surface energy
108      REAL(wp), POINTER, DIMENSION(:) ::   zqfont_bo   ! incoming, bottom energy
109      REAL(wp), POINTER, DIMENSION(:) ::   z_f_surf    ! surface heat for ablation
110      REAL(wp), POINTER, DIMENSION(:) ::   zhgnew      ! new ice thickness
111      REAL(wp), POINTER, DIMENSION(:) ::   zfmass_i    !
112
113      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel   ! snow melt
114      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_pre   ! snow precipitation
115      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_sub   ! snow sublimation
116
117      REAL(wp), POINTER, DIMENSION(:,:) ::   zdeltah
118
119      ! Pathological cases
120      REAL(wp), POINTER, DIMENSION(:) ::   zfdt_init   ! total incoming heat for ice melt
121      REAL(wp), POINTER, DIMENSION(:) ::   zfdt_final  ! total remaing heat for ice melt
122      REAL(wp), POINTER, DIMENSION(:) ::   zqt_i       ! total ice heat content
123      REAL(wp), POINTER, DIMENSION(:) ::   zqt_s       ! total snow heat content
124      REAL(wp), POINTER, DIMENSION(:) ::   zqt_dummy   ! dummy heat content
125
126      REAL(wp), POINTER, DIMENSION(:,:) ::   zqt_i_lay   ! total ice heat content
127
128      ! mass and salt flux (clem)
129      REAL(wp) :: zdvres, zdvsur, zdvbot, zswitch_sal
130      REAL(wp), POINTER, DIMENSION(:) ::   zviold, zvsold   ! old ice volume...
131
132      ! Heat conservation
133      INTEGER  ::   num_iter_max, numce_dh
134      REAL(wp) ::   meance_dh
135      REAL(wp) ::   zinda 
136      REAL(wp), POINTER, DIMENSION(:) ::   zinnermelt
137      REAL(wp), POINTER, DIMENSION(:) ::   zfbase, zdq_i
138      !!------------------------------------------------------------------
139
140      ! Discriminate between varying salinity (num_sal=2) and prescribed cases (other values)
141      SELECT CASE( num_sal )                       ! varying salinity or not
142         CASE( 1, 3, 4 ) ;   zswitch_sal = 0       ! prescribed salinity profile
143         CASE( 2 )       ;   zswitch_sal = 1       ! varying salinity profile
144      END SELECT
145
146      CALL wrk_alloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i )
147      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy )
148      CALL wrk_alloc( jpij, zinnermelt, zfbase, zdq_i )
149      CALL wrk_alloc( jpij, jkmax, zdeltah, zqt_i_lay )
150
151      CALL wrk_alloc( jpij, zviold, zvsold ) ! clem
152     
153      ftotal_fin(:) = 0._wp
154      zfdt_init (:) = 0._wp
155      zfdt_final(:) = 0._wp
156
157      dh_i_surf (:) = 0._wp
158      dh_i_bott (:) = 0._wp
159      dh_snowice(:) = 0._wp
160
161      DO ji = kideb, kiut
162         old_ht_i_b(ji) = ht_i_b(ji)
163         old_ht_s_b(ji) = ht_s_b(ji)
164         zviold(ji) = a_i_b(ji) * ht_i_b(ji) ! clem
165         zvsold(ji) = a_i_b(ji) * ht_s_b(ji) ! clem
166      END DO
167      !
168      !------------------------------------------------------------------------------!
169      !  1) Calculate available heat for surface ablation                            !
170      !------------------------------------------------------------------------------!
171      !
172      DO ji = kideb, kiut
173         isnow         = INT(  1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_s_b(ji) )  )  )
174         ztfs     (ji) = isnow * rtt + ( 1.0 - isnow ) * rtt
175         z_f_surf (ji) = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)
176         z_f_surf (ji) = MAX(  zzero , z_f_surf(ji)  ) * MAX(  zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) )  )
177         zfdt_init(ji) = ( z_f_surf(ji) + MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) ) * rdt_ice
178      END DO ! ji
179
180      zqfont_su  (:) = 0._wp
181      zqfont_bo  (:) = 0._wp
182      dsm_i_se_1d(:) = 0._wp     
183      dsm_i_si_1d(:) = 0._wp   
184      !
185      !------------------------------------------------------------------------------!
186      !  2) Computing layer thicknesses and  snow and sea-ice enthalpies.            !
187      !------------------------------------------------------------------------------!
188      !
189      DO ji = kideb, kiut     ! Layer thickness
190         zh_i(ji) = ht_i_b(ji) / REAL( nlay_i )
191         zh_s(ji) = ht_s_b(ji) / REAL( nlay_s )
192      END DO
193      !
194      zqt_s(:) = 0._wp        ! Total enthalpy of the snow
195      DO jk = 1, nlay_s
196         DO ji = kideb, kiut
197            zqt_s(ji) =  zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s )
198         END DO
199      END DO
200      !
201      zqt_i(:) = 0._wp        ! Total enthalpy of the ice
202      DO jk = 1, nlay_i
203         DO ji = kideb, kiut
204            zzc = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i )
205            zqt_i(ji)        =  zqt_i(ji) + zzc
206            zqt_i_lay(ji,jk) =              zzc
207         END DO
208      END DO
209      !
210      !------------------------------------------------------------------------------|
211      !  3) Surface ablation and sublimation                                         |
212      !------------------------------------------------------------------------------|
213      !
214      !-------------------------
215      ! 3.1 Snow precips / melt
216      !-------------------------
217      ! Snow accumulation in one thermodynamic time step
218      ! snowfall is partitionned between leads and ice
219      ! if snow fall was uniform, a fraction (1-at_i) would fall into leads
220      ! but because of the winds, more snow falls on leads than on sea ice
221      ! and a greater fraction (1-at_i)^beta of the total mass of snow
222      ! (beta < 1) falls in leads.
223      ! In reality, beta depends on wind speed,
224      ! and should decrease with increasing wind speed but here, it is
225      ! considered as a constant. an average value is 0.66
226      ! Martin Vancoppenolle, December 2006
227
228      ! Snow fall
229      DO ji = kideb, kiut
230         zcoeff = ( 1.0 - ( 1.0 - at_i_b(ji) )**betas ) / at_i_b(ji) 
231         zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn
232      END DO
233      zdh_s_mel(:) =  0._wp
234
235      ! Melt of fallen snow
236      DO ji = kideb, kiut
237         ! tatm_ice is now in K
238         zqprec   (ji)   =  rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus ) 
239         zqfont_su(ji)   =  z_f_surf(ji) * rdt_ice
240         zdeltah  (ji,1) =  MIN( 0.e0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) )
241         zqfont_su(ji)   =  MAX( 0.e0 , - zdh_s_pre(ji) - zdeltah(ji,1)              ) * zqprec(ji)
242         zdeltah  (ji,1) =  MAX( - zdh_s_pre(ji) , zdeltah(ji,1) )
243         zdh_s_mel(ji)   =  zdh_s_mel(ji) + zdeltah(ji,1)
244         ! heat conservation
245         qt_s_in(ji,jl)  =  qt_s_in(ji,jl) + zqprec(ji) * zdh_s_pre(ji)
246         zqt_s  (ji)     =  zqt_s  (ji)    + zqprec(ji) * zdh_s_pre(ji)
247         zqt_s  (ji)     =  MAX( zqt_s(ji) - zqfont_su(ji) , 0.e0 ) 
248      END DO
249
250
251      ! Snow melt due to surface heat imbalance
252      DO jk = 1, nlay_s
253         DO ji = kideb, kiut
254            zdeltah  (ji,jk) = - zqfont_su(ji) / q_s_b(ji,jk)
255            zqfont_su(ji)    =  MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * q_s_b(ji,jk) 
256            zdeltah  (ji,jk) =  MAX( zdeltah(ji,jk) , - zh_s(ji) )
257            zdh_s_mel(ji)    =  zdh_s_mel(ji) + zdeltah(ji,jk)        ! resulting melt of snow   
258         END DO
259      END DO
260
261      ! Apply snow melt to snow depth
262      DO ji = kideb, kiut
263         dh_s_tot(ji)   =  zdh_s_mel(ji) + zdh_s_pre(ji)
264         ! Old and new snow depths
265         zhsold(ji)     =  ht_s_b(ji)
266         zhsnew         =  ht_s_b(ji) + dh_s_tot(ji)
267         ! If snow is still present zhn = 1, else zhn = 0
268         zhn            =  1.0 - MAX(  zzero , SIGN( zone , - zhsnew )  )
269         ht_s_b(ji)     =  MAX( zzero , zhsnew )
270         ! we recompute dh_s_tot (clem)
271         dh_s_tot (ji)  =  ht_s_b(ji) - zhsold(ji)
272         ! Volume and mass variations of snow
273         dvsbq_1d  (ji) =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_pre(ji) )
274         dvsbq_1d  (ji) =  MIN( zzero, dvsbq_1d(ji) )
275         !clem rdm_snw_1d(ji) =  rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji)
276      END DO ! ji
277
278      !--------------------------
279      ! 3.2 Surface ice ablation
280      !--------------------------
281      DO ji = kideb, kiut 
282         z_f_surf (ji) =  zqfont_su(ji) * r1_rdtice   ! heat conservation test
283         zdq_i    (ji) =  0._wp
284      END DO ! ji
285
286      DO jk = 1, nlay_i
287         DO ji = kideb, kiut 
288            zEi            = - q_i_b(ji,jk) / rhoic                ! Specific enthalpy of layer k [J/kg, <0]
289
290            ztmelts        = - tmut * s_i_b(ji,jk) + rtt           ! Melting point of layer k [K]
291
292            zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0]
293
294            zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0
295
296            zfmdt          = - zqfont_su(ji) / zdE                 ! Mass flux to the ocean [kg/m2, >0]
297
298            zdeltah(ji,jk) = - zfmdt / rhoic                       ! Melt of layer jk [m, <0]
299
300            zqfont_su(ji)  = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * rhoic * ( - zdE )
301                                                                   ! Energy remaining in case of melting of the full layer [J/m2, >0]
302            zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0]
303
304            dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt
305
306            zfmdt          = - rhoic*zdeltah(ji,jk)                ! Recompute mass flux [kg/m2, >0]
307
308            zQm            = MAX ( zfmdt, 0._wp ) * zEw            ! Energy of the melt water sent to the ocean [J/m2, <0]
309
310            ! Contribution to salt flux
311            sfx_thd_1d(ji)   = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    &
312               &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice
313            ! Contribution to heat flux
314            rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm * a_i_b(ji)
315           
316            ! Conservation test
317            zdq_i(ji)      = zdq_i(ji) + zdeltah(ji,jk) * rhoic * zdE * r1_rdtice  ! ? still don't know
318         END DO
319      END DO
320
321      !                                          !-------------------
322      IF( con_i .AND. jiindex_1d > 0 ) THEN      ! Conservation test
323         !                                       !-------------------
324         numce_dh  = 0
325         meance_dh = 0._wp
326         DO ji = kideb, kiut
327            IF ( ( z_f_surf(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN
328               numce_dh  = numce_dh + 1
329               meance_dh = meance_dh + z_f_surf(ji) + zdq_i(ji)
330            ENDIF
331            IF( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN!
332               WRITE(numout,*) ' ALERTE heat loss for surface melt '
333               WRITE(numout,*) ' ii, ij, jl :', ii, ij, jl
334               WRITE(numout,*) ' ht_i_b       : ', ht_i_b(ji)
335               WRITE(numout,*) ' z_f_surf     : ', z_f_surf(ji)
336               WRITE(numout,*) ' zdq_i        : ', zdq_i(ji)
337               WRITE(numout,*) ' ht_i_b       : ', ht_i_b(ji)
338               WRITE(numout,*) ' fc_bo_i      : ', fc_bo_i(ji)
339               WRITE(numout,*) ' fbif_1d      : ', fbif_1d(ji)
340               WRITE(numout,*) ' qlbbq_1d     : ', qlbbq_1d(ji)
341               WRITE(numout,*) ' s_i_new      : ', s_i_new(ji)
342               WRITE(numout,*) ' sss_m        : ', sss_m(ii,ij)
343            ENDIF
344         END DO
345         !
346         IF( numce_dh > 0 )   meance_dh = meance_dh / numce_dh
347         WRITE(numout,*) ' Error report - Category : ', jl
348         WRITE(numout,*) ' ~~~~~~~~~~~~ '
349         WRITE(numout,*) ' Number of points where there is sur. me. error : ', numce_dh
350         WRITE(numout,*) ' Mean basal growth error on error points : ', meance_dh
351         !
352      ENDIF
353
354      !----------------------
355      ! 3.3 Snow sublimation
356      !----------------------
357
358      DO ji = kideb, kiut
359         ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
360#if defined key_coupled
361         zdh_s_sub(ji)    =  0._wp      ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice)
362#else
363         !                              ! forced  mode: snow thickness change due to sublimation
364         zdh_s_sub(ji)    =  - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice
365#endif
366         dh_s_tot (ji)    =  dh_s_tot(ji) + zdh_s_sub(ji)
367         zdhcf            =  ht_s_b(ji) + zdh_s_sub(ji) 
368         ht_s_b   (ji)    =  MAX( zzero , zdhcf )
369         ! we recompute dh_s_tot
370         dh_s_tot (ji)    =  ht_s_b(ji) - zhsold(ji)
371         qt_s_in  (ji,jl) =  qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1)
372      END DO
373
374      zqt_dummy(:) = 0.e0
375      DO jk = 1, nlay_s
376         DO ji = kideb,kiut
377            q_s_b    (ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
378            zqt_dummy(ji)    =  zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s )            ! heat conservation
379         END DO
380      END DO
381
382      DO jk = 1, nlay_s
383         DO ji = kideb, kiut
384            ! In case of disparition of the snow, we have to update the snow temperatures
385            zhisn  =  MAX(  zzero , SIGN( zone, - ht_s_b(ji) )  )
386            t_s_b(ji,jk) = ( 1.0 - zhisn ) * t_s_b(ji,jk) + zhisn * rtt
387            q_s_b(ji,jk) = ( 1.0 - zhisn ) * q_s_b(ji,jk)
388         END DO
389      END DO
390
391      !
392      !------------------------------------------------------------------------------!
393      ! 4) Basal growth / melt                                                       !
394      !------------------------------------------------------------------------------!
395      !
396      !------------------
397      ! 4.1 Basal growth
398      !------------------
399      ! Basal growth is driven by heat imbalance at the ice-ocean interface,
400      ! between the inner conductive flux  (fc_bo_i), from the open water heat flux
401      ! (qlbbqb) and the turbulent ocean flux (fbif).
402      ! fc_bo_i is positive downwards. fbif and qlbbq are positive to the ice
403
404      ! If salinity varies in time, an iterative procedure is required, because
405      ! the involved quantities are inter-dependent.
406      ! Basal growth (dh_i_bott) depends upon new ice specific enthalpy (zEi),
407      ! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bott)
408      ! -> need for an iterative procedure, which converges quickly
409
410      IF ( num_sal == 2 ) THEN
411         num_iter_max = 5
412      ELSE
413         num_iter_max = 1
414      ENDIF
415
416      ! Initialize dh_i_bott
417      dh_i_bott(:) = 0.0e0
418
419      ! Iterative procedure
420      DO iter = 1, num_iter_max
421         DO ji = kideb, kiut
422
423            IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0  ) THEN
424
425               ! New bottom ice salinity (Cox & Weeks, JGR88 )
426               !--- zswi1  if dh/dt < 2.0e-8
427               !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
428               !--- zswi2  if dh/dt > 3.6e-7
429               zgrr               = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi13 ) )
430               zswi2              = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) )
431               zswi12             = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
432               zswi1              = 1. - zswi2 * zswi12
433               zfracs             = MIN ( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
434                  &               + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 )
435
436               ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
437
438               s_i_new(ji)        = zswitch_sal * zfracs * sss_m(ii,ij)  &  ! New ice salinity
439                                  + ( 1. - zswitch_sal ) * sm_i_b(ji) 
440               ! New ice growth
441               ztmelts            = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K)
442
443               zt_i_new           = zswitch_sal * t_bo_b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)
444               
445               zEi                = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)     
446                  &               - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) )   &
447                  &               + rcp  * ( ztmelts-rtt )         
448
449               zEw                = rcp  * ( t_bo_b(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
450
451               zdE                = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
452
453               dh_i_bott(ji)      = rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / ( zdE * rhoic )
454
455               !!! not sure we still need the next line... useful to keep this in memory ? check limthd_ent...
456               q_i_b(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0)
457
458            ENDIF ! fc_bo_i
459         END DO ! ji
460      END DO ! iter
461
462      ! Contribution to Energy and Salt Fluxes
463      DO ji = kideb, kiut
464
465         zEw            = rcp  * ( t_bo_b(ji) - rt0 )                   ! Specific enthalpy of seawater (J/kg, < 0)
466
467         zfmdt          = - rhoic * dh_i_bott(ji)                       ! Mass flux x time step (kg/m2, < 0)
468
469         ! Contribution to energy flux to the ocean (J/m2)
470         rdq_ice_1d(ji) = rdq_ice_1d(ji) + zEw         * a_i_b(ji) * zfmdt
471
472         ! Contribution to salt flux  ()
473         sfx_thd_1d(ji) = sfx_thd_1d(ji) + s_i_new(ji) * a_i_b(ji) * zfmdt * r1_rdtice
474
475      END DO
476
477      !----------------
478      ! 4.2 Basal melt
479      !----------------
480      meance_dh = 0._wp
481      numce_dh  = 0
482      zinnermelt(:) = 0._wp
483
484      DO ji = kideb, kiut
485         ! heat convergence at the surface > 0
486         IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0._wp  ) THEN
487            s_i_new(ji) = s_i_b(ji,nlay_i) ! is this line still useful ?
488            zqfont_bo(ji) =  rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) )
489            zfbase(ji)    =  zqfont_bo(ji) * r1_rdtice     ! heat conservation test
490            zdq_i(ji)     =  0._wp
491            dh_i_bott(ji) =  0._wp
492         ENDIF
493      END DO
494
495      DO jk = nlay_i, 1, -1
496         DO ji = kideb, kiut
497            IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  >=  0._wp  ) THEN
498
499               ztmelts = - tmut * s_i_b(ji,jk) + rtt  ! Melting point of layer jk (K)
500
501               IF( t_i_b(ji,jk) >= ztmelts ) THEN !!! Internal melting
502
503                  zdeltah   (ji,jk) = - zh_i(ji)  ! internal melting occurs when the internal temperature is above freezing     
504                                                  ! this should normally not happen, but sometimes, heat diffusion leads to this
505
506                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk)
507
508                  zinnermelt(ji)    = 1._wp
509
510                  zQm               = 0.          ! Not sure which specific enthalpy we should use here (MV HC 2014)
511                                                  ! If that happens, heat is probably not well counted
512                                                  ! Put zero by default, but more bug analysis should be done to investigate this case
513
514               ELSE                               !!! Basal melting
515
516                  zEi               = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)
517
518                  zEw               = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0)
519
520                  zdE               = zEi - zEw              ! Specific enthalpy difference   (J/kg, <0)
521
522                  zfmdt             = - zqfont_bo(ji) / zdE  ! Mass flux x time step (kg/m2, >0)
523
524                  zdeltah(ji,jk)    = - zfmdt / rhoic        ! Gross thickness change
525
526                  zqfont_bo(ji)     = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * rhoic * ( - zdE ) ! Update heat available
527
528                  zdeltah(ji,jk)    = MAX( zdeltah(ji,jk), - zh_i(ji) ) ! Update thickness change
529
530                  dh_i_bott(ji)     = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt
531
532                  zfmdt             = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step
533
534                  zQm               = MAX ( zfmdt , 0.0 ) * zEw         ! Heat exchanged with ocean
535
536                  zdq_i(ji)         = zdq_i(ji) + zdeltah(ji,jk) * rhoic * zdE * r1_rdtice ! for heat conservation
537
538               ENDIF
539
540               ! Contribution to salt flux
541               sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    &
542                    &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice
543
544               ! Contribution to energy flux to the ocean (J/m2)
545               rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm * a_i_b(ji)
546
547            ENDIF
548         END DO ! ji
549      END DO ! jk
550
551      !                                          !-------------------
552      IF( con_i .AND. jiindex_1d > 0 ) THEN      ! Conservation test
553      !                                          !-------------------
554         DO ji = kideb, kiut
555            IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0.e0  ) THEN
556               IF( ( zfbase(ji) + zdq_i(ji) ) >= 1.e-3 ) THEN
557                  numce_dh  = numce_dh + 1
558                  meance_dh = meance_dh + zfbase(ji) + zdq_i(ji)
559               ENDIF
560               IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN
561                  WRITE(numout,*) ' ALERTE heat loss for basal melt : ii, ij, jl :', ii, ij, jl
562                  WRITE(numout,*) ' ht_i_b    : ', ht_i_b(ji)
563                  WRITE(numout,*) ' zfbase    : ', zfbase(ji)
564                  WRITE(numout,*) ' zdq_i     : ', zdq_i(ji)
565                  WRITE(numout,*) ' ht_i_b    : ', ht_i_b(ji)
566                  WRITE(numout,*) ' fc_bo_i   : ', fc_bo_i(ji)
567                  WRITE(numout,*) ' fbif_1d   : ', fbif_1d(ji)
568                  WRITE(numout,*) ' qlbbq_1d  : ', qlbbq_1d(ji)
569                  WRITE(numout,*) ' s_i_new   : ', s_i_new(ji)
570                  WRITE(numout,*) ' sss_m     : ', sss_m(ii,ij)
571                  WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
572                  WRITE(numout,*) ' innermelt : ', INT( zinnermelt(ji) )
573               ENDIF
574            ENDIF
575         END DO
576         IF( numce_dh > 0 )   meance_dh = meance_dh / numce_dh
577         WRITE(numout,*) ' Number of points where there is bas. me. error : ', numce_dh
578         WRITE(numout,*) ' Mean basal melt error on error points : ', meance_dh
579         WRITE(numout,*) ' Remaining bottom heat : ', zqfont_bo(jiindex_1d)
580         !
581      ENDIF
582
583      !
584      !------------------------------------------------------------------------------!
585      !  5) Pathological cases                                                       !
586      !------------------------------------------------------------------------------!
587      !
588      !----------------------------------------------
589      ! 5.1 Excessive ablation in a 1-category model
590      !----------------------------------------------
591
592      DO ji = kideb, kiut
593         !                     ! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15)
594         IF( jpl == 1 ) THEN   ;   zdhbf = MAX( hmelt , dh_i_bott(ji) )
595         ELSE                  ;   zdhbf =              dh_i_bott(ji) 
596         ENDIF
597         zdvres        = zdhbf - dh_i_bott(ji)
598         dh_i_bott(ji) = zdhbf
599         sfx_thd_1d(ji)  = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice
600         !                     ! excessive energy is sent to lateral ablation
601         zinda = MAX( 0._wp, SIGN( 1._wp , 1.0 - at_i_b(ji) - epsi10 ) )
602         fsup(ji) =  zinda * rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi10 ) * zdvres * r1_rdtice
603      END DO
604
605      !-----------------------------------
606      ! 5.2 More than available ice melts
607      !-----------------------------------
608      ! then heat applied minus heat content at previous time step should equal heat remaining
609      !
610      DO ji = kideb, kiut
611         ! Adapt the remaining energy if too much ice melts
612         !--------------------------------------------------
613         zdvres     = MAX( 0._wp, - ht_i_b(ji) - dh_i_surf(ji) - dh_i_bott(ji) )
614         zdvsur     = MIN( 0._wp, dh_i_surf(ji) + zdvres ) - dh_i_surf(ji) ! fill the surface first
615         zdvbot     = MAX( 0._wp, zdvres - zdvsur ) ! then the bottom
616         dh_i_surf (ji) = dh_i_surf(ji) + zdvsur ! clem
617         dh_i_bott (ji) = dh_i_bott(ji) + zdvbot ! clem
618
619         ! new ice thickness (clem)
620         zhgnew(ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji)
621         zihgnew    = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice
622         zhgnew(ji) = zihgnew * zhgnew(ji)      ! ice thickness is put to 0
623 
624         !                     !since ice volume is only used for outputs, we keep it global for all categories
625         dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji)
626
627         ! remaining heat
628         zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) +  zqfont_bo(ji) ) 
629
630         ! If snow remains, energy is used to melt snow
631         zhni =  ht_s_b(ji)      ! snow depth at previous time step
632         zihg =  MAX(  zzero , SIGN ( zone , - ht_s_b(ji) )  )   ! =0 if snow
633
634         ! energy of melting of remaining snow
635         zinda = MAX( 0._wp, SIGN( 1._wp , zhni - epsi10 ) )
636         zqt_s(ji) =    ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi10 ) * zinda
637         zdhnm     =  - ( 1. - zihg ) * ( 1. - zihgnew ) * zfdt_final(ji) / MAX( zqt_s(ji) , epsi13 )
638         zhnfi     =  zhni + zdhnm
639         zfdt_final(ji) =  MAX( zfdt_final(ji) + zqt_s(ji) * zdhnm , 0.0 )
640         ht_s_b(ji)     =  MAX( zzero , zhnfi )
641         zqt_s(ji)      =  zqt_s(ji) * ht_s_b(ji)
642         ! we recompute dh_s_tot (clem)
643         dh_s_tot (ji)  =  ht_s_b(ji) - zhsold(ji)
644
645         ! Mass variations of ice and snow
646         !---------------------------------
647         !                                              ! mass variation of the jl category
648         zzfmass_s = - a_i_b(ji) * ( zhni       - ht_s_b(ji) ) * rhosn   ! snow
649         zzfmass_i =   a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic   ! ice 
650         !
651         zfmass_i(ji) = zzfmass_i                       ! ice variation saved to compute salt flux (see below)
652         !
653         !                                              ! mass variation cumulated over category
654         !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s                     ! snow
655         !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i                     ! ice
656
657         ! Remaining heat to the ocean
658         !---------------------------------
659         focea(ji)  = - zfdt_final(ji) * r1_rdtice         ! focea is in W.m-2 * dt
660
661         ! residual salt flux (clem)
662         !--------------------------
663         ! surface
664         sfx_thd_1d(ji)    = sfx_thd_1d(ji) - sm_i_b(ji)  * a_i_b(ji) * zdvsur * rhoic * r1_rdtice
665         ! bottom
666         IF ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) >= 0._wp ) THEN ! melting
667            sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji)  * a_i_b(ji) * zdvbot * rhoic * r1_rdtice
668         ELSE                                                          ! growth
669            sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * zdvbot * rhoic * r1_rdtice
670         ENDIF
671         !
672         ! diagnostic
673         ii = MOD( npb(ji) - 1, jpi ) + 1
674         ij = ( npb(ji) - 1 ) / jpi + 1
675         diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice
676         diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice
677         diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice
678      END DO
679
680      ftotal_fin (:) = zfdt_final(:)  * r1_rdtice
681
682      !---------------------------
683      ! heat fluxes                   
684      !---------------------------
685      DO ji = kideb, kiut
686         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   ! =1 if ice
687         !
688         ! Heat flux
689         ! excessive bottom ablation energy (fsup) - 0 except if jpl = 1
690         ! excessive total  ablation energy (focea) sent to the ocean
691         qfvbq_1d(ji)  = qfvbq_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) * rdt_ice
692
693         zihic   = 1.0 - MAX(  zzero , SIGN( zone , -ht_i_b(ji) )  )      ! equals 0 if ht_i = 0, 1 if ht_i gt 0
694         fscbq_1d(ji) =  a_i_b(ji) * fstbif_1d(ji)
695         qldif_1d(ji)  = qldif_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea   (ji) * a_i_b(ji) * rdt_ice   &
696            &                                    + ( 1.0 - zihic   ) * fscbq_1d(ji)             * rdt_ice
697      END DO  ! ji
698
699      !-------------------------------------------
700      ! Correct temperature, energy and thickness
701      !-------------------------------------------
702      DO ji = kideb, kiut
703         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 
704         t_su_b(ji) =  zihgnew * t_su_b(ji) + ( 1.0 - zihgnew ) * rtt
705      END DO  ! ji
706
707      DO jk = 1, nlay_i
708         DO ji = kideb, kiut
709            zihgnew      =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 
710            t_i_b(ji,jk) =  zihgnew * t_i_b(ji,jk) + ( 1.0 - zihgnew ) * rtt
711            q_i_b(ji,jk) =  zihgnew * q_i_b(ji,jk)
712         END DO
713      END DO  ! ji
714
715      DO ji = kideb, kiut
716         ht_i_b(ji) = zhgnew(ji)
717      END DO ! ji
718     
719      !
720      !------------------------------------------------------------------------------|
721      !  6) Snow-Ice formation                                                       |
722      !------------------------------------------------------------------------------|
723      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
724      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
725      DO ji = kideb, kiut
726         !
727         dh_snowice(ji) = MAX(  zzero , ( rhosn * ht_s_b(ji) + (rhoic-rau0) * ht_i_b(ji) ) / ( rhosn+rau0-rhoic )  )
728         zhgnew(ji)     = MAX(  zhgnew(ji) , zhgnew(ji) + dh_snowice(ji)  )
729         zhnnew         = MIN(  ht_s_b(ji) , ht_s_b(ji) - dh_snowice(ji)  )
730
731         !  Changes in ice volume and ice mass.
732         dvnbq_1d  (ji) =                a_i_b(ji) * ( zhgnew(ji)-ht_i_b(ji) )
733         dmgwi_1d  (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn
734
735         !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic
736         !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew     - ht_s_b(ji) ) * rhosn
737
738         ! Salinity of snow ice
739         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
740
741         zs_snic = zswitch_sal          * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic     &
742                 + ( 1. - zswitch_sal ) * sm_i_b(ji)
743
744         ! entrapment during snow ice formation
745         ! clem: new salinity difference stored (to be used in limthd_ent.F90)
746         IF (  num_sal == 2  ) THEN
747            i_ice_switch = MAX( 0._wp , SIGN( 1._wp , zhgnew(ji) - epsi10 ) )
748            ! salinity dif due to snow-ice formation
749            dsm_i_si_1d(ji) = ( zs_snic - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch     
750            ! salinity dif due to bottom growth
751            IF (  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  < 0._wp ) THEN
752               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch
753            ENDIF
754         ENDIF
755
756         !  Actualize new snow and ice thickness.
757         ht_s_b(ji)  = zhnnew
758         ht_i_b(ji)  = zhgnew(ji)
759
760         ! Total ablation ! new lines added to debug
761         IF( ht_i_b(ji) <= 0._wp )   a_i_b(ji) = 0._wp
762
763         ! diagnostic ( snow ice growth )
764         ii = MOD( npb(ji) - 1, jpi ) + 1 !MV HC 2014 useless
765         ij =    ( npb(ji) - 1 ) / jpi + 1 ! MV HC 2014 useless
766         diag_sni_gr(ii,ij)  = diag_sni_gr(ii,ij) + dh_snowice(ji)*a_i_b(ji) * r1_rdtice
767
768         ! Contribution to energy flux to the ocean [J/m2], >0
769         zfmdt          = ( rhosn - rhoic ) * MAX( dh_snowice(ji), 0.0 )    ! <0
770         zsstK          = sst_m(ii,ij) + rt0                               
771         zEw            = rcp * ( zsstK - rt0 )
772         zQm            = zfmdt * zEw 
773         rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm * a_i_b(ji)
774
775         ! Contribution to salt flux
776         sfx_thd_1d(ji) = sfx_thd_1d(ji) + sss_m(ii,ij) * a_i_b(ji) * zfmdt * r1_rdtice 
777         
778         ! Contribution to mass fluxes
779         rdm_ice_1d(ji) = rdm_ice_1d(ji) + ( a_i_b(ji) * ht_i_b(ji) - zviold(ji) ) * rhoic 
780         rdm_snw_1d(ji) = rdm_snw_1d(ji) + ( a_i_b(ji) * ht_s_b(ji) - zvsold(ji) ) * rhosn 
781
782      END DO !ji
783      !
784      CALL wrk_dealloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i )
785      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy )
786      CALL wrk_dealloc( jpij, zinnermelt, zfbase, zdq_i )
787      CALL wrk_dealloc( jpij, jkmax, zdeltah, zqt_i_lay )
788      !
789      CALL wrk_dealloc( jpij, zviold, zvsold ) ! clem
790      !
791   END SUBROUTINE lim_thd_dh
792   
793#else
794   !!----------------------------------------------------------------------
795   !!   Default option                               NO  LIM3 sea-ice model
796   !!----------------------------------------------------------------------
797CONTAINS
798   SUBROUTINE lim_thd_dh          ! Empty routine
799   END SUBROUTINE lim_thd_dh
800#endif
801
802   !!======================================================================
803END MODULE limthd_dh
Note: See TracBrowser for help on using the repository browser.