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

Last change on this file since 3614 was 3614, checked in by raasch, 3 years ago

unused variables removed, abort renamed inifor_abort to avoid intrinsic problem in Fortran

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