text:
module grid

use domain
use advection, only: NADV, NADVS

implicit none

character(6), parameter :: version = '6.10.9'
character(8), parameter :: version_date = 'Feb 2016'

integer, parameter :: nx = nx_gl/nsubdomains_x
integer, parameter :: ny = ny_gl/nsubdomains_y
integer, parameter :: nz = nz_gl+1
integer, parameter :: nzm = nz-1

integer, parameter :: nsubdomains = nsubdomains_x * nsubdomains_y

logical, parameter :: RUN3D = ny_gl.gt.1
logical, parameter :: RUN2D = .not.RUN3D

integer, parameter :: nxp1 = nx + 1
integer, parameter :: nyp1 = ny + 1 * YES3D
integer, parameter :: nxp2 = nx + 2
integer, parameter :: nyp2 = ny + 2 * YES3D
integer, parameter :: nxp3 = nx + 3
integer, parameter :: nyp3 = ny + 3 * YES3D
integer, parameter :: nxp4 = nx + 4
integer, parameter :: nyp4 = ny + 4 * YES3D

integer, parameter :: dimx1_u = -1                !!-1        -1        -1        -1
integer, parameter :: dimx2_u = nxp3              !!nxp3      nxp3      nxp3      nxp3
integer, parameter :: dimy1_u = 1-(2+NADV)*YES3D  !!1-5*YES3D 1-4*YES3D 1-3*YES3D 1-2*YES3D
integer, parameter :: dimy2_u = nyp2+NADV         !!nyp5      nyp4      nyp3      nyp2
integer, parameter :: dimx1_v = -1-NADV           !!-4        -3        -2        -1
integer, parameter :: dimx2_v = nxp2+NADV         !!nxp5      nxp4      nxp3      nxp2
integer, parameter :: dimy1_v = 1-2*YES3D         !!1-2*YES3D 1-2*YES3D 1-2*YES3D 1-2*YES3D
integer, parameter :: dimy2_v = nyp3              !!nyp3       nyp3      nyp3      nyp3
integer, parameter :: dimx1_w = -1-NADV           !!-4        -3        -2        -1
integer, parameter :: dimx2_w = nxp2+NADV         !!nxp5      nxp4      nxp3      nxp2
integer, parameter :: dimy1_w = 1-(2+NADV)*YES3D  !!1-5*YES3D 1-4*YES3D 1-3*YES3D 1-2*YES3D
integer, parameter :: dimy2_w = nyp2+NADV         !!nyp5      nyp4      nyp3      nyp2
integer, parameter :: dimx1_s = -2-NADVS          !!-4        -3        -2        -2
integer, parameter :: dimx2_s = nxp3+NADVS        !!nxp5      nxp4      nxp3      nxp3
integer, parameter :: dimy1_s = 1-(3+NADVS)*YES3D !!1-5*YES3D 1-4*YES3D 1-3*YES3D 1-3*YES3D
integer, parameter :: dimy2_s = nyp3+NADVS        !!nyp5      nyp4      nyp3      nyp3

integer, parameter :: ncols = nx*ny
integer, parameter :: nadams = 3

! Vertical grid parameters:
real z(nz)      ! height of the pressure levels above surface,m
real pres(nzm)  ! pressure,mb at scalar levels
real zi(nz)     ! height of the interface levels
real presi(nz)  ! pressure,mb at interface levels
real adz(nzm)   ! ratio of the thickness of scalar levels to dz
real adzw(nz)	! ratio of the thinckness of w levels to dz
real pres0      ! Reference surface pressure, Pa

integer:: nstep =0! current number of performed time steps
integer  ncycle  ! number of subcycles over the dynamical timestep
integer icycle  ! current subcycle
real cfl_adv    ! CFL due to advection
real cfl_sgs    ! CFL due to SGS
integer:: na=1, nb=2, nc=3 ! indeces for swapping the rhs arrays for AB scheme
real at, bt, ct ! coefficients for the Adams-Bashforth scheme
real dtn	! current dynamical timestep (can be smaller than dt)
real dt3(3) 	! dynamical timesteps for three most recent time steps
real(8):: time=0.	! current time in sec.
real day	! current day (including fraction)
real dtfactor   ! dtn/dt

!  MPI staff:
integer rank   ! rank of the current subdomain task (default 0)
integer ranknn ! rank of the "northern" subdomain task
integer rankss ! rank of the "southern" subdomain task
integer rankee ! rank of the "eastern"  subdomain task
integer rankww ! rank of the "western"  subdomain task
integer rankne ! rank of the "north-eastern" subdomain task
integer ranknw ! rank of the "north-western" subdomain task
integer rankse ! rank of the "south-eastern" subdomain task
integer ranksw ! rank of the "south-western" subdomain task
logical dompi  ! logical switch to do multitasking
logical masterproc ! .true. if rank.eq.0

character(80) case   ! id-string to identify a case-name(set in CaseName file)

logical dostatis     ! flag to permit the gathering of statistics
logical dostatisrad  ! flag to permit the gathering of radiation statistics
integer nstatis	! the interval between substeps to compute statistics

logical :: compute_reffc = .false.
logical :: compute_reffi = .false.
logical :: compute_reffl = .false. !bloss(2018-02): Include rain/drizzle as radiatively active

logical notopened2D  ! flag to see if the 2D output datafile is opened
logical notopened3D  ! flag to see if the 3D output datafile is opened
logical notopenedmom ! flag to see if the statistical moment file is opened

!-----------------------------------------
! Parameters controled by namelist PARAMETERS

real:: dx =0. 	! grid spacing in x direction
real:: dy =0.	! grid spacing in y direction
real:: dz =0.	! constant grid spacing in z direction (when dz_constant=.true.)
logical:: doconstdz = .false.  ! do constant vertical grid spacing set by dz

integer:: nstop =0   ! time step number to stop the integration
integer:: nelapse =999999999! time step number to elapse before stoping
integer:: nelapsemin=999999999 !bloss: number of wallclock minutes to run before stopping

real:: dt=0.	! dynamical timestep
real:: day0=0.	! starting day (including fraction)

integer:: nrad =1  ! frequency of calling the radiation routines
integer:: nprint =1000	! frequency of printing a listing (steps)
integer:: nrestart =0 ! switch to control starting/restarting of the model
integer:: nstat =1000	! the interval in time steps to compute statistics
integer:: nstatfrq =50 ! frequency of computing statistics

logical:: restart_sep =.false.  ! write separate restart files for sub-domains
integer:: nrestart_skip =0 ! number of skips of writing restart (default 0)
logical:: output_sep =.false.   ! write separate 3D and 2D files for sub-domains

character(80):: caseid =''! id-string to identify a run
character(80):: caseid_restart =''! id-string for branch restart file
character(80):: case_restart =''! id-string for branch restart file

logical:: doisccp = .false.
logical:: domodis = .false.
logical:: domisr = .false.
logical:: dosimfilesout = .false.

logical:: doSAMconditionals = .false. !core updraft,downdraft conditional statistics
logical:: dosatupdnconditionals = .false.!cloudy updrafts,downdrafts and cloud-free
logical:: doscamiopdata = .false.! initialize the case from a SCAM IOP netcdf input file
logical:: dozero_out_day0 = .false.
character(len=120):: iopfile=''
character(256):: rundatadir ='./RUNDATA' ! path to data directory

integer:: nsave3D =1000     ! frequency of writting 3D fields (steps)
integer:: nsave3Dstart =99999999! timestep to start writting 3D fields
integer:: nsave3Dend  =99999999 ! timestep to end writting 3D fields
logical:: save3Dbin =.true.   ! save 3D data in binary format(no 2-byte compression) -- Changed from SAM default of false
logical:: save3Dsep =.false.   ! use separate file for each time point for2-model
real   :: qnsave3D =0.    !threshold manimum cloud water(kg/kg) to save 3D fields
logical:: dogzip3D =.false.    ! gzip compress a 3D output file
logical:: rad3Dout = .false. ! output additional 3D radiation foelds (like reff)

integer:: nsave2D =1000     ! frequency of writting 2D fields (steps)
integer:: nsave2Dstart =99999999! timestep to start writting 2D fields
integer:: nsave2Dend =99999999  ! timestep to end writting 2D fields
logical:: save2Dbin =.true.   ! save 2D data in binary format, rather than compressed -- Changed from SAM default of false
logical:: save2Dsep =.false.   ! write separate file for each time point for 2D output
logical:: save2Davg =.false.   ! flag to time-average 2D output fields (default .false.)
logical:: dogzip2D =.false.    ! gzip compress a 2D output file if save2Dsep=.true.

integer:: nstatmom =1000! frequency of writting statistical moment fields (steps)
integer:: nstatmomstart = 1 ! timestep to start writting statistical moment fields
integer:: nstatmomend = 0  ! timestep to end writting statistical moment fields
logical:: savemomsep =.false.! use one file with stat moments for each time point
logical:: savemombin =.false.! save statistical moment data in binary format

integer:: nmovie =1000! frequency of writting movie fields (steps)
integer:: nmoviestart =99999999! timestep to start writting statistical moment fields
integer:: nmovieend =99999999  ! timestep to end writting statistical moment fields

logical :: isInitialized_scamiopdata = .false.
logical :: wgls_holds_omega = .false.

!bloss: options to change flux divergence in simple radiation scheme.
real :: rad_simple_fluxdiv1 = 70. ! W/m2
real :: rad_simple_fluxdiv2 = 22. ! W/m2
real :: rad_simple_kappa = 85. ! m2/kg

!bloss: option to use fixed divergence and to balance subsidence
!   heating above the inversion with cooling based on a fixed lapse
!   rate (used in Uchida papers)
logical :: dofixdivg = .false.
logical :: do_linear_subsidence = .true.  ! if both true, w_ls=divg_ls*z
real :: divg_ls = 3.75e-6 ! in 1/s, w_ls = -D*z, where [z] = m, [w] = m/s
real :: divg_lapse = 6.e-3 ! lapse rate of dry static energy above inversion

!bloss: option to use extrapolation within upperbound routine, even if
!  dolargescale==.true.  Useful for better consistency of water vapor
!  isotopologues at upper boundary
logical :: doExtrapolate_UpperBound = .false.

!bloss: option to fix wind speed used in calculation of surface fluxes
logical :: doFixedWindSpeedForSurfaceFluxes = .false.
real :: WindSpeedForFluxes = 8.

!bloss: option to run Derbyshire et al (2004) cases
logical :: doDerbyshire = .false.

! heights for bottom of nudging, top of transition layer, top of nudging
real :: Derbyshire_z1 = 1.e3, Derbyshire_z2 = 2.e3, Derbyshire_z3 = 15.e3

! Relative humidity for nudging between z1 and z2 (Low) and z2 and z3 (High)
!    In Derbyshire et al (2004), RelH_Low=0.8 for all cases and
!    RelH_High = 0.25, 0.5, 0.7 and 0.9
real :: Derbyshire_RelH_Low = 0.8, Derbyshire_RelH_High = 0.25

! Theta profile is defined as theta(z) = theta0 + LapseRate*z
!   with theta0 in K, LapseRate for theta in K/m.
real :: Derbyshire_theta0 = 290., Derbyshire_LapseRate = 3.e-3

real :: Derbyshire_tau = 3600. ! nudging timescale is one hour.

! option to require that advective fluxes be computed everywhere.
!   In ADV_UM5b, advective fluxes and updates are not computed
!   if there are chunks of the domain at the top/bottom with
!   only zero values.
logical :: compute_advection_everywhere = .false.

real :: wsub_ref(nz) = 0. ! used to output reference large-scale vertical motion
real :: w_wtg(nz) ! for statistics output

logical :: dowtg_blossey_etal_JAMES2009 = .false. ! alternate implementation
real :: am_wtg  = 2. ! momentum damping rate in 1/d -- note must be non-zero.
            ! default= 2.
real :: am_wtg_exp = 0. ! exponent of p/p0 in momentum damping rate.
real :: lambda_wtg = 650.e3 ! quarter wavelength in m. default = 650.e3 (=650 km).

logical :: dowtg_qnudge = .false. ! if T, nudge q profile to observations on below timescale
real :: itau_wtg_qnudge = 0. ! inverse nudging timescale for q in 1/day.

real :: tauz0_wtg_qnudge = 3.e3
real :: taulz_wtg_qnudge = 1.5e3

logical :: dowtg_tnudge = .false. ! enable temperature nudging above tropopause
real :: itau_wtg_tnudge = 0. ! inverse nudging timescale for temperature, 1/day
real :: taulz_wtg_tnudge = 1.5e3 ! length scale for ramping up t nudging in vertical.

! If true, use Dge computed within microphysics for computing radiative properties
!   of snow and cloud ice
logical :: doDge_SnowAndIce = .false.

! If true, use raw Thompson reff values as Dge for computing radiative properties
!   of snow and cloud ice.  I believe that this is inconsistent with the definition of
!   size in the radiation scheme, but it follows the approach in Thompson
!   et al (2016, Atmos Res).  Ignored if doDge_SnowAndIce==.true.
logical :: doThompsonReffIce = .false.

! Noah additions:
character(len=120) :: initial_condition_netcdf = ''

!-----------------------------------------
end module grid

max_stars_repo_name: kbren/uwnet

max_issues_repo_name: kbren/uwnet

max_forks_repo_name: kbren/uwnet


! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

      MODULE OUTPUT4_MATRICES

! This module defines the characteristics for matrices that are allowed for OUTPUT4 requests in Case Control. The definition of
! the actual matrices (with dimensions) appear in a variety of other modules (with the names listed in ALLOW_OU4_MYSTRAN_NAMES
! below. In addition, all of the output transformation matrices related to Craig-Bampton analyses are defined here. Finally, the
! allowable OUTPUT4 matrix RBM0 is defined here as a 6x6 rigid body mass matrix relative to the basic origin for the 6 DOF's at the
! boundary of a Craig-Bampton model

      USE PENTIUM_II_KIND, ONLY        :  BYTE, LONG, DOUBLE
      USE SCONTR, ONLY                 :  BLNK_SUB_NAM, TSET_CHR_LEN

      IMPLICIT NONE

      SAVE

      INTEGER(LONG), PARAMETER         :: NUM_OU4_VALID_NAMES = 33                   ! Num of OUTPUT4 matrices defined here
      INTEGER(LONG)                    :: NUM_OU4_REQUESTS    = 0                    ! Number of OUTPUT4 requests in a run
      INTEGER(LONG)                    :: NUM_PARTN_REQUESTS  = 0                    ! Number of PARTN requests in a run
      INTEGER(LONG), PRIVATE           :: II

! **********************************************************************************************************************************
! V A R I A B L E S   F O R   O U T P U T 4   M A T R I X   D E F I N I T I O N
! -----------------------------------------------------------------------------

      INTEGER(LONG), PARAMETER         :: TAPE_ACTION_MAX_VAL =  0                   ! See NASTRAN OUTPUT4 ITAPE allow values
      INTEGER(LONG), PARAMETER         :: TAPE_ACTION_MIN_VAL = -3                   ! See NASTRAN OUTPUT4 ITAPE allow values

      INTEGER(LONG)                    :: OU4_FILE_UNITS(NUM_OU4_VALID_NAMES)        ! File unit no's user has req for the mats
      INTEGER(LONG)                    :: OU4_TAPE_ACTION(NUM_OU4_VALID_NAMES)       ! Tape action (0 None, -1 rewind before write
!                                                                                -2 EOF and rewind after write, -3 both)
      CHARACTER(16*BYTE)               :: ACT_OU4_MYSTRAN_NAMES(NUM_OU4_VALID_NAMES) ! MYSTRAN names: act OUTPUT4 requested mats
      CHARACTER(16*BYTE)               :: ACT_OU4_OUTPUT_NAMES(NUM_OU4_VALID_NAMES)  ! Output  names: act OUTPUT4 requested mats

                                                           ! Following are MYSTRAN names that can be requested in Exec Control for
                                                           ! OUTPUT4. Corresponding NASTRAN names (incl for compatibility between
                                                           ! MYSTRAN and NASTRAN) are the ALLOW_OU4_OUTPUT_NAMES listed later
      CHARACTER(16*BYTE), PARAMETER    :: ALLOW_OU4_MYSTRAN_NAMES(NUM_OU4_VALID_NAMES)=                      &
                                                                                (/'CG_LTM          ', &! (  1)
                                                                                'DLR             ', &! (  2)
                                                                                'EIGEN_VAL       ', &! (  3)
                                                                                'EIGEN_VEC       ', &! (  4)
                                                                                'GEN_MASS        ', &! (  5)
                                                                                'IF_LTM          ', &! (  6)
                                                                                'KAA             ', &! (  7)
                                                                                'KGG             ', &! (  8)
                                                                                'KLL             ', &! (  9)
                                                                                'KRL             ', &! ( 10)
                                                                                'KRR             ', &! ( 11)
                                                                                'KRRcb           ', &! ( 12)
                                                                                'KXX             ', &! ( 13)
                                                                                'LTM             ', &! ( 14)
                                                                                'MCG             ', &! ( 15)
                                                                                'MEFFMASS        ', &! ( 16)
                                                                                'MPFACTOR        ', &! ( 17)
                                                                                'MAA             ', &! ( 18)
                                                                                'MGG             ', &! ( 19)
                                                                                'MLL             ', &! ( 20)
                                                                                'MRL             ', &! ( 21)
                                                                                'MRN             ', &! ( 22)
                                                                                'MRR             ', &! ( 23)
                                                                                'MRRcb           ', &! ( 24)
                                                                                'MXX             ', &! ( 25)
                                                                                'PA              ', &! ( 26)
                                                                                'PG              ', &! ( 27)
                                                                                'PL              ', &! ( 28)
                                                                                'PHIXG           ', &! ( 29)
                                                                                'PHIZG           ', &! ( 30)
                                                                                'RBM0            ', &! ( 31)
                                                                                'TR6_0           ', &! ( 32)
                                                                                'TR6_CG          ' /)! ( 33)

                                                           ! The following are allowable names Exec Control OUTPUT4 requests
                                                           ! in order to have compatibility with NASTRAN. The corresponding
                                                           ! MYSTRN names are the ALLOW_OU4_MYSTRAN_NAMES, above
      CHARACTER(16*BYTE), PARAMETER    :: ALLOW_OU4_OUTPUT_NAMES(NUM_OU4_VALID_NAMES)=                       &
                                                                                (/'CG_LTM          ', &! (  1) for CG_LTM
                                                                                'DM              ', &! (  2) for DLR
                                                                                'LAMA            ', &! (  3) for EIGEN_VAL
                                                                                'PHIG            ', &! (  4) for EIGEN_VEC
                                                                                'MI              ', &! (  5) for GEN_MASS
                                                                                'IF_LTM          ', &! (  6) for IF_LTM
                                                                                'KAA             ', &! (  7) for KAA
                                                                                'KGG             ', &! (  8) for KGG
                                                                                'KLL             ', &! (  9) for KLL
                                                                                'KLR(t)          ', &! ( 10) for KRL
                                                                                'KRR             ', &! ( 11) for KRR
                                                                                'KBB             ', &! ( 12) for KRRcb
                                                                                'KRRGN           ', &! ( 13) for KXX
                                                                                'LTM             ', &! ( 14) for LTM
                                                                                'RBMCG           ', &! ( 15) for MCG
                                                                                'MEFFMASS        ', &! ( 16) for MEFFMASS
                                                                                'MPFACTOR        ', &! ( 17) for MPFACTOR
                                                                                'MAA             ', &! ( 18) for MAA
                                                                                'MGG             ', &! ( 19) for MGG
                                                                                'MLL             ', &! ( 20) for MLL
                                                                                'MRL             ', &! ( 21) for MRL
                                                                                'MRN             ', &! ( 22) for MRN
                                                                                'MRR             ', &! ( 23) for MRR
                                                                                'MR              ', &! ( 24) for MRRcb
                                                                                'MRRGN           ', &! ( 25) for MXX
                                                                                'PA              ', &! ( 26) for PA
                                                                                'PG              ', &! ( 27) for PG
                                                                                'PL              ', &! ( 28) for PL
                                                                                'PHIXG           ', &! ( 29) for PHIXG
                                                                                'PHIZG           ', &! ( 30) for PHIZG
                                                                                'RBM0            ', &! ( 31) for RBM0
                                                                                'TR6_0           ', &! ( 32) for TR6_0
                                                                                'RBRCG           ' /)! ( 33) for TR6_CG

      CHARACTER(LEN(BLNK_SUB_NAM)), PARAMETER    :: SUBR_WHEN_TO_WRITE_OU4_MATS(NUM_OU4_VALID_NAMES) =       &
                                                                                (/'LINK6           ', &! (  1) for CG_LTM
                                                                                'LINK6           ', &! (  2) for DLR
                                                                                'LINK4           ', &! (  3) for EIGEN_VAL
                                                                                'LINK5           ', &! (  4) for EIGEN_VEC
                                                                                'LINK4           ', &! (  5) for GEN_MASS
                                                                                'LINK6           ', &! (  6) for IF_LTM
                                                                                'LINK2           ', &! (  7) for KAA
                                                                                'LINK2           ', &! (  8) for KGG
                                                                                'LINK2           ', &! (  9) for KLL
                                                                                'LINK2           ', &! ( 10) for KRL
                                                                                'LINK2           ', &! ( 11) for KRR
                                                                                'LINK6           ', &! ( 12) for KRRcb
                                                                                'LINK6           ', &! ( 13) for KXX
                                                                                'LINK6           ', &! ( 14) for LTM
                                                                                'LINK6           ', &! ( 15) for MCG
                                                                                'LINK9           ', &! ( 16) for MEFFMASS
                                                                                'LINK9           ', &! ( 17) for MPFACTOR
                                                                                'LINK2           ', &! ( 18) for MAA
                                                                                'LINK2           ', &! ( 19) for MGG
                                                                                'LINK2           ', &! ( 20) for MLL
                                                                                'LINK2           ', &! ( 21) for MRL
                                                                                'LINK6           ', &! ( 22) for MRN
                                                                                'LINK2           ', &! ( 23) for MRR
                                                                                'LINK6           ', &! ( 24) for MRRcb
                                                                                'LINK6           ', &! ( 25) for MXX
                                                                                'LINK2           ', &! ( 26) for PA
                                                                                'LINK2           ', &! ( 27) for PG
                                                                                'LINK2           ', &! ( 28) for PL
                                                                                'LINK5           ', &! ( 29) for PHIXG
                                                                                'LINK5           ', &! ( 30) for PHIZG
                                                                                'LINK6           ', &! ( 31) for RBM0
                                                                                'LINK6           ', &! ( 32) for TR6_0
                                                                                'LINK6           ' /)! ( 33) for TR6_CG

      CHARACTER(LEN(TSET_CHR_LEN)), PARAMETER                                                                                &
                                       :: OU4_MAT_COL_SETS(NUM_OU4_VALID_NAMES) =     (/'--'              , &! (  1) for CG_LTM
                                                                                'R '              , &! (  2) for DLR
                                                                                '--'              , &! (  3) for EIGEN_VAL
                                                                                '--'              , &! (  4) for EIGEN_VEC
                                                                                '--'              , &! (  5) for GEN_MASS
                                                                                '--'              , &! (  6) for IF_LTM
                                                                                'A '              , &! (  7) for KAA
                                                                                'G '              , &! (  8) for KGG
                                                                                'L '              , &! (  9) for KLL
                                                                                'L '              , &! ( 10) for KRL
                                                                                'R '              , &! ( 11) for KRR
                                                                                'R '              , &! ( 12) for KRRcb
                                                                                '--'              , &! ( 13) for KXX
                                                                                '--'              , &! ( 14) for LTM
                                                                                '--'              , &! ( 15) for MCG
                                                                                '--'              , &! ( 16) for MEFFMASS
                                                                                'R '              , &! ( 17) for MPFACTOR
                                                                                'A '              , &! ( 18) for MAA
                                                                                'G '              , &! ( 19) for MGG
                                                                                'L '              , &! ( 20) for MLL
                                                                                'L '              , &! ( 21) for MRL
                                                                                '--'              , &! ( 22) for MRN
                                                                                'R '              , &! ( 23) for MRR
                                                                                'R '              , &! ( 24) for MRRcb
                                                                                '--'              , &! ( 25) for MXX
                                                                                '--'              , &! ( 26) for PA
                                                                                '--'              , &! ( 27) for PG
                                                                                '--'              , &! ( 28) for PL
                                                                                '--'              , &! ( 29) for PHIXG
                                                                                '--'              , &! ( 30) for PHIZG
                                                                                '--'              , &! ( 31) for RBM0
                                                                                '--'              , &! ( 32) for TR6_0
                                                                                '--'               /)! ( 33) for TR6_CG

      CHARACTER(LEN(TSET_CHR_LEN)), PARAMETER                                                                                &
                                       :: OU4_MAT_ROW_SETS(NUM_OU4_VALID_NAMES) =     (/'--'              , &! (  1) for CG_LTM
                                                                                'L '              , &! (  2) for DLR
                                                                                '--'              , &! (  3) for EIGEN_VAL
                                                                                'G '              , &! (  4) for EIGEN_VEC
                                                                                '--'              , &! (  5) for GEN_MASS
                                                                                'R '              , &! (  6) for IF_LTM
                                                                                'A '              , &! (  7) for KAA
                                                                                'G '              , &! (  8) for KGG
                                                                                'L '              , &! (  9) for KLL
                                                                                'R '              , &! ( 10) for KRL
                                                                                'R '              , &! ( 11) for KRR
                                                                                'R '              , &! ( 12) for KRRcb
                                                                                '--'              , &! ( 13) for KXX
                                                                                '--'              , &! ( 14) for LTM
                                                                                '--'              , &! ( 15) for MCG
                                                                                '--'              , &! ( 16) for MEFFMASS
                                                                                '--'              , &! ( 17) for MPFACTOR
                                                                                'A '              , &! ( 18) for MAA
                                                                                'G '              , &! ( 19) for MGG
                                                                                'L '              , &! ( 20) for MLL
                                                                                'R '              , &! ( 21) for MRL
                                                                                'R '              , &! ( 22) for MRN
                                                                                'R '              , &! ( 23) for MRR
                                                                                'R '              , &! ( 24) for MRRcb
                                                                                '--'              , &! ( 25) for MXX
                                                                                'A '              , &! ( 26) for PA
                                                                                'G '              , &! ( 27) for PG
                                                                                'L '              , &! ( 28) for PL
                                                                                'G '              , &! ( 29) for PHIXG
                                                                                'G '              , &! ( 30) for PHIZG
                                                                                '--'              , &! ( 31) for RBM0
                                                                                'R '              , &! ( 32) for TR6_0
                                                                                'R '               /)! ( 33) for TR6_CG

      INTEGER(LONG), ALLOCATABLE       :: OU4_MAT_ROW_GRD_COMP(:,:)                 ! Grid/comp nos for rows in the part'd matrix
      INTEGER(LONG), ALLOCATABLE       :: OU4_MAT_COL_GRD_COMP(:,:)                 ! Grid/comp nos for cols in the part'd matrix

! 1st col of HAS_OU4_MAT_BEEN_PROCESSED is whether OU4 output request was processed. 2nd, 3rd col for whether partition request
! was processed for cols, rows

      CHARACTER(1*BYTE)                :: HAS_OU4_MAT_BEEN_PROCESSED(NUM_OU4_VALID_NAMES,3) =                                      &
                                                    RESHAPE ( (/('N', II=1,3*NUM_OU4_VALID_NAMES)/), (/NUM_OU4_VALID_NAMES,3/) )

! **********************************************************************************************************************************
! V A R I A B L E S   F O R   P A R T I T I O N I N G   O U T P U T 4   M A T R I C E S
! -------------------------------------------------------------------------------------

      CHARACTER(16*BYTE)               :: OU4_PART_MAT_NAMES(NUM_OU4_VALID_NAMES,5)  ! Names of matrices (from PARTN entry)

      CHARACTER(16*BYTE)               :: OU4_PART_VEC_NAMES(NUM_OU4_VALID_NAMES,2)  ! Row/col partition vec names for OU4 matrices

      INTEGER(LONG), ALLOCATABLE       :: OU4_PARTVEC_COL(:)                         ! Col partition vector for an OU4 matrix
      INTEGER(LONG), ALLOCATABLE       :: OU4_PARTVEC_ROW(:)                         ! Row partition vector for an OU4 matrix

! **********************************************************************************************************************************
! V A R I A B L E S   F O R   O T M's   F O R   O U T P U T 4   M A T R I C E S

! In a Craig-Bampton model generation run, Case Control output requests for displacements, element forces, stresses, etc are
! interpreted as requests for the Output Transformation Matrix (OTM) for that quantity. See Appendix D to the MYSTRAN User's
! Reference Manual for a discussion of Craig-Bampton model OTM's

      CHARACTER(132*BYTE), ALLOCATABLE :: TXT_ACCE(:)      ! CB accel     text descriptor for rows of OTM_ACCE
      CHARACTER(132*BYTE), ALLOCATABLE :: TXT_DISP(:)      ! CB disp      text descriptor for rows of OTM_DISP
      CHARACTER(132*BYTE), ALLOCATABLE :: TXT_MPCF(:)      ! CB MPC force text descriptor for rows of OTM_MPCF
      CHARACTER(132*BYTE), ALLOCATABLE :: TXT_SPCF(:)      ! CB SPC force text descriptor for rows of OTM_SPCF

      CHARACTER(132*BYTE), ALLOCATABLE :: TXT_ELFE(:)      ! CB elm engr force text descriptor for rows of OTM_ELFE
      CHARACTER(132*BYTE), ALLOCATABLE :: TXT_ELFN(:)      ! CB elm node force text descriptor for rows of OTM_ELFN
      CHARACTER(132*BYTE), ALLOCATABLE :: TXT_STRE(:)      ! CB elm stress     text descriptor for rows of OTM_STRE
      CHARACTER(132*BYTE), ALLOCATABLE :: TXT_STRN(:)      ! CB elm strain     text descriptor for rows of OTM_STRE

      REAL(DOUBLE), ALLOCATABLE        :: OTM_ACCE(:,:)    ! CB accel     OTM for grids req in Case Control ACCE = sid entry
      REAL(DOUBLE), ALLOCATABLE        :: OTM_DISP(:,:)    ! CB disp      OTM for grids req in Case Control DISP = sid entry
      REAL(DOUBLE), ALLOCATABLE        :: OTM_MPCF(:,:)    ! CB MPC force OTM for grids req in Case Control MPCF = sid entry
      REAL(DOUBLE), ALLOCATABLE        :: OTM_SPCF(:,:)    ! CB SPC force OTM for grids req in Case Control SPCF = sid entry

      REAL(DOUBLE), ALLOCATABLE        :: OTM_ELFE(:,:)    ! CB elm engr force OTM for elems req in CC ELFORCE(ENGR) = sid entry
      REAL(DOUBLE), ALLOCATABLE        :: OTM_ELFN(:,:)    ! CB elm node force OTM for elems req in CC ELFORCE(NODE) = sid entry
      REAL(DOUBLE), ALLOCATABLE        :: OTM_STRE(:,:)    ! CB elm stress     OTM for elems req in CC STRESS        = sid entry
      REAL(DOUBLE), ALLOCATABLE        :: OTM_STRN(:,:)    ! CB elm strain     OTM for elems req in CC STRESS        = sid entry

! **********************************************************************************************************************************
! O T H E R

      REAL(DOUBLE)                     :: RBM0(6,6)        ! Rigid body mass matrix rel to basic origin for the 6 DOF's at the
!                                                            boundary of a Craig-Bampton model

      END MODULE OUTPUT4_MATRICES


text:
!> Unit tests for particle types.
module particle_pfunit
  use m_point, only: point, particlearray_to_json
  use json_module
  use pfunit
  use m_rod, only: rod
  use m_particlejson_parser
  implicit none

contains

  subroutine test_rod_json_io
    type(json_value), pointer :: json_val
    type(rod) :: out, in

    out%x = 1.
    out%y = 2.
    out%z = 3.
    out%ux = -1./sqrt(2.)
    out%ux = 1./sqrt(2.)
    out%uz = 0.
    !! write coordinates to json
    call out%coordinates_to_json(json_val)
    !! read coordinates from json
    call in%from_json(json_val)
    call assertTrue(out == in, 'Input and output rod not equal.')

  end subroutine test_rod_json_io

  subroutine test_point_json_io
    type(json_value), pointer :: json_val
    type(point) :: out, in

    out%x = 1.
    out%y = 2.
    out%z = 3.
    !! write coordinates to json
    call out%coordinates_to_json(json_val)

    !! read coordinates from json
    call in%from_json(json_val)
    call assertTrue(out == in, 'Input and output point not equal.')
  end subroutine test_point_json_io

  subroutine test_rodarray_json_io
    type(rod) :: out(2)
    class(point), allocatable :: in(:)
    type(json_value), pointer :: json_val

    call json_create_object(json_val, 'test_particlearray_json_io')
    call particlearray_to_json(json_val, out)
    call particlearray_from_json(json_val, in)

    select type(in)
    type is (rod)
       call assertTrue(all(out(:) == in(1:2)), &
            'Input and output particlearrays not equal.')
    class default
       call assertTrue(.false., 'Input is not type rod.')
    end select

  end subroutine test_rodarray_json_io

  subroutine test_particlearray_json_io
    type(point) :: out(2)
    class(point), allocatable :: in(:)
    type(json_value), pointer :: json_val

    call json_create_object(json_val, 'test_particlearray_json_io')
    call particlearray_to_json(json_val, out)
    call particlearray_from_json(json_val, in)

    select type(in)
    type is (point)
       call assertTrue(all(out(:) == in(1:2)), &
            'Input and output particlearrays not equal.')
    class default
       call assertTrue(.false., 'Input is not of type point.')
    end select

  end subroutine test_particlearray_json_io

  subroutine test_point_downcast
    class(point), allocatable :: somepoint
    class(point), allocatable :: another
    integer :: err
    allocate(somepoint, source=point(x=1.0, y=2.0, z=3.0))
    allocate(another, source=point())
    call another%downcast_assign(somepoint, err)
    call assertEqual(0, err, 'point downcast failed')
    select type (somepoint)
    type is (point)
       select type (another)
       type is (point)
          call assertTrue(somepoint == another, 'assignment did not work')
       class default
          call assertTrue(.false., 'another is not point')
       end select
    class default
       call assertTrue(.false., 'somepoint is not point')
    end select
  end subroutine test_point_downcast

  subroutine test_rod_downcast
    class(point), allocatable :: somerod
    class(point), allocatable :: another
    integer :: err
    allocate(somerod, source=rod(x=1.0, ux=1.0, uz=0.0))
    allocate(another, source=rod())
    call another%downcast_assign(somerod, err)
    call assertEqual(0, err, 'rod downcast failed')
    select type (somerod)
    type is (rod)
       select type (another)
       type is (rod)
          call assertTrue(somerod == another, 'assignment did not work')
       class default
          call assertTrue(.false., 'another is not rod')
       end select
    class default
       call assertTrue(.false., 'somerod is not rod')
    end select

  end subroutine test_rod_downcast

end module particle_pfunit

max_stars_repo_name: jonekoo/CoarseMC

max_issues_repo_name: jonekoo/CoarseMC

max_forks_repo_name: jonekoo/CoarseMC


text:
      SUBROUTINE VER4_27
      CONTINUE
      RETURN
      END

ext: for

max_stars_repo_name: zhumiansanren/FIRI-WHS

max_issues_repo_name: zhumiansanren/FIRI-WHS

max_forks_repo_name: zhumiansanren/FIRI-WHS


text:
      FUNCTION K51R(J,M,K,L,R,JA)
C BY JOHN WILLIAM COLBERT
C TYPED BY JM
      K51R=JA     ! Added 7/32/93 to report own old value
      IF (M .GT. 8) GOTO 5
        K51R = 0
      IF (M .EQ. 6) K51R = 1
        LASTI = 0
        GOTO 10
5     K51R = 0
      LASTI = LASTI - 1
      IF (LASTI .EQ. 3) K51R = 1
      IF (LASTI .GT. 0) GOTO 10
      IF (J .EQ. 1) K51R = 1
      IF (J .EQ. 1) LASTI = 4
10    RETURN
      END

ext: f

max_stars_repo_name: Axelrod-Python/TourExec

max_issues_repo_name: Axelrod-Python/TourExec

max_forks_repo_name: Axelrod-Python/TourExec


text:
The International Festival is an annual festival promoting unity in diversity with free admission for everyone. Its a wonderful and educational event when you can go around the world in one day without leaving Davis!
The festival will be held at the Veterans Memorial Center at 203 E. 14th St. on Saturday, October 12, 2013 from 12:00pm to 6:00pm.

This is a fundraiser to support International House Davis, a community based nonprofit organization that promotes a global community by providing opportunities for crosscultural interaction and exchange.

 What will happen at the International Festival?
  Dancing, music, performances, food, educational lectures and information, and lots of youth activities.
  More than 20 country tables devoted to education about a particular country so you can visit a new culture without leaving the Festival.
  Ethnic marketplace where vendors sell international wares.
  Food for sale!  No alcohol will be sold.
  The youth corner is full of things for young people to do, see, and even create.
  Fun for the whole family!

The 2011 festival coincides with IHouse celebrating its 30th year welcoming international students and scholars to Davis and hosting current events lectures, languages classes and cultural events for community members at 10 College Park.
 To get involved or for more information
Check out our website:  http://www.internationalfestivaldavis.org.

20111001 17:09:34 nbsp The 5 year old and I went to this today. Vets Hall was really too crowded, dangerously so. The performances were terrific. I dont know why they couldnt have performed in the theater though. It would be more comfortable, all the patrons would be able to see, and it would free up room in the hall. Users/NoelBruening

20130814 13:26:14 nbsp Thank you for your feedback Noel. The festival is a multifaceted event with 8 venues to participate in, including the main stage. We will be exploring moving to another site in the future to accommodate the enthusiasm and growing attendance. Last year, we tried having entertainment in the theater and no one attended due to the fact that the weather was beautiful and there were other venues going on outside. Hence our performers in the theatre were very disappointed even though they loved the building. We hope you attend again this year, but we appreciate your feedback. Sincerely, International House. Users/ChardaeJones

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module probdata_module

!     Sod variables
      use amrex_fort_module, only : rt => amrex_real
      real(rt)        , save ::  p_l, u_l, rho_l, p_r, u_r, rho_r, rhoe_l, rhoe_r, frac
      real(rt)        , save :: T_l, T_r

      logical, save :: use_Tinit

!     These help specify which specific problem
      integer        , save ::  probtype,idir

      real(rt)        , save ::  split(3)

end module probdata_module

max_stars_repo_name: yingtchen/Castro

max_issues_repo_name: yingtchen/Castro

max_forks_repo_name: yingtchen/Castro


text: The Margaret B and Michael Harrison Seating Area is on the north side of Shields Library facing the Quad.


text:
program testmpfun

!   This briefly tests most individual MPFUN operations and functions
!   (including all mixed mode arithmetic and comparison operations), by
!   comparing each result with the equivalent operation performed in double
!   precision or double complex as appropriate.  This is clearly not an
!   exhaustive or fully high-precision test, but it often detects software bugs
!   and compiler issues.  All multiprecision results are output to 500 digits,
!   so that the output file can be compared with a reference file.

!   David H Bailey  19 Apr 2017

use mpmodule
implicit none
integer i, ndp, nwds
parameter (ndp = 500, nwds = ndp / mpdpw + 2)
character(500) chr500
character(32) chr32
character(1) chr1(ndp+30), chrpi(ndp+30)
integer i1, i2, i3, i4
logical l1, l2, l3, l4, l5, l6, l7, l8, l9, l10
real (8) d1, d2, d3, d4, e1, e2
complex (8) dc1, dc2, dc3, dc4, ec1, ec2
type (mp_real) t1, t2, t3, t4
type (mp_complex) z1, z2, z3, z4

write (6, '(a)') 'MPUN2015 quick check of operations and functions'

!   Define a few sample data values.

chr500 = &
'3.14159265358979323846264338327950288419716939937510582097494459230781&
&6406286208998628034825342117067982148086513282306647093844609550582231&
&7253594081284811174502841027019385211055596446229489549303819644288109&
&7566593344612847564823378678316527120190914564856692346034861045432664&
&8213393607260249141273724587006606315588174881520920962829254091715364&
&3678925903600113305305488204665213841469519415116094330572703657595919&
&5309218611738193261179310511854807446237996274956735188575272489122793&
&81830119491'

do i = 1, ndp
  chrpi(i) = chr500(i:i)
enddo

t1 = mppi (nwds)
t2 = - mplog2 (nwds)
d1 = t1
d2 = t2
e1 = 3141.d0 / 8192.d0
e2 = 6931.d0 / 8192.d0
z1 = mpcmplx (0.5d0 * mppi (nwds), exp (mpreal (0.5d0, nwds)), nwds)
z2 = mpcmplx (- gamma (mpreal (0.5d0, nwds)), cos (mpreal (1.d0, nwds)), nwds)
dc1 = z1
dc2 = z2
ec1 = cmplx (e1, e2, mprknd)
ec2 = cmplx (-e2, e1, mprknd)
i1 = 5
i2 = -3

write (6, '(/a/)') 'Data arrays:'
write (6, '(a)') 't1 = pi:'
call mpwrite (6, ndp + 20, ndp, t1)
write (6, '(a)') 't2 = -log(2):'
call mpwrite (6, ndp + 20, ndp, t2)
write (6, '(a)') 'z1 = (0.5*pi, exp(0.5)):'
call mpwrite (6, ndp + 20, ndp, z1)
write (6, '(a)') 'z2 = (-Gamma(0.5), Cos(1)):'
call mpwrite (6, ndp + 20, ndp, z2)

write (6, '(a)') 'e1 = 3141/8192:'
write (6, '(1p,d25.15)') e1
write (6, '(a)') 'e2 = 6931 / 8192:'
write (6, '(1p,d25.15)') e2
write (6, '(a)')  'ec1 = (3141/8192, 6931/8192)'
write (6, '(1p,2d25.15)') ec1
write (6, '(a)') 'ec2 = (6931/8192, 3141/8192):'
write (6, '(1p,2d25.15)') ec2

write (6, '(/a/)') 'Real data operations:'

write (6, '(a)') 'addition: t1+t2 ='
call mpwrite (6, ndp + 20, ndp, t1 + t2)
write (6, 1) d1 + d2
1 format (1p,2d25.15)
call checkdp (t1 + t2, d1 + d2, ndp)

write (6, '(a)') 'addition: t1+e2 ='
call mpwrite (6, ndp + 20, ndp, t1 + e2)
write (6, 1) d1 + e2
call checkdp (t1 + e2, d1 + e2, ndp)

write (6, '(a)') 'addition: e1+t2 ='
call mpwrite (6, ndp + 20, ndp, e1 + t2)
write (6, 1) e1 + d2
call checkdp (e1 + t2, e1 + d2, ndp)

write (6, '(a)') 'subtraction: t1-t2 ='
call mpwrite (6, ndp + 20, ndp, t1 - t2)
write (6, 1) d1 - d2
call checkdp (t1 - t2, d1 - d2, ndp)

write (6, '(a)') 'subtraction: t1-e2 ='
call mpwrite (6, ndp + 20, ndp, t1 - e2)
write (6, 1) d1 - e2
call checkdp (t1 - e2, d1 - e2, ndp)

write (6, '(a)') 'subtraction: e1-t2 ='
call mpwrite (6, ndp + 20, ndp, e1 - t2)
write (6, 1) e1 - d2
call checkdp (e1 - t2, e1 - d2, ndp)

write (6, '(a)') 'multiplication: t1*t2 ='
call mpwrite (6, ndp + 20, ndp, t1 * t2)
write (6, 1) d1 * d2
call checkdp (t1 * t2, d1 * d2, ndp)

write (6, '(a)') 'multiplication: t1*e2 ='
call mpwrite (6, ndp + 20, ndp, t1 * e2)
write (6, 1) d1 * e2
call checkdp (t1 * e2, d1 * e2, ndp)

write (6, '(a)') 'multiplication: e1*t2 ='
call mpwrite (6, ndp + 20, ndp, e1 * t2)
write (6, 1) e1 * d2
call checkdp (e1 * t2, e1 * d2, ndp)

write (6, '(a)') 'division: t1/t2 ='
call mpwrite (6, ndp + 20, ndp, t1 / t2)
write (6, 1) d1 / d2
call checkdp (t1 / t2, d1 / d2, ndp)

write (6, '(a)') 'division: t1/e2 ='
call mpwrite (6, ndp + 20, ndp, t1 / e2)
write (6, 1) d1 / e2
call checkdp (t1 / e2, d1 / e2, ndp)

write (6, '(a)') 'division: e1/t2 ='
call mpwrite (6, ndp + 20, ndp, e1 / t2)
write (6, 1) e1 / d2
call checkdp (e1 / t2, e1 / d2, ndp)

write (6, '(a)') 'exponentiation: t1**i1 ='
call mpwrite (6, ndp + 20, ndp, t1 ** i1)
write (6, 1) d1 ** i1
call checkdp (t1 ** i1, d1 ** i1, ndp)

write (6, '(a)') 'exponentiation: t1**t2 ='
call mpwrite (6, ndp + 20, ndp, t1 ** t2)
write (6, 1) d1 ** d2
call checkdp (t1 ** t2, d1 ** d2, ndp)

write (6, '(a)') 'equal test: t1 == t2, e1 == t2, t1 == e2'
l1 = t1 == t2; l2 = d1 == d2; l3 = e1 == t2; l4 = e1 == d2; l5 = t1 == e2; l6 = d1 == e2
write (6, '(6l4)') l1, l2, l3, l4, l5, l6
call checkl (l1, l2); call checkl (l3, l4); call checkl (l5, l6)

write (6, '(a)') 'not-equal test: t1 /= t2, e1 /= t2, t1 =/ e2'
l1 = t1 /= t2; l2 = d1 /= d2; l3 = e1 /= t2; l4 = e1 /= d2; l5 = t1 /= e2; l6 = d1 /= e2
write (6, '(6l4)') l1, l2, l3, l4, l5, l6
call checkl (l1, l2); call checkl (l3, l4); call checkl (l5, l6)

write (6, '(a)') 'less-than-or-equal test: t1 <= t2, e1 <= t2, t1 <= e2'
l1 = t1 <= t2; l2 = d1 <= d2; l3 = e1 <= t2; l4 = e1 <= d2; l5 = t1 <= e2; l6 = d1 <= e2
write (6, '(6l4)') l1, l2, l3, l4, l5, l6
call checkl (l1, l2); call checkl (l3, l4); call checkl (l5, l6)

write (6, '(a)') 'greater-than-or-equal test: t1 >= t2, e1 >= t2, t1 >= e2'
l1 = t1 >= t2; l2 = d1 >= d2; l3 = e1 >= t2; l4 = e1 >= d2; l5 = t1 >= e2; l6 = d1 >= e2
write (6, '(6l4)') l1, l2, l3, l4, l5, l6
call checkl (l1, l2); call checkl (l3, l4); call checkl (l5, l6)

write (6, '(a)') 'less-than test: t1 < t2, e1 < t2, t1 < e2'
l1 = t1 < t2; l2 = d1 < d2; l3 = e1 < t2; l4 = e1 < d2; l5 = t1 < e2; l6 = d1 < e2
write (6, '(6l4)') l1, l2, l3, l4, l5, l6
call checkl (l1, l2); call checkl (l3, l4); call checkl (l5, l6)

write (6, '(a)') 'greater-than test: t1 > t2, e1 > t2, t1 > e2'
l1 = t1 > t2; l2 = d1 > d2; l3 = e1 > t2; l4 = e1 > d2; l5 = t1 > e2; l6 = d1 > e2
write (6, '(6l4)') l1, l2, l3, l4, l5, l6
call checkl (l1, l2); call checkl (l3, l4); call checkl (l5, l6)

write (6, '(a)') 'abs(t2) ='
call mpwrite (6, ndp + 20, ndp, abs (t2))
write (6, '(1p,2d25.15)') abs (d2)
call checkdp (abs (t2), abs (d2), ndp)

write (6, '(a)') 'acos(t2) ='
call mpwrite (6, ndp + 20, ndp, acos (t2))
write (6, '(1p,2d25.15)') acos (d2)
call checkdp (acos (t2), acos (d2), ndp)

write (6, '(a)') 'aint(t1) ='
call mpwrite (6, ndp + 20, ndp, aint (t1))
write (6, '(1p,2d25.15)') aint (d1)
call checkdp (aint (t1), aint (d1), ndp)

write (6, '(a)') 'anint(t1) ='
call mpwrite (6, ndp + 20, ndp, anint (t1))
write (6, '(1p,2d25.15)') anint (d1)
call checkdp (anint (t1), anint (d1), ndp)

write (6, '(a)') 'asin(t2) ='
call mpwrite (6, ndp + 20, ndp, asin (t2))
write (6, '(1p,2d25.15)') asin (d2)
call checkdp (asin (t2), asin (d2), ndp)

write (6, '(a)') 'atan(t1) ='
call mpwrite (6, ndp + 20, ndp, atan (t1))
write (6, '(1p,2d25.15)') atan (d1)
call checkdp (atan (t1) , atan (d1), ndp)

write (6, '(a)') 'atan2(t1,t2) ='
call mpwrite (6, ndp + 20, ndp, atan2 (t1,t2))
write (6, '(1p,2d25.15)') atan2 (d1,d2)
call checkdp (atan2 (t1, t2), atan2 (d1, d2), ndp)

write (6, '(a)') 'cos(t2) ='
call mpwrite (6, ndp + 20, ndp, cos (t2))
write (6, '(1p,2d25.15)') cos (d2)
call checkdp (cos (t2), cos (d2), ndp)

write (6, '(a)') 'cosh(t1) ='
call mpwrite (6, ndp + 20, ndp, cosh (t1))
write (6, '(1p,2d25.15)') cosh (d1)
call checkdp (cosh (t1), cosh (d1), ndp)

write (6, '(a)') 'exp(t1) ='
call mpwrite (6, ndp + 20, ndp, exp (t1))
write (6, '(1p,2d25.15)') exp (d1)
call checkdp (exp (t1), exp (d1), ndp)

write (6, '(a)') 'gamma(t1) ='
call mpwrite (6, ndp + 20, ndp, gamma (t1))
write (6, '(1p,2d25.15)') gamma (d1)
call checkdp (gamma (t1), gamma (d1), ndp)

write (6, '(a)') 'log(t1) ='
call mpwrite (6, ndp + 20, ndp, log (t1))
write (6, '(1p,2d25.15)') log (d1)
call checkdp (log (t1), log (d1), ndp)

write (6, '(a)') 'max(t1,t2) ='
call mpwrite (6, ndp + 20, ndp, max (t1,t2))
write (6, '(1p,2d25.15)') max (d1,d2)
call checkdp (max (t1, t2), max (d1, d2), ndp)

write (6, '(a)') 'min(t1,t2) ='
call mpwrite (6, ndp + 20, ndp, min (t1,t2))
write (6, '(1p,2d25.15)') min (d1,d2)
call checkdp (min (t1, t2), min (d1, d2), ndp)

write (6, '(a)') 'mpcssh(t1) ='
call mpcssh (t1, t3, t4)
call mpwrite (6, ndp + 20, ndp, t3, t4)
d3 = cosh (d1)
d4 = sinh (d1)
write (6, '(1p,2d25.15)') d3, d4
call checkdp (t3, d3, ndp)
call checkdp (t4, d4, ndp)

write (6, '(a)') 'mpcssn(t2) ='
call mpcssn (t2, t3, t4)
call mpwrite (6, ndp + 20, ndp, t3, t4)
d3 = cos (d2)
d4 = sin (d2)
write (6, '(1p,2d25.15)') d3, d4
call checkdp (t3, d3, ndp)
call checkdp (t4, d4, ndp)

write (6, '(a)') 'mpeform(t1) ='
call mpeform (t1, ndp + 20, ndp, chr1)
write (6, '(80a1)') (chr1(i), i = 1, ndp + 20)
write (6, '(1p,d25.15)') d1
call checkdp (mpreal (chr1, ndp, nwds), d1, nwds)

write (6, '(a)') 'mpfform(t1) ='
call mpfform (t1, ndp + 20, ndp, chr1)
write (6, '(80a1)') (chr1(i), i = 1, ndp + 20)
write (6, '(f25.15)') d1
call checkdp (mpreal (chr1, ndp + 20, nwds), d1, nwds)

write (6, '(a)') 'mpmdi(t1) ='
call mpmdi (t1, d4, i4)
write (6, '(1p,d25.15,i8/d25.15)') d4, i4, d4*2**i4
write (6, '(1p,2d25.15)') d1
call checkdp (t1, d4*2**i4, ndp)

write (6, '(a)') 'mpnrt(t1,i1) ='
call mpwrite (6, ndp + 20, ndp, mpnrt (t1,i1))
write (6, '(1p,2d25.15)') d1 ** (1.d0 / i1)
call checkdp (mpnrt (t1, i1), d1 ** (1.d0 / i1), ndp)

write (6, '(a)') 'mpprodd (t1,d2) ='
call mpwrite (6, ndp + 20, ndp, mpprodd (t1,d2))
write (6, '(1p,2d25.15)') d1*d2
call checkdp (mpprodd (t1,d2), d1 * d2, ndp)

write (6, '(a)') 'mpquotd (t1,d2) ='
call mpwrite (6, ndp + 20, ndp, mpquotd (t1,d2))
write (6, '(1p,2d25.15)') d1/d2
call checkdp (mpquotd (t1,d2), d1 / d2, ndp)

write (6, '(a)') 'mpreal (chr500, nwds) ='
call mpwrite (6, ndp + 20, ndp, mpreal (chr500, nwds))
read (chr500, *) d4
write (6, '(1p,2d25.15)') d4
call checkdp (mpreal (chr500, nwds), d4, ndp)

write (6, '(a)') 'mpreal (chrpi, ndp, nwds) ='
call mpwrite (6, ndp + 20, ndp, mpreal (chrpi, ndp, nwds))
read (chr500, *) d4
write (6, '(1p,2d25.15)') d4
call checkdp (mpreal (chrpi, ndp, nwds), d4, ndp)

write (6, '(a)') 'sign(t1,t2) ='
call mpwrite (6, ndp + 20, ndp, sign (t1,t2))
write (6, '(1p,2d25.15)') sign (d1, d2)
call checkdp (sign (t1, t2), sign (d1, d2), ndp)

write (6, '(a)') 'sin(t2) ='
call mpwrite (6, ndp + 20, ndp, sin (t2))
write (6, '(1p,2d25.15)') sin (d2)
call checkdp (sin (t2), sin (d2), ndp)

write (6, '(a)') 'sinh(t1) ='
call mpwrite (6, ndp + 20, ndp, sinh (t1))
write (6, '(1p,2d25.15)') sinh (d1)
call checkdp (sinh (t1), sinh (d1), ndp)

write (6, '(a)') 'sqrt(t1) ='
call mpwrite (6, ndp + 20, ndp, sqrt (t1))
write (6, '(1p,2d25.15)') sqrt (d1)
call checkdp (sqrt (t1), sqrt (d1), ndp)

write (6, '(a)') 'tan(t1) ='
call mpwrite (6, ndp + 20, ndp, tan (t2))
write (6, '(1p,2d25.15)') tan (d2)
call checkdp (tan (t2), tan (d2), ndp)

write (6, '(a)') 'tanh(t1) ='
call mpwrite (6, ndp + 20, ndp, tanh (t1))
write (6, '(1p,2d25.15)') tanh (d1)
call checkdp (tanh (t1), tanh (d1), ndp)

! write (6, '(a)') 'zeta(t1) ='
! call mpwrite (6, ndp + 20, ndp, zeta (t1))
! write (6, '(1p,2d25.15)') zeta (d1)
! call checkdp (zeta (t1), zeta (d1))

write (6, '(/a/)') 'Complex data operations:'

write (6, '(a)') 'addition: z1+z2 ='
call mpwrite (6, ndp + 20, ndp, z1 + z2)
write (6, 1) dc1 + dc2
call checkdc (z1 + z2, dc1 + dc2, ndp)

write (6, '(a)') 'addition: z1+e2 ='
call mpwrite (6, ndp + 20, ndp, z1 + e2)
write (6, 1) dc1 + e2
call checkdc (z1 + e2, dc1 + e2, ndp)

write (6, '(a)') 'addition: e1+z2 ='
call mpwrite (6, ndp + 20, ndp, e1 + z2)
write (6, 1) e1 + dc2
call checkdc (e1 + z2, e1 + dc2, ndp)

call mpwrite (6, ndp + 20, ndp, z1 + ec2)
write (6, 1) dc1 + ec2
call checkdc (z1 + ec2, dc1 + ec2, ndp)

write (6, '(a)') 'addition: ec1+z2 ='
call mpwrite (6, ndp + 20, ndp, ec1 + z2)
write (6, 1) ec1 + dc2
call checkdc (ec1 + z2, ec1 + dc2, ndp)

write (6, '(a)') 'addition: z1+t2 ='
call mpwrite (6, ndp + 20, ndp, z1 + t2)
write (6, 1) dc1 + d2
call checkdc (z1 + t2, dc1 + d2, ndp)

write (6, '(a)') 'addition: t1+z2 ='
call mpwrite (6, ndp + 20, ndp, t1 + z2)
write (6, 1) d1 + dc2
call checkdc (t1 + z2, d1 + dc2, ndp)

write (6, '(a)') 'subtraction: z1-z2 ='
call mpwrite (6, ndp + 20, ndp, z1 - z2)
write (6, 1) dc1 - dc2
call checkdc (z1 - z2, dc1 - dc2, ndp)

write (6, '(a)') 'subtraction: z1-e2 ='
call mpwrite (6, ndp + 20, ndp, z1 - e2)
write (6, 1) dc1 - e2
call checkdc (z1 - e2, dc1 - e2, ndp)

write (6, '(a)') 'subtraction: e1-z2 ='
call mpwrite (6, ndp + 20, ndp, e1 - z2)
write (6, 1) e1 - dc2
call checkdc (e1 - z2, e1 - dc2, ndp)

write (6, '(a)') 'subtraction: z1-ec2 ='
call mpwrite (6, ndp + 20, ndp, z1 - ec2)
write (6, 1) dc1 - ec2
call checkdc (z1 - ec2, dc1 - ec2, ndp)

write (6, '(a)') 'subtraction: ec1-z2 ='
call mpwrite (6, ndp + 20, ndp, ec1 - z2)
write (6, 1) ec1 - dc2
call checkdc (ec1 - z2, ec1 - dc2, ndp)

write (6, '(a)') 'subtraction: z1-t2 ='
call mpwrite (6, ndp + 20, ndp, z1 - t2)
write (6, 1) dc1 - d2
call checkdc (z1 - t2, dc1 - d2, ndp)

write (6, '(a)') 'subtraction: t1-z2 ='
call mpwrite (6, ndp + 20, ndp, t1 - z2)
write (6, 1) d1 - dc2
call checkdc (t1 - z2, d1 - dc2, ndp)

write (6, '(a)') 'multiplication: z1*z2 ='
call mpwrite (6, ndp + 20, ndp, z1 * z2)
write (6, 1) dc1 * dc2
call checkdc (z1 * z2, dc1 * dc2, ndp)

write (6, '(a)') 'multiplication: z1*e2 ='
call mpwrite (6, ndp + 20, ndp, z1 * e2)
write (6, 1) dc1 * e2
call checkdc (z1 * e2, dc1 * e2, ndp)

write (6, '(a)') 'multiplication: e1*z2 ='
call mpwrite (6, ndp + 20, ndp, e1 * z2)
write (6, 1) e1 * dc2
call checkdc (e1 * z2, e1 * dc2, ndp)

write (6, '(a)') 'multiplication: z1*ec2 ='
call mpwrite (6, ndp + 20, ndp, z1 * ec2)
write (6, 1) dc1 * ec2
call checkdc (z1 * ec2, dc1 * ec2, ndp)

write (6, '(a)') 'multiplication: ec1*z2 ='
call mpwrite (6, ndp + 20, ndp, ec1 * z2)
write (6, 1) ec1 * dc2
call checkdc (ec1 * z2, ec1 * dc2, ndp)

write (6, '(a)') 'multiplication: z1*t2 ='
call mpwrite (6, ndp + 20, ndp, z1 * t2)
write (6, 1) dc1 * d2
call checkdc (z1 * t2, dc1 * d2, ndp)

write (6, '(a)') 'multiplication: t1*z2 ='
call mpwrite (6, ndp + 20, ndp, t1 * z2)
write (6, 1) d1 * dc2
call checkdc (t1 * z2, d1 * dc2, ndp)

write (6, '(a)') 'multiplication: z1*e2 ='
call mpwrite (6, ndp + 20, ndp, z1 * e2)
write (6, 1) dc1 * e2
call checkdc (z1 * e2, dc1 * e2, ndp)

write (6, '(a)') 'multiplication: e1*z2 ='
call mpwrite (6, ndp + 20, ndp, e1 * z2)
write (6, 1) e1 * dc2
call checkdc (e1 * z2, e1 * dc2, ndp)

write (6, '(a)') 'division: z1/z2 ='
call mpwrite (6, ndp + 20, ndp, z1 / z2)
write (6, 1) dc1 / dc2
call checkdc (z1 / z2, dc1 / dc2, ndp)

write (6, '(a)') 'division: z1/e2 ='
call mpwrite (6, ndp + 20, ndp, z1 / e2)
write (6, 1) dc1 / e2
call checkdc (z1 / e2, dc1 / e2, ndp)

write (6, '(a)') 'division: e1/z2 ='
call mpwrite (6, ndp + 20, ndp, e1 / z2)
write (6, 1) e1 / dc2
call checkdc (e1 / z2, e1 / dc2, ndp)

write (6, '(a)') 'division: z1/ec2 ='
call mpwrite (6, ndp + 20, ndp, z1 / ec2)
write (6, 1) dc1 / ec2
call checkdc (z1 / ec2, dc1 / ec2, ndp)

write (6, '(a)') 'division: ec1/z2 ='
call mpwrite (6, ndp + 20, ndp, ec1 / z2)
write (6, 1) ec1 / dc2
call checkdc (ec1 / z2, ec1 / dc2, ndp)

write (6, '(a)') 'division: z1/t2 ='
call mpwrite (6, ndp + 20, ndp, z1 / t2)
write (6, 1) dc1 / d2
call checkdc (z1 / t2, dc1 / d2, ndp)

write (6, '(a)') 'division: t1/z2 ='
call mpwrite (6, ndp + 20, ndp, t1 / z2)
write (6, 1) d1 / dc2
call checkdc (t1 / z2, d1 / dc2, ndp)

write (6, '(a)') 'exponentiation: z1**i1 ='
call mpwrite (6, ndp + 20, ndp, z1 ** i1)
write (6, 1) dc1 ** i1
call checkdc (z1 ** i1, dc1 ** i1, ndp)

write (6, '(a)') 'exponentiation: z1**z2 ='
call mpwrite (6, ndp + 20, ndp, z1 ** z2)
write (6, 1) dc1 ** dc2
call checkdc (z1 ** z2, dc1 ** dc2, ndp)

write (6, '(a)') 'exponentiation: t1**z2 ='
call mpwrite (6, ndp + 20, ndp, t1 ** z2)
write (6, 1) d1 ** dc2
call checkdc (t1 ** z2, d1 ** dc2, ndp)

write (6, '(a)') 'exponentiation: z1**t2 ='
call mpwrite (6, ndp + 20, ndp, z1 ** t2)
write (6, 1) dc1 ** d2
call checkdc (z1 ** t2, dc1 ** d2, ndp)

write (6, '(a)') 'equal test: z1 == z2, e1 == z2, z1 == e2, ec1 == z2, z1 == ec2'
l1 = z1 == z2; l2 = dc1 == dc2; l3 = e1 == z2; l4 = e1 == dc2; l5 = z1 == e2
l6 = ec1 == e2; l7 = ec1 == z2; l8 = ec1 == dc2; l9 = z1 == ec2; l10 = dc1 == ec2
write (6, '(10l4)') l1, l2, l3, l4, l5, l6, l7, l8, l9, l10
call checkl (l1, l2); call checkl (l3, l4); call checkl (l5, l6)
call checkl (l7, l8); call checkl (l9, l10)

write (6, '(a)') 'not-equal test: z1 /= z2, e1 /= z2, z1 /= e2, ec1 /= z2, z1 /= ec2'
l1 = z1 /= z2; l2 = dc1 /= dc2; l3 = e1 /= z2; l4 = e1 /= dc2; l5 = z1 /= e2
l6 = ec1 /= e2; l7 = ec1 /= z2; l8 = ec1 /= dc2; l9 = z1 /= ec2; l10 = dc1 /= ec2
write (6, '(10l4)') l1, l2, l3, l4, l5, l6, l7, l8, l9, l10
call checkl (l1, l2); call checkl (l3, l4); call checkl (l5, l6)
call checkl (l7, l8); call checkl (l9, l10)

write (6, '(a)') 'abs(z2) ='
call mpwrite (6, ndp + 20, ndp, abs (z2))
write (6, '(1p,2d25.15)') abs (dc2)
call checkdp (abs (z2), abs (dc2), ndp)

write (6, '(a)') 'aimag(z1) ='
call mpwrite (6, ndp + 20, ndp, aimag (z1))
write (6, '(1p,2d25.15)') aimag (dc1)
call checkdp (aimag (z1), aimag (dc1), ndp)

write (6, '(a)') 'conjg(z1) ='
call mpwrite (6, ndp + 20, ndp, conjg (z1))
write (6, '(1p,2d25.15)') conjg (dc1)
call checkdc (conjg (z1), conjg (dc1), ndp)

write (6, '(a)') 'cos(z2) ='
call mpwrite (6, ndp + 20, ndp, cos (z2))
write (6, '(1p,2d25.15)') cos (dc2)
call checkdc (cos (z2), cos (dc2), ndp)

write (6, '(a)') 'dcmplx(z1) ='
write (6, '(1p,2d25.15)') cmplx (dc1, kind=mprknd)
write (6, '(1p,2d25.15)') dc1
call checkdc (z1, dc1, ndp)

write (6, '(a)') 'exp(z1) ='
call mpwrite (6, ndp + 20, ndp, exp (z1))
write (6, '(1p,2d25.15)') exp (dc1)
call checkdc (exp (z1), exp (dc1), ndp)

write (6, '(a)') 'log(z1) ='
call mpwrite (6, ndp + 20, ndp, log (z1))
write (6, '(1p,2d25.15)') log (dc1)
call checkdc (log (z1), log (dc1), ndp)

write (6, '(a)') 'mpcmplx(t1,t2) ='
call mpwrite (6, ndp + 20, ndp, mpcmplx (t1, t2, nwds))
write (6, '(1p,2d25.15)') cmplx (d1, d2, mprknd)
call checkdc (mpcmplx (t1, t2, nwds), cmplx (d1, d2, mprknd), ndp)

write (6, '(a)') 'mpreal(z1) ='
call mpwrite (6, ndp + 20, ndp, mpreal (z1, nwds))
write (6, '(1p,2d25.15)') dble (dc1)
call checkdp (mpreal (z1, nwds), dble (dc1), ndp)

write (6, '(a)') 'sin(z2) ='
call mpwrite (6, ndp + 20, ndp, sin (z2))
write (6, '(1p,2d25.15)') sin (dc2)
call checkdc (sin (z2), sin (dc2), ndp)

write (6, '(a)') 'sqrt(z1) ='
call mpwrite (6, ndp + 20, ndp, sqrt (z1))
write (6, '(1p,2d25.15)') sqrt (dc1)
call checkdc (sqrt (z1), sqrt (dc1), ndp)

write (6, '(a)') 'sqrt(z2) ='
call mpwrite (6, ndp + 20, ndp, sqrt (z2))
write (6, '(1p,2d25.15)') sqrt (dc2)
call checkdc (sqrt (z2), sqrt (dc2), ndp)

stop
end

subroutine checkdp (t1, d1, ndp)
use mpmodule
implicit none
type (mp_real) t1
real (8) d1, dtol
parameter (dtol = 1d-14)
integer ndp

if (abs ((dble (t1) - d1) / d1) > dtol) then
  write (6, '(a)') 'error:'
  call mpwrite (6, ndp + 20, ndp, t1)
  write (6, '(1p,d30.16)') d1
  stop
endif

return
end

subroutine checkdc (z1, dc1, ndp)
use mpmodule
implicit none
type (mp_complex) z1
complex (8) dc1
real (8) dtol
parameter (dtol = 1d-14)
integer ndp

if (abs ((dcmplx (z1) - dc1) / dc1) > dtol) then
  write (6, '("error")')
  call mpwrite (6, ndp + 20, ndp, z1)
  write (6, '(1p,2d30.16)') dc1
  stop
endif

return
end

subroutine checkl (l1, l2)
logical l1, l2

if (l1 .neqv. l2) then
  write (6, '("error:",2l4)') l1, l2
  stop
endif

return
end

max_stars_repo_name: trcameron/FPML-Comp

max_issues_repo_name: trcameron/FPML-Comp

max_forks_repo_name: trcameron/FPML-Comp


#INCLUDE 'MR_H_ALIGN_PADDING.H'
!***********************************************************************************************************************************
! UNIT:
!
!  (MODULE)
!
! PURPOSE:
!
!   TO
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  MODULE MR_MOD_CALC_EQKD_PRO_XY

    USE MR_KINDS

    USE MR_DEF_CURVED_GEOS
    USE MR_DEF_CONSTS_N_REF_PARS
    USE MR_DEF_FIELD_VARS

    IMPLICIT NONE

    PRIVATE

    PUBLIC :: MR_CALC_EQKD_PRO_XY

    REAL   (FDRD_KIND) , ALLOCATABLE , DIMENSION(:,:,:    ) :: UVT , UUT , VVT

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

  CONTAINS

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!   TO
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_CALC_EQKD_PRO_XY( NI , NJ , NK , EQKD_PRO_XY )

    USE MR_MOD_OPERATOR_UV
    USE MR_MOD_OPERATOR_XUV

    USE MR_MOD_CALC_QADV_N_QDIF_XY
    USE MR_MOD_INTERP_XY

    IMPLICIT NONE

    INTEGER(IJID_KIND) , INTENT(IN ) :: NI , NJ
    INTEGER(KKID_KIND) , INTENT(IN ) :: NK

    REAL   (FDRD_KIND) , INTENT(OUT) , DIMENSION(1:NI1(NI,FDRD_KIND),1:NJ,1:NK) :: EQKD_PRO_XY

    REAL   (FDRD_KIND) , DIMENSION(1:NI1(NI,FDRD_KIND),1:NJ,1:2,1:2) :: HYD_T_XY , HYD_T_XY_EXT
    REAL   (FDRD_KIND) , DIMENSION(0:NI0(NI,FDRD_KIND),1:NJ,1:2    ) :: HYD_T_XY_U
    REAL   (FDRD_KIND) , DIMENSION(1:NI1(NI,FDRD_KIND),0:NJ,1:2    ) :: HYD_T_XY_V

    INTEGER(IJID_KIND) :: I , J
    INTEGER(KKID_KIND) :: K

    DO K = 1 , NK

      ALLOCATE( UVT(1:NI1(NI,FDRD_KIND),1:NJ,1:2) )

        UVT = ( JUV .MRUVTFM. UV(:,:,1:2, K ) )

        ALLOCATE( VVT(1:NI1(NI,FDRD_KIND),0:NJ,1:2) )
          CALL MR_INTERP_XY_UV_V_TO_GET_U_AT_V( NI , NJ , UV(:,:,1:2, K ) , VVT(:,:,1) )
         !CALL MR_INTERP_XY_UV_V_BY_LINEAR( NI , NJ , UV(:,:,1:2, K ) , VVT(:,:,2) )
          VVT(:,:,2) = V(:,:, K )

          VVT = ( JVV .MRUVTFM. VVT )

          CALL MR_CALC_T_XY_UV_U( NI , NJ , UVT , VVT , RB , U(:,:, K ) , EKXY , VXYU(:,:, K ) , HYD_T_XY_U )

        DEALLOCATE( VVT )

        ALLOCATE( UUT(0:NI0(NI,FDRD_KIND),1:NJ,1:2) )
          CALL MR_INTERP_XY_UV_U_TO_GET_V_AT_U( NI , NJ , UV(:,:,1:2, K ) , UUT(:,:,2) )
         !CALL MR_INTERP_XY_UV_U_BY_LINEAR( NI , NJ , UV(:,:,1:2, K ) , UUT(:,:,1) )
          UUT(:,:,1) = U(:,:, K )

          UUT = ( JUU .MRUVTFM. UUT )

          CALL MR_CALC_T_XY_UV_V( NI , NJ , UVT , UUT , RB , V(:,:, K ) , EKXY , VXYV(:,:, K ) , HYD_T_XY_V )

        DEALLOCATE( UUT )

      DEALLOCATE( UVT )

    ! CALCULATE MECHANICAL PRODUCTION
      CALL MR_INV_INTERP_XY_UV_FROM_UU( NI , NJ , HYD_T_XY_U , HYD_T_XY(:,:,1:2,1) )
      CALL MR_INV_INTERP_XY_UV_FROM_VV( NI , NJ , HYD_T_XY_V , HYD_T_XY(:,:,1:2,2) )
      HYD_T_XY = ( HYD_T_XY .MRXUVMAT. ( .MRXUVTPS. JUV ) )
      HYD_T_XY_EXT = ( HYD_T_XY + ( .MRXUVTPS. HYD_T_XY ) )
      DO J = 1 , NJ
       !DIR$ VECTOR ALIGNED
        DO I = 1 , NI
          EQKD_PRO_XY( I , J , K ) = SUM( HYD_T_XY_EXT( I , J ,1:2,1:2) * HYD_T_XY( I , J ,1:2,1:2) ) / ( EKXY * VZWW( I , J , K ) )
        END DO
      END DO

    END DO

  END SUBROUTINE MR_CALC_EQKD_PRO_XY

  END MODULE MR_MOD_CALC_EQKD_PRO_XY


text:
!===============================================================================
! One of Andy Nowacki's Fortran utility modules for dealing with seismic
! anisotropy and other problems.
!
! Andy Nowacki <andy.nowacki@bristol.ac.uk>
!
! See the file LICENCE for licence details.
!===============================================================================
program azi_prog
! Program to take two sets of spherical co-ordinates
! and write out an azimuth from the first to the second.

use spherical_geometry, only: azimuth

implicit none

integer, parameter :: rs = selected_real_kind(15,307)
real(rs) :: lon1,lat1,lon2,lat2
integer :: iostatus

if (command_argument_count() /= 0) then
	write(0,'(a)') 'Usage: azi < [stlon stlat evlon evlat]'
	stop
endif

iostatus = 0
do while (iostatus == 0)
	read(*,*,iostat=iostatus) lon1,lat1,lon2,lat2
	if (iostatus < 0) exit
	if (iostatus > 0) then
		write(0,'(a)') 'azi: Error: problem reading coordinates from stdin.'
		stop
	endif
	write(*,'(f10.6)') azimuth(lon1,lat1,lon2,lat2,degrees=.true.)
enddo

end program azi_prog

max_stars_repo_name: JackWalpole/seismo-fortran-fork

max_issues_repo_name: JackWalpole/seismo-fortran-fork

max_forks_repo_name: JackWalpole/seismo-fortran-fork


text:
	SUBROUTINE MN_3CHC ( report, istart, iprms, rdata, iret )
C************************************************************************
C* MN_3CHC								*
C*									*
C* This subroutine decodes a line of the NGM MOS report that has 1-, 2-,*
C* or 3- letter codes. After the strings are parsed, they are converted	*
C* to reals. 								*
C*									*
C* MN_3CHC ( REPORT, ISTART, IPRMS, RDATA, IRET )			*
C*									*
C* Input parameters:							*
C*	REPORT		CHAR*		NGM MOS report			*
C*	ISTART		INTEGER		Index of line beginnings	*
C*	IPRMS		INTEGER		Position of parameter in list	*
C*									*
C* Output parameters:							*
C*	RDATA		REAL		Forecast data			*
C*	(FCSTTM,MMPARM)							*
C*	IRET		INTEGER		Return Code			*
C**									*
C* Log:									*
C* D. Keiser/GSC	 2/96						*
C* D. Keiser/GSC	 5/96	Declare FCSTTM an integer		*
C* F. J. Yen/NCEP	11/98	Converted for restructuring of DC_NMDC	*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INTEGER		FCSTTM
	PARAMETER	( FCSTTM = 19 )
C*
	CHARACTER*(*)	report
	REAL		rdata(FCSTTM,MMPARM)
C
	CHARACTER	tstr*3
C------------------------------------------------------------------------
	iret = 0
C
C*	Parse the specified line of the report.
C
	IF ( istart .ne. 0 ) THEN
	    ibeg = istart
	    iend = INDEX ( report(istart:), CHCR ) + istart
	    incr = 2
	    DO i = 1, FCSTTM
		IF ( ibeg .lt. iend ) THEN
		    tstr = report (ibeg:ibeg+incr)
C
C*		    Convert each string to a real.
C
C*		    If the CLDS line, convert to SAO code.
C
		    IF ( iprms .eq. 5 ) THEN
		    	IF ( tstr .eq. ' CL' ) tstr = 'CLR'
		    	IF ( tstr .eq. ' SC' ) tstr = 'SCT'
		    	IF ( tstr .eq. ' BK' ) tstr = 'BKN'
		    	IF ( tstr .eq. ' OV' ) tstr = 'OVC'
			rdata(i,iprms) = PT_CCNM ( tstr )
		    END IF
C
C*		    PTYPE line.
C
		    IF ( iprms .eq. 16 ) THEN
			IF ( tstr .eq. '  R' ) rdata(i,iprms) = 0.0
			IF ( tstr .eq. '  S' ) rdata(i,iprms) = 1.0
			IF ( tstr .eq. '  Z' ) rdata(i,iprms) = 2.0
		    END IF
C
C*		    OBVIS line.
C
		    IF ( iprms .eq. 23 ) THEN
			IF ( tstr .eq. '  N' ) rdata(i,iprms) = 0.0
			IF ( tstr .eq. '  F' ) rdata(i,iprms) = 1.0
			IF ( tstr .eq. '  H' ) rdata(i,iprms) = 2.0
		    END IF
		    ibeg = ibeg + incr + 1
		END IF
	    END DO
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
great_rsd(v1,hh1).
great_rsd(c1,x1).
great_rsd(c1,jj1).
great_rsd(v1,g1).
great_rsd(p1,w1).
great_rsd(u1,g1).
great_rsd(u1,ee1).
great_rsd(e1,p1).
great_rsd(f1,ee1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
program test_debug_print
#include "ccs_macros.inc"
  use utils, only : str, debug_print
  use kinds, only : ccs_int, ccs_real
  use mpi
  implicit none

  integer :: ierror
  character(len=:), allocatable :: msg

  call mpi_init(ierror)
  msg = "The numbers are " // str(24_ccs_int) // " and " // &
        str(42_ccs_int, "(I3)") // " and " // &
        str(6.28_ccs_real, "(F5.3)") // " and " // str(3.14_ccs_real)
  call dprint(msg)
  call mpi_finalize(ierror)
end program test_debug_print

max_stars_repo_name: asimovpp/asimov-ccs

max_issues_repo_name: asimovpp/asimov-ccs

max_forks_repo_name: asimovpp/asimov-ccs


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

!>    @brief write a Tecplot header to file "fileid"
!>    @param[in] fileid file handler
!>    @param[in] ni grid size in I0 direction
!>    @param[in] nj grid size in J0 direction
!>    @param[in] nk grid size in K0 direction
!>    @param[in] feld_name array/component name
!>    @param[in] ismpl local sample index
      SUBROUTINE head_tecpl(fileid,ni,nj,nk,feld_name,ismpl)
        use arrays
        use mod_inverse
        use mod_genrl
        use mod_genrlc
        IMPLICIT NONE
        integer :: ismpl
        integer :: i, j
        CHARACTER*25 zeilen_form, zz_form
        INTEGER, PARAMETER :: zbuff_size=16768
        CHARACTER*(zbuff_size) zeilen_buff
        CHARACTER*4 feld_name
        CHARACTER*16 pname
        INTEGER l10, maxp, sollp, fileid, ni, nj, nk

        CALL chln(project,i,j)
        WRITE(fileid,'(3A)') 'title = "', project(i:j), '"'

        maxp = 1
        sollp = 4*15 + 9 + 19*mpara
        IF (zbuff_size<sollp) THEN
          WRITE(*,*) &
            'error, "zeilen_buff" to small in "head_tecpl" !!!'
          STOP
        END IF

        l10 = int(log(dble(sollp))/log(10.0D0)) + 1
        zz_form = '(A13,I' // achar(l10+48) // ',A4)'
        WRITE(zeilen_form,zz_form) '(A', sollp, ')'

        zeilen_buff(maxp:maxp+15) = 'variables = "x"'
        maxp = maxp + 15
        zeilen_buff(maxp:maxp+15) = ',           "y"'
        maxp = maxp + 15
        zeilen_buff(maxp:maxp+15) = ',           "z"'
        maxp = maxp + 15
        zeilen_buff(maxp:maxp+9) = ',"uindex"'
        maxp = maxp + 9
        zeilen_buff(maxp:maxp+15) = ',        "' // feld_name // '"'
        maxp = maxp + 15

100     FORMAT (1A2,1A16,1A1)
        DO i = 1, mpara
          CALL param_name(i,pname,ismpl)
          WRITE(zeilen_buff(maxp:maxp+19),100) ',"', pname, '"'
          maxp = maxp + 19
        END DO

        IF ((maxp-1)/=sollp) THEN
          WRITE(*,*) 'error, "maxP"<>"sollP" in "head_tecpl" !!!'
          STOP
        END IF

        WRITE(fileid,'(A)') zeilen_buff(1:maxp-1)

        WRITE(fileid,'(3(a,i5),a)') 'zone i=', ni, ', j=', nj, &
          ', k=', nk, ', f=point'

        RETURN
      END


subroutine iter_profcv2(btop,bzd,bcf,bsfc,zKuL,zKaL,dr,n1d,eps,imu,&
     dn1d,dm1d,rrate1d,zKuC,zKaSim,epst,piaKu,piaKa,itype,dncv,dnp,&
     dzdn,piaSRTKu,relPIASRTKu,&
     rrate1d_sub, dn_sub, dm_sub, zkuc_sub, piahb_sub, piaKa_sub,zetaS,&
     piaKuS,piaKaS)
  use tables2
  use tableP2
  use ran_mod
  implicit none
  integer :: btop, bzd, bcf, bsfc, n1d, imu
  real :: zKuL(n1d), dr, zKaL(n1d), dnp(n1d)
  !real :: dmCoeff(2)=(/0.027773772993636318,-0.6624076086959071/)
  real,intent(out) :: dm1d(n1d)
  real,intent(out) :: epst,piaKu,piaKa
  real,intent(out) :: dn1d(n1d), rrate1d(n1d), zKuC(n1d), zKaSim(n1d), dzdn(n1d,n1d)
  real :: dnCoeff_new(2)
  real :: dzKa(n1d), dns(n1d)
  real :: eps, dm, pia, dncv
  integer :: it, ik, k, n1, n1H, itype
  real :: ztrueS,ztrueH,ztrue,f, attKuH, attKaH, attKuS, attKaS
  real :: ftran, probs(31), sumprob,  zeta1d(n1d,31), q
  real,intent(out) :: rrate1d_sub(n1d,31), dn_sub(n1d,31), dm_sub(n1d,31), zkuc_sub(n1d,31), piahb_sub(31)
  real,intent(out) :: piaKa_sub(31),zetaS(31)
  real,intent(out) :: piaKuS, piaKaS
  real :: beta, piamax, attKu, attKa, dn, dni
  real :: dnCoeff(2)=(/-0.00570364,  0.13319214/)
  real :: dmCoeff(3)= (/-2.56885571e-04,  7.18909743e-02, -1.60889523e+00/)
  !(/-4.57599678e-04,  7.852247958e-02, -1.78316499e+00/)
  !(/-4.57599678e-04,  8.52247958e-02, -1.78316499e+00/)
  real :: rn, dm_old!, dm_sub(n1d,31)
  real :: zka1,zka2,attka2
  real :: dt1, dt2
  real :: start, finish1, finish2
  integer :: kk, imc, isub
  real :: dmm
  real :: piaSRTKu, dzdn_sub2(n1d,n1d,31), dzdn_sub1(n1d,n1d,31), zkasim_sub(n1d,31)
  integer :: relPIASRTKu
  real :: ddn
  piaKu=0
  !dmCoeff=array([ 0.02893781, -0.69481455])
  dzdn=0.0
  !if(itype.eq.1) then
  !end if
  zKuC=zKuL
  dns=0
  dZKa=0
  rrate1d=0
  call cpu_time(start)
  !print*, n1d, bcf,bzd
  if(bzd<0) bzd=bsfc
  rrate1d=0
  rrate1d_sub=0
  do it=1,2
     piaKu=0.0
     piaKa=0.0
     do k=btop,min(bcf,bzd-1)
        if (zKuC(k+1).gt.10) then
           if (k.lt.bzd) then
              ztrue=zKuC(k+1)
           end if
           if(ztrue.lt.40) then
              ztrueS=ztrue
              ztrueH=0.0
           else
              ztrueS=40
              ztrueH=10*log10(10**(0.1*ztrue)-10**(0.1*39.999))
           endif
           !n1,n2=bisection(zKuSJ,ztrueS-10*dns[k+1])
           call bisection2(zKuSJ,nbinS2,ztrueS-10*dns(k+1), n1)
           attKuS=att13TableS2(n1,imu)*10**dns(k+1)
           attKaS=att35TableS2(n1,imu)*10**dns(k+1)
           !n1H,n2H=bisection(zKuHJ,ztrueH)
           !print*, ztrueS,n1,attKuS
           call bisection2(zKuSJ,nbinS2,ztrueH, n1H)
           attKuH=att13TableH(n1H,imu)
           attKaH=att35TableH(n1H,imu)
           piaKu=piaKu+(attKuS+attKuH)*dr
           piaKa=piaKa+(attKaS+attKaH)*dr
           zKuC(k+1)=zKuL(k+1)+piaKu
           zKaSim(k+1)=10*log10(10**(0.1*z35TableS2(n1,imu)+dns(k+1))+&
                10**(0.1*z35TableH(n1H,imu)))-piaKa
           if(zKaL(k+1).gt.10) then
              dZKa(k+1)=zKaSim(k+1)-zKaL(k+1)
           end if
           piaKu=piaKu+(attKuS+attKuH)*dr
           piaKa=piaKa+(attKaS+attKaH)*dr
           dm1d(k+1)=0.5*(d013TableS2(n1,imu)+d013TableH(n1H,imu))
           dn1d(k+1)=dns(k+1)
           rrate1d(k+1)=pr13TableS2(n1,imu)*10**dns(k+1)+pr13TableH(n1H,imu)
        end if
     end do
     if(it.eq.1) then
        do k=btop,min(bcf,bzd-1)
           dns(k+1)=dns(k+1)-0.6*dZKa(k+1)
           if(dns(k+1).gt.0.975) dns(k+1)=0.975
           if(dns(k+1).lt.-1.0) dns(k+1)=-1.0
        end do
     end if
  end do
  !print*, piaKu,piaKa,btop,bcf, zkuL(btop+1:bcf+1)
  !return
  piaKuS=piaKu+0.0
  piaKaS=piaKa+0.0

  zeta1d=0
  beta=0.71

  piamax=54-zKuL(bcf+1)
  if(piamax<0) piamax=56-zKuL(bcf+1)
  if(piamax<0) piamax=4
  q=0.2*log(10.0)
  zKuC(bzd:bcf+1)=zKuL(bzd:bcf+1)+piaKuS
  epst=eps
  attKu=0.0
  attKa=0.0
  f=1.0
  zeta1d=0
  sumprob=0
  piaKu=0.0
  rrate1d(bzd:bcf+1)=0.0
  dm1d(bzd+1:bcf+1)=0.0
  dn1d(bzd+1:bcf+1)=0.0
  if(bzd.le.bcf) then
     rrate1d_sub(bzd:bcf+1,:)=0
     dm_sub(bzd:bcf+1,:)=0
  endif
  dn_sub=0
  zkasim_sub=0
  dzdn_sub2=0
  dzdn_sub1=0
  zkasim(bzd:bcf+1)=0
  dzdn(bzd:bcf+1,bzd:bcf+1)=0
  call cpu_time(finish1)
  ddn=0.025
!!$OMP PARALLEL DO  PRIVATE(k,dn,n1,attKu,attKa,zKa1,zKa2)
  do isub=1,31
     zetaS(isub)=0
     do k=bzd,bcf
        if (zKuL(k+1).gt.10) then
           call bisection2(zKudN,nbins,zKuL(k+1)+piaKuS-10*(isub-16)*ddn-10*dncv, n1)
           dn_sub(k+1,isub)=(zKudN(n1)-zKuSJ(n1))/10.0+(isub-16)*ddn+dncv
           dn=dn_sub(k+1,isub)
           attKu=att13Table(n1,imu)*10**dn
           attKa=att35Table(n1,imu)*10**dn
           zetaS(isub)=zetaS(isub)+attKu*dr
           if(isub.eq.1.and.(zetaS(isub)*q*beta).gt.0.99) then
              print*, zKuL(k+1)+piaKuS,ddn,dncv,zetaS(isub), piamax
           end if
           zeta1d(k+1,isub)=zetaS(isub)
        else
           attKu=0.0
           attKa=0.0
        endif
     end do
     if(q*beta*zetaS(1).gt.0.995) then
        zetaS(isub)=(1-10**(-piamax*0.1*beta))/(q*beta)
     end if
     if (q*beta*zetaS(isub).lt.0.995) then
        piaKa_sub(isub)=piaKas
        piaHB_sub(isub)=-10/beta*log10(1-q*beta*zetaS(isub))+attKu*(bsfc-bcf)*2*dr
        do k=bzd,bcf
           if(zKuL(k+1).gt.10) then
              zKuC_sub(k+1,isub)=zKuL(k+1)-10/beta*log10(1-q*beta*zeta1d(k+1,isub))
              n1=int((zKuC_sub(k+1,isub)-zmin-10*dn_sub(k+1,isub))/dzbin)
              if(n1.lt.1) n1=1
              if(n1.gt.289) n1=289
              !dn=
              dn=dn_sub(k+1,isub)

              attKa=att35Table(n1,imu)*10**dn
              piaKa_sub(isub)=piaKa_sub(isub)+attKa*dr
              if(n1.lt.1) n1=1
              if(n1.gt.nbins) n1=nbins
              dm_sub(k+1,isub)=d013Table(n1,imu)
              rrate1d_sub(k+1,isub)=pr13Table(n1,imu)*10**dn_sub(k+1,isub)
              zKa1=z35Table(n1,imu)+10*dn
              zkasim_sub(k+1,isub)= zKa1-piaKa_sub(isub)
              piaKa_sub(isub)=piaKa_sub(isub)+attKa*dr
              if(n1.lt.nbins-4) then
                 zKa2=z35Table(n1+4,imu)+10*(dn-0.1)
                 attKa2=att35Table(n1+4,imu)*10**(dn-0.1)
                 dzdn_sub2(k+1,k+1,isub)=dzdn_sub2(k+1,k+1,isub)+(zka2-attKa2*dr)-((zka1-attKa*dr))
                 dzdn_sub1(k+1,k+1,isub)=dzdn_sub1(k+1,k+1,isub)+zka1-attKa*dr
                 do kk=k+1,bcf
                    dzdn_sub2(kk+1,k+1,isub)=dzdn_sub2(kk+1,k+1,isub)-(attKa2*dr-attKa*dr)
                    dzdn_sub1(kk+1,k+1,isub)=dzdn_sub1(kk+1,k+1,isub)-attKa*dr
                 end do
              end if
           else
              attKa=0
              rrate1d_sub(k+1,isub)=0.0
           end if
        end do
        piaKa_sub(isub)=piaKa_sub(isub)+attKa*(bsfc-bcf)*2*dr
        probs(isub)=exp(-((isub-16)*ddn)**2/0.5**2)
        if(rrate1d_sub(bcf+1,isub).gt.150) probs(isub)=probs(isub)*0.9
     else
        probs(isub)=0.0
     end if
  end do
  !print*,zetaS,bzd,bcf
  !!$OMP END PARALLEL DO
  piaKu=0.0
  piaKa=0.0
  sumprob=0.0
  do isub=1,31
     do k=bzd,bcf
        if(zkuL(k+1).gt.10) then
           rrate1d(k+1)= rrate1d(k+1)+probs(isub)*rrate1d_sub(k+1,isub)
           dm1d(k+1)=dm1d(k+1)+probs(isub)*dm_sub(k+1,isub)
           dn1d(k+1)=dn1d(k+1)+probs(isub)*dn_sub(k+1,isub)
           zkasim(k+1)=zkasim(k+1)+10**(0.1*zkasim_sub(k+1,isub))*probs(isub)
           dzdn(k+1,k+1)=dzdn(k+1,k+1)+10**(0.1*zkasim_sub(k+1,isub)+0.1*dzdn_sub2(k+1,k+1,isub))*probs(isub)
           do kk=k+1,bcf
              if(zkuL(kk+1).gt.10) then
                 dzdn(kk+1,k+1)=dzdn(kk+1,k+1)+10**(0.1*zkasim_sub(kk+1,isub)+0.1*dzdn_sub2(kk+1,k+1,isub))*probs(isub)
              end if
           end do
        end if
     end do
     piaKu=piaKu+probs(isub)*piaHB_sub(isub)
     piaKa=piaKa+probs(isub)*piaKa_sub(isub)
     sumprob=sumprob+probs(isub)
  end do

  piaKu=piaKu/(sumprob+1e-10)+piaKuS
  piaKa=piaKa/(sumprob+1e-10)
  do k=bzd,bcf
     if(zkuL(k+1).gt.10) then
        rrate1d(k+1)= rrate1d(k+1)/(sumprob+1e-10)
        dm1d(k+1)=dm1d(k+1)/(sumprob+1e-10)
        dn1d(k+1)=dn1d(k+1)/(sumprob+1e-10)
        zkasim(k+1)=log10(zkasim(k+1)/(sumprob+1e-10))*10.0
        dzdn(k+1,k+1)=(log10(dzdn(k+1,k+1)/(sumprob+1e-10))*10.0-zkasim(k+1))/(-0.1)
        do kk=k+1,bcf
           if(zkuL(kk+1).gt.10) then
              dzdn(kk+1,k+1)=(log10(dzdn(kk+1,k+1)/(sumprob+1e-10))*10.0-log10(zkasim(kk+1)/(sumprob+1e-9))*10)/(-0.1)
           end if
        end do
     end if
  end do
  call cpu_time(finish2)
  dt1=finish1-start
  dt2=finish2-start
  if(rrate1d(bcf+1).gt.150) then
     print*, rrate1d_sub(bcf+1,1:31:3)
     print*, 'probs=',probs(1:31:3)
     print*, 'zetaS=',zetaS(1:31:3)
  end if
end subroutine iter_profcv2

subroutine convRetf90(rrate,dmOut,dm_sub,rrate_sub,bzd,bcf,piaka_sub,piahb_sub,&
     piaSRTKu,dpiaSRT,relFlagKu,relFlag,zetaS)
  implicit none
  integer :: bzd,bcf
  real,intent(out) :: rrate(176), dmOut(176)
  real :: dm_sub(176,31), rrate_sub(176,31), piaka_sub(31),&
       piahb_sub(31), zetaS(31)
  real :: probS,prob1,q,beta
  real :: dpiaSRT, piaSRTKu
  integer  :: relFlag,relFlagKu
  integer :: i, k
  probS=0
  q=0.2*log(10.0)
  beta=0.71
  !print*, bzd,bcf
  !print*, zetaS
  if(bzd.le.bcf) then
     rrate(bzd+1:bcf+1)=0.0
  endif
  !print*, rrate_sub(bcf+1,:), rrate(bcf+1)
  !print*, 'bzd',bzd, bcf
  do i=1,31
     if(q*beta*zetaS(i)<0.995.and.rrate_sub(bcf,i).lt.250.0.and.&
          rrate_sub(bcf+1,i).lt.250.0) then
        if (relFlag==1 .and. dpiaSRT<20) then
           prob1=exp(-0.5*(piaka_sub(i)-piahb_sub(i)-dpiaSRT)**2/2.0)
        else
           prob1=1.0
        endif
        if(relFlagKu==1 .and. piaSRTKu>5.0) then
           prob1=prob1*exp(-0.5*(piahb_sub(i)-piaSRTKu)**2/4.0)
        end if
        prob1=prob1+1e-9
        do k=bzd+1,bcf+1
           rrate(k)=rrate(k)+prob1*rrate_sub(k,i)
           dmOut(k)=dmOut(k)+prob1*dm_sub(k,i)
        end do
        probS=probS+prob1
     endif
  end do
  !print*,rrate(bcf+1), probS
  if(bzd.le.bcf) then
     rrate(bzd+1:bcf+1)=rrate(bzd+1:bcf+1)/probS
     dmOut(bzd+1:bcf+1)=dmOut(bzd+1:bcf+1)/probS
  endif
  !print*,rrate(bcf+1)

end subroutine convRetf90


subroutine cf_synop_check(vor_in, vor_part, &
                        & n_part, nx, ny,   &
                        & lon, lat, mtype_part)

  use types, only: wp
  use constants, only: pi, ra, rkilo, pmax, deg2rad
  use params, only: proj, length_cf, area_ec, theta_cf_min, theta_cf_max, curv_cf, det_coeff_cf
  use utils, only: cosd, sind, great_circle

  implicit none

  real   (wp), intent(in)  :: vor_in  (0:nx, 0:ny)
  integer    , intent(in)  :: vor_part(0:nx, 0:ny)
  integer    , intent(in)  :: n_part
  integer    , intent(in)  :: nx
  integer    , intent(in)  :: ny
  real   (wp), intent(in)  :: lon     (0:nx      )
  real   (wp), intent(in)  :: lat     (      0:ny)
  integer    , intent(out) :: mtype_part
  ! Local variables
  real   (wp)              :: lonin, latin
  integer                  :: pnum !----cold front
  integer    , parameter   :: pnum_thresh = 6
  integer                  :: j_n, j_s, i_n, i_s
  logical                  :: flag_one
  real   (wp)              :: d_cf, theta_d_cf
  real   (wp)              :: one(pmax, 1:2)
  integer                  :: one_num
  real   (wp)              :: a, b, c, k, r2
  real   (wp)              :: size_vor ! synoptic
  integer                  :: i, j

  lonin = lon(1)-lon(0)
  latin = lat(1)-lat(0)

  theta_d_cf = 0. ! to avoid warnings
  d_cf = 0. ! just in case

  i_n = 0
  i_s = nx
  j_n = 0
  j_s = ny

  size_vor = 0.

  one_num = 0
  one(:, :) = 0

  do j = 1, ny-1
    do i = 1, nx-1
      if (vor_part(i, j) == n_part) then
        if (proj == 1) then
          size_vor = size_vor &
                   & + (lonin * deg2rad * ra / rkilo * cosd(lat(j)) &
                   &    * latin * deg2rad * ra / rkilo)
        elseif (proj == 2) then
          size_vor = size_vor + lonin * latin * rkilo**(-2)
        endif

        call median_check(vor_in(i-1:i+1, j-1:j+1), flag_one, pnum)

        if (flag_one) then
          if (pnum >= pnum_thresh) then
            one_num = one_num + 1
            one(one_num, 1) = lon(i)
            one(one_num, 2) = lat(j)
          endif
        endif

        if (j < j_s) then
          j_s = j
          i_s = i
        endif
        if (j >= j_n) then
          j_n = j
          i_n = i
        endif
      endif
    enddo
  enddo

  !if (dbg) then
  !    print*, 'cf_synop_check(99):', 'i_s', i_s, 'i_n', i_n
  !endif

  if (proj == 1) then
    d_cf = great_circle(lon(i_n), lon(i_s), lat(j_n), lat(j_s), ra)
    theta_d_cf = acos((ra * (lat(j_n) - lat(j_s)) * deg2rad) / d_cf)
  elseif(proj==2)then
    d_cf = sqrt((lon(i_s) - lon(i_n))**2 + (lat(j_s) - lat(j_n))**2) * 1.0e-3
    theta_d_cf = (lat(j_s) - lat(j_n)) * 1.0e-3 / d_cf
  endif

  if (d_cf >= length_cf * rkilo) then
    if (lon(i_n) > lon(i_s) .and. theta_d_cf <= theta_cf_max * deg2rad) then
      call quadric_fit(one(1:one_num, 1:2), one_num, a, b, c, k, r2)

      if(r2 >= det_coeff_cf .and. k <= curv_cf) then
        mtype_part=1
      endif

    elseif (theta_d_cf <= theta_cf_min * deg2rad) then
      call quadric_fit(one(1:one_num, 1:2), one_num, a, b, c, k, r2)

      if(r2 >= det_coeff_cf .and. k<= curv_cf) then
        mtype_part=1
      endif
    endif
  endif

!--------------- Check synoptic low ----------------
  if (size_vor >= area_ec) then
    mtype_part = mtype_part + 2
  endif
end subroutine cf_synop_check

subroutine median_check(val, flag, pnum)

  use types, only: wp
  use constants, only: mx, my

  implicit none

  real   (wp), intent (in)  :: val(-1:1, -1:1)
  logical    , intent (out) :: flag
  integer    , intent (out) :: pnum
  ! Local variables
  integer                   :: m
  integer                   :: mmax

  flag = .false.
  pnum = 0
  mmax = size(mx)

  do m = 1, mmax
    if (val (0, 0) > val(mx(m), my(m))) then
      pnum = pnum + 1
    endif
  enddo

  if (pnum >= 5) flag = .true.
end subroutine median_check

subroutine quadric_fit(one, one_num, a, b, c, k, r2)

  use types, only: wp

  implicit none

  integer    , intent (in) :: one_num
  real   (wp), intent (in) :: one(1:one_num, 1:2)
  real   (wp), intent (out):: a, b, c, k, r2
  ! Local variables
  integer                  :: i
  real   (wp)              :: sx, sx2, sx3, sx4, sx2y, sxy, sy, s1
  real   (wp)              :: xbar, ybar
  real   (wp)              :: x  (1:one_num)
  real   (wp)              :: y  (1:one_num)
  real   (wp)              :: fit(1:one_num)

  sx = 0.
  sx2 = 0.
  sx3 = 0.
  sx4 = 0.
  sx2y = 0.
  sxy = 0.
  sy = 0.
  s1 = 0.

  xbar = 0.
  ybar = 0.

  do i = 1, one_num
    xbar = xbar + one(i, 2)
    ybar = ybar + one(i, 1)
  enddo

  xbar = xbar / one_num
  ybar = ybar / one_num

  do i = 1, one_num
    x(i) = one(i, 2) - xbar
    y(i) = one(i, 1) - ybar
  enddo

  do i = 1, one_num
    sx4 = sx4 + x(i)**4
    sx3 = sx3 + x(i)**3
    sx2 = sx2 + x(i)**2
    sx = sx + x(i)**1
    sx2y = sx2y + x(i)**2 * y(i)
    sxy = sxy + x(i) * y(i)
    sy = sy + y(i)
    s1 = s1 + 1.0
  enddo

  a = (s1*sx2*sx2y-sx*sx*sx2y+sx*sx2*sxy-s1*sx3*sxy+sx*sx3*sy-sx2*sx2*sy)&
       &/(2*sx*sx2*sx3+s1*sx2*sx4-sx*sx*sx4-s1*sx3*sx3-sx2*sx2*sx2)

  b = (sx*sx2*sx2y-s1*sx3*sx2y+s1*sx4*sxy-sx2*sx2*sxy+sx2*sx3*sy-sx*sx4*sy)&
       &/(2*sx*sx2*sx3+s1*sx2*sx4-sx*sx*sx4-s1*sx3*sx3-sx2*sx2*sx2)

  c = (-sx2*sx2*sx2y+sx*sx3*sx2y-sx*sx4*sxy+sx2*sx3*sxy-sx3*sx3*sy+sx2*sx4*sy)&
       &/(2*sx*sx2*sx3+s1*sx2*sx4-sx*sx*sx4-s1*sx3*sx3-sx2*sx2*sx2)

  do i = 1, one_num
    fit(i) = a * x(i)**2 + b * x(i) + c
  enddo

  call correlation(y(1:one_num), fit(1:one_num), one_num, r2)

  k = 0.
  do i = 1, one_num
    k = k + a / (1 + (2 * a * x(i) + b)**2.) ** 1.5
  enddo

  k = k / one_num
end subroutine quadric_fit

subroutine correlation(x, y, n, r2)

  use types, only: wp

  implicit none

  integer    , intent (in) :: n
  real   (wp), intent (in) :: x(n)
  real   (wp), intent (in) :: y(n)
  real   (wp), intent (out):: r2
  ! Local variables
  real   (wp)              :: xm, ym, sigx, sigy, cor
  real   (wp)              :: sx, sy, sxy, sx2, sy2
  integer                  :: i

  sx = 0; sy = 0; sxy = 0; sx2 = 0; sy2 = 0

  if (n > 2) then
    do i = 1, n
      sx = sx + x(i)
      sy = sy + y(i)
    enddo
    xm = sx / n
    ym = sy / n
    do i = 1, n
      sx2 = sx2 + (x(i) - xm) * (x(i) - xm)
      sy2 = sy2 + (y(i) - ym) * (y(i) - ym)
      sxy = sxy + (x(i) - xm) * (y(i) - ym)
    enddo

    r2 = sxy * sxy / (sx2 * sy2)

    sigx = sqrt(sx2/n)
    sigy = sqrt(sy2/n)
    if (sigx > 0 .and. sigy > 0) then
      cor = (sxy / n) / (sigx * sigy)
    endif
  endif
end subroutine correlation