! (c) British Crown Copyright 2008, the Met Office.
! All rights reserved.
! $Revision: 88 $, $Date: 2013-11-13 07:08:38 -0700 (Wed, 13 Nov 2013) $
! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/cosp_test.F90 $
! 
! Redistribution and use in source and binary forms, with or without modification, are permitted 
! provided that the following conditions are met:
! 
!     * Redistributions of source code must retain the above copyright notice, this list 
!       of conditions and the following disclaimer.
!     * Redistributions in binary form must reproduce the above copyright notice, this list
!       of conditions and the following disclaimer in the documentation and/or other materials 
!       provided with the distribution.
!     * Neither the name of the Met Office nor the names of its contributors may be used 
!       to endorse or promote products derived from this software without specific prior written 
!       permission.
! 
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 
! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND 
! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 
! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER 
! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT 
! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

!
! History:
! Feb 2008 - A. Bodas-Salcedo - Initial version
! Dec 2010 - A. Bodas-Salcedo - Added capability for processing multiple files
!

#include "cosp_defs.h"
PROGRAM COSPTEST
  USE MOD_COSP_CONSTANTS
  USE MOD_COSP_TYPES
  USE MOD_COSP
  USE MOD_COSP_IO
  IMPLICIT NONE

  ! Local variables
  character(len=64)  :: cosp_input_nl='cosp_input_nl.txt'
