-
Notifications
You must be signed in to change notification settings - Fork 11
/
iri90.f
2994 lines (2951 loc) · 104 KB
/
iri90.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
C Subroutine IRI90
C
C Special-purpose version of IRI90 for use with the GLOW model.
C This is not the original IRI90.
C
C Adapted 7/93 from 10/91 version of IRIS12 by Stan Solomon.
C Replaces 3/92 adaptation of 12/90 version of IRIS12.
C Replaced height range and interval with array of heights ZKM
C Also supply number of heights NZ; unlimited number of heights allowed.
C Added DIRECT argument to specify location of CCIR and URSI files.
C Uses ASCII versions of CCIR and URSI files.
C Subroutine DFP splices directory and filename together.
C Accepts longitudes -180 to +180; maps internal longitudes 0-360.
C Added SAVE statements for volatile memory machines.
C Changed 'file not found' error to output message on unit 6 and stop.
C Changed JF(12)=.true. to write messages on unit 12 instead of 6.
C Changed call to IONCOM to use floating point variable ZMONTH.
C Note: JF(5)=.false. uses recommended (URSI) coefficients.
C Note: JF(4)=.false. uses recommended Gulyeava B0 coefficients.
C All other JF's usually should be .true.
C
C The following is the original comment from IRIS12:
C
C IRIS12.FOR ---------------------------------------- OCTOBER 1991
C
C*****************************************************************
C CHANGES FROM IRIS11.FOR TO IRIS12.FOR:
C - CIRA-1986 INSTEAD OF CIRA-1972 FOR NEUTRAL TEMPERATURE
C - 10/30/91 VNER FOR NIGHTTIME LAY-VERSION: ABS(..)
C - 10/30/91 XNE(..) IN CASE OF LAY-VERSION
C - 10/30/91 CHANGE SSIN=F/T TO IIQU=0,1,2
C - 10/30/91 Te > Ti > Tn ENFORCED IN FINAL PROFILE
C - 10/30/91 SUB ALL NAMES WITH 6 OR MORE CHARACTERS
C - 10/31/91 CORRECTED HF1 IN HST SEARCH: NE(HF1)>NME
C------------- inlcuded on diskette ------------------------------
C - 11/14/91 C1=0 IF NO F1-REGION
C - 11/14/91 CORRECTED HHMIN AND HZ FOR LIN. APP.
C - 1/28/92 RZ12=0 included
C - 1/29/92 NEQV instead of NE between URSIF2 and URSIFO
C - 5/ 1/92 CCIR and URSI input as in IRID12
C
C*****************************************************************
C********* INTERNATIONAL REFERENCE IONOSPHERE (IRI). *************
C*****************************************************************
C**************** OCTOBER 1991 ****************************
C**************** SUBROUTINE ****************************
C*****************************************************************
C
C
C-----------------------------------------------------------------
C INTERNATIONAL REFERENCE IONOSPHERE 1991
C
C INPUT: JMAG=0/1 GEODETIC/GEOMAGNETIC LATITUDE AND LONGITUDE
C ALATI,ALONG LATITUDE NORTH AND LONGITUDE EAST IN DEGREES
C RZ12 (-COV) 12-MONTHS-RUNNING MEAN OF SOLAR SUNSPOT NUMBER
C (OR EQUIVALENT F10.7 SOLAR RADIO FLUX AS
C NEGATIVE NUMBER)
C MMDD (-DDD) DATE (OR DAY OF YEAR AS A NEGATIVE NUMBER)
C DHOUR LOCAL TIME (OR UNIVERSAL TIME + 25) IN DECIMAL
C HOURS
C HEIBEG, BEGIN, END, AND STEPWIDTH OF HEIGHT RANGE
C HEIEND,HEISTP IN KM (MAXIMUM NUMBER OF STEPS IS 50 !!)
C JF(1:12) TRUE/FALSE FLAGS FOR SEVERAL OPTIONS
C JF(1)=.TRUE.[.FALSE.] ELECTRON DENSITY IS [NOT] CALCULATED
C JF(2)=T[F] TEMPERATURES ARE [NOT] CALCULATED
C JF(3)=T[F] ION COMPOSITION IS [NOT] CALCULATED
C JF(4)=T[F] B0 FROM TABLE [FROM GULYEAVA 1987]
C JF(5)=T[F] F2 PEAK FROM CCIR [FROM URSI]
C JF(6)=T[F] ION COMP. STANDARD [DANILOV-YAICHNIKOV-1985]
C JF(7)=T[F] STAND. IRI TOPSIDE [IRI-79]
C JF(8)=T[F] NMF2 PEAK MODEL [INPUT VALUES]
C JF(9)=T[F] HMF2 PEAK MODEL [INPUT VALUES]
C JF(10)=T[F] TE MODEL [TE-NE MODEL WITH NE INPUT]
C JF(11)=T[F] NE STANDARD [LAY-FUNCTIONS VERSION]
C JF(12)=T[F] MESSAGE ARE WRITTEN TO UNIT=12 [=6]
C
C JF(1:11)=.TRUE. GENERATES THE STANDARD IRI-90 PARAMETERS.
C IF YOU SET JF(8)=.FALSE., THAN YOU HAVE TO PROVIDE THE F2 PEAK
C NMF2/M-3 OR FOF2/MHZ IN OARR(1). SIMILARLY, IF YOU SET JF(9)=
C .FALSE., THAN YOU HAVE TO PROVIDE THE F2 PEAK HEIGHT HMF2/KM IN
C OARR(2). IF YOU SET JF(10)=.FALSE., THAN YOU HAVE TO PROVIDE THE
C ELECTRON DENSITY IN M-3 AT 300KM AND/OR 400KM AND/OR 600KM IN
C OARR(3), OARR(4), AND OARR(5). IF YOU WANT TO USE THIS OPTION AT
C ONLY ONE OF THE THREE ALTITUDES, THAN SET THE DENSITIES AT THE
C OTHER TWO TO ZERO.
C
C OUTPUT: OUTF(1:10,1:50) IRI PROFILES
C OUTF(1,*) ELECTRON DENSITY/M-3
C OUTF(2,*) NEUTRAL TEMPERATURE/K
C OUTF(3,*) ION TEMPERATURE/K
C OUTF(4,*) ELECTRON TEMPERATURE/K
C OUTF(5,*) PERCENTAGE OF O+ IONS IN %
C OUTF(6,*) PERCENTAGE OF H+ IONS IN %
C OUTF(7,*) PERCENTAGE OF HE+ IONS IN %
C OUTF(8,*) PERCENTAGE OF O2+ IONS IN %
C OUTF(9,*) PERCENTAGE OF NO+ IONS IN %
C AND, IF JF(6)=.FALSE.:
C OUTF(10,*) PERCENTAGE OF CLUSTER IONS IN %
C OUTF(11,*) PERCENTAGE OF N+ IONS IN %
C
C OARR(1:30) ADDITIONAL OUTPUT PARAMETERS
C OARR(1) = NMF2/M-3 OARR(2) = HMF2/KM
C OARR(3) = NMF1/M-3 OARR(4) = HMF1/KM
C OARR(5) = NME/M-3 OARR(6) = HME/KM
C OARR(7) = NMD/M-3 OARR(8) = HMD/KM
C OARR(9) = HHALF/KM OARR(10) = B0/KM
C OARR(11) =VALLEY-BASE/M-3 OARR(12) = VALLEY-TOP/KM
C OARR(13) = TE-PEAK/K OARR(14) = TE-PEAK HEIGHT/KM
C OARR(15) = TE-MOD(300KM) OARR(16) = TE-MOD(400KM)/K
C OARR(17) = TE-MOD(600KM) OARR(18) = TE-MOD(1400KM)/K
C OARR(19) = TE-MOD(3000KM) OARR(20) = TE(120KM)=TN=TI/K
C OARR(21) = TI-MOD(430KM) OARR(22) = X/KM, WHERE TE=TI
C OARR(23) = SOLAR ZENITH ANGLE/DEG
C OARR(24) = SUN DECLINATION/DEG
C OARR(25) = DIP
C OARR(26) = DIP LATITUDE
C OARR(27) = MODIFIED DIP LATITUDE
C OARR(28:30) FREE
C-------------------------------------------------------------------
C*** THIS PROGRAM PRODUCES PROFILES OF ***
C*** ELECTRON DENSITY ***
C*** NEUTRAL TEMPERATURE (CIRA 86) ***
C*** ELECTRON TEMPERATURE ***
C*** ION TEMPERATURE ***
C*** RELATIVE PERCENTAGE DENSITIES OF THE IONS ***
C*** ATOMIC OXYGEN, HYDROGEN, HELIUM, ***
C*** MOLECULAR OXYGEN AND NITROGEN OXYD (NO+) ***
C*****************************************************************
C*** THE ALTITUDE LIMITS ARE: LOWER (DAY/NIGHT) UPPER ***
C*** ELECTRON DENSITY 60/80 KM 1000 KM ***
C*** TEMPERATURES 120 KM 3000 KM ***
C*** ION DENSITIES 100 KM 1000 KM ***
C*****************************************************************
C* --------------------ADDRESSES------------------------ *
C* I PROF. K. RAWER DR. D. BILITZA I *
C* I HERRENSTR. 43 GSFC/NSSDC CODE 633 I *
C* I D-7801 MARCH GREENBELT MD 20771 I *
C* I F.R.G. USA I *
C* ----------------------------------------------------- *
C*****************************************************************
C*****************************************************************
C*****************************************************************
C********* ALL ANGLES ARE IN DEGREE **************
C********* ALL DENSITIES ARE IN M-3 **************
C********* ALL ALTITUDES ARE IN KM **************
C********* ALL TEMPERATURES ARE IN KELVIN **************
C********* ALL TIMES ARE IN DECIMAL HOURS **************
C*****************************************************************
C******************** OPTIONS **********************************
C*****************************************************************
C* FOR HMF2=0 OR FOF2=0 THE F2 PEAK VALUES ARE CALCULATED WITH *
C* THE CCIR OR URSI MODELS. THE CCIR COEFFICIENT SET FOR THE *
C* MONTH "mm" IS EXPECTED IN THE BINARY FILE "CCIRmm.BIN" AND *
C* THE URSI SET IN "URSImm.BIN". IF YOU USE THE ASCII CODED *
C* FILES "CCIRmm.ASC", YOU HAVE TO INCORPORATE THE CHANGES *
C* INDICTED IN PROGRAM SECTION ENTITLED "READ CCIR COEFFICIENT *
C* SET FOR CHOSEN MONTH." *
C*****************************************************************
C*****************************************************************
C*****************************************************************
C
SUBROUTINE IRI90(JF,JMAG,ALATI,ALONG,RZ12,MMDD,DHOUR,
& ZKM,NZ,DIRECT,OUTF,OARR)
dimension zkm(nz), outf(11,nz), oarr(30)
character*(*) direct
character*50 path
character*10 filename
INTEGER EGNR,AGNR,DAYNR,DDO,DO2,SEASON,SEADAY
REAL LATI,LONGI,MO2,MO,MODIP,NMF2,MAGBR
REAL NMF1,NME,NMD,NEI,MM,MLAT,MLONG,NOBO2
DIMENSION F(3),RIF(4),E(4),XDELS(4),DNDS(4)
DIMENSION FF0(988),XM0(441),F2(13,76,2),FM3(9,49,2)
DIMENSION AMP(4),HXL(4),SCL(4),B0B1(5)
DIMENSION CTN(3),CTNN(3),XSM(4),MM(5),DTI(4)
DIMENSION AHH(7),STTE(6),DTE(5),ATE(7),TEA(6),HOA(3),XNAR(3)
DIMENSION PG1O(80),PG2O(32),PG3O(80),PF1O(12),PF2O(4),PF3O(12)
DIMENSION HO(4),MO(5),DDO(4),HO2(2),MO2(3),DO2(2),DION(7)
LOGICAL EXT,SCHALT,NIGHT,TCON(3)
LOGICAL F1REG,FOF2IN,HMF2IN,URSIF2,LAYVER,DY,GULB0
LOGICAL NODEN,NOTEM,NOION,TENEOP
LOGICAL OLD79,TOPSI,BOTTO,BELOWE,JF(12),URSIFO
COMMON /BLOCK1/HMF2,NMF2,HMF1 /CONST/UMR
& /BLOCK2/B0,B1,C1 /BLOCK3/HZ,T,HST,STR
& /BLOCK4/HME,NME,HEF /BLOCK5/NIGHT,E
& /BLOCK6/HMD,NMD,HDX /BLOCK7/D1,XKK,FP30,FP3U,FP1,FP2
& /BLOCK8/HS,TNHS,XSM,MM,DTI,MXSM
& /BLOTN/XSM1,TEXOS,TLBDH,SIGMA /BLOTE/AHH,ATE1,STTE,DTE
& /BLO10/BETA,ETA,DELTA,ZETA /ARGEXP/ARGMAX
EXTERNAL XE1,XE2,XE3,XE4,XE5,XE6,TEDER
DATA HOA /300.,400.,600./, XNAR /3*0.0/,
& XDELS /3*5.,10./, DNDS /.016,.01,2*.016/,
& DDO /9,5,5,25/, DO2 /5,5/,
& B0B1 /.755566,.778596,.797332,.812928,.826146/
data icalls/0/
C
SAVE EGNR,AGNR,DAYNR,DO2,SEASON,SEADAY,LATI,LONGI,MO2,
& MODIP,MAGBR,NMF1,NEI,MLAT,MLONG,NOBO2,
& F,RIF,FF0,XM0,F2,FM3,AMP,HXL,SCL,
& CTN,CTNN,ATE,TEA,HOA,XNAR,PG1O,
& PG2O,PG3O,PF1O,PF2O,PF3O,HO,MO,DDO,HO2,DION,EXT,
& SCHALT,SSIN,TCON,F1REG,FOF2IN,HMF2IN,URSIF2,LAYVER,
& DY,GULB0,NODEN,NOTEM,NOION,TENEOP,OLD79,TOPSI,BOTTO,BELOWE,
& URSIFO,MONTH,MONTHO,RG,RGO
C
C PROGAM CONSTANTS
C
icalls=icalls+1
HHALF = 0.
! write(6,"('Enter iri90: icalls=',i4,' jf=',12l2,' jmax=',i4,
! | ' alati=',f8.2,' along=',f8.2)")
! | icalls,jf,jmax,alati,along
! write(6,"('rz12=',f8.2,' mmdd=',i4,' dhour=',f8.2,' nz=',i4)")
! | rz12,mmdd,dhour,nz
! write(6,"('DIRECT (file path)=',a)") direct
! write(6,"('zkm=',/,(8f9.2))") zkm
ARGMAX=88.0
UMR=ATAN(1.0)*4./180.
ALOG2=ALOG(2.)
ALG100=ALOG(100.)
ISTART=1
heibeg=zkm(1)
heiend=zkm(nz)
C
C Code inserted to aleviate block data problem for PC version.
C Thus avoiding DATA statement with parameters from COMMON block.
C
AHH(1)=120.
AHH(2)=0.
AHH(3)=300.
AHH(4)=400.
AHH(5)=600.
AHH(6)=1400.
AHH(7)=3000.
DTE(1)=5.
DTE(2)=5.
DTE(3)=10.
DTE(4)=20.
DTE(5)=20.
DTI(1)=10.
DTI(2)=10.
DTI(3)=20.
DTI(4)=20.
C
C FIRST SPECIFY YOUR COMPUTERS CHANNEL NUMBERS ....................
C AGNR=OUTPUT (OUTPUT IS DISPLAYED OR STORED IN FILE OUTPUT.IRI)...
C IUCCIR=UNIT NUMBER FOR CCIR COEFFICIENTS ........................
C
MONITO=6
IUCCIR=10
KONSOL=6
IF (JF(12)) KONSOL=12
c
c selection of density and ion composition options ..................
c
NODEN=(.NOT.JF(1))
NOTEM=(.NOT.JF(2))
NOION=(.NOT.JF(3))
DY=(.NOT.JF(6))
LAYVER=(.NOT.JF(11))
OLD79=(.NOT.JF(7))
GULB0=(.NOT.JF(4))
c
c f peak density ....................................................
c
FOF2IN=(.NOT.JF(8))
IF(FOF2IN) THEN
AFOF2=OARR(1)
IF(AFOF2.GT.100.) AFOF2=SQRT(AFOF2/1.24E10)
ENDIF
URSIF2=(.NOT.JF(5))
c
c f peak altitude ..................................................
c
HMF2IN=(.NOT.JF(9))
IF(HMF2IN) AHMF2=OARR(2)
c
C TE-NE MODEL OPTION ..............................................
C
TENEOP=(.NOT.JF(10))
IF(TENEOP) THEN
DO 8154 JXNAR=1,3
XNAR(JXNAR)=OARR(JXNAR+2)
TCON(JXNAR)=.FALSE.
8154 IF(XNAR(JXNAR).GT.0.) TCON(JXNAR)=.TRUE.
ENDIF
if(icalls.gt.1) goto 8201
write(konsol,*) '*** IRI parameters are being calculated ***'
if(NODEN) goto 2889
if(LAYVER) write(konsol,*) 'Ne, E-F: The LAY-Version is ',
& 'prelimenary. Erroneous profile features can occur.'
if(GULB0) write(konsol,*) 'Ne, B0: Bottomside thickness is ',
& 'obtained with Gulyaeva-1987 model.'
if(OLD79) write(konsol,*) 'Ne: Using IRI-79. Correction',
& ' of equatorial topside is not included.'
if(HMF2IN) write(konsol,*) 'Ne, hmF2: Input values are used.'
if(FOF2IN) then
write(konsol,*) 'Ne, foF2: Input values are used.'
goto 2889
endif
if(URSIF2) then
write(konsol,*) 'Ne, foF2: URSI model is used.'
else
write(konsol,*) 'Ne, foF2: CCIR model is used.'
endif
2889 if((.not.NOION).and.(DY))
& write(konsol,*) 'Ion Com.: Using Danilov-Yaichnikov-1985.'
if((.not.NOTEM).and.(TENEOP))
& write(konsol,*) 'Te: Temperature-density correlation is used'
8201 continue
C
C CALCULATION OF MEAN F10.7CM SOLAR RADIO FLUX (COV)................
C CALCULATION OF RESTRICTED SOLAR ACTIVITIES (RG,COVG)..............
C
IF(RZ12.GE.0.0) THEN
R=RZ12
COV=63.75+R*(0.728+R*0.00089)
ELSE
COV=-RZ12
R=33.52*(SQRT(COV+85.12)-12.2)
ENDIF
RG=R
COVG=COV
IF(R.GT.150.) RG=150.
IF(COV.GT.193.) COVG=193.
C
C CALCULATION OF GEOG. OR GEOM. COORDINATES IN DEG....................
C CALCULATION OF MAGNETIC INCLINATION (DIP), DECLINATION (DEC)........
C DIP LATITUDE (MAGBR) AND MODIFIED DIP (MODIP). ALL IN DEGREE......
C
IF(JMAG.GT.0) THEN
MLAT=ALATI
MLONG=ALONG
if (mlong .lt. 0.) mlong=mlong+360.
ELSE
LATI=ALATI
LONGI=ALONG
if (longi .lt. 0.) longi=longi+360.
ENDIF
CALL GGM(JMAG,LONGI,LATI,MLONG,MLAT)
ABSLAT=ABS(LATI)
CALL FIELDG(LATI,LONGI,300.0,XMA,YMA,ZMA,BET,DIP,DEC,MODIP)
MAGBR=ATAN(0.5*TAN(DIP*UMR))/UMR
ABSMLT=ABS(MLAT)
ABSMDP=ABS(MODIP)
ABSMBR=ABS(MAGBR)
C
C CALCULATION OF SEASON (SUMMER=2, WINTER=4)..........................
C CALCULATION OF DAY OF YEAR AND SUN DECLINATION......................
C
if(MMDD.lt.0) then
DAYNR=-MMDD
call MODA(1,MONTH,IDAY,DAYNR)
else
MONTH=MMDD/100
IDAY=MMDD-MONTH*100
call MODA(0,MONTH,IDAY,DAYNR)
endif
SEASON=INT((DAYNR+45.0)/92.0)
IF(SEASON.LT.1) SEASON=4
NSESON=SEASON
seaday=daynr
IF(LATI.GT.0.0) GOTO 5592
SEASON=SEASON-2
IF(SEASON.LT.1) SEASON=SEASON+4
seaday=daynr+183
if(seaday.gt.366) seaday=seaday-366
C
C CALCULATION OF SOLAR ZENITH ANGLE (XHI/DEG).........................
C NOON VALUE (XHINON).................................................
C
5592 IF(DHOUR.GT.24.1) THEN
UT=DHOUR-25.
HOUR=UT+LONGI/15.
IF(HOUR.GT.24.) HOUR=HOUR-24.
ELSE
HOUR=DHOUR
UT=HOUR-LONGI/15.
IF(UT.LT.0.) UT=UT+24.
ENDIF
CALL SOCO(DAYNR,HOUR,LATI,LONGI,SUNDEC,XHI,SAX,SUX)
CALL SOCO(DAYNR,12.0,LATI,LONGI,SUNDE1,XHINON,SAXNON,SUXNON)
NIGHT=.FALSE.
if(abs(sax).gt.25.0) then
if(sax.lt.0.0) NIGHT=.TRUE.
goto 1334
endif
if(SAX.le.SUX) goto 1386
if((hour.gt.sux).and.(hour.lt.sax)) night=.true.
goto 1334
1386 IF((HOUR.GT.SUX).OR.(HOUR.LT.SAX)) NIGHT=.TRUE.
C
C CALCULATION OF ELECTRON DENSITY PARAMETERS................
C
1334 HNEA=65.
IF(NIGHT) HNEA=80.
HNEE=2000.
IF(NODEN) GOTO 4933
DELA=4.32
IF(ABSMDP.GE.18.) DELA=1.0+EXP(-(ABSMDP-30.0)/10.0)
DELL=1+EXP(-(ABSLAT-20.)/10.)
C!!!!!!! F-REGION PARAMETERS AND E-PEAK !!!!!!!!!!!!!!!!!!!!!!!!!!
FOE=FOEEDI(COV,XHI,XHINON,ABSLAT)
NME=1.24E10*FOE*FOE
HME=105.0
IF((FOF2IN).AND.(HMF2IN)) GOTO 501
IF (URSIF2 .NEQV. URSIFO) GOTO 7797
IF((MONTH.EQ.MONTHO).AND.(RG.EQ.RGO)) GOTO 4292
IF(MONTH.EQ.MONTHO) GOTO 4291
C
C READ CCIR COEFFICIENT SET FOR CHOSEN MONTH....................
C
7797 WRITE(filename,104) MONTH+10
104 FORMAT('ccir',I2,'.asc')
call dfp(direct,filename,path)
OPEN(IUCCIR,FILE=path,STATUS='OLD',ERR=8448)
READ(IUCCIR,4689) F2,FM3
4689 FORMAT(4E15.8)
CLOSE(IUCCIR)
C
C READ URSI COEFFICIENT SET FOR CHOSEN MONTH....................
C
if (URSIF2) then
WRITE(filename,1144) MONTH+10
1144 FORMAT('ursi',I2,'.asc')
call dfp(direct,filename,path)
OPEN(IUCCIR,FILE=path,STATUS='OLD',ERR=8448)
READ(IUCCIR,4689) F2
CLOSE(IUCCIR)
endif
URSIFO=URSIF2
MONTHO=MONTH
GOTO 4291
8448 write(monito,8449) path
8449 format(' IRI90: File ',A50,'not found')
stop
C
C LINEAR INTERPOLATION IN SOLAR ACTIVITY
C
4291 RR2=RG/100.
RR1=1.-RR2
DO 20 I=1,76
DO 20 J=1,13
K=J+13*(I-1)
20 FF0(K)=F2(J,I,1)*RR1+F2(J,I,2)*RR2
DO 30 I=1,49
DO 30 J=1,9
K=J+9*(I-1)
30 XM0(K)=FM3(J,I,1)*RR1+FM3(J,I,2)*RR2
RGO=RG
4292 CALL F2OUT(MODIP,LATI,LONGI,FF0,XM0,UT,YFOF2,XM3000)
501 IF(FOF2IN) THEN
FOF2=AFOF2
ELSE
FOF2=YFOF2
ENDIF
NMF2=1.24E10*FOF2*FOF2
IF(HMF2IN) THEN
HMF2=AHMF2
ELSE
HMF2=HMF2ED(MAGBR,RG,FOF2/FOE,XM3000)
ENDIF
TOPSI=(HEIEND.GT.HMF2)
BOTTO=((HEIEND.GE.HME).AND.(HEIBEG.LE.HMF2))
BELOWE=(HEIBEG.LT.HME)
c
c topside profile parameters .............................
c
IF(.NOT.TOPSI) GOTO 1501
COS2=COS(MLAT*UMR)
COS2=COS2*COS2
FLU=(COVG-40.0)/30.0
IF(OLD79) then
ETA1=-0.0070305*COS2
else
EX=EXP(-MLAT/15.)
EX1=EX+1
EPIN=4.*EX/(EX1*EX1)
ETA1=-0.02*EPIN
endif
ETA=0.058798+ETA1+FLU*(-0.014065+0.0069724*COS2)+
&(0.0024287+0.0042810*COS2-0.00015280*FOF2)*FOF2
ZETA=0.078922-0.0046702*COS2+FLU*(-0.019132+0.0076545*COS2)+
&(0.0032513+0.0060290*COS2-0.00020872*FOF2)*FOF2
BETA=-128.03+20.253*COS2+FLU*(-8.0755-0.65896*COS2)+(0.44041
&+0.71458*COS2-0.042966*FOF2)*FOF2
Z=EXP(94.45/BETA)
Z1=Z+1
Z2=Z/(BETA*Z1*Z1)
DELTA=(ETA/Z1-ZETA/2.0)/(ETA*Z2+ZETA/400.0)
c
c bottomside profile parameters .............................
C
1501 HMF1=HMF2
HZ=HMF2
HEF=HME
IF(.not.BOTTO) GOTO 2727
B1=3.0
C!!!!!!! INTERPOLATION FOR B0 OUT OF ARRAY B0F !!!!!!!!!!!!!!!!!!!!!
if(GULB0) then
call ROGUL(SEADAY,XHI,SEAX,GRAT)
if(NIGHT) GRAT=0.91-HMF2/4000.
B0CNEW=HMF2*(1.-GRAT)
B0=B0CNEW/B0B1(1)
else
B0 = B0POL(HOUR,SAX,SUX,SEASON,RG,DELA)
endif
C!!!!!!! F1-REGION PARAMETERS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
F1REG=.FALSE.
HMF1=0.
PNMF1=0.
C1=0.
IF(NIGHT.OR.(SEASON.EQ.4)) GOTO 150
FOF1=FOF1ED(ABSMBR,R,XHI)
IF(FOF1.LT.1.E-3) GOTO 150
F1REG=.TRUE.
C1=.09+.11/DELA
PNMF1=1.24E10*FOF1*FOF1
150 NMF1=PNMF1
C!!!!!!! PARAMETER FOR E AND VALLEY-REGION !!!!!!!!!!!!!!!!!!!!!
XDEL=XDELS(SEASON)/DELA
DNDHBR=DNDS(SEASON)/DELA
HDEEP=HPOL(HOUR,10.5/DELA,28.,SAX,SUX,1.,1.)
WIDTH=HPOL(HOUR,17.8/DELA,45.+22./DELA,SAX,SUX,1.,1.)
DEPTH=HPOL(HOUR,XDEL,81.,SAX,SUX,1.,1.)
DLNDH=HPOL(HOUR,DNDHBR,.06,SAX,SUX,1.,1.)
IF(DEPTH.LT.1.0) GOTO 600
IF(NIGHT) DEPTH=-DEPTH
CALL TAL(HDEEP,DEPTH,WIDTH,DLNDH,EXT,E)
IF(.NOT.EXT) GOTO 667
WRITE(KONSOL,650)
650 FORMAT(1X,'*NE* E-REGION VALLEY CAN NOT BE MODELLED')
600 WIDTH=.0
667 HEF=HME+WIDTH
VNER = (1. - ABS(DEPTH) / 100.) * NME
c
c Parameters below E .............................
c
2727 IF(.not.BELOWE) GOTO 2726
C!!!!!!!D-REGION PARAMETER!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
NMD=XMDED(XHI,R,4.0E8)
HMD=HPOL(HOUR,81.0,88.0,SAX,SUX,1.,1.)
F(1)=HPOL(HOUR,0.02+0.03/DELA,0.05,SAX,SUX,1.,1.)
F(2)=HPOL(HOUR,4.6,4.5,SAX,SUX,1.,1.)
F(3)=HPOL(HOUR,-11.5,-4.0,SAX,SUX,1.,1.)
FP1=F(1)
FP2=-FP1*FP1/2.0
FP30=(-F(2)*FP2-FP1+1.0/F(2))/(F(2)*F(2))
FP3U=(-F(3)*FP2-FP1-1.0/F(3))/(F(3)*F(3))
HDX=HMD+F(2)
X=HDX-HMD
XDX=NMD*EXP(X*(FP1+X*(FP2+X*FP30)))
DXDX=XDX*(FP1+X*(2.0*FP2+X*3.0*FP30))
X=HME-HDX
XKK=-DXDX*X/(XDX*ALOG(XDX/NME))
D1=DXDX/(XDX*XKK*X**(XKK-1.0))
C
C SEARCH FOR HMF1 ..................................................
C
2726 IF(.not.BOTTO) GOTO 4933
if(LAYVER) goto 6153
924 IF(.not.F1REG) GOTO 380
XE2H=XE2(HEF)
CALL REGFA1(HEF,HMF2,XE2H,NMF2,0.001,NMF1,XE2,SCHALT,HMF1)
IF(.not.SCHALT) GOTO 380
WRITE(KONSOL,11)
11 FORMAT(1X,'*NE* HMF1 IS NOT EVALUATED BY THE FUNCTION XE2')
IREGFA=1
c
c change B1 and try again ..........................................
c
9244 IF(B1.GT.4.5) GOTO (7398,8922) IREGFA
B1=B1+0.5
WRITE(KONSOL,902) B1-0.5,B1
902 FORMAT(6X,'CORR.: B1(OLD)=',F4.1,' B1(NEW)=',F4.1)
IF(GULB0) then
ib1=int(b1*2.-5.)
B0=B0CNEW/b0b1(ib1)
endif
GOTO 924
c
c omit F1 feature ....................................................
c
7398 WRITE(KONSOL,9269)
9269 FORMAT(1X,'CORR.: NO F1 REGION, B1=3, C1=0.0')
HMF1=0.
NMF1=0.
C1=0.0
B1=3.
F1REG=.FALSE.
C
C SEARCH FOR HST [NE3(HST)=NME] ..........................................
C
380 RRRR=0.5
IF(F1REG) then
hf1=hmf1
xf1=nmf1
GOTO 3972
ENDIF
RATHH=0.5
3973 hf1=hef+(hmf2-hef)*RATHH
xf1=xe3(hf1)
IF(XF1.LT.NME) THEN
RATHH=RATHH+.1
GOTO 3973
ENDIF
3972 h=hf1
deh=10.
XXMIN=XF1
HHMIN=HF1
3895 h=h-deh
if(h.lt.HEF) then
h=h+2*deh
deh=deh/10.
if(deh.lt.1.) goto 3885
endif
XE3H=XE3(h)
IF(XE3H.LT.XXMIN) then
XXMIN=XE3H
HHMIN=h
endif
if(XE3H.gt.NME) goto 3895
CALL REGFA1(h,HF1,XE3H,XF1,0.001,NME,XE3,SCHALT,HST)
STR=HST
IF(.not.SCHALT) GOTO 360
3885 WRITE(KONSOL,100)
100 FORMAT(1X,'*NE* HST IS NOT EVALUATED BY THE FUNCTION XE3')
IREGFA=2
IF(XXMIN/NME.LT.1.3) GOTO 9244
c
c assume linear interpolation between HZ and HEF ..................
c
8922 HZ=HHMIN+(HF1-HHMIN)*RRRR
XNEHZ=XE3(HZ)
if(xnehz-nme.lt.0.001) then
RRRR=RRRR+.1
GOTO 8922
endif
WRITE(KONSOL,901) HZ,HEF
901 FORMAT(6X,'CORR.: LIN. APP. BETWEEN HZ=',F5.1,
& ' AND HEF=',F5.1)
T=(XNEHZ-NME)/(HZ-HEF)
HST=-333.
GOTO 4933
c
c calculate HZ, D and T ............................................
c
360 HZ=(HST+HF1)/2.0
D=HZ-HST
T=D*D/(HZ-HEF-D)
GOTO 4933
C
C LAY-functions for middle ionosphere
C
6153 HMF1M=165.+0.6428*XHI
HHALF = GRAT * HMF2
HV1R = HME + WIDTH
HV2R = HME + HDEEP
HHMF2 = HMF2
CALL INILAY(NIGHT,NMF2,NMF1,NME,VNER,HHMF2,HMF1M,HME,
& HV1R,HV2R,HHALF,HXL,SCL,AMP,IIQU)
IF(IIQU.EQ.1) WRITE(KONSOL,7733)
7733 FORMAT('*NE* LAY amplitudes found with 2nd choice of HXL(1).')
IF(IIQU.EQ.2) WRITE(KONSOL,7722)
7722 FORMAT('*NE* LAY amplitudes could not be found.')
C---------- CALCULATION OF NEUTRAL TEMPERATURE PARAMETER-------
4933 HTA=120.0
HTE=3000.0
IF(NOTEM) GOTO 240
SEC=UT*3600.
CALL CIRA86(DAYNR,SEC,LATI,LONGI,HOUR,COV,TEXOS,TN120,SIGMA)
IF(HOUR.NE.0.0) THEN
SECNI=(24.-LONGI/15.)*3600.
CALL CIRA86(DAYNR,SECNI,LATI,LONGI,0.,COV,TEXNI,TN1NI,SIGNI)
ELSE
TEXNI=TEXOS
TN1NI=TN120
SIGNI=SIGMA
ENDIF
TLBDH=TEXOS-TN120
TLBDN=TEXNI-TN1NI
C
C--------- CALCULATION OF ELECTRON TEMPERATURE PARAMETER--------
C
881 CONTINUE
C !!!!!!!!!! TE(120KM)=TN(120KM) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ATE(1)=TN120
C !!!!!!!!!! TE-MAXIMUM (JICAMARCA,ARECIBO) !!!!!!!!!!!!!!!!!!!!
HMAXD=60.*EXP(-(MLAT/22.41)**2)+210.
HMAXN=150.
AHH(2)=HPOL(HOUR,HMAXD,HMAXN,SAX,SUX,1.,1.)
TMAXD=800.*EXP(-(MLAT/33.)**2)+1500.
TMAXN=TN(HMAXN,TEXNI,TLBDN,SIGNI)+20
ATE(2)=HPOL(HOUR,TMAXD,TMAXN,SAX,SUX,1.,1.)
C !!!!!!!!!! TE(300,400KM)=TE-AE-C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C !!!!!!!!!! TE(1400,3000KM)=TE-ISIS !!!!!!!!!!!!!!!!!!!!!!!!!!!
DIPLAT=MAGBR
CALL TEBA(DIPLAT,HOUR,NSESON,TEA)
ATE(3)=TEA(1)
ATE(4)=TEA(2)
ATE(6)=TEA(3)
ATE(7)=TEA(4)
C !!!!!!!!!! TE(600KM)=TE-AEROS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ETT=EXP(-MLAT/11.35)
TET=2900.-5600.*ETT/((ETT+1)**2.)
TEN=839.+1161./(1.+EXP(-(ABSMLT-45.)/5.))
ATE(5)=HPOL(HOUR,TET,TEN,SAX,SUX,1.5,1.5)
C !!!!!!!!!! OPTION TO USE TE-NE-RELATION !!!!!!!!!!!!!!!!!!!!!!
C !!!!!!!!!! AT 300, 400 OR 600 KM !!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF(TENEOP) THEN
DO 3395 I=1,3
3395 IF(TCON(I)) ATE(I+2)=TEDE(HOA(I),XNAR(I),-COV)
ENDIF
C !!!!!!!!!! TE'S ARE CORRECTED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C !!!!!!!!!! ALSO TE > TN ENFORCED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
TNAHH2=TN(AHH(2),TEXOS,TLBDH,SIGMA)
IF(ATE(2).LT.TNAHH2) ATE(2)=TNAHH2
STTE1=(ATE(2)-ATE(1))/(AHH(2)-AHH(1))
DO 1901 I=2,6
TNAHHI=TN(AHH(I+1),TEXOS,TLBDH,SIGMA)
IF(ATE(I+1).LT.TNAHHI) ATE(I+1)=TNAHHI
STTE2=(ATE(I+1)-ATE(I))/(AHH(I+1)-AHH(I))
ATE(I)=ATE(I)-(STTE2-STTE1)*DTE(I-1)*ALOG2
1901 STTE1=STTE2
C !!!!!!!!!! GRADIENTS ARE CALCULATED WITH !!!!!!!!!!!!!!!!!!!!
C !!!!!!!!!! CORRECTED REGION BOUNDARIES !!!!!!!!!!!!!!!!!!!!!!
DO 1902 I=1,6
1902 STTE(I)=(ATE(I+1)-ATE(I))/(AHH(I+1)-AHH(I))
ATE1=ATE(1)
887 CONTINUE
C
C------------ CALCULATION OF ION TEMPERATURE PARAMETERS--------
C
C !!!!!!!!!! TI(430KM,DAY)=TI-AEROS !!!!!!!!!!!!!!!!!!!!!!!!!!!
XSM1=430.0
XSM(1)=XSM1
Z1=EXP(-0.09*MLAT)
Z2=Z1+1.
TID1 = 1240.0 - 1400.0 * Z1 / ( Z2 * Z2 )
MM(2)=HPOL(HOUR,3.0,0.0,SAX,SUX,1.,1.)
C !!!!!!!!!! TI < TE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
TED1=TEA(6)+30.
IF(TID1.GT.TED1) TID1=TED1
C !!!!!!!!!! TI(430KM,NIGHT)=TI-AEROS !!!!!!!!!!!!!!!!!!!!!!!!!
Z1=ABSMLT
Z2=Z1*(0.47+Z1*0.024)*UMR
Z3=COS(Z2)
TIN1=1200.0-300.0*SIGN(1.0,Z3)*SQRT(ABS(Z3))
C !!!!!!!!!! TN < TI < TE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
TEN1=TEA(5)
TNN1=TN(XSM1,TEXNI,TLBDN,SIGNI)
IF(TEN1.LT.TNN1) TEN1=TNN1
IF(TIN1.GT.TEN1) TIN1=TEN1
IF(TIN1.LT.TNN1) TIN1=TNN1
C !!!!!!!!!! TI(430KM,LT) FROM STEP FUNCTION !!!!!!!!!!!!!!!!!!
TI1=TIN1
IF(TID1.GT.TIN1) TI1=HPOL(HOUR,TID1,TIN1,SAX,SUX,1.,1.)
C !!!!!!!!!! TANGENT ON TN DETERMINES HS !!!!!!!!!!!!!!!!!!!!!!
TI13=TEDER(130.)
TI50=TEDER(500.)
CALL REGFA1(130.0,500.0,TI13,TI50,0.01,TI1,TEDER,SCHALT,HS)
IF(SCHALT) HS=200.
TNHS=TN(HS,TEXOS,TLBDH,SIGMA)
MM(1)=DTNDH(HS,TEXOS,TLBDH,SIGMA)
IF(SCHALT) MM(1)=(TI1-TNHS)/(XSM1-HS)
MXSM=2
C !!!!!!!!!! XTETI ALTITTUDE WHERE TE=TI !!!!!!!!!!!!!!!!!!!!!!
2391 XTTS=500.
X=500.
2390 X=X+XTTS
IF(X.GE.AHH(7)) GOTO 240
TEX=ELTE(X)
TIX=TI(X)
IF(TIX.LT.TEX) GOTO 2390
X=X-XTTS
XTTS=XTTS/10.
IF(XTTS.GT.0.1) GOTO 2390
XTETI=X+XTTS*5.
C !!!!!!!!!! TI=TE ABOVE XTETI !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
MXSM=3
MM(3)=STTE(6)
XSM(2)=XTETI
IF(XTETI.GT.AHH(6)) GOTO 240
MXSM=4
MM(3)=STTE(5)
MM(4)=STTE(6)
XSM(3)=AHH(6)
IF(XTETI.GT.AHH(5)) GOTO 240
MXSM=5
DTI(1)=5.
DTI(2)=5.
MM(3)=STTE(4)
MM(4)=STTE(5)
MM(5)=STTE(6)
XSM(3)=AHH(5)
XSM(4)=AHH(6)
C
C CALCULATION OF ION DENSITY PARAMETER..................
C
240 IF(NOION) GOTO 141
HNIA=100.
HNIE=2000.
if(DY) goto 141
C
C INPUT OF THE ION DENSITY PARAMETER ARRAYS PF1O,PF2O AND PF3O......
C
RIF(1)=2.
IF(ABSLAT.LT.30.0) RIF(1)=1.
RIF(2)=2.
IF(COV.LT.100.0) RIF(2)=1.
RIF(3)=SEASON
IF(SEASON.EQ.1) RIF(3)=3.
RIF(4)=1.
IF(NIGHT) RIF(4)=2.
CALL KOEFP1(PG1O)
CALL KOEFP2(PG2O)
CALL KOEFP3(PG3O)
CALL SUFE(PG1O,RIF,12,PF1O)
CALL SUFE(PG2O,RIF, 4,PF2O)
CALL SUFE(PG3O,RIF,12,PF3O)
c
c calculate O+ profile parameters
c
IF(ABS(XHI).LE.90.0) THEN
ZZZ1=COS(XHI*UMR)
ELSE
ZZZ1=0.0
ENDIF
msumo=4
RDOMAX=100.0
MO(1)=EPSTEP(PF1O(1),PF1O(2),PF1O(3),PF1O(4),ZZZ1)
MO(2)=EPSTEP(PF1O(5),PF1O(6),PF1O(7),PF1O(8),ZZZ1)
MO(3)=0.0
HO(1)=EPSTEP(PF1O(9),PF1O(10),PF1O(11),PF1O(12),ZZZ1)
HO(2)=290.0
IF((RIF(2).EQ.2.).AND.(RIF(3).EQ.2.)) HO(2)=237.0
HO(4)=PF2O(1)
ho05=pf2o(4)
MO(4)=PF2O(2)
MO(5)=PF2O(3)
c
c adjust gradient MO(4) of O+ profile segment above F peak
c
7100 HO(3)=(ALG100-MO(5)*(HO(4)-ho05))/MO(4)+HO(4)
IF(HO(3).LE.HO(2)+20.) THEN
MO(4)=MO(4)-0.001
GOTO 7100
endif
hfixo=(ho(2)+ho(3))/2.
c
c find height H0O of maximum O+ relative density
c
DELX=5.0
X=HO(2)
YMAXX=0.0
7102 X=X+DELX
Y=RPID(X,HFIXO,RDOMAX,msumo,MO,DDO,HO)
IF(Y.LE.YMAXX) then
if(delx.le.0.1) GOTO 7104
x=x-delx
delx=delx/5.
ELSE
YMAXX=Y
ENDIF
GOTO 7102
7104 H0O=X-DELX/2.
7101 if(y.lt.100.0) goto 7103
rdomax=rdomax-0.01
y=rpid(h0o,hfixo,rdomax,msumo,mo,ddo,ho)
goto 7101
7103 yo2h0o=100.-y
yoh0o=y
c
c calculate parameters for O2+ profile
c
hfixo2 = pf3o(1)
rdo2mx = pf3o(2)
DO 7105 L=1,2
I = L * 2
HO2(L)=PF3O(1+I)+PF3O(2+I)*ZZZ1
7105 MO2(L+1)=PF3O(7+I)+PF3O(8+I)*ZZZ1
MO2(1)=PF3O(7)+PF3O(8)*ZZZ1
if(hfixo2.gt.ho2(1)) then
ymo2z=mo2(2)
else
ymo2z=mo2(1)
endif
aldo21=alog(rdo2mx)+ymo2z*(ho2(1)-hfixo2)
hfixo2=(ho2(2)+ho2(1))/2.
rdo2mx=exp(aldo21+mo2(2)*(hfixo2-ho2(1)))
c
c make sure that rd(O2+) is less or equal 100-rd(O+) at O+ maximum
c
7106 Y=RPID(H0O,hfixo2,rdo2mx,2,MO2,DO2,HO2)
IF(Y.GT.yo2h0o) then
MO2(3)=MO2(3)-0.02
GOTO 7106
endif
c
C use ratio of NO+ to O2+ density at O+ maximum to calculate
c NO+ density above the O+ maximum (H0O)
c
IF(y.LT.1.) then
NOBO2=0.0
ELSE
NOBO2= (yo2h0o-y)/y
ENDIF
C
C CALCULATION FOR THE REQUIRED HEIGHT RANGE.......................
C
141 IF(.NOT.F1REG) HMF1=HZ
DO 7397 KI=1,11
DO 7397 KK=1,nz
7397 OUTF(KI,KK)=-1.
do 7118 kk=1,nz
height=zkm(kk)
300 IF(NODEN) GOTO 330
IF((HEIGHT.GT.HNEE).OR.(HEIGHT.LT.HNEA)) GOTO 330
IF(LAYVER) THEN
ELEDE=-9.
IF(IIQU.LT.2) ELEDE=XEN(HEIGHT,HMF2,NMF2,HME,4,HXL,SCL,AMP)
ELSE
ELEDE=XE(HEIGHT)
ENDIF
OUTF(1,KK)=ELEDE
330 IF(NOTEM) GOTO 7108
IF((HEIGHT.GT.HTE).OR.(HEIGHT.LT.HTA)) GOTO 7108
TNH=TN(HEIGHT,TEXOS,TLBDH,SIGMA)
TIH=TNH
IF(HEIGHT.GE.HS) TIH=TI(HEIGHT)
TEH=ELTE(HEIGHT)
IF(TIH.LT.TNH) TIH=TNH
IF(TEH.LT.TIH) TEH=TIH
OUTF(2,KK)=TNH
OUTF(3,KK)=TIH
OUTF(4,KK)=TEH
7108 IF(NOION) GOTO 7118
IF((HEIGHT.GT.HNIE).OR.(HEIGHT.LT.HNIA)) GOTO 7118
if(DY) then
zmonth = month + iday / 30.0
call IONCOM(HEIGHT,XHI*UMR,LATI*UMR,COV,ZMONTH,DION)
ROX=DION(1)
RHX=DION(2)
RNX=DION(3)
RHEX=DION(4)
RNOX=DION(5)
RO2X=DION(6)
RCLUST=DION(7)
else
ROX=RPID(HEIGHT,HFIXO,RDOMAX,msumo,MO,DDO,HO)
RO2X=RPID(HEIGHT,HFIXO2,rdo2mx,2,MO2,DO2,HO2)
CALL RDHHE(HEIGHT,H0O,ROX,RO2X,NOBO2,10.,RHX,RHEX)
RNOX=RDNO(HEIGHT,H0O,RO2X,ROX,NOBO2)
RNX=-1.
RCLUST=-1.
endif
OUTF(5,KK)=ROX
OUTF(6,KK)=RHX
OUTF(7,KK)=RHEX
OUTF(8,KK)=RO2X
OUTF(9,KK)=RNOX
OUTF(10,KK)=RNX
OUTF(11,KK)=RCLUST
7118 continue
C
C ADDITIONAL PARAMETER FIELD OARR
C
IF(NODEN) GOTO 6192
OARR(1)=NMF2
OARR(2)=HMF2
OARR(3)=NMF1
OARR(4)=HMF1
OARR(5)=NME