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

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

Bugfix: replace simulated_time by time_since_reference_point where required

  • 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.7 KB
Line 
1!> @file biometeorology_mod.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM-4U.
4!
5! PALM-4U is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM-4U is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2018-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 3646 2018-12-28 17:58:49Z kanani $
29! Remove check for simulated_time > 0, it is not required since biometeorology
30! is only called from time_integration and not during time_integration_spinup
31!
32! 3614 2018-12-10 07:05:46Z raasch
33! unused variables removed
34!
35! 3593 2018-12-03 13:51:13Z kanani
36! Bugfix: additional tmrt_grid allocation in case bio_mrt not selected as ouput,
37! replace degree symbol by degree_C
38!
39! 3582 2018-11-29 19:16:36Z suehring
40! Consistently use bio_fill_value everywhere,
41! move allocation and initialization of output variables to bio_check_data_output
42! and bio_3d_data_averaging,
43! dom_dwd_user, Schrempf:
44! - integration of UVEM module part from r3474 (edited)
45! - split UTCI regression into 6 parts
46! - all data_output_3d is now explicity casted to sp
47!
48! 3525 2018-11-14 16:06:14Z kanani
49! Clean up, renaming from "biom" to "bio", summary of thermal index calculation
50! into one module (dom_dwd_user)
51!
52! 3479 2018-11-01 16:00:30Z gronemeier
53! - reworked check for output quantities
54! - assign new palm-error numbers
55! - set unit according to data standard.
56!
57! 3475 2018-10-30 21:16:31Z kanani
58! Add option for using MRT from RTM instead of simplified global value
59!
60! 3464 2018-10-30 18:08:55Z kanani
61! From branch resler@3462, pavelkrc:
62! make use of basic_constants_and_equations_mod
63!
64! 3448 2018-10-29 18:14:31Z kanani
65! Initial revision
66!
67!
68!
69! Authors:
70! --------
71! @author Dominik Froehlich <dominik.froehlich@dwd.de>, thermal indices
72! @author Jaroslav Resler <resler@cs.cas.cz>, mean radiant temperature
73! @author Michael Schrempf <schrempf@muk.uni-hannover.de>, uv exposure
74!
75!
76! Description:
77! ------------
78!> Biometeorology module consisting of two parts:
79!> 1.: Human thermal comfort module calculating thermal perception of a sample
80!> human being under the current meteorological conditions.
81!> 2.: Calculation of vitamin-D weighted UV exposure
82!>
83!> @todo Alphabetical sorting of "USE ..." lists, "ONLY" list, variable declarations
84!>       (per subroutine: first all CHARACTERs, then INTEGERs, LOGICALs, REALs, )
85!> @todo Comments start with capital letter --> "!-- Include..."
86!> @todo uv_vitd3dose-->new output type necessary (cumulative)
87!> @todo consider upwelling radiation in UV
88!>
89!> @note nothing now
90!>
91!> @bug  no known bugs by now
92!------------------------------------------------------------------------------!
93 MODULE biometeorology_mod
94
95    USE arrays_3d,                                                             &
96       ONLY:  pt, p, u, v, w, q
97
98    USE averaging,                                                             &
99       ONLY:  pt_av, q_av, u_av, v_av, w_av
100
101    USE basic_constants_and_equations_mod,                                     &
102       ONLY:  c_p, degc_to_k, l_v, magnus, sigma_sb, pi
103
104    USE control_parameters,                                                    &
105       ONLY:  average_count_3d, biometeorology, dz, dz_stretch_factor,         &
106              dz_stretch_level, humidity, initializing_actions, nz_do3d,       &
107              surface_pressure
108
109    USE date_and_time_mod,                                                                                                        &
110        ONLY:  calc_date_and_time, day_of_year, time_utc
111
112    USE grid_variables,                                                        &
113       ONLY:  ddx, dx, ddy, dy
114
115    USE indices,                                                               &
116       ONLY:  nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr, nxlg, nxrg,    &
117              nysg, nyng
118
119    USE kinds  !< Set precision of INTEGER and REAL arrays according to PALM
120
121    USE netcdf_data_input_mod,                                                                                                    &
122        ONLY:  netcdf_data_input_uvem, uvem_projarea_f, uvem_radiance_f,                                                          &
123               uvem_irradiance_f, uvem_integration_f, building_obstruction_f
124!
125!-- Import radiation model to obtain input for mean radiant temperature
126    USE radiation_model_mod,                                                   &
127       ONLY: ix, iy, iz, id, mrt_nlevels, mrt_include_sw,                      &
128             mrtinsw, mrtinlw, mrtbl, nmrtbl, radiation,                       &
129             radiation_interactions, rad_sw_in,                                &
130             rad_sw_out, rad_lw_in, rad_lw_out
131
132    USE surface_mod,                                                           &
133       ONLY: get_topography_top_index_ji
134
135    IMPLICIT NONE
136
137    PRIVATE
138
139!
140!-- Declare all global variables within the module (alphabetical order)
141    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tmrt_grid  !< tmrt results (degree_C)
142    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  perct      !< PT results   (degree_C)
143    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci       !< UTCI results (degree_C)
144    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet        !< PET results  (degree_C)
145!
146!-- Grids for averaged thermal indices
147    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  mrt_av_grid   !< time average mean
148    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  perct_av      !< PT results (aver. input)   (degree_C)
149    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci_av       !< UTCI results (aver. input) (degree_C)
150    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet_av        !< PET results (aver. input)  (degree_C)
151
152
153    INTEGER( iwp ) ::  bio_cell_level     !< cell level biom calculates for
154    REAL ( wp )    ::  bio_output_height  !< height output is calculated in m
155    REAL ( wp )    ::  time_bio_results   !< the time, the last set of biom results have been calculated for
156    REAL ( wp ), PARAMETER ::  human_absorb = 0.7_wp  !< SW absorbtivity of a human body (Fanger 1972)
157    REAL ( wp ), PARAMETER ::  human_emiss = 0.97_wp  !< LW emissivity of a human body after (Fanger 1972)
158    REAL ( wp ), PARAMETER ::  bio_fill_value = -9999._wp  !< set module fill value, replace by global fill value as soon as available
159!
160!--
161    LOGICAL ::  aver_perct = .FALSE.  !< switch: do perct averaging in this module? (if .FALSE. this is done globally)
162    LOGICAL ::  aver_q     = .FALSE.  !< switch: do e  averaging in this module?
163    LOGICAL ::  aver_u     = .FALSE.  !< switch: do u  averaging in this module?
164    LOGICAL ::  aver_v     = .FALSE.  !< switch: do v  averaging in this module?
165    LOGICAL ::  aver_w     = .FALSE.  !< switch: do w  averaging in this module?
166    LOGICAL ::  average_trigger_perct = .FALSE.  !< update averaged input on call to bio_perct?
167    LOGICAL ::  average_trigger_utci  = .FALSE.  !< update averaged input on call to bio_utci?
168    LOGICAL ::  average_trigger_pet   = .FALSE.  !< update averaged input on call to bio_pet?
169
170    LOGICAL ::  thermal_comfort = .TRUE.  !< Turn all thermal indices on or off
171    LOGICAL ::  bio_perct     = .TRUE.   !< Turn index PT (instant. input) on or off
172    LOGICAL ::  bio_perct_av  = .TRUE.   !< Turn index PT (averaged input) on or off
173    LOGICAL ::  bio_pet       = .TRUE.   !< Turn index PET (instant. input) on or off
174    LOGICAL ::  bio_pet_av    = .TRUE.   !< Turn index PET (averaged input) on or off
175    LOGICAL ::  bio_utci      = .TRUE.   !< Turn index UTCI (instant. input) on or off
176    LOGICAL ::  bio_utci_av   = .TRUE.   !< Turn index UTCI (averaged input) on or off
177
178!
179!-- UVEM parameters from here
180!
181!-- Declare all global variables within the module (alphabetical order)
182    INTEGER(iwp) ::  ai                      = 0  !< loop index in azimuth direction
183    INTEGER(iwp) ::  bi                      = 0  !< loop index of bit location within an 8bit-integer (one Byte)
184    INTEGER(iwp) ::  clothing                = 1  !< clothing (0=unclothed, 1=Arms,Hands,Face free, 3=Hand,Face free)
185    INTEGER(iwp) ::  iq                      = 0  !< loop index of irradiance quantity
186    INTEGER(iwp) ::  pobi                    = 0  !< loop index of the position of corresponding byte within ibset byte vektor
187    INTEGER(iwp) ::  obstruction_direct_beam = 0  !< Obstruction information for direct beam   
188    INTEGER(iwp) ::  zi                      = 0  !< loop index in zenith direction
189
190    INTEGER(KIND=1), DIMENSION(0:44)  ::  obstruction_temp1 = 0  !< temporary obstruction information stored with ibset
191    INTEGER(iwp),    DIMENSION(0:359) ::  obstruction_temp2 = 0  !< restored temporary obstruction information from ibset file
192
193    INTEGER(iwp), DIMENSION(0:35,0:9) ::  obstruction       = 1  !< final 2D obstruction information array
194
195    LOGICAL ::  consider_obstructions = .TRUE.   !< namelist parameter (see documentation)
196    LOGICAL ::  sun_in_south          = .FALSE.  !< namelist parameter (see documentation)
197    LOGICAL ::  turn_to_sun           = .TRUE.   !< namelist parameter (see documentation)
198    LOGICAL ::  uv_exposure           = .FALSE.  !< namelist parameter (see documentation)
199
200    REAL(wp) ::  diffuse_exposure            =   0.0_wp  !< calculated exposure by diffuse radiation
201    REAL(wp) ::  direct_exposure             =   0.0_wp  !< calculated exposure by direct solar beam   
202    REAL(wp) ::  orientation_angle           =   0.0_wp  !< orientation of front/face of the human model   
203    REAL(wp) ::  projection_area_direct_beam =   0.0_wp  !< projection area for direct solar beam
204    REAL(wp) ::  saa                         = 180.0_wp  !< solar azimuth angle
205    REAL(wp) ::  startpos_human              =   0.0_wp  !< start value for azimuth interpolation of human geometry array
206    REAL(wp) ::  startpos_saa_float          =   0.0_wp  !< start value for azimuth interpolation of radiance array
207    REAL(wp) ::  sza                         =  20.0_wp  !< solar zenith angle
208    REAL(wp) ::  xfactor                     =   0.0_wp  !< relative x-position used for interpolation
209    REAL(wp) ::  yfactor                     =   0.0_wp  !< relative y-position used for interpolation
210
211    REAL(wp), DIMENSION(0:2)  ::  irradiance =   0.0_wp  !< iradiance values extracted from irradiance lookup table 
212
213    REAL(wp), DIMENSION(0:2,0:90) ::  irradiance_lookup_table      = 0.0_wp  !< irradiance lookup table
214    REAL(wp), DIMENSION(0:35,0:9) ::  integration_array            = 0.0_wp  !< solid angle factors for hemispherical integration
215    REAL(wp), DIMENSION(0:35,0:9) ::  projection_area              = 0.0_wp  !< projection areas of a human (all directions)
216    REAL(wp), DIMENSION(0:35,0:9) ::  projection_area_lookup_table = 0.0_wp  !< human geometry lookup table (projection areas)
217    REAL(wp), DIMENSION(0:71,0:9) ::  projection_area_direct_temp  = 0.0_wp  !< temporary projection area for direct solar beam
218    REAL(wp), DIMENSION(0:71,0:9) ::  projection_area_temp         = 0.0_wp  !< temporary projection area for all directions
219    REAL(wp), DIMENSION(0:35,0:9) ::  radiance_array               = 0.0_wp  !< radiance extracted from radiance_lookup_table 
220    REAL(wp), DIMENSION(0:71,0:9) ::  radiance_array_temp          = 0.0_wp  !< temporary radiance data
221
222    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vitd3_exposure     !< result variable for instantaneous vitamin-D weighted exposures
223    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vitd3_exposure_av  !< result variable for summation of vitamin-D weighted exposures
224
225    REAL(wp), DIMENSION(0:35,0:9,0:90) ::  radiance_lookup_table   = 0.0_wp  !< radiance lookup table
226
227!
228!-- INTERFACES that must be available to other modules (alphabetical order)
229
230    PUBLIC bio_3d_data_averaging, bio_check_data_output,                       &
231    bio_calculate_mrt_grid, bio_calculate_thermal_index_maps, bio_calc_ipt,    &
232    bio_check_parameters, bio_data_output_3d, bio_data_output_2d,              &
233    bio_define_netcdf_grid, bio_get_thermal_index_input_ij, bio_header,        &
234    bio_init, bio_parin, bio_perct, bio_perct_av, bio_pet,                     &
235    bio_pet_av, bio_utci, bio_utci_av, thermal_comfort, time_bio_results,      &
236!
237!-- UVEM PUBLIC variables and methods
238    uvem_calc_exposure, uv_exposure
239
240!
241!-- PALM interfaces:
242!
243!-- 3D averaging for HTCM _INPUT_ variables
244    INTERFACE bio_3d_data_averaging
245       MODULE PROCEDURE bio_3d_data_averaging
246    END INTERFACE bio_3d_data_averaging
247!
248!-- Calculate mtr from rtm fluxes and assign into 2D grid
249    INTERFACE bio_calculate_mrt_grid
250       MODULE PROCEDURE bio_calculate_mrt_grid
251    END INTERFACE bio_calculate_mrt_grid
252!
253!-- Calculate static thermal indices PT, UTCI and/or PET
254    INTERFACE bio_calculate_thermal_index_maps
255       MODULE PROCEDURE bio_calculate_thermal_index_maps
256    END INTERFACE bio_calculate_thermal_index_maps
257!
258!-- Calculate the dynamic index iPT (to be caled by the agent model)
259    INTERFACE bio_calc_ipt
260       MODULE PROCEDURE bio_calc_ipt
261    END INTERFACE bio_calc_ipt
262!
263!-- Data output checks for 2D/3D data to be done in check_parameters
264    INTERFACE bio_check_data_output
265       MODULE PROCEDURE bio_check_data_output
266    END INTERFACE bio_check_data_output
267!
268!-- Input parameter checks to be done in check_parameters
269    INTERFACE bio_check_parameters
270       MODULE PROCEDURE bio_check_parameters
271    END INTERFACE bio_check_parameters
272!
273!-- Data output of 2D quantities
274    INTERFACE bio_data_output_2d
275       MODULE PROCEDURE bio_data_output_2d
276    END INTERFACE bio_data_output_2d
277!
278!-- no 3D data, thus, no averaging of 3D data, removed
279    INTERFACE bio_data_output_3d
280       MODULE PROCEDURE bio_data_output_3d
281    END INTERFACE bio_data_output_3d
282!
283!-- Definition of data output quantities
284    INTERFACE bio_define_netcdf_grid
285       MODULE PROCEDURE bio_define_netcdf_grid
286    END INTERFACE bio_define_netcdf_grid
287!
288!-- Obtains all relevant input values to estimate local thermal comfort/stress
289    INTERFACE bio_get_thermal_index_input_ij
290       MODULE PROCEDURE bio_get_thermal_index_input_ij
291    END INTERFACE bio_get_thermal_index_input_ij
292!
293!-- Output of information to the header file
294    INTERFACE bio_header
295       MODULE PROCEDURE bio_header
296    END INTERFACE bio_header
297!
298!-- Initialization actions
299    INTERFACE bio_init
300       MODULE PROCEDURE bio_init
301    END INTERFACE bio_init
302!
303!-- Reading of NAMELIST parameters
304    INTERFACE bio_parin
305       MODULE PROCEDURE bio_parin
306    END INTERFACE bio_parin
307
308!
309!-- Calculate UV exposure grid
310    INTERFACE uvem_calc_exposure
311       MODULE PROCEDURE uvem_calc_exposure
312    END INTERFACE uvem_calc_exposure
313
314 CONTAINS
315
316
317!------------------------------------------------------------------------------!
318! Description:
319! ------------
320!> Sum up and time-average biom input quantities as well as allocate
321!> the array necessary for storing the average.
322!> There is a considerable difference to the 3d_data_averaging subroutines
323!> used by other modules:
324!> For the thermal indices, the module needs to average the input conditions
325!> not the result!
326!------------------------------------------------------------------------------!
327 SUBROUTINE bio_3d_data_averaging( mode, variable )
328
329    IMPLICIT NONE
330
331    CHARACTER (LEN=*) ::  mode     !< averaging mode: allocate, sum, or average
332    CHARACTER (LEN=*) ::  variable !< The variable in question
333
334    INTEGER(iwp) ::  i        !< Running index, x-dir
335    INTEGER(iwp) ::  j        !< Running index, y-dir
336    INTEGER(iwp) ::  k        !< Running index, z-dir
337
338
339    IF ( mode == 'allocate' )  THEN
340
341       SELECT CASE ( TRIM( variable ) )
342
343          CASE ( 'bio_mrt' )
344                IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
345                   ALLOCATE( mrt_av_grid(nmrtbl) )
346                ENDIF
347                mrt_av_grid = 0.0_wp
348
349          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
350
351!
352!--          Averaging, as well as the allocation of the required grids must be
353!            done only once, independent from for how many thermal indices
354!            averaged output is desired.
355!            Therefore wee need to memorize which index is the one that controls
356!            the averaging (what must be the first thermal index called).
357!            Indices are in unknown order as depending on the input file,
358!            determine first index to average und update only once
359
360!--          Only proceed here if this was not done for any index before. This
361!            is done only once during the whole model run.
362             IF ( .NOT. average_trigger_perct .AND.                            &
363                  .NOT. average_trigger_utci  .AND.                            &
364                  .NOT. average_trigger_pet ) THEN
365!
366!--             Allocate the required grids
367                IF ( .NOT. ALLOCATED( perct_av ) ) THEN
368                   ALLOCATE( perct_av (nys:nyn,nxl:nxr) )
369                ENDIF
370                perct_av = REAL( bio_fill_value, KIND = wp )
371
372                IF ( .NOT. ALLOCATED( utci_av ) ) THEN
373                   ALLOCATE( utci_av (nys:nyn,nxl:nxr) )
374                ENDIF
375                utci_av = REAL( bio_fill_value, KIND = wp )
376
377                IF ( .NOT. ALLOCATED( pet_av ) ) THEN
378                   ALLOCATE( pet_av (nys:nyn,nxl:nxr) )
379                ENDIF
380                pet_av = REAL( bio_fill_value, KIND = wp )
381!
382!--             Memorize the first index called to control averaging
383                IF ( TRIM( variable ) == 'bio_perct*' ) THEN
384                    average_trigger_perct = .TRUE.
385                ENDIF
386                IF ( TRIM( variable ) == 'bio_utci*' ) THEN
387                    average_trigger_utci = .TRUE.
388                ENDIF
389                IF ( TRIM( variable ) == 'bio_pet*' ) THEN
390                    average_trigger_pet = .TRUE.
391                ENDIF
392             ENDIF
393!
394!--          Only continue if var is the index, that controls averaging.
395!            Break immediatelly (doing nothing) for the other indices.
396             IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') RETURN
397             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*')   RETURN
398             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'bio_pet*')    RETURN
399!
400!--          Now memorize which of the input grids are not averaged by other
401!            modules. Set averaging switch to .TRUE. in that case.
402             IF ( .NOT. ALLOCATED( pt_av ) )  THEN
403                ALLOCATE( pt_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
404                aver_perct = .TRUE.
405                pt_av = 0.0_wp
406             ENDIF
407
408             IF ( .NOT. ALLOCATED( q_av ) )  THEN
409                ALLOCATE( q_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
410                aver_q = .TRUE.
411                q_av = 0.0_wp
412             ENDIF
413
414             IF ( .NOT. ALLOCATED( u_av ) )  THEN
415                ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
416                aver_u = .TRUE.
417                u_av = 0.0_wp
418             ENDIF
419
420             IF ( .NOT. ALLOCATED( v_av ) )  THEN
421                ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
422                aver_v = .TRUE.
423                v_av = 0.0_wp
424             ENDIF
425
426             IF ( .NOT. ALLOCATED( w_av ) )  THEN
427                ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
428                aver_w = .TRUE.
429                w_av = 0.0_wp
430             ENDIF
431
432          CASE ( 'uvem_vitd3dose*' )
433             IF ( .NOT. ALLOCATED( vitd3_exposure_av ) )  THEN
434                ALLOCATE( vitd3_exposure_av(nysg:nyng,nxlg:nxrg) )
435             ENDIF
436             vitd3_exposure_av = 0.0_wp
437
438          CASE DEFAULT
439             CONTINUE
440
441       END SELECT
442
443    ELSEIF ( mode == 'sum' )  THEN
444
445       SELECT CASE ( TRIM( variable ) )
446
447          CASE ( 'bio_mrt' )
448             IF ( ALLOCATED( mrt_av_grid ) )  THEN
449
450                IF ( mrt_include_sw )  THEN
451                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
452                      (( human_absorb * mrtinsw(:) + human_emiss * mrtinlw(:)) &
453                      / (human_emiss * sigma_sb)) ** .25_wp - degc_to_k
454                ELSE
455                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
456                      (human_emiss * mrtinlw(:) / sigma_sb) ** .25_wp          &
457                      - degc_to_k
458                ENDIF
459             ENDIF
460
461          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
462!
463!--          Only continue if updateindex, see above
464             IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') &
465                RETURN
466             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*')  &
467                RETURN
468             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'bio_pet*')   &
469                RETURN
470
471             IF ( ALLOCATED( pt_av ) .AND. aver_perct ) THEN
472                DO  i = nxl, nxr
473                   DO  j = nys, nyn
474                      DO  k = nzb, nzt+1
475                         pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i)
476                      ENDDO
477                   ENDDO
478                ENDDO
479             ENDIF
480
481             IF ( ALLOCATED( q_av )  .AND. aver_q ) THEN
482                DO  i = nxl, nxr
483                   DO  j = nys, nyn
484                      DO  k = nzb, nzt+1
485                         q_av(k,j,i) = q_av(k,j,i) + q(k,j,i)
486                      ENDDO
487                   ENDDO
488                ENDDO
489             ENDIF
490
491             IF ( ALLOCATED( u_av )  .AND. aver_u ) THEN
492                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
493                   DO  j = nysg, nyng
494                      DO  k = nzb, nzt+1
495                         u_av(k,j,i) = u_av(k,j,i) + u(k,j,i)
496                      ENDDO
497                   ENDDO
498                ENDDO
499             ENDIF
500
501             IF ( ALLOCATED( v_av )  .AND. aver_v ) THEN
502                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
503                   DO  j = nysg, nyng
504                      DO  k = nzb, nzt+1
505                         v_av(k,j,i) = v_av(k,j,i) + v(k,j,i)
506                      ENDDO
507                   ENDDO
508                ENDDO
509             ENDIF
510
511             IF ( ALLOCATED( w_av )  .AND. aver_w ) THEN
512                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
513                   DO  j = nysg, nyng
514                      DO  k = nzb, nzt+1
515                         w_av(k,j,i) = w_av(k,j,i) + w(k,j,i)
516                      ENDDO
517                   ENDDO
518                ENDDO
519             ENDIF
520!
521!--       This is a cumulated dose. No mode == 'average' for this quantity.
522          CASE ( 'uvem_vitd3dose*' )
523             IF ( ALLOCATED( vitd3_exposure_av ) ) THEN
524                DO  i = nxlg, nxrg
525                   DO  j = nysg, nyng
526                      vitd3_exposure_av(j,i) = vitd3_exposure_av(j,i) + vitd3_exposure(j,i)
527                   ENDDO
528                ENDDO
529             ENDIF
530
531          CASE DEFAULT
532             CONTINUE
533
534       END SELECT
535
536    ELSEIF ( mode == 'average' )  THEN
537
538       SELECT CASE ( TRIM( variable ) )
539
540          CASE ( 'bio_mrt' )
541             IF ( ALLOCATED( mrt_av_grid ) )  THEN
542                mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d, KIND=wp )
543             ENDIF
544
545          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
546!
547!--          Only continue if update index, see above
548             IF ( average_trigger_perct .AND.                                  &
549                TRIM( variable ) /= 'bio_perct*' ) RETURN
550             IF ( average_trigger_utci .AND.                                   &
551                TRIM( variable ) /= 'bio_utci*' ) RETURN
552             IF ( average_trigger_pet  .AND.                                   &
553                TRIM( variable ) /= 'bio_pet*' ) RETURN
554
555             IF ( ALLOCATED( pt_av ) .AND. aver_perct ) THEN
556                DO  i = nxl, nxr
557                   DO  j = nys, nyn
558                      DO  k = nzb, nzt+1
559                         pt_av(k,j,i) = pt_av(k,j,i) /                         &
560                            REAL( average_count_3d, KIND=wp )
561                      ENDDO
562                   ENDDO
563                ENDDO
564             ENDIF
565
566             IF ( ALLOCATED( q_av ) .AND. aver_q ) THEN
567                DO  i = nxl, nxr
568                   DO  j = nys, nyn
569                      DO  k = nzb, nzt+1
570                         q_av(k,j,i) = q_av(k,j,i) /                           &
571                            REAL( average_count_3d, KIND=wp )
572                      ENDDO
573                   ENDDO
574                ENDDO
575             ENDIF
576
577             IF ( ALLOCATED( u_av ) .AND. aver_u ) THEN
578                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
579                   DO  j = nysg, nyng
580                      DO  k = nzb, nzt+1
581                         u_av(k,j,i) = u_av(k,j,i) /                           &
582                            REAL( average_count_3d, KIND=wp )
583                      ENDDO
584                   ENDDO
585                ENDDO
586             ENDIF
587
588             IF ( ALLOCATED( v_av ) .AND. aver_v ) THEN
589                DO  i = nxlg, nxrg
590                   DO  j = nysg, nyng
591                      DO  k = nzb, nzt+1
592                         v_av(k,j,i) = v_av(k,j,i) /                           &
593                            REAL( average_count_3d, KIND=wp )
594                      ENDDO
595                   ENDDO
596                ENDDO
597             ENDIF
598
599             IF ( ALLOCATED( w_av ) .AND. aver_w ) THEN
600                DO  i = nxlg, nxrg
601                   DO  j = nysg, nyng
602                      DO  k = nzb, nzt+1
603                         w_av(k,j,i) = w_av(k,j,i) /                           &
604                            REAL( average_count_3d, KIND=wp )
605                      ENDDO
606                   ENDDO
607                ENDDO
608             ENDIF
609!
610!--          Udate all thermal index grids with updated averaged input
611             CALL bio_calculate_thermal_index_maps ( .TRUE. )
612
613!
614!--     No averaging for UVEM since we are calculating a dose (only sum is
615!--     calculated and saved to av.nc file)
616
617        END SELECT
618
619    ENDIF
620
621
622 END SUBROUTINE bio_3d_data_averaging
623
624
625
626!------------------------------------------------------------------------------!
627! Description:
628! ------------
629!> Check data output for biometeorology model
630!------------------------------------------------------------------------------!
631 SUBROUTINE bio_check_data_output( var, unit, i, ilen, k )
632
633    USE control_parameters,                                                    &
634        ONLY: data_output, message_string
635
636    IMPLICIT NONE
637
638    CHARACTER (LEN=*) ::  unit     !< The unit for the variable var
639    CHARACTER (LEN=*) ::  var      !< The variable in question
640
641    INTEGER(iwp) ::  i      !<
642    INTEGER(iwp) ::  ilen   !<   
643    INTEGER(iwp) ::  k      !<
644
645    SELECT CASE ( TRIM( var ) )
646!
647!-- Allocate a temporary array with the desired output dimensions.
648       CASE ( 'bio_mrt' )
649          unit = 'degree_C'
650          IF ( .NOT. ALLOCATED( tmrt_grid ) )  THEN
651             ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) )
652          ENDIF
653          tmrt_grid = REAL( bio_fill_value, KIND = wp )
654
655       CASE ( 'bio_perct*' )
656          unit = 'degree_C'
657          IF ( .NOT. ALLOCATED( perct ) )  THEN
658             ALLOCATE( perct (nys:nyn,nxl:nxr) )
659          ENDIF
660          perct = REAL( bio_fill_value, KIND = wp )
661
662       CASE ( 'bio_utci*' )
663          unit = 'degree_C'
664          IF ( .NOT. ALLOCATED( utci ) )  THEN
665             ALLOCATE( utci (nys:nyn,nxl:nxr) )
666          ENDIF
667          utci = REAL( bio_fill_value, KIND = wp )
668
669       CASE ( 'bio_pet*' )
670          unit = 'degree_C'
671          IF ( .NOT. ALLOCATED( pet ) )  THEN
672             ALLOCATE( pet (nys:nyn,nxl:nxr) )
673          ENDIF
674          pet = REAL( bio_fill_value, KIND = wp )
675
676       CASE ( 'uvem_vitd3*' )
677          IF (  .NOT.  uv_exposure )  THEN
678             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
679                      'res a namelist &uvexposure_par'
680             CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 )
681          ENDIF
682          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
683             message_string = 'illegal value for data_output: "' //            &
684                              TRIM( var ) // '" & only 2d-horizontal ' //      &
685                              'cross sections are allowed for this value'
686             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
687          ENDIF
688          unit = 'IU/s'
689          IF ( .NOT. ALLOCATED( vitd3_exposure ) )  THEN
690             ALLOCATE( vitd3_exposure(nysg:nyng,nxlg:nxrg) )
691          ENDIF
692          vitd3_exposure = 0.0_wp
693
694       CASE ( 'uvem_vitd3dose*' )
695          IF (  .NOT.  uv_exposure )  THEN
696             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
697                      'res  a namelist &uvexposure_par'
698             CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 )
699          ENDIF
700          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
701             message_string = 'illegal value for data_output: "' //            &
702                              TRIM( var ) // '" & only 2d-horizontal ' //      &
703                              'cross sections are allowed for this value'
704             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
705          ENDIF
706          unit = 'IU/av-h'
707          IF ( .NOT. ALLOCATED( vitd3_exposure_av ) )  THEN
708             ALLOCATE( vitd3_exposure_av(nysg:nyng,nxlg:nxrg) )
709          ENDIF
710          vitd3_exposure_av = 0.0_wp
711
712       CASE DEFAULT
713          unit = 'illegal'
714
715    END SELECT
716
717    IF ( thermal_comfort .AND. unit == 'degree_C' )  THEN
718!
719!--    Futher checks if output belongs to biometeorology
720!      Break if required modules "radiation" and "humidity" are not running.
721       IF ( .NOT.  radiation ) THEN
722          message_string = 'output of "' // TRIM( var ) // '" require'         &
723                           // 's radiation = .TRUE.'
724          CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
725          unit = 'illegal'
726       ENDIF
727       IF ( mrt_nlevels == 0 ) THEN
728          message_string = 'output of "' // TRIM( var ) // '" require'         &
729                           // 's mrt_nlevels > 0'
730          CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
731          unit = 'illegal'
732       ENDIF
733
734
735    ENDIF
736
737 END SUBROUTINE bio_check_data_output
738
739!------------------------------------------------------------------------------!
740! Description:
741! ------------
742!> Check parameters routine for biom module
743!> check_parameters 1302
744!------------------------------------------------------------------------------!
745 SUBROUTINE bio_check_parameters
746
747    USE control_parameters,                                                    &
748        ONLY:  message_string
749
750
751    IMPLICIT NONE
752
753!
754!--    Disabled as long as radiation model not available
755
756       IF ( thermal_comfort  .AND.  .NOT.  humidity )  THEN
757          message_string = 'The estimation of thermal comfort requires '    // &
758                           'air humidity information, but humidity module ' // &
759                           'is disabled!'
760          CALL message( 'check_parameters', 'PA0561', 0, 0, 0, 6, 0 )
761       ENDIF
762
763 END SUBROUTINE bio_check_parameters
764
765
766!------------------------------------------------------------------------------!
767!
768! Description:
769! ------------
770!> Subroutine defining 2D output variables
771!> data_output_2d 1188ff
772!------------------------------------------------------------------------------!
773 SUBROUTINE bio_data_output_2d( av, variable, found, grid, local_pf,           &
774                                two_d, nzb_do, nzt_do )
775
776
777    USE kinds
778
779
780    IMPLICIT NONE
781!
782!-- Input variables
783    CHARACTER (LEN=*), INTENT(IN) ::  variable    !< Char identifier to select var for output
784    INTEGER(iwp), INTENT(IN)      ::  av          !< Use averaged data? 0 = no, 1 = yes?
785    INTEGER(iwp), INTENT(IN)      ::  nzb_do      !< Unused. 2D. nz bottom to nz top
786    INTEGER(iwp), INTENT(IN)      ::  nzt_do      !< Unused.
787!
788!-- Output variables
789    CHARACTER (LEN=*), INTENT(OUT) ::  grid   !< Grid type (always "zu1" for biom)
790    LOGICAL, INTENT(OUT)           ::  found  !< Output found?
791    LOGICAL, INTENT(OUT)           ::  two_d  !< Flag parameter that indicates 2D variables, horizontal cross sections, must be .TRUE.
792    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !< Temp. result grid to return
793!
794!-- Internal variables
795    INTEGER(iwp) ::  i        !< Running index, x-dir
796    INTEGER(iwp) ::  j        !< Running index, y-dir
797    INTEGER(iwp) ::  k        !< Running index, z-dir
798    INTEGER(iwp) ::  l        !< Running index, radiation grid
799
800
801    found = .TRUE.
802    local_pf = bio_fill_value
803
804    SELECT CASE ( TRIM( variable ) )
805
806
807        CASE ( 'bio_mrt_xy' )
808           grid = 'zu1'
809           two_d = .FALSE.  !< can be calculated for several levels
810           local_pf = REAL( bio_fill_value, KIND = wp )
811           DO  l = 1, nmrtbl
812              i = mrtbl(ix,l)
813              j = mrtbl(iy,l)
814              k = mrtbl(iz,l)
815              IF ( k < nzb_do .OR. k > nzt_do .OR. j < nys .OR. j > nyn .OR.   &
816                 i < nxl .OR. i > nxr ) CYCLE
817              IF ( av == 0 )  THEN
818                 IF ( mrt_include_sw )  THEN
819                    local_pf(i,j,k) = ((human_absorb * mrtinsw(l) +            &
820                    human_emiss * mrtinlw(l)) /                                &
821                    (human_emiss * sigma_sb)) ** .25_wp - degc_to_k
822                 ELSE
823                    local_pf(i,j,k) = (human_emiss * mrtinlw(l) /              &
824                                       sigma_sb) ** .25_wp - degc_to_k
825                 ENDIF
826              ELSE
827                 local_pf(i,j,k) = mrt_av_grid(l)
828              ENDIF
829           ENDDO
830
831
832        CASE ( 'bio_perct*_xy' )        ! 2d-array
833           grid = 'zu1'
834           two_d = .TRUE.
835           IF ( av == 0 )  THEN
836              DO  i = nxl, nxr
837                 DO  j = nys, nyn
838                    local_pf(i,j,nzb+1) = perct(j,i)
839                 ENDDO
840              ENDDO
841           ELSE
842              DO  i = nxl, nxr
843                 DO  j = nys, nyn
844                    local_pf(i,j,nzb+1) = perct_av(j,i)
845                 ENDDO
846              ENDDO
847           END IF
848
849
850        CASE ( 'bio_utci*_xy' )        ! 2d-array
851           grid = 'zu1'
852           two_d = .TRUE.
853           IF ( av == 0 )  THEN
854              DO  i = nxl, nxr
855                 DO  j = nys, nyn
856                    local_pf(i,j,nzb+1) = utci(j,i)
857                 ENDDO
858              ENDDO
859           ELSE
860              DO  i = nxl, nxr
861                 DO  j = nys, nyn
862                    local_pf(i,j,nzb+1) = utci_av(j,i)
863                 ENDDO
864              ENDDO
865           END IF
866
867
868        CASE ( 'bio_pet*_xy' )        ! 2d-array
869           grid = 'zu1'
870           two_d = .TRUE.
871           IF ( av == 0 )  THEN
872              DO  i = nxl, nxr
873                 DO  j = nys, nyn
874                    local_pf(i,j,nzb+1) = pet(j,i)
875                 ENDDO
876              ENDDO
877           ELSE
878              DO  i = nxl, nxr
879                 DO  j = nys, nyn
880                    local_pf(i,j,nzb+1) = pet_av(j,i)
881                 ENDDO
882              ENDDO
883           END IF
884
885           !
886!--    Before data is transfered to local_pf, transfer is it 2D dummy variable and exchange ghost points therein.
887!--    However, at this point this is only required for instantaneous arrays, time-averaged quantities are already exchanged.
888       CASE ( 'uvem_vitd3*_xy' )        ! 2d-array
889          IF ( av == 0 )  THEN
890             DO  i = nxl, nxr
891                DO  j = nys, nyn
892                   local_pf(i,j,nzb+1) = vitd3_exposure(j,i)
893                ENDDO
894             ENDDO
895          ENDIF
896
897          two_d = .TRUE.
898          grid = 'zu1'
899
900       CASE ( 'uvem_vitd3dose*_xy' )        ! 2d-array
901          IF ( av == 1 )  THEN
902             DO  i = nxl, nxr
903                DO  j = nys, nyn
904                   local_pf(i,j,nzb+1) = vitd3_exposure_av(j,i)
905                ENDDO
906             ENDDO
907          ENDIF
908
909          two_d = .TRUE.
910          grid = 'zu1'
911
912
913       CASE DEFAULT
914          found = .FALSE.
915          grid  = 'none'
916
917    END SELECT
918
919
920 END SUBROUTINE bio_data_output_2d
921
922
923!------------------------------------------------------------------------------!
924!
925! Description:
926! ------------
927!> Subroutine defining 3D output variables (dummy, always 2d!)
928!> data_output_3d 709ff
929!------------------------------------------------------------------------------!
930 SUBROUTINE bio_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
931
932    USE indices
933
934    USE kinds
935
936
937    IMPLICIT NONE
938!
939!-- Input variables
940    CHARACTER (LEN=*), INTENT(IN) ::  variable   !< Char identifier to select var for output
941    INTEGER(iwp), INTENT(IN) ::  av       !< Use averaged data? 0 = no, 1 = yes?
942    INTEGER(iwp), INTENT(IN) ::  nzb_do   !< Unused. 2D. nz bottom to nz top
943    INTEGER(iwp), INTENT(IN) ::  nzt_do   !< Unused.
944!
945!-- Output variables
946    LOGICAL, INTENT(OUT) ::  found   !< Output found?
947    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< Temp. result grid to return
948!
949!-- Internal variables
950    INTEGER(iwp) ::  l    !< Running index, radiation grid
951    INTEGER(iwp) ::  i    !< Running index, x-dir
952    INTEGER(iwp) ::  j    !< Running index, y-dir
953    INTEGER(iwp) ::  k    !< Running index, z-dir
954
955!     REAL(wp) ::  mrt  !< Buffer for mean radiant temperature
956
957    found = .TRUE.
958
959    SELECT CASE ( TRIM( variable ) )
960
961        CASE ( 'bio_mrt' )
962            local_pf = REAL( bio_fill_value, KIND = sp )
963            DO  l = 1, nmrtbl
964               i = mrtbl(ix,l)
965               j = mrtbl(iy,l)
966               k = mrtbl(iz,l)
967               IF ( k < nzb_do .OR. k > nzt_do .OR. j < nys .OR. j > nyn .OR.  &
968                  i < nxl .OR. i > nxr ) CYCLE
969               IF ( av == 0 )  THEN
970                  IF ( mrt_include_sw )  THEN
971                     local_pf(i,j,k) = REAL(((human_absorb * mrtinsw(l) +      &
972                                    human_emiss * mrtinlw(l)) /                &
973                                    (human_emiss * sigma_sb)) ** .25_wp -      &
974                                    degc_to_k, kind=sp )
975                  ELSE
976                     local_pf(i,j,k) = REAL((human_emiss * mrtinlw(l) /        &
977                                    sigma_sb) ** .25_wp - degc_to_k, kind=sp )  !< why not (human_emiss * sigma_sb) as above?
978                  ENDIF
979               ELSE
980                  local_pf(i,j,k) = REAL(mrt_av_grid(l), kind=sp)
981               ENDIF
982            ENDDO
983
984       CASE DEFAULT
985          found = .FALSE.
986
987    END SELECT
988
989 END SUBROUTINE bio_data_output_3d
990
991!------------------------------------------------------------------------------!
992! Description:
993! ------------
994!> Subroutine defining appropriate grid for netcdf variables.
995!> It is called out from subroutine netcdf_interface_mod.
996!> netcdf_interface_mod 918ff
997!------------------------------------------------------------------------------!
998 SUBROUTINE bio_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
999
1000    IMPLICIT NONE
1001!
1002!-- Input variables
1003    CHARACTER (LEN=*), INTENT(IN)  ::  var      !< Name of output variable
1004!
1005!-- Output variables
1006    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x   !< x grid of output variable
1007    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y   !< y grid of output variable
1008    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z   !< z grid of output variable
1009
1010    LOGICAL, INTENT(OUT)           ::  found    !< Flag if output var is found
1011!
1012!-- Local variables
1013    LOGICAL      :: is2d  !< Var is 2d?
1014
1015    INTEGER(iwp) :: l     !< Length of the var array
1016
1017
1018    found  = .FALSE.
1019    grid_x = 'none'
1020    grid_y = 'none'
1021    grid_z = 'none'
1022
1023    l = MAX( 2, LEN_TRIM( var ) )
1024    is2d = ( var(l-1:l) == 'xy' )
1025
1026    IF ( var(1:4) == 'bio_' ) THEN
1027       found  = .TRUE.
1028       grid_x = 'x'
1029       grid_y = 'y'
1030       grid_z = 'zu'
1031       IF ( is2d  .AND.  var(1:7) /= 'bio_mrt' ) grid_z = 'zu1'
1032    ENDIF
1033
1034    IF ( is2d  .AND.  var(1:4) == 'uvem' ) THEN
1035       grid_x = 'x'
1036       grid_y = 'y'
1037       grid_z = 'zu1'
1038    ENDIF
1039
1040 END SUBROUTINE bio_define_netcdf_grid
1041
1042!------------------------------------------------------------------------------!
1043! Description:
1044! ------------
1045!> Header output for biom module
1046!> header 982
1047!------------------------------------------------------------------------------!
1048 SUBROUTINE bio_header( io )
1049
1050    IMPLICIT NONE
1051!
1052!-- Input variables
1053    INTEGER(iwp), INTENT(IN) ::  io           !< Unit of the output file
1054!
1055!-- Internal variables
1056    CHARACTER (LEN=86) ::  output_height_chr  !< String for output height
1057
1058    WRITE( output_height_chr, '(F8.1,7X)' )  bio_output_height
1059!
1060!-- Write biom header
1061    WRITE( io, 1 )
1062    WRITE( io, 2 )  TRIM( output_height_chr )
1063    WRITE( io, 3 )  TRIM( ACHAR( bio_cell_level ) )
1064
10651   FORMAT (//' Human thermal comfort module information:'/                    &
1066              ' ------------------------------'/)
10672   FORMAT ('    --> All indices calculated for a height of (m): ', A )
10683   FORMAT ('    --> This corresponds to cell level : ', A )
1069
1070 END SUBROUTINE bio_header
1071
1072
1073!------------------------------------------------------------------------------!
1074! Description:
1075! ------------
1076!> Initialization of the HTCM
1077!> init_3d_model 1987ff
1078!------------------------------------------------------------------------------!
1079 SUBROUTINE bio_init
1080
1081    USE control_parameters,                                                    &
1082        ONLY: message_string
1083
1084    USE netcdf_data_input_mod,                                                 &
1085        ONLY:  netcdf_data_input_uvem
1086
1087    IMPLICIT NONE
1088!
1089!-- Internal vriables
1090    REAL ( wp )  :: height  !< current height in meters
1091!
1092!-- Determine cell level corresponding to 1.1 m above ground level
1093!   (gravimetric center of sample human)
1094
1095    time_bio_results = 0.0_wp
1096    bio_cell_level = 0_iwp
1097    bio_output_height = 0.5_wp * dz(1)
1098    height = 0.0_wp
1099
1100    bio_cell_level = INT ( 1.099_wp / dz(1) )
1101    bio_output_height = bio_output_height + bio_cell_level * dz(1)
1102
1103    IF ( .NOT. radiation_interactions ) THEN
1104       message_string = 'The mrt calculation requires ' //                     &
1105                        'enabled radiation_interactions but it ' //            &
1106                        'is disabled!'
1107       CALL message( 'check_parameters', 'PAHU03', 0, 0, -1, 6, 0 )
1108    ENDIF
1109
1110!
1111!-- Init UVEM and load lookup tables
1112    CALL netcdf_data_input_uvem
1113
1114 END SUBROUTINE bio_init
1115
1116
1117!------------------------------------------------------------------------------!
1118! Description:
1119! ------------
1120!> Parin for &biometeorology_parameters for reading biomet parameters
1121!------------------------------------------------------------------------------!
1122 SUBROUTINE bio_parin
1123
1124    IMPLICIT NONE
1125
1126!
1127!-- Internal variables
1128    CHARACTER (LEN=80) ::  line  !< Dummy string for current line in parameter file
1129
1130    NAMELIST /biometeorology_parameters/  bio_pet,                             &
1131                                          bio_pet_av,                          &
1132                                          bio_perct,                           &
1133                                          bio_perct_av,                        &
1134                                          bio_utci,                            &
1135                                          bio_utci_av,                         &
1136                                          thermal_comfort,                     &
1137!
1138!-- UVEM namelist parameters
1139                                          clothing,                            &
1140                                          consider_obstructions,               &
1141                                          orientation_angle,                   &
1142                                          sun_in_south,                        &
1143                                          turn_to_sun,                         &
1144                                          uv_exposure
1145
1146
1147!-- Try to find biometeorology_parameters namelist
1148    REWIND ( 11 )
1149    line = ' '
1150    DO WHILE ( INDEX( line, '&biometeorology_parameters' ) == 0 )
1151       READ ( 11, '(A)', END = 20 )  line
1152    ENDDO
1153    BACKSPACE ( 11 )
1154
1155!
1156!-- Read biometeorology_parameters namelist
1157    READ ( 11, biometeorology_parameters, ERR = 10, END = 20 )
1158
1159!
1160!-- Set flag that indicates that the biomet_module is switched on
1161    biometeorology = .TRUE.
1162
1163    GOTO 20
1164
1165!
1166!-- In case of error
1167 10 BACKSPACE( 11 )
1168    READ( 11 , '(A)') line
1169    CALL parin_fail_message( 'biometeorology_parameters', line )
1170
1171!
1172!-- Complete
1173 20 CONTINUE
1174
1175
1176 END SUBROUTINE bio_parin
1177
1178!------------------------------------------------------------------------------!
1179! Description:
1180! ------------
1181!> Calculate biometeorology MRT for all 2D grid
1182!------------------------------------------------------------------------------!
1183 SUBROUTINE bio_calculate_mrt_grid ( av )
1184
1185    IMPLICIT NONE
1186
1187    LOGICAL, INTENT(IN)         ::  av    !< use averaged input?
1188!
1189!-- Internal variables
1190    INTEGER(iwp)                ::  i     !< Running index, x-dir, radiation coordinates
1191    INTEGER(iwp)                ::  j     !< Running index, y-dir, radiation coordinates
1192    INTEGER(iwp)                ::  k     !< Running index, y-dir, radiation coordinates
1193    INTEGER(iwp)                ::  l     !< Running index, radiation coordinates
1194
1195!
1196!-- Calculate biometeorology MRT from local radiation fluxes calculated by RTM and assign
1197!-- into 2D grid. Depending on selected output quantities, tmrt_grid might not have been
1198!-- allocated in bio_check_data_output yet.
1199    IF ( .NOT. ALLOCATED( tmrt_grid ) )  THEN
1200       ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) )
1201    ENDIF
1202    tmrt_grid = REAL( bio_fill_value, KIND = wp )
1203
1204    DO  l = 1, nmrtbl
1205       i = mrtbl(ix,l)
1206       j = mrtbl(iy,l)
1207       k = mrtbl(iz,l)
1208       IF ( k - get_topography_top_index_ji( j, i, 's' ) == bio_cell_level +   &
1209             1_iwp) THEN
1210          IF ( mrt_include_sw )  THEN
1211              tmrt_grid(j,i) = ((human_absorb*mrtinsw(l) +                     &
1212                                human_emiss*mrtinlw(l))  /                     &
1213                               (human_emiss*sigma_sb)) ** .25_wp - degc_to_k
1214          ELSE
1215              tmrt_grid(j,i) = (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp   &
1216                                 - degc_to_k
1217          ENDIF
1218       ENDIF
1219    ENDDO
1220
1221END SUBROUTINE bio_calculate_mrt_grid
1222
1223
1224!------------------------------------------------------------------------------!
1225! Description:
1226! ------------
1227!> Calculate static thermal indices for 2D grid point i, j
1228!------------------------------------------------------------------------------!
1229 SUBROUTINE bio_get_thermal_index_input_ij( average_input, i, j, ta, vp, ws,   &
1230    pair, tmrt )
1231
1232    IMPLICIT NONE
1233!
1234!-- Input variables
1235    LOGICAL,      INTENT ( IN ) ::  average_input  !< Determine averaged input conditions?
1236    INTEGER(iwp), INTENT ( IN ) ::  i     !< Running index, x-dir
1237    INTEGER(iwp), INTENT ( IN ) ::  j     !< Running index, y-dir
1238!
1239!-- Output parameters
1240    REAL(wp), INTENT ( OUT )    ::  tmrt  !< Mean radiant temperature        (degree_C)
1241    REAL(wp), INTENT ( OUT )    ::  ta    !< Air temperature                 (degree_C)
1242    REAL(wp), INTENT ( OUT )    ::  vp    !< Vapour pressure                 (hPa)
1243    REAL(wp), INTENT ( OUT )    ::  ws    !< Wind speed    (local level)     (m/s)
1244    REAL(wp), INTENT ( OUT )    ::  pair  !< Air pressure                    (hPa)
1245!
1246!-- Internal variables
1247    INTEGER(iwp)                ::  k     !< Running index, z-dir
1248!     INTEGER(iwp)                ::  ir    !< Running index, x-dir, radiation coordinates
1249!     INTEGER(iwp)                ::  jr    !< Running index, y-dir, radiation coordinates
1250!     INTEGER(iwp)                ::  kr    !< Running index, y-dir, radiation coordinates
1251    INTEGER(iwp)                ::  k_wind  !< Running index, z-dir, wind speed only
1252!     INTEGER(iwp)                ::  l     !< Running index, radiation coordinates
1253
1254    REAL(wp)                    ::  vp_sat  !< Saturation vapor pressure     (hPa)
1255
1256!
1257!-- Determine cell level closest to 1.1m above ground
1258!   by making use of truncation due to int cast
1259    k = get_topography_top_index_ji(j, i, 's') + bio_cell_level  !< Vertical cell center closest to 1.1m
1260
1261!
1262!-- Avoid non-representative horizontal u and v of 0.0 m/s too close to ground
1263    k_wind = k
1264    IF ( bio_cell_level < 1_iwp ) THEN
1265       k_wind = k + 1_iwp
1266    ENDIF
1267!
1268!-- Determine local values:
1269    IF ( average_input ) THEN
1270!
1271!--    Calculate ta from Tp assuming dry adiabatic laps rate
1272       ta = pt_av(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
1273
1274       vp = bio_fill_value
1275       IF ( humidity .AND. ALLOCATED( q_av ) ) THEN
1276          vp = q_av(k, j, i)
1277       ENDIF
1278
1279       ws = ( 0.5_wp * ABS( u_av(k_wind, j, i) + u_av(k_wind, j, i+1) )  +     &
1280          0.5_wp * ABS( v_av(k_wind, j, i) + v_av(k_wind, j+1, i) )  +         &
1281          0.5_wp * ABS( w_av(k_wind, j, i) + w_av(k_wind+1, j, i) ) )
1282    ELSE
1283!
1284!-- Calculate ta from Tp assuming dry adiabatic laps rate
1285       ta = pt(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
1286
1287       vp = bio_fill_value
1288       IF ( humidity ) THEN
1289          vp = q(k, j, i)
1290       ENDIF
1291
1292       ws = ( 0.5_wp * ABS( u(k_wind, j, i) + u(k_wind, j, i+1) )  +           &
1293          0.5_wp * ABS( v(k_wind, j, i) + v(k_wind, j+1, i) )  +               &
1294          0.5_wp * ABS( w(k_wind, j, i) + w(k_wind+1, j, i) ) )
1295
1296    ENDIF
1297!
1298!-- Local air pressure
1299    pair = surface_pressure
1300!
1301!-- Calculate water vapour pressure at saturation and convert to hPa
1302!-- The magnus formula is limited to temperatures up to 333.15 K to
1303!   avoid negative values of vp_sat
1304    IF ( vp > -998._wp ) THEN
1305       vp_sat = 0.01_wp * magnus( MIN( ta + degc_to_k, 333.15_wp ) )
1306       vp  = vp * pair / ( vp + 0.622_wp )
1307       IF ( vp > vp_sat ) vp = vp_sat
1308    ENDIF
1309!
1310!-- local mtr value at [i,j]
1311    tmrt = bio_fill_value  !< this can be a valid result (e.g. for inside some ostacle)
1312    IF ( radiation ) THEN
1313!
1314!--    Use MRT from RTM precalculated in tmrt_grid
1315       tmrt = tmrt_grid(j,i)
1316    ENDIF
1317
1318 END SUBROUTINE bio_get_thermal_index_input_ij
1319
1320
1321!------------------------------------------------------------------------------!
1322! Description:
1323! ------------
1324!> Calculate static thermal indices for any point within a 2D grid
1325!> time_integration.f90: 1065ff
1326!------------------------------------------------------------------------------!
1327 SUBROUTINE bio_calculate_thermal_index_maps ( av )
1328
1329    IMPLICIT NONE
1330!
1331!-- Input attributes
1332    LOGICAL, INTENT ( IN ) ::  av  !< Calculate based on averaged input conditions?
1333!
1334!-- Internal variables
1335    INTEGER(iwp) ::  i, j     !< Running index
1336
1337    REAL(wp) ::  clo          !< Clothing index                (no dimension)
1338    REAL(wp) ::  ta           !< Air temperature                  (degree_C)
1339    REAL(wp) ::  vp           !< Vapour pressure                  (hPa)
1340    REAL(wp) ::  ws           !< Wind speed    (local level)      (m/s)
1341    REAL(wp) ::  pair         !< Air pressure                     (hPa)
1342    REAL(wp) ::  perct_ij     !< Perceived temperature            (degree_C)
1343    REAL(wp) ::  utci_ij      !< Universal thermal climate index  (degree_C)
1344    REAL(wp) ::  pet_ij       !< Physiologically equivalent temperature  (degree_C)
1345    REAL(wp) ::  tmrt_ij      !< Mean radiant temperature         (degree_C)
1346
1347!
1348!-- fill out the MRT 2D grid from appropriate source (RTM, RRTMG,...)
1349    CALL bio_calculate_mrt_grid ( av )
1350
1351    DO i = nxl, nxr
1352       DO j = nys, nyn
1353!
1354!--       Determine local input conditions
1355          tmrt_ij = bio_fill_value
1356          vp      = bio_fill_value
1357!
1358!--       Determine input only if
1359          CALL bio_get_thermal_index_input_ij ( av, i, j, ta, vp,           &
1360                                                ws, pair, tmrt_ij )
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.