program kcm1d

  use radinc_h,      only: NAERKIND
  use radcommon_h,   only: L_NSPECTI, L_NSPECTV, sigma
  use watercommon_h, only: mH2O
  use ioipsl_getincom 

  implicit none

!==================================================================
!     
!     Purpose
!     -------
!     Run the universal model radiative transfer once in a 1D column.
!     Useful for climate evolution studies etc.
!     
!     It can be compiled with a command like (e.g. for 25 layers):
!                                  "makegcm -p std -d 25 kcm1d"
!     It requires the files "callphys.def", "gases.def"
!     "traceur.def", and "run.def" with a line "INCLUDEDEF=callphys.def"
!
!     Authors
!     -------
!     R. Wordsworth
!
!==================================================================

#include "dimensions.h"
#include "dimphys.h"
#include "callkeys.h"
#include "comcstfi.h"
#include "planete.h"
#include "comsaison.h"
#include "control.h"
#include "advtrac.h"

! --------------------------------------------------------------
!  Declarations
! --------------------------------------------------------------

  integer nlayer,nlevel,nq
  integer ilay,ilev,iq,iw,iter
  real play(nlayermx)     ! Pressure at the middle of the layers [Pa]
  real zlay(nlayermx)     ! Altitude at middle of the layers [km]
  real plev(nlayermx+1)   ! Intermediate pressure levels [Pa]
  real temp(nlayermx)     ! temperature at the middle of the layers [K]
  real q(nlayermx,nqmx)   ! tracer mixing ratio [kg/kg]
  real vmr(nlayermx,nqmx) ! tracer mixing ratio [mol/mol]
  real qsurf(nqmx)        ! tracer surface budget [kg/kg] ????
  real psurf,psurf_n,tsurf      

  real emis, albedo

  real cpp3D(nlayermx)   ! specific heat capacity at const. pressure
  real rcp3D(nlayermx)   ! R / specific heat capacity
  real muvar(nlayermx+1)

  real dtsw(nlayermx) ! heating rate (K/s) due to SW
  real dtlw(nlayermx) ! heating rate (K/s) due to LW
  real fluxsurf_lw   ! incident LW flux to surf (W/m2)
  real fluxtop_lw    ! outgoing LW flux to space (W/m2)
  real fluxsurf_sw   ! incident SW flux to surf (W/m2)
  real fluxabs_sw    ! SW flux absorbed by planet (W/m2)
  real fluxtop_dn    ! incident top of atmosphere SW flux (W/m2)
      
! not used
  real reffrad(nlayermx,naerkind)
  real nueffrad(nlayermx,naerkind)
  real cloudfrac(nlayermx)
  real totcloudfrac
  real tau_col

  real dTstrat
  real aerosol(nlayermx,naerkind) ! aerosol tau (kg/kg)
  real OLR_nu(L_NSPECTI)
  real GSR_nu(L_NSPECTV)
  real Eatmtot

  integer ierr
  logical firstcall,lastcall

! --------------
! Initialisation
! --------------

  pi=2.E+0*asin(1.E+0)

  reffrad(:,:)  = 0.0
  nueffrad(:,:) = 0.0
  cloudfrac(:)  = 0.0
  totcloudfrac  = 0.0

!  Load parameters from "run.def" 
! -------------------------------

