text:
      function find_fl_item(fl_start,command,command_list,nlist,
     &                      label,label_res,backward)

      implicit none

      include 'opdim.h'
      include 'def_contraction.h'
      include 'def_operator.h'
      include 'def_formula_item.h'

      type(formula_item), pointer :: find_fl_item
      type(formula_item), pointer :: fl_start
      integer, intent(in), optional :: command
      integer, intent(in), optional :: command_list(*), nlist
      character(len=len_opname), intent(in), optional :: label
      character(len=len_opname), intent(in), optional :: label_res
      logical, intent(in), optional :: backward

      logical ::
     &    back, ok
      integer, pointer ::
     &    comm(:)
      integer ::
     &    icom, jcom, ncom
      type(formula_item), pointer ::
     &       fl_pnt

      back = .false.
      if (present(backward)) back = backward
      ncom = 1
      if (present(nlist)) ncom = nlist
      allocate(comm(ncom))
      comm(1) = command_end_of_formula
      if (present(command)) comm(1) = command
      if (present(command_list)) comm(1:ncom) = command_list(1:ncom)

      fl_pnt => fl_start

      do
        if (.not.back) fl_pnt => fl_pnt%next
        if (     back) fl_pnt => fl_pnt%prev
        if (.not.associated(fl_pnt)) exit

        ok = .false.
        do icom = 1, ncom
          if (fl_pnt%command.eq.comm(icom)) jcom = icom
          ok = ok.or.(fl_pnt%command.eq.comm(icom))
        end do
        if (.not.ok) cycle

        if (present(label)) then
          ok = .false.
          select case (comm(jcom))
          case (command_bc,command_add_bc,
     &          command_bc_reo,command_add_bc_reo)
            ok = (trim(fl_pnt%bcontr%label_op1).eq.trim(label).or.
     &            trim(fl_pnt%bcontr%label_op2).eq.trim(label).or.
     &            trim(fl_pnt%bcontr%label_res).eq.trim(label))
          case (command_add_intm,command_cp_intm)
            ok = (trim(fl_pnt%bcontr%label_op1).eq.trim(label).or.
     &            trim(fl_pnt%bcontr%label_res).eq.trim(label))
          case (command_new_intermediate)
            ok = (trim(fl_pnt%interm%name).eq.trim(label))
          case (command_del_intermediate)
            ok = (trim(fl_pnt%label).eq.trim(label))
          case (command_reorder)
            ok = (trim(fl_pnt%reo%label_in ).eq.trim(label).or.
     &            trim(fl_pnt%reo%label_out).eq.trim(label))
          end select
          if (ok) exit
        else if (present(label_res)) then
          ok = .false.
          select case (comm(jcom))
          case (command_bc,command_add_bc,
     &          command_bc_reo,command_add_bc_reo)
            ok = trim(fl_pnt%bcontr%label_res).eq.trim(label_res)
          case (command_add_intm,command_cp_intm)
            ok = trim(fl_pnt%bcontr%label_res).eq.trim(label_res)
          case (command_reorder)
            ok = trim(fl_pnt%reo%label_out).eq.trim(label_res)
          end select
          if (ok) exit
        else
          exit
        end if

      end do

      find_fl_item => fl_pnt

      deallocate(comm)

      return
      end function

ext: f

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
      subroutine prtsca(ctb,jflag,jsflag,mrmlra,mwtsp,narn1,narn2,
     $ nbasp,nbaspd,nbt,nbtmax,nelect,nhydr,nhydx,noutpt,no2gaq,
     $ nstmax,qrho,rho,uspec,wfh2o)
c
c     This subroutine computes and prints a table of the sensible
c     composition of the aqueous solution in terms of mass balance
c     totals for component (basis) species that generally correspond to
c     analyzable solutes. This set of basis species is generally not
c     identical to either the active basis set or the 'd' basis set,
c     apart from the fact that it doesn't include the solvent, water.
c
c     This subroutine is called by:
c
c       EQ3NR/scripx.f
c       EQ6/scripz.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c       ctb    = array of molalities of the component (basis) species
c       qrho   = .true. if a solution density value is available
c       rho    = the density of the aqueous solution
c       uspec  = array of names of aqueous species
c
c     Principal output:
c
c       None
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer nbtmax,nstmax
c
      integer noutpt
c
      integer jflag(nstmax),jsflag(nstmax),nbasp(nbtmax),nbaspd(nbtmax)
      integer narn1,narn2,nbt,nelect,nhydr,nhydx,no2gaq
c
      logical qrho
c
      character(len=48) uspec(nstmax)
c
      real(8) ctb(nbtmax),mwtsp(nstmax)
      real(8) mrmlra,rho,wfh2o
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer nb,ns1,ns2
c
      real(8) ctbh,ctboh,molrx,ppmv,ppmw
c
c-----------------------------------------------------------------------
c
      write (noutpt,1000)
 1000 format(/11x,'--- Sensible Composition of the',
     $ ' Aqueous Solution ---')
c
      if (qrho) then
c
c       Have density data, include mg/L and Molarity.
c
        write (noutpt,1010)
 1010   format(/3x,'Species',20x,'mg/L',7x,'mg/kg.sol',4x,
     $  'Molarity',5x,'Molality',/)
c
        do nb = 1,nbt
          ns1 = nbaspd(nb)
          ns2 = nbasp(nb)
          if (ns1.ge.narn1 .and. ns1.le.narn2) then
            if (jsflag(ns2).lt.2 .and. jflag(ns2).ne.30) then
              if (ns1.ne.narn1 .and. ns1.ne.nelect .and.
     $          ns1.ne.no2gaq) then
                if (ns1.eq.nhydr .or. ns1.eq.nhydx) then
                  if (ctb(nb) .ge. 0.) then
                    ppmw = 1000.*ctb(nb)*mwtsp(ns1)*wfh2o
                    ppmv = ppmw*rho
                    molrx = ctb(nb)*mrmlra
                    write (noutpt,1020) uspec(ns1),ppmv,ppmw,
     $              molrx,ctb(nb)
 1020               format(1x,a24,1x,1pe12.5,1x,1pe12.5,1x,1pe12.5,
     $              1x,1pe12.5)
                  elseif (ns1 .eq. nhydr) then
                    ctboh = -ctb(nb)
                    ppmw = 1000.*ctboh*mwtsp(nhydx)*wfh2o
                    ppmv = ppmw*rho
                    molrx = ctb(nb)*mrmlra
                    write (noutpt,1020) uspec(nhydx),ppmv,ppmw,
     $              molrx,ctboh
                  else
                    ctbh = -ctb(nb)
                    ppmw = 1000.*ctbh*mwtsp(nhydr)*wfh2o
                    ppmv = ppmw*rho
                    molrx = ctb(nb)*mrmlra
                    write (noutpt,1020) uspec(nhydr),ppmv,ppmw,
     $              molrx,ctbh
                  endif
                else
                  ppmw = 1000.*ctb(nb)*mwtsp(ns1)*wfh2o
                  ppmv = ppmw*rho
                  molrx = ctb(nb)*mrmlra
                  write (noutpt,1020) uspec(ns1),ppmv,ppmw,
     $            molrx,ctb(nb)
                endif
              endif
            endif
          endif
        enddo
      else
c
c       Have no density data, so do not include mg/L and Molarity.
c
        write (noutpt,1030)
 1030   format(/3x,'Species',18x,'mg/kg.sol',4x,'Molality',/)
c
        do nb = 1,nbt
          ns1 = nbaspd(nb)
          ns2 = nbasp(nb)
          if (ns1.ge.narn1 .and. ns1.le.narn2) then
            if (jsflag(ns2).lt.2 .and. jflag(ns2).ne.30) then
              if (ns1.ne.narn1 .and. ns1.ne.nelect .and.
     $          ns1.ne.no2gaq) then
                if (ns1.eq.nhydr .or. ns1.eq.nhydx) then
                  if (ctb(nb) .ge. 0.) then
                    ppmw = 1000.*ctb(nb)*mwtsp(ns1)*wfh2o
                    write (noutpt,1040) uspec(ns1),ppmw,ctb(nb)
 1040               format(1x,a24,1x,1pe12.5,1x,1pe12.5)
                  elseif (ns1 .eq. nhydr) then
                    ctboh = -ctb(nb)
                    ppmw = 1000.*ctboh*mwtsp(nhydx)*wfh2o
                    write (noutpt,1040) uspec(nhydx),ppmw,ctboh
                  else
                    ctbh = -ctb(nb)
                    ppmw = 1000.*ctbh*mwtsp(nhydr)*wfh2o
                    write (noutpt,1040) uspec(nhydr),ppmw,ctbh
                  endif
                else
                  ppmw = 1000.*ctb(nb)*mwtsp(ns1)*wfh2o
                  write (noutpt,1040) uspec(ns1),ppmw,ctb(nb)
                endif
              endif
            endif
          endif
        enddo
      endif
c
      write (noutpt,1050)
 1050 format(/3x,'The above data have physical significance, but',
     $ ' some may be',/3x,'inconsistent with certain analytical',
     $ ' methods or reporting schemes.',/)
c
      end


text:
!
! Copyright (C) 2001-2011 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!----------------------------------------------------------------------------
PROGRAM pwscf
  !----------------------------------------------------------------------------
  !
  ! ... Plane Wave Self-Consistent Field code
  !
  USE io_global,        ONLY : stdout, ionode, ionode_id
  USE parameters,       ONLY : ntypx, npk, lmaxx
  USE cell_base,        ONLY : fix_volume
  USE control_flags,    ONLY : conv_elec, lscf
  USE control_flags,    ONLY : conv_ions, istep, nstep, restart, lmd, lbfgs
  USE force_mod,        ONLY : lforce, lstres, sigma
  USE environment,      ONLY : environment_start
  USE check_stop,       ONLY : check_stop_init
  USE mp_global,        ONLY : mp_startup, mp_global_end, intra_image_comm
  USE mp_global,        ONLY : nimage, me_image, root_image, my_image_id
  USE io_files,           ONLY : tmp_dir
  USE image_io_routines,  ONLY : io_image_start
  USE read_input,         ONLY : read_input_file
  !
  IMPLICIT NONE
  !
#ifdef __OPENMP
  include 'omp_lib.h'
#endif
  !
  INTEGER :: ierr
  CHARACTER(len=256) :: dirname
  !
#ifdef __HPCTK
  call hpctoolkit_sampling_stop();
#endif
  !
  !
  CALL mp_startup ( )
  ! reset IO nodes
  ! (do this to make each "image head node" an ionode)
  ! Has to be used ONLY to run nimage copies of pwscf
  !
#ifdef __OPENMP
  call dfftw_init_threads( omp_get_max_threads() )
  call dfftw_plan_with_nthreads( omp_get_max_threads() )
#endif
  !
  IF ( nimage > 1 ) CALL io_image_start( )
  CALL environment_start ( 'MiniDFT' )
  !
  IF ( ionode ) WRITE( unit = stdout, FMT = 9010 ) &
         ntypx, npk, lmaxx
  !
  !
  ! ... open, read, close input file
  !call MQEoptions_read()
  CALL read_input_file ('PW')
  !
  ! ... convert to internal variables
  !
  CALL iosys()
  !
  !
  IF( nimage > 1 ) THEN
     !
     ! ... When nimage are used, open a directory for each one
     ! ...It has to be done here in order not to disturb NEB like calculations
     !
     WRITE( dirname, FMT = '( I5.5 )' ) my_image_id
     tmp_dir = TRIM( tmp_dir )//TRIM( dirname )//'/'
     !
  END IF
  !
  !
  CALL check_stop_init()
  !
  !
  CALL setup ()
  !
  !
  CALL init_run()
  !
  CALL MPI_Barrier( intra_image_comm, ierr )
  CALL start_clock( 'Benchmark_Time' )
#ifdef __HPCTK
  call hpctoolkit_sampling_start();
