module MusicBox_main !----------------------------------------------------------- ! use ccpp_kinds, only: r8 => kind_phys !----------------------------------------------------------- use ccpp_kinds, only: kind_phys use read_envConditions, only: read_envConditions_init, read_envConditions_timestep, read_envConditions_update_timestep use json_loader, only: json_loader_read use output_file, only: output_file_type !----------------------------------------------------------- ! MusicBox host model data !----------------------------------------------------------- use MusicBox_mod, only: box_press, box_temp, relhum, box_h2o, photo_lev, nspecies, vmr use MusicBox_mod, only: nbox, ntimes use MusicBox_mod, only: nkRxt, njRxt, TimeStart, TimeEnd use MusicBox_mod, only: nlayer, nlevel, zenith, albedo, press_mid, press_int use MusicBox_mod, only: alt, temp, o2vmrcol, o3vmrcol, so2vmrcol, no2vmrcol use MusicBox_mod, only: dt, density, mbar use MusicBox_mod, only: cnst_info use MusicBox_mod, only: jnames use MusicBox_mod, only: press_top use MusicBox_mod, only: cldfrc, cldwat implicit none public MusicBox_sub contains !> \section arg_table_MusicBox_sub Argument Table !! \htmlinclude arg_table_MusicBox_sub.html !! subroutine MusicBox_sub() !----------------------------------------------------------- ! Main driver routine for MusicBox - The box model of MICM !----------------------------------------------------------- use MusicBox_ccpp_cap, only: MusicBox_ccpp_physics_initialize use MusicBox_ccpp_cap, only: MusicBox_ccpp_physics_timestep_initial use MusicBox_ccpp_cap, only: MusicBox_ccpp_physics_run use MusicBox_ccpp_cap, only: MusicBox_ccpp_physics_timestep_final use MusicBox_ccpp_cap, only: MusicBox_ccpp_physics_finalize use MusicBox_ccpp_cap, only: ccpp_physics_suite_list use MusicBox_ccpp_cap, only: ccpp_physics_suite_part_list use tuv_photolysis, only: tuv_photolysis_readnl implicit none integer :: col_start, col_end integer :: index character(len=128), allocatable :: part_names(:) character(len=512) :: errmsg integer :: errflg integer :: i,n integer, allocatable :: PhotoLevs(:) ! photolysis level index for all boxes real(kind=kind_phys), allocatable :: vmrboxes(:,:) ! vmr for all boxes real(kind=kind_phys), allocatable :: RHboxes(:) ! RH for all boxes real(kind=kind_phys), allocatable :: Zenboxes(:) ! zenith angle for all boxes real(kind=kind_phys), allocatable :: Densityboxes(:) ! atm density for all boxes real(kind=kind_phys), allocatable :: Mbarboxes(:) ! atm mbar for all boxes type(output_file_type) :: outfile integer :: ibox integer :: thrdNdx real(kind_phys) :: sim_beg_time, sim_end_time ! run-time options character(len=120) :: env_conds_file = '../data/env_conditions.nc' character(len=120) :: outfile_name = 'test_output.nc' ! These need to be plain reals for the reading routine real, parameter :: NOT_SET = -huge(1.0) real, allocatable :: env_lat(:) real, allocatable :: env_lon(:) real, allocatable :: env_lev(:) real :: user_begin_time = NOT_SET ! seconds real :: user_end_time = NOT_SET real :: user_dtime = NOT_SET real :: timingStart, timingFinish, elapsedTime logical :: SolverDiags character(len=*), parameter :: photo_opts_file = '../Photolysis_options' character(len=*), parameter :: nml_options = '../MusicBox_options' character(len=120), parameter :: jsonfile = '../molec_info.json' integer :: omp_get_thread_num ! read namelist run-time options namelist /boxcnt/ nbox namelist /options/ outfile_name, env_conds_file namelist /options/ env_lat, env_lon, env_lev namelist /options/ user_begin_time, user_end_time, user_dtime !--------------------------- ! Read in the MusicBox_options file nbox = 1 open(unit=10,file=nml_options) read(unit=10,nml=boxcnt) if( nbox < 1 ) then write(*,*) 'MusicBox: box count = ',nbox,' is invalid; must be > 0' stop endif allocate(env_lat(nbox)) allocate(env_lon(nbox)) allocate(env_lev(nbox)) env_lat(:) = NOT_SET env_lon(:) = NOT_SET env_lev(:) = NOT_SET read(unit=10,nml=options) close(10) call tuv_photolysis_readnl(photo_opts_file, errmsg, errflg) if (errflg /= 0) then write(6, *) trim(errmsg) stop end if !--------------------------- ! error checking of namelist settings if (any(env_lat(:)<-90.) .or. any(env_lat(:) >90.)) then write(*,*) 'Invalid namelist setting: env_lat = ',env_lat(:) write(*,*) 'Must be set between -90 and 90 degrees north' stop end if if (any(env_lon(:)<0.) .or. any(env_lon(:)>360.)) then write(*,*) 'Invalid namelist setting: env_lon = ',env_lon(:) write(*,*) 'Must be set between 0 and 360 degrees east' stop end if if (any(env_lev(:)<0)) then write(*,*) 'Invalid namelist setting: env_lev = ',env_lev(:) write(*,*) 'Must be set to a positive pressure level (hPa)' stop end if !--------------------------- ! Read in the molecular information call json_loader_read( jsonfile, cnst_info, nSpecies, nkRxt, njRxt, jnames ) !--------------------------- ! Create the fields for the output netCDF file call outfile%create(nbox,outfile_name) call outfile%add(cnst_info) call outfile%add('Zenith','solar zenith angle','degrees') call outfile%add('Density','total number density','molecules/cm3') call outfile%add('Mbar','mean molar mass','g/mole') call outfile%add('RelHum','relative humidity','') call outfile%define() ! cannot add more fields after this call !--------------------------- ! Initialize the envrionmental conditions allocate(vmrboxes(nSpecies,nbox)) allocate(RHboxes(nbox),Zenboxes(nbox),Densityboxes(nbox),Mbarboxes(nbox)) allocate(PhotoLevs(nbox)) call read_envConditions_init(nbox, nSpecies, env_conds_file, env_lat, env_lon, env_lev, & user_begin_time, user_end_time, user_dtime, cnst_info, & vmrboxes, dt, sim_beg_time, sim_end_time, nlayer, PhotoLevs) !--------------------------- ! Set up the various dimensions nlevel = nlayer+1 ! number if vertical interface levels !--------------------------- ! allocate host model arrays allocate(vmr(nSpecies)) allocate(alt(nlayer)) allocate(press_mid(nlayer)) allocate(press_int(nlevel)) allocate(temp(nlayer)) allocate(o2vmrcol(nlayer)) allocate(o3vmrcol(nlayer)) allocate(so2vmrcol(nlayer)) allocate(no2vmrcol(nlayer)) allocate(cldwat(nlayer)) allocate(cldfrc(nlayer)) !--------------------------- ! Set the times (note this needs to be set prior to call ccpp_initialize) ! Once Rosenbrock_init is separated into init and time_step_init, this may go ! down right above time_step_init TimeStart = sim_beg_time TimeEnd = TimeStart + dt !--------------------------- ! Use the suite information to setup the run call MusicBox_ccpp_physics_initialize( & 'MusicBox_suite', & box_press, box_temp, relhum, box_h2o, zenith, & albedo, density, photo_lev, alt, temp, cldfrc, & cldwat, o2vmrcol, o3vmrcol, so2vmrcol, no2vmrcol, & press_mid, press_int, vmr, mbar, SolverDiags, & errmsg, errflg) if (errflg /= 0) then write(6, *) trim(errmsg) stop end if ElapsedTime = 0. !----------------------------------------------------------- ! loop over time !----------------------------------------------------------- time_loop: do while (timestart <= sim_end_time) if( timestart == 66150._kind_phys ) then write(*,*) 'At debug time' endif !--------------------------- ! Initialize the timestep call MusicBox_ccpp_physics_timestep_initial( & 'MusicBox_suite', & box_press, box_temp, relhum, box_h2o, zenith, & albedo, density, photo_lev, alt, temp, cldfrc, & cldwat, o2vmrcol, o3vmrcol, so2vmrcol, no2vmrcol, & press_mid, press_int, vmr, mbar, SolverDiags, & errmsg, errflg) if (errflg /= 0) then write(6, *) trim(errmsg) stop end if !--------------------------- ! set the timestep for the output file call outfile%advance(TimeStart) call cpu_time(timingStart) !--------------------------- ! Loop over the boxes !$OMP PARALLEL DO PRIVATE(ibox,box_press,box_temp,relhum,box_h2o,zenith,albedo, & !$OMP density,photo_lev,press_top,alt,temp,cldfrc,cldwat, & !$OMP o2vmrcol,o3vmrcol,so2vmrcol,no2vmrcol, & !$OMP press_mid,press_int,vmr,mbar,SolverDiags) Box_loop: do ibox=1,nbox photo_lev = PhotoLevs(ibox) SolverDiags = ibox == 1 .and. TimeStart == sim_beg_time !--------------------------- ! Read in the environmental conditions at TimeStart call read_envConditions_timestep(TimeStart,ibox, nlayer, photo_lev, vmrboxes, zenith, albedo, & press_mid, press_int, alt, temp, o2vmrcol, o3vmrcol, so2vmrcol, no2vmrcol, vmr, box_h2o, & box_temp, box_press) cldwat(:) = 0._kind_phys cldfrc(:) = 0._kind_phys press_top = press_int(1) !--------------------------- ! Call the schemes for the timestep col_start=1 col_end=1 call MusicBox_ccpp_physics_run( & 'MusicBox_suite', 'physics', & box_press, box_temp, relhum, box_h2o, zenith, & albedo, density, photo_lev, alt, temp, cldfrc, & cldwat, o2vmrcol, o3vmrcol, so2vmrcol, no2vmrcol, & press_mid, press_int, vmr, mbar, SolverDiags, & col_start, col_end, errmsg, errflg) if (errflg /= 0) then write(6, *) trim(errmsg) call ccpp_physics_suite_part_list('MusicBox_suite', part_names, errmsg, errflg) write(6, *) 'Available suite parts are:' do index = 1, size(part_names) write(6, *) trim(part_names(index)) end do stop end if !--------------------------- ! Update the environmental conditions for the timestep call read_envConditions_update_timestep(ibox, vmr, vmrboxes) RHboxes(ibox) = RelHum Zenboxes(ibox) = zenith Densityboxes(ibox) = density Mbarboxes(ibox) = mbar !--------------------------- ! output the timestep variables thrdNdx = omp_get_thread_num() write(*,'(a,1p,g15.7,1x,g0,1x,g0)') ' total density, pressure, temperature :', density, box_press, box_temp write(*,'(a,i3,1x,i3,1x,1p,g0)') 'Concentration @ thrd,box,hour = ',thrdNdx,ibox,TimeStart/3600. write(*,'(1p,5(1x,g0))') vmrboxes(:,ibox),sum(vmrboxes(:,ibox)) end do Box_loop !$OMP END PARALLEL DO call cpu_time(timingFinish) ElapsedTime = ElapsedTime + (timingFinish - timingStart) !--------------------------- ! write out the timestep values call outfile%out( 'RelHum', RHboxes ) call outfile%out( 'Zenith', Zenboxes ) call outfile%out( 'Density', Densityboxes ) call outfile%out( 'Mbar', Mbarboxes ) call outfile%out( cnst_info, vmrboxes ) !--------------------------- ! Advance the timestep TimeStart = TimeEnd TimeEnd = TimeStart + dt !--------------------------- ! call the timestep_final scheme routines call MusicBox_ccpp_physics_timestep_final( & 'MusicBox_suite', & box_press, box_temp, relhum, box_h2o, zenith, & albedo, density, photo_lev, alt, temp, cldfrc, & cldwat, o2vmrcol, o3vmrcol, so2vmrcol, no2vmrcol, & press_mid, press_int, vmr, mbar, SolverDiags, & errmsg, errflg) if (errflg /= 0) then write(6, *) trim(errmsg) write(6,'(a)') 'An error occurred in ccpp_timestep_final, Exiting...' stop end if end do time_loop !--------------------------- ! Finalize all of the schemes call MusicBox_ccpp_physics_finalize( & 'MusicBox_suite', & box_press, box_temp, relhum, box_h2o, zenith, & albedo, density, photo_lev, alt, temp, cldfrc, & cldwat, o2vmrcol, o3vmrcol, so2vmrcol, no2vmrcol, & press_mid, press_int, vmr, mbar, SolverDiags, & errmsg, errflg) if (errflg /= 0) then write(6, *) trim(errmsg) write(6,'(a)') 'An error occurred in ccpp_finalize, Exiting...' stop end if call outfile%close() deallocate(vmr) deallocate(vmrboxes) deallocate(RHboxes,Zenboxes,Densityboxes,Mbarboxes) deallocate(PhotoLevs) deallocate(alt) deallocate(press_mid) deallocate(press_int) deallocate(temp) deallocate(o2vmrcol) deallocate(o3vmrcol) deallocate(so2vmrcol) deallocate(no2vmrcol) deallocate(cldwat) deallocate(cldfrc) write(*,*) ' ' write(*,*) 'Elapsed box loop time = ',ElapsedTime end subroutine MusicBox_sub end module MusicBox_main !> \brief Main SCM program that calls the main SCM subroutine !! !! The Doxygen documentation system cannot handle in-body comments in Fortran main programs, so the "main" program was put in the !! subroutine \ref MusicBox_main_sub above. program MusicBox use MusicBox_main call MusicBox_sub() end program MusicBox