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

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

biometeorology_mod.f90:
(N) Fixed auto-setting of thermal index calculation flags by output as
originally proposed by resler.
(C) removed bio_pet and outher configuration variables.
(C) Updated namelist.
(B) Forcing initialization of tmrt_av_grid to avoid mysterious mrt
values at i==0, j==0

module_interface_mod.f90:
(C) Receiving parameter j (averaging 0==.F./1==.T.) in
module_interface_check_data_output from check_parameters.f90.
(C) Passing j to bio_check_parameters.

check_parameters.f90:
(C) Passing j to module_interface_check_data_output

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