CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC Empirical vertical drifts model for storm time; CC After Fejer and Scherliess, JGR, 102, 24047-24056,1997; CC The equatorial vertical drifts are described in; CC SCHERLIESS AND FEJER, JGR, 104, 6829-6842, 1999; CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PROGRAM MAIN PARAMETER (MAX=366*24*4) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*1 c1 CHARACTER*5 c5 CHARACTER*20 FLname INTEGER FLAG,jF,jL,yymmdd(1:MAX),YMD(1:MAX) REAL*8 UT(1:MAX),AE(1:MAX),time,long,param(2) REAL*8 F107(1:MAX) CC do i=1,max YYMMDD(i)=0.0D0 YMD(i)=0.0D0 AE(i)=0.0D0 UT(i)=0.0D0 F107(i)=0.0D0 end do CC WRITE(*,*) 'Input AE data file name:' READ(*,*) FLname CC c open(20,file=FLname,mode='read') open(20,file=FLname) READ(20,*,END=30) DO I=1,MAX READ(20,*,END=30) yymmdd(I),UT(I),AE(I) jL=I ENDDO 30 CLOSE(20) CC IF ((UT(jL)-UT(jL-1)).GT.0.5D0.AND.(UT(2)-UT(1)).GT.0.5D0) THEN FLAG=1 !1h resolution; jF=49 ELSE FLAG=-1 !15 min resolution; jF=193 ENDIF C Read F107; WRITE(*,*) 'Input F107 data file name:' READ(*,*) FLname c open(30,file=FLname,mode='read') open(30,file=FLname) DO I=1,MAX C2000-07-20 00:00 235.0 "" "" READ(30,35,END=40) iYEAR,c1,MONTH,c1,iDAY,c5,F107(I) YMD(I)=iYEAR*10000+MONTH*100+iDAY ENDDO 35 FORMAT(I4,A1,I2,A1,I2,A7,F5.1) 40 CLOSE(30) CC WRITE(*,*) 'INPUT GLON:' READ(*,*) long ! Geographic Longitude write(*,100) 'SLT','Drift (m/s)','Geogr. Long.','DoY','F10.7cm' CC c OPEN(10,file='drift.dat',MODE='WRITE') OPEN(10,file='drift.dat') DO iP=jF,jL time=UT(iP)+long/15.0D0 IF (time.GT.24.0D0) time=time-24.0D0 C Year, Month and Day; iYEAR=yymmdd(iP)/10000 MONTH=(yymmdd(iP)-iYEAR*10000)/100 iDAY=yymmdd(iP)-(yymmdd(iP)/100)*100 CALL MODA(0,iYEAR,MONTH,iDAY,iYD) C param(1)=iYD ! Day of year (e.g.,January 1st: param(1)=1.) param(2)=0.0 ! F10.7cm solar flux ij=1 45 IF (yymmdd(iP).EQ.YMD(ij)) THEN param(2)=F107(ij) !F10.7cm solar flux; GOTO 50 ELSE ij=ij+1 IF (ij.LT.MAX) GOTO 45 STOP 'F107 data error!' END IF C Vertical drifts during quiet time; 50 CALL vdrift_model(time,long,param,Vd_mean) C Vertical drifts during storm times; CALL StormVd(FLAG,iP,AE,time,PromptVd,DynamoVd,Vd_storm) VdTot=Vd_mean+Vd_storm write(10,200) yymmdd(iP),UT(iP),Vd_mean,Vd_storm,VdTot write(*, 300) yymmdd(iP),UT(iP),Vd_mean,param,Vd_storm,VdTot END DO CLOSE(10) 100 format(5x,a3,6x,a11,2x,a11,6x,a3,8x,a7) 200 format(1X,I8,F8.3,3F10.3) 300 format(1X,I8,F8.3,5F10.3) END C ******************************************************************* SUBROUTINE StormVd(FLAG,iP,AE,SLT,PromptVd,DynamoVd,Vd) C ******************************************************************* C Empirical vertical disturbance drifts model C After Fejer and Scherliess, JGR, 102, 24047-24056,1997 C********************************************************************* C INPUT: C AE: AE(in nT) in 1 hour or 15 minute resolution; C SLT: Local time(in hrs) for wanted Vd; C OUTPUT: C PromptVd: Prompt penetration vertical drifts at given conditions; C DynamoVd: Disturbane dynamo vertical drifts at given conditions; C Vd: PromptVd+DynamoVd; C********************************************************************* IMPLICIT REAL*8(A-H,O-Z) REAL*8 AE(1:366*24*4),Coff1(1:5,1:9),Coff15(1:6,1:9) INTEGER FLAG DATA Coff1/ @ 0.0124,-0.0168,-0.0152,-0.0174,-0.0704, @ -0.0090,-0.0022,-0.0107, 0.0152,-0.0674, @ 0.0275, 0.0051,-0.0132, 0.0020,-0.0110, @ -0.0022, 0.0044, 0.0095, 0.0036,-0.0206, @ 0.0162, 0.0007, 0.0085,-0.0140, 0.0583, @ 0.0181, 0.0185,-0.0109,-0.0031,-0.0427, @ -0.0057, 0.0002, 0.0086, 0.0149, 0.2637, @ -0.0193, 0.0035, 0.0117, 0.0099, 0.3002, @ -0.0492,-0.0201, 0.0338, 0.0099, 0.0746/ DATA Coff15/ @ 0.0177, 0.0118,-0.0006,-0.0152,-0.0174,-0.0704, @ 0.0051,-0.0074,-0.0096,-0.0107, 0.0152,-0.0674, @ 0.0241, 0.0183, 0.0122,-0.0132, 0.0020,-0.0110, @ 0.0019,-0.0010, 0.0001, 0.0095, 0.0036,-0.0206, @ 0.0170, 0.0183, 0.0042, 0.0085,-0.0140, 0.0583, @ 0.0086, 0.0189, 0.0200,-0.0109,-0.0031,-0.0427, @ -0.0070,-0.0053,-0.0090, 0.0086, 0.0149, 0.2637, @ -0.0326,-0.0101, 0.0076, 0.0117, 0.0099, 0.3002, @ -0.0470,-0.0455,-0.0274, 0.0338, 0.0099, 0.0746/ CCCCCCCCCCCCCCCCC**Define to variables**CCCCCCCCCCCCCCCCCCCCC C To 1 h time resolution: C dAEt_30=AE(t)-AE(t-1 hour); C dAEt_90=AE(t-1 hour)-AE(t-2 hour); CC C To 15 MIN time resolution : C dAEt_7P5=AE(t)-AE(t-15min); C dAEt_30=AE(t-15)-AE(t-45min); C dAEt_75=AE(t-45)-AE(t-105min); CC C Following variables are the same to two resolution: C AE1_6=average(AE(1-6hours)); C AE7_12=average(AE(7-12hours)); C AE1_12=average(AE(1-12hours)); C AEd1_6=average(X(AE(1-6hours)-130 nT)); C AEd7_12=average(X(AE(7-12hours)-130 nT)); C AEd1_12=average(X(AE(1-12hours)-130 nT)); C AEd22_28=average(X(AE(22-28hours)-130 nT)); C Here X(a)=a, a>0; =0, a<=0; C Alfa=0, AE1_6<200 nT; C AE1_6/100-2, 200 nT300 nT; C Beta=exp(-AE1_12/90), AE1_12>=70nT; C 0.46, AE1_12<70 nT; CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCccccccc C***************************************************** CC FLAG>0--> 1 h time resolution C**************************************************** IF (FLAG.GT.0) THEN C dAEt_30=AE(iP)-AE(iP-1) dAEt_90=AE(iP-1)-AE(iP-2) C AE1_6=0.0D0 AEd1_6=0.0D0 DO i=-1,-6,-1 AE1_6=AE1_6+AE(iP+i) AEd1_6S=AE(iP+i)-130.0D0 IF (AEd1_6S.LE.0.0D0) AEd1_6S=0.0D0 AEd1_6=AEd1_6+AEd1_6S END DO AE1_6=AE1_6/6.0D0 AEd1_6=AEd1_6/6.0D0 C AEd7_12=0.0D0 DO i=-7,-12,-1 AEd7_12S=AE(iP+i)-130.0D0 IF (AEd7_12S.LE.0.0D0) AE7_12S=0.0D0 AEd7_12=AEd7_12+AEd7_12S END DO AEd7_12=AEd7_12/6.0D0 C AE1_12=0.0D0 DO i=-1,-12,-1 AE1_12=AE1_12+AE(iP+i) END DO AE1_12=AE1_12/12.0D0 C AEd22_28=0.0D0 DO i=-22,-28,-1 AEd22_28S=AE(iP+i)-130.0D0 IF (AED22_28S.LE.0.0D0) AEd22_28S=0.0D0 AEd22_28=AEd22_28+AEd22_28S END DO AEd22_28=AEd22_28/7.0D0 AEd22_28P=AEd22_28-200.0D0 IF (AEd22_28P.LE.0.0D0) AEd22_28P=0.0D0 CC IF (AE1_6.GT.300.0D0) THEN Alfa=1.0D0 ELSE IF (AE1_6.GT.200.0D0) THEN ALfa=AE1_6/100.0D0-2.0D0 ELSE ALfa=0.0D0 ENDIF CC IF (AE1_12.GE.70.0D0) THEN Beta=dexp(-AE1_12/90.0D0) ELSE Beta=0.46D0 END IF PromptVd=0.0D0 DO J=1,9 PromptVd=PromptVd +(Coff1(1,J)*dAEt_30 +Coff1(2,J)*dAEt_90 # )*bspl4_ptime(J,SLT) END DO DynamoVd=0.0D0 DO J=1,9 DynamoVd=DynamoVd+ # (Coff1(3,J)*AEd1_6+Coff1(4,J)*Alfa*AEd7_12 # +Coff1(5,J)*Beta*AEd22_28P)*bspl4_ptime(J,SLT) END DO Vd=PromptVd+DynamoVd RETURN CC C 1 h time resolution end; C******************************************************************** C 15 min time resolution C******************************************************************** ELSE dAEt_7P5=AE(iP)-AE(iP-1) dAEt_30=AE(iP-1)-AE(iP-3) dAEt_75=AE(iP-3)-AE(iP-7) CC AE1_6=0.0D0 AEd1_6=0.0D0 DO i=-4,-24,-1 AE1_6=AE1_6+AE(iP+i) AEd1_6s=AE(iP+i)-130. IF (AEd1_6s.LE.0.0) AEd1_6s=0.0 AEd1_6=AEd1_6+AEd1_6S ENDDO AE1_6=AE1_6/21.0D0 AEd1_6=AEd1_6/21.0D0 CC AEd7_12=0.0D0 DO i=-28,-48,-1 AEd7_12s=AE(iP+i)-130.0 IF (AEd7_12s.LE.0) AEd7_12s=0.0 AEd7_12=AEd7_12+AEd7_12S ENDDO AEd7_12=AEd7_12/21.0D0 CC AE1_12=0.0D0 DO i=-4,-48,-1 AE1_12=AE1_12+AE(iP+i) END DO AE1_12=AE1_12/45.0D0 CC AEd22_28=0.0D0 DO i=-88,-112,-1 AEd22_28s=AE(iP+i)-130. IF (AEd22_28s.LE.0) AEd22_28s=0.0 AEd22_28=AEd22_28+AEd22_28s ENDDO AEd22_28=AEd22_28/25.0D0 AEd22_28P=AEd22_28-200.0D0 IF (AEd22_28P.LE.0.0D0) AEd22_28P=0.0D0 c AE1_6=0.0D0 c AEd1_6=0.0D0 c AEd7_12=0.0D0 c AEd22_28P=0.0D0 c AE1_12=0.0D0 c dAEt_7P5=400.D0 c dAEt_30=0.D0 c dAEt_75=0.D0 CC IF (AE1_6.GT.300.0D0) THEN Alfa=1.0D0 ELSE IF (AE1_6.GT.200.0D0) THEN ALfa=AE1_6/100.0D0-2.0D0 ELSE ALfa=0.0D0 ENDIF CC IF (AE1_12.GE.70.0D0) THEN Beta=dexp(-AE1_12/90.0D0) ELSE Beta=0.46D0 END IF CC PromptVd=0.0D0 DO J=1,9 PromptVd=PromptVd+(Coff15(1,J)*dAEt_7P5+Coff15(2,J)*dAEt_30 # +Coff15(3,J)*dAEt_75)*bspl4_ptime(J,SLT) END DO DynamoVd=0.0D0 print*,AEd1_6,AEd7_12,AEd22_28P,Alfa,Beta DO J=1,9 DynamoVd=DynamoVd +(Coff15(4,J)*AEd1_6+ # Coff15(5,J)*Alfa*AEd7_12+ # Coff15(6,J)*Beta*AEd22_28P # )*bspl4_ptime(J,SLT) END DO Vd=PromptVd+DynamoVd ENDIF RETURN END C ************************************************* real*8 function bspl4_ptime(i,x1) C ************************************************* IMPLICIT REAL*8 (A-H,O-Z) integer i,order/4/,j,k real*8 t_t(0:27) real*8 x,b(20,20),x1 data t_t/0.00,3.00,4.50,6.00,9.00,12.0,15.0,18.0,21.0, * 24.0,27.0,28.5,30.0,33.0,36.0,39.0,42.0,45.0, * 48.0,51.0,52.5,54.0,57.0,60.0,63.0,66.0,69.0,72.0/ C x=x1 if(i.ge.0) then if (x.lt.t_t(i-0)) then x=x+24 end if end if do j=i,i+order-1 if(x.ge.t_t(j).and.x.lt.t_t(j+1)) then b(j,1)=1 else b(j,1)=0 end if end do c do j=2,order do k=i,i+order-j b(k,j)=(x-t_t(k))/(t_t(k+j-1)-t_t(k))*b(k,j-1) b(k,j)=b(k,j)+(t_t(k+j)-x)/(t_t(k+j)-t_t(k+1))*b(k+1,j-1) end do end do bspl4_ptime=b(i,order) return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC CC SUBROUTINE MODA(IN,YEAR,MONTH,IDAY,IDOY) C------------------------------------------------------------------- C CALCULATES DAY OF YEAR (IDOY) FROM MONTH (MONTH) AND DAY (IDAY) C IF IN=0, OR MONTH (MONTH) AND DAY (IDAY) FROM DAY OF C YEAR (IDOY), IF IN>0. C------------------------------------------------------------------- IMPLICIT INTEGER (A-Z) DIMENSION MO(12),MO28(1:12),MO29(1:12) DATA MO28/0,31,59,90,120,151,181,212,243,273,304,334/ DATA MO29/0,31,60,91,121,152,182,213,244,274,305,335/ IMO=0 MOBE=0 IF (mod(YEAR,4).NE.0) THEN do i=1,12 MO(i)=MO28(i) end do ELSE IF (mod(YEAR,100).NE.0) THEN do i=1,12 MO(i)=MO29(i) end do ELSE IF (mod(YEAR,400).NE.0) THEN do i=1,12 MO(i)=MO28(i) end do ELSE do i=1,12 MO(i)=MO29(i) end do END IF IF (IN.GT.0) GOTO 5 IDOY=MO(MONTH)+IDAY RETURN 5 IMO=IMO+1 MOOLD=MOBE IF(IMO.GT.12) GOTO 55 MOBE=MO(IMO) IF(MOBE.LT.IDOY) GOTO 5 55 MONTH=IMO-1 IDAY=IDOY-MOOLD RETURN END C C ************************************************************ subroutine vdrift_model(xt,xl,param,y) C SUBROUTINE CALCULATES EQUATORIAL VERTICAL DRIFT AS DESCRIBED C IN SCHERLIESS AND FEJER, JGR, 104, 6829-6842, 1999 C ************************************************************ C INPUT: XT: SOLAR LOCAL TIME C XL: GEOGRAPHIC LONGITUDE (+ EAST) C PARAM: 2-DIM ARRAY (DOY,F10.7CM) C DOY :Day of Year has to run from 1 to 365 (366) C F10.7cm : F10.7cm solar flux C OUTPUT: Y: EQUATORIAL VERTICAL DRIFT C ************************************************************ IMPLICIT REAL*8 (A-H,O-Z) real*8 param(2),coeff(624),funct(6) real*8 xt,xl,y real*8 bspl4,bspl4_time,bspl4_long integer i,j,ind,il,kk integer index_t/13/ !,dim_t/78/ integer index_l/8/ !,dim_l/48/ C integer index/104/,dim/624/ integer nfunc/6/ data coeff/ * -10.80592, -9.63722,-11.52666, -0.05716, -0.06288, 0.03564, * -5.80962, -7.86988, -8.50888, -0.05194, -0.05798, -0.00138, * 2.09876,-19.99896, -5.11393, -0.05370, -0.06585, 0.03171, * -10.22653, -3.62499,-14.85924, -0.04023, -0.01190, -0.09656, * -4.85180,-26.26264, -6.20501, -0.05342, -0.05174, 0.02419, * -13.98936,-18.10416, -9.30503, -0.01969, -0.03132, -0.01984, * -18.36633,-24.44898,-16.69001, 0.02033, -0.03414, -0.02062, * -20.27621,-16.95623,-36.58234, 0.01445, -0.02044, -0.08297, * 1.44450, 5.53004, 4.55166, -0.02356, -0.04267, 0.05023, * 5.50589, 7.05381, 1.94387, -0.03147, -0.03548, 0.01166, * 3.24165, 10.05002, 4.26218, -0.03419, -0.02651, 0.07456, * 7.02218, 0.06708,-11.31012, -0.03252, -0.01021, -0.09008, * -3.47588, -2.82534, -4.17668, -0.03719, -0.01519, 0.06507, * -4.02607,-11.19563,-10.52923, -0.00592, -0.01286, -0.00477, * -11.47478, -9.57758,-10.36887, 0.04555, -0.02249, 0.00528, * -14.19283, 7.86422, -8.76821, 0.05758, -0.02398, -0.04075, * 14.58890, 36.63322, 27.57497, 0.01358, -0.02316, 0.04723, * 12.53122, 29.38367, 21.40356, -0.00071, -0.00553, 0.01484, * 18.64421, 26.27327, 18.32704, 0.00578, 0.03349, 0.11249, * 4.53014, 6.15099, 7.41935, -0.02860, -0.00395, -0.08394, * 14.29422, 9.77569, 2.85689, -0.00107, 0.04263, 0.10739, * 7.17246, 4.40242, -1.00794, 0.00089, 0.01436, 0.00626, * 7.75487, 5.01928, 4.36908, 0.03952, -0.00614, 0.03039, * 10.25556, 8.82631, 24.21745, 0.05492, -0.02968, 0.00177, * 21.86648, 24.03218, 39.82008, 0.00490, -0.01281, -0.01715, * 19.18547, 23.97403, 34.44242, 0.01978, 0.01564, -0.02434, * 26.30614, 14.22662, 31.16844, 0.06495, 0.19590, 0.05631, * 21.09354, 25.56253, 29.91629, -0.04397, -0.08079, -0.07903, * 28.30202, 16.80567, 38.63945, 0.05864, 0.16407, 0.07622, * 22.68528, 25.91119, 40.45979, -0.03185, -0.01039, -0.01206, * 31.98703, 24.46271, 38.13028, -0.08738, -0.00280, 0.01322, * 46.67387, 16.80171, 22.77190, -0.13643, -0.05277, -0.01982, * 13.87476, 20.52521, 5.22899, 0.00485, -0.04357, 0.09970, * 21.46928, 13.55871, 10.23772, -0.04457, 0.01307, 0.06589, * 16.18181, 16.02960, 9.28661, -0.01225, 0.14623, -0.01570, * 18.16289, -1.58230, 14.54986, -0.00375, -0.00087, 0.04991, * 10.00292, 11.82653, 0.44417, -0.00768, 0.15940, -0.01775, * 12.15362, 5.65843, -1.94855, -0.00689, 0.03851, 0.04851, * -1.25167, 9.05439, 0.74164, 0.01065, 0.03153, 0.02433, * -15.46799, 18.23132, 27.45320, 0.00899, -0.00017, 0.03385, * 2.70396, -0.87077, 6.11476, -0.00081, 0.05167, -0.08932, * 3.21321, -1.06622, 5.43623, 0.01942, 0.05449, -0.03084, * 17.79267, -3.44694, 7.10702, 0.04734, -0.00945, 0.11516, * 0.46435, 6.78467, 4.27231, -0.02122, 0.10922, -0.03331, * 15.31708, 1.70927, 7.99584, 0.07462, 0.07515, 0.08934, * 4.19893, 6.01231, 8.04861, 0.04023, 0.14767, -0.04308, * 9.97541, 5.99412, 5.93588, 0.06611, 0.12144, -0.02124, * 13.02837, 10.29950, -4.86200, 0.04521, 0.10715, -0.05465, * 5.26779, 7.09019, 1.76617, 0.09339, 0.22256, 0.09222, * 9.17810, 5.27558, 5.45022, 0.14749, 0.11616, 0.10418, * 9.26391, 4.19982, 12.66250, 0.11334, 0.02532, 0.18919, * 13.18695, 6.06564, 11.87835, 0.26347, 0.02858, 0.14801, * 10.08476, 6.14899, 17.62618, 0.09331, 0.08832, 0.28208, * 10.75302, 7.09244, 13.90643, 0.09556, 0.16652, 0.22751, * 6.70338, 11.97698, 18.51413, 0.15873, 0.18936, 0.15705, * 5.68102, 23.81606, 20.65174, 0.19930, 0.15645, 0.08151, * 29.61644, 5.49433, 48.90934, 0.70710, 0.40791, 0.26325, * 17.11994, 19.65380, 44.88810, 0.45510, 0.41689, 0.22398, * 8.45700, 34.54442, 27.25364, 0.40867, 0.37223, 0.22374, * -2.30305, 32.00660, 47.75799, 0.02178, 0.43626, 0.30187, * 8.98134, 33.01820, 33.09674, 0.33703, 0.33242, 0.41156, * 14.27619, 20.70858, 50.10005, 0.30115, 0.32570, 0.45061, * 14.44685, 16.14272, 45.40065, 0.37552, 0.31419, 0.30129, * 6.19718, 18.89559, 28.24927, 0.08864, 0.41627, 0.19993, * 7.70847, -2.36281,-21.41381, 0.13766, 0.05113, -0.11631, * -9.07236, 3.76797,-20.49962, 0.03343, 0.08630, 0.00188, * -8.58113, 5.06009, -6.23262, 0.04967, 0.03334, 0.24214, * -27.85742, 8.34615,-27.72532, -0.08935, 0.15905, -0.03655, * 2.77234, 0.14626, -4.01786, 0.22338, -0.04478, 0.18650, * 5.61364, -3.82235,-16.72282, 0.26456, -0.03119, -0.08376, * 13.35847, -6.11518,-16.50327, 0.28957, -0.01345, -0.19223, * -5.37290, -0.09562,-27.27889, 0.00266, 0.22823, -0.35585, * -15.29676,-18.36622,-24.62948, -0.31299, -0.23832, -0.08463, * -23.37099,-13.69954,-26.71177, -0.19654, -0.18522, -0.20679, * -26.33762,-15.96657,-42.51953, -0.13575, -0.00329, -0.28355, * -25.42140,-14.14291,-21.91748, -0.20960, -0.19176, -0.32593, * -23.36042,-23.89895,-46.05270, -0.10336, 0.03030, -0.21839, * -19.46259,-21.27918,-32.38143, -0.17673, -0.15484, -0.11226, * -19.06169,-21.13240,-34.01677, -0.25497, -0.16878, -0.11004, * -18.39463,-16.11516,-19.55804, -0.19834, -0.23271, -0.25699, * -19.93482,-17.56433,-18.58818, 0.06508, -0.18075, 0.02796, * -23.64078,-18.77269,-22.77715, -0.02456, -0.12238, 0.02959, * -12.44508,-21.06941,-19.36011, 0.02746, -0.16329, 0.19792, * -26.34187,-19.78854,-24.06651, -0.07299, -0.03082, -0.03535, * -10.71667,-26.04401,-16.59048, 0.02850, -0.09680, 0.15143, * -18.40481,-23.37770,-16.31450, -0.03989, -0.00729, -0.01688, * -9.68886,-20.59304,-18.46657, 0.01092, -0.07901, 0.03422, * -0.06685,-19.24590,-29.35494, 0.12265, -0.24792, 0.05978, * -15.32341, -9.07320,-13.76101, -0.17018, -0.15122, -0.06144, * -14.68939,-14.82251,-13.65846, -0.11173, -0.14410, -0.07133, * -18.38628,-18.94631,-19.00893, -0.08062, -0.14481, -0.12949, * -16.15328,-17.40999,-14.08705, -0.08485, -0.06896, -0.11583, * -14.50295,-16.91671,-25.25793, -0.06814, -0.13727, -0.12213, * -10.92188,-14.10852,-24.43877, -0.09375, -0.11638, -0.09053, * -11.64716,-14.92020,-19.99063, -0.14792, -0.08681, -0.12085, * -24.09766,-16.14519, -8.05683, -0.24065, -0.05877, -0.23726, * -25.18396,-15.02034,-15.50531, -0.12236, -0.09610, -0.00529, * -15.27905,-19.36708,-12.94046, -0.08571, -0.09560, -0.03544, * -7.48927,-16.00753,-13.02842, -0.07862, -0.10110, -0.05807, * -13.06383,-27.98698,-18.80004, -0.05875, -0.03737, -0.11214, * -13.67370,-16.44925,-16.12632, -0.07228, -0.09322, -0.05652, * -22.61245,-21.24717,-18.09933, -0.05197, -0.07477, -0.05235, * -27.09189,-21.85181,-20.34676, -0.05123, -0.05683, -0.07214, * -27.09561,-22.76383,-25.41151, -0.10272, -0.02058, -0.16720/ call ggg(param,funct,xl) C ********************************** y=0. C ********************************** do i=1,index_t do il=1,index_l kk=index_l*(i-1)+il do j=1,nfunc ind=nfunc*(kk-1)+j bspl4=bspl4_time(i,xt)*bspl4_long(il,xl) y=y+bspl4*funct(j)*coeff(ind) end do end do end do return end C------------------------------------------------------------------ c ************************************************* c ************************************************* real*8 function bspl4_time(i,x1) c ************************************************* implicit REAL*8 (A-H,O-Z) integer i,order/4/,j,k real*8 t_t(0:39) real*8 x,b(20,20),x1 data t_t/ * 0.00,2.75,4.75,5.50,6.25, * 7.25,10.00,14.00,17.25,18.00, * 18.75,19.75,21.00,24.00,26.75, * 28.75,29.50,30.25,31.25,34.00, * 38.00,41.25,42.00,42.75,43.75, * 45.00,48.00,50.75,52.75,53.50, * 54.25,55.25,58.00,62.00,65.25, * 66.00,66.75,67.75,69.00,72.00/ x=x1 if(i.ge.0) then if (x.lt.t_t(i-0)) then x=x+24 end if end if do j=i,i+order-1 if(x.ge.t_t(j).and.x.lt.t_t(j+1)) then b(j,1)=1 else b(j,1)=0 end if end do do j=2,order do k=i,i+order-j b(k,j)=(x-t_t(k))/(t_t(k+j-1)-t_t(k))*b(k,j-1) b(k,j)=b(k,j)+(t_t(k+j)-x)/(t_t(k+j)-t_t(k+1))*b(k+1,j-1) end do end do bspl4_time=b(i,order) return end C------------------------------------------------------------------ c ************************************************* c ************************************************* real*8 function bspl4_long(i,x1) c ************************************************* implicit real*8 (A-H,O-Z) integer i,order/4/,j,k real*8 t_l(0:24) real*8 x,b(20,20),x1 data t_l/ * 0,10,100,190,200,250,280,310, * 360,370,460,550,560,610,640,670, * 720,730,820,910,920,970,1000,1030,1080/ x=x1 if(i.ge.0) then if (x.lt.t_l(i-0)) then x=x+360 end if end if do j=i,i+order-1 if(x.ge.t_l(j).and.x.lt.t_l(j+1)) then b(j,1)=1 else b(j,1)=0 end if end do do j=2,order do k=i,i+order-j b(k,j)=(x-t_l(k))/(t_l(k+j-1)-t_l(k))*b(k,j-1) b(k,j)=b(k,j)+(t_l(k+j)-x)/(t_l(k+j)-t_l(k+1))*b(k+1,j-1) end do end do bspl4_long=b(i,order) return end C------------------------------------------------------------------ c ************************************************* c ************************************************* subroutine ggg(param,funct,x) c ************************************************* implicit real*8 (A-H,O-Z) integer i real*8 param(2),funct(6) real*8 x,a,sigma,gauss,flux,cflux c ************************************************* flux=param(2) if(param(2).le.75) flux=75. if(param(2).ge.230) flux=230. cflux=flux a=0. if((param(1).ge.120).and.(param(1).le.240)) a=170. if((param(1).ge.120).and.(param(1).le.240)) sigma=60 if((param(1).le.60).or.(param(1).ge.300)) a=170. if((param(1).le.60).or.(param(1).ge.300)) sigma=40 if((flux.le.95).and.(a.ne.0)) then gauss=exp(-0.5*((x-a)**2)/sigma**2) cflux=gauss*95.+(1-gauss)*flux end if c ************************************************* c ************************************************* do i=1,6 funct(i)=0. end do c ************************************************* c ************************************************* if((param(1).ge.135).and.(param(1).le.230)) funct(1)=1 if((param(1).le.45).or.(param(1).ge.320)) funct(2)=1 if((param(1).gt.75).and.(param(1).lt.105)) funct(3)=1 if((param(1).gt.260).and.(param(1).lt.290)) funct(3)=1 c ************************************************* if((param(1).ge.45).and.(param(1).le.75)) then ! W-E funct(2)=1.-(param(1)-45.)/30. funct(3)=1-funct(2) end if if((param(1).ge.105).and.(param(1).le.135)) then ! E-S funct(3)=1.-(param(1)-105.)/30. funct(1)=1-funct(3) end if if((param(1).ge.230).and.(param(1).le.260)) then ! S-E funct(1)=1.-(param(1)-230.)/30. funct(3)=1-funct(1) end if if((param(1).ge.290).and.(param(1).le.320)) then ! E-W funct(3)=1.-(param(1)-290.)/30. funct(2)=1-funct(3) end if c ************************************************* funct(4)=(cflux-140)*funct(1) funct(5)=(cflux-140)*funct(2) funct(6)=(flux-140)*funct(3) c ************************************************* return end C ********The End of the Code ************************