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

 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 1552 2015-03-03 14:27:15Z keck $
!
! 1551 2015-03-03 14:18:16Z maronga
! Added support for data output. Various variables have been renamed. Added
! interface for different radiation schemes (currently: clear-sky, constant, and
! RRTM (not yet implemented).
! 
! 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

    USE netcdf_control,                                                        &
        ONLY:  dots_label, dots_num, dots_unit


    IMPLICIT NONE

    CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtm'

    INTEGER(iwp) :: i, j, k


    INTEGER(iwp) :: day_init     = 172,  & !: day of the year at model start (21/06)
                    dots_rad     = 0,    & !: starting index for timeseries output 
                    irad_scheme  = 0

    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 = 0.0_wp,       & !: radiation model timestep
                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 = 43200.0_wp,  & !: UTC time at model start (noon)
                day,                         & !: current day of the year
                lambda = 0.0_wp,             & !: longitude in degrees
                declination,                 & !: solar declination angle
                net_radiation = 0.0_wp,      & !: net radiation at surface
                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
                rad_net,                     & !: net radiation at the surface
                rad_net_av,                  & !: average of rad_net
                rad_lw_in,                   & !: incoming longwave radiation
                rad_lw_out,                  & !: outgoing longwave radiation
                rad_sw_in,                   & !: incoming shortwave radiation
                rad_sw_in_av,                & !: average of rad_sw_in
                rad_sw_out                     !: outgoing shortwave radiation


    INTERFACE init_radiation
       MODULE PROCEDURE init_radiation
    END INTERFACE init_radiation

    INTERFACE radiation_clearsky
       MODULE PROCEDURE radiation_clearsky
    END INTERFACE radiation_clearsky

    INTERFACE radiation_constant
       MODULE PROCEDURE radiation_constant
    END INTERFACE radiation_constant

    INTERFACE radiation_rrtm
       MODULE PROCEDURE radiation_rrtm
    END INTERFACE radiation_rrtm


    SAVE

    PRIVATE

    PUBLIC albedo, day_init, dots_rad, dt_radiation, init_radiation,           &
           irad_scheme, lambda, net_radiation, rad_net, rad_net_av, radiation, &
           radiation_clearsky, radiation_constant, radiation_rrtm,             &
           radiation_scheme, rad_sw_in, rad_sw_in_av, 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 ( rad_net(nysg:nyng,nxlg:nxrg) )
       ALLOCATE ( rad_lw_in(nysg:nyng,nxlg:nxrg) )
       ALLOCATE ( rad_lw_out(nysg:nyng,nxlg:nxrg) )
       ALLOCATE ( rad_sw_in(nysg:nyng,nxlg:nxrg) )
       ALLOCATE ( rad_sw_out(nysg:nyng,nxlg:nxrg) )

       rad_sw_in  = 0.0_wp
       rad_sw_out = 0.0_wp
       rad_lw_in  = 0.0_wp
       rad_lw_out = 0.0_wp
       rad_net    = 0.0_wp

       alpha = albedo

!
!--    Fix net radiation in case of radiation_scheme = 'constant'
       IF ( irad_scheme == 0 )  THEN
          rad_net = net_radiation
!
!--    Calculate radiation scheme constants
       ELSEIF ( irad_scheme == 1 .OR. irad_scheme == 2 )  THEN
          decl_1 = SIN(23.45_wp * pi / 180.0_wp)
          decl_2 = 2.0_wp * 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
       ENDIF
!
!--    Add timeseries for radiation model
       dots_label(dots_num+1) = "rad_net"
       dots_label(dots_num+2) = "rad_sw_in"
       dots_unit(dots_num+1:dots_num+2) = "W/m2"

       dots_rad  = dots_num + 1
       dots_num  = dots_num + 2

       RETURN

    END SUBROUTINE init_radiation


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

       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 (rad_net) for each grid 
!--    point
       DO i = nxlg, nxrg
          DO j = nysg, nyng

             k = nzb_s_inner(j,i)
             rad_sw_in(j,i)  = SW_0 * sky_trans * zenith
             rad_sw_out(j,i) = - alpha(j,i) * rad_sw_in(j,i)
             rad_lw_out(j,i) = - sigma_SB * (pt(k,j,i) * exn)**4
             rad_lw_in(j,i)  = 0.8_wp * sigma_SB * (pt(k+1,j,i) * exn)**4
             rad_net(j,i)    = rad_sw_in(j,i) + rad_sw_out(j,i)                &
                                + rad_lw_in(j,i) + rad_lw_out(j,i)

          ENDDO
       ENDDO

       RETURN

    END SUBROUTINE radiation_clearsky

!------------------------------------------------------------------------------!
! Description:
! ------------
!-- Prescribed net radiation
!------------------------------------------------------------------------------!
    SUBROUTINE radiation_constant

       rad_net = net_radiation

    END SUBROUTINE radiation_constant

!------------------------------------------------------------------------------!
! Description:
! ------------
!-- Implementation of the RRTM radiation_scheme
!------------------------------------------------------------------------------!
    SUBROUTINE radiation_rrtm


    END SUBROUTINE radiation_rrtm

 END MODULE radiation_model_mod