File: C:\NOAA\NEMS_11731\src\atmos\gfs\phys\gbphys_adv_hyb_gc.f

1           subroutine gbphys_adv(im,ix,deltim,told,uold,vold,rqold,ps,
2          &                            tg ,ug ,vg ,rqg ,prsi )
3     c
4     ! version of ppi=ak+bk*psfc+ck*(tv/t0)^(1/kappa)
5     !
6     ! hmhj : this is modified hybrid by finite difference from henry juang
7     !
8     ! input
9     !	im	number of grid point in given latitude to compute
10     !	ix	dimension of grid point number in a given latitude
11     !	deltim	time step in second
12     !	told	(ix,levs) virtual temperature before physics
13     !	uold	(ix,levs) pseudo u-wind before physics
14     !	vold	(ix,levs) pseudo v-wind before physics
15     !	rqold	(ix,levs) tracers before physics
16     !	prsi	(ix,levs+1) presure at interface
17     ! input and output
18     !	tg	(ix,levs) virtual temperature after physics then update here
19     !	ug	(ix,levs) pseudo u-wind after physics
20     !	vg	(ix,levs) pseudo v-wind after physics
21     !	rqg	(ix,levs) tracers after physics
22     !
23           use machine , only : kind_phys
24                                                                                     
25           use resol_def
26           use coordinate_def
27           use physcons , rd => con_rd , cp => con_cp
28     
29           implicit none
30           integer ix,im
31           integer i,k,kk,n,nvcn,ifirst,lat
32           real kappa,rkappa,deltim,dt2,rdt2
33           real det,tkrt0,wmkm1,wmkp1,xvcn
34           real ps(ix),
35          1    tg(ix,levs),told(ix,levs),
36          2    ug(ix,levs),uold(ix,levs),
37          3    vg(ix,levs),vold(ix,levs),
38          4    rqg(ix,levs,ntrac),rqold(ix,levs,ntrac),
39          5    prsi(ix,levs+1)
40           real  dtvdt(ix,levs)
41           real  rtm,rtp,tfm,tfp,afk,afp,tgm,tgp,onethird,twothird
42     c
43           real rdelp(im,levs)
44           real tki,tkci(im,levs+1)
45           real dpp(im,levs),rpp(im,levs)
46           real wflx(im,levs+1),wf(im,levs+1),wf1(im,levs+1)
47           real wml(im,levs),wmm(im,levs),wmu(im,levs)
48           real work(im,levs)
49           real dup(im,levs),dum(im,levs)
50           real alpha(im,levs),betta(im,levs),gamma(im,levs)
51           real delta(im,levs)
52           real ppi(im,levs+1),ppl(im,levs)
53           real zadv(im,levs,3+ntrac)
54           logical adjusted
55     c
56     !     print *,' enter  gbphys_adv_hyb_gc_fd '		! hmhj check
57     c
58           dt2=deltim*2.0
59           rdt2=1./dt2
60     !
61           kappa   = rd / cp
62           rkappa  = cp / rd
63     
64     !
65     ! ----- prepare ck/kappa/t(tv/t0)^(1/kappa)
66     !
67           tkci= 0.0
68           do k=2,levs
69             do i=1,im
70               tkrt0 = (told(i,k-1)+told(i,k))/(thref(k-1)+thref(k))
71               tki   = ck5(k)*tkrt0**rkappa
72               tkci(i,k)=tki*rkappa/(told(i,k-1)+told(i,k))
73             enddo
74           enddo
75           do k=1,levs+1
76             do i=1,im
77               ppi(i,k)=prsi(i,k)	! from hyb2press with overturn adjusted
78             enddo
79           enddo
80           do k=1,levs
81             do i=1,im
82               rpp(i,k)=1./(prsi(i,k)+prsi(i,k+1))
83               dpp(i,k)=prsi(i,k)-prsi(i,k+1)
84               rdelp(i,k) = 0.5/dpp(i,k)
85               ppl(i,k)=0.5*(ppi(i,k)+ppi(i,k+1))
86     ! temperature tendency = temperature change / 2dt
87               dtvdt(i,k)=(tg(i,k)-told(i,k))*rdt2
88             enddo
89           enddo
90     !
91           alpha(:,levs)=0.0
92           betta(:,   1)=0.0
93           do k=2,levs
94             do i=1,im
95               alpha(i,k-1)=(prsi(i,k)+prsi(i,k+1))/(prsi(i,k-1)+prsi(i,k))
96               alpha(i,k-1)=alpha(i,k-1)**kappa
97             enddo
98           enddo
99           do k=1,levs-1
100             do i=1,im
101               betta(i,k+1)=(prsi(i,k)+prsi(i,k+1))/(prsi(i,k+1)+prsi(i,k+2))
102               betta(i,k+1)=betta(i,k+1)**kappa
103             enddo
104           enddo
105           do k=1,levs
106             do i=1,im
107               gamma(i,k)= 1.0-kappa*dpp(i,k)*rpp(i,k)*2.0
108               delta(i,k)= 1.0+kappa*dpp(i,k)*rpp(i,k)*2.0
109             enddo
110           enddo
111     
112     ! ------ hybrid to solve vertical flux * 2dt ----------
113     
114           dup(:,levs)=0.0
115           dum(:,   1)=0.0
116           do k=1,levs-1
117             do i=1,im
118               dup(i,k  )=delta(i,k)*told(i,k)-betta(i,k+1)*told(i,k+1)
119               dum(i,k+1)=alpha(i,k)*told(i,k)-gamma(i,k+1)*told(i,k+1)
120             enddo
121           enddo
122     
123           k=2
124             do i=1,im
125               wmkm1=tkci(i,k)*rdelp(i,k-1)
126               wmkp1=tkci(i,k)*rdelp(i,  k)
127               wmm(i,k-1)=wmkm1*dup(i,k-1)+wmkp1*dum(i,k)-1.0
128               wmu(i,k-1)=wmkp1*dup(i,k)
129             enddo
130           do k=3,levs-1
131             do i=1,im
132               wmkm1=tkci(i,k)*rdelp(i,k-1)
133               wmkp1=tkci(i,k)*rdelp(i,  k)
134               wml(i,k-2)=wmkm1*dum(i,k-1)
135               wmm(i,k-1)=wmkm1*dup(i,k-1)+wmkp1*dum(i,k)-1.0
136               wmu(i,k-1)=wmkp1*dup(i,k)
137             enddo
138           enddo
139           k=levs
140             do i=1,im
141               wmkm1=tkci(i,k)*rdelp(i,k-1)
142               wmkp1=tkci(i,k)*rdelp(i,  k)
143               wml(i,k-2)=wmkm1*dum(i,k-1)
144               wmm(i,k-1)=wmkm1*dup(i,k-1)+wmkp1*dum(i,k)-1.0
145             enddo
146           wf(:,levs:levs+1)=0.0
147           do k=2,levs
148             do i=1,im
149               wf(i,k-1)=tkci(i,k)*(dtvdt(i,k-1)+dtvdt(i,k))
150             enddo
151           enddo
152           call tridim_hyb(im,im,levs-1,levs+1,1,
153          &                wml,wmm,wmu,wf,work,wf1)
154           wflx(:,1)=0.0
155           wflx(:,levs+1)=0.0
156           do k=2,levs
157             do i=1,im
158               wflx(i,k)=wf1(i,k-1)
159             enddo
160           enddo
161     
162     ! ------ vertical advection changes for all --------
163     c
164     c do vertical advection of tt first, since dup and dum are obtained
165     c
166           do k=1,levs
167             do i=1,im
168               zadv(i,k,3)=-rdelp (i,k)*
169          &             (wflx(i,k)*dum(i,k)+wflx(i,k+1)*dup(i,k))
170             enddo
171           enddo
172     c
173     c vertical advection of uu 
174     c
175           do k=1,levs-1
176             do i=1,im
177               dup(i,k  )=uold(i,k)-uold(i,k+1)
178               dum(i,k+1)=uold(i,k)-uold(i,k+1)
179             enddo
180           enddo
181           do k=1,levs
182             do i=1,im
183               zadv(i,k,1)=-rdelp (i,k)*
184          &             (wflx(i,k)*dum(i,k)+wflx(i,k+1)*dup(i,k))
185             enddo
186           enddo
187     c
188     c vertical advection of vv 
189     c
190           do k=1,levs-1
191             do i=1,im
192               dup(i,k  )=vold(i,k)-vold(i,k+1)
193               dum(i,k+1)=vold(i,k)-vold(i,k+1)
194             enddo
195           enddo
196           do k=1,levs
197             do i=1,im
198               zadv(i,k,2)=-rdelp (i,k)*
199          &             (wflx(i,k)*dum(i,k)+wflx(i,k+1)*dup(i,k))
200             enddo
201           enddo
202     c
203     c vertical advection of qq
204     c
205           do n=1,ntrac
206           do k=1,levs-1
207             do i=1,im
208               dup(i,k  )=rqold(i,k,n)-rqold(i,k+1,n)
209               dum(i,k+1)=rqold(i,k,n)-rqold(i,k+1,n)
210             enddo
211           enddo
212           do k=1,levs
213             do i=1,im
214               zadv(i,k,3+n)=-rdelp (i,k)*
215          &               (wflx(i,k)*dum(i,k)+wflx(i,k+1)*dup(i,k))
216             enddo
217           enddo
218           enddo
219     ! hmhj
220     ! do vertical advection filter
221           call vcnhyb(im,levs,3+ntrac,deltim,ppi,ppl,wflx,zadv,nvcn,xvcn)
222           if( nvcn.gt.0 ) print *,' in gbphys_adv nvcn=',nvcn,' xvcn=',xvcn 
223     c
224     ! add vertical filterd advection
225           do k=1,levs
226           do i=1,im
227            ug(i,k)=ug(i,k)+zadv(i,k,1)*dt2
228            vg(i,k)=vg(i,k)+zadv(i,k,2)*dt2
229            tg(i,k)=tg(i,k)+zadv(i,k,3)*dt2
230           enddo
231           enddo
232           do  n=1,ntrac
233            do k=1,levs
234            do i=1,im
235             rqg(i,k,n)=rqg(i,k,n)+zadv(i,k,3+n)*dt2
236            enddo
237            enddo
238           enddo
239     !
240     !     print *,' end of gbphys_adv_hyb_gc_fd. '			! hmhj check
241     !!
242           return
243           end
244     
245