! module ncep_rean_module ! ! Read T,U,V,Z lower boundary (10 mb) data from ncep reanalysis data file. ! The file is specified with namelist input parameter "ncep_reanalysis" ! (see input.F) ! See hao:/home/tgcm/ncep_reanalysis for procedure to acquire the data ! (NCEP/NCAR dataset ds090.0), and make the 10 mb ncep_reanalysis netcdf ! file. ! use params_module,only: nlon,nlonp4,glon1,nlat,dlat,dlon,nlevp1 use cons_module,only: grav,gask,atm_amu use nchist_module,only: handle_ncerr,nc_open use addfld_module,only: addfld implicit none integer :: ncid,ncep_nlev,ncep_ntime #include ! ! Year, day of year, and hour from ncep file: integer,allocatable,dimension(:) :: ncep_year,ncep_day,ncep_hour ! ! Output of this module is T,Z at 10 mb interpolated to current time. real,dimension(nlonp4,nlat) :: z_ncep, t_ncep, u_ncep, v_ncep character(len=16),private :: units_z,units_u,units_v ! ! Use daily ncep reanalysis: logical, parameter :: ncep_rean_daily = .true. contains !----------------------------------------------------------------------- subroutine ncep_rean_bndry(ncfile,istep,mtime,iyear,iday, | isecs) ! ! Called from advance. Return global z_ncep and t_ncep at current time. ! implicit none ! ! Args: character(len=*),intent(in) :: ncfile integer,intent(in) :: | istep, ! model time step (1->nstep) | iyear, ! current model year (4-digit) | iday, ! current model calendar day (1-367) | isecs, ! current model time in seconds | mtime(4) ! model time (day,hr,min,sec), for info only ! ! Local: integer :: i real,dimension(nlon,nlat) :: tncep0,tncep,tncep1, | zncep0,zncep,zncep1, | uncep0,uncep,uncep1, | vncep0,vncep,vncep1 real :: hour ! hour = float(isecs)/3600. ! ! If first step, open and read coord vars of ncep data file: if (istep==1) then write(6,"('ncep_rean_bndry: ncfile=',a)") ncfile call get_ncep_file(ncfile) endif ! ! Get daily ncep reanalysis: if (ncep_rean_daily) then call get_ncep_daily(iyear,iday,hour) ! ! Get ncep data at current model time: else call get_ncep_time(iyear,iday,hour) endif end subroutine ncep_rean_bndry !----------------------------------------------------------------------- subroutine get_ncep_file(ncfile) implicit none ! ! Args: character(len=*),intent(in) :: ncfile ! ! Local: integer :: istat,ndims,idunlim,ngatts,nvars integer :: nlat_rd,nlon_rd integer :: id_lat,id_lon,id_time,id_lev integer :: idv_lev,idv_year,idv_day,idv_hour integer :: istart1(1),icount1(1) character(len=80) :: dskfile real,allocatable :: ncep_levs(:) ! dskfile = ' ' call getfile(ncfile,dskfile) call nc_open(ncid,ncfile,'OLD','READ') istat = nf_inq(ncid,ndims,nvars,ngatts,idunlim) ! write(6,"('get_ncep_file: ndims=',i3,' nvars=',i3, ! | ' ngatts=',i3,' idunlim=',i3)") ndims,nvars,ngatts,idunlim ! ! Get and verify latitude dimension. ! (read nlat_rd, compare with nlat, which is in params.h): istat = nf_inq_dimid(ncid,'lat',id_lat) istat = nf_inq_dimlen(ncid,id_lat,nlat_rd) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'get_ncep_file: error getting latitude dimension') if (nlat_rd /= nlat) then write(6,"(/,'>>> get_ncep_file: bad nlat_rd=',i3, | ' -- should be nlat=',i3)") nlat_rd,nlat call shutdown('get_ncep_file') endif ! ! Get and verify longitude dimension. ! (read nlon_rd, compare with nlon, which is in params.h): istat = nf_inq_dimid(ncid,'lon',id_lon) istat = nf_inq_dimlen(ncid,id_lon,nlon_rd) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'get_ncep_file: error getting longitude dimension') if (nlon_rd /= nlon) then write(6,"(/,'>>> get_ncep_file: bad nlon_rd=',i3, | ' -- should be nlon=',i3)") nlon_rd,nlon call shutdown('get_ncep_file') endif ! ! Get levels dimension (should be 1, i.e., only 10 mb level): istat = nf_inq_dimid(ncid,'lev',id_lev) istat = nf_inq_dimlen(ncid,id_lev,ncep_nlev) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'get_ncep_file: error getting levels dimension') ! ! Verify that first level on file = 10 mb. istat = nf_inq_varid(ncid,"lev",idv_lev) ! day var id allocate(ncep_levs(ncep_nlev),stat=istat) istart1(1) = 1 icount1(1) = ncep_nlev istat = nf_get_vara_double(ncid,idv_lev,istart1,icount1, | ncep_levs) write(6,"('get_ncep_file: ncep_nlev=',i3,' ncep_levs(1)=',f8.2)") | ncep_nlev,ncep_levs(1) if (ncep_levs(1) /= 10.) then write(6,"('>>> ncep_levs(1)=',f8.2,' -- expected 10. mb')") | ncep_levs(1) call shutdown('get_ncep_file') endif deallocate(ncep_levs) ! ! Get number of times (value of unlimited dimension) on the file: istat = nf_inq_unlimdim(ncid,id_time) ! id of unlimited record var if (istat /= NF_NOERR) call handle_ncerr(istat, | 'get_ncep_file: error getting unlimited dimension id') istat = nf_inq_dimlen(ncid,id_time,ncep_ntime) ! length of time var if (istat /= NF_NOERR) then call handle_ncerr(istat, | 'Error getting length of day record variable') call shutdown('get_ncep_file') else ! write(6,"('get_ncep_file: ncep_ntime=',i3)") ncep_ntime endif ! ! Read vars year,day_of_year,hour: ! integer,allocatable,dimension(:) :: ncep_year,ncep_day,ncep_hour ! allocate(ncep_year(ncep_ntime),stat=istat) allocate(ncep_day (ncep_ntime),stat=istat) allocate(ncep_hour(ncep_ntime),stat=istat) istart1(1) = 1 icount1(1) = ncep_ntime ! ! Year: istat = nf_inq_varid(ncid,"year",idv_year) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'get_ncep_file: error getting id of var year') istat = nf_get_vara_int(ncid,idv_year,istart1,icount1, | ncep_year) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'get_ncep_file: error reading var year') ! write(6,"('get_ncep_file: year=',/,(10i6))") ncep_year ! ! Day (of year): istat = nf_inq_varid(ncid,"day",idv_day) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'get_ncep_file: error getting id of var day (day of year)') istat = nf_get_vara_int(ncid,idv_day,istart1,icount1, | ncep_day) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'get_ncep_file: error reading var day (day of year)') ! write(6,"('get_ncep_file: day of year=',/,(10i6))") ncep_day ! ! Hour: istat = nf_inq_varid(ncid,"hour",idv_hour) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'get_ncep_file: error getting id of var hour') istat = nf_get_vara_int(ncid,idv_hour,istart1,icount1, | ncep_hour) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'get_ncep_file: error reading var hour') ! write(6,"('get_ncep_file: hour=',/,(10i6))") ncep_hour end subroutine get_ncep_file !----------------------------------------------------------------------- subroutine get_ncep_time(iyear,iday,hour) implicit none ! ! Args: integer,intent(in) :: iyear,iday real,intent(in) :: hour ! ! Local: integer :: i,itncep0,itncep1,ncephr1 real :: hour0,hour1 ! ! write(6,"('get_ncep_time: searching for iyear=',i4,' iday=',i3, ! | ' hour=',f8.3)") iyear,iday,hour itncep0 = 0 itncep1 = 0 ! ! Check first time (if exact match, interp will not be necessary): if (ncep_year(1)==iyear.and.ncep_day(1)==iday.and. | float(ncep_hour(1))==hour) then itncep0 = 1 endif ! ! Search interior times: if (itncep0==0) then searchloop: do i=1,ncep_ntime-1 ! write(6,"(/,'get_ncep_time: searching for iyear=',i4,' iday=', ! | i4,' hour=',f9.3)") iyear,iday,hour ! write(6,"('i =',i5,': found ncep year=',i4,' day=',i4, ! | ' hour=',f9.3)") i,ncep_year(i),ncep_day(i), ! | float(ncep_hour(i)) ! write(6,"('i+1=',i5,': found ncep year=',i4,' day=',i4, ! | ' hour=',f9.3)") i+1,ncep_year(i+1),ncep_day(i+1), ! | float(ncep_hour(i+1)) ! ! Check for exact match at i: if (ncep_year(i)==iyear.and.ncep_day(i)==iday.and. | float(ncep_hour(i))==hour) then itncep0 = i exit searchloop endif ! ! Check for exact match at i+1: if (ncep_year(i+1)==iyear.and.ncep_day(i+1)==iday.and. | float(ncep_hour(i+1))==hour) then itncep0 = i+1 exit searchloop endif ! ! Check for bracket between i and i+1 (interp will be necessary): ! Note if ncep_hour(i) is 18, then ncep_hour(i+1) is 0, but need ! to test between 18 and 24 (not 18 and 0), see ncephr1: if (itncep0==0) then ncephr1 = ncep_hour(i+1) if (ncep_hour(i)==18) ncephr1 = 24 if (ncep_year(i)==iyear.and.ncep_day(i)==iday.and. | hour > float(ncep_hour(i)).and. | hour < float(ncephr1)) then itncep0 = i itncep1 = i+1 exit searchloop endif endif enddo searchloop endif ! ! Check for exact match at last time: if (itncep0==0) then if (ncep_year(ncep_ntime)==iyear.and. | ncep_day(ncep_ntime)==iday.and. | hour==ncep_hour(ncep_ntime)) then itncep0 = ncep_ntime endif endif ! ! Time not found: if (itncep0==0) then write(6,"('>>> get_ncep_time: could not find iyear=',i4, | ' iday=',i3,' hour=',f8.3)") iyear,iday,hour call shutdown('get_ncep_time') endif ! ! write(6,"(/,'get_ncep_time: found iyear,iday,hour=',i5,i3,f8.3, ! | ' between itncep0,1=',2i4)") iyear,iday,hour,itncep0,itncep1 ! write(6,"('ncep_year,day,hour at itncep0=',3i4)") ! | ncep_year(itncep0),ncep_day(itncep0),ncep_hour(itncep0) ! if (itncep1 > 0) ! | write(6,"('ncep_year,day,hour at itncep1=',3i4)") ! | ncep_year(itncep1),ncep_day(itncep1),ncep_hour(itncep1) ! ! Read ncep data and do interpolation if necessary: ! (4th arg, ncep_hour, is ignored if itncep1==0) if (itncep1==0) then ! exact time match at itncep0 hour0 = float(ncep_hour(itncep0)) hour1 = float(ncep_hour(itncep1)) call get_ncep_data(itncep0,itncep1,hour0,hour1,hour) else ! will interpolate between it0,it1 ncephr1 = ncep_hour(itncep1) if (ncephr1==0) ncephr1 = 24 hour0 = float(ncep_hour(itncep0)) hour1 = float(ncephr1) call get_ncep_data(itncep0,itncep1,hour0,hour1,hour) endif end subroutine get_ncep_time !----------------------------------------------------------------------- subroutine get_ncep_daily(iyear,iday,hour) implicit none ! ! Args: integer,intent(in) :: iyear,iday real,intent(in) :: hour ! ! Local: integer :: i,itncep0,itncep1,ncephr1 real :: hour0,hour1,hour_int ! itncep0 = 0 itncep1 = 0 ! ! Search interior times for model day: if (itncep0==0) then searchloop: do i=1,ncep_ntime if (ncep_year(i)==iyear.and.ncep_day(i)==iday.and. | float(ncep_hour(i))==0.) then itncep0 = i exit searchloop endif enddo searchloop endif ! ! check if exact time (12UT) or need day before or day after this model day ! for interpolation (only for first and last day of NCEP don't interpolate ! ncep_ntime-3 since we get the itncep0 for 0 UT and the file last time is 18 UT if(hour==12.) then itncep1 = 0 elseif(hour > 12. .and. itncep0 < ncep_ntime-3) then itncep1 = itncep0 + 4 hour0 = 0. hour1 = 24. hour_int = hour-12. if(hour_int < 0.or.hour_int >12.) then write(6,*) 'error in subroutine get_ncep_daily: ', | 'hour_int= ',hour_int stop endif elseif(hour <= 12..and. itncep0 > 4) then itncep1 = itncep0 itncep0 = itncep0 - 4 hour0 = 0. hour1 = 24. hour_int = hour+12. if(hour_int > 24.or.hour_int<12.) then write(6,*) 'error in subroutine get_ncep_daily: ', | 'hour_int= ',hour_int stop endif endif ! ! Time not found: if (itncep0==0) then write(6,"('>>> get_ncep_time: could not find iyear=',i4, | ' iday=',i3)") iyear,iday call shutdown('get_ncep_time') endif ! ! write(6,"(/,'get_ncep_time: found iyear,iday,hour=',i5,i3,f8.3, ! | ' between itncep0,1=',2i4)") iyear,iday,hour,itncep0,itncep1 ! write(6,"('ncep_year,day,hour at itncep0=',3i4)") ! | ncep_year(itncep0),ncep_day(itncep0),ncep_hour(itncep0) ! if (itncep1 > 0) ! | write(6,"('ncep_year,day,hour at itncep1=',3i4)") ! | ncep_year(itncep1),ncep_day(itncep1),ncep_hour(itncep1) ! ! Read ncep data and do daily average call get_ncep_data_daily(itncep0,itncep1, | hour0,hour1,hour_int) end subroutine get_ncep_daily !----------------------------------------------------------------------- subroutine get_ncep_data(it0,it1,hour0,hour1,hour) ! ! Read ncep data at time it0. If it1 > 0, then read ncep data ! at both it0 and it1, and interpolate to hour. ! implicit none ! ! Args: integer,intent(in) :: it0,it1 real,intent(in) :: hour0,hour1,hour ! ! Local: integer :: i,j,k real,dimension(nlon,nlat) :: tncep0,tncep1,zncep0,zncep1, | uncep0,uncep1,vncep0,vncep1 real,dimension(nlonp4,nlevp1) :: tncepik,zncepik,uncepik,vncepik ! ! External (util.F): real,external :: time_interp ! ! Read data at it0: call rd_ncep_data(it0,tncep0,zncep0,uncep0,vncep0) ! write(6,"('get_ncep_data: it0,1=',2i4,' hour=',f7.3, ! | ' tncep0 min,max=',2e12.4,' zncep0 min,max=',2e12.4)") ! | it0,it1,hour,minval(tncep0),maxval(tncep0), ! | minval(zncep0),maxval(zncep0) ! ! If it1 == 0, then no need to interpolate: ! (exact time match was found at it0): if (it1==0) then t_ncep(3:nlon+2,:) = tncep0(:,:) z_ncep(3:nlon+2,:) = zncep0(:,:) u_ncep(3:nlon+2,:) = uncep0(:,:) v_ncep(3:nlon+2,:) = vncep0(:,:) ! ! If it1 > 0, then interpolate to needed hour: elseif (it1 > 0) then call rd_ncep_data(it1,tncep1,zncep1,uncep1,vncep1) ! ! Double check interpolates: if (hour0 > hour1) then write(6,"('>>> get_ncep_data: hour0 must be < hour1: hour0=', | f8.3,' hour1=',f8.3)") call shutdown('get_ncep_data') elseif (hour < hour0.or.hour > hour1) then write(6,"('>>> get_ncep_data: hour must be between hour0,1:', | ' hour0,1=',2f8.3,' hour=',f8.3)") hour0,hour1,hour call shutdown('get_ncep_data') endif ! ! Interpolate to hour: ! Util.F: real function time_interp(f0,f1,time0,time1,time) ! do j=1,nlat do i=3,nlon+2 t_ncep(i,j) = time_interp(tncep0(i-2,j),tncep1(i-2,j), | hour0,hour1,hour) z_ncep(i,j) = time_interp(zncep0(i-2,j),zncep1(i-2,j), | hour0,hour1,hour) u_ncep(i,j) = time_interp(uncep0(i-2,j),uncep1(i-2,j), | hour0,hour1,hour) v_ncep(i,j) = time_interp(vncep0(i-2,j),vncep1(i-2,j), | hour0,hour1,hour) enddo enddo endif ! it1 > 0 ! ! 5/24/07: Units of z,u,v were read by rd_ncep_data. Use these to convert ! from m to cm if necessary. See hao:/home/tgcm/ncep_reanalysis/ ! regrid.ncl. ! ! write(6,"('get_ncep_data: units_z=',a,' units_u=',a,' units_v=', ! | a)") trim(units_z),trim(units_u),trim(units_v) if (units_z(1:1)=='m') then z_ncep = z_ncep*100. ! write(6,"('Multiply z_ncep by 100 for conversion from', ! | ' m to cm')") endif if (units_u(1:3)=='m/s') then u_ncep = u_ncep*100. ! write(6,"('Multiply u_ncep by 100 for conversion from', ! | ' m/s to cm/s')") endif if (units_v(1:3)=='m/s') then v_ncep = v_ncep*100. ! write(6,"('Multiply v_ncep by 100 for conversion from', ! | ' m/s to cm/s')") endif ! ! Periodic points: ! lons 1,2 <- nlon+1,nlon+2, and nlonp+3,nlon+4 <- 3,4 ! t_ncep(1:2,:) = t_ncep(nlon+1:nlon+2,:) t_ncep(nlon+3:nlon+4,:) = t_ncep(3:4,:) z_ncep(1:2,:) = z_ncep(nlon+1:nlon+2,:) z_ncep(nlon+3:nlon+4,:) = z_ncep(3:4,:) u_ncep(1:2,:) = u_ncep(nlon+1:nlon+2,:) u_ncep(nlon+3:nlon+4,:) = u_ncep(3:4,:) v_ncep(1:2,:) = v_ncep(nlon+1:nlon+2,:) v_ncep(nlon+3:nlon+4,:) = v_ncep(3:4,:) ! ! write(6,"('get_ncep_data: it0,1=',2i4,' hour0=',f7.3,' hour=', ! | f7.3,' hour1=',f7.3,' tncep0 mnmx=',2e12.4,' tncep1 mnmx=', ! | 2e12.4)") it0,it1,hour0,hour,hour1, ! | minval(tncep0),maxval(tncep0),minval(tncep1),maxval(tncep1) ! write(6,"('get_ncep_data: it0,1=',2i4,' hour0=',f7.3,' hour=', ! | f7.3,' hour1=',f7.3,' t_ncep mnmx=',2e12.4)") ! | it0,it1,hour0,hour,hour1,minval(t_ncep),maxval(t_ncep) ! ! write(6,"('get_ncep_data: it0,1=',2i4,' hour0=',f7.3,' hour=', ! | f7.3,' hour1=',f7.3,' zncep0 mnmx=',2e12.4,' zncep1 mnmx=', ! | 2e12.4)") it0,it1,hour0,hour,hour1, ! | minval(zncep0),maxval(zncep0),minval(zncep1),maxval(zncep1) ! write(6,"('get_ncep_data: it0,1=',2i4,' hour0=',f7.3,' hour=', ! | f7.3,' hour1=',f7.3,' z_ncep mnmx=',2e12.4)") ! | it0,it1,hour0,hour,hour1,minval(z_ncep),maxval(z_ncep) ! ! call addfld('NCEP_T','NCEP reanalysis: TN at 10 mb','K', ! | t_ncep,'lon',1,nlonp4,'lat',1,nlat,0) ! call addfld('NCEP_Z','NCEP reanalysis: Z at 10 mb','cm', ! | z_ncep,'lon',1,nlonp4,'lat',1,nlat,0) ! ! real,dimension(nlonp4,nlevp1) :: tncepik,zncepik ! Save on sec hist redundantly in vertical for tgcmproc_f90: ! do j=1,nlat do i=3,nlon+2 tncepik(i,:) = t_ncep(i,j) zncepik(i,:) = z_ncep(i,j) uncepik(i,:) = u_ncep(i,j) vncepik(i,:) = v_ncep(i,j) enddo do k=1,nlevp1 tncepik(1:2,k) = tncepik(nlon+1:nlon+2,j) tncepik(nlon+3:nlon+4,k) = tncepik(3:4,j) zncepik(1:2,k) = zncepik(nlon+1:nlon+2,j) zncepik(nlon+3:nlon+4,k) = zncepik(3:4,j) uncepik(1:2,k) = uncepik(nlon+1:nlon+2,j) uncepik(nlon+3:nlon+4,k) = uncepik(3:4,j) vncepik(1:2,k) = vncepik(nlon+1:nlon+2,j) vncepik(nlon+3:nlon+4,k) = vncepik(3:4,j) enddo ! call addfld('NCEP_T','NCEP reanalysis: TN at 10 mb','K', | tncepik,'lon',1,nlonp4,'lev',1,nlevp1,j) call addfld('NCEP_Z','NCEP reanalysis: Z at 10 mb','cm', | zncepik,'lon',1,nlonp4,'lev',1,nlevp1,j) call addfld('NCEP_U','NCEP reanalysis: U at 10 mb','cm', | uncepik,'lon',1,nlonp4,'lev',1,nlevp1,j) call addfld('NCEP_V','NCEP reanalysis: V at 10 mb','cm', | vncepik,'lon',1,nlonp4,'lev',1,nlevp1,j) enddo ! end subroutine get_ncep_data !----------------------------------------------------------------------- subroutine get_ncep_data_daily(it0,it1,hour0,hour1,hour) ! ! Read ncep data at time it0, it0+1, it0+2, it0+3; then do a daily average ! implicit none ! ! Args: integer,intent(in) :: it0,it1 real,intent(in) :: hour0,hour1,hour ! ! Local: integer :: i,j,k real,dimension(nlon,nlat) :: tncep0,tncep1, tncep2,tncep3, | zncep0,zncep1,zncep2,zncep3, | uncep0,uncep1,uncep2,uncep3, | vncep0,vncep1,vncep2,vncep3, | tncep_d0,zncep_d0,uncep_d0,vncep_d0, | tncep_d1,zncep_d1,uncep_d1,vncep_d1 real,dimension(nlonp4,nlevp1) :: tncepik,zncepik,uncepik, | vncepik ! ! External (util.F): real,external :: time_interp ! ! ! Read data at it0: call rd_ncep_data(it0,tncep0,zncep0,uncep0,vncep0) call rd_ncep_data(it0+1,tncep1,zncep1,uncep1,vncep1) call rd_ncep_data(it0+2,tncep2,zncep2,uncep2,vncep2) call rd_ncep_data(it0+3,tncep3,zncep3,uncep3,vncep3) ! !If it1 == 0, then no need to interpolate: ! (exact time match was found at it0): if (it1==0) then do j=1,nlat do i=3,nlon+2 t_ncep(i,j) = 0.25*(tncep0(i-2,j)+tncep1(i-2,j)+ | tncep2(i-2,j)+tncep3(i-2,j)) z_ncep(i,j) = 0.25*(zncep0(i-2,j)+zncep1(i-2,j)+ | zncep2(i-2,j)+zncep3(i-2,j)) u_ncep(i,j) = 0.25*(uncep0(i-2,j)+uncep1(i-2,j)+ | uncep2(i-2,j)+uncep3(i-2,j)) v_ncep(i,j) = 0.25*(vncep0(i-2,j)+vncep1(i-2,j)+ | vncep2(i-2,j)+vncep3(i-2,j)) enddo enddo ! ! If it1 > 0, then interpolate to needed hour: elseif (it1 > 0) then ! calculate avarge of it0 do j=1,nlat do i=1,nlon tncep_d0(i,j) = 0.25*(tncep0(i,j)+tncep1(i,j)+ | tncep2(i,j)+tncep3(i,j)) zncep_d0(i,j) = 0.25*(zncep0(i,j)+zncep1(i,j)+ | zncep2(i,j)+zncep3(i,j)) uncep_d0(i,j) = 0.25*(uncep0(i,j)+uncep1(i,j)+ | uncep2(i,j)+uncep3(i,j)) vncep_d0(i,j) = 0.25*(vncep0(i,j)+vncep1(i,j)+ | vncep2(i,j)+vncep3(i,j)) enddo enddo call rd_ncep_data(it1 ,tncep0,zncep0,uncep0,vncep0) call rd_ncep_data(it1+1,tncep1,zncep1,uncep1,vncep1) call rd_ncep_data(it1+2,tncep2,zncep2,uncep2,vncep2) call rd_ncep_data(it1+3,tncep3,zncep3,uncep3,vncep3) ! calculate avarge of it1 do j=1,nlat do i=1,nlon tncep_d1(i,j) = 0.25*(tncep0(i,j)+tncep1(i,j)+ | tncep2(i,j)+tncep3(i,j)) zncep_d1(i,j) = 0.25*(zncep0(i,j)+zncep1(i,j)+ | zncep2(i,j)+zncep3(i,j)) uncep_d1(i,j) = 0.25*(uncep0(i,j)+uncep1(i,j)+ | uncep2(i,j)+uncep3(i,j)) vncep_d1(i,j) = 0.25*(vncep0(i,j)+vncep1(i,j)+ | vncep2(i,j)+vncep3(i,j)) enddo enddo ! ! Interpolate to hour: ! Util.F: real function time_interp(f0,f1,time0,time1,time) do j=1,nlat do i=3,nlon+2 t_ncep(i,j) = time_interp(tncep_d0(i-2,j), | tncep_d1(i-2,j),hour0,hour1,hour) z_ncep(i,j) = time_interp(zncep_d0(i-2,j), | zncep_d1(i-2,j),hour0,hour1,hour) u_ncep(i,j) = time_interp(uncep_d0(i-2,j), | uncep_d1(i-2,j),hour0,hour1,hour) v_ncep(i,j) = time_interp(vncep_d0(i-2,j), | vncep_d1(i-2,j),hour0,hour1,hour) enddo enddo endif ! ! Periodic points: ! lons 1,2 <- nlon+1,nlon+2, and nlonp+3,nlon+4 <- 3,4 ! t_ncep(1:2,:) = t_ncep(nlon+1:nlon+2,:) t_ncep(nlon+3:nlon+4,:) = t_ncep(3:4,:) z_ncep(1:2,:) = z_ncep(nlon+1:nlon+2,:) z_ncep(nlon+3:nlon+4,:) = z_ncep(3:4,:) u_ncep(1:2,:) = u_ncep(nlon+1:nlon+2,:) u_ncep(nlon+3:nlon+4,:) = u_ncep(3:4,:) v_ncep(1:2,:) = v_ncep(nlon+1:nlon+2,:) v_ncep(nlon+3:nlon+4,:) = v_ncep(3:4,:) ! ! am 11/09: convert units like it was done for 6 hourly ncep_rean. data ! 5/24/07: Units of z,u,v were read by rd_ncep_data. Use these to convert ! from m to cm if necessary. See hao:/home/tgcm/ncep_reanalysis/ ! regrid.ncl. ! ! write(6,"('get_ncep_data: units_z=',a,' units_u=',a,' units_v=', ! | a)") trim(units_z),trim(units_u),trim(units_v) if (units_z(1:1)=='m') then z_ncep = z_ncep*100. ! write(6,"('Multiply z_ncep by 100 for conversion from', ! | ' m to cm')") endif if (units_u(1:3)=='m/s') then u_ncep = u_ncep*100. ! write(6,"('Multiply u_ncep by 100 for conversion from', ! | ' m/s to cm/s')") endif if (units_v(1:3)=='m/s') then v_ncep = v_ncep*100. ! write(6,"('Multiply v_ncep by 100 for conversion from', ! | ' m/s to cm/s')") endif ! ! ! real,dimension(nlonp4,nlevp1) :: tncepik,zncepik ! Save on sec hist redundantly in vertical for tgcmproc_f90: ! do j=1,nlat do i=3,nlon+2 tncepik(i,:) = t_ncep(i,j) zncepik(i,:) = z_ncep(i,j) uncepik(i,:) = u_ncep(i,j) vncepik(i,:) = v_ncep(i,j) enddo do k=1,nlevp1 tncepik(1:2,k) = tncepik(nlon+1:nlon+2,j) tncepik(nlon+3:nlon+4,k) = tncepik(3:4,j) zncepik(1:2,k) = zncepik(nlon+1:nlon+2,j) zncepik(nlon+3:nlon+4,k) = zncepik(3:4,j) uncepik(1:2,k) = uncepik(nlon+1:nlon+2,j) uncepik(nlon+3:nlon+4,k) = uncepik(3:4,j) vncepik(1:2,k) = vncepik(nlon+1:nlon+2,j) vncepik(nlon+3:nlon+4,k) = vncepik(3:4,j) enddo ! call addfld('NCEP_T','NCEP: TN at 10 mb','K', | tncepik,'lon',1,nlonp4,'lev',1,nlevp1,j) call addfld('NCEP_Z','NCEP: Z at 10 mb','cm', | zncepik,'lon',1,nlonp4,'lev',1,nlevp1,j) call addfld('NCEP_U','NCEP: U at 10 mb','cm/s', | uncepik,'lon',1,nlonp4,'lev',1,nlevp1,j) call addfld('NCEP_V','NCEP: V at 10 mb','cm/s', | vncepik,'lon',1,nlonp4,'lev',1,nlevp1,j) enddo ! end subroutine get_ncep_data_daily !----------------------------------------------------------------------- subroutine rd_ncep_data(it,tncep,zncep,uncep,vncep) ! ! Read ncep T,Z,U,V at time index it. ! implicit none ! ! Args: integer,intent(in) :: it real,dimension(nlon,nlat),intent(out) :: tncep,zncep,uncep,vncep ! ! Local: integer :: istat,start_4d(4),count_4d(4),idv_t,idv_z,idv_u, | idv_v,lenunits character(len=120) :: char120 ! start_4d(1:3) = 1 start_4d(4) = it count_4d(1) = nlon count_4d(2) = nlat count_4d(3) = ncep_nlev count_4d(4) = 1 ! ! Read T: istat = nf_inq_varid(ncid,"TMP_2_ISBL_10",idv_t) if (istat /= NF_NOERR) | istat = nf_inq_varid(ncid,"TMP_2_ISBL",idv_t) if (istat /= NF_NOERR) then call handle_ncerr(istat, | 'Error getting var id of ncep TMP_2_ISBL_10 or TMP_2_ISBL') call shutdown('rd_ncep_data') endif istat = nf_get_vara_double(ncid,idv_t,start_4d,count_4d,tncep) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_get_vara_double', | ' for field tncep: it=',i4)") it call handle_ncerr(istat,char120) else ! write(6,"('rd_ncep_data: read tncep: it=',i4,' min,max=', ! | 2e12.4)") it,minval(tncep),maxval(tncep) endif ! ! Read Z: istat = nf_inq_varid(ncid,"HGT_2_ISBL_10",idv_z) if (istat /= NF_NOERR) | istat = nf_inq_varid(ncid,"HGT_2_ISBL",idv_z) if (istat /= NF_NOERR) then call handle_ncerr(istat, | 'Error getting var id of ncep HGT_2_ISBL_10 or HGT_2_ISBL') call shutdown('rd_ncep_data') endif istat = nf_get_vara_double(ncid,idv_z,start_4d,count_4d,zncep) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_get_vara_double', + ' for field zncep')") call handle_ncerr(istat,char120) else ! write(6,"('rd_ncep_data: read zncep: it=',i4,' min,max=', ! | 2e12.4)") it,minval(zncep),maxval(zncep) endif units_z = ' ' istat = nf_get_att_text(ncid,idv_z,"units",units_z) if (istat /= NF_NOERR) then call handle_ncerr(istat, | 'Error getting units attribute of ncep HGT_2_ISBL_10') call shutdown('rd_ncep_data') endif ! ! Read U: istat = nf_inq_varid(ncid,"U_GRD_2_ISBL_10",idv_u) if (istat /= NF_NOERR) | istat = nf_inq_varid(ncid,"U_GRD_2_ISBL",idv_u) if (istat /= NF_NOERR) then call handle_ncerr(istat, | 'Error getting var id of ncep U_GRD_2_ISBL_10 or U_GRD_ISBL') call shutdown('rd_ncep_data') endif istat = nf_get_vara_double(ncid,idv_u,start_4d,count_4d,uncep) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_get_vara_double', + ' for field uncep')") call handle_ncerr(istat,char120) else ! write(6,"('rd_ncep_data: read uncep: it=',i4,' min,max=', ! | 2e12.4)") it,minval(uncep),maxval(uncep) endif units_u = ' ' ! istat = nf_inq_attlen(ncid,idv_u,"units",lenunits) istat = nf_get_att_text(ncid,idv_u,"units",units_u) if (istat /= NF_NOERR) then call handle_ncerr(istat, | 'Error getting units attribute of ncep U_GRD_2_ISBL_10') call shutdown('rd_ncep_data') endif ! ! Read V: istat = nf_inq_varid(ncid,"V_GRD_2_ISBL_10",idv_v) if (istat /= NF_NOERR) | istat = nf_inq_varid(ncid,"V_GRD_2_ISBL",idv_v) if (istat /= NF_NOERR) then call handle_ncerr(istat, | 'Error getting var id of ncep V_GRD_2_ISBL_10 or V_GRD_2_ISBL') call shutdown('rd_ncep_data') endif istat = nf_get_vara_double(ncid,idv_v,start_4d,count_4d,vncep) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_get_vara_double', + ' for field vncep')") call handle_ncerr(istat,char120) else ! write(6,"('rd_ncep_data: read vncep: it=',i4,' min,max=', ! | 2e12.4)") it,minval(vncep),maxval(vncep) endif units_v = ' ' istat = nf_get_att_text(ncid,idv_v,"units",units_v) if (istat /= NF_NOERR) then call handle_ncerr(istat, | 'Error getting units attribute of ncep V_GRD_2_ISBL_10') call shutdown('rd_ncep_data') endif end subroutine rd_ncep_data !----------------------------------------------------------------------- end module ncep_rean_module