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

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

biometeorology_mod.f90:
(B) Added safety-meassure to catch the case that 'bio_mrt_av' is stated after 'bio_<index>' in the output section of the p3d file. This could lead to double-averaging and thus to silly results otherwise.

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