Changeset 4529


Ignore:
Timestamp:
May 12, 2020 9:14:57 AM (11 months ago)
Author:
moh.hefny
Message:

New features in coupling RTM-radiation model: consider PC-LW interactions, consider PC emiss. in the effective urban emiss., new algorithm for better performance, minor formatting changes

File:
1 edited

Legend:

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

    r4517 r4529  
    2828! -----------------
    2929! $Id$
     30! - added the following new features to the coupling of RTM-radiation model:
     31!   1) considering the vegetation interaction with LW in the coupling
     32!   2) considering PC emissivity in calculating the effective emissivity
     33!   3) new algorithm for claculating the coupling parameters so that each term
     34!      is calculated within its original line and not at the end.
     35! - minor formatting and comments changes
     36!
     37! 4517 2020-05-03 14:29:30Z raasch
    3038! added restart with MPI-IO for reading local arrays
    3139!
     
    58055813     REAL(wp)                          ::  asrc               !< area of source face
    58065814     REAL(wp)                          ::  pcrad              !< irradiance from plant canopy
    5807      REAL(wp)                          ::  pabsswl  = 0.0_wp  !< total absorbed SW radiation energy in local processor (W)
    5808      REAL(wp)                          ::  pabssw   = 0.0_wp  !< total absorbed SW radiation energy in all processors (W)
    5809      REAL(wp)                          ::  pabslwl  = 0.0_wp  !< total absorbed LW radiation energy in local processor (W)
    5810      REAL(wp)                          ::  pabslw   = 0.0_wp  !< total absorbed LW radiation energy in all processors (W)
    5811      REAL(wp)                          ::  pemitlwl = 0.0_wp  !< total emitted LW radiation energy in all processors (W)
    5812      REAL(wp)                          ::  pemitlw  = 0.0_wp  !< total emitted LW radiation energy in all processors (W)
    5813      REAL(wp)                          ::  pinswl   = 0.0_wp  !< total received SW radiation energy in local processor (W)
    5814      REAL(wp)                          ::  pinsw    = 0.0_wp  !< total received SW radiation energy in all processor (W)
    5815      REAL(wp)                          ::  pinlwl   = 0.0_wp  !< total received LW radiation energy in local processor (W)
    5816      REAL(wp)                          ::  pinlw    = 0.0_wp  !< total received LW radiation energy in all processor (W)
    5817      REAL(wp)                          ::  emiss_sum_surfl    !< sum of emissisivity of surfaces in local processor
    5818      REAL(wp)                          ::  emiss_sum_surf     !< sum of emissisivity of surfaces in all processor
    5819      REAL(wp)                          ::  area_surfl         !< total area of surfaces in local processor
    5820      REAL(wp)                          ::  area_surf          !< total area of surfaces in all processor
    5821      REAL(wp)                          ::  area_hor           !< total horizontal area of domain in all processor
     5815!-   variables for coupling the radiation modle (e.g. RRTMG) and RTM
     5816     REAL(wp)                          ::  pabsswl            !< total absorbed SW radiation energy in local processor (W)
     5817     REAL(wp)                          ::  pabssw             !< total absorbed SW radiation energy in all processors (W)
     5818     REAL(wp)                          ::  pabslwl            !< total absorbed LW radiation energy in local processor (W)
     5819     REAL(wp)                          ::  pabslw             !< total absorbed LW radiation energy in all processors (W)
     5820     REAL(wp)                          ::  pemitlwl           !< total emitted LW radiation energy in all processors (W)
     5821     REAL(wp)                          ::  pemitlw            !< total emitted LW radiation energy in all processors (W)
     5822     REAL(wp)                          ::  pinswl             !< total received SW radiation energy in local processor (W)
     5823     REAL(wp)                          ::  pinsw              !< total received SW radiation energy in all processor (W)
     5824     REAL(wp)                          ::  pinlwl             !< total received LW radiation energy in local processor (W)
     5825     REAL(wp)                          ::  pinlw              !< total received LW radiation energy in all processor (W)
     5826     REAL(wp)                          ::  area_norm          !< reference horizontal area of domain in all processor
     5827     REAL(wp)                          ::  pabs_surf_lwdifl   !< total absorbed LW radiation in surfaces from sky in local processor (W)
     5828     REAL(wp)                          ::  pabs_surf_lwdif    !< total absorbed LW radiation in surfaces from sky in all processors (W)
     5829     REAL(wp)                          ::  pabs_pc_lwdifl     !< total absorbed LW radiation in plant canopy from sky in local processor (W)
     5830     REAL(wp)                          ::  pabs_pc_lwdif      !< total absorbed LW radiation in plant canopy from sky in all processors (W)
     5831
    58225832     REAL(wp)                          ::  sun_direct_factor  !< factor for direct normal radiation from direct horizontal
    58235833#if defined( __parallel )
     
    58885898     ENDIF
    58895899
    5890 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    5891 !--     First pass: direct + diffuse irradiance + thermal
    5892 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    5893      surfinswdir   = 0._wp !nsurfl
    5894      surfins       = 0._wp !nsurfl
    5895      surfinl       = 0._wp !nsurfl
    5896      surfoutsl(:)  = 0.0_wp !start-end
    5897      surfoutll(:)  = 0.0_wp !start-end
     5900
     5901!--  First pass of radiation interaction:
     5902!--   1) direct and diffuse irradiance
     5903!--   2) thermal emissions
     5904
     5905!-   Initialize relavant surface flux arrays and radiation energy sum
     5906!-   surface flux
     5907     surfinswdir   = 0.0_wp
     5908     surfins       = 0.0_wp
     5909     surfinl       = 0.0_wp
     5910     surfoutsl(:)  = 0.0_wp
     5911     surfoutll(:)  = 0.0_wp
    58985912     IF ( nmrtbl > 0 )  THEN
    5899         mrtinsw(:) = 0._wp
    5900         mrtinlw(:) = 0._wp
     5913        mrtinsw(:) = 0.0_wp
     5914        mrtinlw(:) = 0.0_wp
    59015915     ENDIF
    5902      surfinlg(:)  = 0._wp !global
    5903 
     5916     surfinlg(:)  = 0.0_wp
     5917!-   radiation energy sum
     5918     pinlwl   = 0.0_wp
     5919     pinswl   = 0.0_wp
     5920     pemitlwl = 0.0_wp
     5921     pabsswl  = 0.0_wp
     5922     pabslwl  = 0.0_wp
     5923     pabs_surf_lwdifl = 0.0_wp
     5924     pabs_pc_lwdifl   = 0.0_wp
    59045925
    59055926!--  Set up thermal radiation from surfaces
     
    60196040        j = surfl(iy, isurf)
    60206041        i = surfl(ix, isurf)
     6042        d = surfl(id, isurf)
    60216043        surfinswdif(isurf) = rad_sw_in_diff(j,i) * skyvft(isurf)
     6044!-      update received SW energy for RTM coupling
     6045        pinswl = pinswl + surfinswdif(isurf) * facearea(d)
    60226046        IF ( plant_lw_interact )  THEN
    60236047           surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvft(isurf)
     6048!-         update received LW energy for RTM coupling
     6049           pinlwl = pinlwl + surfinlwdif(isurf) * facearea(d)
    60246050        ELSE
    60256051           surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvf(isurf)
     6052!-         update received LW energy for RTM coupling
     6053           pinlwl = pinlwl + surfinlwdif(isurf) * facearea(d)
    60266054        ENDIF
    60276055     ENDDO
     
    60536081           j = surfl(iy, isurf)
    60546082           i = surfl(ix, isurf)
     6083           d = surfl(id, isurf)
    60556084           surfinswdir(isurf) = rad_sw_in_dir(j,i) *                        &
    60566085                costheta(surfl(id, isurf)) * dsitrans(isurf, isd) * sun_direct_factor
     6086!-        update received SW energy for RTM coupling
     6087           pinswl = pinswl + surfinswdir(isurf) * facearea(d)
    60576088        ENDDO
    60586089!
     
    60926123!--             Diffuse radiation from sky
    60936124                pcbinswdif(ipcgb) = csf(1,icsf) * rad_sw_in_diff(j,i)
     6125!--             add to the sum of SW radiation energy
     6126                pinswl = pinswl + pcbinswdif(ipcgb)
    60946127!
    60956128!--             Absorbed diffuse LW radiation from sky minus emitted to sky
     
    60986131                                       * (rad_lw_in_diff(j, i)                   &
    60996132                                          - sigma_sb * (pt(k,j,i)*exner(k))**4)
     6133                   pinlwl = pinlwl + csf(1,icsf) * rad_lw_in_diff(j,i)
     6134                   pabslwl = pabslwl + csf(1,icsf) * rad_lw_in_diff(j,i)
     6135                   pemitlwl = pemitlwl +                                         &
     6136                              csf(1,icsf) * sigma_sb * (pt(k,j,i)*exner(k))**4
     6137                   pabs_pc_lwdifl = pabs_pc_lwdifl +                             &
     6138                                    csf(1,icsf) * rad_lw_in_diff(j,i)
    61006139                ENDIF
    61016140!
     
    61086147                   pcbinswdir(ipcgb) = rad_sw_in_dir(j, i) * pc_box_area &
    61096148                                       * pc_abs_frac * dsitransc(ipcgb, isd)
     6149!--               add to the sum of SW radiation energy
     6150                  pinswl = pinswl + pcbinswdir(ipcgb)
    61106151                ENDIF
    61116152             ELSE
     
    61226163                                       - pcrad)                         & ! Remove emitted heatflux
    61236164                                    * asrc
     6165                   pabslwl = pabslwl + csf(1,icsf) * surfoutl(isurfsrc) * asrc
     6166                   pemitlwl = pemitlwl + pcrad * asrc
    61246167                ENDIF
    61256168             ENDIF
     
    61786221     ENDIF
    61796222
    6180 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    6181 !--     Next passes - reflections
    6182 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     6223!-- Next passes of radiation interactions:
     6224!--  radiation reflections
     6225
    61836226     DO refstep = 1, nrefsteps
    61846227
     
    62736316
    62746317     ENDDO ! refstep
    6275 
    6276 !--  push heat flux absorbed by plant canopy to respective 3D arrays
     6318!
     6319!--  push heat flux absorbed by plant canopy to respective 3D arrays and
     6320!--  add absorbed SW radiation energy for RTM coupling variables
    62776321     IF ( npcbl > 0 )  THEN
    62786322         pcm_heating_rate(:,:,:) = 0.0_wp
     
    62866330             pcm_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) &
    62876331                 * pchf_prep(k) * pt(k, j, i) !-- = dT/dt
     6332!--          add the absorbed SW radiation energy by plant canopy
     6333             pabsswl = pabsswl + pcbinsw(ipcgb)
    62886334         ENDDO
    62896335
     
    62976343                 k = pcbl(iz, ipcgb)
    62986344                 kk = k - topo_top_ind(j,i,0)  !- lad arrays are defined flat
    6299                  CALL pcm_calc_transpiration_rate( i, j, k, kk, pcbinsw(ipcgb), pcbinlw(ipcgb), &
    6300                                                    pcm_transpiration_rate(kk,j,i), pcm_latent_rate(kk,j,i) )
     6345                 CALL pcm_calc_transpiration_rate( i, j, k, kk, pcbinsw(ipcgb),     &
     6346                                                   pcbinlw(ipcgb),                  &
     6347                                                   pcm_transpiration_rate(kk,j,i),  &
     6348                                                   pcm_latent_rate(kk,j,i) )
    63016349              ENDDO
    63026350         ENDIF
     
    63126360     ENDIF
    63136361!
    6314 !--     Transfer radiation arrays required for energy balance to the respective data types
     6362!--  Transfer radiation arrays required for energy balance to the respective data types
     6363!    and claculate relevant radiation model-RTM coupling terms
     6364
    63156365     DO  i = 1, nsurfl
    63166366        m  = surfl(im,i)
     6367        d  = surfl(id, i)
    63176368!
    63186369!--     (1) Urban surfaces
     
    64936544           surf_lsm_v(3)%rad_lw_res(m) = surfinl(i)
    64946545        ENDIF
    6495 
     6546!
     6547!--     RTM coupling terms
     6548!--     sum of absorbed SW & LW radiation energy
     6549        pabsswl = pabsswl +                                                 &
     6550                  (1.0_wp - albedo_surf(i)) * surfinsw(i) * facearea(d)
     6551        pabslwl = pabslwl + emiss_surf(i) * surfinlw(i) * facearea(d)
     6552!--     sum of emitted LW radiation energy
     6553        pemitlwl = pemitlwl + surfemitlwl(i) * facearea(d)
     6554!--     emiss1
     6555        pabs_surf_lwdifl = pabs_surf_lwdifl +                               &
     6556                           emiss_surf(i) * facearea(d) * surfinlwdif(i)
    64966557     ENDDO
    64976558
     
    65276588     ENDDO
    65286589!
    6529 !--  Calculate the average temperature, albedo, and emissivity for urban/land
    6530 !--  domain when using average_radiation in the respective radiation model
    6531 
    6532 !--  calculate horizontal area
    6533 ! !!! ATTENTION!!! uniform grid is assumed here
    6534      area_hor = (nx+1) * (ny+1) * dx * dy
    6535 !
    6536 !--  absorbed/received SW & LW and emitted LW energy of all physical
    6537 !--  surfaces (land and urban) in local processor
    6538      pinswl = 0._wp
    6539      pinlwl = 0._wp
    6540      pabsswl = 0._wp
    6541      pabslwl = 0._wp
    6542      pemitlwl = 0._wp
    6543      emiss_sum_surfl = 0._wp
    6544      area_surfl = 0._wp
    6545      DO  i = 1, nsurfl
    6546         d = surfl(id, i)
    6547 !--  received SW & LW
    6548         pinswl = pinswl + (surfinswdir(i) + surfinswdif(i)) * facearea(d)
    6549         pinlwl = pinlwl + surfinlwdif(i) * facearea(d)
    6550 !--   absorbed SW & LW
    6551         pabsswl = pabsswl + (1._wp - albedo_surf(i)) *                   &
    6552                                                 surfinsw(i) * facearea(d)
    6553         pabslwl = pabslwl + emiss_surf(i) * surfinlw(i) * facearea(d)
    6554 !--   emitted LW
    6555         pemitlwl = pemitlwl + surfemitlwl(i) * facearea(d)
    6556 !--   emissivity and area sum
    6557         emiss_sum_surfl = emiss_sum_surfl + emiss_surf(i) * facearea(d)
    6558         area_surfl = area_surfl + facearea(d)
    6559      END DO
    6560 !
    6561 !--  add the absorbed SW energy by plant canopy
    6562      IF ( npcbl > 0 )  THEN
    6563         pabsswl = pabsswl + SUM(pcbinsw)
    6564         pabslwl = pabslwl + SUM(pcbinlw)
    6565         pinswl  = pinswl + SUM(pcbinswdir) + SUM(pcbinswdif)
    6566      ENDIF
    6567 !
    65686590!--  gather all rad flux energy in all processors. In order to reduce
    65696591!--  the number of MPI calls (to reduce latencies), combine the required
     
    65766598     combine_allreduce_l(4) = pabslwl
    65776599     combine_allreduce_l(5) = pemitlwl
    6578      combine_allreduce_l(6) = emiss_sum_surfl
    6579      combine_allreduce_l(7) = area_surfl
     6600     combine_allreduce_l(6) = pabs_surf_lwdifl
     6601     combine_allreduce_l(7) = pabs_pc_lwdifl
    65806602
    65816603     CALL MPI_ALLREDUCE( combine_allreduce_l,                                  &
     
    65876609                         ierr )
    65886610
    6589      pinsw          = combine_allreduce(1)
    6590      pinlw          = combine_allreduce(2)
    6591      pabssw         = combine_allreduce(3)
    6592      pabslw         = combine_allreduce(4)
    6593      pemitlw        = combine_allreduce(5)
    6594      emiss_sum_surf = combine_allreduce(6)
    6595      area_surf      = combine_allreduce(7)
     6611     pinsw           = combine_allreduce(1)
     6612     pinlw           = combine_allreduce(2)
     6613     pabssw          = combine_allreduce(3)
     6614     pabslw          = combine_allreduce(4)
     6615     pemitlw         = combine_allreduce(5)
     6616     pabs_surf_lwdif = combine_allreduce(6)
     6617     pabs_pc_lwdif   = combine_allreduce(7)
    65966618#else
    6597      pinsw          = pinswl
    6598      pinlw          = pinlwl
    6599      pabssw         = pabsswl
    6600      pabslw         = pabslwl
    6601      pemitlw        = pemitlwl
    6602      emiss_sum_surf = emiss_sum_surfl
    6603      area_surf      = area_surfl
     6619     pinsw           = pinswl
     6620     pinlw           = pinlwl
     6621     pabssw          = pabsswl
     6622     pabslw          = pabslwl
     6623     pemitlw         = pemitlwl
     6624     pabs_surf_lwdif = pabs_surf_lwdifl
     6625     pabs_pc_lwdif   = pabs_pc_lwdif
    66046626#endif
    6605 
    6606 !--  (1) albedo
     6627!
     6628!--  Calculate the effective radiation surface parameters based on
     6629!--  the parameterizations in Krc et al. 2020
     6630
     6631!-   (1) albedo Eq. * in Krc et al. 2020
    66076632     IF ( pinsw /= 0.0_wp )  albedo_urb = ( pinsw - pabssw ) / pinsw
    6608 !--  (2) average emmsivity
    6609      IF ( area_surf /= 0.0_wp )  emissivity_urb = emiss_sum_surf / area_surf
    6610 !
    6611 !--  Temporally comment out calculation of effective radiative temperature.
    6612 !--  See below for more explanation.
    6613 !--  (3) temperature
    6614 !--   first we calculate an effective horizontal area to account for
    6615 !--   the effect of vertical surfaces (which contributes to LW emission)
    6616 !--   We simply use the ratio of the total LW to the incoming LW flux
    6617       area_hor = pinlw / rad_lw_in_diff(nyn,nxl)
    6618       t_rad_urb = ( ( pemitlw - pabslw + emissivity_urb * pinlw ) / &
    6619            (emissivity_urb * sigma_sb * area_hor) )**0.25_wp
     6633
     6634!-   (2) emmsivity Eq. * in Krc et al. 2020
     6635!-    emissivity_urb weighted average of surface and PC emissivity
     6636!-    = absorbed LW in [surfaces + plant canopy] / pinlw
     6637     IF ( pinsw /= 0.0_wp )                                                         &
     6638          emissivity_urb = (pabs_surf_lwdif + pabs_pc_lwdif) / pinlw
     6639
     6640!-    (3) temperature
     6641!-     effective horizontal area to account for the effect of vertical
     6642!-     surfaces, Eq. * in Krc et al. 2020
     6643      area_norm = pinlw / rad_lw_in_diff(nyn,nxl)
     6644!-     temperature, Eq. * in Krc et al. 2020
     6645      t_rad_urb = ( (pemitlw - pabslw + emissivity_urb * pinlw) /                    &
     6646           (emissivity_urb * sigma_sb * area_norm) )**0.25_wp
    66206647
    66216648     IF ( debug_output_timestep )  CALL debug_message( 'radiation_interaction', 'end' )
Note: See TracChangeset for help on using the changeset viewer.