SUBROUTINE TRANSCOM C**** C**** THIS SUBROUTINE ACCUMULATES SEVERAL DIAGNOSTIC ARRAYS C**** FOR TRANSCOM2 EXPERIMENT. C**** ALL ARRAYS ARE INTERPOLATED TO PREFERRED PRESSURE LEVELS C**** BEFORE TIME AVERAGING C**** Pressure Levels for TRANSCOM2: C**** 50,100,150,200,250,300,(350),400,(450),500,(550), C**** (600),(650),700,(750),(800),850, C**** (875),(900),(925),(950),(975),1000,(1010),(1020) C C**** PARAMETER (IM=72,JM=46,LM=9,MTRACE=1) c** the following are defined or computed in tracer model. REAL*4 M COMMON/PASS/ NTRACE,GRAV,IDACC(8),SIGE(LM+1),SIG(LM),DSIG(LM), * DXYP(JM),M(IM,JM,LM),RM(IM,JM,LM,MTRACE),SRC(IM,JM,MTRACE) DIMENSION PU(IM,JM,LM),PV(IM,JM,LM),P(IM,JM),SD(IM,JM,LM-1) c** transcom arrays INCLUDE 'TRANSCOM.COM' COMMON/STIJ/IS(NSTN),JS(NSTN),LS(NSTN) C--------------------------------------------------------------- ENTRY TRANSI C C** one-time initialization of these quantities C** call in subroutine INPUT (eg at line 606.5) c** unit number for station data KSTN=60 if(kstn.le.0)then write(6,*)' resetting unit number to 60' kstn=60 end if c** NDT value from rsf will take precedence NDT=0 C** scale factor to go from kg SF6/kg air to mole fraction C*** 1 mole SF6 = 146e-3 kg C** UNIT=29./146.*1.e12 c** c0 is initial mass mixing ratio of SF6 source c** =1.04E-11 kg SF6/kg air C0=1.04E-11 c** pressure levels plev(1)=1020.0 plev(2)=1010.0 plev(3)=1000.0 plev(4)=975.0 plev(5)=950.0 plev(6)=925.0 plev(7)=900.0 plev(8)=875.0 plev(9)=850.0 plev(10)=800.0 plev(11)=750.0 plev(12)=700.0 plev(13)=650.0 plev(14)=600.0 plev(15)=550.0 plev(16)=500.0 plev(17)=450.0 plev(18)=400.0 plev(19)=350.0 plev(20)=300.0 plev(21)=250.0 plev(22)=200.0 plev(23)=150.0 plev(24)=100.0 plev(25)=50.0 C C** alternate cbkg, background mixing ratio, needs gsrc C** gsrc = globally integrated source since time zero (in kg/s) c** value from rsf will take precedence C write(6,*)' ntrace = ',ntrace DO 20 N=1,NTRACE GSRC(N)=0.0 20 CONTINUE C**** FIND STATION LOCATIONS (igrid=2 for fgrid) IGRID=2 CALL STNIJ(IGRID) SKIP=0.0 RETURN C--------------------------------------------------------------- ENTRY TRANSIM C C** initialize these quantities at start of each month C** call in MAIN (at line 70.5) C IDTRANS=0 DO 30 J=1,JM DO 30 I=1,IM APS(I,J)=0.0 30 CONTINUE DO 50 K=1,NPL DO 50 J=1,JM DO 50 I=1,IM ABETA(I,J,K)=0.0 ABETAC(I,J,K)=0.0 ABETAD(I,J,K)=0.0 AUP(I,J,K)=0.0 AVP(I,J,K)=0.0 AWP(I,J,K)=0.0 DO 40 N=1,NTRACE ACHIP(I,J,K,N)=0.0 ACNP(I,J,K,N)=0.0 AUCNP(I,J,K,N)=0.0 AVCNP(I,J,K,N)=0.0 AWCNP(I,J,K,N)=0.0 ARMCNVP(I,J,K,N)=0.0 ARMDIFP(I,J,K,N)=0.0 40 CONTINUE 50 CONTINUE RETURN C--------------------------------------------------------------- ENTRY GETUVW(p,pu,pv,sd) C C** Tracer model uses PU, PV and SD(=PW). C** get U,V,W from PU,PV,PW at each timestep c** this should be called at end of or right after DYNAM0 (1105.5) c** c usual setup: if p(i,j) is at center of box, c u(i,j) and v(i,j) are at bottom right corner c | | c ---u,v--------------pv--------------u,v-------- c | | | c | | dyp(j) c pu p(i,j) pu(i,j) | c |<- - - - - - dxp(j)- - -- - - ->| | c | | | c ---u,v-------------pv(i,j)----------u,v(i,j)--- c |<- - - - - - dxv(j) - - - - - ->| C c write(6,*)' in getuvw' c** reassign sd to pw array do 90 l=1,lm-1 do 90 j=1,jm imax=im if(j.eq.1.or.j.eq.jm)imax=1 do 80 i=1,imax pw(i,j,l)=sd(i,j,l) 80 continue 90 continue c** temporarily fill pw (l=lm) with zeros do 110 j=1,jm imax=im if(j.eq.1.or.j.eq.jm)imax=1 do 100 i=1,imax pw(i,j,lm)=0.0 100 continue 110 continue C**** FILL IN POLES TO SIMPLIFY CODE DO 120 J=1,JM,JM-1 DO 120 I=2,IM P(I,J)=P(1,J) 120 CONTINUE DO 130 L=1,LM DO 130 J=1,JM,JM-1 DO 130 I=2,IM PU(I,J,L)=PU(1,J,L) PV(I,J,L)=PV(1,J,L) c IF(L.EQ.LM)GO TO 130 PW(I,J,L)=PW(1,J,L) 130 CONTINUE c** determine u at pu location, v at pv location DO 140 L=1,LM DO 140 J=1,JM DO 140 I=1,IM u(i,j,l)=0.0 IP1=I+1 IF(I.EQ.IM) IP1=1 PAVG=0.5*(P(I,J)+P(IP1,J)) if(pavg*dyp(j).ne.0.0)U(I,J,L)=PU(I,J,L)/(PAVG*DYP(J)) 140 CONTINUE DO 150 L=1,LM DO 150 J=1,JM DO 150 I=1,IM v(i,j,l)=0.0 IF (J.EQ.1) THEN PM1=0. ELSE PM1=P(I,J-1) ENDIF c*** unequal weighting for polar pressures: PAVG=0.5*(P(I,J)+PM1) IF (J.EQ.JM) PAVG=(2.*P(I,JM)+4.*P(I,JM-1))/6. IF (J.EQ.2 ) PAVG=(2.*P(I,1 )+4.*P(I,2 ))/6. if(pavg*dxvsav(j).ne.0.0)V(I,J,L)=PV(I,J,L)/(PAVG*DXVSAV(J)) 150 CONTINUE c** PW is at center of box c** W is at center of box, but at different sigma level DO 160 L=1,LM-1 DO 160 J=1,JM DO 160 I=1,IM w(i,j,l)=0.0 if(p(i,j).ne.0.0)W(I,J,L)=PW(I,J,L)/(P(I,J)*dxyp(j)) 160 CONTINUE c** save sfc pressure for later diagnostics (TRANSIA) c** also determine terrain mask (needed in CONVEC, DIFFUS, TRANSIA) DO 180 J=1,JM IMAX = IM IF((J.EQ.1).OR.(J.EQ.JM)) IMAX=1 DO 170 I=1,IMAX PSAVE(I,J)=P(I,J) c** p is surface pressure (at sigma=1, p=ps) DO 170 K=1,NPL IF(PLEV(K).LT.PSAVE(I,J))THEN BETA(I,J,K)=1.0 ELSE BETA(I,J,K)=0.0 END IF 170 CONTINUE 180 CONTINUE c** fill in poles on saved value of P and BETA DO 190 J=1,JM,JM-1 DO 190 I=2,IM PSAVE(I,J)=PSAVE(1,J) DO 185 K=1,NPL BETA(I,J,K)=BETA(1,J,K) 185 CONTINUE 190 CONTINUE RETURN C--------------------------------------------------------------- ENTRY TRANSIA C c** interpolate to pressure levels and accumulate at each timestep c** call right after DIAGA (156.5) IDTRANS=IDTRANS+1 write(6,*)' in transia: idtrans idacc(8) ',idtrans,idacc(8) c c** keep track of global Ps (used to determine cbkg) c** be careful: P is in work1 and lost after DYNAM0 c HBYG=100./GRAV PSUM=0.0 DO 250 J=1,JM IMAX = IM IF((J.EQ.1).OR.(J.EQ.JM)) IMAX=1 DO 240 I=1,IMAX PSUM=PSUM+HBYG*PSAVE(I,J)*DXYP(J) c** global integrated source in kg DO 230 N=1,NTRACE GSRC(N)=GSRC(N)+SRC(I,J,N)*DXYP(J)*NDIAG*3600. 230 CONTINUE 240 CONTINUE 250 CONTINUE write(6,*)' gsrc (kg) ',(gsrc(n),n=1,ntrace) C**** INTERPOLATE U,V,W TO PRESSURE LEVELS C** be careful PU,PV,SD(=PW) are in work1 and lost after DYNAM0 C** get U,V,W from PU,PV,SD inside or right after DYNAM0 C**** CALL PNTRP(IM,JM,LM,NPL,SIGE,SIG,DSIG,PLEV,SKIP,PSAVE,U,UP) CALL PNTRP(IM,JM,LM,NPL,SIGE,SIG,DSIG,PLEV,SKIP,PSAVE,V,VP) c** W only has lm-1 levels. last level filled with zeros for now CALL PNTRP(IM,JM,LM,NPL,SIGE,SIG,DSIG,PLEV,SKIP,PSAVE,W,WP) C**** C**** ACCUMULATE at pressure levels. C*** Points below ground will have a value of SKIP C*** Use terrain mask to skip points below ground C**** DO 270 J=1,JM IMAX = IM IF((J.EQ.1).OR.(J.EQ.JM)) IMAX = 1 DO 260 I=1,IMAX APS(I,J)=APS(I,J)+PSAVE(I,J) DO 260 K=1,NPL ABETA(I,J,K)=ABETA(I,J,K)+BETA(I,J,K) AUP(I,J,K)=AUP(I,J,K)+UP(I,J,K)*BETA(I,J,K) AVP(I,J,K)=AVP(I,J,K)+VP(I,J,K)*BETA(I,J,K) AWP(I,J,K)=AWP(I,J,K)+WP(I,J,K)*BETA(I,J,K) 260 CONTINUE 270 CONTINUE C**** C**** AIR MASS C**** ASUM=0.0 ASUM2=0.0 DO 290 J=1,JM IMAX = IM IF((J.EQ.1).OR.(J.EQ.JM)) IMAX=1 DO 280 L=1,LM DO 280 I=1,IMAX C** be careful P is in work1 and lost after DYNAM0 c** use M instead. M=hbyg*Ps*dxyp*dsig=Mass of air (kg). c** asum=asum2 ASUM=ASUM+M(I,J,L) ASUM2=ASUM2+HBYG*PSAVE(I,J)*DXYP(J)*DSIG(L) 280 CONTINUE 290 CONTINUE C**** C**** INSTANTANEOUS TRACER CONCENTRATION (kg tracer/kg air) C**** DO 400 N=1,NTRACE TSUM=0.0 TSUM2=0.0 DO 320 J=1,JM IMAX = IM IF((J.EQ.1).OR.(J.EQ.JM)) IMAX = 1 DO 310 L=1,LM DO 310 I=1,IMAX ci(i,j,l)=0.0 if(m(i,j,l).ne.0.0)CI(I,J,L)=RM(I,J,L,N)/M(I,J,L) C** be careful P is in work1 and lost after DYNAM0 c** use RM instead : TSUM=TSUM2 TSUM=TSUM+RM(I,J,L,N) TSUM2=TSUM2+CI(I,J,L)*HBYG*PSAVE(I,J)*DXYP(J)*DSIG(L) 310 CONTINUE 320 CONTINUE c** global mass-weighted mean "background" mixing ratio c** cbkg=cbkg1=cbkg2 if(psum.ne.0.0)CBKG=TSUM2/PSUM if(asum.ne.0.0)CBKG1=TSUM/ASUM c** alternate cbkg: (NOTE RM=0, not C0 on initial start) cc if(asum.ne.0.0)CBKG2=C0+GSRC(N)/ASUM if(asum.ne.0.0)CBKG2=GSRC(N)/ASUM write(6,*)' sum(ps) sum(M) sum(hbyg*ps) ',psum,asum,asum2 write(6,*)' sum(RM) sum(hbyg*CI*ps) ',tsum,tsum2 write(6,*)' tau cbkg cbkg1 cbkg2 ',tau,cbkg,cbkg1,cbkg2 C**** INTERPOLATE CI TO PRESSURE LEVELS before subtracting Cbkg C** needed to compute mole fraction of SF6 CALL PNTRP(IM,JM,LM,NPL,SIGE,SIG,DSIG,PLEV,SKIP,PSAVE, * CI,CIP(1,1,1,N)) c** subtract background mixing ratio at sigma levels DO 350 J=1,JM IMAX = IM IF((J.EQ.1).OR.(J.EQ.JM)) IMAX = 1 DO 330 L=1,LM DO 330 I=1,IMAX CI(I,J,L)=CI(I,J,L)-CBKG cc CI(I,J,L)=CI(I,J,L)-CBKG2 330 CONTINUE 350 CONTINUE C**** C**** CI is now normalized C**** INTERPOLATE TO PRESSURE LEVELS after subtracting Cbkg C**** CALL PNTRP(IM,JM,LM,NPL,SIGE,SIG,DSIG,PLEV,SKIP,PSAVE, * CI,CNP(1,1,1,N)) C**** C**** ACCUMULATE at pressure levels C*** Points below ground will have a value of SKIP C*** Use terrain mask to skip points below ground C**** DO 380 J=1,JM IMAX = IM IF((J.EQ.1).OR.(J.EQ.JM)) IMAX = 1 DO 360 I=1,IMAX DO 360 K=1,NPL ACNP(I,J,K,N)=ACNP(I,J,K,N)+CNP(I,J,K,N)*BETA(I,J,K) C** products UC, VC, WC at pressure levels c** check grids! AUCNP(I,J,K,N)=AUCNP(I,J,K,N)+UP(I,J,K)*CNP(I,J,K,N)*BETA(I,J,K) AVCNP(I,J,K,N)=AVCNP(I,J,K,N)+VP(I,J,K)*CNP(I,J,K,N)*BETA(I,J,K) AWCNP(I,J,K,N)=AWCNP(I,J,K,N)+WP(I,J,K)*CNP(I,J,K,N)*BETA(I,J,K) C** accumulate Mole Fraction of SF6 (pptv) (global mean not subtracted) ACHIP(I,J,K,N)=ACHIP(I,J,K,N)+CIP(I,J,K,N)*UNIT*BETA(I,J,K) 360 CONTINUE 380 CONTINUE 400 CONTINUE C C.................................................................... C** save sf6 mixing ratio and meteorological info at station locations c** for last year of run C write(6,*)' tau jyear iyeare ',tau,jyear,iyeare IF(JYEAR.LT.IYEARE)RETURN c** position dataset before writing taus=0.0 do 410 ks=1,ndt read(kstn,end=480)taus 410 continue WRITE (6,*) ' DATA SET POSITIONED BEFORE WRITING STATION DATA ' WRITE (6,*) ' NO OF RECORDS SO FAR ',NDT WRITE (6,*) ' LAST TAU IN DATASET, CURRENT TAU ',TAUS,TAU C**** c** get geopotential height (not implemeted yet) c call getgph c NDT=NDT+1 write(6,*)' ndt = ',ndt TIME=TAU DO 450 NS=1,NSTN PSTN(NS)=PSAVE(IS(NS),JS(NS)) DO 440 L=1,LM c** also need height above ground (m) of each model layer c HSTN(L,NS)=?? c TSTN(L,NS)=T(IS(NS),JS(NS),L) c QSTN(L,NS)=Q(IS(NS),JS(NS),L) VSTN(L,NS)=V(IS(NS),JS(NS),L) c ** boundary layer depth if we have it c BHSTN(L,NS)=??(IS(NS),JS(NS),L) c** mixing ratio DO 430 N=1,NTRACE XVSTN(L,N,NS)=RM(IS(NS),JS(NS),L,N)/M(IS(NS),JS(NS),L) 430 continue 440 continue 450 continue c** write current station info to disk write(kstn)tau,ndt * ,(PSTN(ns),ns=1,nstn) * ,(((XVSTN(l,n,ns),l=1,lm),n=1,ntrace),ns=1,nstn) * ,((HSTN(l,ns),l=1,lm),ns=1,nstn) * ,((TSTN(l,ns),l=1,lm),ns=1,nstn) * ,((QSTN(l,ns),l=1,lm),ns=1,nstn) * ,((VSTN(l,ns),l=1,lm),ns=1,nstn) * ,((BHSTN(l,ns),l=1,lm),ns=1,nstn) * ,tau rewind kstn RETURN 480 WRITE(6,*)' UNABLE TO POSITION DISK BEFORE WRITING STATION DATA' write(6,*)' KSTN,TAU,TAUS=',KSTN,TAU,TAUS STOP 480 RETURN C--------------------------------------------------------------- ENTRY TRANSAC C C** accumulate change of tracer mass due to convection C** add line 2033.5 RMCNV(i,j,l,n)=rmij(l)-rm(i,j,l,n) C** call in subroutine CONVEC (at line 2038.5) c** this should be done every nconv hours, not ndiag hours C write(6,*)' in transac: idacc(4) ',idacc(4) DO 500 J=1,JM IMAX = IM IF((J.EQ.1).OR.(J.EQ.JM)) IMAX = 1 DO 490 K=1,NPL DO 490 I=1,IMAX ABETAC(I,J,K)=ABETAC(I,J,K)+BETA(I,J,K) 490 CONTINUE 500 CONTINUE DO 550 N=1,NTRACE DO 520 J=1,JM IMAX = IM IF((J.EQ.1).OR.(J.EQ.JM)) IMAX = 1 DO 510 L=1,LM DO 510 I=1,IMAX if(m(i,j,l).ne.0.0)RMCNV(I,J,L,N)=RMCNV(I,J,L,N)/M(I,J,L) 510 CONTINUE 520 CONTINUE C** interpolate to pressure levels CALL PNTRP(IM,JM,LM,NPL,SIGE,SIG,DSIG,PLEV,SKIP,PSAVE, * RMCNV(1,1,1,N),RMCNVP(1,1,1,N)) DO 540 J=1,JM IMAX = IM IF((J.EQ.1).OR.(J.EQ.JM)) IMAX = 1 DO 530 K=1,NPL DO 530 I=1,IMAX ARMCNVP(I,J,K,N)=ARMCNVP(I,J,K,N)+RMCNVP(I,J,K,N)*BETA(I,J,K) 530 CONTINUE 540 CONTINUE 550 CONTINUE RETURN C--------------------------------------------------------------- ENTRY TRANSAD C C** accumulate change of tracer mass due to diffusion C** add line 1780.5 RMDIF(i,j,l,n)=(FX(J)-FX(J+1)) (NS only) C** add line 1785.3 RMDIF(1,1,l,n)=RM8S C** add line 1785.4 RMDIF(1,jm,l,n)=RM8N C** call in subroutine DIFFUS (at line 1788.5) c** this should be done every ndfus hours, not ndiag hours C IDACC(3) = IDACC(3)+1 write(6,*)' in transad: idacc(3) ',idacc(3) DO 555 J=1,JM IMAX = IM IF((J.EQ.1).OR.(J.EQ.JM)) IMAX = 1 DO 554 K=1,NPL DO 554 I=1,IMAX ABETAD(I,J,K)=ABETAD(I,J,K)+BETA(I,J,K) 554 CONTINUE 555 CONTINUE DO 600 N=1,NTRACE DO 570 J=1,JM IMAX = IM IF((J.EQ.1).OR.(J.EQ.JM)) IMAX = 1 DO 560 L=1,LM DO 560 I=1,IMAX if(m(i,j,l).ne.0.0)RMDIF(I,J,L,N)=RMDIF(I,J,L,N)/M(I,J,L) 560 CONTINUE 570 CONTINUE C** interpolate to pressure levels CALL PNTRP(IM,JM,LM,NPL,SIGE,SIG,DSIG,PLEV,SKIP,PSAVE, * RMDIF(1,1,1,N),RMDIFP(1,1,1,N)) DO 590 J=1,JM IMAX = IM IF((J.EQ.1).OR.(J.EQ.JM)) IMAX = 1 DO 580 K=1,NPL DO 580 I=1,IMAX ARMDIFP(I,J,K,N)=ARMDIFP(I,J,K,N)+RMDIFP(I,J,K,N)*BETA(I,J,K) 580 CONTINUE 590 CONTINUE 600 CONTINUE C RETURN END C==================================================================== SUBROUTINE PNTRP(IM,JM,LM,NPL,SIGE,SIG,DSIG,PLEV,SKIP,PSURF, * ARRAY,ARRAYP) c c** subroutine to interpolate an array to pressure levels DIMENSION SIGE(LM+1),SIG(LM),DSIG(LM),PLEV(NPL) DIMENSION PSURF(IM,JM),ARRAY(IM,JM,LM),ARRAYP(IM,JM,NPL) c c write(6,*)' interpolating array to pressure levels' c write(6,*)' original array at I=36,J=23' c write(6,*)(array(36,23,l),l=1,lm) c** initialize output array DO 60 K=1,NPL DO 60 J=1,JM IMAX = IM IF((J.EQ.1).OR.(J.EQ.JM)) IMAX = 1 DO 50 I=1,IMAX ARRAYP(I,J,K)=SKIP 50 CONTINUE 60 CONTINUE C C*** INTERPOLATE TO CONSTANT PRESSURE SURFACES C PTOP=10. SP=1000. DO 300 J=1,JM IMAX=IM IF(J.EQ.1.OR.J.EQ.JM)IMAX=1 DO 200 I=1,IMAX PSF=PSURF(I,J) P1 =SIGE(1)*(PSF-PTOP)+PTOP PLM=SIGE(LM+1)*(PSF-PTOP)+PTOP LBOT=0 LTOP=0 C c** determine start/end layers and distance for interpolation DO 190 LP=1,NPL K=LP PINT=PLEV(LP) c** skip out if pressure level > surface pressure IF(PINT.GT.PSF) GO TO 185 IF(PINT.LT.PLM) GO TO 200 c IF(PINT.LT.PTOP) GO TO 200 SIGP=(PINT-PTOP)/(PSF-PTOP) DO 170 L=2,LM LTOP=L IF(SIGP.GE.SIG(L)) GO TO 180 170 CONTINUE 180 LBOT=LTOP-1 DIST=(SIGP-SIG(LTOP))/DSIG(LBOT) ARRAYP(I,J,K)=ARRAY(I,J,LBOT)*DIST+ARRAY(I,J,LTOP)*(1.-DIST) GO TO 190 c 185 CONTINUE c** Note: there may be missing layers in interpolated array c** the terrain mask will take care of these points in the accumulation 190 CONTINUE 200 CONTINUE 300 CONTINUE c write(6,*)' interpolated array at I=36,J=23' c write(6,*)(arrayp(36,23,k),k=1,npl) RETURN END C==================================================================== SUBROUTINE STNIJ(NGRID) C**** THIS PROGRAM FINDS THE GRIDBOX OF THE GMCC MONITORING STATIONS PARAMETER (NSTN=20) CHARACTER STN*4,AGRID*8 COMMON /STATION/ STN(NSTN),BLAT(NSTN),BLON(NSTN),HTM(NSTN) LOGICAL DEBUG/.TRUE./ COMMON/STIJ/IS(NSTN),JS(NSTN),LS(NSTN) DATA IFIRST/1/ C C**** FIND LEVEL IN VERTICAL IF(IFIRST.NE.1) GO TO 70 IFIRST=2 DO 50 NS = 1, NSTN LS(NS) = 1 IF (HTM(NS).GT.1000. .AND. HTM(NS).LT.3000.) LS(NS) = 2 IF (HTM(NS).GE.3000. .AND. HTM(NS).LT.4000.) LS(NS) = 3 IF (HTM(NS).GE.4000.) LS(NS) = 4 50 CONTINUE WRITE(6,*) LS C 70 CONTINUE GOTO (71,72,73),NGRID 71 AGRID = 'MEDIUM' IM = 36 JM = 24 DIVJ=JM-1. OFFJ=-0.5 OFFI=-0.5 GO TO 75 72 AGRID = 'FINE' IM=72 JM=46 DIVJ=JM-1. OFFJ=-0.5 OFFI=-0.5 GO TO 75 73 AGRID= 'N GRID' IM=36 JM=24 DIVJ=22.5 OFFJ=-0.75 OFFI=-0.25 75 CONTINUE DLON = 360./IM DLAT = 180./DIVJ ELAT = -90.+ OFFJ*DLAT ELON = OFFI*DLON IGREN = IM/2 + 1 C DO 100 N = 1, NSTN LAT = BLAT(N) AMIN = 100.*(BLAT(N)-LAT) FLAT = LAT + AMIN/60. JS(N) = (FLAT - ELAT)/DLAT + 1 LON = BLON(N) AMIN = 100.*(BLON(N) - LON) FLON = LON + AMIN/60. IS(N) = (FLON - ELON)/DLON + IGREN 100 CONTINUE 300 CONTINUE IF (.NOT.DEBUG) RETURN WRITE(6,60) AGRID DO 400 N = 1, NSTN RLAT=BLAT(N)*3.14159/180. SLAT=SIN(RLAT) WRITE(6,61) STN(N),BLAT(N),SLAT,BLON(N),HTM(N),IS(N),JS(N),LS(N) 400 CONTINUE C 60 FORMAT(1X,'STN',7X,'LAT',5X,'SIN(LAT)',3X,'LON',10X,A8,/) 61 FORMAT(1X,A4,2X,4F8.2,10X,3I3,10X,I3) RETURN END C==================================================================== BLOCK DATA c** station locations for high-frequency diagnostics C NOTE THAT DECIMALS IN BLAT AND BLON ARE MINUTES c PARAMETER (NSTN=20) CHARACTER STN*4 COMMON /STATION/ STN(NSTN),BLAT(NSTN),BLON(NSTN),HTM(NSTN) DATA STN/ 1 'SPO ','NEU ','TDF ','CGO ','SMO ', 2 'RPB ','GMI ','MLO ','IZO ','BME ', 3 'WITN','QPC ','UTA ','WLEF','HUN ', 4 'FRA ','MHT ','CBA ','SVA ','ALT '/ c c** latitude in degrees and minutes eg smo: 14 degrees 14 mins DATA BLAT/ 1 -90.00,-71.00,-54.52,-41.00,-14.14, 2 13.10, 13.26, 19.32, 28.00, 32.22, 3 35.37, 36.16, 39.54, 45.56, 46.58, 4 50.00, 53.20, 55.12, 78.55, 82.27/ c c** longitude in degrees and minutes DATA BLON/ 1 -102.00, -8.00, -68.29, 145.00, 170.34, 2 -59.26, 144.47,-155.35, -16.00, -64.39, 3 -77.39, 100.55,-113.43, -90.16, 16.23, 4 -82.00, -9.54,-162.43, -11.53, -62.31/ c c** elevation in meters DATA HTM/ 1 2835., 42., 20., 95., 77., 2 3., 2.,3397.,2367., 30., 3 496.,3810.,1320., 868., 240., 4 200., 26., 25., 474., 210./ END