source: trunk/NEMO/LIM_SRC/limthd_zdf.F90 @ 12

Last change on this file since 12 was 12, checked in by opalod, 17 years ago

CT : UPDATE001 : First major NEMO update

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 42.9 KB
Line 
1MODULE limthd_zdf
2#if defined key_ice_lim
3   !!======================================================================
4   !!                       ***  MODULE limthd_zdf ***
5   !!                thermodynamic growth and decay of the ice
6   !!======================================================================
7
8   !!----------------------------------------------------------------------
9   !!   lim_thd_zdf  : vertical accr./abl. and lateral ablation of sea ice
10   !! * Modules used
11   USE par_oce          ! ocean parameters
12   USE phycst           ! ???
13   USE ice_oce          ! ice variables
14   USE thd_ice
15   USE iceini
16   USE limistate
17   USE in_out_manager
18     
19   IMPLICIT NONE
20   PRIVATE
21
22   !! * Routine accessibility
23   PUBLIC lim_thd_zdf        ! called by lim_thd
24
25   !! * Module variables
26   REAL(wp)  ::           &  ! constant values
27      epsi20 = 1e-20   ,  &
28      epsi13 = 1e-13   ,  &
29      zzero  = 0.0     ,  &
30      zone   = 1.0
31   !!----------------------------------------------------------------------
32   !!   LIM 2.0 , UCL-LODYC-IPSL  (2003)
33   !!----------------------------------------------------------------------
34CONTAINS
35
36   SUBROUTINE lim_thd_zdf( kideb , kiut )
37       !!------------------------------------------------------------------
38       !!                ***  ROUTINE lim_thd_zdf  ***
39       !!             
40       !! ** Purpose : This routine determines the time evolution of snow
41       !!      and sea-ice thicknesses, concentration and heat content
42       !!      due to the vertical and lateral thermodynamic accretion-
43       !!      ablation processes. One only treats the case of lat. abl.
44       !!      For lateral accretion, see routine lim_lat_accr
45       !!
46       !! ** Method  : The representation of vertical growth and decay of
47       !!      the sea-ice model is based upon the diffusion of heat
48       !!      through the external and internal boundaries of a
49       !!      three-layer system (two layers of ice and one layer and
50       !!      one layer of snow, if present, on top of the ice).
51       !!
52       !! ** Action  : - Calculation of some intermediates variables
53       !!              - Calculation of surface temperature
54       !!              - Calculation of available heat for surface ablation
55       !!              - Calculation of the changes in internal temperature
56       !!                of the three-layer system, due to vertical diffusion
57       !!                processes
58       !!              - Performs surface ablation and bottom accretion-ablation
59       !!              - Performs snow-ice formation
60       !!              - Performs lateral ablation
61       !!
62       !! References :
63       !!   Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
64       !!   Fichefet T. and M. Maqueda 1999, Clim. Dyn, 15(4), 251-268 
65       !!
66       !! History :
67       !!   original    : 01-04 (LIM)
68       !!   addition    : 02-08 (C. Ethe, G. Madec)
69       !!------------------------------------------------------------------
70       !! * Arguments
71       INTEGER , INTENT (in) ::  &
72          kideb ,  &  ! Start point on which the  the computation is applied
73          kiut        ! End point on which the  the computation is applied
74
75       !! * Local variables
76       INTEGER ::   ji       ! dummy loop indices
77
78       REAL(wp) , DIMENSION(jpij,2) ::  &
79          zqcmlt        ! energy due to surface( /1 ) and bottom melting( /2 )
80
81       REAL(wp), DIMENSION(jpij) ::  &
82          ztsmlt      &    ! snow/ice surface melting temperature
83          ,ztbif      &    ! int. temp. at the mid-point of the 1st layer of the snow/ice sys.
84          ,zksn       &    ! effective conductivity of snow
85          ,zkic       &    ! effective conductivity of ice
86          ,zksndh     &    ! thermal cond. at the mid-point of the 1st layer of the snow/ice sys.
87          , zfcsu     &    ! conductive heat flux at the surface of the snow/ice system
88          , zfcsudt   &    ! = zfcsu * dt
89          , zi0       &    ! frac. of the net SW rad. which is not absorbed at the surface
90          , z1mi0     &    ! fraction of the net SW radiation absorbed at the surface
91          , zqmax     &    ! maximum energy stored in brine pockets
92          , zrcpdt    &    ! h_su*rho_su*cp_su/dt(h_su being the thick. of surf. layer)
93          , zts_old   &    ! previous surface temperature
94          , zidsn , z1midsn , zidsnic ! tempory variables
95
96       REAL(wp), DIMENSION(jpij) :: &
97          zfnet       &  ! net heat flux at the top surface( incl. conductive heat flux)
98          , zsprecip  &    ! snow accumulation
99          , zhsnw_old &    ! previous snow thickness
100          , zdhictop  &    ! change in ice thickness due to top surf ablation/accretion
101          , zdhicbot  &    ! change in ice thickness due to bottom surf abl/acc
102          , zqsup     &    ! energy transmitted to ocean (coming from top surf abl/acc)
103          , zqocea    &    ! energy transmitted to ocean (coming from bottom sur abl/acc)
104          , zfrl_old  &    ! previous sea/ice fraction
105          , zfrld_1d    &    ! new sea/ice fraction
106          , zep            ! internal temperature of the 2nd layer of the snow/ice system
107
108       REAL(wp), DIMENSION(3) :: & 
109          zplediag  &    ! principle diagonal, subdiag. and supdiag. of the
110          , zsubdiag  &    ! tri-diagonal matrix coming from the computation
111          , zsupdiag  &    ! of the temperatures inside the snow-ice system
112          , zsmbr          ! second member
113
114       REAL(wp) :: & 
115          zhsu     &     ! thickness of surface layer
116          , zhe      &     ! effective thickness for compu. of equ. thermal conductivity
117          , zheshth  &     ! = zhe / thth
118          , zghe     &     ! correction factor of the thermal conductivity
119          , zumsb    &     ! parameter for numerical method to solve heat-diffusion eq.
120          , zkhsn    &     ! conductivity at the snow layer
121          , zkhic    &     ! conductivity at the ice layers
122          , zkint    &     ! equivalent conductivity at the snow-ice interface
123          , zkhsnint &     ! = zkint*dt / (hsn*rhosn*cpsn) 
124          , zkhicint &     ! = 2*zkint*dt / (hic*rhoic*cpic)
125          , zpiv1 , zpiv2  &       ! tempory scalars used to solve the tri-diagonal system
126          , zb2 , zd2 , zb3 , zd3 &
127          , ztint          ! equivalent temperature at the snow-ice interface
128
129       REAL(wp) :: & 
130          zexp      &     ! exponential function of the ice thickness
131          , zfsab     &     ! part of solar radiation stored in brine pockets
132          , zfts      &     ! value of energy balance function when the temp. equal surf. temp.
133          , zdfts     &     ! value of derivative of ztfs when the temp. equal surf. temp.
134          , zdts      &     ! surface temperature increment
135          , zqsnw_mlt &     ! energy needed to melt snow
136          , zdhsmlt   &     ! change in snow thickness due to melt
137          , zhsn      &     ! snow thickness (previous+accumulation-melt)
138          , zqsn_mlt_rem &  ! remaining heat coming from snow melting
139          , zqice_top_mlt & ! energy used to melt ice at top surface
140          , zdhssub      &  ! change in snow thick. due to sublimation or evaporation
141          , zdhisub      &  ! change in ice thick. due to sublimation or evaporation   
142          , zdhsn        &  ! snow ice thickness increment
143          , zdtsn        &  ! snow internal temp. increment
144          , zdtic        &  ! ice internal temp. increment
145          , zqnes          ! conductive energy due to ice melting in the first ice layer
146
147       REAL(wp) :: & 
148          ztbot     &      ! temperature at the bottom surface
149          , zfcbot    &      ! conductive heat flux at bottom surface
150          , zqice_bot &      ! energy used for bottom melting/growing
151          , zqice_bot_mlt &  ! energy used for bottom melting
152          , zqstbif_bot  &  ! part of energy stored in brine pockets used for bottom melting
153          , zqstbif_old  &  ! tempory var. for zqstbif_bot
154          , zdhicmlt      &  ! change in ice thickness due to bottom melting
155          , zdhicm        &  ! change in ice thickness var.
156          , zdhsnm        &  ! change in snow thickness var.
157          , zhsnfi        &  ! snow thickness var.
158          , zc1, zpc1, zc2, zpc2, zp1, zp2 & ! tempory variables
159          , ztb2, ztb3
160
161       REAL(wp) :: & 
162          zdrmh         &   ! change in snow/ice thick. after snow-ice formation
163          , zhicnew       &   ! new ice thickness
164          , zhsnnew       &   ! new snow thickness
165          , zquot , ztneq &   ! tempory temp. variables
166          , zqice, zqicetot & ! total heat inside the snow/ice system
167          , zdfrl         &   ! change in ice concentration
168          , zdvsnvol      &   ! change in snow volume
169          , zdrfrl1, zdrfrl2 &  ! tempory scalars
170          , zihsn, zidhb, zihic, zihe, zihq, ziexp, ziqf, zihnf, zibmlt, ziqr, zihgnew, zind
171       !!----------------------------------------------------------------------
172
173       !-----------------------------------------------------------------------
174       !  1. Boundaries conditions for snow/ice system internal temperature
175       !       - If tbif_1d(ji,1) > rt0_snow, tbif_1d(ji,1) = rt0_snow
176       !       - If tbif_1d(ji,2/3) > rt0_ice, tbif_1d(ji,2/3) = rt0_ice
177       !     Computation of energies due to surface and bottom melting
178       !-----------------------------------------------------------------------
179       
180       
181       DO ji = kideb , kiut
182          zihsn = MAX( zzero , SIGN( zone , hsndif - h_snow_1d(ji) ) )
183          zihic = MAX( zzero , SIGN( zone , hicdif - h_ice_1d(ji) ) )
184          !--computation of energy due to surface melting
185          zqcmlt(ji,1) = ( MAX ( zzero ,  &
186             &                   rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn )
187          !--computation of energy due to bottom melting
188          zqcmlt(ji,2) = ( MAX( zzero , &
189             &                  rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) &
190             &           + MAX( zzero , &
191             &                  rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) &
192             &           ) * ( 1.0 - zihic  )
193          !--limitation of  snow/ice system internal temperature
194          tbif_1d(ji,1)   = MIN( rt0_snow, tbif_1d(ji,1) )
195          tbif_1d(ji,2)   = MIN( rt0_ice , tbif_1d(ji,2) )
196          tbif_1d(ji,3)   = MIN( rt0_ice , tbif_1d(ji,3) )
197       END DO
198
199
200       !-------------------------------------------
201       !  2. Calculate some intermediate variables. 
202       !-------------------------------------------
203       
204       ! initialisation of the thickness of surface layer
205       zhsu = hnzst 
206
207       DO ji = kideb , kiut
208          zind   = MAX( zzero , SIGN( zone , zhsu - h_snow_1d(ji) ) )
209          zihsn  = MAX( zzero , SIGN( zone , hsndif - h_snow_1d(ji) ) )
210          zihsn  = MAX( zihsn , zind )
211          zihic  = MAX( zzero , sign( zone , hicdif - h_ice_1d(ji) ) )
212          !     2.1. Computation of surface melting temperature
213          !----------------------------------------------------
214          zind  = MAX( zzero , SIGN( zone , -h_snow_1d(ji) ) )
215          ztsmlt(ji) = ( 1.0 - zind ) * rt0_snow + zind * rt0_ice
216          !
217          !     2.2. Effective conductivity of snow and ice
218          !-----------------------------------------------
219
220          !---computation of the correction factor on the thermal conductivity
221          !-- (Morales Maqueda, 1995 ; Fichefet and Morales Maqueda, 1997)
222          zhe      =  ( rcdsn / ( rcdsn + rcdic ) ) * h_ice_1d(ji)   &
223             &     + ( rcdic / ( rcdsn + rcdic ) ) * h_snow_1d(ji) 
224          zihe     = MAX( zzero , SIGN( zone , 2.0 * zhe - thth ) )
225          zheshth  = zhe / thth
226          zghe     = ( 1.0 - zihe ) * zheshth * ( 2.0 - zheshth )   &
227             &     +         zihe   * 0.5 * ( 1.5 + LOG( 2.0 * zheshth ) )
228#if defined key_lim_cp3
229          zghe = 1.0
230#endif 
231
232          !---effective conductivities
233          zksn(ji)  = zghe * rcdsn 
234          zkic(ji)  = zghe * rcdic
235
236          !
237          !     2.3. Computation of the conductive heat flux from the snow/ice
238          !          system interior toward the top surface
239          !------------------------------------------------------------------
240
241          !---Thermal conductivity at the mid-point of the first snow/ice system layer
242          zksndh(ji) =   ( ( 1.0 - zihsn ) * 2.0 * zksn(ji) + zihsn * 4.0 * zkic(ji) )   &
243             &         / ( ( 1.0 - zihsn ) *  h_snow_1d(ji)                              &
244             &           +        zihsn   *  ( ( 1.0 + 3.0 * zihic ) * h_ice_1d(ji)      &
245             &           + 4.0 * zkic(ji)/zksn(ji) * h_snow_1d(ji) ) )
246
247          !---internal temperature at the mid-point of the first snow/ice system layer
248          ztbif(ji)  = ( 1.0 - zihsn ) * tbif_1d(ji,1)                       &
249             &       +         zihsn   * ( ( 1.0 - zihic ) * tbif_1d(ji,2)   &
250             &       +         zihic   * tfu_1d(ji)   )
251          !---conductive heat flux
252          zfcsu(ji) = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) )
253
254       END DO
255
256       !--------------------------------------------------------------------
257       !  3. Calculate :
258       !     - fstbif_1d, part of solar radiation absorbing inside the ice
259       !       assuming an exponential absorption (Grenfell and Maykut, 1977)
260       !     - zqmax,  maximum energy stored in brine pockets
261       !     - qstbif_1d, total energy stored in brine pockets (updating)
262       !-------------------------------------------------------------------
263
264
265       DO ji = kideb , kiut
266          zihsn  = MAX( zzero , SIGN (zone , -h_snow_1d(ji) ) )
267          zihic  = MAX( zzero , 1.0 - ( h_ice_1d(ji) / zhsu ) )     
268          zind   = MAX( zzero , SIGN (zone , hicdif - h_ice_1d(ji) ) )
269          !--Computation of the fraction of the net shortwave radiation which
270          !--penetrates inside the ice cover ( See Forcat)
271          zi0(ji)  = zihsn * ( fr1_i0_1d(ji) + zihic * fr2_i0_1d(ji) )
272          zexp     = MIN( zone , EXP( -1.5 * ( h_ice_1d(ji) - zhsu ) ) )
273          fstbif_1d(ji) = zi0(ji) * qsr_ice_1d(ji) * zexp
274          !--Computation of maximum energy stored in brine pockets zqmax and update
275          !--the total energy stored in brine pockets, if less than zqmax
276          zqmax(ji) = MAX( zzero , 0.5 * xlic * ( h_ice_1d(ji) - hicmin ) )
277          zfsab   = zi0(ji) * qsr_ice_1d(ji) * ( 1.0 - zexp )
278          zihq    = ( 1.0 - zind ) * MAX(zzero, SIGN( zone , qstbif_1d(ji) - zqmax(ji) ) ) &
279             &    +         zind   * zone
280          qstbif_1d(ji) = ( qstbif_1d(ji) + ( 1.0 - zihq ) * zfsab * rdt_ice ) * swiqst
281          !--fraction of shortwave radiation absorbed at surface
282          ziexp = zihq * zexp + ( 1.0 - zihq ) * ( swiqst + ( 1.0 - swiqst ) * zexp )
283          z1mi0(ji) = 1.0 - zi0(ji) * ziexp
284       END DO
285
286
287       !--------------------------------------------------------------------------------
288       !  4. Computation of the surface temperature : determined by considering the
289       !     budget of a thin layer of thick. zhsu at the top surface (H. Grenier, 1995)
290       !     and based on a surface energy balance :
291       !     hsu * rcp * dT/dt = Fsr + Fnsr(T) + Fcs(T),
292       !     where - Fsr is the net absorbed solar radiation,
293       !           - Fnsr is the total non solar radiation (incoming and outgoing long-wave,
294       !             sensible and latent heat fluxes)
295       !           - Fcs the conductive heat flux at the top of surface
296       !------------------------------------------------------------------------------
297
298       !     4.1. Computation of intermediate values
299       !---------------------------------------------
300       DO ji = kideb, kiut
301          zrcpdt(ji) = ( rcpsn * MIN( h_snow_1d(ji) , zhsu )    &
302             &       + rcpic * MAX( zhsu - h_snow_1d(ji) , zzero ) ) / rdt_ice
303          zts_old(ji) =  sist_1d(ji)
304       END DO
305
306       !     4.2. Computation of surface temperature by expanding the eq. of energy balance
307       !          with Ts = Tp + DT. One obtain , F(Tp) + DT * DF(Tp) = 0
308       !          where  - F(Tp) = Fsr + Fnsr(Tp) + Fcs(Tp)
309       !                 - DF(Tp)= (dFnsr(Tp)/dT) + (dFcs(Tp)/dT) - hsu*rcp/dt
310       !---------------------------------------------------------------------------------
311
312       DO ji = kideb, kiut
313          !---computation of the derivative of energy balance function
314#if defined key_coupled
315#   if defined key_lim_cp2
316          zdfts   =   zksndh(ji)   & ! contribution of the conductive heat flux
317             &      + zrcpdt(ji)   & ! contribution of hsu * rcp / dt
318             &      - dqns_ice_1d(ji)      ! contribution of the total non solar radiation
319#   else
320          zdfts   =   zksndh(ji)   & ! contribution of the conductive heat flux
321             &      + zrcpdt(ji)    ! contribution of hsu * rcp / dt
322#   endif
323
324#else
325          zdfts    =  zksndh(ji)   & ! contribution of the conductive heat flux
326             &      + zrcpdt(ji)   & ! contribution of hsu * rcp / dt
327             &      - dqns_ice_1d (ji)     ! contribution of the total non solar radiation
328#endif
329          !---computation of the energy balance function
330          zfts    = - z1mi0 (ji) * qsr_ice_1d(ji)   & ! net absorbed solar radiation
331             &      - qnsr_ice_1d(ji)                & ! total non solar radiation
332             &      - zfcsu (ji)                  ! conductive heat flux from the surface
333          !---computation of surface temperature increment 
334          zdts    = -zfts / zdfts
335#if defined key_lim_cp3
336          zdts = zdts / 3.0
337#endif
338          !---computation of the new surface temperature
339          sist_1d(ji) = sist_1d(ji) + zdts
340
341       END DO
342
343       !----------------------------------------------------------------------------
344       !  5. Boundary condition at the top surface
345       !--    IF Tsb < Tmelt, Fnet = Fcs (the net heat flux equal the conductive heat flux)
346       !      Otherwise Tsb = Tmelt and Qnet(Tmelt) > 0
347       !      Fnet(Tmelt) is therefore the net surface flux needed for melting
348       !----------------------------------------------------------------------------
349       
350       
351       !     5.1.  Limitation of surface temperature and update total non solar fluxes,
352       !          latent heat flux and conductive flux at the top surface
353       !---------------------------------------------------------------------- 
354                     
355       DO ji = kideb, kiut
356          sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) )
357#if defined key_coupled
358          qnsr_ice_1d(ji) =  qnsr_ice_1d(ji) - ( emic * stefan * sist_1d(ji) * sist_1d(ji) * sist_1d(ji) * sist_1d(ji) )
359#else         
360          qnsr_ice_1d(ji) = qnsr_ice_1d(ji) + dqns_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) )
361          qla_ice_1d (ji) = qla_ice_1d (ji) + dqla_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) )
362#endif
363          zfcsu(ji)  = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) )
364       END DO
365
366     
367
368       !     5.2. Calculate available heat for surface ablation.
369       !---------------------------------------------------------------------
370
371       DO ji = kideb, kiut
372          zfnet(ji) = qnsr_ice_1d(ji) + z1mi0(ji) * qsr_ice_1d(ji) + zfcsu(ji)         
373          zfnet(ji) = MAX( zzero , zfnet(ji) )
374          zfnet(ji) = zfnet(ji) * MAX( zzero , SIGN( zone , sist_1d(ji) - ztsmlt(ji) ) )
375       END DO
376
377       !-------------------------------------------------------------------------
378       !  6. Calculate changes in internal temperature due to vertical diffusion   
379       !     processes. The evolution of this temperature is governed by the one-
380       !     dimensionnal heat-diffusion equation.
381       !     Given the temperature tbif(1/2/3), at time m we solve a set
382       !     of finite difference equations to obtain new tempe. Each tempe is coupled
383       !     to the temp. immediatly above and below by heat conduction terms. Thus
384       !     we have a set of equations of the form A * T = B, where A is a tridiagonal
385       !     matrix, T a vector whose components are the unknown new temp.
386       !-------------------------------------------------------------------------
387       
388       !--parameter for the numerical methode use to solve the heat-diffusion equation
389       !- implicit, explicit or Crank-Nicholson
390       zumsb = 1.0 - sbeta 
391       DO ji = kideb, kiut
392          zidsn(ji)   = MAX ( zzero, SIGN( zone, hsndif - h_snow_1d(ji) ) ) 
393          z1midsn(ji) = 1.0 - zidsn(ji)
394          zihic       = MAX ( zzero, SIGN( zone, hicdif - h_ice_1d(ji) ) ) 
395          zidsnic(ji) = zidsn(ji) *  zihic 
396          zfcsudt(ji) = zfcsu(ji) * rdt_ice 
397       END DO
398   
399       DO ji = kideb, kiut
400
401          !     6.1 Calculate intermediate variables.
402          !----------------------------------------
403
404          !--conductivity at the snow surface
405          zkhsn = 2.0 * zksn(ji) * rdt_ice / rcpsn
406          !--conductivity at the ice surface
407          zkhic = 4.0 * zkic(ji) * rdt_ice / MAX( h_ice_1d(ji) * h_ice_1d(ji) * rcpic , epsi20 )
408          !--conductivity at the snow/ice interface
409          zkint = 4.0 * zksn(ji) * zkic(ji)  &
410             &        / ( zksn(ji) * h_ice_1d(ji) + 2.0 * zkic(ji) * h_snow_1d(ji) * z1midsn(ji)) 
411          zkhsnint = zkint * rdt_ice / rcpsn
412          zkhicint = zkint * 2.0 * rdt_ice / MAX( h_ice_1d(ji) * rcpic , epsi20 )
413         
414          !     6.2. Fulfill the linear system matrix.
415          !-----------------------------------------
416!$$$          zplediag(1) = 1 + sbeta * z1midsn(ji) * ( zkhsn + zkhsnint )       
417          zplediag(1) =   zidsn(ji) + z1midsn(ji) * h_snow_1d(ji)   &
418             &          + sbeta * z1midsn(ji) * zkhsnint 
419          zplediag(2) = 1 + sbeta * ( z1midsn(ji) * zkhicint + zkhic ) 
420          zplediag(3) = 1 + 3.0 * sbeta * zkhic   
421
422          zsubdiag(1) = 0.0             
423          zsubdiag(2) = -1.0 * z1midsn(ji) * sbeta * zkhicint
424          zsubdiag(3) = -1.0 * sbeta * zkhic 
425
426          zsupdiag(1) =  -1.0 * z1midsn(ji) * sbeta * zkhsnint 
427          zsupdiag(2) = zsubdiag(3)
428          zsupdiag(3) = 0.0
429         
430          !     6.3. Fulfill the idependent term vector.
431          !-------------------------------------------
432         
433!$$$          zsmbr(1) = zidsn(ji) * sist_1d(ji) + z1midsn(ji) *   &
434!$$$             &         ( tbif_1d(ji,1) + zkhsn * sist_1d(ji)
435!$$$             &         - zumsb * ( zkhsn * tbif_1d(ji,1)
436!$$$             &                   + zkhsnint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) ) )
437          zsmbr(1) = zidsn(ji) * sist_1d(ji) + z1midsn(ji) *    &
438             &       ( h_snow_1d(ji) * tbif_1d(ji,1) - ( zfcsudt(ji) / rcpsn )  &
439             &       - zumsb * zkhsnint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) )
440
441          zsmbr(2) =  tbif_1d(ji,2)  &
442             &      - zidsn(ji) * ( 1.0 - zidsnic(ji) ) &
443             &        * ( zfcsudt(ji) / MAX( h_ice_1d(ji) * rcpic , epsi20 ) ) &
444             &      + zumsb * ( zkhicint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) &
445             &                   - zkhic * ( tbif_1d(ji,2) - tbif_1d(ji,3) )  )
446
447          zsmbr(3) =  tbif_1d(ji,3)  &
448             &      + zkhic * ( 2.0 * tfu_1d(ji) &
449             &                + zumsb * ( tbif_1d(ji,2) - 3.0 * tbif_1d(ji,3) ) )
450         
451          !     6.4. Solve the system (Gauss elimination method).
452          !----------------------------------------------------
453         
454          zpiv1 = zsubdiag(2) / zplediag(1) 
455          zb2   = zplediag(2) - zpiv1 * zsupdiag(1)
456          zd2   = zsmbr(2) - zpiv1 * zsmbr(1)
457
458          zpiv2 = zsubdiag(3) / zb2
459          zb3   = zplediag(3) - zpiv2 * zsupdiag(2)
460          zd3   = zsmbr(3) - zpiv2 * zd2
461
462          tbif_1d(ji,3) = zd3 / zb3
463          tbif_1d(ji,2) = ( zd2 - zsupdiag(2) * tbif_1d(ji,3) ) / zb2
464          tbif_1d(ji,1) = ( zsmbr(1) - zsupdiag(1) * tbif_1d(ji,2) ) / zplediag(1)           
465
466          !--- taking into account the particular case of  zidsnic(ji) = 1
467          ztint =  (  zkic(ji) * h_snow_1d(ji) * tfu_1d (ji)    &
468             &      + zksn(ji) * h_ice_1d(ji) * sist_1d(ji) )   &
469             &   / ( zkic(ji) * h_snow_1d(ji) + zksn(ji) * h_ice_1d(ji) ) 
470
471          tbif_1d(ji,1) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,1)   &
472             &                + zidsnic(ji)   * ( ztint + sist_1d(ji) ) / 2.0
473          tbif_1d(ji,2) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,2)   &
474             &                + zidsnic(ji)   * ( 3.0 * ztint + tfu_1d(ji) ) / 4.0
475          tbif_1d(ji,3) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,3)   &
476             &                + zidsnic(ji)   * ( ztint + 3.0 * tfu_1d(ji) ) / 4.0     
477       END DO
478 
479       !----------------------------------------------------------------------
480       !  9. Take into account surface ablation and bottom accretion-ablation.|
481       !----------------------------------------------------------------------
482       
483       !---Snow accumulation in one thermodynamic time step
484       zsprecip(kideb:kiut) = sprecip_1d(kideb:kiut) * rdt_ice / rhosn
485
486
487       DO ji = kideb, kiut
488         
489          !      9.1. Surface ablation and update of snow thickness and qstbif_1d
490          !--------------------------------------------------------------------
491         
492          !--------------------------------------------------------------------------
493          !--      Melting snow processes :
494          !--      Melt at the upper surface is computed from the difference between
495          !--      the net heat flux (including the conductive heat flux) at the upper
496          !--      surface and the pre-existing energy due to surface melting
497          !------------------------------------------------------------------------------
498         
499          !-- store the snow thickness
500          zhsnw_old(ji) =  h_snow_1d(ji)
501          !--computation of the energy needed to melt snow
502          zqsnw_mlt  = zfnet(ji) * rdt_ice - zqcmlt(ji,1)
503          !--change in snow thickness due to melt
504          zdhsmlt = - zqsnw_mlt / xlsn
505         
506          !-- compute new snow thickness, taking into account the part of snow accumulation
507          !   (as snow precipitation) and the part of snow lost due to melt
508          zhsn =  h_snow_1d(ji) + zsprecip(ji) + zdhsmlt
509          h_snow_1d(ji) = MAX( zzero , zhsn )
510          !-- compute the volume of snow lost after surface melting and the associated mass
511          dvsbq_1d(ji) =  ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnw_old(ji) - zsprecip(ji) )
512          dvsbq_1d(ji) =  MIN( zzero , dvsbq_1d(ji) )
513          rdmsnif_1d(ji) =  rhosn * dvsbq_1d(ji)
514          !-- If the snow is completely melted the remaining heat is used to melt ice
515          zqsn_mlt_rem  = MAX( zzero , -zhsn ) * xlsn
516          zqice_top_mlt = zqsn_mlt_rem 
517          zqstbif_old   = qstbif_1d(ji)
518
519          !--------------------------------------------------------------------------
520          !--      Melting ice processes at the top surface :
521          !--      The energy used to melt ice, zqice_top_mlt, is taken from the energy
522          !--      stored in brine pockets qstbif_1d and the remaining energy coming
523          !--      from the melting snow process zqsn_mlt_rem.
524          !--      If qstbif_1d > zqsn_mlt_rem then, one uses only a zqsn_mlt_rem part
525          !--      of qstbif_1d to melt ice,
526          !--         zqice_top_mlt = zqice_top_mlt + zqsn_mlt_rem
527          !--         qstbif_1d = qstbif_1d - zqsn_mlt_rem
528          !--      Otherwise one uses all qstbif_1d to melt ice
529          !--         zqice_top_mlt = zqice_top_mlt + qstbif_1d
530          !--         qstbif_1d = 0
531          !------------------------------------------------------
532         
533          ziqf =  MAX ( zzero , SIGN( zone , qstbif_1d(ji) - zqsn_mlt_rem  ) )
534          zqice_top_mlt =         ziqf   * ( zqice_top_mlt + zqsn_mlt_rem )   &
535             &          + ( 1.0 - ziqf ) * ( zqice_top_mlt + qstbif_1d(ji)  )
536
537          qstbif_1d(ji) =         ziqf   * ( qstbif_1d(ji) - zqsn_mlt_rem )   &
538             &        + ( 1.0 - ziqf ) * ( qstbif_1d(ji) - qstbif_1d(ji)  )
539
540          !--    The contribution of the energy stored in brine pockets qstbif_1d to melt
541          !--    ice is taking into account only when qstbif_1d is less than zqmax.
542          !--    Otherwise, only the remaining energy coming from the melting snow
543          !--    process is used
544          zihq =  MAX ( zzero , SIGN( zone , qstbif_1d(ji) - zqmax(ji) ) )
545
546          zqice_top_mlt =         zihq   * zqice_top_mlt   &
547             &          + ( 1.0 - zihq ) * zqsn_mlt_rem
548
549          qstbif_1d(ji) =         zihq   * qstbif_1d(ji)   &
550             &        + ( 1.0 - zihq ) * zqstbif_old
551
552          !--change in ice thickness due to melt at the top surface
553          zdhictop(ji) = -zqice_top_mlt / xlic
554          !--compute the volume formed after surface melting
555          dvsbq_1d(ji) =  zdhictop(ji) * ( 1.0 - frld_1d(ji) )
556
557          !-------------------------------------------------------------------------
558          !--      A small variation at the surface also occurs because of sublimation
559          !--      associated with the latent flux. If qla_ice_1d is negative, snow condensates at
560          !        the surface. Otherwise, snow evaporates
561          !-----------------------------------------------------------------------
562          !----change in snow and ice thicknesses due to sublimation or evaporation
563          zdhssub  = parsub * ( qla_ice_1d(ji) / ( rhosn * xsn ) ) * rdt_ice 
564          zhsn     = h_snow_1d(ji) - zdhssub
565          zdhisub  = MAX( zzero , -zhsn ) * rhosn/rhoic
566          zdhictop(ji) =  zdhictop(ji) - zdhisub
567          h_snow_1d(ji)  =  MAX( zzero , zhsn )
568          !-------------------------------------------------
569          !--  Update Internal temperature and qstbif_1d.
570          !-------------------------------------------
571          zihsn  =  MAX( zzero , SIGN( zone, -h_snow_1d(ji) ) )
572          tbif_1d(ji,1) = ( 1.0 - zihsn ) * tbif_1d(ji,1) + zihsn   * tfu_1d(ji)
573          !--change in snow internal temperature if snow has increased
574          zihnf = MAX( zzero , SIGN( zone , h_snow_1d(ji) - zhsnw_old(ji) ) )
575          zdhsn = 1.0 - zhsnw_old(ji) / MAX( h_snow_1d(ji) , epsi20 )
576          zdtsn = zdhsn * ( sist_1d(ji) - tbif_1d(ji,1) )
577          tbif_1d(ji,1) = tbif_1d(ji,1) + z1midsn(ji) * zihnf * zdtsn
578          !--energy created due to ice melting in the first ice layer
579          zqnes  = ( rt0_ice - tbif_1d(ji,2) ) * rcpic * ( h_ice_1d(ji) / 2. )
580          !--change in first ice layer internal temperature
581          ziqr  = MAX( zzero , SIGN( zone , qstbif_1d(ji) - zqnes ) )
582          zdtic = qstbif_1d(ji) / ( rcpic * ( h_ice_1d(ji) / 2. ) )
583          tbif_1d(ji,2) =  ziqr * rt0_ice + ( 1 - ziqr ) * ( tbif_1d(ji,2) + zdtic )
584          !--update qstbif_1d
585          qstbif_1d(ji) = ziqr * ( qstbif_1d(ji) - zqnes ) * swiqst
586
587
588          !--      9.2. Calculate bottom accretion-ablation and update qstbif_1d.
589          !             Growth and melting at bottom ice surface are governed by 
590          !                 -xlic * Dh = (Fcb - Fbot ) * Dt
591          !             where Fbot is the net downward heat flux from ice to the ocean
592          !            and Fcb is the conductive heat flux at the bottom surface
593          !---------------------------------------------------------------------------
594          ztbot = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,3) + zidsnic(ji) * sist_1d(ji)
595          !---computes conductive heat flux at bottom surface
596          zfcbot =  4.0 * zkic(ji) * ( tfu_1d(ji) - ztbot )   &
597             &   / ( h_ice_1d(ji) + zidsnic(ji) * ( 3. * h_ice_1d(ji) &
598             &   + 4.0 * zkic(ji)/zksn(ji) * h_snow_1d(ji) ) )
599          !---computation of net energy needed for bottom melting/growing
600          zqice_bot = ( zfcbot - ( fbif_1d(ji) + qlbbq_1d(ji) ) ) * rdt_ice
601          zqstbif_bot = qstbif_1d(ji)
602          !---switch to know if bottom surface melts ( = 1 ) or grows ( = 0 )occurs
603          zibmlt = MAX( zzero , SIGN( zone , -zqice_bot ) )
604          !--particular case of melting (in the same way as the top surface)
605          zqice_bot_mlt = zqice_bot 
606          zqstbif_old = zqstbif_bot
607
608          ziqf =  MAX ( zzero , SIGN( zone , qstbif_1d(ji) + zqice_bot_mlt  ) )
609          zqice_bot_mlt =         ziqf   * ( zqice_bot_mlt + zqice_bot_mlt ) &
610             &          + ( 1.0 - ziqf ) * ( zqice_bot_mlt + qstbif_1d(ji)  )         
611          qstbif_1d(ji)   =         ziqf   * ( qstbif_1d(ji) + zqice_bot_mlt ) &
612             &          + ( 1.0 - ziqf ) * ( qstbif_1d(ji) - qstbif_1d(ji)  )
613          !--    The contribution of the energy stored in brine pockets qstbif_1d to melt
614          !--    ice is taking into account only when qstbif_1d is less than zqmax.
615          zihq =  MAX ( zzero , SIGN( zone , qstbif_1d(ji) - zqmax(ji) ) )
616          zqice_bot_mlt =         zihq   * zqice_bot_mlt   &
617             &          + ( 1.0 - zihq ) * zqice_bot
618          qstbif_1d(ji)   =         zihq   * qstbif_1d(ji)   &
619             &             + ( 1.0 - zihq ) * zqstbif_old
620
621          !---treatment of the case of melting/growing
622          zqice_bot   =         zibmlt   * ( zqice_bot_mlt - zqcmlt(ji,2) )   &
623             &        + ( 1.0 - zibmlt ) * ( zqice_bot - zqcmlt(ji,2)  )
624          qstbif_1d(ji) =         zibmlt   * qstbif_1d(ji)   &
625             &           + ( 1.0 - zibmlt ) * zqstbif_bot
626
627          !--computes change in ice thickness due to melt or growth
628          zdhicbot(ji) = zqice_bot / xlic
629          !--limitation of bottom melting if so : hmelt maximum melting at bottom
630          zdhicmlt  = MAX( hmelt , zdhicbot(ji) ) 
631          !--energy after bottom melting/growing
632          zqsup(ji) = ( 1.0 - frld_1d(ji) ) * xlic * ( zdhicmlt - zdhicbot(ji) )
633          !-- compute the new thickness and the newly formed volume after bottom melting/growing
634          zdhicbot(ji)  = zdhicmlt
635          dvbbq_1d(ji) = ( 1.0 - frld_1d(ji) ) * zdhicbot(ji)
636
637
638          !        9.3.  Updating ice thickness after top surface ablation
639          !              and bottom surface accretion/ablation
640          !---------------------------------------------------------------
641          zhicnew  = h_ice_1d(ji) + zdhictop(ji) + zdhicbot(ji)
642
643          !
644          !        9.4. Case of total ablation (ice is gone but snow may be left)
645          !-------------------------------------------------------------------
646          zhsn  = h_snow_1d(ji)
647          zihgnew = 1.0 - MAX( zzero , SIGN( zone , -zhicnew ) )
648          zihsn   = MAX( zzero , SIGN( zone , -zhsn ) )
649          !---convert
650          zdhicm  = ( 1.0 - zihgnew ) * ( zhicnew - qstbif_1d(ji) / xlic )
651          zdhsnm  = ( 1.0 - zihsn ) * zdhicm * rhoic / rhosn
652          !---updating new ice thickness and computing the newly formed ice mass
653          zhicnew   =  zihgnew * zhicnew
654          rdmicif_1d(ji) =  rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) ) * rhoic
655          !---updating new snow thickness and computing the newly formed snow mass
656          zhsnfi   = zhsn + zdhsnm
657          h_snow_1d(ji) = MAX( zzero , zhsnfi )
658          rdmsnif_1d(ji) =  rdmsnif_1d(ji) + ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsn ) * rhosn
659          !--remaining energy in case of total ablation
660          zqocea(ji) = - ( zihsn * xlic * zdhicm + xlsn * ( zhsnfi - h_snow_1d(ji) ) ) * ( 1.0 - frld_1d(ji) )
661          qstbif_1d(ji) = zihgnew * qstbif_1d(ji)
662
663          !
664          !        9.5. Update internal temperature and ice thickness.
665          !-------------------------------------------------------
666          !
667          sist_1d(ji) = zihgnew * sist_1d(ji) + ( 1.0 - zihgnew ) * tfu_1d(ji)
668          zidhb  = MAX( zzero , SIGN( zone , - zdhicbot(ji) ) )
669          zc1    = - zhicnew * 0.5
670          zpc1   = MIN( 0.5 * zone , - h_ice_1d(ji) * 0.5 - zdhictop(ji) )
671          zc2    = - zhicnew
672          zpc2   =  zidhb * zc2 + ( 1.0 - zidhb ) * ( - h_ice_1d(ji) - zdhictop(ji) )
673          zp1    =  MAX( zpc1 , zc1 )
674          zp2    =  MAX( zpc2 , zc1 )
675          zep(ji) =  tbif_1d(ji,2)
676          ztb2  = 2.0 * (         - zp1   * tbif_1d(ji,2)  &
677             &  + ( zp1 - zp2 ) * tbif_1d(ji,3)  &
678             &  + ( zp2 - zc1 ) * tfu_1d(ji) ) / MAX( zhicnew , epsi20 ) 
679          tbif_1d(ji,2) = zihgnew * ztb2 + ( 1.0 - zihgnew ) * tfu_1d(ji)
680          !---
681          zp1  =  MIN( zpc1 , zc1 )
682          zp2  =  MIN( zpc2 , zc1 )
683          zp1  =  MAX( zc2  , zp1 )
684          ztb3 =  2.0 * (   ( 1.0 - zidhb ) * (  ( zc1 - zp2 ) * tbif_1d(ji,3)  &
685             &                                 + ( zp2 - zc2 ) * tfu_1d(ji) )   &
686             &               +      zidhb   * (  ( zc1 - zp1 ) * zep(ji)      &
687             &                                 + ( zp1 - zc2 ) * tbif_1d(ji,3))  ) / MAX( zhicnew , epsi20 )
688          tbif_1d(ji,3) =  zihgnew * ztb3 + ( 1.0 - zihgnew ) * tfu_1d(ji)
689          h_ice_1d(ji)  =  zhicnew
690       END DO
691
692
693       !----------------------------------------------------------------------------
694       !  10. Surface accretion.
695       !      The change of ice thickness after snow/ice formation is such that
696       !      the interface between snow and ice is located at the same height
697       !      as the ocean surface. It is given by (Fichefet and Morales Maqueda 1999)
698       !          D(h_ice) = (- D(hsn)/alph) =  [rhosn*hsn - (rau0 - rhoic)*hic]
699       !                                     / [alph*rhosn+rau0 - rhoic]
700       !----------------------------------------------------------------------------
701       !
702       DO ji = kideb , kiut
703
704          !--  Computation of the change of ice thickness after snow-ice formation
705          zdrmh =  ( rhosn * h_snow_1d(ji) + ( rhoic - rau0 ) * h_ice_1d(ji) )  &
706             &  / ( alphs * rhosn + rau0 - rhoic )
707          zdrmh = MAX( zzero , zdrmh )
708
709          !--New ice and snow thicknesses Fichefet and Morales Maqueda (1999)
710          zhicnew  = MAX( h_ice_1d(ji) , h_ice_1d(ji) + zdrmh )
711          zhsnnew  = MIN( h_snow_1d(ji) , h_snow_1d(ji) - alphs * zdrmh )
712          !---Compute new ice temperatures. snow temperature remains unchanged
713          !   Lepparanta (1983):
714          zihic = 1.0 - MAX( zzero , SIGN( zone , -zhicnew ) )
715          zquot  = ( 1.0 - zihic ) &
716             &   +         zihic * MIN( zone , h_ice_1d(ji) / MAX( zhicnew , epsi20 ) ) 
717          ztneq  =         alphs * cnscg * tbif_1d(ji,1)    &
718             &   + ( 1.0 - alphs * ( rhosn/rhoic ) ) * tfu_1d(ji)
719          zep(ji) = tbif_1d(ji,2)
720          tbif_1d(ji,2) = ztneq - zquot * zquot * ( ztneq - tbif_1d(ji,2) )
721          tbif_1d(ji,3) = 2.0 * ztneq &
722             &        + zquot * ( tbif_1d(ji,3) + zep(ji) - 2.0 * ztneq ) - tbif_1d(ji,2)
723
724          !---  Lepparanta (1983) (latent heat released during white ice formation
725          !     goes to the ocean -for lateral ablation-)
726          qldif_1d(ji)  = qldif_1d(ji) + zdrmh * ( 1.0 - alphs * ( rhosn/rhoic ) ) * xlic * ( 1.0 - frld_1d(ji) )
727          !--   Changes in ice volume and ice mass Lepparanta (1983):
728          dvnbq_1d(ji) = ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) )
729          dmgwi_1d(ji) = dmgwi_1d(ji) + ( 1.0 -frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnnew ) * rhosn
730          !  case of natural freshwater flux
731#if defined key_lim_fdd 
732          rdmicif_1d(ji) = rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhicnew - h_ice_1d(ji) )  * rhoic
733          rdmsnif_1d(ji) = rdmsnif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhsnnew - h_snow_1d(ji) ) * rhosn
734
735#else
736          rdmicif_1d(ji) = rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) ) * (  ( zhicnew - h_ice_1d(ji) ) * rhoic   &
737             &                                                    + ( zhsnnew - h_snow_1d(ji) ) * rhosn )
738#endif
739
740          !---  Actualize new snow and ice thickness.
741
742          h_snow_1d(ji)  = zhsnnew
743          h_ice_1d(ji)  = zhicnew
744
745       END DO
746
747       !----------------------------------------------------
748       !  11. Lateral ablation (Changes in sea/ice fraction)
749       !----------------------------------------------------
750       DO ji = kideb , kiut
751          zfrl_old(ji)  = frld_1d(ji)
752          zihic   = 1.0 - MAX( zzero , SIGN( zone , -h_ice_1d(ji) ) )
753          zihsn   = 1.0 - MAX( zzero , SIGN( zone , -h_snow_1d(ji) ) )
754          !--In the case of total ablation (all the ice ice has melted) frld = 1
755          frld_1d(ji)  = ( 1.0 - zihic ) + zihic * zfrl_old(ji)
756          !--Part of solar radiation absorbing inside the ice and going
757          !--through the ocean
758          fscbq_1d(ji) = ( 1.0 - zfrl_old(ji) ) * ( 1.0 - thcm_1d(ji) ) * fstbif_1d(ji)
759          !--Total remaining energy after bottom melting/growing
760          qfvbq_1d(ji) = zqsup(ji) + ( 1.0 - zihic ) * zqocea(ji)
761          !--Updating of total heat from the ocean
762          qldif_1d(ji)  = qldif_1d(ji) + qfvbq_1d(ji) + ( 1.0 - zihic ) * fscbq_1d(ji) * rdt_ice
763          !--Computation of total heat inside the snow/ice system
764          zqice  = h_snow_1d(ji) * xlsn + h_ice_1d(ji) * xlic
765          zqicetot  = ( 1.0 - frld_1d(ji) ) * zqice
766          !--The concentration of ice is reduced (frld increases) if the heat
767          !--exchange between ice and ocean is positive
768          ziqf = MAX( zzero , SIGN( zone ,  zqicetot - qldif_1d(ji) ) )
769          zdfrl = qldif_1d(ji) / MAX( epsi20 , zqice ) 
770          frld_1d(ji)  = ( 1.0 - ziqf )    &
771             &       +         ziqf * ( frld_1d(ji) + MAX( zzero , zdfrl ) ) 
772          fltbif_1d(ji) = ( ( 1.0 - zfrl_old(ji) ) * qstbif_1d(ji) - zqicetot  ) / rdt_ice
773          !--  Opening of leads: Hakkinen & Mellor, 1992.
774          zdfrl = - ( zdhictop(ji) + zdhicbot(ji) ) * hakspl * ( 1.0 - zfrl_old(ji) ) &
775             &  / MAX( epsi13 , h_ice_1d(ji) + h_snow_1d(ji) * rhosn/rhoic ) 
776          zfrld_1d(ji) =  frld_1d(ji) + MAX( zzero , zdfrl )
777          !--Limitation of sea-ice fraction <= 1
778          zfrld_1d(ji) = ziqf * MIN( 0.99 * zone , zfrld_1d(ji) ) + ( 1 - ziqf )
779          !---Update surface and internal temperature and snow/ice thicknesses.
780          sist_1d(ji)   = sist_1d(ji)   + ( 1.0 - ziqf ) * ( tfu_1d(ji) - sist_1d(ji)   )
781          tbif_1d(ji,1) = tbif_1d(ji,1) + ( 1.0 - ziqf ) * ( tfu_1d(ji) - tbif_1d(ji,1) )
782          tbif_1d(ji,2) = tbif_1d(ji,2) + ( 1.0 - ziqf ) * ( tfu_1d(ji) - tbif_1d(ji,2) )
783          tbif_1d(ji,3) = tbif_1d(ji,3) + ( 1.0 - ziqf ) * ( tfu_1d(ji) - tbif_1d(ji,3) )
784          !--variation of ice volume and ice mass
785          dvlbq_1d(ji)   = zihic * ( zfrl_old(ji) - frld_1d(ji) ) * h_ice_1d(ji)
786          rdmicif_1d(ji) = rdmicif_1d(ji) + dvlbq_1d(ji) * rhoic
787          !--variation of snow volume and snow mass
788          zdvsnvol    = zihsn * ( zfrl_old(ji) - frld_1d(ji) ) * h_snow_1d(ji)
789          rdmsnif_1d(ji) = rdmsnif_1d(ji) + zdvsnvol * rhosn
790          h_snow_1d(ji)  = ziqf * h_snow_1d(ji)
791
792          zdrfrl1 = ziqf * ( 1.0 -  frld_1d(ji) ) / MAX( epsi20 , 1.0 - zfrld_1d(ji) )
793          zdrfrl2 = ziqf * ( 1.0 - zfrl_old(ji) ) / MAX( epsi20 , 1.0 - zfrld_1d(ji) )
794
795          h_snow_1d (ji) = zdrfrl1 * h_snow_1d(ji)
796          h_ice_1d  (ji) = zdrfrl1 * h_ice_1d(ji)
797          qstbif_1d(ji) = zdrfrl2 * qstbif_1d(ji)
798          frld_1d(ji)    = zfrld_1d(ji)
799
800       END DO
801       
802    END SUBROUTINE lim_thd_zdf
803#else
804   !!======================================================================
805   !!                       ***  MODULE limthd_zdf   ***
806   !!                              no sea ice model 
807   !!======================================================================
808CONTAINS
809   SUBROUTINE lim_thd_zdf          ! Empty routine
810   END SUBROUTINE lim_thd_zdf
811#endif
812 END MODULE limthd_zdf
Note: See TracBrowser for help on using the repository browser.