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

Last change on this file since 3739 was 3739, checked in by dom_dwd_user, 2 years ago

biometeorology_mod.f90:
(N) Auto-adjusting thermal_comfort flag if not set by user, but thermal_indices set as output quantities.
(C) Renamed flags "bio_<index>" to "do_calculate_<index>" for better readability
(C) Removed everything related to "time_bio_results" as this is never used.
(C) Moved humidity warning to check_data_output so it will also be triggered if thermal_comofrt flag was auto-set later.
(B) Fixed bug in mrt calculation introduced with my commit yesterday.

time_integration.f90:
(C) Removed everything related to "time_bio_results" as this is never used

module_interface.f90:
(C) Removed call to empty method bio_check_parameters.

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