text:
	SUBROUTINE DA_RFHR ( iflno, fhdnam, mxword,
     +			     rheadr, nword, iret )
C************************************************************************
C* DA_RFHR								*
C*									*
C* This subroutine reads real header info from a non-GEMPAK		*
C* data source.								*
C*									*
C* DA_RFHR  ( IFLNO, FHDNAM, MXWORD, RHEADR, NWORS, IRET )		*
C*									*
C* Input parameters:							*
C*	IFLNO		INTEGER		File number			*
C*	FHDNAM		CHAR*4		Header name			*
C*	MXWORD		INTEGER		Maximum words to return		*
C*									*
C* Output parameters:							*
C*	RHEADR (NWORD)	REAL		File header			*
C*	NWORD		INTEGER		Header length 			*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* S. Jacobs/NCEP	 8/13	Created					*
C* S. Gilbert/NCEP	 6/14	Fixed calling sequence for da_getfilhdr     *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'../dm/dmcmn.cmn'
C*
	CHARACTER*(*)	fhdnam
	REAL		rheadr (*)
C*
	CHARACTER	fhdn*5
C------------------------------------------------------------------------
	iret = 0
C
C*	Request the data
C
	CALL ST_NULL ( fhdnam, fhdn, lenf, ier )
	CALL DA_GETFILHDR ( iflno, fhdn, mxword, rheadr, nword, iret )
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


! chk_lbound_long_error --
!     Check: does the compiler complain if lbound returns a value that is too large?
!
!     Note 1: I came across that when first using long/large array indices
!     Note 2: This doesn't test the dynamic behaviour of lbound()
!
program chk_lbound_long_error
    implicit none

    integer, parameter                         :: long = selected_int_kind(range(1)+1)
    integer(kind=long), parameter              :: bignumber = huge(1_long)

    integer, dimension(bignumber-10:bignumber) :: value
    integer(kind=long)                         :: i
    integer(kind=long), dimension(1)           :: lower, upper

    write(*, '(a)' ) 'Using "long" integers as array indices:'

    write(*, '(a,i0)') 'Lower bound:   ', lbound(value)
    write(*, '(a,i0)') 'Correct value: ', bignumber-10

    write(*, '(/a)') 'Apparently the compiler does not do a check on the range'

end program chk_lbound_long_error


text:
Despite all her names Pennys owners arent really big fans of Dio, Warren Zevon or Prince.

Penny has many hobbies and interests, some of which include the following:

  Chasing / wrestling with cats
  Sleeping
  Making sure nobody is sleeping
  Users/TaylorStreet Coprophagia
  Humping inanimate objects
  Chewing (she loves Beechnut)
  Vampire hunting
  Playfully biting people in inappropriate places
  Creating situations with electronics for Brad to troubleshoot / fix (to make sure he doesnt get rusty).
  Hiding gifts for people to find later
  Licking
  Stealing panties
  Theoretical physics, particularly string theory

Some of these she learned from her owners (no, not the coprophagia) and some she learned from her buddy Users/AlinaAsrian Alina

20080219 11:16:10 nbsp Great, show everybody my underwear! Users/BradBenedict

20080904 22:09:29 nbsp Yknow, any dog that has a Zevon inspired nickname gets my support. Users/JabberWokky

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!==============================================================================!
  subroutine CalcVISt_SPA_ALL(grid, n)
!------------------------------------------------------------------------------!
!   Computes the turbulent viscosity for RANS models.                          !
!------------------------------------------------------------------------------!
!----------------------------------[Modules]-----------------------------------!
  use all_mod
  use pro_mod
  use les_mod
  use rans_mod
  use Grid_Mod
  use Parameters_Mod
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Grid_Type) :: grid
!-----------------------------------[Locals]-----------------------------------!
  integer :: c, n
  real    :: Xrat, Fv1, lf, Cs
!==============================================================================!

  if(SIMULA == DES_SPA) then
    do c = 1, grid % n_cells
      Xrat     = VIS % n(c)/VISc
      Fv1      = Xrat**3/(Xrat**3 + Cvis1**3)
      VISt(c)  = DENc(material(c)) * Fv1 * VIS % n(c)
    end do
  end if

  if(SIMULA == SPA_ALL) then
    do c = 1, grid % n_cells
      Xrat     = VIS % n(c)/VISc
      Fv1      = Xrat**3/(Xrat**3 + Cvis1**3)
      VISt(c)  = DENc(material(c)) * Fv1 * VIS % n(c)
    end do
  end if

  call Exchange(grid, VISt)

  end subroutine

max_stars_repo_name: mhadziabdic/T-FlowS

max_issues_repo_name: mhadziabdic/T-FlowS

max_forks_repo_name: mhadziabdic/T-FlowS


text:
-makelib xcelium_lib/xil_defaultlib -sv \
  "D:/program/Vivado/Vivado/2018.3/data/ip/xpm/xpm_cdc/hdl/xpm_cdc.sv" \
  "D:/program/Vivado/Vivado/2018.3/data/ip/xpm/xpm_memory/hdl/xpm_memory.sv" \
-endlib
-makelib xcelium_lib/xpm \
  "D:/program/Vivado/Vivado/2018.3/data/ip/xpm/xpm_VCOMP.vhd" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  "../../../bd/cm3_core/ipshared/ea75/rtl/cm3_dap_ahb_ap.v" \
  "../../../bd/cm3_core/ipshared/ea75/rtl/cm3_wic.v" \
  "../../../bd/cm3_core/ipshared/ea75/rtl/cm3_nvic.v" \
  "../../../bd/cm3_core/ipshared/ea75/rtl/cm3_dwt.v" \
  "../../../bd/cm3_core/ipshared/ea75/rtl/AhbToAxi.v" \
  "../../../bd/cm3_core/ipshared/ea75/rtl/cm3_mpu.v" \
  "../../../bd/cm3_core/ipshared/ea75/rtl/cm3_bus_matrix.v" \
  "../../../bd/cm3_core/ipshared/ea75/rtl/models.v" \
  "../../../bd/cm3_core/ipshared/ea75/rtl/cmsdk.v" \
  "../../../bd/cm3_core/ipshared/ea75/rtl/cm3_tpiu.v" \
  "../../../bd/cm3_core/ipshared/ea75/rtl/dapswjdp.v" \
  "../../../bd/cm3_core/ipshared/ea75/rtl/cm3_fpb.v" \
  "../../../bd/cm3_core/ipshared/ea75/rtl/cm3_dpu.v" \
  "../../../bd/cm3_core/ipshared/ea75/rtl/cm3_itm.v" \
  "../../../bd/cm3_core/ipshared/ea75/rtl/cm3_etm.v" \
  "../../../bd/cm3_core/ipshared/ea75/rtl/cortexm3.v" \
  "../../../bd/cm3_core/ipshared/ea75/rtl/cortexm3_integration.v" \
  "../../../bd/cm3_core/ip/cm3_core_CORTEXM3_AXI_0_0/sim/cm3_core_CORTEXM3_AXI_0_0.v" \
  "../../../bd/cm3_core/ip/cm3_core_clk_wiz_0_0/cm3_core_clk_wiz_0_0_clk_wiz.v" \
  "../../../bd/cm3_core/ip/cm3_core_clk_wiz_0_0/cm3_core_clk_wiz_0_0.v" \
-endlib
-makelib xcelium_lib/lib_cdc_v1_0_2 \
  "../../../../../bd/cm3_core/ipshared/ef1e/hdl/lib_cdc_v1_0_rfs.vhd" \
-endlib
-makelib xcelium_lib/proc_sys_reset_v5_0_13 \
  "../../../../../bd/cm3_core/ipshared/8842/hdl/proc_sys_reset_v5_0_vh_rfs.vhd" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  "../../../bd/cm3_core/ip/cm3_core_proc_sys_reset_0_0/sim/cm3_core_proc_sys_reset_0_0.vhd" \
-endlib
-makelib xcelium_lib/generic_baseblocks_v2_1_0 \
  "../../../../../bd/cm3_core/ipshared/b752/hdl/generic_baseblocks_v2_1_vl_rfs.v" \
-endlib
-makelib xcelium_lib/axi_infrastructure_v1_1_0 \
  "../../../../../bd/cm3_core/ipshared/ec67/hdl/axi_infrastructure_v1_1_vl_rfs.v" \
-endlib
-makelib xcelium_lib/axi_register_slice_v2_1_18 \
  "../../../../../bd/cm3_core/ipshared/cc23/hdl/axi_register_slice_v2_1_vl_rfs.v" \
-endlib
-makelib xcelium_lib/fifo_generator_v13_2_3 \
  "../../../../../bd/cm3_core/ipshared/64f4/simulation/fifo_generator_vlog_beh.v" \
-endlib
-makelib xcelium_lib/fifo_generator_v13_2_3 \
  "../../../../../bd/cm3_core/ipshared/64f4/hdl/fifo_generator_v13_2_rfs.vhd" \
-endlib
-makelib xcelium_lib/fifo_generator_v13_2_3 \
  "../../../../../bd/cm3_core/ipshared/64f4/hdl/fifo_generator_v13_2_rfs.v" \
-endlib
-makelib xcelium_lib/axi_data_fifo_v2_1_17 \
  "../../../../../bd/cm3_core/ipshared/c4fd/hdl/axi_data_fifo_v2_1_vl_rfs.v" \
-endlib
-makelib xcelium_lib/axi_crossbar_v2_1_19 \
  "../../../../../bd/cm3_core/ipshared/6c9d/hdl/axi_crossbar_v2_1_vl_rfs.v" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  "../../../bd/cm3_core/ip/cm3_core_xbar_0/sim/cm3_core_xbar_0.v" \
-endlib
-makelib xcelium_lib/util_vector_logic_v2_0_1 \
  "../../../../../bd/cm3_core/ipshared/2137/hdl/util_vector_logic_v2_0_vl_rfs.v" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  "../../../bd/cm3_core/ip/cm3_core_util_vector_logic_0_0/sim/cm3_core_util_vector_logic_0_0.v" \
-endlib
-makelib xcelium_lib/xlconstant_v1_1_5 \
  "../../../../../bd/cm3_core/ipshared/4649/hdl/xlconstant_v1_1_vl_rfs.v" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  "../../../bd/cm3_core/ip/cm3_core_xlconstant_0_0/sim/cm3_core_xlconstant_0_0.v" \
  "../../../bd/cm3_core/ip/cm3_core_xlconstant_0_1/sim/cm3_core_xlconstant_0_1.v" \
  "../../../bd/cm3_core/ip/cm3_core_xlconstant_1_0/sim/cm3_core_xlconstant_1_0.v" \
-endlib
-makelib xcelium_lib/axi_lite_ipif_v3_0_4 \
  "../../../../../bd/cm3_core/ipshared/66ea/hdl/axi_lite_ipif_v3_0_vh_rfs.vhd" \
-endlib
-makelib xcelium_lib/interrupt_control_v3_1_4 \
  "../../../../../bd/cm3_core/ipshared/a040/hdl/interrupt_control_v3_1_vh_rfs.vhd" \
-endlib
-makelib xcelium_lib/axi_gpio_v2_0_20 \
  "../../../../../bd/cm3_core/ipshared/a7c9/hdl/axi_gpio_v2_0_vh_rfs.vhd" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  "../../../bd/cm3_core/ip/cm3_core_axi_gpio_0_1/sim/cm3_core_axi_gpio_0_1.vhd" \
-endlib
-makelib xcelium_lib/lib_pkg_v1_0_2 \
  "../../../../../bd/cm3_core/ipshared/0513/hdl/lib_pkg_v1_0_rfs.vhd" \
-endlib
-makelib xcelium_lib/lib_srl_fifo_v1_0_2 \
  "../../../../../bd/cm3_core/ipshared/51ce/hdl/lib_srl_fifo_v1_0_rfs.vhd" \
-endlib
-makelib xcelium_lib/axi_uartlite_v2_0_22 \
  "../../../../../bd/cm3_core/ipshared/7371/hdl/axi_uartlite_v2_0_vh_rfs.vhd" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  "../../../bd/cm3_core/ip/cm3_core_axi_uartlite_0_0/sim/cm3_core_axi_uartlite_0_0.vhd" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  "../../../bd/cm3_core/sim/cm3_core.v" \
  "../../../bd/cm3_core/ip/cm3_core_swdio_tri_buffer_0_0/sim/cm3_core_swdio_tri_buffer_0_0.v" \
-endlib
-makelib xcelium_lib/axi_protocol_converter_v2_1_18 \
  "../../../../../bd/cm3_core/ipshared/7a04/hdl/axi_protocol_converter_v2_1_vl_rfs.v" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  "../../../bd/cm3_core/ip/cm3_core_auto_pc_0/sim/cm3_core_auto_pc_0.v" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  glbl.v
-endlib


