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

Last change on this file since 4517 was 4517, checked in by raasch, 14 months ago

added restart with MPI-IO for reading local arrays

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