source: palm/trunk/SOURCE/radiation_model.f90 @ 1526

 MODULE radiation_model_mod

!--------------------------------------------------------------------------------!
! This file is part of PALM.
!
! PALM is free software: you can redistribute it and/or modify it under the terms
! of the GNU General Public License as published by the Free Software Foundation,
! either version 3 of the License, or (at your option) any later version.
!
! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with
! PALM. If not, see <http://www.gnu.org/licenses/>.
!
! Copyright 1997-2014 Leibniz Universitaet Hannover
!--------------------------------------------------------------------------------!
!
! Current revisions:
! -----------------
! 
! 
! Former revisions:
! -----------------
! $Id: radiation_model.f90 1497 2014-12-02 17:28:07Z keck $
!
! 1496 2014-12-02 17:25:50Z maronga
! Initial revision
! 
!
! Description:
! ------------
! Radiation model(s), to be used e.g. with the land surface scheme
!------------------------------------------------------------------------------!

    USE arrays_3d,                                                             &
        ONLY: pt

    USE control_parameters,                                                    &
        ONLY: phi, surface_pressure, time_since_reference_point

    USE indices,                                                               &
        ONLY:  nxlg, nxrg, nyng, nysg, nzb_s_inner

    USE kinds


    IMPLICIT NONE

    INTEGER(iwp) :: i, j, k


    INTEGER(iwp) :: day_init = 1 !: day of the year at model start

    LOGICAL :: radiation = .FALSE.  !: flag parameter indicating wheather the radiation model is used

    REAL(wp), PARAMETER :: SW_0 = 1368.0, &       !: solar constant  
                           pi = 3.14159265358979323_wp, &
                           sigma_SB  = 5.67E-8_wp !: Stefan-Boltzmann constant
 
    REAL(wp) :: albedo = 0.2_wp,             & !: NAMELIST alpha
                dt_radiation = 9999999.9_wp, &
                exn,                         & !: Exner function
                lon = 0.0_wp,                & !: longitude in radians
                lat = 0.0_wp,                & !: latitude in radians
                decl_1,                      & !: declination coef. 1
                decl_2,                      & !: declination coef. 2
                decl_3,                      & !: declination coef. 3
                time_utc,                    & !: current time in UTC
                time_utc_init = 0.0_wp,      & !: UTC time at model start
                day,                         & !: current day of the year
                lambda = 0.0_wp,             & !: longitude in degrees
                declination,                 & !: solar declination angle
                hour_angle,                  & !: solar hour angle
                time_radiation = 0.0_wp,     &
                zenith,                      & !: solar zenith angle
                sky_trans                      !: sky transmissivity

    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
                alpha,                       & !: surface albedo
                Rn,                          & !: net radiation at the surface
                LW_in,                       & !: incoming longwave radiation
                LW_out,                      & !: outgoing longwave radiation
                SW_in,                       & !: incoming shortwave radiation
                SW_out                         !: outgoing shortwave radiation


    INTERFACE init_radiation
       MODULE PROCEDURE init_radiation
    END INTERFACE init_radiation

    INTERFACE lsm_radiation
       MODULE PROCEDURE lsm_radiation
    END INTERFACE lsm_radiation

    SAVE

    PRIVATE

    PUBLIC albedo, day_init, dt_radiation, init_radiation, lambda,             &
           lsm_radiation, Rn, radiation, SW_in, sigma_SB, time_radiation,      &
           time_utc_init 



 CONTAINS

!------------------------------------------------------------------------------!
! Description:
! ------------
!-- Initialization of the radiation model
!------------------------------------------------------------------------------!
    SUBROUTINE init_radiation
    

       IMPLICIT NONE

       ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
       ALLOCATE ( Rn(nysg:nyng,nxlg:nxrg) )
       ALLOCATE ( LW_in(nysg:nyng,nxlg:nxrg) )
       ALLOCATE ( LW_out(nysg:nyng,nxlg:nxrg) )
       ALLOCATE ( SW_in(nysg:nyng,nxlg:nxrg) )
       ALLOCATE ( SW_out(nysg:nyng,nxlg:nxrg) )

       alpha = albedo

!
!--    Calculate radiation scheme constants
       decl_1 = SIN(23.45_wp * pi / 180.0_wp)
       decl_2 = 2.0 * pi / 365.0_wp
       decl_3 = decl_2 * 81.0_wp

!
!--    Calculate latitude and longitude angles (lat and lon, respectively)
       lat = phi * pi / 180.0_wp
       lon = lambda * pi / 180.0_wp

       RETURN

    END SUBROUTINE init_radiation


!------------------------------------------------------------------------------!
! Description:
! ------------
!-- A simple clear sky radiation model
!------------------------------------------------------------------------------!
    SUBROUTINE lsm_radiation
    

       IMPLICIT NONE

!
!--    Calculate current day and time based on the initial values and simulation
!--    time
       day = day_init + FLOOR( (time_utc_init + time_since_reference_point)    &
                               / 86400.0_wp )
       time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)


!
!--    Calculate solar declination and hour angle   
       declination = ASIN( decl_1 * SIN(decl_2 * day - decl_3) )
       hour_angle  = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi

!
!--    Calculate zenith angle
       zenith = SIN(lat)*SIN(declination) + COS(lat) * COS(declination)        &
                                            * COS(hour_angle)
       zenith = MAX(0.0_wp,zenith)

!
!--    Calculate sky transmissivity
       sky_trans = 0.6_wp + 0.2_wp * zenith

!
!--    Calculate value of the Exner function
       exn = (surface_pressure / 1000.0_wp )**0.286_wp

!
!--    Calculate radiation fluxes and net radiation (Rn) for each grid point
       DO i = nxlg, nxrg
          DO j = nysg, nyng

             k = nzb_s_inner(j,i)
             SW_in(j,i)  = SW_0 * sky_trans * zenith
             SW_out(j,i) = - alpha(j,i) * SW_in(j,i)
             LW_out(j,i) = - sigma_SB * (pt(k,j,i) * exn)**4
             LW_in(j,i)  = 0.8 * sigma_SB * (pt(k+1,j,i) * exn)**4
             Rn(j,i)     = SW_in(j,i) + SW_out(j,i) + LW_in(j,i) + LW_out(j,i)

          ENDDO
       ENDDO

       RETURN

    END SUBROUTINE lsm_radiation



 END MODULE radiation_model_mod