source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90 @ 8486

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

changes in style - part1 - (now the code looks better txs to Gurvan's comments)

File size: 35.3 KB
Line 
1MODULE icethd_dh
2   !!======================================================================
3   !!                       ***  MODULE icethd_dh ***
4   !!  LIM-3 :   thermodynamic growth and decay of the ice
5   !!======================================================================
6   !! History :  LIM  ! 2003-05 (M. Vancoppenolle) Original code in 1D
7   !!                 ! 2005-06 (M. Vancoppenolle) 3D version
8   !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in wfx_snw & wfx_ice
9   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
10   !!            3.5  ! 2012-10 (G. Madec & co) salt flux + bug fixes
11   !!----------------------------------------------------------------------
12#if defined key_lim3
13   !!----------------------------------------------------------------------
14   !!   'key_lim3'                                       LIM3 sea-ice model
15   !!----------------------------------------------------------------------
16   !!   ice_thd_dh    : vertical accr./abl. and lateral ablation of sea ice
17   !!----------------------------------------------------------------------
18   USE par_oce        ! ocean parameters
19   USE phycst         ! physical constants (OCE directory)
20   USE ice            ! LIM variables
21   USE ice1D          ! LIM thermodynamics
22   !
23   USE in_out_manager ! I/O manager
24   USE lib_mpp        ! MPP library
25   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
26   
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   ice_thd_dh      ! called by ice_thd
31   PUBLIC   ice_thd_snwblow ! called in sbcblk/sbcclio/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 (2017)
39   !! $Id: icethd_dh.F90 8420 2017-08-08 12:18:46Z clem $
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE ice_thd_dh
45      !!------------------------------------------------------------------
46      !!                ***  ROUTINE ice_thd_dh  ***
47      !!
48      !! ** Purpose :   determines variations of ice and snow thicknesses.
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      !!                 1) Compute available flux of heat for surface ablation
57      !!                 2) Compute snow and sea ice enthalpies
58      !!                 3) Surface ablation and sublimation
59      !!                 4) Bottom accretion/ablation
60      !!                 5) Case of Total ablation
61      !!                 6) Snow ice formation
62      !!
63      !! References : Bitz and Lipscomb, 1999, J. Geophys. Res.
64      !!              Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
65      !!              Vancoppenolle, Fichefet and Bitz, 2005, Geophys. Res. Let.
66      !!              Vancoppenolle et al.,2009, Ocean Modelling
67      !!------------------------------------------------------------------
68      INTEGER  ::   ji, jk       ! dummy loop indices
69      INTEGER  ::   iter         ! local integer
70
71      REAL(wp) ::   ztmelts      ! local scalar
72      REAL(wp) ::   zdum       
73      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
74      REAL(wp) ::   zswi1        ! switch for computation of bottom salinity
75      REAL(wp) ::   zswi12       ! switch for computation of bottom salinity
76      REAL(wp) ::   zswi2        ! switch for computation of bottom salinity
77      REAL(wp) ::   zgrr         ! bottom growth rate
78      REAL(wp) ::   zt_i_new     ! bottom formation temperature
79
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_i) ::   zdeltah
98      REAL(wp), DIMENSION(jpij,nlay_i) ::   zh_i      ! ice layer thickness
99      INTEGER , DIMENSION(jpij,nlay_i) ::   icount    ! number of layers vanished by melting
100
101      REAL(wp), DIMENSION(jpij) ::   zeh_i       ! total ice heat content  (J.m-2)
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 varying salinity (nn_icesal=2) and prescribed cases (other values)
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      DO ji = 1, nidx
116         icount (ji,:) = 0
117         zdh_s_mel(ji) = 0._wp
118         e_i_1d(ji,nlay_i+1) = 0._wp ! Initialize enthalpy at nlay_i+1
119      END DO
120
121      ! initialize layer thicknesses and enthalpies
122      h_i_old (1:nidx,0:nlay_i+1) = 0._wp
123      eh_i_old(1:nidx,0:nlay_i+1) = 0._wp
124      DO jk = 1, nlay_i
125         DO ji = 1, nidx
126            h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i
127            eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk)
128         END DO
129      END DO
130      !
131      !------------------------------------------------------------------------------!
132      !  1) Calculate available heat for surface and bottom ablation                 !
133      !------------------------------------------------------------------------------!
134      !
135      DO ji = 1, nidx
136         zdum       = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
137         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 
138
139         zq_su (ji) = MAX( 0._wp, zdum      * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
140         zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
141      END DO
142
143      !------------------------------------------------------------------------------!
144      ! If snow temperature is above freezing point, then snow melts
145      ! (should not happen but sometimes it does)
146      !------------------------------------------------------------------------------!
147      DO ji = 1, nidx
148         IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting
149            ! Contribution to heat flux to the ocean [W.m-2], < 0 
150            hfx_res_1d(ji) = hfx_res_1d(ji) + e_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
151            ! Contribution to mass flux
152            wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
153            ! updates
154            ht_s_1d(ji)   = 0._wp
155            e_s_1d (ji,1) = 0._wp
156            t_s_1d (ji,1) = rt0
157         END IF
158      END DO
159
160      !------------------------------------------------------------!
161      !  2) Computing layer thicknesses and enthalpies.            !
162      !------------------------------------------------------------!
163      DO jk = 1, nlay_i
164         DO ji = 1, nidx
165            zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i
166            zeh_i(ji)   = zeh_i(ji) + e_i_1d(ji,jk) * zh_i(ji,jk)
167         END DO
168      END DO
169
170      !------------------------------------------------------------------------------|
171      !  3) Surface ablation and sublimation                                         |
172      !------------------------------------------------------------------------------|
173      !
174      !-------------------------
175      ! 3.1 Snow precips / melt
176      !-------------------------
177      ! Snow accumulation in one thermodynamic time step
178      ! snowfall is partitionned between leads and ice
179      ! if snow fall was uniform, a fraction (1-at_i) would fall into leads
180      ! but because of the winds, more snow falls on leads than on sea ice
181      ! and a greater fraction (1-at_i)^beta of the total mass of snow
182      ! (beta < 1) falls in leads.
183      ! In reality, beta depends on wind speed,
184      ! and should decrease with increasing wind speed but here, it is
185      ! considered as a constant. an average value is 0.66
186      ! Martin Vancoppenolle, December 2006
187
188      CALL ice_thd_snwblow( 1. - at_i_1d(1:nidx), zsnw(1:nidx) ) ! snow distribution over ice after wind blowing
189
190      zdeltah(1:nidx,:) = 0._wp
191      DO ji = 1, nidx
192         !-----------
193         ! Snow fall
194         !-----------
195         ! thickness change
196         zdh_s_pre(ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji)
197         ! enthalpy of the precip (>0, J.m-3)
198         zqprec   (ji) = - qprec_ice_1d(ji)   
199         IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp
200         ! heat flux from snow precip (>0, W.m-2)
201         hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
202         ! mass flux, <0
203         wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice
204
205         !---------------------
206         ! Melt of falling snow
207         !---------------------
208         ! thickness change
209         rswitch        = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )
210         zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 )
211         zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting
212         ! heat used to melt snow (W.m-2, >0)
213         hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
214         ! snow melting only = water into the ocean (then without snow precip), >0
215         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
216         ! updates available heat + precipitations after melting
217         zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )     
218         zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1)
219
220         ! update thickness
221         ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) )
222      END DO
223
224      ! If heat still available (zq_su > 0), then melt more snow
225      zdeltah(1:nidx,:) = 0._wp
226      DO jk = 1, nlay_s
227         DO ji = 1, nidx
228            ! thickness change
229            rswitch          = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 
230            rswitch          = rswitch * ( MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) ) ) 
231            zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( e_s_1d(ji,jk), epsi20 )
232            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - ht_s_1d(ji) ) ! bound melting
233            zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)   
234            ! heat used to melt snow(W.m-2, >0)
235            hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d(ji,jk) * r1_rdtice 
236            ! snow melting only = water into the ocean (then without snow precip)
237            wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
238            ! updates available heat + thickness
239            zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
240            ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) )
241         END DO
242      END DO
243
244      !------------------------------
245      ! 3.2 Sublimation (part1: snow)
246      !------------------------------
247      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
248      ! clem comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean)
249      zdeltah(1:nidx,:) = 0._wp
250      DO ji = 1, nidx
251         zdh_s_sub(ji)  = MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice )
252         ! remaining evap in kg.m-2 (used for ice melting later on)
253         zevap_rema(ji)  = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhosn
254         ! Heat flux by sublimation [W.m-2], < 0 (sublimate first snow that had fallen, then pre-existing snow)
255         zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
256         hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1)  &
257            &                              ) * a_i_1d(ji) * r1_rdtice
258         ! Mass flux by sublimation
259         wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice
260
261         ! new snow thickness
262         ht_s_1d(ji)    =  MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) )
263         ! update precipitations after sublimation and correct sublimation
264         zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
265         zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1)
266      END DO
267     
268      ! --- Update snow diags --- !
269      DO ji = 1, nidx
270         dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
271      END DO
272
273      !-------------------------------------------
274      ! 3.3 Update temperature, energy
275      !-------------------------------------------
276      ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
277      DO jk = 1, nlay_s
278         DO ji = 1,nidx
279            rswitch       = MAX( 0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 ) )
280            e_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *           &
281              &            ( ( zdh_s_pre(ji)               ) * zqprec(ji) +  &
282              &              ( ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
283         END DO
284      END DO
285
286      !--------------------------
287      ! 3.4 Surface ice ablation
288      !--------------------------
289      zdeltah(1:nidx,:) = 0._wp ! important
290      DO jk = 1, nlay_i
291         DO ji = 1, nidx
292            ztmelts           = - tmut * s_i_1d(ji,jk) + rt0          ! Melting point of layer k [K]
293           
294            IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
295
296               zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]       
297               zdE            = 0._wp                                 ! Specific enthalpy difference   (J/kg, <0)
298                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted)
299               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing     
300                                                                      ! this should normally not happen, but sometimes, heat diffusion leads to this
301               zfmdt          = - zdeltah(ji,jk) * rhoic              ! Mass flux x time step > 0
302                         
303               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt
304               
305               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0]
306
307               ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
308               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
309               
310               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
311               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
312               
313               ! Contribution to mass flux
314               wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
315
316            ELSE                                !!! Surface melting
317               
318               zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]
319               zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0]
320               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0
321               
322               zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0]
323               
324               zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0]
325               
326               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]
327               
328               zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat
329               
330               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt
331               
332               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0]
333               
334               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0]
335               
336               ! Contribution to salt flux >0 (clem: using sm_i_1d and not s_i_1d(jk) is ok)
337               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
338               
339               ! Contribution to heat flux [W.m-2], < 0
340               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
341               
342               ! Total heat flux used in this process [W.m-2], > 0 
343               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
344               
345               ! Contribution to mass flux
346               wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
347               
348            END IF
349            ! ----------------------
350            ! Sublimation part2: ice
351            ! ----------------------
352            zdum      = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoic )
353            zdeltah(ji,jk) = zdeltah(ji,jk) + zdum
354            dh_i_sub(ji)  = dh_i_sub(ji) + zdum
355            ! Salt flux > 0 (clem2016: flux is sent to the ocean for simplicity but salt should remain in the ice except if all ice is melted.
356            !                          It must be corrected at some point)
357            sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * sm_i_1d(ji) * r1_rdtice
358            ! Heat flux [W.m-2], < 0
359            hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice
360            ! Mass flux > 0
361            wfx_ice_sub_1d(ji) =  wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice
362
363            ! update remaining mass flux
364            zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoic
365           
366            ! record which layers have disappeared (for bottom melting)
367            !    => icount=0 : no layer has vanished
368            !    => icount=5 : 5 layers have vanished
369            rswitch       = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) 
370            icount(ji,jk) = NINT( rswitch )
371            zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
372                       
373            ! update heat content (J.m-2) and layer thickness
374            eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
375            h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
376         END DO
377      END DO
378      ! update ice thickness
379      DO ji = 1, nidx
380         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) + dh_i_sub(ji) )
381      END DO
382
383      ! remaining "potential" evap is sent to ocean
384      DO ji = 1, nidx
385         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)
386      END DO
387
388
389      !------------------------------------------------------------------------------!
390      ! 4) Basal growth / melt                                                       !
391      !------------------------------------------------------------------------------!
392      !
393      !------------------
394      ! 4.1 Basal growth
395      !------------------
396      ! Basal growth is driven by heat imbalance at the ice-ocean interface,
397      ! between the inner conductive flux  (fc_bo_i), from the open water heat flux
398      ! (fhld) and the turbulent ocean flux (fhtur).
399      ! fc_bo_i is positive downwards. fhtur and fhld are positive to the ice
400
401      ! If salinity varies in time, an iterative procedure is required, because
402      ! the involved quantities are inter-dependent.
403      ! Basal growth (dh_i_bott) depends upon new ice specific enthalpy (zEi),
404      ! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bott)
405      ! -> need for an iterative procedure, which converges quickly
406
407      num_iter_max = 1
408      IF( nn_icesal == 2 ) num_iter_max = 5
409
410      ! Iterative procedure
411      DO ji = 1, nidx
412         IF(  zf_tt(ji) < 0._wp  ) THEN
413            DO iter = 1, num_iter_max
414
415               ! New bottom ice salinity (Cox & Weeks, JGR88 )
416               !--- zswi1  if dh/dt < 2.0e-8
417               !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
418               !--- zswi2  if dh/dt > 3.6e-7
419               zgrr               = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) )
420               zswi2              = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
421               zswi12             = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
422               zswi1              = 1. - zswi2 * zswi12
423               zfracs             = MIN ( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
424                  &               + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 )
425
426               s_i_new(ji)        = zswitch_sal * zfracs * sss_1d(ji)  &  ! New ice salinity
427                                  + ( 1. - zswitch_sal ) * sm_i_1d(ji) 
428               ! New ice growth
429               ztmelts            = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K)
430
431               zt_i_new           = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
432               
433               zEi                = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)     
434                  &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   &
435                  &               + rcp  * ( ztmelts-rt0 )         
436
437               zEw                = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
438
439               zdE                = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
440
441               dh_i_bott(ji)      = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
442
443               e_i_1d(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0)
444               
445            END DO
446            ! Contribution to Energy and Salt Fluxes                                   
447            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0)
448
449            ztmelts        = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K)
450           
451            zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
452           
453            zEi            = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)     
454               &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   &
455               &               + rcp  * ( ztmelts-rt0 )         
456           
457            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
458           
459            zdE            = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
460           
461            ! Contribution to heat flux to the ocean [W.m-2], >0 
462            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
463
464            ! Total heat flux used in this process [W.m-2], <0 
465            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
466           
467            ! Contribution to salt flux, <0
468            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice
469
470            ! Contribution to mass flux, <0
471            wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice
472
473            ! update heat content (J.m-2) and layer thickness
474            eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * e_i_1d(ji,nlay_i+1)
475            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji)
476
477         ENDIF
478
479      END DO
480
481      !----------------
482      ! 4.2 Basal melt
483      !----------------
484      zdeltah(1:nidx,:) = 0._wp ! important
485      DO jk = nlay_i, 1, -1
486         DO ji = 1, nidx
487            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting
488
489               ztmelts = - tmut * s_i_1d(ji,jk) + rt0  ! Melting point of layer jk (K)
490
491               IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
492
493                  zEi               = - e_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0)
494                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0)
495                                                                    ! set up at 0 since no energy is needed to melt water...(it is already melted)
496                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing     
497                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this
498
499                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk)
500
501                  zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0
502
503                  ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
504                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
505
506                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
507                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
508                                   
509                  ! Contribution to mass flux
510                  wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
511
512                  ! update heat content (J.m-2) and layer thickness
513                  eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
514                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
515
516               ELSE                               !!! Basal melting
517
518                  zEi             = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)
519                  zEw             = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0)
520                  zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0)
521
522                  zfmdt           = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0)
523
524                  zdeltah(ji,jk)  = - zfmdt * r1_rhoic         ! Gross thickness change
525
526                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
527                 
528                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors
529
530                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt
531
532                  zfmdt           = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0
533
534                  zQm             = zfmdt * zEw         ! Heat exchanged with ocean
535
536                  ! Contribution to heat flux to the ocean [W.m-2], <0 
537                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
538
539                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
540                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
541                 
542                  ! Total heat flux used in this process [W.m-2], >0 
543                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
544                 
545                  ! Contribution to mass flux
546                  wfx_bom_1d(ji)  =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
547
548                  ! update heat content (J.m-2) and layer thickness
549                  eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
550                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
551               ENDIF
552           
553            ENDIF
554         END DO
555      END DO
556
557      !-------------------------------------------
558      ! Update temperature, energy
559      !-------------------------------------------
560      DO ji = 1, nidx
561         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_bott(ji) )
562      END DO 
563
564      !-------------------------------------------
565      ! 5. What to do with remaining energy
566      !-------------------------------------------
567      ! If heat still available for melting and snow remains, then melt more snow
568      !-------------------------------------------
569      zdeltah(1:nidx,:) = 0._wp ! important
570      DO ji = 1, nidx
571         zq_rema(ji)     = zq_su(ji) + zq_bo(ji) 
572         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )   ! =1 if snow
573         rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) )
574         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 )
575         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting
576         dh_s_tot (ji)   = dh_s_tot(ji)  + zdeltah(ji,1)
577         ht_s_1d   (ji)  = ht_s_1d(ji)   + zdeltah(ji,1)
578       
579         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                ! update available heat (J.m-2)
580         ! heat used to melt snow
581         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0)
582         ! Contribution to mass flux
583         wfx_snw_sum_1d(ji)  =  wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
584         !   
585         ! Remaining heat flux (W.m-2) is sent to the ocean heat budget
586         hfx_out_1d(ji)  = hfx_out_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
587
588         IF( ln_limctl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
589      END DO
590
591      !
592      !------------------------------------------------------------------------------|
593      !  6) Snow-Ice formation                                                       |
594      !------------------------------------------------------------------------------|
595      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
596      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
597      DO ji = 1, nidx
598         !
599         dh_snowice(ji) = MAX(  0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic )  )
600
601         ht_i_1d(ji)    = ht_i_1d(ji) + dh_snowice(ji)
602         ht_s_1d(ji)    = ht_s_1d(ji) - dh_snowice(ji)
603
604         ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
605         zfmdt          = ( rhosn - rhoic ) * dh_snowice(ji)    ! <0
606         zEw            = rcp * sst_1d(ji)
607         zQm            = zfmdt * zEw 
608         
609         ! Contribution to heat flux
610         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
611
612         ! Contribution to salt flux
613         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice 
614
615         ! virtual salt flux to keep salinity constant
616         IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
617            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
618               &                            - sm_i_1d(ji) * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice    ! and get  rn_icesal from the ocean
619         ENDIF
620
621         ! Contribution to mass flux
622         ! All snow is thrown in the ocean, and seawater is taken to replace the volume
623         wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice
624         wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice
625
626         ! update heat content (J.m-2) and layer thickness
627         eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw
628         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
629         
630      END DO
631
632      !
633      !-------------------------------------------
634      ! Update temperature, energy
635      !-------------------------------------------
636      DO ji = 1, nidx
637         rswitch     =  1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 
638         t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0
639      END DO
640
641      DO jk = 1, nlay_s
642         DO ji = 1,nidx
643            ! mask enthalpy
644            rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) )  )
645            e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk)
646            ! recalculate t_s_1d from e_s_1d
647            t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic )
648         END DO
649      END DO
650
651      ! --- ensure that a_i = 0 where ht_i = 0 ---
652      DO ji = 1, nidx
653         IF( ht_i_1d(ji) == 0._wp )   a_i_1d(ji) = 0._wp
654      END DO         
655      !
656   END SUBROUTINE ice_thd_dh
657
658
659   !!--------------------------------------------------------------------------
660   !! INTERFACE ice_thd_snwblow
661   !! ** Purpose :   Compute distribution of precip over the ice
662   !!--------------------------------------------------------------------------
663!!gm  I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE....
664   SUBROUTINE ice_thd_snwblow_2d( pin, pout )
665      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( 1. - a_i_b )
666      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
667      pout = ( 1._wp - ( pin )**rn_betas )
668   END SUBROUTINE ice_thd_snwblow_2d
669
670   SUBROUTINE ice_thd_snwblow_1d( pin, pout )
671      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin
672      REAL(wp), DIMENSION(:), INTENT(inout) :: pout
673      pout = ( 1._wp - ( pin )**rn_betas )
674   END SUBROUTINE ice_thd_snwblow_1d
675
676#else
677   !!----------------------------------------------------------------------
678   !!   Default option                               NO  LIM3 sea-ice model
679   !!----------------------------------------------------------------------
680#endif
681
682   !!======================================================================
683END MODULE icethd_dh
Note: See TracBrowser for help on using the repository browser.