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

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

Changes related to clean-up of biometeorology (by dom_dwd_user)

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