Merge remote-tracking branch 'origin/release-v4.6.1'
[WRF.git] / var / da / da_radiance / da_predictor_rttov.inc
blob067cbbfadbd7c5ac9f622c5608eb843e17bf3804
1 subroutine da_predictor_rttov(inst, pred,npred,nlevels,temp,hum,t_skin)
3    implicit none
5    ! temp - RTM level temperatures (k)
6    ! hum  - RTM level moistures    (kg/kg)
7    ! t-skin - model skin temperature (k)
8    ! nlevels - number of RTM levels (1:nlevels)
9    ! pf   - RTM full level pressure  (hpa)
11    ! pred(1) - 1000-300 thickness
12    ! pred(2) - 200-50 thickness
13    ! pred(3) - t_skin
14    ! pred(4) - total column precipitable water
16    integer, intent(in)  :: inst, npred, nlevels
17    real,    intent(in)  :: temp(nlevels), hum(nlevels), t_skin
18    real,    intent(out) :: pred(npred)
20    real, parameter :: kth = gas_constant/gravity
21    real, parameter :: kpc = 100.0/gravity
22    real :: pf(nlevels)
23    real :: tv(nlevels-1), qm(nlevels-1)
24    real :: dlp(nlevels-1), dp(nlevels-1)
25    real    :: add_thk
26    integer :: itmp(1), k
27    integer :: index1000, index300, index200, index50
29    if (trace_use) call da_trace_entry("da_predictor_rttov")
31    pf = coefs(inst)%coef%ref_prfl_p
33    do k = 1, nlevels-1
34       ! full-level to half-level
35       qm(k)=0.5*(hum(k)+hum(k+1))
36    end do
38    dlp(1:nlevels-1) = log(pf(2:nlevels)) - log(pf(1:nlevels-1))
39     dp(1:nlevels-1) =     pf(2:nlevels)  -     pf(1:nlevels-1)
41    ! 0.0 find the pressure level index that is closest to 
42    ! 1000hPa, 300hPa, 200hPa, and 50hPa respectively
44    itmp(1:1) = minloc(abs(pf(:)-1000.0))
45    index1000 = itmp(1)
46    itmp(1:1) = minloc(abs(pf(:)-300.0))
47    index300  = itmp(1)
48    itmp(1:1) = minloc(abs(pf(:)-200.0))
49    index200  = itmp(1)
50    itmp(1:1) = minloc(abs(pf(:)-50.0))
51    index50   = itmp(1)
53    ! 1.0 convert all temperatures to virtual temperatures
54    !     and full-level to half-level
55    ! ----------------------------------------------------
56    do k = 1, nlevels-1
57       tv(k) = 0.5*(temp(k)+temp(k+1))*(1.0+0.608*qm(k))
58    end do
60    if ( (1000.0 > pf(nlevels)) ) then
61       add_thk = kth*( tv(nlevels-1)*(log(1000.0)-log(pf(nlevels))) )  ! approximation
62    else
63       add_thk = 0.0
64    end if
66    ! 2.0 construct averages for nesdis thick layers
67    ! ----------------------------------------------
69    pred(1) = kth*sum( tv(index300+1:index1000)*dlp(index300+1:index1000) ) + add_thk
70    pred(2) = kth*sum( tv(index50+1:index200)*dlp(index50+1:index200) )
71    pred(3) = t_skin
72    pred(4) = kpc*sum( qm(1:nlevels-1)*dp(1:nlevels-1) )
74    if (trace_use) call da_trace_exit("da_predictor_rttov")
76 end subroutine da_predictor_rttov