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

IMPLICIT NONE

CONTAINS

      SUBROUTINE physiqex(nlon,nlev, &
              debut,lafin,pdtphys, &
              paprs,pplay,pphi,pphis,presnivs, &
              u,v,rot,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 lmdz_geometry, ONLY: latitude
!      USE comcstphy, ONLY: rg
      USE ioipsl, ONLY: ymds2ju
      USE phys_state_var_mod, ONLY: phys_state_var_init
      USE phyetat0_mod, ONLY: phyetat0
      USE output_physiqex_mod, ONLY: output_physiqex

      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) :: rot(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

!    include "clesphys.h"
    INTEGER        length
    PARAMETER    ( length = 100 )
    REAL tabcntr0( length       )
    INTEGER, PARAMETER :: longcles=20
    REAL, SAVE :: clesphy0(longcles)
    !$OMP THREADPRIVATE(clesphy0)


REAL :: temp_newton(klon,klev)
INTEGER :: k
logical, save :: first=.TRUE.
!$OMP THREADPRIVATE(first)

REAL,save :: rg=9.81
!$OMP THREADPRIVATE(rg)

! For I/Os
INTEGER :: itau0
REAL :: zjulian


!------------------------------------------------------------
! Initialisations de la physique au premier pas de temps
!------------------------------------------------------------

PRINT*,'Debut physiqex',debut
! initializations
IF (debut) then ! Things to do only for the first CALL to physics
PRINT*,'Debut physiqex IN'

! load initial conditions for physics (including the grid)
  CALL phys_state_var_init(1) ! some initializations, required before calling phyetat0
  CALL phyetat0("startphy.nc", clesphy0, tabcntr0)

! Initialize outputs:
  itau0=0
  ! 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)

#ifndef CPP_IOIPSL_NO_OUTPUT
  ! Initialize IOIPSL output file
#endif

END IF ! of if (debut)

!------------------------------------------------------------
! Initialisations a chaque pas de temps
!------------------------------------------------------------


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

!------------------------------------------------------------
! Calculs
!------------------------------------------------------------

! 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
END DO


!------------------------------------------------------------
! Entrees sorties
!------------------------------------------------------------


CALL output_physiqex(debut,zjulian,pdtphys,presnivs,paprs,u,v,t,qx,0.*t,0.*t,0.*t,0.*t,0.*t,0.*t)


! if lastcall, then it is time to write "restartphy.nc" file
IF (lafin) THEN
  CALL phyredem("restartphy.nc")
END IF


END SUBROUTINE  physiqex

END MODULE physiqex_mod