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

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

changes in style - part2 -

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