forked from 3rdparty/wrf-python
You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
117 lines
3.5 KiB
117 lines
3.5 KiB
C NCLFORTSTART |
|
subroutine wrfcttcalc(prs,tk,qci,qcw,qvp,ght,ter,ctt, |
|
& haveqci,nz,ns,ew) |
|
|
|
implicit none |
|
integer nz,ns,ew,haveqci |
|
double precision ght(ew,ns,nz) |
|
double precision prs(ew,ns,nz),tk(ew,ns,nz) |
|
double precision qci(ew,ns,nz),qcw(ew,ns,nz) |
|
double precision qvp(ew,ns,nz) |
|
double precision ctt(ew,ns),ter(ew,ns) |
|
c double precision znfac(nz) |
|
C NCLEND |
|
c |
|
c |
|
integer i,j,k,mjx,miy,mkzh,ripk,wrfout |
|
double precision vt,rgas,grav,opdepthu,opdepthd,dp |
|
double precision ratmix,eps,arg1,arg2,agl_hgt,ussalr |
|
double precision abscoefi,abscoef,fac,prsctt,celkel |
|
c double precision ght(ew,ns,nz),stuff(ew,ns) |
|
double precision pf(ns,ew,nz),p1,p2 |
|
c |
|
c |
|
mjx = ew |
|
miy = ns |
|
mkzh = nz |
|
eps = 0.622d0 |
|
ussalr = .0065d0 ! deg C per m |
|
rgas = 287.04d0 !J/K/kg |
|
grav = 9.81d0 |
|
abscoefi = .272d0 ! cloud ice absorption coefficient in m^2/g |
|
abscoef =.145d0 ! cloud water absorption coefficient in m^2/g |
|
celkel = 273.15d0 |
|
wrfout = 1 |
|
|
|
|
|
cCalculate the surface pressure |
|
do j=1,ew |
|
do i=1,ns |
|
ratmix = .001d0*qvp(j,i,1) |
|
arg1 = eps + ratmix |
|
arg2 = eps*(1.+ratmix) |
|
vt = tk(j,i,1) * arg1/arg2 !Virtual temperature |
|
agl_hgt = ght(j,i,nz) - ter(j,i) |
|
arg1 = -grav/(rgas*ussalr) |
|
pf(i,j,nz) = prs(j,i,1)* |
|
& (vt/(vt+ussalr*(agl_hgt)))**(arg1) |
|
enddo |
|
enddo |
|
|
|
|
|
c |
|
do j=1,ew |
|
do i=1,ns |
|
do k=1,nz-1 |
|
ripk = nz-k+1 |
|
pf(i,j,k)=.5d0*(prs(j,i,ripk)+prs(j,i,ripk-1)) |
|
enddo |
|
enddo |
|
enddo |
|
|
|
do 190 j=1,ew |
|
do 190 i=1,ns |
|
opdepthd=0.d0 |
|
k=0 |
|
|
|
c |
|
c Integrate downward from model top, calculating path at full |
|
c model vertical levels. |
|
c |
|
20 opdepthu=opdepthd |
|
k=k+1 |
|
ripk = nz-k+1 |
|
|
|
if (k.eq.1) then |
|
dp=200.d0*(pf(i,j,1)-prs(j,i,nz)) ! should be in Pa |
|
else |
|
dp=100.d0*(pf(i,j,k)-pf(i,j,k-1)) ! should be in Pa |
|
endif |
|
if (haveqci .eq. 0) then |
|
if (tk(i,j,k).lt.celkel) then |
|
c Note: abscoefi is m**2/g, qcw is g/kg, |
|
c so no convrsion needed |
|
opdepthd=opdepthu+abscoefi*qcw(j,i,k)*dp/grav |
|
else |
|
opdepthd=opdepthu+abscoef*qcw(j,i,k)*dp/grav |
|
endif |
|
else |
|
opdepthd=opdepthd+(abscoef*qcw(j,i,ripk)+ |
|
& abscoefi*qci(j,i,ripk))*dp/grav |
|
endif |
|
|
|
if (opdepthd.lt.1..and.k.lt.nz) then |
|
goto 20 |
|
elseif (opdepthd.lt.1..and.k.eq.nz) then |
|
prsctt=prs(j,i,1) |
|
else |
|
fac=(1.-opdepthu)/(opdepthd-opdepthu) |
|
prsctt=pf(i,j,k-1)+fac*(pf(i,j,k)-pf(i,j,k-1)) |
|
prsctt=min(prs(j,i,1),max(prs(j,i,nz),prsctt)) |
|
endif |
|
|
|
do 30 k=2,nz |
|
ripk = nz-k+1 |
|
p1 = prs(j,i,ripk+1) |
|
p2 = prs(j,i,ripk) |
|
if (prsctt .ge. p1 .and. prsctt .le .p2) then |
|
fac=(prsctt-p1)/(p2-p1) |
|
arg1 = fac*(tk(j,i,ripk)-tk(j,i,ripk+1))-celkel |
|
ctt(j,i) = tk(j,i,ripk+1)+ arg1 |
|
goto 40 |
|
endif |
|
30 continue |
|
40 continue |
|
190 continue |
|
return |
|
end
|
|
|