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