text:
! This file is part of DASH.
! SPDX-Identifier: MIT
!
! Copyright 2006 Cambridge Crystallographic Data Centre
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!
!*****************************************************************************
!
      SUBROUTINE Launch_TOPAS(input_file_name)

      USE dash_gui_resources
      USE WINTERACTER
      USE TAVAR

      IMPLICIT NONE

      CHARACTER*(*), INTENT (IN   ) :: input_file_name

      LOGICAL, EXTERNAL :: DASHWDialogGetCheckBoxLogical
      CHARACTER*20, EXTERNAL :: Integer2String
      LOGICAL exists, run_TOPAS_in_background
      INTEGER M, I, tLen
      CHARACTER*20 stage_str
      CHARACTER*(255) tDirName, tFileName
      INTEGER hFile

      ! Launch TOPAS and wait for it to return
      CALL PushActiveWindowID
      CALL SelectDASHDialog(IDD_Configuration)
      run_TOPAS_in_background = DASHWDialogGetCheckBoxLogical(IDC_TOPAS_in_background)
!      CALL DASHWDialogGetString(IDF_TOPASExe, TOPASEXE)
      CALL PopActiveWindowID
      tLen = LEN_TRIM(input_file_name)
      I = LEN_TRIM(TOPASEXE)
      CALL SplitPath(TOPASEXE(1:I), tDirName, tFileName)
      IF ( run_TOPAS_in_background ) THEN
        IF ( I .EQ. 0 ) THEN
          CALL ErrorMessage("DASH could not launch TOPAS. No executable is currently specified."//CHAR(13)//&
                            "This can be changed in the Configuration... window"//CHAR(13)//&
                            "under Options in the menu bar.")
          CALL WCursorShape(CurCrossHair)
          RETURN
        ENDIF
        INQUIRE(FILE = TOPASEXE(1:I),EXIST=exists)
        IF (.NOT. exists) GOTO 998
        M = InfoError(1) ! Clear errors
        CALL IOsDirChange(tDirName)
