text:
      PROGRAM LPTEST

c      Compute the coefficients in the Legendre polynomial expansions
c      of the phase matrix elements.  Compare values from angular
c      integration with values computed directly by MIEV0.

c      NOTES:

c       (0) Do NOT use the NoPMOM version of MIEV0 here.

c       (1) In MIEV0, set the PARAMETER MAXANG to 10001 (in order
c           to do the angular integrations) and set MAXTRM to 150
c           or so (ALL occurrences; to save storage).

c       (2) Large differences between the exact and quadrature
c           values can appear, but these are invariably in cases
c           where the Legendre coefficients are many orders of
c           magnitude down from their maximum values.  These
c           differences can be seen to be mainly due to quadrature
c           error, by comparing the 1001-point to the immediately
c           following 100001-point quadrature (which also shows
c           how hard it would be to get Legendre coefficients by
c           quadrature in some cases).

c       (3) Because of the size of the output, only a few cases
c           are sampled here.  The user is urged to try cases
c           corresponding to his own interests; but remember that
c           computer time escalates rapidly for size parameters
c           over 100 or so.

c      CALLING TREE:

c         LPTEST
c            ERRMSG
c            MIEV0
c            SIMPSO
c               ERRMSG

c ----------------------------------------------------------------------
c ------------  I/O SPECIFICATIONS FOR SUBROUTINE  MIEV0  --------------
c ----------------------------------------------------------------------
      PARAMETER  ( MAXANG = 10001, MOMDIM = 250 )
      LOGICAL  ANYANG, PERFCT, PRNT( 2 )
      INTEGER  IPOLZN, NUMANG, NMOM
      REAL     GQSC, MIMCUT, PMOM( 0:MOMDIM, 4 ), QEXT, QSCA, SPIKE,
     $         XMU( MAXANG ), XX
      COMPLEX  CREFIN, SFORW, SBACK, S1( MAXANG ), S2( MAXANG ),
     $         TFORW( 2 ), TBACK( 2 )
c ----------------------------------------------------------------------

c         *** Definitions of local variables ***

c   PROD      : Product of a Legendre polynomial and a phase quantity

c   PMOMAP    : Moments calculated by angular integration

c   PQ( I,N ) : I-th phase quantity at N-th angle

c   PL( N )   : Legendre polynomial P-sub-MOM  of argument  XMU(N)
c   PLm1( N ) : Legendre Polynomial P-sub-(MOM-1)  of argument  XMU(N)

c   STEP      : Size of integration interval in Simpson quadrature

c     .. Local Scalars ..

      CHARACTER PMINDX*6
      INTEGER   IANG, IPOL, ISZPAR, J, MOM, N, NREFIN
      REAL      FNORM, PMOMAP, POL, STEP
      COMPLEX   CTMP
c     ..
c     .. Local Arrays ..

      REAL      PL( MAXANG ), PLM1( MAXANG ),
     &          PQ( 4, MAXANG ), PROD( MAXANG )
c     ..
c     .. External Functions ..

      REAL      SIMPSO
      EXTERNAL  SIMPSO
c     ..
c     .. External Subroutines ..

      EXTERNAL  ERRMSG, MIEV0
c     ..
c     .. Intrinsic Functions ..

      INTRINSIC ABS, AIMAG, CONJG, REAL
c     ..
c     .. Statement Functions ..

      REAL      SQ
c     ..
c     .. Statement Function definitions ..

      SQ( CTMP ) = REAL( CTMP )**2 + AIMAG( CTMP )**2
c     ..

c                      ** Set some input variables to MIEV0
      PRNT(1) = .FALSE.
      PRNT(2) = .FALSE.
      PERFCT = .FALSE.
      ANYANG = .FALSE.
      MIMCUT = 1.E-6

      WRITE( *, 9000 )

c          ** Loop over values of IPOLZN, CREFIN, XX, and NUMANG

      DO 130 IPOL = 1, 2

         IF( IPOL.EQ.1 ) IPOLZN = -1234
         IF( IPOL.EQ.2 ) IPOLZN = +1234

         DO 120 NREFIN = 1, 2

            IF( NREFIN.EQ.1 ) CREFIN = ( 1.342, 0.0 )
            IF( NREFIN.EQ.2 ) CREFIN = ( 1.342, -1.0 )

            DO 110 ISZPAR = 1, 3

               IF( ISZPAR.EQ.1 ) XX = 1.E-4
               IF( ISZPAR.EQ.2 ) XX = 1.0
               IF( ISZPAR.EQ.3 ) XX = 100.0

c                                           ** Calculate max. possible
c                                           ** no. of Legendre coeffs.
c                                           ** (twice Mie series length)

               NMOM = 2.*( XX + 4.*XX**(1./3.) + 2.)

               DO 100 IANG = 1, 2
c                                             ** WARNING -- NUMANG must
c                                             ** be odd for Simpson Rule
                  IF( IANG.EQ.1 ) NUMANG = 1001
                  IF( IANG.EQ.2 ) NUMANG = 10001

                  IF( NUMANG.GT.MAXANG ) CALL ERRMSG
     &                ( 'PMOMTEST--MAXANG too small', .TRUE.)

                  WRITE( *, 9001 ) NUMANG, XX, IPOLZN, CREFIN

                  STEP   = 2./ ( NUMANG - 1 )

                  DO 10 N = 1, NUMANG
                     XMU( N ) = 1. - ( N - 1 )*STEP
   10             CONTINUE

                  CALL MIEV0( XX, CREFIN, PERFCT, MIMCUT, ANYANG,
     &                        NUMANG, XMU, NMOM, IPOLZN, MOMDIM, PRNT,
     &                        QEXT, QSCA, GQSC, PMOM, SFORW, SBACK, S1,
     &                        S2, TFORW, TBACK, SPIKE )

                  IF( IPOLZN.LT.0 ) THEN
c                                          ** Use S1,S2 to store T1,T2
c                                          ** (see MIEV0 doc for defns)

                     DO 20 N = 2, NUMANG - 1
                        CTMP   = S1( N )
                        S1( N ) = ( S2( N ) - XMU( N )*S1( N ) ) /
     &                            ( 1. - XMU( N )**2 )
                        S2( N ) = ( CTMP - XMU( N )*S2( N ) ) /
     &                            ( 1. - XMU( N )**2 )
   20                CONTINUE

                     S1( 1 ) = TFORW( 1 )
                     S2( 1 ) = TFORW( 2 )
                     S1( NUMANG ) = TBACK( 1 )
                     S2( NUMANG ) = TBACK( 2 )

                  END IF
c                                       ** Calculate phase quantities

                  DO 30 N = 1, NUMANG
                     PQ( 1, N ) = SQ( S1( N ) )
                     PQ( 2, N ) = SQ( S2( N ) )
                     PQ( 3, N ) =  REAL(  S1( N )*CONJG( S2( N ) ) )
                     PQ( 4, N ) = -AIMAG( S1( N )*CONJG( S2( N ) ) )
   30             CONTINUE

                  FNORM  = 4./ ( XX**2 * QSCA )

                  DO 90 MOM = 0, NMOM
c                                         ** Calculate Legendre polys.
                     IF( MOM.EQ.0 ) THEN

                        DO 40 N = 1, NUMANG
                           PL( N ) = 1.0
   40                   CONTINUE

                     ELSE IF( MOM.EQ.1 ) THEN

                        DO 50 N = 1, NUMANG
                           PLM1( N ) = 1.0
                           PL( N ) = XMU( N )
   50                   CONTINUE

                     ELSE

                        DO 60 N = 1, NUMANG
                           POL    = ( ( 2*MOM - 1 )*XMU( N )*PL( N ) -
     &                              ( MOM - 1 )*PLM1( N ) ) / MOM
                           PLM1( N ) = PL( N )
                           PL( N ) = POL
   60                   CONTINUE

                     END IF

                     DO 80 J = 1, 4

                        DO 70 N = 1, NUMANG
                           PROD( N ) = PL( N )*PQ( J, N )
   70                   CONTINUE
c                                            ** Legendre moments by
c                                               angular quadrature

                        PMOMAP = 0.5 * SIMPSO( PROD, NUMANG, STEP )

c                                        ** Put moment number and polar-
c                                           ization flag into a string

                        IF(J.EQ.1) WRITE( PMINDX,'(I4,A,I1)') MOM,',',J
                        IF(J.NE.1) WRITE( PMINDX,'(I6)') J

                        IF( PMOM(MOM,J).NE.0.) THEN

                           WRITE(*,'(A,1P,2E20.5,0P,F24.5)')
     &                           PMINDX, PMOM(MOM,J),
     &                           ( 2*MOM + 1 ) * FNORM * PMOM(MOM,J),
     &                           PMOMAP / PMOM(MOM,J) - 1.

                        ELSE

                           WRITE(*,'(A,F20.1)') PMINDX, 0.0

                        END IF

   80                CONTINUE

   90             CONTINUE

  100          CONTINUE

  110       CONTINUE

  120    CONTINUE

  130 CONTINUE

      STOP

 9000 FORMAT( ////, 30X, '************************************',/,
     &              30X, '* L E G E N D R E    M O M E N T S *',/,
     &              30X, '************************************' )
 9001 FORMAT( //, ' No. Quad Angs=', I5, '   Size Param=', F8.4,
     &  '  IPOLZN=', I6, '   Refr Index=', F6.3, 1P,E9.1, /,'  PMOM',8X,
     &  'unnormalized     Dave normalized    quadrature/exact - 1')
      END

      REAL FUNCTION SIMPSO( F, N, H )

c         Quadrature by Simpsons Rule

c     I N P U T

c          F(I), I=1 TO N : Values of the function being integrated
c          N :  No. of function values (**MUST BE ODD**)
c          H :  Step (in variable of integration) between values of F

c     O U T P U T

c          SIMPSO :  Approximate value of integral  ( H/3 * ( F(1) +
c                        4F(2) + 2F(3) +...+ 4F(N-1) + F(N)) )

c ----------------------------------------------------------------------

c     .. Scalar Arguments ..

      INTEGER   N
      REAL      H
c     ..
c     .. Array Arguments ..

      REAL      F( * )
c     ..
c     .. Local Scalars ..

      INTEGER   I
c     ..
c     .. External Subroutines ..

      EXTERNAL  ERRMSG
c     ..
c     .. Intrinsic Functions ..

      INTRINSIC MOD
c     ..

      IF( N.LT.3 .OR. MOD( N,2 ).NE.1 .OR. H.LE.0.0 )
     &    CALL ERRMSG( 'SIMPSO--Bad input arguments',.TRUE. )

      SIMPSO = F( 1 ) + F( N )

      DO 10 I = 2, N - 1
         SIMPSO = SIMPSO + 2.*( 2 - MOD( I,2 ) ) * F( I )
   10 CONTINUE

      SIMPSO = ( H / 3.) * SIMPSO

      RETURN
      END

max_stars_repo_name: peterkuma/alcf

max_issues_repo_name: peterkuma/alcf

max_forks_repo_name: peterkuma/alcf


text:
C=======================================================================
C  FOR_OPMOB, Subroutine, S. J. Rymph
C-----------------------------------------------------------------------
C  Generates output file for daily growth variables
C-----------------------------------------------------------------------
C  REVISION       HISTORY
C  10/13/2005 SJR  Written
C-----------------------------------------------------------------------
C  Called by: CROPGRO
C  Calls:     None
!=======================================================================
      SUBROUTINE FOR_OPMOB(CONTROL, ISWITCH,
     &  YRPLT, MDATE, DAS, YRDOY, DTX, DXR57, PGAVL, NAVL, PG, PPMFAC,
     &  NMOBR, NMOBSR, MAINR, ASMDOT, RSPNO3, RSPNH4, RPRO,
     &  CNOD, CGRSD, CGRSH, CADVG, CSAVEV, AGRVG, PCH2O, WTLF, WLDOT,
     &  WLIDOT, WLFDOT, WLDOTN, RTWT, WRDOT, WRDOTN, STMWT, WSDOT,
     &  WSIDOT, WSFDOT, WSDOTN,
     &  STRWT, WSRDOT, WSRIDOT, WSRFDOT, WSRDOTN,
     &  SLMDOT, LFSENWT, LTSEN, SLNDOT, SLCADDOT, SLNADDOT,
     &  SLDOT, SRMDOT, SRNDOT, SRCADDOT, SRNADDOT, SRDOT,
     &  SSMDOT, STSENWT, STLTSEN, SSNDOT, SSCADDOT, SSNADDOT,
     &  SSDOT, SSRMDOT, SSRNDOT, SSRCADDOT, SSRNADDOT, SSRDOT,
     &  NDMVEG, NDMNEW, NDMOLD, NDMTOT, NMINEP, NMINEA, SDNPL, TRNU,
     &  NFIXN, NGRSD, NGRSH, TSNMOB,
     &  WTNLF, WNRLF, LFSNMOB, NMINELF, ANMINELF,
     &  NGRLF, NLALL, NADLF, NRUSLF, PCNL,
     &  WTNRT, WNRRT, RTSNMOB, NMINERT, ANMINERT,
     &  NGRRT, NRALL, NADRT, NRUSRT, PCNRT,
     &  WTNST, WNRST, STSNMOB, NMINEST, ANMINEST,
     &  NGRST, NSALL, NADST, NRUSST, PCNST,
     &  WTNSR, WNRSR, SRSNMOB, NMINESR, ANMINESR,
     &  NGRSR, NSRALL, NADSR, NRUSSR, PCNSR,
     &  WTNSH, WNRSH, SHNMINE, NRUSSH,
     &  TSCMOB, WCRLF, WRCLDT, LFSCMOB, CMINELF, CADLF, CRUSLF, RHOL,
     &  WCRRT, WRCRDT, RTSCMOB, CMINERT, CADRT, CRUSRT, RHOR,
     &  WCRST, WRCSDT, STSCMOB, CMINEST, CADST, CRUSST, RHOS,
     &  WCRSR, WRCSRDT, SRSCMOB, CMINESR, CADSR, CRUSSR, RHOSR,
     &  WCRSH, WRCSHD, SHCMINE, CRUSSH, CHORECOVER, NLKSPENT, NLKNUSED,
     &  NLKCHK, TNLKCHK, CMOBSR, LAIMOBR, VNMOBR)

!-----------------------------------------------------------------------
      USE ModuleDefs     !Definitions of constructed variable types,
        ! which contain control information, soil
        ! parameters, hourly weather data.
      IMPLICIT NONE
      SAVE
