source: branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_2/limthd_zdf_2.F90 @ 5206

Last change on this file since 5206 was 5206, checked in by cetlod, 7 years ago

First improvment of diurnal cycle for passive tracers, see ticket #1484

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