MODULE dcmip_initial_conditions_test_1_2_3 !======================================================================= ! ! Functions for setting up initial conditions for the dynamical core tests: ! ! 11 - Deformational Advection Test ! 12 - Hadley Cell Advection Test ! 13 - Orography Advection Test ! 20 - Impact of orography on a steady-state at rest ! 21 and 22 - Non-Hydrostatic Mountain Waves Over A Schaer-Type Mountain without and with vertical wind shear ! 31 - Non-Hydrostatic Gravity Waves ! ! Given a point specified by: ! lon longitude (radians) ! lat latitude (radians) ! p/z pressure/height ! the functions will return: ! u zonal wind (m s^-1) ! v meridional wind (m s^-1) ! w vertical velocity (m s^-1) ! t temperature (K) ! phis surface geopotential (m^2 s^-2) ! ps surface pressure (Pa) ! rho density (kj m^-3) ! q specific humidity (kg/kg) ! qi tracers (kg/kg) ! p pressure if height based (Pa) ! ! ! Authors: James Kent, Paul Ullrich, Christiane Jablonowski ! (University of Michigan, dcmip@ucar.edu) ! version 4 ! July/8/2012 ! ! Change log: (v3, June/8/2012, v4 July/8/2012, v5 July/20/2012) ! ! v2: bug fixes in the tracer initialization for height-based models ! v3: test 3-1: the density is now initialized with the unperturbed background temperature (not the perturbed temperature) ! v3: constants converted to double precision ! v4: modified tracers in test 1-1, now with cutoff altitudes. Outside of the vertical domain all tracers are set to 0 ! v4: modified cos-term in vertical velocity (now cos(2 pi t/tau)) in test 1-1, now completing one full up and down cycle ! v4: added subroutine test1_advection_orography for test 1-3 ! v4: added subroutine test2_steady_state_mountain for test 2-0 ! v4: modified parameter list for tests 2-1 and 2-2 (routine test2_schaer_mountain): addition of parameters hybrid_eta, hyam, hybm ! if the logical flag hybrid_eta is true then the pressure in pressure-based model with hybrid sigma-p (eta) coordinates is ! computed internally. In that case the hybrid coefficients hyam and hybm need to be supplied via the parameter list, ! otherwise they are not used. ! v5: Change in test 11 - change to u and w, cutoff altitudes (introduced in v4) are removed again ! v5: Change in test 12 - velocities multiplies by rho0/rho, different w0 and vertical location of the initial tracer ! ! !======================================================================= USE prec IMPLICIT NONE PRIVATE PUBLIC :: test3_gravity_wave, & ! (lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q) test2_steady_state_mountain, & test2_schaer_mountain !----------------------------------------------------------------------- ! Physical Parameters !----------------------------------------------------------------------- REAL(rstd), parameter :: & a = 6371220.0d0, & ! Earth's Radius (m) Rd = 287.0d0, & ! Ideal gas const dry air (J kg^-1 K^1) g = 9.80616d0, & ! Gravity (m s^2) cp = 1004.5d0, & ! Specific heat capacity (J kg^-1 K^1) pi = 4.d0*atan(1.d0) ! pi !----------------------------------------------------------------------- ! Additional constants !----------------------------------------------------------------------- real(rstd), parameter :: p0 = 100000.d0 ! reference pressure (Pa) CONTAINS !========================================================================================== ! TEST CASE 11 - PURE ADVECTION - 3D DEFORMATIONAL FLOW !========================================================================================== ! The 3D deformational flow test is based on the deformational flow test of Nair and Lauritzen (JCP 2010), ! with a prescribed vertical wind velocity which makes the test truly 3D. An unscaled planet (with scale parameter ! X = 1) is selected. SUBROUTINE test1_advection_deformation (lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q,q1,q2,q3,q4) IMPLICIT NONE !----------------------------------------------------------------------- ! input/output params parameters at given location !----------------------------------------------------------------------- real(rstd), intent(in) :: lon, & ! Longitude (radians) lat, & ! Latitude (radians) z ! Height (m) real(rstd), intent(inout) :: p ! Pressure (Pa) integer, intent(in) :: zcoords ! 0 or 1 see below real(rstd), intent(out) :: u, & ! Zonal wind (m s^-1) v, & ! Meridional wind (m s^-1) w, & ! Vertical Velocity (m s^-1) t, & ! Temperature (K) phis, & ! Surface Geopotential (m^2 s^-2) ps, & ! Surface Pressure (Pa) rho, & ! density (kg m^-3) q, & ! Specific Humidity (kg/kg) q1, & ! Tracer q1 (kg/kg) q2, & ! Tracer q2 (kg/kg) q3, & ! Tracer q3 (kg/kg) q4 ! Tracer q4 (kg/kg) ! if zcoords = 1, then we use z and output p ! if zcoords = 0, then we use p !----------------------------------------------------------------------- ! test case parameters !----------------------------------------------------------------------- real(rstd), parameter :: tau = 12.d0 * 86400.d0, & ! period of motion 12 days u0 = (2.d0*pi*a)/tau, & ! 2 pi a / 12 days k0 = (10.d0*a)/tau, & ! Velocity Magnitude omega0 = (23000.d0*pi)/tau, & ! Velocity Magnitude T0 = 300.d0, & ! temperature H = Rd * T0 / g, & ! scale height RR = 1.d0/2.d0, & ! horizontal half width divided by 'a' ZZ = 1000.d0, & ! vertical half width z0 = 5000.d0, & ! center point in z lambda0 = 5.d0*pi/6.d0, & ! center point in longitudes lambda1 = 7.d0*pi/6.d0, & ! center point in longitudes phi0 = 0.d0, & ! center point in latitudes phi1 = 0.d0 real(rstd) :: height ! The height of the model levels real(rstd) :: ptop ! Model top in p real(rstd) :: sin_tmp, cos_tmp, sin_tmp2, cos_tmp2 ! Calculate great circle distances real(rstd) :: d1, d2, r, r2 ! For tracer calculations real(rstd) :: s, bs ! Shape function, and parameter real(rstd) :: lonp ! Translational longitude, depends on time real(rstd) :: ud ! Divergent part of u real(rstd) :: time ! Initially set to zero seconds, needs ! to be modified when used in dycore !----------------------------------------------------------------------- ! HEIGHT AND PRESSURE !----------------------------------------------------------------------- ! Height and pressure are aligned (p = p0 exp(-z/H)) if (zcoords .eq. 1) then height = z p = p0 * exp(-z/H) else height = H * log(p0/p) endif ! Model top in p ptop = p0*exp(-12000.d0/H) !----------------------------------------------------------------------- ! THE VELOCITIES ARE TIME DEPENDENT AND THEREFORE MUST BE UPDATED ! IN THE DYNAMICAL CORE !----------------------------------------------------------------------- ! These are initial conditions hence time = 0 time = 0.d0 ! Translational longitude = longitude when time = 0 lonp = lon - 2.d0*pi*time/tau ! Shape function !******** ! change in version 5: shape function !******** bs = 0.2 s = 1.0 + exp( (ptop-p0)/(bs*ptop) ) - exp( (p-p0)/(bs*ptop)) - exp( (ptop-p)/(bs*ptop)) ! Zonal Velocity !******** ! change in version 5: ud !******** ud = (omega0*a)/(bs*ptop) * cos(lonp) * (cos(lat)**2.0) * cos(2.0*pi*time/tau) * & ( - exp( (p-p0)/(bs*ptop)) + exp( (ptop-p)/(bs*ptop)) ) u = k0*sin(lonp)*sin(lonp)*sin(2.d0*lat)*cos(pi*time/tau) + u0*cos(lat) + ud ! Meridional Velocity v = k0*sin(2.d0*lonp)*cos(lat)*cos(pi*time/tau) ! Vertical Velocity - can be changed to vertical pressure velocity by ! omega = -(g*p)/(Rd*T0)*w ! !******** ! change in version 4: cos(2.0*pi*time/tau) is now used instead of cos(pi*time/tau) !******** !******** ! change in version 5: shape function in w !******** w = -((Rd*T0)/(g*p))*omega0*sin(lonp)*cos(lat)*cos(2.0*pi*time/tau)*s !----------------------------------------------------------------------- ! TEMPERATURE IS CONSTANT 300 K !----------------------------------------------------------------------- t = T0 !----------------------------------------------------------------------- ! PHIS (surface geopotential) !----------------------------------------------------------------------- phis = 0.d0 !----------------------------------------------------------------------- ! PS (surface pressure) !----------------------------------------------------------------------- ps = p0 !----------------------------------------------------------------------- ! RHO (density) !----------------------------------------------------------------------- rho = p/(Rd*t) !----------------------------------------------------------------------- ! initialize Q, set to zero !----------------------------------------------------------------------- q = 0.d0 !----------------------------------------------------------------------- ! initialize tracers !----------------------------------------------------------------------- ! Tracer 1 - Cosine Bells ! To calculate great circle distance sin_tmp = sin(lat) * sin(phi0) cos_tmp = cos(lat) * cos(phi0) sin_tmp2 = sin(lat) * sin(phi1) cos_tmp2 = cos(lat) * cos(phi1) ! great circle distance without 'a' r = ACOS (sin_tmp + cos_tmp*cos(lon-lambda0)) r2 = ACOS (sin_tmp2 + cos_tmp2*cos(lon-lambda1)) d1 = min( 1., (r/RR)**2 + ((height-z0)/ZZ)**2 ) d2 = min( 1., (r2/RR)**2 + ((height-z0)/ZZ)**2 ) q1 = 0.5d0 * (1.d0 + cos(pi*d1)) + 0.5d0 * (1.d0 + cos(pi*d2)) ! Tracer 2 - Correlated Cosine Bells q2 = 0.9d0 - 0.8d0*q1**2 ! Tracer 3 - Slotted Ellipse ! Make the ellipse if (d1 .le. RR) then q3 = 1.d0 elseif (d2 .le. RR) then q3 = 1.d0 else q3 = 0.1d0 endif ! Put in the slot if (height .gt. z0 .and. abs(lat) .lt. 0.125d0) then q3 = 0.1d0 endif ! Tracer 4: q4 is chosen so that, in combination with the other three tracer ! fields with weight (3/10), the sum is equal to one q4 = 1.d0 - 0.3d0*(q1+q2+q3) !************ ! change in version 4: added cutoff altitudes, tracers are set to zero outside this region ! prevents tracers from being trapped near the bottom and top of the domain !************ ! use a cutoff altitude for ! tracer 2 3 and 4 ! Set them to zero outside `buffer zone' !************ ! change in version 5: change from v4 reversed, no cutoff altitudes due to continuity equation being satisfied ! commented out below !************ !if (height .gt. (z0+1.25*ZZ) .or. height .lt. (z0-1.25*ZZ)) then ! q2 = 0.0 ! q3 = 0.0 ! q4 = 0.0 !endif END SUBROUTINE test1_advection_deformation !========================================================================================== ! TEST CASE 12 - PURE ADVECTION - 3D HADLEY-LIKE FLOW !========================================================================================== SUBROUTINE test1_advection_hadley (lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q,q1) IMPLICIT NONE !----------------------------------------------------------------------- ! input/output params parameters at given location !----------------------------------------------------------------------- real(rstd), intent(in) :: lon, & ! Longitude (radians) lat, & ! Latitude (radians) z ! Height (m) real(rstd), intent(inout) :: p ! Pressure (Pa) integer, intent(in) :: zcoords ! 0 or 1 see below real(rstd), intent(out) :: u, & ! Zonal wind (m s^-1) v, & ! Meridional wind (m s^-1) w, & ! Vertical Velocity (m s^-1) t, & ! Temperature (K) phis, & ! Surface Geopotential (m^2 s^-2) ps, & ! Surface Pressure (Pa) rho, & ! density (kg m^-3) q, & ! Specific Humidity (kg/kg) q1 ! Tracer q1 (kg/kg) ! if zcoords = 1, then we use z and output p ! if zcoords = 0, then we use p !----------------------------------------------------------------------- ! test case parameters !----------------------------------------------------------------------- real(rstd), parameter :: tau = 1.d0 * 86400.d0, & ! period of motion 1 day (in s) u0 = 40.d0, & ! Zonal velocity magnitude (m/s) w0 = 0.15d0, & ! Vertical velocity magnitude (m/s), changed in v5 T0 = 300.d0, & ! temperature (K) H = Rd * T0 / g, & ! scale height K = 5.d0, & ! number of Hadley-like cells z1 = 2000.d0, & ! position of lower tracer bound (m), changed in v5 z2 = 5000.d0, & ! position of upper tracer bound (m), changed in v5 z0 = 0.5d0*(z1+z2), & ! midpoint (m) ztop = 12000.d0 ! model top (m) real(rstd) :: rho0 ! reference density at z=0 m real(rstd) :: height ! Model level heights real(rstd) :: time ! Initially set to zero seconds, needs ! to be modified when used in dycore !----------------------------------------------------------------------- ! HEIGHT AND PRESSURE !----------------------------------------------------------------------- ! Height and pressure are aligned (p = p0 exp(-z/H)) if (zcoords .eq. 1) then height = z p = p0 * exp(-z/H) else height = H * log(p0/p) endif !----------------------------------------------------------------------- ! TEMPERATURE IS CONSTANT 300 K !----------------------------------------------------------------------- t = T0 !----------------------------------------------------------------------- ! PHIS (surface geopotential) !----------------------------------------------------------------------- phis = 0.d0 !----------------------------------------------------------------------- ! PS (surface pressure) !----------------------------------------------------------------------- ps = p0 !----------------------------------------------------------------------- ! RHO (density) !----------------------------------------------------------------------- rho = p/(Rd*t) rho0 = p0/(Rd*t) !----------------------------------------------------------------------- ! THE VELOCITIES ARE TIME DEPENDENT AND THEREFORE MUST BE UPDATED ! IN THE DYNAMICAL CORE !----------------------------------------------------------------------- time = 0.d0 ! Zonal Velocity u = u0*cos(lat) ! Meridional Velocity !************ ! change in version 5: multiply v and w by rho0/rho !************ v = -(rho0/rho) * (a*w0*pi)/(K*ztop) *cos(lat)*sin(K*lat)*cos(pi*height/ztop)*cos(pi*time/tau) ! Vertical Velocity - can be changed to vertical pressure velocity by ! omega = -g*rho*w w = (rho0/rho) *(w0/K)*(-2.d0*sin(K*lat)*sin(lat) + K*cos(lat)*cos(K*lat)) & *sin(pi*height/ztop)*cos(pi*time/tau) !----------------------------------------------------------------------- ! initialize Q, set to zero !----------------------------------------------------------------------- q = 0.d0 !----------------------------------------------------------------------- ! initialize tracers !----------------------------------------------------------------------- ! Tracer 1 - Layer if (height .lt. z2 .and. height .gt. z1) then q1 = 0.5d0 * (1.d0 + cos( 2.d0*pi*(height-z0)/(z2-z1) ) ) else q1 = 0.d0 endif END SUBROUTINE test1_advection_hadley !========================================================================================== ! TEST CASE 13 - HORIZONTAL ADVECTION OF THIN CLOUD-LIKE TRACERS IN THE PRESENCE OF OROGRAPHY !========================================================================================== SUBROUTINE test1_advection_orography (lon,lat,p,z,zcoords,cfv,hybrid_eta,hyam,hybm,gc,u,v,w,t,phis,ps,rho,q,q1,q2,q3,q4) IMPLICIT NONE !----------------------------------------------------------------------- ! input/output params parameters at given location !----------------------------------------------------------------------- real(rstd), intent(in) :: lon, & ! Longitude (radians) lat, & ! Latitude (radians) z, & ! Height (m) hyam, & ! A coefficient for hybrid-eta coordinate, at model level midpoint hybm, & ! B coefficient for hybrid-eta coordinate, at model level midpoint gc ! bar{z} for Gal-Chen coordinate logical, intent(in) :: hybrid_eta ! flag to indicate whether the hybrid sigma-p (eta) coordinate is used ! if set to .true., then the pressure will be computed via the ! hybrid coefficients hyam and hybm, they need to be initialized ! if set to .false. (for pressure-based models): the pressure is already pre-computed ! and is an input value for this routine ! for height-based models: pressure will always be computed based on the height and ! hybrid_eta is not used ! Note that we only use hyam and hybm for the hybrid-eta coordiantes, and we only use ! gc for the Gal-Chen coordinates. If not required then they become dummy variables real(rstd), intent(inout) :: p ! Pressure (Pa) integer, intent(in) :: zcoords ! 0 or 1 see below integer, intent(in) :: cfv ! 0, 1 or 2 see below real(rstd), intent(out) :: u, & ! Zonal wind (m s^-1) v, & ! Meridional wind (m s^-1) w, & ! Vertical Velocity (m s^-1) t, & ! Temperature (K) phis, & ! Surface Geopotential (m^2 s^-2) ps, & ! Surface Pressure (Pa) rho, & ! density (kg m^-3) q, & ! Specific Humidity (kg/kg) q1, & ! Tracer q1 (kg/kg) q2, & ! Tracer q2 (kg/kg) q3, & ! Tracer q3 (kg/kg) q4 ! Tracer q4 (kg/kg) ! if zcoords = 1, then we use z and output p ! if zcoords = 0, then we use p ! if cfv = 0 we assume that our horizontal velocities are not coordinate following ! if cfv = 1 then our velocities follow hybrid eta coordinates and we need to specify w ! if cfv = 2 then our velocities follow Gal-Chen coordinates and we need to specify w ! In hybrid-eta coords: p = hyam p0 + hybm ps ! In Gal-Chen coords: z = zs + (gc/ztop)*(ztop - zs) ! if other orography-following coordinates are used, the w wind needs to be newly derived for them !----------------------------------------------------------------------- ! test case parameters !----------------------------------------------------------------------- real(rstd), parameter :: tau = 12.d0 * 86400.d0, & ! period of motion 12 days (s) u0 = 2.d0*pi*a/tau, & ! Velocity Magnitude (m/s) T0 = 300.d0, & ! temperature (K) H = Rd * T0 / g, & ! scale height (m) alpha = pi/6.d0, & ! rotation angle (radians), 30 degrees lambdam = 3.d0*pi/2.d0, & ! mountain longitude center point (radians) phim = 0.d0, & ! mountain latitude center point (radians) h0 = 2000.d0, & ! peak height of the mountain range (m) Rm = 3.d0*pi/4.d0, & ! mountain radius (radians) zetam = pi/16.d0, & ! mountain oscillation half-width (radians) lambdap = pi/2.d0, & ! cloud-like tracer longitude center point (radians) phip = 0.d0, & ! cloud-like tracer latitude center point (radians) Rp = pi/4.d0, & ! cloud-like tracer radius (radians) zp1 = 3050.d0, & ! midpoint of first (lowermost) tracer (m) zp2 = 5050.d0, & ! midpoint of second tracer (m) zp3 = 8200.d0, & ! midpoint of third (topmost) tracer (m) dzp1 = 1000.d0, & ! thickness of first (lowermost) tracer (m) dzp2 = 1000.d0, & ! thickness of second tracer (m) dzp3 = 400.d0, & ! thickness of third (topmost) tracer (m) ztop = 12000.d0 ! model top (m) real(rstd) :: height ! Model level heights (m) real(rstd) :: r ! Great circle distance (radians) real(rstd) :: rz ! height differences real(rstd) :: zs ! Surface elevation (m) !----------------------------------------------------------------------- ! PHIS (surface geopotential) !----------------------------------------------------------------------- r = acos( sin(phim)*sin(lat) + cos(phim)*cos(lat)*cos(lon - lambdam) ) if (r .lt. Rm) then zs = (h0/2.d0)*(1.d0+cos(pi*r/Rm))*cos(pi*r/zetam)**2.d0 else zs = 0.d0 endif phis = g*zs !----------------------------------------------------------------------- ! PS (surface pressure) !----------------------------------------------------------------------- ps = p0 * exp(-zs/H) !----------------------------------------------------------------------- ! HEIGHT AND PRESSURE !----------------------------------------------------------------------- ! Height and pressure are aligned (p = p0 exp(-z/H)) if (zcoords .eq. 1) then height = z p = p0 * exp(-z/H) else if (hybrid_eta) p = hyam*p0 + hybm*ps ! compute the pressure based on the surface pressure and hybrid coefficients height = H * log(p0/p) endif !----------------------------------------------------------------------- ! THE VELOCITIES ARE TIME INDEPENDENT !----------------------------------------------------------------------- ! Zonal Velocity u = u0*(cos(lat)*cos(alpha)+sin(lat)*cos(lon)*sin(alpha)) ! Meridional Velocity v = -u0*(sin(lon)*sin(alpha)) !----------------------------------------------------------------------- ! TEMPERATURE IS CONSTANT 300 K !----------------------------------------------------------------------- t = T0 !----------------------------------------------------------------------- ! RHO (density) !----------------------------------------------------------------------- rho = p/(Rd*t) !----------------------------------------------------------------------- ! initialize Q, set to zero !----------------------------------------------------------------------- q = 0.d0 !----------------------------------------------------------------------- ! VERTICAL VELOCITY IS TIME INDEPENDENT !----------------------------------------------------------------------- ! Vertical Velocity - can be changed to vertical pressure velocity by ! omega = -(g*p)/(Rd*T0)*w ! NOTE that if orography-following coordinates are used then the vertical ! velocity needs to be translated into the new coordinate system due to ! the variation of the height along coordinate surfaces ! See section 1.3 and the appendix of the test case document if (cfv .eq. 0) then ! if the horizontal velocities do not follow the vertical coordinate w = 0.d0 elseif (cfv .eq. 1) then ! if the horizontal velocities follow hybrid eta coordinates then ! the perceived vertical velocity is call test1_advection_orograph_hybrid_eta_velocity(w) elseif (cfv .eq. 2) then ! if the horizontal velocities follow Gal Chen coordinates then ! the perceived vertical velocity is call test1_advection_orograph_Gal_Chen_velocity(w) ! else ! compute your own vertical velocity if other orography-following ! vertical coordinate is used ! endif !----------------------------------------------------------------------- ! initialize tracers !----------------------------------------------------------------------- ! Tracer 1 - Cloud Layer r = acos( sin(phip)*sin(lat) + cos(phip)*cos(lat)*cos(lon - lambdap) ) rz = abs(height - zp1) if (rz .lt. 0.5d0*dzp1 .and. r .lt. Rp) then q1 = 0.25d0*(1.d0+cos(2.d0*pi*rz/dzp1))*(1.d0+cos(pi*r/Rp)) else q1 = 0.d0 endif rz = abs(height - zp2) if (rz .lt. 0.5d0*dzp2 .and. r .lt. Rp) then q2 = 0.25d0*(1.d0+cos(2.d0*pi*rz/dzp2))*(1.d0+cos(pi*r/Rp)) else q2 = 0.d0 endif rz = abs(height - zp3) if (rz .lt. 0.5d0*dzp3 .and. r .lt. Rp) then q3 = 1.d0 else q3 = 0.d0 endif q4 = q1 + q2 + q3 CONTAINS !----------------------------------------------------------------------- ! SUBROUTINE TO CALCULATE THE PERCEIVED VERTICAL VELOCITY ! UNDER HYBRID-ETA COORDINATES !----------------------------------------------------------------------- SUBROUTINE test1_advection_orograph_hybrid_eta_velocity(w) IMPLICIT NONE real(rstd), intent(out) :: w real(rstd) :: press, & ! hyam *p0 + hybm *ps r, & ! Great Circle Distance dzsdx, & ! Part of surface height derivative dzsdlambda, & ! Derivative of zs w.r.t lambda dzsdphi, & ! Derivative of zs w.r.t phi dzdlambda, & ! Derivative of z w.r.t lambda dzdphi, & ! Derivative of z w.r.t phi dpsdlambda, & ! Derivative of ps w.r.t lambda dpsdphi ! Derivative of ps w.r.t phi ! Calculate pressure and great circle distance to mountain center press = hyam*p0 + hybm*ps r = acos( sin(phim)*sin(lat) + cos(phim)*cos(lat)*cos(lon - lambdam) ) ! Derivatives of surface height if (r .lt. Rm) then dzsdx = -h0*pi/(2.d0*Rm)*sin(pi*r/Rm)*cos(pi*r/zetam)**2 - & (h0*pi/zetam)*(1.d0+cos(pi*r/Rm))*cos(pi*r/zetam)*sin(pi*r/zetam) else dzsdx = 0.d0 endif ! Prevent division by zero if (1.d0-cos(r)**2 .gt. 0.d0) then dzsdlambda = dzsdx * (cos(phim)*cos(lat)*sin(lon-lambdam)) & /sqrt(1.d0-cos(r)**2) dzsdphi = dzsdx * (-sin(phim)*cos(lat) + cos(phim)*sin(lat)*cos(lon-lambdam)) & /sqrt(1.d0-cos(r)**2) else dzsdlambda = 0.d0 dzsdphi = 0.d0 endif ! Derivatives of surface pressure dpsdlambda = -(g*p0/(Rd*T0))*exp(-g*zs/(Rd*T0))*dzsdlambda dpsdphi = -(g*p0/(Rd*T0))*exp(-g*zs/(Rd*T0))*dzsdphi ! Derivatives of coordinate height dzdlambda = -(Rd*T0/(g*press))*hybm*dpsdlambda dzdphi = -(Rd*T0/(g*press))*hybm*dpsdphi ! Prevent division by zero if (abs(lat) .lt. pi/2.d0) then w = - (u/(a*cos(lat)))*dzdlambda - (v/a)*dzdphi else w = 0.d0 endif END SUBROUTINE test1_advection_orograph_hybrid_eta_velocity !----------------------------------------------------------------------- ! SUBROUTINE TO CALCULATE THE PERCEIVED VERTICAL VELOCITY ! UNDER GAL-CHEN COORDINATES !----------------------------------------------------------------------- SUBROUTINE test1_advection_orograph_Gal_Chen_velocity(w) IMPLICIT NONE real(rstd), intent(out) :: w real(rstd) :: r, & ! Great Circle Distance dzsdx, & ! Part of surface height derivative dzsdlambda, & ! Derivative of zs w.r.t lambda dzsdphi, & ! Derivative of zs w.r.t phi dzdlambda, & ! Derivative of z w.r.t lambda dzdphi ! Derivative of z w.r.t phi ! Calculate great circle distance to mountain center r = acos( sin(phim)*sin(lat) + cos(phim)*cos(lat)*cos(lon - lambdam) ) ! Derivatives of surface height if (r .lt. Rm) then dzsdx = -h0*pi/(2.d0*Rm)*sin(pi*r/Rm)*cos(pi*r/zetam)**2 - & (h0*pi/zetam)*(1.d0+cos(pi*r/Rm))*cos(pi*r/zetam)*sin(pi*r/zetam) else dzsdx = 0.d0 endif ! Prevent division by zero if (1.d0-cos(r)**2 .gt. 0.d0) then dzsdlambda = dzsdx * (cos(phim)*cos(lat)*sin(lon-lambdam)) & /sqrt(1.d0-cos(r)**2) dzsdphi = dzsdx * (-sin(phim)*cos(lat) + cos(phim)*sin(lat)*cos(lon-lambdam)) & /sqrt(1.d0-cos(r)**2) else dzsdlambda = 0.d0 dzsdphi = 0.d0 endif ! Derivatives of coordinate height dzdlambda = (1.d0-gc/ztop)*dzsdlambda dzdphi = (1.d0-gc/ztop)*dzsdphi ! Prevent division by zero if (abs(lat) .lt. pi/2.d0) then w = - (u/(a*cos(lat)))*dzdlambda - (v/a)*dzdphi else w = 0.d0 endif END SUBROUTINE test1_advection_orograph_Gal_Chen_velocity END SUBROUTINE test1_advection_orography !========================================================================================== ! TEST CASE 2X - IMPACT OF OROGRAPHY ON A NON-ROTATING PLANET !========================================================================================== ! The tests in section 2-x examine the impact of 3D Schaer-like circular mountain profiles on an ! atmosphere at rest (2-0), and on flow fields with wind shear (2-1) and without vertical wind shear (2-2). ! A non-rotating planet is used for all configurations. Test 2-0 is conducted on an unscaled regular-size ! planet and primarily examines the accuracy of the pressure gradient calculation in a steady-state ! hydrostatically-balanced atmosphere at rest. This test is especially appealing for models with ! orography-following vertical coordinates. It increases the complexity of test 1-3, that investigated ! the impact of the same Schaer-type orographic profile on the accuracy of purely-horizontal passive ! tracer advection. ! ! Tests 2-1 and 2-2 increase the complexity even further since non-zero flow fields are now prescribed ! with and without vertical wind shear. In order to trigger non-hydrostatic responses the two tests are ! conducted on a reduced-size planet with reduction factor $X=500$ which makes the horizontal and ! vertical grid spacing comparable. This test clearly discriminates between non-hydrostatic and hydrostatic ! models since the expected response is in the non-hydrostatic regime. Therefore, the flow response is ! captured differently by hydrostatic models. !========================================================================= ! Test 2-0: Steady-State Atmosphere at Rest in the Presence of Orography !========================================================================= SUBROUTINE test2_steady_state_mountain (lon,lat,p,z,zcoords,hybrid_eta,hyam,hybm,u,v,w,t,phis,ps,rho,q) IMPLICIT NONE !----------------------------------------------------------------------- ! input/output params parameters at given location !----------------------------------------------------------------------- real(rstd), intent(in) :: lon, & ! Longitude (radians) lat, & ! Latitude (radians) z, & ! Height (m) hyam, & ! A coefficient for hybrid-eta coordinate, at model level midpoint hybm ! B coefficient for hybrid-eta coordinate, at model level midpoint logical, intent(in) :: hybrid_eta ! flag to indicate whether the hybrid sigma-p (eta) coordinate is used ! if set to .true., then the pressure will be computed via the ! hybrid coefficients hyam and hybm, they need to be initialized ! if set to .false. (for pressure-based models): the pressure is already pre-computed ! and is an input value for this routine ! for height-based models: pressure will always be computed based on the height and ! hybrid_eta is not used real(rstd), intent(inout) :: p ! Pressure (Pa) integer, intent(in) :: zcoords ! 0 or 1 see below real(rstd), intent(out) :: u, & ! Zonal wind (m s^-1) v, & ! Meridional wind (m s^-1) w, & ! Vertical Velocity (m s^-1) t, & ! Temperature (K) phis, & ! Surface Geopotential (m^2 s^-2) ps, & ! Surface Pressure (Pa) rho, & ! density (kg m^-3) q ! Specific Humidity (kg/kg) ! if zcoords = 1, then we use z and output p ! if zcoords = 0, then we compute or use p ! ! In hybrid-eta coords: p = hyam p0 + hybm ps ! ! The grid-point based initial data are computed in this routine. !----------------------------------------------------------------------- ! test case parameters !----------------------------------------------------------------------- real(rstd), parameter :: T0 = 300.d0, & ! temperature (K) gamma = 0.0065d0, & ! temperature lapse rate (K/m) lambdam = 3.d0*pi/2.d0, & ! mountain longitude center point (radians) phim = 0.d0, & ! mountain latitude center point (radians) h0 = 2000.d0, & ! peak height of the mountain range (m) Rm = 3.d0*pi/4.d0, & ! mountain radius (radians) zetam = pi/16.d0, & ! mountain oscillation half-width (radians) ztop = 12000.d0 ! model top (m) real(rstd) :: height ! Model level heights (m) real(rstd) :: r ! Great circle distance (radians) real(rstd) :: zs ! Surface elevation (m) real(rstd) :: exponent ! exponent: g/(Rd * gamma) real(rstd) :: exponent_rev ! reversed exponent !----------------------------------------------------------------------- ! compute exponents !----------------------------------------------------------------------- exponent = g/(Rd*gamma) exponent_rev = 1.d0/exponent !----------------------------------------------------------------------- ! PHIS (surface geopotential) !----------------------------------------------------------------------- r = acos( sin(phim)*sin(lat) + cos(phim)*cos(lat)*cos(lon - lambdam) ) if (r .lt. Rm) then zs = (h0/2.d0)*(1.d0+cos(pi*r/Rm))*cos(pi*r/zetam)**2.d0 ! mountain height else zs = 0.d0 endif phis = g*zs !----------------------------------------------------------------------- ! PS (surface pressure) !----------------------------------------------------------------------- ps = p0 * (1.d0 - gamma/T0*zs)**exponent !----------------------------------------------------------------------- ! HEIGHT AND PRESSURE !----------------------------------------------------------------------- ! Height and pressure are aligned (p = p0 * (1.d0 - gamma/T0*z)**exponent) if (zcoords .eq. 1) then height = z p = p0 * (1.d0 - gamma/T0*z)**exponent else if (hybrid_eta) p = hyam*p0 + hybm*ps ! compute the pressure based on the surface pressure and hybrid coefficients height = T0/gamma * (1.d0 - (p/p0)**exponent_rev) ! compute the height at this pressure endif !----------------------------------------------------------------------- ! THE VELOCITIES ARE ZERO (STATE AT REST) !----------------------------------------------------------------------- ! Zonal Velocity u = 0.d0 ! Meridional Velocity v = 0.d0 ! Vertical Velocity w = 0.d0 !----------------------------------------------------------------------- ! TEMPERATURE WITH CONSTANT LAPSE RATE !----------------------------------------------------------------------- t = T0 - gamma*height !----------------------------------------------------------------------- ! RHO (density) !----------------------------------------------------------------------- rho = p/(Rd*t) !----------------------------------------------------------------------- ! initialize Q, set to zero !----------------------------------------------------------------------- q = 0.d0 END SUBROUTINE test2_steady_state_mountain !===================================================================================== ! Tests 2-1 and 2-2: Non-hydrostatic Mountain Waves over a Schaer-type Mountain !===================================================================================== SUBROUTINE test2_schaer_mountain (lon,lat,p,z,zcoords,hybrid_eta,hyam,hybm,shear,u,v,w,t,phis,ps,rho,q) IMPLICIT NONE !----------------------------------------------------------------------- ! input/output params parameters at given location !----------------------------------------------------------------------- real(rstd), intent(in) :: lon, & ! Longitude (radians) lat, & ! Latitude (radians) z, & ! Height (m) hyam, & ! A coefficient for hybrid-eta coordinate, at model level midpoint hybm ! B coefficient for hybrid-eta coordinate, at model level midpoint logical, intent(in) :: hybrid_eta ! flag to indicate whether the hybrid sigma-p (eta) coordinate is used ! if set to .true., then the pressure will be computed via the ! hybrid coefficients hyam and hybm, they need to be initialized ! if set to .false. (for pressure-based models): the pressure is already pre-computed ! and is an input value for this routine ! for height-based models: pressure will always be computed based on the height and ! hybrid_eta is not used real(rstd), intent(inout) :: p ! Pressure (Pa) integer, intent(in) :: zcoords, & ! 0 or 1 see below shear ! 0 or 1 see below real(rstd), intent(out) :: u, & ! Zonal wind (m s^-1) v, & ! Meridional wind (m s^-1) w, & ! Vertical Velocity (m s^-1) t, & ! Temperature (K) phis, & ! Surface Geopotential (m^2 s^-2) ps, & ! Surface Pressure (Pa) rho, & ! density (kg m^-3) q ! Specific Humidity (kg/kg) ! if zcoords = 1, then we use z and output p ! if zcoords = 0, then we either compute or use p ! if shear = 1, then we use shear flow ! if shear = 0, then we use constant u !----------------------------------------------------------------------- ! test case parameters !----------------------------------------------------------------------- real(rstd), parameter :: X = 500.d0, & ! Reduced Earth reduction factor Om = 0.d0, & ! Rotation Rate of Earth as = a/X, & ! New Radius of small Earth ueq = 20.d0, & ! Reference Velocity Teq = 300.d0, & ! Temperature at Equator Peq = 100000.d0, & ! Reference PS at Equator ztop = 30000.d0, & ! Model Top lambdac = pi/4.d0, & ! Lon of Schar Mountain Center phic = 0.d0, & ! Lat of Schar Mountain Center h0 = 250.d0, & ! Height of Mountain d = 5000.d0, & ! Mountain Half-Width xi = 4000.d0, & ! Mountain Wavelength cs = 0.00025d0 ! Wind Shear (shear=1) real(rstd) :: height ! Model level heights real(rstd) :: sin_tmp, cos_tmp ! Calculation of great circle distance real(rstd) :: r ! Great circle distance real(rstd) :: zs ! Surface height real(rstd) :: c ! Shear !----------------------------------------------------------------------- ! PHIS (surface geopotential) !----------------------------------------------------------------------- sin_tmp = sin(lat) * sin(phic) cos_tmp = cos(lat) * cos(phic) ! great circle distance with 'a/X' r = as * ACOS (sin_tmp + cos_tmp*cos(lon-lambdac)) zs = h0 * exp(-(r**2)/(d**2))*(cos(pi*r/xi)**2) phis = g*zs !----------------------------------------------------------------------- ! SHEAR FLOW OR CONSTANT FLOW !----------------------------------------------------------------------- if (shear .eq. 1) then c = cs else c = 0.d0 endif !----------------------------------------------------------------------- ! TEMPERATURE !----------------------------------------------------------------------- t = Teq *(1.d0 - (c*ueq*ueq/(g))*(sin(lat)**2) ) !----------------------------------------------------------------------- ! PS (surface pressure) !----------------------------------------------------------------------- ps = peq*exp( -(ueq*ueq/(2.d0*Rd*Teq))*(sin(lat)**2) - phis/(Rd*t) ) !----------------------------------------------------------------------- ! HEIGHT AND PRESSURE !----------------------------------------------------------------------- if (zcoords .eq. 1) then height = z p = peq*exp( -(ueq*ueq/(2.d0*Rd*Teq))*(sin(lat)**2) - g*height/(Rd*t) ) else if (hybrid_eta) p = hyam*p0 + hybm*ps ! compute the pressure based on the surface pressure and hybrid coefficients height = (Rd*t/(g))*log(peq/p) - (t*ueq*ueq/(2.d0*Teq*g))*(sin(lat)**2) endif !----------------------------------------------------------------------- ! THE VELOCITIES !----------------------------------------------------------------------- ! Zonal Velocity u = ueq * cos(lat) * sqrt( (2.d0*Teq/(t))*c*height + t/(Teq) ) ! Meridional Velocity v = 0.d0 ! Vertical Velocity = Vertical Pressure Velocity = 0 w = 0.d0 !----------------------------------------------------------------------- ! RHO (density) !----------------------------------------------------------------------- rho = p/(Rd*t) !----------------------------------------------------------------------- ! initialize Q, set to zero !----------------------------------------------------------------------- q = 0.d0 END SUBROUTINE test2_schaer_mountain !========================================================================================== ! TEST CASE 3 - GRAVITY WAVES !========================================================================================== ! The non-hydrostatic gravity wave test examines the response of models to short time-scale wavemotion ! triggered by a localized perturbation. The formulation presented in this document is new, ! but is based on previous approaches by Skamarock et al. (JAS 1994), Tomita and Satoh (FDR 2004), and ! Jablonowski et al. (NCAR Tech Report 2008) !========== ! Test 3-1 !========== SUBROUTINE test3_gravity_wave (X,lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q) IMPLICIT NONE !----------------------------------------------------------------------- ! input/output params parameters at given location !----------------------------------------------------------------------- real(rstd), intent(in) :: lon, & ! Longitude (radians) lat, & ! Latitude (radians) X ! Reduced Earth reduction factor (DCMIP value = 125) real(rstd), intent(inout) :: p, & ! Pressure (Pa) z ! Height (m) integer, intent(in) :: zcoords ! 0 or 1 see below real(rstd), intent(out) :: u, & ! Zonal wind (m s^-1) v, & ! Meridional wind (m s^-1) w, & ! Vertical Velocity (m s^-1) t, & ! Temperature (K) phis, & ! Surface Geopotential (m^2 s^-2) ps, & ! Surface Pressure (Pa) rho, & ! density (kg m^-3) q ! Specific Humidity (kg/kg) ! if zcoords = 1, then we use z and output z ! if zcoords = 0, then we use p !----------------------------------------------------------------------- ! test case parameters !----------------------------------------------------------------------- real(rstd), parameter :: & Om = 0.d0, & ! Rotation Rate of Earth u0 = 20.d0, & ! Reference Velocity ! u0 = 0.d0, & ! FIXME : no zonal wind for NH tests Teq = 300.d0, & ! Temperature at Equator Peq = 100000.d0, & ! Reference PS at Equator ztop = 10000.d0, & ! Model Top lambdac = 2.d0*pi/3.d0, & ! Lon of Pert Center d = 5000.d0, & ! Width for Pert phic = 0.d0, & ! Lat of Pert Center delta_theta = 1.d0, & ! Max Amplitude of Pert Lz = 20000.d0, & ! Vertical Wavelength of Pert N = 0.01d0, & ! Brunt-Vaisala frequency N2 = N*N, & ! Brunt-Vaisala frequency Squared bigG = (g*g)/(N2*cp) ! Constant real(rstd) :: as ! New Radius of small Earth real(rstd) :: height ! Model level height real(rstd) :: sin_tmp, cos_tmp ! Calculation of great circle distance real(rstd) :: r, s ! Shape of perturbation real(rstd) :: TS ! Surface temperature real(rstd) :: t_mean, t_pert ! Mean and pert parts of temperature real(rstd) :: theta_pert ! Pot-temp perturbation as = a/X !----------------------------------------------------------------------- ! THE VELOCITIES !----------------------------------------------------------------------- ! Zonal Velocity u = u0 * cos(lat) ! Meridional Velocity v = 0.d0 ! Vertical Velocity = Vertical Pressure Velocity = 0 w = 0.d0 !----------------------------------------------------------------------- ! PHIS (surface geopotential) !----------------------------------------------------------------------- phis = 0.d0 !----------------------------------------------------------------------- ! SURFACE TEMPERATURE !----------------------------------------------------------------------- TS = bigG + (Teq-bigG)*exp( -(u0*N2/(4.d0*g*g))*(u0+2.d0*om*as)*(cos(2.d0*lat)-1.d0) ) !----------------------------------------------------------------------- ! PS (surface pressure) !----------------------------------------------------------------------- ps = peq*exp( (u0/(4.0*bigG*Rd))*(u0+2.0*Om*as)*(cos(2.0*lat)-1.0) ) & * (TS/Teq)**(cp/Rd) !----------------------------------------------------------------------- ! HEIGHT AND PRESSURE AND MEAN TEMPERATURE !----------------------------------------------------------------------- if (zcoords .eq. 1) then height = z p = ps*( (bigG/TS)*exp(-N2*height/g)+1.d0 - (bigG/TS) )**(cp/Rd) else height = (-g/N2)*log( (TS/bigG)*( (p/ps)**(Rd/cp) - 1.d0 ) + 1.d0 ) z = height ! modified : initialize z when pressure-based endif t_mean = bigG*(1.d0 - exp(N2*height/g))+ TS*exp(N2*height/g) !----------------------------------------------------------------------- ! rho (density), unperturbed using the background temperature t_mean !----------------------------------------------------------------------- !*********** ! change in version 3: density is now initialized with unperturbed background temperature, ! temperature perturbation is added afterwards !*********** rho = p/(Rd*t_mean) !----------------------------------------------------------------------- ! POTENTIAL TEMPERATURE PERTURBATION, ! here: converted to temperature and added to the temperature field ! models with a prognostic potential temperature field can utilize ! the potential temperature perturbation theta_pert directly and add it ! to the background theta field (not included here) !----------------------------------------------------------------------- sin_tmp = sin(lat) * sin(phic) cos_tmp = cos(lat) * cos(phic) ! great circle distance with 'a/X' r = as * ACOS (sin_tmp + cos_tmp*cos(lon-lambdac)) s = (d**2)/(d**2 + r**2) theta_pert = delta_theta*s*sin(2.d0*pi*height/Lz) t_pert = theta_pert*(p/p0)**(Rd/cp) t = t_mean + t_pert !----------------------------------------------------------------------- ! initialize Q, set to zero !----------------------------------------------------------------------- q = 0.d0 END SUBROUTINE test3_gravity_wave END MODULE dcmip_initial_conditions_test_1_2_3