!-----------------------------------------------------------------------
      CHARACTER*1  IDETG, RNMODE
      CHARACTER*2  CROP
      CHARACTER*6, PARAMETER :: ERRKEY = 'OPMOB1'
!      CHARACTER*8  FNAME
      CHARACTER*12 OUTMOB

      CHARACTER*30 FILEIO

      INTEGER DAP, DAS, DOY, DYNAMIC, ERRNUM, FROP, I
      INTEGER RUN
      INTEGER NOUTMOB
      INTEGER MDATE, TIMDIF, YEAR, YRDOY, YRPLT
      REAL DTX, DXR57, PPMFAC, NMOBR, NMOBSR
      REAL WTLF, WLDOT, WLDOTN, RTWT, WRDOT, WRDOTN, STMWT, WSDOT,
     &  WSDOTN, STRWT, WSRDOT, WSRDOTN,
     &  SLMDOT, LFSENWT, LTSEN, SLNDOT, SLCADDOT, SLNADDOT,
     &  SLDOT, SRMDOT, SRNDOT, SRCADDOT, SRNADDOT, SRDOT,
     &  SSMDOT, STSENWT, STLTSEN, SSNDOT, SSCADDOT, SSNADDOT,
     &  SSDOT, SSRMDOT, SSRNDOT, SSRCADDOT, SSRNADDOT, SSRDOT,
     &  NDMVEG, NDMNEW, NDMOLD, NDMTOT, NMINEP, NMINEA, TRNU, TSNMOB,
     &  WTNLF, WNRLF, LFSNMOB, NMINELF, ANMINELF, NRUSLF, PCNL,
     &  WTNRT, WNRRT, RTSNMOB, NMINERT, ANMINERT, NRUSRT, PCNRT,
     &  WTNST, WNRST, STSNMOB, NMINEST, ANMINEST, NRUSST, PCNST,
     &  WTNSR, WNRSR, SRSNMOB, NMINESR, ANMINESR, NRUSSR, PCNSR,
     &  WTNSH, WNRSH, SHNMINE, NRUSSH,
     &  TSCMOB, WCRLF, WRCLDT, LFSCMOB, CMINELF, CRUSLF, RHOL,
     &  WCRRT, WRCRDT, RTSCMOB, CMINERT, CRUSRT, RHOR,
     &  WCRST, WRCSDT, STSCMOB, CMINEST, CRUSST, RHOS,
     &  WCRSR, WRCSRDT, SRSCMOB, CMINESR, CRUSSR, RHOSR,
     &  WCRSH, WRCSHD, SHCMINE, CRUSSH, WLIDOT, WLFDOT,
     &  WSIDOT, WSFDOT, WSRIDOT, WSRFDOT, PG, PGAVL, NAVL,
     &  MAINR, ASMDOT, RSPNO3, RSPNH4, RPRO, CNOD, CGRSD, CGRSH,
     &  CADVG, CSAVEV, AGRVG, PCH2O, SDNPL, NFIXN,
     &  CADLF, CADRT, CADST, CADSR, NADLF, NADRT, NADST, NADSR,
     &  NGRSD, NGRSH, NGRLF, NGRRT, NGRST, NGRSR,
     &  NLALL, NRALL, NSALL, NSRALL, CHORECOVER, NLKSPENT, NLKNUSED,
     &  NLKCHK, TNLKCHK, CMOBSR, LAIMOBR, VNMOBR

C-------------------------------------------

      LOGICAL FEXIST, FIRST

!-----------------------------------------------------------------------
!     Define constructed variable types based on definitions in
!     ModuleDefs.for.

!     The variable "CONTROL" is of type "ControlType".
      TYPE (ControlType) CONTROL

!     The variable "ISWITCH" is of type "SwitchType".
      TYPE (SwitchType) ISWITCH

!     No output for fallow crop
      CROP    = CONTROL % CROP
      IDETG   = ISWITCH % IDETG
      IF (CROP .EQ. 'FA' .OR. IDETG .EQ. 'N') RETURN

!     Transfer values from constructed data types into local variables.
      DAS     = CONTROL % DAS
      DYNAMIC = CONTROL % DYNAMIC
      FILEIO  = CONTROL % FILEIO
      FROP    = CONTROL % FROP
!      LUNIO   = CONTROL % LUNIO
      RUN     = CONTROL % RUN
      RNMODE  = CONTROL % RNMODE
      YRDOY   = CONTROL % YRDOY
!      YRSIM   = CONTROL % YRSIM

      IDETG = ISWITCH % IDETG

!***********************************************************************
!***********************************************************************
!     Run initialization - run once per simulation
!***********************************************************************
      IF (DYNAMIC .EQ. RUNINIT) THEN
!-----------------------------------------------------------------------
!      IF (IDETG .EQ. 'Y') THEN
        OUTMOB  = 'PlantMob.OUT'
        CALL GETLUN('OUTMOB', NOUTMOB)

!      ENDIF

!***********************************************************************
!***********************************************************************
!     Seasonal initialization - run once per season
!***********************************************************************
      ELSEIF (DYNAMIC .EQ. SEASINIT) THEN
!-----------------------------------------------------------------------
!       Initialize daily growth output file
        INQUIRE (FILE = OUTMOB, EXIST = FEXIST)
        IF (FEXIST) THEN
        OPEN (UNIT = NOUTMOB, FILE = OUTMOB, STATUS = 'OLD',
     &  IOSTAT = ERRNUM, ACCESS = 'APPEND')
        FIRST = .FALSE.
        ELSE
        OPEN (UNIT = NOUTMOB, FILE = OUTMOB, STATUS = 'NEW',
     &  IOSTAT = ERRNUM)
        WRITE(NOUTMOB,'("*SENESCENCE/MOBILIZATION OUTPUT FILE")')
        FIRST = .TRUE.
        ENDIF

        !Write headers
!        CALL HEADER(SEASINIT, FILEIO, NOUTMOB, RUN)
        CALL HEADER(SEASINIT, NOUTMOB, RUN)
        WRITE (NOUTMOB,210)
210   FORMAT('@YEAR DOY   DAS   DAP   ',
     &  'DTX DXR57  PPMF NMOBR NMBSR    ',
     &  'PG PGAVL MAINR ASMDT RSPN3 RSPN4  RPRO  CNOD ',
     &  'CGRSD CGRSH CADVG CSAVE AGRVG PCH2O   ',
     &  'WTLF  WLDOT  WLDTN   RTWT WRDOT  WRDTN  STMWT  WSDOT  ',
     &  'WSDTN   STRWT  WSRDT  WSRDN ',
     &  'SLMDT LSNWT LLTSN SLNDT SLCAD SLNAD ',
     &  'SLDOT SRMDT  SRNDT SRCAD SRNAD  SRDOT ',
     &  'SSMDT SSNWT SLTSN SSNDT SSCAD SSNAD ',
     &  'SSDOT  SSRMD SSRND SSRCA SSRNA  SSRDT ',
     &  'NDMVG NDMNW NDMOD NDMTT NMINP NMINA  TRNU SDNPL NFIXN TSNMO   ',
     &  'NAVL NGRSD NGRSH  ',
     &  'WTNLF  WNRLF NGRLF   WLIDT WLFDT ',
     &  'LFSNM NMNLF ANMLF NLALL NADLF NRULF  PCNL   ',
     &  'WTNRT WNRRT NGRRT RTSNM NMNRT ANMRT NRALL NADRT NRURT PCNRT ',
     &  'WTNST WNRST NGRST  WSIDT WSFDT ',
     &  'STSNM NMNST ANMST NSALL NADST NRUST PCNST   ',
     &  'WTNSR WNRSR NGRSR  WSRID WSRFD ',
     &  'SRSNM NMNSR ANMSR NSRAL NADSR NRUSR PCNSR ',
     &  'WTNSH WNRSH SHNMN NRUSH ',
     &  'TSCMO WCRLF   WCLDT LFSCM CMNLF CADLF CRULF  RHOL ',
     &  'WCRRT  WCRDT RTSCM CMNRT CADRT CRURT  RHOR ',
     &  'WCRST  WCSDT STSCM CMNST CADST CRUST  RHOS ',
     &  'WCRSR WCSRD SRSCM  CMNSR CADSR  CRUSR RHOSR ',
     &  'WCRSH WCSHD SHCMN CRUSH  CHORC  NLKSP  NLKNU  ',
     &  'LKCHK   TLKCK CMBSR LAIMB VNMOB')

!***********************************************************************
!***********************************************************************
!     DAILY OUTPUT
!***********************************************************************
      ELSE IF (DYNAMIC .EQ. OUTPUT) THEN

!-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C   CHECK FOR OUTPUT FREQUENCY
C-----------------------------------------------------------------------
      IF (YRDOY .LT. YRPLT .OR. YRPLT .LT. 0) RETURN
      IF ((MOD(DAS,FROP) .EQ. 0)          !Daily output every FROP days,
     &  .OR. (YRDOY .EQ. YRPLT)           !on planting date, and
     &  .OR. (YRDOY .EQ. MDATE)) THEN     !at harvest maturity

!       Print
        DAP = MAX(0,TIMDIF(YRPLT,YRDOY))
        CALL YR_DOY(YRDOY, YEAR, DOY)

        WRITE (NOUTMOB,310)YEAR, DOY, DAS, DAP,
     &  DTX, DXR57, PPMFAC, NMOBR, NMOBSR,
     &  PG, PGAVL, MAINR, ASMDOT, RSPNO3, RSPNH4, RPRO, CNOD,
     &  CGRSD, CGRSH, CADVG, CSAVEV, AGRVG, PCH2O,
     &  WTLF, WLDOT, WLDOTN, RTWT, WRDOT, WRDOTN, STMWT, WSDOT,
     &  WSDOTN, STRWT, WSRDOT, WSRDOTN,
     &  SLMDOT, LFSENWT, LTSEN, SLNDOT, SLCADDOT, SLNADDOT,
     &  SLDOT, SRMDOT, SRNDOT, SRCADDOT, SRNADDOT, SRDOT,
     &  SSMDOT, STSENWT, STLTSEN, SSNDOT, SSCADDOT, SSNADDOT,
     &  SSDOT, SSRMDOT, SSRNDOT, SSRCADDOT, SSRNADDOT, SSRDOT,
     &  NDMVEG, NDMNEW, NDMOLD, NDMTOT, NMINEP, NMINEA, TRNU,
     &  SDNPL, NFIXN, TSNMOB, NAVL, NGRSD, NGRSH,
     &  WTNLF, WNRLF, NGRLF, WLIDOT, WLFDOT,
     &  LFSNMOB, NMINELF, ANMINELF, NLALL, NADLF, NRUSLF, PCNL,
     &  WTNRT, WNRRT, NGRRT,
     &  RTSNMOB, NMINERT, ANMINERT, NRALL, NADRT, NRUSRT, PCNRT,
     &  WTNST, WNRST, NGRST, WSIDOT, WSFDOT,
     &  STSNMOB, NMINEST, ANMINEST, NSALL, NADST, NRUSST, PCNST,
     &  WTNSR, WNRSR, NGRSR, WSRIDOT, WSRFDOT,
     &  SRSNMOB, NMINESR, ANMINESR, NSRALL, NADSR, NRUSSR, PCNSR,
     &  WTNSH, WNRSH, SHNMINE, NRUSSH,
     &  TSCMOB, WCRLF, WRCLDT, LFSCMOB, CMINELF, CADLF, CRUSLF, RHOL,
     &  WCRRT, WRCRDT, RTSCMOB, CMINERT, CADRT, CRUSRT, RHOR,
     &  WCRST, WRCSDT, STSCMOB, CMINEST, CADST, CRUSST, RHOS,
     &  WCRSR, WRCSRDT, SRSCMOB, CMINESR, CADSR, CRUSSR, RHOSR,
     &  WCRSH, WRCSHD, SHCMINE, CRUSSH, CHORECOVER, NLKSPENT, NLKNUSED,
     &  NLKCHK, TNLKCHK,
     &  CMOBSR, LAIMOBR, VNMOBR
310       FORMAT (1X,I4,1X,I3.3,2(1X,I5),2(1X,F5.1),3(1X,F5.3),
     &    1X,F5.2,1X,F5.3,2(1X,F5.3),2(1X,F5.2),8(1X,F5.3),1X,F6.2,
     &    1X,F6.1,1X,F6.3,1X,F6.1,1X,F5.1,1X,F6.3,1X,F6.2,1X,F6.1,
     &    1X,F6.3,1X,F7.2,1X,F6.2,1X,F6.3,8(1X,F5.3),1X,F6.3,
     &    2(1X,F5.3),1X,F6.3,7(1X,F5.3),1X,F6.3,3(1X,F5.3),1X,F6.3,
     &    10(1X,F5.3),1X,F6.3,2(1X,F5.3),2(1X,F6.3),1X,F5.3,1X,F7.3,
     &    1X,F5.2,7(1X,F5.3),1X,F7.3,1X,F5.2,11(1X,F5.3),1X,F6.2,
     &    1X,F5.2,7(1X,F5.3),1X,F7.3,1X,F5.2,1X,F5.3,1X,F6.3,
     &    13(1X,F5.3),1X,F5.2,1X,F7.3,5(1X,F5.3),1X,F5.1,1X,F6.3,
     &    5(1X,F5.3),1X,F5.2,1X,F6.2,5(1X,F5.3),1X,F5.1, 1X,F5.2,
     &    1X,F5.3,1X,F6.3, 1X,F5.3, 1X,F6.3,5(1X,F5.3), 3(1X,F6.4),
     &    1X,F6.4,1X,F7.3,3(1X,F5.3))

      ENDIF

!***********************************************************************
!***********************************************************************
!     Seasonal Output
!***********************************************************************
      ELSE IF (DYNAMIC .EQ. SEASEND) THEN
C-----------------------------------------------------------------------

        !Close daily output files.
        CLOSE (NOUTMOB)

!***********************************************************************
!***********************************************************************
!     END OF DYNAMIC IF CONSTRUCT
!***********************************************************************
      ENDIF
!***********************************************************************
      RETURN
      END ! SUBROUTINE FOR_OPMOB
!=======================================================================

!=======================================================================
!       Variable definitions for FOR_OPMOB
!-----------------------------------------------------------------------
!***********************************************************************
!       END SUBROUTINE FOR_OPMOB
!=======================================================================

ext: for

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:
program writefloats
  implicit none

  real, dimension(10) :: a, sqrta
  integer :: i
  integer, parameter :: unit = 40

  a = (/ (i, i=1,10) /)
  sqrta = sqrt(a)

  open(unit, file="xydata.txt", status="new", action="write")
  call writexy(unit, a, sqrta)
  close(unit)

