Changeset 4575 for palm

Timestamp:

Jun 24, 2020 4:58:40 PM (4 years ago)

Author:

pavelkrc

Message:

Allow the use of rotation_angle in RTM

File:

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

palm/trunk/SOURCE/radiation_model_mod.f90

@@ -28 +28 @@
 ! -----------------
 ! $Id$
+! Allow the use of rotation_angle in RTM
+! 
+! 4574 2020-06-24 16:33:32Z pavelkrc
 ! - Restructure code in radiation_check_data_output
 ! - Move calculation of MPI global array offsets to a subroutine
@@ -1485 +1488 @@
 
        USE control_parameters,                                                 &
-           ONLY: land_surface, message_string, rotation_angle, urban_surface
+           ONLY: land_surface, message_string, urban_surface
 
        USE netcdf_data_input_mod,                                              &
@@ -1585 +1588 @@
           message_string = 'dt_radiation must be > 0.0'
           CALL message( 'check_parameters', 'PA0591', 1, 2, 0, 6, 0 )
-       ENDIF
-!
-!--    Check rotation angle
-       !> @todo Remove this limitation
-       IF ( rotation_angle /= 0.0 )  THEN
-          message_string = 'rotation of the model domain is not considered in the radiation ' //   &
-                           'model.&Using rotation_angle /= 0.0 is not allowed in combination ' //  &
-                           'with the radiation model at the moment!'
-          CALL message( 'check_parameters', 'PA0675', 1, 2, 0, 6, 0 )
        ENDIF
 
@@ -5847 +5841 @@
 
      REAL(wp)                          ::  sun_direct_factor  !< factor for direct normal radiation from direct horizontal
+     REAL(wp)                          ::  sin_rot            !< sine of rotation_angle
+     REAL(wp)                          ::  cos_rot            !< cosine of rotation_angle
+     REAL(wp)                          ::  solar_azim         !< solar azimuth in rotated model coordinates
 #if defined( __parallel )
      REAL(wp), DIMENSION(1:7)          ::  combine_allreduce   !< dummy array used to combine several MPI_ALLREDUCE calls
@@ -5866 +5863 @@
 
 !
-!--     prepare rotated normal vectors and irradiance factor
+!--  prepare rotated normal vectors and irradiance factor
+     sin_rot = SIN( rotation_angle * pi / 180.0_wp )
+     cos_rot = COS( rotation_angle * pi / 180.0_wp )
      vnorm(1,:) = kdir(:)
      vnorm(2,:) = jdir(:)
      vnorm(3,:) = idir(:)
-     mrot(1, :) = (/ 1._wp,  0._wp,      0._wp      /)
-     mrot(2, :) = (/ 0._wp,  COS(rotation_angle), SIN(rotation_angle) /)
-     mrot(3, :) = (/ 0._wp, -SIN(rotation_angle), COS(rotation_angle) /)
+     mrot(1, :) = (/ 1._wp,  0._wp,   0._wp   /)
+     mrot(2, :) = (/ 0._wp,  cos_rot, sin_rot /)
+     mrot(3, :) = (/ 0._wp, -sin_rot, cos_rot /)
      sunorig = (/ cos_zenith, sun_dir_lat, sun_dir_lon /)
      sunorig = MATMUL(mrot, sunorig)
@@ -6088 +6087 @@
 !
 !--     Identify solar direction vector (discretized number) (1)
+        solar_azim = ATAN2(sun_dir_lon, sun_dir_lat) * (180.0_wp/pi) - rotation_angle
         j = FLOOR(ACOS(cos_zenith) / pi * raytrace_discrete_elevs)
-        i = MODULO(NINT(ATAN2(sun_dir_lon, sun_dir_lat)               &
-                        / (2.0_wp*pi) * raytrace_discrete_azims-0.5_wp, iwp), &
+        i = MODULO(NINT(solar_azim / 360.0_wp *                 &
+                        raytrace_discrete_azims - 0.5_wp, iwp), &
                    raytrace_discrete_azims)
         isd = dsidir_rev(j, i)
@@ -9169 +9169 @@
    SUBROUTINE radiation_presimulate_solar_pos
 
+      USE control_parameters,                                              &
+         ONLY:  rotation_angle
+
       IMPLICIT NONE
 
@@ -9220 +9223 @@
 
          REAL(wp), INTENT(IN) ::  time_since_reference_local  !< local time since reference
+         REAL(wp)             ::  solar_azim                  !< solar azimuth in rotated model coordinates
 !
 !--      Update apparent solar position based on modified t_s_r_p
@@ -9229 +9233 @@
 !--
 !--         Identify solar direction vector (discretized number) 1)
-            i = MODULO(NINT(ATAN2(sun_dir_lon, sun_dir_lat)               &
-                            / (2._wp*pi) * raytrace_discrete_azims-.5_wp, iwp), &
+            solar_azim = ATAN2(sun_dir_lon, sun_dir_lat) * (180.0_wp/pi) - rotation_angle
+            i = MODULO(NINT(solar_azim / 360.0_wp *                &
+                            raytrace_discrete_azims - .5_wp, iwp), &
                        raytrace_discrete_azims)
             j = FLOOR(ACOS(cos_zenith) / pi * raytrace_discrete_elevs)