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

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

source: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 3625

Last change on this file since 3625 was 3625, checked in by acc, 12 years ago

Branch dev_NOC_2012_r3555. #1006. Step 7. Check in code now merged with dev_r3385_NOCS04_HAMF

  • Property svn:keywords set to Id
File size: 34.7 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         ! if qla is positive (upwards), heat goes to the atmosphere, therefore
318         ! snow sublimates, if qla is negative (downwards), snow condensates
319         zdh_s_sub(ji)    =  - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice
320         dh_s_tot (ji)    =  dh_s_tot(ji) + zdh_s_sub(ji)
321         zdhcf            =  ht_s_b(ji) + zdh_s_sub(ji) 
322         ht_s_b   (ji)    =  MAX( zzero , zdhcf )
323         ! we recompute dh_s_tot
324         dh_s_tot (ji)    =  ht_s_b(ji) - zhsold(ji)
325         qt_s_in  (ji,jl) =  qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1)
326      END DO
327
328      zqt_dummy(:) = 0.e0
329      DO jk = 1, nlay_s
330         DO ji = kideb,kiut
331            q_s_b    (ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
332            zqt_dummy(ji)    =  zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s            ! heat conservation
333         END DO
334      END DO
335
336      DO jk = 1, nlay_s
337         DO ji = kideb, kiut
338            ! In case of disparition of the snow, we have to update the snow temperatures
339            zhisn  =  MAX(  zzero , SIGN( zone, - ht_s_b(ji) )  )
340            t_s_b(ji,jk) = ( 1.0 - zhisn ) * t_s_b(ji,jk) + zhisn * rtt
341            q_s_b(ji,jk) = ( 1.0 - zhisn ) * q_s_b(ji,jk)
342         END DO
343      END DO
344
345      !
346      !------------------------------------------------------------------------------!
347      ! 4) Basal growth / melt                                                       !
348      !------------------------------------------------------------------------------!
349      !
350      ! Ice basal growth / melt is given by the ratio of heat budget over basal
351      ! ice heat content.  Basal heat budget is given by the difference between
352      ! the inner conductive flux  (fc_bo_i), from the open water heat flux
353      ! (qlbbqb) and the turbulent ocean flux (fbif).
354      ! fc_bo_i is positive downwards. fbif and qlbbq are positive to the ice
355
356      !-----------------------------------------------------
357      ! 4.1 Basal growth - (a) salinity not varying in time
358      !-----------------------------------------------------
359      IF(  num_sal /= 2  ) THEN   ! ice salinity constant in time
360         DO ji = kideb, kiut
361            IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0._wp  ) THEN
362               s_i_new(ji)         =  sm_i_b(ji)
363               ! Melting point in K
364               ztmelts             =  - tmut * s_i_new(ji) + rtt 
365               ! New ice heat content (Bitz and Lipscomb, 1999)
366               ztform              =  t_i_b(ji,nlay_i)  ! t_bo_b crashes in the
367               ! Baltic
368               q_i_b(ji,nlay_i+1)  = rhoic * (  cpic * ( ztmelts - ztform )                                &
369                  &                           + lfus * (  1.0 - ( ztmelts - rtt ) / ( ztform - rtt )  )    &
370                  &                           - rcp  * ( ztmelts - rtt )                                 )
371               ! Basal growth rate = - F*dt / q
372               dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
373            ENDIF
374         END DO
375      ENDIF
376
377      !-------------------------------------------------
378      ! 4.1 Basal growth - (b) salinity varying in time
379      !-------------------------------------------------
380      IF(  num_sal == 2  ) THEN
381         ! the growth rate (dh_i_bott) is function of the new ice heat content (q_i_b(nlay_i+1)).
382         ! q_i_b depends on the new ice salinity (snewice).
383         ! snewice depends on dh_i_bott ; it converges quickly, so, no problem
384         ! See Vancoppenolle et al., OM08 for more info on this
385
386         ! Initial value (tested 1D, can be anything between 1 and 20)
387         num_iter_max = 4
388         s_i_new(:)   = 4.0
389
390         ! Iterative procedure
391         DO iter = 1, num_iter_max
392            DO ji = kideb, kiut
393               IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0  ) THEN
394                  ii = MOD( npb(ji) - 1, jpi ) + 1
395                  ij = ( npb(ji) - 1 ) / jpi + 1
396                  ! Melting point in K
397                  ztmelts             =   - tmut * s_i_new(ji) + rtt 
398                  ! New ice heat content (Bitz and Lipscomb, 1999)
399                  q_i_b(ji,nlay_i+1)  =  rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             &
400                     &                            + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )   &
401                     &                            - rcp * ( ztmelts-rtt )                                     )
402                  ! Bottom growth rate = - F*dt / q
403                  dh_i_bott(ji) =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)
404                  ! New ice salinity ( Cox and Weeks, JGR, 1988 )
405                  ! zswi2  (1) if dh_i_bott/rdt .GT. 3.6e-7
406                  ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8
407                  ! zswi1  (1) if dh_i_bott/rdt .LT. 2.0e-8
408                  zgrr   = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi13 ) )
409                  zswi2  = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) ) 
410                  zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
411                  zswi1  = 1. - zswi2 * zswi12 
412                  zfracs = zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
413                     &                   + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 
414                  zds         = zfracs * sss_m(ii,ij) - s_i_new(ji)
415                  s_i_new(ji) = zfracs * sss_m(ii,ij)
416               ENDIF ! fc_bo_i
417            END DO ! ji
418         END DO ! iter
419
420         ! Final values
421         DO ji = kideb, kiut
422            IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN
423               ! New ice salinity must not exceed 15 psu
424               s_i_new(ji) = MIN( s_i_new(ji), s_i_max )
425               ! Metling point in K
426               ztmelts     =   - tmut * s_i_new(ji) + rtt 
427               ! New ice heat content (Bitz and Lipscomb, 1999)
428               q_i_b(ji,nlay_i+1)  =  rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                              &
429                  &                            + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )    &
430                  &                            - rcp * ( ztmelts - rtt )                                    )
431               ! Basal growth rate = - F*dt / q
432               dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)
433               ! Salinity update
434               ! entrapment during bottom growth
435               dsm_i_se_1d(ji) = ( s_i_new(ji) * dh_i_bott(ji) + sm_i_b(ji) * ht_i_b(ji) )    &
436                  &            / MAX( ht_i_b(ji) + dh_i_bott(ji) ,epsi13 ) - sm_i_b(ji)
437            ENDIF ! heat budget
438         END DO
439      ENDIF
440
441      !----------------
442      ! 4.2 Basal melt
443      !----------------
444      meance_dh = 0._wp
445      numce_dh  = 0
446      zinnermelt(:) = 0._wp
447
448      DO ji = kideb, kiut
449         ! heat convergence at the surface > 0
450         IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0._wp  ) THEN
451            s_i_new(ji)   =  s_i_b(ji,nlay_i)
452            zqfont_bo(ji) =  rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) )
453            zfbase(ji)    =  zqfont_bo(ji) * r1_rdtice     ! heat conservation test
454            zdq_i(ji)     =  0._wp
455            dh_i_bott(ji) =  0._wp
456         ENDIF
457      END DO
458
459      DO jk = nlay_i, 1, -1
460         DO ji = kideb, kiut
461            IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  >=  0._wp  ) THEN
462               ztmelts = - tmut * s_i_b(ji,jk) + rtt 
463               IF( t_i_b(ji,jk) >= ztmelts ) THEN   !!gm : a comment is needed
464                  zdeltah   (ji,jk) = - zh_i(ji)
465                  dh_i_bott (ji   ) = dh_i_bott(ji) + zdeltah(ji,jk)
466                  zinnermelt(ji   ) = 1._wp
467               ELSE                                  ! normal ablation
468                  zdeltah  (ji,jk) = - zqfont_bo(ji) / q_i_b(ji,jk)
469                  zqfont_bo(ji   ) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk)
470                  zdeltah  (ji,jk) = MAX(zdeltah(ji,jk), - zh_i(ji) )
471                  dh_i_bott(ji   ) = dh_i_bott(ji) + zdeltah(ji,jk)
472                  zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice
473               ENDIF
474               ! contribution to salt flux
475               zsfx_melt(ji) = zsfx_melt(ji) - sm_i_b(ji) * a_i_b(ji) * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 
476            ENDIF
477         END DO ! ji
478      END DO ! jk
479
480      !                     !-------------------
481      IF( con_i ) THEN      ! Conservation test
482      !                     !-------------------
483         DO ji = kideb, kiut
484            IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0.e0  ) THEN
485               IF( ( zfbase(ji) + zdq_i(ji) ) >= 1.e-3 ) THEN
486                  numce_dh  = numce_dh + 1
487                  meance_dh = meance_dh + zfbase(ji) + zdq_i(ji)
488               ENDIF
489               IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN
490                  WRITE(numout,*) ' ALERTE heat loss for basal melt : ii, ij, jl :', ii, ij, jl
491                  WRITE(numout,*) ' ht_i_b    : ', ht_i_b(ji)
492                  WRITE(numout,*) ' zfbase    : ', zfbase(ji)
493                  WRITE(numout,*) ' zdq_i     : ', zdq_i(ji)
494                  WRITE(numout,*) ' ht_i_b    : ', ht_i_b(ji)
495                  WRITE(numout,*) ' fc_bo_i   : ', fc_bo_i(ji)
496                  WRITE(numout,*) ' fbif_1d   : ', fbif_1d(ji)
497                  WRITE(numout,*) ' qlbbq_1d  : ', qlbbq_1d(ji)
498                  WRITE(numout,*) ' s_i_new   : ', s_i_new(ji)
499                  WRITE(numout,*) ' sss_m     : ', sss_m(ii,ij)
500                  WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
501                  WRITE(numout,*) ' innermelt : ', INT( zinnermelt(ji) )
502               ENDIF
503            ENDIF
504         END DO
505         IF( numce_dh > 0 )   meance_dh = meance_dh / numce_dh
506         WRITE(numout,*) ' Number of points where there is bas. me. error : ', numce_dh
507         WRITE(numout,*) ' Mean basal melt error on error points : ', meance_dh
508         WRITE(numout,*) ' Remaining bottom heat : ', zqfont_bo(jiindex_1d)
509         !
510      ENDIF
511
512      !
513      !------------------------------------------------------------------------------!
514      !  5) Pathological cases                                                       !
515      !------------------------------------------------------------------------------!
516      !
517      !----------------------------------------------
518      ! 5.1 Excessive ablation in a 1-category model
519      !----------------------------------------------
520
521      DO ji = kideb, kiut
522         !                     ! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15)
523         IF( jpl == 1 ) THEN   ;   zdhbf = MAX( hmelt , dh_i_bott(ji) )
524         ELSE                  ;   zdhbf =              dh_i_bott(ji) 
525         ENDIF
526         !                     ! excessive energy is sent to lateral ablation
527         fsup     (ji) =  rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 )   &
528            &                          * ( zdhbf - dh_i_bott(ji) ) * r1_rdtice
529         dh_i_bott(ji)  = zdhbf
530         !                     !since ice volume is only used for outputs, we keep it global for all categories
531         dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji)
532         !                     !new ice thickness
533         zhgnew   (ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji)
534         !                     ! diagnostic ( bottom ice growth )
535         ii = MOD( npb(ji) - 1, jpi ) + 1
536         ij = ( npb(ji) - 1 ) / jpi + 1
537         diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice
538         diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice
539         diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice
540      END DO
541
542      !-----------------------------------
543      ! 5.2 More than available ice melts
544      !-----------------------------------
545      ! then heat applied minus heat content at previous time step should equal heat remaining
546      !
547      DO ji = kideb, kiut
548         ! Adapt the remaining energy if too much ice melts
549         !--------------------------------------------------
550         zihgnew =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   ! =1 if ice
551         ! 0 if no more ice
552         zhgnew    (ji) =         zihgnew   * zhgnew(ji)      ! ice thickness is put to 0
553         ! remaining heat
554         zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) +  zqfont_bo(ji) ) 
555
556         ! If snow remains, energy is used to melt snow
557         zhni =  ht_s_b(ji)      ! snow depth at previous time step
558         zihg =  MAX(  zzero , SIGN ( zone , - ht_s_b(ji) )  )   ! =0 if snow
559
560         ! energy of melting of remaining snow
561         zqt_s(ji) =    ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi13 )
562         zdhnm     =  - ( 1. - zihg ) * ( 1. - zihgnew ) * zfdt_final(ji) / MAX( zqt_s(ji) , epsi13 )
563         zhnfi     =  zhni + zdhnm
564         zfdt_final(ji) =  MAX( zfdt_final(ji) + zqt_s(ji) * zdhnm , 0.0 )
565         ht_s_b(ji)     =  MAX( zzero , zhnfi )
566         zqt_s(ji)      =  zqt_s(ji) * ht_s_b(ji)
567
568         ! Mass variations of ice and snow
569         !---------------------------------
570         !                                              ! mass variation of the jl category
571         zzfmass_s = - a_i_b(ji) * ( zhni       - ht_s_b(ji) ) * rhosn   ! snow
572         zzfmass_i =   a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic   ! ice 
573         !
574         zfmass_i(ji) = zzfmass_i                       ! ice variation saved to compute salt flux (see below)
575         !
576         !                                              ! mass variation cumulated over category
577         rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s                     ! snow
578         rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i                     ! ice
579
580         ! Remaining heat to the ocean
581         !---------------------------------
582         focea(ji)  = - zfdt_final(ji) * r1_rdtice         ! focea is in W.m-2 * dt
583
584      END DO
585
586      ftotal_fin (:) = zfdt_final(:)  * r1_rdtice
587
588      !---------------------------
589      ! Salt flux and heat fluxes                   
590      !---------------------------
591      DO ji = kideb, kiut
592         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   ! =1 if ice
593         !
594         ! Salt flux
595         sfx_thd_1d(ji) = sfx_thd_1d(ji) +        zihgnew  * zsfx_melt(ji)               &
596            &                            - (1.0 - zihgnew) * zfmass_i (ji) * sm_i_b(ji)  * r1_rdtice
597         !
598         ! Heat flux
599         ! excessive bottom ablation energy (fsup) - 0 except if jpl = 1
600         ! excessive total  ablation energy (focea) sent to the ocean
601         qfvbq_1d(ji)  = qfvbq_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) * rdt_ice
602
603         zihic   = 1.0 - MAX(  zzero , SIGN( zone , -ht_i_b(ji) )  )      ! equals 0 if ht_i = 0, 1 if ht_i gt 0
604         fscbq_1d(ji) =  a_i_b(ji) * fstbif_1d(ji)
605         qldif_1d(ji)  = qldif_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea   (ji) * a_i_b(ji) * rdt_ice   &
606            &                                    + ( 1.0 - zihic   ) * fscbq_1d(ji)             * rdt_ice
607      END DO  ! ji
608
609      !-------------------------------------------
610      ! Correct temperature, energy and thickness
611      !-------------------------------------------
612      DO ji = kideb, kiut
613         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 
614         t_su_b(ji) =  zihgnew * t_su_b(ji) + ( 1.0 - zihgnew ) * rtt
615      END DO  ! ji
616
617      DO jk = 1, nlay_i
618         DO ji = kideb, kiut
619            zihgnew      =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 
620            t_i_b(ji,jk) =  zihgnew * t_i_b(ji,jk) + ( 1.0 - zihgnew ) * rtt
621            q_i_b(ji,jk) =  zihgnew * q_i_b(ji,jk)
622         END DO
623      END DO  ! ji
624
625      DO ji = kideb, kiut
626         ht_i_b(ji) = zhgnew(ji)
627      END DO  ! ji
628      !
629      !------------------------------------------------------------------------------|
630      !  6) Snow-Ice formation                                                       |
631      !------------------------------------------------------------------------------|
632      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
633      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
634      DO ji = kideb, kiut
635         !
636         dh_snowice(ji) = MAX(  zzero , ( rhosn * ht_s_b(ji) + (rhoic-rau0) * ht_i_b(ji) ) / ( rhosn+rau0-rhoic )  )
637         zhgnew(ji)     = MAX(  zhgnew(ji) , zhgnew(ji) + dh_snowice(ji)  )
638         zhnnew         = MIN(  ht_s_b(ji) , ht_s_b(ji) - dh_snowice(ji)  )
639
640         !  Changes in ice volume and ice mass.
641         dvnbq_1d  (ji) =                a_i_b(ji) * ( zhgnew(ji)-ht_i_b(ji) )
642         dmgwi_1d  (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn
643
644         ! All snow is thrown in the ocean, and seawater is taken to replace the volume
645         rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic * ( 1. - rhosn / rhoic )
646         rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew     - ht_s_b(ji) ) * rhosn
647
648         !        Equivalent salt flux (1) Snow-ice formation component
649         !        -----------------------------------------------------
650         ii = MOD( npb(ji) - 1, jpi ) + 1
651         ij =    ( npb(ji) - 1 ) / jpi + 1
652
653         IF( num_sal == 2 ) THEN   ;   zsm_snowice = sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic
654         ELSE                      ;   zsm_snowice = sm_i_b(ji)   
655         ENDIF
656         sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice
657         !
658         ! entrapment during snow ice formation
659         i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - ht_i_b(ji) + 1.0e-6 ) )
660         isnowic      = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - dh_snowice(ji)      ) ) * i_ice_switch
661         IF(  num_sal == 2  )   &
662            dsm_i_si_1d(ji) = (  zsm_snowice * dh_snowice(ji)    &
663            &                  + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13 )   &
664            &                  - sm_i_b(ji)  ) * isnowic     
665
666         !  Actualize new snow and ice thickness.
667         ht_s_b(ji)  = zhnnew
668         ht_i_b(ji)  = zhgnew(ji)
669
670         ! Total ablation ! new lines added to debug
671         IF( ht_i_b(ji) <= 0._wp )   a_i_b(ji) = 0._wp
672
673         ! diagnostic ( snow ice growth )
674         ii = MOD( npb(ji) - 1, jpi ) + 1
675         ij =    ( npb(ji) - 1 ) / jpi + 1
676         diag_sni_gr(ii,ij)  = diag_sni_gr(ii,ij) + dh_snowice(ji)*a_i_b(ji) * r1_rdtice
677         !
678      END DO !ji
679      !
680      CALL wrk_dealloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i )
681      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 )
682      CALL wrk_dealloc( jpij, zinnermelt, zfbase, zdq_i )
683      CALL wrk_dealloc( jpij, jkmax, zdeltah, zqt_i_lay )
684      !
685   END SUBROUTINE lim_thd_dh
686   
687#else
688   !!----------------------------------------------------------------------
689   !!   Default option                               NO  LIM3 sea-ice model
690   !!----------------------------------------------------------------------
691CONTAINS
692   SUBROUTINE lim_thd_dh          ! Empty routine
693   END SUBROUTINE lim_thd_dh
694#endif
695
696   !!======================================================================
697END MODULE limthd_dh
Note: See TracBrowser for help on using the repository browser.