|
Setting up initial conditions
User defined subroutines tinit_i/tinit_r
c
if(initi.eq.1) call tinit_i(z,x,y,tau,lipps)
if(initi.eq.0) call tinit_r(z,tau,lipps)
c
subroutine tinit_r(z,tau,lipps)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc environmental sounding:
ccc you may enter either pressure or height for a given level
ccc (the other one should be set to zero and this subroutine
ccc will automatically calculate it); if you do not use potential
ccc temperature on input, activate the appropriate code below.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
parameter(npin=23)
dimension press(npin),temp(npin),zin(npin),
1 vap(npin),uu(npin),vv(npin)
cc pressure in hPa:
data press /
1 1008.00, 991.25, 945.50, 893.79, 836.06, 772.82, 705.22,
1 635.05, 564.48, 495.73, 430.71, 370.78, 316.72, 268.82,
1 226.98, 190.82, 159.87, 133.55, 111.29, 92.56, 52.31,
1 22.08, 9.32/
cc height in m:
DATA ZIN /NPIN*0./
cc potential temperature in K:
data temp /
1 25.26, 24.13, 21.04, 18.66, 16.50, 13.41, 9.06,
1 3.73, -1.51, -6.97, -14.09, -22.44, -30.57, -39.60,
1 -48.69, -57.40, -65.21, -72.58, -76.71, -74.98, -74.98,
1 -74.98, -74.98/
ccwater vapor mixing ratio in g/kg:
data vap /
1 0.178E+02, 0.172E+02, 0.156E+02, 0.134E+02, 0.111E+02,
1 0.888E+01, 0.631E+01, 0.487E+01, 0.396E+01, 0.200E+01,
1 0.984E+00, 0.806E+00, 0.370E+00, 0.135E+00, 0.599E-01,
1 0.258E-01, 0.123E-01, 0.582E-02, 0.367E-02, 0.589E-02,
1 0.104E-02, 0.247E-02, 0.585E-02/
cc x and y velocity components in m/s:
DATA UU/npin*0./
DATA VV/npin*0./
#include "vrtstr.fnc"
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
convert from temperature (deg C or K) into potential temperature
do k=1,npin
c temp(k)=temp(k)*(1.e3/press(k))**(rg/cp)
temp(k)=(temp(k)+273.16)*(1.e3/press(k))**(rg/cp)
enddo
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
utr=0.
do ip=1,npin
uu(ip)=uu(ip)-utr
enddo
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
compute approximated height of pressure levels:
if (zin(2).lt.1.e-4) then
zin(1)=0.
do k=2,npin
km=k-1
tempk =temp(k )*(1.e3/press(k ))**(-rg/cp)
. * (1.+.6e-3*vap(k ))
tempkm=temp(km)*(1.e3/press(km))**(-rg/cp)
. * (1.+.6e-3*vap(km))
delt=tempk-tempkm
if (delt.gt.1.e-4) then
tavi=alog(tempk/tempkm)/delt
else
tavi=1./tempk
endif
deltz=-rg/(tavi*g) * alog(press(k)/press(km))
zin(k)=zin(km)+deltz
end do
end if
compute approximate pressure of height levels:
if (press(1).lt.1.e-4) then
press(1)=pr00/100.
do k=2,npin
km=k-1
tempk =temp(k )
tempkm=temp(km)
delt=tempk-tempkm
if (delt.gt.1.e-4) then
tavi=alog(tempk/tempkm)/delt
else
tavi=1./tempk
endif
rc=rg/cp
rci=1./rc
press(k)=( press(km)**rc -
. g*(1.e3)**rc*tavi*(zin(k)-zin(km))/cp )**rci
enddo
endif
compute environmental profiles from sounding assuming no topography:
cc surface data:
iisn=1
the1(1)=temp(iisn)
tme1(1)=the1(1) * (1000./press(iisn))**(-rg/cp)
qve1(1)=vap(iisn)*1.e-3
ue1(1)=uu(1)
ve1(1)=vv(1)
c print*,'DETERMINED SURFACE DATA'
cc higher levels - interpolate:
c print*,'INTERPOLATION TO HIGHER LEVELS'
do k=2,l
c print*,'k=',k
do kk=2,npin
iisn=kk-1
if(zin(kk).ge.z(k)) go to 665
enddo
c print*,'INPUT SOUNDING DOES NOT GO HIGH ENOUGH. STOP.'
stop 'SOUNDING'
665 continue
c print*,'iisn=',iisn
coe2=(z(k)-zin(iisn))/(zin(iisn+1)-zin(iisn))
the1(k)=coe2*temp(iisn+1) + (1.-coe2)*temp(iisn)
qve1(k)=(coe2*vap(iisn+1) + (1.-coe2)*vap(iisn))*1.e-3
presnl=coe2*press(iisn+1) + (1.-coe2)*press(iisn)
tme1(k)=the1(k) * (1000./presnl)**(-rg/cp)
ue1(k)=coe2*uu(iisn+1) + (1.-coe2)*uu(iisn)
ve1(k)=coe2*vv(iisn+1) + (1.-coe2)*vv(iisn)
end do
c smooth environmental profiles
c filter 2 and 4 dz waves
call filtprf(the1,l,1,1)
call filtprf(ue1 ,l,1,1)
call filtprf(ve1 ,l,1,1)
call filtprf(tme1,l,1,1)
call filtprf(qve1,l,1,1)
create environmental profiles; topography is now considered:
do i=1,np
do j=1,mp
do k=1,l
zcr(k)=zstr(k)/gi(i,j)+zs(i,j)
end do
if (moist.eq.1) then
do k=1,l
kk=zcr(k)/dz+1
kk=min(l,max(1,kk))
kkp=min(l,kk+1)
coe=(zcr(k)-z(kk))/dz
tme(i,j,k)=coe*tme1(kkp)+(1.-coe)*tme1(kk)
qve(i,j,k)=coe*qve1(kkp)+(1.-coe)*qve1(kk)
the(i,j,k)=coe*the1(kkp)+(1.-coe)*the1(kk)
ue(i,j,k)=coe*ue1(kkp)+(1.-coe)*ue1(kk)
ve(i,j,k)=coe*ve1(kkp)+(1.-coe)*ve1(kk)
compute reference state vertical profiles for every x,y point
call thprof(tau(1,i,j,1),zcr,l,lipps)
do k=1,l
th0(i,j,k)=tau(k,i,j,1)
end do
call rhprof(tau(1,i,j,1),zcr,l,lipps)
do k=1,l
dnmi=1./(stryy(i,j)*strxx(i,j)-stryx(i,j)*strxy(i,j))
rho(i,j,k)=tau(k,i,j,1)/(gi(i,j)*gmus(k))
+ *((1-icylind)*gmm(i,j,k)**2*cosa(i,j)
+ +icylind*gmm(i,j,k))*dnmi
call absorber(z,x,y,tau)
call theimpl(tau) ! theta_e derivatives for implict model
create flux and roughness fields Z.S.
c call nois2(hfx,zo,1)
do 30 j=1,mp
do 30 i=1,np
hfx(i,j)=hf00
C *(1.+0.01*hfx(i,j)) ! heat flux in K*m/s
30 zo(i,j)=0.0 ! roughness parameter, m
if (moist.eq.1) then
c call nois2(qfx,zo,1)
do 31 j=1,mp
do 31 i=1,np
qfx(i,j)=qf00
C *(1.+0.01*qfx(i,j)) ! spec.hum.flux in kg/kg*m/s
31 zo(i,j)=0.0 ! roughness parameter, m
subroutine tinit_i(z,x,y,tau,lipps)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
compute basic and environmental state profiles for every x,y point
c --- CAUTION: in general these will be time variable!
do 10 j=1,mp
do 10 i=1,np
c --- physical altitude
do 11 k=1,l
11 zcr(k)=zstr(k)/gi(i,j)+zs(i,j)
c------------------------------------------------------------------c
C IDEALIZED BASIC STATE PROFILES c
c c
c NOTE: pressure and temperature profiles c
c are NOT computed unless MOISTMOD > 0 c
c c
c lipps = 0: This is classical BOUSSINESQ approximation. c
c Constant density, linearly decreasing c
c (p,T) -> there exists max. altitude c
c beyond which profiles are nonphysical c
c Th = T. c
c c
c lipps = 1: Ogura and Phillips ANELASTIC: see JAS V. c
c 19, pp. 173-179, 1962. Characterized c
c by single scale height - hr. c
c Th = constant (isentropic atm.) c
c T = linearly decreasing profile. c
c -> there exists a height restriction c
c beyond which profiles are nonphysical c
c c
c lipps = 2: Clark and Farley ANELASTIC: see JAS V.41, c
c pp. 329-350, 1984. Assumes exp. profile c
c for th with scale height - ht. Value of c
c cp (or equivalently, cap) is chosen c
c as second indepenent parameter. c
c (rho,p,T) are all nonlinearly decreasing c
c and go to zero value at finite altitude c
c -> there exists a height restriction c
c beyond which profiles are nonphysical c
c c
c lipps = 3: Bacmeister and Schoeberl ANELASTIC: see c
c JAS V. 46, pp. 2109-2134, 1989. Assumes c
c exp. profiles for rho and th with scale c
c heights hr and ht, respectively. c
c -> T = constant (isothermal atm.) c
c No height restriction. c
c c
c------------------------------------------------------------------c
c --- basic state potential temperature
call thprof(tau(1,i,j,1),zcr,l,lipps)
do 12 k=1,l
12 th0(i,j,k)=tau(k,i,j,1)
c --- basic state density
call rhprof(rhobr,zcr,l,lipps)
do 13 k=1,l
dnmi=1./(stryy(i,j)*strxx(i,j)-stryx(i,j)*strxy(i,j))
rho(i,j,k)=rhobr(k)/(gi(i,j)*gmus(k))
+ *((1-icylind)*gmm(i,j,k)**2*cosa(i,j)
+ +icylind*gmm(i,j,k))*dnmi
rhoi(i,j,k)=1./rho(i,j,k)
13 continue
c------- construct environmental wind
if(u0z.ge.1.e-15) then
do k=1,l
ue0(k)=u00+u0z*zcr(k) ! specified shear
enddo
else
do k=1,l
ue0(k)=u00*gmm(i,j,k) ! ~constant in altitude
enddo
endif
u0s=u00*rdsi
c------- construct geostrophic environmental state
do 14 k=1,l
aa(k)=u0s*(2.*ue0(k)+gmm(i,j,k)*fcr0*rds)+ue0(k)*fcr0
bb(k)=ue0(k)*(fcr0+ue0(k)/gmm(i,j,k)*rdsi)
ue(i,j,k)=ue0(k)*cosa(i,j)
ve(i,j,k)=v00
the0(k)=th0(i,j,k)
if(lipps.eq.0) the0(k)=th00*(1.+ st*zcr(k))
the(i,j,k)=the0(k)+(th0(i,j,k)/g)*(
. isphere*sina(i,j)**2*(bb(k)-.5*aa(k))
. +(1-isphere)*fcr3(i,j)*(v0z*xcr(i,j)-u0z*ycr(i,j)) )
14 continue
ihs=0 ! Held-Suarez flow
if(ihs.eq.1) then
c------------------------------------------------------------------c
c c
c See Held and Suarez, Bulletin of the AMS, V. 75, c
c pp. 1825-1830, 1994 c
c c
c NOTE: th00 = 315K for HS paper c
c c
c------------------------------------------------------------------c
c
c --- construct dry part held-suarez environmental state
hskf=1./( 24.*3600.)
hska=1./(40.*24.*3600.)
hsks=1./(4. *24.*3600.)
sigb=0.7
sigi=1./(1.-sigb)
dthy=60.
dthz=10.
hd=7000. ! density scale height
xs=th00-dthy*sina(i,j)**2
xc=dthz*cosa(i,j)**2
do 5 k=1,l
sgm=exp(-zcr(k)/hd) ! Basic state (p/p_ref) approx.
sg0=sgm**cap ! conversion factor for temp->theta
thstr=200./sg0
dfhs(i,j,k)=amax1(0., sigi*(sgm-sigb) )
theq(i,j,k)=amax1(thstr,xs-xc*alog(sgm)) ! full HS profile
c the(i,j,k)=theq(i,1,k) ! sub HS profile
theqtmp = theq(i,1,k) ! sub HS profile
if(ice.eq.1) then ! COLD RAIN
a=rg/rv
b=hlats/rv
d=hlatv/rv
e=-cp/rg
do k=1,l
coe_l=comb(tme(i,j,k),tdn,tup,0.,1.) ! liquid contribution
thetme=the(i,j,k)/tme(i,j,k)
pre=1.e5*thetme**e ! modify?
tti=tme(i,j,k)
delt=(tti-t00)/(tti*t00)
esw=ee0*exp(d * delt)
esi=ee0*exp(b * delt)
qvsw=a * esw /(pre-esw) ! modify?
qvsi=a * esi /(pre-esi) ! modify?
qvs=coe_l*qvsw + (1.-coe_l)*qvsi
qve(i,j,k)=qve(i,j,k)*qvs
c --- held-suarez test
if(ihs.eq.1) then
tti=theq(i,j,k)/thetme
coe_l=comb(tti,tdn,tup,0.,1.) ! liquid contribution
delt=(tti-t00)/(tti*t00)
esw=ee0*exp(d * delt)
esi=ee0*exp(b * delt)
qvsw=a * esw /(pre-esw) ! modify?
qvsi=a * esi /(pre-esi) ! modify?
qvs=coe_l*qvsw + (1.-coe_l)*qvsi
qveq(i,j,k)=0.9*qvs ! set relative humidity
end if
c --- end HS test
end do
else ! WARM RAIN
a=rg/rv
b=hlat/(rv*t00)
c=hlat/cp
d=hlat/rv
e=-cp/rg
do k=1,l
thetme=the(i,j,k)/tme(i,j,k)
thi=1./the(i,j,k)
c thetme=the0(k)/te0(k) !uniform qv on the planet
c thi=1./the0(k) !uniform qv on the planet
pre=1.e5*thetme**e ! modify?
ylc=b*thetme*t00*thi
ees=ee0*exp(b-ylc)
qvs=a*ees/(pre-ees) ! modify?
qve(i,j,k)=qve(i,j,k)*qvs
c --- held-suarez test
if(ihs.eq.1) then
thi=1./theq(i,j,k)
ylc=b*thetme*t00*thi
ees=ee0*exp(b-ylc)
qvs=a*ees/(pre-ees) ! modify?
qveq(i,j,k)=0.9*qvs ! set relative humidity
end if
c --- end HS test
end do
endif
#endif
#endif
call absorber(z,x,y,tau)
call theimpl(tau) ! theta_e derivatives for implict model
create flux and roughness fields Z.S.
do 30 j=1,mp
do 30 i=1,np
hfx(i,j)=hf00 ! heat flux in K*m/s
30 zo(i,j)=0.1 ! roughness parameter, m
if (moist.eq.1) then
do 31 j=1,mp
do 31 i=1,np
31 qfx(i,j)=qf00 ! spec.hum.flux in kg/kg*m/s
endif
return
end
subroutine thprof(th0,z,l,lipps)
if(lipps.eq.0) then ! Boussinesq profile
do k=1,l
c th0(k)=th00-g*z(k)/rg ! CAUTION: max. alt. exists
th0(k)=th00
end do
endif
if(lipps.eq.1) then ! Ogura-Phillips profile
do k=1,l
th0(k)=th00
end do
endif
if(lipps.eq.2) then ! Clark-Farley profile
do k=1,l
th0(k)=th00*exp(st*z(k))
end do
endif
if(lipps.eq.3) then ! Bacmeister-Schoeberl profile
do k=1,l
th0(k)=th00*exp(st*z(k))
end do
endif
subroutine rhprof(rh0,z,l,lipps)
if(lipps.eq.0) then ! Boussinesq profile
do k=1,l
rh0(k)=rh00
end do
endif
if(lipps.eq.1) then ! Ogura-Phillips profile
capi=1./cap
sdi=1./7000. ! density scale height (m^-1)
do k=1,l
rh0(k)=rh00*(1.-cap*z(k)*sdi)**(capi-1.)
end do ! CAUTION: max. alt. exists
endif
if(lipps.eq.2) then ! Clark-Farley profile
capi=1./cap
cs=g/(cp*tt00*st)
do k=1,l
exs=exp(-st*z(k)) ! CAUTION: max. alt. exists
rh0(k)=rh00*exs*(1.-cs*(1.-exs))**(capi-1.)
end do
endif
if(lipps.eq.3) then ! Bacmeister-Schoeberl profile
sdi=1./6622.51 ! density scale height (m^-1)
do k=1,l
rh0(k)=rh00*exp(-sdi*z(k))
end do
endif
subroutine teprof(pm0,tm0,z,l,lipps)
if(lipps.eq.0) then ! Boussinesq profile
do k=1,l
pm0(k)=pr00-rh00*g*z(k) ! CAUTION: max. alt. exists
tm0(k)=tt00-g*z(k)/rg ! CAUTION: max. alt. exists
end do
end if
if(lipps.eq.1) then ! Ogura-Phillips profile
capi=1./cap
sdi=1./7000. ! density scale height (m^-1)
do k=1,l
pm0(k)=pr00*(1.-cap*z(k)*sdi)**capi ! CAUTION: max. alt. exists
tm0(k)=tt00*(1.-cap*z(k)*sdi) ! CAUTION: max. alt. exists
end do
end if
if(lipps.eq.2) then ! Clark-Farley profile
capi=1./cap
cs=g/(cp*tt00*st)
do k=1,l
exs=exp(-st*z(k))
pm0(k)=pr00*(1.-cs*(1.-exs))**capi ! CAUTION: max. alt. exists
tm0(k)=tt00/exs*(1.-cs*(1.-exs)) ! CAUTION: max. alt. exists
end do
end if
if(lipps.eq.3) then ! Bacmeister-Schoeberl profile
sdi=1./7000. ! density scale height (m^-1)
do k=1,l
pm0(k)=pr00*exp(-sdi*z(k))
tm0(k)=pr00/(rg*rh00)
end do
end if
Initializing flow fields:
do 4 k=1,l
do 4 j=1,mp
do 4 i=1,np
ia = (npos-1)*np + i
ja = (mpos-1)*mp + j
th(i,j,k)=0.
u(i,j,k,0)=ue(i,j,k)
v(i,j,k,0)=ve(i,j,k)
w(i,j,k,0)=0.e-3*fz(i,j,k)
p(i,j,k) =0.
if(ivis.eq.1) then
do k=1,l
do j=1,mp
do i=1,np
tke(i,j,k)=0.
if(ichm.eq.1) then
do ispc=1,nspc
do k=1,l
do j=1,mp
do i=1,np
if(ispc.eq.1) chm(i,j,k,ispc) = zsib(i,j,k)
c chm(i,j,k,ispc) = zsib(i,j,k)
chm(i,j,k,ispc) = 0.
fchm(i,j,k,ispc) = 1000.0
do 40 k=1,l
do 40 j=1,mp
do 40 i=1,np
qv(i,j,k)=qve(i,j,k)
qc(i,j,k)=0.
qr(i,j,k)=0.
compute poisson equation for potential
if(ispcpr.eq.1) call spcpri
call gcrk( p,pfx,pfy,pfz,ox(1-ih,1-ih,1,0),
. oy(1-ih,1-ih,1,0),oz(1-ih,1-ih,1,0),
. fx,fy,fz,ft,itp0,epp0,0)
call prforc(p,pfx,pfy,pfz,ox(1-ih,1-ih,1,0),
. oy(1-ih,1-ih,1,0),oz(1-ih,1-ih,1,0),
. fx,fy,fz,ft)
|
Eulag
