source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90 @ 8554

Last change on this file since 8554 was 8534, checked in by clem, 3 years ago

changes in style - part6 - pure cosmetics

File size: 20.6 KB
Line 
1MODULE icealb
2   !!======================================================================
3   !!                       ***  MODULE  icealb  ***
4   !! Atmospheric forcing:  Albedo over sea ice
5   !!=====================================================================
6   !! History :  4.0  ! 2017-07  (C. Rousset) Split ice and ocean albedos
7   !!----------------------------------------------------------------------
8#if defined key_lim3
9   !!----------------------------------------------------------------------
10   !!   'key_lim3'                                       ESIM sea-ice model
11   !!----------------------------------------------------------------------
12   !!   ice_alb        : albedo for ice (clear and overcast skies)
13   !!   ice_alb_init   : initialisation of albedo computation
14   !!----------------------------------------------------------------------
15   USE ice, ONLY: jpl ! sea-ice: number of categories
16   USE phycst         ! physical constants
17   !
18   USE in_out_manager ! I/O manager
19   USE lib_mpp        ! MPP library
20   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
21   USE timing         ! Timing
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   ice_alb_init   ! called in icestp
27   PUBLIC   ice_alb        ! called in iceforcing.F90 and iceupdate.F90
28
29   REAL(wp), PUBLIC, PARAMETER ::   rn_alb_oce = 0.066   !: ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
30
31   REAL(wp) , PARAMETER ::   rc1    = 0.05    ! snow thickness (only for nn_ice_alb=0)
32   REAL(wp) , PARAMETER ::   rc2    = 0.10    !  "        "
33   REAL(wp) , PARAMETER ::   rcloud = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0)
34   REAL(wp) , PARAMETER ::   r1_c1 = 1. / rc1
35   REAL(wp) , PARAMETER ::   r1_c2 = 1. / rc2
36   !
37   ! ** albedo namelist (namalb)
38   INTEGER  ::   nn_ice_alb       ! type of albedo scheme: 0: Shine & Henderson-Sellers (JGR 1985)
39   !                                      !                         1: "home made" based on Brandt et al. (JClim 2005)
40   !                                      !                             and Grenfell & Perovich (JGR 2004)
41   !                                      !                         2: Same as 1 but with melt ponds
42   REAL(wp) ::   rn_alb_sdry      ! dry snow albedo
43   REAL(wp) ::   rn_alb_smlt      ! melting snow albedo
44   REAL(wp) ::   rn_alb_idry      ! dry ice albedo
45   REAL(wp) ::   rn_alb_imlt      ! bare puddled ice albedo
46   REAL(wp) ::   rn_alb_dpnd      ! ponded ice albedo
47
48   !!----------------------------------------------------------------------
49   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
50   !! $Id: icealb.F90 8268 2017-07-03 15:01:04Z clem $
51   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE ice_alb( pt_ice, ph_ice, ph_snw, pafrac_pnd, ph_pnd, ld_pnd, pa_ice_cs, pa_ice_os )
56      !!----------------------------------------------------------------------
57      !!               ***  ROUTINE ice_alb  ***
58      !!         
59      !! ** Purpose :   Computation of the albedo of the snow/ice system
60      !!       
61      !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb)
62      !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies
63      !!                  1: the scheme is "home made" (for cloudy skies) and based on Brandt et al. (J. Climate 2005)
64      !!                                                                           and Grenfell & Perovich (JGR 2004)
65      !!                  2: fractional surface-based formulation of scheme 1 (NEW)
66      !!                Description of scheme 1:
67      !!                  1) Albedo dependency on ice thickness follows the findings from Brandt et al (2005)
68      !!                     which are an update of Allison et al. (JGR 1993) ; Brandt et al. 1999
69      !!                     0-5cm  : linear function of ice thickness
70      !!                     5-150cm: log    function of ice thickness
71      !!                     > 150cm: constant
72      !!                  2) Albedo dependency on snow thickness follows the findings from Grenfell & Perovich (2004)
73      !!                     i.e. it increases as -EXP(-snw_thick/0.02) during freezing and -EXP(-snw_thick/0.03) during melting
74      !!                  3) Albedo dependency on clouds is speculated from measurements of Grenfell and Perovich (2004)
75      !!                     i.e. cloudy-clear albedo depend on cloudy albedo following a 2d order polynomial law
76      !!                  4) The needed 4 parameters are: dry and melting snow, freezing ice and bare puddled ice
77      !!
78      !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions:
79      !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo
80      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger
81      !!                     under melting conditions than under freezing conditions
82      !!                  3) the evolution of ice albedo as a function of ice thickness shows 
83      !!                     3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic
84      !!
85      !! References :   Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250.
86      !!                Brandt et al. 2005, J. Climate, vol 18
87      !!                Grenfell & Perovich 2004, JGR, vol 109
88      !!----------------------------------------------------------------------
89      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice       !  ice surface temperature (Kelvin)
90      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_ice       !  sea-ice thickness
91      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_snw       !  snow depth
92      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pafrac_pnd   !  melt pond relative fraction (per unit ice area)
93      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_pnd       !  melt pond depth
94      LOGICAL , INTENT(in   )                   ::   ld_pnd       !  melt ponds radiatively active or not
95      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_cs    !  albedo of ice under clear    sky
96      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_os    !  albedo of ice under overcast sky
97      !
98      INTEGER  ::   ji, jj, jl                ! dummy loop indices
99      REAL(wp) ::   zswitch, z1_c1, z1_c2     ! local scalar
100      REAL(wp) ::   z1_href_pnd               !   -      -                 
101      REAL(wp) ::   zalb_sm, zalb_sf, zalb_st ! albedo of snow melting, freezing, total
102      !
103      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb, zalb_it   ! intermediate variable & albedo of ice (snow free)
104!! MV MP
105      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_pnd        ! ponded sea ice albedo
106      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_ice        ! bare sea ice albedo
107      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_snw        ! snow-covered sea ice albedo
108      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_snw      ! relative snow fraction
109      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_ice      ! relative ice fraction
110      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_pnd      ! relative ice fraction (effective)
111      !!---------------------------------------------------------------------
112      !
113      IF( nn_timing == 1 )   CALL timing_start('icealb')
114      !
115      !-----------------------------------------------------
116      !  Snow-free albedo (no ice thickness dependence yet)
117      !-----------------------------------------------------
118      !
119      ! Part common to nn_ice_alb = 0, 1, 2
120      !
121      IF ( .NOT. ld_pnd ) THEN   !--- No melt ponds OR radiatively inactive melt ponds
122         ! Bare ice albedo is prescribed, with implicit assumption on pond fraction and depth
123         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0 )   ;   zalb(:,:,:) = rn_alb_imlt
124         ELSEWHERE                                           ;   zalb(:,:,:) = rn_alb_idry
125         END WHERE
126      ENDIF
127
128      SELECT CASE ( nn_ice_alb )       ! select a parameterization
129      !
130      !              !------------------------------------------
131      CASE( 0 )      !  Shine and Henderson-Sellers (1985)     !   (based on clear sky values)
132         !           !------------------------------------------
133         !
134         !                       ! Thickness-dependent bare ice albedo
135         WHERE    ( 1.5  < ph_ice                     )   ;   zalb_it = zalb
136         ELSEWHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )   ;   zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )
137         ELSEWHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )   ;   zalb_it = 0.2467 + 0.7049 * ph_ice              &
138            &                                                                  - 0.8608 * ph_ice * ph_ice     &
139            &                                                                  + 0.3812 * ph_ice * ph_ice * ph_ice
140         ELSEWHERE                                        ;   zalb_it = 0.1    + 3.6    * ph_ice
141         END WHERE
142         !
143         IF ( ld_pnd ) THEN      ! Depth-dependent ponded ice albedo
144            z1_href_pnd = 1. / 0.05    ! inverse of the characteristic length scale (Lecomte et al. 2015)
145            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd * z1_href_pnd ) 
146            !                          ! Snow-free ice albedo is a function of pond fraction
147            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0 )   
148               zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd 
149            END WHERE
150         ENDIF 
151         !
152         DO jl = 1, jpl
153            DO jj = 1, jpj
154               DO ji = 1, jpi
155                  !                    ! Freezing snow
156                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > rc2
157                  zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - rc1 ) ) )
158                  zalb_sf = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
159                     &                                   + ph_snw(ji,jj,jl) * ( rn_alb_sdry - zalb_it(ji,jj,jl) ) * r1_c1  )   &
160                     &    +           zswitch   * rn_alb_sdry 
161                     !
162                  !                    ! Melting snow
163                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > rc2
164                  zswitch = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - rc2 ) )
165                  zalb_sm = ( 1._wp - zswitch ) * ( rn_alb_imlt + ph_snw(ji,jj,jl) * ( rn_alb_smlt - rn_alb_imlt ) * r1_c2 )   &
166                     &      +         zswitch   *   rn_alb_smlt 
167                     !
168                  !                    ! Snow albedo
169                  zswitch =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
170                  zalb_st =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
171                  !
172                  !                    ! Surface albedo
173                  zswitch             = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
174                  pa_ice_cs(ji,jj,jl) = zswitch * zalb_st + ( 1._wp - zswitch ) * zalb_it(ji,jj,jl)
175                  !
176               END DO
177            END DO
178         END DO
179         !
180         pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rcloud       ! Oberhuber correction for overcast sky
181         !
182         !
183         !           !------------------------------------------
184      CASE( 1 )      !  New parameterization (2016)            !    (based on overcast sky values)
185         !           !------------------------------------------
186         !
187         !                    compilation of values from literature (reference overcast sky values)
188         !                          rn_alb_sdry = 0.85      ! dry snow
189         !                          rn_alb_smlt = 0.75      ! melting snow
190         !                          rn_alb_idry = 0.60      ! bare frozen ice
191         !                          rn_alb_dpnd = 0.36      ! ponded-ice overcast albedo (Lecomte et al, 2015)
192         !                                                  ! early melt pnds 0.27, late melt ponds 0.14 Grenfell & Perovich
193         !                    Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved
194         !                          rn_alb_sdry = 0.85      ! dry snow
195         !                          rn_alb_smlt = 0.72      ! melting snow
196         !                          rn_alb_idry = 0.65      ! bare frozen ice
197         !                    Brandt et al 2005 (East Antarctica)
198         !                          rn_alb_sdry = 0.87      ! dry snow
199         !                          rn_alb_smlt = 0.82      ! melting snow
200         !                          rn_alb_idry = 0.54      ! bare frozen ice
201         !
202         !              !--- Computation of snow-free ice albedo
203         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
204         z1_c2 = 1. / 0.05
205
206         !              !--- Accounting for the ice-thickness dependency
207         WHERE     ( 1.5  < ph_ice                     )   ;   zalb_it = zalb
208         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )   ;   zalb_it = zalb + ( 0.18 - zalb ) * z1_c1 * ( LOG(1.5) - LOG(ph_ice) )
209         ELSE WHERE                                        ;   zalb_it = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice
210         END WHERE
211         !
212         IF ( ld_pnd ) THEN      ! Depth-dependent ponded ice albedo
213            z1_href_pnd = 0.05        ! inverse of the characteristic length scale (Lecomte et al. 2015)
214            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd * z1_href_pnd ) 
215            !
216            !                    ! Snow-free ice albedo is weighted mean of ponded ice and bare ice contributions
217            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0 )
218               zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd
219            END WHERE
220         ENDIF 
221         !
222         !              !--- Overcast sky surface albedo (accounting for snow, ice melt ponds)
223         z1_c1 = 1. / 0.02
224         z1_c2 = 1. / 0.03
225         DO jl = 1, jpl
226            DO jj = 1, jpj
227               DO ji = 1, jpi
228                  ! Snow depth dependence of snow albedo
229                  zalb_sf = rn_alb_sdry - ( rn_alb_sdry - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c1 );
230                  zalb_sm = rn_alb_smlt - ( rn_alb_smlt - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c2 );
231
232                  ! Snow albedo (MV I guess we could use rt0 instead of rt0_snow)
233                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
234                  zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
235
236                  ! Surface albedo   
237                  zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
238                  pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl)
239
240              END DO
241            END DO
242         END DO
243         !
244         !              !--- Clear sky surface albedo
245         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ); 
246         !
247         !
248         !           !------------------------------------------
249      CASE( 2 )      !  Fractional surface-based parameterization (2017)
250         !           !------------------------------------------
251         !              MV: I propose this formulation that is more elegant, and more easy to expand towards
252         !              varying pond and snow fraction.
253         !              Formulation 1 has issues to handle ponds and snow together that can't easily be fixed.
254         !              This one handles it better, I believe.
255         !
256         !                 !==  Snow, bare ice and ponded ice fractions  ==!
257         !
258         !                       ! Specific fractions (zafrac) refer to relative area covered by snow within each ice category
259         !
260         !                       !--- Effective pond fraction (for now, we prevent melt ponds and snow at the same time)
261         zafrac_pnd = 0._wp
262         IF ( ld_pnd ) THEN
263            WHERE( ph_snw == 0._wp )   zafrac_pnd = pafrac_pnd   ! Snow fully "shades" melt ponds
264         ENDIF         
265         !
266         !                       !--- Specific snow fraction (for now, prescribed)
267         WHERE( ph_snw > 0._wp )   ;   zafrac_snw = 1.
268         ELSEWHERE                 ;   zafrac_snw = 0.
269         END WHERE
270         !
271         !                       !--- Specific ice fraction
272         zafrac_ice = 1. - zafrac_snw - zafrac_pnd
273         !
274         !                 !==  Snow-covered, pond-covered, and bare ice albedos  ==!
275         !
276         !                       !--- Bare ice albedo
277         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
278         z1_c2 = 1. / 0.05
279         WHERE    ( 1.5  < ph_ice                     )   ;   zalb_ice = zalb
280         ELSEWHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )   ;   zalb_ice = zalb       + ( 0.18 - zalb       ) * z1_c1 *  &
281            &                                                                       ( LOG(1.5) - LOG(ph_ice) )
282         ELSEWHERE                                        ;   zalb_ice = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice
283         END WHERE
284         !
285         z1_c1 = 1. / 0.02       !--- Snow-covered ice albedo (freezing, melting cases)
286         z1_c2 = 1. / 0.03
287         !
288         WHERE( pt_ice < rt0_snow )   ;   zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw * z1_c1 )
289         ELSEWHERE                    ;   zalb_snw = rn_alb_smlt - ( rn_alb_smlt - zalb_ice ) * EXP( - ph_snw * z1_c2 )
290         END WHERE
291         !
292         IF ( ld_pnd ) THEN      !--- Depth-dependent ponded ice albedo
293            z1_href_pnd = 0.05        ! inverse of the characteristic length scale (Lecomte et al. 2015)
294            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd * z1_href_pnd ) 
295         ELSE
296            zalb_pnd  = rn_alb_dpnd
297         ENDIF
298         !                       !--- Surface albedo is weighted mean of snow, ponds and bare ice contributions
299         pa_ice_os = zafrac_snw * zalb_snw  +  zafrac_pnd * zalb_pnd  +  zafrac_ice * zalb_ice
300         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 )
301
302      END SELECT
303      !
304      IF( nn_timing == 1 )   CALL timing_stop('icealb')
305      !
306   END SUBROUTINE ice_alb
307
308
309   SUBROUTINE ice_alb_init
310      !!----------------------------------------------------------------------
311      !!                 ***  ROUTINE alb_init  ***
312      !!
313      !! ** Purpose :   initializations for the albedo parameters
314      !!
315      !! ** Method  :   Read the namelist namalb
316      !!----------------------------------------------------------------------
317      INTEGER  ::   ios                 ! Local integer output status for namelist read
318      !!
319      NAMELIST/namalb/ nn_ice_alb, rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd
320      !!----------------------------------------------------------------------
321      !
322      REWIND( numnam_ice_ref )              ! Namelist namalb in reference namelist : Albedo parameters
323      READ  ( numnam_ice_ref, namalb, IOSTAT = ios, ERR = 901)
324901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namalb in reference namelist', lwp )
325
326      REWIND( numnam_ice_cfg )              ! Namelist namalb in configuration namelist : Albedo parameters
327      READ  ( numnam_ice_cfg, namalb, IOSTAT = ios, ERR = 902 )
328902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namalb in configuration namelist', lwp )
329      IF(lwm) WRITE ( numoni, namalb )
330      !
331      IF(lwp) THEN                      ! Control print
332         WRITE(numout,*)
333         WRITE(numout,*) 'ice_alb_init: set albedo parameters'
334         WRITE(numout,*) '~~~~~~~~~~~~'
335         WRITE(numout,*) '   Namelist namalb:'
336         WRITE(numout,*) '      choose the albedo parameterization   nn_ice_alb  = ', nn_ice_alb
337         WRITE(numout,*) '      albedo of dry snow                   rn_alb_sdry = ', rn_alb_sdry
338         WRITE(numout,*) '      albedo of melting snow               rn_alb_smlt = ', rn_alb_smlt
339         WRITE(numout,*) '      albedo of dry ice                    rn_alb_idry = ', rn_alb_idry
340         WRITE(numout,*) '      albedo of bare puddled ice           rn_alb_imlt = ', rn_alb_imlt
341         WRITE(numout,*) '      albedo of ponded ice                 rn_alb_dpnd = ', rn_alb_dpnd
342      ENDIF
343      !
344   END SUBROUTINE ice_alb_init
345
346#else
347   !!----------------------------------------------------------------------
348   !!   Default option           Dummy module         NO ESIM sea-ice model
349   !!----------------------------------------------------------------------
350#endif
351
352   !!======================================================================
353END MODULE icealb
Note: See TracBrowser for help on using the repository browser.