! read 20090416/high/height;temper;t-td;uv
! lon 32~160 interval: 4 degree
! lat 12~80 interval: 4 degree
! time 1~12: 2009.04.16.20:00~04.22.08:00 interval: 12 hours
! U,v,height,t: 1000,925,850,700,500,400,300,250,200,150,100 hpa
! BeiJing Time
!: 2010.11.27,
REAL,PARAMETER :: Omega=7.292e-5,R=6371e+3,PI=3.1415926,Delta=4,DeltaT=12,epsilonPhi=1.e3,epsilonPsi=1.e3
INTEGER,PARAMETER :: nx=33,ny=18,nz=11,nt=12
REAL F,EE,AA,BB,sign
INTEGER WtTopLevel(NX,NY,NT),Pn,Qn,Ltotal,Lnum
REAl P(NZ),deltP(10),lat(ny),lon(nx),sigmadeltD,dx,dy,epsilon,aaa,bbb
CHARACTER timefile(12)*12 ,levelfile(11)*4
real,allocatable :: temper(:,:,:,:),q(:,:,:,:),u(:,:,:,:), &
v(:,:,:,:),height(:,:,:,:),w(:,:,:,:), &
qu(:,:,:,:),qv(:,:,:,:),adq(:,:,:,:), &
adqv(:,:,:,:),theta(:,:,:,:),Wt(:,:,:,:), &
vorg(:,:,:,:),voro(:,:,:,:),div(:,:,:,:), &
deltD(:,:,:),ttd(:,:,:,:),thetaP(:,:,:,:), &
phi(:,:,:,:),residual(:,:,:,:), &
psi(:,:,:,:),advoro(:,:,:,:),advorg(:,:,:,:)
allocate (temper(NX,NY,NZ,NT),q(NX,NY,nz,NT),u(NX,NY,NZ,NT), &
v(NX,NY,NZ,NT),height(NX,NY,NZ,NT),w(NX,NY,NZ,NT), &
qu(NX,NY,NZ,NT),qv(NX,NY,NZ,NT),adq(NX,NY,NZ,NT), &
adqv(NX,NY,NZ,NT),theta(NX,NY,NZ,NT),Wt(NX,NY,NZ,NT), &
vorg(NX,NY,NZ,NT),voro(NX,NY,NZ,NT),div(NX,NY,NZ,NT), &
deltD(NX,NY,NT),ttd(NX,NY,NZ,NT),thetaP(NX,NY,NZ,NT), &
phi(NX,NY,NZ,NT),residual(NX,NY,NZ,NT), &
psi(NX,NY,NZ,NT),advoro(NX,NY,NZ,NT),advorg(NX,NY,NZ,NT))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
DATA P /1000,925,850,700,500,400,300,250,200,150,100/
DATA deltP/75, 75, 50, 200,100,100,50, 50, 50, 50/
DATA
timefile/'09041620.000','09041708.000','09041720.000','09041808.000','09041820.000','090419 08.000', &
'09041920.000','09042008.000','09042020.000','09042108.000','09042120.000','09042208.000'/ DATA levelfile/'1000','925','850','700','500','400','300','250','200','150','100'/
do i=1,nx
lon(i)=32+(i-1)*Delta
lon(i)=lon(i)*pi/180.
enddo
do i=1,ny
lat(i)=12+(i-1)*Delta
lat(i)=lat(i)*pi/180.
enddo
!--------------read temperature----------------
DO IZ=1,NZ ! level
DO IT=1,NT ! time
OPEN(11,file='c:\micaps\temper\'//trim(levelfile(iz))//'\'//timefile(it))
do i=1,4
read(11,*)
enddo
do ij=NY,1,-1
read(11,*)(temper(ii,ij,iz,it),ii=1,NX) ! C
enddo
CLOSE(11)
ENDDO
ENDDO
!--------------read u,v -------------------------
DO IZ=1,NZ ! level
DO IT=1,NT ! time
OPEN(11,file='c:\micaps\uv\'//trim(levelfile(iz))//'\'//timefile(it))
do i=1,3
read(11,*)
enddo
do ij=NY,1,-1
read(11,*)(u(ii,ij,iz,it),ii=1,NX) ! m/s
enddo
do ij=NY,1,-1
read(11,*)(v(ii,ij,iz,it),ii=1,NX) ! m/s
enddo
CLOSE(11)
ENDDO
ENDDO
!--------------read t-td -------------------------
DO IZ=1,NZ ! level
DO IT=1,NT ! time
OPEN(11,file='c:\micaps\t-td/'//trim(levelfile(iz))//'\'//timefile(it))
do i=1,4
read(11,*)
enddo
do ij=NY,1,-1
read(11,*)(ttd(ii,ij,iz,it),ii=1,NX) ! C
enddo
CLOSE(11)
ENDDO
ENDDO
!--------------read height -------------------------
DO IZ=1,NZ ! level
DO IT=1,NT ! time
OPEN(11,file='c:\micaps\height\'//trim(levelfile(iz))//'\'//timefile(it))
do i=1,4
read(11,*)
enddo
do ij=NY,1,-1
read(11,*)(height(ii,ij,iz,it),ii=1,NX) ! 10 m
enddo
CLOSE(11)
ENDDO
ENDDO
height=height*10.
!@@@@@@@@@@@@@@@@@@@@ read end @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @@@@@@@@@@
!======== t,td => specific humidity:
DO IT=1,NT
DO IZ=1,NZ
do ix=1,NX
do iy=1,NY
if(temper(ix,iy,iz,it).ge.-15)then
AA=17.2693882
BB=35.86
elseif(temper(ix,iy,iz,it).le.-40)then
AA=21.8745584
BB=7.66
else !!!! ?????????????????????????????????????
AA=21.8745584+(17.2693882-21.8745584)/(-15-(-40.))*(temper(ix,iy,iz,it)-(-40))
BB=7.66+(35.86-7.66)/(-15-(-40.))*(temper(ix,iy,iz,it)-(-40))
endif
EE=6.1078*exp( AA*(temper(ix,iy,iz,it)-ttd(ix,iy,iz,it))/(273.16+(temper(ix,iy,iz,it) & -ttd(ix,iy,iz,it))-BB) )
q(ix,iy,iz,it)=622.*EE/(P(iz)-0.378*EE)
enddo
enddo
ENDDO
ENDDO
!=====t => T:
temper=temper+273.16
!######################## geostrophic wind vorticity ########################### print*,'geostrophic wind vorticity'
vorg=1.e+36
DO it=1,NT
DO iz=1,NZ
do i=1+1,NX-1
do j=1+1,NY-1
F=Omega*2*sin(lat(j))
vorg(i,j,iz,it)=9.8/(F*R**2)*( (height(i+1,j,iz,it)+height(i-1,j,iz,it)- &
2*height(i,j,iz,it))/(Delta*pi/180*cos(lat(j)))**2 + &
(height(i,j+1,iz,it)+height(i,j-1,iz,it)- &
2*height(i,j,iz,it))/(Delta*pi/180)**2 - &
(height(i,j+1,iz,it)-height(i,j-1,iz,it))/(Delta*pi/180*2)*tan(lat(j)) )
enddo
enddo
ENDDO
ENDDO
!######################## vorticity of observed wind ######################## print*,'vorticity of observed wind'
voro=1.e+36
DO it=1,NT
DO iz=1,NZ
do i=1+1,NX-1
do j=1+1,NY-1
voro(i,j,iz,it)=1/(R*2)*( (v(i+1,j,iz,it)-v(i-1,j,iz,it))/(Delta*pi/180*cos(lat(j))) - &
(u(i,j+1,iz,it)-u(i,j-1,iz,it))/(Delta*pi/180) + &
2*u(i,j,iz,it)*tan(lat(j)) )
enddo
enddo
ENDDO
ENDDO
!######################## advection of vorticity ########################
print*,'vorticity of observed wind'
advoro=1.e+36
advorg=1.e+36
DO it=1,NT
DO iz=1,NZ
do i=2+1,NX-2
do j=2+1,NY-2
advoro(i,j,iz,it)=-u(i,j,iz,it)*(voro(i+1,j,iz,it)-voro(i-1,j,iz,it))/(2*R*Delta*pi/180.*cos(lat(j))) - &
v(i,j,iz,it)*(voro(i,j+1,iz,it)-voro(i,j-1,iz,it))/(2*R*Delta*pi/180.)
advorg(i,j,iz,it)=-u(i,j,iz,it)*(vorg(i+1,j,iz,it)-vorg(i-1,j,iz,it))/(2*R*Delta*pi/180.*cos(lat(j))) - &
v(i,j,iz,it)*(vorg(i,j+1,iz,it)-vorg(i,j-1,iz,it))/(2*R*Delta*pi/180.) enddo
enddo
ENDDO
ENDDO
!######################## divergence ########################### print*,'divergence'
div=1.e+36
DO it=1,NT
DO iz=1,NZ
do i=1+1,NX-1
do j=1+1,NY-1
div(i,j,iz,it)=1./(2*R)*( (u(i+1,j,iz,it)-u(i-1,j,iz,it))/(Delta*pi/180*cos(lat(j))) + &
(v(i,j+1,iz,it)-v(i,j-1,iz,it))/(Delta*pi/180) - &
2*v(i,j,iz,it)*tan(lat(j)) )
enddo
enddo
ENDDO
ENDDO
!######################## vertical velocity ########################### print*,'vertical velocity '
w=1.e+36
DO it=1,NT
DO iz=2,NZ
do i=1+2,NX-2
do j=1+2,NY-2
w(i,j,1,it)=0.
w(i,j,iz,it)=w(i,j,iz-1,it)+0.5*(div(i,j,iz,it)+div(i,j,iz-1,it))*(p(iz-1)-p(iz))
enddo
enddo
ENDDO
ENDDO
!-------- correct w --------------
!- - - - - thermaldynamic method for w at the top level - - - - - - -
!- - - - - 100,150hpa: reference to 100,150,200 hpa
theta=1.e+36
wt=1.e+36
DO IT=1,NT
DO IZ=1,NZ
do i=1,NX
do j=1,NY
theta(i,j,iz,it)=temper(i,j,iz,it)*(1000/p(iz))**0.286
enddo
enddo
ENDDO
ENDDO
iz=NZ ! 100hpa
DO it=1,nt-1
do i=1+1,NX-1
sign=0.
thetaP(i,j,iz,it)=1.e+36
IF(theta(i,j,iz-1,it)-theta(i,j,iz,it).ne.0)THEN
sign=1
thetaP(i,j,iz,it)=(theta(i,j,iz-1,it)-theta(i,j,iz,it))/deltP(iz-1) ! 100hpa: forward interpolation
ELSE
iz=NZ-1 ! 150 hpa
if(theta(i,j,iz-1,it)-theta(i,j,iz+1,it).ne.0)then
sign=1
thetaP(i,j,iz,it)=(theta(i,j,iz-1,it)-theta(i,j,iz+1,it))/(deltP(iz-1)+deltP(iz)) ! 150hpa: central difference interpolation
elseif(theta(i,j,iz,it)-theta(i,j,iz+1,it).ne.0)then
sign=1
thetaP(i,j,iz,it)=(theta(i,j,iz,it)-theta(i,j,iz+1,it))/deltP(iz) ! 150hpa: backward interpolation
endif
ENDIF
! if(thetaP(i,j,iz,it).ne.1.e+36)thetaP(i,j,iz,it)=thetaP(i,j,iz,it) ! hpa=>pa: it's unnecessary
Wt(i,j,iz,it)=( (theta(i,j,iz,it+1)-theta(i,j,iz,it))/(DeltaT*60*60)+ &
u(i,j,iz,it)*(theta(i+1,j,iz,it)-theta(i-1,j,iz,it))/(2*R*cos(lat(j))*Delta*pi/180.)+ &
v(i,j,iz,it)*(theta(i,j+1,iz,it)-theta(i,j-1,iz,it))/(2*R*Delta*pi/180.) ) &
/(-1.*thetaP(i,j,iz,it))
if(sign.eq.0)Wt(i,j,iz,it)=1.e+36
WtTopLevel(i,j,it)=iz ! the Top level of Wt
enddo
enddo
ENDDO
!- - - - -
do i=1+2,NX-2
do j=1+2,NY-2
!!!!!!!!!!!!!!: first revised scheme
!sigmadeltP=0.
!do iz=1,WtTopLevel(i,j,it)-1
!sigmadeltP=sigmadeltP+deltP(iz)
!enddo
! IF(Wt(i,j,WtTopLevel(i,j,it),it).ne.1.e+36)THEN
! deltD(i,j,it)=(w(i,j,WtTopLevel(i,j,it),it)-Wt(i,j,WtTopLevel(i,j,it),it))/sigmadeltP ! do iz=1+1,WtTopLevel(i,j,it)
! w(i,j,iz,it)=w(i,j,iz,it)-(P(1)-P(iz))*deltD(i,j,it)
! div(i,j,iz,it)=div(i,j,iz,it)-deltD(i,j,it)
! enddo
! ELSE
! w(i,j,:,it)=1.e+36
! ENDIF
!!!!!!!!!!!!!!: second revised scheme
sigmadeltP=0.
do iz=1,INT((1000.-P(WtTopLevel(i,j,it)))/10)
sigmadeltP=sigmadeltP+iz
enddo
IF(Wt(i,j,WtTopLevel(i,j,it),it).ne.1.e+36)THEN
deltD(i,j,it)=(w(i,j,WtTopLevel(i,j,it),it)-Wt(i,j,WtTopLevel(i,j,it),it))/sigmadeltP
do iz=1+1,WtTopLevel(i,j,it)
w(i,j,iz,it)=w(i,j,iz,it)-((1000-P(iz))/10)*((1000-P(iz))/10+1)/2.0*deltD(i,j,it) div(i,j,iz,it)=div(i,j,iz,it)-((1000-P(iz))/10)/10*deltD(i,j,it)
enddo
ELSE
w(i,j,:,it)=1.e+36
ENDIF
enddo
enddo
ENDDO
!############### vapor flux & vapor flux divergence ##################
print*,'vapor flux & vapor flux divergence'
w(:,:,:,nt)=1.e+36
w(1,:,:,:) =1.e+36
w(:,1,:,:) =1.e+36
w(nx,:,:,:)=1.e+36
w(:,ny,:,:)=1.e+36
adq =1.e+36
adqv=1.e+36
qu=1.e+36
qv=1.e+36
DO IT=1,nt
DO IZ=1,nz
DO i=1+1,NX-1
DO j=1+1,NY-1
qv(i,j,iz,it)=q(i,j,iz,it)*v(i,j,iz,it)/9.8
qu(i,j,iz,it)=q(i,j,iz,it)*u(i,j,iz,it)/9.8
adq(i,j,iz,it)=u(i,j,iz,it)*((q(i+1,j,iz,it)-q(i-1,j,iz,it))/(2*R*cos(lat(j))*Delta*pi/180)) + &
v(i,j,iz,it)*((q(i,j+1,iz,it)-q(i,j-1,iz,it))/(2*R*Delta*pi/180)) adqv(i,j,iz,it)=(q(i,j,iz,it)*div(i,j,iz,it)+adq(i,j,iz,it))*1e5/9.8
ENDDO
ENDDO
ENDDO
ENDDO
!################## stream function #################
!################## velocity potential #################
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
OPEN(8,FILE='c:\dignostic\shixi\out-micaps.grd',form='binary')
do it=1,nt
WRITE(8) (((u(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((v(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((q(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((temper(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((height(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((voro(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((vorg(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((div(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((advoro(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((advorg(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((adqv(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((qu(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((qv(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((adq(i,j,iz,it)*1e5/9.8,i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((q(i,j,iz,it)*div(i,j,iz,it)*1e5/9.8,i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((w(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
! WRITE(8) (((psi(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
! WRITE(8) (((phi(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
enddo
CLOSE(8)
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
deallocate (temper,q,u,v,height,w,qu,qv,adq,adqv,theta,Wt, &
vorg,voro,div,deltD,ttd,thetaP,phi,residual,psi,advoro,advorg)
END