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 @ 12566

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

Move calculation of effective pond fraction to icethd_pnd.F90 (not sbccpl)

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