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

Last change on this file since 3469 was 3464, checked in by kanani, 5 years ago

from branch resler@3462: add MRT shaping function (radiation_model_mod), use basic constants (biometeorology_mod), adjust precision to wp (biometeorology_pt_mod)

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