MODULE module_wrf_to_grads_util USE module_wrf_to_grads_netcdf CONTAINS ! -- MMM version of 6/22/05 modified by RGF so datasets starting off the top of ! the hour will have time handled correctly in ctl file in tdef line ! !--------------------------------------------------------------------------------- subroutine create_grads( grads_file, file_for_time, file_time_index, times, & output_times, variables3d, number_of_3dvar, & variables2d, number_of_2dvar, & desc3d, desc2d, & output_input_grid, use_lowest_heights, grads_levels, & num_grads_levels, case_type, vertical_type, map, & debug1, debug2 ) implicit none include "netcdf.inc" integer :: output_times character (len=80), intent(in) :: grads_file character (len=80), intent(in), dimension(output_times) :: file_for_time, times integer, intent(in), dimension(output_times) :: file_time_index integer :: number_of_3dvar, number_of_2dvar character (len=20), dimension(number_of_3dvar) :: variables3d character (len=20), dimension(number_of_2dvar) :: variables2d character (len=50), dimension(number_of_3dvar) :: desc3d character (len=50), dimension(number_of_2dvar) :: desc2d logical, intent(in) :: output_input_grid, use_lowest_heights integer :: num_grads_levels real, dimension( num_grads_levels ), intent(in) :: grads_levels character (len=6) :: case_type character (len=1) :: vertical_type integer :: map logical, intent(in) :: debug1, debug2 real, allocatable, dimension(:,:,:) :: z, ph, phb real, allocatable, dimension(:,:,:) :: p, pb real, allocatable, dimension(:,:,:) :: data_out, data_out_z character (len=30) :: var_to_get, var_to_plot integer :: length_var, length_plot integer :: it integer :: ivar integer :: number_of_levels integer :: num_vert integer, dimension(2) :: loc_of_min_z real :: vert_args(100) logical :: soil real :: MISSING integer :: ncid, dimid, nf_status integer :: rcode, trcode real :: value_real integer :: nx, ny, nz integer :: nlgen integer :: ndims, dims(4) integer :: i, j, k integer :: ii, jj, kk integer :: ilon character (len=40) :: ctlfile, datfile character (len=38) :: tdef ! rgf integer, dimension(output_times) :: timestamp, datestamp character (len=19) :: wrf_time integer :: irec, file_recl real, allocatable, dimension(:,:) :: xlat, xlon real :: xlat_a(4), xlon_a(4) real :: xlat_n_max, xlat_s_max real :: xlon_w, xlon_e real :: abslatmin, abslatmax real :: abslonmin, abslonmax real :: truelat1, truelat2, temp real :: cen_lat, cen_lon real :: centeri, centerj integer :: ipoints, jpoints integer :: ncref, nrref real :: dx, dy real :: dxll integer :: map_proj parameter (MISSING=1.0E35) soil = .false. irec = 1 !================================================================================== ! need to pull out some data to set up dimensions, etc. ! nf_status = nf_open (file_for_time(1), nf_nowrite, ncid) call handle_err('Error opening file',nf_status) ! nf_status = nf_inq_dimid (ncid, 'west_east_stag', dimid) call handle_err('west_east_stag',nf_status) nf_status = nf_inq_dimlen (ncid, dimid, nx) call handle_err('Get NX',nf_status) nx = nx-1 ! nf_status = nf_inq_dimid (ncid, 'south_north_stag', dimid) call handle_err('south_north_stag',nf_status) nf_status = nf_inq_dimlen (ncid, dimid, ny) call handle_err('Get NY',nf_status) ny = ny-1 ! IF ( case_type /= 'static' ) THEN nf_status = nf_inq_dimid (ncid, 'bottom_top', dimid) call handle_err('bottom_top',nf_status) nf_status = nf_inq_dimlen (ncid, dimid, nz) call handle_err('Get NZ',nf_status) ELSE nz = 100 ENDIF nlgen = nz ! nf_status = nf_close (ncid) call handle_err('Error closing file',nf_status) ! !================================================================================== ! open output files ctlfile = trim(grads_file)//".ctl" datfile = trim(grads_file)//".dat" open (13, file=ctlfile) write (13, '("dset ^",a40)') datfile #ifdef bytesw write (13, '("options byteswapped")') #endif write (13, '("undef 1.e35")') #ifdef RECL4 file_recl = 4 #endif #ifdef RECL1 file_recl = 1 #endif if ( nx == 2 .and. ny /= 2 ) then open (15, file=datfile, form="unformatted",access="direct", & recl=(ny*file_recl)) elseif ( nx /= 2 .and. ny == 2 ) then open (15, file=datfile, form="unformatted",access="direct", & recl=(nx*file_recl)) elseif ( nx == 2 .and. ny == 2 ) then open (15, file=datfile, form="unformatted",access="direct", & recl=file_recl) else open (15, file=datfile, form="unformatted",access="direct", & recl=(nx*ny*file_recl)) endif !================================================================================== ! How will the vertical coordinate look like IF ( (.not. output_input_grid) .and. (.not. use_lowest_heights)) THEN ! we have user supplied vertical levels - CAN WE DO IT? nf_status = nf_open (file_for_time(1), nf_nowrite, ncid) call handle_err('Error opening file',nf_status) if ( vertical_type == 'p' ) then rcode = nf_inq_varid ( ncid, "P", dimid ) trcode = rcode rcode = nf_inq_varid ( ncid, "PB", dimid ) if ( nf_status == 0 ) rcode = trcode else if ( vertical_type == 'z' ) then rcode = nf_inq_varid ( ncid, "PH", dimid ) trcode = rcode rcode = nf_inq_varid ( ncid, "PHB", dimid ) if ( nf_status == 0 ) rcode = trcode endif nf_status = nf_close (ncid) call handle_err('Error closing file',nf_status) if ( rcode == 0 ) then ! we can do it write(6,*) ' ' write(6,*) " Asking to interpolate to ",vertical_type," levels - we can do that" write(6,*) ' ' number_of_levels = num_grads_levels vert_args(1:number_of_levels) = grads_levels(1:number_of_levels) else ! no interp, just put out computational grid write(6,*) ' ' write(6,*) ' FIELDS MISSING TO INTERPOLATE TO USER SPECIFIED VERTICAL LEVELS' write(6,*) ' WILL OUTPUT ON MODEL GRID' write(6,*) ' ' number_of_levels = nz vertical_type = 'n' endif END IF IF ( (output_input_grid) .and. (use_lowest_heights)) THEN ! use lowest column for z heights of grids - CAN WE DO IT? nf_status = nf_open (file_for_time(1), nf_nowrite, ncid) call handle_err('Error opening file',nf_status) rcode = nf_inq_varid ( ncid, "P", dimid ) trcode = rcode rcode = nf_inq_varid ( ncid, "PB", dimid ) if ( nf_status == 0 ) rcode = trcode nf_status = nf_close (ncid) call handle_err('Error closing file',nf_status) if ( rcode == 0 ) then ! we can do it write(6,*) ' ' write(6,*) " Asking to interpolate to lowerst h level - we can do that" write(6,*) ' ' allocate( z(nx,ny,nz) ) allocate( ph(nx,ny,nz+1) ) allocate( phb(nx,ny,nz+1) ) ! get base and perturbation height at full eta levels: call get_var_3d_real_cdf( file_for_time(1),'PH',ph, & nx,ny,nz+1,1,debug2 ) call get_var_3d_real_cdf( file_for_time(1),'PHB',phb, & nx,ny,nz+1,1,debug2 ) ! compute Z at half levels: ph = (ph+phb)/9.81 z = 0.5*(ph(:,:,1:nz)+ph(:,:,2:nz+1)) z = z/1000. ! convert to kilometers number_of_levels = nz vertical_type = 'z' loc_of_min_z = minloc(z(:,:,1)) vert_args(1:number_of_levels) = & z(loc_of_min_z(1),loc_of_min_z(2),1:number_of_levels) vert_args(1) = vert_args(1) + 0.002 vert_args(nz) = vert_args(nz) - 0.002 deallocate( z ) deallocate( ph ) deallocate( phb ) else ! no interp, just put out computational grid write(6,*) ' ' write(6,*) ' FIELDS MISSING TO INTERPOLATE TO HEIGHT LEVELS' write(6,*) ' WILL OUTPUT ON MODEL GRID' write(6,*) ' ' number_of_levels = nz vertical_type = 'n' endif END IF IF ( output_input_grid .and. (.not. use_lowest_heights)) THEN ! no interp, just put out computational grid write(6,*) " Will use model levels for output" number_of_levels = nz ENDIF !================================================================================== if(debug1) then write(6,*) ' number of levels = ',number_of_levels do k=1, number_of_levels write(6,*) ' k, vert_args(k) ',k,vert_args(k) enddo end if !================================================================================== ! work out times and time differences tdef = ' 11 linear 00:00z01jan2000 1hr' ! rgf IF ( case_type /= 'static' ) THEN do it = 1, output_times call time_calc(times(it), timestamp(it), datestamp(it), debug2 , & tdef, it, case_type ) enddo write (tdef(9:10),'(i2)') output_times ENDIF !================================================================================== ! try to get map information IF (map == 1) THEN call get_gl_att_int_cdf( file_for_time(1), 'MAP_PROJ', map_proj, debug2 ) if ( map_proj /= 0 ) then ! get more map parameters first call get_gl_att_real_cdf( file_for_time(1), 'DX', dx, debug2 ) call get_gl_att_real_cdf( file_for_time(1), 'DY', dy, debug2 ) call get_gl_att_real_cdf( file_for_time(1), 'CEN_LAT', cen_lat, debug2 ) call get_gl_att_real_cdf( file_for_time(1), 'TRUELAT1', truelat1, debug2 ) call get_gl_att_real_cdf( file_for_time(1), 'TRUELAT2', truelat2, debug2 ) nf_status = nf_open (file_for_time(1), nf_nowrite, ncid) call handle_err('Error opening file',nf_status) rcode = NF_GET_ATT_REAL(ncid, nf_global, 'STAND_LON', value_real ) nf_status = nf_close (ncid) call handle_err('Error closing file',nf_status) if ( rcode == 0) then call get_gl_att_real_cdf( file_for_time(1), 'STAND_LON', cen_lon, debug2 ) else write(6,*) ' ##### #####' write(6,*) ' ##### NOTE probably dealing with version 1 data #####' write(6,*) ' ##### Using CEN_LON in calculations #####' write(6,*) ' ##### Please check project of GrADS data #####' write(6,*) ' ##### #####' write(6,*) ' ' call get_gl_att_real_cdf( file_for_time(1), 'CEN_LON', cen_lon, debug2 ) endif allocate( xlat(nx,ny) ) allocate( xlon(nx,ny) ) call get_var_2d_real_cdf( file_for_time(1), 'XLAT', xlat, nx,ny, 1, debug2 ) call get_var_2d_real_cdf( file_for_time(1), 'XLONG',xlon, nx,ny, 1, debug2 ) end if if (map_proj == 0 .OR. map_proj == 3) then ! NO or MERCATOR ! check for dateline ilon = 0 if ( abs(xlon(1,1) - xlon(nx,ny)) .GT. 180. ) ilon = 1 IF ( ilon == 1 ) THEN write(13,'("xdef ",i4," linear ",f9.4," ",f8.4)') & nx,xlon(1,1),(abs(xlon(1,1)-(360.+xlon(nx,ny)))/(nx-1)) ELSE write(13,'("xdef ",i4," linear ",f9.4," ",f8.4)') & nx,xlon(1,1),(abs(xlon(1,1)-xlon(nx,ny))/(nx-1)) ENDIF write(13,'("ydef ",i4," linear ",f9.4," ",f8.4)') & ny,xlat(1,1),(abs(xlat(1,1)-xlat(nx,ny))/(ny-1)) if (vertical_type == 'n' ) then write (13, '("zdef ",i3, " linear 1 1")') number_of_levels else write(13,'("zdef ",i3, " levels ")') number_of_levels do k = 1,number_of_levels write(13,'(" ",f10.5)') vert_args(k) enddo endif else if (map_proj == 1) then ! LAMBERT-CONFORMAL ! make sure truelat1 is the larger number if (truelat1 < truelat2) then if (debug2) write (6,*) ' switching true lat values' temp = truelat1 truelat1 = truelat2 truelat2 = temp endif xlat_a(1) = xlat(1,1) xlat_a(2) = xlat(1,ny) xlat_a(3) = xlat(nx,1) xlat_a(4) = xlat(nx,ny) xlon_a(1) = xlon(1,1) xlon_a(2) = xlon(1,ny) xlon_a(3) = xlon(nx,1) xlon_a(4) = xlon(nx,ny) abslatmin = 99999. abslatmax = -99999. abslonmin = 99999. abslonmax = -99999. ! check for dateline ilon = 0 if ( abs(xlon_a(1) - xlon_a(3)) .GT. 180. ) ilon = 1 if ( abs(xlon_a(2) - xlon_a(4)) .GT. 180. ) ilon = 1 do i=1,4 abslatmin=min(abslatmin,xlat_a(i)) abslatmax=max(abslatmax,xlat_a(i)) IF ( xlon_a(i) .lt. 0.0 .AND. ilon .eq. 1 ) THEN abslonmin=min(abslonmin,360.+xlon_a(i)) abslonmax=max(abslonmax,360.+xlon_a(i)) ELSE abslonmin=min(abslonmin,xlon_a(i)) abslonmax=max(abslonmax,xlon_a(i)) ENDIF enddo ! ! xlat_s_max = -90. ! xlat_n_max = -90. ! ! do i = 1, nx ! xlat_s_max = max (xlat_s_max,xlat(i,1)) ! xlat_n_max = max (xlat_n_max,xlat(i,ny)) ! enddo ! xlon_w = xlon(1, ny) ! xlon_e = xlon(nx, ny) ! centeri = int((cen_lon-xlon_w)*((nx-1)/(xlon_e-xlon_w))+1) ! centerj = ((cen_lat-xlat_s_max)*((ny)/(xlat_n_max-xlat_s_max))) dxll = (dx/1000.)/111./2. ipoints = int((abslatmax-abslatmin+2)/dxll) jpoints = int((abslonmax-abslonmin+2)/dxll) write(13,'("pdef ",i3," ",i3," lcc ",f7.3," ",f8.3," ",& & f8.3," ",f8.3," ",f4.0," ",f4.0," ",f8.3," ",& & f7.0," ",f7.0)')& ! nx,ny,cen_lat,cen_lon,centeri,centerj,& nx,ny,xlat(1,1),xlon(1,1),1.0,1.0,& truelat1,truelat2,cen_lon,dx,dy write(13,'("xdef ",i4," linear ",f6.1," ",f12.8)') jpoints, & (abslonmin-1.),dxll write(13,'("ydef ",i4," linear ",f6.1," ",f12.8)') ipoints, & (abslatmin-1.),dxll if (vertical_type == 'n' ) then write (13, '("zdef ",i3, " linear 1 1")') number_of_levels else write(13,'("zdef ",i3, " levels ")') number_of_levels do k = 1,number_of_levels write(13,'(" ",f10.5)') vert_args(k) enddo endif elseif (map_proj == 2) then ! POLAR STEREO xlat_a(1) = xlat(1,1) xlat_a(2) = xlat(1,ny) xlat_a(3) = xlat(nx,1) xlat_a(4) = xlat(nx,ny) xlon_a(1) = xlon(1,1) xlon_a(2) = xlon(1,ny) xlon_a(3) = xlon(nx,1) xlon_a(4) = xlon(nx,ny) abslatmin = 99999. abslatmax = -99999. abslonmin = 99999. abslonmax = -99999. do i=1,4 abslatmin=min(abslatmin,xlat_a(i)) abslonmin=min(abslonmin,xlon_a(i)) abslatmax=max(abslatmax,xlat_a(i)) abslonmax=max(abslonmax,xlon_a(i)) enddo dxll = (dx/1000.)/111./2. ipoints = int((abslatmax-abslatmin+2)/dxll) + 20 jpoints = int((abslonmax-abslonmin+2)/dxll) + 20 ncref = nx/2 nrref = ny/2 write(13,'("pdef ",i3," ",i3," ops ",f8.3," ",f8.3," ",f12.4," ", & & f12.4," ",i4.1," ",i4.1," ",f12.2," ",f12.2)') & nx,ny,xlat(ncref,nrref), xlon(ncref,nrref),dx*0.1,dy*0.1, & ncref,nrref,dx,dy write(13,'("xdef ",i4," linear ",f6.1," ",f12.8)') jpoints, & (abslonmin-1.),dxll write(13,'("ydef ",i4," linear ",f6.1," ",f12.8)') ipoints, & (abslatmin-1.),dxll if (vertical_type == 'n' ) then write (13, '("zdef ",i3, " linear 1 1")') number_of_levels else write(13,'("zdef ",i3, " levels ")') number_of_levels do k = 1,number_of_levels write(13,'(" ",f10.5)') vert_args(k) enddo endif endif ! END of map projections ELSE if ( case_type == 'ideal' .and. nx == 2) then write (13, '("xdef 1 linear 0 0.0001")') else write (13, '("xdef ",i4, " linear 0 0.0001")') nx endif if ( case_type == 'ideal' .and. ny == 2) then write (13, '("ydef 1 linear 0 0.0001")') else write (13, '("ydef ",i4, " linear 0 0.0001")') ny endif if (vertical_type == 'n' ) then write (13, '("zdef ",i3, " linear 1 1")') number_of_levels else write(13,'("zdef ",i3, " levels ")') number_of_levels do k = 1,number_of_levels write(13,'(" ",f10.5)') vert_args(k) enddo endif ENDIF write (13, '("tdef",a38)') tdef ! rgf !================================================================================== ! create the grads file if(debug1) then write(6,*) ' number of 3d variables to process = ',number_of_3dvar write(6,*) ' number of 2d variables to process = ',number_of_2dvar do k=1,number_of_3dvar write(6,*) k, variables3d(k) enddo do k=1,number_of_2dvar write(6,*) k+number_of_3dvar, variables2d(k) enddo write(6,*) ' output times ' do k=1,output_times write(6,*) k, timestamp(k) enddo write(6,*) ' horizontal dims ',nx,ny write(6,*) ' vertical dims nl ' ! do k=1,number_of_3dvar ! write(6,*) k, nl(k) ! enddo ! do k=1+number_of_3dvar,number_of_2dvar+number_of_3dvar ! write(6,*) k, nl(k) ! enddo write(6,*) 'vert coord: ',vert_args(1:number_of_levels) endif !================================================================================== ! First check to see if we have all the variables write(6,*) ' CHECK to make sure we have all the variables in the input file' call check_all_variables ( file_for_time(1), & variables3d, desc3d, number_of_3dvar, & variables2d, desc2d, number_of_2dvar, & debug1 ) write (13, '("vars ",i3)') number_of_3dvar+number_of_2dvar !================================================================================== ! first, get some base-state-stuff if ( vertical_type == 'p' ) then allocate( p (nx, ny, nz) ) allocate( pb(nx, ny, nz) ) call get_var_3d_real_cdf( file_for_time(1),'PB',pb,nx,ny,nz,1,debug2 ) endif if ( vertical_type == 'z' ) then allocate( z (nx, ny, nz) ) allocate( ph (nx, ny, nz+1) ) allocate( phb(nx, ny, nz+1) ) call get_var_3d_real_cdf( file_for_time(1),'PHB',phb,nx,ny,nz+1,1,debug2 ) endif !=============================TIME LOOP STARTS HERE================================ print *,' ' print *,' ############### Begin time loop ############### ' print *,' ' do it = 1, output_times !================================================================================== ! first, get p/z if needed if ( vertical_type == 'p' ) then call get_var_3d_real_cdf( file_for_time(it),'P',p, nx, ny, nz, & file_time_index(it),debug2 ) p = p+pb endif if ( vertical_type == 'z' ) then call get_var_3d_real_cdf( file_for_time(it),'PH',ph, nx, ny, nz+1, & file_time_index(it),debug2 ) ph = (ph+phb)/9.81 z = 0.5*(ph(:,:,1:nz)+ph(:,:,2:nz+1)) ! need to convert to kilometers for coordinate z = z/1000. endif !=============================DEALING WITH 3D VARIABLES============================ do ivar = 1, number_of_3dvar var_to_get = variables3d(ivar) var_to_plot = var_to_get call check_special_variable( var_to_get ) length_var = len_trim(var_to_get) length_plot = len_trim(var_to_plot) if(debug2) write(6,*) ' getting dims for ',var_to_get(1:length_var) call get_dims_cdf( file_for_time(it), var_to_get(1:length_var), & dims, ndims, debug2 ) if ( dims(3) < nlgen ) then num_vert = dims(3) soil = .true. endif if(debug1) write(6,*) ' getting data for ',var_to_get(1:length_var) if (soil) then allocate ( data_out_z (nx, ny, num_vert) ) call get_var_3d_real_cdf( file_for_time(it),var_to_get(1:length_var), & data_out_z, nx, ny, num_vert, & file_time_index(it),debug2 ) else allocate ( data_out (nx, ny, nz) ) allocate ( data_out_z (nx, ny, number_of_levels) ) call g_output_3d (file_for_time(it), file_time_index(it), it, & var_to_plot, length_plot, nx,ny,nz, & data_out, debug2) if(debug2) write(6,*) 'number_of_levels = ', number_of_levels if ( vertical_type == 'p' ) then call interp_to_z( data_out , nx, ny, nz, & data_out_z, nx, ny, number_of_levels, & p/100., vert_args, missing, & vertical_type, debug1 ) else if ( vertical_type == 'z' ) then call interp_to_z( data_out , nx, ny, nz, & data_out_z, nx, ny, number_of_levels, & z, vert_args, missing, & vertical_type, debug1 ) else data_out_z = data_out endif num_vert = number_of_levels deallocate ( data_out ) endif if(debug1) write(6,*) ' writing out variable, time ',ivar,it wrf_time = times(it) write(6,*) ' time ',wrf_time,', output variable ',var_to_plot(1:length_plot) if (it == 1) write(13,'(a15,i3," 0 ",a50)') var_to_plot, num_vert, & desc3d(ivar) if ( case_type == 'ideal' .and. ny == 2 .and. nx == 2 )then do kk=1,num_vert write(15,rec=irec) ((data_out_z(ii,jj,kk),ii=2,2),jj=2,2) irec=irec+1 enddo else if ( case_type == 'ideal' .and. ny == 2 .and. nx .ne. 2 )then do kk=1,num_vert write(15,rec=irec) ((data_out_z(ii,jj,kk),ii=1,nx),jj=2,2) irec=irec+1 enddo else if ( case_type == 'ideal' .and. nx == 2 .and. ny .ne. 2 )then do kk=1,num_vert write(15,rec=irec) ((data_out_z(ii,jj,kk),ii=2,2),jj=1,ny) irec=irec+1 enddo else do kk=1,num_vert write(15,rec=irec) ((data_out_z(ii,jj,kk),ii=1,nx),jj=1,ny) irec=irec+1 enddo endif deallocate ( data_out_z ) soil = .false. enddo ! END of 3D LOOP !=============================DEALING WITH 2D VARIABLES============================ do ivar = 1, number_of_2dvar var_to_get = variables2d(ivar) var_to_plot = var_to_get length_var = len_trim(var_to_get) length_plot = len_trim(var_to_plot) allocate( data_out(nx, ny, 1) ) if(debug1) write(6,*) ' getting data for ',var_to_get(1:length_var) call g_output_2d (file_for_time(it), file_time_index(it), it, & var_to_plot, length_plot, nx,ny,nz, & data_out, debug2) if(debug1) write(6,*) ' writing out variable, time ',ivar+number_of_3dvar,it wrf_time = times(it) write(6,*) ' time ',wrf_time,', output variable ',var_to_plot(1:length_plot) if (it == 1) write(13,'(a15," 0 0 ",a50)') var_to_plot, desc2d(ivar) if ( case_type == 'ideal' .and. ny == 2 .and. nx == 2 )then write(15,rec=irec) ((data_out(ii,jj,1),ii=2,2),jj=2,2) irec=irec+1 elseif ( case_type == 'ideal' .and. ny == 2 .and. nx .ne. 2 )then write(15,rec=irec) ((data_out(ii,jj,1),ii=1,nx),jj=2,2) irec=irec+1 else if ( case_type == 'ideal' .and. nx == 2 .and. ny .ne. 2 )then write(15,rec=irec) ((data_out(ii,jj,1),ii=2,2),jj=1,ny) irec=irec+1 else write(15,rec=irec) ((data_out(ii,jj,1),ii=1,nx),jj=1,ny) irec=irec+1 endif deallocate(data_out) enddo ! END OF 2D VARIABLES print*," " enddo ! END of TIME LOOP !================================================================================== ! we're finished - clean up if ( vertical_type == 'p' ) then deallocate( p ) deallocate( pb ) endif if ( vertical_type == 'z' ) then deallocate( z ) deallocate( ph ) deallocate( phb ) endif write(13, '("endvars")' ) close (15) end subroutine create_grads !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ subroutine read_time_list( io_unit, times, max_times, output_times, & all_times, debug1, debug2 ) implicit none integer, intent(in) :: max_times, io_unit logical, intent(in ) :: debug1, debug2 character (len=80), intent(out) :: times(max_times) integer, intent(out) :: output_times logical, intent(out) :: all_times character (len=80) :: tmp_time integer :: itimes, length read(io_unit, *) output_times if(output_times > 0) then all_times = .false. itimes = 0 do itimes = itimes + 1 read(io_unit,fmt='(a80)') tmp_time length = max(1,index(tmp_time,' ')-1) if ( tmp_time(1:16) == 'end_of_time_list') exit if ( itimes .le. output_times ) then times(itimes) = tmp_time(1:length) if(debug1) write(6,*) ' output time ',itimes,' is ',trim(times(itimes)) endif enddo else all_times = .true. if(debug1) write(6,*) ' All times in file will be processed' do read(io_unit,fmt='(a80)') tmp_time if ( tmp_time(1:16) == 'end_of_time_list') exit enddo output_times = 0 end if end subroutine read_time_list !--------------------------------------------------------------------------------- subroutine read_variable_list( io_unit, & variables3d, desc3d, number_of_3dvar, & variables2d, desc2d, number_of_2dvar, & max_variables, debug1, debug2 ) implicit none integer, intent(in) :: max_variables, io_unit integer, intent(out) :: number_of_3dvar, number_of_2dvar character (len=20), intent(out) :: variables3d(max_variables) character (len=20), intent(out) :: variables2d(max_variables) character (len=50), intent(out) :: desc3d(max_variables) character (len=50), intent(out) :: desc2d(max_variables) character (len=20) :: tmp_var character (len=50) :: tmp_desc integer :: length logical, intent(in) :: debug1, debug2 logical read_complete_3d logical read_complete_2d read_complete_3d = .false. read_complete_2d = .false. number_of_3dvar = 0 number_of_2dvar = 0 do while( .not. read_complete_3d ) read(io_unit, fmt='(a20,2x,a50)') tmp_var, tmp_desc if(tmp_var(1:1) /= ' ' .and. tmp_var(1:17) /= 'end_of_3dvar_list') then number_of_3dvar = number_of_3dvar + 1 length = max(1,index(tmp_var,' ')-1) variables3d(number_of_3dvar) = tmp_var(1:length) desc3d(number_of_3dvar) = tmp_desc if(debug2) write(6,*) ' 3D variable ', number_of_3dvar,' ', & tmp_var(1:15), ': ', & trim(desc3d(number_of_3dvar)) else if(tmp_var(1:17) == 'end_of_3dvar_list') then read_complete_3d = .true. end if enddo do while( .not. read_complete_2d ) read(io_unit, fmt='(a20,2x,a50)') tmp_var, tmp_desc if(tmp_var(1:1) /= ' ' .and. tmp_var(1:17) /= 'end_of_2dvar_list') then number_of_2dvar = number_of_2dvar + 1 length = max(1,index(tmp_var,' ')-1) variables2d(number_of_2dvar) = tmp_var(1:length) desc2d(number_of_2dvar) = tmp_desc if(debug2) write(6,*) ' 2D variable ', number_of_2dvar,' ', & tmp_var(1:15), ': ', & trim(desc2d(number_of_2dvar)) else if(tmp_var(1:17) == 'end_of_2dvar_list') then read_complete_2d = .true. end if enddo end subroutine read_variable_list !--------------------------------------------------------------------------------- subroutine read_file_list( io_unit, files, max_files, & number_of_files, debug1, debug2 ) implicit none integer, intent(in) :: max_files, io_unit logical, intent(in) :: debug1, debug2 integer, intent(out) :: number_of_files character (len=80), intent(out) :: files(max_files) character (len=80) :: file_in integer :: length number_of_files = 0 do read(io_unit, fmt='(a80)') file_in if(file_in(1:16) == 'end_of_file_list') exit if(file_in(1:1) /= ' ') then number_of_files = number_of_files + 1 length = max(1,index(file_in,' ')-1) files(number_of_files) = file_in(1:length) if(debug1) write(6,*) ' file ', number_of_files,' ', & file_in(1:length) end if enddo end subroutine read_file_list !--------------------------------------------------------------------------------- subroutine create_time_file_list( data_files, & file_for_time, & file_time_index, & times, & max_times, & all_times, & output_times, & number_of_files, & debug ) implicit none integer, intent(in) :: max_times, number_of_files integer, intent(out) :: output_times character (len=80) :: data_files(max_times),times(max_times), & file_for_time(max_times) character (len=80) :: times_in_file(max_times) integer, intent(out), dimension(max_times) :: file_time_index logical, intent(in) :: debug, all_times character (len=80) :: input_time character (len=80) :: time1, time2, file_check character (len=19) :: wrf_time integer :: file_index, i, j, n_times, length,length1 logical :: already_on_list ! put together the list of "times" and "file_for_time". ! for each file, find the times in that file and either ! (1) see if it corresponds to a required time - then set the ! file for that time, or ! (2) if all times are desired, add if to the list if it does ! not duplicate a time already in the list if( .not. all_times ) then do i=1,output_times file_for_time(i) = 'no_time_found' enddo end if loop_over_files : do file_index = 1, number_of_files if(debug) write(6,*) ' getting times for file ',data_files(file_index) call get_times_cdf( data_files(file_index), times_in_file, & n_times, max_times, debug ) if( n_times <= 0 ) then write(6,*) ' no output times found in ',data_files(file_index) write(6,*) ' error stop ' stop end if if(debug) then write(6,*) ' times from netcdf file ' do i=1,n_times wrf_time = times_in_file(i) write(6,*) i,wrf_time enddo endif if(all_times) then if(debug) write(6,*) ' sorting all times ' do i=1, n_times ! length = max(1,index(times_in_file(i),' ')-1) length = 19 already_on_list = .false. time1 = times_in_file(i) do j=1, output_times time2 = times(j) if( time1(1:length) == time2(1:length)) already_on_list = .true. enddo if(.not.already_on_list) then output_times = output_times+1 times(output_times) = time1 file_for_time(output_times) = data_files(file_index) file_time_index(output_times) = i endif enddo else do i=1, n_times ! length = max(1,index(times_in_file(i),' ')-1) length = 19 do j=1, output_times time1 = times_in_file(i) time2 = times(j) if( time1(1:length) == time2(1:length)) then file_for_time(j) = data_files(file_index) file_time_index(j) = i end if enddo enddo end if enddo loop_over_files ! check here to see if we have data for all times if ! times were specified as input if( .not. all_times) then do i=1,output_times file_check = file_for_time(i) if( file_check(1:13) == 'no_time_found') then write(6,*) ' no data found for time ',times(i) write(6,*) ' error stop ' stop end if enddo end if if(debug) then write(6,*) ' time and file list ' do i=1,output_times time1 = times(i) file_check = file_for_time(i) length = max(1,index(file_check,' ')-1) length1 = 19 write(6,*) i,time1(1:length1),' ',file_check(1:length),file_time_index(i) enddo end if end subroutine create_time_file_list !------------------------------------------------------------------- subroutine time_calc( time, timestamp, datestamp, debug , tdef,it,case_type) implicit none character (len=80), intent(in) :: time character (len=38) :: tdef ! rgf character (len=5) :: case_type integer, intent(out) :: timestamp, datestamp logical, intent(in) :: debug integer :: hours, minutes, seconds, year, month, day,it,hour1,hourint integer :: mins1,minsint read(time(18:19),*) seconds read(time(15:16),*) minutes read(time(12:13),*) hours read(time(1:4),*) year read(time(6:7),*) month read(time(9:10),*) day if(debug) write(6,*) ' day, month, year, hours, minutes, seconds ' if(debug) write(6,*) day, month, year, hours, minutes, seconds if ( it == 1) then if (case_type(1:4) == 'real')then write (tdef(19:20),'(i2)') hours write (tdef(21:21),'(":")') ! rgf if ( minutes < 10 ) then ! RGF write (tdef(22:23),'("0",i1)') minutes ! RGF else ! RGF write (tdef(22:23),'(i2)') minutes ! RGF endif if ( day < 10 ) then write (tdef(26:26),'(i1)') day ! rgf else write (tdef(25:26),'(i2)') day ! rgf endif write (tdef(30:33),'(i4)') year ! rgf if (month == 1) write (tdef(27:29),'(a3)') 'jan' ! rgf if (month == 2) write (tdef(27:29),'(a3)') 'feb' ! rgf if (month == 3) write (tdef(27:29),'(a3)') 'mar' ! rgf if (month == 4) write (tdef(27:29),'(a3)') 'apr' ! rgf if (month == 5) write (tdef(27:29),'(a3)') 'may' ! rgf if (month == 6) write (tdef(27:29),'(a3)') 'jun' ! rgf if (month == 7) write (tdef(27:29),'(a3)') 'jul' ! rgf if (month == 8) write (tdef(27:29),'(a3)') 'aug' ! rgf if (month == 9) write (tdef(27:29),'(a3)') 'sep' ! rgf if (month ==10) write (tdef(27:29),'(a3)') 'oct' ! rgf if (month ==11) write (tdef(27:29),'(a3)') 'nov' ! rgf if (month ==12) write (tdef(27:29),'(a3)') 'dec' ! rgf endif hour1=hours mins1=minutes elseif ( it == 2) then hourint = abs(hours-hour1) minsint = abs(minutes-mins1) if (hourint == 0 ) then if (minsint == 0 ) minsint = 1 if(debug) write(6,*) "interval is",minsint write (tdef(37:38),'(a2)') "mn" ! rgf write (tdef(35:36),'(i2)') minsint ! rgf if(debug) write(6,*) "TDEF is",tdef else if(debug) write(6,*) "Interval is",hourint write (tdef(35:36),'(i2)') hourint ! rgf if(debug) write(6,*) "TDEF is",tdef endif endif timestamp = seconds+100*minutes+10000*hours if((year > 1800) .and. (year < 2000)) year = year-1900 if((year >= 2000)) year = year-2000 if(month >= 2) day = day+31 ! add january if(month >= 3) day = day+28 ! add february if(month >= 4) day = day+31 ! add march if(month >= 5) day = day+30 ! add april if(month >= 6) day = day+31 ! add may if(month >= 7) day = day+30 ! add june if(month >= 8) day = day+31 ! add july if(month >= 9) day = day+31 ! add august if(month >= 10) day = day+30 ! add september if(month >= 11) day = day+31 ! add october if(month >= 12) day = day+30 ! add november if((month > 2) .and. (mod(year,4) == 0)) day = day+1 ! get leap year day datestamp = (year)*1000 + day ! datestamp = (year+2100)*1000 + day if(debug) then write(6,*) ' time, timestamp, datestamp ',time(1:19),timestamp,datestamp endif end subroutine time_calc !------------------------------------------------------------------------- subroutine g_output_3d (file, file_time_index, it, var, length_var, & nx, ny, nz, data_out, debug) implicit none character (len=80) :: file integer :: it, file_time_index character (len=30) :: var integer :: length_var integer :: nx, ny, nz real, dimension(nx,ny,nz) :: data_out logical :: debug real, allocatable, dimension(:,:,:) :: data_tmp, data_tmp2 real, allocatable, dimension(:,:,:) :: u, v real, allocatable, dimension(:,:) :: xlat, xlon real, allocatable, dimension(:,:,:) :: z, ph, phb real, allocatable, dimension(:,:,:) :: p, pb real, allocatable, dimension(:,:,:) :: t, qv integer :: map_proj real :: cen_lon, truelat1, truelat2 REAL , PARAMETER :: g = 9.81 ! acceleration due to gravity (m {s}^-2) REAL , PARAMETER :: r_d = 287. REAL , PARAMETER :: r_v = 461.6 REAL , PARAMETER :: cp = 7.*r_d/2. REAL , PARAMETER :: cv = cp-r_d REAL , PARAMETER :: cliq = 4190. REAL , PARAMETER :: cice = 2106. REAL , PARAMETER :: psat = 610.78 REAL , PARAMETER :: rcv = r_d/cv REAL , PARAMETER :: rcp = r_d/cp REAL , PARAMETER :: c2 = cp * rcv REAL , PARAMETER :: T0 = 273.16 REAL , PARAMETER :: p1000mb = 100000. REAL , PARAMETER :: cpovcv = cp/(cp-r_d) REAL , PARAMETER :: cvovcp = 1./cpovcv REAL :: pp if(debug) then write(6,*) ' calculations for variable ',var end if if(var == 'U' ) then allocate ( data_tmp(nx+1,ny,nz) ) call get_var_3d_real_cdf( file,"U", data_tmp, nx+1, ny, nz, & file_time_index, debug ) data_out = 0.5*(data_tmp(1:nx,:,:)+data_tmp(2:nx+1,:,:)) deallocate ( data_tmp ) else if(var == 'V' ) then allocate ( data_tmp(nx,ny+1,nz) ) call get_var_3d_real_cdf( file,"V", data_tmp, nx, ny+1, nz, & file_time_index, debug ) data_out = 0.5*(data_tmp(:,1:ny,:)+data_tmp(:,2:ny+1,:)) deallocate ( data_tmp ) else if(var == 'UMET' ) then call get_gl_att_int_cdf ( file, 'MAP_PROJ', map_proj, debug ) IF ( map_proj == 1 .OR. map_proj == 2 ) THEN allocate ( u(nx,ny,nz) ) allocate ( v(nx,ny,nz) ) allocate ( xlat(nx,ny) ) allocate ( xlon(nx,ny) ) allocate ( data_tmp(nx+1,ny,nz) ) call get_var_3d_real_cdf( file,"U", data_tmp, nx+1, ny, nz, & file_time_index, debug ) u = 0.5*(data_tmp(1:nx,:,:)+data_tmp(2:nx+1,:,:)) deallocate ( data_tmp ) allocate ( data_tmp(nx,ny+1,nz) ) call get_var_3d_real_cdf( file,"V", data_tmp, nx, ny+1, nz, & file_time_index, debug ) v = 0.5*(data_tmp(:,1:ny,:)+data_tmp(:,2:ny+1,:)) deallocate ( data_tmp ) call get_gl_att_real_cdf( file, 'STAND_LON', cen_lon, debug ) call get_gl_att_real_cdf( file, 'TRUELAT1', truelat1, debug ) call get_gl_att_real_cdf( file, 'TRUELAT2', truelat2, debug ) call get_var_2d_real_cdf( file, 'XLAT', xlat,nx,ny, 1,debug ) call get_var_2d_real_cdf( file, 'XLONG',xlon,nx,ny, 1,debug ) call rotate_wind (u,v,nx,ny,nz,var, & map_proj,cen_lon,xlat,xlon, & truelat1,truelat2,data_out) deallocate ( xlat ) deallocate ( xlon ) deallocate ( u ) deallocate ( v ) ELSE allocate ( data_tmp(nx+1,ny,nz) ) call get_var_3d_real_cdf( file,"U", data_tmp, nx+1, ny, nz, & file_time_index, debug ) data_out = 0.5*(data_tmp(1:nx,:,:)+data_tmp(2:nx+1,:,:)) deallocate ( data_tmp ) ENDIF else if(var == 'VMET' ) then call get_gl_att_int_cdf ( file, 'MAP_PROJ', map_proj, debug ) IF ( map_proj == 1 .OR. map_proj == 2 ) THEN allocate ( u(nx,ny,nz) ) allocate ( v(nx,ny,nz) ) allocate ( xlat(nx,ny) ) allocate ( xlon(nx,ny) ) allocate ( data_tmp(nx+1,ny,nz) ) call get_var_3d_real_cdf( file,"U", data_tmp, nx+1, ny, nz, & file_time_index, debug ) u = 0.5*(data_tmp(1:nx,:,:)+data_tmp(2:nx+1,:,:)) deallocate ( data_tmp ) allocate ( data_tmp(nx,ny+1,nz) ) call get_var_3d_real_cdf( file,"V", data_tmp, nx, ny+1, nz, & file_time_index, debug ) v = 0.5*(data_tmp(:,1:ny,:)+data_tmp(:,2:ny+1,:)) deallocate ( data_tmp ) call get_gl_att_real_cdf( file, 'STAND_LON', cen_lon, debug ) call get_gl_att_real_cdf( file, 'TRUELAT1', truelat1, debug ) call get_gl_att_real_cdf( file, 'TRUELAT2', truelat2, debug ) call get_var_2d_real_cdf( file, 'XLAT', xlat,nx,ny, 1,debug ) call get_var_2d_real_cdf( file, 'XLONG',xlon,nx,ny, 1,debug ) call rotate_wind (u,v,nx,ny,nz,var, & map_proj,cen_lon,xlat,xlon, & truelat1,truelat2,data_out) deallocate ( xlat ) deallocate ( xlon ) deallocate ( u ) deallocate ( v ) ELSE allocate ( data_tmp(nx,ny+1,nz) ) call get_var_3d_real_cdf( file,"V", data_tmp, nx, ny+1, nz, & file_time_index, debug ) data_out = 0.5*(data_tmp(:,1:ny,:)+data_tmp(:,2:ny+1,:)) deallocate ( data_tmp ) ENDIF else if(var == 'W' ) then allocate ( data_tmp(nx,ny,nz+1) ) call get_var_3d_real_cdf( file,"W", data_tmp, nx, ny, nz+1, & file_time_index, debug ) data_out = 0.5*(data_tmp(:,:,1:nz)+data_tmp(:,:,2:nz+1)) deallocate ( data_tmp ) else if(var == 'P' ) then allocate ( p(nx,ny,nz) ) allocate ( pb(nx,ny,nz) ) call get_var_3d_real_cdf( file,"P", p, nx, ny, nz, & file_time_index, debug ) call get_var_3d_real_cdf( file,"PB", pb, nx, ny, nz, & file_time_index, debug ) data_out = (p+pb)*.01 deallocate ( p ) deallocate ( pb ) else if(var == 'Z' ) then allocate ( ph(nx,ny,nz+1) ) allocate ( phb(nx,ny,nz+1) ) call get_var_3d_real_cdf( file,"PH", ph, nx, ny, nz+1, & file_time_index, debug ) call get_var_3d_real_cdf( file,"PHB", phb, nx, ny, nz+1, & file_time_index, debug ) ph = (ph+phb)/9.81 data_out = 0.5*(ph(:,:,1:nz)+ph(:,:,2:nz+1)) deallocate ( ph ) deallocate ( phb ) else if(var == 'THETA' ) then call get_var_3d_real_cdf( file,"T", data_out, nx, ny, nz, & file_time_index, debug ) data_out = data_out + 300. else if(var == 'TK' ) then allocate ( p(nx,ny,nz) ) allocate ( pb(nx,ny,nz) ) allocate ( data_tmp(nx,ny,nz) ) call get_var_3d_real_cdf( file,"P", p, nx, ny, nz, & file_time_index, debug ) call get_var_3d_real_cdf( file,"PB", pb, nx, ny, nz, & file_time_index, debug ) p = p+pb call get_var_3d_real_cdf( file,"T", data_tmp, nx, ny, nz, & file_time_index, debug ) data_out = (data_tmp+300.)*(p/p1000mb)**rcp deallocate ( p ) deallocate ( pb ) deallocate ( data_tmp ) else if(var == 'TC' ) then allocate ( p(nx,ny,nz) ) allocate ( pb(nx,ny,nz) ) allocate ( data_tmp(nx,ny,nz) ) call get_var_3d_real_cdf( file,"P", p, nx, ny, nz, & file_time_index, debug ) call get_var_3d_real_cdf( file,"PB", pb, nx, ny, nz, & file_time_index, debug ) p = p+pb call get_var_3d_real_cdf( file,"T", data_tmp, nx, ny, nz, & file_time_index, debug ) data_out = (data_tmp+300.)*(p/p1000mb)**rcp -T0 deallocate ( p ) deallocate ( pb ) deallocate ( data_tmp ) else if(var == 'TD' ) then allocate ( p(nx,ny,nz) ) allocate ( pb(nx,ny,nz) ) allocate ( qv(nx,ny,nz) ) allocate ( data_tmp(nx,ny,nz) ) call get_var_3d_real_cdf( file,"P", p, nx, ny, nz, & file_time_index, debug ) call get_var_3d_real_cdf( file,"PB", pb, nx, ny, nz, & file_time_index, debug ) p = p+pb call get_var_3d_real_cdf( file,"QVAPOR", qv, nx, ny, nz, & file_time_index, debug ) data_tmp = qv*(p/100.)/(0.622+qv) data_tmp = AMAX1(data_tmp,0.001) data_out = (243.5*log(data_tmp)-440.8)/(19.48-log(data_tmp)) deallocate ( p ) deallocate ( pb ) deallocate ( qv ) deallocate ( data_tmp ) else if(var == 'RH' ) then allocate ( p(nx,ny,nz) ) allocate ( pb(nx,ny,nz) ) allocate ( qv(nx,ny,nz) ) allocate ( t(nx,ny,nz) ) allocate ( data_tmp(nx,ny,nz) ) allocate ( data_tmp2(nx,ny,nz) ) call get_var_3d_real_cdf( file,"P", p, nx, ny, nz, & file_time_index, debug ) call get_var_3d_real_cdf( file,"PB", pb, nx, ny, nz, & file_time_index, debug ) p = p+pb call get_var_3d_real_cdf( file,"T", t, nx, ny, nz, & file_time_index, debug ) call get_var_3d_real_cdf( file,"QVAPOR", qv, nx, ny, nz, & file_time_index, debug ) t = (t+300.)*(p/p1000mb)**rcp data_tmp2 = 10.*0.6112*exp(17.67*(t-T0)/(t-29.65)) data_tmp = 0.622*data_tmp2/(0.01 * p - (1.-0.622)*data_tmp2) data_out = 100.*AMAX1(AMIN1(qv/data_tmp,1.0),0.0) deallocate ( p ) deallocate ( pb ) deallocate ( qv ) deallocate ( t ) deallocate ( data_tmp ) deallocate ( data_tmp2 ) else call get_var_3d_real_cdf( file,var(1:length_var), & data_out, nx,ny,nz, & file_time_index, debug ) endif end subroutine g_output_3d !------------------------------------------------------------------------- subroutine g_output_2d (file, file_time_index, it, var, length_var, & nx, ny, nz, data_out, debug) implicit none character (len=80) :: file integer :: it, file_time_index character (len=30) :: var integer :: length_var integer :: nx, ny, nz real, dimension(nx,ny,1) :: data_out logical :: debug integer, allocatable, dimension(:,:,:) :: data_int real, allocatable, dimension(:,:,:) :: u10, v10 real, allocatable, dimension(:,:) :: xlat, xlon real, allocatable, dimension(:,:,:) :: z,ph,phb real, allocatable, dimension(:,:,:) :: p,pb real, allocatable, dimension(:,:,:) :: ts,qv integer :: map_proj real :: cen_lon, truelat1, truelat2 if(debug) then write(6,*) ' calculations for variable ',var end if if(var == 'slvl') then allocate ( z(nx,ny,nz) ) allocate ( ph(nx,ny,nz+1) ) allocate ( phb(nx,ny,nz+1) ) allocate ( p(nx,ny,nz) ) allocate ( pb(nx,ny,nz) ) allocate ( ts(nx,ny,nz) ) allocate ( qv(nx,ny,nz) ) call get_var_3d_real_cdf( file,"PH", ph, nx, ny,nz+1, & file_time_index, debug ) call get_var_3d_real_cdf( file,"PHB", phb, nx, ny,nz+1, & file_time_index, debug ) ph = (ph+phb)/9.81 z = 0.5*(ph(:,:,1:nz)+ph(:,:,2:nz+1)) call get_var_3d_real_cdf( file,"P", p, nx, ny,nz, & file_time_index, debug ) call get_var_3d_real_cdf( file,"PB", pb, nx, ny,nz, & file_time_index, debug ) p = p+pb call get_var_3d_real_cdf( file,"T", ts, nx, ny,nz, & file_time_index, debug ) call get_var_3d_real_cdf( file,"QVAPOR", qv, nx, ny,nz, & file_time_index, debug ) call compute_seaprs (nx, ny, nz, z, ts, p, qv, data_out, debug) deallocate ( z ) deallocate ( ph ) deallocate ( phb ) deallocate ( p ) deallocate ( pb ) deallocate ( ts ) deallocate ( qv ) else if(var == 'U10M' ) then call get_gl_att_int_cdf ( file, 'MAP_PROJ', map_proj, debug ) IF ( map_proj == 1 .OR. map_proj == 2 ) THEN allocate ( u10(nx,ny,1) ) allocate ( v10(nx,ny,1) ) allocate ( xlat(nx, ny) ) allocate ( xlon(nx, ny) ) call get_var_2d_real_cdf( file,"U10", u10, nx, ny, & file_time_index, debug ) call get_var_2d_real_cdf( file,"V10", v10, nx, ny, & file_time_index, debug ) call get_gl_att_real_cdf( file, 'STAND_LON', cen_lon, debug ) call get_gl_att_real_cdf( file, 'TRUELAT1', truelat1, debug ) call get_gl_att_real_cdf( file, 'TRUELAT2', truelat2, debug ) call get_var_2d_real_cdf( file, 'XLAT', xlat,nx,ny, 1,debug ) call get_var_2d_real_cdf( file, 'XLONG',xlon,nx,ny, 1,debug ) call rotate_wind (u10,v10,nx,ny,1,var, & map_proj,cen_lon,xlat,xlon, & truelat1,truelat2,data_out) deallocate ( xlat ) deallocate ( xlon ) deallocate ( u10 ) deallocate ( v10 ) ELSE call get_var_2d_real_cdf( file,"U10", data_out, nx, ny, & file_time_index, debug ) ENDIF else if(var == 'V10M' ) then call get_gl_att_int_cdf ( file, 'MAP_PROJ', map_proj, debug ) IF ( map_proj == 1 .OR. map_proj == 2 ) THEN allocate ( u10(nx,ny,1) ) allocate ( v10(nx,ny,1) ) allocate ( xlat(nx, ny) ) allocate ( xlon(nx, ny) ) call get_var_2d_real_cdf( file,"U10", u10, nx, ny, & file_time_index, debug ) call get_var_2d_real_cdf( file,"V10", v10, nx, ny, & file_time_index, debug ) call get_gl_att_real_cdf( file, 'STAND_LON', cen_lon, debug ) call get_gl_att_real_cdf( file, 'TRUELAT1', truelat1, debug ) call get_gl_att_real_cdf( file, 'TRUELAT2', truelat2, debug ) call get_var_2d_real_cdf( file, 'XLAT', xlat,nx,ny, 1,debug ) call get_var_2d_real_cdf( file, 'XLONG',xlon,nx,ny, 1,debug ) call rotate_wind (u10,v10,nx,ny,1,var, & map_proj,cen_lon,xlat,xlon, & truelat1,truelat2,data_out) deallocate ( xlat ) deallocate ( xlon ) deallocate ( u10 ) deallocate ( v10 ) ELSE call get_var_2d_real_cdf( file,"V10", data_out, nx, ny, & file_time_index, debug ) ENDIF else if(var == 'XLONG' ) then call get_var_2d_real_cdf( file,var(1:length_var), & data_out, nx,ny, & file_time_index, debug ) WHERE ( data_out < 0.0 ) data_out = data_out + 360.0 ENDWHERE else if(var == 'IVGTYP' .or. var == 'ISLTYP') then allocate (data_int(nx,ny,1)) call get_var_2d_int_cdf( file,var(1:length_var), & data_int, nx,ny, & file_time_index, debug ) data_out = data_int deallocate (data_int) else call get_var_2d_real_cdf( file,var(1:length_var), & data_out, nx,ny, & file_time_index, debug ) endif end subroutine g_output_2d !------------------------------------------------------------------ subroutine check_special_variable( var_to_get ) implicit none character (len=20), intent(inout) :: var_to_get if(var_to_get(1:6) == 'THETA ' .or. var_to_get(1:6) == 'TC ' & .or. var_to_get(1:6) == 'TK ') then var_to_get(1:6) = 'T ' else if(var_to_get(1:6) == 'TD ' .or. var_to_get(1:6) == 'RH ' ) then var_to_get(1:6) = 'QVAPOR' else if(var_to_get(1:2) == 'Z ') then var_to_get(1:6) = 'PH ' else if(var_to_get(1:4) == 'UMET') then var_to_get(1:6) = 'U ' else if(var_to_get(1:4) == 'VMET') then var_to_get(1:6) = 'V ' end if end subroutine check_special_variable !-------------------------------------------------------- subroutine interp_to_z( data_in , nx_in , ny_in , nz_in , & data_out, nx_out, ny_out, nz_out, & z_in, z_out, missing_value, & vertical_type, debug ) implicit none integer, intent(in) :: nx_in , ny_in , nz_in integer, intent(in) :: nx_out, ny_out, nz_out real, intent(in) :: missing_value real, dimension(nx_in , ny_in , nz_in ), intent(in ) :: data_in, z_in real, dimension(nx_out, ny_out, nz_out), intent(out) :: data_out real, dimension(nz_out), intent(in) :: z_out logical, intent(in) :: debug character (len=1) :: vertical_type real, dimension(nz_in) :: data_in_z, zz_in real, dimension(nz_out) :: data_out_z integer :: i,j,k do i=1,nx_in do j=1,ny_in do k=1,nz_in data_in_z(k) = data_in(i,j,k) zz_in(k) = z_in(i,j,k) enddo do k=1,nz_out data_out_z(k) = data_out(i,j,k) enddo call interp_1d( data_in_z, zz_in, nz_in, & data_out_z, z_out, nz_out, & vertical_type, missing_value ) do k=1,nz_out data_out(i,j,k) = data_out_z(k) enddo enddo enddo end subroutine interp_to_z !---------------------------------------------- subroutine interp_1d( a, xa, na, & b, xb, nb, vertical_type, missing_value ) implicit none integer, intent(in) :: na, nb real, intent(in), dimension(na) :: a, xa real, intent(in), dimension(nb) :: xb real, intent(out), dimension(nb) :: b real, intent(in) :: missing_value integer :: n_in, n_out logical :: interp real :: w1, w2 character (len=1) :: vertical_type if ( vertical_type == 'p' ) then do n_out = 1, nb b(n_out) = missing_value interp = .false. n_in = 1 do while ( (.not.interp) .and. (n_in < na) ) if( (xa(n_in) >= xb(n_out)) .and. & (xa(n_in+1) <= xb(n_out)) ) then interp = .true. w1 = (xa(n_in+1)-xb(n_out))/(xa(n_in+1)-xa(n_in)) w2 = 1. - w1 b(n_out) = w1*a(n_in) + w2*a(n_in+1) end if n_in = n_in +1 enddo enddo else do n_out = 1, nb b(n_out) = missing_value interp = .false. n_in = 1 do while ( (.not.interp) .and. (n_in < na) ) if( (xa(n_in) <= xb(n_out)) .and. & (xa(n_in+1) >= xb(n_out)) ) then interp = .true. w1 = (xa(n_in+1)-xb(n_out))/(xa(n_in+1)-xa(n_in)) w2 = 1. - w1 b(n_out) = w1*a(n_in) + w2*a(n_in+1) end if n_in = n_in +1 enddo enddo endif end subroutine interp_1d !------------------------------------------------------------------------- ! ! This routines has been taken "as is" from wrf_user_fortran_util_0.f ! ! This routine assumes ! index order is (i,j,k) ! wrf staggering ! units: pressure (Pa), temperature(K), height (m), mixing ratio (kg kg{-1}) ! availability of 3d p, t, and qv; 2d terrain; 1d half-level zeta string ! output units of SLP are Pa, but you should divide that by 100 for the ! weather weenies. ! virtual effects are included ! ! Dave subroutine compute_seaprs ( nx , ny , nz , & z, t , p , q , & sea_level_pressure,debug) ! & t_sea_level, t_surf, level ) IMPLICIT NONE ! Estimate sea level pressure. INTEGER nx , ny , nz REAL z(nx,ny,nz) REAL t(nx,ny,nz) , p(nx,ny,nz) , q(nx,ny,nz) ! The output is the 2d sea level pressure. REAL sea_level_pressure(nx,ny) INTEGER level(nx,ny) REAL t_surf(nx,ny) , t_sea_level(nx,ny) LOGICAL debug ! Some required physical constants: REAL R, G, GAMMA PARAMETER (R=287.04, G=9.81, GAMMA=0.0065) ! Specific constants for assumptions made in this routine: REAL TC, PCONST PARAMETER (TC=273.16+17.5, PCONST = 10000) LOGICAL ridiculous_mm5_test PARAMETER (ridiculous_mm5_test = .TRUE.) ! PARAMETER (ridiculous_mm5_test = .false.) ! Local variables: INTEGER i , j , k INTEGER klo , khi REAL plo , phi , tlo, thi , zlo , zhi REAL p_at_pconst , t_at_pconst , z_at_pconst REAL z_half_lowest REAL , PARAMETER :: cp = 7.*R/2. REAL , PARAMETER :: rcp = R/cp REAL , PARAMETER :: p1000mb = 100000. LOGICAL l1 , l2 , l3, found ! Find least zeta level that is PCONST Pa above the surface. We later use this ! level to extrapolate a surface pressure and temperature, which is supposed ! to reduce the effect of the diurnal heating cycle in the pressure field. t(:,:,:) = (t(:,:,:)+300.)*(p(:,:,:)/p1000mb)**rcp DO j = 1 , ny DO i = 1 , nx level(i,j) = -1 k = 1 found = .false. do while( (.not. found) .and. (k.le.nz)) IF ( p(i,j,k) .LT. p(i,j,1)-PCONST ) THEN level(i,j) = k found = .true. END IF k = k+1 END DO IF ( level(i,j) .EQ. -1 ) THEN PRINT '(A,I4,A)','Troubles finding level ', & NINT(PCONST)/100,' above ground.' PRINT '(A,I4,A,I4,A)', & 'Problems first occur at (',i,',',j,')' PRINT '(A,F6.1,A)', & 'Surface pressure = ',p(i,j,1)/100,' hPa.' STOP 'Error_in_finding_100_hPa_up' END IF END DO END DO ! Get temperature PCONST Pa above surface. Use this to extrapolate ! the temperature at the surface and down to sea level. DO j = 1 , ny DO i = 1 , nx klo = MAX ( level(i,j) - 1 , 1 ) khi = MIN ( klo + 1 , nz - 1 ) IF ( klo .EQ. khi ) THEN PRINT '(A)','Trapping levels are weird.' PRINT '(A,I3,A,I3,A)','klo = ',klo,', khi = ',khi, & ': and they should not be equal.' STOP 'Error_trapping_levels' END IF plo = p(i,j,klo) phi = p(i,j,khi) tlo = t(i,j,klo)*(1. + 0.608 * q(i,j,klo) ) thi = t(i,j,khi)*(1. + 0.608 * q(i,j,khi) ) ! zlo = zetahalf(klo)/ztop*(ztop-terrain(i,j))+terrain(i,j) ! zhi = zetahalf(khi)/ztop*(ztop-terrain(i,j))+terrain(i,j) zlo = z(i,j,klo) zhi = z(i,j,khi) p_at_pconst = p(i,j,1) - pconst t_at_pconst = thi-(thi-tlo)*LOG(p_at_pconst/phi)*LOG(plo/phi) z_at_pconst = zhi-(zhi-zlo)*LOG(p_at_pconst/phi)*LOG(plo/phi) t_surf(i,j) = t_at_pconst*(p(i,j,1)/p_at_pconst)**(gamma*R/g) t_sea_level(i,j) = t_at_pconst+gamma*z_at_pconst END DO END DO ! If we follow a traditional computation, there is a correction to the sea level ! temperature if both the surface and sea level temnperatures are *too* hot. IF ( ridiculous_mm5_test ) THEN DO j = 1 , ny DO i = 1 , nx l1 = t_sea_level(i,j) .LT. TC l2 = t_surf (i,j) .LE. TC l3 = .NOT. l1 IF ( l2 .AND. l3 ) THEN t_sea_level(i,j) = TC ELSE t_sea_level(i,j) = TC - 0.005*(t_surf(i,j)-TC)**2 END IF END DO END DO END IF ! The grand finale: ta da! DO j = 1 , ny DO i = 1 , nx ! z_half_lowest=zetahalf(1)/ztop*(ztop-terrain(i,j))+terrain(i,j) z_half_lowest=z(i,j,1) sea_level_pressure(i,j) = p(i,j,1) * & EXP((2.*g*z_half_lowest)/ & (R*(t_sea_level(i,j)+t_surf(i,j)))) sea_level_pressure(i,j) = sea_level_pressure(i,j)*0.01 END DO END DO if (debug) then print *,'sea pres input at weird location i=20,j=1,k=1' print *,'t=',t(20,1,1),t(20,2,1),t(20,3,1) print *,'z=',z(20,1,1),z(20,2,1),z(20,3,1) print *,'p=',p(20,1,1),p(20,2,1),p(20,3,1) print *,'slp=',sea_level_pressure(20,1), & sea_level_pressure(20,2),sea_level_pressure(20,3) endif ! print *,'t=',t(10:15,10:15,1),t(10:15,2,1),t(10:15,3,1) ! print *,'z=',z(10:15,1,1),z(10:15,2,1),z(10:15,3,1) ! print *,'p=',p(10:15,1,1),p(10:15,2,1),p(10:15,3,1) ! print *,'slp=',sea_level_pressure(10:15,10:15), & ! sea_level_pressure(10:15,10:15),sea_level_pressure(20,10:15) end subroutine compute_seaprs !------------------------------------------------------------------ subroutine rotate_wind (u,v,d1,d2,d3,var, & map_proj,cen_lon,xlat,xlon, & truelat1,truelat2,data_out) implicit none integer, intent(in) :: d1, d2, d3 real, dimension(d1,d2,d3) :: data_out integer :: map_proj,i,j,k real :: cen_lon, truelat1, truelat2, cone real, dimension(d1,d2,d3) :: u,v real, dimension(d1,d2) :: xlat, xlon, diff, alpha character (len=10), intent(in) :: var REAL , PARAMETER :: pii = 3.14159265 REAL , PARAMETER :: radians_per_degree = pii/180. cone = 1. ! PS if( map_proj .eq. 1) then ! Lambert Conformal mapping IF (ABS(truelat1-truelat2) .GT. 0.1) THEN cone=(ALOG(COS(truelat1*radians_per_degree))- & ALOG(COS(truelat2*radians_per_degree))) / & (ALOG(TAN((90.-ABS(truelat1))*radians_per_degree*0.5 ))- & ALOG(TAN((90.-ABS(truelat2))*radians_per_degree*0.5 )) ) ELSE cone = SIN(ABS(truelat1)*radians_per_degree ) ENDIF end if diff = xlon - cen_lon do i = 1, d1 do j = 1, d2 if(diff(i,j) .gt. 180.) then diff(i,j) = diff(i,j) - 360. end if if(diff(i,j) .lt. -180.) then diff(i,j) = diff(i,j) + 360. end if end do end do do i = 1, d1 do j = 1, d2 if(xlat(i,j) .lt. 0.) then alpha(i,j) = - diff(i,j) * cone * radians_per_degree else alpha(i,j) = diff(i,j) * cone * radians_per_degree end if end do end do if(var(1:1) .eq. "U") then do k=1,d3 data_out(:,:,k) = v(:,:,k)*sin(alpha) + u(:,:,k)*cos(alpha) end do else if(var(1:1) .eq. "V") then do k=1,d3 data_out(:,:,k) = v(:,:,k)*cos(alpha) - u(:,:,k)*sin(alpha) end do end if end subroutine rotate_wind !------------------------------------------------------------------ subroutine handle_err(rmarker,nf_status) include "netcdf.inc" integer nf_status character*(*) :: rmarker if (nf_status .ne. nf_noerr) then write(*,*) 'NetCDF error : ',rmarker write(*,*) ' ',nf_strerror(nf_status) stop endif end subroutine handle_err !------------------------------------------------------------------ subroutine check_all_variables ( file, & variables3d, desc3d, number_of_3dvar, & variables2d, desc2d, number_of_2dvar, & debug ) include "netcdf.inc" character (len=80) :: file integer :: number_of_3dvar ,number_of_2dvar character (len=20), dimension(number_of_3dvar) :: variables3d character (len=20), dimension(number_of_2dvar) :: variables2d character (len=50), dimension(number_of_3dvar) :: desc3d character (len=50), dimension(number_of_2dvar) :: desc2d logical :: debug integer :: nf_status, ncid, rcode, id_var, trcode integer :: missing3d, missing2d integer :: newi nf_status = nf_open (file, nf_nowrite, ncid) call handle_err('Error opening file',nf_status) missing3d = 0 do i = 1,number_of_3dvar if ( variables3d(i) == 'UMET' ) then rcode = nf_inq_varid ( ncid, "U", id_var ) trcode = rcode rcode = nf_inq_varid ( ncid, "V", id_var ) if ( rcode == 0 ) rcode = trcode else if ( variables3d(i) == 'VMET' ) then rcode = nf_inq_varid ( ncid, "U", id_var ) trcode = rcode rcode = nf_inq_varid ( ncid, "V", id_var ) if ( rcode == 0 ) rcode = trcode else if ( variables3d(i) == 'Z' ) then rcode = nf_inq_varid ( ncid, "PH", id_var ) trcode = rcode rcode = nf_inq_varid ( ncid, "PHB", id_var ) if ( rcode == 0 ) rcode = trcode else if ( variables3d(i) == 'P' ) then rcode = nf_inq_varid ( ncid, "P", id_var ) trcode = rcode rcode = nf_inq_varid ( ncid, "PB", id_var ) if ( rcode == 0 ) rcode = trcode else if ( variables3d(i) == 'THETA' ) then rcode = nf_inq_varid ( ncid, "T", id_var ) else if ( variables3d(i) == 'TK' ) then rcode = nf_inq_varid ( ncid, "P", id_var ) trcode = rcode rcode = nf_inq_varid ( ncid, "PB", id_var ) if ( rcode == 0 ) rcode = trcode trcode = rcode rcode = nf_inq_varid ( ncid, "T", id_var ) if ( rcode == 0 ) rcode = trcode else if ( variables3d(i) == 'TC' ) then rcode = nf_inq_varid ( ncid, "P", id_var ) trcode = rcode rcode = nf_inq_varid ( ncid, "PB", id_var ) if ( rcode == 0 ) rcode = trcode trcode = rcode rcode = nf_inq_varid ( ncid, "T", id_var ) if ( rcode == 0 ) rcode = trcode else if ( variables3d(i) == 'TD' ) then rcode = nf_inq_varid ( ncid, "P", id_var ) trcode = rcode rcode = nf_inq_varid ( ncid, "PB", id_var ) if ( rcode == 0 ) rcode = trcode trcode = rcode rcode = nf_inq_varid ( ncid, "QVAPOR", id_var ) if ( rcode == 0 ) rcode = trcode else if ( variables3d(i) == 'RH' ) then rcode = nf_inq_varid ( ncid, "P", id_var ) trcode = rcode rcode = nf_inq_varid ( ncid, "PB", id_var ) if ( rcode == 0 ) rcode = trcode trcode = rcode rcode = nf_inq_varid ( ncid, "T", id_var ) if ( rcode == 0 ) rcode = trcode trcode = rcode rcode = nf_inq_varid ( ncid, "QVAPOR", id_var ) if ( rcode == 0 ) rcode = trcode else rcode = nf_inq_varid ( ncid, variables3d(i), id_var ) endif if (rcode .ne. 0) then write(6,*) ' Not in file, remove from list: ',trim(variables3d(i)) variables3d(i) = ' ' desc3d(i) = ' ' missing3d = missing3d+1 endif enddo missing2d = 0 do i = 1,number_of_2dvar if ( variables2d(i) == 'U10M' ) then rcode = nf_inq_varid ( ncid, "U10", id_var ) trcode = rcode rcode = nf_inq_varid ( ncid, "V10", id_var ) if ( rcode == 0 ) rcode = trcode else if ( variables2d(i) == 'V10M' ) then rcode = nf_inq_varid ( ncid, "U10", id_var ) trcode = rcode rcode = nf_inq_varid ( ncid, "V10", id_var ) if ( rcode == 0 ) rcode = trcode else if ( variables2d(i) == 'slvl' ) then rcode = nf_inq_varid ( ncid, "P", id_var ) trcode = rcode rcode = nf_inq_varid ( ncid, "PB", id_var ) if ( rcode == 0 ) rcode = trcode trcode = rcode rcode = nf_inq_varid ( ncid, "PH", id_var ) if ( rcode == 0 ) rcode = trcode trcode = rcode rcode = nf_inq_varid ( ncid, "PHB", id_var ) if ( rcode == 0 ) rcode = trcode trcode = rcode rcode = nf_inq_varid ( ncid, "T", id_var ) if ( rcode == 0 ) rcode = trcode trcode = rcode rcode = nf_inq_varid ( ncid, "QVAPOR", id_var ) if ( rcode == 0 ) rcode = trcode else rcode = nf_inq_varid ( ncid, variables2d(i), id_var ) endif if (rcode .ne. 0) then write(6,*) ' Not in file, remove from list: ',trim(variables2d(i)) variables2d(i) = ' ' desc2d(i) = ' ' missing2d = missing2d+1 endif enddo newi = 0 do i = 1,number_of_3dvar if ( variables3d(i) /= ' ' ) then newi = newi+1 variables3d(newi) = variables3d(i) desc3d(newi) = desc3d(i) endif enddo number_of_3dvar = number_of_3dvar - missing3d newi = 0 do i = 1,number_of_2dvar if ( variables2d(i) /= ' ' ) then newi = newi+1 variables2d(newi) = variables2d(i) desc2d(newi) = desc2d(i) endif enddo number_of_2dvar = number_of_2dvar - missing2d nf_status = nf_close (ncid) call handle_err('Error closing file',nf_status) end subroutine check_all_variables !------------------------------------------------------------------ END MODULE module_wrf_to_grads_util