Changeset 3272 for palm


Ignore:
Timestamp:
Sep 24, 2018 10:16:32 AM (6 years ago)
Author:
suehring
Message:

Merge with branch radition - improved representation of diffuse radiation

Location:
palm/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/LIB/rrtmg/rrtmg_sw_rad.nomcica.f90

    r1585 r3272  
    7878      subroutine rrtmg_sw &
    7979            (ncol    ,nlay    ,icld    ,iaer    , &
    80              play    ,plev    ,tlay    ,tlev    ,tsfc    , &
    81              h2ovmr  ,o3vmr   ,co2vmr  ,ch4vmr  ,n2ovmr  ,o2vmr, &
    82              asdir   ,asdif   ,aldir   ,aldif   , &
    83              coszen  ,adjes   ,dyofyr  ,scon    , &
    84              inflgsw ,iceflgsw,liqflgsw,cldfr   , &
    85              taucld  ,ssacld  ,asmcld  ,fsfcld  , &
    86              cicewp  ,cliqwp  ,reice   ,reliq   , &
    87              tauaer  ,ssaaer  ,asmaer  ,ecaer   , &
    88              swuflx  ,swdflx  ,swhr    ,swuflxc ,swdflxc ,swhrc)
     80             play    ,plev    ,tlay    ,tlev    , &
     81             tsfc    ,h2ovmr  ,o3vmr   ,co2vmr  , &
     82             ch4vmr  ,n2ovmr  ,o2vmr   ,asdir   , &
     83             asdif   ,aldir   ,aldif   ,coszen  , &
     84             adjes   ,dyofyr  ,scon    ,inflgsw , &
     85             iceflgsw,liqflgsw,cldfr   ,taucld  , &
     86             ssacld  ,asmcld  ,fsfcld  ,cicewp  , &
     87             cliqwp  ,reice   ,reliq   ,tauaer  , &
     88             ssaaer  ,asmaer  ,ecaer   ,swuflx  , &
     89             swdflx  ,swhr    ,swuflxc ,swdflxc , &
     90             swhrc   ,dirdflux,difdflux)
    8991
    9092! ------- Description -------
     
    169171!   delta scaling based on setting of idelm flag.
    170172!     Dec 2008: M. J. Iacono, AER, Inc.
     173!-- Added the output direct and diffuse fluxes to the argument of rrtmg_sw
     174!     Aug 2018: M. Salim, Humboldt Uni, Berlin, Germany
    171175
    172176! --------- Modules ---------
     
    417421      real(kind=rb) :: swnflx(nlay+2)         ! Total sky shortwave net flux (W/m2)
    418422      real(kind=rb) :: swnflxc(nlay+2)        ! Clear sky shortwave net flux (W/m2)
    419       real(kind=rb) :: dirdflux(nlay+2)       ! Direct downward shortwave surface flux
    420       real(kind=rb) :: difdflux(nlay+2)       ! Diffuse downward shortwave surface flux
     423      real(kind=rb), intent(out) :: dirdflux(nlay+2)       ! Direct downward shortwave surface flux
     424      real(kind=rb), intent(out) :: difdflux(nlay+2)       ! Diffuse downward shortwave surface flux
    421425      real(kind=rb) :: uvdflx(nlay+2)         ! Total sky downward shortwave flux, UV/vis 
    422426      real(kind=rb) :: nidflx(nlay+2)         ! Total sky downward shortwave flux, near-IR
  • palm/trunk/SOURCE/radiation_model_mod.f90

    r3261 r3272  
    2828! -----------------
    2929! $Id$
     30! - split direct and diffusion shortwave radiation using RRTMG rather than using
     31!   calc_diffusion_radiation, in case of RRTMG
     32! - removed the namelist variable split_diffusion_radiation. Now splitting depends
     33!   on the choise of radiation radiation scheme
     34! - removed calculating the rdiation flux for surfaces at the radiation scheme
     35!   in case of using RTM since it will be calculated anyway in the radiation
     36!   interaction routine.
     37! - set SW radiation flux for surfaces to zero at night in case of no RTM is used
     38! - precalculate the unit vector yxdir of ray direction to avoid the temporarly
     39!   array allocation during the subroutine call
     40! - fixed a bug in calculating the max number of boxes ray can cross in the domain
     41!
     42! 3264 2018-09-20 13:54:11Z moh.hefny
    3043! Bugfix in raytrace_2d calls
    3144!
     
    722735                                             rrtm_swuflxc,   & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
    723736                                             rrtm_swhr,      & !< RRTM output of shortwave radiation heating rate (K/d)
    724                                              rrtm_swhrc        !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
    725 
     737                                             rrtm_swhrc,     & !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
     738                                             rrtm_dirdflux,  & !< RRTM output of incoming direct shortwave (W/m²)
     739                                             rrtm_difdflux     !< RRTM output of incoming diffuse shortwave (W/m²)
    726740
    727741    REAL(wp), DIMENSION(1) ::                rrtm_aldif,     & !< surface albedo for longwave diffuse radiation
     
    822836
    823837!-- configuration parameters (they can be setup in PALM config)
    824     LOGICAL                                        ::  split_diffusion_radiation = .TRUE. !< split direct and diffusion dw radiation
    825                                                                                           !< (.F. in case the radiation model already does it)   
    826838    LOGICAL                                        ::  rma_lad_raytrace = .FALSE.         !< use MPI RMA to access LAD for raytracing (instead of global array)
    827839    LOGICAL                                        ::  mrt_factors = .FALSE.              !< whether to generate MRT factor files during init
     
    10801092           skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
    10811093           zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,       &
    1082            split_diffusion_radiation,                                          &
    10831094           nrefsteps, mrt_factors, dist_max_svf, nsvfl, svf,                   &
    10841095           svfsurf, surfinsw, surfinlw, surfins, surfinl, surfinswdir,         &
     
    11531164       SELECT CASE ( TRIM( var ) )
    11541165
    1155          CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr' )
     1166         CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in' )
    11561167             IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
    11571168                message_string = '"output of "' // TRIM( var ) // '" requi' // &
     
    11891200             IF ( TRIM( var ) == 'rad_sw_in*'    ) unit = 'W/m2'
    11901201             IF ( TRIM( var ) == 'rad_sw_out*'   ) unit = 'W/m2'
     1202             IF ( TRIM( var ) == 'rad_sw_in'     ) unit = 'W/m2'
    11911203             IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''   
    11921204             IF ( TRIM( var ) == 'rrtm_aldir*'   ) unit = ''
     
    27692781                                  radiation_scheme, skip_time_do_radiation,    &
    27702782                                  sw_radiation, unscheduled_radiation_calls,   &
    2771                                   split_diffusion_radiation,                   &
    27722783                                  max_raytracing_dist, min_irrf_value,         &
    27732784                                  nrefsteps, mrt_factors, rma_lad_raytrace,    &
     
    27822793                                  radiation_scheme, skip_time_do_radiation,    &
    27832794                                  sw_radiation, unscheduled_radiation_calls,   &
    2784                                   split_diffusion_radiation,                   &
    27852795                                  max_raytracing_dist, min_irrf_value,         &
    27862796                                  nrefsteps, mrt_factors, rma_lad_raytrace,    &
     
    30543064                rad_lw_out(k,:,:) = rrtm_lwuflx(0,k)
    30553065             ENDDO
    3056 
     3066             rad_lw_in_diff(:,:) = rad_lw_in(0,:,:)
    30573067!
    30583068!--          Save heating rates (convert from K/d to K/h).
     
    30703080             ENDDO
    30713081
    3072 !
    3073 !--          Save surface radiative fluxes and change in LW heating rate
    3074 !--          onto respective surface elements
    3075 !--          Horizontal surfaces
    3076              IF ( surf_lsm_h%ns > 0 )  THEN
    3077                 surf_lsm_h%rad_lw_in           = rrtm_lwdflx(0,nzb)
    3078                 surf_lsm_h%rad_lw_out          = rrtm_lwuflx(0,nzb)
    3079                 surf_lsm_h%rad_lw_out_change_0 = rrtm_lwuflx_dt(0,nzb)
    3080              ENDIF             
    3081              IF ( surf_usm_h%ns > 0 )  THEN
    3082                 surf_usm_h%rad_lw_in           = rrtm_lwdflx(0,nzb)
    3083                 surf_usm_h%rad_lw_out          = rrtm_lwuflx(0,nzb)
    3084                 surf_usm_h%rad_lw_out_change_0 = rrtm_lwuflx_dt(0,nzb)
    3085              ENDIF
    3086 !
    3087 !--          Vertical surfaces.
    3088              DO  l = 0, 3
    3089                 IF ( surf_lsm_v(l)%ns > 0 )  THEN
    3090                    surf_lsm_v(l)%rad_lw_in           = rrtm_lwdflx(0,nzb)
    3091                    surf_lsm_v(l)%rad_lw_out          = rrtm_lwuflx(0,nzb)
    3092                    surf_lsm_v(l)%rad_lw_out_change_0 = rrtm_lwuflx_dt(0,nzb)
    3093                 ENDIF
    3094                 IF ( surf_usm_v(l)%ns > 0 )  THEN
    3095                    surf_usm_v(l)%rad_lw_in           = rrtm_lwdflx(0,nzb)
    3096                    surf_usm_v(l)%rad_lw_out          = rrtm_lwuflx(0,nzb)
    3097                    surf_usm_v(l)%rad_lw_out_change_0 = rrtm_lwuflx_dt(0,nzb)
    3098                 ENDIF
    3099              ENDDO
    3100 
    31013082          ENDIF
    31023083
    31033084          IF ( sw_radiation .AND. sun_up )  THEN
    3104              CALL rrtmg_sw( 1, nzt_rad      , rrtm_icld  , rrtm_iaer        ,&
    3105              rrtm_play       , rrtm_plev    , rrtm_tlay  , rrtm_tlev        ,&
    3106              rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr , rrtm_co2vmr      ,&
    3107              rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr , rrtm_asdir       ,&
    3108              rrtm_asdif      , rrtm_aldir   , rrtm_aldif , zenith,           &
    3109              0.0_wp          , day_of_year  , solar_constant,   rrtm_inflgsw,&
    3110              rrtm_iceflgsw   , rrtm_liqflgsw, rrtm_cldfr , rrtm_sw_taucld  ,&
    3111              rrtm_sw_ssacld  , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp  ,&
    3112              rrtm_cliqwp     , rrtm_reice   , rrtm_reliq , rrtm_sw_tauaer  ,&
    3113              rrtm_sw_ssaaer  , rrtm_sw_asmaer  , rrtm_sw_ecaer ,             &
    3114              rrtm_swuflx     , rrtm_swdflx  , rrtm_swhr  ,                   &
    3115              rrtm_swuflxc    , rrtm_swdflxc , rrtm_swhrc )
     3085             CALL rrtmg_sw( 1, nzt_rad      , rrtm_icld     , rrtm_iaer      ,&
     3086             rrtm_play      , rrtm_plev     , rrtm_tlay     , rrtm_tlev      ,&
     3087             rrtm_tsfc      , rrtm_h2ovmr   , rrtm_o3vmr    , rrtm_co2vmr    ,&
     3088             rrtm_ch4vmr    , rrtm_n2ovmr   , rrtm_o2vmr    , rrtm_asdir     ,&
     3089             rrtm_asdif     , rrtm_aldir    , rrtm_aldif    , zenith         ,&
     3090             0.0_wp         , day_of_year   , solar_constant, rrtm_inflgsw   ,&
     3091             rrtm_iceflgsw  , rrtm_liqflgsw , rrtm_cldfr    , rrtm_sw_taucld ,&
     3092             rrtm_sw_ssacld , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp    ,&
     3093             rrtm_cliqwp    , rrtm_reice    , rrtm_reliq    , rrtm_sw_tauaer ,&
     3094             rrtm_sw_ssaaer , rrtm_sw_asmaer, rrtm_sw_ecaer , rrtm_swuflx    ,&
     3095             rrtm_swdflx    , rrtm_swhr     , rrtm_swuflxc  , rrtm_swdflxc   ,&
     3096             rrtm_swhrc     , rrtm_dirdflux , rrtm_difdflux )
    31163097 
    31173098!
    3118 !--          Save fluxes
     3099!--          Save fluxes:
     3100!--          - whole domain
    31193101             DO k = nzb, nzt+1
    31203102                rad_sw_in(k,:,:)  = rrtm_swdflx(0,k)
    31213103                rad_sw_out(k,:,:) = rrtm_swuflx(0,k)
    31223104             ENDDO
     3105!--          - direct and diffuse SW at urban-surface-layer (required by RTM)
     3106             rad_sw_in_dir(:,:) = rrtm_dirdflux(0,nzb)
     3107             rad_sw_in_diff(:,:) = rrtm_difdflux(0,nzb)
    31233108
    31243109!
     
    31283113                rad_sw_cs_hr(k,:,:)  = rrtm_swhrc(0,k) * d_hours_day
    31293114             ENDDO
    3130 
    3131 !
    3132 !--          Save surface radiative fluxes onto respective surface elements
    3133 !--          Horizontal surfaces
    3134              IF ( surf_lsm_h%ns > 0 )  THEN
    3135                    surf_lsm_h%rad_sw_in     = rrtm_swdflx(0,nzb)
    3136                    surf_lsm_h%rad_sw_out    = rrtm_swuflx(0,nzb)
    3137              ENDIF
    3138              IF ( surf_usm_h%ns > 0 )  THEN
    3139                    surf_usm_h%rad_sw_in     = rrtm_swdflx(0,nzb)
    3140                    surf_usm_h%rad_sw_out    = rrtm_swuflx(0,nzb)
    3141              ENDIF
    3142 !
    3143 !--          Vertical surfaces. Fluxes are obtain at respective vertical
    3144 !--          level of the surface element
    3145              DO  l = 0, 3
    3146                 IF ( surf_lsm_v(l)%ns > 0 )  THEN
    3147                       surf_lsm_v(l)%rad_sw_in  = rrtm_swdflx(0,nzb)
    3148                       surf_lsm_v(l)%rad_sw_out = rrtm_swuflx(0,nzb)
    3149                 ENDIF             
    3150                 IF ( surf_usm_v(l)%ns > 0 )  THEN
    3151                       surf_usm_v(l)%rad_sw_in  = rrtm_swdflx(0,nzb)
    3152                       surf_usm_v(l)%rad_sw_out = rrtm_swuflx(0,nzb)
    3153                 ENDIF       
    3154              ENDDO
    31553115!
    31563116!--       Solar radiation is zero during night
     
    31583118             rad_sw_in  = 0.0_wp
    31593119             rad_sw_out = 0.0_wp
     3120             rad_sw_in_dir(:,:) = 0.0_wp
     3121             rad_sw_in_diff(:,:) = 0.0_wp
    31603122          ENDIF
    31613123!
    31623124!--    RRTMG is called for each (j,i) grid point separately, starting at the
    3163 !--    highest topography level
     3125!--    highest topography level. Here no RTM is used since average_radiation is false
    31643126       ELSE
    31653127!
     
    35513513                                  rrtm_swuflxc(:,k_topo:nzt_rad+1),            &
    35523514                                  rrtm_swdflxc(:,k_topo:nzt_rad+1),            &
    3553                                   rrtm_swhrc(:,k_topo+1:nzt_rad+1) )
     3515                                  rrtm_swhrc(:,k_topo+1:nzt_rad+1),            &
     3516                                  rrtm_dirdflux(:,k_topo:nzt_rad+1),           &
     3517                                  rrtm_difdflux(:,k_topo:nzt_rad+1) )
    35543518
    35553519                   DEALLOCATE( rrtm_sw_taucld_dum )
     
    36093573                   rad_sw_in  = 0.0_wp
    36103574                   rad_sw_out = 0.0_wp
     3575!--             !!!!!!!! ATTENSION !!!!!!!!!!!!!!!
     3576!--             Surface radiative fluxes should be also set to zero here                 
     3577!--                Save surface radiative fluxes onto respective surface elements
     3578!--                Horizontal surfaces
     3579                   DO  m = surf_lsm_h%start_index(j,i),                        &
     3580                           surf_lsm_h%end_index(j,i)
     3581                      surf_lsm_h%rad_sw_in(m)     = 0.0_wp
     3582                      surf_lsm_h%rad_sw_out(m)    = 0.0_wp
     3583                   ENDDO             
     3584                   DO  m = surf_usm_h%start_index(j,i),                        &
     3585                           surf_usm_h%end_index(j,i)
     3586                      surf_usm_h%rad_sw_in(m)     = 0.0_wp
     3587                      surf_usm_h%rad_sw_out(m)    = 0.0_wp
     3588                   ENDDO
     3589!
     3590!--                Vertical surfaces. Fluxes are obtain at respective vertical
     3591!--                level of the surface element
     3592                   DO  l = 0, 3
     3593                      DO  m = surf_lsm_v(l)%start_index(j,i),                  &
     3594                              surf_lsm_v(l)%end_index(j,i)
     3595                         k                           = surf_lsm_v(l)%k(m)
     3596                         surf_lsm_v(l)%rad_sw_in(m)  = 0.0_wp
     3597                         surf_lsm_v(l)%rad_sw_out(m) = 0.0_wp
     3598                      ENDDO             
     3599                      DO  m = surf_usm_v(l)%start_index(j,i),                  &
     3600                              surf_usm_v(l)%end_index(j,i)
     3601                         k                           = surf_usm_v(l)%k(m)
     3602                         surf_usm_v(l)%rad_sw_in(m)  = 0.0_wp
     3603                         surf_usm_v(l)%rad_sw_out(m) = 0.0_wp
     3604                      ENDDO
     3605                   ENDDO
    36113606                ENDIF
    36123607
     
    36173612!
    36183613!--    Finally, calculate surface net radiation for surface elements.
    3619 !--    First, for horizontal surfaces   
    3620        DO  m = 1, surf_lsm_h%ns
    3621           surf_lsm_h%rad_net(m) = surf_lsm_h%rad_sw_in(m)                      &
    3622                                 - surf_lsm_h%rad_sw_out(m)                     &
    3623                                 + surf_lsm_h%rad_lw_in(m)                      &
    3624                                 - surf_lsm_h%rad_lw_out(m)
    3625        ENDDO
    3626        DO  m = 1, surf_usm_h%ns
    3627           surf_usm_h%rad_net(m) = surf_usm_h%rad_sw_in(m)                      &
    3628                                 - surf_usm_h%rad_sw_out(m)                     &
    3629                                 + surf_usm_h%rad_lw_in(m)                      &
    3630                                 - surf_usm_h%rad_lw_out(m)
    3631        ENDDO
    3632 !
    3633 !--    Vertical surfaces.
    3634 !--    Todo: weight with azimuth and zenith angle according to their orientation!
    3635        DO  l = 0, 3     
    3636           DO  m = 1, surf_lsm_v(l)%ns
    3637              surf_lsm_v(l)%rad_net(m) = surf_lsm_v(l)%rad_sw_in(m)             &
    3638                                       - surf_lsm_v(l)%rad_sw_out(m)            &
    3639                                       + surf_lsm_v(l)%rad_lw_in(m)             &
    3640                                       - surf_lsm_v(l)%rad_lw_out(m)
     3614       IF (  .NOT.  radiation_interactions  ) THEN
     3615!--       First, for horizontal surfaces   
     3616          DO  m = 1, surf_lsm_h%ns
     3617             surf_lsm_h%rad_net(m) = surf_lsm_h%rad_sw_in(m)                   &
     3618                                   - surf_lsm_h%rad_sw_out(m)                  &
     3619                                   + surf_lsm_h%rad_lw_in(m)                   &
     3620                                   - surf_lsm_h%rad_lw_out(m)
    36413621          ENDDO
    3642           DO  m = 1, surf_usm_v(l)%ns
    3643              surf_usm_v(l)%rad_net(m) = surf_usm_v(l)%rad_sw_in(m)             &
    3644                                       - surf_usm_v(l)%rad_sw_out(m)            &
    3645                                       + surf_usm_v(l)%rad_lw_in(m)             &
    3646                                       - surf_usm_v(l)%rad_lw_out(m)
     3622          DO  m = 1, surf_usm_h%ns
     3623             surf_usm_h%rad_net(m) = surf_usm_h%rad_sw_in(m)                   &
     3624                                   - surf_usm_h%rad_sw_out(m)                  &
     3625                                   + surf_usm_h%rad_lw_in(m)                   &
     3626                                   - surf_usm_h%rad_lw_out(m)
    36473627          ENDDO
    3648        ENDDO
     3628!
     3629!--       Vertical surfaces.
     3630!--       Todo: weight with azimuth and zenith angle according to their orientation!
     3631          DO  l = 0, 3     
     3632             DO  m = 1, surf_lsm_v(l)%ns
     3633                surf_lsm_v(l)%rad_net(m) = surf_lsm_v(l)%rad_sw_in(m)          &
     3634                                         - surf_lsm_v(l)%rad_sw_out(m)         &
     3635                                         + surf_lsm_v(l)%rad_lw_in(m)          &
     3636                                         - surf_lsm_v(l)%rad_lw_out(m)
     3637             ENDDO
     3638             DO  m = 1, surf_usm_v(l)%ns
     3639                surf_usm_v(l)%rad_net(m) = surf_usm_v(l)%rad_sw_in(m)          &
     3640                                         - surf_usm_v(l)%rad_sw_out(m)         &
     3641                                         + surf_usm_v(l)%rad_lw_in(m)          &
     3642                                         - surf_usm_v(l)%rad_lw_out(m)
     3643             ENDDO
     3644          ENDDO
     3645       ENDIF
    36493646
    36503647
     
    39053902          DEALLOCATE ( rrtm_swhr  )
    39063903          DEALLOCATE ( rrtm_swhrc )
     3904          DEALLOCATE ( rrtm_dirdflux )
     3905          DEALLOCATE ( rrtm_difdflux )
    39073906
    39083907       ENDIF
     
    40704069       ALLOCATE ( rrtm_swdflxc(0:0,nzb:nzt_rad+1) )
    40714070       ALLOCATE ( rrtm_swhrc(0:0,nzb+1:nzt_rad+1) )
     4071       ALLOCATE ( rrtm_dirdflux(0:0,nzb:nzt_rad+1) )
     4072       ALLOCATE ( rrtm_difdflux(0:0,nzb:nzt_rad+1) )
    40724073
    40734074       rrtm_swdflx  = 0.0_wp
     
    40774078       rrtm_swdflxc = 0.0_wp
    40784079       rrtm_swhrc   = 0.0_wp
     4080       rrtm_dirdflux = 0.0_wp
     4081       rrtm_difdflux = 0.0_wp
    40794082
    40804083       ALLOCATE ( rrtm_lwdflx(0:0,nzb:nzt_rad+1)  )
     
    44824485
    44834486     IMPLICIT NONE
    4484      INTEGER(iwp)                      :: i, j, k, kk, d, refstep, m, mm, l, ll
     4487
     4488     INTEGER(iwp)                      :: i, j, k, kk, is, js, d, ku, refstep, m, mm, l, ll
    44854489     INTEGER(iwp)                      :: isurf, isurfsrc, isvf, icsf, ipcgb
    44864490     INTEGER(iwp)                      :: isd                !< solar direction number
     
    45194523     REAL(wp)                          :: area_hor           !< total horizontal area of domain in all processor
    45204524
    4521 
    45224525#if ! defined( __nopointer )
    45234526     IF ( plant_canopy )  THEN
     
    45454548
    45464549     IF ( zenith(0) > 0 )  THEN
    4547 !--         now we will "squash" the sunorig vector by grid box size in
    4548 !--         each dimension, so that this new direction vector will allow us
    4549 !--         to traverse the ray path within grid coordinates directly
     4550!--      now we will "squash" the sunorig vector by grid box size in
     4551!--      each dimension, so that this new direction vector will allow us
     4552!--      to traverse the ray path within grid coordinates directly
    45504553         sunorig_grid = (/ sunorig(1)/dz(1), sunorig(2)/dy, sunorig(3)/dx /)
    4551 !--         sunorig_grid = sunorig_grid / norm2(sunorig_grid)
     4554!--      sunorig_grid = sunorig_grid / norm2(sunorig_grid)
    45524555         sunorig_grid = sunorig_grid / SQRT(SUM(sunorig_grid**2))
    45534556
    45544557         IF ( npcbl > 0 )  THEN
    4555 !--            precompute effective box depth with prototype Leaf Area Density
     4558!--         precompute effective box depth with prototype Leaf Area Density
    45564559            pc_box_dimshift = MAXLOC(ABS(sunorig), 1) - 1
    45574560            CALL box_absorb(CSHIFT((/dz(1),dy,dx/), pc_box_dimshift),      &
     
    45644567     ENDIF
    45654568
    4566 !--     split diffusion and direct part of the solar downward radiation
    4567 !--     comming from radiation model and store it in 2D arrays
    4568 !--     rad_sw_in_diff, rad_sw_in_dir and rad_lw_in_diff
    4569      IF ( split_diffusion_radiation )  THEN
    4570          CALL calc_diffusion_radiation
    4571      ELSE
    4572          rad_sw_in_diff = 0.0_wp
    4573          rad_sw_in_dir(:,:)  = rad_sw_in(0,:,:)
    4574          rad_lw_in_diff(:,:) = rad_lw_in(0,:,:)
    4575      ENDIF
     4569!--  if radiation scheme is not RRTMG, split diffusion and direct part of the solar downward radiation
     4570!--  comming from radiation model and store it in 2D arrays
     4571     IF (  radiation_scheme /= 'rrtmg'  )  CALL calc_diffusion_radiation
    45764572
    45774573!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     
    66966692         !--
    66976693         maxboxes = (ntrack + MAX(origin(1) - nzub, nzut - origin(1))) * SIZE(zdirs, 1)
    6698 !          IF ( ncsfl + maxboxes > ncsfla )  THEN
    6699 ! !--         use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)
    6700 ! !--         k = CEILING(grow_factor ** real(CEILING(log(real(ncsfl + maxboxes, kind=wp)) &
    6701 ! !--                                                / log(grow_factor)), kind=wp))
    6702 ! !--         or use this code to simply always keep some extra space after growing
    6703 !             k = CEILING(REAL(ncsfl + maxboxes, kind=wp) * grow_factor)
    6704 !             CALL merge_and_grow_csf(k)
    6705 !          ENDIF
     6694         IF ( ncsfl + maxboxes > ncsfla )  THEN
     6695!--         use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)
     6696!--         k = CEILING(grow_factor ** real(CEILING(log(real(ncsfl + maxboxes, kind=wp)) &
     6697!--                                                / log(grow_factor)), kind=wp))
     6698!--         or use this code to simply always keep some extra space after growing
     6699            k = CEILING(REAL(ncsfl + maxboxes, kind=wp) * grow_factor)
     6700            CALL merge_and_grow_csf(k)
     6701         ENDIF
    67066702
    67076703         !--Calculate transparencies and store new CSFs
     
    67486744                     IF ( rt2_track_lad(iz, i) > 0._wp )  THEN
    67496745                        curtrans = exp(-ext_coef * rt2_track_lad(iz, i) * (rt2_dist(l)-rt2_dist(l-1)))
    6750                         IF (ncsfl+1 > SIZE(acsf) )                             &
    6751                            CALL merge_and_grow_csf(CEILING(REAL(ncsfl, kind=wp) * grow_factor))
    6752 
    67536746                        ncsfl = ncsfl + 1
    67546747                        acsf(ncsfl)%ip = ip
     
    68106803
    68116804                     IF ( create_csf )  THEN
    6812                         IF (ncsfl+1 > SIZE(acsf) )                             &
    6813                            CALL merge_and_grow_csf(CEILING(REAL(ncsfl, kind=wp) * grow_factor))
    68146805                        ncsfl = ncsfl + 1
    68156806                        acsf(ncsfl)%ip = ip
Note: See TracChangeset for help on using the changeset viewer.