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

Last change on this file since 2636 was 2233, checked in by suehring, 7 years ago

last commit documented

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