source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 5134

Last change on this file since 5134 was 5134, checked in by clem, 6 years ago

LIM3: changes to constrain ice thickness after advection. Ice volume and concentration are advected while ice thickness is deduced. This could result in overly high thicknesses associated with very low concentrations (in the 5th category)

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