text:
      PROGRAM CCPPOLE
C
C Define the error file, the Fortran unit number, the workstation type,
C and the workstation ID to be used in calls to GKS routines.
C
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)   ! NCGM
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=8,  IWKID=1)   ! X Windows
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=11, IWKID=1)   ! PDF
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=20, IWKID=1)   ! PostScript
C
      PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)
C
C Parameterize the number of latitudes, the number of longitudes, the
C sizes of the real and integer workspaces, the size of the area map,
C and the size of the arrays used by ARSCAM for X and Y coordinates.
C
      PARAMETER (NLON=361,NLAT=181,LRWK=5000,LIWK=5000,LAMA=600000)
      PARAMETER (NCRA=20000)
C
C Declare the data array.
C
      DIMENSION ZDAT(NLON,NLAT)
C
C Declare the area map array.
C
      DIMENSION IAMA(LAMA)
C
C Declare the real and integer workspace arrays for CONPACK.
C
      DIMENSION RWRK(LRWK),IWRK(LIWK)
C
C Declare arrays for ARSCAM and MAPGRM to use in calls to COLRAM and
C COLRLL, respectively.  XCRA and YCRA hold X and Y coordinates; IAIA
C and IGIA hold area identifiers and group identifiers.
C
      DIMENSION XCRA(NCRA),YCRA(NCRA),IAIA(10),IGIA(10)
C
C Declare an array in which to retrieve aspect source flags.
C
      DIMENSION IASF(13)
C
C Declare a routine to color the areas represented by the area map.
C
      EXTERNAL COLRAM
C
C Declare a routine to draw contour lines over land only.
C
      EXTERNAL COLRCL
C
C Declare a routine to draw lat/lon lines over ocean only.
C
      EXTERNAL COLRLL
C
C Define the values to be used for GKS aspect source flags.
C
      DATA IASF / 13*1 /
C
C Open GKS.
C
      CALL GOPKS (IERRF, ISZDM)
      CALL GOPWK (IWKID, LUNIT, IWTYPE)
      CALL GACWK (IWKID)
C
C Turn off the clipping indicator.
C
      CALL GSCLIP (0)
C
C Set all aspect source flags to "individual".
C
      CALL GSASF (IASF)
C
C Force solid fill.
C
      CALL GSFAIS (1)
C
C Define 16 different color indices, for indices 0 through 15.  The
C color corresponding to index 0 is black and the color corresponding
C to index 1 is white.  Colors 2 and 3 are used for alternate contour
C bands over ocean.  Colors 4 through 16 are used for contour bands
C over land.
C
      CALL GSCR (IWKID, 0,0.000,0.000,0.000)
      CALL GSCR (IWKID, 1,1.000,1.000,1.000)
      CALL GSCR (IWKID, 2,0.000,1.000,1.000)
      CALL GSCR (IWKID, 3,0.000,0.850,0.850)
      CALL GSCR (IWKID, 4,0.700,0.700,0.700)
      CALL GSCR (IWKID, 5,0.750,0.500,1.000)
      CALL GSCR (IWKID, 6,0.500,0.000,1.000)
      CALL GSCR (IWKID, 7,0.000,0.000,1.000)
      CALL GSCR (IWKID, 8,0.000,0.500,1.000)
      CALL GSCR (IWKID, 9,0.000,1.000,0.600)
      CALL GSCR (IWKID,10,0.000,1.000,0.000)
      CALL GSCR (IWKID,11,0.700,1.000,0.000)
      CALL GSCR (IWKID,12,1.000,1.000,0.000)
      CALL GSCR (IWKID,13,1.000,0.750,0.000)
      CALL GSCR (IWKID,14,1.000,0.380,0.380)
      CALL GSCR (IWKID,15,1.000,0.000,0.380)
      CALL GSCR (IWKID,16,1.000,0.000,0.000)
C
C Generate an array of test data.  It is important that the data should
C represent a continous function around the globe.  What is used here
C is a simple trigonometric function of latitude and longitude.
C
      ZMIN= 1.E36
      ZMAX=-1.E36
C
      DO 102 I=1,NLON
         RLON=.017453292519943*(-180.+360.*REAL(I-1)/REAL(NLON-1))
         DO 101 J=1,NLAT
            RLAT=.017453292519943*(-90.+180.*REAL(J-1)/REAL(NLAT-1))
            ZDAT(I,J)=.5*COS(8.*RLAT)+.25*COS(RLAT)*SIN(4.*RLON)
            ZMIN=MIN(ZMIN,ZDAT(I,J))
            ZMAX=MAX(ZMAX,ZDAT(I,J))
 101     CONTINUE
 102  CONTINUE
C
C Reduce the test data to the desired range.
C
      DO 104 I=1,NLON
         DO 103 J=1,NLAT
            ZDAT(I,J)=((ZDAT(I,J)-ZMIN)/(ZMAX-ZMIN))*130.-10.
 103     CONTINUE
 104  CONTINUE
C
C Put a label at the top of the plot.
C
      CALL SET    (0.,1.,0.,1.,0.,1.,0.,1.,1)
      CALL PLCHHQ (.5,.975,
     +     'CONTOUR BANDS ON A STEREOGRAPHIC PROJECTION',
     +     .02,0.,0.)
C
C Initialize the area map.
C
      CALL ARINAM (IAMA,LAMA)
C
C Tell Ezmap where to put the plot on the plotter frame.
C
      CALL MAPPOS (.03,.97,.01,.95)
C
C Tell Ezmap to use a stereographic projection.
C
      CALL MAPROJ ('ST',90.,0.,0.)
C
C Tell Ezmap to use a 90-degree distance to each of the four edges
C of the map.
C
      CALL MAPSET ('AN - ANGLES',90.,90.,90.,90.)
C
C Initialize Ezmap.
C
      CALL MAPINT
C
C Put continental outlines in the area map.
C
      CALL MAPBLA (IAMA)
