!----------------------------------------------------------------------
!  Copyright (C) 2001--2008 SPMODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
! Sample program for gtool_history/gtool5 and ISPACK   2002/08/21 S.Takehiro
!                                                      2004/01/26 M.Odaka
!                                                      2009/02/27 Y.Morikawa
!
! Solving a linear 2-D shallow water system on a sphere
!     with an isolated mountain
!     du/dt + u0/a\cos\phi du/d\lambda + v/a du0/d\phi - u0v\tan\phi/a
!           - 2\Omega\sin\phi v = -g/a\cos\phi dh/d\lambda - \nu\lapla^4 u,
!     dv/dt + u0/a\cos\phi dv/d\lambda + 2 u0 u \tan\phi/a
!           + 2\Omega\sin\phi u  =  -g/a  dh/d\phi - \nu\lapla^4 v,
!     dh/dt + ( 1/\cos\phi d( (H+h0)u+(h-ht)u0 )/d\lambda
!               + 1/\cos\phi d( (H+h0)v\cos\phi)/d\phi ) = - \nu\lapla^4 h.
!
! A setup is similar to the experiment of Grose and Hoskins (1979)
! with a superrotating rigid-rotation zonal wind profile.
!
  program sp_topo_gtauto

    use w_module
    use gtool5
    use gtool_historyauto
    implicit none

  !---- ç©ºé–“è§£åƒåº¦è¨å®š ----
    integer, parameter :: im=32, jm=16             ! æ ¼åç‚¹ã®è¨å®š(X,Y)
    integer, parameter :: nm=10

  !---- å¤‰æ•° ----
    real(8)            :: xy_U(0:im-1,jm)          ! æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿(é€Ÿåº¦çµŒåº¦æˆåˆ†)
    real(8)            :: xy_V(0:im-1,jm)          ! æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿(é€Ÿåº¦ç·¯åº¦æˆåˆ†)
    real(8)            :: xy_H(0:im-1,jm)          ! æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿(å¤‰ä½)

    real(8)            :: w_U((nm+1)*(nm+1))       ! ã‚¹ãƒšã‚¯ãƒˆãƒ«(é€Ÿåº¦çµŒåº¦æˆåˆ†)
    real(8)            :: w_V((nm+1)*(nm+1))       ! ã‚¹ãƒšã‚¯ãƒˆãƒ«(é€Ÿåº¦ç·¯åº¦æˆåˆ†)
    real(8)            :: w_H((nm+1)*(nm+1))       ! ã‚¹ãƒšã‚¯ãƒˆãƒ«(å¤‰ä½)

    real(8)            :: xy_U0(0:im-1,jm)         ! æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿(é€Ÿåº¦çµŒåº¦æˆåˆ†)
    real(8)            :: xy_H0(0:im-1,jm)         ! æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿(å¤‰ä½)
    real(8)            :: w_U0((nm+1)*(nm+1))      ! ã‚¹ãƒšã‚¯ãƒˆãƒ«(é€Ÿåº¦çµŒåº¦æˆåˆ†)

    real(8)            :: xy_Zeta(0:im-1,jm)       ! æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿(æ¸¦åº¦)
    real(8)            :: xy_Htopo(0:im-1,jm)      ! æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿(åœ°å½¢)

  !---- æ™‚é–“ç©åˆ†ãƒ‘ãƒ©ãƒ¡ã‚¿ãƒ¼ ----
    real(8):: dt                                   ! æ™‚é–“é–“éš” [s]
    type(DC_DIFFTIME):: deltime, dispint           ! æ™‚é–“é–“éš”, å‡ºåŠ›æ™‚é–“é–“éš”
    type(DC_DIFFTIME):: endtime                    ! è¨ˆç®—çµ‚äº†æ™‚é–“
    type(DC_DIFFTIME):: ct                         ! ç¾åœ¨æ™‚åˆ»

  !---- ç‰©ç†ãƒ‘ãƒ©ãƒ¡ã‚¿ãƒ¼ ----
    real(8), parameter :: pi = 3.141592653589793   ! å††å‘¨çŽ‡
    real(8), parameter :: R0 = 6.371D6             ! çƒã®åŠå¾„
    real(8), parameter :: H0 = 10D3                ! å¹³å‡å¤‰ä½
    real(8), parameter :: Omega = 7.292D-5         ! å›žè»¢è§’é€Ÿåº¦
    real(8), parameter :: Grav = 9.8               ! é‡åŠ›åŠ é€Ÿåº¦

    real(8), parameter :: U0   = 15                ! å¹³å‡æµ
    real(8), parameter :: Lon0 = pi, Lat0=pi/6     ! å±±ã®ä¸å¿ƒä½ç½®
    real(8), parameter :: W0 = 8                   ! å±±ã®åºƒãŒã‚Šã®é€†æ•°
    real(8), parameter :: Htopo0= 2500             ! å±±ã®é«˜ã•

    integer, parameter :: ndiff=4                  ! è¶…ç²˜æ€§ã®æ¬¡æ•°
    real(8), parameter :: Nu=2.338D16              ! è¶…ç²˜æ€§ä¿‚æ•°

  !------ NAMELIST ãƒ•ã‚¡ã‚¤ãƒ«å ------
    character(256):: nmlfile

  !------ ä½œæ¥å¤‰æ•° ------
    integer            :: i, j
    real(8)            :: Delta

  !---- æ™‚é–“ç©åˆ†ãƒ‘ãƒ©ãƒ¡ã‚¿ãƒ¼ã®è¨å®š ----
    dt = 100.                                      ! æ™‚é–“é–“éš” [s]
    call DCDiffTimeCreate(deltime, dt, 'sec')      ! æ™‚é–“é–“éš”
    dispint = deltime * 500                        ! å‡ºåŠ›æ™‚é–“é–“éš”
    endtime = deltime * 10000                      ! è¨ˆç®—çµ‚äº†æ™‚é–“
    call DCDiffTimeCreate(ct,      0., 'sec')      ! ç¾åœ¨æ™‚åˆ»

  !---------------- åº§æ¨™å€¤ã®è¨å®š ---------------------
    call w_Initial(nm,im,jm)                ! ISPACKåˆæœŸåŒ–

  !------------------- åˆæœŸå€¤è¨å®š ----------------------
    xy_U0  = U0*cos(xy_Lat)
    xy_H0  = ( Omega*R0*U0/(2*Grav) + U0**2/(4*Grav) )*cos(2*xy_Lat)

    xy_U  = 0 ; xy_V  = 0 ; xy_H  = 0

    w_U0 = w_xy(xy_U0) !; w_H0 = w_xy(xy_H0)
    w_U = w_xy(xy_U) ; w_V = w_xy(xy_V) ; w_H = w_xy(xy_H)

  !------------------- åœ°å½¢è¨å®š ----------------------
    do j=1,jm
       do i=0,im-1
          delta = acos( cos(xy_Lat(i,j))*cos(Lat0)*cos(xy_Lon(i,j)-Lon0) &
                       + sin(xy_Lat(i,j))*sin(Lat0) )
          if ( abs(delta) <= pi/W0 ) then
             xy_Htopo(i,j)= htopo0 * (1+cos(W0*Delta))/2
          else
             xy_Htopo(i,j)= 0
          endif
       enddo
    enddo

  !------------------- ãƒ’ã‚¹ãƒˆãƒªãƒ¼åˆæœŸè¨å®š ----------------------
    call output_gtool5_init                        ! ãƒ’ã‚¹ãƒˆãƒªãƒ¼åˆæœŸåŒ–
    call output_gtool5                             ! åˆæœŸå€¤å‡ºåŠ›

  !------------------- æ™‚é–“ç©åˆ† ----------------------
    do while ( ct <= endtime )
       if ( mod(ct, dispint) == 0 ) then
         write(6,*) 'it = ', int( ct / deltime )
       end if

       ct = ct + deltime  ! æ™‚åˆ»ã®é€²è¡Œ

       w_U = ( w_U &
               + dt * w_xy( - xy_U0 * xy_GradLon_w(w_U) / R0   &
                            - xy_V  * xy_GradLat_w(w_U0) / R0  &
                            + xy_U0 * xy_V * tan(xy_Lat) / R0  &
                            + 2 * Omega * sin(xy_Lat) * xy_V   &
                           - Grav * xy_GradLon_w(w_H)/ R0   ) &
             )/(1+Nu*(-rn(:,1)/R0**2)**(ndiff/2)*dt)

       xy_U = xy_w(w_U)

       w_V = ( w_V &
              + dt * w_xy( - xy_U0 * xy_GradLon_w(w_V) / R0   &
                           - xy_U * xy_U0 * tan(xy_Lat) / R0  &
                           - xy_U0 * xy_U * tan(xy_Lat) / R0  &
                           - 2 * Omega * sin(xy_Lat) * xy_U   &
                           - Grav * xy_GradLat_w(w_H) / R0 )  &
             )/(1+Nu*(-rn(:,1)/r0**2)**(ndiff/2)*dt)

       xy_V = xy_w(w_V)

       w_H = ( w_H &
              + dt * ( - w_Div_xy_xy( xy_U*(H0+xy_H0)+xy_U0*(xy_H-xy_Htopo),   &
                                      xy_V*(H0+xy_H0) ) / R0      ) &
              )/(1+Nu*(-rn(:,1)/R0**2)**(ndiff/2)*dt)

       xy_H = xy_w(w_H)

       call output_gtool5  ! å‡ºåŠ›
    enddo

    call output_gtool5_close  ! ãƒ•ã‚¡ã‚¤ãƒ«ã®ã‚¯ãƒãƒ¼ã‚º
    stop

  contains
    subroutine output_gtool5_init
      write(6,'(a)',advance='NO') '  Input NAMELIST file: '
      read (5,'(a)') nmlfile

      call HistoryAutoCreate( &                               ! ãƒ’ã‚¹ãƒˆãƒªãƒ¼ä½œæˆ
           title='Shallow water equation on a sphere',             &
           source='Sample program of gtool_historyauto/gtool5',    &
           institution='GFD_Dennou Club davis/spmodel project',    &
           dims=(/'lon','lat','t  '/), dimsizes=(/im,jm,0/),       &
           longnames=(/'longitude','latitude ','time     '/),      &
           units=(/'degree_east ','degree_north','sec.        '/), &
           origin=ct, interval=dispint, terminus=endtime,          &
           namelist_filename=nmlfile )

      call HistoryAutoPutAxis('lon',x_Lon*180/pi)              ! å¤‰æ•°å‡ºåŠ›
      call HistoryAutoAddAttr('lon','topology','circular')     ! å‘¨æœŸå±žæ€§
      call HistoryAutoAddAttr('lon','modulo',360.0)            ! å‘¨æœŸå±žæ€§
      call HistoryAutoPutAxis('lat',y_Lat*180/pi)              ! å¤‰æ•°å‡ºåŠ›

      call HistoryAutoAddVariable( &                          ! å¤‰æ•°å®šç¾©
           varname='h', dims=(/'lon','lat','t  '/), &
           longname='surface displacement ', units='m')

      call HistoryAutoAddVariable( &                          ! å¤‰æ•°å®šç¾©
           varname='u', dims=(/'lon','lat','t  '/), &
           longname='velocity(longitude) ', units='m/s')

      call HistoryAutoAddVariable( &                          ! å¤‰æ•°å®šç¾©
           varname='v', dims=(/'lon','lat','t  '/), &
           longname='velocity(latitude) ', units='m/s')

      call HistoryAutoAddVariable( &                          ! å¤‰æ•°å®šç¾©
           varname='zeta', dims=(/'lon','lat','t  '/), &
           longname='vorticity', units='1/s')

    end subroutine output_gtool5_init

    subroutine output_gtool5
      call HistoryAutoPut(ct, 'u', xy_U)
      call HistoryAutoPut(ct, 'v', xy_V)
      call HistoryAutoPut(ct, 'h', xy_H)
      xy_Zeta = xy_w(w_Divlon_xy(xy_V) - w_Divlat_xy(xy_U))/r0
      call HistoryAutoPut(ct, 'zeta', xy_Zeta)
    end subroutine output_gtool5

    subroutine output_gtool5_close
      call HistoryAutoClose
    end subroutine output_gtool5_close

  end program sp_topo_gtauto