!tks - program to put raqms and wrf data on 100 layer atmosphere. wrf is used where it can be ! raqms fills the atmosphere above the wrf data. wrf data has 40 layers, raqms has 35 USE Type_Kinds INTEGER, PARAMETER :: nk100=100,nk_wrf=40, nk_raqms=35, nk_use_raqms=20,nchem=15, nmet=2, & nwrfmet=5,n_sections=18 CHARACTER (140) :: gfile,output_dir,output_file,div1*1,div2*2 INTEGER :: n_level,n_layer,n_aerosol,n_profile,i_prof,kk !tks here fp_kind is real*8, fp_kind4 is r*4 integer (kind=4) :: nxx,nyy real (fp_kind4),allocatable, dimension(:,:) :: zsfc,u10,v10,slmsk,landuse,icemask, & tsea,t2,smois,vegfrac,st1,snow,ps real(fp_kind4),allocatable, dimension(:,:,:) :: p3d_wrf, t3d_wrf,q3d_wrf,cld3d_wrf,cldice3d_wrf,o3d_wrf real(fp_kind4),allocatable, dimension(:,:,:,:) :: chem_wrf real(fp_kind4),allocatable, dimension(:,:) :: lats4,lons4,zenith4,azimuth,solzen real(fp_kind4),allocatable, dimension(:,:,:) :: p3d_raqms,t3d_raqms,o3d_raqms real(fp_kind) t100(nk100),q100(nk100),o3100(nk100),cld100(nk100),ch100(nk100,nchem),cldice100(nk100) INTEGER :: stype,sflag,iclouds,icld,npnts,nll,kr,j_loop,i_loop,ichem,it REAL( fp_kind) :: smc,sheleg,soilt,vfrac,prs,pu,ooo,xx,ttt,ppp,qqq,pd,pm REAL( fp_kind ) :: a, a1, a2, alat, alon, height, scover, u, v, total_clw, total_aerosol REAL( fp_kind ) :: scanang,satz,slat,slon,latmax,latmin,lonmax,lonmin,den REAL( fp_kind ) :: earthradius,sathght,solazi,satzen,satazi,aod,szg INTEGER :: Absorber_ID(10), Absorber_UNITS(10) INTEGER :: n_Layers, n_Absorbers, n_Clouds, Climatology real play(nk100),prlvl(nk100+1) character pfile*140,hfile*140,infile*140,outfile*140,ch3*3,cem(nchem)*6 character raqms_3d_met_lab(nmet)*12,odir*140,wrf_3d_chem_map(nchem)*20,wrf_3d_met_lab(nwrfmet)*20 data prlvl/0.005, 0.016, 0.038, 0.077, 0.137, 0.224, 0.345, 0.506, & 0.714, 0.975, 1.297, 1.687, 2.153, 2.701, 3.340, 4.077, & 4.920, 5.878, 6.957, 8.165, 9.512, 11.004, 12.649, 14.456, & 16.432, 18.585, 20.922, 23.453, 26.183, 29.121, 32.274, 35.650, & 39.257, 43.100, 47.188, 51.528, 56.126, 60.990, 66.125, 71.540, & 77.240, 83.231, 89.520, 96.114, 103.017, 110.237, 117.777, 125.646, & 133.846, 142.385, 151.266, 160.496, 170.078, 180.018, 190.320, 200.989, & 212.028, 223.441, 235.234, 247.409, 259.969, 272.919, 286.262, 300.000, & 314.137, 328.675, 343.618, 358.967, 374.724, 390.893, 407.474, 424.470, & 441.882, 459.712, 477.961, 496.630, 515.720, 535.232, 555.167, 575.525, & 596.306, 617.511, 639.140, 661.192, 683.667, 706.565, 729.886, 753.627, & 777.790, 802.371, 827.371, 852.788, 878.620, 904.866, 931.524, 958.591, & 986.067,1013.948,1042.232,1070.917,1100.000/ data play/0.009, 0.026, 0.055, 0.104, 0.177, 0.281, 0.421, 0.604, & 0.838, 1.129, 1.484, 1.910, 2.416, 3.009, 3.696, 4.485, & 5.385, 6.402, 7.545, 8.822, 10.240, 11.807, 13.532, 15.423, & 17.486, 19.730, 22.163, 24.793, 27.626, 30.671, 33.934, 37.425, & 41.148, 45.113, 49.326, 53.794, 58.524, 63.523, 68.797, 74.353, & 80.198, 86.338, 92.778, 99.526, 106.586, 113.965, 121.669, 129.703, & 138.072, 146.781, 155.836, 165.241, 175.001, 185.121, 195.606, 206.459, & 217.685, 229.287, 241.270, 253.637, 266.392, 279.537, 293.077, 307.014, & 321.351, 336.091, 351.236, 366.789, 382.751, 399.126, 415.914, 433.118, & 450.738, 468.777, 487.236, 506.115, 525.416, 545.139, 565.285, 585.854, & 606.847, 628.263, 650.104, 672.367, 695.054, 718.163, 741.693, 765.645, & 790.017, 814.807, 840.016, 865.640, 891.679, 918.130, 944.993, 972.264, & 999.942,1028.025,1056.510,1085.394/ data wrf_3d_chem_map/ 'o3','sulf', 'BC1','BC2','OC1','OC2','DUST_1', & 'DUST_2','DUST_3', 'DUST_4','DUST_5','SEAS_1','SEAS_2','SEAS_3','SEAS_4'/ data raqms_3d_met_lab/'Pressure','Temperature'/ data wrf_3d_met_lab/'Pressure','Temperature','WaterVapor','LiquidWaterContent','IceWaterContent'/ data cem/'o3vmr','iso4aer','ibc1','ibc2','ioc1','ioc2','idu1','idu2','idu3','idu4', & 'idu5','iss1','iss2','iss3','iss4'/ 1001 format(i1) 1002 format(i2) ! LOOP OVER NUMBER OF SECTIONS do 212 it=1,n_sections if(it.lt.10)then write(div1,1001)it div2='0'//div1 else write(div2,1002)it endif !__________READ in DATA____________ !Open Input and output files !odir -directory where input data is located odir='/sierra/raqms_data/Model/RAQMS-WRF-CHEM/Proxy/' output_dir='./' !location of output output_file=trim(output_dir)//'tom_test'//div2//'.dat' !_________________________________ !Read in surface, RAQMS and WRF datasets ! read WRF/Chem 2d sfc data gfile=trim(odir)//'wrf-chem.2d.2006.08.24.06Z.W'//div2//'.dat' open(3,file=gfile,form='unformatted',convert="little_endian") read(3)nxx,nyy !allocate arrays based on nxx,nyy allocate(lats4(nxx,nyy),lons4(nxx,nyy),zenith4(nxx,nyy),azimuth(nxx,nyy),solzen(nxx,nyy)) allocate(zsfc(nxx,nyy),u10(nxx,nyy),v10(nxx,nyy),landuse(nxx,nyy),slmsk(nxx,nyy),icemask(nxx,nyy), & tsea(nxx,nyy),t2(nxx,nyy),smois(nxx,nyy),vegfrac(nxx,nyy),st1(nxx,nyy),snow(nxx,nyy),ps(nxx,nyy)) allocate(p3d_wrf(nxx,nyy,nk_wrf),t3d_wrf(nxx,nyy,nk_wrf),q3d_wrf(nxx,nyy,nk_wrf), & o3d_wrf(nxx,nyy,nk_wrf),cld3d_wrf(nxx,nyy,nk_wrf),cldice3d_wrf(nxx,nyy,nk_wrf)) read(3)lats4,lons4,zenith4,azimuth,solzen read(3)zsfc,u10,v10,landuse,slmsk,icemask,tsea,t2,smois,vegfrac,st1,snow,ps close(3) ! read WRF 3d met do K=1,nwrfmet gfile=trim(odir)//'wrf-chem.'//trim(wrf_3d_met_lab(k))//'.2006.08.24.06Z.W'//div2//'.dat' open(3,file=gfile,form='unformatted',convert="little_endian") read(3)nxx,nyy read(3)lats4,lons4,zenith4,azimuth,solzen read(3)nxx,nyy,nll if(wrf_3d_met_lab(k).eq.'Pressure')read(3)p3d_wrf if(wrf_3d_met_lab(k).eq.'Temperature')read(3)t3d_wrf if(wrf_3d_met_lab(k).eq.'WaterVapor')read(3)q3d_wrf if(wrf_3d_met_lab(k).eq.'LiquidWaterContent')read(3)cld3d_wrf if(wrf_3d_met_lab(k).eq.'IceWaterContent')read(3)cldice3d_wrf close(3) end do print*,'Read WRFf 3d MET ',nxx,nyy,nll,nk_wrf allocate(chem_wrf(nxx,nyy,nk_wrf,nchem)) ! read WRF/Chem 3d chems !k=1 is ozone do K=1,nchem gfile=trim(odir)//'wrf-chem.'//trim(wrf_3d_chem_map(k))//'.2006.08.24.06Z.W'//div2//'.dat' open(4,file=gfile,form='unformatted',convert="little_endian") read(4)nxx,nyy read(4)lats4,lons4,zenith4,azimuth,solzen read(4)nxx,nyy,nll read(4)chem_wrf(:,:,:,k) close(4) end do o3d_wrf(:,:,:)=chem_wrf(:,:,:,1) print*,'Read WRFf 3d chem ',nxx,nyy,nll allocate(p3d_raqms(nxx,nyy,nk_raqms),t3d_raqms(nxx,nyy,nk_raqms),o3d_raqms(nxx,nyy,nk_raqms)) ! first read pressure and temperature !read RAQMS MET DATA do K=1,nmet gfile=trim(odir)//'raqms-on-wrf-chem.'//trim(raqms_3d_met_lab(k))//'.2006.08.24.06Z.W'//div2//'.dat' open(3,file=gfile,form='unformatted',convert="little_endian") read(3)nxx,nyy read(3)lats4,lons4,zenith4,azimuth,solzen read(3)nxx,nyy,nll if(raqms_3d_met_lab(k).eq.'Pressure')read(3)p3d_raqms if(raqms_3d_met_lab(k).eq.'Temperature')read(3)t3d_raqms close(3) end do print*,'Read raqms 3d data' !now read raqms ozone gfile=trim(odir)//'raqms-on-wrf-chem.o3vmr.2006.08.24.06Z.W'//div2//'.dat' open(3,file=gfile,form='unformatted',convert="little_endian") read(3)nxx,nyy read(3)lats4,lons4,zenith4,azimuth,solzen read(3)nxx,nyy,nll read(3)o3d_raqms close(3) print*,'Read raqms o3d' ! output file OPEN(23,file=trim(output_file),status='unknown') write(23,*)nxx,nyy write(23,*)nChannels i_prof=0 !MAIN LOOP!!!!!!!!!! do 81 i_loop=1,nxx do 82 j_loop=1,nyy if(t3d_wrf(i_loop,j_loop,1).lt.100.)go to 82 i_prof=i_prof+1 t100(:)=-99999. o3100(:)=-99999. q100(:)=0. cld100(:)=0. cldice100(:)=0. ch100(:,:)=-99999. !combine wrf and raqms atmospheres into one 100 layer profile !first do 100 layer atm from raqms, T and o3 ! raqms comes in with k=1 at top of atmosphere as needed by crtm !fill top twenty layers with raqms data - any layers above the ! top of raqms will be set to top value of raqms t100(1:20)=t3d_raqms(i_loop,j_loop,1) o3100(1:20)=o3d_raqms(i_loop,j_loop,1) do k=1,nk100 prs=play(k) do kr=1,nk_raqms-1 if(p3d_raqms(i_loop,j_loop,kr).le.prs.and.p3d_raqms(i_loop,j_loop,kr+1).ge.prs)then pm=log(prs) pu=log(p3d_raqms(i_loop,j_loop,kr)) pd=log(p3d_raqms(i_loop,j_loop,kr+1)) ppp=pu-pd ttt=t3d_raqms(i_loop,j_loop,kr)-t3d_raqms(i_loop,j_loop,kr+1) xx=ttt/ppp t100(k)=t3d_raqms(i_loop,j_loop,kr+1)-xx*pd+xx*pm ooo=o3d_raqms(i_loop,j_loop,kr)-o3d_raqms(i_loop,j_loop,kr+1) xx=ooo/ppp o3100(k)=o3d_raqms(i_loop,j_loop,kr+1)-xx*pd+xx*pm exit ! else if(prs.gt.p3d_raqms(i_loop,j_loop,nk_raqms))then t100(k)=t3d_raqms(i_loop,j_loop,nk_raqms) o3100(k)=o3d_raqms(i_loop,j_loop,nk_raqms) exit endif end do end do !now do 100 layer atm from WRF/Chem, replacing RAQMS where WRF available !NOTE wrf comes in with k=1 at earth sfc !!!!!! !k=nk_wrf is TOA !p3d_wrf(kr+1) .lt. p3d_wrf(kr) ! do met first do k=1,nk100 prs=play(k) do kr=1,nk_wrf-1 ! do nothing so raqms values remain if(prs.lt.p3d_wrf(i_loop,j_loop,nk_wrf))then if(p3d_wrf(i_loop,j_loop,kr).ge.prs.and.p3d_wrf(i_loop,j_loop,kr+1).le.prs)then pm=log(prs) pu=log(p3d_wrf(i_loop,j_loop,kr+1)) pd=log(p3d_wrf(i_loop,j_loop,kr)) ppp=pu-pd ttt=t3d_wrf(i_loop,j_loop,kr+1)-t3d_wrf(i_loop,j_loop,kr) xx=ttt/ppp t100(k)=t3d_wrf(i_loop,j_loop,kr)-xx*pd+xx*pm ooo=o3d_wrf(i_loop,j_loop,kr+1)-o3d_wrf(i_loop,j_loop,kr) xx=ooo/ppp o3100(k)=o3d_wrf(i_loop,j_loop,kr)-xx*pd+xx*pm qqq=q3d_wrf(i_loop,j_loop,kr+1)-q3d_wrf(i_loop,j_loop,kr) xx=qqq/ppp q100(k)=q3d_wrf(i_loop,j_loop,kr)-xx*pd+xx*pm ooo=cld3d_wrf(i_loop,j_loop,kr+1)-cld3d_wrf(i_loop,j_loop,kr) xx=ooo/ppp cld100(k)=cld3d_wrf(i_loop,j_loop,kr)-xx*pd+xx*pm ooo=cldice3d_wrf(i_loop,j_loop,kr+1)-cldice3d_wrf(i_loop,j_loop,kr) xx=ooo/ppp cldice100(k)=cldice3d_wrf(i_loop,j_loop,kr)-xx*pd+xx*pm exit else if(prs.gt.p3d_wrf(i_loop,j_loop,1))then t100(k)=t3d_wrf(i_loop,j_loop,1) o3100(k)=o3d_wrf(i_loop,j_loop,1) q100(k)=q3d_wrf(i_loop,j_loop,1) cld100(k)=cld3d_wrf(i_loop,j_loop,1) cldice100(k)=cldice3d_wrf(i_loop,j_loop,1) exit endif end do end do !Now do WRF/Chem aerosols ! do ichem=1,nchem !k=1 is ozone (done in the above loop), loop 2,nchem do ichem=2,nchem ch100(:,ichem)=0. do k=1,nk100 prs=play(k) do kr=1,nk_wrf-1 ! do nothing so values is zero; if(prs.lt.p3d_wrf(i_loop,j_loop,nk_wrf))then if(p3d_wrf(i_loop,j_loop,kr).ge.prs.and.p3d_wrf(i_loop,j_loop,kr+1).le.prs)then pm=log(prs) pu=log(p3d_wrf(i_loop,j_loop,kr+1)) pd=log(p3d_wrf(i_loop,j_loop,kr)) ttt=chem_wrf(i_loop,j_loop,kr+1,ichem)-chem_wrf(i_loop,j_loop,kr,ichem) ppp=pu-pd xx=ttt/ppp ch100(k,ichem)=chem_wrf(i_loop,j_loop,kr,ichem)-xx*pd+xx*pm exit else if(prs.gt.p3d_wrf(i_loop,j_loop,1))then ch100(k,ichem)=chem_wrf(i_loop,j_loop,1,ichem) exit endif end do end do end do if(i_prof.eq.2)then 1121 format(1x,i3,2x,3(f7.2,2x),4(e10.3,2x)) print*,'Final Column met' print*,'top press ',prlvl(1) do k=1,nk100 write(6,1121)k,play(k),prlvl(k+1),t100(k),q100(k),cld100(k),o3100(k),ch100(k,4) end do stop endif 82 continue 81 continue 212 continue ! here T, 03, q, cldwater, cldice and 14 aerosols should be on the 100 pressure layers in play above ! t100, 03100, q100, cld100, cldice100 and ch100(nk100,2:nchem) stop end