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/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 6 years ago

Remove svn keywords

File size: 36.7 KB
RevLine 
[825]1MODULE limthd_dh
[1572]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
[4688]8   !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in wfx_snw & wfx_ice
[3625]9   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
10   !!            3.5  ! 2012-10 (G. Madec & co) salt flux + bug fixes
[1572]11   !!----------------------------------------------------------------------
[825]12#if defined key_lim3
[834]13   !!----------------------------------------------------------------------
14   !!   'key_lim3'                                      LIM3 sea-ice model
15   !!----------------------------------------------------------------------
[3625]16   !!   lim_thd_dh    : vertical accr./abl. and lateral ablation of sea ice
[825]17   !!----------------------------------------------------------------------
[3625]18   USE par_oce        ! ocean parameters
19   USE phycst         ! physical constants (OCE directory)
20   USE sbc_oce        ! Surface boundary condition: ocean fields
21   USE ice            ! LIM variables
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    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
[4688]27   
[825]28   IMPLICIT NONE
29   PRIVATE
30
[5407]31   PUBLIC   lim_thd_dh      ! called by lim_thd
[5487]32   PUBLIC   lim_thd_snwblow ! called in sbcblk/sbcclio/sbccpl and here
[825]33
[5407]34   INTERFACE lim_thd_snwblow
35      MODULE PROCEDURE lim_thd_snwblow_1d, lim_thd_snwblow_2d
36   END INTERFACE
37
[825]38   !!----------------------------------------------------------------------
[4161]39   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
[1156]40   !! $Id$
[2715]41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[825]42   !!----------------------------------------------------------------------
43CONTAINS
44
[4688]45   SUBROUTINE lim_thd_dh( kideb, kiut )
[921]46      !!------------------------------------------------------------------
47      !!                ***  ROUTINE lim_thd_dh  ***
48      !!
[1572]49      !! ** Purpose :   determines variations of ice and snow thicknesses.
[921]50      !!
[1572]51      !! ** Method  :   Ice/Snow surface melting arises from imbalance in surface fluxes
52      !!              Bottom accretion/ablation arises from flux budget
53      !!              Snow thickness can increase by precipitation and decrease by sublimation
54      !!              If snow load excesses Archmiede limit, snow-ice is formed by
55      !!              the flooding of sea-water in the snow
[921]56      !!
[1572]57      !!                 1) Compute available flux of heat for surface ablation
58      !!                 2) Compute snow and sea ice enthalpies
59      !!                 3) Surface ablation and sublimation
60      !!                 4) Bottom accretion/ablation
61      !!                 5) Case of Total ablation
62      !!                 6) Snow ice formation
[921]63      !!
[1572]64      !! References : Bitz and Lipscomb, 1999, J. Geophys. Res.
65      !!              Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
66      !!              Vancoppenolle, Fichefet and Bitz, 2005, Geophys. Res. Let.
67      !!              Vancoppenolle et al.,2009, Ocean Modelling
[921]68      !!------------------------------------------------------------------
[1572]69      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
70      !!
71      INTEGER  ::   ji , jk        ! dummy loop indices
[4161]72      INTEGER  ::   ii, ij         ! 2D corresponding indices to ji
[1572]73      INTEGER  ::   iter
[825]74
[4688]75      REAL(wp) ::   ztmelts             ! local scalar
[5123]76      REAL(wp) ::   zfdum       
[1572]77      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
[5407]78      REAL(wp) ::   zs_snic      ! snow-ice salinity
[1572]79      REAL(wp) ::   zswi1        ! switch for computation of bottom salinity
80      REAL(wp) ::   zswi12       ! switch for computation of bottom salinity
81      REAL(wp) ::   zswi2        ! switch for computation of bottom salinity
82      REAL(wp) ::   zgrr         ! bottom growth rate
[4688]83      REAL(wp) ::   zt_i_new     ! bottom formation temperature
84
85      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
86      REAL(wp) ::   zEi          ! specific enthalpy of sea ice (J/kg)
87      REAL(wp) ::   zEw          ! specific enthalpy of exchanged water (J/kg)
88      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg)
89      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean
90      REAL(wp) ::   zsstK        ! SST in Kelvin
91
92      REAL(wp), POINTER, DIMENSION(:) ::   zqprec      ! energy of fallen snow                       (J.m-3)
93      REAL(wp), POINTER, DIMENSION(:) ::   zq_su       ! heat for surface ablation                   (J.m-2)
94      REAL(wp), POINTER, DIMENSION(:) ::   zq_bo       ! heat for bottom ablation                    (J.m-2)
95      REAL(wp), POINTER, DIMENSION(:) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2)
[5123]96      REAL(wp), POINTER, DIMENSION(:) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2)
[3294]97
[3625]98      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel   ! snow melt
99      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_pre   ! snow precipitation
100      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_sub   ! snow sublimation
[3294]101
102      REAL(wp), POINTER, DIMENSION(:,:) ::   zdeltah
[4688]103      REAL(wp), POINTER, DIMENSION(:,:) ::   zh_i      ! ice layer thickness
[5202]104      INTEGER , POINTER, DIMENSION(:,:) ::   icount    ! number of layers vanished by melting
[3294]105
[4688]106      REAL(wp), POINTER, DIMENSION(:) ::   zqh_i       ! total ice heat content  (J.m-2)
107      REAL(wp), POINTER, DIMENSION(:) ::   zqh_s       ! total snow heat content (J.m-2)
108      REAL(wp), POINTER, DIMENSION(:) ::   zq_s        ! total snow enthalpy     (J.m-3)
[5407]109      REAL(wp), POINTER, DIMENSION(:) ::   zsnw        ! distribution of snow after wind blowing
[3294]110
[5123]111      REAL(wp) :: zswitch_sal
[4161]112
[3294]113      ! Heat conservation
[4688]114      INTEGER  ::   num_iter_max
115
[1572]116      !!------------------------------------------------------------------
[825]117
[5123]118      ! Discriminate between varying salinity (nn_icesal=2) and prescribed cases (other values)
119      SELECT CASE( nn_icesal )                       ! varying salinity or not
[4688]120         CASE( 1, 3, 4 ) ;   zswitch_sal = 0       ! prescribed salinity profile
121         CASE( 2 )       ;   zswitch_sal = 1       ! varying salinity profile
122      END SELECT
[825]123
[5487]124      CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw )
125      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
[5202]126      CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i )
127      CALL wrk_alloc( jpij, nlay_i, icount )
[5487]128       
[4688]129      dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp
130      dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp   
[2715]131
[5487]132      zqprec   (:) = 0._wp ; zq_su    (:) = 0._wp ; zq_bo    (:) = 0._wp ; zf_tt(:) = 0._wp
133      zq_rema  (:) = 0._wp ; zsnw     (:) = 0._wp
134      zdh_s_mel(:) = 0._wp ; zdh_s_pre(:) = 0._wp ; zdh_s_sub(:) = 0._wp ; zqh_i(:) = 0._wp
135      zqh_s    (:) = 0._wp ; zq_s     (:) = 0._wp     
[4161]136
[5487]137      zdeltah(:,:) = 0._wp ; zh_i(:,:) = 0._wp       
138      icount (:,:) = 0
[4688]139
[5487]140
[5167]141      ! Initialize enthalpy at nlay_i+1
142      DO ji = kideb, kiut
143         q_i_1d(ji,nlay_i+1) = 0._wp
144      END DO
145
[4688]146      ! initialize layer thicknesses and enthalpies
147      h_i_old (:,0:nlay_i+1) = 0._wp
148      qh_i_old(:,0:nlay_i+1) = 0._wp
149      DO jk = 1, nlay_i
150         DO ji = kideb, kiut
[5123]151            h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i
[4872]152            qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk)
[4688]153         ENDDO
154      ENDDO
[921]155      !
156      !------------------------------------------------------------------------------!
[4688]157      !  1) Calculate available heat for surface and bottom ablation                 !
[921]158      !------------------------------------------------------------------------------!
159      !
[2715]160      DO ji = kideb, kiut
[4990]161         zfdum      = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
162         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 
[4688]163
[5167]164         zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
[4688]165         zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
166      END DO
167
[921]168      !
169      !------------------------------------------------------------------------------!
[4688]170      ! If snow temperature is above freezing point, then snow melts
171      ! (should not happen but sometimes it does)
[921]172      !------------------------------------------------------------------------------!
[4688]173      DO ji = kideb, kiut
[5123]174         IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting
[4688]175            ! Contribution to heat flux to the ocean [W.m-2], < 0 
[4872]176            hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
[4688]177            ! Contribution to mass flux
[4872]178            wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
[4688]179            ! updates
[4872]180            ht_s_1d(ji)   = 0._wp
181            q_s_1d (ji,1) = 0._wp
[5123]182            t_s_1d (ji,1) = rt0
[4688]183         END IF
184      END DO
185
186      !------------------------------------------------------------!
187      !  2) Computing layer thicknesses and enthalpies.            !
188      !------------------------------------------------------------!
[921]189      !
[825]190      DO jk = 1, nlay_s
[2715]191         DO ji = kideb, kiut
[5167]192            zqh_s(ji) =  zqh_s(ji) + q_s_1d(ji,jk) * ht_s_1d(ji) * r1_nlay_s
[825]193         END DO
194      END DO
[2715]195      !
[825]196      DO jk = 1, nlay_i
[2715]197         DO ji = kideb, kiut
[5123]198            zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i
[4872]199            zqh_i(ji)   = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk)
[825]200         END DO
201      END DO
[921]202      !
203      !------------------------------------------------------------------------------|
204      !  3) Surface ablation and sublimation                                         |
205      !------------------------------------------------------------------------------|
206      !
[834]207      !-------------------------
208      ! 3.1 Snow precips / melt
209      !-------------------------
[825]210      ! Snow accumulation in one thermodynamic time step
211      ! snowfall is partitionned between leads and ice
212      ! if snow fall was uniform, a fraction (1-at_i) would fall into leads
213      ! but because of the winds, more snow falls on leads than on sea ice
214      ! and a greater fraction (1-at_i)^beta of the total mass of snow
[834]215      ! (beta < 1) falls in leads.
[825]216      ! In reality, beta depends on wind speed,
217      ! and should decrease with increasing wind speed but here, it is
[834]218      ! considered as a constant. an average value is 0.66
[825]219      ! Martin Vancoppenolle, December 2006
220
[5487]221      CALL lim_thd_snwblow( 1. - at_i_1d(kideb:kiut), zsnw(kideb:kiut) ) ! snow distribution over ice after wind blowing
222
[5167]223      zdeltah(:,:) = 0._wp
[825]224      DO ji = kideb, kiut
[4688]225         !-----------
226         ! Snow fall
227         !-----------
228         ! thickness change
[5407]229         zdh_s_pre(ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji)
230         ! enthalpy of the precip (>0, J.m-3)
231         zqprec   (ji) = - qprec_ice_1d(ji)   
[4688]232         IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp
233         ! heat flux from snow precip (>0, W.m-2)
[4872]234         hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
[4688]235         ! mass flux, <0
[4872]236         wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice
[825]237
[4688]238         !---------------------
239         ! Melt of falling snow
240         !---------------------
241         ! thickness change
[5134]242         rswitch        = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )
[5167]243         zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 )
244         zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting
[4688]245         ! heat used to melt snow (W.m-2, >0)
[5167]246         hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
[4688]247         ! snow melting only = water into the ocean (then without snow precip), >0
[5167]248         wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice   
249         ! updates available heat + precipitations after melting
250         zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )     
251         zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1)
[4688]252
[5167]253         ! update thickness
254         ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) )
[825]255      END DO
256
[5167]257      ! If heat still available (zq_su > 0), then melt more snow
258      zdeltah(:,:) = 0._wp
[825]259      DO jk = 1, nlay_s
260         DO ji = kideb, kiut
[4688]261            ! thickness change
[4990]262            rswitch          = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 
[5134]263            rswitch          = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) ) 
[4990]264            zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 )
[5167]265            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - ht_s_1d(ji) ) ! bound melting
[4688]266            zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)   
267            ! heat used to melt snow(W.m-2, >0)
[4872]268            hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice 
[4688]269            ! snow melting only = water into the ocean (then without snow precip)
[4872]270            wfx_snw_1d(ji)   = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
[4688]271            ! updates available heat + thickness
[5123]272            zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) )
[4872]273            ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) )
[1572]274         END DO
275      END DO
[825]276
[4688]277      !----------------------
278      ! 3.2 Snow sublimation
279      !----------------------
280      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
[5167]281      ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean)
[4688]282      ! clem comment: ice should also sublimate
[5167]283      zdeltah(:,:) = 0._wp
[5407]284      ! coupled mode: sublimation is set to 0 (evap_ice = 0) until further notice
285      ! forced  mode: snow thickness change due to sublimation
286      DO ji = kideb, kiut
287         zdh_s_sub(ji)  =  MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice )
288         ! Heat flux by sublimation [W.m-2], < 0
289         !      sublimate first snow that had fallen, then pre-existing snow
290         zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
291         hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * q_s_1d(ji,1)  &
292            &                              ) * a_i_1d(ji) * r1_rdtice
293         ! Mass flux by sublimation
294         wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice
295         ! new snow thickness
296         ht_s_1d(ji)    =  MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) )
297         ! update precipitations after sublimation and correct sublimation
298         zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
299         zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1)
300      END DO
301     
[4688]302      ! --- Update snow diags --- !
[825]303      DO ji = kideb, kiut
[5167]304         dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
305      END DO
[825]306
[4688]307      !-------------------------------------------
308      ! 3.3 Update temperature, energy
309      !-------------------------------------------
310      ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
311      zq_s(:) = 0._wp 
[825]312      DO jk = 1, nlay_s
313         DO ji = kideb,kiut
[5134]314            rswitch       = MAX(  0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 )  )
315            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *                          &
[5167]316              &            ( (   zdh_s_pre(ji)             ) * zqprec(ji) +  &
317              &              (   ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
[4872]318            zq_s(ji)     =  zq_s(ji) + q_s_1d(ji,jk)
[825]319         END DO
320      END DO
321
[4688]322      !--------------------------
323      ! 3.4 Surface ice ablation
324      !--------------------------
325      zdeltah(:,:) = 0._wp ! important
326      DO jk = 1, nlay_i
[5167]327         DO ji = kideb, kiut
[5202]328            ztmelts           = - tmut * s_i_1d(ji,jk) + rt0          ! Melting point of layer k [K]
329           
330            IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
[4688]331
[5202]332               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]       
333               zdE            = 0._wp                                 ! Specific enthalpy difference   (J/kg, <0)
334                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted)
335               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing     
336                                                                      ! this should normally not happen, but sometimes, heat diffusion leads to this
337               zfmdt          = - zdeltah(ji,jk) * rhoic              ! Mass flux x time step > 0
338                         
339               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt
340               
341               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0]
[4688]342
[5202]343               ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
344               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
345               
346               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
347               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
348               
349               ! Contribution to mass flux
350               wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
[4688]351
[5202]352            ELSE                                !!! Surface melting
353               
354               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]
355               zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0]
356               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0
357               
358               zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0]
359               
360               zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0]
361               
362               zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0]
363               
364               zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat
365               
366               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt
367               
368               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0]
369               
370               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0]
371               
372               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
373               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
374               
375               ! Contribution to heat flux [W.m-2], < 0
376               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
377               
378               ! Total heat flux used in this process [W.m-2], > 0 
379               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
380               
381               ! Contribution to mass flux
382               wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
383               
[5167]384            END IF
[4688]385            ! record which layers have disappeared (for bottom melting)
386            !    => icount=0 : no layer has vanished
387            !    => icount=5 : 5 layers have vanished
[5202]388            rswitch       = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) 
389            icount(ji,jk) = NINT( rswitch )
390            zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
[4688]391
392            ! update heat content (J.m-2) and layer thickness
[4872]393            qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
[4688]394            h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
[825]395         END DO
[921]396      END DO
[4688]397      ! update ice thickness
398      DO ji = kideb, kiut
[4872]399         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) )
[4688]400      END DO
[825]401
[921]402      !
403      !------------------------------------------------------------------------------!
404      ! 4) Basal growth / melt                                                       !
405      !------------------------------------------------------------------------------!
406      !
[4688]407      !------------------
408      ! 4.1 Basal growth
409      !------------------
410      ! Basal growth is driven by heat imbalance at the ice-ocean interface,
411      ! between the inner conductive flux  (fc_bo_i), from the open water heat flux
412      ! (fhld) and the turbulent ocean flux (fhtur).
413      ! fc_bo_i is positive downwards. fhtur and fhld are positive to the ice
[825]414
[4688]415      ! If salinity varies in time, an iterative procedure is required, because
416      ! the involved quantities are inter-dependent.
417      ! Basal growth (dh_i_bott) depends upon new ice specific enthalpy (zEi),
418      ! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bott)
419      ! -> need for an iterative procedure, which converges quickly
420
[5167]421      num_iter_max = 1
422      IF( nn_icesal == 2 ) num_iter_max = 5
[825]423
[4688]424      ! Iterative procedure
425      DO iter = 1, num_iter_max
426         DO ji = kideb, kiut
427            IF(  zf_tt(ji) < 0._wp  ) THEN
[825]428
[4688]429               ! New bottom ice salinity (Cox & Weeks, JGR88 )
430               !--- zswi1  if dh/dt < 2.0e-8
431               !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
432               !--- zswi2  if dh/dt > 3.6e-7
433               zgrr               = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) )
434               zswi2              = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
435               zswi12             = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
436               zswi1              = 1. - zswi2 * zswi12
437               zfracs             = MIN ( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
438                  &               + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 )
[825]439
[4688]440               ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
[825]441
[4688]442               s_i_new(ji)        = zswitch_sal * zfracs * sss_m(ii,ij)  &  ! New ice salinity
[4872]443                                  + ( 1. - zswitch_sal ) * sm_i_1d(ji) 
[4688]444               ! New ice growth
[5123]445               ztmelts            = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K)
[825]446
[4872]447               zt_i_new           = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
[4688]448               
449               zEi                = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)     
[5123]450                  &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   &
451                  &               + rcp  * ( ztmelts-rt0 )         
[4688]452
[4872]453               zEw                = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
[4688]454
455               zdE                = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
456
457               dh_i_bott(ji)      = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
458
[4872]459               q_i_1d(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0)
[4688]460               
[5134]461            ENDIF
462         END DO
463      END DO
[4688]464
465      ! Contribution to Energy and Salt Fluxes
[825]466      DO ji = kideb, kiut
[4688]467         IF(  zf_tt(ji) < 0._wp  ) THEN
468            ! New ice growth
469                                   
470            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0)
471
[5123]472            ztmelts        = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K)
[4688]473           
[4872]474            zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
[4688]475           
476            zEi            = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)     
[5123]477               &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   &
478               &               + rcp  * ( ztmelts-rt0 )         
[4688]479           
[4872]480            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
[4688]481           
482            zdE            = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
483           
484            ! Contribution to heat flux to the ocean [W.m-2], >0 
[4872]485            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
[4688]486
487            ! Total heat flux used in this process [W.m-2], <0 
[4872]488            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
[4688]489           
490            ! Contribution to salt flux, <0
[5167]491            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice
[4688]492
493            ! Contribution to mass flux, <0
[4872]494            wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice
[4688]495
496            ! update heat content (J.m-2) and layer thickness
[4872]497            qh_i_old(ji,nlay_i+1) = qh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * q_i_1d(ji,nlay_i+1)
[4688]498            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji)
[825]499         ENDIF
500      END DO
501
[4688]502      !----------------
503      ! 4.2 Basal melt
504      !----------------
505      zdeltah(:,:) = 0._wp ! important
[825]506      DO jk = nlay_i, 1, -1
507         DO ji = kideb, kiut
[5202]508            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting
[825]509
[5123]510               ztmelts = - tmut * s_i_1d(ji,jk) + rt0  ! Melting point of layer jk (K)
[825]511
[4872]512               IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
[825]513
[5123]514                  zEi               = - q_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0)
[4688]515                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0)
516                                                                    ! set up at 0 since no energy is needed to melt water...(it is already melted)
[5202]517                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing     
518                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this
[4161]519
[4688]520                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk)
[825]521
[5202]522                  zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0
[825]523
[4688]524                  ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
[4872]525                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
[825]526
[4872]527                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
[5167]528                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
[4688]529                                   
530                  ! Contribution to mass flux
[4872]531                  wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
[825]532
[4688]533                  ! update heat content (J.m-2) and layer thickness
[4872]534                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
[4688]535                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
[825]536
[4688]537               ELSE                               !!! Basal melting
[825]538
[5202]539                  zEi             = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)
540                  zEw             = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0)
541                  zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0)
[825]542
[5202]543                  zfmdt           = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0)
[825]544
[5202]545                  zdeltah(ji,jk)  = - zfmdt * r1_rhoic         ! Gross thickness change
[4688]546
[5202]547                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
[4688]548                 
[5202]549                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors
[4688]550
[5202]551                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt
[4688]552
[5202]553                  zfmdt           = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0
[4688]554
[5202]555                  zQm             = zfmdt * zEw         ! Heat exchanged with ocean
[4688]556
557                  ! Contribution to heat flux to the ocean [W.m-2], <0 
[5202]558                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
[4688]559
[4872]560                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
[5202]561                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
[4688]562                 
563                  ! Total heat flux used in this process [W.m-2], >0 
[5202]564                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
[4688]565                 
566                  ! Contribution to mass flux
[5202]567                  wfx_bom_1d(ji)  =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
[4688]568
569                  ! update heat content (J.m-2) and layer thickness
[4872]570                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
[4688]571                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
572               ENDIF
573           
574            ENDIF
[5123]575         END DO
576      END DO
[4688]577
[825]578      !-------------------------------------------
[4688]579      ! Update temperature, energy
[825]580      !-------------------------------------------
581      DO ji = kideb, kiut
[4872]582         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_bott(ji) )
[4688]583      END DO 
[825]584
[4688]585      !-------------------------------------------
586      ! 5. What to do with remaining energy
587      !-------------------------------------------
588      ! If heat still available for melting and snow remains, then melt more snow
589      !-------------------------------------------
590      zdeltah(:,:) = 0._wp ! important
591      DO ji = kideb, kiut
592         zq_rema(ji)     = zq_su(ji) + zq_bo(ji) 
[5128]593         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )   ! =1 if snow
[5134]594         rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,1) - epsi20 ) )
595         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 )
[5128]596         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting
597         dh_s_tot (ji)   = dh_s_tot(ji)  + zdeltah(ji,1)
598         ht_s_1d   (ji)  = ht_s_1d(ji)   + zdeltah(ji,1)
599       
[5134]600         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * q_s_1d(ji,1)                ! update available heat (J.m-2)
[5128]601         ! heat used to melt snow
[5134]602         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * q_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0)
[5128]603         ! Contribution to mass flux
604         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
605         !   
[4688]606         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
607         ! Remaining heat flux (W.m-2) is sent to the ocean heat budget
[5123]608         hfx_out(ii,ij)  = hfx_out(ii,ij) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
[825]609
[5128]610         IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
[4688]611      END DO
612     
[921]613      !
614      !------------------------------------------------------------------------------|
615      !  6) Snow-Ice formation                                                       |
616      !------------------------------------------------------------------------------|
[1572]617      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
618      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
[825]619      DO ji = kideb, kiut
[1572]620         !
[4872]621         dh_snowice(ji) = MAX(  0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic )  )
[825]622
[5167]623         ht_i_1d(ji)    = ht_i_1d(ji) + dh_snowice(ji)
624         ht_s_1d(ji)    = ht_s_1d(ji) - dh_snowice(ji)
[825]625
[4688]626         ! Salinity of snow ice
627         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
[5123]628         zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) * r1_rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji)
[825]629
630         ! entrapment during snow ice formation
[5134]631         ! new salinity difference stored (to be used in limthd_sal.F90)
[5123]632         IF (  nn_icesal == 2  ) THEN
[5134]633            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )
[4161]634            ! salinity dif due to snow-ice formation
[5134]635            dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch     
[4161]636            ! salinity dif due to bottom growth
[4688]637            IF (  zf_tt(ji)  < 0._wp ) THEN
[5134]638               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch
[4161]639            ENDIF
640         ENDIF
[825]641
[4688]642         ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
643         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
644         zfmdt          = ( rhosn - rhoic ) * MAX( dh_snowice(ji), 0._wp )    ! <0
645         zsstK          = sst_m(ii,ij) + rt0                               
646         zEw            = rcp * ( zsstK - rt0 )
647         zQm            = zfmdt * zEw 
648         
649         ! Contribution to heat flux
[4872]650         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
[825]651
[4688]652         ! Contribution to salt flux
[4872]653         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_1d(ji) * zfmdt * r1_rdtice 
[4688]654         
655         ! Contribution to mass flux
656         ! All snow is thrown in the ocean, and seawater is taken to replace the volume
[4872]657         wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice
658         wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice
[4688]659
660         ! update heat content (J.m-2) and layer thickness
[4872]661         qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_1d(ji,1) + zfmdt * zEw
[4688]662         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
663         
[5134]664      END DO
[825]665
[2715]666      !
[4688]667      !-------------------------------------------
668      ! Update temperature, energy
669      !-------------------------------------------
670      DO ji = kideb, kiut
[4990]671         rswitch     =  1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 
[5123]672         t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0
[5134]673      END DO
[4688]674
675      DO jk = 1, nlay_s
676         DO ji = kideb,kiut
677            ! mask enthalpy
[5134]678            rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) )  )
679            q_s_1d(ji,jk) = rswitch * q_s_1d(ji,jk)
[4872]680            ! recalculate t_s_1d from q_s_1d
[5134]681            t_s_1d(ji,jk) = rt0 + rswitch * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic )
[4688]682         END DO
683      END DO
[5134]684
685      ! --- ensure that a_i = 0 where ht_i = 0 ---
686      WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp
[5123]687     
[5487]688      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw )
689      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
[5202]690      CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i )
691      CALL wrk_dealloc( jpij, nlay_i, icount )
[2715]692      !
[4161]693      !
[921]694   END SUBROUTINE lim_thd_dh
[5407]695
696
697   !!--------------------------------------------------------------------------
698   !! INTERFACE lim_thd_snwblow
699   !! ** Purpose :   Compute distribution of precip over the ice
700   !!--------------------------------------------------------------------------
701   SUBROUTINE lim_thd_snwblow_2d( pin, pout )
[5487]702      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( pfrld or (1. - a_i_b) )
703      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
[5407]704      pout = ( 1._wp - ( pin )**rn_betas )
705   END SUBROUTINE lim_thd_snwblow_2d
706
707   SUBROUTINE lim_thd_snwblow_1d( pin, pout )
[5487]708      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin
709      REAL(wp), DIMENSION(:), INTENT(inout) :: pout
[5407]710      pout = ( 1._wp - ( pin )**rn_betas )
711   END SUBROUTINE lim_thd_snwblow_1d
712
[1572]713   
[825]714#else
[1572]715   !!----------------------------------------------------------------------
716   !!   Default option                               NO  LIM3 sea-ice model
717   !!----------------------------------------------------------------------
[825]718CONTAINS
719   SUBROUTINE lim_thd_dh          ! Empty routine
720   END SUBROUTINE lim_thd_dh
721#endif
[1572]722
723   !!======================================================================
[921]724END MODULE limthd_dh
Note: See TracBrowser for help on using the repository browser.