subroutine garethtrough(u,v,dim1,dim2,ll_lat,ll_lon,ur_lat,ur_lon,& &globflag,smoothpas,vortmin,ucutoff,ucutoff2,advcu,aej,advcuw,aejw,vor,shvor,cuvor) ! dim1 is the fastest varying (fortran issue) - will assume longitude ! this version uses variables in the order latxlon, reversed from orig. ! on mca complie using: ! WRAPIT -g95 -L /opt/local/lib trough.stub sub.trough.f90 implicit none integer :: dim1,dim2 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,amiss integer :: g,h,iret,uncalc integer :: imain, igem1, igem2, icoord, iret1, coord, globflag real :: deltax,deltay,radius,pi,vortmin,ucutoff,ucutoff2 real,dimension (dim1,dim2) :: u,v,vor,p real,dimension (dim1,dim2) :: shvor,cuvor,advcu,aej,fulladv real,dimension (dim1,dim2) :: ddxvor,ddyvor,aejcond,magwnd real,dimension (dim1,dim2) :: ddxcu,ddycu,ddxadv,ddyadv,mask real,dimension (dim1,dim2) :: advcuw,aejw real :: ll_lat,ll_lon,ur_lat,ur_lon,nglat_int,nglon_int ! real :: ll_lat2,ll_lon2,ur_lat2,ur_lon2 integer :: smoothpas ! fixed parameters parameter (radius=6371000) !radius of earth in metres pi = 4.*atan(1.) amiss = -9999.0 dnopts = 1 nglat_int=(ur_lat - ll_lat)/(dim2-1) nglon_int=(ur_lon - ll_lon)/(dim1-1) print*,'lat interval =',nglat_int,'lon interval = ',nglon_int ! PRE-SMOOTHING print*,'DIM1 = ',dim1 print*,'DIM2 = ',dim2 call smooth ( u,dim1,dim2,amiss,smoothpas ) call smooth ( v,dim1,dim2,amiss,smoothpas ) uncalc = int(dnopts+1) ! --------------- CALCULATE VORTICITY --------------- ! note that having to lose another point near poles ! set the non-calculatable values to -9999.9 do g=1,dim1 do h=1,dim2 if ((h.le.uncalc).or.(h.gt.(dim2-uncalc)))then vor(g,h)=amiss else if ((globflag.eq.0).and.((g.le.uncalc).or.& &(g.gt.(dim1-uncalc))))then vor(g,h)=amiss else !calculate deltax, deltay call calc_dxdy(dim1,dim2,g,h,globflag,ll_lat,nglon_int,nglat_& &int,deltax,deltay) ! for points at the edges: m=h+dnopts n=h-dnopts s=g+dnopts t=g-dnopts ! account for split in global datasets (non-global already broken out) if (s.gt.dim1)then s=s-dim1 endif if (t.le.0) then t=dim1+t endif ! vortcity calculation vor(g,h)=((v(s,h)-v(t,h))/((2*dnopts)*deltax))-((u(g,m)-u(g,n))/((2*dnopts)*deltay)) ! multiply vortcity by 10^5 to make it easier to handle: vor(g,h)=vor(g,h)*100000. endif enddo enddo ! --------------- Calculate Shear vorticity ------------- !set the non-calculatable values to -9999.9 do g=1,dim1 do h=1,dim2 if ((h.le.uncalc).or.(h.gt.(dim2-uncalc)))then shvor(g,h)=amiss else if ((globflag.eq.0).and.((g.le.uncalc).or.& &(g.gt.(dim1-uncalc))))then shvor(g,h)=amiss else ! calculate deltax, deltay call calc_dxdy(dim1,dim2,g,h,globflag,ll_lat,nglon_int,nglat_& &int,deltax,deltay) ! for point n of ref ! c ! b a ! d !General points: m=h+dnopts n=h-dnopts s=g+dnopts t=g-dnopts ! account for split in global datasets (non-global already broken out) if (s.gt.dim1)then s=s-dim1 endif if (t.le.0) then t=dim1+t 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(t,h) vee=v(t,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(s,h) vee=v(s,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,m) vee=v(g,m) 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,n) vee=v(g,n) 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 !---------- CURVATURE VORTICITY ------------------ !Calculate curvature vort as residual of full and shear vort !set the non-calculatable values to -9999.9 do g=1,dim1 do h=1,dim2 if ((h.le.uncalc).or.(h.gt.(dim2-uncalc)))then cuvor(g,h)=amiss else if ((globflag.eq.0).and.((g.le.uncalc).or.& &(g.gt.(dim1-uncalc))))then cuvor(g,h)=amiss else ! all the other points: cuvor(g,h)=vor(g,h)-shvor(g,h) endif enddo enddo ! ------------ Advection of curvature vort -------------- ! set the non-calculatable values to -9999.9 do g=1,dim1 do h=1,dim2 if ((h.le.(uncalc+1)).or.(h.gt.(dim2-(uncalc+1))))then advcu(g,h)=amiss advcuw(g,h)=amiss else if ((globflag.eq.0).and.((g.le.(uncalc+1)).or.& &(g.gt.(dim1-(uncalc+1))))) then advcu(g,h)=amiss advcuw(g,h)=amiss else ! calculate deltax, deltay call calc_dxdy(dim1,dim2,g,h,globflag,ll_lat,nglon_int,nglat_& &int,deltax,deltay) ! deriavtives m=h+dnopts n=h-dnopts s=g+dnopts t=g-dnopts ! account for split in global datasets (non-global already broken out) if (s.gt.dim1)then s=s-dim1 endif if (t.le.0) then t=dim1+t endif ! simple centred difference ddxcu(g,h)=(cuvor(s,h)-cuvor(t,h))/(2*dnopts*deltax) ddycu(g,h)=(cuvor(g,m)-cuvor(g,n))/(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 ! MASKING VARIABLE - REMOVE TROUGHS WHERE V(str).(del)advcu not >0 ! set the non-calculatable values to -9999.9 do g=1,dim1 do h=1,dim2 if ((h.le.(uncalc+2)).or.(h.gt.(dim2-(uncalc+2))))then mask(g,h)=amiss else if ((globflag.eq.0).and.((g.le.(uncalc+2)).or.& &(g.gt.(dim1-(uncalc+2))))) then mask(g,h)=amiss else ! calculate deltax, deltay call calc_dxdy(dim1,dim2,g,h,globflag,ll_lat,nglon_int,nglat_& &int,deltax,deltay) ! General points: m=h+dnopts n=h-dnopts s=g+dnopts t=g-dnopts ! account for split in global datasets (non-global already broken out) if (s.gt.dim1)then s=s-dim1 endif if (t.le.0) then t=dim1+t endif ! calc derivatives: ddxadv(g,h)=(advcu(s,h)-advcu(t,h))/(2*dnopts*deltax) ddyadv(g,h)=(advcu(g,m)-advcu(g,n))/(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 ! -------- TROUGH AXES (AT LAST) ---------- ! dump non +ve points if(globflag.eq.1)then do g=1,dim1 do h=(uncalc+2),(dim2-(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),(dim1-(uncalc+2)) do h=(uncalc+2),(dim2-(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 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,dim1 do h=1,dim2 if ((h.le.(uncalc+2)).or.(h.gt.(dim2-(uncalc+2))))then aej(g,h)=amiss aejw(g,h)=amiss else if ((globflag.eq.0).and.(g.le.(uncalc+2)).or.(g.gt.(dim1-(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(dim1,dim2,g,h,globflag,ll_lat,nglon_int,nglat_& &int,deltax,deltay) ! General points: m=h+dnopts n=h-dnopts s=g+dnopts t=g-dnopts ! account for split in global datasets (non-global already broken out) if (s.gt.dim1)then s=s-dim1 endif if (t.le.0) then t=dim1+t endif ! calc derivatives: ddxvor(g,h)=((shvor(s,h)-shvor(t,h))/(2*dnopts*deltax)) ddyvor(g,h)=((shvor(g,m)-shvor(g,n))/(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 999 continue 555 continue return end subroutine garethtrough subroutine calc_dxdy(dim1,dim2,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 :: dim2,dim1,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(dim2-1))*pi deltax=(((sin(xgridpt))*radius)*(2*pi))/float(dim1) deltay=(pi*radius)/float(dim2-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 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_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 smooth(indata,nxdim,nydim,amiss,nsmooth) ! Anantha's smoothing routine implicit none integer :: nxdim,nydim,nsmooth,i,j,ii real, dimension (nxdim,nydim) :: indata real, dimension (:,:), allocatable :: Ldat,Rdat,Tdat,Bdat,Sdat,ic,& iarray real :: amiss allocate ( Ldat(nxdim,nydim) ) ; allocate ( Rdat(nxdim,nydim) ) allocate ( Tdat(nxdim,nydim) ) ; allocate ( Bdat(nxdim,nydim) ) allocate ( Sdat(nxdim,nydim) ) ; allocate ( ic(nxdim,nydim) ) allocate ( iarray(nxdim,nydim)) do ii = 1,nsmooth ! print *, "Now smoothing", nxdim, nydim, amiss Tdat=amiss ; Bdat = amiss ; Ldat=amiss ; Rdat = amiss Ldat(1:nxdim-1,:) = indata(2:nxdim,:) Rdat(2:nxdim,:) = indata(1:nxdim-1,:) Tdat(:,2:nydim) = indata(:,1:nydim-1) Bdat(:,1:nydim-1) = indata(:,2:nydim) 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 ( indata == amiss ) iarray = 0 end where Sdat = Sdat + 2.*indata*iarray ic = ic + 2.*iarray indata = amiss where ( abs(ic) > epsilon(1.) ) indata = sdat/ic end where !!$ !!$ enddo deallocate (Ldat,Rdat,Tdat,Bdat,Sdat,ic,iarray) return end subroutine smooth