| 1 | ! $Id: physiq.F 1565 2011-08-31 12:53:29Z jghattas $ |
|---|
| 2 | MODULE physiq_mod |
|---|
| 3 | |
|---|
| 4 | IMPLICIT NONE |
|---|
| 5 | |
|---|
| 6 | CONTAINS |
|---|
| 7 | |
|---|
| 8 | SUBROUTINE physiq(nlon, nlev, & |
|---|
| 9 | debut, lafin, pdtphys, & |
|---|
| 10 | paprs, pplay, pphi, pphis, presnivs, & |
|---|
| 11 | u, v, t, qx, & |
|---|
| 12 | flxmass_w, & |
|---|
| 13 | d_u, d_v, d_t, d_qx, d_ps) |
|---|
| 14 | |
|---|
| 15 | USE dimphy, ONLY: klon, klev |
|---|
| 16 | USE infotrac_phy, ONLY: nqtot |
|---|
| 17 | USE lmdz_geometry, ONLY: latitude |
|---|
| 18 | USE comcstphy, ONLY: rg |
|---|
| 19 | USE iophy, ONLY: histbeg_phy, histwrite_phy |
|---|
| 20 | USE ioipsl, ONLY: getin, histvert, histdef, histend, ymds2ju |
|---|
| 21 | USE lmdz_phys_para, ONLY: jj_nb |
|---|
| 22 | USE phys_state_var_mod, ONLY: phys_state_var_init |
|---|
| 23 | USE lmdz_grid_phy, ONLY: nbp_lon, nbp_lat |
|---|
| 24 | |
|---|
| 25 | USE lmdz_xios, ONLY: xios_update_calendar, using_xios |
|---|
| 26 | USE lmdz_wxios, ONLY: wxios_add_vaxis, wxios_set_cal, wxios_closedef |
|---|
| 27 | USE iophy, ONLY: histwrite_phy |
|---|
| 28 | |
|---|
| 29 | IMPLICIT NONE |
|---|
| 30 | |
|---|
| 31 | ! Routine argument: |
|---|
| 32 | |
|---|
| 33 | INTEGER, INTENT(IN) :: nlon ! number of atmospheric colums |
|---|
| 34 | INTEGER, INTENT(IN) :: nlev ! number of vertical levels (should be =klev) |
|---|
| 35 | LOGICAL, INTENT(IN) :: debut ! signals first CALL to physics |
|---|
| 36 | LOGICAL, INTENT(IN) :: lafin ! signals last CALL to physics |
|---|
| 37 | REAL, INTENT(IN) :: pdtphys ! physics time step (s) |
|---|
| 38 | REAL, INTENT(IN) :: paprs(klon, klev + 1) ! interlayer pressure (Pa) |
|---|
| 39 | REAL, INTENT(IN) :: pplay(klon, klev) ! mid-layer pressure (Pa) |
|---|
| 40 | REAL, INTENT(IN) :: pphi(klon, klev) ! geopotential at mid-layer |
|---|
| 41 | REAL, INTENT(IN) :: pphis(klon) ! surface geopotential |
|---|
| 42 | REAL, INTENT(IN) :: presnivs(klev) ! pseudo-pressure (Pa) of mid-layers |
|---|
| 43 | REAL, INTENT(IN) :: u(klon, klev) ! eastward zonal wind (m/s) |
|---|
| 44 | REAL, INTENT(IN) :: v(klon, klev) ! northward meridional wind (m/s) |
|---|
| 45 | REAL, INTENT(IN) :: t(klon, klev) ! temperature (K) |
|---|
| 46 | REAL, INTENT(IN) :: qx(klon, klev, nqtot) ! tracers (.../kg_air) |
|---|
| 47 | REAL, INTENT(IN) :: flxmass_w(klon, klev) ! vertical mass flux |
|---|
| 48 | REAL, INTENT(OUT) :: d_u(klon, klev) ! physics tendency on u (m/s/s) |
|---|
| 49 | REAL, INTENT(OUT) :: d_v(klon, klev) ! physics tendency on v (m/s/s) |
|---|
| 50 | REAL, INTENT(OUT) :: d_t(klon, klev) ! physics tendency on t (K/s) |
|---|
| 51 | REAL, INTENT(OUT) :: d_qx(klon, klev, nqtot) ! physics tendency on tracers |
|---|
| 52 | REAL, INTENT(OUT) :: d_ps(klon) ! physics tendency on surface pressure |
|---|
| 53 | |
|---|
| 54 | INTEGER, save :: itau = 0 ! counter to count number of calls to physics |
|---|
| 55 | !$OMP THREADPRIVATE(itau) |
|---|
| 56 | REAL :: temp_newton(klon, klev) |
|---|
| 57 | INTEGER :: k |
|---|
| 58 | logical, save :: first = .TRUE. |
|---|
| 59 | !$OMP THREADPRIVATE(first) |
|---|
| 60 | |
|---|
| 61 | ! For I/Os |
|---|
| 62 | INTEGER :: itau0 |
|---|
| 63 | REAL :: zjulian |
|---|
| 64 | REAL :: dtime |
|---|
| 65 | INTEGER :: nhori ! horizontal coordinate ID |
|---|
| 66 | INTEGER, save :: nid_hist ! output file ID |
|---|
| 67 | !$OMP THREADPRIVATE(nid_hist) |
|---|
| 68 | INTEGER :: zvertid ! vertical coordinate ID |
|---|
| 69 | INTEGER, save :: iwrite_phys ! output every iwrite_phys physics step |
|---|
| 70 | !$OMP THREADPRIVATE(iwrite_phys) |
|---|
| 71 | INTEGER, save :: iwrite_phys_omp ! intermediate variable to read iwrite_phys |
|---|
| 72 | ! (must be shared by all threads) |
|---|
| 73 | REAL :: t_ops ! frequency of the IOIPSL operations (eg average over...) |
|---|
| 74 | REAL :: t_wrt ! frequency of the IOIPSL outputs |
|---|
| 75 | |
|---|
| 76 | ! initializations |
|---|
| 77 | IF (debut) then ! Things to do only for the first CALL to physics |
|---|
| 78 | ! load initial conditions for physics (including the grid) |
|---|
| 79 | CALL phys_state_var_init() ! some initializations, required before calling phyetat0 |
|---|
| 80 | CALL phyetat0("startphy.nc") |
|---|
| 81 | |
|---|
| 82 | ! Initialize outputs: |
|---|
| 83 | itau0 = 0 |
|---|
| 84 | !$OMP MASTER |
|---|
| 85 | iwrite_phys_omp = 1 !default: output every physics timestep |
|---|
| 86 | ! NB: getin() is not threadsafe; only one thread should CALL it. |
|---|
| 87 | CALL getin("iwrite_phys", iwrite_phys_omp) |
|---|
| 88 | !$OMP END MASTER |
|---|
| 89 | !$OMP BARRIER |
|---|
| 90 | iwrite_phys = iwrite_phys_omp |
|---|
| 91 | t_ops = pdtphys * iwrite_phys ! frequency of the IOIPSL operation |
|---|
| 92 | t_wrt = pdtphys * iwrite_phys ! frequency of the outputs in the file |
|---|
| 93 | ! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0 |
|---|
| 94 | !CALL ymds2ju(annee0, month, dayref, hour, zjulian) |
|---|
| 95 | CALL ymds2ju(1979, 1, 1, 0.0, zjulian) |
|---|
| 96 | dtime = pdtphys |
|---|
| 97 | #ifndef CPP_IOIPSL_NO_OUTPUT |
|---|
| 98 | ! Initialize IOIPSL output file |
|---|
| 99 | CALL histbeg_phy("histins.nc", itau0, zjulian, dtime, nhori, nid_hist) |
|---|
| 100 | #endif |
|---|
| 101 | |
|---|
| 102 | !$OMP MASTER |
|---|
| 103 | |
|---|
| 104 | #ifndef CPP_IOIPSL_NO_OUTPUT |
|---|
| 105 | ! IOIPSL |
|---|
| 106 | ! define vertical coordinate |
|---|
| 107 | CALL histvert(nid_hist, "presnivs", "Vertical levels", "Pa", klev, & |
|---|
| 108 | presnivs, zvertid, 'down') |
|---|
| 109 | ! define variables which will be written in "histins.nc" file |
|---|
| 110 | CALL histdef(nid_hist, 'temperature', 'Atmospheric temperature', 'K', & |
|---|
| 111 | nbp_lon, jj_nb, nhori, klev, 1, klev, zvertid, 32, & |
|---|
| 112 | 'inst(X)', t_ops, t_wrt) |
|---|
| 113 | CALL histdef(nid_hist, 'u', 'Eastward Zonal Wind', 'm/s', & |
|---|
| 114 | nbp_lon, jj_nb, nhori, klev, 1, klev, zvertid, 32, & |
|---|
| 115 | 'inst(X)', t_ops, t_wrt) |
|---|
| 116 | CALL histdef(nid_hist, 'v', 'Northward Meridional Wind', 'm/s', & |
|---|
| 117 | nbp_lon, jj_nb, nhori, klev, 1, klev, zvertid, 32, & |
|---|
| 118 | 'inst(X)', t_ops, t_wrt) |
|---|
| 119 | CALL histdef(nid_hist, 'ps', 'Surface Pressure', 'Pa', & |
|---|
| 120 | nbp_lon, jj_nb, nhori, 1, 1, 1, zvertid, 32, & |
|---|
| 121 | 'inst(X)', t_ops, t_wrt) |
|---|
| 122 | ! end definition sequence |
|---|
| 123 | CALL histend(nid_hist) |
|---|
| 124 | #endif |
|---|
| 125 | |
|---|
| 126 | IF (using_xios) THEN |
|---|
| 127 | !XIOS |
|---|
| 128 | ! Declare available vertical axes to be used in output files: |
|---|
| 129 | CALL wxios_add_vaxis("presnivs", klev, presnivs) |
|---|
| 130 | |
|---|
| 131 | ! Declare calendar and time step |
|---|
| 132 | CALL wxios_set_cal(dtime, "earth_360d", 1, 1, 1, 0.0, 1, 1, 1, 0.0) |
|---|
| 133 | |
|---|
| 134 | !Finalize the context: |
|---|
| 135 | CALL wxios_closedef() |
|---|
| 136 | ENDIF |
|---|
| 137 | !$OMP END MASTER |
|---|
| 138 | !$OMP BARRIER |
|---|
| 139 | END IF ! of if (debut) |
|---|
| 140 | |
|---|
| 141 | ! increment local time counter itau |
|---|
| 142 | itau = itau + 1 |
|---|
| 143 | |
|---|
| 144 | ! set all tendencies to zero |
|---|
| 145 | d_u(1:klon, 1:klev) = 0. |
|---|
| 146 | d_v(1:klon, 1:klev) = 0. |
|---|
| 147 | d_t(1:klon, 1:klev) = 0. |
|---|
| 148 | d_qx(1:klon, 1:klev, 1:nqtot) = 0. |
|---|
| 149 | d_ps(1:klon) = 0. |
|---|
| 150 | |
|---|
| 151 | ! compute tendencies to return to the dynamics: |
|---|
| 152 | ! "friction" on the first layer |
|---|
| 153 | d_u(1:klon, 1) = -u(1:klon, 1) / 86400. |
|---|
| 154 | d_v(1:klon, 1) = -v(1:klon, 1) / 86400. |
|---|
| 155 | ! newtonian relaxation towards temp_newton() |
|---|
| 156 | DO k = 1, klev |
|---|
| 157 | temp_newton(1:klon, k) = 280. + cos(latitude(1:klon)) * 40. - pphi(1:klon, k) / rg * 6.e-3 |
|---|
| 158 | d_t(1:klon, k) = (temp_newton(1:klon, k) - t(1:klon, k)) / 1.e5 |
|---|
| 159 | END DO |
|---|
| 160 | |
|---|
| 161 | PRINT*, 'PHYDEV: itau=', itau |
|---|
| 162 | |
|---|
| 163 | ! write some outputs: |
|---|
| 164 | ! IOIPSL |
|---|
| 165 | #ifndef CPP_IOIPSL_NO_OUTPUT |
|---|
| 166 | IF (modulo(itau, iwrite_phys)==0) THEN |
|---|
| 167 | CALL histwrite_phy(nid_hist, .FALSE., "temperature", itau, t) |
|---|
| 168 | CALL histwrite_phy(nid_hist, .FALSE., "u", itau, u) |
|---|
| 169 | CALL histwrite_phy(nid_hist, .FALSE., "v", itau, v) |
|---|
| 170 | CALL histwrite_phy(nid_hist, .FALSE., "ps", itau, paprs(:, 1)) |
|---|
| 171 | END IF |
|---|
| 172 | #endif |
|---|
| 173 | |
|---|
| 174 | !XIOS |
|---|
| 175 | IF (using_xios) THEN |
|---|
| 176 | !$OMP MASTER |
|---|
| 177 | !Increment XIOS time |
|---|
| 178 | CALL xios_update_calendar(itau) |
|---|
| 179 | !$OMP END MASTER |
|---|
| 180 | !$OMP BARRIER |
|---|
| 181 | |
|---|
| 182 | !Send fields to XIOS: (NB these fields must also be defined as |
|---|
| 183 | ! <field id="..." /> in iodef.xml to be correctly used |
|---|
| 184 | CALL histwrite_phy("temperature", t) |
|---|
| 185 | CALL histwrite_phy("temp_newton", temp_newton) |
|---|
| 186 | CALL histwrite_phy("u", u) |
|---|
| 187 | CALL histwrite_phy("v", v) |
|---|
| 188 | CALL histwrite_phy("ps", paprs(:, 1)) |
|---|
| 189 | ENDIF |
|---|
| 190 | |
|---|
| 191 | ! if lastcall, then it is time to write "restartphy.nc" file |
|---|
| 192 | IF (lafin) THEN |
|---|
| 193 | CALL phyredem("restartphy.nc") |
|---|
| 194 | END IF |
|---|
| 195 | |
|---|
| 196 | END SUBROUTINE physiq |
|---|
| 197 | |
|---|
| 198 | END MODULE physiq_mod |
|---|
