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

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

files re-formatted to follow the PALM coding standard

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