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

Last change on this file since 4495 was 4495, checked in by raasch, 16 months ago

restart data handling with MPI-IO added, first part

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