text:
! This file is part of nlopt-f.
! SPDX-Identifier: Apache-2.0 OR MIT
!
! Licensed under either of Apache License, Version 2.0 or MIT license
! at your option; you may not use this file except in compliance with
! the License.
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.

!> Implementation of the callback functionality for the NLopt library.
!>
!> This module provides ISO-C compatible callback functions which are
!> passed to the NLopt library and redirect the computation to a regular
!> Fortran procedure, which in turn can be more generic since it does
!> have to adhere to the constraints placed on ISO-C compatible functions.
!>
!> This is archived by providing small types, for the actual procedure reference
!> and an optional polymorphic data container, which are carried along in
!> NLopt as opaque data pointer and resolved again in the callback wrapper
!> before invoking the actual procedure.
!>
!> While the wrapper types reference the data, they do not own the data,
!> therefore the caller side is responsible to keep the object alive as long
!> it is referenced from the [[nlopt_opt]] instance. One option is to create
!> the wrappers by association:
!>
!>```fortran
!>associate(f => nlopt_func(objective_function, data_container)
!>  call create(opt, algoritm, n)
!>  ! ...
!>  call opt%set_min_objective(f)
!>  ! ...
!>  call opt%optimize(stat)
!>  ! ...
!>  call destroy(opt)
!>end associate
!>```
module nlopt_callback
  use, intrinsic :: iso_c_binding
  implicit none

  abstract interface
    function nlopt_func_interface(x, gradient, func_data) result(f)
      import :: c_int, c_double, c_ptr
      implicit none
      real(c_double), intent(in) :: x(:)
      real(c_double), intent(inout), optional :: gradient(:)
      class(*), intent(in), optional :: func_data
      real(c_double) :: f
    end function

    subroutine nlopt_mfunc_interface(result, x, gradient, func_data)
      import :: c_int, c_double, c_ptr
      implicit none
      real(c_double), intent(inout) :: result(:)
      real(c_double), intent(in) :: x(:)
      real(c_double), intent(inout), optional :: gradient(:)
      class(*), intent(in), optional :: func_data
    end subroutine

    subroutine nlopt_precond_interface(x, v, vpre, func_data)
      import :: c_int, c_double
      implicit none
      real(c_double), intent(in) :: x(:)
      real(c_double), intent(in) :: v(:)
      real(c_double), intent(inout) :: vpre(:)
      class(*), intent(in), optional :: func_data
    end subroutine
  end interface

  type :: nlopt_func
    private
    procedure(nlopt_func_interface), nopass, pointer :: f => null()
    class(*), pointer :: f_data => null()
  end type nlopt_func

  interface nlopt_func
    module procedure :: create_nlopt_func
  end interface nlopt_func

  type :: nlopt_mfunc
    private
    procedure(nlopt_mfunc_interface), nopass, pointer :: f => null()
    class(*), pointer :: f_data => null()
  end type nlopt_mfunc

  interface nlopt_mfunc
    module procedure :: create_nlopt_mfunc
  end interface nlopt_mfunc

  type :: nlopt_precond
    private
    procedure(nlopt_func_interface), nopass, pointer :: f => null()
    procedure(nlopt_precond_interface), nopass, pointer :: pre => null()
    class(*), pointer :: f_data => null()
  end type nlopt_precond

  interface nlopt_precond
    module procedure :: create_nlopt_precond
  end interface nlopt_precond

contains

  function create_nlopt_func(f, f_data) result(self)
    procedure(nlopt_func_interface) :: f
    class(*), intent(in), target, optional :: f_data
    type(nlopt_func) :: self

    self%f => f
    nullify(self%f_data)
    if (present(f_data)) self%f_data => f_data
  end function create_nlopt_func

  function wrap_func(n, x, gradient, func_data) result(f) &
      bind(c, name="nlopt_wrap_func")
    integer(c_int), value :: n
    real(c_double), intent(in) :: x(n)
    real(c_double), intent(inout), optional :: gradient(n)
    type(c_ptr), value :: func_data
    real(c_double) :: f

    type(nlopt_func), pointer :: func
    real(c_double), allocatable :: stub(:)

    call c_f_pointer(func_data, func)
    f = func%f(x, gradient, func%f_data)
  end function

  function create_nlopt_mfunc(f, f_data) result(self)
    procedure(nlopt_mfunc_interface) :: f
    class(*), intent(in), target, optional :: f_data
    type(nlopt_mfunc) :: self

    self%f => f
    nullify(self%f_data)
    if (present(f_data)) self%f_data => f_data
  end function create_nlopt_mfunc

  subroutine wrap_mfunc(m, result, n, x, gradient, func_data) &
      bind(c, name="nlopt_wrap_mfunc")
    integer(c_int), value :: m
    real(c_double), intent(inout) :: result(m)
    integer(c_int), value :: n
    real(c_double), intent(in) :: x(n)
    real(c_double), intent(inout), optional :: gradient(n)
    type(c_ptr), value :: func_data

    type(nlopt_mfunc), pointer :: mfunc

    call c_f_pointer(func_data, mfunc)
    call mfunc%f(result, x, gradient, mfunc%f_data)
  end subroutine

  function create_nlopt_precond(f, pre, f_data) result(self)
     procedure(nlopt_func_interface) :: f
     procedure(nlopt_precond_interface) :: pre
     class(*), intent(in), target, optional :: f_data
     type(nlopt_precond) :: self

     self%f => f
     self%pre => pre
    nullify(self%f_data)
     if (present(f_data)) self%f_data => f_data
  end function create_nlopt_precond

  subroutine wrap_precond(n, x, v, vpre, func_data) &
      bind(c, name="nlopt_wrap_precond")
    integer(c_int), value :: n
    real(c_double), intent(in) :: x(n)
    real(c_double), intent(in) :: v(n)
    real(c_double), intent(inout) :: vpre(n)
    type(c_ptr), value :: func_data

    type(nlopt_precond), pointer :: precond

    call c_f_pointer(func_data, precond)
    call precond%pre(x, v, vpre, precond%f_data)
  end subroutine wrap_precond

end module nlopt_callback

max_stars_repo_name: Nicholaswogan/nlopt-f

max_issues_repo_name: Nicholaswogan/nlopt-f

max_forks_repo_name: Nicholaswogan/nlopt-f


text:
!> @file general_constant.f90
!! @brief Global Constant Definition
!! @author Innocent Souopgui
!!
!<

!>
!!
!<
module general_constant
    use compiler_dependant
implicit none
    !> @brief Unknown integer value
    !<
    integer, parameter :: ip_unknown = -999!!unknown value
    !> @brief size of simple precision real variables
    !<
    integer, parameter :: sp = kind(1.0)
    !> @brief size of double precision real variables
    !<
    integer, parameter :: dp = kind(1d0)
    !> @brief size of the default precision for real variables, stand from computation precision
    !<
    integer, parameter :: cp = dp
    !   !> @brief size of the default precision for graphic, (dislin), stand from graphic precision
    !   !<
    !   integer, parameter :: gp = sp

    !>
    !! \enum noname
    !! @brief String Comparison Member Function status.
    !!
    !! \var ASC explanation 1
    !! \var DESC explanation 2
    !<
    enum, bind(c)
        enumerator :: ASC , DESC !< descending order
        enumerator :: ASCENDING = ASC, DESCENDING = DESC
    end enum

    enum, bind(c)
        enumerator :: MASKED = 1 !< masking mark
        enumerator :: UNMASKED = 0 !<
    end enum

    !> @brief pi constant
    real(cp), parameter :: rp_pi = 3.14159265358979323846_cp
    !> @brief small number for real values comparison, used as zero
    real(cp), parameter :: rp_epsilon = REAL(1d-15, cp)
    !> @brief Strings length, used for short file name or other short string variables
    !<
    integer, parameter :: ip_snl = 80
    !> @brief maximum length of the string representation of common constants
    !<
    integer, parameter :: ip_ccl = 20
    !> @brief Strings length, used for full file name or other short string variables
    !<
    integer, parameter :: ip_fnl = 255
    !> @brief format to print integer with comment
    !<
    character(len=*), parameter :: PRINT_IFORMAT = "(A,I5)"
    !> @brief format to print real with comment
    !<
    character(len=*), parameter :: PRINT_RFORMAT = "(A,E10.2)"
    !> @brief old status for file input/output
    !<
    character(len=*), parameter :: FILE_OLD = 'OLD'
    !> @brief replaced status for file input/output
    !<
    character(len=*), parameter :: FILE_REPLACE = 'REPLACE'
    !> @brief format for file input/output
    !<
    character(len=*), parameter :: FILE_FORMATTED = 'FORMATTED'
    !> @brief string delimiter
    !! @details This delimiter is used when saving string or characters data into namelist
    !! Usefull for special characters in namelist, for example \ and &
    !<
    character(len=*), parameter :: SD= "'"
    !> @brief constant for array vectorization, concatenation of rows(row by row)
    integer, parameter :: ROW_WISE =1341
    !> @brief constant for array vectorization, concatenation of columns(column by column)
    integer, parameter :: COL_WISE =1342

    !> @brief standard output file id
    !<
    integer, parameter :: STDOUT = 6

end module general_constant

max_stars_repo_name: souopgui/da-toolbox

max_issues_repo_name: souopgui/da-toolbox

max_forks_repo_name: souopgui/da-toolbox


text:
SUBROUTINE rad_full()

