New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limthd_dh.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 5167

Last change on this file since 5167 was 5167, checked in by clem, 10 years ago

LIM3 bug fix. see ticket #1492 (forthcoming update) which also solve ticket #1497

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