text:
module type_functional_cache_handler
#include <messenger.h>
    use mod_kinds,                  only: rk, ik
    use type_chidg_data,            only: chidg_data_t
    use type_chidg_worker,          only: chidg_worker_t
    use type_cache_handler,         only: cache_handler_t
    use type_element_info,          only: element_info_t
    use mod_constants,              only: DIAG, ZERO, dQ_DIFF, dX_DIFF, dBC_DIFF, NO_DIFF
    use type_functional_cache,      only: functional_cache_t
    use type_evaluator,             only: evaluator_t
    use type_integral_cache,        only: integral_cache_t
    use type_geometry_cache,        only: geometry_cache_t
    use mod_DNAD_tools,             only: face_compute_seed, element_compute_seed

    use mod_chidg_mpi,              only: IRANK, NRANK, ChiDG_COMM, GLOBAL_MASTER
    use mod_io,                     only: verbosity
    use mpi_f08,                    only: MPI_Barrier, MPI_Bcast, MPI_AllReduce, MPI_Gather,   &
                                          MPI_INTEGER4, MPI_REAL8, MPI_LOGICAL, MPI_CHARACTER, &
                                          MPI_SUM, MPI_LOR
    use DNAD_D
    implicit none

    !>  An object to handle functional computation operations.
    !!
    !!  This allows to temporarily store and handle the data coming from the functional calculation.
    !!  It is meant to avoid the time consuming multiple loops over the elements/faces belonging
    !!  to the functional geometries.
    !!
    !!  It is also useful for easy processor communication in parallel mode.
    !!
    !!  The functional real and derivatives values are store in a sistematic way, independently whether
    !!  it is a surface or a volume type of functional.
    !!
    !!
    !!  @author Matteo Ugolotti
    !!  @date   2/23/2018
    !!
    !!  Restructured
    !!
    !!  @author Matteo Ugolotti
    !!  @date   3/7/2019
    !!
    !---------------------------------------------------------------------------------------------------
    type, public :: functional_cache_handler_t

        type(functional_cache_t),   pointer      :: cache            ! Point to the functional cache to be used for storing info
        type(chidg_data_t),         pointer      :: data             ! ChiDG data
        character(:),               allocatable  :: geom             ! Geometry we are computing integrals on (reference/auxiliary)
        character(:),               allocatable  :: component        ! Component element/face_interior for worker cache update
        integer(ik)                              :: dtype            ! Type of differentiation requested
        integer(ik)                              :: ifunc            ! Functional currently computed

    contains

        procedure,  public   :: update       ! Functional computation update

        ! Private procedures (made public for tests)
        procedure,  public   :: init         ! Initialize functional cache handler
        procedure,  public   :: compute      ! Compute functional
        procedure,  public   :: store        ! Share functional integral with other processors
        procedure,  public   :: destructor   ! Nullify all pointers

    end type functional_cache_handler_t
    !***************************************************************************************************

