! $Id: physiq.F 1565 2011-08-31 12:53:29Z jghattas $
MODULE physiq_mod

IMPLICIT NONE

CONTAINS

      SUBROUTINE physiq (nlon,nlev, &
     &            debut,lafin,pdtphys, &
     &            paprs,pplay,pphi,pphis,presnivs, &
     &            u,v,t,qx, &
     &            flxmass_w, &
     &            d_u, d_v, d_t, d_qx, d_ps)

      USE dimphy, only : klon,klev
      USE infotrac_phy, only : nqtot
      USE geometry_mod, only : latitude
      USE comcstphy, only : rg
      USE iophy, only : histbeg_phy,histwrite_phy
      USE ioipsl, only : getin,histvert,histdef,histend,ymds2ju
      USE mod_phys_lmdz_para, only : jj_nb
      USE phys_state_var_mod, only : phys_state_var_init
      USE mod_grid_phy_lmdz, ONLY: nbp_lon,nbp_lat

      USE lmdz_xios, ONLY: xios_update_calendar, using_xios
      use wxios_mod, only: wxios_add_vaxis, wxios_set_cal, wxios_closedef
      USE iophy, ONLY: histwrite_phy

      IMPLICIT none
!
! Routine argument:
!
      integer,intent(in) :: nlon ! number of atmospheric colums
      integer,intent(in) :: nlev ! number of vertical levels (should be =klev)
      logical,intent(in) :: debut ! signals first call to physics
      logical,intent(in) :: lafin ! signals last call to physics
      real,intent(in) :: pdtphys ! physics time step (s)
      real,intent(in) :: paprs(klon,klev+1) ! interlayer pressure (Pa)
      real,intent(in) :: pplay(klon,klev) ! mid-layer pressure (Pa)
      real,intent(in) :: pphi(klon,klev) ! geopotential at mid-layer
      real,intent(in) :: pphis(klon) ! surface geopotential
      real,intent(in) :: presnivs(klev) ! pseudo-pressure (Pa) of mid-layers
      real,intent(in) :: u(klon,klev) ! eastward zonal wind (m/s)
      real,intent(in) :: v(klon,klev) ! northward meridional wind (m/s)
      real,intent(in) :: t(klon,klev) ! temperature (K)
      real,intent(in) :: qx(klon,klev,nqtot) ! tracers (.../kg_air)
      real,intent(in) :: flxmass_w(klon,klev) ! vertical mass flux
      real,intent(out) :: d_u(klon,klev) ! physics tendency on u (m/s/s)
      real,intent(out) :: d_v(klon,klev) ! physics tendency on v (m/s/s)
      real,intent(out) :: d_t(klon,klev) ! physics tendency on t (K/s)
      real,intent(out) :: d_qx(klon,klev,nqtot) ! physics tendency on tracers
      real,intent(out) :: d_ps(klon) ! physics tendency on surface pressure

integer,save :: itau=0 ! counter to count number of calls to physics
!$OMP THREADPRIVATE(itau)
real :: temp_newton(klon,klev)
integer :: k
logical, save :: first=.true.
!$OMP THREADPRIVATE(first)

! For I/Os
integer :: itau0
real :: zjulian
real :: dtime
integer :: nhori ! horizontal coordinate ID
integer,save :: nid_hist ! output file ID
!$OMP THREADPRIVATE(nid_hist)
integer :: zvertid ! vertical coordinate ID
integer,save :: iwrite_phys ! output every iwrite_phys physics step
!$OMP THREADPRIVATE(iwrite_phys)
integer,save :: iwrite_phys_omp ! intermediate variable to read iwrite_phys
                                ! (must be shared by all threads)
real :: t_ops ! frequency of the IOIPSL operations (eg average over...)
real :: t_wrt ! frequency of the IOIPSL outputs

! initializations
if (debut) then ! Things to do only for the first call to physics 
! load initial conditions for physics (including the grid)
  call phys_state_var_init() ! some initializations, required before calling phyetat0
  call phyetat0("startphy.nc")

! Initialize outputs:
  itau0=0
!$OMP MASTER
  iwrite_phys_omp=1 !default: output every physics timestep
  ! NB: getin() is not threadsafe; only one thread should call it.
  call getin("iwrite_phys",iwrite_phys_omp)
!$OMP END MASTER
!$OMP BARRIER
  iwrite_phys=iwrite_phys_omp
  t_ops=pdtphys*iwrite_phys ! frequency of the IOIPSL operation
  t_wrt=pdtphys*iwrite_phys ! frequency of the outputs in the file
  ! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0
  !CALL ymds2ju(annee0, month, dayref, hour, zjulian)
  call ymds2ju(1979, 1, 1, 0.0, zjulian)
  dtime=pdtphys
#ifndef CPP_IOIPSL_NO_OUTPUT
  ! Initialize IOIPSL output file
  call histbeg_phy("histins.nc",itau0,zjulian,dtime,nhori,nid_hist)
#endif

!$OMP MASTER

#ifndef CPP_IOIPSL_NO_OUTPUT 
! IOIPSL
  ! define vertical coordinate
  call histvert(nid_hist,"presnivs","Vertical levels","Pa",klev, &
                presnivs,zvertid,'down')
  ! define variables which will be written in "histins.nc" file
  call histdef(nid_hist,'temperature','Atmospheric temperature','K', &
               nbp_lon,jj_nb,nhori,klev,1,klev,zvertid,32, &
               'inst(X)',t_ops,t_wrt)
  call histdef(nid_hist,'u','Eastward Zonal Wind','m/s', &
               nbp_lon,jj_nb,nhori,klev,1,klev,zvertid,32, &
               'inst(X)',t_ops,t_wrt)
  call histdef(nid_hist,'v','Northward Meridional Wind','m/s', &
               nbp_lon,jj_nb,nhori,klev,1,klev,zvertid,32, &
               'inst(X)',t_ops,t_wrt)
  call histdef(nid_hist,'ps','Surface Pressure','Pa', &
               nbp_lon,jj_nb,nhori,1,1,1,zvertid,32, &
               'inst(X)',t_ops,t_wrt)
  ! end definition sequence
  call histend(nid_hist)
#endif

    IF (using_xios) THEN
!XIOS
      ! Declare available vertical axes to be used in output files:    
      CALL wxios_add_vaxis("presnivs", klev, presnivs)

      ! Declare calendar and time step
      CALL wxios_set_cal(dtime,"earth_360d",1,1,1,0.0,1,1,1,0.0)
    
      !Finalize the context:
      CALL wxios_closedef()
    ENDIF
!$OMP END MASTER
!$OMP BARRIER
endif ! of if (debut)

! increment local time counter itau
itau=itau+1

! set all tendencies to zero
d_u(1:klon,1:klev)=0.
d_v(1:klon,1:klev)=0.
d_t(1:klon,1:klev)=0.
d_qx(1:klon,1:klev,1:nqtot)=0.
d_ps(1:klon)=0.

! compute tendencies to return to the dynamics:
! "friction" on the first layer
d_u(1:klon,1)=-u(1:klon,1)/86400.
d_v(1:klon,1)=-v(1:klon,1)/86400.
! newtonian relaxation towards temp_newton()
do k=1,klev
  temp_newton(1:klon,k)=280.+cos(latitude(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3
  d_t(1:klon,k)=(temp_newton(1:klon,k)-t(1:klon,k))/1.e5
enddo


print*,'PHYDEV: itau=',itau

! write some outputs:
! IOIPSL
#ifndef CPP_IOIPSL_NO_OUTPUT 
if (modulo(itau,iwrite_phys)==0) then
  call histwrite_phy(nid_hist,.false.,"temperature",itau,t)
  call histwrite_phy(nid_hist,.false.,"u",itau,u)
  call histwrite_phy(nid_hist,.false.,"v",itau,v)
  call histwrite_phy(nid_hist,.false.,"ps",itau,paprs(:,1))
endif
#endif

!XIOS
IF (using_xios) THEN
  !$OMP MASTER
    !Increment XIOS time
    CALL xios_update_calendar(itau)
  !$OMP END MASTER
  !$OMP BARRIER

    !Send fields to XIOS: (NB these fields must also be defined as
    ! <field id="..." /> in iodef.xml to be correctly used
    CALL histwrite_phy("temperature",t)
    CALL histwrite_phy("temp_newton",temp_newton)
    CALL histwrite_phy("u",u)
    CALL histwrite_phy("v",v)
    CALL histwrite_phy("ps",paprs(:,1))
ENDIF

! if lastcall, then it is time to write "restartphy.nc" file
if (lafin) then
  call phyredem("restartphy.nc")
endif

end subroutine physiq

END MODULE physiq_mod