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

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

changes in style - part5 - start changing init routines

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