PROGRAM testphys1d

use comsoil_h,           only: inertiedat, inertiesoil, nsoilmx, tsoil, nqsoil, qsoil
use surfdat_h,           only: albedodat, perennial_co2ice, watercap, tsurf, emis, qsurf
use comslope_mod,        only: def_slope, subslope_dist
use phyredem,            only: physdem0, physdem1
use watersat_mod,        only: watersat
use tracer_mod,          only: igcm_h2o_vap, igcm_h2o_ice, noms
use comcstfi_h,          only: pi, g, rcp, cpp
use time_phylmdz_mod,    only: daysec
use dimradmars_mod,      only: tauvis, totcloudfrac, albedo
use dust_param_mod,      only: tauscaling
use comvert_mod,         only: ap, bp, aps, bps, pa, preff
use physiq_mod,          only: physiq
use turb_mod,            only: q2
use write_output_mod,    only: write_output
use ioipsl_getincom,     only: getin ! To use 'getin'
use init_testphys1d_mod, only: init_testphys1d
use writerestart1D_mod,  only: writerestart1D
use parse_args_mod,      only: parse_args
use callkeys_mod,        only: water
! Mostly for XIOS outputs
use mod_const_mpi,       only: init_const_mpi
use parallel_lmdz,       only: init_parallel

implicit none

!=======================================================================
!   subject:
!   --------
!   PROGRAM useful to run physical part of the martian GCM in a 1D column
!
! Can be compiled with a command like (e.g. for 25 layers)
!   "makegcm -p mars -d 25 testphys1d"
! It requires the files "testphys1d.def" "callphys.def"
!   and a 'run.def' file (containing a "INCLUDEDEF=callphys.def" line)
!   and a file describing the sigma layers (e.g. "z2sig.def")
!
!   author: Frederic Hourdin, R. Fournier, F. Forget
!   -------
!
!   update: 12/06/2003, including chemistry (S. Lebonnois)
!                       and water ice (F. Montmessin)
!           27/09/2023, upgrade to F90 (JB Clément)
!
!=======================================================================

include "dimensions.h"

!--------------------------------------------------------------
! Declarations
!--------------------------------------------------------------
integer, parameter                  :: ngrid = 1     ! (2+(jjm-1)*iim - 1/jjm)
integer, parameter                  :: nlayer = llm
real, parameter                     :: odpref = 610. ! DOD reference pressure (Pa)
integer                             :: ndt, ilayer, idt
logical                             :: firstcall, lastcall
integer                             :: day0          ! initial (sol ; =0 at Ls=0)
real                                :: day           ! date during the run
real                                :: time          ! time (0<time<1 ; time=0.5 a midi)
real                                :: dttestphys    ! testphys1d timestep
real, dimension(nlayer)             :: play          ! Pressure at the middle of the layers (Pa)
real, dimension(nlayer + 1)         :: plev          ! intermediate pressure levels (Pa)
real                                :: psurf         ! Surface pressure
real, dimension(nlayer)             :: u, v          ! zonal, meridional wind
real                                :: gru, grv      ! prescribed "geostrophic" background wind
real, dimension(nlayer)             :: temp          ! temperature at the middle of the layers
real, dimension(:,:,:), allocatable :: q             ! tracer mixing ratio (e.g. kg/kg)
real, dimension(1)                  :: wstar = 0.    ! Thermals vertical velocity

! Physical and dynamical tendencies (e.g. m.s-2, K/s, Pa/s)
real, dimension(nlayer)             :: du, dv, dtemp
real, dimension(1)                  :: dpsurf
real, dimension(:,:,:), allocatable :: dq, dqdyn

! Various intermediate variables
real                    :: zls, pks, ptif
real, dimension(nlayer) :: phi, h, s, w
integer                 :: nq = 1 ! number of tracers
real, dimension(1)      :: latitude, longitude, cell_area

! Use of starting files for 1D
logical :: startfiles_1D, therestart1D, therestartfi, there

! Prescription of atmospheric water/ice profiles
logical                         :: ctrl_h2ovap, ctrl_h2oice
real                            :: relaxtime_h2ovap, relaxtime_h2oice
real, dimension(nlayer)         :: qref_h2ovap, qref_h2oice
real, dimension(:), allocatable :: zqsat
!=======================================================================

