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

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