contains

  subroutine writexy(u, x, y)
    real, dimension(:), intent(in) :: x, y
    integer, intent(in) :: u

    integer :: i

    write(u, "(2F10.4)") (x(i), y(i), i=lbound(x,1), ubound(x,1))
  end subroutine writexy

end program writefloats

ext: f

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
Property Manager, Redwood Grove.  Phone:  5302048669.  Email:  Mailto(changcl AT comcast DOT net).

20100128 19:53:18 nbsp Thanks for editing under your real name. I assume youve already seen Welcome to the Wiki/Business Owner? Users/JasonAller

20110610 16:11:40 nbsp Hey Christina  feel free to just remove information that is incorrect.  No need to just put a line through it if its untrue (giving a reason in the edit comment helps prevent people from misinterpreting or not understanding the edit and restoring the info).  The wiki tries to document Davis, and sometimes that includes history, but functional information is pretty much just intended to be informative. Users/JabberWokky Evan JabberWokky Edwards, jw.dw@timewarp.org

20110610 16:14:33 nbsp Ah!  I see  I didnt realize what that page was for!  You were right, I was wrong.  Take my offered advice as a mistaken but well intentioned attempt to be helpful. Users/JabberWokky

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
Alison Keenan is a Graduate Students graduate student at UC Davis, starting in her PhD in the Nutritional Biology program. She moved here from New York City, via Amtrak train.
Her undergraduate degree is in Forensic Science, specifically toxicology. Not only is she a scientist, but also a writer and artist of sorts. Currently she lives in a studio with her two cats Miller and Sage.

Alison is a GSA representative for the Executive Council of the Academic Senate, as well as the Educational Policy Committee. She has responsibilities as a member of the Student Advisory Committee for the Graduate Group in Nutritional Biology.

She is also a DJ at 95.7 KDRT the local grassroots station in Davis. Her show East Coast mashup can be also found at http://www.kdrt.org/station/archives/398

From California State Fair in Sacramento.


text:
!  *********************************************************************
!  *                                                                   *
!  *                         subroutine dcvit3                         *
!  *                                                                   *
!  *********************************************************************
!  Double Precision Version 2.0
!  Written by Gordon A. Fenton, TUNS, Aug. 2, 1993
!
!  PURPOSE  computes the initial k1*k2*k3 x k1*k2*k3 covariance matrix (and a
!           27x27 submatrix) between local averages. Used by LAS3G.
!
!
!  This routine forms the covariance matrix R0 between local averages arranged
!  on a k1 x k2 x k3 grid using a user provided variance function "dvarfn" that
!  returns the variance of a local average (note that this is the point
!  variance sigma^2 times the traditionally "variance" function, as defined
!  by Vanmarcke in "Random Fields", MIT Press, 1984). The grid and element
!  numbering schemes appear as follows (for k1 = k2 = 5, k3 = 1 plane shown)
!
!          | Dx |
!     ---  --------------------------
!      Dy  | 21 | 22 | 23 | 24 | 25 |
!     ---  --------------------------
!          | 16 | 17 | 18 | 19 | 20 |
!          --------------------------           ----------------
!          | 11 | 12 | 13 | 14 | 15 |           |  7 |  8 |  9 |
!          --------------------------           ----------------
!          |  6 |  7 |  8 |  9 | 10 |           |  4 |  5 |  6 |
!          --------------------------           ----------------
!          |  1 |  2 |  3 |  4 |  5 |           |  1 |  2 |  3 |
!          --------------------------           ----------------
!
!                  R0 Array                          R Array
!
!  with the k1 x k2 array numbering on the left and a basic 3 x 3 subarray
!  shown on the right. In the depth direction (not shown), the numbering
!  continues; ie the element immediately behind element 1 on the left is
!  element number 26, etc. If we call Z_i the local average of a random process
!  over the i'th cell, then for E[Z] = 0, the elements of R0 are defined by
!
!           R0(i,j) = E[Z_i*Z_j]
!
!  Note that the elements of R are simply extracted directly from the
!  appropriate elements of R0 (for example, R(1,5) = R0(1,7)) if k1, k2 and k3
!  are greater than 2. Note also that since both R0 and R are symmetric, ONLY
!  THE UPPER TRIANGULAR VALUES are stored herein. Finally note that
!  the random process is assumed to be quadrant symmetric (that is the
!  covariance function is even in each individual lag component) to reduce
!  the total number of covariance calculations.
!
!  Arguments to this routine are as follows;
!
!    dvarfn     external double precision function which returns the
!               variance of a local average of size X x Y x Z. This function
!               is referenced using the call "V = dvarfn(X,Y,Z)" (other
!               parameters to the function must be passed via common blocks).
!
!    R0          double precision array of size at least k1*k2*k3 x k1*k2*k3
!               which on output will contain the covariance matrix between
!               the local average elements shown above on the left. Note
!               that since R0 is symmetric, only the upper triangular values
!               are actually stored (compatible with DCHOL2). (output)
!
!    iq         leading dimension of the array R0 as specified in the calling
!               routine. (input)
!
!    R          double precision array of size at least 27 x 27 which on output
!               will contain the covariance matrix between the local average
!               elements shown above on the right. Note that since R is
!               symmetric, only the upper triangular values are actually
!               stored (compatible with DSIFA). (output)
!
!    ir         leading dimension of the array R as specified in the calling
!               routine. (input)
!
!    k1         number of local average cells in the x direction
!               (horizontally). (input)
!
!    k2         number of local average cells in the y direction
!               (vertically). (input)
!
!    k3         number of local average cells in the z direction
!               (depth). (input)
!
!    Dx         x-dimension of the local average cells. (input)
!
!    Dy         y-dimension of the local average cells. (input)
!
!    Dz         z-dimension of the local average cells. (input)
!
!  Required:
!   1) from libGAFsim:	DCVAA3
!----------------------------------------------------------------------------
      subroutine dcvit3( dvarfn, R0, iq, R, ir, k1, k2, k3, Dx, Dy, Dz )
      implicit real*8 (a-h,o-z)
      dimension R0(iq,*), R(ir,*)
      external dvarfn
!					first form the essential elements of R0
      ii = 0
      do 30 k = 1, k3
         tz = dble(k-1)
         do 20 j = 1, k2
            ty = dble(j-1)
            do 10 i = 1, k1
               ii = ii + 1
               tx = dble(i-1)
               R0(1,ii) = dcvaa3( dvarfn, Dx, Dy, Dz, tx, ty, tz )
  10        continue
  20     continue
  30  continue
!					now distribute into the upper triangle
      kk  = k1*k2*k3
      k12 = k1*k2
      do 50 n = 2, kk
         n3 = (n-1)/k12
         n2 = (n - k12*n3 - 1)/k1
         do 40 m = 2, n
            m3 = (m-1)/k12
            m2 = (m - k12*m3 - 1)/k1
            i3 = m3 - n3
            i2 = m2 - n2
            i1 = (m-n) - k1*i2 - k12*i3
            k  = 1 + iabs(i1) + k1*iabs(i2) + k12*iabs(i3)
            R0(m,n) = R0(1,k)
  40     continue
  50  continue
!					form R (27 x 27)

      if( k1 .lt. 3 .or. k2 .lt. 3 .or. k3 .lt. 3 ) then
!						find essential elements of R
         ii = 0
         do 80 k = 1, 3
            tz = dble(k-1)
            do 70 j = 1, 3
               ty = dble(j-1)
               do 60 i = 1, 3
                  tx = dble(i-1)
                  ii = ii + 1
                  R(1,ii) = dcvaa3( dvarfn, Dx, Dy, Dz, tx, ty, tz )
  60           continue
  70        continue
  80     continue
!						and distribute into upper tri.
         do 100 n = 2, 27
            n3 = (n-1)/9
            n2 = (n - 9*n3 - 1)/3
            do 90 m = 2, n
               m3 = (m-1)/9
               m2 = (m - 9*m3 - 1)/3
               i3 = m3 - n3
               i2 = m2 - n2
               i1 = (m-n) - 3*i2 - 9*i3
               k  = 1 + iabs(i1) + 3*iabs(i2) + 9*iabs(i3)
               R(m,n) = R(1,k)
  90        continue
 100     continue
      else
!						or extract R from R0
         do 120 n = 1, 27
            n3 = (n-1)/9
            n2 = (n - 9*n3 - 1)/3
            do 110 m = 1, n
               m3 = (m-1)/9
               m2 = (m - 9*m3 - 1)/3
               i3 = m3 - n3
               i2 = m2 - n2
               i1 = (m-n) - 3*i2 - 9*i3
               k  = 1 + iabs(i1) + k1*iabs(i2) + k12*iabs(i3)
               R(m,n) = R0(1,k)
 110        continue
 120     continue
      endif

      return
      end

max_stars_repo_name: leemargetts/PoreFEM

max_issues_repo_name: leemargetts/PoreFEM

max_forks_repo_name: leemargetts/PoreFEM


text:
! To compile and run: make

program q3
  use fitsio
  use integrator
  implicit none

  real, parameter :: pi = 3.1415926535897932384626433832795028841971693993751058
  real :: energy_flip_min

  energy_flip_min = min( &
    energy([pi, 0.0, 0.0, 0.0]), &
    energy([0.0, 0.0, pi, 0.0]))

  call main()

contains

  subroutine main()
    integer, parameter :: steps_per_second = 7
    real, parameter :: time_period = 10000.0 * sqrt(1.0 / 9.806)
    real, parameter :: dt = 1.0 / steps_per_second
    integer, parameter :: steps = steps_per_second * time_period
    real, parameter :: zoom_factor = 3.0
    character(len=64) :: filenames(5)
    real :: th1(2) = [-1, 0], th2(2) = [2, 3]
    integer :: i

    print "(a)", "Q3. Fractal generation"
    print *
    call write_fractal("plot_fractal_16x16.fit",     dt, steps, [-pi, pi], [-pi, pi], 16, 16)
    call write_fractal("plot_fractal_32x32.fit",     dt, steps, [-pi, pi], [-pi, pi], 32, 32)
    call write_fractal("plot_fractal_64x64.fit",     dt, steps, [-pi, pi], [-pi, pi], 64, 64)
    call write_fractal("plot_fractal_128x128.fit",   dt, steps, [-pi, pi], [-pi, pi], 128, 128)
    call write_fractal("plot_fractal_256x256.fit",   dt, steps, [-pi, pi], [-pi, pi], 256, 256)
    ! call write_fractal("plot_fractal_512x512.fit",   dt, steps, [-pi, pi], [-pi, pi], 512, 512)
    ! call write_fractal("plot_fractal_1024x1024.fit", dt, steps, [-pi, pi], [-pi, pi], 1024, 1024)
    ! call write_fractal("plot_fractal_2048x2048.fit", dt, steps, [-pi, pi], [-pi, pi], 2048, 2048)

    ! Recursively zoom into a specific region
    filenames(1) = "plot_fractal_zoom3x_256x256.fit"
    filenames(2) = "plot_fractal_zoom9x_256x256.fit"
    filenames(3) = "plot_fractal_zoom27x_256x256.fit"
    filenames(4) = "plot_fractal_zoom81x_256x256.fit"
    filenames(5) = "plot_fractal_zoom243x_256x256.fit"
    do i = 1, 5
      call write_fractal(filenames(i), dt, steps, th1, th2, 256, 256)
      call zoom_ranges(th1, th2, zoom_factor)
    end do
  end subroutine

  subroutine write_fractal(filename, dt, steps, th1, th2, width, height)
    !! Write phase plot within ranges given by th1 and th2
    character(len=*) :: filename
    character(len=256) :: cmd
    integer, parameter :: iters_since_msg_max = 1024
    integer :: steps, width, height, i, j, iters_since_msg
    real :: dt, th1(2), th2(2), data_(1, width, height)

    print "(a)", "Plotting fractal"
    print "(a, i4, a, i4)", "Dimensions: ", width, " x ", height
    print "(a, 2f6.2)", "theta_1 range: ", th1
    print "(a, 2f6.2)", "theta_2 range: ", th2

    iters_since_msg = iters_since_msg_max

    do j = 1, height
      if (iters_since_msg >= iters_since_msg_max) then
        print "(i3, a)", 100 * (j - 1) / (height - 1), "% complete"
        iters_since_msg = 0
      end if

      !$omp parallel do
      do i = 1, width
        data_(1, i, j) = find_flip(steps, dt, y0=[ &
          th1(1) + (th1(2) - th1(1)) * (i - 1) / (width  - 1), 0.0, &
          th2(1) + (th2(2) - th2(1)) * (j - 1) / (height - 1), 0.0])
      end do

      iters_since_msg = iters_since_msg + width
    end do

    print "(i3, a)", 100, "% complete"
    data_ = log(data_)
    call write2fits(filename, data_, th1, th2, &
      ['magnitude'], '($theta_1$,$theta_2$)')

    write (cmd, "(3a)") "python plot.py ", filename, " --plot-fractal"
    call execute_command_line(cmd)

    print "(2a)", "Written output to ", filename
    print *
  end subroutine

  function find_flip(n, dt, y0) result(i)
    !! Run simulation for given initial condition y0 at time step dt
    !! Returns steps taken until a flip occurs
    real, parameter :: g = 9.806
    real :: y0(4), y(4), dt
    integer :: n, i

    ! Exit if energy of system is too low to allow flip
    if (energy(y0) < energy_flip_min) then
      i = n + 1
      return
    end if

    y = y0
    do i = 1, n
      if ((abs(y(1)) > pi) .or. (abs(y(3)) > pi)) exit
      call gl10(y, dt)
    end do
  end function

  subroutine zoom_ranges(th1, th2, factor)
    real :: th1(2), th2(2), factor, avg(2), diff(2)
    avg  = [(th1(1) + th1(2)), (th2(1) + th2(2))] / 2
    diff = [(th1(2) - th1(1)), (th2(2) - th2(1))] / (factor * 2)
    th1 = [avg(1) - diff(1), avg(1) + diff(1)]
    th2 = [avg(2) - diff(2), avg(2) + diff(2)]
  end subroutine

end program q3

max_stars_repo_name: YodaEmbedding/experiments

max_issues_repo_name: YodaEmbedding/experiments

max_forks_repo_name: YodaEmbedding/experiments


text:
module loadingLib

use funcAux
implicit none

