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

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

source: branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/LIM_SRC_2/limthd_zdf_2.F90 @ 7470

Last change on this file since 7470 was 7470, checked in by jcastill, 7 years ago

Remove svn keys

File size: 47.0 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 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       !--------------------------------------------------------------------------------
276       !  4. Computation of the surface temperature : determined by considering the
277       !     budget of a thin layer of thick. zhsu at the top surface (H. Grenier, 1995)
278       !     and based on a surface energy balance :
279       !     hsu * rcp * dT/dt = Fsr + Fnsr(T) + Fcs(T),
280       !     where - Fsr is the net absorbed solar radiation,
281       !           - Fnsr is the total non solar radiation (incoming and outgoing long-wave,
282       !             sensible and latent heat fluxes)
283       !           - Fcs the conductive heat flux at the top of surface
284       !------------------------------------------------------------------------------
285
286       !     4.1. Computation of intermediate values
287       !---------------------------------------------
288       DO ji = kideb, kiut
289          zrcpdt(ji) = ( rcpsn * MIN( h_snow_1d(ji) , zhsu )    &
290             &       + rcpic * MAX( zhsu - h_snow_1d(ji) , zzero ) ) / rdt_ice
291          zts_old(ji) =  sist_1d(ji)
292       END DO
293
294       !     4.2. Computation of surface temperature by expanding the eq. of energy balance
295       !          with Ts = Tp + DT. One obtain , F(Tp) + DT * DF(Tp) = 0
296       !          where  - F(Tp) = Fsr + Fnsr(Tp) + Fcs(Tp)
297       !                 - DF(Tp)= (dFnsr(Tp)/dT) + (dFcs(Tp)/dT) - hsu*rcp/dt
298       !---------------------------------------------------------------------------------
299
300       DO ji = kideb, kiut
301          !---computation of the derivative of energy balance function
302          zdfts    =  zksndh(ji)   & ! contribution of the conductive heat flux
303             &      + zrcpdt(ji)   & ! contribution of hsu * rcp / dt
304             &      - dqns_ice_1d (ji)     ! contribution of the total non solar radiation
305          !---computation of the energy balance function
306          zfts    = - z1mi0 (ji) * qsr_ice_1d(ji)   & ! net absorbed solar radiation
307             &      - qns_ice_1d(ji)                & ! total non solar radiation
308             &      - zfcsu (ji)                      ! conductive heat flux from the surface
309          !---computation of surface temperature increment 
310          zdts    = -zfts / zdfts
311          !---computation of the new surface temperature
312          sist_1d(ji) = sist_1d(ji) + zdts
313       END DO
314
315       !----------------------------------------------------------------------------
316       !  5. Boundary condition at the top surface
317       !--    IF Tsb < Tmelt, Fnet = Fcs (the net heat flux equal the conductive heat flux)
318       !      Otherwise Tsb = Tmelt and Qnet(Tmelt) > 0
319       !      Fnet(Tmelt) is therefore the net surface flux needed for melting
320       !----------------------------------------------------------------------------
321       
322       
323       !     5.1.  Limitation of surface temperature and update total non solar fluxes,
324       !          latent heat flux and conductive flux at the top surface
325       !---------------------------------------------------------------------- 
326                     
327       IF ( .NOT. ln_cpl ) THEN   ! duplicate the loop for performances issues
328          DO ji = kideb, kiut
329             sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) )
330             qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) )
331             qla_ice_1d(ji) = qla_ice_1d(ji) + dqla_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) )
332             zfcsu(ji)  = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) )
333          END DO
334       ELSE
335          DO ji = kideb, kiut
336             sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) )
337             qla_ice_1d(ji) = -9999.   ! default definition, not used as parsub = 0. in this case
338             zfcsu(ji)  = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) )
339          END DO
340       ENDIF
341
342       !     5.2. Calculate available heat for surface ablation.
343       !---------------------------------------------------------------------
344
345       DO ji = kideb, kiut
346          zfnet(ji) = qns_ice_1d(ji) + z1mi0(ji) * qsr_ice_1d(ji) + zfcsu(ji)         
347          zfnet(ji) = MAX( zzero , zfnet(ji) )
348          zfnet(ji) = zfnet(ji) * MAX( zzero , SIGN( zone , sist_1d(ji) - ztsmlt(ji) ) )
349       END DO
350
351       !-------------------------------------------------------------------------
352       !  6. Calculate changes in internal temperature due to vertical diffusion   
353       !     processes. The evolution of this temperature is governed by the one-
354       !     dimensionnal heat-diffusion equation.
355       !     Given the temperature tbif(1/2/3), at time m we solve a set
356       !     of finite difference equations to obtain new tempe. Each tempe is coupled
357       !     to the temp. immediatly above and below by heat conduction terms. Thus
358       !     we have a set of equations of the form A * T = B, where A is a tridiagonal
359       !     matrix, T a vector whose components are the unknown new temp.
360       !-------------------------------------------------------------------------
361       
362       !--parameter for the numerical methode use to solve the heat-diffusion equation
363       !- implicit, explicit or Crank-Nicholson
364       zumsb = 1.0 - sbeta 
365       DO ji = kideb, kiut
366          zidsn(ji)   = MAX ( zzero, SIGN( zone, hsndif - h_snow_1d(ji) ) ) 
367          z1midsn(ji) = 1.0 - zidsn(ji)
368          zihic       = MAX ( zzero, SIGN( zone, hicdif - h_ice_1d(ji) ) ) 
369          zidsnic(ji) = zidsn(ji) *  zihic 
370          zfcsudt(ji) = zfcsu(ji) * rdt_ice 
371       END DO
372   
373       DO ji = kideb, kiut
374
375          !     6.1 Calculate intermediate variables.
376          !----------------------------------------
377
378          !--conductivity at the snow surface
379          zkhsn = 2.0 * zksn(ji) * rdt_ice / rcpsn
380          !--conductivity at the ice surface
381          zkhic = 4.0 * zkic(ji) * rdt_ice / MAX( h_ice_1d(ji) * h_ice_1d(ji) * rcpic , epsi20 )
382          !--conductivity at the snow/ice interface
383          zkint = 4.0 * zksn(ji) * zkic(ji)  &
384             &        / ( zksn(ji) * h_ice_1d(ji) + 2.0 * zkic(ji) * h_snow_1d(ji) * z1midsn(ji)) 
385          zkhsnint = zkint * rdt_ice / rcpsn
386          zkhicint = zkint * 2.0 * rdt_ice / MAX( h_ice_1d(ji) * rcpic , epsi20 )
387         
388          !     6.2. Fulfill the linear system matrix.
389          !-----------------------------------------
390!$$$          zplediag(1) = 1 + sbeta * z1midsn(ji) * ( zkhsn + zkhsnint )       
391          zplediag(1) =   zidsn(ji) + z1midsn(ji) * h_snow_1d(ji)   &
392             &          + sbeta * z1midsn(ji) * zkhsnint 
393          zplediag(2) = 1 + sbeta * ( z1midsn(ji) * zkhicint + zkhic ) 
394          zplediag(3) = 1 + 3.0 * sbeta * zkhic   
395
396          zsubdiag(1) =  0.e0             
397          zsubdiag(2) = -1.e0 * z1midsn(ji) * sbeta * zkhicint
398          zsubdiag(3) = -1.e0 * sbeta * zkhic 
399
400          zsupdiag(1) = -1.e0 * z1midsn(ji) * sbeta * zkhsnint 
401          zsupdiag(2) = zsubdiag(3)
402          zsupdiag(3) =  0.e0
403         
404          !     6.3. Fulfill the idependent term vector.
405          !-------------------------------------------
406         
407!$$$          zsmbr(1) = zidsn(ji) * sist_1d(ji) + z1midsn(ji) *   &
408!$$$             &         ( tbif_1d(ji,1) + zkhsn * sist_1d(ji)
409!$$$             &         - zumsb * ( zkhsn * tbif_1d(ji,1)
410!$$$             &                   + zkhsnint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) ) )
411          zsmbr(1) = zidsn(ji) * sist_1d(ji) + z1midsn(ji) *    &
412             &       ( h_snow_1d(ji) * tbif_1d(ji,1) - ( zfcsudt(ji) / rcpsn )  &
413             &       - zumsb * zkhsnint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) )
414
415          zsmbr(2) =  tbif_1d(ji,2)  &
416             &      - zidsn(ji) * ( 1.0 - zidsnic(ji) ) &
417             &        * ( zfcsudt(ji) / MAX( h_ice_1d(ji) * rcpic , epsi20 ) ) &
418             &      + zumsb * ( zkhicint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) &
419             &                   - zkhic * ( tbif_1d(ji,2) - tbif_1d(ji,3) )  )
420
421          zsmbr(3) =  tbif_1d(ji,3)  &
422             &      + zkhic * ( 2.0 * tfu_1d(ji) &
423             &                + zumsb * ( tbif_1d(ji,2) - 3.0 * tbif_1d(ji,3) ) )
424         
425          !     6.4. Solve the system (Gauss elimination method).
426          !----------------------------------------------------
427         
428          zpiv1 = zsubdiag(2) / zplediag(1) 
429          zb2   = zplediag(2) - zpiv1 * zsupdiag(1)
430          zd2   = zsmbr(2) - zpiv1 * zsmbr(1)
431
432          zpiv2 = zsubdiag(3) / zb2
433          zb3   = zplediag(3) - zpiv2 * zsupdiag(2)
434          zd3   = zsmbr(3) - zpiv2 * zd2
435
436          tbif_1d(ji,3) = zd3 / zb3
437          tbif_1d(ji,2) = ( zd2 - zsupdiag(2) * tbif_1d(ji,3) ) / zb2
438          tbif_1d(ji,1) = ( zsmbr(1) - zsupdiag(1) * tbif_1d(ji,2) ) / zplediag(1)           
439
440          !--- taking into account the particular case of  zidsnic(ji) = 1
441          ztint =  (  zkic(ji) * h_snow_1d(ji) * tfu_1d (ji)    &
442             &      + zksn(ji) * h_ice_1d(ji) * sist_1d(ji) )   &
443             &   / ( zkic(ji) * h_snow_1d(ji) + zksn(ji) * h_ice_1d(ji) ) 
444
445          tbif_1d(ji,1) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,1)   &
446             &                + zidsnic(ji)   * ( ztint + sist_1d(ji) ) / 2.0
447          tbif_1d(ji,2) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,2)   &
448             &                + zidsnic(ji)   * ( 3.0 * ztint + tfu_1d(ji) ) / 4.0
449          tbif_1d(ji,3) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,3)   &
450             &                + zidsnic(ji)   * ( ztint + 3.0 * tfu_1d(ji) ) / 4.0     
451       END DO
452 
453       !----------------------------------------------------------------------
454       !  9. Take into account surface ablation and bottom accretion-ablation.|
455       !----------------------------------------------------------------------
456       
457       !---Snow accumulation in one thermodynamic time step
458       zsprecip(kideb:kiut) = sprecip_1d(kideb:kiut) * rdt_ice / rhosn
459
460
461       DO ji = kideb, kiut
462         
463          !      9.1. Surface ablation and update of snow thickness and qstbif_1d
464          !--------------------------------------------------------------------
465         
466          !--------------------------------------------------------------------------
467          !--      Melting snow processes :
468          !--      Melt at the upper surface is computed from the difference between
469          !--      the net heat flux (including the conductive heat flux) at the upper
470          !--      surface and the pre-existing energy due to surface melting
471          !------------------------------------------------------------------------------
472         
473          !-- store the snow thickness
474          zhsnw_old(ji) =  h_snow_1d(ji)
475          !--computation of the energy needed to melt snow
476          zqsnw_mlt  = zfnet(ji) * rdt_ice - zqcmlts(ji)
477          !--change in snow thickness due to melt
478          zdhsmlt = - zqsnw_mlt / xlsn
479         
480          !-- compute new snow thickness, taking into account the part of snow accumulation
481          !   (as snow precipitation) and the part of snow lost due to melt
482          zhsn =  h_snow_1d(ji) + zsprecip(ji) + zdhsmlt
483          h_snow_1d(ji) = MAX( zzero , zhsn )
484          !-- compute the volume of snow lost after surface melting and the associated mass
485          dvsbq_1d(ji) =  ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnw_old(ji) - zsprecip(ji) )
486          dvsbq_1d(ji) =  MIN( zzero , dvsbq_1d(ji) )
487          ztmp = rhosn * dvsbq_1d(ji)
488          rdm_snw_1d(ji) =  ztmp
489          !--heat content of the water provided to the ocean (referenced to rt0)
490          rdq_snw_1d(ji) =  cpic * ztmp * ( rt0_snow - rt0 )
491          !-- If the snow is completely melted the remaining heat is used to melt ice
492          zqsn_mlt_rem  = MAX( zzero , -zhsn ) * xlsn
493          zqice_top_mlt = zqsn_mlt_rem 
494          zqstbif_old   = qstbif_1d(ji)
495
496          !--------------------------------------------------------------------------
497          !--      Melting ice processes at the top surface :
498          !--      The energy used to melt ice, zqice_top_mlt, is taken from the energy
499          !--      stored in brine pockets qstbif_1d and the remaining energy coming
500          !--      from the melting snow process zqsn_mlt_rem.
501          !--      If qstbif_1d > zqsn_mlt_rem then, one uses only a zqsn_mlt_rem part
502          !--      of qstbif_1d to melt ice,
503          !--         zqice_top_mlt = zqice_top_mlt + zqsn_mlt_rem
504          !--         qstbif_1d = qstbif_1d - zqsn_mlt_rem
505          !--      Otherwise one uses all qstbif_1d to melt ice
506          !--         zqice_top_mlt = zqice_top_mlt + qstbif_1d
507          !--         qstbif_1d = 0
508          !------------------------------------------------------
509         
510          ziqf =  MAX ( zzero , SIGN( zone , qstbif_1d(ji) - zqsn_mlt_rem  ) )
511          zqice_top_mlt =         ziqf   * ( zqice_top_mlt + zqsn_mlt_rem )   &
512             &          + ( 1.0 - ziqf ) * ( zqice_top_mlt + qstbif_1d(ji)  )
513
514          qstbif_1d(ji) =         ziqf   * ( qstbif_1d(ji) - zqsn_mlt_rem )   &
515             &          + ( 1.0 - ziqf ) * ( qstbif_1d(ji) - qstbif_1d(ji)  )
516
517          !--    The contribution of the energy stored in brine pockets qstbif_1d to melt
518          !--    ice is taking into account only when qstbif_1d is less than zqmax.
519          !--    Otherwise, only the remaining energy coming from the melting snow
520          !--    process is used
521          zihq =  MAX ( zzero , SIGN( zone , qstbif_1d(ji) - zqmax(ji) ) )
522
523          zqice_top_mlt =         zihq   * zqice_top_mlt   &
524             &          + ( 1.0 - zihq ) * zqsn_mlt_rem
525
526          qstbif_1d(ji) =         zihq   * qstbif_1d(ji)   &
527             &          + ( 1.0 - zihq ) * zqstbif_old
528
529          !--change in ice thickness due to melt at the top surface
530          zdhictop(ji) = -zqice_top_mlt / xlic
531          !--compute the volume formed after surface melting
532          dvsbq_1d(ji) =  zdhictop(ji) * ( 1.0 - frld_1d(ji) )
533
534          !-------------------------------------------------------------------------
535          !--      A small variation at the surface also occurs because of sublimation
536          !--      associated with the latent flux. If qla_ice_1d is negative, snow condensates at
537          !        the surface. Otherwise, snow evaporates
538          !-----------------------------------------------------------------------
539          !----change in snow and ice thicknesses due to sublimation or evaporation
540          zdhssub  = parsub * ( qla_ice_1d(ji) / ( rhosn * xsn ) ) * rdt_ice 
541          zhsn     = h_snow_1d(ji) - zdhssub
542          zdhisub  = MAX( zzero , -zhsn ) * rhosn/rhoic
543          zdhictop(ji) =  zdhictop(ji) - zdhisub
544          h_snow_1d(ji)  =  MAX( zzero , zhsn )
545          !-------------------------------------------------
546          !--  Update Internal temperature and qstbif_1d.
547          !-------------------------------------------
548          zihsn  =  MAX( zzero , SIGN( zone, -h_snow_1d(ji) ) )
549          tbif_1d(ji,1) = ( 1.0 - zihsn ) * tbif_1d(ji,1) + zihsn   * tfu_1d(ji)
550          !--change in snow internal temperature if snow has increased
551          zihnf = MAX( zzero , SIGN( zone , h_snow_1d(ji) - zhsnw_old(ji) ) )
552          zdhsn = 1.0 - zhsnw_old(ji) / MAX( h_snow_1d(ji) , epsi20 )
553          zdtsn = zdhsn * ( sist_1d(ji) - tbif_1d(ji,1) )
554          tbif_1d(ji,1) = tbif_1d(ji,1) + z1midsn(ji) * zihnf * zdtsn
555          !--energy created due to ice melting in the first ice layer
556          zqnes  = ( rt0_ice - tbif_1d(ji,2) ) * rcpic * ( h_ice_1d(ji) / 2. )
557          !--change in first ice layer internal temperature
558          ziqr  = MAX( zzero , SIGN( zone , qstbif_1d(ji) - zqnes ) )
559          zdtic = qstbif_1d(ji) / ( rcpic * ( h_ice_1d(ji) / 2. ) )
560          tbif_1d(ji,2) =  ziqr * rt0_ice + ( 1 - ziqr ) * ( tbif_1d(ji,2) + zdtic )
561          !--update qstbif_1d
562          qstbif_1d(ji) = ziqr * ( qstbif_1d(ji) - zqnes ) * swiqst
563
564
565          !--      9.2. Calculate bottom accretion-ablation and update qstbif_1d.
566          !             Growth and melting at bottom ice surface are governed by 
567          !                 -xlic * Dh = (Fcb - Fbot ) * Dt
568          !             where Fbot is the net downward heat flux from ice to the ocean
569          !            and Fcb is the conductive heat flux at the bottom surface
570          !---------------------------------------------------------------------------
571          ztbot = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,3) + zidsnic(ji) * sist_1d(ji)
572          !---computes conductive heat flux at bottom surface
573          zfcbot =  4.0 * zkic(ji) * ( tfu_1d(ji) - ztbot )   &
574             &   / ( h_ice_1d(ji) + zidsnic(ji) * ( 3. * h_ice_1d(ji) &
575             &   + 4.0 * zkic(ji)/zksn(ji) * h_snow_1d(ji) ) )
576          !---computation of net energy needed for bottom melting/growing
577          zqice_bot = ( zfcbot - ( fbif_1d(ji) + qlbbq_1d(ji) ) ) * rdt_ice
578          zqstbif_bot = qstbif_1d(ji)
579          !---switch to know if bottom surface melts ( = 1 ) or grows ( = 0 )occurs
580          zibmlt = MAX( zzero , SIGN( zone , -zqice_bot ) )
581          !--particular case of melting (in the same way as the top surface)
582          zqice_bot_mlt = zqice_bot 
583          zqstbif_old = zqstbif_bot
584
585          ziqf =  MAX ( zzero , SIGN( zone , qstbif_1d(ji) + zqice_bot_mlt  ) )
586          zqice_bot_mlt =         ziqf   * ( zqice_bot_mlt + zqice_bot_mlt ) &
587             &          + ( 1.0 - ziqf ) * ( zqice_bot_mlt + qstbif_1d(ji)  )         
588          qstbif_1d(ji)   =         ziqf   * ( qstbif_1d(ji) + zqice_bot_mlt ) &
589             &          + ( 1.0 - ziqf ) * ( qstbif_1d(ji) - qstbif_1d(ji)  )
590          !--    The contribution of the energy stored in brine pockets qstbif_1d to melt
591          !--    ice is taking into account only when qstbif_1d is less than zqmax.
592          zihq =  MAX ( zzero , SIGN( zone , qstbif_1d(ji) - zqmax(ji) ) )
593          zqice_bot_mlt =         zihq   * zqice_bot_mlt   &
594             &          + ( 1.0 - zihq ) * zqice_bot
595          qstbif_1d(ji)   =         zihq   * qstbif_1d(ji)   &
596             &             + ( 1.0 - zihq ) * zqstbif_old
597
598          !---treatment of the case of melting/growing
599          zqice_bot   =         zibmlt   * ( zqice_bot_mlt - zqcmltb(ji) )   &
600             &        + ( 1.0 - zibmlt ) * ( zqice_bot - zqcmltb(ji)  )
601          qstbif_1d(ji) =         zibmlt   * qstbif_1d(ji)   &
602             &           + ( 1.0 - zibmlt ) * zqstbif_bot
603
604          !--computes change in ice thickness due to melt or growth
605          zdhicbot(ji) = zqice_bot / xlic
606          !--limitation of bottom melting if so : hmelt maximum melting at bottom
607          zdhicmlt  = MAX( hmelt , zdhicbot(ji) ) 
608          !-- output part due to bottom melting only
609          IF( zdhicmlt < 0.e0 ) rdvomif_1d(ji) = ( 1.0 - frld_1d(ji) ) * zdhicmlt
610          !--energy after bottom melting/growing
611          zqsup(ji) = ( 1.0 - frld_1d(ji) ) * xlic * ( zdhicmlt - zdhicbot(ji) )
612          !-- compute the new thickness and the newly formed volume after bottom melting/growing
613          zdhicbot(ji)  = zdhicmlt
614          dvbbq_1d(ji) = ( 1.0 - frld_1d(ji) ) * zdhicbot(ji)
615
616
617          !        9.3.  Updating ice thickness after top surface ablation
618          !              and bottom surface accretion/ablation
619          !---------------------------------------------------------------
620          zhicnew  = h_ice_1d(ji) + zdhictop(ji) + zdhicbot(ji)
621
622          !
623          !        9.4. Case of total ablation (ice is gone but snow may be left)
624          !-------------------------------------------------------------------
625          zhsn  = h_snow_1d(ji)
626          zihgnew = 1.0 - MAX( zzero , SIGN( zone , -zhicnew ) )
627          zihsn   = MAX( zzero , SIGN( zone , -zhsn ) )
628          !---convert
629          zdhicm  = ( 1.0 - zihgnew ) * ( zhicnew - qstbif_1d(ji) / xlic )
630          zdhsnm  = ( 1.0 - zihsn ) * zdhicm * rhoic / rhosn
631          !---updating new ice thickness and computing the newly formed ice mass
632          zhicnew   =  zihgnew * zhicnew
633          ztmp    =  ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) ) * rhoic
634          rdm_ice_1d(ji) =  rdm_ice_1d(ji) + ztmp
635          !---heat content of the water provided to the ocean (referenced to rt0)
636          !   use of rt0_ice is OK for melting ice; in the case of freezing, tfu_1d should be used.
637          !   This is done in 9.5 section (see below)
638          rdq_ice_1d(ji) =  cpic * ztmp * ( rt0_ice - rt0 )
639          !---updating new snow thickness and computing the newly formed snow mass
640          zhsnfi   = zhsn + zdhsnm
641          h_snow_1d(ji) = MAX( zzero , zhsnfi )
642          ztmp = ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsn ) * rhosn
643          rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp
644          !---updating the heat content of the water provided to the ocean (referenced to rt0)
645          rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( rt0_snow - rt0 )
646          !--remaining energy in case of total ablation
647          zqocea(ji) = - ( zihsn * xlic * zdhicm + xlsn * ( zhsnfi - h_snow_1d(ji) ) ) * ( 1.0 - frld_1d(ji) )
648          qstbif_1d(ji) = zihgnew * qstbif_1d(ji)
649
650          !
651          !        9.5. Update internal temperature and ice thickness.
652          !-------------------------------------------------------
653          !
654          sist_1d(ji) = zihgnew * sist_1d(ji) + ( 1.0 - zihgnew ) * tfu_1d(ji)
655          zidhb  = MAX( zzero , SIGN( zone , - zdhicbot(ji) ) )
656          zc1    = - zhicnew * 0.5
657          zpc1   = MIN( 0.5 * zone , - h_ice_1d(ji) * 0.5 - zdhictop(ji) )
658          zc2    = - zhicnew
659          zpc2   =  zidhb * zc2 + ( 1.0 - zidhb ) * ( - h_ice_1d(ji) - zdhictop(ji) )
660          zp1    =  MAX( zpc1 , zc1 )
661          zp2    =  MAX( zpc2 , zc1 )
662          zep(ji) =  tbif_1d(ji,2)
663          ztb2  = 2.0 * (         - zp1   * tbif_1d(ji,2)  &
664             &  + ( zp1 - zp2 ) * tbif_1d(ji,3)  &
665             &  + ( zp2 - zc1 ) * tfu_1d(ji) ) / MAX( zhicnew , epsi20 ) 
666          tbif_1d(ji,2) = zihgnew * ztb2 + ( 1.0 - zihgnew ) * tfu_1d(ji)
667          !---
668          zp1  =  MIN( zpc1 , zc1 )
669          zp2  =  MIN( zpc2 , zc1 )
670          zp1  =  MAX( zc2  , zp1 )
671          ztb3 =  2.0 * (   ( 1.0 - zidhb ) * (  ( zc1 - zp2 ) * tbif_1d(ji,3)  &
672             &                                 + ( zp2 - zc2 ) * tfu_1d(ji) )   &
673             &               +      zidhb   * (  ( zc1 - zp1 ) * zep(ji)      &
674             &                                 + ( zp1 - zc2 ) * tbif_1d(ji,3))  ) / MAX( zhicnew , epsi20 )
675          tbif_1d(ji,3) =  zihgnew * ztb3 + ( 1.0 - zihgnew ) * tfu_1d(ji)
676          h_ice_1d(ji)  =  zhicnew
677          ! update the ice heat content given to the ocean in freezing case
678          ! (part due to difference between rt0_ice and tfu_1d)
679          ztmp = ( 1. - zidhb ) * rhoic * dvbbq_1d(ji)
680          rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * ( tfu_1d(ji) - rt0_ice )
681       END DO
682
683
684       !----------------------------------------------------------------------------
685       !  10. Surface accretion.
686       !      The change of ice thickness after snow/ice formation is such that
687       !      the interface between snow and ice is located at the same height
688       !      as the ocean surface. It is given by (Fichefet and Morales Maqueda 1999)
689       !          D(h_ice) = (- D(hsn)/alph) =  [rhosn*hsn - (rau0 - rhoic)*hic]
690       !                                     / [alph*rhosn+rau0 - rhoic]
691       !----------------------------------------------------------------------------
692       !
693       DO ji = kideb , kiut
694
695          !--  Computation of the change of ice thickness after snow-ice formation
696          zdrmh =  ( rhosn * h_snow_1d(ji) + ( rhoic - rau0 ) * h_ice_1d(ji) )  &
697             &  / ( alphs * rhosn + rau0 - rhoic )
698          zdrmh = MAX( zzero , zdrmh )
699
700          !--New ice and snow thicknesses Fichefet and Morales Maqueda (1999)
701          zhicnew  = MAX( h_ice_1d(ji) , h_ice_1d(ji) + zdrmh )
702          zhsnnew  = MIN( h_snow_1d(ji) , h_snow_1d(ji) - alphs * zdrmh )
703          !---Compute new ice temperatures. snow temperature remains unchanged
704          !   Lepparanta (1983):
705          zihic = 1.0 - MAX( zzero , SIGN( zone , -zhicnew ) )
706          zquot  = ( 1.0 - zihic ) &
707             &   +         zihic * MIN( zone , h_ice_1d(ji) / MAX( zhicnew , epsi20 ) ) 
708          ztneq  =         alphs * cnscg * tbif_1d(ji,1)    &
709             &   + ( 1.0 - alphs * ( rhosn/rhoic ) ) * tfu_1d(ji)
710          zep(ji) = tbif_1d(ji,2)
711          tbif_1d(ji,2) = ztneq - zquot * zquot * ( ztneq - tbif_1d(ji,2) )
712          tbif_1d(ji,3) = 2.0 * ztneq &
713             &        + zquot * ( tbif_1d(ji,3) + zep(ji) - 2.0 * ztneq ) - tbif_1d(ji,2)
714
715          !---  Lepparanta (1983) (latent heat released during white ice formation
716          !     goes to the ocean -for lateral ablation-)
717          qldif_1d(ji)  = qldif_1d(ji) + zdrmh * ( 1.0 - alphs * ( rhosn/rhoic ) ) * xlic * ( 1.0 - frld_1d(ji) )
718          !--   Changes in ice volume and ice mass Lepparanta (1983):
719          dvnbq_1d(ji) = ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) )
720          dmgwi_1d(ji) = dmgwi_1d(ji) + ( 1.0 -frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnnew ) * rhosn
721          !---  volume change of ice and snow (used for ocean-ice freshwater flux computation)
722          ztmp = ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d (ji) ) * rhoic
723          rdm_ice_1d(ji) = rdm_ice_1d(ji) + ztmp
724          rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * ( tfu_1d(ji) - rt0 )
725          !!gm BUG ??   snow ==>  only needed for nn_ice_embd == 0  (standard levitating sea-ice)
726          ztmp = ( 1.0 - frld_1d(ji) )   * ( zhsnnew - h_snow_1d(ji) ) * rhosn
727          rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp
728          rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( rt0_snow - rt0 )
729
730          !---  Actualize new snow and ice thickness.
731          h_snow_1d(ji)  = zhsnnew
732          h_ice_1d (ji)  = zhicnew
733
734       END DO
735
736       !----------------------------------------------------
737       !  11. Lateral ablation (Changes in sea/ice fraction)
738       !----------------------------------------------------
739       DO ji = kideb , kiut
740          zfrl_old(ji)  = frld_1d(ji)
741          zihic   = 1.0 - MAX( zzero , SIGN( zone , -h_ice_1d(ji) ) )
742          zihsn   = 1.0 - MAX( zzero , SIGN( zone , -h_snow_1d(ji) ) )
743          !--In the case of total ablation (all the ice ice has melted) frld = 1
744          frld_1d(ji)  = ( 1.0 - zihic ) + zihic * zfrl_old(ji)
745          !--Part of solar radiation absorbing inside the ice and going
746          !--through the ocean
747          fscbq_1d(ji) = ( 1.0 - zfrl_old(ji) ) * ( 1.0 - thcm_1d(ji) ) * fstbif_1d(ji)
748          !--Total remaining energy after bottom melting/growing
749          qfvbq_1d(ji) = zqsup(ji) + ( 1.0 - zihic ) * zqocea(ji)
750          !--Updating of total heat from the ocean
751          qldif_1d(ji)  = qldif_1d(ji) + qfvbq_1d(ji) + ( 1.0 - zihic ) * fscbq_1d(ji) * rdt_ice
752          !--Computation of total heat inside the snow/ice system
753          zqice  = h_snow_1d(ji) * xlsn + h_ice_1d(ji) * xlic
754          zqicetot  = ( 1.0 - frld_1d(ji) ) * zqice
755          !--The concentration of ice is reduced (frld increases) if the heat
756          !--exchange between ice and ocean is positive
757          ziqf = MAX( zzero , SIGN( zone ,  zqicetot - qldif_1d(ji) ) )
758          zdfrl = qldif_1d(ji) / MAX( epsi20 , zqice ) 
759          frld_1d(ji)  = ( 1.0 - ziqf )    &
760             &       +         ziqf * ( frld_1d(ji) + MAX( zzero , zdfrl ) ) 
761          fltbif_1d(ji) = ( ( 1.0 - zfrl_old(ji) ) * qstbif_1d(ji) - zqicetot  ) / rdt_ice
762          !--  Opening of leads: Hakkinen & Mellor, 1992.
763          zdfrl = - ( zdhictop(ji) + zdhicbot(ji) ) * hakspl * ( 1.0 - zfrl_old(ji) ) &
764             &  / MAX( epsi13 , h_ice_1d(ji) + h_snow_1d(ji) * rhosn/rhoic ) 
765          zfrld_1d(ji) =  frld_1d(ji) + MAX( zzero , zdfrl )
766          !--Limitation of sea-ice fraction <= 1
767          zfrld_1d(ji) = ziqf * MIN( 0.99 * zone , zfrld_1d(ji) ) + ( 1 - ziqf )
768          !---Update surface and internal temperature and snow/ice thicknesses.
769          sist_1d(ji)   = sist_1d(ji)   + ( 1.0 - ziqf ) * ( tfu_1d(ji) - sist_1d(ji)   )
770          tbif_1d(ji,1) = tbif_1d(ji,1) + ( 1.0 - ziqf ) * ( tfu_1d(ji) - tbif_1d(ji,1) )
771          tbif_1d(ji,2) = tbif_1d(ji,2) + ( 1.0 - ziqf ) * ( tfu_1d(ji) - tbif_1d(ji,2) )
772          tbif_1d(ji,3) = tbif_1d(ji,3) + ( 1.0 - ziqf ) * ( tfu_1d(ji) - tbif_1d(ji,3) )
773          !--variation of ice volume and ice mass
774          dvlbq_1d(ji)   = zihic * ( zfrl_old(ji) - frld_1d(ji) ) * h_ice_1d(ji)
775          ztmp = dvlbq_1d(ji) * rhoic
776          rdm_ice_1d(ji) = rdm_ice_1d(ji) + ztmp
777!!gm
778!!gm   This should be split in two parts:
779!!gm         1-  heat required to bring sea-ice to tfu  : this part should be added to the heat flux taken from the ocean
780!!gm                 cpic * ztmp * 0.5 * ( tbif_1d(ji,2) + tbif_1d(ji,3) - 2.* rt0_ice )
781!!gm         2-  heat content of lateral ablation referenced to rt0 : this part only put in rdq_ice_1d
782!!gm                 cpic * ztmp * ( rt0_ice - rt0 )
783!!gm   Currently we put all the heat in rdq_ice_1d
784          rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * 0.5 * ( tbif_1d(ji,2) + tbif_1d(ji,3) - 2.* rt0 )
785          !
786          !--variation of snow volume and snow mass
787          zdvsnvol = zihsn * ( zfrl_old(ji) - frld_1d(ji) ) * h_snow_1d(ji)
788          ztmp     = zdvsnvol * rhosn
789          rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp
790!!gm
791!!gm   This should be split in two parts:
792!!gm         1-  heat required to bring snow to tfu  : this part should be added to the heat flux taken from the ocean
793!!gm                 cpic * ztmp * ( tbif_1d(ji,1) - rt0_snow )
794!!gm         2-  heat content of lateral ablation referenced to rt0 : this part only put in rdq_snw_1d
795!!gm                 cpic * ztmp * ( rt0_snow - rt0 )
796!!gm   Currently we put all the heat in rdq_snw_1d
797          rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( tbif_1d(ji,1) - rt0 )
798
799          h_snow_1d(ji)  = ziqf * h_snow_1d(ji)
800
801          zdrfrl1 = ziqf * ( 1.0 -  frld_1d(ji) ) / MAX( epsi20 , 1.0 - zfrld_1d(ji) )
802          zdrfrl2 = ziqf * ( 1.0 - zfrl_old(ji) ) / MAX( epsi20 , 1.0 - zfrld_1d(ji) )
803
804          h_snow_1d (ji) = zdrfrl1 * h_snow_1d(ji)
805          h_ice_1d  (ji) = zdrfrl1 * h_ice_1d(ji)
806          qstbif_1d(ji) = zdrfrl2 * qstbif_1d(ji)
807          frld_1d(ji)    = zfrld_1d(ji)
808          !
809       END DO
810       !
811       CALL wrk_dealloc( jpij, ztsmlt, ztbif  , zksn    , zkic    , zksndh , zfcsu  , zfcsudt , zi0      , z1mi0   , zqmax    )
812       CALL wrk_dealloc( jpij, zrcpdt, zts_old, zidsn   , z1midsn , zidsnic, zfnet  , zsprecip, zhsnw_old, zdhictop, zdhicbot )
813       CALL wrk_dealloc( jpij, zqsup , zqocea , zfrl_old, zfrld_1d, zep    , zqcmlts, zqcmltb                                 )
814       !
815    END SUBROUTINE lim_thd_zdf_2
816
817#else
818   !!----------------------------------------------------------------------
819   !!   Default Option                                     NO sea-ice model
820   !!----------------------------------------------------------------------
821CONTAINS
822   SUBROUTINE lim_thd_zdf_2          ! Empty routine
823   END SUBROUTINE lim_thd_zdf_2
824#endif
825
826   !!======================================================================
827END MODULE limthd_zdf_2
Note: See TracBrowser for help on using the repository browser.