text:
subroutine rdtddn(lundia    ,lunout    ,lunrd     ,error     ,filout    , &
                & fildis    ,runid     ,eol       ,itstrt    ,itfinish  , &
                & nsrc      ,lstsc     ,rval      ,maxval    ,namsrc    , &
                & disint    ,parnam    ,parunt    ,bubble    ,gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: rdtddn.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/io/src/preprocessor/rdtddn.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: This general purpose routine reads the time depen-
!              dent data from file (new)
! Method used:
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    !
    use globaldata
    use string_module
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    integer                    , pointer :: itdate
    real(fp)                   , pointer :: tstop
    real(fp)                   , pointer :: dt
    real(fp)                   , pointer :: tunit
    character*20, dimension(:) , pointer :: keywrd
    character*39, dimension(:) , pointer :: fmtdis
!
! Global variables
!
    integer                            , intent(in)  :: itfinish !  Description and declaration in inttim.igs
    integer                                          :: itstrt   !  Description and declaration in inttim.igs
    integer                            , intent(in)  :: lstsc    !  Description and declaration in dimens.igs
    integer                                          :: lundia   !  Description and declaration in inout.igs
    integer                                          :: lunout   !  Unit number for the transformed ouput file
    integer                                          :: lunrd    !  Unit number of the attribute file containing the time series
    integer                            , intent(in)  :: maxval   !  Maximum number of values 1+LSTSC+2
    integer                            , intent(in)  :: nsrc     !  Description and declaration in esm_alloc_int.f90
    logical                            , intent(in)  :: bubble   !  Description and declaration in procs.igs
    logical                                          :: error    !  Flag=TRUE if an error is encountered
    real(fp)     , dimension(maxval)                 :: rval     !  Array for the time dependent data
    character(*)                                     :: fildis   !  Name of time dependent input file
    character(*)                                     :: filout   !  Name of the output file to be opened
    character(*)                       , intent(in)  :: runid
    character(1)                       , intent(in)  :: eol      !  ASCII code for End-Of-Line (^J)
    character(1) , dimension(nsrc)                   :: disint   !  Description and declaration in esm_alloc_char.f90
    character(10), dimension(lstsc + 4), intent(in)  :: parunt   !  Unit name fitting the parameter
    character(20), dimension(nsrc)                   :: namsrc   !  Description and declaration in esm_alloc_char.f90
    character(36), dimension(lstsc + 4), intent(in)  :: parnam   !  Names of theparamaters to write to time dependent file DIS
!
! Local variables
!
    integer                      :: ibeg    ! Begin position in the RECORD from where the search for data/record is started
    integer                      :: iend    ! Last position in the RECORD when the searched data/record is finished
    integer                      :: ier     ! =  0 -> end of record encountered =  1 -> real value found = -1 -> length or number of data is larger then specified by the calling routine
    integer                      :: iocond  ! IO status for reading
    integer                      :: irec
    integer                      :: irecrd  ! Counter of records if input file is a direct access file
    integer                      :: itold   ! Help var. to store last read time to test accending order
    integer                      :: ittdep  ! Help var. for the time read (now de- fined as multiples of DT, but in fu- ture it may take any value)
    integer                      :: ix
    integer                      :: j
    integer                      :: l
    integer                      :: lflout  ! Actual length of FILOUT
    integer                      :: lrecrd  ! Length of record read
    integer                      :: mxlrec  ! Actual maximum record length
    integer                      :: n       ! Help var. for the nr. of data to be read (see LENDAT)
    integer                      :: nb
    integer                      :: newlun
    integer                      :: np
    integer                      :: npara   ! NR. of parameter records equal to number of data values in a record to be read
    integer                      :: nparrd  ! NR. of parameter records actual read
    integer                      :: nr
    integer                      :: nrval
    integer                      :: ntimrd
    logical                      :: access  ! Flag to read file as direct access or sequential
    logical                      :: ex
    real(fp)                     :: dtmin   ! Time step dt in minutes
    real(fp)                     :: rdfaul
    real(fp)                     :: timrd   ! Time in minutes read
    real(fp)                     :: timscl  ! Multiple factor to create minutes from read times
    character(1)                 :: quote   ! Apostrophe ASCII-character 39
    character(10)                :: disflg
    character(10)                :: interp  ! Character string containing inter- polation option
    character((4+lstsc)*16)      :: record  ! Record for data to be read. Maximum length is 11 * 24 + 1
    character(36), dimension(:), allocatable :: parrd   ! Parameter names read
    character(40)                :: cntain
    character(42)                :: tablnm  ! Table name specification
    character(300)               :: errmsg
!
!! executable statements -------------------------------------------------------
!
    fmtdis  => gdp%gdfmtdis%fmtdis
    keywrd  => gdp%gdkeywtd%keywrd
    itdate  => gdp%gdexttim%itdate
    tstop   => gdp%gdexttim%tstop
    dt      => gdp%gdexttim%dt
    tunit   => gdp%gdexttim%tunit
    !
    access = .false.
    irecrd = 0
    !
    quote  = char(39)
    !
    npara = 2 + lstsc + 2
    allocate(parrd(npara))
    !
    call remove_leading_spaces(filout    ,lflout    )
    !
    ! Determine the record length for direct access file and the
    ! corresponding format
    !
    mxlrec = max(125, 16 + (lstsc + 3)*14 + 1)
    if (mxlrec>125) then
       do nb = 1, mxfmtd
          ix = index(fmtdis(nb), 't125')
          if (ix/=0) write (fmtdis(nb)(ix + 1:ix + 4), '(i4)') mxlrec
       enddo
    endif
    !
    ! Open output file
    ! and write length direct access file to file first
    !
    lunout = newlun(gdp)
    inquire (file = filout(:lflout), exist = ex)
    if (ex) then
       open (lunout, file = filout(:lflout))
       close (lunout, status = 'delete')
    endif
    open (lunout, file = filout(:lflout), form = 'formatted', access = 'direct',&
        & status = 'unknown', recl = mxlrec)
    irec = 1
    write (lunout, fmtdis(1), rec = irec) '#', mxlrec, eol
    !
    ! Start reading the time dependent data.
    !
    do n = 1, nsrc
       ittdep = -1
       itold  = -1
       read (lunrd, '(a)', iostat = iocond) record
       !
       ! Premature EOR encountered or error (IOCOND <> 0)
       ! Test last time read
       ! Part of error message not fully correct
       !
       if (iocond/=0) then
          error = .true.
          call prterr(lundia    ,'G007'    ,fildis    )
          !
          if (iocond<0) then
             if (itold < itfinish) then
                write(errmsg,'(a,a,a)') 'Last time in file ', trim(fildis), ' <'
                call prterr(lundia    ,'U042'    ,errmsg    )
             elseif (n<nsrc) then
                call prterr(lundia    ,'U021'    ,'not for all discharge locations DIS data defined'    )
             else
             endif
          endif
       endif
       !
       ! Check for EOL in file and find maximum number of characters
       ! to read if EOL is located. Only for UNIX systems where EOL = ^J
       !
       lrecrd = len(record)
       if (eol/=' ' .and. index(record, eol)/=0) lrecrd = index(record, eol)
       !
       ! Read all keywords prior to time dependent data records
       ! Define number for data values (parameter names) to be read
       ! Input NPARA = 2 + LSTSC + 2, output can be -2 depending on
       ! contents of DISFLG
       !
       npara = 2 + lstsc + 2
       call flhnew(lunrd     ,lundia    ,error     ,record(:lrecrd)      ,access    , &
                 & irecrd    ,namsrc(n) ,disflg    ,disint(n) ,itdate    , &
                 & timscl    ,ntimrd    ,parrd     ,npara     ,nparrd    , &
                 & bubble    ,gdp       )
       if (error) then
          exit
       endif
       !
       ! Define NPARA
       !
       npara = nparrd
       nrval = npara - 1
       !
       ! Test string DISFLG
       !
       if (disflg=='regular   ') then
          if (nparrd/=2 + lstsc) then
             call prterr(lundia    ,'V097'    ,' '       )
             error = .true.
             exit
          endif
          cntain = 'regular    # discharge without momentum '
       elseif (disflg=='momentum  ') then
          if (nparrd/=4 + lstsc) then
             call prterr(lundia    ,'V097'    ,' '       )
             error = .true.
             exit
          endif
          cntain = 'momentum   # discharge with momentum'
       elseif (disflg=='walking   ') then
          if (nparrd/=2 + lstsc) then
             call prterr(lundia    ,'V097'    ,' '       )
             error = .true.
             exit
          endif
          cntain = 'walking   # walking discharge      '
       elseif (disflg=='inoutlet  ') then
          if (nparrd/=2 + lstsc) then
             call prterr(lundia    ,'V097'    ,' '       )
             error = .true.
             exit
          endif
          cntain = 'inoutlet  # power station / culvert'
       elseif (disflg=='power     ') then
          if (nparrd/=2 + lstsc) then
             call prterr(lundia    ,'V097'    ,' '       )
             error = .true.
             exit
          endif
          cntain = 'power     # power station in/out   '
       elseif (disflg=='culvert   ') then
          if (nparrd/=2 + lstsc) then
             call prterr(lundia    ,'V097'    ,' '       )
             error = .true.
             exit
          endif
          cntain = 'culvert   # culvert in/out         '
       else
          call prterr(lundia    ,'V096'    ,disflg    )
          error = .true.
          exit
       endif
       !
       ! Test name of constituent conform read parameter name
       ! Only first 20 characters are of significance
       !
       do np = 2, nparrd
          if (parrd(np)(:20)/=parnam(np)(:20)) then
             call prterr(lundia    ,'V096'    ,parrd(np)(:20)       )
             !
             error = .true.
             goto 9999
          endif
       enddo
       !
       ! Define INTERP
       !
       interp = 'linear    '
       if (disint(n)=='N') interp = 'block    '
       !
       ! Define format for number of values NRVAL
       !
       write (fmtdis(13)(8:10), '(i3)') nrval
       !
       ! Define table name
       !
       tablnm = 'T-serie DIS ' // namsrc(n) // ' for run: '
       !
       ! Write first 7 description records to file
       !
       write (lunout, fmtdis(2), rec = irec + 1) &
           & keywrd(1), quote, tablnm, runid, quote, eol
       write (lunout, fmtdis(3), rec = irec + 2) &
           & keywrd(2), quote, cntain(:10), quote, cntain(11:40), eol
       write (lunout, fmtdis(4), rec = irec + 3) &
           & keywrd(3), quote, namsrc(n), quote, eol
       write (lunout, fmtdis(7), rec = irec + 4) &
           & keywrd(8), quote, 'non-equidistant', quote, eol
       write (lunout, fmtdis(8), rec = irec + 5) keywrd(9), itdate, eol
       write (lunout, fmtdis(9), rec = irec + 6) &
           & keywrd(10), quote, 'minutes', quote, eol
       write (lunout, fmtdis(10), rec = irec + 7) &
           & keywrd(12), quote, interp, quote, eol
       !
       ! Write NPARA parameter names
       !
       do l = 1, npara
          write (lunout, fmtdis(11), rec = irec + 7 + l) &
              & keywrd(14), quote, parnam(l), quote, keywrd(15)(:10), quote, parunt(l), quote, eol
       enddo
       !
       ! Write number of time dependent data to file
       !
       write (lunout, fmtdis(12), rec = irec + 8 + npara) &
           & keywrd(16), ntimrd, eol
       !
       ! Re-define IREC 8+NPARA description records forward
       !
       irec = irec + 8 + npara
       !
       ! Read time dependent data for discharge location (N)
       ! for all times NTIMRD on file
       !
       do nr = 1, ntimrd
          !
          ! Read time.
          !
          ibeg = 1
          call read1r(record    ,lrecrd    ,ibeg      ,iend      ,timrd     , &
                    & rdfaul    ,ier       )
          !
          ! Premature EOR (IER = 0) or empty value (IER = -1)
          !
          if (ier<=0) then
             error = .true.
             call prterr(lundia    ,'G007'    ,fildis    )
             write (lundia, '(a,a)') 'RECORD: ', record(:72)
             goto 9999
          endif
          !
          ! Re-defined time read (TIMRD ) in minutes
          !
          timrd = timrd*timscl
          dtmin = dt * tunit / 60.0_fp
          !
          ! Perform some checks on time value read
          !
          call chckit(lundia    ,error     ,fildis    ,timrd     ,dtmin     , &
                    & ittdep    ,itold     ,itstrt    ,nr        ,gdp       )
          if (error) goto 9999
          !
          ! Read discharge value, concentrations for LSTSC and optional
          ! magnitude and direction
          !
          ibeg = iend + 1
          call readnr(record    ,lrecrd    ,ibeg      ,iend      ,nrval     , &
                    & rval      ,rdfaul    ,ier       )
          !
          ! Premature EOR or empty value (IER = -1)
          !
          if (ier<=0) then
             error = .true.
             call prterr(lundia    ,'G007'    ,fildis    )
             write (lundia, '(a,a)') 'RECORD: ', record(:72)
             goto 9999
          endif
          !
          ! If an NaN is read -> error
          !
          do j = 2, nrval
              if ( isnan(rval(j)) ) then
                 write(errmsg,'(a,a)') 'NaN in ', trim(fildis)
                 call prterr(lundia    ,'P004'    ,errmsg      )
                 error = .true.
                 goto 9999
              endif
          enddo
          !
          ! Check for negative values, when not a power station/culvert
          !
          if (disflg/='inoutlet  ' .and. disflg/='power     ' .and.             &
             & disflg/='culvert   ') then
             do j = 2, lstsc + 1
                if (rval(j)<0.0) then
                   errmsg = 'Concentration at discharge ' // namsrc(n)
                   call prterr(lundia    ,'V061'    ,errmsg    )
                   error = .true.
                   goto 9999
                endif
             enddo
          endif
          !
          ! Check for negative magnitude
          !
          if (nrval>lstsc + 1) then
             if (rval(lstsc + 2)<0.0) then
                errmsg = 'Flow magnitude at discharge ' // namsrc(n)
                call prterr(lundia    ,'V061'    ,errmsg    )
                error = .true.
                goto 9999
             endif
          endif
          !
          ! Write NRVAL time dependent data
          !
          write (lunout, fmtdis(13), rec = irec + nr) &
              & timrd, (rval(l), l = 1, nrval), eol
          !
          ! Re-define ITOLD and
          ! As long as NR < NTIMRD read new RECORD
          !
          itold = ittdep
          if (nr<=ntimrd - 1) read (lunrd, '(a)') record(:lrecrd)
       enddo
       !
       ! Re-define IREC NTIMRD records forward
       !
       irec = irec + ntimrd
       !
       ! Define maximum time
       !
       if (itold/= - 1) then
          if (itold < itfinish) then
             write(errmsg,'(a,a,a,a,a)') 'Last time in file ', trim(fildis), ' for discharge #',trim(namsrc(n)),'# <'
             call prterr(lundia    ,'U042'    ,errmsg)
             error = .true.
             goto 9999
          endif
       endif
    enddo
    !
    ! stop reading file
    !
 9999 continue
    deallocate(parrd)
end subroutine rdtddn

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


!
!  Copyright 2018 SALMON developers
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!
!-----------------------------------------------------------------------------------------
module builtin_pbe
  use salmon_math, only: erfc_salmon
  implicit none

  real(8),parameter :: Pi=3.141592653589793d0

contains

  subroutine exc_cor_pbe(nl, rho, grho_s, exc, eexc, vexc, nd, ifdx, ifdy, ifdz, nabx, naby, nabz)
    implicit none
    integer, intent(in) :: nl
    real(8), intent(in) :: rho(nl), grho_s(nl, 3)
    real(8), intent(out) :: exc(nl), eexc(nl), vexc(nl)
    integer, intent(in) :: nd
    integer, intent(in) :: ifdx(-nd:nd, 1:nl)
    integer, intent(in) :: ifdy(-nd:nd, 1:nl)
    integer, intent(in) :: ifdz(-nd:nd, 1:nl)
    real(8), intent(in) :: nabx(-nd:nd)
    real(8), intent(in) :: naby(-nd:nd)
    real(8), intent(in) :: nabz(-nd:nd)

    real(8) :: agrho_s(NL)
    integer :: i
    real(8) :: omega
    real(8) :: FxHSE,drho_FxHSE,ds_FxHSE,ec_unif,drs_ec_unif,H,drs_H,dt_H
    real(8) :: rho_t,rs,kf,ks,grho(NL,3),agrho(NL),srho,trho
    real(8) :: exlda,Ex(NL),Ec(NL),Vx(NL),Vc(NL)
    real(8) :: dFxdrho_tbl(NL,3),dHksrho_tbl(NL,3)

    ! call rho_j_tau(GS_RT,rho_s,tau_s,j_s,grho_s,lrho_s)

    agrho_s(:)=sqrt(grho_s(:,1)**2+grho_s(:,2)**2+grho_s(:,3)**2)
    grho(:,1)=2*grho_s(:,1)
    grho(:,2)=2*grho_s(:,2)
    grho(:,3)=2*grho_s(:,3)
    agrho(:)=2*agrho_s(:)
    omega=0d0

    do i=1,NL
      rho_t=rho(i)
      if(rho(i) < 1d-10) rho_t=rho(i)+1d-10
      rs=(3d0/(4*pi*rho_t))**(1d0/3d0)
      kf=(3*pi**2*rho_t)**(1d0/3d0)
      ks=sqrt(4*kf/pi)
      srho=agrho(i)/(2*kf*rho_t)
      trho=agrho(i)/(2*ks*rho_t)
      call HSEFx(rho_t,srho,omega,FxHSE,drho_FxHSE,ds_FxHSE)
      call PBEc_old(rho_t,trho,ec_unif,drs_ec_unif,H,drs_H,dt_H)
      exlda=-3d0/(4*pi)*kf
      Ex(i)=exlda*FxHSE
      Vx(i)=exlda*(4d0/3d0+rho(i)*drho_FxHSE-4*srho/3d0*ds_FxHSE)
      Ec(i)=ec_unif+H
      Vc(i)=ec_unif+H-rs/3*(drs_ec_unif+drs_H)-7d0/6d0*trho*dt_H
      dFxdrho_tbl(i,:)=ds_FxHSE*grho(i,:)/agrho(i)
      dHksrho_tbl(i,:)=dt_H/(2*ks)*grho(i,:)/agrho(i)
    enddo

    do i=1,NL
      Vx(i)=Vx(i)+3d0/(8*pi)*( &
   &    nabx(1)*(dFxdrho_tbl(ifdx(1,i),1)-dFxdrho_tbl(ifdx(-1,i),1)) &
   &   +nabx(2)*(dFxdrho_tbl(ifdx(2,i),1)-dFxdrho_tbl(ifdx(-2,i),1)) &
   &   +nabx(3)*(dFxdrho_tbl(ifdx(3,i),1)-dFxdrho_tbl(ifdx(-3,i),1)) &
   &   +nabx(4)*(dFxdrho_tbl(ifdx(4,i),1)-dFxdrho_tbl(ifdx(-4,i),1)) &
   &   +naby(1)*(dFxdrho_tbl(ifdy(1,i),2)-dFxdrho_tbl(ifdy(-1,i),2)) &
   &   +naby(2)*(dFxdrho_tbl(ifdy(2,i),2)-dFxdrho_tbl(ifdy(-2,i),2)) &
   &   +naby(3)*(dFxdrho_tbl(ifdy(3,i),2)-dFxdrho_tbl(ifdy(-3,i),2)) &
   &   +naby(4)*(dFxdrho_tbl(ifdy(4,i),2)-dFxdrho_tbl(ifdy(-4,i),2)) &
   &   +nabz(1)*(dFxdrho_tbl(ifdz(1,i),3)-dFxdrho_tbl(ifdz(-1,i),3)) &
   &   +nabz(2)*(dFxdrho_tbl(ifdz(2,i),3)-dFxdrho_tbl(ifdz(-2,i),3)) &
   &   +nabz(3)*(dFxdrho_tbl(ifdz(3,i),3)-dFxdrho_tbl(ifdz(-3,i),3)) &
   &   +nabz(4)*(dFxdrho_tbl(ifdz(4,i),3)-dFxdrho_tbl(ifdz(-4,i),3)) )
      Vc(i)=Vc(i)-( &
   &    nabx(1)*(dHksrho_tbl(ifdx(1,i),1)-dHksrho_tbl(ifdx(-1,i),1)) &
   &   +nabx(2)*(dHksrho_tbl(ifdx(2,i),1)-dHksrho_tbl(ifdx(-2,i),1)) &
   &   +nabx(3)*(dHksrho_tbl(ifdx(3,i),1)-dHksrho_tbl(ifdx(-3,i),1)) &
   &   +nabx(4)*(dHksrho_tbl(ifdx(4,i),1)-dHksrho_tbl(ifdx(-4,i),1)) &
   &   +naby(1)*(dHksrho_tbl(ifdy(1,i),2)-dHksrho_tbl(ifdy(-1,i),2)) &
   &   +naby(2)*(dHksrho_tbl(ifdy(2,i),2)-dHksrho_tbl(ifdy(-2,i),2)) &
   &   +naby(3)*(dHksrho_tbl(ifdy(3,i),2)-dHksrho_tbl(ifdy(-3,i),2)) &
   &   +naby(4)*(dHksrho_tbl(ifdy(4,i),2)-dHksrho_tbl(ifdy(-4,i),2)) &
   &   +nabz(1)*(dHksrho_tbl(ifdz(1,i),3)-dHksrho_tbl(ifdz(-1,i),3)) &
   &   +nabz(2)*(dHksrho_tbl(ifdz(2,i),3)-dHksrho_tbl(ifdz(-2,i),3)) &
   &   +nabz(3)*(dHksrho_tbl(ifdz(3,i),3)-dHksrho_tbl(ifdz(-3,i),3)) &
   &   +nabz(4)*(dHksrho_tbl(ifdz(4,i),3)-dHksrho_tbl(ifdz(-4,i),3)) )
    enddo
    Exc=(Ex+Ec)
    Eexc=rho*(Ex+Ec)
    Vexc=Vx+Vc
    return
  End Subroutine Exc_Cor_PBE

  SUBROUTINE HSEFx(rho,ss,omega,FxHSE,drho_FxHSE,ds_FxHSE)
  implicit none

  real(8),intent(IN) :: rho,ss,omega
  real(8),intent(OUT) :: FxHSE,drho_FxHSE,ds_FxHSE

  !real(8),parameter :: pi=3.1415926535897932d0
  real(8),parameter :: sqpi=1.7724538509055160d0
  real(8),parameter :: eps=1d-10
  real(8),parameter :: A=1.0161144D0,B=-3.7170836D-1
  real(8),parameter :: C=-7.7215461D-2,D=5.7786348D-1
  real(8),parameter :: E=-5.1955731D-2
  real(8),parameter :: Ha1=9.79681d-3,Ha2=4.10834d-2
  real(8),parameter :: Ha3=1.87440d-1,Ha4=1.20824d-3
  real(8),parameter :: Ha5=3.47188d-2
  real(8),parameter :: Fc1=6.4753871d0,Fc2=4.7965830d-1
  real(8),parameter :: ea1=-1.128223946706117d0,ea2=1.452736265762971d0
  real(8),parameter :: ea3=-1.243162299390327d0,ea4=0.971824836115601d0
  real(8),parameter :: ea5=-0.568861079687373d0,ea6=0.246880514820192d0
  real(8),parameter :: ea7=-0.065032363850763d0,ea8=0.008401793031216d0
  real(8),parameter :: wcut=14d0
  real(8),parameter :: EGscut=0.08d0,EGa1=-2.628417880d-2
  real(8),parameter :: EGa2=-7.117647788d-2,EGa3=8.534541323d-2
  real(8),parameter :: expcut=700d0
  real(8),parameter :: exei1=4.03640D0,exei2=1.15198D0
  real(8),parameter :: exei3=5.03627D0,exei4=4.19160D0
  real(8),parameter :: smax=8.572844D0,strans=8.3d0
  real(8),parameter :: sconst=1.879622316D1

  real(8) s,s2,s3,s4,s5,s6,kf,w,w2,w3,w4,w5,w6,w7,w8
  real(8) H,ds_H,ds_Hs,F,ds_CFs,eb1,Hsbw,Hsbw2,Hsbw3,Hsbw4
  real(8) DHsb,DHsb2,DHsb3,DHsb4,DHsb5,DHsb12,DHsb32,DHsb52,DHsb72,DHsb92
  real(8) A94,A942,A943,A944,HA94,HA942,HA943,HA945
  real(8) DHs,DHs2,DHs3,DHs4,DHs52,DHs72
  real(8) DHsw,DHsw12,DHsw32,DHsw52,DHsw72
  real(8) wsqDHsw,wsqDHsw3,wsqDHsw5,wsqDHsw7
  real(8) G_a,G_b,EG,ds_EGs,exer,exei
!  real(8) ei
!  real(8) derfc

  s=ss
! s is modified to enforce Lieb-Oxford bound
  if(s.gt.strans) then
    s=smax-sconst/s**2
  endif

  s2=s**2; s3=s**3; s4=s**4; s5=s**5; s6=s**6

  kf=(3d0*pi**2*rho)**(1d0/3d0)
  w=omega/kf

  w2=w**2; w3=w**3; w4=w**4; w5=w**5; w6=w**6; w7=w**7; w8=w**8

  H=(Ha1*s2+Ha2*s4)/(1d0+Ha3*s4+Ha4*s5+Ha5*s6)
  ds_H=(2*Ha1*s+4*Ha2*s3)/(1+Ha3*s4+Ha4*s5+Ha5*s6) &
 &    -(Ha1*s2+Ha2*s4)/(1+Ha3*s4+Ha4*s5+Ha5*s6)**2 &
 &    *(4*Ha3*s3+5*Ha4*s4+6*Ha5*s5)
  ds_Hs=ds_H*s2+2*s*H

  F=Fc1*H+Fc2
  ds_CFs=C*(2*s*F+s2*Fc1*ds_H)

  if(w.lt.wcut) then
    eb1=1.455915450052607d0
  else
    eb1=2.0d0
  endif

  Hsbw=s2*H+eb1*w2
  Hsbw2=Hsbw**2; Hsbw3=Hsbw**3; Hsbw4=Hsbw**4

  DHsb=D+s2*H+eb1*w2
  DHsb2=DHsb**2; DHsb3=DHsb**3; DHsb4=DHsb**4; DHsb5=DHsb**5
  DHsb12=DHsb**0.5d0; DHsb32=DHsb**1.5d0; DHsb52=DHsb**2.5d0
  DHsb72=DHsb**3.5d0; DHsb92=DHsb**4.5d0

  A94=9d0/(4d0*A)
  A942=A94**2; A943=A94**3; A944=A94**4

  HA94=A94*Hsbw
  HA942=HA94**2; HA943=HA94**3; HA945=HA94**5

  DHs=D+s2*H
  DHs2=DHs**2; DHs3=DHs**3; DHs4=DHs**4
  DHs52=DHs**2.5d0; DHs72=DHs**3.5d0

  DHsw=DHs+w2
  DHsw12=DHsw**0.5d0; DHsw32=DHsw**1.5d0; DHsw52=DHsw**2.5d0; DHsw72=DHsw**3.5d0

  wsqDHsw=w/DHsw12
  wsqDHsw3=wsqDHsw**3; wsqDHsw5=wsqDHsw**5; wsqDHsw7=wsqDHsw**7

  if(s.gt.EGscut) then
    G_a=sqpi*(15d0*E+6d0*C*(1d0+F*s2)*DHs+4d0*B*DHs2+8d0*A*DHs3) &
 &     /(16d0*DHs72) &
 &     -0.75d0*pi*sqrt(A)*exp(A94*H*s2) &
 &     *erfc_salmon(1.5d0*s*sqrt(H/A))
    G_b=(15d0/16d0*sqpi*s2)/DHs72
    EG=-(0.75d0*pi+G_a)/G_b
    ds_EGs=-14d0/5d0*sqpi*DHs52*ds_Hs &
 &         +sqpi*sqrt(A)/5*DHs52  &
 &         *exp(9*H*s2/(4*A))*erfc_salmon(1.5d0*s*sqrt(H/A)) &
 &         *(14+9/A*DHs)*ds_Hs &
 &         -6d0/5d0*DHs72/sqrt(H*s2)*ds_Hs &
 &         -2d0/5d0*DHs*ds_CFs-2d0/5d0*C*(1+F*s2)*ds_Hs &
 &         -8d0/15d0*B*DHs*ds_Hs-8d0/5d0*A*DHs2*ds_Hs
  else
    EG=EGa1+EGa2*s2+EGa3*s4
    ds_EGs=2*s*EG+s2*(2*EGa2*s+4*EGa3*s3)
  endif

  if(HA94 > eps) then
    if(HA94 < expcut) then
      exer=pi*exp(HA94)*erfc_salmon(sqrt(HA94))
      exei=exp(HA94)*ei(-HA94)
    else
      exer=pi*(1d0/(sqpi*sqrt(HA94))-1d0/(2d0*sqrt(pi*HA943)) &
 &    +3d0/(4d0*sqrt(pi*HA945)))
      exei=-(1d0/HA94) &
 &        *(HA942+exei1*HA94+exei2)/(HA942+exei3*HA94+exei4)
    endif
  endif

!  if(w < eps.and.s < eps) then
  if(HA94 < eps) then
    FxHSE=1d0
    drho_FxHSE=0d0
    ds_FxHSE=0d0
    return
  endif

  if(w < eps) then
    FxHSE=-8d0/9d0*(-A/2*(log(DHs)-log(H*s2))-A/2*exei &
 &                  +B/2/DHs+C*(1+F*s2)/2/DHs**2+(E+EG*s2)/DHs3)
    drho_FxHSE=0d0
    ds_FxHSE=-8d0/9d0*(ds_CFs/2/DHs2+ds_EGs/DHs3 &
 &                    +ds_Hs*(A/2*(-1/DHs-A94*exei)-B/2/DHs2-C*(1+F*s2)/DHs3-3*(E+EG*s2)/DHs4))
    return
  endif

  FxHSE=-8d0/9d0*( &
 & -A/2*(log(DHsb)-log(Hsbw))-A/2*exei &
 & +B/2/DHs*(1-wsqDHsw) &
 & +3d0/4d0*C*(1+F*s2)/DHs2*(2d0/3d0-wsqDHsw+wsqDHsw3/3) &
 & +15d0/8d0*(E+EG*s2)/DHs3*(8d0/15d0-wsqDHsw+2d0/3d0*wsqDHsw3-wsqDHsw5/5) &
 & )

  drho_FxHSE= -8d0/9d0*(2*w/3d0/sqpi/rho)*( &
 &   -3d0/4d0*sqrt(A)*pi*exp(A94*(H*s2+w2))*erfc_salmon(sqrt(A94*(H*s2+w2))) &
 &   +sqpi*(A/2/DHsw12+B/4/DHsw32+C*(1+F*s2)*3d0/8d0/DHsw52+(E+EG*s2)*15d0/16d0/DHsw72) )

  ds_FxHSE= -8d0/9d0*( &
 &  3d0/4d0*ds_CFs/DHs2*(2d0/3d0-wsqDHsw+wsqDHsw3/3) &
 & +15d0/8d0*ds_EGs/DHs3*(8d0/15d0-wsqDHsw+2d0/3d0*wsqDHsw3-wsqDHsw5/5) &
 & -ds_Hs*( &
 &         9d0/8d0*exei+A/2/DHsb &
 &        +3d0/4d0*B/DHs2*(2d0/3d0-wsqDHsw+wsqDHsw3/3) &
 &        +15d0/8d0*C*(1+F*s2)/DHs3*(8d0/15d0-wsqDHsw+2d0/3d0*wsqDHsw3-wsqDHsw5/5) &
 &        +105d0/16d0*(E+EG*s2)/DHs4*(16d0/35d0-wsqDHsw+wsqDHsw3 &
 &                                    -3d0/5d0*wsqDHsw5+wsqDHsw7/7) &
 &          ))

  if(w > wcut) return

  FxHSE=FxHSE-8d0/9d0*( &
 &   -9d0/8d0*exei*(-ea2*w2+A94*ea4*w4-A942*ea6*w6+A943*ea8*w8) &
 &   -9d0/8d0*(ea4*w4/Hsbw+ea6*w6*(-A94/Hsbw+1/Hsbw2)+ea8*w8*(A942/Hsbw-A94/Hsbw2+2/Hsbw3)) &
 &   -9d0/8d0*exer/sqrt(A94)*(ea1*w-ea3*w3*A94+ea5*w5*A942-ea7*w7*A943) &
 &   -9d0/8d0*sqpi/sqrt(Hsbw)*(ea3*w3+ea5*w5*(-A94+1/(2*Hsbw)) &
 &                 +ea7*w7*(A942-A94/Hsbw/2+3/(4*Hsbw2))) &
 &   +A/2*(ea2*w2/DHsb+ea4*w4/DHsb2+ea6*w6*2/DHsb3+ea8*w8*6/DHsb4) &
 &   +A*sqpi*(ea1*w/2/DHsb12+ea3*w3/4/DHsb32+ea5*w5*3d0/8d0/DHsb52+ea7*w7*15d0/16d0/DHsb72) )

  ds_FxHSE=ds_FxHSE-8d0/9d0*( &
 &   -ds_Hs*(-9d0/8d0*((ea2*w2*A94-ea4*w4*A942+ea6*w6*A943-ea8*w8*A944)*exei &
 &                     +ea2*w2/Hsbw+ea4*w4*(-A94/Hsbw+1/Hsbw2) &
 &                     +ea6*w6*(A942/Hsbw-A94/Hsbw2+2/Hsbw3) &
 &                     +ea8*w8*(-A943/Hsbw+A942/Hsbw2-2*A94/Hsbw3+6/Hsbw4)) &
 &           -9d0/8d0*(sqrt(A94)*(-ea1*w+ea3*w3*A94-ea5*w5*A942+ea7*w7*A943)*exer &
 &                    +sqpi/sqrt(Hsbw)*(ea1*w+ea3*w3*(-A94+1/Hsbw/2) &
 &                                     +ea5*w5*(A942-A94/2/Hsbw+3d0/4d0/Hsbw2) &
 &                                     +ea7*w7*(-A943+A942/2/Hsbw &
 &                                     -3d0/4d0*A94/Hsbw2+15d0/8d0/Hsbw3))) &
 &           +A/2*(ea2*w2/DHsb2+ea4*w4*2/DHsb3+ea6*w6*6/DHsb4+ea8*w8*24/DHsb5) &
 &           +A*sqpi*(ea1*w/4/DHsb32+ea3*w3*3d0/8d0/DHsb52 &
 &                   +ea5*w5*15d0/16d0/DHsb72+ea7*w7*105d0/32d0/DHsb92) ))
  end subroutine

  SUBROUTINE PBEc_old(rho,t,ec_unif,drs_ec_unif,H,drs_H,dt_H)
  implicit none
  real(8),intent(IN) :: rho,t
  real(8),intent(OUT) :: ec_unif,drs_ec_unif,H,drs_H,dt_H
  !real(8),parameter :: pi=3.141592653589793d0
  integer,parameter :: p=1
  real(8),parameter :: Au=0.031091d0,a1=0.21370d0,b1=7.5957d0
  real(8),parameter :: b2=3.5876d0,b3=1.6382d0,b4=0.49294d0
  real(8),parameter :: beta=0.066725d0,gamma=0.031091d0
  real(8) rs,sqrs,kf,ks,b14,exec,Ah,Ah21,Ah421,dec_Ah,dA_H

  rs=(3d0/(4d0*pi*rho))**(1d0/3d0)
  sqrs=sqrt(rs)
  kf=(3*pi**2*rho)**(1d0/3d0)
  ks=(4*kf/pi)**0.5d0
  b14=b1*sqrs+b2*rs+b3*rs*sqrs+b4*rs**(p+1)
  ec_unif=-2*Au*(1+a1*rs)*log(1+1/(2*Au*b14))
  exec=exp(-ec_unif/gamma)
  Ah=beta/gamma/(exec-1)
  Ah21=1+Ah*t**2
  Ah421=1+Ah*t**2+Ah**2*t**4
  H=gamma*log(1+beta/gamma*t**2*Ah21/Ah421)

  drs_ec_unif=-2*Au*a1*log(1+1/(2*Au*b14)) &
 &            +(1+a1*rs)*(b1/(2*sqrs)+b2+1.5d0*b3*sqrs+(p+1)*b4*rs**p)/b14**2 &
 &            /(1+1/(2*Au*b14))

  dec_Ah=beta/gamma**2*exec/(exec-1)**2
  dA_H=gamma*(beta/gamma)*t**2*(t**2/Ah421-Ah21*(t**2+2*Ah*t**4)/Ah421**2) &
 &     /(1+beta/gamma*t**2*Ah21/Ah421)
  drs_H=drs_ec_unif*dec_Ah*dA_H

  dt_H=gamma*((beta/gamma)*2*t*Ah21/Ah421 &
 &            +beta/gamma*t**2*(2*Ah*t*Ah421-Ah21*(2*Ah*t+4*Ah**2*t**3))/Ah421**2) &
 &    /(1+beta/gamma*t**2*Ah21/Ah421)
  return
  end subroutine

  FUNCTION ei(x)
  implicit none
  real(8) ei,x,xx
  integer,parameter :: one=1
  if(x.ge.0d0) stop 'bad x in ei'
  xx=-x
  ei=-expint(one,xx)
  return
  end function

  FUNCTION expint(n,x)
  implicit none
  INTEGER n,MAXIT
  REAL(8) expint,x,EPS,FPMIN,EULER
  PARAMETER (MAXIT=100,EPS=1.d-7,FPMIN=1.d-30,EULER=.5772156649d0)
  INTEGER i,ii,nm1
  REAL(8) a,b,c,d,del,fact,h,psi
  nm1=n-1
  if(n.lt.0.or.x.lt.0.d0.or.(x.eq.0.d0.and.(n.eq.0.or.n.eq.1))) then
    stop 'bad arguments in expint'
  else if(n.eq.0)then
    expint=exp(-x)/x
  else if(x.eq.0.d0)then
    expint=1.d0/nm1
  else if(x.gt.1.d0)then
    b=x+n
    c=1.d0/FPMIN
    d=1.d0/b
    h=d
    do i=1,MAXIT
      a=-i*(nm1+i)
      b=b+2.d0
      d=1.d0/(a*d+b)
      c=b+a/c
      del=c*d
      h=h*del
      if(abs(del-1.d0).lt.EPS)then
        expint=h*exp(-x)
        return
      endif
    enddo
    stop 'continued fraction failed in expint'
  else
    if(nm1.ne.0)then
      expint=1.d0/nm1
    else
      expint=-log(x)-EULER
    endif
    fact=1.d0
    do i=1,MAXIT
      fact=-fact*x/i
      if(i.ne.nm1) then
        del=-fact/(i-nm1)
      else
        psi=-EULER
        do ii=1,nm1
          psi=psi+1.d0/ii
        enddo
        del=fact*(-log(x)+psi)
      endif
      expint=expint+del
      if(abs(del).lt.abs(expint)*EPS) return
    enddo
    stop 'series failed in expint'
  endif
  return
  END function

  ! Perdew-Burke-Ernzerhof correlation energy
    SUBROUTINE PBEc(rhoa,rhob,grhoa,grhob,ec,dec_drhoa,dec_drhob,dec_dgrhoa,dec_dgrhob)
    implicit none
    real(8),parameter :: Pi=3.141592653589793d0
    real(8),parameter :: beta=0.066725d0,gamma=0.031091d0
    real(8),parameter :: eps=1d-12
    real(8) :: rhoa,rhob,grhoa(3),grhob(3)
    real(8) :: ec,dec_drhoa,dec_drhob,dec_dgrhoa(3),dec_dgrhob(3)
    real(8) :: ec_PW,dec_drhoa_PW,dec_drhob_PW
    real(8) :: rho,grho(3),rs,kf,ks,agrho,zeta,phi,t,exp_ec_PW,A,At2,At21,At421,bgt2A,H
    real(8) :: dH_dphi,dH_dt,dH_dA,dphi_dzeta,dzeta_drhoa,dzeta_drhob,dphi_drhoa,dphi_drhob
    real(8) :: dks_drhoa,dks_drhob,dA_drhoa,dA_drhob,dt_drhoa,dt_drhob
    real(8) :: dH_drhoa,dH_drhob,dH_dgrhoa(3),dH_dgrhob(3)

    call PWc(rhoa,rhob,ec_PW,dec_drhoa_PW,dec_drhob_PW)

    rho=rhoa+rhob
    grho(:)=grhoa(:)+grhob(:)
    rs=(3d0/(4d0*Pi*rho))**(1d0/3d0)
    kf=(3*Pi**2*rho)**(1d0/3d0)
    ks=sqrt(4*kf/Pi)
    agrho=sqrt(grho(1)**2+grho(2)**2+grho(3)**2)
    zeta=(rhoa-rhob)/rho
    if(abs(zeta-1d0) < eps) zeta=1d0-eps
    if(abs(zeta+1d0) < eps) zeta=-1d0+eps
    phi=((1+zeta)**(2d0/3d0)+(1-zeta)**(2d0/3d0))/2d0
    t=agrho/(2*phi*ks*rho)

    exp_ec_PW=exp(-ec_PW/gamma/phi**3)
    A=beta/gamma/(exp_ec_PW-1)
    At2=A*t*t
    At21=1+At2
    At421=At21+At2**2
    bgt2A=1+beta/gamma*t**2*At21/At421
    H=gamma*phi**3*log(bgt2A)

    dH_dphi=3*H/phi
    dH_dt=gamma*phi**3*beta/gamma*2*t*(1+2*At2)/At421**2/bgt2A
    dH_dA=gamma*phi**3*(-beta/gamma)*t**4*At2*(2+At2)/At421**2/bgt2A
    dphi_dzeta=(1/((1+zeta)**2+eps)**(1d0/6d0)-1/((1-zeta)**2+eps)**(1d0/6d0))/3
    dzeta_drhoa=2*rhob/rho**2
    dzeta_drhob=-2*rhoa/rho**2
    dphi_drhoa=dzeta_drhoa*dphi_dzeta
    dphi_drhob=dzeta_drhob*dphi_dzeta

    dks_drhoa=ks/(6*rho)
    dks_drhob=ks/(6*rho)
    dA_drhoa=A**2/beta*exp_ec_PW*(dec_drhoa_PW/phi**3-3*ec_PW/phi**4*dphi_drhoa)
    dA_drhob=A**2/beta*exp_ec_PW*(dec_drhob_PW/phi**3-3*ec_PW/phi**4*dphi_drhob)
    dt_drhoa=-t/phi*dphi_drhoa-t/ks*dks_drhoa-t/rho
    dt_drhob=-t/phi*dphi_drhob-t/ks*dks_drhob-t/rho
    dH_drhoa=dphi_drhoa*dH_dphi+dt_drhoa*dH_dt+dA_drhoa*dH_dA
    dH_drhob=dphi_drhob*dH_dphi+dt_drhob*dH_dt+dA_drhob*dH_dA

    dH_dgrhoa(:)=dH_dt*t*grho(:)/agrho**2
    dH_dgrhob(:)=dH_dgrhoa(:)

    ec=ec_PW+H
    dec_drhoa=dec_drhoa_PW+dH_drhoa
    dec_drhob=dec_drhob_PW+dH_drhob
    dec_dgrhoa(:)=dH_dgrhoa(:)
    dec_dgrhob(:)=dH_dgrhob(:)

    return
    end subroutine

    ! Perdew-Wang correlation energy

      SUBROUTINE PWc(rhoa,rhob,ec,dec_drhoa,dec_drhob)
      implicit none
      real(8),parameter :: Pi=3.141592653589793d0
      real(8),intent(IN) :: rhoa,rhob
      real(8),intent(OUT) :: ec,dec_drhoa,dec_drhob
      real(8),parameter :: pec0=1d0,pec1=1d0,par=1d0
      real(8),parameter :: Aec0=0.031091d0,Aec1=0.015545d0,Aar=0.016887d0
      real(8),parameter :: a1ec0=0.21370d0,a1ec1=0.20548d0,a1ar=0.11125d0
      real(8),parameter :: b1ec0=7.5957d0,b1ec1=14.1189d0,b1ar=10.357d0
      real(8),parameter :: b2ec0=3.5876d0,b2ec1=6.1977d0,b2ar=3.6231d0
      real(8),parameter :: b3ec0=1.6382d0,b3ec1=3.3662d0,b3ar=0.88026d0
      real(8),parameter :: b4ec0=0.49294d0,b4ec1=0.62517d0,b4ar=0.49671d0
      real(8),parameter :: f2d=1.709921d0,eps=1d-15
      real(8) rs,sqrs,zeta,fzeta,dfzeta,Q0ec0,Q0ec1,Q0ar,Q1ec0,Q1ec1,Q1ar
      real(8) Qdec0,Qdec1,Qdar,Gec0,Gec1,Gar,dGec0,dGec1,dGar,decdrs,decdzeta

      rs=(3d0/(4d0*Pi*(rhoa+rhob)))**(1d0/3d0)
      sqrs=sqrt(rs)
      zeta=(rhoa-rhob)/(rhoa+rhob)
      if(abs(zeta-1d0) < eps) zeta=1d0-eps
      if(abs(zeta+1d0) < eps) zeta=-1d0+eps
      fzeta=((1+zeta)**(4d0/3d0)+(1-zeta)**(4d0/3d0)-2)/(2d0**(4d0/3d0)-2d0)
      dfzeta=4d0/3d0*((1+zeta)**(1d0/3d0)-(1-zeta)**(1d0/3d0))/(2d0**(4d0/3d0)-2d0)

      Q0ec0=-2*Aec0*(1+a1ec0*rs)
      Q0ec1=-2*Aec1*(1+a1ec1*rs)
      Q0ar =-2*Aar *(1+a1ar *rs)
      Q1ec0=2*Aec0*(b1ec0*sqrs+b2ec0*rs+b3ec0*rs*sqrs+b4ec0*rs**(pec0+1))
      Q1ec1=2*Aec1*(b1ec1*sqrs+b2ec1*rs+b3ec1*rs*sqrs+b4ec1*rs**(pec1+1))
      Q1ar =2*Aar *(b1ar *sqrs+b2ar *rs+b3ar *rs*sqrs+b4ar *rs**(par +1))
      Qdec0=Aec0*(b1ec0/sqrs+2*b2ec0+3*b3ec0*sqrs+2*(pec0+1)*b4ec0*rs**pec0)
      Qdec1=Aec1*(b1ec1/sqrs+2*b2ec1+3*b3ec1*sqrs+2*(pec1+1)*b4ec1*rs**pec1)
      Qdar =Aar *(b1ar /sqrs+2*b2ar +3*b3ar *sqrs+2*(par +1)*b4ar *rs**par )
      Gec0=Q0ec0*log(1+1d0/Q1ec0)
      Gec1=Q0ec1*log(1+1d0/Q1ec1)
      Gar =Q0ar *log(1+1d0/Q1ar )
      dGec0=-2*Aec0*a1ec0*log(1+1d0/Q1ec0)-Q0ec0*Qdec0/(Q1ec0**2+Q1ec0)
      dGec1=-2*Aec1*a1ec1*log(1+1d0/Q1ec1)-Q0ec1*Qdec1/(Q1ec1**2+Q1ec1)
      dGar =-2*Aar *a1ar *log(1+1d0/Q1ar )-Q0ar *Qdar /(Q1ar **2+Q1ar )

      ec=Gec0-Gar*fzeta/f2d*(1-zeta**4)+(Gec1-Gec0)*fzeta*zeta**4
      decdrs=dGec0*(1-fzeta*zeta**4)+dGec1*fzeta*zeta**4-dGar*fzeta/f2d*(1-zeta**4)
      decdzeta=dfzeta*(-Gar/f2d*(1-zeta**4)+(Gec1-Gec0)*zeta**4) &
    &         +4*zeta**3*(fzeta*(Gec1-Gec0)+Gar*fzeta/f2d)

      dec_drhoa=(-rs/3*decdrs-(zeta-1)*decdzeta)/(rhoa+rhob)
      dec_drhob=(-rs/3*decdrs-(zeta+1)*decdzeta)/(rhoa+rhob)

      return
      end subroutine

end module builtin_pbe


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 05.02/06                              *
C  *    TEST TITLE : Deferral mode BNIL                    *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /GLOBCH/ PIDENT,    GLBERR,    TSTMSG,     FUNCID,
     1                DUMCH
      CHARACTER       PIDENT*40, GLBERR*60, TSTMSG*900, FUNCID*80,
     1                DUMCH(20)*20

C deferral mode
      INTEGER    PASAP,     PBNIG,     PBNIL,     PASTI,     PWAITD
      PARAMETER (PASAP = 0, PBNIG = 1, PBNIL = 2, PASTI = 3, PWAITD = 4)
C modification mode
      INTEGER    PNIVE,     PUWOR,     PUQUM
      PARAMETER (PNIVE = 0, PUWOR = 1, PUQUM = 2)
C dynamic modification
      INTEGER    PIRG,      PIMM,      PCBS
      PARAMETER (PIRG = 0,  PIMM = 1,  PCBS = 2)
C state of visual representation
C                correct    deferred    simulated
      INTEGER    PVROK,     PVRDFR,     PVRSIM
      PARAMETER (PVROK = 0, PVRDFR = 1, PVRSIM = 2)
C clear control flag
C                conditionally  always
      INTEGER    PCONDI,        PALWAY
      PARAMETER (PCONDI = 0,    PALWAY = 1)
C operating mode
C                request  sample    event
      INTEGER    PREQU,   PSAMPL,   PEVENT
      PARAMETER (PREQU=0, PSAMPL=1, PEVENT=2)
C echo switch
      INTEGER    PNECHO,   PECHO
      PARAMETER (PNECHO=0, PECHO=1)

C  Declare program-specific variables
      INTEGER    MODMD, SEVRTY, OLDSVR, INMODE, INCLAS, IND
      INTEGER    ITRIM, PCC, CT, INSTAT, INDEV(6)
      INTEGER    SPECWT, SPECON, ACTSEV, ACTPCC

      INTEGER    UNDWAY,   INACT
      PARAMETER (UNDWAY=0, INACT=1)

      LOGICAL    SVROK

      CHARACTER  PCCHAR*2, SVRLST*3
      CHARACTER  PCCLBL*34, PCCTAB(17)*34
      CHARACTER  MODLBL*4,  MODTAB(PNIVE:PUQUM)*4
      CHARACTER  SEVLBL*3,  SEVTAB(PIRG:PCBS)*3, SEVCHR*2
      CHARACTER  OLDLBL*10, OLDTAB(PVROK:PVRSIM)*10, EXPLBL*40
      CHARACTER             SVRCOD(PVROK:PVRSIM)*1
      CHARACTER  CNDLBL*200

      DATA MODTAB / 'NIVE', 'UWOR', 'UQUM' /
      DATA SEVTAB / 'IRG',  'IMM',  'CBS' /
      DATA OLDTAB / 'CORRECT', 'DEFERRED', 'SIMULATED' /
      DATA SVRCOD / 'c', 'd', 's' /

      CALL INITGL ('05.02/06')

C open PHIGS and workstation
      CALL XPOPPH (ERRFIL, MEMUN)
      CALL POPWK (WKID, CONID, WTYPE)
      CALL PQWKC (WKID, ERRIND, SPECON, SPECWT)
      CALL CHKINQ ('pqwkc', ERRIND)

C  picture change categories are numbered from 1 to 17:
C    workstation changes :  1 to 12
C    structure changes   : 13 to 17
C  set pcctab to names of picture change categories
      CALL SETPCL (PCCTAB)

C  invoke PICENV to set up standard picture environment
      CALL PICENV (WKID)

C  <set display update state> BNIL, NIVE
      CALL PSDUS (WKID, PBNIL, PNIVE)

C  <inquire number of available logical input devices> for primary
C     workstation
C  inclas = input class for which a device is available
      CALL PQLI (SPECWT, ERRIND, INDEV(1), INDEV(2), INDEV(3),
     1                           INDEV(4), INDEV(5), INDEV(6))
      CALL CHKINQ ('pqli', ERRIND)
      DO 100 IND = 1, 6
         IF (INDEV(IND) .GT. 0) THEN
            INCLAS = IND
            GOTO 110
         ENDIF
100   CONTINUE
      CALL UNMSG ('No input devices available for this workstation; ' //
     1            'test requires OUTIN workstation.')
110   CONTINUE

      CT = 0
C  for pcc = 1 to 17 = picture change category
      DO 210 PCC = 1,17
         PCCLBL = PCCTAB(PCC)
         WRITE (PCCHAR, '(I2)') PCC
C        for instat = underway, inactive
         DO 200 INSTAT = UNDWAY, INACT
C           <redraw all structures> to make svr = CORRECT
            CALL PRST (WKID, PCONDI)
            IF (INSTAT .EQ. UNDWAY) THEN
C              activate input:
               CT = CT + 1
               IF (MOD (CT,2) .EQ. 0) THEN
                 INMODE = PSAMPL
               ELSE
                 INMODE = PEVENT
               ENDIF

C              vary order of operation:
               IF (MOD(CT, 3) .EQ. 0) THEN
C                 invoke PICCHG for category = pcc
                  CALL PICCHG (WKID, '0,1,2', PCCHAR, ACTSEV, ACTPCC)
C                 use DOIMOD to set mode of input device to inmode
                  CALL DOIMOD (WKID, INCLAS, INMODE)
               ELSE
                  CALL DOIMOD (WKID, INCLAS, INMODE)
                  CALL PICCHG (WKID, '0,1,2', PCCHAR, ACTSEV, ACTPCC)
               ENDIF

               IF (ACTPCC .NE. PCC) THEN
                  CALL INMSG ('Cannot do picture change, '         //
     1                 '[BNIL, NIVE] test with input underway '    //
     2                 'skipped for category = '                   //
     3                 PCCLBL (1:ITRIM(PCCLBL)) // '.')
                  GOTO 180
               ENDIF

               CALL SETMSG ('3 4 13 30 31', 'A picture change of '    //
     1                      'category ' // PCCLBL (1:ITRIM(PCCLBL))   //
     2                      ' performed under mode [BNIL, NIVE] '     //
     3                      'when input is underway should leave the '//
     4                      'state of visual representation CORRECT.')
               CALL IFPF (SVROK(WKID, 'c'))
C  skip_pcc_act:
180            CONTINUE
C              deactivate local input:
               CALL DOIMOD (WKID, INCLAS, PREQU)
            ELSE
C input is inactive
C              invoke PICCHG for category = pcc
               CALL PICCHG (WKID, '0,1,2', PCCHAR, ACTSEV, ACTPCC)

               IF (ACTPCC .NE. PCC) THEN
                  CALL INMSG ('Cannot do picture change, '         //
     1                 '[BNIL, NIVE] test with input inactive '    //
     2                 'skipped for category = '                   //
     3                 PCCLBL (1:ITRIM(PCCLBL)) // '.')
                  GOTO 185
               ENDIF

               CALL SETMSG ('3 4 13 15 31', 'A picture change of '  //
     1                      'category ' // PCCLBL (1:ITRIM(PCCLBL)) //
     2                      ' performed under mode [BNIL, NIVE] '   //
     3                      'when input is inactive should leave '  //
     4                      'the state of visual representation '   //
     5                      'DEFERRED or CORRECT.')
               CALL IFPF (SVROK(WKID, 'dc'))
C  skip_pcc_inact:
185            CONTINUE
            ENDIF
C  next instat
200      CONTINUE
C  next pcc
210   CONTINUE

      CT = 0
      DO 330 MODMD = PUWOR, PUQUM
         MODLBL = MODTAB(MODMD)
      DO 340 SEVRTY = PIRG, PCBS
         SEVLBL = SEVTAB(SEVRTY)
         WRITE (SEVCHR, '(I2)') SEVRTY
      DO 350 OLDSVR = PVROK, PVRSIM
         OLDLBL = OLDTAB (OLDSVR)
         CNDLBL = 'update mode = [BNIL, ' // MODLBL // '], '     //
     1            'picture change severity = ' // SEVLBL // ', ' //
     2            'old state of visual representation = '        //
     3            OLDLBL
      DO 360 INSTAT = UNDWAY, INACT
C        attempt to set SVR to oldsvr
         CALL SETSVR (WKID, OLDSVR)

C        <set display update state> BNIL, modmd
         CALL PSDUS (WKID, PBNIL, MODMD)
C        if (current_SVR not= oldsvr) goto skip_oldsvr
         IF (.NOT. SVROK (WKID, SVRCOD(OLDSVR))) THEN
            CALL INMSG ('Cannot set old state of visual '           //
     1           'representation = ' // OLDLBL (1:ITRIM(OLDLBL))    //
     2           '; skipping test for ' // CNDLBL(1:ITRIM(CNDLBL))  //
     3           '.')
            GOTO 355
         ENDIF

         IF (INSTAT .EQ. UNDWAY) THEN
            CT = CT + 1
            IF (MOD(CT, 2) .EQ. 0) THEN
               INMODE = PSAMPL
            ELSE
               INMODE = PEVENT
            ENDIF

C           vary order of operation:
            IF (MOD(CT, 3) .EQ. 0) THEN
C              attempt to perform sevrty-type operation:
               CALL PICCHG (WKID, SEVCHR,
     1                      '1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17',
     2                      ACTSEV, ACTPCC)
               CALL DOIMOD (WKID, INCLAS, INMODE)
            ELSE
               CALL DOIMOD (WKID, INCLAS, INMODE)
C              attempt to perform sevrty-type operation:
               CALL PICCHG (WKID, SEVCHR,
     1                      '1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17',
     2                      ACTSEV, ACTPCC)
            ENDIF

            IF (ACTSEV .NE. SEVRTY) THEN
               CALL INMSG ('Cannot do picture change for severity = ' //
     1              SEVLBL // '; skipping test for '                  //
     2              CNDLBL(1:ITRIM(CNDLBL)) // ', and input underway.')
               GOTO 315
            ENDIF

            CALL SETMSG ('3 4 13 30 31', 'If '                        //
     1                   CNDLBL(1:ITRIM(CNDLBL)) // ', and input is ' //
     2                   'underway, then current state of visual '    //
     3                   'representation should be CORRECT.')
            CALL IFPF (SVROK (WKID, 'c'))
C  skip_sev:
315         CONTINUE
C           deactivate local input:
            CALL DOIMOD (WKID, INCLAS, PREQU)
         ELSE
C input is inactive
C           attempt to perform sevrty-type operation:
            CALL PICCHG (WKID, SEVCHR,
     1                   '1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17',
     2                   ACTSEV, ACTPCC)

            IF (ACTSEV .NE. SEVRTY) THEN
               CALL INMSG ('Cannot do picture change for severity = ' //
     1              SEVLBL // '; skipping test for '                  //
     2              CNDLBL(1:ITRIM(CNDLBL)) // ', and input inactive.')
               GOTO 325
            ENDIF

C           calculate expected SVRs and label,
C              based on sevrty, defmod, modmd, and oldsvr
            CALL EXSVRS (SEVRTY, PBNIL, MODMD, OLDSVR, SVRLST, EXPLBL)

            CALL SETMSG ('3 4 13 14 15 31', 'If '                     //
     1                   CNDLBL(1:ITRIM(CNDLBL)) // ', and input is ' //
     2                   'inactive, then current state of visual '    //
     3                   'representation should be '                  //
     4                   EXPLBL (1:ITRIM(EXPLBL)) // '.')
            CALL IFPF (SVROK (WKID, SVRLST))
325         CONTINUE
         ENDIF
C  next instat
360   CONTINUE
C  skip_oldsvr:
355   CONTINUE
C  next oldsvr
350   CONTINUE
C  next sevrty
340   CONTINUE
C  next modmd
330   CONTINUE

666   CONTINUE
      CALL ENDIT
      END

ext: f


text:
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      written by the UFO converter
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      SUBROUTINE MP_COUP3()

      IMPLICIT NONE
      INCLUDE 'model_functions.inc'
      REAL*16 MP__PI, MP__ZERO
      PARAMETER (MP__PI=3.1415926535897932384626433832795E0_16)
      PARAMETER (MP__ZERO=0E0_16)
      INCLUDE 'mp_input.inc'
      INCLUDE 'mp_coupl.inc'

      MP__UV_TMASS_1EPS = MP__MDL_TMASS_UV_1EPS_
      MP__UVWFCT_B_0_1EPS = MP_COND(CMPLX(MP__MDL_MB,KIND=16)
     $ ,CMPLX(0.000000E+00_16,KIND=16),CMPLX(-((MP__MDL_G__EXP__2)
     $ /(2.000000E+00_16*1.600000E+01_16*MP__PI**2))*3.000000E
     $ +00_16*MP__MDL_CF,KIND=16))
      MP__UVWFCT_G_2_1EPS = MP_COND(CMPLX(MP__MDL_MT,KIND=16)
     $ ,CMPLX(0.000000E+00_16,KIND=16),CMPLX(-((MP__MDL_G__EXP__2)
     $ /(2.000000E+00_16*4.800000E+01_16*MP__PI**2))*4.000000E
     $ +00_16*MP__MDL_TF,KIND=16))
      END


text:
	SUBROUTINE TG_ITOF  ( intdtf, iftime, iret )
C************************************************************************
C* TG_ITOF								*
C*									*
C* This subroutine converts three integers containing the date, time	*
C* and forecast field into the two integers stored in a grid file.	*
C*									*
C* TG_ITOF  ( INTDTF, IFTIME, IRET )					*
C*									*
C* Input parameters:							*
C*	INTDTF (3)	INTEGER		Date, time, forecast time	*
C*									*
C* Output parameters:							*
C*	IFTIME (2)	INTEGER		Grid time stored in file	*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 4/89						*
C************************************************************************
	INTEGER		iftime (*), intdtf (*)
C------------------------------------------------------------------------
	iret   = 0
C
C*	If there is no forecast information, the string is stored as
C*	date and time.
C
	IF  ( intdtf (3) .eq. 0 )  THEN
	    iftime (1) = intdtf (1)
	    iftime (2) = intdtf (2)
C
C*	    Otherwise, turn the date and time into a single integer of
C*	    the form MMDDYYHHMM.
C
	  ELSE
	    iyyy = intdtf (1) / 10000
	    mmdd = intdtf (1) - 10000 * iyyy
	    iftime (1) = mmdd * 1000000 + iyyy * 10000 + intdtf (2)
C
C*	    The forecast time remains the same.
C
	    iftime (2) = intdtf (3)
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: sgdecker/pygempak