| 1 | MODULE physics_dcmip_mod |
|---|
| 2 | USE ICOSA |
|---|
| 3 | PRIVATE |
|---|
| 4 | |
|---|
| 5 | INTEGER,SAVE :: testcase |
|---|
| 6 | !$OMP THREADPRIVATE(testcase) |
|---|
| 7 | |
|---|
| 8 | TYPE(t_field),POINTER :: f_out_i(:) |
|---|
| 9 | TYPE(t_field),POINTER :: f_precl(:) |
|---|
| 10 | REAL(rstd),POINTER :: out_i(:,:) |
|---|
| 11 | |
|---|
| 12 | PUBLIC :: compute_phys_wrap, init_physics |
|---|
| 13 | |
|---|
| 14 | CONTAINS |
|---|
| 15 | |
|---|
| 16 | SUBROUTINE init_physics |
|---|
| 17 | IMPLICIT NONE |
|---|
| 18 | testcase=1 ! OK for 4.2 (moist baroclinic instability) |
|---|
| 19 | CALL getin("dcmip_physics",testcase) |
|---|
| 20 | END SUBROUTINE init_physics |
|---|
| 21 | |
|---|
| 22 | SUBROUTINE compute_phys_wrap(args) |
|---|
| 23 | USE physics_interface_mod |
|---|
| 24 | TYPE(physics_inout) :: args |
|---|
| 25 | CALL compute_physics(args%ngrid, args%dt_phys, args%lat, & |
|---|
| 26 | args%p, args%Temp, args%ulon, args%ulat, args%q(:,:,1), & |
|---|
| 27 | args%dTemp, args%dulon, args%dulat, args%dq(:,:,1), args%extra_2D(:,1)) |
|---|
| 28 | END SUBROUTINE compute_phys_wrap |
|---|
| 29 | |
|---|
| 30 | SUBROUTINE compute_physics(ngrid,dt_phys,lat, p,Temp,u,v,q, dTemp,du,dv,dq, precl) |
|---|
| 31 | USE icosa |
|---|
| 32 | IMPLICIT NONE |
|---|
| 33 | INTEGER :: ngrid |
|---|
| 34 | REAL(rstd) :: lat(ngrid) |
|---|
| 35 | REAL(rstd) :: ps(ngrid) |
|---|
| 36 | REAL(rstd) :: precl(ngrid) |
|---|
| 37 | ! arguments with bottom-up indexing (DYNAMICO) |
|---|
| 38 | REAL(rstd) :: p(ngrid,llm+1) |
|---|
| 39 | REAL(rstd) :: Temp(ngrid,llm) |
|---|
| 40 | REAL(rstd) :: u(ngrid,llm) |
|---|
| 41 | REAL(rstd) :: v(ngrid,llm) |
|---|
| 42 | REAL(rstd) :: q(ngrid,llm) |
|---|
| 43 | REAL(rstd) :: dTemp(ngrid,llm) |
|---|
| 44 | REAL(rstd) :: du(ngrid,llm) |
|---|
| 45 | REAL(rstd) :: dv(ngrid,llm) |
|---|
| 46 | REAL(rstd) :: dq(ngrid,llm) |
|---|
| 47 | ! local arrays with top-down vertical indexing (DCMIP) |
|---|
| 48 | REAL(rstd) :: pint(ngrid,llm+1) |
|---|
| 49 | REAL(rstd) :: pmid(ngrid,llm) |
|---|
| 50 | REAL(rstd) :: pdel(ngrid,llm) |
|---|
| 51 | REAL(rstd) :: Tfi(ngrid,llm) |
|---|
| 52 | REAL(rstd) :: ufi(ngrid,llm) |
|---|
| 53 | REAL(rstd) :: vfi(ngrid,llm) |
|---|
| 54 | REAL(rstd) :: qfi(ngrid,llm) |
|---|
| 55 | |
|---|
| 56 | INTEGER :: i,j,l,ll,ij |
|---|
| 57 | REAL(rstd) :: dt_phys, inv_dt |
|---|
| 58 | |
|---|
| 59 | ! prepare input fields and mirror vertical index |
|---|
| 60 | ps(:) = p(:,1) ! surface pressure |
|---|
| 61 | |
|---|
| 62 | DO l=1,llm+1 |
|---|
| 63 | DO j=jj_begin,jj_end |
|---|
| 64 | DO i=ii_begin,ii_end |
|---|
| 65 | ij=(j-1)*iim+i |
|---|
| 66 | pint(ij,l)=p(ij,llm+2-l) |
|---|
| 67 | ENDDO |
|---|
| 68 | ENDDO |
|---|
| 69 | ENDDO |
|---|
| 70 | |
|---|
| 71 | DO l=1,llm |
|---|
| 72 | ll=llm+1-l |
|---|
| 73 | DO j=jj_begin,jj_end |
|---|
| 74 | DO i=ii_begin,ii_end |
|---|
| 75 | ij=(j-1)*iim+i |
|---|
| 76 | pmid(ij,l)=0.5*(pint(ij,l)+pint(ij,l+1)) ! Pressure inside layers |
|---|
| 77 | pdel(ij,l)=pint(ij,l+1)-pint(ij,l) ! Pressure difference between two layers |
|---|
| 78 | ufi(ij,l)=u(ij,ll) |
|---|
| 79 | vfi(ij,l)=v(ij,ll) |
|---|
| 80 | qfi(ij,l)=q(ij,ll) |
|---|
| 81 | Tfi(ij,l)=Temp(ij,ll)/(1+0.608*qfi(ij,l)) |
|---|
| 82 | ENDDO |
|---|
| 83 | ENDDO |
|---|
| 84 | ENDDO |
|---|
| 85 | |
|---|
| 86 | CALL simple_physics(ngrid, llm, dt_phys, lat, tfi, qfi , ufi, vfi, pmid, pint, pdel, 1./pdel, ps, precl, testcase) |
|---|
| 87 | |
|---|
| 88 | ! Mirror vertical index and compute tendencies |
|---|
| 89 | inv_dt = 1./dt_phys |
|---|
| 90 | DO l=1,llm |
|---|
| 91 | ll=llm+1-l |
|---|
| 92 | DO j=jj_begin,jj_end |
|---|
| 93 | DO i=ii_begin,ii_end |
|---|
| 94 | ij=(j-1)*iim+i |
|---|
| 95 | dTemp(ij,l) = inv_dt * ( Tfi(ij,ll)*(1+0.608*qfi(ij,ll)) - Temp(ij,l) ) |
|---|
| 96 | du(ij,l) = inv_dt * (ufi(ij,ll) - u(ij,l)) |
|---|
| 97 | dv(ij,l) = inv_dt * (vfi(ij,ll) - v(ij,l)) |
|---|
| 98 | dq(ij,l) = inv_dt * (qfi(ij,ll) - q(ij,l)) |
|---|
| 99 | ENDDO |
|---|
| 100 | ENDDO |
|---|
| 101 | ENDDO |
|---|
| 102 | |
|---|
| 103 | END SUBROUTINE compute_physics |
|---|
| 104 | |
|---|
| 105 | END MODULE physics_dcmip_mod |
|---|
| 106 | |
|---|
| 107 | subroutine simple_physics (pcols, pver, dtime, lat, t, q, u, v, pmid, pint, pdel, rpdel, ps, precl, test) |
|---|
| 108 | !----------------------------------------------------------------------- |
|---|
| 109 | ! |
|---|
| 110 | ! Purpose: Simple Physics Package |
|---|
| 111 | ! |
|---|
| 112 | ! Author: K. A. Reed (University of Michigan, kareed@umich.edu) |
|---|
| 113 | ! version 5 |
|---|
| 114 | ! July/8/2012 |
|---|
| 115 | ! |
|---|
| 116 | ! Change log: |
|---|
| 117 | ! v2: removal of some NCAR CAM-specific 'use' associations |
|---|
| 118 | ! v3: corrected precl(i) computation, the precipitation rate is now computed via a vertical integral, the previous single-level computation in v2 was a bug |
|---|
| 119 | ! v3: corrected dtdt(i,1) computation, the term '-(i,1)' was missing the temperature variable: '-t(i,1)' |
|---|
| 120 | ! v4: modified and enhanced parameter list to make the routine truly standalone, the number of columns and vertical levels have been added: pcols, pver |
|---|
| 121 | ! v4: 'ncol' has been removed, 'pcols' is used instead |
|---|
| 122 | ! v5: the sea surface temperature (SST) field Tsurf is now an array, the SST now depends on the latitude |
|---|
| 123 | ! v5: addition of the latitude array 'lat' and the flag 'test' in the parameter list |
|---|
| 124 | ! if test = 0: constant SST is used, correct setting for the tropical cyclone test case 5-1 |
|---|
| 125 | ! if test = 1: newly added latitude-dependent SST is used, correct setting for the moist baroclinic wave test with simple-physics (test 4-3) |
|---|
| 126 | ! |
|---|
| 127 | ! Description: Includes large-scale precipitation, surface fluxes and |
|---|
| 128 | ! boundary-leyer mixing. The processes are time-split |
|---|
| 129 | ! in that order. A partially implicit formulation is |
|---|
| 130 | ! used to foster numerical stability. |
|---|
| 131 | ! The routine assumes that the model levels are ordered |
|---|
| 132 | ! in a top-down approach, e.g. level 1 denotes the uppermost |
|---|
| 133 | ! full model level. |
|---|
| 134 | ! |
|---|
| 135 | ! This routine is based on an implementation which was |
|---|
| 136 | ! developed for the NCAR Community Atmosphere Model (CAM). |
|---|
| 137 | ! Adjustments for other models will be necessary. |
|---|
| 138 | ! |
|---|
| 139 | ! The routine provides both updates of the state variables |
|---|
| 140 | ! u, v, T, q (these are local copies of u,v,T,q within this physics |
|---|
| 141 | ! routine) and also collects their time tendencies. |
|---|
| 142 | ! The latter might be used to couple the physics and dynamics |
|---|
| 143 | ! in a process-split way. For a time-split coupling, the final |
|---|
| 144 | ! state should be given to the dynamical core for the next time step. |
|---|
| 145 | ! Test: 0 = Reed and Jablonowski (2011) tropical cyclone test case (test 5-1) |
|---|
| 146 | ! 1 = Moist baroclinic instability test (test 4-3) |
|---|
| 147 | ! |
|---|
| 148 | ! |
|---|
| 149 | ! Reference: Reed, K. A. and C. Jablonowski (2012), Idealized tropical cyclone |
|---|
| 150 | ! simulations of intermediate complexity: A test case for AGCMs, |
|---|
| 151 | ! J. Adv. Model. Earth Syst., Vol. 4, M04001, doi:10.1029/2011MS000099 |
|---|
| 152 | !----------------------------------------------------------------------- |
|---|
| 153 | ! use physics_types , only: physics_dme_adjust ! This is for CESM/CAM |
|---|
| 154 | ! use cam_diagnostics, only: diag_phys_writeout ! This is for CESM/CAM |
|---|
| 155 | |
|---|
| 156 | implicit none |
|---|
| 157 | |
|---|
| 158 | integer, parameter :: r8 = selected_real_kind(12) |
|---|
| 159 | |
|---|
| 160 | ! |
|---|
| 161 | ! Input arguments - MODEL DEPENDENT |
|---|
| 162 | ! |
|---|
| 163 | integer, intent(in) :: pcols ! Set number of atmospheric columns |
|---|
| 164 | integer, intent(in) :: pver ! Set number of model levels |
|---|
| 165 | real(r8), intent(in) :: dtime ! Set model physics timestep |
|---|
| 166 | real(r8), intent(in) :: lat(pcols) ! Latitude |
|---|
| 167 | integer, intent(in) :: test ! Test number |
|---|
| 168 | |
|---|
| 169 | ! |
|---|
| 170 | ! Input/Output arguments |
|---|
| 171 | ! |
|---|
| 172 | ! pcols is the maximum number of vertical columns per 'chunk' of atmosphere |
|---|
| 173 | ! |
|---|
| 174 | real(r8), intent(inout) :: t(pcols,pver) ! Temperature at full-model level (K) |
|---|
| 175 | real(r8), intent(inout) :: q(pcols,pver) ! Specific Humidity at full-model level (kg/kg) |
|---|
| 176 | real(r8), intent(inout) :: u(pcols,pver) ! Zonal wind at full-model level (m/s) |
|---|
| 177 | real(r8), intent(inout) :: v(pcols,pver) ! Meridional wind at full-model level (m/s) |
|---|
| 178 | real(r8), intent(inout) :: pmid(pcols,pver) ! Pressure is full-model level (Pa) |
|---|
| 179 | real(r8), intent(inout) :: pint(pcols,pver+1) ! Pressure at model interfaces (Pa) |
|---|
| 180 | real(r8), intent(inout) :: pdel(pcols,pver) ! Layer thickness (Pa) |
|---|
| 181 | real(r8), intent(in) :: rpdel(pcols,pver) ! Reciprocal of layer thickness (1/Pa) |
|---|
| 182 | real(r8), intent(inout) :: ps(pcols) ! Surface Pressue (Pa) |
|---|
| 183 | |
|---|
| 184 | ! |
|---|
| 185 | ! Output arguments |
|---|
| 186 | ! |
|---|
| 187 | real(r8), intent(out) :: precl(pcols) ! Precipitation rate (m_water / s) |
|---|
| 188 | |
|---|
| 189 | ! |
|---|
| 190 | !---------------------------Local workspace----------------------------- |
|---|
| 191 | ! |
|---|
| 192 | |
|---|
| 193 | ! Integers for loops |
|---|
| 194 | |
|---|
| 195 | integer i,k ! Longitude, level indices |
|---|
| 196 | |
|---|
| 197 | ! Physical Constants - Many of these may be model dependent |
|---|
| 198 | |
|---|
| 199 | real(r8) gravit ! Gravity |
|---|
| 200 | real(r8) rair ! Gas constant for dry air |
|---|
| 201 | real(r8) cpair ! Specific heat of dry air |
|---|
| 202 | real(r8) latvap ! Latent heat of vaporization |
|---|
| 203 | real(r8) rh2o ! Gas constant for water vapor |
|---|
| 204 | real(r8) epsilo ! Ratio of gas constant for dry air to that for vapor |
|---|
| 205 | real(r8) zvir ! Constant for virtual temp. calc. =(rh2o/rair) - 1 |
|---|
| 206 | real(r8) a ! Reference Earth's Radius (m) |
|---|
| 207 | real(r8) omega ! Reference rotation rate of the Earth (s^-1) |
|---|
| 208 | real(r8) pi ! pi |
|---|
| 209 | |
|---|
| 210 | ! Simple Physics Specific Constants |
|---|
| 211 | |
|---|
| 212 | !++++++++ |
|---|
| 213 | real(r8) Tsurf(pcols) ! Sea Surface Temperature (constant for tropical cyclone) |
|---|
| 214 | !++++++++ Tsurf needs to be dependent on latitude for the |
|---|
| 215 | ! moist baroclinic wave test 4-3 with simple-physics, adjust |
|---|
| 216 | |
|---|
| 217 | real(r8) SST_tc ! Sea Surface Temperature for tropical cyclone test |
|---|
| 218 | real(r8) T0 ! Control temp for calculation of qsat |
|---|
| 219 | real(r8) e0 ! Saturation vapor pressure at T0 for calculation of qsat |
|---|
| 220 | real(r8) rhow ! Density of Liquid Water |
|---|
| 221 | real(r8) p0 ! Constant for calculation of potential temperature |
|---|
| 222 | real(r8) Cd0 ! Constant for calculating Cd from Smith and Vogl 2008 |
|---|
| 223 | real(r8) Cd1 ! Constant for calculating Cd from Smith and Vogl 2008 |
|---|
| 224 | real(r8) Cm ! Constant for calculating Cd from Smith and Vogl 2008 |
|---|
| 225 | real(r8) v20 ! Threshold wind speed for calculating Cd from Smith and Vogl 2008 |
|---|
| 226 | real(r8) C ! Drag coefficient for sensible heat and evaporation |
|---|
| 227 | real(r8) T00 ! Horizontal mean T at surface for moist baro test |
|---|
| 228 | real(r8) u0 ! Zonal wind constant for moist baro test |
|---|
| 229 | real(r8) latw ! halfwidth for for baro test |
|---|
| 230 | real(r8) eta0 ! Center of jets (hybrid) for baro test |
|---|
| 231 | real(r8) etav ! Auxiliary variable for baro test |
|---|
| 232 | real(r8) q0 ! Maximum specific humidity for baro test |
|---|
| 233 | |
|---|
| 234 | ! Physics Tendency Arrays |
|---|
| 235 | real(r8) dtdt(pcols,pver) ! Temperature tendency |
|---|
| 236 | real(r8) dqdt(pcols,pver) ! Specific humidity tendency |
|---|
| 237 | real(r8) dudt(pcols,pver) ! Zonal wind tendency |
|---|
| 238 | real(r8) dvdt(pcols,pver) ! Meridional wind tendency |
|---|
| 239 | |
|---|
| 240 | ! Temporary variables for tendency calculations |
|---|
| 241 | |
|---|
| 242 | real(r8) tmp ! Temporary |
|---|
| 243 | real(r8) qsat ! Saturation vapor pressure |
|---|
| 244 | real(r8) qsats ! Saturation vapor pressure of SST |
|---|
| 245 | |
|---|
| 246 | ! Variables for Boundary Layer Calculation |
|---|
| 247 | |
|---|
| 248 | real(r8) wind(pcols) ! Magnitude of Wind |
|---|
| 249 | real(r8) Cd(pcols) ! Drag coefficient for momentum |
|---|
| 250 | real(r8) Km(pcols,pver+1) ! Eddy diffusivity for boundary layer calculations |
|---|
| 251 | real(r8) Ke(pcols,pver+1) ! Eddy diffusivity for boundary layer calculations |
|---|
| 252 | real(r8) rho ! Density at lower/upper interface |
|---|
| 253 | real(r8) za(pcols) ! Heights at midpoints of first model level |
|---|
| 254 | real(r8) dlnpint ! Used for calculation of heights |
|---|
| 255 | real(r8) pbltop ! Top of boundary layer |
|---|
| 256 | real(r8) pblconst ! Constant for the calculation of the decay of diffusivity |
|---|
| 257 | real(r8) CA(pcols,pver) ! Matrix Coefficents for PBL Scheme |
|---|
| 258 | real(r8) CC(pcols,pver) ! Matrix Coefficents for PBL Scheme |
|---|
| 259 | real(r8) CE(pcols,pver+1) ! Matrix Coefficents for PBL Scheme |
|---|
| 260 | real(r8) CAm(pcols,pver) ! Matrix Coefficents for PBL Scheme |
|---|
| 261 | real(r8) CCm(pcols,pver) ! Matrix Coefficents for PBL Scheme |
|---|
| 262 | real(r8) CEm(pcols,pver+1) ! Matrix Coefficents for PBL Scheme |
|---|
| 263 | real(r8) CFu(pcols,pver+1) ! Matrix Coefficents for PBL Scheme |
|---|
| 264 | real(r8) CFv(pcols,pver+1) ! Matrix Coefficents for PBL Scheme |
|---|
| 265 | real(r8) CFt(pcols,pver+1) ! Matrix Coefficents for PBL Scheme |
|---|
| 266 | real(r8) CFq(pcols,pver+1) ! Matrix Coefficents for PBL Scheme |
|---|
| 267 | |
|---|
| 268 | |
|---|
| 269 | ! Variable for Dry Mass Adjustment, this dry air adjustment is necessary to |
|---|
| 270 | ! conserve the mass of the dry air |
|---|
| 271 | |
|---|
| 272 | real(r8) qini(pcols,pver) ! Initial specific humidity |
|---|
| 273 | |
|---|
| 274 | !=============================================================================== |
|---|
| 275 | ! |
|---|
| 276 | ! Physical Constants - MAY BE MODEL DEPENDENT |
|---|
| 277 | ! |
|---|
| 278 | !=============================================================================== |
|---|
| 279 | gravit = 9.80616_r8 ! Gravity (9.80616 m/s^2) |
|---|
| 280 | rair = 287.0_r8 ! Gas constant for dry air: 287 J/(kg K) |
|---|
| 281 | cpair = 1.0045e3_r8 ! Specific heat of dry air: here we use 1004.5 J/(kg K) |
|---|
| 282 | latvap = 2.5e6_r8 ! Latent heat of vaporization (J/kg) |
|---|
| 283 | rh2o = 461.5_r8 ! Gas constant for water vapor: 461.5 J/(kg K) |
|---|
| 284 | epsilo = rair/rh2o ! Ratio of gas constant for dry air to that for vapor |
|---|
| 285 | zvir = (rh2o/rair) - 1._r8 ! Constant for virtual temp. calc. =(rh2o/rair) - 1 is approx. 0.608 |
|---|
| 286 | a = 6371220.0_r8 ! Reference Earth's Radius (m) |
|---|
| 287 | omega = 7.29212d-5 ! Reference rotation rate of the Earth (s^-1) |
|---|
| 288 | pi = 4._r8*atan(1._r8) ! pi |
|---|
| 289 | |
|---|
| 290 | !=============================================================================== |
|---|
| 291 | ! |
|---|
| 292 | ! Local Constants for Simple Physics |
|---|
| 293 | ! |
|---|
| 294 | !=============================================================================== |
|---|
| 295 | C = 0.0011_r8 ! From Smith and Vogl 2008 |
|---|
| 296 | SST_tc = 302.15_r8 ! Constant Value for SST for tropical cyclone test |
|---|
| 297 | T0 = 273.16_r8 ! control temp for calculation of qsat |
|---|
| 298 | e0 = 610.78_r8 ! saturation vapor pressure at T0 for calculation of qsat |
|---|
| 299 | rhow = 1000.0_r8 ! Density of Liquid Water |
|---|
| 300 | Cd0 = 0.0007_r8 ! Constant for Cd calc. Smith and Vogl 2008 |
|---|
| 301 | Cd1 = 0.000065_r8 ! Constant for Cd calc. Smith and Vogl 2008 |
|---|
| 302 | Cm = 0.002_r8 ! Constant for Cd calc. Smith and Vogl 2008 |
|---|
| 303 | v20 = 20.0_r8 ! Threshold wind speed for calculating Cd from Smith and Vogl 2008 |
|---|
| 304 | p0 = 100000.0_r8 ! Constant for potential temp calculation |
|---|
| 305 | pbltop = 85000._r8 ! Top of boundary layer |
|---|
| 306 | pblconst = 10000._r8 ! Constant for the calculation of the decay of diffusivity |
|---|
| 307 | T00 = 288.0_r8 ! Horizontal mean T at surface for moist baro test |
|---|
| 308 | u0 = 35.0_r8 ! Zonal wind constant for moist baro test |
|---|
| 309 | latw = 2.0_r8*pi/9.0_r8 ! Halfwidth for for baro test |
|---|
| 310 | eta0 = 0.252_r8 ! Center of jets (hybrid) for baro test |
|---|
| 311 | etav = (1._r8-eta0)*0.5_r8*pi ! Auxiliary variable for baro test |
|---|
| 312 | q0 = 0.021 ! Maximum specific humidity for baro test |
|---|
| 313 | |
|---|
| 314 | !=============================================================================== |
|---|
| 315 | ! |
|---|
| 316 | ! Definition of local arrays |
|---|
| 317 | ! |
|---|
| 318 | !=============================================================================== |
|---|
| 319 | ! |
|---|
| 320 | ! Calculate hydrostatic height za of the lowest model level |
|---|
| 321 | ! |
|---|
| 322 | do i=1,pcols |
|---|
| 323 | dlnpint = log(ps(i)) - log(pint(i,pver)) ! ps(i) is identical to pint(i,pver+1), note: this is the correct sign (corrects typo in JAMES paper) |
|---|
| 324 | za(i) = rair/gravit*t(i,pver)*(1._r8+zvir*q(i,pver))*0.5_r8*dlnpint |
|---|
| 325 | end do |
|---|
| 326 | ! |
|---|
| 327 | ! Set Initial Specific Humidity - For dry mass adjustment at the end |
|---|
| 328 | ! |
|---|
| 329 | qini(:pcols,:pver) = q(:pcols,:pver) |
|---|
| 330 | !-------------------------------------------------------------- |
|---|
| 331 | ! Set Sea Surface Temperature (constant for tropical cyclone) |
|---|
| 332 | ! Tsurf needs to be dependent on latitude for the |
|---|
| 333 | ! moist baroclinic wave test 4-3 with simple-physics |
|---|
| 334 | !-------------------------------------------------------------- |
|---|
| 335 | if (test .eq. 1) then ! moist baroclinic wave with simple-physics |
|---|
| 336 | do i=1,pcols |
|---|
| 337 | Tsurf(i) = (T00 + pi*u0/rair * 1.5_r8 * sin(etav) * (cos(etav))**0.5_r8 * & |
|---|
| 338 | ((-2._r8*(sin(lat(i)))**6 * ((cos(lat(i)))**2 + 1._r8/3._r8) + 10._r8/63._r8)* & |
|---|
| 339 | u0 * (cos(etav))**1.5_r8 + & |
|---|
| 340 | (8._r8/5._r8*(cos(lat(i)))**3 * ((sin(lat(i)))**2 + 2._r8/3._r8) - pi/4._r8)*a*omega*0.5_r8 ))/ & |
|---|
| 341 | (1._r8+zvir*q0*exp(-(lat(i)/latw)**4)) |
|---|
| 342 | |
|---|
| 343 | end do |
|---|
| 344 | else |
|---|
| 345 | do i=1,pcols ! constant SST for the tropical cyclone test case |
|---|
| 346 | Tsurf(i) = SST_tc |
|---|
| 347 | end do |
|---|
| 348 | end if |
|---|
| 349 | |
|---|
| 350 | !=============================================================================== |
|---|
| 351 | ! |
|---|
| 352 | ! Set initial physics time tendencies and precipitation field to zero |
|---|
| 353 | ! |
|---|
| 354 | !=============================================================================== |
|---|
| 355 | dtdt(:pcols,:pver) = 0._r8 ! initialize temperature tendency with zero |
|---|
| 356 | dqdt(:pcols,:pver) = 0._r8 ! initialize specific humidity tendency with zero |
|---|
| 357 | dudt(:pcols,:pver) = 0._r8 ! initialize zonal wind tendency with zero |
|---|
| 358 | dvdt(:pcols,:pver) = 0._r8 ! initialize meridional wind tendency with zero |
|---|
| 359 | precl(:pcols) = 0._r8 ! initialize precipitation rate with zero |
|---|
| 360 | |
|---|
| 361 | !=============================================================================== |
|---|
| 362 | ! |
|---|
| 363 | ! Large-Scale Condensation and Precipitation Rate |
|---|
| 364 | ! |
|---|
| 365 | !=============================================================================== |
|---|
| 366 | ! |
|---|
| 367 | ! Calculate Tendencies |
|---|
| 368 | ! |
|---|
| 369 | do k=1,pver |
|---|
| 370 | do i=1,pcols |
|---|
| 371 | qsat = epsilo*e0/pmid(i,k)*exp(-latvap/rh2o*((1._r8/t(i,k))-1._r8/T0)) ! saturation specific humidity |
|---|
| 372 | ! out_i(i,llm+1-k)=q(i,k)-qsat |
|---|
| 373 | if (q(i,k) > qsat) then ! saturated? |
|---|
| 374 | tmp = 1._r8/dtime*(q(i,k)-qsat)/(1._r8+(latvap/cpair)*(epsilo*latvap*qsat/(rair*t(i,k)**2))) |
|---|
| 375 | dtdt(i,k) = dtdt(i,k)+latvap/cpair*tmp |
|---|
| 376 | dqdt(i,k) = dqdt(i,k)-tmp |
|---|
| 377 | precl(i) = precl(i) + tmp*pdel(i,k)/(gravit*rhow) ! precipitation rate, computed via a vertical integral |
|---|
| 378 | ! corrected in version 1.3 |
|---|
| 379 | end if |
|---|
| 380 | end do |
|---|
| 381 | end do |
|---|
| 382 | ! |
|---|
| 383 | ! Update moisture and temperature fields from Large-Scale Precipitation Scheme |
|---|
| 384 | ! |
|---|
| 385 | do k=1,pver |
|---|
| 386 | do i=1,pcols |
|---|
| 387 | t(i,k) = t(i,k) + dtdt(i,k)*dtime ! update the state variables T and q |
|---|
| 388 | q(i,k) = q(i,k) + dqdt(i,k)*dtime |
|---|
| 389 | end do |
|---|
| 390 | end do |
|---|
| 391 | |
|---|
| 392 | !=============================================================================== |
|---|
| 393 | ! Send variables to history file - THIS PROCESS WILL BE MODEL SPECIFIC |
|---|
| 394 | ! |
|---|
| 395 | ! note: The variables, as done in many GCMs, are written to the history file |
|---|
| 396 | ! after the moist physics process. This ensures that the moisture fields |
|---|
| 397 | ! are somewhat in equilibrium. |
|---|
| 398 | !=============================================================================== |
|---|
| 399 | ! call diag_phys_writeout(state) ! This is for CESM/CAM |
|---|
| 400 | |
|---|
| 401 | !=============================================================================== |
|---|
| 402 | ! |
|---|
| 403 | ! Turbulent mixing coefficients for the PBL mixing of horizontal momentum, |
|---|
| 404 | ! sensible heat and latent heat |
|---|
| 405 | ! |
|---|
| 406 | ! We are using Simplified Ekman theory to compute the diffusion coefficients |
|---|
| 407 | ! Kx for the boundary-layer mixing. The Kx values are calculated at each time step |
|---|
| 408 | ! and in each column. |
|---|
| 409 | ! |
|---|
| 410 | !=============================================================================== |
|---|
| 411 | ! |
|---|
| 412 | ! Compute magnitude of the wind and drag coeffcients for turbulence scheme: |
|---|
| 413 | ! they depend on the conditions at the lowest model level and stay constant |
|---|
| 414 | ! up to the 850 hPa level. Above this level the coefficients are decreased |
|---|
| 415 | ! and tapered to zero. At the 700 hPa level the strength of the K coefficients |
|---|
| 416 | ! is about 10% of the maximum strength. |
|---|
| 417 | ! |
|---|
| 418 | do i=1,pcols |
|---|
| 419 | wind(i) = sqrt(u(i,pver)**2+v(i,pver)**2) ! wind magnitude at the lowest level |
|---|
| 420 | end do |
|---|
| 421 | do i=1,pcols |
|---|
| 422 | Ke(i,pver+1) = C*wind(i)*za(i) |
|---|
| 423 | if( wind(i) .lt. v20) then |
|---|
| 424 | Cd(i) = Cd0+Cd1*wind(i) |
|---|
| 425 | Km(i,pver+1) = Cd(i)*wind(i)*za(i) |
|---|
| 426 | else |
|---|
| 427 | Cd(i) = Cm |
|---|
| 428 | Km(i,pver+1) = Cm*wind(i)*za(i) |
|---|
| 429 | endif |
|---|
| 430 | end do |
|---|
| 431 | |
|---|
| 432 | do k=1,pver |
|---|
| 433 | do i=1,pcols |
|---|
| 434 | if( pint(i,k) .ge. pbltop) then |
|---|
| 435 | Km(i,k) = Km(i,pver+1) ! constant Km below 850 hPa level |
|---|
| 436 | Ke(i,k) = Ke(i,pver+1) ! constant Ke below 850 hPa level |
|---|
| 437 | else |
|---|
| 438 | Km(i,k) = Km(i,pver+1)*exp(-(pbltop-pint(i,k))**2/(pblconst)**2) ! Km tapered to 0 |
|---|
| 439 | Ke(i,k) = Ke(i,pver+1)*exp(-(pbltop-pint(i,k))**2/(pblconst)**2) ! Ke tapered to 0 |
|---|
| 440 | end if |
|---|
| 441 | end do |
|---|
| 442 | end do |
|---|
| 443 | |
|---|
| 444 | |
|---|
| 445 | !=============================================================================== |
|---|
| 446 | ! Update the state variables u, v, t, q with the surface fluxes at the |
|---|
| 447 | ! lowest model level, this is done with an implicit approach |
|---|
| 448 | ! see Reed and Jablonowski (JAMES, 2012) |
|---|
| 449 | ! |
|---|
| 450 | ! Sea Surface Temperature Tsurf is constant for tropical cyclone test 5-1 |
|---|
| 451 | ! Tsurf needs to be dependent on latitude for the |
|---|
| 452 | ! moist baroclinic wave test 4-3 with simple-physics, adjust |
|---|
| 453 | !=============================================================================== |
|---|
| 454 | |
|---|
| 455 | do i=1,pcols |
|---|
| 456 | qsats = epsilo*e0/ps(i)*exp(-latvap/rh2o*((1._r8/Tsurf(i))-1._r8/T0)) ! saturation specific humidity at the surface |
|---|
| 457 | dudt(i,pver) = dudt(i,pver) + (u(i,pver) & |
|---|
| 458 | /(1._r8+Cd(i)*wind(i)*dtime/za(i))-u(i,pver))/dtime |
|---|
| 459 | dvdt(i,pver) = dvdt(i,pver) + (v(i,pver) & |
|---|
| 460 | /(1._r8+Cd(i)*wind(i)*dtime/za(i))-v(i,pver))/dtime |
|---|
| 461 | u(i,pver) = u(i,pver)/(1._r8+Cd(i)*wind(i)*dtime/za(i)) |
|---|
| 462 | v(i,pver) = v(i,pver)/(1._r8+Cd(i)*wind(i)*dtime/za(i)) |
|---|
| 463 | dtdt(i,pver) = dtdt(i,pver) +((t(i,pver)+C*wind(i)*Tsurf(i)*dtime/za(i)) & |
|---|
| 464 | /(1._r8+C*wind(i)*dtime/za(i))-t(i,pver))/dtime |
|---|
| 465 | t(i,pver) = (t(i,pver)+C*wind(i)*Tsurf(i)*dtime/za(i)) & |
|---|
| 466 | /(1._r8+C*wind(i)*dtime/za(i)) |
|---|
| 467 | dqdt(i,pver) = dqdt(i,pver) +((q(i,pver)+C*wind(i)*qsats*dtime/za(i)) & |
|---|
| 468 | /(1._r8+C*wind(i)*dtime/za(i))-q(i,pver))/dtime |
|---|
| 469 | q(i,pver) = (q(i,pver)+C*wind(i)*qsats*dtime/za(i))/(1._r8+C*wind(i)*dtime/za(i)) |
|---|
| 470 | end do |
|---|
| 471 | !=============================================================================== |
|---|
| 472 | |
|---|
| 473 | |
|---|
| 474 | !=============================================================================== |
|---|
| 475 | ! Boundary layer mixing, see Reed and Jablonowski (JAMES, 2012) |
|---|
| 476 | !=============================================================================== |
|---|
| 477 | ! Calculate Diagonal Variables for Implicit PBL Scheme |
|---|
| 478 | ! |
|---|
| 479 | do k=1,pver-1 |
|---|
| 480 | do i=1,pcols |
|---|
| 481 | rho = (pint(i,k+1)/(rair*(t(i,k+1)+t(i,k))/2.0_r8)) |
|---|
| 482 | CAm(i,k) = rpdel(i,k)*dtime*gravit*gravit*Km(i,k+1)*rho*rho & |
|---|
| 483 | /(pmid(i,k+1)-pmid(i,k)) |
|---|
| 484 | CCm(i,k+1) = rpdel(i,k+1)*dtime*gravit*gravit*Km(i,k+1)*rho*rho & |
|---|
| 485 | /(pmid(i,k+1)-pmid(i,k)) |
|---|
| 486 | CA(i,k) = rpdel(i,k)*dtime*gravit*gravit*Ke(i,k+1)*rho*rho & |
|---|
| 487 | /(pmid(i,k+1)-pmid(i,k)) |
|---|
| 488 | CC(i,k+1) = rpdel(i,k+1)*dtime*gravit*gravit*Ke(i,k+1)*rho*rho & |
|---|
| 489 | /(pmid(i,k+1)-pmid(i,k)) |
|---|
| 490 | end do |
|---|
| 491 | end do |
|---|
| 492 | do i=1,pcols |
|---|
| 493 | CAm(i,pver) = 0._r8 |
|---|
| 494 | CCm(i,1) = 0._r8 |
|---|
| 495 | CEm(i,pver+1) = 0._r8 |
|---|
| 496 | CA(i,pver) = 0._r8 |
|---|
| 497 | CC(i,1) = 0._r8 |
|---|
| 498 | CE(i,pver+1) = 0._r8 |
|---|
| 499 | CFu(i,pver+1) = 0._r8 |
|---|
| 500 | CFv(i,pver+1) = 0._r8 |
|---|
| 501 | CFt(i,pver+1) = 0._r8 |
|---|
| 502 | CFq(i,pver+1) = 0._r8 |
|---|
| 503 | end do |
|---|
| 504 | do i=1,pcols |
|---|
| 505 | do k=pver,1,-1 |
|---|
| 506 | CE(i,k) = CC(i,k)/(1._r8+CA(i,k)+CC(i,k)-CA(i,k)*CE(i,k+1)) |
|---|
| 507 | CEm(i,k) = CCm(i,k)/(1._r8+CAm(i,k)+CCm(i,k)-CAm(i,k)*CEm(i,k+1)) |
|---|
| 508 | CFu(i,k) = (u(i,k)+CAm(i,k)*CFu(i,k+1)) & |
|---|
| 509 | /(1._r8+CAm(i,k)+CCm(i,k)-CAm(i,k)*CEm(i,k+1)) |
|---|
| 510 | CFv(i,k) = (v(i,k)+CAm(i,k)*CFv(i,k+1)) & |
|---|
| 511 | /(1._r8+CAm(i,k)+CCm(i,k)-CAm(i,k)*CEm(i,k+1)) |
|---|
| 512 | CFt(i,k) = ((p0/pmid(i,k))**(rair/cpair)*t(i,k)+CA(i,k)*CFt(i,k+1)) & |
|---|
| 513 | /(1._r8+CA(i,k)+CC(i,k)-CA(i,k)*CE(i,k+1)) |
|---|
| 514 | CFq(i,k) = (q(i,k)+CA(i,k)*CFq(i,k+1)) & |
|---|
| 515 | /(1._r8+CA(i,k)+CC(i,k)-CA(i,k)*CE(i,k+1)) |
|---|
| 516 | end do |
|---|
| 517 | end do |
|---|
| 518 | |
|---|
| 519 | ! |
|---|
| 520 | ! Calculate the updated temperature, specific humidity and horizontal wind |
|---|
| 521 | ! |
|---|
| 522 | ! First we need to calculate the updates at the top model level |
|---|
| 523 | ! |
|---|
| 524 | do i=1,pcols |
|---|
| 525 | dudt(i,1) = dudt(i,1)+(CFu(i,1)-u(i,1))/dtime |
|---|
| 526 | dvdt(i,1) = dvdt(i,1)+(CFv(i,1)-v(i,1))/dtime |
|---|
| 527 | u(i,1) = CFu(i,1) |
|---|
| 528 | v(i,1) = CFv(i,1) |
|---|
| 529 | dtdt(i,1) = dtdt(i,1)+(CFt(i,1)*(pmid(i,1)/p0)**(rair/cpair)-t(i,1))/dtime ! corrected in version 1.3 |
|---|
| 530 | t(i,1) = CFt(i,1)*(pmid(i,1)/p0)**(rair/cpair) |
|---|
| 531 | dqdt(i,1) = dqdt(i,1)+(CFq(i,1)-q(i,1))/dtime |
|---|
| 532 | q(i,1) = CFq(i,1) |
|---|
| 533 | end do |
|---|
| 534 | ! |
|---|
| 535 | ! Loop over the remaining level |
|---|
| 536 | ! |
|---|
| 537 | do i=1,pcols |
|---|
| 538 | do k=2,pver |
|---|
| 539 | dudt(i,k) = dudt(i,k)+(CEm(i,k)*u(i,k-1)+CFu(i,k)-u(i,k))/dtime |
|---|
| 540 | dvdt(i,k) = dvdt(i,k)+(CEm(i,k)*v(i,k-1)+CFv(i,k)-v(i,k))/dtime |
|---|
| 541 | u(i,k) = CEm(i,k)*u(i,k-1)+CFu(i,k) |
|---|
| 542 | v(i,k) = CEm(i,k)*v(i,k-1)+CFv(i,k) |
|---|
| 543 | dtdt(i,k) = dtdt(i,k)+((CE(i,k)*t(i,k-1) & |
|---|
| 544 | *(p0/pmid(i,k-1))**(rair/cpair)+CFt(i,k)) & |
|---|
| 545 | *(pmid(i,k)/p0)**(rair/cpair)-t(i,k))/dtime |
|---|
| 546 | t(i,k) = (CE(i,k)*t(i,k-1)*(p0/pmid(i,k-1))**(rair/cpair)+CFt(i,k)) & |
|---|
| 547 | *(pmid(i,k)/p0)**(rair/cpair) |
|---|
| 548 | dqdt(i,k) = dqdt(i,k)+(CE(i,k)*q(i,k-1)+CFq(i,k)-q(i,k))/dtime |
|---|
| 549 | q(i,k) = CE(i,k)*q(i,k-1)+CFq(i,k) |
|---|
| 550 | end do |
|---|
| 551 | end do |
|---|
| 552 | |
|---|
| 553 | !=============================================================================== |
|---|
| 554 | ! |
|---|
| 555 | ! Dry Mass Adjustment - THIS PROCESS WILL BE MODEL SPECIFIC |
|---|
| 556 | ! |
|---|
| 557 | ! note: Care needs to be taken to ensure that the model conserves the total |
|---|
| 558 | ! dry air mass. Add your own routine here. |
|---|
| 559 | !=============================================================================== |
|---|
| 560 | ! call physics_dme_adjust(state, tend, qini, dtime) ! This is for CESM/CAM |
|---|
| 561 | |
|---|
| 562 | return |
|---|
| 563 | end subroutine simple_physics |
|---|
| 564 | |
|---|
![(please configure the [header_logo] section in trac.ini)](/dynamico/chrome/site/your_project_logo.png){kind=logo}
