Changeset 4205 for palm


Ignore:
Timestamp:
Aug 30, 2019 1:25:00 PM (5 years ago)
Author:
suehring
Message:

plant canopy: Missing working precision + bugfix in calculation of wind speed; surface data output: Correct x,y-coordinates of vertical surfaces in netcdf output; Change definition of azimuth angle, reference is north 0 degree

Location:
palm/trunk/SOURCE
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/plant_canopy_model_mod.f90

    r4188 r4205  
    2727! -----------------
    2828! $Id$
     29! Missing working precision + bugfix in calculation of wind speed
     30!
     31! 4188 2019-08-26 14:15:47Z suehring
    2932! Minor adjustment in error number
    3033!
     
    235238     IMPLICIT NONE
    236239!--  input parameters
    237      INTEGER(iwp), INTENT(IN)          :: i, j, k, kk        !< indices of the pc gridbox
    238      REAL(wp), INTENT(IN)              :: pcbsw              !< sw radiation in gridbox (W)
    239      REAL(wp), INTENT(IN)              :: pcblw              !< lw radiation in gridbox (W)
    240      REAL(wp), INTENT(OUT)             :: pcbtr              !< transpiration rate dq/dt (kg/kg/s)
    241      REAL(wp), INTENT(OUT)             :: pcblh              !< latent heat from transpiration dT/dt (K/s)
     240     INTEGER(iwp), INTENT(IN) :: i, j, k, kk        !< indices of the pc gridbox
     241     REAL(wp), INTENT(IN)     :: pcbsw              !< sw radiation in gridbox (W)
     242     REAL(wp), INTENT(IN)     :: pcblw              !< lw radiation in gridbox (W)
     243     REAL(wp), INTENT(OUT)    :: pcbtr              !< transpiration rate dq/dt (kg/kg/s)
     244     REAL(wp), INTENT(OUT)    :: pcblh              !< latent heat from transpiration dT/dt (K/s)
    242245
    243246!--  variables and parameters for calculation of transpiration rate
    244      REAL(wp)                          :: sat_press, sat_press_d, temp, v_lad
    245      REAL(wp)                          :: d_fact, g_b, g_s, wind_speed, evapor_rate
    246      REAL(wp)                          :: f1, f2, f3, f4, vpd, rswc, e_eq, e_imp, rad
    247      REAL(wp), PARAMETER               :: gama_psychr = 66  !< psychrometric constant (Pa/K)
    248      REAL(wp), PARAMETER               :: g_s_max = 0.01     !< maximum stomatal conductivity (m/s)
    249      REAL(wp), PARAMETER               :: m_soil = 0.4_wp    !< soil water content (needs to adjust or take from LSM)
    250      REAL(wp), PARAMETER               :: m_wilt = 0.01_wp   !< wilting point soil water content (needs to adjust or take from LSM)
    251      REAL(wp), PARAMETER               :: m_sat = 0.51_wp    !< saturation soil water content (needs to adjust or take from LSM)
    252      REAL(wp), PARAMETER               :: t2_min = 0.0_wp    !< minimal temperature for calculation of f2
    253      REAL(wp), PARAMETER               :: t2_max = 40.0_wp   !< maximal temperature for calculation of f2
     247     REAL(wp)             :: sat_press, sat_press_d, temp, v_lad
     248     REAL(wp)             :: d_fact, g_b, g_s, wind_speed, evapor_rate
     249     REAL(wp)             :: f1, f2, f3, f4, vpd, rswc, e_eq, e_imp, rad
     250     REAL(wp), PARAMETER  ::  gama_psychr = 66.0_wp !< psychrometric constant (Pa/K)
     251     REAL(wp), PARAMETER  ::  g_s_max = 0.01        !< maximum stomatal conductivity (m/s)
     252     REAL(wp), PARAMETER  ::  m_soil = 0.4_wp       !< soil water content (needs to adjust or take from LSM)
     253     REAL(wp), PARAMETER  ::  m_wilt = 0.01_wp      !< wilting point soil water content (needs to adjust or take from LSM)
     254     REAL(wp), PARAMETER  ::  m_sat = 0.51_wp       !< saturation soil water content (needs to adjust or take from LSM)
     255     REAL(wp), PARAMETER  ::  t2_min = 0.0_wp       !< minimal temperature for calculation of f2
     256     REAL(wp), PARAMETER  ::  t2_max = 40.0_wp      !< maximal temperature for calculation of f2
    254257
    255258
     
    257260     temp = pt(k,j,i) * exner(k) - degc_to_k
    258261!--  Coefficient for conversion of radiation to grid to radiation to unit leaves surface
    259      v_lad = 1.0_wp / ( MAX( lad_s(kk,j,i), 1.0e-10_wp ) * dx * dy * dz(1) )
     262     v_lad = 1.0_wp / ( MAX( lad_s(kk,j,i), 1.0E-10_wp ) * dx * dy * dz(1) )
    260263!--  Magnus formula for the saturation pressure (see Ngao, Adam and Saudreau (2017) eq. 1)
    261264!--  There are updated formulas available, kept consistent with the rest of the parametrization
     
    265268!--  Wind speed
    266269     wind_speed = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 +            &
    267                         ( 0.5_wp * ( v(k,j,i) + v(k,j,i+1) ) )**2 +            &
    268                         ( 0.5_wp * ( w(k,j,i) + w(k,j,i+1) ) )**2 )
     270                        ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 +            &
     271                        ( 0.5_wp * ( w(k,j,i) + w(k-1,j,i) ) )**2 )
    269272!--  Aerodynamic conductivity (Daudet et al. (1999) eq. 14
    270273     g_b = 0.01_wp * wind_speed + 0.0071_wp
     
    273276!--  First function for calculation of stomatal conductivity (radiation dependency)
    274277!--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 17
    275      f1 = rad * (1000._wp+42.1_wp) / 1000._wp / (rad+42.1_wp)
     278     f1 = rad * (1000.0_wp+42.1_wp) / 1000.0_wp / (rad+42.1_wp)
    276279!--  Second function for calculation of stomatal conductivity (temperature dependency)
    277280!--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 21
     
    286289!--  than the coefficients from Stewart (1988) which correspond to conifer trees.
    287290     vpd = MIN(MAX(vpd,770.0_wp),3820.0_wp)
    288      f3 = -2e-4_wp * vpd + 1.154_wp
     291     f3 = -2E-4_wp * vpd + 1.154_wp
    289292!--  Fourth function for calculation of stomatal conductivity (soil moisture dependency)
    290293!--  Residual soil water content
     
    293296     rswc = ( m_sat - m_soil ) / ( m_sat - m_wilt )
    294297!--  van Wijk et al. (1998; Tree Physiology 20) eq. 5-6 (it is a reformulation of eq. 22-23 of Stewart(1988))
    295      f4 = MAX(0._wp, MIN(1.0_wp - 0.041_wp * EXP(3.2_wp * rswc), 1.0_wp - 0.041_wp))
     298     f4 = MAX(0.0_wp, MIN(1.0_wp - 0.041_wp * EXP(3.2_wp * rswc), 1.0_wp - 0.041_wp))
    296299!--  Stomatal conductivity
    297300!--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 12
    298301!--  (notation according to Ngao, Adam and Saudreau (2017) and others)
    299      g_s = g_s_max * f1 * f2 * f3 * f4 + 1.0e-10_wp
     302     g_s = g_s_max * f1 * f2 * f3 * f4 + 1.0E-10_wp
    300303!--  Decoupling factor
    301304!--  Daudet et al. (1999) eq. 6
    302      d_fact = (sat_press_d / gama_psychr + 2._wp ) /                        &
    303               (sat_press_d / gama_psychr + 2._wp + 2 * g_b / g_s )
     305     d_fact = (sat_press_d / gama_psychr + 2.0_wp ) /                        &
     306              (sat_press_d / gama_psychr + 2.0_wp + 2.0_wp * g_b / g_s )
    304307!--  Equilibrium evaporation rate
    305308!--  Daudet et al. (1999) eq. 4
  • palm/trunk/SOURCE/surface_data_output_mod.f90

    r4182 r4205  
    2525! -----------------
    2626! $Id$
     27! - Correct x,y-coordinates of vertical surfaces in netcdf output
     28! - Change definition of azimuth angle, reference is north 0 degree
     29! - zenith angle is always defined, also for vertical surfaces where it is
     30!   90 degree, while azimuth angle is only defined for vertical surfaces, not
     31!   for horizontal ones
     32!
     33! 4182 2019-08-22 15:20:23Z scharf
    2734! Corrected "Former revisions" section
    2835!
     
    10681075!
    10691076!--       For vertical surfaces, zenith angles are not defined (fill value).
    1070 !--       Azimuth angle: eastward (0), westward (180), northward (270),
    1071 !--       southward (90).
     1077!--       Azimuth angle: northward (0), eastward (90), southward (180),
     1078!--       westward (270).
     1079!--       Note, for vertical surfaces, zenith angles are 90.0_wp.
    10721080          DO  l = 0, 3
    10731081             IF ( l == 0 )  THEN
    1074                 az    = 270.0_wp
     1082                az    = 0.0_wp
    10751083                off_x = 0.5_wp
    10761084                off_y = 0.0_wp
    10771085             ELSEIF ( l == 1 )  THEN
     1086                az    = 180.0_wp
     1087                off_x = 0.5_wp
     1088                off_y = 1.0_wp
     1089             ELSEIF ( l == 2 )  THEN
    10781090                az    = 90.0_wp
    1079                 off_x = 0.5_wp
    1080                 off_y = 0.0_wp
    1081              ELSEIF ( l == 2 )  THEN
    1082                 az    = 0.0_wp
    10831091                off_x = 0.0_wp
    10841092                off_y = 0.5_wp
    10851093             ELSEIF ( l == 3 )  THEN
    1086                 az    = 180.0_wp
    1087                 off_x = 0.0_wp
     1094                az    = 270.0_wp
     1095                off_x = 1.0_wp
    10881096                off_y = 0.5_wp
    10891097             ENDIF
     
    10951103                surfaces%zs(mm)      = zu(surf_def_v(l)%k(m))
    10961104                surfaces%azimuth(mm) = az
    1097                 surfaces%zenith(mm)  = surfaces%fillvalue
     1105                surfaces%zenith(mm)  = 90.0_wp
    10981106                i                    = i + 1
    10991107                mm                   = mm + 1
     
    11051113                surfaces%zs(mm)      = zu(surf_lsm_v(l)%k(m))
    11061114                surfaces%azimuth(mm) = az
    1107                 surfaces%zenith(mm)  = surfaces%fillvalue
     1115                surfaces%zenith(mm)  = 90.0_wp
    11081116                i                    = i + 1
    11091117                mm                   = mm + 1
     
    11151123                surfaces%zs(mm)      = zu(surf_usm_v(l)%k(m))
    11161124                surfaces%azimuth(mm) = az
    1117                 surfaces%zenith(mm)  = surfaces%fillvalue
     1125                surfaces%zenith(mm)  = 90.0_wp
    11181126                i                    = i + 1
    11191127                mm                   = mm + 1
Note: See TracChangeset for help on using the changeset viewer.