source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_2/limthd_zdf_2.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

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