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 ! ********************************************************************** ! ********************************************************************** ! **** **** ! **** Program to create an eta_hdf file from GEOS5 restarts **** ! **** **** ! ********************************************************************** ! ********************************************************************** character*256 dynrst, moistrst, topofile, bkgsfc, bkgeta,tag integer headr1(6) integer headr2(5) integer nymd,nhms integer im,jm,lm real*8 undef,rgas,rvap,eps,grav ! restart variables and topography ! -------------------------------- real*8, allocatable :: dp(:,:,:) real*8, allocatable :: u(:,:,:) real*8, allocatable :: v(:,:,:) real*8, allocatable :: th(:,:,:) real*8, allocatable :: thv(:,:,:) real*8, allocatable :: pk(:,:,:) real*8, allocatable :: ple(:,:,:) real*8, allocatable :: q(:,:,:) real*8, allocatable :: ps(:,:) real*8, allocatable :: ak(:) real*8, allocatable :: bk(:) real*8, allocatable :: phis (:,:) real*8, allocatable :: topo_std (:,:) real*8, allocatable :: oro (:,:) real*8, allocatable :: frland (:,:) real*8, allocatable :: frlandice(:,:) real*8, allocatable :: frlake (:,:) real*8, allocatable :: frocean (:,:) real*8, allocatable :: u10m (:,:) real*8, allocatable :: v10m (:,:) real*8, allocatable :: ts (:,:) real*8, allocatable :: qs (:,:) real*8, allocatable :: snowmass (:,:) real*8, allocatable :: wet1 (:,:) real*8, allocatable :: tsoil1 (:,:) real*8, allocatable :: oz (:,:,:) real*4, allocatable :: dum(:,:) character*256, allocatable :: arg(:) character*5 extn character*8 date character*2 hour integer n,nargs,iargc,i,j,L,nymd0,nhms0 integer precision logical file_exists ! ********************************************************************** ! **** Initialize Filenames **** ! ********************************************************************** undef = 1.0e15 tag = 'geos5' dynrst = 'xxx' moistrst = 'xxx' topofile = 'xxx' nymd0 = -999 nhms0 = -999 precision = 0 ! 32-bit: 0, 64-bit: 1 extn = 'nc4' 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.'-h' ) call usage() if( trim(arg(n)).eq.'-help' ) call usage() if( trim(arg(n)).eq.'-H' ) call usage() if( trim(arg(n)).eq.'-Help' ) call usage() if( trim(arg(n)).eq.'-tag' ) then tag = trim(arg(n+1)) endif if( trim(arg(n)).eq.'-dynrst' ) then dynrst = trim(arg(n+1)) endif if( trim(arg(n)).eq.'-topo' ) then topofile = trim(arg(n+1)) endif if( trim(arg(n)).eq.'-moistrst' ) then moistrst = trim(arg(n+1)) endif 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.'-precision' ) read(arg(n+1),*) precision if( trim(arg(n)).eq.'-ext' ) read(arg(n+1),*) extn enddo endif if( trim(dynrst) .eq.'xxx' .or. . trim(moistrst).eq.'xxx' ) then print * call usage() endif print * print *, ' tag: ',trim(tag) print *, ' dyn restart filename: ',trim(dynrst) print *, 'moist restart filename: ',trim(moistrst) if( trim(topofile).ne.'xxx' ) then print *, ' topo filename: ',trim(topofile) endif ! ********************************************************************** ! **** Read dycore internal Restart **** ! ********************************************************************** inquire ( file=trim(dynrst), exist=file_exists ) if( file_exists ) then open (10,file=trim(dynrst),form='unformatted',access='sequential') read (10) headr1 read (10) headr2 if( nymd0.ne.-999 ) then nymd = nymd0 else nymd = headr1(1)*10000 + headr1(2)*100 + headr1(3) endif if( nhms0.ne.-999 ) then nhms = nhms0 else nhms = headr1(4)*10000 + headr1(5)*100 + headr1(6) endif im = headr2(1) jm = headr2(2) lm = headr2(3) print *, ' resolution: ',im,jm,lm print *, ' date: ',nymd,nhms if( precision.eq.0 ) then print *, ' precision: 0 (32-bit)' endif if( precision.eq.1 ) then print *, ' precision: 1 (64-bit)' endif print * allocate ( u(im,jm,lm) ) allocate ( v(im,jm,lm) ) allocate ( th(im,jm,lm) ) allocate ( thv(im,jm,lm) ) allocate ( dp(im,jm,lm) ) allocate ( pk(im,jm,lm) ) allocate ( ple(im,jm,lm+1) ) allocate ( ps(im,jm) ) allocate ( ak(lm+1) ) allocate ( bk(lm+1) ) read (10) ak read (10) bk do L=1,lm read(10) (( u(i,j,L),i=1,im),j=1,jm) enddo do L=1,lm read(10) (( v(i,j,L),i=1,im),j=1,jm) enddo do L=1,lm read(10) (( th(i,j,L),i=1,im),j=1,jm) ! Note: GEOS-5 variable is DRY potential temperature enddo do L=1,lm+1 read(10) ((ple(i,j,L),i=1,im),j=1,jm) enddo do L=1,lm read(10) (( pk(i,j,L),i=1,im),j=1,jm) enddo close (10) ps(:,:) = ple(:,:,lm+1) do L=lm,1,-1 dp(:,:,L) = ple(:,:,L+1) - ple(:,:,L) enddo else print *, trim(dynrst),' does not exist' stop endif print * write(6,500) 500 format(2x,' k ',' A(k) ',2x,' B(k) ',2x,' Pref ',2x,' DelP',/, . 1x,'----',3x,'----------',2x,'--------',2x,'----------',2x,'---------' ) L=1 write(6,501) L,ak(L)*0.01, bk(L), ak(L)*0.01 + 1000.0*bk(L) do L=2,lm+1 write(6,502) L,ak(L)*0.01, bk(L), ak(L)*0.01 + 1000.0*bk(L), . (ak(L)-ak(L-1))*0.01 + 1000.0*(bk(L)-bk(L-1)) enddo print * 501 format(2x,i3,2x,f10.6,2x,f8.4,2x,f10.4) 502 format(2x,i3,2x,f10.6,2x,f8.4,2x,f10.4,3x,f8.4) ! ********************************************************************** ! **** Read moist internal Restart **** ! ********************************************************************** allocate ( q(im,jm,lm) ) allocate ( dum(im,jm) ) inquire ( file=trim(moistrst), exist=file_exists ) if( file_exists ) then open (10,file=trim(moistrst),form='unformatted',access='sequential') do L=1,lm read(10) dum q(:,:,L) = dum(:,:) ! First moist variable is SPHU enddo close (10) endif ! ********************************************************************** ! **** Read topo file **** ! ********************************************************************** allocate ( phis (im,jm) ) allocate ( topo_std (im,jm) ) allocate ( oro (im,jm) ) allocate ( frland (im,jm) ) allocate ( frlandice(im,jm) ) allocate ( frlake (im,jm) ) allocate ( frocean (im,jm) ) allocate ( u10m (im,jm) ) allocate ( v10m (im,jm) ) allocate ( ts (im,jm) ) allocate ( qs (im,jm) ) allocate ( snowmass (im,jm) ) allocate ( wet1 (im,jm) ) allocate ( tsoil1 (im,jm) ) allocate ( oz (im,jm,lm) ) phis = 0.0 inquire ( file=trim(topofile), exist=file_exists ) if( file_exists ) then open (10,file=trim(topofile),form='unformatted',access='sequential') print *, 'Reading topography from: ',trim(topofile) read(10) dum ; phis = dum close (10) endif ! ********************************************************************** ! **** Adjust Fields **** ! ********************************************************************** grav = 9.8 rgas = 8314.3/28.97 rvap = 8314.3/18.01 eps = rvap/rgas-1.0 thv = th*(1+eps*q) ! Construct virtual potential temperature phis = phis*grav topo_std = 0.0 ! Analysis needs to read it, but doesn't use it snowmass = snowmass/1000 ! Convert millimeter to meter oz = oz/1.655e-6 oro = 1.0 ! Land where ( frocean >= 0.6 .or. frlake >= 0.6 ) oro = 0.0 ! Water where ( oro == 0.0 .and. ts <= 272.81 ) oro = 2.0 ! Ice do j=1,jm do i=1,im wet1(i,j) = ( frland(i,j)*wet1(i,j) + 1.0-frland(i,j) ) if( oro(i,j) == 1.0 ) then tsoil1(i,j) = tsoil1(i,j) + 273.16 else tsoil1(i,j) = undef endif enddo enddo where ( frlandice == 1.0 ) snowmass = 4.0 ! ********************************************************************** ! **** Write Analysis BKG Files **** ! ********************************************************************** write(date,101) nymd write(hour,102) nhms/10000 101 format(i8.8) 102 format(i2.2) bkgeta = trim(tag) // '.bkg.eta.' // date // '_' // hour // 'z.'//trim(extn) print *, 'Creating ',trim(tag),' bkg.eta HDF file: ',trim(bkgeta) print * call hflip ( phis,im,jm, 1 ) call hflip ( ps ,im,jm, 1 ) call hflip ( dp ,im,jm,lm ) call hflip ( u ,im,jm,lm ) call hflip ( v ,im,jm,lm ) call hflip ( thv ,im,jm,lm ) call hflip ( q ,im,jm,lm ) call hflip ( oz ,im,jm,lm ) call hflip ( u10m ,im,jm, 1 ) call hflip ( v10m ,im,jm, 1 ) call hflip ( ts ,im,jm, 1 ) call hflip ( qs ,im,jm, 1 ) call hflip ( snowmass ,im,jm, 1 ) call hflip ( wet1 ,im,jm, 1 ) call hflip ( tsoil1 ,im,jm, 1 ) call hflip ( oro ,im,jm, 1 ) call hflip ( frland ,im,jm, 1 ) call hflip ( frlandice,im,jm, 1 ) call hflip ( frlake ,im,jm, 1 ) call hflip ( frocean ,im,jm, 1 ) call dynhdf ( bkgeta,nymd,nhms,im,jm,lm,ak,bk, . phis,topo_std,ts,oro, . ps,dp,u,v,thv,q,oz,precision ) stop end subroutine minmax (q,im,jm,L) real*8 q(im,jm) qmin = q(1,1) qmax = q(1,1) do j=1,jm do i=1,im qmin = min( qmin,q(i,j) ) qmax = max( qmax,q(i,j) ) enddo enddo print *, 'L: ',L,' qmin: ',qmin,' qmax: ',qmax return end subroutine dynhdf ( filename,nymd,nhms,im,jm,lm,ak,bk, . phis,topo_std,ts,oro,ps,dp,u,v,thv,q,oz,precision ) implicit none integer im,jm,lm,nymd,nhms real*8 sump(im,jm) real*8 phis(im,jm) real*8 slp(im,jm) real*8 ps(im,jm) real*8 ts(im,jm) real*8 oro(im,jm) real*8 topo_std(im,jm) real*8 dp(im,jm,lm) real*8 u(im,jm,lm) real*8 v(im,jm,lm) real*8 t(im,jm,lm) real*8 thv(im,jm,lm) real*8 q(im,jm,lm) real*8 oz(im,jm,lm) real*8 pk(im,jm,lm) real*8 pke(im,jm,lm+1) real*8 ple(im,jm,lm+1) real*8 phi(im,jm,lm+1) real*8 ak(lm+1) real*8 bk(lm+1) real lats(jm) real lons(im) real levs(lm) real*8 latsd(jm) real*8 lonsd(im) real ptop,dlon,dlat,pref,undef,pint integer i,j,L,n,timeinc,rc,ks character*80 filename character*3 ch3 integer nvars,fid,precision character*256 levunits character*256 title character*256 source character*256 contact character*256, allocatable :: vname(:) character*256, allocatable :: vtitle(:) character*256, allocatable :: vunits(:) integer, allocatable :: lmvar(:) real, allocatable :: v_range(:,:) real, allocatable :: p_range(:,:) real, allocatable :: dum1(:) real, allocatable :: dum2(:,:) real, allocatable :: dum3(:,:,:) real rgas,rvap,eps,kappa,grav,cp real dpref dpref(L) = ( ak(L+1)-ak(L) ) + ( bk(L+1)-bk(L) ) * 98400.0 undef = 1.0e15 timeinc = 060000 cp = 1004.16 rgas = 8314.3/28.97 rvap = 8314.3/18.01 eps = rvap/rgas-1.0 kappa = 2.0/7.0 grav = 9.81 ! String and vars settings ! ------------------------ title = 'GEOS5 Dynamics State Vector (Hybrid Coordinates)' source = 'Goddard Modeling and Assimilation Office, NASA/GSFC' contact = 'data@gmao.gsfc.nasa.gov' levunits = 'hPa' nvars = 14 allocate ( vname(nvars) ) allocate ( vtitle(nvars) ) allocate ( vunits(nvars) ) allocate ( lmvar(nvars) ) allocate ( v_range(2,nvars) ) allocate ( p_range(2,nvars) ) vname(01) = 'phis' vtitle(01) = 'Topography geopotential' vunits(01) = 'meter2/sec2' lmvar(01) = 0 vname(02) = 'hs_stdv' vtitle(02) = 'Topography Height Standard Deviation' vunits(02) = 'm' lmvar(02) = 0 vname(03) = 'ts' vtitle(03) = 'Surface temperature' vunits(03) = 'K' lmvar(03) = 0 vname(04) = 'lwi' vtitle(04) = 'Land-water-ice mask' vunits(04) = 'na' lmvar(04) = 0 vname(05) = 'ps' vtitle(05) = 'Surface Pressure' vunits(05) = 'Pa' lmvar(05) = 0 vname(06) = 'delp' vtitle(06) = 'Pressure Thickness' vunits(06) = 'Pa' lmvar(06) = lm vname(07) = 'uwnd' vtitle(07) = 'Zonal Wind' vunits(07) = 'm/s' lmvar(07) = lm vname(08) = 'vwnd' vtitle(08) = 'Meridional Wind' vunits(08) = 'm/s' lmvar(08) = lm vname(09) = 'theta' vtitle(09) = 'Scaled Virtual Potential Temperature' vunits(09) = 'K/Pa^kappa' lmvar(09) = lm vname(10) = 'sphu' vtitle(10) = 'Specific Humidity' vunits(10) = 'kg/kg' lmvar(10) = lm vname(11) = 'ozone' vtitle(11) = 'Ozone Mixing Ratio' vunits(11) = 'kg/kg' lmvar(11) = lm vname(12) = 'slp' vtitle(12) = 'Sea Level Pressure' vunits(12) = 'Pa' lmvar(12) = 0 vname(13) = 'tmpu' vtitle(13) = 'Temperature' vunits(13) = 'K' lmvar(13) = lm vname(14) = 'zle' vtitle(14) = 'Edge Heights' vunits(14) = 'm' lmvar(14) = lm v_range(:,:) = undef p_range(:,:) = undef ! Compute grid ! ------------ dlon = 360.0/ im dlat = 180.0/(jm-1) do j=1,jm latsd(j) = -90.0 + (j-1)*dlat enddo do i=1,im lonsd(i) = 0.0 + (i-1)*dlon enddo do L=1,lm levs(L) = L enddo lats = latsd lons = lonsd sump(:,:) = 0.0 do L=1,lm sump(:,:) = sump(:,:) + dp(:,:,L) enddo sump(:,:) = ps(:,:)-sump(:,:) ! print *, 'ptop values (mb): ',sump/100 ptop = sump(1,1) ptop = ak(1) levs(1) = ptop + 0.5 * dpref(1) do L = 2, lm levs(L) = levs(L-1) + 0.5 * ( dpref(L-1) + dpref(L) ) enddo levs(1:lm) = levs(1:lm) / 100.0 ! Create Dry Temperature and SLP ! ------------------------------ do L=1,lm+1 ple(:,:,L) = ak(L) + ps(:,:)*bk(L) enddo pke(:,:,:) = ple(:,:,:)**kappa phi(:,:,lm+1) = phis(:,:) do L=lm,1,-1 phi(:,:,L) = phi(:,:,L+1) + cp*thv(:,:,L)*( pke(:,:,L+1)-pke(:,:,L) ) pk(:,:,L) = ( pke(:,:,L+1)-pke(:,:,L) ) . / ( kappa*log(ple(:,:,L+1)/ple(:,:,L)) ) enddo t = thv*pk call get_slp ( ps,phis,slp,ple,pk,t,rgas,grav,im,jm,lm ) t = t/(1+eps*q) phi = phi/grav ! Create GFIO file ! ---------------- call GFIO_Create ( filename, title, source, contact, undef, . im, jm, lm, lons, lats, levs, levunits, . nymd, nhms, timeinc, . nvars, vname, vtitle, vunits, lmvar, . v_range, p_range, precision, . fid, rc ) ! Write GFIO data ! --------------- allocate( dum1(lm+1) ) allocate( dum2(im,jm) ) allocate( dum3(im,jm,lm) ) dum2 = phis ; call Gfio_putVar ( fid,vname(01),nymd,nhms,im,jm,0, 1,dum2,rc ) dum2 = topo_std ; call Gfio_putVar ( fid,vname(02),nymd,nhms,im,jm,0, 1,dum2,rc ) dum2 = ts ; call Gfio_putVar ( fid,vname(03),nymd,nhms,im,jm,0, 1,dum2,rc ) dum2 = oro ; call Gfio_putVar ( fid,vname(04),nymd,nhms,im,jm,0, 1,dum2,rc ) dum2 = ps ; call Gfio_putVar ( fid,vname(05),nymd,nhms,im,jm,0, 1,dum2,rc ) dum3 = dp ; call Gfio_putVar ( fid,vname(06),nymd,nhms,im,jm,1,lm,dum3,rc ) dum3 = u ; call Gfio_putVar ( fid,vname(07),nymd,nhms,im,jm,1,lm,dum3,rc ) dum3 = v ; call Gfio_putVar ( fid,vname(08),nymd,nhms,im,jm,1,lm,dum3,rc ) dum3 = thv ; call Gfio_putVar ( fid,vname(09),nymd,nhms,im,jm,1,lm,dum3,rc ) dum3 = q ; call Gfio_putVar ( fid,vname(10),nymd,nhms,im,jm,1,lm,dum3,rc ) dum3 = oz ; call Gfio_putVar ( fid,vname(11),nymd,nhms,im,jm,1,lm,dum3,rc ) dum2 = slp ; call Gfio_putVar ( fid,vname(12),nymd,nhms,im,jm,0, 1,dum2,rc ) dum3 = t ; call Gfio_putVar ( fid,vname(13),nymd,nhms,im,jm,1,lm,dum3,rc ) dum3 = phi(:,:,1:lm) ; call Gfio_putVar ( fid,vname(14),nymd,nhms,im,jm,1,lm,dum3,rc ) ! Write GFIO global attributes ! ---------------------------- ks = 0 pint = ak(ks+1) call GFIO_PutIntAtt ( fid,'nstep', 1,ks ,0 ,rc ) call GFIO_PutRealAtt ( fid,'ptop', 1,ptop ,precision,rc ) call GFIO_PutRealAtt ( fid,'pint', 1,pint ,precision,rc ) call GFIO_PutIntAtt ( fid,'ks', 1,ks ,0 ,rc ) dum1 = ak ; call GFIO_PutRealAtt ( fid,'ak', lm+1,dum1 ,precision,rc ) dum1 = bk ; call GFIO_PutRealAtt ( fid,'bk', lm+1,dum1 ,precision,rc ) call gfio_close ( fid,rc ) return end subroutine sfchdf ( filename,nymd,nhms,im,jm,lm, . phis,u10m,v10m,ts,qs, . snomas,wet1,tsoil1, . frland,frlandice,frlake,frocean,oro,precision ) implicit none integer im,jm,lm,nymd,nhms real*8 phis(im,jm) real*8 u10m(im,jm) real*8 v10m(im,jm) real*8 ts(im,jm) real*8 qs(im,jm) real*8 snomas(im,jm) real*8 wet1(im,jm) real*8 tsoil1(im,jm) real*8 frland(im,jm) real*8 frlandice(im,jm) real*8 frlake(im,jm) real*8 frocean(im,jm) real*8 oro(im,jm) real lats(jm) real lons(im) real levs(lm) real*8 latsd(jm) real*8 lonsd(im) real*8 ptop,dlon,dlat,pref,undef integer i,j,L,n,timeinc,rc character*80 filename character*3 ch3 integer nvars,fid,precision character*256 levunits character*256 title character*256 source character*256 contact character*256, allocatable :: vname(:) character*256, allocatable :: vtitle(:) character*256, allocatable :: vunits(:) integer, allocatable :: lmvar(:) real, allocatable :: v_range(:,:) real, allocatable :: p_range(:,:) undef = 1.0e15 timeinc = 060000 ! String and vars settings ! ------------------------ title = 'GEOS5 Dynamics State Vector (Hybrid Coordinates)' source = 'Goddard Modeling and Assimilation Office, NASA/GSFC' contact = 'data@gmao.gsfc.nasa.gov' levunits = 'hPa' nvars = 13 allocate ( vname(nvars) ) allocate ( vtitle(nvars) ) allocate ( vunits(nvars) ) allocate ( lmvar(nvars) ) allocate ( v_range(2,nvars) ) allocate ( p_range(2,nvars) ) vname(01) = 'PHIS' vtitle(01) = 'Topography geopotential' vunits(01) = 'meter2/sec2' lmvar(01) = 0 vname(02) = 'TSKIN' vtitle(02) = 'Surface skin temperature' vunits(02) = 'K' lmvar(02) = 0 vname(03) = 'QS' vtitle(03) = 'Surface Specific Humidity' vunits(03) = 'kg/kg' lmvar(03) = 0 vname(04) = 'U10M' vtitle(04) = '10 meter U Wind' vunits(04) = 'm/s' lmvar(04) = 0 vname(05) = 'V10M' vtitle(05) = '10 meter V Wind' vunits(05) = 'm/s' lmvar(05) = 0 vname(06) = 'SNOWDP' vtitle(06) = 'Snow Depth Water Equivalent' vunits(06) = 'm' lmvar(06) = 0 vname(07) = 'GWETTOP' vtitle(07) = 'Top soil layer wetness' vunits(07) = '0-1' lmvar(07) = 0 vname(08) = 'TSOIL1' vtitle(08) = 'Top soil layer wetness' vtitle(08) = 'Surface skin temperature' vunits(08) = 'K' lmvar(08) = 0 vname(09) = 'FRLAND' vtitle(09) = 'Fraction of Land' vunits(09) = '0-1' lmvar(09) = 0 vname(10) = 'FRLANDICE' vtitle(10) = 'Fraction of Land Ice' vunits(10) = '0-1' lmvar(10) = 0 vname(11) = 'FRLAKE' vtitle(11) = 'Fraction of Lake' vunits(11) = '0-1' lmvar(11) = 0 vname(12) = 'FROCEAN' vtitle(12) = 'Fraction of Ocean' vunits(12) = '0-1' lmvar(12) = 0 vname(13) = 'ORO' vtitle(13) = 'Water(0) Land(1) Ice(2) Flag' vunits(13) = '0-1-2' lmvar(13) = 0 v_range(:,:) = undef p_range(:,:) = undef ! Compute grid ! ------------ dlon = 360.0/ im dlat = 180.0/(jm-1) do j=1,jm latsd(j) = -90.0 + (j-1)*dlat enddo do i=1,im lonsd(i) = 0.0 + (i-1)*dlon enddo do L=1,lm levs(L) = L enddo lats = latsd lons = lonsd ! Create GFIO file ! ---------------- call GFIO_Create ( filename, title, source, contact, undef, . im, jm, lm, lons, lats, levs, levunits, . nymd, nhms, timeinc, . nvars, vname, vtitle, vunits, lmvar, . v_range, p_range, precision, . fid, rc ) ! Write GFIO data ! --------------- call Gfio_putVar ( fid,vname(01),nymd,nhms,im,jm,0, 1,phis ,rc ) call Gfio_putVar ( fid,vname(02),nymd,nhms,im,jm,0, 1,ts ,rc ) call Gfio_putVar ( fid,vname(03),nymd,nhms,im,jm,0, 1,qs ,rc ) call Gfio_putVar ( fid,vname(04),nymd,nhms,im,jm,0, 1,u10m ,rc ) call Gfio_putVar ( fid,vname(05),nymd,nhms,im,jm,0, 1,v10m ,rc ) call Gfio_putVar ( fid,vname(06),nymd,nhms,im,jm,0, 1,snomas ,rc ) call Gfio_putVar ( fid,vname(07),nymd,nhms,im,jm,0, 1,wet1 ,rc ) call Gfio_putVar ( fid,vname(08),nymd,nhms,im,jm,0, 1,tsoil1 ,rc ) call Gfio_putVar ( fid,vname(09),nymd,nhms,im,jm,0, 1,frland ,rc ) call Gfio_putVar ( fid,vname(10),nymd,nhms,im,jm,0, 1,frlandice,rc ) call Gfio_putVar ( fid,vname(11),nymd,nhms,im,jm,0, 1,frlake ,rc ) call Gfio_putVar ( fid,vname(12),nymd,nhms,im,jm,0, 1,frocean ,rc ) call Gfio_putVar ( fid,vname(13),nymd,nhms,im,jm,0, 1,oro ,rc ) ! Write GFIO global attributes ! ---------------------------- ! call GFIO_PutRealAtt ( fid,'ptop', 1,ptop ,precision,rc ) call gfio_close ( fid,rc ) return end subroutine hflip ( q,im,jm,lm ) implicit none integer im,jm,lm,i,j,L real*8 q(im,jm,lm),dum(im) do L=1,lm do j=1,jm do i=1,im/2 dum(i) = q(i+im/2,j,L) dum(i+im/2) = q(i,j,L) enddo q(:,j,L) = dum(:) enddo enddo return end subroutine get_slp ( ps,phis,slp,pe,pk,tv,rgas,grav,im,jm,km ) implicit none integer im,jm,km real grav real rgas real*8 pk(im,jm,km) ! layer-mean P**kappa real*8 tv(im,jm,km) ! layer-mean virtual Temperature real*8 pe(im,jm,km+1) ! press at layer edges (Pa) real*8 ps(im,jm) ! surface pressure (Pa) real*8 phis(im,jm) ! surface geopotential real*8 slp(im,jm) ! sea-level pressure (hPa) real p_offset real p_bot real tstar ! extrapolated temperature (K) real tref ! Reference virtual temperature (K) real pref ! Reference pressure level (Pa) real pkref ! Reference pressure level (Pa) ** kappa real dp1, dp2 real factor, yfactor real gg real gamma integer k_bot, k, k1, k2, i,j gamma = 6.5e-3 gg = gamma / grav factor = grav / ( Rgas * gamma ) yfactor = Rgas * gg p_offset = 15000. ! 150 hPa above surface do j=1,jm do i=1,im p_bot = ps(i,j) - p_offset k_bot = -1 do k = km, 2, -1 if ( pe(i,j,k+1) .lt. p_bot ) then k_bot = k go to 123 endif enddo 123 continue k1 = k_bot - 1 k2 = k_bot dp1 = pe(i,j,k_bot) - pe(i,j,k_bot-1) dp2 = pe(i,j,k_bot+1) - pe(i,j,k_bot) pkref = ( pk(i,j,k1)*dp1 + pk(i,j,k2)*dp2 ) / (dp1+dp2) tref = ( tv(i,j,k1)*dp1 + tv(i,j,k2)*dp2 ) / (dp1+dp2) pref = 0.5 * ( pe(i,j,k_bot+1) + pe(i,j,k_bot-1) ) tstar = tref*( ps(i,j)/pref )**yfactor slp(i,j) = ps(i,j)*( 1.0+gg*phis(i,j)/tstar )**factor enddo enddo return end subroutine usage() print *, "Usage: " print * print *, " rs2hdf.x [-dynrst dynrst_fname] Default: fvcore_internal_restart" print *, " [-moistrst moistrst_fname] Default: moist_internal_restart" print *, " [-topo topo_fname] Optional Topography File" print *, " [-tag tag] Optional Tag" print *, " [-nymd yyyymmdd] Optional Overriding DateStamp for Output" print *, " [-nhms hhmmss] Optional Overriding TimeStamp for Output" print *, " [-ext filename extension] Optional Overriding Output filename extension" print * print *, "where:" print *, "-----" print *, " -dynrst dynrst_fname: Filename of dynamics internal restart" print *, " -moistrst moistrst_fname: Filename of moist internal restart" print * print *, "creates:" print *, "-------" print *, " tag.bkg.eta.yyyymmdd_hhz.nc4" print * call exit(7) end