subroutine convert_binary_mass !$$$ subprogram documentation block ! . . . . ! subprogram: convert_binary_mass read wrf mass binary restart ! prgmmr: parrish org: np22 date: 2003-09-05 ! ! abstract: using wrf library routines, read a wrf mass core binary ! format restart file. write the result to temporary binary ! file expected by read_wrf_mass_guess. ! ! program history log: ! 2004-09-10 parrish ! 2004-10-26 wu - date info from the field (T) instead of START_DATE, ! current time is necessary for cycling. ! 2004-11-05 wu - add check on input wrf guess file, stop code if problem ! 2004-11-07 parrish - change so byte offset information is written, instead ! of a whole new binary file--this done so mpi-io ! can be used later to read/write directly from ! the wrf mass core binary file. also, do inventory ! of whole file, so offsets are general--not dependent ! on type of file. ! 2004-12-15 treadon - remove get_lun, read guess from file "wrf_inout" ! 2005-07-06 parrish - add read of pint byte address ! 2005-11-29 parrish - add changes to allow earlier reading of surface fields needed ! for intelligent thinning of satellite data. ! 2006-04-06 middlecoff - changed out_unit from 55 to lendian_out ! 2006-09-15 treadon - use nhr_assimilation to build local guess filename ! 2007-04-12 parrish - add modifications to allow any combination of ikj or ijk ! grid ordering for input 3D fields ! ! input argument list: ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ ! COMMENTS ON MPI I-O AS USED TO READ AND WRITE FROM WRF BINARY FORMAT FILES: ! The wrf binary format model files are sequential unformatted files. There is a sequence ! of 2048 byte informational records, which are sometimes followed by a record with a ! variable specified by an identifier in the preceding informational record. ! The first step, carried out in this routine, convert_binary_mass (or convert_binary_nmm), ! is to do an inventory, by reading the file opened for direct access. Sequential files ! have a 4 leading and trailing bytes for each record, containing the length of the record ! in bytes. These are first counted by subroutine count_recs_wrf_binary_file to get ! the total number of records on the file, and then an inventory is created with ! subroutine inventory_wrf_binary_file. This creates a list of field identifiers for ! everything contained on the file, and byte offsets to the beginning of each field ! relative to the start of the file (taking into account the 4-byte leading and trailing ! record markers also). This inventory information is then scanned to select the ! variables required by the analysis, and the corresponding byte offsets are written to ! a temporary file for later access by subroutines read_wrf_mass_binary_guess, ! wrwrf_mass_binary, etc., all of which read and write to the wrf binary file using ! mpi i-o routines. ! This initial inventory is very fast because of the direct access reads, with record length ! set to the optimal value for the disks. Also, not all direct access records are read, but ! only those which contain sequential record marks. ! ! The mpi i-o routines as used here are quite simple. Taking the wrf file as a contiguous ! string of bytes, any piece of the file, characterized by a byte offset from the start of ! the file, and a length, can be read or written to/from any processor, as long as one is ! careful not to write into the same area from different processors at the same time. ! ! The routines used here are ! mpi_file_read_at and mpi_file_write_at ! which are blocking commands, ie, after the call, the buffer used to transfer bytes to/from ! memory is free to be used. See subroutines read_wrf_mass_binary_guess, read_wrf_nmm_binary_guess, ! wrwrfnmma_binary, and wrwrfmassa_binary. ! These routines have been made a bit more complicated due to the ordering of 3-d variables ! on the wrf files. They are currently stored in ikj format, but the analysis requires reading ! and writing in ijk format. To get around this in an efficient manner, 3-d variables are ! read in in ik groups for a range of j-index on each processor, where the j-index is evenly ! distributed across all processors. Once the read is complete, an mpi_alltoallv command is ! used to redistribute so the k-index is distributed evenly across processors, and each processor ! has a subset of full ij fields over its designated k-index range. use kinds, only: r_single,i_llong,i_kind use gsi_4dvar, only: nhr_assimilation use gsi_io, only: lendian_out implicit none ! Declare local parameters integer(i_kind),parameter:: in_unit = 15 real(r_single),parameter:: one_single = 1.0_r_single real(r_single),parameter:: r45 = 45.0_r_single character(6) filename character(9) wrfges integer(i_kind),allocatable:: start_block(:),end_block(:) integer(i_kind),allocatable:: start_byte(:),end_byte(:) integer(i_llong),allocatable:: file_offset(:) integer(i_llong) n_position character(132),allocatable:: datestr_all(:),varname_all(:),memoryorder_all(:) integer(i_kind),allocatable:: domainend_all(:,:) integer(i_kind) nrecs integer(i_kind) status_hdr integer(i_kind) hdrbuf(512) integer(i_kind) iyear,imonth,iday,ihour,iminute,isecond integer(i_kind) nlon_regional,nlat_regional,nsig_regional real(r_single) pt_regional integer(i_kind) k real(r_single),allocatable::field1(:),field1p(:),field2(:,:),field2b(:,:),field2c(:,:) real(r_single) rad2deg_single real(r_single)rdx,rdy integer(i_kind) ksize integer(i_kind) index wrfges = 'wrf_inout' open(in_unit,file=wrfges,form='unformatted') write(6,*)' convert_binary_mass: in_unit,lendian_out=',in_unit,lendian_out write(filename,100) nhr_assimilation 100 format('sigf',i2.2) open(lendian_out,file=filename,form='unformatted') rewind lendian_out ! Check for valid input file read(in_unit,iostat=status_hdr)hdrbuf if(status_hdr /= 0) then write(6,*)'CONVERT_BINARY_MASS: problem with wrfges = ',& trim(wrfges),', Status = ',status_hdr call stop2(74) endif ! reopen for direct access reading of sequential file close(in_unit) call count_recs_wrf_binary_file(in_unit,wrfges,nrecs) write(6,*) ' after count_recs_wrf_binary_file, nrecs=',nrecs allocate(datestr_all(nrecs),varname_all(nrecs),domainend_all(3,nrecs)) allocate(memoryorder_all(nrecs)) allocate(start_block(nrecs),end_block(nrecs),start_byte(nrecs),end_byte(nrecs),file_offset(nrecs)) call inventory_wrf_binary_file(in_unit,wrfges,nrecs, & datestr_all,varname_all,memoryorder_all,domainend_all, & start_block,end_block,start_byte,end_byte,file_offset) ! start with date record for date forecast was started ! y,m,d,h,m,s call retrieve_index(index,'START_DATE',varname_all,nrecs) if(index<0) stop read(datestr_all(index),'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') & iyear,imonth,iday,ihour,iminute,isecond write(6,*)' convert_binary_mass: START_DATE =',& iyear,imonth,iday,ihour,iminute,isecond ! nlon_regional, nlat_regional, nsig_regional call retrieve_index(index,'T',varname_all,nrecs) if(index<0) stop if(trim(memoryorder_all(index))=='XZY') then nlon_regional=domainend_all(1,index) nlat_regional=domainend_all(3,index) nsig_regional=domainend_all(2,index) end if if(trim(memoryorder_all(index))=='XYZ') then nlon_regional=domainend_all(1,index) nlat_regional=domainend_all(2,index) nsig_regional=domainend_all(3,index) end if read(datestr_all(index),'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') & iyear,imonth,iday,ihour,iminute,isecond write(6,*)' convert_binary_mass: iy,m,d,h,m,s=',& iyear,imonth,iday,ihour,iminute,isecond write(6,*)' convert_binary_mass: nlon,lat,sig_regional=',& nlon_regional,nlat_regional,nsig_regional ! pt_regional call retrieve_index(index,'P_TOP',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,pt_regional,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_mass: pt_regional=',pt_regional write(lendian_out) iyear,imonth,iday,ihour,iminute,isecond, & nlon_regional,nlat_regional,nsig_regional,pt_regional allocate(field1(nsig_regional),field1p(nsig_regional+1)) ! znu call retrieve_index(index,'ZNU',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,field1,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) do k=1,nsig_regional write(6,*)' convert_binary_mass: k,znu(k)=',k,field1(k) end do write(lendian_out)field1 ! ZNU ! znw call retrieve_index(index,'ZNW',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,field1p,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) do k=1,nsig_regional+1 write(6,*)' convert_binary_mass: k,znw(k)=',k,field1p(k) end do write(lendian_out)field1p ! ZNW deallocate(field1,field1p) ! rdx call retrieve_index(index,'RDX',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,rdx,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_mass: 1/rdx=',& one_single/rdx ! 1._4 necessary only to get bit reproducibility ! rdy call retrieve_index(index,'RDY',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,rdy,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_mass: 1/rdy=',& one_single/rdy ! 1._4 necessary only to get bit reproducibility allocate(field2(nlon_regional,nlat_regional)) allocate(field2b(nlon_regional,nlat_regional)) allocate(field2c(nlon_regional,nlat_regional)) ! mapfac_m call retrieve_index(index,'MAPFAC_M',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,field2,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_mass: max,min mapfac_m=',maxval(field2),minval(field2) write(6,*)' convert_binary_mass: max,min MAPFAC_M(:,1)=', & maxval(field2(:,1)),minval(field2(:,1)) write(6,*)' convert_binary_mass: max,min MAPFAC_M(1,:)=', & maxval(field2(1,:)),minval(field2(1,:)) write(6,*)' convert_binary_mass: mapfac_m(1,1),mapfac_m(nlon,1)=', & field2(1,1),field2(nlon_regional,1) write(6,*)' convert_binary_mass: mapfac_m(1,nlat),mapfac_m(nlon,nlat)=', & field2(1,nlat_regional),field2(nlon_regional,nlat_regional) field2b=one_single/(field2*rdx) !DX_MC field2c=one_single/(field2*rdy) !DY_MC ! XLAT rad2deg_single=r45/atan(one_single) call retrieve_index(index,'XLAT',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,field2,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_mass: max,min XLAT(:,1)=',& maxval(field2(:,1)),minval(field2(:,1)) write(6,*)' convert_binary_mass: max,min XLAT(1,:)=',& maxval(field2(1,:)),minval(field2(1,:)) write(6,*)' convert_binary_mass: xlat(1,1),xlat(nlon,1)=',& field2(1,1),field2(nlon_regional,1) write(6,*)' convert_binary_mass: xlat(1,nlat),xlat(nlon,nlat)=', & field2(1,nlat_regional),field2(nlon_regional,nlat_regional) field2=field2/rad2deg_single write(lendian_out)field2,field2b ! XLAT,DX_MC ! XLONG call retrieve_index(index,'XLONG',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,field2,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_mass: max,min XLONG(:,1)=',& maxval(field2(:,1)),minval(field2(:,1)) write(6,*)' convert_binary_mass: max,min XLONG(1,:)=',& maxval(field2(1,:)),minval(field2(1,:)) write(6,*)' convert_binary_mass: xlong(1,1),xlong(nlon,1)=',& field2(1,1),field2(nlon_regional,1) write(6,*)' convert_binary_mass: xlong(1,nlat),xlong(nlon,nlat)=', & field2(1,nlat_regional),field2(nlon_regional,nlat_regional) field2=field2/rad2deg_single write(lendian_out)field2,field2c ! XLONG,DY_MC write(lendian_out) wrfges ! index for START_DATE record call retrieve_index(index,'START_DATE',varname_all,nrecs) if(index<0) stop n_position=file_offset(index) write(lendian_out)n_position ! offset for START_DATE record ! MUB call retrieve_index(index,'MUB',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for MUB = ',n_position write(lendian_out)n_position ! offset for mub ! MU call retrieve_index(index,'MU',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for MU = ',n_position write(lendian_out)n_position ! offset for mu ! PHB call retrieve_index(index,'PHB',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset, memoryorder for PHB(',k,') = ',n_position,memoryorder_all(index) write(lendian_out)n_position,memoryorder_all(index) ! offset for PHB ! (zsfc*g is 1st level of this 3d field) ! but more efficient to read in whole 3-d field because of ikj order ! T call retrieve_index(index,'T',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset, memoryorder for T(',k,') = ',n_position,memoryorder_all(index) write(lendian_out)n_position,memoryorder_all(index) ! offset for T(k) ! POT TEMP (sensible) ! QVAPOR call retrieve_index(index,'QVAPOR',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset, memoryorder for QVAPOR(',k,' = ',n_position,memoryorder_all(index) write(lendian_out)n_position,memoryorder_all(index) ! offset for QVAPOR(k) ! U call retrieve_index(index,'U',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset, memoryorder for U(',k,' = ',n_position,memoryorder_all(index) write(lendian_out)n_position,memoryorder_all(index) ! offset for U(k) ! V call retrieve_index(index,'V',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset, memoryorder for V(',k,' = ',n_position,memoryorder_all(index) write(lendian_out)n_position,memoryorder_all(index) ! offset for V(k) ! LANDMASK call retrieve_index(index,'LANDMASK',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for LANDMASK = ',n_position write(lendian_out)n_position ! LANDMASK (1=land, 0=water) ! XICE call retrieve_index(index,'XICE',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for XICE = ',n_position write(lendian_out)n_position ! XICE ! SST call retrieve_index(index,'SST',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for SST = ',n_position write(lendian_out)n_position ! SST ! IVGTYP call retrieve_index(index,'IVGTYP',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for IVGTYP = ',n_position write(lendian_out)n_position ! IVGTYP ! ISLTYP call retrieve_index(index,'ISLTYP',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for ISLTYP = ',n_position write(lendian_out)n_position ! ISLTYP ! VEGFRA call retrieve_index(index,'VEGFRA',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for VEGFRA = ',n_position write(lendian_out)n_position ! VEGFRA ! SNOW call retrieve_index(index,'SNOW',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for SNOW = ',n_position write(lendian_out)n_position ! SNOW ! U10 call retrieve_index(index,'U10',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for U10 = ',n_position write(lendian_out)n_position ! U10 ! V10 call retrieve_index(index,'V10',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for V10 = ',n_position write(lendian_out)n_position ! V10 ! SMOIS call retrieve_index(index,'SMOIS',varname_all,nrecs) if(index<0) stop if(trim(memoryorder_all(index))=='XZY') then ksize=domainend_all(2,index) end if if(trim(memoryorder_all(index))=='XYZ') then ksize=domainend_all(3,index) end if n_position=file_offset(index+1) write(6,*)' byte offset, ksize, memoryorder for SMOIS(',k,') = ', & n_position,ksize,memoryorder_all(index) write(lendian_out)n_position,ksize,memoryorder_all(index) ! SMOIS ! TSLB call retrieve_index(index,'TSLB',varname_all,nrecs) if(index<0) stop if(trim(memoryorder_all(index))=='XZY') then ksize=domainend_all(2,index) end if if(trim(memoryorder_all(index))=='XYZ') then ksize=domainend_all(3,index) end if n_position=file_offset(index+1) write(6,*)' byte offset, ksize, memoryorder for TSLB(',k,') = ', & n_position,ksize,memoryorder_all(index) write(lendian_out)n_position,ksize,memoryorder_all(index) ! TSLB ! TSK call retrieve_index(index,'TSK',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for TSK = ',n_position write(lendian_out)n_position ! TSK !??????????????????/later put in z0 here, but for now just fill with something call retrieve_index(index,'TSK',varname_all,nrecs) call retrieve_field(in_unit,wrfges,field2,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_mass: max,min Z0=', & maxval(field2),minval(field2) write(lendian_out)field2 ! Z0 call retrieve_index(index,'SST',varname_all,nrecs) call retrieve_field(in_unit,wrfges,field2,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_mass: max,min SST=', & maxval(field2),minval(field2) write(lendian_out)field2 ! SST call retrieve_index(index,'TSK',varname_all,nrecs) call retrieve_field(in_unit,wrfges,field2,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_mass: max,min TSK=', & maxval(field2),minval(field2) write(lendian_out)field2 ! TSK call retrieve_index(index,'LANDMASK',varname_all,nrecs) call retrieve_field(in_unit,wrfges,field2,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_mass: max,min LANDMASK=', & maxval(field2),minval(field2) write(lendian_out)field2 ! LANDMASK call retrieve_index(index,'XICE',varname_all,nrecs) call retrieve_field(in_unit,wrfges,field2,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_mass: max,min XICE=', & maxval(field2),minval(field2) write(lendian_out)field2 ! XICE call retrieve_index(index,'SNOW',varname_all,nrecs) call retrieve_field(in_unit,wrfges,field2,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_mass: max,min SNOW=', & maxval(field2),minval(field2) write(lendian_out)field2 ! SNOW deallocate(field2,field2b,field2c) deallocate(datestr_all,varname_all,domainend_all,memoryorder_all) deallocate(start_block,end_block,start_byte,end_byte,file_offset) close(in_unit) close(lendian_out) end subroutine convert_binary_mass subroutine convert_binary_nmm(update_pint,ctph0,stph0,tlm0) !$$$ subprogram documentation block ! . . . . ! subprogram: convert_binary_nmm read wrf nmm binary restart ! prgmmr: parrish org: np22 date: 2003-09-05 ! ! abstract: using wrf library routines, read a wrf nmm binary ! format restart file. write the result to temporary binary ! file expected by read_wrf_nmm_guess. ! ! program history log: ! 2004-09-10 parrish ! 2004-11-05 wu - add check on input wrf guess file, stop code if problem ! 2004-11-11 parrish - change so byte offset information is written, instead ! of a whole new binary file--this done so mpi-io ! can be used later to read/write directly from ! the wrf mass core binary file. also, do inventory ! of whole file, so offsets are general--not dependent ! on type of file. ! 2004-12-15 treadon - remove get_lun, read guess from file "wrf_inout" ! 2005-07-06 parrish - add read of pint byte address ! 2005-10-17 parrish - add ctph0,stph0,tlm0 ! 2006-04-06 middlecoff - changed out_unit from 55 to lendian_out ! 2006-06-19 wu - changes to allow nfldsig=3 (multiple first guess) ! 2007-04-12 parrish - add modifications to allow any combination of ikj or ijk ! grid ordering for input 3D fields ! ! input argument list: ! update_pint: false on input ! ! output argument list: ! update_pint: true on output if field pint (non-hydrostatic pressure in nmm model) ! is available, in which case pint gets updated by analysis increment of pd, ! the nmm hydrostatic pressure thickness variable. ! ctph0,stph0: cos and sin thp0, earth lat of center of nmm grid ! (0 deg lat in rotated nmm coordinate) ! (used by calctends routines) ! tlm0 ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_single,i_llong,r_kind,i_kind use constants, only: zero,half,rad2deg use gsi_4dvar, only: nhr_assimilation use gsi_io, only: lendian_out implicit none integer(i_kind),parameter:: in_unit = 15 logical ,intent(inout) :: update_pint real(r_kind),intent( out) :: ctph0,stph0,tlm0 character(9) wrfges,fileout integer(i_kind),allocatable:: start_block(:),end_block(:) integer(i_kind),allocatable:: start_byte(:),end_byte(:) integer(i_llong),allocatable:: file_offset(:) integer(i_llong) n_position character(132),allocatable:: datestr_all(:),varname_all(:),memoryorder_all(:) integer(i_kind),allocatable:: domainend_all(:,:) integer(i_kind) nrecs integer(i_kind) status_hdr integer(i_kind) hdrbuf(512) integer(i_kind) iyear,imonth,iday,ihour,iminute,isecond integer(i_kind) nlon_regional,nlat_regional,nsig_regional real(r_single) dlmd_regional,dphd_regional,pt_regional,pdtop_regional real(r_single) dy_nmm integer(i_kind) k,n real(r_single),allocatable::field1(:),field1p(:),field2(:,:),field2b(:,:) integer(i_kind) ksize integer(i_kind) index n_loop: do n=1,3 if(n==1)then wrfges = 'wrf_inout' else write(wrfges,'("wrf_inou",i1.1)')n endif open(in_unit,file=trim(wrfges),form='unformatted') write(6,*)' convert_binary_nmm: in_unit,lendian_out=',in_unit,lendian_out ! Check for valid input file read(in_unit,iostat=status_hdr)hdrbuf if(n==1)then if(status_hdr /= 0) then write(6,*)'CONVERT_BINARY_NMM: problem with wrfges = ',& trim(wrfges),', Status = ',status_hdr call stop2(74) endif else if(status_hdr /= 0) then write(6,*)'CONVERT_BINARY_NMM: no off hour guess = ', fileout close(in_unit) cycle n_loop endif endif write(fileout,'("sigf",i2.2)')n+nhr_assimilation-1 write(6,*)' convert_binary_nmm: in_unit,out_unit=',wrfges,',',fileout open(lendian_out,file=fileout,form='unformatted') rewind lendian_out ! reopen for direct access reading of sequential file close(in_unit) call count_recs_wrf_binary_file(in_unit,wrfges,nrecs) write(6,*) ' after count_recs_wrf_binary_file, nrecs=',nrecs allocate(datestr_all(nrecs),varname_all(nrecs),domainend_all(3,nrecs)) allocate(memoryorder_all(nrecs)) allocate(start_block(nrecs),end_block(nrecs),start_byte(nrecs),end_byte(nrecs),file_offset(nrecs)) call inventory_wrf_binary_file(in_unit,wrfges,nrecs, & datestr_all,varname_all,memoryorder_all,domainend_all, & start_block,end_block,start_byte,end_byte,file_offset) ! start with date record for date forecast was started ! y,m,d,h,m,s call retrieve_index(index,'START_DATE',varname_all,nrecs) if(index<0) stop read(datestr_all(index),'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') & iyear,imonth,iday,ihour,iminute,isecond write(6,*)' convert_binary_nmm: START_DATE =',& iyear,imonth,iday,ihour,iminute,isecond ! nlon_regional, nlat_regional, nsig_regional call retrieve_index(index,'T',varname_all,nrecs) if(index<0) stop if(trim(memoryorder_all(index))=='XZY') then nlon_regional=domainend_all(1,index) nlat_regional=domainend_all(3,index) nsig_regional=domainend_all(2,index) end if if(trim(memoryorder_all(index))=='XYZ') then nlon_regional=domainend_all(1,index) nlat_regional=domainend_all(2,index) nsig_regional=domainend_all(3,index) end if read(datestr_all(index),'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') & iyear,imonth,iday,ihour,iminute,isecond write(6,*)' convert_binary_nmm: iy,m,d,h,m,s=',& iyear,imonth,iday,ihour,iminute,isecond write(6,*)' convert_binary_nmm: nlon,lat,sig_regional=',& nlon_regional,nlat_regional,nsig_regional ! dlmd_regional call retrieve_index(index,'DLMD',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,dlmd_regional,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_nmm: dlmd_regional=',dlmd_regional ! dphd_regional call retrieve_index(index,'DPHD',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,dphd_regional,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_nmm: dphd_regional=',dphd_regional ! pt_regional call retrieve_index(index,'PT',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,pt_regional,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_nmm: pt_regional=',pt_regional ! pdtop_regional call retrieve_index(index,'PDTOP',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,pdtop_regional,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_nmm: pdtop_regional=',pdtop_regional write(lendian_out) iyear,imonth,iday,ihour,iminute,isecond, & nlon_regional,nlat_regional,nsig_regional, & dlmd_regional,dphd_regional,pt_regional,pdtop_regional allocate(field1(nsig_regional),field1p(nsig_regional+1)) ! deta1 call retrieve_index(index,'DETA1',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,field1,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) do k=1,nsig_regional write(6,*)' convert_binary_nmm: k,deta1(k)=',k,field1(k) end do write(lendian_out)field1 ! DETA1 ! aeta1 call retrieve_index(index,'AETA1',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,field1,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) do k=1,nsig_regional write(6,*)' convert_binary_nmm: k,aeta1(k)=',k,field1(k) end do write(lendian_out)field1 ! AETA1 ! eta1 call retrieve_index(index,'ETA1',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,field1p,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) do k=1,nsig_regional+1 write(6,*)' convert_binary_nmm: k,eta1(k)=',k,field1p(k) end do write(lendian_out)field1p ! ETA1 ! deta2 call retrieve_index(index,'DETA2',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,field1,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) do k=1,nsig_regional write(6,*)' convert_binary_nmm: k,deta2(k)=',k,field1(k) end do write(lendian_out)field1 ! DETA2 ! aeta2 call retrieve_index(index,'AETA2',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,field1,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) do k=1,nsig_regional write(6,*)' convert_binary_nmm: k,aeta2(k)=',k,field1(k) end do write(lendian_out)field1 ! AETA2 ! eta2 call retrieve_index(index,'ETA2',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,field1p,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) do k=1,nsig_regional+1 write(6,*)' convert_binary_nmm: k,eta2(k)=',k,field1p(k) end do write(lendian_out)field1p ! ETA2 deallocate(field1,field1p) allocate(field2(nlon_regional,nlat_regional)) allocate(field2b(nlon_regional,nlat_regional)) ! GLAT call retrieve_index(index,'GLAT',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,field2,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_nmm: max,min GLAT=', & rad2deg*maxval(field2),rad2deg*minval(field2) write(6,*)' convert_binary_nmm: glat(1,1),glat(nlon,1)=', & rad2deg*field2(1,1),rad2deg*field2(nlon_regional,1) write(6,*)' convert_binary_nmm: glat(1,nlat),glat(nlon,nlat)=', & rad2deg*field2(1,nlat_regional),rad2deg*field2(nlon_regional,nlat_regional) write(6,*)' convert_binary_nmm: my guess at tph0d = ', & rad2deg*field2(1+(nlon_regional-1)/2,1+(nlat_regional-1)/2) ctph0=cos(field2(1+(nlon_regional-1)/2,1+(nlat_regional-1)/2)) stph0=sin(field2(1+(nlon_regional-1)/2,1+(nlat_regional-1)/2)) ! DX_NMM call retrieve_index(index,'DX_NMM',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,field2b,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_nmm: max,min DX_NMM=', & maxval(field2b),minval(field2b) write(6,*)' convert_binary_nmm: dx_nmm(1,1),dx_nmm(nlon,1)=', & field2b(1,1),field2b(nlon_regional,1) write(6,*)' convert_binary_nmm: dx_nmm(1,nlat),dx_nmm(nlon,nlat)=', & field2b(1,nlat_regional),field2b(nlon_regional,nlat_regional) write(lendian_out)field2,field2b ! GLAT,DX_NMM ! GLON call retrieve_index(index,'GLON',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,field2,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_nmm: max,min GLON=', & rad2deg*maxval(field2),rad2deg*minval(field2) write(6,*)' convert_binary_nmm: glon(1,1),glon(nlon,1)=', & rad2deg*field2(1,1),rad2deg*field2(nlon_regional,1) write(6,*)' convert_binary_nmm: glon(1,nlat),glon(nlon,nlat)=', & rad2deg*field2(1,nlat_regional),rad2deg*field2(nlon_regional,nlat_regional) write(6,*)' convert_binary_nmm: my guess at tlm0d = ', & half*rad2deg*(field2(1+(nlon_regional-1)/2,1+(nlat_regional-1)/2)+ & field2(2+(nlon_regional-1)/2,1+(nlat_regional-1)/2)) tlm0=half*(field2(1+(nlon_regional-1)/2,1+(nlat_regional-1)/2)+ & field2(2+(nlon_regional-1)/2,1+(nlat_regional-1)/2)) ! DY_NMM call retrieve_index(index,'DY_NMM',varname_all,nrecs) if(index<0) stop call retrieve_field(in_unit,wrfges,dy_nmm,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_nmm: DY_NMM=',dy_nmm field2b=dy_nmm write(lendian_out)field2,field2b ! GLON,DY_NMM write(lendian_out) wrfges ! index for START_DATE record call retrieve_index(index,'START_DATE',varname_all,nrecs) if(index<0) stop n_position=file_offset(index) write(lendian_out)n_position ! offset for START_DATE record ! PD call retrieve_index(index,'PD',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(lendian_out)n_position ! offset for PD ! FIS call retrieve_index(index,'FIS',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for FIS = ',n_position write(lendian_out)n_position ! offset for FIS ! PINT call retrieve_index(index,'PINT',varname_all,nrecs) update_pint=.false. if(index>0) then update_pint=.true. n_position=file_offset(index+1) write(6,*)' byte offset, memoryorder for PINT = ',n_position,memoryorder_all(index) write(lendian_out)n_position,memoryorder_all(index) ! offset for PINT ! end if ! T call retrieve_index(index,'T',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset, memoryorder for T = ',n_position,memoryorder_all(index) write(lendian_out)n_position,memoryorder_all(index) ! offset for T ! ! Q call retrieve_index(index,'Q',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset, memoryorder for Q = ',n_position,memoryorder_all(index) write(lendian_out)n_position,memoryorder_all(index) ! offset for Q ! ! U call retrieve_index(index,'U',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset, memoryorder for U = ',n_position,memoryorder_all(index) write(lendian_out)n_position,memoryorder_all(index) ! offset for U ! ! V call retrieve_index(index,'V',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset, memoryorder for V = ',n_position,memoryorder_all(index) write(lendian_out)n_position,memoryorder_all(index) ! offset for V ! ! SM call retrieve_index(index,'SM',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for SM = ',n_position write(lendian_out)n_position ! offset for SM ! ! SICE call retrieve_index(index,'SICE',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for SICE = ',n_position write(lendian_out)n_position ! offset for SICE ! ! SST call retrieve_index(index,'SST',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for SST = ',n_position write(lendian_out)n_position ! offset for SST ! ! IVGTYP call retrieve_index(index,'IVGTYP',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for IVGTYP = ',n_position write(lendian_out)n_position ! offset for IVGTYP ! ! ISLTYP call retrieve_index(index,'ISLTYP',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for ISLTYP = ',n_position write(lendian_out)n_position ! offset for ISLTYP ! ! VEGFRC call retrieve_index(index,'VEGFRC',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for VEGFRC = ',n_position write(lendian_out)n_position ! offset for VEGFRC ! ! SNO call retrieve_index(index,'SNO',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for SNO = ',n_position write(lendian_out)n_position ! offset for SNO ! ! U10 call retrieve_index(index,'U10',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for U10 = ',n_position write(lendian_out)n_position ! offset for U10 ! ! V10 call retrieve_index(index,'V10',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for V10 = ',n_position write(lendian_out)n_position ! offset for V10 ! ! SMC call retrieve_index(index,'SMC',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) if(trim(memoryorder_all(index))=='XZY') then ksize=domainend_all(2,index) end if if(trim(memoryorder_all(index))=='XYZ') then ksize=domainend_all(3,index) end if write(6,*)' byte offset, ksize, memoryorder for SMC = ',n_position,ksize,memoryorder_all(index) write(lendian_out)n_position,ksize,memoryorder_all(index) ! offset for SMC ! ! STC call retrieve_index(index,'STC',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) if(trim(memoryorder_all(index))=='XZY') then ksize=domainend_all(2,index) end if if(trim(memoryorder_all(index))=='XYZ') then ksize=domainend_all(3,index) end if write(6,*)' byte offset, ksize, memoryorder for STC = ',n_position,ksize,memoryorder_all(index) write(lendian_out)n_position,ksize,memoryorder_all(index) ! offset for STC ! ! TSK call retrieve_index(index,'TSK',varname_all,nrecs) if(index<0) stop n_position=file_offset(index+1) write(6,*)' byte offset for TSK = ',n_position write(lendian_out)n_position ! offset for TSK ! !????????????????????????????????????????????????????????????????read z0 here to see what it looks like call retrieve_index(index,'Z0',varname_all,nrecs) call retrieve_field(in_unit,wrfges,field2b,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_nmm: max,min Z0=', & maxval(field2b),minval(field2b) write(lendian_out)field2b ! Z0 !????????????????????????????????????????????????????????????????? call retrieve_index(index,'SST',varname_all,nrecs) call retrieve_field(in_unit,wrfges,field2b,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_nmm: max,min SST=', & maxval(field2b),minval(field2b) write(lendian_out)field2b ! SST call retrieve_index(index,'TSK',varname_all,nrecs) call retrieve_field(in_unit,wrfges,field2b,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_nmm: max,min TSK=', & maxval(field2b),minval(field2b) write(lendian_out)field2b ! TSK call retrieve_index(index,'SM',varname_all,nrecs) call retrieve_field(in_unit,wrfges,field2b,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_nmm: max,min SM=', & maxval(field2b),minval(field2b) write(lendian_out)field2b ! SM call retrieve_index(index,'SICE',varname_all,nrecs) call retrieve_field(in_unit,wrfges,field2b,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_nmm: max,min SICE=', & maxval(field2b),minval(field2b) write(lendian_out)field2b ! SICE call retrieve_index(index,'SNO',varname_all,nrecs) call retrieve_field(in_unit,wrfges,field2b,start_block(index+1),end_block(index+1), & start_byte(index+1),end_byte(index+1)) write(6,*)' convert_binary_nmm: max,min SNO=', & maxval(field2b),minval(field2b) write(lendian_out)field2b ! SNO deallocate(field2,field2b) deallocate(datestr_all,varname_all,domainend_all,memoryorder_all) deallocate(start_block,end_block,start_byte,end_byte,file_offset) close(in_unit) close(lendian_out) enddo n_loop end subroutine convert_binary_nmm subroutine convert_nems_nmmb(update_pint,ctph0,stph0,tlm0) !$$$ subprogram documentation block ! . . . . ! subprogram: convert_nems_nmmb read nems nmmb restart ! prgmmr: parrish org: np22 date: 2009-03-17 ! ! abstract: using nemsio library routines, read nems nmmb restart file ! write the result to temporary binary ! file expected by read_wrf_nmm_guess. ! ! program history log: ! 2004-09-10 parrish ! 2004-11-05 wu - add check on input wrf guess file, stop code if problem ! 2004-11-11 parrish - change so byte offset information is written, instead ! of a whole new binary file--this done so mpi-io ! can be used later to read/write directly from ! the wrf mass core binary file. also, do inventory ! of whole file, so offsets are general--not dependent ! on type of file. ! 2004-12-15 treadon - remove get_lun, read guess from file "wrf_inout" ! 2005-07-06 parrish - add read of pint byte address ! 2005-10-17 parrish - add ctph0,stph0,tlm0 ! 2006-04-06 middlecoff - changed out_unit from 55 to lendian_out ! 2006-06-19 wu - changes to allow nfldsig=3 (multiple first guess) ! 2007-04-12 parrish - add modifications to allow any combination of ikj or ijk ! grid ordering for input 3D fields ! ! input argument list: ! update_pint: false on input ! ! output argument list: ! update_pint: true on output if field pint (non-hydrostatic pressure in nmm model) ! is available, in which case pint gets updated by analysis increment of pd, ! the nmm hydrostatic pressure thickness variable. ! ctph0,stph0: cos and sin thp0, earth lat of center of nmm grid ! (0 deg lat in rotated nmm coordinate) ! (used by calctends routines) ! tlm0 ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_single,r_kind,i_kind use constants, only: one_tenth,half,deg2rad,rad2deg use gsi_4dvar, only: nhr_assimilation use gsi_io, only: lendian_out use nemsio_module, only: nemsio_init,nemsio_open,nemsio_close use nemsio_module, only: nemsio_gfile,nemsio_getfilehead,nemsio_getheadvar,nemsio_readrecv implicit none ! integer(i_kind),parameter:: in_unit = 15 real(r_kind),parameter:: r0_01 = 0.01_r_kind real(r_kind),parameter:: r100 = 100.0_r_kind real(r_kind),parameter:: rd_over_cp = 0.285725661955006982_r_kind logical ,intent(inout) :: update_pint real(r_kind),intent( out) :: ctph0,stph0,tlm0 type(nemsio_gfile) :: gfile character(255) wrfges,fileout integer(i_kind) iyear,imonth,iday,ihour,iminute,isecond integer(i_kind) nlon_regional,nlat_regional,nsig_regional real(r_single) dlmd_regional,dphd_regional,pt_regional,pdtop_regional real(r_single) dy_nmm integer(i_kind) i,j,ii,k,n real(r_single),allocatable::field1(:),field1p(:),field2(:),field2b(:),field2c(:) integer(i_kind) idate(7),nrec,iret character(4) gdatatype,modelname character(32) gtype integer(i_kind) nfhour,nfminute,nfsecondn,nfsecondd,nframe,ntrac,nsoil,nmeta logical extrameta real(r_single),allocatable,dimension(:):: dsg1,sgml1,sg1,dsg2,sgml2,sg2 real(r_single),allocatable,dimension(:,:):: glat,glon,dx,dy real(r_single),allocatable,dimension(:):: glata,glona,dxa,dya character(8),allocatable:: recname(:) character(16),allocatable :: reclevtyp(:) integer(i_kind),allocatable:: reclev(:) real(r_kind) date6,date7,second,fhour character(3) nmmb_verttype ! 'OLD' for old vertical coordinate definition ! old def: p = eta1*pdtop+eta2*(psfc-pdtop-ptop)+ptop ! 'NEW' for new vertical coordinate definition ! new def: p = eta1*pdtop+eta2*(psfc-ptop)+ptop real(r_kind) pd,pd_to_ps,psfc_this real(r_kind) factor,ratio call nemsio_init(iret=iret) if(iret/=0) then write(6,*)'CONVERT_NEMS_NMMB: problem with nemsio_init, Status = ',iret call stop2(74) end if n_loop: do n=1,3 if(n==1)then wrfges = 'wrf_inout' else write(wrfges,'("wrf_inou",i1.1)')n endif call nemsio_open(gfile,wrfges,'READ',iret=iret) write(6,*)' convert_nems_nmmb: nemsio_open, file name, iret=',trim(wrfges),iret if(n==1) then if(iret/=0) then write(6,*)'CONVERT_NEMS_NMMB: problem with wrfges = ',& trim(wrfges),', Status = ',iret call stop2(74) end if else if(iret/=0) then write(6,*)'CONVERT_NEMS_NMMB: no off hour guess = ',trim(wrfges) call nemsio_close(gfile,iret=iret) write(6,*)' close nemsio file, iret=',iret cycle n_loop end if end if write(fileout,'("sigf",i2.2)')n+nhr_assimilation-1 write(6,*)' convert_nems_nmmb: in_unit,out_unit=',trim(wrfges),',',trim(fileout) open(lendian_out,file=trim(fileout),form='unformatted') rewind lendian_out ! obtain model info (date, grid dims, etc.) ! idate(1) = year ! idate(2) = month ! idate(3) = day ! idate(4) = hour ! idate(5) = minute ! idate(6) = seconds*idate(7) ! idate(7) = scale factor for seconds call nemsio_getfilehead(gfile,iret=iret,nrec=nrec,dimx=nlon_regional,dimy=nlat_regional, & dimz=nsig_regional,idate=idate,gdatatype=gdatatype,gtype=gtype,modelname=modelname, & nfhour=nfhour,nfminute=nfminute,nfsecondn=nfsecondn,nfsecondd=nfsecondd, & nframe=nframe,ntrac=ntrac,nsoil=nsoil,extrameta=extrameta,nmeta=nmeta) write(6,*)' nemsio_getfilehead, iret=',iret,'nrec=',nrec write(6,*)' nlon_regional=',nlon_regional,' nlat_regional=',nlat_regional,' nsig_regional=',nsig_regional write(6,*)' idate=',idate write(6,*)' gdatatype=',gdatatype,' gtype=',trim(gtype) write(6,*)' nfhour=',nfhour,' nfminute=',nfminute,' nfsecondn=',nfsecondn,' nfsecondd=',nfsecondd write(6,*)' modelname=',modelname,' extrameta=',extrameta,' nframe=',nframe,' nmeta=',nmeta ! start with date record for date forecast was started ! y,m,d,h,m,s iyear= idate(1) imonth= idate(2) iday= idate(3) ihour= idate(4) iminute=idate(5) date6=idate(6) ; date7=idate(7) second=date6/date7 isecond=nint(second) write(6,*)' convert_nems_nmmb: START_DATE =',& iyear,imonth,iday,ihour,iminute,isecond ! dlmd_regional call nemsio_getheadvar(gfile,'DLMD',dlmd_regional,iret) write(6,*)' convert_nems_nmmb: dlmd_regional,iret=',dlmd_regional,iret ! dphd_regional call nemsio_getheadvar(gfile,'DPHD',dphd_regional,iret) write(6,*)' convert_nems_nmmb: dphd_regional,iret=',dphd_regional,iret ! pt_regional call nemsio_getheadvar(gfile,'PT',pt_regional,iret) write(6,*)' convert_nems_nmmb: pt_regional,iret=',pt_regional,iret ! pdtop_regional call nemsio_getheadvar(gfile,'PDTOP',pdtop_regional,iret) write(6,*)' convert_nems_nmmb: pdtop_regional,iret=',pdtop_regional,iret fhour=nfhour write(lendian_out) iyear,imonth,iday,ihour,iminute,isecond,fhour, & nlon_regional,nlat_regional,nsig_regional, & dlmd_regional,dphd_regional,pt_regional,pdtop_regional ! dsg1 (used to be deta1) allocate(dsg1(nsig_regional),field1(nsig_regional)) call nemsio_getheadvar(gfile,'DSG1',dsg1,iret) write(6,*)' convert_nems_nmmb: retrieve dsg1,iret=',iret do k=1,nsig_regional field1(k)=dsg1(nsig_regional+1-k) write(6,*)' convert_nems_nmmb: k,dsg1 (deta1)(k)=',k,field1(k) end do write(lendian_out)field1 ! DETA1 ! sgml1 (used to be aeta1) allocate(sgml1(nsig_regional)) call nemsio_getheadvar(gfile,'SGML1',sgml1,iret) write(6,*)' convert_nems_nmmb: retrieve sgml1,iret=',iret do k=1,nsig_regional field1(k)=sgml1(nsig_regional+1-k) write(6,*)' convert_nems_nmmb: k,sgml1 (aeta1)(k)=',k,field1(k) end do nmmb_verttype='OLD' if(field1(1)<.6_r_single) then nmmb_verttype='NEW' end if write(lendian_out)field1 ! AETA1 ! sg1 (used to be eta1) allocate(sg1(nsig_regional+1),field1p(nsig_regional+1)) call nemsio_getheadvar(gfile,'SG1',sg1,iret) write(6,*)' convert_nems_nmmb: retrieve sg1,iret=',iret do k=1,nsig_regional+1 field1p(k)=sg1(nsig_regional+2-k) write(6,*)' convert_nems_nmmb: k,sg1 (eta1)(k)=',k,field1p(k) end do write(lendian_out)field1p ! ETA1 ! dsg2 (used to be deta2) allocate(dsg2(nsig_regional)) call nemsio_getheadvar(gfile,'DSG2',dsg2,iret) write(6,*)' convert_nems_nmmb: retrieve dsg2,iret=',iret do k=1,nsig_regional field1(k)=dsg2(nsig_regional+1-k) write(6,*)' convert_nems_nmmb: k,dsg2 (deta2)(k)=',k,field1(k) end do write(lendian_out)field1 ! DETA2 ! sgml2 (used to be aeta2) allocate(sgml2(nsig_regional)) call nemsio_getheadvar(gfile,'SGML2',sgml2,iret) write(6,*)' convert_nems_nmmb: retrieve sgml2,iret=',iret do k=1,nsig_regional field1(k)=sgml2(nsig_regional+1-k) write(6,*)' convert_nems_nmmb: k,sgml2 (aeta2)(k)=',k,field1(k) end do write(lendian_out)field1 ! AETA2 ! sg2 (used to be eta2) allocate(sg2(nsig_regional+1)) call nemsio_getheadvar(gfile,'SG2',sg2,iret) write(6,*)' convert_nems_nmmb: retrieve sg2,iret=',iret do k=1,nsig_regional+1 field1p(k)=sg2(nsig_regional+2-k) write(6,*)' convert_nems_nmmb: k,sg2 (eta2)(k)=',k,field1p(k) end do write(lendian_out)field1p ! ETA2 write(6,*)' NEW NMMB EQUIVALENT VERTICAL COORDINATE PARAMETERS FOLLOW:' do k=1,nsig_regional write(6,'(" k,dsg1,sgml1,sg1,dsg2,sgml2,sg2=",i3,6f12.7)') & k,dsg1(k),sgml1(k),sg1(k),dsg2(k),sgml2(k),sg2(k) end do k=nsig_regional+1 write(6,'(" k, sg1, sg2=",i3,24x,f12.7,24x,f12.7)') & k, sg1(k), sg2(k) deallocate(field1,field1p,sg1,sg2,sgml1,sgml2,dsg1,dsg2) allocate(field2(nlon_regional*nlat_regional)) allocate(field2b(nlon_regional*nlat_regional)) allocate(field2c(nlon_regional*nlat_regional)) allocate(recname(nrec),reclevtyp(nrec),reclev(nrec)) allocate(glat(nlon_regional,nlat_regional),glon(nlon_regional,nlat_regional)) allocate( dx(nlon_regional,nlat_regional), dy(nlon_regional,nlat_regional)) allocate(glata(nlon_regional*nlat_regional),glona(nlon_regional*nlat_regional)) allocate( dxa(nlon_regional*nlat_regional), dya(nlon_regional*nlat_regional)) call nemsio_getfilehead(gfile,iret=iret,recname=recname,reclevtyp=reclevtyp, & reclev=reclev,lat=glata(:),lon=glona(:),dx=dxa(:),dy=dya(:)) ii=0 do j=1,nlat_regional do i=1,nlon_regional ii=ii+1 glat(i,j)=glata(ii)*deg2rad ! input lat in degrees glon(i,j)=glona(ii)*deg2rad ! input lon in degrees dx (i,j)= dxa (ii) dy (i,j)= dya (ii) end do end do ! GLAT write(6,*)' convert_nems_nmmb: max,min GLAT=', & rad2deg*maxval(glat),rad2deg*minval(glat) write(6,*)' convert_nems_nmmb: glat(1,1),glat(nlon,1)=', & rad2deg*glat(1,1),rad2deg*glat(nlon_regional,1) write(6,*)' convert_nems_nmmb: glat(1,nlat),glat(nlon,nlat)=', & rad2deg*glat(1,nlat_regional),rad2deg*glat(nlon_regional,nlat_regional) write(6,*)' convert_nems_nmmb: my guess at tph0d = ', & rad2deg*glat(1+(nlon_regional-1)/2,1+(nlat_regional-1)/2) ctph0=cos(glat(1+(nlon_regional-1)/2,1+(nlat_regional-1)/2)) stph0=sin(glat(1+(nlon_regional-1)/2,1+(nlat_regional-1)/2)) ! DX_NMM write(6,*)' convert_nems_nmmb: max,min DX_NMM=', & maxval(dx),minval(dx) write(6,*)' convert_nems_nmmb: dx_nmm(1,1),dx_nmm(nlon,1)=', & dx(1,1), dx(nlon_regional,1) write(6,*)' convert_nems_nmmb: dx_nmm(1,nlat),dx_nmm(nlon,nlat)=', & dx(1,nlat_regional), dx(nlon_regional,nlat_regional) write(lendian_out)glat,dx ! GLAT,DX_NMM !?????????????check to see if dx, dy backwards !?????????????? in existing wrf nmm ???????? ! GLON write(6,*)' convert_nems_nmmb: max,min GLON=',rad2deg*maxval( glon),rad2deg*minval( glon) write(6,*)' convert_nems_nmmb: glon(1,1),glon(nlon,1)=',rad2deg*glon(1,1),rad2deg*glon(nlon_regional,1) write(6,*)' convert_nems_nmmb: glon(1,nlat),glon(nlon,nlat)=', & rad2deg*glon(1,nlat_regional),rad2deg*glon(nlon_regional,nlat_regional) write(6,*)' convert_nems_nmmb: my guess at tlm0d = ', & half*rad2deg*( glon(1+(nlon_regional-1)/2,1+(nlat_regional-1)/2)+ & glon(2+(nlon_regional-1)/2,1+(nlat_regional-1)/2)) tlm0=half*( glon(1+(nlon_regional-1)/2,1+(nlat_regional-1)/2)+ & glon(2+(nlon_regional-1)/2,1+(nlat_regional-1)/2)) ! DY_NMM write(6,*)' convert_nems_nmmb: max,min DY_NMM=', & maxval(dy),minval(dy) write(6,*)' convert_nems_nmmb: dy_nmm(1,1),dy_nmm(nlon,1)=', & dy(1,1), dy(nlon_regional,1) write(6,*)' convert_nems_nmmb: dy_nmm(1,nlat),dy_nmm(nlon,nlat)=', & dy(1,nlat_regional), dy(nlon_regional,nlat_regional) write(lendian_out)glon,dy ! GLON,DY_NMM write(lendian_out) wrfges ! PINT call nemsio_readrecv(gfile,'pres','layer',1,field2(:),iret=iret) update_pint=.false. if(iret==0) update_pint=.true. write(6,*)' convert_nems_nmmb: pint, iret,update_pint=',iret,update_pint !????????????????????????????????????????????????????????????????read z0 here to see what it looks like call nemsio_readrecv(gfile,'zorl','sfc',1,field2b(:),iret=iret) write(6,*)' convert_nems_nmmb: iret,max,min Z0=',iret, & maxval(field2b),minval(field2b) write(lendian_out)field2b ! Z0 (?) ask if zorl is same as z0 !????????????????????????????????????????????????????????????????? call nemsio_readrecv(gfile,'tsea','sfc',1,field2b(:),iret=iret) write(6,*)' convert_nems_nmmb: iret,max,min SST=',iret, & maxval(field2b),minval(field2b) write(lendian_out)field2b ! SST !???????????????????????????????????????????? call nemsio_readrecv(gfile,'tg','sfc',1,field2b(:),iret=iret) write(6,*)' convert_nems_nmmb: iret,max,min TG=',iret, & maxval(field2b),minval(field2b) call nemsio_readrecv(gfile,'ths','sfc',1,field2b(:),iret=iret) write(6,*)' convert_nems_nmmb: iret,max,min THS=',iret, & maxval(field2b),minval(field2b) call nemsio_readrecv(gfile,'dpres','hybrid sig lev',1,field2c(:),iret=iret) write(6,*)' convert_nems_nmmb: iret,max,min PD=',iret, & maxval(field2c),minval(field2c) if(nmmb_verttype=='OLD') then pd_to_ps=r0_01*(pdtop_regional+pt_regional) else pd_to_ps=r0_01*pt_regional end if write(6,*)' pdtop_regional,pt_regional,pd_to_ps=',pdtop_regional,pt_regional,pd_to_ps do i=1,nlon_regional*nlat_regional pd=r0_01*field2c(i) psfc_this=pd+pd_to_ps ratio=(one_tenth*psfc_this/r100) factor=ratio**rd_over_cp field2c(i)=field2b(i)*factor end do write(6,*)' nmmb_verttype=',nmmb_verttype write(6,*)' max diff ths-ts=',maxval(field2b-field2c) write(6,*)' min diff ths-ts=',minval(field2b-field2c) write(6,*)' convert_nems_nmmb: iret,max,min TS=',iret, & maxval(field2c),minval(field2c) write(lendian_out)field2c ! TSK (ths converted to ts) !????????????????????????????????????????sm call nemsio_readrecv(gfile,'sm','sfc',1,field2b(:),iret=iret) write(6,*)' convert_nems_nmmb: iret,max,min SM=',iret, & maxval(field2b),minval(field2b) write(lendian_out)field2b ! SM call nemsio_readrecv(gfile,'sice','sfc',1,field2b(:),iret=iret) write(6,*)' convert_nems_nmmb: iret,max,min SICE=',iret, & maxval(field2b),minval(field2b) write(lendian_out)field2b ! SICE call nemsio_readrecv(gfile,'sno','sfc',1,field2b(:),iret=iret) write(6,*)' convert_nems_nmmb: iret,max,min SNO=',iret, & maxval(field2b),minval(field2b) write(lendian_out)field2b ! SNO deallocate(field2,field2b,field2c,recname,reclevtyp,reclev,glat,glon,dx,dy) call nemsio_close(gfile,iret=iret) write(6,*)' close nemsio file, iret=',iret close(lendian_out) enddo n_loop end subroutine convert_nems_nmmb #ifdef WRF subroutine count_recs_wrf_binary_file(in_unit,wrfges,nrecs) !$$$ subprogram documentation block ! . . . . ! subprogram: count_recs_binary_file count # recs on wrf binary file ! prgmmr: parrish org: np22 date: 2004-11-29 ! ! abstract: count number of sequential records contained in wrf binary ! file. this is done by opening the file in direct access ! mode with block length of 2**20, the size of the physical ! blocks on ibm "blue" and "white" machines. for optimal ! performance, change block length to correspond to the ! physical block length of host machine disk space. ! records are counted by looking for the 4 byte starting ! and ending sequential record markers, which contain the ! record size in bytes. only blocks are read which are known ! by simple calculation to contain these record markers. ! even though this is done on one processor, it is still ! very fast, and the time will always scale by the number of ! sequential records, not their size. this step and the ! following inventory step consistently take less than 0.1 seconds ! to complete. ! ! program history log: ! 2004-11-29 parrish ! 2005-02-17 todling, ifdef'ed wrf code out ! 2006-04-06 middlecoff - replace fortran open,close with openfileread,closefile ! ! input argument list: ! in_unit - fortran unit number where input file is opened through. ! wrfges - filename of input wrf binary restart file ! ! output argument list: ! nrecs - number of sequential records found on input wrf binary restart file. ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ ! do an initial read through of a wrf binary file, and get total number of sequential file records use kinds, only: i_byte,i_long,i_llong,i_kind implicit none integer(i_kind),intent(in ) :: in_unit character(9) ,intent(in ) :: wrfges integer(i_kind),intent( out) :: nrecs character(10) cwrfges integer(i_llong) nextbyte,locbyte,thisblock integer(i_byte) lenrec4(4) integer(i_long) lenrec,lensave equivalence (lenrec4(1),lenrec) integer(i_byte) missing4(4) integer(i_long) missing equivalence (missing,missing4(1)) integer(i_llong),parameter:: lrecl=2**20_i_llong integer(i_byte) buf(lrecl) integer(i_kind) i,loc_count,nreads logical lastbuf integer(i_kind) ierr cwrfges = wrfges cwrfges(10:10) = char(0) call openfileread (in_unit, ierr, cwrfges) ! open(in_unit,file=trim(wrfges),access='direct',recl=lrecl) nrecs=0 missing=-9999_i_long nextbyte=0_i_llong locbyte=lrecl nreads=0 lastbuf=.false. do ! get length of next record do i=1,4 nextbyte=nextbyte+1_i_llong locbyte=locbyte+1_i_llong if(locbyte > lrecl .and. lastbuf) go to 900 if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf) lenrec4(i)=buf(locbyte) end do if(lenrec <= 0_i_long .and. lastbuf) go to 900 if(lenrec <= 0_i_long .and. .not.lastbuf) go to 885 nextbyte=nextbyte+1_i_llong locbyte=locbyte+1_i_llong if(locbyte > lrecl .and. lastbuf) go to 900 if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf) nrecs=nrecs+1 loc_count=1 do i=2,4 if(loc_count>=lenrec) exit loc_count=loc_count+1 nextbyte=nextbyte+1_i_llong locbyte=locbyte+1_i_llong if(locbyte > lrecl .and. lastbuf) go to 900 if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf) end do do i=1,4 if(loc_count>=lenrec) exit loc_count=loc_count+1 nextbyte=nextbyte+1_i_llong locbyte=locbyte+1_i_llong if(locbyte > lrecl .and. lastbuf) go to 900 if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf) end do nextbyte=nextbyte-loc_count+lenrec locbyte=locbyte-loc_count+lenrec if(locbyte > lrecl .and. lastbuf) go to 900 if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf) lensave=lenrec do i=1,4 nextbyte=nextbyte+1_i_llong locbyte=locbyte+1_i_llong if(locbyte > lrecl .and. lastbuf) go to 900 if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf) lenrec4(i)=buf(locbyte) end do if(lenrec /= lensave) go to 890 end do 880 continue write(6,*)' reached impossible place in count_recs_wrf_binary_file' call closefile(in_unit,ierr) return 885 continue write(6,*)' problem in count_recs_wrf_binary_file, lenrec has bad value before end of file' write(6,*)' lenrec =',lenrec call closefile(in_unit,ierr) return 890 continue write(6,*)' problem in count_recs_wrf_binary_file, beginning and ending rec len words unequal' write(6,*)' begining reclen =',lensave write(6,*)' ending reclen =',lenrec write(6,*)' in_unit =',in_unit call closefile(in_unit,ierr) return 900 continue write(6,*)' normal end of file reached in count_recs_wrf_binary_file' write(6,*)' nblocks=',thisblock write(6,*)' nrecs=',nrecs write(6,*)' nreads=',nreads call closefile(in_unit,ierr) end subroutine count_recs_wrf_binary_file subroutine inventory_wrf_binary_file(in_unit,wrfges,nrecs, & datestr_all,varname_all,memoryorder_all,domainend_all, & start_block,end_block,start_byte,end_byte,file_offset) !$$$ subprogram documentation block ! . . . . ! subprogram: inventory_wrf_binary_file get contents of wrf binary file ! prgmmr: parrish org: np22 date: 2004-11-29 ! ! abstract: generate list of contents and map of wrf binary file which can be ! used for reading and writing with mpi io routines. ! same basic routine as count_recs_wrf_binary_file, except ! now wrf unpacking routines are used to decode wrf header ! records, and send back lists of variable mnemonics, dates, ! grid dimensions, and byte addresses relative to start of ! file for each field (this is used by mpi io routines). ! ! program history log: ! 2004-11-29 parrish ! 2006-04-06 middlecoff - replace fortran open,close with openfileread,closefile ! 2007-04-12 parrish - add output variable memoryorder_all to be used with modifications ! which allow any combination of ikj/ijk grid ordering for 3D fields. ! ! input argument list: ! in_unit - fortran unit number where input file is opened through. ! wrfges - filename of input wrf binary restart file ! nrecs - number of sequential records found on input wrf binary restart file. ! (obtained by a previous call to count_recs_wrf_binary_file) ! ! output argument list: (all following dimensioned nrecs) ! datestr_all - date character string for each field, where applicable (or else blanks) ! varname_all - wrf mnemonic for each variable, where applicable (or blank) ! memoryorder_all ! domainend_all - dimensions of each field, where applicable (up to 3 dimensions) ! start_block - direct access block number containing 1st byte of record ! (after 4 byte record mark) ! end_block - direct access block number containing last byte of record ! (before 4 byte record mark) ! start_byte - relative byte address in direct access block of 1st byte of record ! end_byte - relative byte address in direct access block of last byte of record ! file_offset - absolute address of byte before 1st byte of record (used by mpi io) ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: i_byte,i_long,i_llong,i_kind ! use module_internal_header_util, only: int_get_ti_header_char,int_get_write_field_header implicit none integer(i_kind) ,intent(in ) :: in_unit,nrecs character(9) ,intent(in ) :: wrfges character(132) ,intent( out) :: datestr_all(nrecs),varname_all(nrecs),memoryorder_all(nrecs) integer(i_kind) ,intent( out) :: domainend_all(3,nrecs) integer(i_kind) ,intent( out) :: start_block(nrecs),end_block(nrecs) integer(i_kind) ,intent( out) :: start_byte(nrecs),end_byte(nrecs) integer(i_llong),intent( out) :: file_offset(nrecs) character(10) cwrfges integer(i_kind) irecs integer(i_llong) nextbyte,locbyte,thisblock integer(i_byte) lenrec4(4) integer(i_long) lenrec,lensave equivalence (lenrec4(1),lenrec) integer(i_byte) missing4(4) integer(i_long) missing equivalence (missing,missing4(1)) integer(i_llong),parameter:: lrecl=2**20_i_llong integer(i_byte) buf(lrecl) integer(i_kind) i,loc_count,nreads logical lastbuf integer(i_byte) hdrbuf4(2048) integer(i_long) hdrbuf(512) equivalence(hdrbuf(1),hdrbuf4(1)) integer(i_kind),parameter:: int_field = 530 integer(i_kind),parameter:: int_dom_ti_char = 220 integer(i_kind) hdrbufsize integer(i_kind) inttypesize integer(i_kind) datahandle character(128) element,dumstr,strdata integer(i_kind) loccode character(132) blanks integer(i_kind) typesize integer(i_kind) fieldtype integer(i_kind) domaindesc character(132) memoryorder,stagger,dimnames(3) integer(i_kind) domainstart(3),domainend(3) integer(i_kind) patchstart(3),patchend(3) character(132) datestr,varname integer(i_kind) itypesize integer(i_kind) ierr call wrf_sizeof_integer(itypesize) inttypesize=itypesize blanks=trim(' ') cwrfges = wrfges cwrfges(10:10) = char(0) call openfileread (in_unit, ierr, cwrfges) ! open(in_unit,file=trim(wrfges),access='direct',recl=lrecl) irecs=0 missing=-9999_i_long nextbyte=0_i_llong locbyte=lrecl nreads=0 lastbuf=.false. do ! get length of next record do i=1,4 nextbyte=nextbyte+1_i_llong locbyte=locbyte+1_i_llong if(locbyte > lrecl .and. lastbuf) go to 900 if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf) lenrec4(i)=buf(locbyte) end do if(lenrec <= 0_i_long .and. lastbuf) go to 900 if(lenrec <= 0_i_long .and. .not. lastbuf) go to 885 nextbyte=nextbyte+1_i_llong locbyte=locbyte+1_i_llong if(locbyte > lrecl .and. lastbuf) go to 900 if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf) irecs=irecs+1 start_block(irecs)=thisblock start_byte(irecs)=locbyte file_offset(irecs)=nextbyte-1_i_llong hdrbuf4(1)=buf(locbyte) hdrbuf4(2:4)=missing4(2:4) hdrbuf4(5:8)=missing4(1:4) datestr_all(irecs)=blanks varname_all(irecs)=blanks memoryorder_all(irecs)=blanks domainend_all(1:3,irecs)=0 loc_count=1 do i=2,8 if(loc_count>=lenrec) exit loc_count=loc_count+1 nextbyte=nextbyte+1_i_llong locbyte=locbyte+1_i_llong if(locbyte > lrecl .and. lastbuf) go to 900 if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf) hdrbuf4(i)=buf(locbyte) end do if(lenrec==2048_i_long) write(6,*)' irecs,hdrbuf(2),int_dom_ti_char,int_field=', & irecs,hdrbuf(2),int_dom_ti_char,int_field if(lenrec==2048_i_long.and.(hdrbuf(2) == int_dom_ti_char .or. hdrbuf(2) == int_field)) then ! bring in next full record, so we can unpack datestr, varname, and domainend do i=9,lenrec loc_count=loc_count+1 nextbyte=nextbyte+1_i_llong locbyte=locbyte+1_i_llong if(locbyte > lrecl .and. lastbuf) go to 900 if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf) hdrbuf4(i)=buf(locbyte) end do if(hdrbuf(2) == int_dom_ti_char) then call int_get_ti_header_char(hdrbuf,hdrbufsize,inttypesize, & datahandle,element,dumstr,strdata,loccode) varname_all(irecs)=trim(element) datestr_all(irecs)=trim(strdata) write(6,*)' irecs,varname,datestr = ',irecs,trim(varname_all(irecs)),trim(datestr_all(irecs)) else call int_get_write_field_header(hdrbuf,hdrbufsize,typesize, & datahandle,datestr,varname,fieldtype, & domaindesc,memoryorder,stagger,dimnames, & domainstart,domainend,patchstart,patchend) varname_all(irecs)=trim(varname) datestr_all(irecs)=trim(datestr) memoryorder_all(irecs)=trim(memoryorder) domainend_all(1:3,irecs)=domainend(1:3) write(6,*)' irecs,datestr,domend,varname = ', & irecs,trim(datestr_all(irecs)),domainend_all(1:3,irecs),trim(varname_all(irecs)) end if end if nextbyte=nextbyte-loc_count+lenrec locbyte=locbyte-loc_count+lenrec if(locbyte > lrecl .and. lastbuf) go to 900 if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf) end_block(irecs)=thisblock end_byte(irecs)=locbyte lensave=lenrec do i=1,4 nextbyte=nextbyte+1_i_llong locbyte=locbyte+1_i_llong if(locbyte > lrecl .and. lastbuf) go to 900 if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf) lenrec4(i)=buf(locbyte) end do if(lenrec /= lensave) go to 890 end do 880 continue write(6,*)' reached impossible place in inventory_wrf_binary_file' call closefile(in_unit,ierr) return 885 continue write(6,*)' problem in inventory_wrf_binary_file, lenrec has bad value before end of file' write(6,*)' lenrec =',lenrec call closefile(in_unit,ierr) return 890 continue write(6,*)' problem in inventory_wrf_binary_file, beginning and ending rec len words unequal' write(6,*)' begining reclen =',lensave write(6,*)' ending reclen =',lenrec write(6,*)' irecs =',irecs write(6,*)' nrecs =',nrecs call closefile(in_unit,ierr) return 900 continue write(6,*)' normal end of file reached in inventory_wrf_binary_file' write(6,*)' nblocks=',thisblock write(6,*)' irecs,nrecs=',irecs,nrecs write(6,*)' nreads=',nreads call closefile(in_unit,ierr) end subroutine inventory_wrf_binary_file #else /* Start no WRF-library block */ subroutine count_recs_wrf_binary_file() !$$$ subprogram documentation block ! . . . . ! subprogram: count_recs_binary_file count # recs on wrf binary file ! prgmmr: parrish org: np22 date: 2004-11-29 ! ! absract: dummy call, does nothing ! ! program history log: ! 2005-02-17 todling, ifdef'ed wrf code out ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ implicit none write(6,*)'COUNT_RECS_WRF_BINARY_FILE: dummy routine, does nothing' return end subroutine count_recs_wrf_binary_file subroutine inventory_wrf_binary_file() !$$$ subprogram documentation block ! . . . . ! subprogram: inventory_wrf_binary_file get contents of wrf binary file ! prgmmr: parrish org: np22 date: 2004-11-29 ! ! abstract: dummy routine, does nonthing ! ! program history log: ! 2005-02-17 todling, ifdef'ed wrf code out ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ implicit none write(6,*)'INVENTORY_WRF_BINARY_FILE: dummy routine, does nothing' return end subroutine inventory_wrf_binary_file #endif /* End no WRF-library block */ subroutine next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf) !$$$ subprogram documentation block ! . . . . ! subprogram: next_buf bring in next direct access block ! prgmmr: parrish org: np22 date: 2004-11-29 ! ! abstract: bring in next direct access block when needed, as the file is scanned ! from beginning to end during counting and inventory of records. ! (subroutines count_recs_wrf_binary_file and inventory_wrf_binary_file) ! ! program history log: ! 2004-11-29 parrish ! 2006-04-06 middlecoff - replace direct access read with getbytes ! ! input argument list: ! in_unit - fortran unit number where input file is opened through. ! nextbyte - byte number from beginning of file that is desired ! locbyte - byte number from beginning of last block read for desired byte ! lrecl - direct access block length ! nreads - number of blocks read before now (for diagnostic information only) ! lastbuf - logical, if true, then no more blocks, so return ! ! output argument list: ! buf - output array containing contents of next block ! locbyte - byte number from beginning of new block read for desired byte ! thisblock - number of new block being read by this routine ! nreads - number of blocks read now (for diagnostic information only) ! lastbuf - logical, if true, then at end of file. ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: i_byte,i_llong,i_kind implicit none integer(i_llong),intent(in ) :: lrecl integer(i_kind) ,intent(in ) :: in_unit integer(i_llong),intent(in ) :: nextbyte integer(i_byte) ,intent( out) :: buf(lrecl) integer(i_llong),intent( out) :: thisblock integer(i_kind) ,intent(inout) :: nreads integer(i_llong),intent(inout) :: locbyte logical ,intent(inout) :: lastbuf integer(i_kind) ierr if(lastbuf) return ierr=0 nreads=nreads+1 ! compute thisblock: thisblock = 1_i_llong + (nextbyte-1_i_llong)/lrecl locbyte = 1_i_llong+mod(locbyte-1_i_llong,lrecl) ! The Fortran standard does not ! - specify what iostat should be for a DA read past the EOF ! - provide a way to detect end-of-file for a DA file ! - apply the concept end-of-file to direct-access files ! Consequently,the standard does not specify what the contents ! of the buffer will be for locations past the EOF. ! Hence the replacement of the DA read below with the call ! to getbytes ! read(in_unit,rec=thisblock,iostat=ierr)buf ! lastbuf = ierr /= 0 #ifdef WRF call getbytes(in_unit, buf, thisblock, lrecl, ierr) #else buf(:)=-1 #endif /* WRF */ lastbuf = ierr == 1 end subroutine next_buf subroutine retrieve_index(index,string,varname_all,nrecs) !$$$ subprogram documentation block ! . . . . ! subprogram: retrieve_index get record number of desired variable ! prgmmr: parrish org: np22 date: 2004-11-29 ! ! abstract: by examining previously generated inventory of wrf binary restart file, ! find record number that contains the header record for variable ! identified by input character variable "string". ! ! program history log: ! 2004-11-29 parrish ! ! input argument list: ! string - mnemonic for variable desired ! varname_all - list of all mnemonics obtained from inventory of file ! nrecs - total number of sequential records counted in wrf ! binary restart file ! ! output argument list: ! index - desired record number ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: i_kind implicit none integer(i_kind),intent(in ) :: nrecs integer(i_kind),intent( out) :: index character(*) ,intent(in ) :: string character(132) ,intent(in ) :: varname_all(nrecs) integer(i_kind) i do i=1,nrecs if(trim(string) == trim(varname_all(i))) then index=i return end if end do write(6,*)'RETRIEVE_INDEX: ***PROBLEM*** reading wrf nmm binary file, ',& 'rec id "',trim(string),'" not found' index=-1 return end subroutine retrieve_index subroutine retrieve_field(in_unit,wrfges,out,start_block,end_block,start_byte,end_byte) !$$$ subprogram documentation block ! . . . . ! subprogram: retrieve_field retrieve field from wrf binary file ! prgmmr: parrish org: np22 date: 2004-11-29 ! ! abstract: still using direct access, retrieve a field from the wrf binary restart file. ! ! program history log: ! 2004-11-29 parrish ! ! input argument list: ! in_unit - fortran unit number where input file is opened through. ! wrfges - filename of input wrf binary restart file ! start_block - direct access block number containing 1st byte of record ! (after 4 byte record mark) ! end_block - direct access block number containing last byte of record ! (before 4 byte record mark) ! start_byte - relative byte address in direct access block of 1st byte of record ! end_byte - relative byte address in direct access block of last byte of record ! ! output argument list: ! out - output buffer where desired field is deposited ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: i_byte,i_kind implicit none integer(i_kind),intent(in ) :: in_unit character(9) ,intent(in ) :: wrfges integer(i_kind),intent(in ) :: start_block,end_block,start_byte,end_byte integer(i_byte),intent( out) :: out(*) integer(i_kind),parameter:: lrecl=2**20 integer(i_byte) buf(lrecl) integer(i_kind) i,ii,k,ibegin,iend,ierr open(in_unit,file=trim(wrfges),access='direct',recl=lrecl) write(6,*)'RETRIEVE_FIELD: start_block,end_block,s_,e_byte=',& start_block,end_block,start_byte,end_byte ii=0 do k=start_block,end_block read(in_unit,rec=k,iostat=ierr)buf ibegin=1 ; iend=lrecl if(k == start_block) ibegin=start_byte if(k == end_block) iend=end_byte do i=ibegin,iend ii=ii+1 out(ii)=buf(i) end do end do close(in_unit) end subroutine retrieve_field SUBROUTINE int_get_ti_header_char( hdrbuf, hdrbufsize, itypesize, & DataHandle, Element, VarName, Data, code ) !$$$ subprogram documentation block ! . . . . ! subprogram: int_get_ti_header_char ! prgmmr: ! ! abstract: Same as int_gen_ti_header_char except that Data is read from the file. ! ! program history log: ! 2008-03-31 safford - add subroutine doc block ! 2009-01-03 todling - wrapped unavailable routine int_get_ti_header_c within ifdef ! 2009-09-28 guo - flagged uninitialized variable to signal possible conflict. ! ! input argument list: ! hdrbuf - ! itypesize - ! Element - ! Data - ! VarName - ! ! output argument list: ! hdrbuf - ! hdrbufsize - ! Element - ! Data - ! VarName - ! DataHandle - ! code - ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: i_kind IMPLICIT NONE ! INCLUDE 'intio_tags.h' INTEGER(i_kind), INTENT(INOUT) :: hdrbuf(*) INTEGER(i_kind), INTENT( OUT) :: hdrbufsize INTEGER(i_kind), INTENT(IN ) :: itypesize CHARACTER*(*) , INTENT(INOUT) :: Element, Data, VarName INTEGER(i_kind), INTENT( OUT) :: DataHandle, code !Local INTEGER(i_kind) i, n, DummyCount, typesize CHARACTER * 132 dummyData ! logical, external :: debug_foo ! #ifdef WRF CALL int_get_ti_header_c ( hdrbuf, hdrbufsize, n, itypesize, typesize, & DataHandle, dummyData, DummyCount, code ) #else DataHandle=-1 code=-1 #endif /* WRF */ i = n/itypesize+1 ; CALL int_unpack_string ( Element, hdrbuf( i ), n ) ; i = i + n CALL int_unpack_string ( Data , hdrbuf( i ), n ) ; i = i + n CALL int_unpack_string ( VarName , hdrbuf( i ), n ) ; i = i + n hdrbufsize = hdrbuf(1) RETURN END SUBROUTINE int_get_ti_header_char SUBROUTINE int_get_write_field_header ( hdrbuf, hdrbufsize, ftypesize, & DataHandle , DateStr , VarName , FieldType , & DomainDesc , MemoryOrder , Stagger , DimNames , & DomainStart , DomainEnd , & PatchStart , PatchEnd ) !$$$ subprogram documentation block ! . . . . ! subprogram: int_get_write_field_header ! prgmmr: ! ! abstract: See documentation block in int_gen_write_field_header() for ! a description of a "write field" header. ! ! program history log: ! 2008-03-31 safford - add subroutine doc block ! ! input argument list: ! hdrbuf - ! ftypesize - ! DateStr - ! VarName - ! MemoryOrder ! Stagger ! DimNames ! ! output argument list: ! hdrbuf - ! hdrbufsize - ! ftypesize - ! DataHandle - ! DateStr - ! VarName - ! FieldType ! DomainDesc ! MemoryOrder ! Stagger ! DimNames ! DomainStart,DomainEnd ! PatchStart,PatchEnd ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: i_kind IMPLICIT NONE ! INCLUDE 'intio_tags.h' INTEGER(i_kind) , INTENT(INOUT) :: hdrbuf(*) INTEGER(i_kind) , INTENT( OUT) :: hdrbufsize INTEGER(i_kind) , INTENT(INOUT) :: ftypesize INTEGER(i_kind) , INTENT( OUT) :: DataHandle CHARACTER*(*) , INTENT(INOUT) :: DateStr CHARACTER*(*) , INTENT(INOUT) :: VarName INTEGER(i_kind) , INTENT( OUT) :: FieldType INTEGER(i_kind) , INTENT( OUT) :: DomainDesc CHARACTER*(*) , INTENT(INOUT) :: MemoryOrder CHARACTER*(*) , INTENT(INOUT) :: Stagger CHARACTER*(*) ,dimension(*), INTENT(INOUT) :: DimNames INTEGER(i_kind) ,dimension(*), INTENT( OUT) :: DomainStart, DomainEnd INTEGER(i_kind) ,dimension(*), INTENT( OUT) :: PatchStart, PatchEnd !Local integer(i_kind),parameter:: int_field = 530 CHARACTER*132 mess INTEGER(i_kind) i, n hdrbufsize = hdrbuf(1) IF ( hdrbuf(2) /= int_field ) THEN write(mess,*)'int_get_write_field_header: hdrbuf(2) ne int_field ',hdrbuf(2),int_field CALL wrf_error_fatal3 ( "module_internal_header_util.b" , 220 , mess ) ENDIF ftypesize = hdrbuf(3) i = 4 DataHandle = hdrbuf(i) ; i = i+1 call int_unpack_string( DateStr, hdrbuf(i), n ) ; i = i+n call int_unpack_string( VarName, hdrbuf(i), n ) ; i = i+n FieldType = hdrbuf(i) ; i = i+1 call int_unpack_string( MemoryOrder, hdrbuf(i), n ) ; i = i+n call int_unpack_string( Stagger, hdrbuf(i), n ) ; i = i+n call int_unpack_string( DimNames(1), hdrbuf(i), n ) ; i = i+n call int_unpack_string( DimNames(2), hdrbuf(i), n ) ; i = i+n call int_unpack_string( DimNames(3), hdrbuf(i), n ) ; i = i+n DomainStart(1) = hdrbuf(i) ; i = i+1 DomainStart(2) = hdrbuf(i) ; i = i+1 DomainStart(3) = hdrbuf(i) ; i = i+1 DomainEnd(1) = hdrbuf(i) ; i = i+1 DomainEnd(2) = hdrbuf(i) ; i = i+1 DomainEnd(3) = hdrbuf(i) ; i = i+1 PatchStart(1) = hdrbuf(i) ; i = i+1 PatchStart(2) = hdrbuf(i) ; i = i+1 PatchStart(3) = hdrbuf(i) ; i = i+1 PatchEnd(1) = hdrbuf(i) ; i = i+1 PatchEnd(2) = hdrbuf(i) ; i = i+1 PatchEnd(3) = hdrbuf(i) ; i = i+1 DomainDesc = hdrbuf(i) ; i = i+1 RETURN END SUBROUTINE int_get_write_field_header SUBROUTINE int_unpack_string ( str, buf, n ) !$$$ subprogram documentation block ! . . . . ! subprogram: int_unpack_string ! prgmmr: ! ! abstract: This routine is used to extract a string from a sequence of integers. ! The first integer is the string length. ! ! program history log: ! 2008-03-31 safford - add subroutine doc block ! ! input argument list: ! str - ! buf - ! ! output argument list: ! str - ! n - ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: i_kind IMPLICIT NONE CHARACTER*(*) , INTENT( OUT) :: str INTEGER(i_kind) , INTENT( OUT) :: n ! on return, N is the number of ints copied from buf INTEGER(i_kind), DIMENSION(*), INTENT(IN ) :: buf !Local INTEGER(i_kind) i INTEGER(i_kind) strlen strlen = buf(1) str = "" DO i = 1, strlen str(i:i) = char(buf(i+1)) ENDDO n = strlen + 1 END SUBROUTINE int_unpack_string #ifndef MACOS SUBROUTINE wrf_sizeof_integer( retval ) !$$$ subprogram documentation block ! . . . . ! subprogram: wrf_sizeof_integer ! prgmmr: ! ! abstract: ! ! program history log: ! 2008-03-31 safford - add subroutine doc block ! ! input argument list: ! ! output argument list: ! retval - integer size ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: i_kind IMPLICIT NONE INTEGER(i_kind),intent( out) :: retval ! IWORDSIZE is defined by CPP retval = i_kind RETURN END SUBROUTINE wrf_sizeof_integer #endif !WRF:DRIVER_LAYER:UTIL ! MODULE module_wrf_error !$$$ module documentation block ! ! module: module_wrf_error ! ! abstract: ! ! program history log: ! 2008-03-31 safford - add module and subroutine doc blocks ! ! subroutines included: ! wrf_at_debug_level --- ! init_module_wrf_error --- ! ! variable definitions: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: i_kind implicit none ! set default to private private ! set subroutines to public public :: wrf_at_debug_level public :: init_module_wrf_error ! set passed variables to public public :: wrf_debug_level INTEGER(i_kind) :: wrf_debug_level = 0 CONTAINS LOGICAL FUNCTION wrf_at_debug_level ( level ) !$$$ subprogram documentation block ! ! subprogram: wrf_at_debug_level ! ! abstract: ! ! program history log: ! 2008-03-31 safford - add subroutine doc block ! ! input argument list: ! level - debug level ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: i_kind IMPLICIT NONE INTEGER(i_kind) , INTENT(IN ) :: level wrf_at_debug_level = ( level <= wrf_debug_level ) RETURN END FUNCTION wrf_at_debug_level SUBROUTINE init_module_wrf_error !$$$ subprogram documentation block ! ! subprogram: init_module_wrf_error ! ! abstract: ! ! program history log: ! 2008-03-31 safford - add subroutine doc block ! ! input argument list: ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block IMPLICIT NONE END SUBROUTINE init_module_wrf_error END MODULE module_wrf_error SUBROUTINE set_wrf_debug_level ( level ) !$$$ subprogram documentation block ! ! subprogram: set_wrf_debug_level ! ! abstract: ! ! program history log: ! 2008-03-31 safford - add subroutine doc block ! ! input argument list: ! level - debug level ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block USE module_wrf_error use kinds, only: i_kind IMPLICIT NONE INTEGER(i_kind) , INTENT(IN ) :: level wrf_debug_level = level RETURN END SUBROUTINE set_wrf_debug_level SUBROUTINE get_wrf_debug_level ( level ) !$$$ subprogram documentation block ! ! subprogram: get_wrf_debug_level ! ! abstract: ! ! program history log: ! 2008-03-31 safford - add subroutine doc block ! ! input argument list: ! ! output argument list: ! level - debug level ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block USE module_wrf_error use kinds, only: i_kind IMPLICIT NONE INTEGER(i_kind) , INTENT( OUT) :: level level = wrf_debug_level RETURN END SUBROUTINE get_wrf_debug_level SUBROUTINE wrf_debug( level , str ) !$$$ subprogram documentation block ! ! subprogram: wrf_debug ! ! abstract: ! ! program history log: ! 2008-03-31 safford - add subroutine doc block ! ! input argument list: ! level - debug level ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block USE module_wrf_error use kinds, only: i_kind IMPLICIT NONE CHARACTER*(*) str INTEGER(i_kind) , INTENT (IN ) :: level INTEGER(i_kind) :: debug_level CALL get_wrf_debug_level( debug_level ) IF ( level <= debug_level ) THEN ! old behavior CALL wrf_message( str ) ENDIF RETURN END SUBROUTINE wrf_debug SUBROUTINE wrf_message( str ) !$$$ subprogram documentation block ! ! subprogram: wrf_message ! ! abstract: ! ! program history log: ! 2008-03-31 safford - add subroutine doc block ! ! input argument list: ! str ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block USE module_wrf_error IMPLICIT NONE CHARACTER*(*), intent(in ) :: str write(6,*) TRIM(str) print*, TRIM(str) !TBH: Calls to logwrite cut off str in ESMF 2.2.0. !TBH: Restore this call once this ESMF bug is fixed. !TBH#ifdef USE_LOGERR !TBH IF ( WRFU_IsInitialized() ) THEN !TBH CALL WRFU_LogWrite( TRIM(str), WRFU_LOG_INFO ) !TBH ENDIF !TBH#endif END SUBROUTINE wrf_message ! intentionally write to stderr only SUBROUTINE wrf_message2( str ) !$$$ subprogram documentation block ! ! subprogram: wrf_message2 ! ! abstract: ! ! program history log: ! 2008-03-31 safford - add subroutine doc block ! ! input argument list: ! str ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block USE module_wrf_error IMPLICIT NONE CHARACTER*(*), intent(in ) :: str write(6,*) str !TBH: Calls to logwrite cut off str in ESMF 2.2.0. !TBH: Restore this call once this ESMF bug is fixed. !TBH#ifdef USE_LOGERR !TBH IF ( WRFU_IsInitialized() ) THEN !TBH CALL WRFU_LogWrite( str, WRFU_LOG_INFO ) !TBH ENDIF !TBH#endif END SUBROUTINE wrf_message2 SUBROUTINE wrf_error_fatal3( file_str, line, str ) !$$$ subprogram documentation block ! ! subprogram: wrf_error_fatal3 ! ! abstract: ! ! program history log: ! 2008-03-31 safford - add subroutine doc block ! ! input argument list: ! file_str - ! line - ! str - ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block USE module_wrf_error use kinds, only: i_kind IMPLICIT NONE CHARACTER*(*) , intent (in ) :: file_str INTEGER(i_kind) , INTENT (IN ) :: line ! only print file and line if line > 0 CHARACTER*(*) , intent (in ) :: str CHARACTER*256 :: line_str write(line_str,'(i6)') line CALL wrf_message( '-------------- FATAL CALLED ---------------' ) ! only print file and line if line is positive IF ( line > 0 ) THEN CALL wrf_message( 'FATAL CALLED FROM FILE: '//file_str//' LINE: '//TRIM(line_str) ) ENDIF CALL wrf_message( str ) CALL wrf_message( '-------------------------------------------' ) ! CALL wrf_abort call stop2(199) END SUBROUTINE wrf_error_fatal3 SUBROUTINE wrf_error_fatal( str ) !$$$ subprogram documentation block ! ! subprogram: wrf_error_fatal ! ! abstract: ! ! program history log: ! 2008-03-31 safford - add subroutine doc block ! ! input argument list: ! str - ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block USE module_wrf_error IMPLICIT NONE CHARACTER*(*),intent(in ) :: str CALL wrf_error_fatal3 ( ' ', 0, str ) END SUBROUTINE wrf_error_fatal SUBROUTINE wrf_check_error( expected, actual, str, file_str, line ) !$$$ subprogram documentation block ! ! subprogram: wrf_check_error ! ! abstract: Check to see if expected value == actual value ! If not, print message and exit. ! ! program history log: ! 2008-03-31 safford - add subroutine doc block ! ! input argument list: ! expected - ! actual - ! line - ! str ! file_str ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block USE module_wrf_error use kinds, only: i_kind IMPLICIT NONE INTEGER(i_kind) , INTENT (IN ) :: expected INTEGER(i_kind) , INTENT (IN ) :: actual CHARACTER*(*) , intent (in ) :: str CHARACTER*(*) , intent (in ) :: file_str INTEGER(i_kind) , INTENT (IN ) :: line CHARACTER (LEN=512) :: rc_str CHARACTER (LEN=512) :: str_with_rc IF ( expected /= actual ) THEN WRITE (rc_str,*) ' Routine returned error code = ',actual str_with_rc = TRIM(str // rc_str) CALL wrf_error_fatal3 ( file_str, line, str_with_rc ) ENDIF END SUBROUTINE wrf_check_error