! This is the main program of weps driver for air quality ! program weps_driver USE pd_dates_vars USE pd_update_vars USE pd_report_vars USE pd_var_tables USE mandate_vars ! Load shared mandates() array use omp_lib ! load openMP lib include 'p1werm.inc' include 'wpath.inc' include 'p1unconv.inc' include 'm1subr.inc' include 'm1sim.inc' include 'm1geo.inc' include 'm1flag.inc' include 'm1dbug.inc' include 's1layr.inc' include 's1sgeo.inc' include 's1phys.inc' include 's1agg.inc' include 'c1info.inc' include 'c1gen.inc' include 'w1clig.inc' include 'w1wind.inc' include 'w1pavg.inc' include 'file.fi' include 'b1glob.inc' include 'd1glob.inc' include 'c1glob.inc' include 'c1db1.inc' include 'h1hydro.inc' include 'timer.fi' include 'command.inc' !declarations for commandline args include 'precision.inc' !declaration for portable math range checking ! + + + LOCAL COMMON BLOCKS + + + include 'main/main.inc' include 'manage/man.inc' include 'manage/oper.inc' include 'decomp/decomp.inc' include 'erosion/p1erode.inc' !Needs the SURF_UPD_FLG variable include 'erosion/m2geo.inc' !Need tsterode cmdline arg vars(xgdpt,ygdpt) include 'erosion/e2erod.inc' include 'airpact/emit.inc' include 'airpact/spatialGIS.inc' include 'hydro/dvolwparam.inc' ! define local variables integer cd,cm,cy integer julday integer cnt,sd,sm,sy,ed,em,ey integer stjday,endjday real ci integer TID ! thread ID ! local variables integer get_nperiods integer nperiods integer i,j,idx,d,k integer r,s integer m,n,nrow,ncln,srow,scln logical flag_idx(landc_type,soil_comp) ! 8 is manage tpyes and 200 is soil types ! logical initial_flag(landc_type,soil_comp) ! flag for initial ! temp erosion variable to save subregion grid erosion real total_tmp(nday,landc_type,soil_comp) real ssp_tmp(nday,landc_type,soil_comp) real pm10_tmp(nday,landc_type,soil_comp) ! temp erosion variable to save hourly subregion grid erosion real, allocatable ::total_htmp(:,:,:,:) real, allocatable ::ssp_htmp(:,:,:,:) real, allocatable ::pm10_htmp(:,:,:,:) ! This is for 1-km grid hourly output real, allocatable ::total_1kh(:,:,:,:,:,:) real, allocatable ::ssp_1kh(:,:,:,:,:,:) real, allocatable ::pm10_1kh(:,:,:,:,:,:) ! the temporial variables to save the hourly output real :: ehtotal(24,nday,n_grid,n_grid) real :: ehssp(24,nday,n_grid,n_grid) real :: ehpm10(24,nday,n_grid,n_grid) integer*4 start_date(3), start_time(3),end_date(3),end_time(3) allocate (total_htmp(24,nday,landc_type,soil_comp)) allocate (ssp_htmp(24,nday,landc_type,soil_comp)) allocate (pm10_htmp(24,nday,landc_type,soil_comp)) allocate (total_1kh(24,nday,n_grid,n_grid,12,12)) allocate (ssp_1kh(24,nday,n_grid,n_grid,12,12)) allocate (pm10_1kh(24,nday,n_grid,n_grid,12,12)) ! write the cligen and wingen variables into a text file open(4,file="MCIP.txt") id = 1 ! initialization dates, time im = 1 iy = 0006 ld = 24 lm = 08 ly = 0009 sd = 25 !simulation dates, time sm = 08 sy = 0009 ed = 25 em = 08 ey = 0009 call idate(start_date) call itime(start_time) srow = 17 scln = 1 nrow = 34 ncln = 50 max_real = huge(1.0) * 0.999150 max_arg_exp = log(max_real) SURF_UPD_FLG = 1 xgdpt = 0 !use default grid spacing values if ygdpt = 0 !these are not specified on the commandline erod_interval = 0 !default value for updating eroding soil surface !(currently only used in standalone erosion submodel) calib_cycle = 0 max_calib_cycles = 3 ! Default value unless increased via cmdline option calib_done = .false. wc_type = 0 ! default of water content type ntstep = 24 nsubr= 1 nacctr = 1 SoilRockFragments(1) = -1 ! Setting default value to -1 (single subregion only for now!!!) ci = 0.90 ! default confidence interval value ! temporarily initialize old random roughness ! aslrrc(1) = 10. ! as0rrk(1) = 0.9 run_erosion = 1 init_cycle = 1 rootp ='../test1/' wepp_hydro = 1 ! This is a test of darcy function by Jin report_loop = .false. ! if (btest(am0efl,0)) then ! write(*,*) 'Luo_erod',luo_erod call fopenk(luo_erod, rootp(1:len_trim(rootp))// & & 'output/daily_erosion.out','unknown') call fopenk(luo_emit, rootp(1:len_trim(rootp))// & & 'output/erosion_hour.out','unknown') call fopenk(50,rootp(1:len_trim(rootp))//'output/EGT.out', & & 'unknown') call fopenk(55,rootp(1:len_trim(rootp))//'output/EGT_d.out', & & 'unknown') call fopenk(53,rootp(1:len_trim(rootp))//'output/SSP.out', & & 'unknown') call fopenk(66,rootp(1:len_trim(rootp))//'output/SSP_d.out', & & 'unknown') call fopenk(52,rootp(1:len_trim(rootp))//'output/PM10.out', & & 'unknown') call fopenk(67,rootp(1:len_trim(rootp))//'output/PM10_d.out', & & 'unknown') call fopenk(76,rootp(1:len_trim(rootp))//'output/total_1kh.out', & & 'unknown') call fopenk(77,rootp(1:len_trim(rootp))//'output/ssp_1kh.out', & & 'unknown') call fopenk(78,rootp(1:len_trim(rootp))//'output/PM10_1kh.out', & & 'unknown') ! endif if (btest(am0efl,1)) then call fopenk (luo_egrd,rootp//'output/daily_egrd.out', 'unknown') endif if (calibrate_crops > 3) max_calib_cycles = calibrate_crops ijday = julday(id, im, iy) ljday = julday(ld, lm, ly) ! call cliginit() ! open clig file ! call erodinit() cnt = 0 call read_soilGIS() close (47) call read_landcov() call read_cligen_wingenGIS() call cmzlist_read() call mcip_weps ! call from mcip_weps.f95 ! initial flag index and temp erosion variables ! For 10 landcov types and 200 soil components do i=1,landc_type do j=1,soil_comp flag_idx(i,j) = .false. ! initial_flag(i,j) = .false. do d=1,nday total_tmp(d,i,j) = 0.0 ssp_tmp(d,i,j) = 0.0 pm10_tmp(d,i,j) = 0.0 do k=1,24 total_htmp(k,d,i,j) = 0.0 ssp_htmp(k,d,i,j) = 0.0 pm10_htmp(k,d,i,j) = 0.0 do m=1,12 do n=1,12 total_1kh(k,d,i,j,m,n) = 0.0 ssp_1kh(k,d,i,j,m,n) = 0.0 pm10_1kh(k,d,i,j,m,n) = 0.0 end do end do end do end do end do end do ! create an index array for each element in 12x12 km grid ! isr = rdx(i,j) idx = 1 ! soilDIM(i,j) is the unique soil ID ! Assign a temp index for each subregion grid do i=1,12 !row do j=1,12 !col rdx(i,j) = idx idx = idx+1 end do end do cnt = 0 ! run the whole region with 95x95 12-km grids and 12x12 1-km grids (34,50) !!!! parameterize these indices here. !$omp PARALLEL PRIVATE(TID) TID = OMP_GET_THREAD_NUM() do m=srow,nrow do n=scln,ncln ! initialize the erosion daily and hourly output !change day to nday do d=1,nday etotal(d,m,n) = 0.0 essp(d,m,n) = 0.0 epm10(d,m,n) = 0.0 do k=1,24 ehtotal(k,d,m,n) = 0.0 ehssp(k,d,m,n) = 0.0 ehpm10(k,d,m,n) = 0.0 end do end do ! end initialize ! set all flag_idx as false at the beginning of model running at a sub-grid ! in a 12-km grid ! turn off this initial as we use only one climatic data on 5_26_10 !do i=12,12 !do j=12,12 do i=1,12 do j=1,12 if (landcov(m,n,i,j) .gt. 0 .AND. rdx(i,j) .gt. 0 .and. & & soilDim(m,n,i,j) .gt. 0 ) then call get_cmz_idx_name(landcovIDX(m,n,i,j), & & landcov(m,n,i,j)) flag_idx(landcov(m,n,i,j),soilDim(m,n,i,j)) = .FALSE. end if end do end do ! run each subgrid with 1-km grid do i = 1,12 do j = 1,12 ! chang to 1 for first column run if (landcov(m,n,i,j) .gt. 0 .AND. rdx(i,j) .gt. 0 .and. & & soilDim(m,n,i,j) .gt. 0 ) then if (flag_idx(landcov(m,n,i,j),soilDim(m,n,i,j)) & & .neqv. .TRUE.) then call get_cmz_idx_name(landcovIDX(m,n,i,j), & & landcov(m,n,i,j)) call read_manage1(landcov(m,n,i,j)) ! loading manage file for this grid call read_soil(m,n,i,j,rdx(i,j)) !loading soil features call open_cli_win(rdx(i,j),m,n) ! Running the weps model daily ! add day here ! If layer is 0, we will not run the weps model if ( nslay(rdx(i,j)) .gt. 0) then tinfil = rootp(1:len_trim(rootp))//'manage/'//tinfil call mfinit(rdx(i,j), tinfil, maxper) write(*,*) 'tinfil:',tinfil ! check for consistency maxper, n_rot_cycles, number of years to run if( maxper*n_rot_cycles .ne. ly-iy+1 ) then write(*,*) 'Warning: Number of rotations (',n_rot_cycles,') ',& & 'times Years in rotation (',maxper,') ', & & ' does not match Number of simulation years (', & & ly-iy+1,') ' n_rot_cycles = (ly-iy+1) / maxper if( mod( (ly-iy+1), maxper ) .gt. 0 ) then write(*,*) & & 'Warning: Not simulating complete rotations' n_rot_cycles = n_rot_cycles + 1 end if end if if (calibrate_rotcycles .eq. 0) then calibrate_rotcycles = n_rot_cycles end if call mandates(rdx(i,j)) !Get man dates, op names, and crop names nperiods = get_nperiods(maxper) !Get # of periods for reports if( report_debug >= 1 ) then write(*,*) '# rot years', maxper, "nperiods", & & nperiods, '# cycles', n_rot_cycles end if call weps_init(rdx(i,j)) ! calculate first and last Julian dates for simulation (initialization) ijday = julday(id, im, iy) ljday = julday(ld, lm, ly) ! calculate the actually Julian dates for simulation stjday = julday(sd,sm,sy) endjday = julday(ed,em,ey) ! calculate last julian date for initialization cycle am0csr = 1 ! set global current subregion variable? if (init_cycle > 0) then ! to avoid printing it when not being done write(6,*) "Starting initialization phase" else lopyr = 1 endif ! init_file ='N'//rdx(i,j)//'_'//m//n ! begin initialization simulation phase init_loop = .true. ! Signifies that we are in the "initialization" loop ! Running the initialization for a whole period of rotation ! write(*,*) 'Running initialization ... ' ! if (initial_flag(landcov(m,n,i,j),soilDim(m,n,i,j)).eq. .false.) & ! & then do am0jd = ijday,ljday call caldatw (cd, cm, cy) call getcli(cd,cm,cy) call getwin(cd,cm,cy) call submodels(rdx(i,j),cd,cm,cy) rewind (luicli) rewind (luiwin) end do write(*,*) 'Finish initialization!' call plotstate3(m,n,rdx(i,j),1,1,1) ! save the initial variables for later running ! initial_flag(landcov(m,n,i,j),soilDim(m,n,i,j)) = .TRUE. ! end if ! running the date loop d = 0 do am0jd = stjday, endjday d = d+1 call caldatw (cd, cm, cy) ! call cligen and wingen data set here call getcli_win(n+6,83-m,d) ! write(4,*) n+6,83-m,cd,cm,cy,d,awtdmx,awtdmn,awzdpt, & ! & aweirr,awadir,(awu(k),k=1,24) if (am0jd .eq. stjday) then call inpSaveFile0(m,n,rdx(i,j),1,1,1) call submodels(rdx(i,j),cd,cm,cy) amnrotcycle(rdx(i,j)) = 1 ! set here for use in confidence interval calculation (no other use?) am0sif = .false. ! Done with all initialization and calibration phases report_loop = .true. else call caldat(am0jd-1,cd,cm,cy) ! call the data from previous day call inpSaveFile(m,n,rdx(i,j),cd,cm,cy) call caldat(am0jd,cd,cm,cy) call submodels(rdx(i,j),cd,cm,cy) end if call erodinit(rdx(i,j)) if (run_erosion > 0) then if (awudmx .gt. 5.0) then call calcwu call erosion(rdx(i,j),5.0) ! write(*,*)'Erosion Running:',egt(1,1),egtss(1,1),egt10(1,1) !if (egt(1,1) .ne. 0) then ! write(luo_erod,*)cd,cm,cy,m,n,i,j,egt(1,1),egtss(1,1),egt10(1,1) !end if !if (btest(am0efl,0) .or. btest(am0efl,1)) then call daily_erodout (luo_egrd, luo_erod) ! write(*,*) 'done daily on date!',cm,cd,cy !endif end if end if !if (egt(1,rdx(i,j)).gt. 0 .or. egtss(1,rdx(i,j)) .gt. 0 .or. & !& egt10(1,rdx(i,j)) .gt. 0) then ! write(luo_egrd,*)cm, cd,cy,egt(1,rdx(i,j)),egtss(1,rdx(i,j)), & !& egt10(1,rdx(i,j)) !end if ! call plotstate(m,n,rdx(i,j),cd,cm,cy) ! save into a text file for testing call plotstate2(m,n,rdx(i,j),cd,cm,cy) ! write into a temp file which will be !called by inpSaveFile() !if (egtss(1,1) .NE. 0) then ! write(*,*) 'EROSION:',isr,cd,cm,cy,egtss(1,1),egt(1,1),egt10(1,1) !end if ! add the total erosion from subregions do r=1,imax do s=1,jmax total_tmp(d,landcov(m,n,i,j),soilDim(m,n,i,j)) & & = total_tmp(d,landcov(m,n,i,j),soilDim(m,n,i,j)) + egt(r,s) ssp_tmp(d,landcov(m,n,i,j),soilDim(m,n,i,j)) & & = ssp_tmp(d,landcov(m,n,i,j),soilDim(m,n,i,j)) + egtss(r,s) pm10_tmp(d,landcov(m,n,i,j),soilDim(m,n,i,j)) & & = pm10_tmp(d,landcov(m,n,i,j),soilDim(m,n,i,j)) + egt10(r,s) end do end do ! total_tmp(d,landcov(m,n,i,j),soilDim(m,n,i,j)) = egt(1,1) ! ssp_tmp(d,landcov(m,n,i,j),soilDim(m,n,i,j)) = egtss(1,1) ! pm10_tmp(d,landcov(m,n,i,j),soilDim(m,n,i,j)) = egt10(1,1) etotal(d,m,n) = etotal(d,m,n) + total_tmp(d,landcov(m,n,i,j), soilDim(m,n,i,j)) essp(d,m,n) = essp(d,m,n) + ssp_tmp(d,landcov(m,n,i,j), soilDim(m,n,i,j)) epm10(d,m,n) = epm10(d,m,n) + pm10_tmp(d,landcov(m,n,i,j), soilDim(m,n,i,j)) ! if (etotal(d,m,n) .GT. 0) then ! write(luo_erod,*)cd,cm,cy,m,n,i,j,etotal(d,m,n),essp(d,m,n),epm10(d,m,n) ! write(luo_erod,*)'ahgt:',(ahgt(k),k=1,24) ! write(luo_erod,*)'ahssp:',(ahssp(k),k=1,24) ! write(luo_erod,*)'ahpm10:',(ahpm10(k),k=1,24) ! write(luo_erod,*)' ' ! end if ! assign the hourly spatial output do k=1,24 total_htmp(k,d,landcov(m,n,i,j),soilDim(m,n,i,j)) = ahgt(k) ssp_htmp(k,d,landcov(m,n,i,j),soilDim(m,n,i,j)) = ahssp(k) pm10_htmp(k,d,landcov(m,n,i,j),soilDim(m,n,i,j)) = ahpm10(k) ehtotal(k,d,m,n) = ehtotal(k,d,m,n) + ahgt(k) ehssp(k,d,m,n) = ehssp(k,d,m,n) + ahssp(k) ehpm10(k,d,m,n) = ehpm10(k,d,m,n) + ahpm10(k) ! assign the hourly 1-km grid vale here total_1kh(k,d,m,n,i,j) = ahgt(k) ssp_1kh(k,d,m,n,i,j) = ahssp(k) pm10_1kh(k,d,m,n,i,j) = ahpm10(k) end do end do ! end for day loop deallocate (mandate) ! this should be moved inside "mandates" ! call weps_run(m,n,rdx(i,j)) cnt = cnt+1 end if ! end for nslay if rewind (luicli) ! call save_erosion(rdx(i,j)) ! adding codes for this function flag_idx(landcov(m,n,i,j),soilDim(m,n,i,j)) = .TRUE. else ! Modify egt from 2-D into 3-D by adding a time dimension ??? ??? ! change day to nday do d=1,nday ! egt_tp(d)=total_tmp(d,landcov(m,n,i,j),soilDim(m,n,i,j)) ! ssp_tp(d)=ssp_tmp(d,landcov(m,n,i,j),soilDim(m,n,i,j)) ! pm10_tp(d)=pm10_tmp(d,landcov(m,n,i,j),soilDim(m,n,i,j)) etotal(d,m,n) = etotal(d,m,n)+total_tmp(d,landcov(m,n,i,j), soilDim(m,n,i,j)) essp(d,m,n) = essp(d,m,n)+ssp_tmp(d,landcov(m,n,i,j), soilDim(m,n,i,j)) epm10(d,m,n) = epm10(d,m,n)+pm10_tmp(d,landcov(m,n,i,j), soilDim(m,n,i,j)) do k=1,24 ehtotal(k,d,m,n)=ehtotal(k,d,m,n) + total_htmp(k,d,landcov(m,n,i,j),soilDim(m,n,i,j)) ehssp(k,d,m,n) = ehssp(k,d,m,n) + ssp_htmp(k,d,landcov(m,n,i,j),soilDim(m,n,i,j)) ehpm10(k,d,m,n) = ehpm10(k,d,m,n) + pm10_htmp(k,d,landcov(m,n,i,j),soilDim(m,n,i,j)) ! assign the hourly 1-km grid erosion value total_1kh(k,d,m,n,i,j) = total_htmp(k,d,landcov(m,n,i,j),soilDim(m,n,i,j)) ssp_1kh(k,d,m,n,i,j) = ssp_htmp(k,d,landcov(m,n,i,j),soilDim(m,n,i,j)) pm10_1kh(k,d,m,n,i,j) = pm10_htmp(k,d,landcov(m,n,i,j),soilDim(m,n,i,j)) end do end do ! write(*,*) 'The same soil and landcov is already run!' ! get the erosion from the same landcov and soil end if ! add the erosion together in a grid !etotal(d,m,n) = etotal(d,m,n)+egt_tp(d) !essp(d,m,n) = essp(d,m,n)+ssp_tp(d) !epm10(d,m,n) = epm10(d,m,n)+pm10_tp(d) end if end do ! end of 1-km grid loop end do ! end of 1-km grid loop ! write(luo_erod,*)'Eros:',m,n,etotal(m,n),essp(m,n),epm10(m,n) end do ! end of 12-km grid loop end do ! end of 12-km grid loop !OMP END PARALLEL ! writing the out put into a text file write(*,*) 'Total model running is:',cnt do d=1, nday write(55,*) 'Date:',d write(66,*) 'Date:',d write(67,*) 'Date:',d do i=srow,nrow write(55,1010)(etotal(d,i,j)*1000, j=scln,ncln) write(66,1010)(essp(d,i,j)*1000, j=scln,ncln) write(67,1010)(epm10(d,i,j)*1000, j=scln,ncln) end do do k=1,24 write(50,*) 'Date and Hour:',d,k write(53,*) 'Date and Hour:',d,k write(52,*) 'Date and Hour:',d,k do i=srow,nrow write(50,1010)(ehtotal(k,d,i,j)*1000, j=scln,ncln) write(53,1010)(ehssp(k,d,i,j)*1000, j=scln,ncln) write(52,1010)(ehpm10(k,d,i,j)*1000, j=scln,ncln) do j=scln,ncln do m = 1,12 write(76,1020),i,j,m,(total_1kh(k,d,i,j,m,n)*1000,n=1,12) write(77,1020),i,j,m,(ssp_1kh(k,d,i,j,m,n)*1000,n=1,12) write(78,1020),i,j,m,(pm10_1kh(k,d,i,j,m,n)*1000,n=1,12) end do end do end do end do end do ! write(*,*) 'Total model running is:',cnt close (luiwin) close (luicli) close (luo_egrd) close (luo_erod) close (luo_emit) close (luimandate) close (50) close (53) close (52) close (42) close (4) close (55) close (66) close (67) close (76) close (77) close (78) call idate(end_date) call itime(end_time) write(*,1000)start_date(2),start_date(1),start_date(3),start_time write(*,1000)end_date(2),end_date(1),end_date(3),end_time 1000 format('D:',i2.2,'/',i2.2,'/',i4.4,';T:',i2.2,':',i2.2,':',i2.2) 1010 format (50(f13.5,' ')) 1020 format (3(i2,','),12(f13.5,1X)) ! end of subregions luo_erod ! Done with simulation here .................. end