text:
!  ------------------------------------------------------------------------------------------------------
!
!    Program HAMS for the diffraction and radiation of waves
!    by 3D structures.
!
!  License:
!
!    This routine is part of HAMS.
!
!    HAMS is a free software framework: you can redistribute it and/or modify it
!    under the terms of the Apache License, Version 2.0 (the "License"); you may
!    not use these subroutines except in compliance with the License. The software
!    is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND.
!
!    You should have received a copy of the Apache License, Version 2.0, along with
!    HAMS. If not, see <http://www.apache.org/licenses/LICENSE-2.0>.
!
!  Code Original Author:
!
!    Yingyi Liu
!
!  ------------------------------------------------------------------------------------------------------
MODULE BodyIntgr_irr

   USE HAMS_mod
   USE Body_mod
   USE Const_mod
   USE WaveDyn_mod
   USE PanelMesh_mod
   USE Inerfs_mod
   USE PatcVelct
   USE Potentials_mod
   USE omp_lib
   IMPLICIT NONE

      ! ..... Public Subroutines ...................................................................................................

   PUBLIC :: BODINT_LEFT_IRR
   PUBLIC :: RBC_IRR
   PUBLIC :: DBC_IRR

CONTAINS

!   --------------------------------------------------------
!     Integration inside an element, which is not singular
!   --------------------------------------------------------

      SUBROUTINE BODINT_LEFT_IRR(IS,IEL,JEL,TINDP,IRR,FLAG)
      IMPLICIT NONE

      INTEGER,INTENT(IN):: IS,IEL,JEL,IRR,FLAG
      COMPLEX*16,INTENT(OUT):: TINDP(4)

      INTEGER  IP
      COMPLEX*16  GRN(4)
      REAL*8::   RKN(4),ENV(6),TNV(6),EAR,XEC(3)

      IF (IRR.EQ.1) THEN

        RKN(:)=RKBN(IEL,JEL,IS,:)
        GRN(:)=CGRN(IEL,JEL,IS,:)

        EAR=DS(JEL)
        XEC(:)=XYZ_P(JEL,:)
        ENV(:)=DXYZ_P(JEL,:)
        TNV(:)=DXYZ_P(IEL,:)

      ELSEIF (IRR.EQ.3) THEN

        RKN(:)=PKBN(IEL,JEL,IS,:)
        GRN(:)=DGRN(IEL,JEL,IS,:)

        EAR=DS(JEL)
        XEC(:)=XYZ_P(JEL,:)
        ENV(:)=DXYZ_P(JEL,:)
        TNV(:)=iDXYZ_P(IEL,:)

      ENDIF

      IF (FLAG.EQ.1) THEN

        TINDP(IS)=(RKN(2)+GRN(2)*EAR)*ENV(1)   &
                 +(RKN(3)+GRN(3)*EAR)*ENV(2)   &
                 +(RKN(4)+GRN(4)*EAR)*ENV(3)

      ELSE

        TINDP(IS)=(GRN(2)*ENV(1)+GRN(3)*ENV(2)+GRN(4)*ENV(3))*EAR

      ENDIF

      RETURN
      END SUBROUTINE BODINT_LEFT_IRR

!   --------------------------------------------------------
!     Assemble the radiated right-side boundary conditions
!   --------------------------------------------------------

      SUBROUTINE RBC_IRR(IS,IEL,JEL,TINRD,IRR,FLAG)
      IMPLICIT NONE

      INTEGER,INTENT(IN):: IS,IEL,JEL,IRR,FLAG
      COMPLEX*16,INTENT(OUT):: TINRD(4,6,4)

      INTEGER  IP
      REAL*8  XQ,YQ,ZQ,EAR
      REAL*8  RKN(4),ENV(6),TNV(6),XEC(3)
      COMPLEX*16  GRN(4),DUM

      IF (IRR.EQ.1) THEN

        RKN(:)=RKBN(IEL,JEL,IS,:)
        GRN(:)=CGRN(IEL,JEL,IS,:)

        EAR=DS(JEL)
        XEC(:)=XYZ_P(JEL,:)
        ENV(:)=DXYZ_P(JEL,:)
        TNV(:)=DXYZ_P(IEL,:)

      ELSEIF (IRR.EQ.3) THEN

        RKN(:)=PKBN(IEL,JEL,IS,:)
        GRN(:)=DGRN(IEL,JEL,IS,:)

        EAR=DS(JEL)
        XEC(:)=XYZ_P(JEL,:)
        ENV(:)=DXYZ_P(JEL,:)
        TNV(:)=iDXYZ_P(IEL,:)

      ENDIF

!  Above----
!  Field point changes, source element keeps the same as in the first quadrant

      IF (FLAG.EQ.1) THEN
        DUM=RKN(1)+GRN(1)*EAR
      ELSE
        DUM=GRN(1)*EAR
      ENDIF

      DO 200 IP=1, NSYS
       IF (ISX.EQ.1.AND.ISY.EQ.0) THEN
        TINRD(IS,1,IP)=TINRD(IS,1,IP)+DUM*SY(IS,IP)*ENV(1)
        TINRD(IS,2,IP)=TINRD(IS,2,IP)+DUM*SX(IS,IP)*ENV(2)
        TINRD(IS,3,IP)=TINRD(IS,3,IP)+DUM*ENV(3)
        TINRD(IS,4,IP)=TINRD(IS,4,IP)+DUM*SX(IS,IP)*ENV(4)
        TINRD(IS,5,IP)=TINRD(IS,5,IP)+DUM*SY(IS,IP)*ENV(5)
        TINRD(IS,6,IP)=TINRD(IS,6,IP)+DUM*SY(IS,IP)*SX(IS,IP)*ENV(6)
       ELSE
        TINRD(IS,1,IP)=TINRD(IS,1,IP)+DUM*SX(IS,IP)*ENV(1)
        TINRD(IS,2,IP)=TINRD(IS,2,IP)+DUM*SY(IS,IP)*ENV(2)
        TINRD(IS,3,IP)=TINRD(IS,3,IP)+DUM*ENV(3)
        TINRD(IS,4,IP)=TINRD(IS,4,IP)+DUM*SY(IS,IP)*ENV(4)
        TINRD(IS,5,IP)=TINRD(IS,5,IP)+DUM*SX(IS,IP)*ENV(5)
        TINRD(IS,6,IP)=TINRD(IS,6,IP)+DUM*SX(IS,IP)*SY(IS,IP)*ENV(6)
       ENDIF

200   CONTINUE

      RETURN
      END SUBROUTINE RBC_IRR

!   --------------------------------------------------------
!     Assemble the diffracted right-side boundary conditions
!   --------------------------------------------------------

      SUBROUTINE DBC_IRR(IS,IEL,JEL,TINRD,IRR,FLAG)
      IMPLICIT NONE

      INTEGER,INTENT(IN):: IS,IEL,JEL,IRR,FLAG
      COMPLEX*16,INTENT(OUT):: TINRD(4,4)

      INTEGER  IP
      REAL*8  XQ,YQ,ZQ,EAR
      REAL*8  RKN(4),ENV(6),TNV(6),XEC(3)

      COMPLEX*16  GRN(4),DUM(2)
      COMPLEX*16  DPOX,DPOY,DPOZ

      IF (IRR.EQ.1) THEN

        RKN(:)=RKBN(IEL,JEL,IS,:)
        GRN(:)=CGRN(IEL,JEL,IS,:)

        EAR=DS(JEL)
        XEC(:)=XYZ_P(JEL,:)
        ENV(:)=DXYZ_P(JEL,:)
        TNV(:)=DXYZ_P(IEL,:)

      ELSEIF (IRR.EQ.3) THEN

        RKN(:)=PKBN(IEL,JEL,IS,:)
        GRN(:)=DGRN(IEL,JEL,IS,:)

        EAR=DS(JEL)
        XEC(:)=XYZ_P(JEL,:)
        ENV(:)=DXYZ_P(JEL,:)
        TNV(:)=iDXYZ_P(IEL,:)

      ENDIF

!  Above----
!  Field point changes, source element keeps the same as in the first quadrant

      IF (FLAG.EQ.1) THEN
        DUM(1)=RKN(1)+GRN(1)*EAR
      ELSE
        DUM(1)=GRN(1)*EAR
      ENDIF

      DO 200 IP=1, NSYS

       IF (ISX.EQ.1.AND.ISY.EQ.0) THEN
        XQ=SY(IS,IP)*XEC(1)
        YQ=SX(IS,IP)*XEC(2)
        ZQ=           XEC(3)
        CALL DINP(XQ,YQ,ZQ,XW(1),XW(2),BETA,DPOX,DPOY,DPOZ)
        DUM(2)=SY(IS,IP)*DPOX*ENV(1)+SX(IS,IP)*DPOY*ENV(2)+DPOZ*ENV(3)
        TINRD(IS,IP)=TINRD(IS,IP)-DUM(1)*DUM(2)
       ELSE
        XQ=SX(IS,IP)*XEC(1)
        YQ=SY(IS,IP)*XEC(2)
        ZQ=          XEC(3)
        CALL DINP(XQ,YQ,ZQ,XW(1),XW(2),BETA,DPOX,DPOY,DPOZ)
        DUM(2)=SX(IS,IP)*DPOX*ENV(1)+SY(IS,IP)*DPOY*ENV(2)+DPOZ*ENV(3)
        TINRD(IS,IP)=TINRD(IS,IP)-DUM(1)*DUM(2)
       ENDIF

200   CONTINUE

      RETURN
      END SUBROUTINE DBC_IRR
!-------------------------------------------------------------------------------
END MODULE BodyIntgr_irr
!*******************************************************************************

max_stars_repo_name: shousner/pyHAMS

max_issues_repo_name: shousner/pyHAMS

max_forks_repo_name: shousner/pyHAMS


text:
! a simple code to check the analytic Jacobian via numerical
! differencing