C
C Tell CONPACK not to do the SET call (since it's already been done)
C and to use mapping function 1 (EZMAP background).
C
      CALL CPSETI ('SET - DO-SET-CALL FLAG',0)
      CALL CPSETI ('MAP - MAPPING FLAG',1)
C
C Tell CONPACK what ranges of X and Y coordinates to send into the
C mapping function.
C
      CALL CPSETR ('XC1 - X COORDINATE AT I=1',-180.)
      CALL CPSETR ('XCM - X COORDINATE AT I=M',+180.)
      CALL CPSETR ('YC1 - Y COORDINATE AT J=1', -90.)
      CALL CPSETR ('YCN - Y COORDINATE AT J=N', +90.)
C
C Tell CONPACK exactly what contour levels to use.
C
      CALL CPSETI ('CLS - CONTOUR LEVEL SELECTOR',1)
      CALL CPSETR ('CMN - CONTOUR LEVEL MINIMUM',0.)
      CALL CPSETR ('CMX - CONTOUR LEVEL MAXIMUM',110.)
      CALL CPSETR ('CIS - CONTOUR INTERVAL SPECIFIER',10.)
C
C Tell CONPACK what to use as the out-of-range flag.  This is the
C value returned by the Ezmap routine MAPTRA for off-map points.
C
      CALL CPSETR ('ORV - OUT-OF-RANGE VALUE',1.E12)
C
C Initialize the drawing of the contour plot.
C
      CALL CPRECT (ZDAT,NLON,NLON,NLAT,RWRK,LRWK,IWRK,LIWK)
C
C Add contour lines to the area map.
C
      CALL CPCLAM (ZDAT,RWRK,IWRK,IAMA)
C
C Color the map.
C
      CALL ARSCAM (IAMA,XCRA,YCRA,NCRA,IAIA,IGIA,10,COLRAM)
C
C Switch the current polyline color to black.
C
      CALL GSPLCI (0)
C
C Outline the continents and draw lines of latitude and longitude over
C the ocean only.
C
      CALL MAPGRM (IAMA,XCRA,YCRA,NCRA,IAIA,IGIA,10,COLRLL)
      CALL MAPSTC ('OU - OUTLINE DATASET','CO')
      CALL MAPLOT
C
C Draw the contour lines over land only.
C
      CALL CPCLDM (ZDAT,RWRK,IWRK,IAMA,COLRCL)
C
C Advance the frame.
C
      CALL FRAME
C
C Close GKS.
C
      CALL GDAWK (IWKID)
      CALL GCLWK (IWKID)
      CALL GCLKS
C
C Done.
C
      STOP
C
      END

      SUBROUTINE COLRAM (XCRA,YCRA,NCRA,IAIA,IGIA,NAGI)
C
      DIMENSION XCRA(*),YCRA(*),IAIA(*),IGIA(*)
C
C This routine is called to color an area from an area map.  Its
C coordinates are given by the NCRA coordinates in the arrays XCRA and
C YCRA.  For each I from 1 to NAGI, IAIA(I) is the area identifier of
C the area relative to the group whose group identifier is IGIA(I).
C
C Find the area identifier for the area relative to groups 1 and 3.
C The first of these tells us whether the area is over land or water,
C and the second tells us what contour band the area is in.
C
      IAI1=-1
      IAI3=-1
C
      DO 101 I=1,NAGI
         IF (IGIA(I).EQ.1) IAI1=IAIA(I)
         IF (IGIA(I).EQ.3) IAI3=IAIA(I)
 101  CONTINUE
C
C Color-fill the area, using two slightly-different shades of blue
C over water (so that the contour bands will be minimally visible)
C and brighter colors over land.
C
      IF (IAI1.GT.0) THEN
         IF (MAPACI(IAI1).EQ.1) THEN
            CALL GSFACI (2+MOD(IAI3,2))
            CALL GFA (NCRA-1,XCRA,YCRA)
         ELSE
            IF (IAI3.GE.1.AND.IAI3.LE.13) THEN
               CALL GSFACI (IAI3+3)
               CALL GFA (NCRA-1,XCRA,YCRA)
            END IF
         END IF
      END IF
C
C Done.
C
      RETURN
C
      END

      SUBROUTINE COLRCL (XCRA,YCRA,NCRA,IAIA,IGIA,NAGI)
C
      DIMENSION XCRA(*),YCRA(*),IAIA(*),IGIA(*)
C
C This routine is called to draw a portion of a contour line which is
C wholly contained in some area of an area map.  Its coordinates are
C given by the NCRA coordinates in the arrays XCRA and YCRA.  For each
C I from 1 to NAGI, IAIA(I) is the area identifier of the area relative
C to the group whose group identifier is IGIA(I).
C
C Find the area identifier for the area relative to groups 1 and 3.
C The first of these tells us whether the area is over land or water,
C and the second tells us what contour band the area is in.
C
      IAI1=-1
      IAI3=-1
C
      DO 101 I=1,NAGI
         IF (IGIA(I).EQ.1) IAI1=IAIA(I)
         IF (IGIA(I).EQ.3) IAI3=IAIA(I)
 101  CONTINUE
C
C Draw the line only if the area it is in is over land.
C
      IF (IAI1.GT.0.AND.MAPACI(IAI1).NE.1) CALL GPL (NCRA,XCRA,YCRA)
C
C Done.
C
      RETURN
C
      END

      SUBROUTINE COLRLL (XCRA,YCRA,NCRA,IAIA,IGIA,NAGI)
C
      DIMENSION XCRA(*),YCRA(*),IAIA(*),IGIA(*)
C
C This routine is called to draw a portion of a line of latitude or
C longitude which is wholly contained in some area of an area map.  Its
C coordinates are given by the NCRA coordinates in the arrays XCRA and
C YCRA.  For each I from 1 to NAGI, IAIA(I) is the area identifier of
C the area relative to the group whose group identifier is IGIA(I).
C
C Find the area identifier for the area relative to group 1, which will
C tell us whether the area is over land or water.
C
      IAI1=-1
C
      DO 101 I=1,NAGI
         IF (IGIA(I).EQ.1.AND.IAIA(I).GT.0) IAI1=IAIA(I)
 101  CONTINUE
C
C Draw the line only if it is over water.
C
      IF (IAI1.GT.0.AND.MAPACI(IAI1).EQ.1) CALL GPL (NCRA,XCRA,YCRA)
C
C Done.
C
      RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
!>
module parameters

  implicit none

  ! numeric parameters
  ! - <TOL>: double precision tolerance value.
  double precision , parameter :: TOL = 1.0D-10

end module parameters

max_stars_repo_name: dgsaf/scattering

max_issues_repo_name: dgsaf/scattering

max_forks_repo_name: dgsaf/scattering


text:
!##############################################################################
!# ****************************************************************************
!# <name> basicgeometry </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module contains the basic entity of the whole FEAT world: a point,
!# a coordinate system and a set of routines to maintain that.
!#
!# The following routines can be found in this module:
!#
!#  1.) bgeom_initCoordSys2D, bgeom_initCoordSys3D
!#      -> Creates a 2D/3D coordinate system from an optionally given origin,
!#         rotation and scaling factor.
!#
!#  2.) bgeom_transformPoint2D
!#      -> Transforms a 2D point with coordinates relative to a given
!#         coordinate system to world coordinates.
!#
!#  3.) bgeom_transformBackPoint2D
!#      -> Transforms a 2D point with world coordinates to coordinates relative
!#         to a given coordinate system.
!#
!#  4.) bgeom_transformPoint3D
!#      -> Transforms a 3D point with coordinates relative to a given
!#         coordinate system to world coordinates.
!#
!#  5.) bgeom_transformBackPoint3D
!#      -> Transforms a 3D point with world coordinates to coordinates relative
!#         to a given coordinate system.

!# </purpose>
!##############################################################################

module basicgeometry

!$ use omp_lib
  use fsystem

  implicit none

  private

!<constants>

!<constantblock description="Dimension constants">

  ! Dimension constant for 1D triangulations.
  integer, parameter, public :: NDIM1D = 1

  ! Dimension constant for 2D triangulations.
  integer, parameter, public :: NDIM2D = 2

  ! Dimension constant for 3D triangulations.
  integer, parameter, public :: NDIM3D = 3

!</constantblock>

!<constantblock description="Element shape identifiers">
  ! Unknown shape
  integer(I32), parameter, public :: BGEOM_SHAPE_UNKNOWN = 0_I32

  ! Line shape (1D)
  integer(I32), parameter, public :: BGEOM_SHAPE_LINE  = 1000002_I32

  ! Triangle shape (2D)
  integer(I32), parameter, public :: BGEOM_SHAPE_TRIA  = 2000303_I32

  ! Quadrilateral shape (2D)
  integer(I32), parameter, public :: BGEOM_SHAPE_QUAD  = 2000404_I32

  ! Tetrahedron shape (3D)
  integer(I32), parameter, public :: BGEOM_SHAPE_TETRA = 3040604_I32

  ! Hexahedron shape (3D)
  integer(I32), parameter, public :: BGEOM_SHAPE_HEXA  = 3061208_I32

  ! Pyramid shape (3D)
  integer(I32), parameter, public :: BGEOM_SHAPE_PYRA  = 3050805_I32

  ! Prism shape (3D)
  integer(I32), parameter, public :: BGEOM_SHAPE_PRISM = 3050906_I32

!</constantblock>

!</constants>

  !<types>

  !<typeblock>

  ! The point structure for 2D points.
  type t_point2D

    ! X-coordinate
    real(DP) :: X

    ! Y-coordinate
    real(DP) :: Y

  end type

  public :: t_point2D

  !<typeblock>

  !</typeblock>

  ! The point structure for 3D points.
  type t_point3D

    ! X-coordinate
    real(DP) :: X

    ! Y-coordinate
    real(DP) :: Y

    ! Z-coordinate
    real(DP) :: Z

  end type

  public :: t_point3D

  !</typeblock>

  !<typeblock>

  ! The point structure for regular 2D coordinate systems, consisting of
  ! a rotated X- and Y-axes
  type t_coordinateSystem2D

    ! Coordinates of the origin
    real(DP), dimension(2) :: Dorigin = (/0.0_DP,0.0_DP/)

    ! Rotation angle; 0..2*PI
    real(DP) :: drotation = 0.0_DP

    ! precalculated value: sin(rotation); for quicker calculations
    real(DP) :: dsin_rotation = 0.0_DP

    ! precalculated value: cos(rotation); for quicker calculations
    real(DP) :: dcos_rotation = 1.0_DP

    ! scaling factor of the coordinate system; usually = 1.0
    real(DP) :: dscalingFactor = 1.0_DP

  end type

  public :: t_coordinateSystem2D

  !</typeblock>

  !<typeblock>

  ! The point structure for regular 3D coordinate systems, consisting of
  ! a rotated X-, Y-axes and Z-axes
  type t_coordinateSystem3D

    ! Coordinates of the origin
    real(DP), dimension(3) :: Dorigin = (/0.0_dp,0.0_dp,0.0_dp/)

    ! Rotation angle; 0..2*PI
    real(DP) :: drotationX = 0.0_dp

    real(DP) :: drotationY = 0.0_dp

    real(DP) :: drotationZ = 0.0_dp

    ! scaling factor of the coordinate system; usually = 1.0
    real(DP) :: dscalingFactor = 1.0_dp

  end type

    public :: t_coordinateSystem3D

  !</typeblock>

  !</types>

  public :: bgeom_initCoordSys2D
  public :: bgeom_initCoordSys3D
  public :: bgeom_transformPoint2D
  public :: bgeom_transformBackPoint2D
  public :: bgeom_transformPoint3D
  public :: bgeom_transformBackPoint3D

contains

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

!<subroutine>

  pure subroutine bgeom_initCoordSys2D(rcoordSys, Dorigin, drotation, &
                                       dscalingFactor)

!<description>
  ! Creates a 2D coordinate system from an optionally given origin, rotation
  ! and scaling factor.
!</description>

!<input>
  ! OPTIONAL: The origin of the coordinate system.
  ! Set to (/ 0.0_DP, 0.0_DP /) if not given.
  real(DP), dimension(:), optional,    intent(in)  :: Dorigin

  ! OPTIONAL: The rotation angle of the coordinate system.
  ! Angle range: 0..2*PI
  ! Set to 0.0_DP if not given.
  real(DP), optional,                  intent(in)  :: drotation

  ! The scaling factor. Should be != 0.0_DP
  ! Set to 1.0_DP if not given.
  real(DP), optional,                  intent(in)  :: dscalingFactor

!</input>

!<output>
  ! A t_coordinateSystem2D structure to be written.
  type(t_coordinateSystem2D),          intent(out) :: rcoordSys

!</output>

!</subroutine>

    ! Set the origin, if given.
    if (present(Dorigin)) then
      rcoordSys%Dorigin = (/ Dorigin(1), Dorigin(2) /)
    else
      rcoordSys%Dorigin = (/ 0.0_DP, 0.0_DP /)
    end if

    ! Set the rotation, if given.
    if (present(drotation)) then

      rcoordSys%drotation = drotation

      ! Calculate SIN and COS of rotation angle
      rcoordSys%dsin_rotation = sin(drotation)
      rcoordSys%dcos_rotation = cos(drotation)

    else

      rcoordSys%drotation = 0.0_DP

      ! Set SIN and COS values
      rcoordSys%dsin_rotation = 0.0_DP
      rcoordSys%dcos_rotation = 1.0_DP

    end if

    ! Set the scaling factor, if given.
    if (present(dscalingFactor)) then
      rcoordSys%dscalingFactor = dscalingFactor
    else
      rcoordSys%dscalingFactor = 1.0_DP
    endif

    ! That is it

  end subroutine

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

!<subroutine>

  pure subroutine bgeom_transformPoint2D(rcoordSys, DpointIn, DpointOut)

!<description>
  ! Transforms a given 2D point with coordinates relative to the given
  ! coordinate system to world coordinates.
!</description>

!<input>
  ! The 2D coordinate system that is to be applied to the given point.
  type(t_coordinateSystem2D),  intent(in)  :: rcoordSys

  ! The 2D point that is to be transformed.
  real(DP), dimension(:),      intent(in)  :: DpointIn
!</input>

!<output>
  ! The transformed 2D point.
  real(DP), dimension(:),      intent(out) :: DpointOut
!</output>

!</subroutine>

    ! Check if the rotation of the coordinate system is non-zero.
    if (rcoordSys%drotation .ne. 0.0_DP) then

      ! Apply rotation to the input coordinates.
      DpointOut(1) = (rcoordSys%dcos_rotation * DpointIn(1)) - &
                     (rcoordSys%dsin_rotation * DpointIn(2))
      DpointOut(2) = (rcoordSys%dsin_rotation * DpointIn(1)) + &
                     (rcoordSys%dcos_rotation * DpointIn(2))

    else

      ! No rotation in the coordinate system, so simply copy the input coords.
      DpointOut(1) = DpointIn(1)
      DpointOut(2) = DpointIn(2)

    end if

    ! Now scale the given coordinates by the scaling factor and translate them
    ! by the origin of our coordinate system.
    DpointOut(1) = (rcoordSys%dscalingFactor * DpointOut(1)) + &
                    rcoordSys%Dorigin(1)
    DpointOut(2) = (rcoordSys%dscalingFactor * DpointOut(2)) + &
                    rcoordSys%Dorigin(2)

    ! That is it!

  end subroutine

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

!<subroutine>

  pure subroutine bgeom_transformBackPoint2D(rcoordSys, DpointIn, DpointOut)

!<description>
  ! Transform a 2D point given in world coordinates to coordinates relative
  ! to a given 2D coordinate system.
  ! This subroutine is the inverse routine of bgeom_transformPoint2D.
!</description>

!<input>
  ! The 2D coordinate system of our input point.
  type(t_coordinateSystem2D), intent(in)  :: rcoordSys

  ! The 2D point that is to be transformed, relative to the given coordinate
  ! system.
  real(DP), dimension(:),     intent(in)  :: DpointIn
!</input>

!<output>
  ! The transformed 2D point, in world coordinates.
  real(DP), dimension(:),     intent(out) :: DpointOut
!</output>

!</subroutine>

    ! local variables
    real(DP) :: X,Y

    ! Translate the given point by the negatives of our coordinate system
    ! origin.
    X = DpointIn(1) - rcoordSys%Dorigin(1)
    Y = DpointIn(2) - rcoordSys%Dorigin(2)

    ! Now scale the point by the inverse of our coordinate system sclaing
    ! factor, of course only if it is non-zero ...
    if (rcoordSys%dscalingFactor .ne. 0.0_DP) then

      X = X / rcoordSys%dscalingFactor
      Y = Y / rcoordSys%dscalingFactor

    else

      ! ... or set the result to zero and exit subroutine otherwise.
      DpointOut(1) = 0.0_DP
      DpointOut(2) = 0.0_DP
      return

    endif

    ! Check if the rotation of the coordinate system is non-zero.
    if (rcoordSys%drotation .ne. 0.0_DP) then

      ! Apply rotation to the input coordinates.
      DpointOut(1) = ( rcoordSys%dcos_rotation * X) + (rcoordSys%dsin_rotation * Y)
      DpointOut(2) = (-rcoordSys%dsin_rotation * X) + (rcoordSys%dcos_rotation * Y)

    else

      ! No rotation in the coordinate system, so simply copy the input coords.
      DpointOut(1) = X
      DpointOut(2) = Y

    end if

    ! That is it

  end subroutine

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

!<subroutine>

  pure subroutine bgeom_initCoordSys3D(rcoordSys, Dorigin, drotationX, drotationY,&
                                       drotationZ, dscalingFactor)

!<description>
  ! Creates a 3D coordinate system from an optionally given origin, rotation angles
  ! and scaling factor.
!</description>

!<input>
  ! OPTIONAL: The origin of the coordinate system.
  ! Set to (/ 0.0_dp, 0.0_dp, 0.0_dp /) if not given.
  real(DP), dimension(:), optional,    intent(in)  :: Dorigin

  ! OPTIONAL: The rotation angles of the coordinate system.
  ! Angle range: 0..2*PI
  ! Set to 0.0_DP if not given.
  real(DP), optional,                  intent(in)  :: drotationX
  real(DP), optional,                  intent(in)  :: drotationY
  real(DP), optional,                  intent(in)  :: drotationZ

  ! The scaling factor. Should be != 0.0_DP
  ! Set to 1.0_DP if not given.
  real(DP), optional,                  intent(in)  :: dscalingFactor

!</input>

!<output>
  ! A t_coordinateSystem3D structure to be written.
  type(t_coordinateSystem3D),          intent(out) :: rcoordSys

!</output>

!</subroutine>

    ! Set the origin, if given.
    if (present(Dorigin)) then
      rcoordSys%Dorigin = (/ Dorigin(1), Dorigin(2), Dorigin(3) /)
    else
      rcoordSys%Dorigin = (/ 0.0_dp, 0.0_dp, 0.0_dp /)
    end if

    ! Set the rotation, if given.
    if (present(drotationX)) then
      rcoordSys%drotationX = drotationX
    else
      rcoordSys%drotationX = 0.0_DP
    end if

    ! Set the rotation, if given.
    if (present(drotationY)) then
      rcoordSys%drotationY = drotationY
    else
      rcoordSys%drotationY = 0.0_DP
    end if

    ! Set the rotation, if given.
    if (present(drotationZ)) then
      rcoordSys%drotationZ = drotationZ
    else
      rcoordSys%drotationZ = 0.0_DP
    end if

    ! Set the scaling factor, if given.
    if (present(dscalingFactor)) then
      rcoordSys%dscalingFactor = dscalingFactor
    else
      rcoordSys%dscalingFactor = 1.0_DP
    endif

    ! That is it

  end subroutine

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

!<subroutine>

  pure subroutine bgeom_transformPoint3D(rcoordSys, DpointIn, DpointOut)

!<description>
  ! Transforms a 3D point in model coordinates to world coordinates
  ! This routine may be inefficient if a large number of of points is to
  ! be transformed. Then a matrix based transformation would be more efficient.
  ! But this would require a rewriting of the mechanism how routines, that call
  ! this function(i.e. geom_isingeometry) work.
!</description>

!<input>
  ! The 3D coordinate system that is to be applied to the given point.
  type(t_coordinateSystem3D),  intent(in)  :: rcoordSys

  ! The 3D point that is to be transformed.
  real(DP), dimension(:),      intent(in)  :: DpointIn
!</input>

!<output>
  ! The transformed 3D point.
  real(DP), dimension(:),      intent(out) :: DpointOut
!</output>

!</subroutine>

  real(dp) :: sx,sy,sz,cx,cy,cz,x1,y1,z1,x2,y2,z2

  ! Check if the rotation of the coordinate system is non-zero.
  if (rcoordSys%drotationX .ne. 0.0_DP) then

    ! setup
    sx = sin(rcoordSys%drotationX)
    cx = cos(rcoordSys%drotationX)
    sy = sin(rcoordSys%drotationY)
    cy = cos(rcoordSys%drotationY)
    sz = sin(rcoordSys%drotationZ)
    cz = cos(rcoordSys%drotationZ)

    x1 =  DpointIn(1) * cz +  DpointIn(2) * sz
    y1 =  DpointIn(2) * cz -  DpointIn(1) * sz
    z1 =  DpointIn(3)
    x2 = x1 * cy + z1 * sy
    y2 = z1
    z2 = z1 * cy - x1 * sy
    DpointOut(1) = x2
    DpointOut(2) = y2 * cx + z1 * sx
    DpointOut(3) = z2 * cx - x1 * sx

  else

    ! No rotation in the coordinate system, so simply copy the input coords.
    DpointOut(1) = DpointIn(1)
    DpointOut(2) = DpointIn(2)
    DpointOut(3) = DpointIn(3)

  end if

  ! Now scale the given coordinates by the scaling factor and translate them
  ! by the origin of our coordinate system.
  DpointOut(1) = (rcoordSys%dscalingFactor * DpointOut(1)) + &
                  rcoordSys%Dorigin(1)
  DpointOut(2) = (rcoordSys%dscalingFactor * DpointOut(2)) + &
                  rcoordSys%Dorigin(2)
  DpointOut(3) = (rcoordSys%dscalingFactor * DpointOut(3)) + &
                  rcoordSys%Dorigin(3)

  ! That is it!

  end subroutine

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

!<subroutine>

  pure subroutine bgeom_transformBackPoint3D(rcoordSys, DpointIn, DpointOut)

!<description>
  ! Transform a 3d point from world coordinates to model coordinates
!</description>

!<input>
  ! The 3D coordinate system of our input point.
  type(t_coordinateSystem3D), intent(in)  :: rcoordSys

  ! The 3D point that is to be transformed, relative to the given coordinate
  ! system.
  real(DP), dimension(:),     intent(in)  :: DpointIn
!</input>

!<output>
  ! The transformed 3D point, in world coordinates.
  real(DP), dimension(:),     intent(out) :: DpointOut
!</output>

!</subroutine>

  ! local variables
  real(DP) :: X,Y,Z

  real(dp) :: sx,sy,sz,cx,cy,cz,x1,y1,z1,x2,y2,z2

  ! Translate the given point by the negatives of our coordinate system
  ! origin.
  X = DpointIn(1) - rcoordSys%Dorigin(1)
  Y = DpointIn(2) - rcoordSys%Dorigin(2)
  Z = DpointIn(3) - rcoordSys%Dorigin(3)

  ! Now scale the point by the inverse of our coordinate system sclaing
  ! factor, of course only if it is non-zero ...
  if (rcoordSys%dscalingFactor .ne. 0.0_DP) then

    X = X / rcoordSys%dscalingFactor
    Y = Y / rcoordSys%dscalingFactor
    Z = Z / rcoordSys%dscalingFactor

  else

    ! ... or set the result to zero and exit subroutine otherwise.
    DpointOut(1) = 0.0_DP
    DpointOut(2) = 0.0_DP
    DpointOut(3) = 0.0_DP
    return

  endif

  ! Check if the rotation of the coordinate system is non-zero.
  if (rcoordSys%drotationX .ne. 0.0_DP) then

    ! setup
    sx = sin(-rcoordSys%drotationX)
    cx = cos(-rcoordSys%drotationX)
    sy = sin(-rcoordSys%drotationY)
    cy = cos(-rcoordSys%drotationY)
    sz = sin(-rcoordSys%drotationZ)
    cz = cos(-rcoordSys%drotationZ)

    x1 =  X * cz + Y * sz
    y1 =  Y * cz - X * sz
    z1 =  Z
    x2 = x1 * cy + z1 * sy
    y2 = z1
    z2 = z1 * cy - x1 * sy
    DpointOut(1) = x2
    DpointOut(2) = y2 * cx + z1 * sx
    DpointOut(3) = z2 * cx - x1 * sx

  else

    ! No rotation in the coordinate system, so simply copy the input coords.
    DpointOut(1) = X
    DpointOut(2) = Y
    DpointOut(3) = Z

  end if

  end subroutine

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

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


! { dg-do compile }
! { dg-options "-fdec-structure" }
!
! Test a regression where union components could compare equal to structure/map
! components, causing an ICE in gfc_conv_component_ref.
!

implicit none

structure /s1/
  integer(4) i
end structure

structure /s2/
  union
    map
      record /s1/ r
    end map
  end union
end structure

record /s2/ x

x.r.i = 0

end


text:
!-------------------------------------------------------------------------------
!  The @header, @table_section, @table_subsection, @item and @end_table commands
!  are custom defined commands in Doxygen.in. They are defined under ALIASES.
!  For the page created here, the 80 column limit is exceeded. Arguments of
!  aliases are separated by ','. If you intended ',' to be a string you must use
!  an escaped comma '\,'.
!
!>  @page emission_file_sec Soft X-Ray Emission File Manual
!>
!>  @tableofcontents
!>  @section emission_file_intro_sec Introduction
!>  This page documents the NetCDF file format for the Soft X-Ray emission file.
!>  This file contains all the information to compute the Soft X-Ray emission
!>  from physical parameters. To allow for mutiple cameras and  multi-color
!>  systems, transmission and absorption functions need to be defined for each
!>  camera.
!>
!>  @section emission_file_contents_sec File Contents
!>  @header{Variable, Description, Dimension}
!>  @begin_table
!>     @item{te_start,   Initial temperature.,          Scalar}
!>     @item{te_step,    Temperature step.,             Scalar}
!>     @item{emissivity, Soft X-Ray emission function., (Number Cameras\, Number Energies)}
!>  @end_table
!>
!*******************************************************************************
!>  @file emission.f
!>  @brief Contains module @ref emission.
!
!  Note separating the Doxygen comment block here so detailed decription is
!  found in the Module not the file.
!
!>  Defines the base class of the type @ref emission_class. This contains the
!>  X-Ray emission as function of temperature and energy. This needs to still be
!>  multiplied by ne^2.
!*******************************************************************************

      MODULE emission
      USE stel_kinds
      USE profiler
      USE mpi_inc

      IMPLICIT NONE

!*******************************************************************************
!  DERIVED-TYPE DECLARATIONS
!  1) emission class
!
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  Base class representing the soft x-ray emission function.
!-------------------------------------------------------------------------------
      TYPE emission_class
