| 1 |
! $Header: $ |
| 2 |
! $Name: $ |
| 3 |
|
| 4 |
module dargan_bettsmiller_mod |
| 5 |
|
| 6 |
! modified by pog: |
| 7 |
! |
| 8 |
! - changed comments to say r is mixing ratio in capecalcnew |
| 9 |
! - removed unused avg_bl code |
| 10 |
! - removed second order in p code that was commented out |
| 11 |
! - add fatal error if LCL table lookup out of range |
| 12 |
! - changed so mixing ratio has conventional definition |
| 13 |
! - changed lcl table routine call so uses conventional mixing ratio |
| 14 |
! - didn't change approximate saturated adiabatic ascent integration or |
| 15 |
! approximate 'saturate out at constant p if r0>rs' calculation to be |
| 16 |
! consistent with conventional definition of mixing ratio |
| 17 |
! - CAPE/CIN calculation uses virtual temperature if option set |
| 18 |
|
| 19 |
|
| 20 |
!---------------------------------------------------------------------- |
| 21 |
!use fms_mod, only: file_exist, error_mesg, open_file, & |
| 22 |
! check_nml_error, mpp_pe, FATAL, & |
| 23 |
! close_file |
| 24 |
|
| 25 |
use simple_sat_vapor_pres_mod, only: escomp, descomp |
| 26 |
use gcm_params_mod, only: gcm_LEN_MBUF, gcm_SQZ_R, gcm_stdMsgUnit |
| 27 |
use constants_mod, only: HLv,HLs,Cp_air,Grav,rdgas,rvgas, & |
| 28 |
kappa |
| 29 |
|
| 30 |
implicit none |
| 31 |
private |
| 32 |
!--------------------------------------------------------------------- |
| 33 |
! ---- public interfaces ---- |
| 34 |
|
| 35 |
public dargan_bettsmiller, dargan_bettsmiller_init |
| 36 |
|
| 37 |
!----------------------------------------------------------------------- |
| 38 |
! ---- version number ---- |
| 39 |
|
| 40 |
character(len=128) :: version = '$Id: dargan_bettsmiller_mod.F90,v 1.1 2012/09/11 03:53:04 jmc Exp $' |
| 41 |
character(len=128) :: tag = '$Name: $' |
| 42 |
|
| 43 |
!----------------------------------------------------------------------- |
| 44 |
! ---- local/private data ---- |
| 45 |
|
| 46 |
logical :: do_init=.true. |
| 47 |
|
| 48 |
!----------------------------------------------------------------------- |
| 49 |
! --- namelist ---- |
| 50 |
|
| 51 |
real :: tau_bm=7200. |
| 52 |
real :: rhbm = .8 |
| 53 |
logical :: do_virtual = .false. |
| 54 |
logical :: do_shallower = .false. |
| 55 |
logical :: do_changeqref = .false. |
| 56 |
logical :: do_envsat = .false. |
| 57 |
logical :: do_taucape = .false. |
| 58 |
logical :: do_bm_shift = .false. |
| 59 |
real :: capetaubm = 900. |
| 60 |
real :: tau_min = 2400. |
| 61 |
|
| 62 |
namelist /dargan_bettsmiller_nml/ tau_bm, rhbm, & |
| 63 |
do_shallower, do_changeqref, & |
| 64 |
do_envsat, do_taucape, capetaubm, tau_min, & |
| 65 |
do_virtual, do_bm_shift |
| 66 |
|
| 67 |
!----------------------------------------------------------------------- |
| 68 |
! description of namelist variables |
| 69 |
! |
| 70 |
! tau_bm = betts-miller relaxation timescale (seconds) |
| 71 |
! |
| 72 |
! rhbm = relative humidity that you're relaxing towards |
| 73 |
! |
| 74 |
! do_shallower = do the shallow convection scheme where it chooses a smaller |
| 75 |
! depth such that precipitation is zero |
| 76 |
! |
| 77 |
! do_changeqref = do the shallow convection scheme where if changes the |
| 78 |
! profile of both q and T in order make precip zero |
| 79 |
! |
| 80 |
! do_envsat = reference profile is rhbm times saturated wrt environment |
| 81 |
! (if false, it's rhbm times parcel) |
| 82 |
! |
| 83 |
! do_taucape = scheme where taubm is proportional to CAPE**-1/2 |
| 84 |
! |
| 85 |
! capetaubm = for the above scheme, the value of CAPE for which |
| 86 |
! tau = tau_bm |
| 87 |
! |
| 88 |
! tau_min = minimum relaxation time allowed for the above scheme |
| 89 |
! |
| 90 |
! do_virtual = use virtual temperature for CAPE/CIN calculation |
| 91 |
! |
| 92 |
! do_bm_shift = always conserve enthalpy by shifting temperature |
| 93 |
!----------------------------------------------------------------------- |
| 94 |
|
| 95 |
contains |
| 96 |
|
| 97 |
!####################################################################### |
| 98 |
|
| 99 |
subroutine dargan_bettsmiller (dt, tin, qin, pfull, phalf, coldT, & |
| 100 |
rain, snow, tdel, qdel, q_ref, bmflag, & |
| 101 |
klzbs, cape, cin, t_ref,invtau_bm_t,invtau_bm_q, & |
| 102 |
capeflag, bi,bj,myIter, myThid, mask, conv) |
| 103 |
|
| 104 |
!----------------------------------------------------------------------- |
| 105 |
! |
| 106 |
! Betts-Miller Convection Scheme |
| 107 |
! |
| 108 |
!----------------------------------------------------------------------- |
| 109 |
! |
| 110 |
! input: dt time step in seconds |
| 111 |
! tin temperature at full model levels |
| 112 |
! qin specific humidity of water vapor at full |
| 113 |
! model levels |
| 114 |
! pfull pressure at full model levels |
| 115 |
! phalf pressure at half (interface) model levels |
| 116 |
! coldT should precipitation be snow at this point? |
| 117 |
! optional: |
| 118 |
! mask optional mask (0 or 1.) |
| 119 |
! conv logical flag; if true then no betts-miller |
| 120 |
! adjustment is performed at that grid-point or |
| 121 |
! model level |
| 122 |
! |
| 123 |
! output: rain liquid precipitation (kg/m2) |
| 124 |
! snow frozen precipitation (kg/m2) |
| 125 |
! tdel temperature tendency at full model levels |
| 126 |
! qdel specific humidity tendency (of water vapor) at |
| 127 |
! full model levels |
| 128 |
! bmflag flag for which routines you're calling |
| 129 |
! klzbs stored klzb values |
| 130 |
! cape convectively available potential energy |
| 131 |
! cin convective inhibition (this and the above are before the |
| 132 |
! adjustment) |
| 133 |
! invtau_bm_t temperature relaxation timescale |
| 134 |
! invtau_bm_q humidity relaxation timescale |
| 135 |
! capeflag a flag that says why cape=0 |
| 136 |
|
| 137 |
!----------------------------------------------------------------------- |
| 138 |
!--------------------- interface arguments ----------------------------- |
| 139 |
|
| 140 |
real , intent(in) , dimension(:,:,:) :: tin, qin, pfull, phalf |
| 141 |
real , intent(in) :: dt |
| 142 |
logical , intent(in) , dimension(:,:):: coldT |
| 143 |
real , intent(out), dimension(:,:) :: rain,snow, bmflag, klzbs, cape, & |
| 144 |
cin, invtau_bm_t, invtau_bm_q, capeflag |
| 145 |
real , intent(out), dimension(:,:,:) :: tdel, qdel, q_ref, t_ref |
| 146 |
integer, intent(in) :: bi,bj, myIter |
| 147 |
integer, intent(in) :: myThid |
| 148 |
real , intent(in) , dimension(:,:,:), optional :: mask |
| 149 |
logical, intent(in) , dimension(:,:,:), optional :: conv |
| 150 |
!----------------------------------------------------------------------- |
| 151 |
!---------------------- local data ------------------------------------- |
| 152 |
|
| 153 |
logical,dimension(size(tin,1),size(tin,2),size(tin,3)) :: do_adjust |
| 154 |
real,dimension(size(tin,1),size(tin,2),size(tin,3)) :: & |
| 155 |
rin, esat, qsat, desat, dqsat, pmes, pmass |
| 156 |
real,dimension(size(tin,1),size(tin,2)) :: & |
| 157 |
hlcp, precip, precip_t |
| 158 |
real,dimension(size(tin,3)) :: eref, rpc, tpc, & |
| 159 |
tpc1, rpc1 |
| 160 |
|
| 161 |
real :: & |
| 162 |
cape1, cin1, tot, deltak, deltaq, qrefint, deltaqfrac, deltaqfrac2, & |
| 163 |
ptopfrac, es, capeflag1, plzb, plcl, cape2, small |
| 164 |
integer i, j, k, ix, jx, kx, klzb, ktop, klzb2 |
| 165 |
!----------------------------------------------------------------------- |
| 166 |
! computation of precipitation by betts-miller scheme |
| 167 |
!----------------------------------------------------------------------- |
| 168 |
capeflag1 = 0. |
| 169 |
|
| 170 |
! if (do_init) call error_mesg ('dargan_bettsmiller', & |
| 171 |
! 'dargan_bettsmiller_init has not been called.', FATAL) |
| 172 |
|
| 173 |
ix=size(tin,1) |
| 174 |
jx=size(tin,2) |
| 175 |
kx=size(tin,3) |
| 176 |
small = 1.e-10 |
| 177 |
|
| 178 |
! calculate r (where r is the mixing ratio) |
| 179 |
rin = qin/(1.0 - qin) |
| 180 |
|
| 181 |
do i=1,ix |
| 182 |
do j=1,jx |
| 183 |
cape1 = 0. |
| 184 |
cin1 = 0. |
| 185 |
tot = 0. |
| 186 |
klzb=0 |
| 187 |
! the bmflag is written out to show what aspects of the bm scheme is called |
| 188 |
! bmflag = 0 is no cape, no convection |
| 189 |
! bmflag = 1 is shallow conv, the predicted precip is less than zero |
| 190 |
! bmflag = 2 is deep convection |
| 191 |
bmflag(i,j) = 0. |
| 192 |
tpc = tin(i,j,:) |
| 193 |
rpc = rin(i,j,:) |
| 194 |
! calculate cape, cin, level of zero buoyancy, and parcel properties |
| 195 |
! new code (second order in delta ln p and exact LCL calculation) |
| 196 |
call capecalcnew( kx, pfull(i,j,:), phalf(i,j,:),& |
| 197 |
cp_air, rdgas, rvgas, hlv, kappa, tin(i,j,:), & |
| 198 |
rin(i,j,:), cape1, cin1, tpc, & |
| 199 |
rpc, klzb, i,j,bi,bj,myIter,myThid) |
| 200 |
|
| 201 |
! set values for storage |
| 202 |
capeflag(i,j) = capeflag1 |
| 203 |
cape(i,j) = cape1 |
| 204 |
cin(i,j) = cin1 |
| 205 |
klzbs(i,j) = klzb |
| 206 |
if(cape1.gt.0.) then |
| 207 |
! if((tot.gt.0.).and.(cape1.gt.0.)) then |
| 208 |
bmflag(i,j) = 1. |
| 209 |
! reference temperature is just that of the parcel all the way up |
| 210 |
t_ref(i,j,:) = tpc |
| 211 |
do k=klzb,kx |
| 212 |
! sets reference spec hum to a certain relative hum (change to vapor pressure, |
| 213 |
! multiply by rhbm, then back to spec humid) |
| 214 |
if(do_envsat) then |
| 215 |
call escomp(tin(i,j,k),es) |
| 216 |
es = es*rhbm |
| 217 |
rpc(k) = mixing_ratio(es, pfull(i,j,k)) |
| 218 |
q_ref(i,j,k) = rpc(k)/(1 + rpc(k)) |
| 219 |
else |
| 220 |
eref(k) = rhbm*pfull(i,j,k)*rpc(k)/(rdgas/rvgas + rpc(k)) |
| 221 |
rpc(k) = mixing_ratio(eref(k),pfull(i,j,k)) |
| 222 |
q_ref(i,j,k) = rpc(k)/(1 + rpc(k)) |
| 223 |
endif |
| 224 |
end do |
| 225 |
! set the tendencies to zero where you don't adjust |
| 226 |
! set the reference profiles to be the original profiles (for diagnostic |
| 227 |
! purposes only -- you can think of this as what you're relaxing to in |
| 228 |
! areas above the actual convection |
| 229 |
do k=1,max(klzb-1,1) |
| 230 |
qdel(i,j,k) = 0.0 |
| 231 |
tdel(i,j,k) = 0.0 |
| 232 |
q_ref(i,j,k) = qin(i,j,k) |
| 233 |
t_ref(i,j,k) = tin(i,j,k) |
| 234 |
end do |
| 235 |
! initialize p to zero for the loop |
| 236 |
precip(i,j) = 0. |
| 237 |
precip_t(i,j) = 0. |
| 238 |
! makes t_bm prop to (CAPE)**-.5. Gives a relaxation time of tau_bm when |
| 239 |
! CAPE = sqrt(capetaubm) |
| 240 |
if(do_taucape) then |
| 241 |
tau_bm = sqrt(capetaubm)*tau_bm/sqrt(cape1) |
| 242 |
if(tau_bm.lt.tau_min) tau_bm = tau_min |
| 243 |
endif |
| 244 |
do k=klzb, kx |
| 245 |
! relax to reference profiles |
| 246 |
tdel(i,j,k) = - (tin(i,j,k) - t_ref(i,j,k))/tau_bm*dt |
| 247 |
qdel(i,j,k) = - (qin(i,j,k) - q_ref(i,j,k))/tau_bm*dt |
| 248 |
! Precipitation can be calculated already, based on the change in q on the |
| 249 |
! way up (this isn't altered in the energy conservation scheme). |
| 250 |
precip(i,j) = precip(i,j) - qdel(i,j,k)*(phalf(i,j,k+1)- & |
| 251 |
phalf(i,j,k))/grav |
| 252 |
precip_t(i,j)= precip_t(i,j) + cp_air/(hlv+small)*tdel(i,j,k)* & |
| 253 |
(phalf(i,j,k+1)-phalf(i,j,k))/grav |
| 254 |
end do |
| 255 |
if ((precip(i,j).gt.0.).and.(precip_t(i,j).gt.0.)) then |
| 256 |
! If precip > 0, then correct energy. |
| 257 |
bmflag(i,j) = 2. |
| 258 |
! not simple scheme: shift the reference profile of temperature |
| 259 |
! deltak is the energy correction that you make to the temperature reference |
| 260 |
! profile |
| 261 |
if(precip(i,j).gt.precip_t(i,j) .and. (.not. do_bm_shift)) then |
| 262 |
! if the q precip is greater, then lengthen the relaxation timescale on q to |
| 263 |
! conserve energy. qdel is therefore changed. |
| 264 |
invtau_bm_q(i,j) = precip_t(i,j)/precip(i,j)/tau_bm |
| 265 |
qdel(i,j,klzb:kx) = tau_bm*invtau_bm_q(i,j)* & |
| 266 |
qdel(i,j,klzb:kx) |
| 267 |
precip(i,j) = precip_t(i,j) |
| 268 |
invtau_bm_t(i,j) = 1./tau_bm |
| 269 |
else |
| 270 |
|
| 271 |
deltak = 0. |
| 272 |
do k=klzb, kx |
| 273 |
! Calculate the integrated difference in energy change within each level. |
| 274 |
deltak = deltak - (tdel(i,j,k) + hlv/cp_air*& |
| 275 |
qdel(i,j,k))* & |
| 276 |
(phalf(i,j,k+1) - phalf(i,j,k)) |
| 277 |
end do |
| 278 |
! Divide by total pressure. |
| 279 |
deltak = deltak/(phalf(i,j,kx+1) - phalf(i,j,klzb)) |
| 280 |
! Adjust the reference profile (uniformly with height), and correspondingly |
| 281 |
! the temperature change. |
| 282 |
t_ref(i,j,klzb:kx) = t_ref(i,j,klzb:kx)+ & |
| 283 |
deltak*tau_bm/dt |
| 284 |
tdel(i,j,klzb:kx) = tdel(i,j,klzb:kx) + deltak |
| 285 |
endif |
| 286 |
else if(precip_t(i,j).gt.0.) then |
| 287 |
! If precip < 0, then do the shallow conv routine. |
| 288 |
! First option: do_shallower = true |
| 289 |
! This chooses the depth of convection based on choosing the height that |
| 290 |
! it can make precip zero, i.e., subtract off heights until that precip |
| 291 |
! becomes positive. |
| 292 |
if (do_shallower) then |
| 293 |
! ktop is the new top of convection. set this initially to klzb. |
| 294 |
ktop = klzb |
| 295 |
! Work your way down until precip is positive again. |
| 296 |
do while ((precip(i,j).lt.0.).and.(ktop.le.kx)) |
| 297 |
precip(i,j) = precip(i,j) - qdel(i,j,ktop)* & |
| 298 |
(phalf(i,j,ktop) - phalf(i,j,ktop+1))/grav |
| 299 |
ktop = ktop + 1 |
| 300 |
end do |
| 301 |
! since there will be an overshoot (precip is going to be greater than zero |
| 302 |
! once we finish this), the actual new top of convection is somewhere between |
| 303 |
! the current ktop, and one level above this. set ktop to the level above. |
| 304 |
ktop = ktop - 1 |
| 305 |
! Adjust the tendencies in the places above back to zero, and the reference |
| 306 |
! profiles back to the original t,q. |
| 307 |
if (ktop.gt.klzb) then |
| 308 |
qdel(i,j,klzb:ktop-1) = 0. |
| 309 |
q_ref(i,j,klzb:ktop-1) = qin(i,j,klzb:ktop-1) |
| 310 |
tdel(i,j,klzb:ktop-1) = 0. |
| 311 |
t_ref(i,j,klzb:ktop-1) = tin(i,j,klzb:ktop-1) |
| 312 |
end if |
| 313 |
! Then make the change only a fraction of the new top layer so the precip is |
| 314 |
! identically zero. |
| 315 |
! Calculate the fractional penetration of convection through that top layer. |
| 316 |
! This is the amount necessary to make precip identically zero. |
| 317 |
if (precip(i,j).gt.0.) then |
| 318 |
ptopfrac = precip(i,j)/(qdel(i,j,ktop)* & |
| 319 |
(phalf(i,j,ktop+1) - phalf(i,j,ktop)))*grav |
| 320 |
! Reduce qdel in the top layer by this fraction. |
| 321 |
qdel(i,j,ktop) = ptopfrac*qdel(i,j,ktop) |
| 322 |
! Set precip to zero |
| 323 |
precip(i,j) = 0. |
| 324 |
! Now change the reference temperature in such a way to make the net |
| 325 |
! heating zero. |
| 326 |
|
| 327 |
!! modification: pog |
| 328 |
!! Reduce tdel in the top layer |
| 329 |
tdel(i,j,ktop) = ptopfrac*tdel(i,j,ktop) |
| 330 |
!! end modification: pog |
| 331 |
|
| 332 |
deltak = 0. |
| 333 |
if (ktop.lt.kx) then |
| 334 |
!! modification: pog |
| 335 |
!! include the full mass of the top layer when calculating delta_k |
| 336 |
!! also use the entire mass of the column when normalizing |
| 337 |
! |
| 338 |
! Integrate temperature tendency |
| 339 |
do k=ktop,kx |
| 340 |
deltak = deltak + tdel(i,j,k)* & |
| 341 |
(phalf(i,j,k) - phalf(i,j,k+1)) |
| 342 |
end do |
| 343 |
|
| 344 |
! Normalize by the pressure difference. |
| 345 |
deltak = deltak/(phalf(i,j,kx+1) - phalf(i,j,ktop)) |
| 346 |
|
| 347 |
!! end modification: pog |
| 348 |
|
| 349 |
! Subtract this value uniformly from tdel, and make the according change to |
| 350 |
! t_ref. |
| 351 |
do k=ktop,kx |
| 352 |
tdel(i,j,k) = tdel(i,j,k) + deltak |
| 353 |
t_ref(i,j,k) = t_ref(i,j,k) + deltak*tau_bm/dt |
| 354 |
end do |
| 355 |
end if |
| 356 |
else |
| 357 |
precip(i,j) = 0. |
| 358 |
qdel(i,j,kx) = 0. |
| 359 |
q_ref(i,j,kx) = qin(i,j,kx) |
| 360 |
tdel(i,j,kx) = 0. |
| 361 |
t_ref(i,j,kx) = tin(i,j,kx) |
| 362 |
invtau_bm_t(i,j) = 0. |
| 363 |
invtau_bm_q(i,j) = 0. |
| 364 |
end if |
| 365 |
else if(do_changeqref) then |
| 366 |
! Change the reference profile of q by a certain fraction so that precip is |
| 367 |
! zero. This involves calculating the total integrated q_ref dp (this is the |
| 368 |
! quantity intqref), as well as the necessary change in q_ref (this is the |
| 369 |
! quantity deltaq). Then the fractional change in q_ref at each level (the |
| 370 |
! quantity deltaqfrac) is 1-deltaq/intqref. (have to multiply q_ref by |
| 371 |
! 1-deltaq/intqref at every level) Then the change in qdel is |
| 372 |
! -deltaq/intqref*q_ref*dt/tau_bm. |
| 373 |
! Change the reference profile of T by a uniform amount so that precip is zero. |
| 374 |
deltak = 0. |
| 375 |
deltaq = 0. |
| 376 |
qrefint = 0. |
| 377 |
do k=klzb,kx |
| 378 |
! deltaq = a positive quantity (since int qdel is positive). It's how |
| 379 |
! much q_ref must be changed by, in an integrated sense. The requisite |
| 380 |
! change in qdel is this without the factors of tau_bm and dt. |
| 381 |
deltaq = deltaq - qdel(i,j,k)*tau_bm/dt* & |
| 382 |
(phalf(i,j,k) - phalf(i,j,k+1)) |
| 383 |
! deltak = the amount tdel needs to be changed |
| 384 |
deltak = deltak + tdel(i,j,k)* & |
| 385 |
(phalf(i,j,k) - phalf(i,j,k+1)) |
| 386 |
! qrefint = integrated value of qref |
| 387 |
qrefint = qrefint - q_ref(i,j,k)* & |
| 388 |
(phalf(i,j,k) - phalf(i,j,k+1)) |
| 389 |
end do |
| 390 |
! Normalize deltak by total pressure. |
| 391 |
deltak = deltak /(phalf(i,j,kx+1) - phalf(i,j,klzb)) |
| 392 |
! multiplying factor for q_ref is 1 + the ratio |
| 393 |
deltaqfrac = 1. - deltaq/qrefint |
| 394 |
! multiplying factor for qdel adds dt/tau_bm |
| 395 |
deltaqfrac2 = - deltaq/qrefint*dt/tau_bm |
| 396 |
precip(i,j) = 0.0 |
| 397 |
do k=klzb,kx |
| 398 |
qdel(i,j,k) = qdel(i,j,k) + deltaqfrac2*q_ref(i,j,k) |
| 399 |
q_ref(i,j,k) = deltaqfrac*q_ref(i,j,k) |
| 400 |
tdel(i,j,k) = tdel(i,j,k) + deltak |
| 401 |
t_ref(i,j,k) = t_ref(i,j,k) + deltak*tau_bm/dt |
| 402 |
end do |
| 403 |
else |
| 404 |
precip(i,j) = 0. |
| 405 |
tdel(i,j,:) = 0. |
| 406 |
qdel(i,j,:) = 0. |
| 407 |
invtau_bm_t(i,j) = 0. |
| 408 |
invtau_bm_q(i,j) = 0. |
| 409 |
end if |
| 410 |
! for cases where virtual temp predicts CAPE but precip_t < 0. |
| 411 |
! - also note cape and precip_t are different because cape |
| 412 |
! involves integration wrt dlogp and precip_t wrt dp |
| 413 |
else |
| 414 |
tdel(i,j,:) = 0.0 |
| 415 |
qdel(i,j,:) = 0.0 |
| 416 |
precip(i,j) = 0.0 |
| 417 |
q_ref(i,j,:) = qin(i,j,:) |
| 418 |
t_ref(i,j,:) = tin(i,j,:) |
| 419 |
invtau_bm_t(i,j) = 0. |
| 420 |
invtau_bm_q(i,j) = 0. |
| 421 |
end if |
| 422 |
! if no CAPE, set tendencies to zero. |
| 423 |
else |
| 424 |
tdel(i,j,:) = 0.0 |
| 425 |
qdel(i,j,:) = 0.0 |
| 426 |
precip(i,j) = 0.0 |
| 427 |
q_ref(i,j,:) = qin(i,j,:) |
| 428 |
t_ref(i,j,:) = tin(i,j,:) |
| 429 |
invtau_bm_t(i,j) = 0. |
| 430 |
invtau_bm_q(i,j) = 0. |
| 431 |
end if |
| 432 |
end do |
| 433 |
end do |
| 434 |
|
| 435 |
rain = precip |
| 436 |
snow = 0. |
| 437 |
|
| 438 |
|
| 439 |
end subroutine dargan_bettsmiller |
| 440 |
|
| 441 |
!####################################################################### |
| 442 |
|
| 443 |
!all new cape calculation. |
| 444 |
|
| 445 |
subroutine capecalcnew(kx,p,phalf,cp_air,rdgas,rvgas,hlv,kappa,tin,rin,& |
| 446 |
cape,cin,tp,rp,klzb, i,j,bi,bj,myIter,myThid) |
| 447 |
|
| 448 |
! |
| 449 |
! Input: |
| 450 |
! |
| 451 |
! kx number of levels |
| 452 |
! p pressure (index 1 refers to TOA, index kx refers to surface) |
| 453 |
! phalf pressure at half levels |
| 454 |
! cp_air specific heat of dry air |
| 455 |
! rdgas gas constant for dry air |
| 456 |
! rvgas gas constant for water vapor |
| 457 |
! hlv latent heat of vaporization |
| 458 |
! kappa the constant kappa |
| 459 |
! tin temperature of the environment |
| 460 |
! rin mixing ratio of the environment |
| 461 |
! |
| 462 |
! Output: |
| 463 |
! cape Convective available potential energy |
| 464 |
! cin Convective inhibition (if there's no LFC, then this is set |
| 465 |
! to zero) |
| 466 |
! tp Parcel temperature (set to the environmental temperature |
| 467 |
! where no adjustment) |
| 468 |
! rp Parcel mixing ratio (set to the environmental humidity |
| 469 |
! where no adjustment, and set to the saturation humidity at |
| 470 |
! the parcel temperature below the LCL) |
| 471 |
! klzb Level of zero buoyancy |
| 472 |
! |
| 473 |
! Algorithm: |
| 474 |
! Start with surface parcel. |
| 475 |
! Calculate the lifting condensation level (uses an analytic formula and a |
| 476 |
! lookup table). |
| 477 |
! Calculate parcel ascent up to LZB. |
| 478 |
! Calculate CAPE and CIN. |
| 479 |
implicit none |
| 480 |
integer, intent(in) :: kx |
| 481 |
real, intent(in), dimension(:) :: p, phalf, tin, rin |
| 482 |
real, intent(in) :: rdgas, rvgas, hlv, kappa, cp_air |
| 483 |
integer, intent(out) :: klzb |
| 484 |
real, intent(out), dimension(:) :: tp, rp |
| 485 |
real, intent(out) :: cape, cin |
| 486 |
integer, intent(in) :: i,j,bi,bj,myIter |
| 487 |
integer, intent(in) :: myThid |
| 488 |
integer :: k, klcl, klfc, klcl2 |
| 489 |
logical :: nocape |
| 490 |
real, dimension(kx) :: theta, tin_virtual |
| 491 |
real :: t0, r0, es, rs, theta0, pstar, value, tlcl, & |
| 492 |
a, b, dtdlnp, d2tdlnp2, thetam, rm, tlcl2, & |
| 493 |
plcl2, plcl, plzb, small |
| 494 |
|
| 495 |
pstar = 1.e5 |
| 496 |
! so we can run dry limit (one expression involves 1/hlv) |
| 497 |
small = 1.e-10 |
| 498 |
|
| 499 |
nocape = .true. |
| 500 |
cape = 0. |
| 501 |
cin = 0. |
| 502 |
plcl = 0. |
| 503 |
plzb = 0. |
| 504 |
klfc = 0 |
| 505 |
klcl = 0 |
| 506 |
klzb = 0 |
| 507 |
tp(1:kx) = tin(1:kx) |
| 508 |
rp(1:kx) = rin(1:kx) |
| 509 |
|
| 510 |
! calculate the virtual temperature |
| 511 |
do k=1,kx |
| 512 |
tin_virtual(k) = virtual_temp(tin(k), rin(k)) |
| 513 |
enddo |
| 514 |
|
| 515 |
! start with surface parcel |
| 516 |
t0 = tin(kx) |
| 517 |
r0 = rin(kx) |
| 518 |
! calculate the lifting condensation level by the following: |
| 519 |
! are you saturated to begin with? |
| 520 |
call escomp(t0,es) |
| 521 |
rs = mixing_ratio(es, p(kx)) |
| 522 |
if (r0.ge.rs) then |
| 523 |
! if you¹re already saturated, set lcl to be the surface value. |
| 524 |
plcl = p(kx) |
| 525 |
! the first level where you¹re completely saturated. |
| 526 |
klcl = kx |
| 527 |
! saturate out to get the parcel temp and humidity at this level |
| 528 |
! first order (in delta T) accurate expression for change in temp |
| 529 |
! note this is assumed to happen at constant pressure, also the resulting |
| 530 |
! saturation mixing ratio is different from rs and this is accounted for in |
| 531 |
! the formula, but the formula is approximate and should be replaced |
| 532 |
tp(kx) = t0 + (r0 - rs)/(cp_air/(hlv+small) + hlv*rs/rvgas/t0**2.) |
| 533 |
call escomp(tp(kx),es) |
| 534 |
rp(kx) = mixing_ratio(es, p(kx)) |
| 535 |
else |
| 536 |
! if not saturated to begin with, use the analytic expression to calculate the |
| 537 |
! exact pressure and temperature where you¹re saturated. |
| 538 |
theta0 = tin(kx)*(pstar/p(kx))**kappa |
| 539 |
! the expression that we utilize is log(r/(Rd/Rv+r)/theta**(1/kappa)*pstar) = |
| 540 |
! log(es/T**(1/kappa)) |
| 541 |
! The right hand side of this is only a function of temperature, therefore |
| 542 |
! this is put into a lookup table to solve for temperature. |
| 543 |
if (r0.gt.0.) then |
| 544 |
value = log(theta0**(-1/kappa)*pstar*r0/(rdgas/rvgas+r0)) |
| 545 |
! call lcltabl(value,tlcl,myThid) |
| 546 |
call lcltabl(value,tlcl,i,j,bi,bj,myIter,myThid) |
| 547 |
plcl = pstar*(tlcl/theta0)**(1/kappa) |
| 548 |
! just in case plcl is very high up |
| 549 |
if (plcl.lt.p(1)) then |
| 550 |
plcl = p(1) |
| 551 |
tlcl = theta0*(plcl/pstar)**kappa |
| 552 |
write (*,*) 'hi lcl' |
| 553 |
end if |
| 554 |
k = kx |
| 555 |
else |
| 556 |
! if the parcel mixing ratio is zero or negative, set lcl to top level |
| 557 |
plcl = p(1) |
| 558 |
tlcl = theta0*(plcl/pstar)**kappa |
| 559 |
! write (*,*) 'zero r0', r0 |
| 560 |
go to 11 |
| 561 |
end if |
| 562 |
! calculate the parcel temperature (adiabatic ascent) below the LCL. |
| 563 |
! the mixing ratio stays the same |
| 564 |
do while (p(k).gt.plcl) |
| 565 |
tp(k) = theta0*(p(k)/pstar)**kappa |
| 566 |
call escomp(tp(k),es) |
| 567 |
! note rp is not the actual parcel mixing ratio here, but rather will be used |
| 568 |
! to calculate the reference moisture profile using rhbm |
| 569 |
rp(k) = mixing_ratio(es,p(k)) |
| 570 |
! this definition of CIN contains everything below the LCL |
| 571 |
! we use the actual parcel mixing ratio for the virtual temperature effect |
| 572 |
cin = cin + rdgas*(tin_virtual(k)-virtual_temp(tp(k),r0))*log(phalf(k+1)/phalf(k)) |
| 573 |
k = k-1 |
| 574 |
end do |
| 575 |
! first level where you¹re saturated at the level |
| 576 |
klcl = k |
| 577 |
! do a saturated ascent to get the parcel temp at the LCL. |
| 578 |
! use your 2nd order equation up to the pressure above. |
| 579 |
! moist adaibat derivatives: (use the lcl values for temp, humid, and |
| 580 |
! pressure) |
| 581 |
! note moist adiabat derivatives are approximate and should be replaced |
| 582 |
! with more accurate expressions (see e.g. Holton pg 503 for derivation) |
| 583 |
a = kappa*tlcl + hlv/cp_air*r0 |
| 584 |
b = hlv**2.*r0/cp_air/rvgas/tlcl**2. |
| 585 |
dtdlnp = a/(1. + b) |
| 586 |
! first order in p |
| 587 |
! tp(klcl) = tlcl + dtdlnp*log(p(klcl)/plcl) |
| 588 |
! second order in p (RK2) |
| 589 |
! first get temp halfway up |
| 590 |
tp(klcl) = tlcl + dtdlnp*log(p(klcl)/plcl)/2. |
| 591 |
if ((tp(klcl).lt.173.16).and.nocape) go to 11 |
| 592 |
call escomp(tp(klcl),es) |
| 593 |
rp(klcl) = mixing_ratio(es,(p(klcl) + plcl)/2) |
| 594 |
a = kappa*tp(klcl) + hlv/cp_air*rp(klcl) |
| 595 |
b = hlv**2./cp_air/rvgas*rp(klcl)/tp(klcl)**2. |
| 596 |
dtdlnp = a/(1. + b) |
| 597 |
! second half of RK2 |
| 598 |
tp(klcl) = tlcl + dtdlnp*log(p(klcl)/plcl) |
| 599 |
if ((tp(klcl).lt.173.16).and.nocape) go to 11 |
| 600 |
call escomp(tp(klcl),es) |
| 601 |
rp(klcl) = mixing_ratio(es,p(klcl)) |
| 602 |
! write (*,*) 'tp, rp klcl:kx, new', tp(klcl:kx), rp(klcl:kx) |
| 603 |
! CAPE/CIN stuff |
| 604 |
if ((virtual_temp(tp(klcl),rp(klcl)).lt.tin_virtual(klcl)).and.nocape) then |
| 605 |
! if you¹re not yet buoyant, then add to the CIN and continue |
| 606 |
cin = cin + rdgas*(tin_virtual(klcl) - & |
| 607 |
virtual_temp(tp(klcl),rp(klcl)))*log(phalf(klcl+1)/phalf(klcl)) |
| 608 |
else |
| 609 |
! if you¹re buoyant, then add to cape |
| 610 |
cape = cape + rdgas*(virtual_temp(tp(klcl),rp(klcl)) - & |
| 611 |
tin_virtual(klcl))*log(phalf(klcl+1)/phalf(klcl)) |
| 612 |
! if it¹s the first time buoyant, then set the level of free convection to k |
| 613 |
if (nocape) then |
| 614 |
nocape = .false. |
| 615 |
klfc = klcl |
| 616 |
endif |
| 617 |
end if |
| 618 |
end if |
| 619 |
! then average the properties over the boundary layer if so desired. to give |
| 620 |
! a new "parcel". this may not be saturated at the LCL, so make sure you get |
| 621 |
! to a level where it is before moist adiabatic ascent! |
| 622 |
! |
| 623 |
! then, start at the LCL, and do moist adiabatic ascent by the first order |
| 624 |
! scheme -- 2nd order as well |
| 625 |
do k=klcl-1,1,-1 |
| 626 |
! note moist adiabat derivatives are approximate and should be replaced |
| 627 |
! with more accurate expressions (see e.g. Holton pg 503 for derivation) |
| 628 |
a = kappa*tp(k+1) + hlv/cp_air*rp(k+1) |
| 629 |
b = hlv**2./cp_air/rvgas*rp(k+1)/tp(k+1)**2. |
| 630 |
dtdlnp = a/(1. + b) |
| 631 |
! first order in p |
| 632 |
! tp(k) = tp(k+1) + dtdlnp*log(p(k)/p(k+1)) |
| 633 |
! second order in p (RK2) |
| 634 |
! first get temp halfway up |
| 635 |
tp(k) = tp(k+1) + dtdlnp*log(p(k)/p(k+1))/2. |
| 636 |
if ((tp(k).lt.173.16).and.nocape) go to 11 |
| 637 |
call escomp(tp(k),es) |
| 638 |
rp(k) = mixing_ratio(es,(p(k) + p(k+1))/2) |
| 639 |
a = kappa*tp(k) + hlv/cp_air*rp(k) |
| 640 |
b = hlv**2./cp_air/rvgas*rp(k)/tp(k)**2. |
| 641 |
dtdlnp = a/(1. + b) |
| 642 |
! second half of RK2 |
| 643 |
tp(k) = tp(k+1) + dtdlnp*log(p(k)/p(k+1)) |
| 644 |
! if you're below the lookup table value, just presume that there's no way |
| 645 |
! you could have cape and call it quits |
| 646 |
if ((tp(k).lt.173.16).and.nocape) go to 11 |
| 647 |
call escomp(tp(k),es) |
| 648 |
rp(k) = mixing_ratio(es,p(k)) |
| 649 |
if ((virtual_temp(tp(k),rp(k)).lt.tin_virtual(k)).and.nocape) then |
| 650 |
! if you¹re not yet buoyant, then add to the CIN and continue |
| 651 |
cin = cin + rdgas*(tin_virtual(k)-virtual_temp(tp(k),rp(k)))*log(phalf(k+1)/phalf(k)) |
| 652 |
elseif((virtual_temp(tp(k),rp(k)).lt.tin_virtual(k)).and.(.not.nocape)) then |
| 653 |
! if you have CAPE, and it¹s your first time being negatively buoyant, |
| 654 |
! then set the level of zero buoyancy to k+1, and stop the moist ascent |
| 655 |
klzb = k+1 |
| 656 |
go to 11 |
| 657 |
else |
| 658 |
! if you¹re buoyant, then add to cape |
| 659 |
cape = cape + rdgas*(virtual_temp(tp(k),rp(k))-tin_virtual(k))*log(phalf(k+1)/phalf(k)) |
| 660 |
! if it¹s the first time buoyant, then set the level of free convection to k |
| 661 |
if (nocape) then |
| 662 |
nocape = .false. |
| 663 |
klfc = k |
| 664 |
endif |
| 665 |
end if |
| 666 |
end do |
| 667 |
11 if(nocape) then |
| 668 |
! this is if you made it through without having a LZB |
| 669 |
! set LZB to be the top level. |
| 670 |
plzb = p(1) |
| 671 |
klzb = 0 |
| 672 |
klfc = 0 |
| 673 |
cin = 0. |
| 674 |
tp(1:kx) = tin(1:kx) |
| 675 |
rp(1:kx) = rin(1:kx) |
| 676 |
end if |
| 677 |
! write (*,*) 'plcl, klcl, tlcl, r0 new', plcl, klcl, tlcl, r0 |
| 678 |
! write (*,*) 'tp, rp new', tp, rp |
| 679 |
! write (*,*) 'tp, new', tp |
| 680 |
! write (*,*) 'tin new', tin |
| 681 |
! write (*,*) 'klcl, klfc, klzb new', klcl, klfc, klzb |
| 682 |
end subroutine capecalcnew |
| 683 |
|
| 684 |
! lookup table for the analytic evaluation of LCL |
| 685 |
! subroutine lcltabl(value,tlcl,myThid) |
| 686 |
subroutine lcltabl(value,tlcl,i,j,bi,bj,myIter,myThid) |
| 687 |
! |
| 688 |
! Table of values used to compute the temperature of the lifting condensation |
| 689 |
! level. |
| 690 |
! |
| 691 |
! the expression that we utilize is log(r/(Rd/Rv+r)/theta**(1/kappa)*pstar) = |
| 692 |
! log(es/T**(1/kappa)) |
| 693 |
! |
| 694 |
! Gives the values of the temperature for the following range: |
| 695 |
! starts with -23, is uniformly distributed up to -10.4. There are a |
| 696 |
! total of 127 values, and the increment is .1. |
| 697 |
! |
| 698 |
implicit none |
| 699 |
real, intent(in) :: value |
| 700 |
real, intent(out) :: tlcl |
| 701 |
integer, intent(in) :: i,j,bi,bj,myIter |
| 702 |
integer, intent(in) :: myThid |
| 703 |
|
| 704 |
integer :: ival |
| 705 |
real, dimension(127) :: lcltable |
| 706 |
real :: v1, v2 |
| 707 |
|
| 708 |
integer :: errorMessageUnit |
| 709 |
DATA errorMessageUnit / 15 / |
| 710 |
|
| 711 |
data lcltable/ 1.7364512e+02, 1.7427449e+02, 1.7490874e+02, & |
| 712 |
1.7554791e+02, 1.7619208e+02, 1.7684130e+02, 1.7749563e+02, & |
| 713 |
1.7815514e+02, 1.7881989e+02, 1.7948995e+02, 1.8016539e+02, & |
| 714 |
1.8084626e+02, 1.8153265e+02, 1.8222461e+02, 1.8292223e+02, & |
| 715 |
1.8362557e+02, 1.8433471e+02, 1.8504972e+02, 1.8577068e+02, & |
| 716 |
1.8649767e+02, 1.8723077e+02, 1.8797006e+02, 1.8871561e+02, & |
| 717 |
1.8946752e+02, 1.9022587e+02, 1.9099074e+02, 1.9176222e+02, & |
| 718 |
1.9254042e+02, 1.9332540e+02, 1.9411728e+02, 1.9491614e+02, & |
| 719 |
1.9572209e+02, 1.9653521e+02, 1.9735562e+02, 1.9818341e+02, & |
| 720 |
1.9901870e+02, 1.9986158e+02, 2.0071216e+02, 2.0157057e+02, & |
| 721 |
2.0243690e+02, 2.0331128e+02, 2.0419383e+02, 2.0508466e+02, & |
| 722 |
2.0598391e+02, 2.0689168e+02, 2.0780812e+02, 2.0873335e+02, & |
| 723 |
2.0966751e+02, 2.1061074e+02, 2.1156316e+02, 2.1252493e+02, & |
| 724 |
2.1349619e+02, 2.1447709e+02, 2.1546778e+02, 2.1646842e+02, & |
| 725 |
2.1747916e+02, 2.1850016e+02, 2.1953160e+02, 2.2057364e+02, & |
| 726 |
2.2162645e+02, 2.2269022e+02, 2.2376511e+02, 2.2485133e+02, & |
| 727 |
2.2594905e+02, 2.2705847e+02, 2.2817979e+02, 2.2931322e+02, & |
| 728 |
2.3045895e+02, 2.3161721e+02, 2.3278821e+02, 2.3397218e+02, & |
| 729 |
2.3516935e+02, 2.3637994e+02, 2.3760420e+02, 2.3884238e+02, & |
| 730 |
2.4009473e+02, 2.4136150e+02, 2.4264297e+02, 2.4393941e+02, & |
| 731 |
2.4525110e+02, 2.4657831e+02, 2.4792136e+02, 2.4928053e+02, & |
| 732 |
2.5065615e+02, 2.5204853e+02, 2.5345799e+02, 2.5488487e+02, & |
| 733 |
2.5632953e+02, 2.5779231e+02, 2.5927358e+02, 2.6077372e+02, & |
| 734 |
2.6229310e+02, 2.6383214e+02, 2.6539124e+02, 2.6697081e+02, & |
| 735 |
2.6857130e+02, 2.7019315e+02, 2.7183682e+02, 2.7350278e+02, & |
| 736 |
2.7519152e+02, 2.7690354e+02, 2.7863937e+02, 2.8039954e+02, & |
| 737 |
2.8218459e+02, 2.8399511e+02, 2.8583167e+02, 2.8769489e+02, & |
| 738 |
2.8958539e+02, 2.9150383e+02, 2.9345086e+02, 2.9542719e+02, & |
| 739 |
2.9743353e+02, 2.9947061e+02, 3.0153922e+02, 3.0364014e+02, & |
| 740 |
3.0577420e+02, 3.0794224e+02, 3.1014515e+02, 3.1238386e+02, & |
| 741 |
3.1465930e+02, 3.1697246e+02, 3.1932437e+02, 3.2171609e+02, & |
| 742 |
3.2414873e+02, 3.2662343e+02, 3.2914139e+02, 3.3170385e+02 / |
| 743 |
|
| 744 |
! v1 = value |
| 745 |
v1 = MIN( MAX( value, -23.d0 ), -10.4d0 ) |
| 746 |
! if (value.lt.-23.0) call error_mesg ('dargan_bettsmiller', & |
| 747 |
! 'lcltable: value too low', FATAL) |
| 748 |
|
| 749 |
! if (value.gt.-10.4) call error_mesg ('dargan_bettsmiller', & |
| 750 |
! 'lcltable: value too high', FATAL) |
| 751 |
if (value.lt.-23.0) then |
| 752 |
CALL PRINT_ERROR( 'In: dargan_bettsmiller '// & |
| 753 |
'lcltable: value too low', myThid) |
| 754 |
WRITE(errorMessageUnit,'(A,3I4,I6,1PE14.6)') & |
| 755 |
'i,j,bj,myIter,value=',i,j,bj,myIter,value |
| 756 |
! STOP 'ABNORMAL END: S/R LCLTABL (dargan_bettsmiller_mod)' |
| 757 |
endif |
| 758 |
if (value.gt.-10.4) then |
| 759 |
CALL PRINT_ERROR( 'In: dargan_bettsmiller '// & |
| 760 |
'lcltable: value too high', myThid) |
| 761 |
WRITE(errorMessageUnit,'(A,3I4,I6,1PE14.6)') & |
| 762 |
'i,j,bj,myIter,value=',i,j,bj,myIter,value |
| 763 |
! STOP 'ABNORMAL END: S/R LCLTABL (dargan_bettsmiller_mod)' |
| 764 |
endif |
| 765 |
ival = floor(10.*(v1 + 23.0)) |
| 766 |
v2 = -230. + ival |
| 767 |
v1 = 10.*v1 |
| 768 |
tlcl = (v2 + 1.0 - v1)*lcltable(ival+1) + (v1 - v2)*lcltable(ival+2) |
| 769 |
|
| 770 |
|
| 771 |
end subroutine lcltabl |
| 772 |
|
| 773 |
!####################################################################### |
| 774 |
|
| 775 |
subroutine dargan_bettsmiller_init (myThid) |
| 776 |
|
| 777 |
!----------------------------------------------------------------------- |
| 778 |
! |
| 779 |
! initialization for bettsmiller |
| 780 |
! |
| 781 |
!----------------------------------------------------------------------- |
| 782 |
|
| 783 |
integer unit,io,ierr |
| 784 |
integer, intent(in) ::myThid |
| 785 |
!------------------------------------------------------------------------------------- |
| 786 |
integer, dimension(3) :: half = (/1,2,4/) |
| 787 |
!integer :: ierr, io |
| 788 |
integer :: iUnit |
| 789 |
CHARACTER*(gcm_LEN_MBUF) :: msgBuf |
| 790 |
!------------------------------------------------------------------------------- |
| 791 |
|
| 792 |
!----------- read namelist --------------------------------------------- |
| 793 |
|
| 794 |
! if (file_exist('input.nml')) then |
| 795 |
! unit = open_file (file='input.nml', action='read') |
| 796 |
! ierr=1; do while (ierr /= 0) |
| 797 |
! read (unit, nml=dargan_bettsmiller_nml, iostat=io, end=10) |
| 798 |
! ierr = check_nml_error (io,'dargan_bettsmiller_nml') |
| 799 |
! enddo |
| 800 |
! 10 call close_file (unit) |
| 801 |
! endif |
| 802 |
|
| 803 |
!---------- output namelist -------------------------------------------- |
| 804 |
|
| 805 |
! unit = open_file (file='logfile.out', action='append') |
| 806 |
! if ( mpp_pe() == 0 ) then |
| 807 |
! write (unit,'(/,80("="),/(a))') trim(version), trim(tag) |
| 808 |
! write (unit,nml=dargan_bettsmiller_nml) |
| 809 |
! endif |
| 810 |
! call close_file (unit) |
| 811 |
|
| 812 |
! do_init=.false. |
| 813 |
|
| 814 |
! _BARRIER |
| 815 |
! _BEGIN_MASTER(myThid) |
| 816 |
CALL BARRIER(myThid) |
| 817 |
IF ( myThid.EQ.1 ) THEN |
| 818 |
|
| 819 |
WRITE(msgBuf,'(A)') 'DARGAN_BETTSMILLER_init: opening data.atm_gray' |
| 820 |
CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid ) |
| 821 |
CALL OPEN_COPY_DATA_FILE( & |
| 822 |
'data.atm_gray', 'DARGAN_BETTSMILLER_init', & |
| 823 |
iUnit, & |
| 824 |
myThid ) |
| 825 |
! Read parameters from open data file |
| 826 |
READ(UNIT=iUnit,NML=dargan_bettsmiller_nml) |
| 827 |
WRITE(msgBuf,'(A)') & |
| 828 |
'DARGAB_BETTSMILLER_INIT: finished reading data.atm_gray' |
| 829 |
CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid ) |
| 830 |
! Close the open data file |
| 831 |
CLOSE(iUnit) |
| 832 |
|
| 833 |
ENDIF |
| 834 |
CALL BARRIER(myThid) |
| 835 |
|
| 836 |
return |
| 837 |
end subroutine dargan_bettsmiller_init |
| 838 |
|
| 839 |
!####################################################################### |
| 840 |
|
| 841 |
real function mixing_ratio(vapor_pressure, pressure) |
| 842 |
|
| 843 |
! calculates the mixing ratio from the vapor pressure and pressure |
| 844 |
|
| 845 |
implicit none |
| 846 |
real, intent(in) :: vapor_pressure, pressure |
| 847 |
|
| 848 |
mixing_ratio = rdgas*vapor_pressure/rvgas/(pressure-vapor_pressure) |
| 849 |
|
| 850 |
end function mixing_ratio |
| 851 |
|
| 852 |
!####################################################################### |
| 853 |
|
| 854 |
real function virtual_temp(temp, r) |
| 855 |
|
| 856 |
! calculates the virtual temperature from the temperature and mixing ratio |
| 857 |
! consistent with the approximation used in the fms code |
| 858 |
|
| 859 |
implicit none |
| 860 |
real, intent(in) :: temp ! temperature |
| 861 |
real, intent(in) :: r ! mixing ratio |
| 862 |
|
| 863 |
real :: q ! specific humidity |
| 864 |
|
| 865 |
if (do_virtual) then |
| 866 |
q = r/(1.0+r) |
| 867 |
virtual_temp = temp*(1.0+q*(rvgas/rdgas-1.0)) |
| 868 |
else |
| 869 |
virtual_temp = temp |
| 870 |
endif |
| 871 |
|
| 872 |
end function virtual_temp |
| 873 |
|
| 874 |
!####################################################################### |
| 875 |
|
| 876 |
end module dargan_bettsmiller_mod |
| 877 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph