!############################# Change Log ################################## ! 5.0.0 ! !########################################################################### ! Copyright (C) 1990, 1995, 1999, 2000, 2003 - All Rights Reserved ! Regional Atmospheric Modeling System - RAMS !########################################################################### SUBROUTINE input_rawi (olat1,olat2,olon1,olon2) use ModDateUtils use isan_coms !--(DMK-CCATT-INI)---------------------------------------------------------------- use chem_isan_coms !--(DMK-CCATT-FIM)---------------------------------------------------------------- use rconstants implicit none real :: olat1,olat2,olon1,olon2 character(len=256) :: line,line2 character(len=32) :: cflags,tokens(40) real, dimension(maxlev) :: p,t,z,h,d,f,pp,zp,tp,hp,vz,uz,zz,rp,dz,fz integer :: iqflagsp(maxlev,4,3),iqflagsz(maxlev,3,3),iqfl(4) character(len=4), dimension(maxlev) :: qp,q character(len=1) :: iq1,iq2,iq3,iq4 character(len=8) :: idsta character(len=6) :: asta character(len=14) :: obsdate real(kind=8) :: nasecs,obssecs integer :: imarker,iver_up,nobs,k integer :: ntok,jdate,jt,lp,lz,i,j,jyr,jmo,jdy,kl,iqf,iq real :: xlat,xlon,elev call date_abs_secs (iproc_dates(natime),nasecs) call date_unmake_big (iyear,imonth,idate,ihour,iproc_dates(natime)) ihour=ihour/100 print 559,inrawi 559 format(/' Acquiring rawindsonde input file ',A) call rams_f_open(4,inrawi(1:len_trim(inrawi)),'FORMATTED','OLD','READ',0) read(4,*) imarker rewind 4 if(imarker.eq.999999) then read(4,*) imarker,iver_up else iver_up=1 endif do nobs=1,10000 do k=1,maxlev p(k)=0. t(k)=0. z(k)=0. h(k)=0. f(k)=0. pp(k)=0. zp(k)=0. tp(k)=0. hp(k)=0. vz(k)=0. uz(k)=0. zz(k)=0. rp(k)=0. dz(k)=0. fz(k)=0. enddo if(iver_up.eq.1) then read(4,'(a)',end=100,err=101) line call parse(line,tokens,ntok) read(tokens(1),*) jdate read(tokens(2),*) jt read(tokens(3),'(a)') idsta read(tokens(4),*) lp read(tokens(5),*) lz read(tokens(6),*) xlat read(tokens(7),*) xlon read(tokens(8),*) elev IF(lp.GT.maxlev.or. lz.gt.maxlev) then print*,'Rawindsonde read error!' print*,' Number of input levels greater than MAXLEV' print*,' lp,lz,maxlev:',lp,lz,maxlev stop 'rawin-maxlev' endif do k=1,lp read(4,80)p(k),z(k),t(k),h(k),((iqflagsp(k,i,j),j=1,3),i=1,4) enddo 80 format(2f12.3,f10.2,f10.4,2x,4(':',3i1)) do k=1,lz read(4,85)zz(k),fz(k),dz(k),((iqflagsz(k,i,j),j=1,3),i=1,3) enddo 85 format(f12.3,2f10.2,2x,3(':',3i1)) elseif(iver_up.eq.2) then read(4,'(a)',end=100,err=101) line call parse(line,tokens,ntok) read(tokens(1),*) jyr read(tokens(2),*) jmo read(tokens(3),*) jdy read(tokens(4),*) jt read(tokens(5),'(a)') idsta read(tokens(6),*) lp read(tokens(7),*) lz read(tokens(8),*) xlat read(tokens(9),*) xlon read(tokens(10),*) elev IF(lp.GT.maxlev.or. lz.gt.maxlev) then print*,'Rawindsonde read error!' print*,' Number of input levels greater than MAXLEV' print*,' lp,lz,maxlev:',lp,lz,maxlev stop 'rawin-maxlev' endif jdate=jyr*10000+jmo*100+jdy do k=1,lp read(4,*,end=100,err=101) p(k),iqfl(1),z(k),iqfl(2) & ,t(k),iqfl(3),h(k),iqfl(4) do i=1,4 iqflagsp(k,i,1)=iqfl(i)/100 iqflagsp(k,i,2)=mod(iqfl(i)/10,10) iqflagsp(k,i,3)=mod(iqfl(i),10) enddo enddo do k=1,lz read(4,*,end=100,err=101) zz(k),iqfl(1),fz(k),iqfl(2) & ,dz(k),iqfl(3) do i=1,3 iqflagsz(k,i,1)=iqfl(i)/100 iqflagsz(k,i,2)=mod(iqfl(i)/10,10) iqflagsz(k,i,3)=mod(iqfl(i),10) enddo enddo endif ! Check date call date_make_big(jyr,jmo,jdy,jt*100,obsdate) call date_abs_secs(obsdate,obssecs) !print*,'===> ',natime,iobswin !print*,'===> ',iproc_dates(natime),' ',nasecs ,iyear,imonth,idate,ihour !print*,'===> ',obsdate ,' ',obssecs,jyr ,jmo ,jdy ,jt*100 if((iobswin < 0 .and. & (obssecs .lt. nasecs+iobswin .or. obssecs .gt. nasecs)).or. & (iobswin > 0 .and. & (obssecs .lt. nasecs-iobswin .or. obssecs .gt. nasecs+iobswin)).or. & (iobswin .eq. 0 .and. nasecs .ne. obssecs)) then print*,'station ',idsta(1:len_trim(idsta)) & ,' excluded - out of obs window (s)',iobswin goto 987 !else ! print*,'station ',idsta(1:len_trim(idsta)) & ! ,' within obs window (s)',iobswin endif ! --------------------------------------------------------- ! Convert input variables to SI units and appropriate form ! Set missing values to 1E30. ! PP - Pressure (Pascals) ! TP - Temperature (Kelvin) ! RP - Relative humidity (fraction) ! ZP - Height of pressure levels (meters) ! ZZ - Height of wind levels (meters) ! DZ - Wind direction (degrees) ! FZ - Wind speed (m/s) ! XLAT - Station latitude ! XLON - Station longitude ! ELEV - Station elevation (meters) ! LP - Number of pressure levels ! LZ - Number of wind levels in height ! IDSTA- Station ID ! --------------------------------------------------------- do kl=1,lp iqf=1 do iq=1,3 if((iqflagsp(kl,1,iq).ne.5.and.iqflagsp(kl,1,iq).ne.0 & .and.iqflagsp(kl,1,iq).ne.8).or.p(kl).lt.-998.) iqf=0 enddo !print*,'kl:',kl,p(kl),iqf,iqflagsp(kl,1,1:3) if (iqf.eq.1) pp(kl)=p(kl) if (iqf.eq.0) pp(kl)=1.e30 iqf=1 do iq=1,3 if((iqflagsp(kl,2,iq).ne.5.and.iqflagsp(kl,2,iq).ne.0 & .and.iqflagsp(kl,2,iq).ne.8).or.z(kl).lt.-998.) iqf=0 enddo if (iqf.eq.1) zp(kl)=z(kl) if (iqf.eq.0) zp(kl)=1.e30 !MJW Make sure that either pressure or height is non-missing if(pp(kl) > 1.e20 .and. zp(kl) > 1.e20) then tp(kl)=1.e30 rp(kl)=1.e30 cycle endif iqf=1 do iq=1,3 if((iqflagsp(kl,3,iq).ne.5.and.iqflagsp(kl,3,iq).ne.0 & .and.iqflagsp(kl,3,iq).ne.8).or.t(kl).lt.-998.) iqf=0 enddo if (iqf.eq.1) tp(kl)=t(kl)+273.16 if (iqf.eq.0) tp(kl)=1.e30 iqf=1 do iq=1,3 if((iqflagsp(kl,4,iq).ne.5.and.iqflagsp(kl,4,iq).ne.0 & .and.iqflagsp(kl,4,iq).ne.8).or.h(kl).lt.-998.) iqf=0 enddo if (iqf.eq.1) rp(kl)=h(kl) if (iqf.eq.0) rp(kl)=1.e30 enddo do kl=1,lz iqf=1 do iq=1,3 if((iqflagsz(kl,1,iq).ne.5.and.iqflagsz(kl,1,iq).ne.0 & .and.iqflagsz(kl,1,iq).ne.8).or.zz(kl).lt.-998.) iqf=0 enddo if (iqf.eq.1) zz(kl)=zz(kl) if (iqf.eq.0) zz(kl)=1.e30 !MJW Make sure that height is non-missing if(zz(kl) > 1.e20) then fz(kl)=1.e30 dz(kl)=1.e30 cycle endif iqf=1 do iq=1,3 if((iqflagsz(kl,2,iq).ne.5.and.iqflagsz(kl,2,iq).ne.0 & .and.iqflagsz(kl,2,iq).ne.8).or.fz(kl).lt.-998.) iqf=0 enddo if (iqf.eq.1) fz(kl)=fz(kl) if (iqf.eq.0) fz(kl)=1.e30 iqf=1 do iq=1,3 if((iqflagsz(kl,3,iq).ne.5.and.iqflagsz(kl,3,iq).ne.0 & .and.iqflagsz(kl,3,iq).ne.8).or.dz(kl).lt.-998.) iqf=0 enddo if (iqf.eq.1) dz(kl)=dz(kl) if (iqf.eq.0) dz(kl)=1.e30 enddo !---------------------------------- call sndproc (lp,lz,xlat,xlon,elev,idsta & ,tp,zp,pp,rp,dz,fz,zz) 987 continue enddo 100 continue print*,'End of rawindsonde file:',inrawi goto 110 101 continue print*,'Error in rawindsonde file:',inrawi goto 110 102 continue print*,'Error opening rawindsonde file:',inrawi 110 continue close(4) return end SUBROUTINE input_rawi !*************************************************************************** subroutine sndproc (lp,lz,xlat,xlon,elev,idsta,tp,zp,pp,rp,dz,fz,zz) use isan_coms use rconstants implicit none real, dimension(*) :: tp,zp,pp,rp,dz,fz,zz integer :: lp,lz real :: xlat,xlon,elev character(len=*) :: idsta real, dimension(maxlev) :: uz,vz,thz integer :: not,non,llp,k,llz,kk,lbc real :: bchyd,rss,pio,tho,zso real, external :: rs do not=1,notsta if('r'//idsta(1:5).eq.notid(not)(1:6)) then print 31,idsta,xlat,xlon,elev,0,0 return endif enddo 31 format(' ----- Station omitted ',1X,A8,' Lat:',F8.2,' Lon:',F8.2 & ,' Elev:',F10.2,' P Levels:',I5,' Wind levels:',I5) nsta=nsta+1 !-------------------------------------------------------------------- !mjb added the following from 4a as the following fix does not handle ! pressure levels flagged as suspect by qc (501) ! Throw out any levels that have suspect pressure ! ALSO, don't read in any levels above the gridded data, since that leads ! to bulls-eyes around all soundings near model top llp=0 do k=1,lp if(pp(k).lt.1.e20) then if(pp(k) < levpr(nprz)*100.)then print 35,pp(k)/100.,levpr(nprz) 35 format(' discarding ',f6.1, & 'mb level since it is above gridded data top at ',i4,'mb') cycle endif llp=llp+1 tp(llp)=tp(k) zp(llp)=zp(k) pp(llp)=pp(k) rp(llp)=rp(k) endif enddo lp=llp ! Throw out any levels that have suspect winds llz=0 do k=1,lz if(zz(k).lt.1.e20.and.dz(k).lt.1.e20.and.fz(k).lt.1.e20) then llz=llz+1 zz(llz)=zz(K) fz(llz)=fz(K) dz(llz)=dz(K) !print*,'mjb w1',llz,zz(llz),fz(llz),dz(llz) endif enddo !print*,'mjb - wind levels reduced: ',lz,llz lz=llz !-------------------------------------------------------------------- ! arrange winds on height levels do k=1,lz if(dz(k).lt.1e19.and.fz(k).lt.1e19) call winduv (dz(k),fz(k),uz(k),vz(k)) enddo call sort3(zz,uz,vz,lz) up_lat(nsta)=xlat up_lon(nsta)=xlon up_top(nsta)=elev up_lp(nsta)=lp up_lz(nsta)=lz up_chstid(nsta)=idsta(1:8) ! fix missing values if possible do k=1,lp ! missing temperatures are: ! set to valid upper temperature if none below, ! interpolated linearly in logP if between valid values, ! left missing if no valid temperatures are above if(tp(k).gt.1e29.and.k.ne.lp) then !if(k > 1)then ! print*,'missing temp',k,lp,tp(k-1) !else ! print*,'missing temp',k,lp,tp(k) !endif kk=k+1 do while (kk.lt.lp.and.tp(kk).gt.1e29) kk=kk+1 enddo if (k.eq.1.and.tp(kk).lt.1e30) then tp(k)=tp(kk) elseif(tp(kk).lt.1e30) then tp(k)=tp(k-1)+(tp(kk)-tp(k-1))*log(pp(k)/pp(k-1))/log(pp(kk)/pp(k-1)) endif endif ! same for humidities if(rp(k).gt.1e29.and.k.ne.lp) then !if(k > 1)then ! print*,'missing rh',k,lp,rp(k-1) !else ! print*,'missing rh',k,lp,rp(k) !endif kk=k+1 do while (kk.lt.lp.and.rp(kk).gt.1e29) kk=kk+1 enddo if (k.eq.1.and.rp(kk).lt.1e30) then rp(k)=rp(kk) elseif(rp(kk).lt.1e30) then rp(k)=rp(k-1)+(rp(kk)-rp(k-1))*log(pp(k)/pp(k-1))/log(pp(kk)/pp(k-1)) endif endif enddo ! Now that we filled in as many temps and rh's we could, we will recompute heights. ! Find a boundary condition as first level at or below 500 mb. If we do ! not find a good level, we don't use the sounding. bchyd=50000. lbc=0 do k=lp,1,-1 if(pp(k) < 1.e19 .and. pp(k) >= bchyd .and. & tp(k) < 1.e19 .and. zp(k) < 1.e19) then lbc=k go to 52 endif enddo ! Didn't find a good level, so let's look up from 500mb... do k=1,lp if(pp(k) < 1.e19 .and. pp(k) < bchyd .and. & tp(k) < 1.e19 .and. zp(k) < 1.e19) then lbc=k go to 52 endif enddo nsta=nsta-1 print*,' sndproc: Could not find good rawindsonde z boundary level.' print*,' Discarding sounding:',idsta(1:8) return 52 continue ! Compute theta - if rh non-missing, compute vitual theta. do k=1,lp thz(k)=1.e30 if(tp(k) < 1.e19 .and. pp(k) < 1.e19) then thz(k)=tp(k)*(p00/pp(k))**rocp if(rp(k) < 1.e19) then rss=rs(pp(k),tp(k))*rp(k) thz(k)=thz(k)*(1.+.61*rss) endif endif enddo ! Recompute heights pio=cp*(pp(lbc)/p00)**rocp zso=zp(lbc) tho=thz(lbc) do k=lbc+1,lp zp(k)=1.e30 if(pp(k) < 1e19 .and. thz(k) < 1e19)then zp(k)=zso+.5*(thz(k)+tho)*(pio-cp*(pp(k)/p00)**rocp)/g zso=zp(k) pio=cp*(pp(k)/p00)**rocp tho=thz(k) endif enddo zso=zp(lbc) pio=cp*(pp(lbc)/p00)**rocp tho=thz(lbc) do k=lbc-1,1,-1 zp(k)=1.e30 if(pp(k) < 1e19 .and. thz(k) < 1e19)then zp(k)=zso+.5*(thz(k)+tho)*(pio-cp*(pp(k)/p00)**rocp)/g zso=zp(k) pio=cp*(pp(k)/p00)**rocp tho=thz(k) endif enddo ! Do one last filter and fill final arrays. ! If any heights are missing, discard the level. llp=0 do k=1,lp if(zp(k) < 1.e19) then llp=llp+1 up_t(nsta,llp)=tp(k) up_z(nsta,llp)=zp(k) up_p(nsta,llp)=pp(k) up_r(nsta,llp)=rp(k) endif enddo up_lp(nsta)=llp do k=1,lz up_uz(nsta,k)=uz(k) up_vz(nsta,k)=vz(k) up_zz(nsta,k)=zz(k) enddo print 1001,nsta,idsta,up_lat(nsta),up_lon(nsta),up_top(nsta),lp,lz 1001 format(' Rawindsonde station found ',I5,' ID:',1X,A8,' Lat:' & ,F8.2,' Lon:', F8.2,' Elev:',0PF10.2,' P levels:',I5 & ,' Wind levels:',I5) if(nsta.ge.maxsta) then print 92,maxsta 92 format(' Number of stations found is greater than MAXSTA',I6) stop 'st3-maxsta' endif return end subroutine sndproc !*************************************************************************** subroutine staprt (nsta,lp,lz,m1,m2,usndz,vsndz,zsndz,psnd,zsnd,rsnd,tsnd) use rconstants implicit none integer :: nsta,lp,lz,m1,m2 real :: usndz(m1,m2),vsndz(m1,m2),psnd(m1,m2),zsnd(m1,m2) & ,rsnd(m1,m2),tsnd(m1,m2),zsndz(m1,m2) integer :: k print 1003,(psnd(nsta,k),zsnd(nsta,k) & ,tsnd(nsta,k)*(p00/psnd(nsta,k))**rocp & ,tsnd(nsta,k),rsnd(nsta,k) ,k=1,lp) 1003 format(//' Pressure sounding'//,T5,'P(mb)',T15,'Height(m)' & ,T25,'Theta(K)' & ,T35,'Temp(K)',T45,'Rel hum(frac)',// & (1x,-2pf10.2,0p3f10.2,0pf10.3)) print 1004,(zsndz(nsta,k),usndz(nsta,k),vsndz(nsta,k),k=1,lz) 1004 format(//' Height levels',//,T5,'Z(m)',T15,'U(m/s)', & T25,'V(m/s)',//,(1X, 0P3F10.2)) print '(///)' return end subroutine staprt !*************************************************************************** SUBROUTINE input_sfc (olat1,olat2,olon1,olon2) use ModDateUtils use isan_coms use rconstants implicit none real :: olat1,olat2,olon1,olon2 character*1 itflg,iddd*4,idsta*8,line*256,line2*256,cflags*64 character(len=32) :: tokens(40),var_string integer :: iqflags(20,3),iqfl(20) real*8 nasecs integer :: iver_sfc,nvar,imarker integer :: ic,nv,ntok,nsss,k,jyr,jmo,jdy,jd,jt,iqf,iq real :: xlat,xlon,zx,ffx,ddx,tx,tdx,px call date_abs_secs (iproc_dates(natime),nasecs) call date_unmake_big (iyear,imonth,idate,ihour,iproc_dates(natime)) ihour=ihour/100 ! Read surface observations from a file print 559,insrfce 559 format(/' Acquiring surface input file ',A40) call rams_f_open(4,insrfce(1:len_trim(insrfce)),'FORMATTED','OLD','READ',0) !read(4,*) imarker !rewind 4 !if(imarker.eq.999999) then ! read(4,*) imarker,iver_sfc !else ! iver_sfc=1 !endif !read(4,*) nvar !do nv=1,nvar ! read(4,*) var_string !enddo !if (iver_sfc.eq.1) read(4,*) do nsss=1,1000000 read(4,'(a)',end=20,err=21) line ! Check to see if this is a header line. If so, (re)read header. if(line(1:6) == '999999') then !read(4,*) imarker ! rewind 4 ! if(imarker.eq.999999) then ! read(4,*) imarker,iver_sfc ! else ! iver_sfc=1 ! endif iver_sfc=2 read(4,*) nvar do nv=1,nvar read(4,*) var_string enddo !if (iver_sfc.eq.1) read(4,*) read(4,'(a)',end=20,err=21) line endif ! Read a surface obs if(iver_sfc.eq.1) then ! read(4,'(a)',end=20,err=21) line call parse(line,tokens,ntok) read(tokens(1),*) jd read(tokens(2),*) jt read(tokens(3),'(a)') idsta read(tokens(4),*) xlat read(tokens(5),*) xlon read(tokens(6),*) zx read(tokens(7),*) ffx read(tokens(8),*) ddx read(tokens(9),*) tx read(tokens(10),*) tdx read(tokens(11),*) px read(tokens(12),'(a)') cflags ic=1 do nv=1,nvar read(cflags(ic:),'(1x,3i1)') (iqflags(nv,k),k=1,3) ic=ic+4 enddo jyr=jd/10000 jmo=mod(jd,10000)/100 jdy=mod(jd,100) elseif (iver_sfc.eq.2) then ! read(4,'(a)',end=20,err=21) line call parse(line,tokens,ntok) read(tokens(1),*) jyr read(tokens(2),*) jmo read(tokens(3),*) jdy read(tokens(4),*) jt read(tokens(5),'(a)') idsta read(tokens(6),*) xlat read(tokens(7),*) xlon read(tokens(8),*) zx read(tokens(9),*) ffx read(tokens(10),*) iqfl(1) read(tokens(11),*) ddx read(tokens(12),*) iqfl(2) read(tokens(13),*) tx read(tokens(14),*) iqfl(3) read(tokens(15),*) tdx read(tokens(16),*) iqfl(4) read(tokens(17),*) px read(tokens(18),*) iqfl(5) do nv=1,nvar iqflags(nv,1)=iqfl(nv)/100 iqflags(nv,2)=mod(iqfl(nv)/10,10) iqflags(nv,3)=mod(iqfl(nv),10) enddo endif !------------------------------------------------------------------- ! Section for conversion to SI units !------------------------------------------------------------------- ! At the end of this section, the following variables are ! required to be in SI units. ! ! XLAT - STATION LATITUDE ! XLON - STATION LONGITUDE ! ZX - STATION ELEVATION ! DDX - WIND DIRECTION ! FFX - WIND SPEED ! TX - TEMPERATURE IN KELVIN ! TDX - DEWPOINT TEMPERATURE IN KELVIN ! PX - SURFACE PRESSURE ! IDSTA - Station ID (char*5) ! ! Any missing values must be set to 1.E30 !-------------------------------------------------------------------- if(zx.lt.-998.) goto 987 iqf=1 do iq=1,3 if((iqflags(1,iq).ne.5.and.iqflags(1,iq).ne.0 & .and.iqflags(1,iq).ne.8) .or. ffx.lt.-998.) iqf=0 enddo if (iqf.eq.1) ffx=ffx if (iqf.eq.0) ffx=1.e30 iqf=1 do iq=1,3 if((iqflags(2,iq).ne.5.and.iqflags(2,iq).ne.0 & .and.iqflags(2,iq).ne.8) .or. ddx.lt.-998.) iqf=0 enddo if (iqf.eq.1) ddx=ddx if (iqf.eq.0) ddx=1.e30 iqf=1 do iq=1,3 if((iqflags(3,iq).ne.5.and.iqflags(3,iq).ne.0 & .and.iqflags(3,iq).ne.8) .or. tx.lt.-998.) iqf=0 enddo if (iqf.eq.1) tx=tx+273.16 if (iqf.eq.0) tx=1.e30 iqf=1 do iq=1,3 if((iqflags(4,iq).ne.5.and.iqflags(4,iq).ne.0 & .and.iqflags(4,iq).ne.8) .or. tdx.lt.-998.) iqf=0 enddo if (iqf.eq.1) tdx=tdx+273.16 if (iqf.eq.0) tdx=1.e30 iqf=1 do iq=1,3 if((iqflags(5,iq).ne.5.and.iqflags(5,iq).ne.0 & .and.iqflags(5,iq).ne.8) .or. px.lt.-998.) iqf=0 enddo if (iqf.eq.1) px=px if (iqf.eq.0) px=1.e30 !----------------------------------- call sfcproc (nasecs,jyr,jmo,jdy,jt,idsta,xlat,xlon & ,zx,ddx,ffx,tx,tdx,px) 987 continue ENDDO 20 continue print*,'End of surface file:',insrfce goto 110 21 continue print*,'Error reading surface file:', insrfce goto 110 102 continue print*,'Error opening surface obs file',insrfce goto 110 110 continue close(4) return end SUBROUTINE input_sfc !*************************************************************************** subroutine sfcproc (nasecs,jyr,jmo,jdy,jt,idsta,xlat,xlon & ,zx,ddx,ffx,tx,tdx,px) use ModDateUtils use isan_coms use rconstants implicit none integer :: jyr,jmo,jdy,jt real :: xlat,xlon,zx,ddx,ffx,tx,tdx,px real*8 nasecs character(len=*) :: idsta real :: t1000,t900,t850,t700,t600,t500,p1000,p900,p850,p700,p600,p500 & ,z1000,z900,z850,z700,z600,z500,spd real*8 obssecs,obssecsns character(len=14) :: obsdate character(len=8) :: csfc integer :: nsu,ns,not,mis1,mis2,irepl,isfc real, external :: rs integer, save ::ncall=0 call date_make_big(jyr,jmo,jdy,jt*100,obsdate) call date_abs_secs(obsdate,obssecs) isfc=0 nssfc=nssfc+1 nsu=nssfc sf_chstid(nsu)=idsta sf_top(nsu)=zx sf_lat(nsu)=xlat sf_lon(nsu)=xlon if(ddx.lt.1.e19.and.ffx.lt.1.e19) then call winduv(ddx,ffx,sf_u(nsu),sf_v(nsu)) else sf_u(nsu)=1.e30 sf_v(nsu)=1.e30 endif sf_date(nsu)=obsdate ! check whether to omit station do not=1,notsta if('s'//idsta(1:5).eq.notid(not)(1:6)) then isfc=5 nssfc=nssfc-1 goto 110 endif enddo ! Check date against iobswin if((iobswin < 0 .and. & (obssecs .lt. nasecs+iobswin .or. obssecs .gt. nasecs)).or. & (iobswin > 0 .and. & (obssecs .lt. nasecs-iobswin .or. obssecs .gt. nasecs+iobswin)).or. & (iobswin .eq. 0 .and. nasecs .ne. obssecs)) then isfc=3 nssfc=nssfc-1 goto 110 endif t1000=287.6 t900=281.9 t850=278.6 t700=268.6 t600=260.8 t500=252.0 p1000=100000. p900=90000. p850=85000. p700=70000. p600=60000. p500=50000. z1000=111. z900=988. z850=1457. z700=3012. z600=4206. z500=5574. if(px.gt.1.e20.and.tx.lt.1.e20) then if(2.*zx.lt.z1000+z900) then px=p1000*exp(-g*(zx-z1000)/(rgas*.5*(tx+t1000))) elseif(2.*zx.lt.z900+z850) then px=p900*exp(-g*(zx-z900)/(rgas*.5*(tx+t900))) elseif(2.*zx.lt.z850+z700) then px=p850*exp(-g*(zx-z850)/(rgas*.5*(tx+t850))) elseif(2.*zx.lt.z700+z600) then px=p700*exp(-g*(zx-z700)/(rgas*.5*(tx+t700))) elseif(2.*zx.lt.z600+z500) then px=p600*exp(-g*(zx-z600)/(rgas*.5*(tx+t600))) endif endif if(tx.lt.1.e19.and.px.lt.1.e19) then if(tdx.lt.1.e19) then sf_r(nsu)=rs(px,tdx)/rs(px,tx) else sf_r(nsu)=1.e30 endif if(zx.lt.1.e19) then sf_s(nsu)=cp*tx+g*zx else sf_s(nsu)=1.e30 endif sf_t(nsu)=tx*(p00/px)**rocp sf_p(nsu)=px else sf_t(nsu)=1.e30 sf_p(nsu)=1.e30 sf_r(nsu)=1.e30 sf_s(nsu)=1.e30 endif do ns=1,nsu-1 ! determine if current date is closer for idsta if(sf_chstid(ns)==idsta) then call date_abs_secs(sf_date(ns),obssecsns) if(abs(nasecs-obssecsns) > abs(nasecs-obssecs)) then isfc=7 nssfc=nssfc-1 mis1=0 if(sf_u(ns).gt.1.e19) mis1=mis1+1 if(sf_v(ns).gt.1.e19) mis1=mis1+1 if(sf_t(ns).gt.1.e19) mis1=mis1+1 if(sf_p(ns).gt.1.e19) mis1=mis1+1 if(sf_r(ns).gt.1.e19) mis1=mis1+1 if(sf_top(ns).gt.1.e19) mis1=mis1+1 mis2=0 if(sf_u(nsu).gt.1.e19) mis2=mis2+1 if(sf_v(nsu).gt.1.e19) mis2=mis2+1 if(sf_t(nsu).gt.1.e19) mis2=mis2+1 if(sf_p(nsu).gt.1.e19) mis2=mis2+1 if(sf_r(nsu).gt.1.e19) mis2=mis2+1 if(sf_top(nsu).gt.1.e19) mis2=mis2+1 if(mis2.ge.mis1) then sf_u(ns)=sf_u(nsu) sf_v(ns)=sf_v(nsu) sf_t(ns)=sf_t(nsu) sf_s(ns)=sf_s(nsu) sf_p(ns)=sf_p(nsu) sf_top(ns)=sf_top(nsu) sf_r(ns)=sf_r(nsu) sf_lat(ns)=sf_lat(nsu) sf_lon(ns)=sf_lon(nsu) csfc=sf_chstid(ns) sf_chstid(ns)=sf_chstid(nsu) sf_date(ns)=sf_date(nsu) isfc=4 irepl=ns goto 110 endif endif endif enddo do ns=1,nsu-1 ! check for redundant surface obs if(abs(xlon-sf_lon(ns))