!>  Initial temperature.
         REAL (rprec)                          :: te_start = 0.0
!>  Temperature step.
         REAL (rprec)                          :: te_step = 0.0

!>  X-Ray Emission function.
         REAL (rprec), DIMENSION(:,:), POINTER :: emissivity => null()
      END TYPE

!*******************************************************************************
!  INTERFACE BLOCKS
!*******************************************************************************
      INTERFACE emission_construct
         MODULE PROCEDURE emission_construct_netcdf
      END INTERFACE

      CONTAINS
!*******************************************************************************
!  CONSTRUCTION SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  @brief Construct a @ref emission_class object.
!>
!>  Allocates memory and initialized a @ref emission_class object from a netcdf
!>  file.
!>
!>  @param[in] filename File path to the emissivity netcdf file.
!>  @returns A pointer to a constructed @ref emission_class object.
!-------------------------------------------------------------------------------
      FUNCTION emission_construct_netcdf(filename)
      USE ezcdf
      USE v3_utilities

      IMPLICIT NONE

!  Declare Arguments
      TYPE (emission_class), POINTER :: emission_construct_netcdf
      CHARACTER (len=*), INTENT(in)  :: filename

!  local variables
      REAL (rprec)                   :: start_time
      INTEGER                        :: iou
      INTEGER                        :: varid
      INTEGER                        :: status
      INTEGER, DIMENSION(2)          :: dim_lengths

!  Start of executable code
      start_time = profiler_get_start_time()

      ALLOCATE(emission_construct_netcdf)

      status = nf90_open(TRIM(filename), NF90_NOWRITE, iou)
      CALL assert_eq(0, status, 'failed to open emission file')

      status = nf90_inq_varid(iou, 'te_start', varid)
      status = nf90_get_var(iou, varid,                                        &
     &                      emission_construct_netcdf%te_start)
      status = nf90_inq_varid(iou, 'te_step', varid)
      status = nf90_get_var(iou, varid,                                        &
     &                      emission_construct_netcdf%te_step)

!  Start the dimids in the dim_lengths then overwrite them. It is hard coded
!  that emissivity is two dimensional.
      status = nf90_inq_varid(iou, 'emissivity', varid)
      status = nf90_inquire_variable(iou, varid, dimids=dim_lengths)
      status = nf90_inquire_dimension(iou, dim_lengths(1),                     &
     &                                len=dim_lengths(1))
      status = nf90_inquire_dimension(iou, dim_lengths(2),                     &
     &                                len=dim_lengths(2))
      ALLOCATE(emission_construct_netcdf%emissivity(dim_lengths(1),            &
     &                                              dim_lengths(2)))
      status = nf90_get_var(iou, varid,                                        &
     &                      emission_construct_netcdf%emissivity)

      status = nf90_close(iou)

      CALL profiler_set_stop_time('emission_construct_netcdf',                 &
     &                            start_time)

      END FUNCTION

!*******************************************************************************
!  DESTRUCTION SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  @brief Deconstruct a @ref emission_class object.
!>
!>  Deallocates memory and uninitializes a @ref emission_class object.
!>
!>  @param[inout] this A @ref emission_class instance.
!-------------------------------------------------------------------------------
      SUBROUTINE emission_destruct(this)

      IMPLICIT NONE

!  Declare Arguments
      TYPE (emission_class), POINTER :: this

!  Start of executable code
      IF (ASSOCIATED(this%emissivity)) THEN
         DEALLOCATE(this%emissivity)
         this%emissivity => null()
      END IF

      DEALLOCATE(this)

      END SUBROUTINE

!*******************************************************************************
!  GETTER SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  @brief Gets the emission as a function of energy for a fixed temperature.
!>
!>  The emission for a fixed temperature is interpolated from a the temperature
!>  grid points. The emission array columns are the temperature and rows are the
!>  energy. Transmission and absorption function convert to the signal
!>
!>  @param[in] this  A @ref emission_class instance.
!>  @param[in] te    The electron temperature.
!>  @param[in] ne    The electron density.
!>  @param[in] index Filter index of the camera array.
!>  @returns The emissivity function for a fixed temperature.
!-------------------------------------------------------------------------------
      FUNCTION emission_get_emission(this, te, ne, index)

      IMPLICIT NONE

!  Declare Arguments
      REAL (rprec)                      :: emission_get_emission
      TYPE (emission_class), INTENT(in) :: this
      REAL (rprec), INTENT(in)          :: te
      REAL (rprec), INTENT(in)          :: ne
      INTEGER, INTENT(in)               :: index

!  local variables
      INTEGER                           :: i_low, i_high
      REAL (rprec)                      :: w_low, w_high
      INTEGER                           :: num_temp
      REAL (rprec)                      :: start_time

!  Start of executable code
      start_time = profiler_get_start_time()

!  Avoid a divid by zero error.
      IF (te .eq. 0.0) THEN
         emission_get_emission = 0.0
         CALL profiler_set_stop_time('emission_get_emission',                  &
     &                               start_time)
         RETURN
      END IF

!  Determine the high and low indices for interpolating. If the temperature is
!  beyond the first or last temperature, linearly interpolate from the last two
!  temperatures.
      num_temp = SIZE(this%emissivity, 2)

      IF (te .lt. this%te_start) THEN
