source: NEMO/trunk/src/ICE/icethd_dh.F90 @ 10069

Last change on this file since 10069 was 10069, checked in by nicolasmartin, 2 years ago

Fix mistakes of previous commit on SVN keywords property

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