ocean_albedo.F90

Go to the documentation of this file.

!
! $Id$
!

!     #########

subroutine ocean_albedo(knon,zrmu0,knindex,pwind,SFRWL,alb_dir_new,alb_dif_new)


!     ##################################################################
!
!!****  *ALBEDO_RS14*  
!!
!!    PURPOSE
!!    -------
!       computes the direct & diffuse albedo over open water
!
!     
!!**  METHOD
!!    ------
!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!      
!!    AUTHOR
!!    ------
!!              R. Séférian           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    03/2014
!                   05/2014 R. Séférian & B. Decharme :: Adaptation to spectral
!                   computation for diffuse and direct albedo
!                   08/2014 S. Baek :: for wider wavelength range 200-4000nm and
!                   adaptation to LMDZ + whitecap effect by Koepke + chrolophyll
!                   map from climatology file
!       
!-------------------------------------------------------------------------------
!
!*           DECLARATIONS
!            ------------
!
USE ocean_albedo_para
use dimphy
!LF USE PARKIND1  ,ONLY : JPIM     ,JPRB
use phys_state_var_mod, only : chl_con
!
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!              -------------------------
!

include "clesphys.h"

integer, intent(in) :: knon
integer, dimension(klon), intent(in) :: knindex
real, dimension(klon), intent(in) :: zrmu0,pwind
real, dimension(klon,nsw), intent(out) ::  alb_dir_new,alb_dif_new
real, dimension(6),intent(in) :: SFRWL


!=== LOCAL VARIABLES

REAL, parameter :: XPI=4.*atan(1.)

!
!*      0.2    declarations of local variables
!              -------------------------
!
REAL, DIMENSION(klon)               :: ZCHL        ! surface chlorophyll
REAL, DIMENSION(klon,NNWL)               :: ZDIR_ALB    ! direct  ocean surface albedo (spectral)
REAL, DIMENSION(klon,NNWL)               :: ZSCA_ALB    ! diffuse ocean surface albedo (spectral)
!
INTEGER                         :: JI, JWL                  ! indexes
REAL                            :: ZWL                      ! input parameter: wavelength and diffuse/direct fraction of light
REAL:: ZSIG, ZREFM, ZXX2, ZR00, ZRR0, ZRRR                        ! computation variables
REAL:: ZR22, ZUE, ZUE2, ZR11DF, ZALBT, ZFWC                       ! computation variables
REAL:: ZCHLABS, ZAW, ZBW, ZAP, ZYLMD, ZBP550                      ! computation variables
REAL:: ZBBP, ZNU, ZHB                                       ! computation variables
REAL:: ZCOSZEN                                              ! Cosine of the zenith solar angle
REAL:: ZR11, ZRW, ZRWDF, ZRDF                                   ! 4 components of the OSA
! new damping coefficient
REAL:: ZDAMP     

! 
REAL            :: ZWORK                           ! dummy variable
! 
!LF REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
!



!
alb_dir_new(:,:) = 0. 
alb_dif_new(:,:) = 0. 
!
zdir_alb(:,:) = 0. 
zsca_alb(:,:) = 0.
!
!

!ZCHL(:) = CHL_CON!0.05 ! averaged global values for surface chlorophyll
if(ok_chlorophyll)then
  do ji=1,knon
  zchl(ji)=chl_con(knindex(ji))
  enddo
else 
  zchl(:) = 0.05
endif


!
DO jwl=1,nnwl           ! loop over the wavelength
!
  DO ji=1,knon  ! loop over the grid points


    !---------------------------------------------------------------------------------
    ! 0- Compute baseline values
    !---------------------------------------------------------------------------------
    
    ! Get refractive index for the correspoding wavelength
    zwl=xakwl(jwl) !!!----------  wavelength value
    zrefm= xakrefm(jwl) !!!--------- refraction index value
  

    ! compute the cosine of the solar zenith angle
