* PACKAGE P2RADI   
*
*"履歴 90/10/09 沼口  敦
*"     92/04/17 沼口  敦  
*"     92/08/27 竹広真一
*
***********************************************************************
      SUBROUTINE SHTINS            
     O         ( SINS  , SECZET ,
     I           IDATE ,
     C           ALON  , ALAT    )
*
*"   理想的季節・日変化(オプション)
*
*   [PARAM] 




* PACKAGE  ZXT10L32 !" 格子点数，波数 X-Z
*
      INTEGER     IMAX           !" 経度格子数
      INTEGER     JMAX           !" 緯度格子数
      INTEGER     IDIM           !" 経度配列サイズ
      INTEGER     JDIM           !" 緯度配列サイズ
      INTEGER     MMAX           !" 最大東西波数
      INTEGER     LMAX           !" 最大南北波数
      INTEGER     NMAX           !" 最大全波数
      INTEGER     MINT           !" 東西波数間隔
      INTEGER     KMAX           !" 鉛直層数
      INTEGER     KDIM           !" 鉛直配列サイズ
      INTEGER     NMDIM          !" 水平波数成分数
      INTEGER     JMXHF          !" JMAX/2+1
      INTEGER     NWORK          !" ワーク領域大きさ
*
      PARAMETER ( IMAX  = 32,     JMAX =  1,
     &            IDIM  = IMAX+1, JDIM = JMAX,
     &            MMAX  = 10,     LMAX =  0,    NMAX = MMAX, MINT = 1,
     &            KMAX  = 32,     KDIM = KMAX,
     &            NMDIM = (MMAX/MINT+1)*(2*(NMAX+1)-MMAX)
     &                  - (NMAX-LMAX)/MINT*(NMAX-LMAX+1),
     &            JMXHF = JMAX/2+1,
     &            NWORK = IDIM*JDIM*KMAX*12                       )

*
*   [OUTPUT] 
      ReaL*8       SINS  ( IDIM, JDIM )        
      ReaL*8       SECZET( IDIM, JDIM )        
*
*   [INPUT] 
      INTEGER    IDATE ( 6 )                 
      ReaL*8       ALON  ( IDIM )              
      ReaL*8       ALAT  ( JDIM )              
*
*   [INTERNAL PARM] 
      ReaL*8       SOLCON
      LOGICAL    OYRAVR, ODYAVR
      ReaL*8       ATALB
      ReaL*8       EPSORB, EQNORB
      ReaL*8       AINS, BINS
      ReaL*8       AZET, BZET
      DATA       SOLCON     / 1380.   /      
      DATA       OYRAVR     / .TRUE. /       
      DATA       ODYAVR     / .TRUE. /       
      DATA       ATALB      / 0.2     /      
      DATA       EPSORB     / 23.5    /      
      DATA       EQNORB     / 110.    /      
      DATA       AINS       / 0.127   /      
      DATA       BINS       / 0.183   /      
      DATA       AZET       / 0.410   /      
      DATA       BZET       / 0.590   /      
*
      NAMELIST  /NMRADI/ SOLCON, OYRAVR, ODYAVR, ATALB ,
     &                   EPSORB, EQNORB , 
     &                   AINS  , BINS  , AZET  , BZET
      SAVE
*
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
      SAVE       OFIRST
*
*   [INTERNAL WORK] 
      ReaL*8       RDATE ( 2 )                 
*
      INTEGER    I, J
      INTEGER    IFPAR, JFPAR
      ReaL*8       PI, COSZET, ANGHR, SINDEL
*      
*   [ONCE] 
      IF ( OFIRST ) THEN
         WRITE ( 6,* ) ' SOLAR INSOLATION  P2.2 DATE=92/04/17'
         OFIRST = .FALSE.
*
         CALL   REWNML ( IFPAR , JFPAR )
         READ   ( IFPAR, NMRADI, END=190 )
  190    WRITE  ( JFPAR, NMRADI )
      ENDIF
*
*"          < 0. 季節, 時刻を小数に >
*
      CALL RDATES
     O      ( RDATE ,
     I        IDATE  )
*
      PI = ATAN( 1. ) *4.
*
      IF ( OYRAVR .AND. ODYAVR ) THEN
*
*"           < 1. 年・日平均日射 >
*
         DO 1100 I = 1, IDIM
            DO 1110 J = 1, JDIM
               SINS( I,J ) = - SOLCON * ( 1. - ATALB )
     &                         * ( AINS + BINS*COS( ALAT(J) )**2 )
*
               IF ( SINS( I,J ) .LT. 0. ) THEN
                  SECZET( I,J ) = 1.  /( AZET + BZET*COS( ALAT(J) )**2 )
               ELSE
                  SINS  ( I,J ) = 0.  
                  SECZET( I,J ) = 0.  
               ENDIF
 1110       CONTINUE
 1100    CONTINUE
*
      ELSE IF ( (.NOT.OYRAVR).AND.(.NOT.ODYAVR) ) THEN
*
*"          < 2. 季節・日変化あり >  
*"                   92/08/27 bug-fixed by takepiro <-- ここから

         SINDEL =  SIN( EPSORB/180.* PI  )
     &           * SIN( ( RDATE( 1 ) - EQNORB/360. )*2.*PI ) 
*
         DO 2100 I = 1, IDIM
            DO 2110 J = 1, JDIM
               ANGHR  = RDATE( 2 )*2.*PI - PI + ALON( I ) 
               COSZET = SIN( ALAT(J) )*SINDEL
     &                + COS( ALAT(J) )*SQRT(1.-SINDEL**2)*COS( ANGHR )
*
               IF ( COSZET .GT. 0. ) THEN
                  SINS  ( I,J ) = - SOLCON * ( 1. - ATALB ) * COSZET
                  SECZET( I,J ) = 1. / COSZET
               ELSE
                  SINS  ( I,J ) = 0.  
                  SECZET( I,J ) = 0.  
               ENDIF
 2110       CONTINUE
 2100    CONTINUE
*
*"                                                  <--- ここまで
*
      ELSE
*
         CALL MSGDMP( 'E', 'SHTINS',
     &                'THIS TYPE OF INSOLATION IS NOT IMPREMENTED' ) 
*
      ENDIF
*
      RETURN
      END


