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.
albedo.F90 in branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90 @ 8099

Last change on this file since 8099 was 8099, checked in by vancop, 7 years ago

impact of melt ponds on albedo, first run

  • Property svn:keywords set to Id
File size: 20.2 KB
Line 
1MODULE albedo
2   !!======================================================================
3   !!                       ***  MODULE  albedo  ***
4   !! Ocean forcing:  bulk thermohaline forcing of the ocean (or ice)
5   !!=====================================================================
6   !! History :  8.0  ! 2001-04  (LIM 1.0)
7   !!   NEMO     1.0  ! 2003-07  (C. Ethe, G. Madec)  Optimization (old name:shine)
8   !!             -   ! 2004-11  (C. Talandier)  add albedo_init
9   !!             -   ! 2001-06  (M. Vancoppenolle) LIM 3.0
10   !!             -   ! 2006-08  (G. Madec)  cleaning for surface module
11   !!            3.6  ! 2016-01  (C. Rousset) new parameterization for sea ice albedo
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   albedo_ice    : albedo for   ice (clear and overcast skies)
16   !!   albedo_oce    : albedo for ocean (clear and overcast skies)
17   !!   albedo_init   : initialisation of albedo computation
18   !!----------------------------------------------------------------------
19   USE phycst         ! physical constants
20   USE in_out_manager ! I/O manager
21   USE lib_mpp        ! MPP library
22   USE wrk_nemo       ! work arrays
23   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   albedo_ice   ! routine called sbcice_lim.F90
29   PUBLIC   albedo_oce   ! routine called by ???
30
31   INTEGER  ::   albd_init = 0      !: control flag for initialization
32 
33   REAL(wp) ::   rmue     = 0.40    !  cosine of local solar altitude
34   REAL(wp) ::   ralb_oce = 0.066   ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
35   REAL(wp) ::   c1       = 0.05    ! snow thickness (only for nn_ice_alb=0)
36   REAL(wp) ::   c2       = 0.10    !  "        "
37   REAL(wp) ::   rcloud   = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0)
38 
39   !                             !!* namelist namsbc_alb
40   INTEGER  ::   nn_ice_alb
41   REAL(wp) ::   rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd
42
43   !!----------------------------------------------------------------------
44   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
45   !! $Id$
46   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
47   !!----------------------------------------------------------------------
48CONTAINS
49
50   ! MV MP 2016
51   ! SUBROUTINE albedo_ice( pt_ice, ph_ice, ph_snw, pa_ice_cs, pa_ice_os )
52   SUBROUTINE albedo_ice( pt_ice, ph_ice, pa_ice, ph_snw, pa_pnd, pv_pnd, pa_ice_cs, pa_ice_os, ld_pnd )
53      !!----------------------------------------------------------------------
54      !!               ***  ROUTINE albedo_ice  ***
55      !!         
56      !! ** Purpose :   Computation of the albedo of the snow/ice system
57      !!       
58      !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb)
59      !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies
60      !!                  1: the scheme is "home made" (for cloudy skies) and based on Brandt et al. (J. Climate 2005)
61      !!                                                                           and Grenfell & Perovich (JGR 2004)
62      !!                Description of scheme 1:
63      !!                  1) Albedo dependency on ice thickness follows the findings from Brandt et al (2005)
64      !!                     which are an update of Allison et al. (JGR 1993) ; Brandt et al. 1999
65      !!                     0-5cm  : linear function of ice thickness
66      !!                     5-150cm: log    function of ice thickness
67      !!                     > 150cm: constant
68      !!                  2) Albedo dependency on snow thickness follows the findings from Grenfell & Perovich (2004)
69      !!                     i.e. it increases as -EXP(-snw_thick/0.02) during freezing and -EXP(-snw_thick/0.03) during melting
70      !!                  3) Albedo dependency on clouds is speculated from measurements of Grenfell and Perovich (2004)
71      !!                     i.e. cloudy-clear albedo depend on cloudy albedo following a 2d order polynomial law
72      !!                  4) The needed 4 parameters are: dry and melting snow, freezing ice and bare puddled ice
73      !!
74      !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions:
75      !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo
76      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger
77      !!                     under melting conditions than under freezing conditions
78      !!                  3) the evolution of ice albedo as a function of ice thickness shows 
79      !!                     3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic
80      !!
81      !! References :   Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250.
82      !!                Brandt et al. 2005, J. Climate, vol 18
83      !!                Grenfell & Perovich 2004, JGR, vol 109
84      !!----------------------------------------------------------------------
85      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice      !  ice surface temperature (Kelvin)
86      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_ice      !  sea-ice thickness
87      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_snw      !  snow depth
88      ! MV MP 2016
89      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pa_ice      !  ice fraction
90      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pa_pnd      !  melt pond fraction
91      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pv_pnd      !  melt pond volume
92      LOGICAL , INTENT(in   )                   ::   ld_pnd      !  melt ponds radiatively active
93      ! END MV MP 2016
94      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_cs   !  albedo of ice under clear    sky
95      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_os   !  albedo of ice under overcast sky
96      !!
97      INTEGER  ::   ji, jj, jl         ! dummy loop indices
98      INTEGER  ::   ijpl               ! number of ice categories (3rd dim of ice input arrays)
99      REAL(wp)                            ::   ralb_im, ralb_sf, ralb_sm, ralb_if, ralb_dp ! MV MP 2016
100      REAL(wp)                            ::   zswitch, z1_c1, z1_c2
101      REAL(wp)                            ::   zhref_pnd                                   ! MV MP 2016
102      REAL(wp)                            ::   zalb_sm, zalb_sf, zalb_st ! albedo of snow melting, freezing, total
103      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalb_it             ! intermediate variable & albedo of ice (snow free)
104      ! MV MP 2016
105      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zafrac_pnd, zh_pnd        ! melt pond fraction and thickness
106      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb_pnd                  ! ponded sea ice albedo
107      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb_ice                  ! bare sea ice albedo
108      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb_snw                  ! snow-covered sea ice albedo
109      ! END MV MP 2016
110      !!---------------------------------------------------------------------
111      ! TO DO
112      ! Update namelist
113      ! V Change the reference for rn_albice - should be bare ice with no puddles!!!
114      ! ... add 3rd computation
115
116
117      ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories
118     
119      CALL wrk_alloc( jpi,jpj,ijpl, zalb, zalb_it )
120      CALL wrk_alloc( jpi,jpj,ijpl, zafrac_pnd, zh_pnd, zalb_pnd )
121
122      IF( albd_init == 0 )   CALL albedo_init      ! initialization
123
124      ralb_sf = rn_alb_sdry ! dry snow
125      ralb_sm = rn_alb_smlt ! melting snow
126      ralb_if = rn_alb_idry ! bare frozen ice
127      ralb_im = rn_alb_imlt ! bare puddled ice
128      ralb_dp = rn_alb_dpnd ! deep pond albedo   ! MV MP 2016
129     
130      ! MV MP 2016
131      !------------------------------------------------
132      !  Pond fraction and thickness, snow-free albedo
133      !------------------------------------------------
134      ! If melt ponds are radiatively active ( ld_pnd = .TRUE. )
135      ! Then - melt pond fraction and thickness are assumed non-zero in the albedo computation
136      !      - the snow-free albedo is the weighted mean of bare ice (ralb_if) and ponded-ice albedos (zalb_pnd)
137
138      ! Otherwise, we just specify the bare ice albedo
139      ! for cold ice and melting ice, the latter implicitly assuming a specified fraction of melt ponds
140
141      IF ( ld_pnd ) THEN
142
143         ! ponded ice fraction and thickness
144         WHERE ( pa_ice > 1.e-10_wp ) ;   zafrac_pnd = pa_pnd / pa_ice ;   zh_pnd = pv_pnd / pa_ice
145         ELSE WHERE                   ;   zafrac_pnd = 0._wp           ;   zh_pnd = 0._wp 
146         END WHERE
147         
148         ! ponded ice albedo, including thickness dependence:
149         zhref_pnd = 0.05         ! characteristic length scale for albedo dependence of ponded ice, Lecomte et al (2015)
150         zalb_pnd  = ralb_dp - ( ralb_dp - ralb_if ) * EXP( - zh_pnd / zhref_pnd )
151
152         ! weighted bare / ponded ice albedo (no ice thickness correction)
153         WHERE     ( ph_snw == 0._wp )  ;   zalb = ralb_if * ( 1. - zafrac_pnd ) + zalb_pnd * zafrac_pnd
154         END  WHERE
155
156      ELSE ! No melt ponds or radiatively inactive
157     
158         ! ponded ice fraction and thickness
159         zafrac_pnd = 0._wp
160         zh_pnd     = 0._wp
161
162         ! snow-free ice albedo
163         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb = ralb_im    ! prescribed bare ice albedo with assumed contribution of melt ponds
164         ELSE WHERE                                              ;   zalb = ralb_if
165         END  WHERE
166
167      ENDIF
168
169      ! END MV MP 2016
170
171      SELECT CASE ( nn_ice_alb )
172
173      !------------------------------------------
174      !  Shine and Henderson-Sellers (1985)
175      !------------------------------------------
176      CASE( 0 )
177       
178         ! Reference (clear sky) values
179         !ralb_sf = 0.80       ! dry snow
180         !ralb_sm = 0.65       ! melting snow
181         !ralb_if = 0.72       ! bare frozen ice
182         !ralb_im = rn_albice  ! bare ice
183         !ralb_dp = 0.30       ! deep ponded ice albedo, clear sky value, Lecomte et al (2015) MV MP 2016
184         !                     ! early melt ponds 0.24, late melt ponds 0.13 Grenfell & Perovich (2004)
185         
186         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
187         ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )
188         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zalb_it = 0.2467 + 0.7049 * ph_ice              &
189            &                                                                 - 0.8608 * ph_ice * ph_ice     &
190            &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice
191         ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice
192         END WHERE
193     
194         DO jl = 1, ijpl
195            DO jj = 1, jpj
196               DO ji = 1, jpi
197                  ! freezing snow
198                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2
199                  !                                        !  freezing snow       
200                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) )
201                  zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
202                     &                           + ph_snw(ji,jj,jl) * ( ralb_sf - zalb_it(ji,jj,jl) ) / c1  )   &
203                     &        +         zswitch   * ralb_sf 
204
205                  ! melting snow
206                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > c2
207                  zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) )
208                  zalb_sm = ( 1._wp - zswitch ) * ( ralb_im + ph_snw(ji,jj,jl) * ( ralb_sm - ralb_im ) / c2 )   &
209                      &     +         zswitch   *   ralb_sm 
210                  !
211                  ! snow albedo
212                  zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
213                  zalb_st  =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
214               
215                  ! surface albedo
216                  zswitch             = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
217                  pa_ice_cs(ji,jj,jl) = zswitch * zalb_st + ( 1._wp - zswitch ) * zalb_it(ji,jj,jl)
218
219               END DO
220            END DO
221         END DO
222
223         pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rcloud       ! Oberhuber correction for overcast sky
224
225      !------------------------------------------
226      !  New parameterization (2016)
227      !------------------------------------------
228      CASE( 1 ) 
229
230!         ! Reference (overcast sky) values
231!       ralb_im = rn_albice  ! bare ice
232! compilation of values from literature
233!        ralb_sf = 0.85      ! dry snow
234!        ralb_sm = 0.75      ! melting snow
235!        ralb_if = 0.60      ! bare frozen ice
236!        ralb_dp = 0.36      ! ponded-ice overcast albedo (Lecomte et al, 2015)
237!                            ! early melt pnds 0.27, late melt ponds 0.14 Grenfell & Perovich
238! Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved
239!         ralb_sf = 0.85       ! dry snow
240!         ralb_sm = 0.72       ! melting snow
241!         ralb_if = 0.65       ! bare frozen ice
242! Brandt et al 2005 (East Antarctica)
243!         ralb_sf = 0.87      ! dry snow
244!         ralb_sm = 0.82      ! melting snow
245!         ralb_if = 0.54      ! bare frozen ice
246!
247         ! Computation of ice albedo (free of snow)
248         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
249         z1_c2 = 1. / 0.05
250
251         ! Accounting for the ice-thickness dependency (this has not been carefully checked with melt ponds)
252         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
253         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = zalb     + ( 0.18 - zalb     ) * z1_c1 *  &
254            &                                                                     ( LOG(1.5) - LOG(ph_ice) )
255         ELSE WHERE                                       ;  zalb_it = ralb_oce + ( 0.18 - ralb_oce ) * z1_c2 * ph_ice
256         END WHERE
257
258         z1_c1 = 1. / 0.02
259         z1_c2 = 1. / 0.03
260         
261         ! Surface albedo (accounting for snow, ice melt ponds), overcast skies
262         DO jl = 1, ijpl
263            DO jj = 1, jpj
264               DO ji = 1, jpi
265
266                  zalb_sf = ralb_sf - ( ralb_sf - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c1 );
267                  zalb_sm = ralb_sm - ( ralb_sm - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c2 );
268
269                   ! snow albedo
270                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
271                  zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
272
273                  ! Ice/snow albedo   
274                  zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
275                  pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl) 
276
277              END DO
278            END DO
279         END DO
280
281         ! Clear sky albedo
282         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ); 
283
284!     !------------------------------------------
285!     !  Surface-based parameterization (2017)
286!     !------------------------------------------
287!      CASE( 2 )
288!
289!         !-----------------------------------------
290!         ! Snow, bare ice and ponded ice fractions
291!         !-----------------------------------------
292!         ! Specific fractions (zafrac) refer to relative area covered by snow within each ice category
293!
294!         !--- Snow (for now, prescribed)
295!         WHERE     ( pa_ice > epsi10 .OR. ph_snw <= epsi06 )  ;  zafrac_snw = 0.
296!         ELSE WHERE                                           ;  zafrac_snw = 1.
297!         END WHERE
298!
299!         !--- Melt ponds
300!         ! for now, done at the beginning
301!
302!         !--- Ice
303!         WHERE     ( pa_ice > epsi10 )   ;  zafrac_ice = 1. - zafrac_snw - zafrac_pnd
304!         ELSE WHERE                      ;  zafrac_ice = 0.
305!         ENDWHERE
306!
307!         !--------------------------------------------------
308!         ! Snow-covered, pond-covered, and bare ice albedos
309!         !--------------------------------------------------
310!         ! zalb_snw, zalb_pnd, zalb_ice
311!         k
312!
313!
314!         
315!
316!         
317!
318!
319!
320!      ENDch grid cell has specific fractions of snow, ponds and ice
321!      ãs = 1 if hs > 0; 0 otherwise
322!      ãp = a_pnd / a_i
323!      ãi = ( 1 - ãs - ãp )
324!     
325!      2) Compute albedo for each of the three surface types
326!      αi = f(hi,Tsu)
327!      αs = f(hi, hs, Tsu)
328!      αp = f(hp, hi?, Tsu?)
329!     
330!      3) Surface albedo = weighted mean
331!     
332!      α = αs * ãs + αp * ãp + αi * ( 1 - ãs - ãp)
333
334      END SELECT
335     
336      CALL wrk_dealloc( jpi,jpj,ijpl, zalb, zalb_it )
337      CALL wrk_dealloc( jpi,jpj,ijpl, zafrac_pnd, zh_pnd ) ! MV MP 2016
338      !
339   END SUBROUTINE albedo_ice
340
341
342   SUBROUTINE albedo_oce( pa_oce_os , pa_oce_cs )
343      !!----------------------------------------------------------------------
344      !!               ***  ROUTINE albedo_oce  ***
345      !!
346      !! ** Purpose :   Computation of the albedo of the ocean
347      !!----------------------------------------------------------------------
348      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_os   !  albedo of ocean under overcast sky
349      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_cs   !  albedo of ocean under clear sky
350      !!
351      REAL(wp) :: zcoef 
352      !!----------------------------------------------------------------------
353      !
354      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982
355      pa_oce_cs(:,:) = zcoef 
356      pa_oce_os(:,:) = 0.06                       ! Parameterization of Kondratyev, 1969 and Payne, 1972
357      !
358   END SUBROUTINE albedo_oce
359
360
361   SUBROUTINE albedo_init
362      !!----------------------------------------------------------------------
363      !!                 ***  ROUTINE albedo_init  ***
364      !!
365      !! ** Purpose :   initializations for the albedo parameters
366      !!
367      !! ** Method  :   Read the namelist namsbc_alb
368      !!----------------------------------------------------------------------
369      INTEGER  ::   ios                 ! Local integer output status for namelist read
370      NAMELIST/namsbc_alb/ nn_ice_alb, rn_alb_sdry, rn_alb_smlt, rn_alb_idry , rn_alb_imlt, rn_alb_dpnd
371      !!----------------------------------------------------------------------
372      !
373      albd_init = 1                     ! indicate that the initialization has been done
374      !
375      REWIND( numnam_ref )              ! Namelist namsbc_alb in reference namelist : Albedo parameters
376      READ  ( numnam_ref, namsbc_alb, IOSTAT = ios, ERR = 901)
377901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_alb in reference namelist', lwp )
378
379      REWIND( numnam_cfg )              ! Namelist namsbc_alb in configuration namelist : Albedo parameters
380      READ  ( numnam_cfg, namsbc_alb, IOSTAT = ios, ERR = 902 )
381902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_alb in configuration namelist', lwp )
382      IF(lwm) WRITE ( numond, namsbc_alb )
383      !
384      IF(lwp) THEN                      ! Control print
385         WRITE(numout,*)
386         WRITE(numout,*) 'albedo : set albedo parameters'
387         WRITE(numout,*) '~~~~~~~'
388         WRITE(numout,*) '   Namelist namsbc_alb : albedo '
389         WRITE(numout,*) '      choose the albedo parameterization                  nn_ice_alb = ', nn_ice_alb
390         WRITE(numout,*) '      albedo of dry snow                                  rn_alb_sdry = ', rn_alb_sdry
391         WRITE(numout,*) '      albedo of melting snow                              rn_alb_smlt = ', rn_alb_smlt
392         WRITE(numout,*) '      albedo of dry ice                                   rn_alb_idry = ', rn_alb_idry
393         WRITE(numout,*) '      albedo of bare ice                                  rn_alb_imlt = ', rn_alb_imlt
394         WRITE(numout,*) '      albedo of ponded ice                                rn_alb_dpnd = ', rn_alb_dpnd
395      ENDIF
396      !
397   END SUBROUTINE albedo_init
398
399   !!======================================================================
400END MODULE albedo
Note: See TracBrowser for help on using the repository browser.