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

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

source: branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 3963

Last change on this file since 3963 was 3963, checked in by clem, 11 years ago

bugs correction + creation of glob_max and glob_min in lib_fortran.F90, see ticket:#1116

  • Property svn:keywords set to Id
File size: 38.2 KB
Line 
1MODULE limthd_dh
2   !!======================================================================
3   !!                       ***  MODULE limthd_dh ***
4   !!  LIM-3 :   thermodynamic growth and decay of the ice
5   !!======================================================================
6   !! History :  LIM  ! 2003-05 (M. Vancoppenolle) Original code in 1D
7   !!                 ! 2005-06 (M. Vancoppenolle) 3D version
8   !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdmsnif & rdmicif
9   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
10   !!----------------------------------------------------------------------
11#if defined key_lim3
12   !!----------------------------------------------------------------------
13   !!   'key_lim3'                                      LIM3 sea-ice model
14   !!----------------------------------------------------------------------
15   !!   lim_thd_dh  : vertical accr./abl. and lateral ablation of sea ice
16   !!----------------------------------------------------------------------
17   USE par_oce          ! ocean parameters
18   USE phycst           ! physical constants (OCE directory)
19   USE sbc_oce          ! Surface boundary condition: ocean fields
20   USE ice              ! LIM variables
21   USE par_ice          ! LIM parameters
22   USE thd_ice          ! LIM thermodynamics
23   USE in_out_manager   ! I/O manager
24   USE lib_mpp          ! MPP library
25   USE wrk_nemo         ! work arrays
26   USE lib_fortran      ! to use key_nosignedzero
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   lim_thd_dh   ! called by lim_thd
32
33   REAL(wp) ::   epsi20 = 1e-20   ! constant values
34   REAL(wp) ::   epsi13 = 1e-13   !
35   REAL(wp) ::   epsi16 = 1e-16   !
36   REAL(wp) ::   zzero  = 0.e0    !
37   REAL(wp) ::   zone   = 1.e0    !
38
39   !!----------------------------------------------------------------------
40   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
41   !! $Id$
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE lim_thd_dh( kideb, kiut, jl )
47      !!------------------------------------------------------------------
48      !!                ***  ROUTINE lim_thd_dh  ***
49      !!
50      !! ** Purpose :   determines variations of ice and snow thicknesses.
51      !!
52      !! ** Method  :   Ice/Snow surface melting arises from imbalance in surface fluxes
53      !!              Bottom accretion/ablation arises from flux budget
54      !!              Snow thickness can increase by precipitation and decrease by sublimation
55      !!              If snow load excesses Archmiede limit, snow-ice is formed by
56      !!              the flooding of sea-water in the snow
57      !!
58      !!                 1) Compute available flux of heat for surface ablation
59      !!                 2) Compute snow and sea ice enthalpies
60      !!                 3) Surface ablation and sublimation
61      !!                 4) Bottom accretion/ablation
62      !!                 5) Case of Total ablation
63      !!                 6) Snow ice formation
64      !!
65      !! References : Bitz and Lipscomb, 1999, J. Geophys. Res.
66      !!              Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
67      !!              Vancoppenolle, Fichefet and Bitz, 2005, Geophys. Res. Let.
68      !!              Vancoppenolle et al.,2009, Ocean Modelling
69      !!------------------------------------------------------------------
70      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
71      INTEGER , INTENT(in) ::   jl            ! Thickness cateogry number
72      !!
73      INTEGER  ::   ji , jk        ! dummy loop indices
74      INTEGER  ::   zji, zjj       ! 2D corresponding indices to ji
75      INTEGER  ::   iter
76
77      REAL(wp) ::   isnow          ! switch for presence (1) or absence (0) of snow
78      REAL(wp) ::   isnowic        ! snow ice formation not
79      REAL(wp) ::   i_ice_switch   ! ice thickness above a certain treshold or not
80
81      REAL(wp) ::   zzfmass_i, zihgnew                     ! local scalar
82      REAL(wp) ::   zzfmass_s, zhsnew, ztmelts             ! local scalar
83      REAL(wp) ::   zhn, zdhcf, zdhbf, zhni, zhnfi, zihg   !
84      REAL(wp) ::   zdhnm, zhnnew, zhisn, zihic, zzc       !
85      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
86      REAL(wp) ::   zds          ! increment of bottom ice salinity
87      REAL(wp) ::   zcoeff       ! dummy argument for snowfall partitioning over ice and leads
88      REAL(wp) ::   zsm_snowice  ! snow-ice salinity
89      REAL(wp) ::   zswi1        ! switch for computation of bottom salinity
90      REAL(wp) ::   zswi12       ! switch for computation of bottom salinity
91      REAL(wp) ::   zswi2        ! switch for computation of bottom salinity
92      REAL(wp) ::   zgrr         ! bottom growth rate
93      REAL(wp) ::   ztform       ! bottom formation temperature
94      !
95      REAL(wp), POINTER, DIMENSION(:) ::   zh_i        ! ice layer thickness
96      REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness
97      REAL(wp), POINTER, DIMENSION(:) ::   ztfs        ! melting point
98      REAL(wp), POINTER, DIMENSION(:) ::   zhsold      ! old snow thickness
99      REAL(wp), POINTER, DIMENSION(:) ::   zqprec      ! energy of fallen snow
100      REAL(wp), POINTER, DIMENSION(:) ::   zqfont_su   ! incoming, remaining surface energy
101      REAL(wp), POINTER, DIMENSION(:) ::   zqfont_bo   ! incoming, bottom energy
102      REAL(wp), POINTER, DIMENSION(:) ::   z_f_surf    ! surface heat for ablation
103      REAL(wp), POINTER, DIMENSION(:) ::   zhgnew      ! new ice thickness
104      REAL(wp), POINTER, DIMENSION(:) ::   zfmass_i    !
105
106      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel     ! snow melt
107      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_pre     ! snow precipitation
108      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_sub     ! snow sublimation
109      REAL(wp), POINTER, DIMENSION(:) ::   zfsalt_melt   ! salt flux due to ice melt
110
111      REAL(wp), POINTER, DIMENSION(:,:) ::   zdeltah
112
113      ! Pathological cases
114      REAL(wp), POINTER, DIMENSION(:) ::   zfdt_init   ! total incoming heat for ice melt
115      REAL(wp), POINTER, DIMENSION(:) ::   zfdt_final  ! total remaing heat for ice melt
116      REAL(wp), POINTER, DIMENSION(:) ::   zqt_i       ! total ice heat content
117      REAL(wp), POINTER, DIMENSION(:) ::   zqt_s       ! total snow heat content
118      REAL(wp), POINTER, DIMENSION(:) ::   zqt_dummy   ! dummy heat content
119
120      REAL(wp), POINTER, DIMENSION(:,:) ::   zqt_i_lay   ! total ice heat content
121
122      ! mass and salt flux (clem)
123      REAL(wp) :: zdvres, zdvsur, zdvbot
124      REAL(wp), POINTER, DIMENSION(:) ::   zviold, zvsold   ! old ice volume...
125
126      ! Heat conservation
127      INTEGER  ::   num_iter_max, numce_dh
128      REAL(wp) ::   meance_dh
129      REAL(wp), POINTER, DIMENSION(:) ::   zinnermelt
130      REAL(wp), POINTER, DIMENSION(:) ::   zfbase, zdq_i
131      !!------------------------------------------------------------------
132
133      CALL wrk_alloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i )
134      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfsalt_melt, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy )
135      CALL wrk_alloc( jpij, zinnermelt, zfbase, zdq_i )
136      CALL wrk_alloc( jpij, jkmax, zdeltah, zqt_i_lay )
137
138      CALL wrk_alloc( jpij, zviold, zvsold ) ! clem
139     
140      zfsalt_melt(:)  = 0._wp
141      ftotal_fin(:)   = 0._wp
142      zfdt_init(:)    = 0._wp
143      zfdt_final(:)   = 0._wp
144      dh_i_surf(:)    = 0._wp
145      dh_i_bott(:)    = 0._wp
146      dh_snowice(:)    = 0._wp
147
148      DO ji = kideb, kiut
149         old_ht_i_b(ji) = ht_i_b(ji)
150         old_ht_s_b(ji) = ht_s_b(ji)
151         zviold(ji) = a_i_b(ji) * ht_i_b(ji) ! clem
152         zvsold(ji) = a_i_b(ji) * ht_s_b(ji) ! clem
153      END DO
154      !
155      !------------------------------------------------------------------------------!
156      !  1) Calculate available heat for surface ablation                            !
157      !------------------------------------------------------------------------------!
158      !
159      DO ji = kideb, kiut
160         isnow         = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) )
161         ztfs(ji)      = isnow * rtt + ( 1.0 - isnow ) * rtt
162         z_f_surf(ji)  = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)
163         z_f_surf(ji)  = MAX( zzero , z_f_surf(ji) ) * MAX( zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) )
164         zfdt_init(ji) = ( z_f_surf(ji) + MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) ) * rdt_ice
165      END DO ! ji
166
167      zqfont_su  (:) = 0._wp
168      zqfont_bo  (:) = 0._wp
169      dsm_i_se_1d(:) = 0._wp     
170      dsm_i_si_1d(:) = 0._wp   
171      !
172      !------------------------------------------------------------------------------!
173      !  2) Computing layer thicknesses and  snow and sea-ice enthalpies.            !
174      !------------------------------------------------------------------------------!
175      !
176      DO ji = kideb, kiut     ! Layer thickness
177         zh_i(ji) = ht_i_b(ji) / REAL( nlay_i )
178         zh_s(ji) = ht_s_b(ji) / REAL( nlay_s )
179      END DO
180      !
181      zqt_s(:) = 0._wp        ! Total enthalpy of the snow
182      DO jk = 1, nlay_s
183         DO ji = kideb, kiut
184            zqt_s(ji) =  zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s )
185         END DO
186      END DO
187      !
188      zqt_i(:) = 0._wp        ! Total enthalpy of the ice
189      DO jk = 1, nlay_i
190         DO ji = kideb, kiut
191            zzc = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i )
192            zqt_i(ji)        =  zqt_i(ji) + zzc
193            zqt_i_lay(ji,jk) =              zzc
194         END DO
195      END DO
196      !
197      !------------------------------------------------------------------------------|
198      !  3) Surface ablation and sublimation                                         |
199      !------------------------------------------------------------------------------|
200      !
201      !-------------------------
202      ! 3.1 Snow precips / melt
203      !-------------------------
204      ! Snow accumulation in one thermodynamic time step
205      ! snowfall is partitionned between leads and ice
206      ! if snow fall was uniform, a fraction (1-at_i) would fall into leads
207      ! but because of the winds, more snow falls on leads than on sea ice
208      ! and a greater fraction (1-at_i)^beta of the total mass of snow
209      ! (beta < 1) falls in leads.
210      ! In reality, beta depends on wind speed,
211      ! and should decrease with increasing wind speed but here, it is
212      ! considered as a constant. an average value is 0.66
213      ! Martin Vancoppenolle, December 2006
214
215      ! Snow fall
216      DO ji = kideb, kiut
217         zcoeff = ( 1.0 - ( 1.0 - at_i_b(ji) )**betas ) / at_i_b(ji) 
218         zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn
219      END DO
220      zdh_s_mel(:) =  0._wp
221
222      ! Melt of fallen snow
223      DO ji = kideb, kiut
224         ! tatm_ice is now in K
225         zqprec   (ji)   =  rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus ) 
226         zqfont_su(ji)   =  z_f_surf(ji) * rdt_ice
227         zdeltah  (ji,1) =  MIN( 0.e0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) )
228         zqfont_su(ji)   =  MAX( 0.e0 , - zdh_s_pre(ji) - zdeltah(ji,1)              ) * zqprec(ji)
229         zdeltah  (ji,1) =  MAX( - zdh_s_pre(ji) , zdeltah(ji,1) )
230         zdh_s_mel(ji)   =  zdh_s_mel(ji) + zdeltah(ji,1)
231         ! heat conservation
232         qt_s_in(ji,jl)  =  qt_s_in(ji,jl) + zqprec(ji) * zdh_s_pre(ji)
233         zqt_s  (ji)     =  zqt_s  (ji)    + zqprec(ji) * zdh_s_pre(ji)
234         zqt_s  (ji)     =  MAX( zqt_s(ji) - zqfont_su(ji) , 0.e0 ) 
235      END DO
236
237
238      ! Snow melt due to surface heat imbalance
239      DO jk = 1, nlay_s
240         DO ji = kideb, kiut
241            zdeltah  (ji,jk) = - zqfont_su(ji) / q_s_b(ji,jk)
242            zqfont_su(ji)    =  MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * q_s_b(ji,jk) 
243            zdeltah  (ji,jk) =  MAX( zdeltah(ji,jk) , - zh_s(ji) )
244            zdh_s_mel(ji)    =  zdh_s_mel(ji) + zdeltah(ji,jk)        ! resulting melt of snow   
245         END DO
246      END DO
247
248      ! Apply snow melt to snow depth
249      DO ji = kideb, kiut
250         dh_s_tot(ji)   =  zdh_s_mel(ji) + zdh_s_pre(ji)
251         ! Old and new snow depths
252         zhsold(ji)     =  ht_s_b(ji)
253         zhsnew         =  ht_s_b(ji) + dh_s_tot(ji)
254         ! If snow is still present zhn = 1, else zhn = 0
255         zhn            =  1.0 - MAX( zzero , SIGN( zone , - zhsnew ) )
256         ht_s_b(ji)     =  MAX( zzero , zhsnew )
257         ! we recompute dh_s_tot (clem)
258         dh_s_tot (ji)  =  ht_s_b(ji) - zhsold(ji)
259         ! Volume and mass variations of snow
260         ! dvsbq_1d  (ji) =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_mel(ji) )
261         dvsbq_1d  (ji) =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_pre(ji) )   ! IOVINO
262         dvsbq_1d  (ji) =  MIN( zzero, dvsbq_1d(ji) )
263         !clem rdmsnif_1d(ji) =  rdmsnif_1d(ji) + rhosn * dvsbq_1d(ji)
264      END DO ! ji
265
266      !--------------------------
267      ! 3.2 Surface ice ablation
268      !--------------------------
269      DO ji = kideb, kiut 
270         z_f_surf (ji) =  zqfont_su(ji) / rdt_ice ! heat conservation test
271         zdq_i    (ji) =  0._wp
272      END DO ! ji
273
274      DO jk = 1, nlay_i
275         DO ji = kideb, kiut 
276            !                                                    ! melt of layer jk
277            zdeltah  (ji,jk) = - zqfont_su(ji) / q_i_b(ji,jk)
278            !                                                    ! recompute heat available
279            zqfont_su(ji)    = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk) 
280            !                                                    ! melt of layer jk cannot be higher than its thickness
281            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) )
282            !                                                    ! update surface melt
283            dh_i_surf(ji)    = dh_i_surf(ji) + zdeltah(ji,jk) 
284            !                                                    ! for energy conservation
285            zdq_i    (ji)    = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice
286            !
287            ! clem
288            fseqv_1d(ji) = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji)    &
289               &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice
290!            fseqv_1d(ji) = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic / rdt_ice
291         END DO
292      END DO
293
294      !                     !-------------------
295      IF( con_i ) THEN      ! Conservation test
296         !                  !-------------------
297         numce_dh  = 0
298         meance_dh = 0._wp
299         DO ji = kideb, kiut
300            IF ( ( z_f_surf(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN
301               numce_dh  = numce_dh + 1
302               meance_dh = meance_dh + z_f_surf(ji) + zdq_i(ji)
303            ENDIF
304            IF( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN!
305               WRITE(numout,*) ' ALERTE heat loss for surface melt '
306               WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl
307               WRITE(numout,*) ' ht_i_b       : ', ht_i_b(ji)
308               WRITE(numout,*) ' z_f_surf     : ', z_f_surf(ji)
309               WRITE(numout,*) ' zdq_i        : ', zdq_i(ji)
310               WRITE(numout,*) ' ht_i_b       : ', ht_i_b(ji)
311               WRITE(numout,*) ' fc_bo_i      : ', fc_bo_i(ji)
312               WRITE(numout,*) ' fbif_1d      : ', fbif_1d(ji)
313               WRITE(numout,*) ' qlbbq_1d     : ', qlbbq_1d(ji)
314               WRITE(numout,*) ' s_i_new      : ', s_i_new(ji)
315               WRITE(numout,*) ' sss_m        : ', sss_m(zji,zjj)
316            ENDIF
317         END DO
318         !
319         IF( numce_dh > 0 )   meance_dh = meance_dh / numce_dh
320         WRITE(numout,*) ' Error report - Category : ', jl
321         WRITE(numout,*) ' ~~~~~~~~~~~~ '
322         WRITE(numout,*) ' Number of points where there is sur. me. error : ', numce_dh
323         WRITE(numout,*) ' Mean basal growth error on error points : ', meance_dh
324         !
325      ENDIF
326
327      !----------------------
328      ! 3.3 Snow sublimation
329      !----------------------
330
331      DO ji = kideb, kiut
332         ! if qla is positive (upwards), heat goes to the atmosphere, therefore
333         ! snow sublimates, if qla is negative (downwards), snow condensates
334         zdh_s_sub(ji)    =  - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice
335         dh_s_tot (ji)    =  dh_s_tot(ji) + zdh_s_sub(ji)
336         zdhcf            =  ht_s_b(ji) + zdh_s_sub(ji) 
337         ht_s_b   (ji)    =  MAX( zzero , zdhcf )
338         ! we recompute dh_s_tot
339         dh_s_tot (ji)    =  ht_s_b(ji) - zhsold(ji)
340         qt_s_in  (ji,jl) =  qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1)
341      END DO
342
343      zqt_dummy(:) = 0.e0
344      DO jk = 1, nlay_s
345         DO ji = kideb,kiut
346            q_s_b    (ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
347            zqt_dummy(ji)    =  zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s )            ! heat conservation
348         END DO
349      END DO
350
351      DO jk = 1, nlay_s
352         DO ji = kideb, kiut
353            ! In case of disparition of the snow, we have to update the snow temperatures
354            zhisn  =  MAX( zzero , SIGN( zone, - ht_s_b(ji) ) )
355            t_s_b(ji,jk) = ( 1.0 - zhisn ) * t_s_b(ji,jk) + zhisn * rtt
356            q_s_b(ji,jk) = ( 1.0 - zhisn ) * q_s_b(ji,jk)
357         END DO
358      END DO
359
360      !
361      !------------------------------------------------------------------------------!
362      ! 4) Basal growth / melt                                                       !
363      !------------------------------------------------------------------------------!
364      !
365      ! Ice basal growth / melt is given by the ratio of heat budget over basal
366      ! ice heat content.  Basal heat budget is given by the difference between
367      ! the inner conductive flux  (fc_bo_i), from the open water heat flux
368      ! (qlbbqb) and the turbulent ocean flux (fbif).
369      ! fc_bo_i is positive downwards. fbif and qlbbq are positive to the ice
370
371      !-----------------------------------------------------
372      ! 4.1 Basal growth - (a) salinity not varying in time
373      !-----------------------------------------------------
374      IF(  num_sal /= 2  .AND.  num_sal /= 4  ) THEN
375         DO ji = kideb, kiut
376            IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0.0  ) THEN
377               s_i_new(ji)         =  sm_i_b(ji)
378               ! Melting point in K
379               ztmelts             =  - tmut * s_i_new(ji) + rtt 
380               ! New ice heat content (Bitz and Lipscomb, 1999)
381               ztform              =  t_i_b(ji,nlay_i)  ! t_bo_b crashes in the
382               ! Baltic
383               q_i_b(ji,nlay_i+1)  = rhoic * (  cpic * ( ztmelts - ztform )                                &
384                  &                           + lfus * (  1.0 - ( ztmelts - rtt ) / ( ztform - rtt )  )    &
385                  &                           - rcp  * ( ztmelts - rtt )                                 )
386               ! Basal growth rate = - F*dt / q
387               dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
388               fseqv_1d(ji) = fseqv_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic / rdt_ice
389            ENDIF
390         END DO
391      ENDIF
392
393      !-------------------------------------------------
394      ! 4.1 Basal growth - (b) salinity varying in time
395      !-------------------------------------------------
396      IF(  num_sal == 2 .OR.  num_sal == 4  ) THEN
397         ! the growth rate (dh_i_bott) is function of the new ice
398         ! heat content (q_i_b(nlay_i+1)). q_i_b depends on the new ice
399         ! salinity (snewice). snewice depends on dh_i_bott
400         ! it converges quickly, so, no problem
401         ! See Vancoppenolle et al., OM08 for more info on this
402
403         ! Initial value (tested 1D, can be anything between 1 and 20)
404         num_iter_max = 4
405         s_i_new(:)   = 4.0
406
407         ! Iterative procedure
408         DO iter = 1, num_iter_max
409            DO ji = kideb, kiut
410               IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0  ) THEN
411                  zji = MOD( npb(ji) - 1, jpi ) + 1
412                  zjj = ( npb(ji) - 1 ) / jpi + 1
413                  ! Melting point in K
414                  ztmelts             =   - tmut * s_i_new(ji) + rtt 
415                  ! New ice heat content (Bitz and Lipscomb, 1999)
416                  q_i_b(ji,nlay_i+1)  =  rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             &
417                     &                            + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )   &
418                     &                            - rcp * ( ztmelts-rtt )                                     )
419                  ! Bottom growth rate = - F*dt / q
420                  dh_i_bott(ji) =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)
421                  ! New ice salinity ( Cox and Weeks, JGR, 1988 )
422                  ! zswi2  (1) if dh_i_bott/rdt .GT. 3.6e-7
423                  ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8
424                  ! zswi1  (1) if dh_i_bott/rdt .LT. 2.0e-8
425                  zgrr   = MIN( 1.0e-3, MAX ( dh_i_bott(ji) / rdt_ice , epsi13 ) )
426                  zswi2  = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) ) 
427                  zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
428                  zswi1  = 1. - zswi2 * zswi12 
429                  zfracs = zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
430                     &                   + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 
431                  zds         = zfracs * sss_m(zji,zjj) - s_i_new(ji)
432                  s_i_new(ji) = zfracs * sss_m(zji,zjj)
433               ENDIF ! fc_bo_i
434            END DO ! ji
435         END DO ! iter
436
437         ! Final values
438         DO ji = kideb, kiut
439            IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN
440               ! New ice salinity must not exceed 15 psu
441               s_i_new(ji) = MIN( s_i_new(ji), s_i_max )
442               ! Metling point in K
443               ztmelts     =   - tmut * s_i_new(ji) + rtt 
444               ! New ice heat content (Bitz and Lipscomb, 1999)
445               q_i_b(ji,nlay_i+1)  =  rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                              &
446                  &                            + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )    &
447                  &                            - rcp * ( ztmelts - rtt )                                    )
448               ! Basal growth rate = - F*dt / q
449               dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
450               ! Salinity update
451               ! entrapment during bottom growth
452               !clem:movedown dsm_i_se_1d(ji) = ( s_i_new(ji) * dh_i_bott(ji) + sm_i_b(ji) * ht_i_b(ji) )    &
453               !clem:movedown   &            / MAX( ht_i_b(ji) + dh_i_bott(ji) ,epsi13 ) - sm_i_b(ji)
454               fseqv_1d(ji) = fseqv_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic / rdt_ice
455            ENDIF ! heat budget
456         END DO
457      ENDIF
458
459      !----------------
460      ! 4.2 Basal melt
461      !----------------
462      meance_dh = 0._wp
463      numce_dh  = 0
464      zinnermelt(:) = 0._wp
465
466      DO ji = kideb, kiut
467         ! heat convergence at the surface > 0
468         IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0._wp  ) THEN
469            s_i_new(ji)   =  s_i_b(ji,nlay_i)
470            zqfont_bo(ji) =  rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) )
471            zfbase(ji)    =  zqfont_bo(ji) / rdt_ice     ! heat conservation test
472            zdq_i(ji)     =  0._wp
473            dh_i_bott(ji) =  0._wp
474         ENDIF
475      END DO
476
477      DO jk = nlay_i, 1, -1
478         DO ji = kideb, kiut
479            IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN
480               ztmelts            =   - tmut * s_i_b(ji,jk) + rtt 
481               IF( t_i_b(ji,jk) >= ztmelts ) THEN
482                  zdeltah(ji,jk)  = - zh_i(ji)
483                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)
484                  zinnermelt(ji)   = 1._wp
485               ELSE  ! normal ablation
486                  zdeltah(ji,jk)  = - zqfont_bo(ji) / q_i_b(ji,jk)
487                  zqfont_bo(ji)   = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk)
488                  zdeltah(ji,jk)  = MAX(zdeltah(ji,jk), - zh_i(ji) )
489                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)
490                  zdq_i(ji)       = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice
491               ENDIF
492               ! clem
493               fseqv_1d(ji) = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji)    &
494                    &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice
495!               fseqv_1d(ji) = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic / rdt_ice
496            ENDIF
497         END DO ! ji
498      END DO ! jk
499
500      !                     !-------------------
501      IF( con_i ) THEN      ! Conservation test
502      !                     !-------------------
503         DO ji = kideb, kiut
504            IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0.e0  ) THEN
505               IF( ( zfbase(ji) + zdq_i(ji) ) >= 1.e-3 ) THEN
506                  numce_dh  = numce_dh + 1
507                  meance_dh = meance_dh + zfbase(ji) + zdq_i(ji)
508               ENDIF
509               IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN
510                  WRITE(numout,*) ' ALERTE heat loss for basal melt : zji, zjj, jl :', zji, zjj, jl
511                  WRITE(numout,*) ' ht_i_b    : ', ht_i_b(ji)
512                  WRITE(numout,*) ' zfbase    : ', zfbase(ji)
513                  WRITE(numout,*) ' zdq_i     : ', zdq_i(ji)
514                  WRITE(numout,*) ' ht_i_b    : ', ht_i_b(ji)
515                  WRITE(numout,*) ' fc_bo_i   : ', fc_bo_i(ji)
516                  WRITE(numout,*) ' fbif_1d   : ', fbif_1d(ji)
517                  WRITE(numout,*) ' qlbbq_1d  : ', qlbbq_1d(ji)
518                  WRITE(numout,*) ' s_i_new   : ', s_i_new(ji)
519                  WRITE(numout,*) ' sss_m     : ', sss_m(zji,zjj)
520                  WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
521                  WRITE(numout,*) ' innermelt : ', INT( zinnermelt(ji) )
522               ENDIF
523            ENDIF
524         END DO
525         IF( numce_dh > 0 )   meance_dh = meance_dh / numce_dh
526         WRITE(numout,*) ' Number of points where there is bas. me. error : ', numce_dh
527         WRITE(numout,*) ' Mean basal melt error on error points : ', meance_dh
528         WRITE(numout,*) ' Remaining bottom heat : ', zqfont_bo(jiindex_1d)
529         !
530      ENDIF
531
532      !
533      !------------------------------------------------------------------------------!
534      !  5) Pathological cases                                                       !
535      !------------------------------------------------------------------------------!
536      !
537      !----------------------------------------------
538      ! 5.1 Excessive ablation in a 1-category model
539      !----------------------------------------------
540
541      DO ji = kideb, kiut
542         !                     ! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15)
543         IF( jpl == 1 ) THEN   ;   zdhbf = MAX( hmelt , dh_i_bott(ji) )
544         ELSE                  ;   zdhbf =              dh_i_bott(ji) 
545         ENDIF
546         zdvres        = zdhbf - dh_i_bott(ji)
547         dh_i_bott(ji) = zdhbf
548         fseqv_1d(ji)  = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic / rdt_ice
549         !                     ! excessive energy is sent to lateral ablation
550         fsup     (ji) =  rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) * zdvres / rdt_ice
551      END DO
552
553      !-----------------------------------
554      ! 5.2 More than available ice melts
555      !-----------------------------------
556      ! then heat applied minus heat content at previous time step
557      ! should equal heat remaining
558      !
559      DO ji = kideb, kiut
560         ! Adapt the remaining energy if too much ice melts
561         !--------------------------------------------------
562         zdvres     = MAX( 0._wp, - ht_i_b(ji) - dh_i_surf(ji) - dh_i_bott(ji) )
563         zdvsur     = MIN( 0._wp, dh_i_surf(ji) + zdvres ) - dh_i_surf(ji) ! fill the surface first
564         zdvbot     = MAX( 0._wp, zdvres - zdvsur ) ! then the bottom
565         dh_i_surf (ji) = dh_i_surf(ji) + zdvsur ! clem
566         dh_i_bott (ji) = dh_i_bott(ji) + zdvbot ! clem
567
568         ! new ice thickness (clem)
569         zhgnew(ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji)
570         zihgnew    = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice
571         zhgnew(ji) = zihgnew * zhgnew(ji)      ! ice thickness is put to 0
572 
573         !                     !since ice volume is only used for outputs, we keep it global for all categories
574         dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji)
575
576        ! remaining heat
577         zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) +  zqfont_bo(ji) ) 
578
579         ! If snow remains, energy is used to melt snow
580         zhni =  ht_s_b(ji)      ! snow depth at previous time step
581         zihg =  MAX( zzero , SIGN ( zone , - ht_s_b(ji) ) ) ! 0 if snow
582
583         ! energy of melting of remaining snow
584         zqt_s(ji) =    ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi13 )
585         zdhnm     =  - ( 1. - zihg ) * ( 1. - zihgnew ) * zfdt_final(ji) / MAX( zqt_s(ji) , epsi13 )
586         zhnfi          =  zhni + zdhnm
587         zfdt_final(ji) =  MAX( zfdt_final(ji) + zqt_s(ji) * zdhnm , 0.0 )
588         ht_s_b(ji)     =  MAX( zzero , zhnfi )
589         zqt_s(ji)      =  zqt_s(ji) * ht_s_b(ji)
590         ! we recompute dh_s_tot (clem)
591         dh_s_tot (ji)  =  ht_s_b(ji) - zhsold(ji)
592
593         ! Mass variations of ice and snow
594         !---------------------------------
595         !                                              ! mass variation of the jl category
596         zzfmass_s = - a_i_b(ji) * ( zhni       - ht_s_b(ji) ) * rhosn   ! snow
597         zzfmass_i =   a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic   ! ice 
598         !
599         zfmass_i(ji) = zzfmass_i                       ! ice variation saved to compute salt flux (see below)
600         !
601         !                                              ! mass variation cumulated over category
602         !clem rdmsnif_1d(ji) = rdmsnif_1d(ji) + zzfmass_s                     ! snow
603         !clem rdmicif_1d(ji) = rdmicif_1d(ji) + zzfmass_i                     ! ice
604
605         ! Remaining heat to the ocean
606         !---------------------------------
607         focea(ji)  = - zfdt_final(ji) / rdt_ice         ! focea is in W.m-2 * dt
608
609         ! residual salt flux (clem)
610         !--------------------------
611         ! surface
612         fseqv_1d(ji)  = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvsur * rhoic / rdt_ice
613         ! bottom
614         IF     ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0._wp ) THEN ! melting
615            fseqv_1d(ji)  = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvbot * rhoic / rdt_ice
616         ELSE ! growth
617            fseqv_1d(ji)  = fseqv_1d(ji) - s_i_new(ji) * a_i_b(ji) * zdvbot * rhoic / rdt_ice
618         ENDIF
619         !
620         !                     ! diagnostic ( bottom ice growth )
621         zji = MOD( npb(ji) - 1, jpi ) + 1
622         zjj = ( npb(ji) - 1 ) / jpi + 1
623         diag_bot_gr(zji,zjj) = diag_bot_gr(zji,zjj) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) / rdt_ice
624         diag_sur_me(zji,zjj) = diag_sur_me(zji,zjj) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) / rdt_ice
625         diag_bot_me(zji,zjj) = diag_bot_me(zji,zjj) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) / rdt_ice
626      END DO
627
628      ftotal_fin (:) = zfdt_final(:)  / rdt_ice
629
630      !---------------------------
631      ! Salt flux and heat fluxes                   
632      !---------------------------
633      DO ji = kideb, kiut
634         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   !1 if ice
635
636         ! Salt flux
637         zji = MOD( npb(ji) - 1, jpi ) + 1
638         zjj = ( npb(ji) - 1 ) / jpi + 1
639         ! new lines
640         !IF( num_sal == 4 ) THEN
641         !   ! IOVINO
642         !   fseqv_1d(ji) = fseqv_1d(ji) +        zihgnew  * zfsalt_melt(ji)                                &
643         !      &                        - (1.0 - zihgnew) * zfmass_i(ji) * bulk_sal / rdt_ice
644         !ELSE
645         !   ! IOVINO
646         !   fseqv_1d(ji) = fseqv_1d(ji) +        zihgnew  * zfsalt_melt(ji)                                &
647         !      &                        - (1.0 - zihgnew) * zfmass_i(ji) * sm_i_b(ji) / rdt_ice
648         !ENDIF
649         ! Heat flux
650         ! excessive bottom ablation energy (fsup) - 0 except if jpl = 1
651         ! excessive total ablation energy (focea) sent to the ocean
652         qfvbq_1d(ji)  = qfvbq_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) * rdt_ice
653
654         zihic   = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) )
655         ! equals 0 if ht_i = 0, 1 if ht_i gt 0
656         fscbq_1d(ji) =  a_i_b(ji) * fstbif_1d(ji)
657         qldif_1d(ji)  = qldif_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji)    * a_i_b(ji) * rdt_ice   &
658            &                                    + ( 1.0 - zihic   ) * fscbq_1d(ji)             * rdt_ice
659         !IF ( (zji.eq.jiindx).AND.(zjj.eq.jjindx) ) THEN
660         !   !clemclem
661         !   WRITE(numout,*) 'lim_thd_dh : qldif = ', qldif_1d(ji)
662         !   !clemclem
663         !ENDIF
664      END DO  ! ji
665
666      !-------------------------------------------
667      ! Correct temperature, energy and thickness
668      !-------------------------------------------
669      DO ji = kideb, kiut
670         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 
671         t_su_b(ji) =  zihgnew * t_su_b(ji) + ( 1.0 - zihgnew ) * rtt
672      END DO  ! ji
673
674      DO jk = 1, nlay_i
675         DO ji = kideb, kiut
676            zihgnew      =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 
677            t_i_b(ji,jk) =  zihgnew * t_i_b(ji,jk) + ( 1.0 - zihgnew ) * rtt
678            q_i_b(ji,jk) =  zihgnew * q_i_b(ji,jk)
679         END DO
680      END DO  ! ji
681
682      DO ji = kideb, kiut
683         ht_i_b(ji) = zhgnew(ji)
684      END DO  ! ji
685      !
686      !------------------------------------------------------------------------------|
687      !  6) Snow-Ice formation                                                       |
688      !------------------------------------------------------------------------------|
689      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
690      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
691      DO ji = kideb, kiut
692         !
693         dh_snowice(ji) = MAX(  zzero , ( rhosn * ht_s_b(ji) + (rhoic-rau0) * ht_i_b(ji) ) / ( rhosn+rau0-rhoic )  )
694         zhgnew(ji)     = MAX(  zhgnew(ji) , zhgnew(ji) + dh_snowice(ji)  )
695         zhnnew         = MIN(  ht_s_b(ji) , ht_s_b(ji) - dh_snowice(ji)  )
696
697         !  Changes in ice volume and ice mass.
698         dvnbq_1d  (ji) =                a_i_b(ji) * ( zhgnew(ji)-ht_i_b(ji) )
699         dmgwi_1d  (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn
700
701         !clem rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic
702         !clem rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) * ( zhnnew     - ht_s_b(ji) ) * rhosn
703
704         !        Equivalent salt flux (1) Snow-ice formation component
705         !        -----------------------------------------------------
706         zji = MOD( npb(ji) - 1, jpi ) + 1
707         zjj =    ( npb(ji) - 1 ) / jpi + 1
708
709         IF( num_sal /= 2 ) THEN   ;   zsm_snowice = sm_i_b(ji)
710         ELSE                      ;   zsm_snowice = ( rhoic - rhosn ) / rhoic * sss_m(zji,zjj) 
711         ENDIF
712         !IF( num_sal == 4 ) THEN 
713         !   ! IOVINO
714         !   fseqv_1d(ji) = fseqv_1d(ji) - bulk_sal * a_i_b(ji) * dh_snowice(ji) * rhoic / rdt_ice
715         !ELSE
716         !   ! IOVINO
717         !   fseqv_1d(ji) = fseqv_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic / rdt_ice
718         !ENDIF
719         ! entrapment during snow ice formation
720
721         !clem IF(  num_sal == 2  .OR.  num_sal == 4  )   &
722         !clem   dsm_i_si_1d(ji) = ( zsm_snowice*dh_snowice(ji) &
723         !clem   &               + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13)   &
724         !clem   &               - sm_i_b(ji) ) * isnowic
725
726         ! clem: new salinity difference stored (to be used in limthd_ent.F90)
727         IF (  num_sal == 2  .OR.  num_sal == 4  ) THEN
728            i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - zhgnew(ji) + epsi13 ) )
729            ! salinity dif due to snow-ice formation
730            dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch     
731            ! salinity dif due to bottom growth
732            IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0._wp ) THEN
733               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch
734            ENDIF
735         ENDIF
736
737         !  Actualize new snow and ice thickness.
738         ht_s_b(ji)  = zhnnew
739         ht_i_b(ji)  = zhgnew(ji)
740
741         ! Total ablation ! new lines added to debug
742         IF( ht_i_b(ji) <= 0._wp )   a_i_b(ji) = 0._wp
743
744         ! diagnostic ( snow ice growth )
745         zji = MOD( npb(ji) - 1, jpi ) + 1
746         zjj =    ( npb(ji) - 1 ) / jpi + 1
747         diag_sni_gr(zji,zjj)  = diag_sni_gr(zji,zjj) + dh_snowice(ji)*a_i_b(ji) / rdt_ice
748
749         ! salt flux
750         fseqv_1d(ji) = fseqv_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic / rdt_ice
751         !--------------------------------
752         ! Update mass fluxes (clem)
753         !--------------------------------
754         rdmicif_1d(ji) = rdmicif_1d(ji) + ( a_i_b(ji) * ht_i_b(ji) - zviold(ji) ) * rhoic 
755         rdmsnif_1d(ji) = rdmsnif_1d(ji) + ( a_i_b(ji) * ht_s_b(ji) - zvsold(ji) ) * rhosn 
756
757      END DO !ji
758      !
759      CALL wrk_dealloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i )
760      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfsalt_melt, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy )
761      CALL wrk_dealloc( jpij, zinnermelt, zfbase, zdq_i )
762      CALL wrk_dealloc( jpij, jkmax, zdeltah, zqt_i_lay )
763      !
764      CALL wrk_dealloc( jpij, zviold, zvsold ) ! clem
765      !
766   END SUBROUTINE lim_thd_dh
767   
768#else
769   !!----------------------------------------------------------------------
770   !!   Default option                               NO  LIM3 sea-ice model
771   !!----------------------------------------------------------------------
772CONTAINS
773   SUBROUTINE lim_thd_dh          ! Empty routine
774   END SUBROUTINE lim_thd_dh
775#endif
776
777   !!======================================================================
778END MODULE limthd_dh
Note: See TracBrowser for help on using the repository browser.