PROGRAM main use geogrid_t3d ,only: nlev,nlat,nlon use maggrid ,only: nmlat,nmlon,nmlonp1 use edynamo_driver,only: edyn_init,edyn_driver INCLUDE 'EIG.h' COMMON/input/ F107,F107A,AP,Uth,TimeStep,Hindex(6),iyd,IYR,IMO,IDA COMMON/GRIDI/ alti(NLi),Cli(NLi),SIi(NLi),MLATi(lmi+1,NLi),Mloni( & NFi),altH(NHN),ClH(NHN),SIH(NHN),MlatH(NHN,2,NLi), & MlonH(NHN) real*8,allocatable,dimension(:,:,:) :: Z,Tn,cO,cO2,cN2,cH,cHe,cN, & Ws,We,Wu,Qop,QHep,OpL,Qep,QHp,HepL,HpL, & cO2p,cNOp,cN2p,Op,Te,Ti,Vi,Ped,Hall REAL*8 ed1(nmlat,nmlonp1),ed2(nmlat,nmlonp1),zigm11(nmlat,nmlonp1 & ),zigm22(nmlat,nmlonp1),zigm2(nmlat,nmlonp1),zigmc(nmlat, & nmlonp1),rim1(nmlat,nmlonp1),rim2(nmlat,nmlonp1) ! ! Set use_conduct = 'edynamo' to have edynamo do fieldline integrations, ! or set use_conduct = 'time3d' to have TIME3D do fieldline integrations. ! In the latter case, time3d will pass 2d conductances to edynamo. ! ! character(len=16) :: use_conduct='edynamo ' character(len=16) :: use_conduct='time3d ' write(6,"('time3d: use_conduct=',a)") use_conduct CALL PreGrid() !Set parameters for geogrid, remove when coupling with Dynamo allocate(Z(nlat,nlon,Nlev)) allocate(Tn(nlat,nlon,Nlev)) allocate(cO(nlat,nlon,Nlev)) allocate(cO2(nlat,nlon,Nlev)) allocate(cN2(nlat,nlon,Nlev)) allocate(cH(nlat,nlon,Nlev)) allocate(cHe(nlat,nlon,Nlev)) allocate(cN(nlat,nlon,Nlev)) allocate(Ws(nlat,nlon,Nlev)) allocate(We(nlat,nlon,Nlev)) allocate(Wu(nlat,nlon,Nlev)) allocate(Qop(nlat,nlon,Nlev)) allocate(QHep(nlat,nlon,Nlev)) allocate(OpL(nlat,nlon,Nlev)) allocate(Qep(nlat,nlon,Nlev)) allocate(QHp(nlat,nlon,Nlev)) allocate(HepL(nlat,nlon,Nlev)) allocate(HpL(nlat,nlon,Nlev)) allocate(cO2p(nlat,nlon,Nlev)) allocate(cNOp(nlat,nlon,Nlev)) allocate(cN2p(nlat,nlon,Nlev)) allocate(Op(nlat,nlon,Nlev)) allocate(Te(nlat,nlon,Nlev)) allocate(Ti(nlat,nlon,Nlev)) allocate(Vi(nlat,nlon,Nlev)) allocate(Ped(nlat,nlon,Nlev)) allocate(Hall(nlat,nlon,Nlev)) IYR=2001 !year iyd=1 !doy of year ! ! Timestep in seconds (default is 300, or 5 minutes): ! TimeStep=3.0D2 ! time step TimeStep=6.0D2 CALL GCMpara() !Set parameters for GCM, remove when coupling with GCM model call edyn_init CALL TIME3DParaGrid(IYR,iyd,HB) ! Set TIME3D Parameters and Grid System Do ifsn=3,3 C--------------------Parameters for GCM and dynamo, remove when coupling with GCM and dynamo model-------------------- F107=65.0D0 F107A=F107 AP=4.0D0 iyd=int(30.5D0*ifsn)-9 !87 IYR=2001 UTh=0.0D0 Hindex(1)=0.0D0 Hindex(2)=-0.5D0 Hindex(3)=3.0D2 Hindex(4:5)=1.0D1 Hindex(6)=2.0D1 ih=3600.0D0/Timestep Nnum=ih*24 write(*,*)ifsn,F107,iyd CALL DoytoMonDay(IYR,iyd,IMO,IDA) CALL EUVFluxCrossSection(F107,F107A,iyd) CALL Prealt(Z) C----------------------------------------------------------------------------- CALL TIME3DInitialCondition(iyd,F107,F107A,AP,Uth) !TIME3D Initilze CALL TIME3D2GCM(Z,Op,Te,Ti,Vi,Ped,Hall) ! output Op,Te,Ti,Vi for GCM code ed1=0.0D0 ed2=0.0D0 ! DO iday=1,3 DO iday=1,1 DO i=1,Nnum UTh=DBLE(i)*Timestep/3.6D3 IF(UTh.GT.2.4d1)UTh=UTh-2.4d1 ! Instead of GCM code, remove when coupling with GCM ! input: Op: O+ density in m-3; Te,Ti: electron and ion temperature in K ! output: Z: altitude in meter; Tn: neutral temperature in K; ! cO,cO2,cN2,cH,cHe,cN: O, O2, N2, H, He, and N densities in m-3 (in now cersion, [He] is set to be 0 or 1.0D-25) ! Ws,We,Wu: southward, eastward, and upward wind velocities in m/s; ! Qop,QHep,QHp,Qep: production rates of O+,He+,H+, and photoelectron in M-3*s-1 (in now cersion, QHep,QHp are set to be 0 or 1.0D-25) ! OpL,HepL,HpL: loss rates of O+,He+,and H+ in s-1 (in now cersion, HepL,HpL are set to be 0 or 1.0D-25) ! cO2p,cNOp,cN2p: O2+, NO+, and N2+ density in m-3 CALL GCMsim(Z,Tn,cO,cO2,cN2,cH,cHe,cN,Ws,We,Wu,Qop,QHep,OpL, & Qep,QHp,HepL,HpL,cO2p,cNOp,cN2p,Op,Te,Ti) ! ! TIME3D module ! Input: TimeStep: time step; Z: altitude in meter; Tn: neutral temperature in K; ! cO,cO2,cN2,cH,cHe,cN: O, O2, N2, H, He, and N densities in m-3 (in now cersion, [He] is set to be 0 or 1.0D-25) ! Ws,We,Wu: southward, eastward, and upward wind velocities in m/s; ! Qop,QHep,QHp,Qep: production rates of O+,He+,H+, and photoelectron in M-3*s-1 (in now cersion, QHep,QHp are set to be 0 or 1.0D-25) ! OpL,HepL,HpL: loss rates of O+,He+,and H+ in s-1 (in now cersion, HepL,HpL are set to be 0 or 1.0D-25) ! cO2p,cNOp,cN2p: O2+, NO+, and N2+ density in m-3; (above are output from GCM) ! ed1,ed2: electric field index; (output from dynamo) ! output: Op: O+ density in m-3; Te,Ti: electron and ion temperature in K; Vi: anti-field velocity of O+; (input into GCM) ! Ped,Hall: Pedersen and Hall conductivities in s/m; (input into dynamo and GCM) ! zigm11,zigm22,zigm2,zigmc,rim1,rim2; dynamo index in S or S*m/s*Be3 (input into dynamo); the output of [subroutine fieldline_integrals] in dynamo CALL TIME3Dmodule(TimeStep,Z,Tn,cO,cO2,cN2,cH,cHe,cN,Ws,We, & Wu,Qop,QHep,OpL,Qep,QHp,HepL,HpL,cO2p,cNOp & ,cN2p,ed1,ed2,Op,Te,Ti,Vi,Ped,Hall,zigm11, & zigm22,zigm2,zigmc,rim1,rim2,t) ! ! Instead of dynamo code, remove when coupling with dynamo ! output: ed1,ed2;(input into TIME3D) input: HB (need not in dymamo code) ! write(6,"(i12,i12,' 2')")iday,i ! CALL DynamoSim(ed1,ed2,HB) ! c IF((mod(i,ih).eq.0).and.(iday.eq.3))call output(ifsn,i/ih) c IF((mod(i,ih).eq.0).and.(iday.eq.3))call OUTG(ifsn,i/ih,Op) ! WRITE(*,'(3I5,F8.3)')ifsn,iday,i,UTH ! ! If use_conduct = 'edynamo', then edynamo will do the fieldline integrations, ! If use_conduct = 'time3d', then pass the 2d conductances from time3d to edynamo ! (edynamo will not do the fieldline integrations in this case). ! (The last 6 arguments to edyn_driver are optional). ! if (trim(use_conduct)=='edynamo') then call edyn_driver(ed1,ed2) elseif (trim(use_conduct)=='time3d') then call edyn_driver(ed1,ed2,zigm11,zigm22,zigm2,zigmc, | rim1,rim2) else write(6,"('>>> main TIME3D: unknown use_conduct=',a, | ' (must be either ''edynamo'' or ''time3d'')')") | trim(use_conduct) endif write(6,"('Main after edyn_driver: i=',i4,' Nnum=',i4, | ' iday=',i3,' uth=',f9.3)") i,Nnum,iday,uth ! write(6,"(i12,i12,f12.3' 3')")iday,i,Uth ENDDO ! i=1,Nnum ENDDO ! iday=1,3 ENDDO ! ifsn=3,3 WRITE(*,*) 'IGG3D model normally terminated!' END PROGRAM main C----------------------------------- Ootput DATA ------------------------------- SUBROUTINE OUTput(ifn,Iuth) INCLUDE 'EIG.h' character*20 dirn character id1,id2,id3,id4 integer ifn,Iuth real*8 xx(12) COMMON/REIM3/ REIM(lmi,NFi,NLi,12),REIMold(lmi,NFi,NLi,8),REIH(NHN & ,2,NLi,12),REIHold(NHN,2,NLi,8) i1=mod(ifn,10) i2=mod((ifn-i1)/10,10) id1=char(i1+48) id2=char(i2+48) i1=mod(Iuth,10) i2=mod((Iuth-i1)/10,10) id3=char(i1+48) id4=char(i2+48) kx=23 open(100,file='ReI'//id2//id1//'_'//id4//id3//'.dat') DO j=1,12 DO k=1,NLi DO i=k,lmi-k+1 write(100,101)REIM(i,1:NFi,k,j) ENDDO ENDDO ENDDO close(100) open(100,file='ReH'//id2//id1//'_'//id4//id3//'.dat') DO i=1,12 DO j=1,2 DO k=1,NLi write(100,102)REIH(1:NHN,j,k,i) ENDDO ENDDO ENDDO close(100) 101 format(72d20.10) 102 format(150d20.10) END SUBROUTINE Output C----------------------------------- Ootput DATA ------------------------------- SUBROUTINE OUTG(ifn,Iuth,op) use geogrid_t3d,only: nlev,nlat,nlon INCLUDE 'EIG.h' character*20 dirn character id1,id2,id3,id4 integer ifn,Iuth real*8 xy(nlon),Op(nlat,nlon,Nlev) i1=mod(ifn,10) i2=mod((ifn-i1)/10,10) id1=char(i1+48) id2=char(i2+48) i1=mod(Iuth,10) i2=mod((Iuth-i1)/10,10) id3=char(i1+48) id4=char(i2+48) kx=23 open(100,file='outputg\ReG'//id2//id1//'_'//id4//id3//'.dat') DO i=1,nlat DO j=1,nlon xmax=Op(i,j,1) DO k=2,nlev if(xmax.lt.Op(i,j,k))xmax=Op(i,j,k) ENDDO xy(j)=xmax ENDDO write(100,102)xy(1:nlon) c write(100,102)op(i,1:nlon,22) ENDDO 102 format(150d20.10) close(100) END SUBROUTINE OUTG