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, igcm_co2, 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
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
! Mostly for XIOS outputs:
use mod_const_mpi,       only: init_const_mpi
use parallel_lmdz,       only: init_parallel
use version_info_mod,    only: print_version_info

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"
!#include "dimradmars.h"
!#include "comgeomfi.h"
!#include "surfdat.h"
!#include "slope.h"
!#include "comsoil.h"
!#include "comdiurn.h"
include "callkeys.h"
!#include "comsaison.h"
!#include "control.h"
include "netcdf.inc"
!#include "advtrac.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
character(100)                      :: arg           ! To read command-line arguments

! 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

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

! JN & JBC: Force atmospheric water profiles
real                            :: atm_wat_profile, atm_wat_tau
real, dimension(:), allocatable :: zqsat ! useful for (atm_wat_profile=2)
!=======================================================================

!=======================================================================
! INITIALISATION
!=======================================================================
if (command_argument_count() > 0) then ! Get the number of command-line arguments
    call get_command_argument(1,arg) ! Read the argument given to the program
    select case (trim(adjustl(arg)))
        case('version')
            call print_version_info()
            stop
        case default
            error stop 'The argument given to the program is unknown!'
    end select
endif

#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 (startfiles_1D .and. .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

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,atm_wat_profile,atm_wat_tau)

! 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

    ! Modify atmospheric water if asked
    ! Added "atm_wat_profile" flag (JN + JBC)
    ! ---------------------------------------
    if (water) then
        call watersat(nlayer,temp,play,zqsat)
        if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.) then
        ! If atmospheric water is monitored
            if (atm_wat_tau < 0.) then ! Prescribed atm_wat_profile: wet if >0, dry if =0
                q(1,:,igcm_h2o_vap) = min(zqsat,atm_wat_profile*g/psurf)
                q(1,:,igcm_h2o_ice) = 0. ! reset h2o ice
            else ! Relaxation towards the value atm_wat_profile with relaxation time atm_wat_tau
                q(1,:,igcm_h2o_vap) = atm_wat_profile*g/psurf + (q(1,:,igcm_h2o_vap) - atm_wat_profile*g/psurf)*dexp(-dttestphys/atm_wat_tau)
                q(1,:,igcm_h2o_ice) = 0. ! reset h2o ice
            endif
        endif
     endif

    ! 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 compute 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

    ! 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 and temperature for next time step
    ! ------------------------------------------------
    u = u + dttestphys*du
    v = v + dttestphys*dv
    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

    ! Increment tracers
    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,tsurf(1,:),nlayer,size(tsurf,2),temp,u,v,nq,noms,qsurf(1,:,:),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

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