!  Assume no emission at zero eV. Interpolate between zero and the first
!  temperature level.
         i_low = 1
         i_high = 2

         w_high = 0.0
         w_low = te/this%te_start
      ELSE IF (te .gt. this%te_start + (num_temp - 1)*this%te_step) THEN
         i_high = num_temp
         i_low = i_high - 1

         w_high = 1.0
         w_low = 0.0
      ELSE
         i_low = (te - this%te_start)/this%te_step + 1
         i_high = i_low + 1

         w_high = (te - this%te_start)/this%te_step + 1.0 - i_low
         w_low = 1.0 - w_high
      END IF

      emission_get_emission = w_high*this%emissivity(index,i_high)             &
     &                      + w_low*this%emissivity(index,i_low)

      emission_get_emission = ne*ne/SQRT(te)*MAX(emission_get_emission,        &
     &                                           0.0)

      CALL profiler_set_stop_time('emission_get_emission', start_time)

      END FUNCTION

      END MODULE


      module material_behavior

      use linearsolvers

!     =========================element center point
      real*8 :: el_center_coord(dimen)
      real*8 :: qtransform(dimen,dimen) ! this is the transformtion matrix
!     ================================================================
!                      history variables
!     ================================================================
!      real*8,allocatable::gauss_coords(:,:) !gauss_coords(igauss,dimen)
      real*8,allocatable::ele_pol_vector(:,:)  !gauss_coords(ele,dimen)
!      real*8,allocatable::shape_function_weight(:,:) ! shape_function_weight(nem,igauss_element)
!     ================================================================
!                      history variables
!     ================================================================
      real*8,allocatable::q_t(:),q_t_p(:),q_t_b(:)
      real*8,allocatable::                      &
        q_sigma_c_t(:,:,:),q_sigma_e_t(:,:,:),  &
        q_sigma_c_p(:,:,:),q_sigma_e_p(:,:,:),   &
        q_sigma_c_b(:,:,:),q_sigma_e_b(:,:,:),    &
        q_d_e_t(:,:),q_d_k_t(:,:),                 &
        q_d_e_p(:,:),q_d_k_p(:,:),                  &
        q_d_e_b(:,:),q_d_k_b(:,:)
!     ================================================================
!                      polarization history variables
!     ================================================================
      real*8,allocatable::      &
        q_rempol_t(:,:),         &
        q_rempol_p(:,:),          &
        q_rempol_b(:,:)

!     =======================mechanical properties variables
      real*8::cmate_voight(9,9)
      real*8::ctens(3,3,3,3)
      real*8::ctenst(3,3,3,3)
      real*8::cmat_voight(6,6)
      real*8::stense(3,3,3,3)
      real*8::e1,e2,e3,nu12,nu13,nu23,g12,g13,g23
      real*8::lambda,nu,mu,kbulk
!     =======================electric variables
      real*8,parameter::eps_vac_permit=8.58e-12
      real*8::epz(3,3,3)
      real*8,dimension(3,3)::ktense,ktense_trans
      real*8::ktense_inv(3,3) !the material matrix
      real*8::xeta(3,3,3)! the non linear electric permibility tensor

!     ===========================materials tensors
      real*8,dimension(3,3,3):: epz_t,epz_trans,epz_nl !the electro coupling tensor
      real*8,dimension(3,3,3,3)      :: b_tilt,b_tilt_trans ! the nonlinear electro coupling tensor

!     ===================================remanent polarization
      real*8,dimension(dimen):: rem_pol_t,drem_pol_t
      real*8,dimension(dimen):: rem_pol_p,drem_pol_p
      real*8,dimension(dimen):: rem_pol_b
!     ================================================================
!                      polarization history variables local
!     ================================================================
      real*8,dimension(dimen)::  &
        q_rempol_t_local,        &
        q_rempol_p_local,        &
        q_rempol_b_local
!     =======================polarization switching variables
      real*8::e_coercive ! coercive field
      integer::ipol
      real*8::p_per_sat
      real*8::pol_rate
      real*8::der_total(3)
!     ================================================================
!                      time variables
!     ================================================================
      real*8::time(2),dtime
      integer::ntime,itime  ! time
!     ================================================================
      integer::noelem
!     ================================================================
!     ================================================================
!     ================================================================
!     ================================================================
!     ================================================================
      contains

      subroutine material_properties()
      implicit none
!     ================================================================
!                      material variables
!     ================================================================

!     ================================================================
      integer::eye(dimen,dimen)
      integer::i,j,k,l,m,n

!     =======================polarization switching variables
      real*8 ::ecoup_norma
      real*8 ::ecoup_perpe
      real*8 ::ecoup_shear
!     =======================electric polarization (displacement in refrence confiquaration)
      real*8 :: polarization_vec(dimen)
      real*8 :: der_total(dimen),pol(3)
      real*8 :: der_unit(3)
      real*8,dimension(dimen)::  elr_t,elr_p,elr_b
      real*8 :: alpha_de(3,3)
      real*8 :: pol_intesity ,sat_pol
      real*8 :: kappa_0
!     ================================================================
      real*8::ey
      real*8::k00,k01,labda01
      real*8::labda_e_01,k_e_00,k_e_01
      real*8::labda_k_01,k_k_00,k_k_01
      common/visco_elastic/k00,k01,labda01,            &
                          k_e_00,k_e_01,labda_e_01,     &
                          k_k_00,k_k_01,labda_k_01
!     ================================================================
      eye = 0.0; do i = 1,dimen; eye(i,i) = 1;enddo

!      e_coercive=0.67e6  ![v/m] coercive field
!      pol_rate=5.4   !remnant polarization rate
!      p_per_sat= 0.26 ! [c/m^2]
!      kappa_0=0.04e-6
!!
!      e1= 6.06e10; ; ![pa]
!      nu23=0.3;
!     ==============================general polarization
      pol_intesity=norm_vect(polarization_vec);

      e1= 1 ![n/m^2] pzt 5a;
!      e1= 62e9 ![n/m^2] pzt 5a;
      nu23=0.3;

      ktense=0.0d0
      ktense(1,1)=(4.0d0)*1e-9;
      ktense(2,2)=ktense(1,1)
      ktense(3,3)=ktense(1,1)/2.0;

      ecoup_norma =-44.69;   ![n/(v.m)]
      ecoup_perpe =18.93;    ![n/(v.m)]
      ecoup_shear =-23.27;   ![n/(v.m)]

!      k00=1.0 ; k01=0.0; labda01=0.0
!      k_e_00=1.0 ; k_e_01=1.0; labda_e_01=1.0 ;
!      k_k_00=1.0 ; k_k_01=0.0; labda_k_01=0.0 ;

      e2=e1;e3=e2
      nu13=nu23;nu12=nu13;
      g23 = e1/(1+nu12)/2.0;
      g13 = g23;g12=g23;
      ey=e1
      nu=nu12

      mu=ey/(1+nu)/2.0
      lambda=ey*nu/(1+nu)/(1-2.0*nu)
      kbulk=lambda+2.0*mu/3.0

      der_unit=unit_vect(polarization_vec)
!     ==============================test
       der_unit=0.0d0
       der_unit(1)=1.0d0
!     ==============================
      do i=1,dimen
      do j=1,dimen
      alpha_de(i,j)=eye(i,j)-der_unit(i)*der_unit(j)
      enddo;enddo

      epz=0.0d0
      do i=1,dimen
      do j=1,dimen
      do k=1,dimen
      epz(k,i,j)=ecoup_norma*der_unit(k)*der_unit(i)*der_unit(j)  &
               +ecoup_perpe*der_unit(k)*alpha_de(i,j)             &
               +ecoup_shear*der_unit(i)*alpha_de(j,k)*0.5         &
               +ecoup_shear*der_unit(j)*alpha_de(i,k)*0.5
      enddo;enddo;enddo

     ktense_inv=0.0d0
     forall(k=1:3) ktense_inv(k,k)=1/ktense(k,k)

      xeta=0.0d0;
      b_tilt=0.0d0

      do i=1,dimen;do j=1,dimen;do k=1,dimen;do l=1,dimen
      ctens(i,j,k,l)=lambda*eye(i,j)*eye(k,l)+mu*(eye(i,k)*eye(j,l)+eye(i,l)*eye(j,k))
      enddo;enddo;enddo;enddo

      do i=1,dimen;do j=1,dimen;do k=1,dimen;do l=1,dimen
      stense(i,j,k,l)=(1.0/kbulk/9.0)*eye(i,j)*eye(k,l)+(1/mu/4.0)*(eye(i,k)*eye(j,l)+eye(i,l)*eye(j,k)-(2.0/3.0)*eye(i,j)*eye(k,l))
      enddo;enddo;enddo;enddo
!     ==============================x3 polarized
!     epz=0.0d0
!     epz(3,3,3)=45.744;
!     epz(3,1,1)=-19.375;
!     epz(1,1,3)=23.81
!
!     epz(3,2,2)=epz(3,1,1)
!
!     epz(2,2,3)=epz(1,1,3)
!     epz(1,3,1)=epz(1,1,3)
!     epz(2,3,2)=epz(1,1,3)
!     b_tilt(3,3,3,3)=-3.35e-2;
!     b_tilt(3,3,1,1)=3.35e-2;
!     b_tilt(3,3,2,2)=3.35e-2;
!     ==============================x3 polarized
      ktense=0.0d0
      ktense(1,1)=(4.0d0)*1e-9;
      ktense(2,2)=ktense(1,1)
      ktense(3,3)=ktense(1,1)/2.0;

!      call material_transofrmation_tensor_x2()
!      call tensor_transform_rank_3(qtransform,epz,epz_trans)
!      call tensor_transform_rank_4(qtransform,b_tilt,b_tilt_trans)
!      call tensor_transform_rank_2(qtransform,ktense,ktense_trans)
!
!      epz=epz_trans
!      b_tilt=b_tilt_trans
!      ktense=ktense_trans

      epz_t=0.0;
      epz_nl=0.0;

      do 50 i=1,dimen; do 50 j=1,dimen;do 50 k=1,dimen;do 50 l=1,dimen
      epz_nl(k,i,j)=epz(k,i,j)+0.5*b_tilt(k,l,i,j)*abs(elr_t(l))
50    epz_t(k,i,j)=epz(k,i,j)+b_tilt(k,l,i,j)*abs(elr_t(l))

!    1D formulation
ctens=0; ctens(1,1,1,1)=1;
epz=0; epz(1,1,1)=-1;
ctenst=ctens;
epz_t=epz;
ktense=0.0d0  ;ktense(1,1)=1;

      endsubroutine material_properties

      subroutine  update_history(ngauss,dimen)
      integer::ngauss,dimen

      q_sigma_c_b=q_sigma_c_p;
      q_sigma_c_p=q_sigma_c_t;

      q_sigma_e_b=q_sigma_e_p;
      q_sigma_e_p=q_sigma_e_t;

      q_d_e_b=q_d_e_p
      q_d_e_p=q_d_e_t

      q_d_k_b=q_d_k_p
      q_d_k_p=q_d_k_t

      q_rempol_b=q_rempol_p
      q_rempol_p=q_rempol_t

      endsubroutine  update_history

!     ===============================================================
!           forming history variables
!     ===============================================================
      subroutine  form_history(ngauss,dimen,ipdf,nem)
      integer::ngauss,dimen,ipdf,nem
!     ===================allocate history variables
      allocate(                                                                  &
      q_sigma_c_t(dimen,dimen,ngauss),q_sigma_e_t(dimen,dimen,ngauss),            &
      q_sigma_c_p(dimen,dimen,ngauss),q_sigma_e_p(dimen,dimen,ngauss),            &
      q_sigma_c_b(dimen,dimen,ngauss),q_sigma_e_b(dimen,dimen,ngauss),             &
        q_d_e_t(dimen,ngauss),q_d_k_t(dimen,ngauss),                                &
        q_d_e_p(dimen,ngauss),q_d_k_p(dimen,ngauss),                                 &
        q_d_e_b(dimen,ngauss),q_d_k_b(dimen,ngauss),                                  &
      q_t(ngauss),q_t_p(ngauss),q_t_b(ngauss),                                         &
      q_rempol_t(dimen,ngauss),                                                         &
      q_rempol_p(dimen,ngauss),                                                          &
      q_rempol_b(dimen,ngauss),          &
      ele_pol_vector(nem,dimen))

      q_sigma_c_t=0.0d0;q_sigma_e_t=0.0d0;
      q_sigma_c_p=0.0d0;q_sigma_e_p=0.0d0;
      q_sigma_c_b=0.0d0;q_sigma_e_b=0.0d0;

       q_d_e_t=0.0d0;q_d_k_t=0.0d0;
       q_d_e_p=0.0d0;q_d_k_p=0.0d0;
       q_d_e_b=0.0d0;q_d_k_b=0.0d0;

     q_rempol_t=0.0d0;
     q_rempol_p=0.0d0;
     q_rempol_b=0.0d0;
     ele_pol_vector=0.0d0

      q_t=0.0d0 ; q_t_p=0.0d0; q_t_b=0.0d0

      end subroutine  form_history

