Changeset 3464 for palm/trunk

Timestamp:

Oct 30, 2018 6:08:55 PM (5 years ago)

Author:

kanani

Message:

from branch resler@3462: add MRT shaping function (radiation_model_mod), use basic constants (biometeorology_mod), adjust precision to wp (biometeorology_pt_mod)

Location:

palm/trunk/SOURCE

Files:

• biometeorology_mod.f90 (modified) (10 diffs)
• biometeorology_pt_mod.f90 (modified) (20 diffs)
• netcdf_data_input_mod.f90 (modified) (1 diff)
• netcdf_interface_mod.f90 (modified) (2 diffs)
• radiation_model_mod.f90 (modified) (6 diffs)

palm/trunk/SOURCE/biometeorology_mod.f90

@@ -27 +27 @@
 ! -----------------
 ! $Id$
+! From branch resler@3462, pavelkrc:
+! make use of basic_constants_and_equations_mod
+! 
+! 3448 2018-10-29 18:14:31Z kanani
 ! Initial revision
 !
@@ -60 +64 @@
 
     USE basic_constants_and_equations_mod,                                     &
-       ONLY:  magnus
+       ONLY:  degc_to_k, magnus, sigma_sb
 
     USE biometeorology_ipt_mod
@@ -113 +117 @@
     REAL ( wp )    ::  biom_output_height  !< height output is calculated in m
     REAL ( wp )    ::  time_biom_results   !< the time, the last set of biom results have been calculated for
-    REAL ( wp ), PARAMETER ::  cels_offs = 273.15_wp  !< Kelvin-Celsius offset (K)
-    REAL ( wp ), PARAMETER ::  sigma_sb  = 5.67037321E-8_wp  !< Stefan-Boltzmann constant
     REAL ( wp ), PARAMETER ::  human_absorb = 0.7_wp  !< SW absorbtivity of a human body (Fanger 1972)
     REAL ( wp ), PARAMETER ::  human_emiss = 0.97_wp  !< LW emissivity of a human body after (Fanger 1972)
@@ -313 +315 @@
                        biom_mrt_av(:) = biom_mrt_av(:) + &
                           ((human_absorb*mrtinsw(:) + human_emiss*mrtinlw(:))  &
-                          / (human_emiss*sigma_sb)) ** .25_wp - cels_offs
+                          / (human_emiss*sigma_sb)) ** .25_wp - degc_to_k
                     ELSE
                        biom_mrt_av(:) = biom_mrt_av(:) + &
                           (human_emiss * mrtinlw(:) / sigma_sb) ** .25_wp      &
-                          - cels_offs
+                          - degc_to_k
                     ENDIF
                  ENDIF
@@ -1046 +1048 @@
 
 !-- Compute mean radiant temperature
-    tmrt = ( nrfd / (human_emiss * sigma_sb) )**0.25_wp - cels_offs
+    tmrt = ( nrfd / (human_emiss * sigma_sb) )**0.25_wp - degc_to_k
 
  END SUBROUTINE calculate_tmrt_6_directions
@@ -1122 +1124 @@
 
     sw_side = sw_in * 0.125_wp  ! distribute half of upper sw_in to the 4 sides
-    lw_air = ( sigma_sb * 0.95_wp * ( ta + cels_offs )**4 )
+    lw_air = ( sigma_sb * 0.95_wp * ( ta + degc_to_k )**4 )
 
 !-- Compute mean radiation flux density absorbed by rotational symetric human
@@ -1133 +1135 @@
 
 !-- Compute mean radiant temperature
-    tmrt = ( nrfd / (human_emiss * sigma_sb) )**0.25_wp - cels_offs
+    tmrt = ( nrfd / (human_emiss * sigma_sb) )**0.25_wp - degc_to_k
 
  END SUBROUTINE calculate_tmrt_2_directions
@@ -1221 +1223 @@
     IF ( average_input ) THEN
 !--    Calculate ta from Tp assuming dry adiabatic laps rate
