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/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

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