!        CALL IOSCommand('CMD.EXE /C '//TOPASEXE(1:I)//' '//input_file_name(1:tLen), ProcBlocked)
        CALL IOSCommand(TOPASEXE(1:I)//' '//'"'//input_file_name(1:tLen)//'"', ProcBlocked)
        IF ( InfoError(1) .NE. 0 ) GOTO 998
        CALL WCursorShape(CurCrossHair)
      ELSE
        ! Write out the launch_file.txt file
        IF ( ext_RR_stage .EQ. 1 ) THEN
          IF ( I .EQ. 0 ) THEN
            CALL InfoMessage("The launch_file.txt could not be written,"//CHAR(13)//&
                             "because no TOPAS directory has been entered."//CHAR(13)//&
                             "This can be changed in the Configuration... window"//CHAR(13)//&
                             "under Options in the menu bar.")
          ELSE
            hFile = 112
            OPEN(UNIT=hFile, FILE=TRIM(tDirName)//'launch_file.txt', STATUS='unknown', ERR=999)
            WRITE(hFile, '(A)', ERR=999) '"'//input_file_name(1:tLen)//'.inp"'
            GOTO 100
  999       CONTINUE ! Error writing out launch_file.txt: we don't really care.
!           CALL InfoMessage("launch_file.txt could not be written.")
            CALL DebugErrorMessage("launch_file.txt could not be written.")
  100       CLOSE(hFile)
          ENDIF
        ENDIF
!        CALL InfoMessage('TOPAS .inp file for Pawley has been written.')
        stage_str = Integer2String(ext_RR_stage)
        CALL InfoMessage('TOPAS .inp file '//TRIM(stage_str)//' has been written.')
      ENDIF
      RETURN
  998 CALL ErrorMessage("DASH could not launch TOPAS. The executable is currently configured"//CHAR(13)//&
                        "to launch the program "//TOPASEXE(1:I)//CHAR(13)//&
                        "This can be changed in the Configuration... window"//CHAR(13)//&
                        "under Options in the menu bar.")
      CALL WCursorShape(CurCrossHair)
      RETURN

      END SUBROUTINE Launch_TOPAS
!
!*****************************************************************************
!
      INTEGER FUNCTION WriteTOPASFilePawley(TheFileName)

      USE dash_gui_resources
      USE WINTERACTER
      USE VARIABLES
      USE TAVAR

      IMPLICIT NONE

      CHARACTER*(*), INTENT (IN   ) :: TheFileName

      INCLUDE 'params.inc'
      INCLUDE 'Lattice.inc'
      INCLUDE 'GLBVAR.INC'

      INTEGER          NBIN, LBIN
      REAL                         XBIN,       YOBIN,       YCBIN,       YBBIN,       EBIN,       AVGESD
      COMMON /PROFBIN/ NBIN, LBIN, XBIN(MOBS), YOBIN(MOBS), YCBIN(MOBS), YBBIN(MOBS), EBIN(MOBS), AVGESD

      LOGICAL, EXTERNAL :: DASHWDialogGetCheckBoxLogical
      REAL, EXTERNAL :: FnWavelengthOfMenuOption
      INTEGER hFileTOPAS, hTempFile
      INTEGER iLen, iLen_1, i, tLen
      CHARACTER*(1) tChar
      CHARACTER*(MaxPathLength) tDirName, tFileName
      CHARACTER*(40) space_group_str
      INTEGER tStrLen
      INTEGER tIRadSelection
      CHARACTER*(6) bkg(4)
      DATA bkg/' 100.0', '-100.0', ' 10.0', ' 10.0'/

      ! The way this code has curently been written, this routine can only be called
      ! from one of the Wizard windows as part of a "iRietveldMethod" Rietveld refinement
      ! Initialise to failure
      WriteTOPASFilePawley = 1
      hFileTOPAS = 116
      tLen = LEN_TRIM(TheFileName)
      OPEN(UNIT=hFileTOPAS, FILE=TheFileName(1:tLen), STATUS='unknown', ERR=999)
      WRITE(hFileTOPAS, '(A)', ERR=999) '#include "'// &
                                        TRIM(ShareDashDirectory)//'TOPAS.inc"'
      ! Ideally, we need to figure out if penalties_weighting_K1 can be made local to the "str" keyword.
      WRITE(hFileTOPAS, '(A)', ERR=999) 'penalties_weighting_K1 5'

      ! Ideally this next line would account for topas4/5 distinction
      ! but it is commented out so should not hurt either way
      WRITE(hFileTOPAS, '(A)', ERR=999) "'do_errors"
      WRITE(hFileTOPAS, '(A)', ERR=999) 'r_exp 1.0'
      WRITE(hFileTOPAS, '(A)', ERR=999) 'r_exp_dash 1.0'
      WRITE(hFileTOPAS, '(A)', ERR=999) 'r_wp 1.0'
      WRITE(hFileTOPAS, '(A)', ERR=999) 'r_wp_dash 1.0'
      WRITE(hFileTOPAS, '(A)', ERR=999) 'r_p 1.0'
      WRITE(hFileTOPAS, '(A)', ERR=999) 'r_p_dash 1.0'
      WRITE(hFileTOPAS, '(A)', ERR=999) 'weighted_Durbin_Watson 1.0'
      WRITE(hFileTOPAS, '(A)', ERR=999) 'gof 1.0'
      ! Next line contains the file name of the file containing the diffraction data
      ! TOPAS can read .raw files, but not if they contain multiple ranges, and
      ! TOPAS can read .xye files, but not if they contain a wavelength on the first line
      ! We need to write out the data to a temporary file
      ! This also has the advantage that the user can choose a different file
      ! from the one they used for structure solution.
      hTempFile = 117
      CALL SplitPath(TheFileName, tDirName, tFileName)
      iLen = LEN_TRIM(tDirName)
      tFileName = tDirName(1:iLen)//'temp.xye'
      iLen_1 = iLen + 8
      OPEN(UNIT=hTempFile,FILE=tFileName(1:iLen_1),STATUS='unknown',ERR=998)
      DO i = 1, NBIN
        WRITE(hTempFile, '(F10.5,1X,F12.4,1X,F12.5)', ERR=998) XBIN(i), YOBIN(i), EBIN(i)
      ENDDO
      CLOSE(hTempFile)
      WRITE(hFileTOPAS, '(A)', ERR=999) 'xdd "'//tFileName(1:iLen_1)//'" xye_format'
      WRITE(hFileTOPAS, '(A)', ERR=999) '  bkg @ '
      DO i = 1, NumOfBkgTerm
        IF ( i .LE. 4 ) THEN
          WRITE(hFileTOPAS, '(3X,A)', ERR=999) bkg(i)
        ELSE
          WRITE(hFileTOPAS, '(3X,A)', ERR=999) ' 0.0'
        ENDIF
      ENDDO
      WRITE(hFileTOPAS, '(A,F12.4)', ERR=999) '  start_X ', XBIN(1)
      WRITE(hFileTOPAS, '(A,F12.4)', ERR=999) '  finish_X ', XBIN(NBIN)
      ! We do not write out the zero-point error in case the user has decided to use a different data set
      WRITE(hFileTOPAS, '(A)', ERR=999) '  Zero_Error(@, 0.0)'
      IF ( JRadOption .EQ. 1 ) THEN
        WRITE(hFileTOPAS, '(A)', ERR=999) '  LP_Factor( 26)'
      ELSE IF ( JRadOption .EQ. 2 ) THEN
        WRITE(hFileTOPAS, '(A)', ERR=999) '  LP_Factor( 90)'
      ELSE
        CALL ErrorMessage('Unknown radiation type.')
      ENDIF
      WRITE(hFileTOPAS, '(A)', ERR=999) '  Rp 217.5'
      WRITE(hFileTOPAS, '(A)', ERR=999) '  Rs 217.5'
!      WRITE(hFileTOPAS, '(A)', ERR=999) '  Simple_Axial_Model(@, 9.0)'
      WRITE(hFileTOPAS, '(A)', ERR=999) '  axial_conv'
      IF ( JRadOption .EQ. 1 ) THEN
        WRITE(hFileTOPAS, '(A)', ERR=999) '    filament_length   @   6.0'
        WRITE(hFileTOPAS, '(A)', ERR=999) '    sample_length   @   3.5'
        WRITE(hFileTOPAS, '(A)', ERR=999) '    receiving_slit_length   @  12.0'
      ELSE
        WRITE(hFileTOPAS, '(A)', ERR=999) '    filament_length   @   1.0'
        WRITE(hFileTOPAS, '(A)', ERR=999) '    sample_length   @   1.0'
        WRITE(hFileTOPAS, '(A)', ERR=999) '    receiving_slit_length   @  3.0'
      ENDIF
      WRITE(hFileTOPAS, '(A)', ERR=999) '    axial_n_beta  20'
      WRITE(hFileTOPAS, '(A)', ERR=999) '    axial_del  0.0053'
      WRITE(hFileTOPAS, '(A)', ERR=999) '  lam'
      WRITE(hFileTOPAS, '(A)', ERR=999) '    ymin_on_ymax 0.001'
      CALL SelectDASHDialog(IDD_PW_Page4)
      IF ( DASHWDialogGetCheckBoxLogical(IDC_Monochromated) .OR. (JRadOption .NE. 1) ) THEN
        WRITE(hFileTOPAS, '(A,F9.6)', ERR=999) '    la 1 lo ', ALambda
      ELSE
        tIRadSelection = -1
        DO I = 2, 6
          IF (ABS(ALambda - FnWavelengthOfMenuOption(I)) .LT. 0.0003) tIRadSelection = I
        ENDDO
        SELECT CASE (tIRadSelection)
        CASE (2) ! Cu
          WRITE(hFileTOPAS, '(A)', ERR=999) '    la 0.6666667 lo 1.54060'
          WRITE(hFileTOPAS, '(A)', ERR=999) '    la 0.3333333 lo 1.54439'
        CASE (3) ! Mo
          WRITE(hFileTOPAS, '(A)', ERR=999) '    la 0.6666667 lo 0.70930'
          WRITE(hFileTOPAS, '(A)', ERR=999) '    la 0.3333333 lo 0.71359'
        CASE (4) ! Co
          WRITE(hFileTOPAS, '(A)', ERR=999) '    la 0.6666667 lo 1.78897'
          WRITE(hFileTOPAS, '(A)', ERR=999) '    la 0.3333333 lo 1.79285'
        CASE (5) ! Cr
          WRITE(hFileTOPAS, '(A)', ERR=999) '    la 0.6666667 lo 2.28970'
          WRITE(hFileTOPAS, '(A)', ERR=999) '    la 0.3333333 lo 2.29361'
        CASE (6) ! Fe
          WRITE(hFileTOPAS, '(A)', ERR=999) '    la 0.6666667 lo 1.93604'
          WRITE(hFileTOPAS, '(A)', ERR=999) '    la 0.3333333 lo 1.93998'
        CASE DEFAULT
          CALL InfoMessage("Wavelength not recognised as a standard anode material--monochromated assumed.")
          WRITE(hFileTOPAS, '(A,F9.6)', ERR=999) '    la 1 lo ', ALambda
        END SELECT
      ENDIF
      WRITE(hFileTOPAS, '(A,F8.6)', ERR=999) '  x_calculation_step ', (XBIN(NBIN) - XBIN(1)) / (NBIN - 1)
      WRITE(hFileTOPAS, '(A)', ERR=999) '  hkl_Is'
      WRITE(hFileTOPAS, '(A)', ERR=999) '    CS_L(@, 5000.00)' !200.00
      WRITE(hFileTOPAS, '(A)', ERR=999) '    CS_G(@, 200.00)'
      IF ( JRadOption .EQ. 1 ) THEN
        WRITE(hFileTOPAS, '(A)', ERR=999) '    Strain_G(@, 0.2)'
      ELSE
        WRITE(hFileTOPAS, '(A)', ERR=999) '    Strain_G(@, 0.01)'
      ENDIF
      WRITE(hFileTOPAS, '(A)', ERR=999) '    Strain_L(@, 0.1)' !0.01
      IF ( LatBrav .LT. 6 ) THEN ! Triclinic through orthorhombic: no constraints on a, b or c.
        WRITE(hFileTOPAS, '(A,F10.5)', ERR=999) '    a  @ ', CELLPAR(1)
        WRITE(hFileTOPAS, '(A,F10.5)', ERR=999) '    b  @ ', CELLPAR(2)
        WRITE(hFileTOPAS, '(A,F10.5)', ERR=999) '    c  @ ', CELLPAR(3)
      ELSE
        ! Unit-cell parameters b or c or both are constrained:
        ! we must define a in terms of a TOPAS "parameter".
        WRITE(hFileTOPAS, '(A,F10.5)', ERR=999) '    a  uc_prm ', CELLPAR(1)
        WRITE(hFileTOPAS, '(A)', ERR=999) '    b = uc_prm;'
        IF ( CellParConstrained(3) ) THEN
          WRITE(hFileTOPAS, '(A)', ERR=999) '    c = uc_prm;'
        ELSE
          WRITE(hFileTOPAS, '(A,F10.5)', ERR=999) '    c  @ ', CELLPAR(3)
        ENDIF
      ENDIF
      IF ( LatBrav .EQ. 8 ) THEN ! Rhombohedral
        ! Unit-cell parameters beta and gamma are constrained:
        ! we must define alpha in terms of a TOPAS "parameter".
        WRITE(hFileTOPAS, '(A,F10.5)', ERR=999) '    al  uc_ang_prm ', CELLPAR(4)
        WRITE(hFileTOPAS, '(A)', ERR=999) '    be = uc_ang_prm;'
        WRITE(hFileTOPAS, '(A)', ERR=999) '    ga = uc_ang_prm;'
      ELSE
        IF ( CellParConstrained(4) ) THEN
          tChar = ' '
        ELSE
          tChar = '@'
        ENDIF
        WRITE(hFileTOPAS, '(A,F10.5)', ERR=999) '    al '//tChar//' ', CELLPAR(4)
        IF ( CellParConstrained(5) ) THEN
          tChar = ' '
        ELSE
          tChar = '@'
        ENDIF
        WRITE(hFileTOPAS, '(A,F10.5)', ERR=999) '    be '//tChar//' ', CELLPAR(5)
        IF ( CellParConstrained(6) ) THEN
          tChar = ' '
        ELSE
          tChar = '@'
        ENDIF
        WRITE(hFileTOPAS, '(A,F10.5)', ERR=999) '    ga '//tChar//' ', CELLPAR(6)
      ENDIF
!      WRITE(hFileTOPAS, '(A)', ERR=999) '    phase_name dummy'
      ! By writing the space group name last, the extra information on hkl's and intensities
      ! will be appended at the end of the file. That makes it a lot easier to read the .out
      ! file back in again and to just discard everything following after the space group.

      ! TODO ##### for space groups higher than orthorhombic we're almost certainly better off
      ! using the format space_group 222:1 rather than the space-group name.

      space_group_str = TRIM(SGHMaStr(NumberSGTable))
      CALL StrRemoveSpaces(space_group_str, tStrLen)
      WRITE(hFileTOPAS, '(A)', ERR=999) '    space_group '//space_group_str(1:tStrLen)
!      WRITE(hFileTOPAS, '(A,A,A)', ERR=999) '    Out_X_Ycalc("', tDirName(1:iLen), 'temp.txt")'
      WriteTOPASFilePawley = 0
      CLOSE(hFileTOPAS)

!      CALL InfoMessage(TheFileName(1:tLen-4))

      CALL Launch_TOPAS(TheFileName(1:tLen-4))
      RETURN
  999 CALL ErrorMessage("Error writing TOPAS input file (Pawley).")
      CLOSE(hFileTOPAS)
  998 CALL ErrorMessage("Error writing temporary diffraction data file.")
      CLOSE(hTempFile)
      CLOSE(hFileTOPAS)

      END FUNCTION WriteTOPASFilePawley
!
!*****************************************************************************
!
      SUBROUTINE UpdateTOPASCheckBoxes

      USE dash_gui_resources
      USE VARIABLES
      USE TAVAR

      IMPLICIT NONE

      LOGICAL, EXTERNAL :: DASHWDialogGetCheckBoxLogical

      CALL PushActiveWindowID
      CALL SelectDASHDialog(IDD_SAW_Page7_TOPAS)
      IF ( ext_RR_stage .GE. 3 ) THEN
        CALL WDialogFieldState(IDC_Use_anisotropic_broadening,         Disabled)
        CALL WDialogFieldState(IDB_RESTORE,                            Enabled)
      ELSE
        ! Initial state
        CALL WDialogFieldState(IDC_Use_anisotropic_broadening,         Enabled)
        CALL WDialogFieldState(IDC_Anisotropic_broadening,             Disabled)
        CALL WDialogPutCheckBoxLogical(IDC_Use_anisotropic_broadening, .TRUE.)
        CALL WDialogPutCheckBoxLogical(IDC_Anisotropic_broadening,     .TRUE.)
        CALL WDialogFieldState(IDB_RESTORE,                            Disabled)
      ENDIF
      IF ( .NOT. DASHWDialogGetCheckBoxLogical(IDC_UseDASHRecommendation) ) THEN
        IF ( ext_RR_stage .GE. 3 .AND. DASHWDialogGetCheckBoxLogical(IDC_Use_anisotropic_broadening) ) &
            CALL WDialogFieldState(IDC_Anisotropic_broadening,         Enabled)
        IF ( ext_RR_stage .EQ. 3 ) &
          CALL WDialogPutCheckBoxLogical(IDC_Anisotropic_broadening,   .FALSE.)
        RETURN
      ENDIF
      SELECT CASE ( ext_RR_stage )
        CASE ( 2 ) ! Anisotropic Pawley refinement
          CALL WDialogPutCheckBoxLogical(IDC_Anisotropic_broadening,   .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_Scale,       .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_Background,  .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_Biso,        .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_Coordinates, .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_IncludeESDs, .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_SaveProfile, .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_WriteCIF,    .TRUE.)
        CASE ( 3 ) ! First Rietveld refinement: refine scale only
          CALL WDialogPutCheckBoxLogical(IDC_Anisotropic_broadening,   .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_Scale,       .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_Background,  .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_Biso,        .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_Coordinates, .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_IncludeESDs, .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_SaveProfile, .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_WriteCIF,    .TRUE.)
        CASE ( 4 ) ! Second Rietveld refinement: include background
          CALL WDialogPutCheckBoxLogical(IDC_Anisotropic_broadening,   .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_Scale,       .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_Background,  .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_Biso,        .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_Coordinates, .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_IncludeESDs, .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_SaveProfile, .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_WriteCIF,    .TRUE.)
        CASE ( 5 ) ! Third Rietveld refinement: include Biso
          CALL WDialogPutCheckBoxLogical(IDC_Anisotropic_broadening,  .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_Scale,       .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_Background,  .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_Biso,        .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_Coordinates, .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_IncludeESDs, .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_SaveProfile, .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_WriteCIF,    .TRUE.)
        CASE ( 6 ) ! Fourth Rietveld refinement: include coordinates
          CALL WDialogPutCheckBoxLogical(IDC_Anisotropic_broadening,  .FALSE.)
          CALL WDialogPutCheckBoxLogical(IDC_Scale,       .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_Background,  .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_Biso,        .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_Coordinates, .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_IncludeESDs, .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_SaveProfile, .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_WriteCIF,    .TRUE.)
        CASE ( 7 ) ! Fifth Rietveld refinement: reenable+include aniso broadening
          IF ( DASHWDialogGetCheckBoxLogical(IDC_Use_anisotropic_broadening) ) THEN
            CALL WDialogFieldState(IDC_Anisotropic_broadening,         Enabled)
            CALL WDialogPutCheckBoxLogical(IDC_Anisotropic_broadening, .TRUE.)
          ENDIF
          CALL WDialogPutCheckBoxLogical(IDC_Scale,       .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_Background,  .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_Biso,        .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_Coordinates, .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_IncludeESDs, .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_SaveProfile, .TRUE.)
          CALL WDialogPutCheckBoxLogical(IDC_WriteCIF,    .TRUE.)
      END SELECT
      CALL PopActiveWindowID

      END SUBROUTINE UpdateTOPASCheckBoxes
!
!*****************************************************************************
!
! The initialisation of this window is done in DealWithWizardWindowProfileRange()
! in Wizard_routines.f90
      SUBROUTINE DealWithRR_TOPAS

      USE dash_gui_resources
      USE VARIABLES
      USE TAVAR

      IMPLICIT NONE

      LOGICAL, EXTERNAL :: DASHWDialogGetCheckBoxLogical, Confirm
      INTEGER, EXTERNAL :: WriteTOPASFileRietveld2
      INTEGER, EXTERNAL :: WriteTOPASPawleyAnisotropic
      INTEGER, EXTERNAL :: BackupRestorTopasFiles

      INTEGER iStat, iRetCode

      CALL PushActiveWindowID
      CALL SelectDASHDialog(IDD_SAW_Page7_TOPAS)
      SELECT CASE (EventType)
        CASE (PushButton) ! one of the buttons was pushed
          SELECT CASE (EventInfo%VALUE1)
            CASE (IDBACK)
              CALL WizardWindowShow(IDD_RR_External)
            CASE (IDCANCEL, IDCLOSE)
              IF ( iRietveldMethod .NE. INTERNAL_RB ) THEN
                CALL CopyBackup2Pattern()
                iRietveldMethod = INTERNAL_RB
              ENDIF
              CALL EndWizardPastPawley
            CASE (IDB_RESTORE)
              IF (Confirm('This will restore .inp and .out files to the moment'//CHAR(13)// &
                          'immediately after last writting.'//CHAR(13)//CHAR(13)// &
                          'Continue to verwrite current .inp and .out files ?')) THEN
                IF (BackupRestorTopasFiles(ext_RR_input_file_name, .FALSE.) .EQ. 0) THEN
                  CALL InfoMessage('The backuped .inp and .out files have been restored.'//CHAR(13))
                ELSE
                  CALL ErrorMessage("Error while trying to restore .inp and/or .out files."//CHAR(13))
                ENDIF
              ENDIF
            CASE (IDB_WRITE)
              ! ## The way I have programmed it at the moment (inserting an extra stage "2" that
              ! is simply skipped when no aniso) means that we could move the two lines:
              !      ext_RR_stage = ext_RR_stage + 1
              !      CALL UpdateTOPASCheckBoxes()
              ! to Launch_TOPAS()

              IF ( ext_RR_stage .EQ. 2 ) THEN
                ! Anisotropic broadening...
                use_anisotropic_broadening = DASHWDialogGetCheckBoxLogical(IDC_Use_anisotropic_broadening)
                IF ( use_anisotropic_broadening ) THEN
                  IF ( WriteTOPASPawleyAnisotropic(ext_RR_input_file_name) .EQ. 0 ) THEN
                    CALL Launch_TOPAS(ext_RR_input_file_name)
                    ext_RR_stage = ext_RR_stage + 1
                    iRetCode = BackupRestorTopasFiles(ext_RR_input_file_name, .TRUE.)
                    CALL UpdateTOPASCheckBoxes()
                  ENDIF
                  CALL PopActiveWindowID
                  RETURN
                ENDIF
                ext_RR_stage = ext_RR_stage + 1
                CALL UpdateTOPASCheckBoxes()
              ENDIF
              IF ( WriteTOPASFileRietveld2(ext_RR_input_file_name) .EQ. 0 ) THEN
                CALL Launch_TOPAS(ext_RR_input_file_name)
                ext_RR_stage = ext_RR_stage + 1
                iRetCode = BackupRestorTopasFiles(ext_RR_input_file_name, .TRUE.)
                CALL UpdateTOPASCheckBoxes()
              ENDIF
            CASE (IDB_View)
              CALL Launch_Viewer(ext_RR_input_file_name, '.cif')
          END SELECT
        CASE (FieldChanged)
          SELECT CASE (EventInfo%VALUE1)
            CASE (IDC_Use_anisotropic_broadening)
              CALL WDialogPutCheckBoxLogical(IDC_Anisotropic_broadening, &
                       DASHWDialogGetCheckBoxLogical(IDC_Use_anisotropic_broadening))
            CASE (IDC_UseDASHRecommendation)
              iStat = Disabled
              IF (DASHWDialogGetCheckBoxLogical(IDC_UseDASHRecommendation)) THEN
                IF (ext_RR_stage .GE. 7 .AND. DASHWDialogGetCheckBoxLogical(IDC_Use_anisotropic_broadening)) &
                  iStat = Enabled
              ELSE
                IF (ext_RR_stage .GE. 3 .AND. DASHWDialogGetCheckBoxLogical(IDC_Use_anisotropic_broadening)) &
                  iStat = Enabled
              ENDIF
              CALL WDialogFieldState(IDC_Anisotropic_broadening, iStat)
          END SELECT
      END SELECT
      CALL PopActiveWindowID

      END SUBROUTINE DealWithRR_TOPAS
!
!*****************************************************************************
!
! This function reads the Pawley refinement file and and writes the file for
! Pawley refinement with anisotropic broadening.
      INTEGER FUNCTION WriteTOPASPawleyAnisotropic(FileNameBase)

      ! If we have just done the Pawley refinement with isotropic broadening, then the
      ! next step could be Pawley refinement with anisotropic broadening.
      ! To be written out after "Strain_L". All other parameters remain variable.

      USE dash_gui_resources
      USE VARIABLES
      USE ATMVAR
      USE ZMVAR
      USE PO_VAR
      USE RRVAR
      USE SAMVAR

      IMPLICIT NONE

      CHARACTER*(*), INTENT (INOUT) :: FileNameBase

      LOGICAL, EXTERNAL :: DASHWDialogGetCheckBoxLogical
      INTEGER hFileTOPAS, hOutputFile
      INTEGER iLen
      LOGICAL is_last_line
      CHARACTER(512) tLine
      CHARACTER(30) word
      INTEGER       word_len
      CHARACTER(MaxPathLength) FileNameToRead, FileNameToWrite
      CHARACTER*255 tDirName, tFileName, tExtension
      INTEGER ExtLength

      ! Initialise to failure
      WriteTOPASPawleyAnisotropic = 1
      ExtLength = 3
      CALL SplitPath2(FileNameBase, tDirName, tFileName, tExtension, ExtLength)
      FileNameBase = TRIM(tDirName)//tFileName
      FileNameToRead = TRIM(FileNameBase)//'.out'
      FileNameToWrite = TRIM(FileNameBase)//'.inp'
      hFileTOPAS = 116
      hOutputFile = 117 ! This is the file that is being *read*
      OPEN(UNIT=hFileTOPAS, FILE=FileNameToWrite, STATUS='unknown', ERR=999)
      OPEN(UNIT=hOutputFile, FILE=FileNameToRead, STATUS='unknown', ERR=998)
      is_last_line = .FALSE.
   10 CONTINUE
      READ(hOutputFile, '(A)', ERR=998, END=40) tLine
      iLen = LEN_TRIM(tLine)
      word_len = 30
      CALL FirstWord(tLine, word, word_len)
      CALL StrUpperCase(word)
      WRITE(hFileTOPAS, '(A)', ERR=999) tLine(1:iLen)
      IF ( (word_len .EQ. 8) .AND. (word(1:8) .EQ. 'STRAIN_L') ) THEN
        WRITE(hFileTOPAS, '(A)', ERR=999) '    prm    sh_scale 0.01 min 0.0001'
        WRITE(hFileTOPAS, '(A)', ERR=999) '    spherical_harmonics_hkl sh'
        WRITE(hFileTOPAS, '(A)', ERR=999) '      sh_order 6'
        WRITE(hFileTOPAS, '(A)', ERR=999) '    lor_fwhm = sh_scale * sh;'
      ENDIF
      GOTO 10
   40 CLOSE(hOutputFile)

      WriteTOPASPawleyAnisotropic = 0
      CLOSE(hFileTOPAS)
      RETURN
  999 CALL ErrorMessage("Error writing TOPAS input file (Pawley anisotropic)")
      CLOSE(hFileTOPAS)
      RETURN
  998 CALL ErrorMessage("Error reading TOPAS output file (Pawley anisotropic)")
      CLOSE(hOutputFile)
      CLOSE(hFileTOPAS)

      END FUNCTION WriteTOPASPawleyAnisotropic
!
!*****************************************************************************
!
! This function reads and writes the file for Rietveld refinement (as opposed to Pawley).
      INTEGER FUNCTION WriteTOPASFileRietveld2(FileNameBase)

      USE dash_gui_resources
      USE VARIABLES
      USE PO_VAR
      USE RRVAR

      IMPLICIT NONE

      CHARACTER*(*), INTENT (INOUT) :: FileNameBase

      INCLUDE 'Lattice.inc'

      INTEGER, EXTERNAL :: StrFind
      CHARACTER*20, EXTERNAL :: Integer2String
      LOGICAL, EXTERNAL :: DASHWDialogGetCheckBoxLogical
      LOGICAL, EXTERNAL :: PutAtomsForSpecialPosition
      LOGICAL, EXTERNAL :: PutRestraints, TOPASWriteDistance, TOPASWriteAngle, TOPASWritePlane
      INTEGER hFileTOPAS, hOutputFile
      INTEGER iLen, iPos, iStrPos
      LOGICAL is_last_line
      CHARACTER(512) tLine
      CHARACTER(60) word
      INTEGER       word_len
      CHARACTER*20 tStr
      INTEGER K
      ! We need to remember which atoms have been assigned to an assembly
      CHARACTER(MaxPathLength) FileNameToRead, FileNameToWrite
      CHARACTER*255 tDirName, tFileName, tExtension
      INTEGER ExtLength, hSP_out_file
      LOGICAL was_Pawley
      LOGICAL in_spherical_harmonics
      INTEGER wscale_bond, wscale_angle, wscale_flatten

      was_Pawley = .FALSE. ! This is a remnant from when it was still possible
      ! to start from a previously generated .inp file.
      ! In principle, it looks at first glance as if we could now use ext_RR_stage .EQ. 3.
      ! However, it is much better to keep the existing mechanism based on interpreting the
      ! .out file: that way, if the user clicked "Write" twice in a row *without having run TOPAS*,
      ! either by accident or because they feel that scale and background can be refined together,
      ! the correct new .inp will be written out. If we used ext_RR_stage .EQ. 3,
      ! we would inadvertently miss the stage where the atomic coordinates are written out,
      ! which would be a disaster.

      CALL SelectDASHDialog(IDD_SAW_Page7_TOPAS)
      ! Initialise to failure
      WriteTOPASFileRietveld2 = 1
      ExtLength = 3
      CALL SplitPath2(FileNameBase, tDirName, tFileName, tExtension, ExtLength)
      FileNameBase = TRIM(tDirName)//tFileName
      FileNameToRead = TRIM(FileNameBase)//'.out'
      FileNameToWrite = TRIM(FileNameBase)//'.inp'
      hFileTOPAS = 116
      hOutputFile = 117 ! This is the file that is being *read*
      hSP_out_file = 60 ! Special positions output file
      OPEN(UNIT=hFileTOPAS, FILE=FileNameToWrite, STATUS='unknown', ERR=999)
      OPEN(UNIT=hOutputFile, FILE=FileNameToRead, STATUS='unknown', ERR=998)
      is_last_line = .FALSE.
      in_spherical_harmonics = .FALSE.
   10 CONTINUE
      READ(hOutputFile, '(A)', ERR=998, END=40) tLine
      iLen = LEN_TRIM(tLine)
      word_len = LEN(word)
      CALL FirstWord(tLine, word, word_len)
      CALL StrUpperCase(word)
      IF ( in_spherical_harmonics ) THEN
!!!!  prm sh_scale   0.00907`_0.00078 min .0001
!     spherical_harmonics_hkl sh
!       sh_order  6 load sh_Cij_prm {
!          y00   !sh_c00  1.00000
!          y20   sh_c20  -0.65363`_0.15875
!          y22p  sh_c22p -2.09821`_0.24794
!          y40   sh_c40  -1.90493`_0.23022
!          y42p  sh_c42p -0.59773`_0.09433
!          y44p  sh_c44p -1.38602`_0.17287
!          y60   sh_c60   0.11779`_0.02838
!          y62p  sh_c62p  0.03288`_0.18446
!          y64p  sh_c64p -0.25236`_0.14136
!          y66p  sh_c66p  0.19174`_0.15399
!          }
!!!!  lor_fwhm = sh_scale * sh;
        IF ( (word_len .EQ. 8) .AND. (word(1:8) .EQ. 'LOR_FWHM') ) THEN
          in_spherical_harmonics = .FALSE.
        ELSE
          ! Undocumented TOPAS stuff, such as y20 or k61
          IF ( word(1:1) .EQ. 'Y' .OR. word(1:1) .EQ. 'K' ) THEN
            IF ( DASHWDialogGetCheckBoxLogical(IDC_Anisotropic_broadening) ) THEN
              iPos = StrFind(tLine, iLen, ' !sh_', 5)
              IF ( iPos .NE. 0 .AND. tLine(iPos:iPos+7) .NE. ' !sh_c00' ) tLine(iPos+1:iPos+1) = ' '
            ELSE
              iPos = StrFind(tLine, iLen, '  sh_', 5)
              IF ( iPos .NE. 0 ) tLine(iPos+1:iPos+1) = '!'
            ENDIF
          ENDIF
        ENDIF
        WRITE(hFileTOPAS, '(A)', ERR=999) tLine(1:iLen)
      ELSE IF ( (word_len .EQ. 22) .AND. (word(1:22) .EQ. 'PENALTIES_WEIGHTING_K1') ) THEN
        CALL DASHWDialogGetInteger(IDF_K1, K)
        IF ( K .LT. 0 ) K = 1
        tStr = Integer2String(K)
        WRITE(hFileTOPAS, '(A)', ERR=999) 'penalties_weighting_K1 '//TRIM(tStr)
      ELSE IF ( (word_len .EQ. 6) .AND. (word(1:6) .EQ. 'HKL_IS') ) THEN
        was_Pawley = .TRUE.
        WRITE(hFileTOPAS, '(A)', ERR=999) '  str'
      ELSE IF ( (word_len .EQ. 11) .AND. (word(1:11) .EQ. 'SPACE_GROUP') ) THEN
        WRITE(hFileTOPAS, '(A)', ERR=999) tLine(1:iLen)
        ! This is where the very first Rietveld file is written out
        IF ( was_Pawley ) THEN
          ! #@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@
          !                             BEGIN
          ! #@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@
          IF ( DASHWDialogGetCheckBoxLogical(IDC_Scale) ) THEN
            WRITE(hFileTOPAS, '(A)', ERR=999) '    scale @ 0.001'
          ELSE
            WRITE(hFileTOPAS, '(A)', ERR=999) '    scale   0.001'
          ENDIF
! We must call ShowWizardWindowRietveld() here, which will fill
! Xato (and all the RR variables). The Wizard window is suppressed because of the iRietveldMethod flag.
          CALL ShowWizardWindowRietveld(RR_SA_Sol)
! Also need to write out PO if used during SA
          IF ( PrefParExists ) THEN
            WRITE(hFileTOPAS, '(A,F5.3,A,I3,1X,I3,1X,I3,1X,A)', ERR=999) '    PO(@, ', RR_PO, ', , ', PO_Direction(1:3), ')'
          ELSE
            WRITE(hFileTOPAS, '(A)', ERR=999) "'    PO(@, 1.0, , 0 0 1)"
          ENDIF
          IF ( DASHWDialogGetCheckBoxLogical(IDC_Coordinates) ) THEN
            WRITE(hFileTOPAS, '(A)', ERR=999) '    macro ref_flag { @ }'
          ELSE
            WRITE(hFileTOPAS, '(A)', ERR=999) '    macro ref_flag {   }'
          ENDIF
          ! #########################################################################################
          ! TOPAS needs unique atom labels ("site labels"), especially when applying restraints.
          ! This is a problem for DASH where we use user-supplied labels.
          ! This is especially a problem if there is more than one Z-matrix in the asymmetric unit.
          ! (Would it have been possible to solve this using the "scopes" in TOPAS?, I.e. define each
          ! Z-matrix in its own scope and then define its restraints within the same scope?)
          ! #########################################################################################

          IF ( PutAtomsForSpecialPosition(.TRUE.) ) GOTO 997
          OPEN(UNIT=hSP_out_file, FILE="special_positions.out", STATUS='unknown', ERR=996)
       45 CONTINUE
          READ(hSP_out_file, '(A)', ERR=996, END=50) tLine
          WRITE(hFileTOPAS, '(A)', ERR=999) TRIM(tLine)
          GOTO 45
       50 CONTINUE
          CLOSE(hSP_out_file)
          CALL DASHWDialogGetInteger(IDF_INTEGER_TOPAS_WSCALE_DIST, wscale_bond)
          CALL DASHWDialogGetInteger(IDF_INTEGER_TOPAS_WSCALE_ANGLE, wscale_angle)
          CALL DASHWDialogGetInteger(IDF_INTEGER_TOPAS_WSCALE_FLAT, wscale_flatten)
          WRITE(hFileTOPAS, '(A)', ERR=999) '    prm !bond_width 0'
          WRITE(hFileTOPAS, '(A,I10)', ERR=999) '    prm !bond_weight ', wscale_bond !10000
          WRITE(hFileTOPAS, '(A)', ERR=999) '    prm !angle_width 1'
          WRITE(hFileTOPAS, '(A,I10)', ERR=999) '    prm !angle_weight ', wscale_angle !1
          WRITE(hFileTOPAS, '(A)', ERR=999) '    prm !flatten_width 0'
          WRITE(hFileTOPAS, '(A,I10)', ERR=999) '    prm !flatten_weight ', wscale_flatten !1000000
          IF ( PutRestraints(hFileTOPAS, TOPASWriteDistance, TOPASWriteAngle, TOPASWritePlane) ) GOTO 999
          ! #@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@
          !                            END
          ! #@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@

          ! Probably need to improve mechanism to indicate end
          is_last_line = .TRUE.
        ENDIF
      ELSE IF ( (word_len .EQ. 3) .AND. (word(1:3) .EQ. 'BKG') ) THEN
        IF ( DASHWDialogGetCheckBoxLogical(IDC_Background) ) THEN
          WRITE(hFileTOPAS, '(A)', ERR=999) '  bkg @ '
        ELSE
          WRITE(hFileTOPAS, '(A)', ERR=999) '  bkg '
        ENDIF
      ELSE IF ( (word_len .EQ. 5) .AND. (word(1:5) .EQ. 'SCALE') ) THEN
        IF ( DASHWDialogGetCheckBoxLogical(IDC_Scale) ) THEN
          tLine(11:11) = '@'
        ELSE
          tLine(11:11) = ' '
        ENDIF
        WRITE(hFileTOPAS, '(A)', ERR=999) tLine(1:iLen)
      ELSE IF ( (word_len .EQ. 11) .AND. (word(1:11) .EQ. 'OUT_CIF_STR') ) THEN
        ! Do nothing: it will be added at the end if necessary
      ELSE IF ( (word_len .EQ. 16) .AND. (word(1:16) .EQ. 'OUT_CIF_STR_UISO') ) THEN
        ! Do nothing: it will be added at the end if necessary
      ELSE IF ( (word_len .EQ. 11) .AND. (word(1:11) .EQ. 'OUT_X_YCALC') ) THEN
        ! Do nothing: it will be added at the end if necessary
      ELSE IF ( (word_len .EQ. 10) .AND. (word(1:10) .EQ. 'OUT_X_YOBS') ) THEN
        ! Do nothing: it will be added at the end if necessary
      ELSE IF ( (word_len .EQ. 16) .AND. (word(1:16) .EQ. 'OUT_X_DIFFERENCE') ) THEN
        ! Do nothing: it will be added at the end if necessary
      ELSE IF ( (word_len .EQ. 29) .AND. (word(1:29) .EQ. 'OUT_YOBS_YCALC_AND_DIFFERENCE') ) THEN
        ! Do nothing: it will be added at the end if necessary
      ELSE IF ( (word_len .EQ. 11) .AND. (word(1:11) .EQ. 'OUT_PROFILE') ) THEN
        ! Do nothing: it will be added at the end if necessary
      ELSE IF ( (word_len .EQ. 8) .AND. (word(1:8) .EQ. 'OUT_TICK') ) THEN
        ! Do nothing: it will be added at the end if necessary
      ELSE IF ( (word_len .EQ. 19) .AND. (word(1:19) .EQ. 'OUT_POWDERDATABLOCK') ) THEN
        ! Do nothing: it will be added at the end if necessary
      ELSE IF ( ((word_len .EQ. 9) .AND. (word(1:9) .EQ. 'DO_ERRORS')) &
          .OR.  ((word_len .EQ. 10) .AND. (word(1:10) .EQ. "'DO_ERRORS")) &
          .OR.  ((word_len .EQ. 27) .AND. (word(1:27) .EQ. 'DO_ERRORS_INCLUDE_PENALTIES')) &
          .OR.  ((word_len .EQ. 28) .AND. (word(1:28) .EQ. "'DO_ERRORS_INCLUDE_PENALTIES")) ) THEN

        IF ( DASHWDialogGetCheckBoxLogical(IDC_IncludeESDs) ) THEN
          IF ( DASHWDialogGetCheckBoxLogical(IDC_topas5_format_check) ) THEN
             WRITE(hFileTOPAS, '(A)', ERR=999) 'do_errors_include_penalties'
          ELSE
             WRITE(hFileTOPAS, '(A)', ERR=999) 'do_errors'
          ENDIF
        ELSE
          IF ( DASHWDialogGetCheckBoxLogical(IDC_topas5_format_check) ) THEN
             WRITE(hFileTOPAS, '(A)', ERR=999) "'do_errors_include_penalties"
          ELSE
             WRITE(hFileTOPAS, '(A)', ERR=999) "'do_errors"
          ENDIF
        ENDIF

      ELSE IF ( (word_len .EQ. 3) .AND. (word(1:3) .EQ. 'PRM') ) THEN
        ! Need to check if this is the "prm sh_scale 1.0 min 0.0001" line.
        iPos = StrFind(tLine, iLen, 'sh_scale', 8)
        IF ( iPos .NE. 0 ) THEN
          IF ( DASHWDialogGetCheckBoxLogical(IDC_Anisotropic_broadening) ) THEN
            tLine(iPos-1:iPos-1) = ' '
          ELSE
            tLine(iPos-1:iPos-1) = '!'
          ENDIF
          in_spherical_harmonics = .TRUE.
        ELSE
          ! Need to check if this is the "prm bnonh 3.000" line.
          iPos = StrFind(tLine, iLen, 'bnonh', 5)
          IF ( iPos .NE. 0 ) THEN
            IF ( DASHWDialogGetCheckBoxLogical(IDC_Biso) ) THEN
              tLine(iPos-1:iPos-1) = ' '
            ELSE
              tLine(iPos-1:iPos-1) = '!'
            ENDIF
          ELSE
            ! Need to check if the line is for special_position eg "prm dash_2 0.4670".
            iPos = StrFind(tLine, iLen, 'dash_', 5)
            IF ( iPos .NE. 0 ) THEN
              IF ( DASHWDialogGetCheckBoxLogical(IDC_Coordinates) ) THEN
                tLine(iPos-1:iPos-1) = ' '
              ELSE
                tLine(iPos-1:iPos-1) = '!'
              ENDIF
              GOTO 20
            ENDIF
            IF ( INDEX(tLine, '!bond_weight') .GT. 0 ) THEN
              CALL DASHWDialogGetInteger(IDF_INTEGER_TOPAS_WSCALE_DIST, wscale_bond)
              WRITE(tLine, '(A,I10)', ERR=999) '    prm !bond_weight ', wscale_bond
            ELSE IF ( INDEX(tLine, '!angle_weight') .GT. 0 ) THEN
              CALL DASHWDialogGetInteger(IDF_INTEGER_TOPAS_WSCALE_ANGLE, wscale_angle)
              WRITE(tLine, '(A,I10)', ERR=999) '    prm !angle_weight ', wscale_angle
            ELSE IF ( INDEX(tLine, '!flatten_weight') .GT. 0 ) THEN
              CALL DASHWDialogGetInteger(IDF_INTEGER_TOPAS_WSCALE_FLAT, wscale_flatten)
              WRITE(tLine, '(A,I10)', ERR=999) '    prm !flatten_weight ', wscale_flatten
            ENDIF
          ENDIF
        ENDIF
 20     CONTINUE
        WRITE(hFileTOPAS, '(A)', ERR=999) TRIM(tLine)
      ELSE IF ( (word_len .EQ. 4) .AND. (word(1:4) .EQ. 'SITE') ) THEN
        ! We could read the fractional coordinates back in and use them to
        ! overlay the refined solution with the original one.
        WRITE(hFileTOPAS, '(A)', ERR=999) tLine(1:iLen)
      ELSE IF ( (word_len .EQ. 5) .AND. (word(1:5) .EQ. 'MACRO') ) THEN
        ! Check if rest is ref_flag
        IF ( StrFind(tLine, iLen, 'ref_flag', 8) .NE. 0 ) THEN
          IF ( DASHWDialogGetCheckBoxLogical(IDC_Coordinates) ) THEN
            WRITE(hFileTOPAS, '(A)', ERR=999) '    macro ref_flag { @ }'
          ELSE
            WRITE(hFileTOPAS, '(A)', ERR=999) '    macro ref_flag {   }'
          ENDIF
        ENDIF
      ELSE IF ( (word_len .EQ. 2) .AND. (word(1:2) .EQ. 'PO') ) THEN
        WRITE(hFileTOPAS, '(A)', ERR=999) tLine(1:iLen)
      ELSE IF ( (word_len .EQ. 1) .AND. (word(1:1) .EQ. 'A') ) THEN
        iPos = StrFind(tLine, iLen, 'uc_prm', 6)
        IF ( iPos .NE. 0 ) THEN
          IF ( DASHWDialogGetCheckBoxLogical(IDC_IncludeESDs) ) THEN
            tLine(iPos-1:iPos-1) = ' '
          ELSE
            tLine(iPos-1:iPos-1) = '!'
          ENDIF
        ELSE
          iPos = StrFind(tLine, iLen, 'a', 1)
          IF ( DASHWDialogGetCheckBoxLogical(IDC_IncludeESDs) ) THEN
            tLine(iPos+3:iPos+3) = '@'
          ELSE
            tLine(iPos+3:iPos+3) = ' '
          ENDIF
        ENDIF
        WRITE(hFileTOPAS, '(A)', ERR=999) tLine(1:iLen)
      ELSE IF ( (word_len .EQ. 1) .AND. (word(1:1) .EQ. 'B') ) THEN
        IF ( .NOT. CellParConstrained(2) ) THEN
          iPos = StrFind(tLine, iLen, 'b', 1)
          IF ( DASHWDialogGetCheckBoxLogical(IDC_IncludeESDs) ) THEN
            tLine(iPos+3:iPos+3) = '@'
          ELSE
            tLine(iPos+3:iPos+3) = ' '
          ENDIF
        ENDIF
        WRITE(hFileTOPAS, '(A)', ERR=999) tLine(1:iLen)
      ELSE IF ( (word_len .EQ. 1) .AND. (word(1:1) .EQ. 'C') ) THEN
        IF ( .NOT. CellParConstrained(3) ) THEN
          iPos = StrFind(tLine, iLen, 'c', 1)
          IF ( DASHWDialogGetCheckBoxLogical(IDC_IncludeESDs) ) THEN
            tLine(iPos+3:iPos+3) = '@'
          ELSE
            tLine(iPos+3:iPos+3) = ' '
          ENDIF
        ENDIF
        WRITE(hFileTOPAS, '(A)', ERR=999) tLine(1:iLen)
      ELSE IF ( (word_len .EQ. 2) .AND. (word(1:2) .EQ. 'AL') ) THEN
        iPos = StrFind(tLine, iLen, 'uc_ang_prm', 10)
        IF ( iPos .NE. 0 ) THEN
          IF ( DASHWDialogGetCheckBoxLogical(IDC_IncludeESDs) ) THEN
            tLine(iPos-1:iPos-1) = ' '
          ELSE
            tLine(iPos-1:iPos-1) = '!'
          ENDIF
        ELSE
          iPos = StrFind(tLine, iLen, 'al', 2)
          IF ( DASHWDialogGetCheckBoxLogical(IDC_IncludeESDs) .AND. (.NOT. CellParConstrained(4)) ) THEN
            tLine(iPos+3:iPos+3) = '@'
          ELSE
            tLine(iPos+3:iPos+3) = ' '
          ENDIF
        ENDIF
        WRITE(hFileTOPAS, '(A)', ERR=999) tLine(1:iLen)
      ELSE IF ( (word_len .EQ. 2) .AND. (word(1:2) .EQ. 'BE') ) THEN
        iPos = StrFind(tLine, iLen, 'uc_ang_prm', 10)
        IF ( iPos .EQ. 0 ) THEN
          iPos = StrFind(tLine, iLen, 'be', 2)
          IF ( DASHWDialogGetCheckBoxLogical(IDC_IncludeESDs) .AND. (.NOT. CellParConstrained(5)) ) THEN
            tLine(iPos+3:iPos+3) = '@'
          ELSE
            tLine(iPos+3:iPos+3) = ' '
          ENDIF
        ENDIF
        WRITE(hFileTOPAS, '(A)', ERR=999) tLine(1:iLen)
      ELSE IF ( (word_len .EQ. 2) .AND. (word(1:2) .EQ. 'GA') ) THEN
        iPos = StrFind(tLine, iLen, 'uc_ang_prm', 10)
        IF ( iPos .EQ. 0 ) THEN
          iPos = StrFind(tLine, iLen, 'ga', 2)
          IF ( DASHWDialogGetCheckBoxLogical(IDC_IncludeESDs) .AND. (.NOT. CellParConstrained(6)) ) THEN
            tLine(iPos+3:iPos+3) = '@'
          ELSE
            tLine(iPos+3:iPos+3) = ' '
          ENDIF
        ENDIF
        WRITE(hFileTOPAS, '(A)', ERR=999) tLine(1:iLen)
      ELSE IF ( word(1:1) .EQ. "'" ) THEN ! It's a comment
        WRITE(hFileTOPAS, '(A)', ERR=999) tLine(1:iLen)
      ELSE
        ! Change all the @'s to spaces.
        ! This switches off refinement of zero-point error, CS_L, CS_G, Strain_L, Strain_G
        DO iStrPos = 1, iLen
          IF (tLine(iStrPos:iStrPos) .EQ. '@') &
            tLine(iStrPos:iStrPos) = ' '
        ENDDO
        WRITE(hFileTOPAS, '(A)', ERR=999) tLine(1:iLen)
      ENDIF
      IF ( .NOT. is_last_line ) GOTO 10
   40 CLOSE(hOutputFile)
      IF ( DASHWDialogGetCheckBoxLogical(IDC_WriteCIF) ) THEN
!        WRITE(hFileTOPAS, '(A)', ERR=999) '    Out_CIF_STR("'//TRIM(FileNameBase)//'.cif")'
        WRITE(hFileTOPAS, '(A)', ERR=999) '    Out_CIF_STR_Uiso("'//TRIM(FileNameBase)//'.cif")'
      ENDIF
      IF ( DASHWDialogGetCheckBoxLogical(IDC_SaveProfile) ) THEN
        WRITE(hFileTOPAS, '(A)', ERR=999) '    Out_Profile("'//TRIM(FileNameBase)//'_plot.pro")'
        WRITE(hFileTOPAS, '(A)', ERR=999) '    Out_Tick("'//TRIM(FileNameBase)//'_plot.tic")'
        WRITE(hFileTOPAS, '(A)', ERR=999) '    Out_PowderDataBlock("'//TRIM(FileNameBase)//'_plot.pcif")'
      ENDIF
!      WRITE(hFileTOPAS, '(A)', ERR=999) '    Out_X_Ycalc("'//TRIM(FileNameBase)//'.pp")'
!      WRITE(hFileTOPAS, '(A)', ERR=999) '    Out_Yobs_Ycalc_and_Difference("'//TRIM(FileNameBase)//'_1.xco")'
      WriteTOPASFileRietveld2 = 0
      CLOSE(hFileTOPAS)
      RETURN
  999 CALL ErrorMessage("Error writing TOPAS input file (Rietveld).")
      CLOSE(hFileTOPAS)
      RETURN
  998 CALL ErrorMessage("Error reading TOPAS output file (Rietveld).")
      CLOSE(hOutputFile)
      CLOSE(hFileTOPAS)
  997 CALL ErrorMessage("Error writing special positions input file (Rietveld).")
      CLOSE(hFileTOPAS)
  996 CALL ErrorMessage("Error reading special positions output file (Rietveld).")
      CLOSE(hSP_out_file)
      CLOSE(hFileTOPAS)

      END FUNCTION WriteTOPASFileRietveld2
!
!*****************************************************************************
!
!      LOGICAL FUNCTION CheckTOPASFileName (name)
!
!      USE TAVAR
!
!      IMPLICIT NONE
!
!      CHARACTER*(*), INTENT (INOUT) :: name
!
!      CHARACTER(MaxPathLength) :: DirName
!      CHARACTER(MaxPathLength) :: FileName
!      CHARACTER(3) :: Extension
!      INTEGER iLen, dLen, iPos
!      LOGICAL something_changed

      ! It turns out that TOPAS cannot cope with file names that have dots in them,
      ! e.g. the following gives an error message:
      ! "Daresbury9.1_Nov2006.inp"
      ! The error message refers to "Daresbury9", indicating that everything
      ! after the first "." was assumed to be the file extension.
!      CheckTOPASFileName = .FALSE.
!10    CONTINUE
!      ExtLength = 3
!      CALL SplitPath2(name, DirName, FileName, Extension, ExtLength)
!      iLen = LEN_TRIM(FileName)
!      iPos = StrFind(FileName, iLen, '.', 1)
!      IF ( iPos .NE. 0 ) THEN
!        FileName(iPos:iPos) = "_"
!        name = TRIM(DirName)//FileName(1:iLen)//'.inp'
!        CheckTOPASFileName = .TRUE.
!        ! There could be more "." in the file name
!        GOTO 10
!      ENDIF
!
!      END FUNCTION CheckTOPASFileName
!
!*****************************************************************************
!
      LOGICAL FUNCTION TOPASWriteDistance (hFile, iAtom1, iAtom2, iElement1, iElement2, Distance)

      IMPLICIT NONE

      INTEGER, INTENT (IN) :: hFile, iAtom1, iAtom2, iElement1, iElement2
      REAL, INTENT (IN) :: Distance

      CHARACTER(7) label1, label2

      TOPASWriteDistance = .TRUE.
      CALL GenerateAtomLabel(iElement1, iAtom1, label1)
      CALL GenerateAtomLabel(iElement2, iAtom2, label2)
      WRITE(hFile, '(A,A7,1X,A7,A,F9.5,A)', ERR=999) '    Distance_Restrain(', label1, label2, &
                ', ', Distance, ', 1.0, bond_width, bond_weight)'
      TOPASWriteDistance = .FALSE.
  999 RETURN

      END FUNCTION TOPASWriteDistance
!
!*****************************************************************************
!
      LOGICAL FUNCTION TOPASWriteAngle (hFile, iAtom1, iAtom, iAtom2, &
                                        iElement1, iElement, iElement2, Angle)

      IMPLICIT NONE

      INTEGER, INTENT (IN) :: hFile, iAtom1, iAtom, iAtom2, iElement1, iElement, iElement2
      REAL, INTENT (IN) :: Angle

      CHARACTER(7) label1, Label, label2

      TOPASWriteAngle = .TRUE.
      CALL GenerateAtomLabel(iElement1, iAtom1, label1)
      CALL GenerateAtomLabel(iElement, iAtom, Label)
      CALL GenerateAtomLabel(iElement2, iAtom2, label2)
      WRITE(hFile, '(A,A7,1X,A7,1X,A7,A,F9.5,A)', ERR=999) '    Angle_Restrain(', label1, &
                    Label, label2, ', ', Angle, ', 1.0, angle_width, angle_weight)'
      TOPASWriteAngle = .FALSE.
  999 RETURN

      END FUNCTION TOPASWriteAngle
!
!*****************************************************************************
!
      LOGICAL FUNCTION TOPASWritePlane (hFile, n, members, iElements)

      IMPLICIT NONE

      INTEGER, INTENT (IN) :: hFile, n, members(n), iElements(n)

      INTEGER I
      CHARACTER*7 Label

      TOPASWritePlane = .TRUE.
      WRITE(hFile, '(A$)', ERR=999) '    Flatten('
      DO I = 1, n
        CALL GenerateAtomLabel(iElements(I), members(I), Label)
        WRITE(hFile, '(1X,A7$)', ERR=999) Label
      ENDDO
! There are two flatten macros in TOPAS3:
!   Flatten(sites, c, t_calc, tol, wscale)
!   Flatten(sites, t_calc, tol, wscale) and
! The second one is only a wrapper of the first. But it seems the arguments t_calc
! and tol were put in wrong place. Call the first one instead
!      WRITE(hFile, '(A)', ERR=999) ', 0, flatten_width, flatten_weight)'
      WRITE(hFile, '(A)', ERR=999) ',, 0, flatten_width, flatten_weight)'
      TOPASWritePlane = .FALSE.
  999 RETURN

      END FUNCTION TOPASWritePlane
!
!*****************************************************************************
!
! This function backup/restore both .inp and .out files.
      INTEGER FUNCTION BackupRestorTopasFiles(FileName, Backup)

      USE dash_gui_resources
      USE VARIABLES

      IMPLICIT NONE

      CHARACTER*(*), INTENT (IN   ) :: FileName
      LOGICAL,       INTENT (IN   ) :: Backup

!      CHARACTER*20, EXTERNAL :: Integer2String
!      CHARACTER*20 tStr
      CHARACTER*(*), PARAMETER :: cBackSuffix = '.bak'
      CHARACTER*4 tSuffix(2)
      DATA tSuffix/'.inp', '.out'/
      CHARACTER(MaxPathLength) tDirName, tFileName, tName
      INTEGER ExtLength, I, iRetCode

      BackupRestorTopasFiles = 0
      ExtLength = 3
      CALL SplitPath2(FileName, tDirName, tFileName, tName, ExtLength)
      iRetCode = InfoError(1)
      DO I = 1, 2
        tName = TRIM(tDirName)//TRIM(tFileName)//tSuffix(I)
        IF (Backup) THEN
          CALL IOsCopyFile(TRIM(tName), TRIM(tName)//cBackSuffix)
        ELSE
          CALL IOsCopyFile(TRIM(tName)//cBackSuffix, TRIM(tName))
        ENDIF
        IF (InfoError(1) .EQ. ErrOSCommand) BackupRestorTopasFiles = BackupRestorTopasFiles + 1
      ENDDO
      RETURN

      END FUNCTION BackupRestorTopasFiles
!
!*****************************************************************************
!

max_stars_repo_name: ccdc-opensource/dash

max_issues_repo_name: ccdc-opensource/dash

max_forks_repo_name: ccdc-opensource/dash


! RUN: not %flang_fc1 %s 2>&1 | FileCheck %s
! Check errors found in folding
! TODO: test others emitted from flang/lib/Evaluate
module m
 contains
  subroutine s1(a,b)
    real :: a(*), b(:)
    !CHECK: error: DIM=1 dimension is out of range for rank-1 assumed-size array
    integer :: ub1(ubound(a,1))
    !CHECK-NOT: error: DIM=1 dimension is out of range for rank-1 assumed-size array
    integer :: lb1(lbound(a,1))
    !CHECK: error: DIM=0 dimension is out of range for rank-1 array
    integer :: ub2(ubound(a,0))
    !CHECK: error: DIM=2 dimension is out of range for rank-1 array
    integer :: ub3(ubound(a,2))
    !CHECK: error: DIM=0 dimension is out of range for rank-1 array
    integer :: lb2(lbound(b,0))
    !CHECK: error: DIM=2 dimension is out of range for rank-1 array
    integer :: lb3(lbound(b,2))
  end subroutine
end module


text:
MODULE ghex_future_mod
  use iso_c_binding

  ! should be established and defined by cmake, also for request_bind.hpp
  implicit none
#include "ghex_config.h"

  type, bind(c) :: ghex_future
     integer(c_int8_t) :: data(GHEX_FUTURE_SIZE) = 0
  end type ghex_future

  type, bind(c) :: ghex_future_multi
     integer(c_int8_t) :: data(GHEX_FUTURE_MULTI_SIZE) = 0
  end type ghex_future_multi

  interface

     subroutine ghex_future_wait_single(future) bind(c)
       use iso_c_binding
       import ghex_future
       type(ghex_future) :: future
     end subroutine ghex_future_wait_single

     subroutine ghex_future_wait_multi(future) bind(c)
       use iso_c_binding
       import ghex_future_multi
       type(ghex_future_multi) :: future
     end subroutine ghex_future_wait_multi

     logical(c_bool) function ghex_future_ready_single(future) bind(c)
       use iso_c_binding
       import ghex_future
       type(ghex_future) :: future
     end function ghex_future_ready_single

     logical(c_bool) function ghex_future_ready_multi(future) bind(c)
       use iso_c_binding
       import ghex_future_multi
       type(ghex_future_multi) :: future
     end function ghex_future_ready_multi

     integer(c_int) function ghex_future_test_any_wrapped(futures, n_futures) bind(c, name="ghex_future_test_any")
       use iso_c_binding
       type(c_ptr), value :: futures
       integer(c_int), value :: n_futures
     end function ghex_future_test_any_wrapped

     logical(c_bool) function ghex_future_cancel(future) bind(c)
       use iso_c_binding
       import ghex_future
       type(ghex_future) :: future
     end function ghex_future_cancel

  end interface

  interface ghex_future_init
     procedure :: ghex_future_init_single
     procedure :: ghex_future_init_multi
  end interface ghex_future_init

  interface ghex_future_wait
     procedure :: ghex_future_wait_single
     procedure :: ghex_future_wait_multi
  end interface ghex_future_wait

  interface ghex_future_ready
     procedure :: ghex_future_ready_single
     procedure :: ghex_future_ready_multi
  end interface ghex_future_ready

CONTAINS

  subroutine ghex_future_init_single(future)
    type(ghex_future) :: future
    future%data = 0
  end subroutine ghex_future_init_single

  subroutine ghex_future_init_multi(future)
    type(ghex_future_multi) :: future
    future%data = 0
  end subroutine ghex_future_init_multi

  integer(c_int) function ghex_future_test_any(futures)
    type(ghex_future), dimension(:), target :: futures
    ghex_future_test_any = ghex_future_test_any_wrapped(c_loc(futures), size(futures))
  end function ghex_future_test_any

END MODULE ghex_future_mod

max_stars_repo_name: bettiolm/GHEX

max_issues_repo_name: bettiolm/GHEX

max_forks_repo_name: bettiolm/GHEX


text:
cdis    Forecast Systems Laboratory
cdis    NOAA/OAR/ERL/FSL
cdis    325 Broadway
cdis    Boulder, CO     80303
cdis
cdis    Forecast Research Division
cdis    Local Analysis and Prediction Branch
cdis    LAPS
cdis
cdis    This software and its documentation are in the public domain and
cdis    are furnished "as is."  The United States government, its
cdis    instrumentalities, officers, employees, and agents make no
cdis    warranty, express or implied, as to the usefulness of the software
cdis    and documentation for any purpose.  They assume no responsibility
cdis    (1) for the use of the software and documentation; or (2) to provide
cdis     technical support to users.
cdis
cdis    Permission to use, copy, modify, and distribute this software is
cdis    hereby granted, provided that the entire disclaimer notice appears
cdis    in all copies.  All modifications to this software must be clearly
cdis    documented, and are solely the responsibility of the agent making
cdis    the modifications.  If significant modifications or enhancements
cdis    are made to this software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis
cdis
cdis
cdis
cdis
cdis
      subroutine readcdf(csat_id,csat_type,chtype,
     1record,n_elem,n_lines,r4_image,scale_img,
     1latitude,longitude,
     1La1,Lo1,Dx,Dy,Latin,
     1LoV, ivalidTime , ncid, istatus)

c     written by Dan Birkenheuer February 1994
c     J Smart   5/95    Included code to deal with gvar /public files
c     J Smart   6/96    Included code to deal with WFO satellite netcdf files
c     J Smart   3/97    Included subroutine rdblock_line_elem to read only the
c			appropriate sector of visible satellite data. Argument
c                       list lets subroutine know whether logical or I*2 data type.
c     J Smart   3/97    Converted _vis routine to be the only netCDF reader. Works
c			for vis, ir, wv and sounder. Block reading allows this.
c     J Smart   4/97    Adapted code further to work with gvar. NetCDF headers for
c                       raw gvar are different than for fsl-conus (ie., no La1, Lo1,
c                       Lov, or Latin).
c     J Smart   5/07    Added capability to read FMI (Finish Met Inst) data
c=====================================================================================
      implicit none

      Include 'netcdf.inc'

      Integer n_elem,n_lines,record
      INTEGER RCODE
C
      real      r4_image(n_elem,n_lines)
      real      latitude(n_elem,n_lines)
      real      longitude(n_elem,n_lines)
      integer  ib
      integer       bi (2)
      equivalence (ib,bi (1) )
      integer istatus
      Integer   iDx
      Integer   iD
      INTEGER   Nx
      INTEGER   Ny
      INTEGER   Nz
      REAL      La1
      REAL      Lo1
      REAL      Dx
      REAL      Dy
      REAL      Latin
      REAL      LoV
      REAL      poLat
      REAL      La2,Lo2,level
      REAL      scale_img
      real arg1,arg2
      INTEGER START(10)
      INTEGER COUNT(10)
      integer varid,ncid
      integer center_id, process_id,
     +     wmo_sat_id(record),ivalidtime,nav
      Integer imax,jmax,kmax,kdim
      Integer channel_fcinv

      integer lineRes, elemRes
      common /cdfdata/ lineRes, elemRes

      double precision reftime(record), valtime(record)
      character*30  c_valtime
      character*30  c_Lov
      character*30  c_xres
      character*30  c_yres
      Character*3   csat_type
      character*3   chtype
      character*6   csat_id
      character*132 origin_name
      character*132 x_dim
      character*132 y_dim
      character*132 earth_shape
      character*132 wavelength(record)
      character*132 grid_name
      character*132 process_name
      character*132 grid_type
      character*132 channel_comment_
      character*132 asctime
c---------------------------------------------------------
c   code

      print*,'Subroutine readcdf...'
      print*,'n_elem/n_lines = ',n_elem,n_lines

      istatus = 0 ! bad istatus

      call NCPOPT(0)

      if(csat_type.eq.'rll')then ! read lat/lon
         rcode = NF_INQ_VARID(ncid,'latitude',varid)
         if(rcode.ne.NF_NOERR) then
            print *, NF_STRERROR(rcode)
            print *,'in var latitude'
         else
            write(6,*)
            Write(6,*)'Calling rdblock_line_elem - latitude'

            Call rdblock_line_elem(csat_id,csat_type,chtype,
     &          ncid,varid,n_elem,n_lines,latitude,istatus)

            if(istatus .ne. 1)then
               write(6,*)'Error in rdblock_line_elem'
               return
            endif
         endif
         call clean_nan2(latitude,n_elem,n_lines,istatus)
         arg1 = minval(latitude)
         arg2 = maxval(latitude)
         write(6,*)'readcdf latitude range:       ',arg1,arg2
         if(arg1 .eq. 0. .and. arg2 .eq. 0.)then
            write(6,*)' ERROR in latitude range'
            istatus = 0
            return
         endif

         rcode = NF_INQ_VARID(ncid,'longitude',varid)
         if(rcode.ne.NF_NOERR) then
            print *, NF_STRERROR(rcode)
            print *,'in var longitude'
         else
            write(6,*)
            Write(6,*)'Calling rdblock_line_elem - longitude'

            Call rdblock_line_elem(csat_id,csat_type,chtype,
     &          ncid,varid,n_elem,n_lines,longitude,istatus)

            if(istatus .ne. 1)then
               write(6,*)'Error in rdblock_line_elem'
               return
            endif
         endif
         call clean_nan2(longitude,n_elem,n_lines,istatus)
         write(6,*)'readcdf longitude range:      '
     1          ,minval(longitude),maxval(longitude)

      endif

      if(csat_type.eq.'ncp')then
         rcode = NF_INQ_VARID(ncid,'channel',varid)
      else
         rcode = NF_INQ_VARID(ncid,'image',varid)
      endif
      if(rcode.ne.NF_NOERR) then
         print *, NF_STRERROR(rcode)
         print *,'in var image'
      endif
C
C    statements to fill image
C
      write(6,*)
      Write(6,*)'Calling rdblock_line_elem - image read sub'

      Call rdblock_line_elem(csat_id,csat_type,chtype,
     &ncid,varid,n_elem,n_lines,r4_image,istatus)

      if(istatus .ne. 1)then
         write(6,*)'Error in rdblock_line_elem'
         return
      endif

      if(csat_type.eq.'rll')then

!         read scaling attribute
          rcode=NF_GET_ATT_REAL(ncid,varid,'scale_factor'
     1                              ,scale_img)
          if(rcode.ne.NF_NOERR) then
             write(6,*)'Error reading image scaling attribute'
             scale_img = .01
             write(6,*)' Use default value for image scale ',scale_img
          else
             write(6,*)' Successfully read image scale ',scale_img
          endif

          rcode = NF_INQ_VARID(ncid,'lineRes',varid)
          if(rcode.ne.NF_NOERR) then
             print *, NF_STRERROR(rcode)
             print *,'in var lineRes'
          endif
          rcode = NF_GET_VAR_INT(ncid,varid,lineRes)
          if(rcode.ne.NF_NOERR) then
             print *, NF_STRERROR(rcode)
             print *,'in var lineRes'
          endif
          if(rcode.ne.NF_NOERR) then
             write(6,*)'Error reading lineRes'
          else
             write(6,*)' Successfully read lineRes ',lineRes
          endif

          rcode = NF_INQ_VARID(ncid,'elemRes',varid)
          if(rcode.ne.NF_NOERR) then
             print *, NF_STRERROR(rcode)
             print *,'in var elemRes'
          endif
          rcode = NF_GET_VAR_INT(ncid,varid,elemRes)
          if(rcode.ne.NF_NOERR) then
             print *, NF_STRERROR(rcode)
             print *,'in var elemRes'
          endif
          if(rcode.ne.NF_NOERR) then
             write(6,*)'Error reading elemRes'
          else
             write(6,*)' Successfully read elemRes ',elemRes
          endif
C
      endif

      istatus = 0
c
c JSmart: 6-4-97.  WFO satellite netCDF files changed rather dramatically
c                  such that NO header info exists and we must now jump over
c                  the statements to read the header info.
c
      if(csat_type.eq.'cdf')then
         call read_netcdf_sat_head(ncid,record,
     + Nx, Ny, center_id,process_id, wmo_sat_id,Dx,Dy,
     + La1, Latin, Lo1, Lov, reftime, valtime, earth_shape,
     + grid_name, grid_type, origin_name, process_name,
     + wavelength, x_dim, y_dim, istatus)

         ivalidtime = int(valtime(1))
      endif

      write(6,*)'readcdf image range: '
     1          ,minval(r4_image),maxval(r4_image)

      istatus = 1  ! ok!

      Return
      END

ext: f

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:
        SUBROUTINE YCLCST (RESULT, MAXLEN, NUPDAT, CHAN, EQUIV, ISTAT)
*
*  Module number:
*
*  Module name: YCLCST
*
*  Keyphrase:
*  ----------
*       Polarimetry processing
*
*  Description:
*  ------------
*       Routine which recomputes all the Stokes parameters, the linear and
*       circular polarization, the phase angle, and associated errors.
*
*  FORTRAN name: YCLCST
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*       none
*
*  Subroutines Called:
*  -------------------
*  CDBS:
*       yclnrm
*
*  SDAS:
*
*
*  History:
*  --------
*  Version      Date    Author		  Description
*       1       Mar97   M.D. De La Pena   Original Implementation for CALFOS
*                                         based upon code from R. Allen.
*       1.1     Jul97   M.D. De La Pena   Mods for new parameter
*       1.2     Feb98   M.D. De La Pena   Changed I cutoff to allow negative I
*                                         for scaled values.
*-------------------------------------------------------------------------------
*
* INPUTS:
*       RESULT - Array containing Stokes parameters, linear and circular
*                polarization, polarization angle, and error estimates
*       MAXLEN - Maximum number of channels
*       NUPDAT - Number of quantities being updated
*       CHAN   - Current channel being processed
*       EQUIV  - Matrix of BLAH
*       ISTAT  - Return status
*
* OUTPUT:
*       RESULT - Modified input array
*       ISTAT  - Modified return status
*----------------------------------------------------------------------------
C
C UMSPUT DESTINATIONS -- CB, DAO, 4-SEP-87
C
      INTEGER STDOUT
      PARAMETER (STDOUT = 1)
      INTEGER STDERR
      PARAMETER (STDERR = 2)
C
C Passed Parameters
C
      DOUBLE PRECISION EQUIV(4,4)
      INTEGER MAXLEN, NUPDAT, CHAN, ISTAT
      REAL    RESULT(MAXLEN,NUPDAT)
C
C Local Variables
C
      CHARACTER*80     CONTXT
      DOUBLE PRECISION STOKIN(8), STOKOUT(8), P(14), SUM, SCALE
      DOUBLE PRECISION DZERO
      INTEGER I, J
C
      PARAMETER (DZERO = 0.0D0)
C
C----------------------------------------------------------------------------
C
C Initialize the temporary array
C
      DO 10 J = 1, NUPDAT
          P(J) = DZERO
10    CONTINUE
C
C Extract the Stokes parameters I, Q, U, V and the associated errors
C
      DO 20 J = 1, 8
          STOKIN(J) = RESULT(CHAN,J)
20    CONTINUE
C
C If the Stokes I for this channel is > 0.0, then proceed.
C Let the scaled Stokes I be negative - MDD 02/11/98
C
CMDD      IF (STOKIN(1) .GT. DZERO) THEN
C
          SCALE = 1.0D+00
          IF(STOKIN(1) .LT. 1.0D-08) SCALE = 1.0D+12
          IF(STOKIN(1) .LT. 1.0D-16) SCALE = 1.0D+20
          IF(STOKIN(1) .LT. 1.0D-24) SCALE = 1.0D+28
C
C Rescale the data
C
          DO 30 J = 1, 8
              STOKIN(J) = SCALE * STOKIN(J)
30        CONTINUE
C
C Scaled Stokes I can be negative - MDD 02/11/98
C
      IF (ABS(STOKIN(1)) .GT. 1.0D-8) THEN
C
C Recompute the Stokes I, Q, U, and V
C
          DO 50 I = 1, 4
              SUM = DZERO
              DO 40 J = 1, 4
                  SUM = SUM + EQUIV(I,J) * STOKIN(J)
40            CONTINUE
              STOKOUT(I) = SUM
50        CONTINUE
C
C Recompute the errors associated with the Stokes I, Q, U, and V
C
          DO 70 I = 1, 4
              SUM = DZERO
              DO 60 J = 1, 4
                  SUM = SUM + (EQUIV(I,J)*EQUIV(I,J)*
     *                        STOKIN(J+4)*STOKIN(J+4))
60            CONTINUE
              STOKOUT(I+4) = DSQRT(SUM)
70        CONTINUE
C
          DO 80 J = 1, 8
              P(J) = STOKOUT(J) / SCALE
80        CONTINUE
C
C Recompute the linear and circular polarization, polarization angle,
C and errors
C
          CALL YCLNRM (P, NUPDAT, CHAN, ISTAT)
          IF (ISTAT .NE. 0) THEN
              CONTXT = 'ERROR in computing the corrected polarization'//
     *                 ' values.'
              GO TO 999
          END IF
C
C Store the updated P values in the RESULT array
C
          DO 90 I = 1, NUPDAT
              RESULT(CHAN,I) = P(I)
90        CONTINUE
C
      ELSE
C
C When Stokes I is <= zero, set PL, PC, and THETA (and associated errors)
C to zero.
C
          DO 110 I = 9, NUPDAT
             RESULT(CHAN,I) = DZERO
110       CONTINUE
C
      ENDIF
C
      ISTAT = 0
      GO TO 1000
999   CALL YMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
      ISTAT = 1
1000  RETURN
      END

max_stars_repo_name: iraf-community/stsdas


subroutine sw_compute_solve_diagnostics(dt, nCells, nEdges, nVertices,        &
     maxEdges, maxEdges2, nVertLevels, vertexDegree,                          &
     config_monotonic, config_thickness_adv_order, config_apvm_upwinding,     &
     h, u, v, vh, h_edge, h_vertex, circulation,                              &
     vorticity, divergence, ke, pv_edge, pv_vertex, pv_cell, vorticity_cell,  &
     gradPVn, gradPVt, weightsOnEdge, kiteAreasOnVertex, cellsOnEdge,         &
     cellsOnVertex, cellsOnCell, verticesOnEdge, nEdgesOnCell, edgesOnCell,   &
     nEdgesOnEdge, edgesOnEdge, edgesOnVertex, dcEdge, dvEdge, areaCell,      &
     areaTriangle, fVertex, deriv_two, boundaryEdge, boundaryCell)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Compute diagnostic fields used in the tendency computations
  !
  ! Input: grid - grid metadata
  !
  ! Output: s - computed diagnostics
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  implicit none

  integer, parameter :: RKIND  = selected_real_kind(12)
  logical :: config_monotonic
  integer :: config_thickness_adv_order
  real (kind=RKIND) :: config_apvm_upwinding
  real (kind=RKIND), intent(in) :: dt
  integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
  integer :: maxEdges, maxEdges2

  real(kind=RKIND), dimension(nVertices) :: fVertex
  real(kind=RKIND), dimension(nEdges) :: dvEdge
  real(kind=RKIND), dimension(nEdges) :: dcEdge
  real(kind=RKIND), dimension(nCells) :: areaCell
  real(kind=RKIND), dimension(nVertices) :: areaTriangle
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: vh
  real(kind=RKIND), dimension(maxEdges2,nEdges) :: weightsOnEdge
  real(kind=RKIND), dimension(vertexDegree,nVertices) :: kiteAreasOnVertex
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: h_edge
  real(kind=RKIND), dimension(nVertLevels,nCells), intent(in) :: h
  real(kind=RKIND), dimension(nVertLevels,nEdges), intent(in) :: u
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: v
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: circulation
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: vorticity
  real(kind=RKIND), dimension(nVertLevels,nCells) :: ke
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: pv_edge
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: pv_vertex
  real(kind=RKIND), dimension(nVertLevels,nCells) :: pv_cell
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: gradPVn
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: gradPVt
  real(kind=RKIND), dimension(nVertLevels,nCells) :: divergence
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: h_vertex
  real(kind=RKIND), dimension(nVertLevels,nCells) :: vorticity_cell

  integer, dimension(2,nEdges) :: cellsOnEdge
  integer, dimension(vertexDegree,nVertices) :: cellsOnVertex
  integer, dimension(2,nEdges) :: verticesOnEdge
  integer, dimension(maxEdges,nCells) :: edgesOnCell
  integer, dimension(maxEdges2,nEdges) :: edgesOnEdge
  integer, dimension(vertexDegree,nVertices) :: edgesOnVertex
  integer, dimension(nVertLevels,nEdges) :: boundaryEdge
  integer, dimension(nVertLevels,nCells) :: boundaryCell
  integer, dimension(maxEdges,nCells) :: cellsOnCell

  integer, dimension(nCells) :: nEdgesOnCell
  integer, dimension(nEdges) :: nEdgesOnEdge
  real (kind=RKIND), dimension(15,2,nEdges) :: deriv_two

  ! ----- local variables -----
  real (kind=RKIND) :: r, h1, h2
  real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2

  real (kind=RKIND) :: coef_3rd_order

  integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
  real (kind=RKIND) :: flux, vorticity_abs, workpv
  ! ----- end local variables -----

  ! Find those cells that have an edge on the boundary
  boundaryCell(:,:) = 0
  do iEdge = 1, nEdges
     do k = 1, nVertLevels
        if(boundaryEdge(k,iEdge)==1) then
           cell1 = cellsOnEdge(1,iEdge)
           cell2 = cellsOnEdge(2,iEdge)
           boundaryCell(k,cell1) = 1
           boundaryCell(k,cell2) = 1
        endif
     enddo
  enddo

  ! Compute height on cell edges at velocity locations
  !   Namelist options control the order of accuracy of the reconstructed h_edge value
  coef_3rd_order = 0.
  if (config_thickness_adv_order == 3) coef_3rd_order = 1.0
  if (config_thickness_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25

  if (config_thickness_adv_order == 2) then

     do iEdge = 1, nEdges
        cell1 = cellsOnEdge(1,iEdge)
        cell2 = cellsOnEdge(2,iEdge)
        if (cell1 <= nCells .and. cell2 <= nCells) then
           do k = 1, nVertLevels
              h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
           end do
        end if
     end do

  else if (config_thickness_adv_order == 3) then

     do iEdge = 1, nEdges
        cell1 = cellsOnEdge(1,iEdge)
        cell2 = cellsOnEdge(2,iEdge)

        !-- if a cell not on the most outside ring of the halo
        if (cell1 <= nCells .and. cell2 <= nCells) then

           do k = 1, nVertLevels

              d2fdx2_cell1 = 0.0
              d2fdx2_cell2 = 0.0

              !-- if not a boundary cell
              if(boundaryCell(k,cell1) == 0 .and. boundaryCell(k,cell2) == 0) then

                 d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
                 d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)

                 !-- all edges of cell 1
                 do i = 1, nEdgesOnCell(cell1)
                    d2fdx2_cell1 = d2fdx2_cell1 + &
                         deriv_two(i+1,1,iEdge) * h(k, cellsOnCell(i,cell1))
                 end do

                 !-- all edges of cell 2
                 do i = 1, nEdgesOnCell(cell2)
                    d2fdx2_cell2 = d2fdx2_cell2 + &
                         deriv_two(i+1,2,iEdge) * h(k, cellsOnCell(i,cell2))
                 end do

              endif

              !-- if u > 0:
              if (u(k,iEdge) > 0) then
                 h_edge(k,iEdge) =     &
                      0.5*(h(k,cell1) + h(k,cell2))      &
                      -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.          &
                      -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
                 !-- else u <= 0:
              else
                 h_edge(k,iEdge) =     &
                      0.5*(h(k,cell1) + h(k,cell2))      &
                      -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.          &
                      +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
              end if

           end do   ! do k
        end if      ! if (cell1 <=
     end do         ! do iEdge

  else  if (config_thickness_adv_order == 4) then

     do iEdge = 1, nEdges
        cell1 = cellsOnEdge(1,iEdge)
        cell2 = cellsOnEdge(2,iEdge)

        !-- if a cell not on the most outside ring of the halo
        if (cell1 <= nCells .and. cell2 <= nCells) then

           do k = 1, nVertLevels

              d2fdx2_cell1 = 0.0
              d2fdx2_cell2 = 0.0

              !-- if not a boundary cell
              if(boundaryCell(k,cell1) == 0 .and. boundaryCell(k,cell2) == 0) then

                 d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
                 d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)

                 !-- all edges of cell 1
                 do i = 1, nEdgesOnCell(cell1)
                    d2fdx2_cell1 = d2fdx2_cell1 + &
                         deriv_two(i+1,1,iEdge) * h(k, cellsOnCell(i,cell1))
                 end do

                 !-- all edges of cell 2
                 do i = 1, nEdgesOnCell(cell2)
                    d2fdx2_cell2 = d2fdx2_cell2 + &
                         deriv_two(i+1,2,iEdge) * h(k, cellsOnCell(i,cell2))
                 end do

              endif

              h_edge(k,iEdge) =   &
                   0.5*(h(k,cell1) + h(k,cell2))      &
                   -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.

           end do   ! do k
        end if      ! if (cell1 <=
     end do         ! do iEdge

  endif   ! if(config_thickness_adv_order == 2)

  ! set the velocity in the nEdges+1 slot to zero, this is a dummy address
  !    used to when reading for edges that do not exist
  !u(:,nEdges+1) = 0.0

  ! Compute circulation and relative vorticity at each vertex
  circulation(:,:) = 0.0
  do iEdge = 1, nEdges
     do k = 1, nVertLevels
        circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
        circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
     end do
  end do
  do iVertex = 1, nVertices
     do k = 1, nVertLevels
        vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
     end do
  end do

  ! Compute the divergence at each cell center
  divergence(:,:) = 0.0
  do iEdge = 1, nEdges
     cell1 = cellsOnEdge(1,iEdge)
     cell2 = cellsOnEdge(2,iEdge)
     if (cell1 <= nCells) then
        do k = 1, nVertLevels
           divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
        enddo
     endif
     if(cell2 <= nCells) then
        do k = 1, nVertLevels
           divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
        enddo
     end if
  end do
  do iCell = 1, nCells
     r = 1.0 / areaCell(iCell)
     do k = 1, nVertLevels
        divergence(k,iCell) = divergence(k,iCell) * r
     enddo
  enddo

  ! Compute kinetic energy in each cell
  ke(:,:) = 0.0
  do iCell = 1, nCells
     do i = 1, nEdgesOnCell(iCell)
        iEdge = edgesOnCell(i,iCell)
        do k = 1, nVertLevels
           ke(k,iCell) = ke(k,iCell) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
        end do
     end do
     do k = 1, nVertLevels
        ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
     end do
  end do

  ! Compute v (tangential) velocities
  v(:,:) = 0.0
  do iEdge = 1,nEdges
     do i = 1, nEdgesOnEdge(iEdge)
        eoe = edgesOnEdge(i,iEdge)
        do k = 1,nVertLevels
           v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
        end do
     end do
  end do

  ! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
  !vh(:,:) = 0.0
  !do iEdge = 1, nEdges
  !   do j = 1, nEdgesOnEdge(iEdge)
  !      eoe = edgesOnEdge(j,iEdge)
  !      do k = 1, nVertLevels
  !         vh(k,iEdge) = vh(k,iEdge) + weightsOnEdge(j,iEdge) * u(k,eoe) * h_edge(k,eoe)
  !      end do
  !   end do
  !end do

  ! Compute height at vertices, pv at vertices, and average pv to edge locations
  !  ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
  do iVertex = 1,nVertices
     do k = 1, nVertLevels
        h_vertex(k,iVertex) = 0.0
        do i = 1, vertexDegree
           h_vertex(k,iVertex) = h_vertex(k,iVertex) + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
        end do
        h_vertex(k,iVertex) = h_vertex(k,iVertex) / areaTriangle(iVertex)

        pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex(k,iVertex)
     end do
  end do

  ! Compute gradient of PV in the tangent direction
  !   ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
  do iEdge = 1, nEdges
     do k = 1, nVertLevels
        gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
             dvEdge(iEdge)
     enddo
  enddo

  ! Compute pv at the edges
  !   ( this computes pv_edge at all edges bounding real cells )
  pv_edge(:,:) = 0.0
  do iVertex = 1,nVertices
     do i = 1, vertexDegree
        iEdge = edgesOnVertex(i,iVertex)
        do k = 1, nVertLevels
           pv_edge(k,iEdge) =  pv_edge(k,iEdge)  + 0.5 * pv_vertex(k,iVertex)
        end do
     end do
  end do

  ! Modify PV edge with upstream bias.
  do iEdge = 1, nEdges
     do k = 1, nVertLevels
        pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * v(k,iEdge) * dt * gradPVt(k,iEdge)
     enddo
  enddo

  ! Compute pv at cell centers
  !    ( this computes pv_cell for all real cells and distance-1 ghost cells )
  pv_cell(:,:) = 0.0
  vorticity_cell(:,:) = 0.0
  do iVertex = 1, nVertices
     do i = 1, vertexDegree
        iCell = cellsOnVertex(i,iVertex)
        if (iCell <= nCells) then
           do k = 1, nVertLevels
              pv_cell(k,iCell) = pv_cell(k,iCell) &
                   + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) &
                   / areaCell(iCell)
              vorticity_cell(k,iCell) = vorticity_cell(k,iCell) &
                   + kiteAreasOnVertex(i, iVertex) * vorticity(k, iVertex) &
                   / areaCell(iCell)
           enddo
        endif
     enddo
  enddo

  ! Compute gradient of PV in normal direction
  !   ( this computes gradPVn for all edges bounding real cells )
  gradPVn(:,:) = 0.0
  do iEdge = 1, nEdges
     if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
        do k = 1, nVertLevels
           gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
                dcEdge(iEdge)
        enddo
     endif
  enddo

  ! Modify PV edge with upstream bias.
  do iEdge = 1, nEdges
     do k = 1, nVertLevels
        pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * u(k,iEdge) * dt * gradPVn(k,iEdge)
     enddo
  enddo

end subroutine sw_compute_solve_diagnostics

subroutine compute_mesh_scaling(nCells, nEdges, config_h_ScaleWithMesh, &
     cellsOnEdge, meshDensity, meshScalingDel2, meshScalingDel4)
  implicit none

  integer, parameter :: RKIND  = selected_real_kind(12)
  logical :: config_h_ScaleWithMesh
  integer :: nEdges, nCells
  integer, dimension(2,nEdges) :: cellsOnEdge
  real(kind=RKIND), dimension(nCells) :: meshDensity
  real(kind=RKIND), dimension(nEdges) :: meshScalingDel2
  real(kind=RKIND), dimension(nEdges) :: meshScalingDel4

  ! ----- local variable -----
  integer :: iEdge, cell1, cell2
  ! ----- end local variable -----

  ! Compute the scaling factors to be used in the del2 and del4 dissipation
  meshScalingDel2(:) = 1.0
  meshScalingDel4(:) = 1.0
  if (config_h_ScaleWithMesh) then
     do iEdge = 1, nEdges
        cell1 = cellsOnEdge(1,iEdge)
        cell2 = cellsOnEdge(2,iEdge)
        meshScalingDel2(iEdge) = 1.0 &
             / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**(5.0/12.0)
        meshScalingDel4(iEdge) = 1.0 &
             / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**(5.0/6.0)
     end do
  end if

end subroutine compute_mesh_scaling

subroutine sw_compute_tend(nCells, nEdges, nVertices, nVertLevels, &
     maxEdges, maxEdges2, vertexDegree, &
     config_bottom_drag, config_wind_stress, config_h_mom_eddy_visc2, &
     config_h_mom_eddy_visc4, h, u, v, h_edge, circulation, vorticity, &
     divergence, ke, pv_edge, vh, weightsOnEdge, kiteAreasOnVertex, &
     cellsOnEdge, cellsOnVertex, verticesOnEdge, nEdgesOnCell, edgesOnCell, &
     nEdgesOnEdge, edgesOnEdge, edgesOnVertex, dcEdge, dvEdge, &
     areaCell, areaTriangle, h_s, fVertex, fEdge, u_src, &
     meshScalingDel2, meshScalingDel4, tend_h, tend_u)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Compute height and normal wind tendencies, as well as diagnostic variables
  ! Input: s - current model state
  !        grid - grid metadata
  ! Output: tend - computed tendencies for prognostic variables
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  implicit none

  integer, parameter :: RKIND  = selected_real_kind(12)
  real (kind=RKIND), parameter :: gravity = 9.80616
  integer :: nCells, nEdges, nVertices, nVertLevels
  integer :: maxEdges, maxEdges2, vertexDegree
  real(kind=RKIND), dimension(nCells) :: h_s
  real(kind=RKIND), dimension(nVertices) :: fVertex
  real(kind=RKIND), dimension(nEdges) :: fEdge
  real(kind=RKIND), dimension(nEdges) :: dvEdge
  real(kind=RKIND), dimension(nEdges) :: dcEdge
  real(kind=RKIND), dimension(nCells) :: areaCell
  real(kind=RKIND), dimension(nVertices) :: areaTriangle
  real(kind=RKIND), dimension(nEdges) :: meshScalingDel2
  real(kind=RKIND), dimension(nEdges) :: meshScalingDel4
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: vh
  real(kind=RKIND), dimension(maxEdges2,nEdges) :: weightsOnEdge
  real(kind=RKIND), dimension(vertexDegree,nVertices) :: kiteAreasOnVertex
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: h_edge
  real(kind=RKIND), dimension(nVertLevels,nCells) :: h
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: u
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: v
  real(kind=RKIND), dimension(nVertLevels,nCells) :: tend_h
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: tend_u
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: circulation
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: vorticity
  real(kind=RKIND), dimension(nVertLevels,nCells) :: ke
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: pv_edge
  real(kind=RKIND), dimension(nVertLevels,nCells) :: divergence
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: h_vertex

  integer, dimension(2,nEdges) :: cellsOnEdge
  integer, dimension(vertexDegree,nVertices) :: cellsOnVertex
  integer, dimension(2,nEdges) :: verticesOnEdge
  integer, dimension(maxEdges,nCells) :: edgesOnCell
  integer, dimension(maxEdges2,nEdges) :: edgesOnEdge
  integer, dimension(vertexDegree,nVertices) :: edgesOnVertex
  integer, dimension(nCells) :: nEdgesOnCell
  integer, dimension(nEdges) :: nEdgesOnEdge

  real (kind=RKIND), dimension(nVertLevels, nEdges) :: u_src

  logical :: config_wind_stress, config_bottom_drag
  real (kind=RKIND) :: config_h_mom_eddy_visc2, config_h_mom_eddy_visc4

  ! ----- local variables -----
  integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, timeLevel
  real (kind=RKIND) :: flux, vorticity_abs, workpv, q, upstream_bias
  real (kind=RKIND) :: ke_edge
  real (kind=RKIND), parameter :: rho_ref = 1000.0
  real (kind=RKIND) :: r, u_diffusion
  real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
  real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
  real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
  ! ----- end local variables

  ! Compute height tendency for each cell
  tend_h(:,:) = 0.0
  do iEdge = 1, nEdges
     cell1 = cellsOnEdge(1,iEdge)
     cell2 = cellsOnEdge(2,iEdge)
     do k = 1, nVertLevels
        flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
        tend_h(k,cell1) = tend_h(k,cell1) - flux
        tend_h(k,cell2) = tend_h(k,cell2) + flux
     end do
  end do
  do iCell = 1, nCells
     do k = 1, nVertLevels
        tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
     end do
  end do

  ! Compute u (normal) velocity tendency for each edge (cell face)
  tend_u(:,:) = 0.0
  do iEdge = 1, nEdges
     cell1 = cellsOnEdge(1,iEdge)
     cell2 = cellsOnEdge(2,iEdge)
     vertex1 = verticesOnEdge(1,iEdge)
     vertex2 = verticesOnEdge(2,iEdge)

     do k = 1, nVertLevels
        q = 0.0
        do j = 1, nEdgesOnEdge(iEdge)
           eoe = edgesOnEdge(j,iEdge)
           workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
           q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
        end do

        tend_u(k,iEdge) = q     &
             - (   ke(k,cell2) - ke(k,cell1) + &
             gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
             ) / dcEdge(iEdge)
     end do
  end do

  ! Compute diffusion, computed as \nabla divergence - k \times \nabla vorticity
  !                    only valid for visc == constant
  if (config_h_mom_eddy_visc2 > 0.0) then
     do iEdge = 1, nEdges
        cell1 = cellsOnEdge(1,iEdge)
        cell2 = cellsOnEdge(2,iEdge)
        vertex1 = verticesOnEdge(1,iEdge)
        vertex2 = verticesOnEdge(2,iEdge)

        do k = 1, nVertLevels
           u_diffusion =   ( divergence(k,cell2)  -  divergence(k,cell1) ) / dcEdge(iEdge) &
                -(vorticity(k,vertex2)  - vorticity(k,vertex1) ) / dvEdge(iEdge)
           u_diffusion = meshScalingDel2(iEdge) * config_h_mom_eddy_visc2 * u_diffusion
           tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
        end do
     end do
  end if

  ! velocity tendency: del4 dissipation, -\nu_4 \nabla^4 u
  !   computed as \nabla^2 u = \nabla divergence + k \times \nabla vorticity
  !   applied recursively.
  !   strictly only valid for h_mom_eddy_visc4 == constant
  if (config_h_mom_eddy_visc4 > 0.0) then
     allocate(delsq_divergence(nVertLevels, nCells+1))
     allocate(delsq_u(nVertLevels, nEdges+1))
     allocate(delsq_circulation(nVertLevels, nVertices+1))
     allocate(delsq_vorticity(nVertLevels, nVertices+1))

     delsq_u(:,:) = 0.0

     ! Compute \nabla^2 u = \nabla divergence + k \times \nabla vorticity
     do iEdge = 1, nEdges
        cell1 = cellsOnEdge(1,iEdge)
        cell2 = cellsOnEdge(2,iEdge)
        vertex1 = verticesOnEdge(1,iEdge)
        vertex2 = verticesOnEdge(2,iEdge)

        do k = 1, nVertLevels

           delsq_u(k,iEdge) = ( divergence(k,cell2)  - divergence(k,cell1) ) / dcEdge(iEdge)  &
                -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)

        end do
     end do

     ! vorticity using \nabla^2 u
     delsq_circulation(:,:) = 0.0
     do iEdge = 1, nEdges
        vertex1 = verticesOnEdge(1,iEdge)
        vertex2 = verticesOnEdge(2,iEdge)
        do k=1,nVertLevels
           delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) &
                - dcEdge(iEdge) * delsq_u(k,iEdge)
           delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
                + dcEdge(iEdge) * delsq_u(k,iEdge)
        end do
     end do
     do iVertex = 1, nVertices
        r = 1.0 / areaTriangle(iVertex)
        do k = 1, nVertLevels
           delsq_vorticity(k,iVertex) = delsq_circulation(k,iVertex) * r
        end do
     end do

     ! Divergence using \nabla^2 u
     delsq_divergence(:,:) = 0.0
     do iEdge = 1, nEdges
        cell1 = cellsOnEdge(1,iEdge)
        cell2 = cellsOnEdge(2,iEdge)
        do k=1,nVertLevels
           delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &
                + delsq_u(k,iEdge)*dvEdge(iEdge)
           delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
                - delsq_u(k,iEdge)*dvEdge(iEdge)
        end do
     end do
     do iCell = 1,nCells
        r = 1.0 / areaCell(iCell)
        do k = 1, nVertLevels
           delsq_divergence(k,iCell) = delsq_divergence(k,iCell) * r
        end do
     end do

     ! Compute - \kappa \nabla^4 u
     ! as  \nabla div(\nabla^2 u) + k \times \nabla ( k \cross curl(\nabla^2 u) )
     do iEdge = 1, nEdges
        cell1 = cellsOnEdge(1,iEdge)
        cell2 = cellsOnEdge(2,iEdge)
        vertex1 = verticesOnEdge(1,iEdge)
        vertex2 = verticesOnEdge(2,iEdge)

        do k = 1, nVertLevels

           u_diffusion = (  delsq_divergence(k,cell2) &
                - delsq_divergence(k,cell1) ) / dcEdge(iEdge)  &
                -(  delsq_vorticity(k,vertex2) &
                - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)

           u_diffusion = meshScalingDel4(iEdge) * config_h_mom_eddy_visc4 * u_diffusion
           tend_u(k,iEdge) = tend_u(k,iEdge) - u_diffusion

        end do
     end do

     deallocate(delsq_divergence)
     deallocate(delsq_u)
     deallocate(delsq_circulation)
     deallocate(delsq_vorticity)

  end if

  ! Compute u (velocity) tendency from wind stress (u_src)
  if(config_wind_stress) then
     do iEdge = 1, nEdges
        tend_u(1,iEdge) =  tend_u(1,iEdge) &
             + u_src(1,iEdge) / rho_ref / h_edge(1,iEdge)
     end do
  endif

  if (config_bottom_drag) then
     do iEdge = 1,  nEdges
        ! bottom drag is the same as POP:
        ! -c |u| u  where c is unitless and 1.0e-3.
        ! see POP Reference guide, section 3.4.4.
        ke_edge = 0.5 * ( ke(1,cellsOnEdge(1,iEdge)) &
             + ke(1,cellsOnEdge(2,iEdge)))
        tend_u(1,iEdge) = tend_u(1,iEdge) &
             - 1.0e-3*u(1,iEdge)*sqrt(2.0*ke_edge)/h_edge(1,iEdge)
     end do
  endif

end subroutine sw_compute_tend