text:
      subroutine singletopreconstruct(p,failed,costheta,ylight)
c--- Given the usual momentum array, try to reconstruct which
c--- set of momenta reconstruct the top antitop quark;
c--- given those assignments, the routine returns the following quantities:
c---
c---    costheta      the angle between the lepton and the light jet,
c---                    computed in the top rest frame
c---    ylight            the rapidity of the light jet
c---
c--- For real radiation events, the jet algorithm may result in
c--- events for which the invariant mass is not exactly equal to mt,
c--- despite the phase space producing tops exactly on-shell. In that
c--- case, allow |s-mt| <= 'toler' [GeV]
c---
c--- If the event does not contain exactly one b-jet and one light jet
c--- then failed is equal to .true.
      implicit none
      include 'constants.f'
      include 'masses.f'
      include 'plabel.f'
      logical failed
      integer i
      double precision tiny,p(mxpart,4),
     & ptop(4),ptoprest(4),plep(4),plight(4),pleprest(4),plightrest(4),
     & mtcand,costheta,ylight,yrap,ptopalt(4),mtcandalt
      parameter (tiny=1d-8)
c      double precision small,toler
c      parameter (small=1d-4,toler=5d9)

c--- default: reconstruction is okay
      failed=.false.

c--- note: this routine was written for the case 'notag=1' in chooser.f
c--- (the default is notag=0) and inclusive=F, with the result that
c--- only 2 jets are ever accepted in the event
      if (p(5,4) .lt. tiny) then
c        write(6,*) 'No b-quark present'
      failed=.true.
      return
      endif

      if (p(7,4) .lt. tiny) then
c        write(6,*) 'No light jet present'
      failed=.true.
      return
      endif

c--- this requires that two jets are b and light-jet
      do i=1,4
        ptop(i)=p(3,i)+p(4,i)+p(5,i)
        ptopalt(i)=p(3,i)+p(4,i)+p(5,i)+p(7,i)
        plight(i)=p(7,i)
      if (plabel(3) .eq. 'el') then
        plep(i)=p(3,i)
      else
        plep(i)=p(4,i)
      endif
      enddo
      mtcand=dsqrt(ptop(4)**2-ptop(1)**2-ptop(2)**2-ptop(3)**2)
      mtcandalt=dsqrt(ptopalt(4)**2-ptopalt(1)**2
     &               -ptopalt(2)**2-ptopalt(3)**2)
c      write(6,*) 'Event accepted, top candidate mass=',mtcand,mtcandalt

c--- check to see if 345+7 is a better top candidate than 345:
c--- if so, assume light get radiated in decay, so no angle to construct
      if (abs(mtcandalt-mt) .lt. abs(mtcand-mt)) then
c        write(6,*) 'Radiation in decay: ',mtcandalt,' vs ',mtcand
      failed=.true.
      return
      endif

      ptoprest(4)=mtcand
      do i=1,3
      ptoprest(i)=0d0
      enddo

      call boostx(plep,ptop,ptoprest,pleprest)
      call boostx(plight,ptop,ptoprest,plightrest)

c--- cos(theta*)
      costheta=(pleprest(1)*plightrest(1)
     &         +pleprest(2)*plightrest(2)
     &         +pleprest(3)*plightrest(3))
     &        /dsqrt(pleprest(1)**2+pleprest(2)**2+pleprest(3)**2)
     &        /dsqrt(plightrest(1)**2+plightrest(2)**2+plightrest(3)**2)
c      write(6,*) 'cos(theta*)=',costheta

c--- rapidity of light jet
      ylight=yrap(7,p)

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
C***********************************************************************
C*
C*	Copyright 1996, University Corporation for Atmospheric Research.
C*
C*	wprofwnd.f
C*
C*	Called by wind profiler driver to draw wind barbs. Derived from
C*	the GEMPAK program SNSWND.
C*
C*	History:
C*
C*	11/96	COMET	Original copy
C*
C************************************************************************

	SUBROUTINE wprofwnd ( wintyp, nstn, ipsdat, nlvls,
     +			      stndat, idtype, sloc, ybot, ytop, filtfc,
     +			      vcord, iret )
C************************************************************************
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	wintyp, vcord
	REAL		stndat (*), sloc (*)
	INTEGER		ipsdat (*), nlvls (*), idtype (LLMXLV, * )
C*
	CHARACTER	cdata (3)*8
	REAL		data (3), p (2000), d (2000), s (2000),
     +			xx (2000), vx (2000), vy (2000)
csd
        integer set_barb_color, BY_SPEED
        data BY_SPEED/1/
csd        INCLUDE         'DEVSET.CMN'
csd
C*
csd	INCLUDE		'ERMISS.FNC'
C------------------------------------------------------------------------
	iret = 0
C
C*	Get barb size and calculate filter value
C
	CALL GQSYSZ  ( rmkw, rmkh, rtxw, rtxh, brbx, brby, ier )
	brbftr = .6 * brbx * filtfc
C
C*	Extract data for each station.  Don't keep levels which are
C*	outside the plot area; this allows plotting in V coord.,
C*	with barbs extending outside the plot region.  Levels with
C*	missing data are also ignored.  Only winds at mandatory and
C*	significant wind levels are drawn.
C
	DO  i = 1, nstn
	    nlv = 0
	    CALL PC_SSTN  ( '    ',0,0,0,0,ispri,IMISSD, nlvls (i),
     +			iret )
	    DO  j = 1, nlvls (i)
		CALL PC_CMLV (j, stndat (ipsdat (i)), data, cdata, iret)