#endif
#ifdef __IPM
  call MPI_Pcontrol( 1, "Benchmark_Time"//char(0))
#endif
  !
  main_loop: DO
     !
     ! ... electronic self-consistentcy
     !
     CALL electrons()
     !
     IF ( .NOT. conv_elec ) THEN
       CALL punch( 'all' )
       CALL stop_run( conv_elec )
     ENDIF
     !
     ! ... ionic section starts here
     !
     CALL start_clock( 'ions' )
     conv_ions = .TRUE.
     !
     ! ... force calculation
     !
     IF ( lforce ) CALL forces()
     !
     ! ... stress calculation
     !
     IF ( lstres ) CALL stress ( sigma )
     !
     !
     CALL stop_clock( 'ions' )
     !
     ! ... exit condition (ionic convergence) is checked here
     !
     IF ( conv_ions ) EXIT main_loop
     !
     ! ... terms of the hamiltonian depending upon nuclear positions
     ! ... are reinitialized here
     !
     !
  END DO main_loop
  CALL MPI_Barrier( intra_image_comm, ierr )
  CALL stop_clock( 'Benchmark_Time' )
  !
#ifdef __IPM
  call MPI_Pcontrol( -1, "Benchmark_Time"//char(0))
#endif
#ifdef __HPCTK
  call hpctoolkit_sampling_stop();
#endif
  !
  ! ... save final data file
  !
  CALL punch('all')
  CALL stop_run( conv_ions )
#ifdef __OPENMP
  call dfftw_cleanup_threads()
#endif
  !
  !  END IF
  !
  STOP
  !
9010 FORMAT( /,5X,'Current dimensions of program MiniDFT are:', &
           & /,5X,'Max number of different atomic species (ntypx) = ',I2,&
           & /,5X,'Max number of k-points (npk) = ',I6,&
           & /,5X,'Max angular momentum in pseudopotentials (lmaxx) = ',i2)
  !
END PROGRAM pwscf

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
!##############################################################################
!# ****************************************************************************
!# <name> bouss2dmini_mit_2 </name>
!# ****************************************************************************
!#
!# <purpose>
!# Todo
!#
!# System:
!#
!# / A(u) 0    B1    M1   \   / u_1 \
!# | 0    A(u) B2    M2   | * | u_2 | = 0
!# | B1^T B2^T 0     0    |   |  p  |
!# \ 0    0    0     N(u) /   \  T  /
!#
!# u = 0 on the complete boundary
!# T =  1/2 on the left boundary region
!# T = -1/2 on the right boundary region
!# </purpose>
!##############################################################################

module bouss2dmini_mit_2

  use fsystem
  use storage
  use genoutput
  use boundary
  use cubature
  use derivatives
  use matrixfilters
  use vectorfilters
  use bcassembly
  use triangulation
  use element
  use spatialdiscretisation
  use linearalgebra
  use linearsystemscalar
  use linearsystemblock
  use spdiscprojection
  use scalarpde
  use trilinearformevaluation
  use bilinearformevaluation
  use linearformevaluation
  use stdoperators
  use discretebc
  use filtersupport
  use coarsegridcorrection
  use multilevelprojection
  use linearsolver
  use ucd
  use convection
  use statistics

  use bouss2dmini_callback

  implicit none

!<types>

!<typeblock description="Type block defining all information about one level">

  type t_level

    ! An object for saving the triangulation on the domain
    type(t_triangulation) :: rtriangulation

    ! A block discretisation for the primary system (Q1~/Q1~/Q0)
    type(t_blockDiscretisation) :: rdiscrPrimary

    ! A block discretisation for the secondary system (Q1)
    type(t_blockDiscretisation) :: rdiscrSecondary

    ! A primary system matrix which recieves the Navier-Stokes system
    ! for that specific level.
    type(t_matrixBlock) :: rmatrixPrimary

    ! A secondary system matrix which recieves the Convection-Diffusion system
    ! of the temperature for that specific level.
    type(t_matrixBlock) :: rmatrixSecondary

    ! Laplace-matrix of the primary system for that specific level.
    type(t_matrixScalar) :: rmatrixA

    ! B1-matrix for that specific level.
    type(t_matrixScalar) :: rmatrixB1

    ! B2-matrix for that specific level.
    type(t_matrixScalar) :: rmatrixB2

    ! M1- / M2-matrix for that specific level.
    type(t_matrixScalar) :: rmatrixM

    ! Laplace-Matrix of the secondary system for that specific level.
    type(t_matrixScalar) :: rmatrixN

    ! Primary solution vector for that specific level.
    type(t_vectorBlock) :: rvecSolPrimary

    ! Secondary solution vector for that specific level.
    type(t_vectorBlock) :: rvecSolSecondary

    ! A variable describing the discrete boundary conditions for the
    ! primary system.
    type(t_discreteBC) :: rdiscreteBCprimary

    ! A variable describing the discrete boundary conditions for the
    ! secondary system.
    type(t_discreteBC) :: rdiscreteBCsecondary

  end type

!</typeblock>

!</types>

contains

  ! ***************************************************************************

!<subroutine>

  subroutine b2dm_mit_2

!<description>
  ! This is an all-in-one stokes solver for directly solving a stokes
  ! problem without making use of special features like collections
  ! and so on. The routine performs the following tasks:
  !
  ! 1.) Read in parametrisation
  ! 2.) Read in triangulation
  ! 3.) Set up RHS
  ! 4.) Set up matrix
  ! 5.) Create solver structure
  ! 6.) Solve the problem
  ! 7.) Write solution to GMV file
  ! 8.) Release all variables, finish
!</description>

