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

Last change on this file since 3022 was 3014, checked in by maronga, 6 years ago

series of bugfixes

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