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

Last change on this file since 4535 was 4535, checked in by raasch, 2 years ago

bugfix for restart data format query

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