Home

Context Navigation

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

Last change on this file since 2750 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

Rev	Line
[2687]	1	!> @file uv_exposure_model_mod.f90
[2696]	2	!------------------------------------------------------------------------------!
	3	! This file is part of the PALM model system.
[2687]	4	!
[2696]	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
[2687]	8	! version.
	9	!
[2696]	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.
[2687]	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	!
[2718]	17	! Copyright 2017-2018 Leibniz Universitaet Hannover
[2687]	18	!------------------------------------------------------------------------------!
	19	!
	20	! Current revisions:
	21	! -----------------
[2696]	22	!
	23	!
[2687]	24	! Former revisions:
	25	! -----------------
	26	! $Id: uv_exposure_model_mod.f90 2718 2018-01-02 08:49:38Z knoop $
[2696]	27	! Initial revision
[2687]	28	!
[2696]	29	!
	30	!
[2687]	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 |