source: palm/tags/release-4.0/SOURCE/radiation_model.f90 @ 4106

Last change on this file since 4106 was 1497, checked in by maronga, 9 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 6.7 KB
Line 
1 MODULE radiation_model_mod
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: radiation_model.f90 1497 2014-12-02 17:28:07Z gronemeier $
27!
28! 1496 2014-12-02 17:25:50Z maronga
29! Initial revision
30!
31!
32! Description:
33! ------------
34! Radiation model(s), to be used e.g. with the land surface scheme
35!------------------------------------------------------------------------------!
36
37    USE arrays_3d,                                                             &
38        ONLY: pt
39
40    USE control_parameters,                                                    &
41        ONLY: phi, surface_pressure, time_since_reference_point
42
43    USE indices,                                                               &
44        ONLY:  nxlg, nxrg, nyng, nysg, nzb_s_inner
45
46    USE kinds
47
48
49    IMPLICIT NONE
50
51    INTEGER(iwp) :: i, j, k
52
53
54    INTEGER(iwp) :: day_init = 1 !: day of the year at model start
55
56    LOGICAL :: radiation = .FALSE.  !: flag parameter indicating wheather the radiation model is used
57
58    REAL(wp), PARAMETER :: SW_0 = 1368.0, &       !: solar constant 
59                           pi = 3.14159265358979323_wp, &
60                           sigma_SB  = 5.67E-8_wp !: Stefan-Boltzmann constant
61 
62    REAL(wp) :: albedo = 0.2_wp,             & !: NAMELIST alpha
63                dt_radiation = 9999999.9_wp, &
64                exn,                         & !: Exner function
65                lon = 0.0_wp,                & !: longitude in radians
66                lat = 0.0_wp,                & !: latitude in radians
67                decl_1,                      & !: declination coef. 1
68                decl_2,                      & !: declination coef. 2
69                decl_3,                      & !: declination coef. 3
70                time_utc,                    & !: current time in UTC
71                time_utc_init = 0.0_wp,      & !: UTC time at model start
72                day,                         & !: current day of the year
73                lambda = 0.0_wp,             & !: longitude in degrees
74                declination,                 & !: solar declination angle
75                hour_angle,                  & !: solar hour angle
76                time_radiation = 0.0_wp,     &
77                zenith,                      & !: solar zenith angle
78                sky_trans                      !: sky transmissivity
79
80    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
81                alpha,                       & !: surface albedo
82                Rn,                          & !: net radiation at the surface
83                LW_in,                       & !: incoming longwave radiation
84                LW_out,                      & !: outgoing longwave radiation
85                SW_in,                       & !: incoming shortwave radiation
86                SW_out                         !: outgoing shortwave radiation
87
88
89    INTERFACE init_radiation
90       MODULE PROCEDURE init_radiation
91    END INTERFACE init_radiation
92
93    INTERFACE lsm_radiation
94       MODULE PROCEDURE lsm_radiation
95    END INTERFACE lsm_radiation
96
97    SAVE
98
99    PRIVATE
100
101    PUBLIC albedo, day_init, dt_radiation, init_radiation, lambda,             &
102           lsm_radiation, Rn, radiation, SW_in, sigma_SB, time_radiation,      &
103           time_utc_init 
104
105
106
107 CONTAINS
108
109!------------------------------------------------------------------------------!
110! Description:
111! ------------
112!-- Initialization of the radiation model
113!------------------------------------------------------------------------------!
114    SUBROUTINE init_radiation
115   
116
117       IMPLICIT NONE
118
119       ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
120       ALLOCATE ( Rn(nysg:nyng,nxlg:nxrg) )
121       ALLOCATE ( LW_in(nysg:nyng,nxlg:nxrg) )
122       ALLOCATE ( LW_out(nysg:nyng,nxlg:nxrg) )
123       ALLOCATE ( SW_in(nysg:nyng,nxlg:nxrg) )
124       ALLOCATE ( SW_out(nysg:nyng,nxlg:nxrg) )
125
126       alpha = albedo
127
128!
129!--    Calculate radiation scheme constants
130       decl_1 = SIN(23.45_wp * pi / 180.0_wp)
131       decl_2 = 2.0 * pi / 365.0_wp
132       decl_3 = decl_2 * 81.0_wp
133
134!
135!--    Calculate latitude and longitude angles (lat and lon, respectively)
136       lat = phi * pi / 180.0_wp
137       lon = lambda * pi / 180.0_wp
138
139       RETURN
140
141    END SUBROUTINE init_radiation
142
143
144!------------------------------------------------------------------------------!
145! Description:
146! ------------
147!-- A simple clear sky radiation model
148!------------------------------------------------------------------------------!
149    SUBROUTINE lsm_radiation
150   
151
152       IMPLICIT NONE
153
154!
155!--    Calculate current day and time based on the initial values and simulation
156!--    time
157       day = day_init + FLOOR( (time_utc_init + time_since_reference_point)    &
158                               / 86400.0_wp )
159       time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
160
161
162!
163!--    Calculate solar declination and hour angle   
164       declination = ASIN( decl_1 * SIN(decl_2 * day - decl_3) )
165       hour_angle  = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
166
167!
168!--    Calculate zenith angle
169       zenith = SIN(lat)*SIN(declination) + COS(lat) * COS(declination)        &
170                                            * COS(hour_angle)
171       zenith = MAX(0.0_wp,zenith)
172
173!
174!--    Calculate sky transmissivity
175       sky_trans = 0.6_wp + 0.2_wp * zenith
176
177!
178!--    Calculate value of the Exner function
179       exn = (surface_pressure / 1000.0_wp )**0.286_wp
180
181!
182!--    Calculate radiation fluxes and net radiation (Rn) for each grid point
183       DO i = nxlg, nxrg
184          DO j = nysg, nyng
185
186             k = nzb_s_inner(j,i)
187             SW_in(j,i)  = SW_0 * sky_trans * zenith
188             SW_out(j,i) = - alpha(j,i) * SW_in(j,i)
189             LW_out(j,i) = - sigma_SB * (pt(k,j,i) * exn)**4
190             LW_in(j,i)  = 0.8 * sigma_SB * (pt(k+1,j,i) * exn)**4
191             Rn(j,i)     = SW_in(j,i) + SW_out(j,i) + LW_in(j,i) + LW_out(j,i)
192
193          ENDDO
194       ENDDO
195
196       RETURN
197
198    END SUBROUTINE lsm_radiation
199
200
201
202 END MODULE radiation_model_mod
Note: See TracBrowser for help on using the repository browser.