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

Last change on this file since 4018 was 3885, checked in by kanani, 5 years ago

restructure/add location/debug messages

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