Home

Context Navigation

source: palm/trunk/SOURCE/uv_exposure_model_mod.f90 @ 2764

Last change on this file since 2764 was 2718, checked in by maronga, 7 years ago
deleting of deprecated files; headers updated where needed
Property svn:keywords set to `Id`
File size: 31.1 KB

Line
1	!> @file uv_exposure_model_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM 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 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 2017-2018 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: uv_exposure_model_mod.f90 2718 2018-01-02 08:49:38Z gronemeier $
27	! Initial revision
28	!
29	!
30	!
31	! Authors:
32	! --------
33	! @author Michael Schrempf
34	!
35	!
36	! Description:
37	! ------------
38	!> Calculation of vitamin-D weighted UV exposure
39	!>
40	!>
41	!> @todo uv_vitd3dose-->new output type necessary (cumulative)
42	!> @todo consider upwelling radiation
43	!>
44	!> @note <Enter notes on the module>
45	!>
46	!> @bug <Enter known bugs here>
47	!------------------------------------------------------------------------------!
48	MODULE uv_exposure_model_mod
49
50
51	USE kinds
52
53	!
54	!-- Load required variables from existing modules
55	!-- USE modulename, &
56	!-- ONLY: ...
57
58
59	IMPLICIT NONE
60
61
62	!
63	!-- Declare all global variables within the module (alphabetical order)
64	INTEGER(iwp) :: ai = 0 !< running index
65	INTEGER(iwp) :: clothing = 3 !< clothing (0=unclothed, 1=Arms,Hands Face free, 3=Hand, Face free)
66	INTEGER(iwp) :: consider_obstructions = 1 !< 0 = unobstructed scenario,
67	!< 1 = scenario where obstrcutions are considered
68	INTEGER(iwp) :: orientation_angle = 0 !< orientation of front/face of the human model
69	INTEGER(iwp) :: saa_in_south = 1 !< 0 = sun always in south direction, 1 = actual sun position
70	INTEGER(iwp) :: obstruction_direct_beam = 0 !< Obstruction information for direct beam
71	INTEGER(iwp) :: turn_to_sun = 1 !< 0 = orientation of human as in orientation_angle,
72	!< 1 = human is always orientated towards the sun
73	INTEGER(iwp) :: zi = 0 !< running index
74
75	INTEGER(iwp), DIMENSION(0:44) :: obstruction_temp1 = 0 !< temporary obstruction information
76	!< stored with ibset as logical information
77	INTEGER(iwp), DIMENSION(0:359) :: obstruction_temp2 = 0 !< temporary obstruction information restored values
78	!< from logical information, which where stored by ibset
79
80	INTEGER(iwp), DIMENSION(0:35,0:9) :: obstruction = 0 !< final obstruction information array for all
81	!< hemispherical directions
82
83	INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: obstruction_lookup_table !< lookup table of obstruction
84	!< information of entire domain
85
86	REAL(wp) :: diffuse_exposure = 0.0_wp !< calculated exposure by diffuse radiation
87	REAL(wp) :: direct_exposure = 0.0_wp !< calculated exposure by direct beam
88	REAL(wp) :: projection_area_direct_beam = 0.0_wp !< projection area for by direct beam
89	REAL(wp) :: saa = 180.0_wp !< solar azimuth angle
90	REAL(wp) :: startpos_human = 0.0_wp !< start position in azimuth direction for the
91	!< interpolation of the projection area
92	REAL(wp) :: startpos_saa_float = 0.0_wp !< start position in azimuth direction for the
93	!< interpolation of the radiance field
94	REAL(wp) :: sza = 20.0_wp !< solar zenith angle
95	REAL(wp) :: xfactor = 0.0_wp !< relative x-position used for interpolation
96	REAL(wp) :: yfactor = 0.0_wp !< relative y-position used for interpolation
97
98	REAL(wp), DIMENSION(0:2) :: irradiance = 0.0_wp !< iradiance values extracted from irradiance lookup table
99
100	REAL(wp), DIMENSION(0:2,0:90) :: irradiance_lookup_table = 0.0_wp !< irradiance lookup table contains values
101	!< for direct, diffuse and global component
102	REAL(wp), DIMENSION(0:35,0:9) :: integration_array = 0.0_wp
103	REAL(wp), DIMENSION(0:35,0:9) :: projection_area = 0.0_wp
104	REAL(wp), DIMENSION(0:35,0:9) :: projection_area_lookup_table = 0.0_wp
105	REAL(wp), DIMENSION(0:71,0:9) :: projection_area_temp = 0.0_wp
106	REAL(wp), DIMENSION(0:35,0:9) :: radiance_array = 0.0_wp
107	REAL(wp), DIMENSION(0:71,0:9) :: radiance_array_temp = 0.0_wp
108	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vitd3_exposure
109	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vitd3_exposure_av
110
111	REAL(wp), DIMENSION(0:35,0:9,0:90) :: radiance_lookup_table = 0.0_wp
112
113
114	SAVE
115
116	PRIVATE
117
118
119	!
120	!-- Add INTERFACES that must be available to other modules (alphabetical order)
121	PUBLIC uvem_3d_data_averaging, uvem_calc_exposure, uvem_check_data_output, &
122	uvem_data_output_2d, uvem_define_netcdf_grid, uvem_init, &
123	uvem_init_arrays, uvem_parin
124
125	!
126	!-- Add VARIABLES that must be available to other modules (alphabetical order)
127	! PUBLIC
128
129	!
130	!-- Add PROGNOSTIC QUANTITIES that must be available to other modules (alphabetical order)
131	!-- PUBLIC ...
132
133
134	!
135	!-- Default procedures for all new modules (not all are necessarily required,
136	!-- alphabetical order is not essential)
137
138	!
139	!-- PALM interfaces:
140	!-- Data output checks for 2D/3D data to be done in check_parameters
141	INTERFACE uvem_check_data_output
142	MODULE PROCEDURE uvem_check_data_output
143	END INTERFACE uvem_check_data_output
144
145	! !
146	! !-- Data output checks for profile data to be done in check_parameters
147	! INTERFACE uvem_check_data_output_pr
148	! MODULE PROCEDURE uvem_check_data_output_pr
149	! END INTERFACE uvem_check_data_output_pr
150	!
151	! !
152	! !-- Input parameter checks to be done in check_parameters
153	! INTERFACE uvem_check_parameters
154	! MODULE PROCEDURE uvem_check_parameters
155	! END INTERFACE uvem_check_parameters
156	!
157	!
158	!-- Averaging of 3D data for output
159	INTERFACE uvem_3d_data_averaging
160	MODULE PROCEDURE uvem_3d_data_averaging
161	END INTERFACE uvem_3d_data_averaging
162	!
163	! !
164	!-- Data output of 2D quantities
165	INTERFACE uvem_data_output_2d
166	MODULE PROCEDURE uvem_data_output_2d
167	END INTERFACE uvem_data_output_2d
168	!
169	! !
170	! !-- Data output of 3D data
171	! INTERFACE uvem_data_output_3d
172	! MODULE PROCEDURE uvem_data_output_3d
173	! END INTERFACE uvem_data_output_3d
174	!
175	! !
176	! !-- Definition of data output quantities
177	INTERFACE uvem_define_netcdf_grid
178	MODULE PROCEDURE uvem_define_netcdf_grid
179	END INTERFACE uvem_define_netcdf_grid
180	!
181	! !
182	! !-- Output of information to the header file
183	! INTERFACE uvem_header
184	! MODULE PROCEDURE uvem_header
185	! END INTERFACE uvem_header
186
187	!
188	!-- Initialization actions
189	INTERFACE uvem_init
190	MODULE PROCEDURE uvem_init
191	END INTERFACE uvem_init
192
193	!
194	!-- Initialization of arrays
195	INTERFACE uvem_init_arrays
196	MODULE PROCEDURE uvem_init_arrays
197	END INTERFACE uvem_init_arrays
198	!
199	! !
200	! !-- Writing of binary output for restart runs !!! renaming?!
201	! INTERFACE uvem_last_actions
202	! MODULE PROCEDURE uvem_last_actions
203	! END INTERFACE uvem_last_actions
204	!
205	!
206	!-- Reading of NAMELIST parameters
207	INTERFACE uvem_parin
208	MODULE PROCEDURE uvem_parin
209	END INTERFACE uvem_parin
210	!
211	! !
212	! !-- Reading of parameters for restart runs
213	! INTERFACE uvem_read_restart_data
214	! MODULE PROCEDURE uvem_read_restart_data
215	! END INTERFACE uvem_read_restart_data
216	!
217	! !
218	! !-- Writing of parameters for restart runs !!! renaming?!
219	! INTERFACE uvem_skip_var_list
220	! MODULE PROCEDURE uvem_skip_var_list
221	! END INTERFACE uvem_skip_var_list
222	!
223	! !
224	! !-- Swapping of time levels (required for prognostic variables)
225	! INTERFACE uvem_swap_timelevel
226	! MODULE PROCEDURE uvem_swap_timelevel
227	! END INTERFACE uvem_swap_timelevel
228	!
229	! !
230	! !-- New module-specific procedure(s) (alphabetical order):
231	! INTERFACE uvem_newprocedure
232	! MODULE PROCEDURE uvem_newprocedure
233	! END INTERFACE uvem_newprocedure
234
235	CONTAINS
236
237	!------------------------------------------------------------------------------!
238	! Description:
239	! ------------
240	!> Check data output for new module.
241	!------------------------------------------------------------------------------!
242	SUBROUTINE uvem_check_data_output( var, unit, i, ilen, k )
243
244	USE control_parameters, &
245	ONLY: data_output, message_string, uv_exposure
246
247	IMPLICIT NONE
248
249	CHARACTER (LEN=*) :: unit !< string for unit of output quantity
250	CHARACTER (LEN=*) :: var !< string for output quantity
251
252	INTEGER(iwp) :: i !<
253	INTEGER(iwp) :: ilen !<
254	INTEGER(iwp) :: k !<
255
256	SELECT CASE ( TRIM( var ) )
257
258
259	CASE ( 'uvem_vitd3*' )
260	IF ( .NOT. uv_exposure ) THEN
261	message_string = 'output of "' // TRIM( var ) // '" requi' // &
262	'res a namelist &uvexposure_par'
263	CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 )
264	ENDIF
265	IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN
266	message_string = 'illegal value for data_output: "' // &
267	TRIM( var ) // '" & only 2d-horizontal ' // &
268	'cross sections are allowed for this value'
269	CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
270	ENDIF
271	unit = 'IU/s'
272
273	CASE ( 'uvem_vitd3dose*' )
274	IF ( .NOT. uv_exposure ) THEN
275	message_string = 'output of "' // TRIM( var ) // '" requi' // &
276	'res a namelist &uvexposure_par'
277	CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 )
278	ENDIF
279	IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN
280	message_string = 'illegal value for data_output: "' // &
281	TRIM( var ) // '" & only 2d-horizontal ' // &
282	'cross sections are allowed for this value'
283	CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
284	ENDIF
285	unit = 'IU/av-h'
286
287	CASE DEFAULT
288	unit = 'illegal'
289
290	END SELECT
291
292	END SUBROUTINE uvem_check_data_output
293
294
295	!-----------------------------------------------------------------------------!
296	!
297	! Description:
298	! ------------
299	!> Subroutine defining 2D output variables
300	!-----------------------------------------------------------------------------!
301	SUBROUTINE uvem_data_output_2d( av, variable, found, grid, mode, local_pf, &
302	two_d, nzb_do, nzt_do )
303
304	USE indices
305
306	USE kinds
307
308
309	IMPLICIT NONE
310
311	CHARACTER (LEN=*) :: grid !<
312	CHARACTER (LEN=*) :: mode !<
313	CHARACTER (LEN=*) :: variable !<
314
315	INTEGER(iwp) :: av !<
316	INTEGER(iwp) :: i !< running index
317	INTEGER(iwp) :: j !< running index
318	INTEGER(iwp) :: k !< running index
319	INTEGER(iwp) :: m !< running index
320	INTEGER(iwp) :: nzb_do !<
321	INTEGER(iwp) :: nzt_do !<
322
323	LOGICAL :: found !<
324	LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
325
326	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) :: local_pf !<
327
328
329	found = .TRUE.
330
331	SELECT CASE ( TRIM( variable ) )
332	!
333	!-- Before data is transfered to local_pf, transfer is it 2D dummy variable and exchange ghost points therein.
334	!-- However, at this point this is only required for instantaneous arrays, time-averaged quantities are already exchanged.
335	CASE ( 'uvem_vitd3*_xy' ) ! 2d-array
336	IF ( av == 0 ) THEN
337	DO i = nxl, nxr
338	DO j = nys, nyn
339	local_pf(i,j,nzb+1) = vitd3_exposure(j,i)
340	ENDDO
341	ENDDO
342	ENDIF
343
344	two_d = .TRUE.
345	grid = 'zu1'
346
347	CASE ( 'uvem_vitd3dose*_xy' ) ! 2d-array
348	IF ( av == 1 ) THEN
349	DO i = nxl, nxr
350	DO j = nys, nyn
351	local_pf(i,j,nzb+1) = vitd3_exposure_av(j,i)
352	ENDDO
353	ENDDO
354	ENDIF
355
356	two_d = .TRUE.
357	grid = 'zu1'
358
359	CASE DEFAULT
360	found = .FALSE.
361	grid = 'none'
362
363	END SELECT
364
365	END SUBROUTINE uvem_data_output_2d
366
367
368	!------------------------------------------------------------------------------!
369	!
370	! Description:
371	! ------------
372	!> Subroutine defining appropriate grid for netcdf variables.
373	!> It is called out from subroutine netcdf.
374	!------------------------------------------------------------------------------!
375	SUBROUTINE uvem_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
376
377	IMPLICIT NONE
378
379	CHARACTER (LEN=*), INTENT(IN) :: var !<
380	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<
381	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<
382	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<
383
384	LOGICAL, INTENT(OUT) :: found !<
385
386	found = .TRUE.
387
388	!
389	!-- Check for the grid
390	SELECT CASE ( TRIM( var ) )
391
392	CASE ( 'uvem_vitd3_xy', 'uvem_vitd3dose_xy' )
393	grid_x = 'x'
394	grid_y = 'y'
395	grid_z = 'zu1'
396
397	CASE DEFAULT
398	found = .FALSE.
399	grid_x = 'none'
400	grid_y = 'none'
401	grid_z = 'none'
402
403	END SELECT
404
405	END SUBROUTINE uvem_define_netcdf_grid
406
407
408	!------------------------------------------------------------------------------!
409	! Description:
410	! ------------
411	!> Parin for &uvexposure_par for UV exposure model
412	!------------------------------------------------------------------------------!
413	SUBROUTINE uvem_parin
414
415	USE control_parameters, &
416	ONLY: uv_exposure
417
418
419	IMPLICIT NONE
420
421	CHARACTER (LEN=80) :: line !< dummy string for current line in parameter file
422
423	NAMELIST /uvexposure_par/ clothing
424
425	line = ' '
426
427	!
428	!-- Try to find uv exposure model namelist
429	REWIND ( 11 )
430	line = ' '
431	DO WHILE ( INDEX( line, '&uvexposure_par' ) == 0 )
432	READ ( 11, '(A)', END=10 ) line
433	ENDDO
434	BACKSPACE ( 11 )
435
436	!
437	!-- Read user-defined namelist
438	READ ( 11, uvexposure_par )
439
440	!
441	!-- Set flag that indicates that the uv exposure model is switched on
442	uv_exposure = .TRUE.
443
444
445	10 CONTINUE
446
447
448	END SUBROUTINE uvem_parin
449
450
451	!------------------------------------------------------------------------------!
452	!
453	! Description:
454	! ------------
455	!> Subroutine for averaging 3D data
456	!------------------------------------------------------------------------------!
457	SUBROUTINE uvem_3d_data_averaging( mode, variable )
458
459
460	USE control_parameters
461
462	USE indices
463
464	USE kinds
465
466	IMPLICIT NONE
467
468	CHARACTER (LEN=*) :: mode !<
469	CHARACTER (LEN=*) :: variable !<
470
471	INTEGER(iwp) :: i !<
472	INTEGER(iwp) :: j !<
473	INTEGER(iwp) :: k !<
474	INTEGER(iwp) :: m !< running index
475
476	IF ( mode == 'allocate' ) THEN
477
478	SELECT CASE ( TRIM( variable ) )
479
480	CASE ( 'uvem_vitd3dose*' )
481	IF ( .NOT. ALLOCATED( vitd3_exposure_av ) ) THEN
482	ALLOCATE( vitd3_exposure_av(nysg:nyng,nxlg:nxrg) )
483	ENDIF
484	vitd3_exposure_av = 0.0_wp
485
486
487	CASE DEFAULT
488	CONTINUE
489
490	END SELECT
491
492	ELSEIF ( mode == 'sum' ) THEN
493
494	SELECT CASE ( TRIM( variable ) )
495
496	CASE ( 'uvem_vitd3dose*' )
497	DO i = nxlg, nxrg
498	DO j = nysg, nyng
499	vitd3_exposure_av(j,i) = vitd3_exposure_av(j,i) + vitd3_exposure(j,i)
500	ENDDO
501	ENDDO
502
503
504	CASE DEFAULT
505	CONTINUE
506
507	END SELECT
508
509	!
510	!-- No averaging since we are calculating a dose (only sum is calculated and saved to
511	!-- av.nc file)
512	! ELSEIF ( mode == 'average' ) THEN
513
514	ENDIF
515
516	END SUBROUTINE uvem_3d_data_averaging
517
518
519
520
521	!------------------------------------------------------------------------------!
522	! Description:
523	! ------------
524	!> Initialization of the new module
525	!------------------------------------------------------------------------------!
526	SUBROUTINE uvem_init
527
528
529	USE control_parameters, &
530	ONLY: initializing_actions
531
532	IMPLICIT NONE
533
534	!
535	!-- Actions for initial runs
536	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
537	OPEN(90, STATUS='old',FILE=&
538	'RADIANCE', FORM='UNFORMATTED')
539	READ(90) radiance_lookup_table
540	CLOSE(90)
541	OPEN(90, STATUS='old',FILE=&
542	'IRRADIANCE', FORM='UNFORMATTED')
543	READ(90) irradiance_lookup_table
544	CLOSE(90)
545
546	!________________________LOAD Obstruction information _______________________________
547	IF ( consider_obstructions == 1 ) THEN
548	OPEN(90, STATUS='old',FILE=&
549	'OBSTRUCTION', FORM='UNFORMATTED')
550	READ(90) obstruction_lookup_table
551	CLOSE(90)
552	ELSEIF( consider_obstructions == 0 ) THEN
553	obstruction(:,:) = 1
554	ENDIF
555
556	!________________________LOAD integration_array ________________________________________________________________________
557	OPEN(90, STATUS='old',FILE=&
558	'SOLIDANGLE', FORM='UNFORMATTED')
559	READ(90) integration_array
560	CLOSE(90)
561
562	IF (clothing==1) THEN
563	OPEN(90, STATUS='old',FILE='HUMAN', FORM='UNFORMATTED')
564	ENDIF
565	READ(90) projection_area_lookup_table
566	CLOSE(90)
567
568	!
569	!-- Actions for restart runs
570	ELSE
571	! ....
572
573	ENDIF
574
575	END SUBROUTINE uvem_init
576
577
578	!-----------------------------------------------------------------------------!
579	! Description:
580	! ------------
581	!> Allocate new module arrays and define pointers if required
582	!-----------------------------------------------------------------------------!
583	SUBROUTINE uvem_init_arrays
584
585	USE indices, &
586	ONLY: nxlg, nxrg, nyng, nysg
587
588
589	IMPLICIT NONE
590
591	!
592	!-- Allocate arrays
593	ALLOCATE ( vitd3_exposure(nysg:nyng,nxlg:nxrg) )
594	ALLOCATE ( vitd3_exposure_av(nysg:nyng,nxlg:nxrg) )
595	ALLOCATE ( obstruction_lookup_table(nxlg:nxrg,nysg:nyng,0:44) )
596
597	!
598	!-- Initialize arrays
599	vitd3_exposure = 0.0_wp
600	vitd3_exposure_av = 0.0_wp
601	obstruction_lookup_table = 0
602
603
604	END SUBROUTINE uvem_init_arrays
605
606	!------------------------------------------------------------------------------!
607	! Description:
608	! ------------
609	!> Module-specific routine for new module
610	!------------------------------------------------------------------------------!
611	SUBROUTINE uvem_solar_position
612
613	USE constants, &
614	ONLY: pi
615
616	USE date_and_time_mod, &
617	ONLY: calc_date_and_time, day_of_year, time_utc
618
619	USE control_parameters, &
620	ONLY: latitude, longitude
621
622	IMPLICIT NONE
623
624
625	REAL(wp) :: alpha = 0.0_wp !< solar azimuth angle in radiant
626	REAL(wp) :: declination = 0.0_wp !< declination
627	REAL(wp) :: dtor = 0.0_wp !< factor to convert degree to radiant
628	REAL(wp) :: js = 0.0_wp !< parameter for solar position calculation
629	REAL(wp) :: lat = 52.39_wp !< latitude
630	REAL(wp) :: lon = 9.7_wp !< longitude
631	REAL(wp) :: thetar = 0.0_wp !< angle for solar zenith angle calculation
632	REAL(wp) :: thetasr = 0.0_wp !< angle for solar azimuth angle calculation
633	REAL(wp) :: zgl = 0.0_wp !< calculated exposure by direct beam
634	REAL(wp) :: woz = 0.0_wp !< calculated exposure by diffuse radiation
635	REAL(wp) :: wsp = 0.0_wp !< calculated exposure by direct beam
636
637
638	CALL calc_date_and_time
639
640	dtor = pi / 180.0_wp
641	lat = latitude
642	lon = longitude
643	!
644	!-- calculation of js :
645	js= 72.0_wp * ( day_of_year + ( time_utc / 86400.0_wp ) ) / 73.0_wp
646	!
647	!-- calculation of equation of time (zgl):
648	zgl = 0.0066_wp + 7.3525_wp * cos( ( js + 85.9_wp ) * dtor ) + 9.9359_wp * &
649	cos( ( 2.0_wp * js + 108.9_wp ) * dtor ) + 0.3387_wp * cos( ( 3 * js + 105.2_wp ) * dtor )
650	!
651	!-- calculation of apparent solar time woz:
652	woz = ( ( time_utc / 3600.0_wp ) - ( 4.0_wp * ( 15.0_wp - lon ) ) / 60.0_wp ) + ( zgl / 60.0_wp )
653	!
654	!-- calculation of hour angle (wsp):
655	wsp = ( woz - 12.0_wp ) * 15.0_wp
656	!
657	!-- calculation of declination:
658	declination = 0.3948_wp - 23.2559_wp * cos( ( js + 9.1_wp ) * dtor ) - &
659	0.3915_wp * cos( ( 2.0_wp * js + 5.4_wp ) * dtor ) - 0.1764_wp * cos( ( 3.0_wp * js + 26.0_wp ) * dtor )
660	!
661	!-- calculation of solar zenith angle
662	thetar = acos( sin( lat * dtor) * sin( declination * dtor ) + cos( wsp * dtor ) * &
663	cos( lat * dtor ) * cos( declination * dtor ) )
664	thetasr = asin( sin( lat * dtor) * sin( declination * dtor ) + cos( wsp * dtor ) * &
665	cos( lat * dtor ) * cos( declination * dtor ) )
666	sza = thetar / dtor
667	!
668	!-- calculation of solar azimuth angle
669	IF (woz .LE. 12.0_wp) alpha = pi - acos( sin(thetasr) * sin(lat * dtor) - &
670	sin( declination * dtor ) / ( cos(thetasr) * cos( lat * dtor ) ) )
671	IF (woz .GT. 12.0_wp) alpha = pi + acos( sin(thetasr) * sin(lat * dtor) - &
672	sin( declination * dtor ) / ( cos(thetasr) * cos( lat * dtor ) ) )
673	saa = alpha / dtor
674
675	END SUBROUTINE uvem_solar_position
676
677
678	!------------------------------------------------------------------------------!
679	! Description:
680	! ------------
681	!> Module-specific routine for new module
682	!------------------------------------------------------------------------------!
683	SUBROUTINE uvem_calc_exposure
684
685	USE constants, &
686	ONLY: pi
687
688	USE indices, &
689	ONLY: nxlg, nxrg, nyng, nysg
690
691
692	IMPLICIT NONE
693
694	INTEGER(iwp) :: i !< running index
695	INTEGER(iwp) :: j !< running index
696
697
698	CALL uvem_solar_position
699
700	IF (sza .GE. 90) THEN
701	vitd3_exposure(:,:) = 0.0_wp
702	ELSE
703
704	!
705	!-- rotate 3D-Model human to desired direction -----------------------------
706	projection_area_temp( 0:35,:) = projection_area_lookup_table
707	projection_area_temp(36:71,:) = projection_area_lookup_table
708	IF ( Turn_to_Sun .EQ. 0 ) startpos_human = orientation_angle / 10.0_wp
709	IF ( Turn_to_Sun .EQ. 1 ) startpos_human = startpos_saa_float
710	DO ai = 0, 35
711	xfactor = ( startpos_human ) - INT( startpos_human )
712	DO zi = 0, 9
713	projection_area(ai,zi) = ( projection_area_temp( ai + INT( startpos_human ), zi) * &
714	(1.0_wp - xfactor ) ) &
715	+( projection_area_temp( ai + 1 + INT( startpos_human ), zi) * &
716	xfactor)
717	ENDDO
718	ENDDO
719	!
720	!-- calculate Projectionarea for direct beam -----------------------------'
721	projection_area_temp( 0:35,:) = projection_area
722	projection_area_temp(36:71,:) = projection_area
723	yfactor = ( sza / 10.0_wp ) - INT( sza / 10.0_wp )
724	xfactor = ( startpos_saa_float ) - INT( startpos_saa_float )
725	projection_area_direct_beam = ( projection_area_temp( INT(startpos_saa_float) ,INT( sza / 10.0_wp ) ) * &
726	( 1.0_wp - xfactor ) * ( 1.0_wp - yfactor ) ) + &
727	( projection_area_temp( INT(startpos_saa_float) + 1,INT(sza/10.0_wp)) *&
728	( xfactor ) * ( 1.0_wp - yfactor ) ) + &
729	( projection_area_temp( INT(startpos_saa_float) ,INT(sza/10.0_wp)+1)*&
730	( 1.0_wp - xfactor ) * ( yfactor ) ) + &
731	( projection_area_temp( INT(startpos_saa_float) + 1,INT(sza/10.0_wp)+1)*&
732	( xfactor ) * ( yfactor ) )
733
734
735	!
736	!-- interpolate to accurate Solar Zenith Angle ------------------
737	DO ai = 0, 35
738	xfactor = (sza)-INT(sza)
739	DO zi = 0, 9
740	radiance_array(ai,zi) = ( radiance_lookup_table(ai, zi, INT(sza) ) * ( 1.0_wp - xfactor) ) +&
741	( radiance_lookup_table(ai,zi,INT(sza) + 1) * xfactor )
742	ENDDO
743	ENDDO
744	Do ai = 0, 2
745	irradiance(ai) = ( irradiance_lookup_table(ai, INT(sza) ) * ( 1.0_wp - xfactor)) + &
746	(irradiance_lookup_table(ai, INT(sza) + 1) * xfactor )
747	ENDDO
748	!
749	!-- interpolate to accurate Solar Azimuth Angle ------------------
750	IF ( saa_in_south .EQ. 0 ) THEN
751	startpos_saa_float = 180.0_wp / 10.0_wp
752	ELSE
753	startpos_saa_float = saa / 10.0_wp
754	ENDIF
755	radiance_array_temp( 0:35,:) = radiance_array
756	radiance_array_temp(36:71,:) = radiance_array
757	xfactor = (startpos_saa_float) - INT(startpos_saa_float)
758	DO ai = 0, 35
759	DO zi = 0, 9
760	radiance_array(ai,zi) = ( radiance_array_temp( ai + INT( startpos_saa_float ), zi ) * &
761	(1.0_wp - xfactor)) &
762	+ ( radiance_array_temp( ai + 1 + INT( startpos_saa_float ), zi ) &
763	* xfactor)
764	ENDDO
765	ENDDO
766
767
768
769
770	DO i = nxlg, nxrg
771	DO j = nysg, nyng
772	!
773	!-- extract obstrcution from IBSET-Integer_Array ------------------'
774	IF (consider_obstructions == 1) THEN
775	obstruction_temp1 = obstruction_lookup_table(i,j,:)
776	IF (obstruction_temp1(0) .NE. 9) THEN
777	DO zi = 0, 44
778	DO ai = 0, 7
779	IF ( btest( obstruction_temp1(zi), ai ) .EQV. .TRUE.) THEN
780	obstruction_temp2( ( zi * 8 ) + ai ) = 1
781	ELSE
782	obstruction_temp2( ( zi * 8 ) + ai ) = 0
783	ENDIF
784	ENDDO
785	ENDDO
786	DO zi = 0, 9
787	obstruction(:,zi) = obstruction_temp2( zi * 36 :( zi * 36) + 35 )
788	ENDDO
789	ELSE
790	obstruction(:,:) = 0
791	ENDIF
792	ENDIF
793	!
794	!-- calculated human exposure ------------------'
795	diffuse_exposure = SUM( radiance_array * projection_area * integration_array * obstruction )
796
797	obstruction_direct_beam = obstruction( nint(startpos_saa_float), nint( sza / 10.0_wp ) )
798	IF (sza .GE. 89.99_wp) THEN
799	sza = 89.99999_wp
800	ENDIF
801	!
802	!-- calculate direct normal irradiance (direct beam) ------------------'
803	direct_exposure = ( irradiance(1) / cos( pi * sza / 180.0_wp ) ) * &
804	projection_area_direct_beam * obstruction_direct_beam
805
806	vitd3_exposure(j,i) = ( diffuse_exposure + direct_exposure ) / 1000.0_wp * 70.97_wp
807	! unit = international units vitamin D per second
808	ENDDO
809	ENDDO
810	ENDIF
811
812	END SUBROUTINE uvem_calc_exposure
813
814
815
816
817
818	! ------------------------------------------------------------------------------!
819	! Description:
820	! ------------
821	! > Check parameters routine
822	! ------------------------------------------------------------------------------!
823	! SUBROUTINE uvem_check_parameters
824	!
825	! USE control_parameters, &
826	! ONLY: message_string
827	!
828	!
829	! IMPLICIT NONE
830	!
831	!
832	! -- Checks go here (cf. check_parameters.f90).
833	!
834	! END SUBROUTINE uvem_check_parameters
835
836
837
838	! !------------------------------------------------------------------------------!
839	! ! Description:
840	! ! ------------
841	! !> Header output
842	! !------------------------------------------------------------------------------!
843	! SUBROUTINE uvem_header ( io )
844	!
845	!
846	! IMPLICIT NONE
847	!
848	!
849	! INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
850	!
851	! !
852	! !-- Header output goes here
853	! !-- ...
854	!
855	! END SUBROUTINE uvem_header
856
857
858
859	! !------------------------------------------------------------------------------!
860	! ! Description:
861	! ! ------------
862	! !> This routine reads the respective restart data.
863	! !------------------------------------------------------------------------------!
864	! SUBROUTINE uvem_read_restart_data
865	!
866	!
867	! IMPLICIT NONE
868	!
869	! CHARACTER (LEN=30) :: variable_chr !< dummy variable to read string
870	!
871	!
872	! READ ( 13 ) variable_chr
873	! DO WHILE ( TRIM( variable_chr ) /= '* end new module *' )
874	!
875	! SELECT CASE ( TRIM( variable_chr ) )
876	!
877	! ! CASE ( 'param1' )
878	! ! IF ( .NOT. ALLOCATED( param1 ) ) ALLOCATE( param1(...) )
879	! ! READ ( 13 ) param1
880	!
881	! END SELECT
882	!
883	! READ ( 13 ) variable_chr
884	!
885	! ENDDO
886	!
887	! END SUBROUTINE uvem_read_restart_data
888
889
890	! !------------------------------------------------------------------------------!
891	! ! Description:
892	! ! ------------
893	! !> This routine writes the respective restart data.
894	! !------------------------------------------------------------------------------!
895	! SUBROUTINE uvem_last_actions
896	!
897	! IMPLICIT NONE
898	!
899	! ! WRITE ( 14 ) 'param1 '
900	! ! WRITE ( 14 ) param1
901	!
902	!
903	! WRITE ( 14 ) '* end new module * '
904	!
905	! END SUBROUTINE uvem_last_actions
906
907
908	! !------------------------------------------------------------------------------!
909	! ! Description:
910	! ! ------------
911	! !> Skipping of parameters from restart file (binary format).
912	! !------------------------------------------------------------------------------!
913	! SUBROUTINE uvem_skip_var_list
914	!
915	! IMPLICIT NONE
916	!
917	! CHARACTER (LEN=1) :: param2 !<
918	! CHARACTER (LEN=30) :: variable_chr !<
919	!
920	! READ ( 13 ) variable_chr
921	! DO WHILE ( TRIM( variable_chr ) /= '* end new module *' )
922	! READ ( 13 ) param2
923	! ENDDO
924	!
925	! END SUBROUTINE uvem_skip_var_list
926
927
928	END MODULE uv_exposure_model_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |