subroutine f2_trough
a (lunmap,month,uthour,hourlt,idayno,flat,flon,foe_solar,
a tindex,f107a,f107,ap,fof2,ymf2,hmf2,iflag)
c
c fof2 & hmf2 are OVERWRITTEN (by trough values)
c
c This s/r is a revised version of the ICED s/r F2LAYR, with
c corrections, and a change to the use of the IPS world maps
c
c written by leo mcnamara 5-nov-90
c
c MODIFIED 8-OCT-91 to set foe = 0.6 without setting foe_solar
c to 0.6 - that was stupid anyway. Note that this means that
c foe_solar is no longer used
c
c----- last modified 8-oct-91
c
ccc SUBROUTINE F2LAYR(IFLAG,FOF2,HMAX,FOE)
C***********************************************************************
C* SUBROUTINE NAME : F2LAYR CPC: SS/MOD/RTN/F2LAYR *
C* VERSION: 88-II DATE: 16 MAR 88 *
C* *
C* PURPOSE: THIS SUBROUTINE IS THE ELKINS-RUSH MODIFICATION *
C* TO THE URSI MODEL. GIVEN A SSN AND A KP/QE AS *
C* INPUT, IT WILL CALCULATE THE POSITION OF THE *
C* POLEWARD TROUGH AND MAKE APPROPRIATE MODIFIC- *
C* ATIONS OF THE TROUGH AND AURORAL FOF2 VALUES. *
C* *
C* CALLING *
C* ARGUMENTS: IFLAG,FOF2,HMAX,FOE *
C* *
C* ABBREVIATIONS: NONE *
C* *
C* COMMON BLOCKS: F6PARM, GWC1 *
C* *
C* FILES ACCESSED: NONE *
C* *
C* METHOD: 1. CALCULATE THE LOCAL TIME AND THE POSITION *
C* OF THE TROUGH WALL. *
C* 2. DETERMINE THE CORRECTION FACTORS FOR THE *
C* POLAR CAVITY,AURORAL ZONE, AND AURORAL TROUGH*
C* 3. CALCULATE FOF2 AND HMAX AT THE GRIDPOINT *
C* *
C* REMARKS: NONE *
C* *
C* REFERENCES: 1. ICED-III SYSTEM DOCUMENTATION *
C* 2. SPERRY FORTRAN-77 (ASCII) PRM *
C* 3. AFGWC FORM 10 DX-50194 *
C* 4. AFGWC FORM 10 DX-60096 *
C* *
C* GLOBAL VARIABLES: *
C* COMMON /F6PARM/IYR,IMO,IDA,TIME,DAYJUL,EFFSSN,EFFQ,EFFKP, *
C* GGLAT,GGLON,CGLAT,CGLON,TCGM *
C* IYR IS YEAR OF DATA SET *
C* IMO IS MONTH OF DATA SET *
C* IDA IS DAY OF DATA SET *
C* TCGM RS CORRECTED GEOMAGNETIC TIME OF A GRIDPOINT *
C* TIME RS UNIVERSAL TIME IN HH.HH OF DATA SET *
C* DAYJUL RS JULIAN DATE OF DATA SET *
C* EFFSSN RS EFFECTIVE SUNSPOT NUMBER *
C* EFFQ RS EFFECTIVE AURORAL Q VALUE *
C* EFFKP RS EFFECTIVE PLANETARY K INDEX *
C* GGLAT RS GEOGRAPHIC LATITUDE OF A GRIDPOINT *
C* GGLON RS GEOGRAPHIC LONGITUDE OF A GRIDPOINT *
C* CGLAT RS CORRECTED GEOMAGNETIC LATITUDE OF A GRIDPOINT *
C* CGLON RS CORRECTED GEOMAGNETIC LONGITUDE OF A GRIDPOINT *
C* *
C* COMMON /GWC1/ KF(10),UF(13,76),KX(10),UX(13,76) *
C* KF IV HEADER RECORD VALUES FOR FOF2 COEFFICIENTS IN URSI *
C* UF RA FOURIER SERIES COEFFICIENTS FOR FOF2 IN URSI *
C* KX IV HEADER RECORD VALUES FOR M3000 COEFFICIENTS IN URSI *
C* UX RA FOURIER SERIES COEFFICIENTS FOR M3000 IN URSI *
C* *
C* LOCAL VARIABLES: IFLAG IS QUALIFIER AT GRIDPOINT *
C* IERR IS ERROR FLAG FROM COORDINATE TRANSFORM*
C* CHI RS SOLAR ZENITH ANGLE *
C* DELTAN RS SMALL INCREMENT OF TIME *
C* DELTAT RS MLT HOURS FROM 0300 MLT *
C* EB RS EQUATORIAL BOUNDARY VALUE *
C* DUM RS PLACEHOLDER IN CALL TO SOLPOS *
C* FOE RS FREQUENCY OF E LAYER *
C* FOF RS TEMP VALUE OF FREQUENCY OF F2 LAYER *
C* FOF2 RS FREQUENCY OF THE FOF2 LAYER *
C* HMAX RS HEIGHT OF THE F2 LAYER *
C* PB RS POLAR BOUNDARY OF AURORAL ZONE *
C* PHIA RS AURORAL ZONE EQUATORWARD BOUNDARY *
C* PI RS 3.1415926 - CONSTANT *
C* RDUM RS UNUSED RETURN VALUE FROM COORDINATE *
C* TRANSFORMATION *
C* T RS TEMP TIME VARIABLE *
C* T1 RS TEMP TIME VARIABLE *
C* TIMLOC RS LOCAL TIME AT GGLON *
C* TCGM RS CORRECTED GEOMAGNETIC TIME *
C* X1 RS TEMP VARIABLE *
C* XA RS TEMP VARIABLE *
C* TRFPB RS POLEWARE BOUNDARY OF THE TROUGH *
C* TRFMAX RS CENTER OF THE TROUGH *
C* TRFEB RS EQUATORWARD BOUNDARY OF THE TROUGH *
C* HFRM20 RS MLT HOURS PAST 20 MLT *
C* CGLN20 RS CORRECTED GEOMAG LONG. AT 20 MLT *
C* CGLN03 RS CORRECTED GEOMAG LONG. AT 03 MLT *
C* CGLN07 RS CORRECTED GEOMAG LONG. AT 07 MLT *
C* GLAT RS TEMPORARY GEOGRAPHIC LATITUDE *
C* GLON RS TEMPORARY GEOGRAPHIC LONGITUDE *
C* M3000 RS TEMPORARY M(3000) FACTOR *
C* HMAXTL RS F2 HEIGHT AT TOP LEFT OF *
C* INTERPOLATION BOX *
C* HMAXBL RS F2 HEIGHT AT BOTTOM LEFT OF *
C* INTERPOLATION BOX *
C* HMAXTR RS F2 HEIGHT AT TOP RIGHT OF *
C* INTERPOLATION BOX *
C* HMAXBR RS F2 HEIGHT AT BOTTOM RIGHT OF *
C* INTERPOLATION BOX *
C* FRLAT RS FRACTIONAL LAT. FROM BOTTOM TO TOP *
C* OF INTERPOLATION BOX *
C* FRLON RS FRACTIONAL LONG. FROM LEFT TO RIGHT *
C* OF INTERPOLATION BOX *
C* PT1 RS TEMP. BOTTOM F2 HEIGHT IN LONGITUDE *
C* INTERPOLATION BOX *
C* PT2 RS TEMP. TOP F2 HEIGHT IN LONGITUDE *
C* INTERPOLATION BOX *
C* F RS M3000 CALCULATION VARIABLE *
C* XE RS RATIO OF FOF2 TO FOE *
C* DELTAM RS M3000 CORRECTION TERM *
C* (I=INTEGER, R=REAL, C=CHARACTER, S=SCALAR, V=VECTOR, A=ARRAY) *
C* *
C* SETC OPTIONS: NONE USED *
C* *
C* SUBPROGRAMS *
C* CALLED: QOVAL1 - CALCULATES THE AURORAL OVAL POSITION *
C* FLDHR - RETRIEVES FOURIER ANALYSIS COEFFICIENTS *
C* EQBOUN - CALCULATES THE EQ. BOUNDARY OF OVAL *
C* CNVCRD - CONVERTS GEOMAG TO GEOGRAPHIC COORDS *
C* SOLPOS - COMPUTES SOLAR ZENITH ANGLE *
C* *
C* SUBPROGRAM *
C* CALLED BY : FIELD6 *
C* *
C* PROGRAM WRITTEN: JAN 86 - NATIONAL GEOPHYSICAL DATA CENTER *
C* MAY 86 - AFGWC/SDDE - UPGRADED DOCUMENTATION *
C* AUG 87 - NATIONAL GEOPHYSICAL DATA CENTER (87-I) *
C* MAR 88 - NGDC (88-II) NEW SOLAR ZENITH ANGLE *
C* JAN 89 - NGDC (89) *
C* APR 90 - NGDC (90-REF) *
C* *
C***********************************************************************
INTEGER IDA
INTEGER IMO
INTEGER IYR
REAL CGLAT
REAL CGLON
REAL DAYJUL
REAL EFFKP
REAL EFFQ
REAL EFFSSN
REAL GGLAT
REAL GGLON
REAL TCGM
REAL TIME
COMMON/F6PARM/IYR,IMO,IDA,TIME,DAYJUL,EFFSSN,EFFQ,EFFKP,
+ GGLAT,GGLON,CGLAT,CGLON,TCGM
INTEGER KF
INTEGER KX
REAL UF
REAL UX
COMMON/GWC1/KF(10),UF(13,76),KX(10),UX(13,76)
INTEGER IFLAG
INTEGER IERR
REAL CGLN20
REAL CGLN03
REAL CGLN07
REAL CHI
REAL DELTAM
REAL DELTAN
REAL DELTAT
REAL EB
REAL DUM
REAL F
REAL FOE
REAL FOF
REAL FOF2
REAL FRLAT
REAL FRLON
REAL GLAT
REAL GLON
REAL HFRM20
REAL HMAX
REAL HMAXBL
REAL HMAXBR
REAL HMAXTL
REAL HMAXTR
REAL M3000
REAL PB
REAL PHIA
REAL PI
REAL PT1
REAL PT2
REAL RDUM
REAL T
REAL T1
REAL TEMPDM
REAL TIMLOC
REAL TRFEB
REAL TRFMAX
REAL TRFPB
REAL X1
REAL XA
REAL XE
DATA PI/3.141593/
dimension ap(1)
c
c We are going to modify the supplied values of fof2 and hmf2,
c so we need to rename them to avoid confusion.
c Note that the trough and auroral enhancements just tail off
c exponentially, so we end up modifying all supplied values of
c fof2, even if by just a tiny bit
fof2_in = fof2
hmf2_in = hmf2
c We may not in fact modify hmf2 (or fof2 ??)
fof2_mod = fof2_in
hmf2_mod = hmf2_in
c
c We need the solar value of foE to use for the Dudeney
c formula for the calculationof hmf2 in the trough
c Since the trough is a nightime feature, we set foe = 0.6,
c and NOT USE THE SUPPLIED VALUE OF FOE_SOLAR
foe = 0.6
c
C.......................................................................
C AURORAL ZONE BOUNDARIES (CENTER OF AZ = EQ BOUNDARY OF TROUGH) .
C.......................................................................
CALL QOVAL1(EFFQ,TCGM,PB,EB,PHIA)
C.......................................................................
C DETERMINE THE CORRECTION FACTOR IN AURORAL ZONE AND POLAR CAVITY .
C.......................................................................
DELTAT = ABS(TCGM-3.)
IF (TCGM.GT.15.0) DELTAT = 27. - TCGM
IF (CGLAT.LT.PHIA) GO TO 40
XA = PB-PHIA
X1 = (PHIA-CGLAT)/XA
IF (CGLAT.GT.EB) X1 = X1 * 2.0
DELTAN=EXP((-1.0*X1*X1)/2.0)
c mods to soften the discontinuity at phia
x1 = (cglat - phia) / xa
if(x1.ge.0.25) deltan = exp( 3./32.) * exp (-(2.*x1)**2 / 2.)
if(x1.lt.0.25) deltan = exp(-1./32.) * sin(pi/2. * x1/0.25)
GO TO 80
C.......................................................................
C CALCULATE THE CORRECTION FACTOR SOUTH OF AURORAL ZONE .
C.......................................................................
40 T1 = 0.0
XA = 3.7 + (1.3*EFFKP)
C.......................................................................
C COMPUTE SOLAR ZENITH ANGLE .
C.......................................................................
TIMLOC = TIME + GGLON/15.0
IF ( TIMLOC.GE.24.00 ) TIMLOC = TIMLOC - 24.00
IF ( TIMLOC.LT. 0.00 ) TIMLOC = TIMLOC + 24.00
CALL SOLPOS( IYR,IMO,IDA,TIMLOC,GGLAT,GGLON,CHI,DUM,DUM,DUM )
X1=(CGLAT-PHIA)/XA
IF (CHI.LE.90.0.OR.(TCGM.GT.6.0.AND.TCGM.LT.18.0)) T1 = 0.0
IF ((CHI.GT.94.6).AND.((18.0.LE.TCGM).OR.(TCGM.LE.6.0))) T1 =-0.2
IF ( ((CHI.GE.90.0).AND.(CHI.LE.94.6)).AND.
+ ((18.0.GE.TCGM).OR.(TCGM.LE.6.0)) ) T1=-0.2 *(CHI-90.0)/4.6
C.....................................................................
C MILLER AND GIBBS 1975 .
C.....................................................................
c the latitudinal variation within the trough
c Note that XA > 0; X1 < 0.
fac_lat = EXP((X1-(X1*X1))/2.0)
c mod to remove the discontinuity at the poleward edge
if(abs(x1).le.0.5)
a fac_lat = exp(-3./8.) * sin(pi/2. * abs(x1)/0.5)
c
c===================================================================
c deepen the trough by using -x1
fac_lat = EXP((-X1-(X1*X1))/2.0)
c mod to remove the discontinuity at the poleward edge
if(abs(x1).le.0.5)
a fac_lat = exp(1./8.) * sin(pi/2. * abs(x1)/0.5)
c===================================================================
c
T=T1* fac_lat *EXP((-1.0*DELTAT*DELTAT)/12.0)
DELTAN=T*(1.0+COS(PI*(DAYJUL+11.0)/182.5))
c
cccc80 CALL FLDHR(KF,UF,GGLAT,GGLON,FOF,TIME)
c we have already calculated fof2
80 continue
c
c Trough value of fof2
if(cglat.lt.phia) FOF2_mod = FOF2_in * (1.0+DELTAN)
c
C......................................................................
C 'DENSITY' CORRECTION IN AURORAL ZONE AND POLAR CAVITY CONSIDERED.
C CONSTANT IN TIME. CORRECTION BASED ON DENSITY DIFFERENCE BETWEEN.
C 4 AND 5 MHZ. 5 MHZ BEING APPROX 1.22663 TIMES 4 MHZ WHICH IS A .
C TYPICAL NIGHTTIME URSI FOF2 VALUE. .
C......................................................................
ccc IF (CGLAT.GE.PHIA) FOF2 = SQRT( (FOF*FOF) + (9.0+DELTAN) )
c changed lfm 5-nov-90 |
IF (CGLAT.GE.PHIA) FOF2_mod = SQRT((FOF2_in **2) + (9.0 * DELTAN))
c
C......................................................................
C IFLAG >= 100 IMPLIES GRIDPOINT IN THE SUB-AURORAL TROUGH .
c If so, we need to calculate the modified value of hmf2
C......................................................................
IF ( IFLAG.LT.100 ) GO TO 170
c
C......................................................................
C F2 TROUGH HEIGHT CALCULATION BY INTERPOLATION. .
C INTERPOLATION CASES : .
C |-------------------------|-----------------------| TRFPB .
C | | | .
C | 1 | 3 | ~PHIA .
C | | | .
C |-------------------------*-----------------------| TRFMAX .
C | | | .
C | 2 | 4 | .
C | | | .
C |-------------------------|-----------------------| TRFEB .
C 20 MLT 03 MLT 07 MLT .
C .
C * = 450.0 KILOMETERS (MAX HEIGHT IN TROUGH) .
C......................................................................
TRFMAX = PHIA - 1.5
TRFPB = PHIA + 1.5
CALL EQBOUN( TCGM,TRFMAX,DAYJUL,EFFKP,TRFEB )
IF ( TRFEB.GT.TRFMAX-1.5 ) TRFEB = TRFMAX - 1.5
c
c This code plays the same role as that which checks to see
c if flag > 100
ccc if(cglat.lt.trfeb) go to 170
ccc if(cglat.gt.trfpb) go to 170
ccc if(tcgm.gt.7.0.and.tcgm.lt.20.0) go to 170
c
C......................................................................
C HOURS PAST 20 MLT .
C......................................................................
HFRM20 = 4.0 + TCGM
IF ( TCGM.GE.20.0 ) HFRM20 = TCGM - 20.0
C.......................................................................
C CG LONGITUDES OF 20,03,07 MLT .
C.......................................................................
CGLN20 = CGLON - ( HFRM20 * 15.0 )
CGLN03 = CGLN20 + 105.0
CGLN07 = CGLN20 + 165.0
IF ( CGLN20.LE. 0.0 ) CGLN20 = CGLN20 + 360.0
IF ( CGLN03.GE.360.0 ) CGLN03 = CGLN03 - 360.0
IF ( CGLN07.GE.360.0 ) CGLN07 = CGLN07 - 360.0
c
c the height is given by f2layer, so we do not need this
c
ccc XE = FOF2 / FOE
ccc IF ( XE.LT.1.7 ) XE = 1.7
ccc DELTAM = ( 0.253/(XE-1.215) ) - 0.012
c
C......................................................................
C INTERPOLATION CASE 1 .
C......................................................................
IF ( (CGLAT.GE.TRFMAX).AND.(HFRM20.LE.7.) ) THEN
CALL CNVCRD( TRFPB,CGLN20,0.,4,GLAT,GLON,RDUM,2,IERR)
call F2layer (lunmap,month,uthour,hourlt,idayno,glat,glon,
a foe,tindex,f107a,f107,ap,FOF,ymf2,hmf2)
hmaxtl = hmf2
CALL CNVCRD( TRFMAX,CGLN20,0.,4,GLAT,GLON,RDUM,2,IERR)
call F2layer (lunmap,month,uthour,hourlt,idayno,glat,glon,
a foe,tindex,f107a,f107,ap,FOF,ymf2,hmf2)
hmaxbl = hmf2
CALL CNVCRD( TRFPB,CGLN03,0.,4,GLAT,GLON,RDUM,2,IERR)
call F2layer (lunmap,month,uthour,hourlt,idayno,glat,glon,
a foe,tindex,f107a,f107,ap,FOF,ymf2,hmf2)
hmaxtr = hmf2
HMAXBR = 450.0
FRLAT = (CGLAT-TRFMAX) / (TRFPB-TRFMAX)
FRLON = HFRM20 / (27.-20.)
ENDIF
C......................................................................
C INTERPOLATION CASE 2 .
C......................................................................
IF ( (CGLAT.LT.TRFMAX).AND.(HFRM20.LE.7.) ) THEN
CALL CNVCRD( TRFMAX,CGLN20,0.,4,GLAT,GLON,RDUM,2,IERR)
call F2layer (lunmap,month,uthour,hourlt,idayno,glat,glon,
a foe,tindex,f107a,f107,ap,FOF,ymf2,hmf2)
hmaxtl = hmf2
CALL CNVCRD( TRFEB,CGLN20,0.,4,GLAT,GLON,RDUM,2,IERR)
call F2layer (lunmap,month,uthour,hourlt,idayno,glat,glon,
a foe,tindex,f107a,f107,ap,FOF,ymf2,hmf2)
hmaxbl = hmf2
HMAXTR = 450.0
CALL CNVCRD( TRFEB,CGLN03,0.,4,GLAT,GLON,RDUM,2,IERR)
call F2layer (lunmap,month,uthour,hourlt,idayno,glat,glon,
a foe,tindex,f107a,f107,ap,FOF,ymf2,hmf2)
hmaxbr = hmf2
FRLAT = (CGLAT-TRFEB) / (TRFMAX-TRFEB)
FRLON = HFRM20 / (27.-20.)
ENDIF
C......................................................................
C INTERPOLATION CASE 3 .
C......................................................................
IF ( (CGLAT.GE.TRFMAX).AND.(HFRM20.GT.7.) ) THEN
CALL CNVCRD( TRFPB,CGLN03,0.,4,GLAT,GLON,RDUM,2,IERR)
call F2layer (lunmap,month,uthour,hourlt,idayno,glat,glon,
a foe,tindex,f107a,f107,ap,FOF,ymf2,hmf2)
hmaxtl = hmf2
HMAXBL = 450.0
CALL CNVCRD( TRFPB,CGLN07,0.,4,GLAT,GLON,RDUM,2,IERR)
call F2layer (lunmap,month,uthour,hourlt,idayno,glat,glon,
a foe,tindex,f107a,f107,ap,FOF,ymf2,hmf2)
hmaxtr = hmf2
CALL CNVCRD( TRFMAX,CGLN07,0.,4,GLAT,GLON,RDUM,2,IERR)
call F2layer (lunmap,month,uthour,hourlt,idayno,glat,glon,
a foe,tindex,f107a,f107,ap,FOF,ymf2,hmf2)
hmaxbr = hmf2
FRLAT = (CGLAT-TRFMAX) / (TRFPB-TRFMAX)
FRLON = (HFRM20-7.) / (31.-27.)
ENDIF
C......................................................................
C INTERPOLATION CASE 4 .
C......................................................................
IF ( (CGLAT.LT.TRFMAX).AND.(HFRM20.GT.7.) ) THEN
HMAXTL = 450.0
CALL CNVCRD( TRFEB,CGLN03,0.,4,GLAT,GLON,RDUM,2,IERR)
call F2layer (lunmap,month,uthour,hourlt,idayno,glat,glon,
a foe,tindex,f107a,f107,ap,FOF,ymf2,hmf2)
hmaxbl = hmf2
CALL CNVCRD( TRFMAX,CGLN07,0.,4,GLAT,GLON,RDUM,2,IERR)
call F2layer (lunmap,month,uthour,hourlt,idayno,glat,glon,
a foe,tindex,f107a,f107,ap,FOF,ymf2,hmf2)
hmaxtr = hmf2
CALL CNVCRD( TRFEB,CGLN07,0.,4,GLAT,GLON,RDUM,2,IERR)
call F2layer (lunmap,month,uthour,hourlt,idayno,glat,glon,
a foe,tindex,f107a,f107,ap,FOF,ymf2,hmf2)
hmaxbr = hmf2
FRLAT = (CGLAT-TRFEB) / (TRFMAX-TRFEB)
FRLON = (HFRM20-7.) / (31.-27.)
ENDIF
C......................................................................
C TWO-DIMENSIONAL LINEAR INTERPOLATION BETWEEN FOUR CORNERS .
C......................................................................
PT1 = HMAXBL + FRLON*(HMAXBR-HMAXBL)
PT2 = HMAXTL + FRLON*(HMAXTR-HMAXTL)
HMAX = PT1 + FRLAT*(PT2-PT1)
hmf2_mod = hmax
ccc RETURN
C......................................................................
C NON-TROUGH HMF2 COMPUTATION ----- deleted
C......................................................................
c
170 continue
c
fof2 = fof2_mod
hmf2 = hmf2_mod
c
RETURN
END