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

Last change on this file since 2932 was 2932, checked in by maronga, 3 years ago

renamed all Fortran NAMELISTS

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