contains

    !>  This procedure begins the functional computation either on primary or auxiliary geometries
    !!
    !!  @author Matteo Ugolotti
    !!  @date   2/23/2018
    !!
    !!  restructured
    !!
    !!  @author matteo ugolotti
    !!  @date   3/7/2019
    !!
    !!  @param[inout]      worker           chidg worker
    !!  @param[inout]      data             chidg data
    !!  @param[inout]      fcl_cache        functional storage
    !!  @param[in]         ifunc            functional ID
    !!  @param[in]         geom_feature     geometry on which the functional needs to
    !!                                      be compute: reference/auxiliary
    !!  @param[in]         differentiate    flag for functional derivative computation
    !!
    !---------------------------------------------------------------------------------------------------
    subroutine update(self,worker,fcl_cache,data,ifunc,geom_feature,differentiate)
        class(functional_cache_handler_t),  intent(inout)       :: self
        type(chidg_worker_t),               intent(inout)       :: worker
        type(functional_cache_t),           intent(inout)       :: fcl_cache
        type(chidg_data_t),                 intent(inout)       :: data
        integer(ik),                        intent(in)          :: ifunc
        character(*),                       intent(in)          :: geom_feature
        integer(ik),                        intent(in)          :: differentiate

        ! Initialize functional_cache_handler
        call self%init(data,fcl_cache,ifunc,geom_feature,differentiate)

        ! Compute functional
        call self%compute(worker)

        ! Store functional
        call self%store()

        ! Destructor
        call self%destructor()

    end subroutine update
    !***************************************************************************************************

    !>  Here the handler is initialize so that it knows functional and geometry we are currently working on.
    !!  Pointers are created to facilitate the access to the main objects in the handler, that are:
    !!      - ChiDG data (for storage)
    !!      - Functional cache (for integral storage during functional computation)
    !!      - Functional (containing the procedures to actually compute the specific functional)
    !!
    !!  Furthermore, type of differentiation, geometry of current computation (auxiliary/reference) and
    !!  cache component are saved to be easily accessible during the computation.
    !!
    !!  @author Matteo Ugolotti
    !!  @date   3/7/2019
    !!
    !!  @param[inout]      worker           chidg worker
    !!  @param[inout]      data             chidg data
    !!  @param[inout]      fcl_cache        functional storage
    !!  @param[in]         ifunc            functional ID
    !!  @param[in]         geom_feature     geometry on which the functional needs to
    !!                                      be compute: reference/auxiliary
    !!  @param[in]         differentiate    flag for functional derivative computation
    !!
    !---------------------------------------------------------------------------------------------------
    subroutine init(self,data,fcl_cache,ifunc,geom_feature,differentiate)
        class(functional_cache_handler_t),  intent(inout)           :: self
        type(chidg_data_t),                 intent(inout),  target  :: data
        type(functional_cache_t),           intent(inout),  target  :: fcl_cache
        integer(ik),                        intent(in)              :: ifunc
        character(*),                       intent(in)              :: geom_feature
        integer(ik),                        intent(in)              :: differentiate

        ! Currently working on auxiliary geometry
        self%geom  =  trim(geom_feature)
        self%cache => fcl_cache
        self%data  => data

        associate (func => data%functional_group%fcl_entities(ifunc)%func )

            ! Assign active geometry (reference/auxiliary),data, cache, initialize active_cache
            select case (trim(geom_feature))
                case("auxiliary")
                    ! Initialize fuctional auxiliary cache
                    call self%cache%init(self%geom,data%mesh,func%auxiliary_geom,func%intermediate_integrals,func%get_int_type(),differentiate)

                case("reference")
                    ! Initialize fuctional auxiliary cache
                    call self%cache%init(self%geom,data%mesh,func%reference_geom,func%intermediate_integrals,func%get_int_type(),differentiate)

                case default
                    call chidg_signal_two(FATAL,"functional_cache_handler_t%init: wrong geometry string provided to the functional cache handler", trim(geom_feature),".")
            end select

            ! Define type of integral
            ! NOTE: this can be modified to handle different type of integrals on auxiliary and reference geometries
            select case (func%get_int_type())
                case('VOLUME INTEGRAL')
                    self%component = 'element'
                case('FACE INTEGRAL')
                    self%component = 'interior faces'
                case default
                    call chidg_signal_two(FATAL,"functional_cache_handler_t%init: wrong type of integral provided to the worker",func%get_int_type(),".")
            end select

            ! Save locally type of differentiation
            self%dtype = differentiate

            ! Functiona ID
            self%ifunc = ifunc

        end associate

    end subroutine init
    !***************************************************************************************************

    !>  TASKS:
    !!      1) Compute the functional at each element/face by looping through the entities in the cache and
    !!         store them in the functional cache
    !!      2) Carry out the communication among the processor of the integrals values (functional can be made
    !!         off multiple integrals)
    !!      3) Finalize the functional on either the auxiliary or reference geometry (operations done using overall
    !!         integrals.
    !!
    !!  @author Matteo Ugolotti
    !!  @date   2/25/2018
    !!
    !!  Added dX linearization capability
    !!
    !!  @author Matteo Ugolotti
    !!  @date   9/5/2018
    !!
    !!  restructured
    !!
    !!  @author matteo ugolotti
    !!  @date   3/7/2019
    !!
    !!  @param[inout]      worker           chidg worker
    !!
    !---------------------------------------------------------------------------------------------------
    subroutine compute(self,worker)
        class(functional_cache_handler_t),  intent(inout)       :: self
        type(chidg_worker_t),               intent(inout)       :: worker

        type(element_info_t)        :: elem_info
        type(cache_handler_t)       :: worker_cache_handler
        integer(ik)                 :: idepend, ientity, idiff, iface, ierr, idomain_l, ielement_l
        logical                     :: compute_auxiliary, compute_reference

        associate (func => self%data%functional_group%fcl_entities(self%ifunc)%func )

            ! Select the functional computation
            compute_auxiliary = ( self%geom == 'auxiliary' )
            compute_reference = ( self%geom == 'reference' )

            ! Loop through all the entities initialized (elements/faces)
            do ientity = 1,self%cache%nentities(self%geom)

                ! Update the worker
                if (compute_auxiliary) then
                    idomain_l  = self%cache%aux_cache%idomain_g%at(ientity)
                    ielement_l = self%cache%aux_cache%ielement_l%at(ientity)
                    iface      = self%cache%aux_cache%iface%at(ientity)
                else
                    idomain_l  = self%cache%ref_cache%idomain_l%at(ientity)
                    ielement_l = self%cache%ref_cache%ielement_l%at(ientity)
                    iface      = self%cache%ref_cache%iface%at(ientity)
                end if

                elem_info = worker%mesh%get_element_info(idomain_l,ielement_l)
                call worker%set_element(elem_info)

                ! Compute differential interpolator for mesh sensititivites
                call worker%mesh%compute_derivatives_dx(elem_info,self%dtype)

                ! Update worker cache
                call worker_cache_handler%update(worker,self%data%eqnset,       &
                                                 self%data%bc_state_group,      &
                                                 components = self%component,   &
                                                 face = iface,                  &
                                                 differentiate = self%dtype,    &
                                                 lift = .false.)

                ! Set up worker
                idepend = 1
                idiff   = DIAG

                if (self%component == 'interior faces') then
                    worker%function_info%seed = face_compute_seed(worker%mesh,          &
                                                                  elem_info%idomain_l,  &
                                                                  elem_info%ielement_l, &
                                                                  iface,idepend,idiff,worker%itime)
                    call worker%set_face(iface)

                else if (self%component == 'element') then
                    worker%function_info%seed = element_compute_seed(worker%mesh,          &
                                                                     elem_info%idomain_l,  &
                                                                     elem_info%ielement_l, &
                                                                     idepend,idiff,worker%itime)
                end if

                worker%function_info%idepend = idepend
                worker%function_info%idiff   = idiff
                worker%function_info%dtype   = self%dtype

                ! Compute functional on the entity (element/face)
                if (compute_auxiliary) call func%compute_auxiliary(worker,self%cache)
                if (compute_reference) call func%compute_functional(worker,self%cache)

                ! Release differential interpolators allocated memory
                call worker%mesh%release_derivatives_dx(elem_info,self%dtype)

            end do !ientity

            ! Communicate integrals through processors
            call self%cache%comm(self%geom)

            ! Finalize functional on auxiliary/reference geometry
            if (compute_auxiliary) call func%finalize_auxiliary(worker,self%cache)
            if (compute_reference) call func%finalize_functional(worker,self%cache)

        end associate

    end subroutine compute
    !***************************************************************************************************

    !>  This procedure stores the computed functional based on the type of the functional itself:
    !!      - functional derivatives (derivatives of the functional wrt dQ or dX are stored)
    !!      - functional             (only the real value of the functional is of interest and stored)
    !!      - auxiliary              (store the auxiliary value needed for the upcoming actuall functional)
    !!
    !!  - Functional derivatives are stored in the adjoint%Jq vector
    !!  - Functional real values are stored in the sdata%functional%func
    !!  - Auxiliary value is stored in the functional itself
    !!
    !!
    !!  @author Matteo Ugolotti
    !!  @date   2/25/2018
    !!
    !!  Added dX linearization capability
    !!
    !!  @author Matteo Ugolotti
    !!  @date   9/5/2018
    !!
    !!  restructured
    !!
    !!  @author matteo ugolotti
    !!  @date   3/7/2019
    !!
    !!
    !!  TODO: maybe this can be improved
    !!
    !---------------------------------------------------------------------------------------------------
    subroutine store(self)
        class(functional_cache_handler_t),  intent(inout)       :: self

        logical                 :: proceed_storing, store_functional,          &
                                   store_derivatives_dQ, store_derivatives_dX,  &
                                   store_auxiliary, compatibility
        type(integral_cache_t)  :: global_integral
        type(AD_D)              :: finalized_integral

        real(rk)                :: time
        Integer(ik)             :: idomain_l, ielement_l, ientity, i_int, step

        ! Prepare the type of store we need to execute
        if ( (self%geom == 'reference') .and. (self%dtype == dQ_DIFF) ) then
            store_derivatives_dQ = .true.
            store_derivatives_dX = .false.
            store_functional     = .true.
        end if
        if ( (self%geom == 'reference') .and. (self%dtype == dX_DIFF) ) then
            store_derivatives_dQ = .false.
            store_derivatives_dX = .true.
            store_functional     = .false.
        end if
        if ( (self%geom == 'reference') .and. (self%dtype == NO_DIFF) ) then
            store_derivatives_dQ = .false.
            store_derivatives_dX = .false.
            store_functional     = .true.
        end if
        if ( self%geom == 'auxiliary' ) then
            store_derivatives_dQ = .false.
            store_derivatives_dX = .false.
            store_functional     = .false.
        end if

        associate (func => self%data%functional_group%fcl_entities(self%ifunc)%func )

            ! Update local integrals with overall value coming from all ranks and store
            ! dQ derivatives in the sdata%adjoint%Jq
            if (store_derivatives_dQ) then
                self%data%sdata%adjoint%Jq(self%ifunc) = func%store_deriv(self%cache)
                call self%data%sdata%adjoint%Jq(self%ifunc)%assemble()
            end if

            ! Update local integrals with overall value coming from all ranks and store
            ! dX derivatives in the Jx(ifunc) vector
            if (store_derivatives_dX) then
                self%data%sdata%adjointx%Jx(self%ifunc) = func%store_deriv(self%cache)
                call self%data%sdata%adjointx%Jx(self%ifunc)%assemble()
            end if

            ! Store real value of the functional
            if (store_functional) then

                ! Store value
                call self%data%sdata%functional%func(self%ifunc)%push_back(func%store_value(self%cache))

                ! Store time and step
                ! Istep and itime are the same for all the functional
                ! We do not want to keep pushing times and steps, we do it once for the first
                ! functional.
                step = self%data%time_manager%istep
                if (self%ifunc == 1) call self%data%sdata%functional%step%push_back(step)
                time = self%data%time_manager%t
                if (self%ifunc == 1) call self%data%sdata%functional%time%push_back(time)

            end if

        end associate

    end subroutine store
    !***************************************************************************************************

    !>  Nullify all the pointers before moving to the next geometry or functional
    !!
    !!  @author matteo ugolotti
    !!  @date   3/8/2019
    !!
    !---------------------------------------------------------------------------------------------------
    subroutine destructor(self)
        class(functional_cache_handler_t),  intent(inout)       :: self

        ! Deassociate pointers in functional cache handler
        if (associated(self%cache)) nullify(self%cache)
        if (associated(self%data))  nullify(self%data)

    end subroutine destructor
    !***************************************************************************************************

end module type_functional_cache_handler

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text:
c###mpath.for
      SUBROUTINE MPATH(IF)
C-------------------------------
C
C     THIS ROUTINE PERFORMS THE MULTIPATH PROBABILITY AND DETERMINES
C     IF ANY OTHER MODE WILL INTERFERE WITH THE MOST RELIABLE MODE
c**********************************************************************
c          As changed 26 March 2009
c          This return the reliability of the next most probable mode
c          that is within the Maximum tolerable time delay (DMP)
c          and whose signal power is greater than PMP below the
c          signal power of the MRM.
c**********************************************************************
C
      COMMON /DON /ALATD, AMIN, AMIND, BTR, BTRD, DLONG, DMP, ERTR, GCD,
     1 GCDKM, PMP, PWR, TLAT, TLATD, TLONG, TLONGD, RSN, SIGTR, RLAT,
     2 RLATD,RLONG,RLONGD,BRTD,FLUX,ULAT,ULATD,ULONG,ULONGD,SSN,D90R,
     3 D50R,D10R,D90S,D50S,D10S
      COMMON / SON / ANGLE(13), ANGLER(13), CPROB(13), DBLOS(13),
     A DBLOSL(13), DBLOSU(13), DBU(13), DELAY(13), DBW(13), NHP(13),
     B XNYNOIS(13), PROBMP(13), RELIAB(13), SNDB(13), SNPR(13),
     C SNRLW(13), SNRUP(13), SPROB(13), VHIGH(13), RNEFF(13),MDL(13)
       CHARACTER MDL*1
      COMMON / allMODE /ABPS(20),CREL(20),FLDST(20),HN(20),HP(20),
     1PROB(20),RELY(20),RGAIN(20),SIGPOW(20),SN(20),
     2SPRO(20),TGAIN(20),TIMED(20),TLOSS(20),B(20),FSLOS(20),
     3grlos(20),adv(20),obf(20),
     CNMODE(20),TLLOW(20),TLHGH(20),EFF(20),NREL,NMMOD
      COMMON / FILES / LUI, LUO, LU2, LU5, LU6, LU15, LU16,
     1                 LU20, LU25, LU26, LU35
      COMMON / ION / IHRE, IHRO, IHRS, LUFP, METHOD, NOISE, NPSL
      COMMON / CONTRL / IELECT(3), KTOUT(12), MONTHS(12), SUNSP(12),
     A IANTOU, ICARD, INTEG, IRED, ISOUT, ISPROC, ISSN, ITYPE, JDASH,
     B JFREQ, JLONG, KCARD, KRUN, MAPIN, MAXNAM, MONOLD, MOREM, MORES,
     C NUMNAM, NUPROC, MAXMET, MSPEC, M100

      PROBMP(IF) = 0.001
      IF(DMP.le.0. .or. PMP.le.0.) return        !  no multipath requested
c          Is this METHOD=30? If so, not valid>7000km
      if(method.eq.20 .and. mspec.eq.121 .and. gcdkm.gt.7000.) return
      sig_power=sigpow(nrel)         !  signal power of Most-Reliable-Mode
      sig_power_limit=sig_power-PMP  !  ignore signals below this level
      TTIM = timed(nrel)             !  time delay of MRM
ccc      write(luo,4) method,mspec,gcdkm
ccc4     format('method, mspec=',2i5,f10.2)
ccc      write(luo,1) nrel,sig_power,sig_power_limit,ttim,pmp,dmp
ccc1     format('mpath=',i3,5f10.3)

      DO 135 IM = 1,NMMOD   ! DO LOOP up to max. number of modes.
      IF(IM.EQ.NREL) GO TO 135  ! Don't check mode for MRM -- found already by RELBIL.FOR.
      IF(HP(IM).LE.0.) GO TO 135  ! This mode does not exist.
ccc      write(luo,2) im,timed(im),sigpow(im)
ccc2     format(i4,2f10.3)
      IF(ABS(TIMED(IM)-TTIM).LE.DMP) GO TO 135   !  inside acceptable time delay

C        Check on Multipath Power criterion set by operator.
C        The logic of this next part is revised from original code.
c        Ignore signal powers below this limit
      if(sigpow(im).lt.sig_power_limit) go to 135

      PROBMP(IF)=AMAX1(PROBMP(IF),RELY(IM))
ccc      write(luo,3) probmp(if),rely(im)
ccc3     format(20x,2f10.3)
135   CONTINUE
      RETURN
      END
C--------------------------------

ext: for

max_stars_repo_name: jawatson/voacapl

max_issues_repo_name: jawatson/voacapl

max_forks_repo_name: jawatson/voacapl


SUBROUTINE deemat(dee,e,v)
!
! This subroutine returns the elastic dee matrix for ih=3 (plane strain),
! ih=4 (axisymmetry or plane strain elastoplasticity) or ih=6
! (three dimensions).
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::e,v
 REAL(iwp),INTENT(OUT)::dee(:,:)
 REAL(iwp)::v1,v2,c,vv,zero=0.0_iwp,pt5=0.5_iwp,one=1.0_iwp,two=2.0_iwp
 INTEGER::i,ih
 dee=zero
 ih=UBOUND(dee,1)
 v1=one-v
 c=e/((one+v)*(one-two*v))
 SELECT CASE(ih)
 CASE(3)
   dee(1,1)=v1*c
   dee(2,2)=v1*c
   dee(1,2)=v*c
   dee(2,1)=v*c
   dee(3,3)=pt5*c*(one-two*v)
 CASE(4)
   dee(1,1)=v1*c
   dee(2,2)=v1*c
   dee(4,4)=v1*c
   dee(3,3)=pt5*c*(one-two*v)
   dee(1,2)=v*c
   dee(2,1)=v*c
   dee(1,4)=v*c
   dee(4,1)=v*c
   dee(2,4)=v*c
   dee(4,2)=v*c
 CASE(6)
   v2=v/(one-v)
   vv=(one-two*v)/(one-v)*pt5
   DO i=1,3
     dee(i,i)=one
   END DO
   DO i=4,6
     dee(i,i)=vv
   END DO
   dee(1,2)=v2
   dee(2,1)=v2
   dee(1,3)=v2
   dee(3,1)=v2
   dee(2,3)=v2
   dee(3,2)=v2
   dee=dee*e/(two*(one+v)*vv)
 CASE DEFAULT
   WRITE(*,*)'wrong size for dee matrix'
 END SELECT
RETURN
END SUBROUTINE deemat


PROGRAM main

  USE mo_kind, ONLY: i4
  use mo_ansi_colors, only: color, c_red, c_green
  USE mo_nml,  ONLY: open_nml, close_nml, position_nml

  IMPLICIT NONE

  INTEGER, PARAMETER :: nnml = 100
  CHARACTER(len=*), PARAMETER :: nfile = 'namelist.txt'
  INTEGER(i4) :: stat

  LOGICAL :: isgood

  Write(*,*) ''
  Write(*,*) 'Test mo_nml.f90'

  isgood = .true.
  call open_nml(nfile, nnml)
  call position_nml('name1', unit=nnml, first=.true., status=stat)
  if (stat /= 0) isgood =.false.
  call position_nml('name2', unit=nnml, status=stat)
  !call position_nml('kkkk', unit=nnml)
  if (stat /= 0) isgood =.false.
  call position_nml('name3', unit=nnml, first=.true., status=stat)
  call close_nml()

  Write(*,*) ''
  if (isgood) then
     write(*,*) 'mo_nml ', color('o.k.', c_green)
  else
     write(*,*) 'mo_nml ', color('failed!', c_red)
  endif

END PROGRAM main


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

Hi RoseAnne! Is your Gluten/Dairy Free page about Gluten & Dairy Free food? Or is it just a business? It looks like youre creating a page for Natural Food Works which we already have. You can list places to find and buy Gluten/Dairy Free foods, but it shouldnt be a page for one place. But great idea for the additions to the Food page, you making those pages really adds to the wiki! Thank you! Users/MichelleAccurso

Also, look into our style of markup at these pages:
wiki:wikispot:Help with Editing
wiki:wikispot:Help with Formatting
Wiki Sandbox
wiki:wikispot:Help with Macros

Welcome to the wiki. Sorry to tramp through the first pages you created, Ive tried to preserve the information you contributed while moving it around to match the existing layout a little better. Users/JasonAller

20080210 13:08:34 nbsp Hi RoseAnne, how goes NFW? Users/StevenDaubert

20081214 18:35:19 nbsp Please dont completely blank a page. Perhaps reading Welcome to the Wiki/Business Owner will put it in perspective. Users/StevenDaubert

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE init_modular_coils (nvariables, xvariables, nfp)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE modular_coils
      IMPLICIT NONE
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: nc, i, n, nfp
      INTEGER :: nvariables, modes
      REAL(rprec) :: xvariables(*)
!-----------------------------------------------

      nc = nmod_coils_per_period

      nfper = nfp
      nmod_coils = nc*nfper
      nmid = (nc+1)/2
      nodd = MOD(nc,2)
      IF ((nodd.eq.0) .and. (.not.lsymm)) nmid = nmid + 1

      nvariables = 0
      nmod_coeffs = 0

      IF (nc .le. 0) RETURN

!     Initialize the variables to values of unique coil parameters
!     and count the number of variables

      n = 0

!     First consider the CASE with coils on both symmetry planes at
!     phi = 0 and phi = pi/nfp (lsymm = F). This implies that the number
!     of coils per field period must be even (nodd = 0).

      IF ((nodd .eq. 0) .and. (.not.lsymm)) THEN

!     Symmetry coil at phi = 0.

        i = 1
        DO modes = 1, nf_phi
          n = n + 1
          xvariables(n) = phis(i,modes)
        END DO

        DO modes = 1, nf_rho
          n = n + 1
          xvariables(n) = rhos(i,modes)
        END DO

!     Coils 2 through nmid-1.
        DO i = 2, nmid-1
          modes = 0
          n = n + 1
          xvariables(n) = phic(i,modes)
          DO modes = 1,nf_phi
            n = n + 1
            xvariables(n) = phic(i,modes)
            n = n + 1
            xvariables(n) = phis(i,modes)
          END DO

          DO modes = 1,nf_rho
            n = n + 1
            xvariables(n) = rhos(i,modes)
          END DO
        END DO

!     Symmetry coil at phi = pi/nfp.
        i = nmid
        DO modes = 1, nf_phi
          n = n + 1
          xvariables(n) = phis(i,modes)
        END DO

        DO modes = 1, nf_rho
          n = n + 1
          xvariables(n) = rhos(i,modes)
        END DO

!     Next consider the cases with a coil on phi = 0 (lsymm = T), or
!     on phi = pi/nfp (lsymm = F), but not both. Then there may be an
!     even number of coils per period (nodd = 0) or an odd number of
!     coils per period (nodd = 1).

      ELSE

        DO i = 1, nmid-nodd
          modes = 0
          n = n + 1
          xvariables(n) = phic(i,modes)
          DO modes = 1,nf_phi
            n = n + 1
            xvariables(n) = phic(i,modes)
            n = n + 1
            xvariables(n) = phis(i,modes)
          END DO

          DO modes = 1,nf_rho
            n = n + 1
            xvariables(n) = rhos(i,modes)
          END DO
        END DO

        IF (nodd .eq. 1) THEN
          i = nmid
          DO modes = 1, nf_phi
            n = n + 1
            xvariables(n) = phis(i,modes)
          END DO

          DO modes = 1, nf_rho
            n = n + 1
            xvariables(n) = rhos(i,modes)
          END DO
        END IF

      END IF  ! END IF ((nodd .eq. 0) .and. (lsymm .eqv. .false.))

      nmod_coeffs = n
      nvariables = n

      END SUBROUTINE init_modular_coils


text:
PROGRAM MULTIMODEFLOQUET

  USE ATOMIC_PROPERTIES
  USE TYPES
  USE SUBINTERFACE
  USE SUBINTERFACE_LAPACK
  USE FLOQUETINITINTERFACE
  USE ARRAYS

  IMPLICIT NONE
  TYPE(MODE),       DIMENSION(:),   ALLOCATABLE :: FIELDS
  TYPE(ATOM)                                       ID
  INTEGER,          DIMENSION(:),   ALLOCATABLE :: MODES_NUM
  INTEGER                                          TOTAL_FREQUENCIES,D_BARE
  INTEGER                                          INFO,m,INDEX0,r
  DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE :: ENERGY,E_FLOQUET
  COMPLEX*16,       DIMENSION(:,:), ALLOCATABLE :: H__,U_F,U_AUX,U_B2D
  DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: P_AVG
  DOUBLE PRECISION                              :: T1,T2

  INTEGER N_

  OPEN(UNIT=3,FILE="qubit_oscillation_DRIVER.dat",ACTION="WRITE")

  INFO = 0
  CALL FLOQUETINIT(ID,'qubit',INFO)
!  WRITE(*,*) total_states_lsi

  D_BARE = ID%D_BARE
  ALLOCATE(P_AVG(D_BARE,D_BARE))
  ALLOCATE(U_AUX(D_BARE,D_BARE))

  ALLOCATE(MODES_NUM(2))

  MODES_NUM(1) = 1 !(STATIC FIELD)
  MODES_NUM(2) = 1 !(DRIVING BY ONE HARMONIC)

  TOTAL_FREQUENCIES = SUM(MODES_NUM,1)
  ALLOCATE(FIELDS(TOTAL_FREQUENCIES))
  DO m=1,TOTAL_FREQUENCIES
     ALLOCATE(FIELDS(m)%V(ID%D_BARE,ID%D_BARE))
  END DO

  FIELDS(1)%X    = 0.0
  FIELDS(1)%Y    = 0.0
  FIELDS(1)%Z    = 1.0
  FIELDS(1)%phi_x = 0.0
  FIELDS(1)%phi_y = 0.0
  FIELDS(1)%phi_z = 0.0
  FIELDS(1)%omega = 0.0
  FIELDS(1)%N_Floquet = 0

  FIELDS(2)%X     = 2.0
  FIELDS(2)%Y     = 0.0
  FIELDS(2)%Z     = 0.0
  FIELDS(2)%phi_x = 0.0
  FIELDS(2)%phi_y = 0.0
  FIELDS(2)%phi_z = 0.0
  FIELDS(2)%omega = 1.0
  FIELDS(2)%N_Floquet = 20

  N_ = 32!128
  DO m=1,N_

     ! --- SET DRIVING PARAMETERS
     FIELDS(2)%omega = 0.2 + (m-1)*2.0/N_

     !--- EVALUATE TIME-EVOLUTION OPERATOR IN THE BARE BASIS
     T1 = 0.0
     DO r=1,N_
        T2 = r*32.0*4.0*atan(1.0)/N_
        CALL TIMEEVOLUTIONOPERATOR(ID,D_BARE,SIZE(MODES_NUM,1),SIZE(MODES_NUM,1),MODES_NUM,FIELDS,T1,T2,U_AUX,INFO)
        P_AVG = ABS(U_AUX)**2
        WRITE(3,*) t2,FIELDS(2)%OMEGA,ABS(U_AUX)**2
     END DO
     WRITE(3,*)
  END DO

END PROGRAM MULTIMODEFLOQUET

max_stars_repo_name: openMMF/Multimode-Floquet

max_issues_repo_name: openMMF/Multimode-Floquet

max_forks_repo_name: openMMF/Multimode-Floquet


module construct_vion_vh_floc_ffte_module

  use rgrid_variables, only: Ngrid, dV, Igrid
  use ffte_sub_module, only: zwork1_ffte, zwork2_ffte &
                           , comm_fftx, comm_ffty, comm_fftz, npuy, npuz
  use atom_module, only: ki_atom, aa_atom
  use bb_module, only: bb
  use ggrid_module, only: construct_Ggrid, MGL
  use ps_local_variables, only: vqlg
  use parallel_module, only: RSDFT_MPI_COMPLEX16, comm_grid, myrank &
                            ,construct_id_ir_parallel
  use rsdft_mpi_module, only: rsdft_allreduce_sum, rsdft_allgatherv
  use watch_module, only: watchb

  implicit none

  private
  public :: construct_vion_vh_floc_ffte

  logical :: first_time=.true.
  integer,allocatable :: LLG(:,:), LGHS(:,:)
  integer,allocatable :: LGHT(:,:), IGHT(:)
  real(8),allocatable :: GGHT(:)
  integer :: a1b,b1b,a2b,b2b,a3b,b3b,ab1,ab12
  complex(8),allocatable :: fg(:), GGVL(:)
  integer,allocatable :: icnta(:), idisa(:)
  integer :: MI_0, MI_1

contains

  subroutine construct_vion_vh_floc_ffte( rho, Vion, Vh, force_local, Eh )
    implicit none
    real(8),intent(in)  :: rho(:,:)
    real(8),intent(out) :: Vion(:), Vh(:), force_local(:,:)
    real(8),intent(out) :: Eh
    real(8),allocatable :: work(:,:)
    real(8) :: a1,a2,a3,pi2,Gx,Gy,Gz,Gr
    real(8) :: ttmp(2),ttt(2,13)
    complex(8),parameter :: z0=(0.0d0,0.0d0), zi=(0.0d0,1.0d0)
    complex(8) :: zsum1,zsum2,zsum3,ztmpv,ztmpf
    integer :: i1,i2,i3,is,i,ii,j,ik,a,ierr
    integer :: Nspin, Natom
    include 'mpif.h'

    !call watchb( ttmp, barrier='on' ); ttt=0.0d0

    !call write_border( 1, " construct_vion_vh_floc_ffte(start)" )

    Nspin = size( rho, 2 )
    Natom = size( force_local, 2 )
    pi2   = 2.0d0*acos(-1.0d0)

    if ( first_time ) then
      call initial_preparation
      call construct_id_ir_parallel( idisa,icnta,Natom,comm_fftx,MI_0,MI_1 )
      idisa=idisa*3
      icnta=icnta*3
      first_time = .false.
    end if

    !call watchb( ttmp, ttt(:,1), barrier='on' )

! ---

!$OMP parallel workshare
    Vh(:)=0.0d0
    fg(:)=z0
    GGVL(:)=z0
    force_local=0.0d0
    zwork1_ffte=z0
!$OMP end parallel workshare

    !call watchb( ttmp, ttt(:,2), barrier='on' )

!$OMP parallel do collapse(3) private(i)
    do i3=a3b,b3b
    do i2=a2b,b2b
    do i1=a1b,b1b
      i=1+(i1-a1b)+(i2-a2b)*ab1+(i3-a3b)*ab12
      zwork1_ffte(i1,i2,i3)=rho(i,1)
    end do
    end do
    end do
!$OMP end parallel do
    do is=2,Nspin
!$OMP parallel do collapse(3) private(i)
      do i3=a3b,b3b
      do i2=a2b,b2b
      do i1=a1b,b1b
        i=1+(i1-a1b)+(i2-a2b)*ab1+(i3-a3b)*ab12
        zwork1_ffte(i1,i2,i3)=zwork1_ffte(i1,i2,i3)+rho(i,is)
      end do
      end do
      end do
!$OMP end parallel do
    end do

    !call watchb( ttmp, ttt(:,3), barrier='on' )

    call MPI_Allreduce(zwork1_ffte,zwork2_ffte,size(zwork2_ffte) &
         ,RSDFT_MPI_COMPLEX16,MPI_SUM,comm_fftx,ierr)

    !call watchb( ttmp, ttt(:,4), barrier='on' )

    call pzfft3dv(zwork2_ffte,zwork1_ffte,Ngrid(1),Ngrid(2),Ngrid(3) &
         ,comm_ffty,comm_fftz,npuy,npuz,-1)

    !call watchb( ttmp, ttt(:,5), barrier='on' )

!$OMP parallel do private( i1,i2,i3,ii )
    do i=1,size(LGHT,2)
      i1=LGHT(1,i)
      i2=LGHT(2,i)
      i3=LGHT(3,i)
      fg(i)=conjg( zwork1_ffte(i1,i2,i3) )
    end do
!$OMP end parallel do

    !call watchb( ttmp, ttt(:,6), barrier='on' )

! ---

    do a=MI_0,MI_1

      ik=ki_atom(a)
      a1=pi2*aa_atom(1,a)
      a2=pi2*aa_atom(2,a)
      a3=pi2*aa_atom(3,a)

      zsum1=z0
      zsum2=z0
      zsum3=z0

!$OMP parallel do reduction(+:zsum1,zsum2,zsum3) private( i,i1,i2,i3,Gx,Gy,Gz,Gr,j,ztmpv,ztmpf )
    do i=1,size(LGHS,2)

      i1=LGHS(1,i)
      i2=LGHS(2,i)
      i3=LGHS(3,i)

      Gx=bb(1,1)*i1+bb(1,2)*i2+bb(1,3)*i3
      Gy=bb(2,1)*i1+bb(2,2)*i2+bb(2,3)*i3
      Gz=bb(3,1)*i1+bb(3,2)*i2+bb(3,3)*i3
      Gr=a1*i1+a2*i2+a3*i3

      j=MGL( IGHT(i) )

      ztmpv = vqlg(j,ik)*dcmplx( cos(Gr),-sin(Gr) )

      GGVL(i) = GGVL(i) + ztmpv

      ztmpf =-ztmpv*fg(i)*zi
      zsum1 = zsum1 + Gx*ztmpf
      zsum2 = zsum2 + Gy*ztmpf
      zsum3 = zsum3 + Gz*ztmpf

    end do ! i
!$OMP end parallel do

      force_local(1,a) = -zsum1*dV
      force_local(2,a) = -zsum2*dV
      force_local(3,a) = -zsum3*dV

    end do ! a

    !call watchb( ttmp, ttt(:,7), barrier='on' )

    call rsdft_allreduce_sum( force_local, comm_grid )

    !call watchb( ttmp, ttt(:,8), barrier='on' )

    call rsdft_allreduce_sum( GGVL, comm_fftx )

    !call watchb( ttmp, ttt(:,9), barrier='on' )

! ---

    Eh=0.0d0
!$OMP parallel private( i,ii,i1,i2,i3,ztmpv,ztmpf )
!$OMP workshare
    zwork2_ffte(:,:,:)=z0
!$OMP end workshare
!$OMP do reduction(+:Eh)
    do i=1,size(LGHT,2)
      i1=LGHT(1,i)
      i2=LGHT(2,i)
      i3=LGHT(3,i)
      ztmpv=zwork1_ffte(i1,i2,i3)
      ztmpf=ztmpv*GGHT(i)
      Eh=Eh+conjg(ztmpv)*ztmpf
      zwork2_ffte(i1,i2,i3)=ztmpf+GGVL(i)*Ngrid(0)
    end do
!$OMP end do
!$OMP end parallel

    !call watchb( ttmp, ttt(:,10), barrier='on' )

    call rsdft_allreduce_sum( Eh, comm_ffty )
    call rsdft_allreduce_sum( Eh, comm_fftz )

    !call watchb( ttmp, ttt(:,11), barrier='on' )

    Eh=-Eh*0.5d0*dV/Ngrid(0)

    call pzfft3dv(zwork2_ffte,zwork1_ffte,Ngrid(1),Ngrid(2),Ngrid(3) &
         ,comm_ffty,comm_fftz,npuy,npuz,1)

    !call watchb( ttmp, ttt(:,12), barrier='on' )

!$OMP parallel do collapse(3) private(i)
    do i3=a3b,b3b
    do i2=a2b,b2b
    do i1=a1b,b1b
      i = 1 + i1-a1b + (i2-a2b)*ab1 + (i3-a3b)*ab12
      Vion(i)=real( zwork1_ffte(i1,i2,i3) )
    end do
    end do
    end do
!$OMP end parallel do

    !call watchb( ttmp, ttt(:,13), barrier='on' )

    !if ( myrank == 0 ) then
    !   do i=1,size(ttt,2)
    !      write(*,'(1x,"vvf2(",i2,")",2f10.5)') i,ttt(:,i)
    !   end do
    !end if

    !call write_border( 1, " construct_vion_vh_floc_ffte(end)" )

  end subroutine construct_vion_vh_floc_ffte

  subroutine initial_preparation
    implicit none
    integer :: n,i,i1,i2,i3,j1,j2,j3
    real(8) :: pi4, g2

    pi4 = 4.0d0*acos(-1.0d0)

    a1b = Igrid(1,1)
    b1b = Igrid(2,1)
    a2b = Igrid(1,2)
    b2b = Igrid(2,2)
    a3b = Igrid(1,3)
    b3b = Igrid(2,3)
    ab1 = b1b-a1b+1
    ab12= ab1*(b2b-a2b+1)

    call construct_Ggrid( 0, LLG )

    n=0
    do i=1,size(LLG,2)
      i1=mod( Ngrid(1)+LLG(1,i), Ngrid(1) )
      i2=mod( Ngrid(2)+LLG(2,i), Ngrid(2) )
      i3=mod( Ngrid(3)+LLG(3,i), Ngrid(3) )
      if ( a2b <= i2 .and. i2 <= b2b .and. a3b <= i3 .and. i3 <= b3b ) n=n+1
    end do

    allocate( LGHS(3,n) ); LGHS=0
    allocate( LGHT(3,n) ); LGHT=0
    allocate( IGHT(n)   ); IGHT=0
    allocate( GGHT(n)   ); GGHT=0.0d0

    n=0
    do i=1,size(LLG,2)
      i1=LLG(1,i)
      i2=LLG(2,i)
      i3=LLG(3,i)
      j1=mod( Ngrid(1)+i1, Ngrid(1) )
      j2=mod( Ngrid(2)+i2, Ngrid(2) )
      j3=mod( Ngrid(3)+i3, Ngrid(3) )
      if ( a2b<=j2.and.j2<=b2b .and. a3b<=j3.and.j3<=b3b ) then
        n=n+1
        LGHS(1,n)=i1
        LGHS(2,n)=i2
        LGHS(3,n)=i3
        LGHT(1,n)=j1
        LGHT(2,n)=j2
        LGHT(3,n)=j3
        IGHT(n)=i
        if ( i1/=0 .or. i2/=0 .or. i3/=0 ) then
          g2=( bb(1,1)*i1+bb(1,2)*i2+bb(1,3)*i3 )**2 &
            +( bb(2,1)*i1+bb(2,2)*i2+bb(2,3)*i3 )**2 &
            +( bb(3,1)*i1+bb(3,2)*i2+bb(3,3)*i3 )**2
          GGHT(n)=pi4/g2
        end if
      end if
    end do

    deallocate( LLG )
    allocate( fg(n)   ); fg  =(0.0d0,0.0d0)
    allocate( GGVL(n) ); GGVL=(0.0d0,0.0d0)

  end subroutine initial_preparation

end module construct_vion_vh_floc_ffte_module


text:
!***************************************************************************************************************************
! This is the code from the section (DYNAMIC == RATE) lines 4075 - 4215 of the original CSCAS code. The names of the
! dummy arguments are the same as in the original CSCAS code and the call statement and are declared here. The variables
! that are not arguments are declared in module YCA_First_Trans_m. Unless identified as by MF, all comments are those of
! the original CSCAS.FOR code.
!
! Subroutine YCA_Growth_Evapo calculates the water and thermal conditions and water factor for initiation of growth of the
! planting material (called germination).
!***************************************************************************************************************************

    SUBROUTINE YCA_Growth_Evapo ( &
        ALBEDOS     , BRSTAGE     , CLOUDS      , CO2         , DLAYR       , DUL         , EO          , EOP         , &
        ES          , ISWWAT      , KEP         , LL          , NLAYR       , RLV         , RWUMX       , RWUPM       , &
        SAT         , SRAD        , SW          , TAIRHR      , TDEW        , TMAX        , TMIN        , TRWUP       , &
        UH2O        , &
        !WEATHER     ,                                                                                ! MF WEATHER needed for VPD
        WINDSP      , YEAR        , ST          &    !LPM20MAR2016 To consider ST for germination
        )

        USE ModuleDefs
        USE YCA_First_Trans_m
        USE YCA_Control_Environment

        IMPLICIT NONE

        !TYPE (WeatherType) WEATHER                                                            ! MF Defined in ModuleDefs

        INTEGER NLAYR       , YEAR
        INTEGER CSIDLAYR                                                                      ! Integer function call.

        REAL    ALBEDOS     , BRSTAGE     , CLOUDS      , CO2         , DLAYR(NL)   , DUL(NL)     , EO          , EOP
        REAL    ES          , KEP         , LL(NL)      , RLV(NL)     , RWUMX       , RWUPM       , SAT(NL)     , SRAD
        REAL    SW(NL)      , TAIRHR(24)  , TDEW        , TMAX        , TMIN        , TRWUP       , UH2O(NL)    , WINDSP
        REAL    ST(NL)                                  !LPM20MAR2016 To consider ST for germination
        REAL    CSVPSAT     , TFAC4                     ! Real function call.  !LPM 19SEP2017 Added tfac5

        CHARACTER(LEN=1) IDETG       , ISWWAT

          IF (PLYEAR <= 0) PLYEAR = YEAR

!-----------------------------------------------------------------------
!         Calculate potential plant evaporation,water uptake if neeeded
!-----------------------------------------------------------------------

          ! EO is brought into the module. The following calculations
          ! (apart from the root water uptake module) are for
          ! comparative purposes only. The root water uptake module
          ! is not necessary when running in CSM, but is necessary for
          ! CROPSIM.

          ! Co2 effect on stomatal resistances. General for C3 crops
          RLF  = 9.72 + 0.0757 * 330.0 + 10.0
          RLFC = 9.72 + 0.0757 *  CO2  + 10.0
          YEARDOY = YEARDOY                                                 ! MF For WORK.OUT

          IF (FILEIOT == 'DS4'.AND.IDETG /= 'N'.OR.FILEIOT /= 'DS4') THEN

            IF (ISWWAT /= 'N') THEN

                ! Call 1  Basic calculations with rcrop = 0
                CALL EVAPO('M',SRAD,CLOUDS,TMAX,TMIN,TDEW,WINDSP,ALBEDO,RATM,RCROP*0.0, &
                    EO,EOPEN,EOMPEN,EOPT,EOEBUD,TCAN,'M')

                EOPENC = EOPENC + EOPEN
                EOPTC = EOPTC + EOPT
                EOMPENC = EOMPENC + EOMPEN
                EOEBUDC = EOEBUDC + EOEBUD

                ! CSM with LAI=1.0,CROP=WH,TAVG=20.0,WINDSP=86.4 has
                ! RATM = 55  RCROP = 45
                ! Monteith had RATM = 300, RCROP = 150->500

                ! Call 2  Using rcrop as read-in from spe file
                CALL EVAPO('M',SRAD,CLOUDS,TMAX,TMIN,TDEW,WINDSP,ALBEDO,RATM,RCROP, &
                    EO,EOPEN,EOMPCRP,EOPT,EOEBUDCRP,TCAN,'M')
                EOMPCRPC = EOMPCRPC + EOMPCRP
                EOEBUDCRPC = EOEBUDCRPC + EOEBUDCRP

                ! Call 3 Using rcrop adjusted for CO2 effect
                IF (RLF > 0.0)CALL EVAPO('M',SRAD,CLOUDS,TMAX,TMIN,TDEW,WINDSP,ALBEDO,RATM,RCROP*RLFC/RLF, &
                    EO,EOPEN,EOMPCRPCO2,EOPT,EOEBUDCRPCO2,TCAN,'M')
                EOMPCRPCO2C = EOMPCRPCO2C + EOMPCRPCO2
                EOEBUDCRPCO2C = EOEBUDCRPCO2C + EOEBUDCRPCO2

                ! Transpiration ratio (Pot.pl.evap/Pot.soil evap)
                EPSRATIO = 1.0
                IF (EOMPEN > 0.0) EPSRATIO = EOMPCRPCO2 / EOMPEN
                TRATIO = 1.0
                IF (EOMPCRP > 0.0) TRATIO = EOMPCRPCO2 / EOMPCRP
            ENDIF

            IF (fileiot(1:2) /= 'DS') THEN
              ! Calculate plant potential evaporation
              EOP = MAX(0.0,EO/EOMPEN*EOMPCRPCO2 * (1.0-EXP(-LAI*KEP)))
              ! Ratio necessary because EO method may not be Monteith
              CALL CSCRPROOTWU(ISWWAT,NLAYR, DLAYR, LL, SAT, WFEU, MEWNU,EOP, RLV, RWUPM, RLFWU, RWUMX, RTDEP, &
                  SW, WTDEP, uh2o, trwup, trwu)
            ENDIF

            ! Call 4 Using rcrop adjusted for CO2 & H2O effect
            ! NB. Using previous days WFP.
            ! If use this for other than comparison, must check  LAH
            IF (RLF>0.0 .AND. WFP<=1.0 .AND. WFP>0.0) THEN
             CALL EVAPO('M',SRAD,CLOUDS,TMAX,TMIN,TDEW,WINDSP,ALBEDO,RATM,RCROP*RLFC/RLF*(1.0+(1.0-WFP)), &
                 EO,EOPEN,EOMPCRPCO2H2O,EOPT,EOEBUDCRPCO2H2O,TCAN,'M')
            ENDIF
            EOMPCRPCO2H2OC = EOMPCRPCO2H2OC + EOMPCRPCO2H2O
            EOEBUDCRPCO2H2OC = EOEBUDCRPCO2H2OC + EOEBUDCRPCO2H2O

            ! Call 5 to calcuate canopy temperature
            TVR1 = (TRWU*10.0+ES)
            IF (RLF>0.0) THEN
                CALL EVAPO('M',SRAD,CLOUDS,TMAX,TMIN,TDEW,WINDSP,ALBEDO,RATM,RCROP*RLFC/RLF*(1.0+(1.0-WFP)), &
                EO,tvr1,tvr2,tvr3,tvr4,TCAN,'C')
            ENDIF
          ENDIF

          ! Cumulative potential ET as used
          IF (EO > 0.0) EOC = EOC + EO

          ! Cumulative canopy-air temperature difference
          TDIFSUM = TDIFSUM+(TCAN-TMEAN)
          TDIFNUM = TDIFNUM + 1
          TDIFAV = TDIFSUM/TDIFNUM

!!-----------------------------------------------------------------------
!!         Calculate thermal time                       !LPM 21MAR2016 Moved after estimation of water factors to consider ST(LSEED)
!!-----------------------------------------------------------------------
!
!          Tfd = TFAC4(trdv1,tmean,TT)
!          IF (brstage+1.0 < 10.0) &
!           Tfdnext = TFAC4(trdv2,tmean,TTNEXT)
!          IF (trgem(3) > 0.0) THEN
!            Tfgem = TFAC4(trgem,ST(LSEED),TTGEM)
!            !Tfgem = TFAC4(trgem,tmean,TTGEM)  !LPM 20MAR2016 To modify
!          ELSE
!            Ttgem = tt
!          ENDIF
!          !IF (Cfllflife == 'D') THEN
!          !  ! Leaf life read-in as days (eg.7 phyllochrons->7 days)
!          !  Ttlflife = Phints
!          !ELSE
!          !  !Tflflife = TFAC4(trdv1,tmean,TTlflife)
!            Tflflife = TFAC4(trdv3,tmean,TTlflife)                         ! LPM 18MAR15 modified trdv1 to trdv3 to consider the cardinal temperatures for leaf development
!            Tflfsize = TFAC4(trdv4,tmean,TTlfsize)                         ! LPM 18MAR15 modified trdv1 to trdv4 to consider different optimum temperature for leaf size
!          !ENDIF
!
!-----------------------------------------------------------------------
!         Calculate soil water 'status' (Used as a sw 'potential')
!-----------------------------------------------------------------------

          DO L = 1,NLAYR
            SWP(L) = &
             AMIN1(1.0,AMAX1(.0,((SW(L)-LL(L))/(DUL(L)-LL(L)))))
          ENDDO

!-----------------------------------------------------------------------
!         Calculate water factor for germination
!-----------------------------------------------------------------------

          IF (LSEED < 0) LSEED = CSIDLAYR (NLAYR, DLAYR, SDEPTH)
          WFGE = 1.0
          IF (ISWWAT /= 'N') THEN
              !IF (GESTAGE < 1.0) THEN  !LPM 21MAR2015 to estimate WFGE for germination and emergence
              IF (EMRGFR < 1.0) THEN
              !IF (LSEED < 0) LSEED = CSIDLAYR (NLAYR, DLAYR, SDEPTH) !LPM 21MAR2015 To estimate LSEED when ISWWAT == N (linking ST and germination)
                  !IF (LSEED > 1) THEN
                  !  SWPSD = SWP(LSEED)
                  !ELSE
                  ! SWP(0) = AMIN1(1.,AMAX1(.0,(SWP(1)-0.5*(SWP(2)-SWP(1)))))
                  ! SWPSD = SWP(0) + (SDEPTH/DLAYR(1))*(SWP(2)-SWP(0))
                  !ENDIF
                  IF (WFGEM > 0.0) &
                   WFGE = AMAX1(0.0,AMIN1(1.0,(SWP(LSEED)/WFGEM)))
            ENDIF
          ENDIF
          IF (ISWWATCROP == 'N') WFGE = 1.0

!-----------------------------------------------------------------------
!         Calculate thermal time
!-----------------------------------------------------------------------

          Tfd = TFAC4(trdv1,tmean,TT)
          IF (brstage+1.0 < PSX) & !LPM 23JUL19 set as limit PSX instead of 10
           Tfdnext = TFAC4(trdv2,tmean,TTNEXT)
          IF (trgem(3) > 0.0) THEN
            Tfgem = TFAC4(trgem,ST(LSEED),TTGEM)
            !Tfgem = TFAC4(trgem,tmean,TTGEM)  !LPM 20MAR2016 To modify the TT for germination using ST
          ELSE
            Ttgem = tt
          ENDIF
          !IF (Cfllflife == 'D') THEN
          !  ! Leaf life read-in as days (eg.7 phyllochrons->7 days)
          !  Ttlflife = Phints
          !ELSE
          !  !Tflflife = TFAC4(trdv1,tmean,TTlflife)
          !LPM 09OCT2019 Remove TFLfgrowth because it is the same than TFG
          !Tflfgrowth = calculateTemperatureFactor(trdv3,tmean,TTlfgrowth)                         ! LPM 18MAR15 modified trdv1 to trdv3 to consider the cardinal temperatures for leaf development
            Tflfsize = TFAC4(trdv4,tmean,TTlfsize)                         ! LPM 18MAR15 modified trdv1 to trdv4 to consider different optimum temperature for leaf size
            Tflflife = calculateTemperatureFactor(trlfl,tmean,TTlflife)              !LPM 14SEP2017 Added new cardinal temperatures for leaf life and other for leaf growth (trdv3 or trlfg)
            !ENDIF

      END SUBROUTINE YCA_Growth_Evapo

max_stars_repo_name: vstocca/dssat-csm-os

max_issues_repo_name: vstocca/dssat-csm-os

max_forks_repo_name: vstocca/dssat-csm-os


text:
Address(Jerome Street) is a residential street in West Davis.

 Intersecting Streets

  Marina Circle and across the intersection Bermuda Avenue
  Marina Circle again
  Portsmouth Avenue
  Breton Avenue


text:
            program demo_abs
            integer :: i = -1, iout
            real :: x = -1.e0, xout, zout
            complex :: z = (-3.e0,4.e0)
            doubleprecision :: r8 = 45.78D+00, dout
               write(*,*)'INPUTS:',i,x,z,r8
               iout = abs(i)
               xout = abs(x)
               zout = abs(z)
               dout = abs(r8)
               write(*,*)'OUTPUTS:',iout,xout,zout,dout
               write ( *, '(a,f12.4,12x,f12.4)' ) ' Double precision  ', -r8, abs(r8)
               ! COMPLEX
               ! 3 - 4 -5 right triangle test :
               write(*,*)'The abs() of (3.0,4.0) should be 5.0',abs((3.0,4.0))
               ! ELEMENTAL
               write(*,*)'abs is ELEMENTAL: ',abs([-10, 20, 0, -1, -3, 100])
            end program demo_abs

max_stars_repo_name: urbanjost/fortran-intrinsic-manpages

max_issues_repo_name: urbanjost/fortran-intrinsic-manpages

max_forks_repo_name: urbanjost/fortran-intrinsic-manpages


text:
      SUBROUTINE NBINT(I,IKS,IR,FIRR,FD)
*
*
*       Irregular integration.
*       ----------------------
*
      INCLUDE 'common6.h'
      COMMON/CHAINC/  XC(3,NCMAX),UC(3,NCMAX),BODYC(NCMAX),ICH,
     &                LISTC(LMAX)
      COMMON/PREDICT/ TPRED(NMAX)
      REAL*8  XI(3),XIDOT(3),FIRR(3),FREG(3),FD(3),FDUM(3),DX(3),DV(3)
      REAL*8  FCM(3),FCMD(3),FP(6),FPD(6)
      SAVE TCALL,TCALL2,ISTART,TRY
      DATA ISTART /0/
      LOGICAL TRY
*
*
*       Predict body #I (and possible KS components) to order FDOT.
      CALL JPRED(I)
*
*       Initialize binary search time to current TTOT at start/restart.
      IF (ISTART.EQ.0) THEN
          ISTART = 1
          TCALL = TTOT
          TCALL2 = TTOT
          TRY = .FALSE.
      END IF
*
*       Check regularization criterion for single particles.
      IF (STEP(I).LT.DTMIN.AND.I.LE.N) THEN
*       See whether dominant body can be regularized.
          IF (IKS.EQ.0) THEN
              CALL SEARCH(I,IKS)
          END IF
      END IF
*
*       Include close encounter search for low-eccentric massive binaries.
      IF (IKS.EQ.0.AND.TCALL.LT.TTOT.AND.STEP(I).LT.8.0*DTMIN) THEN
*       Consider massive single bodies in absence of subsystems.
          IF (I.LE.N.AND.BODY(I).GT.2.0*BODYM.AND.NSUB.EQ.0) THEN
*
*       Advance TCALL and obtain two-body elements & relative perturbation.
              TCALL = TTOT + 0.01
              JMIN = 0
              CALL ORBIT(I,JMIN,SEMI,ECC,GI)
*
              EB = -0.5*BODY(I)*BODY(JMIN)/SEMI
              IF (EB.LT.EBH.AND.GI.LT.0.001.AND.JMIN.GE.IFIRST) THEN
                  APO = SEMI*(1.0 + ECC)
*       Check eccentricity (cf. max perturbation) and neighbour radius.
                  IF (ECC.LT.0.25.AND.APO.LT.0.02*RS(I)) THEN
*                     WRITE (6,3)  NAME(I), NAME(JMIN), ECC, SEMI, EB
*   3                 FORMAT (' KS TRY:    NAM E A EB ',
*    &                                     2I6,F7.3,1P,2E10.2)
                      IKS = IKS + 1
                      ICOMP = I
                      JCOMP = JMIN
*       Reduce TCALL below TTOT for second call to ensure equal TNEW.
                      TCALL = TTOT - 1.0D-04
                  END IF
              END IF
          END IF
      END IF
*
*       Perform KS search for massive wide binary candidates.
      IF (IKS.EQ.0.AND.TCALL2.LT.TTOT.AND.STEP(I).LT.30.0*DTMIN) THEN
*       Note extra safety tests inside routine SWEEP2.
          IF (BODY(I).GT.20.0*BODYM.OR.TRY) THEN
*       Use a logical variable to catch lighter second KS component.
              TRY = .FALSE.
              TCALL2 = TCALL2 + 0.01
              DTCL = 30.0*DTMIN
              RCL = 12.0*RMIN
              CALL SWEEP2(I,IKS,DTCL,RCL)
*       Reduce TCALL2 after first IKS > 0 to allow second call.
              IF (IKS.GT.0) THEN
                  TCALL2 = TTOT - 1.0D-04
                  TRY = .TRUE.
              END IF
          END IF
      END IF
*
*       Obtain irregular force & first derivative.
      DO 5 K = 1,3
          XI(K) = X(K,I)
          XIDOT(K) = XDOT(K,I)
    5 CONTINUE
*
*       Assume small mass at centre for special case of no neighbours.
      NNB0 = LIST(1,I)
      IF (NNB0.EQ.0) THEN
          RI2 = XI(1)**2 + XI(2)**2 + XI(3)**2
          FIJ = 0.01*BODYM/(RI2*SQRT(RI2))
          RDOT = 3.0*(XI(1)*XIDOT(1) + XI(2)*XIDOT(2) +
     &                                 XI(3)*XIDOT(3))/RI2
          DO 10 K = 1,3
              FIRR(K) = -FIJ*XI(K)
              FD(K) = -(XIDOT(K) - RDOT*XI(K))*FIJ
   10     CONTINUE
          IF (I.GT.N) IPAIR = I - N
          GO TO 70
      END IF
*
*       Choose force loop for single particle or regularized c.m. body.
      IF (I.LE.N) GO TO 20
*
*       Set KS pair index and components.
      IPAIR = I - N
      I2 = 2*IPAIR
      I1 = I2 - 1
*
*       Adopt c.m. approximation for small total perturbation.
      IF (LIST(1,I1).GT.0) THEN
*       Correct irregular force on perturbed c.m. body (including any chain).
          CALL CMFIRR(I,I1,FIRR,FD)
          GO TO 70
      END IF
*
*       Copy c.m. coordinates & velocities for rare unperturbed intruder.
      DO 15 K = 1,3
          X(K,I1) = XI(K)
          X(K,I2) = XI(K)
          XDOT(K,I1) = XIDOT(K)
          XDOT(K,I2) = XIDOT(K)
   15 CONTINUE
*
*       Set neighbour number & list index of the last single particle.
   20 NNB1 = NNB0 + 1
      NNB2 = NNB1
   25 IF (LIST(NNB2,I).LE.N) GO TO 30
      NNB2 = NNB2 - 1
      IF (NNB2.GT.1) GO TO 25
*       Include special case of only c.m. neighbours.
      GO TO 40
*
*       Sum over single particles (unperturbed case included).
*  30 DO 35 L = 2,NNB2
*         K = LIST(L,I)
*         A1 = X(1,K) - XI(1)
*         A2 = X(2,K) - XI(2)
*         A3 = X(3,K) - XI(3)
*         DV(1) = XDOT(1,K) - XIDOT(1)
*         DV(2) = XDOT(2,K) - XIDOT(2)
*         DV(3) = XDOT(3,K) - XIDOT(3)
*         RIJ2 = A1*A1 + A2*A2 + A3*A3
*
*         DR2I = 1.0/RIJ2
*         DR3I = BODY(K)*DR2I*SQRT(DR2I)
*         DRDV = 3.0*(A1*DV(1) + A2*DV(2) + A3*DV(3))*DR2I
*
*         FIRR(1) = FIRR(1) + A1*DR3I
*         FIRR(2) = FIRR(2) + A2*DR3I
*         FIRR(3) = FIRR(3) + A3*DR3I
*         FD(1) = FD(1) + (DV(1) - A1*DRDV)*DR3I
*         FD(2) = FD(2) + (DV(2) - A2*DRDV)*DR3I
*         FD(3) = FD(3) + (DV(3) - A3*DRDV)*DR3I
*  35 CONTINUE
*
*       Replace irregular force loop by fast routine in C++ and OpenMP.
*  30 CALL CNBINT(I,X,XDOT,BODY,NNB2,LIST(2,I),FIRR,FD)
*
*       See whether any c.m. neighbours should be considered.
   30 IF (NNB2.EQ.NNB1) GO TO 60
*
*       Perform differential correction for active KS neighbours.
   40 DO 50 LL = NNB2+1,NNB1
          J = LIST(LL,I)
*       See whether to sum over binary components.
          JPAIR = J - N
          J1 = 2*JPAIR - 1
*       Skip unperturbed binary (treated as single particle).
          IF (LIST(1,J1).EQ.0) THEN
              GO TO 50
          END IF
          IF (TIME - TPRED(J).EQ.0.0D0) THEN
              ZZ = 1.0
              IF (GAMMA(JPAIR).GT.1.0D-04) ZZ = 0.0
              CALL KSRES2(JPAIR,J1,J2,ZZ)
          ELSE
              CALL JPRED(J)
          END IF
          J2 = J1 + 1
*       Obtain individual c.m. force with single particle approximation.
          A1 = X(1,J) - X(1,I)
          A2 = X(2,J) - X(2,I)
          A3 = X(3,J) - X(3,I)
          RIJ2 = A1*A1 + A2*A2 + A3*A3
*
          DV(1) = XDOT(1,J) - XDOT(1,I)
          DV(2) = XDOT(2,J) - XDOT(2,I)
          DV(3) = XDOT(3,J) - XDOT(3,I)
          DR2I = 1.0/RIJ2
          DR3I = BODY(J)*DR2I*SQRT(DR2I)
          DRDV = 3.0*(A1*DV(1) + A2*DV(2) + A3*DV(3))*DR2I
*
          FCM(1) = A1*DR3I
          FCM(2) = A2*DR3I
          FCM(3) = A3*DR3I
          FCMD(1) = (DV(1) - A1*DRDV)*DR3I
          FCMD(2) = (DV(2) - A2*DRDV)*DR3I
          FCMD(3) = (DV(3) - A3*DRDV)*DR3I
*
*       Evaluate perturbation due to first component of body #J.
          dr2 = 0.0
          drdv = 0.0
          DO 42 L = 1,3
              dx(L) = X(L,J1) - X(L,I)
              dv(L) = XDOT(L,J1) - XDOT(L,I)
              dr2 = dr2 + dx(L)**2
              drdv = drdv + dx(L)*dv(L)
   42     CONTINUE
*
          dr2i = 1.0/dr2
          dr3i = BODY(J1)*dr2i*SQRT(dr2i)
          drdv = 3.0*drdv*dr2i
*
          DO 45 L = 1,3
              FP(L) = dx(L)*dr3i
              FPD(L) = (dv(L) - dx(L)*drdv)*dr3i
   45     CONTINUE
*
*       Evaluate perturbation due to second component.
          dr2 = 0.0
          drdv = 0.0
          DO 46 L = 1,3
              dx(L) = X(L,J2) - X(L,I)
              dv(L) = XDOT(L,J2) - XDOT(L,I)
              dr2 = dr2 + dx(L)**2
              drdv = drdv + dx(L)*dv(L)
   46     CONTINUE
*
          dr2i = 1.0/dr2
          dr3i = BODY(J2)*dr2i*SQRT(dr2i)
          drdv = 3.0*drdv*dr2i
*
          DO 48 L = 1,3
              FP(L+3) = dx(L)*dr3i
              FPD(L+3) = (dv(L) - dx(L)*drdv)*dr3i
   48     CONTINUE
*
*       Accumulate individual correction terms together for accuracy.
          DO 49 L = 1,3
              FIRR(L) = FIRR(L) + (FP(L) + FP(L+3) - FCM(L))
              FD(L) = FD(L) + (FPD(L) + FPD(L+3) - FCMD(L))
   49     CONTINUE
   50 CONTINUE
*
*       Include differential force treatment for regularized subsystem.
   60 IF (NCH.GT.0) THEN
*       Distinguish between chain c.m. and any other particle.
          IF (NAME(I).EQ.0) THEN
              CALL CHFIRR(I,0,XI,XIDOT,FIRR,FD)
          ELSE
*       See if chain perturber list contains body #I.
              NP1 = LISTC(1) + 1
              DO 65 L = 2,NP1
                  J = LISTC(L)
                  IF (J.GT.I) GO TO 70
                  IF (J.EQ.I) THEN
                      CALL FCHAIN(I,1,XI,XIDOT,FIRR,FD)
                      GO TO 70
                  END IF
   65         CONTINUE
          END IF
      END IF
*
*       Check option for external tidal field using predicted FREG.
   70 DT = TIME - T0(I)
      IF (KZ(14).GT.0) THEN
          DTR = TIME - T0R(I)
          DO 75 K = 1,3
              FREG(K) = FR(K,I) + FRDOT(K,I)*DTR
   75     CONTINUE
*       Jpetts- added argument I to end of XTRNLF (zero here)
          CALL XTRNLF(XI,XIDOT,FIRR,FREG,FD,FDUM,0,0)
      END IF
*
*       Include the corrector and set new T0, F, FDOT, D1, D2 & D3.
      DTSQ = DT**2
      DT6 = 6.0D0/(DT*DTSQ)
      DT2 = 2.0D0/DTSQ
      DTSQ12 = ONE12*DTSQ
      DT13 = ONE3*DT
      T0(I) = TIME
*
      DO 80 K = 1,3
          DF = FI(K,I) - FIRR(K)
          FID = FIDOT(K,I)
          SUM = FID + FD(K)
          AT3 = 2.0D0*DF + DT*SUM
          BT2 = -3.0D0*DF - DT*(SUM + FID)
*
          X0(K,I) = XI(K) + (0.6D0*AT3 + BT2)*DTSQ12
          X0DOT(K,I) = XIDOT(K) + (0.75D0*AT3 + BT2)*DT13
*
          FI(K,I) = FIRR(K)
          FIDOT(K,I) = FD(K)
*       Use total force for irregular step (cf. Makino & Aarseth PASJ, 1992).
          FDUM(K) = FIRR(K) + FR(K,I)
*
          D0(K,I) = FIRR(K)
          D1(K,I) = FD(K)
          D2(K,I) = (3.0D0*AT3 + BT2)*DT2
          D3(K,I) = AT3*DT6
*       NOTE: These are real derivatives!
   80 CONTINUE
*
*       Specify new time-step by standard criterion (STEPI version not good).
      TTMP = TSTEP(FDUM,FD,D2(1,I),D3(1,I),ETAI)
      DT0 = TTMP
*
*     IF (I.GT.N) THEN
*       Check for hierarchical configuration but exclude small perturbations.
*         IF (H(IPAIR).LT.-ECLOSE.AND.KZ(36).GT.0) THEN
*             IF (GAMMA(IPAIR).GT.1.0E-04) THEN
*                 CALL KEPLER(I,TTMP)
*                 DT0 = TTMP
*             END IF
*         END IF
*     END IF
*
*       Include convergence test for large step (cf. Makino, Ap.J. 369, 200).
      IF (TTMP.GT.STEPJ.AND.N.GT.1000) THEN
         DV2 = 0.0
         F2 = 0.0
         DO 85 K = 1,3
            DV2 = DV2 + (XDOT(K,I) - X0DOT(K,I))**2
            F2 = F2 + FIRR(K)**2
   85    CONTINUE
*       Employ Jun's criterion to avoid over-shooting (cf. Book, 2.16).
         DTJ = STEP(I)*(1.0D-06*STEP(I)**2*F2/DV2)**0.1
         TTMP = MIN(TTMP,DTJ)
      END IF
*
*       Select discrete value (increased by 2, decreased by 2 or unchanged).
      IF (TTMP.GT.2.0*STEP(I)) THEN
          IF (DMOD(TIME,2.0*STEP(I)).EQ.0.0D0) THEN
              TTMP = MIN(2.0*STEP(I),SMAX)
*       Include factor 4 increase for FPOLY initializations with small STEP.
              IF (DT0.GT.10.0*TTMP.AND.
     &            DMOD(TIME,2.0D0*TTMP).EQ.0.0D0) THEN
                  TTMP = MIN(2.0*TTMP,SMAX)
              END IF
          ELSE
              TTMP = STEP(I)
          END IF
      ELSE IF (TTMP.LT.STEP(I)) THEN
          TTMP = 0.5*STEP(I)
          IF (TTMP.GT.DT0) THEN
              TTMP = 0.5*TTMP
          END IF
      ELSE
          TTMP = STEP(I)
      END IF
*
*       Set new block step and update next time.
      STEP(I) = TTMP
*       Jpetts - enforce minimum timestep
  111 IF (STEP(I) < 1.0d-8) THEN
          STEP(I) = 2.0*STEP(I)
          GO TO 111
      END IF
      TNEW(I) = STEP(I) + T0(I)
*
*       See whether any KS candidates are in the same block as body #I.
      IF (IKS.GT.0.AND.I.EQ.ICOMP) THEN
*       Accept same time, otherwise reduce STEP(ICOMP) and/or delay.
          IF (T0(JCOMP).EQ.T0(ICOMP)) THEN
              ICOMP = MIN(ICOMP,JCOMP)
              JCOMP = MAX(I,JCOMP)
          ELSE IF (T0(JCOMP) + STEP(JCOMP).LT.T0(ICOMP)) THEN
              STEP(ICOMP) = 0.5D0*STEP(ICOMP)
              TNEW(ICOMP) = STEP(ICOMP) + T0(ICOMP)
              IKS = 0
          ELSE
              IKS = 0
          END IF
      END IF
*
*       See whether total force & derivative needs updating.
      IF (IR.EQ.0) THEN
*       Combine regular force & first derivatives to obtain F & FDOT.
          DTR = TIME - T0R(I)
          DO 90 K = 1,3
              F(K,I) = 0.5D0*(FRDOT(K,I)*DTR + FR(K,I) + FIRR(K))
              FDOT(K,I) = ONE6*(FRDOT(K,I) + FD(K))
   90     CONTINUE
      END IF
*
*       Increase step counter and count perturbed c.m. steps.
*     NSTEPI = NSTEPI + 1
      IF (I.GT.N) THEN
          IF (LIST(1,2*IPAIR-1).GT.0) NSTEPB = NSTEPB + 1
      END IF
*
      RETURN
*
      END


text:
c----------------------------------------------------------------------
c
        module rigidBodyForce
          integer, allocatable :: rbIndex(:)
          real*8, allocatable  :: rbForce(:,:)
          real*8, allocatable  :: rbForceOld(:,:)
          real*8, allocatable  :: rbTorque(:,:)
          real*8, allocatable  :: rbTorqueOld(:,:)
c
          real*8, allocatable  :: rbDisp(:,:)
          real*8, allocatable  :: rbTotalDisp(:,:)
          real*8, allocatable  :: rbVel(:,:)
          real*8, allocatable  :: rbVelOld(:,:)
          real*8, allocatable  :: rbAcc(:,:)
          real*8, allocatable  :: rbAccOld(:,:)
          real*8, allocatable  :: rbAng(:)
          real*8, allocatable  :: rbAngOld(:)
          real*8, allocatable  :: rbTotalAng(:)
        end module
c
c----------------------------------------------------------------------
c
c----------------------------------------------------------------------
c
        subroutine malloc_rbForce
c
        use rigidbody_m
        use rigidBodyReadData
        use rigidBodyForce
        use number_def_m
c
        allocate( rbForce(numrbs, 3)  )
        allocate( rbForceOld(numrbs, 3)  )
        allocate( rbTorque(numrbs, 3) )
        allocate( rbTorqueOld(numrbs, 3)  )
c
        allocate( rbDisp(numrbs, 3) )
        allocate( rbTotalDisp(numrbs, 3) )
        allocate( rbVel(numrbs, 3)  )
        allocate( rbVelOld(numrbs, 3)  )
        allocate( rbAcc(numrbs, 3)  )
        allocate( rbAccOld(numrbs, 3)  )
        allocate( rbAng(numrbs)  )
        allocate( rbAngOld(numrbs)  )
        allocate( rbTotalAng(numrbs)  )
c
        rbForce  = zero
        rbTorque = zero
        rbDisp = zero
        rbVel = zero
        rbAcc = zero
        rbAng = zero
c
        rbTorqueOld = zero
c
c.... rbParam = rbTotalDisp + rbVelOld + rbAccOld + rbForceOld
c             + rbTotalAng  + rbAngOld
c
        if (rbUseReadData .eq. 1) then
c
          if (rbParamSize .ne. 14) then
            if(myrank .eq. master)
     &        write(*,*) "change rigid body parameter size rbParamSize"
            call error('rbParamSize','change size',rbParamSize)
          endif
c
          rbTotalDisp(:,1:3) = rbParamRead(:,1:3)
          rbVelOld(:,1:3) = rbParamRead(:,4:6)
          rbAccOld(:,1:3) = rbParamRead(:,7:9)
          rbForceOld(:,1:3) = rbParamRead(:,10:12)
          rbTotalAng(:) = rbParamRead(:,13)
          rbAngOld(:) = rbParamRead(:,14)
        else
          rbTotalDisp = zero
          rbVelOld = zero
          rbAccOld = zero
          rbForceOld = zero
          rbTotalAng = zero
          rbAngOld = zero
        endif
c
        if (allocated(rbParamRead))
     &    deallocate( rbParamRead )
c
        return
        end
c
c----------------------------------------------------------------------
c
c----------------------------------------------------------------------
c
      subroutine set_rbBC (x,  iBC,  BC,  BC_flow)
c
        use rigidBodyReadData
        use rigidBodyForce
c
        include "common.h"
c
        real*8    x(numnp,nsd)
        dimension iBC(nshg), BC(nshg,3), BC_flow(nshg,3)
        real*8    t
        real*8, dimension(nsd) :: ra, rc, r_disp
c
c.... loop over mesh vertices
        do i = 1,numnp
          if ( (ibits(iBC(i),14,3) .eq. 7) .and.
     &         (rbFlags(i) .gt. 0) ) then
c.... calculate the displacement due to rotation
            t = rbAngOld(rbFlags(i))/180.0*pi
            ra = rb_prop(rbFlags(i),5:7)
            rc = rb_prop(rbFlags(i),8:10)
            r_disp(1) =
     &               (cos(t)+ra(1)*ra(1)*(1-cos(t))-1)*(x(i,1)-rc(1))
     &         + (ra(1)*ra(2)*(1-cos(t))-ra(3)*sin(t))*(x(i,2)-rc(2))
     &         + (ra(1)*ra(3)*(1-cos(t))+ra(3)*sin(t))*(x(i,3)-rc(3))
            r_disp(2) =
     &           (ra(2)*ra(1)*(1-cos(t))+ra(3)*sin(t))*(x(i,1)-rc(1))
     &             + (cos(t)+ra(2)*ra(2)*(1-cos(t))-1)*(x(i,2)-rc(2))
     &         + (ra(2)*ra(3)*(1-cos(t))-ra(1)*sin(t))*(x(i,3)-rc(3))
             r_disp(2) =
     &           (ra(3)*ra(1)*(1-cos(t))-ra(2)*sin(t))*(x(i,1)-rc(1))
     &         + (ra(3)*ra(2)*(1-cos(t))+ra(1)*sin(t))*(x(i,2)-rc(2))
     &             + (cos(t)+ra(3)*ra(3)*(1-cos(t))-1)*(x(i,3)-rc(3))
c.... update flow BC
            BC_flow(i,1:3) = rbVelOld(rbFlags(i), 1:3)
          endif
        enddo ! end loop numnp
c
      return
      end
c
c----------------------------------------------------------------------
c
c----------------------------------------------------------------------
c
        subroutine init_rbForce
c
        use rigidBodyForce
        use number_def_m
c
        rbForce  = zero
        rbTorque = zero
c
        return
        end
c
c----------------------------------------------------------------------
c
c----------------------------------------------------------------------
c
        subroutine commu_rbForce
c
        use rigidBodyReadData
        use rigidBodyForce
c
        include "common.h"
        include "mpif.h"
        include "auxmpi.h"
c
        real*8, dimension(3) :: Forin, Forout
c
        do i = 1, numrbs
c.... XXX need to add rotation force
c
          if (rbsMM(i) .ne. 1 .and. rbsMM(i) .ne. 2)
     &      call error('rigidBodyBCElas','not support mode',rbsMM(i))
          if (iter .eq. nitr) then
            Forin  = (/ rbForce(i,1), rbForce(i,2), rbForce(i,3) /)
            if (numpe > 1) then
              call MPI_ALLREDUCE ( Forin,  Forout, 3,
     &                             MPI_DOUBLE_PRECISION,
     &                             MPI_SUM, MPI_COMM_WORLD, ierr)
              rbForce(i,1:3) = Forout(1:3)
            endif
          endif
        enddo
c
        return
        end
c
c----------------------------------------------------------------------
c
c----------------------------------------------------------------------
c
        subroutine release_rbForce
c
        use rigidBodyForce
c
        include "common.h"
c
        if (allocated(rbForce))
     &    deallocate( rbForce )
        if (allocated(rbTorque))
     &    deallocate( rbTorque )
        if (allocated(rbDisp))
     &    deallocate( rbDisp )
        if (allocated(rbTotalDisp))
     &    deallocate( rbTotalDisp )
        if (allocated(rbForceOld))
     &    deallocate( rbForceOld )
        if (allocated(rbTorqueOld))
     &    deallocate( rbTorqueOld )
        if (allocated(rbVel))
     &    deallocate( rbVel )
        if (allocated(rbVelOld))
     &    deallocate( rbVelOld )
        if (allocated(rbAcc))
     &    deallocate( rbAcc )
        if (allocated(rbAccOld))
     &    deallocate( rbAccOld )
        if (allocated(rbAng))
     &    deallocate( rbAng )
        if (allocated(rbAngOld))
     &    deallocate( rbAngOld )
        if (allocated(rbTotalAng))
     &    deallocate( rbTotalAng )
c
        return
        end
c
c----------------------------------------------------------------------
c
c----------------------------------------------------------------------
c
        subroutine local_rbIndex ( ien )
c
        use rigidBodyReadData
        use rigidBodyForce
c
        include "common.h"
c
        dimension    ien(npro,nshl)
        dimension    lnode(27)
        integer, dimension(npro, nshl) :: rbFlagsl
        integer   i,  j,  k,  counter
c
        allocate ( rbIndex(npro) )
        rbIndex = 0 ! 0 means it isn't a rigid body
c
c.... localize rigid body flag
c
        do i = 1, nshl
          rbFlagsl(:,i) = rbFlags(ien(:,i))
        enddo
c
c.... get local node index
c
        call getbnodes(lnode)
c
c.... if all boundary nodes share the same flag, set element flag
c
        do i = 1, npro
          counter = 0
          do j = 1, nenbl
            if ( (rbFlagsl(i, lnode(j)) .gt. 0) .and.
     &           (rbFlagsl(i, lnode(j)) .eq. rbFlagsl(i, lnode(1))) ) then
              counter = counter + 1
            endif
          enddo
          if (counter .eq. nenbl) then ! we find a mesh face on rigid body
            rbIndex(i) = rbFlagsl(i, lnode(1))
          endif
        enddo
c
        return
        end
c
c----------------------------------------------------------------------
c
c----------------------------------------------------------------------
c
        subroutine release_rbIndex
c
        use rigidBodyForce
c
        if (allocated(rbIndex))
     &    deallocate( rbIndex )
c
        return
        end
c
c----------------------------------------------------------------------
c
c----------------------------------------------------------------------
c
        subroutine calc_rbMotion
c
        use rigidBodyReadData
        use rigidBodyForce
c
        include "common.h"
        include "mpif.h"
        include "auxmpi.h"
c
        real*8                           t_mag
        real*8, dimension(nsd)        :: t_dir,  tmpAcc
        real*8, dimension(numrbs,nsd) :: t_acc
        real*8, dimension(3)          :: Forin1, Forin2
        real*8, dimension(1)          :: Forin3
c.... generalized-alpha parameters
        real*8  am, af, gm, bt
c
        am = 0.5 * (3.0-rhoinf_rb(1)) / (1.0+rhoinf_rb(1))
        af = 1.0 / (1.0+rhoinf_rb(1))
        gm = 0.5 + am - af
        bt = 0.25 * (1.0+am-af) * (1.0+am-af)
c
c.... loop over rigid body
        do j = 1,numrbs
          if (rbsMM(j) .ne. 1 .and. rbsMM(j) .ne. 2)
     &      call error('rigidBodyBCElas','not support mode',rbsMM(j))
c.... debugging {
c.... XXX need to implement rotation motion calculation
c          call core_get_centroid(rbMTs(j), cent(j,:))
c.... debugging }
c
c.... apply constraint on force
c
          t_dir(:) = rb_prop(j,2:4) ! translation direction
          t_mag = sqrt( t_dir(1)*t_dir(1)
     &                + t_dir(2)*t_dir(2)
     &                + t_dir(3)*t_dir(3) )
          if (t_mag .lt. 1e-12) then ! no constraint
            write(*,*)
     &      "WARNING: translation direction magnitude is too small!"
            t_dir(1:3) = 1.0
          else
            t_dir(1:3) = t_dir(1:3) / t_mag
          endif
c
          rbForce(j,1:3) = rbForce(j,1:3) * t_dir(1:3)
c
c.... set initial acceleration
c
          if ( (istep .eq. 0) .and. (rbUseReadData .eq. 0) ) then
            rbAccOld(j,1:3) = rbForce(j,1:3) / rb_prop(j,1)
            rbForceOld(j,1:3) = rbForce(j,1:3)
          endif
c
c.... prepare variable
c
          tmpAcc(1:3) = af/am/rb_prop(j,1) * rbForce(j,1:3)
     &                + (1-af)/am/rb_prop(j,1) * rbForceOld(j,1:3)
c
c.... get translation displacement
c
          rbDisp(j,1:3) = Delt(1) * rbVelOld(j,1:3)
     &                  + Delt(1) * Delt(1) * (
     &                    (0.5 + bt/am) * rbAccOld(j,1:3)
     &                  - bt * tmpAcc(1:3)    )
c
c.... get translation velocity
c
          rbVel(j,1:3) = rbVelOld(j,1:3)
     &                 + Delt(1) * (1.0 - gm/am) * rbAccOld(j,1:3)
     &                 + Delt(1) * gm * tmpAcc(1:3)
c
c.... get translation acceleration
c
          rbAcc(j,1:3) = tmpAcc(1:3) - (1.0-am)/am * rbAccOld(j,1:3)
c
c
c.... get rotation angle for this time step
c
          if (rbsMM(j) .eq. 2) then ! translation-spinning mode
c.... the angle = spinning ratio * magitude of translation disp
            rbAng(j) = rb_prop(j,11) * sqrt(rbDisp(j,1)*rbDisp(j,1)
     &                                     +rbDisp(j,2)*rbDisp(j,2)
     &                                     +rbDisp(j,3)*rbDisp(j,3))
          endif
c
c.... broadcast rb velocity and displacement
c.... broadcast rotation angle
c
          if ((rb_commuMotion .eq. 1) .and. (numpe .gt. 1)) then
            Forin1  = (/ rbVel(j,1), rbVel(j,2), rbVel(j,3) /)
            call MPI_BCAST ( Forin1(1),  3,  MPI_DOUBLE_PRECISION,
     &                       master,     MPI_COMM_WORLD,  ierr)
            rbVel(j,1:3) = Forin1(1:3)
c
            Forin2  = (/ rbDisp(j,1),   rbDisp(j,2),   rbDisp(j,3) /)
            call MPI_BCAST ( Forin2(1),  3,  MPI_DOUBLE_PRECISION,
     &                       master,     MPI_COMM_WORLD,  ierr)
            rbDisp(j,1:3) = Forin2(1:3)

            Forin3  = (/ rbAng(j) /)
            call MPI_BCAST ( Forin3(1),  1,  MPI_DOUBLE_PRECISION,
     &                       master,     MPI_COMM_WORLD,  ierr)
            rbAng(j) = Forin3(1)
          endif
c
c.... debugging {
c          if (myrank .eq. master) then
c            write(*,*) "rbForce", rbForce(j,1)
c     &                ,"disp:", rbTotalDisp(j,1)
c     &                ,"vel:",rbVel(j,1)
c     &                ,"acc:",rbAcc(j,1)
c     &                ,"ang:",rbTotalAng(j)
c          endif
c.... debugging }
c
        enddo
c
        return
        end
c
c----------------------------------------------------------------------
c
c----------------------------------------------------------------------
c
        subroutine update_rbParam
c
        use rigidBodyForce
c
        include "common.h"
c
        do j = 1,numrbs
          rbTotalDisp(j,1:3) = rbTotalDisp(j,1:3) + rbDisp(j,1:3)
          rbAccOld(j,1:3) = rbAcc(j,1:3)
          rbVelOld(j,1:3) = rbVel(j,1:3)
          rbForceOld(j,1:3) = rbForce(j,1:3)
          rbTotalAng(j) = rbTotalAng(j) + rbAng(j)
          rbAngOld(j) = rbAng(j)
c
c.... debugging {
c          if (myrank .eq. master) then
c            write(*,*) "rbForce", rbForce(j,1)
c     &                ,"disp:", rbTotalDisp(j,1)
c     &                ,"vel:",rbVel(j,1)
c     &                ,"acc:",rbAcc(j,1)
c     &                ,"ang:",rbTotalAng(j)
c          endif
c.... debugging }
c
        enddo
c
        return
        end
c
c----------------------------------------------------------------------
c
c----------------------------------------------------------------------
c
        subroutine synchronize_rbParam
c
        use rigidBodyReadData
        use rigidBodyForce
        use core_rigid_body
c
        include "common.h"
c
        real*8, dimension(numrbs) :: ax, ay, az, px, py, pz, sc
c
c.... XXX need to add rotation motion here
        do j = 1,numrbs
          if (rbsMM(j) .ne. 1 .and. rbsMM(j) .ne. 2)
     &      call error('rigidBodyBCElas','not support mode',rbsMM(j))
        enddo
c
        ax = rb_prop(1:numrbs,5)
        ay = rb_prop(1:numrbs,6)
        az = rb_prop(1:numrbs,7)
        px = rb_prop(1:numrbs,8)
        py = rb_prop(1:numrbs,9)
        pz = rb_prop(1:numrbs,10)
        sc = 0.0
c
        call core_update_rbms(rbTotalDisp(:,1),
     &                        rbTotalDisp(:,2),
     &                        rbTotalDisp(:,3),
     &                        ax, ay, az,
     &                        px, py, pz,
     &                        rbTotalAng(:), sc,
     &                        rbMTs(:), numrbs)
c
        return
        end
c
c----------------------------------------------------------------------
c
c----------------------------------------------------------------------
c
        subroutine write_rbParam
c
c        use rigidBodyReadData
        use rigidBodyForce
c
        include "common.h"
c
c.... rbParam = rbTotalDisp + rbVelOld + rbAccOld + rbForceOld
c             + rbTotalAng  + rbAngOld
c
        real*8, dimension(numrbs, rbParamSize) :: rbParam
c
        if (rbParamSize .ne. 14) then
          if(myrank .eq. master)
     &      write(*,*) "change rigid body parameter size rbParamSize"
        endif
c
        rbParam(:,1:3)   = rbTotalDisp(:,1:3)
        rbParam(:,4:6)   = rbVelOld(:,1:3)
        rbParam(:,7:9)   = rbAccOld(:,1:3)
        rbParam(:,10:12) = rbForceOld(:,1:3)
        rbParam(:,13)    = rbTotalAng(:)
        rbParam(:,14)    = rbAngOld(:)
c
        call write_field(
     &       myrank,  'a'//char(0),'rbParams'//char(0), 8,
     &       rbParam, 'd'//char(0), numrbs, rbParamSize, lstep)
c
        return
        end
c
c----------------------------------------------------------------------
c

max_stars_repo_name: SCOREC/phasta

max_issues_repo_name: SCOREC/phasta

max_forks_repo_name: SCOREC/phasta


text:
! KGEN-generated Fortran source file
!
! Filename    : rrsw_ref.f90
! Generated at: 2015-07-27 00:47:03
! KGEN version: 0.4.13

    MODULE rrsw_ref
        USE kgen_utils_mod, ONLY : kgen_dp, check_t, kgen_init_check, kgen_print_check
        USE shr_kind_mod, ONLY: r8 => shr_kind_r8
        !      use parkind, only : jpim, jprb
        IMPLICIT NONE
        !------------------------------------------------------------------
        ! rrtmg_sw reference atmosphere
        ! Based on standard mid-latitude summer profile
        !
        ! Initial version:  JJMorcrette, ECMWF, jul1998
        ! Revised: MJIacono, AER, jun2006
        !------------------------------------------------------------------
        !  name     type     purpose
        ! -----  :  ----   : ----------------------------------------------
        ! pref   :  real   : Reference pressure levels
        ! preflog:  real   : Reference pressure levels, ln(pref)
        ! tref   :  real   : Reference temperature levels for MLS profile
        !------------------------------------------------------------------
        REAL(KIND=r8), dimension(59) :: preflog
        REAL(KIND=r8), dimension(59) :: tref
        PUBLIC kgen_read_externs_rrsw_ref
    CONTAINS

    ! write subroutines

    ! module extern variables

    SUBROUTINE kgen_read_externs_rrsw_ref(kgen_unit)
        INTEGER, INTENT(IN) :: kgen_unit
        READ(UNIT=kgen_unit) preflog
        READ(UNIT=kgen_unit) tref
    END SUBROUTINE kgen_read_externs_rrsw_ref

    END MODULE rrsw_ref

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang