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 trunk/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMO/LIM_SRC_3/limthd_dh.F90 @ 1557

Last change on this file since 1557 was 1465, checked in by smasson, 15 years ago

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