!> @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