source: LMDZ6/branches/Amaury_dev/libf/phydev/physiq_mod.F90 @ 5186

! $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 lmdz_geometry, ONLY: latitude
    USE comcstphy, ONLY: rg
    USE iophy, ONLY: histbeg_phy, histwrite_phy
    USE ioipsl, ONLY: getin, histvert, histdef, histend, ymds2ju
    USE lmdz_phys_para, ONLY: jj_nb
    USE phys_state_var_mod, ONLY: phys_state_var_init
    USE lmdz_grid_phy, ONLY: nbp_lon, nbp_lat

    USE lmdz_xios, ONLY: xios_update_calendar, using_xios
    USE lmdz_wxios, 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
    END IF ! 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
    END DO

    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))
    END IF
#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")
    END IF

  END SUBROUTINE  physiq

END MODULE physiq_mod