source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 4333

Last change on this file since 4333 was 4333, checked in by clem, 7 years ago

remove remaining bugs in LIM3, so that it can run in both regional and global config

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