Home

Context Navigation

biometeorology_mod.f90 @ 3523

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

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

Property svn:keywords set to Id

Property svn:mergeinfo set to (toggle deleted branches)

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

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

Download in other formats: