text:
#include "rundeck_opts.h"

      MODULE TRACER_DIAG_COM
!@sum Tracer diagnostic arrays
!@+    Mostly tracer independent, but this may deped on applications
!@auth Jean Lerner
!ver   1.0
      USE MODEL_COM, only: im,jm,lm
      USE DAGCOM, only: npts !npts are conservation quantities
      USE TRACER_COM, only: ntm
#ifdef TRACERS_ON
     *     , ntsurfsrcmax, nt3Dsrcmax
#endif
      IMPLICIT NONE
      SAVE

C**** TAIJS  <<<< KTAIJS and IJTS_xx are Tracer-Dependent >>>>
C**** TAJLS  <<<< KTAJLS and JLS_xx are Tracer-Dependent >>>>

!@dbparam to_volume_MixRat: For printout of tracer concentration
!@+   to_volume_MixRat=1: printout is in Volume Mixing Ratio
!@+   to_volume_MixRat=0: printout is in Mass Mixing Ratio
#ifdef TRACERS_SPECIAL_Shindell
      INTEGER, DIMENSION(NTM) :: to_volume_MixRat=1
#else
      INTEGER, DIMENSION(NTM) :: to_volume_MixRat=0
#endif
#if (defined TRACERS_WATER) || (defined TRACERS_OCEAN)
!@dbparam to_per_mil For printout of tracer concentration in permil
      INTEGER, DIMENSION(NTM) :: to_per_mil = 0
#endif
#ifdef TRACERS_ON
!@var MMR_to_VMR: converts tracer mass mixing ratio to volume mr
      REAL*8, DIMENSION(NTM) :: MMR_to_VMR

C**** TAIJLN
!@var TAIJLN 3D tracer diagnostics (all tracers)
      real*8, dimension(im,jm,lm,ntm) :: taijln
!@var SNAME_IJT, UNITS_IJT: Names and units of lat-sigma tracer IJ diags
      character(len=30), dimension(lm,ntm) :: sname_ijt,units_ijt
!@var LNAME_IJT: descriptions of tracer IJ diagnostics
      character(len=80), dimension(lm,ntm) :: lname_ijt = 'unused'
!@var SCALE_IJT: printout scaling factor for tracer IJ diagnostics
      REAL*8, dimension(lm,ntm) :: scale_ijt
!@var IR_IJT: range index of IJ diagnosts
      integer, dimension(ntm) :: ir_ijt
!@var IA_IJT: idacc-number for tracer IJ diags (same for all tracers)
      integer ia_ijt
!@var IJTC_POWER: power of 10 used for tracer IJ concentration diags
      integer, dimension(ntm) :: ijtc_power
!@var IJTM_POWER: power of 10 used for tracer IJ mass unit diags
      integer, dimension(ntm) :: ijtm_power
!@var IJT_XX Names for TAIJLN diagnostics (Not currently used)

C**** TAIJN
!@param KTAIJ number of 3D diagnostics for each tracer
#ifdef TRACERS_WATER
      integer, parameter :: ktaij=11
#else
      integer, parameter :: ktaij=3
#endif
!@var IJT_XX names for taijn diagnostics
      integer tij_conc,tij_surf,tij_mass
#ifdef TRACERS_WATER
!@var IJT_XX names for water-based taijn diagnostics
      integer tij_rvr,tij_seaice,tij_prec,tij_evap,tij_grnd,tij_lk1
     *     ,tij_lk2,tij_soil
#endif
!@var TAIJN lat/lon tracer diagnostics (all tracers)
      real*8, dimension(im,jm,ktaij,ntm) :: taijn
!@var SCALE_TIJ: printout scaling factor for tracer IJK diagnostics
      REAL*8, dimension(ktaij,ntm) :: scale_tij
!@var SNAME_TIJ,UNITS_TIJ: Names and units of lat-sigma tracer diags
      character(len=30), dimension(ktaij,ntm) :: sname_tij,units_tij
!@var LNAME_TIJ: descriptions of tracer IJK diags
      character(len=80), dimension(ktaij,ntm) :: lname_tij = 'unused'

C**** TAIJS  <<<< KTAIJS and IJTS_xx are Tracer-Dependent >>>>
!@parm KTAIJS number of special lat/lon tracer diagnostics
#ifdef TRACERS_SPECIAL_Shindell
      integer, parameter :: ktaijs=26
#else
      integer, parameter :: ktaijs=25
#endif
!@var TAIJS  lat/lon special tracer diagnostics; sources, sinks, etc.
      REAL*8, DIMENSION(IM,JM,ktaijs) :: TAIJS
!@var ijts_source tracer independent array for TAIJS surface src. diags
      INTEGER ijts_source(ntsurfsrcmax,ntm)
!@var SNAME_IJTS, UNITS_IJTS: Names & units of lat-sigma tracer diags
      character(len=30), dimension(ktaijs) :: sname_ijts,units_ijts
!@var LNAME_IJTS: descriptions of tracer IJTS diags
      character(len=80), dimension(ktaijs) :: lname_ijts = 'unused'
!@var SCALE_IJTS: printout scaling factor for tracer IJTS diagnostics
      REAL*8, dimension(ktaijs) :: scale_ijts
!@var IA_IJTS: idacc-number for tracer source/sink IJ diags
      integer ia_ijts(ktaijs)
!@var ijts_power: power of 10 used for tracer IJ source/sink diags
      INTEGER, DIMENSION(ktaijs) :: ijts_power
!@var ijts_index: tracer index associated with a TAIJS diagnostic
      INTEGER, DIMENSION(ktaijs) :: ijts_index

C**** TAJLN
!@parm ktajl,ktajlx number of TAJL tracer diagnostics;
!@+          ktajlx includes composits
      INTEGER, PARAMETER :: ktajl=8, ktajlx=ktajl+2
!@var TAJLN  vertical tracer diagnostics (all tracers)
      REAL*8, DIMENSION(JM,LM,ktajlx,NTM) :: TAJLN
!@var jlnt_xx Names for TAJLN diagnostics
      INTEGER jlnt_conc,jlnt_nt_tot,jlnt_nt_mm,jlnt_vt_tot,
     &  jlnt_vt_mm,jlnt_mc,jlnt_turb,jlnt_lscond
!@var SNAME_JLN: Names of lat-sigma tracer JL diagnostics
      character(len=30), dimension(ktajlx,ntm) :: sname_jln
!@var LNAME_JLN,UNITS_JLN: descriptions/units of tracer JL diagnostics
      character(len=50), dimension(ktajlx,ntm) :: lname_jln = 'unused'
     *     ,units_jln
!@var SCALE_JLQ: printout scaling factor for tracer JL diagnostics
      REAL*8, dimension(ktajlx) :: scale_jlq
!@var IA_JLQ,JGRID_JLQ: idacc-numbers,gridtypes for tracer JL diags
      integer, dimension(ktajlx) :: ia_jlq,jgrid_jlq
!@var JLQ_POWER: power of 10 used for tracer JL diagnostics
!@+        It is associated with a specific physical process
      integer jlq_power(ktajlx)
      integer ktajlp !@var helps keeps track of
!@var scale_jln: Scale for jl maps
      REAL*8, DIMENSION(ntm) :: scale_jln

C**** TAJLS  <<<< KTAJLS and JLS_xx are Tracer-Dependent >>>>
!@parm ktajls number of source/sink TAJLS tracer diagnostics;
#ifdef TRACERS_SPECIAL_Shindell
      INTEGER, PARAMETER :: ktajls=58
#else
      INTEGER, PARAMETER :: ktajls=34
#endif
!@var TAJLS  JL special tracer diagnostics for sources, sinks, etc
      REAL*8, DIMENSION(JM,LM,ktajls) :: TAJLS
!@var jls_source tracer independent array for TAJLS surface src. diags
      INTEGER jls_source(ntsurfsrcmax,ntm)
!@var jls_3Dsource tracer independent array for TAJLS 3D source diags
      INTEGER jls_3Dsource(nt3Dsrcmax,ntm)
!@var jls_decay tracer independent array for radioactive sinks
      INTEGER, DIMENSION(NTM) :: jls_decay
!@var jls_grav tracer independent array for grav. settling sink
      INTEGER, DIMENSION(NTM) :: jls_grav
!@var jls_index: Number of source/sink tracer JL diagnostics/tracer
      integer, dimension(ktajls) :: jls_index
!@var SNAME_JLS: Names of lat-sigma tracer JL sources/sinks
      character(len=30), dimension(ktajls) :: sname_jls
!@var LNAME_JLS,UNITS_JLS: descriptions/units of tracer JLS diags
      character(len=50), dimension(ktajls) :: lname_jls = 'unused'
     *     ,units_jls
!@var SCALE_JLS: printout scaling factor for tracer JLS diagnostics
      REAL*8, dimension(ktajls) :: scale_jls
!@var IA_JLS,JGRID_JLS: idacc-numbers,gridtypes for tracer JL diags
      integer, dimension(ktajls) :: ia_jls,jgrid_jls
!@var JLS_POWER: power of 10 used for tracer JLS diagnostics
      integer jls_power(ktajls)
!@var jls_ltop: Top layer for this diagnostic
      integer jls_ltop(ktajls)

C**** TCONSRV
!@parm NTCONS Maximum Number of special tracer conservation points
      INTEGER, PARAMETER :: ntcons=20
!@parm KTCON total number of conservation diagnostics for tracers
      INTEGER, PARAMETER :: KTCON=npts+ntcons+2
!@var TCONSRV conservation diagnostics for tracers
      REAL*8, DIMENSION(JM,ktcon,ntm) :: TCONSRV
!@var SCALE_TCON scales for tracer conservation diagnostics
      REAL*8, DIMENSION(ktcon,ntm) :: SCALE_TCON
!@var TITLE_TCON titles for tracer conservation diagnostics
      CHARACTER*38, DIMENSION(ktcon,ntm) :: TITLE_TCON
!@var IA_TCON IDACC numbers for tracer conservation diagnostics
      INTEGER, DIMENSION(ktcon,ntm) ::  IA_TCON
!@var NSUM_TCON indices for summation of conservation diagnostics
      INTEGER, DIMENSION(ktcon,ntm) :: NSUM_TCON
!@var NOFMT indices for TCONSRV array
      INTEGER, DIMENSION(ktcon,ntm) :: NOFMT
!@var CONPTS names of special processes for tracer conservation diags
      CHARACTER*16, DIMENSION(ntcons) :: CONPTS
!@var kt_power_inst,kt_power_change: Exponents for tracer conservation
      INTEGER, DIMENSION(ntm) :: kt_power_inst,kt_power_change
