Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4528)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4529)
@@ -28,4 +28,12 @@
 ! -----------------
 ! $Id$
+! - added the following new features to the coupling of RTM-radiation model:
+!   1) considering the vegetation interaction with LW in the coupling
+!   2) considering PC emissivity in calculating the effective emissivity
+!   3) new algorithm for claculating the coupling parameters so that each term
+!      is calculated within its original line and not at the end.
+! - minor formatting and comments changes
+! 
+! 4517 2020-05-03 14:29:30Z raasch
 ! added restart with MPI-IO for reading local arrays
 ! 
@@ -5805,19 +5813,21 @@
      REAL(wp)                          ::  asrc               !< area of source face
      REAL(wp)                          ::  pcrad              !< irradiance from plant canopy
-     REAL(wp)                          ::  pabsswl  = 0.0_wp  !< total absorbed SW radiation energy in local processor (W)
-     REAL(wp)                          ::  pabssw   = 0.0_wp  !< total absorbed SW radiation energy in all processors (W)
-     REAL(wp)                          ::  pabslwl  = 0.0_wp  !< total absorbed LW radiation energy in local processor (W)
-     REAL(wp)                          ::  pabslw   = 0.0_wp  !< total absorbed LW radiation energy in all processors (W)
-     REAL(wp)                          ::  pemitlwl = 0.0_wp  !< total emitted LW radiation energy in all processors (W)
-     REAL(wp)                          ::  pemitlw  = 0.0_wp  !< total emitted LW radiation energy in all processors (W)
-     REAL(wp)                          ::  pinswl   = 0.0_wp  !< total received SW radiation energy in local processor (W)
-     REAL(wp)                          ::  pinsw    = 0.0_wp  !< total received SW radiation energy in all processor (W)
-     REAL(wp)                          ::  pinlwl   = 0.0_wp  !< total received LW radiation energy in local processor (W)
-     REAL(wp)                          ::  pinlw    = 0.0_wp  !< total received LW radiation energy in all processor (W)
-     REAL(wp)                          ::  emiss_sum_surfl    !< sum of emissisivity of surfaces in local processor
-     REAL(wp)                          ::  emiss_sum_surf     !< sum of emissisivity of surfaces in all processor
-     REAL(wp)                          ::  area_surfl         !< total area of surfaces in local processor
-     REAL(wp)                          ::  area_surf          !< total area of surfaces in all processor
-     REAL(wp)                          ::  area_hor           !< total horizontal area of domain in all processor
+!-   variables for coupling the radiation modle (e.g. RRTMG) and RTM
+     REAL(wp)                          ::  pabsswl            !< total absorbed SW radiation energy in local processor (W)
+     REAL(wp)                          ::  pabssw             !< total absorbed SW radiation energy in all processors (W)
+     REAL(wp)                          ::  pabslwl            !< total absorbed LW radiation energy in local processor (W)
+     REAL(wp)                          ::  pabslw             !< total absorbed LW radiation energy in all processors (W)
+     REAL(wp)                          ::  pemitlwl           !< total emitted LW radiation energy in all processors (W)
+     REAL(wp)                          ::  pemitlw            !< total emitted LW radiation energy in all processors (W)
+     REAL(wp)                          ::  pinswl             !< total received SW radiation energy in local processor (W)
+     REAL(wp)                          ::  pinsw              !< total received SW radiation energy in all processor (W)
+     REAL(wp)                          ::  pinlwl             !< total received LW radiation energy in local processor (W)
+     REAL(wp)                          ::  pinlw              !< total received LW radiation energy in all processor (W)
+     REAL(wp)                          ::  area_norm          !< reference horizontal area of domain in all processor
+     REAL(wp)                          ::  pabs_surf_lwdifl   !< total absorbed LW radiation in surfaces from sky in local processor (W)
+     REAL(wp)                          ::  pabs_surf_lwdif    !< total absorbed LW radiation in surfaces from sky in all processors (W)
+     REAL(wp)                          ::  pabs_pc_lwdifl     !< total absorbed LW radiation in plant canopy from sky in local processor (W)
+     REAL(wp)                          ::  pabs_pc_lwdif      !< total absorbed LW radiation in plant canopy from sky in all processors (W)
+
      REAL(wp)                          ::  sun_direct_factor  !< factor for direct normal radiation from direct horizontal
 #if defined( __parallel )
