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 . ! ! Copyright 1997-2014 Leibniz Universitaet Hannover !--------------------------------------------------------------------------------! ! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: radiation_model.f90 1552 2015-03-03 14:27:15Z raasch $ ! ! 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