subroutine gfidi_sig(lon_dim,lons_lat,lat, 1 dg,tg,zg,ug,vg,rqg,dphi,dlam,qg, 1 rcl,del,rdel2,ci,tov,spdmax,deltim,sl,nvcn,xvcn, 2 dtdf,dtdl,drdf,drdl,dudl,dvdl,dudf,dvdf, 2 dqdt,dtdt,drdt,dudt,dvdt) cc USE MACHINE , ONLY : kind_phys cc use resol_def use physcons, rerth => con_rerth, rd => con_rd, cp => con_cp &, omega => con_omega implicit none cc integer j,k,n,nvcn real(kind=kind_evod) fnor,rcl,rcl2,rk,sinra,deltim,xvcn c c input variables c integer lon_dim,lat integer lons_lat c real(kind=kind_evod) 1 dg(lon_dim,levs),tg(lon_dim,levs), zg(lon_dim,levs), 2 ug(lon_dim,levs),vg(lon_dim,levs),rqg(lon_dim,levs,ntrac), 3 dphi(lon_dim),dlam(lon_dim),qg(lon_dim) c real(kind=kind_evod) 1 dtdf(lon_dim,levs),dtdl(lon_dim,levs), 1 drdf(lon_dim,levs,ntrac),drdl(lon_dim,levs,ntrac), 1 dudl(lon_dim,levs),dvdl(lon_dim,levs), 1 dudf(lon_dim,levs),dvdf(lon_dim,levs) c c output variables c real(kind=kind_evod) spdmax(levs), 1 dudt(lon_dim,levs),dvdt(lon_dim,levs), 1 dtdt(lon_dim,levs),drdt(lon_dim,levs,ntrac), 1 dqdt(lon_dim) c c constant arrays c real(kind=kind_evod) 1 sl(levs), 1 del(levs),rdel2(levs), 2 ci(levp1),tov(levs) c c local variables c real(kind=kind_evod) 1 cg (lon_dim,levs), db(lon_dim,levs),cb(lon_dim,levs), 2 dot(lon_dim,levp1),dup(lon_dim,levs),dvp(lon_dim,levs), 3 dum(lon_dim,levs ),dvm(lon_dim,levs), ek(lon_dim,levs), 4 rmu(levs ),rnu(levs),rho(levs),si(levp1), 5 x1(levs ), x2(levs), x3(levs),x4(levs) c real(kind=kind_evod) cons0,cons0p5,cons1,cons2 !constant c cons0 = 0.d0 !constant cons0p5 = 0.5d0 !constant cons1 = 1.d0 !constant cons2 = 2.d0 !constant c rk= rd /cp sinra=sqrt(cons1-cons1/rcl) !constant fnor=cons2*omega*sinra !constant sinra=sinra/rerth if(lat.gt.latg2) then fnor=-fnor sinra=-sinra endif c si(1)=cons1 !constant do 4 k=1,levs si(k+1)=si(k)-del(k) 4 continue c do 1 k=1,levm1 rho(k)= log(si(k)/si(k+1)) 1 continue rho(levs)=cons0 !constant c do 2 k=1,levs rmu(k)=cons1-si(k+1)*rho(k)/del(k) !constant 2 continue c do 3 k=1,levm1 rnu(k+1)=-cons1+si(k)*rho(k)/del(k) !constant 3 continue rnu(1)=cons0 !constant c do 20 k=1,levs x1(k)=rmu(k)*(cons1-rk*rnu(k))/(rmu(k)+rnu(k)) !constant x2(k)=cons1-x1(k) !constant x3(k)=(cons1+rk*rmu(k))/(cons1-rk*rnu(k)) !constant x4(k)=cons1/x3(k) !constant 20 continue c do 1234 k=1,levs spdmax(k)=cons0 !constant 1234 continue c rcl2=cons0p5*rcl !constant c do 140 k=1,levs do 140 j=1,lons_lat ek(j,k)=(ug(j,k)*ug(j,k)+vg(j,k)*vg(j,k))*rcl 140 continue c do 10 k=1,levs do 10 j=1,lons_lat if (ek(j,k) .gt. spdmax(k)) spdmax(k)=ek(j,k) 10 continue c c compute c=v(true)*del(ln(ps)).divide by cos for del, cos for v c do 150 j=1,lons_lat dphi(j)=dphi(j)*rcl dlam(j)=dlam(j)*rcl 150 continue c do 180 k=1,levs do 180 j=1,lons_lat cg(j,k)=ug(j,k)*dlam(j)+vg(j,k)*dphi(j) 180 continue c do 190 j=1,lons_lat db(j,1)=del(1)*dg(j,1) cb(j,1)=del(1)*cg(j,1) 190 continue c do 210 k=1,levm1 do 210 j=1,lons_lat db(j,k+1)=db(j,k)+del(k+1)*dg(j,k+1) cb(j,k+1)=cb(j,k)+del(k+1)*cg(j,k+1) 210 continue c c store integral of cg in dlax c do 220 j=1,lons_lat dqdt(j)= -cb(j,levs) 220 continue c c sigma dot computed only at interior interfaces. c do 230 j=1,lons_lat dot(j,1)=cons0 !constant dvm(j,1)=cons0 !constant dum(j,1)=cons0 !constant dot(j,levp1)=cons0 !constant dvp(j,levs )=cons0 !constant dup(j,levs )=cons0 !constant c 230 continue c do 240 k=1,levm1 do 240 j=1,lons_lat dot(j,k+1)=dot(j,k)+ 1 del(k)*(db(j,levs)+cb(j,levs)- 2 dg(j,k)-cg(j,k)) 240 continue c c c implicitly filter input profiles c if vertical advection would be numerically unstable c call vcnfil(lons_lat,lon_dim,levs,ntrac, x deltim,del,sl,si,tov,dot(1,1), x ug(1,1),vg(1,1),tg(1,1),rqg(1,1,1),nvcn,xvcn) c do 260 k=1,levm1 do 260 j=1,lons_lat dvp(j,k )=vg(j,k+1)-vg(j,k) dup(j,k )=ug(j,k+1)-ug(j,k) dvm(j,k+1)=vg(j,k+1)-vg(j,k) dum(j,k+1)=ug(j,k+1)-ug(j,k) 260 continue c do j=1,lons_lat dphi(j)=dphi(j)/rcl dlam(j)=dlam(j)/rcl enddo c do k=1,levs do j=1,lons_lat dudt(j,k)=-ug(j,k)*dudl(j,k)-vg(j,k)*dudf(j,k) 1 -rdel2(k)*(dot(j,k+1)*dup(j,k)+dot(j,k)*dum(j,k)) 2 -rd*tg(j,k)*dlam(j) c dvdt(j,k)=-ug(j,k)*dvdl(j,k)-vg(j,k)*dvdf(j,k) 1 -rdel2(k)*(dot(j,k+1)*dvp(j,k)+dot(j,k)*dvm(j,k)) 2 -rd*tg(j,k)*dphi(j) enddo enddo c do k=1,levs do j=1,lons_lat dudt(j,k)=dudt(j,k)+vg(j,k)*fnor c dvdt(j,k)=dvdt(j,k)-ug(j,k)*fnor-sinra*ek(j,k) enddo enddo c do k=1,levs do j=1,lons_lat dudt(j,k)=dudt(j,k)*rcl dvdt(j,k)=dvdt(j,k)*rcl enddo enddo c c do 280 k=1,levm1 do 280 j=1,lons_lat cecmwf: dup(j,k )=tg(j,k+1)+tov(k+1)-tg(j,k)-tov(k) x +cons2*rk*rnu(k+1)*(tg(j,k)+tov(k)) !constant cecmwf: dum(j,k+1)=tg(j,k+1)+tov(k+1)-tg(j,k)-tov(k) x +cons2*rk*rmu(k+1)*(tg(j,k+1)+tov(k+1)) !constant cecmwf: cecmwf: 280 continue c c do k=1,levs do j=1,lons_lat c dtdt(j,k)= dtdt(j,k)=-ug(j,k)*dtdl(j,k)-vg(j,k)*dtdf(j,k) 1 -rdel2(k)*(dot(j,k+1)*dup(j,k)+dot(j,k)*dum(j,k)) enddo enddo c do k=1,levs do j=1,lons_lat dtdt(j,k)=dtdt(j,k) 1 +rk*(tov(k)+tg(j,k))*(cg(j,k)-cb(j,levs)-db(j,levs)) enddo enddo c do 330 n=1,ntrac do 300 k=1,levm1 do 300 j=1,lons_lat dup(j,k )=rqg(j,k+1,n)-rqg(j,k,n) dum(j,k+1)=rqg(j,k+1,n)-rqg(j,k,n) 300 continue c do 310 j=1,lons_lat dup(j,levs)=cons0 !constant 310 continue c do 320 k=1,levs do 320 j=1,lons_lat drdt(j,k,n)=-ug(j,k)*drdl(j,k,n)-vg(j,k)*drdf(j,k,n) 1 -rdel2(k)*(dot(j,k+1)*dup(j,k)+dot(j,k)*dum(j,k)) 320 continue 330 continue c return end c----------------------------------------------------------------------- subroutine vcnfil(im,ix,km,nt,deltim,del,sl,si,tov,dot, & u,v,t,q,nvcn,xvcn) c$$$ subprogram documentation block c . . . . c subprogram: vcnfil vertical advection instability filter c prgmmr: iredell org: w/nmc23 date: 91-05-07 c c abstract: prefilters fields appearing in the nonlinear terms c in the dynamics tendency equation in order to ensure stability c when the vertical velocity exceeds the cfl criterion. c the vertical velocity in this case is sigmadot. c for simple second-order centered eulerian advection, c filtering is needed when vcn=sigmadot*deltim/dsigma>1. c the maximum eigenvalue of the linear advection equation c with second-order implicit filtering on the tendencies c is less than one for all resolvable wavenumbers (i.e. stable) c if the nondimensional filter parameter is nu=(vcn**2-1)/4. c c program history log: c 97-07-30 iredell c 98-02-20 iredell increase margin of error in filter c by starting filter when vcn>0.9. c c usage: call vcnfil(im,ix,km,nt,deltim,del,sl,si,tov,dot, c & u,v,t,q,nvcn,xvcn) c c input argument list: c im - integer number of gridpoints to filter c ix - integer first dimension of dot,u,v,t,q c km - integer number of vertical levels c nt - integer number of tracers in q c deltim - real timestep in seconds c del - real (km) sigma thicknesses c sl - real (km) full sigma values c si - real (km+1) interface sigma values c tov - real (km) temperature base c dot - real (ix,km+1) sigmadot in 1/seconds c u - real (ix,km) zonal wind c v - real (ix,km) meridional wind c t - real (ix,km) virtual temperature deviation c q - real (ix,km,nt) tracers c c output argument list: c u - real (ix,km) zonal wind c v - real (ix,km) meridional wind c t - real (ix,km) virtual temperature deviation c q - real (ix,km,nt) tracers c nvcn - integer number of points requiring filtering c xvcn - real maximum vertical courant number c c subprograms called: c tridim - tridiagonal matrix solver c c$$$ cfpp$ noconcur r cc USE MACHINE , ONLY : kind_phys use resol_def use physcons, rerth => con_rerth, rocp => con_rocp implicit none cc integer i,im,ix,j,k,km,n,nt,nvcn real(kind=kind_evod) cvcn,deli,deltim,rnu,t0term,t1term ! real(kind=kind_evod) cvcn,deli,deltim,rnu,rocp,t0term,t1term real(kind=kind_evod) tterm,xvcn cc ! parameter(rocp=rd/cp) parameter(cvcn=0.9d0) !constant real(kind=kind_evod) del(km),sl(km),si(km+1),tov(km) real(kind=kind_evod) dot(ix,km+1) real(kind=kind_evod) u(ix,km),v(ix,km),t(ix,km),q(ix,km,nt) logical lvcn(im) integer ivcn(im) real(kind=kind_evod) vcn(ix,km-1) real(kind=kind_evod) cm(ix,km),cu(ix,km-1),cl(ix,km-1) real(kind=kind_evod) rr(ix,km,3+nt) cc real(kind=kind_evod) cons0,cons0p5,cons1,cons4 !constant cc cons0 = 0.d0 !constant cons0p5 = 0.5d0 !constant cons1 = 1.d0 !constant cons4 = 4.d0 !constant cc c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c compute vertical courant number nvcn=0 xvcn=cons0 !constant do i=1,im lvcn(i)=.false. enddo do k=1,km-1 do i=1,im deli=sl(k)-sl(k+1) vcn(i,k)=abs(dot(i,k+1))*deltim/deli lvcn(i)=lvcn(i).or.vcn(i,k).gt.cvcn xvcn=max(xvcn,vcn(i,k)) enddo enddo c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c determine points requiring filtering if(xvcn.gt.cvcn) then do i=1,im if(lvcn(i)) then ivcn(nvcn+1)=i nvcn=nvcn+1 endif enddo c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c compute tridiagonal matrix do j=1,nvcn cm(j,1)=cons1 !constant enddo do k=1,km-1 deli=sl(k)-sl(k+1) do j=1,nvcn i=ivcn(j) if(vcn(i,k).gt.cvcn) then rnu=(vcn(i,k)**2-cvcn**2)/cons4 !constant cu(j,k)=-rnu*deli/del(k) cl(j,k)=-rnu*deli/del(k+1) cm(j,k)=cm(j,k)-cu(j,k) cm(j,k+1)=cons1-cl(j,k) !constant else cu(j,k)=cons0 !constant cl(j,k)=cons0 !constant cm(j,k+1)=cons1 !constant endif enddo enddo c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c fill fields to be filtered do k=1,km do j=1,nvcn i=ivcn(j) rr(j,k,1)=u(i,k) rr(j,k,2)=v(i,k) rr(j,k,3)=t(i,k)+tov(k) enddo do n=1,nt do j=1,nvcn i=ivcn(j) rr(j,k,3+n)=q(i,k,n) enddo enddo enddo c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c adjust temperature for adiabatic change do k=1,km-1 deli=sl(k)-sl(k+1) t1term=cons0p5*rocp*deli/si(k+1) !constant t0term=t1term*(tov(k)+tov(k+1)) do j=1,nvcn i=ivcn(j) tterm=t0term+t1term*(t(i,k)+t(i,k+1)) rr(j,k,3)=rr(j,k,3)-cu(j,k)*tterm rr(j,k+1,3)=rr(j,k+1,3)+cl(j,k)*tterm enddo enddo c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c solve tridiagonal system call tridim(nvcn,ix,km,km,3+nt,cl,cm,cu,rr,cu,rr) c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c replace filtered fields do k=1,km do j=1,nvcn i=ivcn(j) u(i,k)=rr(j,k,1) v(i,k)=rr(j,k,2) t(i,k)=rr(j,k,3)-tov(k) enddo do n=1,nt do j=1,nvcn i=ivcn(j) q(i,k,n)=rr(j,k,3+n) enddo enddo enddo endif c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end c----------------------------------------------------------------------- subroutine tridim(l,lx,n,nx,m,cl,cm,cu,r,au,a) c$$$ subprogram documentation block c . . . . c subprogram: tridim solves tridiagonal matrix problems. c prgmmr: iredell org: w/nmc23 date: 91-05-07 c c abstract: this routine solves multiple tridiagonal matrix problems c with multiple right-hand-side and solution vectors for every matrix. c the solutions are found by eliminating off-diagonal coefficients, c marching first foreward then backward along the matrix diagonal. c the computations are vectorized around the number of matrices. c no checks are made for zeroes on the diagonal or singularity. c c program history log: c 97-07-30 iredell c c usage: call tridim(l,lx,n,nx,m,cl,cm,cu,r,au,a) c c input argument list: c l - integer number of tridiagonal matrices c lx - integer first dimension (lx>=l) c n - integer order of the matrices c nx - integer second dimension (nx>=n) c m - integer number of vectors for every matrix c cl - real (lx,2:n) lower diagonal matrix elements c cm - real (lx,n) main diagonal matrix elements c cu - real (lx,n-1) upper diagonal matrix elements c (may be equivalent to au if no longer needed) c r - real (lx,nx,m) right-hand-side vector elements c (may be equivalent to a if no longer needed) c c output argument list: c au - real (lx,n-1) work array c a - real (lx,nx,m) solution vector elements c c attributes: c language: fortran 77. c machine: cray. c c$$$ cfpp$ noconcur r cc use machine implicit none c integer i,j,k,l,lx,m,n,nx real(kind=kind_evod) fk cc real(kind=kind_evod) cl(lx,2:n),cm(lx,n),cu(lx,n-1),r(lx,nx,m), & au(lx,n-1),a(lx,nx,m) cc real(kind=kind_evod) cons1 !constant cc cons1 = 1.d0 !constant c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c march up do i=1,l fk=cons1/cm(i,1) !constant au(i,1)=fk*cu(i,1) enddo do j=1,m do i=1,l fk=cons1/cm(i,1) !constant a(i,1,j)=fk*r(i,1,j) enddo enddo do k=2,n-1 do i=1,l fk=cons1/(cm(i,k)-cl(i,k)*au(i,k-1)) !constant au(i,k)=fk*cu(i,k) enddo do j=1,m do i=1,l fk=cons1/(cm(i,k)-cl(i,k)*au(i,k-1)) !constant a(i,k,j)=fk*(r(i,k,j)-cl(i,k)*a(i,k-1,j)) enddo enddo enddo c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c march down do j=1,m do i=1,l fk=cons1/(cm(i,n)-cl(i,n)*au(i,n-1)) !constant a(i,n,j)=fk*(r(i,n,j)-cl(i,n)*a(i,n-1,j)) enddo enddo do k=n-1,1,-1 do j=1,m do i=1,l a(i,k,j)=a(i,k,j)-au(i,k)*a(i,k+1,j) enddo enddo enddo c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end