!++++++++++++++TAK: PBL interface for RRTM based on P. Blossey's GCSS-CFMIP interface+++++++++++++
! (1) The USEd parameters were organized for SAM6.10.4.
! (2) Adjustment for SAM6.10.4: t00(=300K) is added to sstxy (SST perturbation) for tg_slice.
! (3) Following RAD_RRTM/rad.f90, an option to pass the effective radius computed with Get_reffc
!     and Get_reffi defined in microphysics.f90 was added.
! (4) The CALL statment for rad_driver_rrtm is modified to pass computed effective radius.
! (5) Modifications were added to the subroutine, read_patch_background_sounding, so that it
!     patches a default sounding ported from WRF3.3.1 (module_ra_rrtmg_lw.f90) when no input file
!     is specified.
! (6) read_rad & write_rad were created for restart runs.
! (7) Several non-standard options were disabled: use_SCAM_ozon, do_scale_co2, co2factor
!     When ozon exist in SCAM IOP file, it will be used,
!     i.e., use_SCAM_ozon = .TRUE. when ozon is in IOP file, vice versa.
! (8) includePrecip: if true, qpl & qpi are included in radiation calculation
! (9) Further modifications are necessary when trace-gases and above-domain profiles are not
!     constant in time.
!
! 2017/01
! - Deleted the disabled non-standard options for cleaner code
! - Code modifications for time-varying patched soundings
! - THIS IS IMPORTANT: npatch_start, npatch_end, nzpatch are assumed to be constant throughout
!   simulation. This means that the domain top pressure always has to be between
!   psnd(npatch_start-1) and psnd(npatch_start). Initially presi(nz)-10 >= psnd(npatch_start), so
!   that resolution of pressure level for psnd should be larger than 20 mb or so if anticipated
!   presi(nz) change is less than 20 mb, which likely guarantees psnd(npatch_start-1) < presi(nz)
!   < psnd(npatch_start). Ways to mitigate this requirement are
!   - Have deep enough domain so that the domain top pressure does not change much,
!   - Make resolution of the background sounding (IOP file) coarse enough.
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

	USE rad ! nradsteps, qrad, pres_input, tabs_slice, insolation_TOA, lwUp, etc from this module
	USE rad_driver, ONLY: rad_driver_rrtm, initialize_radiation, isInitialized_RadDriver, &
		tracesini, tracesupdate, read_rad, write_rad, land_frac
	USE parkind, ONLY: kind_rb ! RRTM expects reals with this kind parameter (8 byte reals)

	!==========================================================================================
	!============================== BEGIN CHANGES REQUIRED HERE ===============================
	!==========================================================================================
	! PORTABILITY NOTE:  Here, all of the stuff needed to call the radiation is drawn from
	!   various modules in SAM as an example.  You will need to bring all of the associated
	!   variables into this routine for your model (either by using the appropriate modules,
	!   passing them in as arguments or defining them here).
	!==========================================================================================

	! Logicals
	USE grid, ONLY: &
		masterproc, &  ! true if MPI rank==0.
		dostatisrad, & ! accumulate radiation statistics at this time step
		compute_reffc, compute_reffi, & ! get effective radius with Get_reffc, Get_reffi in microphysics
                compute_reffl ! option for radiatively-active drizzle

	USE params, ONLY: &
		dolongwave, &  ! do longwave radiation
		doshortwave, & ! do shortwave radiation
		doperpetual, & ! use diurnally-averaged insolation
		dosolarconstant, & ! specify mean insolation, zenith angle for perpetual insolatoin
		doseasons, &   ! allow diurnally-varying insolation to vary with time of year
		ocean ! if true, run is over ocean

        ! Tak Yamaguchi 2015/10: RADAEROS, RESLAST
	USE params, ONLY: &
		doradaerosimple, &
		dorestart_last_only

	! Characters
	USE grid, ONLY: &
		iopfile, & ! name of SCAM IOP forcing file, e.g. 'ctl_s11.nc'
		case ! used to construct path of SCAM IOP file in SAM

	! Integers
	USE grid, ONLY: &
		nx, ny, nzm, & ! number of grid points (x,y,z)
		nstep, icycle, & ! model step and substep number
		nrad, & ! call radiation every nrad timesteps
		nrestart, & ! switch to control starting/restarting of the model
		nrestart_skip, & ! number of skips of writing restart (default 0)
		nstat, & ! the interval in time steps to compute statistics
		nstop, & ! time step number to stop the integration
		nelapse ! time step number to elapse before stoping

	! Reals or real arrays
	USE params, ONLY: &
		cp, & ! specific heat of dry air at constant pressure, J/kg/K
		ggr, & ! gravitational acceleration, m/s2
		coszrs, & ! cosine of solar zenith angle
		latitude0, longitude0, & ! latitude, longitude in degrees
		solar_constant, & ! modified solar constant doperpetual==dosolarconstant==.true.
		zenith_angle ! zenith angle if doperpetual==dosolarconstant==.true.

	USE grid, ONLY: &
		dtn, & ! time step in seconds
		day, day0, & ! model day (current and at start of run) day=0. for 00Z, Jan 1.
		dz, adz ! vertical grid spacing is dz*adz(k) in this model

	! NOTE: when dosolarconstat==.true, insolation = solar_constant*cos(zenith_angle)

	USE vars, ONLY: &
		daysnd, & ! day in background sounding for SCAM IOP file
		t, tabs, qv, qcl, qci, qpl, qpi, & ! model fields
		sstxy, t00, & ! SST perturbation and 300 K
		pres, presi, & ! pressure at model levels/interfaces
		rho, & ! density profile.  In this anelastic model, rho=rho(z).
		radswup, radswdn, radqrsw, & ! radiation statistics, summed in x and y.
		radlwup, radlwdn, radqrlw, & ! radiation statistics, summed in x and y.
		s_flntoa, s_fsntoa, s_flntoac, s_fsntoac, s_solin, & ! TOA radiative flux statistics
                s_flnt, s_fsnt, & ! top-of-model fluxes
		s_flns, s_fsns, s_flnsc, s_fsnsc, s_fsds, s_flds, & ! surface radiative flux statistics
		lwnt_xy, swnt_xy, lwntc_xy, swntc_xy, solin_xy, & ! time-accumulated x-y fields of TOA radiative fluxes
		lwns_xy, swns_xy, lwnsc_xy, swnsc_xy ! time-accumulated x-y fields of surface radiative fluxes

	USE microphysics, ONLY: reffc, reffi, reffr, & ! These arrays replace get_reffc, get_reffi
                 dorrtm_cloud_optics_from_effrad_LegacyOption, &
                 Get_nca ! aerosol number concentration function (Tak Yamaguchi 2015/10: RADAEROS)

        !bloss(2017/07): compute effective radii in this routine if not provided by microphysics.
        use cam_rad_parameterizations, only : computeRe_Liquid, computeRe_Ice

	!==========================================================================================
	!=============================== END CHANGES REQUIRED HERE ================================
	!==========================================================================================

	IMPLICIT NONE

	! Parameters
	INTEGER, PARAMETER :: iyear = 1999 ! (TAK) RAD_RRTM/rad.f90 uses iyear=1999

	! Local variables
	LOGICAL, SAVE :: isInitialized_BackgroundSounding = .FALSE.
	LOGICAL, SAVE :: isRestart = .TRUE. ! true only at the first step of new and restart run
	CHARACTER(LEN=250) :: SoundingFileName ! input file for background soundings
	REAL(KIND=kind_rb), DIMENSION(nx,nzm) :: &
		swHeatingRate, & ! units: K/s
		lwHeatingRate, & !  units: K/s
		swHeatingRateClearSky, & !  units: K/s
		lwHeatingRateClearSky !  units: K/s
	REAL :: coef
	INTEGER :: ierr, i, j, k, i1, i2
	INTEGER, SAVE :: index_daysnd = 1 ! index for daysnd used for the last update of background snds

	!++++++++++++TAK: Including precipitating condensate in radiation calculation?+++++++++++++
	! Set TRUE if precipitation condensate (qpl, qpi) is included in radiation calculation,
	! e.g., effective radius from microphysics.
        !bloss(2018-02): Change default to false since effective radius may not be consisitent
        !  with including rain in the liquid water content.  Similar for snow/graupel.
	LOGICAL :: includePrecip = .false.
	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

	!==========================================================================================
	! PORTABILITY NOTE: all processors should have same pressure soundings (pres and presi) in
	!   hPa. This is important for the consistency of trace soundings across the different
	!   processors.
	!==========================================================================================

	IF ( icycle == 1 ) THEN ! Skip subcycles (i.e., when icycle /= 1)

		IF ( .NOT.isInitialized_RadDriver ) THEN
			! This IF block is TRUE only once.

			! Call initialization routine for radiation
			CALL initialize_radiation( cp, iyear, day0, latitude0, longitude0, doperpetual, ocean )

			! Allocate arrays that are fixed in size for the whole simulation.
			ALLOCATE( &
				NetlwUpSurface(nx,ny), NetlwUpSurfaceClearSky(nx,ny), &
				NetlwUpToa(nx,ny),     NetlwUpToaClearSky(nx,ny),     &
				NetlwUpToM(nx,ny),     NetswDownToM(nx,ny),     &
				NetswDownSurface(nx,ny), NetswDownSurfaceClearSky(nx,ny), &
				NetswDownToa(nx,ny),     NetswDownToaClearSky(nx,ny),     &
				insolation_TOA(nx,ny), swDownSurface(nx,ny), lwDownSurface(nx,ny), &
				swnsxy(nx,ny), lwnsxy(nx,ny), &
				STAT=ierr )
			IF ( ierr /= 0 ) THEN
				WRITE(*,*) 'Could not allocate fixed size arrays in rad_full'
				CALL rad_error()
			ENDIF

			isInitialized_RadDriver = .TRUE.
		ENDIF ! IF (.NOT.isInitialized_RadDriver)

		IF ( .NOT.isInitialized_BackgroundSounding ) THEN
			! This IF block is TRUE only once.

			! Read in sounding and define set up patching.
			! This will be stored in radiation restart file for restarts.
			SoundingFileName = './'//TRIM(case)//'/'//TRIM(iopfile)

			! Specify SoundingFileName as NoInputFile if there is no input file (TAK)
			IF( LEN('./'//TRIM(case)//'/') == LEN(TRIM(SoundingFileName)) ) &
				SoundingFileName='NoInputFile'

			! Setup background soundings
			! - Allocate psnd, tsnd, qsnd, o3snd, tsndng, qsndng, o3sndng
			! - Set have_o3mmr
			! - Read either default sounding or SCAM IOP file
			CALL read_patch_background_sounding( SoundingFileName, presi(nzm+1), masterproc )

			! nzpatch = number of patched levels (above domain top)
			! nzrad = total number of level for radiation
			IF ( npatch_end /= npatch_start ) THEN
				nzpatch = npatch_end - npatch_start + 1
			ELSE
				nzpatch = 0
			ENDIF
			nzrad = nzm + nzpatch

			IF ( masterproc ) WRITE(*,*) 'In RRTM radiation interface, using ozone profile'
			IF ( have_o3mmr .AND. masterproc ) THEN
				WRITE(*,*) '  from SCAM IOP netcdf file', SoundingFileName
			ELSE IF ( masterproc ) THEN
				WRITE(*,*) '  from default absorber profiles in rrtmg_lw.nc'
			ENDIF

			isInitialized_BackgroundSounding = .TRUE.
		ENDIF ! IF ( .NOT.isInitialized_BackgroundSounding )

		IF ( isAllocated_RadInputsOutputs .AND. nzrad /= nzrad_old ) THEN
			! This IF block is TRUE when nzpatch and nzrad were changed because the domain top
			! pressure becomes either smaller than psnd at previous npatch_start or larger than
			! psnd at previous npatch_start-1.
			!
			! This IF block will be FALSE since nzrad will not be updated. This means that nzrad
			! is assumed to be constant. Have deep enough domain depth.

			! Deallocate old arrays
			! These variables have not been allocated yet at initialization/restart, thus this IF
			! block will be skipped.
			DEALLOCATE( &
				tabs_slice, qv_slice, qcl_slice, qci_slice, rel_slice, rei_slice, tg_slice, &
				nca_slice, &! Tak Yamaguchi 2015/10: RADAEROS
				pres_input, presi_input, lwUp, lwDown, lwUpClearSky, lwDownClearSky, &
				swUp, swDown, swUpClearSky, swDownClearSky, STAT=ierr )
			IF ( ierr /= 0 ) THEN
				WRITE(*,*) 'Could not deallocate input/output arrays in rad_full'
				CALL rad_error()
			ENDIF
			IF ( compute_reffc ) THEN
!bloss				DEALLOCATE( rel_rad, STAT=ierr )
!bloss				IF ( ierr /= 0 ) THEN
!bloss					WRITE(*,*) 'Could not deallocate rel_rad(nx,ny,nzm) in rad_full'
!bloss					CALL rad_error()
!bloss				ENDIF
			ENDIF
			IF ( compute_reffi ) THEN
!bloss				DEALLOCATE( rei_rad, STAT=ierr )
!bloss				IF ( ierr /= 0 ) THEN
!bloss					WRITE(*,*) 'Could not deallocate rei_rad(nx,ny,nzm) in rad_full'
!bloss					CALL rad_error()
!bloss				ENDIF
			ENDIF

                       ! Tak Yamaguchi 2015/10: RADAEROS
			IF ( doradaerosimple ) THEN
				DEALLOCATE( nca_rad, STAT=ierr )
				IF ( ierr /= 0 ) THEN
					WRITE(*,*) 'Could not deallocate nca_rad(nx,ny,nzm) in rad_full'
					CALL rad_error()
				ENDIF
			ENDIF

			isAllocated_RadInputsOutputs = .FALSE.
		ENDIF ! IF ( isAllocated_RadInputsOutputs .AND. nzrad /= nzrad_old )

		IF ( .NOT.isAllocated_RadInputsOutputs ) THEN
			! This IF block is TRUE only once for this version.
			! When nzrad is allowed to be changed, then this block should be TRUE whenever the above
			! IF block is TRUE.

			! Allocate arrays
			ALLOCATE( &
				tabs_slice(nx,nzrad), qv_slice(nx,nzrad), qcl_slice(nx,nzrad), &
				qci_slice(nx,nzrad), rel_slice(nx,nzrad), rei_slice(nx,nzrad), tg_slice(nx), &
				nca_slice(nx,nzrad), & ! Tak Yamaguchi 2015/10: RADAEROS
				pres_input(nzrad), presi_input(nzrad+1), &
				lwUp(nx,nzrad+2), lwDown(nx,nzrad+2), &
				lwUpClearSky(nx,nzrad+2), lwDownClearSky(nx,nzrad+2), &
				swUp(nx,nzrad+2), swDown(nx,nzrad+2), &
				swUpClearSky(nx,nzrad+2), swDownClearSky(nx,nzrad+2), &
				STAT=ierr )
			IF ( ierr /= 0 ) THEN
				WRITE(*,*) 'Could not allocate input/output arrays in rad_full'
				CALL rad_error()
			ENDIF
			IF ( compute_reffc ) THEN
!bloss				ALLOCATE( rel_rad(nx,ny,nzm), STAT=ierr )
!bloss				IF ( ierr /= 0 ) THEN
!bloss					WRITE(*,*) 'Could not allocate rel_rad(nx,ny,nzm) in rad_full'
!bloss					CALL rad_error()
!bloss				ENDIF
				rad_reffc(:,:,:) = 0.0
			ENDIF
			IF ( compute_reffi ) THEN
!bloss				ALLOCATE( rei_rad(nx,ny,nzm), STAT=ierr )
!bloss				IF ( ierr /= 0 ) THEN
!bloss					WRITE(*,*) 'Could not allocate rei_rad(nx,ny,nzm) in rad_full'
!bloss					CALL rad_error()
!bloss				ENDIF
				rad_reffi(:,:,:) = 0.0
			ENDIF

                        ! Tak Yamaguchi 2015/10: RADAEROS
			IF ( doradaerosimple ) THEN
				ALLOCATE( nca_rad(nx,ny,nzm), STAT=ierr )
				IF ( ierr /= 0 ) THEN
					WRITE(*,*) 'Could not allocate nca_rad(nx,ny,nzm) in rad_full'
					CALL rad_error()
				ENDIF
				nca_rad(:,:,:) = 0.0
			ENDIF

			isAllocated_RadInputsOutputs = .TRUE.
			nzrad_old = nzrad
		ENDIF ! IF ( .NOT.isAllocated_RadInputsOutputs )

		IF ( nstep == 1 .OR. isRestart ) THEN
			! Set up pressure inputs to radiation -- needed for initialize_radiation
			pres_input(1:nzm) = pres(1:nzm)
			presi_input(1:nzm+1) = presi(1:nzm+1)
			IF ( nzpatch > 0 ) THEN
				pres_input(nzm+1:nzrad) = psnd(npatch_start:npatch_end) ! layer pressures
				presi_input(nzm+2:nzrad) & ! interface pressures.
					= 0.5*(psnd(npatch_start:npatch_end-1) + psnd(npatch_start+1:npatch_end))
				presi_input(nzrad+1) &
					= MAX(0.5*psnd(npatch_end), 1.5*psnd(npatch_end) - 0.5*psnd(npatch_end-1))
			ENDIF

			! Allocate trace gases arrays used in RRTM.
			! Interpolates standard sounding of trace gas concentrations to grid for radiation.
			! For restart runs, trace gas profiles will be overwritten when they are read from the
			! restart data. Thus no time interpolation for o3snd is necessary.
			CALL tracesini( nzrad, pres_input, presi_input, ggr, nzsnd, psnd, o3snd, have_o3mmr, &
				1.0, masterproc )
		ENDIF ! IF ( nstep == 1 .OR. isRestart )

		! Read restart data at restart
		! This will read nradsteps, npatch_start, npatch_end, nzpatch, psnd, tsnd, qsnd, other
		! allocatables (qrad, etc. declared in vars.d90, rad.f90), and trace gases (o3, co2, etc.
		! declared in rad_driver.f90).
		IF ( isRestart ) THEN
			IF ( nrestart /= 0 ) CALL read_rad() ! Read radiation restart data
			isRestart = .FALSE.
		ENDIF

		! Update radiation
		nradsteps = nradsteps + 1
		IF ( nstep == 1 .OR. nradsteps >= nrad ) THEN
			! Zero out radiation statistics that are summed in x- and y-directions.
			radlwup(:) = 0.0
			radlwdn(:) = 0.0
			radqrlw(:) = 0.0
			radswup(:) = 0.0
			radswdn(:) = 0.0
			radqrsw(:) = 0.0
			radqrclw(:) = 0.0
			radqrcsw(:) = 0.0

			! Check if presi(nz) is within psnd(npatch_start-1) and psnd(npatch_start).
			! This is necessary since npatch_start (so as nzpatch, nzrad) is constant with time.
			! psnd(npatch_start-1) < presi(nz) is perhaps OK as long as psnd(npatch_start-1) ~
			! presi(nz), but I include this is also an error.
                        !bloss(2017-07): Make this only a warning/notice if presi(nz) < psnd(npatch_start-2)
			IF ( nzpatch > 0 ) THEN
				IF ( psnd(npatch_start-1) < presi(nzm+1) .OR. presi(nzm+1) < psnd(npatch_start) ) THEN
					IF ( masterproc ) THEN
                                           PRINT*,'Patching background sounding on top of dynamical model domain'
                                           PRINT*,'  for radiation computation.  '
                                           write(*,992) presi(1), presi(nzm+1), psnd(npatch_start), psnd(npatch_end)
                                           992 format('Dynamical model domain: ',F6.1,' < p (hPa) < ',F6.1, &
                                                    ', Patched sounding above: ',F6.1,' < p (hPa) < ',F6.1)
                                         END IF
                                         IF(abs(presi(nzm+1)-psnd(npatch_start)).gt.100.) THEN
                                           if(masterproc) PRINT*,'Gap seems too big, so stopping...'
                                           CALL task_abort()
                                         END IF
                                       ENDIF
			ENDIF

			! Set up pressure inputs to radiation
			! THIS MAY NOT BE UNNECESSARY WHEN PRES AND PRESI DO NOT CHANGE MUCH.
			pres_input(1:nzm) = pres(1:nzm)
			presi_input(1:nzm+1) = presi(1:nzm+1)
			IF ( nzpatch > 0 ) THEN
				pres_input(nzm+1:nzrad) = psnd(npatch_start:npatch_end) ! layer pressures
				presi_input(nzm+2:nzrad) & ! interface pressures
					= 0.5 * ( psnd(npatch_start:npatch_end-1) + psnd(npatch_start+1:npatch_end) )
				presi_input(nzrad+1) &
					= MAX( 0.5*psnd(npatch_end), 1.5*psnd(npatch_end) - 0.5*psnd(npatch_end-1) )
			ENDIF

			! Time interpolate background soundings
			IF ( ntsnd /= 1 ) THEN
				! Background profiles are from a SCAM IOP file
				! Use daysnd to identify the two closest times to the current time
				DO i = index_daysnd, ntsnd
					IF ( day < daysnd(i) ) THEN
						index_daysnd = i - 1
						EXIT
					ENDIF
				ENDDO
				i1 = index_daysnd
				i2 = i1 + 1
				coef = ( day - daysnd(i1) ) / ( daysnd(i2) - daysnd(i1) )
				tsnd(:) = tsndng(:,i1) + ( tsndng(:,i2) - tsndng(:,i1) ) * coef
				qsnd(:) = qsndng(:,i1) + ( qsndng(:,i2) - qsndng(:,i1) ) * coef

				! Update o3snd, then o3 following rad_driver.f90/tracesini
				IF ( have_o3mmr ) THEN
					o3snd(:) = o3sndng(:,i1) + ( o3sndng(:,i2) - o3sndng(:,i1) ) * coef
					CALL tracesupdate( nzrad, pres_input, presi_input, ggr, nzsnd, psnd, o3snd, &
						have_o3mmr, masterproc )
				ENDIF
			ENDIF

			! Obtain effective radius from microphysics.f90 (TAK)
                        !   or compute it using the CAM parameterizations.
			IF ( compute_reffc ) then
                           rad_reffc(1:nx,1:ny,1:nzm) = reffc(:,:,:)
                        ELSE
                           do k = 1,nzm
                             do j = 1,ny
                               do i = 1,nx
                                 rad_reffc(i,j,k) = computeRe_Liquid( tabs(i,j,k), land_frac, 0.0_kind_rb, 0.0_kind_rb )
                               end do
                             end do
                           end do
                         END IF

			IF ( compute_reffi ) THEN
                           rad_reffi(1:nx,1:ny,1:nzm) = reffi(:,:,:)
                        ELSE
                          rad_reffi(1:nx,1:ny,1:nzm) = computeRe_Ice( tabs(1:nx,1:ny,1:nzm) )
                        END IF

                        !bloss(2018-02): option for drizzle/rain being radiatively active
                        if(compute_reffl.AND.dorrtm_cloud_optics_from_effrad_LegacyOption) then
                          !bloss(2018-02): Combine cloud and rain/drizzle to get a composite species
                          !   with an average optical depth
                          do k = 1,nzm
                            do j = 1,ny
                              do i = 1,nx
                                if(qpl(i,j,k)+qcl(i,j,k).gt.0._kind_rb) then
                                  rad_reffc(i,j,k) = ( qcl(i,j,k) + qpl(i,j,k) ) &
                                       / (qcl(i,j,k)/reffc(i,j,k) + qpl(i,j,k)/reffr(i,j,k) )
                                end if
                              end do
                            end do
                          end do
!bloss(2018-02)                          rad_reffc(1:nx,1:ny,1:nzm) = reffl(1:nx,1:ny,1:nzm)
                          includePrecip = .true. ! include rain as radiatively active
                        end if

                        ! Tak Yamaguchi 2015/10: RADAEROS
			! Obtain aerosol number concentration from microphysics.f90
			IF ( doradaerosimple ) nca_rad(1:nx,1:ny,1:nzm) = Get_nca() ! #/cm3

			!
			DO j = 1, ny
				! Extract a slice from the three-dimensional domain on this processor.
				!   We need absolute temperature (K), mass mixing ratios (kg/kg) of
				!   water vapor, cloud liquid water and cloud ice, along with SST (K).
				tabs_slice(1:nx,1:nzm) = tabs(1:nx,j,1:nzm)
				qv_slice(1:nx,1:nzm) = qv(1:nx,j,1:nzm)
				IF ( .NOT.includePrecip ) THEN
					qcl_slice(1:nx,1:nzm) = qcl(1:nx,j,1:nzm)
					qci_slice(1:nx,1:nzm) = qci(1:nx,j,1:nzm)
				ELSE
					qcl_slice(1:nx,1:nzm) = qcl(1:nx,j,1:nzm) + qpl(1:nx,j,1:nzm)
					qci_slice(1:nx,1:nzm) = qci(1:nx,j,1:nzm) + qpi(1:nx,j,1:nzm)
				ENDIF
				tg_slice(1:nx) = sstxy(1:nx,j) + t00

                                !bloss(2017/07):
                                rel_slice(1:nx,1:nzm) = rad_reffc(1:nx,j,1:nzm) ! (TAK)
				rei_slice(1:nx,1:nzm) = rad_reffi(1:nx,j,1:nzm) ! (TAK)

                                ! Tak Yamaguchi 2015/10: RADAEROS
				! Extract aerosol number concentration
				IF ( doradaerosimple ) nca_slice(1:nx,1:nzm) = nca_rad(1:nx,j,1:nzm)

				IF ( nzpatch > 0)  THEN
					! Patch sounding on top of model sounding for more complete radiation calculation.
					tabs_slice(1:nx,nzm+1:nzrad) &
						= SPREAD( tsnd(npatch_start:npatch_end), DIM=1, NCOPIES=nx )
					qv_slice(1:nx,nzm+1:nzrad) &
						= SPREAD( qsnd(npatch_start:npatch_end), DIM=1, NCOPIES=nx )
					qcl_slice(1:nx,nzm+1:nzrad) = 0.0
					qci_slice(1:nx,nzm+1:nzrad) = 0.0
					IF ( compute_reffc ) rel_slice(1:nx,nzm+1:nzrad) = 0.0 ! (TAK)
					IF ( compute_reffi ) rei_slice(1:nx,nzm+1:nzrad) = 0.0 ! (TAK)

                                        ! Tak Yamaguchi 2015/10: RADAEROS
					! Patch aerosol number concentration
                                        ! FOR THE CURRENT VERSION THIS IS MEANINGLESS SINCE AEROSOL OPTICAL DEPTH IS COMPUTED BELOW NZM
                                        ! SINCE THERE IS NO INFORMATION ABOUT DZ IN THE PATCHED SOUNDING. PATCHED SOUNDING IS BASED ON
                                        ! PRESSURE.
					IF ( doradaerosimple ) THEN
						DO k = nzm+1, nzrad
							nca_slice(1:nx,k) = nca_slice(1:nx,nzm)
						ENDDO
					ENDIF
				ENDIF

				! Make call to wrapper routine for RRTMG (v.4.8 for LW, v.3.8 for SW)
				CALL rad_driver_rrtm( nx, nzrad, nzm, j, pres_input, presi_input, &
					tabs_slice, qv_slice, qcl_slice, qci_slice, rel_slice, rei_slice, tg_slice, &
					nca_slice, & ! Tak Yamaguchi 2015/10: RADAEROS
					dolongwave, doshortwave, doperpetual, doseasons, &
					dosolarconstant, solar_constant, zenith_angle, &
					day, day0, latitude0, longitude0, &
					ocean, ggr, cp, masterproc, &
					lwUp, lwDown, lwUpClearSky, lwDownClearSky, &
					swUp, swDown, swUpClearSky, swDownClearSky, coszrs )

				! Compute heating rates from fluxes using local density in model.
				! Results in heating rates in K/s.
				! PORTABILIY NOTE: CHANGE THERMAL MASS TO THAT USED IN YOUR MODEL.
				!   Units below are cp*rho*deltaz ~ J/kg/K * kg/m3 * m
				!   where delta z = dz*adz(k) in SAM.
				DO k = 1, nzm ! loop over model levels
					swHeatingRate(1:nx,k) &
						= (swDown(:,k+1) - swDown(:,k) + swUp(:,k) - swUp(:,k+1)) / (cp*rho(k)*dz*adz(k))
					lwHeatingRate(1:nx,k) &
						= (lwDown(:,k+1) - lwDown(:,k) + lwUp(:,k) - lwUp(:,k+1)) / (cp*rho(k)*dz*adz(k))
					swHeatingRateClearSky(1:nx,k) &
						= ( swDownClearSky(:,k+1) - swDownClearSky(:,k) &
						+ swUpClearSky(:,k) - swUpClearSky(:,k+1) ) &
						/ (cp*rho(k)*dz*adz(k))
					lwHeatingRateClearSky(1:nx,k) &
						= ( lwDownClearSky(:,k+1) - lwDownClearSky(:,k) &
						+ lwUpClearSky(:,k) - lwUpClearSky(:,k+1) ) &
						/ (cp*rho(k)*dz*adz(k))
				ENDDO

				! Update total radiative heating rate of model
				qrad(1:nx,j,1:nzm) = swHeatingRate(1:nx,1:nzm) + lwHeatingRate(1:nx,1:nzm)

				! Accumulate heating rates and fluxes for horizontally-averaged statistics
				radlwup(:) = radlwup(:) + SUM(lwUp(:, 1:nzm+1),   DIM = 1)
				radlwdn(:) = radlwdn(:) + SUM(lwDown(:, 1:nzm+1), DIM = 1)
				radqrlw(1:nzm) = radqrlw(1:nzm) + SUM(lwHeatingRate(:, 1:nzm), DIM = 1)
				radqrclw(1:nzm) = radqrclw(1:nzm) + sum(lwHeatingRateClearSky(:, 1:nzm), dim = 1)
				radswup(:) = radswup(:) + SUM(swUp(:, 1:nzm+1),   DIM = 1)
				radswdn(:) = radswdn(:) + SUM(swDown(:, 1:nzm+1), DIM = 1)
				radqrsw(1:nzm) = radqrsw(1:nzm) + SUM(swHeatingRate(:, 1:nzm), DIM = 1)
				radqrcsw(:nzm) = radqrcsw(:nzm) + sum(swHeatingRateClearSky(:, 1:nzm), dim = 1)

				! Shortwave fluxes at top-of-atmosphere (TOA) and surface -- NOTE POSITIVE DOWNWARDS
				insolation_TOA(1:nx,j) = swDown(:,nzrad+2) ! shortwave down at TOA
				swDownSurface(1:nx,j) = swDown(:,1) ! shortwave down at surface

				NetswDownToM(1:nx,j) = swDown(:,nzm+1) - swUp(:,nzm+1) ! net shortwave down at Top of Model
				NetswDownToa(1:nx,j) = swDown(:,nzrad+2) - swUp(:,nzrad+2) ! net shortwave down at TOA
				NetswDownToaClearSky(1:nx,j) = swDownClearSky(:,nzrad+2) - swUpClearSky(:,nzrad+2) ! net clearsky shortwave down at TOA

				NetswDownSurface(1:nx,j) = swDown(:,1) - swUp(:,1) ! net shortwave down at surface
				NetswDownSurfaceClearSky(1:nx,j) = swDownClearSky(:,1) - swUpClearSky(:,1) ! net clearsky shortwave down at surface

				! longwave fluxes at top-of-atmosphere (TOA) and surface -- NOTE POSITIVE UPWARDS
				lwDownSurface(1:nx,j) = lwDown(:,1) ! longwave down at surface

				NetlwUpToM(1:nx,j) = lwUp(:,nzm+1) - lwDown(:,nzm+1) ! net longwave up at Top of Model
				NetlwUpToa(1:nx,j) = lwUp(:,nzrad+2) - lwDown(:,nzrad+2) ! net longwave up at TOA
				NetlwUpToaClearSky(1:nx,j) = lwUpClearSky(:,nzrad+2) - lwDownClearSky(:,nzrad+2) ! net clearsky longwave up at TOA

				NetlwUpSurface(1:nx,j) = lwUp(:,1) - lwDown(:,1) ! net longwave up at surface
				NetlwUpSurfaceClearSky(1:nx,j) = lwUpClearSky(:,1) - lwDownClearSky(:,1) ! net clearsky longwave up at surface

			ENDDO ! j = 1, ny
			! Re-initialize nradsteps for the next update
			nradsteps = 0
		ENDIF ! IF (updateRadiation)

		! Update 2d diagnostic fields
		! Net surface and toa fluxes
		! First two for ocean evolution
		lwnsxy(:, :) = NetlwUpSurface(:, :) ! instantaneous
		swnsxy(:, :) = NetswDownSurface(:, :)

		! net full sky radiative fluxes (varying in x and y, time-accumulated)
		lwns_xy(:, :) = lwns_xy(:, :) + NetlwUpSurface(:, :)
		swns_xy(:, :) = swns_xy(:, :) + NetswDownSurface(:, :)
		lwnt_xy(:, :) = lwnt_xy(:, :) + NetlwUpToa(:, :)
		swnt_xy(:, :) = swnt_xy(:, :) + NetswDownToa(:, :)

		! net clear sky radiative fluxes (varying in x and y, time-accumulated)
		lwnsc_xy(:, :) = lwnsc_xy(:, :) + NetlwUpSurfaceClearSky(:, :)
		swnsc_xy(:, :) = swnsc_xy(:, :) + NetswDownSurfaceClearSky(:, :)
		lwntc_xy(:, :) = lwntc_xy(:, :) + NetlwUpToaClearSky(:, :)
		swntc_xy(:, :) = swntc_xy(:, :) + NetswDownToaClearSky(:, :)

		! TOA Insolation
		solin_xy(:, :) = solin_xy(:, :) + insolation_Toa(:, :)

		! Update 1D diagnostics
		IF ( dostatisrad ) THEN
			s_flns = s_flns + sum(NetlwUpSurface(:, :))   ! lwnsxy
			s_fsns = s_fsns + sum(NetswDownSurface(:, :)) ! swnsxy
			s_flnt = s_flnt + sum(NetlwUpToM(:, :))       ! LW Net Up Top of Model
			s_fsnt = s_fsnt + sum(NetswDownToM(:, :))     ! SW Net Dn Top of Model
			s_flntoa = s_flntoa + sum(NetlwUpToa(:, :))       ! lwntxy TOA
			s_fsntoa = s_fsntoa + sum(NetswDownToa(:, :))     ! swntxy TOA
			s_flnsc = s_flnsc + sum(NetlwUpSurfaceClearSky(:, :))   ! lwnscxy
			s_fsnsc = s_fsnsc + sum(NetswDownSurfaceClearSky(:, :)) ! swnscxy
			s_flntoac = s_flntoac + sum(NetlwUpToaClearSky(:, :))       ! lwntcxy
			s_fsntoac = s_fsntoac + sum(NetswDownToaClearSky(:, :))     ! swntcxy
			s_solin = s_solin + sum(insolation_TOA(:, :))           ! solinxy
			!
			s_fsds = s_fsds + sum(swDownSurface(:, :)) ! downwelling SW at surface
			s_flds = s_flds + sum(lwDownSurface(:, :)) ! downwelling LW at surface
		ENDIF ! if(dostatis)

		! Write radiation restart file
		IF ( MOD(nstep,nstat*(1+nrestart_skip)) == 0.OR.nstep == nstop.OR.nelapse == 0 ) THEN

                        ! Tak Yamaguchi 2015/10 RESLAST
			IF ( .NOT.dorestart_last_only &
			     .OR. ( dorestart_last_only .AND. (nstep == nstop .OR. nelapse == 0 ) ) ) &
			CALL write_rad()
		ENDIF

	ENDIF ! IF (icycle == 1)

	! Add radiative heating to model thermodynamic variable.
	! (here t is liquid-ice static energy divided by Cp.)
	! PORTABILITY NOTE: don't forget exner function if your model uses theta.
	t(1:nx, 1:ny, 1:nzm) = t(1:nx, 1:ny, 1:nzm) + qrad(:, :, :) * dtn

END SUBROUTINE rad_full

max_stars_repo_name: kbren/uwnet

max_issues_repo_name: kbren/uwnet

max_forks_repo_name: kbren/uwnet


text:
      SUBROUTINE PW92C( NSPIN, DENS, EC, VC )

C ********************************************************************
C Implements the Perdew-Wang'92 local correlation (beyond RPA).
C Ref: J.P.Perdew & Y.Wang, PRB, 45, 13244 (1992)
C Written by L.C.Balbas and J.M.Soler. Dec'96.  Version 0.5.
C ********* INPUT ****************************************************
C INTEGER NSPIN       : Number of spin polarizations (1 or 2)
C REAL*8  DENS(NSPIN) : Local (spin) density
C ********* OUTPUT ***************************************************
C REAL*8  EC        : Correlation energy density
C REAL*8  VC(NSPIN) : Correlation (spin) potential
C ********* UNITS ****************************************************
C Densities in electrons per Bohr**3
C Energies in Hartrees
C ********* ROUTINES CALLED ******************************************
C None
C ********************************************************************

C Next line is nonstandard but may be supressed
      IMPLICIT          NONE

C Argument types and dimensions
      INTEGER           NSPIN
      DOUBLE PRECISION  DENS(NSPIN), EC, VC(NSPIN)

C Internal variable declarations
      INTEGER           IG
      DOUBLE PRECISION  A(0:2), ALPHA1(0:2), B, BETA(0:2,4), C,
     .                  DBDRS, DECDD(2), DECDRS, DECDZ, DENMIN, DFDZ,
     .                  DGDRS(0:2), DCDRS, DRSDD, DTOT, DZDD(2),
     .                  F, FPP0, FOUTHD, G(0:2), HALF,
     .                  P(0:2), PI, RS, THD, THRHLF, ZETA

C Fix lower bound of density to avoid division by zero
      PARAMETER ( DENMIN = 1.D-12 )

C Fix some numerical constants
      PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
     .            THD=1.D0/3.D0, THRHLF=1.5D0 )

C Parameters from Table I of Perdew & Wang, PRB, 45, 13244 (92)
      DATA P      / 1.00,     1.00,     1.00     /
      DATA A      / 0.031091, 0.015545, 0.016887 /
      DATA ALPHA1 / 0.21370,  0.20548,  0.11125  /
      DATA BETA   / 7.5957,  14.1189,  10.357,
     .              3.5876,   6.1977,   3.6231,
     .              1.6382,   3.3662,   0.88026,
     .              0.49294,  0.62517,  0.49671 /

C Find rs and zeta
      PI = 4 * ATAN(1.D0)
      IF (NSPIN .EQ. 1) THEN
        DTOT = MAX( DENMIN, DENS(1) )
        ZETA = 0
        RS = ( 3 / (4*PI*DTOT) )**THD
C       Find derivatives dRs/dDens and dZeta/dDens
        DRSDD = - RS / DTOT / 3
        DZDD(1) = 0
      ELSE
        DTOT = MAX( DENMIN, DENS(1)+DENS(2) )
        ZETA = ( DENS(1) - DENS(2) ) / DTOT
        RS = ( 3 / (4*PI*DTOT) )**THD
        DRSDD = - RS / DTOT / 3
        DZDD(1) =   1 / DTOT - ZETA / DTOT
        DZDD(2) = - 1 / DTOT - ZETA / DTOT
      ENDIF

C Find eps_c(rs,0)=G(0), eps_c(rs,1)=G(1) and -alpha_c(rs)=G(2)
C using eq.(10) of cited reference (Perdew & Wang, PRB, 45, 13244 (92))
      DO 20 IG = 0,2
        B = BETA(IG,1) * RS**HALF   +
     .      BETA(IG,2) * RS         +
     .      BETA(IG,3) * RS**THRHLF +
     .      BETA(IG,4) * RS**(P(IG)+1)
        DBDRS = BETA(IG,1) * HALF      / RS**HALF +
     .          BETA(IG,2)                         +
     .          BETA(IG,3) * THRHLF    * RS**HALF +
     .          BETA(IG,4) * (P(IG)+1) * RS**P(IG)
        C = 1 + 1 / (2 * A(IG) * B)
        DCDRS = - (C-1) * DBDRS / B
        G(IG) = - 2 * A(IG) * ( 1 + ALPHA1(IG)*RS ) * LOG(C)
        DGDRS(IG) = - 2*A(IG) * ( ALPHA1(IG) * LOG(C) +
     .                            (1+ALPHA1(IG)*RS) * DCDRS / C )
   20 CONTINUE

C Find f''(0) and f(zeta) from eq.(9)
      C = 1 / (2**FOUTHD - 2)
      FPP0 = 8 * C / 9
      F = ( (1+ZETA)**FOUTHD + (1-ZETA)**FOUTHD - 2 ) * C
      DFDZ = FOUTHD * ( (1+ZETA)**THD - (1-ZETA)**THD ) * C

C Find eps_c(rs,zeta) from eq.(8)
      EC = G(0) - G(2) * F / FPP0 * (1-ZETA**4) +
     .    (G(1)-G(0)) * F * ZETA**4
      DECDRS = DGDRS(0) - DGDRS(2) * F / FPP0 * (1-ZETA**4) +
     .        (DGDRS(1)-DGDRS(0)) * F * ZETA**4
      DECDZ = - G(2) / FPP0 * ( DFDZ*(1-ZETA**4) - F*4*ZETA**3 ) +
     .        (G(1)-G(0)) * ( DFDZ*ZETA**4 + F*4*ZETA**3 )

C Find correlation potential
      IF (NSPIN .EQ. 1) THEN
        DECDD(1) = DECDRS * DRSDD
        VC(1) = EC + DTOT * DECDD(1)
      ELSE
        DECDD(1) = DECDRS * DRSDD + DECDZ * DZDD(1)
        DECDD(2) = DECDRS * DRSDD + DECDZ * DZDD(2)
        VC(1) = EC + DTOT * DECDD(1)
        VC(2) = EC + DTOT * DECDD(2)
      ENDIF

      END

ext: f

max_stars_repo_name: qsnake/petot

max_issues_repo_name: qsnake/petot

max_forks_repo_name: qsnake/petot


text:
CommuniCare, formerly known as the Davis Community Clinic, provides primary care medical services to the lowincome and Uninsured uninsured population of Davis. Services include laboratory, pediatric, tobacco prevention, reproductive health services, prenatal care, CHDP exams, childhood immunizations, HIV testing, family planning programs, dental clinic, teen clinic, Youth for Recovery drug treatment program and vision care. Rates are on a sliding scale and insurance plans accepted include MediCal, Healthy Families, Partnership Health Plan, and http://www.yolocounty.org/Index.aspx?page67 YCHIP.

Patients who have no form of health coverage are asked to pay for services provided according to a sliding fee scale. The fee for a doctor visit ranges from $30 to $140.

Volunteer opportunities at the seven CommuniCare health centers include work in medical records, the lab, and language translation.

See also: Medical Services, CommuniCare HC Teen Clinic

20080922 14:30:46 nbsp I used to (and still do on occasion) refer clients here.  They offer quality services.  And the Staff seems to truly care about their patients.   Users/CurlyGirl26

20080922 19:56:03 nbsp I agree with that comment.  Users/TaniaG

20100224 22:50:59 nbsp Great Place! :) gave me help when i really needed it :)  Users/chrisqtan

