New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icethd_dh.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

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

Last change on this file since 8565 was 8565, checked in by clem, 6 years ago

trying to respect naming convention

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