!=======================================================================
! INITIALISATION
!=======================================================================
! Parse command-line options
call parse_args()

#ifdef CPP_XIOS
    call init_const_mpi
    call init_parallel
#endif

! Initialize "serial/parallel" related stuff
!call init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
!call init_phys_lmdz(1,1,llm,1,(/1/))
!call initcomgeomphy

startfiles_1D = .false.
!------------------------------------------------------
! Loading run parameters from "run.def" file
!------------------------------------------------------
! check if 'run.def' file is around. Otherwise reading parameters
! from callphys.def via getin() routine won't work.
inquire(file = 'run.def',exist = there)
if (.not. there) then
    write(*,*) 'Cannot find required file "run.def"'
    write(*,*) '  (which should contain some input parameters along with the following line: INCLUDEDEF=callphys.def)'
    write(*,*) ' ... might as well stop here ...'
    error stop
endif

startfiles_1D = .false. ! Default value
write(*,*) 'Do you want to use starting files and/or to write restarting files?'
call getin('startfiles_1D',startfiles_1D)
write(*,*) 'startfiles_1D =', startfiles_1D

therestart1D = .false. ! Default value
therestartfi = .false. ! Default value
if (startfiles_1D) then
    inquire(file = 'start1D.txt',exist = therestart1D)
    if (.not. therestart1D) then
        write(*,*) 'There is no "start1D.txt" file!'
        write(*,*) 'Initialization is done with default values.'
    endif
    inquire(file = 'startfi.nc',exist = therestartfi)
    if (.not. therestartfi) then
        write(*,*) 'There is no "startfi.nc" file!'
        write(*,*) 'Initialization is done with default values.'
    endif
endif
if (therestart1D .and. .not. therestartfi) error stop 'There is a "start1D.txt" but no "startfi.nc". There might be problems for the initialization (for example ''qsurf'')!'

call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,odpref,nq,q, &
                     time,psurf,u,v,temp,ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, &
                     play,plev,latitude,longitude,cell_area,                                        &
                     ctrl_h2ovap,relaxtime_h2ovap,qref_h2ovap,ctrl_h2oice,relaxtime_h2oice,qref_h2oice)

! Write a "startfi" file
! ----------------------
! This file will be read during the first call to "physiq".
! It is needed to transfert physics variables to "physiq"...
if (.not. therestartfi) then
    call physdem0("startfi.nc",longitude,latitude,nsoilmx,ngrid, &
                  llm,nq,dttestphys,float(day0),0.,cell_area,    &
                  albedodat,inertiedat,def_slope,subslope_dist)
    call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,nqsoil,dttestphys,time,      &
                  tsurf,tsoil,inertiesoil,albedo,emis,q2,qsurf,qsoil,tauscaling, &
                  totcloudfrac,wstar,watercap,perennial_co2ice)
endif !(.not. therestartfi)

!=======================================================================
!  1D MODEL TIME STEPPING LOOP
!=======================================================================
firstcall = .true.
lastcall = .false.
do idt = 1,ndt
    if (idt == ndt) lastcall = .true.
