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

Last change on this file since 3582 was 3582, checked in by suehring, 3 years ago

Merge branch salsa with trunk

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