Changeset 3065 for palm/trunk/SOURCE/radiation_model_mod.f90

Timestamp:

Jun 12, 2018 7:03:02 AM (6 years ago)

Author:

Giersch

Message:

New vertical stretching procedure has been introduced

File:

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

palm/trunk/SOURCE/radiation_model_mod.f90

@@ -28 +28 @@
 ! -----------------
 ! $Id$
+! dz was replaced by dz(1), error message concerning vertical stretching was 
+! added  
+! 
+! 3049 2018-05-29 13:52:36Z Giersch
 ! Error messages revised
 ! 
@@ -330 +334 @@
 ! ------------
 !> Radiation models and interfaces
+!> @todo Replace dz(1) appropriatly to account for grid stretching
 !> @todo move variable definitions used in radiation_init only to the subroutine
 !>       as they are no longer required after initialization.
@@ -4371 +4376 @@
      IF ( plant_canopy )  THEN
          pchf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp &
-                     / (cp * hyp(nzub:nzut) * dx*dy*dz) !< equals to 1 / (rho * c_p * Vbox * T)
+                     / (cp * hyp(nzub:nzut) * dx*dy*dz(1)) !< equals to 1 / (rho * c_p * Vbox * T)
          pctf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp &
                      / (l_v * hyp(nzub:nzut) * dx*dy*dz)
@@ -4396 +4401 @@
 !--         each dimension, so that this new direction vector will allow us
 !--         to traverse the ray path within grid coordinates directly
-         sunorig_grid = (/ sunorig(1)/dz, sunorig(2)/dy, sunorig(3)/dx /)
+         sunorig_grid = (/ sunorig(1)/dz(1), sunorig(2)/dy, sunorig(3)/dx /)
 !--         sunorig_grid = sunorig_grid / norm2(sunorig_grid)
          sunorig_grid = sunorig_grid / SQRT(SUM(sunorig_grid**2))
@@ -4403 +4408 @@
 !--            precompute effective box depth with prototype Leaf Area Density
             pc_box_dimshift = MAXLOC(ABS(sunorig), 1) - 1
-            CALL box_absorb(CSHIFT((/dz,dy,dx/), pc_box_dimshift),      &
+            CALL box_absorb(CSHIFT((/dz(1),dy,dx/), pc_box_dimshift),      &
                                 60, prototype_lad,                          &
                                 CSHIFT(ABS(sunorig), pc_box_dimshift),      &
@@ -4801 +4806 @@
         IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx
         IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy
-        IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz
+        IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz(1)
      ENDDO
 !
@@ -5014 +5019 @@
     SUBROUTINE radiation_interaction_init
 
+       USE control_parameters,                                                 &
+           ONLY:  dz_stretch_level_start
+           
        USE netcdf_data_input_mod,                                              &
            ONLY:  leaf_area_density_f
@@ -5095 +5103 @@
        nzpt = nzptl
 #endif
+!
+!--    Stretching (non-uniform grid spacing) is not considered in the radiation 
+!--    model. Therefore, vertical stretching has to be applied above the area 
+!--    where the parts of the radiation model which assume constant grid spacing
+!--    are active. ABS (...) is required because the default value of 
+!--    dz_stretch_level_start is -9999999.9_wp (negative). 
+       IF ( ABS( dz_stretch_level_start(1) ) <= zw(nzut) ) THEN
+          WRITE( message_string, * ) 'The lowest level where vertical ',       &
+                                     'stretching is applied have to be ',      &
+                                     'greater than ', zw(nzut)
+          CALL message( 'init_grid', 'PA0496', 1, 2, 0, 6, 0 )
+       ENDIF 
 !
 !--    global number of urban and plant layers
@@ -5390 +5410 @@
             IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx
             IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy
-            IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz
+            IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz(1)
         ENDDO
 
@@ -5512 +5532 @@
 
 !--                 unit vector source -> target
-                    uv = (/ (ta(1)-sa(1))*dz, (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
+                    uv = (/ (ta(1)-sa(1))*dz(1), (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
                     sqdist = SUM(uv(:)**2)
                     uv = uv / SQRT(sqdist)
@@ -5677 +5697 @@
 
 !--               unit vector source -> target
-                  uv = (/ (ta(1)-sa(1))*dz, (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
+                  uv = (/ (ta(1)-sa(1))*dz(1), (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
                   sqdist = SUM(uv(:)**2)
                   uv = uv / SQRT(sqdist)
@@ -6127 +6147 @@
         ENDIF
         uvect(:) = delta(:) / distance
-        realdist = SQRT(SUM( (uvect(:)*(/dz,dy,dx/))**2 ))
+        realdist = SQRT(SUM( (uvect(:)*(/dz(1),dy,dx/))**2 ))
 
         lastdist = 0._wp
@@ -6509 +6529 @@
                      zb1 = CEILING(zexit * zsgn - .5_wp) - 1  ! because it must be smaller than exit
                      nz = MAX(zb1 - zb0 + 3, 2)
-                     rt2_dist(nz) = SQRT(((zexit-zorig)*dz)**2 + dxxyy)
+                     rt2_dist(nz) = SQRT(((zexit-zorig)*dz(1))**2 + dxxyy)
                      qdist = rt2_dist(nz) / (zexit-zorig)
                      rt2_dist(2:nz-1) = (/( ((REAL(l, wp) + .5_wp) * zsgn - zorig) * qdist , l = zb0, zb1 )/)
@@ -6567 +6587 @@
                   zb1 = CEILING(zexit * zsgn - .5_wp) - 1  ! because it must be smaller than exit
                   nz = MAX(zb1 - zb0 + 3, 2)
-                  rt2_dist(nz) = SQRT(((zexit-zorig)*dz)**2 + dxxyy)
+                  rt2_dist(nz) = SQRT(((zexit-zorig)*dz(1))**2 + dxxyy)
                   qdist = rt2_dist(nz) / (zexit-zorig)
                   rt2_dist(2:nz-1) = (/( ((REAL(l, wp) + .5_wp) * zsgn - zorig) * qdist , l = zb0, zb1 )/)