20100408 08:20:31 nbsp I make 2,000 dollars/month and apparently with that kind of income they charge 180 for an office visit.  Thats about the same as Sutter. Users/AmandaAbughosh
  Are you sure that Sutters prices arent at least double or triple that?  Its been years since Ive been treated at Sutter, but I think I got charged something like $450 for a fiveminute appointment.  Users/ScottMeehleib

  Yes, I am sure that Sutters prices arent at least double or triple that.  The base office visit is 180.  Of course, if they take tests or anything, Im sure it would be more.  Users/AmandaAbughosh

20101103 16:16:22 nbsp Patients who have no form of health coverage are asked to pay for services provided according to a sliding fee scale. The minimum amount on the sliding scale is $30 for an office visit, not $180. You can confirm this for yourself by calling the clinic. Users/JohnBonner

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


Program final2
	Implicit none
	real*8 :: infty, x, dx, e1, e2, e, tol, phi1, phi2
	real*8, dimension(2) :: y1, y2, y
	real*8, external :: f_x, g_x
	character (len=10) :: filename
	integer :: n=1

! parameters
	tol = epsilon(0.0d0) ! tolerance of 0
	infty = 2.4d0 ! where psi(infty) is considered to be tol
	e1 = -8.0d0  ! lower bound energy level

	do while ( n.le.4 ) ! get first 4 lowest energy

		x = -infty
		dx = infty/1.0d4
		y1(1) = tol
		y1(2) = tol
		do while( x.le.infty ) ! integrate to get psi1(infty)
			call rk4(x, dx, f_x, g_x, y1, e1)
			x = x + dx
		enddo
		phi1 = y1(1)
		phi2 = y1(1)
		e2 = e1

		print *, "e1 and phi1", e1, phi1

		do while ( phi1*phi2.gt.0 ) ! find e2
			x = -infty
			e2 = e2 + 3.d0/7.d0
			dx = infty/1.0d5
			y2(1) = tol
			y2(2) = tol

			do while( x.le.infty ) ! integrate to get psi2(infty)
				call rk4(x, dx, f_x, g_x, y2, e2)
				x = x + dx
			enddo
			phi2 = y2(1)

			print *, "finding e2", e2, phi2
		enddo

		y(1) = 1.0
		do while (abs(phi1-phi2).gt.tol) ! find e by updating e1 & e2
			x = -infty
			e = (e1 + e2)/2.d0
			dx = ( e2 - e1 )/9.d0
			y1(1) = tol
			y2(1) = tol
			y(1) = tol
			y1(2) = tol
			y2(2) = tol
			y(2) = tol

			write(filename , '("res",i1,".dat")') n
			open(10, file=filename,action="write",status="replace")
			do while( x.le.infty )

				call rk4(x, dx, f_x, g_x, y1, e1)
				call rk4(x, dx, f_x, g_x, y2, e2)
				call rk4(x, dx, f_x, g_x, y, e)
				write (10,*) x, y(1), y(2)
				x = x + dx
			enddo
			phi1 = y1(1)
			phi2 = y2(1)

			if (phi1*y(1).lt.0) then
				e2 = e
			else
				e1 = e
			endif

			close(10)

