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

Last change on this file since 3004 was 3004, checked in by Giersch, 3 years ago

precipitation_rate removed, further allocation checks for data output of averaged quantities implemented, double CALL of flow_statistics at the beginning of time_integration removed, further minor bugfixes, comments added

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