! !################################################################## !################################################################## !###### ###### !###### PROGRAM WRF2ARPS ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! PROGRAM wrfextsnd,79 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Extracts soundings from the WRF files and variables ! !----------------------------------------------------------------------- ! ! AUTHOR: Keith Brewster ! 2007 Feb 15 Based on WRF2ARPS by Yunheng Wang, Dan Weber ! and extsnd program from ARPS package. ! ! MODIFICATION HISTORY: ! ! Yunheng Wang (03/12/2007) ! Merged getwrfds with getwrfd so that it has the same interface as WRF2ARPS. ! Added a new parameter "grid_id" for WRF nesting outputs. ! Stopped the program when a station location exceeds the WRF domain. ! ! Yunheng Wang (03/26/2007) ! Removed namelist parameter frames_per_outfile. It is not determined by ! the program automatically. ! !----------------------------------------------------------------------- ! IMPLICIT NONE INTEGER :: nx_wrf,ny_wrf,nz_wrf INTEGER :: nzsoil_wrf,nstyp_wrf INTEGER :: iproj_wrf ! external data map projection REAL :: scale_wrf ! external data map scale factor REAL :: trlon_wrf ! external data true longitude REAL :: trlat1_wrf,trlat2_wrf,latnot_wrf(2) ! external data true latitude(s) REAL :: ctrlat_wrf, ctrlon_wrf REAL :: x0_wrf,y0_wrf ! external data origin REAL :: dx_wrf,dy_wrf,dt_wrf INTEGER :: sf_surface_physics, mp_physics ! !----------------------------------------------------------------------- ! ! WRF forecast variables ! !----------------------------------------------------------------------- ! REAL, ALLOCATABLE :: x_wrf(:) ! external data x-coordinate REAL, ALLOCATABLE :: y_wrf(:) ! external data y-coordinate REAL, ALLOCATABLE :: xu_wrf(:) ! external data u x-coordinate REAL, ALLOCATABLE :: yu_wrf(:) ! external data u y-coordinate REAL, ALLOCATABLE :: xv_wrf(:) ! external data v x-coordinate REAL, ALLOCATABLE :: yv_wrf(:) ! external data v y-coordinate REAL, ALLOCATABLE :: lat_wrf(:,:) ! external data latidude REAL, ALLOCATABLE :: lon_wrf(:,:) ! external data longitude REAL, ALLOCATABLE :: latu_wrf(:,:) ! external data latidude (x-stag) REAL, ALLOCATABLE :: lonu_wrf(:,:) ! external data longitude (x-stag) REAL, ALLOCATABLE :: latv_wrf(:,:) ! external data latidude (y-stag) REAL, ALLOCATABLE :: lonv_wrf(:,:) ! external data longitude (y-stag) REAL, ALLOCATABLE :: zp_wrf(:,:,:) ! external data physical height (m) REAL, ALLOCATABLE :: hgt_wrf(:,:,:) ! Height (m) of scalar points REAL, ALLOCATABLE :: zpsoil_wrf(:,:,:)! Height (m) of soil layers REAL, ALLOCATABLE :: p_wrf(:,:,:) ! Pressure (Pascals) REAL, ALLOCATABLE :: pt_wrf(:,:,:) ! Potential Temperature (K) REAL, ALLOCATABLE :: t_wrf(:,:,:) ! Temperature (K) REAL, ALLOCATABLE :: u_wrf(:,:,:) ! Eastward wind component REAL, ALLOCATABLE :: v_wrf(:,:,:) ! Northward wind component REAL, ALLOCATABLE :: vatu_wrf(:,:,:) ! NOTE: only used when use_wrf_grid /= 1 REAL, ALLOCATABLE :: uatv_wrf(:,:,:) ! REAL, ALLOCATABLE :: w_wrf(:,:,:) ! Vertical wind component REAL, ALLOCATABLE :: qv_wrf(:,:,:) ! Specific humidity (kg/kg) REAL, ALLOCATABLE :: qc_wrf(:,:,:) ! Cloud H2O mixing ratio (kg/kg) REAL, ALLOCATABLE :: qr_wrf(:,:,:) ! Rain H2O mixing ratio (kg/kg) REAL, ALLOCATABLE :: qi_wrf(:,:,:) ! Ice H2O mixing ratio (kg/kg) REAL, ALLOCATABLE :: qs_wrf(:,:,:) ! Snow H2O mixing ratio (kg/kg) REAL, ALLOCATABLE :: qh_wrf(:,:,:) ! Hail H2O mixing ratio (kg/kg) REAL, ALLOCATABLE :: psfc_wrf(:,:) ! Surface Pressure REAL, ALLOCATABLE :: t2m_wrf(:,:) ! Temperature (K) at 2 m AGL REAL, ALLOCATABLE :: th2m_wrf(:,:) ! Potential Temperature (K) at 2 m AGL REAL, ALLOCATABLE :: qv2m_wrf(:,:) ! Specific Humidity (kg/kg) at 2 m AGL REAL, ALLOCATABLE :: u10m_wrf(:,:) ! U wind component at 10 m AGL REAL, ALLOCATABLE :: v10m_wrf(:,:) ! V wind component at 10 m AGL REAL, ALLOCATABLE :: rainf(:,:) ! Total accumulated rainfall REAL, ALLOCATABLE :: cumrain(:,:) ! Total accumulated cumulus rainfall REAL, ALLOCATABLE :: snowf(:,:) ! Total accumulated snow and ice REAL, ALLOCATABLE :: graupel(:,:) ! Total accumulated graupel REAL, ALLOCATABLE :: tsoil_wrf (:,:,:,:) ! Soil temperature (K) REAL, ALLOCATABLE :: qsoil_wrf (:,:,:,:) ! Soil moisture (m3/m3) REAL, ALLOCATABLE :: wetcanp_wrf(:,:,:) ! Canopy water amount REAL, ALLOCATABLE :: snowdpth_wrf(:,:) ! Snow depth (m) REAL, ALLOCATABLE :: trn_wrf (:,:) ! External terrain (m) REAL, ALLOCATABLE :: raing_wrf (:,:) ! PBL time-step total precipitation (mm) REAL, ALLOCATABLE :: rainc_wrf (:,:) ! time-step cumulus precipitation (mm) INTEGER, ALLOCATABLE :: soiltyp_wrf (:,:,:) ! Soil type REAL, ALLOCATABLE :: stypfrct_wrf(:,:,:) ! Soil type fraction INTEGER, ALLOCATABLE :: vegtyp_wrf (:,:) ! Vegetation type REAL, ALLOCATABLE :: veg_wrf (:,:) ! Vegetation fraction !------------------------------------------------------------------- ! ! Working arrays ! !------------------------------------------------------------------- REAL, ALLOCATABLE :: tc(:) REAL, ALLOCATABLE :: tdc(:) REAL, ALLOCATABLE :: usnd(:) REAL, ALLOCATABLE :: vsnd(:) REAL, ALLOCATABLE :: dir(:) REAL, ALLOCATABLE :: spd(:) REAL, ALLOCATABLE :: dxfld(:) ! on WRF grid REAL, ALLOCATABLE :: dyfld(:) REAL, ALLOCATABLE :: rdxfld(:) REAL, ALLOCATABLE :: rdyfld(:) REAL, ALLOCATABLE :: dxfldu(:) REAL, ALLOCATABLE :: dyfldu(:) REAL, ALLOCATABLE :: rdxfldu(:) REAL, ALLOCATABLE :: rdyfldu(:) REAL, ALLOCATABLE :: dxfldv(:) REAL, ALLOCATABLE :: dyfldv(:) REAL, ALLOCATABLE :: rdxfldv(:) REAL, ALLOCATABLE :: rdyfldv(:) REAL, ALLOCATABLE :: tem1_wrf(:,:,:) ! Temporary work array REAL, ALLOCATABLE :: tem2_wrf(:,:,:) ! Temporary work array REAL, ALLOCATABLE :: tem3_wrf(:,:,:) ! Temporary work array REAL, ALLOCATABLE :: tem4_wrf(:,:,:) ! Temporary work array REAL, ALLOCATABLE :: tem5_wrf(:,:,:) ! Temporary work array REAL, ALLOCATABLE :: xa_wrf(:,:) ! WRF x coordinate on ARPS grid REAL, ALLOCATABLE :: ya_wrf(:,:) ! WRF y coordinate on ARPS grid ! !----------------------------------------------------------------------- ! ! include files ! !----------------------------------------------------------------------- ! INCLUDE 'globcst.inc' INCLUDE 'grid.inc' INCLUDE 'phycst.inc' INCLUDE 'bndry.inc' INCLUDE 'adjust.inc' INCLUDE 'mp.inc' INTEGER, PARAMETER :: maxwrffil = 100 INTEGER, PARAMETER :: maxsnd = 500 INTEGER, PARAMETER :: maxtime = 1000 ! !----------------------------------------------------------------------- ! ! NAMELIST parameters ! !----------------------------------------------------------------------- ! INTEGER :: nprocx_in, nprocy_in CHARACTER(LEN=256) :: dir_extd ! directory of external data CHARACTER(LEN=256) :: dir_output INTEGER :: io_form LOGICAL :: multifile CHARACTER(LEN=19) :: init_time_str,start_time_str,end_time_str CHARACTER(LEN=11) :: history_interval INTEGER :: outsnd INTEGER :: outsfc REAL :: ugrid REAL :: vgrid INTEGER :: locopt INTEGER :: nsnd REAL :: xpt(maxsnd) REAL :: ypt(maxsnd) CHARACTER (LEN=512) :: uafile CHARACTER (LEN=512) :: sfcfile CHARACTER (LEN=8) :: stid(maxsnd) REAL :: slat(maxsnd) REAL :: slon(maxsnd) REAL :: selev(maxsnd) INTEGER :: ipt(maxsnd) INTEGER :: jpt(maxsnd) INTEGER :: istnm(maxsnd) LOGICAL :: ptingrid(maxsnd) REAL :: sfcpr(maxtime,maxsnd) REAL :: tsfcf(maxtime,maxsnd) REAL :: tdsfcf(maxtime,maxsnd) REAL :: wdir(maxtime,maxsnd) REAL :: wspd(maxtime,maxsnd) REAL :: prcgp(maxtime,maxsnd) REAL :: prccu(maxtime,maxsnd) REAL :: prctot(maxtime,maxsnd) REAL :: tmax(maxsnd) REAL :: tmin(maxsnd) REAL :: spdmax(maxsnd) REAL :: prc0(maxsnd) REAL :: prcsum(maxsnd) INTEGER :: dmp_out_joined INTEGER :: grid_id NAMELIST /message_passing/ nproc_x, nproc_y, max_fopen, & nprocx_in, nprocy_in NAMELIST /wrfdfile/ dir_extd,init_time_str,io_form,grid_id, & start_time_str,history_interval,end_time_str NAMELIST /sndloc/ ugrid,vgrid,locopt,nsnd, & slat,slon,selev,xpt,ypt,stid,istnm NAMELIST /output/ dir_output,outsnd,outsfc INTEGER :: ncompressx, ncompressy ! compression in x and y direction: ! ncompressx=nprocx_in/nproc_x ! ncompressy=nprocy_in/nproc_y !----------------------------------------------------------------------- ! ! Non-dimensional smoothing coefficient ! !----------------------------------------------------------------------- ! REAL, PARAMETER :: smfct1 = 0.5 REAL, PARAMETER :: smfct2 = -0.5 REAL, PARAMETER :: rhmax = 1.0 ! !----------------------------------------------------------------------- ! ! Misc internal variables ! !----------------------------------------------------------------------- ! REAL, PARAMETER :: ms2kts =1.9438 REAL, PARAMETER :: mm2inch=(1.0/25.4) INTEGER, PARAMETER :: isec24h = 86400 CHARACTER (LEN=300) :: systring INTEGER :: isnow,jsnow,ii,jj INTEGER :: is, idummy INTEGER :: kbot,ktop REAL :: amin,amax REAL :: xumin,xumax,yvmin,yvmax REAL :: qvmin,qvmax,qvval REAL :: scrhgt,annhgt,csconst,pconst,tv_wrf,tsfc,tf,tdf,tdsfc REAL :: tcsoil,delt,delz,spdkts,raing,rainc,raint REAL :: pres,temp,qvsat,rh,tvbar,qvprt,qtot,dirsfc,spdsfc REAL :: tctof,tktof CHARACTER(LEN=256) :: extdname(MAXWRFFIL) CHARACTER(LEN=256) :: tmpstr INTEGER, ALLOCATABLE :: fHndl(:,:) CHARACTER(LEN=19) :: timestr CHARACTER(LEN=15) :: timefil(maxtime) CHARACTER(LEN=17) :: timewrf(maxtime) REAL :: latnot(2) REAL :: deltaz,omegasnd INTEGER :: i,j,k,ksmth INTEGER :: strlen INTEGER :: istatus INTEGER :: nextdfil INTEGER :: ifile,itime,jtime,ntime,isnd INTEGER :: iniotfu INTEGER :: iextmn,iextmx,jextmn,jextmx INTEGER :: iyr,imo,iday,ihr,imin,isec,jldy INTEGER :: myr,modyr,initsec,iabssec,istverif,kftime LOGICAL :: rewindyr LOGICAL :: fexist LOGICAL :: first_time LOGICAL :: in_verif LOGICAL :: wait_verif LOGICAL :: exit_early CHARACTER(LEN=1) :: ach INTEGER :: idist DOUBLE PRECISION :: ntmergeinv, mfac REAL :: ppasc,pmb,theta,smix,e,bige,alge REAL :: dd,dmin,latd,lond INTEGER :: abstimes,abstimee,abstimei INTEGER :: iloc, jloc INTEGER :: frames_per_outfile(MAXWRFFIL) ! !----------------------------------------------------------------------- ! ! roufns & lai convert table (see src/arpssfc/sfc_winter.tbl) ! !----------------------------------------------------------------------- REAL, PARAMETER :: rfnstbl(14) = (/0.002, 0.02, 0.01, 0.03, 0.10, & 0.20, 0.40, 2.00, 0.005,0.01, & 0.02, 0.06, 0.04, 0.002 /) REAL, PARAMETER :: laitbl(14) = (/0.5, 0.5, 0.02, 0.02, 0.02, & 0.05, 0.5, 0.5, 0.0, 0.02, & 0.0, 0.5, 0.5, 0.0 /) ! !----------------------------------------------------------------------- ! ! Function f_qvsat and inline directive for Cray PVP ! !----------------------------------------------------------------------- ! REAL :: f_qvsat ! compute saturation specific humidity defined in ! thermolib3d.f90 !fpp$ expand (f_qvsat) !!dir$ inline always f_qvsat !*$* inline routine (f_qvsat) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! !----------------------------------------------------------------------- ! ! MIsc Initializations ! !----------------------------------------------------------------------- ! jtime=0 ntime=0 tmax=-999.0 tmin=999.0 spdmax=-999.0 prc0=-99.0 prcsum=0.0 mp_physics=2 ! !----------------------------------------------------------------------- ! ! Initialize message passing processors. ! !----------------------------------------------------------------------- ! ! Non-MPI defaults: All others initialized in mpinit_var mp_opt = 0 myproc = 0 nprocx_in = 1 nprocy_in = 1 dumpstride = 1 readstride = 1 CALL mpinit_proc IF(myproc == 0) WRITE(6,'(10(/5x,a),/)') & '###################################################################',& '###################################################################',& '#### ####',& '#### Welcome to WRFEXTSND ####',& '#### ####',& '#### Program that reads in output from the WRF model ####',& '#### and produces text soundings. ####',& '#### ####',& '###################################################################',& '###################################################################' mgrid = 1 nestgrd = 0 ! !----------------------------------------------------------------------- ! ! Read in message passing options. ! !----------------------------------------------------------------------- ! IF (myproc == 0) THEN READ (5,message_passing) WRITE(6,'(a)')' Namelist block message_passing sucessfully read.' END IF CALL mpupdatei(nproc_x,1) CALL mpupdatei(nproc_y,1) CALL mpupdatei(max_fopen,1) CALL mpupdatei(nprocx_in,1) CALL mpupdatei(nprocy_in,1) ! !----------------------------------------------------------------------- ! ! Initialize message passing variables. ! !----------------------------------------------------------------------- ! CALL mpinit_var ! !----------------------------------------------------------------------- ! ! Read in namelist &wrfdfile ! !----------------------------------------------------------------------- ! dir_extd = './' init_time_str = '0000-00-00_00:00:00' start_time_str = '0000-00-00_00:00:00' history_interval = '00_00:00:00' end_time_str = '0000-00-00_00:00:00' io_form = 7 grid_id = 1 multifile = .FALSE. ! not in namelist IF (myproc == 0) THEN READ(5,wrfdfile) WRITE(6,'(a)') ' Namelist wrfdfile read in successfully.' strlen = LEN_TRIM(dir_extd) IF(strlen > 0) THEN IF(dir_extd(strlen:strlen) /= '/') THEN dir_extd(strlen+1:strlen+1) = '/' strlen = strlen + 1 END IF ELSE dir_extd = './' END IF IF (io_form > 100) THEN io_form = MOD(io_form,100) IF (mp_opt > 0) multifile = .TRUE. END IF END IF ! myproc == 0 CALL mpupdatec(dir_extd,256) CALL mpupdatei(io_form,1) CALL mpupdatel(multifile,1) CALL mpupdatec(init_time_str,19) CALL mpupdatec(start_time_str,19) CALL mpupdatec(end_time_str,19) CALL mpupdatec(history_interval,11) CALL mpupdatei(grid_id, 1) ! !----------------------------------------------------------------------- ! ! Read in namelist &wrfdfile ! !----------------------------------------------------------------------- ! nsnd=1 omegasnd=0.0 xpt=-999.0 ypt=-999.0 slat=-999.0 slon=-999.0 ipt=-99 jpt=-99 IF (myproc == 0 ) THEN READ(5,sndloc) IF( nsnd > maxsnd ) then WRITE(6,'(a,/a,i5)') & 'The number of sounding locations to be extracted exceeded maximum ',& 'allowed. nsnd is reset to ', maxsnd nsnd = maxsnd ENDIF DO isnd=1,nsnd xpt(isnd)=1000.*xpt(isnd) ypt(isnd)=1000.*ypt(isnd) END DO ENDIF CALL mpupdater(ugrid,1) CALL mpupdater(vgrid,1) CALL mpupdatei(locopt,1) CALL mpupdatei(nsnd,1) CALL mpupdater(xpt,maxsnd) CALL mpupdater(ypt,maxsnd) CALL mpupdater(slat,maxsnd) CALL mpupdater(slon,maxsnd) CALL mpupdatec(stid,8*maxsnd) CALL mpupdatei(istnm,maxsnd) IF (myproc == 0 ) THEN READ(5,output) ENDIF CALL mpupdatec(dir_output,256) CALL mpupdatei(outsnd,1) CALL mpupdatei(outsfc,1) IF (myproc == 0 ) THEN WRITE(6,'(4x,a,a)') ' Output directory: ',TRIM(dir_output) WRITE(6,'(2(4x,a,i3),/)') ' outsnd: ',outsnd,' outsfc: ',outsfc END IF !======================================================================= ! ! NAMELIST readings are done ! !======================================================================= rewindyr = .FALSE. READ(start_time_str,'(I4.4,5(a,I2.2))') & year,ach,month,ach,day,ach,hour,ach,minute,ach,second IF (year < 1960) THEN ! maybe ideal case myr = year year = 1960 rewindyr = .TRUE. END IF CALL ctim2abss(year,month,day,hour,minute,second,abstimes) READ(end_time_str,'(I4.4,a,I2.2,a,I2.2,a,I2.2,a,I2.2,a,I2.2)') & year,ach,month,ach,day,ach,hour,ach,minute,ach,second IF (rewindyr) year = 1960 CALL ctim2abss(year,month,day,hour,minute,second,abstimee) READ(history_interval,'(I2.2,a,I2.2,a,I2.2,a,I2.2)') & day,ach,hour,ach,minute,ach,second abstimei = day*24*3600+hour*3600+minute*60+second IF (multifile) THEN IF (MOD(nprocx_in,nproc_x) /= 0 .OR. MOD(nprocy_in,nproc_y) /= 0) THEN WRITE(6,*) 'nprocx_in (nprocy_in) must be dividable by nproc_x (nproc_y).' CALL arpsstop('WRONG message passing parameter.',1) END IF ncompressx = nprocx_in/nproc_x ncompressy = nprocy_in/nproc_y ELSE ! non-mpi or mpi with one file ncompressx = 1 ncompressy = 1 END IF ALLOCATE(fHndl(ncompressx,ncompressy), STAT = istatus) !----------------------------------------------------------------------- ! ! Check the availability of files and get parameter frames_per_outfile ! !----------------------------------------------------------------------- frames_per_outfile(:) = 1 nextdfil = 0 CALL check_wrf_files(multifile,MAXWRFFIL,grid_id,io_form,nprocx_in, & ncompressx,ncompressy,abstimes,abstimei,abstimee,rewindyr,myr, & dir_extd,extdname,nextdfil,frames_per_outfile,istatus) IF (istatus /= 0) CALL arpsstop('ERROR in check_wrf_files, See STDOUT for details',1) IF (mp_opt > 0) THEN ! should moved into mpinit_var later dumpstride = max_fopen readstride = nprocs IF (joindmp > 0) dumpstride = nprocs ! join and dump IF (multifile) readstride = max_fopen END IF IF (ANY(frames_per_outfile > 1) .AND. readstride < nprocs) THEN WRITE(6,'(3(a,/))') 'WARNING: WRF2ARPS does not support multi-frame in ',& ' one file for split-form WRF history files.', & ' The program is stopping ... ...' CALL arpsstop('frames_per_outfile should be 1.',1) END IF !----------------------------------------------------------------------- ! ! Get dimension parameters from the first input file ! !----------------------------------------------------------------------- ! blocking inserted for ordering i/o for message passing DO i=0,nprocs-1,readstride IF(myproc >= i .AND. myproc <= i+readstride-1) THEN CALL open_wrf_file(TRIM(extdname(1)),io_form,multifile,.TRUE., & ncompressx,ncompressy,fHndl) CALL get_wrf_metadata(fHndl,io_form,multifile,.TRUE.,1,1, & nx_wrf,ny_wrf,nz_wrf,nzsoil_wrf, & iproj_wrf,trlat1_wrf,trlat2_wrf,trlon_wrf, & ctrlat_wrf,ctrlon_wrf,dx_wrf,dy_wrf,dt_wrf, & sf_surface_physics,mp_physics,istatus) CALL close_wrf_file(fHndl,io_form,multifile,.TRUE.,ncompressx,ncompressy) END IF IF (mp_opt > 0) CALL mpbarrier END DO scale_wrf = 1.0 nstyp_wrf = 1 IF(myproc == 0) WRITE(6,'(/a,4(a,I4),/)') ' WRF grid dimensions:', & ' nx_wrf = ',nx_wrf, ' ny_wrf = ',ny_wrf, & ' nz_wrf = ',nz_wrf, ' nzsoil_wrf = ', nzsoil_wrf ! WRF forecast started time READ(init_time_str,'(I4.4,a,I2.2,a,I2.2,a,I2.2,a,I2.2,a,I2.2)') & year,ach,month,ach,day,ach,hour,ach,minute,ach,second CALL ctim2abss(year,month,day,hour,minute,second,initsec) thisdmp = abstimei tstart = 0.0 tstop = abstimee-initsec latitud = ctrlat !----------------------------------------------------------------------- ! ! Allocate and initalize arrays based on dimension parameters ! read in from the input file ! !----------------------------------------------------------------------- ! Note that for MP version nx & ny here are global values. They will ! be reassigned to their per-processor value below. IF (mp_opt > 0) THEN nx_wrf = (nx_wrf - 1)/nproc_x + 1 ! fake zone for WRF is 1 ny_wrf = (ny_wrf - 1)/nproc_y + 1 ELSE nproc_x = 1 nproc_y = 1 nprocs = 1 joindmp = 0 END IF ! !----------------------------------------------------------------------- ! ! Allocate and initialize external grid variables ! !----------------------------------------------------------------------- ! ALLOCATE(x_wrf(nx_wrf),stat=istatus) ALLOCATE(y_wrf(ny_wrf),stat=istatus) ALLOCATE(xu_wrf(nx_wrf),stat=istatus) ALLOCATE(yu_wrf(ny_wrf),stat=istatus) ALLOCATE(xv_wrf(nx_wrf),stat=istatus) ALLOCATE(yv_wrf(ny_wrf),stat=istatus) ALLOCATE(lat_wrf(nx_wrf,ny_wrf),stat=istatus) ALLOCATE(lon_wrf(nx_wrf,ny_wrf),stat=istatus) ALLOCATE(latu_wrf(nx_wrf,ny_wrf),stat=istatus) ALLOCATE(lonu_wrf(nx_wrf,ny_wrf),stat=istatus) ALLOCATE(latv_wrf(nx_wrf,ny_wrf),stat=istatus) ALLOCATE(lonv_wrf(nx_wrf,ny_wrf),stat=istatus) ALLOCATE(zp_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(hgt_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(zpsoil_wrf(nx_wrf,ny_wrf,nzsoil_wrf),stat=istatus) ALLOCATE(p_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(t_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(u_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(v_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(vatu_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(uatv_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(w_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(pt_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(qv_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(qc_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(qr_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(qi_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(qs_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(qh_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(psfc_wrf(nx_wrf,ny_wrf),stat=istatus) ALLOCATE(t2m_wrf(nx_wrf,ny_wrf),stat=istatus) ALLOCATE(th2m_wrf(nx_wrf,ny_wrf),stat=istatus) ALLOCATE(qv2m_wrf(nx_wrf,ny_wrf),stat=istatus) ALLOCATE(u10m_wrf(nx_wrf,ny_wrf),stat=istatus) ALLOCATE(v10m_wrf(nx_wrf,ny_wrf),stat=istatus) ALLOCATE(raing_wrf(nx_wrf,ny_wrf),stat=istatus) ALLOCATE(rainc_wrf(nx_wrf,ny_wrf),stat=istatus) ALLOCATE(tsoil_wrf (nx_wrf,ny_wrf,nzsoil_wrf,0:nstyp_wrf),stat=istatus) ALLOCATE(qsoil_wrf (nx_wrf,ny_wrf,nzsoil_wrf,0:nstyp_wrf),stat=istatus) ALLOCATE(wetcanp_wrf (nx_wrf,ny_wrf,0:nstyp_wrf),stat=istatus) ALLOCATE(soiltyp_wrf (nx_wrf,ny_wrf,nstyp_wrf),stat=istatus) ALLOCATE(stypfrct_wrf(nx_wrf,ny_wrf,nstyp_wrf),stat=istatus) ALLOCATE(vegtyp_wrf (nx_wrf,ny_wrf),stat=istatus) ALLOCATE(veg_wrf (nx_wrf,ny_wrf),stat=istatus) ALLOCATE(snowdpth_wrf(nx_wrf,ny_wrf),stat=istatus) ALLOCATE(trn_wrf (nx_wrf,ny_wrf),stat=istatus) ALLOCATE(tc(nz_wrf)) ALLOCATE(tdc(nz_wrf)) ALLOCATE(usnd(nz_wrf)) ALLOCATE(vsnd(nz_wrf)) ALLOCATE(dir(nz_wrf)) ALLOCATE(spd(nz_wrf)) x_wrf=0.0 y_wrf=0.0 xu_wrf=0.0 yu_wrf=0.0 xv_wrf=0.0 yv_wrf=0.0 lat_wrf=0.0 lon_wrf=0.0 latu_wrf=0.0 lonu_wrf=0.0 latv_wrf=0.0 lonv_wrf=0.0 zp_wrf=0.0 hgt_wrf=0.0 zpsoil_wrf=0.0 p_wrf=0.0 t_wrf=0.0 u_wrf=0.0 v_wrf=0.0 vatu_wrf=0.0 uatv_wrf=0.0 w_wrf=0.0 pt_wrf=0.0 qv_wrf=0.0 qc_wrf=0.0 qr_wrf=0.0 qi_wrf=0.0 qs_wrf=0.0 qh_wrf=0.0 raing_wrf = 0.0 rainc_wrf = 0.0 trn_wrf =0.0 tsoil_wrf =-999.0 qsoil_wrf =-999.0 wetcanp_wrf =-999.0 snowdpth_wrf=-999.0 soiltyp_wrf = -999 stypfrct_wrf= -999.0 vegtyp_wrf = -999 veg_wrf = -999.0 !----------------------------------------------------------------------- ! ! Allocate working arrays ! !----------------------------------------------------------------------- ALLOCATE(dxfld(nx_wrf),stat=istatus) ALLOCATE(dyfld(ny_wrf),stat=istatus) ALLOCATE(rdxfld(nx_wrf),stat=istatus) ALLOCATE(rdyfld(ny_wrf),stat=istatus) ALLOCATE(dxfldu(nx_wrf),stat=istatus) ALLOCATE(dyfldu(ny_wrf),stat=istatus) ALLOCATE(rdxfldu(nx_wrf),stat=istatus) ALLOCATE(rdyfldu(ny_wrf),stat=istatus) ALLOCATE(dxfldv(nx_wrf),stat=istatus) ALLOCATE(dyfldv(ny_wrf),stat=istatus) ALLOCATE(rdxfldv(nx_wrf),stat=istatus) ALLOCATE(rdyfldv(ny_wrf),stat=istatus) ALLOCATE(tem1_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(tem2_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(tem3_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(tem4_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(tem5_wrf(nx_wrf,ny_wrf,nz_wrf),stat=istatus) ALLOCATE(xa_wrf(nx_wrf,ny_wrf),stat=istatus) ALLOCATE(ya_wrf(nx_wrf,ny_wrf),stat=istatus) dxfld=0.0 dyfld=0.0 rdxfld=0.0 rdyfld=0.0 dxfldu=0.0 dyfldu=0.0 rdxfldu=0.0 rdyfldu=0.0 dxfldv=0.0 dyfldv=0.0 rdxfldv=0.0 rdyfldv=0.0 tem1_wrf=0.0 tem2_wrf=0.0 tem3_wrf=0.0 tem4_wrf=0.0 tem5_wrf=0.0 xa_wrf=0.0 ya_wrf=0.0 ! !----------------------------------------------------------------------- ! ! Loop through the data times provided via NAMELIST. ! !----------------------------------------------------------------------- ! iniotfu = 21 ! FORTRAN unit number used for data output first_time = .TRUE. wait_verif = .TRUE. in_verif = .FALSE. readstride = readstride/(ncompressx*ncompressy) IF (readstride < 1) THEN WRITE(6,*) 'ERROR: readstride < 1, please check max_fopen in namelist file.' WRITE(6,*) ' Please remember that readstride = max_fopen/(ncompressx*ncompressy)' CALL arpsstop('max_fopen too small',1) END IF exit_early = .FALSE. DO ifile = 1,nextdfil IF (myproc == 0) WRITE(6,'(2x,2a,/)') 'Reading file ',TRIM(extdname(ifile)) IF (frames_per_outfile(ifile) == 1) THEN ! Finish read in the following block ! ! blocking inserted for ordering i/o for message passing ! DO i=0,nprocs-1,readstride IF(myproc >= i .AND. myproc <= i+readstride-1) THEN CALL open_wrf_file(TRIM(extdname(ifile)),io_form,multifile, & .FALSE.,ncompressx,ncompressy,fHndl) !----------------------------------------------------------------------- ! ! Retrieve global attributes from the file ! !----------------------------------------------------------------------- CALL get_wrf_metadata(fHndl,io_form,multifile,.FALSE., & ncompressx,ncompressy, & nx_wrf,ny_wrf,nz_wrf,nzsoil_wrf, & iproj_wrf,trlat1_wrf,trlat2_wrf,trlon_wrf, & ctrlat_wrf,ctrlon_wrf,dx_wrf,dy_wrf,dt_wrf, & sf_surface_physics,mp_physics,istatus) scale_wrf = 1.0 latnot_wrf(1) = trlat1_wrf latnot_wrf(2) = trlat2_wrf IF (mp_opt > 0) THEN nx_wrf = (nx_wrf - 1)/nproc_x + 1 ny_wrf = (ny_wrf - 1)/nproc_y + 1 END IF itime = 1 CALL get_wrf_Times(fHndl,io_form,multifile,ncompressx,ncompressy, & itime,timestr) ! !----------------------------------------------------------------------- ! ! Call getwrfd to reads and converts data to ARPS units ! ! NOTE: u_wrf, v_wrf are just values extracted from data files. It may ! need to be rotated or extend to be MPI valid. ! !----------------------------------------------------------------------- ! CALL getwrfd(fHndl,io_form,multifile,ncompressx,ncompressy,itime, & timestr,nx_wrf,ny_wrf,nz_wrf,nzsoil_wrf, & iproj_wrf,scale_wrf,trlon_wrf,latnot_wrf, & ctrlat_wrf,ctrlon_wrf,dx_wrf,dy_wrf,x0_wrf,y0_wrf, & sf_surface_physics,mp_physics, & lat_wrf,lon_wrf,latu_wrf,lonu_wrf,latv_wrf,lonv_wrf, & zp_wrf,hgt_wrf,zpsoil_wrf, p_wrf,t_wrf, & u_wrf,v_wrf, w_wrf, & qv_wrf,qc_wrf,qr_wrf,qi_wrf,qs_wrf,qh_wrf, & psfc_wrf,t2m_wrf,th2m_wrf,qv2m_wrf,u10m_wrf,v10m_wrf,& tsoil_wrf(:,:,:,0),qsoil_wrf(:,:,:,0), & wetcanp_wrf(:,:,0),snowdpth_wrf,trn_wrf, & soiltyp_wrf(:,:,1),vegtyp_wrf,veg_wrf, & raing_wrf,rainc_wrf, & tem1_wrf,tem2_wrf,tem3_wrf,tem4_wrf,tem5_wrf, & istatus) CALL close_wrf_file(fHndl,io_form,multifile,.FALSE., & ncompressx,ncompressy) END IF IF (mp_opt > 0) CALL mpbarrier END DO ELSE ! frames_per_outfile > 1, only open the file and read meta CALL open_wrf_file(TRIM(extdname(ifile)),io_form,multifile, & .FALSE.,ncompressx,ncompressy,fHndl) !----------------------------------------------------------------------- ! ! Retrieve global attributes from the file ! !----------------------------------------------------------------------- CALL get_wrf_metadata(fHndl,io_form,multifile,.FALSE., & ncompressx,ncompressy, & nx_wrf,ny_wrf,nz_wrf,nzsoil_wrf, & iproj_wrf,trlat1_wrf,trlat2_wrf,trlon_wrf, & ctrlat_wrf,ctrlon_wrf,dx_wrf,dy_wrf,dt_wrf, & sf_surface_physics,mp_physics,istatus) scale_wrf = 1.0 latnot_wrf(1) = trlat1_wrf latnot_wrf(2) = trlat2_wrf IF (mp_opt > 0) THEN nx_wrf = (nx_wrf - 1)/nproc_x + 1 ny_wrf = (ny_wrf - 1)/nproc_y + 1 END IF END IF DO itime = 1, frames_per_outfile(ifile) IF (frames_per_outfile(ifile) > 1) THEN ! read data only when frames_per_outfile >1 ! otherwise, it has been read before. CALL get_wrf_Times(fHndl,io_form,multifile,ncompressx,ncompressy,& itime,timestr) ! !----------------------------------------------------------------------- ! ! Call getwrfds to reads and converts data to ARPS units ! ! NOTE: u_wrf, v_wrf are just values extracted from data files. It may ! need to be rotated or extend to be MPI valid. ! !----------------------------------------------------------------------- ! CALL getwrfd(fHndl,io_form,multifile,ncompressx,ncompressy,itime,& timestr,nx_wrf,ny_wrf,nz_wrf,nzsoil_wrf, & iproj_wrf,scale_wrf,trlon_wrf,latnot_wrf, & ctrlat_wrf,ctrlon_wrf,dx_wrf,dy_wrf,x0_wrf,y0_wrf, & sf_surface_physics,mp_physics, & lat_wrf,lon_wrf,latu_wrf,lonu_wrf,latv_wrf,lonv_wrf, & zp_wrf,hgt_wrf,zpsoil_wrf, p_wrf,t_wrf, & u_wrf,v_wrf, w_wrf, & qv_wrf,qc_wrf,qr_wrf,qi_wrf,qs_wrf,qh_wrf, & psfc_wrf,t2m_wrf,th2m_wrf,qv2m_wrf,u10m_wrf,v10m_wrf,& tsoil_wrf(:,:,:,0),qsoil_wrf(:,:,:,0), & wetcanp_wrf(:,:,0),snowdpth_wrf,trn_wrf, & soiltyp_wrf(:,:,1),vegtyp_wrf,veg_wrf, & raing_wrf,rainc_wrf, & tem1_wrf,tem2_wrf,tem3_wrf,tem4_wrf,tem5_wrf, & istatus) END IF !----------------------------------------------------------------------- ! ! Post-reading processing, most of this block is inside getwrfd before ! It must be moved here because of ordering I/O for message passing ! !----------------------------------------------------------------------- CALL adj_wrfuv(multifile,1,nx_wrf,ny_wrf,nz_wrf, & iproj_wrf,scale_wrf,trlon_wrf,latnot_wrf,x0_wrf,y0_wrf, & lonu_wrf,lonv_wrf,u_wrf,vatu_wrf,uatv_wrf,v_wrf, & tem1_wrf,tem2_wrf,istatus) pt_wrf(:,:,:) = t_wrf(:,:,:)*((p0/p_wrf(:,:,:))**rddcp) stypfrct_wrf(:,:,1) = 1. tsoil_wrf(:,:,:,1) = tsoil_wrf(:,:,:,0) qsoil_wrf(:,:,:,1) = qsoil_wrf(:,:,:,0) wetcanp_wrf(:,:,1) = wetcanp_wrf(:,:,0) ! !----------------------------------------------------------------------- ! ! Time conversions. ! Formats: timestr='1998-05-25_18:00:00 ! !----------------------------------------------------------------------- ! jtime=jtime+1 READ(timestr,'(I4.4,a1,I2.2,a1,I2.2,a1,I2.2,a1,I2.2,a1,I2.2)') & iyr,ach,imo,ach,iday,ach,ihr,ach,imin,ach,isec WRITE(timefil(jtime),'(i4.4,i2.2,i2.2,a1,3(i2.2))') & iyr,imo,iday,'-',ihr,imin,isec WRITE(timewrf(jtime),'(i4.4,i2.2,i2.2,1x,i2.2,2(a1,i2.2))') & iyr,imo,iday,ihr,':',imin,':',isec CALL julday(iyr,imo,iday,jldy) CALL ctim2abss(iyr,imo,iday,ihr,imin,isec,iabssec) IF(wait_verif) THEN IF(ihr == 6) THEN in_verif=.true. wait_verif=.false. CALL ctim2abss(iyr,imo,iday,6,0,0,istverif) END IF ELSE IF (in_verif) THEN IF((iabssec - istverif) > isec24h) in_verif=.false. END IF kftime=iabssec - initsec curtim=FLOAT(kftime) IF(myproc == 0) WRITE(6,'(/,a,a19,/,a,i12,a,i12,a,/)') & ' Subroutine getwrfd was called for data valid at ',timestr, & ' Which is ',iabssec,' abs seconds or ',kftime, & ' seconds from the WRF/ARPS initial time.' !----------------------------------------------------------------------- CALL a3dmax0(lat_wrf,1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,1,1,1,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'lat_wrf_min= ', amin,', lat_wrf_max=',amax CALL a3dmax0(lon_wrf,1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,1,1,1,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'lon_wrf_min= ', amin,', lon_wrf_max=',amax CALL a3dmax0(p_wrf ,1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,nz_wrf,1,nz_wrf,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'p_wrf_min = ', amin,', p_wrf_max =',amax CALL a3dmax0(hgt_wrf,1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,nz_wrf,1,nz_wrf,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'hgt_wrf_min= ', amin,', hgt_wrf_max=',amax CALL a3dmax0(t_wrf ,1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,nz_wrf,1,nz_wrf,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 't_wrf_min = ', amin,', t_wrf_max =',amax CALL a3dmax0(u_wrf ,1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,nz_wrf,1,nz_wrf,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'u_wrf_min = ', amin,', u_wrf_max =',amax CALL a3dmax0(v_wrf ,1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,nz_wrf,1,nz_wrf,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'v_wrf_min = ', amin,', v_wrf_max =',amax CALL a3dmax0(w_wrf ,1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,nz_wrf,1,nz_wrf,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'w_wrf_min = ', amin,', w_wrf_max =',amax CALL a3dmax0(qv_wrf ,1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,nz_wrf,1,nz_wrf,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'qv_wrf_min = ', amin,', qv_wrf_max =',amax CALL a3dmax0(qc_wrf ,1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,nz_wrf,1,nz_wrf,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'qc_wrf_min = ', amin,', qc_wrf_max =',amax CALL a3dmax0(qr_wrf ,1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,nz_wrf,1,nz_wrf,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'qr_wrf_min = ', amin,', qr_wrf_max =',amax CALL a3dmax0(qi_wrf ,1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,nz_wrf,1,nz_wrf,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'qi_wrf_min = ', amin,', qi_wrf_max =',amax CALL a3dmax0(qs_wrf ,1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,nz_wrf,1,nz_wrf,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'qs_wrf_min = ', amin,', qs_wrf_max =',amax CALL a3dmax0(qh_wrf ,1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,nz_wrf,1,nz_wrf,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'qh_wrf_min = ', amin,', qh_wrf_max =',amax CALL a3dmax0(tsoil_wrf(1,1,1,0),1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,1,1,1,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'tsoil_sfc_wrf_min= ', amin,', tsoil_sfc_wrf_max=',amax CALL a3dmax0(tsoil_wrf(1,1,2,0),1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,1,1,1,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'tsoil_1st_wrf_min= ', amin,', tsoil_1st_wrf_max=',amax CALL a3dmax0(qsoil_wrf(1,1,1,0),1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,1,1,1,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'qsoil_sfc_wrf_min= ', amin,', qsoil_sfc_wrf_max=',amax CALL a3dmax0(qsoil_wrf(1,1,2,0),1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,1,1,1,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'qsoil_1st_wrf_min= ', amin,', qsoil_1st_wrf_max=',amax CALL a3dmax0(wetcanp_wrf,1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,1,1,1,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'wetcanp_wrf_min= ', amin,', wetcanp_wrf_max=',amax CALL a3dmax0(snowdpth_wrf,1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,1,1,1,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'snowd_wrf_min= ', amin,', snow_wrf_max=',amax CALL a3dmax0(trn_wrf ,1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,1,1,1,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'trn_wrf_min= ', amin,', trn_wrf_max=',amax CALL a3dmax0(veg_wrf,1,nx_wrf,1,nx_wrf,1,ny_wrf,1,ny_wrf, & 1,1,1,1,amax,amin) IF(myproc == 0) WRITE(6,'(1x,2(a,e13.6))') & 'veg_wrf_min= ', amin,', veg_wrf_max=',amax IF(myproc == 0) PRINT*,' ' ! !----------------------------------------------------------------------- ! ! First time through the time loop, calculate grid ! transformation info. ! !----------------------------------------------------------------------- ! IF(first_time) THEN ! !----------------------------------------------------------------------- ! ! Find x,y locations of WRF grid. ! !----------------------------------------------------------------------- ! CALL setmapr(iproj_wrf,scale_wrf,latnot_wrf,trlon_wrf) CALL setorig(1,x0_wrf,y0_wrf) DO j=1,ny_wrf CALL lltoxy(1,1,lat_wrf(1,j),lon_wrf(1,j), & x_wrf(1),y_wrf(j)) END DO DO i=1,nx_wrf CALL lltoxy(1,1,lat_wrf(i,1),lon_wrf(i,1), & x_wrf(i),y_wrf(1)) END DO DO j=1,ny_wrf CALL lltoxy(1,1,latu_wrf(1,j),lonu_wrf(1,j), & xu_wrf(1),yu_wrf(j)) END DO DO i=1,nx_wrf CALL lltoxy(1,1,latu_wrf(i,1),lonu_wrf(i,1), & xu_wrf(i),yu_wrf(1)) END DO DO j=1,ny_wrf CALL lltoxy(1,1,latv_wrf(1,j),lonv_wrf(1,j), & xv_wrf(1),yv_wrf(j)) END DO DO i=1,nx_wrf CALL lltoxy(1,1,latv_wrf(i,1),lonv_wrf(i,1), & xv_wrf(i),yv_wrf(1)) END DO ! !----------------------------------------------------------------------- ! ! Find x,y of sounding locations in external grid. ! !----------------------------------------------------------------------- ! DO isnd=1,nsnd IF(locopt == 1) THEN CALL lltoxy(1,1,slat(isnd),slon(isnd),xpt(isnd),ypt(isnd)) END IF IF (xpt(isnd) > x_wrf(1) .AND. xpt(isnd) < x_wrf(nx_wrf) .AND. & ypt(isnd) > y_wrf(1) .AND. ypt(isnd) < y_wrf(ny_wrf) ) THEN ptingrid(isnd)=.TRUE. dmin=((xpt(isnd)-x_wrf(1))*(xpt(isnd)-x_wrf(1))+ & (ypt(isnd)-y_wrf(1))*(ypt(isnd)-y_wrf(1))) ipt(isnd)=1 jpt(isnd)=1 DO j=1,ny_wrf-1 DO i=1,nx_wrf-1 dd=((xpt(isnd)-x_wrf(i))*(xpt(isnd)-x_wrf(i))+ & (ypt(isnd)-y_wrf(j))*(ypt(isnd)-y_wrf(j))) IF(dd < dmin) THEN dmin=dd ipt(isnd)=i jpt(isnd)=j END IF END DO END DO WRITE(6,'(a,f10.4,f10.4,/a,i5,i5,a,f10.4,f10.4)') & ' Nearest WRF pt to diagnostic lat,lon: ', & slat(isnd),slon(isnd), & ' Diagnostic i,j: ', ipt(isnd),jpt(isnd), & ' lat,lon= ', lat_wrf(ipt(isnd),jpt(isnd)), & lon_wrf(ipt(isnd),jpt(isnd)) WRITE(6,'(1x,a,f10.2,a)') 'Terrain at WRF pt: ', & trn_wrf(ipt(isnd),jpt(isnd)),' meters.' ELSE ptingrid(isnd)=.FALSE. ipt(isnd)=1 jpt(isnd)=1 ! WRITE(6,'(1x,a,I2,a)') 'ERROR: Station ',isnd, & ! ' location is outside of the WRF domain.' ! CALL arpsstop('Station location exceeds the WRF domain',1) END IF END DO ! isnd first_time = .FALSE. END IF ! first_time DO isnd=1,nsnd IF(ptingrid(isnd)) THEN iloc=ipt(isnd) jloc=jpt(isnd) WRITE(6,'(///2(a,I4),/2x,a)') & ' External sounding at iloc,jloc = ',iloc,', ',jloc, & 'k pres hgt temp theta dewp u v dir spd' ! !----------------------------------------------------------------------- ! ! Output sounding to look like GEMPAK SNLIST.FIL ! to use the Skew-T and Hodograph programs ! e.g., those by Rich Carpenter and NSHARP ! ! Example of what the header looks like: ! !23456789012345678901234567890123456789012345678901234567890 !SNPARM = PRES;HGHT;TMPC;DWPC;DRCT;SPED;OMEG ! !STID = SEP STNM = 72260 TIME = 920308/1500 !SLAT = 32.21 SLON = -98.18 SELEV = 399.0 ! ! PRES HGHT TMPC DWPC DRCT SPED OMEG ! !----------------------------------------------------------------------- ! IF(outsnd > 0) THEN WRITE(uafile,'(6a)') TRIM(dir_output),'/',TRIM(stid(isnd)),& '-',timefil(jtime),'.gemsnd' OPEN(31,FILE=TRIM(uafile),STATUS='unknown') modyr=mod(iyr,100) WRITE(31,'(/a/)') ' SNPARM = PRES;HGHT;TMPC;DWPC;DRCT;SPED;OMEG' WRITE(31,'(a,a8,3x,a,i8,3x,a,i2.2,i2.2,i2.2,a1,i2.2,i2.2)') & ' STID = ',stid(isnd), & 'STNM = ',istnm(isnd), & 'TIME = ',modyr,imo,iday,'/',ihr,imin WRITE(31,'(a,f8.2,3x,a,f8.2,3x,a,f7.1/)') & ' SLAT = ',slat(isnd),'SLON = ',slon(isnd), & 'SELV = ',selev(isnd) WRITE(31,'(6x,a)') & 'PRES HGHT TMPC DWPC DRCT SPED OMEG' END IF ! ! Convert units of external data and write as a sounding. ! DO k=1,nz_wrf pmb=.01*p_wrf(iloc,jloc,k) tc(k)=t_wrf(iloc,jloc,k)-273.15 theta=pt_wrf(iloc,jloc,k) IF( qv_wrf(iloc,jloc,k) > 0.) THEN smix=qv_wrf(iloc,jloc,k)/(1.-qv_wrf(iloc,jloc,k)) e=(pmb*smix)/(0.62197 + smix) bige=e/( 1.001 + ( (pmb - 100.) / 900.) * 0.0034) alge = ALOG(bige/6.112) tdc(k) = (alge * 243.5) / (17.67 - alge) ELSE tdc(k) = tc(k)-30. END IF usnd=0.5*(u_wrf(iloc,jloc,k)+ & u_wrf((iloc+1),jloc,k)) vsnd=0.5*(v_wrf(iloc,jloc,k)+ & v_wrf(iloc,(jloc+1),k)) CALL uvrotdd(1,1,lon_wrf(iloc,jloc), & usnd(k),vsnd(k),dir(k),spd(k)) WRITE(6,'(i4,f6.0,f9.0,f7.1,f7.1,f7.1,f7.1,f7.1,f7.1,f7.1)') & k,pmb, & hgt_wrf(iloc,jloc,k), & tc(k),theta,tdc(k), & usnd(k),vsnd(k),dir(k),spd(k) IF(outsnd > 0 ) WRITE(31,'(1x,7(F9.2))') & pmb,hgt_wrf(iloc,jloc,k), & tc(k),tdc(k),dir(k),spd(k),omegasnd END DO IF(outsnd > 0) CLOSE(31) ELSE ! pt not in grid ! ! Create files with an error message ! IF(outsnd > 0) THEN WRITE(uafile,'(6a)') TRIM(dir_output),'/',TRIM(stid(isnd)),& '-',timefil(jtime),'.gemsnd' OPEN(31,FILE=TRIM(uafile),STATUS='unknown') WRITE(31,'(3a)') ' Station ',stid(isnd),' not in model domain' CLOSE(31) END IF END IF ! scrhgt=selev(isnd)+2.0 annhgt=selev(isnd)+10.0 tcsoil=tsoil_wrf(iloc,jloc,1,0)-273.15 IF(scrhgt > trn_wrf(iloc,jloc)) THEN DO k=1,nz_wrf IF(hgt_wrf(iloc,jloc,k) > scrhgt) THEN ktop=k EXIT END IF END DO IF(k > 1) THEN kbot=k-1 delt=tc(ktop)-tc(kbot) delz=hgt_wrf(iloc,jloc,ktop)-hgt_wrf(iloc,jloc,kbot) tsfc=tc(kbot)+((scrhgt-hgt_wrf(iloc,jloc,kbot))*(delt/delz)) ELSE delt=tc(1)-tcsoil delz=hgt_wrf(iloc,jloc,1)-trn_wrf(iloc,jloc) tsfc=tcsoil+((scrhgt-trn_wrf(iloc,jloc))*(delt/delz)) END IF ELSE ! below terrain tsfc=tc(1)+0.0065*(scrhgt-hgt_wrf(iloc,jloc,1)) END IF print *, ' terrain, scrhgt, hgt(1):', & trn_wrf(iloc,jloc),scrhgt,hgt_wrf(iloc,jloc,1) print *, ' tcsoil,tsfc,tc(1):',tcsoil,tsfc,tc(1) tf=tctof(tsfc) tdf=tctof(tdc(1)) spdkts=spd(1)*ms2kts raing=raing_wrf(iloc,jloc)*mm2inch rainc=rainc_wrf(iloc,jloc)*mm2inch raint=raing+rainc ! WRITE(isfcunit(isnd),'(1x,a8,1x,i4.4,2i2.2,1x,i2.2,2(a,i2.2),4f6.0,3f8.3)') & ! stid(isnd),iyr,imo,iday,ihr,':',imin,':',isec, & ! tf,tdf,dir(1),spdkts,raing,rainc,raint pmb=0.01*psfc_wrf(iloc,jloc) tf=tktof(t2m_wrf(iloc,jloc)) IF( qv2m_wrf(iloc,jloc) > 0.) THEN smix=qv2m_wrf(iloc,jloc)/(1.-qv2m_wrf(iloc,jloc)) e=(pmb*smix)/(0.62197 + smix) bige=e/( 1.001 + ( (pmb - 100.) / 900.) * 0.0034) alge = ALOG(bige/6.112) tdsfc = (alge * 243.5) / (17.67 - alge) ELSE tdsfc = t2m_wrf(iloc,jloc)-303.15 END IF tdf=tctof(tdsfc) CALL uvrotdd(1,1,lon_wrf(iloc,jloc), & u10m_wrf(iloc,jloc),v10m_wrf(iloc,jloc),dirsfc,spdsfc) spdkts=spdsfc*ms2kts sfcpr(jtime,isnd)=pmb tsfcf(jtime,isnd)=tf tdsfcf(jtime,isnd)=tdf wdir(jtime,isnd)=dirsfc wspd(jtime,isnd)=spdkts prcgp(jtime,isnd)=raing prccu(jtime,isnd)=rainc prctot(jtime,isnd)=raint IF(in_verif) THEN tmax(isnd)=max(tmax(isnd),tf) tmin(isnd)=min(tmin(isnd),tf) spdmax(isnd)=max(spdmax(isnd),spdkts) IF(prc0(isnd) < 0.) prc0(isnd)=raint prcsum(isnd)=raint-prc0(isnd) END IF END DO ! isnd IF (iabssec >= abstimee) THEN ! exit if time exceeds required in namelist exit_early = .TRUE. EXIT END IF END DO ! itime IF (frames_per_outfile(ifile) > 1) THEN CALL close_wrf_file(fHndl,io_form,multifile,.FALSE., & ncompressx,ncompressy) END IF IF (exit_early) EXIT END DO ! ifile ntime=jtime IF(outsfc > 0) THEN DO isnd=1, nsnd WRITE(sfcfile,'(6a)') TRIM(dir_output),'/', & TRIM(stid(isnd)),'-',timefil(1),'.sfc' OPEN(32,FILE=TRIM(sfcfile),STATUS='unknown') IF(ptingrid(isnd)) THEN WRITE(32,'(2a,2(a,f7.0))') & ' Station: ',stid(isnd),' Elev(m): ',selev(isnd), & ' WRF Terrain(m): ', trn_wrf(ipt(isnd),jpt(isnd)) WRITE(32,'(a,a)') ' Date Time(UTC) SfcPr Tsfc Tdsfc', & ' WDir WSpd PrecGP PrecCu PrecTot' DO jtime=1,ntime WRITE(32,'(1x,a,5f6.0,3f8.3)') timewrf(jtime), & sfcpr(jtime,isnd),tsfcf(jtime,isnd),tdsfcf(jtime,isnd), & wdir(jtime,isnd),wspd(jtime,isnd),prcgp(jtime,isnd), & prccu(jtime,isnd),prctot(jtime,isnd) END DO WRITE(32,'(a,f6.0,a,f6.0,a,f6.0,a,f8.3)') & ' Tmin=',tmin(isnd),' Tmax=',tmax(isnd),& ' WSpd Max=',spdmax(isnd),' Precip=',prcsum(isnd) ELSE WRITE(6,'(3a)') ' Station: ',stid(isnd),' is outside WRF domain' WRITE(32,'(3a)') ' Station: ',stid(isnd),' is outside WRF domain' END IF CLOSE(32) END DO END IF CALL io_shutdown(io_form) ! !----------------------------------------------------------------------- ! ! Friendly exit message. ! !----------------------------------------------------------------------- ! IF(myproc == 0) WRITE(6,'(a/a,i4,a)') & ' ==== Normal succesful completion of WRFEXTSND ====', & ' Processed',nextdfil,' file(s)' IF (mp_opt > 0) CALL mpexit(0) STOP ! !----------------------------------------------------------------------- ! ! Error status returned from getwrfd ! !----------------------------------------------------------------------- ! IF(myproc == 0) WRITE(6,'(a,i6)') & ' Aborting, error reading external file. istatus=',istatus IF (mp_opt > 0) CALL mpexit(0) STOP END PROGRAM wrfextsnd FUNCTION tctof(tc) IMPLICIT NONE REAL :: tctof REAL :: tc REAL, PARAMETER :: fconst = 32.0 REAL, PARAMETER :: fcratio = (9./5.) tctof = (tc*fcratio) + fconst RETURN END FUNCTION tktof(tk) IMPLICIT NONE REAL :: tktof REAL :: tk REAL, PARAMETER :: fconst = 32.0 REAL, PARAMETER :: fcratio = (9./5.) tktof = ((tk-273.15)*fcratio) + fconst RETURN END