#include "cfx5ext.h"
dllexport(usr_pmvppd)
      SUBROUTINE USR_PMVPPD(NLOC,NRET,NARG,RET,ARGS,CRESLT,
     &                      CZ,DZ,IZ,LZ,RZ)
      IMPLICIT NONE
CC
CD User CEL routine to calculate PMV
CC
CC --------------------
CC        Input
CC --------------------
CC
CI NLOC    : Number of locations in space
CI NARG    : Number of CEL arguments
CI NRET    : Number of return variables
CI ARGS    : Array containing argument field(s)
CC
CC --------------------
CC      Modified
CC --------------------
CC
CV Stacks
CC
CC --------------------
CC        Output
CC --------------------
CC
CO RET    : Array containing return field(s)
CO CRESLT : Result flag
CC
CC --------------------
CC       Details
CC --------------------
CC
CC Although the calculation routine calculates both PMV & PPD
CC we only need calculate PMV, since PPD can be calculated by
CC a CEL expression.
CC
CC --------------------
CC Local variables & parameters
CC --------------------
CC
CC
CC --------------------
CC Modification History
CC --------------------
CC
CC Jun 2002 CJS  Subroutine created.
CC Dec 2002 YLS  Added switches from ccl, radiation, & ISO 7730 tables.
CC Nov 2006 CJS  Tidy up. Bring up to CFX source code standard.
CC Nov 2006 YLS  Tidy up.
CC Jan 2008 NSW  Concatenated with USR_CAL_PMVPPD.
CC
C=======================================================================
C
C ---------------------------
C    Preprocessor includes
C ---------------------------
C
#include "cfd_constants.h"
C
C ------------------------------
C        Argument list
C ------------------------------
C
      CHARACTER*(4)      CRESLT
C
      INTEGER            NLOC, NARG, NRET
C
      REAL               ARGS(NLOC,NARG), RET(NLOC,NRET)
C
      INTEGER IZ(*)
      CHARACTER CZ(*)*(1)
      DOUBLE PRECISION DZ(*)
      LOGICAL LZ(*)
      REAL RZ(*)
C
C ------------------------------
C        External functions
C ------------------------------
C
C ------------------------------
C        Local Parameters
C ------------------------------
C
      CHARACTER*(*)      ROUTIN
      PARAMETER         (ROUTIN = 'USR_PMVPPD')
C
C ---------------------------
C        Local Variables
C ---------------------------
C
      INTEGER            I, IERR, JGM, JGICL, JGTRAD
      REAL               GM, W, GICL, PH2O, TAA, TRA, VAIR, PMV,
     &                   GTRAD, GINTM
C
C ---------------------------
C        Stack pointers
C ---------------------------
C
C ------------------------------
C  Directory on entry:- ANY
C ------------------------------
C
C---------------------------
C Executable Statements
C---------------------------
C
C-----------------------------------------------------------------------
C     Pick up user parameters
C-----------------------------------------------------------------------
C
      CALL USER_PEEKR('GM',IONE,GM,'STOP',CRESLT,CZ)
      CALL USER_PEEKR('WE',IONE,W,'STOP',CRESLT,CZ)
      CALL USER_PEEKR('GICL',IONE,GICL,'STOP',CRESLT,CZ)
      CALL USER_PEEKR('PH2O',IONE,PH2O,'STOP',CRESLT,CZ)
      CALL USER_PEEKR('GTRAD',IONE,GTRAD,'STOP',CRESLT,CZ)
      CALL USER_PEEKI('JGM',IONE,JGM,'STOP',CRESLT,CZ)
      CALL USER_PEEKI('JGICL',IONE,JGICL,'STOP',CRESLT,CZ)
      CALL USER_PEEKI('JGTRAD',IONE,JGTRAD,'STOP',CRESLT,CZ)
C
C ----------------------------------------------------------------------
C     Calculate PMV at all locations in this locale
C ----------------------------------------------------------------------
C
      TRA = GTRAD
      DO I = 1,NLOC
        TAA = ARGS(I,4)
        IF(JGTRAD.NE.0) THEN
         GINTM = MAX(1.0E-10, ABS(ARGS(I,5)))
         TRA = (5.5402E7*GINTM)**0.25
        ENDIF
        VAIR = SQRT (ARGS(I,1)**2 + ARGS(I,2)**2 + ARGS(I,3)**2)
