Changeset 3824 for palm


Ignore:
Timestamp:
Mar 27, 2019 3:56:16 PM (5 years ago)
Author:
pavelkrc
Message:

Code review of radiation_model_mod.f90

Location:
palm/trunk/SOURCE
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/chem_photolysis_mod.f90

    r3614 r3824  
    234234
    235235       USE radiation_model_mod,                                                   &
    236            ONLY:  calc_zenith, zenith
     236           ONLY:  calc_zenith, cos_zenith
    237237
    238238       IMPLICIT NONE
     
    247247       CALL calc_zenith
    248248
    249        IF ( zenith(0) > 0.0_wp ) THEN
    250           coszi = 1. / zenith(0)
     249       IF ( cos_zenith > 0.0_wp ) THEN
     250          coszi = 1. / cos_zenith
    251251
    252252          DO iphot = 1,nphot
     
    254254                IF ( TRIM( names_s(iav) ) == TRIM( phot_names(iphot) ) ) then
    255255                   phot_frequen(iphot)%freq(nzb+1:nzt,:,:) =    &
    256                                 par_l(iav) * zenith(0)**par_m(iav) *  EXP( -par_n(iav) * coszi )
     256                                par_l(iav) * cos_zenith**par_m(iav) *  EXP( -par_n(iav) * coszi )
    257257                ENDIF
    258258             ENDDO
  • palm/trunk/SOURCE/chemistry_model_mod.f90

    r3821 r3824  
    27992799
    28002800    USE radiation_model_mod,                                                &
    2801          ONLY:  zenith
     2801         ONLY:  cos_zenith
    28022802
    28032803
     
    34193419                !
    34203420                !-- Get Rc_tot
    3421                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
     3421                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
    34223422                     relh, lai, sai, nwet, lu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
    34233423                     r_aero_lu , Rb )
     
    35663566                !
    35673567                !-- Get Rc_tot
    3568                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
     3568                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
    35693569                     relh, lai, sai, nwet, lup_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
    35703570                     r_aero_lu , Rb )
     
    37133713                !
    37143714                !-- Get Rc_tot
    3715                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
     3715                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
    37163716                     relh, lai, sai, nwet, luw_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
    37173717                     r_aero_lu , Rb )
     
    41094109                !
    41104110                !-- Get Rc_tot
    4111                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
     4111                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
    41124112                     relh, lai, sai, nwet, luu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
    41134113                     r_aero_lu , Rb )
     
    42494249                !
    42504250                !-- Get Rc_tot
    4251                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
     4251                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
    42524252                     relh, lai, sai, nwet, lug_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
    42534253                     r_aero_lu , Rb )
     
    43954395                !
    43964396                !-- Get Rc_tot
    4397                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
     4397                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
    43984398                     relh, lai, sai, nwet, lud_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
    43994399                     r_aero_lu , Rb )
  • palm/trunk/SOURCE/radiation_model_mod.f90

    r3814 r3824  
    2828! -----------------
    2929! $Id$
     30! Change zenith(0:0) and others to scalar.
     31! Code review.
     32!
     33! 3814 2019-03-26 08:40:31Z pavelkrc
    3034! Rename exported nzu, nzp and related variables due to name conflict
    3135!
     
    747751
    748752
    749     REAL(wp), DIMENSION(0:0) ::  zenith,         & !< cosine of solar zenith angle
    750                                  sun_dir_lat,    & !< solar directional vector in latitudes
    751                                  sun_dir_lon       !< solar directional vector in longitudes
     753    REAL(wp) ::  cos_zenith        !< cosine of solar zenith angle, also z-coordinate of solar unit vector
     754    REAL(wp) ::  sun_dir_lat       !< y-coordinate of solar unit vector
     755    REAL(wp) ::  sun_dir_lon       !< x-coordinate of solar unit vector
    752756
    753757    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_net_av       !< average of net radiation (rad_net) at surface
     
    12941298           rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,  &
    12951299           rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr,              &
    1296            rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb, solar_constant, &
     1300           rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, solar_constant,          &
    12971301           skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
    1298            zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,       &
     1302           cos_zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,   &
    12991303           idir, jdir, kdir, id, iz, iy, ix,                                   &
    13001304           iup_u, inorth_u, isouth_u, ieast_u, iwest_u,                        &
     
    28342838!
    28352839!--    Calculate sky transmissivity
    2836        sky_trans = 0.6_wp + 0.2_wp * zenith(0)
     2840       sky_trans = 0.6_wp + 0.2_wp * cos_zenith
    28372841
    28382842!
     
    29092913                k = nz_urban_t
    29102914
    2911                 surf%rad_sw_in  = solar_constant * sky_trans * zenith(0)
     2915                surf%rad_sw_in  = solar_constant * sky_trans * cos_zenith
    29122916                surf%rad_sw_out = albedo_urb * surf%rad_sw_in
    29132917               
     
    29332937                   k = surf%k(m)
    29342938
    2935                    surf%rad_sw_in(m) = solar_constant * sky_trans * zenith(0)
     2939                   surf%rad_sw_in(m) = solar_constant * sky_trans * cos_zenith
    29362940
    29372941!
     
    33793383
    33803384       REAL(wp), DIMENSION(0:nzt+1) :: pt_av, q_av, ql_av
     3385       REAL(wp), DIMENSION(0:0)     :: zenith   !< to provide indexed array
    33813386!
    33823387!--    Just dummy arguments
     
    33953400!--    Calculate current (cosine of) zenith angle and whether the sun is up
    33963401       CALL calc_zenith     
     3402       zenith(0) = cos_zenith
    33973403!
    33983404!--    Calculate surface albedo. In case average radiation is applied,
     
    41724178!
    41734179!--    Calculate cosine of solar zenith angle
    4174        zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)   &
     4180       cos_zenith = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)   &
    41754181                                            * COS(hour_angle)
    4176        zenith(0) = MAX(0.0_wp,zenith(0))
     4182       cos_zenith = MAX(0.0_wp,cos_zenith)
    41774183
    41784184!
     
    41824188!
    41834189!--       Direction in longitudes equals to sin(solar_azimuth) * sin(zenith)
    4184           sun_dir_lon(0) = -SIN(hour_angle) * COS(declination)
     4190          sun_dir_lon = -SIN(hour_angle) * COS(declination)
    41854191
    41864192!
    41874193!--       Direction in latitues equals to cos(solar_azimuth) * sin(zenith)
    4188           sun_dir_lat(0) = SIN(declination) * COS(lat) - COS(hour_angle) &
     4194          sun_dir_lat = SIN(declination) * COS(lat) - COS(hour_angle) &
    41894195                              * COS(declination) * SIN(lat)
    41904196       ENDIF
     
    41924198!
    41934199!--    Check if the sun is up (otheriwse shortwave calculations can be skipped)
    4194        IF ( zenith(0) > 0.0_wp )  THEN
     4200       IF ( cos_zenith > 0.0_wp )  THEN
    41954201          sun_up = .TRUE.
    41964202       ELSE
     
    42274233                 IF ( surf%albedo_type(ind_type,m) == 1 )  THEN
    42284234                    surf%rrtm_aldir(ind_type,m) = 0.026_wp /                    &
    4229                                                 ( zenith(0)**1.7_wp + 0.065_wp )&
    4230                                      + 0.15_wp * ( zenith(0) - 0.1_wp )         &
    4231                                                * ( zenith(0) - 0.5_wp )         &
    4232                                                * ( zenith(0) - 1.0_wp )
     4235                                                ( cos_zenith**1.7_wp + 0.065_wp )&
     4236                                     + 0.15_wp * ( cos_zenith - 0.1_wp )         &
     4237                                               * ( cos_zenith - 0.5_wp )         &
     4238                                               * ( cos_zenith - 1.0_wp )
    42334239                    surf%rrtm_asdir(ind_type,m) = surf%rrtm_aldir(ind_type,m)
    42344240!
    42354241!--              Snow
    42364242                 ELSEIF ( surf%albedo_type(ind_type,m) == 16 )  THEN
    4237                     IF ( zenith(0) < 0.5_wp )  THEN
     4243                    IF ( cos_zenith < 0.5_wp )  THEN
    42384244                       surf%rrtm_aldir(ind_type,m) =                           &
    42394245                                 0.5_wp * ( 1.0_wp - surf%aldif(ind_type,m) )  &
    42404246                                        * ( 3.0_wp / ( 1.0_wp + 4.0_wp         &
    4241                                         * zenith(0) ) ) - 1.0_wp
     4247                                        * cos_zenith ) ) - 1.0_wp
    42424248                       surf%rrtm_asdir(ind_type,m) =                           &
    42434249                                 0.5_wp * ( 1.0_wp - surf%asdif(ind_type,m) )  &
    42444250                                        * ( 3.0_wp / ( 1.0_wp + 4.0_wp         &
    4245                                         * zenith(0) ) ) - 1.0_wp
     4251                                        * cos_zenith ) ) - 1.0_wp
    42464252
    42474253                       surf%rrtm_aldir(ind_type,m) =                           &
     
    42824288                          surf%rrtm_aldir(ind_type,m) =                        &
    42834289                                surf%aldif(ind_type,m) * 1.4_wp /              &
    4284                                            ( 1.0_wp + 0.8_wp * zenith(0) )
     4290                                           ( 1.0_wp + 0.8_wp * cos_zenith )
    42854291                          surf%rrtm_asdir(ind_type,m) =                        &
    42864292                                surf%asdif(ind_type,m) * 1.4_wp /              &
    4287                                            ( 1.0_wp + 0.8_wp * zenith(0) )
     4293                                           ( 1.0_wp + 0.8_wp * cos_zenith )
    42884294!
    42894295!--                    Surface types with weak zenith dependence
     
    42914297                          surf%rrtm_aldir(ind_type,m) =                        &
    42924298                                surf%aldif(ind_type,m) * 1.1_wp /              &
    4293                                            ( 1.0_wp + 0.2_wp * zenith(0) )
     4299                                           ( 1.0_wp + 0.2_wp * cos_zenith )
    42944300                          surf%rrtm_asdir(ind_type,m) =                        &
    42954301                                surf%asdif(ind_type,m) * 1.1_wp /              &
    4296                                            ( 1.0_wp + 0.2_wp * zenith(0) )
     4302                                           ( 1.0_wp + 0.2_wp * cos_zenith )
    42974303
    42984304                       CASE DEFAULT
     
    50185024     mrot(2, :) = (/ 0._wp,  COS(alpha), SIN(alpha) /)
    50195025     mrot(3, :) = (/ 0._wp, -SIN(alpha), COS(alpha) /)
    5020      sunorig = (/ zenith(0), sun_dir_lat, sun_dir_lon /)
     5026     sunorig = (/ cos_zenith, sun_dir_lat, sun_dir_lon /)
    50215027     sunorig = MATMUL(mrot, sunorig)
    50225028     DO d = 0, nsurf_type
     
    50245030     ENDDO
    50255031
    5026      IF ( zenith(0) > 0 )  THEN
     5032     IF ( cos_zenith > 0 )  THEN
    50275033!--      now we will "squash" the sunorig vector by grid box size in
    50285034!--      each dimension, so that this new direction vector will allow us
     
    51865192
    51875193     !-- direct radiation
    5188      IF ( zenith(0) > 0 )  THEN
     5194     IF ( cos_zenith > 0 )  THEN
    51895195        !--Identify solar direction vector (discretized number) 1)
    51905196        !--
    5191         j = FLOOR(ACOS(zenith(0)) / pi * raytrace_discrete_elevs)
    5192         i = MODULO(NINT(ATAN2(sun_dir_lon(0), sun_dir_lat(0))               &
     5197        j = FLOOR(ACOS(cos_zenith) / pi * raytrace_discrete_elevs)
     5198        i = MODULO(NINT(ATAN2(sun_dir_lon, sun_dir_lat)               &
    51935199                        / (2._wp*pi) * raytrace_discrete_azims-.5_wp, iwp), &
    51945200                   raytrace_discrete_azims)
     
    51995205           i = surfl(ix, isurf)
    52005206           surfinswdir(isurf) = rad_sw_in_dir(j,i) *                        &
    5201                 costheta(surfl(id, isurf)) * dsitrans(isurf, isd) / zenith(0)
     5207                costheta(surfl(id, isurf)) * dsitrans(isurf, isd) / cos_zenith
    52025208        ENDDO
    52035209!
     
    52075213           i = mrtbl(ix, imrt)
    52085214           mrtinsw(imrt) = mrtinsw(imrt) + mrtdsit(imrt, isd) * rad_sw_in_dir(j,i) &
    5209                                      / zenith(0) / 4._wp ! normal to sphere
     5215                                     / cos_zenith / 4._wp ! normal to sphere
    52105216        ENDDO
    52115217     ENDIF
     
    52175223        mrtinlw(imrt) = mrtinlw(imrt) + mrtf(imrtf) * surfoutl(isurfsrc)
    52185224     ENDDO
    5219 
     5225!
     5226!--  Absorption in each local plant canopy grid box from the first atmospheric
     5227!--  pass of radiation
    52205228     IF ( npcbl > 0 )  THEN
    52215229
     
    52235231         pcbinswdif(:) = 0._wp
    52245232         pcbinlw(:) = 0._wp
    5225 !
    5226 !--      pcsf first pass
     5233
    52275234         DO icsf = 1, ncsfl
    52285235             ipcgb = csfsurf(1, icsf)
     
    52345241             IF ( isurfsrc == -1 )  THEN
    52355242!
    5236 !--             Diffuse rad from sky.
     5243!--             Diffuse radiation from sky
    52375244                pcbinswdif(ipcgb) = csf(1,icsf) * rad_sw_in_diff(j,i)
    52385245!
    5239 !--             Absorbed diffuse LW from sky minus emitted to sky
     5246!--             Absorbed diffuse LW radiation from sky minus emitted to sky
    52405247                IF ( plant_lw_interact )  THEN
    52415248                   pcbinlw(ipcgb) = csf(1,icsf)                                  &
     
    52445251                ENDIF
    52455252!
    5246 !--             Direct rad
    5247                 IF ( zenith(0) > 0 )  THEN
     5253!--             Direct solar radiation
     5254                IF ( cos_zenith > 0 )  THEN
    52485255!--                Estimate directed box absorption
    52495256                   pc_abs_frac = 1._wp - exp(pc_abs_eff * lad_s(k,j,i))
     
    53125319
    53135320         surfoutsl = albedo_surf * surfins
    5314 !--         for non-transparent surfaces, longwave albedo is 1 - emissivity
     5321!
     5322!--      for non-transparent surfaces, longwave albedo is 1 - emissivity
    53155323         surfoutll = (1._wp - emiss_surf) * surfinl
    53165324
     
    53185326         CALL MPI_AllGatherv(surfoutsl, nsurfl, MPI_REAL, &
    53195327             surfouts, nsurfs, surfstart, MPI_REAL, comm2d, ierr)
    5320          IF ( ierr /= 0 ) THEN
     5328         IF ( ierr /= 0 )  THEN
    53215329             WRITE(9,*) 'Error MPI_AllGatherv2:', ierr, SIZE(surfoutsl), nsurfl, &
    53225330                        SIZE(surfouts), nsurfs, surfstart
     
    53265334         CALL MPI_AllGatherv(surfoutll, nsurfl, MPI_REAL, &
    53275335             surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr)
    5328          IF ( ierr /= 0 ) THEN
     5336         IF ( ierr /= 0 )  THEN
    53295337             WRITE(9,*) 'Error MPI_AllGatherv3:', ierr, SIZE(surfoutll), nsurfl, &
    53305338                        SIZE(surfoutl), nsurfs, surfstart
     
    53365344         surfoutl = surfoutll
    53375345#endif
    5338 
    5339 !--         reset for next pass input
     5346!
     5347!--      Reset for the input from next reflective pass
    53405348         surfins = 0._wp
    53415349         surfinl = 0._wp
    5342 
    5343 !--         reflected radiation
     5350!
     5351!--      Reflected radiation
    53445352         DO isvf = 1, nsvfl
    53455353             isurf = svfsurf(1, isvf)
     
    58255833!> This subroutine splits direct and diffusion dw radiation
    58265834!> It sould not be called in case the radiation model already does it
    5827 !> It follows <CITATION>
     5835!> It follows Boland, Ridley & Brown (2008)
    58285836!------------------------------------------------------------------------------!
    58295837    SUBROUTINE calc_diffusion_radiation
     
    58555863!--     towards 0 while keeping full continuity.
    58565864!--   
    5857         IF ( zenith(0) <= lowest_solarUp )  THEN
     5865        IF ( cos_zenith <= lowest_solarUp )  THEN
    58585866            corrected_solarUp = lowest_solarUp
    58595867        ELSE
    5860             corrected_solarUp = zenith(0)
     5868            corrected_solarUp = cos_zenith
    58615869        ENDIF
    58625870       
     
    82258233!--      Update apparent solar position based on modified t_s_r_p
    82268234         CALL calc_zenith
    8227          IF ( zenith(0) > 0 )  THEN
     8235         IF ( cos_zenith > 0 )  THEN
    82288236!--         
    82298237!--         Identify solar direction vector (discretized number) 1)
    8230             i = MODULO(NINT(ATAN2(sun_dir_lon(0), sun_dir_lat(0))               &
     8238            i = MODULO(NINT(ATAN2(sun_dir_lon, sun_dir_lat)               &
    82318239                            / (2._wp*pi) * raytrace_discrete_azims-.5_wp, iwp), &
    82328240                       raytrace_discrete_azims)
    8233             j = FLOOR(ACOS(zenith(0)) / pi * raytrace_discrete_elevs)
     8241            j = FLOOR(ACOS(cos_zenith) / pi * raytrace_discrete_elevs)
    82348242            IF ( dsidir_rev(j, i) == -1 )  THEN
    82358243               ndsidir = ndsidir + 1
  • palm/trunk/SOURCE/urban_surface_mod.f90

    r3814 r3824  
    2828! -----------------
    2929! $Id$
     30! Remove unused imports
     31!
     32!
     33! 3814 2019-03-26 08:40:31Z pavelkrc
    3034! unused subroutine commented out
    3135!
     
    438442
    439443    USE basic_constants_and_equations_mod,                                     &
    440         ONLY:  c_p, g, kappa, pi, r_d, rho_l, l_v
     444        ONLY:  c_p, g, kappa, pi, r_d, rho_l, l_v, sigma_sb
    441445
    442446    USE control_parameters,                                                    &
     
    473477       
    474478    USE radiation_model_mod,                                                   &
    475         ONLY:  albedo_type, radiation_interaction, calc_zenith, zenith,        &
     479        ONLY:  albedo_type, radiation_interaction,                             &
    476480               radiation, rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out,        &
    477                sigma_sb, sun_direction, sun_dir_lat, sun_dir_lon,              &
    478481               force_radiation_call, iup_u, inorth_u, isouth_u, ieast_u,       &
    479482               iwest_u, iup_l, inorth_l, isouth_l, ieast_l, iwest_l, id,       &
Note: See TracChangeset for help on using the changeset viewer.