feature/chgres_cube_grib2_atmos:

This commit references #6.

chgres_cube - process gfs grib2 data on its native isobaric
coordinate instead of converting it to a hybrid coordinate.
Remove routines used for this conversion.

sorc/chgres_cube.fd/grib2_util.F90

@@ -9,26 +9,18 @@ module grib2_util
 ! -------------------
 ! read_vcoord           Reads vertical coordinate data
 !
-! iso2sig               Main code for conversion of isobaric to sigma coordinates
-!
-! p2hyo, p2hyb          Routines for pressure to sigma conversion; used by iso2sig 
 !--------------------------------------------------------------------------
 
 use esmf
-use netcdf
-
-use program_setup, only   : tracers_input,num_tracers, external_model, base_install_dir
 
-use model_grid, only      : i_input,j_input, ip1_input, jp1_input
+use program_setup, only   : external_model, base_install_dir
 
-use atmdata_type
+use model_grid, only      : i_input, j_input
 
 implicit none
 
-
 contains 
 
-
 subroutine read_vcoord(isnative,vcoordi,vcoordo,lev_input,levp1_input,pt,metadata,iret)
 
   implicit none
@@ -106,231 +98,6 @@ subroutine read_vcoord(isnative,vcoordi,vcoordo,lev_input,levp1_input,pt,metadat
 
 end subroutine read_vcoord
 
-subroutine iso2sig(pi,sigma,lev_input,levp1_input,psptr,atm,clb,cub,nvars,iret)
-
-  implicit none
-
-  real(esmf_kind_r8), intent(inout)                       :: pi(lev_input), & 
-                                                           sigma(levp1_input,2)
-  integer, intent(inout)                                  :: lev_input, levp1_input, &
-                                                              clb(3),cub(3)
-  type(atmdata), intent(inout)   :: atm(:)
-  real(esmf_kind_r8), pointer, intent(inout)               :: psptr(:,:)
-  integer, intent(out)                                  :: iret
-
-
-  real                                      :: msg
-  integer                                   :: nvars
-
-
-  print*, "in iso2sig"
-  call p2hyo(pi,atm,psptr,100000.0,sigma(2:levp1_input,1),sigma(2:levp1_input,2), &
-              lev_input,clb,cub,nvars,msg,4,iret)
-  if (iret /= 0) call error_handler(" CONVERTING TO SIGMA COORDINATES. ONE OR BOTH PRESSURE &
-                                    &ARRAYS ARE NOT IN TOP TO BOTTOM ORDER.", iret)
-  iret = 0
-
-end subroutine iso2sig
-
- subroutine p2hyo(pi,xi,psfc,p0,hyao,hybo,lev_input,clb,cub,nvars,xmsg,kflag,iret) 
-      implicit none
-
-!      this routine interploates constant pres levels to hybrid
-!     the formula for the pressure of a hybrid surface is;
-!          phy(k) = pout(k) = hya(k)*p0 + hyb(k)*psfc
-!
-!     input  ["i" input ... "o" output]
-!          pi     - pressure level                     [input]
-!          psfc   - is the surface pressure pa         [input]
-!          mlon   - longitude dimension
-!          nlat   - latitude  dimension
-!          hyao   - is the "a" or pressure hybrid coef
-!          hybo   - is the "b" or sigma coeficient
-!          klevo  - number of output levels
-!          kflag  - specify how values outside the "pi" will be handled
-!                   by "outside" i mean [pout < pi(1)] or [pout > pi(klevi)]
-!                   extrapolation is via log-linear extrapolation.
-!                   =0   no extrapolation. values set to _fillvalue
-!                   =1   values set to nearest valid value
-!                   =2   values at pout less    than pi(1) set to nearest value
-!                        values at pout greater than pi(klevi) are extrapolated 
-!                   =3   values at pout less    than pi(1) are extrapolated
-!                        values at pout greater than pi(klevi) set to nearest value 
-!                   =4   values at pout less    than pi(1) are extrapolated
-!                        values at pout greater than pi(klevi) are extrapolated
-!               
-!          iret    - error code  [=0 no error detected]
-!                               [.ne.0 error detected: one or both
-!                                      pressure arrays not top->bot order]
-!     output
-!          xo     - pressure at hybrid levels [pa]
-
-      integer, intent(inout)       :: lev_input
-      integer                      :: kflag,clb(3),cub(3)
-      real, intent(in)             :: p0,xmsg
-      real, intent(inout)          :: pi(lev_input), hyao(lev_input), hybo(lev_input)
-      type(atmdata), intent(inout):: xi(:)
-      real(esmf_kind_r8), pointer, intent(inout)       :: psfc(:,:)
-
-      integer, intent(out)      :: iret
-      integer                   :: iflag, klevo, nvars
-
-      real          :: po(lev_input)
-      print*, "in p2hyo"
-      iflag = 0
-!                                                 ! ? input asending order     
-      iret   = 0
-      klevo = lev_input
-
-      !print *, hyao(klevo), ' ', hybo(klevo), ' ', psfc(1,1), ' ', p0
-      ! print *, pi
-
-      if (pi(1).gt.pi(lev_input)) then
-          iret = 1
-          return
-      end if
-
-
-      po(1)     = hyao(1)*p0     + hybo(1)*psfc(clb(1),clb(2))
-
-      po(klevo) = hyao(klevo)*p0 + hybo(klevo)*psfc(clb(1),clb(2))
-!                                                 ! ? output ascending order     
-      if (po(1).gt.po(klevo)) then
-          iret = 20 + iret
-      end if
-
-      if (iret /= 0) return
-
-      call p2hyb(pi,xi,psfc,p0,hyao,hybo,po, lev_input,clb,cub,nvars,iflag, kflag, xmsg)
-
-      return
-end subroutine p2hyo
-
-subroutine p2hyb(pi,xi,psfc,p0,hyao,hybo,po,lev_input,clb,cub,nvar,iflag, kflag, xmsg)
-  implicit none
-
-  integer, intent(in)         :: lev_input,clb(3),cub(3)
-  integer, intent(inout)      :: kflag
-  real, intent(in)          :: p0,pi(lev_input), hyao(lev_input), hybo(lev_input), xmsg
-  real, intent(inout)       :: po(lev_input)
-  type(atmdata), intent(inout):: xi(:)
-  real(esmf_kind_r8), pointer, intent(inout)       :: psfc(:,:)
-
-  integer, intent(out)      :: iflag
-
-  integer     :: nl,nlo,ml,mlo,ki,ko,nv,  dims(3),nvar
-  real        :: pimin, pimax, pomin, pomax, dxdp
-  real, allocatable :: xo(:,:,:,:)
-
-  print*, "in p2hyb"
-  dims = shape(xi(1)%var)
-  allocate(xo(dims(1),dims(2),dims(3),nvar+1))
-
-  pimin = pi(1)   
-  pimax = pi(lev_input)
-
-  do nl = clb(2),cub(2)
-    nlo = nl-clb(2)+1
-
-    do ml = clb(1),cub(1)
-      mlo = ml-clb(1)+1
-
-      do ko = clb(3),cub(3)
-        po(ko) = hyao(ko)*p0 + hybo(ko)*psfc(ml,nl)
-        xo(mlo,nlo,ko,1) = po(ko)     
-      end do
-
-      pomin = po(1)
-      pomax = po(lev_input)
-
-
-      do ko = 1,lev_input
-        xo(mlo,nlo,ko,2:nvar) = xmsg
-
-        do ki = 1,lev_input-1
-          if (po(ko) >= pimin .and. po(ko) <= pimax ) then     
-            if (po(ko).ge.pi(ki) .and. po(ko).lt.pi(ki+1)) then
-              do nv = 2,nvar
-                xo(mlo,nlo,ko,nv) = xi(nv)%var(ml,nl,ki)            &
-                               +(xi(nv)%var(ml,nl,ki+1) -xi(nv)%var(ml,nl,ki) )*    &
-                                (log(po(ko))  -log(pi(ki)))/    &
-                                (log(pi(ki+1))-log(pi(ki)))
-              enddo 
-            end if
-           else
-            if (kflag == 0) then
-              iflag = 1
-            elseif (kflag == 1) then
-              if (po(ko) < pimin) then
-                do nv = 2,nvar
-                  xo(mlo,nlo,ko,nv) = xi(nv)%var(ml,nl,1)
-                enddo
-              elseif (po(ko) > pimax) then
-                do nv = 2,nvar
-                  xo(mlo,nlo,ko,nv-1) = xi(nv)%var(ml,nl,lev_input)
-                enddo
-              end if
-            elseif (kflag == 2) then
-              if (po(ko) < pimin) then
-                do nv = 2,nvar
-                  xo(mlo,nlo,ko,nv) = xi(nv)%var(ml,nl,1)
-                enddo
-              elseif (po(ko) > pimax) then
-                do nv = 2,nvar
-                  dxdp = (xi(nv)%var(ml,nl,lev_input) -xi(nv)%var(ml,nl,lev_input-1))*   &
-                     (log(pi(lev_input))-log(pi(lev_input-1)))
-                  xo(mlo,nlo,ko,nv) = xi(nv)%var(ml,nl,lev_input)       &
-                     + (log(po(ko))-log(pi(lev_input)))*dxdp
-                enddo
-
-              end if
-            elseif (kflag == 3) then
-              if (po(ko) < pimin) then
-                do nv = 2,nvar
-                  dxdp = (xi(nv)%var(ml,nl,2) -xi(nv)%var(ml,nl,1))*      &
-                       (log(pi(2))-log(pi(1)))
-                  xo(mlo,nlo,ko,nv) = xi(nv)%var(ml,nl,1)          &
-                       + (log(po(ko))-log(pi(1)))*dxdp
-                enddo
-
-              elseif (po(ko) > pimax) then
-                do nv = 2,nvar       
-                  xo(mlo,nlo,ko,nv-1) = xi(nv)%var(ml,nl,lev_input)
-                enddo   
-              end if
-            elseif (kflag == 4) then
-              if (po(ko) < pimin) then
-                do nv = 2,nvar
-
-                  dxdp = (xi(nv)%var(ml,nl,2) -xi(nv)%var(ml,nl,1))*       &
-                         (log(pi(2))-log(pi(1)))
-                  xo(mlo,nlo,ko,nv) = xi(nv)%var(ml,nl,1)         &
-                         + (log(po(ko))-log(pi(1)))*dxdp
-                enddo
-              elseif (po(ko) > pimax) then
-                do nv = 2,nvar
-                  dxdp = (xi(nv)%var(ml,nl,lev_input) -xi(nv)%var(ml,nl,lev_input-1))*   &
-                         (log(pi(lev_input))-log(pi(lev_input-1)))
-                  xo(mlo,nlo,ko,nv) = xi(nv)%var(ml,nl,lev_input)       &
-                         + (log(po(ko))-log(pi(lev_input)))*dxdp   
-                enddo
-              end if
-            end if !kflag 
-          endif ! outside pi bounds
-         end do ! loop over klevs input
-      end do  !loop over klevs output
-    end do !loop over lon
-  end do  !loop over lat
-
-  do nv = 1,nvar
-    xi(nv)%var(clb(1):cub(1),clb(2):cub(2),clb(3):cub(3)) = xo(:,:,lev_input:1:-1,nv)
-  enddo
-
-
-  deallocate(xo)
-
- end subroutine p2hyb
-
  subroutine rh2spfh(rh_sphum,p,t)
 
   implicit none

sorc/chgres_cube.fd/input_data.F90

@@ -41,8 +41,7 @@ module input_data
                                     tracers_input, num_tracers, &
                                     input_type, num_tracers_input, &
                                     input_type, external_model, &
-                                    get_var_cond, read_from_input, tracers, &
-                                    convert_sfc
+                                    get_var_cond, read_from_input, tracers
 
  use model_grid, only             : input_grid,        &
                                     i_input, j_input,  &
@@ -2167,9 +2166,7 @@ subroutine read_input_atm_grib2_file(localpet)
 
  use wgrib2api
 
- use grib2_util, only                   : read_vcoord, iso2sig, rh2spfh, convert_omega
- use model_grid, only                   : file_is_converted
-
+ use grib2_util, only                   : read_vcoord, rh2spfh, convert_omega
 
  implicit none
 
@@ -2221,10 +2218,8 @@ subroutine read_input_atm_grib2_file(localpet)
                       "rainwat", "snowwat", "graupel", "cld_amt", "ice_nc", &
                       "rain_nc","water_nc","liq_aero","ice_aero", &
                       "sgs_tke"/)
+
  the_file = trim(data_dir_input_grid) // "/" // trim(grib2_file_input_grid)
- !if (file_is_converted) then
- !  the_file = "./test.grib2"
- !endif
 
  print*,"- READ ATMOS DATA FROM GRIB2 FILE: ", trim(the_file)
  print*,"- USE INVENTORY FILE ", inv_file
@@ -2628,7 +2623,6 @@ subroutine read_input_atm_grib2_file(localpet)
     call error_handler("IN FieldGet", rc) 
   end do
 
-  !call iso2sig(rlevs,vcoord,lev_input,levp1_input,psptr,atm,clb,cub,5+num_tracers, iret)
   do i = clb(1),cub(1)
     do j = clb(2),cub(2)
       atm(1)%var(i,j,:) = rlevs(lev_input:1:-1)
@@ -2737,9 +2731,6 @@ subroutine read_input_atm_grib2_file(localpet)
 
  endif
 
- if (localpet == 0 .and. file_is_converted .and. .not. convert_sfc) &
-      call system("rm "//trim(the_file))
-
  end subroutine read_input_atm_grib2_file
 
 !---------------------------------------------------------------------------

sorc/chgres_cube.fd/makefile

@@ -29,7 +29,7 @@ utils.o:  utils.f90
 atmdata_type.o: atmdata_type.F90
 	$(FCOMP) $(FFLAGS) -c $(ESMF_F90COMPILEPATHS) atmdata_type.F90
 
-grib2_util.o:  model_grid.o program_setup.o atmdata_type.o grib2_util.F90
+grib2_util.o:  model_grid.o program_setup.o grib2_util.F90
 	$(FCOMP) $(FFLAGS) -c $(ESMF_F90COMPILEPATHS) grib2_util.F90
 
 program_setup.o:  program_setup.f90

sorc/chgres_cube.fd/model_grid.F90

@@ -72,7 +72,6 @@ module model_grid
  integer, public                        :: ip1_target, jp1_target
  integer, public                        :: num_tiles_input_grid
  integer, public                        :: num_tiles_target_grid
- integer, public                        :: file_is_converted = 0
 
  type(esmf_grid),  public               :: input_grid
  type(esmf_grid),  public               :: target_grid