source: palm/trunk/SOURCE/biometeorology_mod.f90 @ 3358

!> @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 <http://www.gnu.org/licenses/>.
!
! 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 <resler@cs.cas.cz>
! @author Dominik Froehlich <dominik.froehlich@dwd.de>
!
!------------------------------------------------------------------------------!

MODULE biometeorology_mod
!
!-- Load required variables from existing modules
    USE arrays_3d,                                                             &
        ONLY:  hyp, p, pt, q, u, v, w

    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) + 0.622_wp ) / 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) + 0.622_wp ) / 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