subroutine test_jacobian() bind(C)

  use network
  use eos_module
  use burner_module
  use actual_burner_module
  use actual_rhs_module
  use burn_type_module

  use amrex_fort_module, only : rt => amrex_real
  implicit none

  type (burn_t) :: state, statep, statem

  type (eos_t) :: eos_state

  character (len=32) :: probin_file
  integer :: probin_pass(32)
  integer :: i, j

  real(rt)        , parameter :: delta = 1.e-6_rt
  real(rt)        , parameter :: SMALL = 1.e-12_rt
  real(rt)         :: num_jac

  character(len=16) :: namei,namej

  probin_file = "probin"
  do i = 1, len(trim(probin_file))
     probin_pass(i) = ichar(probin_file(i:i))
  enddo

  call runtime_init(probin_pass(1:len(trim(probin_file))), len(trim(probin_file)))

  call network_init()
  call burner_init()
  call eos_init()

  ! Set up state

  state % rho   = 2.0e7_rt
  state % T     = 8.0e9_rt

  state % xn(:) = ONE / nspec

  call burn_to_eos(state, eos_state)
  call normalize_abundances(eos_state)
  call eos(eos_input_rt, eos_state)
  call eos_to_burn(eos_state, state)

  state % self_heat = .true.

  ! Evaluate the analytical Jacobian. Note that we
  ! need to call f_rhs first because that will fill
  ! the state with the rates that the Jacobian needs.

  call actual_rhs(state)
  call actual_jac(state)

888 format(a,"-derivatives that don't match:")
999 format(5x, "df(",a,")/dy(",a,")", g18.10, g18.10, g18.10)

  ! Now evaluate a numerical estimate of the Jacobian
  ! using the RHS.

  do j = 1, neqs

     statep = state
     statem = state

     if (j <= nspec) then
        statep % xn(j) = (ONE + delta) * state % xn(j)
     else if (j == net_itemp) then
        statep % T     = (ONE + delta) * state % T
     else if (j == net_ienuc) then
        statep % e     = (ONE + delta) * state % e
     endif

     call actual_rhs(statep)

     if (j <= nspec) then
        statem % xn(j) = (ONE - delta) * state % xn(j)
     else if (j == net_itemp) then
        statem % T     = (ONE - delta) * state % T
     else if (j == net_ienuc) then
        statem % e     = (ONE - delta) * state % e
     endif

     call actual_rhs(statem)

     if (j <= nspec) then
        namej = short_spec_names(j)
     else if (j==net_ienuc) then
        namej = "e"
     else if (j==net_itemp) then
        namej = "T"
     endif

     write(*,888) trim(namej)

     do i = 2, neqs

        if (j <= nspec) then
           num_jac = (statep % ydot(i) - statem % ydot(i))/(statep % xn(j) - statem % xn(j) + SMALL)
        else if (j == net_itemp) then
           num_jac = (statep % ydot(i) - statem % ydot(i))/(statep % T - statem % T + SMALL)
        else if (j == net_ienuc) then
           num_jac = (statep % ydot(i) - statem % ydot(i))/(statep % e - statem % e + SMALL)
        endif

        if (i <= nspec) then
           namei = short_spec_names(i)
        else if (i==net_ienuc) then
           namei = "e"
        else if (i==net_itemp) then
           namei = "T"
        endif

        ! only dump the ones that don't match
        if (num_jac /= state % jac(i,j)) then
           if (num_jac /= ZERO) then
              write (*,999) trim(namei), &
                   trim(namej), num_jac, state % jac(i,j), (num_jac-state % jac(i,j))/num_jac
           else
              write (*,999) trim(namei), &
                   trim(namej), num_jac, state % jac(i,j)
           endif
        endif

     enddo
  enddo

end subroutine test_jacobian

max_stars_repo_name: doreenfan/Microphysics

max_issues_repo_name: doreenfan/Microphysics

max_forks_repo_name: doreenfan/Microphysics


text:
      SUBROUTINE DV7CPY(P, Y, X)
C
C  ***  SET Y = X, WHERE X AND Y ARE P-VECTORS  ***
C
      INTEGER P
      DOUBLE PRECISION X(P), Y(P)
C
      INTEGER I
C
      DO 10 I = 1, P
 10      Y(I) = X(I)
      RETURN
      END

ext: f

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
The Pavilion at ARC is part of the Activities and Recreation Center.  It was formerly a separate building known as the Rec Hall, which was completed in 1977. Its construction was funded entirely by private and student support. Prior to the Pavilion, there was the now demolished Activities and Recreation Center.

The Pavilion was and still is a decent place to work out,  but the university added a large extension with a store, study lounge, and all, and renamed it the ARC.

The Pavilion hosts both mens and womens Aggie basketball games.

Entertainment events are often held in the Pavilion. WWE World Wrestling Federation/Entertainment house shows have been held there, but sadly rival wrestling promotion NWATNA has never come to Davis.   Turbulence used to be held in the Pavilion annually; its replacement The Buzz use to be held at the Pavilion and ARC.

Notable guests include Bill Clinton, the 42nd President of the United States, who appeared January 15, 2008 to an overflow crowd of 7800 people (according to the press) to campaign for his wife, presidential hopeful Hillary Clinton.

