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

Last change on this file since 4127 was 4127, checked in by suehring, 23 months ago

Merge with branch resler: biomet- output of bio_mrt added; plant_canopy - separate vertical dimension for 3D output (to save disk space); radiation - remove unused plant canopy variables; urban-surface model - do not add anthropogenic heat during wall spin-up

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