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

Last change on this file since 3516 was 3474, checked in by kanani, 5 years ago

Add netcdf input for uv exposure model (netcdf_data_input_mod, Makefile, .palm.iofiles), plus bugfixes (uv_exposure_model_mod)

  • Property svn:keywords set to Id
File size: 28.5 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 3474 2018-10-30 21:07:39Z gronemeier $
27! Schrempf:
28! Bugfix in rotate 3D-Model human to desired direction
29! Bugfix in interpolate to accurate Solar Azimuth Angle
30! Further changes in variables
31!
32! 3448 2018-10-29 18:14:31Z kanani
33! Temporary rename of namelist until uv model moves to biometeorology module
34!
35! 3274 2018-09-24 15:42:55Z knoop
36! Modularization of all bulk cloud physics code components
37!
38! 3248 2018-09-14 09:42:06Z sward
39! Minor formating changes
40!
41! 3246 2018-09-13 15:14:50Z sward
42! Added error handling for input namelist via parin_fail_message and small
43! typo bugfix
44!
45! 3241 2018-09-12 15:02:00Z raasch
46! unused variables removed
47!
48! 3014 2018-05-09 08:42:38Z maronga
49! Bugfix: domain bounds of local_pf corrected
50!
51! 3004 2018-04-27 12:33:25Z Giersch
52! Further allocation checks implemented (averaged data will be assigned to fill
53! values if no allocation happened so far)
54!
55! 2932 2018-03-26 09:39:22Z maronga
56! renamed uvexposure_par to biometeorology_parameters
57!
58! 2894 2018-03-15 09:17:58Z Giersch
59! Routine for skipping global restart data has been removed, uvem_last_actions
60! has been renamed to uvem_wrd_global and uvem_read_restart_data has been
61! renamed to uvem_rrd_global, variable named found has been introduced for
62! checking if restart data was found, reading of restart strings has been moved
63! completely to read_restart_data_mod, marker *** end new module *** has been
64! removed, strings and their respective lengths are written out and read now
65! in case of restart runs to get rid of prescribed character lengths, CASE
66! DEFAULT was added if restart data is read
67!
68! 2848 2018-03-05 10:11:18Z Giersch
69! Initial revision
70!
71!
72!
73! Authors:
74! --------
75! @author Michael Schrempf
76!
77!
78! Description:
79! ------------
80!> Calculation of vitamin-D weighted UV exposure
81!>
82!>
83!> @todo uv_vitd3dose-->new output type necessary (cumulative)
84!> @todo consider upwelling radiation
85!>
86!> @note <Enter notes on the module>
87!>
88!> @bug  <Enter known bugs here>
89!------------------------------------------------------------------------------!
90 MODULE uv_exposure_model_mod
91
92
93    USE basic_constants_and_equations_mod,                                     &
94        ONLY:  pi
95
96    USE kinds
97
98    USE date_and_time_mod,                                                     &
99       ONLY:  calc_date_and_time, day_of_year, time_utc
100
101    USE netcdf_data_input_mod,                                                 &
102       ONLY:  netcdf_data_input_uvem, uvem_projarea_f, uvem_radiance_f,        &
103               uvem_irradiance_f, uvem_integration_f, building_obstruction_f 
104
105
106
107    IMPLICIT NONE
108
109
110!
111!-- Declare all global variables within the module (alphabetical order)
112    INTEGER(iwp) ::  ai                      = 0 !< running index
113    INTEGER(iwp) ::  clothing                = 1 !< clothing (0=unclothed, 1=Arms,Hands Face free, 3=Hand, Face free)
114    INTEGER(iwp) ::  obstruction_direct_beam = 0 !< Obstruction information for direct beam   
115    INTEGER(iwp) ::  zi                      = 0 !< running index
116
117    INTEGER(KIND=1), DIMENSION(0:44)  ::  obstruction_temp1 = 0 !< temporary obstruction information
118                                                             !< stored with ibset as logical information
119    INTEGER(iwp), DIMENSION(0:359) ::  obstruction_temp2 = 0 !< temporary obstruction information restored values
120                                                             !< from logical information, which where stored by ibset 
121   
122    INTEGER(iwp), DIMENSION(0:35,0:9) ::  obstruction = 1 !< final obstruction information array for all
123                                                          !< hemispherical directions
124   
125    LOGICAL ::  consider_obstructions = .TRUE.   !< namelist parameter (see documentation)
126    LOGICAL ::  sun_in_south          = .FALSE.  !< namelist parameter (see documentation)
127    LOGICAL ::  turn_to_sun           = .TRUE.   !< namelist parameter (see documentation)
128
129    REAL(wp) ::  diffuse_exposure            =   0.0_wp !< calculated exposure by diffuse radiation
130    REAL(wp) ::  direct_exposure             =   0.0_wp !< calculated exposure by direct beam   
131    REAL(wp) ::  orientation_angle           =   0.0_wp !< orientation of front/face of the human model   
132    REAL(wp) ::  projection_area_direct_beam =   0.0_wp !< projection area for  by direct beam
133    REAL(wp) ::  saa                         = 180.0_wp !< solar azimuth angle
134    REAL(wp) ::  startpos_human              =   0.0_wp !< start position in azimuth direction for the
135                                                        !< interpolation of the projection area             
136    REAL(wp) ::  startpos_saa_float          =   0.0_wp !< start position in azimuth direction for the
137                                                        !< interpolation of the radiance field             
138    REAL(wp) ::  sza                         =  20.0_wp !< solar zenith angle
139    REAL(wp) ::  xfactor                     =   0.0_wp !< relative x-position used for interpolation
140    REAL(wp) ::  yfactor                     =   0.0_wp !< relative y-position used for interpolation
141   
142    REAL(wp), DIMENSION(0:2)  ::  irradiance =   0.0_wp !< iradiance values extracted from irradiance lookup table   
143    REAL(wp), DIMENSION(0:2,0:90) ::  irradiance_lookup_table = 0.0_wp !< irradiance lookup table contains values
144                                                                       !< for direct, diffuse and global component
145
146    REAL(wp), DIMENSION(0:35,0:9) ::  integration_array            = 0.0_wp
147    REAL(wp), DIMENSION(0:35,0:9) ::  projection_area              = 0.0_wp
148    REAL(wp), DIMENSION(0:35,0:9) ::  projection_area_lookup_table = 0.0_wp
149    REAL(wp), DIMENSION(0:71,0:9) ::  projection_area_direct_temp  = 0.0_wp
150    REAL(wp), DIMENSION(0:71,0:9) ::  projection_area_temp         = 0.0_wp
151    REAL(wp), DIMENSION(0:35,0:9) ::  radiance_array               = 0.0_wp
152    REAL(wp), DIMENSION(0:71,0:9) ::  radiance_array_temp          = 0.0_wp
153    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vitd3_exposure
154    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vitd3_exposure_av
155
156    REAL(wp), DIMENSION(0:35,0:9,0:90) ::  radiance_lookup_table = 0.0_wp
157
158
159    SAVE
160
161    PRIVATE
162
163   
164!
165!-- Add INTERFACES that must be available to other modules (alphabetical order)
166    PUBLIC uvem_3d_data_averaging, uvem_calc_exposure, uvem_check_data_output, &
167           uvem_data_output_2d, uvem_define_netcdf_grid, uvem_init,            &
168           uvem_init_arrays, uvem_parin
169
170!
171!-- Add VARIABLES that must be available to other modules (alphabetical order)
172!     PUBLIC
173
174!
175!-- Add PROGNOSTIC QUANTITIES that must be available to other modules (alphabetical order)
176!-- PUBLIC ...
177
178
179!
180!-- Default procedures for all new modules (not all are necessarily required,
181!-- alphabetical order is not essential)
182
183!
184!-- PALM interfaces:
185!-- Data output checks for 2D/3D data to be done in check_parameters
186    INTERFACE uvem_check_data_output
187       MODULE PROCEDURE uvem_check_data_output
188    END INTERFACE uvem_check_data_output
189!
190!
191!-- Averaging of 3D data for output
192    INTERFACE uvem_3d_data_averaging
193       MODULE PROCEDURE uvem_3d_data_averaging
194    END INTERFACE uvem_3d_data_averaging
195!
196!
197!-- Data output of 2D quantities
198    INTERFACE uvem_data_output_2d
199       MODULE PROCEDURE uvem_data_output_2d
200    END INTERFACE uvem_data_output_2d
201!
202!
203! !-- Definition of data output quantities
204    INTERFACE uvem_define_netcdf_grid
205       MODULE PROCEDURE uvem_define_netcdf_grid
206    END INTERFACE uvem_define_netcdf_grid
207
208!
209!-- Initialization actions 
210    INTERFACE uvem_init
211       MODULE PROCEDURE uvem_init
212    END INTERFACE uvem_init
213!
214!
215! !-- Initialization of arrays
216    INTERFACE uvem_init_arrays
217       MODULE PROCEDURE uvem_init_arrays
218    END INTERFACE uvem_init_arrays
219!     
220!
221! !-- Reading of NAMELIST parameters
222    INTERFACE uvem_parin
223       MODULE PROCEDURE uvem_parin
224    END INTERFACE uvem_parin
225
226 CONTAINS
227
228!------------------------------------------------------------------------------!
229! Description:
230! ------------
231!> Check data output for new module.
232!------------------------------------------------------------------------------!
233 SUBROUTINE uvem_check_data_output( var, unit, i, ilen, k )
234 
235    USE control_parameters,                                                    &
236        ONLY:  data_output, message_string, uv_exposure
237
238    IMPLICIT NONE
239
240    CHARACTER (LEN=*) ::  unit     !< string for unit of output quantity
241    CHARACTER (LEN=*) ::  var      !< string for output quantity
242
243    INTEGER(iwp) ::  i      !<
244    INTEGER(iwp) ::  ilen   !<   
245    INTEGER(iwp) ::  k      !<
246
247    SELECT CASE ( TRIM( var ) )
248
249
250       CASE ( 'uvem_vitd3*' )
251          IF (  .NOT.  uv_exposure )  THEN
252             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
253                      'res a namelist &uvexposure_par'
254             CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 )
255          ENDIF
256          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
257             message_string = 'illegal value for data_output: "' //            &
258                              TRIM( var ) // '" & only 2d-horizontal ' //      &
259                              'cross sections are allowed for this value'
260             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
261          ENDIF
262          unit = 'IU/s'
263
264       CASE ( 'uvem_vitd3dose*' )
265          IF (  .NOT.  uv_exposure )  THEN
266             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
267                      'res  a namelist &uvexposure_par'
268             CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 )
269          ENDIF
270          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
271             message_string = 'illegal value for data_output: "' //            &
272                              TRIM( var ) // '" & only 2d-horizontal ' //      &
273                              'cross sections are allowed for this value'
274             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
275          ENDIF
276          unit = 'IU/av-h'
277             
278       CASE DEFAULT
279          unit = 'illegal'
280
281    END SELECT
282
283 END SUBROUTINE uvem_check_data_output
284
285
286!-----------------------------------------------------------------------------!
287!
288! Description:
289! ------------
290!> Subroutine defining 2D output variables
291!-----------------------------------------------------------------------------!
292 SUBROUTINE uvem_data_output_2d( av, variable, found, grid, local_pf, two_d,   &
293                                 nzb_do, nzt_do )
294 
295    USE indices
296
297    USE kinds
298
299
300    IMPLICIT NONE
301
302    CHARACTER (LEN=*) ::  grid     !<
303    CHARACTER (LEN=*) ::  variable !<
304
305    INTEGER(iwp) ::  av      !<
306    INTEGER(iwp) ::  i       !< running index
307    INTEGER(iwp) ::  j       !< running index
308    INTEGER(iwp) ::  nzb_do  !<
309    INTEGER(iwp) ::  nzt_do  !<
310
311    LOGICAL      ::  found !<
312    LOGICAL      ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
313
314    REAL(wp) ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
315
316    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
317
318
319    found = .TRUE.
320
321    SELECT CASE ( TRIM( variable ) )
322!
323!--    Before data is transfered to local_pf, transfer is it 2D dummy variable and exchange ghost points therein.
324!--    However, at this point this is only required for instantaneous arrays, time-averaged quantities are already exchanged.
325       CASE ( 'uvem_vitd3*_xy' )        ! 2d-array
326          IF ( av == 0 )  THEN
327             DO  i = nxl, nxr
328                DO  j = nys, nyn
329                   local_pf(i,j,nzb+1) = vitd3_exposure(j,i)
330                ENDDO
331             ENDDO
332          ENDIF
333
334          two_d = .TRUE.
335          grid = 'zu1'
336
337       CASE ( 'uvem_vitd3dose*_xy' )        ! 2d-array
338          IF ( .NOT. ALLOCATED( vitd3_exposure_av ) ) THEN
339             ALLOCATE( vitd3_exposure_av(nysg:nyng,nxlg:nxrg) )
340             vitd3_exposure_av = REAL( fill_value, KIND = wp )
341          ENDIF
342          IF ( av == 1 )  THEN
343             DO  i = nxl, nxr
344                DO  j = nys, nyn
345                   local_pf(i,j,nzb+1) = vitd3_exposure_av(j,i)
346                ENDDO
347             ENDDO
348          ENDIF
349
350          two_d = .TRUE.
351          grid = 'zu1'
352
353       CASE DEFAULT
354          found = .FALSE.
355          grid  = 'none'
356
357    END SELECT
358 
359 END SUBROUTINE uvem_data_output_2d
360
361
362!------------------------------------------------------------------------------!
363!
364! Description:
365! ------------
366!> Subroutine defining appropriate grid for netcdf variables.
367!> It is called out from subroutine netcdf.
368!------------------------------------------------------------------------------!
369 SUBROUTINE uvem_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
370   
371    IMPLICIT NONE
372
373    CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
374    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
375    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
376    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
377
378    LOGICAL, INTENT(OUT)           ::  found       !<
379
380    found  = .TRUE.
381
382!
383!-- Check for the grid
384    SELECT CASE ( TRIM( var ) )
385
386       CASE ( 'uvem_vitd3*_xy', 'uvem_vitd3dose*_xy' )
387          grid_x = 'x'
388          grid_y = 'y'
389          grid_z = 'zu1'
390
391       CASE DEFAULT
392          found  = .FALSE.
393          grid_x = 'none'
394          grid_y = 'none'
395          grid_z = 'none'
396
397    END SELECT
398
399 END SUBROUTINE uvem_define_netcdf_grid
400
401
402!------------------------------------------------------------------------------!
403! Description:
404! ------------
405!> Parin for &uvexposure_par for UV exposure model
406!------------------------------------------------------------------------------!
407 SUBROUTINE uvem_parin
408
409    USE control_parameters,                                                   &
410        ONLY:  uv_exposure
411
412
413    IMPLICIT NONE
414
415    CHARACTER (LEN=80) ::  line  !< dummy string for current line in parameter file
416   
417    NAMELIST /biometeorology_parameters_uv/  clothing,                            &
418                                             consider_obstructions,               &
419                                             orientation_angle,                   &
420                                             sun_in_south,                        &
421                                             turn_to_sun
422   
423!    line = ' ' 
424!
425!-- Try to find uv exposure model namelist
426    REWIND ( 11 )
427    line = ' '
428    DO WHILE ( INDEX( line, '&biometeorology_parameters_uv' ) == 0 )
429       READ ( 11, '(A)', END=20 )  line
430    ENDDO
431    BACKSPACE ( 11 )
432
433!
434!-- Read user-defined namelist
435    READ ( 11, biometeorology_parameters_uv, ERR = 10, END = 20 )
436
437!
438!-- Set flag that indicates that the uv exposure model is switched on
439    uv_exposure = .TRUE.
440    GOTO 20
441
442 10 BACKSPACE( 11 )
443    READ( 11 , '(A)') line
444    CALL parin_fail_message( 'biometeorology_parameters_uv', line )
445
446
447 20 CONTINUE
448       
449
450 END SUBROUTINE uvem_parin
451
452
453!------------------------------------------------------------------------------!
454!
455! Description:
456! ------------
457!> Subroutine for averaging 3D data
458!------------------------------------------------------------------------------!
459 SUBROUTINE uvem_3d_data_averaging( mode, variable )
460 
461
462    USE control_parameters
463
464    USE indices
465
466    USE kinds
467
468    IMPLICIT NONE
469
470    CHARACTER (LEN=*) ::  mode    !<
471    CHARACTER (LEN=*) :: variable !<
472
473    INTEGER(iwp) ::  i       !<
474    INTEGER(iwp) ::  j       !<
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             IF ( ALLOCATED( vitd3_exposure_av ) ) THEN
498                DO  i = nxlg, nxrg
499                   DO  j = nysg, nyng
500                      vitd3_exposure_av(j,i) = vitd3_exposure_av(j,i) + vitd3_exposure(j,i)
501                   ENDDO
502                ENDDO
503             ENDIF
504
505
506          CASE DEFAULT
507             CONTINUE
508
509       END SELECT
510
511!
512!-- No averaging since we are calculating a dose (only sum is calculated and saved to
513!-- av.nc file)
514!    ELSEIF ( mode == 'average' )  THEN
515
516    ENDIF
517
518 END SUBROUTINE uvem_3d_data_averaging
519
520   
521
522
523!------------------------------------------------------------------------------!
524! Description:
525! ------------
526!> Initialization of the new module
527!------------------------------------------------------------------------------!
528 SUBROUTINE uvem_init
529   
530
531    USE control_parameters,                                                   &
532        ONLY:  initializing_actions
533
534    USE netcdf_data_input_mod,                                                &
535        ONLY:  netcdf_data_input_uvem, uvem_projarea_f, uvem_radiance_f,      &
536               uvem_irradiance_f, uvem_integration_f, building_obstruction_f 
537
538    IMPLICIT NONE
539
540    CALL netcdf_data_input_uvem
541
542 END SUBROUTINE uvem_init
543
544
545!-----------------------------------------------------------------------------!
546! Description:
547! ------------
548!> Allocate new module arrays and define pointers if required
549!-----------------------------------------------------------------------------!
550 SUBROUTINE uvem_init_arrays
551   
552    USE indices,                                                              &
553        ONLY:  nxlg, nxrg, nyng, nysg
554
555
556
557    IMPLICIT NONE
558
559!
560!-- Allocate arrays
561    ALLOCATE ( vitd3_exposure(nysg:nyng,nxlg:nxrg) )
562    ALLOCATE ( vitd3_exposure_av(nysg:nyng,nxlg:nxrg) )
563
564!
565!-- Initialize arrays
566    vitd3_exposure = 0.0_wp
567    vitd3_exposure_av = 0.0_wp
568
569
570 END SUBROUTINE uvem_init_arrays
571
572!------------------------------------------------------------------------------!
573! Description:
574! ------------
575!> Module-specific routine for new module
576!------------------------------------------------------------------------------!
577 SUBROUTINE uvem_solar_position
578   
579    USE date_and_time_mod,                                                     &
580       ONLY:  calc_date_and_time, day_of_year, time_utc 
581   
582    USE control_parameters,                                                    &
583       ONLY:  latitude, longitude   
584
585    IMPLICIT NONE
586   
587   
588    REAL(wp) ::  alpha       = 0.0_wp   !< solar azimuth angle in radiant   
589    REAL(wp) ::  doy_r       = 0.0_wp   !< real format of day_of_year           
590    REAL(wp) ::  declination = 0.0_wp   !< declination
591    REAL(wp) ::  dtor        = 0.0_wp   !< factor to convert degree to radiant
592    REAL(wp) ::  js          = 0.0_wp   !< parameter for solar position calculation
593    REAL(wp) ::  lat         = 52.39_wp !< latitude
594    REAL(wp) ::  lon         = 9.7_wp   !< longitude       
595    REAL(wp) ::  thetar      = 0.0_wp   !< angle for solar zenith angle calculation
596    REAL(wp) ::  thetasr     = 0.0_wp   !< angle for solar azimuth angle calculation   
597    REAL(wp) ::  zgl         = 0.0_wp   !< calculated exposure by direct beam   
598    REAL(wp) ::  woz         = 0.0_wp   !< calculated exposure by diffuse radiation
599    REAL(wp) ::  wsp         = 0.0_wp   !< calculated exposure by direct beam   
600   
601
602    CALL calc_date_and_time
603    doy_r = real(day_of_year)   
604    dtor = pi / 180.0_wp
605    lat = latitude
606    lon = longitude
607!
608!-- calculation of js :
609    js=  72.0_wp * ( doy_r + ( time_utc / 86400.0_wp ) ) / 73.0_wp 
610!
611!-- calculation of equation of time (zgl):
612    zgl = 0.0066_wp + 7.3525_wp * cos( ( js + 85.9_wp ) * dtor ) + 9.9359_wp *                     &
613    cos( ( 2.0_wp * js + 108.9_wp ) * dtor ) + 0.3387_wp * cos( ( 3 * js + 105.2_wp ) * dtor )
614!
615!-- calculation of apparent solar time woz:
616    woz = ( ( time_utc / 3600.0_wp ) - ( 4.0_wp * ( 15.0_wp - lon ) ) / 60.0_wp ) + ( zgl / 60.0_wp )
617!
618!-- calculation of hour angle (wsp):
619    wsp = ( woz - 12.0_wp ) * 15.0_wp
620!
621!-- calculation of declination:
622    declination = 0.3948_wp - 23.2559_wp * cos( ( js + 9.1_wp ) * dtor ) -                         &
623    0.3915_wp * cos( ( 2.0_wp * js + 5.4_wp ) * dtor ) - 0.1764_wp * cos( ( 3.0_wp * js + 26.0_wp ) * dtor )
624!
625!-- calculation of solar zenith angle
626    thetar  = acos( sin( lat * dtor) * sin( declination * dtor ) + cos( wsp * dtor ) *             &
627    cos( lat * dtor ) * cos( declination * dtor ) )
628    thetasr = asin( sin( lat * dtor) * sin( declination * dtor ) + cos( wsp * dtor ) *             & 
629    cos( lat * dtor ) * cos( declination * dtor ) )
630    sza = thetar / dtor
631!
632!-- calculation of solar azimuth angle
633    IF (woz .LE. 12.0_wp) alpha = pi - acos( (sin(thetasr) * sin(lat * dtor) -                     & 
634    sin( declination * dtor )) / ( cos(thetasr) * cos( lat * dtor ) ) )   
635    IF (woz .GT. 12.0_wp) alpha = pi + acos( (sin(thetasr) * sin(lat * dtor) -                     &
636    sin( declination * dtor )) / ( cos(thetasr) * cos( lat * dtor ) ) )   
637    saa = alpha / dtor
638
639 END SUBROUTINE uvem_solar_position
640
641
642!------------------------------------------------------------------------------!
643! Description:
644! ------------
645!> Module-specific routine for new module
646!------------------------------------------------------------------------------!
647 SUBROUTINE uvem_calc_exposure
648
649    USE indices,                                                               &
650        ONLY:  nxlg, nxrg, nyng, nysg, nys, nyn, nxl, nxr 
651   
652   
653    IMPLICIT NONE   
654   
655    INTEGER(iwp) ::  i   !< running index
656    INTEGER(iwp) ::  j   !< running index
657    INTEGER(iwp) ::  k   !< running index
658
659    CALL uvem_solar_position
660     
661    IF (sza .GE. 90) THEN
662       vitd3_exposure(:,:) = 0.0_wp
663    ELSE
664       
665       DO  i = 0, 35
666          DO  j = 0, 9
667                projection_area_lookup_table(i,j) = uvem_projarea_f%var(clothing,j,i)
668          ENDDO
669       ENDDO
670       DO  i = 0, 35
671          DO  j = 0, 9
672                integration_array(i,j) = uvem_integration_f%var(j,i)
673          ENDDO
674       ENDDO
675       DO  i = 0, 2
676          DO  j = 0, 90
677                irradiance_lookup_table(i,j) = uvem_irradiance_f%var(j,i)
678          ENDDO
679       ENDDO
680       DO  i = 0, 35
681          DO  j = 0, 9
682             DO  k = 0, 90
683                radiance_lookup_table(i,j,k) = uvem_radiance_f%var(k,j,i)
684             ENDDO
685          ENDDO
686       ENDDO
687       
688       
689       
690!--    rotate 3D-Model human to desired direction  -----------------------------
691       projection_area_temp( 0:35,:) = projection_area_lookup_table
692       projection_area_temp(36:71,:) = projection_area_lookup_table               
693       IF ( .NOT. turn_to_sun )  startpos_human = orientation_angle / 10.0_wp
694       IF (       turn_to_sun )  startpos_human = saa / 10.0_wp       
695       DO  ai = 0, 35
696           xfactor = ( startpos_human ) - INT( startpos_human )
697          DO  zi = 0, 9
698              projection_area(ai,zi) = ( projection_area_temp( 36 - INT( startpos_human ) - 1 + ai , zi) * &
699                                       ( xfactor ) ) &
700                                      +( projection_area_temp( 36 - INT( startpos_human ) + ai , zi) * &
701                                        ( 1.0_wp - xfactor ) )
702          ENDDO
703       ENDDO           
704!             
705!           
706!--    interpolate to accurate Solar Zenith Angle  ------------------         
707       DO  ai = 0, 35
708          xfactor = (sza)-INT(sza)
709          DO  zi = 0, 9
710              radiance_array(ai,zi) = ( radiance_lookup_table(ai, zi, INT(sza) ) * ( 1.0_wp - xfactor) ) +&
711              ( radiance_lookup_table(ai,zi,INT(sza) + 1) * xfactor )
712          ENDDO
713       ENDDO
714       Do  ai = 0, 2
715           irradiance(ai) = ( irradiance_lookup_table(ai, INT(sza) ) * ( 1.0_wp - xfactor)) + &
716           (irradiance_lookup_table(ai, INT(sza) + 1) * xfactor )
717       ENDDO   
718!         
719!--    interpolate to accurate Solar Azimuth Angle ------------------
720       IF ( sun_in_south )  THEN
721          startpos_saa_float = 180.0_wp / 10.0_wp
722       ELSE
723          startpos_saa_float = saa / 10.0_wp
724       ENDIF
725       radiance_array_temp( 0:35,:) = radiance_array
726       radiance_array_temp(36:71,:) = radiance_array
727       xfactor = (startpos_saa_float) - INT(startpos_saa_float)
728              DO  ai = 0, 35
729          DO  zi = 0, 9
730             radiance_array(ai,zi) = ( radiance_array_temp( 36 - INT( startpos_saa_float ) - 1 + ai , zi ) * &
731                                     ( xfactor ) ) &
732                                     + ( radiance_array_temp( 36 - INT( startpos_saa_float ) + ai , zi ) &
733                                     * ( 1.0_wp - xfactor ) )
734          ENDDO
735       ENDDO 
736!       
737!     
738!--    calculate Projectionarea for direct beam -----------------------------'
739       projection_area_direct_temp( 0:35,:) = projection_area
740       projection_area_direct_temp(36:71,:) = projection_area
741       yfactor = ( sza / 10.0_wp ) - INT( sza / 10.0_wp )
742       xfactor = ( startpos_saa_float ) - INT( startpos_saa_float )
743       projection_area_direct_beam = ( projection_area_direct_temp( INT(startpos_saa_float)    ,INT(sza/10.0_wp)  ) *&
744                                     ( 1.0_wp - xfactor ) * ( 1.0_wp - yfactor ) ) + &
745                                     ( projection_area_direct_temp( INT(startpos_saa_float) + 1,INT(sza/10.0_wp)  ) *&
746                                     (          xfactor ) * ( 1.0_wp - yfactor ) ) + &
747                                     ( projection_area_direct_temp( INT(startpos_saa_float)    ,INT(sza/10.0_wp)+1) *&
748                                     ( 1.0_wp - xfactor ) * (          yfactor ) ) + &
749                                     ( projection_area_direct_temp( INT(startpos_saa_float) + 1,INT(sza/10.0_wp)+1) *&
750                                     (          xfactor ) * (          yfactor ) )
751!                                               
752!                                               
753!                                               
754       DO  i = nxl, nxr    !nxlg, nxrg
755          DO  j = nys, nyn    !nysg, nyng
756!                   
757! !--        extract obstruction from IBSET-Integer_Array ------------------'
758             IF (consider_obstructions ) THEN
759                obstruction_temp1 = building_obstruction_f%var_3d(:,j,i)
760                IF (obstruction_temp1(0) .NE. 9) THEN
761                   DO  zi = 0, 44 
762                      DO  ai = 0, 7 
763                         IF ( btest( obstruction_temp1(zi), ai ) .EQV. .TRUE.) THEN
764                            obstruction_temp2( ( zi * 8 ) + ai ) = 1
765                         ELSE
766                            obstruction_temp2( ( zi * 8 ) + ai ) = 0
767                         ENDIF
768                      ENDDO
769                   ENDDO       
770                   DO  zi = 0, 9                                         
771                      obstruction(:,zi) = obstruction_temp2( zi * 36 :( zi * 36) + 35 )
772                   ENDDO
773                ELSE
774                   obstruction(:,:) = 0
775                ENDIF
776             ENDIF
777!             
778! !--        calculated human exposure ------------------' 
779             diffuse_exposure = SUM( radiance_array * projection_area * integration_array * obstruction )     
780         
781             obstruction_direct_beam = obstruction( nint(startpos_saa_float), nint( sza / 10.0_wp ) ) 
782             IF (sza .GE. 89.99_wp) THEN
783                sza = 89.99999_wp
784             ENDIF
785!             
786! !--        calculate direct normal irradiance (direct beam) ------------------'
787             direct_exposure = ( irradiance(1) / cos( pi * sza / 180.0_wp ) ) * &
788             projection_area_direct_beam * obstruction_direct_beam 
789               
790             vitd3_exposure(j,i) = ( diffuse_exposure + direct_exposure ) / 1000.0_wp * 70.97_wp 
791             ! unit = international units vitamin D per second             
792          ENDDO
793       ENDDO
794    ENDIF
795
796 END SUBROUTINE uvem_calc_exposure
797
798 END MODULE uv_exposure_model_mod
Note: See TracBrowser for help on using the repository browser.