private  LinearPressure, ConstantPressure, SinusoidalPressure ,  SinusoidalPressureZero, adaptative, LinearPeriodic
public pressureLoad, chooseLoad

abstract interface
	subroutine PressureLoadLike(p,time,dt,loadPar) ! flexible load
		implicit none
		real*8 , intent(in) :: time,dt,loadPar(:)
		real*8 , intent(out) :: p
	end subroutine
end interface

procedure (PressureLoadLike), pointer :: pressureLoad => null ()

contains

subroutine setF(F,LoadPar,Time,dt,LoadProg,LoadType)
	Real*8 , intent(out):: F(:,:)
	Real*8 , intent(in):: Time , dt, LoadPar(:)
	integer , intent(in) :: LoadProg, LoadType
	Real*8 :: beta(2), normal(2), tangent(2) ,  theta, alpha , lamb0, lamb1
	Real*8 :: lamb , eps
	integer :: k,l, dirShear

	call chooseLoad(LoadProg)
	call pressureLoad(lamb,time,dt,loadPar)

	F = 0.0d0

	selectcase(LoadType)
		case(1)
			F(1,1) = lamb
!~ 			write(0,*) "lxx=",lamb
		case(2)
			F(2,2) = lamb
!~ 			write(0,*) "lyy=",lamb

		case(3)
			F(1,2) = lamb
!~ 			write(0,*) "lxy=",lamb
		case(4)
			F(2,1) = lamb
!~ 			write(0,*) "lyx=",lamb
		case(5)
			F(1,2) = 0.5d0*lamb
			F(2,1) = 0.5d0*lamb
		case(6) ! biaxial
			F(1,1) = 0.5d0*lamb
			F(2,2) = 0.5d0*lamb
		case(7) !! incompressible !!
			F(1,1) = lamb
			F(2,2) = -lamb/(1.0d0+lamb)
		case(8) !! incompressible
			F(1,1) = lamb
			F(1,2) = dsqrt(lamb)
			F(2,1) = F(1,2)
		case(9) !! incompressible
			F(1,1) = lamb
			F(2,1) = 0.5*lamb
		case(10) !! prescribed
			F(1,1) = 1.05d0
			F(1,2) = 0.01d0
			F(2,1) = -0.02d0
			F(2,2) = 0.99d0

		case(11) !! beta \otimes normal , where beta = lamb * normal
			theta = datan(0.0d0/1.0d0)
			normal(1) = -dsin(theta)
			normal(2) = dcos(theta)

			F(1,1) = lamb * normal(1)* normal(1)
			F(1,2) = lamb * normal(1)* normal(2)
			F(2,1) = lamb * normal(2)* normal(1)
			F(2,2) = lamb * normal(2)* normal(2)

		case(12) !! beta \otimes normal , where beta = lamb * tangent
			theta = datan(0.0d0/1.0d0)
			normal(1) = -dsin(theta)
			normal(2) = dcos(theta)
			tangent(1) = dcos(theta)
			tangent(2) = dsin(theta)

			F(1,1) = lamb * tangent(1)* normal(1)
			F(1,2) = lamb * tangent(1)* normal(2)
			F(2,1) = lamb * tangent(2)* normal(1)
			F(2,2) = lamb * tangent(2)* normal(2)

		case(13) !! beta \otimes normal , where beta = lamb * ( b_t * tangent + b_n * normal)
!~ 			alpha = datan(0.0d0) + PI!!
			alpha = 0.5d0*PI
			theta = datan(0.3d0/1.0d0)
			normal(1) = -dsin(theta)
			normal(2) = dcos(theta)
			beta(1) = dcos(theta + alpha)
			beta(2) = dsin(theta + alpha)

			F(1,1) = lamb * beta(1)* normal(1)
			F(1,2) = lamb * beta(1)* normal(2)
			F(2,1) = lamb * beta(2)* normal(1)
			F(2,2) = lamb * beta(2)* normal(2)

		case(14) !! beta \otimes normal , where beta angle and normal angle is given by LoadPar
			theta = LoadPar(5)*PI
			alpha = LoadPar(6)*PI
			normal(1) = dcos(theta)
			normal(2) = dsin(theta)
			beta(1) = dcos(alpha)
			beta(2) = dsin(alpha)

			F(1,1) = lamb * beta(1)* normal(1)
			F(1,2) = lamb * beta(1)* normal(2)
			F(2,1) = lamb * beta(2)* normal(1)
			F(2,2) = lamb * beta(2)* normal(2)

		case(15)
			lamb0 = LoadPar(5)

			F = 0.0d0

			if(lamb<lamb0) then
				F(1,1) = lamb
				F(2,1) = 0.0d0
			else
				F(1,1) = lamb0
				F(2,1) = lamb - lamb0
			end if

		case(16)
			eps = LoadPar(4)
			k = nint(LoadPar(5))
			l = nint(LoadPar(6))

			F(1,1) = lamb

			if(k>0 .and. l> 0) then
				F(k,l) = F(k,l) + eps
			end if

		case(17)
			lamb0 = LoadPar(5)
			lamb1 = LoadPar(6)

			F = 0.0d0

			theta = (lamb - lamb0)/(lamb1 - lamb0)

			if(theta<0.0d0) then
				F(1,1) = lamb
				F(2,1) = 0.0d0
			else if(theta < 1.0d0) then
				F(1,1) = lamb0 + theta*lamb1

			else ! > 1.0d0
				F(1,1) = lamb1
				F(2,1) = lamb - lamb1
			end if

			F(1,1) = (1.0d0 - theta)*lamb
			F(2,1) = theta*(lamb - lamb1)

		case default
			write(0,*) 'Loadtype not recognised ' , LoadType
			stop

	end select

!~ 	call numprint(F)

end subroutine

subroutine chooseLoad(LoadProg)
	integer, intent(in) :: LoadProg

	selectcase(LoadProg)
		case(1)
			pressureLoad => LinearPressure
		case(2)
			pressureLoad => ConstantPressure
		case(3)
			pressureLoad => SinusoidalPressure
		case(4)
			pressureLoad => SinusoidalPressureZero
		case(5)
			pressureLoad => adaptative
		case(6)
			pressureLoad => LinearRampPressure
		case(7)
			pressureLoad => LinearPeriodic
		case default
			write(0,*) 'LoadProg not recognised ' , LoadProg
			stop
	end select

end subroutine

subroutine LinearPeriodic(p,time,dt,loadPar) ! periodic load, linear periodic
	implicit none
	real*8 , intent(in) :: time,dt,loadPar(:)
	real*8 , intent(out) :: p
	real*8 :: T,pmax0, dpmax, pmin0, dpmin

	real*8 :: di, trel, p0, p1
	integer :: i

	T = loadPar(1)
	pmax0 = loadPar(2)
	dpmax = loadPar(3)
	pmin0 = loadPar(4)
	dpmin = loadPar(5)

	di = dint(time/T)
	trel = time - di*T
	i = int(di)

	if(mod(i,2) == 0) then ! even number
		p0 = pmin0 + dpmin*di
		p1 = pmax0 + dpmax*(di+1.0d0)
	else ! odd number
		p0 = pmax0 + dpmax*di
		p1 = pmin0 + dpmin*(di+1.0d0)
	end if

	p = p0 + (p1 - p0)*trel/T

	write(0,*) "pressure=", p

end subroutine

subroutine LinearPressure(p,time,dt,loadPar) ! flexible load
	! p(t0) = 0 , p(t1) = pmax, linearly
	implicit none
	real*8 , intent(in) :: time,dt,loadPar(:)
	real*8 , intent(out) :: p
	real*8 :: t0,t1, pmax

	pmax = loadPar(1)
	t0 = loadPar(2)
	t1 = loadPar(3)

	if(time<t1) then
		p = (time-t0)*pmax/(t1-t0)
	else
		p = pmax
	end if

	write(0,*) "pressure=", p

end subroutine

subroutine ConstantPressure(p,time,dt,loadPar) ! flexible load
	implicit none
	real*8 , intent(in) :: time,dt,loadPar(:)
	real*8 , intent(out) :: p

	p = loadPar(1)

end subroutine

subroutine SinusoidalPressure(p,time,dt,loadPar) ! flexible load
	implicit none
	real*8 , intent(in) :: time,dt,loadPar(:)
	real*8 , intent(out) :: p
	real*8 :: t0,t1, p_amp,p_ave, ncycles,phase

	p_ave = loadPar(1)
	p_amp = loadPar(2)
	t0 = loadPar(3)
	t1 = loadPar(4)
	ncycles = loadPar(5)
 	phase = loadPar(6)*2.0d0*PI

	t1 = (time-t0)/(t1 - t0) ! t1 in [0,1]

 	p = p_ave + p_amp*dsin(ncycles*2.0d0*PI*t1+phase)

end subroutine

subroutine SinusoidalPressureZero(p,time,dt,loadPar) ! flexible load
	implicit none
	real*8 , intent(in) :: time,dt,loadPar(:)
	real*8 , intent(out) :: p
	real*8 :: t0,t1, p_amp,p_ave, ncycles,phase

	p_ave = loadPar(1)
	p_amp = loadPar(2)
	t0 = loadPar(3)
	t1 = loadPar(4)
	ncycles = loadPar(5)

 	phase = dasin(-p_ave/p_amp) + t0*ncycles*2.0d0*PI/(t1-t0)

	t1 = (time-t0)/(t1 - t0) ! t1 in [0,1]

 	p = p_ave + p_amp*dsin(ncycles*2.0d0*PI*t1+phase)

end subroutine

subroutine adaptative(p,time,dt,loadPar)
	implicit none
	real*8 , intent(in) :: time,dt,loadPar(:)
	real*8 , intent(out) :: p
	real*8 :: tf,lf,tm,kmin,kini , den, a1, a2, a3

	tf = loadPar(1)
	lf = loadPar(2)
	tm = loadPar(3)
	kmin = loadPar(4)
	kini = loadPar(5)

	den = 3.0d0*(tm**2.0d0)*tf**2.0d0 - 2.0d0*tm*tf**3.0d0
	a2 = (3.0d0*(tm**2.0d0)*(lf-kini*tf) - (tf**3.0d0)*(kmin - kini))/den
	a3 = (-2.0d0*tm*(lf-kini*tf) + (tf**2.0d0)*(kmin - kini))/den
	a1 = kini

	p = a1*time + a2*time**2.0 + a3*time**3.0

end subroutine

subroutine LinearRampPressure(p,time,dt,loadPar) ! flexible load
	! p(t0) = p0 , p(t1) = p1, linearly
	implicit none
	real*8 , intent(in) :: time,dt,loadPar(:)
	real*8 , intent(out) :: p
	real*8 :: t0,t1, p0, p1

	p0 = loadPar(1)
	p1 = loadPar(2)
	t0 = loadPar(3)
	t1 = loadPar(4)

	if(time<t0) then
		p = p0
	else if(time<t1) then
		p = p0 + (time-t0)*(p1-p0)/(t1-t0)
	else
		p = p1
	end if

	write(0,*) "pressure=", p, p0, p1, t0, t1

end subroutine

end module

max_stars_repo_name: felipefr/Piola

max_issues_repo_name: felipefr/Piola

max_forks_repo_name: felipefr/Piola


text:
wiki:woodland:Judys Grinder Judys Grinders is a small chain of restaurant restaurants with locations in Downtown Davis and Woodland. The original Woodland restaurant has been a local favorite since 1977 & the Davis location shows lots of promise

 Serves
  Fish & Chips
  Burgers
  Grinders
  Chicken
  Sandwiches
  Teriyaki Bowl
  Salads
  Shakes
  Ice Cream

20070303 11:48:03 nbsp Ive eaten here a couple of times.  The selection of food is amazing, and everything Ive tried has been decent (not great, but decent).  The ambience of the place is kind of depressing though.  Last time I went they had old pictures of food that were dusty and looked like they were from the 80s.  Not appetizing.  But Ill go back occasionally just to try something different on the extensive menu. Users/AlexPomeranz

20070304 11:30:48 nbsp I went there once when I was a kid and had the best damn pastrami burger Ive ever had. Users/BradBenedict

20070518 12:15:17 nbsp I eat in woodlant every day since I work at Seminis now, and they have the largest menu selection I have ever seen.  They have the seafood menu of London Fish and Chips, the fallafel and gyro menu of Ali babas, they also have tons of grinders (subs), burgers, cheesesteaks and Chicken items.  Its absolutely insane, but prepare yourself for a 20 minute wait on your food. Users/AllanRae

20131106 13:06:04 nbsp I just had lunch at the Judys Grinders in Davis and that is the last time.  The burger was lousy and two tables were served before me.  I had to ask a server to check on my order and she said that they couldnt find it.   There were only five tables in the restaurant and there was no reason to make this kind of error.  Users/DanielPadron
 Seriously? You go to a restaurant on the first (or second?) day after its opening and you post a bad review? All restaurants deserve at least a couple of weeks to work out the kinks. Users/ScottLay

20131110 21:08:40 nbsp I ate here last night at around 5 oclock. Restaurant wasnt very busy. I didnt have any problems with the service, besides the guy at the register forgetting to give my friend his beer right when he ordered (He went up and asked and got it then). The food itself was solid and came out pretty quickly. I had the pastrami burger and enjoyed it. The fries werent spectacular, but came out  hot and crispy, which is better than a lot of places can do. Ill definitely come back to try out some different things on the menu. Also, be warned, the burgers, sandwiches, subs, etc dont come with fries unless theyre one of the specials on the board; thats definitely a downside, especially with a modest portion of fries costing $3 extra. Users/jrb

20131216 11:21:03 nbsp I visited the E Street location and was unimpressed. First, the cashier got my order wrong. I ordered with my girlfriend, who got a sandwich (by itself, with water). I then ordered a bleu cheese and guacamole burger, also just with water. The cashier got my girlfriends order right, but strangely charged me for a #4 burger combo (which comes with fries and a drink, for an extra $2). He also carelessly overcharged two other friends who were there, assuming (without asking) that they wanted the combo deal. Second, the cashier didnt ask how I wanted the burger cooked. It came out mediumwell. Maybe thats just the default at Judys, but most places charging $7+ for a burger cook it to order.

The place wasnt busy and the location has been open and running for more than a month, so they dont have an excuse for poor service. I expect it to be out of business within 6 months (just like Huong Lan Sandwiches and Luigis before them  this storefront is a true revolving door). My burger was good, but Davis already has several good burger options.  Users/AndrewDay

