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 on Ticket #2257 – Attachment – NEMO

Ticket #2257: icethd_dh.F90

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