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

Last change on this file since 3650 was 3650, checked in by kanani, 2 years ago

Bugfix/additions to enable restarts with biometeorology (biometeorology_mod, module_interface)

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