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

Last change on this file since 3716 was 3711, checked in by knoop, 5 years ago

Introduced module_interface_init_checks for post-init checks

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