text:
      SUBROUTINE rkqs(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs)
      character*100 id,revision,revdate
      data id/'$RCSfile: rkqs.f,v $'/
      data revision/'$Revision: 1.4 $'/
      data revdate/'$Date: 1994/03/12 17:59:01 $'/
      INTEGER n
      include "arrays.h"
      REAL*8 eps,hdid,hnext,htry,x,dydx(n),y(n),yscal(n)
      EXTERNAL derivs
CU    USES derivs,rkck
      INTEGER i
      REAL*8 errmax,h,xnew,yerr(IDIMEN),ytemp(IDIMEN),SAFETY,
     $  PGROW,PSHRNK,ERRCON
      PARAMETER (SAFETY=0.9,PGROW=-.2,PSHRNK=-.25,ERRCON=1.89e-4)
      h=htry
1     call rkck(y,dydx,n,x,h,ytemp,yerr,derivs)
      errmax=0.
      do 11 i=1,n
        errmax=max(errmax,abs(yerr(i)/yscal(i)))
11    continue
      errmax=errmax/eps
      if(errmax.gt.1.)then
        h=SAFETY*h*(errmax**PSHRNK)
        if(h.lt.0.1*h)then
          h=.1*h
        endif
        xnew=x+h
        if(xnew.eq.x)pause 'stepsize underflow in rkqs'
        goto 1
      else
        if(errmax.gt.ERRCON)then
          hnext=SAFETY*h*(errmax**PGROW)
        else
          hnext=5.*h
        endif
        hdid=h
        x=x+h
        do 12 i=1,n
          y(i)=ytemp(i)
	if((i.le.dmax).and.(y(i).le.0.))y(i)=1.d-8
12      continue
        return
      endif
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 5"#@-130Rk#3#)KB.

ext: f

max_stars_repo_name: phaustin/empm

max_issues_repo_name: phaustin/empm

max_forks_repo_name: phaustin/empm


text:
      subroutine kmprsnirr(n,z,iz,x,ix,y,iy,count)
C
C
C #####################################################################
C
C     PURPOSE -
C
C        Compress a vector x into vector y where z is not zero.
C
C     INPUT ARGUMENTS -
C
C        n   - (integer) the length of the input vector.
C        z   - (integer) the mask vector.
C        iz  - (integer) the stride in z.
C        x   - (integer) the source vector to be compressed.
C        ix  - (integer) the stride in x.
C        iy    - (integer) the stride in y.
C
C     OUTPUT ARGUMENTS -
C
C        y     - (integer) the compressed output vector.
C        count - (integer) the number of elements in y.
C
C     CHANGE HISTORY -
C
C        $Log: kmprsnirr.f,v $
C        Revision 2.00  2007/11/03 00:49:12  spchu
C        Import to CVS
C
CPVCS
CPVCS       Rev 1.0   21 Mar 2002 10:19:42   dcg
CPVCS    Initial revision.
CPVCS
CPVCS       Rev 1.0   21 Mar 2002 10:19:00   dcg
CPVCS    Initial revision.
CPVCS
CPVCS       Rev 1.2   03/15/95 15:15:46   ejl
CPVCS    Cleaned up, Implicit none.
CPVCS
CPVCS       Rev 1.1   01/04/95 21:55:52   llt
CPVCS    unicos changes (made by het)
CPVCS
CPVCS       Rev 1.0   11/10/94 12:42:44   pvcs
CPVCS    Original version.
C
C ######################################################################
C
      implicit none
C
C ######################################################################
C
      integer n, iz, ix, iy, count
      integer z(*)
      real*8  x(*), y(*)
C
C ######################################################################
C
      integer imax, i, jx, jy
C
C ######################################################################
C
C
C
      imax=iz*n
C
      jx=1
      jy=1
C
      count=0
C
      do i=1,imax,iz
C
         if(z(i).ne.0) then
            count=count+1
            y(jy)=x(jx)
            jy=jy+iy
         endif
C
         jx=jx+ix
C
      enddo
C
      return
      end

max_forks_repo_name: daniellivingston/LaGriT


text:
	PROGRAM  GDPOINT
C************************************************************************
C* PROGRAM GDPOINT							*
C*									*
C* Interpolates any gridded quantity to the location specified and      *
C* prints the output.                                                   *
C*									*
C**									*
C* Log:									*
C* T.W. Barker/WR/MSO	 9/96	Created from GDTHGT				*
C* S. Chiswell/Unidata	11/03	Fixed typo in gr_sscl call			*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(LLMXLN)	gdfile, gpoint, gdatim, scale, gfunc,
     +                  glevel, gvcord, pfunc, parm
C
	REAL	        rgx(1), rgy(1), rlat(1), rlon(1),
     +			hgrd(LLMXGD), yy(1), rnvblk ( LLNNAV )
	CHARACTER	timfnd*36,
     +			time(2)*20, firstm*20, lastim*20
	INTEGER		lev(2)
	LOGICAL		respnd, done, proces, first, gintp, gottm

	INCLUDE		'ERMISS.FNC'
C*
C-----------------------------------------------------------------------
C*	Initialize TAE and GEMPLT.
C
	CALL IP_INIT  ( respnd, iperr )
	CALL IP_IDNT  ( 'GDPOINT', ier )
	IF ( iperr .eq. 0 )  THEN
	   mode = 1
	   CALL GG_INIT  ( mode , iperr )
	   IF (iperr.eq.0) THEN
	      done = .false.
	   ELSE
	      done = .true.
	   ENDIF
	ELSE
	   done = .true.
	END IF
C
C*	Initialize the DG library.
C
	CALL DG_INTL ( ier )
C
C*	Main loop to read in TAE parameters
C
	DO WHILE  ( .not. done )
C
C*	    Set flag to indicate processing will be done.
C
	    proces = .true.
C
C*	    Read in the variables from the TAE.
C
	    CALL GDPTIN  ( gdatim, gpoint, gdfile, scale,
     +                     glevel, gvcord, gfunc, iperr )
C
C*	    Exit if there is an error.
C
	    IF ( iperr .ne. 0 )  THEN
		done = .true.
	    ELSE
C
C*		Process the GDFILE input.
C
		CALL DG_NFIL ( gdfile, ' ', ier )
		IF ( ier .ne. 0 ) THEN
                        CALL ER_WMSG ( 'DG', ier, ' ', irr )
                        proces = .false.
                END IF
C
C*		Process the GDATTIM input; setup the time server.
C
		CALL DG_NDTM ( gdatim, ier )
		IF ( ier .ne. 0 ) THEN
                    CALL ER_WMSG ( 'DG', ier, gdatim, irr )
                    proces = .false.
                END IF
C
C*		Loop over times.
C
		itime = 1
		gottm = proces
		first=.true.
		DO WHILE ( gottm )
C
C*		   Get the next time to process from the time server.
C
		   CALL DG_NTIM ( .true., .true., time, gottm, iret )
		   proces = ( iret .eq. 0 .and. gottm )
		   IF ( iret .ne. 0 ) THEN
		      write (*,*) 'Time not found ',time(1)
		   END IF

		   CALL TG_DUAL ( time, timfnd, iret )
		   IF ( proces .and. first ) THEN
                       CALL DG_FLNO  ( gfunc, iflnos, ier )
                       CALL DG_QDTM ( iflnos, firstm, lastim, ier )
                       IF ( firstm .eq. ' ' ) THEN
                           proces = .false.
                       END IF
                   END IF
C
C*		   Find plotting location
C
		   IF ( proces ) THEN
		       CALL GR_PLOC ( gpoint, rgx (1), rgy (1),
     +				rlat (1), rlon (1), iret )
		       IF ( iret .ne. 0) THEN
			   CALL ER_WMSG ( 'GR', iret, ' ', ier )
			   proces = .false.
		       END IF
		   END IF

C
C*      Set a new output reference navigation using gpoint as the center
C*      of a 9 point grid. This allows vector components to be interpolated
C*      to the desired grid point first within DG_GRID before the quantity is
C*      calculated, rather than interpolated from surrounding points after
C*      the computation is done. This will improve calculations involving
C*      vector components.
C*
C*	Note: This isn't good for functions such as HIGH(), LOWS(), GWFS(),
C*	SM5S(), SM9S(), which use an area of points
C
	           CALL GR_MNAV ( 'CED', 3, 3, rlat(1) - .01,
     +			rlon(1) - .01, rlat(1) + .01,
     +			rlon(1) + .01, 0., 0., 0., .true.,
     +          	rnvblk, iret )
        	   IF ( iret .eq. 0 ) CALL DGCONAV ( rnvblk, iret )
        	   IF ( iret .ne. 0 ) THEN
		       write(*,*) 'Using interpolation of grid results.'
            	       gintp = .true.
                     ELSE
		       gintp = .false.
        	   END IF
C
		   CALL CWRAP_GCALC ( timfnd, glevel, gvcord,
     +			gfunc, scale, gintp, rgy, rgy, yy, iret )
C
c		   CALL DG_GRID(timfnd, glevel, gvcord, gfunc,
c     +                          pfunc, hgrd, kx, ky, time,
c     +                          lev, ivc, parm, iret)
C
                   IF (iret .eq. 0) THEN
c                      CALL IN_SCAL(scale,iscals,iscalv,iret)
c                      CALL GR_SSCL(iscals,kx,ky,1,1,kx,ky,hgrd,
c     +				rmin, rmax, iret)
C
c		      IF ( gintp ) THEN
c                          CALL GR_INTP(1,rgx,rgy,1,kx,ky,hgrd,yy, iret)
c			else
c			  yy(1) = hgrd(5)
c		      END IF
                      IF (first) THEN
                          print 200
                          print 300,gdfile
                          print 301,gpoint
                          print 302,gvcord
                          print 303,glevel
                          print 304,gfunc
                          print 305,scale
 200                      format(' ')
 300                      format('GDFILE: ',a50)
 301                      format('GPOINT: ',a50)
 302                      format('GVCORD: ',a50)
 303                      format('GLEVEL: ',a50)
 304                      format('GFUNC : ',a50)
 305                      format('SCALE : ',a50)
                          first=.false.
                       END IF
                       print 201,time(1),yy(1)
 201                   format('     ',a15,' : ',f12.5)
                   END IF
                END DO

		CALL GDPTUP ( gdatim, gpoint, gdfile, scale,
     +                        glevel, gvcord, gfunc, iret)
		CALL IP_DYNM ( done, ier )
	    END IF
	END DO
	CALL GENDP ( 1, iret )
	CALL IP_EXIT ( iret )

	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      subroutine dpodi(a,lda,n,det,job)
      integer lda,n,job
      double precision a(lda,1)
      double precision det(2)
c
c     dpodi computes the determinant and inverse of a certain
c     double precision symmetric positive definite matrix (see below)
c     using the factors computed by dpoco, dpofa or dqrdc.
c
c     on entry
c
c        a       double precision(lda, n)
c                the output  a  from dpoco or dpofa
c                or the output  x  from dqrdc.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c        job     integer
c                = 11   both determinant and inverse.
c                = 01   inverse only.
c                = 10   determinant only.
c
c     on return
c
c        a       if dpoco or dpofa was used to factor  a  then
c                dpodi produces the upper half of inverse(a) .
c                if dqrdc was used to decompose  x  then
c                dpodi produces the upper half of inverse(trans(x)*x)
c                where trans(x) is the transpose.
c                elements of  a  below the diagonal are unchanged.
c                if the units digit of job is zero,  a  is unchanged.
c
c        det     double precision(2)
c                determinant of  a  or of  trans(x)*x  if requested.
c                otherwise not referenced.
c                determinant = det(1) * 10.0**det(2)
c                with  1.0 .le. det(1) .lt. 10.0
c                or  det(1) .eq. 0.0 .
c
c     error condition
c
c        a division by zero will occur if the input factor contains
c        a zero on the diagonal and the inverse is requested.
c        it will not occur if the subroutines are called correctly
c        and if dpoco or dpofa has set info .eq. 0 .
c
c     linpack.  this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,dscal
c     fortran mod
c
c     internal variables
c
      double precision t
      double precision s
      integer i,j,jm1,k,kp1
