source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 4161

Last change on this file since 4161 was 4161, checked in by cetlod, 7 years ago

dev_LOCEAN_2013 : merge in the 3rd dev branch dev_r4028_CNRS_LIM3, see ticket #1169

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