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

Last change on this file since 9767 was 9767, checked in by clem, 2 years ago

remove cp_ice_msh option in NEMO because the ice model is in C-grid and not anymore in B-grid (LIM2 type)

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