WPS/ungrib/src/output.F

   1 subroutine output(hdate, nlvl, maxlvl, plvl, interval, iflag, out_format, prefix, debug_level)
   2 !                                                                             !
   3 !*****************************************************************************!
   4 !  Write output to a file.
   5 !                                                                             !
   6 !    hdate:
   7 !    nlvl:
   8 !    maxlvl:
   9 !    plvl:
  10 !    interval:
  11 !    iflag:  1 = output for ingest into rrpr ; 2 = output for hinterp
  12 !    out_format: requested output format (WPS, SI, or MM5)
  13 !                                                                             !
  14 !*****************************************************************************!
  15
  16   use table
  17   use gridinfo
  18   use storage_module
  19   use filelist
  20   use module_debug
  21   use misc_definitions_module
  22   use stringutil
  23
  24   implicit none
  25
  26   character(LEN=19) :: hdate
  27   character(LEN=24) :: hdate_output
  28   character(LEN=3)  :: out_format
  29   character(LEN=MAX_FILENAME_LEN)  :: prefix
  30   integer :: iunit = 13
  31
  32   real, pointer, dimension(:,:) :: scr2d
  33
  34   integer :: maxlvl
  35   integer nlvl, debug_level
  36   real , dimension(maxlvl) :: plvl
  37   character (LEN=9) :: field
  38   real :: level
  39   integer :: sunit = 14
  40   integer :: interval
  41   integer :: iflag
  42 ! Local Miscellaneous
  43   integer :: k, n, mm, ilev
  44   integer :: ii, jj
  45   real :: maxv, minv
  46   real :: xplv
  47   real :: xfcst = 0.
  48   character (LEN=25) :: units
  49   character (LEN=46) :: Desc
  50   logical lopen
  51
  52 ! DATELEN:  length of date strings to use for our output file names.
  53   integer :: datelen
  54
  55 ! Decide the length of date strings to use for output file names.
  56 ! DATELEN is 13 for hours, 16 for minutes, and 19 for seconds.
  57   if (mod(interval,3600) == 0) then
  58      datelen = 13
  59   elseif (mod(interval,60) == 0) then
  60      datelen = 16
  61   else
  62      datelen = 19
  63   endif
  64
  65   call get_plvls(plvl, maxlvl, nlvl)
  66
  67   if ( debug_level .ge. 0 ) then
  68   write(*,119) hdate(1:10), hdate(12:19)
  69 119 format(/,79('#'),//,'Inventory for date = ', A10,1x,A8,/)
  70
  71   write(*,advance='NO', fmt='("PRES", 2x)')
  72   WRTLOOP : do n = 1, maxvar
  73      do k = 1, n-1
  74         if (namvar(k).eq.namvar(n)) cycle WRTLOOP
  75      enddo
  76      write(*,advance='NO', fmt='(1x,A8)') namvar(n)
  77   enddo WRTLOOP
  78   write(*,advance='YES', fmt='(1x)')
  79
  80   write(*,FMT='(79("-"))')
  81   end if
  82   KLOOP : do k = 1, nlvl
  83      if ((iflag.eq.2).and.(plvl(k).gt.200100) .and. (plvl(k).lt.200200)) then
  84         cycle KLOOP
  85      endif
  86      ilev = nint(plvl(k))
  87      if ( debug_level .ge. 0 ) then
  88      write(*, advance='NO', FMT='(F6.1)') plvl(k)/100.
  89      end if
  90      MLOOP : do mm = 1, maxvar
  91         do n = 1, mm-1
  92            if (namvar(mm).eq.namvar(n)) cycle MLOOP
  93         enddo
  94         if ( debug_level .ge. 0 ) then
  95         if (is_there(ilev,namvar(mm))) then
  96            write(*, advance='NO', FMT='("  X      ")')
  97         else
  98            if ( plvl(k).gt.200000 ) then
  99              write(*, advance='NO', FMT='("  O      ")')
 100            else
 101              write(*, advance='NO', FMT='("         ")')
 102            endif
 103         endif
 104         endif
 105      enddo MLOOP
 106      if ( debug_level .ge. 0 ) then
 107      write(*,advance='YES', fmt='(1x)')
 108      endif
 109   enddo KLOOP
 110   if ( debug_level .ge. 0 ) then
 111   write(*,FMT='(79("-"))')
 112   endif
 113
 114   if (iflag.eq.1) then
 115      if (nfiles.eq.0) then
 116         open(iunit, file=trim(get_path(prefix))//'PFILE:'//HDATE(1:datelen), form='unformatted', &
 117              position='REWIND')
 118         nfiles = nfiles + 1
 119         filedates(nfiles)(1:datelen) = hdate(1:datelen)
 120      else
 121         DOFILES : do k = 1, nfiles
 122            if (hdate(1:datelen).eq.filedates(k)(1:datelen)) then
 123               open(iunit, file=trim(get_path(prefix))//'PFILE:'//HDATE(1:datelen), form='unformatted',&
 124                    position='APPEND')
 125            endif
 126         enddo DOFILES
 127         inquire (iunit, OPENED=LOPEN)
 128         if (.not. LOPEN) then
 129            open(iunit, file=trim(get_path(prefix))//'PFILE:'//HDATE(1:datelen), form='unformatted', &
 130                 position='REWIND')
 131            nfiles = nfiles + 1
 132            filedates(nfiles)(1:datelen) = hdate(1:datelen)
 133         endif
 134      endif
 135   else if (iflag.eq.2) then
 136      open(iunit, file=trim(prefix)//':'//HDATE(1:datelen), form='unformatted', &
 137           position='REWIND')
 138   endif
 139
 140 !MGD  if ( debug_level .gt. 100 ) then
 141 !MGD     write(6,*) 'begin nloop'
 142 !MGD  end if
 143   NLOOP : do n = 1, nlvl
 144
 145 !MGD  if ( debug_level .gt. 100 ) then
 146 !MGD     write(6,*) 'begin outloop'
 147 !MGD  end if
 148      OUTLOOP : do mm = 1, maxvar
 149         field = namvar(mm)
 150         do k = 1, mm-1
 151            if (field.eq.namvar(k)) cycle OUTLOOP
 152         enddo
 153         level = plvl(n)
 154         if ((iflag.eq.2).and.(level.gt.200100) .and. (level.lt.200200)) then
 155            cycle NLOOP
 156         endif
 157         ilev = nint(level)
 158         desc = ddesc(mm)
 159         if (iflag.eq.2) then
 160            if (desc.eq.' ') cycle OUTLOOP
 161         endif
 162         units = dunits(mm)
 163         if ((iflag.eq.1).or.(iflag.eq.2.and.desc(1:1).ne.' ')) then
 164           if (is_there(ilev,field)) then
 165             call get_dims(ilev, field)
 166
 167 !MGD            if ( debug_level .gt. 100 ) then
 168 !MGD               write(6,*) 'call refr_storage'
 169 !MGD            end if
 170             call refr_storage(ilev, field, scr2d, map%nx, map%ny)
 171
 172 !MGD            if ( debug_level .gt. 100 ) then
 173 !MGD               write(6,*) 'back from refr'
 174 !MGD               write(6,*) 'out_format = ',out_format
 175 !MGD            end if
 176
 177             if (out_format(1:2) .eq. 'SI') then
 178 !MGD              if ( debug_level .gt. 100 ) then
 179 !MGD                 write(6,*) 'writing in SI format'
 180 !MGD              end if
 181               write(iunit) 4
 182               hdate_output = hdate
 183               write (iunit) hdate_output, xfcst, map%source, field, units, &
 184                    Desc, level, map%nx, map%ny, map%igrid
 185               if (map%igrid.eq.3) then ! lamcon
 186                  write (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 187                       map%lov, map%truelat1, map%truelat2
 188               elseif (map%igrid.eq.5) then ! Polar Stereographic
 189                  write (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 190                       map%lov, map%truelat1
 191               elseif (map%igrid.eq.0 .or. map%igrid.eq.4)then ! lat/lon
 192                  write (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx
 193               elseif (map%igrid.eq.1)then ! Mercator
 194                  write (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
 195                       map%truelat1
 196               else
 197                  call mprintf(.true.,ERROR, &
 198                 "Unrecognized map%%igrid: %i in subroutine output 1",i1=map%igrid)
 199               endif
 200               write (iunit) scr2d
 201             else if (out_format(1:2) .eq. 'WP') then
 202                 call mprintf(.true.,DEBUG, &
 203          "writing in WPS format  iunit = %i, map%%igrid = %i",i1=iunit,i2=map%igrid)
 204               write(iunit) 5
 205               hdate_output = hdate
 206               write (iunit) hdate_output, xfcst, map%source, field, units, &
 207                    Desc, level, map%nx, map%ny, map%igrid
 208               if (map%igrid.eq.3) then ! lamcon
 209                  write (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 210                       map%lov, map%truelat1, map%truelat2, map%r_earth
 211               elseif (map%igrid.eq.5) then ! Polar Stereographic
 212                  write (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 213                       map%lov, map%truelat1, map%r_earth
 214               elseif (map%igrid.eq.0 .or. map%igrid.eq.4)then ! lat/lon
 215                  write (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
 216                       map%r_earth
 217               elseif (map%igrid.eq.1)then ! Mercator
 218                  write (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
 219                       map%truelat1, map%r_earth
 220               else
 221                  call mprintf(.true.,ERROR, &
 222                 "Unrecognized map%%igrid: %i in subroutine output 1",i1=map%igrid)
 223               endif
 224               write (iunit) map%grid_wind
 225               write (iunit) scr2d
 226             else if (out_format(1:2) .eq. 'MM') then
 227 !MGD              if ( debug_level .gt. 100 ) then
 228 !MGD             write(6,*) 'writing in MM5 format'
 229 !MGD              end if
 230               write(iunit) 3
 231               hdate_output = hdate
 232               write (iunit) hdate_output, xfcst, field, units, Desc, level,&
 233                    map%nx, map%ny, map%igrid
 234               if (map%igrid.eq.3) then ! lamcon
 235                  write (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
 236                       map%truelat1, map%truelat2
 237               elseif (map%igrid.eq.5) then ! Polar Stereographic
 238                  write (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
 239                       map%truelat1
 240               elseif (map%igrid.eq.0 .or. map%igrid.eq.4)then ! lat/lon
 241                  write (iunit) map%lat1, map%lon1, map%dy, map%dx
 242               elseif (map%igrid.eq.1)then ! Mercator
 243                  write (iunit) map%lat1, map%lon1, map%dy, map%dx, map%truelat1
 244               else
 245                  call mprintf(.true.,ERROR, &
 246                 "Unrecognized map%%igrid: %i in subroutine output 1",i1=map%igrid)
 247               endif
 248               write (iunit) scr2d
 249             endif
 250               if ( debug_level .gt. 100 ) then
 251                 call mprintf(.true.,DEBUG, &
 252                 "hdate = %s,  xfcst = %f ",s1=hdate_output,f1=xfcst)
 253                 call mprintf(.true.,DEBUG, &
 254            "map%%source = %s, field = %s, units = %s",s1=map%source,s2=field,s3=units)
 255                 call mprintf(.true.,DEBUG, &
 256                     "Desc = %s, level = %f",s1=Desc,f1=level)
 257                 call mprintf(.true.,DEBUG, &
 258                     "map%%nx = %i, map%%ny = %i",i1=map%nx,i2=map%ny)
 259               else if ( debug_level .gt. 0 ) then
 260                 call mprintf(.true.,STDOUT, &
 261               " field = %s, level = %f",s1=field,f1=level)
 262                 call mprintf(.true.,LOGFILE, &
 263               " field = %s, level = %f",s1=field,f1=level)
 264               end if
 265               if ( debug_level .gt. 100 ) then
 266               maxv = -99999.
 267               minv = 999999.
 268               do jj = 1, map%ny
 269               do ii = 1, map%nx
 270                 if (scr2d(ii,jj) .gt. maxv) maxv = scr2d(ii,jj)
 271                 if (scr2d(ii,jj) .lt. minv) minv = scr2d(ii,jj)
 272               enddo
 273               enddo
 274               call mprintf(.true.,DEBUG, &
 275                  "max value = %f , min value = %f",f1=maxv,f2=minv)
 276               end if
 277
 278               nullify(scr2d)
 279
 280            endif
 281         endif
 282      enddo OUTLOOP
 283   enddo NLOOP
 284
 285   close(iunit)
 286
 287 end subroutine output
 288