!    ZCOSZEN = COS(XPI/2 - PZENITH(JI))
     zcoszen = zrmu0(knindex(ji))
    ! Compute sigma derived from wind speed (Cox & Munk reflectance model)
    zsig=sqrt(0.003+0.00512*pwind(ji))

    
    !---------------------------------------------------------------------------------
    ! 1- Compute direct surface albedo (ZR11)
    !---------------------------------------------------------------------------------
    !
    zxx2=sqrt(1.0-(1.0-zcoszen**2)/zrefm**2)
    zrr0=0.50*(((zxx2-zrefm*zcoszen)/(zxx2+zrefm*zcoszen))**2 +((zcoszen-zrefm*zxx2)/(zcoszen+zrefm*zxx2))**2)
    zrrr=0.50*(((zxx2-1.34*zcoszen)/(zxx2+1.34*zcoszen))**2 +((zcoszen-1.34*zxx2)/(zcoszen+1.34*zxx2))**2)
    zr11=zrr0-(0.0152-1.7873*zcoszen+6.8972*zcoszen**2-8.5778*zcoszen**3+4.071*zsig-7.6446*zcoszen*zsig) &
        & * exp(0.1643-7.8409*zcoszen-3.5639*zcoszen**2-2.3588*zsig+10.0538*zcoszen*zsig)*zrr0/zrrr
    ! 
    !---------------------------------------------------------------------------------
    ! 2- Compute surface diffuse albedo (ZRDF)
    !---------------------------------------------------------------------------------
    ! Diffuse albedo from Jin et al., 2006 + estimation from diffuse fraction of
    ! light (relying later on AOD)
    zrdf=-0.1482-0.012*zsig+0.1609*zrefm-0.0244*zsig*zrefm ! surface diffuse (Eq 5a-5b)
   
    !---------------------------------------------------------------------------------
    ! *- Determine absorption and backscattering
    ! coefficients to determine reflectance below the surface (Ro) once for all
    !
    ! *.1- Absorption by chlorophyll
    zchlabs= xakachl(jwl) 
    ! *.2- Absorption by seawater 
    zaw= xakaw3(jwl) 
    ! *.3- Backscattering by seawater
    zbw= xakbw(jwl) 
    ! *.4- Backscattering by chlorophyll
    zylmd = exp(0.014*(440.0-zwl))
    zwork= exp(log(zchl(ji))*0.65)
    zap = 0.06*zchlabs*zwork +0.2*(xaw440+0.06*zwork)*zylmd
    zbp550 = 0.416 * exp(log(zchl(ji))*0.766)
   
    IF ( zchl(ji) > 2. ) THEN
      znu=0.
    ELSE
      IF ( zchl(ji) > 0.02 ) THEN
        zwork=log10(zchl(ji))
        znu=0.5*(zwork-0.3)
        zbbp=(0.002+0.01*(0.5-0.25*zwork)*(zwl/550.)**znu)*zbp550
      ELSE
        zbbp=0.019*(550./zwl)*zbp550       !ZBBPf=0.0113 at chl<=0.02
      ENDIF
    ENDIF
   
    ! Morel-Gentili(1991), Eq (12)
    ! ZHB=h/(h+2*ZBBPf*(1.-h))        
    zhb=0.5*zbw/(0.5*zbw+zbbp)
   
    !---------------------------------------------------------------------------------
    ! 3- Compute direct water-leaving albedo (ZRW)
    !---------------------------------------------------------------------------------
    ! Based on Morel & Gentilli 1991 parametrization
    zr22=0.48168549-0.014894708*zsig-0.20703885*zsig**2
    ! Use Morel 91 formula to compute the direct reflectance
    ! below the surface
    zr00=(0.5*zbw+zbbp)/(zaw+zap) *(0.6279-0.2227*zhb-0.0513*zhb**2 + (-0.3119+0.2465*zhb)*zcoszen)
    zrw=zr00*(1.-zr22)*(1.-zr11)/(1.-zr00*zr22)

    zrw=zr00*(1.-zr22)/(1.-zr00*zr22)
    !---------------------------------------------------------------------------------
    ! 4- Compute diffuse water-leaving albedo (ZRWDF)
    !---------------------------------------------------------------------------------
    ! as previous water-leaving computation but assumes a uniform incidence of
    ! shortwave at surface (ue)
    zue=0.676               ! equivalent u_unif for diffuse incidence
    zue2=sqrt(1.0-(1.0-zue**2)/zrefm**2)
    zrr0=0.50*(((zue2-zrefm*zue)/(zue2+zrefm*zue))**2 +((zue-zrefm*zue2)/(zue+zrefm*zue2))**2)
    zrrr=0.50*(((zue2-1.34*zue)/(zue2+1.34*zue))**2 +((zue-1.34*zue2)/(zue+1.34*zue2))**2)
    zr11df=zrr0-(0.0152-1.7873*zue+6.8972*zue**2-8.5778*zue**3+4.071*zsig-7.6446*zue*zsig) &
          & * exp(0.1643-7.8409*zue-3.5639*zue**2-2.3588*zsig+10.0538*zue*zsig)*zrr0/zrrr
    ! Use Morel 91 formula to compute the diffuse
    ! reflectance below the surface
    zr00=(0.5*zbw+zbbp)/(zaw+zap) *(0.6279-0.2227*zhb-0.0513*zhb**2 + (-0.3119+0.2465*zhb)*zue)
    zrwdf=zr00*(1.-zr22)*(1.-zr11df)/(1.-zr00*zr22)
   
    ! original : correction for foam (Eq 16-17)
    zfwc=3.97e-4*pwind(ji)**(1.59) ! Salisbury 2014 eq(2) at 37GHz, value in fraction
    ! has to be update once we have information from wave model (discussion with G. Madec)
     
    ! --------------------------------------------------------------------
    !  *- OSA estimation
    ! --------------------------------------------------------------------
    ! partitionning direct and diffuse albedo
    !

    ! excluding diffuse albedo ZRW on ZDIR_ALB
      zdir_alb(ji,jwl) =  xfrwl(jwl) *((1.-zfwc) * (zr11+zrw) +zfwc*xrwc(jwl))
      zsca_alb(ji,jwl) =  xfrwl(jwl) *((1.-zfwc) * (zrdf+zrwdf) + zfwc*xrwc(jwl))

    !  print*,ji,ZFWC,ZDIR_ALB(JI,JWL),ZSCA_ALB(JI,JWL),pwind(ji)
    ENDDO ! end of the loop over grid points

