C +-======-+ C Copyright (c) 2003-2007 United States Government as represented by C the Admistrator of the National Aeronautics and Space Administration. C All Rights Reserved. C C THIS OPEN SOURCE AGREEMENT ("AGREEMENT") DEFINES THE RIGHTS OF USE, C REPRODUCTION, DISTRIBUTION, MODIFICATION AND REDISTRIBUTION OF CERTAIN C COMPUTER SOFTWARE ORIGINALLY RELEASED BY THE UNITED STATES GOVERNMENT AS C REPRESENTED BY THE GOVERNMENT AGENCY LISTED BELOW ("GOVERNMENT AGENCY"). C THE UNITED STATES GOVERNMENT, AS REPRESENTED BY GOVERNMENT AGENCY, IS AN C INTENDED THIRD-PARTY BENEFICIARY OF ALL SUBSEQUENT DISTRIBUTIONS OR C REDISTRIBUTIONS OF THE SUBJECT SOFTWARE. ANYONE WHO USES, REPRODUCES, C DISTRIBUTES, MODIFIES OR REDISTRIBUTES THE SUBJECT SOFTWARE, AS DEFINED C HEREIN, OR ANY PART THEREOF, IS, BY THAT ACTION, ACCEPTING IN FULL THE C RESPONSIBILITIES AND OBLIGATIONS CONTAINED IN THIS AGREEMENT. C C Government Agency: National Aeronautics and Space Administration C Government Agency Original Software Designation: GSC-15354-1 C Government Agency Original Software Title: GEOS-5 GCM Modeling Software C User Registration Requested. Please Visit http://opensource.gsfc.nasa.gov C Government Agency Point of Contact for Original Software: C Dale Hithon, SRA Assistant, (301) 286-2691 C C +-======-+ program main implicit none include 'alias.h' c ********************************************************************** c ********************************************************************** c **** **** c **** Program to create prs output from an eta file **** c **** **** c ********************************************************************** c ********************************************************************** integer im,jm,lm,nt integer nymd,nhms integer im_out, jm_out integer nymd0,nhms0,hour,day,month,year integer nymdb,nhmsb c Generic Model Variables c ----------------------- real, allocatable :: ps(:,:) real, allocatable :: dp(:,:,:) real, allocatable :: q2d(:,:,:) real, allocatable :: q3d(:,:,:,:) c HDF and other Local Variables c ----------------------------- logical, pointer :: Lsurf (:) character*128, pointer :: names (:) character*128, pointer :: name2d(:), name3d(:) character*128, pointer :: titl2d(:), titl3d(:) integer id,rc,fid,nhmsf,n2d,n3d integer nvars,ngatts,ntime,ntimes,gfrc real, allocatable :: plevs(:) character*128, allocatable :: arg(:) character*128, allocatable :: fname(:) character*128 dummy, name character*128 output, hdfile, ctlfile, output2 character*8 date,date0 character*2 time,time0 character*1 char integer n,m,nargs,iargc,L,nbeg,nfiles,mlev real*8 lonbeg real undef, getcon real, allocatable :: lat(:) real, allocatable :: lon(:) real, allocatable :: lev(:) real, allocatable :: vrange(:,:) real, allocatable :: prange(:,:) integer, allocatable :: yymmdd(:) integer, allocatable :: hhmmss(:) integer, allocatable :: kmvar(:) character*256 title character*256 source character*256 contact character*256 levunits character*256, allocatable :: vname(:) character*256, allocatable :: vtitle(:) character*256, allocatable :: vunits(:) integer i,j,ndt integer imax,jmax logical hdf, quad logical hdfcreate logical ctl_exists logical edges character*8 cdate interface subroutine read_ctl ( ctlfile,im,jm,lm,n2d,n3d,lonbeg,undef, . nvars,names,Lsurf,name2d,titl2d,name3d,titl3d ) logical, pointer :: Lsurf (:) character*128, pointer :: names (:) character*128, pointer :: name2d(:), name3d(:) character*128, pointer :: titl2d(:), titl3d(:) character*128 ctlfile integer im,jm,lm,n2d,n3d,nvars real undef real*8 lonbeg end subroutine read_ctl subroutine read_hdf ( hdffile,im,jm,lm,n2d,n3d,lonbeg,undef,id, . nymdb,nhmsb,ndt,ntimes, . nvars,names,Lsurf,name2d,titl2d,name3d,titl3d ) logical, pointer :: Lsurf (:) character*128, pointer :: names (:) character*128, pointer :: name2d(:), name3d(:) character*128, pointer :: titl2d(:), titl3d(:) character*128 hdffile integer id,im,jm,lm,n2d,n3d,nvars integer nymdb,nhmsb,ndt,ntimes real undef real*8 lonbeg end subroutine read_hdf end interface C ********************************************************************** C **** Initialization **** C ********************************************************************** ctlfile = 'xxx' im_out = -999 jm_out = -999 nymd0 = -999 nhms0 = -999 nt = 1 hdf = .true. quad = .false. output ='eta2prs' output2 ='NONE' nargs = iargc() if( nargs.eq.0 ) then call usage() else allocate ( arg(nargs) ) do n=1,nargs call getarg(n,arg(n)) enddo do n=1,nargs if( trim(arg(n)).eq.'-im' ) read(arg(n+1),*) im_out if( trim(arg(n)).eq.'-jm' ) read(arg(n+1),*) jm_out if( trim(arg(n)).eq.'-nymd' ) read(arg(n+1),*) nymd0 if( trim(arg(n)).eq.'-nhms' ) read(arg(n+1),*) nhms0 if( trim(arg(n)).eq.'-ndt' ) read(arg(n+1),*) ndt if( trim(arg(n)).eq.'-hdf' ) read(arg(n+1),*) hdf if( trim(arg(n)).eq.'-quad' ) quad = .true. if( trim(arg(n)).eq.'-flat' ) hdf = .false. if( trim(arg(n)).eq.'-tag' ) output = arg(n+1) if( trim(arg(n)).eq.'-ctl' ) ctlfile = arg(n+1) if( trim(arg(n)).eq.'-o' ) output2 = arg(n+1) if( trim(arg(n)).eq.'-levs' ) then mlev = 1 read(arg(n+mlev),fmt='(a1)') char do while (char.ne.'-' .and. n+mlev.lt.nargs ) mlev = mlev+1 read(arg(n+mlev),fmt='(a1)') char enddo if( char.eq.'-' ) mlev = mlev-1 allocate ( plevs(mlev) ) do m=1,mlev read(arg(n+m),*) plevs(m) enddo endif if( trim(arg(n)).eq.'-eta' ) then nfiles = 1 read(arg(n+nfiles),fmt='(a1)') char do while (char.ne.'-' .and. n+nfiles.ne.nargs ) nfiles = nfiles+1 read(arg(n+nfiles),fmt='(a1)') char enddo if( char.eq.'-' ) nfiles = nfiles-1 allocate ( fname(nfiles) ) do m=1,nfiles fname(m) = arg(n+m) enddo endif enddo endif C ********************************************************************** C **** Read Grads CLT or HDF Meta Data **** C ********************************************************************** ! Check whether ctl file exists ! ----------------------------- inquire ( file=trim(ctlfile), exist=ctl_exists ) if( ctl_exists ) then call read_ctl ( ctlfile,im,jm,lm,n2d,n3d,lonbeg,undef, . nvars,names,Lsurf,name2d,titl2d,name3d,titl3d ) else call read_hdf ( fname(1),im,jm,lm,n2d,n3d,lonbeg,undef,id, . nymdb,nhmsb,ndt,ntimes, . nvars,names,Lsurf,name2d,titl2d,name3d,titl3d ) endif C ********************************************************************** C **** Summarize Input Variables **** C ********************************************************************** allocate ( ps(im,jm) ) allocate ( dp(im,jm,lm) ) allocate ( q2d(im,jm, n2d) ) allocate ( q3d(im,jm,lm,n3d) ) if( im_out.eq.-999 ) im_out = im if( jm_out.eq.-999 ) jm_out = jm if( nymd0 == -999 ) nymd0 = nymdb if( nhms0 == -999 ) nhms0 = nhmsb print * print *, ' im: ',im_out print *, ' jm: ',jm_out print *, ' lm: ',lm print *, 'Beginning Date: ',nymd0 print *, 'Beginning Time: ',nhms0 print *, 'Time Increment: ',nhmsf(ndt),' (',ndt,' seconds)' print *, ' plevs: ',(plevs(L),L=1,mlev) print * print *, '2-D Fields:' do n=1,n2d print *, trim(name2d(n)),' ',trim(titl2d(n)) enddo print * print *, '3-D Fields:' do n=1,n3d print *, trim(name3d(n)),' ',trim(titl3d(n)) enddo print * print *, 'Files: ' do n=1,nfiles print *, n,trim(fname(n)) enddo print * plevs(1:mlev) = plevs(mlev:1:-1)*100 C ********************************************************************** C **** Read and Interpolate Eta File **** C ********************************************************************** nymd = nymd0 nhms = nhms0 edges = .false. do n=1,nfiles print *, 'Opening: ',trim(fname(n)) write(date0,1000) nymd write(time0,2000) nhms/10000 1000 format(i8.8) 2000 format(i2.2) if(hdf) then hdfile = trim(output) // "." // trim(date0) // "_" // trim(time0) // "z.nc4" else hdfile = trim(output) // "." // trim(date0) // "_" // trim(time0) // "z.bin" endif if ( trim(output2) .ne. 'NONE' ) hdfile = trim(output2) ! use user-specified fname if( ctl_exists ) then close(10) open (10,file=trim(fname(n)),form='unformatted',access='sequential') else call gfio_close ( id,rc ) call gfio_open ( trim(fname(n)),1,id,rc ) endif rc = 0 ntime = 0 hdfcreate = .true. print * dowhile (rc.eq.0) ntime = ntime + 1 print *, 'nymd: ',nymd,' nhms: ',nhms if( ctl_exists ) then call read_flat_data ( 10,dp,ps,q2d,q3d,names,nvars,Lsurf,im,jm,lm,nt,nymd,nhms,rc ) else if( nymd.eq.nymd0 .and. nhms.eq.nhms0 ) then call read_hdf_data ( id,dp,ps,q2d,q3d,n2d,n3d,name2d,name3d, . im,jm,lm,nymdb,nhmsb,rc,ntime,ntimes,edges ) else call read_hdf_data ( id,dp,ps,q2d,q3d,n2d,n3d,name2d,name3d, . im,jm,lm,nymd,nhms,rc,ntime,ntimes,edges ) endif endif if( rc.eq.0 ) then call eta2prs ( dp,ps,q2d,q3d,name2d,name3d,titl2d,titl3d,n2d,n3d,undef, . im,jm,lm,nt,plevs,mlev,im_out,jm_out,lonbeg,nymd,nhms,ndt, . fid,hdf,hdfcreate,hdfile,quad,edges ) call tick (nymd,nhms,ndt) hdfcreate = .false. else if(hdf) then call gfio_close ( fid,gfrc ) else close(55) endif print * print *, 'Created: ',trim(hdfile) print * endif enddo enddo deallocate ( dp,ps,arg ) stop end subroutine read_hdf_data ( id,dp,ps,q2d,q3d,n2d,n3d,name2d,name3d, . im,jm,lm,nymd,nhms,rc,ntime,ntimes,edges ) implicit none include 'alias.h' integer im,jm,lm,nymd,nhms,id,rc integer n2d,n3d,ntime,ntimes real ps(im,jm) real dp(im,jm,lm) real q2d(im,jm ,n2d) real q3d(im,jm,lm,n3d) character*128 name2d(n2d) character*128 name3d(n3d) logical edges integer i,j,L,n real, allocatable :: pl(:,:,:) ps = -999 dp = -999 rc = 0 if( ntime <= ntimes ) then do n=1,n2d if( match( c_ps,name2d(n) ) ) then call gfio_getvar ( id,trim(name2d(n)),nymd,nhms,im,jm,0, 1,ps,rc ) endif enddo edges = .false. if( rc.eq.0 ) then do n=1,n3d if( match( c_dp,name3d(n) ) ) then call gfio_getvar ( id,trim(name3d(n)),nymd,nhms,im,jm,1,lm,dp,rc ) endif if( match( c_pl,name3d(n) ) ) then allocate( pl(im,jm,lm) ) call gfio_getvar ( id,trim(name3d(n)),nymd,nhms,im,jm,1,lm,pl,rc ) dp(:,:,lm) = 2*( ps(:,:)-pl(:,:,lm) ) do L=lm-1,1,-1 dp(:,:,L) = 2*( pl(:,:,L+1)-0.5*dp(:,:,L+1)-pl(:,:,L) ) enddo deallocate( pl ) endif if( match( c_ple,name3d(n) ) ) then edges = .true. call gfio_getvar ( id,trim(name3d(n)),nymd,nhms,im,jm,1,lm,dp,rc ) ps(:,:) = dp(:,:,lm) endif enddo do n=1,n2d call gfio_getvar ( id,trim(name2d(n)),nymd,nhms,im,jm,0,1,q2d(1,1,n),rc ) enddo do n=1,n3d call gfio_getvar ( id,trim(name3d(n)),nymd,nhms,im,jm,1,lm,q3d(1,1,1,n),rc ) enddo do L=1,lm do j=1,jm do i=1,im if( ps(i,j).eq.-999 .or. dp(i,j,L).eq.-999 ) then print *, 'PS and DelP were not created!' stop endif enddo enddo enddo endif else rc = 1 ! No more time periods in file endif return end subroutine read_hdf_data subroutine read_flat_data ( ku,dp,ps,q2d,q3d,names,nvars, . Lsurf,im,jm,lm,nt,nymd,nhms,rc ) implicit none include 'alias.h' integer im,jm,lm,nymd,nhms,ku,rc integer nt,nvars real ps(im,jm) real dp(im,jm,lm) real q2d(im,jm ,1) real q3d(im,jm,lm,1) real ple(im,jm,lm+1) logical Lsurf(nvars), found_ps, found_dp, found_ple real sum,grav,getcon,ptop,kappa,pk0,rgas,rvap,eps character*8 date,date0 character*128 filename,names(nvars) save date0 integer i,j,L,k,m,n c Test for End of File c -------------------- print *, 'Checking for More Data' rc = 0 read(ku,err=999,end=999) backspace(ku) goto 1000 999 continue print *, 'End of File reached' rc = -1 return 1000 continue print *, ' More Data Found' C ********************************************************************** C **** Read Diagnostic Data **** C ********************************************************************** found_ps = .false. found_dp = .false. found_ple = .false. m = 0 n = 0 do k=1,nvars if( match( c_ps,names(k) ) ) then call readit (ps,im,jm, 1,ku) found_ps = .true. m = m+1 q2d(:,:,m) = ps(:,:) else if( match( c_dp,names(k) ) ) then call readit (dp,im,jm,lm,ku) found_dp = .true. n = n+1 q3d(:,:,:,n) = dp(:,:,:) else if( match( c_ple,names(k) ) ) then call readit (ple,im,jm,lm,ku) found_ple = .true. n = n+1 q3d(:,:,:,n) = ple(:,:,:) else if( Lsurf(k) ) then m = m+1 call readit (q2d(1,1,m),im,jm,1,ku) else n = n+1 call readit (q3d(1,1,1,n),im,jm,lm,ku) endif endif enddo if( found_ps .and. found_ple .and. .not.found_dp ) then ple(:,:,lm+1) = ps(:,:) do L=1,lm dp(:,:,L) = ple(:,:,L+1)-ple(:,:,L) enddo endif if( .not.found_ps .or. ( .not.found_dp .and. .not.found_ple ) ) then print *, 'Error!' print *, 'Program requires PLE or PS and DELP to be present' stop endif return end subroutine read_flat_data subroutine readit (q,im,jm,lm,ku) implicit none integer im,jm,lm,ku,L real q(im,jm,lm) real*4 dum(im,jm) do L=1,lm read(ku) dum q(:,:,L) = dum(:,:) enddo return end subroutine readit subroutine eta2prs ( dp,ps,q2d,q3d,name2d,name3d,titl2d,titl3d,n2d,n3d,undef, . im,jm,lm,nt,plev,mlev,im_out,jm_out,lonbeg,nymd,nhms,ninc, . id,hdf,create,filename,quad,edges ) implicit none include 'alias.h' c Input Variables c --------------- integer im,jm,lm,nt,mlev,im_out,jm_out,nymd,nhms,ninc,n2d,n3d integer nymd0,nhms0 real dp(im,jm,lm) real ps(im,jm) real q2d(im,jm, n2d) real q3d(im,jm,lm,n3d) character*128 name2d(n2d), titl2d(n2d) character*128 name3d(n3d), titl3d(n3d) character*128 filename real plev(mlev) ! Target Pressures for Output real*8 lonbeg, dlon, dlat real*8 lat(jm_out),lon(im_out) logical hdf, create, quad, edges c Local Variables c --------------- integer i,j,k,L,n,m real logpl(im,jm,lm) real logps(im,jm) real logpt(im,jm) real, allocatable :: pk (:,:,:) ! Model pressure at edge-levels real, allocatable :: pke(:,:,:) ! Model pressure at edge-levels real, allocatable :: pe(:,:,:) ! Model pressure at edge-levels real, allocatable :: qprs(:,:,:) ! Interpolated Quantity real lats(jm_out),lons(im_out),levs(mlev) real undef,getcon,eps integer precision,id,timeinc,rc,nhmsf character*128 levunits character*128 title character*128 source character*128 contact integer nvars integer i_u, i_v, i_t, i_q, i_th, i_thv, i_phis character*256, allocatable :: vname(:) , vname2(:) character*256, allocatable :: vtitle(:) , vtitle2(:) character*256, allocatable :: vunits(:) , vunits2(:) integer, allocatable :: lmvar(:) , lmvar2(:) real, allocatable :: vrange(:,:) real, allocatable :: prange(:,:) real pref,tref,pkref,dum,kappa c Value-Added Products c -------------------- real, allocatable :: up(:,:,:), uprime(:,:,:) real, allocatable :: vp(:,:,:), vprime(:,:,:) real, allocatable :: tp(:,:,:), tprime(:,:,:) real, allocatable :: qp(:,:,:) real, allocatable :: hp(:,:,:) real, allocatable :: sphu(:,:,:) real, allocatable :: tmpu(:,:,:) real, allocatable :: th(:,:,:) real, allocatable :: thv(:,:,:) real, allocatable :: phis(:,:) integer nuu, nvv, ntt, nqq, nuv, nut, nuq integer nvt, nvq, nupvp, nvptp, nh C ********************************************************************** C **** Initialize Constants And Local Arrays **** C ********************************************************************** nvars = n2d + n3d kappa = getcon('KAPPA') eps = getcon('EPS') allocate ( qprs(im,jm,mlev) ) allocate ( pe(im,jm,lm+1) ) allocate ( pk(im,jm,lm) ) allocate ( pke(im,jm,lm+1) ) C ********************************************************************** C **** Initialize GFIO File **** C ********************************************************************** allocate ( vname(nvars) ) allocate ( vtitle(nvars) ) allocate ( vunits(nvars) ) allocate ( lmvar(nvars) ) timeinc = nhmsf(ninc) precision = 1 ! 64-bit precision = 0 ! 32-bit title = 'GEOS-5 GCM (Pressure Coordinates)' source = 'Goddard Modeling and Assimilation Office, NASA/GSFC' contact = 'data@gmao.gsfc.nasa.gov' levunits = 'mb' c Defined Fields c -------------- do m=1,n2d n = m vname(n) = name2d(m) vtitle(n) = trim(titl2d(m)) vunits(n) = 'unknown' lmvar(n) = 0 enddo do m=1,n3d n = n2d+m vname(n) = name3d(m) vtitle(n) = trim(titl3d(m)) vunits(n) = 'unknown' lmvar(n) = mlev enddo C ********************************************************************** C **** Value Added Products **** C ********************************************************************** i_phis = 0 i_u = 0 i_v = 0 i_t = 0 i_q = 0 i_th = 0 i_thv = 0 do n=1,n2d if( match( c_phis,name2d(n) ) ) i_phis = n enddo do n=1,n3d if( match( c_u ,name3d(n) ) ) i_u = n if( match( c_v ,name3d(n) ) ) i_v = n if( match( c_t ,name3d(n) ) ) i_t = n if( match( c_q ,name3d(n) ) ) i_q = n if( match( c_th ,name3d(n) ) ) i_th = n if( match( c_thv ,name3d(n) ) ) i_thv = n enddo c uwnd*uwnd c --------- if( quad .and. i_u.ne.0 ) then allocate ( up(im,jm,mlev), uprime(im,jm,mlev) ) allocate ( vname2(nvars), vtitle2(nvars), vunits2(nvars), lmvar2(nvars) ) vname2 = vname vtitle2 = vtitle vunits2 = vunits lmvar2 = lmvar nvars = nvars + 1 nuu = nvars deallocate ( vname, vtitle, vunits, lmvar ) allocate ( vname(nvars), vtitle(nvars), vunits(nvars), lmvar(nvars) ) vname(1:nvars-1) = vname2 vtitle(1:nvars-1) = vtitle2 vunits(1:nvars-1) = vunits2 lmvar(1:nvars-1) = lmvar2 vname( nvars ) = 'uu' vtitle( nvars ) = 'Quadratic UWND*UWND' vunits( nvars ) = 'm**2/s**2' lmvar( nvars ) = mlev deallocate ( vname2, vtitle2, vunits2, lmvar2 ) endif c vwnd*vwnd c --------- if( quad .and. i_v.ne.0 ) then allocate ( vp(im,jm,mlev), vprime(im,jm,mlev) ) allocate ( vname2(nvars), vtitle2(nvars), vunits2(nvars), lmvar2(nvars) ) vname2 = vname vtitle2 = vtitle vunits2 = vunits lmvar2 = lmvar nvars = nvars + 1 nvv = nvars deallocate ( vname, vtitle, vunits, lmvar ) allocate ( vname(nvars), vtitle(nvars), vunits(nvars), lmvar(nvars) ) vname(1:nvars-1) = vname2 vtitle(1:nvars-1) = vtitle2 vunits(1:nvars-1) = vunits2 lmvar(1:nvars-1) = lmvar2 vname( nvars ) = 'vv' vtitle( nvars ) = 'Quadratic VWND*VWND' vunits( nvars ) = 'm**2/s**2' lmvar( nvars ) = mlev deallocate ( vname2, vtitle2, vunits2, lmvar2 ) endif c tmpu*tmpu c --------- if( quad .and. i_t.ne.0 ) then allocate ( tp(im,jm,mlev), tprime(im,jm,mlev) ) allocate ( vname2(nvars), vtitle2(nvars), vunits2(nvars), lmvar2(nvars) ) vname2 = vname vtitle2 = vtitle vunits2 = vunits lmvar2 = lmvar nvars = nvars + 1 ntt = nvars deallocate ( vname, vtitle, vunits, lmvar ) allocate ( vname(nvars), vtitle(nvars), vunits(nvars), lmvar(nvars) ) vname(1:nvars-1) = vname2 vtitle(1:nvars-1) = vtitle2 vunits(1:nvars-1) = vunits2 lmvar(1:nvars-1) = lmvar2 vname( nvars ) = 'tt' vtitle( nvars ) = 'Quadratic TMPU*TMPU' vunits( nvars ) = 'K**2' lmvar( nvars ) = mlev deallocate ( vname2, vtitle2, vunits2, lmvar2 ) endif c sphu*sphu c --------- if( quad .and. i_q.ne.0 ) then allocate ( qp(im,jm,mlev) ) allocate ( vname2(nvars), vtitle2(nvars), vunits2(nvars), lmvar2(nvars) ) vname2 = vname vtitle2 = vtitle vunits2 = vunits lmvar2 = lmvar nvars = nvars + 1 nqq = nvars deallocate ( vname, vtitle, vunits, lmvar ) allocate ( vname(nvars), vtitle(nvars), vunits(nvars), lmvar(nvars) ) vname(1:nvars-1) = vname2 vtitle(1:nvars-1) = vtitle2 vunits(1:nvars-1) = vunits2 lmvar(1:nvars-1) = lmvar2 vname( nvars ) = 'qq' vtitle( nvars ) = 'Quadratic SPHU*SPHU' vunits( nvars ) = 'nondim' lmvar( nvars ) = mlev deallocate ( vname2, vtitle2, vunits2, lmvar2 ) endif c uwnd*vwnd c --------- if( quad .and. i_u.ne.0 .and. i_v.ne.0 ) then allocate ( vname2(nvars), vtitle2(nvars), vunits2(nvars), lmvar2(nvars) ) vname2 = vname vtitle2 = vtitle vunits2 = vunits lmvar2 = lmvar nvars = nvars + 1 nuv = nvars deallocate ( vname, vtitle, vunits, lmvar ) allocate ( vname(nvars), vtitle(nvars), vunits(nvars), lmvar(nvars) ) vname(1:nvars-1) = vname2 vtitle(1:nvars-1) = vtitle2 vunits(1:nvars-1) = vunits2 lmvar(1:nvars-1) = lmvar2 vname( nvars ) = 'uv' vtitle( nvars ) = 'Quadratic UWND*VWND' vunits( nvars ) = 'm**2/sec**2' lmvar( nvars ) = mlev deallocate ( vname2, vtitle2, vunits2, lmvar2 ) endif c uwnd*tmpu c --------- if( quad .and. i_u.ne.0 .and. i_t.ne.0 ) then allocate ( vname2(nvars), vtitle2(nvars), vunits2(nvars), lmvar2(nvars) ) vname2 = vname vtitle2 = vtitle vunits2 = vunits lmvar2 = lmvar nvars = nvars + 1 nut = nvars deallocate ( vname, vtitle, vunits, lmvar ) allocate ( vname(nvars), vtitle(nvars), vunits(nvars), lmvar(nvars) ) vname(1:nvars-1) = vname2 vtitle(1:nvars-1) = vtitle2 vunits(1:nvars-1) = vunits2 lmvar(1:nvars-1) = lmvar2 vname( nvars ) = 'ut' vtitle( nvars ) = 'Quadratic UWND*TMPU' vunits( nvars ) = 'K m/sec' lmvar( nvars ) = mlev deallocate ( vname2, vtitle2, vunits2, lmvar2 ) endif c uwnd*sphu c --------- if( quad .and. i_u.ne.0 .and. i_q.ne.0 ) then allocate ( vname2(nvars), vtitle2(nvars), vunits2(nvars), lmvar2(nvars) ) vname2 = vname vtitle2 = vtitle vunits2 = vunits lmvar2 = lmvar nvars = nvars + 1 nuq = nvars deallocate ( vname, vtitle, vunits, lmvar ) allocate ( vname(nvars), vtitle(nvars), vunits(nvars), lmvar(nvars) ) vname(1:nvars-1) = vname2 vtitle(1:nvars-1) = vtitle2 vunits(1:nvars-1) = vunits2 lmvar(1:nvars-1) = lmvar2 vname( nvars ) = 'uq' vtitle( nvars ) = 'Quadratic UWND*SPHU' vunits( nvars ) = 'm/sec' lmvar( nvars ) = mlev deallocate ( vname2, vtitle2, vunits2, lmvar2 ) endif c vwnd*tmpu c --------- if( quad .and. i_v.ne.0 .and. i_t.ne.0 ) then allocate ( vname2(nvars), vtitle2(nvars), vunits2(nvars), lmvar2(nvars) ) vname2 = vname vtitle2 = vtitle vunits2 = vunits lmvar2 = lmvar nvars = nvars + 1 nvt = nvars deallocate ( vname, vtitle, vunits, lmvar ) allocate ( vname(nvars), vtitle(nvars), vunits(nvars), lmvar(nvars) ) vname(1:nvars-1) = vname2 vtitle(1:nvars-1) = vtitle2 vunits(1:nvars-1) = vunits2 lmvar(1:nvars-1) = lmvar2 vname( nvars ) = 'vt' vtitle( nvars ) = 'Quadratic VWND*TMPU' vunits( nvars ) = 'K m/sec' lmvar( nvars ) = mlev deallocate ( vname2, vtitle2, vunits2, lmvar2 ) endif c vwnd*sphu c --------- if( quad .and. i_v.ne.0 .and. i_q.ne.0 ) then allocate ( vname2(nvars), vtitle2(nvars), vunits2(nvars), lmvar2(nvars) ) vname2 = vname vtitle2 = vtitle vunits2 = vunits lmvar2 = lmvar nvars = nvars + 1 nvq = nvars deallocate ( vname, vtitle, vunits, lmvar ) allocate ( vname(nvars), vtitle(nvars), vunits(nvars), lmvar(nvars) ) vname(1:nvars-1) = vname2 vtitle(1:nvars-1) = vtitle2 vunits(1:nvars-1) = vunits2 lmvar(1:nvars-1) = lmvar2 vname( nvars ) = 'vq' vtitle( nvars ) = 'Quadratic VWND*SPHU' vunits( nvars ) = 'm/sec' lmvar( nvars ) = mlev deallocate ( vname2, vtitle2, vunits2, lmvar2 ) endif c uprime*vprime c ------------- if( quad .and. i_u.ne.0 .and. i_v.ne.0 ) then allocate ( vname2(nvars), vtitle2(nvars), vunits2(nvars), lmvar2(nvars) ) vname2 = vname vtitle2 = vtitle vunits2 = vunits lmvar2 = lmvar nvars = nvars + 1 nupvp = nvars deallocate ( vname, vtitle, vunits, lmvar ) allocate ( vname(nvars), vtitle(nvars), vunits(nvars), lmvar(nvars) ) vname(1:nvars-1) = vname2 vtitle(1:nvars-1) = vtitle2 vunits(1:nvars-1) = vunits2 lmvar(1:nvars-1) = lmvar2 vname( nvars ) = 'upvp' vtitle( nvars ) = 'Quadratic UPRIME*VPRIME (Departure from zonal mean)' vunits( nvars ) = 'm**2/sec**2' lmvar( nvars ) = mlev deallocate ( vname2, vtitle2, vunits2, lmvar2 ) endif c vprime*tprime c ------------- if( quad .and. i_v.ne.0 .and. i_t.ne.0 ) then allocate ( vname2(nvars), vtitle2(nvars), vunits2(nvars), lmvar2(nvars) ) vname2 = vname vtitle2 = vtitle vunits2 = vunits lmvar2 = lmvar nvars = nvars + 1 nvptp = nvars deallocate ( vname, vtitle, vunits, lmvar ) allocate ( vname(nvars), vtitle(nvars), vunits(nvars), lmvar(nvars) ) vname(1:nvars-1) = vname2 vtitle(1:nvars-1) = vtitle2 vunits(1:nvars-1) = vunits2 lmvar(1:nvars-1) = lmvar2 vname( nvars ) = 'vptp' vtitle( nvars ) = 'Quadratic VPRIME*TPRIME (Departure from zonal mean)' vunits( nvars ) = 'K m/sec' lmvar( nvars ) = mlev deallocate ( vname2, vtitle2, vunits2, lmvar2 ) endif c heights c ------- if( i_phis.ne.0 .and. . ( (i_t .ne.0 .and. i_q.ne.0) .or. . (i_th.ne.0 .and. i_q.ne.0) .or. i_thv.ne.0 ) ) then allocate ( phis(im,jm), hp(im,jm,mlev) ) allocate ( vname2(nvars), vtitle2(nvars), vunits2(nvars), lmvar2(nvars) ) vname2 = vname vtitle2 = vtitle vunits2 = vunits lmvar2 = lmvar nvars = nvars + 1 nh = nvars deallocate ( vname, vtitle, vunits, lmvar ) allocate ( vname(nvars), vtitle(nvars), vunits(nvars), lmvar(nvars) ) vname(1:nvars-1) = vname2 vtitle(1:nvars-1) = vtitle2 vunits(1:nvars-1) = vunits2 lmvar(1:nvars-1) = lmvar2 vname( nvars ) = 'hght' vtitle( nvars ) = 'Heights' vunits( nvars ) = 'm' lmvar( nvars ) = mlev deallocate ( vname2, vtitle2, vunits2, lmvar2 ) endif c Compute grid c ------------ dlon = 360.0/ im_out dlat = 180.0/(jm_out-1) do j=1,jm_out lat(j) = -90.0 + (j-1)*dlat enddo do i=1,im_out lon(i) = lonbeg + (i-1)*dlon enddo lons = lon lats = lat levs(:) = plev(:)/100 c Create GFIO file c ---------------- allocate ( vrange(2,nvars) ) allocate ( prange(2,nvars) ) vrange(:,:) = undef prange(:,:) = undef if (create) then if (hdf) then call GFIO_Create ( trim(filename), title, source, contact, undef, . im_out, jm_out, mlev, lons, lats, levs, levunits, . nymd, nhms, timeinc, . nvars, vname, vtitle, vunits, lmvar, . vrange, prange, precision, . id, rc ) else open (55,file=trim(filename),form='unformatted',access='sequential') endif endif C ********************************************************************** C **** Compute edge-level pressures **** C ********************************************************************** if( .not.edges) then if( ps(1,1).ne.0.0 ) then pe (:,:,lm+1) = ps(:,:) do L=lm,1,-1 pe (:,:,L) = pe(:,:,L+1)-dp(:,:,L) enddo pke = pe**kappa c Compute mid-level pressures c --------------------------- do L=1,lm pk(:,:,L) = ( pke(:,:,L+1)-pke(:,:,L) )/( kappa*log(pe(:,:,L+1)/pe(:,:,L)) ) enddo endif endif C ********************************************************************** C **** Write Defined Fields **** C ********************************************************************** if( n3d.ne.0 ) then if( edges ) then logpl(:,:,:) = log( dp(:,:,: ) ) logps(:,:) = log( dp(:,:,lm) ) logpt(:,:) = log( dp(:,:,1 ) ) else do L=1,lm do j=1,jm do i=1,im logpl(i,j,L) = log( 0.5*( pe(i,j,L+1)+pe(i,j,L) ) ) enddo enddo enddo logps(:,:) = log( pe(:,:,lm+1) ) logpt(:,:) = log( pe(:,:,1) ) endif endif do n=1,n2d call writit( q2d(1,1,n),im,jm,im_out,jm_out,1,id,name2d(n),nymd,nhms,undef,hdf ) if( n.eq.i_phis .and. allocated(phis) ) phis = q2d(:,:,n) enddo do n=1,n3d do L=1,mlev call sigtopl( qprs(1,1,L),q3d(1,1,1,n),logpl,logps,logpt,log(plev(L)),im,jm,lm,undef ) enddo call writit( qprs,im,jm,im_out,jm_out,mlev,id,name3d(n),nymd,nhms,undef,hdf ) if( n.eq.i_u .and. allocated(up) ) up = qprs if( n.eq.i_v .and. allocated(vp) ) vp = qprs if( n.eq.i_t .and. allocated(tp) ) tp = qprs if( n.eq.i_q .and. allocated(qp) ) qp = qprs enddo C ********************************************************************** C **** Compute Heights **** C ********************************************************************** do n=1,n3d if( n.eq.i_t ) then allocate( tmpu(im,jm,lm) ) tmpu = q3d(:,:,:,n) endif if( n.eq.i_q ) then allocate( sphu(im,jm,lm) ) sphu = q3d(:,:,:,n) endif if( n.eq.i_th ) then allocate( th(im,jm,lm) ) th = q3d(:,:,:,n) endif if( n.eq.i_thv ) then allocate( thv(im,jm,lm) ) thv = q3d(:,:,:,n) endif enddo if( allocated(hp) ) then if( allocated(thv) ) then do L=1,mlev call intgeop ( hp(1,1,L),thv,pke(1,1,2),pk,phis,plev(L),undef,im,jm,lm ) enddo else if( allocated(th) .and. allocated(sphu) ) then allocate( thv(im,jm,lm) ) thv = th*(1+eps*sphu) do L=1,mlev call intgeop ( hp(1,1,L),thv,pke(1,1,2),pk,phis,plev(L),undef,im,jm,lm ) enddo else if( allocated(tmpu) .and. allocated(sphu) ) then allocate( thv(im,jm,lm) ) thv = tmpu*(1+eps*sphu)/pk do L=1,mlev call intgeop ( hp(1,1,L),thv,pke(1,1,2),pk,phis,plev(L),undef,im,jm,lm ) enddo endif endif C ********************************************************************** C **** Write value Added Products **** C ********************************************************************** if( quad .and. allocated(up) ) call zprime ( up,uprime,im,jm,mlev,undef ) if( quad .and. allocated(vp) ) call zprime ( vp,vprime,im,jm,mlev,undef ) if( quad .and. allocated(tp) ) call zprime ( tp,tprime,im,jm,mlev,undef ) if( quad .and. allocated(up) ) . call writit2( up, up, im,jm,im_out,jm_out,mlev,id,vname(nuu ),nymd,nhms,undef,hdf ) if( quad .and. allocated(vp) ) . call writit2( vp, vp, im,jm,im_out,jm_out,mlev,id,vname(nvv ),nymd,nhms,undef,hdf ) if( quad .and. allocated(tp) ) . call writit2( tp, tp, im,jm,im_out,jm_out,mlev,id,vname(ntt ),nymd,nhms,undef,hdf ) if( quad .and. allocated(qp) ) . call writit2( qp, qp, im,jm,im_out,jm_out,mlev,id,vname(nqq ),nymd,nhms,undef,hdf ) if( quad .and. allocated(up) .and. . allocated(vp) ) . call writit2( up, vp, im,jm,im_out,jm_out,mlev,id,vname(nuv ),nymd,nhms,undef,hdf ) if( quad .and. allocated(up) .and. . allocated(tp) ) . call writit2( up, tp, im,jm,im_out,jm_out,mlev,id,vname(nut ),nymd,nhms,undef,hdf ) if( quad .and. allocated(up) .and. . allocated(qp) ) . call writit2( up, qp, im,jm,im_out,jm_out,mlev,id,vname(nuq ),nymd,nhms,undef,hdf ) if( quad .and. allocated(vp) .and. . allocated(tp) ) . call writit2( vp, tp, im,jm,im_out,jm_out,mlev,id,vname(nvt ),nymd,nhms,undef,hdf ) if( quad .and. allocated(vp) .and. . allocated(qp) ) . call writit2( vp, qp, im,jm,im_out,jm_out,mlev,id,vname(nvq ),nymd,nhms,undef,hdf ) if( quad .and. allocated(up) .and. . allocated(vp) ) . call writit2( uprime,vprime,im,jm,im_out,jm_out,mlev,id,vname(nupvp),nymd,nhms,undef,hdf ) if( quad .and. allocated(vp) .and. . allocated(tp) ) . call writit2( vprime,tprime,im,jm,im_out,jm_out,mlev,id,vname(nvptp),nymd,nhms,undef,hdf ) if( allocated(hp) ) call writit ( hp,im,jm,im_out,jm_out,mlev,id,vname(nh),nymd,nhms,undef,hdf ) C ********************************************************************** C **** De-Allocate Dynamics Arrays **** C ********************************************************************** if( allocated(up) ) deallocate (up, uprime) if( allocated(vp) ) deallocate (vp, vprime) if( allocated(tp) ) deallocate (tp, tprime) if( allocated(hp) ) deallocate (hp, phis ) if( allocated(tmpu) ) deallocate (tmpu ) if( allocated(sphu) ) deallocate (sphu ) if( allocated(th ) ) deallocate (th ) if( allocated(thv ) ) deallocate (thv ) deallocate (pk,pe,pke) deallocate ( qprs ) deallocate ( vname ) deallocate ( vtitle ) deallocate ( vunits ) deallocate ( lmvar ) deallocate ( vrange ) deallocate ( prange ) return end subroutine eta2prs function defined ( q,undef ) implicit none logical defined real q,undef defined = abs(q-undef).gt.0.1*undef return end function defined function match (name,alias) character*128 name,alias logical match match = .false. if( trim(name) == 'u' ) then if( trim(alias) == 'u' ) match = .true. if( trim(alias) == 'U' ) match = .true. if( trim(alias) == 'uwnd' ) match = .true. if( trim(alias) == 'UWND' ) match = .true. endif if( trim(name) == 'v' ) then if( trim(alias) == 'v' ) match = .true. if( trim(alias) == 'V' ) match = .true. if( trim(alias) == 'vwnd' ) match = .true. if( trim(alias) == 'VWND' ) match = .true. endif if( trim(name) == 't' ) then if( trim(alias) == 't' ) match = .true. if( trim(alias) == 'T' ) match = .true. if( trim(alias) == 'tmpu' ) match = .true. if( trim(alias) == 'TMPU' ) match = .true. endif if( trim(name) == 'q' ) then if( trim(alias) == 'q' ) match = .true. if( trim(alias) == 'Q' ) match = .true. if( trim(alias) == 'qv' ) match = .true. if( trim(alias) == 'QV' ) match = .true. if( trim(alias) == 'sphu' ) match = .true. if( trim(alias) == 'SPHU' ) match = .true. endif if( trim(name) == 'h' ) then if( trim(alias) == 'h' ) match = .true. if( trim(alias) == 'H' ) match = .true. if( trim(alias) == 'hght' ) match = .true. if( trim(alias) == 'HGHT' ) match = .true. endif if( trim(name) == 'th' ) then if( trim(alias) == 'th' ) match = .true. if( trim(alias) == 'TH' ) match = .true. if( trim(alias) == 'theta') match = .true. if( trim(alias) == 'THETA') match = .true. endif if( trim(name) == 'thv' ) then if( trim(alias) == 'thv' ) match = .true. if( trim(alias) == 'THV' ) match = .true. if( trim(alias) == 'thetav') match = .true. if( trim(alias) == 'THETAV') match = .true. endif if( trim(name) == 'ps' ) then if( trim(alias) == 'ps' ) match = .true. if( trim(alias) == 'PS' ) match = .true. endif if( trim(name) == 'dp' ) then if( trim(alias) == 'dp' ) match = .true. if( trim(alias) == 'DP' ) match = .true. if( trim(alias) == 'delp' ) match = .true. if( trim(alias) == 'DELP' ) match = .true. endif if( trim(name) == 'pl' ) then if( trim(alias) == 'pl' ) match = .true. if( trim(alias) == 'PL' ) match = .true. endif if( trim(name) == 'ple' ) then if( trim(alias) == 'ple' ) match = .true. if( trim(alias) == 'PLE' ) match = .true. endif if( trim(name) == 'phis' ) then if( trim(alias) == 'phis' ) match = .true. if( trim(alias) == 'PHIS' ) match = .true. if( trim(alias) == 'gz' ) match = .true. if( trim(alias) == 'GZ' ) match = .true. endif return end FUNCTION GETCON (NAME) C*********************************************************************** C C GENERIC FUNCTION GETCON IS A REPOSITORY OF GLOBAL VARIABLES, I.E. C A MEMORY FOR SCALAR VALUES NEEDED THROUGHOUT A LARGE PROGRAM. C C THIS SPECIFIC FUNCTION, GETCON, REMEMBERS FLOATING POINT VALUES. C THE FUNCTION IS CALLED WITH A CHARACTER NAME TO INTERROGATE A VALUE. C C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C*********************************************************************** integer, PARAMETER :: MAXCON=40 CHARACTER*(*) NAME CHARACTER*16 ANAME(MAXCON) real ACON (MAXCON) integer i C COMPUTATIONAL CONSTANTS C ----------------------- real, PARAMETER :: VECMAX = 65535.5 real, PARAMETER :: CALTOJ = 4184. real, PARAMETER :: UNDEF = -999 C ASTRONOMICAL CONSTANTS C ---------------------- real, PARAMETER :: OB = 23.45 real, PARAMETER :: PER = 102. real, PARAMETER :: ECC = .0167 real, PARAMETER :: AE = 6371E3 real, PARAMETER :: EQNX = 80.5 real, PARAMETER :: SOLS = 176.5 real, PARAMETER :: S0 = 1365.0 C TERRESTRIAL CONSTANTS C --------------------- c PARAMETER ( GRAV = 9.81 ) ! GEOS real, PARAMETER :: GRAV = 9.80616 ! fv real, PARAMETER :: SRFPRS = 984.7 real, PARAMETER :: PIMEAN = 984.7 real, PARAMETER :: PSTD = 1000.0 real, PARAMETER :: TSTD = 280.0 real, PARAMETER :: SDAY = 86400.0 real, PARAMETER :: SSALB = 0.99 real, PARAMETER :: CO2 = 355.0 real, PARAMETER :: CFC11 = 0.3 real, PARAMETER :: CFC12 = 0.5 real, PARAMETER :: CFC22 = 0.2 C THERMODYNAMIC CONSTANTS C ----------------------- c PARAMETER ( CPD = 1004.16 ) ! GEOS real, PARAMETER :: CPD = 1004.64 ! fv real, PARAMETER :: CPV = 1869.46 real, PARAMETER :: ALHL = 2.499E6 real, PARAMETER :: ALHS = 2.845E6 real, PARAMETER :: STFBOL = 5.67E-8 real, PARAMETER :: AIRMW = 28.97 real, PARAMETER :: H2OMW = 18.01 real, PARAMETER :: RUNIV = 8314.3 c PARAMETER ( RGAS = RUNIV/AIRMW) ! GEOS c PARAMETER ( RVAP = RUNIV/H2OMW) ! GEOS real, PARAMETER :: RGAS = 287.04 ! fv real, PARAMETER :: RVAP = 461.00 ! fv real, PARAMETER :: RKAP = RGAS/CPD real, PARAMETER :: HEATW = 597.2 real, PARAMETER :: HEATI = 680.0 real, PARAMETER :: TICE = 273.16 C TURBULENCE CONSTANTS C -------------------- real, PARAMETER :: VKRM = 0.4 C MOISTURE CONSTANTS C ------------------ real, PARAMETER :: EPS = 0.622 real, PARAMETER :: VIRTCON= 0.609 real, PARAMETER :: EPSFAC = EPS*HEATW/RGAS*CALTOJ DATA ANAME(1 ),ACON(1 ) / 'CP ', CPD / DATA ANAME(2 ),ACON(2 ) / 'RGAS ', RGAS / DATA ANAME(3 ),ACON(3 ) / 'KAPPA ', RKAP / DATA ANAME(4 ),ACON(4 ) / 'LATENT HEAT COND', ALHL / DATA ANAME(5 ),ACON(5 ) / 'GRAVITY ', GRAV / DATA ANAME(6 ),ACON(6 ) / 'STEFAN-BOLTZMAN ', STFBOL / DATA ANAME(7 ),ACON(7 ) / 'VON KARMAN ', VKRM / DATA ANAME(8 ),ACON(8 ) / 'EARTH RADIUS ', AE / DATA ANAME(9 ),ACON(9 ) / 'OBLIQUITY ', OB / DATA ANAME(10),ACON(10) / 'ECCENTRICITY ', ECC / DATA ANAME(11),ACON(11) / 'PERIHELION ', PER / DATA ANAME(12),ACON(12) / 'VERNAL EQUINOX ', EQNX / DATA ANAME(13),ACON(13) / 'SUMMER SOLSTICE ', SOLS / DATA ANAME(14),ACON(14) / 'MAX VECT LENGTH ', VECMAX / DATA ANAME(15),ACON(15) / 'MOL WT H2O ', H2OMW / DATA ANAME(16),ACON(16) / 'MOL WT AIR ', AIRMW / DATA ANAME(17),ACON(17) / 'CPV ', CPV / DATA ANAME(18),ACON(18) / 'CPD ', CPD / DATA ANAME(19),ACON(19) / 'UNIV GAS CONST ', RUNIV / DATA ANAME(20),ACON(20) / 'LATENT HEAT SBLM', ALHS / DATA ANAME(21),ACON(21) / 'FREEZING-POINT ', TICE / DATA ANAME(23),ACON(23) / 'CALTOJ ', CALTOJ / DATA ANAME(24),ACON(24) / 'EPS ', EPS / DATA ANAME(25),ACON(25) / 'HEATW ', HEATW / DATA ANAME(26),ACON(26) / 'EPSFAC ', EPSFAC / DATA ANAME(27),ACON(27) / 'VIRTCON ', VIRTCON/ DATA ANAME(28),ACON(28) / 'PIMEAN ', PIMEAN / DATA ANAME(29),ACON(29) / 'SDAY ', SDAY / DATA ANAME(30),ACON(30) / 'HEATI ', HEATI / DATA ANAME(31),ACON(31) / 'S0 ', S0 / DATA ANAME(32),ACON(32) / 'PSTD ', PSTD / DATA ANAME(33),ACON(33) / 'TSTD ', TSTD / DATA ANAME(34),ACON(34) / 'SSALB ', SSALB / DATA ANAME(35),ACON(35) / 'UNDEF ', UNDEF / DATA ANAME(36),ACON(36) / 'CO2 ', CO2 / DATA ANAME(37),ACON(37) / 'RVAP ', RVAP / DATA ANAME(38),ACON(38) / 'CFC11 ', CFC11 / DATA ANAME(39),ACON(39) / 'CFC12 ', CFC12 / DATA ANAME(40),ACON(40) / 'CFC22 ', CFC22 / DO 10 I=1,MAXCON IF(NAME.EQ.ANAME(I)) THEN GETCON = ACON(I) RETURN ENDIF 10 CONTINUE 900 PRINT *,' CANNOT FIND FLOATING POINT CONSTANT - ',NAME PRINT *,' GETCON - CANNOT FIND CONSTANT REQUESTED' end FUNCTION GETCON function nsecf (nhms) C*********************************************************************** C Purpose C Converts NHMS format to Total Seconds C C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C*********************************************************************** implicit none integer nhms, nsecf nsecf = nhms/10000*3600 + mod(nhms,10000)/100*60 + mod(nhms,100) return end function nsecf function nhmsf (nsec) C*********************************************************************** C Purpose C Converts Total Seconds to NHMS format C C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C*********************************************************************** implicit none integer nhmsf, nsec nhmsf = nsec/3600*10000 + mod(nsec,3600)/60*100 + mod(nsec,60) return end function nhmsf subroutine tick (nymd,nhms,ndt) C*********************************************************************** C Purpose C Tick the Date (nymd) and Time (nhms) by NDT (seconds) C C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C*********************************************************************** IF(NDT.NE.0) THEN NSEC = NSECF(NHMS) + NDT IF (NSEC.GT.86400) THEN DO WHILE (NSEC.GT.86400) NSEC = NSEC - 86400 NYMD = INCYMD (NYMD,1) ENDDO ENDIF IF (NSEC.EQ.86400) THEN NSEC = 0 NYMD = INCYMD (NYMD,1) ENDIF IF (NSEC.LT.00000) THEN DO WHILE (NSEC.LT.0) NSEC = 86400 + NSEC NYMD = INCYMD (NYMD,-1) ENDDO ENDIF NHMS = NHMSF (NSEC) ENDIF RETURN end subroutine tick function incymd (NYMD,M) C*********************************************************************** C PURPOSE C INCYMD: NYMD CHANGED BY ONE DAY C MODYMD: NYMD CONVERTED TO JULIAN DATE C DESCRIPTION OF PARAMETERS C NYMD CURRENT DATE IN YYMMDD FORMAT C M +/- 1 (DAY ADJUSTMENT) C C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C*********************************************************************** INTEGER NDPM(12) DATA NDPM /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/ LOGICAL LEAP DATA NY00 / 1900 / LEAP(NY) = MOD(NY,4).EQ.0 .AND. (NY.NE.0 .OR. MOD(NY00,400).EQ.0) C*********************************************************************** C NY = NYMD / 10000 NM = MOD(NYMD,10000) / 100 ND = MOD(NYMD,100) + M IF (ND.EQ.0) THEN NM = NM - 1 IF (NM.EQ.0) THEN NM = 12 NY = NY - 1 ENDIF ND = NDPM(NM) IF (NM.EQ.2 .AND. LEAP(NY)) ND = 29 ENDIF IF (ND.EQ.29 .AND. NM.EQ.2 .AND. LEAP(NY)) GO TO 20 IF (ND.GT.NDPM(NM)) THEN ND = 1 NM = NM + 1 IF (NM.GT.12) THEN NM = 1 NY = NY + 1 ENDIF ENDIF 20 CONTINUE INCYMD = NY*10000 + NM*100 + ND RETURN C*********************************************************************** C E N T R Y M O D Y M D C*********************************************************************** ENTRY MODYMD (NYMD) NY = NYMD / 10000 NM = MOD(NYMD,10000) / 100 ND = MOD(NYMD,100) 40 CONTINUE IF (NM.LE.1) GO TO 60 NM = NM - 1 ND = ND + NDPM(NM) IF (NM.EQ.2 .AND. LEAP(NY)) ND = ND + 1 GO TO 40 60 CONTINUE MODYMD = ND RETURN end function incymd subroutine sigtopl ( qprs,q,logpl,logps,logpt,logp,im,jm,lm,undef ) C*********************************************************************** C C PURPOSE C To interpolate an arbitrary quantity from Model Vertical Grid to Pressure C C INPUT C Q ..... Q (im,jm,lm) Arbitrary Quantity on Model Grid C PKZ ... PKZ (im,jm,lm) Pressure to the Kappa at Model Levels (From Phillips) C PKSRF . PKSRF(im,jm) Surface Pressure to the Kappa C PTOP .. Pressure at Model Top C P ..... Output Pressure Level (mb) C IM .... Longitude Dimension of Input C JM .... Latitude Dimension of Input C LM .... Vertical Dimension of Input C C OUTPUT C QPRS .. QPRS (im,jm) Arbitrary Quantity at Pressure p C C NOTE C Quantity is interpolated Linear in P**Kappa. C Between PTOP**Kappa and PKZ(1), quantity is extrapolated. C Between PKSRF**Kappa and PKZ(LM), quantity is extrapolated. C Undefined Model-Level quantities are not used. C C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C*********************************************************************** C implicit none integer i,j,l,im,jm,lm real qprs(im,jm) real q (im,jm,lm) real logpl(im,jm,lm) real logps(im,jm) real logpt(im,jm) real p,undef,kappa,getcon real logp,logpmin,logpmax,temp kappa = getcon('KAPPA') c Initialize to UNDEFINED c ----------------------- do i=1,im*jm qprs(i,1) = undef enddo c Interpolate to Pressure Between Model Levels c -------------------------------------------- #if (openmp) !$omp parallel do !$omp& default (shared) !$omp& private (i,j,L,logpmin,logpmax,temp) #endif do L=1,lm-1 logpmin = logpl(1,1,L) logpmax = logpl(1,1,L+1) do i=2,im*jm if( logpl(i,1,L) .lt.logpmin ) logpmin = logpl(i,1,L) if( logpl(i,1,L+1).gt.logpmax ) logpmax = logpl(i,1,L+1) enddo if( logp.le.logpmax .and. logp.ge.logpmin ) then do j=1,jm do i=1,im if( logp.le.logpl(i,j,L+1) .and. logp.ge.logpl(i,j,L) ) then temp = ( logpl(i,j,L)-logp ) / ( logpl(i,j,L)-logpl(i,j,L+1) ) if( q(i,j,L) .ne.undef .and. . q(i,j,L+1).ne.undef ) then qprs(i,j) = q(i,j,L+1)*temp + q(i,j,L)*(1.-temp) else if( q(i,j,L+1).ne.undef .and. temp.ge.0.5 ) then qprs(i,j) = q(i,j,L+1) else if( q(i,j,L) .ne.undef .and. temp.le.0.5 ) then qprs(i,j) = q(i,j,L) endif endif enddo enddo endif enddo do i=1,im*jm c Extrapolate to Pressure between Ptop and First Model Level c ---------------------------------------------------------- if( logp.le.logpl(i,1,1) .and. logp.ge.logpt(i,1) ) then temp = ( logpl(i,1,1)-logp ) / ( logpl(i,1,1)-logpl(i,1,2) ) if( q(i,1,1).ne.undef .and. . q(i,1,2).ne.undef ) then qprs(i,1) = q(i,1,2)*temp + q(i,1,1)*(1.-temp) else if( q(i,1,1).ne.undef ) then qprs(i,1) = q(i,1,1) endif endif c Extrapolate to Pressure between Psurf and Lowest Model Level c ------------------------------------------------------------ if( logp.le.logps(i,1) .and. logp.ge.logpl(i,1,lm ) ) then temp = ( logpl(i,1,lm)-logp ) / ( logpl(i,1,lm)-logpl(i,1,lm-1) ) if( q(i,1,lm) .ne.undef .and. . q(i,1,lm-1).ne.undef ) then qprs(i,1) = q(i,1,lm-1)*temp + q(i,1,lm)*(1.-temp) else if( q(i,1,lm) .ne.undef ) then qprs(i,1) = q(i,1,lm) endif endif c Set Values below Ground to Lowest Model Layer c --------------------------------------------- #if (underg) if( logp.gt.logps(i,1) ) qprs(i,1) = q(i,1,lm) #endif enddo return end subroutine sigtopl subroutine read_hdf ( hdffile,im,jm,lm,n2d,n3d,lonbeg,undef,id, . nymd0,nhms0,ndt,ntime, . nvars,names,Lsurf,name2d,titl2d,name3d,titl3d ) implicit none logical, pointer :: Lsurf (:) character*128, pointer :: names (:) character*128, pointer :: name2d(:), name3d(:) character*128, pointer :: titl2d(:), titl3d(:) character*128 hdffile integer id,im,jm,lm,n2d,n3d,nvars,nsecf integer ntime,ngatts,rc,timinc,nymd0,nhms0,ndt real undef real*8 lonbeg integer L,m,n character*128 dummy,name character*128 title character*128 source character*128 contact character*128 levunits character*128, allocatable :: vname(:) character*128, allocatable :: vtitle(:) character*128, allocatable :: vunits(:) real, allocatable :: lat(:) real, allocatable :: lon(:) real, allocatable :: lev(:) real, allocatable :: vrange(:,:) real, allocatable :: prange(:,:) integer, allocatable :: yymmdd(:) integer, allocatable :: hhmmss(:) integer, allocatable :: kmvar(:) integer, allocatable :: loc(:) C ********************************************************************** C **** Read HDF File for Meta Data **** C ********************************************************************** call gfio_open ( trim(hdffile),1,id,rc ) call gfio_diminquire ( id,im,jm,lm,ntime,nvars,ngatts,rc ) allocate ( lon(im) ) allocate ( lat(jm) ) allocate ( lev(lm) ) allocate ( yymmdd(ntime) ) allocate ( hhmmss(ntime) ) allocate ( vname(nvars) ) allocate ( names(nvars) ) allocate ( vtitle(nvars) ) allocate ( vunits(nvars) ) allocate ( kmvar(nvars) ) allocate ( vrange(2,nvars) ) allocate ( prange(2,nvars) ) timinc = 0 call gfio_inquire ( id,im,jm,lm,ntime,nvars, . title,source,contact,undef, . lon,lat,lev,levunits, . yymmdd,hhmmss,timinc, . vname,vtitle,vunits,kmvar, . vrange,prange,rc ) if( timinc == 0 ) timinc = 060000 nymd0 = yymmdd(1) nhms0 = hhmmss(1) ndt = nsecf (timinc) names = vname lonbeg = lon(1) n2d = 0 n3d = 0 do n=1,nvars if( kmvar(n).eq.0 ) then n2d = n2d + 1 else n3d = n3d + 1 endif enddo allocate( Lsurf(nvars) ) allocate( name2d(n2d) ) allocate( titl2d(n2d) ) allocate( name3d(n3d) ) allocate( titl3d(n3d) ) n2d = 0 n3d = 0 do n=1,nvars if( kmvar(n).eq.0 ) then n2d = n2d + 1 name2d(n2d) = vname (n) titl2d(n2d) = vtitle(n) else n3d = n3d + 1 name3d(n3d) = vname (n) titl3d(n3d) = vtitle(n) endif enddo return end subroutine read_hdf subroutine read_ctl ( ctlfile,im,jm,lm,n2d,n3d,lonbeg,undef, . nvars,names,Lsurf,name2d,titl2d,name3d,titl3d ) implicit none logical, pointer :: Lsurf (:) character*128, pointer :: names (:) character*128, pointer :: name2d(:), name3d(:) character*128, pointer :: titl2d(:), titl3d(:) character*128 ctlfile integer im,jm,lm,n2d,n3d,nvars real undef real*8 lonbeg integer L,m,n character*128 dummy,name C ********************************************************************** C **** Read Grads CLT File for Meta Data **** C ********************************************************************** open (10,file=trim(ctlfile),form='formatted') n2d = 0 n3d = 0 do read(10,*,end=500) dummy if( trim(dummy).eq.'xdef' ) then backspace(10) read(10,*) dummy,im,dummy,lonbeg endif if( trim(dummy).eq.'ydef' ) then backspace(10) read(10,*) dummy,jm endif if( trim(dummy).eq.'zdef' ) then backspace(10) read(10,*) dummy,lm endif if( trim(dummy).eq.'undef' ) then backspace(10) read(10,*) dummy,undef endif if( trim(dummy).eq.'vars' ) then backspace(10) read(10,*) dummy,nvars allocate( names(nvars) ) do n=1,nvars read(10,*) names(n),L if( L.eq.0 ) then n2d = n2d + 1 else n3d = n3d + 1 endif enddo endif enddo 500 continue rewind(10) if( n2d.eq.0 .and. n3d.eq.0 ) then print *, 'Warning, n2d = n3d = 0!' stop endif allocate( Lsurf(nvars) ) allocate( name2d(n2d) ) allocate( titl2d(n2d) ) allocate( name3d(n3d) ) allocate( titl3d(n3d) ) n2d = 0 n3d = 0 do read(10,*,end=501) dummy if( trim(dummy).eq.'vars' ) then backspace(10) read(10,*) dummy,nvars do n=1,nvars read(10,*) name,L backspace(10) if( L.eq.0 ) then Lsurf(n) = .true. n2d = n2d + 1 read(10,*) name2d(n2d),L,m,titl2d(n2d) else Lsurf(n) = .false. n3d = n3d + 1 read(10,*) name3d(n3d),L,m,titl3d(n3d) endif enddo endif enddo 501 continue return end subroutine read_ctl subroutine zprime ( a,ap,im,jnp,lm,undef ) real a(im,jnp,lm), ap(im,jnp,lm), az(jnp) DO L=1,LM DO J=1,JNP AZ(J) = 0.0 IC = 0 DO I = 1, IM IF( abs(A(I,J,L)-UNDEF).gt.0.1 ) THEN AZ(J) = AZ(J) + A(I,J,L) IC = IC + 1 ENDIF ENDDO IF(IC.NE.0) AZ(J) = AZ(J) / IC IF(IC.eq.0) AZ(J) = UNDEF ENDDO do j=1,jnp do i=1,im if( abs(a(i,j,L)-undef).gt.0.1 ) then ap(i,j,L) = a(i,j,L)-az(j) else ap(i,j,L) = undef endif enddo enddo ENDDO return end subroutine zprime subroutine writit (q,im,jm,im_out,jm_out,lm,id,name,nymd,nhms,undef,hdf) integer im,jm,lm,L integer id,nymd,nhms,rc,lbeg character*128 name logical hdf real q (im,jm,lm), qout(im_out,jm_out,lm) real*4 qout4(im_out,jm_out) real undef lbeg = 1 if( lm.eq.1 ) lbeg = 0 print *, ' Writing variable: ',trim(name) if( im.eq.im_out .and. jm.eq.jm_out ) then qout = q else call hinterp ( q,im,jm,qout,im_out,jm_out,lm,undef ) endif if(hdf) then call Gfio_putVar ( id,trim(name),nymd,nhms,im_out,jm_out,lbeg,lm,qout,rc ) else do L=lm,1,-1 qout4(:,:) = reshape(qout(:,:,L),(/im_out,jm_out/)) write(55) qout4 enddo endif return end subroutine writit subroutine writit2 (q1,q2,im,jm,im_out,jm_out,lm,id,name,nymd,nhms,undef,hdf) integer im,jm,lm,L integer id,nymd,nhms,rc,lbeg character*128 name logical hdf real q1(im,jm,lm),q2(im,jm,lm),qout(im_out,jm_out,lm) real*4 qout4(im_out,jm_out) real qq(im,jm,lm) real undef do L=1,lm do j=1,jm do i=1,im if( q1(i,j,L).ne.undef .and. q2(i,j,L).ne.undef ) then qq(i,j,L) = q1(i,j,L)*q2(i,j,L) else qq(i,j,L) = undef endif enddo enddo enddo lbeg = 1 if( lm.eq.1 ) lbeg = 0 print *, ' Writing variable: ',trim(name) if( im.eq.im_out .and. jm.eq.jm_out ) then qout = qq else call hinterp ( qq,im,jm,qout,im_out,jm_out,lm,undef ) endif if(hdf) then call Gfio_putVar ( id,trim(name),nymd,nhms,im_out,jm_out,lbeg,lm,qout,rc ) else do L=lm,1,-1 qout4(:,:) = reshape(qout(:,:,L),(/im_out,jm_out/)) write(55) qout4 enddo endif return end subroutine writit2 subroutine intgeop (PHI,THB,PKE,PKO,PHISS,PSTR,UNDEF,im,jm,lm) C*********************************************************************** C C INPUT C THB .... POTENTIAL TEMPERATURE (T/PKO) AT FULL LEVELS C PKE .... P TO THE KAPPA AT THE BASE OF EACH LAYER C PKO .... P TO THE KAPPA AT FULL LEVELS C PHISS .. SURFACE GEOPOTENTIAL C PSTR ... OUTPUT PRESSURE LEVEL (SCALAR) IN MB C UNDEF .. UNDEFINED VALUE C TEM .... SCRATCH ARRAY C C OUTPUT C PHI .... GEOPOTENTIAL HEIGHT (M) AT PSTR C C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C*********************************************************************** PARAMETER (HALF=.5) PARAMETER (ONE=1.) PARAMETER (TWO=2.) real THB(im*jm,lm), PKE(im*jm,lm), PKO(im*jm,lm) real PHI(im*jm) , PHISS(im*jm) , TEM(im*jm,5) G = GETCON('GRAVITY') CP = GETCON('CP') AKP = GETCON('KAPPA') PKSTR = PSTR ** AKP DO 10 I=1,im*jm PHI(I) = UNDEF 10 CONTINUE PKMIN = PKO(1,lm) DO 20 I=2,im*jm IF( PKO(I,lm).LT.PKMIN ) PKMIN = PKO(I,lm) 20 CONTINUE IF( PKSTR.GE.PKMIN ) THEN DO 30 I=1,im*jm IF (PKSTR .LE. PKE(I,lm) .AND. PKSTR .GE. PKO(I,lm)) THEN PHI(I) = PHISS(I) + (PKE(I,lm)-PKSTR) * THB(I,lm)*CP ENDIF 30 CONTINUE ENDIF DO 40 I=1,im*jm TEM(I,5) = PHISS(I) . + (PKE(I,lm)-PKO(I,lm)) * THB(I,lm)*CP 40 CONTINUE DO 310 L=lm,2,-1 PKMAX = PKO(1,L) DO 50 I=2,im*jm IF( PKO(I,L).GT.PKMAX ) PKMAX = PKO(I,L) 50 CONTINUE IF( PKSTR.GE.PKMAX ) GO TO 320 DO 60 I=1,im*jm TEM(I,1) = THB(I,L-1)*(PKE(I,L-1)-PKO(I,L-1)) . + THB(I,L )*(PKO(I,L )-PKE(I,L-1)) IF( PKSTR.GE.PKE(I,L-1) ) THEN TEM(I,3) = PKO(I,L) TEM(I,2) = TEM(I,5) TEM(I,4) = THB(I,L) ENDIF TEM(I,5) = TEM(I,5) + TEM(I,1) * CP IF( PKSTR.LT.PKE(I,L-1) ) THEN TEM(I,3) = PKO(I,L-1) TEM(I,2) = TEM(I,5 ) TEM(I,4) = THB(I,L-1) ENDIF IF( PKSTR.GE.PKO(I,L-1) .AND. PKSTR.LT.PKO(I,L) ) THEN PHI(I) = TEM(I,2) ENDIF IF(L.EQ.2) THEN IF( PKSTR.LT.PKO(I,1) ) PHI(I) = TEM(I,2) ENDIF TEM(I,1) = TEM(I,1) / (PKO(I,L)-PKO(I,L-1)) TEM(I,2) = (TEM(I,3)-PKSTR) / (TEM(I,3)-PKE(I,L-1)) TEM(I,1) = TEM(I,1) * TEM(I,2) . + TEM(I,4) * (TWO-TEM(I,2)) TEM(I,1) = TEM(I,1) * (TEM(I,3)-PKSTR) IF( PKSTR.GE.PKO(I,L-1) .AND. PKSTR.LT.PKO(I,L)) THEN PHI(I) = PHI(I) + TEM(I,1) * (CP*HALF) ENDIF IF(L.EQ.2) THEN IF( PKSTR.LT.PKO(I,1) ) PHI(I) = PHI(I) + TEM(I,1)*(CP*HALF) ENDIF 60 CONTINUE 310 CONTINUE 320 CONTINUE DO 70 I=1,im*jm IF( PHI(I).NE.UNDEF ) PHI(I) = PHI(I) * (ONE/G) 70 CONTINUE RETURN end subroutine intgeop subroutine hinterp ( qin,iin,jin,qout,iout,jout,mlev,undef ) implicit none integer iin,jin, iout,jout, mlev real qin(iin,jin,mlev), qout(iout,jout,mlev) real undef,pi,dlin,dpin,dlout,dpout real dlam(iin), lons(iout*jout), lon real dphi(jin), lats(iout*jout), lat integer i,j,loc pi = 4.0*atan(1.0) dlin = 2*pi/iin dpin = pi/(jin-1) dlam(:) = dlin dphi(:) = dpin dlout = 2*pi/iout dpout = pi/(jout-1) loc = 0 do j=1,jout do i=1,iout loc = loc + 1 lon = -pi + (i-1)*dlout lons(loc) = lon enddo enddo loc = 0 do j=1,jout lat = -pi/2.0 + (j-1)*dpout do i=1,iout loc = loc + 1 lats(loc) = lat enddo enddo call interp_h ( qin,iin,jin,mlev,dlam,dphi, . qout,iout*jout,lons,lats,undef ) return end subroutine hinterp subroutine interp_h ( q_cmp,im,jm,lm,dlam,dphi, . q_geo,irun,lon_geo,lat_geo,undef ) C*********************************************************************** C C PURPOSE: C ======== C Performs a horizontal interpolation from a field on a computational grid C to arbitrary locations. C C INPUT: C ====== C q_cmp ...... Field q_cmp(im,jm,lm) on the computational grid C im ......... Longitudinal dimension of q_cmp C jm ......... Latitudinal dimension of q_cmp C lm ......... Vertical dimension of q_cmp C dlam ....... Computational Grid Delta Lambda C dphi ....... Computational Grid Delta Phi C irun ....... Number of Output Locations C lon_geo .... Longitude Location of Output C lat_geo .... Latitude Location of Output C C OUTPUT: C ======= C q_geo ...... Field q_geo(irun,lm) at arbitrary locations C C C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C*********************************************************************** implicit none c Input Variables c --------------- integer im,jm,lm,irun real q_geo(irun,lm) real lon_geo(irun) real lat_geo(irun) real q_cmp(im,jm,lm) real dlam(im) real dphi(jm) c Local Variables c --------------- integer i,j,l,m,n integer, allocatable :: ip1(:), ip0(:), im1(:), im2(:) integer, allocatable :: jp1(:), jp0(:), jm1(:), jm2(:) integer ip1_for_jp1, ip0_for_jp1, im1_for_jp1, im2_for_jp1 integer ip1_for_jm2, ip0_for_jm2, im1_for_jm2, im2_for_jm2 integer jm2_for_jm2, jp1_for_jp1 c Bi-Linear Weights c ----------------- real, allocatable :: wl_ip0jp0 (:) real, allocatable :: wl_im1jp0 (:) real, allocatable :: wl_ip0jm1 (:) real, allocatable :: wl_im1jm1 (:) c Bi-Cubic Weights c ---------------- real, allocatable :: wc_ip1jp1 (:) real, allocatable :: wc_ip0jp1 (:) real, allocatable :: wc_im1jp1 (:) real, allocatable :: wc_im2jp1 (:) real, allocatable :: wc_ip1jp0 (:) real, allocatable :: wc_ip0jp0 (:) real, allocatable :: wc_im1jp0 (:) real, allocatable :: wc_im2jp0 (:) real, allocatable :: wc_ip1jm1 (:) real, allocatable :: wc_ip0jm1 (:) real, allocatable :: wc_im1jm1 (:) real, allocatable :: wc_im2jm1 (:) real, allocatable :: wc_ip1jm2 (:) real, allocatable :: wc_ip0jm2 (:) real, allocatable :: wc_im1jm2 (:) real, allocatable :: wc_im2jm2 (:) real ux, ap1, ap0, am1, am2 real uy, bp1, bp0, bm1, bm2 real lon_cmp(im) real lat_cmp(jm) real q_tmp(irun) real pi,d real lam,lam_ip1,lam_ip0,lam_im1,lam_im2 real phi,phi_jp1,phi_jp0,phi_jm1,phi_jm2 real dl,dp,phi_np,lam_0 real lam_geo, lam_cmp real phi_geo, phi_cmp real undef integer im1_cmp,icmp integer jm1_cmp,jcmp c Initialization c -------------- pi = 4.*atan(1.) dl = 2*pi/ im ! Uniform Grid Delta Lambda dp = pi/(jm-1) ! Uniform Grid Delta Phi c Allocate Memory for Weights and Index Locations c ----------------------------------------------- allocate ( wl_ip0jp0(irun) , wl_im1jp0(irun) ) allocate ( wl_ip0jm1(irun) , wl_im1jm1(irun) ) allocate ( wc_ip1jp1(irun) , wc_ip0jp1(irun) , wc_im1jp1(irun) , wc_im2jp1(irun) ) allocate ( wc_ip1jp0(irun) , wc_ip0jp0(irun) , wc_im1jp0(irun) , wc_im2jp0(irun) ) allocate ( wc_ip1jm1(irun) , wc_ip0jm1(irun) , wc_im1jm1(irun) , wc_im2jm1(irun) ) allocate ( wc_ip1jm2(irun) , wc_ip0jm2(irun) , wc_im1jm2(irun) , wc_im2jm2(irun) ) allocate ( ip1(irun) , ip0(irun) , im1(irun) , im2(irun) ) allocate ( jp1(irun) , jp0(irun) , jm1(irun) , jm2(irun) ) c Compute Input Computational-Grid Latitude and Longitude Locations c ----------------------------------------------------------------- lon_cmp(1) = -pi do i=2,im lon_cmp(i) = lon_cmp(i-1) + dlam(i-1) enddo lat_cmp(1) = -pi*0.5 do j=2,jm-1 lat_cmp(j) = lat_cmp(j-1) + dphi(j-1) enddo lat_cmp(jm) = pi*0.5 c Compute Weights for Computational to Geophysical Grid Interpolation c ------------------------------------------------------------------- do i=1,irun lam_cmp = lon_geo(i) phi_cmp = lat_geo(i) c Determine Indexing Based on Computational Grid c ---------------------------------------------- im1_cmp = 1 do icmp = 2,im if( lon_cmp(icmp).lt.lam_cmp ) im1_cmp = icmp enddo jm1_cmp = 1 do jcmp = 2,jm if( lat_cmp(jcmp).lt.phi_cmp ) jm1_cmp = jcmp enddo im1(i) = im1_cmp ip0(i) = im1(i) + 1 ip1(i) = ip0(i) + 1 im2(i) = im1(i) - 1 jm1(i) = jm1_cmp jp0(i) = jm1(i) + 1 jp1(i) = jp0(i) + 1 jm2(i) = jm1(i) - 1 c Fix Longitude Index Boundaries c ------------------------------ if(im1(i).eq.im) then ip0(i) = 1 ip1(i) = 2 endif if(im1(i).eq.1) then im2(i) = im endif if(ip0(i).eq.im) then ip1(i) = 1 endif c Compute Immediate Surrounding Coordinates c ----------------------------------------- lam = lam_cmp phi = phi_cmp c Compute and Adjust Longitude Weights c ------------------------------------ lam_im2 = lon_cmp(im2(i)) lam_im1 = lon_cmp(im1(i)) lam_ip0 = lon_cmp(ip0(i)) lam_ip1 = lon_cmp(ip1(i)) if( lam_im2.gt.lam_im1 ) lam_im2 = lam_im2 - 2*pi if( lam_im1.gt.lam_ip0 ) lam_ip0 = lam_ip0 + 2*pi if( lam_im1.gt.lam_ip1 ) lam_ip1 = lam_ip1 + 2*pi if( lam_ip0.gt.lam_ip1 ) lam_ip1 = lam_ip1 + 2*pi c Compute and Adjust Latitude Weights c Note: Latitude Index Boundaries are Adjusted during Interpolation c ------------------------------------------------------------------ phi_jm2 = lat_cmp(jm2(i)) phi_jm1 = lat_cmp(jm1(i)) phi_jp0 = lat_cmp(jp0(i)) phi_jp1 = lat_cmp(jp1(i)) if( jm2(i).eq.0 ) phi_jm2 = phi_jm1 - dphi(1) if( jm1(i).eq.jm ) then phi_jp0 = phi_jm1 + dphi(jm-1) phi_jp1 = phi_jp0 + dphi(jm-2) endif if( jp1(i).eq.jm+1 ) phi_jp1 = phi_jp0 + dphi(jm-1) c Bi-Linear Weights c ----------------- d = (lam_ip0-lam_im1)*(phi_jp0-phi_jm1) wl_im1jm1(i) = (lam_ip0-lam )*(phi_jp0-phi )/d wl_ip0jm1(i) = (lam -lam_im1)*(phi_jp0-phi )/d wl_im1jp0(i) = (lam_ip0-lam )*(phi -phi_jm1)/d wl_ip0jp0(i) = (lam -lam_im1)*(phi -phi_jm1)/d c Bi-Cubic Weights c ---------------- ap1 = ( (lam -lam_ip0)*(lam -lam_im1)*(lam -lam_im2) ) . / ( (lam_ip1-lam_ip0)*(lam_ip1-lam_im1)*(lam_ip1-lam_im2) ) ap0 = ( (lam_ip1-lam )*(lam -lam_im1)*(lam -lam_im2) ) . / ( (lam_ip1-lam_ip0)*(lam_ip0-lam_im1)*(lam_ip0-lam_im2) ) am1 = ( (lam_ip1-lam )*(lam_ip0-lam )*(lam -lam_im2) ) . / ( (lam_ip1-lam_im1)*(lam_ip0-lam_im1)*(lam_im1-lam_im2) ) am2 = ( (lam_ip1-lam )*(lam_ip0-lam )*(lam_im1-lam ) ) . / ( (lam_ip1-lam_im2)*(lam_ip0-lam_im2)*(lam_im1-lam_im2) ) bp1 = ( (phi -phi_jp0)*(phi -phi_jm1)*(phi -phi_jm2) ) . / ( (phi_jp1-phi_jp0)*(phi_jp1-phi_jm1)*(phi_jp1-phi_jm2) ) bp0 = ( (phi_jp1-phi )*(phi -phi_jm1)*(phi -phi_jm2) ) . / ( (phi_jp1-phi_jp0)*(phi_jp0-phi_jm1)*(phi_jp0-phi_jm2) ) bm1 = ( (phi_jp1-phi )*(phi_jp0-phi )*(phi -phi_jm2) ) . / ( (phi_jp1-phi_jm1)*(phi_jp0-phi_jm1)*(phi_jm1-phi_jm2) ) bm2 = ( (phi_jp1-phi )*(phi_jp0-phi )*(phi_jm1-phi ) ) . / ( (phi_jp1-phi_jm2)*(phi_jp0-phi_jm2)*(phi_jm1-phi_jm2) ) wc_ip1jp1(i) = bp1*ap1 wc_ip0jp1(i) = bp1*ap0 wc_im1jp1(i) = bp1*am1 wc_im2jp1(i) = bp1*am2 wc_ip1jp0(i) = bp0*ap1 wc_ip0jp0(i) = bp0*ap0 wc_im1jp0(i) = bp0*am1 wc_im2jp0(i) = bp0*am2 wc_ip1jm1(i) = bm1*ap1 wc_ip0jm1(i) = bm1*ap0 wc_im1jm1(i) = bm1*am1 wc_im2jm1(i) = bm1*am2 wc_ip1jm2(i) = bm2*ap1 wc_ip0jm2(i) = bm2*ap0 wc_im1jm2(i) = bm2*am1 wc_im2jm2(i) = bm2*am2 enddo c Interpolate Computational-Grid Quantities to Geophysical Grid c ------------------------------------------------------------- do L=1,lm do i=1,irun if( lat_geo(i).le.lat_cmp(2) .or. . lat_geo(i).ge.lat_cmp(jm-1) ) then c 1st Order Interpolation at Poles c -------------------------------- if( q_cmp( im1(i),jm1(i),L ).ne.undef .and. . q_cmp( ip0(i),jm1(i),L ).ne.undef .and. . q_cmp( im1(i),jp0(i),L ).ne.undef .and. . q_cmp( ip0(i),jp0(i),L ).ne.undef ) then q_tmp(i) = wl_im1jm1(i) * q_cmp( im1(i),jm1(i),L ) . + wl_ip0jm1(i) * q_cmp( ip0(i),jm1(i),L ) . + wl_im1jp0(i) * q_cmp( im1(i),jp0(i),L ) . + wl_ip0jp0(i) * q_cmp( ip0(i),jp0(i),L ) else q_tmp(i) = undef endif else c Cubic Interpolation away from Poles c ----------------------------------- if( q_cmp( ip1(i),jp0(i),L ).ne.undef .and. . q_cmp( ip0(i),jp0(i),L ).ne.undef .and. . q_cmp( im1(i),jp0(i),L ).ne.undef .and. . q_cmp( im2(i),jp0(i),L ).ne.undef .and. . q_cmp( ip1(i),jm1(i),L ).ne.undef .and. . q_cmp( ip0(i),jm1(i),L ).ne.undef .and. . q_cmp( im1(i),jm1(i),L ).ne.undef .and. . q_cmp( im2(i),jm1(i),L ).ne.undef .and. . q_cmp( ip1(i),jp1(i),L ).ne.undef .and. . q_cmp( ip0(i),jp1(i),L ).ne.undef .and. . q_cmp( im1(i),jp1(i),L ).ne.undef .and. . q_cmp( im2(i),jp1(i),L ).ne.undef .and. . q_cmp( ip1(i),jm2(i),L ).ne.undef .and. . q_cmp( ip0(i),jm2(i),L ).ne.undef .and. . q_cmp( im1(i),jm2(i),L ).ne.undef .and. . q_cmp( im2(i),jm2(i),L ).ne.undef ) then q_tmp(i) = wc_ip1jp1(i) * q_cmp( ip1(i),jp1(i),L ) . + wc_ip0jp1(i) * q_cmp( ip0(i),jp1(i),L ) . + wc_im1jp1(i) * q_cmp( im1(i),jp1(i),L ) . + wc_im2jp1(i) * q_cmp( im2(i),jp1(i),L ) . + wc_ip1jp0(i) * q_cmp( ip1(i),jp0(i),L ) . + wc_ip0jp0(i) * q_cmp( ip0(i),jp0(i),L ) . + wc_im1jp0(i) * q_cmp( im1(i),jp0(i),L ) . + wc_im2jp0(i) * q_cmp( im2(i),jp0(i),L ) . + wc_ip1jm1(i) * q_cmp( ip1(i),jm1(i),L ) . + wc_ip0jm1(i) * q_cmp( ip0(i),jm1(i),L ) . + wc_im1jm1(i) * q_cmp( im1(i),jm1(i),L ) . + wc_im2jm1(i) * q_cmp( im2(i),jm1(i),L ) . + wc_ip1jm2(i) * q_cmp( ip1(i),jm2(i),L ) . + wc_ip0jm2(i) * q_cmp( ip0(i),jm2(i),L ) . + wc_im1jm2(i) * q_cmp( im1(i),jm2(i),L ) . + wc_im2jm2(i) * q_cmp( im2(i),jm2(i),L ) elseif( q_cmp( im1(i),jm1(i),L ).ne.undef .and. . q_cmp( ip0(i),jm1(i),L ).ne.undef .and. . q_cmp( im1(i),jp0(i),L ).ne.undef .and. . q_cmp( ip0(i),jp0(i),L ).ne.undef ) then q_tmp(i) = wl_im1jm1(i) * q_cmp( im1(i),jm1(i),L ) . + wl_ip0jm1(i) * q_cmp( ip0(i),jm1(i),L ) . + wl_im1jp0(i) * q_cmp( im1(i),jp0(i),L ) . + wl_ip0jp0(i) * q_cmp( ip0(i),jp0(i),L ) else q_tmp(i) = undef endif endif enddo c Load Temp array into Output array c --------------------------------- do i=1,irun q_geo(i,L) = q_tmp(i) enddo enddo deallocate ( wl_ip0jp0 , wl_im1jp0 ) deallocate ( wl_ip0jm1 , wl_im1jm1 ) deallocate ( wc_ip1jp1 , wc_ip0jp1 , wc_im1jp1 , wc_im2jp1 ) deallocate ( wc_ip1jp0 , wc_ip0jp0 , wc_im1jp0 , wc_im2jp0 ) deallocate ( wc_ip1jm1 , wc_ip0jm1 , wc_im1jm1 , wc_im2jm1 ) deallocate ( wc_ip1jm2 , wc_ip0jm2 , wc_im1jm2 , wc_im2jm2 ) deallocate ( ip1 , ip0 , im1 , im2 ) deallocate ( jp1 , jp0 , jm1 , jm2 ) return end subroutine interp_h subroutine usage() print *, "Usage: " print * print *, " eta2prs_$ARCH.x -levs plev(s)" print *, " -eta eta_fname(s)" print * print *, " [-nymd nymd]" print *, " [-nhms nhms]" print *, " [-ndt ndt]" print *, " [-im im_out]" print *, " [-jm jm_out]" print *, " [-hdf Logical HDF Flag]" print *, " [-tag output_tag] " print *, " [-ctl ctl_fname]" print * print *, "where:" print * print *, " -levs plev(s): Pressure Level(s) for output file (bottom to top, i.e: 1000 850 700 etc)" print *, " -eta eta_fname(s): Filename(s) in eta format (HDF or Flat)" print * print *, "Optional Args:" print * print *, " -nymd nymd: Beginning YYYYMMDD (Required for Flat Files)" print *, " -nhms nhms: Beginning HHMMSS (Required for Flat Files)" print *, " -ndt ndt: Time Increment (secs)(Required for Flat Files)" print *, " -tag output_tag: Optional Filename Tag for output file (default: eta2prs.nc4)" print *, " -o outfilename: Output filename (-tag opt ignored; default eta2prs.nc4)" print *, " -im im_out: Output Resolution in X (default: Input Resolution)" print *, " -jm jm_out: Output Resolution in Y (default: Input Resolution)" print *, " -ctl ctl_fname: CTL Filename for Flat Files" print * print *, " -quad : Implies Compute Quadratics (default: no quadratics)" print *, " -flat : Implies Flat Binary Output (default: hdf output)" print * call exit(7) end subroutine usage