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 @ 8623

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

rewrite melt ponds

File size: 35.4 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            dh_s_mlt(ji)       = dh_s_mlt(ji) - h_s_1d(ji)
151            ! updates
152            h_s_1d(ji)    = 0._wp
153            e_s_1d (ji,1) = 0._wp
154            t_s_1d (ji,1) = rt0
155         END IF
156      END DO
157
158      !------------------------------------------------------------!
159      !  2) Computing layer thicknesses and enthalpies.            !
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            zeh_i(ji)   = zeh_i(ji) + e_i_1d(ji,jk) * zh_i(ji,jk)
165         END DO
166      END DO
167
168      !------------------------------------------------------------------------------|
169      !  3) Surface ablation and sublimation                                         |
170      !------------------------------------------------------------------------------|
171      !
172      !-------------------------
173      ! 3.1 Snow precips / melt
174      !-------------------------
175      ! Snow accumulation in one thermodynamic time step
176      ! snowfall is partitionned between leads and ice
177      ! if snow fall was uniform, a fraction (1-at_i) would fall into leads
178      ! but because of the winds, more snow falls on leads than on sea ice
179      ! and a greater fraction (1-at_i)^beta of the total mass of snow
180      ! (beta < 1) falls in leads.
181      ! In reality, beta depends on wind speed,
182      ! and should decrease with increasing wind speed but here, it is
183      ! considered as a constant. an average value is 0.66
184      ! Martin Vancoppenolle, December 2006
185
186      CALL ice_thd_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing
187
188      zdeltah(1:npti,:) = 0._wp
189      DO ji = 1, npti
190         !-----------
191         ! Snow fall
192         !-----------
193         ! thickness change
194         zdh_s_pre(ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji)
195         ! enthalpy of the precip (>0, J.m-3)
196         zqprec   (ji) = - qprec_ice_1d(ji)   
197         IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp
198         ! heat flux from snow precip (>0, W.m-2)
199         hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
200         ! mass flux, <0
201         wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice
202
203         !---------------------
204         ! Melt of falling snow
205         !---------------------
206         ! thickness change
207         rswitch        = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )
208         zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 )
209         zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting
210         ! heat used to melt snow (W.m-2, >0)
211         hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
212         ! snow melting only = water into the ocean (then without snow precip), >0
213         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
214         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1)
215         ! updates available heat + precipitations after melting
216         zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )     
217         zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1)
218
219         ! update thickness
220         h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) )
221      END DO
222
223      ! If heat still available (zq_su > 0), then melt more snow
224      zdeltah(1:npti,:) = 0._wp
225      zdh_s_mel(1:npti) = 0._wp
226      DO jk = 1, nlay_s
227         DO ji = 1, npti
228            ! thickness change
229            rswitch          = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_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) , - h_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            dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,jk)
239            ! updates available heat + thickness
240            zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
241            h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) )
242         END DO
243      END DO
244
245      !------------------------------
246      ! 3.2 Sublimation (part1: snow)
247      !------------------------------
248      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
249      ! clem comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean)
250      zdeltah(1:npti,:) = 0._wp
251      DO ji = 1, npti
252         zdh_s_sub(ji)  = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice )
253         ! remaining evap in kg.m-2 (used for ice melting later on)
254         zevap_rema(ji)  = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhosn
255         ! Heat flux by sublimation [W.m-2], < 0 (sublimate first snow that had fallen, then pre-existing snow)
256         zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
257         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)  &
258            &                              ) * a_i_1d(ji) * r1_rdtice
259         ! Mass flux by sublimation
260         wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice
261
262         ! new snow thickness
263         h_s_1d(ji)    =  MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) )
264         ! update precipitations after sublimation and correct sublimation
265         zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
266         zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1)
267      END DO
268     
269      ! --- Update snow diags --- !
270      DO ji = 1, npti
271         dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
272      END DO
273
274      !-------------------------------------------
275      ! 3.3 Update temperature, energy
276      !-------------------------------------------
277      ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
278      DO jk = 1, nlay_s
279         DO ji = 1,npti
280            rswitch       = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) )
281            e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) *           &
282              &            ( ( zdh_s_pre(ji)               ) * zqprec(ji) +  &
283              &              ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
284         END DO
285      END DO
286
287      !--------------------------
288      ! 3.4 Surface ice ablation
289      !--------------------------
290      zdeltah(1:npti,:) = 0._wp ! important
291      DO jk = 1, nlay_i
292         DO ji = 1, npti
293            ztmelts           = - tmut * sz_i_1d(ji,jk)          ! Melting point of layer k [C]
294           
295            IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !!! Internal melting
296
297               zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]       
298               zdE            = 0._wp                                 ! Specific enthalpy difference   (J/kg, <0)
299                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted)
300               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing     
301                                                                      ! this should normally not happen, but sometimes, heat diffusion leads to this
302               zfmdt          = - zdeltah(ji,jk) * rhoic              ! Mass flux x time step > 0
303                         
304               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt
305               
306               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0]
307
308               ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
309               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
310               
311               ! Contribution to salt flux (clem: using s_i_1d and not sz_i_1d(jk) is ok)
312               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice
313               
314               ! Contribution to mass flux
315               wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
316
317            ELSE                                !!! Surface melting
318               
319               zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]
320               zEw            =    rcp * ztmelts                      ! Specific enthalpy of resulting meltwater [J/kg, <0]
321               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0
322               
323               zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0]
324               
325               zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0]
326               
327               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]
328               
329               zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat
330               
331               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt
332               
333               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0]
334               
335               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0]
336               
337               ! Contribution to salt flux >0 (clem: using s_i_1d and not sz_i_1d(jk) is ok)
338               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice
339               
340               ! Contribution to heat flux [W.m-2], < 0
341               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
342               
343               ! Total heat flux used in this process [W.m-2], > 0 
344               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
345               
346               ! Contribution to mass flux
347               wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
348               
349            END IF
350            ! ----------------------
351            ! Sublimation part2: ice
352            ! ----------------------
353            zdum      = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoic )
354            zdeltah(ji,jk) = zdeltah(ji,jk) + zdum
355            dh_i_sub(ji)  = dh_i_sub(ji) + zdum
356            ! 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.
357            !                          It must be corrected at some point)
358            sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_rdtice
359            ! Heat flux [W.m-2], < 0
360            hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice
361            ! Mass flux > 0
362            wfx_ice_sub_1d(ji) =  wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice
363
364            ! update remaining mass flux
365            zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoic
366           
367            ! record which layers have disappeared (for bottom melting)
368            !    => icount=0 : no layer has vanished
369            !    => icount=5 : 5 layers have vanished
370            rswitch       = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) 
371            icount(ji,jk) = NINT( rswitch )
372            zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
373                       
374            ! update heat content (J.m-2) and layer thickness
375            eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
376            h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
377         END DO
378      END DO
379      ! update ice thickness
380      DO ji = 1, npti
381         h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_surf(ji) + dh_i_sub(ji) )
382      END DO
383
384      ! remaining "potential" evap is sent to ocean
385      DO ji = 1, npti
386         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)
387      END DO
388
389
390      !------------------------------------------------------------------------------!
391      ! 4) Basal growth / melt                                                       !
392      !------------------------------------------------------------------------------!
393      !
394      !------------------
395      ! 4.1 Basal growth
396      !------------------
397      ! Basal growth is driven by heat imbalance at the ice-ocean interface,
398      ! between the inner conductive flux  (fc_bo_i), from the open water heat flux
399      ! (fhld) and the turbulent ocean flux (fhtur).
400      ! fc_bo_i is positive downwards. fhtur and fhld are positive to the ice
401
402      ! If salinity varies in time, an iterative procedure is required, because
403      ! the involved quantities are inter-dependent.
404      ! Basal growth (dh_i_bott) depends upon new ice specific enthalpy (zEi),
405      ! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bott)
406      ! -> need for an iterative procedure, which converges quickly
407
408      num_iter_max = 1
409      IF( nn_icesal == 2 )   num_iter_max = 5  ! salinity varying in time
410     
411      ! Iterative procedure
412      DO ji = 1, npti
413         IF(  zf_tt(ji) < 0._wp  ) THEN
414            DO iter = 1, num_iter_max
415
416               ! New bottom ice salinity (Cox & Weeks, JGR88 )
417               !--- zswi1  if dh/dt < 2.0e-8
418               !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
419               !--- zswi2  if dh/dt > 3.6e-7
420               zgrr               = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) )
421               zswi2              = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
422               zswi12             = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
423               zswi1              = 1. - zswi2 * zswi12
424               zfracs             = MIN ( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
425                  &               + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 )
426
427               s_i_new(ji)        = zswitch_sal * zfracs * sss_1d(ji)  &  ! New ice salinity
428                                  + ( 1. - zswitch_sal ) * s_i_1d(ji) 
429               ! New ice growth
430               ztmelts            = - tmut * s_i_new(ji)                 ! New ice melting point (C)
431
432               zt_i_new           = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
433               
434               zEi                = cpic * ( zt_i_new - (ztmelts+rt0) ) &     ! Specific enthalpy of forming ice (J/kg, <0)     
435                  &               - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) )   &
436                  &               + rcp  * ztmelts         
437
438               zEw                = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
439
440               zdE                = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
441
442               dh_i_bott(ji)      = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
443               
444            END DO
445            ! Contribution to Energy and Salt Fluxes                                   
446            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0)
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 ) )   &
454               &               + rcp  * ztmelts         
455           
456            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
457           
458            zdE            = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
459           
460            ! Contribution to heat flux to the ocean [W.m-2], >0 
461            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
462
463            ! Total heat flux used in this process [W.m-2], <0 
464            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
465           
466            ! Contribution to salt flux, <0
467            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice
468
469            ! Contribution to mass flux, <0
470            wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice
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_bott(ji) * (-zEi * rhoic)
474            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji)
475
476         ENDIF
477
478      END DO
479
480      !----------------
481      ! 4.2 Basal melt
482      !----------------
483      zdeltah(1:npti,:) = 0._wp ! important
484      DO jk = nlay_i, 1, -1
485         DO ji = 1, npti
486            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting
487
488               ztmelts = - tmut * sz_i_1d(ji,jk)  ! Melting point of layer jk (C)
489
490               IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !!! Internal melting
491
492                  zEi               = - e_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0)
493                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0)
494                                                                    ! set up at 0 since no energy is needed to melt water...(it is already melted)
495                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing     
496                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this
497
498                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk)
499
500                  zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0
501
502                  ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
503                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
504
505                  ! Contribution to salt flux (clem: using s_i_1d and not sz_i_1d(jk) is ok)
506                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice
507                                   
508                  ! Contribution to mass flux
509                  wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
510
511                  ! update heat content (J.m-2) and layer thickness
512                  eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
513                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
514
515               ELSE                               !!! Basal melting
516
517                  zEi             = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)
518                  zEw             = rcp * ztmelts              ! Specific enthalpy of meltwater (J/kg, <0)
519                  zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0)
520
521                  zfmdt           = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0)
522
523                  zdeltah(ji,jk)  = - zfmdt * r1_rhoic         ! Gross thickness change
524
525                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
526                 
527                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors
528
529                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt
530
531                  zfmdt           = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0
532
533                  zQm             = zfmdt * zEw         ! Heat exchanged with ocean
534
535                  ! Contribution to heat flux to the ocean [W.m-2], <0 
536                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
537
538                  ! Contribution to salt flux (clem: using s_i_1d and not sz_i_1d(jk) is ok)
539                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice
540                 
541                  ! Total heat flux used in this process [W.m-2], >0 
542                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
543                 
544                  ! Contribution to mass flux
545                  wfx_bom_1d(ji)  =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
546
547                  ! update heat content (J.m-2) and layer thickness
548                  eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
549                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
550               ENDIF
551           
552            ENDIF
553         END DO
554      END DO
555
556      !-------------------------------------------
557      ! Update temperature, energy
558      !-------------------------------------------
559      DO ji = 1, npti
560         h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_bott(ji) )
561      END DO 
562
563      !-------------------------------------------
564      ! 5. What to do with remaining energy
565      !-------------------------------------------
566      ! If heat still available for melting and snow remains, then melt more snow
567      !-------------------------------------------
568      zdeltah(1:npti,:) = 0._wp ! important
569      DO ji = 1, npti
570         zq_rema(ji)     = zq_su(ji) + zq_bo(ji) 
571         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) )   ! =1 if snow
572         rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) )
573         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 )
574         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting
575         dh_s_tot (ji)   = dh_s_tot(ji)  + zdeltah(ji,1)
576         h_s_1d   (ji)   = h_s_1d(ji)   + zdeltah(ji,1)
577       
578         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                ! update available heat (J.m-2)
579         ! heat used to melt snow
580         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)
581         ! Contribution to mass flux
582         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
583         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1)
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_icectl .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      z1_rho = 1._wp / ( rhosn+rau0-rhoic )
598      DO ji = 1, npti
599         !
600         dh_snowice(ji) = MAX(  0._wp , ( rhosn * h_s_1d(ji) + (rhoic-rau0) * h_i_1d(ji) ) * z1_rho )
601
602         h_i_1d(ji)    = h_i_1d(ji) + dh_snowice(ji)
603         h_s_1d(ji)    = h_s_1d(ji) - dh_snowice(ji)
604
605         ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
606         zfmdt          = ( rhosn - rhoic ) * dh_snowice(ji)    ! <0
607         zEw            = rcp * sst_1d(ji)
608         zQm            = zfmdt * zEw 
609         
610         ! Contribution to heat flux
611         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
612
613         ! Contribution to salt flux
614         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice 
615
616         ! Case constant salinity in time: virtual salt flux to keep salinity constant
617         IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
618            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
619               &                            - s_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice    ! and get  rn_icesal from the ocean
620         ENDIF
621
622         ! Contribution to mass flux
623         ! All snow is thrown in the ocean, and seawater is taken to replace the volume
624         wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice
625         wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice
626
627         ! update heat content (J.m-2) and layer thickness
628         eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw
629         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
630         
631      END DO
632
633      !
634      !-------------------------------------------
635      ! Update temperature, energy
636      !-------------------------------------------
637      DO ji = 1, npti
638         rswitch     =  1.0 - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 
639         t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0
640      END DO
641
642      DO jk = 1, nlay_s
643         DO ji = 1,npti
644            ! mask enthalpy
645            rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp, - h_s_1d(ji) )  )
646            e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk)
647            ! recalculate t_s_1d from e_s_1d
648            t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhosn * r1_cpic + lfus * r1_cpic )
649         END DO
650      END DO
651
652      ! --- ensure that a_i = 0 where h_i = 0 ---
653      WHERE( h_i_1d(1:npti) == 0._wp )   a_i_1d(1:npti) = 0._wp
654      !
655   END SUBROUTINE ice_thd_dh
656
657
658   !!--------------------------------------------------------------------------
659   !! INTERFACE ice_thd_snwblow
660   !! ** Purpose :   Compute distribution of precip over the ice
661   !!--------------------------------------------------------------------------
662!!gm  I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE....
663   SUBROUTINE ice_thd_snwblow_2d( pin, pout )
664      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( 1. - a_i_b )
665      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
666      pout = ( 1._wp - ( pin )**rn_blow_s )
667   END SUBROUTINE ice_thd_snwblow_2d
668
669   SUBROUTINE ice_thd_snwblow_1d( pin, pout )
670      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin
671      REAL(wp), DIMENSION(:), INTENT(inout) :: pout
672      pout = ( 1._wp - ( pin )**rn_blow_s )
673   END SUBROUTINE ice_thd_snwblow_1d
674
675#else
676   !!----------------------------------------------------------------------
677   !!   Default option                                NO ESIM sea-ice model
678   !!----------------------------------------------------------------------
679#endif
680
681   !!======================================================================
682END MODULE icethd_dh
Note: See TracBrowser for help on using the repository browser.