@@ -5888,18 +5898,29 @@
      ENDIF
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!--     First pass: direct + diffuse irradiance + thermal
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-     surfinswdir   = 0._wp !nsurfl
-     surfins       = 0._wp !nsurfl
-     surfinl       = 0._wp !nsurfl
-     surfoutsl(:)  = 0.0_wp !start-end
-     surfoutll(:)  = 0.0_wp !start-end
+
+!--  First pass of radiation interaction:
+!--   1) direct and diffuse irradiance 
+!--   2) thermal emissions
+
+!-   Initialize relavant surface flux arrays and radiation energy sum
+!-   surface flux
+     surfinswdir   = 0.0_wp
+     surfins       = 0.0_wp
+     surfinl       = 0.0_wp
+     surfoutsl(:)  = 0.0_wp
+     surfoutll(:)  = 0.0_wp
      IF ( nmrtbl > 0 )  THEN
-        mrtinsw(:) = 0._wp
-        mrtinlw(:) = 0._wp
+        mrtinsw(:) = 0.0_wp
+        mrtinlw(:) = 0.0_wp
      ENDIF
-     surfinlg(:)  = 0._wp !global
-
+     surfinlg(:)  = 0.0_wp
+!-   radiation energy sum
+     pinlwl   = 0.0_wp
+     pinswl   = 0.0_wp
+     pemitlwl = 0.0_wp
+     pabsswl  = 0.0_wp
+     pabslwl  = 0.0_wp
+     pabs_surf_lwdifl = 0.0_wp
+     pabs_pc_lwdifl   = 0.0_wp
 
 !--  Set up thermal radiation from surfaces
@@ -6019,9 +6040,16 @@
         j = surfl(iy, isurf)
         i = surfl(ix, isurf)
+        d = surfl(id, isurf)
         surfinswdif(isurf) = rad_sw_in_diff(j,i) * skyvft(isurf)
+!-      update received SW energy for RTM coupling
+        pinswl = pinswl + surfinswdif(isurf) * facearea(d)
         IF ( plant_lw_interact )  THEN
            surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvft(isurf)
+!-         update received LW energy for RTM coupling
+           pinlwl = pinlwl + surfinlwdif(isurf) * facearea(d)
         ELSE
            surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvf(isurf)
+!-         update received LW energy for RTM coupling
+           pinlwl = pinlwl + surfinlwdif(isurf) * facearea(d)
         ENDIF
      ENDDO
@@ -6053,6 +6081,9 @@
            j = surfl(iy, isurf)
            i = surfl(ix, isurf)
+           d = surfl(id, isurf)
            surfinswdir(isurf) = rad_sw_in_dir(j,i) *                        &
                 costheta(surfl(id, isurf)) * dsitrans(isurf, isd) * sun_direct_factor
+!-        update received SW energy for RTM coupling
+           pinswl = pinswl + surfinswdir(isurf) * facearea(d)
         ENDDO
 !
@@ -6092,4 +6123,6 @@
 !--             Diffuse radiation from sky
                 pcbinswdif(ipcgb) = csf(1,icsf) * rad_sw_in_diff(j,i)
+!--             add to the sum of SW radiation energy
+                pinswl = pinswl + pcbinswdif(ipcgb)
 !
 !--             Absorbed diffuse LW radiation from sky minus emitted to sky
@@ -6098,4 +6131,10 @@
                                        * (rad_lw_in_diff(j, i)                   &
                                           - sigma_sb * (pt(k,j,i)*exner(k))**4)
+                   pinlwl = pinlwl + csf(1,icsf) * rad_lw_in_diff(j,i)
+                   pabslwl = pabslwl + csf(1,icsf) * rad_lw_in_diff(j,i)
+                   pemitlwl = pemitlwl +                                         &
+                              csf(1,icsf) * sigma_sb * (pt(k,j,i)*exner(k))**4
+                   pabs_pc_lwdifl = pabs_pc_lwdifl +                             &
+                                    csf(1,icsf) * rad_lw_in_diff(j,i)
                 ENDIF
 !
@@ -6108,4 +6147,6 @@
                    pcbinswdir(ipcgb) = rad_sw_in_dir(j, i) * pc_box_area &
                                        * pc_abs_frac * dsitransc(ipcgb, isd)
+!--               add to the sum of SW radiation energy
+                  pinswl = pinswl + pcbinswdir(ipcgb)
                 ENDIF
              ELSE
@@ -6122,4 +6163,6 @@
                                        - pcrad)                         & ! Remove emitted heatflux
                                     * asrc
+                   pabslwl = pabslwl + csf(1,icsf) * surfoutl(isurfsrc) * asrc
+                   pemitlwl = pemitlwl + pcrad * asrc
                 ENDIF
              ENDIF
