source: NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd_dh.F90 @ 12484

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