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

Last change on this file since 3246 was 3246, checked in by sward, 3 years ago

Added error handling for wrong input parameters

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