text:
c     ************************************************************
c     * This package provides input routines to read the wind    *
c     * and other fields from IVE necdf files. The routines are  *
c     *                                                          *
c     * 1) input_open  : to open a data file                     *
c     * 2) input_grid  : to read the grid information, including *
c     *                  the vertical levels                     *
c     * 3) input_wind  : to read the wind components             *
c     * 4) input_close : to close an input file                  *
c     *                                                          *
c     * The file is characterised by an filename <filename> and  *
c     * a file identifier <fid>. The horizontal grid is given by *
c     * <xmin,xmax,ymin,ymax,dx,dy,nx,ny> where the pole of the  *
c     * rotated grid is given by <pollon,pollat>. The vertical   *
c     * grid is characterised by the surface pressure <ps> and   *
c     * the pressure at staggered <slev> and unstaggered <ulev>  *
c     * levels. The number of levels is given by <nz>. Finally,  *
c     * the retrieval of the wind <field> with name <fieldname>  *
c     * is characterised by a <time> and a missing data value    *
c     * <mdv>.                                                   *
c     *                                                          *
c     * Author: Michael Sprenger, Autumn 2008                    *
c     ************************************************************

c     ------------------------------------------------------------
c     Open input file
c     ------------------------------------------------------------

      subroutine input_open (fid,filename)

c     Open the input file with filename <filename> and return the
c     file identifier <fid> for further reference.

      implicit none

c     Declaration of subroutine parameters
      integer      fid              ! File identifier
      character*80 filename         ! Filename

c     Declaration of auxiliary variables
      integer      ierr

c     Open IVE netcdf file
      call cdfopn(filename,fid,ierr)
      if (ierr.ne.0) goto 900

c     Exception handling
      return

 900  print*,'Cannot open input file ',trim(filename)
      stop

      end

c     ------------------------------------------------------------
c     Read information about the grid
c     ------------------------------------------------------------

      subroutine input_grid
     >                   (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >                    time,pollon,pollat,p3,ps,nz,ak,bk,stagz,
     >                    timecheck)

c     Read grid information at <time> from file with identifier <fid>.
c     The horizontal grid is characterized by <xmin,xmax,ymin,ymax,dx,dy>
c     with pole position at <pollon,pollat> and grid dimension <nx,ny>.
c     The 3d arrays <p3(nx,ny,nz)> gives the vertical coordinates, either
c     on the staggered or unstaggered grid (with <stagz> as the flag).
c     The surface pressure is given in <ps(nx,ny)>. If <fid> is negative,
c     only the grid dimensions and grid parameters (xmin...pollat,nz) are
c     determined and returned (this is needed for dynamical allocation of
c     memory).

      implicit none

c     Declaration of subroutine parameters
      integer      fid                 ! File identifier
      real         xmin,xmax,ymin,ymax ! Domain size
      real         dx,dy               ! Horizontal resolution
      integer      nx,ny,nz            ! Grid dimensions
      real         pollon,pollat       ! Longitude and latitude of pole
      real         p3(nx,ny,nz)        ! Staggered levels
      real         ps(nx,ny)           ! Surface pressure
      real         time                ! Time of the grid information
      real         ak(nz),bk(nz)       ! Ak and Bk for layers or levels
      real         stagz               ! Vertical staggering (0 or -0.5)
      character*80 fieldname           ! Variable from which to take grid info
      character*80 timecheck           ! Either 'yes' or 'no'

c     Numerical and physical parameters
      real          eps                 ! Numerical epsilon
      parameter    (eps=0.001)

c     Netcdf variables
      integer      vardim(4)
      real         varmin(4),varmax(4)
      real         mdv
      real         stag(4)
      integer      ndim
      character*80 cstfile
      integer      cstid
      real         times(10)
      integer      ntimes
      real         aklay(nz),bklay(nz),aklev(nz),bklev(nz)
      integer      nvars
      character*80 vars(100)

c     Auxiliary varaibles
      integer      ierr
      integer      i,j,k
      integer      isok
      real         tmp(200)
      character*80 varname
      real         rtime
      integer      is2d
      integer      plev

c     Init the flag for 2D variables - assume a 3D field
      is2d = 0

c     Init the flag for pressure levels (PS is not needed)
      plev = 0

c     Inquire dimensions and grid constants if <fid> is negative
      if (fid.lt.0) then

c        Get grid info for the selected variable
         if ( fieldname.eq.'PLEV' ) then
            varname = 'PS'
            stagz   = 0.
            call getdef(-fid,varname,ndim,mdv,vardim,
     >                  varmin,varmax,stag,ierr)
            if (ierr.ne.0) goto 900
            call getcfn(-fid,cstfile,ierr)
            if (ierr.ne.0) goto 903
            call cdfopn(cstfile,cstid,ierr)
            if (ierr.ne.0) goto 903
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
            if (ierr.ne.0) goto 903
            call clscdf(cstid,ierr)
            if (ierr.ne.0) goto 903

         elseif ( ( fieldname.eq.'PLAY' ).or.( fieldname.eq.'P' ) ) then
            varname = 'PS'
            stagz   = -0.5
            call getdef(-fid,varname,ndim,mdv,vardim,
     >                  varmin,varmax,stag,ierr)
            if (ierr.ne.0) goto 900
            call getcfn(-fid,cstfile,ierr)
            if (ierr.ne.0) goto 903
            call cdfopn(cstfile,cstid,ierr)
            if (ierr.ne.0) goto 903
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
            if (ierr.ne.0) goto 903
            call clscdf(cstid,ierr)
            if (ierr.ne.0) goto 903

         else
            varname = fieldname
            call getdef(-fid,varname,ndim,mdv,vardim,
     >                  varmin,varmax,stag,ierr)
            if (ierr.ne.0) goto 900

         endif

c        Set the grid dimensions and constants - vardim(3) is taken from constants file
         nx   = vardim(1)
         ny   = vardim(2)
         nz   = vardim(3)
         xmin = varmin(1)
         ymin = varmin(2)
         xmax = varmax(1)
         ymax = varmax(2)
         dx   = (xmax-xmin)/real(nx-1)
         dy   = (ymax-ymin)/real(ny-1)

c        Get pole position - if no constants file available, set default pole
         call getcfn(-fid,cstfile,ierr)
         if (ierr.eq.0) then
            call cdfopn(cstfile,cstid,ierr)
            if (ierr.ne.0) goto 903
            call getpole(cstid,pollon,pollat,ierr)
            if (ierr.ne.0) goto 903
            call clscdf(cstid,ierr)
            if (ierr.ne.0) goto 903
         else
            pollon = 0.
            pollat = 90.
         endif

c     Get non-constant grid parameters (surface pressure and vertical grid)
      else

c        Special handling for fieldname 'P.ML' -> in this case the fields
c        P and PS are available on the data file and can be read in. There
c        is no need to reconstruct it from PS,AK and BK. This mode is
c        used for model-level (2D) trajectories
         if ( fieldname.eq.'P.ML' ) then

c           Get the right time to read
            call gettimes(fid,times,ntimes,ierr)
            if (ierr.ne.0) goto 901
            isok=0
            do i=1,ntimes
               if (abs(time-times(i)).lt.eps) then
                  isok = 1
                  rtime = times(i)
               elseif (timecheck.eq.'no') then
                  isok = 1
                  rtime = times(1)
               endif
            enddo

c           Read surface pressure and 3D pressure
            varname='PS'
            call getdat(fid,varname,rtime,0,ps,ierr)
            if (ierr.ne.0) goto 904
            varname='P'
            call getdat(fid,varname,rtime,0,p3,ierr)
            if (ierr.ne.0) goto 904

c           Set MDV to 1050. - otherwise interpolation routines don't work
            do i=1,nx
              do j=1,ny
                do k=1,nz
                   if ( p3(i,j,k).lt.0. ) p3(i,j,k) = 1050.
                enddo
              enddo
            enddo

c           Don't care about other stuff - finish subroutine
            goto 800

         endif

c        Get full grid info - in particular staggering flag; set flag for 2D tracing
         if ( fieldname.eq.'PLEV' ) then
            varname = 'PS'
            stagz   = 0.
            call getdef(fid,varname,ndim,mdv,vardim,
     >                  varmin,varmax,stag,ierr)
            if (ierr.ne.0) goto 900
            call getcfn(fid,cstfile,ierr)
            if (ierr.ne.0) goto 903
            call cdfopn(cstfile,cstid,ierr)
            if (ierr.ne.0) goto 903
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
            if (ierr.ne.0) goto 903
            call clscdf(cstid,ierr)
            if (ierr.ne.0) goto 903

         elseif ( ( fieldname.eq.'PLAY' ).or.( fieldname.eq.'P' ) ) then
            varname = 'PS'
            stagz   = -0.5
            call getdef(fid,varname,ndim,mdv,vardim,
     >                  varmin,varmax,stag,ierr)
            if (ierr.ne.0) goto 900
            call getcfn(fid,cstfile,ierr)
            if (ierr.ne.0) goto 903
            call cdfopn(cstfile,cstid,ierr)
            if (ierr.ne.0) goto 903
            call getlevs(cstid,vardim(3),tmp,tmp,tmp,tmp,ierr)
            if (ierr.ne.0) goto 903
            call clscdf(cstid,ierr)
            if (ierr.ne.0) goto 903

         else
            varname=fieldname
            call getdef(fid,varname,ndim,mdv,vardim,
     >                  varmin,varmax,stag,ierr)
            if (ierr.ne.0) goto 900
            if (vardim(3).eq.1) is2d = 1
         endif

c        Get time information (check if time is correct)
         call gettimes(fid,times,ntimes,ierr)
         if (ierr.ne.0) goto 901
         isok=0
         do i=1,ntimes
            if (abs(time-times(i)).lt.eps) then
               isok = 1
               rtime = times(i)
            elseif (timecheck.eq.'no') then
               isok = 1
               rtime = times(1)
            endif
         enddo
         if ( isok.eq.0) goto 905

c        If 2D tracing requested: take dummay values for PS, AKLEV,BKLEV,AKLAY,BKLAY
         if ( is2d.eq.1 ) then

            do i=1,nx
               do j=1,ny
                  ps(i,j) = 1050.
               enddo
            enddo

            do k=1,nz
               aklev(k) = 0.
               bklev(k) = real(nz-k)/real(nz-1)
               aklay(k) = 0.
               bklay(k) = real(nz-k)/real(nz-1)
            enddo

c        3D tracing - read PS, AKLEV,BKLEV,AKLAY;BKLAY
         else

c           Read the level coefficients from the constants file
            call getcfn(fid,cstfile,ierr)
            if (ierr.ne.0) goto 903
            call cdfopn(cstfile,cstid,ierr)
            if (ierr.ne.0) goto 903
            call getlevs(cstid,vardim(3),aklev,bklev,aklay,bklay,ierr)
            if (ierr.ne.0) goto 903
            call clscdf(cstid,ierr)
            if (ierr.ne.0) goto 903

c           Check whether PS is needed to get the 3d pressure field
            plev = 1
            do i=1,nz
              if ( (abs(stagz).lt.eps).and.(abs(bklev(i)).gt.eps) ) then
                plev = 0
              endif
              if ( (abs(stagz).gt.eps).and.(abs(bklay(i)).gt.eps) ) then
                plev = 0
              endif
            enddo

c           Read surface pressure if needed
            if ( plev.ne.1 ) then
              varname='PS'
              call getdat(fid,varname,rtime,0,ps,ierr)
              if (ierr.ne.0) goto 904
            else
              do i=1,nx
                do j=1,ny
                    ps(i,j) = 1000.
                enddo
              enddo
            endif

         endif

c        Calculate layer and level pressures
         do i=1,nx
            do j=1,ny
               do k=1,nz
                  if ( abs(stagz).lt.eps ) then
                     p3(i,j,k)=aklev(k)+bklev(k)*ps(i,j)
                  else
                     p3(i,j,k)=aklay(k)+bklay(k)*ps(i,j)
                  endif
               enddo
            enddo
         enddo

c        Set the ak and bk for the vertical grid
         do k=1,nz
            if ( abs(stagz).lt.eps ) then
               ak(k)=aklev(k)
               bk(k)=bklev(k)
            else
               ak(k)=aklay(k)
               bk(k)=bklay(k)
            endif
         enddo

      endif

c     Exit point for subroutine
 800  continue
      return

