program phase c$$$ main program documentation block c c main program: gfdc_chgvortx c prgmmr: marchok org: np22 date: 2001-03-29 c c abstract: this program adds the bogus storm to the environmental c field. it first removes the analyzed vortex from the global c analysis using the variable rnot filter. we input the forecasted c asymmetries from the previous 12-hour forecast if they exist. c finally the symmetric and asymmetric part of the vortex is added back c into the environmental field at the correct storm location. c c program history log: c 95-05-02 bender c 98-05-07 bender c c input files: c unit 9 the axi-symmetric vortex determined from the axi-sym. mo c unit 88 the various inputs to bogus (rnot is inputed in this ste c unit 55 the inital condition of the nested analysis c c output files: c the 77 the initial condition with the specified vortex added in c c attributes: c language: fortran c c$$$ c cc*********************************************************************** cc cc this subroutine obtains the basic flow field for u,v,t,r,p-star cc field for each nest. the value is first obtained at the cc coarse resolution. values for the finer meshes are then cc obtained by interpolation cc cc cc*********************************************************************** parameter (imx=150,jmx=150,kmax=42,lgi=16,kmaxp=kmax+1) parameter (n = 40, iblkmx = lgi*imx+4*kmax*imx, nd = 20, * jblkmx = imx + kmax*imx*14, nnst = 3, ngdfd = 19) parameter (ixm = 99, iym = 99 , kmaxm=kmax-1, * nxm = 99, nym = 99 , nh = 30,nmx=24) cc cc cc*********************************************************************** cc cc cc ixm = greater then compound columns of the inner mesh cc iym = greater then compound rows of the inner mesh cc nxm = greater then single columns of the inner mesh cc nym = greater then single rows of the inner mesh cc cc*********************************************************************** cc cc parameter (ijm=ixm*iym, njm=nxm*nym,ikx=ijm*1000) common /vect/ro(nmx),xvect(nmx),yvect(nmx) common /consts/ rgas,cappa,tdt,q(kmax),dq(kmax),odq(kmax), * odq1(kmax),dqroq(kmax),qhsm(kmax),qmh(kmaxp),dqph(kmaxm), * dkmh(kmaxp),og,o2g,ag,omag,cpg,cgc,ctque common /cqlgcc/ qlgc1(kmax),qlgc2(kmax),qlgc3(kmax),qlgc4(kmax), * qp(4),nqp(4),drqhcp(kmax),dqrqk(kmax),rgash * ,aorpad5(7) common /posit/ xold,yold,xcorn,ycorn,rox,xv,yv,rb,ihour common /bkinfo/ nstep,nnest,iblock common /xxx/ dummy(imx,jmx),xxcc,yycc,ddxx,ddyy common /hegt/ phi(kmax),ktop common /filec/ typ1c(imx,lgi),typ2c(kmax,imx,4) common /gdinf/ ngd,ngr,ntr,dt,js,jn,ie,iw,iimax,imax,jjmax, * jmax,nstflg,icx,icy,ihx,ihy,dftx,dfty common /valu/ wti(nxm),wtj(nym),idi(nxm),idj(nym) common /compnd/ curr(iblkmx,nnst),diag(jblkmx,nnst) common /ibnd/ jsg(nnst),jng(nnst),iwg(nnst),ieg(nnst), * imag(nnst),jmag(nnst),jjmag(nnst),ngrg(nnst) common / vtemp/ u(imx,jmx,kmax),v(imx,jmx,kmax), * t(imx,jmx,kmax),r(imx,jmx,kmax), * del(imx,jmx,nnst),tha(imx,jmx,nnst), * pstar(imx,jmx),xyd(2,ijm), * xd(ijm),yd(ijm),zd(ijm),zim(njm), * xim(nxm),yim(nym),xps(imx,nh),yps(jmx,nh), * xzs(imx,nh),yzs(jmx,nh), * ytu(kmax,jmx,nh),ytv(kmax,jmx,nh),ytt(kmax,jmx,nh), * xtu(kmax,imx,nh),xtv(kmax,imx,nh),xtt(kmax,imx,nh), * ytr(kmax,jmx,nh),xtr(kmax,imx,nh),pmm(imx,jmx,2) common /wnda/ ua(imx,jmx),va(imx,jmx) common /tang/ tanw(kmax,imx) common /ang/ delg(imx,nnst),thag(jmx,nnst),ddelg(nnst), * dthag(nnst) common /total/ uh(imx,jmx),vh(imx,jmx),ds(imx,jmx),dl(imx,jmx), * th(imx,jmx),f(imx,jmx),ddel,dtha,dzbdx,dzbdy dimension phih(kmaxp),tbrr(kmax) dimension qnw2(kmax),qnw1(kmax) dimension hgc(kmax+1,imx),gc(kmax,imx,imx) dimension vortx(kmax,imx,4),psvr(imx) dimension uu(kmax,jmx),vv(kmax,jmx),aswnd(kmax,imx,2) dimension tdinf(19),filepc(iblkmx),nnrg(nnst),chc(imx,3) dimension cpa(imx,jmx,kmax,4), tk(nh) , ampf(25),pp(imx,jmx) dimension uvp(imx,jmx,kmax,4),xxd(imx,jmx),prime(kmax,imx,2) dimension frmn(imx,4,kmax,4),frme(jmx,4,kmax,4),pspr(imx) dimension fpsn(imx,4), fpse(jmx,4), pert(imx,3) dimension fact1(imx),fact2(imx) dimension u1(imx),v1(imx),t1(imx),r1(imx) cc common /worksp/ rwksp real rwksp cc data tk/.25,.3333333,.5,.7236068,1.0,1.3279853, cc * 1.7071068,2.1371581/ equivalence (filepc(1),typ1c(1,1)),(u(1,1,1),cpa(1,1,1,1)) character*4 iblock cc cc cc cc***************************************************************** cc * cc unit of the original history tape field is : 55 * cc unit of the final filtered history tape field is : 77 * cc be sure to always have enough disc reserved for these units * cc * cc * cc**************************************************************** cc cc cc***************************************************************** cc * cc * external star cc cc determine the positions of the old vortex (xold,yold) * cc in degrees * cc dexf = 1./12. pi180 = 4.*atan(1.0)/180. cc dggg=0.0 ifric=0 cc************************************************************* c c c modification of moisture paramater as a function of intensity c and change in storm intensity c read(21,*)ipold read(22,*)ipcur ipcurt = ipcur pcurt = float(ipcurt) if(ipold .lt. 0) then pold=1000. pcur=1000. else pold = float(ipold) pcur = float(ipcur) print* print*,'previous minimum pressure: ',pold print*,'current minimum pressure: ',pcur print* endif c hcor=.10 c if(pcur .ge. 960. .or. pcur .le. 985.) then corre = (985. - pcur)/25. tints = hcor * corre endif if(pcur .gt. 985.) then tints = 0.0 endif if(pcur .lt. 960.) then tints = hcor endif pbase = .85 + tints alt = .04 pdiff = pold - pcur c fctf = pbase + pdiff * alt if(fctf .lt. .60) then fctf = .60 endif if(fctf .gt. 1.0) then fctf = 1.0 endif c c fctf = 1.0 c c pcur = pcurt print*,'current pressure of storm: ',pcurt c print*,'6 hour pressure change (hpa): ',pdiff print* print*,'moisture factor over water: ',fctf c read(88)xv,yv write(6,*)' write the values of xv,yv',xv,yv read(88)iqd if(iqd.eq.4)then read(88)wne,wse,wsw,wnw endif read(88)tjunk1,tjunk2,tjunk3 cc read(88)hurrr,deep cc read(88)pmax read(88)ihour print* print*,'this is ihour: ',ihour print* read(88)tsize cc read(88)xold,yold,xcorn,ycorn print* print*,'(xold,yold) from disc: ',xold,yold print*,'(xcorn,ycorn) from disc: ',xcorn,ycorn print* read(88)ro,rb c cosft = cos(yold*pi180) print* print*,'(ro,rb) from disc: ',ro,rb print* c rewind 88 cc delc = xv*pi180 thac = yv*pi180 cc******************************************************************* cc******************************************************************* cc cc note: this version is assuming that the outer gird has cc a grid resolution of 1 degree and the innermost grid cc has a grid resolution of 1/6 degree. this is manifested cc in the routine ...init... in the variable ....cord... cc which currently is defined as 1./12. - .5. cc this will have to be changed when the resolutions are cc changed. in the future this will be made more general!!!!!! cc cc********************************************************* cc ccxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx tlap = 6.7e-5/(980.6) cc call const cc do 3346 id =1 , 4 do 3346 k = 1 , kmax do 3346 j = 1 , jmx do 3346 i = 1 , imx cpa(i,j,k,id) = 0.0 3346 continue cc qnw1(1) = ALOG(qmh(2)) - ALOG(q(1)) do k = 2,kmax qnw1(k) = ALOG(qmh(k+1)) - ALOG(q(k)) qnw2(k) = ALOG(qmh(k+1)) - ALOG(qmh(k)) enddo cc do 3349 j = 1 , jmx do 3349 i = 1 , imx pstar(i,j) = 0.0 3349 continue cc cc read(55)nstep,nnest,iblock cc cc tn = float(n) ccc pi = 4.*atan(1.0) cosf = cos(2.*pi/tn) - 1. cc nnrg(1) = 0 cc do 10 ng = 1 , nnest cc read(55)ngd,ngr,ntr,dt,js,jn,ie,iw,iimax,imax,jjmax,jmax, * nstflg,ickx,icky,ihx,ihy,dftx,dfty ccc jsg(ng) = js jng(ng) = jn ieg(ng) = ie iwg(ng) = iw imag(ng) = imax jmag(ng) = jmax jjmag(ng) = jjmax ngrg(ng) = ngr cc cc imaxm = imax - 1 jmaxm = jmax - 1 lngth1 = iblkmx*ngr j = 0 do 15 jj = 1 , jjmax cc read(55)diggg read(55)((curr(i,jrt),i=1,iblkmx),jrt=1,ngr) cc cc njrow = ngr if(jj.le.2.or.jj.ge.jjmax-1)njrow = 1 do 18 jr = 1 , njrow j = j + 1 do 22 i = 1 , iblkmx 22 filepc(i) = curr(i,jr) do 23 i = 1 , imax del(i,j,ng) = typ1c(i,5) tha(i,j,ng) = typ1c(i,6) 23 continue cc thag(j,ng) = typ1c(3,6) cc cc*************************************************************** cc if(ng.eq.1)then cc do 29 i = 1 , imx dl(i,j) = del(i,j,1) th(i,j) = tha(i,j,1) ds(i,j) = typ1c(i,4) f(i,j) = typ1c(i,9) 29 continue cc endif cc******************************* cccc if(j.lt.6.and.jr.eq.njrow)then do 2600 i = 1 , imx delg(i,ng) = typ1c(i,5) 2600 continue dthag(ng) = typ1c(3,3) ddelg(ng) = typ1c(3,2) cc cc cc define other constants for assymetric wind if(ng.eq.1)then dtha = typ1c(3,3) ddel = typ1c(3,2) cc ddxx = cosft*ddel/pi180 ddyy = dtha/pi180 xxcc = (xold-xcorn)*ddxx yycc = (yold-ycorn)*ddyy cc write(6,3111)dtha,ddel,ddxx,ddyy,pi180,cosft 3111 format(2x,6e11.5) cc endif cc********************************* cc endif cc cc**************************************************************** cc if(ng.eq.1)then do i = 1 , imax hgc(kmax+1,i) = typ1c(i,11) do k = kmax,2,-1 hgc(k,i) = hgc(k+1,i) + rgas*typ2c(k,i,3)*qnw2(k) enddo do k = kmax,1,-1 gc(k,i,j) = hgc(k+1,i) + rgas*typ2c(k,i,3)*qnw1(k) enddo enddo do 2311 i = 1 , imax tbar = alog(1.0+tlap*typ1c(i,11)/typ2c(kmax,i,3)) epsi = tbar/(rgas*tlap) psea = typ1c(i,lgi)*exp(epsi) typ1c(i,lgi) = psea 2311 continue cc do 740 id = 1 , 4 do 740 i = 1 , imx do 740 k = 1 , kmax uvp(i,j,k,id) = typ2c(k,i,id) 740 continue do i = 1 , imx do k = 1 , kmax t(i,j,k) = typ2c(k,i,3) + tlap*gc(k,i,j) enddo enddo do 741 i = 1 , imx pp(i,j) = typ1c(i,lgi) 741 continue endif cc 18 continue 15 continue 10 continue cc imax = imag(1) jmax = jmag(1) imaxm = imax - 1 jmaxm = jmax - 1 ccc print*,'imax for smoothing: ',imax,jmax,ktop rewind 15 read(15,*)rfavg print*,'this is rfavg: ',rfavg rewind 55 rfav2 = 5.0 rfav3 = 5.6 cc begin to do the smoothing ....first in the east-west direction cc if(pcur .lt. 1000. )then bbot = .120 bkt2 = .590 ccc kbot =12 ccc kt2 = 25 else bbot = .170 bkt2 = .675 ccc kbot =14 ccc kt2 = 28 endif c print* do 2284 k=1,kmax if(q(k).gt.bbot)go to 2011 2284 continue 2011 continue kbot = k c do 2294 k=1,kmax if(q(k).gt.bkt2)go to 2091 2294 continue 2091 continue kt2 = k print*,'smoothing upper bounds: ',kbot,kt2 c do k = 1 , kmax c do j = 1 , jmx c do i = 1 , imx c t(i,j,k) = uvp(i,j,k,3) c enddo c enddo c enddo do k = kbot + 1 , kmax if(k .ge. kbot .and. k .le. kt2)then ifl = 2 endif if(k .ge. kt2+1 .and. k .le. ktop-1)then ifl = 2 if(rfavg .gt. rfav3 .or. pcur .lt. 981.0) ifl = 3 endif if(k .ge. ktop)then ifl = 2 if(rfavg .gt. rfav2 .or. pcur .lt. 993.0) then ifl = 3 endif if(rfavg .gt. rfav3 .or. pcur .lt. 971.0) then ifl = 4 endif endif c c sandy change ifl = 4 c print*,'smoothing stength: ',k,ifl cc chg = 0.0 kt = 0 if(ifl.eq.1)nty = 15 c 2009 test if(ifl.eq.3)nty = 19 c 2009 test if(ifl.eq.4)nty = 21 if(ifl.eq.2)nty = 17 if(ifl.eq.3)nty = 22 if(ifl.eq.4)nty = 25 print*,'this is nty: ',nty do 802 kty = 1 , nty if((kty.eq.4.or.kty.eq.8.or.kty.ge.12) * .and.ifl.eq.1)chg = 1.0 if((kty.eq.4.or.kty.eq.8.or.kty.eq.12 * .or.kty.ge.14).and.ifl.eq.2)chg = 1.0 c c 2009 test if((kty.eq.4.or.kty.eq.8.or.kty.eq.12 c * .or.kty.ge.15).and.ifl.eq.3)chg = 1.0 c 2009 test if((kty.eq.4.or.kty.eq.8.or.kty.eq.12 c * .or.kty.ge.16).and.ifl.eq.4)chg = 1.0 c if((kty.eq.4.or.kty.eq.8.or.kty.eq.12 * .or.kty.ge.16).and.ifl.eq.3)chg = 1.0 if((kty.eq.4.or.kty.eq.8.or.kty.eq.12 * .or.kty.eq.16 .or. kty .ge. 18).and.ifl.eq.4) * chg = 1.0 cc if(chg.eq.0)kt = kt + 1 if(chg.eq.1.0)tk(kty) = .25 if(chg.eq.1.0)go to 801 fact = 2.0*pi/(float(kt) + 1.0) tk(kty) = -.5/(cos(fact) - 1.0) do 679 na = 2 , 25 ampf(na) = 1 + 2.*tk(kty)*(cos(2.*pi/float(na)) - 1.0) 679 continue c 801 continue chg = 0.0 802 continue write(6,815) (tk(kk),kk = 1 , nty) 815 format(2x,'this is tk:',e12.6) do j = 1 , jmax cc do 70 i = 1 , imax u1(i) = uvp(i,j,k,1) v1(i) = uvp(i,j,k,2) t1(i) = t(i,j,k) cccc t1(i) = uvp(i,j,k,3) r1(i) = uvp(i,j,k,4) 70 continue cc do 588 nn = 1 , nty xps(1,nn) = pp(1,j) xps(imax,nn) = pp(imax,j) 588 continue cc do 58 nn = 1 , nty xtu(k,1,nn) = u1(1) xtu(k,imax,nn) = u1(imax) xtv(k,1,nn) = v1(1) xtv(k,imax,nn) = v1(imax) xtt(k,1,nn) = t1(1) xtt(k,imax,nn) = t1(imax) xtr(k,1,nn) = r1(1) xtr(k,imax,nn) = r1(imax) 58 continue ccc do 660 i = 2 , imaxm xps(i,1) = pp(i,j) + tk(1)*(pp(i-1,j) + * pp(i+1,j) - 2.*pp(i,j)) 660 continue cc do 60 i = 2 , imaxm xtu(k,i,1) = u1(i) + tk(1)*(u1(i-1) + * u1(i+1) - 2.*u1(i)) xtv(k,i,1) = v1(i) + tk(1)*(v1(i-1) + * v1(i+1) - 2.*v1(i)) xtt(k,i,1) = t1(i) + tk(1)*(t1(i-1) + * t1(i+1) - 2.*t1(i)) xtr(k,i,1) = r1(i) + tk(1)*(r1(i-1) + * r1(i+1) - 2.*r1(i)) 60 continue cc do 650 nn = 2 , nty do 620 i = 2 , imaxm xps(i,nn) = xps(i,nn-1) + tk(nn)*(xps(i-1,nn-1) + * xps(i+1,nn-1) - 2.*xps(i,nn-1)) 620 continue 650 continue cc do 65 nn = 2 , nty do 62 i = 2 , imaxm xtu(k,i,nn) = xtu(k,i,nn-1) + tk(nn)*(xtu(k,i-1,nn-1) + *xtu(k,i+1,nn-1) - 2.*xtu(k,i,nn-1)) xtv(k,i,nn) = xtv(k,i,nn-1) + tk(nn)*(xtv(k,i-1,nn-1) + *xtv(k,i+1,nn-1) - 2.*xtv(k,i,nn-1)) xtt(k,i,nn) = xtt(k,i,nn-1) + tk(nn)*(xtt(k,i-1,nn-1) + *xtt(k,i+1,nn-1) - 2.*xtt(k,i,nn-1)) xtr(k,i,nn) = xtr(k,i,nn-1) + tk(nn)*(xtr(k,i-1,nn-1) + *xtr(k,i+1,nn-1) - 2.*xtr(k,i,nn-1)) 62 continue 65 continue cc do i = 1 , imax u(i,j,k) = xtu(k,i,nty) v(i,j,k) = xtv(k,i,nty) t(i,j,k) = xtt(k,i,nty) r(i,j,k) = xtr(k,i,nty) pstar(i,j) = xps(i,nty) enddo cc cc end of j loop cc enddo cc now do the smoothing in the meridional direction cc do 85 i = 1 , imax do 80 nn = 1 , nty yps(1,nn) = pstar(i,1) yps(jmax,nn) = pstar(i,jmax) ytu(k,1,nn) = u(i,1,k) ytu(k,jmax,nn) = u(i,jmax,k) ytv(k,1,nn) = v(i,1,k) ytv(k,jmax,nn) = v(i,jmax,k) ytt(k,1,nn) = t(i,1,k) ytt(k,jmax,nn) = t(i,jmax,k) ytr(k,1,nn) = r(i,1,k) ytr(k,jmax,nn) = r(i,jmax,k) 80 continue cc cc cc do 900 j = 2 , jmaxm yps(j,1) = pstar(i,j) + tk(1)*(pstar(i,j-1) + pstar(i,j+1) * -2.*pstar(i,j)) 900 continue cc cc cc do 90 j = 2 , jmaxm ytu(k,j,1) = u(i,j,k) + tk(1)*(u(i,j-1,k) + u(i,j+1,k) * -2.*u(i,j,k)) ytv(k,j,1) = v(i,j,k) + tk(1)*(v(i,j-1,k) + v(i,j+1,k) * -2.*v(i,j,k)) ytt(k,j,1) = t(i,j,k) + tk(1)*(t(i,j-1,k) + t(i,j+1,k) * -2.*t(i,j,k)) ytr(k,j,1) = r(i,j,k) + tk(1)*(r(i,j-1,k) + r(i,j+1,k) * -2.*r(i,j,k)) 90 continue cc cc do 950 nn = 2 , nty do 950 j = 2 , jmaxm yps(j,nn) = yps(j,nn-1) + tk(nn)*(yps(j-1,nn-1) + * yps(j+1,nn-1) - 2.*yps(j,nn-1)) 950 continue cc cc do 95 nn = 2 , nty do 95 j = 2 , jmaxm ytu(k,j,nn) = ytu(k,j,nn-1) + tk(nn)*(ytu(k,j-1,nn-1) + * ytu(k,j+1,nn-1) - 2.*ytu(k,j,nn-1)) ytv(k,j,nn) = ytv(k,j,nn-1) + tk(nn)*(ytv(k,j-1,nn-1) + * ytv(k,j+1,nn-1) - 2.*ytv(k,j,nn-1)) ytt(k,j,nn) = ytt(k,j,nn-1) + tk(nn)*(ytt(k,j-1,nn-1) + * ytt(k,j+1,nn-1) - 2.*ytt(k,j,nn-1)) ytr(k,j,nn) = ytr(k,j,nn-1) + tk(nn)*(ytr(k,j-1,nn-1) + * ytr(k,j+1,nn-1) - 2.*ytr(k,j,nn-1)) 95 continue cc do 99 j = 1 , jmax pstar(i,j) = yps(j,nty) u(i,j,k) = ytu(k,j,nty) v(i,j,k) = ytv(k,j,nty) t(i,j,k) = ytt(k,j,nty) r(i,j,k) = ytr(k,j,nty) 99 continue cc 85 continue cc cc end of k loop cc enddo cc cc cc pstar(i,j) is the smoothed surface pressure of basic field of nst cc u(i,j,k) is the smoothed u component of the basic field of nest cc v(i,j,k) is the smoothed v component of the basic field of nest cc t(i,j,k) is the smoothed t component of the basic field of nest cc r(i,j,k) is the smoothed r component of the basic field of nest cc this is the estimated phase speed for nest 1 so write it out to di cc cc cc write the smoothed values out to disc files to be later written to cc the history tape..... cc tmx = float(imx)*float(jmx) do 2377 k = 1 , kmax tbrr(k) = 0.0 2377 continue do 2301 j = 1 , jmax do 2301 i = 1 , imax do 2301 k = 1 , kmax t(i,j,k) = t(i,j,k) - tlap*gc(k,i,j) tbrr(k) = tbrr(k) + t(i,j,k)/tmx 2301 continue phih(kmaxp) = 0.0 do 2302 kk = 1 , kmax k = kmax - kk + 1 phih(k) = phih(k+1) + qlgc3(k)*tbrr(k) 2302 continue do 2303 k = 1 , kmax phi(k) = phih(k+1) + rgas*tbrr(k)*(qlgc4(k) -alog(q(k))) 2303 continue cc do 2399 k = 1 , kmax phi(k) = phi(k)/98000. 2399 continue cc write(6,2305)(tbrr(k), k = 1 , kmax) write(6,2305)( phi(k), k = 1 , kmax) 2305 format(2x,'phi:',(/,2x,10e11.5)) cc cc c c set up matrix a c call rodist call amatrix c do 871 id = 1 , 4 do 870 k = kbot + 1 , kmax do 880 j = 1 , jmx do 880 i = 1 , imx xxd(i,j) = uvp(i,j,k,id) - cpa(i,j,k,id) 880 continue call separ(xxd) do 890 j = 1 , jmx do 890 i = 1 , imx cpa(i,j,k,id) = cpa(i,j,k,id) + xxd(i,j) 890 continue 870 continue 871 continue cc do 893 j = 1 , jmx do 893 i = 1 , imx xxd(i,j) = pp(i,j) - pstar(i,j) 893 continue call separ(xxd) do 895 j = 1 , jmx do 895 i = 1 , imx pstar(i,j) = pstar(i,j) + xxd(i,j) 895 continue cc if(kbot .ne. 0) then do id = 1 , 4 do k = 1 , kbot do j = 1 , jmx do i = 1 , imx cpa(i,j,k,id) = uvp(i,j,k,id) enddo enddo enddo enddo endif npa = 36 ncs = 33 nps = 57 do 101 j = 1 , jmax write(npa)(((cpa(i,j,k,id), k = 1 , kmax) , i = 1 , imx) * , id = 1 , 4) cc write(nps)(pstar(i,j) , i = 1 , imx) cc cc c 101 continue cc cc fill in the frame array with data that will be in the window cc frame of the inner nests.... this is taken from smoothed data of cc the outer nest cc cc first for u,v,t,r cc if(nnest.ne.1)call frmin(1,frmn,frme) cc cc next for p-surface cc cc if(nnest.ne.1)call psfin(1,fpsn,fpse) cc cc cc cc cc************ this ends the nest 1 block ****************************** cc cc cc cc now set up for the interpolation of the derived basic field cc to the fine mesh resolutions cc cc if(nnest.ne.1)then cc nest1 = nnest - 1 do 110 ng = 1 , nest1 cc cc fill in the window frame data with values obtained from the outer cc nest....... cc cc cc cc cc start setting up for akima interpolation cc cc first do the pstar interpolation cc jb = jsg(ng) + 2 je = jng(ng) - 2 ib = iwg(ng) + 2 ie = ieg(ng) - 2 cc iold = ie - ib + 1 jold = je - jb + 1 cc cc ii = 0 do 3490 i = ib , ie ii = ii + 1 xd(ii) = del(i,3,ng) 3490 continue jj = 0 do 3491 j = jb , je jj = jj + 1 yd(jj) = tha(3,j,ng) 3491 continue ii = 0 do 990 j = jb , je do 990 i = ib , ie ii = ii + 1 zd(ii) = pstar(i,j) 990 continue ccc nxi = imag(ng+1) - 4 cc nyi = jmag(ng+1) - 4 jb = 3 ib = 3 je = jmag(ng+1) - 2 ie = imag(ng+1) - 2 izm = nxi do 933 i = 3 , ie xim(i-2) = del(i,3,ng+1) 933 continue do 935 j = 3 , je yim(j-2) = tha(3,j,ng+1) 935 continue cc cc do the interpolation to the finer mesh with the akima scheme cc call wghts(iold,jold,nxi,nyi,xd,yd,xim,yim,wti,wtj,idi,idj) call interp(zd,zim,idi,idj,wti,wtj,iold,jold,nxi,nyi) cc ii = 0 do 995 j = 1 , nyi do 995 i = 1 , nxi ii = ii + 1 pstar(i+2,j+2) = zim(ii) 995 continue cc cc 2070 continue cc cc first deal with the u component of the wind at each k level cc second deal with the v component of the wind at each k level cc thrid deal with the t component of the wind at each k level cc fourth deal with the r component of the wind at each k level cc c cc cc c c do 130 k = 1 , kmax cc do 1130 id = 1 , 4 ccc jb = jsg(ng) + 2 je = jng(ng) - 2 ib = iwg(ng) + 2 ie = ieg(ng) - 2 cc iold = ie - ib + 1 jold = je - jb + 1 cc ii = 0 do 3420 i = ib , ie ii = ii + 1 xd(ii) = del(i,3,ng) 3420 continue jj = 0 do 3421 j = jb , je jj = jj + 1 yd(jj) = tha(3,j,ng) 3421 continue cc cc ii = 0 do 120 j = jb , je do 120 i = ib , ie ii = ii + 1 zd(ii) = cpa(i,j,k,id) 120 continue cc cc nxi = imag(ng+1) - 4 nyi = jmag(ng+1) - 4 jb = 3 ib = 3 je = jmag(ng+1) - 2 ie = imag(ng+1) - 2 izm = nxi do 1933 i = 3 , ie xim(i-2) = del(i,3,ng+1) 1933 continue do 1935 j = 3 , je yim(j-2) = tha(3,j,ng+1) 1935 continue cc cc do the interpolation to the finer mesh with the akima scheme cc call wghts(iold,jold,nxi,nyi,xd,yd,xim,yim,wti,wtj,idi,idj) cc call interp(zd,zim,idi,idj,wti,wtj,iold,jold,nxi,nyi) cc cc place the interpolated values into the phase speed array ----cpa-- cc to be written out to tape cc ii = 0 do 150 j = 1 , nyi do 150 i = 1 , nxi ii = ii + 1 cpa(i+2,j+2,k,id) = zim(ii) 150 continue cc 1130 continue cc cc 130 continue cc nps = nps + 1 npa = npa + 1 ncs = ncs + 1 cc cc place the window frame array of the next inner grid into the speci cc arrays ---frmn---- and ---frme---. these will be placed on ---cpa- cc in the routine -------frmout cc do the same for the pstar field call frmout(ng+1,frmn,frme) cc call psfout(ng+1,fpsn,fpse) cc cc if(ng.ne.nest1)call frmin(ng+1,frmn,frme) cc if(ng.ne.nest1)call psfin(ng+1,fpsn,fpse) cc cc cc iml = nxi + 4 jml = nyi + 4 cc cc write the phase speeds out to dics cc cc do 170 j = 1 , jml write(npa)(((cpa(i,j,k,id), k = 1 , kmax) , i = 1 , imx) * ,id = 1 , 4) cc write(nps)(pstar(i,j) , i = 1 , imx) cc cc cc ccc cc 170 continue cc 110 continue cc cc***********end of the inner nest block******************* cc endif cc cc cc end of dealing with the basic flow field............. cc cc cc*********************************************************************** cc cc cc we will compute the symmetric and asym. vortex now cc this will interface with the code of r.ross 12/199o...... cc the symmetric vortex was determined by an earlier cc step of the initialization package........ cc where we ran a axis-symmetirc version of the hurricance cc model....................... cc routine ..init..... is the driver routine for this code cc cc********************* call init cc********************* cc cc the full bogus vortex has ben written out to disc in the routine cc nest which was called from init!!!!!!!!!!!!!!!!!! cc cc******************************************* cc cc cc write out the final history tape..........unit 77 cc cc read(55)nstep,nnest,iblock cc cc write(77)nstep,nnest,iblock do 3166 ip = 1 , ngr do 3166 jp = 1 , jblkmx diag(jp,ip) = 0.0 3166 continue ijaz = 0 do 3156 ir = 1 , nnest npa = 35 + ir nps = 56 + ir ncs = 32 + ir rewind npa rewind nps rewind ncs cc read(55)ngd,ngr,ntr,dt,js,jn,ie,iw,iimax,imax,jjmax,jmax, * nstflg,ickx,icky,ihx,ihy,dftx,dfty cc cc write(77)ngd,ngr,ntr,dt,js,jn,ie,iw,iimax,imax,jjmax,jmax, * nstflg,ickx,icky,ihx,ihy,dftx,dfty cc cc cc do 3244 j = 1 , jmax read(npa)(((cpa(i,j,k,id), k = 1 , kmax) , i = 1 , imx) * , id = 1 , 4) cc read(nps)(pstar(i,j), i = 1 , imx) cc 3244 continue cc nvx = 40 + ir nvp = 45 + ir j = 0 do 3177 jj = 1 , jjmax cc read(55)diggg read(55) ((curr(iq,jq),iq = 1 , iblkmx) , jq = 1 , ngr) njrow = ngr if(jj.le.2.or.jj.ge.jjmax-1)njrow = 1 do 3180 jt = 1 , njrow j = j + 1 do 3187 i = 1 , iblkmx filepc(i) = curr(i,jt) 3187 continue c cc cc cc finally add in the 3m vortex that is stored out on disc cc cc********************************************************* cc read(nvx)(((vortx(k,i,id), k = 1 , kmax), i = 1 , imx) * ,id = 1 , 4) c read(nvp)(psvr(i) , i = 1 , imx) cc cc********************************************************* cc do 3186 id = 1 , 3 do 3186 i = 1 , imax do 3186 k = 1 , kmax typ2c(k,i,id) = cpa(i,j,k,id) + vortx(k,i,id) 3186 continue ccc cc do 4146 id = 1 , 2 cc do 4146 i = 1 , imax cc do 4146 k = 1 , kmax cc typ2c(k,i,id) = vortx(k,i,id) cc4146 continue c do 4186 i = 1 , imax fctml = fctf if(typ1c(i,15).ge.0.0)fctml = amin1(fctf,.5) do 4186 k = 1 , kmax typ2c(k,i,4) = cpa(i,j,k,id) + vortx(k,i,4)*fctml 4186 continue do 3968 i = 1 , imax do 3968 k = 1 , kmax if(typ2c(k,i,4).lt.0.0)ijaz = ijaz + 1 if(typ2c(k,i,4).lt.0.0)typ2c(k,i,4) = 1.0e-20 3968 continue cc do 3183 i = 1 , imax typ1c(i,lgi) = pstar(i,j) + psvr(i) 3183 continue cc cc********************************************************* cc cc do 3186 i = 1 , imax cc typ1c(i,lgi) = pspr(i) cc 3186 continue cc do 2711 i = 1 , imax tbar = alog(1.0+tlap*typ1c(i,11)/typ2c(kmax,i,3)) epsi = tbar/(rgas*tlap) pst = typ1c(i,lgi)/exp(epsi) typ1c(i,lgi) = pst 2711 continue cc cc do 3188 i = 1 , iblkmx curr(i,jt) = filepc(i) 3188 continue ccc 3180 continue write(77) dggg write(77) ((curr(jk,ik),jk = 1 , iblkmx), ik = 1 , ngr) cc 3177 continue cc rewind npa cc 3156 continue ccc rewind 55 cc write(6,3903)ijaz 3903 format(2x,'ijaz',i6) cc stop cc end