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