20131218 21:05:48 nbsp I tried their fish and chips tonight.  Two pieces of fish were about half an inch square by six inches long, battered and fried to dark brown.  The fish had almost no moisture left, which made me sad.  The platter came with fries, three battered/fried green beans, and a side salad.  Salad: spring mix, raw yellow summer squash, red onions, raw mushrooms (cut thick and browning badly), sliced pickled peperoncini, shredded carrots, black olives, underripe tomato, and a couple bits of bell pepper.  It was trying really hard to look nice, but the middle of December is just not the time for raw summerseason salad veggies.  Coleslaw would have been more successful, or Japanesestyle pickles.

De Veres has much nicer fish and chips.  Dumpling House has reasonable fish and chips, too.  Plutos has WAY nicer salad. Users/JudithTruman

20131221 22:54:42 nbsp Service was fine for me, food was ok, and I dont like writing neighsayer reviews, but Im just not sure how this place is going to survive in this location unless they offer some serious discounts. Their Cheeseburger with fries and tax but excluding tip and drink is around $11 making it one of the more expensive in Yolo County. I cant see a situation where I want to pay top dollar for a welldone burger. Maybe Plainfield Station but thats a little cheaper and also provides a unique experience. Users/MikeyCrews

I just went to Judys grinder for the first time. I got one of the larger burgers so it was $9 but it was really good. They look huge & are pretty tasty. Also multiple friends of mine noted how cheap the bottle beers were, with decent selection as well! Daubert

An update. Went again. Found out the shakes are amazing. The Reuben was good, a little small but good. Daubert

20140120 15:14:58 nbsp Just ate there.  I was a bit apprehensive after reading comments disparaging the Davis location, but I had already made plans to meet there with a friend.  All I can say is, it was delicious.  Had the Texas Gold Burger, friend had a bunless burger, his son had a reuben.  We all left happy and satisfied.  Service was quick, everything cooked perfectly.  It was worth the price.  I really dont know what people are complaining about...if its related to a new restaurant getting its footing, thats one thing, you need to cut some slack, but my experience was not what others have been posting here.  So many businesses have come and gone at that location, Id like to see someone do well there, and I think this place can make it happen.
 Users/RobertArosteguy

20140124 13:27:39 nbsp Its a bit pricey, but its a lot better than Burgers and Brew with their stale buns. The fries are a lot better here, too. Theyre not overcooked like B&Bs. INNOUT is a better than these burgers, but sometimes you want something heftier or you want different flavors (they have tons of options). 4/5 stars, minus 1 star for cost. Users/UCD4Me

20140209 15:01:13 nbsp I have been here several times, once to watch the 49ers get trounced by the Seahawks. This place is a secret place to watch sports when the graduate is packed). We have tried a variety of items, essentially all of which are distinctly mediocre, but large portions are a plus for those who want quantity over quality (which counts me out). Beer selection is OK. I predict this restaurant will follow the flock of doomed enterprises that have opened here over the past few years. I give it a year.   Users/MichaelPlotkin

20140214 10:23:57 nbsp Visited here this week for lunch.  It was about 11:30 and there were three of us waiting for service for about 5 minutes.  I ordered the Philly, Fries and a Soda which came to $11.15.  A bit on the high side but Id say the portion was large and it was well made.  Inside was nice and the TVs and beer selection were also good.  They didnt seem to have any happy hour specials which would make this much more attractive to hang out in.  The food is basically burgers, cold sandwiches and fried food.    The philly was really tasty, but it used bell peppers instead of cherry peppers and was missing the spice and acid normally found in a true Philly. The meat was good and again it was well prepared.  The fries were cooked perfect for my taste.

It wasnt very busy during the entire time I was there eating. Thats a bad sign. They need to bring folks into this place if they hope to survive.  Coupons. 2 for 1, Happy Hour specials are going to be key to drumming up business and conquering the curse of this spot.  The food was solid enough but again it was on the heavy side and I know I wont be eating here more than 1 or 2 times a month.

The pictures of food is always odd to me and especially since they have the same burger presented 2 or 3 times around the place.  Additionally they have an electronic picture of the same burger but it doesnt rotate and just stays on the same photo.  Not sure thats a good thing.  Customer service needs to improve a bit.  The incredibly noisy hand cart they use to truck food into from their shopping trips needs to be oiled.  Maybe turn the TVs to a few different shows instead of all showing the same thing.  Id give it a 6.5 out of 10 probably overall.   Users/loneshark


text:
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C   FFTPACK 5.0
C
C   Authors:  Paul N. Swarztrauber and Richard A. Valent
C
C   $Id$
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE MSNTF1(LOT,JUMP,N,INC,X,WSAVE,DSUM,XH,WORK,IER)
      REAL       X(INC,*)       ,WSAVE(*)   ,XH(LOT,*)
      DOUBLE PRECISION           DSUM(*)
      IER = 0
      LJ = (LOT-1)*JUMP+1
      IF (N-2) 101,102,103
 102  SSQRT3 = 1./SQRT(3.)
      DO 112 M=1,LJ,JUMP
         XHOLD = SSQRT3*(X(M,1)+X(M,2))
         X(M,2) = SSQRT3*(X(M,1)-X(M,2))
         X(M,1) = XHOLD
  112 CONTINUE
  101  GO TO 200
  103 NP1 = N+1
      NS2 = N/2
      DO 104 K=1,NS2
         KC = NP1-K
         M1 = 0
         DO 114 M=1,LJ,JUMP
         M1 = M1 + 1
         T1 = X(M,K)-X(M,KC)
         T2 = WSAVE(K)*(X(M,K)+X(M,KC))
         XH(M1,K+1) = T1+T2
         XH(M1,KC+1) = T2-T1
  114    CONTINUE
  104 CONTINUE
      MODN = MOD(N,2)
      IF (MODN .EQ. 0) GO TO 124
      M1 = 0
      DO 123 M=1,LJ,JUMP
         M1 = M1 + 1
         XH(M1,NS2+2) = 4.*X(M,NS2+1)
  123 CONTINUE
  124 DO 127 M=1,LOT
         XH(M,1) = 0.
  127 CONTINUE
      LNXH = LOT-1 + LOT*(NP1-1) + 1
      LNSV = NP1 + INT(LOG(REAL(NP1))) + 4
      LNWK = LOT*NP1
C
      CALL RFFTMF(LOT,1,NP1,LOT,XH,LNXH,WSAVE(NS2+1),LNSV,WORK,
     1            LNWK,IER1)
      IF (IER1 .NE. 0) THEN
        IER = 20
        CALL XERFFT ('MSNTF1',-5)
        GO TO 200
      ENDIF
C
      IF(MOD(NP1,2) .NE. 0) GO TO 30
      DO 20 M=1,LOT
      XH(M,NP1) = XH(M,NP1)+XH(M,NP1)
   20 CONTINUE
   30 SFNP1 = 1./FLOAT(NP1)
      M1 = 0
      DO 125 M=1,LJ,JUMP
         M1 = M1+1
         X(M,1) = .5*XH(M1,1)
         DSUM(M1) = X(M,1)
  125 CONTINUE
      DO 105 I=3,N,2
         M1 = 0
         DO 115 M=1,LJ,JUMP
            M1 = M1+1
            X(M,I-1) = .5*XH(M1,I)
            DSUM(M1) = DSUM(M1)+.5*XH(M1,I-1)
            X(M,I) = DSUM(M1)
  115    CONTINUE
  105 CONTINUE
      IF (MODN .NE. 0) GO TO 200
      M1 = 0
      DO 116 M=1,LJ,JUMP
         M1 = M1+1
         X(M,N) = .5*XH(M1,N+1)
  116 CONTINUE
  200 CONTINUE
      RETURN
      END

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
Kerr Hall is located on the west side of the Quad, between Wellman Hall and Wickson Hall.  It was built in 1969. This building formerly housed the offices of the Davis Honors Challenge which moved to Surge III, and the Mathematics and Statistics departments, which moved to the new Mathematical Sciences Building. The question of the new occupants of Kerr Hall existed for some time, however the latest word is that Kerr Hall will have the Political Science Department bestowed upon it. Kerr Hall is named for wiki:WikiPedia:Clark_Kerr Clark Kerr, 12th President of the UC System University of California.

  Anyone have any idea whats moving into Kerr after the math department moves? By the way its a fitting move, coming to join the engineers and physicists on south campus. Users/JamesDawe
     I heard that History was moving in. But that could just be a vicious rumor. Users/EricTalevich
     According to one of my Freshman Seminar professors (Hill), it will be Political Sciences. Users/JosephBleckman

Following the retreat of the Math departments the majority of the 2nd6th floors were abandoned for a year, leading to the facilities decaying to the point where it looked as if WW3 had hit it. In the summer of 2007 it began to be restored, with the DHC offices moving from the 1st to the 5th floor. This means that the DHC labs no longer suck, but are rather now in a classy setting with a cool view. And much more comfertable temperaturez.

  The DHC labs have been bunked down to the second floor, but the DHC facilities are as roomy as ever. No more free printing, and the door lock is broken though. Users/ChiakiHirai

During its disuse, the building was a prime spot for: the cleanest bathrooms and Urban Exploration

Most recently, as of Fall 2009, the Communication department has left Everson Hall, and relocated to the third floor of Kerr. The World History Project and a few communications professors share the first floor. Half the second floor seems to be taken up by the DHC. The fourth through sixth floors are under construction.

During the renewal of the building, each floor is planned to be color coded for no particular reason:
First floor: Yellow  Moving chaos
Second floor: Brown  Linguistics
Third floor: Orange  Communications
Fourth: Olive Green  Political Science and Political Science Computer Lab
Fifth: Baby Blue  Political Science
Sixth: Brighter Yellow  Political Science TAs

At seven stories high, it is one of the tallest buildings on campus and in Yolo County.  Kerr Hall also proudly sports the antenna for KDVS on its rooftop.

If youre ever in need to fight a medieval war on campus, Kerr would be an ideal location to set up a defensive front for the following reasons:

A cool thing about Kerr Hall is the Tunnels on Campus tunnel that connects it to the basement of Wellman Hall.  There is, however, a lock on the door between them.  The trick is that this lock is oneway, and you can always enter Wellman from Kerr, if you happen to be in Kerr when Wellman is locked up.


text:
Mo Sto is a wiki user

20060906 20:15:40 It is possible that they are just as bad.  I am currently renting from them, and everything has been fine until we had an issue come up.  Our landscaping service has not been keeping up, and after 3 months of contacting Acadian Properties Acadian we still dont have a resolution.  So much of my wasted time!  And one of the girls there, April, actually hung up on me cause she didnt want to help me!  If you can possibly rent elsewhere I say do it.  But if you really like the property, at least go in with your eyes open.  At the first hint of trouble, contact Community Mediation Service at 7575623.  They are VERY helpful!  Good luck!  ShelleyParker


text:
!!!!!
! This subroutine implements the full two-stream model. For
! information on the logic and theoretical basis for the calculations,
! see Pinty et al. (2006).
!
!  Pinty, B., Lavergne, T., Dickinson, R.E., Widlowski, J.-L., Gobron,
!  N., Verstraete, M.M., 2006. Simplifying the interaction of land
!  surfaces with radiation for relating remote sensing products to
!  climate models. J. Geophys. Res. 111, D02116.
!  doi:10.1029/2005JD005952
!!!!!

subroutine two_stream(rl, tl, solar_zenith, L, soil_moisture, &     ! Inputs
        alpha_c, alpha_i, Tc, Ti, Ac, Ai)       ! Outputs

    use mod_dataspec_wavelength
    use mod_hapke_soil
    use mod_biggamma

    implicit none
    real*8, intent(in), dimension(nw) :: rl, tl
    real*8, intent(in) :: L, solar_zenith, soil_moisture

    real*8, intent(out), dimension(nw) :: alpha_c, alpha_i, Tc, Ti, Ac, Ai

    !! TODO: More intelligent treatment of zeta_c and zeta_i
    !! TODO: Scale leaf optics by forward/backward scattering efficiency (depends on leaf orientation function)
    !! TODO: Implement alternative G functions (as callable functions of mu)

    !! Fixed parameters
    real*8, parameter :: mui = 0.5d0 / 0.705d0      ! Isotropic cosine constant
    real*8, parameter :: G = 0.5d0                ! Structure factor for spherically (randomly) distributed leaf angles
    real*8, parameter :: zeta_c = 0.67d0          ! Structure factor for collimated radiation (see Pinty 2006)
    real*8, parameter :: zeta_i_factor = 1.0d0    ! Structure factor scaling term for isotropic radiation (see Pinty 2006)

    !! Derived quantities
    real*8 :: mu0                       ! Solar zenith angle (in radians)
    real*8 :: KL                        ! Vegetation optical depth (K * L, where K = G / mu0)
    real*8 :: zeta_i                    ! Structure factor for isotropic radiation
    real*8 :: tau_c, tau_i, exptaui             ! Optical depths for collimated (c) and isotropic (i) radiation
    real*8, dimension(nw, 4) :: gammas_c, gammas_i      ! gamma terms for Meador & Weaver (1980) solutions
    real*8, dimension(nw) :: w0, rho_soil               ! Canopy albedo, soil reflectance
    real*8, dimension(nw) :: RcBB, TcBB, RiBB, TiBB     ! Black background reflectance and transmittance
    real*8 :: TcBC, TiBC                                ! Black canopy transmittances
    real*8, dimension(nw) :: Rstar                      ! Infinite series reflectance term
    real*8, dimension(nw) :: TcT, TiT                   ! Total one-way transmittances
    real*8, dimension(nw) :: rhoxTiT

    mu0 = cos(solar_zenith)
    zeta_i = zeta_c * zeta_i_factor

    w0 = rl + tl
    KL = L * G / mu0
    tau_c = KL * zeta_c
    tau_i = KL * zeta_i

    !! Calculate soil reflectance from soil moisture (simple Hapke model)
    rho_soil = hapke_soil(soil_moisture)

    !! Gamma coefficients
    call pinty(rl, tl, mu0, gammas_c)
    call pinty(rl, tl, mui, gammas_i)

    !! Black background solutions (from Meador & Weaver 1980, plane-parallel atmospheres)
    call meador_plane_parallel(w0, mu0, tau_c, &
        gammas_c(:,1), gammas_c(:,2), gammas_c(:,3), gammas_c(:,4), &
        RcBB, TcBB)
    call meador_plane_parallel(w0, mui, tau_i, &
        gammas_i(:,1), gammas_i(:,2), gammas_i(:,3), gammas_i(:,4), &
        RiBB, TiBB)

    !! Black canopy transmittances
    TcBC = exp(-tau_c / mu0)
    exptaui = exp(tau_i)
    TiBC = (1.0d0 - tau_i + tau_i * tau_i * BigGamma(tau_i)) / exptaui

    !! Infinite series reflectance term
    Rstar = 1.0d0 / (1.0d0 - RiBB * rho_soil)

    !! Total one-way transmittance
    TcT = TcBB + TcBC
    TiT = TiBB + TiBC

    !! Results: Transmitance to background layer
    Tc = TcT * Rstar
    Ti = TiT * Rstar

    !! Results: Canopy reflectance
    !! Full formulae for reference
    !! alpha_c = RcBB + TcT * rho_soil * TiT * Rstar
    !! alpha_i = RiBB + TiT * rho_soil * TiT * Rstar
    rhoxTiT = rho_soil * TiT
    alpha_c = RcBB + Tc * rhoxTiT
    alpha_i = RiBB + Ti * rhoxTit

    !! Results: Canopy absorption
    !! Full formulae for reference
    !! Ac = 1 - (alpha_c + Tc) + rho_soil * TiT
    !! Ai = 1 - (alpha_i + Ti) + rho_soil * TiT
    Ac = 1 - (alpha_c + Tc) + rhoxTiT
    Ai = 1 - (alpha_i + Ti) + rhoxTiT

