program troughjet_rawwind_v4 ! This program generates objective trough lines and jet axes ! from streamfunction or wind field. Method used is detailed in ! Berry, Thorncroft and Hewson (2006), MWR. ! ---------------------------: CAUTION :------------------------------ ! ONLY set up to use cylindrical equidistant (CED) grids ! Mercator seems to work fine near the equator, but you've been warned. ! --------------------------------------------------------------------- ! ---- UPDATE LOG-------- !October 28 2005 - GB V2 !Allows tuning of lengthscale (# of gridpts) used in finite !differences - to cater for noisy psi field associated with !some hi-res datasets. two controls - one for wind calc (nopts) !one for derived quantities (dnopts) e.g. Vorticity !Minor changes to make code neater. you think this is a mess... !November 8 2005 - GB V2.1 !code modified to use 'raw' wind rather than non-divergent wind !Added flag to use non-global raw wind field !Added flag to use forecast times (FXXX) or 'real' times (F000) ! December 28 2005 - GB V3 ! UPGRADE-O-RAMA: ! Code transported to f90 standard. ! Automated reading of gempak grid navigation and determination of ! global or limited area datasets. ! option to use either streamfunction or 'raw wind' ! implementation of various subrountines to reduce prog size, ! particularly shear vorticity calculation. ! included rewritten f90 versions of gem_rd and gem_wrt as subrountines ! redundant arrays removed !------------------------------- ! Feb 06 - GB V3.03 ! Bug fix - GEMPAK output vertical coordinate was always pressure... ! GB V3.04 ! Anantha's smoothing algortihm added ! Oct 07 - GB V4.00 ! Attempt at incorperating tracking from Anantha's code. ! 4.01 - narrow down text output region !----------------------------------------------------------------------- ! NOTE: on DEAS suns compile with: ! f90 -e -lF77 ***.f90 $GEMLIB/gemlib.a $GEMLIB/cgemlib.a ! on aster: ! pgf90 -pgf77libs -g77libs -pgf90libs troughjetv4.00.f90 $GEMLIB/gemlib.a $GEMLIB/cgemlib.a !-------------------------------------------------------------------------- implicit none integer :: lev, nlon,nlat,yy,mm,dd,hh, date, dtemp integer ::m,n,s,t,dnopts,nopts,d,f real :: maga,dira,anglea,dotproa,ura,vra real :: magb,dirb,angleb,dotprob,urb,vrb real :: magc,dirc,anglec,dotproc,urc,vrc real :: magd,dird,angled,dotprod,urd,vrd ,x1 real :: refdir,you,vee, globetest,uncalc,amiss integer :: ik, sdate, edate, ntimes, level,g,h,iret,interval character (len=20) :: datec character (len=6) :: vvar,uvar,var,streamvar character (len=70) :: ifile, ifile2 integer :: imain, igem1, igem2, icoord, iret1, coord, globflag real :: deltax,deltay,radius,pi,vortmin,ucutoff,ucutoff2 real,dimension (:,:),allocatable :: u,v,vor,p real,dimension (:,:),allocatable :: shvor,cuvor,advcu,aej,fulladv real,dimension (:,:),allocatable :: ddxvor,ddyvor,aejcond,magwnd real,dimension (:,:),allocatable :: ddxcu,ddycu,ddxadv,ddyadv,mask real,dimension (:,:),allocatable :: advcuw,aejw,data1,data2,data3,data4 real,dimension (:,:),allocatable :: Tloce,Tlocw,Jloce,Jlocw integer :: vyes,uyes,styes,shyes,cuyes,fayes,mayes,fcflag real :: ll_lat,ll_lon,ur_lat,ur_lon,nglat_int,nglon_int integer :: fcstart,fcend,fcincr,fchr,streamflag,smvor,smcuv,gp1 logical :: wrtflg1=.false.,file_open integer :: igdfln, navsz, ihdrsz, maxgrd,ianlsz,tcounte,tcountw integer :: jcounte,jcountw,troughyes,req real,dimension(128) :: anlblk real,dimension(256) :: rnvblk ! for trough finder real, dimension(:,:), allocatable :: latArray, lonArray, vorArray real, dimension(:,:), allocatable :: xLat,xLon,xVor integer, dimension(:), allocatable :: axisCount integer, dimension(:,:), allocatable ::pos, look integer :: i,j,k,ii,jj,is,js,jj1,ii1,i1,j1 integer :: na,ic,nn,etrough,ejet real :: lat0,lat1,lon0,lon1 real :: latB,latT,lonL,LonR,minlength,crosslon,crosslat integer :: hB,hT,gL,gR,mingrid,meridflag,grange,hrange real :: searchrad real, dimension(:,:), allocatable :: latArray2, lonArray2, vorArray2 real, dimension(:,:), allocatable :: xLat2,xLon2,xVor2 integer, dimension(:), allocatable :: axisCount2 integer, dimension(:,:), allocatable :: pos2, look2 ! fixed parameters parameter (radius=6371000) !radius of earth in metres pi = 4.*atan(1.) amiss = -9999.0 !-------------------------------------------------------------------------- write(*,*)'---------------------------------------------------' write(*,*)' GARETHS TROUGH PROGRAM: V4.01 ' write(*,*)' ' write(*,*)'---------------------------------------------------' write(*,*)'ENTER INPUT FILENAME' read(*,'(a70)')ifile write(*,*)'ENTER OUTPUT FILENAME' read(*,'(a70)')ifile2 write(*,*) 'ENTER VERTICAL COORD 1=pres, 2=thta, 3=hght' read(*,*) coord write(*,*) 'ENTER LEVEL' read(*,*) level write(*,*) 'ENTER (1)REAL DATE (F000) or (2) FCAST (FXXX)' read(*,*) fcflag write(*,*) 'ENTER STARTING DATE: yymmddtt' read(*,*) sdate if (fcflag.eq.1)then write(*,*) 'ENTER ENDING DATE: yymmddtt' read(*,*) edate write(*,*)'ENTER INTERVAL (HRS)' read(*,*) interval else if(fcflag.eq.2)then write(*,*) 'ENTER FC START,END,INCREMENT (HRS)' read(*,*) fcstart,fcend,fcincr else write(*,*) 'CHECK YOUR DATES BOYO' stop endif write(*,*) 'USE STREAMFUCTION? 1=y, 0=n' read(*,*)streamflag if(streamflag.eq.1)then WRITE(*,*)'ENTER STREAMFUNCTION VARIABLE NAME (CAPS)' read(*,*) streamvar ! write(*,*) 'ENTER # of points used for wind finite diff' ! read(*,*) nopts nopts = 1 else if(streamflag.eq.0)then write(*,*) 'ENTER ZONAL (U) WIND VARIABLE NAME (CAPS)' read(*,*) uvar write(*,*) 'ENTER MERIDIONAL (V) WIND VARIABLE NAME (CAPS)' read(*,*) vvar else write(*,*) 'SOMETHING WRONG HERE OLD BEAN' stop endif ! write(*,*) 'ENTER # of points used for vort etc finite diff' ! read(*,*) dnopts dnopts = 1 write(*,*) '# of smoothing passings on vort' read(*,*) smvor write(*,*) 'ENTER MIN CURVATURE VORT FOR TROUGHS (*10-5)' read(*,*) vortmin write(*,*) 'OUTPUT U,V,VOR,SHVOR,CUVOR,FULLADV,mask? 1=Y 0=N' read(*,*) uyes,vyes,styes,shyes,cuyes,fayes,mayes write(*,*) 'ENTER WIND MAGN CUTOFF FOR AEJ AXIS' read(*,*) ucutoff write(*,*) 'ENTER ZONAL(U) WIND CUTOFF FOR AEWS' read(*,*) ucutoff2 write(*,*) 'OUTPUT TROUGH/JET LOC TO FILE? 1=y,0=n' read(*,*) troughyes if (troughyes.eq.1) then write(*,*) 'ENTER LAT BOUNDS FOR TROUGHFINDER' write(*,*)' FORMAT: lat1,lat2' read(*,*) latB,latT write(*,*) 'ENTER LON BOUNDS FOR TROUGHFINDER' write(*,*)' FORMAT: lonmin,lonmax' read(*,*) lonL,lonR write(*,*) 'ENTER MINIMUM TROUGH LENGTH (km)(Min jet length is 2x this)' read(*,*)minlength write(*,*) 'ENTER SEARCH RADIUS (km)' read(*,*)searchrad endif !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ! OPEN GEMPAK FILE, FIND NAVIGATION !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: call in_bdta ( iret ) if (iret /= 0) then print *, "Gempak initialization faliure...terminating" stop endif !!! open the file and get some information regarding the grids call gd_opnf ( ifile, wrtflg1, igdfln, navsz, rnvblk,& &ianlsz, anlblk, ihdrsz, maxgrd, iret ) file_open = .false. nlon = rnvblk(5) ; nlat = rnvblk(6) !assign the grid sizes ll_lat = rnvblk(7) ; ll_lon = rnvblk(8) ur_lat = rnvblk(9) ; ur_lon = rnvblk(10) if (ur_lon < 0.) then ur_lon = ur_lon+360. endif nglat_int=(ur_lat - ll_lat)/(nlat-1) nglon_int=(ur_lon - ll_lon)/(nlon-1) print*,'' print*,'------------------------------------------------------' print*,'' print*,'nx :',nlon,' ny :',nlat print*,'' print*,'LL lat =',ll_lat print*,'LL lon =',ll_lon print*,'UR lat =',ur_lat print*,'UR lon =',ur_lon print*,'' print*,'lon int:',nglon_int,'lat int:',nglat_int ! check for global grid.... globetest = (ur_lon + nglon_int) if (globetest /= (ll_lon+(360.))) then globflag = 0 print*,'ASSUMING NON-GLOBAL DATA' else globflag = 1 print*,'ASSUMING GLOBAL DATA' endif ! get ready to rotate data if required if((ll_lon.eq.0.).and.(globflag.eq.1))then meridflag = 1 endif ! -------- Outer date loop start ------------------------ date = sdate !set initial date as the starting date call de_date(date,yy,mm,dd,hh) ik = 0 if(fcflag.eq.1)then ntimes = 365*4*60 else if(fcflag.eq.2)then ntimes = ((fcend-fcstart)/fcincr)+1 print*,'Number of times to do =',ntimes endif do 555 imain = 1, ntimes if(fcflag.eq.1)then if(date.gt.edate) goto 600 endif dtemp = int(float(date)/1.e8) dtemp = (date - int( float(dtemp)*1.e8) )/100 hh = date - (date/100 ) *100 ! code the gempak date if(fcflag.eq.1)then write(datec,'(i6.6,a1,i2.2,a6)') dtemp,'/',hh,'00F000' endif if(fcflag.eq.2)then fchr = fcstart+(ik*fcincr) write(datec,'(i6.6,a1,i2.2,a3,i3.3)') dtemp,'/',hh,'00F',fchr endif print*,'GEMPAK DATE IS:',datec ! Allocate all the arrays here allocate (p(nlon,nlat)) allocate (u(nlon,nlat)) allocate (v(nlon,nlat)) allocate (vor(nlon,nlat)) allocate (shvor(nlon,nlat)) allocate (cuvor(nlon,nlat)) allocate (advcu(nlon,nlat)) allocate (aej(nlon,nlat)) allocate (advcuw(nlon,nlat)) allocate (aejw(nlon,nlat)) allocate (fulladv(nlon,nlat)) allocate (ddxvor(nlon,nlat)) allocate (ddyvor(nlon,nlat)) allocate (aejcond(nlon,nlat)) allocate (magwnd(nlon,nlat)) allocate (ddxcu(nlon,nlat)) allocate (ddycu(nlon,nlat)) allocate (ddxadv(nlon,nlat)) allocate (data1(nlon,nlat)) allocate (data2(nlon,nlat)) allocate (data3(nlon,nlat)) allocate (data4(nlon,nlat)) allocate (ddyadv(nlon,nlat)) allocate (mask(nlon,nlat)) allocate (Tloce((nlon*nlat),3)) allocate (Tlocw((nlon*nlat),3)) allocate (Jloce((nlon*nlat),3)) allocate (Jlocw((nlon*nlat),3)) !............................................................. !-------reading part !-------read the gempak file for the given date/variable ! if (streamflag.eq.1) then var = streamvar igem1 = 0 igem2 = 2 icoord = coord !1 pres, 2=thta, 3=hght lev = level call gem_rd(p,nlon,nlat,lev,datec,icoord,ifile,var,igem1,iret1) if(iret1.eq.0) then write(*,'(a24,a20)') 'Read streamfunction for:', datec else call gem_err(datec,var,ifile,lev,igem1,iret1) endif ! ---------- CALC U,V(psi) FROM Streamfunction ------------------ do g=1,nlon do h=1,nlat !calculate deltax, deltay call calc_dxdy(nlon,nlat,g,h,globflag,ll_lat,nglon_int,nglat_& &int,deltax,deltay) ! calculate v component if ((h.le.nopts).or.(h.gt.(nlat-nopts))) then v(g,h)=amiss else m=g+nopts n=g-nopts ! account for split in dataset if (m.gt.nlon)then m=m-nlon endif if (n.le.0) then n=nlon+n endif ! simple centred difference v(g,h)=(((p(m,h))-(p(n,h)))/(2*nopts*deltax)) endif ! calc u component ! Set upper parts to missing if ((h.le.nopts).or.(h.gt.(nlat-nopts))) then u(g,h)=amiss else ! note since dif is in merid dir, data wrap not needed m=h+nopts n=h-nopts ! centred difference: u(g,h)=(((p(g,m))-(p(g,n)))/(-2*nopts*deltay)) endif enddo enddo ! Write data to file if (uyes.eq.1) then var='UWND' igem1=0 igem2=2 icoord=1 lev=level CALL GEMWRT(u,nlon,nlat,lev,datec,icoord,ifile2,var,igem1) endif if (vyes.eq.1) then var='VWND' igem1=0 igem2=2 icoord=1 lev=level CALL GEMWRT(v,nlon,nlat,lev,datec,icoord,ifile2,var,igem1) endif print*,'WRITTEN U,V',datec else if(streamflag.eq.0) then !-------- UWND -------- var = uvar igem1 = 0 igem2 = 2 icoord = coord !1 pres, 2=thta, 3=hght lev = level call gem_rd(u,nlon,nlat,lev,datec,icoord,ifile,var,igem1,iret1) if(iret1.eq.0) then write(*,'(a13,a20)') 'Read U ', datec else call gem_err(datec,var,ifile,lev,igem1,iret1) endif !-------- VWND ------- var = vvar igem1 = 0 igem2 = 2 icoord = coord !1 pres, 2=thta, 3=hght lev = level call gem_rd(v,nlon,nlat,lev,datec,icoord,ifile,var,igem1,iret1) if(iret1.eq.0) then write(*,'(a13,a20)') 'Read V ', datec else write(*,*) 'ERROR OCCURED AT THE FOLLOWING DATE' write(*,'(a20)') datec write(*,'(a6)') var write(*,'(a45)') ifile write(*,*) lev, igem1,iret1 write(*,*) 'EITHER FILE NOT PRESENT OR DATA FOR' write(*,*) 'SPECIFIED DATE/TIME NOT PRESENT' write(*,*) 'CHECK INPUT FILE FOR DATA AVAILABILTY' write(*,*) 'TERMINATING PROGRAM' stop endif ! ATTEMPT SMOOTHING call smoothDat ( u,nlon,nlat,amiss,smvor ) call smoothDat ( v,nlon,nlat,amiss,smvor ) !Write data to file if (uyes.eq.1) then var='UWND' igem1=0 igem2=2 icoord=coord lev=level CALL GEMWRT(u,nlon,nlat,lev,datec,icoord,ifile2,var,igem1) endif if (vyes.eq.1) then var='VWND' igem1=0 igem2=2 icoord=coord lev=level CALL GEMWRT(v,nlon,nlat,lev,datec,icoord,ifile2,var,igem1) endif print*,'WRITTEN U,V ',datec endif ik = ik + 1 ! Set number of uncalculatable grid points - can be avoided if I get round ! to implementing a more exotic finite differencing scheme if(streamflag.eq.1) then uncalc = int(dnopts+nopts) else if(streamflag.eq.0) then uncalc = int(dnopts+1) endif ! --------------- CALCULATE VORTICITY --------------- ! note that having to lose another point near poles ! set the non-calculatable values to -9999.9 do g=1,nlon do h=1,nlat if ((h.le.uncalc).or.(h.gt.(nlat-uncalc)))then vor(g,h)=amiss else if ((globflag.eq.0).and.((g.le.uncalc).or.& &(g.gt.(nlon-uncalc))))then vor(g,h)=amiss else !calculate deltax, deltay call calc_dxdy(nlon,nlat,g,h,globflag,ll_lat,nglon_int,nglat_& &int,deltax,deltay) ! for points at the edges: m=g+dnopts n=g-dnopts s=h+dnopts t=h-dnopts ! account for split in dataset if (m.gt.nlon)then m=m-nlon endif if (n.le.0) then n=nlon+n endif ! vortcity calculation vor(g,h)=((v(m,h)-v(n,h))/((2*dnopts)*deltax))& &-((u(g,s)-u(g,t))/((2*dnopts)*deltay)) ! multiply vortcity by 10^5 to make it easier to handle: vor(g,h)=vor(g,h)*100000. endif enddo enddo ! Write vorticity to file if (styes.eq.1) then var='VORT' igem1=0 igem2=2 icoord=coord lev=level CALL GEMWRT(vor,nlon,nlat,lev,datec,icoord,ifile2,var,igem1) endif print*,'VORT DONE' ! --------------- Calculate Shear vorticity ------------- !set the non-calculatable values to -9999.9 do g=1,nlon do h=1,nlat if ((h.le.uncalc).or.(h.gt.(nlat-uncalc)))then shvor(g,h)=amiss else if ((globflag.eq.0).and.((g.le.uncalc).or.& &(g.gt.(nlon-uncalc))))then shvor(g,h)=amiss else ! calculate deltax, deltay call calc_dxdy(nlon,nlat,g,h,globflag,ll_lat,nglon_int,nglat_& &int,deltax,deltay) ! for point n of ref ! c ! b a ! d !General points: m=g+dnopts n=g-dnopts s=h+dnopts t=h-dnopts ! account for split in dataset if (m.gt.nlon)then m=m-nlon endif if (n.le.0) then n=nlon+n endif ! set reference vector at central point ! Subroutine determine which qudrant vectors is in ! and calculate angle in radians. call calc_refdir(u(g,h),v(g,h),refdir) ! point a you=u(m,h) vee=v(m,h) call pt_dir(you,vee,refdir,dotproa) ! recast dot product in grid relative coordinates ura=dotproa*sin(refdir) vra=dotproa*cos(refdir) ! point b you=u(n,h) vee=v(n,h) call pt_dir(you,vee,refdir,dotproa) ! recast dot product in grid relative coordinates urb=dotproa*sin(refdir) vrb=dotproa*cos(refdir) ! point c you=u(g,s) vee=v(g,s) call pt_dir(you,vee,refdir,dotproa) ! recast dot product in grid relative coordinates urc=dotproa*sin(refdir) vrc=dotproa*cos(refdir) ! point d you=u(g,t) vee=v(g,t) call pt_dir(you,vee,refdir,dotproa) ! recast dot product in grid relative coordinates urd=dotproa*sin(refdir) vrd=dotproa*cos(refdir) !now Calculate the Shear vorticity shvor(g,h)=(((vra-vrb)/((2*dnopts)*deltax))& &-((urc-urd)/((2*dnopts)*deltay))) ! scale by 10^5 to match full vort shvor(g,h)=shvor(g,h)*100000 endif enddo enddo ! Write data to file if (shyes.eq.1) then var='SHVOR' igem1=0 igem2=2 icoord=coord lev=level CALL GEMWRT(shvor,nlon,nlat,lev,datec,icoord,ifile2,var,igem1) endif !---------- CURVATURE VORTICITY ------------------ !Calculate curvature vort as residual of full and shear vort !set the non-calculatable values to -9999.9 do g=1,nlon do h=1,nlat if ((h.le.uncalc).or.(h.gt.(nlat-uncalc)))then cuvor(g,h)=amiss else if ((globflag.eq.0).and.((g.le.uncalc).or.& &(g.gt.(nlon-uncalc))))then cuvor(g,h)=amiss else ! all the other points: cuvor(g,h)=vor(g,h)-shvor(g,h) endif enddo enddo if(cuyes.eq.1)then var='CUVOR' igem1=0 igem2=2 icoord=coord lev=level CALL GEMWRT(cuvor,nlon,nlat,lev,datec,icoord,ifile2,var,igem1) endif print*,'CURV VORT DONE' ! ------------ Advection of curvature vort -------------- ! set the non-calculatable values to -9999.9 do g=1,nlon do h=1,nlat if ((h.le.(uncalc+1)).or.(h.gt.(nlat-(uncalc+1))))then advcu(g,h)=amiss advcuw(g,h)=amiss else if ((globflag.eq.0).and.((g.le.(uncalc+1)).or.& &(g.gt.(nlon-(uncalc+1))))) then advcu(g,h)=amiss advcuw(g,h)=amiss else ! calculate deltax, deltay call calc_dxdy(nlon,nlat,g,h,globflag,ll_lat,nglon_int,nglat_& &int,deltax,deltay) ! deriavtives m=g+dnopts n=g-dnopts s=h+dnopts t=h-dnopts ! account for split in dataset if (m.gt.nlon)then m=m-nlon endif if (n.le.0) then n=nlon+n endif ! simple centred difference ddxcu(g,h)=(cuvor(m,h)-cuvor(n,h))/(2*dnopts*deltax) ddycu(g,h)=(cuvor(g,s)-cuvor(g,t))/(2*dnopts*deltay) advcu(g,h)=-(u(g,h)*ddxcu(g,h))-(v(g,h)*ddycu(g,h)) advcuw(g,h)= advcu(g,h) ! spit out unmasked advcu if you want to do graphically fulladv(g,h)=advcu(g,h) endif enddo enddo if(fayes.eq.1)then var='FADV' igem1=0 igem2=2 icoord=1 lev=level CALL GEMWRT(fulladv,nlon,nlat,lev,datec,icoord,ifile2,var,igem1) endif print*,'FULLADV DONE' ! MASKING VARIABLE - REMOVE TROUGHS WHERE V(str).(del)advcu not >0 ! set the non-calculatable values to -9999.9 do g=1,nlon do h=1,nlat if ((h.le.(uncalc+2)).or.(h.gt.(nlat-(uncalc+2))))then mask(g,h)=amiss else if ((globflag.eq.0).and.((g.le.(uncalc+2)).or.& &(g.gt.(nlon-(uncalc+2))))) then mask(g,h)=amiss else ! calculate deltax, deltay call calc_dxdy(nlon,nlat,g,h,globflag,ll_lat,nglon_int,nglat_& &int,deltax,deltay) ! General points: m=g+dnopts n=g-dnopts s=h+dnopts t=h-dnopts ! account for split in dataset if (m.gt.nlon)then m=m-nlon endif if (n.le.0) then n=nlon+n endif ! calc derivatives: ddxadv(g,h)=(advcu(m,h)-advcu(n,h))/(2*dnopts*deltax) ddyadv(g,h)=(advcu(g,s)-advcu(g,t))/(2*dnopts*deltay) ! calc V.(del(advcu)) mask(g,h)=(u(g,h)*ddxadv(g,h))+(v(g,h)*ddyadv(g,h)) ! mask(g,h)=mask(g,h)*100 endif enddo enddo if(mayes.eq.1)then var='MASK' igem1=0 igem2=2 icoord=1 lev=level CALL GEMWRT(mask,nlon,nlat,lev,datec,icoord,ifile2,var,igem1) endif print*,'MASK DONE' ! -------- TROUGH AXES (AT LAST) ---------- ! dump non +ve points if(globflag.eq.1)then do g=1,nlon do h=(uncalc+2),(nlat-(uncalc+2)) if((cuvor(g,h).lt.vortmin).or.(u(g,h).gt.ucutoff2).& &or.(mask(g,h).lt.0.))then advcu(g,h)=amiss endif if((cuvor(g,h).lt.vortmin).or.(u(g,h).lt.ucutoff2).& &or.(mask(g,h).lt.0.))then advcuw(g,h)=amiss endif enddo enddo else if (globflag.eq.0) then do g=(uncalc+2),(nlon-(uncalc+2)) do h=(uncalc+2),(nlat-(uncalc+2)) if((cuvor(g,h).lt.vortmin).or.(u(g,h).gt.ucutoff2).& &or.(mask(g,h).lt.0.))then advcu(g,h)=amiss endif if((cuvor(g,h).lt.vortmin).or.(u(g,h).lt.ucutoff2).& &or.(mask(g,h).lt.0.))then advcuw(g,h)=amiss endif enddo enddo endif var='TROFE' igem1=0 igem2=2 icoord=coord lev=level CALL GEMWRT(advcu,nlon,nlat,lev,datec,icoord,ifile2,var,igem1) var='TROFW' igem1=0 igem2=2 icoord=coord lev=level CALL GEMWRT(advcuw,nlon,nlat,lev,datec,icoord,ifile2,var,igem1) print*,'TROUGHS DONE' ! ------------ JET AXIS -------------- ! set where shear vort = 0, in regions where ! V.(del(shearvort) x K) > 0 ! set the non-calculatable values to -9999.9 do g=1,nlon do h=1,nlat if ((h.le.(uncalc+2)).or.(h.gt.(nlat-(uncalc+2))))then aej(g,h)=amiss aejw(g,h)=amiss else if ((globflag.eq.0).and.((g.le.(uncalc+2)).or.& &(g.gt.(nlon-(uncalc+2))))) then aej(g,h)=amiss aejw(g,h)=amiss else ! calc x, y derivs of shear vort for all the other points: ! calculate deltax, deltay call calc_dxdy(nlon,nlat,g,h,globflag,ll_lat,nglon_int,nglat_& &int,deltax,deltay) ! General points: m=g+dnopts n=g-dnopts s=h+dnopts t=h-dnopts ! account for split in dataset if (m.gt.nlon)then m=m-nlon endif if (n.le.0) then n=nlon+n endif ! calc derivatives: ddxvor(g,h)=((shvor(m,h)-shvor(n,h))/(2*dnopts*deltax)) ddyvor(g,h)=((shvor(g,s)-shvor(g,t))/(2*dnopts*deltay)) magwnd(g,h)=sqrt(((u(g,h)**2)+(v(g,h)**2))) ! calc V.(del(shvort) x K) aejcond(g,h)=-((v(g,h)*ddxvor(g,h))-(u(g,h)*ddyvor(g,h))) if((aejcond(g,h).ge.0.).and.(magwnd(g,h).ge.ucutoff).and.& &(u(g,h).lt.0.)) then aej(g,h)=shvor(g,h) else aej(g,h)=amiss endif if((aejcond(g,h).ge.0.).and.(magwnd(g,h).ge.ucutoff).and.& &(u(g,h).gt.0.)) then aejw(g,h)=shvor(g,h) else aejw(g,h)=amiss endif endif enddo enddo var='JETE' igem1=0 igem2=2 icoord=coord lev=level CALL GEMWRT(aej,nlon,nlat,lev,datec,icoord,ifile2,var,igem1) var='JETW' igem1=0 igem2=2 icoord=coord lev=level CALL GEMWRT(aejw,nlon,nlat,lev,datec,icoord,ifile2,var,igem1) print*,'JETS DONE' if(troughyes.eq.1) then call trough_e(ll_lat,ll_lon,nglon_int,nglat_int,advcu,nlon,nlat,& &latB,latT,lonL,lonR,cuvor,minlength,datec,searchrad) call jet_e(ll_lat,ll_lon,nglon_int,nglat_int,aej,nlon,nlat,& &latB,latT,lonL,lonR,magwnd,minlength,datec,searchrad) 999 continue else print*,'NO TROUGHS FOR YOU!' endif ! Clean up arrays if(fcflag.eq.1)then call compute_date(yy,mm,dd,hh,date,interval) endif 555 continue deallocate (p,u,v,vor,shvor,cuvor,advcu,aej,advcuw,aejw) deallocate (fulladv,ddxvor,ddyvor,aejcond,magwnd) deallocate (ddxcu,ddycu,ddxadv,ddyadv,mask) deallocate (data1,data2,data3,data4) !--------------------- 600 continue write (*,*) '----------DONE------------' write (*,*) '----------DONE------------' write (*,*) write (*,*) 'Total times processed = ', ik write (*,*) write (*,*) '----------DONE------------' write (*,*) '----------DONE------------' stop end !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: subroutine de_date(date,yy,mm,dd,hh) ! Subroutine to split date block into year,month,day etc. ! written by Anantha. integer :: yy,mm,dd,hh,date yy = int(float(date)/1.e6) mm = int(float(date-yy*1000000)/1.e4) dd = int(float(date - yy*1000000 - mm*10000)/1.e2) hh = date - yy*1000000 - mm*10000 - dd*100 return end subroutine de_date !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: subroutine compute_date(yy,mm,dd,hh,date,interval) ! Subroutine to generate next gempak date ! written by Anantha. integer :: yy,mm,dd,hh,date,interval integer :: md1(12),md2(12) data md1 /31,28,31,30,31,30,31,31,30,31,30,31/ !non leap year data md2 /31,29,31,30,31,30,31,31,30,31,30,31/ !leap year ileap = 0 if (mod(float(yy),4.).eq.0) ileap = 1 if (yy.eq.2000) ileap = 0 hh = hh + interval if(hh.eq.24) then hh = 0 dd = dd + 1 if(ileap.eq.1) then !leap year if(dd.gt.md2(mm)) then dd = 1 mm = mm + 1 if (mm.eq.13) then mm = 1 yy = yy+1 endif endif else if(dd.gt.md1(mm)) then !non leap year dd = 1 mm = mm + 1 if (mm.eq.13) then mm = 1 yy = yy+1 endif endif endif endif date=(yy*1000000)+(mm*10000)+(dd*100)+hh return end subroutine compute_date !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::; subroutine GEMWRT(XDUM,IX,IY,IPRES,ITIME,IVCORD,FILE,PARM,IRET) ! Subroutine to write array to gempak file as a gempak variable ! original by G.Lackman modified by Anantha, converted to f90 by GB character (len=6) :: parm character (len=20) :: gdattm(2) character (len=70) ::file character (len=15) ::ITIME dimension :: XDUM(IX,IY) real :: rnvblk(256),anlblk(128) integer :: level(2),ighdr(2) data ipktyp /1/, nbits/16/ logical :: wrtflg=.true. ,rewrit=.true. ! print*,'TIME =',ITIME igx = ix igy = iy level(1)=IPRES level(2)=-1 gdattm(1) = ITIME gdattm(2) = ' ' call in_bdta(iret) call gd_opnf(file,wrtflg,igdfln,navsz,rnvblk,ianlsz,anlblk,& &ihdrsz,maxgrd,iret) if (iret .ne. 0) print*, 'output open error: ',iret call gd_wpgd(igdfln,XDUM,igx,igy,ighdr,gdattm,level,ivcord,& &parm,wrtflg,ipktyp,nbits,iret) if (iret .ne. 0) print*, 'write error: ',iret call gd_clos(igdfln,iret) RETURN !here i am returning to main program end subroutine GEMWRT !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: subroutine GEM_RD(XDUM,IX,IY,IPRES,ITIME,IVCORD,FILE,PARM,IGDFLN,IRET) ! subroutine to read gempak variable into an array. ! original by G.Lackman modified by Anantha, converted to f90 by GB character (len=6) :: parm character (len=20) :: gdattm(2) character (len=70) ::file character (len=15) ::ITIME dimension :: XDUM(IX,IY) real :: rnvblk(256),anlblk(128) integer :: level(2),ighdr(2) data ipktyp /1/, nbits/16/ logical :: wrtflg=.true. ,rewrit=.true. igx = ix igy = iy level(1)=IPRES level(2)=-1 gdattm(1) = ITIME gdattm(2) = ' ' if (igdfln.eq.0) then call in_bdta(iret) call gd_opnf(file,wrtflg,igdfln,navsz,rnvblk,ianlsz,anlblk,& &ihdrsz,maxgrd,iret) if (iret .ne. 0) print*, 'open error: ',iret endif call gd_rdat (igdfln, gdattm, level, ivcord, parm,& & XDUM, igx, igy, ighdr, iret) if (iret .ne. 0) print*, 'read error: ', iret if (iret .ne. 0) write(99,*)'Read error on this date ',gdattm,parm,iret RETURN end subroutine gem_rd !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: subroutine calc_refdir(u,v,refdir) ! Calculates the direction of wind vector ! used in computation of shear vort as the central 'reference vector'. real :: refdir,u,v,pi,u2,v2 ! print*,'U = ',u,' V = ',v u2=u**2 v2=v**2 pi = 4.*atan(1.) ! case of non-zero components if ((u.ge.(epsilon(1.))).and.(v.ge.(epsilon(1.))))then refdir=atan(u/v) elseif ((u.ge.(epsilon(1.))).and.(v.le.(-1*epsilon(1.))))then refdir=pi+(atan(u/v)) elseif ((u.le.(-1*epsilon(1.))).and.(v.le.(-1*epsilon(1.))))then refdir=(pi)+(atan(u/v)) elseif ((u.le.(-1*epsilon(1.))).and.(v.ge.(epsilon(1.))))then refdir=(2*pi)+(atan(u/v)) endif ! in case one or both components equal zero if((u2.lt.(epsilon(1.))).and.(v.ge.(epsilon(1.))))then refdir=0. elseif((u2.lt.(epsilon(1.))).and.(v.le.(-1*epsilon(1.))))then refdir=pi endif if((v2.lt.(epsilon(1.))).and.(u.ge.(epsilon(1.))))then refdir=(pi/2.) elseif((v2.lt.(epsilon(1.))).and.(u.le.(-1*epsilon(1.))))then refdir=((3*pi)/2.) endif ! case both components are zero: ! Assume ref vector is zero radians - may require workaound if((v2.lt.(epsilon(1.)).and.(u2.lt.(epsilon(1.)))))then refdir=0. endif RETURN end subroutine calc_refdir !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: subroutine pt_dir(you,vee,refdir,dotproa) ! subrountine to calculate the dot product of wind vector with ! reference vector (see other subroutine). ! used in calculation of shear vorticity. real :: you,vee,pi,maga,dira,you2,vee2,dotpro pi = 4.*atan(1.) maga=sqrt((you**2)+(vee**2)) you2 = you**2 vee2 = vee**2 ! for non-zero components if ((you.ge.(epsilon(1.))).and.(vee.ge.(epsilon(1.))))then dira=atan(you/vee) elseif ((you.ge.(epsilon(1.))).and.(vee.le.(-1*epsilon(1.))))then dira=pi+(atan(you/vee)) elseif ((you.le.(-1*epsilon(1.))).and.(vee.le.(-1*epsilon(1.))))then dira=pi+(atan(you/vee)) elseif ((you.le.(-1*epsilon(1.))).and.(vee.ge.(epsilon(1.))))then dira=(2*pi)+(atan(you/vee)) endif ! in case one component is equal to zero if((you2.lt.(epsilon(1.))).and.(vee.ge.(epsilon(1.))))then dira=0. elseif((you2.lt.(epsilon(1.))).and.(vee.le.(-1*epsilon(1.))))then dira=pi endif if((vee2.lt.(epsilon(1.))).and.(you.ge.(epsilon(1.))))then dira=(pi/2.) elseif((vee2.lt.(epsilon(1.))).and.(you.le.(-1*epsilon(1.))))then dira=((3*pi)/2.) endif ! Compute the dot product dotproa=maga*(cos(refdir-dira)) ! case both components are zero: if((vee2.eq.0.).and.(you2.eq.0.))then dotproa=0. endif return end subroutine pt_dir !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: subroutine calc_dxdy(nlon,nlat,g,h,globflag,ll_lat,nglon_int,nglat_& &int,deltax,deltay) ! Calculates grid spacing on a CED grid ! Gets grid dimension data from gempak file via main program integer :: nlon,nlat,g,h,globflag real :: pi,deltax,deltay,xgridpt,radius real :: nglat_int,nglon_int,latofptdeg,latofptrad,ll_lat radius=6371000. !radius of earth in metres pi = 4.*atan(1.) if (globflag.eq.1)then xgridpt = (float(h-1)/float(nlat-1))*pi deltax=(((sin(xgridpt))*radius)*(2*pi))/float(nlon) deltay=(pi*radius)/float(nlat-1) else if (globflag.eq.0) then latofptdeg = (ll_lat)+((h-1)*nglat_int) latofptrad = ((latofptdeg+90)/180.)*pi deltax=((((pi/180)*radius)*(sin(latofptrad))))*nglon_int deltay=((pi/180)*radius)*nglat_int endif return end subroutine calc_dxdy ! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: subroutine gem_err(datec,var,ifile,lev,igem1,iret1) character (len=20) :: datec,var character (len=70) :: ifile integer :: lev,igem1,iret1 ! just outputs a message when gempak data cannot be found write(*,*) '' write(*,*) '' write(*,*) ' E R R O R ' write(*,*) '' write(*,*) 'OCCURED AT DATE:' write(*,'(a20)') datec write(*,*) 'LOOKING FOR VARIABLE:' write(*,'(a6)') var write(*,*) 'IN FILE:' write(*,'(a45)') ifile write(*,*) 'LEV,IGEM1,IRET:' write(*,*) lev, igem1,iret1 write(*,*) '*******************************************' write(*,*) 'EITHER FILE NOT PRESENT OR DATA FOR' write(*,*) 'SPECIFIED DATE/TIME NOT PRESENT' write(*,*) 'CHECK INPUT FILE FOR DATA AVAILABILTY' write(*,*) 'TERMINATING PROGRAM' write(*,*) '*******************************************' write(*,*) '' write(*,*) '' write(*,*) '' stop end subroutine gem_err !------------------------------------------------------- subroutine smoothDat ( data,nx,ny,amiss,nsmooth ) ! Anantha's smoothing routine implicit none real, dimension (nx,ny) :: data real, dimension (:,:), allocatable :: Ldat,Rdat,Tdat,Bdat,Sdat,ic,& iarray real :: amiss integer :: nx,ny,nsmooth,i,j,ii allocate ( Ldat(nx,ny) ) ; allocate ( Rdat(nx,ny) ) allocate ( Tdat(nx,ny) ) ; allocate ( Bdat(nx,ny) ) allocate ( Sdat(nx,ny) ) ; allocate ( ic(nx,ny) ) allocate ( iarray(nx,ny)) do ii = 1,nsmooth ! print *, "Now smoothing", nx, ny, amiss Tdat=amiss ; Bdat = amiss ; Ldat=amiss ; Rdat = amiss Ldat(1:nx-1,:) = data(2:nx,:) Rdat(2:nx,:) = data(1:nx-1,:) Tdat(:,2:ny) = data(:,1:ny-1) Bdat(:,1:ny-1) = data(:,2:ny) ic = 0 iarray = 1 where ( Ldat == amiss ) iarray = 0 end where Sdat = Ldat*iarray ic = ic + iarray iarray = 1 where ( Rdat == amiss ) iarray = 0 end where Sdat = Sdat + Rdat*iarray ic = ic + iarray iarray = 1 where ( Tdat == amiss ) iarray = 0 end where Sdat = Sdat + Tdat*iarray ic = ic + iarray iarray = 1 where ( Bdat == amiss ) iarray = 0 end where Sdat = Sdat + Bdat*iarray ic = ic + iarray iarray = 1 where ( data == amiss ) iarray = 0 end where Sdat = Sdat + 2.*data*iarray ic = ic + 2.*iarray data = amiss where ( abs(ic) > epsilon(1.) ) data = sdat/ic end where !!$ !!$ enddo deallocate (Ldat,Rdat,Tdat,Bdat,Sdat,ic,iarray) return end subroutine smoothDat ! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: !-------------------------------------------------------- subroutine trough_e(ll_lat,ll_lon,nglon_int,nglat_int,axisdata,nlon,nlat,& &latB,latT,lonL,lonR,data,minlength,datec,searchrad) !!$! TROUGH TRACKING !!$! adapted from Anantha's code. sub.TroughAxis.f90 !!$! Looks along latitude lines for where advcu = 0 real :: ll_lat,ll_lon,nglon_int,nglat_int,latB,latT,lonL,lonR,x1 integer :: nlon,nlat,tcounte,hB,hT,gL,gR,grange,hrange,req integer :: g,h real, dimension (nlon,nlat) :: axisdata,data real, dimension (nlon*nlat,3) :: Tloce character (len=20) :: datec real, dimension(:,:), allocatable :: latArray, lonArray, vorArray real, dimension(:,:), allocatable :: xLat,xLon,xVor integer, dimension(:), allocatable :: axisCount,loccheck real, dimension(:,:), allocatable :: pos, look integer :: i,j,k,ii,jj,is,js,jj1,ii1,i1,j1,ip1 integer :: na,ic,nn,etrough,ejet real :: lat0,lat1,lon0,lon1 real :: minlength,crosslon,crosslat,searchrad integer :: mingrid,meridflag,searchradi tcounte=0 hB = ((nint((latB-ll_lat)/nglat_int))+1) hT = ((nint((latT-ll_lat)/nglat_int))+1) hrange = hT-hB do h=hB,hT do g=1,nlon if(axisdata(g,h).ne.amiss)then if (abs(axisdata(g,h)) <= epsilon(1.)) then ! zero line exactly at grid pt tcounte = tcounte+1 Tloce(tcounte,1)=ll_lon + (g-1)*nglon_int !lon Tloce(tcounte,2)=ll_lat + (h-1)*nglat_int !lat Tloce(tcounte,3)=data(g,h) else ! interpolate in lon to find zero line ip1 = g+1 ; if ( ip1 > nlon ) ip1 = ip1 - nlon if (axisdata(g,h) > 0. .and. axisdata(ip1,h) < 0.) then ! for easterlies x1 = ll_lon + (g-1)*nglon_int tcounte = tcounte+1 Tloce(tcounte,1)=(x1+ (nglon_int*(axisdata(g,h)/(axisdata(g,h)-axisdata(g+1,h))))) !lon Tloce(tcounte,2)=ll_lat + (h-1)*nglat_int !lat Tloce(tcounte,3)=max(data(g,h),data(g+1,h)) ! print*,tcounte,Tloce(tcounte,1),Tloce(tcounte,2),Tloce(tcounte,3) endif endif endif enddo enddo print*,'NUMBER OF TROUGH POINTS:' print*,'' print*,tcounte,'EASTERLY' ! now the axis finder routine - easterly open (unit=12, file='trough_east.out') allocate ( xlat(Tcounte,1000) ) allocate ( xlon(Tcounte,1000) ) allocate ( xvor(Tcounte,1000) ) allocate ( axisCount(Tcounte) ) allocate ( loccheck(Tcounte) ) allocate (pos(nlon,nlat) ) allocate (look(nlon,nlat) ) allocate (latArray(nlon,nlat)) allocate (lonArray(nlon,nlat)) allocate (vorArray(nlon,nlat)) searchradi=nint(searchrad/((((nglon_int+nglat_int)*0.5)*111))) do g=1,nlon do h=1,nlat pos(g,h) = 0 look(g,h) = 0 enddo enddo do k = 1, Tcounte h = 1 + nint( (Tloce(k,2)-ll_lat) / nglat_int ) g = 1 + nint( (Tloce(k,1)-ll_lon) / nglon_int ) pos(g,h) = 1 latArray(g,h) = Tloce(k,2) lonArray(g,h) = Tloce(k,1) vorArray(g,h) = Tloce(k,3) enddo do g=1,nlon do h=1,nlat if(pos(g,h).eq.0.)then look(g,h) = 1 endif enddo enddo ic = 0 ; na = 0 !counter to track the individual trough axes ! Concentrate within the domain of interest do g=1,nlon do h=hB,hT nn = 0 !counter to track grids under a given trough if(look(g,h) == 0 .and. pos(g,h) == 1) then !!$!---we have found a fresh trough axis! na = na + 1 !troughs nn = nn + 1 !points ic=ic+1 xlat(na,nn) = latArray(g,h) xlon(na,nn) = lonArray(g,h) xvor(na,nn) = vorArray(g,h) look(g,h) = 1 ii = g ; jj = h !temporarily store the location of the axis begin pt. 330 continue !the previous point now is the seed point to look for the rest of the axis is = ii js = jj axisCount(na) = nn !we look in the vicinity of the seed point for the remainder of the trough ! do ii1 = is-1, is+1 ! do jj1 = js, js+1 do ii1 = is-searchradi, is+searchradi do jj1 = js, js+searchradi if (jj1 >= 1 .and.jj1 <= nlat ) then i1 = ii1 ; j1 = jj1 if (i1 > nlon) i1 = i1-nlon if (i1 < 1) i1 = nlon + i1 if(i1 == is .and. j1 == js) then !---do nothing else ii = i1 jj = j1 if(look(i1,j1) == 0 .and. pos(i1,j1) == 1) then !here we've found the next grid pt. under the trough axis nn = nn + 1 ; ic=ic+1 ; look(i1,j1) = 1 ! print*,na,nn,ii,jj xlat(na,nn) = latArray(ii,jj) xlon(na,nn) = lonArray(ii,jj) xvor(na,nn) = vorArray(ii,jj) goto 330 !go back to look for additional grids under this trough axis endif !do that reasigning the seed point to the current trough pt. endif endif enddo !do jj1 = js, js+1 enddo !do ii1 = is-1, is+1 ! curtail troughs less than required points long mingrid=nint(minlength/((((nglon_int+nglat_int)*0.5)*111))) if(nn.lt.mingrid) then na=na-1 endif endif !if(look(g,h) == 0 .and. pos(g,h) == 1) enddo !do g=gL,gR enddo !do h=hB,hT print *, "Number of troughs found = ", na, ic write(12,*) datec,'Number of troughs: ',na ! determine which trough have some part in out box of interest do k=1,na loccheck(k)=0 enddo do k=1,na do i=1,axiscount(k) if(xlon(k,i) >= 180.) then xlon(k,i)=xlon(k,i)-360. endif if((xlon(k,i) > lonL).and.(xlon(k,i) < lonR)) then loccheck(k) = 1 endif enddo enddo etrough=na do k=1,na if (loccheck(k).eq.1) then write(12,*) 'Trough number',k,'lat/lon/int' ! print*, 'Trough number',k,'lat/lon/int' do i=1,axiscount(k) write(12,967)xlat(k,i),xlon(k,i),xvor(k,i) ! write(*,967)xlat(k,i),xlon(k,i),xvor(k,i) enddo write(12,*)'' endif enddo 967 format (f10.3,1x,f10.3,1x,f12.7) print*,'TROUGHS DONE ',datec deallocate (pos,look,latArray,lonArray,vorArray) deallocate (xlat,xlon,xvor,axisCount,loccheck) return end subroutine trough_e ! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: !-------------------------------------------------------- subroutine jet_e(ll_lat,ll_lon,nglon_int,nglat_int,axisdata,nlon,nlat,& &latB,latT,lonL,lonR,data,minlength,datec,searchrad) !!$ Same ideas as above ! Code adapted from Anantha's original real :: ll_lat,ll_lon,nglon_int,nglat_int,latB,latT,lonL,lonR,x1 integer :: nlon,nlat,jcounte,hB,hT,gL,gR,grange,hrange,req integer :: g,h real, dimension (nlon,nlat) :: axisdata,data real, dimension (nlon*nlat,3) :: Jloce character (len=20) :: datec real, dimension(:,:), allocatable :: latArray, lonArray, vorArray real, dimension(:,:), allocatable :: xLat,xLon,xVor integer, dimension(:), allocatable :: axisCount,loccheck real, dimension(:,:), allocatable :: pos, look integer :: i,j,k,ii,jj,is,js,jj1,ii1,i1,j1,ip1 integer :: na,ic,nn,etrough,ejet real :: lat0,lat1,lon0,lon1 real :: minlength,crosslon,crosslat integer :: mingrid,meridflag tcounte=0 hB = ((nint((latB-ll_lat)/nglat_int))+1) hT = ((nint((latT-ll_lat)/nglat_int))+1) hrange = hT-hB do g=1,nlon do h=hB,hT if(axisdata(g,h).ne.amiss)then if (abs(axisdata(g,h)) <= epsilon(1.)) then ! zero line exactly at grid pt jcounte = jcounte+1 Jloce(jcounte,1)=ll_lon + (g-1)*nglon_int !lon Jloce(jcounte,2)=ll_lat + (h-1)*nglat_int !lat Jloce(jcounte,3)=data(g,h) else ! interpolate in lon to find zero line if (axisdata(g,h) > 0. .and. axisdata(g,h+1) < 0.) then x1 = ll_lat + (h-1)*nglat_int jcounte = jcounte+1 Jloce(jcounte,1)=ll_lon + (g-1)*nglon_int !lon Jloce(jcounte,2)=(x1 + (nglat_int*(axisdata(g,h)/(axisdata(g,h)-axisdata(g,h+1))))) !lat Jloce(jcounte,3)=max(data(g,h),data(g+1,h)) endif endif endif enddo enddo print*,'NUMBER OF JET POINTS:' print*,'' print*,jcounte,'EASTERLY' ! now the axis finder routine - easterly open (unit=12, file='jet_east.out') allocate ( xlat(jcounte,1000) ) allocate ( xlon(jcounte,1000) ) allocate ( xvor(jcounte,1000) ) allocate ( axisCount(jcounte) ) allocate ( loccheck(jcounte) ) allocate (pos(nlon,nlat) ) allocate (look(nlon,nlat) ) allocate (latArray(nlon,nlat)) allocate (lonArray(nlon,nlat)) allocate (vorArray(nlon,nlat)) searchradi=nint(searchrad/((((nglon_int+nglat_int)*0.5)*111))) print*,'SEARCH',searchrad,searchradi do g=1,nlon do h=1,nlat pos(g,h) = 0 look(g,h) = 0 enddo enddo do k = 1, jcounte h = 1 + nint( (jloce(k,2)-ll_lat) / nglat_int ) g = 1 + nint( (jloce(k,1)-ll_lon) / nglon_int ) pos(g,h) = 1 latArray(g,h) = jloce(k,2) lonArray(g,h) = jloce(k,1) vorArray(g,h) = jloce(k,3) enddo do g=1,nlon do h=1,nlat if(pos(g,h).eq.0.)then look(g,h) = 1 endif enddo enddo ic = 0 ; na = 0 !counter to track the individual trough axes ! Concentrate within the domain of interest do g=1,nlon do h=hB,hT nn = 0 !counter to track grids under a given trough if(look(g,h) == 0 .and. pos(g,h) == 1) then !!$!---we have found a fresh trough axis! na = na + 1 !troughs nn = nn + 1 !points ic=ic+1 xlat(na,nn) = latArray(g,h) xlon(na,nn) = lonArray(g,h) xvor(na,nn) = vorArray(g,h) look(g,h) = 1 ii = g ; jj = h !temporarily store the location of the axis begin pt. 330 continue !the previous point now is the seed point to look for the rest of the axis is = ii js = jj axisCount(na) = nn !we look in the vicinity of the seed point for the remainder of the trough ! do ii1 = is-1, is+1 ! do jj1 = js, js+1 do ii1 = is-searchradi, is+searchradi do jj1 = js, js+searchradi if (jj1 >= 1 .and.jj1 <= nlat ) then i1 = ii1 ; j1 = jj1 if (i1 > nlon) i1 = i1-nlon if (i1 < 1) i1 = nlon + i1 if(i1 == is .and. j1 == js) then !---do nothing else ii = i1 jj = j1 if(look(i1,j1) == 0 .and. pos(i1,j1) == 1) then !here we've found the next grid pt. under the trough axis nn = nn + 1 ; ic=ic+1 ; look(i1,j1) = 1 ! print*,na,nn,ii,jj xlat(na,nn) = latArray(ii,jj) xlon(na,nn) = lonArray(ii,jj) xvor(na,nn) = vorArray(ii,jj) goto 330 !go back to look for additional grids under this trough axis endif !do that reasigning the seed point to the current trough pt. endif endif enddo !do jj1 = js, js+1 enddo !do ii1 = is-1, is+1 ! curtail troughs less than required points long mingrid=nint(minlength/((((nglon_int+nglat_int)*0.5)*111)))*2 if(nn.lt.mingrid) then na=na-1 endif endif !if(look(g,h) == 0 .and. pos(g,h) == 1) enddo !do g=gL,gR enddo !do h=hB,hT print *, "Number of jets found = ", na, ic do k=1,na loccheck(k)=0 enddo do k=1,na do i=1,axiscount(k) if(xlon(k,i) >= 180.) then xlon(k,i)=xlon(k,i)-360. endif if((xlon(k,i) > lonL).and.(xlon(k,i) < lonR)) then loccheck(k) = 1 endif enddo enddo write(12,*) datec,'Number of jets: ',na etrough=na do k=1,na if (loccheck(k).eq.1) then write(12,*) 'Jet number',k,'lat/lon/int' ! print*, 'Jet number',k,'lat/lon/int' do i=1,axiscount(k) write(12,967)xlat(k,i),xlon(k,i),xvor(k,i) ! write(*,967)xlat(k,i),xlon(k,i),xvor(k,i) enddo write(12,*)'' endif enddo 967 format (f10.3,1x,f10.3,1x,f12.7) print*,'JETS DONE ',datec deallocate (pos,look,latArray,lonArray,vorArray) deallocate (xlat,xlon,xvor,axisCount) return end subroutine jet_e