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

Last change on this file since 3448 was 3448, checked in by kanani, 4 years ago

Implementation of human thermal indices (from branch biomet_p2 at r3444)

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