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

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

changes in style - part1 - (now the code looks better txs to Gurvan's comments)

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