WPS/ungrib/src/datint.F

   1 subroutine datint(fuldates, nful, hstart, ntimes, interval, out_format, prefix)
   2 !                                                                             !
   3 !*****************************************************************************!
   4 !                                                                             !
   5 !   interpolate missing data in time
   6 !    out_format: requested output format
   7 !                                                                             !
   8 !*****************************************************************************!
   9
  10   use gridinfo
  11   use storage_module
  12   use module_debug
  13   use misc_definitions_module
  14
  15   implicit none
  16   integer :: nful
  17   integer :: interval
  18   character(len=*), dimension(nful) :: fuldates
  19   character(len=*) :: hstart
  20   integer :: ntimes
  21
  22   character(len=24) :: hdate = "0000-00-00_00:00:00.0000"
  23   character(len=24) :: hdate_output, jdate
  24   character(len=9) :: field
  25   character(len=25) :: units
  26   character(len=46) :: desc
  27   character(LEN=3)  :: out_format
  28   character(LEN=MAX_FILENAME_LEN)  :: prefix
  29   real :: xfcst
  30
  31   real :: level
  32   real, allocatable, dimension(:,:) :: scr2d, bfr2d
  33   integer :: iful, intervala, intervalb, ifv
  34   real :: awt
  35   integer :: itime
  36
  37 ! DATELEN:  length of date strings to use for our output file names.
  38   integer :: datelen
  39
  40 ! Decide the length of date strings to use for output file names.
  41 ! DATELEN is 13 for hours, 16 for minutes, and 19 for seconds.
  42   if (mod(interval,3600) == 0) then
  43      datelen = 13
  44   else if (mod(interval, 60) == 0) then
  45      datelen = 16
  46   else
  47      datelen = 19
  48   end if
  49
  50   call mprintf(.true.,STDOUT,"Subroutine DATINT: Interpolating 3-d files to fill in any missing data...")
  51   call mprintf(.true.,LOGFILE,"Subroutine DATINT: Interpolating 3-d files to fill in any missing data...")
  52
  53   TIMELOOP : do itime = 1, ntimes
  54      call geth_newdate(hdate(1:19), hstart(1:19), (itime-1)*interval)
  55      call mprintf(.true.,STDOUT,"Looking for data at time %s",s1=hdate(1:datelen))
  56      call mprintf(.true.,LOGFILE,"Looking for data at time %s",s1=hdate(1:datelen))
  57      do iful = 1, nful
  58         if (fuldates(iful).eq.hdate) then
  59            call mprintf(.true.,STDOUT,"Found file:      %s:%s", &
  60                  s1=trim(prefix),s2=hdate(1:datelen))
  61            call mprintf(.true.,LOGFILE,"Found file:      %s:%s", &
  62                  s1=trim(prefix),s2=hdate(1:datelen))
  63            cycle TIMELOOP
  64         else if ((fuldates(iful).lt.hdate) .and. &
  65              (fuldates(iful+1).gt.hdate) )then
  66
  67            call mprintf(.true.,STDOUT,"Found surrounding files:      %s: %s  and %s: %s", &
  68                s1=trim(prefix),s2=fuldates(iful)(1:datelen), &
  69                s3=trim(prefix),s4=fuldates(iful+1)(1:datelen))
  70            call mprintf(.true.,LOGFILE,"Found surrounding files:      %s: %s  and %s: %s", &
  71                s1=trim(prefix),s2=fuldates(iful)(1:datelen), &
  72                s3=trim(prefix),s4=fuldates(iful+1)(1:datelen))
  73            call mprintf(.true.,STDOUT,"Interpolating to create file:      %s: %s", &
  74                  s1=trim(prefix),s2=hdate(1:datelen))
  75            call mprintf(.true.,LOGFILE,"Interpolating to create file:      %s: %s", &
  76                  s1=trim(prefix),s2=hdate(1:datelen))
  77            call geth_idts(hdate(1:19), fuldates(iful)(1:19), intervalA)
  78            call mprintf(.true.,STDOUT,"A Time Difference = %f",f1=float(intervalA) / 3600.)
  79            call mprintf(.true.,LOGFILE,"A Time Difference = %f",f1=float(intervalA) / 3600.)
  80            call geth_idts(fuldates(iful+1)(1:19), hdate(1:19), intervalB)
  81            call mprintf(.true.,STDOUT,"B Time Difference = %f",f1=float(intervalB) / 3600.)
  82            call mprintf(.true.,LOGFILE,"B Time Difference = %f",f1=float(intervalB) / 3600.)
  83            AWT = 1. - (float(intervalA)/float(intervalA+intervalB))
  84
  85            open(10, file=trim(prefix)//':'//fuldates(iful)(1:datelen), form='unformatted', &
  86                 status='old')
  87            call clear_storage
  88            READLOOP1 : do
  89               read(10, end=44) ifv
  90               if ( ifv .eq. 5) then     ! WPS
  91                 read (10) jdate, xfcst, map%source, field, units, desc, level, &
  92                      map%nx, map%ny, map%igrid
  93                 select case (map%igrid)
  94                 case (0, 4)
  95                    read(10) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%r_earth
  96                 case (3)
  97                    read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
  98                         map%lov, map%truelat1, map%truelat2, map%r_earth
  99                 case (5)
 100                    read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 101                         map%lov, map%truelat1, map%r_earth
 102                 case default
 103                   call mprintf(.true.,ERROR, &
 104                     "Unrecognized map%%igrid: %i in DATINT 1",i1=map%igrid)
 105                 end select
 106                 read (10) map%grid_wind
 107               else if ( ifv .eq. 4 ) then          ! SI
 108                 read (10) jdate, xfcst, map%source, field, units, desc, level, &
 109                       map%nx, map%ny, map%igrid
 110                 select case (map%igrid)
 111                 case (0, 4)
 112                    read(10) map%startloc, map%lat1, map%lon1, map%dy, map%dx
 113                 case (3)
 114                    read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 115                         map%lov, map%truelat1, map%truelat2
 116                 case (5)
 117                    read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 118                         map%lov, map%truelat1
 119                 case default
 120                   call mprintf(.true.,ERROR, &
 121                     "Unrecognized map%%igrid: %i in DATINT 2",i1=map%igrid)
 122                 end select
 123               else if ( ifv .eq. 3 ) then          ! MM5
 124                 read(10) jdate, xfcst, field, units, desc, level,&
 125                      map%nx, map%ny, map%igrid
 126                 select case (map%igrid)
 127                 case (3)      ! lamcon
 128                    read (10) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
 129                         map%truelat1, map%truelat2
 130                 case (5)      ! Polar Stereographic
 131                    read (10) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
 132                         map%truelat1
 133                 case (0, 4)      ! lat/lon
 134                    read (10) map%lat1, map%lon1, map%dy, map%dx
 135                 case (1)      ! Mercator
 136                    read (10) map%lat1, map%lon1, map%dy, map%dx, map%truelat1
 137                 case default
 138                   call mprintf(.true.,ERROR, &
 139                     "Unrecognized map%%igrid: %i in DATINT 3",i1=map%igrid)
 140                 end select
 141               else
 142                 call mprintf(.true.,ERROR, &
 143                   "Unknown out_format: %i in DATINT ",i1=ifv)
 144               endif
 145               allocate(scr2d(map%nx, map%ny))
 146               read (10) scr2d
 147               call put_storage(nint(level), field, scr2d, map%nx, map%ny)
 148               deallocate(scr2d)
 149            enddo READLOOP1
 150 44         close(10)
 151
 152            open(10, file=trim(prefix)//':'//fuldates(iful+1)(1:datelen), status='old', &
 153                 form = 'unformatted')
 154            open(11, file=trim(prefix)//':'//hdate(1:datelen), status='new', form='unformatted')
 155            READLOOP2 : do
 156               read (10,END=45) ifv
 157               if ( ifv .eq. 5) then     ! WPS
 158                 read (10) jdate, xfcst, map%source, field, units, desc, level, &
 159                       map%nx, map%ny, map%igrid
 160                 select case (map%igrid)
 161                 case (0, 4)
 162                    read(10) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%r_earth
 163                 case (3)
 164                    read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 165                         map%lov, map%truelat1, map%truelat2, map%r_earth
 166                 case (5)
 167                    read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 168                       map%lov, map%truelat1, map%r_earth
 169                 case default
 170                    call mprintf(.true.,ERROR, &
 171                      "Unrecognized map%%igrid: %i in DATINT ",i1=map%igrid)
 172                 end select
 173                 read (10) map%grid_wind
 174               else if ( ifv .eq. 4 ) then          ! SI
 175                 read (10) jdate, xfcst, map%source, field, units, desc, level, &
 176                       map%nx, map%ny, map%igrid
 177                 select case (map%igrid)
 178                 case (0, 4)
 179                    read(10) map%startloc, map%lat1, map%lon1, map%dy, map%dx
 180                 case (1)
 181                    read(10) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%truelat1
 182                 case (3)
 183                    read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 184                       map%lov, map%truelat1, map%truelat2
 185                 case (5)
 186                    read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 187                       map%lov, map%truelat1
 188                 case default
 189                    call mprintf(.true.,ERROR, &
 190                      "Unrecognized map%%igrid: %i in DATINT ",i1=map%igrid)
 191                 end select
 192
 193               else if ( ifv .eq. 3 ) then          ! MM5
 194                 read(10) jdate, xfcst, field, units, desc, level,&
 195                      map%nx, map%ny, map%igrid
 196                 select case (map%igrid)
 197                 case (3)    ! lamcon
 198                    read (10) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
 199                         map%truelat1, map%truelat2
 200                 case (5)    ! Polar Stereographic
 201                    read (10) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
 202                       map%truelat1
 203                 case (0, 4)    ! lat/lon
 204                    read (10) map%lat1, map%lon1, map%dy, map%dx
 205                 case (1)    ! Mercator
 206                    read (10) map%lat1, map%lon1, map%dy, map%dx, map%truelat1
 207                 case default
 208                    call mprintf(.true.,ERROR, &
 209                      "Unrecognized map%%igrid: %i in DATINT ",i1=map%igrid)
 210                 end select
 211
 212               else
 213                 call mprintf(.true.,ERROR, &
 214                   "Unknown out_format: %i in DATINT ",i1=ifv)
 215               endif
 216
 217               allocate(scr2d(map%nx, map%ny))
 218               read (10) scr2d
 219               if (is_there(nint(level), field)) then
 220                  allocate(bfr2d(map%nx,map%ny))
 221                  call get_storage(nint(level), field, bfr2d, map%nx, map%ny)
 222                  scr2d = bfr2d * (AWT) + scr2d * (1.-AWT)
 223                  hdate_output = hdate
 224
 225                  if (out_format(1:2) .eq. 'SI') then
 226                  write(11) ifv
 227                  write(11) hdate_output, xfcst, map%source, field, units, desc, &
 228                       level, map%nx, map%ny, map%igrid
 229                  if (map%igrid == 0 .or. map%igrid == 4) then
 230                     write(11) map%startloc, map%lat1, map%lon1, map%dy, map%dx
 231                  elseif (map%igrid == 1) then
 232                     write(11) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%truelat1
 233                  elseif (map%igrid == 3) then
 234                     write (11) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 235                          map%lov, map%truelat1, map%truelat2
 236                  elseif (map%igrid == 5) then
 237                     write (11) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 238                          map%lov, map%truelat1
 239                  else
 240                    call mprintf(.true.,ERROR, &
 241                      "Unrecognized map%%igrid: %i in DATINT ",i1=map%igrid)
 242                  endif
 243
 244                  else if (out_format(1:2) .eq. 'WP') then
 245                  write(11) ifv
 246                  write(11) hdate_output, xfcst, map%source, field, units, desc, &
 247                       level, map%nx, map%ny, map%igrid
 248                  if (map%igrid == 0 .or. map%igrid == 4) then
 249                     write(11) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
 250                               map%r_earth
 251                  elseif (map%igrid == 1) then
 252                     write(11) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
 253                               map%truelat1, map%r_earth
 254                  elseif (map%igrid == 3) then
 255                     write (11) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 256                          map%lov, map%truelat1, map%truelat2, map%r_earth
 257                  elseif (map%igrid == 5) then
 258                     write (11) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 259                          map%lov, map%truelat1, map%r_earth
 260                  else
 261                    call mprintf(.true.,ERROR, &
 262                      "Unrecognized map%%igrid: %i in DATINT ",i1=map%igrid)
 263                  endif
 264                  write(11) map%grid_wind
 265
 266                  else if (out_format(1:2) .eq. 'MM') then
 267                  write (11) ifv
 268                  write (11) hdate_output, xfcst, field, units, Desc, level,&
 269                    map%nx, map%ny, map%igrid
 270                  if (map%igrid .eq. 3) then ! lamcon
 271                    write (11) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
 272                       map%truelat1, map%truelat2
 273                  elseif (map%igrid .eq. 5) then ! Polar Stereographic
 274                    write (11) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
 275                       map%truelat1
 276                  elseif (map%igrid .eq. 0 .or. map%igrid .eq. 4)then ! lat/lon
 277                    write (11) map%lat1, map%lon1, map%dy, map%dx
 278                  elseif (map%igrid.eq.1)then ! Mercator
 279                    write (11) map%lat1, map%lon1, map%dy, map%dx, map%truelat1
 280                  else
 281                    call mprintf(.true.,ERROR, &
 282                      "Unrecognized map%%igrid: %i in DATINT ",i1=map%igrid)
 283                  endif
 284                  endif
 285                  write(11) scr2d
 286               else
 287                  call mprintf(.true.,ERROR, &
 288                   "hdate = %s , fuldates = %s %s, Field = %s",s1=hdate,s2=fuldates(iful),s3=fuldates(iful+1),s4=field)
 289               endif
 290               deallocate(scr2d, bfr2d)
 291            enddo READLOOP2
 292 45         close(10)
 293            close(11)
 294            cycle TIMELOOP
 295         endif
 296      enddo
 297
 298      call mprintf(.true.,ERROR, &
 299         "Data not found: %s",s1=hdate)
 300
 301   enddo TIMELOOP
 302
 303   call mprintf(.true.,STDOUT, &
 304    "End Subroutine DATINT.")
 305   call mprintf(.true.,LOGFILE, &
 306    "End Subroutine DATINT.")
 307
 308 end subroutine datint