! check if 'run1d.def' file is around (otherwise reading parameters
! from callphys.def via getin() routine won't work.)
  open(90,file='run.def',status='old',form='formatted',&
       iostat=ierr)
  if (ierr.ne.0) then
     write(*,*) 'Cannot find required file "run.def"'
     write(*,*) '  (which should contain some input parameters'
     write(*,*) '   along with the following line:'
     write(*,*) '   INCLUDEDEF=callphys.def'
     write(*,*) '   )'
     write(*,*) ' ... might as well stop here ...'
     stop
  else
     close(90)
  endif

  nlayer=nlayermx
  nlevel=nlayer+1

  psurf_n=100000. ! default value for psurf
  print*,'Dry surface pressure (Pa)?'
  call getin("psurf",psurf_n)
  write(*,*) " psurf = ",psurf_n
      
  tsurf=300.0 ! default value for tsurf
  print*,'Surface temperature (K)?'
  call getin("tref",tsurf)
  write(*,*) " tsurf = ",tsurf

  g=10.0 ! default value for g
  print*,'Gravity ?'
  call getin("g",g)
  write(*,*) " g = ",g
  
  periastr = 1.0
  apoastr  = 1.0
  print*,'Periastron (AU)?'
  call getin("periastr",periastr)
  write(*,*) "periastron = ",periastr
  dist_star = periastr 
  fract     = 0.5
  print*,'Apoastron (AU)?'
  call getin("apoastr",apoastr)
  write(*,*) "apoastron = ",apoastr
 
  periastr = periastr*149.598 ! AU to Mkm
  apoastr  = apoastr*149.598  ! AU to Mk

  albedo=0.2 ! default value for albedo
  print*,'Albedo of bare ground?'
  call getin("albedo",albedo)
  write(*,*) " albedo = ",albedo

  emis=1.0 ! default value for emissivity
  PRINT *,'Emissivity of bare ground ?'
  call getin("emis",emis)
  write(*,*) " emis = ",emis

  pceil=100.0 ! Pascals
  PRINT *,'Ceiling pressure (Pa) ?'
  call getin("pceil",pceil)
  write(*,*) " pceil = ", pceil

  mugaz=0.
  cpp=0.

  call su_gases
  call calc_cpp_mugaz


  call inifis(1,llm,0,86400.0,1.0,0.0,0.0,1.0,rad,g,r,cpp)

! Tracer initialisation
! ---------------------
  if (tracer) then
     ! load tracer names from file 'traceur.def'
     open(90,file='traceur.def',status='old',form='formatted',&
          iostat=ierr)
     if (ierr.eq.0) then
        write(*,*) "kcm1d: Reading file traceur.def"
        ! read number of tracers:
        read(90,*,iostat=ierr) nq
        if (ierr.ne.0) then
           write(*,*) "kcm1d: error reading number of tracers"
           write(*,*) "   (first line of traceur.def) "
           stop
        else
           ! check that the number of tracers is indeed nqmx
           if (nq.ne.nqmx) then
              write(*,*) "kcm1d: error, wrong number of tracers:"
              write(*,*) "nq=",nq," whereas nqmx=",nqmx
              stop
           endif
        endif
        
        do iq=1,nq
           ! minimal version, just read in the tracer names, 1 per line
           read(90,*,iostat=ierr) tnom(iq)
           if (ierr.ne.0) then
              write(*,*) 'kcm1d: error reading tracer names...'
              stop
           endif
        enddo !of do iq=1,nq
     endif
  !endif
  
  call initracer()

  endif


  do iq=1,nqmx
     do ilay=1,nlayer
        q(ilay,iq) = 0.
     enddo
  enddo
  
  do iq=1,nqmx
     qsurf(iq) = 0.
  enddo

  firstcall = .true.
  lastcall  = .false.

  iter    = 1
  Tstrat  = 60.0
  dTstrat = 1000.0

! ---------
! Run model
! ---------
  do
     psurf = psurf_n

!    Create vertical profiles
     call kcmprof_fn(psurf,qsurf(1),tsurf,    &
          tstrat,play,plev,zlay,temp,q(:,1),muvar(1))

!    Run radiative transfer
     call callcorrk(ngridmx,nlayer,q,nqmx,qsurf,      &
          albedo,emis,mu0,plev,play,temp,                    &
          tsurf,fract,dist_star,aerosol,cpp3D,muvar,         &
          dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
          fluxabs_sw,fluxtop_dn,reffrad,tau_col,  &
          cloudfrac,totcloudfrac,.false.,firstcall,lastcall)
 
!    Iterate stratospheric temperature
     print*,'Tstrat = ',Tstrat
     dTstrat = Tstrat
     !Tstrat  = Tsurf*(0.2786*(psurf/100000.)**(-1.123))**0.25 
     ! skin temperature (gray approx.) using analytic pure H2 expression
     !Tstrat  = (fluxabs_sw/(2*sigma))**0.25 ! skin temperature (gray approx.)
     Tstrat  = (fluxtop_lw/(2*sigma))**0.25 ! skin temperature (gray approx.)
     dTstrat = dTstrat-Tstrat

     if(abs(dTstrat).lt.1.0)then
        print*,'dTstrat = ',dTstrat
        exit
     endif

     iter=iter+1
     if(iter.eq.100)then
        print*,'Stratosphere failed to converge after'
        print*,'100 iteration, aborting run.'
        call abort
     endif

  end do

! Calculate total atmospheric energy
  Eatmtot=0.0
!  call calcenergy_kcm(tsurf,temp,play,plev,qsurf,&
!     q(:,1),muvar,Eatmtot)

! ------------------------
! Save data to ascii files
! ------------------------

  print*,'Saving profiles...'
  open(115,file='profpres.out',form='formatted')
  open(116,file='proftemp.out',form='formatted')
  open(117,file='profztab.out',form='formatted')
  open(118,file='profqvar.out',form='formatted')
  open(119,file='profvvar.out',form='formatted')

  write(115,*) psurf
  write(116,*) tsurf
  write(117,*) 0.0
  write(118,*) qsurf(1)
  write(119,*) qsurf(1)*(muvar(1)/mH2O)
  do ilay=1,nlayer
     vmr(ilay,1) = q(ilay,1)*(muvar(ilay+1)/mH2O)
     write(115,*) play(ilay)
     write(116,*) temp(ilay)
     write(117,*) zlay(ilay)
     write(118,*) q(ilay,1)
     write(119,*) vmr(ilay,1)
  enddo
  close(115)
  close(116)
  close(117)
  close(118)
  close(119)

  print*,'Saving scalars...'
  open(116,file='surf_vals.out')
  write(116,*) tsurf,psurf,fluxtop_dn,         &
       fluxabs_sw,fluxtop_lw 
  close(116)
  open(111,file='ene_vals.out')
  write(111,*) tsurf,psurf,Eatmtot,Tstrat
  close(111)

  print*,'Saving spectra...'
  open(117,file='OLRnu.out')
  do iw=1,L_NSPECTI
     write(117,*) OLR_nu(iw)
  enddo
  close(117)
  
  open(127,file='GSRnu.out')
  do iw=1,L_NSPECTV
     write(127,*) GSR_nu(iw)
  enddo
  close(127)  

end program kcm1d