c
c     compute determinant
c
      if (job/10 .eq. 0) go to 70
         det(1) = 1.0d0
         det(2) = 0.0d0
         s = 10.0d0
         do 50 i = 1, n
            det(1) = a(i,i)**2*det(1)
c        ...exit
            if (det(1) .eq. 0.0d0) go to 60
   10       if (det(1) .ge. 1.0d0) go to 20
               det(1) = s*det(1)
               det(2) = det(2) - 1.0d0
            go to 10
   20       continue
   30       if (det(1) .lt. s) go to 40
               det(1) = det(1)/s
               det(2) = det(2) + 1.0d0
            go to 30
   40       continue
   50    continue
   60    continue
   70 continue
c
c     compute inverse(r)
c
      if (mod(job,10) .eq. 0) go to 140
         do 100 k = 1, n
            a(k,k) = 1.0d0/a(k,k)
            t = -a(k,k)
            call dscal(k-1,t,a(1,k),1)
            kp1 = k + 1
            if (n .lt. kp1) go to 90
            do 80 j = kp1, n
               t = a(k,j)
               a(k,j) = 0.0d0
               call daxpy(k,t,a(1,k),1,a(1,j),1)
   80       continue
   90       continue
  100    continue
c
c        form  inverse(r) * trans(inverse(r))
c
         do 130 j = 1, n
            jm1 = j - 1
            if (jm1 .lt. 1) go to 120
            do 110 k = 1, jm1
               t = a(k,j)
               call daxpy(k,t,a(1,j),1,a(1,k),1)
  110       continue
  120       continue
            t = a(j,j)
            call dscal(j,t,a(1,j),1)
  130    continue
  140 continue
      return
      end


