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

Last change on this file since 3693 was 3693, checked in by dom_dwd_user, 5 years ago

biometeorology_mod.f90:
(C) renamed averageing switches from e.g. 'aver_q' to 'do_average_q' for better readability
(N) introduced a tmrt_av_grid to store time-averaged mean radiant temperature in a discrete 2d grid (analoge to tmrt_grid). Added mrt_av_grid to restart methods.
(B) replaced mis-named do_average_perct by do_average_theta, as confusing otherwise.
(C) improved general code commenting
(C) 'thermal_comfort' now defaults to .FALSE. (analogue to 'uv_exposure').

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