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_add_pond_lids_prints/src/ICE – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_pnd.F90 @ 12449

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