source: branches/2014/dev_r5134_UKMO4_CF_compliance/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 5350

Last change on this file since 5350 was 5350, checked in by hadcv, 6 years ago

Update to head of the trunk (r5344).

  • Property svn:keywords set to Id
File size: 35.8 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(:) ::   zqprec      ! energy of fallen snow                       (J.m-3)
89      REAL(wp), POINTER, DIMENSION(:) ::   zq_su       ! heat for surface ablation                   (J.m-2)
90      REAL(wp), POINTER, DIMENSION(:) ::   zq_bo       ! heat for bottom ablation                    (J.m-2)
91      REAL(wp), POINTER, DIMENSION(:) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2)
92      REAL(wp), POINTER, DIMENSION(:) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2)
93
94      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel   ! snow melt
95      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_pre   ! snow precipitation
96      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_sub   ! snow sublimation
97
98      REAL(wp), POINTER, DIMENSION(:,:) ::   zdeltah
99      REAL(wp), POINTER, DIMENSION(:,:) ::   zh_i      ! ice layer thickness
100      INTEGER , POINTER, DIMENSION(:,:) ::   icount    ! number of layers vanished by melting
101
102      REAL(wp), POINTER, DIMENSION(:) ::   zqh_i       ! total ice heat content  (J.m-2)
103      REAL(wp), POINTER, DIMENSION(:) ::   zqh_s       ! total snow heat content (J.m-2)
104      REAL(wp), POINTER, DIMENSION(:) ::   zq_s        ! total snow enthalpy     (J.m-3)
105
106      REAL(wp) :: zswitch_sal
107
108      ! Heat conservation
109      INTEGER  ::   num_iter_max
110
111      !!------------------------------------------------------------------
112
113      ! Discriminate between varying salinity (nn_icesal=2) and prescribed cases (other values)
114      SELECT CASE( nn_icesal )                       ! varying salinity or not
115         CASE( 1, 3, 4 ) ;   zswitch_sal = 0       ! prescribed salinity profile
116         CASE( 2 )       ;   zswitch_sal = 1       ! varying salinity profile
117      END SELECT
118
119      CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema )
120      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
121      CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i )
122      CALL wrk_alloc( jpij, nlay_i, icount )
123     
124      dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp
125      dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp   
126 
127      zqprec (:) = 0._wp ; zq_su  (:) = 0._wp ; zq_bo  (:) = 0._wp ; zf_tt  (:) = 0._wp
128      zq_rema(:) = 0._wp
129
130      zdh_s_pre(:) = 0._wp
131      zdh_s_mel(:) = 0._wp
132      zdh_s_sub(:) = 0._wp
133      zqh_s    (:) = 0._wp     
134      zqh_i    (:) = 0._wp   
135
136      zh_i      (:,:) = 0._wp       
137      zdeltah   (:,:) = 0._wp       
138      icount    (:,:) = 0
139
140      ! Initialize enthalpy at nlay_i+1
141      DO ji = kideb, kiut
142         q_i_1d(ji,nlay_i+1) = 0._wp
143      END DO
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) * r1_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         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) - rt0 ) )
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 jk = 1, nlay_s
190         DO ji = kideb, kiut
191            zqh_s(ji) =  zqh_s(ji) + q_s_1d(ji,jk) * ht_s_1d(ji) * r1_nlay_s
192         END DO
193      END DO
194      !
195      DO jk = 1, nlay_i
196         DO ji = kideb, kiut
197            zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i
198            zqh_i(ji)   = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk)
199         END DO
200      END DO
201      !
202      !------------------------------------------------------------------------------|
203      !  3) Surface ablation and sublimation                                         |
204      !------------------------------------------------------------------------------|
205      !
206      !-------------------------
207      ! 3.1 Snow precips / melt
208      !-------------------------
209      ! Snow accumulation in one thermodynamic time step
210      ! snowfall is partitionned between leads and ice
211      ! if snow fall was uniform, a fraction (1-at_i) would fall into leads
212      ! but because of the winds, more snow falls on leads than on sea ice
213      ! and a greater fraction (1-at_i)^beta of the total mass of snow
214      ! (beta < 1) falls in leads.
215      ! In reality, beta depends on wind speed,
216      ! and should decrease with increasing wind speed but here, it is
217      ! considered as a constant. an average value is 0.66
218      ! Martin Vancoppenolle, December 2006
219
220      zdeltah(:,:) = 0._wp
221      DO ji = kideb, kiut
222         !-----------
223         ! Snow fall
224         !-----------
225         ! thickness change
226         zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**rn_betas ) / at_i_1d(ji) 
227         zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice * r1_rhosn
228         ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K)
229         zqprec   (ji) = rhosn * ( cpic * ( rt0 - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )   
230         IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp
231         ! heat flux from snow precip (>0, W.m-2)
232         hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
233         ! mass flux, <0
234         wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice
235
236         !---------------------
237         ! Melt of falling snow
238         !---------------------
239         ! thickness change
240         rswitch        = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )
241         zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 )
242         zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting
243         ! heat used to melt snow (W.m-2, >0)
244         hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
245         ! snow melting only = water into the ocean (then without snow precip), >0
246         wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice   
247         ! updates available heat + precipitations after melting
248         zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )     
249         zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1)
250
251         ! update thickness
252         ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) )
253      END DO
254
255      ! If heat still available (zq_su > 0), then melt more snow
256      zdeltah(:,:) = 0._wp
257      DO jk = 1, nlay_s
258         DO ji = kideb, kiut
259            ! thickness change
260            rswitch          = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 
261            rswitch          = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) ) 
262            zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 )
263            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - ht_s_1d(ji) ) ! bound melting
264            zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)   
265            ! heat used to melt snow(W.m-2, >0)
266            hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice 
267            ! snow melting only = water into the ocean (then without snow precip)
268            wfx_snw_1d(ji)   = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
269            ! updates available heat + thickness
270            zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) )
271            ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) )
272         END DO
273      END DO
274
275      !----------------------
276      ! 3.2 Snow sublimation
277      !----------------------
278      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
279      ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean)
280      ! clem comment: ice should also sublimate
281      zdeltah(:,:) = 0._wp
282      IF( lk_cpl ) THEN
283         ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice)
284         zdh_s_sub(:)      =  0._wp 
285      ELSE
286         ! forced  mode: snow thickness change due to sublimation
287         DO ji = kideb, kiut
288            zdh_s_sub(ji)  =  MAX( - ht_s_1d(ji) , - qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice )
289            ! Heat flux by sublimation [W.m-2], < 0
290            !      sublimate first snow that had fallen, then pre-existing snow
291            zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
292            hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * q_s_1d(ji,1)  &
293               &                              ) * a_i_1d(ji) * r1_rdtice
294            ! Mass flux by sublimation
295            wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice
296            ! new snow thickness
297            ht_s_1d(ji)    =  MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) )
298            ! update precipitations after sublimation and correct sublimation
299            zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
300            zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1)
301         END DO
302      ENDIF
303
304      ! --- Update snow diags --- !
305      DO ji = kideb, kiut
306         dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
307      END DO
308
309      !-------------------------------------------
310      ! 3.3 Update temperature, energy
311      !-------------------------------------------
312      ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
313      zq_s(:) = 0._wp 
314      DO jk = 1, nlay_s
315         DO ji = kideb,kiut
316            rswitch       = MAX(  0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 )  )
317            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *                          &
318              &            ( (   zdh_s_pre(ji)             ) * zqprec(ji) +  &
319              &              (   ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
320            zq_s(ji)     =  zq_s(ji) + q_s_1d(ji,jk)
321         END DO
322      END DO
323
324      !--------------------------
325      ! 3.4 Surface ice ablation
326      !--------------------------
327      zdeltah(:,:) = 0._wp ! important
328      DO jk = 1, nlay_i
329         DO ji = kideb, kiut
330            ztmelts           = - tmut * s_i_1d(ji,jk) + rt0          ! Melting point of layer k [K]
331           
332            IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
333
334               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]       
335               zdE            = 0._wp                                 ! Specific enthalpy difference   (J/kg, <0)
336                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted)
337               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing     
338                                                                      ! this should normally not happen, but sometimes, heat diffusion leads to this
339               zfmdt          = - zdeltah(ji,jk) * rhoic              ! Mass flux x time step > 0
340                         
341               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt
342               
343               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0]
344
345               ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
346               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
347               
348               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
349               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
350               
351               ! Contribution to mass flux
352               wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
353
354            ELSE                                !!! Surface melting
355               
356               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]
357               zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0]
358               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0
359               
360               zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0]
361               
362               zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0]
363               
364               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]
365               
366               zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat
367               
368               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt
369               
370               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0]
371               
372               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0]
373               
374               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
375               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
376               
377               ! Contribution to heat flux [W.m-2], < 0
378               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
379               
380               ! Total heat flux used in this process [W.m-2], > 0 
381               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
382               
383               ! Contribution to mass flux
384               wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
385               
386            END IF
387            ! record which layers have disappeared (for bottom melting)
388            !    => icount=0 : no layer has vanished
389            !    => icount=5 : 5 layers have vanished
390            rswitch       = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) 
391            icount(ji,jk) = NINT( rswitch )
392            zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
393
394            ! update heat content (J.m-2) and layer thickness
395            qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
396            h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
397         END DO
398      END DO
399      ! update ice thickness
400      DO ji = kideb, kiut
401         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) )
402      END DO
403
404      !
405      !------------------------------------------------------------------------------!
406      ! 4) Basal growth / melt                                                       !
407      !------------------------------------------------------------------------------!
408      !
409      !------------------
410      ! 4.1 Basal growth
411      !------------------
412      ! Basal growth is driven by heat imbalance at the ice-ocean interface,
413      ! between the inner conductive flux  (fc_bo_i), from the open water heat flux
414      ! (fhld) and the turbulent ocean flux (fhtur).
415      ! fc_bo_i is positive downwards. fhtur and fhld are positive to the ice
416
417      ! If salinity varies in time, an iterative procedure is required, because
418      ! the involved quantities are inter-dependent.
419      ! Basal growth (dh_i_bott) depends upon new ice specific enthalpy (zEi),
420      ! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bott)
421      ! -> need for an iterative procedure, which converges quickly
422
423      num_iter_max = 1
424      IF( nn_icesal == 2 ) num_iter_max = 5
425
426      ! Iterative procedure
427      DO iter = 1, num_iter_max
428         DO ji = kideb, kiut
429            IF(  zf_tt(ji) < 0._wp  ) THEN
430
431               ! New bottom ice salinity (Cox & Weeks, JGR88 )
432               !--- zswi1  if dh/dt < 2.0e-8
433               !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
434               !--- zswi2  if dh/dt > 3.6e-7
435               zgrr               = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) )
436               zswi2              = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
437               zswi12             = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
438               zswi1              = 1. - zswi2 * zswi12
439               zfracs             = MIN ( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
440                  &               + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 )
441
442               ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
443
444               s_i_new(ji)        = zswitch_sal * zfracs * sss_m(ii,ij)  &  ! New ice salinity
445                                  + ( 1. - zswitch_sal ) * sm_i_1d(ji) 
446               ! New ice growth
447               ztmelts            = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K)
448
449               zt_i_new           = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
450               
451               zEi                = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)     
452                  &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   &
453                  &               + rcp  * ( ztmelts-rt0 )         
454
455               zEw                = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
456
457               zdE                = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
458
459               dh_i_bott(ji)      = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
460
461               q_i_1d(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0)
462               
463            ENDIF
464         END DO
465      END DO
466
467      ! Contribution to Energy and Salt Fluxes
468      DO ji = kideb, kiut
469         IF(  zf_tt(ji) < 0._wp  ) THEN
470            ! New ice growth
471                                   
472            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0)
473
474            ztmelts        = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K)
475           
476            zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
477           
478            zEi            = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)     
479               &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   &
480               &               + rcp  * ( ztmelts-rt0 )         
481           
482            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
483           
484            zdE            = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
485           
486            ! Contribution to heat flux to the ocean [W.m-2], >0 
487            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
488
489            ! Total heat flux used in this process [W.m-2], <0 
490            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
491           
492            ! Contribution to salt flux, <0
493            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice
494
495            ! Contribution to mass flux, <0
496            wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice
497
498            ! update heat content (J.m-2) and layer thickness
499            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)
500            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji)
501         ENDIF
502      END DO
503
504      !----------------
505      ! 4.2 Basal melt
506      !----------------
507      zdeltah(:,:) = 0._wp ! important
508      DO jk = nlay_i, 1, -1
509         DO ji = kideb, kiut
510            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting
511
512               ztmelts = - tmut * s_i_1d(ji,jk) + rt0  ! Melting point of layer jk (K)
513
514               IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
515
516                  zEi               = - q_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0)
517                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0)
518                                                                    ! set up at 0 since no energy is needed to melt water...(it is already melted)
519                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing     
520                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this
521
522                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk)
523
524                  zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0
525
526                  ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
527                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
528
529                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
530                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
531                                   
532                  ! Contribution to mass flux
533                  wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
534
535                  ! update heat content (J.m-2) and layer thickness
536                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
537                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
538
539               ELSE                               !!! Basal melting
540
541                  zEi             = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)
542                  zEw             = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0)
543                  zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0)
544
545                  zfmdt           = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0)
546
547                  zdeltah(ji,jk)  = - zfmdt * r1_rhoic         ! Gross thickness change
548
549                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
550                 
551                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors
552
553                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt
554
555                  zfmdt           = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0
556
557                  zQm             = zfmdt * zEw         ! Heat exchanged with ocean
558
559                  ! Contribution to heat flux to the ocean [W.m-2], <0 
560                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
561
562                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
563                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
564                 
565                  ! Total heat flux used in this process [W.m-2], >0 
566                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
567                 
568                  ! Contribution to mass flux
569                  wfx_bom_1d(ji)  =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
570
571                  ! update heat content (J.m-2) and layer thickness
572                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
573                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
574               ENDIF
575           
576            ENDIF
577         END DO
578      END DO
579
580      !-------------------------------------------
581      ! Update temperature, energy
582      !-------------------------------------------
583      DO ji = kideb, kiut
584         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_bott(ji) )
585      END DO 
586
587      !-------------------------------------------
588      ! 5. What to do with remaining energy
589      !-------------------------------------------
590      ! If heat still available for melting and snow remains, then melt more snow
591      !-------------------------------------------
592      zdeltah(:,:) = 0._wp ! important
593      DO ji = kideb, kiut
594         zq_rema(ji)     = zq_su(ji) + zq_bo(ji) 
595         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )   ! =1 if snow
596         rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,1) - epsi20 ) )
597         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 )
598         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting
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) * q_s_1d(ji,1)                ! 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) * q_s_1d(ji,1) * 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_icectl .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 ) * r1_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_sal.F90)
634         IF (  nn_icesal == 2  ) THEN
635            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )
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), epsi20 ) * 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), epsi20 ) * 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      END DO
667
668      !
669      !-------------------------------------------
670      ! Update temperature, energy
671      !-------------------------------------------
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, 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, zdeltah, zh_i )
693      CALL wrk_dealloc( jpij, nlay_i, 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.