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/2019/dev_r11842_SI3-10_EAP/src/ICE – NEMO

source: NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icethd_dh.F90 @ 13662

Last change on this file since 13662 was 13662, checked in by clem, 3 years ago

update to almost r4.0.4

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