!   character(len=64)  :: cosp_output_nl='cfmip2/cosp_output_cfmip2_short_offline.txt'
  character(len=64)  :: cosp_output_nl='cosp_output_nl.txt'
  character(len=512) :: dinput ! Directory with input files
  integer,parameter :: N_MAX_INPUT_FILES = 10000 ! Maximum number of input files
  character(len=64),dimension(N_MAX_INPUT_FILES) :: finput ! File names
  character(len=600) :: dfinput ! Input file
  character(len=512) :: cmor_nl
  character(len=8)  :: wmode ! Writing mode 'replace' or 'append'
  integer :: overlap   !  overlap type: 1=max, 2=rand, 3=max/rand
  integer :: isccp_topheight,isccp_topheight_direction
  integer :: Ncolumns     ! Number of subcolumns in SCOPS
  integer :: Npoints      ! Number of gridpoints
  integer :: Nlevels      ! Number of levels
  integer :: Nlr          ! Number of levels in statistical outputs
  integer :: Npoints_it   ! Max number of gridpoints to be processed in one iteration
  integer :: Nfiles       ! Number of files to be processed
  integer,parameter :: ntsteps=5 
  integer :: i,k
  type(cosp_config) :: cfg   ! Configuration options
  type(cosp_gridbox) :: gbx ! Gridbox information. Input for COSP
  type(cosp_subgrid) :: sgx     ! Subgrid outputs
  type(cosp_sgradar) :: sgradar ! Output from radar simulator
  type(cosp_sglidar) :: sglidar ! Output from lidar simulator
  type(cosp_isccp)   :: isccp   ! Output from ISCCP simulator
  type(cosp_modis)   :: modis   ! Output from MODIS simulator
  type(cosp_misr)    :: misr    ! Output from MISR simulator
  type(cosp_rttov)   :: rttov   ! Output from RTTOV 
  type(cosp_vgrid)   :: vgrid   ! Information on vertical grid of stats
  type(cosp_radarstats) :: stradar ! Summary statistics from radar simulator
  type(cosp_lidarstats) :: stlidar ! Summary statistics from lidar simulator
  type(var1d) :: v1d(N1D+1) ! Structures needed by output routines for 1D variables
  type(var2d) :: v2d(N2D) ! Structures needed by output routines for 2D variables
  type(var3d) :: v3d(N3D) ! Structures needed by output routines for 3D variables
  integer ::  grid_id,latvar_id,lonvar_id,lon_axid,lat_axid,time_axid,height_axid,height_mlev_axid,column_axid,sza_axid, &
              temp_axid,channel_axid,dbze_axid,sratio_axid,MISR_CTH_axid,tau_axid,pressure2_axid,ReffLiq_axid,ReffIce_axid
  real,dimension(:),allocatable :: lon,lat
  real,dimension(:,:),allocatable,target :: p,ph,zlev,zlev_half,T,sh,rh,tca,cca, &
                    mr_lsliq,mr_lsice,mr_ccliq,mr_ccice,fl_lsrain,fl_lssnow,fl_lsgrpl, &
                    fl_ccrain,fl_ccsnow,dtau_s,dtau_c,dem_s,dem_c,mr_ozone
  real,dimension(:,:,:),allocatable :: Reff
  real,dimension(:),allocatable :: skt,landmask,u_wind,v_wind,sunlit
  integer :: t0,t1,t2,t3,count_rate,count_max
  integer :: Nlon,Nlat,geomode
  real :: radar_freq,k2,ZenAng,co2,ch4,n2o,co,emsfc_lw
  integer,dimension(RTTOV_MAX_CHANNELS) :: Channels
  real,dimension(RTTOV_MAX_CHANNELS) :: Surfem
  integer :: surface_radar,use_mie_tables,use_gas_abs,do_ray,melt_lay
  integer :: Nprmts_max_hydro,Naero,Nprmts_max_aero,lidar_ice_type
  integer :: platform,satellite,Instrument,Nchannels,N1
  logical :: use_vgrid,csat_vgrid,use_precipitation_fluxes,use_reff
  double precision :: time,time_bnds(2),time_step
  real :: toffset_step,half_time_step
  namelist/COSP_INPUT/cmor_nl,overlap,isccp_topheight,isccp_topheight_direction, &
              npoints,npoints_it,ncolumns,nlevels,use_vgrid,nlr,csat_vgrid,dinput,finput, &
              radar_freq,surface_radar,use_mie_tables, &
              use_gas_abs,do_ray,melt_lay,k2,Nprmts_max_hydro,Naero,Nprmts_max_aero, &
              lidar_ice_type,use_precipitation_fluxes,use_reff, &
              platform,satellite,Instrument,Nchannels, &
              Channels,Surfem,ZenAng,co2,ch4,n2o,co

  !---------------- End of declaration of variables --------------
   

  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  ! Read COSP namelists
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  finput(:) = ''
  open(10,file=cosp_input_nl,status='old')
  read(10,nml=cosp_input)
  close(10)
  call read_cosp_output_nl(cosp_output_nl,cfg)
  
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  ! Find number of input files
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  i=1
  do while (i <= N_MAX_INPUT_FILES)
     if (len_trim(finput(i)) < 1) exit
     i=i+1
  enddo
  Nfiles = i-1
  if (Nfiles < 1) call cosp_error('cosp_test','Number of files < 1')
  if (Nfiles > N_MAX_INPUT_FILES) call cosp_error('cosp_test','Number of files > N_MAX_INPUT_FILES')

  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  ! Allocate local arrays
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  allocate(lon(Npoints),lat(Npoints), &
           p(Npoints,Nlevels),ph(Npoints,Nlevels), &
           zlev(Npoints,Nlevels),zlev_half(Npoints,Nlevels),T(Npoints,Nlevels), &
           sh(Npoints,Nlevels),rh(Npoints,Nlevels), &
           tca(Npoints,Nlevels),cca(Npoints,Nlevels),mr_lsliq(Npoints,Nlevels), &
           mr_lsice(Npoints,Nlevels),mr_ccliq(Npoints,Nlevels),mr_ccice(Npoints,Nlevels), &
           fl_lsrain(Npoints,Nlevels),fl_lssnow(Npoints,Nlevels),fl_lsgrpl(Npoints,Nlevels),fl_ccrain(Npoints,Nlevels), &
           fl_ccsnow(Npoints,Nlevels),Reff(Npoints,Nlevels,N_HYDRO),dtau_s(Npoints,Nlevels),dtau_c(Npoints,Nlevels), &
           dem_s(Npoints,Nlevels),dem_c(Npoints,Nlevels),skt(Npoints),landmask(Npoints), &
           mr_ozone(Npoints,Nlevels),u_wind(Npoints),v_wind(Npoints),sunlit(Npoints))

  call system_clock(t0,count_rate,count_max) !!! Only for testing purposes

  ! Example that processes ntsteps. It always uses the same input data
  time_step      = 3.D0/24.D0
  time           = 8*1.D0/8.D0 ! First time step
  toffset_step   = time_step/Npoints
  half_time_step = 0.5*time_step
  do i=1,Nfiles
        dfinput=trim(dinput)//trim(finput(i))
        time_bnds = (/time-half_time_step,time+half_time_step/) ! This may need to be adjusted, 
                                                                ! depending on the approx_interval in the MIP table
        print *, 'Processing file: ', trim(dfinput)
        print *, 'Time: ',time
        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        ! Read input geophysical variables from NetCDF file
        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        ! input : surface to top
        call nc_read_input_file(dfinput,Npoints,Nlevels,N_HYDRO,lon,lat,p,ph,zlev,zlev_half,T,sh,rh,tca,cca, &
                mr_lsliq,mr_lsice,mr_ccliq,mr_ccice,fl_lsrain,fl_lssnow,fl_lsgrpl,fl_ccrain,fl_ccsnow,Reff, &
                dtau_s,dtau_c,dem_s,dem_c,skt,landmask,mr_ozone,u_wind,v_wind,sunlit, &
                emsfc_lw,geomode,Nlon,Nlat)
                ! geomode = 2 for (lon,lat) mode.
                ! geomode = 3 for (lat,lon) mode.
                ! In those modes it returns Nlon and Nlat with the correct values

        call system_clock(t1,count_rate,count_max)
        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        ! Allocate memory for gridbox type
        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        print *, 'Allocating memory for gridbox type...'
        call construct_cosp_gridbox(time,time_bnds,radar_freq,surface_radar,use_mie_tables,use_gas_abs, &
                                    do_ray,melt_lay,k2, &
                                    Npoints,Nlevels,Ncolumns,N_HYDRO,Nprmts_max_hydro,Naero,Nprmts_max_aero,Npoints_it, &
                                    lidar_ice_type,isccp_topheight,isccp_topheight_direction,overlap,emsfc_lw, &
                                    use_precipitation_fluxes,use_reff, &
                                    Platform,Satellite,Instrument,Nchannels,ZenAng, &
                                    channels(1:Nchannels),surfem(1:Nchannels),co2,ch4,n2o,co,gbx)

        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        ! Here code to populate input structure
        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        print *, 'Populating input structure...'
        gbx%longitude = lon
        gbx%latitude = lat
        ! Toffset. This assumes that time is the mid-point of the interval.
        do k=1,Npoints
          gbx%toffset(k) = -half_time_step + toffset_step*(k-0.5)
        enddo
        gbx%p = p
        gbx%ph = ph
        gbx%zlev = zlev
        gbx%zlev_half = zlev_half
        gbx%T = T
        gbx%q = rh
        gbx%sh = sh
        gbx%cca = cca
        gbx%tca = tca
        gbx%psfc = ph(:,1)
        gbx%skt  = skt
        gbx%land = landmask
        gbx%mr_ozone  = mr_ozone
        gbx%u_wind  = u_wind
        gbx%v_wind  = v_wind
        gbx%sunlit  = sunlit

        gbx%mr_hydro(:,:,I_LSCLIQ) = mr_lsliq
        gbx%mr_hydro(:,:,I_LSCICE) = mr_lsice
        gbx%mr_hydro(:,:,I_CVCLIQ) = mr_ccliq
        gbx%mr_hydro(:,:,I_CVCICE) = mr_ccice
        gbx%rain_ls = fl_lsrain
        gbx%snow_ls = fl_lssnow
        gbx%grpl_ls = fl_lsgrpl
        gbx%rain_cv = fl_ccrain
        gbx%snow_cv = fl_ccsnow

        gbx%Reff = Reff
        gbx%Reff(:,:,I_LSRAIN) = 0.0

        ! ISCCP simulator
        gbx%dtau_s   = dtau_s
        gbx%dtau_c   = dtau_c
        gbx%dem_s    = dem_s
        gbx%dem_c    = dem_c

        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        ! Define new vertical grid
        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        print *, 'Defining new vertical grid...'
        call construct_cosp_vgrid(gbx,Nlr,use_vgrid,csat_vgrid,vgrid)

        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        ! Allocate memory for other types
        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        print *, 'Allocating memory for other types...'
        call construct_cosp_subgrid(Npoints, Ncolumns, Nlevels, sgx)
        call construct_cosp_sgradar(cfg,Npoints,Ncolumns,Nlevels,N_HYDRO,sgradar)
        call construct_cosp_radarstats(cfg,Npoints,Ncolumns,vgrid%Nlvgrid,N_HYDRO,stradar)
        call construct_cosp_sglidar(cfg,Npoints,Ncolumns,Nlevels,N_HYDRO,PARASOL_NREFL,sglidar)
        call construct_cosp_lidarstats(cfg,Npoints,Ncolumns,vgrid%Nlvgrid,N_HYDRO,PARASOL_NREFL,stlidar)
        call construct_cosp_isccp(cfg,Npoints,Ncolumns,Nlevels,isccp)
        call construct_cosp_modis(cfg,Npoints,modis)
        call construct_cosp_misr(cfg,Npoints,misr)
        call construct_cosp_rttov(cfg,Npoints,Nchannels,rttov)

        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        ! Call simulator
        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        print *, 'Calling simulator...'