text:
great(h1,kk1).
great(z1,w1).
great(k1,l1).
great(h1,u1).
great(q1,b1).
great(e1,z1).
great(ii1,d1).
great(ee1,n1).
great(r1,p1).
great(z1,x1).
great(i1,o1).
great(p1,w1).
great(j1,jj1).
great(h1,jj1).
great(k1,kk1).
great(j1,i1).
great(g1,q1).
great(l1,ff1).
great(l1,q1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
	SUBROUTINE DSCLR2  ( icolr, icolr2, jcolr, jcolr2, iret )
C************************************************************************
C* DSCLR2								*
C* 									*
C* This subroutine sets two color numbers.  Color numbers larger than	*
C* the number of valid colors will be converted, via modular 		*
C* arithmetic to a valid color number.  If the color number is 		*
C* negative or zero, no change will be made. 				*
C*									*
C* DSCLR2  ( ICOLR, ICOLR2, JCOLR, JCOLR2, IRET )			*
C* 									*
C* Input parameters:							*
C*	ICOLR		INTEGER		Color number 1			*
C*	ICOLR2		INTEGER		Color number 2			*
C* 									*
C* Output parameters:							*
C*	JCOLR		INTEGER		Color number set 1		*
C*	JCOLR2		INTEGER		Color number set 2		*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* I. Durham/GSC	 3/98						*
C************************************************************************
	INCLUDE		'ERROR.PRM'
	INCLUDE		'DEVCHR.CMN'
	INCLUDE		'DEVACT.CMN'
C------------------------------------------------------------------------
	iret  = NORMAL
	jcolr  = mcolr
	jcolr2 = mcolr2
C
C*	Exclude negative and zero colors.
C
	IF  ( icolr .eq. 101 )  THEN
	    jcolr  = icolr
	  ELSE IF  ( icolr .gt. 0 ) THEN
	    jcolr  = MOD ( icolr, nncolr )
	    IF  ( jcolr .eq. 0 ) jcolr = nncolr
	  ELSE
	    RETURN
	END IF
C
C*	Exclude negative and zero colors.
C
	IF  ( icolr2 .eq. 101 )  THEN
	    jcolr2 = icolr2
	  ELSE IF  ( icolr2 .gt. 0 ) THEN
	    jcolr2 = MOD ( icolr2, nncolr )
	    IF  ( jcolr2 .eq. 0 ) jcolr2 = nncolr
	  ELSE
	    RETURN
	END IF
C
C*	Check if this is the active color.
C
	IF  ( jcolr .ne. mcolr ) THEN
C
C*	    Set this color in the device.
C
	    CALL HSCOLR  ( jcolr, iret )
C
C*	    Store color in active buffer.
C
	    mcolr = jcolr
	END IF
C
	IF  ( jcolr2 .ne. mcolr2 ) THEN
C
C*	    Store color in active buffer.
C
	    mcolr2 = jcolr2
	END IF
C
	RETURN
	END


text:
      function mapad(idiag,n,mxc)
c
c     compute coefficient matrix profile and set diagonal pointers
c
      dimension idiag(*)
c
      mxc=0
      lk=mxc
      do 100 i=1,n
      mxc=max(mxc,idiag(i))
      lk=lk+idiag(i)+1
      idiag(i)=lk
  100 continue
      mxc=mxc+1
      mapad=lk
      return
      end

max_stars_repo_name: imohame/LabCode

max_issues_repo_name: imohame/LabCode

max_forks_repo_name: imohame/LabCode


module AirClass
    use fem
    implicit none

    type :: Air_
        type(FEMDomain_) :: femdomain
        real(real64) :: temp = 303.0d0 ! K
        real(real64) :: CO2 = 380.0d0 ! ppm
        real(real64) :: O2 = 202000.0d0 ! ppm
    contains
        procedure, public :: init => initAir
    end type
contains

subroutine initAir(obj,temp,CO2,O2)
    class(Air_),intent(inout) :: obj
    real(real64),optional,intent(in) :: temp,CO2,O2

    obj%temp = input(default=304.0d0, option=temp)
    obj%CO2 = input(default=380.0d0, option=CO2)
    obj%O2 = input(default=202000.0d0, option=O2)

end subroutine

end module


text:
c     ****************************************************************
c     *                                                              *
c     *                      subroutine patch_data                   *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 12/19/2019 rhd             *
c     *                                                              *
c     *     patch command to support development mods of data during *
c     *     execution                                                *
c     *                                                              *
c     ****************************************************************
c
c
c
      subroutine patch_data( sbflg1, sbflg2 )
      use global_data ! old common.main
c
      use main_data
c
      implicit none
c
      logical ::sbflg1, sbflg2
      logical, external :: matchs, numd, integr, endcrd, true, label,
     &                     string, realn
c
c
      if( matchs('props',4) ) then
         call patch_data_props
         return
      else
         write(out,9000)
         call die_abort
      end if
c
 9000 format(2x,"Fatal error: invalid patch command...")
c
      contains
c     ========
c
      subroutine patch_data_props
      implicit none
c
      real :: realval
      integer :: intval, row, col, itype
c

      if( .not. matchs('row',3) ) then
          write(out,9000)
          call die_abort
      end if
      if( .not. integr( row ) ) then
          write(out,9000)
          call die_abort
      end if
      if( .not. matchs('column',3) ) then
          write(out,9000)
          call die_abort
      end if
      if( .not. integr( col ) ) then
          write(out,9000)
          call die_abort
      end if
c
      itype = 1
      if( .not. integr( intval ) ) then
          if( .not. realn( realval ) ) then
            write(out,9000)
            call die_abort
          end if
          itype = 2
      end if
c
      if( itype == 1 ) iprops(row,col) = intval
      if( itype == 2 ) props(row,col)  = realval
c
      return
c
 9000 format(2x,"Fatal error: invalid patch command...")
c
      end subroutine patch_data_props
      end subroutine patch_data


text:
      subroutine findtic (end1,end2,bigtic,smltic)
c******************************************************************************
c     this routine makes tic marks at nice round intervals.
c******************************************************************************

      implicit real*8(a-h,o-z)
      real*4 end1, end2, bigtic, smltic, span, spanlog

      span = end2 - end1
      if (span .lt. 0.) span = - span
      spanlog = alog10(span)
      size = ifix(spanlog)
      if (spanlog .lt. 0.) size = size - 1
      chop = spanlog - size
      if (chop .lt. 0.31) then
         bigtic = 10.**(size)/2
         smltic = bigtic/5
      elseif (chop .lt. 0.71) then
         bigtic = 10.**(size)
         smltic = bigtic/5
      else
         bigtic = 10.**(size)*2
         smltic = bigtic/4
      endif
      return

      end


text:
module mod_fft
  use iso_c_binding , only: C_INT
  use mod_common_mpi, only: ierr
  use mod_fftw_param
  use mod_types
  !$ use omp_lib
  private
  public fftini,fftend,fft
  contains
  subroutine fftini(n_x,n_y,bcxy,c_or_f,arrplan,normfft)
    implicit none
    integer, intent(in), dimension(3) :: n_x,n_y
    character(len=1), intent(in), dimension(0:1,2) :: bcxy
    character(len=1), intent(in), dimension(2) :: c_or_f
    type(C_PTR), intent(out), dimension(2,2) :: arrplan
    real(rp), intent(out) :: normfft
    real(rp), dimension(n_x(1),n_x(2),n_x(3))  :: arrx
    real(rp), dimension(n_y(1),n_y(2),n_y(3))  :: arry
    type(C_PTR) :: plan_fwd_x,plan_bwd_x, &
                   plan_fwd_y,plan_bwd_y
    type(fftw_iodim), dimension(1) :: iodim
    type(fftw_iodim), dimension(2) :: iodim_howmany
    integer :: kind_fwd,kind_bwd
    real(rp), dimension(2) :: norm
    integer(C_INT) :: nx_x,ny_x,nz_x, &
                      nx_y,ny_y,nz_y
    integer :: ix,iy
#if defined(_SINGLE_PRECISION)
    !$ call sfftw_init_threads(ierr)
    !$ call sfftw_plan_with_nthreads(omp_get_max_threads())
#else
    !$ call dfftw_init_threads(ierr)
    !$ call dfftw_plan_with_nthreads(omp_get_max_threads())
#endif
    !
    ! fft in x
    !
    ! prepare plans with guru interface
    !
    nx_x = n_x(1)
    ny_x = n_x(2)
    nz_x = n_x(3)
    nx_y = n_y(1)
    ny_y = n_y(2)
    nz_y = n_y(3)
    !
    normfft = 1.
    ix = 0
    ! size of transform reduced by 1 point with Dirichlet BC in face
    if(bcxy(0,1)//bcxy(1,1) == 'DD'.and.c_or_f(1) == 'f') ix = 1
    iodim(1)%n  = nx_x-ix
    iodim(1)%is = 1
    iodim(1)%os = 1
    iodim_howmany(1)%n  = ny_x
    iodim_howmany(1)%is = nx_x
    iodim_howmany(1)%os = nx_x
    iodim_howmany(2)%n  = nz_x
    iodim_howmany(2)%is = nx_x*ny_x
    iodim_howmany(2)%os = nx_x*ny_x
    call find_fft(bcxy(:,1),c_or_f(1),kind_fwd,kind_bwd,norm)
    plan_fwd_x=fftw_plan_guru_r2r(1,iodim,2,iodim_howmany,arrx,arrx,kind_fwd,FFTW_ESTIMATE)
    plan_bwd_x=fftw_plan_guru_r2r(1,iodim,2,iodim_howmany,arrx,arrx,kind_bwd,FFTW_ESTIMATE)
    normfft = normfft*norm(1)*(nx_x+norm(2)-ix)
    !
    ! fft in y
    !
    ! prepare plans with guru interface
    !
    iy = 0
    ! size of transform reduced by 1 point with Dirichlet BC in face
    if(bcxy(0,2)//bcxy(1,2) == 'DD'.and.c_or_f(2) == 'f') iy = 1
    iodim(1)%n  = ny_y-iy
    iodim(1)%is = nx_y
    iodim(1)%os = nx_y
    iodim_howmany(1)%n  = nx_y
    iodim_howmany(1)%is = 1
    iodim_howmany(1)%os = 1
    iodim_howmany(2)%n  = nz_y
    iodim_howmany(2)%is = nx_y*ny_y
    iodim_howmany(2)%os = nx_y*ny_y
    call find_fft(bcxy(:,2),c_or_f(2),kind_fwd,kind_bwd,norm)
    plan_fwd_y=fftw_plan_guru_r2r(1,iodim,2,iodim_howmany,arry,arry,kind_fwd,FFTW_ESTIMATE)
    plan_bwd_y=fftw_plan_guru_r2r(1,iodim,2,iodim_howmany,arry,arry,kind_bwd,FFTW_ESTIMATE)
    normfft = normfft*norm(1)*(ny_y+norm(2)-iy)
    !
    normfft = normfft**(-1)
    arrplan(1,1) = plan_fwd_x
    arrplan(2,1) = plan_bwd_x
    arrplan(1,2) = plan_fwd_y
    arrplan(2,2) = plan_bwd_y
  end subroutine fftini
  !
  subroutine fftend(arrplan)
    implicit none
    type(C_PTR), intent(in), dimension(2,2) :: arrplan
#if defined(_SINGLE_PRECISION)
    call sfftw_destroy_plan(arrplan(1,1))
    call sfftw_destroy_plan(arrplan(1,2))
    call sfftw_destroy_plan(arrplan(2,1))
    call sfftw_destroy_plan(arrplan(2,2))
    !$call sfftw_cleanup_threads(ierr)
#else
    call dfftw_destroy_plan(arrplan(1,1))
    call dfftw_destroy_plan(arrplan(1,2))
    call dfftw_destroy_plan(arrplan(2,1))
    call dfftw_destroy_plan(arrplan(2,2))
    !$call dfftw_cleanup_threads(ierr)
#endif
  end subroutine fftend
  !
  subroutine fft(plan,arr)
    implicit none
    type(C_PTR), intent(in) :: plan
    real(rp), intent(inout), dimension(:,:,:) :: arr
#if defined(_SINGLE_PRECISION)
    call sfftw_execute_r2r(plan,arr,arr)
#else
    call dfftw_execute_r2r(plan,arr,arr)
#endif
  end subroutine fft
  !
  subroutine find_fft(bc,c_or_f,kind_fwd,kind_bwd,norm)
  implicit none
  character(len=1), intent(in), dimension(0:1) :: bc
  character(len=1), intent(in) :: c_or_f
  integer , intent(out) :: kind_fwd,kind_bwd
  real(rp), intent(out), dimension(2) :: norm
  if(c_or_f == 'c') then
    select case(bc(0)//bc(1))
    case('PP')
      kind_fwd = FFTW_R2HC
      kind_bwd = FFTW_HC2R
      norm = [1.,0.]
    case('NN')
      kind_fwd = FFTW_REDFT10
      kind_bwd = FFTW_REDFT01
      norm = [2.,0.]
    case('DD')
      kind_fwd = FFTW_RODFT10
      kind_bwd = FFTW_RODFT01
      norm = [2.,0.]
    case('ND')
      kind_fwd = FFTW_REDFT11
      kind_bwd = FFTW_REDFT11
      norm = [2.,0.]
    case('DN')
      kind_fwd = FFTW_RODFT11
      kind_bwd = FFTW_RODFT11
      norm = [2.,0.]
    end select
  else if(c_or_f == 'f') then
    select case(bc(0)//bc(1))
    case('PP')
      kind_fwd = FFTW_R2HC
      kind_bwd = FFTW_HC2R
      norm = [1.,0.]
    case('NN')
      kind_fwd = FFTW_REDFT00
      kind_bwd = FFTW_REDFT00
      norm = [2.,-1.]
    case('DD')
      kind_fwd = FFTW_RODFT00
      kind_bwd = FFTW_RODFT00
      norm = [2.,1.]
    case('ND')
      kind_fwd = FFTW_REDFT10
      kind_bwd = FFTW_REDFT01
      norm = [2.,0.]
    case('DN')
      kind_fwd = FFTW_RODFT01
      kind_bwd = FFTW_RODFT10
      norm = [2.,0.]
    end select
  end if
  end subroutine find_fft
end module mod_fft

max_stars_repo_name: GabrieleBoga/CaNS

max_issues_repo_name: GabrieleBoga/CaNS

max_forks_repo_name: GabrieleBoga/CaNS


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2021 Guido Dhondt
!
!     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 2);
!
!
!     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, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine precfdcyc(nelemface,sideface,nface,ipoface,nodface,
     &  ne,ipkon,kon,lakon,ikboun,ilboun,xboun,nboun,nk,isolidsurf,
     &  nsolidsurf,ifreestream,nfreestream,neighsolidsurf,iponoel,inoel,
     &  inoelfree,nef,co,ipompc,nodempc,ikmpc,ilmpc,nmpc,set,istartset,
     &  iendset,ialset,nset,iturbulent)
!
!     cyclic update for cfd applicatons:
!     - determining the external faces of the mesh and storing
!       them in fields nelemface and sideface
!     - determining the fluid elements belonging to a given node
!       and storing them in fields iponoel and inoel
!
      implicit none
!
      logical solidboun,solidmpc,freestream
!
      character*1 sideface(*)
      character*8 lakon(*)
      character*81 set(*),noset
!
      integer nelemface(*),nface,ipoface(*),nodface(5,*),nodes(4),
     &  ne,ipkon(*),kon(*),indexe,ifaceq(8,6),ifacet(7,4),index,
     &  ifacew(8,5),ithree,ifour,iaux,kflag,nnodes,ikboun(*),
     &  ilboun(*),nboun,isolidsurf(*),nsolidsurf,ifreestream(*),
     &  nfreestream,id,nk,node,idof,i,j,k,l,m,neighsolidsurf(*),
     &  iponoel(*),noden,idn,nope,nodemin,ifree,nef,indexold,
     &  inoel(3,*),ifreenew,inoelfree,ikmpc(*),nmpc,indexi,
     &  nodempc(3,*),ipompc(*),ilmpc(*),kmax,impc,idofi,idi,
     &  iturbulent,istartset(*),iendset(*),ialset(*),nset
!
      real*8 xboun(*),dist,distmin,co(3,*)
!
!     nodes belonging to the element faces
!
      data ifaceq /4,3,2,1,11,10,9,12,
     &            5,6,7,8,13,14,15,16,
     &            1,2,6,5,9,18,13,17,
     &            2,3,7,6,10,19,14,18,
     &            3,4,8,7,11,20,15,19,
     &            4,1,5,8,12,17,16,20/
      data ifacet /1,3,2,7,6,5,11,
     &             1,2,4,5,9,8,12,
     &             2,3,4,6,10,9,13,
     &             1,4,3,8,10,7,14/
      data ifacew /1,3,2,9,8,7,0,0,
     &             4,5,6,10,11,12,0,0,
     &             1,2,5,4,7,14,10,13,
     &             2,3,6,5,8,15,11,14,
     &             4,6,3,1,12,15,9,13/
!
      kflag=1
      ithree=3
      ifour=4
!
!     determining the external element faces of the fluid mesh
!     the faces are catalogued by the three lowes nodes numbers
!     in ascending order. ipoface(i) points to a face for which
!     node i is the lowest node and nodface(1,ipoface(i)) and
!     nodface(2,ipoface(i)) are the next lower ones.
!     nodface(3,ipoface(i)) contains the element number,
!     nodface(4,ipoface(i)) the face number and nodface(5,ipoface(i))
!     is a pointer to the next surface for which node i is the
!     lowest node; if there are no more such surfaces the pointer
!     has the value zero
!     An external element face is one which belongs to one element
!     only
!
      ifree=1
      do i=1,6*nef-1
         nodface(5,i)=i+1
      enddo
      do i=1,ne
         if(ipkon(i).lt.0) cycle
         if(lakon(i)(1:1).ne.'F') cycle
         indexe=ipkon(i)
         if((lakon(i)(4:4).eq.'2').or.(lakon(i)(4:4).eq.'8')) then
            do j=1,6
               do k=1,4
                  nodes(k)=kon(indexe+ifaceq(k,j))
               enddo
               call isortii(nodes,iaux,ifour,kflag)
               indexold=0
               index=ipoface(nodes(1))
               do
!
!                 adding a surface which has not been
!                 catalogued so far
!
                  if(index.eq.0) then
                     ifreenew=nodface(5,ifree)
                     nodface(1,ifree)=nodes(2)
                     nodface(2,ifree)=nodes(3)
                     nodface(3,ifree)=i
                     nodface(4,ifree)=j
                     nodface(5,ifree)=ipoface(nodes(1))
                     ipoface(nodes(1))=ifree
                     ifree=ifreenew
                     exit
                  endif
!
!                 removing a surface which has already
!                 been catalogued
!
                  if((nodface(1,index).eq.nodes(2)).and.
     &               (nodface(2,index).eq.nodes(3))) then
                     if(indexold.eq.0) then
                        ipoface(nodes(1))=nodface(5,index)
                     else
                        nodface(5,indexold)=nodface(5,index)
                     endif
                     nodface(5,index)=ifree
                     ifree=index
                     exit
                  endif
                  indexold=index
                  index=nodface(5,index)
               enddo
            enddo
         elseif((lakon(i)(4:4).eq.'4').or.(lakon(i)(4:5).eq.'10')) then
            do j=1,4
               do k=1,3
                  nodes(k)=kon(indexe+ifacet(k,j))
               enddo
               call isortii(nodes,iaux,ithree,kflag)
               indexold=0
               index=ipoface(nodes(1))
               do
!
!                 adding a surface which has not been
!                 catalogues so far
!
                  if(index.eq.0) then
                     ifreenew=nodface(5,ifree)
                     nodface(1,ifree)=nodes(2)
                     nodface(2,ifree)=nodes(3)
                     nodface(3,ifree)=i
                     nodface(4,ifree)=j
                     nodface(5,ifree)=ipoface(nodes(1))
                     ipoface(nodes(1))=ifree
                     ifree=ifreenew
                     exit
                  endif
!
!                 removing a surface which has already
!                 been catalogued
!
                  if((nodface(1,index).eq.nodes(2)).and.
     &               (nodface(2,index).eq.nodes(3))) then
                     if(indexold.eq.0) then
                        ipoface(nodes(1))=nodface(5,index)
                     else
                        nodface(5,indexold)=nodface(5,index)
                     endif
                     nodface(5,index)=ifree
                     ifree=index
                     exit
                  endif
                  indexold=index
                  index=nodface(5,index)
               enddo
            enddo
         else
            do j=1,5
               if(j.le.2) then
                  do k=1,3
                     nodes(k)=kon(indexe+ifacew(k,j))
                  enddo
                  call isortii(nodes,iaux,ithree,kflag)
               else
                  do k=1,4
                     nodes(k)=kon(indexe+ifacew(k,j))
                  enddo
                  call isortii(nodes,iaux,ifour,kflag)
               endif
               indexold=0
               index=ipoface(nodes(1))
               do
!
!                 adding a surface which has not been
!                 catalogues so far
!
                  if(index.eq.0) then
                     ifreenew=nodface(5,ifree)
                     nodface(1,ifree)=nodes(2)
                     nodface(2,ifree)=nodes(3)
                     nodface(3,ifree)=i
                     nodface(4,ifree)=j
                     nodface(5,ifree)=ipoface(nodes(1))
                     ipoface(nodes(1))=ifree
                     ifree=ifreenew
                     exit
                  endif
!
!                 removing a surface which has already
!                 been catalogued
!
                  if((nodface(1,index).eq.nodes(2)).and.
     &               (nodface(2,index).eq.nodes(3))) then
                     if(indexold.eq.0) then
                        ipoface(nodes(1))=nodface(5,index)
                     else
                        nodface(5,indexold)=nodface(5,index)
                     endif
                     nodface(5,index)=ifree
                     ifree=index
                     exit
                  endif
                  indexold=index
                  index=nodface(5,index)
               enddo
            enddo
         endif
      enddo
!
!     storing the external faces in nelemface and sideface
!
      nface=0
!
      do m=1,nk
         index=ipoface(m)
         do
            if(index.eq.0) exit
            nface=nface+1
            i=nodface(3,index)
            j=nodface(4,index)
!
            nelemface(nface)=i
            write(sideface(nface)(1:1),'(i1)') j
!
            index=nodface(5,index)
         enddo
      enddo
!
!     determining the fluid elements belonging to edge nodes of
!     the elements
!
      inoelfree=1
      do i=1,ne
         if(ipkon(i).lt.0) cycle
         if(lakon(i)(1:1).ne.'F') cycle
         if(lakon(i)(4:4).eq.'2') then
            nope=20
         elseif(lakon(i)(4:4).eq.'8') then
            nope=8
         elseif(lakon(i)(4:4).eq.'4') then
            nope=4
         elseif(lakon(i)(4:5).eq.'10') then
            nope=10
         elseif(lakon(i)(4:4).eq.'6') then
            nope=6
         else
            nope=15
         endif
         indexe=ipkon(i)
         do j=1,nope
            node=kon(indexe+j)
            inoel(1,inoelfree)=i
            inoel(2,inoelfree)=j
            inoel(3,inoelfree)=iponoel(node)
            iponoel(node)=inoelfree
            inoelfree=inoelfree+1
         enddo
      enddo
!
!     sorting nelemface
!
      kflag=2
      call isortic(nelemface,sideface,nface,kflag)
c      write(*,*) 'external faces'
c      do i=1,nface
c         write(*,*) nelemface(i),sideface(i)
c      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
      subroutine lattc(rb,qb)
c  sets up the real and reciprocal space lattice vectors
      implicit none
      real*8 rb(3,3), qb(3,3)
      real*8 vol
      integer m,k
      real*8 tripl
      external tripl

      call cross(rb(1,2), rb(1,3), qb(1,1))
      call cross(rb(1,3), rb(1,1), qb(1,2))
      call cross(rb(1,1), rb(1,2), qb(1,3))
      vol = tripl(rb(1,1), rb(1,2), rb(1,3))
      qb = qb * (1.d0/vol)

      return
      end


text:
!-----------------------------------------------------------------------
!
! MODULE: geom_module
!> @brief
!> This module contains the variables that relate to the geometry of the
!> problem and the subroutines necessary to allocate and deallocate
!> geometry related data as necessary.
!
!-----------------------------------------------------------------------

MODULE geom_module

  USE global_module, ONLY: i_knd, r_knd, zero, one, two, allocator

  use iso_c_binding
  use sicm_f90, only: sf_alloc, sf_free

  IMPLICIT NONE

  PUBLIC

  SAVE
!_______________________________________________________________________
!
! Module Input Variables
!
! ndimen   - number of spatial dimensions 1/2/3
! nx       - number of x-dir spatial cells (global)
! ny       - number of y-dir spatial cells (global on input, reset to
!            per PE in setup)
! nz       - number of z-dir spatial cells (global on input, reset to
!            per PE in setup)
! lx       - total length of x domain
! ly       - total length of y domain
! lz       - total length of z domain
!_______________________________________________________________________

  INTEGER(i_knd) :: ndimen=1, nx=4, ny=1, nz=1

  REAL(r_knd) :: lx=one, ly=one, lz=one
!_______________________________________________________________________
!
! Run-time variables
!
! ny_gl    - global number of y-dir spatial cells
! nz_gl    - global number of z-dir spatial cells
! jlb      - global index of local lower y bound
! jub      - global index of local upper y bound
! klb      - global index of local lower z bound
! kub      - global index of local upper z bound
!
! dx       - x width of spatial cell
! dy       - y width of spatial cell
! dz       - z width of spatial cell
!
! nc       - number of i-chunks, nx/ichunk
!
! jdim     - number of directions to sweep in j-direction: 1D=1, else=2
! kdim     - number of directions to sweep in z-direction: 3D=2, else=1
!
! hi       - Spatial DD x-coefficient
! hj       - Spatial DD y-coefficient
! hk       - Spatial DD z-coefficient
!
! dinv(nang,ichunk,ny,nz,nc,ng) - Sweep denominator pre-computed
!
! ndiag    - number of diagonals of mini-KBA sweeps in nested threading
!_______________________________________________________________________

  INTEGER(i_knd) :: ny_gl, nz_gl, jlb, jub, klb, kub, nc, ndiag=0,     &
    jdim, kdim

  REAL(r_knd) :: dx, dy, dz, hi, hj, hk

  REAL(r_knd), pointer, DIMENSION(:,:,:,:,:,:) :: dinv
  type(c_ptr) :: c_tmp
  integer(c_size_t) :: dinvsize
!_______________________________________________________________________
!
! Derived data types for mini-KBA diagonals
!
! cell_id_type      - type for holding ijk indices of a cell on a given
!                     diagonal
! ic, j, k          - ijk integer indices of cell_id_type
! diag_type         - type for holding diagonal information
! len               - number of cells on a diagonal line/plane
! cell_id(len)      - array of cell ijk indices
! diag(ndiag)       - array of diagonal lines/planes information
!_______________________________________________________________________

  TYPE cell_id_type
    INTEGER(i_knd) :: ic, j, k
  END TYPE cell_id_type

  TYPE diag_type
    INTEGER(i_knd) :: len
    TYPE(cell_id_type), ALLOCATABLE, DIMENSION(:) :: cell_id
  END TYPE diag_type

  TYPE(diag_type), ALLOCATABLE, DIMENSION(:) :: diag

  CONTAINS

  SUBROUTINE geom_allocate ( nang, ng, swp_typ, ichunk, ierr )

!-----------------------------------------------------------------------
!
! Allocate the geometry-related solution arrays. Setup the diag_type
! data structures for mini-KBA if requested.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: nang, ng, swp_typ, ichunk

    INTEGER(i_knd), INTENT(OUT) :: ierr
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: i, j, k, nn, ing

    INTEGER(i_knd), ALLOCATABLE, DIMENSION(:) :: indx
!_______________________________________________________________________

    ierr = 0

    dinvsize = int(nang, c_size_t) * int(ichunk, c_size_t) * int(ny, c_size_t)&
      * int(nz, c_size_t) * int(nc, c_size_t) * int(ng, c_size_t) * int(r_knd, c_size_t)
    c_tmp = sf_alloc(allocator, dinvsize)
    call c_f_pointer(c_tmp, dinv, shape=[nang,ichunk,ny,nz,nc,ng])
    !ALLOCATE( dinv(nang,ichunk,ny,nz,nc,ng), STAT=ierr )
    IF ( ierr /= 0 ) RETURN

    hi = zero
    hj = zero
    hk = zero
    dinv = zero

!_______________________________________________________________________
!
!   If mini-KBA on spatial chunks is selected, set up the diagonal
!   'diag' derived data type, allocate many diagonals, and set indices
!   for cells. Because each chunk is same size, can create just one.
!_______________________________________________________________________

    IF ( swp_typ == 1 ) THEN

      ndiag = ichunk + ny + nz - 2

      ALLOCATE( diag(ndiag), indx(ndiag), STAT=ierr )
      IF ( ierr /= 0 ) RETURN

      diag%len = 0
      indx = 0
!_______________________________________________________________________
!
!     Cells of same diagonal all have same value according to i+j+k-2
!     formula. Use that to compute len for each diagonal. Use ichunk.
!_______________________________________________________________________

      DO k = 1, nz
      DO j = 1, ny
      DO i = 1, ichunk
        nn = i + j + k - 2
        diag(nn)%len = diag(nn)%len + 1
      END DO
      END DO
      END DO
!_______________________________________________________________________
!
!     Next allocate cell_id array within diag type according to len
!_______________________________________________________________________

      DO nn = 1, ndiag
        ing = diag(nn)%len
        ALLOCATE( diag(nn)%cell_id(ing), STAT=ierr )
        IF ( ierr /= 0 ) RETURN
      END DO
!_______________________________________________________________________
!
!     Lastly, set each cell's actual ijk indices in this diagonal map
!_______________________________________________________________________

      DO k = 1, nz
      DO j = 1, ny
      DO i = 1, ichunk
        nn = i + j + k - 2
        indx(nn) = indx(nn) + 1
        ing = indx(nn)
        diag(nn)%cell_id(ing)%ic = i
        diag(nn)%cell_id(ing)%j  = j
        diag(nn)%cell_id(ing)%k  = k
      END DO
      END DO
      END DO

      DEALLOCATE( indx )

    END IF
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE geom_allocate

  SUBROUTINE geom_deallocate ( swp_typ )

!-----------------------------------------------------------------------
!
! Deallocate the geometry-related solution arrays.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: swp_typ
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: i
!_______________________________________________________________________
!
!   Deallocate the sweep parameters
!_______________________________________________________________________

    call sf_free(allocator, c_loc(dinv), dinvsize)
    !DEALLOCATE( dinv )
!_______________________________________________________________________
!
!   Deallocate the diagonal related arrays with swp_typ
!_______________________________________________________________________

    IF ( swp_typ == 1 ) THEN
      DO i = 1, ndiag
        DEALLOCATE( diag(i)%cell_id )
      END DO
      DEALLOCATE( diag )
    END IF
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE geom_deallocate

  SUBROUTINE geom_param_calc ( nang, ichunk, mu, eta, xi, cs, vd, d )

!-----------------------------------------------------------------------
!
! Calculate the DD spatial coefficients hi, hj, hk for all angles at
! the start of each time step. Compute the pre-computed/inverted dinv.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: nang, ichunk

    REAL(r_knd), INTENT(IN) :: vd

    REAL(r_knd), DIMENSION(nang), INTENT(IN) :: mu, eta, xi

    REAL(r_knd), DIMENSION(nx,ny,nz), INTENT(IN) :: cs

    REAL(r_knd), DIMENSION(nang,ichunk,ny,nz,nc), INTENT(OUT) :: d
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: k, j, ii, ich, i, m
!_______________________________________________________________________
!
!   Set the DD coefficients
!_______________________________________________________________________

    hi = two/dx
    IF ( ndimen > 1 ) THEN
      hj = two/dy
      IF ( ndimen > 2 ) THEN
        hk = two/dz
      END IF
    END IF
!_______________________________________________________________________
!
!   Compute the inverted denominator, saved for sweep
!_______________________________________________________________________

    DO k = 1, nz
      DO j = 1, ny
        ii = 0
        ich = 1
        DO i = 1, nx
          ii = ii + 1
          DO m = 1, nang
            d(m,ii,j,k,ich) = one / ( cs(i,j,k) + vd + hi*mu(m) +      &
                                      hj*eta(m) + hk*xi(m) )
          END DO
          IF ( ii == ichunk ) THEN
            ii = 0
            ich = ich + 1
          END IF
        END DO
      END DO
    END DO
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE geom_param_calc

END MODULE geom_module

max_stars_repo_name: gvallee/SICM

max_issues_repo_name: gvallee/SICM

max_forks_repo_name: gvallee/SICM


text:
!
! MWwaterLUT_Define
!
! Module defining the MWwaterLUT object containing the
! Look-Up Table (LUT) for the microWave (MW) sea surface emissivity
! model.
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 10-Nov-2011
!                       paul.vandelst@noaa.gov
!

MODULE MWwaterLUT_Define

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module use
  USE Type_Kinds           , ONLY: fp, Long, Double
  USE Message_Handler      , ONLY: SUCCESS, FAILURE, INFORMATION, Display_Message
  USE Compare_Float_Numbers, ONLY: OPERATOR(.EqualTo.)
  USE File_Utility         , ONLY: File_Open, File_Exists
  USE Binary_File_Utility  , ONLY: Open_Binary_File      , &
                                   WriteGAtts_Binary_File, &
                                   ReadGAtts_Binary_File
  ! Disable implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  ! Everything private by default
  PRIVATE
  ! Datatypes
  PUBLIC :: MWwaterLUT_type
  ! Operators
  PUBLIC :: OPERATOR(==)
  ! Procedures
  PUBLIC :: MWwaterLUT_Associated
  PUBLIC :: MWwaterLUT_Destroy
  PUBLIC :: MWwaterLUT_Create
  PUBLIC :: MWwaterLUT_Inspect
  PUBLIC :: MWwaterLUT_ValidRelease
  PUBLIC :: MWwaterLUT_Info
  PUBLIC :: MWwaterLUT_DefineVersion
  PUBLIC :: MWwaterLUT_InquireFile
  PUBLIC :: MWwaterLUT_ReadFile
  PUBLIC :: MWwaterLUT_WriteFile

  ! ---------------------
  ! Procedure overloading
  ! ---------------------
  INTERFACE OPERATOR(==)
    MODULE PROCEDURE MWwaterLUT_Equal
  END INTERFACE OPERATOR(==)

  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
    '$Id: MWwaterLUT_Define.f90 99117 2017-11-27 18:37:14Z tong.zhu@noaa.gov $'
  ! Release and version
  INTEGER, PARAMETER :: MWWATERLUT_RELEASE = 1  ! This determines structure and file formats.
  INTEGER, PARAMETER :: MWWATERLUT_VERSION = 1  ! This is just the default data version.
  ! Close status for write errors
  CHARACTER(*), PARAMETER :: WRITE_ERROR_STATUS = 'DELETE'
  ! Literal constants
  REAL(fp), PARAMETER :: ZERO = 0.0_fp
  REAL(fp), PARAMETER :: ONE  = 1.0_fp
  ! String lengths
  INTEGER,  PARAMETER :: ML = 256 ! Message length
  INTEGER,  PARAMETER :: SL =  80 ! String length

  ! ----------------------------------
  ! MWwaterLUT data type definitions
  ! ----------------------------------
  !:tdoc+:
  TYPE :: MWwaterLUT_type
    ! Allocation indicator
    LOGICAL :: Is_Allocated = .FALSE.
    ! Release and version information
    INTEGER(Long) :: Release = MWwaterLUT_RELEASE
    INTEGER(Long) :: Version = MWwaterLUT_VERSION
    ! Dimensions
    INTEGER(Long) :: n_Angles       = 0   ! I1 dimension
    INTEGER(Long) :: n_Frequencies  = 0   ! I2 dimension
    INTEGER(Long) :: n_Temperatures = 0   ! I3 dimension
    INTEGER(Long) :: n_Wind_Speeds  = 0   ! I4 dimension
    ! Dimensional vectors
    REAL(Double),  ALLOCATABLE :: Angle(:)        ! I1
    REAL(Double),  ALLOCATABLE :: Frequency(:)    ! I2
    REAL(Double),  ALLOCATABLE :: Temperature(:)  ! I3
    REAL(Double),  ALLOCATABLE :: Wind_Speed(:)   ! I4
    ! Large-scale correction emissivity data
    REAL(Double),  ALLOCATABLE :: ev(:,:,:,:)  ! I1 x I2 x I3 x I4
    REAL(Double),  ALLOCATABLE :: eh(:,:,:,:)  ! I1 x I2 x I3 x I4
  END TYPE MWwaterLUT_type
  !:tdoc-:

CONTAINS

!################################################################################
!################################################################################
!##                                                                            ##
!##                           ## PUBLIC PROCEDURES ##                          ##
!##                                                                            ##
!################################################################################
!################################################################################

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       MWwaterLUT_Associated
!
! PURPOSE:
!       Pure function to test the status of the allocatable components
!       of the MWwaterLUT structure.
!
! CALLING SEQUENCE:
!       Status = MWwaterLUT_Associated( MWwaterLUT )
!
! OBJECTS:
!       MWwaterLUT:    Structure which is to have its member's
!                      status tested.
!                      UNITS:      N/A
!                      TYPE:       MWwaterLUT_type
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!       Status:        The return value is a logical value indicating the
!                      status of the components.
!                       .TRUE.  - if ANY of the MWwaterLUT allocatable members
!                                 are in use.
!                       .FALSE. - if ALL of the MWwaterLUT allocatable members
!                                 are not in use.
!                      UNITS:      N/A
!                      TYPE:       LOGICAL
!                      DIMENSION:  Same as input
!
!:sdoc-:
!--------------------------------------------------------------------------------

  PURE FUNCTION MWwaterLUT_Associated( self ) RESULT( Status )
    TYPE(MWwaterLUT_type), INTENT(IN) :: self
    LOGICAL :: Status
    Status = self%Is_Allocated
  END FUNCTION MWwaterLUT_Associated

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       MWwaterLUT_Destroy
!
! PURPOSE:
!       Pure subroutine to re-initialize MWwaterLUT objects.
!
! CALLING SEQUENCE:
!       CALL MWwaterLUT_Destroy( MWwaterLUT )
!
! OBJECTS:
!       MWwaterLUT:    Re-initialized MWwaterLUT structure.
!                      UNITS:      N/A
!                      TYPE:       MWwaterLUT_type
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: INTENT(OUT)
!
!:sdoc-:
!--------------------------------------------------------------------------------

  PURE SUBROUTINE MWwaterLUT_Destroy( self )
    TYPE(MWwaterLUT_type), INTENT(OUT) :: self
    self%Is_Allocated = .FALSE.
    self%n_Angles       = 0
    self%n_Frequencies  = 0
    self%n_Temperatures = 0
    self%n_Wind_Speeds  = 0
  END SUBROUTINE MWwaterLUT_Destroy

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       MWwaterLUT_Create
!
! PURPOSE:
!       Pure subroutine to create an instance of an MWwaterLUT object.
!
! CALLING SEQUENCE:
!       CALL MWwaterLUT_Create( MWwaterLUT    , &
!                               n_Angles      , &
!                               n_Frequencies , &
!                               n_Temperatures, &
!                               n_Wind_Speeds   )
!
! OBJECTS:
!       MWwaterLUT:         MWwaterLUT object structure.
!                           UNITS:      N/A
!                           TYPE:       MWwaterLUT_type
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(OUT)
!
! INPUTS:
!       n_Angles:           Number of zenith angles for which is are data.
!                           Must be > 0.
!                           UNITS:      N/A
!                           TYPE:       INTEGER
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(IN)
!
!       n_Frequencies:      Number of spectral frequencies for which there are
!                           data.
!                           Must be > 0.
!                           UNITS:      N/A
!                           TYPE:       INTEGER
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(IN)
!
!       n_Temperatures:     Number of surface temperatures for which there are
!                           data.
!                           Must be > 0.
!                           UNITS:      N/A
!                           TYPE:       INTEGER
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(IN)
!
!       n_Wind_Speeds:      Number of surface wind speeds for which there are
!                           data.
!                           Must be > 0.
!                           UNITS:      N/A
!                           TYPE:       INTEGER
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(IN)
!
!:sdoc-:
!--------------------------------------------------------------------------------

  PURE SUBROUTINE MWwaterLUT_Create( &
    self          , &  ! Output
    n_Angles      , &  ! Input
    n_Frequencies , &  ! Input
    n_Temperatures, &  ! Input
    n_Wind_Speeds   )  ! Input
    ! Arguments
    TYPE(MWwaterLUT_type), INTENT(OUT) :: self
    INTEGER              , INTENT(IN)  :: n_Angles
    INTEGER              , INTENT(IN)  :: n_Frequencies
    INTEGER              , INTENT(IN)  :: n_Temperatures
    INTEGER              , INTENT(IN)  :: n_Wind_Speeds
    ! Local variables
    INTEGER :: alloc_stat

    ! Check input
    IF ( n_Angles       < 1 .OR. &
         n_Frequencies  < 1 .OR. &
         n_Temperatures < 1 .OR. &
         n_Wind_Speeds  < 1 ) RETURN

    ! Perform the allocation
    ALLOCATE( self%Angle( n_Angles ), &
              self%Frequency( n_Frequencies ), &
              self%Temperature( n_Temperatures ), &
              self%Wind_Speed( n_Wind_Speeds ), &
              self%ev( n_Angles, n_Frequencies, n_Temperatures, n_Wind_Speeds ), &
              self%eh( n_Angles, n_Frequencies, n_Temperatures, n_Wind_Speeds ), &
              STAT = alloc_stat )
    IF ( alloc_stat /= 0 ) RETURN

    ! Initialise
    ! ...Dimensions
    self%n_Angles       = n_Angles
    self%n_Frequencies  = n_Frequencies
    self%n_Temperatures = n_Temperatures
    self%n_Wind_Speeds  = n_Wind_Speeds
    ! ...Arrays
    self%Angle       = ZERO
    self%Frequency   = ZERO
    self%Temperature = ZERO
    self%Wind_Speed  = ZERO
    self%ev          = ZERO
    self%eh          = ZERO

    ! Set allocation indicator
    self%Is_Allocated = .TRUE.

  END SUBROUTINE MWwaterLUT_Create

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       MWwaterLUT_Inspect
!
! PURPOSE:
!       Subroutine to print the contents of a MWwaterLUT object to stdout.
!
! CALLING SEQUENCE:
!       CALL MWwaterLUT_Inspect( MWwaterLUT )
!
! OBJECTS:
!       MWwaterLUT:     MWwaterLUT object to display.
!                       UNITS:      N/A
!                       TYPE:       MWwaterLUT_type
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN)
!
!:sdoc-:
!--------------------------------------------------------------------------------

  SUBROUTINE MWwaterLUT_Inspect( self, pause )
    TYPE(MWwaterLUT_type), INTENT(IN) :: self
    LOGICAL,     OPTIONAL, INTENT(IN) :: pause
    LOGICAL :: wait
    INTEGER :: i2, i3, i4

    wait = .FALSE.
    IF ( PRESENT(pause) ) wait = pause

    WRITE(*,'(1x,"MWwaterLUT OBJECT")')
    ! Release/version info
    WRITE(*,'(3x,"Release.Version : ",i0,".",i0)') self%Release, self%Version
    ! Dimensions
    WRITE(*,'(3x,"n_Angles        : ",i0)') self%n_Angles
    WRITE(*,'(3x,"n_Frequencies   : ",i0)') self%n_Frequencies
    WRITE(*,'(3x,"n_Temperatures  : ",i0)') self%n_Temperatures
    WRITE(*,'(3x,"n_Wind_Speeds   : ",i0)') self%n_Wind_Speeds
    IF ( .NOT. MWwaterLUT_Associated(self) ) RETURN
    ! Dimension arrays
    WRITE(*,'(3x,"Angle :")')
    WRITE(*,'(5(1x,es13.6,:))') self%Angle
    WRITE(*,'(3x,"Frequency :")')
    WRITE(*,'(5(1x,es13.6,:))') self%Frequency
    WRITE(*,'(3x,"Temperature :")')
    WRITE(*,'(5(1x,es13.6,:))') self%Temperature
    WRITE(*,'(3x,"Wind_Speed :")')
    WRITE(*,'(5(1x,es13.6,:))') self%Wind_Speed

    ! Emissivity arrays
    WRITE(*,'(/3x,"Emissivity(vertical polarisation) :")')
    IF ( wait ) THEN
      WRITE(*,FMT='(/1x,"Paused. Press <ENTER> to continue...")',ADVANCE='NO')
      READ(*,*)
    END IF

    DO i4 = 1, self%n_Wind_Speeds
      WRITE(*,'(5x,"WIND_SPEED  :",es13.6)') self%Wind_Speed(i4)
      DO i3 = 1, self%n_Temperatures
        WRITE(*,'(5x,"TEMPERATURE :",es13.6)') self%Temperature(i3)
        DO i2 = 1, self%n_Frequencies
          WRITE(*,'(5x,"FREQUENCY   :",es13.6)') self%Frequency(i2)
          WRITE(*,'(5(1x,es13.6,:))') self%ev(:,i2,i3,i4)
        END DO
      END DO
    END DO

    WRITE(*,'(/3x,"Emissivity(horizontal polarisation) :")')
    IF ( wait ) THEN
      WRITE(*,FMT='(/1x,"Paused. Press <ENTER> to continue...")',ADVANCE='NO')
      READ(*,*)
    END IF

    DO i4 = 1, self%n_Wind_Speeds
      WRITE(*,'(5x,"WIND_SPEED  :",es13.6)') self%Wind_Speed(i4)
      DO i3 = 1, self%n_Temperatures
        WRITE(*,'(5x,"TEMPERATURE :",es13.6)') self%Temperature(i3)
        DO i2 = 1, self%n_Frequencies
          WRITE(*,'(5x,"FREQUENCY   :",es13.6)') self%Frequency(i2)
          WRITE(*,'(5(1x,es13.6,:))') self%eh(:,i2,i3,i4)
        END DO
      END DO
    END DO
  END SUBROUTINE MWwaterLUT_Inspect

!----------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       MWwaterLUT_ValidRelease
!
! PURPOSE:
!       Function to check the MWwaterLUT Release value.
!
! CALLING SEQUENCE:
!       IsValid = MWwaterLUT_ValidRelease( MWwaterLUT )
!
! INPUTS:
!       MWwaterLUT:     MWwaterLUT object for which the Release component
!                       is to be checked.
!                       UNITS:      N/A
!                       TYPE:       MWwaterLUT_type
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!       IsValid:        Logical value defining the release validity.
!                       UNITS:      N/A
!                       TYPE:       LOGICAL
!                       DIMENSION:  Scalar
!
!:sdoc-:
!----------------------------------------------------------------------------------

  FUNCTION MWwaterLUT_ValidRelease( self ) RESULT( IsValid )
    ! Arguments
    TYPE(MWwaterLUT_type), INTENT(IN) :: self
    ! Function result
    LOGICAL :: IsValid
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'MWwaterLUT_ValidRelease'
    ! Local variables
    CHARACTER(ML) :: msg

    ! Set up
    IsValid = .TRUE.

    ! Check release is not too old
    IF ( self%Release < MWWATERLUT_RELEASE ) THEN
      IsValid = .FALSE.
      WRITE( msg,'("An MWwaterLUT data update is needed. ", &
                  &"MWwaterLUT release is ",i0,". Valid release is ",i0,"." )' ) &
                  self%Release, MWWATERLUT_RELEASE
      CALL Display_Message( ROUTINE_NAME, msg, INFORMATION ); RETURN
    END IF

    ! Check release is not too new
    IF ( self%Release > MWWATERLUT_RELEASE ) THEN
      IsValid = .FALSE.
      WRITE( msg,'("An MWwaterLUT software update is needed. ", &
                  &"MWwaterLUT release is ",i0,". Valid release is ",i0,"." )' ) &
                  self%Release, MWWATERLUT_RELEASE
      CALL Display_Message( ROUTINE_NAME, msg, INFORMATION ); RETURN
    END IF

  END FUNCTION MWwaterLUT_ValidRelease

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       MWwaterLUT_Info
!
! PURPOSE:
!       Subroutine to return a string containing version and dimension
!       information about a MWwaterLUT object.
!
! CALLING SEQUENCE:
!       CALL MWwaterLUT_Info( MWwaterLUT, Info )
!
! OBJECTS:
!       MWwaterLUT:      MWwaterLUT object about which info is required.
!                        UNITS:      N/A
!                        TYPE:       MWwaterLUT_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Info:            String containing version and dimension information
!                        about the MWwaterLUT object.
!                        UNITS:      N/A
!                        TYPE:       CHARACTER(*)
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(OUT)
!
!:sdoc-:
!--------------------------------------------------------------------------------

  SUBROUTINE MWwaterLUT_Info( self, Info )
    ! Arguments
    TYPE(MWwaterLUT_type), INTENT(IN)  :: self
    CHARACTER(*)         , INTENT(OUT) :: Info
    ! Parameters
    INTEGER, PARAMETER :: CARRIAGE_RETURN = 13
    INTEGER, PARAMETER :: LINEFEED = 10
    ! Local variables
    CHARACTER(2000) :: Long_String

    ! Write the required data to the local string
    WRITE( Long_String, &
           '(a,1x,"MWwaterLUT RELEASE.VERSION: ",i2,".",i2.2,a,3x, &
           &"N_ANGLES=",i0,2x,&
           &"N_FREQUENCIES=",i0,2x,&
           &"N_TEMPERATURES=",i0,2x,&
           &"N_WIND_SPEEDS=",i0 )' ) &
           ACHAR(CARRIAGE_RETURN)//ACHAR(LINEFEED), &
           self%Release, self%Version, &
           ACHAR(CARRIAGE_RETURN)//ACHAR(LINEFEED), &
           self%n_Angles      , &
           self%n_Frequencies , &
           self%n_Temperatures, &
           self%n_Wind_Speeds

    ! Trim the output based on the
    ! dummy argument string length
    Info = Long_String(1:MIN(LEN(Info), LEN_TRIM(Long_String)))

  END SUBROUTINE MWwaterLUT_Info

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       MWwaterLUT_DefineVersion
!
! PURPOSE:
!       Subroutine to return the module version information.
!
! CALLING SEQUENCE:
!       CALL MWwaterLUT_DefineVersion( Id )
!
! OUTPUTS:
!       Id:    Character string containing the version Id information
!              for the module.
!              UNITS:      N/A
!              TYPE:       CHARACTER(*)
!              DIMENSION:  Scalar
!              ATTRIBUTES: INTENT(OUT)
!
!:sdoc-:
!--------------------------------------------------------------------------------

  SUBROUTINE MWwaterLUT_DefineVersion( Id )
    CHARACTER(*), INTENT(OUT) :: Id
    Id = MODULE_VERSION_ID
  END SUBROUTINE MWwaterLUT_DefineVersion

!------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       MWwaterLUT_InquireFile
!
! PURPOSE:
!       Function to inquire MWwaterLUT object files.
!
! CALLING SEQUENCE:
!       Error_Status = MWwaterLUT_InquireFile( &
!                        Filename                       , &
!                        n_Angles       = n_Angles      , &
!                        n_Frequencies  = n_Frequencies , &
!                        n_Temperatures = n_Temperatures, &
!                        n_Wind_Speeds  = n_Wind_Speeds , &
!                        Release        = Release       , &
!                        Version        = Version       , &
!                        Title          = Title         , &
!                        History        = History       , &
!                        Comment        = Comment         )
!
! INPUTS:
!       Filename:           Character string specifying the name of the
!                           data file to inquire.
!                           UNITS:      N/A
!                           TYPE:       CHARACTER(*)
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(IN)
!
! OPTIONAL OUTPUTS:
!       n_Angles:           Number of zenith angles for which is are data.
!                           UNITS:      N/A
!                           TYPE:       INTEGER
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       n_Frequencies:      Number of spectral frequencies for which there are
!                           data.
!                           UNITS:      N/A
!                           TYPE:       INTEGER
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       n_Temperatures:     Number of surface temperatures for which there are
!                           data.
!                           UNITS:      N/A
!                           TYPE:       INTEGER
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       n_Wind_Speeds:      Number of surface wind speeds for which there are
!                           data.
!                           UNITS:      N/A
!                           TYPE:       INTEGER
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       Release:            The data/file release number. Used to check
!                           for data/software mismatch.
!                           UNITS:      N/A
!                           TYPE:       INTEGER
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       Version:            The data/file version number. Used for
!                           purposes only in identifying the dataset for
!                           a particular release.
!                           UNITS:      N/A
!                           TYPE:       INTEGER
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       Title:              Character string containing a succinct description
!                           of what is in the dataset.
!                           UNITS:      N/A
!                           TYPE:       CHARACTER(*)
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       History:            Character string containing dataset creation
!                           history.
!                           UNITS:      N/A
!                           TYPE:       CHARACTER(*)
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       Comment:            Character string containing any comments about
!                           the dataset.
!                           UNITS:      N/A
!                           TYPE:       CHARACTER(*)
!                           DIMENSION:  Scalar
!                           ATTRIBUTES: INTENT(OUT), OPTIONAL
!
! FUNCTION RESULT:
!       Error_Status:       The return value is an integer defining the error
!                           status. The error codes are defined in the
!                           Message_Handler module.
!                           If == SUCCESS the file inquire was successful
!                              == FAILURE an unrecoverable error occurred.
!                           UNITS:      N/A
!                           TYPE:       INTEGER
!                           DIMENSION:  Scalar
!
!:sdoc-:
!------------------------------------------------------------------------------

  FUNCTION MWwaterLUT_InquireFile( &
    Filename      , &  ! Input
    n_Angles      , &  ! Optional output
    n_Frequencies , &  ! Optional output
    n_Temperatures, &  ! Optional output
    n_Wind_Speeds , &  ! Optional output
    Release       , &  ! Optional output
    Version       , &  ! Optional output
    Title         , &  ! Optional output
    History       , &  ! Optional output
    Comment       ) &  ! Optional output
  RESULT( err_stat )
    ! Arguments
    CHARACTER(*),           INTENT(IN)  :: Filename
    INTEGER     , OPTIONAL, INTENT(OUT) :: n_Angles
    INTEGER     , OPTIONAL, INTENT(OUT) :: n_Frequencies
    INTEGER     , OPTIONAL, INTENT(OUT) :: n_Temperatures
    INTEGER     , OPTIONAL, INTENT(OUT) :: n_Wind_Speeds
    INTEGER     , OPTIONAL, INTENT(OUT) :: Release
    INTEGER     , OPTIONAL, INTENT(OUT) :: Version
    CHARACTER(*), OPTIONAL, INTENT(OUT) :: Title
    CHARACTER(*), OPTIONAL, INTENT(OUT) :: History
    CHARACTER(*), OPTIONAL, INTENT(OUT) :: Comment
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'MWwaterLUT_InquireFile'
    ! Function variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: io_msg
    INTEGER :: io_stat
    INTEGER :: fid
    TYPE(MWwaterLUT_type) :: MWwaterLUT

    ! Setup
    err_stat = SUCCESS
    ! ...Check that the file exists
    IF ( .NOT. File_Exists( Filename ) ) THEN
      msg = 'File '//TRIM(Filename)//' not found.'
      CALL Inquire_Cleanup(); RETURN
    END IF

    ! Open the file
    err_stat = Open_Binary_File( Filename, fid )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error opening '//TRIM(Filename)
      CALL Inquire_Cleanup(); RETURN
    END IF

    ! Read the release and version
    READ( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      MWwaterLUT%Release, &
      MWwaterLUT%Version
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading Release/Version - '//TRIM(io_msg)
      CALL Inquire_Cleanup(); RETURN
    END IF
    IF ( .NOT. MWwaterLUT_ValidRelease( MWwaterLUT ) ) THEN
      msg = 'MWwaterLUT Release check failed.'
      CALL Inquire_Cleanup(); RETURN
    END IF

    ! Read the dimensions
    READ( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      MWwaterLUT%n_Angles      , &
      MWwaterLUT%n_Frequencies , &
      MWwaterLUT%n_Temperatures, &
      MWwaterLUT%n_Wind_Speeds
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading dimension values from '//TRIM(Filename)//' - '//TRIM(io_msg)
      CALL Inquire_Cleanup(); RETURN
    END IF

    ! Read the global attributes
    err_stat = ReadGAtts_Binary_File( &
                 fid, &
                 Title   = Title  , &
                 History = History, &
                 Comment = Comment  )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error reading global attributes'
      CALL Inquire_Cleanup(); RETURN
    END IF

    ! Close the file
    CLOSE( fid, IOSTAT=io_stat, IOMSG=io_msg )
    IF ( io_stat /= 0 ) THEN
      msg = 'Error closing '//TRIM(Filename)//' - '//TRIM(io_msg)
      CALL Inquire_Cleanup(); RETURN
    END IF

    ! Assign the return arguments
    IF ( PRESENT(n_Angles      ) ) n_Angles       = MWwaterLUT%n_Angles
    IF ( PRESENT(n_Frequencies ) ) n_Frequencies  = MWwaterLUT%n_Frequencies
    IF ( PRESENT(n_Temperatures) ) n_Temperatures = MWwaterLUT%n_Temperatures
    IF ( PRESENT(n_Wind_Speeds ) ) n_Wind_Speeds  = MWwaterLUT%n_Wind_Speeds
    IF ( PRESENT(Release       ) ) Release        = MWwaterLUT%Release
    IF ( PRESENT(Version       ) ) Version        = MWwaterLUT%Version

  CONTAINS

    SUBROUTINE Inquire_CleanUp()
      ! Close file if necessary
      IF ( File_Open(fid) ) THEN
        CLOSE( fid, IOSTAT=io_stat, IOMSG=io_msg )
        IF ( io_stat /= 0 ) &
          msg = TRIM(msg)//'; Error closing input file during error cleanup - '//TRIM(io_msg)
      END IF
      ! Set error status and print error message
      err_stat = FAILURE
      CALL Display_Message( ROUTINE_NAME, msg, err_stat )
    END SUBROUTINE Inquire_CleanUp

  END FUNCTION MWwaterLUT_InquireFile

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       MWwaterLUT_ReadFile
!
! PURPOSE:
!       Function to read MWwaterLUT object files.
!
! CALLING SEQUENCE:
!       Error_Status = MWwaterLUT_ReadFile( &
!                        MWwaterLUT         , &
!                        Filename           , &
!                        No_Close = No_Close, &
!                        Quiet    = Quiet   , &
!                        Title    = Title   , &
!                        History  = History , &
!                        Comment  = Comment   )
!
! OBJECTS:
!       MWwaterLUT:     MWwaterLUT object containing the data read from file.
!                       UNITS:      N/A
!                       TYPE:       MWwaterLUT_type
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(OUT)
!
! INPUTS:
!       Filename:       Character string specifying the name of a
!                       MWwaterLUT data file to read.
!                       UNITS:      N/A
!                       TYPE:       CHARACTER(*)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN)
!
! OPTIONAL INPUTS:
!       No_Close:       Set this logical argument to *NOT* close the datafile
!                       upon exiting this routine. This option is required if
!                       the MWwaterLUT data is embedded within another file.
!                       If == .FALSE., File is closed upon function exit [DEFAULT].
!                          == .TRUE.,  File is NOT closed upon function exit
!                       If not specified, default is .FALSE.
!                       UNITS:      N/A
!                       TYPE:       LOGICAL
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN), OPTIONAL
!
!       Quiet:          Set this logical argument to suppress INFORMATION
!                       messages being printed to stdout
!                       If == .FALSE., INFORMATION messages are OUTPUT [DEFAULT].
!                          == .TRUE.,  INFORMATION messages are SUPPRESSED.
!                       If not specified, default is .FALSE.
!                       UNITS:      N/A
!                       TYPE:       LOGICAL
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN), OPTIONAL
!
! OPTIONAL OUTPUTS:
!       Title:          Character string containing a succinct description
!                       of what is in the dataset.
!                       UNITS:      N/A
!                       TYPE:       CHARACTER(*)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       History:        Character string containing dataset creation
!                       history.
!                       UNITS:      N/A
!                       TYPE:       CHARACTER(*)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       Comment:        Character string containing any comments about
!                       the dataset.
!                       UNITS:      N/A
!                       TYPE:       CHARACTER(*)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(OUT), OPTIONAL
!
! FUNCTION RESULT:
!       Error_Status:   The return value is an integer defining the error status.
!                       The error codes are defined in the Message_Handler module.
!                       If == SUCCESS, the file read was successful
!                          == FAILURE, an unrecoverable error occurred.
!                       UNITS:      N/A
!                       TYPE:       INTEGER
!                       DIMENSION:  Scalar
!
!:sdoc-:
!------------------------------------------------------------------------------

  FUNCTION MWwaterLUT_ReadFile( &
    MWwaterLUT, &  ! Output
    Filename  , &  ! Input
    No_Close  , &  ! Optional input
    Quiet     , &  ! Optional input
    Title     , &  ! Optional output
    History   , &  ! Optional output
    Comment   , &  ! Optional output
    Debug     ) &  ! Optional input (Debug output control)
  RESULT( err_stat )
    ! Arguments
    TYPE(MWwaterLUT_type) , INTENT(OUT) :: MWwaterLUT
    CHARACTER(*),           INTENT(IN)  :: Filename
    LOGICAL     , OPTIONAL, INTENT(IN)  :: No_Close
    LOGICAL     , OPTIONAL, INTENT(IN)  :: Quiet
    CHARACTER(*), OPTIONAL, INTENT(OUT) :: Title
    CHARACTER(*), OPTIONAL, INTENT(OUT) :: History
    CHARACTER(*), OPTIONAL, INTENT(OUT) :: Comment
    LOGICAL     , OPTIONAL, INTENT(IN)  :: Debug
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'MWwaterLUT_ReadFile'
    ! Function variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: io_msg
    LOGICAL :: close_file
    LOGICAL :: noisy
    INTEGER :: io_stat
    INTEGER :: fid
    TYPE(MWwaterLUT_type) :: dummy

    ! Setup
    err_stat = SUCCESS
    ! ...Check No_Close argument
    close_file = .TRUE.
    IF ( PRESENT(No_Close) ) close_file = .NOT. No_Close
    ! ...Check Quiet argument
    noisy = .TRUE.
    IF ( PRESENT(Quiet) ) noisy = .NOT. Quiet
    ! ...Override Quiet settings if debug set.
    IF ( PRESENT(Debug) ) THEN
      IF ( Debug ) noisy = .TRUE.
    END IF

    ! Check if the file is open.
    IF ( File_Open( Filename ) ) THEN
      ! ...Inquire for the logical unit number
      INQUIRE( FILE=Filename, NUMBER=fid )
      ! ...Ensure it's valid
      IF ( fid < 0 ) THEN
        msg = 'Error inquiring '//TRIM(Filename)//' for its FileID'
        CALL Read_CleanUp(); RETURN
      END IF
    ELSE
      ! ...Open the file if it exists
      IF ( File_Exists( Filename ) ) THEN
        err_stat = Open_Binary_File( Filename, fid )
        IF ( err_Stat /= SUCCESS ) THEN
          msg = 'Error opening '//TRIM(Filename)
          CALL Read_CleanUp(); RETURN
        END IF
      ELSE
        msg = 'File '//TRIM(Filename)//' not found.'
        CALL Read_CleanUp(); RETURN
      END IF
    END IF

    ! Read and check the release and version
    READ( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      dummy%Release, &
      dummy%Version
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading Release/Version - '//TRIM(io_msg)
      CALL Read_Cleanup(); RETURN
    END IF
    IF ( .NOT. MWwaterLUT_ValidRelease( dummy ) ) THEN
      msg = 'MWwaterLUT Release check failed.'
      CALL Read_Cleanup(); RETURN
    END IF

    ! Read the dimensions
    READ( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      dummy%n_Angles      , &
      dummy%n_Frequencies , &
      dummy%n_Temperatures, &
      dummy%n_Wind_Speeds
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading data dimensions - '//TRIM(io_msg)
      CALL Read_Cleanup(); RETURN
    END IF
    ! ...Allocate the object
    CALL MWwaterLUT_Create( &
           MWwaterLUT       , &
           dummy%n_Angles      , &
           dummy%n_Frequencies , &
           dummy%n_Temperatures, &
           dummy%n_Wind_Speeds   )
    IF ( .NOT. MWwaterLUT_Associated( MWwaterLUT ) ) THEN
      msg = 'MWwaterLUT object allocation failed.'
      CALL Read_Cleanup(); RETURN
    END IF
    ! ...Explicitly assign the version number
    MWwaterLUT%Version = dummy%Version

    ! Read the global attributes
    err_stat = ReadGAtts_Binary_File( &
                 fid, &
                 Title   = Title  , &
                 History = History, &
                 Comment = Comment  )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error reading global attributes'
      CALL Read_Cleanup(); RETURN
    END IF

    ! Read the coefficient data
    ! ...Read the dimensional vectors
    READ( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      MWwaterLUT%Angle      , &
      MWwaterLUT%Frequency  , &
      MWwaterLUT%Temperature, &
      MWwaterLUT%Wind_Speed
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading dimensional vectors - '//TRIM(io_msg)
      CALL Read_Cleanup(); RETURN
    END IF
    ! ...Read the emissivity data
    READ( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      MWwaterLUT%ev, &
      MWwaterLUT%eh
    IF ( io_stat /= 0 ) THEN
      msg = 'Error reading emissivity data - '//TRIM(io_msg)
      CALL Read_Cleanup(); RETURN
    END IF

    ! Close the file
    IF ( close_file ) THEN
      CLOSE( fid, IOSTAT=io_stat, IOMSG=io_msg )
      IF ( io_stat /= 0 ) THEN
        msg = 'Error closing '//TRIM(Filename)//' - '//TRIM(io_msg)
        CALL Read_Cleanup(); RETURN
      END IF
    END IF

    ! Output an info message
    IF ( noisy ) THEN
      CALL MWwaterLUT_Info( MWwaterLUT, msg )
      CALL Display_Message( ROUTINE_NAME, 'FILE: '//TRIM(Filename)//'; '//TRIM(msg), INFORMATION )
    END IF

  CONTAINS

    SUBROUTINE Read_CleanUp()
      IF ( File_Open(Filename) ) THEN
        CLOSE( fid, IOSTAT=io_stat, IOMSG=io_msg )
        IF ( io_stat /= 0 ) &
          msg = TRIM(msg)//'; Error closing input file during error cleanup - '//TRIM(io_msg)
      END IF
      CALL MWwaterLUT_Destroy( MWwaterLUT )
      err_stat = FAILURE
      CALL Display_Message( ROUTINE_NAME, msg, err_stat )
    END SUBROUTINE Read_CleanUp

  END FUNCTION MWwaterLUT_ReadFile

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       MWwaterLUT_WriteFile
!
! PURPOSE:
!       Function to write MWwaterLUT object files.
!
! CALLING SEQUENCE:
!       Error_Status = MWwaterLUT_WriteFile( &
!                        MWwaterLUT         , &
!                        Filename           , &
!                        No_Close = No_Close, &
!                        Quiet    = Quiet     )
!
! OBJECTS:
!       MWwaterLUT:     Object containing the data to write to file.
!                       UNITS:      N/A
!                       TYPE:       MWwaterLUT_type
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN)
!
! INPUTS:
!       Filename:       Character string specifying the name of a
!                       MWwaterLUT format data file to write.
!                       UNITS:      N/A
!                       TYPE:       CHARACTER(*)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN)
!
! OPTIONAL INPUTS:
!       No_Close:       Set this logical argument to *NOT* close the datafile
!                       upon exiting this routine. This option is required if
!                       the MWwaterLUT data is to be embedded within another file.
!                       If == .FALSE., File is closed upon function exit [DEFAULT].
!                          == .TRUE.,  File is NOT closed upon function exit
!                       If not specified, default is .FALSE.
!                       UNITS:      N/A
!                       TYPE:       LOGICAL
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN), OPTIONAL
!
!       Quiet:          Set this logical argument to suppress INFORMATION
!                       messages being printed to stdout
!                       If == .FALSE., INFORMATION messages are OUTPUT [DEFAULT].
!                          == .TRUE.,  INFORMATION messages are SUPPRESSED.
!                       If not specified, default is .FALSE.
!                       UNITS:      N/A
!                       TYPE:       LOGICAL
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN), OPTIONAL
!
!       Title:          Character string containing a succinct description
!                       of what is in the dataset.
!                       UNITS:      N/A
!                       TYPE:       CHARACTER(*)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN), OPTIONAL
!
!       History:        Character string containing dataset creation
!                       history.
!                       UNITS:      N/A
!                       TYPE:       CHARACTER(*)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN), OPTIONAL
!
!       Comment:        Character string containing any comments about
!                       the dataset.
!                       UNITS:      N/A
!                       TYPE:       CHARACTER(*)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN), OPTIONAL
!
! FUNCTION RESULT:
!       Error_Status:   The return value is an integer defining the error status.
!                       The error codes are defined in the Message_Handler module.
!                       If == SUCCESS, the file write was successful
!                          == FAILURE, an unrecoverable error occurred.
!                       UNITS:      N/A
!                       TYPE:       INTEGER
!                       DIMENSION:  Scalar
!
!:sdoc-:
!------------------------------------------------------------------------------

  FUNCTION MWwaterLUT_WriteFile( &
    MWwaterLUT, &  ! Input
    Filename  , &  ! Input
    No_Close  , &  ! Optional input
    Quiet     , &  ! Optional input
    Title     , &  ! Optional input
    History   , &  ! Optional input
    Comment   , &  ! Optional input
    Debug     ) &  ! Optional input (Debug output control)
  RESULT( err_stat )
    ! Arguments
    TYPE(MWwaterLUT_type) , INTENT(IN) :: MWwaterLUT
    CHARACTER(*),           INTENT(IN) :: Filename
    LOGICAL     , OPTIONAL, INTENT(IN) :: No_Close
    LOGICAL     , OPTIONAL, INTENT(IN) :: Quiet
    CHARACTER(*), OPTIONAL, INTENT(IN) :: Title
    CHARACTER(*), OPTIONAL, INTENT(IN) :: History
    CHARACTER(*), OPTIONAL, INTENT(IN) :: Comment
    LOGICAL     , OPTIONAL, INTENT(IN) :: Debug
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'MWwaterLUT_WriteFile'
    ! Function variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: io_msg
    LOGICAL :: close_file
    LOGICAL :: noisy
    INTEGER :: io_stat
    INTEGER :: fid

    ! Setup
    err_stat = SUCCESS
    ! ...Check No_Close argument
    close_file = .TRUE.
    IF ( PRESENT(No_Close) ) close_file = .NOT. No_Close
    ! ...Check Quiet argument
    noisy = .TRUE.
    IF ( PRESENT(Quiet) ) noisy = .NOT. Quiet
    ! ...Override Quiet settings if debug set.
    IF ( PRESENT(Debug) ) THEN
      IF ( Debug ) noisy = .TRUE.
    END IF
    ! ...Check there is data to write
    IF ( .NOT. MWwaterLUT_Associated( MWwaterLUT ) ) THEN
      msg = 'MWwaterLUT object is empty.'
      CALL Write_Cleanup(); RETURN
    END IF

    ! Check if the file is open.
    IF ( File_Open( FileName ) ) THEN
      ! ...Inquire for the logical unit number
      INQUIRE( FILE=Filename, NUMBER=fid )
      ! ...Ensure it's valid
      IF ( fid < 0 ) THEN
        msg = 'Error inquiring '//TRIM(Filename)//' for its FileID'
        CALL Write_CleanUp(); RETURN
      END IF
    ELSE
      ! ...Open the file for output
      err_stat = Open_Binary_File( Filename, fid, For_Output=.TRUE. )
      IF ( err_Stat /= SUCCESS ) THEN
        msg = 'Error opening '//TRIM(Filename)
        CALL Write_CleanUp(); RETURN
      END IF
    END IF

    ! Write the release and version
    WRITE( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      MWwaterLUT%Release, &
      MWwaterLUT%Version
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing Release/Version - '//TRIM(io_msg)
      CALL Write_Cleanup(); RETURN
    END IF

    ! Write the dimensions
    WRITE( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      MWwaterLUT%n_Angles      , &
      MWwaterLUT%n_Frequencies , &
      MWwaterLUT%n_Temperatures, &
      MWwaterLUT%n_Wind_Speeds
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing data dimensions - '//TRIM(io_msg)
      CALL Write_Cleanup(); RETURN
    END IF

    ! Write the global attributes
    err_stat = WriteGAtts_Binary_File( &
                 fid, &
                 Write_Module = MODULE_VERSION_ID, &
                 Title        = Title  , &
                 History      = History, &
                 Comment      = Comment  )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error writing global attributes'
      CALL Write_Cleanup(); RETURN
    END IF

    ! Write the coefficient data
    ! ...Write the dimensional vectors
    WRITE( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      MWwaterLUT%Angle      , &
      MWwaterLUT%Frequency  , &
      MWwaterLUT%Temperature, &
      MWwaterLUT%Wind_Speed
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing dimensional vectors - '//TRIM(io_msg)
      CALL Write_Cleanup(); RETURN
    END IF
    ! ...Write the emissivity data
    WRITE( fid, IOSTAT=io_stat, IOMSG=io_msg ) &
      MWwaterLUT%ev, &
      MWwaterLUT%eh
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing the emissivity data - '//TRIM(io_msg)
      CALL Write_Cleanup(); RETURN
    END IF

    ! Close the file
    IF ( close_file ) THEN
      CLOSE( fid, IOSTAT=io_stat, IOMSG=io_msg )
      IF ( io_stat /= 0 ) THEN
        msg = 'Error closing '//TRIM(Filename)//' - '//TRIM(io_msg)
        CALL Write_Cleanup(); RETURN
      END IF
    END IF

    ! Output an info message
     IF ( noisy ) THEN
       CALL MWwaterLUT_Info( MWwaterLUT, msg )
       CALL Display_Message( ROUTINE_NAME, 'FILE: '//TRIM(Filename)//'; '//TRIM(msg), INFORMATION )
     END IF

   CONTAINS

     SUBROUTINE Write_Cleanup()
       IF ( File_Open(Filename) ) THEN
         CLOSE( fid, STATUS=WRITE_ERROR_STATUS, IOSTAT=io_stat, IOMSG=io_msg )
         IF ( io_stat /= 0 ) &
           msg = TRIM(msg)//'; Error closing output file during error cleanup - '//TRIM(io_msg)
       END IF
       err_stat = FAILURE
       CALL Display_Message( ROUTINE_NAME, msg, err_stat )
     END SUBROUTINE Write_Cleanup

  END FUNCTION MWwaterLUT_WriteFile

!################################################################################
!################################################################################
!##                                                                            ##
!##                          ## PRIVATE PROCEDURES ##                          ##
!##                                                                            ##
!################################################################################
!################################################################################

!--------------------------------------------------------------------------------
!
! NAME:
!       MWwaterLUT_Equal
!
! PURPOSE:
!       Pure function to test the equality of two MWwaterLUT objects.
!       Used in OPERATOR(==) interface block.
!
! CALLING SEQUENCE:
!       is_equal = MWwaterLUT_Equal( x, y )
!
!         or
!
!       IF ( x == y ) THEN
!         ...
!       END IF
!
! OBJECTS:
!       x, y:          Two MWwaterLUT objects to be compared.
!                      UNITS:      N/A
!                      TYPE:       MWwaterLUT_type
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!       is_equal:      Logical value indicating whether the inputs are equal.
!                      UNITS:      N/A
!                      TYPE:       LOGICAL
!                      DIMENSION:  Same as inputs.
!
!--------------------------------------------------------------------------------

  PURE FUNCTION MWwaterLUT_Equal( x, y ) RESULT( is_equal )
    TYPE(MWwaterLUT_type), INTENT(IN) :: x, y
    LOGICAL :: is_equal

    ! Set up
    is_equal = .FALSE.

    ! Check the object association status
    IF ( (.NOT. MWwaterLUT_Associated(x)) .OR. &
         (.NOT. MWwaterLUT_Associated(y))      ) RETURN

   ! Check contents
    ! ...Release/version info
    IF ( (x%Release /= y%Release) .OR. &
         (x%Version /= y%Version) ) RETURN
    ! ...Dimensions
    IF ( (x%n_Angles       /= y%n_Angles       ) .OR. &
         (x%n_Frequencies  /= y%n_Frequencies  ) .OR. &
         (x%n_Temperatures /= y%n_Temperatures ) .OR. &
         (x%n_Wind_Speeds  /= y%n_Wind_Speeds  ) ) RETURN
    ! ...Arrays
    IF ( ALL(x%Angle       .EqualTo. y%Angle       ) .AND. &
         ALL(x%Frequency   .EqualTo. y%Frequency   ) .AND. &
         ALL(x%Temperature .EqualTo. y%Temperature ) .AND. &
         ALL(x%Wind_Speed  .EqualTo. y%Wind_Speed  ) .AND. &
         ALL(x%ev          .EqualTo. y%ev          ) .AND. &
         ALL(x%eh          .EqualTo. y%eh          ) ) &
      is_equal = .TRUE.

  END FUNCTION MWwaterLUT_Equal

END MODULE MWwaterLUT_Define

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb

max_forks_repo_name: teb-model/wrf-teb


text:
*DECK R9LGMC
      FUNCTION R9LGMC (X)
C***BEGIN PROLOGUE  R9LGMC
C***SUBSIDIARY
C***PURPOSE  Compute the log Gamma correction factor so that
C            LOG(GAMMA(X)) = LOG(SQRT(2*PI)) + (X-.5)*LOG(X) - X
C            + R9LGMC(X).
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7E
C***TYPE      SINGLE PRECISION (R9LGMC-S, D9LGMC-D, C9LGMC-C)
C***KEYWORDS  COMPLETE GAMMA FUNCTION, CORRECTION TERM, FNLIB,
C             LOG GAMMA, LOGARITHM, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C Compute the log gamma correction factor for X .GE. 10.0 so that
C  LOG (GAMMA(X)) = LOG(SQRT(2*PI)) + (X-.5)*LOG(X) - X + R9LGMC(X)
C
C Series for ALGM       on the interval  0.          to  1.00000D-02
C                                        with weighted error   3.40E-16
C                                         log weighted error  15.47
C                               significant figures required  14.39
C                                    decimal places required  15.86
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  CSEVL, INITS, R1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770801  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900720  Routine changed from user-callable to subsidiary.  (WRB)
C***END PROLOGUE  R9LGMC
      DIMENSION ALGMCS(6)
      LOGICAL FIRST
      SAVE ALGMCS, NALGM, XBIG, XMAX, FIRST
      DATA ALGMCS( 1) /    .1666389480 45186E0 /
      DATA ALGMCS( 2) /   -.0000138494 817606E0 /
      DATA ALGMCS( 3) /    .0000000098 108256E0 /
      DATA ALGMCS( 4) /   -.0000000000 180912E0 /
      DATA ALGMCS( 5) /    .0000000000 000622E0 /
      DATA ALGMCS( 6) /   -.0000000000 000003E0 /
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  R9LGMC
      IF (FIRST) THEN
         NALGM = INITS (ALGMCS, 6, R1MACH(3))
         XBIG = 1.0/SQRT(R1MACH(3))
         XMAX = EXP (MIN(LOG(R1MACH(2)/12.0), -LOG(12.0*R1MACH(1))) )
      ENDIF
      FIRST = .FALSE.
C
      IF (X .LT. 10.0) CALL XERMSG ('SLATEC', 'R9LGMC',
     +   'X MUST BE GE 10', 1, 2)
      IF (X.GE.XMAX) GO TO 20
C
      R9LGMC = 1.0/(12.0*X)
      IF (X.LT.XBIG) R9LGMC = CSEVL (2.0*(10./X)**2-1., ALGMCS, NALGM)/X
      RETURN
C
 20   R9LGMC = 0.0
      CALL XERMSG ('SLATEC', 'R9LGMC', 'X SO BIG R9LGMC UNDERFLOWS', 2,
     +   1)
      RETURN
C
      END