text:
      SUBROUTINE HRT1 ( RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,TBOT,
     +                 ZBOT, PSISAT, SH2O, DT, B, SH, SR, SF,
C
C12/02  V.K.  ADDED NEW FORMAL PARAMETER - IVERS
CC     +                 STYPE, QUARTZ, SMCMR )
     +                 STYPE, QUARTZ, SMCMR, IVERS, CKSOIL)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY
CC    =======   TERM OF THE SOIL THERMAL DIFFUSION EQUATION.  ALSO TO
CC              COMPUTE ( PREPARE ) THE MATRIX COEFFICIENTS FOR THE
CC              TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
CC     TSS - SNOW-BARE SOIL WEGHTED TEMPERATURE DURING SNOWMELT
CC     SOILT - SOIL SURFACE TEMPERATURE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER NSOIL, K
      REAL CAIR, CH2O, CSOIL, DDZ, DDZ2, DENOM, DFMAX1, DFMAXN,
     +     DTSDZ, DTSDZ2, HCPCT, MXICE, S, SMCMAX, STYPE,
     +     TBOT, ZBOT, YY, SH, SR, SF, SOILT, TSS
      REAL SMC   ( NSOIL )
      REAL SH2O  ( NSOIL )
      REAL STC   ( NSOIL )
      REAL ZSOIL ( NSOIL )
      REAL RHSTS ( NSOIL )

      PARAMETER ( NSZ_ABC = 10 )
      REAL AI    ( NSZ_ABC )
      REAL BI    ( NSZ_ABC )
      REAL CI    ( NSZ_ABC )
      COMMON /ABC/ AI, BI, CI

c  test only
      common /test1/ iday

      PARAMETER( CAIR=1004.0,CH2O=4.2 E6,CSOIL=1.26E6 )
      PARAMETER ( CICE = 2.106 E6, T0 = 273.16 )

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     RETRIEVE THE THERMAL DIFFUSIVITY
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C  ACCOUNT FOR RESIDUE EFFECT ON CONDUCTIVITY
      SICE = 0.
      DO I = 1, NSOIL
        SICE = SICE + SMC(I) - SH2O(I)
      ENDDO

c3/10 surfwater
      smc0=smcmr*0.15
      sh2o0=smc0
      MXICE = SMC0 - SH2O0
      dfmax1 = cnd_johns(stype,sh2o0,mxice,smcmr,0.)
c3/10 surfwater end

C  CALCULATE SNOW/SOIL FLUX AND SOIL SURFACE TEMPERATURE IF SNOW
      SOILT=YY
      TSS=-99.
      IF(SH .NE. 0.) THEN
       TSS=YY
       IF(YY .GT. T0) TSS=SF*T0+(1-SF)*YY
       DFSNOW=CSNOW(SR)

C  CALCULATION OF SOIL TEMPERATURE IS BASED ON Qsnow=Qsoil
       ALP=1./(1. + DFMAX1*SH/(-DFSNOW*0.5*ZSOIL(1)))
       SOILT=ALP*TSS+(1-ALP)*STC(1)
      ENDIF
      S = DFMAX1 * ( STC(1) - SOILT ) / ( 0.5 * ZSOIL(1) )
      S = MIN ( MAX ( S,-100.0 ), 100.0 )

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC THE HEAT CAPACITY OF THE TOP SOIL LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      HCPCT = SH2O0*CH2O + (1.0-SMCMR)*CSOIL + (SMCMR-SMC0)*CAIR
     +        + ( SMC0 - SH2O0 ) * CICE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DDZ = 1.0 / ( -0.5 * ZSOIL(2) )
      AI(1) = 0.0
      CI(1) =  ( DFMAX1 * DDZ ) / ( ZSOIL(1) * HCPCT )
      BI(1) = -CI(1)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC THE VERTICAL SOIL TEMP GRADIENT BTWN THE TOP AND 2ND SOIL
C     LAYERS.  RECALC/ADJUST THE SOIL HEAT FLUX.  USE THE GRADIENT
C     AND FLUX TO CALC RHSTS FOR THE TOP SOIL LAYER.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DTSDZ = ( STC(1) - STC(2) ) / ( -0.5 * ZSOIL(2) )
      RHSTS(1) = ( DFMAX1 * DTSDZ - S ) / ( ZSOIL(1) * HCPCT )

C   CALCULATING OF LAYER BOUNDARY TEMPERATURE USING TBND SUBR.
C   AND SINK/SOURCE TERM OF DIFFUSION EQUATION BY SNKSRC SUBR.

      TBK = TBND ( STC(1), STC(2), ZSOIL, ZBOT, 1, NSOIL )
      IF ( SICE .NE. 0. .OR. SOILT .LT. T0 ) THEN
       TSNSR = SNKSRC ( SOILT, STC(1),TBK, SMC(1), SH2O(1), HCPCT,
C
C12/02  V.K.  ADDED NEW FORMAL PARAMETER - IVERS
CC     +             ZSOIL, SMCMR, PSISAT, B, DT, 1, STYPE, 0.)
c 3/10 surface water
     +     ZSOIL, SMCMR, PSISAT, B, DT, 1, STYPE, 0.,IVERS,CKSOIL,smc0)

       RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
      ENDIF
      DDZ2 = 0.0

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DO 40 K = 2, NSOIL

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CALC THIS SOIL LAYER'S HEAT CAPACITY
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        HCPCT = SH2O(K)*CH2O +(1.0-SMCMAX)*CSOIL +(SMCMAX-SMC(K))*CAIR
     +        + ( SMC(K) - SH2O(K) ) * CICE

        IF ( K .NE. NSOIL ) THEN

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         RETRIEVE THE THERMAL DIFFUSIVITY
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        MXICE = SMC(K) - SH2O(K)
        dfmaxn=cnd_johns(stype,sh2o(k),mxice,smcmax,quartz)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

          DENOM = 0.5 * ( ZSOIL(K-1) - ZSOIL(K+1) )
          DTSDZ2 = ( STC(K) - STC(K+1) ) / DENOM

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

          DDZ2 = 2. / (ZSOIL(K-1) - ZSOIL(K+1))
          CI(K) = -DFMAXN * DDZ2 / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)

C   CALCULATING OF LAYER BOUNDARY TEMPERATURE USING TBND SUBR.
          TBK1 = TBND ( STC(K), STC(K+1), ZSOIL, ZBOT, K, NSOIL )

        ELSE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         RETRIEVE THE THERMAL DIFFUSIVITY FOR THE LOWEST SOIL LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C   SOIL THERMAL DIFFUSIVITY ADJUSTMENT
C   AND  CALCULATING OF LAYER BOUNDARY TEMPERATURE USING TBND SUBR.
C
          MXICE = SMC(K) - SH2O(K)
          dfmaxn=cnd_johns(stype,sh2o(k),mxice,smcmax,quartz)
          TBK1 = TBND ( STC(K), TBOT, ZSOIL, ZBOT, K, NSOIL )

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         CALC THE VERTICAL SOIL TEMP GRADIENT THRU THE LOWEST LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

          DTSDZ2 = (STC(K)-TBOT) / (.5 * (ZSOIL(K-1) + ZSOIL(K))-ZBOT)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         SET MATRIX COEF, CI TO ZERO
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

          CI(K) = 0.
        END IF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CALC RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        DENOM = ( ZSOIL(K) - ZSOIL(K-1) ) * HCPCT
        RHSTS(K) = ( DFMAXN * DTSDZ2 - DFMAX1 * DTSDZ ) / DENOM

C   SINK/SOURCE TERM OF DIFFUSION EQUATION BY SNKSRC SUBR.
      IF ( SICE .NE. 0. .OR. TBK .LT. T0 ) THEN
       TSNSR = SNKSRC ( TBK, STC(K),TBK1, SMC(K), SH2O(K), HCPCT,
C
C12/02  V.K.  ADDED NEW FORMAL PARAMETER - IVERS
CC     +             ZSOIL,SMCMAX, PSISAT, B, DT, K,STYPE,QUARTZ)
c 3/10 surface water
     +   ZSOIL,SMCMAX, PSISAT, B, DT, K,STYPE,QUARTZ,IVERS,CKSOIL,smc0)

       RHSTS(K) = RHSTS(K) - TSNSR / DENOM
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        AI(K) = - DFMAX1 * DDZ / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
        BI(K) = -(AI(K) + CI(K))

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       RESET VALUES OF DFMAX1, DTSDZ, AND DDZ FOR LOOP TO NEXT SOIL LYR
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        TBK = TBK1
        DFMAX1 = DFMAXN
        DTSDZ = DTSDZ2
        DDZ = DDZ2
   40   CONTINUE

      RETURN
      END

      SUBROUTINE HSTEP ( STCOUT, STCIN, RHSTS, DT, NSOIL )

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER NSOIL, K
      REAL DT
      REAL RHSTS ( NSOIL )
      REAL STCOUT   ( NSOIL )
      REAL STCIN  ( NSOIL )

      PARAMETER ( NSZ_ABC = 10 )
      REAL AI    ( NSZ_ABC )
      REAL BI    ( NSZ_ABC )
      REAL CI    ( NSZ_ABC )
      COMMON /ABC/ AI, BI, CI

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CREATE FINITE DIFFERENCE VALUES FOR USE IN ROSR12 ROUTINE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DO 10 K = 1 , NSOIL
        RHSTS(K) = RHSTS(K) * DT
        AI(K) = AI(K) * DT
        BI(K) = 1. + BI(K) * DT
        CI(K) = CI(K) * DT
 10   CONTINUE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     SOLVE THE TRI-DIAGONAL MATRIX EQUATION
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CALL ROSR12 ( CI,AI,BI,CI,RHSTS,RHSTS,NSOIL )

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     CALC/UPDATE THE SOIL TEMPS USING MATRIX SOLUTION
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DO 20 K = 1 , NSOIL
        STCOUT(K) = STCIN(K) + CI(K)
  20  CONTINUE
      RETURN
      END

      SUBROUTINE ROSR12 ( P, A, B, C, D, DELTA, NSOIL )

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC    PURPOSE:  TO INVERT (SOLVE) THE TRI-DIAGONAL MATRIX PROBLEM SHOWN
CC    =======   BELOW:
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER K, KK, NSOIL
      REAL P(NSOIL)
      REAL A(NSOIL)
      REAL B(NSOIL)
      REAL C(NSOIL)
      REAL D(NSOIL)
      REAL DELTA(NSOIL)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     INITIALIZE EQN COEF C FOR THE LOWEST SOIL LAYER.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      C(NSOIL) = 0.0

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     SOLVE THE COEFS FOR THE 1ST SOIL LAYER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      P(1) = -C(1) / B(1)
      DELTA(1) = D(1) / B(1)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     SOLVE THE COEFS FOR SOIL LAYERS 2 THRU NSOIL
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DO 10 K = 2 , NSOIL
        P(K) = -C(K) * ( 1.0 / (B(K) + A (K) * P(K-1)) )
        DELTA(K) = (D(K)-A(K)*DELTA(K-1))*(1.0/(B(K)+A(K)*P(K-1)))
 10   CONTINUE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     SET P TO DELTA FOR LOWEST SOIL LAYER.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      P(NSOIL) = DELTA(NSOIL)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     ADJUST P FOR SOIL LAYERS 2 THRU NSOIL
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DO 20 K = 2 , NSOIL
        KK = NSOIL - K + 1
        P(KK) = P(KK) * P(KK+1) + DELTA(KK)
 20   CONTINUE

      RETURN
      END

