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

1           subroutine gfidi_sig(lon_dim,lons_lat,lat,
2          1 dg,tg,zg,ug,vg,rqg,dphi,dlam,qg,zsphi,zslam,
3          1 rcl,del,rdel2,ci,tov,spdmax,deltim,sl,nvcn,xvcn,
4          2 dtdf,dtdl,drdf,drdl,dudl,dvdl,dudf,dvdf,
5          2 dqdt,dtdt,drdt,dudt,dvdt)
6     cc
7           USE gfs_dyn_MACHINE , ONLY : kind_grid
8     
9     cc
10           use gfs_dyn_resol_def
11           use gfs_dyn_physcons, rerth => con_rerth
12          &,             rd => con_rd, cp => con_cp
13          &,             omega => con_omega
14          &,             grav  => con_g
15           implicit none
16     cc
17           integer              j,k,n,nvcn
18           real(kind=kind_evod) fnor,rcl,rcl2,rk,sinra,deltim,xvcn
19     c
20     c input variables
21     c
22           integer              lon_dim,lat
23           integer              lons_lat
24     c
25           real(kind=kind_evod)
26          1  dg(lon_dim,levs),tg(lon_dim,levs), zg(lon_dim,levs),
27          2  ug(lon_dim,levs),vg(lon_dim,levs),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     c
31           real(kind=kind_evod)
32          1  dtdf(lon_dim,levs),dtdl(lon_dim,levs),
33          1  drdf(lon_dim,levs,ntrac),drdl(lon_dim,levs,ntrac),
34          1  dudl(lon_dim,levs),dvdl(lon_dim,levs),
35          1  dudf(lon_dim,levs),dvdf(lon_dim,levs)
36     c
37     c output variables
38     c
39           real(kind=kind_evod) spdmax(levs),
40          1  dudt(lon_dim,levs),dvdt(lon_dim,levs),
41          1  dtdt(lon_dim,levs),drdt(lon_dim,levs,ntrac),
42          1  dqdt(lon_dim)
43     c
44     c constant arrays
45     c
46            real(kind=kind_evod)
47          1 sl(levs),
48          1 del(levs),rdel2(levs),
49          2 ci(levp1),tov(levs)
50     c
51     c local variables
52     c
53           real(kind=kind_evod)
54          1     cg (lon_dim,levs), db(lon_dim,levs),cb(lon_dim,levs),
55          2     dot(lon_dim,levp1),dup(lon_dim,levs),dvp(lon_dim,levs),
56          3     dum(lon_dim,levs ),dvm(lon_dim,levs), ek(lon_dim,levs),
57          4     rmu(levs ),rnu(levs),rho(levs),si(levp1),
58          5      x1(levs ), x2(levs), x3(levs),x4(levs)
59     c
60           real(kind=kind_evod) cons0,cons0p5,cons1,cons2     !constant
61     c
62           cons0   = 0.d0      !constant
63           cons0p5 = 0.5d0     !constant
64           cons1   = 1.d0      !constant
65           cons2   = 2.d0      !constant
66     c
67           rk= rd /cp
68           sinra=sqrt(cons1-cons1/rcl)     !constant
69           fnor=cons2*omega*sinra          !constant
70           sinra=sinra/rerth
71      
72           if(lat.gt.latg2) then
73           fnor=-fnor
74           sinra=-sinra
75           endif
76     c
77           si(1)=cons1     !constant
78           do 4 k=1,levs
79           si(k+1)=si(k)-del(k)
80     4     continue
81     c
82           do 1 k=1,levm1
83           rho(k)= log(si(k)/si(k+1))
84     1     continue
85           rho(levs)=cons0     !constant
86     c
87           do 2 k=1,levs
88           rmu(k)=cons1-si(k+1)*rho(k)/del(k)     !constant
89     2     continue
90     c
91           do 3 k=1,levm1
92           rnu(k+1)=-cons1+si(k)*rho(k)/del(k)     !constant
93     3     continue
94           rnu(1)=cons0     !constant
95     c
96           do 20 k=1,levs
97           x1(k)=rmu(k)*(cons1-rk*rnu(k))/(rmu(k)+rnu(k))     !constant
98           x2(k)=cons1-x1(k)                                  !constant
99           x3(k)=(cons1+rk*rmu(k))/(cons1-rk*rnu(k))          !constant
100           x4(k)=cons1/x3(k)                                  !constant
101     20    continue
102     c
103           do 1234 k=1,levs
104           spdmax(k)=cons0      !constant
105     1234  continue
106     c
107           rcl2=cons0p5*rcl     !constant
108     c
109           do 140 k=1,levs
110           do 140 j=1,lons_lat
111           ek(j,k)=(ug(j,k)*ug(j,k)+vg(j,k)*vg(j,k))*rcl
112       140 continue
113     c
114           do 10 k=1,levs
115           do 10 j=1,lons_lat
116           if (ek(j,k) .gt. spdmax(k))  spdmax(k)=ek(j,k)
117        10 continue
118     c
119     c     compute c=v(true)*del(ln(ps)).divide by cos for del, cos for v
120     c
121           do 150 j=1,lons_lat
122           dphi(j)=dphi(j)*rcl
123           dlam(j)=dlam(j)*rcl
124       150 continue
125     c
126           do 180 k=1,levs
127           do 180 j=1,lons_lat
128           cg(j,k)=ug(j,k)*dlam(j)+vg(j,k)*dphi(j)
129       180 continue
130     c
131           do 190 j=1,lons_lat
132           db(j,1)=del(1)*dg(j,1)
133           cb(j,1)=del(1)*cg(j,1)
134       190 continue
135     c
136           do 210 k=1,levm1
137           do 210 j=1,lons_lat
138           db(j,k+1)=db(j,k)+del(k+1)*dg(j,k+1)
139           cb(j,k+1)=cb(j,k)+del(k+1)*cg(j,k+1)
140       210 continue
141     c
142     c   store integral of cg in dlax
143     c
144           do 220 j=1,lons_lat
145           dqdt(j)= -cb(j,levs)
146       220 continue
147     c
148     c   sigma dot computed only at interior interfaces.
149     c
150           do 230 j=1,lons_lat
151           dot(j,1)=cons0         !constant
152           dvm(j,1)=cons0         !constant
153           dum(j,1)=cons0         !constant
154           dot(j,levp1)=cons0     !constant
155           dvp(j,levs )=cons0     !constant
156           dup(j,levs )=cons0     !constant
157     c
158       230 continue
159     c
160           do 240 k=1,levm1
161           do 240 j=1,lons_lat
162           dot(j,k+1)=dot(j,k)+
163          1                 del(k)*(db(j,levs)+cb(j,levs)-
164          2                 dg(j,k)-cg(j,k))
165       240 continue
166     c
167     c
168     c  implicitly filter input profiles
169     c  if vertical advection would be numerically unstable
170     c
171           call vcnfil(lons_lat,lon_dim,levs,ntrac,
172          x            deltim,del,sl,si,tov,dot(1,1),
173          x            ug(1,1),vg(1,1),tg(1,1),rqg(1,1,1),nvcn,xvcn)
174     c
175           do 260 k=1,levm1
176           do 260 j=1,lons_lat
177           dvp(j,k  )=vg(j,k+1)-vg(j,k)
178           dup(j,k  )=ug(j,k+1)-ug(j,k)
179           dvm(j,k+1)=vg(j,k+1)-vg(j,k)
180           dum(j,k+1)=ug(j,k+1)-ug(j,k)
181       260 continue
182     c
183           do j=1,lons_lat
184            dphi(j)=dphi(j)/rcl
185            dlam(j)=dlam(j)/rcl
186           enddo
187     c
188           do k=1,levs
189            do j=1,lons_lat
190             dudt(j,k)=-ug(j,k)*dudl(j,k)-vg(j,k)*dudf(j,k)
191          1 -rdel2(k)*(dot(j,k+1)*dup(j,k)+dot(j,k)*dum(j,k))
192          2 -rd*tg(j,k)*dlam(j)
193     c
194             dvdt(j,k)=-ug(j,k)*dvdl(j,k)-vg(j,k)*dvdf(j,k)
195          1 -rdel2(k)*(dot(j,k+1)*dvp(j,k)+dot(j,k)*dvm(j,k))
196          2 -rd*tg(j,k)*dphi(j)
197            enddo
198           enddo
199     c
200           do k=1,levs
201            do j=1,lons_lat
202             dudt(j,k)=dudt(j,k)+vg(j,k)*fnor
203     c
204             dvdt(j,k)=dvdt(j,k)-ug(j,k)*fnor-sinra*ek(j,k)
205            enddo
206           enddo
207     ! hmhj add
208           do k=1,levs
209            do j=1,lons_lat
210             dudt(j,levs+1-k)=dudt(j,levs+1-k)-grav*zslam(j)
211             dvdt(j,levs+1-k)=dvdt(j,levs+1-k)-grav*zsphi(j)
212            enddo
213           enddo
214     c
215     !     do k=1,levs
216     !      do j=1,lons_lat
217     !       dudt(j,k)=dudt(j,k)*rcl
218     !       dvdt(j,k)=dvdt(j,k)*rcl
219     !      enddo
220     !     enddo
221     c
222     c
223           do 280 k=1,levm1
224           do 280 j=1,lons_lat
225     cecmwf:
226           dup(j,k  )=tg(j,k+1)+tov(k+1)-tg(j,k)-tov(k)
227          x          +cons2*rk*rnu(k+1)*(tg(j,k)+tov(k))       !constant
228     cecmwf:
229           dum(j,k+1)=tg(j,k+1)+tov(k+1)-tg(j,k)-tov(k)
230          x        +cons2*rk*rmu(k+1)*(tg(j,k+1)+tov(k+1))     !constant
231     cecmwf:
232     cecmwf:
233       280 continue
234     c
235     c
236           do k=1,levs
237            do j=1,lons_lat
238     c       dtdt(j,k)=
239             dtdt(j,k)=-ug(j,k)*dtdl(j,k)-vg(j,k)*dtdf(j,k)
240          1 -rdel2(k)*(dot(j,k+1)*dup(j,k)+dot(j,k)*dum(j,k))
241            enddo
242           enddo
243     c
244           do k=1,levs
245            do j=1,lons_lat
246             dtdt(j,k)=dtdt(j,k)
247          1  +rk*(tov(k)+tg(j,k))*(cg(j,k)-cb(j,levs)-db(j,levs))
248            enddo
249           enddo
250     c
251           do 330 n=1,ntrac
252           do 300 k=1,levm1
253           do 300 j=1,lons_lat
254           dup(j,k  )=rqg(j,k+1,n)-rqg(j,k,n)
255           dum(j,k+1)=rqg(j,k+1,n)-rqg(j,k,n)
256       300 continue
257     c
258           do 310 j=1,lons_lat
259           dup(j,levs)=cons0     !constant
260       310 continue
261     c
262           do 320 k=1,levs
263           do 320 j=1,lons_lat
264           drdt(j,k,n)=-ug(j,k)*drdl(j,k,n)-vg(j,k)*drdf(j,k,n)
265          1 -rdel2(k)*(dot(j,k+1)*dup(j,k)+dot(j,k)*dum(j,k))
266       320 continue
267       330 continue
268     c
269           return
270           end
271     c-----------------------------------------------------------------------
272           subroutine vcnfil(im,ix,km,nt,deltim,del,sl,si,tov,dot,
273          &                  u,v,t,q,nvcn,xvcn)
274     c$$$  subprogram documentation block
275     c                .      .    .                                       .
276     c subprogram:    vcnfil      vertical advection instability filter
277     c   prgmmr: iredell          org: w/nmc23    date: 91-05-07
278     c
279     c abstract: prefilters fields appearing in the nonlinear terms
280     c   in the dynamics tendency equation in order to ensure stability
281     c   when the vertical velocity exceeds the cfl criterion.
282     c   the vertical velocity in this case is sigmadot.
283     c   for simple second-order centered eulerian advection,
284     c   filtering is needed when vcn=sigmadot*deltim/dsigma>1.
285     c   the maximum eigenvalue of the linear advection equation
286     c   with second-order implicit filtering on the tendencies
287     c   is less than one for all resolvable wavenumbers (i.e. stable)
288     c   if the nondimensional filter parameter is nu=(vcn**2-1)/4.
289     c
290     c program history log:
291     c   97-07-30  iredell
292     c   98-02-20  iredell  increase margin of error in filter
293     c                      by starting filter when vcn>0.9.
294     c
295     c usage:  call vcnfil(im,ix,km,nt,deltim,del,sl,si,tov,dot,
296     c    &                u,v,t,q,nvcn,xvcn)
297     c
298     c   input argument list:
299     c     im       - integer number of gridpoints to filter
300     c     ix       - integer first dimension of dot,u,v,t,q
301     c     km       - integer number of vertical levels
302     c     nt       - integer number of tracers in q
303     c     deltim   - real timestep in seconds
304     c     del      - real (km) sigma thicknesses
305     c     sl       - real (km) full sigma values
306     c     si       - real (km+1) interface sigma values
307     c     tov      - real (km) temperature base
308     c     dot      - real (ix,km+1) sigmadot in 1/seconds
309     c     u        - real (ix,km) zonal wind
310     c     v        - real (ix,km) meridional wind
311     c     t        - real (ix,km) virtual temperature deviation
312     c     q        - real (ix,km,nt) tracers
313     c
314     c   output argument list:
315     c     u        - real (ix,km) zonal wind
316     c     v        - real (ix,km) meridional wind
317     c     t        - real (ix,km) virtual temperature deviation
318     c     q        - real (ix,km,nt) tracers
319     c     nvcn     - integer number of points requiring filtering
320     c     xvcn     - real maximum vertical courant number
321     c
322     c   subprograms called:
323     c     tridim   - tridiagonal matrix solver
324     c
325     c$$$
326     cfpp$ noconcur r
327     cc
328           USE gfs_dyn_MACHINE , ONLY : kind_grid
329     
330           use gfs_dyn_resol_def
331           use gfs_dyn_physcons, rerth => con_rerth, rocp => con_rocp
332           implicit none
333     cc
334           integer              i,im,ix,j,k,km,n,nt,nvcn
335           real(kind=kind_evod) cvcn,deli,deltim,rnu,t0term,t1term
336     !     real(kind=kind_evod) cvcn,deli,deltim,rnu,rocp,t0term,t1term
337           real(kind=kind_evod) tterm,xvcn
338     cc
339     !     parameter(rocp=rd/cp)
340           parameter(cvcn=0.9d0)     !constant
341           real(kind=kind_evod) del(km),sl(km),si(km+1),tov(km)
342           real(kind=kind_evod) dot(ix,km+1)
343           real(kind=kind_evod) u(ix,km),v(ix,km),t(ix,km),q(ix,km,nt)
344           logical              lvcn(im)
345           integer              ivcn(im)
346           real(kind=kind_evod) vcn(ix,km-1)
347           real(kind=kind_evod) cm(ix,km),cu(ix,km-1),cl(ix,km-1)
348           real(kind=kind_evod) rr(ix,km,3+nt)
349     cc
350           real(kind=kind_evod) cons0,cons0p5,cons1,cons4     !constant
351     cc
352           cons0   =  0.d0         !constant
353           cons0p5 = 0.5d0         !constant
354           cons1   =  1.d0         !constant
355           cons4   =  4.d0         !constant
356     cc
357     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
358     c  compute vertical courant number
359           nvcn=0
360           xvcn=cons0     !constant
361           do i=1,im
362             lvcn(i)=.false.
363           enddo
364           do k=1,km-1
365             do i=1,im
366               deli=sl(k)-sl(k+1)
367               vcn(i,k)=abs(dot(i,k+1))*deltim/deli
368               lvcn(i)=lvcn(i).or.vcn(i,k).gt.cvcn
369               xvcn=max(xvcn,vcn(i,k))
370             enddo
371           enddo
372     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
373     c  determine points requiring filtering
374           if(xvcn.gt.cvcn) then
375             do i=1,im
376               if(lvcn(i)) then
377                 ivcn(nvcn+1)=i
378                 nvcn=nvcn+1
379               endif
380             enddo
381     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
382     c  compute tridiagonal matrix
383             do j=1,nvcn
384               cm(j,1)=cons1                           !constant
385             enddo
386             do k=1,km-1
387               deli=sl(k)-sl(k+1)
388               do j=1,nvcn
389                 i=ivcn(j)
390                 if(vcn(i,k).gt.cvcn) then
391                   rnu=(vcn(i,k)**2-cvcn**2)/cons4     !constant
392                   cu(j,k)=-rnu*deli/del(k)
393                   cl(j,k)=-rnu*deli/del(k+1)
394                   cm(j,k)=cm(j,k)-cu(j,k)
395                   cm(j,k+1)=cons1-cl(j,k)             !constant
396                 else
397                   cu(j,k)=cons0                       !constant
398                   cl(j,k)=cons0                       !constant
399                   cm(j,k+1)=cons1                     !constant
400                 endif
401               enddo
402             enddo
403     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
404     c  fill fields to be filtered
405             do k=1,km
406               do j=1,nvcn
407                 i=ivcn(j)
408                 rr(j,k,1)=u(i,k)
409                 rr(j,k,2)=v(i,k)
410                 rr(j,k,3)=t(i,k)+tov(k)
411               enddo
412               do n=1,nt
413                 do j=1,nvcn
414                   i=ivcn(j)
415                   rr(j,k,3+n)=q(i,k,n)
416                 enddo
417               enddo
418             enddo
419     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
420     c  adjust temperature for adiabatic change
421             do k=1,km-1
422               deli=sl(k)-sl(k+1)
423               t1term=cons0p5*rocp*deli/si(k+1)     !constant
424               t0term=t1term*(tov(k)+tov(k+1))
425               do j=1,nvcn
426                 i=ivcn(j)
427                 tterm=t0term+t1term*(t(i,k)+t(i,k+1))
428                 rr(j,k,3)=rr(j,k,3)-cu(j,k)*tterm
429                 rr(j,k+1,3)=rr(j,k+1,3)+cl(j,k)*tterm
430               enddo
431             enddo
432     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
433     c  solve tridiagonal system
434             call tridim(nvcn,ix,km,km,3+nt,cl,cm,cu,rr,cu,rr)
435     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
436     c  replace filtered fields
437             do k=1,km
438               do j=1,nvcn
439                 i=ivcn(j)
440                 u(i,k)=rr(j,k,1)
441                 v(i,k)=rr(j,k,2)
442                 t(i,k)=rr(j,k,3)-tov(k)
443               enddo
444               do n=1,nt
445                 do j=1,nvcn
446                   i=ivcn(j)
447                   q(i,k,n)=rr(j,k,3+n)
448                 enddo
449               enddo
450             enddo
451           endif
452     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
453           end
454     c-----------------------------------------------------------------------
455           subroutine tridim(l,lx,n,nx,m,cl,cm,cu,r,au,a)
456     c$$$  subprogram documentation block
457     c                .      .    .                                       .
458     c subprogram:    tridim      solves tridiagonal matrix problems.
459     c   prgmmr: iredell          org: w/nmc23    date: 91-05-07
460     c
461     c abstract: this routine solves multiple tridiagonal matrix problems
462     c   with multiple right-hand-side and solution vectors for every matrix.
463     c   the solutions are found by eliminating off-diagonal coefficients,
464     c   marching first foreward then backward along the matrix diagonal.
465     c   the computations are vectorized around the number of matrices.
466     c   no checks are made for zeroes on the diagonal or singularity.
467     c
468     c program history log:
469     c   97-07-30  iredell
470     c
471     c usage:    call tridim(l,lx,n,nx,m,cl,cm,cu,r,au,a)
472     c
473     c   input argument list:
474     c     l        - integer number of tridiagonal matrices
475     c     lx       - integer first dimension (lx>=l)
476     c     n        - integer order of the matrices
477     c     nx       - integer second dimension (nx>=n)
478     c     m        - integer number of vectors for every matrix
479     c     cl       - real (lx,2:n) lower diagonal matrix elements
480     c     cm       - real (lx,n) main diagonal matrix elements
481     c     cu       - real (lx,n-1) upper diagonal matrix elements
482     c                (may be equivalent to au if no longer needed)
483     c     r        - real (lx,nx,m) right-hand-side vector elements
484     c                (may be equivalent to a if no longer needed)
485     c
486     c   output argument list:
487     c     au       - real (lx,n-1) work array
488     c     a        - real (lx,nx,m) solution vector elements
489     c
490     c attributes:
491     c   language: fortran 77.
492     c   machine:  cray.
493     c
494     c$$$
495     cfpp$ noconcur r
496     cc
497           use gfs_dyn_machine
498           implicit none
499     c
500           integer              i,j,k,l,lx,m,n,nx
501           real(kind=kind_evod) fk
502     cc
503           real(kind=kind_evod) cl(lx,2:n),cm(lx,n),cu(lx,n-1),r(lx,nx,m),
504          &                                         au(lx,n-1),a(lx,nx,m)
505     cc
506           real(kind=kind_evod) cons1                        !constant
507     cc
508           cons1   =  1.d0                                   !constant
509     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
510     c  march up
511           do i=1,l
512             fk=cons1/cm(i,1)                                !constant
513             au(i,1)=fk*cu(i,1)
514           enddo
515           do j=1,m
516             do i=1,l
517               fk=cons1/cm(i,1)                              !constant
518               a(i,1,j)=fk*r(i,1,j)
519             enddo
520           enddo
521           do k=2,n-1
522             do i=1,l
523               fk=cons1/(cm(i,k)-cl(i,k)*au(i,k-1))          !constant
524               au(i,k)=fk*cu(i,k)
525             enddo
526             do j=1,m
527               do i=1,l
528                 fk=cons1/(cm(i,k)-cl(i,k)*au(i,k-1))        !constant
529                 a(i,k,j)=fk*(r(i,k,j)-cl(i,k)*a(i,k-1,j))
530               enddo
531             enddo
532           enddo
533     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
534     c  march down
535           do j=1,m
536             do i=1,l
537               fk=cons1/(cm(i,n)-cl(i,n)*au(i,n-1))          !constant
538               a(i,n,j)=fk*(r(i,n,j)-cl(i,n)*a(i,n-1,j))
539             enddo
540           enddo
541           do k=n-1,1,-1
542             do j=1,m
543               do i=1,l
544                 a(i,k,j)=a(i,k,j)-au(i,k)*a(i,k+1,j)
545               enddo
546             enddo
547           enddo
548     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
549           end
550