Some weight lifting P.E. classes are still held in the Rec Hall in the old equipment rooms, although the ARC weight rooms are nicer. Some students still use the Rec Halls running track because it is cooler (in temperature), and its also easier to keep track of your miles (7.5 laps at the Rec Hall vs. 7.65 at the ARC. Laps are also easier to run in the Pavilion.  There are 4 corners, which turn in the same direction, was opposed to the ARCs inandout/alongtheside track.   Other hardcore exercisers still use the Rec Hall because it has the option of  stair laps, in which one does a lap around the track via running up and down the rows of stairs. Some of the studios on the second floor are still used by various sports and clubs such and the Ping Pong Club. It is also nice to go to the Rec Hall because it is quieter and not as much of a meat market.

The main floor of the Rec Hall was replaced with a multiuse floor in the summer of 2008, which helped with the nonsporting events which are becoming more common uses.  There is a portable wood floor, so that the arena can still be used for basketball, volleyball, etc.  The old wooden bleachers were also replaced.  This construction dispersed the high school cheerleading camps that summer to places like Freeborn Hall for maximum jailbait potential.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! { dg-do compile }
! PR 79312 - assigning a logical value to a real
! is invalid.
! Test case by John Harper.
program emptyarray5
  implicit none
  real a(0)
  a = [logical::] ! { dg-error "Can't convert LOGICAL" }
  print *,size(a)
end program emptyarray5


text:
      SUBROUTINE MB04QB( TRANC, TRAND, TRANQ, STOREV, STOREW, M, N, K,
     $                   V, LDV, W, LDW, C, LDC, D, LDD, CS, TAU, DWORK,
     $                   LDWORK, INFO )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To overwrite general real m-by-n matrices C and D, or their
C     transposes, with
C
C               [ op(C) ]
C         Q  *  [       ]   if TRANQ = 'N', or
C               [ op(D) ]
C
C          T    [ op(C) ]
C         Q  *  [       ]   if TRANQ = 'T',
C               [ op(D) ]
C
C     where Q is defined as the product of symplectic reflectors and
C     Givens rotations,
C
C         Q = diag( H(1),H(1) ) G(1) diag( F(1),F(1) )
C             diag( H(2),H(2) ) G(2) diag( F(2),F(2) )
C                               ....
C             diag( H(k),H(k) ) G(k) diag( F(k),F(k) ).
C
C     Blocked version.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     TRANC   CHARACTER*1
C             Specifies the form of op( C ) as follows:
C             = 'N':  op( C ) = C;
C             = 'T':  op( C ) = C';
C             = 'C':  op( C ) = C'.
C
C     TRAND   CHARACTER*1
C             Specifies the form of op( D ) as follows:
C             = 'N':  op( D ) = D;
C             = 'T':  op( D ) = D';
C             = 'C':  op( D ) = D'.
C
C     TRANQ   CHARACTER*1
C             = 'N':  apply Q;
C             = 'T':  apply Q'.
C
C     STOREV  CHARACTER*1
C             Specifies how the vectors which define the concatenated
C             Householder reflectors contained in V are stored:
C             = 'C':  columnwise;
C             = 'R':  rowwise.
C
C     STOREW  CHARACTER*1
C             Specifies how the vectors which define the concatenated
C             Householder reflectors contained in W are stored:
C             = 'C':  columnwise;
C             = 'R':  rowwise.
C
C     Input/Output Parameters
C
C     M       (input) INTEGER
C             The number of rows of the matrices op(C) and op(D).
C             M >= 0.
C
C     N       (input) INTEGER
C             The number of columns of the matrices op(C) and op(D).
C             N >= 0.
C
C     K       (input) INTEGER
C             The number of elementary reflectors whose product defines
C             the matrix Q.  M >= K >= 0.
C
C     V       (input) DOUBLE PRECISION array, dimension
C                     (LDV,K) if STOREV = 'C',
C                     (LDV,M) if STOREV = 'R'
C             On entry with STOREV = 'C', the leading M-by-K part of
C             this array must contain in its columns the vectors which
C             define the elementary reflectors F(i).
C             On entry with STOREV = 'R', the leading K-by-M part of
C             this array must contain in its rows the vectors which
C             define the elementary reflectors F(i).
C
C     LDV     INTEGER
C             The leading dimension of the array V.
C             LDV >= MAX(1,M),  if STOREV = 'C';
C             LDV >= MAX(1,K),  if STOREV = 'R'.
C
C     W       (input) DOUBLE PRECISION array, dimension
C                     (LDW,K) if STOREW = 'C',
C                     (LDW,M) if STOREW = 'R'
C             On entry with STOREW = 'C', the leading M-by-K part of
C             this array must contain in its columns the vectors which
C             define the elementary reflectors H(i).
C             On entry with STOREW = 'R', the leading K-by-M part of
C             this array must contain in its rows the vectors which
C             define the elementary reflectors H(i).
C
C     LDW     INTEGER
C             The leading dimension of the array W.
C             LDW >= MAX(1,M),  if STOREW = 'C';
C             LDW >= MAX(1,K),  if STOREW = 'R'.
C
C     C       (input/output) DOUBLE PRECISION array, dimension
C                     (LDC,N) if TRANC = 'N',
C                     (LDC,M) if TRANC = 'T' or TRANC = 'C'
C             On entry with TRANC = 'N', the leading M-by-N part of
C             this array must contain the matrix C.
C             On entry with TRANC = 'C' or TRANC = 'T', the leading
C             N-by-M part of this array must contain the transpose of
C             the matrix C.
C             On exit with TRANC = 'N', the leading M-by-N part of
C             this array contains the updated matrix C.
C             On exit with TRANC = 'C' or TRANC = 'T', the leading
C             N-by-M part of this array contains the transpose of the
C             updated matrix C.
C
C     LDC     INTEGER
C             The leading dimension of the array C.
C             LDC >= MAX(1,M),  if TRANC = 'N';
C             LDC >= MAX(1,N),  if TRANC = 'T' or TRANC = 'C'.
C
C     D       (input/output) DOUBLE PRECISION array, dimension
C                     (LDD,N) if TRAND = 'N',
C                     (LDD,M) if TRAND = 'T' or TRAND = 'C'
C             On entry with TRAND = 'N', the leading M-by-N part of
C             this array must contain the matrix D.
C             On entry with TRAND = 'C' or TRAND = 'T', the leading
C             N-by-M part of this array must contain the transpose of
C             the matrix D.
C             On exit with TRAND = 'N', the leading M-by-N part of
C             this array contains the updated matrix D.
C             On exit with TRAND = 'C' or TRAND = 'T', the leading
C             N-by-M part of this array contains the transpose of the
C             updated matrix D.
C
C     LDD     INTEGER
C             The leading dimension of the array D.
C             LDD >= MAX(1,M),  if TRAND = 'N';
C             LDD >= MAX(1,N),  if TRAND = 'T' or TRAND = 'C'.
C
C     CS      (input) DOUBLE PRECISION array, dimension (2*K)
C             On entry, the first 2*K elements of this array must
C             contain the cosines and sines of the symplectic Givens
C             rotations G(i).
C
C     TAU     (input) DOUBLE PRECISION array, dimension (K)
C             On entry, the first K elements of this array must
C             contain the scalar factors of the elementary reflectors
C             F(i).
C
C     Workspace
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if INFO = 0,  DWORK(1)  returns the optimal
C             value of LDWORK.
C             On exit, if  INFO = -20,  DWORK(1)  returns the minimum
C             value of LDWORK.
C
C     LDWORK  INTEGER
C             The length of the array DWORK.  LDWORK >= MAX(1,N).
C
C             If LDWORK = -1, then a workspace query is assumed;
C             the routine only calculates the optimal size of the
C             DWORK array, returns this value as the first entry of
C             the DWORK array, and no error message related to LDWORK
C             is issued by XERBLA.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value.
C
C     REFERENCES
C
C     [1] Kressner, D.
C         Block algorithms for orthogonal symplectic factorizations.
C         BIT, 43 (4), pp. 775-790, 2003.
C
C     CONTRIBUTORS
C
C     D. Kressner, Technical Univ. Berlin, Germany, and
C     P. Benner, Technical Univ. Chemnitz, Germany, December 2003.
C
C     REVISIONS
C
C     V. Sima, June 2008 (SLICOT version of the HAPACK routine DOSMSB).
C     V. Sima, Aug. 2011, Nov. 2011.
C
C     KEYWORDS
C
C     Elementary matrix operations, orthogonal symplectic matrix.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION  ONE
      PARAMETER         ( ONE = 1.0D+0 )
C     .. Scalar Arguments ..
      CHARACTER         STOREV, STOREW, TRANC, TRAND, TRANQ
      INTEGER           INFO, K, LDC, LDD, LDV, LDW, LDWORK, M, N
C     .. Array Arguments ..
      DOUBLE PRECISION  C(LDC,*), CS(*), D(LDD,*), DWORK(*), TAU(*),
     $                  V(LDV,*), W(LDW,*)
C     .. Local Scalars ..
      CHARACTER*1       SIDE
      LOGICAL           LCOLV, LCOLW, LQUERY, LTRC, LTRD, LTRQ
      INTEGER           I, IB, IC, ID, IERR, JC, JD, KI, KK, MINWRK, NB,
     $                  NBMIN, NX, PDRS, PDT, PDW, WRKOPT
C     .. External Functions ..
      INTEGER           UE01MD
      LOGICAL           LSAME
      EXTERNAL          LSAME, UE01MD
C     .. External Subroutines ..
      EXTERNAL          DORMQR, MB04QC, MB04QF, MB04QU, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC         DBLE, INT, MAX, MIN, SQRT
C
C     .. Executable Statements ..
C
C     Decode the scalar input parameters.
C
      INFO  = 0
      LCOLV = LSAME( STOREV, 'C' )
      LCOLW = LSAME( STOREW, 'C' )
      LTRC  = LSAME( TRANC, 'T' ) .OR. LSAME( TRANC, 'C' )
      LTRD  = LSAME( TRAND, 'T' ) .OR. LSAME( TRAND, 'C' )
      LTRQ  = LSAME( TRANQ, 'T' )
C
C     Check the scalar input parameters.
C
      IF ( .NOT.( LTRC .OR. LSAME( TRANC, 'N' ) ) ) THEN
         INFO = -1
      ELSE IF ( .NOT.( LTRD  .OR. LSAME( TRAND,  'N' ) ) ) THEN
         INFO = -2
      ELSE IF ( .NOT.( LTRQ  .OR. LSAME( TRANQ,  'N' ) ) ) THEN
         INFO = -3
      ELSE IF ( .NOT.( LCOLV .OR. LSAME( STOREV, 'R' ) ) ) THEN
         INFO = -4
      ELSE IF ( .NOT.( LCOLW .OR. LSAME( STOREW, 'R' ) ) ) THEN
         INFO = -5
      ELSE IF ( M.LT.0 ) THEN
         INFO = -6
      ELSE IF ( N.LT.0 ) THEN
         INFO = -7
      ELSE IF ( K.LT.0 .OR. K.GT.M ) THEN
         INFO = -8
      ELSE IF ( ( LCOLV .AND. LDV.LT.MAX( 1, M ) ) .OR.
     $     ( .NOT.LCOLV .AND. LDV.LT.MAX( 1, K ) ) ) THEN
         INFO = -10
      ELSE IF ( ( LCOLW .AND. LDW.LT.MAX( 1, M ) ) .OR.
     $     ( .NOT.LCOLW .AND. LDW.LT.MAX( 1, K ) ) ) THEN
         INFO = -12
      ELSE IF ( ( LTRC  .AND. LDC.LT.MAX( 1, N ) ) .OR.
     $     ( .NOT.LTRC  .AND. LDC.LT.MAX( 1, M ) ) ) THEN
         INFO = -14
      ELSE IF ( ( LTRD  .AND. LDD.LT.MAX( 1, N ) ) .OR.
     $     ( .NOT.LTRD  .AND. LDD.LT.MAX( 1, M ) ) ) THEN
         INFO = -16
      ELSE
         LQUERY = LDWORK.EQ.-1
         MINWRK = MAX( 1, N )
         IF ( LDWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
            DWORK(1) = DBLE( MINWRK )
            INFO = -20
         ELSE
            IF ( N.EQ.0 ) THEN
               WRKOPT = ONE
            ELSE
               I = MAX( 1, M, N )
               IF ( LTRC ) THEN
                  SIDE = 'R'
                  IC   = N
                  JC   = M
               ELSE
                  SIDE = 'L'
                  IC   = M
                  JC   = N
               END IF
               CALL DORMQR( SIDE, TRANC, IC, JC, K, DWORK, I, DWORK,
     $                      DWORK, I, DWORK, -1, INFO )
               WRKOPT = MAX( MINWRK, INT( DWORK(1) ) )
               NB     = MIN( INT( WRKOPT/N ), N )
               WRKOPT = MAX( WRKOPT, 9*N*NB + 15*NB*NB )
            END IF
            IF ( LQUERY ) THEN
               DWORK(1) = DBLE( WRKOPT )
               RETURN
            END IF
         END IF
      END IF
C
C     Return if there were illegal values.
C
      IF ( INFO.NE.0 ) THEN
         CALL XERBLA( 'MB04QB', -INFO )
         RETURN
      END IF
C
C     Quick return if possible.
C
      IF ( MIN( K, M, N ).EQ.0 ) THEN
         DWORK(1) = ONE
         RETURN
      END IF
C
      NBMIN = 2
      NX = 0
      IF ( NB.GT.1 .AND. NB.LT.K ) THEN
C
C        Determine when to cross over from blocked to unblocked code.
C
         NX = MAX( 0, UE01MD( 3, 'MB04QB', TRANC // TRAND // TRANQ, M,
     $                        N, K ) )
         IF ( NX.LT.K ) THEN
C
C           Determine if workspace is large enough for blocked code.
C
            IF ( LDWORK.LT.WRKOPT ) THEN
C
C              Not enough workspace to use optimal NB:  reduce NB and
C              determine the minimum value of NB.
C
               NB = INT( ( SQRT( DBLE( 81*N*N + 60*LDWORK ) )
     $                     - DBLE( 9*N ) ) / 30.0D0 )
               NBMIN = MAX( 2, UE01MD( 2, 'MB04QB', TRANC // TRAND //
     $                                 TRANQ, M, N, K ) )
            END IF
         END IF
      END IF
C
      PDRS = 1
      PDT  = PDRS + 6*NB*NB
      PDW  = PDT  + 9*NB*NB
      IC = 1
      JC = 1
      ID = 1
      JD = 1
C
      IF ( LTRQ ) THEN
C
C        Use blocked code initially.
C
         IF ( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
            DO 10 I = 1, K - NX, NB
               IB = MIN( K-I+1, NB )
C
C              Form the triangular factors of the symplectic block
C              reflector SH.
C
               CALL MB04QF( 'Forward', STOREV, STOREW, M-I+1, IB,
     $                      V(I,I), LDV, W(I,I), LDW, CS(2*I-1), TAU(I),
     $                      DWORK(PDRS), NB, DWORK(PDT), NB,
     $                      DWORK(PDW) )
C
C              Apply SH' to [ op(C)(i:m,:); op(D)(i:m,:) ] from the
C              left.
C
               IF ( LTRC ) THEN
                  JC = I
               ELSE
                  IC = I
               END IF
               IF ( LTRD ) THEN
                  JD = I
               ELSE
                  ID = I
               END IF
               CALL MB04QC( 'No Structure', TRANC, TRAND, TRANQ,
     $                      'Forward', STOREV, STOREW, M-I+1, N, IB,
     $                      V(I,I), LDV, W(I,I), LDW, DWORK(PDRS), NB,
     $                      DWORK(PDT), NB, C(IC,JC), LDC, D(ID,JD),
     $                      LDD, DWORK(PDW) )
   10       CONTINUE
         ELSE
            I = 1
         END IF
C
C        Use unblocked code to update last or only block.
C
         IF ( I.LE.K ) THEN
            IF ( LTRC ) THEN
               JC = I
            ELSE
               IC = I
            END IF
            IF ( LTRD ) THEN
               JD = I
            ELSE
               ID = I
            END IF
            CALL MB04QU( TRANC, TRAND, TRANQ, STOREV, STOREW, M-I+1, N,
     $                   K-I+1, V(I,I), LDV, W(I,I), LDW, C(IC,JC), LDC,
     $                   D(ID,JD), LDD, CS(2*I-1), TAU(I), DWORK,
     $                   LDWORK, IERR )
         END IF
      ELSE
         IF ( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
C
C           Use blocked code after the last block.
C           The first kk columns are handled by the block method.
C
            KI = ( ( K-NX-1 ) / NB )*NB
            KK = MIN( K, KI+NB )
         ELSE
            KK = 0
         END IF
C
C        Use unblocked code for the last or only block.
C
         IF ( KK.LT.K ) THEN
            IF ( LTRC ) THEN
               JC = KK + 1
            ELSE
               IC = KK + 1
            END IF
            IF ( LTRD ) THEN
               JD = KK + 1
            ELSE
               ID = KK + 1
            END IF
            CALL MB04QU( TRANC, TRAND, TRANQ, STOREV, STOREW, M-KK, N,
     $                   K-KK, V(KK+1,KK+1), LDV, W(KK+1,KK+1), LDW,
     $                   C(IC,JC), LDC, D(ID,JD), LDD, CS(2*KK+1),
     $                   TAU(KK+1), DWORK, LDWORK, IERR )
         END IF
C
C        Blocked code.
C
         IF ( KK.GT.0 ) THEN
            DO 20 I = KI + 1, 1, -NB
               IB = MIN( NB, K-I+1 )
C
C              Form the triangular factors of the symplectic block
C              reflector SH.
C
               CALL MB04QF( 'Forward', STOREV, STOREW, M-I+1, IB,
     $                      V(I,I), LDV, W(I,I), LDW, CS(2*I-1), TAU(I),
     $                      DWORK(PDRS), NB, DWORK(PDT), NB,
     $                      DWORK(PDW) )
C
C              Apply SH to [ op(C)(i:m,:); op(D)(i:m,:) ] from
C              the left.
C
               IF ( LTRC ) THEN
                  JC = I
               ELSE
                  IC = I
               END IF
               IF ( LTRD ) THEN
                  JD = I
               ELSE
                  ID = I
               END IF
               CALL MB04QC( 'No Structure', TRANC, TRAND, TRANQ,
     $                      'Forward', STOREV, STOREW, M-I+1, N, IB,
     $                      V(I,I), LDV, W(I,I), LDW, DWORK(PDRS), NB,
     $                      DWORK(PDT), NB, C(IC,JC), LDC, D(ID,JD),
     $                      LDD, DWORK(PDW) )
   20       CONTINUE
         END IF
      END IF
      DWORK(1) = DBLE( WRKOPT )
C
      RETURN
C *** Last line of MB04QB ***
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: RJHKnight/slicotr

max_forks_repo_name: RJHKnight/slicotr


text:
	SUBROUTINE LV_CCRD  ( ivcord, vcoord, iret )
C************************************************************************
C* LV_CCRD								*
C*									*
C* This subroutine translates a numeric value for IVCORD into its 	*
C* character value in VCOORD.						*
C*									*
C* LV_CCRD  ( IVCORD, VCOORD, IRET )					*
C*									*
C* Input parameters:							*
C*	IVCORD		INTEGER		Numeric vertical coordinate	*
C*									*
C* Output parameters:							*
C*	VCOORD		CHAR*		Vertical coordinate 		*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -3 = invalid coordinate	*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 4/89	From LV_CORD				*
C* K. Brill/NMC		 4/91	Added SGMA and DPTH			*
C* M. desJardins/NMC	10/91	Generalize vertical coordinate		*
C* F. J. Yen/NCEP 	 7/08	Added HYBL				*
C* M. James/Unidata 	 1/13	Added PVBL and PVAB			*
C************************************************************************
	CHARACTER*(*)	 vcoord
C*
	CHARACTER	v*1
C*
	PARAMETER	( NVERT= 9 )
	CHARACTER	vert (NVERT)*4
C*
	DATA 		vert  / 'NONE', 'PRES', 'THTA', 'HGHT',
     +				'SGMA', 'DPTH', 'HYBL',
     +				'PVBL', 'PVAB' /
C------------------------------------------------------------------------
	iret  = 0
C
C*	Translate known vertical coordinates or look for parameter name.
C
	ier = 0
C
C*	Check for numeric vertical coordinates.
C
	IF  ( ( ivcord .ge. 0 ) .and. ( ivcord .le. NVERT - 1 ) )  THEN
	    vcoord = vert ( ivcord + 1 )
C
C*	    Check for character name as vertical coordinate.  Check that
C*	    each character is an alphanumeric character.
C
	  ELSE IF  ( ivcord .gt. 100 )  THEN
	    CALL ST_ITOC  ( ivcord, 1, vcoord, ier )
	    IF  ( ier .eq. 0 )  THEN
		DO  i = 1, 4
		    v = vcoord (i:i)
		    IF  ( ( ( v .lt. 'A' ) .or. ( v .gt. 'Z' ) ) .and.
     +			  ( ( v .lt. '0' ) .or. ( v .gt. '9' ) ) )  THEN
			ier = -1
		    END IF
		END DO
	    END IF
	  ELSE
	    ier = -1
	END IF
C
C*	Check for errors.
C
	IF  ( ier .ne. 0 )  THEN
	    vcoord = ' '
	    iret   = -3
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: sgdecker/pygempak


text:
!======================================================================!
SUBROUTINE Calc1
    !----------------------------------------------------------------------!
    !   Calculate node coordinates inside the domain, block by block.      !
    !----------------------------------------------------------------------!
    !------------------------------[Modules]-------------------------------!
    USE allp_mod, only : huge
    USE all_mod
    USE gen_mod
    !----------------------------------------------------------------------!
    IMPLICIT NONE
    !------------------------------[Calling]-------------------------------!
    INTEGER :: IsLine
    !-----------------------------[Interface]------------------------------!
    INTERFACE
        LOGICAL FUNCTION Approx(A,B,tol)
            IMPLICIT NONE
            REAL          :: A,B
            REAL,OPTIONAL :: tol
        END FUNCTION Approx
    END INTERFACE
    !-------------------------------[Locals]-------------------------------!
    INTEGER :: b, bl, i, j, k, n, c, ig
    INTEGER :: l, l1, l2
    INTEGER :: is, js, ks, ie, je, ke, face
    INTEGER :: NI, NJ, NK, CI, CJ, CK
    INTEGER :: trans(3,2)
    !--------------------------------[CVS]---------------------------------!
    !  $Id: Calc1.f90,v 1.10 2002/10/30 16:29:21 niceno Exp $
    !  $Source: /home/muhamed/.CVSROOT/T-Rex/Generate/Calc1.f90,v $
    !======================================================================!

    do n=1,MAXN
        x(n)=HUGE
    end do

    !++++++++++++++++++++++++++++++++++++++++!
    !     calculate the node coordinates     !
    !++++++++++++++++++++++++++++++++++++++++!
    NN = 0                          ! Initialize n.o.n.
    NC = 0                          ! Initialize n.o.v.
    do b=1,Nbloc
        write(*,*) 'Generating block: ', b
        NI=BlkRes(b, 1)
        NJ=BlkRes(b, 2)
        NK=BlkRes(b, 3)

        !-( 1 )-!
        n = NN + ( 1-1)*NI*NJ + ( 1-1)*NI +  1
        x(n)=xp(BlkPnt(b, 1))
        y(n)=yp(BlkPnt(b, 1))
        z(n)=zp(BlkPnt(b, 1))

        !-( 2 )-!
        n = NN + ( 1-1)*NI*NJ + ( 1-1)*NI + NI
        x(n)=xp(BlkPnt(b, 2))
        y(n)=yp(BlkPnt(b, 2))
        z(n)=zp(BlkPnt(b, 2))

        !-( 3 )-!
        n = NN + ( 1-1)*NI*NJ + (NJ-1)*NI +  1
        x(n)=xp(BlkPnt(b, 3))
        y(n)=yp(BlkPnt(b, 3))
        z(n)=zp(BlkPnt(b, 3))

        !-( 4 )-!
        n = NN + ( 1-1)*NI*NJ + (NJ-1)*NI + NI
        x(n)=xp(BlkPnt(b, 4))
        y(n)=yp(BlkPnt(b, 4))
        z(n)=zp(BlkPnt(b, 4))

        !-( 5 )-!
        n = NN + (NK-1)*NI*NJ + ( 1-1)*NI +  1
        x(n)=xp(BlkPnt(b, 5))
        y(n)=yp(BlkPnt(b, 5))
        z(n)=zp(BlkPnt(b, 5))

        !-( 6 )-!
        n = NN + (NK-1)*NI*NJ + ( 1-1)*NI + NI
        x(n)=xp(BlkPnt(b, 6))
        y(n)=yp(BlkPnt(b, 6))
        z(n)=zp(BlkPnt(b, 6))

        !-( 7 )-!
        n = NN + (NK-1)*NI*NJ + (NJ-1)*NI +  1
        x(n)=xp(BlkPnt(b, 7))
        y(n)=yp(BlkPnt(b, 7))
        z(n)=zp(BlkPnt(b, 7))

        !-( 8 )-!
        n = NN + (NK-1)*NI*NJ + (NJ-1)*NI + NI
        x(n)=xp(BlkPnt(b, 8))
        y(n)=yp(BlkPnt(b, 8))
        z(n)=zp(BlkPnt(b, 8))

        !------------------------------!
        !       First on the lines     !
        !     defined point by point   !
        !------------------------------!
        do l=1, Nline

            bl = IsLine(LinPnt(l,1),LinPnt(l,2),b)

            if(bl == b) then

                do n=1,8
                    if( LinPnt(l,1) == BlkPnt(b,n)  ) l1=n
                    if( LinPnt(l,2) == BlkPnt(b,n)  ) l2=n
                end do

                !..... line is defined in the +i direction
                if     ( (l2-l1) == +1 ) then
                    trans(1,1)= 0
                    trans(1,2)=+1
                    trans(2,2)= 0
                    trans(3,2)= 0
                    if( (l1 == 1).or.(l1 == 5) ) then
                        trans(2,1)=1
                    else
                        trans(2,1)=BlkRes(b,2)
                    endif
                    if( (l1 == 1).or.(l1 == 3) ) then
                        trans(3,1)=1
                    else
                        trans(3,1)=BlkRes(b,3)
                    endif

                !..... line is defined in the -i direction
                else if( (l2-l1) == -1 ) then
                    trans(1,1)=BlkRes(b,1)+1   ! NI from the block + 1
                    trans(1,2)=-1
                    trans(2,2)= 0
                    trans(3,2)= 0
                    if( (l1 == 2).or.(l1 == 6) ) then
                        trans(2,1)=1
                    else
                        trans(2,1)=BlkRes(b,2)
                    endif
                    if( (l1 == 2).or.(l1 == 4) ) then
                        trans(3,1)=1
                    else
                        trans(3,1)=BlkRes(b,3)
                    endif

                !..... line is defined in the +j direction
                else if( (l2-l1) == +2 ) then
                    trans(2,1)= 0
                    trans(2,2)=+1
                    trans(1,2)= 0
                    trans(3,2)= 0
                    if( (l1 == 1).or.(l1 == 5) ) then
                        trans(1,1)=1
                    else
                        trans(1,1)=BlkRes(b,1)
                    endif
                    if( (l1 == 1).or.(l1 == 2) ) then
                        trans(3,1)=1
                    else
                        trans(3,1)=BlkRes(b,3)
                    endif

                !..... line is defined in the -j direction
                else if( (l2-l1) == -2 ) then
                    trans(2,1)=BlkRes(b,2)+1   ! NJ from the block + 1
                    trans(2,2)=-1
                    trans(1,2)= 0
                    trans(3,2)= 0
                    if( (l1 == 3).or.(l1 == 7) ) then
                        trans(1,1)=1
                    else
                        trans(1,1)=BlkRes(b,1)
                    endif
                    if( (l1 == 3).or.(l1 == 4) ) then
                        trans(3,1)=1
                    else
                        trans(3,1)=BlkRes(b,3)
                    endif

                !..... line is defined in the +k direction
                else if( (l2-l1) == +4 ) then
                    trans(3,1)= 0
                    trans(3,2)=+1
                    trans(1,2)= 0
                    trans(2,2)= 0
                    if( (l1 == 1).or.(l1 == 3) ) then
                        trans(1,1)=1
                    else
                        trans(1,1)=BlkRes(b,1)
                    end if
                    if( (l1 == 1).or.(l1 == 2) ) then
                        trans(2,1)=1
                    else
                        trans(2,1)=BlkRes(b,2)
                    endif

                !..... line is defined in the -k direction
                else if( (l2-l1) == -4 ) then
                    trans(3,1)=BlkRes(b,3) + 1  ! NK from the block + 1
                    trans(3,2)=-1
                    trans(1,2)= 0
                    trans(2,2)= 0
                    if( (l1 == 5).or.(l1 == 7) ) then
                        trans(1,1)=1
                    else
                        trans(1,1)=BlkRes(b,1)
                    endif
                    if( (l1 == 5).or.(l1 == 6) ) then
                        trans(2,1)=1
                    else
                        trans(2,1)=BlkRes(b,2)
                    endif

                endif ! l1-l2

                !----- linija je zadan tocka po tocka; stara fora
                if(LinWgt(l) ==  0.0) then
                    write(6,*) 'LINIJA: ', l
                    write(6,*) 'l1= ', l1
                    write(6,*) 'l2= ', l2
                    do ig=1,LinRes(l)
                        i=trans(1,1)+trans(1,2)*ig
                        j=trans(2,1)+trans(2,2)*ig
                        k=trans(3,1)+trans(3,2)*ig

                        n = NN + (k-1)*NI*NJ + (j-1)*NI + i
                        x(n) = xl(l,ig)
                        y(n) = yl(l,ig)
                        z(n) = zl(l,ig)
                    end do

                !----- linija je zadana tezinskim faktorom; nova fora
                else
                    is=trans(1,1)+trans(1,2)
                    js=trans(2,1)+trans(2,2)
                    ks=trans(3,1)+trans(3,2)
                    ie=trans(1,1)+trans(1,2)*LinRes(l)
                    je=trans(2,1)+trans(2,2)*LinRes(l)
                    ke=trans(3,1)+trans(3,2)*LinRes(l)
                    call linija( b, LinWgt(l), is, js, ks, ie, je, ke)
                endif

            endif ! if the block contains

        end do ! for the lines

        !------------------------!
        !       lines...         !
        !                        !
        !------------------------!
        do k=1,NK,NK-1
            do j=1,NJ,NJ-1
                call linija( b, BlkWgt(b,1), 1,j,k,NI,j,k)
            end do
        end do

        do k=1,NK,NK-1
            do i=1,NI,NI-1
                call linija( b, BlkWgt(b,2), i,1,k,i,NJ,k)
            end do
        end do

        do j=1,NJ,NJ-1
            do i=1,NI,NI-1
                call linija( b, BlkWgt(b,3), i,j,1,i,j,NK)
            end do
        end do

        !------------------------------------------------------------!
        !       surfaces...                                          !
        !                                                            !
        !   I think this is the propper way to calculate surfaces:   !
        !   it spans the lines in the direction of higher weigh      !
        !                                                            !
        !------------------------------------------------------------!

        !----------------------------------------------------------------------
        ! I (k=1)
        n = (b-1)*6 + 1   ! face index
        k=1
        if( .NOT. Approx(BlfaWt(n,1),1.0 ) ) then
            do j=1,NJ
                call linija( b, BlfaWt(n,1), 1,j,k,NI,j,k)
            end do
        else ! lines in the j direction
            do i=1,NI
                call linija( b, BlfaWt(n,2), i,1,k,i,NJ,k)
            end do
        endif

        ! VI (k=NK)
        n = (b-1)*6 + 6   ! face index
        k=NK
        if( .NOT. Approx(BlfaWt(n,1),1.0 ) ) then
            do j=1,NJ
                call linija( b, BlFaWt(n,1), 1,j,k,NI,j,k)
            end do
        else ! lines in the j direction
            do i=1,NI
                call linija( b, BlfaWt(n,2), i,1,k,i,NJ,k)
            end do
        endif
        !----------------------------------------------------------------------
        ! V (i=1)
        n = (b-1)*6 + 5   ! face index
        i=1
        if( .NOT. Approx(BlfaWt(n,3),1.0 ) ) then
            do j=1,NJ
                call linija( b, BlFaWt(n,3), i,j,1,i,j,NK)
            end do
        else ! lines in the j direction
            do k=1,NK
                call linija( b, BlFaWt(n,2), i,1,k,i,NJ,k)
            end do
        end if
        ! III (i=NI)
        n = (b-1)*6 + 3   ! face index
        i=NI
        if( .NOT. Approx(BlfaWt(n,3),1.0 ) ) then
            do j=1,NJ
                call linija( b, BlFaWt(n,3), i,j,1,i,j,NK)
            end do
        else ! lines in the j direction
            do k=1,NK
                call linija( b, BlFaWt(n,2), i,1,k,i,NJ,k)
            end do
        end if
        !----------------------------------------------------------------------
        ! II (j=1)
        n = (b-1)*6 + 2   ! face index
        j=1
        if( .NOT. Approx(BlfaWt(n,3),1.0 ) ) then
            do i=1,NI
                call linija( b, BlFaWt(n,3), i,j,1,i,j,NK)
            end do
        else ! lines in the i direction
            do k=1,NK
                call linija( b, BlFaWt(n,1), 1,j,k,NI,j,k)
            end do
        endif
        ! IV (j=NJ)
        n = (b-1)*6 + 4   ! face index
        j=NJ
        if( .NOT. Approx(BlfaWt(n,3),1.0 ) ) then
            do i=1,NI
                call linija( b, BlFaWt(n,3), i,j,1,i,j,NK)
            end do
        else ! lines in the i direction
            do k=1,NK
                call linija( b, BlFaWt(n,1), 1,j,k,NI,j,k)
            end do
        endif
        !----------------------------------------------------------------------

        !------------------------!
        !       volumes...       !
        !                        !
        !------------------------!
        if( .NOT. Approx( BlkWgt(b,3), 1.0 ) ) then
            do i=1,NI
                do j=1,NJ
                    call linija( b, BlkWgt(b,3), i,j,1,i,j,NK)
                end do
            end do
        else if( .NOT. Approx( BlkWgt(b,1), 1.0 ) ) then
            do k=1,NK
                do j=1,NJ
                    call linija( b, BlkWgt(b,1), 1,j,k,NI,j,k)
                end do
            end do
        else if( .NOT. Approx( BlkWgt(b,2), 1.0 ) ) then
            do k=1,NK
                do i=1,NI
                    call linija( b, BlkWgt(b,2), i,1,k,i,NJ,k)
                end do
            end do
        else

            do i=1,NI
                do j=1,NJ
                    do k=1,NK
                        n = NN+(k-1)*NI*NJ + (j-1)*NI + i
                        call Laplac(b, i, j, k, 0.333, 0.333, 0.334,           &
                            0.333, 0.333, 0.334,           &
                            0.333, 0.333, 0.334)
                    end do
                end do
            end do
        endif

        !--------------------------------------------!
        !     set the control volume nodes (CellN)  !
        !     and  the control volume neighbours     !
        !--------------------------------------------!
        CI = NI-1
        CJ = NJ-1
        CK = NK-1

        do k=1,CK
            do j=1,CJ
                do i=1,CI
                    c = NC + (k-1)*CI*CJ + (j-1)*CI + i ! cell
                    n = NN + (k-1)*NI*NJ + (j-1)*NI + i ! 1st node

                    !----- nodes
                    CellN(c,1)=n
                    CellN(c,2)=n+1
                    CellN(c,3)=n+NI
                    CellN(c,4)=n+NI+1
                    CellN(c,5)=CellN(c,1)+NI*NJ
                    CellN(c,6)=CellN(c,2)+NI*NJ
                    CellN(c,7)=CellN(c,3)+NI*NJ
                    CellN(c,8)=CellN(c,4)+NI*NJ

                    !----- neighbours
                    CellC(c,1)= c-CI*CJ
                    CellC(c,2)= c-CI
                    CellC(c,3)= c+1
                    CellC(c,4)= c+CI
                    CellC(c,5)= c-1
                    CellC(c,6)= c+CI*CJ

                    !----- this value (-1) is also the default boundary marker
                    if(i == 1)  CellC(c,5)=-1
                    if(i == CI) CellC(c,3)=-1
                    if(j == 1)  CellC(c,2)=-1
                    if(j == CJ) CellC(c,4)=-1
                    if(k == 1)  CellC(c,1)=-1
                    if(k == CK) CellC(c,6)=-1

                end do
            end do
        end do

        BlkRes(b,4) = NI*NJ*NK ! is this needed ???
        BlkRes(b,5) = NN       ! old number of nodes, for fuzion
        BlkRes(b,6) = NC       ! old number of volumes, for fuzion
        NN = NN + NI*NJ*NK
        NC = NC + CI*CJ*CK

    end do   ! through blocks

    !---------------------------------------------!
    !     Insertion of the boundary condition     !
    !           and materials information         !
    !---------------------------------------------!

    !----- initialize all the material markers to 1
    do c=1,NC
        material(c)=1
    end do

    do n=1,Nboun

        b  = abs(Bound(n,7))   ! block

        !---- block resolution
        CI=BlkRes(b,1)-1
        CJ=BlkRes(b,2)-1
        CK=BlkRes(b,3)-1

        !---- default values
        is=1
        ie=CI
        js=1
        je=CJ
        ks=1
        ke=CK

        !---- boundary conditions prescribed with mnemonics
        if(BndFac(n) == 'IMIN') then
            ie=1
            face = 5
        else if(BndFac(n) == 'IMAX') then
            is=CI
            face = 3
        else if(BndFac(n) == 'JMIN') then
            je=1
            face = 2
        else if(BndFac(n) == 'JMAX') then
            js=CJ
            face = 4
        else if(BndFac(n) == 'KMIN') then
            ke=1
            face = 1
        else if(BndFac(n) == 'KMAX') then
            ks=CK
            face = 6
        !---- boundary conditions (materials) prescribed explicitly
        !     (error prone and difficult, but might be usefull)
        else
            is = Bound(n,1)
            js = Bound(n,2)
            ks = Bound(n,3)
            ie = Bound(n,4)
            je = Bound(n,5)
            ke = Bound(n,6)
            face = 0
            if( (is == ie).and.(is ==  1) ) face=5
            if( (is == ie).and.(is == CI) ) face=3
            if( (js == je).and.(js ==  1) ) face=2
            if( (js == je).and.(js == CJ) ) face=4
            if( (ks == ke).and.(ks ==  1) ) face=1
            if( (ks == ke).and.(ks == CK) ) face=6
        end if

        do i=is,ie
            do j=js,je
                do k=ks,ke
                    c = BlkRes(b,6) + (k-1)*CI*CJ + (j-1)*CI + i
                    if(face /= 0) then
                        CellC(c,face)= -Bound(n,8) ! marker
                    else
                        material(c)= Bound(n,8)    ! material
                    end if
                end do
            end do
        end do

    end do  !  Nboun

END SUBROUTINE Calc1

max_stars_repo_name: palkinev/T-FlowS-old-compressible

max_issues_repo_name: palkinev/T-FlowS-old-compressible

max_forks_repo_name: palkinev/T-FlowS-old-compressible


text:
program scaleTest_coretran

!! Test program for the coretran library
  use iso_fortran_env
  use variableKind
  use m_array1D
  use m_errors
  use m_fileIO, only: openFile, closeFile
  use m_indexing
  use m_random
  use m_readLine
  use m_strings
  use m_writeLine
  use m_time
  use m_allocate
  use m_reallocate
  use m_copy
  ! Sorting routines
  use m_Sort
  !use m_PartialQuicksort
  use m_Select
  use m_maths

  use m_searching, only: binarySearch

  use Stopwatch_Class
  use ProgressBar_Class
  use m_KdTree

  implicit none

  character(len=100) :: fName
  character(len=128) :: sa
  character(len=:), allocatable :: cTest
  logical :: lTest
  integer(i32) :: iTest,istat,N, nIterations
  real(r32) :: ar, br, cr
  real(r32),allocatable :: ar1D(:),br1D(:),cr1D(:)
  real(r32),allocatable :: ar2D(:,:),br2D(:,:)
  real(r32),allocatable :: ar3D(:,:,:),br3D(:,:,:)
  real(r64) :: a,b,c
  real(r64),allocatable :: a1D(:),b1D(:),c1D(:)
  real(r64),allocatable :: a2D(:,:),b2D(:,:)
  real(r64),allocatable :: a3D(:,:,:),b3D(:,:,:)
  integer(i32) :: ia,ib,ic,id,ie
  integer(i32), allocatable :: ia1D(:),ib1D(:),ic1D(:)
  integer(i32), allocatable :: ia2D(:,:),ib2D(:,:)
  integer(i32), allocatable :: ia3D(:,:,:),ib3D(:,:,:)
  integer(i64) :: iad,ibd,icd
  integer(i64), allocatable :: iad1D(:),ibd1D(:),icd1D(:)
  integer(i64), allocatable :: iad2D(:,:),ibd2D(:,:)
  integer(i64), allocatable :: iad3D(:,:,:),ibd3D(:,:,:)
  complex(r32) :: x,y,z
  complex(r32), allocatable :: za1D(:),zb1D(:)
  complex(r32), allocatable :: za2D(:,:),zb2D(:,:)
  complex(r32), allocatable :: za3D(:,:,:),zb3D(:,:,:)
  complex(r64), allocatable :: zza1D(:),zzb1D(:)
  complex(r64), allocatable :: zza2D(:,:),zzb2D(:,:)
  complex(r64), allocatable :: zza3D(:,:,:),zzb3D(:,:,:)
  logical :: la,lb,lc
  logical, allocatable :: la1D(:),lb1D(:)
  logical, allocatable :: la2D(:,:),lb2D(:,:)
  logical, allocatable :: la3D(:,:,:),lb3D(:,:,:)
  integer(i32) :: passCount = 0
  integer(i32) :: testTotal = 0
  type(Stopwatch) :: clk
  type(ProgressBar) :: P
  type(KdTree) :: tree
  type(KdTreeSearch) :: search

  integer(i32) :: iunit
  integer(i32) :: maxSize = 25
  integer(i32), parameter :: nSizes = 6
  integer(i32) :: sizes(nSizes)
  real(r64) :: times(nSizes)

  ! Open a file for python plotting
  call openFile('scale_results.txt', iunit, 'unknown', istat)

  ! Set the max length of the tests
  N = 2**maxSize
  ! Set a repetition for very quick algorithms
  ic = 1e6

  ! Set the sizes
  do ia = 0, nSizes-1
    sizes(ia+1) = 2**(maxSize-ia)
  enddo

  write(iunit, '(a)') str(sizes)

  ! Set a fixed seed
  ! call random_seed(size = ia)
  ! call allocate(ia1D, ia)
  ! ia1D = 546420601
  call setPrng(big = .true., display=.true.)

    ! Initial setup for testing
  call allocate(ar1D, N)
  call allocate(br1D, N)
  call allocate(a1D, N)
  call allocate(b1D, N)
  call allocate(c1D, N)
  call allocate(ia1D,N)

  call rngNormal(a1D)
  ar1D = real(a1D)

  call Msg('==========================')
  call Msg('Testing : Sorting')
  call Msg('==========================')

  call Msg('Size  Time(s)')
  call Msg('Double precision introspection sort', iunit)

  b1D = a1D
  call clk%restart()
  call sort(b1D)
  call clk%stop()
  write(*,*) clk%elapsedInSeconds()

  do ia = 0, nSizes-1
    ib = sizes(ia+1)
    b1D(1:ib) = a1D(1:ib)

    call clk%restart()
    call sort(b1D(1:ib))
    call clk%stop()
    times(ia+1) = clk%elapsedInSeconds()
    call msg(str(ib)//' '//str(times(ia+1)))

    if (ia == 0) c1D=b1D
  enddo

  ! Write the times
  write(iunit, '(a)') str(times)

  call Msg('Already sorted double precision introspection sort', iunit)

  b1D = c1D

  do ia = 0, 5
    ib = 2**(24-ia)
    call clk%restart()
    call sort(b1D(1:ib))
    call clk%stop()
    times(ia+1) = clk%elapsedInSeconds()
    call msg(str(ib)//' '//str(times(ia+1)))
  enddo

  ! Write the times
  write(iunit, '(a)') str(times)

  call Msg('Double precision merge sort', iunit)
  do ia = 0, 5
    ib = 2**(24-ia)
    b1D(1:ib) = a1D(1:ib)

    call clk%restart()
    call sort(b1D(1:ib), stable=.true.)
    call clk%stop()
    times(ia+1) = clk%elapsedInSeconds()
    call msg(str(ib)//' '//str(times(ia+1)))
  enddo

  ! Write the times
  write(iunit, '(a)') str(times)

  call Msg('Already sorted double precision merge sort', iunit)
  b1D = c1D
  do ia = 0, 5
    ib = 2**(24-ia)

    call clk%restart()
    call sort(b1D(1:ib), stable=.true.)
    call clk%stop()
    times(ia+1) = clk%elapsedInSeconds()
    call msg(str(ib)//' '//str(times(ia+1)))
  enddo

  ! Write the times
  write(iunit, '(a)') str(times)

  call Msg('==========================')
  call Msg('Testing : Binary Search')
  call Msg('==========================')

  call allocate(ia1D, ic)
  b1D = c1D

  call Msg('Double precision binary search', iunit)
  do ia = 0, 5
    ib = sizes(ia+1)
    call rngInteger(ia1D, 1, ib)

    call clk%restart()
    do id = 1, ic
      a = b1D(ia1D(id))
      ie = binarySearch(b1D(1:ib), a, 1, ib)
    enddo
    call clk%stop()
    times(ia+1) = clk%elapsedInSeconds()/dble(ic)
    call msg(str(ib)//' '//str(times(ia+1)))

  enddo

  write(iunit,'(a)') str(times)

  call allocate(ia1D,N)

  call Msg('==========================')
  call Msg('Testing : Selection')
  call Msg('==========================')

  call Msg('Double precision quick select', iunit)

  ic = 1
  do ia = 0, 5
    ib = 2**(24-ia)

    ie = (ib+1)/2

    call clk%restart()
    do id = 1, ic
      b1D(1:ib) = a1D(1:ib)
      call select(b1D(1:ib), ie, a)
    enddo
    call clk%stop()
    times(ia+1) = clk%elapsedInSeconds()/dble(ic)
    call msg(str(ib)//' '//str(times(ia+1)))
  enddo
  write(iunit,'(a)') str(times)

  call Msg('Already sorted double precision quick select', iunit)
  b1D = c1D
  do ia = 0, 5
    ib = 2**(24-ia)

    ic = (ib+1)/2

    call clk%restart()
    call select(b1D(1:ib), ic, a)
    call clk%stop()
    times(ia+1) = clk%elapsedInSeconds()/dble(ic)
    call msg(str(ib)//' '//str(times(ia+1)))
  enddo
  write(iunit,'(a)') str(times)

  call Msg('==========================')
  call Msg('Testing : Spatial')
  call Msg('==========================')
  call Msg('---  Timing the 2D KdTree ---')
  call rngNormal(a1D)
  call rngNormal(b1D)
  do ia = 0, 5
    ib = 2**(24-ia)

    call clk%restart()
    tree = KdTree(a1D(1:ib), b1D(1:ib))
    call clk%stop()
    write(*,'(i10,1x,f0.3)') ib,clk%elapsedInSeconds()
  enddo
  call tree%deallocate()

  call Msg('---  Timing the 3D KdTree ---')
  call rngNormal(c1D)
  do ia = 0, 5
    ib = 2**(24-ia)

    call clk%restart()
    tree = KdTree(a1D(1:ib), b1D(1:ib), c1D(1:ib))
    call clk%stop()
    write(*,'(i10,1x,f0.3)') ib,clk%elapsedInSeconds()
  enddo
  call tree%deallocate()

  call allocate(a2D, [N, 5])
  call Msg('---  Timing the ND KdTree ---')
  call rngNormal(a2D)
  do ia = 0, 5
    ib = 2**(24-ia)

    call clk%restart()
    tree = KdTree(a2D(1:ib,:))
    call clk%stop()
    write(*,'(i10,1x,f0.3)') ib,clk%elapsedInSeconds()
  enddo
  call tree%deallocate()

  call closeFile('scale_results.txt',iunit,'',istat)

  stop
1 format(a)
contains

  subroutine test(l,msg)
  logical :: l
  character(len=*) :: msg
    testTotal = testTotal + 1
    if (l) then
      write(output_unit,1) msg//' Passed!'
      passCount = passCount + 1
    else
      write(output_unit,1) msg//' Failed!'
    endif
  1 format(a)
  end subroutine

end program

max_stars_repo_name: leonfoks/coretran

max_issues_repo_name: leonfoks/coretran

max_forks_repo_name: leonfoks/coretran


text:
c*********************************************************************c
c c
      subroutine dbzcalc(qvp,qra,qsn,qgr,tmk,prs,dbz,miy,mjx,mkzh,
     &                   in0r,in0s,in0g,iliqskin)
c
c This routine computes equivalent reflectivity factor (in dBZ) at
c each model grid point. In calculating Ze, the RIP algorithm makes
c assumptions consistent with those made in an early version
c (ca. 1996) of the bulk mixed-phase microphysical scheme in the MM5
c model (i.e., the scheme known as "Resiner-2"). For each species:
c
c 1. Particles are assumed to be spheres of constant density. The
c densities of rain drops, snow particles, and graupel particles are
c taken to be rho_r = rho_l = 1000 kg m^-3, rho_s = 100 kg m^-3, and
c rho_g = 400 kg m^-3, respectively. (l refers to the density of
c liquid water.)
c
c 2. The size distribution (in terms of the actual diameter of the
c particles, rather than the melted diameter or the equivalent solid
c ice sphere diameter) is assumed to follow an exponential
c distribution of the form N(D) = N_0 * exp( lambda*D ).
c
c 3. If in0X=0, the intercept parameter is assumed constant (as in
c early Reisner-2), with values of 8x10^6, 2x10^7, and 4x10^6 m^-4,
c for rain, snow, and graupel, respectively. Various choices of
c in0X are available (or can be added). Currently, in0X=1 gives the
c variable intercept for each species that is consistent with
c Thompson, Rasmussen, and Manning (2004, Monthly Weather Review,
c Vol. 132, No. 2, pp. 519-542.)
c
c 4. If iliqskin=1, frozen particles that are at a temperature above
c freezing are assumed to scatter as a liquid particle.
c
c More information on the derivation of simulated reflectivity in RIP
c can be found in Stoelinga (2005, unpublished write-up). Contact
c Mark Stoelinga (stoeling@atmos.washington.edu) for a copy.
c
      dimension qra(miy,mjx,mkzh),qvp(miy,mjx,mkzh),
     & qsn(miy,mjx,mkzh),qgr(miy,mjx,mkzh),tmk(miy,mjx,mkzh),
     & prs(miy,mjx,mkzh),dbz(miy,mjx,mkzh)
      character meth*1
c

!include 'comconst'

c
c Constant intercepts
c
      rn0_r = 8.e6 ! m^-4
      rn0_s = 2.e7 ! m^-4
      rn0_g = 4.e6 ! m^-4
c
c Constants used to calculate variable intercepts
c
      r1=1.e-15
      ron=8.e6
      ron2=1.e10
      son=2.e7
      gon=5.e7
      ron_min = 8.e6
      ron_qr0 = 0.00010
      ron_delqr0 = 0.25*ron_qr0
      ron_const1r = (ron2-ron_min)*0.5
      ron_const2r = (ron2+ron_min)*0.5
c
c Other constants:
c
      gamma_seven = 720.
      rho_r = rhowat ! 1000. kg m^-3
      rho_s = 100. ! kg m^-3
      rho_g = 400. ! kg m^-3
      alpha = 0.224
      factor_r = gamma_seven * 1.e18 * (1./(pi*rho_r))**1.75
      factor_s = gamma_seven * 1.e18 * (1./(pi*rho_s))**1.75
     & * (rho_s/rhowat)**2 * alpha
      factor_g = gamma_seven * 1.e18 * (1./(pi*rho_g))**1.75
     & * (rho_g/rhowat)**2 * alpha
c
      do k=1,mkzh
      do j=1,mjx-1
      do i=1,miy-1
c
      rhoair=prs(i,j,k)*100./
     & (rgas*virtual(tmk(i,j,k),qvp(i,j,k))) ! air density
c
c Adjust factor for brightband, where snow or graupel particle
c scatters like liquid water (alpha=1.0) because it is assumed to
c have a liquid skin.
c
      if (iliqskin.eq.1.and.tmk(i,j,k).gt.celkel) then
        factorb_s=factor_s/alpha
        factorb_g=factor_g/alpha
      else
        factorb_s=factor_s
        factorb_g=factor_g
      endif
c
c Calculate variable intercept parameters if wanted
c
      if (in0s.eq.1) then ! N_0s as in Thompson et al.
        temp_c = amin1(-0.001, tmk(i,j,k)-celkel)
        sonv = amin1(2.0e8, 2.0e6*exp(-0.12*temp_c))
      else
        sonv = rn0_s
      endif
c
      if (in0g.eq.1) then ! N_0g as in Thompson et al.
        gonv = gon
        if (qgr(i,j,k).gt.r1) then
          gonv = 2.38*(pi*rho_g/
     +          (rhoair*qgr(i,j,k)))**0.92
          gonv = max(1.e4, min(gonv,gon))
        endif
      else
        gonv = rn0_g
      endif
c
      if (in0r.eq.1) then ! N_0r as in Thompson et al.
        ronv = ron2
        if (qra(i,j,k).gt. r1) then
          ronv = ron_const1r*tanh((ron_qr0-qra(i,j,k))
     +          /ron_delqr0) + ron_const2r
        endif
c if (i.eq.97.and.j.eq.100) then
c print*,'k,ron2,qra(i,j,k),r1,ronv,rn0_r='
c print*,' ',k,ron2,qra(i,j,k),r1,ronv,rn0_r
c endif
      else
        ronv = rn0_r
      endif
c
c Total equivalent reflectivity factor (z_e, in mm^6 m^-3) is
c the sum of z_e for each hydrometeor species:
c
      z_e = factor_r * (rhoair*qra(i,j,k))**1.75 / ! rain
     & ronv**.75
     & + factorb_s * (rhoair*qsn(i,j,k))**1.75 / ! snow
     & sonv**.75
     & + factorb_g * (rhoair*qgr(i,j,k))**1.75 / ! graupel
     & gonv**.75
c
c Adjust small values of Z_e so that dBZ is no lower than -20
c
      z_e = max(z_e,.01)
c
c Convert to dBZ
c
      dbz(i,j,k) = 10. * log10(z_e)
c
      enddo
      enddo
      enddo
c
      return
      end

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
module semi_implicit_variables
   use real_kind
   implicit none

   ! Expansion coefficients for the radius in spherical harmonics, for rotation and tides
   real(double), allocatable, save :: a_rot2(:,:), a_tide(:, :, :)

   ! Ratio between the average radius r0 and the volume radius rv
   real(double), allocatable, save :: r0_over_rv(:, :)

   ! fp and ft correction terms for distorted stars
   real(double), allocatable, save :: fp_star(:, :)
   real(double), allocatable, save :: ft_star(:, :)

   ! Vertical component of the meridional circulation
   real(double), allocatable, save :: Vmc2_star(:, :)

   ! Difference in rotation rate, ~domega / dk
   real(double), allocatable, save :: diff_omega(:, :)

   ! Total mass loss rate due to RLOF
   real(double), save :: mdot_rlof(2)
   real(double), save :: mdot_rlof0(2)
end module semi_implicit_variables

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
MODULE MWwaterCoeff_FASTEM6

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module use
  USE Type_Kinds         , ONLY: fp
  USE Message_Handler    , ONLY: SUCCESS, Display_Message
  USE FitCoeff_Define    , ONLY: FitCoeff_SetValue
  USE MWwaterCoeff_Define, ONLY: MWwaterCoeff_type, &
                                 MWwaterCoeff_Create
  ! Disable implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  PUBLIC :: MWwaterCoeff_LoadCoeffs
  PUBLIC :: MWwaterCoeff_LoadVersion

  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &

  ! No. of Stokes components
  INTEGER , PARAMETER :: N_STOKES = 4
  INTEGER , PARAMETER :: N_2STOKES = 2  ! For data with just V and H polarisations

  ! Foam coverage coefficients from equation (5) of
  !
  !   Monahan, E.C., and O'Muircheartaigh, I.G., (1986)
  !     Whitecaps and the passive remote sensing of the ocean surface,
  !     International Journal of Remote Sensing, 7, pp627-642.
  !
  ! but assuming a neutral stability condition, i.e. dT=0, no
  ! difference between the skin and air temperature.
  INTEGER , PARAMETER :: N_FCCOEFFS = 2
  REAL(fp), PARAMETER :: FCCOEFF(N_FCCOEFFS) = (/1.95E-05_fp, 2.55_fp/)

  ! Foam reflectivity coefficients from section d of
  !
  !   Kazumori, M. et al. (2008) Impact Study of AMSR-E Radiances
  !     in the NCEP Global Data Assimilation System,
  !     Monthly Weather Review, 136, pp541-559
  !
  ! A constant value, 0.93, is used as the nadir emissivity, as
  ! well as for all incident angles for vertical polarisation.
  ! The relationship for horizontal polarisation is as listed
  ! in equation (15), same as the Stogryn(1972) model, but for
  ! the new nadir emissivity.
  !
  ! The current frequency dependent term is empirical, and undocumented.
  INTEGER , PARAMETER :: N_FRCOEFFS = 6
  REAL(fp), PARAMETER :: FRCOEFF(N_FRCOEFFS) = &
  (/ 0.93_fp, &                                 ! Nadir emissivity
    -1.748e-3_fp, -7.336e-5_fp, 1.044e-7_fp, &  ! Horizontal polarisation coefficients
     0.4_fp, -0.05_fp /)                        ! Frequency correction coefficients

  ! Small-scale correction coefficients for equation (A4) of
  !
  !   Liu, Q. et al. (2011) An Improved Fast Microwave Water
  !     Emissivity Model, TGRSS, 49, pp1238-1250
  !
  ! where the small-scale correction formulation is given in
  ! equation (17a,b) of
  !
  !   Liu, Q. et al. (1998) Monte Carlo simulations of the microwave
  !     emissivity of the sea surface, JGR, 103, pp24983-24989
  !
  ! and originally in equation (30) of
  !
  !   Guissard,A. and P.Sobieski (1987) An approximate model
  !     for the microwave brightness temperature of the sea,
  !     Int.J.Rem.Sens., 8, pp1607-1627.
  INTEGER , PARAMETER :: N_SSCCOEFFS = 8
  REAL(fp), PARAMETER :: SSCCOEFF(N_SSCCOEFFS) = &
  (/ -5.0208480E-06_fp,   2.3297951E-08_fp,   4.6625726E-08_fp,  -1.9765665E-09_fp, &
     -7.0469823E-04_fp,   7.5061193E-04_fp,   9.8103876E-04_fp,   1.5489504E-04_fp /)

  ! Large-scale correction coefficients for equations (A5a,b) of
  !
  !   Liu, Q. et al. (2011) An Improved Fast Microwave Water
  !     Emissivity Model, TGRSS, 49, pp1238-1250
  INTEGER , PARAMETER :: N_FCOEFFS   = 3 ! No. of frequency polynomial coefficients
  INTEGER , PARAMETER :: N_LSCCOEFFS = 6 ! No. of regression coefficients
  REAL(fp), PARAMETER :: LSCCOEFF(N_FCOEFFS, N_LSCCOEFFS, N_2STOKES) = RESHAPE( &
  (/ -5.994667E-02_fp, 9.341346E-04_fp,-9.566110E-07_fp, 8.360313E-02_fp,-1.085991E-03_fp, 6.735338E-07_fp, &
     -2.617296E-02_fp, 2.864495E-04_fp,-1.429979E-07_fp,-5.265879E-04_fp, 6.880275E-05_fp,-2.916657E-07_fp, &
     -1.671574E-05_fp, 1.086405E-06_fp,-3.632227E-09_fp, 1.161940E-04_fp,-6.349418E-05_fp, 2.466556E-07_fp, &

     -2.431811E-02_fp,-1.031810E-03_fp, 4.519513E-06_fp, 2.868236E-02_fp, 1.186478E-03_fp,-5.257096E-06_fp, &
     -7.933390E-03_fp,-2.422303E-04_fp, 1.089605E-06_fp,-1.083452E-03_fp,-1.788509E-05_fp, 5.464239E-09_fp, &
     -3.855673E-05_fp, 9.360072E-07_fp,-2.639362E-09_fp, 1.101309E-03_fp, 3.599147E-05_fp,-1.043146E-07_fp /), &
  (/N_FCOEFFS, N_LSCCOEFFS, N_2STOKES/) )

  ! Reflection correction coefficients
  INTEGER , PARAMETER :: N_ODCOEFFS = 3 ! No. of optical depth coefficients
  INTEGER , PARAMETER :: N_RCCOEFFS = 7 ! No. of regression coefficients
  REAL(fp), PARAMETER :: RCCOEFF(N_ODCOEFFS, N_RCCOEFFS, N_2STOKES) = RESHAPE( &
  (/  0.388242E-01_fp, 0.194901E+00_fp,-0.425093E-01_fp, &
      0.607698E+01_fp,-0.313861E+01_fp,-0.103383E+01_fp, &
     -0.377867E+01_fp, 0.180284E+01_fp, 0.699556E+00_fp, &
     -0.506455E-01_fp,-0.262822E+00_fp, 0.703056E-01_fp, &
      0.362055E+01_fp,-0.120318E+01_fp,-0.124971E+01_fp, &
      0.154014E-01_fp, 0.759848E-01_fp,-0.268604E-01_fp, &
     -0.802073E+01_fp, 0.324658E+01_fp, 0.304165E+01_fp, &

      0.199277E+00_fp, 0.166155E+00_fp, 0.153272E-01_fp, &
      0.399234E+01_fp,-0.130968E+01_fp,-0.874716E+00_fp, &
     -0.169403E+01_fp,-0.260998E-01_fp, 0.540443E+00_fp, &
     -0.282483E+00_fp,-0.219994E+00_fp,-0.203438E-01_fp, &
      0.351731E+00_fp, 0.208641E+01_fp,-0.693299E+00_fp, &
      0.867861E-01_fp, 0.619020E-01_fp, 0.595251E-02_fp, &
     -0.475191E+01_fp,-0.430134E-01_fp, 0.248524E+01_fp /), &
  (/N_ODCOEFFS, N_RCCOEFFS, N_2STOKES/) )

  ! Coefficients for M.Kazumori azimuth model function
  INTEGER , PARAMETER :: N_AZCOEFFS      = 6 ! No. of azimuth emissivity coefficients
  INTEGER , PARAMETER :: N_AZFREQUENCIES = 6 ! No. of fitting frequencies
  REAL(fp), PARAMETER :: AZCOEFF(N_AZCOEFFS, N_AZFREQUENCIES, N_2STOKES) = RESHAPE( &
  [ 4.401E-02_fp, -1.636E+01_fp,  1.478E+00_fp, -4.800E-02_fp,  3.202E-06_fp, -6.002E-05_fp, &     ! 06V OK
    4.379E-02_fp, -1.633E+01_fp,  1.453E+00_fp, -4.176E-02_fp,  5.561E-06_fp, -4.644E-05_fp, &     ! 10V OK
    5.009E-02_fp, -1.638E+01_fp,  1.520E+00_fp, -3.994E-02_fp,  1.330E-05_fp,  1.113E-05_fp, &     ! 19V OK
    5.165E-02_fp, -1.638E+01_fp,  1.543E+00_fp, -4.066E-02_fp,  1.494E-05_fp,  1.010E-05_fp, &     ! 23V interpolated
    5.553E-02_fp, -1.638E+01_fp,  1.602E+00_fp, -4.246E-02_fp,  1.903E-05_fp,  7.524E-06_fp, &     ! 37V OK
   -9.131E-05_fp,  1.251E+00_fp,  6.769E-01_fp, -2.913E-02_fp,  1.092E+00_fp, -1.806E-04_fp, &     ! 89V OK revised
   -1.234E-07_fp, -8.179E-03_fp, -1.040E+01_fp,  4.477E-01_fp,  0.000E+00_fp,  3.390E-05_fp, &     ! 06H OK
   -1.938E-05_fp, -8.007E-03_fp, -1.039E+01_fp,  4.610E-01_fp,  0.000E+00_fp,  4.419E-05_fp, &     ! 10H OK
    1.362E-04_fp, -1.013E-03_fp, -9.235E+00_fp,  3.844E-01_fp,  0.000E+00_fp,  2.891E-04_fp, &     ! 19H OK
    1.519E-04_fp, -7.865E-04_fp, -9.234E+00_fp,  3.884E-01_fp,  0.000E+00_fp,  6.856E-04_fp, &     ! 23H Interpolated
    1.910E-04_fp, -2.224E-04_fp, -9.232E+00_fp,  3.982E-01_fp,  0.000E+00_fp,  1.673E-03_fp, &     ! 37H OK
    3.554E-04_fp,  5.226E-04_fp,  9.816E-01_fp, -7.783E-03_fp,  0.000E+00_fp,  2.437E+01_fp  ], &  ! 89H OK revised
    [N_AZCOEFFS, N_AZFREQUENCIES, N_2STOKES] )

CONTAINS

  SUBROUTINE MWwaterCoeff_LoadCoeffs( MWwaterCoeff )
    TYPE(MWwaterCoeff_type), INTENT(IN OUT) :: MWwaterCoeff
    INTEGER :: version

    version = 6

    ! Create an instance of the structure
    CALL MWwaterCoeff_Create(MWwaterCoeff)
    MWwaterCoeff%Version = version

    ! Load the foam coverage coefficients
    CALL FitCoeff_SetValue(MWwaterCoeff%FCCoeff, FCCOEFF, Version=version)

    ! Load the foam reflectivity coefficients
    CALL FitCoeff_SetValue(MWwaterCoeff%FRCoeff, FRCOEFF, Version=version)

    ! Load the reflection correction coefficients
    CALL FitCoeff_SetValue(MWwaterCoeff%RCCoeff, RCCOEFF, Version=version)

    ! Load the azimuth emissivity coefficients
    CALL FitCoeff_SetValue(MWwaterCoeff%AZCoeff, AZCOEFF, Version=version)

    ! Load the small-scale correction coefficients
    CALL FitCoeff_SetValue(MWwaterCoeff%SSCCoeff, SSCCOEFF, Version=version)

    ! Load the large-scale correction coefficients
    CALL FitCoeff_SetValue(MWwaterCoeff%LSCCoeff, LSCCOEFF, Version=version)

  END SUBROUTINE MWwaterCoeff_LoadCoeffs

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

END MODULE MWwaterCoeff_FASTEM6

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm