* PACKAGE XMKSSTP2  !" SST の作成 ( patch ) 矩形状
*
**************************************************************************
      PROGRAM MKSST
*
*   [PARAM] 
#if   SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
      INCLUDE   (ZHDIM)                      !" 文字列文字数
#else
#include        "zcdim.F"                    !" 格子点数, 波数
#include        "zhdim.F"                    !" 文字列文字数
#endif
*
*"  [VAR] 
      REAL        GSST  ( IDIM , JDIM  )
      REAL        GPAT  ( IDIM , JDIM  )
      REAL        ALON  ( IDIM )             !" 経度 単位 rad
      REAL        ALAT  ( IDIM )             !" 緯度 単位 rad
      REAL        SIG   ( KDIM )             !" σ
      REAL        SIGM  ( KDIM+1 )
      REAL        DLON  ( IDIM )             !" 経度加重
      REAL        DLAT  ( IDIM )             !" 緯度加重
      REAL        DSIG  ( KDIM )
      REAL        DSIGM ( KDIM+1 )
*
*   [INTERNAL PARM] 
      CHARACTER  HITEM*(NCC)
      CHARACTER  HTITL*(NCC)
      CHARACTER  HUNIT*(NCC)
      DATA       HITEM / 'GRSST' /
      DATA       HTITL / 'Sea Surface Temp' /
      DATA       HUNIT / 'K' /
      INTEGER    IT
      INTEGER    IDATE(6)
      INTEGER    ISTEP
      INTEGER    ITDUR
      INTEGER    IOMODE
      INTEGER    NOEND
      DATA       IT    / 0 /
      DATA       IDATE / 1992, 3, 21, 0, 0, 0 /
      DATA       ISTEP  / 0 /
      DATA       ITDUR  / 1 /
      DATA       IOMODE / 1 /
      DATA       NOEND  / 0 /
*
      NAMELIST  /NMDATE/ IT, IDATE, ISTEP
*
      CHARACTER  DDSET *(NCC)
      DATA       DDSET / ' ' /
      CHARACTER  HDFMT *(NCC)
      DATA       HDFMT  / 'UR4' /
      REAL       TEQ                 !" 赤道上の SST
      REAL       BLAT0               !" 最高温度緯度 単位度
      REAL       BLAT1               
      REAL       BETA1
      REAL       BETA2
      REAL       BETA3
      REAL       BLACON              !" 大体一定緯度幅 単位度
*
      REAL       CENTERX, CENTERY    !" パッチ中心 単位 度
      REAL       RENGEX,  RENGEY     !" パッチ大きさ 単位 度
      REAL       PATAMP              !" パッチ最大振幅
      INTEGER    IPTYP           !" 0 なし 1 楕円 2 矩形 3 波
      INTEGER    JPTYP
*
      DATA       TEQ   , BLAT0, BLAT1, BETA1,  BETA2, BETA3
     &         / 300.  ,  0.  , 30.  , 60.  , 32.  , 0.     /
*      DATA       TEQ   , BLAT0, BLAT1, BETA1,  BETA2, BETA3
*     &         / 302.  ,  0.  , 30.  , 60.  , 32.  , 0.     /
      DATA       BLACON / 7.   /
*
      NAMELIST  /NMSST/ TEQ, BLAT0, BLAT1, BETA1, BETA2, BETA3, BLACON
*
      DATA  CENTERX, CENTERY, RENGEX, RENGEY, PATAMP, IPTYP,JPTYP 
     &    /   180. ,   0.   ,   30.  ,  10.   ,  2. ,  2  , 1  /
*
      NAMELIST  /PAT/  CENTERX, CENTERY, RENGEX, RENGEY, 
     &                 PATAMP, IPTYP , JPTYP
*
*   [INTERNAL WORK] 
      INTEGER    JFILE
*
      CALL YPREP
      CALL PCONST
      CALL GTOPEN
*
      WRITE ( 6,*      ) ' DATASET NAME? '
      READ  ( 5,'(A)'  ) DDSET
      CALL GFWOPQ ( JFILE, 'SST OUTPUT' )
*
      WRITE ( 6,NMSST  )
      WRITE ( 6,*      ) ' ENTER NAMELST NMSST?'
      READ  ( 5,NMSST,  END=1900  )
 1900 CONTINUE
*
      WRITE ( 6,NMDATE )
      WRITE ( 6,*      ) ' ENTER NAMELST NMDATE?'
      READ  ( 5,NMDATE, END=1910 )
 1910 CONTINUE
*
      WRITE ( 6,PAT )
      WRITE ( 6,*      ) ' ENTER NAMELST PAT?'
      READ  ( 5,PAT, END=1920 )
 1920 CONTINUE
*
      CALL GTCSET( 'DDSET' , DDSET    )
*
      CALL SETCOR
     O         ( ALON  , DLON  ,
     O           ALAT  , DLAT  ,
     O           SIG   , DSIG  ,
     O           SIGM  , DSIGM  )
*
*" 東西一様 SST
*
      CALL CALSST
     O         ( GSST  ,
     I           ALAT  ,
     I           TEQ   , BLAT0 , BLAT1 , BETA1 , BETA2  , BETA3 ,
     I           BLACON,
     I           IDIM  , JDIM                                   )
*
*" パッチ
*      
      IF ( IPTYP .EQ. 1 ) THEN
        CALL SSTP1
     M         ( GSST  ,
     I           ALAT  , 
     I           CENTERX, CENTERY, RENGEX, RENGEY, PATAMP, JPTYP, 
     I           IDIM  , JDIM , 
     W           GPAT            )
*
      ELSEIF ( IPTYP .EQ. 2 ) THEN
        CALL SSTP2
     M         ( GSST  ,
     I           ALAT  , 
     I           CENTERX, CENTERY, RENGEX, RENGEY, PATAMP, JPTYP, 
     I           IDIM  , JDIM , 
     W           GPAT            )
*
      ELSEIF ( IPTYP .EQ. 3 ) THEN
        CALL SSTP3
     M         ( GSST  ,
     I           ALAT  , 
     I           CENTERX, CENTERY, RENGEX, RENGEY, PATAMP, JPTYP, 
     I           IDIM  , JDIM , 
     W           GPAT            )
*
      ENDIF
*
      WRITE ( 6,* ) (GSST(1,J),J=1,JMAX)
*
      CALL GDWRIT
     O         ( GSST  ,
     I           HITEM , HTITL , HUNIT ,
     I           IT    , IDATE , ISTEP , ITDUR ,
     I           JFILE , IOMODE, NOEND , 'ASFC', HDFMT  )
*
      CALL YFINE
*
      STOP
      END
***********************************************************************
      SUBROUTINE CALSST        !" SST 計算
     O         ( GSST  ,
     I           ALAT  ,
     I           TEQ   , BLAT0 , BLAT1 , BETA1 , BETA2  , BETA3 ,
     I           BLACON,
     I           IDIM  , JDIM                                   )
*
*     SST = SST ( TEQ-BLAT0, BLAT1, BETA1, BETA2, BETA3 )
*         -> BLAT0 +/- BLACON の間, 平坦化
*
*   [PARAM] 
      INTEGER    IDIM, JDIM
*
*   [OUTPUT] 
      REAL       GSST  ( IDIM, JDIM )
*
*   [INPUT] 
      REAL       ALAT  ( JDIM )
      REAL       TEQ
      REAL       BLAT0
      REAL       BLAT1
      REAL       BETA1
      REAL       BETA2
      REAL       BETA3
      REAL       BLACON
*
*   [INTERNAL WORK]       
      REAL       PI, PHI1, AB4, PHI, GSSTP, ALATP, ALATM, GSSTMX
      INTEGER    I, J, JP, JM
*
*" << CRSST : SSTの計算作成 >> 
*
      PI   = ATAN( 1. ) * 4.  
      PHI1 = ABS( BLAT1 * PI/180.   )
      AB4  = 2. *( PHI1**3 )*BETA2/BETA1
*
      DO 2100 J = 1, JDIM
*
         PHI   = ABS( ALAT( J ) - BLAT0 *PI/180. )
         GSSTP = TEQ
     &         - BETA1/2.  
     &           * ( PHI - MAX(   SQRT(  PHI1**2     +AB4 )
     &                           -SQRT( (PHI-PHI1)**2+AB4 ), 0.   )  )
     &         + BETA3 *( PHI**3 )
*
         DO 2110 I = 1, IDIM
            GSST( I,J ) = GSSTP
 2110    CONTINUE
*
 2100 CONTINUE
*
*" << FIXSST : 中心+/-BLACONを平坦に >> 
*
      ALATP = ( BLAT0 + BLACON )*PI/180.  
      ALATM = ( BLAT0 - BLACON )*PI/180.  
*
      JP = 1
      JM = JDIM
      DO 2200 J = 1, JDIM
         IF ( ALAT( J ) .GT. ALATP ) THEN
            JP = J
         ENDIF
         IF ( ALAT( J ) .GE. ALATM ) THEN
            JM = J
         ENDIF
 2200 CONTINUE
*
      GSSTMX = (     GSST(1,JP) * ( ALATP - ALAT(JP+1))
     &            +  GSST(1,JP+1) * ( ALAT(JP) - ALATP ))
     &        /( ALAT(JP) - ALAT(JP+1) )
*
      DO 2500 J = JP+1, JM
         DO 2510 I = 1, IDIM
            GSST( I,J ) = GSSTMX
 2510    CONTINUE
 2500 CONTINUE
*
*
      RETURN
      END
***********************************************
      SUBROUTINE  SSTP1             !" 楕円状パッチ
     M         ( GSST  ,
     I           ALAT  , 
     I           CENTERX, CENTERY, RENGEX, RENGEY, PATAMP, JPTYP, 
     I           IDIM  , JDIM , 
     W           GPAT            )

*   [PARAM] 
      INTEGER    IDIM, JDIM
*
*   [OUTPUT] 
      REAL       GSST  ( IDIM, JDIM )
*
*   [INPUT] 
      REAL       ALAT  ( JDIM )
      REAL       CENTERX, CENTERY  
      REAL       RENGEX , RENGEY 	
      REAL       PATAMP
      INTEGER    JPTYP  
*
*   [INTERNAL WORK]       
      REAL       PI
      REAL       ACX , ACY 
      REAL       ARX , ARY
      INTEGER    I, J
*
      REAL       GPAT ( IDIM, JDIM ) 
      REAL       ALON, RAD
*
*" 定数設定
*
      PI   = ATAN( 1. ) * 4.  
*
      ACX   = CENTERX * PI / 180.
      ACY   = CENTERY * PI / 180.
      ARX   = RENGEX  * PI / 180.
      ARY   = RENGEY  * PI / 180.
*
      CALL RESET ( GPAT , IDIM * JDIM )     
*
      DO 1000 I = 1, IDIM
        ALON = 1. / IDIM * ( I - 1 ) * 2. * PI 
        DO 1000 J = 1, JDIM
*
          RAD =   ( ( ALON      - ACX ) / ARX )**2 
     &          + ( ( ALAT( J ) - ACY ) / ARY )**2 


        IF ( JPTYP .EQ. 1 ) THEN
           IF ( RAD .LE. 1. ) THEN
              GPAT ( I , J ) = PATAMP
           ENDIF
        ELSEIF ( JPTYP .EQ. 2 ) THEN
           IF ( RAD .LE. 1. ) THEN
              GPAT ( I , J ) = PATAMP * ( 1. -SQRT ( RAD ))
           ENDIF
        ELSEIF ( JPTYP .EQ. 3 ) THEN
           IF ( RAD .LE. 1. ) THEN
              GPAT ( I , J ) = PATAMP * ( 1. - RAD )
           ENDIF
        ELSE 
           WRITE ( 6 , * ) 'NO PATCH'
        ENDIF
*
        GSST ( I , J ) = GSST ( I , J ) + GPAT ( I , J )
*      
 1000 CONTINUE
*
      RETURN
      END
***********************************************
      SUBROUTINE  SSTP2             !" 矩形パッチ
     M         ( GSST  ,
     I           ALAT  , 
     I           CENTERX, CENTERY, RENGEX, RENGEY, PATAMP, JPTYP, 
     I           IDIM  , JDIM , 
     W           GPAT            )

*   [PARAM] 
      INTEGER    IDIM, JDIM
*
*   [OUTPUT] 
      REAL       GSST  ( IDIM, JDIM )
*
*   [INPUT] 
      REAL       ALAT  ( JDIM )
      REAL       CENTERX, CENTERY  
      REAL       RENGEX , RENGEY 	
      REAL       PATAMP
      INTEGER    JPTYP   !" 1 一様 2 錐 3 x 方向斜面 4 y方向斜面
*
*   [INTERNAL WORK]       
      REAL       PI
      REAL       ACX , ACY 
      REAL       ARX , ARY
      INTEGER    I, J
*
      REAL       GPAT ( IDIM, JDIM ) 
      REAL       ALON, RAD
*
*" 定数設定
*
      PI   = ATAN( 1. ) * 4.  
*
      ACX   = CENTERX * PI / 180.
      ACY   = CENTERY * PI / 180.
      ARX   = RENGEX  * PI / 180.
      ARY   = RENGEY  * PI / 180.
*
      CALL RESET ( GPAT , IDIM * JDIM )     
*
      DO 1000 I = 1, IDIM
        ALON = 1. / IDIM * ( I - 1 ) * 2. * PI 
        DO 1000 J = 1, JDIM
*
        IF ( JPTYP .EQ. 1 ) THEN
           IF ( ABS ( ALAT(J)- ACY ) .LE. ARY 
     &          .AND. ABS ( ALON - ACX ) .LE. ARX ) THEN


              GPAT ( I , J ) = PATAMP
           ENDIF
        ELSEIF ( JPTYP .EQ. 2 ) THEN
           IF ( ABS ( ALAT(J) - ACY ) .LE. ARY 
     &          .AND. ABS ( ALON - ACX ) .LE. ARX ) THEN
              RAD = MAX ( ABS ( ALON - ACX ) / ARX ,
     &                    ABS ( ALAT(J) - ACY ) / ARY )
              GPAT ( I , J ) = PATAMP * ( 1. - RAD ) 
           ENDIF
        ELSEIF ( JPTYP .EQ. 3 ) THEN
           IF ( ABS ( ALAT(J) - ACY ) .LE. RENGEY 
     &          .AND. ABS ( ALON - ACX ) .LE. ARX ) THEN
              RAD =  ABS ( ALON - ACX ) / ARX 
              GPAT ( I , J ) = PATAMP * ( 1. - RAD )
           ENDIF
        ELSEIF ( JPTYP .EQ. 4 ) THEN
           IF ( ABS ( ALAT(J) - ACY ) .LE. ARY 
     &          .AND. ABS ( ALON - ACX ) .LE. ARX ) THEN
              RAD = ABS ( ALAT(J) - ACY ) / ARY 
              GPAT ( I , J ) = PATAMP * ( 1. - RAD ) 
           ENDIF
        ELSE 
           WRITE ( 6 , * ) 'NO PATCH'
        ENDIF
*
        GSST ( I , J ) = GSST ( I , J ) + GPAT ( I , J )
*      
 1000 CONTINUE
*
      RETURN
      END
***********************************************
      SUBROUTINE  SSTP3             !" 波上ノイズ
     M         ( GSST  ,
     I           ALAT  , 
     I           CENTERX, CENTERY, RENGEX, RENGEY, PATAMP, JPTYP, 
     I           IDIM  , JDIM , 
     W           GPAT            )

*   [PARAM] 
      INTEGER    IDIM, JDIM
*
*   [OUTPUT] 
      REAL       GSST  ( IDIM, JDIM )
*
*   [INPUT] 
      REAL       ALAT  ( JDIM )
      REAL       CENTERX, CENTERY  
      REAL       RENGEX , RENGEY 	
      REAL       PATAMP
      INTEGER    JPTYP   !" 1 一様 2 y 方向に直線減衰
*
*   [INTERNAL WORK]       
      REAL       PI
      REAL       ACX , ACY 
      REAL       ARX , ARY    !" ARX は周期
      INTEGER    I, J
*
      REAL       GPAT ( IDIM, JDIM ) 
      REAL       ALON, RAD
*
*" 定数設定
*
      PI   = ATAN( 1. ) * 4.  
*
      ACX   = CENTERX * PI / 180.
      ACY   = CENTERY * PI / 180.
      ARX   = RENGEX  * PI / 180.
      ARY   = RENGEY  * PI / 180.
*
      CALL RESET ( GPAT , IDIM * JDIM )     
*
      DO 1000 I = 1, IDIM
        ALON = 1. / IDIM * ( I - 1 ) * 2. * PI 
        DO 1000 J = 1, JDIM
*
        RAD = ABS ( ALAT (J) - ACY ) / ARY
*
        IF ( JPTYP .EQ. 1 ) THEN
           IF ( RAD .LE. 1 ) then
              GPAT ( I , J ) 
     &          = PATAMP * COS ( 2. * PI *( ALON - ACX ) / ARX )
           ENDIF
        ELSEIF ( JPTYP .EQ. 2 ) THEN
           IF ( RAD .LE. 1  ) then
              GPAT ( I , J ) 
     &          = PATAMP * COS ( 2. * PI *( ALON - ACX ) / ARX )
     &                   * ( 1. - RAD )
           ENDIF
        ELSE 
           WRITE ( 6 , * ) 'NO PATCH'
        ENDIF
*
        GSST ( I , J ) = GSST ( I , J ) + GPAT ( I , J )
*      
 1000 CONTINUE
*
      RETURN
      END
