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

Last change on this file since 8514 was 8514, checked in by clem, 7 years ago

changes in style - part5 - almost done

File size: 35.4 KB
Line 
1MODULE icethd_dh
2   !!======================================================================
3   !!                       ***  MODULE icethd_dh ***
4   !!  LIM-3 :   thermodynamic growth and decay of the ice
5   !!======================================================================
6   !! History :  LIM  ! 2003-05 (M. Vancoppenolle) Original code in 1D
7   !!                 ! 2005-06 (M. Vancoppenolle) 3D version
8   !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in wfx_snw & wfx_ice
9   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
10   !!            3.5  ! 2012-10 (G. Madec & co) salt flux + bug fixes
11   !!----------------------------------------------------------------------
12#if defined key_lim3
13   !!----------------------------------------------------------------------
14   !!   'key_lim3'                                       LIM3 sea-ice model
15   !!----------------------------------------------------------------------
16   !!   ice_thd_dh    : vertical accr./abl. and lateral ablation of sea ice
17   !!----------------------------------------------------------------------
18   USE par_oce        ! ocean parameters
19   USE phycst         ! physical constants (OCE directory)
20   USE ice            ! LIM variables
21   USE ice1D          ! LIM thermodynamics
22   USE icethd_sal     ! sea-ice: salinity profiles
23   !
24   USE in_out_manager ! I/O manager
25   USE lib_mpp        ! MPP library
26   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
27   
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   ice_thd_dh      ! called by ice_thd
32   PUBLIC   ice_thd_snwblow ! called in sbcblk/sbcclio/sbccpl and here
33
34   INTERFACE ice_thd_snwblow
35      MODULE PROCEDURE ice_thd_snwblow_1d, ice_thd_snwblow_2d
36   END INTERFACE
37
38   !!----------------------------------------------------------------------
39   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
40   !! $Id: icethd_dh.F90 8420 2017-08-08 12:18:46Z clem $
41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE ice_thd_dh
46      !!------------------------------------------------------------------
47      !!                ***  ROUTINE ice_thd_dh  ***
48      !!
49      !! ** Purpose :   determines variations of ice and snow thicknesses.
50      !!
51      !! ** Method  :   Ice/Snow surface melting arises from imbalance in surface fluxes
52      !!              Bottom accretion/ablation arises from flux budget
53      !!              Snow thickness can increase by precipitation and decrease by sublimation
54      !!              If snow load excesses Archmiede limit, snow-ice is formed by
55      !!              the flooding of sea-water in the snow
56      !!
57      !!                 1) Compute available flux of heat for surface ablation
58      !!                 2) Compute snow and sea ice enthalpies
59      !!                 3) Surface ablation and sublimation
60      !!                 4) Bottom accretion/ablation
61      !!                 5) Case of Total ablation
62      !!                 6) Snow ice formation
63      !!
64      !! References : Bitz and Lipscomb, 1999, J. Geophys. Res.
65      !!              Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
66      !!              Vancoppenolle, Fichefet and Bitz, 2005, Geophys. Res. Let.
67      !!              Vancoppenolle et al.,2009, Ocean Modelling
68      !!------------------------------------------------------------------
69      INTEGER  ::   ji, jk       ! dummy loop indices
70      INTEGER  ::   iter         ! local integer
71
72      REAL(wp) ::   ztmelts      ! local scalar
73      REAL(wp) ::   zdum       
74      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
75      REAL(wp) ::   zswi1        ! switch for computation of bottom salinity
76      REAL(wp) ::   zswi12       ! switch for computation of bottom salinity
77      REAL(wp) ::   zswi2        ! switch for computation of bottom salinity
78      REAL(wp) ::   zgrr         ! bottom growth rate
79      REAL(wp) ::   zt_i_new     ! bottom formation temperature
80
81      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
82      REAL(wp) ::   zEi          ! specific enthalpy of sea ice (J/kg)
83      REAL(wp) ::   zEw          ! specific enthalpy of exchanged water (J/kg)
84      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg)
85      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean
86
87      REAL(wp), DIMENSION(jpij) ::   zqprec      ! energy of fallen snow                       (J.m-3)
88      REAL(wp), DIMENSION(jpij) ::   zq_su       ! heat for surface ablation                   (J.m-2)
89      REAL(wp), DIMENSION(jpij) ::   zq_bo       ! heat for bottom ablation                    (J.m-2)
90      REAL(wp), DIMENSION(jpij) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2)
91      REAL(wp), DIMENSION(jpij) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2)
92      REAL(wp), DIMENSION(jpij) ::   zevap_rema  ! remaining mass flux from sublimation        (kg.m-2)
93
94      REAL(wp), DIMENSION(jpij) ::   zdh_s_mel   ! snow melt
95      REAL(wp), DIMENSION(jpij) ::   zdh_s_pre   ! snow precipitation
96      REAL(wp), DIMENSION(jpij) ::   zdh_s_sub   ! snow sublimation
97
98      REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah
99      REAL(wp), DIMENSION(jpij,nlay_i) ::   zh_i      ! ice layer thickness
100      INTEGER , DIMENSION(jpij,nlay_i) ::   icount    ! number of layers vanished by melting
101
102      REAL(wp), DIMENSION(jpij) ::   zeh_i       ! total ice heat content  (J.m-2)
103      REAL(wp), DIMENSION(jpij) ::   zsnw        ! distribution of snow after wind blowing
104
105      REAL(wp) ::   zswitch_sal
106
107      INTEGER  ::   num_iter_max      ! Heat conservation
108      !!------------------------------------------------------------------
109
110      ! Discriminate between time varying salinity and constant
111      SELECT CASE( nn_icesal )                  ! varying salinity or not
112         CASE( 1 , 3 )   ;   zswitch_sal = 0._wp   ! prescribed salinity profile
113         CASE( 2 )       ;   zswitch_sal = 1._wp   ! varying salinity profile
114      END SELECT
115
116      DO ji = 1, nidx
117         icount (ji,:) = 0
118         zdh_s_mel(ji) = 0._wp
119         e_i_1d(ji,nlay_i+1) = 0._wp ! Initialize enthalpy at nlay_i+1
120      END DO
121
122      ! initialize layer thicknesses and enthalpies
123      h_i_old (1:nidx,0:nlay_i+1) = 0._wp
124      eh_i_old(1:nidx,0:nlay_i+1) = 0._wp
125      DO jk = 1, nlay_i
126         DO ji = 1, nidx
127            h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i
128            eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk)
129         END DO
130      END DO
131      !
132      !------------------------------------------------------------------------------!
133      !  1) Calculate available heat for surface and bottom ablation                 !
134      !------------------------------------------------------------------------------!
135      !
136      DO ji = 1, nidx
137         zdum       = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
138         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 
139
140         zq_su (ji) = MAX( 0._wp, zdum      * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
141         zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
142      END DO
143
144      !------------------------------------------------------------------------------!
145      ! If snow temperature is above freezing point, then snow melts
146      ! (should not happen but sometimes it does)
147      !------------------------------------------------------------------------------!
148      DO ji = 1, nidx
149         IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting
150            ! Contribution to heat flux to the ocean [W.m-2], < 0 
151            hfx_res_1d(ji) = hfx_res_1d(ji) + e_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
152            ! Contribution to mass flux
153            wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
154            ! updates
155            ht_s_1d(ji)   = 0._wp
156            e_s_1d (ji,1) = 0._wp
157            t_s_1d (ji,1) = rt0
158         END IF
159      END DO
160
161      !------------------------------------------------------------!
162      !  2) Computing layer thicknesses and enthalpies.            !
163      !------------------------------------------------------------!
164      DO jk = 1, nlay_i
165         DO ji = 1, nidx
166            zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i
167            zeh_i(ji)   = zeh_i(ji) + e_i_1d(ji,jk) * zh_i(ji,jk)
168         END DO
169      END DO
170
171      !------------------------------------------------------------------------------|
172      !  3) Surface ablation and sublimation                                         |
173      !------------------------------------------------------------------------------|
174      !
175      !-------------------------
176      ! 3.1 Snow precips / melt
177      !-------------------------
178      ! Snow accumulation in one thermodynamic time step
179      ! snowfall is partitionned between leads and ice
180      ! if snow fall was uniform, a fraction (1-at_i) would fall into leads
181      ! but because of the winds, more snow falls on leads than on sea ice
182      ! and a greater fraction (1-at_i)^beta of the total mass of snow
183      ! (beta < 1) falls in leads.
184      ! In reality, beta depends on wind speed,
185      ! and should decrease with increasing wind speed but here, it is
186      ! considered as a constant. an average value is 0.66
187      ! Martin Vancoppenolle, December 2006
188
189      CALL ice_thd_snwblow( 1. - at_i_1d(1:nidx), zsnw(1:nidx) ) ! snow distribution over ice after wind blowing
190
191      zdeltah(1:nidx,:) = 0._wp
192      DO ji = 1, nidx
193         !-----------
194         ! Snow fall
195         !-----------
196         ! thickness change
197         zdh_s_pre(ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji)
198         ! enthalpy of the precip (>0, J.m-3)
199         zqprec   (ji) = - qprec_ice_1d(ji)   
200         IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp
201         ! heat flux from snow precip (>0, W.m-2)
202         hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
203         ! mass flux, <0
204         wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice
205
206         !---------------------
207         ! Melt of falling snow
208         !---------------------
209         ! thickness change
210         rswitch        = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )
211         zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 )
212         zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting
213         ! heat used to melt snow (W.m-2, >0)
214         hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
215         ! snow melting only = water into the ocean (then without snow precip), >0
216         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
217         ! updates available heat + precipitations after melting
218         zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )     
219         zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1)
220
221         ! update thickness
222         ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) )
223      END DO
224
225      ! If heat still available (zq_su > 0), then melt more snow
226      zdeltah(1:nidx,:) = 0._wp
227      DO jk = 1, nlay_s
228         DO ji = 1, nidx
229            ! thickness change
230            rswitch          = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 
231            rswitch          = rswitch * ( MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) ) ) 
232            zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( e_s_1d(ji,jk), epsi20 )
233            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - ht_s_1d(ji) ) ! bound melting
234            zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)   
235            ! heat used to melt snow(W.m-2, >0)
236            hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d(ji,jk) * r1_rdtice 
237            ! snow melting only = water into the ocean (then without snow precip)
238            wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
239            ! updates available heat + thickness
240            zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
241            ht_s_1d(ji) = MAX( 0._wp , ht_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:nidx,:) = 0._wp
251      DO ji = 1, nidx
252         zdh_s_sub(ji)  = MAX( - ht_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         ht_s_1d(ji)    =  MAX( 0._wp , ht_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, nidx
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,nidx
280            rswitch       = MAX( 0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 ) )
281            e_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *           &
282              &            ( ( zdh_s_pre(ji)               ) * zqprec(ji) +  &
283              &              ( ht_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:nidx,:) = 0._wp ! important
291      DO jk = 1, nlay_i
292         DO ji = 1, nidx
293            ztmelts           = - tmut * s_i_1d(ji,jk) + rt0          ! Melting point of layer k [K]
294           
295            IF( t_i_1d(ji,jk) >= ztmelts ) 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 sm_i_1d and not s_i_1d(jk) is ok)
312               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_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 - rt0 )            ! 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 sm_i_1d and not s_i_1d(jk) is ok)
338               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_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 * sm_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, nidx
381         ht_i_1d(ji) =  MAX( 0._wp , ht_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, nidx
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, nidx
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 ) * sm_i_1d(ji) 
429               ! New ice growth
430               ztmelts            = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K)
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 ) &     ! Specific enthalpy of forming ice (J/kg, <0)     
435                  &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   &
436                  &               + rcp  * ( ztmelts-rt0 )         
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               e_i_1d(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0)
445               
446            END DO
447            ! Contribution to Energy and Salt Fluxes                                   
448            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0)
449
450            ztmelts        = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K)
451           
452            zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
453           
454            zEi            = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)     
455               &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   &
456               &               + rcp  * ( ztmelts-rt0 )         
457           
458            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
459           
460            zdE            = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
461           
462            ! Contribution to heat flux to the ocean [W.m-2], >0 
463            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
464
465            ! Total heat flux used in this process [W.m-2], <0 
466            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
467           
468            ! Contribution to salt flux, <0
469            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice
470
471            ! Contribution to mass flux, <0
472            wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice
473
474            ! update heat content (J.m-2) and layer thickness
475            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)
476            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji)
477
478         ENDIF
479
480      END DO
481
482      !----------------
483      ! 4.2 Basal melt
484      !----------------
485      zdeltah(1:nidx,:) = 0._wp ! important
486      DO jk = nlay_i, 1, -1
487         DO ji = 1, nidx
488            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting
489
490               ztmelts = - tmut * s_i_1d(ji,jk) + rt0  ! Melting point of layer jk (K)
491
492               IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
493
494                  zEi               = - e_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0)
495                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0)
496                                                                    ! set up at 0 since no energy is needed to melt water...(it is already melted)
497                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing     
498                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this
499
500                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk)
501
502                  zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0
503
504                  ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
505                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
506
507                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
508                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
509                                   
510                  ! Contribution to mass flux
511                  wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
512
513                  ! update heat content (J.m-2) and layer thickness
514                  eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
515                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
516
517               ELSE                               !!! Basal melting
518
519                  zEi             = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)
520                  zEw             = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0)
521                  zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0)
522
523                  zfmdt           = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0)
524
525                  zdeltah(ji,jk)  = - zfmdt * r1_rhoic         ! Gross thickness change
526
527                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
528                 
529                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors
530
531                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt
532
533                  zfmdt           = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0
534
535                  zQm             = zfmdt * zEw         ! Heat exchanged with ocean
536
537                  ! Contribution to heat flux to the ocean [W.m-2], <0 
538                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
539
540                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
541                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
542                 
543                  ! Total heat flux used in this process [W.m-2], >0 
544                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
545                 
546                  ! Contribution to mass flux
547                  wfx_bom_1d(ji)  =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
548
549                  ! update heat content (J.m-2) and layer thickness
550                  eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
551                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
552               ENDIF
553           
554            ENDIF
555         END DO
556      END DO
557
558      !-------------------------------------------
559      ! Update temperature, energy
560      !-------------------------------------------
561      DO ji = 1, nidx
562         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_bott(ji) )
563      END DO 
564
565      !-------------------------------------------
566      ! 5. What to do with remaining energy
567      !-------------------------------------------
568      ! If heat still available for melting and snow remains, then melt more snow
569      !-------------------------------------------
570      zdeltah(1:nidx,:) = 0._wp ! important
571      DO ji = 1, nidx
572         zq_rema(ji)     = zq_su(ji) + zq_bo(ji) 
573         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )   ! =1 if snow
574         rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) )
575         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 )
576         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting
577         dh_s_tot (ji)   = dh_s_tot(ji)  + zdeltah(ji,1)
578         ht_s_1d   (ji)  = ht_s_1d(ji)   + zdeltah(ji,1)
579       
580         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                ! update available heat (J.m-2)
581         ! heat used to melt snow
582         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)
583         ! Contribution to mass flux
584         wfx_snw_sum_1d(ji)  =  wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
585         !   
586         ! Remaining heat flux (W.m-2) is sent to the ocean heat budget
587         hfx_out_1d(ji)  = hfx_out_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
588
589         IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
590      END DO
591
592      !
593      !------------------------------------------------------------------------------|
594      !  6) Snow-Ice formation                                                       |
595      !------------------------------------------------------------------------------|
596      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
597      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
598      DO ji = 1, nidx
599         !
600         dh_snowice(ji) = MAX(  0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic )  )
601
602         ht_i_1d(ji)    = ht_i_1d(ji) + dh_snowice(ji)
603         ht_s_1d(ji)    = ht_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               &                            - sm_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, nidx
638         rswitch     =  1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_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,nidx
644            ! mask enthalpy
645            rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp, - ht_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) / ( rhosn * cpic ) + lfus / cpic )
649         END DO
650      END DO
651
652      ! --- ensure that a_i = 0 where ht_i = 0 ---
653      DO ji = 1, nidx
654         IF( ht_i_1d(ji) == 0._wp )   a_i_1d(ji) = 0._wp
655      END DO         
656      !
657   END SUBROUTINE ice_thd_dh
658
659
660   !!--------------------------------------------------------------------------
661   !! INTERFACE ice_thd_snwblow
662   !! ** Purpose :   Compute distribution of precip over the ice
663   !!--------------------------------------------------------------------------
664!!gm  I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE....
665   SUBROUTINE ice_thd_snwblow_2d( pin, pout )
666      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( 1. - a_i_b )
667      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
668      pout = ( 1._wp - ( pin )**rn_blow_s )
669   END SUBROUTINE ice_thd_snwblow_2d
670
671   SUBROUTINE ice_thd_snwblow_1d( pin, pout )
672      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin
673      REAL(wp), DIMENSION(:), INTENT(inout) :: pout
674      pout = ( 1._wp - ( pin )**rn_blow_s )
675   END SUBROUTINE ice_thd_snwblow_1d
676
677#else
678   !!----------------------------------------------------------------------
679   !!   Default option                               NO  LIM3 sea-ice model
680   !!----------------------------------------------------------------------
681#endif
682
683   !!======================================================================
684END MODULE icethd_dh
Note: See TracBrowser for help on using the repository browser.