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

Last change on this file since 4144 was 4144, checked in by raasch, 23 months ago

relational operators .EQ., .NE., etc. replaced by ==, /=, etc.

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