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/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 5048

Last change on this file since 5048 was 5048, checked in by vancop, 9 years ago

new itd boundaries, namelist changes, mono-category and comments

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