external/io_netcdf/vort.F90

   1 ! on linux, compile wrf then compile as:
   2 ! pgf90 -Mfree -I ../../main -I ../../inc -I /usr/local/netcdf-pgi/include vort.F90 libwrfio_nf.a /usr/local/netcdf-pgi/lib/libnetcdf.a ../../main/libwrflib.a
   3 ! on AIX, compile wrf then compile as:
   4 ! /lib/cpp -C -P vort.F90 > vort.f
   5 ! mpxlf -qfree=f90 -I ../../share -I ../../main -I ../../inc -I /usr/local/netcdf/include vort.f libwrfio_nf.a /usr/local/netcdf/lib/libnetcdf.a ../../main/libwrflib.a
   6
   7 module read_util_module
   8
   9 contains
  10
  11    subroutine arguments(v2file, lmore)
  12      implicit none
  13      character(len=*) :: v2file
  14      character(len=120) :: harg
  15      logical :: lmore
  16
  17      integer :: ierr, i, numarg
  18
  19      numarg = command_argument_count()
  20
  21      i = 1
  22      lmore = .false.
  23
  24      do while ( i < numarg)
  25         call get_command_argument(number=i, value=harg)
  26         print*, 'harg = ', trim(harg)
  27
  28         if (harg == "-v") then
  29            i = i + 1
  30            lmore = .true.
  31         elseif (harg == "-h") then
  32            call help
  33         endif
  34
  35      enddo
  36
  37      call get_command_argument(number=i, value=harg)
  38      v2file = harg
  39    end subroutine arguments
  40
  41    subroutine help
  42      implicit none
  43      character(len=120) :: cmd
  44      call get_command_argument(number=0, value=cmd)
  45
  46      write(*,'(/,"Usage: ", A, " [-v] v2file ")') trim(cmd)
  47      write(*,'(8x, "-v     : Print extra info")')
  48      write(*,'(8x, "v3file : MM5v3 file name to read.")')
  49      write(*,'(8x, "-h     : print this help message and exit.",/)')
  50      stop
  51    end subroutine help
  52 end module read_util_module
  53
  54
  55
  56  program readv3
  57   use wrf_data
  58   use read_util_module
  59   use module_compute_geop
  60
  61
  62   implicit none
  63 #include "wrf_status_codes.h"
  64 #include <netcdf.inc>
  65   character(len=120) :: flnm
  66   character(len=120) :: flnm2
  67   character(len=120) :: arg3
  68   character(len=19) :: DateStr
  69   character(len=19) :: DateStr2
  70   character(len=31) :: VarName
  71   character(len=31) :: VarName2
  72   integer dh1, dh2
  73
  74   integer :: flag, flag2
  75   integer :: iunit, iunit2
  76
  77   integer :: i,j,k
  78   integer :: levlim
  79   integer :: cross
  80   integer :: ndim, ndim2
  81   integer :: WrfType, WrfType2
  82   real :: time, time2
  83   real*8 :: a, b
  84   real*8 :: sum1, sum2, diff1, diff2, serr, perr, rms
  85   integer, dimension(4) :: start_index, end_index, start_index2, end_index2, end_index_u, end_index_uz
  86   integer , Dimension(3) :: MemS,MemE,PatS,PatE
  87   character (len= 4) :: staggering,   staggering2
  88   character (len= 3) :: ordering,     ordering2, ord
  89   character (len=24) :: start_date,   start_date2
  90   character (len=24) :: current_date, current_date2
  91   character (len=31) :: name,         name2,  tmpname
  92   character (len=25) :: units,        units2
  93   character (len=46) :: description,  description2
  94
  95   real, allocatable, dimension(:,:,:) :: ph, phb, p, pb
  96   real, allocatable, dimension(:,:)   :: height
  97
  98   integer ::  ids, ide, jds, jde, kds, kde,    &
  99               ims, ime, jms, jme, kms, kme,    &
 100               its, ite, jts, jte, kts, kte
 101   integer outcount
 102
 103
 104   character (len=80), dimension(3)  ::  dimnames
 105   character (len=80) :: SysDepInfo
 106
 107   integer :: l, n
 108   integer :: ikdiffs, ifdiffs
 109
 110   real, allocatable, dimension(:,:,:,:) :: data,data2
 111
 112   integer :: ierr, ierr2, ier, ier2, Status, Status_next_time, Status_next_time2, Status_next_var, Status_next_var_2
 113
 114   logical :: newtime = .TRUE.
 115   logical :: justplot, efound
 116
 117   integer, external :: iargc
 118   logical, external :: iveceq
 119
 120   levlim = -1
 121
 122   call ext_ncd_ioinit(SysDepInfo,Status)
 123   call set_wrf_debug_level ( 1 )
 124
 125
 126   Justplot = .true.
 127
 128 ! get arguments
 129 !  if ( iargc() .ge. 2 ) then
 130     call get_command_argument(number=1, value=flnm)
 131 !    call get_command_argument(number=2, value=flnm2)
 132     ierr = 0
 133     call ext_ncd_open_for_read( trim(flnm), 0, 0, "", dh1, Status)
 134     if ( Status /= 0 ) then
 135       print*,'error opening ',flnm, ' Status = ', Status ; stop
 136     endif
 137 !    call ext_ncd_open_for_read( trim(flnm2), 0, 0, "", dh2, Status)
 138 !    if ( Status /= 0 ) go to 923
 139 !    goto 924
 140 !923    continue
 141 !
 142 !! bounce here if second name is not openable -- this would mean that
 143 !! it is a field name instead.
 144 !
 145 !    print*,'could not open ',flnm2
 146 !    name = flnm2
 147 !    Justplot = .true.
 148 !924    continue
 149 !  if ( iargc() .eq. 3 ) then
 150 !    call get_command_argument(number=3, value=arg3)
 151 !    read(arg3,*)levlim
 152 !    print*,'LEVLIM = ',LEVLIM
 153 !  endif
 154 !  else
 155 !     print*,'Usage: command file1 file2'
 156 !     stop
 157 !  endif
 158
 159 !print*,'Just plot ',Justplot
 160
 161 start_index = 1
 162 end_index = 0
 163
 164 CALL ext_ncd_get_dom_ti_integer(dh1,'WEST-EAST_GRID_DIMENSION',end_index(1),1,OutCount,Status)
 165 CALL ext_ncd_get_dom_ti_integer(dh1,'BOTTOM-TOP_GRID_DIMENSION',end_index(2),1,OutCount,Status)
 166 CALL ext_ncd_get_dom_ti_integer(dh1,'SOUTH-NORTH_GRID_DIMENSION',end_index(3),1,OutCount,Status)
 167 ord = 'XZY'
 168 staggering = ' '
 169
 170
 171
 172 allocate(ph(end_index(1),end_index(2),end_index(3)))
 173 allocate(phb(end_index(1),end_index(2),end_index(3)))
 174 allocate(p(end_index(1),end_index(2),end_index(3)))
 175 allocate(pb(end_index(1),end_index(2),end_index(3)))
 176 allocate(height(end_index(1),end_index(3)))
 177
 178 ids=start_index(1); ide=end_index(1); jds=start_index(3); jde=end_index(3); kds=start_index(2); kde=end_index(2)
 179 ims=start_index(1); ime=end_index(1);   jms=start_index(3); jme=end_index(3);   kms=start_index(2); kme=end_index(2)
 180 its=start_index(1); ite=end_index(1)-1; jts=start_index(3); jte=end_index(3)-1; kts=start_index(2); kte=end_index(2)-1
 181
 182 end_index_u  = end_index - 1
 183 end_index_uz = end_index - 1
 184 end_index_uz(2) = end_index_uz(2) + 1
 185
 186
 187
 188 if ( Justplot ) then
 189   print*, 'flnm = ', trim(flnm)
 190
 191   call ext_ncd_get_next_time(dh1, DateStr, Status_next_time)
 192
 193   DO WHILE ( Status_next_time .eq. 0 )
 194     write(*,*)'Next Time ',TRIM(Datestr)
 195
 196     staggering = 'Z'
 197     name = 'PH'
 198     call ext_ncd_read_field(dh1,DateStr,TRIM(name),ph,WRF_REAL,0,0,0,ord, &
 199                             staggering, dimnames ,                      &
 200                             start_index,end_index_uz,                      & !dom
 201                             start_index,end_index,                      & !mem
 202                             start_index,end_index_uz,                      & !pat
 203                             ierr)
 204     name = 'PHB'
 205     call ext_ncd_read_field(dh1,DateStr,TRIM(name),phb,WRF_REAL,0,0,0,ord, &
 206                             staggering, dimnames ,                      &
 207                             start_index,end_index_uz,                      & !dom
 208                             start_index,end_index,                      & !mem
 209                             start_index,end_index_uz,                      & !pat
 210                             ierr)
 211     staggering = ' '
 212     name = 'P'
 213     call ext_ncd_read_field(dh1,DateStr,TRIM(name),p,WRF_REAL,0,0,0,ord, &
 214                             staggering, dimnames ,                      &
 215                             start_index,end_index_u,                      & !dom
 216                             start_index,end_index,                      & !mem
 217                             start_index,end_index_u,                      & !pat
 218                             ierr)
 219     name = 'PB'
 220     call ext_ncd_read_field(dh1,DateStr,TRIM(name),pb,WRF_REAL,0,0,0,ord, &
 221                             staggering, dimnames ,                      &
 222                             start_index,end_index_u,                      & !dom
 223                             start_index,end_index,                      & !mem
 224                             start_index,end_index_u,                      & !pat
 225                             ierr)
 226
 227     CALL compute_500mb_height  ( ph, phb, p, pb,                  &
 228                                    height,                          &
 229                                    ids, ide, jds, jde, kds, kde,    &
 230                                    ims, ime, jms, jme, kms, kme,    &
 231                                    its, ite, jts, jte, kts, kte    )
 232
 233     write(88,*)end_index_u(1),end_index_u(3),' height ',trim(Datestr)
 234     do j = 1, end_index_u(3)
 235       do i = 1, end_index_u(1)
 236         write(88,*) height(i,j)
 237       enddo
 238     enddo
 239
 240     call ext_ncd_get_next_time(dh1, DateStr, Status_next_time)
 241   enddo
 242 endif
 243
 244 end program readv3
 245
 246 ! stub for routine called by module_wrf_error (used by netcdf implementation of IO api)
 247 SUBROUTINE wrf_abort
 248   STOP
 249 END SUBROUTINE wrf_abort