max_stars_repo_name: KathyNie/LISF

max_issues_repo_name: dmocko/LISF

max_forks_repo_name: dmocko/LISF


text:
*----------------------------------------------------------------------*
      subroutine add_unity_drv(label,fac,ssign,init,minblk_in,maxblk_in,
     &     op_info,orb_info,str_info)
*----------------------------------------------------------------------*
      implicit none

      include 'opdim.h'
      include 'stdunit.h'
      include 'mdef_operator_info.h'
      include 'def_graph.h'
      include 'def_strinf.h'
      include 'def_orbinf.h'

      character(len=*), intent(in) ::
     &     label
      real(8), intent(in) ::
     &     fac
      logical, intent(in) ::
     &     init
      integer, intent(in) ::
     &     minblk_in, maxblk_in, ssign
      type(operator_info) ::
     &     op_info
      type(strinf) ::
     &     str_info
      type(orbinf) ::
     &     orb_info

      integer ::
     &     idx, jdx, minblk, maxblk

      integer, external ::
     &     idx_mel_list

      type(me_list), pointer ::
     &     mel_pnt

      idx = idx_mel_list(label,op_info)
      if (idx.lt.0)
     &       call quit(1,'add_unity_drv','Label not on list: "'//
     &       trim(label)//'"')

      mel_pnt => op_info%mel_arr(idx)%mel

      minblk = 1
      if (minblk_in.gt.0) minblk = minblk_in
      maxblk = mel_pnt%op%n_occ_cls
      if (maxblk_in.gt.0) maxblk = maxblk_in

      if (init) call zeroop(mel_pnt)

      do idx = minblk, maxblk
        call add_unity(fac,ssign,mel_pnt,idx,orb_info,str_info)
      end do

      end

ext: f

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


module vec2los

character(len=512) :: input_file
character(len=512) :: output_file
double precision :: azimuth
double precision :: inclination
character(len=512) :: look_geometry_file

end module

!==================================================================================================!

program main
!----
! Read three-dimensional E-N-Z displacements and compute line-of-sight displacements, given the
! azimuth (clockwise from north) and inclination (upward from horizontal) from the viewer to the
! location of interest.
!----

use io, only: stderr, stdin, stdout, fileExists, line_count
use trig, only: d2r

use vec2los, only: input_file, &
                   output_file, &
                   azimuth, &
                   inclination, &
                   look_geometry_file

implicit none

! Local variables
double precision :: az, inc, cosaz, cosinc, sinaz, sininc, ue, un, uz, ulos
double precision, allocatable :: vector(:,:), look(:,:)
integer :: i, luin, luout, ios, nvec, nlook
integer, parameter :: NVECMAX = 10000
character(len=256) :: input_line

! Parse command line
call gcmdln()

! Check that look geometry is defined
if ((azimuth.lt.-1.0d50 .or. inclination.lt.-1.0d50) .and. look_geometry_file.eq.'') then
    call usage('vec2los: no look geometry defined')
endif

! Check that input file exists, if defined
if (input_file.ne.'') then
    if (.not.fileExists(input_file)) then
        call usage('vec2los: no input vector displacement file found named '//trim(input_file))
    endif
endif

! Open input and output streams
if (input_file.eq.'') then
    luin = stdin
    allocate(vector(NVECMAX,6))
else
    nvec = line_count(input_file)
    luin = 11
    open(unit=luin,file=input_file,status='old')
    allocate(vector(nvec,6))
endif

if (output_file.eq.'') then
    luout = stdout
else
    luout = 12
    open(unit=luout,file=output_file,status='unknown')
endif

! Read input displacements
nvec = 0
do
    read(luin,'(A)',iostat=ios) input_line
    if (ios.ne.0) then
        exit
    else
        nvec = nvec + 1
        read(input_line,*) (vector(nvec,i),i=1,6)
    endif

    ! Maximum number of inputs is limited when reading from stdin
    if (input_file.eq.''.and.nvec.gt.NVECMAX) then
        write(stderr,*) 'vec2los: number of inputs from stdin exceeds NVECMAX'
        call usage('read input vectors from file for large numbers of observations')
    endif
enddo

if (look_geometry_file.ne.'') then

    ! If there is a look file, read it

    if (.not.fileExists(look_geometry_file)) then
        call usage('vec2los: no look geometry file found named '//trim(look_geometry_file))
    endif

    ! Number of look values must equal number of inputs
    nlook = line_count(look_geometry_file)
    if (nlook.ne.nvec) then
        call usage('vec2los: number of look geometries is not equal to number of vector inputs')
    endif

    ! Allocate memory for reading inputs
    allocate(look(nlook,2))

    ! Read look geometry from file
    open(unit=21,file=look_geometry_file,status='old')
    do i = 1,nlook
        read(21,*) look(i,1),look(i,2)
    enddo
    close(21)

else

    ! Otherwise set the look geometry to be the same for all displacements

    allocate(look(nvec,2))
    look(:,1) = azimuth
    look(:,2) = inclination

endif

! Compute the line-of-sight displacement, from the viewer to the location
! Positive implies displacement away from viewer
do i = 1,nvec
    az = look(i,1)
    inc = look(i,2)
    cosaz = dcos(look(i,1)*d2r)
    sinaz = dsin(look(i,1)*d2r)
    cosinc = dcos(look(i,2)*d2r)
    sininc = dsin(look(i,2)*d2r)

    ue = vector(i,4)
    un = vector(i,5)
    uz = vector(i,6)
    ulos = un*cosaz*cosinc + ue*sinaz*cosinc - uz*sininc
    write(luout,*) vector(i,1),vector(i,2),vector(i,3),uLOS
enddo

! FINISH HIM...IT...THE PROGRAM
if (input_file.ne.'') then
    close(luin)
endif
if (output_file.ne.'') then
    close(luout)
endif

end

!--------------------------------------------------------------------------------------------------!

subroutine gcmdln()

use vec2los, only: input_file, &
                   output_file, &
                   azimuth, &
                   inclination, &
                   look_geometry_file

implicit none

! Local variables
integer :: i, narg
character(len=512) :: tag

! Initialize control parameters
input_file = ''
output_file = ''
azimuth = -1d100
inclination = -1d100
look_geometry_file = ''

! Number of arguments
narg = command_argument_count()
if (narg.eq.0) then
  call usage('')
endif

! Read arguments
i = 1
do while (i.le.narg)

    call get_command_argument(i,tag)

    if (trim(tag).eq.'-a') then
        i = i + 1
        call get_command_argument(i,tag)
        read (tag,*) azimuth
    elseif (trim(tag).eq.'-i') then
        i = i + 1
        call get_command_argument(i,tag)
        read (tag,*) inclination
    elseif (trim(tag).eq.'-look') then
        i = i + 1
        call get_command_argument(i,look_geometry_file)
    elseif (trim(tag).eq.'-f') then
        i = i + 1
        call get_command_argument(i,input_file)
    elseif (trim(tag).eq.'-o') then
        i = i + 1
        call get_command_argument(i,output_file)
    endif

    i = i + 1
enddo

return
end subroutine

!--------------------------------------------------------------------------------------------------!

subroutine usage(str)

use io, only: stderr
implicit none

character(len=*) :: str

if (str.ne.'') then
    write(stderr,*) trim(str)
    write(stderr,*)
endif

write(stderr,*) 'Usage: vec2los -a AZ -i INC | -look FILE [-f IFILE] [-o OFILE]'
write(stderr,*)
write(stderr,*) '-a AZ         Azimuth from viewer to ground (CW from N)'
write(stderr,*) '-i INC        Inclination from viewer to ground (horizontal=0, vertical=90)'
write(stderr,*) '-look FILE    Read azimuth and inclination from FILE'
write(stderr,*) '-f FILE       Input vector displacements (default: stdin) (lo la dp E N Z)'
write(stderr,*) '-o FILE       Output LOS displacements (default: stdout) (lo la dp LOS)'
write(stderr,*)

call error_exit(1)
end subroutine


RECURSIVE SUBROUTINE PDEFusedSrcNCP(Src_BgradQ,Q,gradQin)
    !USE MainVariables, ONLY : nVar, nParam, d, EQN, nDim
    IMPLICIT NONE
#if defined(Dim3)
    INTEGER, PARAMETER             :: nDim = 3                   ! The number of space dimensions
#elif defined(Dim2)
    INTEGER, PARAMETER             :: nDim = 2                   ! The number of space dimensions
#endif
    INTEGER, PARAMETER             :: nParam=0
    INTEGER, PARAMETER             :: d=3
    INTEGER, PARAMETER             :: nVar = 59                           ! The number of variables of the PDE system
  TYPE, bind(C) :: tEquations
      REAL(8)    :: gamma, Pi, c0, g = 9.81, friction = 1.0
      REAL(8)    :: CCZ4k1, CCZ4k2, CCZ4k3, CCZ4eta, CCZ4itau, CCZ4f, CCZ4g, CCZ4xi, CCZ4e, CCZ4c, CCZ4mu, CCZ4ds, CCZ4sk, CCZ4bs
      REAL(8)    :: CCZ4GLMc0 = 0.5, CCZ4GLMc = 0.75, CCZ4GLMd = 0.75, CCZ4GLMepsD = 1e-2, CCZ4GLMepsA = 1e-2, CCZ4GLMepsP = 1e-2, cs, alpha, beta, lambda, cv, rho0, p0, tau1, tau2, mu, kappa ,tau
      INTEGER :: CCZ4LapseType, EinsteinAutoAux = 0, ReferenceDepth = 1.0
      REAL(8)    :: DivCleaning_a = 1.0
  END TYPE tEquations
  TYPE(tEquations) :: EQN
    ! Argument list
    REAL, INTENT(IN)  :: Q(nVar), gradQin(nVar,d)
    REAL, INTENT(OUT) :: Src_BgradQ(nVar)
    ! Local variables
    REAL :: gradQ(nVar,d)
    REAL :: par(nParam), time=0.0
    INTEGER :: i,j,k,l,m,n,ip,iq,ii,jj,kk,ll,mm,iErr,count
    REAL :: p, irho, lam, mu
    REAL :: Qx(nVar), Qy(nVar), Qz(nVar), BgradQ(nVar), src(nVar), V(nVar)
    REAL :: k1,k2,k3,fff,ggg,e,c,ds,xi,sk,sknl,bs,g_cov(3,3),g_contr(3,3),dgup(3,3,3)
    REAL :: det, alpha, fa, k0, dk0(3), beta0(3), b0(3), u(3), ialpha
    REAL :: s1, s2, s3, s4, s5, s6, s7, s8, s9, s10
    REAL :: AQx(nVar), BQy(nVar), CQz(nVar), Smom(3), iphi2, ddm(3,3,3)
    REAL :: lapse, shift(3), gammaij(6), delta(3,3), bv(3), vxb(3), vxb_contr(3), psi, qv_contr(3), qb_contr(3), bv_contr(3)
    REAL :: v2,vf(3),uem,b2,e2,gp,gm,vc(nVar),lf,w,ww,gamma1,rho,vf_cov(3), s_contr(3), w_ij, wim
    REAL :: v_cov(3), v_contr(3), SijTF(3,3), phi2, EE

#ifdef VECTOR
#ifdef AVX512
  INTEGER, PARAMETER :: nVarGRMHD = 24                           ! The number of variables of the PDE system
#else
  INTEGER, PARAMETER :: nVarGRMHD = 20                           ! The number of variables of the PDE system
#endif
#else
  INTEGER, PARAMETER :: nVarGRMHD = 19                           ! The number of variables of the PDE system
#endif
    REAL :: VGRMHD(nVarGRMHD), QGRMHD(nVarGRMHD), gradQGRMHD(nVarGRMHD,d), BgradQGRMHD(nVarGRMHD), eta, itau, Ham, Mom(3), dKex(3,3,3), ov(3)
#ifdef VECTOR
      INTEGER, PARAMETER :: nVarGRGPR = 32                           ! The number of variables of the PDE system
#else
      INTEGER, PARAMETER :: nVarGRGPR = 30                           ! The number of variables of the PDE system
#endif
    REAL :: SGRGPR(nVarGRGPR),QGRGPR(nVarGRGPR), VGRGPR(nVarGRGPR), gradQGRGPR(nVarGRGPR,d), BgradQGRGPR(nVarGRGPR)
    !
    REAL :: Christoffel(3,3,3), RiemannNCP(3,3,3,3), RiemannSrc(3,3,3,3), dChristoffelNCP(3,3,3,3), dChristoffelSrc(3,3,3,3), DD(3,3,3), dDD(3,3,3,3), dChristoffelOrd(3,3,3,3)
    REAL :: dChristoffel_tildeNCP(3,3,3,3), dChristoffel_tildeSrc(3,3,3,3)
    REAL :: R, AA(3), dAA(3,3), BB(3,3), dBB(3,3,3), beta(3), Kex(3,3), Kmix(3,3), Kup(3,3), Z(3), dZ(3,3), nablaZNCP(3,3), nablaZSrc(3,3), RplusNablaZNCP, RplusNablaZSrc
    REAL :: Theta, dTheta(3), nablaijalphaNCP(3,3), nablaijalphaSrc(3,3), Ricci(3,3), RicciNCP(3,3), RicciSrc(3,3), dtraceK(3), dtraceKNCP(3), dKtempNCP(3), dZNCP(3,3), dZSrc(3,3)
    REAL :: dtgamma(3,3), dtK(3,3), dK(3,3,3), dtTheta, dtZ(3), dtalpha, dtGhat(3), dtbeta(3), dtbb(3), dtA(3), dtB(3,3), dtD(3,3,3)
    REAL :: Aupdown, Aex(3,3), dAex(3,3,3), Amix(3,3), Aup(3,3), Ghat(3), Gtilde(3), dGhat(3,3), traceK, Kupdown, phi, PP(3), dPP(3,3), Pup(3), DDcontr(3)
    REAL :: dGtildeSrc(3,3), dGtildeNCP(3,3), RiccitildeNCP(3,3), RicciphiNCP(3,3), RiccitildeSrc(3,3), RicciphiSrc(3,3)
    REAL :: Christoffel_tilde(3,3,3), Christoffel_kind1(3,3,3), Zup(3), RicciPlusNablaZNCP(3,3), RicciPlusNablaZSrc(3,3), traceA, traceB, QG(3), b(3), faa, temp
    REAL :: SecondOrderTermsNCP(3,3), SecondOrderTermsSrc(3,3), traceNCP, traceSrc, dtphi, dtTraceK, dtP(3)
    REAL :: RNCP, RSrc, RPlusTwoNablaZNCP, RPlusTwoNablaZSrc, nablanablaalpha, nablanablaalphaNCP, nablanablaalphaSrc, Riemann(3,3,3,3), dChristoffel(3,3,3,3)
    REAL :: TwoNablaZNCP, TwoNablaZSrc,  divAupNCP(3), divAupSrc(3), XX(3), dXX(3,3), nablaXNCP(3,3), nablaXSrc(3,3), dtX(3)
    ! Matter contributions
    REAL :: sm(3), Sij(3,3), Sij_contr(3,3), sm_contr(3), S, tau, dens, bv_cov(3), sqrdet
    REAL :: srctraceK, srcTheta
    REAL :: SrcK(3,3), SrcGhat(3),mytmp1(3,3,3,3), mytmp2(3,3,3,3)

    REAL, PARAMETER :: Pi = ACOS(-1.0)
    !
    BgradQ = 0.0
    !
#if defined(CCZ4EINSTEIN) || defined(CCZ4GRHD) || defined(CCZ4GRMHD) || defined(CCZ4GRGPR)
    !
    Qx = gradQin(:,1)
    Qy = gradQin(:,2)
    IF(nDim.eq.2) THEN
        Qz = 0.0
        gradQ(:,1:2)=gradQin(:,1:2)
        gradQ(:,3)	=0.0
    ELSE
        Qz = gradQin(:,3)
        gradQ=gradQin
    ENDIF

    k1   = EQN%CCZ4k1
    k2   = EQN%CCZ4k2
    k3   = EQN%CCZ4k3
    fff  = EQN%CCZ4f
    ggg  = EQN%CCZ4g
    eta  = EQN%CCZ4eta
    itau = EQN%CCZ4itau
    e    = EQN%CCZ4e
    c    = EQN%CCZ4c
    mu   = EQN%CCZ4mu
    ds   = EQN%CCZ4ds
    bs   = EQN%CCZ4bs
    xi   = EQN%CCZ4xi
    sk   = EQN%CCZ4sk
    !
    ! These are the tilde quantities, so be careful !
    g_cov(1,1) = Q(1)
    g_cov(1,2) = Q(2)
    g_cov(1,3) = Q(3)
    g_cov(2,1) = Q(2)
    g_cov(2,2) = Q(4)
    g_cov(2,3) = Q(5)
    g_cov(3,1) = Q(3)
    g_cov(3,2) = Q(5)
    g_cov(3,3) = Q(6)
    ! This determinant should be unity, since we use the conformal decomposition
    det = (Q(1)*Q(4)*Q(6)-Q(1)*Q(5)**2-Q(2)**2*Q(6)+2*Q(2)*Q(3)*Q(5)-Q(3)**2*Q(4))

    g_contr(1,1) =  (g_cov(2,2)*g_cov(3,3)-g_cov(2,3)*g_cov(3,2)) / det
    g_contr(1,2) = -(g_cov(1,2)*g_cov(3,3)-g_cov(1,3)*g_cov(3,2)) / det
    g_contr(1,3) = -(-g_cov(1,2)*g_cov(2,3)+g_cov(1,3)*g_cov(2,2))/ det
    g_contr(2,1) = -(g_cov(2,1)*g_cov(3,3)-g_cov(2,3)*g_cov(3,1)) / det
    g_contr(2,2) = (g_cov(1,1)*g_cov(3,3)-g_cov(1,3)*g_cov(3,1))  / det
    g_contr(2,3) = -(g_cov(1,1)*g_cov(2,3)-g_cov(1,3)*g_cov(2,1)) / det
    g_contr(3,1) = -(-g_cov(2,1)*g_cov(3,2)+g_cov(2,2)*g_cov(3,1))/ det
    g_contr(3,2) = -(g_cov(1,1)*g_cov(3,2)-g_cov(1,2)*g_cov(3,1)) / det
    g_contr(3,3) = (g_cov(1,1)*g_cov(2,2)-g_cov(1,2)*g_cov(2,1))  / det

    alpha = EXP(MAX(-20.,MIN(20.,Q(17))))
    SELECT CASE(EQN%CCZ4LapseType)
    CASE(0)  ! harmonic
        fa = 1.0
        faa = 0.0
    CASE DEFAULT  ! 1 + log
        fa = 2.0/alpha
        faa = -2.0/alpha**2
    END SELECT
    !
    K0    = Q(59)
    dK0   = sk*gradQ(59,:)
    !
    Aex(1,1) = Q(7)
    Aex(1,2) = Q(8)
    Aex(1,3) = Q(9)
    Aex(2,1) = Q(8)
    Aex(2,2) = Q(10)
    Aex(2,3) = Q(11)
    Aex(3,1) = Q(9)
    Aex(3,2) = Q(11)
    Aex(3,3) = Q(12)
    !
    traceA = SUM(g_contr*Aex)
    Aex = Aex - 1./3.*g_cov*traceA
    !
    dAex(:,1,1) = gradQ(7,:)
    dAex(:,1,2) = gradQ(8,:)
    dAex(:,1,3) = gradQ(9,:)
    dAex(:,2,1) = gradQ(8,:)
    dAex(:,2,2) = gradQ(10,:)
    dAex(:,2,3) = gradQ(11,:)
    dAex(:,3,1) = gradQ(9,:)
    dAex(:,3,2) = gradQ(11,:)
    dAex(:,3,3) = gradQ(12,:)
    !
    Amix = MATMUL(g_contr, Aex)
    Aup  = MATMUL(g_contr, TRANSPOSE(Amix))
    !
    Theta = Q(13)
    dTheta = gradQ(13,:)
    !
    Ghat = (/ Q(14), Q(15), Q(16) /)
    dGhat(:,1) = gradQ(14,:)
    dGhat(:,2) = gradQ(15,:)
    dGhat(:,3) = gradQ(16,:)
    !
    b = Q(21:23)
    !
    AA    = (/ Q(24), Q(25), Q(26) /)
    dAA(:,1) = gradQ(24,:)
    dAA(:,2) = gradQ(25,:)
    dAA(:,3) = gradQ(26,:)
    !
    traceK = Q(54)
    dtraceK = gradQ(54,:)
    !
    phi   =  EXP(MAX(-20.,MIN(20.,Q(55))))

    PP    = Q(56:58)
    dPP(:,1) = gradQ(56,:)
    dPP(:,2) = gradQ(57,:)
    dPP(:,3) = gradQ(58,:)
    !
    beta = (/ Q(18), Q(19), Q(20) /)
    BB(1,1) = Q(27)
    BB(2,1) = Q(28)
    BB(3,1) = Q(29)
    BB(1,2) = Q(30)
    BB(2,2) = Q(31)
    BB(3,2) = Q(32)
    BB(1,3) = Q(33)
    BB(2,3) = Q(34)
    BB(3,3) = Q(35)
    !
    dBB(:,1,1) = gradQ(27,:)
    dBB(:,2,1) = gradQ(28,:)
    dBB(:,3,1) = gradQ(29,:)
    dBB(:,1,2) = gradQ(30,:)
    dBB(:,2,2) = gradQ(31,:)
    dBB(:,3,2) = gradQ(32,:)
    dBB(:,1,3) = gradQ(33,:)
    dBB(:,2,3) = gradQ(34,:)
    dBB(:,3,3) = gradQ(35,:)
    !
    dBB = dBB*sk
    !
    DD(1,1,1)=Q(36)
    DD(1,1,2)=Q(37)
    DD(1,1,3)=Q(38)
    DD(1,2,1)=Q(37)
    DD(1,2,2)=Q(39)
    DD(1,2,3)=Q(40)
    DD(1,3,1)=Q(38)
    DD(1,3,2)=Q(40)
    DD(1,3,3)=Q(41)
    !
    DD(2,1,1)=Q(42)
    DD(2,1,2)=Q(43)
    DD(2,1,3)=Q(44)
    DD(2,2,1)=Q(43)
    DD(2,2,2)=Q(45)
    DD(2,2,3)=Q(46)
    DD(2,3,1)=Q(44)
    DD(2,3,2)=Q(46)
    DD(2,3,3)=Q(47)
    !
    DD(3,1,1)=Q(48)
    DD(3,1,2)=Q(49)
    DD(3,1,3)=Q(50)
    DD(3,2,1)=Q(49)
    DD(3,2,2)=Q(51)
    DD(3,2,3)=Q(52)
    DD(3,3,1)=Q(50)
    DD(3,3,2)=Q(52)
    DD(3,3,3)=Q(53)
    !
    dDD(:,1,1,1)=gradQ(36,:)
    dDD(:,1,1,2)=gradQ(37,:)
    dDD(:,1,1,3)=gradQ(38,:)
    dDD(:,1,2,1)=gradQ(37,:)
    dDD(:,1,2,2)=gradQ(39,:)
    dDD(:,1,2,3)=gradQ(40,:)
    dDD(:,1,3,1)=gradQ(38,:)
    dDD(:,1,3,2)=gradQ(40,:)
    dDD(:,1,3,3)=gradQ(41,:)
    dDD(:,2,1,1)=gradQ(42,:)
    dDD(:,2,1,2)=gradQ(43,:)
    dDD(:,2,1,3)=gradQ(44,:)
    dDD(:,2,2,1)=gradQ(43,:)
    dDD(:,2,2,2)=gradQ(45,:)
    dDD(:,2,2,3)=gradQ(46,:)
    dDD(:,2,3,1)=gradQ(44,:)
    dDD(:,2,3,2)=gradQ(46,:)
    dDD(:,2,3,3)=gradQ(47,:)
    dDD(:,3,1,1)=gradQ(48,:)
    dDD(:,3,1,2)=gradQ(49,:)
    dDD(:,3,1,3)=gradQ(50,:)
    dDD(:,3,2,1)=gradQ(49,:)
    dDD(:,3,2,2)=gradQ(51,:)
    dDD(:,3,2,3)=gradQ(52,:)
    dDD(:,3,3,1)=gradQ(50,:)
    dDD(:,3,3,2)=gradQ(52,:)
    dDD(:,3,3,3)=gradQ(53,:)
    !
    dgup = 0.0
    DO n = 1, 3
     DO j = 1, 3
      DO l = 1, 3
       DO m = 1, 3
        DO k = 1, 3
           dgup(k,m,l) = dgup(k,m,l)-g_contr(m,n)*g_contr(j,l)*2*DD(k,n,j)
        ENDDO
       ENDDO
      ENDDO
     ENDDO
    ENDDO
    !
    Kex  = Aex/phi**2 + 1./3.*traceK*g_cov/phi**2
    Kmix = MATMUL( phi**2*g_contr, Kex  )
    Kup  = MATMUL( phi**2*g_contr, TRANSPOSE(Kmix) )
    !
    Christoffel_tilde = 0.0
    Christoffel       = 0.0
    Gtilde = 0.0
    !
    DO k = 1, 3
     DO j = 1, 3
      DO i = 1, 3
       Christoffel_kind1(i,j,k) = DD(k,i,j)+DD(j,i,k)-DD(i,j,k)      ! this definition seems to work !
       DO l = 1, 3
          Christoffel_tilde(i,j,k) = Christoffel_tilde(i,j,k) + g_contr(k,l)*( DD(i,j,l)+DD(j,i,l)-DD(l,i,j) )
          Christoffel(i,j,k)       = Christoffel(i,j,k)       + g_contr(k,l)*(DD(i,j,l)+DD(j,i,l)-DD(l,i,j))-g_contr(k,l)*( g_cov(j,l)*PP(i)+g_cov(i,l)*PP(j)-g_cov(i,j)*PP(l) )
        ENDDO
      ENDDO
     ENDDO
    ENDDO

    DO l = 1, 3
     DO j = 1, 3
      DO i = 1, 3
          Gtilde(i) = Gtilde(i) + g_contr(j,l)*Christoffel_tilde(j,l,i)
      ENDDO
     ENDDO
    ENDDO
    !
    Z   = ds*0.5*MATMUL( g_cov, Ghat - Gtilde )
    Zup = MATMUL(phi**2*g_contr, Z)
    !
    dChristoffelNCP = 0.0
    dChristoffelSrc = 0.0
    dChristoffel_tildeNCP = 0.0
    dChristoffel_tildeSrc = 0.0
    DO l = 1, 3
     DO m = 1, 3
      DO j = 1, 3
       DO i = 1, 3
        DO k = 1, 3
            dChristoffelNCP(k,i,j,m) = dChristoffelNCP(k,i,j,m) + g_contr(m,l)*( 0.5*(dDD(k,i,j,l)+dDD(i,k,j,l))+0.5*(dDD(k,j,i,l)+dDD(j,k,i,l))-0.5*(dDD(k,l,i,j)+dDD(l,k,i,j)) )         &
                                                                  - g_contr(m,l)*( g_cov(j,l)*0.5*(dPP(k,i)+dPP(i,k))+g_cov(i,l)*0.5*(dPP(k,j)+dPP(j,k))-g_cov(i,j)*0.5*(dPP(k,l)+dPP(l,k)) )
            !
            dChristoffel_tildeNCP(k,i,j,m) = dChristoffel_tildeNCP(k,i,j,m) + g_contr(m,l)*( 0.5*(dDD(k,i,j,l)+dDD(i,k,j,l))+0.5*(dDD(k,j,i,l)+dDD(j,k,i,l))-0.5*(dDD(k,l,i,j)+dDD(l,k,i,j)) )
            !
            dChristoffelSrc(k,i,j,m) = dChristoffelSrc(k,i,j,m) + dgup(k,m,l)*(DD(i,j,l)+DD(j,i,l)-DD(l,i,j)) - dgup(k,m,l)*(g_cov(j,l)*PP(i)+g_cov(i,l)*PP(j)-g_cov(i,j)*PP(l)) - g_contr(m,l)*( 2*DD(k,j,l)*PP(i)+2*DD(k,i,l)*PP(j)-2*DD(k,i,j)*PP(l) )
            !
            dChristoffel_tildeSrc(k,i,j,m) = dChristoffel_tildeSrc(k,i,j,m) + dgup(k,m,l)*(DD(i,j,l)+DD(j,i,l)-DD(l,i,j))
          ENDDO
       ENDDO
      ENDDO
     ENDDO
    ENDDO
    !
    RiemannSrc = 0.0
    RiemannNCP = 0.0
    DO m = 1, 3
     DO j = 1, 3
      DO k = 1, 3
       DO i = 1, 3
          RiemannNCP(i,k,j,m) = dChristoffelNCP(k,i,j,m)-dChristoffelNCP(j,i,k,m)
          RiemannSrc(i,k,j,m) = dChristoffelSrc(k,i,j,m)-dChristoffelSrc(j,i,k,m)
          DO l = 1, 3
           RiemannSrc(i,k,j,m) = RiemannSrc(i,k,j,m) + Christoffel(i,j,l)*Christoffel(l,k,m) - Christoffel(i,k,l)*Christoffel(l,j,m)
          ENDDO
       ENDDO
      ENDDO
     ENDDO
    ENDDO
    !
    RicciNCP = 0.0
    RicciSrc = 0.0
    DO l = 1, 3
     DO n = 1, 3
      DO m = 1, 3
         RicciNCP(m,n) = RicciNCP(m,n) + RiemannNCP(m,l,n,l)
         RicciSrc(m,n) = RicciSrc(m,n) + RiemannSrc(m,l,n,l)
      ENDDO
     ENDDO
    ENDDO
    !
    RNCP = phi**2*SUM(g_contr*RicciNCP)
    RSrc = phi**2*SUM(g_contr*RicciSrc)
    !
    ! Here we directly compute the derivative of Gtilde from its original definition as contracted Christoffel symbol,
    ! without assuming unit determinant of the conformal metric. Back to the roots, and as few assumptions as possible...
    !
    dGtildeNCP = 0.0
    dGtildeSrc = 0.0
    DO l = 1, 3
     DO j = 1, 3
      DO i = 1, 3
       DO k = 1, 3
           dGtildeSrc(k,i) = dGtildeSrc(k,i) + dgup(k,j,l)*Christoffel_tilde(j,l,i) + g_contr(j,l)*dChristoffel_tildeSrc(k,j,l,i)
           dGtildeNCP(k,i) = dGtildeNCP(k,i)                                        + g_contr(j,l)*dChristoffel_tildeNCP(k,j,l,i)
       ENDDO
      ENDDO
     ENDDO
    ENDDO
    !
    dZNCP = 0.0
    dZSrc = 0.0
    DO j = 1, 3
     DO i = 1, 3
      DO k = 1, 3
        dZNCP(k,i) = dZNCP(k,i) + ds*0.5*g_cov(i,j)*( dGhat(k,j)-dGtildeNCP(k,j) )
        dZSrc(k,i) = dZSrc(k,i) + ds*(DD(k,i,j)*(Ghat(j)-Gtilde(j)) + 0.5*g_cov(i,j)*( -dGtildeSrc(k,j) ) )
       ENDDO
      ENDDO
    ENDDO
    !
    nablaZNCP = dZNCP
    nablaZSrc = 0.0
    DO j = 1, 3
     DO i = 1, 3
      nablaZSrc(i,j) = dZSrc(i,j)
      DO k = 1, 3
        nablaZSrc(i,j) = nablaZSrc(i,j) - Christoffel(i,j,k)*Z(k)
      ENDDO
     ENDDO
    ENDDO
    !
    RicciPlusNablaZNCP = RicciNCP + ( nablaZNCP + TRANSPOSE(nablaZNCP) )
    RicciPlusNablaZSrc = RicciSrc + ( nablaZSrc + TRANSPOSE(nablaZSrc) )
    !
    RPlusTwoNablaZNCP = phi**2*SUM(g_contr*RicciPlusNablaZNCP)
    RPlusTwoNablaZSrc = phi**2*SUM(g_contr*RicciPlusNablaZSrc)
    !
    nablaijalphaNCP = 0.0
    nablaijalphaSrc = 0.0
    DO j = 1, 3
     DO i = 1, 3
       nablaijalphaNCP(i,j) = alpha*0.5*( dAA(i,j)+dAA(j,i) )
       nablaijalphaSrc(i,j) = alpha*AA(j)*AA(i)
       DO k = 1, 3
         nablaijalphaSrc(i,j) = nablaijalphaSrc(i,j) - Christoffel(i,j,k)*alpha*AA(k)
       ENDDO
     ENDDO
    ENDDO
    nablanablaalphaNCP = phi**2*SUM( g_contr*nablaijalphaNCP )
    nablanablaalphaSrc = phi**2*SUM( g_contr*nablaijalphaSrc )
    !
    SecondOrderTermsNCP = -nablaijalphaNCP + alpha*RicciPlusNablaZNCP
    SecondOrderTermsSrc = -nablaijalphaSrc + alpha*RicciPlusNablaZSrc
    traceNCP = SUM( g_contr*SecondOrderTermsNCP )
    SecondOrderTermsNCP = SecondOrderTermsNCP - 1./3.*g_cov*traceNCP
    traceSrc = SUM( g_contr*SecondOrderTermsSrc )
    SecondOrderTermsSrc = SecondOrderTermsSrc - 1./3.*g_cov*traceSrc
    !
    ! Now assemble all this terrible stuff...
    !
    dtgamma = - 2*alpha*Aex - itau*(det-1.0)*g_cov
    DO k = 1, 3
     DO j = 1, 3
      DO i = 1, 3
          dtgamma(i,j) = dtgamma(i,j) + g_cov(k,i)*BB(j,k) + g_cov(k,j)*BB(i,k) - 2./3.*g_cov(i,j)*BB(k,k) + beta(k)*2*DD(k,i,j)
      ENDDO
     ENDDO
    ENDDO
    !
    ! Main variables of the CCZ4 system
    dtK = phi**2*SecondOrderTermsNCP + beta(1)*dAex(1,:,:) + beta(2)*dAex(2,:,:) + beta(3)*dAex(3,:,:)      ! extrinsic curvature
    dtK = dtK + phi**2*SecondOrderTermsSrc + alpha*Aex*(traceK-2*Theta) - 2*alpha*MATMUL(Aex,Amix) - itau*g_cov*traceA
    DO j = 1, 3
     DO i = 1, 3
      DO k = 1, 3
         dtK(i,j) = dtK(i,j) + Aex(k,i)*BB(j,k) + Aex(k,j)*BB(i,k) - 2./3.*Aex(i,j)*BB(k,k)
      ENDDO
     ENDDO
    ENDDO
    !
    dtTraceK = - nablanablaalphaNCP - nablanablaalphaSrc + alpha*( RPlusTwoNablaZNCP + RPlusTwoNablaZSrc + traceK**2 - 2*Theta*traceK ) - 3*alpha*k1*(1+k2)*Theta + SUM(beta(:)*dtraceK(:))
    !
    traceB = BB(1,1) + BB(2,2) + BB(3,3)
    dtphi   = beta(1)*PP(1) + beta(2)*PP(2) + beta(3)*PP(3) + 1./3.*alpha*traceK - 1./3.*traceB
    dtalpha = -alpha*fa*(traceK-K0-c*2*Theta) + beta(1)*AA(1) + beta(2)*AA(2) + beta(3)*AA(3)

    Aupdown = SUM(Aex*Aup)
    ! *** original
    dtTheta =  0.5*alpha*e**2*(RplusTwoNablaZNCP + RplusTwoNablaZSrc) + beta(1)*dTheta(1) + beta(2)*dTheta(2) + beta(3)*dTheta(3)    &            ! temporal Z
             + 0.5*alpha*e**2*( - Aupdown + 2./3.*traceK**2 ) - alpha*Theta*traceK - SUM(Zup*alpha*AA) - alpha*k1*(2+k2)*Theta
    !
    divAupNCP = 0.0
    divAupSrc = 0.0
    DO k = 1, 3
        DO j = 1, 3
         DO l = 1, 3
          DO i = 1, 3
            divAupNCP(i) = divAupNCP(i) + g_contr(i,l)*g_contr(j,k)*dAex(j,l,k)
            divAupSrc(i) = divAupSrc(i) + ( dgup(j,i,l)*g_contr(j,k) + g_contr(i,l)*dgup(j,j,k) )*Aex(l,k)
          ENDDO
         ENDDO
        ENDDO
    ENDDO
    !
    dtGhat = 0.0
    DO i = 1, 3
          dtGhat(i) = dtGhat(i) +        &
                      + 2*alpha*( SUM(Christoffel_tilde(:,:,i)*Aup(:,:)) - 3.0*SUM(Aup(i,:)*PP(:)) - 2./3.*SUM(g_contr(i,:)*dtraceK(:)) )    &
                      + 2*alpha*SUM( g_contr(:,i)*( dTheta(:) - Theta*AA(:) - 2./3.*traceK*Z(:)  ) )  &
                      - 2*SUM( Aup(i,:)*alpha*AA(:) ) - 2*alpha*k1*SUM(g_contr(i,:)*Z(:)) - SUM(Gtilde(:)*BB(:,i))   &
                    + beta(1)*dGhat(1,i) + beta(2)*dGhat(2,i) + beta(3)*dGhat(3,i) + 2./3.*Gtilde(i)*traceB
        DO l = 1, 3
         DO k = 1, 3
             dtGhat(i) = dtGhat(i) + g_contr(k,l)*0.5*(dBB(k,l,i)+dBB(l,k,i)) + 1./3*g_contr(i,k)*0.5*(dBB(k,l,l)+dBB(l,k,l)) + &
                                     2*k3*( 2./3.*g_contr(i,l)*Z(l)*BB(k,k) - g_contr(l,k)*Z(l)*BB(k,i) )
         ENDDO
        ENDDO
    ENDDO
    DO k = 1, 3
        ov(k) = + 2*alpha*( SUM(g_contr(:,:)*dAex(k,:,:)) + SUM(dgup(k,:,:)*Aex(:,:)) )   ! here we can use the constraint that trace A tilde = 0.
    ENDDO
    dtGhat = dtGhat + sk*MATMUL(g_contr,ov)                                               ! the above ordering constraint is "down", so we must raise the index via g_contr.
    !
    dtbb = xi*dtGhat - eta*b                                                              !  <= be careful, this damping term -eta*b may be dangerous for the gamma driver, since it may kill the waves that you want !
    ! Add the following terms if you want shift convection in the PDE for b^i
    dtbb = dtbb + bs*( - beta(1)*gradQ(14:16,1) - beta(2)*gradQ(14:16,2) - beta(3)*gradQ(14:16,3)  + beta(1)*gradQ(21:23,1) + beta(2)*gradQ(21:23,2) + beta(3)*gradQ(21:23,3)  )   !
    dtbb = sk*dtbb
    !
    dtbeta  = + fff*b
    ! Add the following term if you want to have shift convection in the PDE for beta^i
    ! Do not add it if you want a real Lie derivative for beta. In this case, the advection term cancels out.
    dtbeta = dtbeta + bs*( beta(1)*BB(1,:) + beta(2)*BB(2,:) + beta(3)*BB(3,:) )
    dtbeta = sk*dtbeta
    !
    ! Auxiliary variables
    dtA = -alpha*fa*( dtraceK(:) -dK0(:) - c*2*dTheta(:) ) + beta(1)*dAA(1,:) + beta(2)*dAA(2,:) + beta(3)*dAA(3,:) - alpha*AA*(fa+alpha*faa)*(traceK-K0-c*2*Theta) + MATMUL(BB, AA)
    DO k = 1, 3
        dtA(k) = dtA(k) - sk*alpha*fa*( SUM(g_contr(:,:)*dAex(k,:,:)) + SUM(dgup(k,:,:)*Aex(:,:)) )   ! here we can use the constraint that trace A tilde = 0.
    ENDDO
    !
    ! In CCZ4 we have completely removed all the conservative fluxes.
    dtB(:,1) = fff*gradQ(21,:)
    dtB(:,2) = fff*gradQ(22,:)
    dtB(:,3) = fff*gradQ(23,:)
    ! #ordB2#
    ! for the ordering constraints, we have to check whether they should be multiplied by alpha**2, or not...
    DO j = 1, 3
     DO i = 1, 3
      DO k = 1, 3
         dtB(k,i) = dtB(k,i) + mu*alpha**2*g_contr(i,j)*( (dPP(k,j)-dPP(j,k)) )
         DO n = 1, 3
          DO l = 1, 3
            dtB(k,i) = dtB(k,i) - mu*alpha**2*g_contr(i,j)*g_contr(n,l)*( dDD(k,l,j,n)-dDD(l,k,j,n) )
          ENDDO
         ENDDO
      ENDDO
      !
      ENDDO
    ENDDO
    !
    dtB = dtB + bs*( beta(1)*dBB(1,:,:) + beta(2)*dBB(2,:,:) + beta(3)*dBB(3,:,:) + MATMUL(BB,BB) )
    dtB = dtB*sk
    !
    dtD = -alpha*dAex
    DO m = 1, 3
     DO j = 1, 3
      DO i = 1, 3
        DO k = 1, 3
            dtD(k,i,j) = dtD(k,i,j) + ( 0.5*(g_cov(m,i)*0.5*(dBB(k,j,m)+dBB(j,k,m))+g_cov(m,j)*0.5*(dBB(k,i,m)+dBB(i,k,m)) ) - 1./3.*g_cov(i,j)*0.5*(dBB(k,m,m)+dBB(m,k,m)) )
            DO n = 1, 3
                dtD(k,i,j) = dtD(k,i,j) + 1./3*alpha*g_cov(i,j)*g_contr(n,m)*dAex(k,n,m) + 1./3.*alpha*g_cov(i,j)*dgup(k,n,m)*Aex(n,m)   ! explicitly remove the trace of tilde A again
            ENDDO
        ENDDO
       ENDDO
      ENDDO
    ENDDO
    !
    DO j = 1, 3
     DO i = 1, 3
      DO k = 1, 3
        dtD(k,i,j) = dtD(k,i,j) - alpha*AA(k)*Aex(i,j)
        DO l = 1, 3
          dtD(k,i,j) = dtD(k,i,j) + BB(k,l)*DD(l,i,j) + DD(k,l,i)*BB(j,l) + DD(k,l,j)*BB(i,l) - 2./3.*DD(k,i,j)*BB(l,l)
        ENDDO
      ENDDO
     ENDDO
    ENDDO
    !
    dtD = dtD + beta(1)*dDD(1,:,:,:) + beta(2)*dDD(2,:,:,:) + beta(3)*dDD(3,:,:,:)
    !
    dtP = MATMUL(BB,PP) + beta(1)*dPP(1,:) + beta(2)*dPP(2,:) + beta(3)*dPP(3,:)
    DO k = 1, 3
     dtP(k) = dtP(k) + 1./3.*alpha*dtraceK(k) + 1./3.*alpha*AA(k)*traceK + sk*1./3.*alpha*( SUM(g_contr(:,:)*dAex(k,:,:)) + SUM(dgup(k,:,:)*Aex(:,:)) )
     DO i = 1, 3
          dtP(k) = dtP(k) - 1./3.*0.5*(dBB(k,i,i)+dBB(i,k,i))
     ENDDO
    ENDDO
    !
    BgradQ(1:6)    = (/ dtgamma(1,1), dtgamma(1,2), dtgamma(1,3), dtgamma(2,2), dtgamma(2,3), dtgamma(3,3) /)          ! \tilde \gamma_ij
    BgradQ(7:12)   = (/ dtK(1,1), dtK(1,2), dtK(1,3), dtK(2,2), dtK(2,3), dtK(3,3) /)                                  ! \tilde A_ij
    BgradQ(13)     = dtTheta                                                                                ! Theta
    BgradQ(14:16)  = dtGhat(1:3)                                                                                ! \hat \Gamma^i
    BgradQ(17)     = dtalpha                                                                                ! log alpha
    BgradQ(18:20)  = dtbeta                                                                                ! beta^i
    BgradQ(21:23)  = dtbb                                                                                ! b^i
    BgradQ(24:26)  = dtA(1:3)                                                                                ! A_k
    BgradQ(27:35)  = (/ dtB(1,1), dtB(2,1), dtB(3,1), dtB(1,2), dtB(2,2), dtB(3,2), dtB(1,3), dtB(2,3), dtB(3,3) /)    ! B_k^i
    BgradQ(36:41)  = (/ dtD(1,1,1), dtD(1,1,2), dtD(1,1,3), dtD(1,2,2), dtD(1,2,3), dtD(1,3,3) /)                      ! D_kij
    BgradQ(42:47)  = (/ dtD(2,1,1), dtD(2,1,2), dtD(2,1,3), dtD(2,2,2), dtD(2,2,3), dtD(2,3,3) /)                      ! D_kij
    BgradQ(48:53)  = (/ dtD(3,1,1), dtD(3,1,2), dtD(3,1,3), dtD(3,2,2), dtD(3,2,3), dtD(3,3,3) /)                      ! D_kij
    BgradQ(54)     = dtTraceK                                                                                ! traceK
    BgradQ(55)     = dtphi                                                                                ! log phi
    BgradQ(56:58)  = dtP                                                                                ! P_k
    BgradQ(59)     = 0
    !
    Src_BgradQ = BgradQ    ! here, we do not have to change sign, since we work on the right hand side in the fused subroutine
    !
    !
    RETURN
    !
#endif
    !
END SUBROUTINE PDEFusedSrcNCP


text:
! addw0 adds or subtracts the base state velocity, w0, to the MAC
! velocities (which are originally constructed from Utilde).

module addw0_module

  use bl_types, only: dp_t
  use box_module, only: upb, lwb
  use multifab_module, only: multifab, multifab_fill_boundary, &
                             nboxes, nghost, dataptr, get_box, nfabs

  implicit none

  private

  public :: addw0, addw0_3d_sphr

contains

  subroutine addw0(umac,the_bc_level,mla,w0,w0mac,mult)

    use bl_prof_module, only: bl_prof_timer, build, destroy
    use create_umac_grown_module, only: create_umac_grown
    use geometry, only: spherical
    use define_bc_module, only: bc_level
    use ml_layout_module, only: ml_layout
    use ml_cc_restriction_module, only: ml_edge_restriction
    use multifab_physbc_edgevel_module, only: multifab_physbc_edgevel

    type(multifab) , intent(inout) :: umac(:,:)
    type(bc_level) , intent(in   ) :: the_bc_level(:)
    type(ml_layout), intent(inout) :: mla
    real(kind=dp_t), intent(in   ) :: w0(:,0:)
    type(multifab) , intent(in   ) :: w0mac(:,:)
    real(kind=dp_t), intent(in   ) :: mult

    ! Local variables
    integer :: i,lo(mla%dim),hi(mla%dim),n,ng_um,ng_w0,dm,nlevs
    real(kind=dp_t), pointer :: ump(:,:,:,:)
    real(kind=dp_t), pointer :: vmp(:,:,:,:)
    real(kind=dp_t), pointer :: wmp(:,:,:,:)
    real(kind=dp_t), pointer :: w0xp(:,:,:,:)
    real(kind=dp_t), pointer :: w0yp(:,:,:,:)
    real(kind=dp_t), pointer :: w0zp(:,:,:,:)

    type(bl_prof_timer), save :: bpt

    call build(bpt, "addw0")

    dm = mla%dim
    nlevs = mla%nlevel

    ng_um = nghost(umac(1,1)) ! here we are assuming all components have the
                              ! same # of ghostcells

    do n = 1, nlevs
       do i = 1, nfabs(umac(n,dm))
          wmp  => dataptr(umac(n,dm), i)
          lo =  lwb(get_box(umac(n,dm), i))
          hi =  upb(get_box(umac(n,dm), i))
          select case(dm)
          case(1)
             call addw0_1d(wmp(:,1,1,1),ng_um,w0(n,:),lo,hi,mult)
          case(2)
             call addw0_2d(wmp(:,:,1,1),ng_um,w0(n,:),lo,hi,mult)
          case(3)
             if (spherical .eq. 0) then
                call addw0_3d(wmp(:,:,:,1),ng_um,w0(n,:),lo,hi,mult)
             else
                ump  => dataptr(umac(n,1), i)
                vmp  => dataptr(umac(n,2), i)
                wmp  => dataptr(umac(n,3), i)
                w0xp  => dataptr(w0mac(n,1), i)
                w0yp  => dataptr(w0mac(n,2), i)
                w0zp  => dataptr(w0mac(n,3), i)
                ng_w0 = nghost(w0mac(n,1))
                call addw0_3d_sphr(ump(:,:,:,1),vmp(:,:,:,1),wmp(:,:,:,1),ng_um, &
                                   w0xp(:,:,:,1),w0yp(:,:,:,1),w0zp(:,:,:,1),ng_w0, &
                                   lo,hi,mult)
             end if
          end select
       end do
    end do

    if (nlevs .eq. 1) then

       ! fill ghost cells for two adjacent grids at the same level
       ! this includes periodic domain boundary ghost cells
      do i=1,dm
          call multifab_fill_boundary(umac(nlevs,i))
       enddo

       ! fill non-periodic domain boundary ghost cells
       call multifab_physbc_edgevel(umac(nlevs,:),the_bc_level(nlevs))

    else

       ! the loop over nlevs must count backwards to make sure the finer grids are done first
       do n=nlevs,2,-1
          ! set level n-1 data to be the average of the level n data covering it
          do i=1,dm
             call ml_edge_restriction(umac(n-1,i),umac(n,i),mla%mba%rr(n-1,:),i)
          enddo
       end do

       do n=2,nlevs
          ! fill level n ghost cells using interpolation from level n-1 data
          ! note that multifab_fill_boundary and multifab_physbc_edgevel are called.
          call create_umac_grown(umac(n,:),umac(n-1,:),the_bc_level(n-1),the_bc_level(n), &
               n.eq.nlevs)
       end do

    end if

    call destroy(bpt)

  end subroutine addw0

  subroutine addw0_1d(umac,ng_um,w0,lo,hi,mult)

    integer        , intent(in   ) :: lo(:),hi(:),ng_um
    real(kind=dp_t), intent(inout) :: umac(lo(1)-ng_um:)
    real(kind=dp_t), intent(in   ) ::   w0(0:)
    real(kind=dp_t), intent(in   ) :: mult

    integer :: i

    do i = lo(1),hi(1)+1
       umac(i) = umac(i) + mult * w0(i)
    end do

  end subroutine addw0_1d

  subroutine addw0_2d(vmac,ng_um,w0,lo,hi,mult)

    integer        , intent(in   ) :: lo(:),hi(:),ng_um
    real(kind=dp_t), intent(inout) :: vmac(lo(1)-ng_um:,lo(2)-ng_um:)
    real(kind=dp_t), intent(in   ) ::   w0(0:)
    real(kind=dp_t), intent(in   ) :: mult

    integer :: i,j

    do j = lo(2),hi(2)+1
       do i = lo(1)-1,hi(1)+1
          vmac(i,j) = vmac(i,j) + mult * w0(j)
       end do
    end do

  end subroutine addw0_2d

  subroutine addw0_3d(wmac,ng_um,w0,lo,hi,mult)

    integer        , intent(in   ) :: lo(:),hi(:),ng_um
    real(kind=dp_t), intent(inout) :: wmac(lo(1)-ng_um:,lo(2)-ng_um:,lo(3)-ng_um:)
    real(kind=dp_t), intent(in   ) ::   w0(0:)
    real(kind=dp_t), intent(in   ) :: mult

    integer :: i,j,k

    !$OMP PARALLEL DO PRIVATE(i,j,k)
    do k = lo(3),hi(3)+1
       do j = lo(2)-1,hi(2)+1
          do i = lo(1)-1,hi(1)+1
             wmac(i,j,k) = wmac(i,j,k) + mult * w0(k)
          end do
       end do
    end do
    !$OMP END PARALLEL DO

  end subroutine addw0_3d

  subroutine addw0_3d_sphr(umac,vmac,wmac,ng_um,w0macx,w0macy,w0macz,ng_w0,lo,hi,mult)

    integer        , intent(in   ) :: lo(:),hi(:),ng_um,ng_w0
    real(kind=dp_t), intent(inout) ::    umac(lo(1)-ng_um:,lo(2)-ng_um:,lo(3)-ng_um:)
    real(kind=dp_t), intent(inout) ::    vmac(lo(1)-ng_um:,lo(2)-ng_um:,lo(3)-ng_um:)
    real(kind=dp_t), intent(inout) ::    wmac(lo(1)-ng_um:,lo(2)-ng_um:,lo(3)-ng_um:)
    real(kind=dp_t), intent(in   ) ::  w0macx(lo(1)-ng_w0:,lo(2)-ng_w0:,lo(3)-ng_w0:)
    real(kind=dp_t), intent(in   ) ::  w0macy(lo(1)-ng_w0:,lo(2)-ng_w0:,lo(3)-ng_w0:)
    real(kind=dp_t), intent(in   ) ::  w0macz(lo(1)-ng_w0:,lo(2)-ng_w0:,lo(3)-ng_w0:)
    real(kind=dp_t), intent(in   ) :: mult

    integer :: i,j,k

    !$OMP PARALLEL PRIVATE(i,j,k)
    !$OMP DO
    do k = lo(3),hi(3)
       do j = lo(2),hi(2)
          do i = lo(1),hi(1)+1
             umac(i,j,k) = umac(i,j,k) + mult * w0macx(i,j,k)
          end do
       end do
    end do
    !$OMP END DO NOWAIT
    !$OMP DO
    do k = lo(3),hi(3)
       do j = lo(2),hi(2)+1
          do i = lo(1),hi(1)
             vmac(i,j,k) = vmac(i,j,k) + mult * w0macy(i,j,k)
          end do
       end do
    end do
    !$OMP END DO NOWAIT
    !$OMP DO
    do k = lo(3),hi(3)+1
       do j = lo(2),hi(2)
          do i = lo(1),hi(1)
             wmac(i,j,k) = wmac(i,j,k) + mult * w0macz(i,j,k)
          end do
       end do
    end do
    !$OMP END DO
    !$OMP END PARALLEL

  end subroutine addw0_3d_sphr

end module addw0_module

max_stars_repo_name: sailoridy/MAESTRO

max_issues_repo_name: sailoridy/MAESTRO

max_forks_repo_name: sailoridy/MAESTRO


module m
  interface a
    subroutine s(x)
      real :: x
    end subroutine
    !ERROR: 's' is already declared in this scoping unit
    subroutine s(x)
      integer :: x
    end subroutine
  end interface
end module

module m2
  interface s
    subroutine s(x)
      real :: x
    end subroutine
    !ERROR: 's' is already declared in this scoping unit
    subroutine s(x)
      integer :: x
    end subroutine
  end interface
end module


!-------------------------------------------------------------------------------
! Copyright (c) 2016 The University of Tokyo
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
!> \brief Adaptive Mesh Refinement

!C
!C***
!C*** hecmw_adapt_NEW_CELL_341
!C***
!C
!C    create new TET.
!C
subroutine hecmw_adapt_NEW_CELL_341 (hecMESH, icouN)

  use hecmw_util
  implicit real*8 (A-H,O-Z)

  dimension NDIV(6), NNTYP(0:11)

  integer(kind=kint)        :: PAR_CEL_TYP
  type (hecmwST_local_mesh) :: hecMESH

  !C
  !C +----------------+
  !C | embedding TET. |
  !C +----------------+
  !C===
  do i= 0, 11
    NNTYP(i)= 0
  enddo

  do icel0= 1, hecMESH%n_adapt_act_elem_341
    icel_par= hecMESH%adapt_act_elem_341(icel0)
    NPAR= icel_par

    if (hecMESH%elem_ID(2*icel_par).eq.hecMESH%my_rank) then
      inc= 1
    else
      inc= 0
    endif

    is= hecMESH%elem_node_index(icel_par-1)
    n1= hecMESH%elem_node_item (is+1)
    n2= hecMESH%elem_node_item (is+2)
    n3= hecMESH%elem_node_item (is+3)
    n4= hecMESH%elem_node_item (is+4)
    n5= 0
    n6= 0
    n7= 0
    n8= 0
    n9= 0
    n0= 0

    nnn= hecMESH%n_adapt_edge
    call hecmw_adapt_EDGE_INFO ( hecMESH, n1, n2, ie1, 1 )
    call hecmw_adapt_EDGE_INFO ( hecMESH, n1, n3, ie2, 1 )
    call hecmw_adapt_EDGE_INFO ( hecMESH, n1, n4, ie3, 1 )
    call hecmw_adapt_EDGE_INFO ( hecMESH, n2, n3, ie4, 1 )
    call hecmw_adapt_EDGE_INFO ( hecMESH, n2, n4, ie5, 1 )
    call hecmw_adapt_EDGE_INFO ( hecMESH, n3, n4, ie6, 1 )

    NDIV(1)= hecMESH%adapt_iemb(ie1)
    NDIV(2)= hecMESH%adapt_iemb(ie2)
    NDIV(3)= hecMESH%adapt_iemb(ie3)
    NDIV(4)= hecMESH%adapt_iemb(ie4)
    NDIV(5)= hecMESH%adapt_iemb(ie5)
    NDIV(6)= hecMESH%adapt_iemb(ie6)

    NDIVSUM= NDIV(1)+NDIV(2)+NDIV(3)+NDIV(4)+NDIV(5)+NDIV(6)
    PAR_CEL_TYP= hecMESH%elem_type(NPAR)
    !C
    !C-- init. CHILD. cell array
    is= hecMESH%adapt_children_index(NPAR-1)
    iS1= is + 1
    iS2= is + 2
    iS3= is + 3
    iS4= is + 4
    iS5= is + 5
    iS6= is + 6
    iS7= is + 7
    iS8= is + 8

    hecMESH%adapt_children_item(2*iS1)= -1
    hecMESH%adapt_children_item(2*iS2)= -1
    hecMESH%adapt_children_item(2*iS3)= -1
    hecMESH%adapt_children_item(2*iS4)= -1
    hecMESH%adapt_children_item(2*iS5)= -1
    hecMESH%adapt_children_item(2*iS6)= -1
    hecMESH%adapt_children_item(2*iS7)= -1
    hecMESH%adapt_children_item(2*iS8)= -1

    hecMESH%adapt_children_item(2*iS1-1)= 0
    hecMESH%adapt_children_item(2*iS2-1)= 0
    hecMESH%adapt_children_item(2*iS3-1)= 0
    hecMESH%adapt_children_item(2*iS4-1)= 0
    hecMESH%adapt_children_item(2*iS5-1)= 0
    hecMESH%adapt_children_item(2*iS6-1)= 0
    hecMESH%adapt_children_item(2*iS7-1)= 0
    hecMESH%adapt_children_item(2*iS8-1)= 0

    !C
    !C== embedding TYPE
    if (NDIVSUM.eq.0) NTYP= 0
    !C
    !C-- TYP-1
    if (NDIVSUM.eq.1 .and. NDIV(1).eq.1) then
      NTYP= 1
      n5  = hecMESH%adapt_IWK(ie1)

      hecMESH%adapt_type(NPAR)= NTYP
      call hecmw_adapt_CREATE_NEW_TETRA (n1, n5, n3, n4, 1)
      call hecmw_adapt_CREATE_NEW_TETRA (n5, n2, n3, n4, 2)
    endif
    !C
    !C-- TYP-2
    if (NDIVSUM.eq.1 .and. NDIV(2).eq.1) then
      NTYP= 2
      n5  = hecMESH%adapt_IWK(ie2)

      hecMESH%adapt_type(NPAR)= NTYP
      call hecmw_adapt_CREATE_NEW_TETRA (n1, n2, n5, n4, 1)
      call hecmw_adapt_CREATE_NEW_TETRA (n5, n2, n3, n4, 2)
    endif
    !C
    !C-- TYP-3
    if (NDIVSUM.eq.1 .and. NDIV(3).eq.1) then
      NTYP= 3
      n5  = hecMESH%adapt_IWK(ie3)

      hecMESH%adapt_type(NPAR)= NTYP
      call hecmw_adapt_CREATE_NEW_TETRA (n1, n2, n3, n5, 1)
      call hecmw_adapt_CREATE_NEW_TETRA (n5, n2, n3, n4, 2)
    endif
    !C
    !C-- TYP-4
    if (NDIVSUM.eq.1 .and. NDIV(4).eq.1) then
      NTYP= 4
      n5  = hecMESH%adapt_IWK(ie4)

      hecMESH%adapt_type(NPAR)= NTYP
      call hecmw_adapt_CREATE_NEW_TETRA (n1, n2, n5, n4, 1)
      call hecmw_adapt_CREATE_NEW_TETRA (n1, n5, n3, n4, 2)
    endif
    !C
    !C-- TYP-5
    if (NDIVSUM.eq.1 .and. NDIV(5).eq.1) then
      NTYP= 5
      n5  = hecMESH%adapt_IWK(ie5)

      hecMESH%adapt_type(NPAR)= NTYP
      call hecmw_adapt_CREATE_NEW_TETRA (n1, n2, n3, n5, 1)
      call hecmw_adapt_CREATE_NEW_TETRA (n1, n5, n3, n4, 2)
    endif
    !C
    !C-- TYP-6
    if (NDIVSUM.eq.1 .and. NDIV(6).eq.1) then
      NTYP= 6
      n5  = hecMESH%adapt_IWK(ie6)

      hecMESH%adapt_type(NPAR)= NTYP
      call hecmw_adapt_CREATE_NEW_TETRA (n1, n2, n3, n5, 1)
      call hecmw_adapt_CREATE_NEW_TETRA (n1, n2, n5, n4, 2)
    endif
    !C
    !C-- TYP-7
    if (NDIVSUM.eq.3 .and. NDIV(1).eq.1 .and.                       &
        &      NDIV(3).eq.1 .and. NDIV(5).eq.1) then
      NTYP= 7
      n5  = hecMESH%adapt_IWK(ie1)
      n6  = hecMESH%adapt_IWK(ie3)
      n7  = hecMESH%adapt_IWK(ie5)

      hecMESH%adapt_type(NPAR)= NTYP
      call hecmw_adapt_CREATE_NEW_TETRA (n1, n5, n3, n6, 1)
      call hecmw_adapt_CREATE_NEW_TETRA (n5, n2, n3, n7, 2)
      call hecmw_adapt_CREATE_NEW_TETRA (n6, n7, n3, n4, 3)
      call hecmw_adapt_CREATE_NEW_TETRA (n6, n5, n3, n7, 4)
    endif
    !C
    !C-- TYP-8
    if (NDIVSUM.eq.3 .and. NDIV(2).eq.1 .and.                       &
        &      NDIV(3).eq.1 .and. NDIV(6).eq.1) then
      NTYP= 8
      n5  = hecMESH%adapt_IWK(ie2)
      n6  = hecMESH%adapt_IWK(ie3)
      n7  = hecMESH%adapt_IWK(ie6)

      hecMESH%adapt_type(NPAR)= NTYP
      call hecmw_adapt_CREATE_NEW_TETRA (n1, n2, n5, n6, 1)
      call hecmw_adapt_CREATE_NEW_TETRA (n5, n2, n3, n7, 2)
      call hecmw_adapt_CREATE_NEW_TETRA (n6, n2, n7, n4, 3)
      call hecmw_adapt_CREATE_NEW_TETRA (n6, n2, n5, n7, 4)
    endif
    !C
    !C-- TYP-9
    if (NDIVSUM.eq.3 .and. NDIV(1).eq.1 .and.                       &
        &      NDIV(2).eq.1 .and. NDIV(4).eq.1) then
      NTYP= 9
      n5  = hecMESH%adapt_IWK(ie1)
      n6  = hecMESH%adapt_IWK(ie2)
      n7  = hecMESH%adapt_IWK(ie4)

      hecMESH%adapt_type(NPAR)= NTYP
      call hecmw_adapt_CREATE_NEW_TETRA (n1, n5, n6, n4, 1)
      call hecmw_adapt_CREATE_NEW_TETRA (n5, n2, n7, n4, 2)
      call hecmw_adapt_CREATE_NEW_TETRA (n6, n7, n3, n4, 3)
      call hecmw_adapt_CREATE_NEW_TETRA (n5, n7, n6, n4, 4)
    endif
    !C
    !C-- TYP-10
    if (NDIVSUM.eq.3 .and. NDIV(4).eq.1 .and.                       &
        &      NDIV(5).eq.1 .and. NDIV(6).eq.1) then
      NTYP= 10
      n5  = hecMESH%adapt_IWK(ie4)
      n6  = hecMESH%adapt_IWK(ie5)
      n7  = hecMESH%adapt_IWK(ie6)

      hecMESH%adapt_type(NPAR)= NTYP
      call hecmw_adapt_CREATE_NEW_TETRA (n1, n2, n5, n6, 1)
      call hecmw_adapt_CREATE_NEW_TETRA (n1, n5, n3, n7, 2)
      call hecmw_adapt_CREATE_NEW_TETRA (n1, n6, n7, n4, 3)
      call hecmw_adapt_CREATE_NEW_TETRA (n1, n5, n7, n6, 4)
    endif
    !C
    !C-- TYP-11
    if (NDIVSUM.eq.6) then
      NTYP= 11
      n5  = hecMESH%adapt_IWK(ie1)
      n6  = hecMESH%adapt_IWK(ie2)
      n7  = hecMESH%adapt_IWK(ie3)
      n8  = hecMESH%adapt_IWK(ie4)
      n9  = hecMESH%adapt_IWK(ie5)
      n0  = hecMESH%adapt_IWK(ie6)

      hecMESH%adapt_type(NPAR)= NTYP
      call hecmw_adapt_CREATE_NEW_TETRA (n1, n5, n6, n7, 1)
      call hecmw_adapt_CREATE_NEW_TETRA (n5, n2, n8, n9, 2)
      call hecmw_adapt_CREATE_NEW_TETRA (n6, n8, n3, n0, 3)
      call hecmw_adapt_CREATE_NEW_TETRA (n7, n9, n0, n4, 4)
      call hecmw_adapt_CREATE_NEW_TETRA (n5, n8, n6, n7, 5)
      call hecmw_adapt_CREATE_NEW_TETRA (n5, n8, n7, n9, 6)
      call hecmw_adapt_CREATE_NEW_TETRA (n0, n8, n7, n6, 7)
      call hecmw_adapt_CREATE_NEW_TETRA (n0, n8, n9, n7, 8)
    endif
    !C==

    !C
    !C-- TYPE of EMBEDDING
    NNTYP             (NTYP)= NNTYP(NTYP) + 1
  enddo
  !C===

  return

contains
  subroutine hecmw_adapt_CREATE_NEW_TETRA (in1, in2, in3, in4, IDchi)

    hecMESH%n_adapt_elem_341_cur= hecMESH%n_adapt_elem_341_cur + 1
    hecMESH%n_adapt_elem_cur    = hecMESH%n_adapt_elem_cur     + 1

    icel = hecMESH%n_adapt_elem_cur
    icouN= icouN + inc

    if (icel.gt.hecMESH%ne_array) then
      call hecmw_adapt_ERROR_EXIT (hecMESH, 61)
    endif

    hecMESH%when_i_was_refined_elem(icel)= hecMESH%n_adapt
    hecMESH%elem_node_index(icel)= hecMESH%elem_node_index(icel-1) + 4

    is= hecMESH%elem_node_index(icel-1)
    hecMESH%elem_node_item(is+1)= in1
    hecMESH%elem_node_item(is+2)= in2
    hecMESH%elem_node_item(is+3)= in3
    hecMESH%elem_node_item(is+4)= in4

    hecMESH%adapt_parent(2*icel-1)= hecMESH%elem_ID(2*NPAR-1)
    hecMESH%adapt_parent(2*icel  )= hecMESH%elem_ID(2*NPAR  )

    hecMESH%elem_ID(2*icel-1)= icouN + hecMESH%ne_internal
    hecMESH%elem_ID(2*icel  )= hecMESH%elem_ID(2*NPAR  )

    hecMESH%elem_mat_ID_item(icel)= hecMESH%elem_mat_ID_item(NPAR)
    hecMESH%section_ID      (icel)= hecMESH%section_ID      (NPAR)

    hecMESH%adapt_type(icel)= 0

    if (NDIVSUM.eq.6) then
      hecMESH%adapt_level(icel)= hecMESH%adapt_level(NPAR) + 2
    else
      hecMESH%adapt_level(icel)= hecMESH%adapt_level(NPAR) + 1
    endif

    is= hecMESH%adapt_children_index(NPAR-1)
    hecMESH%adapt_children_item(2*(is+IDchi)-1)= icel
    hecMESH%adapt_children_item(2*(is+IDchi)-1)= icouN + hecMESH%ne_internal
    hecMESH%adapt_children_item(2*(is+IDchi)  )= hecMESH%my_rank

    hecMESH%adapt_children_local(is+IDchi)= icel

    is= hecMESH%adapt_children_index(icel-1)
    hecMESH%adapt_children_index(icel)= is + 8

    iS1= is + 1
    iS2= is + 2
    iS3= is + 3
    iS4= is + 4
    iS5= is + 5
    iS6= is + 6
    iS7= is + 7
    iS8= is + 8

    hecMESH%adapt_children_item(2*iS1)= -1
    hecMESH%adapt_children_item(2*iS2)= -1
    hecMESH%adapt_children_item(2*iS3)= -1
    hecMESH%adapt_children_item(2*iS4)= -1
    hecMESH%adapt_children_item(2*iS5)= -1
    hecMESH%adapt_children_item(2*iS6)= -1
    hecMESH%adapt_children_item(2*iS7)= -1
    hecMESH%adapt_children_item(2*iS8)= -1

    hecMESH%adapt_children_item(2*iS1-1)=  0
    hecMESH%adapt_children_item(2*iS2-1)=  0
    hecMESH%adapt_children_item(2*iS3-1)=  0
    hecMESH%adapt_children_item(2*iS4-1)=  0
    hecMESH%adapt_children_item(2*iS5-1)=  0
    hecMESH%adapt_children_item(2*iS6-1)=  0
    hecMESH%adapt_children_item(2*iS7-1)=  0
    hecMESH%adapt_children_item(2*iS8-1)=  0

    hecMESH%elem_type        (icel)= PAR_CEL_TYP
    hecMESH%adapt_parent_type(icel)= hecMESH%adapt_type(NPAR)

  end subroutine hecmw_adapt_CREATE_NEW_TETRA
end subroutine hecmw_adapt_NEW_CELL_341


text:
      SUBROUTINE MESGBF(LUNIT,MESGTYP)

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    MESGBF
C   PRGMMR: WOOLLEN          ORG: NP20       DATE: 1994-01-06
C
C ABSTRACT: THIS SUBROUTINE READS AND EXAMINES SECTION 1 OF MESSAGES
C  IN A BUFR FILE IN SEQUENCE UNTIL IT FINDS THE FIRST MESSAGE THAT
C  IS NOT A BUFR TABLE (DICTIONARY) (I.E., NOT MESSAGE TYPE 11).  IT
C  THEN RETURNS THE MESSAGE TYPE FOR THIS FIRST NON-DICTIONARY MESSAGE.
C  THE BUFR FILE SHOULD NOT BE OPEN VIA BUFR ARCHIVE LIBRARY SUBROUTINE
C  OPENBF PRIOR TO CALLING THIS SUBROUTINE HOWEVER THE BUFR FILE MUST
C  BE CONNECTED TO UNIT LUNIT.  THIS SUBROUTINE IS IDENTICAL TO BUFR
C  ARCHIVE LIBRARY SUBROUTINE MESGBC EXCEPT MESGBC RETURNS THE MESSAGE
C  TYPE FOR THE FIRST NON-DICTIONARY MESSAGE THAT ACTUALLY CONTAINS
C  REPORT DATA (WHEREAS MESGBF WOULD RETURN THE REPORT TYPE OF A DUMMY
C  MESSAGE CONTAINING THE CENTER TIME FOR DUMP FILES), AND MESGBC ALSO
C  INDICATES WHETHER OR NOT THE FIRST REPORT DATA MESSAGE CONTAINS
C  REPORT DATA IS BUFR COMPRESSED.  (MESGBC ALSO HAS AN OPTION TO
C  OPERATE ON THE CURRENT BUFR STORED IN MEMORY, SOMETHING MESGBF
C  CANNOT DO.)
C
C PROGRAM HISTORY LOG:
C 1994-01-06  J. WOOLLEN -- ORIGINAL AUTHOR
C 2000-09-19  J. WOOLLEN -- MAXIMUM MESSAGE LENGTH INCREASED FROM
C                           10,000 TO 20,000 BYTES
C 2003-11-04  S. BENDER  -- ADDED REMARKS/BUFRLIB ROUTINE
C                           INTERDEPENDENCIES
C 2003-11-04  D. KEYSER  -- UNIFIED/PORTABLE FOR WRF; ADDED
C                           DOCUMENTATION (INCLUDING HISTORY)
C 2004-08-09  J. ATOR    -- MAXIMUM MESSAGE LENGTH INCREASED FROM
C                           20,000 TO 50,000 BYTES
C 2005-11-29  J. ATOR    -- USE IUPBS01 AND RDMSGW
C
C USAGE:    CALL MESGBF (LUNIT, MESGTYP)
C   INPUT ARGUMENT LIST:
C     LUNIT    - INTEGER: FORTRAN LOGICAL UNIT NUMBER FOR BUFR FILE
C
C   OUTPUT ARGUMENT LIST:
C     MESGTYP  - INTEGER: BUFR MESSAGE TYPE FOR FIRST NON-DICTIONARY
C                MESSAGE
C                      -1 = no messages read or error
C                      11 = if only BUFR table messages in BUFR file
C
C   INPUT FILES:
C     UNIT "LUNIT" - BUFR FILE
C
C REMARKS:
C    THIS ROUTINE CALLS:        IUPBS01  RDMSGW   WRDLEN
C    THIS ROUTINE IS CALLED BY: None
C                               Normally called only by application
C                               programs.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

      INCLUDE 'bufrlib.prm'

      DIMENSION    MBAY(MXMSGLD4)

C-----------------------------------------------------------------------
C-----------------------------------------------------------------------

      MESGTYP = -1

C  SINCE OPENBF HAS NOT YET BEEN CALLED, MUST CALL WRDLEN TO GET
C  MACHINE INFO NEEDED LATER
C  -------------------------------------------------------------

      CALL WRDLEN

C  READ PAST ANY BUFR TABLES AND RETURN THE FIRST MESSAGE TYPE FOUND
C  -----------------------------------------------------------------

      REWIND LUNIT

1     CALL RDMSGW(LUNIT,MBAY,IER)
      IF(IER.LT.0) GOTO 100

      MESGTYP = IUPBS01(MBAY,'MTYP')
      IF(MESGTYP.EQ.11) GOTO 1

      REWIND LUNIT

C  EXIT
C  ----

100   RETURN
      END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak