var/da/da_tools/da_intpsfc_tem.inc

   1 subroutine da_intpsfc_tem (val, ho, po, to, height, pressure, temp, kts, kte)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: Correct temperature between two levels.
   5    !
   6    ! Reference:
   7    ! ---------
   8    ! The use of surface observations in four dimensional data assimilation
   9    !  Using a mesoscale model.
  10    !  Ruggiero et al., 1996, Monthly Weather Review, Volume 124, 1018-1033
  11    !
  12    !----------------------------------------------------------------------------
  13
  14    implicit none
  15
  16    real,    intent (out) :: val
  17    integer, intent (in)  :: kts, kte
  18    real,    intent (in)  :: ho, po, to
  19    real,    intent (in)  :: height(kts:kte)
  20    real,    intent (in)  :: pressure(kts:kte)
  21    real,    intent (in)  :: temp(kts:kte)
  22
  23    real    :: prs_mb(kts:kte)
  24    ! calculated but never used
  25    !      real    :: dth_12, dth_21, dth_sfc, dth_obs
  26    !     real    :: dhe_12, dhe_21, dhe_sfc1, dhe_obs1, dhe_sfc2, dhe_obs2
  27    real    :: dth_21, dth_sfc, dth_obs
  28    real    :: dhe_12, dhe_sfc1, dhe_obs1, dhe_sfc2, dhe_obs2
  29    real    :: th_100mb, th_200mb, th_obs, th_sfc
  30    real    :: th_obs_int, th_sfc_int
  31    real    :: pdif, rcp
  32    integer :: k_100mb, k_200mb, jk
  33    integer :: inc_100mb, inc_200mb
  34
  35    if (trace_use_dull) call da_trace_entry("da_intpsfc_tem")
  36
  37    rcp = gas_constant/cp
  38
  39    ! 1.Find levels: model surface + 100hpa and model surface + 200hpa ar obs loc
  40    ! ===========================================================================
  41
  42    ! 1.1 Convert model pressure profile from Pa to hPa
  43
  44    prs_mb = pressure / 100.0
  45
  46    ! 1.2  Find levels surface + 100hPA
  47
  48    inc_100mb = 100.0
  49    k_100mb   = kts
  50
  51    do jk =  kts+1, kte
  52       pdif = prs_mb (kts) - prs_mb (jk)
  53       if (pdif .GE. inc_100mb) then
  54          k_100mb = jk
  55          exit
  56       end if
  57    end do
  58
  59    ! 1.2  Find levels surface + 200hPA
  60
  61    inc_200mb = 200.0
  62    k_200mb   = kts
  63
  64    do jk =  kts+1, kte
  65       pdif = prs_mb (kts) - prs_mb (jk)
  66       if (pdif .GE. inc_200mb) then
  67          k_200mb = jk
  68          exit
  69       end if
  70    end do
  71
  72    ! 1.3  Check consistency
  73
  74    if ((k_100mb .LE. kts) .OR. (k_200mb .LE. kts) .OR. &
  75         (k_200mb .LT. k_100mb)) then
  76       write (unit=message(1),fmt='(A)') ' Cannot find sfc + 200hPa and sfc + 100hPa'
  77       write (unit=message(2),fmt='(A,I2,A,F10.3)') ' P (',k_200mb,') = ',prs_mb (k_200mb)
  78       write (unit=message(3),fmt='(A,I2,A,F10.3)') ' P (',k_100mb,') = ',prs_mb (k_100mb)
  79       write (unit=message(4),fmt='(A,F10.3)')      ' P_SFC  = ',         prs_mb (kts)
  80       call da_warning(__FILE__,__LINE__,message(1:4))
  81       val = missing_r
  82       if (trace_use_dull) call da_trace_exit("da_intpsfc_tem")
  83       return
  84    end if
  85
  86    ! 2.  potential temperature
  87    ! =========================
  88
  89    ! 2.1 Potential temperature at 100hPa above model surface
  90
  91    th_100mb = temp (k_100mb) * (1000.0 / prs_mb (k_100mb))**rcp
  92
  93    ! 2.2 Potential temperature at 200hPa above model surface
  94
  95    th_200mb = temp (K_200mb) * (1000.0 / prs_mb (k_200mb))**rcp
  96
  97    ! 2.3 Potential temperature at observation location
  98
  99    th_obs   = to * (1000.0 / (po/100.0)) ** rcp
 100
 101
 102    ! 3.  lapse rate between surface+200hpa and surface+100hpa
 103    ! =========================================================
 104
 105    ! 3.1 Potential temperature increment
 106
 107    dth_21 = th_100mb - th_200mb
 108    ! never used
 109    ! dth_12 = th_200mb - th_100mb
 110
 111    ! 3.1 Height increments
 112
 113    ! never used
 114    ! dhe_21   = height (k_100mb)- height (k_200mb)
 115    dhe_sfc1 = height (k_100mb)- height (kts)
 116    dhe_obs1 = height (k_100mb)- ho
 117
 118    dhe_12   = height (k_200mb)- height (k_100mb)
 119    dhe_sfc2 = height (k_200mb)- height (kts)
 120    dhe_obs2 = height (k_200mb)- ho
 121
 122    ! 3.2 Extrapolated potential temperature at model surface and observation loc
 123
 124    th_sfc_int = th_100mb + (dth_21/dhe_12) * dhe_sfc1
 125    th_obs_int = th_100mb + (dth_21/dhe_12) * dhe_obs1
 126
 127
 128    ! 4.  Bring temperature onto model surface
 129    ! ========================================
 130
 131    ! 4.1 Difference at observation locations
 132
 133    dth_obs = th_obs_int - th_obs
 134
 135    ! 4.2 Difference at model surface
 136
 137    dth_sfc = (dhe_sfc2/dhe_obs2) * dth_obs
 138
 139    ! 4.3 Potentiel temperature brought to model surface
 140
 141    th_sfc = th_sfc_int - dth_sfc
 142
 143    ! 4.3 Corresponding Temperature
 144
 145    val  = th_sfc * (prs_mb (kts) / 1000.0)**rcp
 146
 147    if (trace_use_dull) call da_trace_exit("da_intpsfc_tem")
 148
 149 end subroutine da_intpsfc_tem
 150
 151