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.
icethd_dh.F90 in NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/ICE – NEMO

source: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/ICE/icethd_dh.F90 @ 14021

Last change on this file since 14021 was 14021, checked in by laurent, 3 years ago

Caught up with trunk rev 14020...

  • Property svn:keywords set to Id
File size: 35.9 KB
Line 
1MODULE icethd_dh
2   !!======================================================================
3   !!                       ***  MODULE icethd_dh ***
4   !!   seaice : thermodynamic growth and melt
5   !!======================================================================
6   !! History :       !  2003-05  (M. Vancoppenolle) Original code in 1D
7   !!                 !  2005-06  (M. Vancoppenolle) 3D version
8   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
9   !!----------------------------------------------------------------------
10#if defined key_si3
11   !!----------------------------------------------------------------------
12   !!   'key_si3'                                       SI3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   ice_thd_dh        : vertical sea-ice growth and melt
15   !!----------------------------------------------------------------------
16   USE dom_oce        ! ocean space and time domain
17   USE phycst         ! physical constants
18   USE ice            ! sea-ice: variables
19   USE ice1D          ! sea-ice: thermodynamics variables
20   USE icethd_sal     ! sea-ice: salinity profiles
21   USE icevar         ! for CALL ice_var_snwblow
22   !
23   USE in_out_manager ! I/O manager
24   USE lib_mpp        ! MPP library
25   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   ice_thd_dh        ! called by ice_thd
31
32   !!----------------------------------------------------------------------
33   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
34   !! $Id$
35   !! Software governed by the CeCILL license (see ./LICENSE)
36   !!----------------------------------------------------------------------
37CONTAINS
38
39   SUBROUTINE ice_thd_dh
40      !!------------------------------------------------------------------
41      !!                ***  ROUTINE ice_thd_dh  ***
42      !!
43      !! ** Purpose :   compute ice and snow thickness changes due to growth/melting
44      !!
45      !! ** Method  :   Ice/Snow surface melting arises from imbalance in surface fluxes
46      !!                Bottom accretion/ablation arises from flux budget
47      !!                Snow thickness can increase by precipitation and decrease by sublimation
48      !!                If snow load excesses Archmiede limit, snow-ice is formed by
49      !!                the flooding of sea-water in the snow
50      !!
51      !!                - Compute available flux of heat for surface ablation
52      !!                - Compute snow and sea ice enthalpies
53      !!                - Surface ablation and sublimation
54      !!                - Bottom accretion/ablation
55      !!                - Snow ice formation
56      !!
57      !! ** Note     :  h=max(0,h+dh) are often used to ensure positivity of h.
58      !!                very small negative values can occur otherwise (e.g. -1.e-20)
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  ::   ji, jk       ! dummy loop indices
66      INTEGER  ::   iter         ! local integer
67
68      REAL(wp) ::   ztmelts      ! local scalar
69      REAL(wp) ::   zdum
70      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
71      REAL(wp) ::   zswi1        ! switch for computation of bottom salinity
72      REAL(wp) ::   zswi12       ! switch for computation of bottom salinity
73      REAL(wp) ::   zswi2        ! switch for computation of bottom salinity
74      REAL(wp) ::   zgrr         ! bottom growth rate
75      REAL(wp) ::   zt_i_new     ! bottom formation temperature
76      REAL(wp) ::   z1_rho       ! 1/(rhos+rho0-rhoi)
77
78      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
79      REAL(wp) ::   zEi          ! specific enthalpy of sea ice (J/kg)
80      REAL(wp) ::   zEw          ! specific enthalpy of exchanged water (J/kg)
81      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg)
82      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean
83
84      REAL(wp), DIMENSION(jpij) ::   zq_top      ! heat for surface ablation                   (J.m-2)
85      REAL(wp), DIMENSION(jpij) ::   zq_bot      ! heat for bottom ablation                    (J.m-2)
86      REAL(wp), DIMENSION(jpij) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2)
87      REAL(wp), DIMENSION(jpij) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2)
88      REAL(wp), DIMENSION(jpij) ::   zevap_rema  ! remaining mass flux from sublimation        (kg.m-2)
89      REAL(wp), DIMENSION(jpij) ::   zdeltah
90      REAL(wp), DIMENSION(jpij) ::   zsnw        ! distribution of snow after wind blowing
91
92      INTEGER , DIMENSION(jpij,nlay_i)     ::   icount    ! number of layers vanishing by melting
93      REAL(wp), DIMENSION(jpij,0:nlay_i+1) ::   zh_i      ! ice layer thickness (m)
94      REAL(wp), DIMENSION(jpij,0:nlay_s  ) ::   zh_s      ! snw layer thickness (m)
95      REAL(wp), DIMENSION(jpij,0:nlay_s  ) ::   ze_s      ! snw layer enthalpy (J.m-3)
96
97      REAL(wp) ::   zswitch_sal
98
99      INTEGER  ::   num_iter_max      ! Heat conservation
100      !!------------------------------------------------------------------
101
102      ! Discriminate between time varying salinity and constant
103      SELECT CASE( nn_icesal )                  ! varying salinity or not
104         CASE( 1 , 3 )   ;   zswitch_sal = 0._wp   ! prescribed salinity profile
105         CASE( 2 )       ;   zswitch_sal = 1._wp   ! varying salinity profile
106      END SELECT
107
108      ! initialize ice layer thicknesses and enthalpies
109      eh_i_old(1:npti,0:nlay_i+1) = 0._wp
110      h_i_old (1:npti,0:nlay_i+1) = 0._wp
111      zh_i    (1:npti,0:nlay_i+1) = 0._wp
112      DO jk = 1, nlay_i
113         DO ji = 1, npti
114            eh_i_old(ji,jk) = h_i_1d(ji) * r1_nlay_i * e_i_1d(ji,jk)
115            h_i_old (ji,jk) = h_i_1d(ji) * r1_nlay_i
116            zh_i    (ji,jk) = h_i_1d(ji) * r1_nlay_i
117         END DO
118      END DO
119      !
120      ! initialize snw layer thicknesses and enthalpies
121      zh_s(1:npti,0) = 0._wp
122      ze_s(1:npti,0) = 0._wp
123      DO jk = 1, nlay_s
124         DO ji = 1, npti
125            zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s
126            ze_s(ji,jk) = e_s_1d(ji,jk)
127         END DO
128      END DO
129      !
130      !                       ! ============================================== !
131      !                       ! Available heat for surface and bottom ablation !
132      !                       ! ============================================== !
133      !
134      IF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN
135         !
136         DO ji = 1, npti
137            zq_top(ji)     = MAX( 0._wp, qml_ice_1d(ji) * rDt_ice )
138         END DO
139         !
140      ELSE
141         !
142         DO ji = 1, npti
143            zdum           = qns_ice_1d(ji) + qsr_ice_1d(ji) - qtr_ice_top_1d(ji) - qcn_ice_top_1d(ji)
144            qml_ice_1d(ji) = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
145            zq_top(ji)     = MAX( 0._wp, qml_ice_1d(ji) * rDt_ice )
146         END DO
147         !
148      ENDIF
149      !
150      DO ji = 1, npti
151         zf_tt(ji)         = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) + qtr_ice_bot_1d(ji) * frq_m_1d(ji)
152         zq_bot(ji)        = MAX( 0._wp, zf_tt(ji) * rDt_ice )
153      END DO
154
155      !                       ! ============ !
156      !                       !     Snow     !
157      !                       ! ============ !
158      !
159      ! Internal melting
160      ! ----------------
161      ! IF snow temperature is above freezing point, THEN snow melts (should not happen but sometimes it does)
162      DO jk = 1, nlay_s
163         DO ji = 1, npti
164            IF( t_s_1d(ji,jk) > rt0 ) THEN
165               hfx_res_1d    (ji) = hfx_res_1d    (ji) - ze_s(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice   ! heat flux to the ocean [W.m-2], < 0
166               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos        * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice   ! mass flux
167               ! updates
168               dh_s_mlt(ji)    =             dh_s_mlt(ji) - zh_s(ji,jk)
169               h_s_1d  (ji)    = MAX( 0._wp, h_s_1d  (ji) - zh_s(ji,jk) )
170               zh_s    (ji,jk) = 0._wp
171               ze_s    (ji,jk) = 0._wp
172            END IF
173         END DO
174      END DO
175
176      ! Snow precipitation
177      !-------------------
178      CALL ice_var_snwblow( 1._wp - at_i_1d(1:npti), zsnw(1:npti) )   ! snow distribution over ice after wind blowing
179
180      DO ji = 1, npti
181         IF( sprecip_1d(ji) > 0._wp ) THEN
182            zh_s(ji,0) = zsnw(ji) * sprecip_1d(ji) * rDt_ice * r1_rhos / at_i_1d(ji)   ! thickness of precip
183            ze_s(ji,0) = MAX( 0._wp, - qprec_ice_1d(ji) )                              ! enthalpy of the precip (>0, J.m-3)
184            !
185            hfx_spr_1d(ji) = hfx_spr_1d(ji) + ze_s(ji,0) * zh_s(ji,0) * a_i_1d(ji) * r1_Dt_ice   ! heat flux from snow precip (>0, W.m-2)
186            wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos       * zh_s(ji,0) * a_i_1d(ji) * r1_Dt_ice   ! mass flux, <0
187            !
188            ! update thickness
189            h_s_1d(ji) = h_s_1d(ji) + zh_s(ji,0)
190         ENDIF
191      END DO
192
193      ! Snow melting
194      ! ------------
195      ! If heat still available (zq_top > 0)
196      ! then all snw precip has been melted and we need to melt more snow
197      DO jk = 0, nlay_s
198         DO ji = 1, npti
199            IF( zh_s(ji,jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN
200               !
201               rswitch = MAX( 0._wp , SIGN( 1._wp , ze_s(ji,jk) - epsi20 ) )
202               zdum    = - rswitch * zq_top(ji) / MAX( ze_s(ji,jk), epsi20 )   ! thickness change
203               zdum    = MAX( zdum , - zh_s(ji,jk) )                           ! bound melting
204
205               hfx_snw_1d    (ji) = hfx_snw_1d    (ji) - ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice   ! heat used to melt snow(W.m-2, >0)
206               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos        * zdum * a_i_1d(ji) * r1_Dt_ice   ! snow melting only = water into the ocean
207
208               ! updates available heat + thickness
209               dh_s_mlt(ji)    =              dh_s_mlt(ji)    + zdum
210               zq_top  (ji)    = MAX( 0._wp , zq_top  (ji)    + zdum * ze_s(ji,jk) )
211               h_s_1d  (ji)    = MAX( 0._wp , h_s_1d  (ji)    + zdum )
212               zh_s    (ji,jk) = MAX( 0._wp , zh_s    (ji,jk) + zdum )
213!!$               IF( zh_s(ji,jk) == 0._wp )   ze_s(ji,jk) = 0._wp
214               !
215            ENDIF
216         END DO
217      END DO
218
219      ! Snow sublimation
220      !-----------------
221      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
222      !    comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean)
223      zdeltah   (1:npti) = 0._wp ! total snow thickness that sublimates, < 0
224      zevap_rema(1:npti) = 0._wp
225      DO ji = 1, npti
226         IF( evap_ice_1d(ji) > 0._wp ) THEN
227            zdeltah   (ji) = MAX( - evap_ice_1d(ji) * r1_rhos * rDt_ice, - h_s_1d(ji) )   ! amount of snw that sublimates, < 0
228            zevap_rema(ji) = MAX( 0._wp, evap_ice_1d(ji) * rDt_ice + zdeltah(ji) * rhos ) ! remaining evap in kg.m-2 (used for ice sublimation later on)
229         ENDIF
230      END DO
231
232      DO jk = 0, nlay_s
233         DO ji = 1, npti
234            zdum = MAX( -zh_s(ji,jk), zdeltah(ji) ) ! snow layer thickness that sublimates, < 0
235            !
236            hfx_sub_1d    (ji) = hfx_sub_1d    (ji) + ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice  ! Heat flux of snw that sublimates [W.m-2], < 0
237            wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos        * zdum * a_i_1d(ji) * r1_Dt_ice  ! Mass flux by sublimation
238
239            ! update thickness
240            h_s_1d(ji)    = MAX( 0._wp , h_s_1d(ji)    + zdum )
241            zh_s  (ji,jk) = MAX( 0._wp , zh_s  (ji,jk) + zdum )
242!!$            IF( zh_s(ji,jk) == 0._wp )   ze_s(ji,jk) = 0._wp
243
244            ! update sublimation left
245            zdeltah(ji) = MIN( zdeltah(ji) - zdum, 0._wp )
246         END DO
247      END DO
248
249      !
250      !                       ! ============ !
251      !                       !     Ice      !
252      !                       ! ============ !
253
254      ! Surface ice melting
255      !--------------------
256      DO jk = 1, nlay_i
257         DO ji = 1, npti
258            ztmelts = - rTmlt * sz_i_1d(ji,jk)   ! Melting point of layer k [C]
259
260            IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting
261
262               zEi            = - e_i_1d(ji,jk) * r1_rhoi             ! Specific enthalpy of layer k [J/kg, <0]
263               zdE            =   0._wp                               ! Specific enthalpy difference (J/kg, <0)
264               !                                                          set up at 0 since no energy is needed to melt water...(it is already melted)
265               zdum           = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing
266               !                                                          this should normally not happen, but sometimes, heat diffusion leads to this
267               zfmdt          = - zdum * rhoi                         ! Recompute mass flux [kg/m2, >0]
268               !
269               dh_i_itm(ji)   = dh_i_itm(ji) + zdum                   ! Cumulate internal melting
270               !
271               hfx_res_1d(ji) = hfx_res_1d(ji) + zEi  * zfmdt             * a_i_1d(ji) * r1_Dt_ice    ! Heat flux to the ocean [W.m-2], <0
272               !                                                                                          ice enthalpy zEi is "sent" to the ocean
273               wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * zdum              * a_i_1d(ji) * r1_Dt_ice    ! Mass flux
274               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice    ! Salt flux
275               !                                                                                          using s_i_1d and not sz_i_1d(jk) is ok
276            ELSE                                        !-- Surface melting
277
278               zEi            = - e_i_1d(ji,jk) * r1_rhoi             ! Specific enthalpy of layer k [J/kg, <0]
279               zEw            =    rcp * ztmelts                      ! Specific enthalpy of resulting meltwater [J/kg, <0]
280               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0
281
282               zfmdt          = - zq_top(ji) / zdE                    ! Mass flux to the ocean [kg/m2, >0]
283
284               zdum           = - zfmdt * r1_rhoi                     ! Melt of layer jk [m, <0]
285
286               zdum           = MIN( 0._wp , MAX( zdum , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0]
287
288               zq_top(ji)     = MAX( 0._wp , zq_top(ji) - zdum * rhoi * zdE ) ! update available heat
289
290               dh_i_sum(ji)   = dh_i_sum(ji) + zdum                   ! Cumulate surface melt
291
292               zfmdt          = - rhoi * zdum                         ! Recompute mass flux [kg/m2, >0]
293
294               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0]
295
296               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw  * zfmdt             * a_i_1d(ji) * r1_Dt_ice    ! Heat flux [W.m-2], < 0
297               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zdE  * zfmdt             * a_i_1d(ji) * r1_Dt_ice    ! Heat flux used in this process [W.m-2], > 0
298               wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * zdum              * a_i_1d(ji) * r1_Dt_ice    ! Mass flux
299               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice    ! Salt flux >0
300               !                                                                                          using s_i_1d and not sz_i_1d(jk) is ok)
301            END IF
302            ! update thickness
303            zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum )
304            h_i_1d(ji)  = MAX( 0._wp, h_i_1d(ji)  + zdum )
305            !
306            ! update heat content (J.m-2) and layer thickness
307            eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk)
308            h_i_old (ji,jk) = h_i_old (ji,jk) + zdum
309            !
310            !
311            ! Ice sublimation
312            ! ---------------
313            zdum               = MAX( - zh_i(ji,jk) , - zevap_rema(ji) * r1_rhoi )
314            !
315            hfx_sub_1d(ji)     = hfx_sub_1d(ji)     + e_i_1d(ji,jk) * zdum              * a_i_1d(ji) * r1_Dt_ice ! Heat flux [W.m-2], < 0
316            wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi          * zdum              * a_i_1d(ji) * r1_Dt_ice ! Mass flux > 0
317            sfx_sub_1d(ji)     = sfx_sub_1d(ji)     - rhoi          * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux >0
318            !                                                                                                      clem: flux is sent to the ocean for simplicity
319            !                                                                                                            but salt should remain in the ice except
320            !                                                                                                            if all ice is melted. => must be corrected
321            ! update remaining mass flux and thickness
322            zevap_rema(ji) = zevap_rema(ji) + zdum * rhoi
323            zh_i(ji,jk)    = MAX( 0._wp, zh_i(ji,jk) + zdum )
324            h_i_1d(ji)     = MAX( 0._wp, h_i_1d(ji)  + zdum )
325            dh_i_sub(ji)   = dh_i_sub(ji) + zdum
326
327            ! update heat content (J.m-2) and layer thickness
328            eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk)
329            h_i_old (ji,jk) = h_i_old (ji,jk) + zdum
330
331            ! record which layers have disappeared (for bottom melting)
332            !    => icount=0 : no layer has vanished
333            !    => icount=5 : 5 layers have vanished
334            rswitch       = MAX( 0._wp , SIGN( 1._wp , - zh_i(ji,jk) ) )
335            icount(ji,jk) = NINT( rswitch )
336
337         END DO
338      END DO
339
340      ! remaining "potential" evap is sent to ocean
341      DO ji = 1, npti
342         wfx_err_sub_1d(ji) = wfx_err_sub_1d(ji) - zevap_rema(ji) * a_i_1d(ji) * r1_Dt_ice  ! <=0 (net evap for the ocean in kg.m-2.s-1)
343      END DO
344
345
346      ! Ice Basal growth
347      !------------------
348      ! Basal growth is driven by heat imbalance at the ice-ocean interface,
349      ! between the inner conductive flux  (qcn_ice_bot), from the open water heat flux
350      ! (fhld) and the sensible ice-ocean flux (qsb_ice_bot).
351      ! qcn_ice_bot is positive downwards. qsb_ice_bot and fhld are positive to the ice
352
353      ! If salinity varies in time, an iterative procedure is required, because
354      ! the involved quantities are inter-dependent.
355      ! Basal growth (dh_i_bog) depends upon new ice specific enthalpy (zEi),
356      ! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bog)
357      ! -> need for an iterative procedure, which converges quickly
358
359      num_iter_max = 1
360      IF( nn_icesal == 2 )   num_iter_max = 5  ! salinity varying in time
361
362      DO ji = 1, npti
363         IF(  zf_tt(ji) < 0._wp  ) THEN
364            DO iter = 1, num_iter_max   ! iterations
365
366               ! New bottom ice salinity (Cox & Weeks, JGR88 )
367               !--- zswi1  if dh/dt < 2.0e-8
368               !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
369               !--- zswi2  if dh/dt > 3.6e-7
370               zgrr     = MIN( 1.0e-3, MAX ( dh_i_bog(ji) * r1_Dt_ice , epsi10 ) )
371               zswi2    = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
372               zswi12   = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
373               zswi1    = 1. - zswi2 * zswi12
374               zfracs   = MIN( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
375                  &          + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 )
376
377               s_i_new(ji)    = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji)  ! New ice salinity
378
379               ztmelts        = - rTmlt * s_i_new(ji)                                                  ! New ice melting point (C)
380
381               zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
382
383               zEi            = rcpi * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0)
384                  &             - rLfus * ( 1.0 - ztmelts / ( MIN( zt_i_new - rt0, -epsi10 ) ) ) + rcp * ztmelts
385
386               zEw            = rcp  * ( t_bo_1d(ji) - rt0 )                                           ! Specific enthalpy of seawater (J/kg, < 0)
387
388               zdE            = zEi - zEw                                                              ! Specific enthalpy difference (J/kg, <0)
389
390               dh_i_bog(ji)   = rDt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) )
391
392            END DO
393            ! Contribution to Energy and Salt Fluxes
394            zfmdt = - rhoi * dh_i_bog(ji)                                                              ! Mass flux x time step (kg/m2, < 0)
395
396            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw  * zfmdt                      * a_i_1d(ji) * r1_Dt_ice   ! Heat flux to the ocean [W.m-2], >0
397            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zdE  * zfmdt                      * a_i_1d(ji) * r1_Dt_ice   ! Heat flux used in this process [W.m-2], <0
398            wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * dh_i_bog(ji)               * a_i_1d(ji) * r1_Dt_ice   ! Mass flux, <0
399            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * dh_i_bog(ji) * s_i_new(ji) * a_i_1d(ji) * r1_Dt_ice   ! Salt flux, <0
400
401            ! update thickness
402            zh_i(ji,nlay_i+1) = zh_i(ji,nlay_i+1) + dh_i_bog(ji)
403            h_i_1d(ji)        = h_i_1d(ji)        + dh_i_bog(ji)
404
405            ! update heat content (J.m-2) and layer thickness
406            eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoi)
407            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bog(ji)
408
409         ENDIF
410
411      END DO
412
413      ! Ice Basal melt
414      !---------------
415      DO jk = nlay_i, 1, -1
416         DO ji = 1, npti
417            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting
418
419               ztmelts = - rTmlt * sz_i_1d(ji,jk)  ! Melting point of layer jk (C)
420
421               IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting
422
423                  zEi            = - e_i_1d(ji,jk) * r1_rhoi     ! Specific enthalpy of melting ice (J/kg, <0)
424                  zdE            = 0._wp                         ! Specific enthalpy difference   (J/kg, <0)
425                  !                                                  set up at 0 since no energy is needed to melt water...(it is already melted)
426                  zdum           = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing
427                  !                                                  this should normally not happen, but sometimes, heat diffusion leads to this
428                  dh_i_itm (ji)  = dh_i_itm(ji) + zdum
429                  !
430                  zfmdt          = - zdum * rhoi                 ! Mass flux x time step > 0
431                  !
432                  hfx_res_1d(ji) = hfx_res_1d(ji) + zEi  * zfmdt             * a_i_1d(ji) * r1_Dt_ice   ! Heat flux to the ocean [W.m-2], <0
433                  !                                                                                         ice enthalpy zEi is "sent" to the ocean
434                  wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * zdum              * a_i_1d(ji) * r1_Dt_ice   ! Mass flux
435                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice   ! Salt flux
436                  !                                                                                         using s_i_1d and not sz_i_1d(jk) is ok
437               ELSE                                        !-- Basal melting
438
439                  zEi            = - e_i_1d(ji,jk) * r1_rhoi                       ! Specific enthalpy of melting ice (J/kg, <0)
440                  zEw            = rcp * ztmelts                                   ! Specific enthalpy of meltwater (J/kg, <0)
441                  zdE            = zEi - zEw                                       ! Specific enthalpy difference   (J/kg, <0)
442
443                  zfmdt          = - zq_bot(ji) / zdE                              ! Mass flux x time step (kg/m2, >0)
444
445                  zdum           = - zfmdt * r1_rhoi                               ! Gross thickness change
446
447                  zdum           = MIN( 0._wp , MAX( zdum, - zh_i(ji,jk) ) )       ! bound thickness change
448
449                  zq_bot(ji)     = MAX( 0._wp , zq_bot(ji) - zdum * rhoi * zdE )   ! update available heat. MAX is necessary for roundup errors
450
451                  dh_i_bom(ji)   = dh_i_bom(ji) + zdum                             ! Update basal melt
452
453                  zfmdt          = - zdum * rhoi                                   ! Mass flux x time step > 0
454
455                  zQm            = zfmdt * zEw                                     ! Heat exchanged with ocean
456
457                  hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw  * zfmdt             * a_i_1d(ji) * r1_Dt_ice   ! Heat flux to the ocean [W.m-2], <0
458                  hfx_bom_1d(ji) = hfx_bom_1d(ji) - zdE  * zfmdt             * a_i_1d(ji) * r1_Dt_ice   ! Heat used in this process [W.m-2], >0
459                  wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoi * zdum              * a_i_1d(ji) * r1_Dt_ice   ! Mass flux
460                  sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice   ! Salt flux
461                  !                                                                                         using s_i_1d and not sz_i_1d(jk) is ok
462               ENDIF
463               ! update thickness
464               zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum )
465               h_i_1d(ji)  = MAX( 0._wp, h_i_1d(ji)  + zdum )
466               !
467               ! update heat content (J.m-2) and layer thickness
468               eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk)
469               h_i_old (ji,jk) = h_i_old (ji,jk) + zdum
470            ENDIF
471         END DO
472      END DO
473
474      ! Remove snow if ice has melted entirely
475      ! --------------------------------------
476      DO jk = 0, nlay_s
477         DO ji = 1,npti
478            IF( h_i_1d(ji) == 0._wp ) THEN
479               ! mass & energy loss to the ocean
480               hfx_res_1d(ji) = hfx_res_1d(ji) - ze_s(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice  ! heat flux to the ocean [W.m-2], < 0
481               wfx_res_1d(ji) = wfx_res_1d(ji) + rhos        * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice  ! mass flux
482
483               ! update thickness and energy
484               h_s_1d(ji)    = 0._wp
485               ze_s  (ji,jk) = 0._wp
486               zh_s  (ji,jk) = 0._wp
487            ENDIF
488         END DO
489      END DO
490
491      ! Snow load on ice
492      ! -----------------
493      ! When snow load exceeds Archimede's limit and sst is positive,
494      ! snow-ice formation (next bloc) can lead to negative ice enthalpy.
495      ! Therefore we consider here that this excess of snow falls into the ocean
496      zdeltah(1:npti) = h_s_1d(1:npti) + h_i_1d(1:npti) * (rhoi-rho0) * r1_rhos
497      DO jk = 0, nlay_s
498         DO ji = 1, npti
499            IF( zdeltah(ji) > 0._wp .AND. sst_1d(ji) > 0._wp ) THEN
500               ! snow layer thickness that falls into the ocean
501               zdum = MIN( zdeltah(ji) , zh_s(ji,jk) )
502               ! mass & energy loss to the ocean
503               hfx_res_1d(ji) = hfx_res_1d(ji) - ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice  ! heat flux to the ocean [W.m-2], < 0
504               wfx_res_1d(ji) = wfx_res_1d(ji) + rhos        * zdum * a_i_1d(ji) * r1_Dt_ice  ! mass flux
505               ! update thickness and energy
506               h_s_1d(ji)    = MAX( 0._wp, h_s_1d(ji)  - zdum )
507               zh_s  (ji,jk) = MAX( 0._wp, zh_s(ji,jk) - zdum )
508               ! update snow thickness that still has to fall
509               zdeltah(ji)   = MAX( 0._wp, zdeltah(ji) - zdum )
510            ENDIF
511         END DO
512      END DO
513
514      ! Snow-Ice formation
515      ! ------------------
516      ! When snow load exceeds Archimede's limit, snow-ice interface goes down under sea-level,
517      ! flooding of seawater transforms snow into ice. Thickness that is transformed is dh_snowice (positive for the ice)
518      z1_rho = 1._wp / ( rhos+rho0-rhoi )
519      zdeltah(1:npti) = 0._wp
520      DO ji = 1, npti
521         !
522         dh_snowice(ji) = MAX( 0._wp , ( rhos * h_s_1d(ji) + (rhoi-rho0) * h_i_1d(ji) ) * z1_rho )
523
524         h_i_1d(ji)    = h_i_1d(ji) + dh_snowice(ji)
525         h_s_1d(ji)    = h_s_1d(ji) - dh_snowice(ji)
526
527         ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
528         zfmdt          = ( rhos - rhoi ) * dh_snowice(ji)    ! <0
529         zEw            = rcp * sst_1d(ji)
530         zQm            = zfmdt * zEw
531
532         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw        * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux
533         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Salt flux
534
535         ! Case constant salinity in time: virtual salt flux to keep salinity constant
536         IF( nn_icesal /= 2 )  THEN
537            sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d(ji) * zfmdt                 * a_i_1d(ji) * r1_Dt_ice  &  ! put back sss_m     into the ocean
538               &                            - s_i_1d(ji) * dh_snowice(ji) * rhoi * a_i_1d(ji) * r1_Dt_ice     ! and get  rn_icesal from the ocean
539         ENDIF
540
541         ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume
542         wfx_sni_1d    (ji) = wfx_sni_1d    (ji) - dh_snowice(ji) * rhoi * a_i_1d(ji) * r1_Dt_ice
543         wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + dh_snowice(ji) * rhos * a_i_1d(ji) * r1_Dt_ice
544
545         ! update thickness
546         zh_i(ji,0)  = zh_i(ji,0) + dh_snowice(ji)
547         zdeltah(ji) =              dh_snowice(ji)
548
549         ! update heat content (J.m-2) and layer thickness
550         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
551         eh_i_old(ji,0) = eh_i_old(ji,0) + zfmdt * zEw           ! 1st part (sea water enthalpy)
552
553      END DO
554      !
555      DO jk = nlay_s, 0, -1   ! flooding of snow starts from the base
556         DO ji = 1, npti
557            zdum           = MIN( zdeltah(ji), zh_s(ji,jk) )     ! amount of snw that floods, > 0
558            zh_s(ji,jk)    = MAX( 0._wp, zh_s(ji,jk) - zdum )    ! remove some snow thickness
559            eh_i_old(ji,0) = eh_i_old(ji,0) + zdum * ze_s(ji,jk) ! 2nd part (snow enthalpy)
560            ! update dh_snowice
561            zdeltah(ji)    = MAX( 0._wp, zdeltah(ji) - zdum )
562         END DO
563      END DO
564      !
565      !
566!!$      ! --- Update snow diags --- !
567!!$      !!clem: this is wrong. dh_s_tot is not used anyway
568!!$      DO ji = 1, npti
569!!$         dh_s_tot(ji) = dh_s_tot(ji) + dh_s_mlt(ji) + zdeltah(ji) + zdh_s_sub(ji) - dh_snowice(ji)
570!!$      END DO
571      !
572      !
573      ! Remapping of snw enthalpy on a regular grid
574      !--------------------------------------------
575      CALL snw_ent( zh_s, ze_s, e_s_1d )
576
577      ! recalculate t_s_1d from e_s_1d
578      DO jk = 1, nlay_s
579         DO ji = 1,npti
580            IF( h_s_1d(ji) > 0._wp ) THEN
581               t_s_1d(ji,jk) = rt0 + ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi )
582            ELSE
583               t_s_1d(ji,jk) = rt0
584            ENDIF
585         END DO
586      END DO
587
588      ! Note: remapping of ice enthalpy is done in icethd.F90
589
590      ! --- ensure that a_i = 0 & h_s = 0 where h_i = 0 ---
591      WHERE( h_i_1d(1:npti) == 0._wp )
592         a_i_1d (1:npti) = 0._wp
593         h_s_1d (1:npti) = 0._wp
594         t_su_1d(1:npti) = rt0
595      END WHERE
596
597   END SUBROUTINE ice_thd_dh
598
599   SUBROUTINE snw_ent( ph_old, pe_old, pe_new )
600      !!-------------------------------------------------------------------
601      !!               ***   ROUTINE snw_ent  ***
602      !!
603      !! ** Purpose :
604      !!           This routine computes new vertical grids in the snow,
605      !!           and consistently redistributes temperatures.
606      !!           Redistribution is made so as to ensure to energy conservation
607      !!
608      !!
609      !! ** Method  : linear conservative remapping
610      !!
611      !! ** Steps : 1) cumulative integrals of old enthalpies/thicknesses
612      !!            2) linear remapping on the new layers
613      !!
614      !! ------------ cum0(0)                        ------------- cum1(0)
615      !!                                    NEW      -------------
616      !! ------------ cum0(1)               ==>      -------------
617      !!     ...                                     -------------
618      !! ------------                                -------------
619      !! ------------ cum0(nlay_s+1)                 ------------- cum1(nlay_s)
620      !!
621      !!
622      !! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005
623      !!-------------------------------------------------------------------
624      REAL(wp), DIMENSION(jpij,0:nlay_s), INTENT(in   ) ::   ph_old             ! old thicknesses (m)
625      REAL(wp), DIMENSION(jpij,0:nlay_s), INTENT(in   ) ::   pe_old             ! old enthlapies (J.m-3)
626      REAL(wp), DIMENSION(jpij,1:nlay_s), INTENT(inout) ::   pe_new             ! new enthlapies (J.m-3, remapped)
627      !
628      INTEGER  :: ji         !  dummy loop indices
629      INTEGER  :: jk0, jk1   !  old/new layer indices
630      !
631      REAL(wp), DIMENSION(jpij,0:nlay_s+1) ::   zeh_cum0, zh_cum0   ! old cumulative enthlapies and layers interfaces
632      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeh_cum1, zh_cum1   ! new cumulative enthlapies and layers interfaces
633      REAL(wp), DIMENSION(jpij)            ::   zhnew               ! new layers thicknesses
634      !!-------------------------------------------------------------------
635
636      !--------------------------------------------------------------------------
637      !  1) Cumulative integral of old enthalpy * thickness and layers interfaces
638      !--------------------------------------------------------------------------
639      zeh_cum0(1:npti,0) = 0._wp
640      zh_cum0 (1:npti,0) = 0._wp
641      DO jk0 = 1, nlay_s+1
642         DO ji = 1, npti
643            zeh_cum0(ji,jk0) = zeh_cum0(ji,jk0-1) + pe_old(ji,jk0-1) * ph_old(ji,jk0-1)
644            zh_cum0 (ji,jk0) = zh_cum0 (ji,jk0-1) + ph_old(ji,jk0-1)
645         END DO
646      END DO
647
648      !------------------------------------
649      !  2) Interpolation on the new layers
650      !------------------------------------
651      ! new layer thickesses
652      DO ji = 1, npti
653         zhnew(ji) = SUM( ph_old(ji,0:nlay_s) ) * r1_nlay_s
654      END DO
655
656      ! new layers interfaces
657      zh_cum1(1:npti,0) = 0._wp
658      DO jk1 = 1, nlay_s
659         DO ji = 1, npti
660            zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + zhnew(ji)
661         END DO
662      END DO
663
664      zeh_cum1(1:npti,0:nlay_s) = 0._wp
665      ! new cumulative q*h => linear interpolation
666      DO jk0 = 1, nlay_s+1
667         DO jk1 = 1, nlay_s-1
668            DO ji = 1, npti
669               IF( zh_cum1(ji,jk1) <= zh_cum0(ji,jk0) .AND. zh_cum1(ji,jk1) > zh_cum0(ji,jk0-1) ) THEN
670                  zeh_cum1(ji,jk1) = ( zeh_cum0(ji,jk0-1) * ( zh_cum0(ji,jk0) - zh_cum1(ji,jk1  ) ) +  &
671                     &                 zeh_cum0(ji,jk0  ) * ( zh_cum1(ji,jk1) - zh_cum0(ji,jk0-1) ) )  &
672                     &             / ( zh_cum0(ji,jk0) - zh_cum0(ji,jk0-1) )
673               ENDIF
674            END DO
675         END DO
676      END DO
677      ! to ensure that total heat content is strictly conserved, set:
678      zeh_cum1(1:npti,nlay_s) = zeh_cum0(1:npti,nlay_s+1)
679
680      ! new enthalpies
681      DO jk1 = 1, nlay_s
682         DO ji = 1, npti
683            rswitch      = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) )
684            pe_new(ji,jk1) = rswitch * ( zeh_cum1(ji,jk1) - zeh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 )
685         END DO
686      END DO
687
688   END SUBROUTINE snw_ent
689
690
691#else
692   !!----------------------------------------------------------------------
693   !!   Default option                                NO SI3 sea-ice model
694   !!----------------------------------------------------------------------
695#endif
696
697   !!======================================================================
698END MODULE icethd_dh
Note: See TracBrowser for help on using the repository browser.