@@ -6178,7 +6221,7 @@
      ENDIF
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!--     Next passes - reflections
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!-- Next passes of radiation interactions:
+!--  radiation reflections
+
      DO refstep = 1, nrefsteps
 
@@ -6273,6 +6316,7 @@
 
      ENDDO ! refstep
-
-!--  push heat flux absorbed by plant canopy to respective 3D arrays
+!
+!--  push heat flux absorbed by plant canopy to respective 3D arrays and
+!--  add absorbed SW radiation energy for RTM coupling variables
      IF ( npcbl > 0 )  THEN
          pcm_heating_rate(:,:,:) = 0.0_wp
@@ -6286,4 +6330,6 @@
              pcm_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) &
                  * pchf_prep(k) * pt(k, j, i) !-- = dT/dt
+!--          add the absorbed SW radiation energy by plant canopy
+             pabsswl = pabsswl + pcbinsw(ipcgb)
          ENDDO
 
@@ -6297,6 +6343,8 @@
                  k = pcbl(iz, ipcgb)
                  kk = k - topo_top_ind(j,i,0)  !- lad arrays are defined flat
-                 CALL pcm_calc_transpiration_rate( i, j, k, kk, pcbinsw(ipcgb), pcbinlw(ipcgb), &
-                                                   pcm_transpiration_rate(kk,j,i), pcm_latent_rate(kk,j,i) )
+                 CALL pcm_calc_transpiration_rate( i, j, k, kk, pcbinsw(ipcgb),     &
+                                                   pcbinlw(ipcgb),                  &
+                                                   pcm_transpiration_rate(kk,j,i),  &
+                                                   pcm_latent_rate(kk,j,i) )
               ENDDO
          ENDIF
@@ -6312,7 +6360,10 @@
      ENDIF
 !
-!--     Transfer radiation arrays required for energy balance to the respective data types
+!--  Transfer radiation arrays required for energy balance to the respective data types
+!    and claculate relevant radiation model-RTM coupling terms
+
      DO  i = 1, nsurfl
         m  = surfl(im,i)
+        d  = surfl(id, i)
 !
 !--     (1) Urban surfaces
@@ -6493,5 +6544,15 @@
            surf_lsm_v(3)%rad_lw_res(m) = surfinl(i)
         ENDIF
-
+!
+!--     RTM coupling terms
+!--     sum of absorbed SW & LW radiation energy
+        pabsswl = pabsswl +                                                 &
+                  (1.0_wp - albedo_surf(i)) * surfinsw(i) * facearea(d)
+        pabslwl = pabslwl + emiss_surf(i) * surfinlw(i) * facearea(d)
+!--     sum of emitted LW radiation energy
+        pemitlwl = pemitlwl + surfemitlwl(i) * facearea(d)
+!--     emiss1
+        pabs_surf_lwdifl = pabs_surf_lwdifl +                               &
+                           emiss_surf(i) * facearea(d) * surfinlwdif(i) 
      ENDDO
 
@@ -6527,43 +6588,4 @@
      ENDDO
 !
-!--  Calculate the average temperature, albedo, and emissivity for urban/land
-!--  domain when using average_radiation in the respective radiation model
-
-!--  calculate horizontal area
-! !!! ATTENTION!!! uniform grid is assumed here
-     area_hor = (nx+1) * (ny+1) * dx * dy
-!
-!--  absorbed/received SW & LW and emitted LW energy of all physical
-!--  surfaces (land and urban) in local processor
-     pinswl = 0._wp
-     pinlwl = 0._wp
-     pabsswl = 0._wp
-     pabslwl = 0._wp
-     pemitlwl = 0._wp
-     emiss_sum_surfl = 0._wp
-     area_surfl = 0._wp
-     DO  i = 1, nsurfl
-        d = surfl(id, i)
-!--  received SW & LW
-        pinswl = pinswl + (surfinswdir(i) + surfinswdif(i)) * facearea(d)
-        pinlwl = pinlwl + surfinlwdif(i) * facearea(d)
-!--   absorbed SW & LW
-        pabsswl = pabsswl + (1._wp - albedo_surf(i)) *                   &
-                                                surfinsw(i) * facearea(d)
-        pabslwl = pabslwl + emiss_surf(i) * surfinlw(i) * facearea(d)
-!--   emitted LW
-        pemitlwl = pemitlwl + surfemitlwl(i) * facearea(d)
-!--   emissivity and area sum
-        emiss_sum_surfl = emiss_sum_surfl + emiss_surf(i) * facearea(d)
-        area_surfl = area_surfl + facearea(d)
-     END DO
-!
-!--  add the absorbed SW energy by plant canopy
-     IF ( npcbl > 0 )  THEN
-        pabsswl = pabsswl + SUM(pcbinsw)
-        pabslwl = pabslwl + SUM(pcbinlw)
-        pinswl  = pinswl + SUM(pcbinswdir) + SUM(pcbinswdif)
-     ENDIF
-!
 !--  gather all rad flux energy in all processors. In order to reduce
 !--  the number of MPI calls (to reduce latencies), combine the required