end subroutine two_stream

max_stars_repo_name: ankurdesai/pecan

max_stars_repo_licenses: NCSA Unlicense

max_issues_repo_name: ankurdesai/pecan

max_issues_repo_licenses: NCSA Unlicense

max_forks_repo_name: ankurdesai/pecan

max_forks_repo_licenses: NCSA Unlicense


text:
!
! SPDX-License-Identifier: MIT
!
module mod_funcs
  !
  use mod_param, only: pi
  use mod_types, only: rp
  !
  implicit none
  !
  private
  public ssign,heaviside,dirac,efun,interp_d,interp_g
  !
  contains
  !
  function ssign(delta,phi)
    !
    ! smooth sign function
    !
    implicit none
    !
    real(rp), intent(in) :: delta,phi
    real(rp)             :: ssign
    !
    ssign = sign(1._rp,phi)
    if(abs(phi).le.delta) then
      ssign = phi/(sqrt(phi**2+delta**2))
    endif
    !
    return
  end function ssign
  !
  function heaviside(r,eps)
    !
    ! smooth step function based
    ! on cosine/sine
    !
    implicit none
    !
    real(rp), intent(in) :: r,eps
    real(rp)             :: heaviside
    !
    if(r.lt.-eps) then
      heaviside = 0._rp
    elseif(r.lt.eps) then
      heaviside = 0.5_rp + 0.5_rp*r/eps + 0.5_rp/pi*sin(pi*r/eps)
    else
      heaviside = 1._rp
    endif
    !
    return
  end function heaviside
  !
  function dirac(r,eps)
    !
    ! smooth impulse function
    ! (note: the derivative of heaviside function w.r.t. eps)
    !
    implicit none
    !
    real(rp), intent(in) :: r,eps
    real(rp)             :: dirac
    !
    if(abs(r).ge.eps) then
      dirac = 0._rp
    else
      dirac = 0.5_rp/eps + 0.5_rp/eps*cos(pi*r/eps)
    endif
    !
    return
  end function dirac
  !
  function efun(r,eps)
    !
    ! smooth step function based
    ! on the error function
    !
    implicit none
    !
    real(rp), intent(in) :: r,eps
    real(rp)             :: efun
    !
    efun = 0.5_rp*( 1._rp+erf(r/eps) )
    !
    return
  end function efun
  !
  function dir_efun(r,eps)
    !
    ! smooth impulse function based
    ! on the error function
    ! (note: the derivative of error function w.r.t. eps)
    !
    implicit none
    !
    real(rp), intent(in) :: r,eps
    real(rp)             :: dir_efun
    !
    dir_efun = (1._rp/(eps*pi))*exp(-(r/eps)**2._rp)
    !
    return
  end function dir_efun
  !
  function interp_d(vec,vel,q) result(flux) ! divergence form
    !
    ! QUICK discretization (by Leonard) in divergence form
    ! (see, e.g.: Numerical simulation of incompressible flows, pag. 112)
    !
    implicit none
    !
    real(rp), intent(in), dimension(-2:2) :: vec
    real(rp), intent(in) :: vel ! the transporting velocity (in d(uu)/dx is u, in d(uv)/dy is v)
    integer , intent(in) :: q   ! (q=-1 if f_i, q=0 if f_i+1 with f is the flux)
    !
    real(rp) :: flux
    !
    !flux = 0.5_rp*vel*(vec(0+q)+vec(1+q)) + 0.5_rp*abs(vel)*(+vec(0+q)-vec(1+q))
    !
    flux =      0.0625_rp*(vel)*(-vec(-1+q)+9.0_rp*vec(0+q)+9.0_rp*vec(1+q)-vec(2+q)) &
           + 0.0625_rp*abs(vel)*(-vec(-1+q)+3.0_rp*vec(0+q)-3.0_rp*vec(1+q)+vec(2+q))
    !
    return
  end function interp_d
  !
  function interp_g(vec,vel,dli) result(vel_grad) ! gradient form
    !
    ! II order (Kawamura and Kuwahara equation) discretization in gradient form
    ! (see, e.g.: Numerical simulation of incompressible flows, pag. 113)
    !
    implicit none
    !
    real(rp), intent(in), dimension(-2:2) :: vec
    real(rp), intent(in)                  :: vel ! the transporting velocity (in ud(f)/dx is u, in vd(f)/dy is v)
    real(rp), intent(in)                  :: dli
    !
    real(rp) :: vel_grad
    !
    !vel_grad = 0.5_rp*(vel+abs(vel))*(vec(0)-vec(-1))*dli + 0.5_rp*(vel-abs(vel))*(vec(1)-vec(0))*dli
    !
    vel_grad =             (vel)*(vec(-2)-8.0_rp*vec(-1)               +8.0_rp*vec(+1)-vec(+2))*(1.0_rp/12.0_rp)*dli &
               + 3.0_rp*abs(vel)*(vec(-2)-4.0_rp*vec(-1)+6.0_rp*vec(+0)-4.0_rp*vec(+1)+vec(+2))*(1.0_rp/12.0_rp)*dli
    !
    return
  end function interp_g
  !
end module mod_funcs

max_stars_repo_name: fspiga/FluTAS

max_issues_repo_name: fspiga/FluTAS

max_forks_repo_name: fspiga/FluTAS


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine mechmodel(elconloc,elas,emec,kode,emec0,ithermal,
     &     icmd,beta,stre,xkl,ckl,vj,xikl,vij,plconloc,xstate,xstateini,
     &     ielas,amat,t1l,dtime,time,ttime,iel,iint,nstate_,mi,
     &     iorien,pgauss,orab,eloc,mattyp,pnewdt,istep,iinc,ipkon,
     &     nmethod,iperturb)
!
!     kode=-1: Arruda-Boyce
!          -2: Mooney-Rivlin
!          -3: Neo-Hooke
!          -4: Ogden (N=1)
!          -5: Ogden (N=2)
!          -6: Ogden (N=3)
!          -7: Polynomial (N=1)
!          -8: Polynomial (N=2)
!          -9: Polynomial (N=3)
!          -10: Reduced Polynomial (N=1)
!          -11: Reduced Polynomial (N=2)
!          -12: Reduced Polynomial (N=3)
!          -13: Van der Waals (not implemented yet)
!          -14: Yeoh
!          -15: Hyperfoam (N=1)
!          -16: Hyperfoam (N=2)
!          -17: Hyperfoam (N=3)
!          -50: deformation plasticity
!          -51: incremental plasticity (no viscosity)
!          -52: viscoplasticity
!       < -100: user material routine with -kode-100 user
!               defined constants with keyword *USER MATERIAL
!
      implicit none
!
      character*80 amat
!
      integer kode,ithermal,icmd,ielas,iel,iint,nstate_,mi(*),iorien,
     &  mattyp,istep,iinc,ipkon(*),nmethod,iperturb(*)
!
      real*8 elconloc(*),elas(21),emec(*),emec0(*),beta(*),stre(*),
     &  ckl(*),vj,plconloc(*),t1l,xkl(*),xikl(*),vij,
     &  dtime,didc(27),d2idc2(243),dibdc(27),d2ibdc2(243),
     &  dudc(9),d2udc2(81),dldc(27),d2ldc2(243),dlbdc(27),d2lbdc2(243),
     &  pgauss(3),orab(7,*),time,ttime,eloc(6),pnewdt
!
      real*8 xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*)
!
      if(kode.gt.0) then
         call linel(kode,mattyp,beta,emec,stre,elas,elconloc,
     &  iorien,orab,pgauss)
      elseif(kode.gt.-50) then
         mattyp=3
         call rubber(elconloc,elas,emec,kode,didc,d2idc2,
     &     dibdc,d2ibdc2,dudc,d2udc2,dldc,d2ldc2,dlbdc,d2lbdc2,
     &     ithermal,icmd,beta,stre)
      elseif(kode.eq.-50) then
         mattyp=3
         call defplas(elconloc,elas,emec,ithermal,icmd,beta,stre,
     &     ckl,vj)
      elseif(kode.gt.-100) then
         mattyp=3
         if(iperturb(2).eq.1) then
            call incplas(elconloc,plconloc,xstate,xstateini,elas,emec,
     &           ithermal,icmd,beta,stre,vj,kode,ielas,amat,t1l,dtime,
     &           time,ttime,iel,iint,nstate_,mi(1),eloc,pgauss,nmethod,
     &           pnewdt)
         else
            call incplas_lin(elconloc,plconloc,xstate,xstateini,elas,
     &           emec,
     &           ithermal,icmd,beta,stre,vj,kode,ielas,amat,t1l,dtime,
     &           time,ttime,iel,iint,nstate_,mi(1),eloc,pgauss,nmethod,
     &           pnewdt)
         endif
      else
         mattyp=3
         call umat_main(amat,iel,iint,kode,elconloc,emec,emec0,beta,
     &        xikl,vij,xkl,vj,ithermal,t1l,dtime,time,ttime,icmd,ielas,
     &        mi(1),nstate_,xstateini,xstate,stre,elas,iorien,pgauss,
     &        orab,pnewdt,istep,iinc,ipkon,nmethod,iperturb)
      endif
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
      SUBROUTINE PHO_DSIGDT(Ee,Xta,Nfill)
C*********************************************************************
C
C     calculation of unitarized amplitude
C                    and differential cross section
C
C     input:   EE       cm energy (GeV)
C              XTA(1,*) t values (GeV**2)
C              NFILL    entries in t table
C
C     output:  XTA(2,*)  DSIG/DT  g p --> g h/V (mub/GeV**2)
C              XTA(3,*)  DSIG/DT  g p --> rho0 h/V
C              XTA(4,*)  DSIG/DT  g p --> omega0 h/V
C              XTA(5,*)  DSIG/DT  g p --> phi h/V
C              XTA(6,*)  DSIG/DT  g p --> pi+ pi- h/V (continuum)
C
C*********************************************************************
      IMPLICIT NONE
      DOUBLE PRECISION bmax , DEPS , Ee , etmp , fac , PHO_BESSJ0 ,
     &                 THOUS , wg , wght , xpnt , Xta , xx
      INTEGER i , i1 , i2 , iamp , ITHREE , ITWO , j , j1 , j2 , k ,
     &        k1 , k2 , l , l1 , l2 , Nfill
      SAVE

      PARAMETER (ITWO=2,ITHREE=3,THOUS=1.D3,DEPS=1.D-20)

      DIMENSION Xta(6,Nfill)

C  input/output channels
      INCLUDE 'inc/poinou'
C  some constants
      INCLUDE 'inc/pocons'
C  integration precision for hard cross sections (obsolete)
      INCLUDE 'inc/pogaup'
C  event debugging information
      INCLUDE 'inc/podebg'
C  global event kinematics and particle IDs
      INCLUDE 'inc/pogcms'
C  complex Born graph amplitudes used for unitarization
      INCLUDE 'inc/point4'

      COMPLEX*16 xt , amp , czero
      DIMENSION amp(5) , xpnt(96) , wght(96) , xt(5,100)
      CHARACTER*12 fna

      czero = DCMPLX(0.D0,0.D0)

      etmp = ECM
      ECM = Ee

      IF ( Nfill.GT.100 ) THEN
         IF ( LPRi.GT.4 ) WRITE (LO,'(1X,A,I4)')
     &         'PHO_DSIGDT:ERROR: too many entries in table' , Nfill
         STOP
      END IF
C
      DO k = 1 , Nfill
         DO l = 1 , 5
            xt(l,k) = czero
         END DO
      END DO
C
C  impact parameter integration
C     BMAX=12.D0*SQRT(MAX(BPOM,BREG))
      bmax = 10.D0
      CALL PHO_GAUSET(0.D0,bmax,NGAuso,xpnt,wght)
      iamp = 5
      IF ( (IFPap(1).EQ.22) .AND. (IFPap(2).NE.22) ) THEN
         i1 = 1
         i2 = 0
      ELSE IF ( (IFPap(1).NE.22) .AND. (IFPap(2).EQ.22) ) THEN
         i1 = 0
         i2 = 1
      ELSE IF ( (IFPap(1).EQ.22) .AND. (IFPap(2).EQ.22) ) THEN
         i1 = 1
         i2 = 1
      ELSE
         i1 = 0
         i2 = 0
         iamp = 1
      END IF
      j1 = i1*2
      k1 = i1*3
      l1 = i1*4
      j2 = i2*2
      k2 = i2*3
      l2 = i2*4
C
      DO i = 1 , NGAuso
         wg = wght(i)*xpnt(i)
C  calculate amplitudes
         IF ( i.EQ.1 ) THEN
            CALL PHO_EIKON(1,-1,xpnt(i))
         ELSE
            CALL PHO_EIKON(1,1,xpnt(i))
         END IF
         amp(1) = AMPel
         amp(2) = AMPvm(i1,i2)
         amp(3) = AMPvm(j1,j2)
         amp(4) = AMPvm(k1,k2)
         amp(5) = AMPvm(l1,l2)
C
         DO j = 1 , Nfill
            xx = xpnt(i)*SQRT(Xta(1,j)/GEV2mb)
            fac = PHO_BESSJ0(xx)*wg
            DO k = 1 , iamp
               xt(1,j) = xt(1,j) + amp(k)*fac
            END DO
         END DO
      END DO
