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

Timestamp:

Oct 30, 2019 4:01:14 PM (4 years ago)

Author:

resler

Message:

Merge branch resler into trunk

Location:

palm/trunk

Files:

• . (modified) (1 prop)
• SOURCE (modified) (1 prop)
• SOURCE/radiation_model_mod.f90 (modified) (9 diffs, 1 prop)

palm/trunk/SOURCE/radiation_model_mod.f90

@@ -28 +28 @@
 ! -----------------
 ! $Id$
+! - Fix wrong treating of time_rad during interpolation in radiation_model_mod
+! - Fix wrong checks of time_rad from dynamic driver in radiation_model_mod
+! - Add new directional model of human body for MRT: ellipsoid
+! 
+! 4271 2019-10-23 10:46:41Z maronga
 ! Bugfix: missing parentheses in calculation of snow albedo
 ! 
@@ -670 +675 @@
     LOGICAL                                        ::  mrt_skip_roof = .TRUE.             !< do not calculate MRT above roof surfaces
     LOGICAL                                        ::  mrt_include_sw = .TRUE.            !< should MRT calculation include SW radiation as well?
-    LOGICAL                                        ::  mrt_geom_human = .TRUE.            !< MRT direction weights simulate human body instead of a sphere
+    INTEGER(wp)                                    ::  mrt_geom = 1                       !< method for MRT direction weights simulating a sphere or a human body
+    REAL(wp), DIMENSION(2)                         ::  mrt_geom_params = (/ .12_wp, .88_wp /)   !< parameters for the selected method
     INTEGER(iwp)                                   ::  nrefsteps = 3                      !< number of reflection steps to perform
     REAL(wp), PARAMETER                            ::  ext_coef = 0.6_wp                  !< extinction coefficient (a.k.a. alpha)
@@ -967 +973 @@
 !-- Public variables and constants / NEEDS SORTING
     PUBLIC albedo, albedo_type, decl_1, decl_2, decl_3, dots_rad, dt_radiation,&
-           emissivity, force_radiation_call, lat, lon, mrt_geom_human,         &
+           emissivity, force_radiation_call, lat, lon, mrt_geom,               &
+           mrt_geom_params,                                                    &
            mrt_include_sw, mrt_nlevels, mrtbl, mrtinsw, mrtinlw, nmrtbl,       &
            rad_net_av, radiation, radiation_scheme, rad_lw_in,                 &
@@ -2949 +2956 @@
           CALL get_date_time( 0.0_wp, second_of_day=second_of_day )
 
-          IF ( ABS( time_rad_f%var1d(0) - second_of_day ) > 1E-6_wp )  THEN
-             message_string = 'External radiation forcing: first point in ' // &
-                              'time is /= origin_date_time.'
-             CALL message( 'radiation_init', 'PA0313', 1, 2, 0, 6, 0 )
-          ENDIF
-          
-          IF ( end_time - spinup_time > time_rad_f%var1d(ntime-1)              &
-                                      - second_of_day )  THEN
+          IF ( end_time - spinup_time > time_rad_f%var1d(ntime-1) )  THEN
              message_string = 'External radiation forcing does not cover ' //  &
                               'the entire simulation time.'
@@ -3147 +3147 @@
           CALL get_date_time( 0.0_wp, second_of_day=second_of_day_init )
           t = 0
-          DO WHILE ( time_rad_f%var1d(t) <=                                    &
-                     time_since_reference_point + second_of_day_init )
+          DO WHILE ( time_rad_f%var1d(t) <= time_since_reference_point )
              t = t + 1
           ENDDO
@@ -3154 +3153 @@
           tm = MAX( t-1, 0 )
           
-          fac_dt = ( time_since_reference_point + second_of_day_init           &
+          fac_dt = ( time_since_reference_point                                &
                    - time_rad_f%var1d(tm) + dt_3d )                            &
                  / ( time_rad_f%var1d(t)  - time_rad_f%var1d(tm) )
@@ -3844 +3843 @@
                                   constant_albedo, dt_radiation, emissivity,    &
                                   lw_radiation, max_raytracing_dist,            &
-                                  min_irrf_value, mrt_geom_human,               &
+                                  min_irrf_value, mrt_geom, mrt_geom_params,    &
                                   mrt_include_sw, mrt_nlevels,                  &
                                   mrt_skip_roof, net_radiation, nrefsteps,      &
@@ -3860 +3859 @@
                                   constant_albedo, dt_radiation, emissivity,    &
                                   lw_radiation, max_raytracing_dist,            &
-                                  min_irrf_value, mrt_geom_human,               &
+                                  min_irrf_value, mrt_geom, mrt_geom_params,    &
                                   mrt_include_sw, mrt_nlevels,                  &
                                   mrt_skip_roof, net_radiation, nrefsteps,      &
@@ -7720 +7719 @@
            zbdry(:) = (/( zn0+REAL(izn,wp)*zns, izn=0, nzn )/)
            vffrac0(:) = (COS(zbdry(0:nzn-1)) - COS(zbdry(1:nzn))) / 2._wp / REAL(naz, wp)
-           !
-           !--Modify direction weights to simulate human body (lower weight for top-down)
-           IF ( mrt_geom_human )  THEN
-              vffrac0(:) = vffrac0(:) * MAX(0._wp, SIN(zcent(:))*0.88_wp + COS(zcent(:))*0.12_wp)
+!
+!--        Modify direction weights to simulate human body (lower weight for
+!--        irradiance from zenith, higher from sides) depending on selection.
+!--        For mrt_geom=0, no weighting is done (simulates spherical globe
+!--        thermometer).
+           SELECT CASE ( mrt_geom )
+
+           CASE ( 1 )
+              vffrac0(:) = vffrac0(:) * MAX(0._wp, SIN(zcent(:))*mrt_geom_params(2) &
+                                                   + COS(zcent(:))*mrt_geom_params(1))
               vffrac0(:) = vffrac0(:) / (SUM(vffrac0) * REAL(naz, wp))
-           ENDIF
+
+           CASE ( 2 )
+              vffrac0(:) = vffrac0(:)                                          &
+                           * SQRT( ( mrt_geom_params(1) * COS(zcent(:)) ) ** 2 &
+                                 + ( mrt_geom_params(2) * SIN(zcent(:)) ) ** 2 )
+              vffrac0(:) = vffrac0(:) / (SUM(vffrac0) * REAL(naz, wp))
+
+           END SELECT
 
            DO  imrt = 1, nmrtbl