@@ -6576,6 +6598,6 @@
      combine_allreduce_l(4) = pabslwl
      combine_allreduce_l(5) = pemitlwl
-     combine_allreduce_l(6) = emiss_sum_surfl
-     combine_allreduce_l(7) = area_surfl
+     combine_allreduce_l(6) = pabs_surf_lwdifl
+     combine_allreduce_l(7) = pabs_pc_lwdifl
 
      CALL MPI_ALLREDUCE( combine_allreduce_l,                                  &
@@ -6587,35 +6609,40 @@
                          ierr )
 
-     pinsw          = combine_allreduce(1)
-     pinlw          = combine_allreduce(2)
-     pabssw         = combine_allreduce(3)
-     pabslw         = combine_allreduce(4)
-     pemitlw        = combine_allreduce(5)
-     emiss_sum_surf = combine_allreduce(6)
-     area_surf      = combine_allreduce(7)
+     pinsw           = combine_allreduce(1)
+     pinlw           = combine_allreduce(2)
+     pabssw          = combine_allreduce(3)
+     pabslw          = combine_allreduce(4)
+     pemitlw         = combine_allreduce(5)
+     pabs_surf_lwdif = combine_allreduce(6)
+     pabs_pc_lwdif   = combine_allreduce(7)
 #else
-     pinsw          = pinswl
-     pinlw          = pinlwl
-     pabssw         = pabsswl
-     pabslw         = pabslwl
-     pemitlw        = pemitlwl
-     emiss_sum_surf = emiss_sum_surfl
-     area_surf      = area_surfl
+     pinsw           = pinswl
+     pinlw           = pinlwl
+     pabssw          = pabsswl
+     pabslw          = pabslwl
+     pemitlw         = pemitlwl
+     pabs_surf_lwdif = pabs_surf_lwdifl
+     pabs_pc_lwdif   = pabs_pc_lwdif
 #endif
-
-!--  (1) albedo
+!
+!--  Calculate the effective radiation surface parameters based on
+!--  the parameterizations in Krc et al. 2020
+
+!-   (1) albedo Eq. * in Krc et al. 2020
      IF ( pinsw /= 0.0_wp )  albedo_urb = ( pinsw - pabssw ) / pinsw
-!--  (2) average emmsivity
-     IF ( area_surf /= 0.0_wp )  emissivity_urb = emiss_sum_surf / area_surf
-!
-!--  Temporally comment out calculation of effective radiative temperature.
-!--  See below for more explanation.
-!--  (3) temperature
-!--   first we calculate an effective horizontal area to account for
-!--   the effect of vertical surfaces (which contributes to LW emission)
-!--   We simply use the ratio of the total LW to the incoming LW flux
-      area_hor = pinlw / rad_lw_in_diff(nyn,nxl)
-      t_rad_urb = ( ( pemitlw - pabslw + emissivity_urb * pinlw ) / &
-           (emissivity_urb * sigma_sb * area_hor) )**0.25_wp
+
+!-   (2) emmsivity Eq. * in Krc et al. 2020
+!-    emissivity_urb weighted average of surface and PC emissivity
+!-    = absorbed LW in [surfaces + plant canopy] / pinlw
+     IF ( pinsw /= 0.0_wp )                                                         &
+          emissivity_urb = (pabs_surf_lwdif + pabs_pc_lwdif) / pinlw
+
+!-    (3) temperature
+!-     effective horizontal area to account for the effect of vertical 
+!-     surfaces, Eq. * in Krc et al. 2020
+      area_norm = pinlw / rad_lw_in_diff(nyn,nxl)
+!-     temperature, Eq. * in Krc et al. 2020
+      t_rad_urb = ( (pemitlw - pabslw + emissivity_urb * pinlw) /                    &
+           (emissivity_urb * sigma_sb * area_norm) )**0.25_wp
 
      IF ( debug_output_timestep )  CALL debug_message( 'radiation_interaction', 'end' )
