var/da/da_varbc_tamdar/da_varbc_tamdar_pred.inc

   1    subroutine da_varbc_tamdar_pred(iv, be, ob)
   2
   3    !--------------------------------------!
   4    !  Calculate Predictors in Bias Model  !
   5    !--------------------------------------!
   6
   7    implicit none
   8
   9    type (iv_type), intent(inout) :: iv
  10    type (be_type), intent(inout) :: be
  11    type (y_type),  intent(inout) :: ob
  12
  13    integer                       :: jt_start,size_jt
  14    integer                       :: n,nphase,iflt,iobs,ipred,isn,iqc,iphase
  15    real                          :: dhdt,dtdt,p1,p2,p3,p4
  16
  17    integer, allocatable          :: tdiff(:),osn(:)
  18    real   , allocatable          :: op(:),ot(:),oh(:),mach(:)
  19
  20    character(len=24)             :: otime,stime
  21    character(len=40)             :: stringn
  22
  23
  24    if (trace_use) call da_trace_entry("da_varbc_tamdar_pred")
  25
  26    if (rootproc) write(unit=varbc_tamdar_unit,fmt='(//5X,A/)') &
  27                       'Calculating predictors'
  28
  29    size_jt  = 0
  30    jt_start = be%cv%size_jb + be%cv%size_je + be%cv%size_jp + be%cv%size_js + be%cv%size_jl
  31
  32    nphase = iv%varbc_tamdar%nphase
  33
  34    allocate( op    (iv%varbc_tamdar%nmaxobs) )
  35    allocate( ot    (iv%varbc_tamdar%nmaxobs) )
  36    allocate( oh    (iv%varbc_tamdar%nmaxobs) )
  37    allocate( osn   (iv%varbc_tamdar%nmaxobs) )
  38    allocate( mach  (iv%varbc_tamdar%nmaxobs) )
  39    allocate( tdiff (iv%varbc_tamdar%nmaxobs) )
  40
  41    do iflt = 1, iv%varbc_tamdar%nair
  42       if (iv%varbc_tamdar%nobs_sum(nphase+1,iflt) >= varbc_tamdar_nobsmin) then
  43
  44           n = 0
  45           op(:) = 0.
  46           ot(:) = 0.
  47           oh(:) = 0.
  48           osn(:) = 0
  49           mach(:) = 0.
  50           tdiff(:) = 0
  51
  52           if (iv%varbc_tamdar%nobs(nphase+1,iflt) > 1) then
  53               do iobs = 1, iv%varbc_tamdar%nobs(nphase+1,iflt)
  54
  55                  isn = iv%varbc_tamdar%obs_sn(iobs,nphase+1,iflt)
  56                  iqc = iv%tamdar(isn)%t(1)%qc
  57
  58                  if (iv%tamdar(isn)%p(1) > missing_r .and. iqc >= 0) then
  59                      n = n + 1
  60                      osn(n) = isn
  61                      op(n) = iv%tamdar(isn)%p(1)
  62                      ot(n) = ob%tamdar(isn)%t(1)
  63                      oh(n) = (1-(op(n)/101325.0)**0.190284)*145366.45*0.3048    ! pressure altitude
  64
  65                      if (n == 1) stime = iv%info(tamdar)%date_char(isn)
  66                      otime = iv%info(tamdar)%date_char(isn)
  67                      tdiff(n) = int(da_diff_seconds(otime, stime))
  68
  69                    ! May assign Mach number on the position of 'pw' in ob.ascii
  70                      if (op(n) >= 95000.0) mach(n) = 0.208
  71                      if (op(n) < 95000.0 .and. op(n) >= 90000.0) mach(n) = 0.265
  72                      if (op(n) < 90000.0 .and. op(n) >= 85000.0) mach(n) = 0.288
  73                      if (op(n) < 85000.0 .and. op(n) >= 80000.0) mach(n) = 0.319
  74                      if (op(n) < 80000.0 .and. op(n) >= 75000.0) mach(n) = 0.351
  75                      if (op(n) < 75000.0 .and. op(n) >= 70000.0) mach(n) = 0.387
  76                      if (op(n) < 70000.0 .and. op(n) >= 65000.0) mach(n) = 0.419
  77                      if (op(n) < 65000.0 .and. op(n) >= 60000.0) mach(n) = 0.446
  78                      if (op(n) < 60000.0 .and. op(n) >= 55000.0) mach(n) = 0.474
  79                      if (op(n) < 55000.0 .and. op(n) >= 50000.0) mach(n) = 0.492
  80                      if (op(n) < 50000.0 .and. op(n) >= 45000.0) mach(n) = 0.505
  81                      if (op(n) < 45000.0 .and. op(n) >= 40000.0) mach(n) = 0.529
  82                      if (op(n) < 40000.0 .and. op(n) >= 30000.0) mach(n) = 0.552
  83                      if (op(n) < 30000.0) mach(n) = 0.580
  84                  end if
  85               end do
  86           end if
  87
  88           if (n > 1) then
  89               do iobs = 2, n
  90
  91                  if (tdiff(iobs) == tdiff(iobs-1) .or. abs(tdiff(iobs) - tdiff(iobs-1)) > 360) cycle
  92
  93                  dhdt = (oh(iobs)-oh(iobs-1))/10.0/(tdiff(iobs)-tdiff(iobs-1))
  94                  dtdt = -1.0*(ot(iobs)-ot(iobs-1))/(tdiff(iobs)-tdiff(iobs-1))
  95
  96                  p1 = (dhdt + 2.0) / 5.0
  97                  p2 = 1.0 / (10.0 * mach(iobs))
  98                  p3 = dtdt * 5.0
  99                  p4 = ot(iobs) / 2000.0
 100
 101                  if (dhdt < -0.2) then
 102
 103                      iv%varbc_tamdar%nobs(1,iflt) = iv%varbc_tamdar%nobs(1,iflt) + 1
 104                      iv%varbc_tamdar%obs_sn(iv%varbc_tamdar%nobs(1,iflt),1,iflt) = osn(iobs)
 105
 106                      if (varbc_tamdar_bm == 1) then
 107                          iv%varbc_tamdar%pred(2,1,iflt) = iv%varbc_tamdar%pred(2,1,iflt) + dhdt
 108                      elseif (varbc_tamdar_bm == 2) then
 109                          iv%varbc_tamdar%pred(2,1,iflt) = iv%varbc_tamdar%pred(2,1,iflt) + p1
 110                          iv%varbc_tamdar%pred(3,1,iflt) = iv%varbc_tamdar%pred(3,1,iflt) + p2
 111                          iv%varbc_tamdar%pred(4,1,iflt) = iv%varbc_tamdar%pred(4,1,iflt) + p3
 112                          iv%varbc_tamdar%pred(5,1,iflt) = iv%varbc_tamdar%pred(5,1,iflt) + p4
 113                      end if
 114
 115                  elseif (dhdt > 0.3) then
 116
 117                      iv%varbc_tamdar%nobs(2,iflt) = iv%varbc_tamdar%nobs(2,iflt) + 1
 118                      iv%varbc_tamdar%obs_sn(iv%varbc_tamdar%nobs(2,iflt),2,iflt) = osn(iobs)
 119
 120                      if (varbc_tamdar_bm == 1) then
 121                          iv%varbc_tamdar%pred(2,2,iflt) = iv%varbc_tamdar%pred(2,2,iflt) + dhdt
 122                      elseif (varbc_tamdar_bm == 2) then
 123                          iv%varbc_tamdar%pred(2,2,iflt) = iv%varbc_tamdar%pred(2,2,iflt) + p1
 124                          iv%varbc_tamdar%pred(3,2,iflt) = iv%varbc_tamdar%pred(3,2,iflt) + p2
 125                          iv%varbc_tamdar%pred(4,2,iflt) = iv%varbc_tamdar%pred(4,2,iflt) + p3
 126                          iv%varbc_tamdar%pred(5,2,iflt) = iv%varbc_tamdar%pred(5,2,iflt) + p4
 127                      end if
 128
 129                  else
 130
 131                      iv%varbc_tamdar%nobs(3,iflt) = iv%varbc_tamdar%nobs(3,iflt) + 1
 132                      iv%varbc_tamdar%obs_sn(iv%varbc_tamdar%nobs(3,iflt),3,iflt) = osn(iobs)
 133
 134                      if (varbc_tamdar_bm == 1) then
 135                          iv%varbc_tamdar%pred(2,3,iflt) = iv%varbc_tamdar%pred(2,3,iflt) + dhdt
 136                      elseif (varbc_tamdar_bm == 2) then
 137                          iv%varbc_tamdar%pred(2,3,iflt) = iv%varbc_tamdar%pred(2,3,iflt) + p1
 138                          iv%varbc_tamdar%pred(3,3,iflt) = iv%varbc_tamdar%pred(3,3,iflt) + p2
 139                          iv%varbc_tamdar%pred(4,3,iflt) = iv%varbc_tamdar%pred(4,3,iflt) + p3
 140                          iv%varbc_tamdar%pred(5,3,iflt) = iv%varbc_tamdar%pred(5,3,iflt) + p4
 141                      end if
 142
 143                  end if
 144
 145               end do
 146           end if
 147
 148           do iphase = 1, nphase
 149
 150              iv%varbc_tamdar%nobs_sum(iphase,iflt) = wrf_dm_sum_integer(iv%varbc_tamdar%nobs(iphase,iflt))
 151
 152              if (iv%varbc_tamdar%nobs_sum(iphase,iflt) >= varbc_tamdar_nobsmin) then
 153
 154                  do ipred = 1, iv%varbc_tamdar%npred
 155                     iv%varbc_tamdar%pred(ipred,iphase,iflt) = &
 156                        wrf_dm_sum_real(iv%varbc_tamdar%pred(ipred,iphase,iflt)) / &
 157                        iv%varbc_tamdar%nobs_sum(iphase,iflt)
 158                     if (ipred == 1) iv%varbc_tamdar%pred(ipred,iphase,iflt) = varbc_tamdar_pred0
 159                  end do
 160
 161                  if (iv%varbc_tamdar%nobs(iphase,iflt) > 0 .and. iv%varbc_tamdar%ifuse(iphase,iflt) > 0) then
 162                      do ipred = 1, iv%varbc_tamdar%npred
 163                         size_jt = size_jt + 1
 164                         iv%varbc_tamdar%index(ipred,iphase,iflt) = jt_start + size_jt
 165                         iv%varbc_tamdar%vtox(ipred,ipred,iphase,iflt) = 1.0
 166                      end do
 167                  end if
 168
 169               !   if (iv%varbc_tamdar%nobs(iphase,iflt) > 0 .and. iv%varbc_tamdar%ifuse(iphase,iflt) <= 0) then
 170               !       do iobs = 1, iv%varbc_tamdar%nobs(iphase,iflt)
 171               !          isn = iv%varbc_tamdar%obs_sn(iobs,iphase,iflt)
 172               !          if (iv%tamdar(isn)%t(1)%qc >= 0) iv%tamdar(isn)%t(1)%qc = fail_varbc_aircraft
 173               !       end do
 174               !   end if
 175
 176              else
 177
 178                  iv%varbc_tamdar%ifuse(iphase,iflt) = 0
 179
 180                  if (iv%varbc_tamdar%nobs(iphase,iflt) > 0) then
 181                      do iobs = 1, iv%varbc_tamdar%nobs(iphase,iflt)
 182                         isn = iv%varbc_tamdar%obs_sn(iobs,iphase,iflt)
 183                         if (iv%tamdar(isn)%t(1)%qc >= 0) iv%tamdar(isn)%t(1)%qc = fail_varbc_aircraft
 184                      end do
 185                  end if
 186
 187              end if
 188           end do
 189
 190       else
 191
 192           iv%varbc_tamdar%ifuse(:,iflt) = 0
 193
 194           if (iv%varbc_tamdar%nobs(nphase+1,iflt) > 0) then
 195               do iobs = 1, iv%varbc_tamdar%nobs(nphase+1,iflt)
 196                  isn = iv%varbc_tamdar%obs_sn(iobs,nphase+1,iflt)
 197                  if (iv%tamdar(isn)%t(1)%qc >= 0) iv%tamdar(isn)%t(1)%qc = fail_varbc_aircraft
 198               end do
 199           end if
 200
 201           if (rootproc) &
 202               write(unit=varbc_tamdar_unit,fmt='(10X,A,I4,A,I4)') &
 203                    'Measurement from ', iv%varbc_tamdar%tail_id(iflt),' is insufficient: ', &
 204                     iv%varbc_tamdar%nobs_sum(nphase+1,iflt)
 205
 206       end if
 207    end do
 208
 209    be % cv % size_jt = size_jt
 210
 211    if (rootproc .and. ANY(iv%varbc_tamdar%nobs_sum(nphase+1,:) >= varbc_tamdar_nobsmin)) &
 212        write(unit=varbc_tamdar_unit,fmt='(/10X,A)') " ID    Predictors(npred,nphase)  &  Obs Count(nphase)"
 213
 214    write(stringn,'(I3)') iv%varbc_tamdar%npred
 215    stringn = '(10X,I4,3('//trim(ADJUSTL(stringn))//'f8.3),3I5)'
 216    stringn = trim(adjustl(stringn))
 217
 218    do iflt = 1, iv%varbc_tamdar%nair
 219       if (rootproc .and. iv%varbc_tamdar%nobs_sum(nphase+1,iflt) >= varbc_tamdar_nobsmin) then
 220           write(unit=varbc_tamdar_unit,fmt=stringn) &
 221                 iv%varbc_tamdar%tail_id(iflt), &
 222                 iv%varbc_tamdar%pred(1:iv%varbc_tamdar%npred,1:nphase,iflt), &
 223                 iv%varbc_tamdar%nobs_sum(1:nphase,iflt)
 224       end if
 225    end do
 226
 227    deallocate(op,ot,osn,mach,tdiff)
 228
 229    if (trace_use) call da_trace_exit("da_varbc_tamdar_pred")
 230
 231    end subroutine da_varbc_tamdar_pred