source: LMDZ5/trunk/libf/phylmd/ocean_albedo.F90 @ 2483

!
! $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