C
        CALL USR_CAL_PMVPPD(GM,W,GICL,TAA,TRA,VAIR,PH2O,JGM,
     &                  JGICL,PMV,IERR)
C
        RET(I,1) = PMV
        IF (IERR .NE. 0) CRESLT = 'FAIL'
      END DO
C
C ----------------------------------------------------------------------
C     Issue error message if appropriate
C ----------------------------------------------------------------------
C
      IF (CRESLT .NE. 'GOOD') THEN
        CALL ERRMSG('Exceeded max iterations calculating PMV')
        CALL CFXSTP(ROUTIN)
      ENDIF
C
      END

      SUBROUTINE USR_CAL_PMVPPD(GM,GW,GICL,TAA,TRA,VAIR,PH2O,JGM,
     &                          JGICL,PMV,IERR)
      IMPLICIT NONE
CC
CD Calculation of the PMV Index
CC
CC --------------------
CC        Input
CC --------------------
CC
CI GM   = Metabolic rate (W/m^2)
CI GW    = External work (W/m^2) (normally 0.0)
CI GICL = Clothing thermal resistance (m^2K/W)
CI TAA  = Local gas temperature (K)
CI TRA  = Local mean radiant temperature (K)
CI VAIR = Local air speed (m/s)
CI PH2O = Water vapour partial pressure (Pa)
CC
CI JGM = Switch for GM:
CC
CC  0 will use value supplied as input from ccl.
CC  will use one of 7 values defined in Annex B/Table 2 of ISO 7730:
CC  1 Reclining (46 W/m^2 = 0.8 met)
CC  2 Seated, relaxed (58 W/m^2 = 1.0 met)
CC  3 Standing, relaxed (70 W/m^2 = 1.2 met)
CC  4 Sedentary activity (office, dwelling, school, laboratory, as in 3)
CC  5 Standing activity (shopping, laboratory, light industry; 93 W/m^2
CC    = 1.6 met)
CC  6 Standing activity (shop assistant, domestic work, machine work;
CC    116 W/m^2 = 2.0 met)
CC  7 Medium activity (heavy machine work, garage work; 165 W/m^2
CC     = 2.8 met).
CC
CI JGICL = Switch for GICL:
CC  0 will use value supplied as input from ccl.
CC  1 Nude (0 m^2C/W = 0 clo).
CC  2 Shorts (0.015 m^2C/W = 0.1 clo).
CC  3 Typical tropical clothing ensemble (briefs, shorts,open-neck
CC    short-sleeved shirt, light socks + sandals;
CC    0.045 m^2C/W = 0.3 clo).
CC  4 Light summer clothing (briefs, long light-weight trousers,
CC    open-neck short-sleeved shirt, light socks + shoes;
CC    0.08 m^2C/W = 0.5 clo).
CC  5 Light working ensemble (light underwear, cotton work shirt
CC    with long sleeves, work trousers, woollen socks + shoes;
CC    0.11 m^2C/W = 0.7 clo).
CC  6 Typical indoor winter clothing ensemble (underwear, shirt with
CC    long sleeves, trousers, jacket or sweater with long sleeves,
CC    heavy socks + shoes; 0.16 m^2C/W = 1.0 clo).
CC  7 Heavy traditional European business suit (cotton underwear
CC    with long legs and sleeves, shirt, suit including trousers,
CC    jacket and waiscoat, woollen socks + heavy shoes; 0.23 m^2C/W
CC    = 1.5 clo).
CC
CC --------------------
CC      Modified
CC --------------------
CC
CC --------------------
CC        Output
CC --------------------
CC
CO  PMV  = Predicted Mean Vote
CO  IERR = Error flag (0 -> o.k.)
CC
CC --------------------
CC       Details
CC --------------------
CC
CC Based on ISO 7730-1984(E)
CC
CC --------------------
CC Local variables & parameters
CC --------------------
CC
CC
CC --------------------
CC Modification History
CC --------------------
CC