!			print *, "updating e1, e2", e1, e2, y(1)
		enddo
		e1 = e + 0.1d0/7.d0

		print *, "eigenvalue",n,"is:", e
		n = n + 1
	enddo

end program final2

! 4th-order Runge-Kutta subroutine
subroutine rk4(x, dx, df, dg, y, e)
	implicit none
	real*8, external :: df, dg
	real*8, intent(in) :: x, dx, e
	real*8, intent(inout), dimension(2) :: y
	real*8 :: h, k0, k1, k2, k3, l0, l1, l2, l3

	h=dx/2.0D0

	k0 = dx * df(y(1),y(2))
	l0 = dx * dg(y(1),y(2),x,e)
	k1 = dx * df(y(1)+h, y(2)+0.5d0*l0)
	l1 = dx * dg(y(1)+0.5d0*k0,y(2)+0.5d0*l0,x,e)
	k2 = dx * df(y(1)+0.5d0*k1,y(2)+0.5d0*l1)
	l2 = dx * dg(y(1)+0.5d0*k1,y(2)+0.5d0*l1,x,e)
	k3 = dx * df(y(1)+k2,y(2)+l2)
	l3 = dx * dg(y(1)+k2,y(2)+l2,x,e)
	y(1) = y(1) + (k0+2*k1+2*k2+k3)/6.0d0
	y(2) = y(2) + (l0+2*l1+2*l2+l3)/6.0d0
	Return
End subroutine rk4

! function which returns the derivatives (RHS)
real*8 function f_x(y, yp)
! dpsi/dx = psi'
	Implicit none
	real*8 ,intent(in) :: y, yp
	f_x = yp
	Return
End function f_x

real*8 function g_x(y, yp, x, e)
! d^2y/dx^2 = 2(V(x)-e)*y = (4x^4-16x^2-2e)y
	implicit none
	real*8, intent(in) :: y, yp, x, e
!	g_x = ( -4.0d0 - e) * y
	g_x = ( 4.0d0 * x**4.0d0 - 16.0d0 * x**2.0d0 - 2.0d0*e) * y
	Return
End function g_x


!
!
!                             MLD2P4  version 2.0
!    MultiLevel Domain Decomposition Parallel Preconditioners Package
!               based on PSBLAS (Parallel Sparse BLAS version 3.0)
!
!    (C) Copyright 2008,2009,2010
!
!                        Salvatore Filippone
!                        Alfredo Buttari
!                        Pasqua D'Ambra
!                        Daniela di Serafino
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the MLD2P4 group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
!
module data_input

  interface read_data
    module procedure read_char, read_int,&
         & read_double, read_single, read_logical,&
         & string_read_char, string_read_int,&
         & string_read_double, string_read_single, &
         & string_read_logical
  end interface read_data
  interface trim_string
    module procedure trim_string
  end interface

  character(len=4096), private  :: charbuf
  character, private, parameter :: def_marker="!"

contains

  subroutine read_logical(val,file,marker)
    logical, intent(out) :: val
    integer, intent(in)           :: file
    character(len=1), optional, intent(in) :: marker

    read(file,'(a)')charbuf
    call read_data(val,charbuf,marker)

  end subroutine read_logical

  subroutine read_char(val,file,marker)
    character(len=*), intent(out) :: val
    integer, intent(in)           :: file
    character(len=1), optional, intent(in) :: marker

    read(file,'(a)')charbuf
    call read_data(val,charbuf,marker)

  end subroutine read_char

  subroutine read_int(val,file,marker)
    integer, intent(out) :: val
    integer, intent(in)  :: file
    character(len=1), optional, intent(in) :: marker

    read(file,'(a)')charbuf
    call read_data(val,charbuf,marker)

  end subroutine read_int
  subroutine read_single(val,file,marker)
    use psb_base_mod
    real(psb_spk_), intent(out) :: val
    integer, intent(in)         :: file
    character(len=1), optional, intent(in) :: marker

    read(file,'(a)')charbuf
    call read_data(val,charbuf,marker)

  end subroutine read_single
  subroutine read_double(val,file,marker)
    use psb_base_mod
    real(psb_dpk_), intent(out) :: val
    integer, intent(in)         :: file
    character(len=1), optional, intent(in) :: marker

    read(file,'(a)')charbuf
    call read_data(val,charbuf,marker)

  end subroutine read_double

  subroutine string_read_char(val,file,marker)
    character(len=*), intent(out) :: val
    character(len=*), intent(in)  :: file
    character(len=1), optional, intent(in) :: marker
    character(len=1)    :: marker_
    character(len=1024) :: charbuf
    integer :: idx
    if (present(marker)) then
      marker_ = marker
    else
      marker_ = def_marker
    end if
    read(file,'(a)')charbuf
    charbuf = adjustl(charbuf)
    idx=index(charbuf,marker_)
    if (idx == 0) idx = len(charbuf)+1
    read(charbuf(1:idx-1),'(a)') val
  end subroutine string_read_char

  subroutine string_read_int(val,file,marker)
    integer, intent(out) :: val
    character(len=*), intent(in)  :: file
    character(len=1), optional, intent(in) :: marker
    character(len=1)    :: marker_
    character(len=1024) :: charbuf
    integer :: idx
    if (present(marker)) then
      marker_ = marker
    else
      marker_ = def_marker
    end if
    read(file,'(a)')charbuf
    charbuf = adjustl(charbuf)
    idx=index(charbuf,marker_)
    if (idx == 0) idx = len(charbuf)+1
    read(charbuf(1:idx-1),*) val
  end subroutine string_read_int

  subroutine string_read_single(val,file,marker)
    use psb_base_mod
    real(psb_spk_), intent(out) :: val
    character(len=*), intent(in)         :: file
    character(len=1), optional, intent(in) :: marker
    character(len=1)    :: marker_
    character(len=1024) :: charbuf
    integer :: idx
    if (present(marker)) then
      marker_ = marker
    else
      marker_ = def_marker
    end if
    read(file,'(a)')charbuf
    charbuf = adjustl(charbuf)
    idx=index(charbuf,marker_)
    if (idx == 0) idx = len(charbuf)+1
    read(charbuf(1:idx-1),*) val
  end subroutine string_read_single

  subroutine string_read_double(val,file,marker)
    use psb_base_mod
    real(psb_dpk_), intent(out) :: val
    character(len=*), intent(in)         :: file
    character(len=1), optional, intent(in) :: marker
    character(len=1)    :: marker_
    character(len=1024) :: charbuf
    integer :: idx
    if (present(marker)) then
      marker_ = marker
    else
      marker_ = def_marker
    end if
    read(file,'(a)')charbuf
    charbuf = adjustl(charbuf)
    idx=index(charbuf,marker_)
    if (idx == 0) idx = len(charbuf)+1
    read(charbuf(1:idx-1),*) val
  end subroutine string_read_double

  subroutine string_read_logical(val,file,marker)
    use psb_base_mod
    logical, intent(out) :: val
    character(len=*), intent(in)         :: file
    character(len=1), optional, intent(in) :: marker
    character(len=1)    :: marker_
    character(len=1024) :: charbuf
    integer :: idx
    if (present(marker)) then
      marker_ = marker
    else
      marker_ = def_marker
    end if
    read(file,'(a)')charbuf
    charbuf = adjustl(charbuf)
    idx=index(charbuf,marker_)
    if (idx == 0) idx = len(charbuf)+1
    read(charbuf(1:idx-1),*) val
  end subroutine string_read_logical

  function  trim_string(string,marker)
    character(len=*), intent(in) :: string
    character(len=1), optional, intent(in) :: marker
    character(len=len(string))    :: trim_string
    character(len=1)              :: marker_
    integer :: idx
    if (present(marker)) then
      marker_ = marker
    else
      marker_ = def_marker
    end if
    idx=index(string,marker_)
    trim_string = adjustl(string(idx:))
  end function trim_string
end module data_input


text:
Address(Crystal Grove Lane) is a residential street in West Davis. The Bike Overpasses bike overcrossing from Sycamore Park crosses Highway 113 here.

 Intersecting Streets

  Santa Rosa Street
  Cape Cod Street
  Hampton Drive


text:
program associate_03
    implicit none

    integer, pointer :: p1
    integer, target :: t1 = 2, t2 = 1
    integer :: i

    print *, t1, t2
    ! Runtime pointer association
    if (t1 > t2) then
        p1 => t1
    else
        p1 => t2
    end if
    print *, p1
    ! Does not work yet:
    !if (p1 == t2) error stop
    i = p1
    if (i == t2) error stop

end program

max_stars_repo_name: Thirumalai-Shaktivel/lfortran

max_issues_repo_name: Thirumalai-Shaktivel/lfortran

max_forks_repo_name: Thirumalai-Shaktivel/lfortran


text:
module mod_monolis_neighbor_search
  use mod_monolis_prm
  use mod_monolis_stdlib

  private

  public :: monolis_neighbor_search_init
  public :: monolis_neighbor_search_finalize
  public :: monolis_neighbor_search_get_by_position
  public :: monolis_neighbor_search_get_by_bb
  public :: monolis_neighbor_search_push
  public :: type_monolis_neighbor_search

  type type_monolis_neighbor_search_main
    integer(kint) :: nid = 0
    integer(kint), allocatable :: id(:)
  end type type_monolis_neighbor_search_main

  type type_monolis_neighbor_search
    integer(kint) :: div(3)
    real(kdouble) :: BB(6) !> xmin, xmax, ymin, ymax, zmin, zmax
    real(kdouble) :: dx(3) !> dx, dy, dz
    type(type_monolis_neighbor_search_main), allocatable :: cell(:)
  end type type_monolis_neighbor_search

  real(kdouble), save :: ths = 1.0d-3

contains

  subroutine monolis_neighbor_search_init(monolis_nbsearch, BB, div)
    implicit none
    type(type_monolis_neighbor_search) :: monolis_nbsearch
    integer(kint) :: div(3)
    real(kdouble) :: BB(6), minbb
    monolis_nbsearch%div = div
    monolis_nbsearch%BB = BB
    monolis_nbsearch%dx(1) = (BB(2)-BB(1))/dble(div(1))
    monolis_nbsearch%dx(2) = (BB(4)-BB(3))/dble(div(2))
    monolis_nbsearch%dx(3) = (BB(6)-BB(5))/dble(div(3))
    allocate(monolis_nbsearch%cell(div(1)*div(2)*div(3)))

    minbb = min(BB(2)-BB(1), BB(4)-BB(3), BB(6)-BB(5))
    ths = 1.0d-3*minbb
  end subroutine monolis_neighbor_search_init

  subroutine monolis_neighbor_search_push(monolis_nbsearch, BB, id)
    implicit none
    type(type_monolis_neighbor_search) :: monolis_nbsearch
    integer(kint) :: id, x, y, z, in, imin(3), imax(3), div(3)
    real(kdouble) :: BB(6), pos(3)

    pos(1) = BB(1) - ths
    pos(2) = BB(3) - ths
    pos(3) = BB(5) - ths
    call get_int_coordinate(monolis_nbsearch, pos, imin)

    pos(1) = BB(2) + ths
    pos(2) = BB(4) + ths
    pos(3) = BB(6) + ths
    call get_int_coordinate(monolis_nbsearch, pos, imax)

    div = monolis_nbsearch%div
    do z = imin(3), imax(3)
    do y = imin(2), imax(2)
    do x = imin(1), imax(1)
      in = get_index(monolis_nbsearch%div, x, y, z)
      call monolis_neighbor_search_push_main(monolis_nbsearch, in, id)
    enddo
    enddo
    enddo
  end subroutine monolis_neighbor_search_push

  subroutine monolis_neighbor_search_push_main(monolis_nbsearch, in, id)
    implicit none
    type(type_monolis_neighbor_search) :: monolis_nbsearch
    integer(kint) :: in, id, add(1), nid
    nid = monolis_nbsearch%cell(in)%nid
    add = id
    call monolis_reallocate_integer(monolis_nbsearch%cell(in)%id, nid, 1, add)
    monolis_nbsearch%cell(in)%nid = nid + 1
  end subroutine monolis_neighbor_search_push_main

  subroutine monolis_neighbor_search_get_by_position(monolis_nbsearch, pos, nid, id)
    implicit none
    type(type_monolis_neighbor_search) :: monolis_nbsearch
    integer(kint) :: nid, index
    integer(kint), allocatable :: id(:)
    real(kdouble) :: pos(3)
    logical :: is_in

    nid = 0
    if(allocated(id)) deallocate(id)
    call BB_check(monolis_nbsearch, pos, is_in)
    if(.not. is_in) return

    call get_index_by_position(monolis_nbsearch, pos, index)

    nid = monolis_nbsearch%cell(index)%nid
    if(nid > 0)then
      allocate(id(nid))
      id = monolis_nbsearch%cell(index)%id
    endif
  end subroutine monolis_neighbor_search_get_by_position

  subroutine monolis_neighbor_search_get_by_bb(monolis_nbsearch, BB, nid, id)
    implicit none
    type(type_monolis_neighbor_search) :: monolis_nbsearch
    integer(kint) :: nid, morton_id, imin(3), imax(3), div(3)
    integer(kint) :: in, x, y, z, newlen
    integer(kint), allocatable :: id(:), tmp(:)
    real(kdouble) :: BB(6), pos(3)
    logical :: is_in

    nid = 0
    if(allocated(id)) deallocate(id)
    call BB_modify(monolis_nbsearch, BB, is_in)
    if(.not. is_in) return

    pos(1) = BB(1) - ths
    pos(2) = BB(3) - ths
    pos(3) = BB(5) - ths
    call get_int_coordinate(monolis_nbsearch, pos, imin)

    pos(1) = BB(2) + ths
    pos(2) = BB(4) + ths
    pos(3) = BB(6) + ths
    call get_int_coordinate(monolis_nbsearch, pos, imax)

    div = monolis_nbsearch%div
    do z = imin(3), imax(3)
    do y = imin(2), imax(2)
    do x = imin(1), imax(1)
      in = get_index(monolis_nbsearch%div, x, y, z)
      call monolis_neighbor_search_get_by_bb_main(monolis_nbsearch, nid, in, tmp)
    enddo
    enddo
    enddo

    if(.not. allocated(tmp)) return

    call monolis_qsort_int(tmp, 1, nid)
    call monolis_uniq_int(tmp, nid, newlen)

    allocate(id(newlen))
    id = tmp(1:newlen)
    deallocate(tmp)
    nid = newlen
  end subroutine monolis_neighbor_search_get_by_bb

  subroutine monolis_neighbor_search_get_by_bb_main(monolis_nbsearch, nid, eid, id)
    implicit none
    type(type_monolis_neighbor_search) :: monolis_nbsearch
    integer(kint) :: eid, ntmp, nid
    integer(kint), allocatable :: id(:)

    ntmp = monolis_nbsearch%cell(eid)%nid
    if(ntmp > 0)then
      call monolis_reallocate_integer(id, nid, ntmp, monolis_nbsearch%cell(eid)%id)
      nid = nid + ntmp
    endif
  end subroutine monolis_neighbor_search_get_by_bb_main

  subroutine monolis_neighbor_search_finalize(monolis_nbsearch)
    implicit none
    type(type_monolis_neighbor_search) :: monolis_nbsearch
    deallocate(monolis_nbsearch%cell)
  end subroutine monolis_neighbor_search_finalize

  subroutine BB_check(monolis_nbsearch, pos, is_in)
    implicit none
    type(type_monolis_neighbor_search) :: monolis_nbsearch
    real(kdouble) :: pos(3), BB(6)
    logical :: is_in

    is_in = .false.
    BB = monolis_nbsearch%BB
    if( BB(1)-ths <= pos(1) .and. pos(1) <= BB(2)+ths .and. &
      & BB(3)-ths <= pos(2) .and. pos(2) <= BB(4)+ths .and. &
      & BB(5)-ths <= pos(3) .and. pos(3) <= BB(6)+ths )then
      is_in = .true.
    endif
  end subroutine BB_check

  subroutine BB_modify(monolis_nbsearch, BB_in, is_in)
    implicit none
    type(type_monolis_neighbor_search) :: monolis_nbsearch
    integer(kint) :: i
    real(kdouble) :: pos(3), BB(6), BB_in(6)
    logical :: is_in
    integer(kint), parameter :: perm(3,8) = reshape([ &
      1, 3, 5, &
      2, 3, 5, &
      1, 4, 5, &
      2, 4, 5, &
      1, 3, 6, &
      2, 3, 6, &
      1, 4, 6, &
      2, 4, 6  &
      ], [3,8])

    is_in = .true.
    BB = monolis_nbsearch%BB

    if(BB_in(2) < BB(1)) is_in = .false.
    if(BB_in(4) < BB(3)) is_in = .false.
    if(BB_in(6) < BB(5)) is_in = .false.
    if(BB_in(1) > BB(2)) is_in = .false.
    if(BB_in(3) > BB(4)) is_in = .false.
    if(BB_in(5) > BB(6)) is_in = .false.

    if(.not. is_in) return

    if(BB_in(1) < BB(1)) BB_in(1) = BB(1)
    if(BB_in(3) < BB(3)) BB_in(3) = BB(3)
    if(BB_in(5) < BB(5)) BB_in(5) = BB(5)
    if(BB_in(2) > BB(2)) BB_in(2) = BB(2)
    if(BB_in(4) > BB(4)) BB_in(4) = BB(4)
    if(BB_in(6) > BB(6)) BB_in(6) = BB(6)
  end subroutine BB_modify

  subroutine get_index_by_position(monolis_nbsearch, pos, index)
    implicit none
    type(type_monolis_neighbor_search) :: monolis_nbsearch
    integer(kint) :: index, id(3)
    real(kdouble) :: pos(3)

    call get_int_coordinate(monolis_nbsearch, pos, id)
    index = get_index(monolis_nbsearch%div, id(1), id(2), id(3))
  end subroutine get_index_by_position

  function get_index(div, x, y, z)
    implicit none
    integer(kint) :: get_index, x, y, z, div(3)
    get_index = x + (y-1)*div(1) + (z-1)*div(1)*div(2)
  end function get_index

  subroutine get_int_coordinate(monolis_nbsearch, pos, id)
    implicit none
    type(type_monolis_neighbor_search) :: monolis_nbsearch
    integer(kint) :: id(3)
    real(kdouble) :: pos(3)

    id(1) = int((pos(1) - monolis_nbsearch%BB(1))/monolis_nbsearch%dx(1)) + 1
    id(2) = int((pos(2) - monolis_nbsearch%BB(3))/monolis_nbsearch%dx(2)) + 1
    id(3) = int((pos(3) - monolis_nbsearch%BB(5))/monolis_nbsearch%dx(3)) + 1
    if(id(1) < 1) id(1) = 1
    if(id(2) < 1) id(2) = 1
    if(id(3) < 1) id(3) = 1
    if(id(1) > monolis_nbsearch%div(1)) id(1) = monolis_nbsearch%div(1)
    if(id(2) > monolis_nbsearch%div(2)) id(2) = monolis_nbsearch%div(2)
    if(id(3) > monolis_nbsearch%div(3)) id(3) = monolis_nbsearch%div(3)
  end subroutine get_int_coordinate
end module mod_monolis_neighbor_search

max_stars_repo_name: nqomorita/monolis

max_issues_repo_name: nqomorita/monolis

max_forks_repo_name: nqomorita/monolis


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! NAME:         bestpred_fmt4.f90
! VERSION:      2.0 beta
! RELEASED:     01 AUGUST 2007
! AUTHORS:      Paul M. VanRaden (paul@aipl.arsusda.gov)
!               John B. Cole (john.cole@ars.usda.gov)
! DESCRIPTION:  Prepares format 4 data for subroutine aiplDCR. This program
!               is part of the bestpred package from AIPL.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    subroutine bestpred_fmt4( &
!                                                                 Input
          format4,doprev,maxprnt &
!                                                                Output
          ,DCRm,DCRc,DCRs,YLDvec,PERSvec,RELyld,RELpers &
!                                                              Optional
          ,Yvec,herd305,bump,BLUPout,BLUPn              &
          ,bpfreqin,GRAFplot,DEBUGmsgs,ONscreen         &
          ,mtrait,use3X                                 &
          ,laclen,dailyfreq,INTmethod,maxlen            &
          ,DAILYbp,DAILYherd,WRITEcurve,CURVEfile       &
          ,WRITEdata,DATAfile,plotfreq, INTmethodSCS    &
          ,READparms,UNITSin, UNITSout, breedUNK, dim0  &
          ,dim0flag, LOGon, LOGfile, LOGfreq,maxshow    &
          ,CURVEsmall, CURVEsingle, region, season      &
          )
! .....................................................................
!   Definitions (MT=multi-trait,ST=single-trait,m=milk,f=fat,p=protein)
!                                                                 Input
!       format4 = AIPL format 4 record (710 bytes) as defined
!                 on aipl.arsusda.gov web site
!       doprev  = days not pregnant in previous lactation
!       maxprnt = number of graphs to display (must not be negative)
!               = -1 displays only the input variables
!                                                                Output
!       DCRm   = MT m DCR, DCRc = MT f,p DCR,  DCRs = MT scs DCR
!       YLDvec = MT m,f,p,scs, ST m,f,p,scs, non-305 m,f,p,scs
!       PERSvec= MT m,f,p,scs persistency, ST m,f,p,scs persistency
!       RELyld = MT m,f,p,scs Reliability, ST m,f,p,scs Reliability
!       RELpers= MT m,f,p,scs Rel(persist),ST m,f,p,scs Rel(persist)
!       length = days from fresh date to last test or cull date
!                                                              Optional
!       GRAFplot = Array of switches toggling lactation curve graphs
!                  on (1) or off (0); the order of traits is: M,F,P,S.
!       DEBUGmsgs = Used to display debugging messages.
!       ONscreen = toggle output to the screen on (1) and off (0).
!       laclen specifies default lactation length
!       dailyfreq specifies how often actual dailies are calculated
!       INTmethod specifies the inerpolation method to be used for
!           calculating lactation curves.
!       READparms is a flag passed downstream from maindcr to fmt4dcr
!           and then to aipldcr that indicates which parms should be
!           used by that particular subroutine -- those passed in from
!           maindcr (0) or those stored in the file maindcr.par (1).
! .....................................................................
      integer    mtrait,  use3X,  maxtd, READparms
      parameter (maxtd=50)
!                                            aiplDCR control variables:
!                         mtrait = 1 does only ST to save CPU time
!                                = 3 does MT m,f,p; 4 does MT m,f,p,scs
!                          use3X = 0 doesn't adjust for 3X milking,
!                                  1 uses old factors, 2 new factors,
!                                  3 uses phased-in factors over time
      integer*2 agebase
!                        agebase = 0 uses mature equivalent age base
!                                =36 uses 36 month age base, etc.
      integer   showage
!                        showage = 1 displays age factors at startup
!                                = 0 doesn't
! ....................................................................
      !character*710 :: format4
      character*1400 :: format4
      character*100 :: BLUPout(10)
      character*23  :: segment(maxtd)
      character*17  :: cowid
      character*8   :: fresh,herd,birth
      character     ::  brd
      character(len=*)  :: CURVEfile, DATAfile, LOGfile
      real*8 DCRm,DCRc,DCRs,YLDvec(2,16),PERSvec(4),RELyld(4),RELpers(4), &
             yield(2,4,maxtd),Yvec(2,4),bump(4),agefac(4),avgfac(4),        &
             herd305(2,4),me305(3),dev305(3),zero
      real*8 oldfac(4)
      integer ntd,parity,maxprnt,length,BLUPn(10),                      &
              dim(maxtd),Xmilk(maxtd),weigh(maxtd),sample(maxtd),       &
              MRD(maxtd),super(maxtd),status(maxtd),                    &
              adjscs4,i,WRITEcurve,WRITEdata,LOGon,LOGfreq,maxshow
      integer   count, CURVEsmall, CURVEsingle
      integer*2 pday,age,doprev,docur,lacno,frmo,fryr,state         &
               ,fill2,yld2(3),adjyld2(3),i305,scs305
!                                                 Simple test day means
      real :: ybar(4)
!                               User-specified multiple trait parameter
      integer bpfreqin, GRAFplot(4), DEBUGmsgs, ONscreen, laclen        &
        , dailyfreq, maxlen, plotfreq, breedUNK
      character INTmethod, INTmethodSCS, UNITSin, UNITSout

      ! DAILYbp and DAILYherd will contain daily BP of yield for cows
      ! and herds, respectively, unless dailyfreq=0. If dailyfreq is
      ! 0 then they will contain only zeroes.
      real*8 DAILYbp(4,999)
      real*8 DAILYherd(4,999)
      character*2 :: herdstate

      ! integer :: breed11
      integer :: dim0(8), dim0flag, j, k
      integer :: region, season

      external pday

      data ybar    /70.,2.5,2.2,3.2/
      data count   /0/
      data docur   /90/
      data fill2   /0/
      data i305    /305/
      data agebase /0/
      data showage /0/
!      data showage /1/

!      print *, '[bestpred_fmt4]: CURVEsingle = ', CURVEsingle

