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

Last change on this file since 3479 was 3479, checked in by gronemeier, 3 years ago

reworked check for output quantities; assign new palm-error numbers, set unit according to data standard

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