csd		IF  ( ( .not. ERMISS (data (1)) ) .and.
csd     *		      ( .not. ERMISS (data (2)) ) .and.
csd     *		      ( .not. ERMISS (data (3)) ) .and.
csd     *                ( idtype ( j, i ) .ne. 2  ) .and.
		IF (
     *		      ( ( ( vcord .eq. 'PRES' )   .and.
     *		          ( data (1) .le. ybot )  .and.
     *			  ( data (1) .ge. ytop ) ) .or.
     *			( ( vcord .ne. 'PRES' )   .and.
     *			  ( data (1) .ge. ybot )  .and.
     *			  ( data (1) .le. ytop ) ) ) ) THEN
		    nlv = nlv + 1
		    p (nlv) = data (1)
		    d (nlv) = data (2)
		    s (nlv) = data (3)
		    xx (nlv) = sloc (i)
		END IF
	    END DO
C
C*	    Transform pressures and x locations to V coordinates.
C
	    CALL GTRANS  ( 'M', 'V', nlv, xx, p, vx, vy, iret )
C
C*	    If wind arrows, draw them.
C
	    IF  ( wintyp .eq. 'A' )  THEN
		CALL GARRW   ( 'V', nlv, vx, vy, s, d, iret )
	      ELSE
C
C*		If no filter then draw the barbs.
C
		IF  ( filtfc .eq. 0. )  THEN
csd
csd  Call GBARB one barb at a time to reset color if necessary.
csd
                    do m=1, nlv
                      iwnclr = set_barb_color(s(m),BY_SPEED)
                      CALL GSCOLR  ( iwnclr, ier )
                      CALL GBARB  ( 'V', 1, vx (m), vy (m), s (m),
     *                               d (m), iret )
                    enddo

		  ELSE
C
C*		    Draw the first barb
C
		    yold = vy (1)
                    iwnclr = set_barb_color(s(1),1)
                    CALL GSCOLR  ( iwnclr, ier )
		    CALL GBARB  ( 'V', 1, vx (1), vy (1), s (1), d (1),
     +				  iret )
C
C*		    Check for overlap and draw if ok.
C
		    DO  m = 2, nlv
			IF  ( vy (m) - yold .ge. brbftr )  THEN
                            iwnclr = set_barb_color(s(m),1)
                            CALL GSCOLR  ( iwnclr, ier )
			    CALL GBARB  ( 'V', 1, vx (m), vy (m),
     +					  s (m), d (m), iret )
			    yold = vy (m)
			END IF
		    END DO
		END IF
	    END IF
	END DO
C*
	RETURN
	END

ext: f

max_stars_repo_name: Unidata/Garp

max_issues_repo_name: Unidata/Garp

max_forks_repo_name: Unidata/Garp


! { dg-do compile }
! { dg-options "-Wunused" }
!
! PR fortran/31461
!
! Contributed by Vivek Rao.
!

module util_mod
  integer :: i,j
end module util_mod

program main
  use util_mod, only: i,j ! { dg-warning "Unused module variable .i. which has been explicitly imported" }
  j = 1
  print*,"j=",j
end program main


text:
module message_test
    use iso_varying_string, only: VARYING_STRING, operator(//), var_str
    use erloff_debug_m, only: debug_t, DEBUG
    use erloff_debug_level_m, only: GENERAL
    use erloff_error_m, only: ERROR
    use erloff_fatal_m, only: fatal_t, FATAL
    use erloff_info_m, only: info_t, INFO
    use erloff_internal_m, only: internal_t, INTERNAL
    use erloff_message_m, only: Message_t
    use erloff_message_type_m, only: &
            INPUTS, &
            OUTSIDE_NORMAL_RANGE, &
            UNEQUAL_ARRAY_SIZES, &
            UNKNOWN_TYPE
    use erloff_module_m, only: module_t
    use erloff_procedure_m, only: procedure_t
    use erloff_warning_m, only: warning_t, WARNING
    use veggies, only: &
            Result_t, &
            Test_Item_t, &
            Describe, &
            assert_Not, &
            assert_That, &
            It

    implicit none
    private

    public :: test_message
contains
    function test_message() result(tests)
        type(Test_Item_t) :: tests

        type(Test_Item_t) :: individual_tests(8)

        individual_tests(1) = It( &
                "Can tell whether it is of a given type", checkType)
        individual_tests(2) = It( &
                "Can tell which module it originated from", &
                checkOriginatingModule)
        individual_tests(3) = It( &
                "Can tell which procedure it originated from", &
                checkOriginatingProcedure)
        individual_tests(4) = It( &
                "Can tell which modules it came from", &
                checkFromModule)
        individual_tests(5) = It( &
                "Can tell which procedures it came from", &
                checkFromProcedure)
        individual_tests(6) = It( &
                "Can tell which modules it came through", &
                checkThroughModule)
        individual_tests(7) = It( &
                "Can tell which procedures it came through", &
                checkThroughProcedure)
        individual_tests(8) = It( &
                "Can tell lots about its contents", &
                checkContents)
        tests = Describe("Message_t", individual_tests)
    end function test_message

    pure function checkType() result(result_)
        type(Result_t) :: result_

        class(Message_t), allocatable :: debug_message
        class(Message_t), allocatable :: info_message
        class(Message_t), allocatable :: warning_message
        class(Message_t), allocatable :: fatal_message
        class(Message_t), allocatable :: internal_message

        debug_message = debug_t( &
                INPUTS, &
                module_t("Some_m"), &
                procedure_t("some"), &
                GENERAL, &
                "Test Message")

        info_message = info_t( &
                UNEQUAL_ARRAY_SIZES, &
                module_t("Some_m"), &
                procedure_t("some"), &
                "Test Message")

        warning_message = warning_t( &
                OUTSIDE_NORMAL_RANGE, &
                module_t("Some_m"), &
                procedure_t("some"), &
                "Test Message")

        fatal_message = fatal_t( &
                UNEQUAL_ARRAY_SIZES, &
                module_t("Some_m"), &
                procedure_t("some"), &
                "Test Message")

        internal_message = internal_t( &
                UNKNOWN_TYPE, &
                module_t("Some_m"), &
                procedure_t("some"), &
                "Test Message")

        result_ = &
                assert_That( &
                        debug_message.isType.DEBUG, &
                        debug_message%repr() // ".isType." // DEBUG%repr()) &
                .and.assert_That( &
                        debug_message.isType.INPUTS, &
                        debug_message%repr() // ".isType." // INPUTS%repr()) &
                .and.assert_Not( &
                        debug_message.isType.INFO, &
                        debug_message%repr() // ".isType." // INFO%repr()) &
                .and.assert_Not( &
                        debug_message.isType.ERROR, &
                        debug_message%repr() // ".isType." // ERROR%repr()) &
                .and.assert_That( &
                        warning_message.isType.WARNING, &
                        warning_message%repr() // ".isType." // WARNING%repr()) &
                .and.assert_That( &
                        warning_message.isType.OUTSIDE_NORMAL_RANGE, &
                        warning_message%repr() // ".isType." // OUTSIDE_NORMAL_RANGE%repr()) &
                .and.assert_Not( &
                        warning_message.isType.INFO, &
                        warning_message%repr() // ".isType." // INFO%repr()) &
                .and.assert_Not( &
                        warning_message.isType.ERROR, &
                        warning_message%repr() // ".isType." // ERROR%repr()) &
                .and.assert_That( &
                        fatal_message.isType.ERROR, &
                        fatal_message%repr() // ".isType." // ERROR%repr()) &
                .and.assert_That( &
                        fatal_message.isType.FATAL, &
                        fatal_message%repr() // ".isType." // FATAL%repr()) &
                .and.assert_That( &
                        fatal_message.isType.UNEQUAL_ARRAY_SIZES, &
                        fatal_message%repr() // ".isType." // UNEQUAL_ARRAY_SIZES%repr()) &
                .and.assert_Not( &
                        fatal_message.isType.INTERNAL, &
                        fatal_message%repr() // ".isType." // INTERNAL%repr()) &
                .and.assert_That( &
                        internal_message.isType.ERROR, &
                        internal_message%repr() // ".isType." // ERROR%repr()) &
                .and.assert_That( &
                        internal_message.isType.INTERNAL, &
                        internal_message%repr() // ".isType." // INTERNAL%repr()) &
                .and.assert_That( &
                        internal_message.isType.UNKNOWN_TYPE, &
                        internal_message%repr() // ".isType." // UNKNOWN_TYPE%repr())
    end function checkType

    function checkOriginatingModule() result(result_)
        type(Result_t) :: result_

        type(Module_t) :: another_module
        type(procedure_t) :: another_procedure
        class(Message_t), allocatable :: message
        class(Message_t), allocatable :: message2
        type(Module_t) :: other_module
        type(Module_t) :: the_module
        type(procedure_t) :: the_procedure

        the_module = module_t("Some_m")
        the_procedure = procedure_t("some")
        another_module = module_t("Another_m")
        another_procedure = procedure_t("another")
        other_module = module_t("Other_m")
        message2 = info_t(the_module, the_procedure, "Test Message")
        message = message2%with_names_prepended(another_module, another_procedure)

        result_ = &
                assert_That( &
                        message.originatedFrom.the_module, &
                        message%repr() // '.originatedFrom.' // the_module%repr()) &
                .and.assert_Not( &
                        message.originatedFrom.another_module, &
                        message%repr() // '.originatedFrom.' // another_module%repr()) &
                .and.assert_Not( &
                        message.originatedFrom.other_module, &
                        message%repr() // '.originatedFrom.' // other_module%repr())
    end function checkOriginatingModule

    function checkOriginatingProcedure() result(result_)
        type(Result_t) :: result_

        type(module_t) :: another_module
        type(procedure_t) :: another_procedure
        class(Message_t), allocatable :: message
        class(Message_t), allocatable :: message2
        type(procedure_t) :: other_procedure
        type(module_t) :: the_module
        type(procedure_t) :: the_procedure

        the_module = module_t("Some_m")
        the_procedure = procedure_t("some")
        another_module = module_t("Another_m")
        another_procedure = procedure_t("another")
        other_procedure = procedure_t("other")
        message2 = info_t(the_module, the_procedure, "Test Message")
        message = message2%with_names_prepended(another_module, another_procedure)

        result_ = &
                assert_That( &
                        message.originatedFrom.the_procedure, &
                        message%repr() // '.originatedFrom.' // the_procedure%repr()) &
                .and.assert_Not( &
                        message.originatedFrom.another_procedure, &
                        message%repr() // '.originatedFrom.' // another_procedure%repr()) &
                .and.assert_Not( &
                        message.originatedFrom.other_procedure, &
                        message%repr() // '.originatedFrom.' // other_procedure%repr())
    end function checkOriginatingProcedure

    function checkFromModule() result(result_)
        type(Result_t) :: result_

        type(Module_t) :: another_module
        type(procedure_t) :: another_procedure
        class(Message_t), allocatable :: message
        class(Message_t), allocatable :: message2
        type(Module_t) :: other_module
        type(Module_t) :: the_module
        type(procedure_t) :: the_procedure

        the_module = module_t("Some_m")
        the_procedure = procedure_t("some")
        another_module = module_t("Another_m")
        another_procedure = procedure_t("another")
        other_module = module_t("Other_m")
        message2 = info_t(the_module, the_procedure, "Test Message")
        message = message2%with_names_prepended(another_module, another_procedure)

        result_ = &
                assert_That( &
                        message.isFrom.the_module, &
                        message%repr() // '.isFrom.' // the_module%repr()) &
                .and.assert_That( &
                        message.isFrom.another_module, &
                        message%repr() // '.isFrom.' // another_module%repr()) &
                .and.assert_Not( &
                        message.isFrom.other_module, &
                        message%repr() // '.isFrom.' // other_module%repr())
    end function checkFromModule

    function checkFromProcedure() result(result_)
        type(Result_t) :: result_

        type(Module_t) :: another_module
        type(procedure_t) :: another_procedure
        class(Message_t), allocatable :: message
        class(Message_t), allocatable :: message2
        type(procedure_t) :: other_procedure
        type(Module_t) :: the_module
        type(procedure_t) :: the_procedure

        the_module = module_t("Some_m")
        the_procedure = procedure_t("some")
        another_module = module_t("Another_m")
        another_procedure = procedure_t("another")
        other_procedure = procedure_t("other")
        message2 = info_t(the_module, the_procedure, "Test Message")
        message = message2%with_names_prepended(another_module, another_procedure)

        result_ = &
                assert_That( &
                        message.isFrom.the_procedure, &
                        message%repr() // '.isFrom.' // the_procedure%repr()) &
                .and.assert_That( &
                        message.isFrom.another_procedure, &
                        message%repr() // '.isFrom.' // another_procedure%repr()) &
                .and.assert_Not( &
                        message.isFrom.other_procedure, &
                        message%repr() // '.isFrom.' // other_procedure%repr())
    end function checkFromProcedure

    function checkThroughModule() result(result_)
        type(Result_t) :: result_

        type(Module_t) :: another_module
        type(procedure_t) :: another_procedure
        class(Message_t), allocatable :: message
        class(Message_t), allocatable :: message2
        type(Module_t) :: other_module
        type(Module_t) :: the_module
        type(procedure_t) :: the_procedure

        the_module = module_t("Some_m")
        the_procedure = procedure_t("some")
        another_module = module_t("Another_m")
        another_procedure = procedure_t("another")
        other_module = module_t("Other_m")
        message2 = info_t(the_module, the_procedure, "Test Message")
        message = message2%with_names_prepended(another_module, another_procedure)

        result_ = &
                assert_Not( &
                        message.cameThrough.the_module, &
                        message%repr() // '.cameThrough.' // the_module%repr()) &
                .and.assert_That( &
                        message.cameThrough.another_module, &
                        message%repr() // '.cameThrough.' // another_module%repr()) &
                .and.assert_Not( &
                        message.cameThrough.other_module, &
                        message%repr() // '.cameThrough.' // other_module%repr())
    end function checkThroughModule

    function checkThroughProcedure() result(result_)
        type(Result_t) :: result_

        type(Module_t) :: another_module
        type(procedure_t) :: another_procedure
        class(Message_t), allocatable :: message
        class(Message_t), allocatable :: message2
        type(procedure_t) :: other_procedure
        type(Module_t) :: the_module
        type(procedure_t) :: the_procedure

        the_module = module_t("Some_m")
        the_procedure = procedure_t("some")
        another_module = module_t("Another_m")
        another_procedure = procedure_t("another")
        other_procedure = procedure_t("other")
        message2 = info_t(the_module, the_procedure, "Test Message")
        message = message2%with_names_prepended(another_module, another_procedure)

        result_ = &
                assert_Not( &
                        message.cameThrough.the_procedure, &
                        message%repr() // '.cameThrough.' // the_procedure%repr()) &
                .and.assert_That( &
                        message.cameThrough.another_procedure, &
                        message%repr() // '.cameThrough.' // another_procedure%repr()) &
                .and.assert_Not( &
                        message.cameThrough.other_procedure, &
                        message%repr() // '.cameThrough.' // other_procedure%repr())
    end function checkThroughProcedure

    function checkContents() result(result_)
        type(Result_t) :: result_

        class(Message_t), allocatable :: message
        type(VARYING_STRING) :: includesAnyStrings1(2)
        type(VARYING_STRING) :: includesAnyStrings2(2)
        type(VARYING_STRING) :: includesAnyStrings3(2)
        type(VARYING_STRING) :: includesAllStrings1(2)
        type(VARYING_STRING) :: includesAllStrings2(2)
        type(VARYING_STRING) :: includesAllStrings3(2)

        message = info_t( &
                module_t("Some_m"), procedure_t("some"), "Test Message Content")
        includesAnyStrings1(1) = var_str("Test")
        includesAnyStrings1(2) = var_str("else")
        includesAnyStrings2(1) = var_str("Test")
        includesAnyStrings2(2) = var_str("Content")
        includesAnyStrings3(1) = var_str("test")
        includesAnyStrings3(2) = var_str("else")
        includesAllStrings1(1) = var_str("Test")
        includesAllStrings1(2) = var_str("Message")
        includesAllStrings2(1) = var_str("test")
        includesAllStrings2(2) = var_str("Message")
        includesAllStrings3(1) = var_str("Test")
        includesAllStrings3(2) = var_str("message")

        result_ = &
                assert_That( &
                        message.includes."Test", &
                        message%repr() // '.includes."Test"') &
                .and.assert_Not( &
                        message.includes."test", &
                        message%repr() // '.includes."test"') &
                .and.assert_That( &
                        message.includesAnyOf.includesAnyStrings1, &
                        message%repr() // '.includesAnyOf.[var_str("Test"), var_str("else")]') &
                .and.assert_That( &
                        message.includesAnyOf.includesAnyStrings2, &
                        message%repr() // '.includesAnyOf.[var_str("Test"), var_str("Content")]') &
                .and.assert_Not( &
                        message.includesAnyOf.includesAnyStrings3, &
                        message%repr() // '.includesAnyOf.[var_str("test"), var_str("else")]') &
                .and.assert_That( &
                        message.includesAllOf.includesAllStrings1, &
                        message%repr() // '.includesAllOf.[var_str("Test"), var_str("Message")]') &
                .and.assert_Not( &
                        message.includesAllOf.includesAllStrings2, &
                        message%repr() // '.includesAllOf.[var_str("test"), var_str("Message")]') &
                .and.assert_Not( &
                        message.includesAllOf.includesAllStrings3, &
                        message%repr() // '.includesAllOf.[var_str("Test"), var_str("message")]')
    end function checkContents
end module message_test

max_stars_repo_name: everythingfunctional/erloff

max_issues_repo_name: everythingfunctional/erloff

max_forks_repo_name: everythingfunctional/erloff


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 04.03.02.02/10                        *
C  *    TEST TITLE : Foreground colour                     *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /GLOBCH/ PIDENT,    GLBERR,    TSTMSG,     FUNCID,
     1                DUMCH
      CHARACTER       PIDENT*40, GLBERR*60, TSTMSG*900, FUNCID*80,
     1                DUMCH(20)*20

      COMMON /DIALOG/ DOUTYP, DINTYP, DSTDNR, DSTRID, PSTRID, DTCLIM,
     1                SCRMOD, DTXCI, SPECWT,
     2                DSIZE, EFRAC, DYXRAT, SYXRAT, MTRPDC, WCPDC, QVIS
      INTEGER         DOUTYP, DINTYP, DSTDNR, DSTRID, PSTRID, DTCLIM,
     1                SCRMOD, DTXCI, SPECWT
      REAL            DSIZE, EFRAC, DYXRAT, SYXRAT, MTRPDC, WCPDC, QVIS

C aspect source
C                bundled     individual
      INTEGER    PBUNDL,     PINDIV
      PARAMETER (PBUNDL = 0, PINDIV = 1)
C marker type
      INTEGER    PPOINT,   PPLUS,   PAST,   POMARK,   PXMARK
      PARAMETER (PPOINT=1, PPLUS=2, PAST=3, POMARK=4, PXMARK=5)
C  type of returned value
      INTEGER    PSET,       PREALI
      PARAMETER (PSET = 0,   PREALI = 1)
C  composition type
      INTEGER    PCPRE,      PCPOST,    PCREPL
      PARAMETER (PCPRE = 0,  PCPOST = 1,PCREPL = 2)
C  colour model
      INTEGER    PRGB,       PCIE,      PHSV,     PHLS
      PARAMETER (PRGB = 1,   PCIE = 2,  PHSV = 3, PHLS = 4)

      INTEGER    PICSTR, TXCI, IX, BGI, NUMRET, COLSIZ, MAXCI, SQSIDE
      INTEGER    NGSQ, WINDOW, NXTCOL, PERM(6), RNDINT
      INTEGER    IDUM1, IDUM2, IDUM3, IDUM4, IDUM5, IDUM6, IDUM7

      REAL       GBM(3,3), FG(3), BG(3),BCKCOL(3),NOMMS
      REAL       SIDE1, SIDE2, SQINC, MRKSIZ, XSIZ, YSIZ
      REAL       XLOC, YLOC, Z, U, H, XWINLO(6), YWINLO(6), XF(3,3)
      REAL       XA(100), YA(100), RDUM2, RDUM3

      PARAMETER  (Z = 0.0, U = 1.0, H = 0.5)

      DATA GBM / 0.,1.,0.,0.,0.,1.,1.,0.,1. /

      CALL INITGL ('04.03.02.02/10')

C open PHIGS
      CALL XPOPPH (ERRFIL, MEMUN)
C set-up of workstation and dialogue area
      PICSTR = 101
      TXCI = 1
      CALL SETDLG (PICSTR, 801,TXCI)

      CALL POPST (PICSTR)
C  by convention, view #1 is for picture
      CALL PSVWI (1)
C  use individual attributes
      CALL SETASF (PINDIV)
C  set colour model = RGB
      CALL PSCMD (WKID, PRGB)

      CALL PEXST (106)
      CALL PEXST (102)
      CALL PCLST

C  divide screen up into six labelled square areas
      CALL WIN6 (106, 1, XSIZ, YSIZ, XWINLO, YWINLO)

      CALL POPST (102)

C  bg = requested background colour = randomly: green, blue or magenta
      BGI = RNDINT (1, 3)
      BG(1) = GBM(1,BGI)
      BG(2) = GBM(2,BGI)
      BG(3) = GBM(3,BGI)
C  set color representation #0 to bg
      CALL PSCR (WKID, 0, 3, BG)

C  inquire color representation #0 as realized to determine
C    bckcol = actual background color
      CALL PQCR (WKID, 0, 3, PREALI, ERRIND, NUMRET, BCKCOL)
      CALL CHKINQ ('pqcr', ERRIND)

C  inquire workstation state table lengths to determine
C    colsiz = number of settable entries in color table
      CALL PQWKSL (SPECWT, ERRIND, IDUM1, IDUM2, IDUM3, IDUM4,
     1             IDUM5, IDUM6, COLSIZ, IDUM7)
      CALL CHKINQ ('pqwksl', ERRIND)
C  maxci = maximum color index = colsiz-1
      MAXCI = COLSIZ - 1

C  fg = foreground color farthest from bckcol in color cube
      DO 50 IX = 1, 3
         IF (BCKCOL(IX) .LT. 0.5) THEN
            FG(IX) = 1.0
         ELSE
            FG(IX) = 0.0
         ENDIF
50    CONTINUE

C  set color representations from #1 to #maxci to fg
      DO 100 IX = 1, MAXCI
         CALL PSCR (WKID, IX, 3, FG)
100   CONTINUE

      CALL SETMSG ('30 31 32', 'All the positive entries within the ' //
     1             'reported size of the color table should control ' //
     2             'a foreground colour.')

C  sqside = integer number of polymarkers per side of square array
C         = sqrt(colsiz/4 + 1) + 1
C  (this ensures all color slots will be used)
C  side1 = 0.1
C  side2 = 0.9
C  sqinc = increment for square = (side2-side1) / (sqside-1)
C  mrksiz = polymarker size = min(0.05, 0.2 * sqinc)
      SQSIDE = SQRT(COLSIZ/4. + 1.) + 1
      SIDE1 = 0.1
      SIDE2 = 0.9
      SQINC = (SIDE2 - SIDE1) / (SQSIDE - 1)
      MRKSIZ = MIN(0.05, 0.2*SQINC)

      CALL PQPMF (SPECWT, 0, ERRIND, IDUM1, IDUM2, IDUM3, NOMMS,
     1            RDUM2, RDUM3, IDUM4)
      CALL CHKINQ ('pqpmf', ERRIND)

C  set polymarker style = plus (+)
C  set polymarker size = mrksiz (in WC)
      CALL PSMK (PPLUS)
      CALL PSMKSC (MRKSIZ/(NOMMS * WCPDC))

C  perm = randomize order from 1 to 6
      CALL RNPERM (6, PERM)
C  nxtcol = next color index to use = maxci
      NXTCOL = MAXCI

C  draw square grids of polymarkers
      DO 200 IX = 1, 6
         WINDOW = PERM(IX)
C     set tranformation to scale unit square into this window
         CALL EBLTM (Z, Z, XWINLO(IX), YWINLO(IX), Z, XSIZ, XSIZ, XF)
         CALL PSLMT (XF, PCREPL)

         IF (WINDOW .EQ. 1) THEN
C  simulate incorrect display
            NGSQ = IX
            CALL PSPMCI (1)

            CALL SQGRMK (SIDE1, SIDE2+0.5*SQINC, SQINC)
C  put a polymarker at a non-aligned position
            XA(1) = SIDE1 + SQINC * (RNDINT(0, SQSIDE-2) + 0.5)
            YA(1) = SIDE1 + SQINC * (RNDINT(0, SQSIDE-2) + 0.5)

            CALL PPM (1, XA, YA)
         ELSEIF (WINDOW .EQ. 2) THEN
C  simulate correct display
            CALL PSPMCI (1)
            CALL SQGRMK (SIDE1, SIDE2+0.5*SQINC, SQINC)
         ELSE
C  draw an actual grid
            DO 230 XLOC = SIDE1, SIDE2 + 0.5*SQINC, SQINC
            DO 235 YLOC = SIDE1, SIDE2 + 0.5*SQINC, SQINC
               IF (NXTCOL .LT. MAXCI) THEN
                  NXTCOL = NXTCOL + 1
               ELSE
                  NXTCOL = 1
               ENDIF
C  display single polymarker: xa,ya
               CALL PSPMCI (NXTCOL)
               XA(1) = XLOC
               YA(1) = YLOC
               CALL PPM (1, XA, YA)
235         CONTINUE
230         CONTINUE
         ENDIF
200   CONTINUE

      CALL DCHPFV ('FOREGROUND COLOUR: Which polymarker is ' //
     1             'different in either color or geometry?',
     2             6, NGSQ)
      CALL PEMST (102)

666   CONTINUE
C  wrap it up.
      CALL ENDIT
      END


text:
C-----------------------------------------------------------------------
!
!  revision history:
!
!  Aug, 2012  Jun Wang    bafrio for big and little endian files
!
!  note:
!     This version of bafrio.f is revised to have byteswap in FORTRAN
!  data file control words. It is designed to be run on
!  on WCOSS(little endian machine) and to generate big endian files.
!     It does byteswap on fortran record control words(4 byte integer
!  before and after data field), not on data field itself. Users need
!  to byteswap their data after(for reading)/before(for writing)
!  calling subroutines this file. This is considered to be the best
!  way to keep subroutine inerfaces intact for backward compatible.
!
C-----------------------------------------------------------------------
      SUBROUTINE BAFRINDEX(LU,IB,LX,IX)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: BAFRINDEX      BYTE-ADDRESSABLE FORTRAN RECORD INDEX
C   PRGMMR: IREDELL          ORG: W/NMC23     DATE: 1999-01-21
C
C ABSTRACT: THIS SUBPROGRAM IS CALLING BAFRINDEXL TO EITHER READ AN
C   UNFORMATTED FORTRAN RECORD
C   AND RETURN ITS LENGTH AND START BYTE OF THE NEXT FORTRAN RECORD;
C   OR GIVEN THE RECORD LENGTH, WITHOUT I/O IT DETERMINES THE START BYTE
C   OF THE NEXT FORTRAN RECORD. THE DIFFERENCE BETWEEN BAFRINDEX AND
C   BAFRINDEXL IS THE KIND TYPE OF INTERGERS IN THE ARGUMENT LIST
C
C PROGRAM HISTORY LOG:
C   1999-01-21  IREDELL
C   2009-04-20  J. WANG
C
C USAGE:    CALL BAFRINDEX(LU,IB,LX,IX)
C   INPUT ARGUMENTS:
C     LU           INTEGER LOGICAL UNIT TO READ
C                  IF LU<=0, THEN DETERMINE IX FROM LX
C     IB           INTEGER FORTRAN RECORD START BYTE
C                  (FOR THE FIRST FORTRAN RECORD, IB SHOULD BE 0)
C     LX           INTEGER RECORD LENGTH IN BYTES IF LU<=0
C
C   OUTPUT ARGUMENTS:
C     LX           INTEGER RECORD LENGTH IN BYTES IF LU>0,
C                  OR LX=-1 FOR I/O ERROR (PROBABLE END OF FILE),
C                  OR LX=-2 FOR I/O ERROR (INVALID FORTRAN RECORD)
C     IX           INTEGER START BYTE FOR THE NEXT FORTRAN RECORD
C                  (COMPUTED ONLY IF LX>=0)
C
C SUBPROGRAMS CALLED:
C   BAFRINDEXL     BYTE-ADDRESSABLE FORTRAN RECORD INDEX
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C
C$$$
!
      IMPLICIT NONE
      INTEGER,INTENT(IN):: LU,IB
      INTEGER,INTENT(INOUT):: LX
      INTEGER,INTENT(OUT):: IX
      integer(kind=8) :: LONG_IB,LONG_LX ,LONG_IX
!
      LONG_IB=IB
      LONG_LX=LX
      call BAFRINDEXL(LU,LONG_IB,LONG_LX,LONG_IX)
      LX=LONG_LX
      IX=LONG_IX

      return
      end SUBROUTINE BAFRINDEX
C-----------------------------------------------------------------------
      SUBROUTINE BAFRINDEXL(LU,IB,LX,IX)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: BAFRINDEXL      BYTE-ADDRESSABLE FORTRAN RECORD INDEX
C   PRGMMR: IREDELL          ORG: W/NMC23     DATE: 1999-01-21
C
C ABSTRACT: THIS SUBPROGRAM EITHER READS AN UNFORMATTED FORTRAN RECORD
C   AND RETURN ITS LENGTH AND START BYTE OF THE NEXT FORTRAN RECORD;
C   OR GIVEN THE RECORD LENGTH, WITHOUT I/O IT DETERMINES THE START BYTE
C   OF THE NEXT FORTRAN RECORD.
C
C PROGRAM HISTORY LOG:
C   1999-01-21  IREDELL
C   2009-04-20  J. WANG
C
C USAGE:    CALL BAFRINDEXL(LU,IB,LX,IX)
C   INPUT ARGUMENTS:
C     LU           INTEGER LOGICAL UNIT TO READ
C                  IF LU<=0, THEN DETERMINE IX FROM LX
C     IB           INTEGER(8) FORTRAN RECORD START BYTE
C                  (FOR THE FIRST FORTRAN RECORD, IB SHOULD BE 0)
C     LX           INTEGER(8) RECORD LENGTH IN BYTES IF LU<=0
C
C   OUTPUT ARGUMENTS:
C     LX           INTEGER(8) RECORD LENGTH IN BYTES IF LU>0,
C                  OR LX=-1 FOR I/O ERROR (PROBABLE END OF FILE),
C                  OR LX=-2 FOR I/O ERROR (INVALID FORTRAN RECORD)
C     IX           INTEGER(8) START BYTE FOR THE NEXT FORTRAN RECORD
C                  (COMPUTED ONLY IF LX>=0)
C
C SUBPROGRAMS CALLED:
C   BAREADL        BYTE-ADDRESSABLE READ
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C
C$$$
      IMPLICIT NONE
      INTEGER,INTENT(IN):: LU
      INTEGER(KIND=8),INTENT(IN):: IB
      INTEGER(KIND=8),INTENT(INOUT):: LX
      INTEGER(KIND=8),INTENT(OUT):: IX
      INTEGER(KIND=8),PARAMETER:: LBCW=4
      INTEGER(KIND=LBCW):: BCW1,BCW2
      INTEGER(KIND=8):: KR
      CHARACTER(16) :: MACHINE_ENDIAN
      LOGICAL :: DO_BYTESWAP
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  COMPARE FIRST BLOCK CONTROL WORD AND TRAILING BLOCK CONTROL WORD
      IF(LU.GT.0) THEN
!
!-- set do_byteswap from machine endianness and file endianness
        CALL CHK_ENDIANC(MACHINE_ENDIAN)
        IF( LU<=999) THEN
          IF( trim(MACHINE_ENDIAN)=="big_endian") THEN
            DO_BYTESWAP=.false.
          ELSEIF( trim(MACHINE_ENDIAN)=="little_endian") THEN
            DO_BYTESWAP=.true.
          ENDIF
        ELSEIF(LU<=1999) THEN
          IF( trim(MACHINE_ENDIAN)=="big_endian") THEN
            DO_BYTESWAP=.true.
          ELSEIF( trim(MACHINE_ENDIAN)=="little_endian") THEN
            DO_BYTESWAP=.false.
          ENDIF
        ENDIF
!
!
!-- read out control word
        CALL BAREADL(LU,IB,LBCW,KR,BCW1)
        IF(DO_BYTESWAP) CALL Byteswap(BCW1,LBCW,1)
!
        IF(KR.NE.LBCW) THEN
          LX=-1
        ELSE
          CALL BAREADL(LU,IB+LBCW+BCW1,LBCW,KR,BCW2)
          IF(DO_BYTESWAP) CALL Byteswap(BCW2,LBCW,1)
!
          IF(KR.NE.LBCW.OR.BCW1.NE.BCW2) THEN
            LX=-2
          ELSE
            LX=BCW1
          ENDIF
        ENDIF
!
!end luif
      ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  COMPUTE START BYTE FOR THE NEXT FORTRAN RECORD
      IF(LX.GE.0) IX=IB+LBCW+LX+LBCW
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      END SUBROUTINE BAFRINDEXL
C-----------------------------------------------------------------------
      SUBROUTINE BAFRREAD(LU,IB,NB,KA,A)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: BAFRREAD       BYTE-ADDRESSABLE FORTRAN RECORD READ
C   PRGMMR: IREDELL          ORG: W/NMC23     DATE: 1999-01-21
C
C ABSTRACT: THIS SUBPROGRAM IS CALLING BAFREAD TO REAS AN UNFORMATTED
C   FORTRAN RECORD. THE DIFFERENCE BETWEEN BAFRREAD AND BAFRREADL IS
C    THE KIND TYPE OF INTERGERS IN THE ARGUMENT LIST
C
C PROGRAM HISTORY LOG:
C   1999-01-21  IREDELL
C   2009-04-20  J. WANG
C
C USAGE:    CALL BAFRREAD(LU,IB,NB,KA,A)
C   INPUT ARGUMENTS:
C     LU           INTEGER LOGICAL UNIT TO READ
C     IB           INTEGER FORTRAN RECORD START BYTE
C                  (FOR THE FIRST FORTRAN RECORD, IB SHOULD BE 0)
C     NB           INTEGER NUMBER OF BYTES TO READ
C
C   OUTPUT ARGUMENTS:
C     KA           INTEGER NUMBER OF BYTES IN FORTRAN RECORD
C                  (IN WHICH CASE THE NEXT FORTRAN RECORD
C                  SHOULD HAVE A START BYTE OF IB+KA),
C                  OR KA=-1 FOR I/O ERROR (PROBABLE END OF FILE),
C                  OR KA=-2 FOR I/O ERROR (INVALID FORTRAN RECORD),
C                  OR KA=-3 FOR I/O ERROR (REQUEST LONGER THAN RECORD)
C     A            CHARACTER*1 (NB) DATA READ
C
C SUBPROGRAMS CALLED:
C   BAFRREADL      BYTE-ADDRESSABLE READ
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C
C$$$
!
      IMPLICIT NONE
      INTEGER,INTENT(IN):: LU,IB,NB
      INTEGER,INTENT(OUT):: KA
      CHARACTER,INTENT(OUT):: A(NB)
      INTEGER(KIND=8) :: LONG_IB,LONG_NB,LONG_KA
!
        if((IB<0.and.IB/=-1) .or. NB<0 ) THEN
          print *,'WRONG: in BAFRREAD starting postion IB or read '//    &
     & 'data size NB < 0, STOP! Consider use BAFREADL and long integer'
          KA=0
          return
        ENDIF
        LONG_IB=IB
        LONG_NB=NB
        CALL BAFRREADL(LU,LONG_IB,LONG_NB,LONG_KA,A)
        KA=LONG_KA
      END SUBROUTINE BAFRREAD
C-----------------------------------------------------------------------
      SUBROUTINE BAFRREADL(LU,IB,NB,KA,A)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: BAFRREADL      BYTE-ADDRESSABLE FORTRAN RECORD READ
C   PRGMMR: IREDELL          ORG: W/NMC23     DATE: 1999-01-21
C
C ABSTRACT: THIS SUBPROGRAM READS AN UNFORMATTED FORTRAN RECORD
C
C PROGRAM HISTORY LOG:
C   1999-01-21  IREDELL
C   2009-04-20  J. WANG
C
C USAGE:    CALL BAFRREADL(LU,IB,NB,KA,A)
C   INPUT ARGUMENTS:
C     LU           INTEGER LOGICAL UNIT TO READ
C     IB           INTEGER(8) FORTRAN RECORD START BYTE
C                  (FOR THE FIRST FORTRAN RECORD, IB SHOULD BE 0)
C     NB           INTEGER(8) NUMBER OF BYTES TO READ
C
C   OUTPUT ARGUMENTS:
C     KA           INTEGER(8) NUMBER OF BYTES IN FORTRAN RECORD
C                  (IN WHICH CASE THE NEXT FORTRAN RECORD
C                  SHOULD HAVE A START BYTE OF IB+KA),
C                  OR KA=-1 FOR I/O ERROR (PROBABLE END OF FILE),
C                  OR KA=-2 FOR I/O ERROR (INVALID FORTRAN RECORD),
C                  OR KA=-3 FOR I/O ERROR (REQUEST LONGER THAN RECORD)
C     A            CHARACTER*1 (NB) DATA READ
C
C SUBPROGRAMS CALLED:
C   BAFRINDEXL     BYTE-ADDRESSABLE FORTRAN RECORD INDEX
C   BAREADL        BYTE-ADDRESSABLE READ
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C
C$$$
      IMPLICIT NONE
      INTEGER,INTENT(IN):: LU
      INTEGER(kind=8),INTENT(IN):: IB,NB
      INTEGER(kind=8),INTENT(OUT):: KA
      CHARACTER,INTENT(OUT):: A(NB)
      INTEGER(kind=8),PARAMETER:: LBCW=4
      INTEGER(kind=8):: LX,IX
      INTEGER(kind=8):: KR
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  VALIDATE FORTRAN RECORD
      CALL BAFRINDEXL(LU,IB,LX,IX)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  READ IF VALID
      IF(LX.LT.0) THEN
        KA=LX
      ELSEIF(LX.LT.NB) THEN
        KA=-3
      ELSE
        CALL BAREADL(LU,IB+LBCW,NB,KR,A)
        IF(KR.NE.NB) THEN
          KA=-1
        ELSE
          KA=LBCW+LX+LBCW
        ENDIF
      ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      END SUBROUTINE BAFRREADL
C-----------------------------------------------------------------------
      SUBROUTINE BAFRWRITE(LU,IB,NB,KA,A)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: BAFRWRITE      BYTE-ADDRESSABLE FORTRAN RECORD WRITE
C   PRGMMR: IREDELL          ORG: W/NMC23     DATE: 1999-01-21
C
C ABSTRACT: THIS SUBPROGRAM IS CALLING BAFRWRITE TO WRITE AN UNFORMATTED
C    FORTRAN RECORD. THE DIFFERENCE BETWEEN BAFRWRITE AND BAFRWRITEL IS
C    THE KIND TYPE OF INTERGERS IN THE ARGUMENT LIST
C
C PROGRAM HISTORY LOG:
C   1999-01-21  IREDELL
C   2009-04-20  J. WANG
C
C USAGE:    CALL BAFRWRITE(LU,IB,NB,KA,A)
C   INPUT ARGUMENTS:
C     LU           INTEGER LOGICAL UNIT TO WRITE
C     IB           INTEGER FORTRAN RECORD START BYTE
C                  (FOR THE FIRST FORTRAN RECORD, IB SHOULD BE 0)
C     NB           INTEGER NUMBER OF BYTES TO WRITE
C     A            CHARACTER*1 (NB) DATA TO WRITE
C
C   OUTPUT ARGUMENTS:
C     KA           INTEGER NUMBER OF BYTES IN FORTRAN RECORD
C                  (IN WHICH CASE THE NEXT FORTRAN RECORD
C                  SHOULD HAVE A START BYTE OF IB+KA),
C                  OR KA=-1 FOR I/O ERROR
C
C SUBPROGRAMS CALLED:
C   BAWRITEL       BYTE-ADDRESSABLE WRITE
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C
C$$$

!
      IMPLICIT NONE
      INTEGER,INTENT(IN):: LU,IB,NB
      INTEGER,INTENT(OUT):: KA
      CHARACTER,INTENT(IN):: A(NB)
      INTEGER(KIND=8) :: LONG_IB,LONG_NB,LONG_KA
!
        if((IB<0.and.IB/=-1) .or. NB<0 ) THEN
          print *,'WRONG: in BAFRWRITE starting postion IB or read '//   &
     &   'data size NB <0, STOP! ' //                                    &
     &   'Consider use BAFRRWRITEL and long integer'
          KA=0
          return
        ENDIF
        LONG_IB=IB
        LONG_NB=NB
        CALL BAFRWRITEL(LU,LONG_IB,LONG_NB,LONG_KA,A)
        KA=LONG_KA
!
      END SUBROUTINE BAFRWRITE
C-----------------------------------------------------------------------
      SUBROUTINE BAFRWRITEL(LU,IB,NB,KA,A)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: BAFRWRITEL     BYTE-ADDRESSABLE FORTRAN RECORD WRITE
C   PRGMMR: IREDELL          ORG: W/NMC23     DATE: 1999-01-21
C
C ABSTRACT: THIS SUBPROGRAM WRITES AN UNFORMATTED FORTRAN RECORD
C
C PROGRAM HISTORY LOG:
C   1999-01-21  IREDELL
C   2009-04-20  J. WANG
C
C USAGE:    CALL BAFRWRITEL(LU,IB,NB,KA,A)
C   INPUT ARGUMENTS:
C     LU           INTEGER LOGICAL UNIT TO WRITE
C     IB           INTEGER(8) FORTRAN RECORD START BYTE
C                  (FOR THE FIRST FORTRAN RECORD, IB SHOULD BE 0)
C     NB           INTEGER(8) NUMBER OF BYTES TO WRITE
C     A            CHARACTER*1 (NB) DATA TO WRITE
C
C   OUTPUT ARGUMENTS:
C     KA           INTEGER(8) NUMBER OF BYTES IN FORTRAN RECORD
C                  (IN WHICH CASE THE NEXT FORTRAN RECORD
C                  SHOULD HAVE A START BYTE OF IB+KA),
C                  OR KA=-1 FOR I/O ERROR
C
C SUBPROGRAMS CALLED:
C   BAWRITEL       BYTE-ADDRESSABLE WRITE
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C
C$$$
      IMPLICIT NONE
      INTEGER,INTENT(IN):: LU
      INTEGER(KIND=8),INTENT(IN):: IB,NB
      INTEGER(kind=8),INTENT(OUT):: KA
      CHARACTER,INTENT(IN):: A(NB)
!
      INTEGER(kind=8),PARAMETER:: LBCW=4
      INTEGER(kind=LBCW):: BCW
      INTEGER(kind=8):: KR
      INTEGER(LBCW):: BCW2,LBCW2
      CHARACTER(16) :: MACHINE_ENDIAN
      LOGICAL :: DO_BYTESWAP
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  WRITE DATA BRACKETED BY BLOCK CONTROL WORDS
!
!-- set do_byteswap from machine endianness and file endianness
      CALL CHK_ENDIANC(MACHINE_ENDIAN)
      IF( LU<=999) THEN
        IF( trim(MACHINE_ENDIAN)=="big_endian") THEN
          DO_BYTESWAP=.false.
        ELSEIF( trim(MACHINE_ENDIAN)=="little_endian") THEN
          DO_BYTESWAP=.true.
        ENDIF
      ELSEIF(LU<=1999) THEN
        IF( trim(MACHINE_ENDIAN)=="big_endian") THEN
          DO_BYTESWAP=.true.
        ELSEIF( trim(MACHINE_ENDIAN)=="little_endian") THEN
          DO_BYTESWAP=.false.
        ENDIF
      ENDIF
!
!
      BCW=NB
      IF(DO_BYTESWAP) CALL Byteswap(BCW,LBCW,1)
      CALL BAWRITEL(LU,IB,LBCW,KR,BCW)
      IF(KR.NE.LBCW) THEN
        KA=-1
      ELSE
        CALL BAWRITEL(LU,IB+LBCW,NB,KR,A)
        IF(KR.NE.NB) THEN
          KA=-1
        ELSE
          CALL BAWRITEL(LU,IB+LBCW+NB,LBCW,KR,BCW)
          IF(KR.NE.LBCW) THEN
            KA=-1
          ELSE
            KA=LBCW+NB+LBCW
          ENDIF
        ENDIF
      ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      END SUBROUTINE  BAFRWRITEL

max_stars_repo_name: GMAO-SI-Team/jedi-stack

max_issues_repo_name: GMAO-SI-Team/jedi-stack


program main
    use sim
    implicit none

    type(FEMDomain_),target :: domain
    type(MaterialProp_) :: YoungsModulus, PoissonRatio,Density,Cohesion,Phi,psi
    type(Boundary_) :: disp_x,disp_y,disp_z
    type(IO_) :: f
    type(FiniteDeform_) :: Solver
    integer(int32) :: i

    call domain%create(meshtype="rectangular3D",x_num=10,y_num=10,z_num=10,&
        x_len=50.0d0,y_len=20.0d0,z_len=10.0d0)

    !call domain%create(meshtype="Sphere3D",x_num=10,y_num=10,z_num=10,&
    !    x_len=50.0d0,y_len=20.0d0,z_len=10.0d0)
    !call domain%gmsh(Name="test")

    call YoungsModulus%create(Name="YoungsModulus",ParaValue=100.0d0,Layer=1)
    call PoissonRatio%create(Name="PoissonRatio",ParaValue=0.20d0,Layer=2)
    call Density%create(Name="Density",ParaValue=0.00d0,Layer=3)
    call Cohesion%create(Name="Cohesion",ParaValue=100000000.0d0,Layer=4)
    call Phi%create(Name="Psi",ParaValue=0.0d0,Layer=5)
    call psi%create(Name="psi",ParaValue=0.0d0,Layer=6)

    call disp_x%create(Category="Dirichlet",BoundValue=0.0d0,x_max=1.0d0,x_min=-1.0d0,Layer=1)!x
    call disp_y%create(Category="Dirichlet",BoundValue=0.0d0,x_max=1.0d0,x_min=-1.0d0,Layer=2)!y
    call disp_z%create(Category="Dirichlet",BoundValue=0.0d0,x_max=1.0d0,x_min=-1.0d0,Layer=3)!z
    call disp_x%create(Category="Dirichlet",BoundValue=5.0d0,x_max=51.0d0,x_min=40.0d0,Layer=1)

    call domain%import(Materials=.true., Material=YoungsModulus)
    call domain%import(Materials=.true., Material=PoissonRatio)
    call domain%import(Materials=.true., Material=Density)
    call domain%import(Materials=.true., Material=Cohesion)
    call domain%import(Materials=.true., Material=Phi)
    call domain%import(Materials=.true., Material=Psi)

    call domain%import(Boundaries=.true. , Boundary=disp_x)
    call domain%import(Boundaries=.true. , Boundary=disp_y)
    call domain%import(Boundaries=.true. , Boundary=disp_z)

    call domain%bake(template="FiniteDeform_")

    solver%FEMDomain => domain

    solver%dt=1.0d0
    solver%infinitesimal = .true.

    call solver%DivideBC()
    call solver%solve(solvertype="BiCGSTAB")
    do i=1,10
        call solver%UpdateInitConfig()
        call solver%UpdateBC()
        call solver%solve(solvertype="BiCGSTAB")
        call solver%display(DisplayMode="gmsh",OptionalStep=i,Name="deform")
    enddo

end program main


text:
      subroutine pzeroi(a,n)
      implicit none
      integer n, a(*), i
      do i=1, n
        a(i) = 0
      enddo
      return
      end

      subroutine pzero(a,n)
      implicit none
      integer          n, i
      double precision a(*), r0
      r0 = 0.d0
      do i=1, n
        a(i) = r0
      enddo
      return
      end

      subroutine pmovei(a,b,n)
      implicit none
      integer a(*), b(*), n, i
      do i=1, n
        b(i) = a(i)
      enddo
      return
      end

      subroutine pmove(a,b,n)
      implicit none
      integer          n, i
      double precision a(*), b(*)
      do i=1, n
        b(i) = a(i)
      enddo
      return
      end

      subroutine paddi(a,b,n)
      implicit none
      integer a(*), b(*), n, i
      do i=1, n
        b(i) = b(i) + a(i)
      enddo
      return
      end

      subroutine padd(a,b,n)
      implicit none
      integer          n, i
      double precision a(*), b(*)
      do i=1, n
        b(i) = b(i) + a(i)
      enddo
      return
      end

      subroutine psubi(a,b,n)
      implicit none
      integer a(*), b(*), n, i
      do i=1, n
        b(i) = b(i) - a(i)
      enddo
      return
      end

      subroutine psub(a,b,n)
      implicit none
      integer          n, i
      double precision a(*), b(*)
      do i=1, n
        b(i) = b(i) - a(i)
      enddo
      return
      end

ext: for

max_stars_repo_name: chennachaos/stabfem

max_issues_repo_name: chennachaos/stabfem

max_forks_repo_name: chennachaos/mpap


text:
module timeutils
use phys_consts, only: wp
use, intrinsic :: iso_fortran_env, only: sp=>real32, dp=>real64, int32, int64

implicit none (type, external)
private
public :: ymd2doy, sza, dateinc, utsec2filestem, date_filename, day_wrap, find_lastdate, find_time_elapsed

real(wp), parameter :: pi = 4._wp*atan(1._wp)

contains

elemental integer function daysmonth(year, month) result(days)

integer, intent(in) :: year, month
integer :: monthdays(12)

if ((year < 1600) .or. (year > 2500)) error stop 'is year specified correctly?'
if ((month < 1) .or. (month > 12)) error stop 'impossible month'

monthdays = [31,28,31,30,31,30,31,31,30,31,30,31]

if (mod(year, 4)==0 .and. mod(year, 100)/=0 .or. mod(year, 400) == 0) monthdays(2)=29

days = monthdays(month)

end function daysmonth

elemental integer function ymd2doy(year, month, day) result(doy)

integer, intent(in) :: year, month, day
integer :: i

if ((day < 1) .or. (day > daysmonth(year, month))) error stop 'impossible day'

doy = 0
do i = 1, month-1
  doy = doy + daysmonth(year, i)
enddo

doy = doy + day

end function ymd2doy

elemental function sza(year, month, day, UTsec,glat,glon)
!! computes sza in radians
!! CALCULATE SOLAR ZENITH ANGLE OVER A GIVEN GET OF LAT/LON

integer, intent(in) :: year, month, day
real(wp), intent(in) :: UTsec
real(wp), intent(in) :: glat,glon
real(wp) :: sza

real(wp), parameter :: tau = 2._wp*pi

real(wp) :: doy,soldecrad
real(wp) :: lonrad,LThrs,latrad,hrang

!> SOLAR DECLINATION ANGLE
doy = ymd2doy(year, month, day)
soldecrad = -23.44_wp*cos(tau/365._wp*(doy+10)) * pi/180

!> HOUR ANGLE
lonrad=glon*pi/180
lonrad = modulo(lonrad, 2*pi)

LThrs=UTsec/3600._wp+lonrad/(pi/12._wp)
hrang=(12-LThrs)*(pi/12._wp)

!> SOLAR ZENITH ANGLE
latrad=glat*pi/180._wp
sza=acos(sin(soldecrad)*sin(latrad)+cos(soldecrad)*cos(latrad)*cos(hrang))

end function sza

pure subroutine dateinc(dtsec, ymd, UTsec)
!! increment datetime by dtsec

real(wp), intent(in) :: dtsec
!! seconds to increment
integer, intent(inout) :: ymd(3)
!! year, month, day of month
real(wp), intent(inout) :: UTsec
!! seconds since midnight UTC

integer :: year,month,day
integer :: monthinc          !< we incremented the month

year=ymd(1); month=ymd(2); day=ymd(3);

if (day < 1) error stop 'temporal:timeutils:dateinc(): days are positive integers'
if (utsec < 0) error stop 'negative UTsec, simulation should go forward in time only!'
if (dtsec < 0) error stop 'negative dtsec, simulation should go forward in time only!'
if (dtsec > 86400) error stop 'excessively large dtsec > 86400, simulation step should be small enough!'

UTsec = UTsec + dtsec
do while (UTsec >= 86400)
  UTsec = UTsec - 86400._wp
  day = day+1
  call day_wrap(year, month, day)
end do

ymd(1)=year; ymd(2)=month; ymd(3)=day;    !replace input array with new date

end subroutine dateinc

recursive pure subroutine day_wrap(year, month, day)
!! increment date if needed, according to day
!! that is, if day is beyond month, increment month and year if needed
integer, intent(inout) :: year, month, day

if (month < 1 .or. day < 1) error stop 'day_wrap: months and days must be positive'

!> wrap months
do while (month > 12)
  month = month - 12
  year = year + 1
end do

!> wrap days
do while (day > daysmonth(year, month))
  day = day - daysmonth(year, month)
  month = month + 1
  call day_wrap(year, month, day)
end do

end subroutine day_wrap

pure function date_filename(outdir, ymd, UTsec)
!! GENERATE A FILENAME stem OUT OF A GIVEN DATE/TIME
!! (does not include suffix like .h5)

character(*), intent(in) :: outdir
integer, intent(in) :: ymd(3)
class(*), intent(in) :: UTsec
character(:), allocatable :: date_filename

character(len=21) :: stem

stem = utsec2filestem(ymd, UTsec)

date_filename = outdir // '/' // stem

end function date_filename

pure character(21) function utsec2filestem(ymd, UTsec)
!! file stem is exactly 21 characters long, per Matt Z's de facto spec.
!! we keep microsecond dummy precision filenames to be legacy compatible
!! true filename resolution is 10 milliseconds due to real32 7 digits of precision vis 86400 second day.
integer, intent(in) :: ymd(3)
class(*), intent(in) :: UTsec
!! UTC second: real [0.0 .. 86400.0)

character(12) :: sec_str
integer :: year, month, day, seconds, millisec, frac

year = ymd(1)
month = ymd(2)
day = ymd(3)

select type(UTsec)
  type is (real(sp))
    !! round to nearest ten milliseconds
    millisec = nint(UTsec * 100) * 10
  type is (real(dp))
    !! round to nearest ten milliseconds
    millisec = nint(UTsec * 100) * 10
  type is (integer(int32))
    millisec = UTsec * 1000
  type is (integer(int64))
    millisec = int(UTsec) * 1000
  class default
    error stop "timeutils.f90:utsec2filestem unknown UTsec type"
end select

seconds = millisec / 1000 !< truncate fractional second
if (seconds < 0 .or. seconds > 86400) error stop 'timeutils.f90::utsec2filestem did NOT satisfy 0 <= seconds < 86400'
if (seconds == 86400) then
  !> FIXME: This corner case is from not using integers for microseconds
  ! write(stderr,*) 'utsec2filestem: FIXME: patching UTsec=86400 to next day midnight'
  day = day+1
  seconds = 0
  millisec = 0
  call day_wrap(year, month, day)
endif

frac = millisec*1000 - seconds * 1000000  !< microseconds

write(sec_str, '(I5.5, A1, I6.6)') seconds, '.', frac

write(utsec2filestem, '(i4,I2.2,I2.2,a13)') year, month, day, '_' // sec_str

end function utsec2filestem

pure subroutine find_lastdate(ymd0,UTsec0,ymdtarget,UTsectarget,cadence,ymd,UTsec)

!> Compute the last date before the target date based on a start date and date cadence.  The
!  file is assumed to exist and programmer needs to check for existence outside this
!  procedure.
! FIXME:  this will occasionally hang and may need to be re-examined at some point...

integer, dimension(3), intent(in) :: ymd0
real(wp), intent(in) :: UTsec0
integer, dimension(3), intent(in) :: ymdtarget
real(wp), intent(in) :: UTsectarget
real(wp), intent(in) :: cadence
integer, dimension(3), intent(out) :: ymd
real(wp), intent(out) :: UTsec

integer, dimension(3) :: ymdnext
real(wp) :: UTsecnext
logical :: flagend

ymd=ymd0
UTsec=UTsec0
ymdnext=ymd0
UTsecnext=UTsec0
flagend=ymdnext(1)>=ymdtarget(1) .and. ymdnext(2)>=ymdtarget(2) .and. ymdnext(3)>=ymdtarget(3) .and. UTsecnext>UTsectarget &
          .or. ymdnext(1)>ymdtarget(1) .or. ymdnext(2)>ymdtarget(2) .or. ymdnext(3)>ymdtarget(3) ! in case the first time step is the last before target
do while ( .not.(flagend) )
  ymd=ymdnext
  UTsec=UTsecnext
  call dateinc(cadence,ymdnext,UTsecnext)
  flagend=ymdnext(1)>=ymdtarget(1) .and. ymdnext(2)>=ymdtarget(2) .and. ymdnext(3)>=ymdtarget(3) .and. UTsecnext>UTsectarget &
            .or. ymdnext(1)>ymdtarget(1) .or. ymdnext(2)>ymdtarget(2) .or. ymdnext(3)>ymdtarget(3)
end do
! When the loops exits ymd,UTsec have the date of the last output file before the given target time OR the first output file in the case that the target date is before the begin date...

end subroutine find_lastdate

pure function find_time_elapsed(ymdstart,UTsecstart,ymdend,UTsecend,dt) result(telapsed)

! finds the amount of time that has elapsed between a given start and end date, using a given dt increment
! the resulting elapsed time will be the smallest multiple of dt that is >= true elapsed time

!! inputs
integer, dimension(3), intent(in) :: ymdstart,ymdend
real(wp), intent(in) :: UTsecstart,UTsecend
reaL(wp), intent(in) :: dt

integer, dimension(3) :: ymdnow
real(wp) :: UTsecnow

real(wp) :: telapsed  !! output

telapsed=0._wp; ymdnow=ymdstart; UTsecnow=UTsecstart;
do while (.not. (all(ymdend==ymdnow) .and. UTsecnow>UTsecend) )
  call dateinc(dt,ymdnow,UTsecnow)
  telapsed=telapsed+dt
end do
telapsed=telapsed-dt

end function find_time_elapsed

end module timeutils

max_forks_repo_name: mattzett/GEMINI


text:
great_rsd(a1,m1).
great_rsd(w1,hh1).
great_rsd(z1,dd1).
great_rsd(g1,h1).
great_rsd(a1,w1).
great_rsd(u1,aa1).
great_rsd(e1,t1).
great_rsd(q1,h1).
great_rsd(dd1,x1).
great_rsd(k1,h1).
great_rsd(w1,aa1).
great_rsd(w1,h1).
great_rsd(p1,ff1).
great_rsd(w1,t1).
great_rsd(e1,ff1).
great_rsd(w1,g1).
great_rsd(v1,z1).
great_rsd(q1,c1).
great_rsd(cc1,q1).
great_rsd(p1,cc1).
great_rsd(e1,u1).
great_rsd(f1,t1).
great_rsd(ff1,ii1).
great_rsd(a1,aa1).
great_rsd(c1,dd1).
great_rsd(f1,q1).
great_rsd(v1,hh1).
great_rsd(w1,q1).
great_rsd(t1,kk1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
c234567890

      integer i,j,k,l,tmpline,tmpint,tmpint1,tmpint2,ttt,tta,ttb,ttu
      integer countat,begin,endd,countchgrp,pdbnCGbeads
      integer atom1,atom2,atom3,atom4
      real tmpdist,tmpdistx,tmpdisty,tmpdistz

c DATA for the pdb file
      integer pdbnres,pdbnatom,pdbatnumb(10000),pdbresnum(10000)
      real ax(10000),ay(10000),az(10000)
      character*4 tmp
      character*3 pdbatname(10000),pdbresname(10000)
c CG FF data
      character*3 resname(20)
      integer nbeads(20),nbond(20),nangle(20),ndihe(20),nconst(20)
      character*3 tmpchar,atomname(20,7)
      character*3 beadtype(20,7)
      integer chargrp(20,7)
      real charbead(20,7)
      integer bondatom1(20,10),bondatom2(20,10),bondtype(20,10)
      real bondist(20,10),bondk(20,10)
      integer constatom1(20,10),constatom2(20,10),consttype(20,10)
      real constdist(20,10)
      integer angleatom1(20,10),angleatom2(20,10),angleatom3(20,10)
      integer angletype(20,10)
      real angle(20,10),anglek(20,10)
      integer diheatom1(20,10),diheatom2(20,10),diheatom3(20,10)
      integer diheatom4(20,10),dihetype(20,10)
      real dihe(20,10),dihek(20,10)
c DATA for topology (top)
      character*3 topatname(10000),tpresname(1000),CGname(1000)
      real atcharge(1000,7)
      integer tprestype(1000),tpnbeads(1000),beginres(1000)
      character*3 tpbeadtype(1000,7),tpatomname(1000,7)
      integer tpatnumb(1000,7),tpchrgrp(1000,7)
      real tpatcharge(1000,7)
      real*8 vector1(3),vector2(3),angleCACACA,angleCACACB,angleCBCACA
c AA FF data
      integer tmp2,natomresAA(20),pdbnatomAA,pdbnresAA
      integer natomresCG(20)
      character*3 nameresAA(20)
      character*120 tmpbigchar
c Stuff for the COMs
      integer countatcom
      integer countatCG,countatAA,countatSC
      real tmpcomx,tmpcomy,tmpcomz
      real tmpcom1x,tmpcom1y,tmpcom1z
      real tmpcom2x,tmpcom2y,tmpcom2z
      real tmpcom3x,tmpcom3y,tmpcom3z
      real tmpcom4x,tmpcom4y,tmpcom4z
      real comBBx(1000),comBBy(1000),comBBz(1000)
      real comSCx(1000,4),comSCy(1000,4),comSCz(1000,4)
c Paramters
      parameter     (pi=3.141592654)

C--------------------------------------------------------------
C This program is ment to read a pdb file and generate a new
C pdb file corresponding to the topology of FFCG-2.1.
C The main idea is to use the pdb structure to contruct the
C beads of the CG model.
C All comments to x.periole@rug.nl; NYC-CSICUNY-03-05-2007
C--------------------------------------------------------------
c  Important Notes for Correct Functioning:
c  0- OF COURSE the use of this program is at your own risks
c  1- the pdb file should be cleaned of hydrogens, if any,
c     HEADER/REMARK etc, and Nterm and Cterm atoms.
c  3- the first line of the pdb file must contain the number
c     of residues and the number of atoms.
c  4- The pdb file should not contain a CHAIN-ID (between the
c     resname and the resnumber)
C--------------------------------------------------------------

C  Define stuff for the CG model

C  Integrate the info of the CG FF
c  Note : the resname are red in cg.dat and then they are
c  refered as the index as the order of appearance
C

      open(unit=11,name="cg-2.1.dat",status='old')
c     read(11,*) tmpline

      do i=1,20
      read(11,222) resname(i),nbeads(i),nbond(i),nangle(i),ndihe(i)
     .,nconst(i)
c     write(6,222) resname(i),nbeads(i),nbond(i),nangle(i),ndihe(i)
c    .,nconst(i)
222   format(A3,4X,I1,3X,I1,3X,I1,3x,I1,3x,I1)
      do j=1,nbeads(i)
      read(11,*) tmpint1,beadtype(i,j),tmpint2,tmpchar,atomname(i,j)
     .,chargrp(i,j),charbead(i,j)
c     write(6,*) tmpint1,beadtype(i,j),tmpint2,tmpchar,atomname(i,j)
c    .,chargrp(i,j),charbead(i,j)
      enddo
      do j=1,nbond(i)
      read(11,*) bondatom1(i,j),bondatom2(i,j),bondtype(i,j)
     .,bondist(i,j),bondk(i,j)
c     write(6,*) bondatom1(i,j),bondatom2(i,j),bondtype(i,j)
c    ,,bondist(i,j),bondk(i,j)
      enddo
      do j=1,nconst(i)
      read(11,*) constatom1(i,j),constatom2(i,j),consttype(i,j)
     .,constdist(i,j)
c     write(6,*) constatom1(i,j),constatom2(i,j),consttype(i,j)
c    .,constdist(i,j)
      enddo
      do j=1,nangle(i)
      read(11,*) angleatom1(i,j),angleatom2(i,j),angleatom3(i,j)
     .,angletype(i,j),angle(i,j),anglek(i,j)
c     write(6,*) angleatom1(i,j),angleatom2(i,j),angleatom3(i,j)
c    .,angletype(i,j),angle(i,j),anglek(i,j)
      enddo
      do j=1,ndihe(i)
      read(11,*) diheatom1(i,j),diheatom2(i,j),diheatom3(i,j)
     .,diheatom4(i,j),dihetype(i,j),dihe(i,j),dihek(i,j)
c     write(6,*) diheatom1(i,j),diheatom2(i,j),diheatom3(i,j)
c    .,diheatom4(i,j),dihetype(i,j),dihe(i,j),dihek(i,j)
      enddo
      enddo

C  Read the info concerning residues

      open(unit=22,name="AA-2.1.dat",status='old')

      read (22,*) tmpbigchar
      read (22,*) tmpbigchar
      read (22,*) tmpbigchar
      read (22,*) tmp2
      do i = 1, tmp2
      read (22,*) nameresAA(i),natomresAA(i),natomresCG(i)
c     write (6,*) nameresAA(i),natomresAA(i),natomresCG(i)
      enddo

C  Read the pdb file -

c first read the number of residues and atoms

      read (5,*)  pdbnresAA,pdbnatomAA
      if (pdbnresAA .le. 0 .or. pdbnatomAA .le. 0 ) then
      write(6,*) "The PDB file you give as an input is not formatted
     .correctly"
      write(6,*) "Tip: The first line should be the number of residues
     .and atoms."
      stop
      endif

      do i = 1,pdbnatomAA
c     read (5,*) tmp,pdbatnumb(i),pdbatname(i),pdbresname(i),pdbresnum(i),ax(i),
c    .ay(i),az(i)
      read (5,662) pdbatnumb(i),pdbatname(i),pdbresname(i),pdbresnum(i),
     .ax(i),ay(i),az(i)
c     write (6,662) pdbatnumb(i),pdbatname(i),pdbresname(i),pdbresnum(i)
c    .,ax(i),ay(i),az(i)
      enddo

C  getting the info for each residue

c  here the counter of beads for the CG
      pdbnCGbeads=0
      countat = 1
      do i = 1,pdbnresAA

c define the residue type

       beginres(i)=countat
       tpresname(i)=pdbresname(beginres(i))

        do j =1,20
        if (tpresname(i) .eq. resname(j)) then
        tprestype(i) = j
c       write(6,*) i,tpresname(i),resname(j),tprestype(i),countat
        endif
        enddo

        if ((tprestype(i) .lt. 1) .and. (tprestype(i) .gt. 20)) then
        write (6,*) "Something is wrong with the attribution of the
     .residue index corresponding to its name in the pdb file!"
        stop
        endif

c here adding the number of beads for the current residue

        pdbnCGbeads=pdbnCGbeads+nbeads(tprestype(i))

c initialize current COM for CA

        tmpcomx = 0.0
        tmpcomy = 0.0
        tmpcomz = 0.0
        countatcom = 0.0

c get the COM  for BB

        do k=1,natomresAA(tprestype(i))
c        write(6,663) pdbatname(beginres(i)+k-1)
         if((pdbatname(beginres(i)+k-1) .eq. "CA ")) then
c        write(6,663) pdbatname(beginres(i)+k-1)
         tmpcomx = ax(beginres(i)+k-1)
         tmpcomy = ay(beginres(i)+k-1)
         tmpcomz = az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
         countat = countat + 1
        enddo
663   format(A4)

c calculate the COM for CA ; note the unusual way to calculate(old stuff)

         if(countatcom .eq. 1) then
         comBBx(i) = tmpcomx
         comBBy(i) = tmpcomy
         comBBz(i) = tmpcomz
         endif

C Sidechains
c initialize the COM for the SC

        tmpcomx = 0.0
        tmpcomy = 0.0
        tmpcomz = 0.0
        tmpcom1x = 0.0
        tmpcom1y = 0.0
        tmpcom1z = 0.0
        tmpcom2x = 0.0
        tmpcom2y = 0.0
        tmpcom2z = 0.0
        tmpcom3x = 0.0
        tmpcom3y = 0.0
        tmpcom3z = 0.0
        tmpcom4x = 0.0
        tmpcom4y = 0.0
        tmpcom4z = 0.0
        countatcom = 0.0
        countattmp = 0.0

c GLY; ALA no COM for SC

c CYS
       if(tprestype(i) .eq. 3) then
        do k=1,natomresAA(tprestype(i))
         if(pdbatname(beginres(i)+k-1) .eq. "CB "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "SG "  ) then
         tmpcomx = tmpcomx + ax(beginres(i)+k-1)
         tmpcomy = tmpcomy + ay(beginres(i)+k-1)
         tmpcomz = tmpcomz + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
        enddo

c calculate the COM for SC

         if(countatcom .eq. 2) then
         comSCx(i,1) = tmpcomx/2.0
         comSCy(i,1) = tmpcomy/2.0
         comSCz(i,1) = tmpcomz/2.0
         else
         write(6,*) "something wrong in the Cys ",i
         stop
         endif
         countattmp = countattmp + 1
       endif

c VAL
       if(tprestype(i) .eq. 4) then
        do k=1,natomresAA(tprestype(i))
         if(pdbatname(beginres(i)+k-1) .eq. "CB "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CG1"  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CG2"  ) then
         tmpcomx = tmpcomx + ax(beginres(i)+k-1)
         tmpcomy = tmpcomy + ay(beginres(i)+k-1)
         tmpcomz = tmpcomz + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
         countattmp = countattmp + 1
        enddo

c calculate the COM for SC

         if(countatcom .eq. 3) then
         comSCx(i,1) = tmpcomx/3.0
         comSCy(i,1) = tmpcomy/3.0
         comSCz(i,1) = tmpcomz/3.0
         else
         write(6,*) "something wrong in the Val ",i
         stop
         endif
       endif

c LEU

       if(tprestype(i) .eq. 5) then
        do k=1,natomresAA(tprestype(i))
         if(pdbatname(beginres(i)+k-1) .eq. "CB "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CG "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CD1"  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CD2"  ) then
         tmpcomx = tmpcomx + ax(beginres(i)+k-1)
         tmpcomy = tmpcomy + ay(beginres(i)+k-1)
         tmpcomz = tmpcomz + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
         countattmp = countattmp + 1
        enddo

c calculate the COM for SC

         if(countatcom .eq. 4) then
         comSCx(i,1) = tmpcomx/4.0
         comSCy(i,1) = tmpcomy/4.0
         comSCz(i,1) = tmpcomz/4.0
         else
         write(6,*) "something wrong in the Leu ",i
         stop
         endif
       endif

c ILE
       if(tprestype(i) .eq. 6) then
        do k=1,natomresAA(tprestype(i))
c        write(6,*) pdbatname(beginres(i)+k-1)
         if(pdbatname(beginres(i)+k-1) .eq. "CB "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CG1"  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CG2"  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CD1"  ) then
c        write(6,*) pdbatname(beginres(i)+k-1)
         tmpcomx = tmpcomx + ax(beginres(i)+k-1)
         tmpcomy = tmpcomy + ay(beginres(i)+k-1)
         tmpcomz = tmpcomz + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
         countattmp = countattmp + 1
        enddo

c calculate the COM for SC

         if(countatcom .eq. 4) then
         comSCx(i,1) = tmpcomx/4.0
         comSCy(i,1) = tmpcomy/4.0
         comSCz(i,1) = tmpcomz/4.0
         else
         write(6,*) "something wrong in the Ile ",i
         stop
         endif
       endif

c MET
       if(tprestype(i) .eq. 7) then
        do k=1,natomresAA(tprestype(i))
c        write(6,*) pdbatname(beginres(i)+k-1)
c        write(6,*) tmpcomx,tmpcomy,tmpcomz
         if(pdbatname(beginres(i)+k-1) .eq. "CB "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CG "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "SD "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CE "  ) then
         tmpcomx = tmpcomx + ax(beginres(i)+k-1)
         tmpcomy = tmpcomy + ay(beginres(i)+k-1)
         tmpcomz = tmpcomz + az(beginres(i)+k-1)
c        write(6,*) tmpcomx,tmpcomy,tmpcomz
         countatcom = countatcom+1
         endif
         countattmp = countattmp + 1
        enddo

c calculate the COM for SC

         if(countatcom .eq. 4) then
         comSCx(i,1) = tmpcomx/4.0
         comSCy(i,1) = tmpcomy/4.0
         comSCz(i,1) = tmpcomz/4.0
         else
         write(6,*) "something wrong in the Met ",i
         stop
         endif
       endif

c PRO
       if(tprestype(i) .eq. 8) then
        do k=1,natomresAA(tprestype(i))
         if(pdbatname(beginres(i)+k-1) .eq. "CB "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CG "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CD "  ) then
         tmpcomx = tmpcomx + ax(beginres(i)+k-1)
         tmpcomy = tmpcomy + ay(beginres(i)+k-1)
         tmpcomz = tmpcomz + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
         countattmp = countattmp + 1
        enddo

c calculate the COM for SC

         if(countatcom .eq. 3) then
         comSCx(i,1) = tmpcomx/3.0
         comSCy(i,1) = tmpcomy/3.0
         comSCz(i,1) = tmpcomz/3.0
         else
         write(6,*) "something wrong in the Pro ",i
         stop
         endif
       endif

c ASN
       if(tprestype(i) .eq. 9) then
        do k=1,natomresAA(tprestype(i))
         if(pdbatname(beginres(i)+k-1) .eq. "CB "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CG "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "OD1"  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "ND2"  ) then
         tmpcomx = tmpcomx + ax(beginres(i)+k-1)
         tmpcomy = tmpcomy + ay(beginres(i)+k-1)
         tmpcomz = tmpcomz + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
         countattmp = countattmp + 1
        enddo

c calculate the COM for SC

         if(countatcom .eq. 4) then
         comSCx(i,1) = tmpcomx/4.0
         comSCy(i,1) = tmpcomy/4.0
         comSCz(i,1) = tmpcomz/4.0
         else
         write(6,*) "something wrong in the Asn ",i
         stop
         endif
       endif

c GLN
       if(tprestype(i) .eq. 10) then
        do k=1,natomresAA(tprestype(i))
         if(pdbatname(beginres(i)+k-1) .eq. "CB "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CG "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CD "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "OE1"  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "NE2"  ) then
         tmpcomx = tmpcomx + ax(beginres(i)+k-1)
         tmpcomy = tmpcomy + ay(beginres(i)+k-1)
         tmpcomz = tmpcomz + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
         countattmp = countattmp + 1
        enddo

c calculate the COM for SC

         if(countatcom .eq. 5) then
         comSCx(i,1) = tmpcomx/5.0
         comSCy(i,1) = tmpcomy/5.0
         comSCz(i,1) = tmpcomz/5.0
         else
         write(6,*) "something wrong in the Gln ",i
         stop
         endif
       endif

c LYS ; here two beads!
       if(tprestype(i) .eq. 11) then
        do k=1,natomresAA(tprestype(i))
c COM1
         if(pdbatname(beginres(i)+k-1) .eq. "CB "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CG "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CD "  ) then
         tmpcom1x = tmpcom1x + ax(beginres(i)+k-1)
         tmpcom1y = tmpcom1y + ay(beginres(i)+k-1)
         tmpcom1z = tmpcom1z + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
c COM2
         if(pdbatname(beginres(i)+k-1) .eq. "CE "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "NZ "  ) then
         tmpcom2x = tmpcom2x + ax(beginres(i)+k-1)
         tmpcom2y = tmpcom2y + ay(beginres(i)+k-1)
         tmpcom2z = tmpcom2z + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
         countattmp = countattmp + 1
        enddo

c calculate the COM for SC

         if(countatcom .eq. 5) then
         comSCx(i,1) = tmpcom1x/3.0
         comSCy(i,1) = tmpcom1y/3.0
         comSCz(i,1) = tmpcom1z/3.0
         comSCx(i,2) = tmpcom2x/2.0
         comSCy(i,2) = tmpcom2y/2.0
         comSCz(i,2) = tmpcom2z/2.0
         else
         write(6,*) "something wrong in the Lys ",i
         stop
         endif
       endif

c ASP
       if(tprestype(i) .eq. 12) then
        do k=1,natomresAA(tprestype(i))
         if(pdbatname(beginres(i)+k-1) .eq. "CB "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CG "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "OD1"  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "OD2"  ) then
         tmpcomx = tmpcomx + ax(beginres(i)+k-1)
         tmpcomy = tmpcomy + ay(beginres(i)+k-1)
         tmpcomz = tmpcomz + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
         countattmp = countattmp + 1
        enddo

c calculate the COM for SC

         if(countatcom .eq. 4) then
         comSCx(i,1) = tmpcomx/4.0
         comSCy(i,1) = tmpcomy/4.0
         comSCz(i,1) = tmpcomz/4.0
         else
         write(6,*) "something wrong in the Asp ",i
         stop
         endif
       endif

c GLU
       if(tprestype(i) .eq. 13) then
        do k=1,natomresAA(tprestype(i))
         if(pdbatname(beginres(i)+k-1) .eq. "CB "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CG "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CD "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "OE1"  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "OE2"  ) then
         tmpcomx = tmpcomx + ax(beginres(i)+k-1)
         tmpcomy = tmpcomy + ay(beginres(i)+k-1)
         tmpcomz = tmpcomz + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
         countattmp = countattmp + 1
        enddo

c calculate the COM for SC

         if(countatcom .eq. 5) then
         comSCx(i,1) = tmpcomx/5.0
         comSCy(i,1) = tmpcomy/5.0
         comSCz(i,1) = tmpcomz/5.0
         else
         write(6,*) "something wrong in the Glu ",i
         stop
         endif
       endif

c THR
       if(tprestype(i) .eq. 14) then
        do k=1,natomresAA(tprestype(i))
         if(pdbatname(beginres(i)+k-1) .eq. "CB "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "OG1"  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CG2"  ) then
         tmpcomx = tmpcomx + ax(beginres(i)+k-1)
         tmpcomy = tmpcomy + ay(beginres(i)+k-1)
         tmpcomz = tmpcomz + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
         countattmp = countattmp + 1
        enddo

c calculate the COM for SC

         if(countatcom .eq. 3) then
         comSCx(i,1) = tmpcomx/3.0
         comSCy(i,1) = tmpcomy/3.0
         comSCz(i,1) = tmpcomz/3.0
         else
         write(6,*) "something wrong in the Thr ",i
         stop
         endif
       endif

c SER
       if(tprestype(i) .eq. 15) then
        do k=1,natomresAA(tprestype(i))
         if(pdbatname(beginres(i)+k-1) .eq. "CB "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "OG "  ) then
         tmpcomx = tmpcomx + ax(beginres(i)+k-1)
         tmpcomy = tmpcomy + ay(beginres(i)+k-1)
         tmpcomz = tmpcomz + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
         countattmp = countattmp + 1
        enddo

c calculate the COM for SC

         if(countatcom .eq. 2) then
         comSCx(i,1) = tmpcomx/2.0
         comSCy(i,1) = tmpcomy/2.0
         comSCz(i,1) = tmpcomz/2.0
         else
         write(6,*) "something wrong in the Ser ",i
         stop
         endif
       endif

c ARG ; here two beads!
       if(tprestype(i) .eq. 16) then
        do k=1,natomresAA(tprestype(i))
c COM1
         if(pdbatname(beginres(i)+k-1) .eq. "CB "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CG "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CD "  ) then
         tmpcom1x = tmpcom1x + ax(beginres(i)+k-1)
         tmpcom1y = tmpcom1y + ay(beginres(i)+k-1)
         tmpcom1z = tmpcom1z + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
c COM2
         if(pdbatname(beginres(i)+k-1) .eq. "NE "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CZ "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "NH1"  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "NH2"  ) then
         tmpcom2x = tmpcom2x + ax(beginres(i)+k-1)
         tmpcom2y = tmpcom2y + ay(beginres(i)+k-1)
         tmpcom2z = tmpcom2z + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
         countattmp = countattmp + 1
        enddo

c calculate the COM for SC

         if(countatcom .eq. 7) then
         comSCx(i,1) = tmpcom1x/3.0
         comSCy(i,1) = tmpcom1y/3.0
         comSCz(i,1) = tmpcom1z/3.0
         comSCx(i,2) = tmpcom2x/4.0
         comSCy(i,2) = tmpcom2y/4.0
         comSCz(i,2) = tmpcom2z/4.0
         else
         write(6,*) "something wrong in the Arg ",i
         stop
         endif
       endif

c HIS
       if(tprestype(i) .eq. 17) then
        do k=1,natomresAA(tprestype(i))
         if(pdbatname(beginres(i)+k-1) .eq. "CB "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "CG "  ) then
         tmpcom1x = tmpcom1x + ax(beginres(i)+k-1)
         tmpcom1y = tmpcom1y + ay(beginres(i)+k-1)
         tmpcom1z = tmpcom1z + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif

         if(pdbatname(beginres(i)+k-1) .eq. "ND1"       ) then
         tmpcom2x = tmpcom2x + ax(beginres(i)+k-1)
         tmpcom2y = tmpcom2y + ay(beginres(i)+k-1)
         tmpcom2z = tmpcom2z + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif

         if(pdbatname(beginres(i)+k-1) .eq. "NE2"       ) then
         tmpcom3x = tmpcom3x + ax(beginres(i)+k-1)
         tmpcom3y = tmpcom3y + ay(beginres(i)+k-1)
         tmpcom3z = tmpcom3z + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
         countattmp = countattmp + 1
        enddo

c calculate the COM for SC

         if(countatcom .eq. 4) then
         comSCx(i,1) = tmpcom1x/2.0
         comSCy(i,1) = tmpcom1y/2.0
         comSCz(i,1) = tmpcom1z/2.0
         comSCx(i,2) = tmpcom2x/1.0
         comSCy(i,2) = tmpcom2y/1.0
         comSCz(i,2) = tmpcom2z/1.0
         comSCx(i,3) = tmpcom3x/1.0
         comSCy(i,3) = tmpcom3y/1.0
         comSCz(i,3) = tmpcom3z/1.0
         else
         write(6,*) "something wrong in the His ",i
         stop
         endif
       endif

c PHE
       if(tprestype(i) .eq. 18) then
        do k=1,natomresAA(tprestype(i))
         if(pdbatname(beginres(i)+k-1) .eq. "CD1"       ) then
         tmpcom1x = tmpcom1x + ax(beginres(i)+k-1)
         tmpcom1y = tmpcom1y + ay(beginres(i)+k-1)
         tmpcom1z = tmpcom1z + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif

         if(pdbatname(beginres(i)+k-1) .eq. "CD2"       ) then
         tmpcom2x = tmpcom2x + ax(beginres(i)+k-1)
         tmpcom2y = tmpcom2y + ay(beginres(i)+k-1)
         tmpcom2z = tmpcom2z + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif

         if(pdbatname(beginres(i)+k-1) .eq. "CZ "       ) then
         tmpcom3x = tmpcom3x + ax(beginres(i)+k-1)
         tmpcom3y = tmpcom3y + ay(beginres(i)+k-1)
         tmpcom3z = tmpcom3z + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
         countattmp = countattmp + 1
        enddo

c calculate the COM for SC

         if(countatcom .eq. 3) then
         comSCx(i,1) = tmpcom1x/1.0
         comSCy(i,1) = tmpcom1y/1.0
         comSCz(i,1) = tmpcom1z/1.0
         comSCx(i,2) = tmpcom2x/1.0
         comSCy(i,2) = tmpcom2y/1.0
         comSCz(i,2) = tmpcom2z/1.0
         comSCx(i,3) = tmpcom3x/1.0
         comSCy(i,3) = tmpcom3y/1.0
         comSCz(i,3) = tmpcom3z/1.0
         else
         write(6,*) "something wrong in the Phe ",i
         stop
         endif
       endif

c TYR
       if(tprestype(i) .eq. 19) then
        do k=1,natomresAA(tprestype(i))
         if(pdbatname(beginres(i)+k-1) .eq. "CD1"       ) then
         tmpcom1x = tmpcom1x + ax(beginres(i)+k-1)
         tmpcom1y = tmpcom1y + ay(beginres(i)+k-1)
         tmpcom1z = tmpcom1z + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif

         if(pdbatname(beginres(i)+k-1) .eq. "CD2"       ) then
         tmpcom2x = tmpcom2x + ax(beginres(i)+k-1)
         tmpcom2y = tmpcom2y + ay(beginres(i)+k-1)
         tmpcom2z = tmpcom2z + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif

         if(pdbatname(beginres(i)+k-1) .eq. "CZ "  .or.
     .      pdbatname(beginres(i)+k-1) .eq. "OH "  ) then
         tmpcom3x = tmpcom3x + ax(beginres(i)+k-1)
         tmpcom3y = tmpcom3y + ay(beginres(i)+k-1)
         tmpcom3z = tmpcom3z + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
         countattmp = countattmp + 1
        enddo

c calculate the COM for SC

         if(countatcom .eq. 4) then
         comSCx(i,1) = tmpcom1x/1.0
         comSCy(i,1) = tmpcom1y/1.0
         comSCz(i,1) = tmpcom1z/1.0
         comSCx(i,2) = tmpcom2x/1.0
         comSCy(i,2) = tmpcom2y/1.0
         comSCz(i,2) = tmpcom2z/1.0
         comSCx(i,3) = tmpcom3x/2.0
         comSCy(i,3) = tmpcom3y/2.0
         comSCz(i,3) = tmpcom3z/2.0
         else
         write(6,*) "something wrong in the Tyr ",i
         stop
         endif
       endif

c TRP
       if(tprestype(i) .eq. 20) then
        do k=1,natomresAA(tprestype(i))
         if(pdbatname(beginres(i)+k-1) .eq. "CG "       ) then
         tmpcom1x = tmpcom1x + ax(beginres(i)+k-1)
         tmpcom1y = tmpcom1y + ay(beginres(i)+k-1)
         tmpcom1z = tmpcom1z + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif

         if(pdbatname(beginres(i)+k-1) .eq. "NE1"       ) then
         tmpcom2x = tmpcom2x + ax(beginres(i)+k-1)
         tmpcom2y = tmpcom2y + ay(beginres(i)+k-1)
         tmpcom2z = tmpcom2z + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif

         if(pdbatname(beginres(i)+k-1) .eq. "CZ2"       ) then
         tmpcom3x = tmpcom3x + ax(beginres(i)+k-1)
         tmpcom3y = tmpcom3y + ay(beginres(i)+k-1)
         tmpcom3z = tmpcom3z + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif

         if(pdbatname(beginres(i)+k-1) .eq. "CE3"       ) then
         tmpcom4x = tmpcom4x + ax(beginres(i)+k-1)
         tmpcom4y = tmpcom4y + ay(beginres(i)+k-1)
         tmpcom4z = tmpcom4z + az(beginres(i)+k-1)
         countatcom = countatcom+1
         endif
         countattmp = countattmp + 1
        enddo

c calculate the COM for SC

         if(countatcom .eq. 4) then
         comSCx(i,1) = tmpcom1x/1.0
         comSCy(i,1) = tmpcom1y/1.0
         comSCz(i,1) = tmpcom1z/1.0
         comSCx(i,2) = tmpcom2x/1.0
         comSCy(i,2) = tmpcom2y/1.0
         comSCz(i,2) = tmpcom2z/1.0
         comSCx(i,3) = tmpcom3x/1.0
         comSCy(i,3) = tmpcom3y/1.0
         comSCz(i,3) = tmpcom3z/1.0
         comSCx(i,4) = tmpcom4x/1.0
         comSCy(i,4) = tmpcom4y/1.0
         comSCz(i,4) = tmpcom4z/1.0
         else
         write(6,*) "something wrong in the Trp ",i
         stop
         endif
       endif

      enddo

C Here we write the pdb file according to the FFCG-2.1

c First the number of residues and the number of beads are written in a HEADER!

      write(6,109) pdbnresAA,pdbnCGbeads
109   format('HEADER number of residues: ',I4," number of beads: ",I4)
      countatCG = 1
      countatAA = 1
      do i = 1,pdbnresAA
       ttt = countatCG
       tta = countatAA
       ttu = beginres(i)
c      write(6,662) countatCG,pdbatname(tta),pdbresname(ttu),
       write(6,662) countatCG," CA ",pdbresname(ttu),
     .pdbresnum(ttu),comBBx(i),comBBy(i),comBBz(i)
       countatSC = 1
       countatCG = countatCG + 1
        do j = 1, (natomresCG(tprestype(i))-1)
        ttb = countatSC
        write(6,662) countatCG,"SC ",pdbresname(ttu),
     .pdbresnum(ttu),comSCx(i,ttb),comSCy(i,ttb),comSCz(i,ttb)
        countatSC = countatSC + 1
        countatCG = countatCG + 1
        enddo
c      countat = countat + natomresCG(tprestype(i))
       countatAA = countatAA + natomresAA(tprestype(i))
      enddo

662   format('ATOM',I7,1X,A4,1X,A3,2X,I4,4X,3F8.3)

      end

max_stars_repo_name: ErikSchutte/Structural-Bioinformatics

max_issues_repo_name: ErikSchutte/Structural-Bioinformatics

max_forks_repo_name: ErikSchutte/Structural-Bioinformatics


text:
      SUBROUTINE XSAVE
C     THE PURPOSE OF THIS ROUTINE IS TO PERFORM THE FUNCTIONS ASSIGNED
C     TO THE SAVE DMAP INSTRUCTION.
C
      COMMON/XVPS/ IVPS(1)
      COMMON/BLANK/ IPAR(1)
      COMMON /OSCENT/ IOSCR(7)
C     GET NUMBER OF PARAMETERS FROM OSCAR
      N = IOSCR(7)*2 + 6
      DO 20 I1 = 8,N,2
C     GET VPS POINTER AND POINTER TO VALUE IN BLANK COMMON.
      J = IOSCR(I1)
      K = IOSCR(I1+1)
C     GET LENGTH OF VALUE FROM VPS
      L = IVPS(J-1)
C     TRANSFER VALUE FROM BLANK COMMON TO VPS
      DO 10 I2 = 1,L
      IVPS(J) = IPAR(K)
      J = J + 1
   10 K = K + 1
   20 CONTINUE
      RETURN
      END


text:
program ch0401
    ! Simple Text I/O
    implicit none

    character(20) :: first_name

    print *, ' type in your first name.'
    print *, ' up to 20 characters'
    read *, first_name
    print *, first_name
end program

max_stars_repo_name: hechengda/fortran

max_issues_repo_name: hechengda/fortran

max_forks_repo_name: hechengda/fortran


text:
module sw_utils
    implicit none
    private
    logical, save :: sw_allocated = .false.
    ! - double precision
    integer, parameter, public :: idp = kind(1.0d0)

    real(idp), parameter, public :: PI = 4.0_idp*atan(1.0)

    ! -------------------------------------
    ! number of point in x and y directions
    ! -------------------------------------
    integer, parameter, public :: NX = 101
    integer, parameter, public :: NY = 101

    ! --------------------
    ! number of time steps
    ! --------------------
    integer, parameter, public :: NT = 10001

    ! ----------------------------------------------
    ! upper and lower bounds for x and y dimensions
    ! ----------------------------------------------
    real(idp), parameter, public :: XUPPER = 5.0_idp     ! xend
    real(idp), parameter, public :: XLOWER = -5.0_idp    ! xstart
    real(idp), parameter, public :: YUPPER = 6.0_idp     ! ystart
    real(idp), parameter, public :: YLOWER = -6.0_idp    ! yend

    ! -----------------------------
    ! start and end point for time
    ! -----------------------------
    real(idp), parameter, public :: TSTART = 0.0_idp
    real(idp), parameter, public :: TEND = 1000.0_idp

    ! ----------------------------
    ! length of x-axis and y-axis
    ! ----------------------------
    real(idp), parameter, public :: XLENGTH = (XUPPER - XLOWER)
    real(idp), parameter, public :: YLENGTH = (YUPPER - YLOWER)

    ! -----------------------
    ! total time or duration
    ! -----------------------
    real(idp), parameter, public :: TLENGTH = (TEND - TSTART)

    ! ---------------------------
    ! spatial grid step dx and dy
    ! ---------------------------
    real(idp), parameter, public :: DX = XLENGTH/dble(NX-1)
    real(idp), parameter, public :: DY = YLENGTH/dble(NY-1)

    ! ----------
    ! time step
    ! ----------
    real(idp), parameter, public :: DT = TLENGTH/dble(NT-1)

    ! ---------------------
    ! gravitation constant
    ! ---------------------
    real(idp), parameter, public :: G = 9.8065_idp

    ! ------------
    ! base height
    ! ------------
    real(idp), parameter, public :: H0 = 1.0_idp

    ! --------------------------
    ! rotation parameter, omega
    ! --------------------------
    real(idp), parameter, public :: OMEGA = 0.1_idp

    ! -------------
    ! viscous drag
    ! -------------
    real(idp), parameter, public :: vd = 0.001_idp

    ! -----------------------
    ! viscosity of the fluid
    ! -----------------------
    real(idp), parameter, public :: nu = 0.000015_idp

    ! -------------
    ! perturbation
    ! -------------
    real(idp), parameter, public :: PERTURB = 0.01_idp

    ! ----------------------
    ! output data file name
    ! ----------------------
    character(len=*), parameter, public :: swfname = "sw_data.nc"

    ! =====================================
    ! public variables' name and constants
    ! =====================================
    real(idp), dimension(1:NX), public :: xdata
    real(idp), dimension(1:NY), public :: ydata
    real(idp), allocatable, public :: tdata(:)
    real(idp), allocatable, public :: xxdata(:, :)
    real(idp), allocatable, public :: yydata(:, :)
    real(idp), allocatable, public :: fdata(:, :)
    real(idp), allocatable, public :: udata(:, :, :)
    real(idp), allocatable, public :: vdata(:, :, :)
    real(idp), allocatable, public :: hdata(:, :, :)

    ! ---------------------------------
    ! public subroutines and functions
    ! ---------------------------------
    public swMeshgrid
    public swAllocateData
    public swDeallocateData
    public swSaveData

contains
    subroutine swAllocateData()
        implicit none
        if ( .not. sw_allocated ) then
            allocate(tdata(NT))
            allocate(xxdata(NX, NY))
            allocate(yydata(NX, NY))
            allocate(fdata(NX, NY))
            allocate(udata(NX, NY, NT))
            allocate(vdata(NX, NY, NT))
            allocate(hdata(NX, NY, NT))
            sw_allocated = .true.
        endif
    end subroutine swAllocateData

    subroutine swDeallocateData()
        implicit none
        if (sw_allocated) then
            deallocate(tdata)
            deallocate(xxdata)
            deallocate(yydata)
            deallocate(fdata)
            deallocate(udata)
            deallocate(vdata)
            deallocate(hdata)
            sw_allocated = .false.
        endif
    end subroutine swDeallocateData

    ! ---
    subroutine swMeshgrid(x, y, xx, yy, nx, ny)
        implicit none
        integer, intent(in) :: nx, ny
        real(kind=8), dimension(1:nx), intent(in) :: x
        real(kind=8), dimension(1:ny), intent(in) :: y
        real(kind=8), dimension(1:nx, 1:ny), intent(out) :: xx, yy
        xx = spread(x, 1, size(y))
        yy = spread(y, 2, size(x))
    end subroutine swMeshgrid

    subroutine swSaveData()
        use netcdf
        implicit none

        ! -----------------------------------------
        ! number of dimensions for large variables
        ! -----------------------------------------
        integer, parameter :: NDIMS3 = 3
        integer, parameter :: NDIMS2 = 2

        ! ------------------------
        ! id for netcdf variables
        ! ------------------------
        integer :: ncid

        ! ----------------
        ! dimension names
        ! ----------------
        character(len=*), parameter :: xname = "X Distance"
        character(len=*), parameter :: yname = "Y Distance"
        character(len=*), parameter :: tname = "T Time"

        ! --------------
        ! dimension ids
        ! --------------
        integer :: xdimID, ydimID, tdimID

        ! --------------------
        ! shape of dimensions
        ! --------------------
        integer :: dimID3(NDIMS3)
        integer :: dimID2(NDIMS2)

        ! ---------------
        ! variable names
        ! ---------------
        character(len=*), parameter :: uname = "U_Velocity"
        character(len=*), parameter :: vname = "V_Velocity"
        character(len=*), parameter :: hname = "H_Level"
        character(len=*), parameter :: fname = "F_Coriolis"

        ! ----------------------------------
        ! Units for the dimension variables
        ! ----------------------------------
        character(len=*), parameter :: units = "units"
        character(len=*), parameter :: u_units = "m/s"
        character(len=*), parameter :: v_units = "m/s"
        character(len=*), parameter :: h_units = "m"
        character(len=*), parameter :: f_units = "1/s"
        character(len=*), parameter :: x_units = "meters_east"
        character(len=*), parameter :: y_units = "meters_north"
        character(len=*), parameter :: t_units = "seconds_forward"

        ! ---------------------------
        ! variable IDs for variables
        ! ---------------------------
        integer :: uvarID, vvarID, hvarID, fvarID

        ! -----------------------------------------
        ! variable IDs for the dimension variables
        ! -----------------------------------------
        integer :: xvarID, yvarID, tvarID

        write(*, *) "Defined all data variables and their IDs."

        ! ----------------
        ! create the file
        ! ----------------
        call check(nf90_create(                                     &
            swfname, or(nf90_clobber, nf90_64bit_offset), ncid)   &
        )

        ! ----------------------
        ! define the dimensions
        ! ----------------------
        call check(nf90_def_dim(ncid, xname, NX, xdimID))
        call check(nf90_def_dim(ncid, yname, NY, ydimID))
        call check(nf90_def_dim(ncid, tname, NT, tdimID))

        ! ---------------------------
        ! create dimension variables
        ! ---------------------------
        call check(nf90_def_var(ncid, xname, NF90_DOUBLE, xdimID, xvarID))
        call check(nf90_def_var(ncid, yname, NF90_DOUBLE, ydimID, yvarID))
        call check(nf90_def_var(ncid, tname, NF90_DOUBLE, tdimID, tvarID))

        ! ----------------------------------
        ! set attributes for the dimensions
        ! ----------------------------------
        call check(nf90_put_att(ncid, xvarID, units, x_units))
        call check(nf90_put_att(ncid, yvarID, units, y_units))
        call check(nf90_put_att(ncid, tvarID, units, t_units))

        write(*, *) "Dimensions' unit attribute added"

        ! ----------------------
        ! create variable slots
        ! ----------------------
        dimID3 = (/ xdimID, ydimID, tdimID /)
        dimID2 = (/ xdimID, ydimID /)

        call check(nf90_def_var(ncid, uname, NF90_DOUBLE, dimID3, uvarID))
        call check(nf90_def_var(ncid, vname, NF90_DOUBLE, dimID3, vvarID))
        call check(nf90_def_var(ncid, hname, NF90_DOUBLE, dimID3, hvarID))
        call check(nf90_def_var(ncid, fname, NF90_DOUBLE, dimID2, fvarID))

        write(*, *) "Data variables created"

        ! ----------------------------------
        ! set attribute for data variabales
        ! ----------------------------------
        call check(nf90_put_att(ncid, uvarID, units, u_units))
        call check(nf90_put_att(ncid, vvarID, units, v_units))
        call check(nf90_put_att(ncid, hvarID, units, h_units))
        call check(nf90_put_att(ncid, fvarID, units, f_units))

        write(*, *) "Data variable unit attribute added"

        ! --------------------------
        ! clode variable definition
        ! --------------------------
        call check(nf90_enddef(ncid))

        ! ------------------------------
        ! write the dimension variables
        ! ------------------------------
        call check(nf90_put_var(ncid, xvarID, xdata))
        write(*, *) "Wrote x dimension variable"
        call check(nf90_put_var(ncid, yvarID, ydata))
        write(*, *) "Wrote y dimension variable"
        call check(nf90_put_var(ncid, tvarID, tdata))
        write(*, *) "Wrote t dimension variable"

        ! -------------------------
        ! write the data variables
        ! -------------------------
        call check(nf90_put_var(ncid, uvarID, udata))
        write(*, *) "Wrote u data variable"
        call check(nf90_put_var(ncid, vvarID, vdata))
        write(*, *) "Wrote v data variable"
        call check(nf90_put_var(ncid, hvarID, hdata))
        write(*, *) "Wrote h (i.e height) data variable"
        call check(nf90_put_var(ncid, fvarID, fdata))
        write(*, *) "Wrote f (i.e coriolis) data variable"

        ! --------------
        ! close the file
        ! ---------------
        call check(nf90_close(ncid))

    contains
        subroutine check(status)
            implicit none
            integer, intent(in) :: status
            if(status /= nf90_noerr) then
                write(*, *) trim(nf90_strerror(status))
                stop "Stopped"
            endif
        end subroutine check
    end subroutine swSaveData
end module sw_utils

max_stars_repo_name: apetcho/scicomp

max_issues_repo_name: apetcho/scicomp

max_forks_repo_name: apetcho/scicomp