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

Last change on this file since 2699 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

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