c     Exception handling
 900  print*,'Cannot retrieve grid dimension from ',fid
      stop
 901  print*,'Cannot retrieve grid parameters from ',fid
      stop
 902  print*,'Grid inconsistency detected for ',fid
      stop
 903  print*,'Problem with level coefficients from ',trim(cstfile)
      stop
 904  print*,'Cannot read surface pressure from ',trim(cstfile)
      stop
 905  print*,'Cannot find time ',time,' on ',fid
      stop
 906  print*,'Unable to get grid info ',fid
      stop

      end

c     ------------------------------------------------------------
c     Read wind information
c     ------------------------------------------------------------

      subroutine input_wind (fid,fieldname,field,time,stagz,mdv,
     >                       xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,
     >                       timecheck)

c     Read the wind component <fieldname> from the file with identifier
c     <fid> and save it in the 3d array <field>. The vertical staggering
c     information is provided in <stagz> and gives the reference to either
c     the layer or level field from <input_grid>. A consistency check is
c     performed to have an agreement with the grid specified by <xmin,xmax,
c     ymin,ymax,dx,dy,nx,ny,nz>.

      implicit none

c     Declaration of variables and parameters
      integer      fid                 ! File identifier
      character*80 fieldname           ! Name of the wind field
      integer      nx,ny,nz            ! Dimension of fields
      real         field(nx,ny,nz)     ! 3d wind field
      real         stagz               ! Staggering in the z direction
      real         mdv                 ! Missing data flag
      real         xmin,xmax,ymin,ymax ! Domain size
      real         dx,dy               ! Horizontal resolution
      real         time                ! Time
      character*80 timecheck           ! Either 'yes' or 'no'

c     Numerical and physical parameters
      real        eps                 ! Numerical epsilon
      parameter  (eps=0.001)
      real        notimecheck         ! 'Flag' for no time check
      parameter  (notimecheck=7.26537)

c     Netcdf variables
      integer      ierr
      character*80 varname
      integer      vardim(4)
      real         varmin(4),varmax(4)
      real         stag(4)
      integer      ndim
      real         times(10)
      integer      ntimes
      character*80 cstfile
      integer      cstid
      real         aklay(200),bklay(200),aklev(200),bklev(200)
      real         ps(nx,ny)

c     Auxiliary variables
      integer      isok
      integer      i,j,k
      integer      nz1
      real         rtime

c     Read variable definition - for P, PLEV and PLAY: load also ak,bk
      if ( ( fieldname.eq.'PLEV' ).or.
     >     ( fieldname.eq.'PLAY' ).or.
     >     ( fieldname.eq.'P'    ) )
     >then
         call getcfn(fid,cstfile,ierr)
         if (ierr.ne.0) goto 905
         call cdfopn(cstfile,cstid,ierr)
         if (ierr.ne.0) goto 905
         call getlevs(cstid,nz1,aklev,bklev,aklay,bklay,ierr)
         if (ierr.ne.0) goto 905
         call clscdf(cstid,ierr)
         if (ierr.ne.0) goto 905
         varname = 'PS'
         call getdef(fid,varname,ndim,mdv,vardim,
     >               varmin,varmax,stag,ierr)
         vardim(3) = nz1
         if (ierr.ne.0) goto 900

      else
         varname = fieldname
         call getdef(fid,varname,ndim,mdv,vardim,
     >               varmin,varmax,stag,ierr)
         if (ierr.ne.0) goto 900
         stagz=stag(3)
      endif

c     Get time information (set time to first one in the file)
      call gettimes(fid,times,ntimes,ierr)
      if (ierr.ne.0) goto 902
      isok=0
      do i=1,ntimes
         if (abs(time-times(i)).lt.eps) then
            isok = 1
            rtime = times(i)
         elseif (timecheck.eq.'no') then
            isok = 1
            rtime = times(1)
         endif
      enddo
      if ( isok.eq.0 ) goto 904

c     Read  field
      if ( ( fieldname.eq.'P' ).or.(fieldname.eq.'PLAY') )  then       ! P, PLAY
         stagz   = -0.5
         varname = 'PS'
         call getdat(fid,varname,rtime,0,ps,ierr)
         if (ierr.ne.0) goto 903
         do i=1,nx
            do j=1,ny
               do k=1,nz
                  field(i,j,k)=aklay(k)+bklay(k)*ps(i,j)
               enddo
            enddo
         enddo

      elseif ( fieldname.eq.'PLEV' )  then                             ! PLEV
         stagz   = 0.
         varname = 'PS'
         call getdat(fid,varname,rtime,0,ps,ierr)
         if (ierr.ne.0) goto 903
         do i=1,nx
            do j=1,ny
               do k=1,nz
                  field(i,j,k)=aklev(k)+bklev(k)*ps(i,j)
               enddo
            enddo
         enddo

      else                                                             ! Other fields
         varname=fieldname
         call getdat(fid,varname,rtime,0,field,ierr)
         if (ierr.ne.0) goto 903

      endif

c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
      if ( vardim(3).eq.1 ) then
         do i=1,nx
            do j=1,ny
               do k=1,nz
                  field(i,j,k) = field(i,j,1)
               enddo
            enddo
         enddo
      endif

c     Exception handling
      return

 900  print*,'Cannot retrieve definition for ',trim(varname),'  ',fid
      stop
 901  print*,'Grid inconsistency detected for ',trim(varname),'  ',fid
      stop
 902  print*,'Cannot retrieve time for ',trim(varname),'  ',fid
      stop
 903  print*,'Cannot load wind component ',trim(varname),'  ',fid
      stop
 904  print*,'Cannot load time ',time,' for ',trim(varname),'  ',fid
      stop
 905  print*,'Cannot load time vertical grid AK, BK from file  ',fid
      stop

      end

c     ------------------------------------------------------------
c     Close input file
c     ------------------------------------------------------------

      subroutine input_close(fid)

c     Close the input file with file identifier <fid>.

      implicit none

c     Declaration of subroutine parameters
      integer fid

c     Auxiliary variables
      integer ierr

c     Close file
      call clscdf(fid,ierr)

      end

c     ------------------------------------------------------------
c     Get a list of variables on netCDF file
c     ------------------------------------------------------------

      subroutine input_getvars(fid,vnam,nvars)

c     List of variables on netCDF file

      implicit none

c     Declaration of subroutine parameters
      integer      fid
      integer      nvars
      character*80 vnam(200)

c     Auxiliary variables
      integer ierr

c     Get list and save
      call getvars(fid,nvars,vnam,ierr)

      end

max_stars_repo_name: GEOS-ESM/AeroApps

max_issues_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


text:
C Copyright(C) 2009 Sandia Corporation. Under the terms of Contract
C DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
C certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C     * Neither the name of Sandia Corporation nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C=======================================================================
      SUBROUTINE DBOQA (NDB, NQAREC, QAREC, NINFO, INFO)
C=======================================================================
C$Id: dboqa.f,v 1.2 2009/03/25 12:46:01 gdsjaar Exp $
C$Log: dboqa.f,v $
CRevision 1.2  2009/03/25 12:46:01  gdsjaar
CAdd copyright and license notice to all files.
C
CRevision 1.1.1.1  1990/08/14 16:13:37  gdsjaar
CTesting
C
c Revision 1.1  90/08/14  16:13:36  gdsjaar
c Initial revision
c
c Revision 1.1  90/08/09  13:39:17  gdsjaar
c Initial revision
c

C   --*** DBOQA *** (EXOLIB) Write QA and information records
C   --   Written by Amy Gilkey - revised 02/08/88
C   --
C   --DBOQA writes the QA records and the information records.
C   --
C   --Parameters:
C   --   NDB - IN - the database number
C   --   NQAREC - IN - the number of QA records (must be at least one);
C   --      written only if >= 0
C   --   QAREC - IN - the QA records containing:
C   --      (1) - the analysis code name
C   --      (2) - the analysis code QA descriptor
C   --      (3) - the analysis date
C   --      (4) - the analysis time
C   --   NINFO - IN - the number of information records; written only if >= 0
C   --   INFO - IN - the information records
C   --
C   --Database must be positioned at start of QA records upon entry;
C   --upon exit positioned at end of information records.

      INTEGER NDB
      INTEGER NQAREC
      CHARACTER*8 QAREC(4,*)
      INTEGER NINFO
      CHARACTER*80 INFO(*)

      IF (NQAREC .LT. 0) GOTO 120

      WRITE (NDB) NQAREC

      DO 100 IQA = 1, NQAREC
         WRITE (NDB) (QAREC(I,IQA), I=1,4)
  100 CONTINUE

      IF (NINFO .LT. 0) GOTO 120

      WRITE (NDB) NINFO

      DO 110 I = 1, NINFO
         WRITE (NDB) INFO(I)
  110 CONTINUE

  120 CONTINUE
      RETURN
      END

ext: f

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas


text:
! Function for matrix product C=A*B
! Optimized memory access pattern is considered.
! In Fortran matrices should be accessed column by column.
! Igor Lopes, February 2015
SUBROUTINE MATPRD(A,B,C,NROWA,NCOLA,NCOLB)
    IMPLICIT NONE
    ! Parameter
    REAL(8) R0 /0.0D0/
    ! Arguments
    INTEGER :: NROWA , NCOLA , NCOLB
    REAL(8), DIMENSION(NROWA,NCOLA)  :: A
    REAL(8), DIMENSION(NCOLA, NCOLB) :: B
    REAL(8), DIMENSION(NROWA,NCOLB)  :: C
    ! Locals
    INTEGER I, J, K
    ! Begin algorithm
    C=R0
    DO J=1,NCOLB
	    DO K=1,NCOLA
		    DO I=1,NROWA
			    C(I,J)=C(I,J)+A(I,K)*B(K,J)
		    ENDDO
	    ENDDO
    ENDDO
END SUBROUTINE

max_stars_repo_name: iarlopes/GPROPT

max_issues_repo_name: iarlopes/GPROPT

max_forks_repo_name: iarlopes/GPROPT


text:
C Copyright(C) 2011-2017 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C * Redistributions of source code must retain the above copyright
C    notice, this list of conditions and the following disclaimer.
C
C * Redistributions in binary form must reproduce the above
C   copyright notice, this list of conditions and the following
C   disclaimer in the documentation and/or other materials provided
C   with the distribution.
C
C * Neither the name of NTESS nor the names of its
C   contributors may be used to endorse or promote products derived
C   from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C=======================================================================
      SUBROUTINE NEWMAP (MAPEL, MAPEL3, IXEL, INCEL, NREL, IELCOL)
C=======================================================================

C   --*** NEWMAP *** (GEN3D) Calculate 3D element order map
C   --   Written by Amy Gilkey - revised 05/05/86
C   --
C   --NEWMAP calculates the element order map for the 3D database.
C   --
C   --Parameters:
C   --   MAPEL - IN - the 2D element order map
C   --   MAPEL3 - OUT - the 3D element order map
C   --   IXEL - IN - the new index for each element
C   --   INCEL - IN - the increment for each element, needed for blocks
C   --      that become multiple blocks
C   --   NREL - IN - the number of new elements generated for each element
C   --   IELCOL - IN - the row number for each element, 0 if not needed
C   --
C   --Common Variables:
C   --   Uses NUMEL of /DBNUMS/
C   --   Uses NUMNP3 of /DBNUM3/
C   --   Uses NEREPL of /PARAMS/

      INCLUDE 'g3_dbnums.blk'
      INCLUDE 'g3_dbnum3.blk'
      INCLUDE 'g3_params.blk'

      INTEGER MAPEL(NUMEL), MAPEL3(NUMEL3)
      INTEGER IXEL(*), INCEL(*), NREL(*), IELCOL(*)

C   --Element order map - add on elements for each plate/slice
C ... To avoid potential numbering conflicts, we throw away the old map and just create a
C     dummy map from 1..numel3

      DO 20 I = 1, NUMEL3
        MAPEL3(I) = I
   20 CONTINUE

      RETURN
      END

max_stars_repo_name: brian-kelley/seacas

max_issues_repo_name: brian-kelley/seacas

max_forks_repo_name: brian-kelley/seacas


text:
      SUBROUTINE CALIMP2T
