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