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.
icethd_pnd.F90 in NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd_pnd.F90 @ 12500

Last change on this file since 12500 was 12500, checked in by dancopsey, 4 years ago

Add namelist variables to control new functionality.
File size: 21.3 KB
Line 
1MODULE icethd_pnd 
2   !!======================================================================
3   !!                     ***  MODULE  icethd_pnd   ***
4   !!   sea-ice: Melt ponds on top of sea ice
5   !!======================================================================
6   !! history :       !  2012     (O. Lecomte)       Adaptation from Flocco and Turner
7   !!                 !  2017     (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation
8   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
9   !!----------------------------------------------------------------------
10#if defined key_si3
11   !!----------------------------------------------------------------------
12   !!   'key_si3' :                                     SI3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   ice_thd_pnd_init : some initialization and namelist read
15   !!   ice_thd_pnd      : main calling routine
16   !!----------------------------------------------------------------------
17   USE phycst         ! physical constants
18   USE dom_oce        ! ocean space and time domain
19   USE ice            ! sea-ice: variables
20   USE ice1D          ! sea-ice: thermodynamics variables
21   USE icetab         ! sea-ice: 1D <==> 2D transformation
22   !
23   USE in_out_manager ! I/O manager
24   USE lib_mpp        ! MPP library
25   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
26   USE timing         ! Timing
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   ice_thd_pnd_init    ! routine called by icestp.F90
32   PUBLIC   ice_thd_pnd         ! routine called by icestp.F90
33
34   INTEGER ::              nice_pnd    ! choice of the type of pond scheme
35   !                                   ! associated indices:
36   INTEGER, PARAMETER ::   np_pndNO  = 0   ! No pond scheme
37   INTEGER, PARAMETER ::   np_pndCST = 1   ! Constant pond scheme
38   INTEGER, PARAMETER ::   np_pndH12 = 2   ! Evolutive pond scheme (Holland et al. 2012)
39
40   REAL(wp), PUBLIC, PARAMETER :: pnd_lid_max = 0.015_wp !  pond lid thickness above which the ponds disappear from the albedo calculation
41   REAL(wp), PUBLIC, PARAMETER :: pnd_lid_min = 0.005_wp !  pond lid thickness below which the full pond area is used in the albedo calculation
42
43   REAL(wp), PARAMETER :: viscosity_dyn = 1.79e-3_wp
44
45   !! * Substitutions
46#  include "vectopt_loop_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
49   !! $Id$
50   !! Software governed by the CeCILL license (see ./LICENSE)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE ice_thd_pnd
55      !!-------------------------------------------------------------------
56      !!               ***  ROUTINE ice_thd_pnd   ***
57      !!               
58      !! ** Purpose :   change melt pond fraction
59      !!               
60      !! ** Method  :   brut force
61      !!-------------------------------------------------------------------
62      !
63      SELECT CASE ( nice_pnd )
64      !
65      CASE (np_pndCST)   ;   CALL pnd_CST    !==  Constant melt ponds  ==!
66         !
67      CASE (np_pndH12)   ;   CALL pnd_H12    !==  Holland et al 2012 melt ponds  ==!
68         !
69      END SELECT
70      !
71   END SUBROUTINE ice_thd_pnd 
72
73
74   SUBROUTINE pnd_CST 
75      !!-------------------------------------------------------------------
76      !!                ***  ROUTINE pnd_CST  ***
77      !!
78      !! ** Purpose :   Compute melt pond evolution
79      !!
80      !! ** Method  :   Melt pond fraction and thickness are prescribed
81      !!                to non-zero values when t_su = 0C
82      !!
83      !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd)
84      !!               
85      !! ** Note   : Coupling with such melt ponds is only radiative
86      !!             Advection, ridging, rafting... are bypassed
87      !!
88      !! ** References : Bush, G.W., and Trump, D.J. (2017)
89      !!-------------------------------------------------------------------
90      INTEGER  ::   ji        ! loop indices
91      !!-------------------------------------------------------------------
92      DO ji = 1, npti
93         !
94         IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN
95            a_ip_frac_1d(ji) = rn_apnd
96            h_ip_1d(ji)      = rn_hpnd   
97            a_ip_1d(ji)      = a_ip_frac_1d(ji) * a_i_1d(ji)
98         ELSE
99            a_ip_frac_1d(ji) = 0._wp
100            h_ip_1d(ji)      = 0._wp   
101            a_ip_1d(ji)      = 0._wp
102         ENDIF
103         !
104      END DO
105      !
106   END SUBROUTINE pnd_CST
107
108
109   SUBROUTINE pnd_H12
110      !!-------------------------------------------------------------------
111      !!                ***  ROUTINE pnd_H12  ***
112      !!
113      !! ** Purpose    : Compute melt pond evolution
114      !!
115      !! ** Method     : Empirical method. A fraction of meltwater is accumulated in ponds
116      !!                 and sent to ocean when surface is freezing
117      !!
118      !!                 pond growth:      Vp = Vp + dVmelt
119      !!                    with dVmelt = R/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i
120      !!                 pond contraction: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)
121      !!                    with Tp = -2degC
122      !! 
123      !! ** Tunable parameters : (no real expertise yet, ideas?)
124      !!
125      !! ** Note       : Stolen from CICE for quick test of the melt pond
126      !!                 radiation and freshwater interfaces
127      !!                 Coupling can be radiative AND freshwater
128      !!                 Advection, ridging, rafting are called
129      !!
130      !! ** References : Holland, M. M. et al (J Clim 2012)
131      !!-------------------------------------------------------------------
132      REAL(wp), PARAMETER ::   zrmin       = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding
133      REAL(wp), PARAMETER ::   zrmax       = 0.70_wp  ! maximum     -           -         -         -            -
134      REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp   ! pond aspect ratio
135      REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature
136      REAL(wp), PARAMETER ::   max_h_diff_s = -1.0E-6  ! Maximum meltpond depth change due to leaking or overflow (m s-1)
137      !
138      REAL(wp) ::   tot_mlt          ! Total ice and snow surface melt (some goes into ponds, some into the ocean)
139      REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding
140      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding (equivalent volume change)
141      REAL(wp) ::   z1_Tp            ! inverse reference temperature
142      REAL(wp) ::   z1_rhow          ! inverse freshwater density
143      REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio
144      REAL(wp) ::   z1_rhoi          ! inverse ice density
145      REAL(wp) ::   zfac, zdum
146      REAL(wp) ::   t_grad           ! Temperature deficit for refreezing
147      REAL(wp) ::   omega_dt         ! Time independent accumulated variables used for freezing
148      REAL(wp) ::   lh_ip_end        ! Lid thickness at end of timestep (temporary variable)
149      REAL(wp) ::   zdh_frz          ! Amount of melt pond that freezes (m)
150      REAL(wp) ::   v_ip_old         ! Pond volume before leaking back to the ocean
151      REAL(wp) ::   dh_ip_over       ! Pond thickness change due to overflow or leaking
152      REAL(wp) ::   dh_i_pndleak     ! Grid box mean change in water depth due to leaking back to the ocean
153      REAL(wp) ::   h_gravity_head   ! Height above sea level of the top of the melt pond
154      REAL(wp) ::   h_percolation    ! Distance between the base of the melt pond and the base of the sea ice
155      REAL(wp) ::   Sbr              ! Brine salinity
156      REAL(wp), DIMENSION(nlay_i) ::   phi     ! liquid fraction
157      REAL(wp) ::   perm             ! Permeability of the sea ice
158      REAL(wp) ::   za_ip            ! Temporary pond fraction
159      REAL(wp) ::   max_h_diff_ts    ! Maximum meltpond depth change due to leaking or overflow (m per ts)
160     
161      !
162      INTEGER  ::   ji, jk   ! loop indices
163      !!-------------------------------------------------------------------
164      z1_rhow        = 1._wp / rhow 
165      z1_zpnd_aspect = 1._wp / zpnd_aspect
166      z1_Tp          = 1._wp / zTp 
167      z1_rhoi        = 1._wp / rhoi
168      max_h_diff_ts  = max_h_diff_s * rdt_ice
169
170      ! Define time-independent field for use in refreezing
171      omega_dt = 2.0_wp * rcnd_i * rdt_ice / (rLfus * rhow)
172
173      DO ji = 1, npti
174
175         !                                                            !----------------------------------------------------!
176         IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction !
177            !                                                         !----------------------------------------------------!
178            !--- Remove ponds on thin ice or tiny ice fractions
179            a_ip_1d(ji)      = 0._wp
180            a_ip_frac_1d(ji) = 0._wp
181            h_ip_1d(ji)      = 0._wp
182            lh_ip_1d(ji)     = 0._wp
183            !                                                     !--------------------------------!
184         ELSE                                                     ! Case ice thickness >= rn_himin !
185            !                                                     !--------------------------------!
186            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step
187
188            ! To avoid divide by zero errors in some calculations we will use a temporary pond fraction variable that has a minimum value of epsi06
189            za_ip = a_ip_1d(ji)
190            IF ( za_ip < epsi06 ) za_ip = epsi06
191            !
192            ! available meltwater for melt ponding [m, >0] and fraction
193            tot_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow 
194            IF ( ln_pnd_totfrac ) THEN
195               zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * at_i_1d(ji) ! Use total ice fraction
196            ELSE
197               zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * a_i_1d(ji)  ! Use category ice fraction
198            ENDIF
199            zdv_mlt = zfr_mlt * tot_mlt 
200            !
201            !--- Pond gowth ---!
202            v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt
203            !
204            !--- Lid shrinking. ---!
205            IF ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip
206            !
207            ! melt pond mass flux (<0)
208            IF( ln_use_pndmass .AND. zdv_mlt > 0._wp ) THEN
209               zfac = zdv_mlt * rhow * r1_rdtice
210               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
211               !
212               ! adjust ice/snow melting flux to balance melt pond flux (>0)
213               zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )
214               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)
215               wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum)
216            ENDIF
217            !
218            !--- Pond contraction (due to refreezing) ---!
219            IF ( ln_pnd_lids ) THEN
220
221               ! Code to use if using melt pond lids
222               IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN
223                  t_grad = (zTp+rt0) - t_su_1d(ji)
224
225                  ! The following equation is a rearranged form of:
226                  ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow)
227                  ! where: lid_thickness_start = lh_ip_1d(ji)
228                  !        lid_thickness_end = lh_ip_end
229                  !        omega_dt is a bunch of terms in the equation that do not change
230                  ! Note the use of rhow instead of rhoi as we are working with volumes and it is mathematically easier
231                  ! if the water and ice specific volumes (for the lid and the pond) are the same (have the same density).
232
233                  lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2)
234                  zdh_frz = lh_ip_end - lh_ip_1d(ji)
235
236                  ! Pond shrinking
237                  v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji)
238
239                  ! Lid growing
240                  F ( ln_pnd_lids )  lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz
241               ELSE
242                  zdh_frz = 0._wp
243               END IF
244
245            ELSE
246               ! Code to use if not using melt pond lids
247               v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp )
248            ENDIF
249
250            !
251            ! Make sure pond volume or lid thickness has not gone negative
252            IF ( v_ip_1d(ji) < 0._wp ) v_ip_1d(ji) = 0._wp 
253            IF ( lh_ip_1d(ji) < 0._wp ) lh_ip_1d(ji) = 0._wp
254            !
255            ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac
256            !    h_ip = zpnd_aspect * a_ip_frac = zpnd_aspect * a_ip/a_i
257            a_ip_1d(ji)      = SQRT( v_ip_1d(ji) * z1_zpnd_aspect * a_i_1d(ji) )
258            a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
259            h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji)
260           
261            !--- Pond overflow and vertical flushing ---!
262            IF ( ln_pnd_overflow ) THEN
263
264               ! If pond area exceeds a_pnd_avail_1d(ji) * a_i_1d(ji) then reduce the pond volume
265               IF ( a_ip_1d(ji) > a_pnd_avail_1d(ji) * a_i_1d(ji) ) THEN
266                   v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
267                   dh_ip_over = zpnd_aspect * a_pnd_avail_1d(ji) - h_ip_1d(ji)   ! This will be a negative number
268                   dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit
269                   h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over)
270                   a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect
271                   a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji)
272                   v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
273               ENDIF
274
275               ! If pond depth exceeds half the ice thickness then reduce the pond volume
276               IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN
277                   v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
278                   dh_ip_over = 0.5_wp * h_i_1d(ji) - h_ip_1d(ji)                ! This will be a negative number
279                   dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit
280                   h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over)
281                   a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect
282                   a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji)
283                   v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
284               ENDIF
285
286               !-- Vertical flushing of pond water --!
287               ! The height above sea level of the top of the melt pond is the ratios of density of ice and water times the ice depth.
288               ! This assumes the pond is sitting on top of the ice.
289               h_gravity_head = h_i_1d(ji) * (rhow - rhoi) * z1_rhow
290
291               ! The depth through which water percolates is the distance under the melt pond
292               h_percolation = h_i_1d(ji) - h_ip_1d(ji)
293
294               ! Calculate the permeability of the ice (Assur 1958)
295               DO jk = 1, nlay_i
296                   Sbr = - 1.2_wp                         &
297                         - 21.8_wp     * (t_i_1d(ji,jk)-rt0)    &
298                         - 0.919_wp    * (t_i_1d(ji,jk)-rt0)**2 &
299                         - 0.01878_wp  * (t_i_1d(ji,jk)-rt0)**3
300                   phi(jk) = sz_i_1d(ji,jk)/Sbr
301               END DO
302               perm = 3.0e-08_wp * (minval(phi))**3
303
304               ! Do the drainage using Darcy's law
305               IF ( perm > 0._wp ) THEN
306                   dh_ip_over = -1.0_wp*perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation)  ! This should be a negative number
307                   dh_ip_over = MIN(dh_ip_over, 0._wp)      ! Make sure it is negative
308
309                   v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
310                   dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit
311                   h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over)
312                   a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect
313                   a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji)
314                   v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
315               ENDIF
316            ENDIF
317
318            ! If lid thickness is ten times greater than pond thickness then remove pond
319            IF ( ln_pnd_lids ) THEN           
320               IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN
321                   a_ip_1d(ji)      = 0._wp
322                   a_ip_frac_1d(ji) = 0._wp
323                   h_ip_1d(ji)      = 0._wp
324                   lh_ip_1d(ji)     = 0._wp
325                   v_ip_1d(ji)      = 0._wp
326               ENDIF
327            ENDIF
328
329            ! If any of the previous changes has removed all the ice thickness then remove ice area.
330            IF ( h_i_1d(ji) == 0._wp ) THEN
331                a_i_1d(ji) = 0._wp
332                h_s_1d(ji) = 0._wp
333            ENDIF
334
335            !
336         ENDIF
337
338      END DO
339      !
340   END SUBROUTINE pnd_H12
341
342
343   SUBROUTINE ice_thd_pnd_init 
344      !!-------------------------------------------------------------------
345      !!                  ***  ROUTINE ice_thd_pnd_init   ***
346      !!
347      !! ** Purpose : Physical constants and parameters linked to melt ponds
348      !!              over sea ice
349      !!
350      !! ** Method  :  Read the namthd_pnd  namelist and check the melt pond 
351      !!               parameter values called at the first timestep (nit000)
352      !!
353      !! ** input   :   Namelist namthd_pnd 
354      !!-------------------------------------------------------------------
355      INTEGER  ::   ios, ioptio   ! Local integer
356      !!
357      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb,          &
358                            rn_pnd_min, rn_pnd_max, ln_pnd_overflow, ln_pnd_lids, ln_pnd_totfrac,  &
359                            ln_use_pndmass
360      !!-------------------------------------------------------------------
361      !
362      REWIND( numnam_ice_ref )              ! Namelist namthd_pnd  in reference namelist : Melt Ponds 
363      READ  ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901)
364901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namthd_pnd  in reference namelist' )
365      REWIND( numnam_ice_cfg )              ! Namelist namthd_pnd  in configuration namelist : Melt Ponds
366      READ  ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 )
367902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist' )
368      IF(lwm) WRITE ( numoni, namthd_pnd )
369      !
370      IF(lwp) THEN                        ! control print
371         WRITE(numout,*)
372         WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
373         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
374         WRITE(numout,*) '   Namelist namicethd_pnd:'
375         WRITE(numout,*) '      Melt ponds activated or not                                     ln_pnd     = ', ln_pnd
376         WRITE(numout,*) '         Evolutive  melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12
377         WRITE(numout,*) '            Minimum ice fraction that contributes to melt ponds       rn_pnd_min = ', rn_pnd_min
378         WRITE(numout,*) '            Maximum ice fraction that contributes to melt ponds       rn_pnd_max = ', rn_pnd_max
379         WRITE(numout,*) '            Use total ice fraction instead of category ice fraction   ln_pnd_totfrac  = ',ln_pnd_totfrac
380         WRITE(numout,*) '            Allow ponds to overflow and have vertical flushing        ln_pnd_overflow = ',ln_pnd_overflow
381         WRITE(numout,*) '            Melt ponds can have frozen lids                           ln_pnd_lids     = ',ln_pnd_lids
382         WRITE(numout,*) '            Use melt pond mass flux diagnostic, passing to ocean      ln_use_pndmass  = ',ln_use_pndmass
383         WRITE(numout,*) '         Prescribed melt pond fraction and depth                      ln_pnd_CST = ', ln_pnd_CST
384         WRITE(numout,*) '            Prescribed pond fraction                                  rn_apnd    = ', rn_apnd
385         WRITE(numout,*) '            Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd
386         WRITE(numout,*) '         Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb
387      ENDIF
388      !
389      !                             !== set the choice of ice pond scheme ==!
390      ioptio = 0
391      IF( .NOT.ln_pnd ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndNO     ;   ENDIF
392      IF( ln_pnd_CST  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
393      IF( ln_pnd_H12  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF
394      IF( ioptio /= 1 )   &
395         & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' )
396      !
397      SELECT CASE( nice_pnd )
398      CASE( np_pndNO )         
399         IF( ln_pnd_alb ) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF
400      END SELECT
401      !
402   END SUBROUTINE ice_thd_pnd_init
403   
404#else
405   !!----------------------------------------------------------------------
406   !!   Default option          Empty module          NO SI3 sea-ice model
407   !!----------------------------------------------------------------------
408#endif 
409
410   !!======================================================================
411END MODULE icethd_pnd 
Note: See TracBrowser for help on using the repository browser.