*----------------------------------------------------------------------- * uvVECT *----------------------------------------------------------------------- SUBROUTINE uvVECT(U,MU,V,MV,NX,NY) parameter ( iundef = 9999 ) REAL U(MU,*),V(MV,*) real uxp(1000), uyp(1000) CHARACTER CMSG*80 LOGICAL LMISS,LNRMAL,LSETX,LSETY,LMSG,LEQRAT, + LMISSP,LSMALL,LUNIT,LUMSG logical ldeg logical linvis equivalence ( vrcnst, ivcnst ) SAVE *----------------------------------------------------------------------- * / GET INTERNAL PARAMETERS / *----------------------------------------------------------------------- CALL GLPGET('LMISS ',LMISS ) CALL GLPGET('RMISS ',RMISS ) call glpget('REALMAX ',REMAX ) call glpget('REALMIN ',REMIN ) CALL uvPGET('INDEX ',INDEX ) CALL uvPGET('LNRMAL ',LNRMAL) CALL uvPGET('LEQRAT ',LEQRAT) CALL uvPGET('LMSG ',LMSG ) CALL uvPGET('LUNIT ',LUNIT ) CALL uvPGET('LUMSG ',LUMSG ) CALL uvPGET('ICENT ',ICENT ) CALL uvPGET('LMISSP ',LMISSP) CALL uvPGET('ITYPE1 ',ITYPE1) CALL uvPGET('LSMALL ',LSMALL) CALL uvPGET('RSMALL ',RSMALL) CALL uvPGET('ITYPE2 ',ITYPE2) CALL uvPGET('RSIZEM ',RSIZEM) CALL uvPGET('RSIZET ',RSIZET) CALL uvPGET('XTTL ',XTTL ) CALL uvPGET('IXINT ',IXINT ) CALL uvPGET('IYINT ',IYINT ) call uvpget('linvis ',linvis ) call uvpget('vvpmax ',vvpmax ) call uvpget('iarlen ',iarlen ) call uvpget('vrcnst ',vrcnst ) *----------------------------------------------------------------------- * / get info. on coordinate and mapping / *----------------------------------------------------------------------- call sgpget ( 'LDEG', ldeg ) if ( ldeg ) then cp = 3.14159265 / 180.0 else cp = 1.0 end if call sgqtrn ( itr ) call sgqwnd ( uxmn, uxmx, uymn, uymx ) call sgqvpt ( vxmn, vxmx, vymn, vymx ) if ( itr .ge. 5 .and. itr .le. 7 ) then call sgqsim ( aotg, vxoff, vyoff ) else if ( itr .ge. 10 .and. itr .le. 13 ) then call sgqsim ( radius, vxoff, vyoff ) call sgqmpl ( xlon0, ylat0, rot0 ) else if ( itr .ge. 20 .and. itr .le. 23 ) then call sgqsim ( radius, vxoff, vyoff ) call sgqmpl ( xlon0, ylat0, rot0 ) else if ( itr .ge. 30 .and. itr .le. 33 ) then call sgqsim ( radius, vxoff, vyoff ) call sgqmpl ( xlon0, ylat0, rot0 ) end if if ( itr .ge. 2 .and. itr .le. 4 ) then CMSG='LOG coodinate not supported' CALL MSGDMP('E','uvVECT',CMSG) end if if ( itr .eq. 1 .and. iarlen .eq. 3 ) then CMSG=' iarlen=3 is not good with itr=1 ' CALL MSGDMP('W','uvVECT',CMSG) end if *----------------------------------------------------------------------- * / get ( and set if necessary ) info. on grid points / *----------------------------------------------------------------------- * / SET GRID ATTRIBUTE IF IT HAS NOT BEEN SET YET / CALL UWQGXZ(LSETX) CALL UWQGYZ(LSETY) IF (.NOT.LSETX) THEN CALL UWSGXB(UXMN,UXMX,NX) CALL UWSGXZ(.FALSE.) END IF IF (.NOT.LSETY) THEN CALL UWSGYB(UYMN,UYMX,NY) CALL UWSGYZ(.FALSE.) END IF * / get grid position on WC / call uwqgxa ( uxp, nx ) call uwqgya ( uyp, ny ) *----------------------------------------------------------------------- * / CHECK INPUT DATA ( missing, maximum, minimum ) / *----------------------------------------------------------------------- ndata = 0 rminu = REMAX rminv = REMAX rmaxu = REMIN rmaxv = REMIN DO 10 J=1,NY,IYINT DO 10 I=1,NX,IXINT if( .NOT. ( LMISS .AND. & (U(I,J).EQ.RMISS .OR. V(I,J).EQ.RMISS)) ) then ndata = ndata + 1 RMAXU = max ( U(I,J), rmaxu ) RMAXV = max ( V(I,J), rmaxv ) RMINU = min ( U(I,J), rminu ) RMINV = min ( V(I,J), rminv ) end if 10 CONTINUE * / WRITE MESSAGE and return IF MISSING FIELD / IF ( ndata .eq. 0 ) then CMSG='MISSING FIELD.' CALL MSGDMP('W','uvVECT',CMSG) IF (LMSG) THEN call uvwmsg( cmsg ) END IF RETURN * -------> exit END IF * / WRITE MESSAGE and return ZERO FIELD / IF ( RMINU.EQ.0 .AND. RMAXU.EQ.0 + .AND. RMINV.EQ.0 .AND. RMAXV.EQ.0 ) THEN CMSG='ZERO FIELD.' CALL MSGDMP('W','uvVECT',CMSG) IF (LMSG) THEN call uvwmsg( cmsg ) END IF RETURN * -------> exit END IF *----------------------------------------------------------------------- * / CALCULATE NORMALIZATION FACTOR if necessary / *----------------------------------------------------------------------- IF (LNRMAL) THEN UE=(VXMX-VXMN)/(NX/IXINT) VE=(VYMX-VYMN)/(NY/IYINT) if ( itr .le. 4 ) then AU=RGNLE(UE/MAX(ABS(RMINU),ABS(RMAXU))) AV=RGNLE(VE/MAX(ABS(RMINV),ABS(RMAXV))) IF (LEQRAT) THEN AU=MIN(AU,AV) AV=au END IF else if ( itr.ge.5 .and. itr.le.7 ) then if ( leqrat ) then rumax = max( max( abs(rminu), abs(rmaxu) ), & max( abs(rminv), abs(rmaxv) ) ) aa = sqrt ( ue * ve ) / rumax au = rgnle(aa) av = au else cmsg = 'LEQRAT must be TRUE if 5 <= ITR <= 7 ' call msgdmp('W','uvVECT',cmsg) return * -------> exit end if else if ( itr.ge.10 .and. itr.le.33 ) then if ( leqrat ) then rumax = max( max( abs(rminu), abs(rmaxu) ), & max( abs(rminv), abs(rmaxv) ) ) aa = 3.1416 * sqrt ( ue * ve ) * radius / rumax au = rgnle(aa) av = au else cmsg = 'LEQRAT must be TRUE if 10<=ITR<=33 ' call msgdmp('W','uvVECT',cmsg) return * -------> exit end if end if ELSE CALL uvPGET('XFACT1 ',XFACT1) CALL uvPGET('YFACT1 ',YFACT1) AU=XFACT1 AV=YFACT1 END IF CALL uvPSET('XFACT2 ',AU) CALL uvPSET('YFACT2 ',AV) *----------------------------------------------------------------------- * / output info. on scaling / *----------------------------------------------------------------------- * / determine length of the UNIT VECTOR and draw them / IF (LUNIT) THEN CALL uvUNIT * / WRITE UNIT VALUE / IF (LUMSG) THEN if ( iarlen .eq. 1 .or. iarlen .eq. 3 ) then CALL uvPGET('UXUNIT ',UXUNIT) CALL uvPGET('UYUNIT ',UYUNIT) CMSG='XUNIT =##########, YUNIT =##########' WRITE(CMSG( 8:17),'(1P,E10.3)') UXUNIT WRITE(CMSG(27:36),'(1P,E10.3)') UYUNIT else if ( iarlen .eq. 2 ) then cmsg = ' XUNIT = null, YUNIT = null ' end if call uvwmsg( cmsg ) END IF * / DRAW TITLE FOR UNIT VECTOR / CALL uvDUT END IF * / WRITE SCALING FACTOR / IF (LMSG) THEN if ( iarlen .eq. 1 .or. iarlen .eq. 3 ) then CMSG='XFACT =##########, YFACT =##########' WRITE(CMSG( 8:17),'(1P,E10.3)') AU WRITE(CMSG(27:36),'(1P,E10.3)') AV else if ( iarlen .eq. 2 ) then cmsg = 'XFACT = NULL , YFACT = NULL' end if call uvwmsg ( cmsg ) END IF *----------------------------------------------------------------------- * / DRAW VECTORS / *----------------------------------------------------------------------- * / SET line type of arrow / * / SET the length of arrow for direction-only-plot / if ( ivcnst .eq. iundef ) then vrcnst = 0.5 / sqrt( (real(nx)/ixint)*(real(ny)/iyint) ) else call uvpget ( 'VRCNST', vrcnst ) end if * DO 20 J=1,NY,IYINT DO 20 I=1,NX,IXINT * CALL STFTRF( uxp(i), uyp(j), VXSTRT, VYSTRT ) * / skip missing value and mark the point if required / if ( lmiss .and. & (u(i,j).eq.rmiss .or. v(i,j).eq.rmiss) ) then IF (LMISSP) THEN if ( vxstrt .ne. rmiss .and. & vystrt .ne. rmiss ) then CALL SGPMZV & (1,VXstrt,VYstrt,ITYPE1,INDEX,RSIZEM) end if END IF go to 20 end if * / skip if the grid is on singular point of mapping / IF ( LMISS )THEN CALL STITRF( VXSTRT, VYSTRT, UXSTRT, UYSTRT ) IF ( ( UXSTRT .EQ. RMISS ) .OR. ( UYSTRT .EQ. RMISS ) )THEN go to 20 ENDIF end if * / skip far-side point if orthographic projection / if ( itr.eq.30 ) then xradc = xlon0 * cp yradc = 3.14159265/2.0 - ylat0 * cp xradp = uxp(i) * cp yradp = 3.14159265/2.0 - uyp(j) * cp rtemp = 1.0 call ct3sc( rtemp, yradc, xradc, xc, yc, zc ) call ct3sc( rtemp, yradp, xradp, xp, yp, zp ) zz = xc*xp + yc*yp + zc*zp if ( zz .lt. 0.0 ) then go to 20 end if end if * / get unit (x,y)-tangent vector in NDC / call uvtang I ( uxp(i), uyp(j), itr, cp, O veux, veuy, vevx, vevy ) * / skip singular point / if ( lmiss .and. & ( veux .eq. rmiss .or. veuy .eq. rmiss .or. & vevx .eq. rmiss .or. vevy .eq. rmiss ) ) then go to 20 end if * / vector arrow component in NDC / if ( itr.ge.10 .and. itr.le.13 ) then vdx = ( u(i,j) * au * veux + v(i,j) * av * vevx ) / radius vdy = ( u(i,j) * au * veuy + v(i,j) * av * vevy ) / radius else if ( itr.ge.20 .and. itr.le.23 ) then vdx = ( u(i,j) * au * veux + v(i,j) * av * vevx ) / radius vdy = ( u(i,j) * au * veuy + v(i,j) * av * vevy ) / radius else if ( itr.ge.30 .and. itr.le.33 ) then vdx = ( u(i,j) * au * veux + v(i,j) * av * vevx ) / radius vdy = ( u(i,j) * au * veuy + v(i,j) * av * vevy ) / radius else if ( itr.ge.5 .and. itr.le.7 ) then vdx = ( u(i,j) * au * veux + v(i,j) * av * vevx ) / aotg vdy = ( u(i,j) * au * veuy + v(i,j) * av * vevy ) / aotg else vdx = ( u(i,j) * au * veux + v(i,j) * av * vevx ) vdy = ( u(i,j) * au * veuy + v(i,j) * av * vevy ) end if if ( iarlen .eq. 1 ) then tabs = sqrt( (au*u(i,j))**2 + (av*v(i,j))**2 ) vabs = sqrt( vdx**2 + vdy**2 ) vdx = tabs * vdx / vabs vdy = tabs * vdy / vabs else if ( iarlen .eq. 2 ) then vabs = sqrt(vdx**2 + vdy**2) if ( vabs .ne. 0 ) then vdx = vrcnst * vdx / vabs vdy = vrcnst * vdy / vabs end if else if ( iarlen .eq. 3 ) then continue else CMSG='iarlen not implemented' CALL MSGDMP('E','uvvect',CMSG) end if * / skip if vector is too long or too short / rlvec2 = vdx**2 + vdy**2 if ( rlvec2 .gt. vvpmax**2 ) then CALL SGPMZV(1,VXstrt,VYstrt,ITYPE2,INDEX,RSIZEM) go to 20 else IF (LSMALL .AND. rlvec2.LE.RSMALL**2) THEN CALL SGPMZV(1,VXstrt,VYstrt,ITYPE2,INDEX,RSIZEM) go to 20 END IF * / do centering of arrow, and plot it ! / VX1 = VXstrt - vdx*(ICENT+1)/2.0 VX2 = VXstrt - vdx*(ICENT-1)/2.0 VY1 = VYstrt - vdy*(ICENT+1)/2.0 VY2 = VYstrt - vdy*(ICENT-1)/2.0 c$$$ CALL UFLNZV(VX1,VY1,VX2,VY2,INDEX) CALL SGLAZV(VX1,VY1,VX2,VY2,1,INDEX) 20 CONTINUE return END