source: palm/trunk/SOURCE/biometeorology_mod.f90 @ 3742

Last change on this file since 3742 was 3742, checked in by dom_dwd_user, 2 years ago

biometeorology_mod.f90:
(C) Allocation of the input _av grids was moved to the "sum" section of bio_3d_data_averaging to make sure averaging is only done once!
(B) Moved call of bio_calculate_thermal_index_maps from biometeorology module to time_integration (ll 1712) to make sure averaged input is updated before calculating.

time_integration.f90:
(B) Moved call of bio_calculate_thermal_index_maps from biometeorology module to time_integration (ll 1712) to make sure averaged input is updated before calculating.

  • Property svn:keywords set to Id
  • Property svn:mergeinfo set to (toggle deleted branches)
    /palm/branches/chemistry/SOURCE/biometeorology_mod.f902047-3190,​3218-3297
    /palm/branches/resler/SOURCE/biometeorology_mod.f902023-3320,​3337-3474
    /palm/branches/salsa/SOURCE/biometeorology_mod.f902503-3581
    /palm/trunk/SOURCE/biometeorology_mod.f90mergedeligible
    /palm/branches/forwind/SOURCE/biometeorology_mod.f901564-1913
    /palm/branches/fricke/SOURCE/biometeorology_mod.f90942-977
    /palm/branches/hoffmann/SOURCE/biometeorology_mod.f90989-1052
    /palm/branches/letzel/masked_output/SOURCE/biometeorology_mod.f90296-409
    /palm/branches/palm4u/SOURCE/biometeorology_mod.f902540-2692
    /palm/branches/rans/SOURCE/biometeorology_mod.f902078-3128
    /palm/branches/suehring/biometeorology_mod.f90423-666