!@var name_tconsrv,lname_tconsrv,units_tconsrv: for tracer conservation
      character(len=20), dimension(ktcon,ntm) ::
     &   name_tconsrv,units_tconsrv
      character(len=80), dimension(ktcon,ntm) :: lname_tconsrv
!@var SCALE_INST,SCALE_CHANGE: Scale factors for tracer conservation
      REAL*8, dimension(ntm) :: SCALE_INST,SCALE_CHANGE
!@var itcon_surf Index array for surface source/sink conservation diags
      INTEGER, DIMENSION(ntsurfsrcmax,NTM) :: itcon_surf
!@var itcon_3D Index array for 3D source/sink conservation diags
      INTEGER, DIMENSION(nt3Dsrcmax,NTM) :: itcon_3Dsrc
!@var itcon_decay Index array for decay conservation diags
      INTEGER, DIMENSION(NTM) :: itcon_decay
!@var itcon_grav Index array for gravitational settling conserv. diags
      INTEGER, DIMENSION(NTM) :: itcon_grav
!@var itcon_mc Index array for moist convection conserv. diags
      INTEGER, DIMENSION(NTM) :: itcon_mc
!@var itcon_grav Index array for large-scale condensation conserv. diags
      INTEGER, DIMENSION(NTM) :: itcon_ss

C----------------------------------------------------
!@param KTACC total number of tracer diagnostic words
!@var TACC: Contains all tracer diagnostic accumulations
      INTEGER, PARAMETER ::
     *  ktacc=IM*JM*LM*NTM + IM*JM*ktaij*NTM + IM*JM*ktaijs +
     *        JM*LM*ktajlx*NTM + JM*LM*ktajls + JM*NTM*ktcon
      COMMON /TACCUM/ TAIJLN,TAIJN,TAIJS,TAJLN,TAJLS,TCONSRV
      REAL*8, DIMENSION(KTACC) :: TACC
      EQUIVALENCE (TACC,TAIJLN)
C----------------------------------------------------
#endif
      END MODULE TRACER_DIAG_COM

#ifdef TRACERS_ON
      SUBROUTINE SET_TCON(QCON,NAME_CON,QSUM,INST_UNIT,SUM_UNIT
     *     ,INST_SC,CHNG_SC, itr,CONPTs)
!@sum  SET_TCON assigns conservation diagnostic array indices
!@auth Gavin Schmidt
!@ver  1.0
      USE CONSTANT, only: sday
      USE MODEL_COM, only: dtsrc,nfiltr
      USE DAGCOM, only: npts,ia_d5d,ia_d5s,ia_filt,ia_12hr,ia_src,conpt0
      USE TRACER_DIAG_COM, only: ktcon,title_tcon,scale_tcon,nsum_tcon
     *     ,nofmt,ia_tcon,name_tconsrv,lname_tconsrv,units_tconsrv
     *     ,ntcons
      IMPLICIT NONE
!@var QCON denotes at which points conservation diags are saved
      LOGICAL, INTENT(IN),DIMENSION(ktcon-1) :: QCON
!@var QSUM sets whether each diag is included in final sum
!@+   should be zero for diags set using DIAGTCB (i.e. using difference)
      LOGICAL, INTENT(IN),DIMENSION(ktcon-1) :: QSUM
      LOGICAL, DIMENSION(ktcon-1) :: QSUM_CON   ! local version
!@var INST_SC scale for instantaneous value
      REAL*8, INTENT(IN) :: INST_SC
!@var CHNG_SC scale for changes
      REAL*8, INTENT(IN) :: CHNG_SC
!@var NAME_CON name of conservation quantity
      CHARACTER*8, INTENT(IN) :: NAME_CON
!@var sname name of conservation quantity (no spaces)
      CHARACTER*8 :: sname
!@var INST_UNIT string for unit for instant. values
      CHARACTER*20, INTENT(IN) :: INST_UNIT
!@var SUM_UNIT string for unit for summed changes
      CHARACTER*20, INTENT(IN) :: SUM_UNIT
!@var ITR index for the tracer
      INTEGER, INTENT(IN) :: ITR
!@var CONPTS conservation diag points for special tracer diags
      CHARACTER*16, INTENT(IN), DIMENSION(ntcons) :: CONPTS
      CHARACTER*11 CHGSTR
      INTEGER NI,NM,NS,N,k

C**** remove spaces in NAME_CON for netcdf names
      sname=TRIM(NAME_CON)
      do k=1,len_trim(NAME_CON)
        if (sname(k:k).eq." ") sname(k:k)="_"
      end do
