!<CONTACT EMAIL="petri@pik-potsdam.de"> Stefan Petri
!</CONTACT>
!<REVIEWER EMAIL="">
!
!</REVIEWER>
!-----------------------------------------------------------------------
!<OVERVIEW>
!  Standalone main program for testing the interface between FMS/F90
!  and AtmosphereIII/C++
!</OVERVIEW>

!<DESCRIPTION>

!</DESCRIPTION>

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                   !!
!!                   GNU General Public License                      !!
!!                                                                   !!
!! This file is part of the Flexible Modeling System (FMS).          !!
!!                                                                   !!
!! FMS is free software; you can redistribute it and/or modify       !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation.             !!
!!                                                                   !!
!! FMS is distributed in the hope that it will be useful,            !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of    !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the     !!
!! GNU General Public License for more details.                      !!
!!                                                                   !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to:                                 !!
!!          Free Software Foundation, Inc.                           !!
!!          59 Temple Place, Suite 330                               !!
!!          Boston, MA  02111-1307  USA                              !!
!! or see:                                                           !!
!!          http://www.gnu.org/licenses/gpl.txt                      !!
!!                                                                   !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program testaeolus

  use sinsol_mod, only: sinsol_mod_init

  !-----------------------------------------------------------------------

  implicit none

  !-----------------------------------------------------------------------

  !<PUBLICTYPE>
  type surf_diff_type
     real, pointer, dimension(:,:) :: dtmass  => NULL() ! dt/mass, where dt = atmospheric time step (sec)
                                                        ! mass = mass of lowest atmospheric layer (Kg/m2)
     real, pointer, dimension(:,:) :: dflux_t => NULL() ! derivative of the temperature flux at the top of the lowest
                                                        ! atmospheric layer with respect to the temperature
                                                        ! of that layer  (J/(m2 K))
     real, pointer, dimension(:,:,:) :: dflux_tr => NULL() ! derivative of the flux of specific humidity
                                                        ! at the top of the lowest atmospheric layer with respect to
                                                        ! the specific humidity of that layer  (--/(m2 K))
     real, pointer, dimension(:,:) :: delta_t => NULL() ! the increment in temperature in the lowest atmospheric
                                                        ! layer (((i+1)-(i-1) if atmos model is leapfrog) (K)
                                                        ! (defined in  gcm_vert_diff_down as the increment computed up
                                                        ! to this point in model, including effect of vertical
                                                        ! diffusive flux at top of lowest model level, presumed
                                                        ! to be modified to include effects of surface fluxes
                                                        ! outside of this module, then used to start upward
                                                        ! tridiagonal sweep,
     real, pointer, dimension(:,:,:):: delta_tr => NULL() ! similarly for the increment in specific humidity
                                                        ! (non-dimensional  = Kg/Kg)
     real, pointer, dimension(:,:) :: delta_u => NULL()
     real, pointer, dimension(:,:) :: delta_v => NULL()
     real, pointer, dimension(:,:) :: sst_miz => NULL()
  end type surf_diff_type
  !</PUBLICTYPE>

  !-----------------------------------------------------------------------

  character(len=128) :: version = '$Id$'
  character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $'

  !---- namelist (saved in file input.nml) ----
  !<NAMELIST NAME="atmosphere_nml">
  ! <DATA NAME="NoNx" TYPE="integer">
  !  Longitudinal resolution of atmosphere grid (at the equator).
  !  Must match the data in the grid specification file.
  ! </DATA>
  ! <DATA NAME="NoNy" TYPE="integer">
  !  Latitudinal resolution of atmosphere grid.
  !  Must match the data in the grid specification file.
  ! </DATA>
  ! <DATA NAME="NoNz" TYPE="integer">
  !  Number of vertical layers of atmosphere grid.
  ! </DATA>
  ! <DATA NAME="polar_merging" TYPE="logical">
  !  If true, longitudinal adjacent grid cells are merged towards the poles,
  !  so that cell sizes, and with them the time steps, dont get too small.
  !  If false, simply use a regular lat/lon grid.
  ! </DATA>
  !</NAMELIST>
  namelist /atmosphere_nml/ &
       NoNx, NoNy, NoNz, NoNz_Strato, & ! Number of grid cells in lon, lat, height (Tropo+Strato)
       polar_merging,       & ! merge cells longitudinal towards pole, or regular lat/lon grid?
       lon_0_360              ! whether longitudes run 0 to 360 in the grid spec file

  !-----------------------------------------------------------------------
  !---- private data ----

  integer :: NoNx = 96, NoNy = 48, NoNz = 5, NoNz_Strato = 0  ! values from namelist
  logical :: polar_merging = .true. , lon_0_360 = .true.

  integer :: nr_tracers = 1  ! total number of tracers in surf_diff_type fields
  integer :: q_ind = 1       ! index of humidity in tracers array - for now hardcoded!

  integer :: is, ie, js, je  ! local domain boundary indices

  !------------------------- axis names ----------------------------------
  character(len=8) :: axiset   = 'aeolus'
  character(len=8) :: mod_name = 'aeolus'

  !-----------------------------------------------------------------------

  type(surf_diff_type) :: Surf_diff

  integer :: Time_init_days, Time_init_seconds
  integer :: Time_days, Time_seconds
  integer :: Time_step_days, Time_step_seconds
  integer :: unit, ierr, io
  integer :: log_unit

  integer :: pe, root_pe, npes, layout(2)

  real(kind=8), dimension(:), allocatable :: lonb, latb, lon, lat ! for the axes
  integer :: id_lonb, id_latb, id_lon, id_lat
  integer, dimension(1:2) :: axes_2d

  character(len=64) :: filename

  ! some magical values for testing the interface between Fortran90 and C++ code.
  ! we pass them from Fortran90 to C++ functions and test in there if they were
  ! received correctly.
  integer :: deadbeef, fse;
  real(8) :: onetwothree, pitest
  logical :: true = .true., false = .false.
  integer :: kind_i = KIND(deadbeef), kind_r = KIND(pitest), kind_l = KIND(true)
  real(8), dimension(:,:,:), allocatable :: testarray3d

  real(8), dimension(:,:), allocatable :: &
       solarm, & ! daily solar radiation at the top of the atmosphere
       coszm     ! daily averaged solar zenit angle
  integer :: DaysPerYear

  integer :: i, j, k
  
  ! start of program code

  deadbeef = X'eadbeef' ! X'deadbeef' yields overflow error on 32bit machines
  fse = 4711
  onetwothree = 1234567809.0987654321D0
  pitest = 3.1415926535897932384626433832795029D0 ! from linux <math.h>
  
  is = 1
  ie = NoNx
  js = 1
  je = NoNy
 

  !----- read namelist -----
    
  open(file='input.nml', unit=42)
  ierr=1;
  do while (ierr /= 0)
     read (unit=42, nml=atmosphere_nml, iostat=io, end=10)
     !ierr = check_nml_error (io, 'atmosphere_nml')
  enddo
10 close(42)


  !----- write version and namelist to log file -----
  write (*, nml=atmosphere_nml)


  !---- compute physics/atmos time step in seconds ----

  ! initialize domain decomposition
  
  pe = 0
  root_pe = 0
  npes = 1

  !call mpp_get_compute_domain(grid_domain, is, ie, js, je)

  ! --- now we can allocate variables ---

  ! hack: in case Aeolus uses more than 1 tracer other than water vapor, revamp this code.
  allocate (Surf_diff%dtmass  (is:ie, js:je) )
  allocate (Surf_diff%dflux_t (is:ie, js:je) )
  allocate (Surf_diff%dflux_tr (is:ie, js:je, q_ind) )
  allocate (Surf_diff%delta_t (is:ie, js:je) )
  allocate (Surf_diff%delta_tr (is:ie, js:je, q_ind) )
  allocate (Surf_diff%delta_u (is:ie, js:je) )
  allocate (Surf_diff%delta_v (is:ie, js:je) )
    
  Surf_diff%dtmass  = 0.0
  Surf_diff%dflux_t = 0.0
  Surf_diff%dflux_tr = 0.0
  Surf_diff%delta_t = 0.0
  Surf_diff%delta_tr = 0.0
  Surf_diff%delta_u = 0.0
  Surf_diff%delta_v = 0.0

  Time_init_seconds = 0
  Time_init_days = 0
  Time_seconds = 7200
  Time_days = 42
  Time_step_seconds = 3600
  Time_step_days = 0

  allocate(testarray3d(2,3,4))
  do i=1,2
     do j=1,3
        do k=1,4
           testarray3d(i,j,k) = i*100.0+j+k/10.0
           write(*,*) 'testarray3d', i, j, k, testarray3d(i, j, k)
        end do
     end do
  end do

  ! now call out into the C++ - implemented heart of Aeolus
  write(*,*) 'MAXEXPONENT ', MAXEXPONENT(pitest)
  call aeolus_init(     &
       deadbeef, pitest,  &
       true, false,       &
       testarray3d,       &
       kind_i, kind_r, kind_l,            &
       Time_init_days, Time_init_seconds, &
       Time_days, Time_seconds,           &
       Time_step_days, Time_step_seconds, &
       NoNx, NoNy, NoNz, NoNz_Strato, &
       polar_merging,     &
       lon_0_360,         &
       pe, root_pe, npes, &
       is, ie, js, je,    &
       nr_tracers, q_ind, &
       onetwothree, fse   &
       )


  !----- initialize atmos_axes and diagnostics
    
  allocate(lonb(1:NoNx+1), latb(1:NoNy+1), lon(1:NoNx), lat(1:NoNy))

  call aeolus_get_axes_coords(lonb, latb, lon, lat)
  write (*,*) 'lonb', lonb
  write (*,*) 'latb', latb
  write (*,*) 'lon', lon
  write (*,*) 'lat', lat


  ! now init the tables of cosine zenith angle and daily insolation
  DaysPerYear = 365
  allocate(solarm(NoNy, DaysPerYear))
  allocate(coszm(NoNy, DaysPerYear))
  
  call sinsol_mod_init(real(Time_init_days / DaysPerYear, 8), DaysPerYear, NoNy, lat, solarm, coszm)

  deallocate(lonb)
  deallocate(latb)
  deallocate(lon)
  deallocate(lat)
    
  call test_inner(is,ie,js,je, Surf_diff)


  call aeolus_finish(Time_seconds, Time_days)

contains

  subroutine test_inner(is, ie, js, je, Surf_diff)
    integer, intent(in) :: is, ie, js, je
    type(surf_diff_type), intent(inout) :: Surf_diff

    integer :: i, j

    ! for ``down'' test
    real(8),  dimension(is:ie,js:je) :: frac_land, t_surf, albedo
    real(8),  dimension(is:ie,js:je) :: albedo_vis_dir, albedo_nir_dir
    real(8),  dimension(is:ie,js:je) :: albedo_vis_dif, albedo_nir_dif
    real(8),  dimension(is:ie,js:je) :: rough_mom
    real(8),  dimension(is:ie,js:je) :: u_star, b_star, q_star, dtau_du, dtau_dv
    real(8),  dimension(is:ie,js:je) :: u_flux, v_flux
    real(8),  dimension(is:ie,js:je) :: gust, coszen
    real(8),  dimension(is:ie,js:je) :: flux_sw
    real(8),  dimension(is:ie,js:je) :: flux_sw_dir, flux_sw_dif, flux_sw_down_vis_dir
    real(8),  dimension(is:ie,js:je) :: flux_sw_down_vis_dif, flux_sw_down_total_dir
    real(8),  dimension(is:ie,js:je) :: flux_sw_down_total_dif, flux_sw_vis
    real(8),  dimension(is:ie,js:je) :: flux_sw_vis_dir, flux_sw_vis_dif
    real(8),  dimension(is:ie,js:je) :: flux_lw


    ! for ``up'' test
    !real,                 dimension(is:ie,js:je) :: frac_land
    !type(surf_diff_type)                         :: Surf_diff
    real(8),                 dimension(is:ie,js:je) :: lprec, fprec !, gust
    !real,                 dimension(:,:)         :: u_star, b_star, q_star

    ! for cell_area test
    real(8), dimension(is:ie,js:je) :: area_out

    ! for boundary test
    real(8)         :: lon_boundaries(is:ie+1,js:je+1), lat_boundaries(is:ie+1,js:je+1)   ! radians !!

    ! for bottom_mass test
    real(8), dimension(is:ie,js:je) :: t_bot, p_bot, z_bot, p_surf, slp
    real(8), dimension(is:ie,js:je,3) :: tr_bot

    ! for wind test
    real(8), dimension(is:ie,js:je) :: u_bot, v_bot
  
    ! for stock test
    integer :: idx ! stock type index
    real(8) :: val

    do i=is,ie
       do j=js,je
          frac_land(i,j) = i+j/100.0
       end do
    end do

    ! get info from SphericalGrid
    call aeolus_get_boundary(lon_boundaries, lat_boundaries, .true., &
         lbound(lon_boundaries, dim=1), ubound(lon_boundaries, dim=1), &
         lbound(lon_boundaries, dim=2), ubound(lon_boundaries, dim=2), &
         lbound(lat_boundaries, dim=1), ubound(lat_boundaries, dim=1), &
         lbound(lat_boundaries, dim=2), ubound(lat_boundaries, dim=2))

    write(*,*) 'lon_boundaries ', lon_boundaries
    write(*,*) 'lat_boundaries ', lat_boundaries
    write(*,*) 'end boundaries'

    ! get area information from SphericalGrid
    call aeolus_cell_area(area_out)

    !do xsize=is,ie
    !   do ysize=js,je
    !      write(*,*) 'i ', xsize, ' j ', ysize, ' area ', area_out(xsize, ysize)
    !   end do
    !end do

    call aeolus_test_var_access(         &
         frac_land,                        &
         t_bot, tr_bot, p_bot)


    do i=is,ie
       do j=js,je
          write(*,*) 't_bot ', i, j, t_bot(i,j)
          write(*,*) 'p_bot ', i, j, p_bot(i,j)
          write(*,*) 'tr_bot', i, j, tr_bot(i,j,1)
       end do
    end do
    call flush(6)

    !call aeolus_atmosphere_down (        &
    !     Time_seconds, Time_days,          &
    !     frac_land, t_surf,  albedo,       &
    !     albedo_vis_dir, albedo_nir_dir,   &
    !     albedo_vis_dif, albedo_nir_dif,   &
    !     rough_mom,                        &
    !     u_star,  b_star, q_star,          &
    !     dtau_du, dtau_dv, u_flux, v_flux, &
    !     gust, coszen, flux_sw,            &
    !     flux_sw_dir, flux_sw_dif,         &
    !     flux_sw_down_vis_dir,             &
    !     flux_sw_down_vis_dif,             &
    !     flux_sw_down_total_dir,           &
    !     flux_sw_down_total_dif,           &
    !     flux_sw_vis,                      &
    !     flux_sw_vis_dir,                  &
    !     flux_sw_vis_dif,                  &
    !     flux_lw,                          &
    !     Surf_diff%delta_tr,               &
    !     Surf_diff%dflux_tr,               &
    !     Surf_diff%delta_t,                &
    !     Surf_diff%dflux_t,                &
    !     Surf_diff%dtmass,                 &
    !     Surf_diff%delta_u,                &
    !     Surf_diff%delta_v,                &
    !     )



    !call aeolus_atmosphere_up ( &
    !     Time_seconds, Time_days, &
    !     frac_land,               &
    !     Surf_diff%delta_tr,      &
    !     Surf_diff%dflux_tr,      &
    !     lprec, fprec, gust,      &
    !     u_star, b_star, q_star )

    !call aeolus_get_bottom_mass(t_bot, tr_bot, p_bot, z_bot, p_surf, slp)

    !call aeolus_get_bottom_wind(u_bot, v_bot)
    
    !call aeolus_get_stock_pe(idx, val)

  end subroutine test_inner


end program testaeolus

