! module advance_module implicit none ! ! Advance the model in time. ! contains !----------------------------------------------------------------------- subroutine advance ! ! Advance the model nstep time steps. ! use fields_module use params_module,only: nlon,nlat,nlev,nlevp1 use hist_module,only: nstep,modeltime,nsource use init_module,only: istep,uthr,iter,secs,iday,iyear use input_module,only: start,step,iaurora=>aurora,idynamo=>dynamo, | f107,f107a,ctpoten,power,calendar_advance,start_year,start_day, | potential_model,ntimes_ctpoten, ctpoten_time, ntimes_power, | power_time,ntimes_f107, f107_time, ntimes_f107a, f107a_time, | meped_file,gpi_ncfile,step_compqrj use cons_module,only: dt use magfield_module,only: sunloc_apex use timing_module,only: timer,timing use gpi_module,only: rdgpi,getgpi use solgar_module,only: solgar_bndry,solgar_import use aurora_module,only: aurora_cons,echar_meped, epower_meped use qrj_module,only: init_sflux use efield_module,only: efield use hdif_module,only: hdif1,hdif2 use mgw_module,only: backgrnd use dynamo_module,only: prep_dynamo,prep_dynamo_dyn0,dynamo, | nodynamo,phihm,nmlat0,zpotenm3d use heelis_module,only: heelis,colath use weimer_module,only: weimer01 use output,only: outhist use meped,only: rdmeped use addfld_module,only: addfld use lbc,only: tuvz_lbc ! ! current due to plasma pressure & gravity -> add to RHS of dynamo use magpres_g_module, only: magpres_grav,j_pg ! #ifdef MPI use mpi_module,only: mytid,mp_gather2root,ntask,lat0,lat1,lon0, | lon1,mp_polelats,mp_bndlats,mp_bndlons,mp_periodic_f4d, | mp_updatephi,sum_mpi_time implicit none #include #else use mpi_module,only: lat0,lat1,lon0,lon1,mytid,ntask implicit none #endif ! ! Local: integer(kind=8) :: nsecs ! total model time in seconds (includes day) integer :: itmp,ier,i,j,iprint real :: fmin,fmax,dday logical :: time2write,wrprim,wrsech, | newseries_sech,newseries_prim,iseries_sech,iseries_prim ! ! External: logical,external :: time2print ! real :: secs_per_step, | time0, time1, | time0_step, time1_step, | time0_dynamics, time1_dynamics, | time0_dynamo, time1_dynamo, | time0_phist, time1_phist, | time0_shist, time1_shist, | time0_prep, time1_prep, | time0_init, time1_init character(len=3) :: timing_type ! ! logical :: debug=.true. logical :: debug=.false. ! ! External: logical,external :: wrhist ! write(6,"(/,'Enter advance: iter=',i6,' nstep=',i6)") iter,nstep timing_type = 'sys' if (timing%rtc) timing_type = 'rtc' call timer(time0_init,time1_init,'INIT',0,0) ! ! Init timestep index, modeltime, and nsecs: istep = 0 modeltime(1:3) = start(:,1) modeltime(4) = 0 call modeltime_to_nsecs(modeltime,nsecs) ! sets nsecs, includes day ! ! Source was read into data(:,:,:,itp), and output always writes from ! data(:,:,:,itc), so swap indices to echo itp source data: itmp = itp itp = itc itc = itmp ! ! Start phist timing (source hist may or may not be written): call timer(time0_phist,time1_phist,'PHIST',0,0) ! start phist ! ! Echo source history to primary history output if an initial run, time2write = wrhist(istep,modeltime, | wrprim, newseries_prim, iseries_prim, | wrsech, newseries_sech, iseries_sech) if (time2write) then #ifdef MPI ! ! If an initial run, echo source history to output before first time step. ! Because only subdomains are read by rdhist, and only the root task ! outputs histories, the subdomain data must be gathered to the root. ! call mp_gather2root(itc,'prim') if (mytid==0) call outhist(istep,modeltime) #else do i=1,nf4d foutput(:,lon0:lon1,lat0:lat1,i) = ! foutput is in fields.F | f4d(i)%data(:,lon0:lon1,lat0:lat1,itc) enddo tlbc_glb(lon0:lon1,lat0:lat1) = tlbc(lon0:lon1,lat0:lat1) ulbc_glb(lon0:lon1,lat0:lat1) = ulbc(lon0:lon1,lat0:lat1) vlbc_glb(lon0:lon1,lat0:lat1) = vlbc(lon0:lon1,lat0:lat1) call outhist(istep,modeltime) #endif ! ! nstep may be zero if user only wants to copy source history: if (nstep==0) then write(6,"(/,'Advance: model is not advanced in time because ', | ' start==stop.')") return endif endif ! time2write call timer(time0_phist,time1,'PHIST',1,0) ! suspend phist timing time1_phist = time1 ! ! Reswap indices, so model reads itp data and updates itc data: itmp = itp itp = itc itc = itmp ! call timer(time0_init,time1_init,'INIT',1,0) ! end init timing if (timing%level >= 2) | write(6,"('Time in INIT = ', | f12.2,' Advance: step ',i5)") time1_init,istep ! ! Main time loop: ! 100 continue ! write(6,"('advance call getoverstepped: begin istep=',i5)") istep ! call getoverstepped() call timer(time0_step,time1_step,'STEP',0,0) ! start step timing call timer(time0_prep,time1_prep,'PREP',0,0) ! start prep timing iter = iter+1 istep = istep+1 iprint = 0 if (istep==1) iprint = 1 nsecs = nsecs+step ! ! Note model time is incremented even when calendar_advance==0 call nsecs_to_modeltime(nsecs,modeltime) ! increments modeltime(4) secs = amod(float(iter)*dt,86400.) uthr = secs/3600. ! ! Report to stdout (data is actually from previous step): if (time2print(nstep,istep)) then if (istep > 1) then write(6,"('Step ',i6,' of ',i6,' mtime=',3i3, | ' secs/step (',a,') =',f6.2)") istep,nstep,modeltime(1:3), | timing_type,secs_per_step else write(6,"('Step ',i6,' of ',i6,' mtime=',3i3)") | istep,nstep,modeltime(1:3) endif endif ! time2print ! ! Update calendar day and sun's location and declination: ! (sub advance_day also updates sfeps) ! (note iday is incremented only if calendar_advance > 0) ! if (calendar_advance > 0) call advance_day call sunloc_apex(iyear,iday,secs)! locate sun's longitude ! ! Update lower boundaries of t,u,v,z. ! call tuvz_lbc ! ! Get Solomon-Garcia lower boundaries and imports: ! call solgar_bndry(iyear,iday,int(secs)) call solgar_import(iyear,iday,int(secs)) ! ! Get gpi data if necessary: if (len_trim(gpi_ncfile) > 0) then if (istep==1) call rdgpi ! read gpi data call getgpi(iyear,iday,int(secs),iprint) if (time2print(nstep,istep)) then write(6,"(' GPI: f107=',f7.2,' f107a=',f7.2, | ' ctpoten=',f7.2,' power=',f7.2)") | f107,f107a,ctpoten,power endif endif ! ! Read MEPED data if this is first timestep (see meped.F): if (istep==1.and.len_trim(meped_file) > 0) | call rdmeped(meped_file) ! ! Interpolate power and/or ctpoten to current model time, if time-dependent ! values were read from namelist input: ! if (ntimes_f107 > 0) then call set_index(f107_time,ntimes_f107,nsecs,f107, | 'f107') if (time2print(nstep,istep)) | write(6,"(' User provided f107 =',f8.2)") f107 endif if (ntimes_f107a > 0) then call set_index(f107a_time,ntimes_f107a,nsecs,f107a, | 'f107a') if (time2print(nstep,istep)) | write(6,"(' User provided f107a =',f8.2)") f107a endif if (ntimes_ctpoten > 0) then call set_index(ctpoten_time,ntimes_ctpoten,nsecs,ctpoten, | 'ctpoten') if (time2print(nstep,istep)) | write(6,"(' User provided ctpoten=',f8.2)") ctpoten endif if (ntimes_power > 0) then call set_index(power_time,ntimes_power,nsecs,power,'power') if (time2print(nstep,istep)) | write(6,"(' User provided power =',f8.2)") power endif ! ! Update auroral parameters: ! call aurora_cons(iprint,iyear,iday,secs) if (len_trim(meped_file) > 0.and.time2print(nstep,istep)) | write(6,"(' MEPED epower=',f10.2,' echar=',f10.2)") | epower_meped,echar_meped ! ! Init gravity wave parameterization (mgw.F, see also mgwinit): ! call backgrnd ! ! init_sflux calls ssflux and euvac (in qrj_module, qrj.F): ! if (istep==1.or.mod(istep,step_compqrj/step)==0) | call init_sflux(iprint) ! ! Calculate electric field (ex,ey,ez in efield module, efield.F): ! call efield(1,nlevp1,lon0,lon1,lat0,lat1,iprint) ! ! Calculate vc, barm, xnmbar, and z (tn,o2,ox,vn are input). ! call addiag( | tn (levd0,lond0,latd0,itp), ! in | o2 (levd0,lond0,latd0,itp), ! in | ox (levd0,lond0,latd0,itp), ! in | vn (levd0,lond0,latd0,itp), ! in | vc (levd0,lond0,latd0,itp), ! out | barm (levd0,lond0,latd0,itp), ! out | xnmbar (levd0,lond0,latd0), ! out | xnmbari(levd0,lond0,latd0), ! out | xnmbarm(levd0,lond0,latd0), ! out | z (levd0,lond0,latd0,itp), ! out (note z output is at itp) | zg (levd0,lond0,latd0), ! out (z with varying gravity) | lon0,lon1,1,nlevp1,lat0,lat1) ! write(6,"('advance call getoverstepped: after addiag')") ! call getoverstepped() ! ! do j=lat0,lat1 ! call addfld('un',' ',' ',un(:,lon0:lon1,j,itp), ! | 'lev',1,nlevp1,'lon',lon0,lon1,j) ! call addfld('vc',' ',' ',vc(1:nlev,lon0:lon1,j,itp), ! | 'lev',1,nlev,'lon',lon0,lon1,j) ! call addfld('advbarm',' ',' ',barm(:,lon0:lon1,j,itp), ! | 'lev',1,nlevp1,'lon',lon0,lon1,j) ! enddo ! #ifdef MPI ! ! Define boundary latitudes -1,0 across the south pole, and ! nlat+1,nlat+2 across the north pole: ! call mp_polelats(itp) ! ! Update subdomain boundary latitudes, longitudes: f4d(:)%mpi = .true. call mp_bndlats(f4d,nf4d,itp) call mp_bndlons(f4d,nf4d,itp) ! write(6,"('advance call getoverstepped: after mp_bndlons')") ! call getoverstepped() #else call mk_polelat( 0,1,itp) call mk_polelat(-1,2,itp) call mk_polelat(nlat+1,nlat ,itp) call mk_polelat(nlat+2,nlat-1,itp) #endif ! ! For Z, itc==itp (itp was set by addiag): z(:,:,:,itc) = z(:,:,:,itp) ! ! Horizontal diffusion, first step: ! hdif1 saves nrh and kmh at lats -2->nlat in nrh and kmh using ! fields un_nm and vn_nm at j+1 and j+2, and tn_nm and barm at j+1. ! hdif1 calls mp_bnd[lat,lons]_kmh ! call hdif1( | tn_nm(levd0,lond0,latd0,itp), | un_nm(levd0,lond0,latd0,itp), | vn_nm(levd0,lond0,latd0,itp), | barm (levd0,lond0,latd0,itp), | 1,nlevp1,lon0,lon1,lat0,lat1) ! write(6,"('advance call getoverstepped: after hdif1')") ! call getoverstepped() ! ! Horizontal diffusion, second step: ! hdif2 saves 3d diagnostics f3d (fkldt,u,v,o2,o) at j+1 (0->37) ! hdif2: WRITE (fkldt,u,v,o2,o1) <- READ (tn_nm,un_nm,vn_nm,o2_nm,ox_nm) ! The kld terms will be used later in hdif3. ! call hdif2( | tn_nm (levd0,lond0,latd0,itp), ! 4d input | un_nm (levd0,lond0,latd0,itp), ! 4d input | vn_nm (levd0,lond0,latd0,itp), ! 4d input | o2_nm (levd0,lond0,latd0,itp), ! 4d input | ox_nm (levd0,lond0,latd0,itp), ! 4d input | n4s_nm(levd0,lond0,latd0,itp), ! 4d input | noz_nm(levd0,lond0,latd0,itp), ! 4d input | co_nm (levd0,lond0,latd0,itp), ! 4d input | co2_nm(levd0,lond0,latd0,itp), ! 4d input | h2o_nm(levd0,lond0,latd0,itp), ! 4d input | h2_nm (levd0,lond0,latd0,itp), ! 4d input | hox_nm(levd0,lond0,latd0,itp), ! 4d input | ch4_nm(levd0,lond0,latd0,itp), ! 4d input | ar_nm (levd0,lond0,latd0,itp), ! 4d input | he_nm (levd0,lond0,latd0,itp), ! 4d input | nat_nm(levd0,lond0,latd0,itp), ! 4d input | kldt , ! 3d output | kldu , ! 3d output | kldv , ! 3d output | kldo2, ! 3d output | kldox, ! 3d output | kldn4s, ! 3d output | kldnoz, ! 3d output | kldco, ! 3d output | kldco2, ! 3d output | kldh2o, ! 3d output | kldh2, ! 3d output | kldhox, ! 3d output | kldch4, ! 3d output | kldar, ! 3d output | kldhe, ! 3d output | kldnat, ! 3d output | 1,nlevp1,lon0,lon1,lat0,lat1) ! write(6,"('advance call getoverstepped: after hdif2')") ! call getoverstepped() ! ! Dynamo calls Heelis (heelis.F), Weimer (wei01gcm.F), or neither ! for high latitude electric potential, according to user-provided ! namelist parameter "model_potential". ! Get high latitude (Heelis or other) colatitudes, NH pfrac, and poten phihm. ! If Weimer is used, then theta0,phid etc is changed before use in aurora ! in dynamics. ! if (potential_model == 'WEIMER') then if (debug) write(6,"('advance calling weimer: istep=',i5)") | istep call weimer01 elseif (potential_model == 'HEELIS') then if (debug) write(6,"('advance calling heelis: istep=',i5)") | istep call heelis else ! potential_model='NONE' do j=1,nmlat0 do i=1,nmlonp1 phihm(i,j) = 0. enddo ! i=1,nmlonp1 enddo ! j=1,nmlat0 call colath endif call timer(time0_prep,time1,'PREP',1,0) ! suspend step-prep timing time1_prep = time1 ! ! Model dynamics. Sub dynamics calls physics, chemistry, and dynamics ! routines for current time step (dynamics.F). ! if (debug) write(6,"('advance calling dynamics: istep=',i5)") | istep call timer(time0_dynamics,time1_dynamics,'DYNAMICS',0,0) ! write(6,"('advance call getoverstepped: before dynamics')") ! call getoverstepped() call dynamics(nstep,istep) ! write(6,"('advance call getoverstepped: after dynamics')") ! call getoverstepped() call timer(time0_dynamics,time1_dynamics,'DYNAMICS',1,0) if (debug) write(6,"('advance after dynamics: istep=',i5)") | istep if (timing%level >= 2.and.time2print(nstep,istep)) | write(6,"('Time in DYNAMICS = ',f12.2, | ' Advance: step ',i5)") time1_dynamics,istep ! ! Call dynamo if namelist input flag is set: call timer(time0_dynamo,time1_dynamo,'DYNAMO',0,0) ! ! zpotenm3d is 3-d geopotential height in mag coords. This field is forced ! on the secondary histories. Here init to zero: ! zpotenm3d = 0. ! whole array init (this field forced onto secondary hist) select case (idynamo) ! ! idynamo==0 -> Do not call dynamo. However, if potential_model is ! requested, set electric potential accordingly. case (0) if (trim(potential_model) /= 'NONE') then if (debug) write(6,"('advance calling prep_dynamo_dyn0', | ' istep=',i5)") istep call prep_dynamo_dyn0(z(levd0,lond0,latd0,itp), | 1,nlevp1,lon0,lon1,lat0,lat1) if (debug) write(6,"('advance after prep_dynamo_dyn0')") if (mytid==0) then if (debug) write(6,"('advance calling nodynamo: istep=', | i5)") istep call nodynamo if (debug) write(6,"('advance after nodynamo: istep=', | i5)") istep endif endif ! potential_model /= NONE ! ! idynamo==1 -> All tasks call prep_dynamo, and master task calls dynamo: ! case (1) ! ! calculate addition to electrodynamo RHS (current due to plasma ! pressure and gravity) ! if(j_pg) then call magpres_grav ( | z (levd0,lond0,latd0,itp), ! geopotential input | te (levd0,lond0,latd0,itp), ! electron temperature K | ti (levd0,lond0,latd0,itp), ! ion temperature K | ne (levd0,lond0,latd0,itp), ! electron density 1/cm^3 | op (levd0,lond0,latd0,itp), ! O+ 1/cm^3 | nplus (levd0,lond0,latd0), ! N+ 1/cm^3 | n2p (levd0,lond0,latd0), ! N2+ 1/cm^3 | nop (levd0,lond0,latd0), ! NO+ 1/cm^3 | o2p (levd0,lond0,latd0,itp), ! O2+ 1/cm^3 | 1,nlevp1,lon0,lon1,lat0,lat1) if (debug) write(6,"('advance after magpres_grav: istep=', | i4)") istep endif ! write(6,"('advance call getoverstepped: after magpres_grav')") ! call getoverstepped() call prep_dynamo( | tn (levd0,lond0,latd0,itp), | un (levd0,lond0,latd0,itp), | vn (levd0,lond0,latd0,itp), | w (levd0,lond0,latd0,itp), | z (levd0,lond0,latd0,itp), | barm (levd0,lond0,latd0,itp), | ped (levd0,lond0,latd0), ! from lamdas.F | hall (levd0,lond0,latd0), ! from lamdas.F | 1,nlevp1,lon0,lon1,lat0,lat1) ! write(6,"('advance call getoverstepped: after prep_dynamo')") ! call getoverstepped() if (debug) write(6,"('advance after prep_dynamo: istep=',i5)") | istep if (mytid==0) call dynamo ! write(6,"('advance call getoverstepped: after dynamo')") ! call getoverstepped() if (debug) write(6,"('advance after dynamo: istep=',i5)") | istep ! ! Unknown idynamo flag: ! case default write(6,"('>>> advance: unknown idynamo = ',i3)") idynamo end select ! ! End timing for dynamo and step (before i/o): call timer(time0_dynamo,time1_dynamo,'DYNAMO',1,0) if (timing%level >= 2.and.time2print(nstep,istep)) | write(6,"('Time in DYNAMO = ', | f12.2,' Advance: step ',i5)") time1_dynamo,istep ! ! Update slave tasks with new phim3d from the dynamo. Make this call ! after timing, so that timing will show the difference between master ! and slave tasks as a consequence of the serial dynamo. (slaves will ! be blocked here until the master is done with the dynamo) ! #ifdef MPI if (idynamo > 0 .or. | (idynamo <= 0 .and. trim(potential_model) /= 'NONE')) | call mp_updatephi ! write(6,"('advance call getoverstepped: after mp_updatephi')") ! call getoverstepped() #endif #ifdef VT ! code = 127 ; state = 'timestep' ; activity='ModelCode' call vtend(127,ier) #endif ! ! Write output history if necessary (master task writes data%(:,:,:,itc)). ! time2write = wrhist(istep,modeltime, | wrprim, newseries_prim, iseries_prim, | wrsech, newseries_sech, iseries_sech) ! ! History i/o timing is not included in STEP segment: call timer(time0_step,time1,'STEP',1,0) ! suspend step timing for i/o time1_step = time1 ! ! If its time to write a history, root task must gather subdomain ! data from slave tasks before writing the history. This is done ! by sub mp_gather2root. Only root task writes to the history. ! if (time2write) then if (wrprim) call timer(time0_phist,time1_phist,'PHIST',0,0) ! resume phist if (wrsech) call timer(time0_shist,time1_shist,'SHIST',0,0) #ifdef MPI if (ntask > 1) then call mp_gather2root(itc,'prim') if (wrsech) call mp_gather2root(itc,'sech') elseif (ntask==1) then do i=1,nf4d foutput(:,lon0:lon1,lat0:lat1,i) = | f4d(i)%data(:,lon0:lon1,lat0:lat1,itc) enddo tlbc_glb(lon0:lon1,lat0:lat1) = tlbc(lon0:lon1,lat0:lat1) ulbc_glb(lon0:lon1,lat0:lat1) = ulbc(lon0:lon1,lat0:lat1) vlbc_glb(lon0:lon1,lat0:lat1) = vlbc(lon0:lon1,lat0:lat1) tlbc_nm_glb(lon0:lon1,lat0:lat1)=tlbc_nm(lon0:lon1,lat0:lat1) ulbc_nm_glb(lon0:lon1,lat0:lat1)=ulbc_nm(lon0:lon1,lat0:lat1) vlbc_nm_glb(lon0:lon1,lat0:lat1)=vlbc_nm(lon0:lon1,lat0:lat1) endif #else ! ! Update foutput if serial non-MPI run: do i=1,nf4d foutput(:,lon0:lon1,lat0:lat1,i) = | f4d(i)%data(:,lon0:lon1,lat0:lat1,itc) enddo tlbc_glb(lon0:lon1,lat0:lat1) = tlbc(lon0:lon1,lat0:lat1) ulbc_glb(lon0:lon1,lat0:lat1) = ulbc(lon0:lon1,lat0:lat1) vlbc_glb(lon0:lon1,lat0:lat1) = vlbc(lon0:lon1,lat0:lat1) tlbc_nm_glb(lon0:lon1,lat0:lat1) = tlbc_nm(lon0:lon1,lat0:lat1) ulbc_nm_glb(lon0:lon1,lat0:lat1) = ulbc_nm(lon0:lon1,lat0:lat1) vlbc_nm_glb(lon0:lon1,lat0:lat1) = vlbc_nm(lon0:lon1,lat0:lat1) #endif ! ! Write to history files: if (mytid==0) call outhist(istep,modeltime) ! ! Update i/o timing: if (wrprim) then call timer(time0_phist,time1,'PHIST',1,0) ! end phist timing time1_phist = time1_phist+time1 endif if (wrsech) call timer(time0_shist,time1_shist,'SHIST',1,0) endif ! time2write call timer(time0_step,time1,'STEP',0,0) ! resume step timing #ifdef MPI ! ! Each mpi task must receive the 2 lats before its first ! updated lat (lat0-1,lat0-2), and the 2 lats after its last ! updated lat (lat1+1,lat2+2). Each task must also send its ! first 2 (lat0,lat0+1) and last 2 (lat1,lat1-1) updated lats ! to the appropriate "adjacent" tasks. ! call timer(time0_prep,time1,'PREP',0,0) ! resume prep timing f4d(:)%mpi = .true. call mp_bndlats(f4d,nf4d,itc) ! ! Same for longitude: call mp_bndlons(f4d,nf4d,itc) ! ! Periodic points for all updated fields: call mp_periodic_f4d(itc) #else ! Non-mpi serial run: call mk_polelat(0 ,1 ,itc) call mk_polelat(-1 ,2 ,itc) call mk_polelat(lat1+1,lat1 ,itc) call mk_polelat(lat1+2,lat1-1,itc) call set_periodic_f4d(itc) #endif ! ! Swap field data time indices, so current updated data becomes previous ! step data for next step: itmp = itp itp = itc itc = itmp ! call timer(time0_prep,time1,'PREP',1,0) ! end prep timing time1_prep = time1_prep+time1 call timer(time0_step,time1,'STEP',1,0) ! end step timing time1_step = time1_step+time1 if (timing%level >= 2.and.time2print(nstep,istep)) then write(6,"('Time in PREP = ', | f12.2,' Advance: step ',i5)") time1_prep,istep write(6,"('Time in STEP = ', | f12.2,' Advance: step ',i5)") time1_step,istep endif secs_per_step = time1_step #ifdef MPI timing%mpi = timing%mpi + sum_mpi_time() ! accumulate mpi timing #endif ! ! Return for next time step: if (istep < nstep) then ! write(6,"('advance call getoverstepped: end istep=',i5)") istep ! call getoverstepped() goto 100 endif end subroutine advance !----------------------------------------------------------------------- subroutine advance_day ! ! Advance calendar day if needed. Also update sfeps. ! use init_module,only: iter,iyear,iday,sfeps,sundec,sin_sundec, | cos_sundec use cons_module,only: dt,pi implicit none ! ! Local: integer :: idayit,idayprev,iyearprev,iyr4,iyr100,lpyr,ienda real :: theta0 ! idayit = iter*int(dt)/86400 if (idayit*86400 == iter*int(dt)) then idayprev = iday iyearprev = iyear iday = iday + 1 ! ! lpyr = 1(0) if is (not) a leap year iyr4 = iyear/4 iyr100 = iyear/100 lpyr = 0 if (iyr4*4 == iyear .and. iyr100*100 /= iyear) lpyr=1 ienda = 365 + lpyr if (iday > ienda) then iyear = iyear + 1 iday = iday - ienda endif ! for past year's end ! ! Recalculate sun's declination ! ! 1/27/05 btf: Changed local delta to sundec (sundec is in init_module, ! where it is initialized at beginning of run). This should not change ! any results, because only sin_sundec and cos_sundec are actually used ! (in chapman), but am updating sundec here so it can be used in future ! code development. ! sundec = atan(tan(23.5*pi/180.)*sin(2.*pi*float(iday-80)/365.)) sin_sundec = sin(sundec) ! C(95) cos_sundec = cos(sundec) ! C(96) ! ! Update sfeps: ! sfeps is 6% variation in solar output over a year ! caused by the orbital eccentricity. ! theta0 = 2.*pi*float(iday)/365. sfeps = 1.000110+0.034221*cos(theta0)+0.001280*sin(theta0) 1 +0.000719*cos(2.*theta0)+0.000077*sin(2.*theta0) ! write(6,"('Advancing day (previous,present)=',4i5,' sfeps=', | e12.4,' sundec=',e12.4)") idayprev,iyearprev,iday,iyear,sfeps, | sundec endif end subroutine advance_day !----------------------------------------------------------------------- subroutine set_index(rindex,ntimes,msecs,outindex,name) ! ! User has provided time-dependent ctpoten (ctpoten_time) and/or ! power (power_time) via namelist input. This routine interpolates ! these inputs to the current model time msecs, returning outindex. ! This routine is called separately (from advance) for ctpoten and power, ! (i.e., rindex will be either ctpoten_time or power_time from input). ! Note times in seconds are 8-byte integers. ! use params_module,only: | mxind_time ! max number of time-dependent solar index points implicit none ! ! Args: real,intent(in) :: | rindex(4,mxind_time) ! user input times and values (day,hour,min,value) integer,intent(in) :: | ntimes ! number of valid time/values in rindex(:,1:ntimes) integer(kind=8),intent(in) :: | msecs ! current model time in seconds real,intent(out) :: outindex ! output interpolated value character(len=*),intent(in) :: name ! ! Local: integer :: i integer(kind=8) :: nsec0,nsec1 ! ! External: integer(kind=8),external :: mtime_to_nsec real,external :: finterp_bigints ! ! If model time is beyond last rindex time, use last rindex time: nsec1 = mtime_to_nsec(int(rindex(1:3,ntimes))) if (msecs > nsec1) then outindex = rindex(4,ntimes) goto 100 endif ! ! Bracket model time: do i=1,ntimes-1 nsec0 = mtime_to_nsec(int(rindex(1:3,i))) nsec1 = mtime_to_nsec(int(rindex(1:3,i+1))) ! ! If model time is at a provided time, interpolation is not necessary: if (nsec0 == msecs) then outindex = rindex(4,i) goto 100 endif if (nsec1 == msecs) then outindex = rindex(4,i+1) goto 100 endif ! ! Interpolate to model time msecs: if (msecs >= nsec0 .and. msecs <= nsec1) then outindex = finterp_bigints(rindex(4,i),rindex(4,i+1),nsec0, | nsec1,msecs) goto 100 endif enddo ! i=1,ntimes-1 ! ! Error if model time could not be bracketed. This should not happen, ! but you never know... write(6,"(/,'>>> set_index: could not bracket model time ', | i10)") msecs write(6,"('ntimes=',i3)") ntimes do i=1,ntimes write(6,"(' i=',i3,' ntimes=',i3,' day,hr,min=',3f7.2, | ' value=',e12.4)") i,ntimes,rindex(1:3,i),rindex(4,i) enddo call shutdown('set_index') ! ! Report to stdout: 100 continue ! write(6,"('set_index: ',a,' = ',e12.4)") name,outindex end subroutine set_index !----------------------------------------------------------------------- end module advance_module