source: palm/tags/release-6.0/SOURCE/biometeorology_mod.f90 @ 4141

Last change on this file since 4141 was 3475, checked in by kanani, 5 years ago

Add option for using MRT from RTM instead of simplified global value

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