subroutine stress_elect_displacement(ur,der,sigma)
!     ================================================================
!                          input variables
!     ================================================================
real*8,intent(in) :: ur(:,:)! global coefficient
!     ================================================================
!                         output variables
!     ================================================================
integer :: i !integer counters
real*8  :: der(:),sigma(:,:)
!     ================================================================
real*8::strain(3,3);
real*8::electric_feild(3);
!     ================================================================
strain=0.0d0;electric_feild=0.0d0;
strain(1:dimen,1:dimen)=0.5d0*(   ur(1:dimen,:)+transpose(ur(1:dimen,:))+      &
                              matmul(  transpose(   ur(1:dimen,:)  ),ur(1:dimen,:))    );

electric_feild(1:dimen)=-ur(dimen+1,:)

!write(1,*)'strain=',strain(1,1)
!write(1,*)'electric_feild=',electric_feild(1)

sigma =  0.0d0 ;
der   =  0.0d0 ;

sigma(1,1)=ctens(1,1,1,1)*strain(1,1)-epz(1,1,1)*electric_feild(1)
der(1)=epz(1,1,1)*strain(1,1)+ktense(1,1)*electric_feild(1)

!write(1,*)'ktense(1,1)=',ktense(1,1)
!write(1,*)'sigma=',sigma(1,1)
!write(1,*)'der=',der(1)

endsubroutine stress_elect_displacement

      end module material_behavior


text:
program parabola
	real*8 t0,t1,t2,t,e,a1,a2
	t1=0.0;t0=1.0;t2=2.0;e=0.001;t=0.0;a1=0.0;a2=0.0
	open(14,file="Parabola")
	write(14,"(3(a13))")"t1","t0","t2"
	do
		write(14,"(f13.6,f13.6,f13.6)")t1,t0,t2
		a1=(t0**2.0-t2**2.0)*(8.0*t1**3-2.0*t1**2.0-7.0*t1+3.0)&
		+(t2**2.0-t1**2.0)*(8.0*t0**3.0-2.0*t0**2.0-7.0*t0+3.0)&
		+(t1**2.0-t0**2.0)*(8.0*t2**3.0-2.0*t2**2.0-7.0*t2+3.0)
		a1=a1/((t1-t0)*(t0-t2)*(t2-t1))
		a2=(t0-t2)*(8.0*t1**3-2.0*t1**2-7.0*t1+3.0)&
		+(t2-t1)*(8.0*t0**3-2.0*t0**2-7.0*t0+3.0)&
                +(t1-t0)*(8.0*t2**3-2.0*t2**2-7.0*t2+3)
		a2=a2/((t1-t0)*(t0-t2)*(t2-t1))
		t=a1/(2*a2)
		if((((t-t0)**2)**0.5)<e)exit
		if(t-t0)10,10,20
		10 if((8*t**3-2*t**2-7*t+3)-(8*t0**3-2*t0**2-7*t0+3))101,101,102
		20 if((8*t**3-2*t**2-7*t+3)-(8*t0**3-2*t0**2-7*t0+3))201,201,202
		101 t2=t0;t0=t; goto 30
		102 t1=t; goto 30
		201 t1=t0;t0=t; goto 30
		202 t2=t; goto 30
!程序编到此，继续编写两个if语句，
!记得拷地震学的电子书,以及还矩阵理论

		30 cycle
	end do
	write(14,"(f13.6,f13.6,f13.6,f13.6)")t1,t,t2,(8*t**3-2*t**2-7*t+3)
	close(14)
end program parabola

max_stars_repo_name: DearDon/Optimization

max_issues_repo_name: DearDon/Optimization

max_forks_repo_name: DearDon/Optimization


text:
C Copyright (C) 2002, Carnegie Mellon University and others.
C All Rights Reserved.
C This code is published under the Common Public License.

      subroutine MULT_N(TRANS, NZ, NY, NU, NP, NCOL, NAC, NELE, TI,
     1     NNZFU, NNZAU, NNZFPP, NNZCU, NNZCUU, IRN, JCN, N, VIN, VOUT)
C
C     $Id: mult_n.f,v 1.1 2002/05/02 18:52:15 andreasw Exp $
C
C     Multiply vector with overall nonbasis matrix
C
C     Author:  Andreas Waechter
C              c/o Group of Larry Biegler
C              Department of Chemical Engineering
C              Carnegie Mellon University
C              Pittsburgh, PA
C
      implicit none

      character*(*) TRANS       !if = 't' or 'T', compute vout += CW^T*vin
                                !otherwise vout += CW*vin
                                !NOTE: VOUT is NOT INITIALIZED here!!!
      integer NZ, NY, NU, NP, NCOL, NAC, NELE
      double precision TI(NELE+1)
      integer NNZFU
      integer NNZAU
      integer NNZFPP
      integer NNZCU
      integer NNZCUU(NCOL)
C      integer IRN(NNZFU+NNZAU+NCOL*NNZFPP+NNZCU)
C      integer JCN(NNZFU+NNZAU+NCOL*NNZFPP+NNZCU)
      integer IRN(*)
      integer JCN(*)
C      double precision N(NELE*(NNZFU+NNZAU+NCOL*NNZFPP))
      double precision N(*)
      double precision VIN(*)
      double precision VOUT(*)

      integer lfu, lfp, lcu, lout, lin, ln, lu, lp, iele
      double precision h

      lfu = 1
      lfp = lfu + NNZFU + NNZAU
      lcu = lfp + NCOL*NNZFPP

      if( TRANS(1:1).eq.'t' .or. TRANS(1:1).eq.'T' ) then

         lin = 1 + NZ
         lu  = 1
         lp  = 1 + NELE*(NCOL*NU-NAC)
         ln  = 1

         call DCOPY(NELE*(NCOL*NU-NAC) + NP, 0d0, 0, VOUT, 1)

         do iele = 1, NELE
C
C     FU^T  and  AU^T
C
            call MULTMAT(TRANS, .true., NNZFU+NNZAU, IRN(lfu), JCN(lfu),
     1           N(ln), VIN(lin), VOUT(lu))
            ln  = ln  + NNZFU + NNZAU
C
C     FP^T
C
            call MULTMAT(TRANS, .true., NCOL*NNZFPP, IRN(lfp), JCN(lfp),
     1           N(ln), VIN(lin), VOUT(lp))
            ln  = ln  + NCOL*NNZFPP
            lin = lin + NCOL*(NZ+NY) + NAC
C
C     CU^T
C
            h = TI(iele+1) - TI(iele)
            call MULTCW(TRANS, .true., NZ, NY, NU, NP, NCOL, h,
     1           NNZCU, NNZCUU, IRN(lcu), JCN(lcu), VIN(lin), VOUT(lu))
C
            lin = lin + NZ
            lu  = lu  + NCOL*NU - NAC

         enddo

      else

         call DCOPY(NZ + NELE*(NCOL*(NZ+NY) + NZ + NAC), 0d0, 0,
     1        VOUT, 1)
C
C     loop over all elements
C
         lout = 1 + NZ
         lu   = 1
         lp   = 1 + NELE*(NCOL*NU-NAC)
         ln   = 1

         do iele = 1, NELE
C
C     FU  and  AU
C
            call MULTMAT(TRANS, .true., NNZFU+NNZAU, IRN(lfu), JCN(lfu),
     1           N(ln), VIN(lu), VOUT(lout))
            ln   = ln   + NNZFU + NNZAU
C
C     FP
C
            call MULTMAT(TRANS, .true., NCOL*NNZFPP, IRN(lfp), JCN(lfp),
     1           N(ln), VIN(lp), VOUT(lout))
            ln   = ln   + NCOL*NNZFPP
            lout = lout + NCOL*(NZ+NY) + NAC
C
C     Continuity equations
C
            h = TI(iele+1) - TI(iele)
            call MULTCW(TRANS, .true., NZ, NY, NU, NP, NCOL, h,
     1           NNZCU, NNZCUU, IRN(lcu), JCN(lcu), VIN(lu), VOUT(lout))
C
            lout = lout + NZ
            lu   = lu   + NCOL*NU - NAC

         enddo

      endif

      return
      end

max_stars_repo_name: kant/CMMPPT

max_issues_repo_name: kant/CMMPPT

max_forks_repo_name: kant/CMMPPT


text:
*     -----------------
*     PROGRAMME EXAMPLE
*     -----------------
*
*     Reference : Bureau des Longitudes - PBGF9502.
*
*     ------
*     Object
*     ------
*
*     This program is an example of use of VSOP87 subroutine.
*     It computes the coordinates of a planet for a version of VSOP87.
*     VSOP87 (main version) : elliptic coordinates J2000
*     VSOP87A : heliocentric rectangular coordinates J2000
*     VSOP87B : heliocentric spherical   coordinates J2000
*     VSOP87C : heliocentric rectangular coordinates of date
*     VSOP87D : heliocentric spherical   coordinates of date
*     VSOP87E : barycentric  rectangular coordinates J2000
*     The frame is referred to the dynamical equinox and ecliptic.
*     The time scale is the dynamical barycentric time TDB.
*
*     ------------
*     Declarations
*     ------------
*
      implicit double precision (a-h,o-z)
      logical fexist
      integer day,month,year,d0,h0,dt,datmod
      character rep,a(0:5)
      character*4 ext(0:9)
      character*6 prefix
      character*7 vers
      character*21 body(0:9)
      character*40 fich
      character*37 coord(4)
      character*39 reper(2)
      dimension r(6),kbody(0:8,0:5),maxbody(0:5),prec0(0:9)
      data a/' ','A','B','C','D','E'/
      data ext/'.sun','.mer','.ven','.ear','.mar','.jup',
     .         '.sat','.ura','.nep','.emb'/
      data body/'SUN','MERCURY','VENUS','EARTH','MARS','JUPITER',
     .          'SATURN','URANUS','NEPTUNE','EARTH-MOON BARYCENTER'/
      data coord/'Heliocentric elliptic coordinates.',
     .           'Heliocentric rectangular coordinates.',
     .           'Heliocentric spherical coordinates.',
     .           'Barycentric rectangular coordinates.'/
      data reper/'Dynamical equinox and ecliptic J2000.',
     .           'Dynamical equinox and ecliptic of date.'/
      data prec0/1.d-8,0.6d-8,2.5d-8,2.5d-8,10.d-8,35.d-8,70.d-8,
     .           8.d-8,42.d-8,2.5d-8/
      data kbody/1,2,4,5,6,7,8,9,10,
     .           1,2,3,4,5,6,7,8,9,
     .           1,2,3,4,5,6,7,8,10,
     .           1,2,3,4,5,6,7,8,10,
     .           1,2,3,4,5,6,7,8,10,
     .           0,1,2,3,4,5,6,7,8/
      data maxbody/7,8,7,7,7,8/
      data prefix/'VSOP87'/
      menu=4
      lu0=0
*
*     ----------------
*     Version and body
*     ----------------
*
100   continue
      if (menu.ne.1.and.menu.ne.4) goto 200
      call CLRSCR
      print 1000
      print 2001
      read (*,3001,iostat=nerr) ivers
      if (nerr.ne.0) goto 100
      if (ivers.lt.0.or.ivers.gt.5) goto 100
      vers='VSOP87'//a(ivers)
      call CLRSCR
      print 1000
      print 2002, vers
      do i=0,maxbody(ivers)
         nbody=kbody(i,ivers)
         print 2003, nbody,body(nbody)
      enddo
      print 2004
      read (*,3004,iostat=nerr) ibody
      if (nerr.ne.0) goto 100
      irep=0
      do i=0,maxbody(ivers)
         if (kbody(i,ivers).eq.ibody) irep=1
      enddo
      if (irep.eq.0) goto 100
      print 1001
      read 3000, rep
      if (rep.eq.'N'.or.rep.eq.'n') goto 100
*
*     ----------------
*     Open file VSOP87
*     ----------------
*
      lu=ivers*10+10+ibody
      if (lu0.ne.0.and.lu0.ne.lu) close (lu0)
      if (ivers.eq.0) fich=prefix//ext(ibody)
      if (ivers.ne.0) fich=prefix//a(ivers)//ext(ibody)
      inquire (file=fich,exist=fexist)
      if (fexist) then
         if (lu.ne.lu0) then
            open (lu,file=fich,status='old',iostat=nerr)
            if (nerr.ne.0) goto 610
            lu0=lu
         endif
      else
         goto 620
      endif
*
*     -------------
*     Datation mode
*     -------------
*
200   continue
      if (menu.ne.2.and.menu.ne.4) goto 300
      call CLRSCR
      print 1000
      print 2005
      read (*,3005,iostat=nerr) datmod
      if (nerr.ne.0) goto 200
      if (datmod.lt.1.or.datmod.gt.2) goto 200
      if (datmod.eq.1) then
         print 2006
         read (*,3006,iostat=nerr) t0
         if (nerr.ne.0) goto 200
         if (t0.lt.1.d0) goto 200
      else
         print 2007
         read (*,3007,iostat=nerr) d0
         if (nerr.ne.0) goto 200
         call DATEJD (d0,0,td)
         if (td.eq.0.d0) goto 200
         print 2008
         read (*,3008,iostat=nerr) h0
         if (nerr.ne.0) goto 200
         call DATEJD (d0,h0,t0)
         if (t0.eq.0.d0) goto 200
      endif
      t=t0
      print 2009
      read (*,3009,iostat=nerr) dt
      if (nerr.ne.0.or.dt.lt.0) goto 200
      nj=dt/1000000
      nh=mod(dt,1000000)/10000
      nm=mod(dt,10000)/100
      ns=mod(dt,100)
      if (nj.gt.100.or.nh.gt.24.or.nm.gt.60.or.ns.gt.60) goto 200
      pas=nj+nh/24.d0+nm/1440.d0+ns/86400.d0
      print 1001
      read 3000, rep
      if (rep.eq.'N'.or.rep.eq.'n') goto 200
*
*     --------------------------
*     Precision level (relative)
*     --------------------------
*
300   continue
      if (menu.ne.3.and.menu.ne.4) goto 400
      call CLRSCR
      print 1000
      print 2010
      read (*,3010,iostat=nerr) iprec
      if (nerr.ne.0) goto 300
      if (iprec.lt.1.or.iprec.gt.2) goto 300
      if (iprec.eq.1) then
         prec=0.d0
      else
         print 2011, prec0(ibody)
         read (*,3011,iostat=nerr) prec
         if (nerr.ne.0) goto 300
         if (prec.lt.prec0(ibody).or.prec.gt.1.d-2) goto 300
         print 1001
         read 3000, rep
         if (rep.eq.'N'.or.rep.eq.'n') goto 300
      endif
*
*     ----------------------------
*     Computation of the ephemeris
*     ----------------------------
*
400   continue
      print 2012
*
      call VSOP87 (t,ivers,ibody,prec,lu,r,ierr)
*
      if (ierr.ne.0) goto 600
      call JDDATE (t,idate,iheur)
      jdate=abs(idate)
      day=mod(jdate,100)
      month=mod(jdate/100,100)
      year=idate/10000
      nh=iheur/10000
      nm=mod(iheur/100,100)
      ns=mod(iheur,100)
      ic=1
      ir=1
      if (ivers.eq.1.or.ivers.eq.3) ic=2
      if (ivers.eq.2.or.ivers.eq.4) ic=3
      if (ivers.eq.5) ic=4
      if (ivers.eq.3.or.ivers.eq.4) ir=2
*
410   continue
      call CLRSCR
      if (prec.lt.prec0(ibody)) then
         print 4001, vers,body(ibody),coord(ic),reper(ir),
     .               day,month,year,t,nh,nm,ns
      else
         print 4002, vers,body(ibody),coord(ic),reper(ir),
     .               day,month,year,t,nh,nm,ns,prec
      endif
      if (ivers.eq.0) print 4003,r
      if (ivers.eq.1.or.ivers.eq.3.or.ivers.eq.5)
     .   print 4004, r(1),r(4),r(2),r(5),r(3),r(6)
      if (ivers.eq.2.or.ivers.eq.4)
     .   print 4005, r(1),r(4),r(2),r(5),r(3),r(6)
      print 2013
      read 3000, rep
      if (rep.ne.'+'.and.rep.ne.'-') goto 500
      if (dt.eq.0) goto 410
      if (rep.eq.'+') t=t+pas
      if (rep.eq.'-') t=t-pas
      goto 400
*
*     -------------------
*     Another computation
*     -------------------
*
500   continue
      call CLRSCR
      print 1000
      print 2014
      read (*,3014,iostat=nerr) menu
      if (nerr.ne.0) goto 500
      if (menu.lt.0.or.menu.gt.4) goto 500
      if (menu.ne.0) goto 100
      call CLRSCR
      print 1000
      print 2015
      pause ' Program terminated : Hit Enter'
      stop
*
*     ---------------
*     Errors Messages
*     ---------------
*
600   continue
      call CLRSCR
      print 1000
      if (ierr.eq.1) print 5001, ivers
      if (ierr.eq.2) print 5002, ibody
      if (ierr.eq.3) print 5003, prec
      if (ierr.eq.4) print 5004, fich
      pause ' Program terminated : Hit Enter'
      stop
610   continue
      call CLRSCR
      print 1000
      print 5005, fich
      pause ' Program terminated : Hit Enter'
      stop
620   continue
      call CLRSCR
      print 1000
      print 5006, fich
      pause ' Program terminated : Hit Enter'
      stop
*
*     -------
*     Formats
*     -------
*
1000  format (2x,'PLANETARY SOLUTION VSOP87'/
     .        2x,'-------------------------')
1001  format (/2x,'Confirm your choice (y/n) ? ',$)
*
2001  format (/2x,'Choice of the version'/
     .        /2x,'0 : VSOP87 (main solution), elliptic coordinates'
     .        /2x,'    dynamical equinox and ecliptic J2000'/
     .        /2x,'1 : VSOP87A, heliocentric rectangular coordinates'
     .        /2x,'    dynamical equinox and ecliptic J2000'/
     .        /2x,'2 : VSOP87B, heliocentric spherical coordinates'
     .        /2x,'    dynamical equinox and ecliptic J2000'/
     .        /2x,'3 : VSOP87C, heliocentric rectangular coordinates'
     .        /2x,'    dynamical equinox and ecliptic of the date'/
     .        /2x,'4 : VSOP87D, heliocentric spherical coordinates'
     .        /2x,'    dynamical equinox and ecliptic of the date'/
     .        /2x,'5 : VSOP87E, barycentric rectangular coordinates'
     .        /2x,'    dynamical equinox and ecliptic J2000'/
     .        /2x,'Enter your choice ? ',$)
2002  format (/2x,'Choice of the body in version ',a7/)
2003  format (2x,i1,' : ',a21)
2004  format (/2x,'Enter your choice ? ',$)
2005  format (/2x,'Datation mode'/
     .        /2x,'1 : computation with julian dates'
     .        /2x,'2 : computation with calendar dates'/
     .        /2x,'Enter your choice ? ',$)
2006  format (/2x,'Enter the initial julian date'
     .        /2x,'don''t forget the decimal point ---> ',$)
2007  format (/2x,'Enter the initial date'
     .        /2x,'format : +/-yyyymmdd -------------> ',$)
2008  format (/2x,'Enter the initial time'
     .        /2x,'format : hhmmss ------------------> ',$)
2009  format (/2x,'Enter the time step (<100 days)'
     .        /2x,'format : ddhhmmss ----------------> ',$)
2010  format (/2x,'Precision level'/
     .        /2x,'1 : computation with the complete solution'
     .        /2x,'2 : computation with a truncated solution'/
     .        /2x,'Enter your choice ? ',$)
2011  format (/2x,'Enter the relative precision level'
     .        /2x,'between ',f11.9,' and 0.01'
     .        /2x,'don''t forget the decimal point -----> ',$)
2012  format (/2x,'Wait Please . . .')
2013  format (/2x,77('-')//2x,'Choice of time :      next time -> +',
     .        '      previous time -> -      stop -> 0'/
     .        /2x,'Enter your choice ? ',$)
2014  format (/2x,'Another computation'/
     .        /2x,'0 : stop'/
     .        /2x,'1 : with new version or body'
     .        /2x,'2 : with new datation mode'
     .        /2x,'3 : with new precision level'/
     .        /2x,'4 : with new parameters'/
     .        /2x,'Enter your choice ? ',$)
2015  format (///)
*
3000  format (a)
3001  format (i2)
3004  format (i2)
3005  format (i2)
3006  format (f14.6)
3007  format (i10)
3008  format (i7)
3009  format (i9)
3010  format (i2)
3011  format (f12.9)
3014  format (i2)
*
4001  format (2x,'PLANETARY SOLUTION ',a7,5x,'BODY : ',a21
     .        /2x,77('-')//2x,a37,1x,a39/
     .        /2x,'Date : ',i2.2,' / ',i2.2,' / ',i5,
     .        t41,'Julian date : ',f14.6
     .        /2x,'Time : ',i2.2,'h ',i2.2,'m ',i2.2,'s TDB',
     .        t41,'Precision of the complete solution')
4002  format (2x,'PLANETARY SOLUTION ',a7,5x,'BODY : ',a21
     .        /2x,77('-')//2x,a37,1x,a39/
     .        /2x,'Date : ',i2.2,' / ',i2.2,' / ',i5,
     .        t41,'Julian date : ',f14.6
     .        /2x,'Time : ',i2.2,'h ',i2.2,'m ',i2.2,'s TDB',
     .        t41,'Relative precision : ',f10.9)
4003  format (/2x,'Semi-major axis :  ',f13.10,' au',
     .        t41,'Mean longitude :   ',f13.10,' rad'
     .        /2x,'k = e*cos(p) :     ',f13.10,' rad',
     .        t41,'h = e*sin(p) :     ',f13.10,' rad'
     .        /2x,'q = sin(g)*cos(G) :',f13.10,' rad',
     .        t41,'p = sin(g)*sin(G) :',f13.10,' rad'/
     .        /2x,'e : eccentricity',
     .        t41,'p : perihelion longitude'
     .        /2x,'g : semi-iclination',
     .        t41,'G : ascending node longitude')
4004  format (/2x,'X :',f14.10,' au',t41,'X'' :',f14.10,' ua/d'
     .        /2x,'Y :',f14.10,' au',t41,'Y'' :',f14.10,' ua/d'
     .        /2x,'Z :',f14.10,' au',t41,'Z'' :',f14.10,' ua/d')
4005  format (/2x,'Longitude :',f14.10,' rad',
     .        t41,'vitesse  :',f14.10,' rad/d'
     .        /2x,'Latitude : ',f14.10,' rad',
     .        t41,'vitesse  :',f14.10,' rad/d'
     .        /2x,'Radius :   ',f14.10,' au',
     .        t41,'vitesse  :',f14.10,' au/d')
*
5001  format (/2x,'*****  Error *****'/
     .        /2x,'Wrong version index : ',i1/)
5002  format (/2x,'*****  Error *****'/
     .        /2x,'Wrong body index : ',i1/)
5003  format (/2x,'*****  Error *****'/
     .        /2x,'Wrong precision level : ',d9.2/)
5004  format (/2x,'*****  Error *****'/
     .        /2x,'Fatal reading file error : ',a40/)
5005  format (/2x,'*****  Error *****'/
     .        /2x,'Fatal open file error : ',a40/)
5006  format (/2x,'***** File Error *****'/
     .        /2x,'Missing file : ',a40/)
*
      end
*
*
*
      subroutine JDDATE (tjd,idate,ihour)
*-----------------------------------------------------------------------
*
*     Object :
*     Conversion   Julian date ---> Calendar date    (Meeus formular).
*
*     Input :
*     tjd         julian date (real double precision).
*
*     Ouput :
*     idate       calendar date (integer).
*                 julian calendar before 1582 october 15
*                 gregorian calendar after.
*                 code: *yyyymmdd (* sign).
*     ihour       hour (integer).
*                 code: hhmmss.
*
*-----------------------------------------------------------------------
      implicit double precision (a-h,o-z)
      integer day,month,year
      idate=0
      ihour=0
      if (tjd.lt.0.d0) return
      t=tjd+0.5d0/86400.d0+0.5d0
      z=dint(t)
      f=t-z
      if (z.lt.2299161.d0) then
         a=z
      else
         x=dint((z-1867216.25d0)/36524.25d0)
         a=z+1.d0+x-dint(x/4.d0)
      endif
      b=a+1524.d0
      c=dint((b-122.1d0)/365.25d0)
      d=dint(365.25d0*c)
      e=dint((b-d)/30.6001d0)
      day=b-d-dint(30.6001d0*e)
      month=e-1.d0
      if (e.lt.13.5d0) then
         month=e-1.d0
      else
         month=e-13.d0
      endif
      if (month.lt.3) then
         year=c-4715.d0
      else
         year=c-4716.d0
      endif
      is=+1
      if (year.lt.0) is=-1
      idate=((iabs(year)*100+month)*100+day)*is
      f=f*24.d0
      ih=f
      f=(f-ih)*60.d0
      im=f
      f=(f-im)*60.d0
      is=f
      ihour=(ih*100+im)*100+is
      return
      end
*
*
*
      subroutine DATEJD (idate,ihour,tjd)
*-----------------------------------------------------------------------
*
*     Object :
*     Conversion   Calendar date -> Julian date    (Meeus formular).
*
*     Input :
*     idate       calendar date gregorian/julian (integer).
*                 julian calendar before 1582 october 15
*                 gregorian calendar after.
*                 code: *yyyymmdd (* sign).
*     ihour       hour (integer).
*                 code: hhmmss.
*
*     Output
*     tjd         julian date (real double precision).
*
*-----------------------------------------------------------------------
      implicit double precision (a-h,o-z)
      integer day,month,year
      dimension lm(12)
      data lm/31,28,31,30,31,30,31,31,30,31,30,31/
      tjd=0.d0
      year=idate/10000
      if (year.lt.-4713.or.year.gt.5000) return
      kdate=iabs(idate)-iabs(year)*10000
      month=kdate/100
      if (month.lt.1.or.month.gt.12) return
      day=kdate-month*100
      lm(2)=28
      if (mod(year,4).eq.0) then
         lm(2)=29
         if (year.gt.1582) then
            ncent=year/100
            if (mod(year,100).eq.0.and.mod(ncent,4).ne.0) lm(2)=28
         endif
      endif
      if (day.lt.1.or.day.gt.lm(month)) return
      is=ihour
      ih=ihour/10000
      if (ih.lt.0.or.ih.gt.24) return
      is=is-ih*10000
      im=is/100
      if (im.lt.0.or.im.gt.60) return
      is=is-im*100
      if (is.lt.0.or.is.gt.60) return
      a=0.d0
      b=0.d0
      c=0.d0
      if (month.gt.2) then
         y=year
         m=month
      else
         y=year-1
         m=month+12
      endif
      if (y.lt.0.d0) then
         c=-0.75d0
      else
         if (idate.ge.15821015) then
            a=dint(y/100.d0)
            b=2.d0-a+dint(a/4.d0)
         endif
      endif
      tjd=dint(365.25d0*y+c)+dint(30.6001d0*(m+1))+day+
     .    dfloat(ih)/24.d0+dfloat(im)/1440.d0+dfloat(is)/86400.d0+
     .    1720994.5d0+b
      return
      end
*
*
*
      subroutine CLRSCR
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9412
*
*     Object :
*     Clear the screen.
*
*     Remark :
*     *DOS  for DOS system.
*     *UNX  for UNIX system.
*
*-----------------------------------------------------------------------
*
*DOS  character*7 escscr
*DOS  escscr=char(27)//char(91)//'2J'//char(27)//char(91)//'H'
*DOS  write (*,'(2x,a)') escscr
*
*UNX  call system ('clear')
*
      call system ('clear')
*
      return
      end
*
*
*
      subroutine VSOP87 (tdj,ivers,ibody,prec,lu,r,ierr)
*-----------------------------------------------------------------------
*
*     Reference : Bureau des Longitudes - PBGF9502
*
*     Object :
*
*     Substitution of time in VSOP87 solution written on a file.
*     The file corresponds to a version of VSOP87 theory and to a body.
*
*     Input :
*
*     tdj      julian date (real double precision).
*              time scale : dynamical time TDB.
*
*     ivers    version index (integer).
*              0: VSOP87 (initial solution).
*                 elliptic coordinates
*                 dynamical equinox and ecliptic J2000.
*              1: VSOP87A.
*                 rectangular coordinates
*                 heliocentric positions and velocities
*                 dynamical equinox and ecliptic J2000.
*              2: VSOP87A.
*                 spherical coordinates
*                 heliocentric positions and velocities
*                 dynamical equinox and ecliptic J2000.
*              3: VSOP87C.
*                 rectangular coordinates
*                 heliocentric positions and velocities
*                 dynamical equinox and ecliptic of the date.
*              4: VSOP87D.
*                 spherical coordinates
*                 heliocentric positions and velocities
*                 dynamical equinox and ecliptic of the date.
*              5: VSOP87E.
*                 rectangular coordinates
*                 barycentric positions and velocities
*                 dynamical equinox and ecliptic J2000.
*
*     ibody    body index (integer).
*              0: Sun
*              1: Mercury
*              2: Venus
*              3: Earth
*              4: Mars
*              5: Jupiter
*              6: Saturn
*              7: Uranus
*              8: Neptune
*              9: Earth-Moon barycenter
*
*     prec     relative precision (real double precision).
*
*              if prec is equal to 0 then the precision is the precision
*                 p0 of the complete solution VSOP87.
*                 Mercury    p0 =  0.6 10**-8
*                 Venus      p0 =  2.5 10**-8
*                 Earth      p0 =  2.5 10**-8
*                 Mars       p0 = 10.0 10**-8
*                 Jupiter    p0 = 35.0 10**-8
*                 Saturn     p0 = 70.0 10**-8
*                 Uranus     p0 =  8.0 10**-8
*                 Neptune    p0 = 42.0 10**-8
*
*              if prec is not equal to 0, let us say in between p0 and
*              10**-2, the precision is :
*                 for the positions :
*                 - prec*a0 au for the distances.
*                 - prec rd for the other variables.
*                 for the velocities :
*                 - prec*a0 au/day for the distances.
*                 - prec rd/day for the other variables.
*                   a0 is semi-major axis of the body.
*                 Mercury    a0 =  0.3871 ua
*                 Venus      a0 =  0.7233 ua
*                 Earth      a0 =  1.0000 ua
*                 Mars       a0 =  1.5237 ua
*                 Jupiter    a0 =  5.2026 ua
*                 Saturn     a0 =  9.5547 ua
*                 Uranus     a0 = 19.2181 ua
*                 Neptune    a0 = 30.1096 ua
*
*     lu       logical unit index of the file (integer).
*              The file corresponds to a version of VSOP87 theory and
*              a body, and it must be defined and opened before the
*              first call to subroutine VSOP87.
*
*     Output :
*
*     r(6)     array of the results (real double precision).
*
*              for elliptic coordinates :
*                  1: semi-major axis (au)
*                  2: mean longitude (rd)
*                  3: k = e*cos(pi) (rd)
*                  4: h = e*sin(pi) (rd)
*                  5: q = sin(i/2)*cos(omega) (rd)
*                  6: p = sin(i/2)*sin(omega) (rd)
*                     e:     eccentricity
*                     pi:    perihelion longitude
*                     i:     inclination
*                     omega: ascending node longitude
*
*              for rectangular coordinates :
*                  1: position x (au)
*                  2: position y (au)
*                  3: position z (au)
*                  4: velocity x (au/day)
*                  5: velocity y (au/day)
*                  6: velocity z (au/day)
*
*              for spherical coordinates :
*                  1: longitude (rd)
*                  2: latitude (rd)
*                  3: radius (au)
*                  4: longitude velocity (rd/day)
*                  5: latitude velocity (rd/day)
*                  6: radius velocity (au/day)
*
*     ierr     error index (integer).
*                  0: no error.
*                  1: file error (check up ivers index).
*                  2: file error (check up ibody index).
*                  3: precision error (check up prec parameter).
*                  4: reading file error.
*
*-----------------------------------------------------------------------
*
*     --------------------------------
*     Declarations and initializations
*     --------------------------------
*
      implicit double precision (a-h,o-z)
      character*7 bo,body(0:9)
      dimension r(6),t(-1:5),a0(0:9)
      data body/'SUN','MERCURY','VENUS','EARTH','MARS','JUPITER',
     .          'SATURN','URANUS','NEPTUNE','EMB'/
      data a0/0.01d0,0.3871d0,0.7233d0,1.d0,1.5237d0,5.2026d0,
     .        9.5547d0,19.2181d0,30.1096d0,1.d0/
      data dpi/6.283185307179586d0/
      data t/0.d0,1.d0,5*0.d0/
      data t2000/2451545.d0/
      data a1000/365250.d0/
      k=0
      ideb=0
*
      rewind (lu,err=500)
      do i=1,6
         r(i)=0.d0
      enddo
*
      t(1)=(tdj-t2000)/a1000
      do i=2,5
         t(i)=t(1)*t(i-1)
      enddo
*
      ierr=3
      if (prec.lt.0.d0.or.prec.gt.1.d-2) return
      q=dmax1(3.d0,-dlog10(prec+1.d-50))
*
*     -------------------------------------
*     File reading and substitution of time
*     -------------------------------------
*
100   continue
      read (lu,1001,end=400,err=500) iv,bo,ic,it,in
*
      if (ideb.eq.0) then
         ideb=1
         ierr=1
         if (iv.ne.ivers) return
         ierr=2
         if (bo.ne.body(ibody)) return
         ierr=0
         if (iv.eq.0) k=2
         if (iv.eq.2.or.iv.eq.4) k=1
      endif
*
      if (in.eq.0) goto 100
*
      p=prec/10.d0/(q-2)/(dabs(t(it))+it*dabs(t(it-1))*1.d-4+1.d-50)
      if (k.eq.0.or.(k.ne.0.and.ic.eq.5-2*k)) p=p*a0(ibody)
*
      do 200 n=1,in
         nn=n
         read (lu,1002,err=500) a,b,c
         if (dabs(a).lt.p) goto 300
         u=b+c*t(1)
         cu=dcos(u)
         r(ic)=r(ic)+a*cu*t(it)
         if (iv.eq.0) goto 200
         su=dsin(u)
         r(ic+3)=r(ic+3)+t(it-1)*it*a*cu-t(it)*a*c*su
200   continue
*
      goto 100
300   continue
*
      if (nn.eq.in) goto 100
*
      do n=nn+1,in
         read (lu,1002,err=500)
      enddo
      goto 100
*
400   continue
      if (iv.ne.0) then
         do i=4,6
            r(i)=r(i)/a1000
         enddo
      endif
*
      if (k.eq.0) return
*
      r(k)=dmod(r(k),dpi)
      if (r(k).lt.0.d0) r(k)=r(k)+dpi
      return
*
500   continue
      ierr=4
      return
*
*     -------
*     Formats
*     -------
*
1001  format (17x,i1,4x,a7,12x,i1,17x,i1,i7)
1002  format (79x,f18.11,f14.11,f20.11)
*
      end

max_issues_repo_name: benitocm/my-orbits


text:
subroutine radiat(daynr,hour,dayl,sinld,cosld,avrad,sinb,pardir,pardif)
! ----------------------------------------------------------------------
! --- Authors: this routine Astro is based on Sastro (Daniel van Kraalingen)
! --- calculates the fluxes of diffuse and direct photosynthetically
! --- active radiation from the total daily shortwave radiation actually
! --- received (avrad) for a given day of the year and hour of the day.
! --- the input variables dayl, sinld and cosld are calculated in astro.
! --- for more information: see spitters et al. (1988).
! ----------------------------------------------------------------------
      implicit none

      integer daynr

      real  aob,atmtr,avrad,cosld,dayl,dsinb,dsinbe,dso,frdif,hour
      real  par,pardif,pardir,pi,sc,sinb,sinld

      data    pi /3.1415926d0/
! ----------------------------------------------------------------------
! --- calculations on solar elevation
! --- sine of solar elevation sinb
      aob = sinld/cosld
      sinb = max (0.0d0,sinld+cosld*cos(2.0*pi*(hour+12.0d0)/24.0))
! --- integral of sinb
      dsinb = 3600.*(dayl*sinld+24.*cosld*sqrt(1.0d0-aob*aob)/pi)
! --- integral of sinb, corrected for lower atmospheric transmission
! --- at low solar elevations
      dsinbe = 3600.*(dayl*(sinld+0.4*(sinld*sinld+cosld*cosld*0.5))+12.0*cosld*(2.0d0+3.0*0.4*sinld)*sqrt(1.0d0-aob*aob)/pi)

! --- solar constant and daily extraterrestrial radiation
      sc = 1370.*(1.0d0+0.033*cos(2.0*pi*daynr/365.))
      dso = sc*dsinb

! --- diffuse light fraction from atmospheric transmission
      atmtr = avrad/dso
      if (atmtr.gt.0.75d0)  frdif = 0.23d0
      if (atmtr.le.0.75d0.and.atmtr.gt.0.35d0) frdif = 1.33d0-1.46*atmtr
      if (atmtr.le.0.35d0.and.atmtr.gt.0.07d0) frdif = 1.0d0-2.3*(atmtr-0.07d0)**2
      if (atmtr.le.0.07d0) frdif = 1.0d0

! --- photosynthetic active radiation, diffuse and direct
      par = 0.5*avrad*sinb*(1.0d0+0.4*sinb)/dsinbe
      pardif = min (par,sinb*frdif*atmtr*0.5*sc)
      pardir = par-pardif

      return

end subroutine radiat

max_stars_repo_name: vstocca/dssat-csm-os

max_issues_repo_name: vstocca/dssat-csm-os

max_forks_repo_name: vstocca/dssat-csm-os