source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 8414

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

changing names

  • Property svn:keywords set to Id
File size: 35.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   !
23   USE in_out_manager ! I/O manager
24   USE lib_mpp        ! MPP library
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
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  ::   ji , jk        ! dummy loop indices
69      INTEGER  ::   iter
70
71      REAL(wp) ::   ztmelts             ! local scalar
72      REAL(wp) ::   zdum       
73      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
74      REAL(wp) ::   zswi1        ! switch for computation of bottom salinity
75      REAL(wp) ::   zswi12       ! switch for computation of bottom salinity
76      REAL(wp) ::   zswi2        ! switch for computation of bottom salinity
77      REAL(wp) ::   zgrr         ! bottom growth rate
78      REAL(wp) ::   zt_i_new     ! bottom formation temperature
79
80      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
81      REAL(wp) ::   zEi          ! specific enthalpy of sea ice (J/kg)
82      REAL(wp) ::   zEw          ! specific enthalpy of exchanged water (J/kg)
83      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg)
84      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean
85
86      REAL(wp), DIMENSION(jpij) ::   zqprec      ! energy of fallen snow                       (J.m-3)
87      REAL(wp), DIMENSION(jpij) ::   zq_su       ! heat for surface ablation                   (J.m-2)
88      REAL(wp), DIMENSION(jpij) ::   zq_bo       ! heat for bottom ablation                    (J.m-2)
89      REAL(wp), DIMENSION(jpij) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2)
90      REAL(wp), DIMENSION(jpij) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2)
91      REAL(wp), DIMENSION(jpij) ::   zevap_rema  ! remaining mass flux from sublimation        (kg.m-2)
92
93      REAL(wp), DIMENSION(jpij) ::   zdh_s_mel   ! snow melt
94      REAL(wp), DIMENSION(jpij) ::   zdh_s_pre   ! snow precipitation
95      REAL(wp), DIMENSION(jpij) ::   zdh_s_sub   ! snow sublimation
96
97      REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah
98      REAL(wp), DIMENSION(jpij,nlay_i) ::   zh_i      ! ice layer thickness
99      INTEGER , DIMENSION(jpij,nlay_i) ::   icount    ! number of layers vanished by melting
100
101      REAL(wp), DIMENSION(jpij) ::   zeh_i       ! total ice heat content  (J.m-2)
102      REAL(wp), DIMENSION(jpij) ::   zsnw        ! distribution of snow after wind blowing
103
104      REAL(wp) :: zswitch_sal
105
106      ! Heat conservation
107      INTEGER  ::   num_iter_max
108      !!------------------------------------------------------------------
109
110      ! Discriminate between varying salinity (nn_icesal=2) and prescribed cases (other values)
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         ENDDO
130      ENDDO
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      !------------------------------------------------------------------------------!
146      ! If snow temperature is above freezing point, then snow melts
147      ! (should not happen but sometimes it does)
148      !------------------------------------------------------------------------------!
149      DO ji = 1, nidx
150         IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting
151            ! Contribution to heat flux to the ocean [W.m-2], < 0 
152            hfx_res_1d(ji) = hfx_res_1d(ji) + e_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
153            ! Contribution to mass flux
154            wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
155            ! updates
156            ht_s_1d(ji)   = 0._wp
157            e_s_1d (ji,1) = 0._wp
158            t_s_1d (ji,1) = rt0
159         END IF
160      END DO
161
162      !------------------------------------------------------------!
163      !  2) Computing layer thicknesses and enthalpies.            !
164      !------------------------------------------------------------!
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 lim_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
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_limctl .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      DO ji = 1, nidx
601         !
602         dh_snowice(ji) = MAX(  0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic )  )
603
604         ht_i_1d(ji)    = ht_i_1d(ji) + dh_snowice(ji)
605         ht_s_1d(ji)    = ht_s_1d(ji) - dh_snowice(ji)
606
607         ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
608         zfmdt          = ( rhosn - rhoic ) * dh_snowice(ji)    ! <0
609         zEw            = rcp * sst_1d(ji)
610         zQm            = zfmdt * zEw 
611         
612         ! Contribution to heat flux
613         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
614
615         ! Contribution to salt flux
616         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice 
617
618         ! virtual salt flux to keep salinity constant
619         IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
620            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
621               &                            - sm_i_1d(ji) * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice    ! and get  rn_icesal from the ocean
622         ENDIF
623
624         ! Contribution to mass flux
625         ! All snow is thrown in the ocean, and seawater is taken to replace the volume
626         wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice
627         wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice
628
629         ! update heat content (J.m-2) and layer thickness
630         eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw
631         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
632         
633      END DO
634
635      !
636      !-------------------------------------------
637      ! Update temperature, energy
638      !-------------------------------------------
639      DO ji = 1, nidx
640         rswitch     =  1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 
641         t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0
642      END DO
643
644      DO jk = 1, nlay_s
645         DO ji = 1,nidx
646            ! mask enthalpy
647            rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) )  )
648            e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk)
649            ! recalculate t_s_1d from e_s_1d
650            t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic )
651         END DO
652      END DO
653
654      ! --- ensure that a_i = 0 where ht_i = 0 ---
655      DO ji = 1, nidx
656         IF( ht_i_1d(ji) == 0._wp )   a_i_1d(ji) = 0._wp
657      END DO         
658      !
659   END SUBROUTINE lim_thd_dh
660
661
662   !!--------------------------------------------------------------------------
663   !! INTERFACE lim_thd_snwblow
664   !! ** Purpose :   Compute distribution of precip over the ice
665   !!--------------------------------------------------------------------------
666   SUBROUTINE lim_thd_snwblow_2d( pin, pout )
667      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( 1. - a_i_b )
668      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
669      pout = ( 1._wp - ( pin )**rn_betas )
670   END SUBROUTINE lim_thd_snwblow_2d
671
672   SUBROUTINE lim_thd_snwblow_1d( pin, pout )
673      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin
674      REAL(wp), DIMENSION(:), INTENT(inout) :: pout
675      pout = ( 1._wp - ( pin )**rn_betas )
676   END SUBROUTINE lim_thd_snwblow_1d
677
678   
679#else
680   !!----------------------------------------------------------------------
681   !!   Default option                               NO  LIM3 sea-ice model
682   !!----------------------------------------------------------------------
683CONTAINS
684   SUBROUTINE lim_thd_dh          ! Empty routine
685   END SUBROUTINE lim_thd_dh
686#endif
687
688   !!======================================================================
689END MODULE limthd_dh
Note: See TracBrowser for help on using the repository browser.