!> @file biometeorology.f90
!------------------------------------------------------------------------------!
! This file is part of PALM-4U.
!
! PALM-4U is free software: you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software
! Foundation, either version 3 of the License, or (at your option) any later
! version.
!
! PALM-4U is distributed in the hope that it will be useful, but WITHOUT ANY
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with
! PALM. If not, see .
!
! Copyright 2018, Institute of Computer Science,Academy of Sciences, Prague
!
! Calculation of PET:
! Copyright 2016, Deutscher Wetterdienst (DWD) /
! German Meteorological Service (DWD)
!------------------------------------------------------------------------------!
!
! Current revisions:
! -----------------
!
!
! Former revisions:
! -----------------
! $Id$
! Initial revision
!
!
!
! Authors:
! --------
! @author Jaroslav Resler
! @author Dominik Froehlich
!
!------------------------------------------------------------------------------!
MODULE biometeorology_mod
!
!-- Load required variables from existing modules
USE arrays_3d, &
ONLY: hyp, p, pt, q, u, v, w
USE basic_constants_and_equations_mod, &
ONLY: rd_d_rv
USE kinds !< to set precision of INTEGER and REAL arrays according to PALM
USE radiation_model_mod, &
ONLY: ix, iy, iz, id, mrt_nlevels, mrt_include_sw, &
mrtinsw, mrtinlw, mrtbl, nmrtbl, radiation
IMPLICIT NONE
REAL(wp), DIMENSION(:), ALLOCATABLE :: bio_mrt !< biometeorology mean radiant temperature for each MRT box
REAL(wp), DIMENSION(:), ALLOCATABLE :: bio_mrt_av !< time average mean radiant temperature for each MRT box
REAL(wp), DIMENSION(:), ALLOCATABLE :: bio_pet !< PhysiologiCALLy Equivalent Temperature (PET) for each MRT box
REAL(wp), DIMENSION(:), ALLOCATABLE :: bio_pet_av !< time average PhysiologiCALLy Equivalent Temperature (PET) for each MRT box
REAL(wp), PARAMETER :: sigma_sb = 5.67037321E-8_wp, & !< Stefan-Boltzmann constant
t_zero = -273.15_wp, & !< temperature 0K in Celsius
human_absorb = 0.7_wp, & !< SW absorbtivity of human body
human_emiss = 0.97_wp !< emissivity of human body (Lindberg 2008)
INTERFACE biometeorology_3d_data_averaging
MODULE PROCEDURE biometeorology_3d_data_averaging
END INTERFACE biometeorology_3d_data_averaging
INTERFACE biometeorology_check_data_output
MODULE PROCEDURE biometeorology_check_data_output
END INTERFACE biometeorology_check_data_output
INTERFACE biometeorology_data_output_3d
MODULE PROCEDURE biometeorology_data_output_3d
END INTERFACE biometeorology_data_output_3d
INTERFACE biometeorology_define_netcdf_grid
MODULE PROCEDURE biometeorology_define_netcdf_grid
END INTERFACE biometeorology_define_netcdf_grid
SAVE
PRIVATE
!
!-- Public functions
PUBLIC biometeorology_3d_data_averaging, biometeorology_check_data_output, &
biometeorology_data_output_3d, biometeorology_define_netcdf_grid
!
!-- Public variables and constants / NEEDS SORTING
! PUBLIC
CONTAINS
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Check data output for biometeorology model
!------------------------------------------------------------------------------!
SUBROUTINE biometeorology_check_data_output( var, unit, i, ilen, k )
USE control_parameters, &
ONLY: data_output, message_string
IMPLICIT NONE
CHARACTER (LEN=*) :: unit !<
CHARACTER (LEN=*) :: var !<
INTEGER(iwp) :: i
INTEGER(iwp) :: ilen
INTEGER(iwp) :: k
SELECT CASE ( TRIM( var ) )
CASE ( 'bio_mrt', 'bio_pet' )
IF ( .NOT. radiation ) THEN
message_string = 'output of "' // TRIM( var ) // '" require'&
// 's radiation = .TRUE.'
CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
ENDIF
IF ( mrt_nlevels == 0 ) THEN
message_string = 'output of "' // TRIM( var ) // '" require'&
// 's mrt_nlevels > 0'
CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
ENDIF
IF ( TRIM( var ) == 'bio_mrt' ) THEN
unit = 'K'
ELSE
unit = 'W m-2'
ENDIF
CASE DEFAULT
unit = 'illegal'
END SELECT
END SUBROUTINE biometeorology_check_data_output
!------------------------------------------------------------------------------!
!
! Description:
! ------------
!> Subroutine for averaging 3D data
!------------------------------------------------------------------------------!
SUBROUTINE biometeorology_3d_data_averaging( mode, variable )
USE control_parameters
USE indices
USE kinds
IMPLICIT NONE
CHARACTER (LEN=*) :: mode !<
CHARACTER (LEN=*) :: variable !<
INTEGER(iwp) :: i !<
INTEGER(iwp) :: j !<
INTEGER(iwp) :: k !<
INTEGER(iwp) :: l, m !< index of current surface element
REAL(wp) :: mrt, pet
IF ( mode == 'allocate' ) THEN
SELECT CASE ( TRIM( variable ) )
CASE ( 'bio_mrt' )
IF ( .NOT. ALLOCATED( bio_mrt_av ) ) THEN
ALLOCATE( bio_mrt_av(nmrtbl) )
ENDIF
bio_mrt_av = 0.0_wp
CASE ( 'bio_pet' )
IF ( .NOT. ALLOCATED( bio_pet_av ) ) THEN
ALLOCATE( bio_pet_av(nmrtbl) )
ENDIF
bio_pet_av = 0.0_wp
CASE DEFAULT
CONTINUE
END SELECT
ELSEIF ( mode == 'sum' ) THEN
SELECT CASE ( TRIM( variable ) )
CASE ( 'bio_mrt' )
IF ( ALLOCATED( bio_mrt_av ) ) THEN
IF ( nmrtbl > 0 ) THEN
IF ( mrt_include_sw ) THEN
bio_mrt_av(:) = bio_mrt_av(:) + ((human_absorb*mrtinsw(:) &
+ human_emiss*mrtinlw(:)) / (human_emiss*sigma_sb)) ** .25_wp
ELSE
bio_mrt_av(:) = bio_mrt_av(:) + &
(human_emiss * mrtinlw(:) / sigma_sb) ** .25_wp
ENDIF
ENDIF
ENDIF
CASE ( 'bio_pet' )
IF ( ALLOCATED( bio_pet_av ) ) THEN
DO l = 1, nmrtbl
i = mrtbl(ix,l)
j = mrtbl(iy,l)
k = mrtbl(iz,l)
IF ( mrt_include_sw ) THEN
mrt = ((human_absorb*mrtinsw(l) + human_emiss*mrtinlw(l)) &
/ (human_emiss*sigma_sb)) ** .25_wp
ELSE
mrt = mrt + (human_emiss * mrtinlw(l) / sigma_sb) ** .25_wp
ENDIF
CALL calculate_pet_static( &
pt(k,j,i) * (hyp(k) / 100000.0_wp )**0.286_wp + t_zero, & !< Air temperature (°C)
q(k,j,i) * hyp(k) / ( q(k,j,i) + rd_d_rv ) / 100._wp, & !< Vapor pressure (hPa)
SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + &
( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + &
( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2, &
0.01_wp ) ), & !< Wind speed (at scalar gridpoint) (m/s)
mrt + t_zero, & !< Mean radiant temperature (°C)
(hyp(k) + p(k,j,i)) * 0.01_wp, & !< Air pressure (hPa)
pet ) !< output variable of PET
bio_pet_av(l) = bio_pet_av(l) + pet
ENDDO
ENDIF
CASE DEFAULT
CONTINUE
END SELECT
ELSEIF ( mode == 'average' ) THEN
SELECT CASE ( TRIM( variable ) )
CASE ( 'bio_mrt' )
IF ( ALLOCATED( bio_mrt_av ) ) THEN
bio_mrt_av(:) = bio_mrt_av(:) / REAL( average_count_3d, KIND=wp )
ENDIF
CASE ( 'bio_pet' )
IF ( ALLOCATED( bio_pet_av ) ) THEN
bio_pet_av(:) = bio_pet_av(:) / REAL( average_count_3d, KIND=wp )
ENDIF
END SELECT
ENDIF
END SUBROUTINE biometeorology_3d_data_averaging
!------------------------------------------------------------------------------!
!
! Description:
! ------------
!> Subroutine defining appropriate grid for netcdf variables.
!> It is called out from subroutine netcdf.
!------------------------------------------------------------------------------!
SUBROUTINE biometeorology_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
IMPLICIT NONE
CHARACTER (LEN=*), INTENT(IN) :: var !<
LOGICAL, INTENT(OUT) :: found !<
CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<
CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<
CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<
found = .TRUE.
!
!-- Check for the grid
SELECT CASE ( TRIM( var ) )
CASE ( 'bio_mrt', 'bio_pet' )
grid_x = 'x'
grid_y = 'y'
grid_z = 'zu'
CASE DEFAULT
found = .FALSE.
grid_x = 'none'
grid_y = 'none'
grid_z = 'none'
END SELECT
END SUBROUTINE biometeorology_define_netcdf_grid
!------------------------------------------------------------------------------!
!
! Description:
! ------------
!> Subroutine defining 3D output variables
!------------------------------------------------------------------------------!
SUBROUTINE biometeorology_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
USE indices
USE kinds
IMPLICIT NONE
CHARACTER (LEN=*) :: variable !<
INTEGER(iwp) :: av !<
INTEGER(iwp) :: i, j, k, l !<
INTEGER(iwp) :: nzb_do !<
INTEGER(iwp) :: nzt_do !<
LOGICAL :: found !<
REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute
REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !<
REAL(wp) :: mrt, pet
found = .TRUE.
SELECT CASE ( TRIM( variable ) )
CASE ( 'bio_mrt' )
local_pf = REAL( fill_value, KIND = wp )
DO l = 1, nmrtbl
i = mrtbl(ix,l)
j = mrtbl(iy,l)
k = mrtbl(iz,l)
IF ( mrt_include_sw ) THEN
mrt = ((human_absorb*mrtinsw(l) + human_emiss*mrtinlw(l)) &
/ (human_emiss*sigma_sb)) ** .25_wp
ELSE
mrt = (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp
ENDIF
local_pf(i,j,k) = mrt
ENDDO
CASE ( 'bio_pet' )
local_pf = REAL( fill_value, KIND = wp )
IF ( av == 0 ) THEN
DO l = 1, nmrtbl
i = mrtbl(ix,l)
j = mrtbl(iy,l)
k = mrtbl(iz,l)
IF ( mrt_include_sw ) THEN
mrt = ((human_absorb*mrtinsw(l) + human_emiss*mrtinlw(l)) &
/ (human_emiss*sigma_sb)) ** .25_wp
ELSE
mrt = mrt + (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp
ENDIF
CALL calculate_pet_static( &
pt(k,j,i) * (hyp(k) / 100000.0_wp )**0.286_wp + t_zero, & !< Air temperature (°C)
q(k,j,i) * hyp(k) / ( q(k,j,i) + rd_d_rv ) / 100.0_wp, & !< Vapor pressure (hPa)
SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + &
( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + &
( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2, &
0.01_wp ) ), & !< Wind speed (at scalar gridpoint) (m/s)
mrt + t_zero, & !< Mean radiant temperature (°C)
(hyp(k) + p(k,j,i)) * 0.01_wp, & !< Air pressure (hPa)
pet ) !< output variable of PET
local_pf(i,j,k) = pet
ENDDO
ELSE
IF ( ALLOCATED( bio_pet_av ) ) THEN
DO l = 1, nmrtbl
i = mrtbl(ix,l)
j = mrtbl(iy,l)
k = mrtbl(iz,l)
local_pf(i,j,k) = bio_pet_av(l)
ENDDO
ENDIF
ENDIF
CASE DEFAULT
found = .FALSE.
END SELECT
END SUBROUTINE biometeorology_data_output_3d
!------------------------------------------------------------------------------!
!
! Description:
! ------------
!> PhysiologiCALLy Equivalent Temperature (PET),
! stationary (calculated based on MEMI),
! Subroutine based on PETBER vers. 1.5.1996 by P. Hoeppe
!
! Input arguments:
! ----------------
! - ta : Air temperature (°C) REAL(wp)
! - tmrt : Mean radiant temperature (°C) REAL(wp)
! - v : Wind speed (m/s) REAL(wp)
! - vpa : Vapor pressure (hPa) REAL(wp)
! - p : Air pressure (hPa) REAL(wp)
!
! Output arguments:
! ----------------
! - tx : PET (°C) REAL(wp)
!------------------------------------------------------------------------------!
SUBROUTINE calculate_pet_static( &
!-- Meteorological input
ta, vpa, v, tmrt, p, &
!-- Output variables
tx ) !, &
!-- Configure sample person (optional)
! age, mbody, ht, work, eta, icl, fcl, pos, sex )
IMPLICIT NONE
REAL(wp), INTENT( IN ) :: ta, tmrt, v, vpa, p
REAL(wp), INTENT ( OUT ) :: tx
! REAL(wp), INTENT ( in ), optional :: age, mbody, ht, work, eta, icl, fcl
REAL(wp) :: age, mbody, ht, work, eta, icl, fcl
! INTEGER(iwp), INTENT ( in ), optional :: pos, sex
INTEGER(iwp) :: pos, sex
REAL(wp) :: acl, adu, aeff, cair, cb, emcl, emsk, ere, erel, esw, &
evap, facl, feff, food, h, po, rdcl, rdsk, rob, rtv, &
sigm, vpts, &
! former intent (out)
! - tsk : Skin temperature (°C) real
! - tcl : Clothing temperature (°C) real
! - ws : real
! - wetsk : Fraction of wet skin real
tsk, tcl, wetsk
!-- Person data
! IF ( .NOT. present( age ) ) age = 35.
! IF ( .NOT. present( mbody ) ) mbody = 75.
! IF ( .NOT. present( ht ) ) ht = 1.75
! IF ( .NOT. present( work ) ) work = 80.
! IF ( .NOT. present( eta ) ) eta = 0.
! IF ( .NOT. present( icl ) ) icl = 0.9
! IF ( .NOT. present( fcl ) ) fcl = 1.15
! IF ( .NOT. present( pos ) ) pos = 1
! IF ( .NOT. present( sex ) ) sex = 1
age = 35.
mbody = 75.
ht = 1.75
work = 80.
eta = 0.
icl = 0.9
fcl = 1.15
pos = 1
sex = 1
!-- constants
po = 1013.25 !< preassure at sea level
! p = 1013.25 !< local preassure (hPa), now defined as input variable
rob = 1.06
cb = 3.64 * 1000.
food = 0.
emsk = 0.99
emcl = 0.95
evap = 2.42 * 10. ** 6.
sigm = 5.67 * 10. **(-8.)
! write(9,*) 'Call calculate_pet_static(ta=', ta, ', vpa=', vpa, &
! ', v=', v, ', tmrt=', tmrt, ', p=', p
! flush(9)
!-- call subfunctions
CALL INKOERP ( age, cair, eta, ere, erel, evap, h, ht, mbody, &
p, rtv, sex, ta, vpa, work )
CALL BERECH ( acl, adu, aeff, cair, cb, emcl, emsk, &
ere, erel, esw, evap, facl, fcl, feff, food, h, ht, icl, &
mbody, p, po, rdcl, rdsk, rob, sex, sigm, &
ta, tcl, tmrt, tsk, v, vpa, vpts, wetsk )
CALL PET ( acl, adu, aeff, cair, emcl, emsk, esw, evap, &
facl, feff, h, p, po, rdcl, rdsk, &
rtv, sigm, ta, tcl, tsk, tx, vpts, wetsk )
END SUBROUTINE calculate_pet_static
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Calculate internal energy ballance
!------------------------------------------------------------------------------!
SUBROUTINE inkoerp( age, cair, eta, ere, erel, evap, h, ht, mbody, &
& p, rtv, sex, ta, vpa, work )
REAL(wp) :: age, cair, eta, ere, erel, eres, evap, h, ht, mbody, &
& met, p, rtv, ta, tex, vpa, vpex, work
INTEGER(iwp) :: sex
IF ( sex .EQ. 1 ) THEN
met = 3.45 * mbody ** ( 3. / 4. ) * (1. + 0.004 * &
( 30. - age) + 0.010 * ( ( ht * 100. / &
( mbody ** ( 1. / 3. ) ) ) - 43.4 ) )
ELSE IF ( sex .EQ. 2 ) THEN
met = 3.19 * mbody ** ( 3. / 4. ) * ( 1. + 0.004 * &
( 30. - age ) + 0.018 * ( ( ht * 100. / ( mbody ** &
( 1. / 3. ) ) ) - 42.1 ) )
END IF
met = work + met
h = met * (1. - eta)
!-- SENSIBLE RESPIRATION ENERGY
cair = 1.01 * 1000.
tex = 0.47 * ta + 21.0
rtv = 1.44 * 10. ** (-6.) * met
eres = cair * (ta - tex) * rtv
!-- LATENT RESPIRATION ENERGY
vpex = 6.11 * 10. ** ( 7.45 * tex / ( 235. + tex ) )
erel = 0.623 * evap / p * ( vpa - vpex ) * rtv
!-- SUM OF RESULTS
ere = eres + erel
RETURN
END SUBROUTINE inkoerp
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Calculate heat gain or loss
!------------------------------------------------------------------------------!
SUBROUTINE BERECH( acl, adu, aeff, cair, cb, emcl, emsk, &
ere, erel, esw, evap, facl, fcl, feff, food, h, ht, icl, &
mbody, p, po, rdcl, rdsk, rob, sex, sigm, &
ta, tcl, tmrt, tsk, v, vpa, vpts, wetsk )
REAL(wp) :: acl, adu, aeff, c(0:10), cair, cb, cbare, &
cclo, csum, di, ed, emcl, emsk, enbal, &
enbal2, ere, erel, esw, eswdif, eswphy, eswpot, &
evap, facl, fcl, fec, feff, food, h, hc, he, ht, htcl, icl, &
mbody, p, po, r1, r2, rbare, rcl, &
rclo, rclo2, rdcl, rdsk, rob, rsum, sigm, sw, swf, swm, &
ta, tbody, tcl, tcore(1:7), tmrt, tsk, v, vb, vb1, vb2, &
vpa, vpts, wetsk, wd, wr, ws, wsum, xx, y
INTEGER(iwp) :: count1, count3, j, sex
logical :: skipIncreaseCount
wetsk = 0.
adu = 0.203 * mbody ** 0.425 * ht ** 0.725
hc = 2.67 + ( 6.5 * v ** 0.67)
hc = hc * (p /po) ** 0.55
feff = 0.725 !< Posture: 0.725 for stading, 0.696 for sitting
facl = (- 2.36 + 173.51 * icl - 100.76 * icl * icl + 19.28 &
* (icl ** 3.)) / 100.
IF ( facl .GT. 1. ) facl = 1.
rcl = ( icl / 6.45) / facl
IF ( icl .GE. 2. ) y = 1.
IF ( ( icl .GT. 0.6 ) .AND. ( icl .LT. 2. ) ) y = ( ht - 0.2 ) / ht
IF ( ( icl .LE. 0.6 ) .AND. ( icl .GT. 0.3 ) ) y = 0.5
IF ( ( icl .LE. 0.3 ) .AND. ( icl .GT. 0. ) ) y = 0.1
r2 = adu * (fcl - 1. + facl) / (2. * 3.14 * ht * y)
r1 = facl * adu / (2. * 3.14 * ht * y)
di = r2 - r1
!-- SKIN TEMPERATURE
DO j = 1,7
tsk = 34.
count1 = 0
tcl = ( ta + tmrt + tsk ) / 3.
count3 = 1
enbal2 = 0.
DO
acl = adu * facl + adu * ( fcl - 1. )
rclo2 = emcl * sigm * ( ( tcl + 273.2 )**4. - &
( tmrt + 273.2 )** 4. ) * feff
htcl = 6.28 * ht * y * di / ( rcl * LOG( r2 / r1 ) * acl )
tsk = 1. / htcl * ( hc * ( tcl - ta ) + rclo2 ) + tcl
!-- RADIATION SALDO
aeff = adu * feff
rbare = aeff * ( 1. - facl ) * emsk * sigm * &
( ( tmrt + 273.2 )** 4. - ( tsk + 273.2 )** 4. )
rclo = feff * acl * emcl * sigm * &
( ( tmrt + 273.2 )** 4. - ( tcl + 273.2 )** 4. )
rsum = rbare + rclo
!-- CONVECTION
cbare = hc * ( ta - tsk ) * adu * ( 1. - facl )
cclo = hc * ( ta - tcl ) * acl
csum = cbare + cclo
!-- CORE TEMPERATUR
c(0) = h + ere
c(1) = adu * rob * cb
c(2) = 18. - 0.5 * tsk
c(3) = 5.28 * adu * c(2)
c(4) = 0.0208 * c(1)
c(5) = 0.76075 * c(1)
c(6) = c(3) - c(5) - tsk * c(4)
c(7) = - c(0) * c(2) - tsk * c(3) + tsk * c(5)
c(8) = c(6) * c(6) - 4. * c(4) * c(7)
c(9) = 5.28 * adu - c(5) - c(4) * tsk
c(10) = c(9) * c(9) - 4. * c(4) * &
( c(5) * tsk - c(0) - 5.28 * adu * tsk )
IF ( tsk .EQ. 36. ) tsk = 36.01
tcore(7) = c(0) / ( 5.28 * adu + c(1) * 6.3 / 3600. ) + tsk
tcore(3) = c(0) / ( 5.28 * adu + ( c(1) * 6.3 / 3600. ) / &
( 1. + 0.5 * ( 34. - tsk ) ) ) + tsk
IF ( c(10) .GE. 0.) THEN
tcore(6) = ( - c(9) - c(10)**0.5 ) / ( 2. * c(4) )
tcore(1) = ( - c(9) + c(10)**0.5 ) / ( 2. * c(4) )
END IF
! 22
IF ( c(8) .GE. 0. ) THEN
tcore(2) = ( - c(6) + ABS( c(8) ) ** 0.5 ) / ( 2. * c(4) )
tcore(5) = ( - c(6) - ABS( c(8) ) ** 0.5 ) / ( 2. * c(4) )
tcore(4) = c(0) / ( 5.28 * adu + c(1) * 1. / 40. ) + tsk
END IF
! 24
!-- TRANSPIRATION
tbody = 0.1 * tsk + 0.9 * tcore(j)
swm = 304.94 * ( tbody - 36.6 ) * adu / 3600000.
vpts = 6.11 * 10.**( 7.45 * tsk / ( 235. + tsk ) )
IF ( tbody .LE. 36.6 ) swm = 0.
! swf = 0.7 * swm
IF ( sex .EQ. 1 ) sw = swm
IF ( sex .EQ. 2 ) sw = 0.7 * swm ! swf
eswphy = - sw * evap
he = 0.633 * hc / ( p * cair )
fec = 1. / ( 1. + 0.92 * hc * rcl )
eswpot = he * ( vpa - vpts ) * adu * evap * fec
wetsk = eswphy / eswpot
IF ( wetsk .GT. 1. ) wetsk = 1.
eswdif = eswphy - eswpot
IF ( eswdif .LE. 0. ) esw = eswpot
IF ( eswdif .GT. 0. ) esw = eswphy
IF ( esw .GT. 0. ) esw = 0.
!-- DIFFUSION
rdsk = 0.79 * 10. ** 7.
rdcl = 0.
ed = evap / ( rdsk + rdcl ) * adu * ( 1. - wetsk ) * ( vpa - vpts )
!-- MAX VB
vb1 = 34. - tsk
vb2 = tcore(j) - 36.6
IF ( vb2 .LT. 0. ) vb2 = 0.
IF ( vb1 .LT. 0. ) vb1 = 0.
vb = ( 6.3 + 75. * vb2 ) / ( 1. + 0.5 * vb1 )
!-- ENERGY BALLANCE
enbal = h + ed + ere + esw + csum + rsum + food
!-- CLOTHING TEMPERATURE
xx = 0.001
IF ( count1 .EQ. 0 ) xx = 1.
IF ( count1 .EQ. 1 ) xx = 0.1
IF ( count1 .EQ. 2 ) xx = 0.01
IF ( count1 .EQ. 3 ) xx = 0.001
IF ( enbal .GT. 0. ) tcl = tcl + xx
IF ( enbal .LT. 0. ) tcl = tcl - xx
skipIncreaseCount = .FALSE.
IF ( ( (enbal .LE. 0.) .AND. (enbal2 .GT. 0.) ) .OR. &
( ( enbal .GE. 0. ) .AND. ( enbal2 .LT. 0. ) ) ) THEN
skipIncreaseCount = .TRUE.
ELSE
enbal2 = enbal
count3 = count3 + 1
END IF
IF ( ( count3 .GT. 200 ) .OR. skipIncreaseCount ) THEN
IF ( count1 .LT. 3 ) THEN
count1 = count1 + 1
enbal2 = 0.
ELSE
EXIT
END IF
END IF
END DO
IF ( count1 .EQ. 3 ) THEN
SELECT CASE ( j )
CASE ( 2, 5)
IF ( .NOT. ( ( tcore(j) .GE. 36.6 ) .AND. &
( tsk .LE. 34.050 ) ) ) CYCLE
CASE ( 6, 1 )
IF ( c(10) .LT. 0. ) CYCLE
IF ( .NOT. ( ( tcore(j) .GE. 36.6 ) .AND. &
( tsk .GT. 33.850 ) ) ) CYCLE
CASE ( 3 )
IF ( .NOT. ( ( tcore(j) .LT. 36.6 ) .AND. &
( tsk .LE. 34.000 ) ) ) CYCLE
CASE ( 7 )
IF ( .NOT. ( ( tcore(j) .LT. 36.6 ) .AND. &
( tsk .GT. 34.000 ) ) ) CYCLE
CASE default ! = CASE ( 4 ), does actually nothing
END SELECT
END IF
IF ( ( j .NE. 4 ) .AND. ( vb .GE. 91. ) ) CYCLE
IF ( ( j .EQ. 4 ) .AND. ( vb .LT. 89. ) ) CYCLE
IF ( vb .GT. 90.) vb = 90.
!-- LOSSES BY WATER
ws = sw * 3600. * 1000.
IF ( ws .GT. 2000. ) ws = 2000.
wd = ed / evap * 3600. * ( -1000. )
wr = erel / evap * 3600. * ( -1000. )
wsum = ws + wr + wd
RETURN
END DO
RETURN
END SUBROUTINE Berech
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Calculate PET
!------------------------------------------------------------------------------!
SUBROUTINE PET ( acl, adu, aeff, cair, emcl, emsk, esw, evap, &
facl, feff, h, p, po, rdcl, rdsk, rtv, sigm, ta, tcl, tsk, tx, vpts, wetsk)
REAL ( wp ) :: acl, adu, aeff, cair, cbare, cclo, csum, ed, &
emcl, emsk, enbal, enbal2, ere, erel, eres, esw, evap, &
facl, feff, h, hc, p, po, rbare, rclo, rdcl, rdsk, rsum, &
rtv, sigm, ta, tcl, tex, tsk, tx, vpex, vpts, wetsk, xx
INTEGER ( iwp ) :: count1
tx = ta
enbal2 = 0.
DO count1 = 0, 3
DO
hc = 2.67 + 6.5 * 0.1 ** 0.67
hc = hc * ( p / po ) ** 0.55
!-- Radiation
aeff = adu * feff
rbare = aeff * ( 1. - facl ) * emsk * sigm * &
( ( tx + 273.2 ) ** 4. - ( tsk + 273.2 ) ** 4. )
rclo = feff * acl * emcl * sigm * &
( ( tx + 273.2 ) ** 4. - ( tcl + 273.2 ) ** 4. )
rsum = rbare + rclo
!-- Covection
cbare = hc * ( tx - tsk ) * adu * ( 1. - facl )
cclo = hc * ( tx - tcl ) * acl
csum = cbare + cclo
!-- Diffusion
ed = evap / ( rdsk + rdcl ) * adu * ( 1. - wetsk ) * ( 12. - vpts )
!-- Respiration
tex = 0.47 * tx + 21.
eres = cair * ( tx - tex ) * rtv
vpex = 6.11 * 10. ** ( 7.45 * tex / ( 235. + tex ) )
erel = 0.623 * evap / p * ( 12. - vpex ) * rtv
ere = eres + erel
!-- Energy ballance
enbal = h + ed + ere + esw + csum + rsum
!-- Iteration concerning ta
IF ( count1 .EQ. 0 ) xx = 1.
IF ( count1 .EQ. 1 ) xx = 0.1
IF ( count1 .EQ. 2 ) xx = 0.01
IF ( count1 .EQ. 3 ) xx = 0.001
IF ( enbal .GT. 0. ) tx = tx - xx
IF ( enbal .LT. 0. ) tx = tx + xx
IF ( ( enbal .LE. 0. ) .AND. ( enbal2 .GT. 0. ) ) EXIT
IF ( ( enbal .GE. 0. ) .AND. ( enbal2 .LT. 0. ) ) EXIT
enbal2 = enbal
END DO
END DO
RETURN
END SUBROUTINE PET
END MODULE biometeorology_mod