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

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

source: branches/dev_003_CPL/NEMO/LIM_SRC_3/limthd_dh.F90 @ 991

Last change on this file since 991 was 991, checked in by smasson, 16 years ago

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