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

Last change on this file since 2894 was 2894, checked in by Giersch, 4 years ago

Reading/Writing? data in case of restart runs revised

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