C
C **  SUBROUTINE CALEXP CALCULATES IMPLICIT MOMENTUM EQUATION
C **  CORIOLIS AND CURVATURE TERMS FOR 1/2 STEP PREDICTOR
C CHANGE RECORD
C
      USE GLOBAL
      IF(ISDYNSTP.EQ.0)THEN
        DELT=DT
        DELTD2=0.5*DT
        DELTI=1./DELT
      ELSE
        DELT=DTDYN
        DELTD2=0.5*DTDYN
        DELTI=1./DELT
      ENDIF
      DO K=1,KC
        TVAR1W(1,K)=0.
        TVAR1S(1,K)=0.
        TVAR1E(1,K)=0.
        TVAR1N(1,K)=0.
        TVAR2W(1,K)=0.
        TVAR2S(1,K)=0.
        TVAR2E(1,K)=0.
        TVAR2N(1,K)=0.
        TVAR1W(LC,K)=0.
        TVAR1S(LC,K)=0.
        TVAR1E(LC,K)=0.
        TVAR1N(LC,K)=0.
        TVAR2W(LC,K)=0.
        TVAR2S(LC,K)=0.
        TVAR2E(LC,K)=0.
        TVAR2N(LC,K)=0.
      ENDDO
      DO K=1,KC
        DO L=2,LA
          TVAR1W(L,K)=0.
          TVAR1S(L,K)=0.
          TVAR1E(L,K)=0.
          TVAR1N(L,K)=0.
          TVAR2W(L,K)=0.
          TVAR2S(L,K)=0.
          TVAR2E(L,K)=0.
          TVAR2N(L,K)=0.
        ENDDO
      ENDDO
      TVAR3W(1)=0.
      TVAR3S(1)=0.
      TVAR3E(1)=0.
      TVAR3N(1)=0.
      TVAR3W(LC)=0.
      TVAR3S(LC)=0.
      TVAR3E(LC)=0.
      TVAR3N(LC)=0.
      DO ND=1,NDM
        LF=2+(ND-1)*LDM
        LL=LF+LDM-1
        DO L=LF,LL
          TVAR3W(L)=0.
          TVAR3S(L)=0.
          TVAR3E(L)=0.
          TVAR3N(L)=0.
        ENDDO
      ENDDO
C
C **  INITIALIZE EXTERNAL CORIOLIS-CURVATURE AND ADVECTIVE FLUX TERMS
C
      DO L=1,LC
        FCAXE(L)=0.
        FCAYE(L)=0.
C PMC          FCAX1E(L)=0.
C PMC          FCAY1E(L)=0.
        FXE(L)=0.
        FYE(L)=0.
      ENDDO
C
C **  CALCULATE TIME LEVEL TIME LEVEL N INTERNAL STRESSES
C **  AND WAVE INDUCED CURRENT FORCINGS
C **  CALCULATE FRACTIONAL STEP WITH STRESSES AND WV FORCES
C
      ! *** Surface & Bottom Stresses
      DO ND=1,NDM
        LF=2+(ND-1)*LDM
        LL=LF+LDM-1
        DO L=LF,LL
          TVAR2W(L,1)=TBX(L)
          TVAR2S(L,1)=TBY(L)
          TVAR2E(L,KC)=TSX(L)
          TVAR2N(L,KC)=TSY(L)
        ENDDO
      ENDDO
      ! *** Vertical Eddy Viscosity Impacts
      DO K=1,KS
        TMPVAL=DZIG(K)
        DO L=2,LA
          TVAR2E(L,K)=TMPVAL*(U(L,K+1)-U(L,K))/AVUI(L,K)
          TVAR2N(L,K)=TMPVAL*(V(L,K+1)-V(L,K))/AVVI(L,K)
        ENDDO
      ENDDO

      DO K=2,KC
        DO L=2,LA
          TVAR2W(L,K)=TVAR2E(L,K-1)
          TVAR2S(L,K)=TVAR2N(L,K-1)
        ENDDO
      ENDDO

      DO K=1,KC
        DO L=2,LA
          TVAR1W(L,K)=DYIU(L)*UHDY(L,K)
          TVAR1S(L,K)=DXIV(L)*VHDX(L,K)
        ENDDO
      ENDDO
      DO K=1,KC
        TMPVAL=DELTD2*DZIC(K)
        DO L=2,LA
          ! *** Vertical Shear in the X direction
          TVAR1W(L,K)=TVAR1W(L,K)+TMPVAL*(TVAR2E(L,K)-TVAR2W(L,K))
          ! *** Vertical Shear in the Y direction
          TVAR1S(L,K)=TVAR1S(L,K)+TMPVAL*(TVAR2N(L,K)-TVAR2S(L,K))
        ENDDO
      ENDDO
C
C *** DSLLC BEGIN BLOCK
C
      IF(ISWAVE.EQ.2)THEN
        IF(N.LT.NTSWV)THEN
          TMPVAL=FLOAT(N)/FLOAT(NTSWV)
          WVFACT=0.5-0.5*COS(PI*TMPVAL)
        ELSE
          WVFACT=1.0
        ENDIF
        WVFACT=DELTD2*WVFACT
C
        DO K=1,KC
          DO L=2,LA
            TVAR1W(L,K)=TVAR1W(L,K)+WVFACT*SAAX(L)*DXYIU(L)*FXWAVE(L,K)
            TVAR1S(L,K)=TVAR1S(L,K)+WVFACT*SAAY(L)*DXYIV(L)*FYWAVE(L,K)
          ENDDO
        ENDDO
      ENDIF
C
C *** DSLLC END BLOCK
C
      DO K=1,KC
        DO L=2,LA
          TVAR1W(L,K)=SUB(L)*TVAR1W(L,K)
          TVAR1S(L,K)=SVB(L)*TVAR1S(L,K)
        ENDDO
      ENDDO
C
C **  ADD PRESSURE GRADIENT TERMS
C **  HRU=SUB*HMU*DYU/DXU & HRV=SVB*HMV*DXV/DYV
C
      DO ND=1,NDM
        LF=2+(ND-1)*LDM
        LL=LF+LDM-1
        DO L=LF,LL
          TVAR1E(L,KC)=DZC(KC)*B(L,KC)
        ENDDO
      ENDDO
      DO K=KS,1,-1
        DO L=2,LA
          TVAR1E(L,K)=TVAR1E(L,K+1)+DZC(K)*B(L,K)
        ENDDO
      ENDDO
      DO K=1,KC
        DO L=2,LA
          TVAR1E(L,K)=TVAR1E(L,K)-0.5*DZC(K)*B(L,K)
        ENDDO
      ENDDO
      DO K=1,KC
        DO L=2,LA
          TVAR1N(L,K)=TVAR1E(L,K)+0.5*(Z(K)+Z(K-1))*B(L,K)
        ENDDO
      ENDDO
      DO K=1,KC
        DO L=2,LA
          TVAR1N(L,K)=GP*TVAR1N(L,K)
          TVAR1E(L,K)=GP*TVAR1E(L,K)
        ENDDO
      ENDDO
      DO ND=1,NDM
        LF=2+(ND-1)*LDM
        LL=LF+LDM-1
        DO L=LF,LL
          TVAR3W(L)=DELTD2*SUB(L)*HU(L)*DXIU(L)
          TVAR3S(L)=DELTD2*SVB(L)*HV(L)*DYIV(L)
        ENDDO
      ENDDO
      DO K=1,KC
        DO L=2,LA
          LS=LSC(L)
          TVAR1W(L,K)=TVAR1W(L,K)-TVAR3W(L)*( P(L)-P(LWEST(L)) )
          TVAR1S(L,K)=TVAR1S(L,K)-TVAR3S(L)*( P(L)-P(LS ) )
        ENDDO
      ENDDO
      DO ND=1,NDM
        LF=2+(ND-1)*LDM
        LL=LF+LDM-1
        DO L=LF,LL
          TVAR3W(L)=DELT*SUB(L)*HU(L)*DXYIU(L)*SBX(L)  ! *** PMC H1U --> HU
          TVAR3S(L)=DELT*SVB(L)*HV(L)*DXYIV(L)*SBY(L)  ! *** PMC H1V --> HV
        ENDDO
      ENDDO
      DO K=1,KC
        DO L=2,LA
          LS=LSC(L)
          TVAR1W(L,K)=TVAR1W(L,K)-TVAR3W(L)*(
     &        +HU(L)*(TVAR1E(L,K)-TVAR1E(LWEST(L),K))  ! *** PMC H1U --> HU
     &        +0.5*GP*(B(L,K)+B(LWEST(L),K))*(BELV(L)-BELV(LWEST(L)))
     &        +0.5*(TVAR1N(L,K)+TVAR1N(LWEST(L),K))*(HP(L)-HP(LWEST(L))) )
          TVAR1S(L,K)=TVAR1S(L,K)-TVAR3S(L)*(
     &        +HV(L)*(TVAR1E(L,K)-TVAR1E(LS ,K))   ! *** PMC H1V --> HV
     &        +0.5*GP*(B(L,K)+B(LS ,K))*(BELV(L)-BELV(LS ))
     &        +0.5*(TVAR1N(L,K)+TVAR1N(LS ,K))*(HP(L)-HP(LS )) )
        ENDDO
      ENDDO
C
C **  ADVECTION STEP FOR X MOMENTUM EQUATION
C **  SELECT ADVECTIVE FLUX FORM
C
      DO K=1,KC
        UHDY2(1,K)=0.
        VHDX2(1,K)=0.
        UHDY2(LC,K)=0.
        VHDX2(LC,K)=0.
        FUHU(1,K)=0.
        FUHU(LC,K)=0.
        FVHU(1,K)=0.
        FVHU(LC,K)=0.
      ENDDO
      DO K=0,KC
        W2(1,K)=0.
        FWU(1,K)=0.
        W2(LC,K)=0.
        FWU(LC,K)=0.
      ENDDO
      DO ND=1,NDM
        LF=2+(ND-1)*LDM
        LL=LF+LDM-1
        DO L=LF,LL
          W2(L,0)=0.
          FWU(L,0)=0.
        ENDDO
      ENDDO
      DO K=1,KC
        DO L=2,LC
          UHDY2(L,K)=0.
          VHDX2(L,K)=0.
          W2(L,K)=0.
          FUHU(L,K)=0.
          FVHU(L,K)=0.
          FWU(L,K)=0.
        ENDDO
      ENDDO
C
C **  ADVECTIVE FLUX HALF STEP
C **  CALCULATE ADVECTIVE FLUXES BY UPWIND DIFFERENCE WITH ADVECTION
C **  AT N AND ADVECTED FIELD AT N
C **  UHDY2(L,K) IS DYU*U ON POS X U(L) MOM CV FACE (AT B OR P POINT)
C **  VHDX2(L,K) IS DXV*V ON NEG Y U(L) MOM CV FACE (AT CORNER POINT)
C **  W2(L,K)    IS DXYP*W/H ON TOP U(L) MOM CV FACE
C ** SMOLARKIEZC AND MARGOLIN FORM
C ** STANDARD FORM
C
      DO K=1,KC
        DO L=2,LA
          LE=LEAST(L)
          LW=LWEST(L)
