program plevel ! ! Compilation on GFDL Altix platform: ! ifort -r8 -I/usr/local/include plevel.f90 -L/usr/local/lib -lnetcdf ! use netcdf implicit none ! number of pressure levels in the output file integer, parameter :: nplev = 17 ! pressure levels in Pascals real, dimension(nplev) :: plev = (/ 1000.e2, 925.e2, 850.e2, 700.e2, 600.e2, 500.e2, & 400.e2, 300.e2, 250.e2, 200.e2, 150.e2, 100.e2, & 70.e2, 50.e2, 30.e2, 20.e2, 10.e2 /) ! field name character(len=NF90_MAX_NAME) :: fieldname = "ta" ! input file name character(len=NF90_MAX_NAME) :: file_in = "ta_gfdl_2070120100.nc" ! output file name character(len=NF90_MAX_NAME) :: file_out = "ta_pres_2070120100.nc" logical :: xtop = .true. ! extrapolate above top level logical :: xbot = .false. ! extrapolate below bottom level !---------- probably no changes needed below here ----------- integer, parameter :: nlon = 217, nlat = 131, nlev = 24 real :: lon(nlon), lat(nlat) real :: a_bnds(2,nlev), b_bnds(2,nlev) real :: ps(nlon,nlat), data_in(nlon,nlat,nlev) real :: p_bnds(2), lnp(nlev) real :: data_out(nlon,nlat,nplev), lnpout(nplev) real :: time integer :: i, j, k, kp, km, n integer :: istat, ncid_in, ncid_out, recdim, numtime, varid_in(7), dimid_out(4), varid_out(5) integer :: ntop, nbot real :: fact character(len=NF90_MAX_NAME) :: name_recdim character(len=8) :: pct !----------------------------------- ! open the input file istat = NF90_OPEN (trim(file_in), NF90_NOWRITE, ncid_in) if (istat /= NF90_NOERR) call error_handler (ncode=istat) ! get record dimension information istat = NF90_INQUIRE (ncid_in, unlimiteddimid=recdim) istat = NF90_INQUIRE_DIMENSION (ncid_in, recdim, name=name_recdim, len=numtime) ! get model grid istat = NF90_INQ_VARID (ncid_in, 'lon', varid_in(1)) istat = NF90_INQ_VARID (ncid_in, 'lat', varid_in(2)) ! get model vertical grid istat = NF90_INQ_VARID (ncid_in, 'a_bnds', varid_in(3)) istat = NF90_INQ_VARID (ncid_in, 'b_bnds', varid_in(4)) ! record dimension istat = NF90_INQ_VARID (ncid_in, trim(name_recdim), varid_in(5)) ! surface pressure istat = NF90_INQ_VARID (ncid_in, 'ps', varid_in(6)) ! 3d field istat = NF90_INQ_VARID (ncid_in, trim(fieldname), varid_in(7)) ! read static axis information istat = NF90_GET_VAR (ncid_in, varid_in(1), lon) if (istat /= NF90_NOERR) call error_handler (ncode=istat) istat = NF90_GET_VAR (ncid_in, varid_in(2), lat) if (istat /= NF90_NOERR) call error_handler (ncode=istat) istat = NF90_GET_VAR (ncid_in, varid_in(3), a_bnds) if (istat /= NF90_NOERR) call error_handler (ncode=istat) istat = NF90_GET_VAR (ncid_in, varid_in(4), b_bnds) if (istat /= NF90_NOERR) call error_handler (ncode=istat) !----------------------------------- ! setup output file istat = NF90_CREATE (trim(file_out), NF90_CLOBBER, ncid_out) if (istat /= NF90_NOERR) call error_handler (ncode=istat) istat = NF90_DEF_DIM (ncid_out, 'lon', nlon, dimid_out(1)) istat = NF90_DEF_DIM (ncid_out, 'lat', nlat, dimid_out(2)) istat = NF90_DEF_DIM (ncid_out, 'plev', nplev, dimid_out(3)) istat = NF90_DEF_DIM (ncid_out, 'time', NF90_UNLIMITED, dimid_out(4)) ! output variables istat = NF90_DEF_VAR (ncid_out, 'lon', NF90_DOUBLE, dimid_out(1:1), varid_out(1)) istat = NF90_COPY_ATT (ncid_in, varid_in(1), 'units', ncid_out, varid_out(1)) istat = NF90_COPY_ATT (ncid_in, varid_in(1), 'long_name', ncid_out, varid_out(1)) istat = NF90_COPY_ATT (ncid_in, varid_in(1), 'standard_name', ncid_out, varid_out(1)) istat = NF90_COPY_ATT (ncid_in, varid_in(1), 'axis', ncid_out, varid_out(1)) istat = NF90_DEF_VAR (ncid_out, 'lat', NF90_DOUBLE, dimid_out(2:2), varid_out(2)) istat = NF90_COPY_ATT (ncid_in, varid_in(2), 'units', ncid_out, varid_out(2)) istat = NF90_COPY_ATT (ncid_in, varid_in(2), 'long_name', ncid_out, varid_out(2)) istat = NF90_COPY_ATT (ncid_in, varid_in(2), 'standard_name', ncid_out, varid_out(2)) istat = NF90_COPY_ATT (ncid_in, varid_in(2), 'axis', ncid_out, varid_out(2)) istat = NF90_DEF_VAR (ncid_out, 'plev', NF90_DOUBLE, dimid_out(3:3), varid_out(3)) istat = NF90_PUT_ATT (ncid_out, varid_out(3), 'units', 'Pa') istat = NF90_PUT_ATT (ncid_out, varid_out(3), 'long_name', 'pressure') istat = NF90_PUT_ATT (ncid_out, varid_out(3), 'standard_name', 'air_pressure') istat = NF90_PUT_ATT (ncid_out, varid_out(3), 'axis', 'Z') istat = NF90_DEF_VAR (ncid_out, 'time', NF90_DOUBLE, dimid_out(4:4), varid_out(4)) istat = NF90_COPY_ATT (ncid_in, varid_in(5), 'units', ncid_out, varid_out(4)) istat = NF90_COPY_ATT (ncid_in, varid_in(5), 'long_name', ncid_out, varid_out(4)) istat = NF90_COPY_ATT (ncid_in, varid_in(5), 'standard_name', ncid_out, varid_out(4)) istat = NF90_COPY_ATT (ncid_in, varid_in(5), 'axis', ncid_out, varid_out(4)) istat = NF90_COPY_ATT (ncid_in, varid_in(5), 'calendar', ncid_out, varid_out(4)) istat = NF90_DEF_VAR (ncid_out, trim(fieldname), NF90_REAL4, dimid_out, varid_out(5)) istat = NF90_COPY_ATT (ncid_in, varid_in(7), 'units', ncid_out, varid_out(5)) istat = NF90_COPY_ATT (ncid_in, varid_in(7), 'long_name', ncid_out, varid_out(5)) istat = NF90_COPY_ATT (ncid_in, varid_in(7), 'standard_name', ncid_out, varid_out(5)) istat = NF90_COPY_ATT (ncid_in, varid_in(7), 'cell_methods', ncid_out, varid_out(5)) istat = NF90_PUT_ATT (ncid_out, varid_out(5), '_FillValue', NF90_FILL_REAL4) istat = NF90_ENDDEF (ncid_out) ! write static data istat = NF90_PUT_VAR (ncid_out, varid_out(1), lon) if (istat /= NF90_NOERR) call error_handler (ncode=istat) istat = NF90_PUT_VAR (ncid_out, varid_out(2), lat) if (istat /= NF90_NOERR) call error_handler (ncode=istat) istat = NF90_PUT_VAR (ncid_out, varid_out(3), plev) if (istat /= NF90_NOERR) call error_handler (ncode=istat) ! preparatory calculation do kp = 1, nplev lnpout(kp) = log(plev(kp)) enddo ntop = 0 nbot = 0 !----------------------------------- ! loop through the record dimension do n = 1, numtime istat = NF90_GET_VAR (ncid_in, varid_in(5), time, start=(/n/)) ! read surface pressure istat = NF90_GET_VAR (ncid_in, varid_in(6), ps, start=(/1,1,n/)) ! read 3d field istat = NF90_GET_VAR (ncid_in, varid_in(7), data_in, start=(/1,1,1,n/)) ! interpolate do j = 1, nlat do i = 1, nlon do k = 1, nlev p_bnds(1:2) = a_bnds(1:2,k) + b_bnds(1:2,k) * ps(i,j) lnp(k) = log((p_bnds(2)-p_bnds(1))/log(p_bnds(2)/p_bnds(1))) enddo ! loop through pressure levels do kp = 1, nplev ! find model levels km = nlev do k = 2, nlev if (lnpout(kp) <= lnp(k)) then km = k !mask = .false. exit endif enddo ! k fact = (lnpout(kp)-lnp(km))/ (lnp(km-1)-lnp(km)) if (fact >= 0.0 .and. fact <= 1.0) then ! interpolation data_out(i,j,kp) = data_in(i,j,km) + fact * (data_in(i,j,km-1)-data_in(i,j,km)) else if (fact > 1.0) then ! extrapolation above top level if (xtop) then data_out(i,j,kp) = data_in(i,j,km) + fact * (data_in(i,j,km-1)-data_in(i,j,km)) else data_out(i,j,kp) = NF90_FILL_REAL4 endif ntop = ntop+1 else if (fact < 0.0) then ! extrapolation below bottom level if (xbot) then data_out(i,j,kp) = data_in(i,j,km) + fact * (data_in(i,j,km-1)-data_in(i,j,km)) else data_out(i,j,kp) = NF90_FILL_REAL4 endif nbot = nbot+1 endif enddo ! kp enddo ! i enddo ! j ! write output field istat = NF90_PUT_VAR (ncid_out, varid_out(4), time, start=(/n/)) istat = NF90_PUT_VAR (ncid_out, varid_out(5), data_out, start=(/1,1,1,n/)) enddo istat = NF90_CLOSE (ncid_in) istat = NF90_CLOSE (ncid_out) write(pct,'(f6.2)') real(100*nbot)/real(nlon*nlat*nplev*numtime) print *, 'number of points above top level = ', ntop print *, 'number of points below bottom level = ', nbot, '('//trim(adjustl(pct))//'%)' !--------- END OF PROGRAM --------- contains subroutine error_handler (string, ncode) character(len=*), intent(in), optional :: string integer , intent(in), optional :: ncode character(len=80) :: errstrg if (present(ncode)) then print *, 'NETCDF ERROR' else print *, 'ERROR' endif if (present(string)) print *, trim(string) if (present(ncode)) then errstrg = NF90_STRERROR (ncode) print *, trim(errstrg) endif stop 111 end subroutine error_handler end program plevel