Changeset 4584

Timestamp:

Jun 29, 2020 1:16:14 PM (4 years ago)

Author:

pavelkrc

Message:

Consider only boxes with LAD>0 as plant canopy

File:

• palm/trunk/SOURCE/radiation_model_mod.f90 (modified) (4 diffs)

palm/trunk/SOURCE/radiation_model_mod.f90

@@ -28 +28 @@
 ! -----------------
 ! $Id$
+! Consider only boxes with LAD>0 as plant canopy (credit: S. Schubert)
+! 
+! 4576 2020-06-24 17:58:55Z pavelkrc
 ! Allow the use of rotation_angle in RTM
 ! 
@@ -6918 +6921 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) :: i, j, k, l, m, d
-       INTEGER(iwp) :: k_topo     !< vertical index indicating topography top for given (j,i)
-       INTEGER(iwp) :: nzptl, nzubl, nzutl, isurf, ipcgb, imrt
-       REAL(wp)     :: mrl
+       INTEGER(iwp)        ::  i, j, k, l, m, d
+       INTEGER(iwp)        ::  k_topo               !< vertical index indicating
+                                                    !< topography top for given (j,i)
+       INTEGER(iwp)        ::  isurf, ipcgb, imrt
+       INTEGER(iwp)        ::  nzptl, nzubl, nzutl
+       REAL(wp)            ::  mrl
+       REAL(wp), PARAMETER ::  eps_lad = 1E-10_wp   !< minimum considered nonzero
 #if defined( __parallel )
        INTEGER(iwp), DIMENSION(:), POINTER, SAVE ::  gridsurf_rma   !< fortran pointer, but lower bounds are 1
@@ -6960 +6966 @@
                    k_topo = topo_top_ind(j,i,0)
 
-                   DO k = nzt+1, 0, -1
-                       IF ( lad_s(k,j,i) /= 0.0_wp )  THEN
+                   DO k = nzt+1, 1, -1
+                       IF ( lad_s(k,j,i) > eps_lad )  THEN
 !--                        we are at the top of the pcs
                            pct(j,i) = k + k_topo
                            pch(j,i) = k
-                           npcbl = npcbl + pch(j,i)
+                           npcbl = npcbl + 1 + COUNT(lad_s(1:k-1,j,i) > eps_lad)
                            EXIT
                        ENDIF
@@ -7084 +7090 @@
 
                    DO k = k_topo + 1, pct(j,i)
-                       ipcgb = ipcgb + 1
-                       gridpcbl(k,j,i) = ipcgb
-                       pcbl(:,ipcgb) = (/ k, j, i /)
+                       IF ( lad_s(k-k_topo,j,i) > eps_lad )  THEN
+                          ipcgb = ipcgb + 1
+                          gridpcbl(k,j,i) = ipcgb
+                          pcbl(:,ipcgb) = (/ k, j, i /)
+                       ENDIF
                    ENDDO
                ENDDO