C****
      NI=1
      NOFMT(1,itr) = NI
      TITLE_TCON(NI,itr) =
     *     "0INSTANT "//TRIM(NAME_CON)//" "//TRIM(INST_UNIT)
      SCALE_TCON(NI,itr) = INST_SC
      name_tconsrv(NI,itr) ="inst_"//sname
      lname_tconsrv(NI,itr) = "INSTANT "//TRIM(NAME_CON)
      units_tconsrv(NI,itr) = INST_UNIT
      NSUM_TCON(NI,itr) = -1
      IA_TCON(NI,itr) = 12
      NM=NI
      DO N=2,ktcon-1
        IF (QCON(N)) THEN
          NM = NM + 1
          NOFMT(N,itr) = NM
          QSUM_CON(NM)=.FALSE.
          IF (QSUM(N)) QSUM_CON(NM)=.TRUE.
          CHGSTR=" CHANGE OF "
          if (n.le.npts+1) then
            TITLE_TCON(NM,itr) = CHGSTR//TRIM(NAME_CON)//" BY "//
     *         CONPT0(N-1)
            name_tconsrv(NM,itr) =
     *           "chg_"//trim(sname)//"_"//TRIM(CONPT0(N-1)(1:3))
          else
            IF (.not. QSUM(N)) CHGSTR="     DELTA "
            TITLE_TCON(NM,itr) = CHGSTR//TRIM(NAME_CON)//" BY "//
     *           CONPTs(N-npts-1)
            name_tconsrv(NM,itr) =
     *           "chg_"//trim(sname)//"_"//TRIM(CONPTs(N-npts-1)(1:3))
          end if
          lname_tconsrv(NM,itr) = TITLE_TCON(NM,itr)
          units_tconsrv(NM,itr) = SUM_UNIT
          SELECT CASE (N)
          CASE (2)
            SCALE_TCON(NM,itr) = CHNG_SC/DTSRC
            IA_TCON(NM,itr) = ia_d5d
          CASE (3,4,5,6,7,9,11,12)
            SCALE_TCON(NM,itr) = CHNG_SC/DTSRC
            IA_TCON(NM,itr) = ia_d5s
          CASE (8)
            SCALE_TCON(NM,itr) = CHNG_SC/(NFILTR*DTSRC)
            IA_TCON(NM,itr) = ia_filt
          CASE (10)
            SCALE_TCON(NM,itr) = CHNG_SC*2./SDAY
            IA_TCON(NM,itr) = ia_12hr
          CASE (13:)   ! special tracer sources
            SCALE_TCON(NM,itr) = CHNG_SC/DTSRC
            IA_TCON(NM,itr) = ia_src
          END SELECT
        ELSE
          NOFMT(N,itr) = 0
        END IF
      END DO
      NS=NM+1
      IF (NS.gt.KTCON) THEN
        WRITE(6,*) "KTCON not large enough for extra conserv diags",
     *       KTCON,NI,NM,NS,NAME_CON
        call stop_model(
     &       "Change KTCON in tracer diagnostic common block",255)
      END IF
      DO NM=NI+1,NS-1
        NSUM_TCON(NM,itr) = -1
        IF (QSUM_CON(NM)) NSUM_TCON(NM,itr) = NS
      END DO
      TITLE_TCON(NS,itr) = " SUM OF CHANGES "//TRIM(SUM_UNIT)
      name_Tconsrv(NS,itr) ="sum_chg_"//trim(sname)
      lname_Tconsrv(NS,itr) = " SUM OF CHANGES OF "//TRIM(NAME_CON)
      units_Tconsrv(NS,itr) = SUM_UNIT
      SCALE_TCON(NS,itr) = 1.
      IA_TCON(NS,itr) = 12
      NSUM_TCON(NS,itr) = 0
      RETURN
      END SUBROUTINE set_tcon

#endif


text:
      double precision function Li2(x)
      implicit none
C     returns Li_2(x) for real x, minf < x < 1
      double precision x
      double complex wgplg

      if (x .gt. 1d0) then
        write(6,*) 'x>1 in Li2 function, src/Lib/Li2.f'
        stop
      endif

      Li2 = dble(wgplg(1,1,x))
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
!===============================================================================
!> non-ideal MHD module
!===============================================================================
MODULE non_ideal_mod
  USE io_mod
  USE trace_mod
  USE io_unit_mod
  USE scalar_mod
  USE stagger_mod
  implicit none
  private
  type, public:: non_ideal_t
    real:: eta_Ohm, gamma_AD, rho_ions
    logical:: mhd_Ohm, mhd_AD, is_used
  contains
    procedure:: init
    procedure:: update
  end type
  type(non_ideal_t), public:: non_ideal
CONTAINS

!===============================================================================
!> Read parameters
!===============================================================================
SUBROUTINE init (self)
  class (non_ideal_t):: self
  integer:: iostat
  real:: gamma_AD=75., rho_ions=1.0, eta_Ohm=0.0
  logical:: mhd_AD=.false., mhd_Ohm=.false.
  namelist /non_ideal_params/mhd_AD, mhd_Ohm, gamma_AD, eta_Ohm, rho_ions
  !---------------------------------------------------------------------------
  rewind (io_unit%input)
  read (io_unit%input, non_ideal_params,iostat=iostat)
  if (io%master) write (*, non_ideal_params)
  self%mhd_Ohm  = mhd_Ohm
  self%eta_Ohm  = eta_Ohm
  self%mhd_AD   = mhd_AD
  self%gamma_AD = gamma_AD
  self%rho_ions = rho_ions
  self%is_used  = mhd_Ohm .or. mhd_AD
END SUBROUTINE init

