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

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

last routine names to be changed

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