!    if (idt == ndt - day_step - 1) then ! test
!        lastcall = .true.
!        call solarlong(day*1.0,zls)
!        write(103,*) 'Ls=',zls*180./pi
!        write(103,*) 'Lat=', latitude(1)*180./pi
!        write(103,*) 'Tau=', tauvis/odpref*psurf
!        write(103,*) 'RunEnd - Atmos. Temp. File'
!        write(103,*) 'RunEnd - Atmos. Temp. File'
!        write(104,*) 'Ls=',zls*180./pi
!        write(104,*) 'Lat=', latitude(1)
!        write(104,*) 'Tau=', tauvis/odpref*psurf
!        write(104,*) 'RunEnd - Atmos. Temp. File'
!    endif

    ! Compute geopotential
    ! --------------------
    s = (aps/psurf + bps)**rcp
    h = cpp*temp/(pks*s)

    phi(1) = pks*h(1)*(1. - s(1))
    do ilayer = 2,nlayer
        phi(ilayer) = phi(ilayer - 1) + pks*(h(ilayer - 1) + h(ilayer))*.5*(s(ilayer - 1) - s(ilayer))
    enddo

    ! Call physics
    ! ------------
    call physiq(1,llm,nq,firstcall,lastcall,day,time,dttestphys,plev,play,phi,u,v,temp,q,w,du,dv,dtemp,dq,dpsurf)
    !                                                                                      ^----- outputs -----^

    ! Wind increment: specific for 1D
    ! -------------------------------
    ! The physics computes the tendencies on u and v, here we just add Coriolos effect
    !do ilayer = 1,nlayer
    !    du(ilayer) = du(ilayer) + ptif*(v(ilayer) - grv)
    !    dv(ilayer) = dv(ilayer) + ptif*(-u(ilayer) + gru)
    !enddo
    ! For some tests: No coriolis force at equator
    !if (latitude(1) == 0.) then
    du = du + (gru - u)/1.e4
    dv = dv + (grv - v)/1.e4
    !endif

    ! Water tracer increment: specific for 1D
    ! ---------------------------------------
    ! The physics computes the tendency on q for igcm_h2o_vap/igcm_h2o_ice, here we can mimic an effect
    ! of the "3D dynamics" either by forcing a profile or by relaxing towards a prescribed profile
    if (water) then
        if (ctrl_h2ovap) then ! If atmospheric water vapour profile is prescribed
            if (relaxtime_h2ovap < 0.) then ! Forcing
                ! For some tests: unless it reaches saturation (maximal value)
                !call watersat(nlayer,temp,play,zqsat)
                !dq(1,:,igcm_h2o_vap) = (min(zqsat(:),qref_h2ovap(:)) - q(1,:,igcm_h2o_vap))/dttestphys
                dq(1,:,igcm_h2o_vap) = (qref_h2ovap(:) - q(1,:,igcm_h2o_vap))/dttestphys
            else ! Relaxation
                dq(1,:,igcm_h2o_vap) = dq(1,:,igcm_h2o_vap) - (q(1,:,igcm_h2o_vap) - qref_h2ovap(:))/relaxtime_h2ovap
            endif
        endif
        if (ctrl_h2oice) then ! If atmospheric water ice profile is prescribed
            if (relaxtime_h2oice < 0.) then ! Forcing
                dq(1,:,igcm_h2o_ice) = (qref_h2oice(:) - q(1,:,igcm_h2o_ice))/dttestphys
            else ! Relaxation
                dq(1,:,igcm_h2o_ice) = dq(1,:,igcm_h2o_ice) - (q(1,:,igcm_h2o_ice) - qref_h2oice(:))/relaxtime_h2oice
            endif
        endif
     endif

    ! Compute time for next time step
    ! -------------------------------
    firstcall = .false.
    time = time + dttestphys/daysec
    if (time > 1.) then
        time = time - 1.
        day = day + 1
    endif

    ! Compute winds for next time step
    ! --------------------------------
    u = u + dttestphys*du
    v = v + dttestphys*dv

    ! Compute temperature for next time step
    ! --------------------------------------
    temp = temp + dttestphys*dtemp

    ! Compute pressure for next time step
    ! -----------------------------------
    psurf = psurf + dttestphys*dpsurf(1) ! Surface pressure change
    plev = ap + psurf*bp
    play = aps + psurf*bps

    ! Compute tracers for next time step
    ! ----------------------------------
    q = q + dttestphys*dq
enddo ! End of time stepping loop (idt=1,ndt)

! Writing the "restart1D.txt" file for the next run
if (startfiles_1D) call writerestart1D('restart1D.txt',psurf,pa,preff,tsurf(1,:),nlayer,size(tsurf,2),temp,u,v,nq,noms,q)

write(*,*) "testphys1d: everything is cool!"

END PROGRAM testphys1d

!***********************************************************************
!***********************************************************************
! Dummy subroutine used only in 3D, but required to
! compile testphys1d (to cleanly use writediagfi)
SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)

implicit none

integer                       :: im, jm, ngrid, nfield
real, dimension(im,jm,nfield) :: pdyn
real, dimension(ngrid,nfield) :: pfi

if (ngrid /= 1) error stop 'gr_fi_dyn error: in 1D ngrid should be 1!!!'

pdyn(1,1,1:nfield) = pfi(1,1:nfield)

END SUBROUTINE gr_fi_dyn

!***********************************************************************
!***********************************************************************