!===============================================================================
!> Update heating Q, electric field E, and Courant velocity u_max
!===============================================================================
SUBROUTINE update (self, gn, ds, d, Jx, Jy, Jz, Bx, By, Bz, Q, E, u_max)
  class (non_ideal_t):: self
  integer:: gn(3)
  real(8):: ds
  real:: u_max, gamma_AD
  logical::mhd_AD, mhd_Ohm
  real, dimension(:,:,:,:):: E
  real, dimension(:,:,:):: d, Jx, Jy, Jz, Bx, By, Bz, Q
  real, dimension(:,:,:), allocatable:: Ud_x, Ud_y, Ud_z
  real, dimension(:,:,:), allocatable:: ambi
  integer, save:: itimer=0
  !---------------------------------------------------------------------------
  ! Below is the actual ambipolar drift parameter as given in equation 20 in
  ! Padoan, Zweibel and Nordlund (2000).
  ! To only write ambi, the definition in Padoan et al. is already multiplied
  ! by density**(-3/2) as it appears in the calculation of the drift velocity.
  ! The term d**1.5 is a cosnequence of the single fluid approximation and the
  ! relation n_i = K(n_n/10^5cm^-3)^k + K' (n_n/10^3 cm^-3)^-2
  ! K = 3*10^-3cm^-3, k=0.5 and K'=4.64*10^-4cm^-3.
  ! based on assuming that dust grains do NOT affect ionization balance significantly
  ! => second term falls off quickly => n_i~n_n^0.5
  ! In order to compute this parameter correctly, we need to know TEMPERATURE,
  ! RESOLUTION (not directly applicable here because we use AMR), BOX LENGTH,
  ! AVERAGE DENSITY, COSMIC-RAY IONISATION (K_i parameter:
  ! K_i = n_i*n_H**(-0.5) = x_i*n_H**(-0.5) ) and AVERAGE VOLUMETRIC FLOW RATE av_sigv:
  !---------------------------------------------------------------------------
  ! ambi = 0.28 * (0.1*T)**0.5 * (N/128d0) * (L/10d0)**(-1) * (av_n/200d0)**(-0.5)
  !      * (K_i/1e-5)**(-1) * (av_sigv/2e-9)**(-1) * d**(-1.5)
  !
  ! Alternative definition in Duffin and Pudritz (2008) assumes constant parameter
  ! for gamma_AD
  ! gamma_AD = av_sigv / (m_i + m_n) = 3.28e13 g^-1cm^3s^-1
  ! ambi = 1.4 / (gamma_AD*d**1.5)
  ! rho_i = 1.
  !---------------------------------------------------------------------------
  if (.not.self%is_used) return
  call trace%begin ('non_ideal_t%update', itimer=itimer)
  if (self%mhd_AD) then
    call allocate_scalars_a (gn, ambi, Ud_x, Ud_y, Ud_z)
    ambi = 1./(self%gamma_AD*d*self%rho_ions)
    Ud_x = xdn(ambi) * (zup(Jy)*xdn(zup(Bz)) - yup(Jz)*xdn(yup(By)))
    Ud_y = ydn(ambi) * (xup(Jz)*ydn(xup(Bx)) - zup(Jx)*ydn(zup(Bz)))
    Ud_z = zdn(ambi) * (yup(Jx)*zdn(yup(By)) - xup(Jy)*zdn(xup(Bx)))
    E(:,:,:,1) = E(:,:,:,1) - (Ud_y*Bz - Ud_z*By)
    E(:,:,:,2) = E(:,:,:,2) - (Ud_z*Bx - Ud_x*Bz)
    E(:,:,:,3) = E(:,:,:,3) - (Ud_x*By - Ud_y*Bx)
    Q = Q + self%rho_ions*d*self%gamma_AD*(xup(Ud_x)**2+yup(Ud_y)**2+zup(Ud_z)**2)
    call deallocate_scalars_a (ambi, Ud_x, Ud_y, Ud_z)
    u_max = max(u_max,maxval((xup(Bx)**2+yup(By)**2+zup(Bz)**2)/(self%gamma_AD*self%rho_ions*d*ds)))
  end if
  if (self%mhd_Ohm) then
    E(:,:,:,1) = E(:,:,:,1) + self%eta_Ohm*Jx
    E(:,:,:,2) = E(:,:,:,2) + self%eta_Ohm*Jy
    E(:,:,:,3) = E(:,:,:,3) + self%eta_Ohm*Jz
    Q = Q + yup(zup(self%eta_Ohm*Jx**2)) + zup(xup(self%eta_Ohm*Jy**2)) + xup(yup(self%eta_Ohm*Jz**2))
    u_max = max(u_max,self%eta_Ohm/ds)
  end if
  call trace%end (itimer)
END SUBROUTINE update

END MODULE non_ideal_mod

max_stars_repo_name: applejwjcat/dispatch

max_issues_repo_name: applejwjcat/dispatch

max_forks_repo_name: applejwjcat/dispatch


text:
C Copyright(C) 2009-2017 National Technology & Engineering Solutions of
C Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C     * Neither the name of NTESS nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C=======================================================================
      SUBROUTINE QALAB (DBORD0, DVIEW0, CHLSIZ, DOQA, DOAXIS, DOLAB,
     &   CAPTN, TITLE, CREATE, MODIFY, DRAW, DLEGND, BLKCOL,  *)
C=======================================================================

