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.
albedoice.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/albedoice.F90 @ 8321

Last change on this file since 8321 was 8321, checked in by clem, 7 years ago

STEP3 (1): clean separation between sea-ice and ocean

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