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/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 4733

Last change on this file since 4733 was 4733, checked in by vancop, 10 years ago

Fix energy budget in coupled case

  • Property svn:keywords set to Id
File size: 36.4 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 wfx_snw & wfx_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   USE cpl_oasis3, ONLY : lk_cpl
29   
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   lim_thd_dh   ! called by lim_thd
34
35   REAL(wp) ::   epsi20 = 1.e-20   ! constant values
36   REAL(wp) ::   epsi10 = 1.e-10   !
37
38   !!----------------------------------------------------------------------
39   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
40   !! $Id$
41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE lim_thd_dh( kideb, kiut )
46      !!------------------------------------------------------------------
47      !!                ***  ROUTINE lim_thd_dh  ***
48      !!
49      !! ** Purpose :   determines variations of ice and snow thicknesses.
50      !!
51      !! ** Method  :   Ice/Snow surface melting arises from imbalance in surface fluxes
52      !!              Bottom accretion/ablation arises from flux budget
53      !!              Snow thickness can increase by precipitation and decrease by sublimation
54      !!              If snow load excesses Archmiede limit, snow-ice is formed by
55      !!              the flooding of sea-water in the snow
56      !!
57      !!                 1) Compute available flux of heat for surface ablation
58      !!                 2) Compute snow and sea ice enthalpies
59      !!                 3) Surface ablation and sublimation
60      !!                 4) Bottom accretion/ablation
61      !!                 5) Case of Total ablation
62      !!                 6) Snow ice formation
63      !!
64      !! References : Bitz and Lipscomb, 1999, J. Geophys. Res.
65      !!              Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
66      !!              Vancoppenolle, Fichefet and Bitz, 2005, Geophys. Res. Let.
67      !!              Vancoppenolle et al.,2009, Ocean Modelling
68      !!------------------------------------------------------------------
69      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
70      !!
71      INTEGER  ::   ji , jk        ! dummy loop indices
72      INTEGER  ::   ii, ij         ! 2D corresponding indices to ji
73      INTEGER  ::   iter
74
75      REAL(wp) ::   ztmelts             ! local scalar
76      REAL(wp) ::   zdh, zfdum  !
77      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
78      REAL(wp) ::   zcoeff       ! dummy argument for snowfall partitioning over ice and leads
79      REAL(wp) ::   zs_snic  ! snow-ice salinity
80      REAL(wp) ::   zswi1        ! switch for computation of bottom salinity
81      REAL(wp) ::   zswi12       ! switch for computation of bottom salinity
82      REAL(wp) ::   zswi2        ! switch for computation of bottom salinity
83      REAL(wp) ::   zgrr         ! bottom growth rate
84      REAL(wp) ::   zt_i_new     ! bottom formation temperature
85
86      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
87      REAL(wp) ::   zEi          ! specific enthalpy of sea ice (J/kg)
88      REAL(wp) ::   zEw          ! specific enthalpy of exchanged water (J/kg)
89      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg)
90      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean
91      REAL(wp) ::   zsstK        ! SST in Kelvin
92
93      REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness
94      REAL(wp), POINTER, DIMENSION(:) ::   zqprec      ! energy of fallen snow                       (J.m-3)
95      REAL(wp), POINTER, DIMENSION(:) ::   zq_su       ! heat for surface ablation                   (J.m-2)
96      REAL(wp), POINTER, DIMENSION(:) ::   zq_bo       ! heat for bottom ablation                    (J.m-2)
97      REAL(wp), POINTER, DIMENSION(:) ::   zq_1cat     ! corrected heat in case 1-cat and hmelt>15cm (J.m-2)
98      REAL(wp), POINTER, DIMENSION(:) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2)
99      REAL(wp), POINTER, DIMENSION(:) ::   zf_tt     ! Heat budget to determine melting or freezing(W.m-2)
100      INTEGER , POINTER, DIMENSION(:) ::   icount      ! number of layers vanished by melting
101
102      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel   ! snow melt
103      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_pre   ! snow precipitation
104      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_sub   ! snow sublimation
105
106      REAL(wp), POINTER, DIMENSION(:,:) ::   zdeltah
107      REAL(wp), POINTER, DIMENSION(:,:) ::   zh_i      ! ice layer thickness
108
109      REAL(wp), POINTER, DIMENSION(:) ::   zqh_i       ! total ice heat content  (J.m-2)
110      REAL(wp), POINTER, DIMENSION(:) ::   zqh_s       ! total snow heat content (J.m-2)
111      REAL(wp), POINTER, DIMENSION(:) ::   zq_s        ! total snow enthalpy     (J.m-3)
112
113      ! mass and salt flux (clem)
114      REAL(wp) :: zdvres, zswitch_sal, zswitch
115
116      ! Heat conservation
117      INTEGER  ::   num_iter_max
118      REAL(wp) ::   zinda, zindq, zindh 
119      REAL(wp), POINTER, DIMENSION(:) ::   zintermelt   ! debug
120
121      !!------------------------------------------------------------------
122
123      ! Discriminate between varying salinity (num_sal=2) and prescribed cases (other values)
124      SELECT CASE( num_sal )                       ! varying salinity or not
125         CASE( 1, 3, 4 ) ;   zswitch_sal = 0       ! prescribed salinity profile
126         CASE( 2 )       ;   zswitch_sal = 1       ! varying salinity profile
127      END SELECT
128
129      CALL wrk_alloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema )
130      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
131      CALL wrk_alloc( jpij, zintermelt )
132      CALL wrk_alloc( jpij, jkmax, zdeltah, zh_i )
133      CALL wrk_alloc( jpij, icount )
134     
135      dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp
136      dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp   
137 
138      zqprec (:) = 0._wp ; zq_su  (:) = 0._wp ; zq_bo  (:) = 0._wp ; zf_tt  (:) = 0._wp
139      zq_1cat(:) = 0._wp ; zq_rema(:) = 0._wp
140
141      zh_s     (:) = 0._wp       
142      zdh_s_pre(:) = 0._wp
143      zdh_s_mel(:) = 0._wp
144      zdh_s_sub(:) = 0._wp
145      zqh_s    (:) = 0._wp     
146      zqh_i    (:) = 0._wp   
147
148      zh_i      (:,:) = 0._wp       
149      zdeltah   (:,:) = 0._wp       
150      zintermelt(:)   = 0._wp
151      icount    (:)   = 0
152
153      ! initialize layer thicknesses and enthalpies
154      h_i_old (:,0:nlay_i+1) = 0._wp
155      qh_i_old(:,0:nlay_i+1) = 0._wp
156      DO jk = 1, nlay_i
157         DO ji = kideb, kiut
158            h_i_old (ji,jk) = ht_i_b(ji) / REAL( nlay_i )
159            qh_i_old(ji,jk) = q_i_b(ji,jk) * h_i_old(ji,jk)
160         ENDDO
161      ENDDO
162      !
163      !------------------------------------------------------------------------------!
164      !  1) Calculate available heat for surface and bottom ablation                 !
165      !------------------------------------------------------------------------------!
166      !
167      DO ji = kideb, kiut
168         zinda      = 1._wp - MAX(  0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )
169         ztmelts    = zinda * rtt + ( 1._wp - zinda ) * rtt
170
171         zfdum      = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
172         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 
173
174         zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_b(ji) - ztmelts ) )
175         zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
176      END DO
177
178      !
179      !------------------------------------------------------------------------------!
180      ! If snow temperature is above freezing point, then snow melts
181      ! (should not happen but sometimes it does)
182      !------------------------------------------------------------------------------!
183      DO ji = kideb, kiut
184         IF( t_s_b(ji,1) > rtt ) THEN !!! Internal melting
185            ! Contribution to heat flux to the ocean [W.m-2], < 0 
186            hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_b(ji,1) * ht_s_b(ji) * a_i_b(ji) * r1_rdtice
187            ! Contribution to mass flux
188            wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_b(ji) * a_i_b(ji) * r1_rdtice
189            ! updates
190            ht_s_b(ji)   = 0._wp
191            q_s_b (ji,1) = 0._wp
192            t_s_b (ji,1) = rtt
193         END IF
194      END DO
195
196      !------------------------------------------------------------!
197      !  2) Computing layer thicknesses and enthalpies.            !
198      !------------------------------------------------------------!
199      !
200      DO ji = kideb, kiut     
201         zh_s(ji) = ht_s_b(ji) / REAL( nlay_s )
202      END DO
203      !
204      DO jk = 1, nlay_s
205         DO ji = kideb, kiut
206            zqh_s(ji) =  zqh_s(ji) + q_s_b(ji,jk) * zh_s(ji)
207         END DO
208      END DO
209      !
210      DO jk = 1, nlay_i
211         DO ji = kideb, kiut
212            zh_i(ji,jk) = ht_i_b(ji) / REAL( nlay_i )
213            zqh_i(ji)   = zqh_i(ji) + q_i_b(ji,jk) * zh_i(ji,jk)
214         END DO
215      END DO
216      !
217      !------------------------------------------------------------------------------|
218      !  3) Surface ablation and sublimation                                         |
219      !------------------------------------------------------------------------------|
220      !
221      !-------------------------
222      ! 3.1 Snow precips / melt
223      !-------------------------
224      ! Snow accumulation in one thermodynamic time step
225      ! snowfall is partitionned between leads and ice
226      ! if snow fall was uniform, a fraction (1-at_i) would fall into leads
227      ! but because of the winds, more snow falls on leads than on sea ice
228      ! and a greater fraction (1-at_i)^beta of the total mass of snow
229      ! (beta < 1) falls in leads.
230      ! In reality, beta depends on wind speed,
231      ! and should decrease with increasing wind speed but here, it is
232      ! considered as a constant. an average value is 0.66
233      ! Martin Vancoppenolle, December 2006
234
235      DO ji = kideb, kiut
236         !-----------
237         ! Snow fall
238         !-----------
239         ! thickness change
240         zcoeff = ( 1._wp - ( 1._wp - at_i_b(ji) )**betas ) / at_i_b(ji) 
241         zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn
242         ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K)
243         zqprec   (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )   
244         IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp
245         ! heat flux from snow precip (>0, W.m-2)
246         hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice
247         ! mass flux, <0
248         wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_b(ji) * zdh_s_pre(ji) * r1_rdtice
249         ! update thickness
250         ht_s_b    (ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_pre(ji) )
251
252         !---------------------
253         ! Melt of falling snow
254         !---------------------
255         ! thickness change
256         IF( zdh_s_pre(ji) > 0._wp ) THEN
257         zindq          = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zqprec(ji) + epsi20 ) )
258         zdh_s_mel (ji) = - zindq * zq_su(ji) / MAX( zqprec(ji) , epsi20 )
259         zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting
260         ! heat used to melt snow (W.m-2, >0)
261         hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdh_s_mel(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice
262         ! snow melting only = water into the ocean (then without snow precip), >0
263         wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdh_s_mel(ji) * r1_rdtice
264         
265         ! updates available heat + thickness
266         zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) )     
267         ht_s_b(ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_mel(ji) )
268         zh_s  (ji) = ht_s_b(ji) / REAL( nlay_s )
269
270         ENDIF
271      END DO
272
273      ! If heat still available, then melt more snow
274      zdeltah(:,:) = 0._wp ! important
275      DO jk = 1, nlay_s
276         DO ji = kideb, kiut
277            ! thickness change
278            zindh            = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_b(ji) ) ) 
279            zindq            = 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_b(ji,jk) + epsi20 ) ) 
280            zdeltah  (ji,jk) = - zindh * zindq * zq_su(ji) / MAX( q_s_b(ji,jk), epsi20 )
281            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting
282            zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)   
283            ! heat used to melt snow(W.m-2, >0)
284            hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_b(ji) * q_s_b(ji,jk) * r1_rdtice 
285            ! snow melting only = water into the ocean (then without snow precip)
286            wfx_snw_1d(ji)   = wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
287
288            ! updates available heat + thickness
289            zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_b(ji,jk) )
290            ht_s_b(ji) = MAX( 0._wp , ht_s_b(ji) + zdeltah(ji,jk) )
291
292         END DO
293      END DO
294
295      !----------------------
296      ! 3.2 Snow sublimation
297      !----------------------
298      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
299      ! clem comment: not counted in mass exchange in limsbc since this is an exchange with atm. (not ocean)
300      ! clem comment: ice should also sublimate
301      IF( lk_cpl ) THEN
302         ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice)
303         zdh_s_sub(:)      =  0._wp 
304      ELSE
305         ! forced  mode: snow thickness change due to sublimation
306         DO ji = kideb, kiut
307            zdh_s_sub(ji)  =  MAX( - ht_s_b(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice )
308            ! Heat flux by sublimation [W.m-2], < 0
309            !      sublimate first snow that had fallen, then pre-existing snow
310            zcoeff         =      ( MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) )   * zqprec(ji) +   &
311               &  ( zdh_s_sub(ji) - MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) ) * q_s_b(ji,1) )  &
312               &  * a_i_b(ji) * r1_rdtice
313            hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff
314            ! Mass flux by sublimation
315            wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_b(ji) * zdh_s_sub(ji) * r1_rdtice
316            ! new snow thickness
317            ht_s_b(ji)     =  MAX( 0._wp , ht_s_b(ji) + zdh_s_sub(ji) )
318         END DO
319      ENDIF
320
321      ! --- Update snow diags --- !
322      DO ji = kideb, kiut
323         dh_s_tot(ji)   = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
324         zh_s(ji)       = ht_s_b(ji) / REAL( nlay_s )
325      END DO ! ji
326
327      !-------------------------------------------
328      ! 3.3 Update temperature, energy
329      !-------------------------------------------
330      ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
331      zq_s(:) = 0._wp 
332      DO jk = 1, nlay_s
333         DO ji = kideb,kiut
334            zindh  =  MAX(  0._wp , SIGN( 1._wp, - ht_s_b(ji) + epsi20 )  )
335            q_s_b(ji,jk) = ( 1._wp - zindh ) / MAX( ht_s_b(ji), epsi20 ) *             &
336              &            ( (   MAX( 0._wp, dh_s_tot(ji) )              ) * zqprec(ji) +  &
337              &              ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_b(ji) ) * rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) )
338            zq_s(ji)     =  zq_s(ji) + q_s_b(ji,jk)
339         END DO
340      END DO
341
342      !--------------------------
343      ! 3.4 Surface ice ablation
344      !--------------------------
345      zdeltah(:,:) = 0._wp ! important
346      DO jk = 1, nlay_i
347         DO ji = kideb, kiut 
348            zEi            = - q_i_b(ji,jk) / rhoic                ! Specific enthalpy of layer k [J/kg, <0]
349
350            ztmelts        = - tmut * s_i_b(ji,jk) + rtt           ! Melting point of layer k [K]
351
352            zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0]
353
354            zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0
355
356            zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0]
357
358            zdeltah(ji,jk) = - zfmdt / rhoic                       ! Melt of layer jk [m, <0]
359
360            zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0]
361
362            zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat
363
364            dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt
365
366            zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0]
367
368            zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0]
369
370            ! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok)
371            sfx_sum_1d(ji)   = sfx_sum_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
372
373            ! Contribution to heat flux [W.m-2], < 0
374            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
375
376            ! Total heat flux used in this process [W.m-2], > 0 
377            hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
378
379            ! Contribution to mass flux
380            wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
381           
382            ! record which layers have disappeared (for bottom melting)
383            !    => icount=0 : no layer has vanished
384            !    => icount=5 : 5 layers have vanished
385            zindh       = NINT( MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) ) 
386            icount(ji)  = icount(ji) + zindh
387            zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
388
389            ! update heat content (J.m-2) and layer thickness
390            qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_b(ji,jk)
391            h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
392         END DO
393      END DO
394      ! update ice thickness
395      DO ji = kideb, kiut
396         ht_i_b(ji) =  MAX( 0._wp , ht_i_b(ji) + dh_i_surf(ji) )
397      END DO
398
399      !
400      !------------------------------------------------------------------------------!
401      ! 4) Basal growth / melt                                                       !
402      !------------------------------------------------------------------------------!
403      !
404      !------------------
405      ! 4.1 Basal growth
406      !------------------
407      ! Basal growth is driven by heat imbalance at the ice-ocean interface,
408      ! between the inner conductive flux  (fc_bo_i), from the open water heat flux
409      ! (fhld) and the turbulent ocean flux (fhtur).
410      ! fc_bo_i is positive downwards. fhtur and fhld are positive to the ice
411
412      ! If salinity varies in time, an iterative procedure is required, because
413      ! the involved quantities are inter-dependent.
414      ! Basal growth (dh_i_bott) depends upon new ice specific enthalpy (zEi),
415      ! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bott)
416      ! -> need for an iterative procedure, which converges quickly
417
418      IF ( num_sal == 2 ) THEN
419         num_iter_max = 5
420      ELSE
421         num_iter_max = 1
422      ENDIF
423
424      !clem debug. Just to be sure that enthalpy at nlay_i+1 is null
425      DO ji = kideb, kiut
426         q_i_b(ji,nlay_i+1) = 0._wp
427      END DO
428
429      ! Iterative procedure
430      DO iter = 1, num_iter_max
431         DO ji = kideb, kiut
432            IF(  zf_tt(ji) < 0._wp  ) THEN
433
434               ! New bottom ice salinity (Cox & Weeks, JGR88 )
435               !--- zswi1  if dh/dt < 2.0e-8
436               !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
437               !--- zswi2  if dh/dt > 3.6e-7
438               zgrr               = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) )
439               zswi2              = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
440               zswi12             = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
441               zswi1              = 1. - zswi2 * zswi12
442               zfracs             = MIN ( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
443                  &               + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 )
444
445               ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
446
447               s_i_new(ji)        = zswitch_sal * zfracs * sss_m(ii,ij)  &  ! New ice salinity
448                                  + ( 1. - zswitch_sal ) * sm_i_b(ji) 
449               ! New ice growth
450               ztmelts            = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K)
451
452               zt_i_new           = zswitch_sal * t_bo_b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)
453               
454               zEi                = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)     
455                  &               - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) )   &
456                  &               + rcp  * ( ztmelts-rtt )         
457
458               zEw                = rcp  * ( t_bo_b(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
459
460               zdE                = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
461
462               dh_i_bott(ji)      = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
463
464               q_i_b(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0)
465               
466            ENDIF ! fc_bo_i
467         END DO ! ji
468      END DO ! iter
469
470      ! Contribution to Energy and Salt Fluxes
471      DO ji = kideb, kiut
472         IF(  zf_tt(ji) < 0._wp  ) THEN
473            ! New ice growth
474                                   
475            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0)
476
477            ztmelts        = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K)
478           
479            zt_i_new       = zswitch_sal * t_bo_b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)
480           
481            zEi            = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)     
482               &               - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) )   &
483               &               + rcp  * ( ztmelts-rtt )         
484           
485            zEw            = rcp  * ( t_bo_b(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
486           
487            zdE            = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
488           
489            ! Contribution to heat flux to the ocean [W.m-2], >0 
490            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
491
492            ! Total heat flux used in this process [W.m-2], <0 
493            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
494           
495            ! Contribution to salt flux, <0
496            sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_b(ji) * zfmdt * r1_rdtice
497
498            ! Contribution to mass flux, <0
499            wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_b(ji) * dh_i_bott(ji) * r1_rdtice
500
501            ! update heat content (J.m-2) and layer thickness
502            qh_i_old(ji,nlay_i+1) = qh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * q_i_b(ji,nlay_i+1)
503            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji)
504         ENDIF
505      END DO
506
507      !----------------
508      ! 4.2 Basal melt
509      !----------------
510      zdeltah(:,:) = 0._wp ! important
511      DO jk = nlay_i, 1, -1
512         DO ji = kideb, kiut
513            IF(  zf_tt(ji)  >=  0._wp  .AND. jk > icount(ji) ) THEN   ! do not calculate where layer has already disappeared from surface melting
514
515               ztmelts = - tmut * s_i_b(ji,jk) + rtt  ! Melting point of layer jk (K)
516
517               IF( t_i_b(ji,jk) >= ztmelts ) THEN !!! Internal melting
518                  zintermelt(ji)    = 1._wp
519
520                  zEi               = - q_i_b(ji,jk) / rhoic        ! Specific enthalpy of melting ice (J/kg, <0)
521
522                  !!zEw               = rcp * ( t_i_b(ji,jk) - rtt )  ! Specific enthalpy of meltwater at T = t_i_b (J/kg, <0)
523
524                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0)
525                                                                    ! set up at 0 since no energy is needed to melt water...(it is already melted)
526
527                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing     
528                                                                   ! this should normally not happen, but sometimes, heat diffusion leads to this
529
530                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk)
531
532                  zfmdt             = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0
533
534                  ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
535                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_b(ji) * zEi * r1_rdtice
536
537                  ! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok)
538                  sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
539                                   
540                  ! Contribution to mass flux
541                  wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
542
543                  ! update heat content (J.m-2) and layer thickness
544                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_b(ji,jk)
545                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
546
547               ELSE                               !!! Basal melting
548
549                  zEi               = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)
550
551                  zEw               = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0)
552
553                  zdE               = zEi - zEw              ! Specific enthalpy difference   (J/kg, <0)
554
555                  zfmdt             = - zq_bo(ji) / zdE  ! Mass flux x time step (kg/m2, >0)
556
557                  zdeltah(ji,jk)    = - zfmdt / rhoic        ! Gross thickness change
558
559                  zdeltah(ji,jk)    = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
560                 
561                  zq_bo(ji)         = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors
562
563                  dh_i_bott(ji)     = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt
564
565                  zfmdt             = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0
566
567                  zQm               = zfmdt * zEw         ! Heat exchanged with ocean
568
569                  ! Contribution to heat flux to the ocean [W.m-2], <0 
570                  hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
571
572                  ! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok)
573                  sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
574                 
575                  ! Total heat flux used in this process [W.m-2], >0 
576                  hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
577                 
578                  ! Contribution to mass flux
579                  wfx_bom_1d(ji) =  wfx_bom_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
580
581                  ! update heat content (J.m-2) and layer thickness
582                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_b(ji,jk)
583                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
584               ENDIF
585           
586            ENDIF
587         END DO ! ji
588      END DO ! jk
589
590      !------------------------------------------------------------------------------!
591      ! Excessive ablation in a 1-category model
592      !     in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15)
593      !------------------------------------------------------------------------------!
594      ! ??? keep ???
595      ! clem bug: I think this should be included above, so we would not have to
596      !           track heat/salt/mass fluxes backwards
597!      IF( jpl == 1 ) THEN
598!         DO ji = kideb, kiut
599!            IF(  zf_tt(ji)  >=  0._wp  ) THEN
600!               zdh            = MAX( hmelt , dh_i_bott(ji) )
601!               zdvres         = zdh - dh_i_bott(ji) ! >=0
602!               dh_i_bott(ji)  = zdh
603!
604!               ! excessive energy is sent to lateral ablation
605!               zinda = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_b(ji) - epsi20 ) )
606!               zq_1cat(ji) =  zinda * rhoic * lfus * at_i_b(ji) / MAX( 1._wp - at_i_b(ji) , epsi20 ) * zdvres ! J.m-2 >=0
607!
608!               ! correct salt and mass fluxes
609!               sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice ! this is only a raw approximation
610!               wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_b(ji) * zdvres * r1_rdtice
611!            ENDIF
612!         END DO
613!      ENDIF
614
615      !-------------------------------------------
616      ! Update temperature, energy
617      !-------------------------------------------
618      DO ji = kideb, kiut
619         ht_i_b(ji) =  MAX( 0._wp , ht_i_b(ji) + dh_i_bott(ji) )
620      END DO 
621
622      !-------------------------------------------
623      ! 5. What to do with remaining energy
624      !-------------------------------------------
625      ! If heat still available for melting and snow remains, then melt more snow
626      !-------------------------------------------
627      zdeltah(:,:) = 0._wp ! important
628      DO ji = kideb, kiut
629         zq_rema(ji)     = zq_su(ji) + zq_bo(ji) 
630!         zindh           = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_b(ji) ) )   ! =1 if snow
631!         zindq           = 1._wp - MAX( 0._wp, SIGN( 1._wp, - zq_s(ji) + epsi20 ) )
632!         zdeltah  (ji,1) = - zindh * zindq * zq_rema(ji) / MAX( zq_s(ji), epsi20 )
633!         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_b(ji) ) ) ! bound melting
634!         zdh_s_mel(ji)   = zdh_s_mel(ji) + zdeltah(ji,1)   
635!         dh_s_tot (ji)   = dh_s_tot(ji) + zdeltah(ji,1)
636!         ht_s_b   (ji)   = ht_s_b(ji)   + zdeltah(ji,1)
637!       
638!         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji)                ! update available heat (J.m-2)
639!         ! heat used to melt snow
640!         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_b(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0)
641!         ! Contribution to mass flux
642!         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdeltah(ji,1) * r1_rdtice
643!   
644         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
645         ! Remaining heat flux (W.m-2) is sent to the ocean heat budget
646         hfx_out(ii,ij)  = hfx_out(ii,ij) + ( zq_1cat(ji) + zq_rema(ji) * a_i_b(ji) ) * r1_rdtice
647
648         IF( ln_nicep .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
649      END DO
650     
651      !
652      !------------------------------------------------------------------------------|
653      !  6) Snow-Ice formation                                                       |
654      !------------------------------------------------------------------------------|
655      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
656      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
657      DO ji = kideb, kiut
658         !
659         dh_snowice(ji) = MAX(  0._wp , ( rhosn * ht_s_b(ji) + (rhoic-rau0) * ht_i_b(ji) ) / ( rhosn+rau0-rhoic )  )
660
661         ht_i_b(ji)     = ht_i_b(ji) + dh_snowice(ji)
662         ht_s_b(ji)     = ht_s_b(ji) - dh_snowice(ji)
663
664         ! Salinity of snow ice
665         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
666         zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_b(ji)
667
668         ! entrapment during snow ice formation
669         ! new salinity difference stored (to be used in limthd_ent.F90)
670         IF (  num_sal == 2  ) THEN
671            zswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_b(ji) - epsi10 ) )
672            ! salinity dif due to snow-ice formation
673            dsm_i_si_1d(ji) = ( zs_snic - sm_i_b(ji) ) * dh_snowice(ji) / MAX( ht_i_b(ji), epsi10 ) * zswitch     
674            ! salinity dif due to bottom growth
675            IF (  zf_tt(ji)  < 0._wp ) THEN
676               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( ht_i_b(ji), epsi10 ) * zswitch
677            ENDIF
678         ENDIF
679
680         ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
681         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
682         zfmdt          = ( rhosn - rhoic ) * MAX( dh_snowice(ji), 0._wp )    ! <0
683         zsstK          = sst_m(ii,ij) + rt0                               
684         zEw            = rcp * ( zsstK - rt0 )
685         zQm            = zfmdt * zEw 
686         
687         ! Contribution to heat flux
688         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice 
689
690         ! Contribution to salt flux
691         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_b(ji) * zfmdt * r1_rdtice 
692         
693         ! Contribution to mass flux
694         ! All snow is thrown in the ocean, and seawater is taken to replace the volume
695         wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice
696         wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_b(ji) * dh_snowice(ji) * rhosn * r1_rdtice
697
698         ! update heat content (J.m-2) and layer thickness
699         qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_b(ji,1) + zfmdt * zEw
700         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
701         
702         ! Total ablation (to debug)
703         IF( ht_i_b(ji) <= 0._wp )   a_i_b(ji) = 0._wp
704
705      END DO !ji
706
707      !
708      !-------------------------------------------
709      ! Update temperature, energy
710      !-------------------------------------------
711      !clem bug: we should take snow into account here
712      DO ji = kideb, kiut
713         zindh    =  1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_b(ji) ) ) 
714         t_su_b(ji) =  zindh * t_su_b(ji) + ( 1.0 - zindh ) * rtt
715      END DO  ! ji
716
717      DO jk = 1, nlay_s
718         DO ji = kideb,kiut
719            ! mask enthalpy
720            zinda        =  MAX(  0._wp , SIGN( 1._wp, - ht_s_b(ji) )  )
721            q_s_b(ji,jk) = ( 1.0 - zinda ) * q_s_b(ji,jk)
722            ! recalculate t_s_b from q_s_b
723            t_s_b(ji,jk) = rtt + ( 1._wp - zinda ) * ( - q_s_b(ji,jk) / ( rhosn * cpic ) + lfus / cpic )
724         END DO
725      END DO
726
727      CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema )
728      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
729      CALL wrk_dealloc( jpij, zintermelt )
730      CALL wrk_dealloc( jpij, jkmax, zdeltah, zh_i )
731      CALL wrk_dealloc( jpij, icount )
732      !
733      !
734   END SUBROUTINE lim_thd_dh
735   
736#else
737   !!----------------------------------------------------------------------
738   !!   Default option                               NO  LIM3 sea-ice model
739   !!----------------------------------------------------------------------
740CONTAINS
741   SUBROUTINE lim_thd_dh          ! Empty routine
742   END SUBROUTINE lim_thd_dh
743#endif
744
745   !!======================================================================
746END MODULE limthd_dh
Note: See TracBrowser for help on using the repository browser.