File size: 178.7 KB
Line 
1!> @file biometeorology_mod.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM-4U.
4!
5! PALM-4U 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.
9!
10! PALM-4U 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 2018-2019 Deutscher Wetterdienst (DWD)
18! Copyright 2018-2019 Institute of Computer Science, Academy of Sciences, Prague
19! Copyright 2018-2019 Leibniz Universitaet Hannover
20!--------------------------------------------------------------------------------!
21!
22! Current revisions:
23! ------------------
24!
25!
26! Former revisions:
27! -----------------
28! $Id: biometeorology_mod.f90 3742 2019-02-14 11:25:22Z dom_dwd_user $
29! - Allocation of the input _av grids was moved to the "sum" section of
30! bio_3d_data_averaging to make sure averaging is only done once!
31! - Moved call of bio_calculate_thermal_index_maps from biometeorology module to
32! time_integration to make sure averaged input is updated before calculating.
33!
34! 3740 2019-02-13 12:35:12Z dom_dwd_user
35! - Added safety-meassure to catch the case that 'bio_mrt_av' is stated after
36! 'bio_<index>' in the output section of the p3d file.
37!
38! 3739 2019-02-13 08:05:17Z dom_dwd_user
39! - Auto-adjusting thermal_comfort flag if not set by user, but thermal_indices
40! set as output quantities.
41! - Renamed flags "bio_<index>" to "do_calculate_<index>" for better readability
42! - Removed everything related to "time_bio_results" as this is never used.
43! - Moved humidity warning to check_data_output
44! - Fixed bug in mrt calculation introduced with my commit yesterday.
45!
46! 3735 2019-02-12 09:52:40Z dom_dwd_user
47! - Fixed auto-setting of thermal index calculation flags by output
48!  as originally proposed by resler.
49! - removed bio_pet and outher configuration variables.
50! - Updated namelist.
51!
52! 3711 2019-01-31 13:44:26Z knoop
53! Introduced interface routine bio_init_checks + small error message changes
54!
55! 3693 2019-01-23 15:20:53Z dom_dwd_user
56! Added usage of time_averaged mean radiant temperature, together with calculation,
57! grid and restart routines. General cleanup and commenting.
58!
59! 3685 2019-01-21 01:02:11Z knoop
60! Some interface calls moved to module_interface + cleanup
61!
62! 3650 2019-01-04 13:01:33Z kanani
63! Bugfixes and additions for enabling restarts with biometeorology
64!
65! 3646 2018-12-28 17:58:49Z kanani
66! Remove check for simulated_time > 0, it is not required since biometeorology
67! is only called from time_integration and not during time_integration_spinup
68!
69! 3614 2018-12-10 07:05:46Z raasch
70! unused variables removed
71!
72! 3593 2018-12-03 13:51:13Z kanani
73! Bugfix: additional tmrt_grid allocation in case bio_mrt not selected as ouput,
74! replace degree symbol by degree_C
75!
76! 3582 2018-11-29 19:16:36Z suehring
77! Consistently use bio_fill_value everywhere,
78! move allocation and initialization of output variables to bio_check_data_output
79! and bio_3d_data_averaging,
80! dom_dwd_user, Schrempf:
81! - integration of UVEM module part from r3474 (edited)
82! - split UTCI regression into 6 parts
83! - all data_output_3d is now explicity casted to sp
84!
85! 3525 2018-11-14 16:06:14Z kanani
86! Clean up, renaming from "biom" to "bio", summary of thermal index calculation
87! into one module (dom_dwd_user)
88!
89! 3479 2018-11-01 16:00:30Z gronemeier
90! - reworked check for output quantities
91! - assign new palm-error numbers
92! - set unit according to data standard.
93!
94! 3475 2018-10-30 21:16:31Z kanani
95! Add option for using MRT from RTM instead of simplified global value
96!
97! 3464 2018-10-30 18:08:55Z kanani
98! From branch resler@3462, pavelkrc:
99! make use of basic_constants_and_equations_mod
100!
101! 3448 2018-10-29 18:14:31Z kanani
102! Initial revision
103!
104!
105!
106! Authors:
107! --------
108! @author Dominik Froehlich <dominik.froehlich@dwd.de>, thermal indices
109! @author Jaroslav Resler <resler@cs.cas.cz>, mean radiant temperature
110! @author Michael Schrempf <schrempf@muk.uni-hannover.de>, uv exposure
111!
112!
113! Description:
114! ------------
115!> Biometeorology module consisting of two parts:
116!> 1.: Human thermal comfort module calculating thermal perception of a sample
117!> human being under the current meteorological conditions.
118!> 2.: Calculation of vitamin-D weighted UV exposure
119!>
120!> @todo Alphabetical sorting of "USE ..." lists, "ONLY" list, variable declarations
121!>       (per subroutine: first all CHARACTERs, then INTEGERs, LOGICALs, REALs, )
122!> @todo Comments start with capital letter --> "!-- Include..."
123!> @todo uv_vitd3dose-->new output type necessary (cumulative)
124!> @todo consider upwelling radiation in UV
125!>
126!> @note nothing now
127!>
128!> @bug  no known bugs by now
129!------------------------------------------------------------------------------!
130 MODULE biometeorology_mod
131
132    USE arrays_3d,                                                             &
133        ONLY:  pt, p, u, v, w, q
134
135    USE averaging,                                                             &
136        ONLY:  pt_av, q_av, u_av, v_av, w_av
137
138    USE basic_constants_and_equations_mod,                                     &
139        ONLY:  c_p, degc_to_k, l_v, magnus, sigma_sb, pi
140
141    USE control_parameters,                                                    &
142        ONLY:  average_count_3d, biometeorology, dz, dz_stretch_factor,        &
143               dz_stretch_level, humidity, initializing_actions, nz_do3d,      &
144               surface_pressure
145
146    USE date_and_time_mod,                                                     &
147        ONLY:  calc_date_and_time, day_of_year, time_utc
148
149    USE grid_variables,                                                        &
150        ONLY:  ddx, dx, ddy, dy
151
152    USE indices,                                                               &
153        ONLY:  nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr, nxlg, nxrg,   &
154               nysg, nyng
155
156    USE kinds  !< Set precision of INTEGER and REAL arrays according to PALM
157
158    USE netcdf_data_input_mod,                                                 &
159        ONLY:  netcdf_data_input_uvem, uvem_projarea_f, uvem_radiance_f,       &
160               uvem_irradiance_f, uvem_integration_f, building_obstruction_f
161!
162!-- Import radiation model to obtain input for mean radiant temperature
163    USE radiation_model_mod,                                                   &
164        ONLY:  ix, iy, iz, id, mrt_nlevels, mrt_include_sw,                    &
165               mrtinsw, mrtinlw, mrtbl, nmrtbl, radiation,                     &
166               radiation_interactions, rad_sw_in,                              &
167               rad_sw_out, rad_lw_in, rad_lw_out
168
169    USE surface_mod,                                                           &
170        ONLY:  get_topography_top_index_ji
171
172    IMPLICIT NONE
173
174    PRIVATE
175
176!
177!-- Declare all global variables within the module (alphabetical order)
178    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tmrt_grid  !< tmrt results (degree_C)
179    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  perct      !< PT results   (degree_C)
180    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci       !< UTCI results (degree_C)
181    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet        !< PET results  (degree_C)
182!
183!-- Grids for averaged thermal indices
184    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  mrt_av_grid   !< time average mean
185    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tmrt_av_grid  !< tmrt results (degree_C)
186    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  perct_av      !< PT results (aver. input)   (degree_C)
187    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci_av       !< UTCI results (aver. input) (degree_C)
188    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet_av        !< PET results (aver. input)  (degree_C)
189
190
191    INTEGER( iwp ) ::  bio_cell_level     !< cell level biom calculates for
192    REAL ( wp )    ::  bio_output_height  !< height output is calculated in m
193    REAL ( wp ), PARAMETER ::  human_absorb = 0.7_wp  !< SW absorbtivity of a human body (Fanger 1972)
194    REAL ( wp ), PARAMETER ::  human_emiss = 0.97_wp  !< LW emissivity of a human body after (Fanger 1972)
195    REAL ( wp ), PARAMETER ::  bio_fill_value = -9999._wp  !< set module fill value, replace by global fill value as soon as available
196!
197!--
198    LOGICAL ::  thermal_comfort  = .FALSE.  !< Enables or disables the entire thermal comfort part
199    LOGICAL ::  do_average_theta = .FALSE.  !< switch: do theta averaging in this module? (if .FALSE. this is done globally)
200    LOGICAL ::  do_average_q     = .FALSE.  !< switch: do e averaging in this module?
201    LOGICAL ::  do_average_u     = .FALSE.  !< switch: do u averaging in this module?
202    LOGICAL ::  do_average_v     = .FALSE.  !< switch: do v averaging in this module?
203    LOGICAL ::  do_average_w     = .FALSE.  !< switch: do w averaging in this module?
204    LOGICAL ::  do_average_mrt   = .FALSE.  !< switch: do mrt averaging in this module?
205    LOGICAL ::  average_trigger_perct  = .FALSE.  !< update averaged input on call to bio_perct?
206    LOGICAL ::  average_trigger_utci   = .FALSE.  !< update averaged input on call to bio_utci?
207    LOGICAL ::  average_trigger_pet    = .FALSE.  !< update averaged input on call to bio_pet?
208    LOGICAL ::  do_calculate_perct     = .FALSE.  !< Turn index PT (instant. input) on or off
209    LOGICAL ::  do_calculate_perct_av  = .FALSE.  !< Turn index PT (averaged input) on or off
210    LOGICAL ::  do_calculate_pet       = .FALSE.  !< Turn index PET (instant. input) on or off
211    LOGICAL ::  do_calculate_pet_av    = .FALSE.  !< Turn index PET (averaged input) on or off
212    LOGICAL ::  do_calculate_utci      = .FALSE.  !< Turn index UTCI (instant. input) on or off
213    LOGICAL ::  do_calculate_utci_av   = .FALSE.  !< Turn index UTCI (averaged input) on or off
214
215!
216!-- UVEM parameters from here
217!
218!-- Declare all global variables within the module (alphabetical order)
219    INTEGER(iwp) ::  bio_nmrtbl
220    INTEGER(iwp) ::  ai                      = 0  !< loop index in azimuth direction
221    INTEGER(iwp) ::  bi                      = 0  !< loop index of bit location within an 8bit-integer (one Byte)
222    INTEGER(iwp) ::  clothing                = 1  !< clothing (0=unclothed, 1=Arms,Hands,Face free, 3=Hand,Face free)
223    INTEGER(iwp) ::  iq                      = 0  !< loop index of irradiance quantity
224    INTEGER(iwp) ::  pobi                    = 0  !< loop index of the position of corresponding byte within ibset byte vektor
225    INTEGER(iwp) ::  obstruction_direct_beam = 0  !< Obstruction information for direct beam   
226    INTEGER(iwp) ::  zi                      = 0  !< loop index in zenith direction
227
228    INTEGER(KIND=1), DIMENSION(0:44)  ::  obstruction_temp1 = 0  !< temporary obstruction information stored with ibset
229    INTEGER(iwp),    DIMENSION(0:359) ::  obstruction_temp2 = 0  !< restored temporary obstruction information from ibset file
230
231    INTEGER(iwp), DIMENSION(0:35,0:9) ::  obstruction       = 1  !< final 2D obstruction information array
232
233    LOGICAL ::  consider_obstructions = .TRUE.   !< namelist parameter (see documentation)
234    LOGICAL ::  sun_in_south          = .FALSE.  !< namelist parameter (see documentation)
235    LOGICAL ::  turn_to_sun           = .TRUE.   !< namelist parameter (see documentation)
236    LOGICAL ::  uv_exposure           = .FALSE.  !< namelist parameter (see documentation)
237
238    REAL(wp) ::  diffuse_exposure            =   0.0_wp  !< calculated exposure by diffuse radiation
239    REAL(wp) ::  direct_exposure             =   0.0_wp  !< calculated exposure by direct solar beam   
240    REAL(wp) ::  orientation_angle           =   0.0_wp  !< orientation of front/face of the human model   
241    REAL(wp) ::  projection_area_direct_beam =   0.0_wp  !< projection area for direct solar beam
242    REAL(wp) ::  saa                         = 180.0_wp  !< solar azimuth angle
243    REAL(wp) ::  startpos_human              =   0.0_wp  !< start value for azimuth interpolation of human geometry array
244    REAL(wp) ::  startpos_saa_float          =   0.0_wp  !< start value for azimuth interpolation of radiance array
245    REAL(wp) ::  sza                         =  20.0_wp  !< solar zenith angle
246    REAL(wp) ::  xfactor                     =   0.0_wp  !< relative x-position used for interpolation
247    REAL(wp) ::  yfactor                     =   0.0_wp  !< relative y-position used for interpolation
248
249    REAL(wp), DIMENSION(0:2)  ::  irradiance =   0.0_wp  !< iradiance values extracted from irradiance lookup table 
250
251    REAL(wp), DIMENSION(0:2,0:90) ::  irradiance_lookup_table      = 0.0_wp  !< irradiance lookup table
252    REAL(wp), DIMENSION(0:35,0:9) ::  integration_array            = 0.0_wp  !< solid angle factors for hemispherical integration
253    REAL(wp), DIMENSION(0:35,0:9) ::  projection_area              = 0.0_wp  !< projection areas of a human (all directions)
254    REAL(wp), DIMENSION(0:35,0:9) ::  projection_area_lookup_table = 0.0_wp  !< human geometry lookup table (projection areas)
255    REAL(wp), DIMENSION(0:71,0:9) ::  projection_area_direct_temp  = 0.0_wp  !< temporary projection area for direct solar beam
256    REAL(wp), DIMENSION(0:71,0:9) ::  projection_area_temp         = 0.0_wp  !< temporary projection area for all directions
257    REAL(wp), DIMENSION(0:35,0:9) ::  radiance_array               = 0.0_wp  !< radiance extracted from radiance_lookup_table 
258    REAL(wp), DIMENSION(0:71,0:9) ::  radiance_array_temp          = 0.0_wp  !< temporary radiance data
259
260    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vitd3_exposure     !< result variable for instantaneous vitamin-D weighted exposures
261    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vitd3_exposure_av  !< result variable for summation of vitamin-D weighted exposures
262
263    REAL(wp), DIMENSION(0:35,0:9,0:90) ::  radiance_lookup_table   = 0.0_wp  !< radiance lookup table
264
265!
266!-- INTERFACES that must be available to other modules (alphabetical order)
267
268    PUBLIC bio_3d_data_averaging, bio_check_data_output,                       &
269    bio_calculate_mrt_grid, bio_calculate_thermal_index_maps, bio_calc_ipt,    &
270    bio_check_parameters, bio_data_output_3d, bio_data_output_2d,              &
271    bio_define_netcdf_grid, bio_get_thermal_index_input_ij, bio_header,        &
272    bio_init, bio_init_checks, bio_parin, thermal_comfort,                     &
273    bio_nmrtbl, bio_wrd_local, bio_rrd_local, bio_wrd_global, bio_rrd_global
274!
275!-- UVEM PUBLIC variables and methods
276    PUBLIC uvem_calc_exposure, uv_exposure
277
278!
279!-- PALM interfaces:
280!
281!-- 3D averaging for HTCM _INPUT_ variables
282    INTERFACE bio_3d_data_averaging
283       MODULE PROCEDURE bio_3d_data_averaging
284    END INTERFACE bio_3d_data_averaging
285!
286!-- Calculate mtr from rtm fluxes and assign into 2D grid
287    INTERFACE bio_calculate_mrt_grid
288       MODULE PROCEDURE bio_calculate_mrt_grid
289    END INTERFACE bio_calculate_mrt_grid
290!
291!-- Calculate static thermal indices PT, UTCI and/or PET
292    INTERFACE bio_calculate_thermal_index_maps
293       MODULE PROCEDURE bio_calculate_thermal_index_maps
294    END INTERFACE bio_calculate_thermal_index_maps
295!
296!-- Calculate the dynamic index iPT (to be caled by the agent model)
297    INTERFACE bio_calc_ipt
298       MODULE PROCEDURE bio_calc_ipt
299    END INTERFACE bio_calc_ipt
300!
301!-- Data output checks for 2D/3D data to be done in check_parameters
302    INTERFACE bio_check_data_output
303       MODULE PROCEDURE bio_check_data_output
304    END INTERFACE bio_check_data_output
305!
306!-- Input parameter checks to be done in check_parameters
307    INTERFACE bio_check_parameters
308       MODULE PROCEDURE bio_check_parameters
309    END INTERFACE bio_check_parameters
310!
311!-- Data output of 2D quantities
312    INTERFACE bio_data_output_2d
313       MODULE PROCEDURE bio_data_output_2d
314    END INTERFACE bio_data_output_2d
315!
316!-- no 3D data, thus, no averaging of 3D data, removed
317    INTERFACE bio_data_output_3d
318       MODULE PROCEDURE bio_data_output_3d
319    END INTERFACE bio_data_output_3d
320!
321!-- Definition of data output quantities
322    INTERFACE bio_define_netcdf_grid
323       MODULE PROCEDURE bio_define_netcdf_grid
324    END INTERFACE bio_define_netcdf_grid
325!
326!-- Obtains all relevant input values to estimate local thermal comfort/stress
327    INTERFACE bio_get_thermal_index_input_ij
328       MODULE PROCEDURE bio_get_thermal_index_input_ij
329    END INTERFACE bio_get_thermal_index_input_ij
330!
331!-- Output of information to the header file
332    INTERFACE bio_header
333       MODULE PROCEDURE bio_header
334    END INTERFACE bio_header
335!
336!-- Initialization actions
337    INTERFACE bio_init
338       MODULE PROCEDURE bio_init
339    END INTERFACE bio_init
340!
341!-- Initialization checks
342    INTERFACE bio_init_checks
343       MODULE PROCEDURE bio_init_checks
344    END INTERFACE bio_init_checks
345!
346!-- Reading of NAMELIST parameters
347    INTERFACE bio_parin
348       MODULE PROCEDURE bio_parin
349    END INTERFACE bio_parin
350!
351!-- Read global restart parameters
352    INTERFACE bio_rrd_global
353       MODULE PROCEDURE bio_rrd_global
354    END INTERFACE bio_rrd_global
355!
356!-- Read local restart parameters
357    INTERFACE bio_rrd_local
358       MODULE PROCEDURE bio_rrd_local
359    END INTERFACE bio_rrd_local
360!
361!-- Write global restart parameters
362    INTERFACE bio_wrd_global
363       MODULE PROCEDURE bio_wrd_global
364    END INTERFACE bio_wrd_global
365!
366!-- Write local restart parameters
367    INTERFACE bio_wrd_local
368       MODULE PROCEDURE bio_wrd_local
369    END INTERFACE bio_wrd_local
370!
371!-- Calculate UV exposure grid
372    INTERFACE uvem_calc_exposure
373       MODULE PROCEDURE uvem_calc_exposure
374    END INTERFACE uvem_calc_exposure
375
376 CONTAINS
377
378
379!------------------------------------------------------------------------------!
380! Description:
381! ------------
382!> Sum up and time-average biom input quantities as well as allocate
383!> the array necessary for storing the average.
384!> There is a considerable difference to the 3d_data_averaging subroutines
385!> used by other modules:
386!> For the thermal indices, the module needs to average the input conditions
387!> not the result!
388!------------------------------------------------------------------------------!
389 SUBROUTINE bio_3d_data_averaging( mode, variable )
390
391    IMPLICIT NONE
392
393    CHARACTER (LEN=*) ::  mode     !< averaging mode: allocate, sum, or average
394    CHARACTER (LEN=*) ::  variable !< The variable in question
395
396    INTEGER(iwp) ::  i        !< Running index, x-dir
397    INTEGER(iwp) ::  j        !< Running index, y-dir
398    INTEGER(iwp) ::  k        !< Running index, z-dir
399
400
401    IF ( mode == 'allocate' )  THEN
402
403       SELECT CASE ( TRIM( variable ) )
404
405          CASE ( 'bio_mrt' )
406
407                IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
408                   ALLOCATE( mrt_av_grid(nmrtbl) )
409                ENDIF
410                mrt_av_grid = 0.0_wp
411                do_average_mrt = .FALSE.  !< overwrite if that was enabled somehow
412
413
414          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
415
416!
417!--          Averaging, as well as the allocation of the required grids must be
418!--          done only once, independent from for how many thermal indices
419!--          averaged output is desired.
420!--          Therefore wee need to memorize which index is the one that controls
421!--          the averaging (what must be the first thermal index called).
422!--          Indices are in unknown order as depending on the input file,
423!--          determine first index to average und update only once
424!
425!--          Only proceed here if this was not done for any index before. This
426!--          is done only once during the whole model run.
427             IF ( .NOT. average_trigger_perct  .AND.                           &
428                  .NOT. average_trigger_utci   .AND.                           &
429                  .NOT. average_trigger_pet )  THEN
430!
431!--             Memorize the first index called to control averaging
432                IF ( TRIM( variable ) == 'bio_perct*' )  THEN
433                    average_trigger_perct = .TRUE.
434                ENDIF
435                IF ( TRIM( variable ) == 'bio_utci*' )  THEN
436                    average_trigger_utci = .TRUE.
437                ENDIF
438                IF ( TRIM( variable ) == 'bio_pet*' )  THEN
439                    average_trigger_pet = .TRUE.
440                ENDIF
441             ENDIF
442!
443!--          Allocation of the input _av grids was moved to the "sum" section to
444!--          make sure averaging is only done once!
445
446
447          CASE ( 'uvem_vitd3dose*' )
448             IF ( .NOT. ALLOCATED( vitd3_exposure_av ) )  THEN
449                ALLOCATE( vitd3_exposure_av(nysg:nyng,nxlg:nxrg) )
450             ENDIF
451             vitd3_exposure_av = 0.0_wp
452
453          CASE DEFAULT
454             CONTINUE
455
456       END SELECT
457
458    ELSEIF ( mode == 'sum' )  THEN
459
460       SELECT CASE ( TRIM( variable ) )
461
462          CASE ( 'bio_mrt' )
463!
464!--          Consider the case 'bio_mrt' is called after some thermal index. In
465!--          that case do_average_mrt will be .TRUE. leading to a double-
466!--          averaging.
467             IF ( .NOT. do_average_mrt  .AND.  ALLOCATED( mrt_av_grid ) )  THEN
468
469                IF ( mrt_include_sw )  THEN
470                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
471                      ( ( human_absorb * mrtinsw(:) +                          &
472                      mrtinlw(:) ) /                                           &
473                      ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
474                ELSE
475                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
476                      ( mrtinlw(:) /                                           &
477                      ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
478                ENDIF
479             ENDIF
480
481          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
482!
483!--          Only continue if the current index is the one to trigger the input
484!--          averaging, see above
485             IF ( average_trigger_perct  .AND.  TRIM( variable ) /=            &
486                'bio_perct*')  RETURN
487             IF ( average_trigger_utci   .AND.  TRIM( variable ) /=            &
488                'bio_utci*')   RETURN
489             IF ( average_trigger_pet    .AND.  TRIM( variable ) /=            &
490                'bio_pet*')    RETURN
491!
492!--          Now memorize which of the input grids are not averaged by other
493!--          modules. Set averaging switch to .TRUE. and allocate the respective
494!--          grid in that case.
495             IF ( .NOT. ALLOCATED( pt_av ) )  THEN  !< if not averaged by other module
496                ALLOCATE( pt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
497                do_average_theta = .TRUE.  !< memorize, that bio is responsible
498                pt_av = 0.0_wp
499             ENDIF
500             IF ( ALLOCATED( pt_av )  .AND.  do_average_theta )  THEN
501                DO  i = nxl, nxr
502                   DO  j = nys, nyn
503                      DO  k = nzb, nzt+1
504                         pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i)
505                      ENDDO
506                   ENDDO
507                ENDDO
508             ENDIF
509
510             IF ( .NOT. ALLOCATED( q_av ) )  THEN
511                ALLOCATE( q_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
512                do_average_q = .TRUE.
513                q_av = 0.0_wp
514             ENDIF
515             IF ( ALLOCATED( q_av )  .AND.  do_average_q )  THEN
516                DO  i = nxl, nxr
517                   DO  j = nys, nyn
518                      DO  k = nzb, nzt+1
519                         q_av(k,j,i) = q_av(k,j,i) + q(k,j,i)
520                      ENDDO
521                   ENDDO
522                ENDDO
523             ENDIF
524
525             IF ( .NOT. ALLOCATED( u_av ) )  THEN
526                ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
527                do_average_u = .TRUE.
528                u_av = 0.0_wp
529             ENDIF
530             IF ( ALLOCATED( u_av )  .AND.  do_average_u )  THEN
531                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
532                   DO  j = nysg, nyng
533                      DO  k = nzb, nzt+1
534                         u_av(k,j,i) = u_av(k,j,i) + u(k,j,i)
535                      ENDDO
536                   ENDDO
537                ENDDO
538             ENDIF
539
540             IF ( .NOT. ALLOCATED( v_av ) )  THEN
541                ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
542                do_average_v = .TRUE.
543                v_av = 0.0_wp
544             ENDIF
545             IF ( ALLOCATED( v_av )  .AND.  do_average_v )  THEN
546                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
547                   DO  j = nysg, nyng
548                      DO  k = nzb, nzt+1
549                         v_av(k,j,i) = v_av(k,j,i) + v(k,j,i)
550                      ENDDO
551                   ENDDO
552                ENDDO
553             ENDIF
554
555             IF ( .NOT. ALLOCATED( w_av ) )  THEN
556                ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
557                do_average_w = .TRUE.
558                w_av = 0.0_wp
559             ENDIF
560             IF ( ALLOCATED( w_av )  .AND.  do_average_w )  THEN
561                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
562                   DO  j = nysg, nyng
563                      DO  k = nzb, nzt+1
564                         w_av(k,j,i) = w_av(k,j,i) + w(k,j,i)
565                      ENDDO
566                   ENDDO
567                ENDDO
568             ENDIF
569
570             IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
571                ALLOCATE( mrt_av_grid(nmrtbl) )
572                do_average_mrt = .TRUE.
573                mrt_av_grid = 0.0_wp
574             ENDIF
575             IF ( ALLOCATED( mrt_av_grid )  .AND.  do_average_mrt )  THEN
576
577                IF ( mrt_include_sw )  THEN
578                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
579                      ( ( human_absorb * mrtinsw(:) +                          &
580                      mrtinlw(:) ) /                                           &
581                      ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
582                ELSE
583                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
584                      ( mrtinlw(:) /                                           &
585                      ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
586                ENDIF
587             ENDIF
588!
589!--       This is a cumulated dose. No mode == 'average' for this quantity.
590          CASE ( 'uvem_vitd3dose*' )
591             IF ( ALLOCATED( vitd3_exposure_av ) )  THEN
592                DO  i = nxlg, nxrg
593                   DO  j = nysg, nyng
594                      vitd3_exposure_av(j,i) = vitd3_exposure_av(j,i) + vitd3_exposure(j,i)
595                   ENDDO
596                ENDDO
597             ENDIF
598
599          CASE DEFAULT
600             CONTINUE
601
602       END SELECT
603
604    ELSEIF ( mode == 'average' )  THEN
605
606       SELECT CASE ( TRIM( variable ) )
607
608          CASE ( 'bio_mrt' )
609!
610!--          Consider the case 'bio_mrt' is called after some thermal index. In
611!--          that case do_average_mrt will be .TRUE. leading to a double-
612!--          averaging.
613             IF ( .NOT. do_average_mrt  .AND.  ALLOCATED( mrt_av_grid ) )  THEN
614                mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d, KIND=wp )
615             ENDIF
616
617          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
618!
619!--          Only continue if update index, see above
620             IF ( average_trigger_perct  .AND.                                 &
621                TRIM( variable ) /= 'bio_perct*' )  RETURN
622             IF ( average_trigger_utci  .AND.                                  &
623                TRIM( variable ) /= 'bio_utci*' )  RETURN
624             IF ( average_trigger_pet   .AND.                                  &
625                TRIM( variable ) /= 'bio_pet*' )  RETURN
626
627             IF ( ALLOCATED( pt_av )  .AND.  do_average_theta )  THEN
628                DO  i = nxl, nxr
629                   DO  j = nys, nyn
630                      DO  k = nzb, nzt+1
631                         pt_av(k,j,i) = pt_av(k,j,i) /                         &
632                            REAL( average_count_3d, KIND=wp )
633                      ENDDO
634                   ENDDO
635                ENDDO
636             ENDIF
637
638             IF ( ALLOCATED( q_av )  .AND.  do_average_q )  THEN
639                DO  i = nxl, nxr
640                   DO  j = nys, nyn
641                      DO  k = nzb, nzt+1
642                         q_av(k,j,i) = q_av(k,j,i) /                           &
643                            REAL( average_count_3d, KIND=wp )
644                      ENDDO
645                   ENDDO
646                ENDDO
647             ENDIF
648
649             IF ( ALLOCATED( u_av )  .AND.  do_average_u )  THEN
650                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
651                   DO  j = nysg, nyng
652                      DO  k = nzb, nzt+1
653                         u_av(k,j,i) = u_av(k,j,i) /                           &
654                            REAL( average_count_3d, KIND=wp )
655                      ENDDO
656                   ENDDO
657                ENDDO
658             ENDIF
659
660             IF ( ALLOCATED( v_av )  .AND.  do_average_v )  THEN
661                DO  i = nxlg, nxrg
662                   DO  j = nysg, nyng
663                      DO  k = nzb, nzt+1
664                         v_av(k,j,i) = v_av(k,j,i) /                           &
665                            REAL( average_count_3d, KIND=wp )
666                      ENDDO
667                   ENDDO
668                ENDDO
669             ENDIF
670
671             IF ( ALLOCATED( w_av )  .AND.  do_average_w )  THEN
672                DO  i = nxlg, nxrg
673                   DO  j = nysg, nyng
674                      DO  k = nzb, nzt+1
675                         w_av(k,j,i) = w_av(k,j,i) /                           &
676                            REAL( average_count_3d, KIND=wp )
677                      ENDDO
678                   ENDDO
679                ENDDO
680             ENDIF
681
682             IF ( ALLOCATED( mrt_av_grid )  .AND.  do_average_mrt )  THEN
683                mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d,      &
684                KIND=wp )
685             ENDIF
686
687!
688!--     No averaging for UVEM since we are calculating a dose (only sum is
689!--     calculated and saved to av.nc file)
690
691        END SELECT
692
693    ENDIF
694
695
696 END SUBROUTINE bio_3d_data_averaging
697
698
699
700!------------------------------------------------------------------------------!
701! Description:
702! ------------
703!> Check data output for biometeorology model
704!------------------------------------------------------------------------------!
705 SUBROUTINE bio_check_data_output( var, unit, i, j, ilen, k )
706
707    USE control_parameters,                                                    &
708        ONLY: data_output, message_string
709
710    IMPLICIT NONE
711
712    CHARACTER (LEN=*) ::  unit     !< The unit for the variable var
713    CHARACTER (LEN=*) ::  var      !< The variable in question
714
715    INTEGER(iwp), INTENT(IN) ::  i     !< Current element of data_output
716    INTEGER(iwp), INTENT(IN) ::  j     !< Average quantity? 0 = no, 1 = yes
717    INTEGER(iwp), INTENT(IN) ::  ilen  !< Length of current entry in data_output
718    INTEGER(iwp), INTENT(IN) ::  k     !< Output is xy mode? 0 = no, 1 = yes
719
720    SELECT CASE ( TRIM( var ) )
721!
722!--    Allocate a temporary array with the desired output dimensions.
723!--    Arrays for time-averaged thermal indices are also allocated here because
724!--    they are not running through the standard averaging procedure in
725!--    bio_3d_data_averaging as the values of the averaged thermal indices are
726!--    derived in a single step based on priorly averaged arrays (see
727!--    bio_calculate_thermal_index_maps).
728       CASE ( 'bio_mrt' )
729          unit = 'degree_C'
730          thermal_comfort = .TRUE.  !< enable thermal_comfort if user forgot to do so
731          IF ( .NOT. ALLOCATED( tmrt_grid ) )  THEN
732             ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) )
733             tmrt_grid = REAL( bio_fill_value, KIND = wp )
734          ENDIF
735
736       CASE ( 'bio_perct*' )
737          unit = 'degree_C'
738          thermal_comfort = .TRUE.
739          IF ( j == 0 )  THEN                !< if instantaneous input
740             do_calculate_perct = .TRUE.
741             IF ( .NOT. ALLOCATED( perct ) )  THEN
742                ALLOCATE( perct (nys:nyn,nxl:nxr) )
743                perct = REAL( bio_fill_value, KIND = wp )
744             ENDIF
745          ELSE                              !< if averaged input
746             do_calculate_perct_av = .TRUE.
747             IF ( .NOT. ALLOCATED( perct_av ) )  THEN
748                ALLOCATE( perct_av (nys:nyn,nxl:nxr) )
749                perct_av = REAL( bio_fill_value, KIND = wp )
750             ENDIF
751          ENDIF
752
753       CASE ( 'bio_utci*' )
754          unit = 'degree_C'
755          thermal_comfort = .TRUE.
756          IF ( j == 0 )  THEN
757             do_calculate_utci = .TRUE.
758             IF ( .NOT. ALLOCATED( utci ) )  THEN
759                ALLOCATE( utci (nys:nyn,nxl:nxr) )
760                utci = REAL( bio_fill_value, KIND = wp )
761             ENDIF
762          ELSE
763             do_calculate_utci_av = .TRUE.
764             IF ( .NOT. ALLOCATED( utci_av ) )  THEN
765                ALLOCATE( utci_av (nys:nyn,nxl:nxr) )
766                utci_av = REAL( bio_fill_value, KIND = wp )
767             ENDIF
768          ENDIF
769
770       CASE ( 'bio_pet*' )
771          unit = 'degree_C'
772          thermal_comfort = .TRUE.
773          IF ( j == 0 )  THEN
774             do_calculate_pet = .TRUE.
775             IF ( .NOT. ALLOCATED( pet ) )  THEN
776                ALLOCATE( pet (nys:nyn,nxl:nxr) )
777                pet = REAL( bio_fill_value, KIND = wp )
778             ENDIF
779          ELSE
780             do_calculate_pet_av = .TRUE.
781             IF ( .NOT. ALLOCATED( pet_av ) )  THEN
782                ALLOCATE( pet_av (nys:nyn,nxl:nxr) )
783                pet_av = REAL( bio_fill_value, KIND = wp )
784             ENDIF
785          ENDIF
786
787
788       CASE ( 'uvem_vitd3*' )
789          IF (  .NOT.  uv_exposure )  THEN
790             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
791                      'res a namelist &uvexposure_par'
792             CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 )
793          ENDIF
794          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
795             message_string = 'illegal value for data_output: "' //            &
796                              TRIM( var ) // '" & only 2d-horizontal ' //      &
797                              'cross sections are allowed for this value'
798             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
799          ENDIF
800          unit = 'IU/s'
801          IF ( .NOT. ALLOCATED( vitd3_exposure ) )  THEN
802             ALLOCATE( vitd3_exposure(nysg:nyng,nxlg:nxrg) )
803          ENDIF
804          vitd3_exposure = 0.0_wp
805
806       CASE ( 'uvem_vitd3dose*' )
807          IF (  .NOT.  uv_exposure )  THEN
808             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
809                      'res  a namelist &uvexposure_par'
810             CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 )
811          ENDIF
812          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
813             message_string = 'illegal value for data_output: "' //            &
814                              TRIM( var ) // '" & only 2d-horizontal ' //      &
815                              'cross sections are allowed for this value'
816             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
817          ENDIF
818          unit = 'IU/av-h'
819          IF ( .NOT. ALLOCATED( vitd3_exposure_av ) )  THEN
820             ALLOCATE( vitd3_exposure_av(nysg:nyng,nxlg:nxrg) )
821          ENDIF
822          vitd3_exposure_av = 0.0_wp
823
824       CASE DEFAULT
825          unit = 'illegal'
826
827    END SELECT
828
829!
830!--    Further checks if thermal comfort output is desired.
831    IF ( thermal_comfort  .AND.  unit == 'degree_C' )  THEN
832
833!
834!--    Break if required modules "radiation" and "humidity" are not avalable.
835       IF ( .NOT.  radiation )  THEN
836          message_string = 'output of "' // TRIM( var ) // '" require'         &
837                           // 's radiation = .TRUE.'
838          CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
839          unit = 'illegal'
840       ENDIF
841       IF ( .NOT.  humidity )  THEN
842          message_string = 'The estimation of thermal comfort requires '    // &
843                           'air humidity information, but humidity module ' // &
844                           'is disabled!'
845          CALL message( 'check_parameters', 'PA0561', 1, 2, 0, 6, 0 )
846          unit = 'illegal'
847       ENDIF
848       IF ( mrt_nlevels == 0 )  THEN
849          message_string = 'output of "' // TRIM( var ) // '" require'         &
850                           // 's mrt_nlevels > 0'
851          CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
852          unit = 'illegal'
853       ENDIF
854
855
856    ENDIF
857
858 END SUBROUTINE bio_check_data_output
859
860!------------------------------------------------------------------------------!
861! Description:
862! ------------
863!> Check parameters routine for biom module
864!> Currently unused but might come in handy for future checks?
865!------------------------------------------------------------------------------!
866 SUBROUTINE bio_check_parameters
867
868    USE control_parameters,                                                    &
869        ONLY:  message_string
870
871    IMPLICIT NONE
872
873
874
875 END SUBROUTINE bio_check_parameters
876
877
878!------------------------------------------------------------------------------!
879! Description:
880! ------------
881!> Subroutine defining 2D output variables
882!> data_output_2d 1188ff
883!------------------------------------------------------------------------------!
884 SUBROUTINE bio_data_output_2d( av, variable, found, grid, local_pf,           &
885                                two_d, nzb_do, nzt_do )
886
887
888    USE kinds
889
890
891    IMPLICIT NONE
892!
893!-- Input variables
894    CHARACTER (LEN=*), INTENT(IN) ::  variable    !< Char identifier to select var for output
895    INTEGER(iwp), INTENT(IN)      ::  av          !< Use averaged data? 0 = no, 1 = yes?
896    INTEGER(iwp), INTENT(IN)      ::  nzb_do      !< Unused. 2D. nz bottom to nz top
897    INTEGER(iwp), INTENT(IN)      ::  nzt_do      !< Unused.
898!
899!-- Output variables
900    CHARACTER (LEN=*), INTENT(OUT) ::  grid   !< Grid type (always "zu1" for biom)
901    LOGICAL, INTENT(OUT)           ::  found  !< Output found?
902    LOGICAL, INTENT(OUT)           ::  two_d  !< Flag parameter that indicates 2D variables,
903                                              !< horizontal cross sections, must be .TRUE. for thermal indices and uv
904    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !< Temp. result grid to return
905!
906!-- Internal variables
907    INTEGER(iwp) ::  i        !< Running index, x-dir
908    INTEGER(iwp) ::  j        !< Running index, y-dir
909    INTEGER(iwp) ::  k        !< Running index, z-dir
910    INTEGER(iwp) ::  l        !< Running index, radiation grid
911
912
913    found = .TRUE.
914    local_pf = bio_fill_value
915
916    SELECT CASE ( TRIM( variable ) )
917
918
919        CASE ( 'bio_mrt_xy' )
920           grid = 'zu1'
921           two_d = .FALSE.  !< can be calculated for several levels
922           local_pf = REAL( bio_fill_value, KIND = wp )
923           DO  l = 1, nmrtbl
924              i = mrtbl(ix,l)
925              j = mrtbl(iy,l)
926              k = mrtbl(iz,l)
927              IF ( k < nzb_do  .OR.  k > nzt_do  .OR.  j < nys  .OR.           &
928                 j > nyn  .OR.  i < nxl  .OR.  i > nxr )  CYCLE
929              IF ( av == 0 )  THEN
930                 IF ( mrt_include_sw )  THEN
931                    local_pf(i,j,k) = ( ( human_absorb * mrtinsw(l) +          &
932                                    mrtinlw(l) ) /                             &
933                                    ( human_emiss * sigma_sb ) )**.25_wp -     &
934                                    degc_to_k
935                 ELSE
936                    local_pf(i,j,k) = ( mrtinlw(l)  /                          &
937                                    ( human_emiss * sigma_sb ) )**.25_wp -     &
938                                    degc_to_k
939                 ENDIF
940              ELSE
941                 local_pf(i,j,k) = mrt_av_grid(l)
942              ENDIF
943           ENDDO
944
945
946        CASE ( 'bio_perct*_xy' )        ! 2d-array
947           grid = 'zu1'
948           two_d = .TRUE.
949           IF ( av == 0 )  THEN
950              DO  i = nxl, nxr
951                 DO  j = nys, nyn
952                    local_pf(i,j,nzb+1) = perct(j,i)
953                 ENDDO
954              ENDDO
955           ELSE
956              DO  i = nxl, nxr
957                 DO  j = nys, nyn
958                    local_pf(i,j,nzb+1) = perct_av(j,i)
959                 ENDDO
960              ENDDO
961           ENDIF
962
963
964        CASE ( 'bio_utci*_xy' )        ! 2d-array
965           grid = 'zu1'
966           two_d = .TRUE.
967           IF ( av == 0 )  THEN
968              DO  i = nxl, nxr
969                 DO  j = nys, nyn
970                    local_pf(i,j,nzb+1) = utci(j,i)
971                 ENDDO
972              ENDDO
973           ELSE
974              DO  i = nxl, nxr
975                 DO  j = nys, nyn
976                    local_pf(i,j,nzb+1) = utci_av(j,i)
977                 ENDDO
978              ENDDO
979           ENDIF
980
981
982        CASE ( 'bio_pet*_xy' )        ! 2d-array
983           grid = 'zu1'
984           two_d = .TRUE.
985           IF ( av == 0 )  THEN
986              DO  i = nxl, nxr
987                 DO  j = nys, nyn
988                    local_pf(i,j,nzb+1) = pet(j,i)
989                 ENDDO
990              ENDDO
991           ELSE
992              DO  i = nxl, nxr
993                 DO  j = nys, nyn
994                    local_pf(i,j,nzb+1) = pet_av(j,i)
995                 ENDDO
996              ENDDO
997           ENDIF
998
999!
1000!--    Before data is transfered to local_pf, transfer is it 2D dummy variable and exchange ghost points therein.
1001!--    However, at this point this is only required for instantaneous arrays, time-averaged quantities are already exchanged.
1002       CASE ( 'uvem_vitd3*_xy' )        ! 2d-array
1003          IF ( av == 0 )  THEN
1004             DO  i = nxl, nxr
1005                DO  j = nys, nyn
1006                   local_pf(i,j,nzb+1) = vitd3_exposure(j,i)
1007                ENDDO
1008             ENDDO
1009          ENDIF
1010
1011          two_d = .TRUE.
1012          grid = 'zu1'
1013
1014       CASE ( 'uvem_vitd3dose*_xy' )        ! 2d-array
1015          IF ( av == 1 )  THEN
1016             DO  i = nxl, nxr
1017                DO  j = nys, nyn
1018                   local_pf(i,j,nzb+1) = vitd3_exposure_av(j,i)
1019                ENDDO
1020             ENDDO
1021          ENDIF
1022
1023          two_d = .TRUE.
1024          grid = 'zu1'
1025
1026
1027       CASE DEFAULT
1028          found = .FALSE.
1029          grid  = 'none'
1030
1031    END SELECT
1032
1033
1034 END SUBROUTINE bio_data_output_2d
1035
1036
1037!------------------------------------------------------------------------------!
1038! Description:
1039! ------------
1040!> Subroutine defining 3D output variables (dummy, always 2d!)
1041!> data_output_3d 709ff
1042!------------------------------------------------------------------------------!
1043 SUBROUTINE bio_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
1044
1045    USE indices
1046
1047    USE kinds
1048
1049
1050    IMPLICIT NONE
1051!
1052!-- Input variables
1053    CHARACTER (LEN=*), INTENT(IN) ::  variable   !< Char identifier to select var for output
1054    INTEGER(iwp), INTENT(IN) ::  av       !< Use averaged data? 0 = no, 1 = yes?
1055    INTEGER(iwp), INTENT(IN) ::  nzb_do   !< Unused. 2D. nz bottom to nz top
1056    INTEGER(iwp), INTENT(IN) ::  nzt_do   !< Unused.
1057!
1058!-- Output variables
1059    LOGICAL, INTENT(OUT) ::  found   !< Output found?
1060    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< Temp. result grid to return
1061!
1062!-- Internal variables
1063    INTEGER(iwp) ::  l    !< Running index, radiation grid
1064    INTEGER(iwp) ::  i    !< Running index, x-dir
1065    INTEGER(iwp) ::  j    !< Running index, y-dir
1066    INTEGER(iwp) ::  k    !< Running index, z-dir
1067
1068!     REAL(wp) ::  mrt  !< Buffer for mean radiant temperature
1069
1070    found = .TRUE.
1071
1072    SELECT CASE ( TRIM( variable ) )
1073
1074        CASE ( 'bio_mrt' )
1075            local_pf = REAL( bio_fill_value, KIND = sp )
1076            DO  l = 1, nmrtbl
1077               i = mrtbl(ix,l)
1078               j = mrtbl(iy,l)
1079               k = mrtbl(iz,l)
1080               IF ( k < nzb_do  .OR.  k > nzt_do  .OR.  j < nys  .OR.          &
1081                  j > nyn  .OR.  i < nxl  .OR.  i > nxr )  CYCLE
1082               IF ( av == 0 )  THEN
1083                  IF ( mrt_include_sw )  THEN
1084                     local_pf(i,j,k) = REAL( ( ( human_absorb * mrtinsw(l) +   &
1085                                    mrtinlw(l) ) /                             &
1086                                    ( human_emiss * sigma_sb ) )**.25_wp -     &
1087                                    degc_to_k, KIND = sp )
1088                  ELSE
1089                     local_pf(i,j,k) = REAL( ( mrtinlw(l) /                    &
1090                                    ( human_emiss * sigma_sb ) )**.25_wp -     &
1091                                    degc_to_k, KIND = sp )
1092                  ENDIF
1093               ELSE
1094                  local_pf(i,j,k) = REAL( mrt_av_grid(l), KIND = sp )
1095               ENDIF
1096            ENDDO
1097
1098       CASE DEFAULT
1099          found = .FALSE.
1100
1101    END SELECT
1102
1103 END SUBROUTINE bio_data_output_3d
1104
1105!------------------------------------------------------------------------------!
1106! Description:
1107! ------------
1108!> Subroutine defining appropriate grid for netcdf variables.
1109!> It is called out from subroutine netcdf_interface_mod.
1110!> netcdf_interface_mod 918ff
1111!------------------------------------------------------------------------------!
1112 SUBROUTINE bio_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
1113
1114    IMPLICIT NONE
1115!
1116!-- Input variables
1117    CHARACTER (LEN=*), INTENT(IN)  ::  var      !< Name of output variable
1118!
1119!-- Output variables
1120    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x   !< x grid of output variable
1121    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y   !< y grid of output variable
1122    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z   !< z grid of output variable
1123
1124    LOGICAL, INTENT(OUT)           ::  found    !< Flag if output var is found
1125!
1126!-- Local variables
1127    LOGICAL      :: is2d  !< Var is 2d?
1128
1129    INTEGER(iwp) :: l     !< Length of the var array
1130
1131
1132    found  = .FALSE.
1133    grid_x = 'none'
1134    grid_y = 'none'
1135    grid_z = 'none'
1136
1137    l = MAX( 2, LEN_TRIM( var ) )
1138    is2d = ( var(l-1:l) == 'xy' )
1139
1140    IF ( var(1:4) == 'bio_' )  THEN
1141       found  = .TRUE.
1142       grid_x = 'x'
1143       grid_y = 'y'
1144       grid_z = 'zu'
1145       IF ( is2d  .AND.  var(1:7) /= 'bio_mrt' )  grid_z = 'zu1'
1146    ENDIF
1147
1148    IF ( is2d  .AND.  var(1:4) == 'uvem' )  THEN
1149       grid_x = 'x'
1150       grid_y = 'y'
1151       grid_z = 'zu1'
1152    ENDIF
1153
1154 END SUBROUTINE bio_define_netcdf_grid
1155
1156!------------------------------------------------------------------------------!
1157! Description:
1158! ------------
1159!> Header output for biom module
1160!> header 982
1161!------------------------------------------------------------------------------!
1162 SUBROUTINE bio_header( io )
1163
1164    IMPLICIT NONE
1165!
1166!-- Input variables
1167    INTEGER(iwp), INTENT(IN) ::  io           !< Unit of the output file
1168!
1169!-- Internal variables
1170    CHARACTER (LEN=86) ::  output_height_chr  !< String for output height
1171
1172    WRITE( output_height_chr, '(F8.1,7X)' )  bio_output_height
1173!
1174!-- Write biom header
1175    WRITE( io, 1 )
1176    WRITE( io, 2 )  TRIM( output_height_chr )
1177    WRITE( io, 3 )  TRIM( ACHAR( bio_cell_level ) )
1178
11791   FORMAT (//' Human thermal comfort module information:'/                    &
1180              ' ------------------------------'/)
11812   FORMAT ('    --> All indices calculated for a height of (m): ', A )
11823   FORMAT ('    --> This corresponds to cell level : ', A )
1183
1184 END SUBROUTINE bio_header
1185
1186
1187!------------------------------------------------------------------------------!
1188! Description:
1189! ------------
1190!> Initialization of the HTCM
1191!> init_3d_model 1987ff
1192!------------------------------------------------------------------------------!
1193 SUBROUTINE bio_init
1194
1195    USE netcdf_data_input_mod,                                                 &
1196        ONLY:  netcdf_data_input_uvem
1197
1198    IMPLICIT NONE
1199!
1200!-- Internal vriables
1201    REAL ( wp )  :: height  !< current height in meters
1202
1203    CALL location_message( 'initializing biometeorology module', .FALSE. )
1204!
1205!-- Determine cell level corresponding to 1.1 m above ground level
1206!   (gravimetric center of sample human)
1207
1208    bio_cell_level = 0_iwp
1209    bio_output_height = 0.5_wp * dz(1)
1210    height = 0.0_wp
1211
1212    bio_cell_level = INT ( 1.099_wp / dz(1) )
1213    bio_output_height = bio_output_height + bio_cell_level * dz(1)
1214
1215!
1216!-- Init UVEM and load lookup tables
1217    IF ( uv_exposure )  CALL netcdf_data_input_uvem
1218
1219    CALL location_message( 'finished', .TRUE. )
1220
1221 END SUBROUTINE bio_init
1222
1223
1224!------------------------------------------------------------------------------!
1225! Description:
1226! ------------
1227!> Checks done after the Initialization
1228!------------------------------------------------------------------------------!
1229 SUBROUTINE bio_init_checks
1230
1231    USE control_parameters,                                                    &
1232        ONLY: message_string
1233
1234    IF ( .NOT. radiation_interactions )  THEN
1235       message_string = 'The mrt calculation requires ' //                     &
1236                        'enabled radiation_interactions but it ' //            &
1237                        'is disabled!'
1238       CALL message( 'bio_init_checks', 'PAHU03', 1, 2, 0, 6, 0 )
1239    ENDIF
1240
1241
1242 END SUBROUTINE bio_init_checks
1243
1244
1245!------------------------------------------------------------------------------!
1246! Description:
1247! ------------
1248!> Parin for &biometeorology_parameters for reading biomet parameters
1249!------------------------------------------------------------------------------!
1250 SUBROUTINE bio_parin
1251
1252    IMPLICIT NONE
1253
1254!
1255!-- Internal variables
1256    CHARACTER (LEN=80) ::  line  !< Dummy string for current line in parameter file
1257
1258    NAMELIST /biometeorology_parameters/  thermal_comfort,                     &
1259
1260!
1261!-- UVEM namelist parameters
1262                                          clothing,                            &
1263                                          consider_obstructions,               &
1264                                          orientation_angle,                   &
1265                                          sun_in_south,                        &
1266                                          turn_to_sun,                         &
1267                                          uv_exposure
1268
1269
1270!-- Try to find biometeorology_parameters namelist
1271    REWIND ( 11 )
1272    line = ' '
1273    DO WHILE ( INDEX( line, '&biometeorology_parameters' ) == 0 )
1274       READ ( 11, '(A)', END = 20 )  line
1275    ENDDO
1276    BACKSPACE ( 11 )
1277
1278!
1279!-- Read biometeorology_parameters namelist
1280    READ ( 11, biometeorology_parameters, ERR = 10, END = 20 )
1281
1282!
1283!-- Set flag that indicates that the biomet_module is switched on
1284    biometeorology = .TRUE.
1285
1286    GOTO 20
1287
1288!
1289!-- In case of error
1290 10 BACKSPACE( 11 )
1291    READ( 11 , '(A)') line
1292    CALL parin_fail_message( 'biometeorology_parameters', line )
1293
1294!
1295!-- Complete
1296 20 CONTINUE
1297
1298
1299 END SUBROUTINE bio_parin
1300
1301!------------------------------------------------------------------------------!
1302! Description:
1303! ------------
1304!> Soubroutine reads global biometeorology configuration from restart file(s)
1305!------------------------------------------------------------------------------!
1306 SUBROUTINE bio_rrd_global( found )
1307
1308    USE control_parameters,                                                    &
1309        ONLY:  length, restart_string
1310
1311
1312    IMPLICIT NONE
1313
1314    LOGICAL, INTENT(OUT) ::  found      !< variable found? yes = .T., no = .F.
1315
1316    found = .TRUE.
1317
1318
1319    SELECT CASE ( restart_string(1:length) )
1320
1321!
1322!--    read control flags to determine if input grids need to be averaged
1323       CASE ( 'do_average_theta' )
1324          READ ( 13 )  do_average_theta
1325
1326       CASE ( 'do_average_q' )
1327          READ ( 13 )  do_average_q
1328
1329       CASE ( 'do_average_u' )
1330          READ ( 13 )  do_average_u
1331
1332       CASE ( 'do_average_v' )
1333          READ ( 13 )  do_average_v
1334
1335       CASE ( 'do_average_w' )
1336          READ ( 13 )  do_average_w
1337
1338       CASE ( 'do_average_mrt' )
1339          READ ( 13 )  do_average_mrt
1340
1341!
1342!--    read control flags to determine which thermal index needs to trigger averaging
1343       CASE ( 'average_trigger_perct' )
1344          READ ( 13 )  average_trigger_perct
1345
1346       CASE ( 'average_trigger_utci' )
1347          READ ( 13 )  average_trigger_utci
1348
1349       CASE ( 'average_trigger_pet' )
1350          READ ( 13 )  average_trigger_pet
1351
1352
1353       CASE DEFAULT
1354
1355          found = .FALSE.
1356
1357    END SELECT
1358
1359
1360 END SUBROUTINE bio_rrd_global
1361
1362
1363!------------------------------------------------------------------------------!
1364! Description:
1365! ------------
1366!> Soubroutine reads local biometeorology configuration from restart file(s)
1367!------------------------------------------------------------------------------!
1368 SUBROUTINE bio_rrd_local( found )
1369
1370
1371    USE control_parameters,                                                    &
1372        ONLY:  length, restart_string
1373
1374
1375    IMPLICIT NONE
1376
1377
1378    LOGICAL, INTENT(OUT) ::  found      !< variable found? yes = .T., no = .F.
1379
1380    found = .TRUE.
1381
1382
1383    SELECT CASE ( restart_string(1:length) )
1384
1385       CASE ( 'nmrtbl' )
1386          READ ( 13 )  bio_nmrtbl
1387
1388       CASE ( 'mrt_av_grid' )
1389          IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
1390             ALLOCATE( mrt_av_grid(bio_nmrtbl) )
1391          ENDIF
1392          READ ( 13 )  mrt_av_grid
1393
1394
1395       CASE DEFAULT
1396
1397          found = .FALSE.
1398
1399    END SELECT
1400
1401
1402 END SUBROUTINE bio_rrd_local
1403
1404!------------------------------------------------------------------------------!
1405! Description:
1406! ------------
1407!> Write global restart data for the biometeorology module.
1408!------------------------------------------------------------------------------!
1409 SUBROUTINE bio_wrd_global
1410
1411    IMPLICIT NONE
1412
1413    CALL wrd_write_string( 'do_average_theta' )
1414    WRITE ( 14 )  do_average_theta
1415    CALL wrd_write_string( 'do_average_q' )
1416    WRITE ( 14 )  do_average_q
1417    CALL wrd_write_string( 'do_average_u' )
1418    WRITE ( 14 )  do_average_u
1419    CALL wrd_write_string( 'do_average_v' )
1420    WRITE ( 14 )  do_average_v
1421    CALL wrd_write_string( 'do_average_w' )
1422    WRITE ( 14 )  do_average_w
1423    CALL wrd_write_string( 'do_average_mrt' )
1424    WRITE ( 14 )  do_average_mrt
1425    CALL wrd_write_string( 'average_trigger_perct' )
1426    WRITE ( 14 )  average_trigger_perct
1427    CALL wrd_write_string( 'average_trigger_utci' )
1428    WRITE ( 14 )  average_trigger_utci
1429    CALL wrd_write_string( 'average_trigger_pet' )
1430    WRITE ( 14 )  average_trigger_pet
1431
1432 END SUBROUTINE bio_wrd_global
1433
1434
1435!------------------------------------------------------------------------------!
1436! Description:
1437! ------------
1438!> Write local restart data for the biometeorology module.
1439!------------------------------------------------------------------------------!
1440 SUBROUTINE bio_wrd_local
1441
1442    IMPLICIT NONE
1443
1444!
1445!-- First nmrtbl has to be written/read, because it is the dimension of mrt_av_grid
1446    CALL wrd_write_string( 'nmrtbl' )
1447    WRITE ( 14 )  nmrtbl
1448
1449    IF ( ALLOCATED( mrt_av_grid ) )  THEN
1450       CALL wrd_write_string( 'mrt_av_grid' )
1451       WRITE ( 14 )  mrt_av_grid
1452    ENDIF
1453
1454
1455 END SUBROUTINE bio_wrd_local
1456
1457
1458!------------------------------------------------------------------------------!
1459! Description:
1460! ------------
1461!> Calculate biometeorology MRT for all 2D grid
1462!------------------------------------------------------------------------------!
1463 SUBROUTINE bio_calculate_mrt_grid ( av )
1464
1465    IMPLICIT NONE
1466
1467    LOGICAL, INTENT(IN)         ::  av    !< use averaged input?
1468!
1469!-- Internal variables
1470    INTEGER(iwp)                ::  i     !< Running index, x-dir, radiation coordinates
1471    INTEGER(iwp)                ::  j     !< Running index, y-dir, radiation coordinates
1472    INTEGER(iwp)                ::  k     !< Running index, y-dir, radiation coordinates
1473    INTEGER(iwp)                ::  l     !< Running index, radiation coordinates
1474
1475
1476!
1477!-- We need to differentiate if averaged input is desired (av == .TRUE.) or not.
1478    IF ( av )  THEN
1479!
1480!-- Make sure tmrt_av_grid is present and initialize with the fill value
1481       IF ( .NOT. ALLOCATED( tmrt_av_grid ) )  THEN
1482          ALLOCATE( tmrt_av_grid (nys:nyn,nxl:nxr) )
1483       ENDIF
1484       tmrt_av_grid = REAL( bio_fill_value, KIND = wp )
1485
1486!
1487!-- mrt_av_grid should always be allcoated here, but better make sure ist actually is.
1488       IF ( ALLOCATED( mrt_av_grid ) )  THEN
1489!
1490!-- Iterate over the radiation grid (radiation coordinates) and fill the
1491!-- tmrt_av_grid (x, y coordinates) where appropriate:
1492!-- tmrt_av_grid is written for all i / j if level (k) matches output height.
1493          DO  l = 1, nmrtbl
1494             i = mrtbl(ix,l)
1495             j = mrtbl(iy,l)
1496             k = mrtbl(iz,l)
1497             IF ( k - get_topography_top_index_ji( j, i, 's' ) ==              &
1498                   bio_cell_level + 1_iwp)  THEN
1499!
1500!-- Averaging was done before, so we can just copy the result here
1501                tmrt_av_grid(j,i) = mrt_av_grid(l)
1502
1503             ENDIF
1504          ENDDO
1505       ENDIF
1506
1507!
1508!-- In case instantaneous input is desired, mrt values will be re-calculated.
1509    ELSE
1510!
1511!-- Calculate biometeorology MRT from local radiation fluxes calculated by RTM and assign
1512!-- into 2D grid. Depending on selected output quantities, tmrt_grid might not have been
1513!-- allocated in bio_check_data_output yet.
1514       IF ( .NOT. ALLOCATED( tmrt_grid ) )  THEN
1515          ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) )
1516       ENDIF
1517       tmrt_grid = REAL( bio_fill_value, KIND = wp )
1518
1519       DO  l = 1, nmrtbl
1520          i = mrtbl(ix,l)
1521          j = mrtbl(iy,l)
1522          k = mrtbl(iz,l)
1523          IF ( k - get_topography_top_index_ji( j, i, 's' ) ==                 &
1524                bio_cell_level + 1_iwp)  THEN
1525             IF ( mrt_include_sw )  THEN
1526                 tmrt_grid(j,i) = ( ( human_absorb * mrtinsw(l) +              &
1527                                  mrtinlw(l) )  /                              &
1528                                  ( human_emiss * sigma_sb ) )**.25_wp -       &
1529                                  degc_to_k
1530             ELSE
1531                 tmrt_grid(j,i) = ( mrtinlw(l)  /                              &
1532                                  ( human_emiss * sigma_sb ) )**.25_wp -       &
1533                                  degc_to_k
1534             ENDIF
1535          ENDIF
1536       ENDDO
1537    ENDIF
1538
1539END SUBROUTINE bio_calculate_mrt_grid
1540
1541
1542!------------------------------------------------------------------------------!
1543! Description:
1544! ------------
1545!> Calculate static thermal indices for 2D grid point i, j
1546!------------------------------------------------------------------------------!
1547 SUBROUTINE bio_get_thermal_index_input_ij( average_input, i, j, ta, vp, ws,   &
1548    pair, tmrt )
1549
1550    IMPLICIT NONE
1551!
1552!-- Input variables
1553    LOGICAL,      INTENT ( IN ) ::  average_input  !< Determine averaged input conditions?
1554    INTEGER(iwp), INTENT ( IN ) ::  i     !< Running index, x-dir
1555    INTEGER(iwp), INTENT ( IN ) ::  j     !< Running index, y-dir
1556!
1557!-- Output parameters
1558    REAL(wp), INTENT ( OUT )    ::  tmrt  !< Mean radiant temperature        (degree_C)
1559    REAL(wp), INTENT ( OUT )    ::  ta    !< Air temperature                 (degree_C)
1560    REAL(wp), INTENT ( OUT )    ::  vp    !< Vapour pressure                 (hPa)
1561    REAL(wp), INTENT ( OUT )    ::  ws    !< Wind speed    (local level)     (m/s)
1562    REAL(wp), INTENT ( OUT )    ::  pair  !< Air pressure                    (hPa)
1563!
1564!-- Internal variables
1565    INTEGER(iwp)                ::  k     !< Running index, z-dir
1566    INTEGER(iwp)                ::  k_wind  !< Running index, z-dir, wind speed only
1567
1568    REAL(wp)                    ::  vp_sat  !< Saturation vapor pressure     (hPa)
1569
1570!
1571!-- Determine cell level closest to 1.1m above ground
1572!   by making use of truncation due to int cast
1573    k = get_topography_top_index_ji(j, i, 's') + bio_cell_level  !< Vertical cell center closest to 1.1m
1574
1575!
1576!-- Avoid non-representative horizontal u and v of 0.0 m/s too close to ground
1577    k_wind = k
1578    IF ( bio_cell_level < 1_iwp )  THEN
1579       k_wind = k + 1_iwp
1580    ENDIF
1581!
1582!-- Determine local values:
1583    IF ( average_input )  THEN
1584!
1585!--    Calculate ta from Tp assuming dry adiabatic laps rate
1586       ta = pt_av(k,j,i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
1587
1588       vp = bio_fill_value
1589       IF ( humidity  .AND.  ALLOCATED( q_av ) )  THEN
1590          vp = q_av(k,j,i)
1591       ENDIF
1592
1593       ws = ( 0.5_wp * ABS( u_av(k_wind, j, i) + u_av(k_wind, j, i+1) )  +     &
1594          0.5_wp * ABS( v_av(k_wind, j, i) + v_av(k_wind, j+1, i) )  +         &
1595          0.5_wp * ABS( w_av(k_wind, j, i) + w_av(k_wind+1, j, i) ) )
1596    ELSE
1597!
1598!-- Calculate ta from Tp assuming dry adiabatic laps rate
1599       ta = pt(k,j,i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
1600
1601       vp = bio_fill_value
1602       IF ( humidity )  THEN
1603          vp = q(k,j,i)
1604       ENDIF
1605
1606       ws = ( 0.5_wp * ABS( u(k_wind, j, i) + u(k_wind, j, i+1) )  +           &
1607          0.5_wp * ABS( v(k_wind, j, i) + v(k_wind, j+1, i) )  +               &
1608          0.5_wp * ABS( w(k_wind, j, i) + w(k_wind+1, j, i) ) )
1609
1610    ENDIF
1611!
1612!-- Local air pressure
1613    pair = surface_pressure
1614!
1615!-- Calculate water vapour pressure at saturation and convert to hPa
1616!-- The magnus formula is limited to temperatures up to 333.15 K to
1617!   avoid negative values of vp_sat
1618    IF ( vp > -998._wp )  THEN
1619       vp_sat = 0.01_wp * magnus( MIN( ta + degc_to_k, 333.15_wp ) )
1620       vp  = vp * pair / ( vp + 0.622_wp )
1621       IF ( vp > vp_sat )  vp = vp_sat
1622    ENDIF
1623!
1624!-- local mtr value at [i,j]
1625    tmrt = bio_fill_value  !< this can be a valid result (e.g. for inside some ostacle)
1626    IF ( .NOT. average_input )  THEN
1627!
1628!--    Use MRT from RTM precalculated in tmrt_grid
1629       tmrt = tmrt_grid(j,i)
1630    ELSE
1631       tmrt = tmrt_av_grid(j,i)
1632    ENDIF
1633
1634 END SUBROUTINE bio_get_thermal_index_input_ij
1635
1636
1637!------------------------------------------------------------------------------!
1638! Description:
1639! ------------
1640!> Calculate static thermal indices for any point within a 2D grid
1641!> time_integration.f90: 1065ff
1642!------------------------------------------------------------------------------!
1643 SUBROUTINE bio_calculate_thermal_index_maps ( av )
1644
1645    IMPLICIT NONE
1646!
1647!-- Input attributes
1648    LOGICAL, INTENT ( IN ) ::  av  !< Calculate based on averaged input conditions?
1649!
1650!-- Internal variables
1651    INTEGER(iwp) ::  i, j     !< Running index
1652
1653    REAL(wp) ::  clo          !< Clothing index                (no dimension)
1654    REAL(wp) ::  ta           !< Air temperature                  (degree_C)
1655    REAL(wp) ::  vp           !< Vapour pressure                  (hPa)
1656    REAL(wp) ::  ws           !< Wind speed    (local level)      (m/s)
1657    REAL(wp) ::  pair         !< Air pressure                     (hPa)
1658    REAL(wp) ::  perct_ij     !< Perceived temperature            (degree_C)
1659    REAL(wp) ::  utci_ij      !< Universal thermal climate index  (degree_C)
1660    REAL(wp) ::  pet_ij       !< Physiologically equivalent temperature  (degree_C)
1661    REAL(wp) ::  tmrt_ij      !< Mean radiant temperature         (degree_C)
1662
1663!
1664!-- fill out the MRT 2D grid from appropriate source (RTM, RRTMG,...)
1665    CALL bio_calculate_mrt_grid ( av )
1666
1667    DO  i = nxl, nxr
1668       DO  j = nys, nyn
1669!
1670!--       Determine local input conditions
1671          tmrt_ij = bio_fill_value
1672          vp      = bio_fill_value
1673!
1674!--       Determine input only if
1675          CALL bio_get_thermal_index_input_ij ( av, i, j, ta, vp,              &
1676                                                ws, pair, tmrt_ij )
1677!
1678!--       Only proceed if input is available
1679          pet_ij   = bio_fill_value   !< set fail value, e.g. valid for within
1680          perct_ij = bio_fill_value   !< some obstacle
1681          utci_ij  = bio_fill_value
1682          IF ( .NOT. ( tmrt_ij <= -998._wp  .OR.  vp <= -998._wp ) )  THEN
1683!
1684!--          Calculate static thermal indices based on local tmrt
1685             clo = bio_fill_value
1686
1687             IF ( do_calculate_perct  .OR.  do_calculate_perct_av )  THEN
1688!
1689!--          Estimate local perceived temperature
1690                CALL calculate_perct_static( ta, vp, ws, tmrt_ij, pair, clo,   &
1691                   perct_ij )
1692             ENDIF
1693
1694             IF ( do_calculate_utci  .OR.  do_calculate_utci_av )  THEN
1695!
1696!--          Estimate local universal thermal climate index
1697                CALL calculate_utci_static( ta, vp, ws, tmrt_ij,               &
1698                   bio_output_height, utci_ij )
1699             ENDIF
1700
1701             IF ( do_calculate_pet  .OR.  do_calculate_pet_av )  THEN
1702!
1703!--          Estimate local physiologically equivalent temperature
1704                CALL calculate_pet_static( ta, vp, ws, tmrt_ij, pair, pet_ij )
1705             ENDIF
1706          ENDIF
1707
1708
1709          IF ( av )  THEN
1710!
1711!--          Write results for selected averaged indices
1712             IF ( do_calculate_perct_av )  THEN
1713                perct_av(j, i) = perct_ij
1714             ENDIF
1715             IF ( do_calculate_utci_av )  THEN
1716                utci_av(j, i) = utci_ij
1717             ENDIF
1718             IF ( do_calculate_pet_av )  THEN
1719                pet_av(j, i)  = pet_ij
1720             ENDIF
1721          ELSE
1722!
1723!--          Write result for selected indices
1724             IF ( do_calculate_perct )  THEN
1725                perct(j, i) = perct_ij
1726             ENDIF
1727             IF ( do_calculate_utci )  THEN
1728                utci(j, i) = utci_ij
1729             ENDIF
1730             IF ( do_calculate_pet )  THEN
1731                pet(j, i)  = pet_ij
1732             ENDIF
1733          ENDIF
1734
1735       ENDDO
1736    ENDDO
1737
1738 END SUBROUTINE bio_calculate_thermal_index_maps
1739
1740!------------------------------------------------------------------------------!
1741! Description:
1742! ------------
1743!> Calculate dynamic thermal indices (currently only iPT, but expandable)
1744!------------------------------------------------------------------------------!
1745 SUBROUTINE bio_calc_ipt( ta, vp, ws, pair, tmrt, dt, energy_storage,          &
1746    t_clo, clo, actlev, age, weight, height, work, sex, ipt )
1747
1748    IMPLICIT NONE
1749!
1750!-- Input parameters
1751    REAL(wp), INTENT ( IN )  ::  ta   !< Air temperature                  (degree_C)
1752    REAL(wp), INTENT ( IN )  ::  vp   !< Vapour pressure                  (hPa)
1753    REAL(wp), INTENT ( IN )  ::  ws   !< Wind speed    (local level)      (m/s)
1754    REAL(wp), INTENT ( IN )  ::  pair !< Air pressure                     (hPa)
1755    REAL(wp), INTENT ( IN )  ::  tmrt !< Mean radiant temperature         (degree_C)
1756    REAL(wp), INTENT ( IN )  ::  dt   !< Time past since last calculation (s)
1757    REAL(wp), INTENT ( IN )  ::  age  !< Age of agent                     (y)
1758    REAL(wp), INTENT ( IN )  ::  weight  !< Weight of agent               (Kg)
1759    REAL(wp), INTENT ( IN )  ::  height  !< Height of agent               (m)
1760    REAL(wp), INTENT ( IN )  ::  work    !< Mechanical workload of agent
1761                                         !  (without metabolism!)         (W)
1762    INTEGER(iwp), INTENT ( IN ) ::  sex  !< Sex of agent (1 = male, 2 = female)
1763!
1764!-- Both, input and output parameters
1765    Real(wp), INTENT ( INOUT )  ::  energy_storage    !< Energy storage   (W/m²)
1766    Real(wp), INTENT ( INOUT )  ::  t_clo   !< Clothing temperature       (degree_C)
1767    Real(wp), INTENT ( INOUT )  ::  clo     !< Current clothing in sulation
1768    Real(wp), INTENT ( INOUT )  ::  actlev  !< Individuals activity level
1769                                            !  per unit surface area      (W/m²)
1770!
1771!-- Output parameters
1772    REAL(wp), INTENT ( OUT ) ::  ipt    !< Instationary perceived temp.   (degree_C)
1773!
1774!-- return immediatelly if nothing to do!
1775    IF ( .NOT. thermal_comfort )  THEN
1776        RETURN
1777    ENDIF
1778!
1779!-- If clo equals the initial value, this is the initial call
1780    IF ( clo <= -998._wp )  THEN
1781!
1782!--    Initialize instationary perceived temperature with personalized
1783!      PT as an initial guess, set actlev and clo
1784       CALL ipt_init ( age, weight, height, sex, work, actlev, clo,            &
1785          ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo,                   &
1786          ipt )
1787    ELSE
1788!
1789!--    Estimate local instatinoary perceived temperature
1790       CALL ipt_cycle ( ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo,     &
1791          clo, actlev, work, ipt )
1792    ENDIF
1793
1794 END SUBROUTINE bio_calc_ipt
1795
1796
1797
1798!------------------------------------------------------------------------------!
1799! Description:
1800! ------------
1801!> SUBROUTINE for calculating UTCI Temperature (UTCI)
1802!> computed by a 6th order approximation
1803!>
1804!> UTCI regression equation after
1805!> Bröde P, Fiala D, Blazejczyk K, Holmér I, Jendritzky G, Kampmann B, Tinz B,
1806!> Havenith G (2012) Deriving the operational procedure for the Universal Thermal
1807!> Climate Index (UTCI). International Journal of Biometeorology 56 (3):481-494.
1808!> doi:10.1007/s00484-011-0454-1
1809!>
1810!> original source available at:
1811!> www.utci.org
1812!------------------------------------------------------------------------------!
1813 SUBROUTINE calculate_utci_static( ta_in, vp, ws_hag, tmrt, hag, utci_ij )
1814
1815    IMPLICIT NONE
1816!
1817!-- Type of input of the argument list
1818    REAL(WP), INTENT ( IN )  ::  ta_in    !< Local air temperature (degree_C)
1819    REAL(WP), INTENT ( IN )  ::  vp       !< Loacl vapour pressure (hPa)
1820    REAL(WP), INTENT ( IN )  ::  ws_hag   !< Incident wind speed (m/s)
1821    REAL(WP), INTENT ( IN )  ::  tmrt     !< Local mean radiant temperature (degree_C)
1822    REAL(WP), INTENT ( IN )  ::  hag      !< Height of wind speed input (m)
1823!
1824!-- Type of output of the argument list
1825    REAL(wp), INTENT ( OUT ) ::  utci_ij  !< Universal Thermal Climate Index (degree_C)
1826
1827    REAL(WP) ::  ta           !< air temperature modified by offset (degree_C)
1828    REAL(WP) ::  pa           !< air pressure in kPa      (kPa)
1829    REAL(WP) ::  d_tmrt       !< delta-tmrt               (degree_C)
1830    REAL(WP) ::  va           !< wind speed at 10 m above ground level    (m/s)
1831    REAL(WP) ::  offset       !< utci deviation by ta cond. exceeded      (degree_C)
1832    REAL(WP) ::  part_ta      !< Air temperature related part of the regression
1833    REAL(WP) ::  ta2          !< 2 times ta
1834    REAL(WP) ::  ta3          !< 3 times ta
1835    REAL(WP) ::  ta4          !< 4 times ta
1836    REAL(WP) ::  ta5          !< 5 times ta
1837    REAL(WP) ::  ta6          !< 6 times ta
1838    REAL(WP) ::  part_va      !< Vapour pressure related part of the regression
1839    REAL(WP) ::  va2          !< 2 times va
1840    REAL(WP) ::  va3          !< 3 times va
1841    REAL(WP) ::  va4          !< 4 times va
1842    REAL(WP) ::  va5          !< 5 times va
1843    REAL(WP) ::  va6          !< 6 times va
1844    REAL(WP) ::  part_d_tmrt  !< Mean radiant temp. related part of the reg.
1845    REAL(WP) ::  d_tmrt2      !< 2 times d_tmrt
1846    REAL(WP) ::  d_tmrt3      !< 3 times d_tmrt
1847    REAL(WP) ::  d_tmrt4      !< 4 times d_tmrt
1848    REAL(WP) ::  d_tmrt5      !< 5 times d_tmrt
1849    REAL(WP) ::  d_tmrt6      !< 6 times d_tmrt
1850    REAL(WP) ::  part_pa      !< Air pressure related part of the regression
1851    REAL(WP) ::  pa2          !< 2 times pa
1852    REAL(WP) ::  pa3          !< 3 times pa
1853    REAL(WP) ::  pa4          !< 4 times pa
1854    REAL(WP) ::  pa5          !< 5 times pa
1855    REAL(WP) ::  pa6          !< 6 times pa
1856    REAL(WP) ::  part_pa2     !< Air pressure^2 related part of the regression
1857    REAL(WP) ::  part_pa3     !< Air pressure^3 related part of the regression
1858    REAL(WP) ::  part_pa46    !< Air pressure^4-6 related part of the regression
1859
1860!
1861!-- Initialize
1862    offset = 0._wp
1863    ta = ta_in
1864    d_tmrt = tmrt - ta_in
1865!
1866!-- Use vapour pressure in kpa
1867    pa = vp / 10.0_wp
1868!
1869!-- Wind altitude correction from hag to 10m after Broede et al. (2012), eq.3
1870!-- z(0) is set to 0.01 according to UTCI profile definition
1871    va = ws_hag *  log ( 10.0_wp / 0.01_wp ) / log ( hag / 0.01_wp )
1872!
1873!-- Check if input values in range after Broede et al. (2012)
1874    IF ( ( d_tmrt > 70._wp )  .OR.  ( d_tmrt < -30._wp )  .OR.                 &
1875       ( vp >= 50._wp ) )  THEN
1876       utci_ij = bio_fill_value
1877       RETURN
1878    ENDIF
1879!
1880!-- Apply eq. 2 in Broede et al. (2012) for ta out of bounds
1881    IF ( ta > 50._wp )  THEN
1882       offset = ta - 50._wp
1883       ta = 50._wp
1884    ENDIF
1885    IF ( ta < -50._wp )  THEN
1886       offset = ta + 50._wp
1887       ta = -50._wp
1888    ENDIF
1889!
1890!-- For routine application. For wind speeds and relative
1891!-- humidity values below 0.5 m/s or 5%, respectively, the
1892!-- user is advised to use the lower bounds for the calculations.
1893    IF ( va < 0.5_wp )  va = 0.5_wp
1894    IF ( va > 17._wp )  va = 17._wp
1895
1896!
1897!-- Pre-calculate multiples of input parameters to save time later
1898
1899    ta2 = ta  * ta
1900    ta3 = ta2 * ta
1901    ta4 = ta3 * ta
1902    ta5 = ta4 * ta
1903    ta6 = ta5 * ta
1904
1905    va2 = va  * va
1906    va3 = va2 * va
1907    va4 = va3 * va
1908    va5 = va4 * va
1909    va6 = va5 * va
1910
1911    d_tmrt2 = d_tmrt  * d_tmrt
1912    d_tmrt3 = d_tmrt2 * d_tmrt
1913    d_tmrt4 = d_tmrt3 * d_tmrt
1914    d_tmrt5 = d_tmrt4 * d_tmrt
1915    d_tmrt6 = d_tmrt5 * d_tmrt
1916
1917    pa2 = pa  * pa
1918    pa3 = pa2 * pa
1919    pa4 = pa3 * pa
1920    pa5 = pa4 * pa
1921    pa6 = pa5 * pa
1922
1923!
1924!-- Pre-calculate parts of the regression equation
1925    part_ta = (  6.07562052e-01_wp )       +                                   &
1926              ( -2.27712343e-02_wp ) * ta  +                                   &
1927              (  8.06470249e-04_wp ) * ta2 +                                   &
1928              ( -1.54271372e-04_wp ) * ta3 +                                   &
1929              ( -3.24651735e-06_wp ) * ta4 +                                   &
1930              (  7.32602852e-08_wp ) * ta5 +                                   &
1931              (  1.35959073e-09_wp ) * ta6
1932
1933    part_va = ( -2.25836520e+00_wp ) * va +                                    &
1934        (  8.80326035e-02_wp ) * ta  * va +                                    &
1935        (  2.16844454e-03_wp ) * ta2 * va +                                    &
1936        ( -1.53347087e-05_wp ) * ta3 * va +                                    &
1937        ( -5.72983704e-07_wp ) * ta4 * va +                                    &
1938        ( -2.55090145e-09_wp ) * ta5 * va +                                    &
1939        ( -7.51269505e-01_wp ) *       va2 +                                   &
1940        ( -4.08350271e-03_wp ) * ta  * va2 +                                   &
1941        ( -5.21670675e-05_wp ) * ta2 * va2 +                                   &
1942        (  1.94544667e-06_wp ) * ta3 * va2 +                                   &
1943        (  1.14099531e-08_wp ) * ta4 * va2 +                                   &
1944        (  1.58137256e-01_wp ) *       va3 +                                   &
1945        ( -6.57263143e-05_wp ) * ta  * va3 +                                   &
1946        (  2.22697524e-07_wp ) * ta2 * va3 +                                   &
1947        ( -4.16117031e-08_wp ) * ta3 * va3 +                                   &
1948        ( -1.27762753e-02_wp ) *       va4 +                                   &
1949        (  9.66891875e-06_wp ) * ta  * va4 +                                   &
1950        (  2.52785852e-09_wp ) * ta2 * va4 +                                   &
1951        (  4.56306672e-04_wp ) *       va5 +                                   &
1952        ( -1.74202546e-07_wp ) * ta  * va5 +                                   &
1953        ( -5.91491269e-06_wp ) * va6
1954
1955    part_d_tmrt = (  3.98374029e-01_wp ) *       d_tmrt +                      &
1956            (  1.83945314e-04_wp ) * ta  *       d_tmrt +                      &
1957            ( -1.73754510e-04_wp ) * ta2 *       d_tmrt +                      &
1958            ( -7.60781159e-07_wp ) * ta3 *       d_tmrt +                      &
1959            (  3.77830287e-08_wp ) * ta4 *       d_tmrt +                      &
1960            (  5.43079673e-10_wp ) * ta5 *       d_tmrt +                      &
1961            ( -2.00518269e-02_wp ) *       va  * d_tmrt +                      &
1962            (  8.92859837e-04_wp ) * ta  * va  * d_tmrt +                      &
1963            (  3.45433048e-06_wp ) * ta2 * va  * d_tmrt +                      &
1964            ( -3.77925774e-07_wp ) * ta3 * va  * d_tmrt +                      &
1965            ( -1.69699377e-09_wp ) * ta4 * va  * d_tmrt +                      &
1966            (  1.69992415e-04_wp ) *       va2 * d_tmrt +                      &
1967            ( -4.99204314e-05_wp ) * ta  * va2 * d_tmrt +                      &
1968            (  2.47417178e-07_wp ) * ta2 * va2 * d_tmrt +                      &
1969            (  1.07596466e-08_wp ) * ta3 * va2 * d_tmrt +                      &
1970            (  8.49242932e-05_wp ) *       va3 * d_tmrt +                      &
1971            (  1.35191328e-06_wp ) * ta  * va3 * d_tmrt +                      &
1972            ( -6.21531254e-09_wp ) * ta2 * va3 * d_tmrt +                      &
1973            ( -4.99410301e-06_wp ) * va4 *       d_tmrt +                      &
1974            ( -1.89489258e-08_wp ) * ta  * va4 * d_tmrt +                      &
1975            (  8.15300114e-08_wp ) *       va5 * d_tmrt +                      &
1976            (  7.55043090e-04_wp ) *             d_tmrt2 +                     &
1977            ( -5.65095215e-05_wp ) * ta  *       d_tmrt2 +                     &
1978            ( -4.52166564e-07_wp ) * ta2 *       d_tmrt2 +                     &
1979            (  2.46688878e-08_wp ) * ta3 *       d_tmrt2 +                     &
1980            (  2.42674348e-10_wp ) * ta4 *       d_tmrt2 +                     &
1981            (  1.54547250e-04_wp ) *       va  * d_tmrt2 +                     &
1982            (  5.24110970e-06_wp ) * ta  * va  * d_tmrt2 +                     &
1983            ( -8.75874982e-08_wp ) * ta2 * va  * d_tmrt2 +                     &
1984            ( -1.50743064e-09_wp ) * ta3 * va  * d_tmrt2 +                     &
1985            ( -1.56236307e-05_wp ) *       va2 * d_tmrt2 +                     &
1986            ( -1.33895614e-07_wp ) * ta  * va2 * d_tmrt2 +                     &
1987            (  2.49709824e-09_wp ) * ta2 * va2 * d_tmrt2 +                     &
1988            (  6.51711721e-07_wp ) *       va3 * d_tmrt2 +                     &
1989            (  1.94960053e-09_wp ) * ta  * va3 * d_tmrt2 +                     &
1990            ( -1.00361113e-08_wp ) *       va4 * d_tmrt2 +                     &
1991            ( -1.21206673e-05_wp ) *             d_tmrt3 +                     &
1992            ( -2.18203660e-07_wp ) * ta  *       d_tmrt3 +                     &
1993            (  7.51269482e-09_wp ) * ta2 *       d_tmrt3 +                     &
1994            (  9.79063848e-11_wp ) * ta3 *       d_tmrt3 +                     &
1995            (  1.25006734e-06_wp ) *       va  * d_tmrt3 +                     &
1996            ( -1.81584736e-09_wp ) * ta  * va  * d_tmrt3 +                     &
1997            ( -3.52197671e-10_wp ) * ta2 * va  * d_tmrt3 +                     &
1998            ( -3.36514630e-08_wp ) *       va2 * d_tmrt3 +                     &
1999            (  1.35908359e-10_wp ) * ta  * va2 * d_tmrt3 +                     &
2000            (  4.17032620e-10_wp ) *       va3 * d_tmrt3 +                     &
2001            ( -1.30369025e-09_wp ) *             d_tmrt4 +                     &
2002            (  4.13908461e-10_wp ) * ta  *       d_tmrt4 +                     &
2003            (  9.22652254e-12_wp ) * ta2 *       d_tmrt4 +                     &
2004            ( -5.08220384e-09_wp ) *       va  * d_tmrt4 +                     &
2005            ( -2.24730961e-11_wp ) * ta  * va  * d_tmrt4 +                     &
2006            (  1.17139133e-10_wp ) *       va2 * d_tmrt4 +                     &
2007            (  6.62154879e-10_wp ) *             d_tmrt5 +                     &
2008            (  4.03863260e-13_wp ) * ta  *       d_tmrt5 +                     &
2009            (  1.95087203e-12_wp ) *       va  * d_tmrt5 +                     &
2010            ( -4.73602469e-12_wp ) *             d_tmrt6
2011
2012    part_pa = (  5.12733497e+00_wp ) *                pa +                     &
2013       ( -3.12788561e-01_wp ) * ta  *                 pa +                     &
2014       ( -1.96701861e-02_wp ) * ta2 *                 pa +                     &
2015       (  9.99690870e-04_wp ) * ta3 *                 pa +                     &
2016       (  9.51738512e-06_wp ) * ta4 *                 pa +                     &
2017       ( -4.66426341e-07_wp ) * ta5 *                 pa +                     &
2018       (  5.48050612e-01_wp ) *       va  *           pa +                     &
2019       ( -3.30552823e-03_wp ) * ta  * va  *           pa +                     &
2020       ( -1.64119440e-03_wp ) * ta2 * va  *           pa +                     &
2021       ( -5.16670694e-06_wp ) * ta3 * va  *           pa +                     &
2022       (  9.52692432e-07_wp ) * ta4 * va  *           pa +                     &
2023       ( -4.29223622e-02_wp ) *       va2 *           pa +                     &
2024       (  5.00845667e-03_wp ) * ta  * va2 *           pa +                     &
2025       (  1.00601257e-06_wp ) * ta2 * va2 *           pa +                     &
2026       ( -1.81748644e-06_wp ) * ta3 * va2 *           pa +                     &
2027       ( -1.25813502e-03_wp ) *       va3 *           pa +                     &
2028       ( -1.79330391e-04_wp ) * ta  * va3 *           pa +                     &
2029       (  2.34994441e-06_wp ) * ta2 * va3 *           pa +                     &
2030       (  1.29735808e-04_wp ) *       va4 *           pa +                     &
2031       (  1.29064870e-06_wp ) * ta  * va4 *           pa +                     &
2032       ( -2.28558686e-06_wp ) *       va5 *           pa +                     &
2033       ( -3.69476348e-02_wp ) *             d_tmrt  * pa +                     &
2034       (  1.62325322e-03_wp ) * ta  *       d_tmrt  * pa +                     &
2035       ( -3.14279680e-05_wp ) * ta2 *       d_tmrt  * pa +                     &
2036       (  2.59835559e-06_wp ) * ta3 *       d_tmrt  * pa +                     &
2037       ( -4.77136523e-08_wp ) * ta4 *       d_tmrt  * pa +                     &
2038       (  8.64203390e-03_wp ) *       va  * d_tmrt  * pa +                     &
2039       ( -6.87405181e-04_wp ) * ta  * va  * d_tmrt  * pa +                     &
2040       ( -9.13863872e-06_wp ) * ta2 * va  * d_tmrt  * pa +                     &
2041       (  5.15916806e-07_wp ) * ta3 * va  * d_tmrt  * pa +                     &
2042       ( -3.59217476e-05_wp ) *       va2 * d_tmrt  * pa +                     &
2043       (  3.28696511e-05_wp ) * ta  * va2 * d_tmrt  * pa +                     &
2044       ( -7.10542454e-07_wp ) * ta2 * va2 * d_tmrt  * pa +                     &
2045       ( -1.24382300e-05_wp ) *       va3 * d_tmrt  * pa +                     &
2046       ( -7.38584400e-09_wp ) * ta  * va3 * d_tmrt  * pa +                     &
2047       (  2.20609296e-07_wp ) *       va4 * d_tmrt  * pa +                     &
2048       ( -7.32469180e-04_wp ) *             d_tmrt2 * pa +                     &
2049       ( -1.87381964e-05_wp ) * ta  *       d_tmrt2 * pa +                     &
2050       (  4.80925239e-06_wp ) * ta2 *       d_tmrt2 * pa +                     &
2051       ( -8.75492040e-08_wp ) * ta3 *       d_tmrt2 * pa +                     &
2052       (  2.77862930e-05_wp ) *       va  * d_tmrt2 * pa +                     &
2053       ( -5.06004592e-06_wp ) * ta  * va  * d_tmrt2 * pa +                     &
2054       (  1.14325367e-07_wp ) * ta2 * va  * d_tmrt2 * pa +                     &
2055       (  2.53016723e-06_wp ) *       va2 * d_tmrt2 * pa +                     &
2056       ( -1.72857035e-08_wp ) * ta  * va2 * d_tmrt2 * pa +                     &
2057       ( -3.95079398e-08_wp ) *       va3 * d_tmrt2 * pa +                     &
2058       ( -3.59413173e-07_wp ) *             d_tmrt3 * pa +                     &
2059       (  7.04388046e-07_wp ) * ta  *       d_tmrt3 * pa +                     &
2060       ( -1.89309167e-08_wp ) * ta2 *       d_tmrt3 * pa +                     &
2061       ( -4.79768731e-07_wp ) *       va  * d_tmrt3 * pa +                     &
2062       (  7.96079978e-09_wp ) * ta  * va  * d_tmrt3 * pa +                     &
2063       (  1.62897058e-09_wp ) *       va2 * d_tmrt3 * pa +                     &
2064       (  3.94367674e-08_wp ) *             d_tmrt4 * pa +                     &
2065       ( -1.18566247e-09_wp ) * ta *        d_tmrt4 * pa +                     &
2066       (  3.34678041e-10_wp ) *       va  * d_tmrt4 * pa +                     &
2067       ( -1.15606447e-10_wp ) *             d_tmrt5 * pa
2068
2069    part_pa2 = ( -2.80626406e+00_wp ) *               pa2 +                    &
2070       (  5.48712484e-01_wp ) * ta  *                 pa2 +                    &
2071       ( -3.99428410e-03_wp ) * ta2 *                 pa2 +                    &
2072       ( -9.54009191e-04_wp ) * ta3 *                 pa2 +                    &
2073       (  1.93090978e-05_wp ) * ta4 *                 pa2 +                    &
2074       ( -3.08806365e-01_wp ) *       va *            pa2 +                    &
2075       (  1.16952364e-02_wp ) * ta  * va *            pa2 +                    &
2076       (  4.95271903e-04_wp ) * ta2 * va *            pa2 +                    &
2077       ( -1.90710882e-05_wp ) * ta3 * va *            pa2 +                    &
2078       (  2.10787756e-03_wp ) *       va2 *           pa2 +                    &
2079       ( -6.98445738e-04_wp ) * ta  * va2 *           pa2 +                    &
2080       (  2.30109073e-05_wp ) * ta2 * va2 *           pa2 +                    &
2081       (  4.17856590e-04_wp ) *       va3 *           pa2 +                    &
2082       ( -1.27043871e-05_wp ) * ta  * va3 *           pa2 +                    &
2083       ( -3.04620472e-06_wp ) *       va4 *           pa2 +                    &
2084       (  5.14507424e-02_wp ) *             d_tmrt  * pa2 +                    &
2085       ( -4.32510997e-03_wp ) * ta  *       d_tmrt  * pa2 +                    &
2086       (  8.99281156e-05_wp ) * ta2 *       d_tmrt  * pa2 +                    &
2087       ( -7.14663943e-07_wp ) * ta3 *       d_tmrt  * pa2 +                    &
2088       ( -2.66016305e-04_wp ) *       va  * d_tmrt  * pa2 +                    &
2089       (  2.63789586e-04_wp ) * ta  * va  * d_tmrt  * pa2 +                    &
2090       ( -7.01199003e-06_wp ) * ta2 * va  * d_tmrt  * pa2 +                    &
2091       ( -1.06823306e-04_wp ) *       va2 * d_tmrt  * pa2 +                    &
2092       (  3.61341136e-06_wp ) * ta  * va2 * d_tmrt  * pa2 +                    &
2093       (  2.29748967e-07_wp ) *       va3 * d_tmrt  * pa2 +                    &
2094       (  3.04788893e-04_wp ) *             d_tmrt2 * pa2 +                    &
2095       ( -6.42070836e-05_wp ) * ta  *       d_tmrt2 * pa2 +                    &
2096       (  1.16257971e-06_wp ) * ta2 *       d_tmrt2 * pa2 +                    &
2097       (  7.68023384e-06_wp ) *       va  * d_tmrt2 * pa2 +                    &
2098       ( -5.47446896e-07_wp ) * ta  * va  * d_tmrt2 * pa2 +                    &
2099       ( -3.59937910e-08_wp ) *       va2 * d_tmrt2 * pa2 +                    &
2100       ( -4.36497725e-06_wp ) *             d_tmrt3 * pa2 +                    &
2101       (  1.68737969e-07_wp ) * ta  *       d_tmrt3 * pa2 +                    &
2102       (  2.67489271e-08_wp ) *       va  * d_tmrt3 * pa2 +                    &
2103       (  3.23926897e-09_wp ) *             d_tmrt4 * pa2
2104
2105    part_pa3 = ( -3.53874123e-02_wp ) *               pa3 +                    &
2106       ( -2.21201190e-01_wp ) * ta  *                 pa3 +                    &
2107       (  1.55126038e-02_wp ) * ta2 *                 pa3 +                    &
2108       ( -2.63917279e-04_wp ) * ta3 *                 pa3 +                    &
2109       (  4.53433455e-02_wp ) *       va  *           pa3 +                    &
2110       ( -4.32943862e-03_wp ) * ta  * va  *           pa3 +                    &
2111       (  1.45389826e-04_wp ) * ta2 * va  *           pa3 +                    &
2112       (  2.17508610e-04_wp ) *       va2 *           pa3 +                    &
2113       ( -6.66724702e-05_wp ) * ta  * va2 *           pa3 +                    &
2114       (  3.33217140e-05_wp ) *       va3 *           pa3 +                    &
2115       ( -2.26921615e-03_wp ) *             d_tmrt  * pa3 +                    &
2116       (  3.80261982e-04_wp ) * ta  *       d_tmrt  * pa3 +                    &
2117       ( -5.45314314e-09_wp ) * ta2 *       d_tmrt  * pa3 +                    &
2118       ( -7.96355448e-04_wp ) *       va  * d_tmrt  * pa3 +                    &
2119       (  2.53458034e-05_wp ) * ta  * va  * d_tmrt  * pa3 +                    &
2120       ( -6.31223658e-06_wp ) *       va2 * d_tmrt  * pa3 +                    &
2121       (  3.02122035e-04_wp ) *             d_tmrt2 * pa3 +                    &
2122       ( -4.77403547e-06_wp ) * ta  *       d_tmrt2 * pa3 +                    &
2123       (  1.73825715e-06_wp ) *       va  * d_tmrt2 * pa3 +                    &
2124       ( -4.09087898e-07_wp ) *             d_tmrt3 * pa3
2125
2126    part_pa46 = (  6.14155345e-01_wp ) *              pa4 +                    &
2127       ( -6.16755931e-02_wp ) * ta  *                 pa4 +                    &
2128       (  1.33374846e-03_wp ) * ta2 *                 pa4 +                    &
2129       (  3.55375387e-03_wp ) *       va  *           pa4 +                    &
2130       ( -5.13027851e-04_wp ) * ta  * va  *           pa4 +                    &
2131       (  1.02449757e-04_wp ) *       va2 *           pa4 +                    &
2132       ( -1.48526421e-03_wp ) *             d_tmrt  * pa4 +                    &
2133       ( -4.11469183e-05_wp ) * ta  *       d_tmrt  * pa4 +                    &
2134       ( -6.80434415e-06_wp ) *       va  * d_tmrt  * pa4 +                    &
2135       ( -9.77675906e-06_wp ) *             d_tmrt2 * pa4 +                    &
2136       (  8.82773108e-02_wp ) *                       pa5 +                    &
2137       ( -3.01859306e-03_wp ) * ta  *                 pa5 +                    &
2138       (  1.04452989e-03_wp ) *       va  *           pa5 +                    &
2139       (  2.47090539e-04_wp ) *             d_tmrt  * pa5 +                    &
2140       (  1.48348065e-03_wp ) *                       pa6
2141!
2142!-- Calculate 6th order polynomial as approximation
2143    utci_ij = ta + part_ta + part_va + part_d_tmrt + part_pa + part_pa2 +      &
2144        part_pa3 + part_pa46
2145!
2146!-- Consider offset in result
2147    utci_ij = utci_ij + offset
2148
2149 END SUBROUTINE calculate_utci_static
2150
2151
2152
2153
2154!------------------------------------------------------------------------------!
2155! Description:
2156! ------------
2157!> calculate_perct_static: Estimation of perceived temperature (PT, degree_C)
2158!> Value of perct is the Perceived Temperature, degree centigrade
2159!------------------------------------------------------------------------------!
2160 SUBROUTINE calculate_perct_static( ta, vp, ws, tmrt, pair, clo, perct_ij )
2161
2162    IMPLICIT NONE
2163!
2164!-- Type of input of the argument list
2165    REAL(wp), INTENT ( IN )  :: ta   !< Local air temperature (degC)
2166    REAL(wp), INTENT ( IN )  :: vp   !< Local vapour pressure (hPa)
2167    REAL(wp), INTENT ( IN )  :: tmrt !< Local mean radiant temperature (degC)
2168    REAL(wp), INTENT ( IN )  :: ws   !< Local wind velocitry (m/s)
2169    REAL(wp), INTENT ( IN )  :: pair !< Local barometric air pressure (hPa)
2170!
2171!-- Type of output of the argument list
2172    REAL(wp), INTENT ( OUT ) :: perct_ij  !< Perceived temperature (degC)
2173    REAL(wp), INTENT ( OUT ) :: clo       !< Clothing index (dimensionless)
2174!
2175!-- Parameters for standard "Klima-Michel"
2176    REAL(wp), PARAMETER :: eta = 0._wp  !< Mechanical work efficiency for walking on flat ground
2177                                        !< (compare to Fanger (1972) pp 24f)
2178    REAL(wp), PARAMETER :: actlev = 134.6862_wp  !< Workload by activity per standardized surface (A_Du)
2179!
2180!-- Type of program variables
2181    REAL(wp), PARAMETER :: eps = 0.0005  !< Accuracy in clothing insulation (clo) for evaluation the root of Fanger's PMV (pmva=0)
2182    REAL(wp) ::  sclo           !< summer clothing insulation
2183    REAL(wp) ::  wclo           !< winter clothing insulation
2184    REAL(wp) ::  d_pmv          !< PMV deviation (dimensionless --> PMV)
2185    REAL(wp) ::  svp_ta         !< saturation vapor pressure    (hPa)
2186    REAL(wp) ::  sult_lim       !< threshold for sultrieness    (hPa)
2187    REAL(wp) ::  dgtcm          !< Mean deviation dependent on perct
2188    REAL(wp) ::  dgtcstd        !< Mean deviation plus its standard deviation
2189    REAL(wp) ::  clon           !< clo for neutral conditions   (clo)
2190    REAL(wp) ::  ireq_minimal   !< Minimal required clothing insulation (clo)
2191!     REAL(wp) ::  clo_fanger     !< clo for fanger subroutine, unused
2192    REAL(wp) ::  pmv_w          !< Fangers predicted mean vote for winter clothing
2193    REAL(wp) ::  pmv_s          !< Fangers predicted mean vote for summer clothing
2194    REAL(wp) ::  pmva           !< adjusted predicted mean vote
2195    REAL(wp) ::  ptc            !< perceived temp. for cold conditions (degree_C)
2196    REAL(wp) ::  d_std          !< factor to threshold for sultriness
2197    REAL(wp) ::  pmvs           !< pred. mean vote considering sultrieness
2198    REAL(wp) ::  top            !< Gagge's operative temperatures (degree_C)
2199
2200    INTEGER(iwp) :: ncount      !< running index
2201    INTEGER(iwp) :: nerr_cold   !< error number (cold conditions)
2202    INTEGER(iwp) :: nerr        !< error number
2203
2204    LOGICAL :: sultrieness
2205!
2206!-- Initialise
2207    perct_ij = bio_fill_value
2208
2209    nerr     = 0_iwp
2210    ncount   = 0_iwp
2211    sultrieness  = .FALSE.
2212!
2213!-- Tresholds: clothing insulation (account for model inaccuracies)
2214!
2215!-- summer clothing
2216    sclo     = 0.44453_wp
2217!
2218!-- winter clothing
2219    wclo     = 1.76267_wp
2220!
2221!-- decision: firstly calculate for winter or summer clothing
2222    IF ( ta <= 10._wp )  THEN
2223!
2224!--    First guess: winter clothing insulation: cold stress
2225       clo = wclo
2226       CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva, top )
2227       pmv_w = pmva
2228
2229       IF ( pmva > 0._wp )  THEN
2230!
2231!--       Case summer clothing insulation: heat load ?
2232          clo = sclo
2233          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
2234             top )
2235          pmv_s = pmva
2236          IF ( pmva <= 0._wp )  THEN
2237!
2238!--          Case: comfort achievable by varying clothing insulation
2239!--          Between winter and summer set values
2240             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
2241                pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
2242             IF ( ncount < 0_iwp )  THEN
2243                nerr = -1_iwp
2244                RETURN
2245             ENDIF
2246          ELSE IF ( pmva > 0.06_wp )  THEN
2247             clo = 0.5_wp
2248             CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,           &
2249                           pmva, top )
2250          ENDIF
2251       ELSE IF ( pmva < -0.11_wp )  THEN
2252          clo = 1.75_wp
2253          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
2254             top )
2255       ENDIF
2256    ELSE
2257!
2258!--    First guess: summer clothing insulation: heat load
2259       clo = sclo
2260       CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva, top )
2261       pmv_s = pmva
2262
2263       IF ( pmva < 0._wp )  THEN
2264!
2265!--       Case winter clothing insulation: cold stress ?
2266          clo = wclo
2267          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
2268             top )
2269          pmv_w = pmva
2270
2271          IF ( pmva >= 0._wp )  THEN
2272!
2273!--          Case: comfort achievable by varying clothing insulation
2274!--          between winter and summer set values
2275             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
2276                               pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
2277             IF ( ncount < 0_iwp )  THEN
2278                nerr = -1_iwp
2279                RETURN
2280             ENDIF
2281          ELSE IF ( pmva < -0.11_wp )  THEN
2282             clo = 1.75_wp
2283             CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,           &
2284                           pmva, top )
2285          ENDIF
2286       ELSE IF ( pmva > 0.06_wp )  THEN
2287          clo = 0.5_wp
2288          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
2289             top )
2290       ENDIF
2291
2292    ENDIF
2293!
2294!-- Determine perceived temperature by regression equation + adjustments
2295    pmvs = pmva
2296    CALL perct_regression ( pmva, clo, perct_ij )
2297    ptc = perct_ij
2298    IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp )  THEN
2299!
2300!--    Adjust for cold conditions according to Gagge 1986
2301       CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
2302       IF ( nerr_cold > 0_iwp )  nerr = -5_iwp
2303       pmvs = pmva - d_pmv
2304       IF ( pmvs > -0.11_wp )  THEN
2305          d_pmv  = 0._wp
2306          pmvs   = -0.11_wp
2307       ENDIF
2308       CALL perct_regression ( pmvs, clo, perct_ij )
2309    ENDIF
2310!     clo_fanger = clo
2311    clon = clo
2312    IF ( clo > 0.5_wp  .AND.  perct_ij <= 8.73_wp )  THEN
2313!
2314!--    Required clothing insulation (ireq) is exclusively defined for
2315!--    operative temperatures (top) less 10 (C) for a
2316!--    reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
2317       clon = ireq_neutral ( perct_ij, ireq_minimal, nerr )
2318       clo = clon
2319    ENDIF
2320    CALL calc_sultr ( ptc, dgtcm, dgtcstd, sult_lim )
2321    sultrieness    = .FALSE.
2322    d_std = -99._wp
2323    IF ( pmva > 0.06_wp  .AND.  clo <= 0.5_wp )  THEN
2324!
2325!--    Adjust for warm/humid conditions according to Gagge 1986
2326       CALL saturation_vapor_pressure ( ta, svp_ta )
2327       d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
2328       pmvs   = pmva + d_pmv
2329       CALL perct_regression ( pmvs, clo, perct_ij )
2330       IF ( sult_lim < 99._wp )  THEN
2331          IF ( (perct_ij - ptc) > sult_lim )  sultrieness = .TRUE.
2332!
2333!--       Set factor to threshold for sultriness
2334          IF ( dgtcstd /= 0_iwp )  THEN
2335             d_std = ( ( perct_ij - ptc ) - dgtcm ) / dgtcstd
2336          ENDIF
2337       ENDIF
2338    ENDIF
2339
2340 END SUBROUTINE calculate_perct_static
2341
2342!------------------------------------------------------------------------------!
2343! Description:
2344! ------------
2345!> The SUBROUTINE calculates the (saturation) water vapour pressure
2346!> (hPa = hecto Pascal) for a given temperature ta (degC).
2347!> 'ta' can be the air temperature or the dew point temperature. The first will
2348!> result in the current vapor pressure (hPa), the latter will calulate the
2349!> saturation vapor pressure (hPa).
2350!------------------------------------------------------------------------------!
2351 SUBROUTINE saturation_vapor_pressure( ta, svp_ta )
2352
2353    IMPLICIT NONE
2354
2355    REAL(wp), INTENT ( IN )  ::  ta     !< ambient air temperature (degC)
2356    REAL(wp), INTENT ( OUT ) ::  svp_ta !< water vapour pressure (hPa)
2357
2358    REAL(wp)      ::  b
2359    REAL(wp)      ::  c
2360
2361
2362    IF ( ta < 0._wp )  THEN
2363!
2364!--    ta  < 0 (degC): water vapour pressure over ice
2365       b = 17.84362_wp
2366       c = 245.425_wp
2367    ELSE
2368!
2369!--    ta >= 0 (degC): water vapour pressure over water
2370       b = 17.08085_wp
2371       c = 234.175_wp
2372    ENDIF
2373!
2374!-- Saturation water vapour pressure
2375    svp_ta = 6.1078_wp * EXP ( b * ta / ( c + ta ) )
2376
2377 END SUBROUTINE saturation_vapor_pressure
2378
2379!------------------------------------------------------------------------------!
2380! Description:
2381! ------------
2382!> Find the clothing insulation value clo_res (clo) to make Fanger's Predicted
2383!> Mean Vote (PMV) equal comfort (pmva=0) for actual meteorological conditions
2384!> (ta,tmrt, vp, ws, pair) and values of individual's activity level
2385!------------------------------------------------------------------------------!
2386 SUBROUTINE iso_ridder( ta, tmrt, vp, ws, pair, actlev, eta, sclo,             &
2387                       pmv_s, wclo, pmv_w, eps, pmva, top, nerr,               &
2388                       clo_res )
2389
2390    IMPLICIT NONE
2391!
2392!-- Input variables of argument list:
2393    REAL(wp), INTENT ( IN )  :: ta       !< Ambient temperature (degC)
2394    REAL(wp), INTENT ( IN )  :: tmrt     !< Mean radiant temperature (degC)
2395    REAL(wp), INTENT ( IN )  :: vp       !< Water vapour pressure (hPa)
2396    REAL(wp), INTENT ( IN )  :: ws       !< Wind speed (m/s) 1 m above ground
2397    REAL(wp), INTENT ( IN )  :: pair     !< Barometric air pressure (hPa)
2398    REAL(wp), INTENT ( IN )  :: actlev   !< Individuals activity level per unit surface area (W/m2)
2399    REAL(wp), INTENT ( IN )  :: eta      !< Individuals work efficiency (dimensionless)
2400    REAL(wp), INTENT ( IN )  :: sclo     !< Lower threshold of bracketing clothing insulation (clo)
2401    REAL(wp), INTENT ( IN )  :: wclo     !< Upper threshold of bracketing clothing insulation (clo)
2402    REAL(wp), INTENT ( IN )  :: eps      !< (0.05) accuracy in clothing insulation (clo) for
2403!                                          evaluation the root of Fanger's PMV (pmva=0)
2404    REAL(wp), INTENT ( IN )  :: pmv_w    !< Fanger's PMV corresponding to wclo
2405    REAL(wp), INTENT ( IN )  :: pmv_s    !< Fanger's PMV corresponding to sclo
2406!
2407!-- Output variables of argument list:
2408    REAL(wp), INTENT ( OUT ) :: pmva     !< 0 (set to zero, because clo is evaluated for comfort)
2409    REAL(wp), INTENT ( OUT ) :: top      !< Operative temperature (degC) at found root of Fanger's PMV
2410    REAL(wp), INTENT ( OUT ) :: clo_res  !< Resulting clothing insulation value (clo)
2411    INTEGER(iwp), INTENT ( OUT ) :: nerr !< Error status / quality flag
2412                                         !< nerr >= 0, o.k., and nerr is the number of iterations for convergence
2413                                         !< nerr = -1: error = malfunction of Ridder's convergence method
2414                                         !< nerr = -2: error = maximum iterations (max_iteration) exceeded
2415                                         !< nerr = -3: error = root not bracketed between sclo and wclo
2416!
2417!-- Type of program variables
2418    INTEGER(iwp), PARAMETER  ::  max_iteration = 15_iwp       !< max number of iterations
2419    REAL(wp),     PARAMETER  ::  guess_0       = -1.11e30_wp  !< initial guess
2420    REAL(wp) ::  x_ridder    !< current guess for clothing insulation   (clo)
2421    REAL(wp) ::  clo_lower   !< lower limit of clothing insulation      (clo)
2422    REAL(wp) ::  clo_upper   !< upper limit of clothing insulation      (clo)
2423    REAL(wp) ::  x_lower     !< lower guess for clothing insulation     (clo)
2424    REAL(wp) ::  x_upper     !< upper guess for clothing insulation     (clo)
2425    REAL(wp) ::  x_average   !< average of x_lower and x_upper          (clo)
2426    REAL(wp) ::  x_new       !< preliminary result for clothing insulation (clo)
2427    REAL(wp) ::  y_lower     !< predicted mean vote for summer clothing
2428    REAL(wp) ::  y_upper     !< predicted mean vote for winter clothing
2429    REAL(wp) ::  y_average   !< average of y_lower and y_upper
2430    REAL(wp) ::  y_new       !< preliminary result for pred. mean vote
2431    REAL(wp) ::  sroot       !< sqrt of PMV-guess
2432    INTEGER(iwp) ::  j       !< running index
2433!
2434!-- Initialise
2435    nerr    = 0_iwp
2436!
2437!-- Set pmva = 0 (comfort): Root of PMV depending on clothing insulation
2438    x_ridder    = bio_fill_value
2439    pmva        = 0._wp
2440    clo_lower   = sclo
2441    y_lower     = pmv_s
2442    clo_upper   = wclo
2443    y_upper     = pmv_w
2444    IF ( ( y_lower > 0._wp  .AND.  y_upper < 0._wp )  .OR.                     &
2445         ( y_lower < 0._wp  .AND.  y_upper > 0._wp ) )  THEN
2446       x_lower  = clo_lower
2447       x_upper  = clo_upper
2448       x_ridder = guess_0
2449
2450       DO  j = 1_iwp, max_iteration
2451          x_average = 0.5_wp * ( x_lower + x_upper )
2452          CALL fanger ( ta, tmrt, vp, ws, pair, x_average, actlev, eta,        &
2453                        y_average, top )
2454          sroot = SQRT ( y_average**2 - y_lower * y_upper )
2455          IF ( sroot == 0._wp )  THEN
2456             clo_res = x_average
2457             nerr = j
2458             RETURN
2459          ENDIF
2460          x_new = x_average + ( x_average - x_lower ) *                        &
2461                      ( SIGN ( 1._wp, y_lower - y_upper ) * y_average / sroot )
2462          IF ( ABS ( x_new - x_ridder ) <= eps )  THEN
2463             clo_res = x_ridder
2464             nerr       = j
2465             RETURN
2466          ENDIF
2467          x_ridder = x_new
2468          CALL fanger ( ta, tmrt, vp, ws, pair, x_ridder, actlev, eta,         &
2469                        y_new, top )
2470          IF ( y_new == 0._wp )  THEN
2471             clo_res = x_ridder
2472             nerr       = j
2473             RETURN
2474          ENDIF
2475          IF ( SIGN ( y_average, y_new ) /= y_average )  THEN
2476             x_lower = x_average
2477             y_lower = y_average
2478             x_upper  = x_ridder
2479             y_upper  = y_new
2480          ELSE IF ( SIGN ( y_lower, y_new ) /= y_lower )  THEN
2481             x_upper  = x_ridder
2482             y_upper  = y_new
2483          ELSE IF ( SIGN ( y_upper, y_new ) /= y_upper )  THEN
2484             x_lower = x_ridder
2485             y_lower = y_new
2486          ELSE
2487!
2488!--          Never get here in x_ridder: singularity in y
2489             nerr       = -1_iwp
2490             clo_res = x_ridder
2491             RETURN
2492          ENDIF
2493          IF ( ABS ( x_upper - x_lower ) <= eps )  THEN
2494             clo_res = x_ridder
2495             nerr       = j
2496             RETURN
2497          ENDIF
2498       ENDDO
2499!
2500!--    x_ridder exceed maximum iterations
2501       nerr       = -2_iwp
2502       clo_res = y_new
2503       RETURN
2504    ELSE IF ( y_lower == 0. )  THEN
2505       x_ridder = clo_lower
2506    ELSE IF ( y_upper == 0. )  THEN
2507       x_ridder = clo_upper
2508    ELSE
2509!
2510!--    x_ridder not bracketed by u_clo and o_clo
2511       nerr = -3_iwp
2512       clo_res = x_ridder
2513       RETURN
2514    ENDIF
2515
2516 END SUBROUTINE iso_ridder
2517
2518!------------------------------------------------------------------------------!
2519! Description:
2520! ------------
2521!> Regression relations between perceived temperature (perct) and (adjusted)
2522!> PMV. The regression presumes the Klima-Michel settings for reference
2523!> individual and reference environment.
2524!------------------------------------------------------------------------------!
2525 SUBROUTINE perct_regression( pmv, clo, perct_ij )
2526
2527    IMPLICIT NONE
2528
2529    REAL(wp), INTENT ( IN ) ::  pmv   !< Fangers predicted mean vote (dimensionless)
2530    REAL(wp), INTENT ( IN ) ::  clo   !< clothing insulation index (clo)
2531
2532    REAL(wp), INTENT ( OUT ) ::  perct_ij   !< perct (degC) corresponding to given PMV / clo
2533
2534    IF ( pmv <= -0.11_wp )  THEN
2535       perct_ij = 5.805_wp + 12.6784_wp * pmv
2536    ELSE
2537       IF ( pmv >= + 0.01_wp )  THEN
2538          perct_ij = 16.826_wp + 6.163_wp * pmv
2539       ELSE
2540          perct_ij = 21.258_wp - 9.558_wp * clo
2541       ENDIF
2542    ENDIF
2543
2544 END SUBROUTINE perct_regression
2545
2546!------------------------------------------------------------------------------!
2547! Description:
2548! ------------
2549!> FANGER.F90
2550!>
2551!> SI-VERSION: ACTLEV W m-2, DAMPFDRUCK hPa
2552!> Berechnet das aktuelle Predicted Mean Vote nach Fanger
2553!>
2554!> The case of free convection (ws < 0.1 m/s) is dealt with ws = 0.1 m/s
2555!------------------------------------------------------------------------------!
2556 SUBROUTINE fanger( ta, tmrt, pa, in_ws, pair, in_clo, actlev, eta, pmva, top )
2557
2558    IMPLICIT NONE
2559!
2560!-- Input variables of argument list:
2561    REAL(wp), INTENT ( IN ) ::  ta       !< Ambient air temperature (degC)
2562    REAL(wp), INTENT ( IN ) ::  tmrt     !< Mean radiant temperature (degC)
2563    REAL(wp), INTENT ( IN ) ::  pa       !< Water vapour pressure (hPa)
2564    REAL(wp), INTENT ( IN ) ::  pair     !< Barometric pressure (hPa) at site
2565    REAL(wp), INTENT ( IN ) ::  in_ws    !< Wind speed (m/s) 1 m above ground
2566    REAL(wp), INTENT ( IN ) ::  in_clo   !< Clothing insulation (clo)
2567    REAL(wp), INTENT ( IN ) ::  actlev   !< Individuals activity level per unit surface area (W/m2)
2568    REAL(wp), INTENT ( IN ) ::  eta      !< Individuals mechanical work efficiency (dimensionless)
2569!
2570!-- Output variables of argument list:
2571    REAL(wp), INTENT ( OUT ) ::  pmva    !< Actual Predicted Mean Vote (PMV,
2572                                         !< dimensionless) according to Fanger corresponding to meteorological
2573                                         !< (ta,tmrt,pa,ws,pair) and individual variables (clo, actlev, eta)
2574    REAL(wp), INTENT ( OUT ) ::  top     !< operative temperature (degC)
2575!
2576!-- Internal variables
2577    REAL(wp) ::  f_cl         !< Increase in surface due to clothing    (factor)
2578    REAL(wp) ::  heat_convection  !< energy loss by autocnvection       (W)
2579    REAL(wp) ::  activity     !< persons activity  (must stay == actlev, W)
2580    REAL(wp) ::  t_skin_aver  !< average skin temperature               (degree_C)
2581    REAL(wp) ::  bc           !< preliminary result storage
2582    REAL(wp) ::  cc           !< preliminary result storage
2583    REAL(wp) ::  dc           !< preliminary result storage
2584    REAL(wp) ::  ec           !< preliminary result storage
2585    REAL(wp) ::  gc           !< preliminary result storage
2586    REAL(wp) ::  t_clothing   !< clothing temperature                   (degree_C)
2587    REAL(wp) ::  hr           !< radiational heat resistence
2588    REAL(wp) ::  clo          !< clothing insulation index              (clo)
2589    REAL(wp) ::  ws           !< wind speed                             (m/s)
2590    REAL(wp) ::  z1           !< Empiric factor for the adaption of the heat
2591                              !< ballance equation to the psycho-physical scale (Equ. 40 in FANGER)
2592    REAL(wp) ::  z2           !< Water vapour diffution through the skin
2593    REAL(wp) ::  z3           !< Sweat evaporation from the skin surface
2594    REAL(wp) ::  z4           !< Loss of latent heat through respiration
2595    REAL(wp) ::  z5           !< Loss of radiational heat
2596    REAL(wp) ::  z6           !< Heat loss through forced convection
2597    INTEGER(iwp) :: i         !< running index
2598!
2599!-- Clo must be > 0. to avoid div. by 0!
2600    clo = in_clo
2601    IF ( clo <= 0._wp )  clo = .001_wp
2602!
2603!-- f_cl = Increase in surface due to clothing
2604    f_cl = 1._wp + .15_wp * clo
2605!
2606!-- Case of free convection (ws < 0.1 m/s ) not considered
2607    ws = in_ws
2608    IF ( ws < .1_wp )  THEN
2609       ws = .1_wp
2610    ENDIF
2611!
2612!-- Heat_convection = forced convection
2613    heat_convection = 12.1_wp * SQRT ( ws * pair / 1013.25_wp )
2614!
2615!-- Activity = inner heat produktion per standardized surface
2616    activity = actlev * ( 1._wp - eta )
2617!
2618!-- T_skin_aver = average skin temperature
2619    t_skin_aver = 35.7_wp - .0275_wp * activity
2620!
2621!-- Calculation of constants for evaluation below
2622    bc = .155_wp * clo * 3.96_wp * 10._wp**( -8 ) * f_cl
2623    cc = f_cl * heat_convection
2624    ec = .155_wp * clo
2625    dc = ( 1._wp + ec * cc ) / bc
2626    gc = ( t_skin_aver + bc * ( tmrt + degc_to_k )**4 + ec * cc * ta ) / bc
2627!
2628!-- Calculation of clothing surface temperature (t_clothing) based on
2629!-- Newton-approximation with air temperature as initial guess
2630    t_clothing = ta
2631    DO  i = 1, 3
2632       t_clothing = t_clothing - ( ( t_clothing + degc_to_k )**4 + t_clothing  &
2633          * dc - gc ) / ( 4._wp * ( t_clothing + degc_to_k )**3 + dc )
2634    ENDDO
2635!
2636!-- Empiric factor for the adaption of the heat ballance equation
2637!-- to the psycho-physical scale (Equ. 40 in FANGER)
2638    z1 = ( .303_wp * EXP ( -.036_wp * actlev ) + .0275_wp )
2639!
2640!-- Water vapour diffution through the skin
2641    z2 = .31_wp * ( 57.3_wp - .07_wp * activity-pa )
2642!
2643!-- Sweat evaporation from the skin surface
2644    z3 = .42_wp * ( activity - 58._wp )
2645!
2646!-- Loss of latent heat through respiration
2647    z4 = .0017_wp * actlev * ( 58.7_wp - pa ) + .0014_wp * actlev *            &
2648      ( 34._wp - ta )
2649!
2650!-- Loss of radiational heat
2651    z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + degc_to_k )**4 - ( tmrt +        &
2652       degc_to_k )**4 )
2653    IF ( ABS ( t_clothing - tmrt ) > 0._wp )  THEN
2654       hr = z5 / f_cl / ( t_clothing - tmrt )
2655    ELSE
2656       hr = 0._wp
2657    ENDIF
2658!
2659!-- Heat loss through forced convection cc*(t_clothing-TT)
2660    z6 = cc * ( t_clothing - ta )
2661!
2662!-- Predicted Mean Vote
2663    pmva = z1 * ( activity - z2 - z3 - z4 - z5 - z6 )
2664!
2665!-- Operative temperatur
2666    top = ( hr * tmrt + heat_convection * ta ) / ( hr + heat_convection )
2667
2668 END SUBROUTINE fanger
2669
2670!------------------------------------------------------------------------------!
2671! Description:
2672! ------------
2673!> For pmva > 0 and clo =0.5 the increment (deltapmv) is calculated
2674!> that converts pmva into Gagge's et al. (1986) PMV*.
2675!------------------------------------------------------------------------------!
2676 REAL(wp) FUNCTION deltapmv( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
2677
2678    IMPLICIT NONE
2679
2680!
2681!-- Input variables of argument list:
2682    REAL(wp),     INTENT ( IN )  :: pmva     !< Actual Predicted Mean Vote (PMV) according to Fanger
2683    REAL(wp),     INTENT ( IN )  :: ta       !< Ambient temperature (degC) at screen level
2684    REAL(wp),     INTENT ( IN )  :: vp       !< Water vapour pressure (hPa) at screen level
2685    REAL(wp),     INTENT ( IN )  :: svp_ta   !< Saturation water vapour pressure (hPa) at ta
2686    REAL(wp),     INTENT ( IN )  :: tmrt     !< Mean radiant temperature (degC) at screen level
2687    REAL(wp),     INTENT ( IN )  :: ws       !< Wind speed (m/s) 1 m above ground
2688
2689!
2690!-- Output variables of argument list:
2691    INTEGER(iwp), INTENT ( OUT ) :: nerr     !< Error status / quality flag
2692                                             !<  0 = o.k.
2693                                             !< -2 = pmva outside valid regression range
2694                                             !< -3 = rel. humidity set to 5 % or 95 %, respectively
2695                                             !< -4 = deltapmv set to avoid pmvs < 0
2696
2697!
2698!-- Internal variables:
2699    REAL(wp) ::  pmv          !< temp storage og predicted mean vote
2700    REAL(wp) ::  pa_p50       !< ratio actual water vapour pressure to that of relative humidity of 50 %
2701    REAL(wp) ::  pa           !< vapor pressure (hPa) with hard bounds
2702    REAL(wp) ::  apa          !< natural logarithm of pa (with hard lower border)
2703    REAL(wp) ::  dapa         !< difference of apa and pa_p50
2704    REAL(wp) ::  sqvel        !< square root of local wind velocity
2705    REAL(wp) ::  dtmrt        !< difference mean radiation to air temperature
2706    REAL(wp) ::  p10          !< lower bound for pa
2707    REAL(wp) ::  p95          !< upper bound for pa
2708    REAL(wp) ::  weight       !<
2709    REAL(wp) ::  weight2      !<
2710    REAL(wp) ::  dpmv_1       !<
2711    REAL(wp) ::  dpmv_2       !<
2712    REAL(wp) ::  pmvs         !<
2713    REAL(wp) ::  bpmv(0:7)    !<
2714    REAL(wp) ::  bpa_p50(0:7) !<
2715    REAL(wp) ::  bpa(0:7)     !<
2716    REAL(wp) ::  bapa(0:7)    !<
2717    REAL(wp) ::  bdapa(0:7)   !<
2718    REAL(wp) ::  bsqvel(0:7)  !<
2719    REAL(wp) ::  bta(0:7)     !<
2720    REAL(wp) ::  bdtmrt(0:7)  !<
2721    REAL(wp) ::  aconst(0:7)  !<
2722    INTEGER(iwp) :: nreg      !<
2723
2724    DATA bpmv     /                                                            &
2725     -0.0556602_wp, -0.1528680_wp, -0.2336104_wp, -0.2789387_wp, -0.3551048_wp,&
2726     -0.4304076_wp, -0.4884961_wp, -0.4897495_wp /
2727    DATA bpa_p50 /                                                             &
2728     -0.1607154_wp, -0.4177296_wp, -0.4120541_wp, -0.0886564_wp, +0.4285938_wp,&
2729     +0.6281256_wp, +0.5067361_wp, +0.3965169_wp /
2730    DATA bpa     /                                                             &
2731     +0.0580284_wp, +0.0836264_wp, +0.1009919_wp, +0.1020777_wp, +0.0898681_wp,&
2732     +0.0839116_wp, +0.0853258_wp, +0.0866589_wp /
2733    DATA bapa    /                                                             &
2734     -1.7838788_wp, -2.9306231_wp, -1.6350334_wp, +0.6211547_wp, +3.3918083_wp,&
2735     +5.5521025_wp, +8.4897418_wp, +16.6265851_wp /
2736    DATA bdapa   /                                                             &
2737     +1.6752720_wp, +2.7379504_wp, +1.2940526_wp, -1.0985759_wp, -3.9054732_wp,&
2738     -6.0403012_wp, -8.9437119_wp, -17.0671201_wp /
2739    DATA bsqvel  /                                                             &
2740     -0.0315598_wp, -0.0286272_wp, -0.0009228_wp, +0.0483344_wp, +0.0992366_wp,&
2741     +0.1491379_wp, +0.1951452_wp, +0.2133949_wp /
2742    DATA bta     /                                                             &
2743     +0.0953986_wp, +0.1524760_wp, +0.0564241_wp, -0.0893253_wp, -0.2398868_wp,&
2744     -0.3515237_wp, -0.5095144_wp, -0.9469258_wp /
2745    DATA bdtmrt  /                                                             &
2746     -0.0004672_wp, -0.0000514_wp, -0.0018037_wp, -0.0049440_wp, -0.0069036_wp,&
2747     -0.0075844_wp, -0.0079602_wp, -0.0089439_wp /
2748    DATA aconst  /                                                             &
2749     +1.8686215_wp, +3.4260713_wp, +2.0116185_wp, -0.7777552_wp, -4.6715853_wp,&
2750     -7.7314281_wp, -11.7602578_wp, -23.5934198_wp /
2751!
2752!-- Test for compliance with regression range
2753    IF ( pmva < -1.0_wp  .OR.  pmva > 7.0_wp )  THEN
2754       nerr = -2_iwp
2755    ELSE
2756       nerr = 0_iwp
2757    ENDIF
2758!
2759!-- Initialise classic PMV
2760    pmv  = pmva
2761!
2762!-- Water vapour pressure of air
2763    p10  = 0.05_wp * svp_ta
2764    p95  = 1.00_wp * svp_ta
2765    IF ( vp >= p10  .AND.  vp <= p95 )  THEN
2766       pa = vp
2767    ELSE
2768       nerr = -3_iwp
2769       IF ( vp < p10 )  THEN
2770!
2771!--       Due to conditions of regression: r.H. >= 5 %
2772          pa = p10
2773       ELSE
2774!
2775!--       Due to conditions of regression: r.H. <= 95 %
2776          pa = p95
2777       ENDIF
2778    ENDIF
2779    IF ( pa > 0._wp )  THEN
2780!
2781!--    Natural logarithm of pa
2782       apa = LOG ( pa )
2783    ELSE
2784       apa = -5._wp
2785    ENDIF
2786!
2787!-- Ratio actual water vapour pressure to that of a r.H. of 50 %
2788    pa_p50   = 0.5_wp * svp_ta
2789    IF ( pa_p50 > 0._wp  .AND.  pa > 0._wp )  THEN
2790       dapa   = apa - LOG ( pa_p50 )
2791       pa_p50 = pa / pa_p50
2792    ELSE
2793       dapa   = -5._wp
2794       pa_p50 = 0._wp
2795    ENDIF
2796!
2797!-- Square root of wind velocity
2798    IF ( ws >= 0._wp )  THEN
2799       sqvel = SQRT ( ws )
2800    ELSE
2801       sqvel = 0._wp
2802    ENDIF
2803!
2804!-- Difference mean radiation to air temperature
2805    dtmrt = tmrt - ta
2806!
2807!-- Select the valid regression coefficients
2808    nreg = INT ( pmv )
2809    IF ( nreg < 0_iwp )  THEN
2810!
2811!--    value of the FUNCTION in the case pmv <= -1
2812       deltapmv = 0._wp
2813       RETURN
2814    ENDIF
2815    weight = MOD ( pmv, 1._wp )
2816    IF ( weight < 0._wp )  weight = 0._wp
2817    IF ( nreg > 5_iwp )  THEN
2818       ! nreg=6
2819       nreg  = 5_iwp
2820       weight   = pmv - 5._wp
2821       weight2  = pmv - 6._wp
2822       IF ( weight2 > 0_iwp )  THEN
2823          weight = ( weight - weight2 ) / weight
2824       ENDIF
2825    ENDIF
2826!
2827!-- Regression valid for 0. <= pmv <= 6.
2828    dpmv_1 =                                                                   &
2829       + bpa(nreg) * pa                                                        &
2830       + bpmv(nreg) * pmv                                                      &
2831       + bapa(nreg) * apa                                                      &
2832       + bta(nreg) * ta                                                        &
2833       + bdtmrt(nreg) * dtmrt                                                  &
2834       + bdapa(nreg) * dapa                                                    &
2835       + bsqvel(nreg) * sqvel                                                  &
2836       + bpa_p50(nreg) * pa_p50                                                &
2837       + aconst(nreg)
2838
2839    dpmv_2 = 0._wp
2840    IF ( nreg < 6_iwp )  THEN
2841       dpmv_2 =                                                                &
2842          + bpa(nreg+1_iwp)     * pa                                           &
2843          + bpmv(nreg+1_iwp)    * pmv                                          &
2844          + bapa(nreg+1_iwp)    * apa                                          &
2845          + bta(nreg+1_iwp)     * ta                                           &
2846          + bdtmrt(nreg+1_iwp)  * dtmrt                                        &
2847          + bdapa(nreg+1_iwp)   * dapa                                         &
2848          + bsqvel(nreg+1_iwp)  * sqvel                                        &
2849          + bpa_p50(nreg+1_iwp) * pa_p50                                       &
2850          + aconst(nreg+1_iwp)
2851    ENDIF
2852!
2853!-- Calculate pmv modification
2854    deltapmv = ( 1._wp - weight ) * dpmv_1 + weight * dpmv_2
2855    pmvs = pmva + deltapmv
2856    IF ( ( pmvs ) < 0._wp )  THEN
2857!
2858!--    Prevent negative pmv* due to problems with clothing insulation
2859       nerr = -4_iwp
2860       IF ( pmvs > -0.11_wp )  THEN
2861!
2862!--       Threshold from FUNCTION perct_regression for winter clothing insulation
2863          deltapmv = deltapmv + 0.11_wp
2864       ELSE
2865!
2866!--       Set pmvs to "0" for compliance with summer clothing insulation
2867          deltapmv = -1._wp * pmva
2868       ENDIF
2869    ENDIF
2870
2871 END FUNCTION deltapmv
2872
2873!------------------------------------------------------------------------------!
2874! Description:
2875! ------------
2876!> The subroutine "calc_sultr" returns a threshold value to perceived
2877!> temperature allowing to decide whether the actual perceived temperature
2878!> is linked to perecption of sultriness. The threshold values depends
2879!> on the Fanger's classical PMV, expressed here as perceived temperature
2880!> perct.
2881!------------------------------------------------------------------------------!
2882 SUBROUTINE calc_sultr( perct_ij, dperctm, dperctstd, sultr_res )
2883
2884    IMPLICIT NONE
2885!
2886!-- Input of the argument list:
2887    REAL(wp), INTENT ( IN )  ::  perct_ij      !< Classical perceived temperature: Base is Fanger's PMV
2888!
2889!-- Additional output variables of argument list:
2890    REAL(wp), INTENT ( OUT ) ::  dperctm    !< Mean deviation perct (classical gt) to gt* (rational gt
2891                                            !< calculated based on Gagge's rational PMV*)
2892    REAL(wp), INTENT ( OUT ) ::  dperctstd  !< dperctm plus its standard deviation times a factor
2893                                            !< determining the significance to perceive sultriness
2894    REAL(wp), INTENT ( OUT ) ::  sultr_res
2895!
2896!-- Types of coefficients mean deviation: third order polynomial
2897    REAL(wp), PARAMETER ::  dperctka = +7.5776086_wp
2898    REAL(wp), PARAMETER ::  dperctkb = -0.740603_wp
2899    REAL(wp), PARAMETER ::  dperctkc = +0.0213324_wp
2900    REAL(wp), PARAMETER ::  dperctkd = -0.00027797237_wp
2901!
2902!-- Types of coefficients mean deviation plus standard deviation
2903!-- regression coefficients: third order polynomial
2904    REAL(wp), PARAMETER ::  dperctsa = +0.0268918_wp
2905    REAL(wp), PARAMETER ::  dperctsb = +0.0465957_wp
2906    REAL(wp), PARAMETER ::  dperctsc = -0.00054709752_wp
2907    REAL(wp), PARAMETER ::  dperctsd = +0.0000063714823_wp
2908!
2909!-- Factor to mean standard deviation defining SIGNificance for
2910!-- sultriness
2911    REAL(wp), PARAMETER :: faktor = 1._wp
2912!
2913!-- Initialise
2914    sultr_res = +99._wp
2915    dperctm   = 0._wp
2916    dperctstd = 999999._wp
2917
2918    IF ( perct_ij < 16.826_wp  .OR.  perct_ij > 56._wp )  THEN
2919!
2920!--    Unallowed value of classical perct!
2921       RETURN
2922    ENDIF
2923!
2924!-- Mean deviation dependent on perct
2925    dperctm = dperctka + dperctkb * perct_ij + dperctkc * perct_ij**2._wp +    &
2926       dperctkd * perct_ij**3._wp
2927!
2928!-- Mean deviation plus its standard deviation
2929    dperctstd = dperctsa + dperctsb * perct_ij + dperctsc * perct_ij**2._wp +  &
2930       dperctsd * perct_ij**3._wp
2931!
2932!-- Value of the FUNCTION
2933    sultr_res = dperctm + faktor * dperctstd
2934    IF ( ABS ( sultr_res ) > 99._wp )  sultr_res = +99._wp
2935
2936 END SUBROUTINE calc_sultr
2937
2938!------------------------------------------------------------------------------!
2939! Description:
2940! ------------
2941!> Multiple linear regression to calculate an increment delta_cold,
2942!> to adjust Fanger's classical PMV (pmva) by Gagge's 2 node model,
2943!> applying Fanger's convective heat transfer coefficient, hcf.
2944!> Wind velocitiy of the reference environment is 0.10 m/s
2945!------------------------------------------------------------------------------!
2946 SUBROUTINE dpmv_cold( pmva, ta, ws, tmrt, nerr, dpmv_cold_res )
2947
2948    IMPLICIT NONE
2949!
2950!-- Type of input arguments
2951    REAL(wp), INTENT ( IN ) ::  pmva   !< Fanger's classical predicted mean vote
2952    REAL(wp), INTENT ( IN ) ::  ta     !< Air temperature 2 m above ground (degC)
2953    REAL(wp), INTENT ( IN ) ::  ws     !< Relative wind velocity 1 m above ground (m/s)
2954    REAL(wp), INTENT ( IN ) ::  tmrt   !< Mean radiant temperature (degC)
2955!
2956!-- Type of output argument
2957    INTEGER(iwp), INTENT ( OUT ) ::  nerr !< Error indicator: 0 = o.k., +1 = denominator for intersection = 0
2958    REAL(wp),     INTENT ( OUT ) ::  dpmv_cold_res    !< Increment to adjust pmva according to the results of Gagge's
2959                                                      !< 2 node model depending on the input
2960!
2961!-- Type of program variables
2962    REAL(wp) ::  delta_cold(3)
2963    REAL(wp) ::  pmv_cross(2)
2964    REAL(wp) ::  reg_a(3)
2965    REAL(wp) ::  coeff(3,5)
2966    REAL(wp) ::  r_nenner
2967    REAL(wp) ::  pmvc
2968    REAL(wp) ::  dtmrt
2969    REAL(wp) ::  sqrt_ws
2970    INTEGER(iwp) ::  i
2971    INTEGER(iwp) ::  i_bin
2972!
2973!-- Coefficient of the 3 regression lines
2974!     1:const   2:*pmvc  3:*ta      4:*sqrt_ws  5:*dtmrt
2975    coeff(1,1:5) =                                                             &
2976       (/ +0.161_wp, +0.130_wp, -1.125E-03_wp, +1.106E-03_wp, -4.570E-04_wp /)
2977    coeff(2,1:5) =                                                             &
2978       (/ 0.795_wp, 0.713_wp, -8.880E-03_wp, -1.803E-03_wp, -2.816E-03_wp/)
2979    coeff(3,1:5) =                                                             &
2980       (/ +0.05761_wp, +0.458_wp, -1.829E-02_wp, -5.577E-03_wp, -1.970E-03_wp /)
2981!
2982!-- Initialise
2983    nerr           = 0_iwp
2984    dpmv_cold_res  = 0._wp
2985    pmvc           = pmva
2986    dtmrt          = tmrt - ta
2987    sqrt_ws        = ws
2988    IF ( sqrt_ws < 0.10_wp )  THEN
2989       sqrt_ws = 0.10_wp
2990    ELSE
2991       sqrt_ws = SQRT ( sqrt_ws )
2992    ENDIF
2993
2994    DO  i = 1, 3
2995       delta_cold(i) = 0._wp
2996       IF ( i < 3 )  THEN
2997          pmv_cross(i) = pmvc
2998       ENDIF
2999    ENDDO
3000
3001    DO  i = 1, 3
3002!
3003!--    Regression constant for given meteorological variables
3004       reg_a(i) = coeff(i,1) + coeff(i,3) * ta + coeff(i,4) *                  &
3005                  sqrt_ws + coeff(i,5)*dtmrt
3006       delta_cold(i) = reg_a(i) + coeff(i,2) * pmvc
3007    ENDDO
3008!
3009!-- Intersection points of regression lines in terms of Fanger's PMV
3010    DO  i = 1, 2
3011       r_nenner = coeff(i,2) - coeff(i+1,2)
3012       IF ( ABS ( r_nenner ) > 0.00001_wp )  THEN
3013          pmv_cross(i) = ( reg_a(i+1) - reg_a(i) ) / r_nenner
3014       ELSE
3015          nerr = 1_iwp
3016          RETURN
3017       ENDIF
3018    ENDDO
3019
3020    i_bin = 3
3021    DO  i = 1, 2
3022       IF ( pmva > pmv_cross(i) )  THEN
3023          i_bin = i
3024          EXIT
3025       ENDIF
3026    ENDDO
3027!
3028!-- Adjust to operative temperature scaled according
3029!-- to classical PMV (Fanger)
3030    dpmv_cold_res = delta_cold(i_bin) - dpmv_adj(pmva)
3031
3032 END SUBROUTINE dpmv_cold
3033
3034!------------------------------------------------------------------------------!
3035! Description:
3036! ------------
3037!> Calculates the summand dpmv_adj adjusting to the operative temperature
3038!> scaled according to classical PMV (Fanger)
3039!> Reference environment: v_1m = 0.10 m/s, dTMRT = 0 K, r.h. = 50 %
3040!------------------------------------------------------------------------------!
3041 REAL(wp) FUNCTION dpmv_adj( pmva )
3042
3043    IMPLICIT NONE
3044
3045    REAL(wp), INTENT ( IN ) ::  pmva
3046    INTEGER(iwp), PARAMETER ::  n_bin = 3
3047    INTEGER(iwp), PARAMETER ::  n_ca = 0
3048    INTEGER(iwp), PARAMETER ::  n_ce = 3
3049    REAL(wp), dimension (n_bin, n_ca:n_ce) ::  coef
3050
3051    REAL(wp)      ::  pmv
3052    INTEGER(iwp)  ::  i, i_bin
3053!
3054!--                                 range_1        range_2        range_3
3055    DATA (coef(i,0), i = 1, n_bin) /0.0941540_wp, -0.1506620_wp, -0.0871439_wp/
3056    DATA (coef(i,1), i = 1, n_bin) /0.0783162_wp, -1.0612651_wp,  0.1695040_wp/
3057    DATA (coef(i,2), i = 1, n_bin) /0.1350144_wp, -1.0049144_wp, -0.0167627_wp/
3058    DATA (coef(i,3), i = 1, n_bin) /0.1104037_wp, -0.2005277_wp, -0.0003230_wp/
3059
3060    IF ( pmva <= -2.1226_wp )  THEN
3061       i_bin = 3_iwp
3062    ELSE IF ( pmva <= -1.28_wp )  THEN
3063       i_bin = 2_iwp
3064    ELSE
3065       i_bin = 1_iwp
3066    ENDIF
3067
3068    dpmv_adj   = coef(i_bin,n_ca)
3069    pmv        = 1._wp
3070
3071    DO  i = n_ca + 1, n_ce
3072       pmv      = pmv * pmva
3073       dpmv_adj = dpmv_adj + coef(i_bin,i) * pmv
3074    ENDDO
3075    RETURN
3076 END FUNCTION dpmv_adj
3077
3078!------------------------------------------------------------------------------!
3079! Description:
3080! ------------
3081!> Based on perceived temperature (perct) as input, ireq_neutral determines
3082!> the required clothing insulation (clo) for thermally neutral conditions
3083!> (neither body cooling nor body heating). It is related to the Klima-
3084!> Michel activity level (134.682 W/m2). IREQ_neutral is only defined
3085!> for perct < 10 (degC)
3086!------------------------------------------------------------------------------!
3087 REAL(wp) FUNCTION ireq_neutral( perct_ij, ireq_minimal, nerr )
3088
3089    IMPLICIT NONE
3090!
3091!-- Type declaration of arguments
3092    REAL(wp),     INTENT ( IN )  ::  perct_ij
3093    REAL(wp),     INTENT ( OUT ) ::  ireq_minimal
3094    INTEGER(iwp), INTENT ( OUT ) ::  nerr
3095!
3096!-- Type declaration for internal varables
3097    REAL(wp)                     ::  top02
3098!
3099!-- Initialise
3100    nerr = 0_iwp
3101!
3102!-- Convert perceived temperature from basis 0.1 m/s to basis 0.2 m/s
3103    top02 = 1.8788_wp + 0.9296_wp * perct_ij
3104!
3105!-- IREQ neutral conditions (thermal comfort)
3106    ireq_neutral = 1.62_wp - 0.0564_wp * top02
3107!
3108!-- Regression only defined for perct <= 10 (degC)
3109    IF ( ireq_neutral < 0.5_wp )  THEN
3110       IF ( ireq_neutral < 0.48_wp )  THEN
3111          nerr = 1_iwp
3112       ENDIF
3113       ireq_neutral = 0.5_wp
3114    ENDIF
3115!
3116!-- Minimal required clothing insulation: maximal acceptable body cooling
3117    ireq_minimal = 1.26_wp - 0.0588_wp * top02
3118    IF ( nerr > 0_iwp )  THEN
3119       ireq_minimal = ireq_neutral
3120    ENDIF
3121
3122    RETURN
3123 END FUNCTION ireq_neutral
3124
3125
3126!------------------------------------------------------------------------------!
3127! Description:
3128! ------------
3129!> The SUBROUTINE surface area calculates the surface area of the individual
3130!> according to its height (m), weight (kg), and age (y)
3131!------------------------------------------------------------------------------!
3132 SUBROUTINE surface_area ( height_cm, weight, age, surf )
3133
3134    IMPLICIT NONE
3135
3136    REAL(wp)    , INTENT(in)  ::  weight
3137    REAL(wp)    , INTENT(in)  ::  height_cm
3138    INTEGER(iwp), INTENT(in)  ::  age
3139    REAL(wp)    , INTENT(out) ::  surf
3140    REAL(wp)                  ::  height
3141
3142    height = height_cm * 100._wp
3143!
3144!-- According to Gehan-George, for children
3145    IF ( age < 19_iwp )  THEN
3146       IF ( age < 5_iwp )  THEN
3147          surf = 0.02667_wp * height**0.42246_wp * weight**0.51456_wp
3148          RETURN
3149       ENDIF
3150       surf = 0.03050_wp * height**0.35129_wp * weight**0.54375_wp
3151       RETURN
3152    ENDIF
3153!
3154!-- DuBois D, DuBois EF: A formula to estimate the approximate surface area if
3155!-- height and weight be known. In: Arch. Int. Med.. 17, 1916, S. 863?871.
3156    surf = 0.007184_wp * height**0.725_wp * weight**0.425_wp
3157    RETURN
3158
3159 END SUBROUTINE surface_area
3160
3161!------------------------------------------------------------------------------!
3162! Description:
3163! ------------
3164!> The SUBROUTINE persdat calculates
3165!>  - the total internal energy production = metabolic + workload,
3166!>  - the total internal energy production for a standardized surface (actlev)
3167!>  - the DuBois - area (a_surf [m2])
3168!> from
3169!>  - the persons age (years),
3170!>  - weight (kg),
3171!>  - height (m),
3172!>  - sex (1 = male, 2 = female),
3173!>  - work load (W)
3174!> for a sample human.
3175!------------------------------------------------------------------------------!
3176 SUBROUTINE persdat ( age, weight, height, sex, work, a_surf, actlev )
3177
3178    IMPLICIT NONE
3179
3180    REAL(wp), INTENT(in) ::  age
3181    REAL(wp), INTENT(in) ::  weight
3182    REAL(wp), INTENT(in) ::  height
3183    REAL(wp), INTENT(in) ::  work
3184    INTEGER(iwp), INTENT(in) ::  sex
3185    REAL(wp), INTENT(out) ::  actlev
3186    REAL(wp) ::  a_surf
3187    REAL(wp) ::  energy_prod
3188    REAL(wp) ::  s
3189    REAL(wp) ::  factor
3190    REAL(wp) ::  basic_heat_prod
3191
3192    CALL surface_area ( height, weight, INT( age ), a_surf )
3193    s = height * 100._wp / ( weight**( 1._wp / 3._wp ) )
3194    factor = 1._wp + .004_wp  * ( 30._wp - age )
3195    basic_heat_prod = 0.
3196    IF ( sex == 1_iwp )  THEN
3197       basic_heat_prod = 3.45_wp * weight**( 3._wp / 4._wp ) * ( factor +      &
3198                     .01_wp  * ( s - 43.4_wp ) )
3199    ELSE IF ( sex == 2_iwp )  THEN
3200       basic_heat_prod = 3.19_wp * weight**( 3._wp / 4._wp ) * ( factor +      &
3201                    .018_wp * ( s - 42.1_wp ) )
3202    ENDIF
3203
3204    energy_prod = work + basic_heat_prod
3205    actlev = energy_prod / a_surf
3206
3207 END SUBROUTINE persdat
3208
3209
3210!------------------------------------------------------------------------------!
3211! Description:
3212! ------------
3213!> SUBROUTINE ipt_init
3214!> initializes the instationary perceived temperature
3215!------------------------------------------------------------------------------!
3216
3217 SUBROUTINE ipt_init (age, weight, height, sex, work, actlev, clo,             &
3218     ta, vp, ws, tmrt, pair, dt, storage, t_clothing,                          &
3219     ipt )
3220
3221    IMPLICIT NONE
3222!
3223!-- Input parameters
3224    REAL(wp), INTENT(in) ::  age        !< Persons age          (years)
3225    REAL(wp), INTENT(in) ::  weight     !< Persons weight       (kg)
3226    REAL(wp), INTENT(in) ::  height     !< Persons height       (m)
3227    REAL(wp), INTENT(in) ::  work       !< Current workload     (W)
3228    REAL(wp), INTENT(in) ::  ta         !< Air Temperature      (degree_C)
3229    REAL(wp), INTENT(in) ::  vp         !< Vapor pressure       (hPa)
3230    REAL(wp), INTENT(in) ::  ws         !< Wind speed in approx. 1.1m (m/s)
3231    REAL(wp), INTENT(in) ::  tmrt       !< Mean radiant temperature   (degree_C)
3232    REAL(wp), INTENT(in) ::  pair       !< Air pressure         (hPa)
3233    REAL(wp), INTENT(in) ::  dt         !< Timestep             (s)
3234    INTEGER(iwp), INTENT(in)  :: sex    !< Persons sex (1 = male, 2 = female)
3235!
3236!-- Output parameters
3237    REAL(wp), INTENT(out) ::  actlev
3238    REAL(wp), INTENT(out) ::  clo
3239    REAL(wp), INTENT(out) ::  storage
3240    REAL(wp), INTENT(out) ::  t_clothing
3241    REAL(wp), INTENT(out) ::  ipt
3242!
3243!-- Internal variables
3244    REAL(wp), PARAMETER :: eps = 0.0005_wp
3245    REAL(wp), PARAMETER :: eta = 0._wp
3246    REAL(wp) ::  sclo
3247    REAL(wp) ::  wclo
3248    REAL(wp) ::  d_pmv
3249    REAL(wp) ::  svp_ta
3250    REAL(wp) ::  sult_lim
3251    REAL(wp) ::  dgtcm
3252    REAL(wp) ::  dgtcstd
3253    REAL(wp) ::  clon
3254    REAL(wp) ::  ireq_minimal
3255!     REAL(wp) ::  clo_fanger
3256    REAL(wp) ::  pmv_w
3257    REAL(wp) ::  pmv_s
3258    REAL(wp) ::  pmva
3259    REAL(wp) ::  ptc
3260    REAL(wp) ::  d_std
3261    REAL(wp) ::  pmvs
3262    REAL(wp) ::  top
3263    REAL(wp) ::  a_surf
3264!     REAL(wp) ::  acti
3265    INTEGER(iwp) ::  ncount
3266    INTEGER(iwp) ::  nerr_cold
3267    INTEGER(iwp) ::  nerr
3268
3269    LOGICAL ::  sultrieness
3270
3271    storage = 0._wp
3272    CALL persdat ( age, weight, height, sex, work, a_surf, actlev )
3273!
3274!-- Initialise
3275    t_clothing = bio_fill_value
3276    ipt        = bio_fill_value
3277    nerr       = 0_wp
3278    ncount     = 0_wp
3279    sultrieness    = .FALSE.
3280!
3281!-- Tresholds: clothing insulation (account for model inaccuracies)
3282!-- Summer clothing
3283    sclo     = 0.44453_wp
3284!-- Winter clothing
3285    wclo     = 1.76267_wp
3286!
3287!-- Decision: firstly calculate for winter or summer clothing
3288    IF ( ta <= 10._wp )  THEN
3289!
3290!--    First guess: winter clothing insulation: cold stress
3291       clo = wclo
3292       t_clothing = bio_fill_value  ! force initial run
3293       CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,         &
3294          t_clothing, storage, dt, pmva )
3295       pmv_w = pmva
3296
3297       IF ( pmva > 0._wp )  THEN
3298!
3299!--       Case summer clothing insulation: heat load ?           
3300          clo = sclo
3301          t_clothing = bio_fill_value  ! force initial run
3302          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,      &
3303             t_clothing, storage, dt, pmva )
3304          pmv_s = pmva
3305          IF ( pmva <= 0._wp )  THEN
3306!
3307!--          Case: comfort achievable by varying clothing insulation
3308!--          between winter and summer set values
3309             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta , sclo,     &
3310                            pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
3311             IF ( ncount < 0_iwp )  THEN
3312                nerr = -1_iwp
3313                RETURN
3314             ENDIF
3315          ELSE IF ( pmva > 0.06_wp )  THEN
3316             clo = 0.5_wp
3317             t_clothing = bio_fill_value
3318             CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,   &
3319                t_clothing, storage, dt, pmva )
3320          ENDIF
3321       ELSE IF ( pmva < -0.11_wp )  THEN
3322          clo = 1.75_wp
3323          t_clothing = bio_fill_value
3324          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,      &
3325             t_clothing, storage, dt, pmva )
3326       ENDIF
3327
3328    ELSE
3329!
3330!--    First guess: summer clothing insulation: heat load
3331       clo = sclo
3332       t_clothing = bio_fill_value
3333       CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,         &
3334          t_clothing, storage, dt, pmva )
3335       pmv_s = pmva
3336
3337       IF ( pmva < 0._wp )  THEN
3338!
3339!--       Case winter clothing insulation: cold stress ?
3340          clo = wclo
3341          t_clothing = bio_fill_value
3342          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,      &
3343             t_clothing, storage, dt, pmva )
3344          pmv_w = pmva
3345
3346          IF ( pmva >= 0._wp )  THEN
3347!
3348!--          Case: comfort achievable by varying clothing insulation
3349!--          between winter and summer set values
3350             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
3351                               pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
3352             IF ( ncount < 0_wp )  THEN
3353                nerr = -1_iwp
3354                RETURN
3355             ENDIF
3356          ELSE IF ( pmva < -0.11_wp )  THEN
3357             clo = 1.75_wp
3358             t_clothing = bio_fill_value
3359             CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,   &
3360                t_clothing, storage, dt, pmva )
3361          ENDIF
3362       ELSE IF ( pmva > 0.06_wp )  THEN
3363          clo = 0.5_wp
3364          t_clothing = bio_fill_value
3365          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,      &
3366             t_clothing, storage, dt, pmva )
3367       ENDIF
3368
3369    ENDIF
3370!
3371!-- Determine perceived temperature by regression equation + adjustments
3372    pmvs = pmva
3373    CALL perct_regression ( pmva, clo, ipt )
3374    ptc = ipt
3375    IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp )  THEN
3376!
3377!--    Adjust for cold conditions according to Gagge 1986
3378       CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
3379       IF ( nerr_cold > 0_iwp )  nerr = -5_iwp
3380       pmvs = pmva - d_pmv
3381       IF ( pmvs > -0.11_wp )  THEN
3382          d_pmv  = 0._wp
3383          pmvs   = -0.11_wp
3384       ENDIF
3385       CALL perct_regression ( pmvs, clo, ipt )
3386    ENDIF
3387!     clo_fanger = clo
3388    clon = clo
3389    IF ( clo > 0.5_wp  .AND.  ipt <= 8.73_wp )  THEN
3390!
3391!--    Required clothing insulation (ireq) is exclusively defined for
3392!--    operative temperatures (top) less 10 (C) for a
3393!--    reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
3394       clon = ireq_neutral ( ipt, ireq_minimal, nerr )
3395       clo = clon
3396    ENDIF
3397    CALL calc_sultr ( ptc, dgtcm, dgtcstd, sult_lim )
3398    sultrieness    = .FALSE.
3399    d_std      = -99._wp
3400    IF ( pmva > 0.06_wp  .AND.  clo <= 0.5_wp )  THEN
3401!
3402!--    Adjust for warm/humid conditions according to Gagge 1986
3403       CALL saturation_vapor_pressure ( ta, svp_ta )
3404       d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
3405       pmvs   = pmva + d_pmv
3406       CALL perct_regression ( pmvs, clo, ipt )
3407       IF ( sult_lim < 99._wp )  THEN
3408          IF ( (ipt - ptc) > sult_lim )  sultrieness = .TRUE.
3409       ENDIF
3410    ENDIF
3411
3412 
3413 END SUBROUTINE ipt_init
3414 
3415!------------------------------------------------------------------------------!
3416! Description:
3417! ------------
3418!> SUBROUTINE ipt_cycle
3419!> Calculates one timestep for the instationary version of perceived
3420!> temperature (iPT, degree_C) for
3421!>  - standard measured/predicted meteorological values and TMRT
3422!>    as input;
3423!>  - regressions for determination of PT;
3424!>  - adjustment to Gagge's PMV* (2-node-model, 1986) as base of PT
3425!>    under warm/humid conditions (Icl= 0.50 clo) and under cold
3426!>    conditions (Icl= 1.75 clo)
3427!>
3428!------------------------------------------------------------------------------!
3429 SUBROUTINE ipt_cycle( ta, vp, ws, tmrt, pair, dt, storage, t_clothing, clo,   &
3430     actlev, work, ipt )
3431
3432    IMPLICIT NONE
3433!
3434!-- Type of input of the argument list
3435    REAL(wp), INTENT ( IN )  ::  ta      !< Air temperature             (degree_C)
3436    REAL(wp), INTENT ( IN )  ::  vp      !< Vapor pressure              (hPa)
3437    REAL(wp), INTENT ( IN )  ::  tmrt    !< Mean radiant temperature    (degree_C)
3438    REAL(wp), INTENT ( IN )  ::  ws      !< Wind speed                  (m/s)
3439    REAL(wp), INTENT ( IN )  ::  pair    !< Air pressure                (hPa)
3440    REAL(wp), INTENT ( IN )  ::  dt      !< Timestep                    (s)
3441    REAL(wp), INTENT ( IN )  ::  clo     !< Clothing index              (no dim)
3442    REAL(wp), INTENT ( IN )  ::  actlev  !< Internal heat production    (W)
3443    REAL(wp), INTENT ( IN )  ::  work    !< Mechanical work load        (W)
3444!
3445!-- In and output parameters
3446    REAL(wp), INTENT (INOUT) ::  storage     !< Heat storage            (W)
3447    REAL(wp), INTENT (INOUT) ::  t_clothing  !< Clothig temperature     (degree_C)
3448!
3449!-- Type of output of the argument list
3450    REAL(wp), INTENT ( OUT ) ::  ipt  !< Instationary perceived temperature (degree_C)
3451!
3452!-- Type of internal variables
3453    REAL(wp) ::  d_pmv
3454    REAL(wp) ::  svp_ta
3455    REAL(wp) ::  sult_lim
3456    REAL(wp) ::  dgtcm
3457    REAL(wp) ::  dgtcstd
3458    REAL(wp) ::  pmva
3459    REAL(wp) ::  ptc
3460    REAL(wp) ::  d_std
3461    REAL(wp) ::  pmvs
3462    INTEGER(iwp) ::  nerr_cold
3463    INTEGER(iwp) ::  nerr
3464
3465    LOGICAL ::  sultrieness
3466!
3467!-- Initialise
3468    ipt = bio_fill_value
3469
3470    nerr     = 0_iwp
3471    sultrieness  = .FALSE.
3472!
3473!-- Determine pmv_adjusted for current conditions
3474    CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,            &
3475       t_clothing, storage, dt, pmva )
3476!
3477!-- Determine perceived temperature by regression equation + adjustments
3478    CALL perct_regression ( pmva, clo, ipt )
3479
3480    IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp )  THEN
3481!
3482!--    Adjust for cold conditions according to Gagge 1986
3483       CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
3484       IF ( nerr_cold > 0_iwp )  nerr = -5_iwp
3485       pmvs = pmva - d_pmv
3486       IF ( pmvs > -0.11_wp )  THEN
3487          d_pmv  = 0._wp
3488          pmvs   = -0.11_wp
3489       ENDIF
3490       CALL perct_regression ( pmvs, clo, ipt )
3491    ENDIF
3492!
3493!-- Consider sultriness if appropriate
3494    ptc = ipt
3495    CALL calc_sultr ( ptc, dgtcm, dgtcstd, sult_lim )
3496    sultrieness    = .FALSE.
3497    d_std      = -99._wp
3498    IF ( pmva > 0.06_wp  .AND.  clo <= 0.5_wp )  THEN
3499!
3500!--    Adjust for warm/humid conditions according to Gagge 1986
3501       CALL saturation_vapor_pressure ( ta, svp_ta )
3502       d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
3503       pmvs   = pmva + d_pmv
3504       CALL perct_regression ( pmvs, clo, ipt )
3505       IF ( sult_lim < 99._wp )  THEN
3506          IF ( (ipt - ptc) > sult_lim )  sultrieness = .TRUE.
3507       ENDIF
3508    ENDIF
3509
3510 END SUBROUTINE ipt_cycle
3511
3512!------------------------------------------------------------------------------!
3513! Description:
3514! ------------
3515!> SUBROUTINE fanger_s calculates the
3516!> actual Predicted Mean Vote (dimensionless) according
3517!> to Fanger corresponding to meteorological (ta,tmrt,pa,ws,pair)
3518!> and individual variables (clo, actlev, eta) considering a storage
3519!> and clothing temperature for a given timestep.
3520!------------------------------------------------------------------------------!
3521 SUBROUTINE fanger_s_acti( ta, tmrt, pa, in_ws, pair, in_clo, actlev,          &
3522    activity, t_cloth, s, dt, pmva )
3523
3524    IMPLICIT NONE
3525!
3526!--  Input argument types
3527    REAL(wp), INTENT ( IN )  ::  ta       !< Air temperature          (degree_C)
3528    REAL(wp), INTENT ( IN )  ::  tmrt     !< Mean radiant temperature (degree_C)
3529    REAL(wp), INTENT ( IN )  ::  pa       !< Vapour pressure          (hPa)
3530    REAL(wp), INTENT ( IN )  ::  pair     !< Air pressure             (hPa)
3531    REAL(wp), INTENT ( IN )  ::  in_ws    !< Wind speed               (m/s)
3532    REAL(wp), INTENT ( IN )  ::  actlev   !< Metabolic + work energy  (W/m²)
3533    REAL(wp), INTENT ( IN )  ::  dt       !< Timestep                 (s)
3534    REAL(wp), INTENT ( IN )  ::  activity !< Work load                (W/m²)
3535    REAL(wp), INTENT ( IN )  ::  in_clo   !< Clothing index (clo)     (no dim)
3536!
3537!-- Output argument types
3538    REAL(wp), INTENT ( OUT ) ::  pmva  !< actual Predicted Mean Vote (no dim)
3539
3540    REAL(wp), INTENT (INOUT) ::  s  !< storage var. of energy balance (W/m2)
3541    REAL(wp), INTENT (INOUT) ::  t_cloth  !< clothing temperature (degree_C)
3542!
3543!-- Internal variables
3544    REAL(wp), PARAMETER  ::  time_equil = 7200._wp
3545
3546    REAL(wp) ::  f_cl         !< Increase in surface due to clothing    (factor)
3547    REAL(wp) ::  heat_convection  !< energy loss by autocnvection       (W)
3548    REAL(wp) ::  t_skin_aver  !< average skin temperature               (degree_C)
3549    REAL(wp) ::  bc           !< preliminary result storage
3550    REAL(wp) ::  cc           !< preliminary result storage
3551    REAL(wp) ::  dc           !< preliminary result storage
3552    REAL(wp) ::  ec           !< preliminary result storage
3553    REAL(wp) ::  gc           !< preliminary result storage
3554    REAL(wp) ::  t_clothing   !< clothing temperature                   (degree_C)
3555!     REAL(wp) ::  hr           !< radiational heat resistence
3556    REAL(wp) ::  clo          !< clothing insulation index              (clo)
3557    REAL(wp) ::  ws           !< wind speed                             (m/s)
3558    REAL(wp) ::  z1           !< Empiric factor for the adaption of the heat
3559                              !< ballance equation to the psycho-physical scale (Equ. 40 in FANGER)
3560    REAL(wp) ::  z2           !< Water vapour diffution through the skin
3561    REAL(wp) ::  z3           !< Sweat evaporation from the skin surface
3562    REAL(wp) ::  z4           !< Loss of latent heat through respiration
3563    REAL(wp) ::  z5           !< Loss of radiational heat
3564    REAL(wp) ::  z6           !< Heat loss through forced convection
3565    REAL(wp) ::  en           !< Energy ballance                        (W)
3566    REAL(wp) ::  d_s          !< Storage delta                          (W)
3567    REAL(wp) ::  adjustrate   !< Max storage adjustment rate
3568    REAL(wp) ::  adjustrate_cloth  !< max clothing temp. adjustment rate
3569
3570    INTEGER(iwp) :: i         !< running index
3571    INTEGER(iwp) ::  niter    !< Running index
3572
3573!
3574!-- Clo must be > 0. to avoid div. by 0!
3575    clo = in_clo
3576    IF ( clo < 001._wp )  clo = .001_wp
3577!
3578!-- Increase in surface due to clothing
3579    f_cl = 1._wp + .15_wp * clo
3580!
3581!-- Case of free convection (ws < 0.1 m/s ) not considered
3582    ws = in_ws
3583    IF ( ws < .1_wp )  THEN
3584       ws = .1_wp
3585    ENDIF
3586!
3587!-- Heat_convection = forced convection
3588    heat_convection = 12.1_wp * SQRT ( ws * pair / 1013.25_wp )
3589!
3590!-- Average skin temperature
3591    t_skin_aver = 35.7_wp - .0275_wp * activity
3592!
3593!-- Calculation of constants for evaluation below
3594    bc = .155_wp * clo * 3.96_wp * 10._wp**( -8._wp ) * f_cl
3595    cc = f_cl * heat_convection
3596    ec = .155_wp * clo
3597    dc = ( 1._wp + ec * cc ) / bc
3598    gc = ( t_skin_aver + bc * ( tmrt + 273.2_wp )**4._wp + ec * cc * ta ) / bc
3599!
3600!-- Calculation of clothing surface temperature (t_clothing) based on
3601!-- newton-approximation with air temperature as initial guess
3602    niter = INT( dt * 10._wp, KIND=iwp )
3603    IF ( niter < 1 )  niter = 1_iwp
3604    adjustrate = 1._wp - EXP ( -1._wp * ( 10._wp / time_equil ) * dt )
3605    IF ( adjustrate >= 1._wp )  adjustrate = 1._wp
3606    adjustrate_cloth = adjustrate * 30._wp
3607    t_clothing = t_cloth
3608
3609    IF ( t_cloth <= -998.0_wp )  THEN  ! If initial run
3610       niter = 3_iwp
3611       adjustrate = 1._wp
3612       adjustrate_cloth = 1._wp
3613       t_clothing = ta
3614    ENDIF
3615
3616    DO  i = 1, niter
3617       t_clothing = t_clothing - adjustrate_cloth * ( ( t_clothing +           &
3618          273.2_wp )**4._wp + t_clothing *                                     &
3619          dc - gc ) / ( 4._wp * ( t_clothing + 273.2_wp )**3._wp + dc )
3620    ENDDO
3621!
3622!-- Empiric factor for the adaption of the heat ballance equation
3623!-- to the psycho-physical scale (Equ. 40 in FANGER)
3624    z1 = ( .303_wp * EXP ( -.036_wp * actlev ) + .0275_wp )
3625!
3626!-- Water vapour diffution through the skin
3627    z2 = .31_wp * ( 57.3_wp - .07_wp * activity-pa )
3628!
3629!-- Sweat evaporation from the skin surface
3630    z3 = .42_wp * ( activity - 58._wp )
3631!
3632!-- Loss of latent heat through respiration
3633    z4 = .0017_wp * actlev * ( 58.7_wp - pa ) + .0014_wp * actlev *            &
3634      ( 34._wp - ta )
3635!
3636!-- Loss of radiational heat
3637    z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + 273.2_wp )**4 - ( tmrt +         &
3638       273.2_wp )**4 )
3639!
3640!-- Heat loss through forced convection
3641    z6 = cc * ( t_clothing - ta )
3642!
3643!-- Write together as energy ballance
3644    en = activity - z2 - z3 - z4 - z5 - z6
3645!
3646!-- Manage storage
3647    d_s = adjustrate * en + ( 1._wp - adjustrate ) * s
3648!
3649!-- Predicted Mean Vote
3650    pmva = z1 * d_s
3651!
3652!-- Update storage
3653    s = d_s
3654    t_cloth = t_clothing
3655
3656 END SUBROUTINE fanger_s_acti
3657
3658
3659
3660!------------------------------------------------------------------------------!
3661!
3662! Description:
3663! ------------
3664!> Physiologically Equivalent Temperature (PET),
3665!> stationary (calculated based on MEMI),
3666!> Subroutine based on PETBER vers. 1.5.1996 by P. Hoeppe
3667!------------------------------------------------------------------------------!
3668
3669 SUBROUTINE calculate_pet_static( ta, vpa, v, tmrt, pair, pet_ij )
3670
3671    IMPLICIT NONE
3672!
3673!-- Input arguments:
3674    REAL(wp), INTENT( IN ) ::  ta    !< Air temperature             (degree_C)
3675    REAL(wp), INTENT( IN ) ::  tmrt  !< Mean radiant temperature    (degree_C)
3676    REAL(wp), INTENT( IN ) ::  v     !< Wind speed                  (m/s)
3677    REAL(wp), INTENT( IN ) ::  vpa   !< Vapor pressure              (hPa)
3678    REAL(wp), INTENT( IN ) ::  pair  !< Air pressure                (hPa)
3679!
3680!-- Output arguments:
3681    REAL(wp), INTENT ( OUT ) ::  pet_ij  !< PET                     (degree_C)
3682!
3683!-- Internal variables:
3684    REAL(wp) ::  acl        !< clothing area                        (m²)
3685    REAL(wp) ::  adu        !< Du Bois area                         (m²)
3686    REAL(wp) ::  aeff       !< effective area                       (m²)
3687    REAL(wp) ::  ere        !< energy ballance                      (W)
3688    REAL(wp) ::  erel       !< latent energy ballance               (W)
3689    REAL(wp) ::  esw        !< Energy-loss through sweat evap.      (W)
3690    REAL(wp) ::  facl       !< Surface area extension through clothing (factor)
3691    REAL(wp) ::  feff       !< Surface modification by posture      (factor)
3692    REAL(wp) ::  rdcl       !< Diffusion resistence of clothing     (factor)
3693    REAL(wp) ::  rdsk       !< Diffusion resistence of skin         (factor)
3694    REAL(wp) ::  rtv
3695    REAL(wp) ::  vpts       !< Sat. vapor pressure over skin        (hPa)
3696    REAL(wp) ::  tsk        !< Skin temperature                     (degree_C)
3697    REAL(wp) ::  tcl        !< Clothing temperature                 (degree_C)
3698    REAL(wp) ::  wetsk      !< Fraction of wet skin                 (factor)
3699!
3700!-- Variables:
3701    REAL(wp) :: int_heat    !< Internal heat        (W)
3702!
3703!-- MEMI configuration
3704    REAL(wp) :: age         !< Persons age          (a)
3705    REAL(wp) :: mbody       !< Persons body mass    (kg)
3706    REAL(wp) :: ht          !< Persons height       (m)
3707    REAL(wp) :: work        !< Work load            (W)
3708    REAL(wp) :: eta         !< Work efficiency      (dimensionless)
3709    REAL(wp) :: clo         !< Clothing insulation index (clo)
3710    REAL(wp) :: fcl         !< Surface area modification by clothing (factor)
3711!     INTEGER(iwp) :: pos     !< Posture: 1 = standing, 2 = sitting
3712!     INTEGER(iwp) :: sex     !< Sex: 1 = male, 2 = female
3713!
3714!-- Configuration, keep standard parameters!
3715    age   = 35._wp
3716    mbody = 75._wp
3717    ht    =  1.75_wp
3718    work  = 80._wp
3719    eta   =  0._wp
3720    clo   =  0.9_wp
3721    fcl   =  1.15_wp
3722!
3723!-- Call subfunctions
3724    CALL in_body ( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta,    &
3725            vpa, work )
3726
3727    CALL heat_exch ( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff, ht, &
3728            int_heat,mbody, pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa,      &
3729            vpts, wetsk )
3730
3731    CALL pet_iteration ( acl, adu, aeff, esw, facl, feff, int_heat, pair,      &
3732            rdcl, rdsk, rtv, ta, tcl, tsk, pet_ij, vpts, wetsk )
3733
3734
3735 END SUBROUTINE calculate_pet_static
3736
3737
3738!------------------------------------------------------------------------------!
3739! Description:
3740! ------------
3741!> Calculate internal energy ballance
3742!------------------------------------------------------------------------------!
3743 SUBROUTINE in_body ( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta, &
3744    vpa, work )
3745!
3746!-- Input arguments:
3747    REAL(wp), INTENT( IN )  ::  pair      !< air pressure             (hPa)
3748    REAL(wp), INTENT( IN )  ::  ta        !< air temperature          (degree_C)
3749    REAL(wp), INTENT( IN )  ::  vpa       !< vapor pressure           (hPa)
3750    REAL(wp), INTENT( IN )  ::  age       !< Persons age              (a)
3751    REAL(wp), INTENT( IN )  ::  mbody     !< Persons body mass        (kg)
3752    REAL(wp), INTENT( IN )  ::  ht        !< Persons height           (m)
3753    REAL(wp), INTENT( IN )  ::  work      !< Work load                (W)
3754    REAL(wp), INTENT( IN )  ::  eta       !< Work efficiency      (dimensionless)
3755!
3756!-- Output arguments:
3757    REAL(wp), INTENT( OUT ) ::  ere       !< energy ballance          (W)
3758    REAL(wp), INTENT( OUT ) ::  erel      !< latent energy ballance   (W)
3759    REAL(wp), INTENT( OUT ) ::  int_heat  !< internal heat production (W)
3760    REAL(wp), INTENT( OUT ) ::  rtv       !< respiratory volume
3761!
3762!-- Internal variables:
3763    REAL(wp) ::  eres                     !< Sensible respiratory heat flux (W)
3764    REAL(wp) ::  met
3765    REAL(wp) ::  tex
3766    REAL(wp) ::  vpex
3767
3768!
3769!-- Metabolic heat production
3770    met = 3.45_wp * mbody**( 3._wp / 4._wp ) * (1._wp + 0.004_wp *             &
3771          ( 30._wp - age) + 0.010_wp * ( ( ht * 100._wp /                      &
3772          ( mbody**( 1._wp / 3._wp ) ) ) - 43.4_wp ) )
3773    met = work + met
3774    int_heat = met * (1._wp - eta)
3775!
3776!-- Sensible respiration energy
3777    tex  = 0.47_wp * ta + 21.0_wp
3778    rtv  = 1.44_wp * 10._wp**(-6._wp) * met
3779    eres = c_p * (ta - tex) * rtv
3780!
3781!-- Latent respiration energy
3782    vpex = 6.11_wp * 10._wp**( 7.45_wp * tex / ( 235._wp + tex ) )
3783    erel = 0.623_wp * l_v / pair * ( vpa - vpex ) * rtv
3784!
3785!-- Sum of the results
3786    ere = eres + erel
3787
3788 END SUBROUTINE in_body
3789
3790
3791!------------------------------------------------------------------------------!
3792! Description:
3793! ------------
3794!> Calculate heat gain or loss
3795!------------------------------------------------------------------------------!
3796 SUBROUTINE heat_exch ( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff,  &
3797        ht, int_heat, mbody, pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa,     &
3798        vpts, wetsk )
3799
3800!
3801!-- Input arguments:
3802    REAL(wp), INTENT( IN )  ::  ere    !< Energy ballance          (W)
3803    REAL(wp), INTENT( IN )  ::  erel   !< Latent energy ballance   (W)
3804    REAL(wp), INTENT( IN )  ::  int_heat  !< internal heat production (W)
3805    REAL(wp), INTENT( IN )  ::  pair   !< Air pressure             (hPa)
3806    REAL(wp), INTENT( IN )  ::  ta     !< Air temperature          (degree_C)
3807    REAL(wp), INTENT( IN )  ::  tmrt   !< Mean radiant temperature (degree_C)
3808    REAL(wp), INTENT( IN )  ::  v      !< Wind speed               (m/s)
3809    REAL(wp), INTENT( IN )  ::  vpa    !< Vapor pressure           (hPa)
3810    REAL(wp), INTENT( IN )  ::  mbody  !< body mass                (kg)
3811    REAL(wp), INTENT( IN )  ::  ht     !< height                   (m)
3812    REAL(wp), INTENT( IN )  ::  clo    !< clothing insulation      (clo)
3813    REAL(wp), INTENT( IN )  ::  fcl    !< factor for surface area increase by clothing
3814!
3815!-- Output arguments:
3816    REAL(wp), INTENT( OUT ) ::  acl    !< Clothing surface area        (m²)
3817    REAL(wp), INTENT( OUT ) ::  adu    !< Du-Bois area                 (m²)
3818    REAL(wp), INTENT( OUT ) ::  aeff   !< Effective surface area       (m²)
3819    REAL(wp), INTENT( OUT ) ::  esw    !< Energy-loss through sweat evap. (W)
3820    REAL(wp), INTENT( OUT ) ::  facl   !< Surface area extension through clothing (factor)
3821    REAL(wp), INTENT( OUT ) ::  feff   !< Surface modification by posture (factor)
3822    REAL(wp), INTENT( OUT ) ::  rdcl   !< Diffusion resistence of clothing (factor)
3823    REAL(wp), INTENT( OUT ) ::  rdsk   !< Diffusion resistence of skin (factor)
3824    REAL(wp), INTENT( OUT ) ::  tcl    !< Clothing temperature         (degree_C)
3825    REAL(wp), INTENT( OUT ) ::  tsk    !< Skin temperature             (degree_C)
3826    REAL(wp), INTENT( OUT ) ::  vpts   !< Sat. vapor pressure over skin (hPa)
3827    REAL(wp), INTENT( OUT ) ::  wetsk  !< Fraction of wet skin (dimensionless)
3828!
3829!-- Cconstants:
3830!     REAL(wp), PARAMETER :: cair = 1010._wp      !< replaced by c_p
3831    REAL(wp), PARAMETER :: cb   = 3640._wp        !<
3832    REAL(wp), PARAMETER :: emcl =    0.95_wp      !< Longwave emission coef. of cloth
3833    REAL(wp), PARAMETER :: emsk =    0.99_wp      !< Longwave emission coef. of skin
3834!     REAL(wp), PARAMETER :: evap = 2.42_wp * 10._wp **6._wp  !< replaced by l_v
3835    REAL(wp), PARAMETER :: food =    0._wp        !< Heat gain by food        (W)
3836    REAL(wp), PARAMETER :: po   = 1013.25_wp      !< Air pressure at sea level (hPa)
3837    REAL(wp), PARAMETER :: rob  =    1.06_wp      !<
3838!
3839!-- Internal variables
3840    REAL(wp) ::  c(0:10)        !< Core temperature array           (degree_C)
3841    REAL(wp) ::  cbare          !< Convection through bare skin
3842    REAL(wp) ::  cclo           !< Convection through clothing
3843    REAL(wp) ::  csum           !< Convection in total
3844    REAL(wp) ::  di             !< difference between r1 and r2
3845    REAL(wp) ::  ed             !< energy transfer by diffusion     (W)
3846    REAL(wp) ::  enbal          !< energy ballance                  (W)
3847    REAL(wp) ::  enbal2         !< energy ballance (storage, last cycle)
3848    REAL(wp) ::  eswdif         !< difference between sweat production and evaporation potential
3849    REAL(wp) ::  eswphy         !< sweat created by physiology
3850    REAL(wp) ::  eswpot         !< potential sweat evaporation
3851    REAL(wp) ::  fec            !<
3852    REAL(wp) ::  hc             !<
3853    REAL(wp) ::  he             !<
3854    REAL(wp) ::  htcl           !<
3855    REAL(wp) ::  r1             !<
3856    REAL(wp) ::  r2             !<
3857    REAL(wp) ::  rbare          !< Radiational loss of bare skin    (W/m²)
3858    REAL(wp) ::  rcl            !<
3859    REAL(wp) ::  rclo           !< Radiational loss of clothing     (W/m²)
3860    REAL(wp) ::  rclo2          !< Longwave radiation gain or loss  (W/m²)
3861    REAL(wp) ::  rsum           !< Radiational loss or gain         (W/m²)
3862    REAL(wp) ::  sw             !<
3863!     REAL(wp) ::  swf            !< female factor, currently unused
3864    REAL(wp) ::  swm            !<
3865    REAL(wp) ::  tbody          !<
3866    REAL(wp) ::  tcore(1:7)     !<
3867    REAL(wp) ::  vb             !<
3868    REAL(wp) ::  vb1            !<
3869    REAL(wp) ::  vb2            !<
3870    REAL(wp) ::  wd             !<
3871    REAL(wp) ::  wr             !<
3872    REAL(wp) ::  ws             !<
3873    REAL(wp) ::  wsum           !<
3874    REAL(wp) ::  xx             !< modification step                (K)
3875    REAL(wp) ::  y              !< fraction of bare skin
3876    INTEGER(iwp) ::  count1     !< running index
3877    INTEGER(iwp) ::  count3     !< running index
3878    INTEGER(iwp) ::  j          !< running index
3879    INTEGER(iwp) ::  i          !< running index
3880    LOGICAL ::  skipIncreaseCount   !< iteration control flag
3881
3882!
3883!-- Initialize
3884    wetsk = 0._wp  !< skin is dry everywhere on init (no non-evaporated sweat)
3885!
3886!-- Set Du Bois Area for the sample person
3887    adu = 0.203_wp * mbody**0.425_wp * ht**0.725_wp
3888!
3889!-- Initialize convective heat considering local air preassure
3890    hc = 2.67_wp + ( 6.5_wp * v**0.67_wp )
3891    hc = hc * ( pair / po )**0.55_wp
3892!
3893!-- Set surface modification by posture (the person will always stand)
3894    feff = 0.725_wp                     !< Posture: 0.725 for stading
3895!
3896!-- Set surface modification by clothing
3897    facl = ( - 2.36_wp + 173.51_wp * clo - 100.76_wp * clo * clo + 19.28_wp    &
3898          * ( clo**3._wp ) ) / 100._wp
3899    IF ( facl > 1._wp )  facl = 1._wp
3900!
3901!-- Initialize heat resistences
3902    rcl = ( clo / 6.45_wp ) / facl
3903    IF ( clo >= 2._wp )  y  = 1._wp
3904    IF ( ( clo > 0.6_wp )   .AND.  ( clo < 2._wp ) )   y = ( ht - 0.2_wp ) / ht
3905    IF ( ( clo <= 0.6_wp )  .AND.  ( clo > 0.3_wp ) )  y = 0.5_wp
3906    IF ( ( clo <= 0.3_wp )  .AND.  ( clo > 0._wp ) )   y = 0.1_wp
3907    r2   = adu * ( fcl - 1._wp + facl ) / ( 2._wp * 3.14_wp * ht * y )
3908    r1   = facl * adu / ( 2._wp * 3.14_wp * ht * y )
3909    di   = r2 - r1
3910
3911!
3912!-- Estimate skin temperatur
3913    DO  j = 1, 7
3914
3915       tsk    = 34._wp
3916       count1 = 0_iwp
3917       tcl    = ( ta + tmrt + tsk ) / 3._wp
3918       count3 = 1_iwp
3919       enbal2 = 0._wp
3920
3921       DO  i = 1, 100  ! allow for 100 iterations max
3922          acl   = adu * facl + adu * ( fcl - 1._wp )
3923          rclo2 = emcl * sigma_sb * ( ( tcl + degc_to_k )**4._wp -             &
3924            ( tmrt + degc_to_k )**4._wp ) * feff
3925          htcl  = 6.28_wp * ht * y * di / ( rcl * LOG( r2 / r1 ) * acl )
3926          tsk   = 1._wp / htcl * ( hc * ( tcl - ta ) + rclo2 ) + tcl
3927!
3928!--       Radiation saldo
3929          aeff  = adu * feff
3930          rbare = aeff * ( 1._wp - facl ) * emsk * sigma_sb *                  &
3931            ( ( tmrt + degc_to_k )**4._wp - ( tsk + degc_to_k )**4._wp )
3932          rclo  = feff * acl * emcl * sigma_sb *                               &
3933            ( ( tmrt + degc_to_k )**4._wp - ( tcl + degc_to_k )**4._wp )
3934          rsum  = rbare + rclo
3935!
3936!--       Convection
3937          cbare = hc * ( ta - tsk ) * adu * ( 1._wp - facl )
3938          cclo  = hc * ( ta - tcl ) * acl
3939          csum  = cbare + cclo
3940!
3941!--       Core temperature
3942          c(0)  = int_heat + ere
3943          c(1)  = adu * rob * cb
3944          c(2)  = 18._wp - 0.5_wp * tsk
3945          c(3)  = 5.28_wp * adu * c(2)
3946          c(4)  = 0.0208_wp * c(1)
3947          c(5)  = 0.76075_wp * c(1)
3948          c(6)  = c(3) - c(5) - tsk * c(4)
3949          c(7)  = - c(0) * c(2) - tsk * c(3) + tsk * c(5)
3950          c(8)  = c(6) * c(6) - 4._wp * c(4) * c(7)
3951          c(9)  = 5.28_wp * adu - c(5) - c(4) * tsk
3952          c(10) = c(9) * c(9) - 4._wp * c(4) *                                 &
3953                  ( c(5) * tsk - c(0) - 5.28_wp * adu * tsk )
3954
3955          IF ( tsk == 36._wp )  tsk = 36.01_wp
3956          tcore(7) = c(0) / ( 5.28_wp * adu + c(1) * 6.3_wp / 3600._wp ) + tsk
3957          tcore(3) = c(0) / ( 5.28_wp * adu + ( c(1) * 6.3_wp / 3600._wp ) /   &
3958            ( 1._wp + 0.5_wp * ( 34._wp - tsk ) ) ) + tsk
3959          IF ( c(10) >= 0._wp )  THEN
3960             tcore(6) = ( - c(9) - c(10)**0.5_wp ) / ( 2._wp * c(4) )
3961             tcore(1) = ( - c(9) + c(10)**0.5_wp ) / ( 2._wp * c(4) )
3962          ENDIF
3963
3964          IF ( c(8) >= 0._wp )  THEN
3965             tcore(2) = ( - c(6) + ABS( c(8) )**0.5_wp ) / ( 2._wp * c(4) )
3966             tcore(5) = ( - c(6) - ABS( c(8) )**0.5_wp ) / ( 2._wp * c(4) )
3967             tcore(4) = c(0) / ( 5.28_wp * adu + c(1) * 1._wp / 40._wp ) + tsk
3968          ENDIF
3969!
3970!--       Transpiration
3971          tbody = 0.1_wp * tsk + 0.9_wp * tcore(j)
3972          swm   = 304.94_wp * ( tbody - 36.6_wp ) * adu / 3600000._wp
3973          vpts  = 6.11_wp * 10._wp**( 7.45_wp * tsk / ( 235._wp + tsk ) )
3974
3975          IF ( tbody <= 36.6_wp )  swm = 0._wp  !< no need for sweating
3976
3977          sw = swm
3978          eswphy = - sw * l_v
3979          he     = 0.633_wp * hc / ( pair * c_p )
3980          fec    = 1._wp / ( 1._wp + 0.92_wp * hc * rcl )
3981          eswpot = he * ( vpa - vpts ) * adu * l_v * fec
3982          wetsk  = eswphy / eswpot
3983
3984          IF ( wetsk > 1._wp )  wetsk = 1._wp
3985!
3986!--       Sweat production > evaporation?
3987          eswdif = eswphy - eswpot
3988
3989          IF ( eswdif <= 0._wp )  esw = eswpot  !< Limit is evaporation
3990          IF ( eswdif > 0._wp )   esw = eswphy  !< Limit is sweat production
3991          IF ( esw  > 0._wp )     esw = 0._wp   !< Sweat can't be evaporated, no more cooling effect
3992!
3993!--       Diffusion
3994          rdsk = 0.79_wp * 10._wp**7._wp
3995          rdcl = 0._wp
3996          ed   = l_v / ( rdsk + rdcl ) * adu * ( 1._wp - wetsk ) * ( vpa -     &
3997             vpts )
3998!
3999!--       Max vb
4000          vb1 = 34._wp - tsk
4001          vb2 = tcore(j) - 36.6_wp
4002
4003          IF ( vb2 < 0._wp )  vb2 = 0._wp
4004          IF ( vb1 < 0._wp )  vb1 = 0._wp
4005          vb = ( 6.3_wp + 75._wp * vb2 ) / ( 1._wp + 0.5_wp * vb1 )
4006!
4007!--       Energy ballence
4008          enbal = int_heat + ed + ere + esw + csum + rsum + food
4009!
4010!--       Clothing temperature
4011          xx = 0.001_wp
4012          IF ( count1 == 0_iwp )  xx = 1._wp
4013          IF ( count1 == 1_iwp )  xx = 0.1_wp
4014          IF ( count1 == 2_iwp )  xx = 0.01_wp
4015          IF ( count1 == 3_iwp )  xx = 0.001_wp
4016
4017          IF ( enbal > 0._wp )  tcl = tcl + xx
4018          IF ( enbal < 0._wp )  tcl = tcl - xx
4019
4020          skipIncreaseCount = .FALSE.
4021          IF ( ( (enbal <= 0._wp )  .AND.  (enbal2 > 0._wp ) )  .OR.           &
4022             ( ( enbal >= 0._wp )  .AND.  ( enbal2 < 0._wp ) ) )  THEN
4023             skipIncreaseCount = .TRUE.
4024          ELSE
4025             enbal2 = enbal
4026             count3 = count3 + 1_iwp
4027          ENDIF
4028
4029          IF ( ( count3 > 200_iwp )  .OR.  skipIncreaseCount )  THEN
4030             IF ( count1 < 3_iwp )  THEN
4031                count1 = count1 + 1_iwp
4032                enbal2 = 0._wp
4033             ELSE
4034                EXIT
4035             ENDIF
4036          ENDIF
4037       ENDDO
4038
4039       IF ( count1 == 3_iwp )  THEN
4040          SELECT CASE ( j )
4041             CASE ( 2, 5) 
4042                IF ( .NOT. ( ( tcore(j) >= 36.6_wp )  .AND.                    &
4043                   ( tsk <= 34.050_wp ) ) )  CYCLE
4044             CASE ( 6, 1 )
4045                IF ( c(10) < 0._wp ) CYCLE
4046                IF ( .NOT. ( ( tcore(j) >= 36.6_wp )  .AND.                    &
4047                   ( tsk > 33.850_wp ) ) )  CYCLE
4048             CASE ( 3 )
4049                IF ( .NOT. ( ( tcore(j) < 36.6_wp )  .AND.                     &
4050                   ( tsk <= 34.000_wp ) ) )  CYCLE
4051             CASE ( 7 )
4052                IF ( .NOT. ( ( tcore(j) < 36.6_wp )  .AND.                     &
4053                   ( tsk > 34.000_wp ) ) )  CYCLE
4054             CASE default
4055          END SELECT
4056       ENDIF
4057
4058       IF ( ( j /= 4_iwp )  .AND.  ( vb >= 91._wp ) )  CYCLE
4059       IF ( ( j == 4_iwp )  .AND.  ( vb < 89._wp ) )  CYCLE
4060       IF ( vb > 90._wp ) vb = 90._wp
4061!
4062!--    Loses by water
4063       ws = sw * 3600._wp * 1000._wp
4064       IF ( ws > 2000._wp )  ws = 2000._wp
4065       wd = ed / l_v * 3600._wp * ( -1000._wp )
4066       wr = erel / l_v * 3600._wp * ( -1000._wp )
4067
4068       wsum = ws + wr + wd
4069
4070       RETURN
4071    ENDDO
4072 END SUBROUTINE heat_exch
4073
4074!------------------------------------------------------------------------------!
4075! Description:
4076! ------------
4077!> Calculate PET
4078!------------------------------------------------------------------------------!
4079 SUBROUTINE pet_iteration ( acl, adu, aeff, esw, facl, feff, int_heat, pair,   &
4080        rdcl, rdsk, rtv, ta, tcl, tsk, pet_ij, vpts, wetsk )
4081!
4082!-- Input arguments:
4083    REAL(wp), INTENT( IN ) ::  acl   !< clothing surface area        (m²)
4084    REAL(wp), INTENT( IN ) ::  adu   !< Du-Bois area                 (m²)
4085    REAL(wp), INTENT( IN ) ::  esw   !< energy-loss through sweat evap. (W)
4086    REAL(wp), INTENT( IN ) ::  facl  !< surface area extension through clothing (factor)
4087    REAL(wp), INTENT( IN ) ::  feff  !< surface modification by posture (factor)
4088    REAL(wp), INTENT( IN ) ::  int_heat  !< internal heat production (W)
4089    REAL(wp), INTENT( IN ) ::  pair  !< air pressure                 (hPa)
4090    REAL(wp), INTENT( IN ) ::  rdcl  !< diffusion resistence of clothing (factor)
4091    REAL(wp), INTENT( IN ) ::  rdsk  !< diffusion resistence of skin (factor)
4092    REAL(wp), INTENT( IN ) ::  rtv   !< respiratory volume
4093    REAL(wp), INTENT( IN ) ::  ta    !< air temperature              (degree_C)
4094    REAL(wp), INTENT( IN ) ::  tcl   !< clothing temperature         (degree_C)
4095    REAL(wp), INTENT( IN ) ::  tsk   !< skin temperature             (degree_C)
4096    REAL(wp), INTENT( IN ) ::  vpts  !< sat. vapor pressure over skin (hPa)
4097    REAL(wp), INTENT( IN ) ::  wetsk !< fraction of wet skin (dimensionless)
4098!
4099!-- Output arguments:
4100    REAL(wp), INTENT( OUT ) ::  aeff     !< effective surface area       (m²)
4101    REAL(wp), INTENT( OUT ) ::  pet_ij   !< PET                          (degree_C)
4102!
4103!-- Cconstants:
4104    REAL(wp), PARAMETER :: emcl =    0.95_wp      !< Longwave emission coef. of cloth
4105    REAL(wp), PARAMETER :: emsk =    0.99_wp      !< Longwave emission coef. of skin
4106    REAL(wp), PARAMETER :: po   = 1013.25_wp      !< Air pressure at sea level (hPa)
4107!
4108!-- Internal variables
4109    REAL ( wp ) ::  cbare             !< Convection through bare skin
4110    REAL ( wp ) ::  cclo              !< Convection through clothing
4111    REAL ( wp ) ::  csum              !< Convection in total
4112    REAL ( wp ) ::  ed                !< Diffusion                      (W)
4113    REAL ( wp ) ::  enbal             !< Energy ballance                (W)
4114    REAL ( wp ) ::  enbal2            !< Energy ballance (last iteration cycle)
4115    REAL ( wp ) ::  ere               !< Energy ballance result         (W)
4116    REAL ( wp ) ::  erel              !< Latent energy ballance         (W)
4117    REAL ( wp ) ::  eres              !< Sensible respiratory heat flux (W)
4118    REAL ( wp ) ::  hc                !<
4119    REAL ( wp ) ::  rbare             !< Radiational loss of bare skin  (W/m²)
4120    REAL ( wp ) ::  rclo              !< Radiational loss of clothing   (W/m²)
4121    REAL ( wp ) ::  rsum              !< Radiational loss or gain       (W/m²)
4122    REAL ( wp ) ::  tex               !< Temperat. of exhaled air       (degree_C)
4123    REAL ( wp ) ::  vpex              !< Vapor pressure of exhaled air  (hPa)
4124    REAL ( wp ) ::  xx                !< Delta PET per iteration        (K)
4125
4126    INTEGER ( iwp ) ::  count1        !< running index
4127    INTEGER ( iwp ) ::  i             !< running index
4128
4129    pet_ij = ta
4130    enbal2 = 0._wp
4131
4132    DO  count1 = 0, 3
4133       DO  i = 1, 125  ! 500 / 4
4134          hc = 2.67_wp + 6.5_wp * 0.1_wp**0.67_wp
4135          hc = hc * ( pair / po )**0.55_wp
4136!
4137!--       Radiation
4138          aeff  = adu * feff
4139          rbare = aeff * ( 1._wp - facl ) * emsk * sigma_sb *                  &
4140              ( ( pet_ij + degc_to_k )**4._wp - ( tsk + degc_to_k )**4._wp )
4141          rclo  = feff * acl * emcl * sigma_sb *                               &
4142              ( ( pet_ij + degc_to_k )**4._wp - ( tcl + degc_to_k )**4._wp )
4143          rsum  = rbare + rclo
4144!
4145!--       Covection
4146          cbare = hc * ( pet_ij - tsk ) * adu * ( 1._wp - facl )
4147          cclo  = hc * ( pet_ij - tcl ) * acl
4148          csum  = cbare + cclo
4149!
4150!--       Diffusion
4151          ed = l_v / ( rdsk + rdcl ) * adu * ( 1._wp - wetsk ) * ( 12._wp -    &
4152             vpts )
4153!
4154!--       Respiration
4155          tex  = 0.47_wp * pet_ij + 21._wp
4156          eres = c_p * ( pet_ij - tex ) * rtv
4157          vpex = 6.11_wp * 10._wp**( 7.45_wp * tex / ( 235._wp + tex ) )
4158          erel = 0.623_wp * l_v / pair * ( 12._wp - vpex ) * rtv
4159          ere  = eres + erel
4160!
4161!--       Energy ballance
4162          enbal = int_heat + ed + ere + esw + csum + rsum
4163!
4164!--       Iteration concerning ta
4165          xx = 0.001_wp
4166          IF ( count1 == 0_iwp )  xx = 1._wp
4167          IF ( count1 == 1_iwp )  xx = 0.1_wp
4168          IF ( count1 == 2_iwp )  xx = 0.01_wp
4169!           IF ( count1 == 3_iwp )  xx = 0.001_wp
4170          IF ( enbal > 0._wp )  pet_ij = pet_ij - xx
4171          IF ( enbal < 0._wp )  pet_ij = pet_ij + xx
4172          IF ( ( enbal <= 0._wp )  .AND.  ( enbal2 > 0._wp ) )  EXIT
4173          IF ( ( enbal >= 0._wp )  .AND.  ( enbal2 < 0._wp ) )  EXIT
4174
4175          enbal2 = enbal
4176       ENDDO
4177    ENDDO
4178 END SUBROUTINE pet_iteration
4179
4180!
4181!-- UVEM specific subroutines
4182
4183!---------------------------------------------------------------------------------------------------------------------!
4184! Description:
4185! ------------
4186!> Module-specific routine for new module
4187!---------------------------------------------------------------------------------------------------------------------!
4188 SUBROUTINE uvem_solar_position
4189   
4190    USE date_and_time_mod,                                                                                            &
4191       ONLY:  calc_date_and_time, day_of_year, time_utc
4192   
4193    USE control_parameters,                                                                                           &
4194       ONLY:  latitude, longitude   
4195
4196    IMPLICIT NONE
4197   
4198   
4199    REAL(wp) ::  alpha       = 0.0_wp   !< solar azimuth angle in radiant   
4200    REAL(wp) ::  doy_r       = 0.0_wp   !< real format of day_of_year           
4201    REAL(wp) ::  declination = 0.0_wp   !< declination
4202    REAL(wp) ::  dtor        = 0.0_wp   !< factor to convert degree to radiant
4203    REAL(wp) ::  js          = 0.0_wp   !< parameter for solar position calculation
4204    REAL(wp) ::  lat         = 52.39_wp !< latitude
4205    REAL(wp) ::  lon         = 9.7_wp   !< longitude       
4206    REAL(wp) ::  thetar      = 0.0_wp   !< angle for solar zenith angle calculation
4207    REAL(wp) ::  thetasr     = 0.0_wp   !< angle for solar azimuth angle calculation   
4208    REAL(wp) ::  zgl         = 0.0_wp   !< calculated exposure by direct beam   
4209    REAL(wp) ::  woz         = 0.0_wp   !< calculated exposure by diffuse radiation
4210    REAL(wp) ::  wsp         = 0.0_wp   !< calculated exposure by direct beam   
4211   
4212
4213    CALL calc_date_and_time
4214    doy_r = real(day_of_year)   
4215    dtor = pi / 180.0_wp
4216    lat = latitude
4217    lon = longitude
4218!
4219!-- calculation of js, necessary for calculation of equation of time (zgl) :
4220    js=  72.0_wp * ( doy_r + ( time_utc / 86400.0_wp ) ) / 73.0_wp 
4221!
4222!-- calculation of equation of time (zgl):
4223    zgl = 0.0066_wp + 7.3525_wp * cos( ( js + 85.9_wp ) * dtor ) + 9.9359_wp *                                        &
4224    cos( ( 2.0_wp * js + 108.9_wp ) * dtor ) + 0.3387_wp * cos( ( 3 * js + 105.2_wp ) * dtor )
4225!
4226!-- calculation of apparent solar time woz:
4227    woz = ( ( time_utc / 3600.0_wp ) - ( 4.0_wp * ( 15.0_wp - lon ) ) / 60.0_wp ) + ( zgl / 60.0_wp )
4228!
4229!-- calculation of hour angle (wsp):
4230    wsp = ( woz - 12.0_wp ) * 15.0_wp
4231!
4232!-- calculation of declination:
4233    declination = 0.3948_wp - 23.2559_wp * cos( ( js + 9.1_wp ) * dtor ) -                                            &
4234    0.3915_wp * cos( ( 2.0_wp * js + 5.4_wp ) * dtor ) - 0.1764_wp * cos( ( 3.0_wp * js + 26.0_wp ) * dtor )
4235!
4236!-- calculation of solar zenith angle
4237    thetar  = acos( sin( lat * dtor) * sin( declination * dtor ) + cos( wsp * dtor ) *                                &
4238    cos( lat * dtor ) * cos( declination * dtor ) )
4239    thetasr = asin( sin( lat * dtor) * sin( declination * dtor ) + cos( wsp * dtor ) *                                & 
4240    cos( lat * dtor ) * cos( declination * dtor ) )
4241    sza = thetar / dtor
4242!
4243!-- calculation of solar azimuth angle
4244    IF (woz .LE. 12.0_wp) alpha = pi - acos( ( sin(thetasr) * sin( lat * dtor ) -                                     &
4245    sin( declination * dtor ) ) / ( cos(thetasr) * cos( lat * dtor ) ) )   
4246    IF (woz .GT. 12.0_wp) alpha = pi + acos( ( sin(thetasr) * sin( lat * dtor ) -                                     &
4247    sin( declination * dtor ) ) / ( cos(thetasr) * cos( lat * dtor ) ) )   
4248    saa = alpha / dtor
4249
4250 END SUBROUTINE uvem_solar_position
4251
4252
4253!------------------------------------------------------------------------------!
4254! Description:
4255! ------------
4256!> Module-specific routine for new module
4257!---------------------------------------------------------------------------------------------------------------------!
4258 SUBROUTINE uvem_calc_exposure
4259
4260    USE indices,                                                                                                      &
4261        ONLY:  nys, nyn, nxl, nxr
4262   
4263   
4264    IMPLICIT NONE   
4265   
4266    INTEGER(iwp) ::  i     !< loop index in x direction
4267    INTEGER(iwp) ::  j     !< loop index in y direction
4268    INTEGER(iwp) ::  szai  !< loop index for different sza values
4269
4270    CALL uvem_solar_position
4271     
4272    IF (sza  .GE.  90)  THEN
4273       vitd3_exposure(:,:) = 0.0_wp
4274    ELSE
4275       
4276       DO  ai = 0, 35
4277          DO  zi = 0, 9
4278                projection_area_lookup_table(ai,zi) = uvem_projarea_f%var(clothing,zi,ai)
4279          ENDDO
4280       ENDDO
4281       DO  ai = 0, 35
4282          DO  zi = 0, 9
4283                integration_array(ai,zi) = uvem_integration_f%var(zi,ai)
4284          ENDDO
4285       ENDDO
4286       DO  ai = 0, 2
4287          DO  zi = 0, 90
4288                irradiance_lookup_table(ai,zi) = uvem_irradiance_f%var(zi,ai)
4289          ENDDO
4290       ENDDO
4291       DO  ai = 0, 35
4292          DO  zi = 0, 9
4293             DO  szai = 0, 90
4294                radiance_lookup_table(ai,zi,szai) = uvem_radiance_f%var(szai,zi,ai)
4295             ENDDO
4296          ENDDO
4297       ENDDO
4298       
4299       
4300       
4301!--    rotate 3D-Model human to desired direction  -----------------------------
4302       projection_area_temp( 0:35,:) = projection_area_lookup_table
4303       projection_area_temp(36:71,:) = projection_area_lookup_table               
4304       IF (  .NOT.  turn_to_sun ) startpos_human = orientation_angle / 10.0_wp
4305       IF (       turn_to_sun ) startpos_human = saa / 10.0_wp       
4306       DO  ai = 0, 35
4307          xfactor = ( startpos_human ) - INT( startpos_human )
4308          DO  zi = 0, 9
4309             projection_area(ai,zi) = ( projection_area_temp( 36 - INT( startpos_human ) - 1 + ai , zi) *             &
4310                                      ( xfactor ) )                                                                   &
4311                                      +( projection_area_temp( 36 - INT( startpos_human ) + ai , zi) *                &
4312                                      ( 1.0_wp - xfactor ) )
4313          ENDDO
4314       ENDDO           
4315!             
4316!           
4317!--    interpolate to accurate Solar Zenith Angle  ------------------         
4318       DO  ai = 0, 35
4319          xfactor = (sza)-INT(sza)
4320          DO  zi = 0, 9
4321             radiance_array(ai,zi) = ( radiance_lookup_table(ai, zi, INT(sza) ) * ( 1.0_wp - xfactor) ) +             &
4322             ( radiance_lookup_table(ai,zi,INT(sza) + 1) * xfactor )
4323          ENDDO
4324       ENDDO
4325       DO  iq = 0, 2
4326          irradiance(iq) = ( irradiance_lookup_table(iq, INT(sza) ) * ( 1.0_wp - xfactor)) +                          &
4327          (irradiance_lookup_table(iq, INT(sza) + 1) * xfactor )
4328       ENDDO   
4329!         
4330!--    interpolate to accurate Solar Azimuth Angle ------------------
4331       IF ( sun_in_south )  THEN
4332          startpos_saa_float = 180.0_wp / 10.0_wp
4333       ELSE
4334          startpos_saa_float = saa / 10.0_wp
4335       ENDIF
4336       radiance_array_temp( 0:35,:) = radiance_array
4337       radiance_array_temp(36:71,:) = radiance_array
4338       xfactor = (startpos_saa_float) - INT(startpos_saa_float)
4339       DO  ai = 0, 35
4340          DO  zi = 0, 9
4341             radiance_array(ai,zi) = ( radiance_array_temp( 36 - INT( startpos_saa_float ) - 1 + ai , zi ) *          &
4342                                     ( xfactor ) )                                                                    &
4343                                     + ( radiance_array_temp( 36 - INT( startpos_saa_float ) + ai , zi )              &
4344                                     * ( 1.0_wp - xfactor ) )
4345          ENDDO
4346       ENDDO
4347!       
4348!     
4349!--    calculate Projectionarea for direct beam -----------------------------'
4350       projection_area_direct_temp( 0:35,:) = projection_area
4351       projection_area_direct_temp(36:71,:) = projection_area
4352       yfactor = ( sza / 10.0_wp ) - INT( sza / 10.0_wp )
4353       xfactor = ( startpos_saa_float ) - INT( startpos_saa_float )
4354       projection_area_direct_beam = ( projection_area_direct_temp( INT(startpos_saa_float)    ,INT(sza/10.0_wp)  ) * &
4355                                     ( 1.0_wp - xfactor ) * ( 1.0_wp - yfactor ) ) +                                  &
4356                                     ( projection_area_direct_temp( INT(startpos_saa_float) + 1,INT(sza/10.0_wp)  ) * &
4357                                     (          xfactor ) * ( 1.0_wp - yfactor ) ) +                                  &
4358                                     ( projection_area_direct_temp( INT(startpos_saa_float)    ,INT(sza/10.0_wp)+1) * &
4359                                     ( 1.0_wp - xfactor ) * (          yfactor ) ) +                                  &
4360                                     ( projection_area_direct_temp( INT(startpos_saa_float) + 1,INT(sza/10.0_wp)+1) * &
4361                                     (          xfactor ) * (          yfactor ) )
4362!                                               
4363!                                               
4364!                                               
4365       DO  i = nxl, nxr
4366          DO  j = nys, nyn
4367!                   
4368! !--        extract obstruction from IBSET-Integer_Array ------------------'
4369             IF (consider_obstructions )  THEN
4370                obstruction_temp1 = building_obstruction_f%var_3d(:,j,i)
4371                IF (obstruction_temp1(0)  .NE.  9)  THEN
4372                   DO  pobi = 0, 44 
4373                      DO  bi = 0, 7 
4374                         IF ( btest( obstruction_temp1(pobi), bi )  .EQV.  .TRUE.)  THEN
4375                            obstruction_temp2( ( pobi * 8 ) + bi ) = 1
4376                         ELSE
4377                            obstruction_temp2( ( pobi * 8 ) + bi ) = 0
4378                         ENDIF
4379                      ENDDO
4380                   ENDDO       
4381                   DO  zi = 0, 9                                         
4382                      obstruction(:,zi) = obstruction_temp2( zi * 36 :( zi * 36) + 35 )
4383                   ENDDO
4384                ELSE
4385                   obstruction(:,:) = 0
4386                ENDIF
4387             ENDIF
4388!             
4389! !--        calculated human exposure ------------------' 
4390             diffuse_exposure = SUM( radiance_array * projection_area * integration_array * obstruction )     
4391         
4392             obstruction_direct_beam = obstruction( nint(startpos_saa_float), nint( sza / 10.0_wp ) ) 
4393             IF (sza  .GE.  89.99_wp)  THEN
4394                sza = 89.99999_wp
4395             ENDIF
4396!             
4397!--          calculate direct normal irradiance (direct beam) ------------------'
4398             direct_exposure = ( irradiance(1) / cos( pi * sza / 180.0_wp ) ) * &
4399             projection_area_direct_beam * obstruction_direct_beam 
4400               
4401             vitd3_exposure(j,i) = ( diffuse_exposure + direct_exposure ) / 1000.0_wp * 70.97_wp 
4402             ! unit = international units vitamin D per second             
4403          ENDDO
4404       ENDDO
4405    ENDIF
4406
4407 END SUBROUTINE uvem_calc_exposure
4408
4409 END MODULE biometeorology_mod
Note: See TracBrowser for help on using the repository browser.