Merge pull request #119 from GeorgeGayno-NOAA/feature/gfsv16b

enkf_chgres_recenter_nc.fd - Remove hard-coded vertical levels.
NOAA-EMC · Aug 26, 2020 · 4641a72 · 4641a72
2 parents a1b08c9 + 25ae489
commit 4641a72
Show file tree

Hide file tree

Showing 6 changed files with 213 additions and 149 deletions.
diff --git a/modulefiles/fv3gfs/enkf_chgres_recenter_nc.wcoss_cray b/modulefiles/fv3gfs/enkf_chgres_recenter_nc.wcoss_cray
@@ -4,15 +4,22 @@
 
 module load ncep/1.0
 module load PrgEnv-intel/5.2.56
-module rm intel
-module load intel/16.3.210
-module load cray-mpich/7.2.0
+module unload intel/15.0.3.187
+module load intel/18.1.163
 module load craype-haswell
-module load cray-netcdf
+module load cray-mpich/7.2.0
 
-module load bacio-intel/2.0.2
-module load w3nco-intel/2.0.6
+# Loading netcdf modules
+module use /usrx/local/dev/modulefiles
+module load NetCDF-intel-sandybridge/4.7.4
+module load HDF5-parallel-intel-sandybridge/1.10.6
+
+module use /gpfs/hps/nco/ops/nwprod/lib/modulefiles
 module load ip-intel/3.0.0
 module load sp-intel/2.0.2
+module load w3nco-intel/2.0.6
+
+module use /usrx/local/nceplibs/modulefiles
+module load bacio-intel/2.0.2
 
 export FC=ftn
diff --git a/sorc/enkf_chgres_recenter_nc.fd/input_data.f90 b/sorc/enkf_chgres_recenter_nc.fd/input_data.f90
@@ -9,16 +9,17 @@ module input_data
  private
 
  integer, public                              :: idvc, idsl, idvm, nvcoord
- integer, public                              :: ntrac, ncldt
+ integer, public                              :: nvcoord_input, ntrac, ncldt
  integer, public                              :: ij_input, kgds_input(200)
- integer, public                              :: i_input, j_input, lev
+ integer, public                              :: i_input, j_input, lev, lev_output
  integer, public                              :: idate(6)
  integer, public                              :: icldamt, iicmr,  &
                                                  idelz,idpres,idzdt, &
                                                  irwmr,isnmr,igrle
 
 
  real, allocatable, public                    :: vcoord(:,:)
+ real, allocatable, public                    :: vcoord_input(:,:)
  real, allocatable, public                    :: clwmr_input(:,:)
  real, allocatable, public                    :: dzdt_input(:,:)
  real, allocatable, public                    :: grle_input(:,:)
@@ -52,7 +53,8 @@ subroutine read_input_data
  type(Dataset) :: indset
  type(Dimension) :: ncdim
  real, allocatable                            :: work2d(:,:),work3d(:,:,:)
- integer iret
+ integer iret, k, kk
+ real, allocatable :: ak(:), bk(:)
 
  ! hard code these values that are the same for GFS
  idvc=2
@@ -76,6 +78,19 @@ subroutine read_input_data
 
  ij_input = i_input * j_input
 
+ call read_attribute(indset, 'ak', ak)
+ call read_attribute(indset, 'bk', bk)
+
+ nvcoord_input = 2
+ allocate(vcoord_input(lev+1,nvcoord_input))
+ do k = 1, lev+1
+   kk = lev+2-k
+   vcoord_input(k,1) = ak(kk)
+   vcoord_input(k,2) = bk(kk)
+   print*,'VCOORD OF INPUT DATA ',k,vcoord_input(k,:)
+ enddo
+
+ deallocate(ak, bk)
 
  print*
  print*,"READ SURFACE PRESSURE"
@@ -297,40 +312,33 @@ subroutine read_vcoord_info
 
  implicit none
 
- integer                    :: istat, levs_vcoord, n, k
-
- print*
- print*,"OPEN VERTICAL COORD FILE: ", trim(vcoord_file)
- open(14, file=trim(vcoord_file), form='formatted', iostat=istat)
- if (istat /= 0) then
-   print*,"FATAL ERROR OPENING FILE. ISTAT IS: ", istat
-   call errexit(4)
- endif
-
- read(14, *, iostat=istat) nvcoord, levs_vcoord
- if (istat /= 0) then
-   print*,"FATAL ERROR READING FILE HEADER. ISTAT IS: ",istat
-   call errexit(5)
- endif
+ integer                    :: istat, n, k, k2
 
-!---------------------------------------------------------------------------------
-! The last value in the file is not used for the fv3 core.  Only read the first 
-! (lev + 1) values.
-!---------------------------------------------------------------------------------
+ real, allocatable :: ak(:), bk(:)
 
- allocate(vcoord(lev+1, nvcoord))
- read(14, *, iostat=istat) ((vcoord(n,k), k=1,nvcoord), n=1,lev+1)
- if (istat /= 0) then
-   print*,"FATAL ERROR READING FILE. ISTAT IS: ",istat
-   call errexit(6)
- endif
+ type(Dataset) :: refdset
 
  print*
- do k = 1, (lev+1)
-   print*,'VCOORD FOR LEV ', k, 'IS: ', vcoord(k,:)
+ print*,"READ OUTPUT VERT COORDINATE INFO FROM REFERENCE FILE: ",trim(ref_file)
+
+ refdset = open_dataset(ref_file)
+ call read_attribute(refdset, 'ak', ak)
+ call read_attribute(refdset, 'bk', bk)
+ call close_dataset(refdset)
+
+ lev_output = size(bk) - 1
+
+ nvcoord=2
+ allocate(vcoord(lev_output+1, nvcoord))
+
+ do k = 1, (lev_output+1)
+   k2 = lev_output+2 - k
+   vcoord(k,1) = ak(k2)
+   vcoord(k,2) = bk(k2)
+   print*,'VCOORD OF OUTPUT GRID ',k,vcoord(k,:)
  enddo
 
- close(14)
+ deallocate (ak, bk)
 
  end subroutine read_vcoord_info