var/da/da_gpseph/da_mdl_ref_to_eph.inc

   1 SUBROUTINE da_mdl_ref_to_eph (gps_rays, ref_mean_h, model_eph)
   2
   3 !---------------------------------------------------------------------------------------
   4 !  Purpose: calculate model excess phase
   5 !           (Chen et al., 2009, doi:10.3319/TAO.2007.11.29.01(F3C))
   6 !           It is called by da_get_innov_vector_gpseph.
   7 !---------------------------------------------------------------------------------------
   8
   9    implicit none
  10    ! input : gps_rays: gps rays information
  11    !         ref_mean_h: model refractivity on mean height
  12    ! output : model_eph
  13
  14    type(gpsrays_type),                     intent(in)  :: gps_rays
  15    real,dimension(ids:ide,jds:jde,kds:kde),intent(in)  :: ref_mean_h
  16    real,dimension(kds:kde),                intent(out) :: model_eph
  17
  18    real,dimension(kds:kde) :: mean_h
  19    integer :: i, j, k, l, m, nn, i1, i2, i3, nbot, ntop
  20    integer :: is, ie, js, je, ks, ke
  21    integer :: i1l,i2l,i1r,i2r
  22    real :: step, h, tmp_ref
  23    real,dimension(2)   :: w1
  24    real,dimension(2,2) :: w2
  25    integer :: ip1,ip2,ip3,je2
  26    real :: refp
  27
  28    step = gps_ray_path_step
  29    ks = kds
  30    ke = kde
  31    nbot = gps_rays%nbot
  32    ntop = gps_rays%ntop
  33    mean_h = global_h_mean(:) !km
  34
  35    do i = nbot, ntop
  36
  37       ! use the ilocal decided in obs_ref_to_eph for consistency
  38       ! between obs and model eph calculation
  39
  40       if ( gps_rays%ilocal(i) == 0 ) then ! local
  41
  42          ip1 = gps_rays%ip123(i)%i1(0)
  43          ip2 = gps_rays%ip123(i)%i2(0)
  44          ip3 = gps_rays%ip123(i)%i3(0)
  45          w1  = gps_rays%ip123(i)%w1(1:2,0)
  46          w2  = gps_rays%ip123(i)%w2(1:2,1:2,0)
  47          refp=0.
  48          do l=1,2
  49             do m=1,2
  50                do nn=1,2
  51                   refp=refp+ref_mean_h(ip1+l-1,ip2+m-1,ip3+nn-1)*w2(l,m)*w1(nn)
  52                end do
  53             end do
  54          end do
  55
  56          model_eph(i)=step*refp
  57
  58       else ! not local
  59
  60          model_eph(i)=0.
  61          je2 = gps_rays%je2(i)
  62          do j=1,je2
  63             !* calculate S from TP point and integrate to different direction (WRF)
  64             !* S is asymmetric
  65             do k=-1,+1,2
  66                ! transform coordiante from cartesian to sphere coordinate
  67                h = gps_rays%ip123(i)%h(k*j)
  68                if (h <= mean_h(ke-1)) then
  69                   i1 = gps_rays%ip123(i)%i1(k*j)
  70                   i2 = gps_rays%ip123(i)%i2(k*j)
  71                   i3 = gps_rays%ip123(i)%i3(k*j)
  72                   w1 = gps_rays%ip123(i)%w1(1:2,k*j)
  73                   w2 = gps_rays%ip123(i)%w2(1:2,1:2,k*j)
  74                   tmp_ref=0.
  75                   do l=1,2
  76                      do m=1,2
  77                         do nn=1,2
  78                            tmp_ref=tmp_ref+ref_mean_h(i1+l-1,i2+m-1,i3+nn-1)*w2(l,m)*w1(nn)
  79                         end do
  80                      end do
  81                   end do
  82                   model_eph(i)=model_eph(i)+step*tmp_ref
  83                end if
  84             end do !k= -1, 1 loop
  85          end do !je2 loop
  86
  87       end if
  88
  89    end do !kbot, ntop loop
  90
  91 END SUBROUTINE da_mdl_ref_to_eph