SUBROUTINE VSH(OUTPUT, USING, GLON, GLAT, ALT, UT, JDAY, 1 F107, AP) C++++++++++++++++++++++++ VSH VERSION 2.0 +++++++++++++++++++++++++++++++++ C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C TIM KILLEEN UNIVERSITY OF MICHIGAN AND NCAR C ROB RASKIN UNIVERSITY OF MICHIGAN C ****************** GENERAL DESCRIPTION: ************************ C The geophysical fields are based on output from the NCAR-TIGCM. C The routine uses a stored set of coefficients for each of C several NCAR-TGCM runs. The coefficients contain a C representation of the global thermospheric wind, temperature C and density fields using a set of Fourier C coefficients in time, spline coefficients in altitude C and vector spherical harmonics in space. C A linear interpolation is carried out over the coefficient space C for the specified Ap, F10.7, and Julian day. MSIS is called for H C densities, and for altitudes too high or too low. C Modifications to the mass density are made for some of the TIGCM C runs, based on MSIS global means. These correction C factors are read in from the file NORMALIZ.DAT C The ion drifts are based on the Heelis model used as input C to the TIGCM and contains the effects of collisions at altitudes C where the ion collision frequency is significant. C ***************** OUTPUT VARIABLES: *************************** C Name Type Description C ---- ---- ----------- C OUTPUT(25) REAL the array of output variables: C 1) Zonal neutral wind (U) (m/s) C (positive eastward) C 2) Meridional neutral wind (V) (m/s) C (positive northward) C 3) Neutral temperature (K) C 4) Zonal ion drift (Ui) (m/s) C 5) Meridional ion drift (Vi) (m/s) C 6) Pressure surfaace Z (5X10-4 * EXP(-Z) microbars) C 7) Vertical neutral wind (m/s) C 8) Atomic oxygen mass mixing ratio C 9) Mass density (g/cm3) C 10) Electron density (#/cm3) C 11) O+ density (#/cm3) C 12) Ion temperature (K) C 13) Molecular oxygen mass mixing ratio C 20) Atomic oxygen density (#/cm3) C 21) Molecular oxygen density (#/cm3) C 22) Molecular nitrogen density (#/cm3) C 23) Hydrogen density (#/cm3) C 24) Molecular nitrogen mass mixing ratio C (VARIABLES 14-19 AND 25 ARE RESERVED FOR FUTURE USE) C C ******************** INPUT VARIABLES: *************************** C Name Type Description C ---- ------ ------------- C USING(25) LOGICAL an array of 25 logical variables (.true. or .false.) C corresponding to each of the above 25 variables. C Set a value to .false. to eliminate calculation of C that variable (to reduce computational time). C Set all desired corresponding values to .true. C GLON REAL geographic longitude GLON (degrees) C GLAT REAL geographic latitude GLAT (degrees) C ALT REAL altitude (km) (between 110 and 1500 km) C UT REAL Universal Time UT (hours) C JDAY INTEGER Julian Day (between 1 and 365, inclusive) C F107 REAL F10.7 value of solar flux C AP REAL 3-hr ap index of geomagnetic activity C INPUT OPTIONS C Three input options are temporarily available for validation purposes: C 1. In lieu of specifying geophysical parameters (JDAY, F107, and AP), a C specific TIGCM run can be selected by setting JDAY= -run#, where run# C is a number 1 - 38. See documentation for the run descriptions. With this C option, input values of F107 and AP are not used. On output, JDAY, C F107, and AP retur the actual values for that TIGCM run. C 2. In lieu of specifying an altitude, a constant pressure surface C can be selected by setting ALT to a pressure surface (-7 to +5). C If this option is used, OUTPUT(6) returns the altitude. C C 3. In lieu of specifying an altitude, a natural log of total mass density C can be specified. This option will be used for real-time density updates. C For this option, OUTPUT(6) returns the TIGCM altitude of the specified C density and OUTPUT(9) returns the pressure. C *************** Parameter Declaratons ********************* INTEGER TMAX, VARMAX, FLDMAX, RUNMAX, ZMAX PARAMETER (MMAX=5, NMAX=25, TMAX=7, ZMAX=7, 1 FLDMAX=13, VARMAX=25, RUNMAX=60, LOWALT= 110) C ********************** Variable Declarations ********************* C VSHCOF contains Fourier coefficients as read in from the TIGCM history files C VSH2 is a temporary work space for the coefficients C VSHT contains the coefficients following Fourier time series transformation C VSHTZ contains the coefficients following time and vertical transformation C Prefix FLD represents the number of variables exclusive of densities C Prefix VAR represents the number of variables including densities DIMENSION VSHCOF (FLDMAX, ZMAX, MMAX, NMAX, TMAX, 2) DIMENSION VSH2 (FLDMAX, ZMAX, MMAX, NMAX, TMAX, 2) DIMENSION VSHT (MMAX, NMAX, 2, ZMAX, FLDMAX) DIMENSION VSHTZ (MMAX, NMAX, 2, FLDMAX), OUTPUT(VARMAX) LOGICAL USING(VARMAX), NEWTIME, NEWLAT, NEWREAD, NEWDATA, GETDENS LOGICAL NEWLONG, ALTSPEC, NEWALT, MSIS, NEWRUN, VECTORS, DENSPEC INTEGER ZTRUNC(FLDMAX), MTRUNC(FLDMAX), NTRUNC(FLDMAX), ZLOW DIMENSION PBAR(MMAX+1,NMAX+1), VBAR(MMAX+1, NMAX+1) DIMENSION WBAR(MMAX+1,NMAX+1), RUNWT(RUNMAX), DENSMOD(2) DIMENSION COSPHI(MMAX-1), SINPHI(MMAX-1), DENSFAC(2), TEMPFAC(2) DATA IPRINT/0/, NEWDATA/.FALSE./, DENSMOD/0,0/, DENSFAC/2*1/ C Validate input - What input values have changed since last call? CALL VALIDATE (USING, VARMAX, GETDENS, ALT, F107,JDAY, 1 AP, ALTSPEC, MSIS, VECTORS, DENSPEC, LAYER, LOWALT) IF (LAYER .GT. 1) THEN CALL WHATSNEW (USING, VARMAX, GLON, GLAT, UT, ALT, ALTSPEC, 1 JDAY,F107, AP, NEWLAT, NEWTIME, NEWRUN, NEWLONG, NEWALT,NEWREAD) C Determine which coefficient libraries to read IF (NEWRUN) CALL COEFFSEL (F107,JDAY,AP,RUNWT,RUNMAX) C Read spectral coefficients IF (NEWREAD) CALL COEFFIN(VSHCOF, USING, RUNWT, RUNMAX, DENSMOD, 1 FLDMAX, VARMAX, MMAX, NMAX, ZMAX, TMAX, 2 DENSFAC, MTRUNC, NTRUNC, ZTRUNC, VSH2, TEMPFAC, LOWALT,ZLOW) C Get coefficient modifiers if this is a new run or there is new data C IF (NEWREAD .OR. NEWDATA) C 1 CALL COEFFMOD (VSHCOF, FLDMAX, ZMAX, MMAX, NMAX, TMAX, C 2 USING, VARMAX) C Carry out temporal synthesis IF(NEWTIME) 1 CALL TIMETR(UT, VSHCOF, VSHT, USING, FLDMAX, VARMAX, 2 MMAX, NMAX, ZMAX, TMAX, MTRUNC, NTRUNC, ZTRUNC, ZLOW) C Get basis functions for spatial synthesis IF (NEWLAT .OR. NEWLONG) 1 CALL BASIS (PBAR, WBAR, VBAR, MMAX, NMAX, GLAT,GLON, 2 SINPHI, COSPHI, NEWLAT, NEWLONG, VECTORS) C Carry out altitude synthesis IF (NEWALT .OR. NEWLONG .OR. NEWLAT) # CALL ALTTR(VSHT,ZMAX,MMAX,NMAX,FLDMAX,MTRUNC, 1 NTRUNC,ZTRUNC,VSHTZ, ALTSPEC, COSPHI, SINPHI,PBAR, 2 PRESS,ALT,USING,VARMAX,DENSMOD, OUTPUT, 3 ALTTOP, LAYER, DENSPEC, ZLOW, TEMPMAX) C Carry out spatial synthesis IF (NEWLAT .OR. NEWLONG .OR. NEWALT) 1 CALL SPACTR (VSHTZ,MMAX,NMAX,FLDMAX,COSPHI,SINPHI,PBAR, 2 VBAR, WBAR, OUTPUT, VARMAX, MTRUNC, NTRUNC, USING, ALTSPEC) C Compute DERIVED FIELDS IF (USING(3)) OUTPUT(3)= TEMPMAX - EXP(OUTPUT(3)) IF (GETDENS .AND. LAYER .LE. 3) 1 CALL DENSITY (OUTPUT, USING, VARMAX, PRESS, MSIS, DENSPEC,ALT) IF (USING(5) .OR. USING(10) .OR. USING(11)) 1 CALL IONELEC (OUTPUT,VARMAX,GLAT, PRESS, LAYER, USING) ENDIF C Normalize to MSIS mean IF ( (LAYER .EQ. 2 .OR. LAYER .EQ. 3) .AND. 1 (USING(3) .OR. USING(9)) ) CALL NORMLZE 1 (GLAT, LAYER, USING, OUTPUT, VARMAX, TEMPFAC, DENSFAC) C Get parameters of run if run# was specified IF (JDAY .LT. 0) CALL RUNPARAM(JDAY, AP, F107) C Too high, too low, and for H density call MSIS. IF ( (MSIS .OR. LAYER .NE. 2) .AND. .NOT. DENSPEC) CALL GETMSIS 1 (OUTPUT, USING, VARMAX, F107, JDAY, AP, UT, ALT, GLAT, GLON, 2 ALTTOP, LAYER) RETURN END SUBROUTINE VALIDATE(USING, VARMAX, GETDENS, ALT, F107, JDAY, AP, 1 ALTSPEC, MSIS, VECTORS, DENSPEC, LAYER, LOWALT) INTEGER VARMAX LOGICAL USING(VARMAX), GETDENS, ALTSPEC, MSIS, VECTORS, DENSPEC DENSPEC= .FALSE. ALTSPEC= .TRUE. MSIS= .FALSE. LAYER= 2 C CHECK FOR NEED FOR ALTITUDES FOR INTERPOLATION AND MSIS CALL IF (ALT .GE. -5.51 .AND. ALT .LE. 4.51) THEN C Pressure has been specified ALTSPEC= .FALSE. ELSEIF (ALT .GE. LOWALT .AND. ALT .LE. 800) THEN C Altitude is in range USING(6)= .TRUE. ELSEIF (ALT .GE. 800 .AND. ALT .LE. 1500.001) THEN C Altitude is too high USING(6)= .TRUE. MSIS= .TRUE. LAYER= 4 ELSEIF (ALT .GT. 90 .AND. ALT .LT. LOWALT-.001) THEN C Altitude is too low MSIS= .TRUE. LAYER= 1 ELSEIF (ALT .GE. -41 .AND. ALT .LE. -18) THEN C Log density has been specified ALTSPEC= .FALSE. DENSPEC= .TRUE. ELSE WRITE (*,*) 'ALTITUDE: ',ALT,' OUT OF RANGE' STOP ENDIF IF (USING(23)) THEN C Hydsrogen needs MSIS MSIS= .TRUE. USING(6) = .TRUE. ENDIF C Check for consistent specification of USING C Winds need each other; total density needs constituent densities C Densities need mixing ratios and temperatures; interpolation needs heights VECTORS= .FALSE. IF (USING(1)) USING(2)= .TRUE. IF (USING(2)) THEN USING(1)= .TRUE. VECTORS= .TRUE. ENDIF IF (USING(4)) USING(5)= .TRUE. IF (USING(5)) THEN USING(4)= .TRUE. USING(3)= .TRUE. IF (ALTSPEC) USING(6)= .TRUE. VECTORS= .TRUE. ENDIF IF (USING(22)) USING(24)= .TRUE. IF (USING(20) .OR. USING(24)) USING(8) =.TRUE. IF (USING(21) .OR. USING(24)) USING(13)=.TRUE. IF (USING(20) .OR. USING(21) .OR. USING(22)) USING(9) =.TRUE. IF (USING(9) .OR. USING(20) .OR. USING(21) .OR. USING(24)) 1 GETDENS= .TRUE. C Other input checks to be put here RETURN END SUBROUTINE WHATSNEW (USING, VARMAX, GLON, GLAT, UT, ALT, ALTSPEC, 1 JDAY,F107, AP, NEWLAT, NEWTIME, NEWRUN, NEWLONG, NEWALT,NEWREAD) C Determine what requested values have changed from last call INTEGER VARMAX, VAR LOGICAL USING(VARMAX), LASTUSI(25), NEWLAT,NEWTIME,NEWRUN LOGICAL NEWLONG, NEWALT, NEWREAD, ALTSPEC REAL LASTLAT, LASTUT, LASTAP, LASTF10, LASTLON, LASTALT DATA LASTLAT/-9999/, LASTUT/-9999/, LASTALT/-999/ DATA LASTF10/-999/, LASTAP/-999/, LASTDAY/-999/, LASTLON/-9999/ DATA LASTUSI/25*.FALSE./ SAVE LASTLAT, LASTUT, LASTAP, LASTF10, LASTLON, LASTALT SAVE LASTUSI, LASTDAY C Is this is a new run? IF (ABS(AP - LASTAP) .GT. 1.5 .OR. ABS(F107 - LASTF10) .GT. 5 1 .OR. ABS(JDAY - LASTDAY) .GT. 2.5) THEN NEWRUN= .TRUE. LASTF10= F107 LASTDAY= JDAY LASTAP = AP ELSE NEWRUN= .FALSE. ENDIF C Determine if there is a change in the variables requested NEWREAD= .FALSE. DO 44 VAR= 1,VARMAX IF (USING(VAR) .NE. LASTUSI(VAR)) THEN NEWREAD= .TRUE. LASTUSI(VAR)= USING(VAR) ENDIF 44 CONTINUE IF (NEWRUN) NEWREAD= .TRUE. C Is this is a new time? IF (ABS(UT - LASTUT) .GT. 0.3 .OR. NEWREAD) THEN NEWTIME= .TRUE. LASTUT= UT ELSE NEWTIME= .FALSE. ENDIF C Is this is a new latitude? IF (ABS(GLAT - LASTLAT) .GT. 1) THEN NEWLAT = .TRUE. LASTLAT= GLAT ELSE NEWLAT= .FALSE. ENDIF C Is this a new alt? IF (ABS(ALT-LASTALT) .GT. 1E-3 .OR. NEWTIME) THEN NEWALT= .TRUE. LASTALT= ALT ELSE NEWALT= .FALSE. ENDIF C Is this a new longitude? IF (ABS(GLON - LASTLON) .GT. 1.5) THEN LASTLON= GLON NEWLONG= .TRUE. ELSE NEWLONG= .FALSE. ENDIF C Does altitude need to be recomputed from pressure? IF ((NEWLAT .OR. NEWLON) .AND. ALTSPEC) NEWALT= .TRUE. RETURN END SUBROUTINE COEFFSEL (F10PT7,JDAYREQ,AP1,RUNWT,RUNMAX) C Coefficient Selector C Finds the weights that are used to construct the weighted average of C TIGCM runs for given geophysical conditions C Inputs: F10.7 Flux, Julian day, AP INTEGER F107MAX, RUNMAX, RUN, APMAX PARAMETER (APMAX= 3, JDAYMAX= 3, F107MAX=2, HALFPI=1.570795) REAL F107WT(F107MAX), APWT(APMAX), JDAYWT(JDAYMAX) DIMENSION F107(F107MAX), AP(APMAX), JDAY(JDAYMAX+1) DIMENSION RUNWT(RUNMAX) DATA AP/5,11,32/, F107/67,243/, JDAY/80,182,263,355/ DATA APINC/.4/, F107INC/3.5/, JDAYINC/5/ C Initialize DO 20 J= 1,RUNMAX RUNWT(J)= 0 20 CONTINUE DO 30 J= 1,JDAYMAX JDAYWT(J)= 0 30 CONTINUE DO 40 K= 1,F107MAX F107WT(K)= 0 40 CONTINUE DO 50 L= 1,APMAX APWT(L)= 0 50 CONTINUE C JULIAN DAY IF (JDAYREQ .LT. 0) THEN RUNWT(ABS(JDAYREQ))= 1. RETURN ELSEIF (JDAYREQ .GT. 365) THEN JDAYREQ= MOD(365,JDAYREQ) ENDIF c Within threshold of a TIGCM run? IF (ABS(JDAYREQ-JDAY(1)) .LT. JDAYINC) THEN JDAYWT(1)= 1 ELSEIF (ABS(JDAYREQ-JDAY(2)) .LT. JDAYINC) THEN JDAYWT(2)= 1 ELSEIF (ABS(JDAYREQ-JDAY(3)) .LT. JDAYINC) THEN JDAYWT(1)= 1 ELSEIF (ABS(JDAYREQ-JDAY(4)) .LT. JDAYINC) THEN JDAYWT(3)= 1 ELSEIF (JDAYREQ .LE. JDAY(1)) THEN C Winter ALPHA= ( JDAYREQ - JDAY(4) + 365 ) / 1 FLOAT( JDAY(1) - JDAY(4) + 365 ) JDAYWT(3)= SIN (HALFPI * (1 + ALPHA)) JDAYWT(1)= 1 - JDAYWT(3) ELSEIF (JDAYREQ .LE. JDAY(2)) THEN C Spring ALPHA= ( JDAYREQ - JDAY(1) ) / 1 FLOAT( JDAY(2) - JDAY(1) ) JDAYWT(2)= SIN (HALFPI * ALPHA) JDAYWT(1)= 1 - JDAYWT(2) ELSEIF (JDAYREQ .LE. JDAY(3)) THEN C Summer ALPHA= ( JDAYREQ - JDAY(2) ) / 1 FLOAT( JDAY(3) - JDAY(2) ) JDAYWT(2)= SIN (HALFPI * (1 + ALPHA)) JDAYWT(1)= 1 - JDAYWT(2) ELSEIF (JDAYREQ .LE. JDAY(4)) THEN C Fall ALPHA= ( JDAYREQ - JDAY(3) ) / 1 FLOAT( JDAY(4) - JDAY(3) ) JDAYWT(3)= SIN (HALFPI * ALPHA) JDAYWT(1)= 1 - JDAYWT(3) ELSE C Winter ALPHA= ( JDAYREQ - JDAY(4) ) / 1 FLOAT( 365 + JDAY(1) - JDAY(4) ) JDAYWT(3)= SIN (HALFPI * (1 + ALPHA)) JDAYWT(1)= 1 - JDAYWT(3) ENDIF C F10.7 SOLAR FLUX IF (F10PT7 .LT. 50) THEN F107REQ = 50 ELSEIF (F10PT7 .GT. 500) THEN F107REQ= 500 ELSE F107REQ= F10PT7 ENDIF c See if within threshold of a TIGCM run value of F10.7 IF (ABS(F107REQ-F107(1)) .LT. F107INC) THEN F107WT(2)= 0 ELSEIF (ABS(F107REQ-F107(2)) .LT. F107INC) THEN F107WT(2)= 1 ELSE c Take weighted average. Extrapolate if needed F107WT(2)= ( F107REQ - F107(1) ) / ( F107(2) - F107(1) ) ENDIF F107WT(1)= 1 - F107WT(2) C AP IF (AP1 .LE. 2) THEN APREQ= 2 ELSEIF (AP1 .GE. 60) THEN APREQ= 60 ELSE APREQ= AP1 ENDIF c See if within threshold of a TIGCM run value of AP IF (ABS(APREQ-AP(1)) .LT. APINC) THEN APWT(1)= 1 ELSEIF (ABS(APREQ-AP(2)) .LT. APINC) THEN APWT(2)= 1 ELSEIF (ABS(APREQ-AP(3)) .LT. APINC) THEN APWT(3)= 1 c Below formulas allow for extrapolation ELSEIF (APREQ .LT. AP(2)) THEN APWT(2)= ( APREQ - AP(1) ) / ( AP(2) - AP(1) ) APWT(1)= 1 - APWT(2) ELSE APWT(3)= ( APREQ - AP(2) ) / ( AP(3) - AP(2) ) APWT(2)= 1 - APWT(3) ENDIF DO 500 J= 1,JDAYMAX DO 400 K= 1,F107MAX DO 300 L= 1,APMAX C Convert to RUN # RUN= 11 + L + 6*J + 3*K RUNWT(RUN)= JDAYWT(J) * F107WT(K) * APWT(L) 300 CONTINUE 400 CONTINUE 500 CONTINUE RETURN END SUBROUTINE COEFFIN(VSHCOF, USING, RUNWT, RUNMAX, 1 DENSMOD, FLDMAX, VARMAX, MMAX, NMAX, ZMAX, TMAX, 2 DENSFAC, MTRUNC, NTRUNC, ZTRUNC, VSH2, TEMPFAC,LOWALT,ZLOW) C Read in the TIGCM spectral coefficients INTEGER ZMAX, Z, TMAX, T, FLDMAX INTEGER ZTRUNC(FLDMAX), NTRUNC(FLDMAX), TTRUNC, PARITY INTEGER VAR, VARMAX, RUN, RUNMAX, RUNMOD INTEGER MTRUNC(FLDMAX), ZLOW DIMENSION DEN1FAC(60), DEN2FAC(60), DENSFAC(2), TEMPFAC(2) DIMENSION RUNWT(RUNMAX), DENSMOD(2), DENSADJ(60,2), TEMPADJ(60,2) DIMENSION VSHCOF (FLDMAX, ZMAX, MMAX, NMAX, TMAX, 2) DIMENSION VSH2 (FLDMAX, ZMAX, MMAX, NMAX, TMAX, 2) PARAMETER (SQRT2=1.414) LOGICAL USING(VARMAX), FIRST CHARACTER*2 CHARRUN CHARACTER*35 FINAME DATA DENSADJ /120*1./, TEMPADJ/120*1/, FIRST /.TRUE./ DATA DEN1FAC /60*1/, DEN2FAC/60*1/ SAVE FIRST, DENSADJ, TEMPADJ C Get density modifiers for Northern & Southern Northern Hemispheres IF (FIRST .AND. (USING(9) .OR. USING(3))) THEN OPEN (3, STATUS='OLD', . FILE='/u/tgcm/vsh13/normaliz.dat') DO 95 RUN= 1,RUNMAX READ (3,950) RUNMOD, DENSADJ(RUN,1), TEMPADJ(RUN,1), . DENSADJ(RUN,2), TEMPADJ(RUN,2) 95 CONTINUE CLOSE (3) FIRST= .FALSE. ENDIF C Initialize DO 100 VAR= 1,FLDMAX DO 90 Z= 1,ZMAX DO 80 M= 1,MMAX DO 70 T= 1,TMAX DO 60 N= M,NMAX DO 50 PARITY= 1,2 VSHCOF(VAR,Z,M,N,T,PARITY)= 0. 50 CONTINUE 60 CONTINUE 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE DO 110 I= 1,2 DENSFAC(I)= 0 TEMPFAC(I)= 0 110 CONTINUE DO 2000 RUN= 1,RUNMAX WEIGHT= RUNWT(RUN) IF (ABS(WEIGHT) .GT. 0.001) THEN C Calculate weighted average of density and temperature modifiers IF (USING(9)) THEN DENSFAC(1)= DENSFAC(1) + WEIGHT * DENSADJ(RUN,1) DENSFAC(2)= DENSFAC(2) + WEIGHT * DENSADJ(RUN,2) c DENSMOD(1)= DENSMOD(1) + WEIGHT * DEN1FAC(RUN) c DENSMOD(2)= DENSMOD(2) + WEIGHT * DEN2FAC(RUN) ENDIF IF (USING(3)) THEN TEMPFAC(1)= TEMPFAC(1) + WEIGHT * TEMPADJ(RUN,1) TEMPFAC(2)= TEMPFAC(2) + WEIGHT * TEMPADJ(RUN,2) ENDIF C Convert run# to a character string CALL RUNTOCH( RUN, CHARRUN ) FINAME='/u/tgcm/vsh13/run0' //CHARRUN// '.new' OPEN (2, STATUS='OLD', FILE=FINAME, 1 FORM='UNFORMATTED') READ (2) ZLOW DO 400 VAR= 1,FLDMAX IF (USING(VAR)) THEN READ(2) IVAR,MTRUNC(VAR), NTRUNC(VAR), 1 ZTRUNC(VAR), TTRUNC READ (2) (((((VSH2(VAR,Z,M,N,T,PARITY), 1 N=M,NTRUNC(VAR)), M= PARITY, MTRUNC(VAR)), 2 PARITY= 1,2), T= 1,TTRUNC), Z=ZLOW,ZTRUNC(VAR)) C Take a weighted average of the coefficients DO 900 PARITY= 1,2 DO 800 M= PARITY,MTRUNC(VAR) DO 700 N= M,NTRUNC(VAR) DO 600 Z= ZLOW,ZTRUNC(VAR) DO 500 T= 1,TMAX IF (VAR .EQ. 3 .AND. Z .LE. ZTRUNC(3)-2) THEN C Temperatures need to be converted out of log form for this interpolation IF (N .EQ. 1 .AND. T .EQ. 1) THEN VSHCOF(3,Z,M,N,T,PARITY)= VSHCOF(3,Z,M,N,T,PARITY) 1 + WEIGHT * EXP(VSH2(3,Z,M,N,T,PARITY)/2/SQRT2) ELSE VSHCOF(3,Z,M,N,T,PARITY)= VSHCOF(3,Z,M,N,T,PARITY) 1 + WEIGHT * EXP(VSH2(3,Z,M,N,T,PARITY)) ENDIF ELSE VSHCOF(VAR,Z,M,N,T,PARITY)= VSHCOF(VAR,Z,M,N,T,PARITY) 1 + WEIGHT * VSH2(VAR,Z,M,N,T,PARITY) ENDIF 500 CONTINUE 600 CONTINUE 700 CONTINUE 800 CONTINUE 900 CONTINUE ELSE READ (2) READ (2) ENDIF 400 CONTINUE CLOSE(2) ENDIF 2000 CONTINUE IF (USING(3)) THEN DO 1900 PARITY= 1,2 DO 1800 M= PARITY,MTRUNC(3) DO 1700 N= M,NTRUNC(3) DO 1600 Z= ZLOW,ZTRUNC(3)-2 DO 1500 T= 1,TMAX IF (N .EQ. 1 .AND. T .EQ. 1) THEN VSHCOF(3,Z,M,N,T,PARITY)= 1 ALOG(VSHCOF(3,Z,M,N,T,PARITY))*2*SQRT2 ELSE VSHCOF(3,Z,M,N,T,PARITY)= 1 ALOG(VSHCOF(3,Z,M,N,T,PARITY)) ENDIF 1500 CONTINUE 1600 CONTINUE 1700 CONTINUE 1800 CONTINUE 1900 CONTINUE ENDIF 950 FORMAT (1X,I3,2E10.3,4X,2E10.3) RETURN END SUBROUTINE COEFFMOD (VSHCOF, FLDMAX, ZMAX, MMAX, NMAX, TMAX, 2 USING, VARMAX) C Will provide modification of coefficients based on satellite data C Currently only operational on the Cray INTEGER FLDMAX, ZMAX, TMAX, VARMAX LOGICAL USING(VARMAX) DIMENSION VSHCOF(FLDMAX, ZMAX, MMAX, NMAX, TMAX, 2) RETURN END SUBROUTINE TIMETR(UT, VSHCOF, VSHT, USING, FLDMAX, VARMAX, 1 MMAX, NMAX, ZMAX, TMAX, MTRUNC, NTRUNC, ZTRUNC, ZLOW) C CARRIES OUT A FOURIER SUMMATION IN TIME TO ELIMINATE EXPLICIT TIME C DEPENDENCE. THE ARRAY VSHCOF IS CONVERTED INTO VSHT INTEGER ZMAX, TMAX, VARMAX, Z, T, FLDMAX,HR1,VAR, ZLOW INTEGER ZTRUNC(FLDMAX), PARITY, MTRUNC(FLDMAX), NTRUNC(FLDMAX) DIMENSION VSHCOF(FLDMAX, ZMAX, MMAX, NMAX, TMAX, 2) DIMENSION VSHT (MMAX, NMAX, 2, ZMAX, FLDMAX) DIMENSION SINPHA(12), COSPHA(12) LOGICAL USING(VARMAX) PARAMETER (HR1=1) C CALCULATE PHASE ANGLE TWOPI = 8 * ATAN(1.) PHASE = TWOPI * (UT - HR1) / 24 C PRECOMPUTE SIN AND COSINE FUNCTIONS DO 10 K= 1, (TMAX-1)/2 SINPHA(K)= SIN(K*PHASE) COSPHA(K)= COS(K*PHASE) 10 CONTINUE DO 60 VAR= 1,FLDMAX IF (USING(VAR)) THEN DO 55 PARITY= 1,2 DO 50 M= PARITY,MTRUNC(VAR) DO 40 N= M,NTRUNC(VAR) DO 30 Z= ZLOW,ZTRUNC(VAR) SUM= 0 TMEAN= VSHCOF(VAR,Z,M,N,1,PARITY) DO 20 T= 1,(TMAX-1)/2 SUM = SUM + 1 VSHCOF(VAR,Z,M,N,2*T,PARITY) * COSPHA(T) 2 - VSHCOF(VAR,Z,M,N,2*T+1,PARITY) * SINPHA(T) 20 CONTINUE VSHT(M,N,PARITY,Z,VAR) = TMEAN + 2* SUM 30 CONTINUE 40 CONTINUE 50 CONTINUE 55 CONTINUE ENDIF 60 CONTINUE RETURN END SUBROUTINE BASIS(PBAR, WBAR, VBAR, MMAX, NMAX, GLAT, GLON, 1 SINPHI, COSPHI, NEWLAT, NEWLONG, VECTORS) C Determines the components of the vector spherical harmonics C as a function of order M and degree N C Inputs: MMAX, NMAX (truncation levels); GLAT (geographic latitude in degrees) C Outputs: PBAR, WBAR, VBAR (the three VSH coefficient matricies) DIMENSION PBAR(MMAX+1,NMAX+1), WBAR(MMAX+1,NMAX+1) DIMENSION VBAR(MMAX+1,NMAX+1), W(125), PB(26) DIMENSION COSPHI(MMAX-1), SINPHI(MMAX-1) LOGICAL NEWLAT, NEWLONG, VECTORS PI= 4 * ATAN(1.) IF (NEWLAT) THEN THETA= (90-GLAT) * PI/180 C CALCULATE NORMALISED ASSOCIATED LEGENDRE FUNCTIONS (PBAR): DO 15 MP1= 1,MMAX+1 CALL LFNC(0,MP1-1,NMAX,THETA,PB,W) CALL LFNC(1,MP1-1,NMAX,THETA,PB,W) DO 10 NP1= MP1,NMAX PBAR(MP1,NP1)= PB(NP1) 10 CONTINUE 15 CONTINUE IF (VECTORS) THEN C Calculate vector basis functions (VBAR and WBAR): C VBAR= dP/dlat WBAR= P * m / sin (colat) DO 40 MP1= 1,MMAX M= MP1 - 1 DO 35 NP1= MP1,NMAX N= NP1 - 1 IF(MP1.GT.1) THEN VBAR(MP1,NP1) = (SQRT((N+M)*(N-M+1.)) * PBAR(M,NP1) 1 - SQRT((N-M)*(N+M+1.)) * PBAR(M+2,NP1)) 2 / (2*SQRT((N*(N+1.)))) WBAR(MP1,NP1) = (SQRT((N+M)*(N+M-1.)) * PBAR(M,N) 1 + SQRT((N-M)*(N-M-1.)) * PBAR(M+2,N) ) 2 * SQRT((2*N+1.)/(2*N-1.)) /(2*SQRT(N*(N+1.))) ELSE VBAR(1,NP1)= -PBAR(2,NP1) WBAR(1,NP1)= 0 ENDIF 35 CONTINUE 40 CONTINUE ENDIF ENDIF IF (NEWLONG) THEN PHI = GLON * PI/180. + PI C PRECOMPUTE COS AND SIN FUNCTIONS DO 80 M= 1,MMAX-1 SINPHI(M)= SIN(M*PHI) COSPHI(M)= COS(M*PHI) 80 CONTINUE ENDIF RETURN END SUBROUTINE ALTTR(VSHT,ZMAX,MMAX,NMAX,FLDMAX,MTRUNC,NTRUNC,ZTRUNC, 1 VSHTZ, ALTSPEC, COSPHI, SINPHI,PBAR,PRESS,ALT,USING, 3 VARMAX,DENSMOD,OUTPUT,TOP,LAYER, DENSPEC, ZLOW, TEMPMAX) INTEGER FLDMAX, ZMAX, Z, VARMAX, PARITY, VAR, ZLOW INTEGER MTRUNC(FLDMAX), NTRUNC(FLDMAX), ZTRUNC(FLDMAX) PARAMETER (NCOEFF=7, KORDER=4) LOGICAL ALTSPEC, USING(VARMAX), DENSPEC DIMENSION COSPHI(MMAX-1), SINPHI(MMAX-1), HT(NCOEFF) DIMENSION PBAR(MMAX+1,NMAX+1), OUTPUT(VARMAX) DIMENSION VSHT (MMAX, NMAX, 2, ZMAX, FLDMAX) DIMENSION VSHTZ (MMAX, NMAX, 2, FLDMAX), BSPLINE(KORDER) DIMENSION DENSMOD(2), WORK(100), COEFF(NCOEFF) REAL KNOT(NCOEFF+KORDER), LASTP DATA LASTP /-999/ DATA KNOT /4*-5.51,-3.5,-1.5,1.5,4*4.51/, COEFF/7*0/, HT/7*0/ SAVE BSPLINE, LASTP, KNOTLO DO 77 I= 1,VARMAX OUTPUT(I)= 0. 77 CONTINUE LAYER= 2 IF (ALTSPEC .OR. DENSPEC) THEN C Altitude or density was specified in ALT IF (ALTSPEC) THEN VAR= 6 ELSE VAR= 9 ENDIF C An inversion of the spline is required to obtain pressure from altitude C Take spatial transform for each spline coefficient DO 100 Z= ZLOW,ZTRUNC(VAR) CALL SCATRAN (VSHT(1,1,1,Z,VAR), COSPHI, SINPHI, PBAR, 1 MMAX, NMAX, MTRUNC(VAR), NTRUNC(VAR), HT(Z)) 100 CONTINUE TOP= HT(ZTRUNC(VAR)) IF (ALTSPEC .AND. ALT .GT. TOP + 100) THEN C Altitude is above the 100 km blend layer. Use MSIS LAYER= 4 PRESS= KNOT(NCOEFF+KORDER) - .01 ELSEIF (ALTSPEC .AND. ALT .GT. TOP) THEN C Altitude is too high but within 100 km of the top. Use blend of VSH & MSIS LAYER= 3 PRESS= KNOT(NCOEFF+KORDER) - .01 ELSEIF (ALTSPEC .AND. ALT .LT. HT(1)) THEN C Altitude is too low LAYER= 1 PRESS= KNOT(1) + .01 ELSE C Find root of the polynomial PRESS= GETROOT(-5.5,4.5,IER,KNOT,HT,KORDER,NCOEFF,ALT,VAR) ENDIF OUTPUT(VAR)= PRESS ELSE C Pressure was specified in ALT PRESS= ALT ENDIF C Get temperature max IF (USING(3)) CALL SCATRAN (VSHT(1,1,1,ZTRUNC(3),3), COSPHI, 1 SINPHI, PBAR, MMAX, NMAX, MTRUNC(3), NTRUNC(3), TEMPMAX) IF (ABS(PRESS-LASTP) .GT. 1E-3) THEN C Precompute spline basis functions CALL BSPVN(KNOT,NCOEFF,KORDER,1,PRESS,KNOTLO,BSPLINE,WORK,IWORK) LASTP= PRESS ENDIF C Altitude transform DO 950 VAR= 1, FLDMAX IF (USING(VAR) .AND. (VAR .NE. 6 .OR. .NOT. ALTSPEC)) THEN IF (VAR .EQ. 3 .AND. KNOTLO .EQ. ZTRUNC(3)) THEN KTOP= KORDER-1 ELSE KTOP= KORDER ENDIF DO 900 PARITY= 1,2 DO 800 M= PARITY,MTRUNC(VAR) DO 700 N= M,NTRUNC(VAR) IF (VAR .NE. 4 .AND. VAR .NE. 5) THEN SUM= 0 DO 650 K= 1, KTOP SUM= SUM + VSHT(M,N,PARITY,K+KNOTLO-KORDER,VAR)*BSPLINE(K) 650 CONTINUE VSHTZ(M,N,PARITY,VAR)= SUM ELSE C Ion winds use only one pressure surface VSHTZ(M,N,PARITY,VAR)= 1 VSHT(M,N,PARITY,ZLOW,VAR) ENDIF 700 CONTINUE 800 CONTINUE 900 CONTINUE ENDIF 950 CONTINUE RETURN END SUBROUTINE SPACTR (VSHTZ,MMAX,NMAX,FLDMAX,COSPHI,SINPHI,PBAR, 1 VBAR, WBAR, OUTPUT, VARMAX, MTRUNC, NTRUNC, USING, ALTSPEC) INTEGER VAR, VARMAX, FLDMAX LOGICAL USING(VARMAX), ALTSPEC DIMENSION VSHTZ(MMAX, NMAX, 2, FLDMAX), OUTPUT(VARMAX) DIMENSION COSPHI(MMAX-1), SINPHI(MMAX-1) DIMENSION PBAR(MMAX+1,NMAX+1), VBAR(MMAX+1,NMAX+1) DIMENSION WBAR(MMAX+1,NMAX+1) DIMENSION MTRUNC(FLDMAX), NTRUNC(FLDMAX) C Spatial synthesis DO 950 VAR= 1,FLDMAX IF (USING(VAR) .AND. (VAR .NE. 6 .OR. .NOT. ALTSPEC)) THEN IF (VAR .EQ. 1 .OR. VAR .EQ. 4) THEN CALL VECTRAN(VSHTZ(1,1,1,VAR), VSHTZ(1,1,1,VAR+1), 1 MMAX, NMAX, COSPHI, SINPHI, VBAR, WBAR, 2 MTRUNC(VAR), NTRUNC(VAR), OUTPUT(VAR),OUTPUT(VAR+1)) ELSEIF (VAR .NE. 2 .AND. VAR .NE. 5) THEN CALL SCATRAN (VSHTZ(1,1,1,VAR), COSPHI, SINPHI, PBAR, 1 MMAX, NMAX, MTRUNC(VAR), NTRUNC(VAR), OUTPUT(VAR)) ENDIF ENDIF 950 CONTINUE RETURN END SUBROUTINE VECTRAN(DIV, ROT, MMAX, NMAX, COSPHI, SINPHI, 1 VBAR, WBAR, MTRUNC, NTRUNC, UWIND, VWIND) C Carry out vector spherical harmonic synthesis DIMENSION DIV (MMAX, NMAX, 2) DIMENSION ROT (MMAX, NMAX, 2) DIMENSION COSPHI(MMAX-1), SINPHI(MMAX-1) DIMENSION VBAR(MMAX+1,NMAX+1), WBAR(MMAX+1,NMAX+1) UWIND= 0. VWIND= 0. DO 75 MP1=1,MTRUNC DO 70 NP1= MP1,NTRUNC W= WBAR(MP1,NP1) V= VBAR(MP1,NP1) IF (MP1 .NE. 1) THEN UWIND= UWIND 1 - (ROT(MP1,NP1,1) * V + DIV(MP1,NP1,2) * W) * COSPHI(MP1-1) 2 + (ROT(MP1,NP1,2) * V - DIV(MP1,NP1,1) * W) * SINPHI(MP1-1) VWIND= VWIND 1 + (DIV(MP1,NP1,1) * V - ROT(MP1,NP1,2) * W) * COSPHI(MP1-1) 2 - (DIV(MP1,NP1,2) * V + ROT(MP1,NP1,1) * W) * SINPHI(MP1-1) ELSE UWIND = UWIND - (ROT(1,NP1,1) * V) / 2 VWIND = VWIND + (DIV(1,NP1,1) * V) / 2 ENDIF 70 CONTINUE 75 CONTINUE C Convert from spherical coordinates (V positive southwards) to C meteorological convention (V positive northwards) VWIND= -VWIND RETURN END SUBROUTINE SCATRAN (VSHTZ, COSPHI, SINPHI, PBAR, 1 MMAX, NMAX, MTRUNC, NTRUNC, SCALAR) C Carry out scalar spherical harmonic synthesis DIMENSION VSHTZ (MMAX, NMAX, 2) DIMENSION COSPHI(MMAX-1), SINPHI(MMAX-1) DIMENSION PBAR(MMAX+1,NMAX+1) SCALAR= 0. DO 75 MP1=1,MTRUNC DO 70 NP1= MP1,NTRUNC IF (MP1 .NE. 1) THEN SCALAR = SCALAR + PBAR(MP1,NP1) 1 * ( VSHTZ(MP1,NP1,1) * COSPHI(MP1-1) 2 - VSHTZ(MP1,NP1,2) * SINPHI(MP1-1)) ELSE SCALAR = SCALAR + PBAR(1,NP1) 1 * VSHTZ(1,NP1,1) / 2 ENDIF 70 CONTINUE 75 CONTINUE RETURN END SUBROUTINE IONELEC (OUTPUT,VARMAX, GLAT, PRESS, LAYER, USING) INTEGER VARMAX DIMENSION OUTPUT(VARMAX) LOGICAL USING(VARMAX) IF (USING(4)) THEN IF (LAYER .EQ. 1 .OR. ABS(GLAT) .LT. 1) THEN C Ion drifts: Correct lower levels for collisions OUTPUT(4)= 0 OUTPUT(5)= 0 ELSEIF (LAYER .EQ. 2 .AND. PRESS .LT. 0) THEN ALP= 2.6 * 25. * EXP(- PRESS) / 1 (1.38 * 9.57946 * 0.47 * OUTPUT(3) * SIN(GLAT/57.2958)) UI= OUTPUT(4) VI= OUTPUT(5) OUTPUT(4)= (UI + ALP * VI) / ( 1. + ALP**2 ) OUTPUT(5)= (VI - ALP * UI) / ( 1. + ALP**2 ) ENDIF ENDIF IF (USING(10)) OUTPUT(10)= EXP(OUTPUT(10)) IF (USING(11)) OUTPUT(11)= EXP(OUTPUT(11)) RETURN END SUBROUTINE DENSITY(OUTPUT, USING, VARMAX, PRESS, MSIS,DENSPEC,ALT) C Computes densities and mixing ratios INTEGER VARMAX LOGICAL USING(VARMAX), MSIS, DENSPEC DIMENSION OUTPUT(VARMAX) PARAMETER (PFAC= 6.0139E-12, PDIVBOL=3.62187E+12) C PFAC= (PNOUGHT=5.E-4)/(BOLTZ=1.3805E-16)/(AVOGAD=6.0225E+23) C PNOUGHT is reference pressure level C 3=Temp 8= O mixing ratio 13= O2 mixing ratio C 9= mass density 20= O density 21= 02 density 22= N2 density C 24= N2 mixing ratio IF (OUTPUT(13) .LT. 0) OUTPUT(13)= 0 IF (USING(13) .AND. OUTPUT(13) .EQ. 0) MSIS= .TRUE. IF (OUTPUT(8) .LT. 0) OUTPUT(8)= 0 IF (OUTPUT(8) .GT. 1) OUTPUT(8)= 1 IF (USING(24)) OUTPUT(24)= 1 -OUTPUT(8)-OUTPUT(13) IF (USING(9) .AND. (.NOT. DENSPEC)) THEN TDIVAMU= OUTPUT(9) DENFAC= EXP( -PRESS) / TDIVAMU OUTPUT(9)= PFAC * DENFAC ENDIF IF (DENSPEC) DENFAC= ALT/PFAC IF (USING(20)) OUTPUT(20)= DENFAC * OUTPUT(8) /16 * PDIVBOL IF (USING(21)) OUTPUT(21)= DENFAC * OUTPUT(13)/16 * PDIVBOL IF (USING(22)) OUTPUT(22)= DENFAC * OUTPUT(24)/32 * PDIVBOL RETURN END SUBROUTINE NORMLZE (GLAT, LAYER, USING, OUTPUT, VARMAX, 1 TEMPFAC, DENSFAC) INTEGER VARMAX DIMENSION OUTPUT(VARMAX), DENSFAC(2), TEMPFAC(2) LOGICAL USING(VARMAX) C Apply output field modifier to mass density and temperature IF (GLAT .GE. 10) THEN WEIGHT= 1 ELSEIF (GLAT .LE. -10) THEN WEIGHT= 0 ELSE C If within 10 degrees of equator, take weighted average of factors (to smooth) WEIGHT= (GLAT + 10) / 20 ENDIF IF (USING(9)) OUTPUT(9)= OUTPUT(9) * 1 (WEIGHT*DENSFAC(1) + (1-WEIGHT)*DENSFAC(2)) IF (USING(3)) OUTPUT(3)= OUTPUT(3) * 1 (WEIGHT*TEMPFAC(1) + (1-WEIGHT)*TEMPFAC(2)) END SUBROUTINE RUNPARAM(JDAY, AP, F107) INTEGER RUNMAX PARAMETER (RUNMAX=60) DIMENSION JDAYRUN(RUNMAX), APRUN(RUNMAX), F107RUN(RUNMAX) DATA JDAYRUN /21,52,111,141,202,233,264,294, 1 325,21,52,111,141,202,233,264,294,325,2*0,6*80,6*182,6*355, 2 355,80,180,19*80/ DATA F107RUN /9*67,9*243,2*0,3*67,3*243, 1 3*67,3*243,3*67,3*243,2*150,176,19*67/ c DATA F107A /9*72,73,212,212,8*0,3*72,3*215, c 1 3*72,3*215,3*72,3*215/ DATA APRUN /18*5,2*0,5,11,32,5,11,32,5,11,32, 1 5,11,32,5,11,32,5,11,32,22*5/ AP= APRUN(-JDAY) F107= F107RUN(-JDAY) JDAY= JDAYRUN(-JDAY) RETURN END SUBROUTINE GETMSIS(OUTPUT, USING, VARMAX, F107, JDAY, AP, 1 UT, ALT, GLAT, GLONG, ALTTOP, LAYER) C Get MSIS data if altitude is too high or if H density is requested INTEGER VARMAX, RUNMAX PARAMETER (RUNMAX=60) LOGICAL USING(VARMAX) DIMENSION OUTPUT(VARMAX),D(100),T(100),DTOP(100),TTOP(100) PARAMETER (AVOGAD= 6.0225E+23) IYD= 79000 + JDAY SEC= UT*3600 STL= GLONG/15 + UT IF (ALT .GT. 75) THEN HT= ALT ELSE HT= OUTPUT(6) ENDIF CALL GTS5(IYD, SEC,HT,GLAT,GLONG,STL,F107, F107, AP, 48, D, T) IF (LAYER .EQ. 3) THEN C Mixed layer WEIGHT= 1 - (HT - ALTTOP)/100 CALL GTS5(IYD, SEC,ALTTOP,GLAT,GLONG,STL,F107, F107, AP, 48, 1 DTOP, TTOP) IF (USING(9)) OUTPUT(9) = EXP( ALOG(D(6)) + 1 WEIGHT * ALOG( OUTPUT(9)/DTOP(6))) IF (USING(20)) OUTPUT(20) = EXP( ALOG(D(2)) + 1 WEIGHT * ALOG( OUTPUT(20)/DTOP(2))) IF (USING(21)) THEN IF (OUTPUT(21) .GT. 1E-6) THEN OUTPUT(21) = EXP( ALOG(D(4)) + 1 WEIGHT * ALOG( OUTPUT(21)/DTOP(4))) ELSE OUTPUT(21) = D(4) ENDIF ENDIF IF (USING(22)) OUTPUT(22) = EXP( ALOG(D(3)) + 1 WEIGHT * ALOG( OUTPUT(22)/DTOP(3))) ELSEIF (LAYER .EQ. 4 .OR. LAYER .EQ. 1) THEN C Above or below VSH altitudes OUTPUT(9) = D(6) OUTPUT(20)= D(2) OUTPUT(21)= D(4) OUTPUT(22)= D(3) IF (LAYER .EQ. 1) THEN DO 111 I= 1,7 OUTPUT(I)= 0 111 CONTINUE OUTPUT(3)= T(2) ENDIF ELSEIF (USING(21)) THEN IF (OUTPUT(21) .LT. 1E-6) OUTPUT(21)= D(4) ENDIF C Mixing ratios IF (USING(8) .AND. LAYER .NE. 2) THEN OUTPUT(8)= 16 * D(2) / D(6) / AVOGAD IF (OUTPUT(8) .LT. 0) OUTPUT(8) = 0 IF (OUTPUT(8) .GT. 1) OUTPUT(8) = 1 ENDIF IF (USING(13) .AND. LAYER .NE. 2) THEN OUTPUT(13)= 28 * D(3) / D(6) / AVOGAD IF (OUTPUT(13) .LT. 0) OUTPUT(13) = 0 ENDIF C Hydrogen OUTPUT(23)= D(7) RETURN END