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

Last change on this file since 3685 was 3685, checked in by knoop, 2 years ago

Some interface calls moved to module_interface + cleanup

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