source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 6486

Last change on this file since 6486 was 6486, checked in by davestorkey, 5 years ago

Remove SVN keywords from UKMO/dev_r5518_GO6_package branch.

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