CC May 2002 YLS  v1.
CC Jun 2002 CJS  Modified for CFX-5, v2.
CC Dec 2002 YLS  Added ISO standard GM & GICL tables, with switches
CC               from ccl, v3. Also TRA from the radiation solution
CC               rather than prescribed (see pmvppd.F).
CC Nov 2006 YLS  Bug corrected (definition of TA). v4.
CC Nov 2006 CJS  Tidy up. Bring up to CFX source code standard.
CC Jan 2008 NSW  Concatenated with USR_PMVPPD.
CC
C=======================================================================
C
C ---------------------------
C    Preprocessor includes
C ---------------------------
C
C ------------------------------
C        Argument list
C ------------------------------
C
      INTEGER            IERR
C
      REAL               GM, GW, GICL, TAA, TRA, VAIR, PH2O, PMV
C
C ------------------------------
C        External functions
C ------------------------------
C
C ------------------------------
C        Local Parameters
C ------------------------------
C
C ---------------------------
C        Local Variables
C ---------------------------
C
      INTEGER            NOI, JJ, JGM, JGICL
C
      REAL               EPS, GMW, FCL, FCIC, P1, P2, P3, P4,
     &                   TA, TR, TCLA, XN, XF, HCF, HCN, HC, TCL,
     &                   PM1, PM2, PM3, PM4, GMM(7), GICLL(7)
C
C---- Data statements
C
      DATA GMM(1), GMM(2), GMM(3), GMM(4), GMM(5), GMM(6)
     &    / 46.0,   58.0,   70.0,   70.0,   93.0,   116.0 /
      DATA GMM(7) / 165.0 /
      DATA GICLL(1), GICLL(2), GICLL(3), GICLL(4), GICLL(5), GICLL(6)
     &    /  0.0,      0.015,    0.045,    0.08,     0.11,     0.16 /
      DATA GICLL(7) / 0.23 /
C
C ---------------------------
C        Stack pointers
C ---------------------------
C
C ------------------------------
C  Directory on entry:- N/A
C ------------------------------
C
C---------------------------
C Executable Statements
C---------------------------
C
      IERR = 0
C
C---- Modify GM is JGM is not zero.
C
      IF(JGM.NE.0) THEN
       GM = GMM(JGM)
      ENDIF
C
C     Modify GICL is JGICL is not zero.
C
      IF(JGICL.NE.0) THEN
       GICL = GICLL(JGICL)
      ENDIF
C
      EPS=0.00015
      GMW = GM-GW
C
C---- Compute the corresponding FCL value
C
      FCL = 1.05+0.645*GICL
      IF(GICL.LT.0.078) FCL=FCL-0.05+0.645*GICL
      FCIC = GICL*FCL
      P2 = 3.96*FCIC
      P3 = 100*FCIC
      TA = TAA - 273.15
      TR = TRA - 273.15
      P1 = FCIC*TAA
      P4 = 308.7-0.028*GMW+P2*(TRA/100)**4
C
C---- First guess for surface temperature
C
      TCLA = TAA+(35.5-TA)/(3.5*(6.45*GICL+0.1))
      XN = TCLA/100
      XF = XN
      HCF = 12.1*SQRT(MAX(1.0E-20,ABS(VAIR)))
      NOI = 0
C
C---- Iterate
C
100   XF = 0.5*(XF+XN)
      HCN = 2.38*ABS(100*XF-TAA)**0.25
      HC = AMAX1(HCF,HCN)
      XN = (P4+P1*HC-P2*XF**4)/(100+P3*HC)
      NOI =NOI+1
      IF(NOI.GT.150) THEN
        IERR = 1
        GOTO 999
      ENDIF
      IF(ABS(XN-XF).GT.EPS) GO TO 100
      TCL = 100*XN-273.15
C
C---- Predicted mean vote
C
      PM1 = 3.96*FCL*(XN**4-(TRA/100)**4)
      PM2 = FCL*HC*(TCL-TA)
      PM3 = 0.303*EXP(-0.036*GM)+0.028
      PM4 = 0
      IF(GMW.GT.58.15) PM4=0.42*(GMW-58.15)
      PMV = PM3*(GMW-3.05*0.001*(5733-6.99*GMW-PH2O)
     &      -PM4-1.7*0.00001*GM*(5867-PH2O)
     &      -0.0014*GM*(34-TA)-PM1-PM2)
C
C---- Impose clipping
C
      PMV = SIGN(MIN(ABS(PMV), 3.0), PMV)
C
C---- Skip point
C
999   CONTINUE
C
      END
