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

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

Implementation of human thermal indices (from branch biomet_p2 at r3444)

  • Property svn:keywords set to Id
File size: 49.6 KB
RevLine 
[3448]1!> @file biometeorology_mod.f90
2!--------------------------------------------------------------------------------!
[3321]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!
[3448]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!--------------------------------------------------------------------------------!
[3321]21!
[3337]22! Current revisions:
[3321]23! -----------------
[3337]24!
25!
26! Former revisions:
[3321]27! -----------------
28! $Id: biometeorology_mod.f90 3448 2018-10-29 18:14:31Z kanani $
[3337]29! Initial revision
[3321]30!
[3337]31!
32!
[3321]33! Authors:
34! --------
[3448]35! @author Dominik Froehlich <dominik.froehlich@dwd.de>
[3321]36! @author Jaroslav Resler <resler@cs.cas.cz>
37!
[3448]38!
39! Description:
40! ------------
41!> Human thermal comfort module calculating thermal perception of a sample
42!> human being under the current meteorological conditions.
43!>
44!> @todo Alphabetical sorting of "USE ..." lists, "ONLY" list, variable declarations
45!>       (per subroutine: first all CHARACTERs, then INTEGERs, LOGICALs, REALs, )
46!> @todo Comments start with capital letter --> "!-- Include..."
47!> @todo Variable and routine names strictly in lowercase letters and in English
48!>
49!> @note nothing now
50!>
51!> @bug  no known bugs by now
[3321]52!------------------------------------------------------------------------------!
[3448]53 MODULE biometeorology_mod
[3321]54
55    USE arrays_3d,                                                             &
[3448]56       ONLY:  pt, p, u, v, w, q
[3321]57
[3448]58    USE averaging,                                                             &
59       ONLY:  pt_av, q_av, u_av, v_av, w_av
60
[3361]61    USE basic_constants_and_equations_mod,                                     &
[3448]62       ONLY:  magnus
[3361]63
[3448]64    USE biometeorology_ipt_mod
[3321]65
[3448]66    USE biometeorology_pet_mod
[3321]67
[3448]68    USE biometeorology_pt_mod,                                                 &
69       ONLY: calculate_pt_static
[3321]70
[3448]71    USE biometeorology_utci_mod
[3321]72
[3448]73    USE control_parameters,                                                    &
74       ONLY:  average_count_3d, biometeorology, dz, dz_stretch_factor,         &
75              dz_stretch_level, humidity, initializing_actions, nz_do3d,       &
76              simulated_time, surface_pressure
[3321]77
[3448]78    USE grid_variables,                                                        &
79       ONLY:  ddx, dx, ddy, dy
[3321]80
[3448]81    USE indices,                                                               &
82       ONLY:  nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr
[3321]83
[3448]84    USE kinds  !< Set precision of INTEGER and REAL arrays according to PALM
[3321]85
[3448]86!-- Import radiation model to obtain input for mean radiant temperature
87    USE radiation_model_mod,                                                   &
88       ONLY: ix, iy, iz, id, mrt_nlevels, mrt_include_sw,                      &
89             mrtinsw, mrtinlw, mrtbl, nmrtbl, radiation, rad_sw_in,            &
90             rad_sw_out, rad_lw_in, rad_lw_out
[3321]91
[3448]92    USE surface_mod,                                                           &
93       ONLY: get_topography_top_index_ji
[3321]94
[3448]95    IMPLICIT NONE
[3321]96
[3448]97    PRIVATE
[3321]98
[3448]99!-- Declare all global variables within the module (alphabetical order)
100    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tmrt_grid  !< tmrt results (°C)
101    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pt_grid    !< PT results   (°C)
102    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci_grid  !< UTCI results (°C)
103    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet_grid   !< PET results  (°C)
104!-- Grids for averaged thermal indices       
105    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pt_av_grid    !< PT results (aver. input)   (°C)
106    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci_av_grid  !< UTCI results (aver. input) (°C)
107    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet_av_grid   !< PET results (aver. input)  (°C)
108!-- Grids for radiation_model
109    REAL(wp), DIMENSION(:), ALLOCATABLE ::  biom_mrt     !< biometeorology mean radiant temperature for each MRT box
110    REAL(wp), DIMENSION(:), ALLOCATABLE ::  biom_mrt_av  !< time average mean
[3321]111
[3448]112    INTEGER( iwp ) ::  biom_cell_level     !< cell level biom calculates for
113    REAL ( wp )    ::  biom_output_height  !< height output is calculated in m
114    REAL ( wp )    ::  time_biom_results   !< the time, the last set of biom results have been calculated for
115    REAL ( wp ), PARAMETER ::  cels_offs = 273.15_wp  !< Kelvin-Celsius offset (K)
116    REAL ( wp ), PARAMETER ::  sigma_sb  = 5.67037321E-8_wp  !< Stefan-Boltzmann constant
117    REAL ( wp ), PARAMETER ::  human_absorb = 0.7_wp  !< SW absorbtivity of a human body (Fanger 1972)
118    REAL ( wp ), PARAMETER ::  human_emiss = 0.97_wp  !< LW emissivity of a human body after (Fanger 1972)
119!--
[3321]120
[3448]121    LOGICAL ::  aver_pt  = .FALSE.  !< switch: do pt averaging in this module? (if .FALSE. this is done globally)
122    LOGICAL ::  aver_q   = .FALSE.  !< switch: do e  averaging in this module?
123    LOGICAL ::  aver_u   = .FALSE.  !< switch: do u  averaging in this module?
124    LOGICAL ::  aver_v   = .FALSE.  !< switch: do v  averaging in this module?
125    LOGICAL ::  aver_w   = .FALSE.  !< switch: do w  averaging in this module?
126    LOGICAL ::  average_trigger_pt   = .FALSE.  !< update averaged input on call to biom_pt?
127    LOGICAL ::  average_trigger_utci = .FALSE.  !< update averaged input on call to biom_utci?
128    LOGICAL ::  average_trigger_pet  = .FALSE.  !< update averaged input on call to biom_pet?
[3321]129
[3448]130    LOGICAL ::  biom_pt        = .TRUE.   !< Turn index PT (instant. input) on or off
131    LOGICAL ::  biom_pt_av     = .TRUE.   !< Turn index PT (averaged input) on or off
132    LOGICAL ::  biom_pet       = .TRUE.   !< Turn index PET (instant. input) on or off
133    LOGICAL ::  biom_pet_av    = .TRUE.   !< Turn index PET (averaged input) on or off
134    LOGICAL ::  biom_utci      = .TRUE.   !< Turn index UTCI (instant. input) on or off
135    LOGICAL ::  biom_utci_av   = .TRUE.   !< Turn index UTCI (averaged input) on or off
[3321]136
[3448]137!-- Add INTERFACES that must be available to other modules (alphabetical order)
[3321]138
[3448]139    PUBLIC biom_3d_data_averaging, biom_check_data_output,                     &
140    biom_calculate_static_grid, biom_calc_ipt,                                 &
141    biom_check_parameters, biom_data_output_3d, biom_data_output_2d,           &
142    biom_define_netcdf_grid, biom_determine_input_at, biom_header, biom_init,  &
143    biom_init_arrays, biom_parin, biom_pt, biom_pt_av, biom_pet, biom_pet_av,  &
144    biom_utci, biom_utci_av, time_biom_results
[3321]145
[3448]146!
147!-- PALM interfaces:
148!
149!-- 3D averaging for HTCM _INPUT_ variables
150    INTERFACE biom_3d_data_averaging
151       MODULE PROCEDURE biom_3d_data_averaging
152    END INTERFACE biom_3d_data_averaging
[3321]153
[3448]154!-- Calculate static thermal indices PT, UTCI and/or PET
155    INTERFACE biom_calculate_static_grid
156       MODULE PROCEDURE biom_calculate_static_grid
157    END INTERFACE biom_calculate_static_grid
[3321]158
[3448]159!-- Calculate the dynamic index iPT (to be caled by the agent model)
160    INTERFACE biom_calc_ipt
161       MODULE PROCEDURE biom_calc_ipt
162    END INTERFACE biom_calc_ipt
[3321]163
[3448]164!-- Data output checks for 2D/3D data to be done in check_parameters
165     INTERFACE biom_check_data_output
166        MODULE PROCEDURE biom_check_data_output
167     END INTERFACE biom_check_data_output
[3321]168
[3448]169!-- Input parameter checks to be done in check_parameters
170    INTERFACE biom_check_parameters
171       MODULE PROCEDURE biom_check_parameters
172    END INTERFACE biom_check_parameters
[3321]173
[3448]174!-- Data output of 2D quantities
175    INTERFACE biom_data_output_2d
176       MODULE PROCEDURE biom_data_output_2d
177    END INTERFACE biom_data_output_2d
[3321]178
[3448]179!-- no 3D data, thus, no averaging of 3D data, removed
180    INTERFACE biom_data_output_3d
181       MODULE PROCEDURE biom_data_output_3d
182    END INTERFACE biom_data_output_3d
[3321]183
[3448]184!-- Definition of data output quantities
185    INTERFACE biom_define_netcdf_grid
186       MODULE PROCEDURE biom_define_netcdf_grid
187    END INTERFACE biom_define_netcdf_grid
[3321]188
[3448]189!-- Obtains all relevant input values to estimate local thermal comfort/stress
190    INTERFACE biom_determine_input_at
191       MODULE PROCEDURE biom_determine_input_at
192    END INTERFACE biom_determine_input_at
[3321]193
[3448]194!-- Output of information to the header file
195    INTERFACE biom_header
196       MODULE PROCEDURE biom_header
197    END INTERFACE biom_header
[3321]198
[3448]199!-- Initialization actions
200    INTERFACE biom_init
201       MODULE PROCEDURE biom_init
202    END INTERFACE biom_init
[3321]203
[3448]204!-- Initialization of arrays
205    INTERFACE biom_init_arrays
206       MODULE PROCEDURE biom_init_arrays
207    END INTERFACE biom_init_arrays
[3321]208
[3448]209!-- Reading of NAMELIST parameters
210    INTERFACE biom_parin
211       MODULE PROCEDURE biom_parin
212    END INTERFACE biom_parin
[3321]213
[3448]214 
215 CONTAINS
216 
217 
[3321]218!------------------------------------------------------------------------------!
219! Description:
220! ------------
[3448]221!> Sum up and time-average biom input quantities as well as allocate
222!> the array necessary for storing the average.
[3321]223!------------------------------------------------------------------------------!
[3448]224 SUBROUTINE biom_3d_data_averaging( mode, variable )
[3321]225
226    IMPLICIT NONE
227
[3448]228    CHARACTER (LEN=*) ::  mode     !<
229    CHARACTER (LEN=*) ::  variable !<
[3321]230
[3448]231    INTEGER(iwp) ::  i  !<
232    INTEGER(iwp) ::  j  !<
233    INTEGER(iwp) ::  k  !<
[3321]234
235
236    IF ( mode == 'allocate' )  THEN
237
238       SELECT CASE ( TRIM( variable ) )
[3448]239
240          CASE ( 'biom_mrt' )
241                IF ( .NOT. ALLOCATED( biom_mrt_av ) )  THEN
242                   ALLOCATE( biom_mrt_av(nmrtbl) )
[3321]243                ENDIF
[3448]244                biom_mrt_av = 0.0_wp
[3321]245
[3448]246          CASE ( 'biom_pt', 'biom_utci', 'biom_pet' )
247
248!--          Indices in unknown order as depending on input file, determine
249!            first index to average und update only once
250             IF ( .NOT. average_trigger_pt .AND. .NOT. average_trigger_utci    &
251                .AND. .NOT. average_trigger_pet ) THEN
252                IF ( TRIM( variable ) == 'biom_pt' ) THEN
253                    average_trigger_pt = .TRUE.
[3321]254                ENDIF
[3448]255                IF ( TRIM( variable ) == 'biom_utci' ) THEN
256                    average_trigger_utci = .TRUE.
257                ENDIF
258                IF ( TRIM( variable ) == 'biom_pet' ) THEN
259                    average_trigger_pet = .TRUE.
260                ENDIF
261             ENDIF
[3321]262
[3448]263!--          Only continue if updateindex
264             IF ( average_trigger_pt   .AND. TRIM( variable ) /= 'biom_pt')   RETURN
265             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'biom_utci') RETURN
266             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'biom_pet')  RETURN
267
268!--          Set averaging switch to .TRUE. if not done by other module before!
269             IF ( .NOT. ALLOCATED( pt_av ) )  THEN
270                ALLOCATE( pt_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
271                aver_pt = .TRUE.
272             ENDIF
273             pt_av = 0.0_wp
274
275             IF ( .NOT. ALLOCATED( q_av ) )  THEN
276                ALLOCATE( q_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
277                aver_q = .TRUE.
278             ENDIF
279             q_av = 0.0_wp
280
281             IF ( .NOT. ALLOCATED( u_av ) )  THEN
282                ALLOCATE( u_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
283                aver_u = .TRUE.
284             ENDIF
285             u_av = 0.0_wp
286
287             IF ( .NOT. ALLOCATED( v_av ) )  THEN
288                ALLOCATE( v_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
289                aver_v = .TRUE.
290             ENDIF
291             v_av = 0.0_wp
292
293             IF ( .NOT. ALLOCATED( w_av ) )  THEN
294                ALLOCATE( w_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
295                aver_w = .TRUE.
296             ENDIF
297             w_av = 0.0_wp
298
[3321]299          CASE DEFAULT
300             CONTINUE
301
302       END SELECT
303
304    ELSEIF ( mode == 'sum' )  THEN
305
306       SELECT CASE ( TRIM( variable ) )
307
[3448]308          CASE ( 'biom_mrt' )
309             IF ( ALLOCATED( biom_mrt_av ) )  THEN
[3321]310
311                 IF ( nmrtbl > 0 )  THEN
312                    IF ( mrt_include_sw )  THEN
[3448]313                       biom_mrt_av(:) = biom_mrt_av(:) + &
314                          ((human_absorb*mrtinsw(:) + human_emiss*mrtinlw(:))  &
315                          / (human_emiss*sigma_sb)) ** .25_wp - cels_offs
[3321]316                    ELSE
[3448]317                       biom_mrt_av(:) = biom_mrt_av(:) + &
318                          (human_emiss * mrtinlw(:) / sigma_sb) ** .25_wp      &
319                          - cels_offs
[3321]320                    ENDIF
321                 ENDIF
322             ENDIF
323
[3448]324          CASE ( 'biom_pt', 'biom_utci', 'biom_pet' )
[3321]325
[3448]326!--          Only continue if updateindex
327             IF ( average_trigger_pt   .AND. TRIM( variable ) /= 'biom_pt')    &
328                RETURN
329             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'biom_utci')  &
330                RETURN
331             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'biom_pet')   &
332                RETURN
[3321]333
[3448]334             IF ( ALLOCATED( pt_av ) .AND. aver_pt ) THEN
335                DO  i = nxl, nxr
336                   DO  j = nys, nyn
337                      DO  k = nzb, nzt+1
338                         pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i)
339                      ENDDO
340                   ENDDO
341                ENDDO
342             ENDIF
[3321]343
[3448]344             IF ( ALLOCATED( q_av )  .AND. aver_q ) THEN
345                DO  i = nxl, nxr
346                   DO  j = nys, nyn
347                      DO  k = nzb, nzt+1
348                         q_av(k,j,i) = q_av(k,j,i) + q(k,j,i)
349                      ENDDO
350                   ENDDO
[3321]351                ENDDO
352             ENDIF
353
[3448]354             IF ( ALLOCATED( u_av )  .AND. aver_u ) THEN
355                DO  i = nxl, nxr
356                   DO  j = nys, nyn
357                      DO  k = nzb, nzt+1
358                         u_av(k,j,i) = u_av(k,j,i) + u(k,j,i)
359                      ENDDO
360                   ENDDO
361                ENDDO
362             ENDIF
363
364             IF ( ALLOCATED( v_av )  .AND. aver_v ) THEN
365                DO  i = nxl, nxr
366                   DO  j = nys, nyn
367                      DO  k = nzb, nzt+1
368                         v_av(k,j,i) = v_av(k,j,i) + v(k,j,i)
369                      ENDDO
370                   ENDDO
371                ENDDO
372             ENDIF
373
374             IF ( ALLOCATED( w_av )  .AND. aver_w ) THEN
375                DO  i = nxl, nxr
376                   DO  j = nys, nyn
377                      DO  k = nzb, nzt+1
378                         w_av(k,j,i) = w_av(k,j,i) + w(k,j,i)
379                      ENDDO
380                   ENDDO
381                ENDDO
382             ENDIF
383
384           CASE DEFAULT
[3321]385             CONTINUE
386
387       END SELECT
388
389    ELSEIF ( mode == 'average' )  THEN
390
391       SELECT CASE ( TRIM( variable ) )
392
[3448]393          CASE ( 'biom_mrt' )
394             IF ( ALLOCATED( biom_mrt_av ) )  THEN
395                biom_mrt_av(:) = biom_mrt_av(:) / REAL( average_count_3d, KIND=wp )
[3321]396             ENDIF
397
[3448]398          CASE ( 'biom_pt', 'biom_utci', 'biom_pet' )
399
400!--          Only continue if update index
401             IF ( average_trigger_pt   .AND. TRIM( variable ) /= 'biom_pt')    &
402                RETURN
403             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'biom_utci')  &
404                RETURN
405             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'biom_pet')   &
406                RETURN
407
408             IF ( ALLOCATED( pt_av ) .AND. aver_pt ) THEN
409                DO  i = nxl, nxr
410                   DO  j = nys, nyn
411                      DO  k = nzb, nzt+1
412                         pt_av(k,j,i) = pt_av(k,j,i) / REAL( average_count_3d, KIND=wp )
413                      ENDDO
414                   ENDDO
415                ENDDO
[3321]416             ENDIF
417
[3448]418             IF ( ALLOCATED( q_av ) .AND. aver_q ) THEN
419                DO  i = nxl, nxr
420                   DO  j = nys, nyn
421                      DO  k = nzb, nzt+1
422                         q_av(k,j,i) = q_av(k,j,i) / REAL( average_count_3d, KIND=wp )
423                      ENDDO
424                   ENDDO
425                ENDDO
426             ENDIF
[3321]427
[3448]428             IF ( ALLOCATED( u_av ) .AND. aver_u ) THEN
429                DO  i = nxl, nxr
430                   DO  j = nys, nyn
431                      DO  k = nzb, nzt+1
432                         u_av(k,j,i) = u_av(k,j,i) / REAL( average_count_3d, KIND=wp )
433                      ENDDO
434                   ENDDO
435                ENDDO
436             ENDIF
437
438             IF ( ALLOCATED( v_av ) .AND. aver_v ) THEN
439                DO  i = nxl, nxr
440                   DO  j = nys, nyn
441                      DO  k = nzb, nzt+1
442                         v_av(k,j,i) = v_av(k,j,i) / REAL( average_count_3d, KIND=wp )
443                      ENDDO
444                   ENDDO
445                ENDDO
446             ENDIF
447
448             IF ( ALLOCATED( w_av ) .AND. aver_w ) THEN
449                DO  i = nxl, nxr
450                   DO  j = nys, nyn
451                      DO  k = nzb, nzt+1
452                         w_av(k,j,i) = w_av(k,j,i) / REAL( average_count_3d, KIND=wp )
453                      ENDDO
454                   ENDDO
455                ENDDO
456             ENDIF
457
458!--          Udate thermal indices with derived averages
459             CALL biom_calculate_static_grid ( .TRUE. )
460
461        END SELECT
462
[3321]463    ENDIF
464
465
[3448]466 END SUBROUTINE biom_3d_data_averaging
[3321]467
[3448]468
469
[3321]470!------------------------------------------------------------------------------!
471! Description:
472! ------------
[3448]473!> Check data output for biometeorology model
[3321]474!------------------------------------------------------------------------------!
[3448]475 SUBROUTINE biom_check_data_output( var, unit )
[3321]476
[3448]477       USE control_parameters,                                                 &
478           ONLY: data_output, message_string
[3321]479
[3448]480       IMPLICIT NONE
[3321]481
[3448]482       CHARACTER (LEN=*) ::  unit     !< The unit for the variable var
483       CHARACTER (LEN=*) ::  var      !< The variable in question
[3321]484
485
[3448]486                unit = '°C'
487                IF ( .NOT. biometeorology ) THEN
488                  message_string = 'output of "' // TRIM( var ) // '" req' //  &
489                        'uires biometeorology = .TRUE. in initialisati' &
490                        // 'on_parameters'
491                  CALL message( 'check_parameters', 'PA0114', 1, 2, 0, 6, 0 )
492                  unit = 'illegal'
493                ENDIF
494                IF ( .NOT.  radiation ) THEN
495                   message_string = 'output of "' // TRIM( var ) // '" require'&
496                                    // 's radiation = .TRUE.'
497                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
498                   unit = 'illegal'
499                ENDIF
500                IF ( mrt_nlevels == 0 ) THEN
501                   message_string = 'output of "' // TRIM( var ) // '" require'&
502                                    // 's mrt_nlevels > 0'
503                   CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
504                   unit = 'illegal'
505                ENDIF
[3321]506
[3448]507 END SUBROUTINE biom_check_data_output
[3321]508
[3448]509!------------------------------------------------------------------------------!
510! Description:
511! ------------
512!> Check parameters routine for biom module
513!> check_parameters 1302
514!------------------------------------------------------------------------------!
515 SUBROUTINE biom_check_parameters
[3321]516
[3448]517    USE control_parameters,                                                    &
518        ONLY:  message_string
[3321]519
520
[3448]521    IMPLICIT NONE
[3321]522
[3448]523
524!--    Disabled as long as radiation model not available
525       IF ( .NOT. radiation )  THEN
526          message_string = 'The human thermal comfort module does require ' // &
527                           'radiation information in terms of the mean ' //    &
528                           'radiant temperature, but radiation is not enabled!'
529          CALL message( 'check_parameters', 'PAHU01', 0, 0, 0, 6, 0 )
530       ENDIF
531
532       IF ( .NOT. humidity )  THEN
533          message_string = 'The human thermal comfort module does require ' // &
534                           'air humidity information, but humidity module ' // &
535                           'is disabled!'
536          CALL message( 'check_parameters', 'PAHU02', 0, 0, 0, 6, 0 )
537       ENDIF
538
539
540 END SUBROUTINE biom_check_parameters
541
542
[3321]543!------------------------------------------------------------------------------!
544!
545! Description:
546! ------------
[3448]547!> Subroutine defining 3D output variables (dummy, always 2d!)
548!> data_output_3d 709ff
[3321]549!------------------------------------------------------------------------------!
[3448]550 SUBROUTINE biom_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
[3321]551
552    USE indices
553
554    USE kinds
555
556
557    IMPLICIT NONE
558
[3448]559!-- Input variables
560    CHARACTER (LEN=*), INTENT(IN) ::  variable   !< Char identifier to select var for output
561    INTEGER(iwp), INTENT(IN) ::  av       !< Use averaged data? 0 = no, 1 = yes?
562    INTEGER(iwp), INTENT(IN) ::  nzb_do   !< Unused. 2D. nz bottom to nz top
563    INTEGER(iwp), INTENT(IN) ::  nzt_do   !< Unused.
[3321]564
[3448]565!-- Output variables
566    LOGICAL, INTENT(OUT) ::  found   !< Output found?
567    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< Temp. result grid to return
[3321]568
[3448]569!-- Internal variables
570    INTEGER(iwp) ::  l    !< Running index, radiation grid
571    INTEGER(iwp) ::  i    !< Running index, x-dir
572    INTEGER(iwp) ::  j    !< Running index, y-dir
573    INTEGER(iwp) ::  k    !< Running index, z-dir
[3321]574
[3448]575    CHARACTER (LEN=:), allocatable ::  variable_short  !< Trimmed version of char variable
[3321]576
[3448]577    REAL(wp), PARAMETER ::  fill_value = -999._wp
578    REAL(wp) ::  mrt  !< Buffer for mean radiant temperature
[3321]579
580    found = .TRUE.
[3448]581    variable_short = TRIM( variable )
[3321]582
[3448]583    IF ( variable_short(1:4) /= 'biom' ) THEN
584!--   Nothing to do, set found to FALSE and return immediatelly
585      found = .FALSE.
586      RETURN
587    ENDIF
[3321]588
[3448]589    SELECT CASE ( variable_short )
[3321]590
[3448]591        CASE ( 'biom_mrt' )
[3321]592
593            local_pf = REAL( fill_value, KIND = wp )
594            DO  l = 1, nmrtbl
595                i = mrtbl(ix,l)
596                j = mrtbl(iy,l)
597                k = mrtbl(iz,l)
598                IF ( mrt_include_sw )  THEN
599                    mrt = ((human_absorb*mrtinsw(l) + human_emiss*mrtinlw(l)) &
600                           / (human_emiss*sigma_sb)) ** .25_wp
601                ELSE
602                    mrt = (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp
603                ENDIF
604                local_pf(i,j,k) = mrt
605            ENDDO
606
[3448]607        CASE ( 'biom_tmrt' )        ! 2d-array
608              DO  i = nxl, nxr
609                 DO  j = nys, nyn
610                    local_pf(i,j,nzb_do) = REAL( tmrt_grid(j,i), sp )
611                 ENDDO
612              ENDDO
613
614        CASE ( 'biom_pt' )        ! 2d-array
[3321]615            IF ( av == 0 )  THEN
[3448]616              DO  i = nxl, nxr
617                 DO  j = nys, nyn
618                    local_pf(i,j,nzb_do) = REAL( pt_grid(j,i), sp )
619                 ENDDO
620              ENDDO
621            ELSE
622              DO  i = nxl, nxr
623                 DO  j = nys, nyn
624                    local_pf(i,j,nzb_do) = REAL( pt_av_grid(j,i), sp )
625                 ENDDO
626              ENDDO
627            END IF
[3321]628
[3448]629        CASE ( 'biom_utci' )        ! 2d-array
630            IF ( av == 0 )  THEN
631              DO  i = nxl, nxr
632                 DO  j = nys, nyn
633                    local_pf(i,j,nzb_do) = REAL( utci_grid(j,i), sp )
634                 ENDDO
635              ENDDO
636            ELSE
637              DO  i = nxl, nxr
638                 DO  j = nys, nyn
639                    local_pf(i,j,nzb_do) = REAL( utci_av_grid(j,i), sp )
640                 ENDDO
641              ENDDO
642            END IF
[3321]643
[3448]644        CASE ( 'biom_pet' )        ! 2d-array
645            IF ( av == 0 )  THEN
646              DO  i = nxl, nxr
647                 DO  j = nys, nyn
648                    local_pf(i,j,nzb_do) = REAL( pet_grid(j,i), sp )
649                 ENDDO
650              ENDDO
651            ELSE
652              DO  i = nxl, nxr
653                 DO  j = nys, nyn
654                    local_pf(i,j,nzb_do) = REAL( pet_av_grid(j,i), sp )
655                 ENDDO
656              ENDDO
657            END IF
[3321]658
[3448]659       CASE DEFAULT
660          found = .FALSE.
[3321]661
[3448]662    END SELECT
[3321]663
[3448]664 END SUBROUTINE biom_data_output_3d
665
666!------------------------------------------------------------------------------!
667!
668! Description:
669! ------------
670!> Subroutine defining 2D output variables
671!> data_output_2d 1188ff
672!------------------------------------------------------------------------------!
673 SUBROUTINE biom_data_output_2d( av, variable, found, grid, local_pf,          &
674                                      two_d, nzb_do, nzt_do, fill_value )
675
676    USE indices,                                                               &
677       ONLY: nxl, nxl, nxr, nxr, nyn, nyn, nys, nys, nzb, nzt
678
679    USE kinds
680
681
682    IMPLICIT NONE
683
684!-- Input variables
685    CHARACTER (LEN=*), INTENT(IN) ::  variable    !< Char identifier to select var for output
686    INTEGER(iwp), INTENT(IN)      ::  av          !< Use averaged data? 0 = no, 1 = yes?
687    INTEGER(iwp), INTENT(IN)      ::  nzb_do      !< Unused. 2D. nz bottom to nz top
688    INTEGER(iwp), INTENT(IN)      ::  nzt_do      !< Unused.
689    REAL(wp), INTENT(in)          ::  fill_value  !< Fill value for unassigned locations
690
691!-- Output variables
692    CHARACTER (LEN=*), INTENT(OUT) ::  grid   !< Grid type (always "zu1" for biom)
693    LOGICAL, INTENT(OUT)           ::  found  !< Output found?
694    LOGICAL, INTENT(OUT)           ::  two_d  !< Flag parameter that indicates 2D variables, horizontal cross sections, must be .TRUE.
695    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !< Temp. result grid to return
696
697!-- Internal variables
698    CHARACTER (LEN=:), allocatable ::  variable_short  !< Trimmed version of char variable
699    INTEGER(iwp) ::  i        !< Running index, x-dir
700    INTEGER(iwp) ::  j        !< Running index, y-dir
701    INTEGER(iwp) ::  k        !< Running index, z-dir
702
703
704    found = .TRUE.
705    variable_short = TRIM( variable )
706    IF ( variable_short(1:4) == 'biom' ) THEN
707       two_d = .TRUE.
708       grid = 'zu1'
709    ELSE
710       found = .FALSE.
711       grid  = 'none'
712       RETURN
713    ENDIF
714
715    local_pf = fill_value
716
717    SELECT CASE ( variable_short )
718
719
720        CASE ( 'biom_tmrt_xy' )        ! 2d-array
721              DO  i = nxl, nxr
722                 DO  j = nys, nyn
723                    local_pf(i,j,1) = tmrt_grid(j,i)
724                 ENDDO
725              ENDDO
726
727
728        CASE ( 'biom_pt_xy' )        ! 2d-array
729            IF ( av == 0 )  THEN
730              DO  i = nxl, nxr
731                 DO  j = nys, nyn
732                    local_pf(i,j,nzb+1) = pt_grid(j,i)
733                 ENDDO
734              ENDDO
735            ELSE
736              DO  i = nxl, nxr
737                 DO  j = nys, nyn
738                    local_pf(i,j,nzb+1) = pt_av_grid(j,i)
739                 ENDDO
740              ENDDO
741            END IF
742
743
744        CASE ( 'biom_utci_xy' )        ! 2d-array
745            IF ( av == 0 )  THEN
746              DO  i = nxl, nxr
747                 DO  j = nys, nyn
748                    local_pf(i,j,nzb+1) = utci_grid(j,i)
749                 ENDDO
750              ENDDO
751            ELSE
752              DO  i = nxl, nxr
753                 DO  j = nys, nyn
754                    local_pf(i,j,nzb+1) = utci_av_grid(j,i)
755                 ENDDO
756              ENDDO
757            END IF
758
759
760        CASE ( 'biom_pet_xy' )        ! 2d-array
761            IF ( av == 0 )  THEN
762              DO  i = nxl, nxr
763                 DO  j = nys, nyn
764                    local_pf(i,j,nzb+1) = pet_grid(j,i)
765                 ENDDO
766              ENDDO
767            ELSE
768              DO  i = nxl, nxr
769                 DO  j = nys, nyn
770                    local_pf(i,j,nzb+1) = pet_av_grid(j,i)
771                 ENDDO
772              ENDDO
773            END IF
774
775
[3321]776       CASE DEFAULT
777          found = .FALSE.
[3448]778          grid  = 'none'
[3321]779
780    END SELECT
781
782
[3448]783 END SUBROUTINE biom_data_output_2d
[3321]784
785
[3448]786!------------------------------------------------------------------------------!
787! Description:
788! ------------
789!> Subroutine defining appropriate grid for netcdf variables.
790!> It is called out from subroutine netcdf_interface_mod.
791!> netcdf_interface_mod 918ff
792!------------------------------------------------------------------------------!
793 SUBROUTINE biom_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
[3321]794
[3448]795    IMPLICIT NONE
796
797!-- Input variables
798    CHARACTER (LEN=*), INTENT(IN)  ::  var      !< Name of output variable
799
800!-- Output variables
801    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x   !< x grid of output variable
802    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y   !< y grid of output variable
803    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z   !< z grid of output variable
804
805    LOGICAL, INTENT(OUT)           ::  found    !< Flag if output var is found
806
807!-- Local variables
808    LOGICAL      :: is2d  !< Var is 2d?
809
810    INTEGER(iwp) :: l     !< Length of the var array
811
812
813    found  = .FALSE.
814    grid_x = 'none'
815    grid_y = 'none'
816    grid_z = 'none'
817
818    l = MAX( 2, LEN_TRIM( var ) )
819    is2d = ( var(l-1:l) == 'xy' )
820
821
822    IF ( var(1:4) == 'biom' ) THEN
823       found  = .TRUE.
824       grid_x = 'x'
825       grid_y = 'y'
826       grid_z = 'zu'
827       IF ( is2d ) grid_z = 'zu1'
828    ENDIF
829
830 END SUBROUTINE biom_define_netcdf_grid
831
[3321]832!------------------------------------------------------------------------------!
833! Description:
834! ------------
[3448]835!> Header output for biom module
836!> header 982
837!------------------------------------------------------------------------------!
838 SUBROUTINE biom_header( io )
839
840    IMPLICIT NONE
841
842!-- Input variables
843    INTEGER(iwp), INTENT(IN) ::  io           !< Unit of the output file
844
845!-- Internal variables
846    CHARACTER (LEN=86) ::  output_height_chr  !< String for output height
847
848    WRITE( output_height_chr, '(F8.1,7X)' )  biom_output_height
[3321]849!
[3448]850!-- Write biom header
851    WRITE( io, 1 )
852    WRITE( io, 2 )  TRIM( output_height_chr )
853    WRITE( io, 3 )  TRIM( ACHAR( biom_cell_level ) )
854
8551   FORMAT (//' Human thermal comfort module information:'/                    &
856              ' ------------------------------'/)
8572   FORMAT ('    --> All indices calculated for a height of (m): ', A )
8583   FORMAT ('    --> This corresponds to cell level : ', A )
859
860 END SUBROUTINE biom_header
861
862
[3321]863!------------------------------------------------------------------------------!
[3448]864! Description:
865! ------------
866!> Initialization of the HTCM
867!> init_3d_model 1987ff
868!------------------------------------------------------------------------------!
869 SUBROUTINE biom_init
[3321]870
[3448]871    IMPLICIT NONE
[3321]872
[3448]873!-- Internal vriables
874    REAL ( wp )  :: height  !< current height in meters
[3321]875
[3448]876    INTEGER ( iwp )  :: i  !< iteration index
[3321]877
[3448]878!-- Determine cell level corresponding to 1.1 m above ground level
879!   (gravimetric center of sample human)
[3321]880
[3448]881    time_biom_results = 0.0_wp
882    biom_cell_level = 0_iwp
883    biom_output_height = 0.5_wp * dz(1)
884    height = 0.0_wp
[3321]885
[3448]886    biom_cell_level = INT ( 1.099_wp / dz(1) )
887    biom_output_height = biom_output_height + biom_cell_level * dz(1)
[3321]888
[3448]889 END SUBROUTINE biom_init
[3321]890
891
[3448]892!------------------------------------------------------------------------------!
893! Description:
894! ------------
895!> Allocate biom arrays and define pointers if required
896!> init_3d_model 1047ff
897!------------------------------------------------------------------------------!
898 SUBROUTINE biom_init_arrays
[3321]899
[3448]900    IMPLICIT NONE
[3321]901
[3448]902!-- Allocate a temporary array with the desired output dimensions.
903!   Initialization omitted for performance, will be overwritten anyway
904    IF ( .NOT. ALLOCATED( tmrt_grid ) ) THEN
905      ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) )
906    ENDIF
[3321]907
[3448]908    IF ( biom_pt ) THEN
909      IF ( .NOT. ALLOCATED( pt_grid ) ) THEN
910        ALLOCATE( pt_grid (nys:nyn,nxl:nxr) )
911      ENDIF
912    ENDIF
[3321]913
[3448]914    IF ( biom_utci ) THEN
915      IF ( .NOT. ALLOCATED( utci_grid ) ) THEN
916        ALLOCATE( utci_grid (nys:nyn,nxl:nxr) )
917      ENDIF
918    ENDIF
[3321]919
[3448]920    IF ( biom_pet ) THEN
921      IF ( .NOT. ALLOCATED( pet_grid ) ) THEN
922        ALLOCATE( pet_grid (nys:nyn,nxl:nxr) )
923      ENDIF
924    END IF
[3321]925
[3448]926    IF ( biom_pt_av ) THEN
927      IF ( .NOT. ALLOCATED( pt_av_grid ) ) THEN
928        ALLOCATE( pt_av_grid (nys:nyn,nxl:nxr) )
929      ENDIF
930    ENDIF
[3321]931
[3448]932    IF ( biom_utci_av ) THEN
933      IF ( .NOT. ALLOCATED( utci_av_grid ) ) THEN
934        ALLOCATE( utci_av_grid (nys:nyn,nxl:nxr) )
935      ENDIF
936    ENDIF
[3321]937
[3448]938    IF ( biom_pet_av ) THEN
939      IF ( .NOT. ALLOCATED( pet_av_grid ) ) THEN
940        ALLOCATE( pet_av_grid (nys:nyn,nxl:nxr) )
941      ENDIF
[3321]942
[3448]943    END IF
[3321]944
[3448]945 END SUBROUTINE biom_init_arrays
946
947
[3321]948!------------------------------------------------------------------------------!
949! Description:
950! ------------
[3448]951!> Parin for &biometeorology_parameters for reading biomet parameters
[3321]952!------------------------------------------------------------------------------!
[3448]953 SUBROUTINE biom_parin
[3321]954
[3448]955    IMPLICIT NONE
[3321]956
[3448]957!
958!-- Internal variables
959    CHARACTER (LEN=80) ::  line  !< Dummy string for current line in parameter file
[3321]960
[3448]961    NAMELIST /biometeorology_parameters/  biom_pet,                            &
962                                          biom_pet_av,                         &
963                                          biom_pt,                             &
964                                          biom_pt_av,                          &
965                                          biom_utci,                           &
966                                          biom_utci_av
[3321]967
968
[3448]969!-- Try to find biometeorology_parameters namelist
970    REWIND ( 11 )
971    line = ' '
972    DO WHILE ( INDEX( line, '&biometeorology_parameters' ) == 0 )
973       READ ( 11, '(A)', END = 20 )  line
974    ENDDO
975    BACKSPACE ( 11 )
[3321]976
[3448]977!
978!-- Read biometeorology_parameters namelist
979    READ ( 11, biometeorology_parameters, ERR = 10, END = 20 )
[3321]980
[3448]981!
982!-- Set flag that indicates that the biomet_module is switched on
983    biometeorology = .TRUE.
[3321]984
[3448]985    GOTO 20
[3321]986
[3448]987!
988!-- In case of error
989 10 BACKSPACE( 11 )
990    READ( 11 , '(A)') line
991    CALL parin_fail_message( 'biometeorology_parameters', line )
[3321]992
[3448]993!
994!-- Complete
995 20 CONTINUE
[3321]996
997
[3448]998 END SUBROUTINE biom_parin
[3321]999
1000!------------------------------------------------------------------------------!
1001! Description:
1002! ------------
[3448]1003!> Calculates the mean radiant temperature (tmrt) based on the Six-directions
1004!> method according to VDI 3787 2.
[3321]1005!------------------------------------------------------------------------------!
[3448]1006 SUBROUTINE calculate_tmrt_6_directions( SW_N, SW_E, SW_S, SW_W,               &
1007    SW_U, SW_D, LW_N, LW_E, LW_S, LW_W, LW_U, LW_D, tmrt )
[3321]1008
[3448]1009    IMPLICIT NONE
[3321]1010
[3448]1011!-- Type of input of the argument list
1012!   Short- (SW_) and longwave (LW_) radiation fluxes from the six directions
1013!   North (N), East (E), South (S), West (W), up (U) and down (D)
1014    REAL(wp), INTENT ( IN )  ::  SW_N  !< Sw radflux density from N    (W/m²)
1015    REAL(wp), INTENT ( IN )  ::  SW_E  !< Sw radflux density from E    (W/m²)
1016    REAL(wp), INTENT ( IN )  ::  SW_S  !< Sw radflux density from S    (W/m²)
1017    REAL(wp), INTENT ( IN )  ::  SW_W  !< Sw radflux density from W    (W/m²)
1018    REAL(wp), INTENT ( IN )  ::  SW_U  !< Sw radflux density from U    (W/m²)
1019    REAL(wp), INTENT ( IN )  ::  SW_D  !< Sw radflux density from D    (W/m²)
1020    REAL(wp), INTENT ( IN )  ::  LW_N  !< Lw radflux density from N    (W/m²)
1021    REAL(wp), INTENT ( IN )  ::  LW_E  !< Lw radflux density from E    (W/m²)
1022    REAL(wp), INTENT ( IN )  ::  LW_S  !< Lw radflux density from S    (W/m²)
1023    REAL(wp), INTENT ( IN )  ::  LW_W  !< Lw radflux density from W    (W/m²)
1024    REAL(wp), INTENT ( IN )  ::  LW_U  !< Lw radflux density from U    (W/m²)
1025    REAL(wp), INTENT ( IN )  ::  LW_D  !< Lw radflux density from D    (W/m²)
[3321]1026
[3448]1027!-- Type of output of the argument list
1028    REAL(wp), INTENT ( OUT ) ::  tmrt  !< Mean radiant temperature (°C)
[3321]1029
[3448]1030!-- Directional weighting factors
1031    REAL(wp), PARAMETER      ::  weight_h  = 0.22_wp
1032    REAL(wp), PARAMETER      ::  weight_v  = 0.06_wp
[3321]1033
[3448]1034    REAL(wp) ::  nrfd           !<  Net radiation flux density (W/m²)
[3321]1035
[3448]1036!-- Initialise
1037    nrfd = 0._wp
[3321]1038
[3448]1039!-- Compute mean radiation flux density absorbed by rotational symetric human
1040    nrfd = ( weight_h * ( human_absorb * SW_N + human_emiss * LW_N ) ) +       &
1041       ( weight_h * ( human_absorb * SW_E + human_emiss * LW_E ) ) +           &
1042       ( weight_h * ( human_absorb * SW_S + human_emiss * LW_S ) ) +           &
1043       ( weight_h * ( human_absorb * SW_W + human_emiss * LW_W ) ) +           &
1044       ( weight_v * ( human_absorb * SW_U + human_emiss * LW_U ) ) +           &
1045       ( weight_v * ( human_absorb * SW_D + human_emiss * LW_D ) )
[3321]1046
[3448]1047!-- Compute mean radiant temperature
1048    tmrt = ( nrfd / (human_emiss * sigma_sb) )**0.25_wp - cels_offs
[3321]1049
[3448]1050 END SUBROUTINE calculate_tmrt_6_directions
[3321]1051
[3448]1052!------------------------------------------------------------------------------!
1053! Description:
1054! ------------
1055!> Very crude approximation of mean radiant temperature based on upwards and
1056!> downwards radiation fluxes only (other directions curently not available,
1057!> replace as soon as possible!)
1058!------------------------------------------------------------------------------!
1059 SUBROUTINE calculate_tmrt_2_directions( sw_u, sw_d, lw_u, lw_d, ta, tmrt )
[3321]1060
[3448]1061    IMPLICIT NONE
[3321]1062
[3448]1063!-- Type of input of the argument list
1064    REAL(wp), INTENT ( IN )  ::  sw_u  !< Shortwave radiation flux density from upper direction (W/m²)
1065    REAL(wp), INTENT ( IN )  ::  sw_d  !< Shortwave radiation flux density from lower direction (W/m²)
1066    REAL(wp), INTENT ( IN )  ::  lw_u  !< Longwave radiation flux density from upper direction (W/m²)
1067    REAL(wp), INTENT ( IN )  ::  lw_d  !< Longwave radiation flux density from lower direction (W/m²)
1068    REAL(wp), INTENT ( IN )  ::  ta    !< Air temperature (°C)
[3321]1069
[3448]1070!-- Type of output of the argument list
1071    REAL(wp), INTENT ( OUT ) ::  tmrt  !< mean radiant temperature, (°C)
[3321]1072
[3448]1073!-- Directional weighting factors and parameters
1074    REAL(wp), PARAMETER ::  weight_h  = 0.22_wp     !< Weight for horizontal radiational gain after Fanger (1972)
1075    REAL(wp), PARAMETER ::  weight_v  = 0.06_wp     !< Weight for vertical radiational gain after Fanger (1972)
[3321]1076
[3448]1077!-- Other internal variables
1078    REAL(wp) ::  sw_in
1079    REAL(wp) ::  sw_out
1080    REAL(wp) ::  lw_in
1081    REAL(wp) ::  lw_out
1082    REAL(wp) ::  nrfd     !< Net radiation flux density (W/m²)
1083    REAL(wp) ::  lw_air   !< Longwave emission by surrounding air volume (W/m²)
1084    REAL(wp) ::  sw_side  !< Shortwave immission from the sides (W/m²)
[3321]1085
[3448]1086    INTEGER(iwp) ::  no_input  !< Count missing input radiation fluxes
[3321]1087
[3448]1088!-- initialise
1089    sw_in    = sw_u
1090    sw_out   = sw_d
1091    lw_in    = lw_u
1092    lw_out   = lw_d
1093    nrfd     = 0._wp
1094    no_input = 0_iwp
[3321]1095
[3448]1096!-- test for missing input data
1097    IF ( sw_in <= -998._wp .OR. sw_out <= -998._wp .OR. lw_in <= -998._wp .OR. &
1098       lw_out <= -998._wp .OR. ta <= -998._wp ) THEN
1099       IF ( sw_in <= -998._wp ) THEN
1100          sw_in = 0.
1101          no_input = no_input + 1
1102       ENDIF
1103       IF ( sw_out <= -998._wp ) THEN
1104          sw_out = 0.
1105          no_input = no_input + 1
1106       ENDIF
1107       IF ( lw_in <= -998._wp ) THEN
1108          lw_in = 0.
1109          no_input = no_input + 1
1110       ENDIF
1111       IF ( lw_out <= -998._wp ) THEN
1112          lw_out = 0.
1113          no_input = no_input + 1
1114       ENDIF
[3321]1115
[3448]1116!-- Accept two missing radiation flux directions, fail otherwise as error might become too large
1117       IF ( ta <= -998._wp .OR. no_input >= 3 ) THEN
1118          tmrt = -999._wp
1119          RETURN
1120       ENDIF
1121    ENDIF
[3321]1122
[3448]1123    sw_side = sw_in * 0.125_wp  ! distribute half of upper sw_in to the 4 sides
1124    lw_air = ( sigma_sb * 0.95_wp * ( ta + cels_offs )**4 )
[3321]1125
[3448]1126!-- Compute mean radiation flux density absorbed by rotational symetric human
1127    nrfd = ( weight_h * ( human_absorb * sw_side + human_emiss * lw_air ) ) +         &
1128       ( weight_h * ( human_absorb * sw_side + human_emiss * lw_air ) ) +             &
1129       ( weight_h * ( human_absorb * sw_side + human_emiss * lw_air ) ) +             &
1130       ( weight_h * ( human_absorb * sw_side + human_emiss * lw_air ) ) +             &
1131       ( weight_v * ( human_absorb * (sw_in * 0.5_wp) + human_emiss * lw_in ) ) +     &
1132       ( weight_v * ( human_absorb * sw_out + human_emiss * lw_out ) )
[3321]1133
[3448]1134!-- Compute mean radiant temperature
1135    tmrt = ( nrfd / (human_emiss * sigma_sb) )**0.25_wp - cels_offs
[3321]1136
[3448]1137 END SUBROUTINE calculate_tmrt_2_directions
[3321]1138
[3448]1139!------------------------------------------------------------------------------!
1140! Description:
1141! ------------
1142!> Calculate static thermal indices for given meteorological conditions
1143!------------------------------------------------------------------------------!
1144 SUBROUTINE calculate_static_thermal_indices ( ta, vp, ws, pair, tmrt,         &
1145    pt_static, utci_static, pet_static )
[3321]1146
[3448]1147    IMPLICIT NONE
[3321]1148
[3448]1149!-- Input parameters
1150    REAL(wp), INTENT ( IN )  ::  ta           !< Air temperature                  (°C)
1151    REAL(wp), INTENT ( IN )  ::  vp           !< Vapour pressure                  (hPa)
1152    REAL(wp), INTENT ( IN )  ::  ws           !< Wind speed    (local level)      (m/s)
1153    REAL(wp), INTENT ( IN )  ::  pair         !< Air pressure                     (hPa)
1154    REAL(wp), INTENT ( IN )  ::  tmrt         !< Mean radiant temperature         (°C)
1155!-- Output parameters
1156    REAL(wp), INTENT ( OUT ) ::  pt_static    !< Perceived temperature            (°C)
1157    REAL(wp), INTENT ( OUT ) ::  utci_static  !< Universal thermal climate index  (°C)
1158    REAL(wp), INTENT ( OUT ) ::  pet_static   !< Physiologically equivalent temp. (°C)
1159!-- Temporary field, not used here
1160    REAL(wp)                 ::  clo          !< Clothing index (no dim.)
[3321]1161
[3448]1162    clo = -999._wp
[3321]1163
[3448]1164    IF ( biom_pt ) THEN
1165!-- Estimate local perceived temperature
1166       CALL calculate_pt_static( ta, vp, ws, tmrt, pair, clo, pt_static )
1167    ENDIF
[3321]1168
[3448]1169    IF ( biom_utci ) THEN
1170!-- Estimate local universal thermal climate index
1171       CALL calculate_utci_static( ta, vp, ws, tmrt, biom_output_height,       &
1172          utci_static )
1173    ENDIF
[3321]1174
[3448]1175    IF ( biom_pet ) THEN
1176!-- Estimate local physiologically equivalent temperature
1177       CALL calculate_pet_static( ta, vp, ws, tmrt, pair, pet_static )
1178    ENDIF
[3321]1179
[3448]1180 END SUBROUTINE calculate_static_thermal_indices
[3321]1181
1182
[3448]1183!------------------------------------------------------------------------------!
1184! Description:
1185! ------------
1186!> Calculate static thermal indices for 2D grid point i, j
1187!------------------------------------------------------------------------------!
1188 SUBROUTINE biom_determine_input_at( average_input, i, j, ta, vp, ws, pair,    &
1189    tmrt )
[3321]1190
[3448]1191    IMPLICIT NONE
[3321]1192
[3448]1193!-- Input variables
1194    LOGICAL,      INTENT ( IN ) ::  average_input  !< Determine averaged input conditions?
1195    INTEGER(iwp), INTENT ( IN ) ::  i     !< Running index, x-dir
1196    INTEGER(iwp), INTENT ( IN ) ::  j     !< Running index, y-dir
[3321]1197
[3448]1198!-- Output parameters
1199    REAL(wp), INTENT ( OUT )    ::  tmrt  !< Mean radiant temperature        (°C)
1200    REAL(wp), INTENT ( OUT )    ::  ta    !< Air temperature                 (°C)
1201    REAL(wp), INTENT ( OUT )    ::  vp    !< Vapour pressure                 (hPa)
1202    REAL(wp), INTENT ( OUT )    ::  ws    !< Wind speed    (local level)     (m/s)
1203    REAL(wp), INTENT ( OUT )    ::  pair  !< Air pressure                    (hPa)
[3321]1204
[3448]1205!-- Internal variables
1206    INTEGER(iwp)                ::  k     !< Running index, z-dir
1207    INTEGER(iwp)                ::  k_wind  !< Running index, z-dir, wind speed only
[3321]1208
[3448]1209    REAL(wp)                    ::  vp_sat  !< Saturation vapor pressure     (hPa)
[3321]1210
1211
[3448]1212!-- Determine cell level closest to 1.1m above ground
1213!   by making use of truncation due to int cast
1214    k = get_topography_top_index_ji(j, i, 's') + biom_cell_level  !< Vertical cell center closest to 1.1m
1215    k_wind = k
1216    IF( k_wind < 1_iwp ) THEN  ! Avoid horizontal u and v of 0.0 m/s close to ground
1217       k_wind = 1_iwp
1218    ENDIF
[3321]1219
[3448]1220!-- Determine local values:
1221    IF ( average_input ) THEN
1222!--    Calculate ta from Tp assuming dry adiabatic laps rate
1223       ta = pt_av(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - cels_offs
[3321]1224
[3448]1225       vp = 0.034_wp
1226       IF ( humidity .AND. ALLOCATED( q_av ) ) THEN
1227          vp = q_av(k, j, i)
1228       ENDIF
[3321]1229
[3448]1230       ws = ( 0.5_wp * ABS( u_av(k_wind, j, i) + u_av(k_wind, j, i+1) )  +     &
1231          0.5_wp * ABS( v_av(k_wind, j, i) + v_av(k_wind, j+1, i) )  +         &
1232          0.5_wp * ABS( w_av(k_wind, j, i) + w_av(k_wind+1, j, i) ) )
1233    ELSE
1234!-- Calculate ta from Tp assuming dry adiabatic laps rate
1235       ta = pt(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - cels_offs
[3321]1236
[3448]1237       vp = q(k, j, i)
[3321]1238
[3448]1239       ws = ( 0.5_wp * ABS( u(k_wind, j, i) + u(k_wind, j, i+1) )  +           &
1240          0.5_wp * ABS( v(k_wind, j, i) + v(k_wind, j+1, i) )  +               &
1241          0.5_wp * ABS( w(k_wind, j, i) + w(k_wind+1, j, i) ) )
[3321]1242
[3448]1243    ENDIF
[3321]1244
[3448]1245!-- Local air pressure
1246    pair = surface_pressure
1247!
1248!-- Calculate water vapour pressure at saturation and convert to hPa
1249!-- The magnus formula is limited to temperatures up to 333.15 K to
1250!   avoid negative values of vp_sat
1251    vp_sat = 0.01_wp * magnus( MIN( ta + cels_offs, 333.15_wp ) )
1252    vp  = vp * pair / ( vp + 0.622_wp )
1253    IF ( vp > vp_sat ) vp = vp_sat
1254
1255    tmrt = ta
1256    IF ( radiation ) THEN
1257       CALL calculate_tmrt_2_directions (rad_sw_in(0, j, i),                   &
1258          rad_sw_out(0, j, i), rad_lw_in(0, j, i), rad_lw_out(0, j, i), ta,    &
1259          tmrt )
1260    ENDIF
1261
1262 END SUBROUTINE biom_determine_input_at
1263
1264
[3321]1265!------------------------------------------------------------------------------!
1266! Description:
1267! ------------
[3448]1268!> Calculate static thermal indices for any point within a 2D grid
1269!> time_integration.f90: 1065ff
[3321]1270!------------------------------------------------------------------------------!
[3448]1271 SUBROUTINE biom_calculate_static_grid ( average_input )
[3321]1272
[3448]1273    IMPLICIT NONE
[3321]1274
[3448]1275!-- Input attributes
1276    LOGICAL, INTENT ( IN ) ::  average_input  !< Calculate based on averaged input? conditions?
[3321]1277
[3448]1278!-- Internal variables
1279    INTEGER(iwp) ::  i, j      !< Running index
[3321]1280
[3448]1281    REAL(wp) ::  ta            !< Air temperature                  (°C)
1282    REAL(wp) ::  vp            !< Vapour pressure                  (hPa)
1283    REAL(wp) ::  ws            !< Wind speed    (local level)      (m/s)
1284    REAL(wp) ::  pair          !< Air pressure                     (hPa)
1285    REAL(wp) ::  tmrt_tmp      !< Mean radiant temperature
1286    REAL(wp) ::  pt_tmp        !< Perceived temperature
1287    REAL(wp) ::  utci_tmp      !< Universal thermal climate index
1288    REAL(wp) ::  pet_tmp       !< Physiologically equivalent temperature
[3321]1289
1290
[3448]1291    CALL biom_init_arrays
[3321]1292
[3448]1293    DO i = nxl, nxr
1294       DO j = nys, nyn
1295!--       Determine local input conditions
1296          CALL biom_determine_input_at ( average_input, i, j, ta, vp, ws,      &
1297             pair, tmrt_tmp ) 
1298          tmrt_grid(j, i) = tmrt_tmp
[3321]1299
[3448]1300!--       Only proceed if tmrt is available
1301          IF ( tmrt_tmp <= -998._wp ) THEN
1302             pt_tmp   = -999._wp
1303             utci_tmp = -999._wp
1304             pet_tmp  = -999._wp
1305             CYCLE
1306          END IF
[3321]1307
[3448]1308!--       Calculate static thermal indices based on local tmrt
1309          CALL calculate_static_thermal_indices ( ta, vp, ws,                  &
1310             pair, tmrt_tmp, pt_tmp, utci_tmp, pet_tmp )
[3321]1311
[3448]1312          IF ( average_input ) THEN
1313!--          Write results for selected averaged indices
1314             IF ( biom_pt_av )  THEN
1315                pt_av_grid(j, i)   = pt_tmp
1316             END IF
1317             IF ( biom_utci_av ) THEN
1318                utci_av_grid(j, i) = utci_tmp
1319             END IF
1320             IF ( biom_pet_av ) THEN
1321                pet_av_grid(j, i)  = pet_tmp
1322             END IF
1323          ELSE
1324!--          Write result for selected indices
1325             IF ( biom_pt )  THEN
1326                pt_grid(j, i)   = pt_tmp
1327             END IF
1328             IF ( biom_utci ) THEN
1329                utci_grid(j, i) = utci_tmp
1330             END IF
1331             IF ( biom_pet ) THEN
1332                pet_grid(j, i)  = pet_tmp
1333             END IF
1334          END IF
[3321]1335
[3448]1336       END DO
1337    END DO
[3321]1338
[3448]1339 END SUBROUTINE biom_calculate_static_grid
[3321]1340
[3448]1341!------------------------------------------------------------------------------!
1342! Description:
1343! ------------
1344!> Calculate dynamic thermal indices (currently only iPT, but expandable)
1345!------------------------------------------------------------------------------!
1346 SUBROUTINE biom_calc_ipt( ta, vp, ws, pair, tmrt, dt, energy_storage,         &
1347    t_clo, clo, actlev, age, weight, height, work, sex, ipt )
[3321]1348
[3448]1349    IMPLICIT NONE
[3321]1350
[3448]1351!-- Input parameters
1352    REAL(wp), INTENT ( IN )  ::  ta   !< Air temperature                  (°C)
1353    REAL(wp), INTENT ( IN )  ::  vp   !< Vapour pressure                  (hPa)
1354    REAL(wp), INTENT ( IN )  ::  ws   !< Wind speed    (local level)      (m/s)
1355    REAL(wp), INTENT ( IN )  ::  pair !< Air pressure                     (hPa)
1356    REAL(wp), INTENT ( IN )  ::  tmrt !< Mean radiant temperature         (°C)
1357    REAL(wp), INTENT ( IN )  ::  dt   !< Time past since last calculation (s)
1358    REAL(wp), INTENT ( IN )  ::  age  !< Age of agent                     (y)
1359    REAL(wp), INTENT ( IN )  ::  weight  !< Weight of agent               (Kg)
1360    REAL(wp), INTENT ( IN )  ::  height  !< Height of agent               (m)
1361    REAL(wp), INTENT ( IN )  ::  work    !< Mechanical workload of agent
1362                                         !  (without metabolism!)         (W)
1363    INTEGER(iwp), INTENT ( IN ) ::  sex  !< Sex of agent (1 = male, 2 = female)
[3321]1364
[3448]1365!-- Both, input and output parameters
1366    Real(wp), INTENT ( INOUT )  ::  energy_storage    !< Energy storage   (W/m²)
1367    Real(wp), INTENT ( INOUT )  ::  t_clo   !< Clothing temperature       (°C)
1368    Real(wp), INTENT ( INOUT )  ::  clo     !< Current clothing in sulation
1369    Real(wp), INTENT ( INOUT )  ::  actlev  !< Individuals activity level
1370                                            !  per unit surface area      (W/m²)
1371!-- Output parameters
1372    REAL(wp), INTENT ( OUT ) ::  ipt    !< Instationary perceived temp.   (°C)
[3321]1373
[3448]1374!-- If clo equals the initial value, this is the initial call
1375    IF ( clo <= -998._wp ) THEN
1376!--    Initialize instationary perceived temperature with personalized
1377!      PT as an initial guess, set actlev and clo
1378       CALL ipt_init ( age, weight, height, sex, work, actlev, clo,            &
1379          ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo,                   &
1380          ipt )
1381    ELSE
1382!--    Estimate local instatinoary perceived temperature
1383       CALL ipt_cycle ( ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo,     &
1384          clo, actlev, work, ipt )
1385    ENDIF
[3321]1386
[3448]1387 END SUBROUTINE biom_calc_ipt
1388
1389 END MODULE biometeorology_mod
Note: See TracBrowser for help on using the repository browser.