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

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

biometeorology_mod.f90:
(B) Added addittional safety meassures to bio_calculate_thermal_index_maps.
(B) Replaced several REAL (un-)equality comparisons.

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