C
C  change units to mb/GeV**2
      fac = 4.D0*PI/GEV2mb
      fna = '(mb/GeV**2) '
      IF ( i1+i2.EQ.1 ) THEN
         fac = fac*THOUS
         fna = '(mub/GeV**2)'
      ELSE IF ( i1+i2.EQ.2 ) THEN
         fac = fac*THOUS*THOUS
         fna = '(nb/GeV**2) '
      END IF
      IF ( IDEb(56).GE.5 ) THEN
         IF ( LPRi.GT.4 ) WRITE (LO,'(1X,A,A12,/1X,A)')
     &         'table:  -T (GeV**2)   DSIG/DT ' , fna ,
     &        '------------------------------------------'
      END IF
      DO j = 1 , Nfill
         DO k = 1 , iamp
            Xta(k+1,j) = ABS(xt(k,j))**2*fac
         END DO
         IF ( IDEb(56).GE.5 ) THEN
            IF ( LPRi.GT.4 ) WRITE (LO,'(5X,6E12.3)')
     &           (Xta(i,j),i=1,iamp+1)
         END IF
      END DO

      ECM = etmp
      END SUBROUTINE

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
SUBROUTINE dple(rowk, n, b, c, ierr)

! Code converted using TO_F90 by Alan Miller
! Date: 2003-06-16  Time: 12:26:32

!     ******************************************************************
!     SOLUTION OF LINEAR EQUATIONS WITH REDUCED STORAGE
!     ******************************************************************

! Uses the Henderson-Wassyng partial pivot algorithm.
! Wassyng, A. 'Solving Ax = b: A method with reduced storage requirements',
! SIAM J. Numerical Analysis, vol.19 (1982), pp. 197-204.

! The user must provide a routine ROWK to return the requested row of the
! matrix A.

! N.B. Arguments D and IP have been removed.

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(14, 60)

INTEGER, INTENT(IN)     :: n
REAL (dp), INTENT(IN)   :: b(n)
REAL (dp), INTENT(OUT)  :: c(n)
INTEGER, INTENT(OUT)    :: ierr

! EXTERNAL rowk
INTERFACE
  SUBROUTINE rowk(n, k, r)
    IMPLICIT NONE
    INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(14, 60)
    INTEGER, INTENT(IN)     :: n, k
    REAL (dp), INTENT(OUT)  :: r(:)
  END SUBROUTINE rowk
END INTERFACE

! Local variables

REAL (dp)  :: bk, cj, ck, c1, dkj
REAL (dp), PARAMETER  :: zero = 0.0_dp
REAL (dp)  :: wk(n*n/4 + n + 3)
INTEGER    :: i, iflag, ij, ijold, ik, iwk(n), j, k, kjold, km1, kp1,   &
              last, lastm1, lcol, lcolp1, m, maxwk, mjold, nm1, np1

!     SET THE NECESSARY CONSTANTS

ierr = 0
maxwk = n * n / 4 + n + 3
np1 = n + 1
k = 1
iflag = -1

!     GET THE FIRST COLUMN OF THE TRANSPOSED SYSTEM

CALL rowk(n, 1, c)
bk = b(1)

IF (n <= 1) THEN
  IF (c(1) == zero) GO TO 130
  c(1) = bk / c(1)
  RETURN
END IF

!     FIND THE PIVOT FOR COLUMN 1

m = 1
DO  i = 2, n
  IF (ABS(c(m)) < ABS(c(i))) m = i
END DO

iwk(1) = m
c1 = c(m)
c(m) = c(1)
c(1) = c1
IF (c(1) /= zero) THEN

!     FIND THE FIRST ELEMENTARY MATRIX AND STORE IT IN D

  DO  i = 2, n
    wk(i-1) = -c(i) / c(1)
  END DO
  wk(n) = bk / c(1)

!     K LOOP - EACH K FOR A NEW COLUMN OF THE TRANSPOSED SYSTEM

  DO  k = 2, n
    kp1 = k + 1
    km1 = k - 1

!       GET COLUMN K

    CALL rowk(n, k, c)
    DO  j = 1, km1
      m = iwk(j)
      cj = c(j)
      c(j) = c(m)
      c(m) = cj
    END DO
    bk = b(k)

    iflag = -iflag
    lcol = np1 - k
    lcolp1 = lcol + 1
    lastm1 = 1
    last = maxwk - n + k
    IF (k /= 2) THEN

      lastm1 = maxwk - n + km1
      IF (iflag < 0) last = last - n + k - 2
      IF (iflag > 0) lastm1 = lastm1 - n + k - 3
    END IF

!     J LOOP - EFFECT OF COLUMNS 1 TO K-1 OF L-INVERSE

    DO  j = 1, km1
      cj = c(j)
      ij = (j-1) * lcolp1
      IF (j == km1) ij = lastm1 - 1

!     I LOOP - EFFECT OF L-INVERSE ON ROWS K TO N+1

      DO  i = k, n
        ij = ij + 1
        c(i) = c(i) + wk(ij) * cj
      END DO
      bk = bk - wk(ij+1) * cj
    END DO

!       K=N CASE

    m = k
    IF (k >= n) THEN
      IF (c(k) == zero) GO TO 130
      wk(last) = bk / c(k)
    ELSE

!     FIND THE PIVOT

      DO  i = kp1, n
        IF (ABS(c(m)) < ABS(c(i))) m = i
      END DO

      iwk(k) = m
      ck = c(m)
      c(m) = c(k)
      c(k) = ck
      IF (c(k) == zero) GO TO 130

!     FIND THE K-TH ELEMENTARY MATRIX

      ik = last
      DO  i = kp1, n
        wk(ik) = -c(i) / c(k)
        ik = ik + 1
      END DO
      wk(ik) = bk / c(k)
    END IF

!     FORM THE PRODUCT OF THE ELEMENTARY MATRICES

    DO  j = 1, km1
      kjold = j * lcolp1 + k - np1
      mjold = kjold + m - k
      ij = (j-1) * lcol
      ijold = ij + j
      IF (j == km1) THEN

        kjold = lastm1
        mjold = lastm1 + m - k
        ijold = lastm1
      END IF

      ik = last - 1
      dkj = wk(mjold)
      wk(mjold) = wk(kjold)
      DO  i = kp1, np1
        ij = ij + 1
        ijold = ijold + 1
        ik = ik + 1
        wk(ij) = wk(ijold) + wk(ik) * dkj
      END DO
    END DO
  END DO

  last = maxwk
  IF (iflag < 0) last = maxwk - 2
  wk(n) = wk(last)

!     INSERT THE SOLUTION IN C

  c(1:n) = wk(1:n)

  nm1 = n - 1
  DO  i = 1, nm1
    k = n - i
    m = iwk(k)
    ck = c(k)
    c(k) = c(m)
    c(m) = ck
  END DO
  RETURN
END IF

!     THE SYSTEM IS SINGULAR

130 ierr = k
RETURN
END SUBROUTINE dple

SUBROUTINE rowk(n, k, r)

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(14, 60)
INTEGER, INTENT(IN)     :: n, k
REAL (dp), INTENT(OUT)  :: r(:)

INTEGER    :: i
REAL (dp)  :: temp

r(k) = 1.0_dp
IF (k > 1) THEN
  temp = 1.0_dp
  DO i = 1, k-1
    temp = 0.99_dp * temp
    r(k-i) = temp
  END DO
END IF

IF (k < n) THEN
  temp = 1.0_dp
  DO i = 1, n-k
    temp = 0.99_dp * temp
    r(k+i) = temp
  END DO
END IF

RETURN
END SUBROUTINE rowk

PROGRAM Test_DPLE
IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(14, 60)

REAL (dp)  :: a(100), x(100), b(100), error, temp, soln(100)
INTEGER    :: i, ierr, row

INTERFACE
  SUBROUTINE dple(rowk, n, b, c, ierr)
    IMPLICIT NONE
    INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(14, 60)
    INTEGER, INTENT(IN)     :: n
    REAL (dp), INTENT(IN)   :: b(n)
    REAL (dp), INTENT(OUT)  :: c(n)
    INTEGER, INTENT(OUT)    :: ierr
    INTERFACE
      SUBROUTINE rowk(n, k, r)
        IMPLICIT NONE
        INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(14, 60)
        INTEGER, INTENT(IN)     :: n, k
        REAL (dp), INTENT(OUT)  :: r(:)
      END SUBROUTINE rowk
    END INTERFACE
  END SUBROUTINE dple

  SUBROUTINE rowk(n, k, r)
    IMPLICIT NONE
    INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(14, 60)
    INTEGER, INTENT(IN)     :: n, k
    REAL (dp), INTENT(OUT)  :: r(:)
  END SUBROUTINE rowk
END INTERFACE

! Generate the solution using the random number generator
CALL RANDOM_NUMBER(x)
! Calculate the right-hand side
DO row = 1, 100
  CALL rowk(100, row, a)
  b(row) = DOT_PRODUCT(a, x)
END DO

CALL dple(rowk, 100, b, soln, ierr)
IF (ierr == 0) THEN
  temp = 0.0
  DO i = 1, 100
    error = ABS(soln(i) - x(i))
    temp = MAX(temp, error)
  END DO
  WRITE(*, '(a, g13.5)') ' Max. error = ', temp
ELSE
  WRITE(*, *) ' Error = ', ierr
END IF

STOP
END PROGRAM Test_DPLE

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


subroutine ensctl2state_ad(eval,mval,grad)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    ensctl2state_ad
!   prgmmr: kleist
!
! abstract:  Contribution from state space to ensemble control vector
!
! program history log:
!   2011-11-17  kleist - initial code
!   2013-10-28  todling - rename p3d to prse
!   2013-11-22  kleist - add option for q perturbations
!   2014-12-03  derber   - introduce parallel regions for optimization
!   2017-05-12  Y. Wang and X. Wang - add w as state variable for rw DA, POC: xuguang.wang@ou.edu
!   2019-07-11  Todling - there should be no need to check on the existence of w and dw
!
!   input argument list:
!     eval - Ensemble state variable variable
!     grad - Control variable
!
!   output argument list:
!     grad - Control variable
!
!$$$ end documentation block

use kinds, only: r_kind,i_kind
use control_vectors, only: control_vector,cvars3d
use gsi_4dvar, only: ibin_anl
use hybrid_ensemble_parameters, only: uv_hyb_ens,dual_res,ntlevs_ens,q_hyb_ens
use hybrid_ensemble_isotropic, only: ensemble_forward_model_ad
use hybrid_ensemble_isotropic, only: ensemble_forward_model_ad_dual_res
use balmod, only: strong_bk_ad
use gsi_bundlemod, only: gsi_bundlecreate
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
use gsi_bundlemod, only: gsi_bundlegetvar
use gsi_bundlemod, only: gsi_bundleputvar
use gsi_bundlemod, only: gsi_bundledestroy
use gsi_bundlemod, only: assignment(=)
use gsi_bundlemod, only : self_add
use constants, only: zero,max_varname_length
use mpeu_util, only: getindex
use gsi_metguess_mod, only: gsi_metguess_get
use mod_strong, only: tlnmc_option
use cwhydromod, only: cw2hydro_ad
use cwhydromod, only: cw2hydro_ad_hwrf
use timermod, only: timer_ini,timer_fnl
use gridmod, only: nems_nmmb_regional
implicit none

! Declare passed variables
type(control_vector), intent(inout) :: grad
type(gsi_bundle)    , intent(inout) :: mval
type(gsi_bundle)    , intent(in   ) :: eval(ntlevs_ens)

! Declare local variables
character(len=*),parameter::myname='ensctl2state_ad'
character(len=max_varname_length),allocatable,dimension(:) :: clouds
integer(i_kind) :: jj,ic,id,istatus,nclouds

integer(i_kind), parameter :: ncvars = 8
integer(i_kind) :: icps(ncvars)
type(gsi_bundle):: wbundle_c ! work bundle
character(len=3), parameter :: mycvars(ncvars) = (/  &  ! vars from CV needed here
                               'sf ', 'vp ', 'ps ', 't  ',    &
                               'q  ', 'cw ', 'w  ', 'dw '/)
logical :: lc_sf,lc_vp,lc_ps,lc_t,lc_rh,lc_cw
logical :: lc_w,lc_dw
real(r_kind),pointer,dimension(:,:,:) :: cv_sf,cv_vp,cv_rh
! Declare required local state variables
integer(i_kind), parameter :: nsvars = 13
integer(i_kind) :: isps(nsvars)
character(len=4), parameter :: mysvars(nsvars) = (/  &  ! vars from ST needed here
                               'u   ', 'v   ', 'prse', 'q   ', 'tsen','ql  ','qi  ', &
                               'qr  ', 'qs  ', 'qg  ', 'qh  ', 'w   ','dw  ' /)
logical :: ls_u,ls_v,ls_prse,ls_q,ls_tsen,ls_ql,ls_qi
logical :: ls_qr,ls_qs,ls_qg,ls_qh
logical :: ls_w,ls_dw
real(r_kind),pointer,dimension(:,:)   :: rv_ps,rv_sst
real(r_kind),pointer,dimension(:,:,:) :: rv_u,rv_v,rv_prse,rv_q,rv_tsen,rv_tv,rv_oz
real(r_kind),pointer,dimension(:,:,:) :: rv_rank3,rv_w,rv_dw

logical :: do_getuv,do_tv_to_tsen_ad,do_normal_rh_to_q_ad,do_getprs_ad,lstrong_bk_vars
logical :: do_tlnmc,do_q_copy
logical :: do_cw_to_hydro_ad
logical :: do_cw_to_hydro_ad_hwrf
logical :: wdw_exist

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

! Initialize timer
!call timer_ini(trim(myname))

! Inquire about chemistry
call gsi_metguess_get('clouds::3d',nclouds,istatus)
if (nclouds>0) then
    allocate(clouds(nclouds))
    call gsi_metguess_get('clouds::3d',clouds,istatus)
endif

! Since each internal vector of grad has the same structure, pointers are
! the same independent of the subwindow jj
call gsi_bundlegetpointer (grad%step(1),mycvars,icps,istatus)
lc_sf =icps(1)>0; lc_vp =icps(2)>0; lc_ps =icps(3)>0
lc_t  =icps(4)>0; lc_rh =icps(5)>0; lc_cw =icps(6)>0
lc_w  =icps(7)>0; lc_dw =icps(8)>0