#ifdef RTTOV
        call cosp(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,rttov,stradar,stlidar)
#else
        call cosp(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,stradar,stlidar)
#endif
        call system_clock(t2,count_rate,count_max)

        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        ! Write outputs to CMOR-compliant NetCDF
        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        ! Initialise CMOR interface
        gbx%time = time
        ! Write one time step to file
        if (cfg%Lwrite_output) then
            N1 = N1D
            if (geomode == 1) N1 = N1D+1
            print *, 'Writing outputs...'
            if (i == 1) call nc_cmor_init(cmor_nl,'replace',cfg,vgrid,gbx,sgx,sgradar,sglidar, &
                                  isccp,misr,modis,rttov,stradar,stlidar,geomode,Nlon,Nlat,N1,N2D,N3D, &
                                  lon_axid,lat_axid,time_axid,height_axid,height_mlev_axid,grid_id,lonvar_id,latvar_id, &
                                  column_axid,sza_axid,temp_axid,channel_axid,dbze_axid, &
                                  sratio_axid,MISR_CTH_axid,tau_axid,pressure2_axid,ReffLiq_axid,ReffIce_axid, &
                                  v1d(1:N1),v2d,v3d)
            if (geomode == 1) then
               print *, 'Associate'
               call nc_cmor_associate_1d(grid_id,time_axid,height_axid,height_mlev_axid,column_axid,sza_axid, &
                         temp_axid,channel_axid,dbze_axid,sratio_axid,MISR_CTH_axid,tau_axid,pressure2_axid,ReffLiq_axid,ReffIce_axid, &
                         Nlon,Nlat,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,rttov,stradar,stlidar, &
                         v1d(1:N1),v2d,v3d)
               print *, 'Write'
               call nc_cmor_write_1d(gbx,time_bnds,lonvar_id,latvar_id,N1,N2D,N3D,v1d(1:N1),v2d,v3d)
            elseif (geomode >  1) then
               call nc_cmor_associate_2d(lon_axid,lat_axid,time_axid,height_axid,height_mlev_axid,column_axid,sza_axid, &
                         temp_axid,channel_axid,dbze_axid,sratio_axid,MISR_CTH_axid,tau_axid,pressure2_axid,ReffLiq_axid,ReffIce_axid, &
                         Nlon,Nlat,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,rttov,stradar,stlidar, &
                         v1d(1:N1),v2d,v3d)
               call nc_cmor_write_2d(time,time_bnds,geomode,Nlon,Nlat,N1,N2D,N3D,v1d(1:N1),v2d,v3d)

            endif
            if (i == Nfiles) then
              call nc_cmor_close()
            endif
        endif

        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        ! Deallocate memory in derived types
        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        print *, 'Deallocating memory...'
        call free_cosp_gridbox(gbx)
        call free_cosp_subgrid(sgx)
        call free_cosp_sgradar(sgradar)
        call free_cosp_radarstats(stradar)
        call free_cosp_sglidar(sglidar)
        call free_cosp_lidarstats(stlidar)
        call free_cosp_isccp(isccp)
        call free_cosp_misr(misr)
        call free_cosp_modis(modis)
        call free_cosp_rttov(rttov)
        call free_cosp_vgrid(vgrid)
        ! Update time
        time = time + time_step
  enddo
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  ! Deallocate memory in local arrays
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  deallocate(lon,lat,p,ph,zlev,zlev_half,T,sh,rh,tca,cca, mr_lsliq,mr_lsice,mr_ccliq,mr_ccice, &
           fl_lsrain,fl_lssnow,fl_lsgrpl,fl_ccrain,fl_ccsnow,Reff,dtau_s,dtau_c,dem_s,dem_c,skt, &
           landmask,mr_ozone,u_wind,v_wind,sunlit)

  ! Time in s. Only for testing purposes
  call system_clock(t3,count_rate,count_max)
  do i=1,N_SIMULATORS
      print *,'=== '//trim(SIM_NAME(i))//': ', float(tsim(i))/count_rate
  enddo
  print *,'=== COSP: ', (t2-t1)*1.0/count_rate
  print *,'=== TOTAL: ', (t3-t0)*1.0/count_rate
    
END PROGRAM COSPTEST