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

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

last commit documented

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