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

Last change on this file since 4990 was 4990, checked in by timgraham, 6 years ago

Merged branches/2014/dev_MERGE_2014 back onto the trunk as follows:

In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
1 conflict in LIM_SRC_3/limdiahsb.F90
Resolved by keeping the version from dev_MERGE_2014 branch
and commited at r4989

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge —reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2014/dev_MERGE_2014
to merge the branch into the trunk - no conflicts at this stage.

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