[388] | 1 | MODULE dcmip2016_supercell_mod |
---|
| 2 | |
---|
| 3 | !======================================================================= |
---|
| 4 | ! |
---|
| 5 | ! Date: April 22, 2016 |
---|
| 6 | ! |
---|
| 7 | ! Functions for setting up idealized initial conditions for the |
---|
| 8 | ! Klemp et al. supercell test. Before sampling the result, |
---|
| 9 | ! supercell_test_init() must be called. |
---|
| 10 | ! |
---|
| 11 | ! SUBROUTINE supercell_test( |
---|
| 12 | ! lon,lat,p,z,zcoords,u,v,t,thetav,ps,rho,q,pert) |
---|
| 13 | ! |
---|
| 14 | ! Given a point specified by: |
---|
| 15 | ! lon longitude (radians) |
---|
| 16 | ! lat latitude (radians) |
---|
| 17 | ! p/z pressure (Pa) / height (m) |
---|
| 18 | ! zcoords 1 if z is specified, 0 if p is specified |
---|
| 19 | ! pert 1 if thermal perturbation included, 0 if not |
---|
| 20 | ! |
---|
| 21 | ! the functions will return: |
---|
| 22 | ! p pressure if z is specified (Pa) |
---|
| 23 | ! z geopotential height if p is specified (m) |
---|
| 24 | ! u zonal wind (m s^-1) |
---|
| 25 | ! v meridional wind (m s^-1) |
---|
| 26 | ! t temperature (K) |
---|
| 27 | ! thetav virtual potential temperature (K) |
---|
| 28 | ! ps surface pressure (Pa) |
---|
| 29 | ! rho density (kj m^-3) |
---|
| 30 | ! q water vapor mixing ratio (kg/kg) |
---|
| 31 | ! |
---|
| 32 | ! Author: Paul Ullrich |
---|
| 33 | ! University of California, Davis |
---|
| 34 | ! Email: paullrich@ucdavis.edu |
---|
| 35 | ! |
---|
| 36 | ! Based on a code by Joseph Klemp |
---|
| 37 | ! (National Center for Atmospheric Research) |
---|
| 38 | ! |
---|
| 39 | !======================================================================= |
---|
| 40 | |
---|
| 41 | IMPLICIT NONE |
---|
| 42 | |
---|
| 43 | !======================================================================= |
---|
| 44 | ! Physical constants |
---|
| 45 | !======================================================================= |
---|
| 46 | |
---|
| 47 | REAL(8), PARAMETER :: & |
---|
| 48 | a = 6371220.0d0, & ! Reference Earth's Radius (m) |
---|
| 49 | Rd = 287.0d0, & ! Ideal gas const dry air (J kg^-1 K^1) |
---|
| 50 | g = 9.80616d0, & ! Gravity (m s^2) |
---|
| 51 | cp = 1004.5d0, & ! Specific heat capacity (J kg^-1 K^1) |
---|
| 52 | Lvap = 2.5d6, & ! Latent heat of vaporization of water |
---|
| 53 | Rvap = 461.5d0, & ! Ideal gas constnat for water vapor |
---|
| 54 | Mvap = 0.608d0, & ! Ratio of molar mass of dry air/water |
---|
| 55 | pi = 3.14159265358979d0, & ! pi |
---|
| 56 | p0 = 100000.0d0, & ! surface pressure (Pa) |
---|
| 57 | kappa = 2.d0/7.d0, & ! Ratio of Rd to cp |
---|
| 58 | omega = 7.29212d-5, & ! Reference rotation rate of the Earth (s^-1) |
---|
| 59 | deg2rad = pi/180.d0 ! Conversion factor of degrees to radians |
---|
| 60 | |
---|
| 61 | !======================================================================= |
---|
| 62 | ! Test case parameters |
---|
| 63 | !======================================================================= |
---|
| 64 | INTEGER(4), PARAMETER :: & |
---|
| 65 | nz = 30 , & ! number of vertical levels in init |
---|
| 66 | nphi = 16 ! number of meridional points in init |
---|
| 67 | |
---|
| 68 | REAL(8), PARAMETER :: & |
---|
| 69 | z1 = 0.0d0 , & ! lower sample altitude |
---|
| 70 | z2 = 50000.0d0 ! upper sample altitude |
---|
| 71 | |
---|
| 72 | REAL(8), PARAMETER :: & |
---|
| 73 | X = 120.d0 , & ! Earth reduction factor |
---|
| 74 | theta0 = 300.d0 , & ! theta at the equatorial surface |
---|
| 75 | theta_tr = 343.d0 , & ! theta at the tropopause |
---|
| 76 | z_tr = 12000.d0 , & ! altitude at the tropopause |
---|
| 77 | T_tr = 213.d0 , & ! temperature at the tropopause |
---|
| 78 | pseq = 100000.0d0 ! surface pressure at equator (Pa) |
---|
| 79 | !pseq = 95690.0d0 ! surface pressure at equator (Pa) |
---|
| 80 | |
---|
| 81 | REAL(8), PARAMETER :: & |
---|
| 82 | us = 30.d0 , & ! maximum zonal wind velocity |
---|
| 83 | uc = 15.d0 , & ! coordinate reference velocity |
---|
| 84 | zs = 5000.d0 , & ! lower altitude of maximum velocity |
---|
| 85 | zt = 1000.d0 ! transition distance of velocity |
---|
| 86 | |
---|
| 87 | REAL(8), PARAMETER :: & |
---|
| 88 | pert_dtheta = 3.d0 , & ! perturbation magnitude |
---|
| 89 | pert_lonc = 0.d0 , & ! perturbation longitude |
---|
| 90 | pert_latc = 0.d0 , & ! perturbation latitude |
---|
| 91 | pert_rh = 10000.d0 * X , & ! perturbation horiz. halfwidth |
---|
| 92 | pert_zc = 1500.d0 , & ! perturbation center altitude |
---|
| 93 | pert_rz = 1500.d0 ! perturbation vert. halfwidth |
---|
| 94 | |
---|
| 95 | !----------------------------------------------------------------------- |
---|
| 96 | ! Coefficients computed from initialization |
---|
| 97 | !----------------------------------------------------------------------- |
---|
| 98 | INTEGER(4) :: initialized = 0 |
---|
| 99 | |
---|
| 100 | REAL(8), DIMENSION(nphi) :: phicoord |
---|
| 101 | REAL(8), DIMENSION(nz) :: zcoord |
---|
| 102 | REAL(8), DIMENSION(nphi,nz) :: thetavyz |
---|
| 103 | REAL(8), DIMENSION(nphi,nz) :: exneryz |
---|
| 104 | REAL(8), DIMENSION(nz) :: qveq |
---|
| 105 | |
---|
| 106 | CONTAINS |
---|
| 107 | |
---|
| 108 | !======================================================================= |
---|
| 109 | ! Generate the supercell initial conditions |
---|
| 110 | !======================================================================= |
---|
| 111 | SUBROUTINE supercell_init() & |
---|
| 112 | BIND(c, name = "supercell_init") |
---|
| 113 | |
---|
| 114 | IMPLICIT NONE |
---|
| 115 | |
---|
| 116 | ! d/dphi and int(dphi) operators |
---|
| 117 | REAL(8), DIMENSION(nphi,nphi) :: ddphi, intphi |
---|
| 118 | |
---|
| 119 | ! d/dz and int(dz) operators |
---|
| 120 | REAL(8), DIMENSION(nz, nz) :: ddz, intz |
---|
| 121 | |
---|
| 122 | ! Buffer matrices for computing SVD of d/dphi operator |
---|
| 123 | REAL(8), DIMENSION(nphi,nphi) :: ddphibak |
---|
| 124 | REAL(8), DIMENSION(nphi,nphi) :: svdpu, svdpvt |
---|
| 125 | REAL(8), DIMENSION(nphi) :: svdps |
---|
| 126 | REAL(8), DIMENSION(5*nphi) :: pwork |
---|
| 127 | |
---|
| 128 | ! Buffer matrices for computing SVD of d/dz operator |
---|
| 129 | REAL(8), DIMENSION(nz, nz) :: ddzbak |
---|
| 130 | REAL(8), DIMENSION(nz, nz) :: svdzu, svdzvt |
---|
| 131 | REAL(8), DIMENSION(nz) :: svdzs |
---|
| 132 | REAL(8), DIMENSION(5*nz) :: zwork |
---|
| 133 | |
---|
| 134 | ! Buffer data for calculation of SVD |
---|
| 135 | INTEGER(4) :: lwork, info |
---|
| 136 | |
---|
| 137 | ! Sampled values of ueq**2 and d/dz(ueq**2) |
---|
| 138 | REAL(8), DIMENSION(nphi, nz) :: ueq2, dueq2 |
---|
| 139 | |
---|
| 140 | ! Buffer matrices for iteration |
---|
| 141 | REAL(8), DIMENSION(nphi, nz) :: phicoordmat, dztheta, rhs, irhs |
---|
| 142 | |
---|
| 143 | ! Buffer for sampled potential temperature at equator |
---|
| 144 | REAL(8), DIMENSION(nz) :: thetaeq |
---|
| 145 | |
---|
| 146 | ! Buffer for computed equatorial Exner pressure and relative humidity |
---|
| 147 | REAL(8), DIMENSION(nz) :: exnereq, H |
---|
| 148 | |
---|
| 149 | ! Variables for calculation of equatorial profile |
---|
| 150 | REAL(8) :: exnereqs, p, T, qvs, qv |
---|
| 151 | |
---|
| 152 | ! Error metric |
---|
| 153 | REAL(8) :: err |
---|
| 154 | |
---|
| 155 | ! Loop indices |
---|
| 156 | INTEGER(4) :: i, k, iter |
---|
| 157 | |
---|
| 158 | ! Chebyshev nodes in the phi direction |
---|
| 159 | do i = 1, nphi |
---|
| 160 | phicoord(i) = - cos(dble(i-1) * pi / dble(nphi-1)) |
---|
| 161 | phicoord(i) = 0.25d0 * pi * (phicoord(i) + 1.0d0) |
---|
| 162 | end do |
---|
| 163 | |
---|
| 164 | ! Matrix of phis |
---|
| 165 | do k = 1, nz |
---|
| 166 | phicoordmat(:,k) = phicoord |
---|
| 167 | end do |
---|
| 168 | |
---|
| 169 | ! Chebyshev nodes in the z direction |
---|
| 170 | do k = 1, nz |
---|
| 171 | zcoord(k) = - cos(dble(k-1) * pi / dble(nz-1)) |
---|
| 172 | zcoord(k) = z1 + 0.5d0*(z2-z1)*(zcoord(k)+1.0d0) |
---|
| 173 | end do |
---|
| 174 | |
---|
| 175 | ! Compute the d/dphi operator |
---|
| 176 | do i = 1, nphi |
---|
| 177 | call diff_lagrangian_polynomial_coeffs( & |
---|
| 178 | nphi, phicoord, ddphi(:,i), phicoord(i)) |
---|
| 179 | end do |
---|
| 180 | |
---|
| 181 | ! Zero derivative at pole |
---|
| 182 | ddphi(:,nphi) = 0.0d0 |
---|
| 183 | |
---|
| 184 | ! Compute the d/dz operator |
---|
| 185 | do k = 1, nz |
---|
| 186 | call diff_lagrangian_polynomial_coeffs( & |
---|
| 187 | nz, zcoord, ddz(:,k), zcoord(k)) |
---|
| 188 | end do |
---|
| 189 | |
---|
| 190 | ! Compute the int(dphi) operator via pseudoinverse |
---|
| 191 | lwork = 5*nphi |
---|
| 192 | |
---|
| 193 | ddphibak = ddphi |
---|
| 194 | call DGESVD('A', 'A', & |
---|
| 195 | nphi, nphi, ddphibak, nphi, & |
---|
| 196 | svdps, svdpu, nphi, svdpvt, nphi, & |
---|
| 197 | pwork, lwork, info) |
---|
| 198 | |
---|
| 199 | if (info .ne. 0) then |
---|
| 200 | write(*,*) 'Unable to compute SVD of d/dphi matrix' |
---|
| 201 | stop |
---|
| 202 | end if |
---|
| 203 | |
---|
| 204 | do i = 1, nphi |
---|
| 205 | if (abs(svdps(i)) .le. 1.0d-12) then |
---|
| 206 | call DSCAL(nphi, 0.0d0, svdpu(1,i), 1) |
---|
| 207 | else |
---|
| 208 | call DSCAL(nphi, 1.0d0 / svdps(i), svdpu(1,i), 1) |
---|
| 209 | end if |
---|
| 210 | end do |
---|
| 211 | call DGEMM('T', 'T', & |
---|
| 212 | nphi, nphi, nphi, 1.0d0, svdpvt, nphi, svdpu, nphi, 0.0d0, & |
---|
| 213 | intphi, nphi) |
---|
| 214 | |
---|
| 215 | ! Compute the int(dz) operator via pseudoinverse |
---|
| 216 | lwork = 5*nz |
---|
| 217 | |
---|
| 218 | ddzbak = ddz |
---|
| 219 | call DGESVD('A', 'A', & |
---|
| 220 | nz, nz, ddzbak, nz, & |
---|
| 221 | svdzs, svdzu, nz, svdzvt, nz, & |
---|
| 222 | zwork, lwork, info) |
---|
| 223 | |
---|
| 224 | if (info .ne. 0) then |
---|
| 225 | write(*,*) 'Unable to compute SVD of d/dz matrix' |
---|
| 226 | stop |
---|
| 227 | end if |
---|
| 228 | |
---|
| 229 | do i = 1, nz |
---|
| 230 | if (abs(svdzs(i)) .le. 1.0d-12) then |
---|
| 231 | call DSCAL(nz, 0.0d0, svdzu(1,i), 1) |
---|
| 232 | else |
---|
| 233 | call DSCAL(nz, 1.0d0 / svdzs(i), svdzu(1,i), 1) |
---|
| 234 | end if |
---|
| 235 | end do |
---|
| 236 | call DGEMM('T', 'T', & |
---|
| 237 | nz, nz, nz, 1.0d0, svdzvt, nz, svdzu, nz, 0.0d0, & |
---|
| 238 | intz, nz) |
---|
| 239 | |
---|
| 240 | ! Sample the equatorial velocity field and its derivative |
---|
| 241 | do k = 1, nz |
---|
| 242 | ueq2(1,k) = zonal_velocity(zcoord(k), 0.0d0) |
---|
| 243 | ueq2(1,k) = ueq2(1,k)**2 |
---|
| 244 | end do |
---|
| 245 | do k = 1, nz |
---|
| 246 | dueq2(1,k) = dot_product(ddz(:,k), ueq2(1,:)) |
---|
| 247 | end do |
---|
| 248 | do i = 2, nphi |
---|
| 249 | ueq2(i,:) = ueq2(1,:) |
---|
| 250 | dueq2(i,:) = dueq2(1,:) |
---|
| 251 | end do |
---|
| 252 | |
---|
| 253 | ! Initialize potential temperature at equator |
---|
| 254 | do k = 1, nz |
---|
| 255 | thetaeq(k) = equator_theta(zcoord(k)) |
---|
| 256 | H(k) = equator_relative_humidity(zcoord(k)) |
---|
| 257 | end do |
---|
| 258 | thetavyz(1,:) = thetaeq |
---|
| 259 | |
---|
| 260 | ! Exner pressure at the equatorial surface |
---|
| 261 | exnereqs = (pseq / p0)**(Rd/cp) |
---|
| 262 | |
---|
| 263 | ! Iterate on equatorial profile |
---|
| 264 | do iter = 1, 12 |
---|
| 265 | |
---|
| 266 | ! Calculate Exner pressure in equatorial column (p0 at surface) |
---|
| 267 | rhs(1,:) = - g / cp / thetavyz(1,:) |
---|
| 268 | do k = 1, nz |
---|
| 269 | exnereq(k) = dot_product(intz(:,k), rhs(1,:)) |
---|
| 270 | end do |
---|
| 271 | do k = 2, nz |
---|
| 272 | exnereq(k) = exnereq(k) + (exnereqs - exnereq(1)) |
---|
| 273 | end do |
---|
| 274 | exnereq(1) = exnereqs |
---|
| 275 | |
---|
| 276 | ! Calculate new pressure and temperature |
---|
| 277 | do k = 1, nz |
---|
| 278 | p = p0 * exnereq(k)**(cp/Rd) |
---|
| 279 | T = thetaeq(k) * exnereq(k) |
---|
| 280 | |
---|
| 281 | qvs = saturation_mixing_ratio(p, T) |
---|
| 282 | qveq(k) = qvs * H(k) |
---|
| 283 | |
---|
| 284 | thetavyz(1,k) = thetaeq(k) * (1.d0 + 0.61d0 * qveq(k)) |
---|
| 285 | end do |
---|
| 286 | end do |
---|
| 287 | |
---|
| 288 | !do k = 1, nz |
---|
| 289 | ! write(*,*) exnereq(k) * thetaeq(k) |
---|
| 290 | !end do |
---|
| 291 | |
---|
| 292 | ! Iterate on remainder of domain |
---|
| 293 | do iter = 1, 12 |
---|
| 294 | |
---|
| 295 | ! Compute d/dz(theta) |
---|
| 296 | do i = 1, nphi |
---|
| 297 | do k = 1, nz |
---|
| 298 | dztheta(i,k) = dot_product(ddz(:,k), thetavyz(i,:)) |
---|
| 299 | end do |
---|
| 300 | end do |
---|
| 301 | |
---|
| 302 | ! Compute rhs |
---|
| 303 | rhs = sin(2.0d0*phicoordmat)/(2.0d0*g) & |
---|
| 304 | * (ueq2 * dztheta - thetavyz * dueq2) |
---|
| 305 | |
---|
| 306 | ! Integrate |
---|
| 307 | do k = 1, nz |
---|
| 308 | do i = 1, nphi |
---|
| 309 | irhs(i,k) = dot_product(intphi(:,i), rhs(:,k)) |
---|
| 310 | end do |
---|
| 311 | end do |
---|
| 312 | |
---|
| 313 | ! Apply boundary conditions (fixed Dirichlet condition at equator) |
---|
| 314 | do i = 2, nphi |
---|
| 315 | irhs(i,:) = irhs(i,:) + (thetavyz(1,:) - irhs(1,:)) |
---|
| 316 | end do |
---|
| 317 | irhs(1,:) = thetavyz(1,:) |
---|
| 318 | |
---|
| 319 | ! Compute difference after iteration |
---|
| 320 | !err = sum(irhs - thetavyz) |
---|
| 321 | !write(*,*) iter, err |
---|
| 322 | |
---|
| 323 | ! Update iteration |
---|
| 324 | thetavyz = irhs |
---|
| 325 | end do |
---|
| 326 | |
---|
| 327 | ! Calculate pressure through remainder of domain |
---|
| 328 | rhs = - ueq2 * sin(phicoordmat) * cos(phicoordmat) / cp / thetavyz |
---|
| 329 | |
---|
| 330 | do k = 1, nz |
---|
| 331 | do i = 1, nphi |
---|
| 332 | exneryz(i,k) = dot_product(intphi(:,i), rhs(:,k)) |
---|
| 333 | end do |
---|
| 334 | do i = 2, nphi |
---|
| 335 | exneryz(i,k) = exneryz(i,k) + (exnereq(k) - exneryz(1,k)) |
---|
| 336 | end do |
---|
| 337 | |
---|
| 338 | exneryz(1,k) = exnereq(k) |
---|
| 339 | end do |
---|
| 340 | |
---|
| 341 | ! Initialization successful |
---|
| 342 | initialized = 1 |
---|
| 343 | |
---|
| 344 | END SUBROUTINE supercell_init |
---|
| 345 | |
---|
| 346 | !----------------------------------------------------------------------- |
---|
| 347 | ! Evaluate the supercell initial conditions |
---|
| 348 | !----------------------------------------------------------------------- |
---|
| 349 | SUBROUTINE supercell_test(lon,lat,p,z,zcoords,u,v,t,thetav,ps,rho,q,pert) & |
---|
| 350 | BIND(c, name = "supercell_test") |
---|
| 351 | |
---|
| 352 | IMPLICIT NONE |
---|
| 353 | |
---|
| 354 | !------------------------------------------------ |
---|
| 355 | ! Input / output parameters |
---|
| 356 | !------------------------------------------------ |
---|
| 357 | REAL(8), INTENT(IN) :: & |
---|
| 358 | lon, & ! Longitude (radians) |
---|
| 359 | lat ! Latitude (radians) |
---|
| 360 | |
---|
| 361 | REAL(8), INTENT(INOUT) :: & |
---|
| 362 | p, & ! Pressure (Pa) |
---|
| 363 | z ! Altitude (m) |
---|
| 364 | |
---|
| 365 | INTEGER, INTENT(IN) :: zcoords ! 1 if z coordinates are specified |
---|
| 366 | ! 0 if p coordinates are specified |
---|
| 367 | |
---|
| 368 | REAL(8), INTENT(OUT) :: & |
---|
| 369 | u, & ! Zonal wind (m s^-1) |
---|
| 370 | v, & ! Meridional wind (m s^-1) |
---|
| 371 | t, & ! Temperature (K) |
---|
| 372 | thetav, & ! Virtual potential Temperature (K) |
---|
| 373 | ps, & ! Surface Pressure (Pa) |
---|
| 374 | rho, & ! density (kg m^-3) |
---|
| 375 | q ! water vapor mixing ratio (kg/kg) |
---|
| 376 | |
---|
| 377 | INTEGER, INTENT(IN) :: pert ! 1 if perturbation should be included |
---|
| 378 | ! 0 if no perturbation should be included |
---|
| 379 | |
---|
| 380 | !------------------------------------------------ |
---|
| 381 | ! Local variables |
---|
| 382 | !------------------------------------------------ |
---|
| 383 | |
---|
| 384 | ! Absolute latitude |
---|
| 385 | REAL(8) :: nh_lat |
---|
| 386 | |
---|
| 387 | ! Check that we are initialized |
---|
| 388 | if (initialized .ne. 1) then |
---|
| 389 | write(*,*) 'supercell_init() has not been called' |
---|
| 390 | stop |
---|
| 391 | end if |
---|
| 392 | |
---|
| 393 | !------------------------------------------------ |
---|
| 394 | ! Begin sampling |
---|
| 395 | !------------------------------------------------ |
---|
| 396 | |
---|
| 397 | ! Northern hemisphere latitude |
---|
| 398 | if (lat .le. 0.0d0) then |
---|
| 399 | nh_lat = -lat |
---|
| 400 | else |
---|
| 401 | nh_lat = lat |
---|
| 402 | end if |
---|
| 403 | |
---|
| 404 | ! Sample surface pressure |
---|
| 405 | CALL supercell_z(lon, lat, 0.d0, ps, thetav, rho, q, pert) |
---|
| 406 | |
---|
| 407 | ! Calculate dependent variables |
---|
| 408 | if (zcoords .eq. 1) then |
---|
| 409 | CALL supercell_z(lon, lat, z, p, thetav, rho, q, pert) |
---|
| 410 | else |
---|
| 411 | CALL supercell_p(lon, lat, p, z, thetav, rho, q, pert) |
---|
| 412 | end if |
---|
| 413 | |
---|
| 414 | ! Sample the zonal velocity |
---|
| 415 | u = zonal_velocity(z, lat) |
---|
| 416 | |
---|
| 417 | ! Zero meridional velocity |
---|
| 418 | v = 0.d0 |
---|
| 419 | |
---|
| 420 | ! Temperature |
---|
| 421 | t = thetav / (1.d0 + 0.61d0 * q) * (p / p0)**(Rd/cp) |
---|
| 422 | |
---|
| 423 | END SUBROUTINE supercell_test |
---|
| 424 | |
---|
| 425 | !----------------------------------------------------------------------- |
---|
| 426 | ! Calculate pointwise pressure and temperature |
---|
| 427 | !----------------------------------------------------------------------- |
---|
| 428 | SUBROUTINE supercell_z(lon, lat, z, p, thetav, rho, q, pert) |
---|
| 429 | |
---|
| 430 | REAL(8), INTENT(IN) :: & |
---|
| 431 | lon, & ! Longitude (radians) |
---|
| 432 | lat, & ! Latitude (radians) |
---|
| 433 | z ! Altitude (m) |
---|
| 434 | |
---|
| 435 | INTEGER, INTENT(IN) :: pert ! 1 if perturbation should be included |
---|
| 436 | ! 0 if no perturbation should be included |
---|
| 437 | |
---|
| 438 | ! Evaluated variables |
---|
| 439 | REAL(8), INTENT(OUT) :: p, thetav, rho, q |
---|
| 440 | |
---|
| 441 | ! Northern hemisphere latitude |
---|
| 442 | REAL(8) :: nh_lat |
---|
| 443 | |
---|
| 444 | ! Pointwise Exner pressure |
---|
| 445 | REAL(8) :: exner |
---|
| 446 | |
---|
| 447 | ! Assembled variable values in a column |
---|
| 448 | REAL(8), DIMENSION(nz) :: varcol |
---|
| 449 | |
---|
| 450 | ! Coefficients for computing a polynomial fit in each coordinate |
---|
| 451 | REAL(8), DIMENSION(nphi) :: fitphi |
---|
| 452 | REAL(8), DIMENSION(nz) :: fitz |
---|
| 453 | |
---|
| 454 | ! Loop indices |
---|
| 455 | INTEGER(4) :: k |
---|
| 456 | |
---|
| 457 | ! Northern hemisphere latitude |
---|
| 458 | if (lat .le. 0.0d0) then |
---|
| 459 | nh_lat = -lat |
---|
| 460 | else |
---|
| 461 | nh_lat = lat |
---|
| 462 | end if |
---|
| 463 | |
---|
| 464 | ! Perform fit |
---|
| 465 | CALL lagrangian_polynomial_coeffs(nz, zcoord, fitz, z) |
---|
| 466 | CALL lagrangian_polynomial_coeffs(nphi, phicoord, fitphi, nh_lat) |
---|
| 467 | |
---|
| 468 | ! Obtain exner pressure of background state |
---|
| 469 | do k = 1, nz |
---|
| 470 | varcol(k) = dot_product(fitphi, exneryz(:,k)) |
---|
| 471 | end do |
---|
| 472 | exner = dot_product(fitz, varcol) |
---|
| 473 | p = p0 * exner**(cp/Rd) |
---|
| 474 | |
---|
| 475 | ! Sample the initialized fit at this point for theta_v |
---|
| 476 | do k = 1, nz |
---|
| 477 | varcol(k) = dot_product(fitphi, thetavyz(:,k)) |
---|
| 478 | end do |
---|
| 479 | thetav = dot_product(fitz, varcol) |
---|
| 480 | |
---|
| 481 | ! Sample water vapor mixing ratio |
---|
| 482 | q = dot_product(fitz, qveq) |
---|
| 483 | |
---|
| 484 | ! Fixed density |
---|
| 485 | rho = p / (Rd * exner * thetav) |
---|
| 486 | |
---|
| 487 | ! Modified virtual potential temperature |
---|
| 488 | if (pert .ne. 0) then |
---|
| 489 | thetav = thetav & |
---|
| 490 | + thermal_perturbation(lon, lat, z) * (1.d0 + 0.61d0 * q) |
---|
| 491 | end if |
---|
| 492 | |
---|
| 493 | ! Updated pressure |
---|
| 494 | p = p0 * (rho * Rd * thetav / p0)**(cp/(cp-Rd)) |
---|
| 495 | |
---|
| 496 | END SUBROUTINE supercell_z |
---|
| 497 | |
---|
| 498 | !----------------------------------------------------------------------- |
---|
| 499 | ! Calculate pointwise z and temperature given pressure |
---|
| 500 | !----------------------------------------------------------------------- |
---|
| 501 | SUBROUTINE supercell_p(lon, lat, p, z, thetav, rho, q, pert) |
---|
| 502 | |
---|
| 503 | REAL(8), INTENT(IN) :: & |
---|
| 504 | lon, & ! Longitude (radians) |
---|
| 505 | lat, & ! Latitude (radians) |
---|
| 506 | p ! Pressure (Pa) |
---|
| 507 | |
---|
| 508 | INTEGER, INTENT(IN) :: pert ! 1 if perturbation should be included |
---|
| 509 | ! 0 if no perturbation should be included |
---|
| 510 | |
---|
| 511 | ! Evaluated variables |
---|
| 512 | REAL(8), INTENT(OUT) :: z, thetav, rho, q |
---|
| 513 | |
---|
| 514 | ! Bounding interval and sampled values |
---|
| 515 | REAL(8) :: za, zb, zc, pa, pb, pc |
---|
| 516 | |
---|
| 517 | ! Iterate |
---|
| 518 | INTEGER(4) :: iter |
---|
| 519 | |
---|
| 520 | za = z |
---|
| 521 | zb = z2 |
---|
| 522 | |
---|
| 523 | CALL supercell_z(lon, lat, za, pa, thetav, rho, q, pert) |
---|
| 524 | CALL supercell_z(lon, lat, zb, pb, thetav, rho, q, pert) |
---|
| 525 | |
---|
| 526 | if (pa .lt. p) then |
---|
| 527 | write(*,*) 'Requested pressure out of range on bottom, adjust sample interval' |
---|
| 528 | write(*,*) pa, p |
---|
| 529 | stop |
---|
| 530 | end if |
---|
| 531 | if (pb .gt. p) then |
---|
| 532 | write(*,*) 'Requested pressure out of range on top, adjust sample interval' |
---|
| 533 | write(*,*) pb, p |
---|
| 534 | stop |
---|
| 535 | end if |
---|
| 536 | |
---|
| 537 | ! Iterate using fixed point method |
---|
| 538 | do iter = 1, 100 |
---|
| 539 | |
---|
| 540 | zc = (za * (pb - p) - zb * (pa - p)) / (pb - pa) |
---|
| 541 | |
---|
| 542 | CALL supercell_z(lon, lat, zc, pc, thetav, rho, q, pert) |
---|
| 543 | |
---|
| 544 | !write(*,*) pc |
---|
| 545 | |
---|
| 546 | if (abs((pc - p) / p) .lt. 1.d-10) then |
---|
| 547 | exit |
---|
| 548 | end if |
---|
| 549 | |
---|
| 550 | if (pc .gt. p) then |
---|
| 551 | za = zc |
---|
| 552 | pa = pc |
---|
| 553 | else |
---|
| 554 | zb = zc |
---|
| 555 | pb = pc |
---|
| 556 | end if |
---|
| 557 | end do |
---|
| 558 | |
---|
| 559 | if (iter .eq. 101) then |
---|
| 560 | write(*,*) 'Iteration failed to converge' |
---|
| 561 | stop |
---|
| 562 | end if |
---|
| 563 | |
---|
| 564 | z = zc |
---|
| 565 | |
---|
| 566 | END SUBROUTINE supercell_p |
---|
| 567 | |
---|
| 568 | !----------------------------------------------------------------------- |
---|
| 569 | ! Calculate pointwise z and temperature given pressure |
---|
| 570 | !----------------------------------------------------------------------- |
---|
| 571 | REAL(8) FUNCTION thermal_perturbation(lon, lat, z) |
---|
| 572 | |
---|
| 573 | REAL(8), INTENT(IN) :: & |
---|
| 574 | lon, & ! Longitude (radians) |
---|
| 575 | lat, & ! Latitude (radians) |
---|
| 576 | z ! Altitude (m) |
---|
| 577 | |
---|
| 578 | ! Great circle radius from the perturbation centerpoint |
---|
| 579 | REAL(8) :: gr |
---|
| 580 | |
---|
| 581 | ! Approximately spherical radius from the perturbation centerpoint |
---|
| 582 | REAL(8) :: Rtheta |
---|
| 583 | |
---|
| 584 | gr = a*acos(sin(pert_latc*deg2rad)*sin(lat) + & |
---|
| 585 | (cos(pert_latc*deg2rad)*cos(lat)*cos(lon-pert_lonc*deg2rad))) |
---|
| 586 | |
---|
| 587 | Rtheta = sqrt((gr/pert_rh)**2 + ((z - pert_zc) / pert_rz)**2) |
---|
| 588 | |
---|
| 589 | if (Rtheta .le. 1.d0) then |
---|
| 590 | thermal_perturbation = pert_dtheta * (cos(0.5d0 * pi * Rtheta))**2 |
---|
| 591 | else |
---|
| 592 | thermal_perturbation = 0.0d0 |
---|
| 593 | end if |
---|
| 594 | |
---|
| 595 | END FUNCTION thermal_perturbation |
---|
| 596 | |
---|
| 597 | !----------------------------------------------------------------------- |
---|
| 598 | ! Calculate the reference zonal velocity |
---|
| 599 | !----------------------------------------------------------------------- |
---|
| 600 | REAL(8) FUNCTION zonal_velocity(z, lat) |
---|
| 601 | |
---|
| 602 | IMPLICIT NONE |
---|
| 603 | |
---|
| 604 | REAL(8), INTENT(IN) :: z, lat |
---|
| 605 | |
---|
| 606 | if (z .le. zs - zt) then |
---|
| 607 | zonal_velocity = us * (z / zs) - uc |
---|
| 608 | elseif (abs(z - zs) .le. zt) then |
---|
| 609 | zonal_velocity = & |
---|
| 610 | (-4.0d0/5.0d0 + 3.0d0*z/zs - 5.0d0/4.0d0*(z**2)/(zs**2)) * us - uc |
---|
| 611 | else |
---|
| 612 | zonal_velocity = us - uc |
---|
| 613 | end if |
---|
| 614 | |
---|
| 615 | zonal_velocity = zonal_velocity * cos(lat) |
---|
| 616 | |
---|
| 617 | END FUNCTION zonal_velocity |
---|
| 618 | |
---|
| 619 | !----------------------------------------------------------------------- |
---|
| 620 | ! Calculate pointwise theta at the equator at the given altitude |
---|
| 621 | !----------------------------------------------------------------------- |
---|
| 622 | REAL(8) FUNCTION equator_theta(z) |
---|
| 623 | |
---|
| 624 | IMPLICIT NONE |
---|
| 625 | |
---|
| 626 | REAL(8), INTENT(IN) :: z |
---|
| 627 | |
---|
| 628 | if (z .le. z_tr) then |
---|
| 629 | equator_theta = & |
---|
| 630 | theta0 + (theta_tr - theta0) * (z / z_tr)**(1.25d0) |
---|
| 631 | else |
---|
| 632 | equator_theta = & |
---|
| 633 | theta_tr * exp(g/cp/T_tr * (z - z_tr)) |
---|
| 634 | end if |
---|
| 635 | |
---|
| 636 | END FUNCTION equator_theta |
---|
| 637 | |
---|
| 638 | !----------------------------------------------------------------------- |
---|
| 639 | ! Calculate pointwise relative humidity (in %) at the equator at the |
---|
| 640 | ! given altitude |
---|
| 641 | !----------------------------------------------------------------------- |
---|
| 642 | REAL(8) FUNCTION equator_relative_humidity(z) |
---|
| 643 | |
---|
| 644 | IMPLICIT NONE |
---|
| 645 | |
---|
| 646 | REAL(8), INTENT(IN) :: z |
---|
| 647 | |
---|
| 648 | if (z .le. z_tr) then |
---|
| 649 | equator_relative_humidity = 1.0d0 - 0.75d0 * (z / z_tr)**(1.25d0) |
---|
| 650 | else |
---|
| 651 | equator_relative_humidity = 0.25d0 |
---|
| 652 | end if |
---|
| 653 | |
---|
| 654 | END FUNCTION equator_relative_humidity |
---|
| 655 | |
---|
| 656 | !----------------------------------------------------------------------- |
---|
| 657 | ! Calculate saturation mixing ratio (in kg/kg) in terms of pressure |
---|
| 658 | ! (in Pa) and temperature (in K) |
---|
| 659 | !----------------------------------------------------------------------- |
---|
| 660 | REAL(8) FUNCTION saturation_mixing_ratio(p, T) |
---|
| 661 | |
---|
| 662 | IMPLICIT NONE |
---|
| 663 | |
---|
| 664 | REAL(8), INTENT(IN) :: & |
---|
| 665 | p, & ! Pressure in Pa |
---|
| 666 | T ! Temperature |
---|
| 667 | |
---|
| 668 | saturation_mixing_ratio = & |
---|
| 669 | 380.d0 / p * exp(17.27d0 * (T - 273.d0) / (T - 36.d0)) |
---|
| 670 | |
---|
| 671 | if (saturation_mixing_ratio > 0.014) then |
---|
| 672 | saturation_mixing_ratio = 0.014 |
---|
| 673 | end if |
---|
| 674 | |
---|
| 675 | END FUNCTION saturation_mixing_ratio |
---|
| 676 | |
---|
| 677 | !----------------------------------------------------------------------- |
---|
| 678 | ! Calculate coefficients for a Lagrangian polynomial |
---|
| 679 | !----------------------------------------------------------------------- |
---|
| 680 | SUBROUTINE lagrangian_polynomial_coeffs(npts, x, coeffs, xs) |
---|
| 681 | |
---|
| 682 | IMPLICIT NONE |
---|
| 683 | |
---|
| 684 | ! Number of points to fit |
---|
| 685 | INTEGER(4), INTENT(IN) :: npts |
---|
| 686 | |
---|
| 687 | ! Sample points to fit |
---|
| 688 | REAL(8), DIMENSION(npts), INTENT(IN) :: x |
---|
| 689 | |
---|
| 690 | ! Computed coefficients |
---|
| 691 | REAL(8), DIMENSION(npts), INTENT(OUT) :: coeffs |
---|
| 692 | |
---|
| 693 | ! Point at which sample is taken |
---|
| 694 | REAL(8), INTENT(IN) :: xs |
---|
| 695 | |
---|
| 696 | ! Loop indices |
---|
| 697 | INTEGER(4) :: i, j |
---|
| 698 | |
---|
| 699 | ! Compute the Lagrangian polynomial coefficients |
---|
| 700 | do i = 1, npts |
---|
| 701 | coeffs(i) = 1.0d0 |
---|
| 702 | do j = 1, npts |
---|
| 703 | if (i .eq. j) then |
---|
| 704 | cycle |
---|
| 705 | end if |
---|
| 706 | coeffs(i) = coeffs(i) * (xs - x(j)) / (x(i) - x(j)) |
---|
| 707 | end do |
---|
| 708 | end do |
---|
| 709 | |
---|
| 710 | END SUBROUTINE lagrangian_polynomial_coeffs |
---|
| 711 | |
---|
| 712 | !----------------------------------------------------------------------- |
---|
| 713 | ! Calculate coefficients of the derivative of a Lagrangian polynomial |
---|
| 714 | !----------------------------------------------------------------------- |
---|
| 715 | SUBROUTINE diff_lagrangian_polynomial_coeffs(npts, x, coeffs, xs) |
---|
| 716 | |
---|
| 717 | IMPLICIT NONE |
---|
| 718 | |
---|
| 719 | ! Number of points to fit |
---|
| 720 | INTEGER(4), INTENT(IN) :: npts |
---|
| 721 | |
---|
| 722 | ! Sample points to fit |
---|
| 723 | REAL(8), DIMENSION(npts), INTENT(IN) :: x |
---|
| 724 | |
---|
| 725 | ! Computed coefficients |
---|
| 726 | REAL(8), DIMENSION(npts), INTENT(OUT) :: coeffs |
---|
| 727 | |
---|
| 728 | ! Point at which sample is taken |
---|
| 729 | REAL(8), INTENT(IN) :: xs |
---|
| 730 | |
---|
| 731 | ! Loop indices |
---|
| 732 | INTEGER(4) :: i, j, imatch |
---|
| 733 | |
---|
| 734 | ! Buffer sum |
---|
| 735 | REAL(8) :: coeffsum, differential |
---|
| 736 | |
---|
| 737 | ! Check if xs is equivalent to one of the values of x |
---|
| 738 | imatch = (-1) |
---|
| 739 | do i = 1, npts |
---|
| 740 | if (abs(xs - x(i)) < 1.0d-14) then |
---|
| 741 | imatch = i |
---|
| 742 | exit |
---|
| 743 | end if |
---|
| 744 | end do |
---|
| 745 | |
---|
| 746 | ! Equivalence detected; special treatment required |
---|
| 747 | if (imatch .ne. (-1)) then |
---|
| 748 | do i = 1, npts |
---|
| 749 | coeffs(i) = 1.0d0 |
---|
| 750 | coeffsum = 0.0d0 |
---|
| 751 | |
---|
| 752 | do j = 1, npts |
---|
| 753 | if ((j .eq. i) .or. (j .eq. imatch)) then |
---|
| 754 | cycle |
---|
| 755 | end if |
---|
| 756 | |
---|
| 757 | coeffs(i) = coeffs(i) * (xs - x(j)) / (x(i) - x(j)) |
---|
| 758 | coeffsum = coeffsum + 1.0 / (xs - x(j)) |
---|
| 759 | end do |
---|
| 760 | |
---|
| 761 | if (i .ne. imatch) then |
---|
| 762 | coeffs(i) = coeffs(i) & |
---|
| 763 | * (1.0 + (xs - x(imatch)) * coeffsum) & |
---|
| 764 | / (x(i) - x(imatch)) |
---|
| 765 | else |
---|
| 766 | coeffs(i) = coeffs(i) * coeffsum |
---|
| 767 | end if |
---|
| 768 | end do |
---|
| 769 | |
---|
| 770 | ! No equivalence; simply differentiate Lagrangian fit |
---|
| 771 | else |
---|
| 772 | call lagrangian_polynomial_coeffs(npts, x, coeffs, xs) |
---|
| 773 | |
---|
| 774 | do i = 1, npts |
---|
| 775 | differential = 0.0d0 |
---|
| 776 | do j = 1, npts |
---|
| 777 | if (i .eq. j) then |
---|
| 778 | cycle |
---|
| 779 | end if |
---|
| 780 | differential = differential + 1.0 / (xs - x(j)) |
---|
| 781 | end do |
---|
| 782 | coeffs(i) = coeffs(i) * differential |
---|
| 783 | end do |
---|
| 784 | end if |
---|
| 785 | |
---|
| 786 | END SUBROUTINE |
---|
| 787 | |
---|
| 788 | END MODULE dcmip2016_supercell_mod |
---|