-       ta = pt_av(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - cels_offs
+       ta = pt_av(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
 
        vp = 0.034_wp
@@ -1233 +1235 @@
     ELSE
 !-- Calculate ta from Tp assuming dry adiabatic laps rate
-       ta = pt(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - cels_offs
+       ta = pt(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
 
        vp = q(k, j, i)
@@ -1249 +1251 @@
 !-- The magnus formula is limited to temperatures up to 333.15 K to
 !   avoid negative values of vp_sat
-    vp_sat = 0.01_wp * magnus( MIN( ta + cels_offs, 333.15_wp ) )
+    vp_sat = 0.01_wp * magnus( MIN( ta + degc_to_k, 333.15_wp ) )
     vp  = vp * pair / ( vp + 0.622_wp )
     IF ( vp > vp_sat ) vp = vp_sat

palm/trunk/SOURCE/biometeorology_pt_mod.f90

@@ -20 +20 @@
 !--------------------------------------------------------------------------------!
 !
-! Current revisions: 001
-! -----------------
+! Current revisions:
+! ------------------
 ! 
 ! 
-! Former revisions: 001
+! Former revisions:
 ! -----------------
 ! $Id$
-! Initial revision 001
+! From branch resler@3462, pavelkrc:
+! replace specific precision DLOG, DMOD intrinsics with generic precision for
+! compatibility with 'wp'
+! 
+! 
+! Initial revision
 !
 ! 
@@ -116 +121 @@
 ! Description:
 ! ------------
-!> PT_BASIC.F90   Version of perceived temperature (PT, °C) for
+!> PT_BASIC.F90   Version of perceived temperature (PT, ï¿œC) for
 !> - standard measured/predicted meteorological values and TMRT 
 !> as input;
@@ -131 +136 @@
 
 !-- Type of input of the argument list
-    REAL(wp), INTENT ( IN )  :: tt2m   !< Local air temperature (°C)
+    REAL(wp), INTENT ( IN )  :: tt2m   !< Local air temperature (ï¿œC)
     REAL(wp), INTENT ( IN )  :: el2m   !< Local vapour pressure (hPa)
-    REAL(wp), INTENT ( IN )  :: tmrt   !< Local mean radiant temperature (°C)
+    REAL(wp), INTENT ( IN )  :: tmrt   !< Local mean radiant temperature (ï¿œC)
     REAL(wp), INTENT ( IN )  :: vau1m  !< Local wind velocitry (m/s)
     REAL(wp), INTENT ( IN )  :: pb     !< Local barometric air pressure (hPa)
 
 !-- Type of output of the argument list
-    REAL(wp), INTENT ( OUT ) :: pt_basic  !< Perceived temperature (°C)
+    REAL(wp), INTENT ( OUT ) :: pt_basic  !< Perceived temperature (ï¿œC)
     REAL(wp), INTENT ( OUT ) :: clo       !< Clothing index (dimensionless)
 
@@ -301 +306 @@
 ! ------------
 !> The SUBROUTINE calculates the saturation water vapour pressure 
-!> (hPa = hecto Pascal) for a given temperature tt (°C).
+!> (hPa = hecto Pascal) for a given temperature tt (ï¿œC).
 !> For example, tt can be the air temperature or the dew point temperature.
 !------------------------------------------------------------------------------!
@@ -308 +313 @@
     IMPLICIT NONE
 
-    REAL(wp), INTENT ( IN )  ::  tt   !< ambient air temperature (°C)
+    REAL(wp), INTENT ( IN )  ::  tt   !< ambient air temperature (ï¿œC)
     REAL(wp), INTENT ( OUT ) ::  p_st !< saturation water vapour pressure (hPa)
 
@@ -316 +321 @@
 
     IF ( tt < 0._wp ) THEN
-!--    tt  < 0 (°C): saturation water vapour pressure over ice
+!--    tt  < 0 (ï¿œC): saturation water vapour pressure over ice
        b = 17.84362_wp
        c = 245.425_wp
     ELSE
-!--    tt >= 0 (°C): saturation water vapour pressure over water
+!--    tt >= 0 (ï¿œC): saturation water vapour pressure over water
        b = 17.08085_wp
        c = 234.175_wp
@@ -346 +351 @@
 
 !-- Input variables of argument list:
-    REAL(wp), INTENT ( IN )  :: tt2m     !< Ambient temperature (°C)
-    REAL(wp), INTENT ( IN )  :: tmrt     !< Mean radiant temperature (°C)
+    REAL(wp), INTENT ( IN )  :: tt2m     !< Ambient temperature (ï¿œC)
+    REAL(wp), INTENT ( IN )  :: tmrt     !< Mean radiant temperature (ï¿œC)
     REAL(wp), INTENT ( IN )  :: el2m     !< Water vapour pressure (hPa)
     REAL(wp), INTENT ( IN )  :: vau1m    !< Wind speed (m/s) 1 m above ground
@@ -362 +367 @@
 ! Output variables of argument list:
     REAL(wp), INTENT ( OUT ) :: pmva     !< 0 (set to zero, because clo is evaluated for comfort)
-    REAL(wp), INTENT ( OUT ) :: top      !< Operative temperature (°C) at found root of Fanger's PMV
+    REAL(wp), INTENT ( OUT ) :: top      !< Operative temperature (ï¿œC) at found root of Fanger's PMV
     REAL(wp), INTENT ( OUT ) :: clo_res  !< Resulting clothing insulation value (clo)
     INTEGER(iwp), INTENT ( OUT ) :: nerr !< Error status / quality flag
@@ -488 +493 @@
     REAL(wp), INTENT ( IN ) ::  pmv   !< Fangers predicted mean vote (dimensionless)
     REAL(wp), INTENT ( IN ) ::  clo   !< clothing insulation index (clo)
-    REAL(wp), INTENT ( OUT ) ::  pt   !< pt (°C) corresponding to given PMV / clo
+    REAL(wp), INTENT ( OUT ) ::  pt   !< pt (ï¿œC) corresponding to given PMV / clo
 
     IF ( pmv <= -0.11_wp ) THEN
@@ -517 +522 @@
 
 !-- Input variables of argument list:
-    REAL(wp), INTENT ( IN ) ::  tt       !< Ambient air temperature (°C)
-    REAL(wp), INTENT ( IN ) ::  tmrt     !< Mean radiant temperature (°C)
+    REAL(wp), INTENT ( IN ) ::  tt       !< Ambient air temperature (ï¿œC)
+    REAL(wp), INTENT ( IN ) ::  tmrt     !< Mean radiant temperature (ï¿œC)
     REAL(wp), INTENT ( IN ) ::  pa       !< Water vapour pressure (hPa)
     REAL(wp), INTENT ( IN ) ::  pb       !< Barometric pressure (hPa) at site
@@ -530 +535 @@
 !            to Fanger corresponding to meteorological (tt,tmrt,pa,vau,pb) 
 !            and individual variables (clo, actlev, eta)
-    REAL(wp), INTENT ( OUT ) ::  top      !< operative temperature (°C)
+    REAL(wp), INTENT ( OUT ) ::  top      !< operative temperature (ï¿œC)
 
 !-- Internal variables
@@ -638 +643 @@
 !-- Input variables of argument list:
     REAL(wp),     INTENT ( IN )  :: pmva     !< Actual Predicted Mean Vote (PMV) according to Fanger
-    REAL(wp),     INTENT ( IN )  :: tt2m     !< Ambient temperature (°C) at screen level
+    REAL(wp),     INTENT ( IN )  :: tt2m     !< Ambient temperature (ï¿œC) at screen level
     REAL(wp),     INTENT ( IN )  :: el2m     !< Water vapour pressure (hPa) at screen level
     REAL(wp),     INTENT ( IN )  :: svp_tt   !< Saturation water vapour pressure (hPa) at tt2m
-    REAL(wp),     INTENT ( IN )  :: tmrt     !< Mean radiant temperature (°C) at screen level
+    REAL(wp),     INTENT ( IN )  :: tmrt     !< Mean radiant temperature (ï¿œC) at screen level
     REAL(wp),     INTENT ( IN )  :: vau1m    !< Wind speed (m/s) 1 m above ground
 
@@ -724 +729 @@
        nerr = -3_iwp
        IF ( el2m < p10 ) THEN
-!--       Due to conditions of regressíon: r.H. >= 5 %
+!--       Due to conditions of regressï¿œon: r.H. >= 5 %
           pa = p10
        ELSE
-!--       Due to conditions of regressíon: r.H. <= 95 %
+!--       Due to conditions of regressï¿œon: r.H. <= 95 %
           pa = p95
        ENDIF
@@ -733 +738 @@
     IF ( pa > 0._wp ) THEN
 !--    Natural logarithm of pa
-       apa = DLOG ( pa )
+       apa = LOG ( pa )
     ELSE
        apa = -5._wp
@@ -741 +746 @@
     pa_p50   = 0.5_wp * svp_tt
     IF ( pa_p50 > 0._wp .AND. pa > 0._wp ) THEN
-       dapa   = apa - DLOG ( pa_p50 )
+       dapa   = apa - LOG ( pa_p50 )
        pa_p50 = pa / pa_p50
     ELSE
@@ -765 +770 @@
        RETURN
     ENDIF
-    gew = DMOD ( pmv, 1._wp )
+    gew = MOD ( pmv, 1._wp )
     IF ( gew < 0._wp ) gew = 0._wp
     IF ( nreg > 5_iwp ) THEN
@@ -834 +839 @@
 !> velocity:
 !> - defined only under warm conditions: gtc based on 0.5 (clo) and
-!>                                       gtc > 16.826 (°C)
+!>                                       gtc > 16.826 (ï¿œC)
 !> - undefined: sultr_res set at +99.
 !------------------------------------------------------------------------------!
@@ -904 +909 @@
 !-- Type of input arguments
     REAL(wp), INTENT ( IN ) :: pmva      !< Fanger's classical predicted mean vote
-    REAL(wp), INTENT ( IN ) :: tt2m      !< Air temperature (°C) 2 m above ground
+    REAL(wp), INTENT ( IN ) :: tt2m      !< Air temperature (ï¿œC) 2 m above ground
     REAL(wp), INTENT ( IN ) :: vau1m     !< Relative wind velocity 1 m above ground (m/s)
-    REAL(wp), INTENT ( IN ) :: tmrt      !< Mean radiant temperature (°C)
+    REAL(wp), INTENT ( IN ) :: tmrt      !< Mean radiant temperature (ï¿œC)
 
 !-- Type of output argument
@@ -1036 +1041 @@
 !> (neither body cooling nor body heating). It is related to the Klima- 
 !> Michel activity level (134.682 W/m2). IREQ_neutral is only defined 
-!> for gt < 10 (°C)
+!> for gt < 10 (ï¿œC)
 !------------------------------------------------------------------------------!
  REAL(wp) FUNCTION ireq_neutral( gt, ireq_minimal, nerr )
@@ -1059 +1064 @@
     ireq_neutral = 1.62_wp - 0.0564_wp * top02
 
-!-- Regression only defined for gt <= 10 (°C)
+!-- Regression only defined for gt <= 10 (ï¿œC)
     IF ( ireq_neutral < 0.5_wp ) THEN
        IF ( ireq_neutral < 0.48_wp ) THEN

palm/trunk/SOURCE/netcdf_data_input_mod.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! Define coordinate reference system (crs) and read from input dataset
-! Revise default values for reference coordinates
+! 
 ! 
 ! Former revisions:
 ! -----------------
 ! $Id$
+! Define coordinate reference system (crs) and read from input dataset
+! Revise default values for reference coordinates
+! 
+! 3459 2018-10-30 15:04:11Z gronemeier
 ! from chemistry branch r3443, banzhafs, Russo:
 ! Uncommented lines on dimension of surface_fractions

palm/trunk/SOURCE/netcdf_interface_mod.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
+! 
+! 
+! Former revisions:
+! -----------------
+! $Id$
 ! - Add variable crs to output files
 ! - Add long_name to UTM coordinates
@@ -25 +30 @@
 !   only written if parallel output is used.
 ! 
-! Former revisions:
-! -----------------
-! $Id$
+! 3459 2018-10-30 15:04:11Z gronemeier
 ! Adjustment of biometeorology calls
 ! 

palm/trunk/SOURCE/radiation_model_mod.f90

@@ -28 +28 @@
 ! -----------------
 ! $Id$
+! From branch resler@3462, pavelkrc:
+! add MRT shaping function for human
+! 
+! 3449 2018-10-29 19:36:56Z suehring
 ! New RTM version 3.0: (Pavel Krc, Jaroslav Resler, ICS, Prague)
 !   - Interaction of plant canopy with LW radiation
@@ -902 +906 @@
     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(iwp)                                   ::  nrefsteps = 3                      !< number of reflection steps to perform
     REAL(wp), PARAMETER                            ::  ext_coef = 0.6_wp                  !< extinction coefficient (a.k.a. alpha)
@@ -1172 +1177 @@
 !-- 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,                         &
+           emissivity, force_radiation_call, lat, lon, mrt_geom_human,         &
            mrt_include_sw, mrt_nlevels, mrtbl, mrtinsw, mrtinlw, nmrtbl,       &
            rad_net_av, radiation, radiation_scheme, rad_lw_in,                 &
@@ -2942 +2947 @@
                                   constant_albedo, dt_radiation, emissivity,    &
                                   lw_radiation, max_raytracing_dist,            &
-                                  min_irrf_value, mrt_include_sw, mrt_nlevels,  &
+                                  min_irrf_value, mrt_geom_human,               &
+                                  mrt_include_sw, mrt_nlevels,                  &
                                   mrt_skip_roof, net_radiation, nrefsteps,      &
                                   plant_lw_interact, rad_angular_discretization,&
@@ -2957 +2963 @@
                                   constant_albedo, dt_radiation, emissivity,    &
                                   lw_radiation, max_raytracing_dist,            &
-                                  min_irrf_value, mrt_include_sw, mrt_nlevels,  &
+                                  min_irrf_value, mrt_geom_human,               &
+                                  mrt_include_sw, mrt_nlevels,                  &
                                   mrt_skip_roof, net_radiation, nrefsteps,      &
                                   plant_lw_interact, rad_angular_discretization,&
@@ -6553 +6560 @@
            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)
+              vffrac0(:) = vffrac0(:) / (SUM(vffrac0) * REAL(naz, wp))
+           ENDIF
 
            DO  imrt = 1, nmrtbl