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.
limthd_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/limthd_dh.F90 @ 8341

Last change on this file since 8341 was 8341, checked in by clem, 8 years ago

correct the heat conservation (all fine except limthd_da for a reason I do not understand)

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