C   --*** QALAB *** (BLOT) Draw standard label (PLT)
C   --   Written by Amy Gilkey - revised 02/02/88
C   --
C   --QALAB starts a "standard" plot by:
C   --   o starting a new plot,
C   --   o outlining the display area (if DOQA true),
C   --   o displaying the database title and creation, modification and
C   --        plot information (if DOQA true),
C   --   o displaying the plot caption, and
C   --   o returning the legend coordinates.
C   --
C   --Parameters:
C   --   DBORD0 - IN - the plot boundary, including legend
C   --      (left, right, bottom, top)
C   --   DVIEW0 - IN - the plot view boundary (left, right, bottom, top)
C   --   CHLSIZ - IN - the size of a character line
C   --   DOQA - IN - true iff QA information is to be included in label
C   --   DOAXIS - IN - true iff axis is to be numbered (not in this routine)
C   --   DOLAB - IN - true iff axis is to be labeled (not in this routine)
C   --   CAPTN - IN - the three-line plot caption
C   --   TITLE - IN - the database title
C   --   CREATE - IN - the database creation code name, version, date, time
C   --   MODIFY - IN - the database modification code name, version, date, time
C   --   DRAW - IN - the database plot code name, version, date, time
C   --   DLEGND - OUT - the location of the legend (device units)
C   --      (left, right, bottom, top)
C   --   BLKCOL - IN/OUT - the user selected colors of the element blocks.
C   --                    BLKCOL(0) = 1 if the user defined material
C   --                                colors should be used in mesh plots.
C   --                              = -1 if program selected colors should
C   --                                be used.
C   --                    BLKCOL(i) = the user selected color of element
C   --                               block i:
C   --                                  -2 - no color selected by user.
C   --                                  -1 - black
C   --                                   0 - white
C   --                                   1 - red
C   --                                   2 - green
C   --                                   3 - yellow
C   --                                   4 - blue
C   --                                   5 - cyan
C   --                                   6 - magenta
C   --   * - the return statement if the cancel function is active

C   --Routines Called:
C   --   GRABRT - (GRPLIB) Check for plot set abort
C   --   GRBOX  - (GRPLIB) Draw box
C   --   GRPBEG - (PLTLIB) Begin a new plot
C   --   GRTEXT - (GRPLIB) Display a software/hardware string
C   --   GRYCEN - (GRPLIB) Find center lines of text area
C   --   LENSTR - (STRLIB) Find string length

      PARAMETER (KLFT=1, KRGT=2, KBOT=3, KTOP=4)

      include 'params.blk'
      include 'dbnums.blk'

      REAL DBORD0(KTOP), DVIEW0(KTOP)
      REAL CHLSIZ
      LOGICAL DOQA, DOAXIS, DOLAB
      CHARACTER*80 CAPTN(3)
      CHARACTER*80 TITLE
      CHARACTER*(MXSTLN) CREATE(4), MODIFY(4), DRAW(4)
      REAL DLEGND(KTOP)
      CHARACTER*256 CDUM1(5)
      CHARACTER*8 CDUM2(10)
      INTEGER IDUM2, IDUM3(10)
      INTEGER BLKCOL(0:NELBLK)

      LOGICAL GRABRT

C   --Set background color

      CALL SETBCK (2, CDUM1, IDUM2, IDUM3, CDUM2, *100)
      GOTO 110
  100 CONTINUE
      WRITE (*, *) 'Problem setting background color'

  110 CONTINUE

C   --Begin a new plot

      CALL GRPBEG

C   --Set foreground color

      CALL UGRCOL (0, BLKCOL)

C   --Set up layout

      DLEGND(KLFT) = DVIEW0(KRGT) + CHLSIZ
      DLEGND(KRGT) = DBORD0(KRGT)
      DLEGND(KTOP) = DVIEW0(KTOP)
      DLEGND(KBOT) = DVIEW0(KBOT)

      IF (DOQA) THEN

C      --Outline the display area

         IF (GRABRT ()) RETURN 1
         CALL GRBOX ('L',
     &      DBORD0(KLFT), DBORD0(KRGT), DBORD0(KBOT), DBORD0(KTOP))

C *** Database title and dates ***

C      --Title

         IF (GRABRT ()) RETURN 1
         DXTITL = 0.5 * (DBORD0(KRGT) - DBORD0(KLFT))
         DYTITL = DBORD0(KTOP) -
     &      0.5 * (DBORD0(KTOP) - DVIEW0(KTOP) + CHLSIZ)

C      --Left-justify title
         DO 120 N1 = 1, 79
            IF (TITLE(N1:N1) .NE. ' ') GO TO 130
  120    CONTINUE
  130    CONTINUE

         CALL GRTEXC (DXTITL, DYTITL, TITLE(N1:))

C      --Creator
         call qaleg(dlegnd, chlsiz, create, 1.0, 'Cre')

C      --Modifier
         call qaleg(dlegnd, chlsiz, modify, 1.5, 'Mod')

C      --User
         call qaleg(dlegnd, chlsiz, draw,   1.5, 'Drw')

      end if
C *** Plot caption ***

      DO 140 IEND = 3, 1, -1
         IF (CAPTN(IEND) .NE. ' ') GOTO 150
  140 CONTINUE
  150 CONTINUE
      IF (IEND .GT. 0) THEN
         DXCAPT = DVIEW0(KLFT) + 0.5 * (DVIEW0(KRGT) - DVIEW0(KLFT))
         DTOP = DVIEW0(KBOT) - CHLSIZ
         IF (DOAXIS) THEN
            DTOP = DTOP - 3.0*CHLSIZ
         ELSE IF (DOLAB) THEN
            DTOP = DTOP - 1.5*CHLSIZ
         END IF
         DBOT = DBORD0(KBOT)
         CALL GRYCEN (CHLSIZ, DTOP, DBOT, IEND, IDUM)
         DO 160 I = 1, IEND
            IF (GRABRT ()) RETURN 1
            CALL GRTEXC (DXCAPT, DTOP, CAPTN(I))
            DTOP = DTOP - CHLSIZ
  160    CONTINUE
      END IF

