var/da/da_gpseph/da_transform_xtoy_gpseph.inc

   1 subroutine da_transform_xtoy_gpseph (iv, y)
   2
   3    !-------------------------------------------------------------------------
   4    ! Purpose: y = H(grid%xa)
   5    !-------------------------------------------------------------------------
   6
   7    implicit none
   8
   9    type (iv_type), intent(in)    :: iv       ! Innovation vector (O-B).
  10    type (y_type),  intent(inout) :: y        ! y = h (grid%xa)
  11
  12    integer :: is, ie, js, je, ks, ke
  13    real, dimension(kds:kde)   :: mean_h   ! mean altitude in each level
  14    real, dimension(kds:kde)   :: mdl_z
  15    real, dimension(kds:kde)   :: mdl_ref, mdl_ref9
  16    real, dimension(kds:kde)   :: refm, refm9
  17    real, dimension(3,kds:kde) :: cc, cc9
  18    real, dimension(ids:ide,jds:jde,kds:kde) :: ref_mean_h_tl
  19    integer :: i, j, k, l, m, n, i1, i2, i3, nbot, ntop, nk, nn
  20    real :: step, h, tmp_ref
  21    real,dimension(2) :: w1
  22    real,dimension(2,2) :: w2
  23    integer :: ip1,ip2,ip3,je2
  24    real :: refp
  25
  26    if (trace_use_dull) call da_trace_entry("da_transform_xtoy_gpseph")
  27
  28    is = ids
  29    ie = ide
  30    js = jds
  31    je = jde
  32    ks = kds
  33    ke = kde
  34    nk = kde-kds+1
  35
  36    step = gps_ray_path_step
  37    mean_h(:) = global_h_mean(:) !km
  38
  39    ! calculate ref increment on mean_h from ref increment on model grid
  40    do j=js,je
  41       do i=is,ie
  42
  43          mdl_z(kds:kde)    = global_h(i,j,kds:kde)
  44          mdl_ref(kds:kde)  = global_xa_ref(i,j,kds:kde)*(1./global_ref(i,j,kds:kde))
  45          mdl_ref9(kds:kde) = log(global_ref(i,j,kds:kde))
  46
  47          call da_splinx_lin(nk,mdl_z,mdl_ref9,mdl_ref,cc9,cc,nk,mean_h*1000.,refm9,refm)
  48
  49          ref_mean_h_tl(i,j,kds:kde) = refm(kds:kde)*exp(refm9(kds:kde))
  50
  51       enddo
  52    enddo
  53
  54    do n=iv%info(gpseph)%n1,iv%info(gpseph)%n2
  55
  56       nbot = gps_rays(n)%nbot
  57       ntop = gps_rays(n)%ntop
  58       if ( nbot == 0 .and. ntop == 0 ) cycle
  59
  60       do i=nbot,ntop
  61          if ( iv%gpseph(n)%eph(i)%qc < obs_qc_pointer ) cycle
  62          y%gpseph(n)%eph(i)=0.0
  63
  64          if ( gps_rays(n)%ilocal(i)==1 ) then
  65
  66             je2 = gps_rays(n)%je2(i)
  67             do j=1,je2
  68                do k=-1, 1, 2
  69                   h = gps_rays(n)%ip123(i)%h(k*j)
  70                   if ( h <= mean_h(ke-1) ) then
  71                      i1 = gps_rays(n)%ip123(i)%i1(k*j)
  72                      i2 = gps_rays(n)%ip123(i)%i2(k*j)
  73                      i3 = gps_rays(n)%ip123(i)%i3(k*j)
  74                      w1 = gps_rays(n)%ip123(i)%w1(1:2,k*j)
  75                      w2 = gps_rays(n)%ip123(i)%w2(1:2,1:2,k*j)
  76                      tmp_ref=0.0
  77                      do l=1,2
  78                         do m=1,2
  79                            do nn=1,2
  80                               tmp_ref=tmp_ref+ref_mean_h_tl(i1+l-1,i2+m-1,i3+nn-1)*w2(l,m)*w1(nn)
  81                            end do
  82                         end do
  83                      end do
  84                      y%gpseph(n)%eph(i)=y%gpseph(n)%eph(i)+step*tmp_ref
  85                   end if
  86                end do
  87             end do
  88
  89          else if ( gps_rays(n)%ilocal(i)==0 ) then
  90
  91             ip1 = gps_rays(n)%ip123(i)%i1(0)
  92             ip2 = gps_rays(n)%ip123(i)%i2(0)
  93             ip3 = gps_rays(n)%ip123(i)%i3(0)
  94             w1 = gps_rays(n)%ip123(i)%w1(1:2,0)
  95             w2 = gps_rays(n)%ip123(i)%w2(1:2,1:2,0)
  96
  97             refp=0.
  98             do l=1,2
  99                do m=1,2
 100                   do nn=1,2
 101                      refp=refp+ref_mean_h_tl(ip1+l-1,ip2+m-1,ip3+nn-1)*w2(l,m)*w1(nn)
 102                   end do
 103                end do
 104             end do
 105             y%gpseph(n)%eph(i)=step*refp
 106
 107          end if
 108
 109       end do !nbot-ntop loop
 110    end do !n1-n2 loop
 111
 112    if (trace_use_dull) call da_trace_exit("da_transform_xtoy_gpseph")
 113
 114 end subroutine da_transform_xtoy_gpseph
 115