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

Last change on this file since 3253 was 3248, checked in by sward, 6 years ago

Minor format changes

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