! ! For details, see http://www.rap.ucar.edu/~gthompsn/workshop2008/ ! This F90 program will generate a vertical profile of initial conditions for ! Case#1's "idealized" experiment with bell-shaped mountain. Below in comments ! is the output of the program. If you use different vertical levels, you ! should enter those into the "zvco_z60" array below. To simulate all ! experiments, you should alter variable "tsl" below according to instructions ! available at the web site above. ! PROGRAM PROFILE IMPLICIT NONE INTEGER :: ke ! domain index in z-direction INTEGER :: ke1 ! KE+1 REAL, PARAMETER :: cp_d = 1005.0 ! specific heat of dry air at const. pressure REAL, PARAMETER :: r_d = 287.05 ! gas constant for dry air REAL, PARAMETER :: r_v = 461.51 ! gas constant for water vapor REAL, PARAMETER :: g = 9.80665 ! acceleration due to gravity REAL :: rdocp ! r_d / cp_d REAL :: rdv ! r_d / r_v REAL :: o_m_rdv ! 1 - r_d/r_v REAL, DIMENSION(61):: zvco_z60 REAL, DIMENSION(:), ALLOCATABLE :: zml, hhl, t, p, rh, qv, u, v ! Variables for temperature and pressure profile REAL :: n_bv, tsl, psl ! Variables for moisture profile REAL :: za, zb, zc, zzc, zes, zp, zqs ! Variables for wind profile REAL :: zub, zut, zzb, zzt ! Loop indices INTEGER :: k, i, j ! === User specifications: ! Number of vertical levels ke = 60 ! User specification for the vertical model half-levels [meters] ! e.g. IF ( ke .EQ. 60 ) THEN DATA zvco_z60 / 23588.50, 22141.60, 20764.47, 19454.77, 18210.22, & 17028.57, 15907.61, 14845.17, 13839.12, 12887.38, & 11987.90, 11138.66, 10337.69, 9583.07, 8872.91, & 8205.35, 7578.58, 6990.84, 6440.39, 5925.54, & 5444.64, 4996.08, 4578.29, 4189.73, 3828.91, & 3494.38, 3184.73, 2898.58, 2634.61, 2391.51, & 2168.03, 1962.97, 1775.14, 1603.41, 1446.70, & 1303.93, 1174.11, 1056.24, 949.41, 852.71, & 765.29, 686.33, 615.06, 550.73, 492.66, & 440.19, 392.70, 349.61, 310.38, 274.53, & 241.58, 211.13, 182.79, 156.23, 131.15, & 107.29, 84.42, 62.38, 41.03, 20.25, 0.00 / ELSE STOP ! User specific model levels ENDIF ! Dry Brunt-Vaiasala frequency, ! sea level temperature, sea level pressure n_bv = 0.011 ! [s-1] tsl = 280.15 ! [K] psl = 1.0E5 ! [Pa] ! Parameters for moisture profile za = 0.95 ! [Rel. hum. surface] zb = 0.03 ! [Rel. hum. TOA] zc = 0.0015 ! Parameter for vertical decay zzc= 6000.0 ! Parameter for vertical decay ! Parameters for wind profile zub = 15.0 ! Wind speed below zzb zut = 40.0 ! Wind speed above zzb zzb = 10000.0 ! Height beyond which windspeed increases linearly zzt = zvco_z60(1) ! Top of domain ! === Initializations: ke1 = ke + 1 rdocp = r_d/cp_d rdv = r_d / r_v o_m_rdv = 1.0 - rdv ! === Allocate some fields ALLOCATE(zml(ke), hhl(ke1), t(ke), p(ke), rh(ke), qv(ke), u(ke), v(ke)) ! === Model half-levels hhl IF (ke .EQ. 60) THEN hhl(ke1) = 0.0 DO k = 1, ke1 hhl(k) = zvco_z60(k) ENDDO ENDIF ! === Compute height of full model levels zml DO k = 1, ke zml(k) = 0.5 * ( hhl(k) + hhl(k+1) ) ENDDO ! === Compute initial profile of temperature and pressure ! === on full model levels CALL CLARK_FARLEY_1984(n_bv,tsl,psl,zml,t,p) ! === Compute moisture profile CALL FERMI_PROFILE(za,zb,zc,zzc,t,p,zml,rh,qv) ! === Compute wind profile CALL WIND_PROFILE(zub,zut,zzb,zzt,zml,u) v(:) = 0.0 ! === Write sounding to IO and to file OPEN(unit=34,file="sounding.dat",status="replace",action="write") DO k = 1, ke ! Level index, height main levels [m], temp. [K], pressure [Pa], ! rh [0-1], qv [kg/kg], u [m/s], v [m/s] WRITE(*,'(i6,tr2,f8.2,tr2,f9.2,tr2,f8.2,tr2,f10.8,tr2,f6.4,tr2,f8.2,tr2,f8.2)') & k, zml(k), p(k), t(k), qv(k), rh(k), u(k), v(k) WRITE(34,'(i6,tr2,f8.2,tr2,f9.2,tr2,f8.2,tr2,f10.8,tr2,f6.4,tr2,f8.2,tr2,f8.2)') & k, zml(k), p(k), t(k), qv(k), rh(k), u(k), v(k) ENDDO CLOSE(34,status="keep") ! === Deallocation of fields DEALLOCATE(zml, hhl, t, p, rh, qv, u, v) CONTAINS SUBROUTINE CLARK_FARLEY_1984(n_bv,tsl,psl,zml,t,p) IMPLICIT NONE REAL, INTENT(IN) :: n_bv, tsl, psl REAL, DIMENSION(ke), INTENT(IN) :: zml REAL, DIMENSION(ke), INTENT(OUT) :: t, p REAL :: ztmp, zS, zth0 INTEGER :: i, j, k zS = n_bv*n_bv/g zth0 = tsl * ( 1.E5 / psl )**(rdocp) DO k = 1, ke ztmp = 1.0 - g/(cp_d*zth0*zS) * ( 1.0 - exp(-zS*zml(k)) ) t(k) = zth0 * exp(zS*zml(k)) * ztmp p(k)= psl * ztmp**(cp_d/r_d) ENDDO END SUBROUTINE CLARK_FARLEY_1984 SUBROUTINE FERMI_PROFILE(za,zb,zc,zzc,t,p,hml,rh,qv) IMPLICIT NONE REAL, INTENT(IN) :: za, zb, zc, zzc REAL, DIMENSION(ke), INTENT(IN) :: t, p, hml REAL, DIMENSION(ke), INTENT(INOUT) :: rh, qv REAL, PARAMETER :: b1 = 610.78 REAL, PARAMETER :: b2w= 17.2693882 REAL, PARAMETER :: b3 = 273.16 REAL, PARAMETER :: b4w= 35.86 REAL :: zes, zp, zqs DO k = 1, ke rh(k) = za + (zb-za)/(exp(-zc*(hml(k)-zzc))+1.0) ENDDO DO k = 1, ke ! Saturation vapor pressure over water zes = b1*EXP( b2w*(t(k)-b3) / (t(k)-b4w) ) ! Total pressure zp = p(k) ! Sat. spec. humidity zqs= rdv*zes / ( zp - o_m_rdv * zes) ! Actual spec. humidity qv(k) = rh(k) * zqs ENDDO END SUBROUTINE FERMI_PROFILE SUBROUTINE WIND_PROFILE(zub,zut,zzb,zzt,zml,u) IMPLICIT NONE REAL, INTENT(IN) :: zub, zut, zzb, zzt REAL, DIMENSION(ke), INTENT(IN) :: zml REAL, DIMENSION(ke), INTENT(OUT) :: u DO k = 1, ke IF ( zml(k) .LE. zzb ) THEN u(k) = zub ELSEIF ( zml(k) .GE. zzt ) THEN u(k) = zut ELSE u(k) = zub + (zut-zub)/(zzt-zzb) * (zml(k)-zzb) ENDIF ENDDO END SUBROUTINE WIND_PROFILE END PROGRAM PROFILE ! ! The program above produces the following vertical profile ! K Z(m) P(Pa) T(K) qv (kg/kg) RH u(m/s) v(m/s) ! 1 22865.05 1584.19 113.69 0.00000000 0.0300 38.67 0.00 ! 2 21453.04 2372.33 125.39 0.00000000 0.0300 36.07 0.00 ! 3 20109.62 3369.50 136.33 0.00000000 0.0300 33.60 0.00 ! 4 18832.50 4587.49 146.56 0.00000000 0.0300 31.25 0.00 ! 5 17619.39 6032.94 156.13 0.00000000 0.0300 29.02 0.00 ! 6 16468.09 7707.37 165.09 0.00000000 0.0300 26.90 0.00 ! 7 15376.39 9607.48 173.46 0.00000000 0.0300 24.89 0.00 ! 8 14342.14 11725.54 181.29 0.00000002 0.0300 22.99 0.00 ! 9 13363.25 14049.80 188.61 0.00000006 0.0300 21.19 0.00 ! 10 12437.64 16565.09 195.44 0.00000015 0.0301 19.48 0.00 ! 11 11563.28 19253.37 201.83 0.00000036 0.0302 17.88 0.00 ! 12 10738.18 22094.32 207.80 0.00000075 0.0308 16.36 0.00 ! 13 9960.38 25065.89 213.37 0.00000146 0.0324 15.00 0.00 ! 14 9227.99 28144.98 218.56 0.00000288 0.0372 15.00 0.00 ! 15 8539.13 31307.80 223.41 0.00000621 0.0500 15.00 0.00 ! 16 7891.96 34530.55 227.92 0.00001523 0.0809 15.00 0.00 ! 17 7284.71 37789.77 232.12 0.00003992 0.1469 15.00 0.00 ! 18 6715.62 41062.68 236.03 0.00009940 0.2644 15.00 0.00 ! 19 6182.96 44327.63 239.67 0.00021440 0.4273 15.00 0.00 ! 20 5685.09 47564.31 243.04 0.00038722 0.5967 15.00 0.00 ! 21 5220.36 50754.03 246.17 0.00059776 0.7320 15.00 0.00 ! 22 4787.19 53879.75 249.08 0.00082418 0.8216 15.00 0.00 ! 23 4384.01 56926.33 251.76 0.00105566 0.8751 15.00 0.00 ! 24 4009.32 59880.57 254.25 0.00128945 0.9058 15.00 0.00 ! 25 3661.65 62731.19 256.55 0.00152536 0.9232 15.00 0.00 ! 26 3339.55 65468.86 258.67 0.00176319 0.9333 15.00 0.00 ! 27 3041.66 68086.04 260.62 0.00200212 0.9392 15.00 0.00 ! 28 2766.60 70577.02 262.42 0.00224079 0.9429 15.00 0.00 ! 29 2513.06 72937.88 264.07 0.00247756 0.9451 15.00 0.00 ! 30 2279.77 75166.38 265.58 0.00271071 0.9465 15.00 0.00 ! 31 2065.50 77261.48 266.97 0.00293855 0.9475 15.00 0.00 ! 32 1869.05 79223.63 268.24 0.00315954 0.9481 15.00 0.00 ! 33 1689.28 81054.60 269.39 0.00337240 0.9486 15.00 0.00 ! 34 1525.05 82756.95 270.45 0.00357602 0.9489 15.00 0.00 ! 35 1375.31 84334.38 271.41 0.00376958 0.9491 15.00 0.00 ! 36 1239.02 85791.37 272.28 0.00395253 0.9493 15.00 0.00 ! 37 1115.18 87132.91 273.07 0.00412446 0.9494 15.00 0.00 ! 38 1002.82 88364.64 273.79 0.00428524 0.9495 15.00 0.00 ! 39 901.06 89492.56 274.44 0.00443490 0.9496 15.00 0.00 ! 40 809.00 90522.95 275.03 0.00457364 0.9496 15.00 0.00 ! 41 725.81 91462.28 275.56 0.00470178 0.9497 15.00 0.00 ! 42 650.70 92317.30 276.03 0.00481979 0.9497 15.00 0.00 ! 43 582.90 93094.65 276.46 0.00492821 0.9497 15.00 0.00 ! 44 521.70 93800.82 276.85 0.00502761 0.9498 15.00 0.00 ! 45 466.42 94442.36 277.20 0.00511868 0.9498 15.00 0.00 ! 46 416.45 95025.56 277.52 0.00520210 0.9498 15.00 0.00 ! 47 371.15 95556.58 277.81 0.00527855 0.9498 15.00 0.00 ! 48 329.99 96041.29 278.07 0.00534877 0.9498 15.00 0.00 ! 49 292.46 96485.06 278.30 0.00541342 0.9498 15.00 0.00 ! 50 258.05 96893.20 278.52 0.00547316 0.9498 15.00 0.00 ! 51 226.36 97270.53 278.72 0.00552863 0.9498 15.00 0.00 ! 52 196.96 97621.50 278.91 0.00558047 0.9498 15.00 0.00 ! 53 169.51 97950.22 279.08 0.00562921 0.9499 15.00 0.00 ! 54 143.69 98260.23 279.24 0.00567535 0.9499 15.00 0.00 ! 55 119.22 98554.75 279.40 0.00571931 0.9499 15.00 0.00 ! 56 95.85 98836.63 279.55 0.00576153 0.9499 15.00 0.00 ! 57 73.40 99108.19 279.69 0.00580234 0.9499 15.00 0.00 ! 58 51.71 99371.09 279.82 0.00584194 0.9499 15.00 0.00 ! 59 30.64 99626.89 279.96 0.00588060 0.9499 15.00 0.00 ! 60 10.13 99876.58 280.09 0.00591845 0.9499 15.00 0.00