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