File: C:\NOAA\NEMS_11731\src\atmos\gfs\dyn\gfidi_hyb.f

1           subroutine gfidi_hyb(lon_dim,lons_lat,lat,
2          &  dg,tg,zg,ug,vg,rqg,dphi,dlam,qg,zsphi,zslam,
3          &  rcl,spdmax,deltim,nvcn,xvcn,
4          &  dtdf,dtdl,drdf,drdl,dudl,dvdl,dudf,dvdf,
5          &  dqdt,dtdt,drdt,dudt,dvdt,szdrdt,zfirst)
6      
7     ! sela  add calculation of spdlat as in gfidiu
8     ! Fanglin Yang, June 2007, use flux-limited scheme for vertical advection of tracers
9      
10           use gfs_dyn_machine , only : kind_grid
11      
12           use namelist_dynamics_def
13           use gfs_dyn_resol_def
14           use gfs_dyn_coordinate_def	! hmhj
15           use gfs_dyn_physcons, rerth => con_rerth
16          &,             rd => con_rd, cp => con_cp
17          &,             omega => con_omega, cvap => con_cvap
18          &,             grav  => con_g
19           implicit none
20     
21           integer lon_dim,lons_lat
22           integer j,k,n,nvcn,ifirst,lat
23           real coriol,rcl,rk,sinra,deltim,xvcn,sinlat
24           real
25          1    dg(lon_dim,levs), tg(lon_dim,levs),  zg(lon_dim,levs),
26          2    ug(lon_dim,levs), vg(lon_dim,levs),
27          2   rqg(lon_dim,levs,ntrac),
28          3  dphi(lon_dim), dlam(lon_dim), qg(lon_dim),
29          3 zsphi(lon_dim),zslam(lon_dim)
30           real
31          1  dtdf(lon_dim,levs),       dtdl(lon_dim,levs),
32          1  dudf(lon_dim,levs),       dudl(lon_dim,levs),
33          1  dvdf(lon_dim,levs),       dvdl(lon_dim,levs),
34          1  drdf(lon_dim,levs,ntrac), drdl(lon_dim,levs,ntrac)
35           real
36          1  dudt(lon_dim,levs),       dvdt(lon_dim,levs),
37          1  dqdt(lon_dim),
38          1  dtdt(lon_dim,levs),
39          1  drdt(lon_dim,levs,ntrac), spdmax(levs),
40          1  szdrdt(lon_dim,levs,ntrac)    !vertical advection of tracers saved from time step n-1
41     
42           real rqg_half(lons_lat,0:levs,ntrac), rqg_d(lons_lat,0:levs,ntrac)
43           real rrkp,rrk1m,phkp,phk1m,bb,cc,tmpdrdt
44           logical zfirst
45      
46           real  pk5(lons_lat,levp1), dpk(lons_lat,levs)
47           real zadv(lons_lat,levs,3+ntrac)
48      
49           real
50          &     dot(lons_lat,levp1), dotinv(lons_lat,levp1),
51          &      ek(lons_lat,levs),      cg(lons_lat,levs),
52          &      cb(lons_lat,levs),      db(lons_lat,levs),
53          &    zlam(lons_lat,levs),    zphi(lons_lat,levs),
54          &   worka(lons_lat,levs),   workb(lons_lat,levs),
55          &   workc(lons_lat,levs),
56          &    phiu(lons_lat,levs),    phiv(lons_lat,levs),
57          &    uprs(lons_lat,levs),    vprs(lons_lat,levs),
58          &    cofa(lons_lat,levs),    cofb(lons_lat,levs),
59          &    alfa(lons_lat,levs),    rlnp(lons_lat,levs),
60          &    px1u(lons_lat,levs),    px1v(lons_lat,levs),
61          &    px2u(lons_lat,levs),    px2v(lons_lat,levs),
62          &     px2(lons_lat,levs),
63          &    px3u(lons_lat,levs),    px3v(lons_lat,levs),
64          &    px4u(lons_lat,levs),    px4v(lons_lat,levs),
65          &    px5u(lons_lat,levs),    px5v(lons_lat,levs),
66          &    uphi(lons_lat,levs),    vphi(lons_lat,levs),
67          &    expq(lons_lat),
68          &    rdel(lons_lat,levs),   rdel2(lons_lat,levs),
69          &  sumdel(lons_lat),          del(lons_lat,levs),
70          &      si(lons_lat,levp1),     sl(lons_lat,levs),
71          &    rk1,rkr,                 dif(lons_lat)
72      
73           real cons0,cons0p5,cons1,cons2,clog2   !constant
74           real rmin,rmax,delta,delta1
75           save clog2,ifirst,delta,delta1
76           data ifirst /1/
77      
78     !     print *,' enter gfidi_hyb_fd ' 		! hmhj
79     
80           cons0   = 0.d0      !constant
81           cons0p5 = 0.5d0     !constant
82           cons1   = 1.d0      !constant
83           cons2   = 2.d0      !constant
84           rk= rd /cp
85      
86     !-------------------------------------------------------
87           sinra=sqrt(cons1-cons1/rcl)     !constant
88           coriol=cons2*omega*sinra          !constant
89           sinra=sinra/rerth
90      
91           if(lat.gt.latg2) then
92           coriol=-coriol
93           sinra=-sinra
94           endif
95     !-------------------------------------------------------
96      
97     !!sela     rcl = cons1/(cons1-sinlat*sinlat)  !constant
98      
99           if(ifirst.eq.1)then
100      
101           clog2=log(cons2)     ! constant
102           delta=cvap/cp  ! check these cpv cpd (at const p for vapor and dry
103           delta1=delta-cons1
104           rk1 = rk + 1.e0
105           rkr = 1.0/rk
106      
107           ifirst=0
108           endif
109      
110           do j=1,lons_lat
111            expq(j)=exp(qg(j))
112           enddo
113     !.............................................................
114     ! get vertical coordinate for vcnhyb going bot. to top.
115           do k=1,levp1
116           do j=1,lons_lat
117            si(j,levs+2-k)= ak5(k)+bk5(k)*expq(j) !ak(k) bk(k) go top to bottom
118           enddo
119           enddo
120      
121     !sela sumdel=0.
122     !sela do k=1,levs
123     !sela   del(k)= si(k)-si(k+1)
124     !sela   sumdel=sumdel+del(k)
125     !sela enddo
126      
127           do  k=1,levs
128           do j=1,lons_lat
129     !sela  dif = si(k)**rk1 - si(k+1)**rk1
130     !sela  dif = dif / (rk1*(si(k)-si(k+1)))
131     !sela  sl(k) = dif**rkr
132            sl(j,k) = cons0p5*(si(j,k)+si(j,k+1))
133           enddo
134           enddo
135     !.............................................................
136      
137           do k=1,levp1
138           do j=1,lons_lat
139             pk5(j,k)=ak5(k) + bk5(k)*expq(j)
140     !sela  print*,'sela pk5=',pk5(k),'k=',k
141           enddo
142           enddo
143      
144           do k=1,levs
145           do j=1,lons_lat
146              dpk(j,k)=    pk5(j,k+1) - pk5(j,k)
147             rdel(j,k)=    cons1/dpk(j,k)            ! constant
148            rdel2(j,k)=  cons0p5/dpk(j,k)            ! constant
149           enddo
150           enddo
151      
152           k=1
153           do j=1,lons_lat
154            alfa(j,1)=clog2                          ! constant
155     !sela  alfa(j,1)=cons1                          ! constant
156           enddo
157      
158           do j=1,lons_lat
159             rlnp(j,1)= 99999.99
160           enddo
161           do  k=2,levs
162           do j=1,lons_lat
163             rlnp(j,k)= log( pk5(j,k+1)/pk5(j,k) )
164             alfa(j,k)= cons1-( pk5(j,k)/dpk(j,k) )*rlnp(j,k)
165           enddo
166           enddo
167      
168           spdmax=0.
169           do  k=1,levs
170           do j=1,lons_lat
171             ek(j,k)=(ug(j,levs+1-k)*ug(j,levs+1-k)+
172          1           vg(j,levs+1-k)*vg(j,levs+1-k))*rcl
173           if (ek(j,k) .gt. spdmax(levs+1-k))  spdmax(levs+1-k)=ek(j,k)
174           enddo
175           enddo
176           do k=1,levs
177           do j=1,lons_lat
178            cg(j,k)=(ug(j,levs+1-k)*dlam(j)+vg(j,levs+1-k)*dphi(j))*rcl
179           enddo
180           enddo
181           k=1
182           do j=1,lons_lat
183            db(j,1)=dg(j,levs)*dpk(j,1)
184            cb(j,1)=cg(j,1)*dbk(1)
185           enddo
186      
187           do k=1,levm1
188           do j=1,lons_lat
189            db(j,k+1)=db(j,k)+dg(j,levs-k)*dpk(j,k+1)
190            cb(j,k+1)=cb(j,k)+cg(j,k+1)*dbk(k+1)
191           enddo
192           enddo
193      
194           do j=1,lons_lat
195            dqdt(j)= -db(j,levs)/expq(j)-cb(j,levs)
196            dot(j,    1)=cons0                    !constant
197            dot(j,levp1)=cons0                    !constant
198           enddo
199           do k=1,levs-1
200           do j=1,lons_lat
201             dot(j,k+1)=-expq(j)*(bk5(k+1)*dqdt(j)+cb(j,k)) -db(j,k)
202           enddo
203           enddo
204           do k=1,levp1
205           do j=1,lons_lat
206            dotinv(j,k)=dot(j,levp1+1-k)
207           enddo
208           enddo
209      
210      
211     ! variables are in bottom to top order  !!!!!!!!!!!!!!!!!
212     ! do horizontal advection.
213      
214           k=1
215           do j=1,lons_lat
216             dudt(j,levs+1-k)=-ug(j,levs+1-k)*dudl(j,levs+1-k)
217          1                   -vg(j,levs+1-k)*dudf(j,levs+1-k)
218             dvdt(j,levs+1-k)=-ug(j,levs+1-k)*dvdl(j,levs+1-k)
219          1                   -vg(j,levs+1-k)*dvdf(j,levs+1-k)
220             dtdt(j,levs+1-k)=-ug(j,levs+1-k)*dtdl(j,levs+1-k)
221          1                   -vg(j,levs+1-k)*dtdf(j,levs+1-k)
222           enddo
223      
224           k=levs
225           do j=1,lons_lat
226             dudt(j,levs+1-k)=-ug(j,levs+1-k)*dudl(j,levs+1-k)
227          1                   -vg(j,levs+1-k)*dudf(j,levs+1-k)
228             dvdt(j,levs+1-k)=-ug(j,levs+1-k)*dvdl(j,levs+1-k)
229          1                   -vg(j,levs+1-k)*dvdf(j,levs+1-k)
230             dtdt(j,levs+1-k)=-ug(j,levs+1-k)*dtdl(j,levs+1-k)
231          1                   -vg(j,levs+1-k)*dtdf(j,levs+1-k)
232           enddo
233      
234           do k=2,levm1
235           do j=1,lons_lat
236             dudt(j,levs+1-k)=-ug(j,levs+1-k)*dudl(j,levs+1-k)
237          1                   -vg(j,levs+1-k)*dudf(j,levs+1-k)
238             dvdt(j,levs+1-k)=-ug(j,levs+1-k)*dvdl(j,levs+1-k)
239          1                   -vg(j,levs+1-k)*dvdf(j,levs+1-k)
240             dtdt(j,levs+1-k)=-ug(j,levs+1-k)*dtdl(j,levs+1-k)
241          1                   -vg(j,levs+1-k)*dtdf(j,levs+1-k)
242           enddo
243           enddo
244     ! add coriolis,deformation note  coriolis sign for s.hemi
245           do k=1,levs
246           do j=1,lons_lat
247             dudt(j,levs+1-k)=dudt(j,levs+1-k)+vg(j,levs+1-k)*coriol
248      
249             dvdt(j,levs+1-k)=dvdt(j,levs+1-k)-ug(j,levs+1-k)*coriol
250          &                   -sinra*ek(j,k)
251           enddo
252           enddo
253     !.................................................................
254     ! calculate pressure force:
255            k=1
256           do j=1,lons_lat
257             cofb(j,k)=-rdel(j,k)*(                 alfa(j,k)*dbk(k))
258           enddo
259      
260           do k=2,levs
261           do j=1,lons_lat
262             cofb(j,k)=-rdel(j,k)*(bk5(k)*rlnp(j,k)+alfa(j,k)*dbk(k))
263           enddo
264           enddo
265      
266           do k=1,levs
267           do j=1,lons_lat
268             uprs(j,k)=cofb(j,k)*rd*tg(j,levs+1-k)*expq(j)*dlam(j)
269             vprs(j,k)=cofb(j,k)*rd*tg(j,levs+1-k)*expq(j)*dphi(j)
270           enddo
271           enddo
272           do k=1,levs
273           do j=1,lons_lat
274             cofa(j,k)=-rdel(j,k)*(
275          &   bk5(k+1)*pk5(j,k)/pk5(j,k+1) - bk5(k)
276          &  +rlnp(j,k)*( bk5(k)-pk5(j,k)*dbk(k)*rdel(j,k) )  )
277           enddo
278           enddo
279     !.................................................................
280           do k=1,levs
281           do j=1,lons_lat
282             px1u(j,k)=cons0              ! grid topography =0 for testing
283             px1v(j,k)=cons0              ! grid tpopgraphy =0 for testing
284           enddo
285           enddo
286     ! see programming notes for looping in  calculating  px2 and px3
287           do j=1,lons_lat
288             px2(j,levs)=cons0                             ! constant
289             px2(j,levs-1)=
290          &  -rd*( bk5(levp1)/pk5(j,levp1)-bk5(levs)/pk5(j,levs) )*tg(j,1)
291           enddo
292      
293      
294           do k=2,levs-1
295           do j=1,lons_lat
296           px2(j,levs-k)=px2(j,levp1-k)
297          & -rd*(bk5(levs+2-k)/pk5(j,levs+2-k)-bk5(levp1-k)/pk5(j,levp1-k))*
298          &                                                  tg(j,k)
299           enddo
300           enddo
301      
302           do k=1,levs
303           do j=1,lons_lat
304             px2u(j,k)=px2(j,k)*expq(j)*dlam(j)
305             px2v(j,k)=px2(j,k)*expq(j)*dphi(j)
306           enddo
307           enddo
308           do j=1,lons_lat
309            px3u(j,levs)=cons0 ! constant
310            px3v(j,levs)=cons0 ! constant
311           enddo
312      
313           do j=1,lons_lat
314            px3u(j,levs-1)=-rd*rlnp(j,levs)*dtdl(j,1)
315            px3v(j,levs-1)=-rd*rlnp(j,levs)*dtdf(j,1)
316           enddo
317      
318           do k=2,levs-1
319           do j=1,lons_lat
320            px3u(j,levs-k)=px3u(j,levp1-k)-rd*rlnp(j,levp1-k)*dtdl(j,k)
321            px3v(j,levs-k)=px3v(j,levp1-k)-rd*rlnp(j,levp1-k)*dtdf(j,k)
322           enddo
323           enddo
324           do k=1,levs
325           do j=1,lons_lat
326            px3u(j,k)=px3u(j,k)/rcl
327            px3v(j,k)=px3v(j,k)/rcl
328           enddo
329           enddo
330           do k=1,levs
331           do j=1,lons_lat
332             px4u(j,k)=-rd*alfa(j,k)*dtdl(j,levs+1-k)/rcl
333             px4v(j,k)=-rd*alfa(j,k)*dtdf(j,levs+1-k)/rcl
334           enddo
335           enddo
336      
337           do k=1,levs
338           do j=1,lons_lat
339             px5u(j,k)=-cofa(j,k)*rd*tg(j,levs+1-k)*expq(j)*dlam(j)
340             px5v(j,k)=-cofa(j,k)*rd*tg(j,levs+1-k)*expq(j)*dphi(j)
341           enddo
342           enddo
343      
344           do k=1,levs
345           do j=1,lons_lat
346             uphi(j,k)=px1u(j,k)+px2u(j,k)+px3u(j,k)+px4u(j,k)+px5u(j,k)
347             vphi(j,k)=px1v(j,k)+px2v(j,k)+px3v(j,k)+px4v(j,k)+px5v(j,k)
348           enddo
349           enddo
350           do k=1,levs
351           do j=1,lons_lat
352             dudt(j,levs+1-k)=dudt(j,levs+1-k)+uphi(j,k)+uprs(j,k)
353             dvdt(j,levs+1-k)=dvdt(j,levs+1-k)+vphi(j,k)+vprs(j,k)
354           enddo
355           enddo
356     ! hmhj add
357           do k=1,levs
358           do j=1,lons_lat
359             dudt(j,levs+1-k)=dudt(j,levs+1-k)-grav*zslam(j)
360             dvdt(j,levs+1-k)=dvdt(j,levs+1-k)-grav*zsphi(j)
361           enddo
362           enddo
363     !..............................................................
364           do k=1,levs
365           do j=1,lons_lat
366             worka(j,k)=rk*tg(j,levs+1-k)/(cons1+(delta1)*
367          &   rqg(j,levs+1-k,1)) * rdel(j,k)
368           enddo
369           enddo
370      
371           k=1
372           do j=1,lons_lat
373            workb(j,1)=
374          &  alfa(j,1)*( dg(j,levs)*dpk(j,1)+expq(j)*cb(j,1)*dbk(1) )
375           enddo
376      
377           do k=2,levs
378           do j=1,lons_lat
379             workb(j,k)=rlnp(j,k)*( db(j,k-1)+expq(j)*cb(j,k-1) )
380          &  +alfa(j,k)*( dg(j,levs+1-k)*dpk(j,k)+expq(j)*cg(j,k)*dbk(k) )
381           enddo
382           enddo
383      
384           k=1
385           do j=1,lons_lat
386            workc(j,1)=expq(j)*cg(j,1)*dbk(1)
387           enddo
388      
389           do k=2,levs
390           do j=1,lons_lat
391             workc(j,k)=expq(j)*cg(j,k)*( dbk(k)+ck(k)*rlnp(j,k)*rdel(j,k) )
392           enddo
393           enddo
394      
395           do k=1,levs
396           do j=1,lons_lat
397            dtdt(j,levs+1-k)=
398          & dtdt(j,levs+1-k)+worka(j,k)*( -workb(j,k) + workc(j,k))
399           enddo
400           enddo
401     !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
402           do 330 n=1,ntrac
403           k=1
404           do j=1,lons_lat
405             drdt(j,levs+1-k,n)=-ug(j,levs+1-k)*drdl(j,levs+1-k,n)
406          1                     -vg(j,levs+1-k)*drdf(j,levs+1-k,n)
407           enddo
408           k=levs
409           do j=1,lons_lat
410             drdt(j,levs+1-k,n)=-ug(j,levs+1-k)*drdl(j,levs+1-k,n)
411          1                     -vg(j,levs+1-k)*drdf(j,levs+1-k,n)
412           enddo
413           do k=2,levm1
414           do j=1,lons_lat
415             drdt(j,levs+1-k,n)=-ug(j,levs+1-k)*drdl(j,levs+1-k,n)
416          1                     -vg(j,levs+1-k)*drdf(j,levs+1-k,n)
417           enddo
418           enddo
419     330   continue
420     !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
421      
422     ! do vertical advection
423      
424           k=1
425           do j=1,lons_lat
426             zadv(j,levs+1-k,1)=
427          1 -rdel2(j,k)*dot(j,k+1)*( ug(j,levs-k)-ug(j,levs+1-k))
428      
429             zadv(j,levs+1-k,2)=
430          1 -rdel2(j,k)*dot(j,k+1)*( vg(j,levs-k)-vg(j,levs+1-k))
431      
432             zadv(j,levs+1-k,3)=
433          1 -rdel2(j,k)*dot(j,k+1)*( tg(j,levs-k)-tg(j,levs+1-k))
434      
435           enddo
436      
437           k=levs
438           do j=1,lons_lat
439             zadv(j,levs+1-k,1)=
440          1 -rdel2(j,k)*dot(j,k)*( ug(j,levs+1-k)-ug(j,levs+2-k) )
441      
442             zadv(j,levs+1-k,2)=
443          1 -rdel2(j,k)*dot(j,k)*( vg(j,levs+1-k)-vg(j,levs+2-k) )
444      
445             zadv(j,levs+1-k,3)=
446          1 -rdel2(j,k)*dot(j,k)*( tg(j,levs+1-k)-tg(j,levs+2-k) )
447      
448           enddo
449      
450           do k=2,levm1
451           do j=1,lons_lat
452             zadv(j,levs+1-k,1)=
453          1 -rdel2(j,k)*( dot(j,k+1)*( ug(j,levs  -k)-ug(j,levs+1-k) ) +
454          2               dot(j,k  )*( ug(j,levs+1-k)-ug(j,levs+2-k) ) )
455      
456             zadv(j,levs+1-k,2)=
457          1 -rdel2(j,k)*( dot(j,k+1)*( vg(j,levs  -k)-vg(j,levs+1-k) ) +
458          2               dot(j,k  )*( vg(j,levs+1-k)-vg(j,levs+2-k) ) )
459      
460             zadv(j,levs+1-k,3)=
461          1 -rdel2(j,k)*( dot(j,k+1)*( tg(j,levs  -k)-tg(j,levs+1-k) ) +
462          2               dot(j,k  )*( tg(j,levs+1-k)-tg(j,levs+2-k) ) )
463      
464           enddo
465           enddo
466     !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
467     ! Fanglin Yang, June 2007
468     ! 1. use Total Variation Diminishing (TVD) flux-limited scheme
469     !    for vertical advection of tracers ( J. Thuburn, QJRMS, 1993, 469-487)
470     !    Vertical advection  dQ/dt = AA = -W*dQ/dP = -[d(Q*W)/dP-Q*dW/dP]
471     !    let BB=d(Q*W)/dP and CC=-Q*dW/dP, AA=-(BB+CC), then Q(n+1)=Q(n)+AA*dt
472     ! 2. The current scheme is central in space and central in time.  To use
473     !    the TVD scheme for vertical advection, the time differencing must be
474     !    forward in time otherwise it is unstable.  However, for horizonatl
475     !    advection which is center in space, the forward-in-time scheme is
476     !    always unstable.  To overcome this conflict, the vertical adevtion
477     !    from time step n-1 is used to get mean advection at current time
478     !    step.  Then, the central-in-time scheme is applied to both the
479     !    vertical and horizontal advections of tracers.
480     !
481           if(zflxtvd) then    !flux-limited vertical advection
482     
483           do 339 n=1,ntrac
484     
485             do j=1,lons_lat
486             do k=1,levm1            !k=1, top
487              rqg_half(j,k,n)=cons0p5*(rqg(j,levs-k,n)+rqg(j,levs+1-k,n))
488             enddo
489              rqg_half(j,0,n)=rqg(j,levs,n)
490              rqg_half(j,levs,n)=rqg(j,1,n)
491     
492             do k=1,levm1            !k=1, top
493              rqg_d(j,k,n)=rqg(j,levs-k,n)-rqg(j,levs+1-k,n)
494             enddo
495             if(rqg(j,levs,n).ge.cons0) then
496              rqg_d(j,0,n)=rqg(j,levs,n)-
497          1     max(cons0,cons2*rqg(j,levs,n)-rqg(j,levs-1,n))
498             else
499              rqg_d(j,0,n)=rqg(j,levs,n)-
500          1     min(cons0,cons2*rqg(j,levs,n)-rqg(j,levs-1,n))
501             endif
502             if(rqg(j,1,n).ge.cons0) then
503              rqg_d(j,levs,n)=max(cons0,cons2*rqg(j,1,n)-rqg(j,2,n))-
504          1     rqg(j,1,n)
505             else
506              rqg_d(j,levs,n)=min(cons0,cons2*rqg(j,1,n)-rqg(j,2,n))-
507          1     rqg(j,1,n)
508             endif
509             enddo
510     
511     ! --update tracers at half-integer layers using Van Leer (1974) limiter
512     !   (without this update, the scheme is the same as that in loop 340)
513             do j=1,lons_lat
514             do k=1,levs-1
515             if(dot(j,k+1).gt.cons0) then            !dot is from top to bottom
516              rrkp=cons0
517              if(rqg_d(j,k,n).ne.cons0) rrkp=rqg_d(j,k-1,n)/rqg_d(j,k,n)
518              phkp=(rrkp+abs(rrkp))/(1+abs(rrkp))
519              rqg_half(j,k,n)=rqg(j,levs+1-k,n)+
520          1                   phkp*(rqg_half(j,k,n)-rqg(j,levs+1-k,n))
521             else
522              rrk1m=cons0
523              if(rqg_d(j,k,n).ne.cons0) rrk1m=rqg_d(j,k+1,n)/rqg_d(j,k,n)
524              phk1m=(rrk1m+abs(rrk1m))/(1+abs(rrk1m))
525              rqg_half(j,k,n)=rqg(j,levs-k,n)+
526          1                   phk1m*(rqg_half(j,k,n)-rqg(j,levs-k,n))
527             endif
528             enddo
529             enddo
530     
531             do j=1,lons_lat
532             do k=1,levs
533              bb=rqg_half(j,k,n)*dot(j,k+1)-rqg_half(j,k-1,n)*dot(j,k)
534              cc=-rqg(j,levs+1-k,n)*(dot(j,k+1)-dot(j,k))
535              tmpdrdt=-rdel(j,k)*(bb+cc)
536              if(zfirst) then
537               zadv(j,levs+1-k,3+n)=tmpdrdt
538              else
539               zadv(j,levs+1-k,3+n)=cons0p5*(tmpdrdt+szdrdt(j,levs+1-k,n))
540              endif
541              szdrdt(j,levs+1-k,n)=tmpdrdt
542             enddo
543             enddo
544     339   continue
545     !     --------------------
546           else           !existing central-differencing
547     
548           do 340 n=1,ntrac
549           k=1
550           do j=1,lons_lat
551             zadv(j,levs+1-k,3+n)=
552          1 -rdel2(j,k)*dot(j,k+1)*( rqg(j,levs-k,n)-rqg(j,levs+1-k,n) )
553           enddo
554      
555           k=levs
556           do j=1,lons_lat
557             zadv(j,levs+1-k,3+n)=
558          1 -rdel2(j,k)*dot(j,k)*( rqg(j,levs+1-k,n)-rqg(j,levs+2-k,n) )
559           enddo
560      
561           do k=2,levm1
562           do j=1,lons_lat
563             zadv(j,levs+1-k,3+n)=-rdel2(j,k)*
564          1 ( dot(j,k+1)*( rqg(j,levs  -k,n)-rqg(j,levs+1-k,n) ) +
565          2   dot(j,k  )*( rqg(j,levs+1-k,n)-rqg(j,levs+2-k,n) ) )
566           enddo
567           enddo
568     340   continue
569     !     --------------------
570           endif
571     
572     !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
573           call vcnhyb(lons_lat,levs,3+ntrac,deltim,
574          &            si,sl,dotinv,zadv,nvcn,xvcn)
575     !sela if(xvcn.ne.0.) print*,'xvcn=',xvcn,' nvcn=',nvcn
576     !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
577      
578     ! add vertical filterd advection
579           do k=1,levs
580           do j=1,lons_lat
581            dudt(j,k)=dudt(j,k)+zadv(j,k,1)
582            dvdt(j,k)=dvdt(j,k)+zadv(j,k,2)
583            dtdt(j,k)=dtdt(j,k)+zadv(j,k,3)
584           enddo
585           enddo
586           do  n=1,ntrac
587            do k=1,levs
588            do j=1,lons_lat
589             drdt(j,k,n)=drdt(j,k,n)+zadv(j,k,3+n)
590            enddo
591            enddo
592           enddo
593      
594     !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
595     ! this multiplication must be on  completed tendencies.
596     !     do k=1,levs
597     !     do j=1,lons_lat
598     !       dudt(j,levs+1-k)=dudt(j,levs+1-k)*rcl
599     !       dvdt(j,levs+1-k)=dvdt(j,levs+1-k)*rcl
600     !     enddo
601     !     enddo
602     
603     !     print *,' leave gfidi_hyb_fd ' 		! hmhj
604     
605           return
606           end
607      
608           subroutine vcnhyb(im,km,nm,dt,zint,zmid,zdot,zadv,nvcn,xvcn)
609     c                .      .    .                                       .
610     c subprogram:    vcnhyb      vertical advection instability filter
611     c   prgmmr: iredell          org: w/nmc23    date: 91-05-07
612     c
613     c abstract: filters vertical advection tendencies
614     c   in the dynamics tendency equation in order to ensure stability
615     c   when the vertical velocity exceeds the cfl criterion.
616     c   the vertical velocity in this case is sigmadot.
617     c   for simple second-order centered eulerian advection,
618     c   filtering is needed when vcn=zdot*dt/dz>1.
619     c   the maximum eigenvalue of the linear advection equation
620     c   with second-order implicit filtering on the tendencies
621     c   is less than one for all resolvable wavenumbers (i.e. stable)
622     c   if the nondimensional filter parameter is nu=(vcn**2-1)/4.
623     c
624     c program history log:
625     c   97-07-30  iredell
626     c
627     c usage:    call vcnhyb(im,km,nm,dt,zint,zmid,zdot,zadv,nvcn,xvcn)
628     c
629     c   input argument list:
630     c     im       - integer number of gridpoints to filter
631     c     km       - integer number of vertical levels
632     c     nm       - integer number of fields
633     c     dt       - real timestep in seconds
634     c     zint     - real (im,km+1) interface vertical coordinate values
635     c     zmid     - real (im,km) midlayer vertical coordinate values
636     c     zdot     - real (im,km+1) vertical coordinate velocity
637     c     zadv     - real (im,km,nm) vertical advection tendencies
638     c
639     c   output argument list:
640     c     zadv     - real (im,km,nm) vertical advection tendencies
641     c     nvcn     - integer number of points requiring filtering
642     c     xvcn     - real maximum vertical courant number
643     c
644     c   subprograms called:
645     c     tridim_hyb   - tridiagonal matrix solver
646     c
647           implicit none
648           integer,intent(in):: im,km,nm
649           real,intent(in):: dt,zint(im,km+1),zmid(im,km),zdot(im,km+1)
650           real,intent(inout):: zadv(im,km,nm)
651           integer,intent(out):: nvcn
652           real,intent(out):: xvcn
653           integer i,j,k,n,ivcn(im)
654           logical lvcn(im)
655           real zdm,zda,zdb,vcn(im,km-1)
656           real rnu,cm(im,km),cu(im,km-1),cl(im,km-1)
657           real rr(im,km,nm)
658     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
659     c  compute vertical courant number
660     c  increase by 10% for safety
661           nvcn=0
662           xvcn=0.
663           lvcn=.false.
664           do k=1,km-1
665             do i=1,im
666               zdm=zmid(i,k)-zmid(i,k+1)
667               vcn(i,k)=abs(zdot(i,k+1)*dt/zdm)*1.1
668               lvcn(i)=lvcn(i).or.vcn(i,k).gt.1
669               xvcn=max(xvcn,vcn(i,k))
670             enddo
671           enddo
672     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
673     c  determine points requiring filtering
674           if(xvcn.gt.1) then
675             do i=1,im
676               if(lvcn(i)) then
677                 ivcn(nvcn+1)=i
678                 nvcn=nvcn+1
679               endif
680             enddo
681     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
682     c  compute tridiagonal matrim
683             do j=1,nvcn
684               cm(j,1)=1
685             enddo
686             do k=1,km-1
687               do j=1,nvcn
688                 i=ivcn(j)
689                 if(vcn(i,k).gt.1) then
690                  zdm=zmid(i,k)-zmid(i,k+1)
691                  zda=zint(i,k+1)-zint(i,k+2)
692                  zdb=zint(i,k)-zint(i,k+1)
693                   rnu=(vcn(i,k)**2-1)/4
694                   cu(j,k)=-rnu*zdm/zdb
695                   cl(j,k)=-rnu*zdm/zda
696                   cm(j,k)=cm(j,k)-cu(j,k)
697                   cm(j,k+1)=1-cl(j,k)
698                 else
699                   cu(j,k)=0
700                   cl(j,k)=0
701                   cm(j,k+1)=1
702                 endif
703               enddo
704             enddo
705     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
706     c  fill fields to be filtered
707             do n=1,nm
708               do k=1,km
709                 do j=1,nvcn
710                   i=ivcn(j)
711                   rr(j,k,n)=zadv(i,k,n)
712                 enddo
713               enddo
714             enddo
715     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
716     c  solve tridiagonal system
717             call tridim_hyb(nvcn,im,km,km,nm,cl,cm,cu,rr,cu,rr)
718     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
719     c  replace filtered fields
720             do n=1,nm
721               do k=1,km
722                 do j=1,nvcn
723                   i=ivcn(j)
724                   zadv(i,k,n)=rr(j,k,n)
725                 enddo
726               enddo
727             enddo
728           endif
729     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
730           end
731     c-----------------------------------------------------------------------
732           subroutine tridim_hyb(l,lx,n,nx,m,cl,cm,cu,r,au,a)
733     c                .      .    .                                       .
734     c subprogram:    tridim_hyb      solves tridiagonal matrix problems.
735     c   prgmmr: iredell          org: w/nmc23    date: 91-05-07
736     c
737     c abstract: this routine solves multiple tridiagonal matrix problems
738     c   with multiple right-hand-side and solution vectors for every matrix.
739     c   the solutions are found by eliminating off-diagonal coefficients,
740     c   marching first foreward then backward along the matrix diagonal.
741     c   the computations are vectorized around the number of matrices.
742     c   no checks are made for zeroes on the diagonal or singularity.
743     c
744     c program history log:
745     c   97-07-30  iredell
746     c
747     c usage:    call tridim_hyb(l,lx,n,nx,m,cl,cm,cu,r,au,a)
748     c
749     c   input argument list:
750     c     l        - integer number of tridiagonal matrices
751     c     lx       - integer first dimension (lx>=l)
752     c     n        - integer order of the matrices
753     c     nx       - integer second dimension (nx>=n)
754     c     m        - integer number of vectors for every matrix
755     c     cl       - real (lx,2:n) lower diagonal matrix elements
756     c     cm       - real (lx,n) main diagonal matrix elements
757     c     cu       - real (lx,n-1) upper diagonal matrix elements
758     c                (may be equivalent to au if no longer needed)
759     c     r        - real (lx,nx,m) right-hand-side vector elements
760     c                (may be equivalent to a if no longer needed)
761     c
762     c   output argument list:
763     c     au       - real (lx,n-1) work array
764     c     a        - real (lx,nx,m) solution vector elements
765     c
766     c attributes:
767     c   language: fortran 77.
768     c   machine:  cray.
769     c
770           real cl(lx,2:n),cm(lx,n),cu(lx,n-1),r(lx,nx,m),
771          &                         au(lx,n-1),a(lx,nx,m)
772     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
773     c  march up
774           do i=1,l
775             fk=1./cm(i,1)
776             au(i,1)=fk*cu(i,1)
777           enddo
778           do j=1,m
779             do i=1,l
780               fk=1./cm(i,1)
781               a(i,1,j)=fk*r(i,1,j)
782             enddo
783           enddo
784           do k=2,n-1
785             do i=1,l
786               fk=1./(cm(i,k)-cl(i,k)*au(i,k-1))
787               au(i,k)=fk*cu(i,k)
788             enddo
789             do j=1,m
790               do i=1,l
791                 fk=1./(cm(i,k)-cl(i,k)*au(i,k-1))
792                 a(i,k,j)=fk*(r(i,k,j)-cl(i,k)*a(i,k-1,j))
793               enddo
794             enddo
795           enddo
796     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
797     c  march down
798           do j=1,m
799             do i=1,l
800               fk=1./(cm(i,n)-cl(i,n)*au(i,n-1))
801               a(i,n,j)=fk*(r(i,n,j)-cl(i,n)*a(i,n-1,j))
802             enddo
803           enddo
804           do k=n-1,1,-1
805             do j=1,m
806               do i=1,l
807                 a(i,k,j)=a(i,k,j)-au(i,k)*a(i,k+1,j)
808               enddo
809             enddo
810           enddo
811     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
812           end
813