!</subroutine>

    ! Definitions of variables.
    !
    ! We need a couple of variables for this problem. Let's see...
    !
    ! An array of problem levels for the multigrid solver
    type(t_level), dimension(:), pointer :: Rlevels

    ! An object for saving the domain:
    type(t_boundary) :: rboundary

    ! A trilinear form describing the convection.
    type(t_trilinearform) :: rtriform1, rtriform2

    ! A streamline diffusion object for the convection.
    type(t_convStreamlineDiffusion) :: rsd

    ! A block matrix and a couple of block vectors. These will be filled
    ! with data for the linear solver.
    type(t_vectorBlock) :: rrhsPrimary,rrhsSecondary,&
        rvecTmpPrimary,rvecTmpSecondary,rvecDefPrimary,rvecDefSecondary

    ! A set of variables describing the analytic and discrete boundary
    ! conditions.
    type(t_boundaryRegion) :: rboundaryRegion

    ! A solver node that accepts parameters for the linear solver
    type(t_linsolNode), pointer :: p_rsolverPrimary,p_rsolverSecondary,&
        p_rcoarseGridSolver,p_rsmoother,p_rpreconditioner

    ! An array for the system matrix(matrices) during the initialisation of
    ! the linear solver.
    type(t_matrixBlock), dimension(:), pointer :: Rmatrices

    ! A filter chain that describes how to filter the matrix/vector
    ! before/during the solution process. The filters usually implement
    ! boundary conditions.
    type(t_filterChain), dimension(1), target :: RfilterChain
    type(t_filterChain), dimension(:), pointer :: p_RfilterChain

    ! An interlevel projection structure for changing levels
    type(t_interlevelProjectionBlock) :: rprojPrimary, rprojSecondary

    ! One level of multigrid
    type(t_linsolMGLevelInfo), pointer :: p_rlevelInfo

    ! NLMIN receives the level of the coarse grid.
    integer :: NLMIN

    ! NLMAX receives the level where we want to solve.
    integer :: NLMAX

    ! Viscosity parameter nu = 1/Re
    real(DP) :: dnu

    ! Thermal conductivity eta
    real(DP) :: deta

    ! Gravity direction (dgrav1, dgrav2)
    real(DP) :: dgrav1, dgrav2

    ! Non-linear damping parameters
    real(DP) :: dnlDampPrimary, dnlDampSecondary

    ! Error indicator during initialisation of the solver
    integer :: ierror

    ! A counter variable for the non-linear loop
    integer :: NLiter

    ! Non-linear residuals
    real(DP) :: dnlResPrimary, dnlResSecondary, dnlResTotal !, dnlResScale

    ! Output block for UCD output to GMV file
    type(t_ucdExport) :: rexport
    real(DP), dimension(:), pointer :: p_Ddata,p_Ddata2

    ! A counter variable
    integer :: i
    logical :: bError

    ! A timer object for time measuring
    type(t_timer) :: rtimer

    ! Ok, let's start.
    !
    ! We want to solve our Stokes problem on level...
    NLMIN = 3
    NLMAX = 6

    ! Viscosity parameter:
    dnu = 0.04_DP

    ! Thermal conductivity parameter:
    deta = 0.01_DP

    ! Gravity direction:
    dgrav1 =  0.0_DP
    dgrav2 = -1.0_DP

    ! Set the damping parameters for the systems.
    dnlDampPrimary = 1.0_DP
    dnlDampSecondary = 1.0_DP

    ! Allocate memory for all levels
    allocate(Rlevels(NLMIN:NLMAX))

    ! At first, read in the parametrisation of the boundary and save
    ! it to rboundary.
    call boundary_read_prm(rboundary, './pre/mit1.prm')

    ! Now read in the basic triangulation.
    call tria_readTriFile2D (Rlevels(NLMIN)%rtriangulation, &
                             './pre/mit1.tri', rboundary)

    ! Refine the mesh up to the minimum level
    call tria_quickRefine2LevelOrdering (NLMIN-1,&
        Rlevels(NLMIN)%rtriangulation,rboundary)

    ! And create information about adjacencies and everything one needs from
    ! a triangulation.
    call tria_initStandardMeshFromRaw (Rlevels(NLMIN)%rtriangulation,&
        rboundary)

    ! Now refine the grid for the fine levels.
    do i = NLMIN+1, NLMAX

      ! Refine the grid using the 2-Level-Ordering algorithm
      call tria_refine2LevelOrdering(Rlevels(i-1)%rtriangulation,&
          Rlevels(i)%rtriangulation,rboundary)

      ! Create a standard mesh
      call tria_initStandardMeshFromRaw(Rlevels(i)%rtriangulation,&
        rboundary)

    end do

    ! Now we can start to initialise the discretisation.
    do i = NLMIN, NLMAX

      ! Our primary discretisation has 3 components (Q1~/Q1~/Q0)
      call spdiscr_initBlockDiscr (Rlevels(i)%rdiscrPrimary, 3, &
                                   Rlevels(i)%rtriangulation, rboundary)

      ! For simplicity, we set up one discretisation structure for the
      ! velocity...
      call spdiscr_initDiscr_simple (Rlevels(i)%rdiscrPrimary%RspatialDiscr(1),&
          EL_EM30, CUB_G2X2, Rlevels(i)%rtriangulation, rboundary)

      ! ...and copy this structure also to the discretisation structure
      ! of the 2nd component (Y-velocity). This needs no additional memory,
      ! as both structures will share the same dynamic information afterwards.
      call spdiscr_duplicateDiscrSc (&
          Rlevels(i)%rdiscrPrimary%RspatialDiscr(1),&
          Rlevels(i)%rdiscrPrimary%RspatialDiscr(2))

      ! For the pressure (3rd component), we set up a separate discretisation
      ! structure, as this uses different finite elements for trial and test
      ! functions.
      call spdiscr_deriveSimpleDiscrSc (Rlevels(i)%rdiscrPrimary%RspatialDiscr(1), &
          EL_Q0, CUB_G2X2, Rlevels(i)%rdiscrPrimary%RspatialDiscr(3))

      ! Our secondary discretisation has 1 component (Q1)
      call spdiscr_initBlockDiscr (Rlevels(i)%rdiscrSecondary, 1, &
                                   Rlevels(i)%rtriangulation, rboundary)

      call spdiscr_initDiscr_simple (Rlevels(i)%rdiscrSecondary%RspatialDiscr(1),&
          EL_Q1, CUB_G2X2, Rlevels(i)%rtriangulation, rboundary)

    end do

    do i = NLMIN, NLMAX

      ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
      ! Initialise primary system for this level
      ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

      ! Initialise a block matrix based on the primary discretisation.
      call lsysbl_createMatBlockByDiscr (Rlevels(i)%rdiscrPrimary,&
                                         Rlevels(i)%rmatrixPrimary)

      ! Inform the matrix that we build a saddle-point problem.
      Rlevels(i)%rmatrixPrimary%imatrixSpec = LSYSBS_MSPEC_SADDLEPOINT

      ! Create the matrix structure of the A-matrices of the primary system.
      call bilf_createMatrixStructure (Rlevels(i)%rdiscrPrimary%RspatialDiscr(1),&
          LSYSSC_MATRIX9, Rlevels(i)%rmatrixA)

      ! Create the matrix structure of the B-matrices
      call bilf_createMatrixStructure (Rlevels(i)%rdiscrPrimary%RspatialDiscr(3),&
          LSYSSC_MATRIX9, Rlevels(i)%rmatrixB1,&
          Rlevels(i)%rdiscrPrimary%RspatialDiscr(1))

      ! Since B2 has the same structure as B1, we'll simply get a copy of it.
      call lsyssc_duplicateMatrix (Rlevels(i)%rmatrixB1, Rlevels(i)%rmatrixB2,&
                                   LSYSSC_DUP_SHARE, LSYSSC_DUP_REMOVE)

      ! Build the Laplace-Matrix of the primary system.
      call stdop_assembleLaplaceMatrix(Rlevels(i)%rmatrixA, .true., dnu)

      ! Build the B1-matrix entries
      call stdop_assembleSimpleMatrix(Rlevels(i)%rmatrixB1, DER_FUNC, DER_DERIV_X,&
                                      -1.0_DP, .true.)

      ! Build the B2-matrix entries
      call stdop_assembleSimpleMatrix(Rlevels(i)%rmatrixB2, DER_FUNC, DER_DERIV_Y,&
                                      -1.0_DP, .true.)

      ! Now let's set up our primary system block matrix.

      ! Copy the A-matrix into our primary system block matrix.
      call lsyssc_duplicateMatrix (Rlevels(i)%rmatrixA, &
          Rlevels(i)%rmatrixPrimary%RmatrixBlock(1,1),LSYSSC_DUP_SHARE,LSYSSC_DUP_COPY)
      call lsyssc_duplicateMatrix (Rlevels(i)%rmatrixPrimary%RmatrixBlock(1,1), &
          Rlevels(i)%rmatrixPrimary%RmatrixBlock(2,2),LSYSSC_DUP_SHARE,LSYSSC_DUP_SHARE)

      ! Copy the B1-/B2-matrices into our primary system block matrix.
      call lsyssc_duplicateMatrix (Rlevels(i)%rmatrixB1, &
          Rlevels(i)%rmatrixPrimary%RmatrixBlock(1,3),LSYSSC_DUP_SHARE,LSYSSC_DUP_COPY)
      call lsyssc_duplicateMatrix (Rlevels(i)%rmatrixB2, &
          Rlevels(i)%rmatrixPrimary%RmatrixBlock(2,3),LSYSSC_DUP_SHARE,LSYSSC_DUP_COPY)

      ! Furthermore, put B1^T and B2^T to the block matrix.
      call lsyssc_transposeMatrix (Rlevels(i)%rmatrixB1, &
          Rlevels(i)%rmatrixPrimary%RmatrixBlock(3,1),LSYSSC_TR_VIRTUAL)

      call lsyssc_transposeMatrix (Rlevels(i)%rmatrixB2, &
          Rlevels(i)%rmatrixPrimary%RmatrixBlock(3,2),LSYSSC_TR_VIRTUAL)

      ! Update the structural information of the block matrix.
      call lsysbl_updateMatStrucInfo (Rlevels(i)%rmatrixPrimary)

      ! Finally, create a solution vector based on the matrix.
      call lsysbl_createVecBlockIndMat (Rlevels(i)%rmatrixPrimary,&
                                        Rlevels(i)%rvecSolPrimary,.false.)

      ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
      ! Initialise secondary system for this level
      ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

      ! Initialise a block matrix based on the secondary discretisation.
      call lsysbl_createMatBlockByDiscr (Rlevels(i)%rdiscrSecondary,&
                                         Rlevels(i)%rmatrixSecondary)

      ! Create the matrix structure of the secondary system.
      call bilf_createMatrixStructure (Rlevels(i)%rdiscrSecondary%RspatialDiscr(1),&
          LSYSSC_MATRIX9, Rlevels(i)%rmatrixN)

      ! Build the Laplace-Matrix of the secondary system.
      call stdop_assembleLaplaceMatrix(Rlevels(i)%rmatrixN, .true., deta)

      ! Copy the N-matrix into our secondary system block matrix.
      call lsyssc_duplicateMatrix (Rlevels(i)%rmatrixN, &
          Rlevels(i)%rmatrixSecondary%RmatrixBlock(1,1),LSYSSC_DUP_SHARE,LSYSSC_DUP_COPY)

      ! Update the structural information of the block matrix.
      call lsysbl_updateMatStrucInfo (Rlevels(i)%rmatrixSecondary)

      ! Finally, create a solution vector based on the matrix.
      call lsysbl_createVecBlockIndMat (Rlevels(i)%rmatrixSecondary,&
                                        Rlevels(i)%rvecSolSecondary,.false.)

      ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
      ! Initialise coupling matrices for this level
      ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

      ! Create the matrix structure of the M-matrices
      call bilf_createMatrixStructure (Rlevels(i)%rdiscrSecondary%RspatialDiscr(1),&
          LSYSSC_MATRIX9, Rlevels(i)%rmatrixM,&
          Rlevels(i)%rdiscrPrimary%RspatialDiscr(1))

      ! Build a mass matrix.
      call stdop_assembleSimpleMatrix(Rlevels(i)%rmatrixM, DER_FUNC, DER_FUNC,&
                                      1.0_DP, .true.)

      ! That's it for now.

    end do

    ! Create the RHS vectors on the finest level.
    call lsysbl_createVecBlockIndMat (Rlevels(NLMAX)%rmatrixPrimary,&
                                      rrhsPrimary, .false.)
    call lsysbl_createVecBlockIndMat (Rlevels(NLMAX)%rmatrixPrimary,&
                                      rvecTmpPrimary, .false.)
    call lsysbl_createVecBlockIndMat (Rlevels(NLMAX)%rmatrixPrimary,&
                                      rvecDefPrimary, .false.)
    call lsysbl_createVecBlockIndMat (Rlevels(NLMAX)%rmatrixSecondary,&
                                      rrhsSecondary, .false.)
    call lsysbl_createVecBlockIndMat (Rlevels(NLMAX)%rmatrixSecondary,&
                                      rvecTmpSecondary, .false.)
    call lsysbl_createVecBlockIndMat (Rlevels(NLMAX)%rmatrixSecondary,&
                                      rvecDefSecondary, .false.)

    ! Clear both RHS vectors and the solution vectors on the finest level.
    call lsysbl_clearVector(rrhsPrimary)
    call lsysbl_clearVector(rrhsSecondary)
    call lsysbl_clearVector(Rlevels(NLMAX)%rvecSolPrimary)
    call lsysbl_clearVector(Rlevels(NLMAX)%rvecSolSecondary)

    do i = NLMIN, NLMAX

      ! Now implement the boundary conditions.
      call bcasm_initDiscreteBC(Rlevels(i)%rdiscreteBCprimary)
      call bcasm_initDiscreteBC(Rlevels(i)%rdiscreteBCsecondary)

      ! Get the bottom edge - Dirichlet(0) BCs for X- and Y-velocities, but
      ! Neumann for the temperature.
      call boundary_createRegion(rboundary,1,1,rboundaryRegion)

      ! X-/Y-velocity
      call bcasm_newDirichletBConRealBD (Rlevels(i)%rdiscrPrimary,1,&
          rboundaryRegion, Rlevels(i)%rdiscreteBCprimary, getBoundaryValuesPrimary)
      call bcasm_newDirichletBConRealBD (Rlevels(i)%rdiscrPrimary,2,&
          rboundaryRegion, Rlevels(i)%rdiscreteBCprimary, getBoundaryValuesPrimary)

      ! Get the left edge - Dirichlet(0) BC for the X- and Y-velocities and
      ! Dirichlet(-1/2) for the temperature.
      call boundary_createRegion(rboundary,1,2,rboundaryRegion)
      rboundaryRegion%iproperties = BDR_PROP_WITHSTART + BDR_PROP_WITHEND

      ! X-/Y-velocity
      call bcasm_newDirichletBConRealBD (Rlevels(i)%rdiscrPrimary,1,&
          rboundaryRegion, Rlevels(i)%rdiscreteBCprimary, getBoundaryValuesPrimary)
      call bcasm_newDirichletBConRealBD (Rlevels(i)%rdiscrPrimary,2,&
          rboundaryRegion, Rlevels(i)%rdiscreteBCprimary, getBoundaryValuesPrimary)

      ! temperature
      call bcasm_newDirichletBConRealBD (Rlevels(i)%rdiscrSecondary,1,&
          rboundaryRegion, Rlevels(i)%rdiscreteBCsecondary, getBoundaryValuesSecondary)

      ! Get the top edge - Dirichlet(0) BCs for X- and Y-velocities, but
      ! Neumann for the temperature.
      call boundary_createRegion(rboundary,1,3,rboundaryRegion)

      ! X-/Y-velocity
      call bcasm_newDirichletBConRealBD (Rlevels(i)%rdiscrPrimary,1,&
          rboundaryRegion, Rlevels(i)%rdiscreteBCprimary, getBoundaryValuesPrimary)
      call bcasm_newDirichletBConRealBD (Rlevels(i)%rdiscrPrimary,2,&
          rboundaryRegion, Rlevels(i)%rdiscreteBCprimary, getBoundaryValuesPrimary)

      ! Get the right edge - Dirichlet(0) BC for the X- and Y-velocities and
      ! Dirichlet(1/2) for the temperature.
      call boundary_createRegion(rboundary,1,4,rboundaryRegion)
      rboundaryRegion%iproperties = BDR_PROP_WITHSTART + BDR_PROP_WITHEND

      ! X-/Y-velocity
      call bcasm_newDirichletBConRealBD (Rlevels(i)%rdiscrPrimary,1,&
          rboundaryRegion, Rlevels(i)%rdiscreteBCprimary, getBoundaryValuesPrimary)
      call bcasm_newDirichletBConRealBD (Rlevels(i)%rdiscrPrimary,2,&
          rboundaryRegion, Rlevels(i)%rdiscreteBCprimary, getBoundaryValuesPrimary)

      ! temperature
      call bcasm_newDirichletBConRealBD (Rlevels(i)%rdiscrSecondary,1,&
          rboundaryRegion, Rlevels(i)%rdiscreteBCsecondary, getBoundaryValuesSecondary)

      ! Hang in the discrete BCs into the system matrices and solution vectors.
      Rlevels(i)%rmatrixPrimary%p_rdiscreteBC => Rlevels(i)%rdiscreteBCprimary
      Rlevels(i)%rmatrixSecondary%p_rdiscreteBC => Rlevels(i)%rdiscreteBCsecondary
      Rlevels(i)%rvecSolPrimary%p_rdiscreteBC => Rlevels(i)%rdiscreteBCprimary
      Rlevels(i)%rvecSolSecondary%p_rdiscreteBC => Rlevels(i)%rdiscreteBCsecondary

    end do

    ! Hang in the BCs into the RHS vectors.
    rrhsPrimary%p_rdiscreteBC => Rlevels(NLMAX)%rdiscreteBCprimary
    rvecDefPrimary%p_rdiscreteBC => Rlevels(NLMAX)%rdiscreteBCprimary
    rrhsSecondary%p_rdiscreteBC => Rlevels(NLMAX)%rdiscreteBCsecondary
    rvecDefSecondary%p_rdiscreteBC => Rlevels(NLMAX)%rdiscreteBCsecondary

    ! And filter the RHS and solution vectors on the finest level.
    call vecfil_discreteBCrhs (rrhsPrimary)
    call vecfil_discreteBCrhs (rrhsSecondary)
    call vecfil_discreteBCsol (Rlevels(NLMAX)%rvecSolPrimary)
    call vecfil_discreteBCsol (Rlevels(NLMAX)%rvecSolSecondary)

    ! Initialise a filter chain.
    RfilterChain(1)%ifilterType = FILTER_DISCBCDEFREAL

    ! Set up two inter-level-projection structures - one for the primary and
    ! one for the secondary system.
    call mlprj_initProjectionMat (rprojPrimary,Rlevels(NLMAX)%rmatrixPrimary)
    call mlprj_initProjectionMat (rprojSecondary,Rlevels(NLMAX)%rmatrixSecondary)

    ! Create two multi-grid solvers.
    p_RfilterChain => RfilterChain
    call linsol_initMultigrid (p_rsolverPrimary, p_RfilterChain)
    call linsol_initMultigrid (p_rsolverSecondary, p_RfilterChain)

    ! Set up a BiCGStab solver with VANKA preconditioner as coarse grid solver
    ! for the primary system.
    !CALL linsol_initVANKA (p_rpreconditioner,1.0_DP,LINSOL_VANKA_2DNAVST)
    !CALL linsol_initBiCGStab (p_rcoarseGridSolver,p_rpreconditioner,p_RfilterChain)
    call linsol_initUMFPACK4(p_rcoarseGridSolver)
    call linsol_addMultiGridLevel(p_rlevelInfo,p_rsolverPrimary,rprojPrimary,&
                                  null(), null(), p_rcoarseGridSolver)

    ! Set up a BiCGStab solver with SSOR preconditioner as coarse grid solver
    ! for the secondary system.
    !CALL linsol_initSSOR (p_rpreconditioner)
    !CALL linsol_initBiCGStab (p_rcoarseGridSolver,p_rpreconditioner,p_RfilterChain)
    call linsol_initUMFPACK4(p_rcoarseGridSolver)
    call linsol_addMultiGridLevel(p_rlevelInfo,p_rsolverSecondary,rprojSecondary,&
                                  null(), null(), p_rcoarseGridSolver)

    ! Now set up the other levels...
    do i = NLMIN+1, NLMAX

      ! Set up a VANKA smoother for the primary system.
      call linsol_initVANKA (p_rsmoother,1.0_DP,LINSOL_VANKA_2DNAVST)
      call linsol_convertToSmoother(p_rsmoother, 4, 1.0_DP)
      call linsol_addMultiGridLevel(p_rlevelInfo,p_rsolverPrimary,rprojPrimary,&
                                    p_rsmoother, p_rsmoother, null())

      ! Set up a SSOR smoother for the secondary system.
      call linsol_initSSOR (p_rsmoother)
      call linsol_convertToSmoother(p_rsmoother, 4, 1.0_DP)
      call linsol_addMultiGridLevel(p_rlevelInfo,p_rsolverSecondary,rprojSecondary,&
                                    p_rsmoother, p_rsmoother, null())

    end do

    ! Set the output level of the solvers to 2 for some output
    p_rsolverPrimary%ioutputLevel = -1
    p_rsolverSecondary%ioutputLevel = -1

    ! And set the correct tolerance parameters
    p_rsolverPrimary%iresNorm = 0
    p_rsolverPrimary%depsRel = 1E-5_DP
    p_rsolverPrimary%depsAbs = 0.0_DP !1E-9_DP
    p_rsolverPrimary%istoppingCriterion = LINSOL_STOP_ONEOF
    p_rsolverPrimary%nminIterations = 3
    p_rsolverPrimary%nmaxIterations = 5

    p_rsolverSecondary%iresNorm = 0
    p_rsolverSecondary%depsRel = 1E-5_DP
    p_rsolverSecondary%depsAbs = 0.0_DP !1E-9_DP
    p_rsolverSecondary%istoppingCriterion = LINSOL_STOP_ONEOF
    p_rsolverSecondary%nminIterations = 3
    p_rsolverSecondary%nmaxIterations = 5

    ! Attach the system matrices to the solvers.
    allocate(Rmatrices(NLMIN:NLMAX))
    ! Primary System
    do i = NLMIN, NLMAX
      call lsysbl_duplicateMatrix (Rlevels(i)%rmatrixPrimary,&
          Rmatrices(i),LSYSSC_DUP_SHARE,LSYSSC_DUP_SHARE)
    end do
    call linsol_setMatrices(p_RsolverPrimary,Rmatrices(NLMIN:NLMAX))
    do i=NLMIN,NLMAX
      call lsysbl_releaseMatrix (Rmatrices(i))
    end do

    ! Secondary System
    do i = NLMIN, NLMAX
      call lsysbl_duplicateMatrix (Rlevels(i)%rmatrixSecondary,&
          Rmatrices(i),LSYSSC_DUP_SHARE,LSYSSC_DUP_SHARE)
    end do
    call linsol_setMatrices(p_RsolverSecondary,Rmatrices(NLMIN:NLMAX))
    do i=NLMIN,NLMAX
      call lsysbl_releaseMatrix (Rmatrices(i))
    end do

    deallocate(Rmatrices)

    ! Initialise only the structure of the solvers - the data is initialised
    ! in the non-linear loop below.
    call linsol_initStructure (p_rsolverPrimary, ierror)
    if (ierror .ne. LINSOL_ERR_NOERROR) stop
    call linsol_initStructure (p_rsolverSecondary, ierror)
    if (ierror .ne. LINSOL_ERR_NOERROR) stop

    ! Now prepare the trilinear forms for the convection in the secondary system.
    ! u1 * d_x * T
    rtriform1%itermCount = 1
    rtriform1%ballCoeffConstant = .true.
    rtriform1%BconstantCoeff = .true.
    rtriform1%Dcoefficients = 1.0_DP
    rtriform1%Idescriptors(1,1) = DER_FUNC
    rtriform1%Idescriptors(2,1) = DER_DERIV_X
    rtriform1%Idescriptors(3,1) = DER_FUNC
    ! u2 * d_y * T
    rtriform2%itermCount = 1
    rtriform2%ballCoeffConstant = .true.
    rtriform2%BconstantCoeff = .true.
    rtriform2%Dcoefficients = 1.0_DP
    rtriform2%Idescriptors(1,1) = DER_FUNC
    rtriform2%Idescriptors(2,1) = DER_DERIV_Y
    rtriform2%Idescriptors(3,1) = DER_FUNC

    ! And prepare the streamline diffusion object for the convection in the
    ! primary system.
    rsd%dnu = dnu

    ! Start the timer
    call stat_clearTimer(rtimer)
    call stat_startTimer(rtimer)

    ! Now comes the fun part - the non-linear loop.
    do NLiter = 0, 50

      ! The very first thing is we need to update the non-linear parts of
      ! the system matrices, but to do this, we need the solution vectors
      ! on all levels - so let's restrict the solution vectors.
      do i = NLMAX, NLMIN+1, -1

        ! Restrict primary solution vector.
        call mlprj_performInterpolation(rprojPrimary, Rlevels(i-1)%rvecSolPrimary,&
                         Rlevels(i)%rvecSolPrimary, rvecTmpPrimary%RvectorBlock(1))

        ! Restrict secondary solution vector.
        call mlprj_performInterpolation(rprojSecondary, Rlevels(i-1)%rvecSolSecondary,&
                         Rlevels(i)%rvecSolSecondary, rvecTmpSecondary%RvectorBlock(1))

        ! And filter the restricted vectors.
        call vecfil_discreteBCsol (Rlevels(i-1)%rvecSolPrimary)
        call vecfil_discreteBCsol (Rlevels(i-1)%rvecSolSecondary)

      end do

      ! Now we can update the matrices - so go for it.
      do i = NLMIN, NLMAX

        ! Copy A-matrix into the X-velocity block of the primary system.
        call lsyssc_duplicateMatrix(Rlevels(i)%rmatrixA,&
            Rlevels(i)%rmatrixPrimary%RmatrixBlock(1,1),&
            LSYSSC_DUP_IGNORE,LSYSSC_DUP_COPY)

        ! Assemble the convection of the primary system using streamline diffusion.
        call conv_streamlineDiffusion2d(Rlevels(i)%rvecSolPrimary, &
            Rlevels(i)%rvecSolPrimary, 1.0_DP, 0.0_DP, rsd, CONV_MODMATRIX, &
            Rlevels(i)%rmatrixPrimary%RmatrixBlock(1,1))

        ! And filter the primary system matrix.
        call matfil_discreteBC (Rlevels(i)%rmatrixPrimary)

        ! The Y-velocity block is automatically updated, as it is just a shared
        ! copy of the X-velocity block.

        ! Copy N-matrix into the secondary system block matrix.
        call lsyssc_duplicateMatrix(Rlevels(i)%rmatrixN,&
            Rlevels(i)%rmatrixSecondary%RmatrixBlock(1,1),&
            LSYSSC_DUP_IGNORE,LSYSSC_DUP_COPY)

        ! Assemble the convection of the secondary system.
        ! Assemble u_1 * d_x T
        call trilf_buildMatrixScalar(rtriform1,.false.,&
             Rlevels(i)%rmatrixSecondary%RmatrixBlock(1,1), &
             Rlevels(i)%rvecSolPrimary%RvectorBlock(1))
        ! Assemble u_2 * d_y T
        call trilf_buildMatrixScalar(rtriform2,.false.,&
             Rlevels(i)%rmatrixSecondary%RmatrixBlock(1,1), &
             Rlevels(i)%rvecSolPrimary%RvectorBlock(2))

        ! And filter the secondary system matrix.
        call matfil_discreteBC (Rlevels(i)%rmatrixSecondary)

      end do

      ! Okay, the system matrices are up-to-date. Now we need to calculate
      ! the non-linear defect. But as we do not have a global system matrix,
      ! we need to calculate the defect of both the primary and secondary
      ! systems.
      ! So let's begin with the secondary system.
      call lsysbl_copyVector(rrhsSecondary, rvecDefSecondary)
      call lsysbl_blockMatVec(Rlevels(NLMAX)%rmatrixSecondary, &
          Rlevels(NLMAX)%rvecSolSecondary, rvecDefSecondary, -1.0_DP, 1.0_DP)

      ! Filter the defect vector.
      call vecfil_discreteBCdef(rvecDefSecondary)

      ! Calculate the secondary residual.
      dnlResSecondary = lsysbl_vectorNorm(rvecDefSecondary, LINALG_NORMEUCLID)

      ! Calculate the defect of the primary system.
      call lsysbl_copyVector(rrhsPrimary, rvecDefPrimary)
      call lsysbl_blockMatVec(Rlevels(NLMAX)%rmatrixPrimary, &
          Rlevels(NLMAX)%rvecSolPrimary, rvecDefPrimary, -1.0_DP, 1.0_DP)

      ! Don't forget the M-matrices which couple our primary system to
      ! the secondary system.
      call lsyssc_scalarMatVec(Rlevels(NLMAX)%rmatrixM,&
          Rlevels(NLMAX)%rvecSolSecondary%RvectorBlock(1),&
          rvecDefPrimary%RvectorBlock(1), -dgrav1, 1.0_DP)
      call lsyssc_scalarMatVec(Rlevels(NLMAX)%rmatrixM,&
          Rlevels(NLMAX)%rvecSolSecondary%RvectorBlock(1),&
          rvecDefPrimary%RvectorBlock(2), -dgrav2, 1.0_DP)

      ! Filter the defect vector
      call vecfil_discreteBCdef(rvecDefPrimary)

      ! Calculate the primary residual.
      dnlResPrimary = lsysbl_vectorNorm(rvecDefPrimary, LINALG_NORMEUCLID)

      ! And calculate the total residual.
      dnlResTotal = sqrt(dnlResPrimary**2 + dnlResSecondary**2)

      ! Print the residual
      !CALL output_separator(OU_SEP_MINUS)
      call output_line('NL-Iteration: ' // trim(sys_si(NLiter,2)) // &
                       '   |RES| = ' // trim(sys_sdEP(dnlResTotal,20,12)))
      !CALL output_separator(OU_SEP_MINUS)

      ! Is this precise enough?
      if(dnlResTotal .le. 1E-8_DP) exit

      ! Okay, we need to perform another non-linear iteration.

      !IF(dnlResSecondary .GT. 1E-10_DP) THEN
      if(.true.) then

        ! So let's first solve the secondary system to get a new temperature.

        ! Initialise the data of the secondary solver.
        call linsol_initData(p_rsolverSecondary, ierror)
        if(ierror .ne. LINSOL_ERR_NOERROR) stop

        ! And call the secondary solver
        !CALL output_line('Solving secondary system...')
        !CALL output_line('---------------------------')
        call linsol_precondDefect(p_rsolverSecondary, rvecDefSecondary)
        bError = (p_rsolverSecondary%iresult .ne. 0)

        ! Release the secondary solver data.
        call linsol_doneData(p_rsolverSecondary)

        ! Did the secondary solver break down?
        if(bError) then

          ! Print an error
          call output_separator(OU_SEP_STAR)
          call output_line('NL-Iteration: ERROR: secondary solver broke down')
          call output_separator(OU_SEP_STAR)

          ! Exit the non-linear loop
          exit

        end if

        ! Now update the secondary solution vector - which is the temperature.
        call lsysbl_vectorLinearComb(rvecDefSecondary, &
            Rlevels(NLMAX)%rvecSolSecondary, dnlDampSecondary, 1.0_DP)

        ! Now recalculate the defect of the primary system as we have updated
        ! the temperature.
        call lsysbl_copyVector(rrhsPrimary, rvecDefPrimary)
        call lsysbl_blockMatVec(Rlevels(NLMAX)%rmatrixPrimary, &
            Rlevels(NLMAX)%rvecSolPrimary, rvecDefPrimary, -1.0_DP, 1.0_DP)
        call lsyssc_scalarMatVec(Rlevels(NLMAX)%rmatrixM,&
            Rlevels(NLMAX)%rvecSolSecondary%RvectorBlock(1),&
            rvecDefPrimary%RvectorBlock(1), -dgrav1, 1.0_DP)
        call lsyssc_scalarMatVec(Rlevels(NLMAX)%rmatrixM,&
            Rlevels(NLMAX)%rvecSolSecondary%RvectorBlock(1),&
            rvecDefPrimary%RvectorBlock(2), -dgrav2, 1.0_DP)

        ! Filter the defect vector
        call vecfil_discreteBCdef(rvecDefPrimary)

        ! Calculate the primary residual.
        dnlResPrimary = lsysbl_vectorNorm(rvecDefPrimary, LINALG_NORMEUCLID)

      end if

      !IF (dnlResPrimary .GT. 1E-10_DP) THEN
      if(.true.) then

        ! Initialise the data of the primary solver.
        call linsol_initData(p_rsolverPrimary, ierror)
        if(ierror .ne. LINSOL_ERR_NOERROR) stop

        ! And call the primary solver
        !CALL output_line('Solving primary system...')
        !CALL output_line('---------------------------')
        call linsol_precondDefect(p_rsolverPrimary, rvecDefPrimary)
        bError = (p_rsolverPrimary%iresult .ne. 0)

        ! Release the primary solver data.
        call linsol_doneData(p_rsolverPrimary)

        ! Did the primary solver break down?
        if(bError) then

          ! Print an error
          call output_separator(OU_SEP_STAR)
          call output_line('NL-Iteration: ERROR: Primary solver broke down')
          call output_separator(OU_SEP_STAR)

          ! Exit the non-linear loop
          exit

        end if

        ! Now update the primary solution vector.
        call lsysbl_vectorLinearComb(rvecDefPrimary, &
            Rlevels(NLMAX)%rvecSolPrimary, dnlDampPrimary, 1.0_DP)

      end if

      ! That's it for this non-linear iteration...

    end do ! NLiter

    ! Stop the timer and print the elapsed time
    call stat_stopTimer(rtimer)
    call output_line('Total time for non-linear loop: ' // stat_sgetTime_byTimer(rtimer))

    ! We can now start the postprocessing.
    ! Start UCD export to GMV file:
    call ucd_startGMV (rexport,UCD_FLAG_STANDARD,&
        Rlevels(NLMAX)%rtriangulation,'gmv/u_mit_2.gmv')

    ! Project X- and Y-velocity to the vertices of the mesh.
    nullify(p_Ddata)
    nullify(p_Ddata2)
    call spdp_projectToVertices(Rlevels(NLMAX)%rvecSolPrimary%RvectorBlock(1), p_Ddata)
    call spdp_projectToVertices(Rlevels(NLMAX)%rvecSolPrimary%RvectorBlock(2), p_Ddata2)

    ! In case we use the VTK exporter, which supports vector output, we will
    ! pass the X- and Y-velocity at once to the ucd module.
    call ucd_addVarVertBasedVec(rexport,'velocity',p_Ddata,p_Ddata2)

    ! Release the temporary memory.
    deallocate(p_Ddata)
    deallocate(p_Ddata2)

    ! Write pressure
    call lsyssc_getbase_double (Rlevels(NLMAX)%rvecSolPrimary%RvectorBlock(3),p_Ddata)
    call ucd_addVariableElementBased (rexport,'pressure',UCD_VAR_STANDARD, p_Ddata)

    ! Write temperature
    call lsyssc_getbase_double (Rlevels(NLMAX)%rvecSolSecondary%RvectorBlock(1),p_Ddata)
    call ucd_addVariableVertexBased (rexport,'temperature',UCD_VAR_STANDARD, p_Ddata)

    ! Write the file to disc, that's it.
    call ucd_write (rexport)
    call ucd_release (rexport)

    ! We are finished - but not completely!
    ! Now, clean up so that all the memory is available again.
    !
    ! Release solver data and structure
   ! CALL linsol_doneData (p_rsolverNode)
    call linsol_doneStructure (p_rsolverSecondary)
    call linsol_doneStructure (p_rsolverPrimary)

    ! Release the interlevel projection structure
    call mlprj_doneProjection (rprojSecondary)
    call mlprj_doneProjection (rprojPrimary)

    ! Release the solver node and all subnodes attached to it (if at all):
    call linsol_releaseSolver (p_rsolverSecondary)
    call linsol_releaseSolver (p_rsolverPrimary)

    ! Release the block matrix/vectors
    call lsysbl_releaseVector (rvecDefPrimary)
    call lsysbl_releaseVector (rvecDefSecondary)
    call lsysbl_releaseVector (rvecTmpPrimary)
    call lsysbl_releaseVector (rvecTmpSecondary)
    call lsysbl_releaseVector (rrhsPrimary)
    call lsysbl_releaseVector (rrhsSecondary)
    do i = NLMAX, NLMIN, -1
      call lsysbl_releaseVector (Rlevels(i)%rvecSolSecondary)
      call lsysbl_releaseVector (Rlevels(i)%rvecSolPrimary)
      call lsysbl_releaseMatrix (Rlevels(i)%rmatrixSecondary)
      call lsysbl_releaseMatrix (Rlevels(i)%rmatrixPrimary)
    end do

    ! Release scalar matrices
    do i = NLMAX, NLMIN, -1
      call lsyssc_releaseMatrix (Rlevels(i)%rmatrixM)
      call lsyssc_releaseMatrix (Rlevels(i)%rmatrixN)
      call lsyssc_releaseMatrix (Rlevels(i)%rmatrixB2)
      call lsyssc_releaseMatrix (Rlevels(i)%rmatrixB1)
      call lsyssc_releaseMatrix (Rlevels(i)%rmatrixA)
    end do

    ! Release our discrete version of the boundary conditions
    do i = NLMAX, NLMIN, -1
      call bcasm_releaseDiscreteBC (Rlevels(i)%rdiscreteBCsecondary)
      call bcasm_releaseDiscreteBC (Rlevels(i)%rdiscreteBCprimary)
    end do

    ! Release the discretisation structure and all spatial discretisation
    ! structures in it.
    do i = NLMAX, NLMIN, -1
      call spdiscr_releaseBlockDiscr(Rlevels(i)%rdiscrSecondary)
      call spdiscr_releaseBlockDiscr(Rlevels(i)%rdiscrPrimary)
    end do

    ! Release the triangulation.
    do i = NLMAX, NLMIN, -1
      call tria_done (Rlevels(i)%rtriangulation)
    end do

    deallocate(Rlevels)

    ! Finally release the domain, that's it.
    call boundary_release (rboundary)

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
C$Procedure      SIGDGT ( Retain significant digits )

      SUBROUTINE SIGDGT ( IN, OUT )

C$ Abstract
C
C      Retain only the significant digits in a numeric string.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C      None.
C
C$ Keywords
C
C      CHARACTER,  PARSING
C
C$ Declarations

      CHARACTER*(*)    IN
      CHARACTER*(*)    OUT

C$ Brief_I/O
C
C      VARIABLE  I/O  DESCRIPTION
C      --------  ---  --------------------------------------------------
C      IN         I   Input numeric string.
C      OUT        O   Numeric string, with insignificant digits removed.
C
C$ Detailed_Input
C
C      IN          is a numeric string.
C
C$ Detailed_Output
C
C      OUT         is the same numeric string with insignificant
C                  zeros and spaces removed. The special case '.000...'
C                  becomes just '0'. OUT may overwrite IN. If the
C                  output string is too long, it is truncated on the
C                  right.
C
C$ Parameters
C
C     None.
C
C$ Particulars
C
C      There are only two interesting cases:
C
C         1) There is a decimal point and an exponent immediately
C            preceded by zero ('...0E', '...0D', '...0e', '...0d')
C            or by a space ('... E', '... D', '... e', '... d').
C
C         2) There is a decimal point and no exponent, and the last non-
C            blank character is a zero ('...0').
C
C      In each of these cases, go to the zero in question, and step
C      backwards until you find something other than a blank or a zero.
C
C      Finally, remove all leading spaces, and all occurrences of more
C      than one consecutive space within the string.
C
C$ Examples
C
C      The following examples illustrate the use of SIGDGT.
C
C      '0.123456000000D-04'        becomes     '0.123456D-04'
C      '  -9.2100000000000'                    '-9.21'
C      '       13'                             '13'
C      '    00013'                             '00013'
C      ' .314 159 265 300 000 e1'              '.314 159 265 3e1'
C      '   123    45     6'                    '123 45 6'
C      '  .000000000'                          '0'
C
C$ Restrictions
C
C      None.
C
C$ Exceptions
C
C      Error free.
C
C      If IN is a non-numeric string, the contents of OUT are
C      unpredictable.
C
C$ Files
C
C      None.
C
C$ Author_and_Institution
C
C      N.J. Bachman    (JPL)
C      H.A. Neilan     (JPL)
C      W.L. Taber      (JPL)
C      I.M. Underwood  (JPL)
C
C$ Literature_References
C
C      None.
C
C$ Version
C
C-     SPICELIB Version 1.0.1, 10-MAR-1992 (WLT)
C
C         Comment section for permuted index source lines was added
C         following the header.
C
C-     SPICELIB Version 1.0.0, 31-JAN-1990 (IMU)
C
C-&

C$ Index_Entries
C
C     retain significant digits
C
C-&

C$ Revisions
C
C-     Beta Version 1.3.0, 21-MAR-1989 (WLT)
C
C         Previous fix was unbelievably bad, very buggy.  This
C         has been fixed along with other bugs and non-standard
C         code has been removed.
C
C-     Beta Version 1.2.0, 28-FEB-1989 (WLT)
C
C         Reference to INSSUB replaced by SUFFIX
C
C-     Beta Version 1.1.1, 17-FEB-1989 (HAN)  (NJB)
C
C         Declaration of the unused function ISRCHC removed.
C
C
C-&

C
C     SPICELIB functions
C
      INTEGER               FRSTNB
      INTEGER               LASTNB
      INTEGER               CPOS

C
C     Local variables
C
      INTEGER               BEGIN
      INTEGER               END
      CHARACTER*1           LCHAR
      INTEGER               ZERO
      INTEGER               I
      INTEGER               J
      INTEGER               K
      INTEGER               L

C
C     Find the first and last non-blank characters in the string.
C
      BEGIN  = MAX ( 1, FRSTNB (IN) )
      END    = MAX ( 1, LASTNB (IN) )
      LCHAR  = ' '
C
C     Trivial case.
C
      IF ( BEGIN .EQ. END ) THEN

         OUT(1:1) = IN(BEGIN:BEGIN)

         IF ( LEN(OUT) .GT. 1 ) OUT(2:) = ' '