! Since each internal vector of grad has the same structure, pointers are
! the same independent of the subwindow jj
call gsi_bundlegetpointer (eval(1),mysvars,isps,istatus)
ls_u  =isps(1)>0; ls_v   =isps(2)>0; ls_prse=isps(3)>0
ls_q  =isps(4)>0; ls_tsen=isps(5)>0; ls_ql =isps(6)>0; ls_qi =isps(7)>0
ls_qr  =isps(8)>0; ls_qs  =isps(9)>0
ls_qg  =isps(10)>0; ls_qh =isps(11)>0
ls_w   =isps(12)>0; ls_dw =isps(13)>0

! Define what to do depending on what's in CV and SV
lstrong_bk_vars     =lc_sf.and.lc_vp.and.lc_ps .and.lc_t
do_getuv            =lc_sf.and.lc_vp.and.ls_u  .and.ls_v
do_tv_to_tsen_ad    =lc_t .and.ls_q .and.ls_tsen
do_normal_rh_to_q_ad=(.not.q_hyb_ens).and.&
                     lc_t .and.lc_rh.and.ls_prse.and.ls_q
do_q_copy=.false.
if(.not. do_normal_rh_to_q_ad) then
  do_q_copy = lc_rh.and.lc_t .and.ls_prse.and.ls_q.and.q_hyb_ens
end if
do_getprs_ad        =lc_t .and.lc_ps.and.ls_prse

do_cw_to_hydro_ad=.false.
do_cw_to_hydro_ad=lc_cw.and.ls_ql.and.ls_qi
do_cw_to_hydro_ad_hwrf=.false.
do_cw_to_hydro_ad_hwrf=lc_cw.and.ls_ql.and.ls_qi.and.ls_qr.and.ls_qs.and.ls_qg.and.ls_qh

wdw_exist = lc_w.and.lc_dw.and.ls_w.and.ls_dw

! Initialize
mval%values=zero
!  Create a temporary bundle similar to grad, and copy contents of grad into it
call gsi_bundlecreate ( wbundle_c, grad%step(1), 'ensctl2state_ad work', istatus )
if(istatus/=0) then
   write(6,*) trim(myname), ': trouble creating work bundle'
   call stop2(999)
endif

do jj=1,ntlevs_ens

! If calling TLNMC, already have u,v (so set last argument to true)
   do_tlnmc = lstrong_bk_vars .and. ( (tlnmc_option==3) .or. &
            (jj==ibin_anl .and. tlnmc_option==2))

   wbundle_c%values=zero

! Get sv pointers here
!  Get pointers to required state variables
   call gsi_bundlegetpointer (eval(jj),'u'   ,rv_u,   istatus)
   call gsi_bundlegetpointer (eval(jj),'v'   ,rv_v,   istatus)
   call gsi_bundlegetpointer (eval(jj),'ps'  ,rv_ps,  istatus)
   call gsi_bundlegetpointer (eval(jj),'prse',rv_prse,istatus)
   call gsi_bundlegetpointer (eval(jj),'tv'  ,rv_tv,  istatus)
   call gsi_bundlegetpointer (eval(jj),'tsen',rv_tsen,istatus)
   call gsi_bundlegetpointer (eval(jj),'q'   ,rv_q ,  istatus)
   call gsi_bundlegetpointer (wbundle_c,'q'  ,cv_rh ,istatus)

!  Adjoint of consistency for sensible temperature, calculate sensible temperature
   if(do_tv_to_tsen_ad) call tv_to_tsen_ad(rv_tv,rv_q,rv_tsen)

   if(do_tlnmc) then

      ! Adjoint to convert ps to 3-d pressure
      if(do_getprs_ad) call getprs_ad(rv_ps,rv_tv,rv_prse)
      rv_prse=zero

      ! Adjoint of strong_bk
      call strong_bk_ad(rv_u,rv_v,rv_ps,rv_tv,.true.)

   end if

   call self_add(mval,eval(jj))

!$omp parallel sections private(ic,id,istatus)

!$omp section

!  Convert RHS calculations for u,v to st/vp
   if (do_getuv) then
      if(uv_hyb_ens) then
         call gsi_bundleputvar ( wbundle_c, 'sf', rv_u, istatus )
         call gsi_bundleputvar ( wbundle_c, 'vp', rv_v, istatus )
      else
         call gsi_bundlegetpointer (wbundle_c,'sf' ,cv_sf ,istatus)
         call gsi_bundlegetpointer (wbundle_c,'vp' ,cv_vp ,istatus)
         call getuv(rv_u,rv_v,cv_sf,cv_vp,1)
      end if
   end if

!$omp section

   call gsi_bundlegetpointer (eval(jj),'oz'  ,rv_oz , istatus)
   call gsi_bundlegetpointer (eval(jj),'sst' ,rv_sst, istatus)
   call gsi_bundleputvar ( wbundle_c, 'oz',  rv_oz,  istatus )
   call gsi_bundleputvar ( wbundle_c, 'sst', rv_sst, istatus )
   if(wdw_exist)then
     call gsi_bundlegetpointer (eval(jj),'w' ,rv_w, istatus)
     call gsi_bundleputvar ( wbundle_c, 'w', rv_w, istatus )
     if(nems_nmmb_regional)then
       call gsi_bundlegetpointer (eval(jj),'dw' ,rv_dw, istatus)
       call gsi_bundleputvar ( wbundle_c, 'dw', rv_dw, istatus )
     end if
   end if

!$omp section

   if (do_cw_to_hydro_ad .and. .not.do_cw_to_hydro_ad_hwrf) then
!     Case when cloud-vars do not map one-to-one
!     e.g. cw-to-ql&qi
      call cw2hydro_ad(eval(jj),wbundle_c,clouds,nclouds)
   elseif (do_cw_to_hydro_ad_hwrf) then
!!     Case when cloud-vars do not map one-to-one
!!     e.g. cw-to-ql&qi&qr&qs&qg&qh
      call cw2hydro_ad_hwrf(eval(jj),wbundle_c,rv_tsen)
   else
!  Since cloud-vars map one-to-one, take care of them together
      do ic=1,nclouds
         id=getindex(cvars3d,clouds(ic))
         if (id>0) then
            call gsi_bundlegetpointer (eval(jj),clouds(ic),rv_rank3,istatus)
            call gsi_bundleputvar     (wbundle_c, clouds(ic),rv_rank3,istatus)
         endif
      enddo
   endif

!  Calculate sensible temperature
   if(do_q_copy) then
      call gsi_bundleputvar (wbundle_c, 'q', rv_q, istatus )
   else

!     Adjoint of convert input normalized RH to q to add contribution of moisture
!     to t, p , and normalized rh
      if(do_normal_rh_to_q_ad) call normal_rh_to_q_ad(cv_rh,rv_tv,rv_prse,rv_q)

!     Adjoint to convert ps to 3-d pressure
      if(do_getprs_ad) call getprs_ad(rv_ps,rv_tv,rv_prse)
   end if

!  Adjoint of control to initial state
   call gsi_bundleputvar ( wbundle_c, 't' ,  rv_tv,  istatus )
   call gsi_bundleputvar ( wbundle_c, 'ps',  rv_ps,  istatus )
!  call gsi_bundleputvar ( wbundle_c, 'q' ,  zero,   istatus )
!$omp end parallel sections

   if(dual_res) then
      call ensemble_forward_model_ad_dual_res(wbundle_c,grad%aens(1,:),jj)
   else
      call ensemble_forward_model_ad(wbundle_c,grad%aens(1,:),jj)
   end if

end do

call gsi_bundledestroy(wbundle_c,istatus)
if (istatus/=0) then
   write(6,*) trim(myname),': trouble destroying work bundle'
   call stop2(999)
endif

if (nclouds>0) deallocate(clouds)

! Finalize timer
!call timer_fnl(trim(myname))

return
end subroutine ensctl2state_ad


text:
! RUN: %flang -fopenmp -Hy,69,0x1000 -S -emit-llvm %s -o - | FileCheck %s
! RUN: %flang -i8 -fopenmp -Hy,69,0x1000 -S -emit-llvm %s -o - | FileCheck %s

subroutine reduce()
integer :: j
logical(kind=8) :: error_status = .FALSE.
!$omp parallel do reduction(.or.: error_status)
do j=1,100
end do
!$omp end parallel do
end subroutine
! //CHECK: atomicrmw

subroutine atomic_logical_1(n,val)
logical(kind=1) :: lg = .FALSE.
integer :: n, val, i
!$omp parallel
do i=1,n
!$omp atomic
lg = lg .or. val==n
end do
!$omp end parallel
end subroutine
! //CHECK: atomicrmw

subroutine atomic_logical_8(n,val)
logical(kind=8) :: lg = .FALSE.
integer :: n, val, i
!$omp parallel
do i=1,n
!$omp atomic
lg = lg .or. val==n
end do
!$omp end parallel
end subroutine
! //CHECK: atomicrmw

subroutine atomic_integer_1(n,val)
integer(kind=1) :: lg, val
integer :: n, i
!$omp parallel
do i=1,n
!$omp atomic
lg = lg + val
end do
!$omp end parallel
end subroutine
! //CHECK: atomicrmw

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
      SUBROUTINE Q4GMGS (MID,FACTOR,G)
C     &    ENTRY Q4GMGD (MID,FACTOD,D)
C
C
C     MATERIAL PROPERTY MATRIX GENERATOR ROUTINE FOR QUAD4 ELEMENT
C
C     THIS ROUTINE BUILDS THE MATERIAL PROPERTY MATRIX, G, USING THE
C     OUTPUT OF SUBROUTINE 'MAT' (/MATOUT/).
C
C     ALL THE MATERIAL OPTIONS, ISOTROPIC, ORTHOTROPIC, AND ANISOTROPIC,
C     ARE AVAILABLE.
C
C     OUTPUT WILL BE G(9) OR D(9) FOR MID1, MID2 AND MID4.
C     FOR MID3,  G(4) OR D(4) IS SENT BACK.
C
      REAL             FACTOR,G(9),CONST,MTYPE,NU12,NU21
      DOUBLE PRECISION FACTOD,D(9),DONST
      COMMON  /MATOUT/ RMTOUT(25)
      EQUIVALENCE      (RMTOUT(1),E1),(RMTOUT(2),NU12),(RMTOUT(3),E2)
C
C     SINGLE PRECISION SECTION -
C
      DO 10 I=1,9
   10 G(I) = 0.0
      MTYPE= RMTOUT(25)
      MTYP = IFIX(MTYPE+.05) - 2
      IF (MTYP) 20,30,80
C
C     ISOTROPIC MATERIALS (MAT1)
C
   20 IF (MID .NE. 3) GO TO 40
      G(1) = RMTOUT(6)
      G(4) = G(1)
      GO TO 100
C
C     ANISOTROPIC MATERIALS (MAT2)
C
   30 IF (MID .EQ. 3) GO TO 60
   40 DO 50 I=1,3
   50 G(I) = RMTOUT(I)
      G(4) = G(2)
      G(5) = RMTOUT(4)
      G(6) = RMTOUT(5)
      G(7) = G(3)
      G(8) = G(6)
      G(9) = RMTOUT(6)
      GO TO 100
C
   60 DO 70 I=1,4
   70 G(I) = RMTOUT(I)
      G(3) = G(2)
      GO TO 100
C
C     ORTHOTROPIC MATERIALS (MAT8)
C
   80 IF (MID .EQ. 3) GO TO 90
      NU21 = NU12 * E2 / E1
      CONST= 1.0 - (NU21*NU12)
      G(1) = E1 / CONST
      G(2) = NU12 * E2 / CONST
      G(4) = G(2)
      G(5) = E2 / CONST
      G(9) = RMTOUT(4)
      GO TO 100
C
   90 G(1) = RMTOUT(6)
      G(4) = RMTOUT(5)
      IF (G(1).EQ.0.0 .AND. G(4).EQ.0.0) GO TO 120
C
C     STANDARD RETURN
C
  100 DO 110 I=1,9
  110 G(I) = G(I)*FACTOR
      GO TO 310
C
C     FATAL RETURN
C
  120 MID = -MID
      GO TO 310
C
      ENTRY Q4GMGD (MID,FACTOD,D)
C     ===========================
C
      DO 200 I=1,9
  200 D(I) = 0.0D0
      MTYPE= RMTOUT(25)
      MTYP = IFIX(MTYPE+.05) - 2
      IF (MTYP) 210,220,270
C
C     ISOTROPIC MATERIALS (MAT1)
C
  210 IF (MID .NE. 3) GO TO 230
      D(1) = RMTOUT(6)
      D(4) = D(1)
      GO TO 290
C
C     ANISOTROPIC MATERIALS (MAT2)
C
  220 IF (MID .EQ. 3) GO TO 250
  230 DO 240 I=1,3
  240 D(I) = RMTOUT(I)
      D(4) = D(2)
      D(5) = RMTOUT(4)
      D(6) = RMTOUT(5)
      D(7) = D(3)
      D(8) = D(6)
      D(9) = RMTOUT(6)
      GO TO 290
C
  250 DO 260 I=1,4
  260 D(I) = RMTOUT(I)
      D(3) = D(2)
      GO TO 290
C
C     ORTHOTROPIC MATERIALS (MAT8)
C
  270 IF (MID .EQ. 3) GO TO 280
      NU21 = NU12 * E2 / E1
      DONST= 1.0D0 - DBLE(NU21*NU12)
      D(1) = E1 / DONST
      D(2) = NU12 * E2 / DONST
      D(4) = D(2)
      D(5) = E2 / DONST
      D(9) = RMTOUT(4)
      GO TO 290
C
  280 D(1) = RMTOUT(6)
      D(4) = RMTOUT(5)
      IF (D(1).EQ.0.0D0 .AND. D(4).EQ.0.0D0) GO TO 120
C
C     STANDARD RETURN
C
  290 DO 300 I=1,9
  300 D(I) = D(I)*FACTOD
C
  310 RETURN
      END


text:
      MODULE GWFHFBMODULE
        INTEGER, POINTER  ::MXHFB,NHFB,IPRHFB,NHFBNP,NPHFB,IHFBPB,
     &                      NumHFBs
        REAL,    DIMENSION(:,:), POINTER   ::HFB
      TYPE GWFHFBTYPE
        INTEGER, POINTER  ::MXHFB,NHFB,IPRHFB,NHFBNP,NPHFB,IHFBPB,
     &                      NumHFBs
        REAL,    DIMENSION(:,:), POINTER   ::HFB
      END TYPE
      TYPE(GWFHFBTYPE), SAVE    ::GWFHFBDAT(10)
      END MODULE GWFHFBMODULE