! ....................................................................
      ! Make sure that mtrait has a valid value. 1 is used as a defauly
      ! because that's its value in PVR's original code, where mtrait
      ! had a hard-coded value.
      if (mtrait == 1 .or. mtrait == 3 .or. mtrait == 4) then
        continue
      else
        mtrait = 1
      end if
!      print *,'Frequency of BP: ', bpfreqin

      zero = 0.

      !herd305(:,:) = 0.d0

      do 5 j=1,3
 5      yld2(j) = 10000
!                                                 Display age factors
      count = count + 1
      if(count == 1 .and. showage == 1) then
        do i=1,6
          if (i == 1) brd = 'H'
          if (i == 2) brd = 'J'
          if (i == 3) brd = 'B'
          if (i == 4) brd = 'G'
          if (i == 5) brd = 'A'
          if (i == 6) brd = 'M'
          if (i == 1 ) print '(26x,a)','Milk    Fat     Prot    SCS   prevDO'
          do fryr = 1990,1966,-6
            !print *,'Brd = ',brd,', Fresh year = ',fryr
            do age=22,99,14
              lacno = (age - 6)/14
              parity = lacno
              if(doprev < 1) doprev = 140
              if(lacno == 1) doprev = 0
              avgfac = 0.
              do frmo = 1,12
                do k=1,5
                  if (k == 1) state = 21
                  if (k == 2) state = 35
                  if (k == 3) state = 42
                  if (k == 4) state = 74
                  if (k == 5) state = 93
                  scs305 = 329
                  call aiplage(brd,age,fryr,frmo,lacno,state,doprev, &
                       agebase,agefac(1),agefac(2),agefac(3))
                  !print *, "Calling adjscs4()"
                  agefac(4) = adjscs4(brd,lacno,frmo,state, &
                       age,i305,scs305)/(scs305*1.d0)
                  do j=1,4
                    if(j < 4) oldfac(j) = adjyld2(j) / (yld2(j) * 1.d0)
                    avgfac(j) = avgfac(j) + agefac(j) / 60.
                    end do
                  end do
               print '(a,2i3,a,4f8.2,i8)','Age-par',age,lacno,' factors ',agefac,doprev
               print '(a,2i3,a,4f8.2,i8)','Age-par',age,lacno,' factors ',oldfac,doprev
!               print *
                end do
              print '(a,2i3,a,4f8.2,i8)','Age-par',age,lacno,' factors ',avgfac,doprev
              end do
            end do
          end do
        end if
!.......................................................... Read record
      read(format4,11) cowid,birth,herd,fresh,length,lacno,doprev,ntd, &
                      (segment(i),i=1,ntd)
! 11   format(t3,a17,t71,a8,t107,a8,t128,a8 &
!            ,t136,i3,t159,i2,t249,i2,50a23)
! Added to support previous days open (thanks to BJH).
 11   format(t3,a17,t71,a8,t107,a8,t128,a8 &
        ,t136,i3,t159,i2,t246,i3,t249,i2,50a23)
      !print *, cowid, " ", birth, " ", herd, " ", fresh, " ", length, " ", lacno, " ", ntd, " ", doprev
!                                                Read test day segments
      do 16 i=1,ntd
        read(segment(i),13) dim(i),super(i),status(i), &
                   Xmilk(i),weigh(i),sample(i),MRD(i), &
                   (yield(1,k,i),k=1,4)
 13     format(i3,5i1,i2,3x,f4.1,3f2.1)
        !print *, "    ", dim(i), super(i), status(i), Xmilk(i), weigh(i), sample(i), MRD(i), yield(1,1,i) &
        !  , yield(1,2,i), yield(1,3,i), yield(1,4,i)
        if ( status(i) == 3 ) then
          do 14 k=1,4
 14         yield(1,k,i) = zero
        end if
        if(status(i) == 2) then
          yield(1,2,i) = zero
          yield(1,3,i) = zero
          ! JBC 08/11/2009 Commented this out so that records with SCS but no
          ! fat or protein test make it through to bestpred().
          !yield(1,4,i) = zero
        end if
!        if ( yield(1,4,i) == 0.d0 ) then
!            yield(1,4,i) = -0.01
!        end if
 16     continue

! Read the state code stored in bytes 107-108 of Format 4
      herdstate = herd(1:2)
      scs305 = 329
      if ( DEBUGmsgs > 0 ) print *, '[bestpred_fmt4] herdstate: ', herdstate
!                                       Read cow's deviation and me305
!      print *, "Format 4: ", format4
      read(format4,17) me305(1),me305(2),me305(3), &
                       dev305(1),dev305(2),dev305(3)
 17   format(t188,f5.0,2f4.0,f6.0,2f5.0)
      if ( DEBUGmsgs > 0 ) print *, '[bestpred_fmt4]: herd305 before bestpred_fmt4()', herd305
!                                    Backsolve to get mature herd mean
!      print *, '[bestpred_fmt4]: herd305 before back-calculation in fmt4', herd305
!      print *, '[bestpred_fmt4]: herd305 before zero check in bestpred_fmt4()', herd305(1,:)
      do j = 1,4
!        herd305(j) = me305(j) - dev305(j)
! I added this IF to help with Albert's work on using some sort of PPA by
! increasing the herd average for a later-lactation cow to reflect the
! repeatability of her earlier performance. The herd average for a trait
! is only back-calculated from me305 and dev305 if the herd average in
! herd305 is <= 0.
        if ( herd305(1,j) <= 0 ) then
          if ( j < 4 ) then
            herd305(1,j) = me305(j) - dev305(j)
            if ( me305(j) == 0.d0 ) herd305(1,j) = 0.d0
          else
            herd305(1,j) = 0.d0
            !herd305(1,4) = scs305
          end if
        end if
      end do
      herd305(2,:) = herd305(1,:)
!      print *, '[bestpred_fmt4]: herd305 after back-calculation in fmt4', herd305

!      print *, '[bestpred_fmt4]: herd305 after zero check in bestpred_fmt4()', herd305(1,:)
!      herd305(4) = 0.d0
      if(birth == '        ') birth = '19500101'
      if(birth  < '19500101') birth = '19500101'
      age = (pday(fresh) - pday(birth))/30.5
      if(fresh == '        ') age = 84
!                                         Assign missing lactation num
      if(lacno <= 0) lacno = max(1,(age-6)/13)
      if(birth == '19500101') age = 12 + 15*lacno
      read(fresh(1:4),'(i4)') fryr
      read(fresh(5:6),'(i2)') frmo
      read(herd(1:2),'(i2)') state
      if(frmo < 1) frmo = 1
      brd = cowid(1:1)
      if(index('ABGHJM',brd) == 0) brd = 'H'
      parity = lacno
      if(doprev < 1) doprev = 140
      if(lacno == 1) doprev = 0
      !scs305 = 329
!                                    Get factor to adjust 305-d records
!          ageadj4                               is official C function
!          ageadjF                       is substitute Fortran function
!          aiplage             also adjusts pre-1990 data by C function
!     call ageadj4(brd,lacno,frmo,state,age, &
!                  fill2,doprev,docur,yld2,adjyld2)

      if ( DEBUGmsgs > 1 ) then
          print *, 'Age factors before call to aiplage(): ', agefac
          print *, 'aiplage() parameters'
          print *, '--------------------'
          print *, '    brd     : ', brd
          print *, '    age     : ', age
          print *, '    fryr    : ', fryr
          print *, '    frmo    : ', frmo
          print *, '    lacno   : ', lacno
          print *, '    state   : ', state
          print *, '    doprev  : ', doprev
          print *, '    agebase : ', agebase
      end if

      ! If the state code is empty/zero then default to 35, which is the
      ! code for Wisconsin. This prevents a nasty little side-effect with
      ! source 10 records (partial Format 4 records).
      if ( state == 0 ) then
          print *, '[bestpred_fmt4]: Changing state from ', &
              state, ' to 35 (Wisconsin).'
          state = 35
      end if

      call aiplage(brd,age,fryr,frmo,lacno,state,doprev, &
                   agebase,agefac(1),agefac(2),agefac(3))
!                                        Get factor to adjust 305-d SCS
!                 adjscs4                        is official C function
!                 adjscsF                is substitute Fortran function
      !agefac(4) = adjscs4(brd,lacno,frmo,state, &
      !            age,i305,scs305)/(scs305*1.d0)
      !print *, "Calling adjscs4()"
      agefac(4) = adjscs4(brd,lacno,frmo,state, &
                  age,i305,scs305)/(scs305*1.d0)
!      print *, 'Age factors after call to aiplage() :', agefac

      ! JBC 08/11/2008 -- We need to divide herd average SCS by 100 so that it is on the proper
      ! scale before we call bestpred(). The edits system in the routine run processing sends
      ! in SCS in X.XX format, not (X.XX*100). If we do the dividion in bestpred() we end
      ! up with herd average SCS near 0, which results in a systematic underestimation of
      ! ME SCS, and is particularly pronounced for early RIPs.
      herd305(:,4) = herd305(:,4) / 100

!      print *, '[bestpred_fmt4]: length in fmt4 before entry to bestpred: ', length
!      print *, '[bestpred_fmt4: maxprnt: ', maxprnt
!      print *, '[bestpred_fmt4: maxshow: ', maxshow
      !print *, "ntd before calling bestpred(): ", ntd
!      print *, '[bestpred_fmt4]: yield going into bestpred(): ', yield

      call bestpred(mtrait,use3X,maxprnt                        &
                  ,ntd,dim,super,Xmilk,weigh,sample,MRD,yield   &
                  ,herd305,agefac,cowid,fresh,parity,length     &
                  ,DCRm,DCRc,DCRs,YLDvec,PERSvec,RELyld,RELpers &
!                                                              Optional
                  ,Yvec,bump,BLUPout,BLUPn,GRAFplot             &
                  ,DEBUGmsgs,ONscreen,laclen,dailyfreq          &
                  ,INTmethod,maxlen,DAILYbp,DAILYherd           &
                  ,WRITEcurve,CURVEfile,WRITEdata,DATAfile      &
                  ,plotfreq,INTmethodSCS,READparms,UNITSin      &
                  ,UNITSout, breedUNK, dim0, dim0flag           &
                  ,LOGon,LOGfile,LOGfreq,maxshow                &
                  ,herdstate,CURVEsmall, CURVEsingle            &
                  ,region, season                               &
      )
!	  print *, '[bestpred_fmt4]: Leaving bestpred()'
      return
    end subroutine bestpred_fmt4
!----------------------------------------------------------------------
!                                 Fortran routine to standardize yield
!                                !!! Example factors for age,season !!!
!                                !!! Poor substitute for C function !!
!
    subroutine ageadjF(brd,lacno,frmo,state,age, &
                        fill2,doprev,docur,yld2,adjyld2)
!
      character brd*1
      integer*2 lacno,frmo,state,age,fill2,doprev,docur,yld2(3) &
               ,adjyld2(3)
      real agefac
      integer begin, j
      data begin   /0/

      if(begin == 0) then
        print *,'Age factors for example only, use C function instead'
        begin = 1
        end if
!                                                     Age adjustment
      if(lacno == 0) lacno = max(1,(age-6)/13)
      agefac = 1.0 + .0056*(72 - age)
      if(age > 72) agefac = 1.0
!                                                  Season adjustment
      agefac = agefac*(1.0 + .02*sin((frmo - 3)*3.1416/6.0))
!                                             Convert actual to mature
      do 20 j=1,3
 20     adjyld2(j) = yld2(j)*agefac
      return
    end subroutine ageadjF
!----------------------------------------------------------------------
!                              !!Example factors to standardize SCS !!!
!                                !!! Poor substitute for C function !!
!
    integer function adjscsF(brd,lacno,frmo,state,age,dim,scs305)
!
      character brd*1
      integer*2 lacno,frmo,state,age,dim,scs305
      real agefac
      integer begin
      data begin   /0/
      if(begin == 0) then
        print *,'SCS factors for example only, use C function instead'
        begin = 1
        end if
!                                     Adjustment to average age (45 mo)
      if(lacno == 0) lacno = max(1,(age)/13)
      agefac = 1.0 + .0094*(45 - age)
      if(age > 72) agefac = .75
!                                                     Season adjustment
      agefac = agefac*(1.0 + .04*sin((frmo - 9)*3.1416/6.0))
!                                                     Adjust input data
      adjscsF = agefac*scs305
      return
    end function adjscsF
! -------------------------- Compute perpetual days --------------------
!                                                  by Chris Schanefelter
!                                     Year 2000 changes by Paul VanRaden
    integer*2 function pday(date8)
      implicit none
      character  date8*8
      integer modays(12,2) &
        ,i,year,month,day,leap,balance_years        &
        ,number_of_years,year_groups

      data modays /31,28,31,30,31,30,31,31,30,31,30,31  &
                  ,31,29,31,30,31,30,31,31,30,31,30,31/
      data year   /0/
      data month  /0/
!                                                    Accepts 8-byte date
      read(date8,'(i4,2i2)') year,month,day
      year = year - 1900
!                                           0 represents an unknown date
!                                               set to December 31, 1889
      if (year <= 0) then
        pday = -21916
        return
      endif
      leap = year/4                  !     leap years are divisible by 4
      number_of_years = year - 60    !       number years from base year
      year_groups = abs(number_of_years) / 4    !          4 year groups
!                        Calculate number of years not in a 4 year group
      balance_years = abs(number_of_years) - (year_groups * 4)
      if(year .eq. (leap*4)) then
        leap = 2                                            !  leap year
      else
        leap = 1                                       ! not a leap year
      end if
!                                        Establish minimums and maximums
      if (month <= 0) month = 1
      if (month > 12) month = 12
      if (day .eq. 0) day = 1
      if (day .gt. modays(month,leap)) day = modays(month,leap)
      do i = 1, month-1
        day = day + modays(i,leap)
      end do
      if (number_of_years .le. -1 ) then !  before 1/1/1960
!         Adjust days when going away from zero in a negative direction
        day = (365 - day)
!                                              Calculate perpetual days
        pday = -((year_groups * 1461) +  day &
                + (balance_years - 1) * 365)
      else                               ! after and including 1/1/1960
!                For years not in a year group, add 1 leap and 365 days
!                                               for each year completed
        if (balance_years > 0) &
          balance_years = (((balance_years - 1)* 365)+ 366)

!                                              Calculate perpetual days
        pday = (year_groups * 1461) + day + balance_years
      end if
      pday = pday - 1
      return
    end function pday

max_stars_repo_name: stuvet/bestpred

max_issues_repo_name: stuvet/bestpred

max_forks_repo_name: stuvet/bestpred


        !COMPILER-GENERATED INTERFACE MODULE: Fri May 22 17:00:06 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE COMPUTE_ADD_SOURCES_VISCOELASTIC_GPU__genmod
          INTERFACE
            SUBROUTINE COMPUTE_ADD_SOURCES_VISCOELASTIC_GPU
            END SUBROUTINE COMPUTE_ADD_SOURCES_VISCOELASTIC_GPU
          END INTERFACE
        END MODULE COMPUTE_ADD_SOURCES_VISCOELASTIC_GPU__genmod


text:
C
C-----------------------------------------------------------------------------
C
      double precision function alfa(alfa0,qsq )
C
C-----------------------------------------------------------------------------
C
C	This function returns the 1-loop value of alpha.
C
C	INPUT:
C		qsq   = Q^2
C
C-----------------------------------------------------------------------------
C
      implicit none
      double precision  qsq,alfa0
c
c constants
c
      double precision  One, Three, Pi,zmass
      parameter( One = 1.0d0, Three = 3.0d0 )
      parameter( Pi = 3.14159265358979323846d0 )
      parameter( zmass = 91.188d0 )
cc
      alfa = alfa0 / ( 1.0d0 - alfa0*dlog( qsq/zmass**2 ) /Three /Pi )
ccc
      return
      end

C
C-----------------------------------------------------------------------------
C
      double precision function alfaw(alfaw0,qsq,nh )
C
C-----------------------------------------------------------------------------
C
C	This function returns the 1-loop value of alpha_w.
C
C	INPUT:
C		qsq = Q^2
C               nh  = # of Higgs doublets
C
C-----------------------------------------------------------------------------
C
      implicit none
      double precision  qsq, alphaw, dum,alfaw0
      integer  nh, nq
c
c	  include
c

c
c constants
c
      double precision  Two, Four, Pi, Twpi, zmass,tmass
      parameter( Two = 2.0d0, Four = 4.0d0 )
      parameter( Pi = 3.14159265358979323846d0 )
      parameter( Twpi = 3.0d0*Four*Pi )
      parameter( zmass = 91.188d0,tmass=174d0 )
cc
      if ( qsq.ge.tmass**2 ) then
         nq = 6
      else
         nq = 5
      end if
      dum = ( 22.0d0 - Four*nq - nh/Two ) / Twpi
      alfaw = alfaw0 / ( 1.0d0 + dum*alfaw0*dlog( qsq/zmass**2 ) )
ccc
      return
      end

      DOUBLE PRECISION FUNCTION ALPHAS(Q)

      DOUBLE PRECISION Q,ALPHAS_not_timed
      external ALPHAS_not_timed

c     timing statistics
      include "timing_variables.inc"

c     This function takes roughly 0.6 micro-seconds and the function
c     cpu_time takes 0.3 so it does not make much sense to always leave
c     the timing profile. But it can be uncommented for dedicated study.
c      call cpu_time(tbefore)
      ALPHAS=ALPHAS_not_timed(Q)
c      call cpu_time(tAfter)

c      tPDF = tPDF + (tAfter-tBefore)
      RETURN
      end

C-----------------------------------------------------------------------------
C
      DOUBLE PRECISION FUNCTION ALPHAS_not_timed(Q)
c
c     Evaluation of strong coupling constant alpha_S
c     Author: R.K. Ellis
c
c     q -- scale at which alpha_s is to be evaluated
c
c-- common block alfas.inc
c     asmz -- value of alpha_s at the mass of the Z-boson
c     nloop -- the number of loops (1,2, or 3) at which beta
c
c     function is evaluated to determine running.
c     the values of the cmass and the bmass should be set
c     in common block qmass.
C-----------------------------------------------------------------------------

      IMPLICIT NONE
c
      include 'alfas.inc'
      DOUBLE PRECISION Q,T,AMZ0,AMB,AMC
      DOUBLE PRECISION AS_OUT
      INTEGER NLOOP0,NF3,NF4,NF5
      PARAMETER(NF5=5,NF4=4,NF3=3)
C
      REAL*8       CMASS,BMASS
      COMMON/QMASS/CMASS,BMASS
      DATA CMASS,BMASS/1.42D0,4.7D0/  ! HEAVY QUARK MASSES FOR THRESHOLDS
C
      REAL*8 ZMASS
      DATA ZMASS/91.188D0/
C
      SAVE AMZ0,NLOOP0,AMB,AMC
      DATA AMZ0,NLOOP0/0D0,0/
      IF (Q .LE. 0D0) THEN
         WRITE(6,*) 'q .le. 0 in alphas'
         WRITE(6,*) 'q= ',Q
         STOP
      ENDIF
      IF (asmz .LE. 0D0) THEN
         WRITE(6,*) 'asmz .le. 0 in alphas',asmz
c         WRITE(6,*) 'continue with asmz=0.1185'
         STOP
         asmz=0.1185D0
      ENDIF
      IF (CMASS .LE. 0.3D0) THEN
         WRITE(6,*) 'cmass .le. 0.3GeV in alphas',CMASS
c         WRITE(6,*) 'continue with cmass=1.5GeV?'
         STOP
         CMASS=1.42D0
      ENDIF
      IF (BMASS .LE. 0D0) THEN
         WRITE(6,*) 'bmass .le. 0 in alphas',BMASS
         WRITE(6,*) 'COMMON/QMASS/CMASS,BMASS'
         STOP
         BMASS=4.7D0
      ENDIF
c--- establish value of coupling at b- and c-mass and save
      IF ((asmz .NE. AMZ0) .OR. (NLOOP .NE. NLOOP0)) THEN
         AMZ0=asmz
         NLOOP0=NLOOP
         T=2D0*DLOG(BMASS/ZMASS)
         CALL NEWTON1(T,asmz,AMB,NLOOP,NF5)
         T=2D0*DLOG(CMASS/BMASS)
         CALL NEWTON1(T,AMB,AMC,NLOOP,NF4)
      ENDIF

c--- evaluate strong coupling at scale q
      IF (Q  .LT. BMASS) THEN
           IF (Q  .LT. CMASS) THEN
             T=2D0*DLOG(Q/CMASS)
             CALL NEWTON1(T,AMC,AS_OUT,NLOOP,NF3)
           ELSE
             T=2D0*DLOG(Q/BMASS)
             CALL NEWTON1(T,AMB,AS_OUT,NLOOP,NF4)
           ENDIF
      ELSE
      T=2D0*DLOG(Q/ZMASS)
      CALL NEWTON1(T,asmz,AS_OUT,NLOOP,NF5)
      ENDIF
      ALPHAS_not_timed=AS_OUT
      RETURN
      END

      SUBROUTINE NEWTON1(T,A_IN,A_OUT,NLOOP,NF)
C     Author: R.K. Ellis

c---  calculate a_out using nloop beta-function evolution
c---  with nf flavours, given starting value as-in
c---  given as_in and logarithmic separation between
c---  input scale and output scale t.
c---  Evolution is performed using Newton's method,
c---  with a precision given by tol.

      IMPLICIT NONE
      INTEGER NLOOP,NF
      REAL*8 T,A_IN,A_OUT,AS,TOL,F2,F3,F,FP,DELTA
      REAL*8 B0(3:5),C1(3:5),C2(3:5),DEL(3:5)
      PARAMETER(TOL=5.D-4)
C---     B0=(11.-2.*NF/3.)/4./PI
      DATA B0/0.716197243913527D0,0.66314559621623D0,0.61009394851893D0/
C---     C1=(102.D0-38.D0/3.D0*NF)/4.D0/PI/(11.D0-2.D0/3.D0*NF)
      DATA C1/.565884242104515D0,0.49019722472304D0,0.40134724779695D0/
C---     C2=(2857.D0/2.D0-5033*NF/18.D0+325*NF**2/54)
C---     /16.D0/PI**2/(11.D0-2.D0/3.D0*NF)
      DATA C2/0.453013579178645D0,0.30879037953664D0,0.14942733137107D0/
C---     DEL=SQRT(4*C2-C1**2)
      DATA DEL/1.22140465909230D0,0.99743079911360D0,0.66077962451190D0/
      F2(AS)=1D0/AS+C1(NF)*LOG((C1(NF)*AS)/(1D0+C1(NF)*AS))
      F3(AS)=1D0/AS+0.5D0*C1(NF)
     & *LOG((C2(NF)*AS**2)/(1D0+C1(NF)*AS+C2(NF)*AS**2))
     & -(C1(NF)**2-2D0*C2(NF))/DEL(NF)
     & *ATAN((2D0*C2(NF)*AS+C1(NF))/DEL(NF))

      A_OUT=A_IN/(1D0+A_IN*B0(NF)*T)
      IF (NLOOP .EQ. 1) RETURN
      A_OUT=A_IN/(1D0+B0(NF)*A_IN*T+C1(NF)*A_IN*LOG(1D0+A_IN*B0(NF)*T))
      IF (A_OUT .LT. 0D0) AS=0.3D0
 30   AS=A_OUT

      IF (NLOOP .EQ. 2) THEN
      F=B0(NF)*T+F2(A_IN)-F2(AS)
      FP=1D0/(AS**2*(1D0+C1(NF)*AS))
      ENDIF
      IF (NLOOP .EQ. 3) THEN
      F=B0(NF)*T+F3(A_IN)-F3(AS)
      FP=1D0/(AS**2*(1D0+C1(NF)*AS+C2(NF)*AS**2))
      ENDIF
      A_OUT=AS-F/FP
      DELTA=ABS(F/FP/AS)
      IF (DELTA .GT. TOL) GO TO 30
      RETURN
      END

C-----------------------------------------------------------------------------
C
      double precision function mfrun(mf,scale,asmz,nloop)
C
C-----------------------------------------------------------------------------
C
C	This function returns the 2-loop value of a MSbar fermion mass
C       at a given scale.
C
C	INPUT: mf    = MSbar mass of fermion at MSbar fermion mass scale
C	       scale = scale at which the running mass is evaluated
C	       asmz  = AS(MZ) : this is passed to alphas(scale,asmz,nloop)
C              nloop = # of loops in the evolution
C
C
C
C	EXTERNAL:      double precision alphas(scale,asmz,nloop)
C
C-----------------------------------------------------------------------------
C
      implicit none
C
C     ARGUMENTS
C
      double precision  mf,scale,asmz
      integer           nloop
C
C     LOCAL
C
      double precision  beta0, beta1,gamma0,gamma1
      double precision  A1,as,asmf,l2
      integer  nf
C
C     EXTERNAL
C
      double precision  alphas
      external          alphas
c
c     CONSTANTS
c
      double precision  One, Two, Three, Pi
      parameter( One = 1.0d0, Two = 2.0d0, Three = 3.0d0 )
      parameter( Pi = 3.14159265358979323846d0)
      double precision tmass
      parameter(tmass=174d0)
cc
C
C
      if ( mf.gt.tmass ) then
         nf = 6
      else
         nf = 5
      end if

      beta0 = ( 11.0d0 - Two/Three *nf )/4d0
      beta1 = ( 102d0  - 38d0/Three*nf )/16d0
      gamma0= 1d0
      gamma1= ( 202d0/3d0  - 20d0/9d0*nf )/16d0
      A1    = -beta1*gamma0/beta0**2+gamma1/beta0
      as    = alphas(scale)
      asmf  = alphas(mf)
      l2    = (1+ A1*as/Pi)/(1+ A1*asmf/Pi)

      mfrun = mf * (as/asmf)**(gamma0/beta0)

      if(nloop.eq.2) mfrun =mfrun*l2
ccc
      return
      end

max_stars_repo_name: JaySandesara/madminer

max_issues_repo_name: JaySandesara/madminer

max_forks_repo_name: JaySandesara/madminer


! { dg-do compile }
! { dg-options "-ffree-form -std=f95" }
!
! Test errors for DEC-style PARAMETER statements with a real standard.
!

subroutine sub()
  implicit real(8) (A-Z)
  parameter pi = 3.1415926535d0 ! { dg-error "Legacy Extension: PARAMETER" }
  print *, pi
end subroutine

end


text:
C
C $Id: isgeti.f,v 1.5 2008-07-27 00:17:16 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE ISGETI (IPN,IVL)
C
      CHARACTER*(*) IPN
C
C This subroutine is called to retrieve the integer value of a specified
C parameter.
C
C IPN is the name of the parameter whose value is to be retrieved.
C
C IVL is an integer variable in which the desired value is to be
C returned by ISGETI.
C
C Get the real value of the parameter and return the integer equivalent.
C
      CALL ISGETR (IPN,RVL)
      IVL=INT(RVL)
C
C Done.
C
      RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl