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

Last change on this file since 3593 was 3593, checked in by kanani, 3 years ago

Bugfix for missing array allocation (biometeorology_mod), remove degree symbol (biometeorology_mod, indoor_model_mod, multi_agent_system_mod, surface_mod, wind_turbine_model_mod)

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