module module_tracker implicit none private #if ( HWRF == 1 ) public :: ncep_tracker_center, ncep_tracker_init, update_tracker_post_move real, parameter :: invE=0.36787944117 ! 1/e ! Copied from tracker: real,parameter :: searchrad_6=250.0 ! km - ignore data more than this far from domain center real,parameter :: searchrad_7=200.0 ! km - ignore data more than this far from domain center integer, parameter :: maxtp=11 ! number of tracker parameters real, parameter :: uverrmax = 225.0 ! For use in get_uv_guess real, parameter :: ecircum = 40030.2 ! Earth's circumference ! (km) using erad=6371.e3 real, parameter :: rads_vmag=120.0 ! max search radius for wind minimum real, parameter :: err_reg_init=300.0 ! max err at initial time (km) real, parameter :: err_reg_max=225.0 ! max err at other times (km) real, parameter :: errpmax=485.0 ! max stddev of track parameters real, parameter :: errpgro=1.25 ! stddev multiplier real, parameter :: max_wind_search_radius=searchrad_7 ! max radius for vmax search real, parameter :: min_mlsp_search_radius=searchrad_7 ! max radius for pmin search ! Also used: real, parameter :: km2nmi = 0.539957, kn2mps=0.514444, mps2kn=1./kn2mps, pi180=0.01745329251 contains !---------------------------------------------------------------------------------- ! These two simple routines return an N, S, E or W for the ! hemisphere of a latitude or longitude. They are copied from ! module_HIFREQ to avoid a relatively pointless compiler dependency. character(1) function get_lat_ns(lat) ! This could be written simply as merge('N','S',lat>=0) if WRF allowed F95 implicit none ; real lat if(lat>=0) then get_lat_ns='N' else get_lat_ns='S' endif end function get_lat_ns character(1) function get_lon_ew(lon) ! This could be written simply as merge('E','W',lon>=0) if WRF allowed F95 implicit none ; real lon if(lon>=0) then get_lon_ew='E' else get_lon_ew='W' endif end function get_lon_ew subroutine ncep_tracker_init(grid) ! Initialize tracker variables in the grid structure. use module_domain, only: domain implicit none type(domain), intent(inout) :: grid call wrf_message('ncep_tracker_init') grid%track_stderr_m1=-99.9 grid%track_stderr_m2=-99.9 grid%track_stderr_m3=-99.9 grid%track_n_old=0 grid%track_old_lon=0 grid%track_old_lat=0 grid%track_old_ntsd=0 grid%tracker_angle=0 grid%tracker_fixlon=-999.0 grid%tracker_fixlat=-999.0 grid%tracker_ifix=-99 grid%tracker_jfix=-99 grid%tracker_havefix=.false. grid%tracker_gave_up=.false. grid%tracker_pmin=-99999. grid%tracker_vmax=-99. grid%tracker_rmw=-99. grid%track_have_guess=.false. grid%track_guess_lat=-999.0 grid%track_guess_lon=-999.0 end subroutine ncep_tracker_init subroutine ncep_tracker_center(grid) ! Top-level entry to the inline ncep tracker. Finds the center of ! the storm in the specified grid and updates the grid variables. ! Will do nothing and return immediately if ! grid%tracker_gave_up=.true. USE MODULE_DOMAIN, ONLY : domain,get_ijk_from_grid implicit none type(domain), intent(inout) :: grid character*255 :: message integer :: IDS,IDE,JDS,JDE,KDS,KDE integer :: IMS,IME,JMS,JME,KMS,KME integer :: IPS,IPE,JPS,JPE,KPS,KPE CALL get_ijk_from_grid ( grid , & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe ) call ntc_impl(grid, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe ) end subroutine ncep_tracker_center subroutine ntc_impl(grid, & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE) ! This is the main entry point to the tracker. It is most similar ! to the function "tracker" in the GFDL/NCEP vortex tracker. USE MODULE_DOMAIN, ONLY : domain,get_ijk_from_grid #ifdef DM_PARALLEL use module_dm, only: wrf_dm_sum_real #endif implicit none logical, external :: wrf_dm_on_monitor type(domain), intent(inout) :: grid integer, intent(in) :: IDS,IDE,JDS,JDE,KDS,KDE integer, intent(in) :: IMS,IME,JMS,JME,KMS,KME integer, intent(in) :: IPS,IPE,JPS,JPE,KPS,KPE real :: dxdymean, sum integer :: i,j, iweights,ip integer :: iguess, jguess ! first guess location real :: latguess, longuess ! same, but in lat & lon integer :: iuvguess, juvguess ! "second guess" location using everything except wind maxima real :: srsq integer :: ifinal,jfinal real :: latfinal,lonfinal integer :: ierr integer :: icen(maxtp), jcen(maxtp) ! center locations for each parameter real :: loncen(maxtp), latcen(maxtp) ! lat, lon locations in degrees logical :: calcparm(maxtp) ! do we have a valid center location for this parameter? real :: max_wind,min_pres ! for ATCF output real :: rcen(maxtp) ! center value (max wind, min mslp, etc.) character*255 :: message logical :: north_hemi ! true = northern hemisphere logical :: have_guess ! first guess is available real :: guessdist,guessdeg ! first guess distance to nearest point on grid real :: latnear, lonnear ! nearest point in grid to first guess ! icen,jcen: Same meaning as clon, clat in tracker, but uses i and ! j indexes of the center instead of lat/lon. Tracker comment: ! Holds the coordinates for the center positions for ! all storms at all times for all parameters. ! (max_#_storms, max_fcst_times, max_#_parms). ! For the third position (max_#_parms), here they are: ! 1: Relative vorticity at 850 mb ! 2: Relative vorticity at 700 mb ! 3: Vector wind magnitude at 850 mb ! 4: NOT CURRENTLY USED ! 5: Vector wind magnitude at 700 mb ! 6: NOT CURRENTLY USED ! 7: Geopotential height at 850 mb ! 8: Geopotential height at 700 mb ! 9: Mean Sea Level Pressure ! 10: Vector wind magnitude at 10 m ! 11: Relative vorticity at 10 m call wrf_message('ncep_tracker_center') ! Initialize center information to invalid values for all centers: icen=-99 jcen=-99 latcen=9e9 loncen=9e9 rcen=9e9 calcparm=.false. if(grid%vortex_tracker==6) then srsq=searchrad_6*searchrad_6*1e6 else srsq=searchrad_7*searchrad_7*1e6 endif ! Get the first guess from the prior nest motion timestep: have_guess=grid%track_have_guess if(have_guess) then ! We have a first guess center. We have to translate it to gridpoint space. longuess=grid%track_guess_lon latguess=grid%track_guess_lat call get_nearest_lonlat(grid,iguess,jguess,ierr,longuess,latguess, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe, lonnear, latnear) if(ierr==0) then call calcdist(longuess,latguess, lonnear,latnear, guessdist,guessdeg) if(guessdist*1e3>3*grid%dy) then 108 format('WARNING: guess lon=',F0.3,',lat=',F0.3, & ' too far (',F0.3,'km) from nearest point lon=',F0.3,',lat=',F0.3, & '. Will use domain center as first guess.') write(message,108) grid%track_guess_lon,grid%track_guess_lat, & guessdist,lonnear,latnear call wrf_message(message) have_guess=.false. ! indicate that the first guess is unusable else latguess=latnear longuess=lonnear endif else have_guess=.false. ! indicate that the first guess is unusable. 109 format('WARNING: guess lon=',F0.3,',lat=',F0.3, & ' does not exist in this domain. Will use domain center as first guess.') write(message,109) grid%track_guess_lon,grid%track_guess_lat call wrf_message(message) endif endif ! If we could not get the first guess from the prior nest motion ! timestep, then use the default first guess: the domain center. if(grid%vortex_tracker==6 .or. .not.have_guess) then ! vt=6: hard coded first-guess center is domain center: ! vt=7: first guess comes from prior timestep ! Initial first guess is domain center. ! Backup first guess is domain center if first guess is unusable. iguess=ide/2 jguess=jde/2 if(grid%vortex_tracker==7) then call wrf_message('Using domain center as first guess since no valid first guess is available.') endif call get_lonlat(grid,iguess,jguess,longuess,latguess,ierr, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe) if(ierr/=0) then call wrf_error_fatal("ERROR: center of domain is not inside the domain") endif have_guess=.true. endif if(.not.have_guess) then call wrf_error_fatal("INTERNAL ERROR: No first guess is available (should never happen).") endif north_hemi = latguess>0.0 ! Get the mean V-to-H point-to-point distance: sum=0 do j=jps,min(jde-1,jpe) do i=ips,min(ide-1,ipe) sum=sum+grid%dx_nmm(i,j) enddo enddo #ifdef DM_PARALLEL sum=wrf_dm_sum_real(sum) #endif dxdymean=0.5*(grid%dy_nmm + sum/( (ide-ids) * (jde-jds) ))/1000.0 33 format ('dxdymean=',F0.3,' dx=',F0.3,' dy=',F0.3,' sum=',F0.3,' count=',I0) !write(0,33) dxdymean,grid%dx_nmm((ips+ipe)/2,(jps+jpe)/2),grid%dy_nmm, & ! sum,(ide-ids) * (jde-jds) ! Find the centers of all fields except the wind minima: call find_center(grid,grid%p850rv,grid%sp850rv,srsq, & icen(1),jcen(1),rcen(1),calcparm(1),loncen(1),latcen(1),dxdymean,'zeta', & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE, north_hemi=north_hemi) call find_center(grid,grid%p700rv,grid%sp700rv,srsq, & icen(2),jcen(2),rcen(2),calcparm(2),loncen(2),latcen(2),dxdymean,'zeta', & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE, north_hemi=north_hemi) call find_center(grid,grid%p850z,grid%sp850z,srsq, & icen(7),jcen(7),rcen(7),calcparm(7),loncen(7),latcen(7),dxdymean,'hgt', & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE) call find_center(grid,grid%p700z,grid%sp700z,srsq, & icen(8),jcen(8),rcen(8),calcparm(8),loncen(8),latcen(8),dxdymean,'hgt', & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE) call find_center(grid,grid%membrane_mslp,grid%smslp,srsq, & icen(9),jcen(9),rcen(9),calcparm(9),loncen(9),latcen(9),dxdymean,'slp', & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE) call find_center(grid,grid%m10rv,grid%sm10rv,srsq, & icen(11),jcen(11),rcen(11),calcparm(11),loncen(11),latcen(11),dxdymean,'zeta', & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE, north_hemi=north_hemi) ! Get a guess center location for the wind minimum searches: call get_uv_guess(grid,icen,jcen,loncen,latcen,calcparm, & iguess,jguess,longuess,latguess,iuvguess,juvguess, & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE) ! Find wind minima. Requires a first guess center: windmin: if(grid%vortex_tracker==6) then call find_center(grid,grid%p850wind,grid%sp850wind,srsq, & icen(3),jcen(3),rcen(3),calcparm(3),loncen(3),latcen(3),dxdymean,'wind', & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE, & iuvguess=iuvguess, juvguess=juvguess) call find_center(grid,grid%p700wind,grid%sp700wind,srsq, & icen(5),jcen(5),rcen(5),calcparm(5),loncen(5),latcen(5),dxdymean,'wind', & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE, & iuvguess=iuvguess, juvguess=juvguess) call find_center(grid,grid%m10wind,grid%sm10wind,srsq, & icen(10),jcen(10),rcen(10),calcparm(10),loncen(10),latcen(10),dxdymean,'wind', & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE, & iuvguess=iuvguess, juvguess=juvguess) else call get_uv_center(grid,grid%p850wind, & icen(3),jcen(3),rcen(3),calcparm(3),loncen(3),latcen(3),dxdymean,'wind', & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE, & iuvguess=iuvguess, juvguess=juvguess) call get_uv_center(grid,grid%p700wind, & icen(5),jcen(5),rcen(5),calcparm(5),loncen(5),latcen(5),dxdymean,'wind', & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE, & iuvguess=iuvguess, juvguess=juvguess) call get_uv_center(grid,grid%m10wind, & icen(10),jcen(10),rcen(10),calcparm(10),loncen(10),latcen(10),dxdymean,'wind', & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE, & iuvguess=iuvguess, juvguess=juvguess) endif windmin ! Get a final guess center location: call fixcenter(grid,icen,jcen,calcparm,loncen,latcen, & iguess,jguess,longuess,latguess, & ifinal,jfinal,lonfinal,latfinal, & north_hemi, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe) grid%tracker_fixes=0 do ip=1,maxtp if(calcparm(ip)) then 300 format('Parameter ',I0,': i=',I0,' j=',I0,' lon=',F0.2,' lat=',F0.2) !write(0,300) ip,icen(ip),jcen(ip),loncen(ip),latcen(ip) if(icen(ip)>=ips .and. icen(ip)<=ipe & .and. jcen(ip)>=jps .and. jcen(ip)<=jpe) then grid%tracker_fixes(icen(ip),jcen(ip))=ip endif else 301 format('Parameter ',I0,' invalid') !write(0,301) ip endif enddo if(iguess>=ips .and. iguess<=ipe .and. jguess>=jps .and. jguess<=jpe) then grid%tracker_fixes(iguess,jguess)=-1 201 format('First guess: i=',I0,' j=',I0,' lon=',F0.2,' lat=',F0.2) !write(0,201) iguess,jguess,longuess,latguess endif if(iuvguess>=ips .and. iuvguess<=ipe .and. juvguess>=jps .and. juvguess<=jpe) then grid%tracker_fixes(iuvguess,juvguess)=-2 202 format('UV guess: i=',I0,' j=',I0) !write(0,202) iguess,jguess endif 1000 format('Back with final lat/lon at i=',I0,' j=',I0,' lon=',F0.3,' lat=',F0.3) !write(0,1000) ifinal,jfinal,lonfinal,latfinal if(ifinal>=ips .and. ifinal<=ipe .and. jfinal>=jps .and. jfinal<=jpe) then grid%tracker_fixes(ifinal,jfinal)=-3 203 format('Final fix: i=',I0,' j=',I0,' lon=',F0.2,' lat=',F0.2) !write(0,201) ifinal,jfinal,lonfinal,latfinal endif call get_tracker_distsq(grid, & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE) call get_wind_pres_intensity(grid, & grid%tracker_pmin,grid%tracker_vmax,grid%tracker_rmw, & max_wind_search_radius, min_mlsp_search_radius, & lonfinal,latfinal, & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE) #ifdef DM_PARALLEL if(wrf_dm_on_monitor()) then #endif call output_partial_atcfunix(grid, & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE) #ifdef DM_PARALLEL endif #endif ! if(grid%vortex_tracker==7) then ! call get_first_ges(grid,iguess,jguess,longuess,latguess, & ! IDS,IDE,JDS,JDE,KDS,KDE, & ! IMS,IME,JMS,JME,KMS,KME, & ! IPS,IPE,JPS,JPE,KPS,KPE) ! call store_old_fixes(grid, & ! IDS,IDE,JDS,JDE,KDS,KDE, & ! IMS,IME,JMS,JME,KMS,KME, & ! IPS,IPE,JPS,JPE,KPS,KPE) ! ! Store the first guess: ! grid%track_have_guess=.true. ! grid%track_guess_lat=latguess ! grid%track_guess_lon=longuess ! 3011 format('First guess: lon=',F0.3,' lat=',F0.3) ! write(message,3011) grid%track_guess_lon,grid%track_guess_lat ! call wrf_debug(1,message) ! endif end subroutine ntc_impl subroutine get_first_ges(grid, & iguess,jguess,longuess,latguess, & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE) ! This replicates the functionality of the tracker get_first_ges ! routine, whose purpose is to analyze the storm and guess where ! it will be at the next nest motion timestep. It does that using ! two different methods, similar to the GFDL/NCEP Tracker's ! methods: ! ! 1. Use the present, and past few, fix locations and extrapolate ! to the next location. ! ! 2. Calculate the mean motion and extrapolate to get the ! location at the next nest motion timestep. ! ! The average of the two results is used. #ifdef DM_PARALLEL use module_dm, only: wrf_dm_maxval_real #endif USE MODULE_DOMAIN, ONLY : domain,get_ijk_from_grid use module_wrf_error, only: wrf_at_debug_level implicit none type(domain), intent(inout) :: grid integer, intent(in) :: IDS,IDE,JDS,JDE,KDS,KDE integer, intent(in) :: IMS,IME,JMS,JME,KMS,KME integer, intent(in) :: IPS,IPE,JPS,JPE,KPS,KPE integer, intent(out) :: iguess,jguess real, intent(out) :: longuess,latguess character*255 message integer :: iold, inew, jold, jnew integer :: ifix,jfix,jrot,irot,ierr, pinky,brain, n, tsum, ntsd_plus_1, i, told real :: motion_grideast, motion_gridnorth, fixdx real :: dxeast,dynorth, xeast, ynorth real :: dxrot, dyrot, tracker_dt, xsum, ysum, ytsum, xtsum, xxsum, yysum, ttsum real :: mx, my, bx, by ! x=mx*t+bx ; y=my*t+by real :: xrot,yrot logical :: have_motion_guess, have_line_guess have_motion_guess=.false. have_line_guess=.false. if(grid%tracker_havefix) then ifix=grid%tracker_ifix jfix=grid%tracker_jfix call mean_motion(grid, motion_grideast, motion_gridnorth, & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE) #ifdef DM_PARALLEL fixdx=0 if(ifix>=ips .and. ifix<=ipe .and. jfix>=jps .and. jfix<=jpe) then fixdx=grid%dx_nmm(ifix,jfix) endif pinky=2 ; brain=308 call wrf_dm_maxval_real(fixdx,pinky,brain) #else fixdx=grid%dx_nmm(ifix,jfix) #endif ! Rotated east and north motion in gridpoints per second, on the combined H+V grid: tracker_dt=grid%dt*grid%nphs*grid%ntrack dxeast = motion_grideast * tracker_dt / fixdx dynorth = motion_gridnorth * tracker_dt / grid%dy_nmm ! Combine the H & V coordinate systems and rotate 45 degrees. ! This puts the H points on a rectangular grid. Add storm motion ! to the rotated coordinates and round to nearest H point xeast=ifix*2 ynorth=jfix if(mod(jfix,2)==0) xeast=xeast+1 jrot=nint((xeast+ynorth)/2 + (dxeast+dynorth)/2) irot=nint((ynorth-xeast)/2 + ((jde-1)/2) + (dynorth-dxeast)/2) ! Translate back to usual E grid H points iguess=irot-jrot+((jde-1)/2) jguess=irot+jrot-((jde-1)/2) if(mod(jguess,2)==0) then iguess=(iguess-1)/2 else iguess=iguess/2 endif ! This last step should not be necessary but done just in case: have_motion_guess = .not.(iguesside*3/4 .or. jguessjde*3/4) !print *,'got have_motion_guess=',have_motion_guess endif if(.not.have_motion_guess) then ! Could not find the storm, so give the domain center as the ! next first guess location. iguess=ide/2 jguess=jde/2 !print *,'cannot find storm, so using domain center for motion guess' endif if(grid%track_n_old>0) then !print *,'line guess: have old' n=1 call to_rot45_grid(grid%tracker_ifix,grid%tracker_jfix,jde,xrot,yrot) xsum=xrot ysum=yrot tsum=grid%ntsd xtsum=xsum*tsum xxsum=xsum*xsum yysum=ysum*ysum ytsum=ysum*tsum ttsum=tsum*tsum do i=1,grid%track_n_old call get_nearest_lonlat(grid,iold,jold,ierr, & grid%track_old_lon(i),grid%track_old_lat(i), & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe) if(ierr==0) then !print *,'insert: i=',iold,' j=',jold,' lon=',grid%track_old_lon(i),' lat=',grid%track_old_lat(i),' t=',grid%track_old_ntsd(i) call to_rot45_grid(iold,jold,jde,xrot,yrot) n=n+1 xsum=xsum+xrot ysum=ysum+yrot told=grid%track_old_ntsd(i) tsum=tsum+told xtsum=xtsum+xrot*told xxsum=xxsum+xrot*xrot ytsum=ytsum+yrot*told yysum=xxsum+yrot*yrot ttsum=ttsum+told*told endif enddo !print *,'line guess: n=',n if(n>1) then ntsd_plus_1 = grid%ntsd + grid%ntrack*grid%nphs mx=(xtsum-(xsum*tsum)/real(n))/(ttsum-(tsum*tsum)/real(n)) my=(ytsum-(ysum*tsum)/real(n))/(ttsum-(tsum*tsum)/real(n)) bx=(xsum-mx*tsum)/real(n) by=(ysum-my*tsum)/real(n) !print *,'mx=',mx,' my=',my,' bx=',bx,' by=',by,' t+1=',ntsd_plus_1 xrot=nint(mx*ntsd_plus_1+bx) yrot=nint(my*ntsd_plus_1+by) call from_rot45_grid(inew,jnew,jde,xrot,yrot) !print *,'inew=',inew,' jnew=',jnew,' xrot=',xrot,' yrot=',yrot have_line_guess=.not.(inewide*3/4 & .or. jnewjde*3/4) else have_line_guess=.false. endif endif print_locs: if(wrf_at_debug_level(2)) then call get_lonlat(grid,iguess,jguess,longuess,latguess,ierr, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe) if(ierr==0) then if(have_motion_guess) then 3088 format('Motion Guess: lon=',F0.3,' lat=',F0.3) write(message,3088) longuess,latguess call wrf_debug(2,message) else 3089 format('Motion Guess failed; use domain center: lon=',F0.3,' lat=',F0.3) write(message,3089) longuess,latguess call wrf_debug(2,message) endif else 3090 format('Motion guess ierr=',I0) write(message,3090) ierr call wrf_debug(2,message) endif if(have_line_guess) then call get_lonlat(grid,inew,jnew,longuess,latguess,ierr, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe) if(ierr==0) then 3091 format('Line guess: lon=',F0.3,' lat=',F0.3) write(message,3091) longuess,latguess call wrf_debug(2,message) else 3092 format('Line guess ierr=',I0) write(message,3092) ierr call wrf_debug(2,message) endif endif endif print_locs if(have_line_guess) then if(have_motion_guess) then call wrf_debug(1,'get_first_ges: have MOTION and LINE guesses') iguess=(iguess+inew)/2 jguess=(jguess+jnew)/2 else call wrf_debug(1,'get_first_ges: have LINE guess only') iguess=inew jguess=jnew endif elseif(have_motion_guess) then call wrf_debug(1,'get_first_ges: have MOTION guess only') else call wrf_debug(1,'get_first_ges: have no guesses; will use domain center') endif ! Now get lats & lons: call get_lonlat(grid,iguess,jguess,longuess,latguess,ierr, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe) if(ierr/=0) then ! Should never get here due to max/min check before. call wrf_error_fatal("ERROR: domain is not inside the domain in get_first_ges (!?)") endif 38 format('First guess: i=',I0,' j=',I0,' lat=',F8.3,' lon=',F8.3) write(message,38) iguess,jguess,latguess,longuess call wrf_message(message) end subroutine get_first_ges subroutine store_old_fixes(grid, & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE) ! This stores old fix locations for later use in the get_first_ges ! routine's line of best fit. USE MODULE_DOMAIN, ONLY : domain,get_ijk_from_grid implicit none type(domain), intent(inout) :: grid integer, intent(in) :: IDS,IDE,JDS,JDE,KDS,KDE integer, intent(in) :: IMS,IME,JMS,JME,KMS,KME integer, intent(in) :: IPS,IPE,JPS,JPE,KPS,KPE integer i if(grid%tracker_havefix) then !print *,'in store old, have fix' if(grid%track_n_old>0) then !print *,'in store old, shifting old' do i=1,grid%num_old_fixes-1 grid%track_old_lon(i+1)=grid%track_old_lon(i) grid%track_old_lat(i+1)=grid%track_old_lat(i) grid%track_old_ntsd(i+1)=grid%track_old_ntsd(i) enddo endif grid%track_old_lon(1)=grid%tracker_fixlon grid%track_old_lat(1)=grid%tracker_fixlat grid%track_old_ntsd(1)=grid%ntsd grid%track_n_old=min(grid%num_old_fixes,grid%track_n_old+1) !print *,'in store old, now have ',grid%track_n_old endif end subroutine store_old_fixes subroutine to_rot45_grid(i,j,jde,x,y) implicit none integer, intent(in) :: i,j,jde real, intent(inout) :: x,y real :: a,b a=i*2 b=j if(mod(j,2)==0) a=a+1 x=(a+b)/2 y=(b-a)/2+((jde-1)/2) end subroutine to_rot45_grid subroutine from_rot45_grid(i,j,jde,x,y) implicit none integer, intent(inout) :: i,j integer, intent(in) :: jde real, intent(in) :: x,y i=x-y+((jde-1)/2) j=x+y-((jde-1)/2) if(mod(j,2)==0) then i=(i-1)/2 else i=i/2 endif end subroutine from_rot45_grid subroutine get_nearest_lonlat(grid,iloc,jloc,ierr,lon,lat, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe, & lonnear, latnear) ! Finds the nearest point in the domain to the specified lon,lat ! location. #ifdef DM_PARALLEL use module_dm, only: wrf_dm_minloc_real #endif USE MODULE_DOMAIN, ONLY : domain,get_ijk_from_grid implicit none type(domain), intent(inout) :: grid integer, intent(in) :: IDS,IDE,JDS,JDE,KDS,KDE integer, intent(in) :: IMS,IME,JMS,JME,KMS,KME integer, intent(in) :: IPS,IPE,JPS,JPE,KPS,KPE integer, intent(out) :: iloc,jloc,ierr real, intent(in) :: lon,lat real :: dx,dy,d,dmin, zdummy, latmin,lonmin integer :: i,j,imin,jmin real, intent(out), optional :: latnear, lonnear zdummy=42 dmin=9e9 imin=-99 jmin=-99 latmin=9e9 lonmin=9e9 ierr=0 do j=jps,min(jde-1,jpe) do i=ips,min(ide-1,ipe) dy=abs(lat-grid%hlat(i,j)) dx=abs(mod(3600.+180.+(lon-grid%hlon(i,j)),360.)-180.) d=dx*dx+dy*dy if(dlocalextreme) then localextreme=windsq locali=i localj=j endif endif enddo enddo if(localextreme>0) localextreme=sqrt(localextreme) globalextreme=localextreme globali=locali globalj=localj #ifdef DM_PARALLEL call wrf_dm_maxval_real(globalextreme,globali,globalj) #endif call get_lonlat(grid,globali,globalj,globallon,globallat,ierr, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte) if(ierr/=0) then call wrf_message("WARNING: Unable to find location of wind maximum.") rmw=-99 else call calcdist(clon,clat,globallon,globallat,rmw,degrees) end if ! Get the guess location for the next time: max_wind=globalextreme if(globali<0 .or. globalj<0) then call wrf_message("WARNING: No wind values found that were greater than -9*10^9.") min_mslp=-999 endif end subroutine get_wind_pres_intensity subroutine mean_motion(grid,motion_grideast,motion_gridnorth, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte) ! This calculates the mean motion of the storm by calculating the ! average wind vector at 850, 700 and 500 mbars. #ifdef DM_PARALLEL use module_dm, only: wrf_dm_sum_real8, wrf_dm_sum_integer #endif use module_wrf_error USE MODULE_DOMAIN, ONLY : domain, domain_clock_get implicit none integer, intent(in) :: & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte type(domain), intent(in) :: grid real, intent(out) :: motion_grideast,motion_gridnorth integer :: count,i,j,ierr real :: distsq, dist real*8 :: e,n e=0 ; n=0 ; count=0 ! east sum, north sum, count dist = min(grid%tracker_edge_dist, max(50e3, 3e3*grid%tracker_rmw)) distsq = dist * dist ! print *,'motion search radius (m) = ',dist ! print *,' considered edge dist = ',grid%tracker_edge_dist ! print *,' considered 3e3*rmw = ',3e3*grid%tracker_rmw ! print *,' considered 50e3.' do j=jts,min(jte,jde-1) do i=its,min(ite,ide-1) if(grid%tracker_distsq(i,j) 2 .and. ip < 7) .or. ip == 10) then cycle ! because 3-6 are for 850 & 700 u & v and 10 is ! for surface wind magnitude. elseif(calcparm(ip)) then call calcdist (longuess,latguess,loncen(ip),latcen(ip),dist,degrees) if(dist90.) then ylat1=180.-ylat1 xlon1=mod(xlon1+360.,360.)-180. elseif(ylat1<-90.) then ylat1=-180. - ylat1 xlon1=mod(xlon1+360.,360.)-180. endif end subroutine clean_lon_lat subroutine calcdist(rlonb,rlatb,rlonc,rlatc,xdist,degrees) ! Copied from gettrk_main.f ! ! ABSTRACT: This subroutine computes the distance between two ! lat/lon points by using spherical coordinates to ! calculate the great circle distance between the points. ! Figure out the angle (a) between pt.B and pt.C, ! N. Pole then figure out how much of a % of a great ! x circle distance that angle represents. ! / \ ! b/ \ cos(a) = (cos b)(cos c) + (sin b)(sin c)(cos A) ! / \ . ! pt./<--A-->\c NOTE: The latitude arguments passed to the ! B / \ subr are the actual lat vals, but in ! \ the calculation we use 90-lat. ! a \ . ! \pt. NOTE: You may get strange results if you: ! C (1) use positive values for SH lats AND ! you try computing distances across the ! equator, or (2) use lon values of 0 to ! -180 for WH lons AND you try computing ! distances across the 180E meridian. ! ! NOTE: In the diagram above, (a) is the angle between pt. B and ! pt. C (with pt. x as the vertex), and (A) is the difference in ! longitude (in degrees, absolute value) between pt. B and pt. C. ! ! !!! NOTE !!! -- THE PARAMETER ecircum IS DEFINED (AS OF THE ! ORIGINAL WRITING OF THIS SYSTEM) IN KM, NOT M, SO BE AWARE THAT ! THE DISTANCE RETURNED FROM THIS SUBROUTINE IS ALSO IN KM. ! implicit none real, intent(inout) :: degrees real, intent(out) :: xdist real, intent(in) :: rlonb,rlatb,rlonc,rlatc real, parameter :: dtr = 0.0174532925199433 real :: distlatb,distlatc,pole,difflon,cosanga,circ_fract ! if (rlatb < 0.0 .or. rlatc < 0.0) then pole = -90. else pole = 90. endif ! distlatb = (pole - rlatb) * dtr distlatc = (pole - rlatc) * dtr difflon = abs( (rlonb - rlonc)*dtr ) ! cosanga = ( cos(distlatb) * cos(distlatc) + & sin(distlatb) * sin(distlatc) * cos(difflon)) ! This next check of cosanga is needed since I have had ACOS crash ! when calculating the distance between 2 identical points (should ! = 0), but the input for ACOS was just slightly over 1 ! (e.g., 1.00000000007), due to (I'm guessing) rounding errors. if (cosanga > 1.0) then cosanga = 1.0 endif degrees = acos(cosanga) / dtr circ_fract = degrees / 360. xdist = circ_fract * ecircum ! ! NOTE: whether this subroutine returns the value of the distance ! in km or m depends on the scale of the parameter ecircum. ! At the original writing of this subroutine (7/97), ecircum ! was given in km. ! return end subroutine calcdist ! subroutine get_lonlat(grid,iguess,jguess,longuess,latguess, & ! ids,ide, jds,jde, kds,kde, & ! ims,ime, jms,jme, kms,kme, & ! ips,ipe, jps,jpe, kps,kpe) ! ! Returns the latitude (latguess) and longitude (longuess) of the ! ! specified location (iguess,jguess) in the specified grid. ! USE MODULE_DOMAIN, ONLY : domain,get_ijk_from_grid ! USE MODULE_DM, ONLY: wrf_dm_at_ij_real ! implicit none ! integer, intent(in) :: & ! ids,ide, jds,jde, kds,kde, & ! ims,ime, jms,jme, kms,kme, & ! ips,ipe, jps,jpe, kps,kpe ! type(domain), intent(inout) :: grid ! integer, intent(in) :: iguess,jguess ! real, intent(inout) :: longuess,latguess ! call wrf_dm_at_ij_real(grid,iguess,jguess,ims,ime, jms,jme, & ! longuess,grid%hlon, & ! val2=latguess,field2=grid%hlat) ! end subroutine get_lonlat subroutine get_lonlat(grid,iguess,jguess,longuess,latguess,ierr, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe) ! Returns the latitude (latguess) and longitude (longuess) of the ! specified location (iguess,jguess) in the specified grid. USE MODULE_DOMAIN, ONLY : domain,get_ijk_from_grid USE MODULE_DM, ONLY: wrf_dm_maxloc_real implicit none integer, intent(in) :: & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe integer, intent(out) :: ierr type(domain), intent(inout) :: grid integer, intent(in) :: iguess,jguess real, intent(inout) :: longuess,latguess real :: weight,zjunk integer :: itemp,jtemp ierr=0 zjunk=1 if(iguess>=ips .and. iguess<=ipe .and. jguess>=jps .and. jguess<=jpe) then weight=1 longuess=grid%hlon(iguess,jguess) latguess=grid%hlat(iguess,jguess) itemp=iguess jtemp=jguess else weight=0 longuess=-999.9 latguess=-999.9 itemp=-99 jtemp=-99 endif #ifdef DM_PARALLEL call wrf_dm_maxloc_real(weight,latguess,longuess,zjunk,itemp,jtemp) #endif if(itemp==-99 .and. jtemp==-99) then ierr=95 endif end subroutine get_lonlat subroutine update_tracker_post_move(grid) ! This updates the tracker i/j fix location and square of the ! distance to the tracker center after a nest move. USE MODULE_DOMAIN, ONLY : domain,get_ijk_from_grid type(domain), intent(inout) :: grid integer :: ierr, & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE ! Get the grid bounds: CALL get_ijk_from_grid ( grid , & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe ) ! Get the i/j center location from the fix location: ierr=0 call get_nearest_lonlat(grid,grid%tracker_ifix,grid%tracker_jfix, & ierr,grid%tracker_fixlon,grid%tracker_fixlat, & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE) ! Get the square of the approximate distance to the tracker center ! at all points: if(ierr==0) & call get_tracker_distsq(grid, & IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & IPS,IPE,JPS,JPE,KPS,KPE) end subroutine update_tracker_post_move #endif end module module_tracker #if (HWRF == 1) subroutine nmm_med_tracker_post_move(grid) ! This updates the tracker i/j fix location and square of the ! distance to the tracker center after a nest move. use module_tracker, only: update_tracker_post_move use module_domain, only : domain type(domain), intent(inout) :: grid call update_tracker_post_move(grid) end subroutine nmm_med_tracker_post_move #endif