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

Last change on this file since 3381 was 3274, checked in by knoop, 6 years ago

Modularization of all bulk cloud physics code components

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