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

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

changes in style - part6 - pure cosmetics

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