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

Last change on this file since 4442 was 4392, checked in by pavelkrc, 4 years ago

Merge branch resler (updated documentation, optional flux tracing, bugfix)

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