C   --Flush buffer, so label is complete at this point
      CALL PLTFLU

      RETURN
      END

      subroutine qaleg(dlegnd, chlsiz, entry, factor, label)
      include 'params.blk'
      PARAMETER (KLFT=1, KRGT=2, KBOT=3, KTOP=4)
      REAL DLEGND(KTOP)
      REAL CHLSIZ
      CHARACTER*(MXSTLN) entry(4)
      character*3 label

      DLEGND(KTOP) = DLEGND(KTOP) - factor * CHLSIZ
C ... Limit creator, modified by, and drawn by strings to 20 characters
C     Larger strings cause problems with metafiles
      L1 = MIN(20, LENSTR(ENTRY(1)))
      CALL GRTEXT (DLEGND(KLFT), DLEGND(KTOP),
     &     label//': '//ENTRY(1)(:L1))
      DLEGND(KTOP) = DLEGND(KTOP) - CHLSIZ

      l3 = lenstr(entry(3))
      l4 = lenstr(entry(4))
      CALL GRTEXT (DLEGND(KLFT), DLEGND(KTOP),
     &     '  '//ENTRY(3)(:l3) //
     *     '  '//ENTRY(4)(:l4))

      return
      end

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


program main
  use parent, only : child_function,grandchild_subroutine
  use parent, only : sibling_function,GreatGrandChild_subroutine
  implicit none
  if (child_function()) call grandchild_subroutine
  if (sibling_function()) call GreatGrandChild_subroutine
end program


module set_velocity_module

  use ml_layout_module
  use bl_constants_module

  implicit none

  private

  public :: set_velocity_on_level, set_velocity

contains

  subroutine set_velocity_on_level(velocity,dx,time)

    type(multifab) , intent(inout) :: velocity(:)
    real(kind=dp_t), intent(in   ) :: dx,time

    ! local
    integer :: i,dm,ng
    integer :: lo(velocity(1)%dim), hi(velocity(1)%dim)

    real(kind=dp_t), pointer :: dp1(:,:,:,:)
    real(kind=dp_t), pointer :: dp2(:,:,:,:)
    real(kind=dp_t), pointer :: dp3(:,:,:,:)

    do i=1,nfabs(velocity(1))
       dp1 => dataptr(velocity(1),i)
       dp2 => dataptr(velocity(2),i)
       lo = lwb(get_box(velocity(1),i))
       hi = upb(get_box(velocity(1),i))
       select case(dm)
       case (2)
          call set_velocity_2d(dp1(:,:,1,1), dp2(:,:,1,1), ng, &
                                   lo, hi, dx, time)
       case (3)
          dp3 => dataptr(velocity(3),i)
          call set_velocity_3d(dp1(:,:,:,1), dp2(:,:,:,1), dp3(:,:,:,1), ng, &
                                   lo, hi, dx, time)
       end select
    end do

    do i=1,dm
       call multifab_fill_boundary(velocity(i))
    end do

  end subroutine set_velocity_on_level

  subroutine set_velocity(mla,velocity,dx,time)

    type(ml_layout), intent(in   ) :: mla
    type(multifab) , intent(inout) :: velocity(:,:)
    real(kind=dp_t), intent(in   ) :: dx(:),time

    ! local variables
    integer :: lo(mla%dim), hi(mla%dim)
    integer :: nlevs, dm, ng, i, n

    real(kind=dp_t), pointer :: dp1(:,:,:,:)
    real(kind=dp_t), pointer :: dp2(:,:,:,:)
    real(kind=dp_t), pointer :: dp3(:,:,:,:)

    ng = velocity(1,1)%ng
    dm = mla%dim
    nlevs = mla%nlevel

    do n=1,nlevs
       do i=1,nfabs(velocity(n,1))
          dp1 => dataptr(velocity(n,1),i)
          dp2 => dataptr(velocity(n,2),i)
          lo = lwb(get_box(velocity(n,1),i))
          hi = upb(get_box(velocity(n,1),i))
          select case(dm)
          case (2)
             call set_velocity_2d(dp1(:,:,1,1), dp2(:,:,1,1), ng, &
                                      lo, hi, dx(n), time)
          case (3)
             dp3 => dataptr(velocity(n,3),i)
             call set_velocity_3d(dp1(:,:,:,1), dp2(:,:,:,1), dp3(:,:,:,1), ng, &
                                      lo, hi, dx(n), time)
          end select
       end do
    end do

    do n=1,nlevs
       do i=1,dm
          call multifab_fill_boundary(velocity(n,i))
       end do
    end do

  end subroutine set_velocity

  subroutine set_velocity_2d(velx, vely, ng, lo, hi, dx, time)

    integer          :: lo(2), hi(2), ng
    double precision :: velx(lo(1)-ng:,lo(2)-ng:)
    double precision :: vely(lo(1)-ng:,lo(2)-ng:)
    double precision :: dx,time

    ! local varables
    integer          :: i,j
    double precision :: x,y
    double precision :: psi(lo(1)-2:hi(1)+2,lo(2)-2:hi(2)+2)

    ! streamfunction psi
    do j = lo(2)-2, hi(2)+2
       y = (dble(j)+0.5d0)*dx
       do i = lo(1)-2, hi(1)+2
          x = (dble(i)+0.5d0)*dx
          psi(i,j) =  sin(M_PI*x)**2 * sin(M_PI*y)**2 * cos (M_PI*time/2.d0) / M_PI
       end do
    end do

    ! x velocity
    do j=lo(2)-1,hi(2)+1
       y = (dble(j)+0.5d0) * dx
       do i=lo(1)-1,hi(1)+2
          x = dble(i) * dx
          velx(i,j) =  -( (psi(i,j+1)+psi(i-1,j+1)) - (psi(i,j-1)+psi(i-1,j-1)) ) / (4.d0*dx)
       end do
    end do

    ! y velocity
    do j=lo(2)-1,hi(2)+2
       y = dble(j) * dx
       do i=lo(1)-1,hi(1)+1
          x = (dble(i)+0.5d0) * dx
          vely(i,j) = ( (psi(i+1,j)+psi(i+1,j-1)) - (psi(i-1,j)+psi(i-1,j-1)) ) / (4.d0*dx)
       end do
    end do

  end subroutine set_velocity_2d

  subroutine set_velocity_3d(velx, vely, velz, ng, lo, hi, dx, time)

    integer          :: lo(3), hi(3), ng
    double precision :: velx(lo(1)-ng:,lo(2)-ng:,lo(3)-ng:)
    double precision :: vely(lo(1)-ng:,lo(2)-ng:,lo(3)-ng:)
    double precision :: velz(lo(1)-ng:,lo(2)-ng:,lo(3)-ng:)
    double precision :: dx,time

    ! local varables
    integer          :: i,j,k
    double precision :: x,y
    double precision :: psi(lo(1)-2:hi(1)+2,lo(2)-2:hi(2)+2,lo(3)-1:hi(3)+1)

    ! streamfunction psi
    do k = lo(3)-1, hi(3)+1
       do j = lo(2)-2, hi(2)+2
          y = (dble(j)+0.5d0)*dx
          do i = lo(1)-2, hi(1)+2
             x = (dble(i)+0.5d0)*dx
             psi(i,j,k) =  sin(M_PI*x)**2 * sin(M_PI*y)**2 * cos (M_PI*time/2.d0) / M_PI
          end do
       end do
    end do

    ! x velocity
    do k=lo(3)-1,hi(3)+1
       do j=lo(2)-1,hi(2)+1
          y = (dble(j)+0.5d0) * dx
          do i=lo(1)-1,hi(1)+2
             x = dble(i) * dx
             velx(i,j,k) =  -( (psi(i,j+1,k)+psi(i-1,j+1,k)) - (psi(i,j-1,k)+psi(i-1,j-1,k)) ) / (4.d0*dx)
          end do
       end do
    end do

    ! y velocity
    do k=lo(3)-1,hi(3)+1
       do j=lo(2)-1,hi(2)+2
          y = dble(j) * dx
          do i=lo(1)-1,hi(1)+1
             x = (dble(i)+0.5d0) * dx
             vely(i,j,k) = ( (psi(i+1,j,k)+psi(i+1,j-1,k)) - (psi(i-1,j,k)+psi(i-1,j-1,k)) ) / (4.d0*dx)
          end do
       end do
    end do

    ! z velocity
    velz = 1.d0

  end subroutine set_velocity_3d

end module set_velocity_module


text:
      subroutine pkill(word,wordp,latt,mult,imon,nmon,istr,nstr,
     1                 lfno)
      use ffs
      use tffitcode
      character*(*) word,wordp
      logical mhogal,exist
      integer*8 latt(nlat)
      dimension mult(*)
      dimension imon(nmona,4),istr(nstra,4)
      external pack
      include 'inc/common.inc'
c
      call getwdl2(word,wordp)
      if(word.eq.'M') then
        if(nmon.eq.0) return
        call getwdl2(word,wordp)
        ni=ielm(wordp,exist)
        if(exist) then
          call getwdl2(word,wordp)
          nf=ielm(wordp,exist)
          if(exist) then
            call getwdl2(word,wordp)
          else
            nf=ni
          endif
          do 10 i=1,nmon
            if(mhogal(ni,nf,imon(imon(i,2),1))) then
              imon(imon(i,2),3)=1
            endif
   10     continue
          nmo=nmon
          call pack(imon(1,2),imon(1,3),nmon,nmon)
          write(lfno,'(2(A,I4))')
     1              '  BPM available:',nmon,' in ',nmo
          nmonact=nmon
        else
          return
        endif
      elseif(word.eq.'C') then
        if(nmon.eq.0) return
        call getwdl2(word,wordp)
        ni=ielm(wordp,exist)
        if(exist) then
          call getwdl2(word,wordp)
          nf=ielm(wordp,exist)
          if(exist) then
            call getwdl2(word,wordp)
          else
            nf=ni
          endif
          do 20 i=1,nstr
            if(mhogal(ni,nf,istr(istr(i,2),1))) then
              istr(istr(i,2),3)=1
            endif
   20     continue
          nst=nstr
          call pack(istr(1,2),istr(1,3),nstr,nst)
          write(lfno,'(2(A,I4))')
     1             '  Correction dipole available :',nstr,' in ',nst
        else
          return
        endif
      endif
      return
      end

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD