program create_firstguess implicit none integer(4),parameter:: nt=10! # of horiz. slabs from wgrib integer(4),parameter:: nx=1073,ny=689,nsig=1 real(4),parameter::ds=5079.406 real(4),parameter::elonv=265.000000 real(4),parameter:: alatan=25.000000 real(4),parameter:: gravity=9.81 integer(4) iyear,imonth,iday,ihour,iminute,isecond integer(4) ihdrbuf(512) integer(4) ifield(nx,ny) integer(4) ivgtyp(nx,ny) integer(4) isltyp(nx,ny) integer(4) i,j,n real(4) glat(nx,ny) !Earth latitude real(4) glon(nx,ny) ! Earth longitude real(4) field(nx,ny) real(4) sst(nx,ny) real(4) tslb(nx,ny) real(4) tsk(nx,ny) real(4) vegfra(nx,ny) real(4) smois(nx,ny) real(4) qgrid(nx,ny) real(4) tdgrid(nx,ny) real(4) psfcgrid(nx,ny) real(4) ugrid(nx,ny) real(4) vgrid(nx,ny) real(4) uearth(nx,ny),ue,u real(4) vearth(nx,ny),ve,v real(4) udiff(nx,ny) real(4) vdiff(nx,ny) real(4) landmask(nx,ny) real(4) dx(nx,ny),dy(nx,ny) real(4) mapfac(nx,ny) real(4) xn,dg2rad,ylon real(4) rmax,rmin namelist/timeinfo/iyear,imonth,iday,ihour,iminute,isecond data iyear/2005/ data imonth/6/ data iday/17/ data ihour/12/ data iminute/0/ data isecond/0/ read(9,timeinfo) print*,'in create_firstguess,iyear,imonth,iday,ihour,iminute,isecond=',& iyear,imonth,iday,ihour,iminute,isecond !********************************************************************* print*,'in create_firstguess, nx,ny=',nx,ny !=> generate file in "first-guess format" for surface analysis rewind(30) read(30) glat read(30) glon read(30) mapfac read(40) landmask dx=ds*mapfac dy=ds*mapfac do i=1,512 ihdrbuf(i)=0 enddo write(88) ihdrbuf write(88) iyear,imonth,iday,ihour,iminute,isecond,nx,ny,nsig write(88) dx,dy write(88) glat write(88) glon rewind(20) do 300 n=1,nt print*,'read slab# ',n,' in read_write_ndfd_bckg' read (20) field print*,'rmin,rmax=',minval(field),maxval(field) if (n.eq.1) then psfcgrid=field write(88) psfcgrid ! psfc0 endif if (n.eq.2) then do j=1,ny do i=1,nx if (landmask(i,j).le.0.5) field(i,j)=field(i,j)-500. enddo enddo endif if (n.eq.2) write(88) field*gravity ! PHB (zsfc*g) if (n.eq.3) then write(88) field ! T(k) ! TEMP (sensible) sst=field tslb=field tsk=field do j=1,ny do i=1,nx if (landmask(i,j).gt.0.5) sst(i,j)=0. enddo enddo endif if (n.eq.4) tdgrid=field if (n.eq.5) ugrid=field if (n.eq.6) vgrid=field if (n.eq.7) then qgrid=field call td_to_q(tdgrid,psfcgrid,qgrid,nx,ny) !when in use, makes preceding statement redundant write(88) qgrid ! Q(k) endif 300 continue !==> rotate winds dg2rad=atan(1.)/45. xn=sin(alatan*dg2rad) do 400 j=1,ny do 400 i=1,nx u=ugrid(i,j) v=vgrid(i,j) ylon=glon(i,j) call get_earth_winds(ue,ve,u,v,xn,elonv-360.,ylon,dg2rad) uearth(i,j)=ue vearth(i,j)=ve udiff(i,j)=uearth(i,j)-ugrid(i,j) vdiff(i,j)=vearth(i,j)-vgrid(i,j) 400 continue print*,'for udiff/ rmin,rmax=',minval(udiff),maxval(udiff) print*,'for vdiff/ rmin,rmax=',minval(vdiff),maxval(vdiff) field=0. ifield=0 write(88)ugrid ! U(K) write(88)vgrid ! V(K) write(88)landmask ! LANDMASK (0=water and >0.5 for land) write(88)field ! XICE write(88)sst ! SST write(88)ifield ! IVGTYP write(88)ifield ! ISLTYP write(88)field ! VEGFRA write(88)field ! SNOW write(88)ugrid ! U10 write(88)vgrid ! V10 write(88)field ! SMOIS write(88)tslb ! TSLB write(88)tsk ! TSK !==> generate slabs_new.dat ! rewind(20) ! do 500 n=1,nt ! print*,'read slab# ',n,' in read_write_ndfd_bckg' ! read (20) field ! if (n.eq.5) field=uearth ! if (n.eq.6) field=vearth ! write(50) field !500 continue end ! ************************************************************** subroutine get_earth_winds(un,vn,ut,vt,ROTCON_P,LON_XX_P, & olon,dg2rad) ! ! Abstract: Rotates wind relative to grid on Lambert Conformal ! or polar streo projection to get winds with respect ! to the true north. ! Adapted from Stan Benjamin's code snippet implicit none real(4) un,vn,ut,vt,ROTCON_P,LON_XX_P,olon, & dg2rad,angle2,sinx2,cosx2 !** ROTCON_P R WIND ROTATION CONSTANT, = 1 FOR POLAR STEREO !** AND SIN(LAT_TAN_P) FOR LAMBERT CONFORMAL !** LAT_TAN_P R LATITUDE AT LAMBERT CONFORMAL PROJECTION !** IS TRUE (DEG) !** LON_XX_P R MERIDIAN ALIGNED WITH CARTESIAN X-AXIS(DEG) angle2 = rotcon_p*(olon-lon_xx_p)*dg2rad sinx2 = sin(angle2) cosx2 = cos(angle2) un = cosx2*ut+sinx2*vt vn =-sinx2*ut+cosx2*vt return end !******************************************************************** subroutine td_to_q(td,p,q,nx,ny) ! ! abstract: given dewpt in K and pressure in Pa, compute ! specific humidity in kg/kg. this is done by inverting: ! ! qv=q(i,j)/(1.-q(i,j)) ! e=p(i,j)/100.*qv/(eps+qv) ! eln=alog(e) ! td(i,j) = (243.5*eln-440.8)/(19.48-eln)+273.15 !******************************************************************** implicit none real(4),parameter::alpha=243.5 real(4),parameter::beta=440.8 real(4),parameter::gamma=19.48 real(4),parameter::eta=273.15 real(4),parameter::eps=0.62197 !=Rd/Rv integer(4) nx,ny,i,j real(4),dimension(nx,ny)::td,p,q real(4) loge,e,qv do j=1,ny do i=1,nx loge=(gamma*(td(i,j)-eta)+beta)/(alpha+td(i,j)-eta) e=exp(loge) qv=100.*e*eps/(p(i,j)-100.*e) q(i,j)=qv/(1.+qv) enddo enddo return end