source: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90 @ 11442

Last change on this file since 11442 was 11442, checked in by mattmartin, 21 months ago

Introduction of stochastic physics in NEMO, based on Andrea Storto's code.
For details, see ticket https://code.metoffice.gov.uk/trac/utils/ticket/251.

File size: 14.9 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   USE stopack
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   albedo_ice   ! routine called sbcice_lim.F90
30   PUBLIC   albedo_oce   ! routine called by ???
31
32   INTEGER  ::   albd_init = 0      !: control flag for initialization
33 
34   REAL(wp) ::   rmue     = 0.40    !  cosine of local solar altitude
35   REAL(wp) ::   ralb_oce = 0.066   ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
36   REAL(wp) ::   c1       = 0.05    ! snow thickness (only for nn_ice_alb=0)
37   REAL(wp) ::   c2       = 0.10    !  "        "
38   REAL(wp) ::   rcloud   = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0)
39 
40   !                             !!* namelist namsbc_alb
41   INTEGER  ::   nn_ice_alb
42   REAL(wp) ::   rn_albice
43
44   !!----------------------------------------------------------------------
45   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
46   !! $Id$
47   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
48   !!----------------------------------------------------------------------
49CONTAINS
50
51   SUBROUTINE albedo_ice( pt_ice, ph_ice, ph_snw, pa_ice_cs, pa_ice_os )
52      !!----------------------------------------------------------------------
53      !!               ***  ROUTINE albedo_ice  ***
54      !!         
55      !! ** Purpose :   Computation of the albedo of the snow/ice system
56      !!       
57      !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb)
58      !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies
59      !!                  1: the scheme is "home made" (for cloudy skies) and based on Brandt et al. (J. Climate 2005)
60      !!                                                                           and Grenfell & Perovich (JGR 2004)
61      !!                Description of scheme 1:
62      !!                  1) Albedo dependency on ice thickness follows the findings from Brandt et al (2005)
63      !!                     which are an update of Allison et al. (JGR 1993) ; Brandt et al. 1999
64      !!                     0-5cm  : linear function of ice thickness
65      !!                     5-150cm: log    function of ice thickness
66      !!                     > 150cm: constant
67      !!                  2) Albedo dependency on snow thickness follows the findings from Grenfell & Perovich (2004)
68      !!                     i.e. it increases as -EXP(-snw_thick/0.02) during freezing and -EXP(-snw_thick/0.03) during melting
69      !!                  3) Albedo dependency on clouds is speculated from measurements of Grenfell and Perovich (2004)
70      !!                     i.e. cloudy-clear albedo depend on cloudy albedo following a 2d order polynomial law
71      !!                  4) The needed 4 parameters are: dry and melting snow, freezing ice and bare puddled ice
72      !!
73      !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions:
74      !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo
75      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger
76      !!                     under melting conditions than under freezing conditions
77      !!                  3) the evolution of ice albedo as a function of ice thickness shows 
78      !!                     3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic
79      !!
80      !! References :   Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250.
81      !!                Brandt et al. 2005, J. Climate, vol 18
82      !!                Grenfell & Perovich 2004, JGR, vol 109
83      !!----------------------------------------------------------------------
84      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice      !  ice surface temperature (Kelvin)
85      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_ice      !  sea-ice thickness
86      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_snw      !  snow thickness
87      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_cs   !  albedo of ice under clear    sky
88      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_os   !  albedo of ice under overcast sky
89      !!
90      INTEGER  ::   ji, jj, jl         ! dummy loop indices
91      INTEGER  ::   ijpl               ! number of ice categories (3rd dim of ice input arrays)
92      REAL(wp)            ::   ralb_im, ralb_sf, ralb_sm, ralb_if
93      REAL(wp)            ::   zswitch, z1_c1, z1_c2
94      REAL(wp)                            ::   zalb_sm, zalb_sf, zalb_st ! albedo of snow melting, freezing, total
95      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalb_it             ! intermediate variable & albedo of ice (snow free)
96      !!---------------------------------------------------------------------
97
98      ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories
99     
100      CALL wrk_alloc( jpi,jpj,ijpl, zalb, zalb_it )
101
102      IF( albd_init == 0 )   CALL albedo_init      ! initialization
103
104     
105      SELECT CASE ( nn_ice_alb )
106
107      !------------------------------------------
108      !  Shine and Henderson-Sellers (1985)
109      !------------------------------------------
110      CASE( 0 )
111       
112         ralb_sf = 0.80       ! dry snow
113         ralb_sm = 0.65       ! melting snow
114         ralb_if = 0.72       ! bare frozen ice
115         ralb_im = rn_albice  ! bare puddled ice
116         
117         !  Computation of ice albedo (free of snow)
118         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im
119         ELSE WHERE                                              ;   zalb(:,:,:) = ralb_if
120         END  WHERE
121     
122         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
123         ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )
124         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zalb_it = 0.2467 + 0.7049 * ph_ice              &
125            &                                                                 - 0.8608 * ph_ice * ph_ice     &
126            &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice
127         ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice
128         END WHERE
129     
130         DO jl = 1, ijpl
131            DO jj = 1, jpj
132               DO ji = 1, jpi
133                  ! freezing snow
134                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2
135                  !                                        !  freezing snow       
136                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) )
137                  zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
138                     &                           + ph_snw(ji,jj,jl) * ( ralb_sf - zalb_it(ji,jj,jl) ) / c1  )   &
139                     &        +         zswitch   * ralb_sf 
140
141                  ! melting snow
142                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > c2
143                  zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) )
144                  zalb_sm = ( 1._wp - zswitch ) * ( ralb_im + ph_snw(ji,jj,jl) * ( ralb_sm - ralb_im ) / c2 )   &
145                      &     +         zswitch   *   ralb_sm 
146                  !
147                  ! snow albedo
148                  zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
149                  zalb_st  =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
150               
151                  ! Ice/snow albedo
152                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
153                  pa_ice_cs(ji,jj,jl) =  zswitch * zalb_st + ( 1._wp - zswitch ) * zalb_it(ji,jj,jl)
154                  !
155               END DO
156            END DO
157           
158            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) &
159               & CALL spp_gen( 1, pa_ice_cs(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl )
160                       
161         END DO
162
163         pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rcloud       ! Oberhuber correction for overcast sky
164
165      !------------------------------------------
166      !  New parameterization (2016)
167      !------------------------------------------
168      CASE( 1 ) 
169
170         ralb_im = rn_albice  ! bare puddled ice
171! compilation of values from literature
172         ralb_sf = 0.85      ! dry snow
173         ralb_sm = 0.75      ! melting snow
174         ralb_if = 0.60      ! bare frozen ice
175! Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved
176!         ralb_sf = 0.85       ! dry snow
177!         ralb_sm = 0.72       ! melting snow
178!         ralb_if = 0.65       ! bare frozen ice
179! Brandt et al 2005 (East Antarctica)
180!         ralb_sf = 0.87      ! dry snow
181!         ralb_sm = 0.82      ! melting snow
182!         ralb_if = 0.54      ! bare frozen ice
183!
184         !  Computation of ice albedo (free of snow)
185         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
186         z1_c2 = 1. / 0.05
187         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb = ralb_im
188         ELSE WHERE                                              ;   zalb = ralb_if
189         END  WHERE
190         
191         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
192         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = zalb     + ( 0.18 - zalb     ) * z1_c1 *  &
193            &                                                                     ( LOG(1.5) - LOG(ph_ice) )
194         ELSE WHERE                                       ;  zalb_it = ralb_oce + ( 0.18 - ralb_oce ) * z1_c2 * ph_ice
195         END WHERE
196
197         z1_c1 = 1. / 0.02
198         z1_c2 = 1. / 0.03
199         !  Computation of the snow/ice albedo
200         DO jl = 1, ijpl
201            DO jj = 1, jpj
202               DO ji = 1, jpi
203                  zalb_sf = ralb_sf - ( ralb_sf - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c1 );
204                  zalb_sm = ralb_sm - ( ralb_sm - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c2 );
205
206                   ! snow albedo
207                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
208                  zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
209
210                  ! Ice/snow albedo   
211                  zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
212                  pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl)
213
214              END DO
215            END DO
216           
217            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) &
218               & CALL spp_gen( 1, pa_ice_os(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl )
219           
220         END DO
221         ! Effect of the clouds (2d order polynomial)
222         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ); 
223
224      END SELECT
225     
226      CALL wrk_dealloc( jpi,jpj,ijpl, zalb, zalb_it )
227      !
228   END SUBROUTINE albedo_ice
229
230
231   SUBROUTINE albedo_oce( pa_oce_os , pa_oce_cs )
232      !!----------------------------------------------------------------------
233      !!               ***  ROUTINE albedo_oce  ***
234      !!
235      !! ** Purpose :   Computation of the albedo of the ocean
236      !!----------------------------------------------------------------------
237      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_os   !  albedo of ocean under overcast sky
238      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_cs   !  albedo of ocean under clear sky
239      !!
240      REAL(wp) :: zcoef 
241      !!----------------------------------------------------------------------
242      !
243      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982
244      pa_oce_cs(:,:) = zcoef 
245      pa_oce_os(:,:) = 0.06                       ! Parameterization of Kondratyev, 1969 and Payne, 1972
246      !
247   END SUBROUTINE albedo_oce
248
249
250   SUBROUTINE albedo_init
251      !!----------------------------------------------------------------------
252      !!                 ***  ROUTINE albedo_init  ***
253      !!
254      !! ** Purpose :   initializations for the albedo parameters
255      !!
256      !! ** Method  :   Read the namelist namsbc_alb
257      !!----------------------------------------------------------------------
258      INTEGER  ::   ios                 ! Local integer output status for namelist read
259      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice 
260      !!----------------------------------------------------------------------
261      !
262      albd_init = 1                     ! indicate that the initialization has been done
263      !
264      REWIND( numnam_ref )              ! Namelist namsbc_alb in reference namelist : Albedo parameters
265      READ  ( numnam_ref, namsbc_alb, IOSTAT = ios, ERR = 901)
266901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_alb in reference namelist', lwp )
267
268      REWIND( numnam_cfg )              ! Namelist namsbc_alb in configuration namelist : Albedo parameters
269      READ  ( numnam_cfg, namsbc_alb, IOSTAT = ios, ERR = 902 )
270902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_alb in configuration namelist', lwp )
271      IF(lwm) WRITE ( numond, namsbc_alb )
272      !
273      IF(lwp) THEN                      ! Control print
274         WRITE(numout,*)
275         WRITE(numout,*) 'albedo : set albedo parameters'
276         WRITE(numout,*) '~~~~~~~'
277         WRITE(numout,*) '   Namelist namsbc_alb : albedo '
278         WRITE(numout,*) '      choose the albedo parameterization                  nn_ice_alb = ', nn_ice_alb
279         WRITE(numout,*) '      albedo of bare puddled ice                          rn_albice  = ', rn_albice
280      ENDIF
281      !
282   END SUBROUTINE albedo_init
283
284   !!======================================================================
285END MODULE albedo
Note: See TracBrowser for help on using the repository browser.