source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 5443

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

Update 2015/dev_r5021_UKMO1_CICE_coupling branch to revision 5442 of the trunk.

File size: 36.6 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/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 )
125      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s, zsnw )
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
134
135      zdh_s_pre(:) = 0._wp
136      zdh_s_mel(:) = 0._wp
137      zdh_s_sub(:) = 0._wp
138      zqh_s    (:) = 0._wp     
139      zqh_i    (:) = 0._wp   
140
141      zh_i      (:,:) = 0._wp       
142      zdeltah   (:,:) = 0._wp       
143      icount    (:,:) = 0
144
145      ! Initialize enthalpy at nlay_i+1
146      DO ji = kideb, kiut
147         q_i_1d(ji,nlay_i+1) = 0._wp
148      END DO
149
150      ! initialize layer thicknesses and enthalpies
151      h_i_old (:,0:nlay_i+1) = 0._wp
152      qh_i_old(:,0:nlay_i+1) = 0._wp
153      DO jk = 1, nlay_i
154         DO ji = kideb, kiut
155            h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i
156            qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk)
157         ENDDO
158      ENDDO
159      !
160      !------------------------------------------------------------------------------!
161      !  1) Calculate available heat for surface and bottom ablation                 !
162      !------------------------------------------------------------------------------!
163      !
164      DO ji = kideb, kiut
165         zfdum      = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
166         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 
167
168         zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
169         zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
170      END DO
171
172      !
173      !------------------------------------------------------------------------------!
174      ! If snow temperature is above freezing point, then snow melts
175      ! (should not happen but sometimes it does)
176      !------------------------------------------------------------------------------!
177      DO ji = kideb, kiut
178         IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting
179            ! Contribution to heat flux to the ocean [W.m-2], < 0 
180            hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
181            ! Contribution to mass flux
182            wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
183            ! updates
184            ht_s_1d(ji)   = 0._wp
185            q_s_1d (ji,1) = 0._wp
186            t_s_1d (ji,1) = rt0
187         END IF
188      END DO
189
190      !------------------------------------------------------------!
191      !  2) Computing layer thicknesses and enthalpies.            !
192      !------------------------------------------------------------!
193      !
194      DO jk = 1, nlay_s
195         DO ji = kideb, kiut
196            zqh_s(ji) =  zqh_s(ji) + q_s_1d(ji,jk) * ht_s_1d(ji) * r1_nlay_s
197         END DO
198      END DO
199      !
200      DO jk = 1, nlay_i
201         DO ji = kideb, kiut
202            zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i
203            zqh_i(ji)   = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk)
204         END DO
205      END DO
206      !
207      !------------------------------------------------------------------------------|
208      !  3) Surface ablation and sublimation                                         |
209      !------------------------------------------------------------------------------|
210      !
211      !-------------------------
212      ! 3.1 Snow precips / melt
213      !-------------------------
214      ! Snow accumulation in one thermodynamic time step
215      ! snowfall is partitionned between leads and ice
216      ! if snow fall was uniform, a fraction (1-at_i) would fall into leads
217      ! but because of the winds, more snow falls on leads than on sea ice
218      ! and a greater fraction (1-at_i)^beta of the total mass of snow
219      ! (beta < 1) falls in leads.
220      ! In reality, beta depends on wind speed,
221      ! and should decrease with increasing wind speed but here, it is
222      ! considered as a constant. an average value is 0.66
223      ! Martin Vancoppenolle, December 2006
224
225      zdeltah(:,:) = 0._wp
226      CALL lim_thd_snwblow( 1. - at_i_1d, zsnw ) ! snow distribution over ice after wind blowing
227      DO ji = kideb, kiut
228         !-----------
229         ! Snow fall
230         !-----------
231         ! thickness change
232         zdh_s_pre(ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji)
233         ! enthalpy of the precip (>0, J.m-3)
234         zqprec   (ji) = - qprec_ice_1d(ji)   
235         IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp
236         ! heat flux from snow precip (>0, W.m-2)
237         hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
238         ! mass flux, <0
239         wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice
240
241         !---------------------
242         ! Melt of falling snow
243         !---------------------
244         ! thickness change
245         rswitch        = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )
246         zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 )
247         zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting
248         ! heat used to melt snow (W.m-2, >0)
249         hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
250         ! snow melting only = water into the ocean (then without snow precip), >0
251         wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice   
252         ! updates available heat + precipitations after melting
253         zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )     
254         zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1)
255
256         ! update thickness
257         ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) )
258      END DO
259
260      ! If heat still available (zq_su > 0), then melt more snow
261      zdeltah(:,:) = 0._wp
262      DO jk = 1, nlay_s
263         DO ji = kideb, kiut
264            ! thickness change
265            rswitch          = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 
266            rswitch          = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) ) 
267            zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 )
268            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - ht_s_1d(ji) ) ! bound melting
269            zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)   
270            ! heat used to melt snow(W.m-2, >0)
271            hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice 
272            ! snow melting only = water into the ocean (then without snow precip)
273            wfx_snw_1d(ji)   = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
274            ! updates available heat + thickness
275            zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) )
276            ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) )
277         END DO
278      END DO
279
280      !----------------------
281      ! 3.2 Snow sublimation
282      !----------------------
283      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
284      ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean)
285      ! clem comment: ice should also sublimate
286      zdeltah(:,:) = 0._wp
287      ! coupled mode: sublimation is set to 0 (evap_ice = 0) until further notice
288      ! forced  mode: snow thickness change due to sublimation
289      DO ji = kideb, kiut
290         zdh_s_sub(ji)  =  MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice )
291         ! Heat flux by sublimation [W.m-2], < 0
292         !      sublimate first snow that had fallen, then pre-existing snow
293         zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
294         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)  &
295            &                              ) * a_i_1d(ji) * r1_rdtice
296         ! Mass flux by sublimation
297         wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice
298         ! new snow thickness
299         ht_s_1d(ji)    =  MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) )
300         ! update precipitations after sublimation and correct sublimation
301         zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
302         zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1)
303      END DO
304     
305      ! --- Update snow diags --- !
306      DO ji = kideb, kiut
307         dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
308      END DO
309
310      !-------------------------------------------
311      ! 3.3 Update temperature, energy
312      !-------------------------------------------
313      ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
314      zq_s(:) = 0._wp 
315      DO jk = 1, nlay_s
316         DO ji = kideb,kiut
317            rswitch       = MAX(  0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 )  )
318            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *                          &
319              &            ( (   zdh_s_pre(ji)             ) * zqprec(ji) +  &
320              &              (   ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
321            zq_s(ji)     =  zq_s(ji) + q_s_1d(ji,jk)
322         END DO
323      END DO
324
325      !--------------------------
326      ! 3.4 Surface ice ablation
327      !--------------------------
328      zdeltah(:,:) = 0._wp ! important
329      DO jk = 1, nlay_i
330         DO ji = kideb, kiut
331            ztmelts           = - tmut * s_i_1d(ji,jk) + rt0          ! Melting point of layer k [K]
332           
333            IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
334
335               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]       
336               zdE            = 0._wp                                 ! Specific enthalpy difference   (J/kg, <0)
337                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted)
338               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing     
339                                                                      ! this should normally not happen, but sometimes, heat diffusion leads to this
340               zfmdt          = - zdeltah(ji,jk) * rhoic              ! Mass flux x time step > 0
341                         
342               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt
343               
344               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0]
345
346               ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
347               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
348               
349               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
350               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
351               
352               ! Contribution to mass flux
353               wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
354
355            ELSE                                !!! Surface melting
356               
357               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]
358               zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0]
359               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0
360               
361               zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0]
362               
363               zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0]
364               
365               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]
366               
367               zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat
368               
369               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt
370               
371               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0]
372               
373               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0]
374               
375               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
376               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
377               
378               ! Contribution to heat flux [W.m-2], < 0
379               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
380               
381               ! Total heat flux used in this process [W.m-2], > 0 
382               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
383               
384               ! Contribution to mass flux
385               wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
386               
387            END IF
388            ! record which layers have disappeared (for bottom melting)
389            !    => icount=0 : no layer has vanished
390            !    => icount=5 : 5 layers have vanished
391            rswitch       = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) 
392            icount(ji,jk) = NINT( rswitch )
393            zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
394
395            ! update heat content (J.m-2) and layer thickness
396            qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
397            h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
398         END DO
399      END DO
400      ! update ice thickness
401      DO ji = kideb, kiut
402         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) )
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 iter = 1, num_iter_max
429         DO ji = kideb, kiut
430            IF(  zf_tt(ji) < 0._wp  ) THEN
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               ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
444
445               s_i_new(ji)        = zswitch_sal * zfracs * sss_m(ii,ij)  &  ! New ice salinity
446                                  + ( 1. - zswitch_sal ) * sm_i_1d(ji) 
447               ! New ice growth
448               ztmelts            = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K)
449
450               zt_i_new           = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
451               
452               zEi                = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)     
453                  &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   &
454                  &               + rcp  * ( ztmelts-rt0 )         
455
456               zEw                = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
457
458               zdE                = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
459
460               dh_i_bott(ji)      = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
461
462               q_i_1d(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0)
463               
464            ENDIF
465         END DO
466      END DO
467
468      ! Contribution to Energy and Salt Fluxes
469      DO ji = kideb, kiut
470         IF(  zf_tt(ji) < 0._wp  ) THEN
471            ! New ice growth
472                                   
473            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0)
474
475            ztmelts        = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K)
476           
477            zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
478           
479            zEi            = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)     
480               &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   &
481               &               + rcp  * ( ztmelts-rt0 )         
482           
483            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
484           
485            zdE            = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
486           
487            ! Contribution to heat flux to the ocean [W.m-2], >0 
488            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
489
490            ! Total heat flux used in this process [W.m-2], <0 
491            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
492           
493            ! Contribution to salt flux, <0
494            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice
495
496            ! Contribution to mass flux, <0
497            wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice
498
499            ! update heat content (J.m-2) and layer thickness
500            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)
501            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji)
502         ENDIF
503      END DO
504
505      !----------------
506      ! 4.2 Basal melt
507      !----------------
508      zdeltah(:,:) = 0._wp ! important
509      DO jk = nlay_i, 1, -1
510         DO ji = kideb, kiut
511            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting
512
513               ztmelts = - tmut * s_i_1d(ji,jk) + rt0  ! Melting point of layer jk (K)
514
515               IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
516
517                  zEi               = - q_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0)
518                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0)
519                                                                    ! set up at 0 since no energy is needed to melt water...(it is already melted)
520                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing     
521                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this
522
523                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk)
524
525                  zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0
526
527                  ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
528                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
529
530                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
531                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
532                                   
533                  ! Contribution to mass flux
534                  wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
535
536                  ! update heat content (J.m-2) and layer thickness
537                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
538                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
539
540               ELSE                               !!! Basal melting
541
542                  zEi             = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)
543                  zEw             = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0)
544                  zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0)
545
546                  zfmdt           = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0)
547
548                  zdeltah(ji,jk)  = - zfmdt * r1_rhoic         ! Gross thickness change
549
550                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
551                 
552                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors
553
554                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt
555
556                  zfmdt           = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0
557
558                  zQm             = zfmdt * zEw         ! Heat exchanged with ocean
559
560                  ! Contribution to heat flux to the ocean [W.m-2], <0 
561                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
562
563                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
564                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
565                 
566                  ! Total heat flux used in this process [W.m-2], >0 
567                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
568                 
569                  ! Contribution to mass flux
570                  wfx_bom_1d(ji)  =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
571
572                  ! update heat content (J.m-2) and layer thickness
573                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
574                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
575               ENDIF
576           
577            ENDIF
578         END DO
579      END DO
580
581      !-------------------------------------------
582      ! Update temperature, energy
583      !-------------------------------------------
584      DO ji = kideb, kiut
585         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_bott(ji) )
586      END DO 
587
588      !-------------------------------------------
589      ! 5. What to do with remaining energy
590      !-------------------------------------------
591      ! If heat still available for melting and snow remains, then melt more snow
592      !-------------------------------------------
593      zdeltah(:,:) = 0._wp ! important
594      DO ji = kideb, kiut
595         zq_rema(ji)     = zq_su(ji) + zq_bo(ji) 
596         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )   ! =1 if snow
597         rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,1) - epsi20 ) )
598         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 )
599         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting
600         dh_s_tot (ji)   = dh_s_tot(ji)  + zdeltah(ji,1)
601         ht_s_1d   (ji)  = ht_s_1d(ji)   + zdeltah(ji,1)
602       
603         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * q_s_1d(ji,1)                ! update available heat (J.m-2)
604         ! heat used to melt snow
605         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)
606         ! Contribution to mass flux
607         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
608         !   
609         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
610         ! Remaining heat flux (W.m-2) is sent to the ocean heat budget
611         hfx_out(ii,ij)  = hfx_out(ii,ij) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
612
613         IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
614      END DO
615     
616      !
617      !------------------------------------------------------------------------------|
618      !  6) Snow-Ice formation                                                       |
619      !------------------------------------------------------------------------------|
620      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
621      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
622      DO ji = kideb, kiut
623         !
624         dh_snowice(ji) = MAX(  0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic )  )
625
626         ht_i_1d(ji)    = ht_i_1d(ji) + dh_snowice(ji)
627         ht_s_1d(ji)    = ht_s_1d(ji) - dh_snowice(ji)
628
629         ! Salinity of snow ice
630         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
631         zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) * r1_rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji)
632
633         ! entrapment during snow ice formation
634         ! new salinity difference stored (to be used in limthd_sal.F90)
635         IF (  nn_icesal == 2  ) THEN
636            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )
637            ! salinity dif due to snow-ice formation
638            dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch     
639            ! salinity dif due to bottom growth
640            IF (  zf_tt(ji)  < 0._wp ) THEN
641               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch
642            ENDIF
643         ENDIF
644
645         ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
646         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
647         zfmdt          = ( rhosn - rhoic ) * MAX( dh_snowice(ji), 0._wp )    ! <0
648         zsstK          = sst_m(ii,ij) + rt0                               
649         zEw            = rcp * ( zsstK - rt0 )
650         zQm            = zfmdt * zEw 
651         
652         ! Contribution to heat flux
653         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
654
655         ! Contribution to salt flux
656         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_1d(ji) * zfmdt * r1_rdtice 
657         
658         ! Contribution to mass flux
659         ! All snow is thrown in the ocean, and seawater is taken to replace the volume
660         wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice
661         wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice
662
663         ! update heat content (J.m-2) and layer thickness
664         qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_1d(ji,1) + zfmdt * zEw
665         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
666         
667      END DO
668
669      !
670      !-------------------------------------------
671      ! Update temperature, energy
672      !-------------------------------------------
673      DO ji = kideb, kiut
674         rswitch     =  1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 
675         t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0
676      END DO
677
678      DO jk = 1, nlay_s
679         DO ji = kideb,kiut
680            ! mask enthalpy
681            rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) )  )
682            q_s_1d(ji,jk) = rswitch * q_s_1d(ji,jk)
683            ! recalculate t_s_1d from q_s_1d
684            t_s_1d(ji,jk) = rt0 + rswitch * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic )
685         END DO
686      END DO
687
688      ! --- ensure that a_i = 0 where ht_i = 0 ---
689      WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp
690     
691      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema )
692      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s, zsnw )
693      CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i )
694      CALL wrk_dealloc( jpij, nlay_i, icount )
695      !
696      !
697   END SUBROUTINE lim_thd_dh
698
699
700   !!--------------------------------------------------------------------------
701   !! INTERFACE lim_thd_snwblow
702   !! ** Purpose :   Compute distribution of precip over the ice
703   !!--------------------------------------------------------------------------
704   SUBROUTINE lim_thd_snwblow_2d( pin, pout )
705      REAL(wp), DIMENSION(:,:), INTENT(in)  :: pin   ! previous fraction lead ( pfrld or (1. - a_i_b) )
706      REAL(wp), DIMENSION(:,:), INTENT(out) :: pout
707      pout = ( 1._wp - ( pin )**rn_betas )
708   END SUBROUTINE lim_thd_snwblow_2d
709
710   SUBROUTINE lim_thd_snwblow_1d( pin, pout )
711      REAL(wp), DIMENSION(:), INTENT(in)  :: pin
712      REAL(wp), DIMENSION(:), INTENT(out) :: pout
713      pout = ( 1._wp - ( pin )**rn_betas )
714   END SUBROUTINE lim_thd_snwblow_1d
715
716   
717#else
718   !!----------------------------------------------------------------------
719   !!   Default option                               NO  LIM3 sea-ice model
720   !!----------------------------------------------------------------------
721CONTAINS
722   SUBROUTINE lim_thd_dh          ! Empty routine
723   END SUBROUTINE lim_thd_dh
724#endif
725
726   !!======================================================================
727END MODULE limthd_dh
Note: See TracBrowser for help on using the repository browser.