C
C     If there is no decimal point, all zeros are significant.
C
      ELSE IF ( INDEX (IN, '.') .EQ. 0 ) THEN

         L      = 1
         K      = BEGIN

         DO WHILE (       ( L .LE. LEN(OUT) )
     .              .AND. ( K .LE. END      ) )

            OUT(L:L) = IN(K:K)

C
C           Don't increment L if the last item copied was a space
C           (we don't want to copy extra spaces).
C
            IF ( IN(K:K) .NE. ' ' .OR. LCHAR .NE. ' ' )
     .      L        = L+1
            LCHAR    = IN(K:K)
            K        = K+1

         END DO

         IF ( L .LE. LEN(OUT) ) OUT(L:) = ' '

      ELSE

C
C        Is there is a decimal point and an exponent immediately
C        preceded by zero ('...0E', '...0D', '...0e', '...0d') or
C        by a space ('... E', '... D', '... e', '... d')?
C
         ZERO = INDEX ( IN, '0E' )

         IF ( ZERO .EQ. 0 )    ZERO = INDEX ( IN, '0D' )
         IF ( ZERO .EQ. 0 )    ZERO = INDEX ( IN, '0e' )
         IF ( ZERO .EQ. 0 )    ZERO = INDEX ( IN, '0d' )
         IF ( ZERO .EQ. 0 )    ZERO = INDEX ( IN, ' E' )
         IF ( ZERO .EQ. 0 )    ZERO = INDEX ( IN, ' D' )
         IF ( ZERO .EQ. 0 )    ZERO = INDEX ( IN, ' e' )
         IF ( ZERO .EQ. 0 )    ZERO = INDEX ( IN, ' d' )

C
C        Begin there, and move toward the front of the string until
C        something other than a blank or a zero is encountered. Then
C        remove the superfluous characters.
C
         IF ( ZERO .GT. 0 ) THEN

            J = ZERO + 1
            I = ZERO

            DO WHILE ( IN(I:I) .EQ. '0'  .OR. IN(I:I) .EQ. ' ' )
               I = I - 1
            END DO

            L = 1
            K = BEGIN

            DO WHILE (       ( L .LE. LEN(OUT) )
     .                 .AND. ( K .LE. I        ) )

               OUT(L:L) = IN(K:K)

C
C              Don't increment L if the last item copied was a space.
C
               IF ( IN(K:K) .NE. ' ' .OR. LCHAR .NE. ' ' )
     .         L        = L+1
               LCHAR    = IN(K:K)
               K        = K+1

            END DO

            K = J

            DO WHILE (       ( L .LE. LEN(OUT) )
     .                 .AND. ( K .LE. END      ) )

               OUT(L:L) = IN(K:K)

C
C              Increment L only if we don't have two consecutive
C              spaces.
C
               IF ( IN(K:K) .NE. ' ' .OR. LCHAR .NE. ' ' )
     .         L        = L+1
               LCHAR    = IN(K:K)
               K        = K+1

            END DO

            IF ( L .LE. LEN(OUT) ) OUT(L:) = ' '
C
C
C        Is there is a decimal point and no exponent, and is the last
C        non-blank character a zero ('...0')? Then truncate the string
C        after the last character that is neither a blank nor a zero.
C
         ELSE IF (       IN(END:END)            .EQ. '0'
     .             .AND. CPOS ( IN, 'EeDd', 1 ) .eq.  0 ) THEN

            I = END

            DO WHILE ( IN(I:I) .EQ. '0' .OR. IN(I:I) .EQ. ' ' )
               I = I - 1
            END DO

            L = 1
            K = BEGIN

            DO WHILE (       ( L .LE. LEN(OUT) )
     .                 .AND. ( K .LE. I        ) )

               OUT(L:L) = IN(K:K)

C
C              Increment L only if we don't have two consecutive
C              spaces.
C
               IF ( IN(K:K) .NE. ' ' .OR. LCHAR .NE. ' ' )
     .         L        = L+1
               LCHAR    = IN(K:K)
               K        = K+1

            END DO

            IF ( L .LE. LEN(OUT) ) OUT(L:) = ' '

         ELSE

            L = 1
            K = BEGIN

            DO WHILE (       ( L .LE. LEN(OUT) )
     .                 .AND. ( K .LE. END      ) )

               OUT(L:L) = IN(K:K)

C
C              Increment L only if we don't have two consecutive spaces.
C
               IF ( IN(K:K) .NE. ' ' .OR. LCHAR .NE. ' ' )
     .         L        = L+1
               LCHAR    = IN(K:K)
               K        = K+1

            END DO

            IF ( L .LE. LEN(OUT) ) OUT(L:) = ' '

         END IF

      END IF

C
C     Special case. The string '.0000....' reduces to '.' after the
C     zeros are removed.
C
      IF ( OUT .EQ. '.' ) THEN
         OUT = '0'
      END IF

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
! hgpt.f90
!
!Mateus, P.; Catalão, J.; Mendes, V.B.; Nico, G. An ERA5-Based Hourly Global Pressure and Temperature (HGPT) Model.
!Remote Sens. 2020, 12, 1098; https://doi.org/10.3390/rs12071098
!
! This routine determines the surface pressure (P), surface air temperature (T),
! weighed mean temperature (Tm), and zenith hydrostatic delay (ZHD) from binary coefficient files
! As available from:
! https://github.com/pjmateus/hgpt_model (release v1.0)
! press_grid.bin; temp_grid.bin; and tm_grid.bin
!
! It is admitted that the binary files with the coefficients are in the same directory as this script
! If you don't want it you can change it in the open statement
!
! The epoch can be an real*8 array of size 1, and in this case is the Modified Julian Date (MJD)
! or can be an integer array of size 6, with the Gregorian Calendar in the following format (year, month, day, hour, min, sec)
!
! All parameters are bilinear interpolated to the input ellipsoidal longitude and latitude
!
! Reference for HGPT:
! An ERA5-based hourly global temperature and pressure (HGTP) model (submitted to Remote Sensing, MDPI)
!
! INPUT:
!              dt : if size(dt)=1 => modified julian date
!                   if size(dt)=6 => year, month, day, hour, min, sec
!             ndt : size of dt
!              x0 : ellipsoidal longitude (degrees)
!              y0 : ellipsoidal latitude (degrees)
!              z0 : height (m)
!         z0_type : ‘orth’ for orthometric height or ‘elli’ for ellipsoidal height
!
! OUTPUT:
!               P : surface pressure valid at (x0, y0, z0), in hPa
!               T : surface air temperature valid at (x0, y0, z0), in Kelvins
!              Tm : weighed mean temperature valid at (x0, y0, z0), in Kelvins
!             ZHD : zenith hydrostatic delay, valid at (x0, y0, z0), in meters
!
! Set the following variables in the script from which you call hgpt.f90:
!--------------------------------------------------------------------------!
!       real :: x0, y0, z0, P, T, Tm, ZHD                                  !
!	real*8, dimension(6) :: dt           ! Gregorian Calendar          !
!                                                                          !
! Define variables values and call hgtp:                                   !
! 	call hgpt(dt, 6, x0, y0, z0, 'orth', P, T, Tm, ZHD)                !
! or                                                                       !
!       real :: x0, y0, z0, P, T, Tm, ZHD                                  !
!	real*8, dimension(1) :: dt           ! Modified Julian Date (MJD)  !
!                                                                          !
! Define variables values and call hgtp in the same way:                   !
! 	call hgpt(dt, 1, x0, y0, z0, 'orth', P, T, Tm, ZHD)                !
!--------------------------------------------------------------------------!
!Example:                                                                  !
!program call_hgpt                                                         !
!	implicit None                                                      !
!       real :: x0, y0, z0, P, T, Tm, ZHD                                  !
!	real*8, dimension(1) :: dt                                         !
!	y0 = 38.5519                                                       !
!	x0 = -9.0147                                                       !
!	z0 = 25                                                            !
!	dt(1) = 58119.5                                                    !
!	call hgpt(dt, size(dt), x0, y0, z0, 'orth', P, T, Tm, ZHD)         !
!	write(*,*) P, T, Tm, ZHD                                           !
!end program call_hgpt                                                     !
!--------------------------------------------------------------------------!
! written by Pedro Mateus (2020/01/15)
! Instituto Dom Luiz (IDL), Faculdade de Ciências, Universidade de Lisboa, 1749-016 Lisboa, Portugal
! pjmateus@fc.ul.pt
!
!
subroutine hgpt(dt, ndt, x0, y0, z0, z0_type, P, T, Tm, ZHD)
	implicit None
	real, intent(in) :: x0, y0, z0
	integer, intent(in) :: ndt
	real*8, intent(in) :: dt(ndt)
	character(len=4), intent(in) :: z0_type
	real, intent(out) :: T, P, Tm, ZHD
	real*8  :: mjd
	integer :: ios, i, j, year, month, day, hour, minute, sec
	integer, parameter :: row = 721, col = 1440
	real, parameter :: pi = 4*atan(1.0), p1 = 365.25, p2 = 182.6250, p3 = 91.3125, deg2rad = pi/180.0
	real, dimension(:) :: lon(col), lat(row)
	real, dimension(:,:) :: y_intercept(row, col), slope(row, col), a1(row, col), a2(row, col), a3(row, col), &
	                        orography(row, col), undu(row, col)
	integer*2, dimension(:,:) :: f1(row, col), f2(row, col), f3(row, col)
	real :: a, b, amp1, pha1, amp2, pha2, amp3, pha3, geo_height, N, H_orth, H_ellip

	! Input datetime format
	if ( ndt == 6)  then
		year = idint(dt(1)); month = idint(dt(2)); day = idint(dt(3))
		hour = idint(dt(4)); minute = idint(dt(5)); sec = idint(dt(6))
		call greg2modjul(year, month, day, hour, minute, sec, mjd)
	else if ( ndt == 1 ) then
		mjd = dt(1)
		call modjul2greg(mjd, year, month, day, hour, minute, sec)
	else
		stop 'Use 1) Modified Julian Date (MJD) or 2) Gregorian date (y,m,d,HH,MM,SS).'
	end if

	! Geographic latitude and longitude (equal to ERA5 grid)
	lon = (/((  -180+0.25*i), i=1, col)/)
	lat = (/((-90.25+0.25*i), i=1, row)/)

	! Locate i, j indices
	call locate_ij(row, lat, y0, i)
	call locate_ij(col, lon, x0, j)

	! Open and read the surface air temperature coefficients file
	open(10, file='temp_grid.bin', status='old', action='read', &
	     form='unformatted', access='stream', iostat=ios)
	if ( ios /= 0 )	stop 'Error opening file temp_grid.bin'

	! hour varies from 0 to 23 UTC
	call fseek(10, (row*col*26)*hour, 0)
	read(10) y_intercept
	read(10) slope
	read(10) a1
	read(10) f1
	read(10) a2
	read(10) f2
	read(10) a3
	read(10) f3
	close(10)

	! Bilinear interpolation
	call interpolate(lon, lat, y_intercept, x0, y0, i, j, a)
	call interpolate(lon, lat, slope, x0, y0, i, j, b)
	call interpolate(lon, lat, a1, x0, y0, i, j, amp1)
	call interpolate(lon, lat, float(f1)/10000.0, x0, y0, i, j, pha1)
	call interpolate(lon, lat, a2, x0, y0, i, j, amp2)
	call interpolate(lon, lat, float(f2)/10000.0, x0, y0, i, j, pha2)
	call interpolate(lon, lat, a3, x0, y0, i, j, amp3)
	call interpolate(lon, lat, float(f3)/10000.0, x0, y0, i, j, pha3)

	! Surface air temperature model
	T = a + b*(mjd - 51178) + amp1*cos(2*pi*(mjd - 51178)/p1+pha1) + &
             		          amp2*cos(2*pi*(mjd - 51178)/p2+pha2) + &
             	   	          amp3*cos(2*pi*(mjd - 51178)/p3+pha3)

	! Open and read the surface pressure coefficients file
	open(11, file='press_grid.bin', status='old', action='read', &
	     form='unformatted', access='stream', iostat=ios)
	if ( ios /= 0 ) stop 'Error opening file press_grid.bin'

	! hour varies from 0 to 23 UTC
	call fseek(11, (row*col*20)*hour, 0)
	read(11) y_intercept
	read(11) slope
	read(11) a1
	read(11) f1
	read(11) a2
	read(11) f2
	close(11)

	! Bilinear interpolation
	call interpolate(lon, lat, y_intercept, x0, y0, i, j, a)
	call interpolate(lon, lat, slope, x0, y0, i, j, b)
	call interpolate(lon, lat, a1, x0, y0, i, j, amp1)
	call interpolate(lon, lat, float(f1)/10000.0, x0, y0, i, j, pha1)
	call interpolate(lon, lat, a2, x0, y0, i, j, amp2)
	call interpolate(lon, lat, float(f2)/10000.0, x0, y0, i, j, pha2)

	! Surface pressure model
	P = a + b*(mjd - 51178) + amp1*cos(2*pi*(mjd - 51178)/p1+pha1) + &
                	          amp2*cos(2*pi*(mjd - 51178)/p2+pha2)

	! Open and read the weight mean temperature and undulation coefficients file
	open(12, file='tm_grid.bin', status='old', action='read', &
             form='unformatted', access='stream', iostat=ios)
	if ( ios /= 0 ) stop 'Error opening file tm_grid.bin'

	read(12) y_intercept
	read(12) slope
	read(12) orography
	read(12) undu
	close(12)

	! Bilinear interpolation
	call interpolate(lon, lat, y_intercept, x0, y0, i, j, a)
	call interpolate(lon, lat, slope, x0, y0, i, j, b)
	call interpolate(lon, lat, orography, x0, y0, i, j, geo_height)
	call interpolate(lon, lat, undu, x0, y0, i, j, N)

	! Zenith hydrostatic delay (ZHD), Saastamoinen model
	if ( z0_type == 'orth' ) then
		H_orth = z0
	else if ( z0_type == 'elli' ) then
		H_orth = z0 - N
	else
		stop 'Use 1) <<orth>> for Orthometric height or 2) <<elli>> for Ellipsoidal height (in m).'
	end if

	! Correction to P and T (see Guochanf Xu, GPS Theory, Algorithms and Applications, 2nd Edition, page 56)
	P = (P*100.0 * (1.0 - 0.0065/T * (H_orth - geo_height))**5.2559)/100.0
	T = T - 0.0065*(H_orth - geo_height)

	! Weight mean temperature (Tm) linear model
	Tm = a + b * T

	! ZHD using the Saastamoinen Model (see Saastamoinen, 1973)
	ZHD = (0.0022768 * P)/(1 - 0.0026*cos(2*deg2rad*y0)-0.00000028*H_orth)

	return
end subroutine hgpt

subroutine locate_ij(length, array, value, idx)
        ! Given an array and a value, returns the index of the element that
        ! is closest to, but less than, the given value.
        ! Uses a binary search algorithm.
        implicit none
        integer, intent(in) :: length
        real, dimension(length), intent(in) :: array
        real, intent(in) :: value
	integer, intent(out) :: idx
        integer :: left, middle, right

        left = 1
        right = length
        do
            if (left > right) then
                exit
            end if
            middle = nint((left+right) / 2.0)
            if ( abs(array(middle) - value) <= 1e-9) then
                idx = middle
                return
            else if (array(middle) > value) then
                right = middle - 1
            else
                left = middle + 1
            end if
        end do
        idx = right

	return
end subroutine locate_ij

subroutine interpolate(lon, lat, f, x0, y0, i, j, value)
        ! This function uses bilinear interpolation to estimate the value
        ! of a function f at point (x0,y0) using the point grid (i, j), (i+1,j), (i,j+1)
        ! and (i+1,j+1), see locate_ij subroutine
	! In the limits (borders) the nearest-neighbor interpolation is used
        ! f is assumed to be sampled on a regular grid, with the grid x values specified
        ! by lon array and the grid y values specified by lat array
        implicit none
	integer, parameter :: col = 1440, row = 721
        real, dimension(:), intent(in) :: lon(col), lat(row)
        real, dimension(row, col), intent(in) :: f
        integer, intent(in) :: i, j
	real, intent(in) :: x0, y0
	integer :: ix, jy
	real, intent(out) :: value
        real :: x1, x2, y1, y2, xy1, xy2

	ix = i
	jy = j
	if ( ix == 0 .or. ix >= row .or. jy == 0 .or. jy >= col ) then
		! Nearest-neighbor interpolation
		if ( ix == 0 ) ix = ix + 1
		if ( ix >= row ) ix = row
		if ( jy == 0 ) jy = jy + 1
		if ( jy >= col ) jy = col
		value = f(ix, jy)
	else
		! Bilinear interpolation
		x1 = lon( jy )
		x2 = lon( jy + 1 )
		y1 = lat( ix )
		y2 = lat( ix + 1 )
		xy1 = (x2 - x0)/(x2 - x1)*f(ix, jy) + (x0 - x1)/(x2 - x1)*f(ix, jy+1)
		xy2 = (x2 - x0)/(x2 - x1)*f(ix+1, jy) + (x0 - x1)/(x2 - x1)*f(ix+1, jy+1)
                value = (y2 - y0)/(y2 - y1)*xy1 + (y0 - y1)/(y2 - y1)*xy2
	end if

	return
end subroutine interpolate

subroutine greg2modjul(year, month, day, hour, minute, sec, mjd)
	! greg2modjul(..., mjd) returns the
	! Modified Julian Date (MJD) corresponding to the
	! Gregorian calendar date (year, month, day, hour, minute, and second).
	implicit none
        integer, intent(in) :: year, month, day, hour, minute, sec
	real*8, intent(out) :: mjd
        real*8              :: yyyy, mm, jd

        ! January & February
	if ( month <= 2 ) then
        	yyyy = dfloat(year) - 1_8
        	mm = dfloat(month) + 12_8
	else
        	yyyy = dfloat(year)
        	mm = dfloat(month)
        end if

	jd = floor( 365.25_8*(yyyy + 4716.0_8)) + floor( 30.6001_8*( mm + 1.0_8)) + 2.0_8 - &
	     floor( yyyy/100.0_8 ) + floor( floor( yyyy/100.0_8 )/4.0_8 ) + dfloat(day) - 1524.5_8 + &
             (dfloat(hour) + dfloat(minute)/60_8 + dfloat(sec)/3600_8)/24_8
        mjd = jd - 2400000.5_8

	return
end subroutine greg2modjul

subroutine modjul2greg(mjd, year, month, day, hour, minute, sec)
	! modjul2greg(mjd, ...) returns the
	! Gregorian calendar date (year, month, day, hour, minute, and second)
	! corresponding to the Modified Julian Date (MJD).
	implicit none
        real*8, intent(in)  :: mjd
	real*8              :: jd, ijd, fjd, a, b, c, d, e, m
	integer, intent(out):: year, month, day, hour, minute, sec

	jd = mjd + 2400000.5_8
	ijd = floor(jd + 0.5_8)
        fjd = jd - ijd + 0.5_8

        sec = 86400_8 * fjd
        hour = ifix( sec / 3600.0 )
        sec = sec - 3600_8 * hour
        minute = ifix( sec / 60.0 )
        sec = sec - 60_8 * minute

        a = ijd + 32044_8
        b = floor((4_8 * a + 3_8) / 146097_8)
        c = a - floor((b * 146097_8) / 4_8)
	d = floor((4_8 * c + 3_8) / 1461_8)
        e = c - floor((1461_8 * d) / 4_8)
        m = floor((5_8 * e + 2_8) / 153_8)

        day   = e - floor((153_8 * m + 2_8) / 5_8) + 1_8
        month = m + 3_8 - 12_8 * floor(m / 10_8)
        year  = b * 100_8 + d - 4800_8 + floor(m / 10_8)

	return
end subroutine modjul2greg

max_stars_repo_name: pjmateus/hgpt

max_issues_repo_name: pjmateus/hgpt

max_forks_repo_name: pjmateus/hgpt


text:
!===============================================================================
program OHC_adv
!===============================================================================
implicit none

!===============================================================================
!
!  This script calculates OHC budget from the original POP output files.
!  The area is a box around the  Weddel Sea deep water formation spot:
!  [78,60S]x[35W,80E]
!
!  components:
!  (1) advection:
!      a. meridional advection through section at specific latitude
!      b. zonal advection through section at certain longitude
!  (2) d(OHC)/dt
!  (3) surface heat fluxes
!
!  generated output
!
!  for curvilinear grid:
!  longitude (i) 0=250E, 1100=0E, 1500=40E
!  latitude  (j) 0=78.5S, 518=55S, 677=45S, 866=30S, 1181=0S
!  depth     (k) 1=5m, 9=96.9m, 10=112m, 14=198m, 17=318m, 35=4125, 42=5875m
!
!===============================================================================

!===============================================================================
!  variables
!===============================================================================

character*120 ::                                                               &
  LEC_folder,file_name,file_base,                                              &
  geometry1_file, geometry2_file
character*4   ::                                                               &
  year, yr
character*2   ::                                                               &
  month
integer                             ::                                         &
  imt,jmt,km,y,start_year,end_year,n,rec_length,i,j,k,m,t,                     &
  imin,imax,jmin,jmax,                                                         &
  nrec_UVEL,nrec_VVEL,nrec_TEMP
real, dimension(:),     allocatable ::                                         &
  dz,tdepth,area!,                                                              &

real, dimension(:,:),   allocatable ::                                         &
  DXT,DYT,TAREA,DXU,DYU,UAREA,geometry2!,                                       &

real, dimension(:,:,:), allocatable ::                                         &
  DZT,DZU,UVEL,VVEL,TEMP!,                                                      &

double precision, parameter         ::                                         &
c0 = 0., p125 = 0.125, p25 = 0.25, p5 = 0.5, c1 = 1.

imt        = 3600
jmt        = 2400
km         =   42

start_year =   75
end_year   =  326

imin       =  750 ! 35W
imax       = 1900 ! 80E
jmin       =    2 ! ca. 78S
jmax       =  600 ! ca. 60S

nrec_UVEL  =    1
nrec_VVEL  =   43

write (*,*) ''
write (*,*) '--- OHC advection along transects ---'
write (*,*) ''

!===============================================================================
!  FILES
!===============================================================================

LEC_folder      = '/home/dijkbio/andre/LEC/'
geometry1_file  = trim(LEC_folder)//'input/geometry1'
geometry2_file  = trim(LEC_folder)//'input/geometry2'

ADV_E_j_file = trim(LEC_folder)//'results/OHC_adv/OHC_adv_E_mint'

open(10,file=ADV_E_j_file,access='direct',form='unformatted',                  &
       recl=km,status='unknown')
!===============================================================================
!  GEOMETRY
!===============================================================================

! read 2D geometry fields
allocate( DXT(imt,jmt), DYT(imt,jmt), TAREA(imt,jmt), DZT(imt,jmt,km),         &
          DXU(imt,jmt), DYU(imt,jmt), UAREA(imt,jmt), DZU(imt,jmt,km) )
open(1,file=geometry1_file,access='direct',form='unformatted',recl=imt*jmt,    &
       status='old')
read(1,rec=1) DXT ! [m]
read(1,rec=2) DXU ! [m]
read(1,rec=3) DYT ! [m]read(1,rec=4) DYU ! [m]
read(1,rec=4) DYU ! [m]
read(1,rec=5) TAREA ! [m^2]
read(1,rec=6) UAREA ! [m^2]
do k=1,km
  read(1,rec=6+   k) DZT(:,:,k) ! [m]
  read(1,rec=6+km+k) DZU(:,:,k) ! [m]
enddo
close(1)

! read 1D geometry fields
allocate( geometry2(6,km),dz(km),tdepth(km),area(km) )
open(1,file=geometry2_file,access='direct',form='unformatted',recl=6,     &
       status='old')
do k=1,km
  read(1,rec=k) geometry2(:,k) ! all real: k,dz[m],tdepth[m],area[m^2],p[bar]
enddo
dz(:)     = geometry2(2,:)
tdepth(:) = geometry2(3,:)
area(:)   = geometry2(4,:)
close(1)

!===============================================================================
!  CALCULATIONS
!===============================================================================

write (*,*) ''
write (*,*) 'CALCULATIONS'

allocate( UVEL(imt,jmt,km),      VVEL(imt,jmt,km),     TEMP(imt,jmt,km),       &

! loop over years
do y=start_year,end_year

  if ( y<276 ) then
    file_base = '/projects/0/samoc/pop/tx0.1/output/run_henk_mixedbc/tavg/t.t0.1_42l_nccs01.'
    nrec_TEMP  =  127 ! for 12/14 GB files
  elseif ( y>=276 ) then
    file_base = '/projects/0/samoc/pop/tx0.1/output/run_henk_mixedbc_extravars_viebahn/tavg/t.t0.1_42l_nccs01.'
    nrec_TEMP  =  222 ! in ev case
  endif

  ! loop over months
  do m=1,12
    if ( m<12 ) then
      write(year,"(I0.4)") y
      write(month,"(I0.2)") m+1
    elseif ( m==12 ) then
      write(year,"(I0.4)") y+1
      write(month,"(I0.2)") 1
    endif

    ! file 27601 and 28101 missing!
    if ( y==275 .and. m==12 ) then
      write(year,"(I0.4)") y
      write(month,"(I0.2)") m
    elseif ( y==280 .and. m==12 ) then
      write(year,"(I0.4)") y
      write(month,"(I0.2)") m
    endif

    file_name = trim(file_base)//year//month
    open(1,file=file_name,access='direct',form='unformatted',                  &
           recl=imt*jmt,status='unknown')

    write(*,*) y, file_name

    do k=1,km
      read(1,rec=nrec_UVEL+k-1) UVEL(:,:,k)
      read(1,rec=nrec_VVEL+k-1) VVEL(:,:,k)
      read(1,rec=nrec_TEMP+k-1) TEMP(:,:,k)
    enddo

    !===========================================================================
    ! OHC advection
    !===========================================================================

    !===========================================================================
    ! d(OHC)/dt
    !===========================================================================

    !===========================================================================
    ! surface heat flux
    !===========================================================================

    !===========================================================================
    ! temporary storage
    ADV_N_year(:,:,m) = ADV_N(:,:)

    ! yearly averages
    if ( m==12 ) then
      ADV_E_j(:) = 0.0
      do t=1,12
        ADV_E_j(:) = ADV_E_j(:) + ADV_E_j_year(:,t)/12.0
      enddo !t
      !=========================================================================
      !  OUTPUT
      !=========================================================================
      write(10,rec=y-start_year+1) ADV_E_j(:)
    endif

    close (1)
  enddo !m
enddo !y

contains

subroutine OHC_merid_adv(imin,imax,lat,TEMP,VVEL,DZT,DXT,ADV,ADV_k,ADV_i)
! merdional advection at latitude lat
! northward flux through Northern face of T-cell(i,lat)
! ADV       in [J/m^2/s]
! ADV_k/i   in [J/m/s]

integer                                         :: i,k
integer,                            intent(in)  :: imin,imax,lat
integer, parameter                              :: imt=3600,jmt=2400,km=42
real,    dimension(imt,jmt),        intent(in)  :: DXT
real,    dimension(imt,jmt,km),     intent(in)  :: TEMP, VVEL, DZT
real,    dimension(imax-imin+1,km), intent(out) :: ADV
real,    dimension(km),             intent(out) :: ADV_i ! zonal-int
real,    dimension(imax-imin+1),    intent(out) :: ADV_k ! depth-int
real,    parameter                              :: rho0=4.1/3.996*1000, c=3996.0

do k=1,km
  do i=imin,imax
    ADV(i-imin+1,k) = (TEMP(i,lat+1,k)-TEMP(i,lat,k))                          &
                    * (VVEL(i-1,lat,k)+VVEL(i,lat,k))/2                        &
                    * rho0 * c
  enddo !i
enddo !k

ADV_k(:) = 0.0
do k=1,km
  ADV_k(:) = ADV_k(:) +  ADV(:,k)*DZT(imin:imax,lat,k)
enddo !k

ADV_i(:) = 0.0
do i=imin,imax
  ADV_i(:) = ADV_i(:) + ADV(i-imin+1,:)*DXT(i,lat)
enddo !i

end subroutine OHC_merid_adv

subroutine OHC_zonal_adv(jmin,jmax,lon,TEMP,UVEL,DZT,DYT,ADV,ADV_k,ADV_j)
! zonal advection at longitude imin/imax
! eastward flux through E face of T-cell(lon,i)
! ADV       in [J/m^2/s]
! ADV_k/j   in [J/m/s]

integer                                         :: j,k
integer,                            intent(in)  :: jmin,jmax,lon
integer, parameter                              :: imt=3600,jmt=2400,km=42
real,    dimension(imt,jmt),        intent(in)  :: DYT
real,    dimension(imt,jmt,km),     intent(in)  :: TEMP, UVEL, DZT
real,    dimension(jmax-jmin+1,km), intent(out) :: ADV
real,    dimension(km),             intent(out) :: ADV_j ! merid-int
real,    dimension(jmax-jmin+1),    intent(out) :: ADV_k ! depth-int
real,    parameter                              :: rho0=4.1/3.996*1000, c=3996.0

do k=1,km
  do j=jmin,jmax
    ADV(j-jmin+1,k) = (TEMP(lon+1,j,k)-TEMP(lon,j,k))                          &
                    * (UVEL(lon,j-1,k)+UVEL(lon,j,k))/2                        &
                    * rho0 * c
  enddo !j
enddo !k

ADV_k(:) = 0.0
do k=1,km
  ADV_k(:) = ADV_k(:) + ADV(:,k)*DZT(lon,jmin:jmax,k)
enddo !k

ADV_j(:) = 0.0
do j=jmin,jmax
  ADV_j(:) = ADV_j(:) + ADV(j-jmin+1,:)*DYT(lon,j)
enddo !i

end subroutine OHC_zonal_adv
!===============================================================================
end program OHC_adv
!===============================================================================

max_stars_repo_name: AJueling/LEC

max_issues_repo_name: AJueling/LEC

max_forks_repo_name: AJueling/LEC


text:
c     ****************************************************************
c     *                                                              *
c     *                      subroutine rktstf                       *
c     *                                                              *
c     *                       written by : bh                        *
c     *                                                              *
c     *                   last modified : 9/5/2017  rhd              *
c     *                                                              *
c     *     drive computation of tangent stiffness matrices for a    *
c     *     block of similar elements.                               *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine rktstf( props, iprops, lprops, glb_ek_blk, nrow_ek,
     &                   ispan, local_work )
      use main_data, only : matprp, lmtprp, asymmetric_assembly
      implicit none
      include 'param_def'
      include 'include_tan_ek'
c
c                 parameter declarations
c
      integer :: nrow_ek  ! set for symmetric or asymmetric
      integer :: ispan
      double precision ::  glb_ek_blk(nrow_ek,ispan)
      real    :: props(mxelpr,*)   !  all 3 are same. read-only here
      integer :: iprops(mxelpr,*)  !  1st col is for felem
      logical :: lprops(mxelpr,*)
c
c                 local
c
      integer :: span, felem, utsz, iter, type, order, ngp, nnode,
     &           ndof, totdof, mat_type, surf, local_iout, gpn,
     &           nlength, i, i_local

      logical :: geonl, local_debug, bbar, first, qbar_flag,
     &           compute_shape, cohes_mirror, dummy_logic,
     &           average, symmetric_assembly
c
      double precision ::
     &  xi, eta, zeta, beta_fact, eps_bbar, zero, one, dummy,
     &  temp_ref, d_temp, temp_np1
c
      data local_debug, dummy_logic / .false., .true. /
      data zero, one, dummy / 0.d0, 1.d0, 1.d0 /
c
c
      local_iout = local_work%iout
      if( local_debug ) write(local_iout,9100)
c
c              set common properties for all elements in block
c
      span      = local_work%span
      felem     = local_work%felem
      utsz      = local_work%utsz
      beta_fact = local_work%beta_fact
      eps_bbar  = local_work%eps_bbar
      first     = local_work%first
      iter      = local_work%iter
      qbar_flag = local_work%qbar_flag
c
      type      = local_work%elem_type
      order     = local_work%int_order
      ngp       = local_work%num_int_points
      nnode     = local_work%num_enodes
      ndof      = local_work%num_enode_dof
      geonl     = local_work%geo_non_flg
      totdof    = nnode * ndof
      bbar      = local_work%bbar_flg
      mat_type  = local_work%mat_type
      surf       = local_work%surface
      compute_shape =  local_work%is_cohes_elem .or.
     &                 local_work%is_axisymm_elem .or.
     &                 local_work%fgm_enode_props
      symmetric_assembly = .not. asymmetric_assembly
c
      if( symmetric_assembly ) then
            allocate( local_work%ek_symm(span,nrow_ek) )
            local_work%ek_symm = zero
            allocate( local_work%ek_full(1,1) )  ! for safe passing
      else
            allocate( local_work%ek_full(span,nrow_ek) )
            local_work%ek_full = zero
            allocate( local_work%ek_symm(1,1) )
      end if
c
c      if( asymmetric_assembly ) ek = zero
c
c               compute all the shape function derivatives,
c               inverse coordinate jacobians and determinants.
c               get weight function values for gauss integration.
c               for bbar elements (8-node bricks) we also
c               the additional volume terms. for geonl,
c               we use the nodal coordinates updated to
c               n+1 to compute jacobians, bbar terms,etc.
c
      if( bbar .and. type .eq. 2 ) then
         call rktstf_zero_vol( local_work%vol_block,
     &                  local_work%volume_block, span, mxvl )
      end if
c
c             for cohesive elements the element coordinates and
c             the global displacements are rotated to a coordinate
c             system in which the normal axis (Z rotated) is
c             perpendicular to the referernce surface of the
c             cohesive element
c
      if( local_work%is_cohes_elem )  then
            call cohes_rot_mat( span, felem, nnode, type,
     &                          local_work%ce,
     &                          local_work%cohes_rot_block )
          if ( geonl )
     &      call cohes_mirror_refsurf( span, mxvl, totdof, nnode,
     &                                 local_work%ce )
      end if
c
      do gpn = 1, ngp
         call getgpts( type, order, gpn, xi, eta, zeta,
     &                 local_work%weights(gpn) )
         call derivs( type, xi, eta, zeta, local_work%nxi(1,gpn),
     &                local_work%neta(1,gpn), local_work%nzeta(1,gpn) )
         if( compute_shape )
     &      call shapef( type, xi, eta, zeta, local_work%shape(1,gpn) )
         call jacob1( type, span, felem, gpn, local_work%jac_block,
     &                local_work%det_jac_block(1,gpn),
     &                local_work%gama_block(1,1,1,gpn),
     &                local_work%cohes_rot_block,
     &                local_work%nxi(1,gpn), local_work%neta(1,gpn),
     &                local_work%nzeta(1,gpn), local_work%ce,
     &                nnode )
         if ( bbar ) then
           call vol_terms( local_work%gama_block(1,1,1,gpn),
     &                     local_work%det_jac_block(1,gpn),
     &                     local_work%vol_block, local_work%nxi(1,gpn),
     &                     local_work%neta(1,gpn),
     &                     local_work%nzeta(1,gpn),
     &                     local_work%volume_block,
     &                     span, mxvl )
         end if
      end do
c
      if( bbar ) call vol_avg( local_work%vol_block,
     &                          local_work%volume_block, span, mxvl )
c
c               generate tangent stiffnesses for elements
c               in block, one gauss point at a time. zero
c               [b] matrices used in computations for block -
c               lower routines fill in only non-zero terms.
c
      nlength = mxvl * mxedof * nstr
      call rktstf_zero_vec( local_work%b_block,  nlength )
      call rktstf_zero_vec( local_work%bd_block, nlength )
c
      do gpn = 1, ngp
         call gptns1( local_work%cp, local_work%icp, gpn, props,
     &                iprops, glb_ek_blk, local_work )
       end do
c
      if( symmetric_assembly ) then
          call rktstf_do_transpose( local_work%ek_symm, glb_ek_blk,
     &                              span, nrow_ek )
      else
        call rktstf_do_transpose(  local_work%ek_full, glb_ek_blk,
     &                            span, nrow_ek )
      end if
      deallocate( local_work%ek_symm, local_work%ek_full )
c
c               modify element stiffness matrix by thickness
c               factor for plane strain analysis
c
      if( beta_fact .ne. one ) glb_ek_blk = beta_fact * glb_ek_blk
c
c               transform element stiffnesses to constraint
c               compatible coordinates if required. skip
c               whole block or specfic element to reduce work.
c
      if( local_work%trn_e_block ) then
        do i = 1, span
           i_local = i
           if( local_work%trn_e_flags(i) )
     &         call trnmtx( glb_ek_blk, local_work%cp,
     &                      local_work%trnmte, local_work%trne, ndof,
     &                      nnode, totdof, i_local, nrow_ek )
         end do
      end if
c
      return
c
 9100 format(8x,'>> entered rktstf...' )
 9200 format(12x,'>> upper-triangular stiffnesses for elements',
     &   ' in block...')
 9210 format(3x,i8,1000(6e15.6,/11x))
 9400 format(1x,//,'>> Upper triangluar matrix:')
 9405 format(1x,'   element #',i7)
 9410 format(1x,'j=',i3,2x,'row=',i3,2x,'col=',i3,2x,'ek=',f25.8)
 9300 format(1x,////,
     & '  ELEMENT TESTING, values in rktstf.f',//,
     & '       span = ',i7,/,
     & '      nnode = ',i7,/,
     & '     totdof = ',i7,/,
     & '  elem_type = ',i7,/,
     & '      felem = ',i7,/,
     & '       utsz = ',i7,/)
 9500 format(1x,'>> Fatal Error: rktstf. invalid material type..',
     &    /, 1x,'                job terminated' )
c
      end
c
c     ****************************************************************
c     *                                                              *
c     *                   subroutine rktstf_do_transpose             *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 09/19/2016rhd              *
c     *                                                              *
c     ****************************************************************
c

      subroutine rktstf_do_transpose( local_ek, ek, span, nrow_ek )
      implicit none
c
      integer :: span, nrow_ek
      integer :: i, j
      double precision ::
     & local_ek(span,nrow_ek), ek(nrow_ek,span)
c
c                 time 4.59 secs, tanstf time 12.64
c
c       ek = transpose( local_ek )
c       return
c
c
c           time 2.22 secs,  for 1830  tanstf.f 12.73395100
c             when nrow_ek is used, 2.78 secs

c
      select case ( nrow_ek ) ! careful. type could be integer constant
c
      case( 300 )
        do j = 1, span
!DIR$ VECTOR ALIGNED
              do i = 1, 300
                    ek(i,j) = local_ek(j,i)
              end do
        end do
c
      case( 465 )
        do j = 1, span
!DIR$ VECTOR ALIGNED
              do i = 1, 465
                    ek(i,j) = local_ek(j,i)
              end do
        end do
c
      case( 666 )
        do j = 1, span
!DIR$ VECTOR ALIGNED
              do i = 1, 666
                    ek(i,j) = local_ek(j,i)
              end do
        end do
c
      case( 1830 )
        do j = 1, span
!DIR$ VECTOR ALIGNED
              do i = 1, 1830
                    ek(i,j) = local_ek(j,i)
              end do
        end do
c
      case default
        do j = 1, span
!DIR$ VECTOR ALIGNED
              do i = 1, nrow_ek
                    ek(i,j) = local_ek(j,i)
              end do
        end do
      end select
c
c
c
c                  time 5.16 secs
c
c      do i = 1, 1830
c            do j = 1, span
c                  ek(i,j) = local_ek(j,i)
c            end do
c      end do

      return
      end
c
c     ****************************************************************
c     *                                                              *
c     *                   subroutine rktstf_zero_vec                 *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 6/8/2016 rhd               *
c     *                                                              *
c     ****************************************************************
c
      subroutine rktstf_zero_vec( vec, nterms )
      implicit none
c
      integer :: nterms
      double precision :: vec(nterms), zero
      data zero / 0.0d00 /
c
!DIR$ VECTOR ALIGNED
      vec = zero
c
      return
      end
c     ****************************************************************
c     *                                                              *
c     *                      subroutine zero_vol                     *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 09/5/2017 rhd              *
c     *                                   meant to be inlined        *
c     *                                                              *
c     ****************************************************************
c
      subroutine rktstf_zero_vol ( vol, volume, span, mxvl )
      implicit none
      integer :: span, mxvl
      double precision  :: vol(mxvl,8,3), volume(span)
      double precision, parameter :: zero = 0.0d0
      vol = zero
      volume = zero
      return
      end