source: trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90 @ 3067

PROGRAM testphys1d

use comsoil_h,           only: inertiedat, inertiesoil, nsoilmx, tsoil
use surfdat_h,           only: albedodat, perenial_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, rad, omeg, g, mugaz, rcp, r, cpp
use time_phylmdz_mod,    only: daysec, day_step
use dimradmars_mod,      only: tauvis, totcloudfrac, albedo
use dust_param_mod,      only: tauscaling
use comvert_mod,         only: ap, bp, aps, bps, pa, preff, sig
use physiq_mod,          only: physiq
use turb_mod,            only: q2
use write_output_mod,    only: write_output
use init_testphys1d_mod, only: init_testphys1d
! 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"
!#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                             :: unitstart     ! unite d'ecriture de "startfi"
integer                             :: ndt, ilayer, ilevel, isoil, idt, iq
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 tandencies (e.g. m.s-2, K/s, Pa/s)
real, dimension(nlayer)             :: du, dv, dtemp, dudyn, dvdyn, dtempdyn
real, dimension(1)                  :: dpsurf
real, dimension(:,:,:), allocatable :: dq, dqdyn

! Various intermediate variables
integer                 :: aslun
real                    :: zls, pks, ptif, qtotinit, qtot
real, dimension(nlayer) :: phi, h, s, w
integer                 :: nq = 1 ! number of tracers
real, dimension(1)      :: latitude, longitude, cell_area
logical                 :: there
character(len = 2)      :: str2
character(len = 7)      :: str7
character(len = 44)     :: txt

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

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

!=======================================================================
! INITIALISATION
!=======================================================================
#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

call init_testphys1d(.false.,ngrid,nlayer,odpref,nq,q,time,psurf,u,v,temp,startfiles_1D,therestart1D, &
                     therestartfi,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,dttestphys,time,       &
                  tsurf,tsoil,inertiesoil,albedo,emis,q2,qsurf,tauscaling, &
                  totcloudfrac,wstar,watercap,perenial_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_vap) = min(zqsat(:),q(1,:,igcm_h2o_vap))
                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(1,:,:) = q(1,:,:) + dttestphys*dq(1,:,:)
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,temp,u,v,nq,noms,qsurf(1,:,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) then
    write(*,*) "gr_fi_dyn error: in 1D ngrid should be 1!!!"
    stop
endif

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

END SUBROUTINE gr_fi_dyn

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