source: palm/trunk/SOURCE/calc_radiation.f90 @ 1682

Last change on this file since 1682 was 1682, checked in by knoop, 9 years ago

Code annotations made doxygen readable

  • Property svn:keywords set to Id
File size: 12.9 KB
RevLine 
[1682]1!> @file calc_radiation.f90
[1036]2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
[1310]16! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]17!--------------------------------------------------------------------------------!
18!
[484]19! Current revisions:
[1]20! -----------------
[1682]21! Code annotations made doxygen readable
[1354]22!
[1321]23! Former revisions:
24! -----------------
25! $Id: calc_radiation.f90 1682 2015-10-07 23:56:08Z knoop $
26!
[1354]27! 1353 2014-04-08 15:21:23Z heinze
28! REAL constants provided with KIND-attribute
29!
[1323]30! 1322 2014-03-20 16:38:49Z raasch
31! exponent 4.0 changed to integer
32!
[1321]33! 1320 2014-03-20 08:40:49Z raasch
[1320]34! ONLY-attribute added to USE-statements,
35! kind-parameters added to all INTEGER and REAL declaration statements,
36! kinds are defined in new module kinds,
37! revision history before 2012 removed,
38! comment fields (!:) to be used for variable explanations added to
39! all variable declaration statements
[1]40!
[1037]41! 1036 2012-10-22 13:43:42Z raasch
42! code put under GPL (PALM 3.9)
43!
[1]44! Revision 1.1  2000/04/13 14:42:45  schroeter
45! Initial revision
46!
47!
48! Description:
49! -------------
[1682]50!> Calculation of the vertical divergences of the long-wave radiation-fluxes
51!> based on the parameterization of the cloud effective emissivity
[1]52!------------------------------------------------------------------------------!
[1682]53 MODULE calc_radiation_mod
54 
[1320]55    USE kinds
56   
[1]57    PRIVATE
58    PUBLIC calc_radiation
59   
[1682]60    LOGICAL, SAVE ::  first_call = .TRUE. !<
61    REAL(wp), SAVE ::  sigma = 5.67E-08_wp   !<
[1]62
[1682]63    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  lwp_ground         !<
64    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  lwp_top            !<
65    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  blackbody_emission !<
[1]66
67    INTERFACE calc_radiation
68       MODULE PROCEDURE calc_radiation
69       MODULE PROCEDURE calc_radiation_ij
70    END INTERFACE calc_radiation
71 
72 CONTAINS
73
74
75!------------------------------------------------------------------------------!
[1682]76! Description:
77! ------------
78!> Call for all grid points
[1]79!------------------------------------------------------------------------------!
80    SUBROUTINE calc_radiation
81
[1320]82       USE arrays_3d,                                                          &
83           ONLY:  dzw, pt, ql, tend
84
85       USE cloud_parameters,                                                   &
86           ONLY:  cp, l_d_cp, pt_d_t, t_d_pt
87
88       USE control_parameters,                                                 &
89           ONLY:  rho_surface
90
91       USE indices,                                                            &
92           ONLY:  nxl, nxr, nyn, nys, nzb, nzb_2d, nzt
93
94       USE kinds
95
[1]96       USE pegrid
97
[1320]98
[1]99       IMPLICIT NONE
100
[1682]101       INTEGER(iwp) ::  i      !<
102       INTEGER(iwp) ::  j      !<
103       INTEGER(iwp) ::  k      !<
104       INTEGER(iwp) ::  k_help !<
[1]105 
[1682]106       REAL(wp) :: df_p                      !<
107       REAL(wp) :: df_m                      !<
108       REAL(wp) :: effective_emission_up_m   !<
109       REAL(wp) :: effective_emission_up_p   !<
110       REAL(wp) :: effective_emission_down_m !<
111       REAL(wp) :: effective_emission_down_p !<
112       REAL(wp) :: f_up_m                    !<
113       REAL(wp) :: f_up_p                    !<
114       REAL(wp) :: f_down_m                  !<
115       REAL(wp) :: f_down_p                  !<
116       REAL(wp) :: impinging_flux_at_top     !<
117       REAL(wp) :: temperature               !<
[1]118
119
120!
121!--    On first call, allocate temporary arrays
122       IF ( first_call )  THEN
[1320]123          ALLOCATE( blackbody_emission(nzb:nzt+1), lwp_ground(nzb:nzt+1),      &
[1]124                    lwp_top(nzb:nzt+1) )
125          first_call = .FALSE.
126       ENDIF
127
128
129       DO  i = nxl, nxr
130          DO  j = nys, nyn
131!
132!--          Compute the liquid water path (LWP) and blackbody_emission
133!--          at all vertical levels
[1353]134             lwp_ground(nzb) = 0.0_wp
[1]135             lwp_top(nzt+1)  = rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
136
137             temperature     = pt(nzb,j,i) * t_d_pt(nzb) + l_d_cp * ql(nzb,j,i)
[1322]138             blackbody_emission(nzb) = sigma * temperature**4
[1]139
140             DO  k = nzb_2d(j,i)+1, nzt
141
142                k_help = ( nzt+nzb+1 ) - k
[1320]143                lwp_ground(k)   = lwp_ground(k-1) + rho_surface * ql(k,j,i) *  &
[1]144                                  dzw(k)
145
[1320]146                lwp_top(k_help) = lwp_top(k_help+1) +                          &
[1]147                                  rho_surface * ql(k_help,j,i) * dzw(k_help)
148
149                temperature     = pt(k,j,i) * t_d_pt(k) + l_d_cp * ql(k,j,i)
[1322]150                blackbody_emission(k) = sigma * temperature**4
[1]151
152             ENDDO
153
[1320]154             lwp_ground(nzt+1) = lwp_ground(nzt) +                             &
[1]155                                 rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
156             lwp_top(nzb)      = lwp_top(nzb+1)
157
[1320]158             temperature       = pt(nzt+1,j,i) * t_d_pt(nzt+1) + l_d_cp *      &
[1]159                                 ql(nzt+1,j,i)
[1322]160             blackbody_emission(nzt+1) = sigma * temperature**4
[1]161
162!
163!--          See Chlond '92, this is just a first guess
[1353]164             impinging_flux_at_top = blackbody_emission(nzb) - 100.0_wp
[1]165
166             DO  k = nzb_2d(j,i)+1, nzt
167!
168!--             Save some computational time, but this may cause load
169!--             imbalances if ql is not distributed uniformly
[1353]170                IF ( ql(k,j,i) /= 0.0_wp ) THEN
[1]171!
172!--                Compute effective emissivities
[1353]173                   effective_emission_up_p   = 1.0_wp -                        &
174                                               EXP( -130.0_wp * lwp_ground(k+1) )
175                   effective_emission_up_m   = 1.0_wp -                        &
176                                               EXP( -130.0_wp * lwp_ground(k-1) )
177                   effective_emission_down_p = 1.0_wp -                        &
178                                               EXP( -158.0_wp * lwp_top(k+1) )
179                   effective_emission_down_m = 1.0_wp -                        &
180                                               EXP( -158.0_wp * lwp_top(k-1) ) 
[1]181
182!
183!--                Compute vertical long wave radiation fluxes
[1320]184                   f_up_p = blackbody_emission(nzb) +                          &
185                            effective_emission_up_p *                          &
[1]186                           ( blackbody_emission(k) - blackbody_emission(nzb) )
187
[1320]188                   f_up_m = blackbody_emission(nzb) +                          &
189                            effective_emission_up_m *                          &
[1]190                           ( blackbody_emission(k-1) - blackbody_emission(nzb) )
191
[1320]192                   f_down_p = impinging_flux_at_top +                          &
193                              effective_emission_down_p *                      &
[1]194                             ( blackbody_emission(k) - impinging_flux_at_top )
195
[1320]196                   f_down_m = impinging_flux_at_top +                          &
197                              effective_emission_down_m *                      &
[1]198                             ( blackbody_emission(k-1) - impinging_flux_at_top )
199
200!
201!--                Divergence of vertical long wave radiation fluxes
202                   df_p = f_up_p - f_down_p
203                   df_m = f_up_m - f_down_m
204
205!
206!--                Compute tendency term         
[1320]207                   tend(k,j,i) = tend(k,j,i) -                                 &
208                                ( pt_d_t(k) / ( rho_surface * cp ) *           &
[1]209                                  ( df_p - df_m ) / dzw(k) )
210
211                ENDIF
212
213             ENDDO
214          ENDDO
215       ENDDO
216
217    END SUBROUTINE calc_radiation
218
219
220!------------------------------------------------------------------------------!
[1682]221! Description:
222! ------------
223!> Call for grid point i,j
[1]224!------------------------------------------------------------------------------!
225    SUBROUTINE calc_radiation_ij( i, j )
226
[1320]227       USE arrays_3d,                                                          &
228           ONLY:  dzw, pt, ql, tend
229
230       USE cloud_parameters,                                                   &
231           ONLY:  cp, l_d_cp, pt_d_t, t_d_pt
232
233       USE control_parameters,                                                 &
234           ONLY:  rho_surface
235
236       USE indices,                                                            &
237           ONLY:  nzb, nzb_2d, nzt
238
239       USE kinds
240
[1]241       USE pegrid
[1320]242
[1]243   
244       IMPLICIT NONE
245
[1682]246       INTEGER(iwp) ::  i      !<
247       INTEGER(iwp) ::  j      !<
248       INTEGER(iwp) ::  k      !<
249       INTEGER(iwp) ::  k_help !<
[1]250
[1682]251       REAL(wp) :: df_p                      !<
252       REAL(wp) :: df_m                      !<
253       REAL(wp) :: effective_emission_up_m   !<
254       REAL(wp) :: effective_emission_up_p   !<
255       REAL(wp) :: effective_emission_down_m !<
256       REAL(wp) :: effective_emission_down_p !<
257       REAL(wp) :: f_up_m                    !<
258       REAL(wp) :: f_up_p                    !<
259       REAL(wp) :: f_down_m                  !<
260       REAL(wp) :: f_down_p                  !<
261       REAL(wp) :: impinging_flux_at_top     !<
262       REAL(wp) :: temperature               !<
[1]263
[1320]264       
[1]265!
266!--    On first call, allocate temporary arrays
267       IF ( first_call )  THEN
[1320]268          ALLOCATE( blackbody_emission(nzb:nzt+1), lwp_ground(nzb:nzt+1),      &
[1]269                    lwp_top(nzb:nzt+1) )
270          first_call = .FALSE.
271       ENDIF
272
273!
274!--    Compute the liquid water path (LWP) and blackbody_emission
275!--    at all vertical levels
[1353]276       lwp_ground(nzb) = 0.0_wp
[1]277       lwp_top(nzt+1)  = rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
278
279       temperature     = pt(nzb,j,i) * t_d_pt(nzb) + l_d_cp * ql(nzb,j,i)
[1322]280       blackbody_emission(nzb) = sigma * temperature**4
[1]281
282       DO  k = nzb_2d(j,i)+1, nzt
283          k_help = ( nzt+nzb+1 ) - k
284          lwp_ground(k)   = lwp_ground(k-1) + rho_surface * ql(k,j,i) * dzw(k)
285
[1320]286          lwp_top(k_help) = lwp_top(k_help+1) +                                &
[1]287                            rho_surface * ql(k_help,j,i) * dzw(k_help)
288
289          temperature     = pt(k,j,i) * t_d_pt(k) + l_d_cp * ql(k,j,i)
[1322]290          blackbody_emission(k) = sigma * temperature**4
[1]291
292       ENDDO
[1320]293       lwp_ground(nzt+1) = lwp_ground(nzt) +                                   &
[1]294                           rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
295       lwp_top(nzb)      = lwp_top(nzb+1)
296
[1320]297       temperature       = pt(nzt+1,j,i) * t_d_pt(nzt+1) + l_d_cp *            &
[1]298                           ql(nzt+1,j,i)
[1322]299       blackbody_emission(nzt+1) = sigma * temperature**4
[1]300
301!
302!--    See Chlond '92, this is just a first guess
[1353]303       impinging_flux_at_top = blackbody_emission(nzb) - 100.0_wp
[1]304
305       DO  k = nzb_2d(j,i)+1, nzt
306!
307!--       Store some computational time,
308!--       this may cause load imbalances if ql is not distributed uniformly
[1353]309          IF ( ql(k,j,i) /= 0.0_wp ) THEN
[1]310!
311!--          Compute effective emissivities
[1353]312             effective_emission_up_p   = 1.0_wp -                              &
313                                         EXP( -130.0_wp * lwp_ground(k+1) )
314             effective_emission_up_m   = 1.0_wp -                              &
315                                         EXP( -130.0_wp * lwp_ground(k-1) )
316             effective_emission_down_p = 1.0_wp -                              &
317                                         EXP( -158.0_wp * lwp_top(k+1) )
318             effective_emission_down_m = 1.0_wp -                              &
319                                         EXP( -158.0_wp * lwp_top(k-1) ) 
[1]320             
321!
322!--          Compute vertical long wave radiation fluxes
[1320]323             f_up_p = blackbody_emission(nzb) + effective_emission_up_p *      &
[1]324                     ( blackbody_emission(k) - blackbody_emission(nzb) )
325
[1320]326             f_up_m = blackbody_emission(nzb) + effective_emission_up_m *      &
[1]327                     ( blackbody_emission(k-1) - blackbody_emission(nzb) )
328
[1320]329             f_down_p = impinging_flux_at_top + effective_emission_down_p *    &
[1]330                       ( blackbody_emission(k) - impinging_flux_at_top )
331
[1320]332             f_down_m = impinging_flux_at_top + effective_emission_down_m *    &
[1]333                       ( blackbody_emission(k-1) - impinging_flux_at_top )
334
335!
336!-           Divergence of vertical long wave radiation fluxes
337             df_p = f_up_p - f_down_p
338             df_m = f_up_m - f_down_m
339
340!
341!--          Compute tendency term         
[1320]342             tend(k,j,i) = tend(k,j,i) - ( pt_d_t(k) / ( rho_surface * cp ) *  &
[1]343                                         ( df_p - df_m ) / dzw(k) )
344
345          ENDIF
346
347       ENDDO
348
349    END SUBROUTINE calc_radiation_ij
350
351 END MODULE calc_radiation_mod
Note: See TracBrowser for help on using the repository browser.