ENDDO ! ending loop over wavelengths


! integral for each nsw band

select case(nsw)
case(2)
  do ji=1,knon
   alb_dir_new(ji,1)=sum(zdir_alb(ji,1:49))/sfrwl(1) ! from 200nm to 680nm 
   alb_dir_new(ji,2)=sum(zdir_alb(ji,50:381))/sfrwl(2) ! from 690nm to 4000 nm

   alb_dif_new(ji,1)=sum(zsca_alb(ji,1:49))/sfrwl(1) ! from 200nm to 680nm 
   alb_dif_new(ji,2)=sum(zsca_alb(ji,50:381))/sfrwl(2) ! from 690nm to 4000 nm
  enddo
case(4)
  do ji=1,knon
   alb_dir_new(ji,1)=sum(zdir_alb(ji,1:49))/sfrwl(1) ! from 200nm to 680nm 
   alb_dir_new(ji,2)=sum(zdir_alb(ji,50:99))/sfrwl(2) ! from 690nm to 1180 nm
   alb_dir_new(ji,3)=sum(zdir_alb(ji,100:218))/sfrwl(3) ! from 1190nm to 2370 nm
   alb_dir_new(ji,4)=sum(zdir_alb(ji,219:381))/sfrwl(4) ! from 2380nm to 4000 nm

   alb_dif_new(ji,1)=sum(zsca_alb(ji,1:49))/sfrwl(1) ! from 200nm to 680nm 
   alb_dif_new(ji,2)=sum(zsca_alb(ji,50:99))/sfrwl(2) ! from 690nm to 1180 nm
   alb_dif_new(ji,3)=sum(zsca_alb(ji,100:218))/sfrwl(3) ! from 1190nm to 2370 nm
   alb_dif_new(ji,4)=sum(zsca_alb(ji,219:381))/sfrwl(4) ! from 2380nm to 4000 nm
  enddo
case(6)
 do ji=1,knon
   alb_dir_new(ji,1)=sum(zdir_alb(ji,1:5))/sfrwl(1) ! from 200nm to 240nm 
   alb_dir_new(ji,2)=sum(zdir_alb(ji,6:24))/sfrwl(2) ! from 250nm to 430 nm
   alb_dir_new(ji,3)=sum(zdir_alb(ji,25:49))/sfrwl(3) ! from 440nm to 680 nm
   alb_dir_new(ji,4)=sum(zdir_alb(ji,50:99))/sfrwl(4) ! from 690nm to 1180 nm
   alb_dir_new(ji,5)=sum(zdir_alb(ji,100:218))/sfrwl(5) ! from 1190nm to 2370 nm
   alb_dir_new(ji,6)=sum(zdir_alb(ji,219:381))/sfrwl(6) ! from 2380nm to 4000 nm

   alb_dif_new(ji,1)=sum(zsca_alb(ji,1:5))/sfrwl(1) ! from 200nm to 240nm 
   alb_dif_new(ji,2)=sum(zsca_alb(ji,6:24))/sfrwl(2) ! from 250nm to 430 nm
   alb_dif_new(ji,3)=sum(zsca_alb(ji,25:49))/sfrwl(3) ! from 440nm to 680 nm
   alb_dif_new(ji,4)=sum(zsca_alb(ji,50:99))/sfrwl(4) ! from 690nm to 1180 nm
   alb_dif_new(ji,5)=sum(zsca_alb(ji,100:218))/sfrwl(5) ! from 1190nm to 2370 nm
   alb_dif_new(ji,6)=sum(zsca_alb(ji,219:381))/sfrwl(6) ! from 2380nm to 4000 nm
 enddo
end select



END subroutine ocean_albedo