C PMC         UHDY2(L,K)=0.25*STCUV(L)*( (UHDY(L  ,K)+UHDY(L  ,K)
          UHDY2(L,K)=0.25*( (UHDY(L  ,K)+UHDY(L  ,K))
     &        +(UHDY(LE,K)+UHDY(LE,K)) )
          VHDX2(L,K)=0.25*( (VHDX(L  ,K)+VHDX(L  ,K))
     &        +(VHDX(LW,K)+VHDX(LW,K)) )
        ENDDO
      ENDDO
      ! *** DSLLC BEGIN BLOCK
      IF(ITRICELL.GT.0)THEN
        DO K=1,KC
          DO L=2,LA
            UHDY2(L,K)=UHDY2(L,K)*STCUV(L)
            VHDX2(L,K)=VHDX2(L,K)*STCUV(L)
          ENDDO
        ENDDO
      ENDIF
      ! *** DSLLC END BLOCK
      DO K=1,KS
        DO L=2,LA
          LW=LWEST(L)
          W2(L,K)=0.25*( DXYP(L )*(W(L ,K)+W(L ,K))
     &                  +DXYP(LW)*(W(LW,K)+W(LW,K)) )
        ENDDO
      ENDDO
C
C ** SMOLARKIEZC AND MARGOLIN FORM
C ** STANDARD FORM
C
      DO K=1,KC
        DO L=2,LA
          LS=LSC(L)
          LE=LEAST(L)
          FUHU(L,K)=MAX(UHDY2(L,K),0.)*U(L ,K)
     &             +MIN(UHDY2(L,K),0.)*U(LE,K)
          FVHU(L,K)=MAX(VHDX2(L,K),0.)*U(LS,K)
     &             +MIN(VHDX2(L,K),0.)*U(L ,K)
        ENDDO
      ENDDO
      DO K=1,KS
        DO L=2,LA
          FWU(L,K)=MAX(W2(L,K),0.)*U(L,K  )
     &            +MIN(W2(L,K),0.)*U(L,K+1)
        ENDDO
      ENDDO
C
C **  ADVANCE U ADVECTION 1/2 TIME STEP
C
      DO K=1,KC
        DO L=2,LA
          FX(L,K)=FUHU(L,K)-FUHU(LWEST(L),K)+FVHU(LNC(L),K)-FVHU(L,K)
        ENDDO
      ENDDO
      DO K=1,KC
        DO L=2,LA
          FX(L,K)=FX(L,K)+(FWU(L,K)-FWU(L,K-1))*DZIC(K)
        ENDDO
      ENDDO
      TMPVAL=DELTD2*SNLT
      DO K=1,KC
        DO L=2,LA
          TVAR1W(L,K)=TVAR1W(L,K)-TMPVAL*SUB(L)*DXYIU(L)*SAAX(L)*FX(L,K)
        ENDDO
      ENDDO
C
C **  ADVECTION STEP FOR Y MOMENTUM EQUATION
C **  SELECT ADVECTIVE FLUX FORM
C
      DO K=1,KC
        UHDY2(1,K)=0.
        VHDX2(1,K)=0.
        UHDY2(LC,K)=0.
        VHDX2(LC,K)=0.
        FUHV(1,K)=0.
        FUHV(LC,K)=0.
        FVHV(1,K)=0.
        FVHV(LC,K)=0.
      ENDDO
      DO K=0,KC
        W2(1,K)=0.
        FWV(1,K)=0.
        W2(LC,K)=0.
        FWV(LC,K)=0.
      ENDDO
      DO ND=1,NDM
        LF=2+(ND-1)*LDM
        LL=LF+LDM-1
        DO L=LF,LL
          W2(L,0)=0.
          FWV(L,0)=0.
        ENDDO
      ENDDO
      DO K=1,KC
        DO L=2,LC
          UHDY2(L,K)=0.
          VHDX2(L,K)=0.
          W2(L,K)=0.
          FUHV(L,K)=0.
          FVHV(L,K)=0.
          FWV(L,K)=0.
        ENDDO
      ENDDO
C
C **  TWO TIME LEVEL STEP
C **  CALCULATE ADVECTIVE FLUXES BY UPWIND DIFFERENCE WITH ADVECTION
C **  AVERAGED BETWEEN (N) AND (N+1) AND ADVECTED FIELD AT N
C **  UHDY2(L,K) IS DYU*U ON NEG X V(L) MOM CV FACE (CORNER POINT)
C **  VHDX2(L,K) IS DXV*V ON POS Y V(L) MOM CV FACE (B 0R P POINT)
C **  W2(L,K)    IS DXYP*W/H ON TOP V MOM CV FACE
C ** SMOLARKIEWCZ AND MARGOLIN FORM
C ** STANDARD FORM
C
      DO K=1,KC
        DO L=2,LA
          LS=LSC(L)
          LN=LNC(L)
          UHDY2(L,K)=0.25*( (UHDY(L ,K)+UHDY(L ,K))
     &        +(UHDY(LS,K)+UHDY(LS,K)) )
          VHDX2(L,K)=0.25*STCUV(L)*( (VHDX(L ,K)+VHDX(L ,K))
     &        +(VHDX(LN,K)+VHDX(LN,K)) )
        ENDDO
      ENDDO
      DO K=1,KS
        DO L=2,LA
          LS=LSC(L)
          W2(L,K)=0.25*( DXYP(L )*(W(L ,K)+W(L ,K))
     &        +DXYP(LS)*(W(LS,K)+W(LS,K)) )
        ENDDO
      ENDDO
C
C ** SMOLARKIEWCZ AND MARGOLIN FORM
C ** STANDARD FORM
C
      DO K=1,KC
        DO L=2,LA
          LW=LWEST(L)
          LN=LNC(L)
          FUHV(L,K)=MAX(UHDY2(L,K),0.)*V(LW,K)
     &             +MIN(UHDY2(L,K),0.)*V(L ,K)
          FVHV(L,K)=MAX(VHDX2(L,K),0.)*V(L ,K)
     &             +MIN(VHDX2(L,K),0.)*V(LN,K)
        ENDDO
      ENDDO
      DO K=1,KS
        DO L=2,LA
          FWV(L,K)=MAX(W2(L,K),0.)*V(L,K)
     &        +MIN(W2(L,K),0.)*V(L,K+1)
        ENDDO
      ENDDO
C
C **  ADVANCE V ADVECTION 1/2 TIME STEP
C
      DO K=1,KC
        DO L=2,LA
          LE=LEAST(L)
          LS=LSC(L)
          FY(L,K)=FUHV(LE,K)-FUHV(L,K)+FVHV(L,K)-FVHV(LS,K)
        ENDDO
      ENDDO
      DO K=1,KC
        DO L=2,LA
          FY(L,K)=FY(L,K)+(FWV(L,K)-FWV(L,K-1))*DZIC(K)
        ENDDO
      ENDDO
      TMPVAL=DELTD2*SNLT
      DO K=1,KC
        DO L=2,LA
          TVAR1S(L,K)=TVAR1S(L,K)-TMPVAL*SVB(L)*DXYIV(L)*SAAY(L)*FY(L,K)
        ENDDO
      ENDDO
C
C **  ADVANCE CONTINUITY EQUATION 1/2 STEP
C **  INITIALIZE MID TIME LEVEL VELOCITY
C **  BEGIN CORIOLIS AND CURVATURE ACCELERATION ITERATION
C **  CALCULATE CORIOLIS AND CURVATURE PARAMETER
C
      DO K=1,KC
        DO L=2,LA
          LN=LNC(L)
          LE=LEAST(L)
          CAC(L,K)=( FCORC(L)*DXYP(L)
     &        +0.5*SNLT*(V(LN ,K)+V(L,K))*DYDI(L)
     &        -0.5*SNLT*(U(LE,K)+U(L,K))*DXDJ(L) )*HP(L)
        ENDDO
      ENDDO
C
C **  CALCULATE CORIOLIS-CURVATURE AND ADVECTIVE ACCELERATIONS
C
      DO K=1,KC
        DO L=2,LA
          LN=LNC(L)
          LS=LSC(L)
          LNW=LNWC(L)
          LSE=LSEC(L)
          LE=LEAST(L)
          LW=LWEST(L)
          FCAX(L,K)=0.25*SCAX(L)*(CAC(L,K)*(V(LN,K)+V(L,K))
     &        +CAC(LW,K)*(V(LNW,K)+V(LW,K)))
          FCAY(L,K)=0.25*SCAY(L)*(CAC(L,K)*(U(LE,K)+U(L,K))
     &        +CAC(LS,K)*(U(LSE,K)+U(LS,K)))
        ENDDO
      ENDDO
C
C **  CALCULATE MID TIME LEVEL TRANSPORT AND VELOCTIY
C
      DO K=1,KC
        DO L=2,LA
          TVAR1W(L,K)=TVAR1W(L,K)+DELTD2*SUB(L)*DXYIU(L)*FCAX(L,K)
          TVAR1S(L,K)=TVAR1S(L,K)-DELTD2*SVB(L)*DXYIV(L)*FCAY(L,K)
        ENDDO
      ENDDO
C
C **  ADD HORIZONTAL MOMENTUN DIFFUSION TO ADVECTIVE ACCELERATIONS
C **  CALCULATE EXTERNAL ACCELERATIONS
C
      DO K=1,KC
        DO L=2,LA
          FXE(L)=FXE(L)+DYU(L)*TVAR1W(L,K)*DZC(K)
          FYE(L)=FYE(L)+DXV(L)*TVAR1S(L,K)*DZC(K)
        ENDDO
      ENDDO
C
C **  ADD REMAINING SURFACE AND BOTTOM STRESS
C
      DO L=2,LA
        FXE(L)=FXE(L)+DELTD2*SUB(L)*DYU(L)*(TSX(L)-TBX1(L))  ! *** PMC WHY TBX1 NOT TBX ?
        FYE(L)=FYE(L)+DELTD2*SVB(L)*DXV(L)*(TSY(L)-TBY1(L))
      ENDDO
C
C **  CALCULATE EXPLICIT INTERNAL U AND V SHEAR EQUATION TERMS
C **  CALCULATE EXPLICIT INTERNAL BUOYANCY FORCINGS CENTERED AT N FOR
C **  THREE TIME LEVEL STEP AND AT (N+1/2) FOR TWO TIME LEVEL STEP
C **  SBX=SBX*0.5*DYU & SBY=SBY*0.5*DXV
C
      DO K=1,KS
        DO L=2,LA
          LS=LSC(L)
          FBBX(L,K)=SBX(L)*GP*HU(L)*
     &        ( HU(L)*( (B(L,K+1)-B(LWEST(L),K+1))*DZC(K+1)
     &        +(B(L,K)-B(LWEST(L),K))*DZC(K) )
     &        -(B(L,K+1)-B(L,K)+B(LWEST(L),K+1)-B(LWEST(L),K))*
     &        (BELV(L)-BELV(LWEST(L))+Z(K)*(HP(L)-HP(LWEST(L)))) )
          FBBY(L,K)=SBY(L)*GP*HV(L)*
     &        ( HV(L)*( (B(L,K+1)-B(LS,K+1))*DZC(K+1)
     &        +(B(L,K)-B(LS,K))*DZC(K) )
     &        -(B(L,K+1)-B(L,K)+B(LS,K+1)-B(LS,K))*
     &        (BELV(L)-BELV(LS)+Z(K)*(HP(L)-HP(LS))) )
        ENDDO
      ENDDO
 1111 FORMAT(2I5,2X,8E12.4)
C
C **  CALCULATE EXPLICIT INTERNAL U AND V SHEAR EQUATION TERMS
C      ELSE
C
      DO K=1,KS
        RCDZF=CDZF(K)
        DO L=2,LA
          DU(L,K)=RCDZF*( 2.*(TVAR1W(L,K+1)-TVAR1W(L,K))*DELTI
     &        +DXYIU(L)*FBBX(L,K) )
          DV(L,K)=RCDZF*( 2.*(TVAR1S(L,K+1)-TVAR1S(L,K))*DELTI
     &        +DXYIV(L)*FBBY(L,K) )
        ENDDO
      ENDDO

      IF(NWSER.GT.0)THEN
        DO L=2,LA
          DU(L,KS)=DU(L,KS)-CDZU(KS)*TSX(L)
          DV(L,KS)=DV(L,KS)-CDZU(KS)*TSY(L)
        ENDDO
      ENDIF
C
      RETURN
      END


text:
!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------

!C
!C***
!C*** module hecmw_precond_ML_33
!C***
!C
module hecmw_precond_ML_33
  use  hecmw_util

  private

  public:: hecmw_precond_ML_33_setup
  public:: hecmw_precond_ML_33_apply
  public:: hecmw_precond_ML_33_clear

  integer(kind=kint), save :: id

  logical, save :: INITIALIZED = .false.

contains

  subroutine hecmw_precond_ML_33_setup(hecMAT, hecMESH, sym)
    use hecmw_matrix_misc
    use hecmw_mat_id
    implicit none
    type(hecmwST_matrix), intent(inout) :: hecMAT
    type(hecmwST_local_mesh), intent(in) :: hecMESH
    integer(kind=kint), intent(in) :: sym
    integer(kind=kint) :: ierr
    integer(kind=kint), save :: n_recycle = 0
    if (INITIALIZED) then
      if (hecMAT%Iarray(98) == 1) then      ! need symbolic and numerical setup
        call hecmw_precond_ML_33_clear()
      else if (hecMAT%Iarray(97) == 1) then ! need numerical setup only
        call hecmw_precond_ML_33_clear()
      else                                  ! no need to setup or skip setup
        call hecmw_mat_id_clear(id)
        call hecmw_mat_id_set(hecMAT, hecMESH, id)
        return
      endif
    endif
    call hecmw_mat_id_set(hecMAT, hecMESH, id)
    call hecmw_ML_wrapper_setup_33(id, sym, ierr)
    INITIALIZED = .true.
    hecMAT%Iarray(98) = 0 ! symbolic setup done
    hecMAT%Iarray(97) = 0 ! numerical setup done
    n_recycle = 0
  end subroutine hecmw_precond_ML_33_setup

  subroutine hecmw_precond_ML_33_apply(WW)
    implicit none
    real(kind=kreal), intent(inout) :: WW(:)
    integer(kind=kint) :: ierr
    call hecmw_ML_wrapper_apply_33(id, WW, ierr)
  end subroutine hecmw_precond_ML_33_apply

  subroutine hecmw_precond_ML_33_clear()
    use hecmw_mat_id
    implicit none
    integer(kind=kint) :: ierr
    call hecmw_ML_wrapper_clear_33(id, ierr)
    call hecmw_mat_id_clear(id)
    INITIALIZED = .false.
  end subroutine hecmw_precond_ML_33_clear

end module     hecmw_precond_ML_33

max_stars_repo_name: masae-hayashi/FrontISTR

max_issues_repo_name: masae-hayashi/FrontISTR

max_forks_repo_name: masae-hayashi/FrontISTR


text:
Davis Municipal Code Davis Municipal Code/1 Chapter 1
 GENERAL PROVISIONS
 1.01.050 Effect on past actions and obligations.

Neither the adoption of this Code nor the repeal hereby of any ordinance of this city shall in any manner affect the prosecution for violations of ordinances, which violations were committed prior to its adoption, nor be construed as a waiver of any license or penalty due and unpaid under such ordinances, nor construed as affecting any of the provisions of such ordinances relating to the collection of any such license or penalty or penal provisions applicable to any violation thereof, nor to affect the validity of any bond or cash deposit in lieu thereof, required to be posted, filed or deposited pursuant to any ordinance, and all rights and obligations thereunder appertaining shall continue in full force and effect.

The repeal of an ordinance shall not revive any ordinances in force before or at the time the ordinance repealed took effect. (Code 1964, § 11.03.)

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      module if_global_m
        implicit none
        real*8, dimension(:,:), pointer :: if_kappa
      end module if_global_m

max_stars_repo_name: SCOREC/phasta

max_issues_repo_name: SCOREC/phasta

max_forks_repo_name: SCOREC/phasta


text:
      subroutine gf_unpack1(cgrib,lcgrib,iofst,ids,idslen,ierr)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .                                       .
! SUBPROGRAM:    gf_unpack1
!   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2000-05-26
!
! ABSTRACT: This subroutine unpacks Section 1 (Identification Section)
!   starting at octet 6 of that Section.
!
! PROGRAM HISTORY LOG:
! 2000-05-26  Gilbert
! 2002-01-24  Gilbert  - Changed to dynamically allocate arrays
!                        and to pass pointers to those arrays through
!                        the argument list.
!
! USAGE:    CALL gf_unpack1(cgrib,lcgrib,iofst,ids,idslen,ierr)
!   INPUT ARGUMENT LIST:
!     cgrib    - Character array containing Section 1 of the GRIB2 message
!     lcgrib   - Length (in bytes) of GRIB message array cgrib.
!     iofst    - Bit offset of the beginning of Section 1.
!
!   OUTPUT ARGUMENT LIST:
!     iofst    - Bit offset at the end of Section 1, returned.
!     ids      - Pointer to integer array containing information read from
!                Section 1, the Identification section.
!            ids(1)  = Identification of originating Centre
!                                 ( see Common Code Table C-1 )
!            ids(2)  = Identification of originating Sub-centre
!            ids(3)  = GRIB Master Tables Version Number
!                                 ( see Code Table 1.0 )
!            ids(4)  = GRIB Local Tables Version Number
!                                 ( see Code Table 1.1 )
!            ids(5)  = Significance of Reference Time (Code Table 1.2)
!            ids(6)  = Year ( 4 digits )
!            ids(7)  = Month
!            ids(8)  = Day
!            ids(9)  = Hour
!            ids(10)  = Minute
!            ids(11)  = Second
!            ids(12)  = Production status of processed data
!                                 ( see Code Table 1.3 )
!            ids(13)  = Type of processed data ( see Code Table 1.4 )
!     idslen   - Number of elements in ids().
!     ierr     - Error return code.
!                0 = no error
!                6 = memory allocation error
!
! REMARKS:
!
! ATTRIBUTES:
!   LANGUAGE: Fortran 90
!   MACHINE:  IBM SP
!
!$$$

      character(len=1),intent(in) :: cgrib(lcgrib)
      integer,intent(in) :: lcgrib
      integer,intent(inout) :: iofst
      integer,pointer,dimension(:) :: ids
      integer,intent(out) :: ierr,idslen

      integer,dimension(:) :: mapid(13)

      data mapid /2,2,1,1,1,2,1,1,1,1,1,1,1/

      ierr=0
      idslen=13
      nullify(ids)

      call gbyte(cgrib,lensec,iofst,32)        ! Get Length of Section
      iofst=iofst+32
      iofst=iofst+8     ! skip section number
      !
      !   Unpack each value into array ids from the
      !   the appropriate number of octets, which are specified in
      !   corresponding entries in array mapid.
      !
      istat=0
      allocate(ids(idslen),stat=istat)
      if (istat.ne.0) then
         ierr=6
         nullify(ids)
         return
      endif

      do i=1,idslen
        nbits=mapid(i)*8
        call gbyte(cgrib,ids(i),iofst,nbits)
        iofst=iofst+nbits
      enddo

      return    ! End of Section 1 processing
      end

max_stars_repo_name: martinremy/wps

max_issues_repo_name: martinremy/wps

max_forks_repo_name: WRF-CMake/WPS


text:
! +-----------------------------------------------------------------------------+
! | TRNSYS Type3254: Variable capacity air-air heat pump with performance files |
! +-----------------------------------------------------------------------------+

! This routine implements an air-air heat pump with variable speed compressor.

! Inputs
! --------------------------------------------------------------------------------------------------
!  # | Variable     | Description                                   | Input  Units  | Internal Units
! --------------------------------------------------------------------------------------------------
!  1 | Tr           | Inlet (return) air temperature                | °C            | °C
!  2 | wr           | Inlet (return) air humidity ratio             | -             | -
!  3 | RHr          | Inlet (return) air relative humidity          | % (base 100)  | -
!  4 | pr           | Inlet (return) air pressure                   | atm           | atm
!  5 | mDot         | Inlet (return) air mass flow rate             | kg/h          | kg/h
!  6 | Toa          | Outdoor air dry bulb temperature              | °C            | °C
!  7 | woa          | Outdoor air humidity ratio                    | -             | -
!  8 | RHoa         | Outdoor air realtive humidity                 | % (base 100)  | -
!  9 | poa          | Outdoor air pressure                          | atm           | atm
! 10 | freq         | Normalized frequency signal                   | -             | -
! 11 | AFR          | Inlet (return) normalized air flow rate       | -             | -
! 12 | mode         | 0 = cooling mode                              | -             | -
!                   | 1 = heating mode                              |               |
! 13 | defrost_mode | 0 = defrost (off) mode                        | -             | -
!                   | 1 = recovery mode (transient)                 |               |
!                   | 2 = steady-state mode                         |               |
! 14 | recov_penalty| Penalty factor for defrost recovery mode      | -             | -
! --------------------------------------------------------------------------------------------------

! Parameters
! --------------------------------------------------------------------------------------------------
!  # | Variable     | Description                                   | Param. Units  | Internal Units
! --------------------------------------------------------------------------------------------------
!  1 | psymode      | 2 = Humidity ratio as humidity input          | -             | -
!                   | 4 = Relative humidity as humidity input       |               |
!  2 | PelcRated    | Rated total cooling power                     | kJ/h          | kJ/h
!  3 | QcRated      | Rated cooling capacity                        | kJ/h          | kJ/h
!  4 | PelhRated    | Rated total heating power                     | kJ/h          | kJ/h
!  5 | QhRated      | Rated heating capacity                        | kJ/h          | kJ/h
!  6 | AFRrated     | Rated inlet air mass flow rate                | kg/h          | kg/h
!  7 | freqRatedCool| Rated cooling frequency                       | 1/s           | 1/s
!  8 | freqRatedHeat| Rated heating frequency                       | 1/s           | 1/s
!  9 | PfanRatedIn  | Indoor fan rated power                        | kJ/h          | kJ/h
! 10 | PfanRatedOut | Outdoor fan rated power                       | kJ/h          | kJ/h
! 11 | backupHeat   | Backup heater capacity                        | kJ/h          | kJ/h
! 12 | LUcool       | Logical Unit - cooling mode                   | -             | -
! 13 | LUheat       | Logical Unit - heating mode                   | -             | -
! --------------------------------------------------------------------------------------------------

! Outputs
! --------------------------------------------------------------------------------------------------
!  # | Variable     | Description                                   | Output  Units | Internal Units
! --------------------------------------------------------------------------------------------------
!  1 | Ts           | Outlet (supply) air temperature               | °C            | °C
!  2 | ws           | Outlet (supply) air humidity ratio            | -             | -
!  3 | RHs          | Outlet (supply) air % RH                      | % (base 100)  | % (base 100)
!  4 | ps           | Outlet (supply) air pressure                  | atm           | atm
!  5 | mDot         | Outlet (supply) air mass flow rate            | kg/h          | kg/h
!  6 | Qc           | Total cooling rate                            | kJ/h          | kJ/h
!  7 | Qcs          | Sensible cooling rate                         | kJ/h          | kJ/h
!  8 | Qcl          | Latent cooling rate                           | kJ/h          | kJ/h
!  9 | Qrej         | Heat rejection rate                           | kJ/h          | kJ/h
! 10 | Qh           | Total heating rate                            | kJ/h          | kJ/h
! 11 | Qabs         | Heat absorption rate                          | kJ/h          | kJ/h
! 12 | Pel          | Total power consumption                       | kJ/h          | kJ/h
! 13 | COP          | Coefficient of performance                    | -             | -
! 14 | EER          | Energy efficiency rating                      | -             | -
! 15 | PfanIn       | Indoor fan power                              | kJ/h          | kJ/h
! 16 | PfanOut      | Outdoor fan power                             | kJ/h          | kJ/h
! 17 | Pcomp        | Compressor power                              | kJ/h          | kJ/h
! 18 | fComp        | Compressor frequency                          | 1/s           | 1/s
! 19 | Tc           | Condensate temperature                        | °C            | °C
! 20 | cmfr         | Condensate mass flow rate                     | kg/h          | kg/h
! 21 | defrost_mode | 0 = defrost (off) mode                        | -             | -
!                   | 1 = Recovery mode (transient)                 |               |
!                   | 2 = Steady-state mode                         |               |
! 22 | shutdown     | 0 = No (false)                                | -             | -
!                   | 1 = Yes (true)                                |               |
! --------------------------------------------------------------------------------------------------

module Type3254Data

use, intrinsic :: iso_fortran_env, only : wp=>real64    ! Defines a constant "wp" (working precision) that can be used in real numbers, e.g. 1.0_wp, and sets it to real64 (double precision)
implicit none

type Type3254DataStruct

    ! Parameters
    real(wp), allocatable :: entries(:, :, :), lbounds(:, :), ubounds(:, :)
    integer, allocatable :: extents(:, :)
    integer :: PMClength, PMHlength  ! Length of the flattened performance map

    ! Performance matrices
    real(wp), allocatable :: PelcMap(:, :, :, :, :)
    real(wp), allocatable :: QcsMap(:, :, :, :, :)
    real(wp), allocatable :: QclMap(:, :, :, :, :)
    real(wp), allocatable :: PelhMap(:, :, :, :)
    real(wp), allocatable :: QhMap(:, :, :, :)

end type Type3254DataStruct

type(Type3254DataStruct), allocatable, save :: s(:)

end module Type3254Data

subroutine Type3254
!export this subroutine for its use in external DLLs
!DEC$Attributes DLLexport :: Type3254

use, intrinsic :: iso_fortran_env, only : wp=>real64    ! Defines a constant "wp" (working precision) that can be used in real numbers, e.g. 1.0_wp, and sets it to real64 (double precision)

use TrnsysConstants
use TrnsysFunctions
use Type3254Data

implicit none

integer :: thisUnit, thisType  ! unit and type numbers
real(wp) :: time, dt  ! TRNSYS time and timestep

! Proforma variables
real(wp) :: Tr, wr, RHr, mDot, AFR, pr, Toa, woa, RHoa, poa, freq  ! Inputs
integer :: mode, defrost_mode = 1  ! assume that the heat pump starts up at the beginning -> start with transient state
integer :: psymode, LUcool, LUheat  ! Parameters
real(wp) :: PelcRated, QcRated, PelhRated, QhRated, AFRrated ! Parameters (rated values)
real(wp) :: freqRatedCool, freqRatedHeat, PfanRatedIn, PfanRatedOut, backupHeat
real(wp) :: Ts, ws, RHs, ps  ! Outputs (supply conditions)
real(wp) :: Pel, Qc, Qcs, Qcl, Qrej, Qh, Qabs, Pcomp, PfanIn, PfanOut  ! Outputs (heat and power)
real(wp) :: fComp, COP, EER, Tc, cmfr, recov_penalty  ! Outputs (misc)

! Local variables
real(wp) :: psydat(9), Twbr, Twboa, hr, hx, hs, dr, shutdown_fraction, extrapolation_fraction
integer :: status, shutdowns, extrapolations, shutdown_alt
integer, parameter :: Ninstances = 1  ! Number of units (should be provided by a function)
integer :: Ni = 1  ! temporary, should use a kernel function to get the actual instance number.
real(wp) :: defrost_corr(2) = 0.0_wp  ! defrost correction factors
logical :: shutdown, previous_shutdown, defrost_cooling
Character (len=maxMessageLength) aString, bString

! Performance map reading variables
integer, parameter :: Nc = 5, Nh = 4  ! Number of interpolation variables
integer, parameter :: Nmax = max(Nc, Nh)
integer :: i, j, N, Nd, Noutc = 3, Nouth = 2, Nout

! Interpolation variables
real(wp), allocatable :: interpolation_results(:), point(:)

! Set the version number for this Type
if ( GetIsVersionSigningTime() ) then
    call SetTypeVersion(18)
    return
endif

time = GetSimulationTime()
dt = GetSimulationTimeStep()
thisUnit = GetCurrentUnit()
thisType = GetCurrentType()

! All the stuff that must be done once at the beginning
if ( GetIsFirstCallofSimulation() ) then
    call ExecuteFirstCallOfSimulation()
    return
endif

! Parameters must be re-read - indicates another unit of this Type
if ( GetIsReReadParameters() ) call ReadParameters()

! Start of the first timestep: no iterations, outputs initial conditions
if ( GetIsStartTime() ) then
    call ExecuteStartTime()
	return
endif

! End of timestep call (after convergence or too many iterations)
if ( GetIsEndOfTimestep() ) then
    call ExecuteEndOfTimestep()
    return
endif

if ( GetIsLastCallofSimulation() ) then
    call ExecuteLastCallOfSimulation()
    return
endif

call GetInputValues()
fComp = freq * ((1-mode) * freqRatedCool + mode * freqRatedHeat)
shutdown = .false.
if (ErrorFound()) return

! Ni = GetCurrentUnit()

defrost_cooling = (mode == 1) .and. (defrost_mode == 0)
N = Nc * (1 - mode) + Nh * mode  ! Number of interpolation variables
if ( mode == 0 .or. defrost_cooling ) then
    Nout = Noutc  ! Number of output values for the interpolation
    Nd = Nc  ! For allocations in defrost mode
else
    Nout = Nouth
    Nd = Nh
end if

! Determine return air state
psydat(1) = pr
psydat(2) = Tr
psydat(4) = RHr/100.0_wp
psydat(6) = wr
if ( mode==0 .or. defrost_cooling ) then
    call Psychro(thisUnit, thisType, 1, psymode, 1, psydat, 1, status)
    ! (unit, type, si units used, psych inputs, Twb computed, inputs, warning mgmt, warning occurences)
else
    call Psychro(thisUnit, thisType, 1, psymode, 0, psydat, 1, status)  ! Twb not computed
end if
pr = psydat(1)
Tr = psydat(2)
Twbr = psydat(3)
RHr = psydat(4)  ! RHr between 0 and 1 (not 0 and 100)
wr = psydat(6)
hr = psydat(7)
dr = psydat(9)

if (AFR >= 0.0_wp) then
    mDot = AFR * AFRrated * dr  ! Use normalized AFR as input if it is positive
else
    AFR = mDot / (dr * AFRrated)
endif

PfanIn = PfanRatedIn * AFR**3  ! Adjust the fan power according to the flowrate
PfanOut = PfanRatedOut  ! Assume constant speed for outdoor fan

if (freq > 0.0_wp) then
    ! Interpolate using wet bulb in cooling
    allocate(point(N))
    if (mode == 1) then
        point = (/Tr, Toa, AFR, freq/)
        shutdown = CheckShutdown(point, 1, Nh)
        if ( defrost_cooling .and. .not. shutdown ) then
            deallocate(point)
            allocate(point(Nc))
        end if
        defrost_corr = Correction(defrost_mode, recov_penalty)
    end if

    if ( (mode == 0) .or. (defrost_cooling .and. .not. shutdown) ) then
        point = (/Tr, Twbr, Toa, AFR, freq/)
        shutdown = CheckShutdown(point, 0, Nc)
    end if

    previous_shutdown = GetOutputValue(22) == 1.0_wp
    if ( previous_shutdown .neqv. shutdown ) shutdown_alt = GetOutputValue(23) + 1
    if (shutdown_alt > 4) shutdown = previous_shutdown

    if (shutdown) then
        call ForceZeroPerformance()
        call SetDynamicArrayValueThisIteration(1, GetDynamicArrayValueLastTimestep(1) + 1.0_wp)
        fcomp = 0.0_wp  ! Force shutdown -> set output frequency to zero
        if ( (mode == 1) .and. (freq > 0.0_wp) ) then  ! Only activate auxiliary when required by the controller
            Qh = backupHeat
            Pel = backupHeat
        end if
    else
        allocate(interpolation_results(Noutc + Nouth - 1))

        if (defrost_cooling) then
            interpolation_results = Interpolate(point, 0, Nout)
            if (OutOfRange(point, 0)) call SetDynamicArrayValueThisIteration(2, GetDynamicArrayValueLastTimestep(2) + 1.0_wp)
        else
            interpolation_results = Interpolate(point, mode, Nout)
            if (OutOfRange(point, mode)) call SetDynamicArrayValueThisIteration(2, GetDynamicArrayValueLastTimestep(2) + 1.0_wp)
        end if

        if ( mode == 0 .or. defrost_cooling ) then
            Pel = interpolation_results(1) * PelcRated
        else
            Pel = interpolation_results(1) * PelhRated * defrost_corr(1)
        end if

        Qcs = interpolation_results(2) * QcRated
        Qcl = interpolation_results(3) * QcRated
        Qh = interpolation_results(4) * QhRated * defrost_corr(2)
        Qc = Qcs + Qcl

        deallocate(interpolation_results)
    end if

    deallocate(point)
else
    call ForceZeroPerformance()
endif

! Determine supply air state

ps = pr  ! Fan pressure drop neglected

! Moist air state
if (Qc < Qcs) then
    Qc = Qcs
    ! Add warning
endif
if (mDot /= 0.0_wp) then
    ws = wr
    if (mode == 0) then
        hs = hr - Qc/mDot
        hx = hr  ! useful when the following if clause is not true
        if (Qcl > 0.0_wp) then  ! compute humidity after condensation
            psydat(1) = pr
            psydat(2) = Tr
            hx = hr - Qcl/mDot
            psydat(7) = hx  ! enthalpy of the state (Tr, ws)
            call Psychro(thisUnit, thisType, 1, 5, 0, psydat, 1, status)  ! dry-bulb and enthalpy as inputs
            if (ErrorFound()) return
            ws = psydat(6)
        endif
    else
        hs = hr + Qh/mDot
    end if
else
    hs = hr
    ws = wr
    hx = hr
endif

psydat(1) = ps
psydat(6) = ws
psydat(7) = hs
call Psychro(thisUnit, thisType, 1, 7, 0, psydat, 1, status)  ! humidity ratio and enthalpy as inputs
if (ErrorFound()) return
ps = psydat(1)
Ts = psydat(2)
RHs = psydat(4)
ws = psydat(6)
hs = psydat(7)

! Re-calculate heat transfer whose value is modified if saturation occurs
if ( mode == 0 .and. freq > 0.0_wp ) then
    Qcs = mDot * (hx - hs)  ! Sensible cooling rate
    Qcl = mDot * (hr - hx)  ! Latent cooling rate
    Qc = Qcs + Qcl  ! Total cooling rate
    Qrej = Qc + Pel  ! Heat rejection in ambient air
    Qabs = 0.0_wp
else if (freq > 0.0_wp) then
    Qrej = 0.0_wp
    Qabs = Qh - Pel  ! Heat absorption in ambient air
end if

if (fcomp > 0.0_wp) then
    Pcomp = max(0.0_wp, Pel - PfanIn - PfanOut)  ! Compressor power
else
    Pcomp = 0.0_wp
end if

! COP, EER and condensate
if (Pel /= 0.0_wp) then
    COP = (Qc + Qh) / Pel
else
    COP = 0.0_wp
endif
EER = 3.413_wp * COP
Tc = Ts
cmfr = mDot * (wr - ws)  ! Condensate flow rate - water balance

call SetOutputValues()

return

    contains

    subroutine ReadPermap(LUc, LUh)
        integer, intent(in) :: LUc, LUh
        character (len=maxPathLength) :: permapCoolPath
        character (len=maxPathLength) :: permapHeatPath
        integer :: i, j, LUs(2), LUcool(1), LUheat(1)
        integer :: nTr, nTwbr, nToa, nAFR, nfreq  ! number of entries for each variable
        real(wp) :: filler(Nmax)
        LUcool(1) = LUc
        LUheat(1) = LUh

        ! Ni = GetCurrentUnit()
        LUs = (/LUc, LUh/)

        permapCoolPath = GetLUfileName(LUc)
        permapHeatPath = GetLUfileName(LUh)
        call CheckPMfile(permapCoolPath)
        call CheckPMfile(permapHeatPath)
        if (ErrorFound()) return

        open(LUc, file=permapCoolPath, status='old')
        open(LUh, file=permapHeatPath, status='old')

            call SkipLines(LUs, 6)
        allocate(s(Ni)%extents(Nmax, 0:1))
        allocate(s(Ni)%lbounds(Nmax, 0:1))
        allocate(s(Ni)%ubounds(Nmax, 0:1))
        do i = 1, Nc
            call SkipLines(LUcool, 1)
            read(LUc, *) s(Ni)%extents(i, 0), s(Ni)%lbounds(i, 0), s(Ni)%ubounds(i, 0)
        end do
        do i = 1, Nh
            call SkipLines(LUheat, 1)
            read(LUh, *) s(Ni)%extents(i, 1), s(Ni)%lbounds(i, 1), s(Ni)%ubounds(i, 1)
        end do
        s(Ni)%PMClength = product(s(Ni)%extents(1:Nc, 0))
        s(Ni)%PMHlength = product(s(Ni)%extents(1:Nh, 1))
        allocate(s(Ni)%entries(maxval(s(Ni)%extents), Nmax, 0:1))
        do i = 1, Nc
            call SkipLines(LUcool, 1)
            read(LUc, *) (s(Ni)%entries(j, i, 0), j = 1, s(Ni)%extents(i, 0))
        end do
        do i = 1, Nh
            call SkipLines(LUheat, 1)
            read(LUh, *) (s(Ni)%entries(j, i, 1), j = 1, s(Ni)%extents(i, 1))
        end do
            call SkipLines(LUs, 4)

        nTr = s(Ni)%extents(1, 0)
        nTwbr = s(Ni)%extents(2, 0)
        nToa = s(Ni)%extents(3, 0)
        nAFR = s(Ni)%extents(4, 0)
        nfreq = s(Ni)%extents(5, 0)
        allocate(s(Ni)%PelcMap(nTr, nTwbr, nToa, nAFR, nfreq))
        allocate(s(Ni)%QcsMap(nTr, nTwbr, nToa, nAFR, nfreq))
        allocate(s(Ni)%QclMap(nTr, nTwbr, nToa, nAFR, nfreq))
        do i = 1, s(Ni)%PMClength
            read(LUc, *) (filler(j), j = 1, Nc), Pel, Qcs, Qcl
            call SetPMvalue(s(Ni)%PelcMap, RowToColMajorOrder(i, s(Ni)%extents(1:Nc, 0)), Pel, s(Ni)%PMClength)
            call SetPMvalue(s(Ni)%QcsMap, RowToColMajorOrder(i, s(Ni)%extents(1:Nc, 0)), Qcs, s(Ni)%PMClength)
            call SetPMvalue(s(Ni)%QclMap, RowToColMajorOrder(i, s(Ni)%extents(1:Nc, 0)), Qcl, s(Ni)%PMClength)
        end do

        close(LUc)

        nTr = s(Ni)%extents(1, 1)
        nToa = s(Ni)%extents(2, 1)
        nAFR = s(Ni)%extents(3, 1)
        nfreq = s(Ni)%extents(4, 1)
        allocate(s(Ni)%PelhMap(nTr, nToa, nAFR, nfreq))
        allocate(s(Ni)%QhMap(nTr, nToa, nAFR, nfreq))
        do i = 1, s(Ni)%PMHlength
            read(LUh, *) (filler(j), j = 1, Nh), Pel, Qh
            call SetPMvalue(s(Ni)%PelhMap, RowToColMajorOrder(i, s(Ni)%extents(1:Nh, 1)), Pel, s(Ni)%PMHlength)
            call SetPMvalue(s(Ni)%QhMap, RowToColMajorOrder(i, s(Ni)%extents(1:Nh, 1)), Qh, s(Ni)%PMHlength)
        end do

        close(LUh)

    end subroutine ReadPermap

    subroutine CheckPMfile(permapPath)
        logical :: permapFileFound = .false.
        character (len=maxPathLength) :: permapPath
        character (len=maxMessageLength) :: msg
        inquire(file=trim(permapPath), exist=permapFileFound)
        if ( .not. permapFileFound ) then
            write(msg, '("""",a,"""")') trim(permapPath)
            msg = "Could not find the specified performance map file. Searched for: " // trim(msg)
            call Messages(-1, msg, 'fatal', thisUnit, thisType)
            return
        end if
    end subroutine CheckPMfile

    subroutine SkipLines(LUs, N)
    ! Skip lines in several files at once
    !
    ! Inputs
    !   LUs (integer array) : logical unit of each file
    !                         where lines must be skipped.
    !   N (integer) : number of lines to skip.
        integer, intent(in) :: LUs(:), N
        integer :: i, j
        do i = 1, size(LUs)
            do j = 1, N
                read(LUs(i), *)
            end do
        end do
    end subroutine SkipLines

    function RowToColMajorOrder(rowIndex, extents) result(colIndex)
    ! Transform a row-major order index corresponding to a given array shape
    ! into a column-major index.
    !
    ! Inputs
    !   rowIndex (integer) : one-based index of an array in row-major order.
    !   extents (integer array) : shape of the array that rowIndex is indexing.
    !
    ! Output
    !   colIndex (integer) : one-based index corresponding to the array element
    !                        indexed by rowIndex, but with a column-major order.
        integer, intent(in) :: extents(:), rowIndex
        integer :: i, colIndex, j, p

        i = rowIndex - 1  ! rowIndex is one-based
        colIndex = 1  ! colIndex is one-based
        p = product(extents)
        do j = size(extents), 1, -1
            p = p / extents(j)  ! ("/" performs integer division)
            colIndex = colIndex + p * modulo(i, extents(j))
            i = i / extents(j)
        end do
    end function RowToColMajorOrder

    function OutOfRange(conditions, mode)
    ! Check whether the input conditions are within the operating ranges
    ! specified in the performance map.
    !
    ! Inputs
    !   conditions (real(wp) array) : input conditions.
    !   mode (integer) : operating mode, cooling (0) or heating (1)
    !
    ! Output
    !   inRange (logical) : true if all conditions are within their operating range,
    !                       false otherwise.
        real(wp), intent(in) :: conditions(N)
        integer, intent(in) :: mode
        integer :: i
        logical :: outOfRange

        outOfRange = .false.
        do i = 1, N
            if ( (point(i) > s(Ni)%ubounds(i, mode)) .or. (point(i) < s(Ni)%lbounds(i, mode)) ) then
                outOfRange = .true.
                exit
            end if
        end do
    end function OutOfRange

    function CheckShutdown(conditions, mode, N) result(shutdown)
    ! Check whether the heat pump should be forced to shut down
    ! because of invalid input conditions.
    !
    ! Inputs
    !   conditions (real(wp) array) : input conditions.
    !   mode (integer) : operating mode, cooling (0) or heating (1)

        integer, intent(in) :: mode, N
        real(wp), intent(in) :: conditions(N)
        integer :: i
        logical :: increasing(N), shutdown

        if (mode == 0) increasing = (/.true., .false., .false., .true., .true./)
        if (mode == 1) increasing = (/.false., .true., .true., .true./)
        shutdown = .false.
        do i = 1, N
            if ( increasing(i) .and. (point(i) < s(Ni)%lbounds(i, mode)) .or. &
            (.not. increasing(i)) .and. (point(i) > s(Ni)%ubounds(i, mode)) ) then
                shutdown = .true.
                exit
            end if
        end do
    end function CheckShutdown

    function Interpolate(point, mode, Nout) result(interpolation)
    ! Perform a piecewise multilinear interpolation on the tables in
    ! the structure Type3254DataStruct.
    !
    ! Inputs
    !   point (real(wp) array) : point where the interpolation has to be performed.
    !   mode (integer) : operating mode, cooling (0) or heating (1)
    !   Nout (integer) : number of tables for the interpolation,
    !                    also the number of required output values.
    !
    ! Output
    !   interpolation (real(wp) array) : results of the interpolation.
        real(wp), intent(in) :: point(Nd)
        real(wp) :: scaled_point(Nd), left, right, sp
        real(wp), allocatable :: hypercube(:, :)
        integer, intent(in) :: mode
        integer, dimension(Nd) :: idx, lb_idx, counter_int
        integer :: i
        logical :: counter_bool(Nd)
        integer, intent(in) :: Nout
        real(wp) :: interpolation(Noutc + Nouth - 1)

        ! Map the point to the unit N-hypercube with the table values bounding the point
        do i = 1, Nd
            ! Find the index of the lower bound for the ith component of point
            j = Findlb(s(Ni)%entries(:, i, mode), point(i), s(Ni)%extents(i, mode))
            lb_idx(i) = j
            left = s(Ni)%entries(j, i, mode)
            right = s(Ni)%entries(j+1, i, mode)
            ! Compute the scaled point, and keep it between 0 and 1
            scaled_point(i) = min(1.0_wp, max(0.0_wp, (point(i) - left) / (right - left)))
        end do
        allocate(hypercube(Nout, 2**Nd))
        counter_bool = .true.  ! All N values of the binary counter are initialized at 1
        do i = 1, 2**Nd
            call Increment(counter_bool)
            counter_int = merge(1, 0, counter_bool)  ! convert logical to int
            idx = lb_idx + counter_int  ! index of each of the 2**N vertices
            hypercube(:, i) = Vertex(idx, mode)  ! Get table values associated with the vertices
        end do

        ! Interpolate using the scaled point and the table values
        do i = 1, Nd
            sp = scaled_point(i)
            j = Nd - i
            hypercube(:, :2**j) = (1-sp) * hypercube(:, :2**j) &
                                    + sp * hypercube(:, 2**j+1:2**(j+1))
        end do

        if (mode == 0) then  ! cooling -> fill the first Noutc values, the rest are zeros
            do i = 1, Noutc
                interpolation(i) = hypercube(i, 1)
            end do
            do i = Noutc+1, Noutc+Nouth-1
                interpolation(i) = 0.0_wp
            end do
        else  ! heating -> first value is still the input power, values 2 through Noutc are zeros
            interpolation(1) = hypercube(1, 1)
            do i = 2, Noutc
                interpolation(i) = 0.0_wp
            end do
            do i = Noutc+1, Noutc+Nouth-1
                interpolation(i) = hypercube(i-Noutc+1, 1)
            end do
        end if

        deallocate(hypercube)
    end function Interpolate

    function Vertex(idx, mode)
    ! Provide the performance map values associated with a certain index,
    ! taking the operating mode into account.
    !
    ! Inputs
    !   idx (integer array) : index of the table values to be retrieved.
    !   mode (integer) : operating mode, cooling (0) or heating (1).
    !
    ! Output
    !   vertex (real(wp) array) : array with the table values, given in the form
    !                             (power, sensible capacity, latent capacity) in
    !                             cooling and (power, total capacity) in heating.
        integer, intent(in) :: idx(:), mode
        real(wp) :: Pel, Qcs, Qcl, Qh, vertex(Nout)
        if (mode == 0) then
            Pel = GetPMvalue(mode, s(Ni)%PelcMap, idx)
            Qcs = GetPMvalue(mode, s(Ni)%QcsMap, idx)
            Qcl = GetPMvalue(mode, s(Ni)%QclMap, idx)
            vertex = (/Pel, Qcs, Qcl/)
        else
            Pel = GetPMvalue(mode, s(Ni)%PelhMap, idx)
            Qh = GetPMvalue(mode, s(Ni)%QhMap, idx)
            vertex = (/Pel, Qh/)
        end if
    end function Vertex

    function GetPMvalue(mode, array, idx) result(value)
    ! Retrieve the value at a given index in a performance map.
    !
    ! Inputs
    !   mode (integer, 0 or 1) : mode (heating or cooling) associated with the performance map.
    !   array : performance map, contained in the structure Type3254DataStruct.
    !   idx (integer array with the same length as the array shape) : index of the element to retrieve.
    !
    ! Output
    !   value (real(wp)) : element retireved at index idx.
        integer, intent(in) :: idx(:)
        integer :: mode, i, array_idx
        ! array has a length equal to the number of elements in the appropriate performance map
        real(wp) :: array((1-mode)*s(Ni)%PMClength + mode*s(Ni)%PMHlength)
        real(wp) :: value
        array_idx = idx(1)  ! begin with the leftmost index, since arrays are stored in column-major order
        do i = 2, Nd
            array_idx = array_idx + product(s(Ni)%extents(1:i-1, mode)) * (idx(i) - 1)
        end do
        value = array(array_idx)
    end function GetPMvalue

    subroutine SetPMvalue(array, idx, value, PMlength)
    ! Set a value in a performance map at a given index.
    !
    ! Inputs
    !   array : performance map, contained in the structure Type3254DataStruct.
    !   idx (integer) : one-dimensional index of the element to set, in column-major order.
    !   value (real(wp)) : the value to set at index idx.
    !   PMlength (integer) : length of the performance map.
        integer, intent(in) :: idx, PMlength
        real(wp) :: array(PMlength)
        real(wp), intent(in) :: value
        array(idx) = value
    end subroutine SetPMvalue

    function Findlb(array, value, extent) result(lowerBoundIndex)
    ! Find the index of the element of an ordered array which is smaller
    ! and closest to a given value.
    ! Special cases: If the value is exactly equal to an element of the array
    ! with index i, the function returns i-1. If the value is smaller than the first
    ! (and lowest) value of the array, the function returns 1. If the value is higher
    ! than the last (and largest) value of the array, the function returns size(array) - 1.
    !
    ! Inputs
    !   array (1-D real(wp) array) : array to be searched.
    !   value (real(wp)) : value to search.
    !   extent (integer) : size of the array
    !
    ! Output
    !   lowerBoundIndex (integer) : index of the lower bound of the value in the array.
        real(wp), intent(in) :: array(:)
        real(wp), intent(in) :: value
        integer, intent(in) :: extent
        integer :: L, R, mid, lowerBoundIndex
        L = 1
        R = extent
        do while (L < R)
            mid = (L + R) / 2  ! L & R are integers -> automatic floor
            if (array(mid) < value) then
                L = mid + 1
            else
                R = mid
            end if
        end do
        lowerBoundIndex = L - 1
        if (lowerBoundIndex == 0) lowerBoundIndex = 1
    end function Findlb

    function FullAdder(a, b, carry_in) result(results)
    ! Implement a full adder, to perform addition on binary numbers.
    !
    ! Inputs
    !   a, b (logicals) : operands of the addition
    !   carry_in (logical) : carry from the previous stage.
    !
    ! Outputs (given as a logical array of size 2)
    !   sum (logical) : a + b + carry_in (+ being exclusive or).
    !   carry_out (logical) : carry for the next stage.
        implicit none
        logical, intent(in) :: a, b, carry_in
        logical :: sum, carry_out, results(2)
        sum = a .neqv. b .neqv. carry_in
        carry_out = a .and. b .or. carry_in .and. (a .neqv. b)
        results = (/sum, carry_out /)
    end function FullAdder

    subroutine Increment(C)
    ! Increment a binary counter.
    !
    ! Input
    !   C (logical array) : the counter to increment,
    !                       with the rightmost bit being the lower level.
        implicit none
        logical, intent(inout) :: C(:)
        logical :: sumcarry(2)
        integer :: N, k, i
        N = size(C)
        if ( all(C) ) then  ! counter has reached max capacity
            C = .false.  ! reset counter
        else  ! increment counter using binary addition
            sumcarry = FullAdder(C(N), .true., .false.)
            C(N) = sumcarry(1)
            k = N - 1
            do while (sumcarry(2))
                sumcarry = FullAdder(C(k), .false., sumcarry(2))
                C(k) = sumcarry(1)
                k = k-1
            end do
        end if
    end subroutine Increment

    function Correction(defrost_mode, recov_penalty)
    ! Provide the correction factor to apply to the input power
    ! and heating capacity depending on the defrost mode.
    !
    ! Inputs
    !   defrost_mode (integer) : the defrost mode: defrost(0), recovery (1)
    !                            or steady-state (2).
    !   recov_penalty (real(wp)) : the correction factor for recovery mode.
    !
    ! Output
    !   correction (real(wp) array) : array containing the power correction
    !                                 factor at (1) and the capacity correction
    !                                 factor at (2).
        integer, intent(in) :: defrost_mode
        real(wp), intent(in) :: recov_penalty
        real(wp) :: correction(2)
        if (defrost_mode == 0) then
            correction = (/0.6_wp, 0.0_wp/)
        else if (defrost_mode == 1) then
            correction = (/1.0_wp, recov_penalty/)
        else if (defrost_mode == 2) then
            correction = (/1.0_wp, 1.0_wp/)
        else
            ! add warning
        end if
    end function Correction

    subroutine ForceZeroPerformance
        Qcs = 0.0_wp
        Qcl = 0.0_wp
        Qh = 0.0_wp
        Qc = 0.0_wp
        Pel = 0.0_wp
    end subroutine ForceZeroPerformance

    subroutine Psychro(thisUnit, thisType, iUnits, psymode, wbmode, psydat, eMode, stat)
    ! Wrapper for the routine MoistAirProperties (see documentation for inputs/outputs).
    ! If the moist air state is determined from the humidity ratio and enthalpy (i.e. psymode = 7),
    ! Psychro checks that the state is not over the saturation curve, as that can cause the
    ! MoistAirProperties to loop indefinitely. If it is, humidity ratio is reset to saturation
    ! humidity ratio at the given value of enthalpy.
        integer :: thisUnit, thisType, iUnits, psymode, wbmode, eMode, stat
        real(wp) :: psydat(9), satdat(9), hsat
        if (psymode == 7) then  ! Air state determined with humidity ratio (#6) and enthalpy (#7)
            satdat = psydat
            ! Setting RH too close to 1 may cause convergence problems in MoistAirProperties,
            ! e.g. when RH = 0.999 and h = 9.4303316807275.
            satdat(4) = 0.99_wp  ! Go to saturation curve following isenthalpic line (set RH=100%)
            call MoistAirProperties(thisUnit, thisType, iUnits, 8, 0, satdat, eMode, status)
            if (satdat(6) < psydat(6)) then  ! Over saturation curve
                psydat = satdat  ! Keep saturated state with same enthalpy.
            else
                call MoistAirProperties(thisUnit, thisType, iUnits, psymode, wbmode, psydat, eMode, status)
            end if
        else
            call MoistAirProperties(thisUnit, thisType, iUnits, psymode, wbmode, psydat, eMode, status)
        end if
    end subroutine Psychro

    subroutine ExecuteFirstCallOfSimulation
  	    call SetNumberofParameters(13)
  	    call SetNumberofInputs(14)
  	    call SetNumberofDerivatives(0)
  	    call SetNumberofOutputs(23)
  	    call SetIterationMode(1)
  	    call SetNumberStoredVariables(0, 2)
  	    call SetNumberofDiscreteControls(0)

        ! Allocate stored data structure
        if ( .not. allocated(s) ) then
            allocate(s(Ninstances))
        endif

        call ReadParameters()  ! required to get LUcool and LUheat
        call ReadPermap(LUcool, LUheat)

    end subroutine ExecuteFirstCallOfSimulation

    subroutine ExecuteStartTime
        call ReadParameters()
        call GetInputValues()
        call SetOutputValue(1, 0.0_wp)  ! Outlet air temperature
        call SetOutputValue(2, 0.0_wp)  ! Outlet air humidity ratio
        call SetOutputValue(3, 0.0_wp)  ! Outlet air % RH
        call SetOutputValue(4, 0.0_wp)  ! Outlet air pressure
        call SetOutputValue(5, 0.0_wp)  ! Outlet air flow rate
        call SetOutputValue(6, 0.0_wp)  ! Total cooling rate
        call SetOutputValue(7, 0.0_wp)  ! Sensible cooling rate
        call SetOutputValue(8, 0.0_wp)  ! Latent cooling rate
        call SetOutputValue(9, 0.0_wp)  ! Heat rejection rate
        call SetOutputValue(10, 0.0_wp)  ! Total heating rate
        call SetOutputValue(11, 0.0_wp)  ! Heat absorption rate
        call SetOutputValue(12, 0.0_wp)  ! Total power consumption
        call SetOutputValue(13, 0.0_wp)  ! COP
        call SetOutputValue(14, 0.0_wp)  ! EER
        call SetOutputValue(15, 0.0_wp)  ! Indoor fan power
        call SetOutputValue(16, 0.0_wp)  ! Outdoor fan power
        call SetOutputValue(17, 0.0_wp)  ! Compressor power
        call SetOutputValue(18, 0.0_wp)  ! Compressor frequency
        call SetOutputValue(19, 0.0_wp)  ! Condensate temperature
        call SetOutputValue(20, 0.0_wp)  ! Condensate flow rate
        call SetOutputValue(21, 0.0_wp)  ! Defrost mode
        call SetOutputValue(22, 0.0_wp)  ! Force shutdown
        call SetOutputValue(23, 0.0_wp)  ! Force shutdown alternations

        call SetDynamicArrayInitialValue(1, 0.0_wp)  ! Force shutdowns counter
        call SetDynamicArrayInitialValue(2, 0.0_wp)  ! Constant extrapolations counter
    end subroutine ExecuteStartTime

    subroutine ExecuteEndOfTimestep
        shutdown_alt = 0
    end subroutine ExecuteEndOfTimestep

    subroutine ExecuteLastCallOfSimulation
        deallocate(s(Ni)%extents)
        deallocate(s(Ni)%entries)
        deallocate(s(Ni)%lbounds)
        deallocate(s(Ni)%ubounds)
        deallocate(s(Ni)%PelcMap)
        deallocate(s(Ni)%PelhMap)
        deallocate(s(Ni)%QcsMap)
        deallocate(s(Ni)%QclMap)
        deallocate(s(Ni)%QhMap)

        ! Force shutdowns warning
        shutdowns = int(GetDynamicArrayValueLastTimestep(1))
        shutdown_fraction = 100.0_wp * shutdowns / real(GetNTimeSteps() - 1, wp)
        if (shutdowns > 0) then
            write(aString,'(g)') shutdowns
            write(bString,'(g)') shutdown_fraction
            aString = 'Heat pump was forced to shut down  for'&
                //trim(adjustl(aString))//' time steps ( '//trim(adjustl(bString))//'% of the simulation)'

            if (shutdown_fraction < 5) Then
                call Messages(-1, aString, 'Notice', thisUnit, thisType)
            else
                call Messages(-1, aString, 'Warning', thisUnit, thisType)
            end if
        end if

        ! Constant extrapolations warning
        extrapolations = int(GetDynamicArrayValueLastTimestep(2))
        extrapolation_fraction = 100.0_wp * extrapolations / real(GetNTimeSteps() - 1, wp)
        if (extrapolations > 0) then
            write(aString,'(g)') extrapolations
            write(bString,'(g)') extrapolation_fraction
            aString = 'Performance was extrapolated using the closest values in the performance map for '&
                //trim(adjustl(aString))//' time steps ('//trim(adjustl(bString))//'% of the simulation)'

            if (extrapolation_fraction < 5.0_wp) then
                call Messages(-1, aString, 'Notice', thisUnit, thisType)
            else
                call Messages(-1, aString, 'Warning', thisUnit, thisType)
            end if
        end if

    end subroutine ExecuteLastCallOfSimulation

    subroutine ReadParameters
        psymode = GetParameterValue(1)
        PelcRated = GetParameterValue(2)
        QcRated = GetParameterValue(3)
        PelhRated = GetParameterValue(4)
        QhRated = GetParameterValue(5)
        AFRrated = GetParameterValue(6)
        freqRatedCool = GetParameterValue(7)
        freqRatedHeat = GetParameterValue(8)
        PfanRatedIn = GetParameterValue(9)
        PfanRatedOut = GetParameterValue(10)
        backupHeat = GetParameterValue(11)
        if (backupHeat < 0.0_wp) backupHeat = QhRated
        LUcool = GetParameterValue(12)
        LUheat = GetParameterValue(13)
    end subroutine ReadParameters

    subroutine GetInputValues
        Tr = GetInputValue(1)
        wr = GetInputValue(2)
        RHr = GetInputValue(3)
        pr = GetInputValue(4)
        mDot = GetInputValue(5)
        Toa = GetInputValue(6)
        woa = GetInputValue(7)
        RHoa = GetInputValue(8)
        poa = GetInputValue(9)
        freq = GetInputValue(10)
        AFR = GetInputValue(11)
        mode = GetInputValue(12)
        defrost_mode = GetInputValue(13)
        recov_penalty = GetInputValue(14)
    end subroutine GetInputValues

    subroutine SetOutputValues
        call SetOutputValue(1, Ts)  ! Outlet air temperature
        call SetOutputValue(2, ws)  ! Outlet air humidity ratio
        call SetOutputValue(3, RHs*100.0_wp)  ! Outlet air % RH
        call SetOutputValue(4, ps)  ! Outlet air pressure
        call SetOutputValue(5, mDot)  ! Outlet air flow rate
        call SetOutputValue(6, Qc)  ! Total cooling rate
        call SetOutputValue(7, Qcs)  ! Sensible cooling rate
        call SetOutputValue(8, Qcl)  ! Latent cooling rate
        call SetOutputValue(9, Qrej)  ! Heat rejection rate
        call SetOutputValue(10, Qh)  ! Total heating rate
        call SetOutputValue(11, Qabs)  ! Heat absorption rate
        call SetOutputValue(12, Pel)  ! Total power consumption
        call SetOutputValue(13, COP)  ! COP
        call SetOutputValue(14, EER)  ! EER
        call SetOutputValue(15, PfanIn)  ! Indoor fan power
        call SetOutputValue(16, PfanOut)  ! Outdoor fan power
        call SetOutputValue(17, Pcomp)  ! Compressor power
        call SetOutputValue(18, fComp)  ! Compressor frequency
        call SetOutputValue(19, Tc)  ! Condensate temperature
        call SetOutputValue(20, cmfr)  ! Condensate flow rate
        call SetOutputValue(21, real(defrost_mode, wp))  ! Defrost mode
        call SetOutputValue(22, merge(1.0_wp, 0.0_wp, shutdown))  ! Force shutdown
        call SetOutputValue(22, merge(1.0_wp, 0.0_wp, shutdown))  ! Force shutdown
        call SetOutputValue(23, shutdown_alt)  ! Force shutdown alternations
    end subroutine SetOutputValues

end subroutine Type3254

max_stars_repo_name: polymtl-bee/vcaahp-model

max_issues_repo_name: polymtl-bee/vcaahp-model

max_forks_repo_name: polymtl-bee/vcaahp-model


text:
ncpu 1

compactalignmentformat 1
coocurrencefile en_fr.cooc
corpusfile en_fr.snt
sourcevocabularyfile en.vcb
targetvocabularyfile fr.vcb

hmmiterations 5
model1iterations 5
model2iterations 0
model3iterations 0
model4iterations 5
model5iterations 0
model6iterations 0

hmmdumpfrequency 5
model1dumpfrequency 0
model2dumpfrequency 0
model345dumpfrequency 0

max_stars_repo_name: plandes/amrlib

max_issues_repo_name: plandes/amrlib

max_forks_repo_name: plandes/amrlib