Merge pull request ufs-community#96 from chan-hoo/feature/ic_netcdf

Add 'netcdf' and 'grib2' to input 'source' options in regional IC/LBC routines

model/fv_regional_bc.F90

@@ -248,7 +248,7 @@ module fv_regional_mod
 
  integer :: a_step, p_step, k_step, n_step
 !
- character(len=80) :: data_source
+ logical :: data_source_fv3gfs
 contains
 
 !-----------------------------------------------------------------------
@@ -1254,7 +1254,7 @@ subroutine start_regional_cold_start(Atm, dt_atmos, ak, bk, levp &
 !*** Get the source of the input data
 !-----------------------------------------------------------------------
 !
- call get_data_source(data_source,Atm%flagstruct%regional)
+ call get_data_source(data_source_fv3gfs,Atm%flagstruct%regional)
 !
  call setup_regional_BC(Atm, dt_atmos &
  ,isd, ied, jsd, jed &
@@ -1367,7 +1367,7 @@ subroutine start_regional_restart(Atm, dt_atmos &
 !*** Get the source of the input data.
 !-----------------------------------------------------------------------
 !
- call get_data_source(data_source,Atm%flagstruct%regional)
+ call get_data_source(data_source_fv3gfs,Atm%flagstruct%regional)
 !
 !-----------------------------------------------------------------------
 !*** Preliminary setup for the forecast.
@@ -1723,7 +1723,7 @@ subroutine regional_bc_data(Atm,bc_hour &
 !*** Sensible temperature
 !--------------------------
 !
- if (data_source == 'FV3GFS GAUSSIAN NEMSIO FILE') then
+ if ( data_source_fv3gfs ) then
  nlev=klev_in
  var_name_root='t'
  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
@@ -3676,7 +3676,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm &
 
 ! Compute true temperature using hydrostatic balance if not read from input.
 
- if (data_source /= 'FV3GFS GAUSSIAN NEMSIO FILE') then
+ if ( .not. data_source_fv3gfs ) then
  do k=1,npz
  BC_side%pt_BC(i,j,k) = (gz_fv(k)-gz_fv(k+1))/( rdgas*(pn1(i,k+1)-pn1(i,k))*(1.+zvir*BC_side%q_BC(i,j,k,sphum)) )
  enddo
@@ -3696,13 +3696,13 @@ subroutine remap_scalar_nggps_regional_bc(Atm &
 ! From Jan-Huey Chen's HiRAM code
 !-----------------------------------------------------------------------
 !
-! If the source is FV3GFS GAUSSIAN NEMSIO FILE then all the tracers are in the boundary files
+! If the source is FV3GFS GAUSSIAN NEMSIO/NETCDF and GRIB2 FILE then all the tracers are in the boundary files
 ! and will be read in.
 ! If the source is from old GFS or operational GSM then the tracers will be fixed in the boundaries
 ! and may not provide a very good result
 !
  if (cld_amt .gt. 0) BC_side%q_BC(:,:,:,cld_amt) = 0.
- if (trim(data_source) /= 'FV3GFS GAUSSIAN NEMSIO FILE') then
+ if ( .not. data_source_fv3gfs ) then
  if ( Atm%flagstruct%nwat .eq. 6 ) then
  do k=1,npz
  do i=is,ie
@@ -3745,7 +3745,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm &
  enddo
  enddo
  endif
- endif ! data source /= FV3GFS GAUSSIAN NEMSIO FILE
+ endif ! data source /= FV3GFS GAUSSIAN NEMSIO/NETCDF and GRIB2 FILE
 !
 ! For GFS spectral input, omega in pa/sec is stored as w in the input data so actual w(m/s) is calculated
 ! For GFS nemsio input, omega is 0, so best not to use for input since boundary data will not exist for w
@@ -3762,7 +3762,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm &
 
  call mappm(km, pe0, qp, npz, pe1, qn1, is,ie, -1, 4, Atm%ptop)
 
- if (data_source == 'FV3GFS GAUSSIAN NEMSIO FILE') then
+ if ( data_source_fv3gfs ) then
  do k=1,npz
  do i=is,ie
  BC_side%w_BC(i,j,k) = qn1(i,k)
@@ -3785,7 +3785,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm &
  enddo
  enddo
 
- else !<-- datasource /= 'FV3GFS GAUSSIAN NEMSIO FILE'
+ else !<-- datasource /= 'FV3GFS GAUSSIAN NEMSIO/NETCDF and GRIB2 FILE'
  do k=1,npz
  do i=is,ie
  BC_side%w_BC(i,j,k) = qn1(i,k)/BC_side%delp_BC(i,j,k)*BC_side%delz_BC(i,j,k)
@@ -6632,15 +6632,18 @@ end subroutine exch_uv
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 !---------------------------------------------------------------------
 
- subroutine get_data_source(source,regional)
+ subroutine get_data_source(data_source_fv3gfs,regional)
 !
 ! This routine extracts the data source information if it is present in the datafile.
 !
- character (len = 80) :: source
- integer :: ncids,sourceLength
- logical :: lstatus,regional
+ logical, intent(in):: regional
+ logical, intent(out):: data_source_fv3gfs
+
+ character (len=80) :: source 
+ logical :: lstatus
 !
 ! Use the fms call here so we can actually get the return code value.
+! The term 'source' is specified by 'chgres_cube'
 !
  if (regional) then
  lstatus = get_global_att_value('INPUT/gfs_data.nc',"source", source)
@@ -6651,6 +6654,18 @@ subroutine get_data_source(source,regional)
  if (mpp_pe() == 0) write(0,*) 'INPUT source not found ',lstatus,' set source=No Source Attribute'
  source='No Source Attribute'
  endif
+ if (mpp_pe()==0) write(*,*) 'INPUT gfs_data source string=',source
+
+! Logical flag for fv3gfs nemsio/netcdf/grib2 --------
+ if ( trim(source)=='FV3GFS GAUSSIAN NEMSIO FILE' .or. &
+ trim(source)=='FV3GFS GAUSSIAN NETCDF FILE' .or. &
+ trim(source)=='FV3GFS GRIB2 FILE' ) then
+ data_source_fv3gfs = .TRUE.
+ else
+ data_source_fv3gfs = .FALSE.
+ endif
+ if (mpp_pe()==0) write(*,*) 'data_source_fv3gfs=',data_source_fv3gfs
+
  end subroutine get_data_source
 
 !---------------------------------------------------------------------
@@ -6683,7 +6698,7 @@ subroutine set_delp_and_tracers(BC_side,npz,nwat)
  graupel = get_tracer_index(MODEL_ATMOS, 'graupel')
  cld_amt = get_tracer_index(MODEL_ATMOS, 'cld_amt')
 !
- source: if (trim(data_source) == 'FV3GFS GAUSSIAN NEMSIO FILE') then
+ source: if ( data_source_fv3gfs ) then
 !
 ! if (cld_amt > 0) BC_side%q_BC(:,:,:,cld_amt) = 0.0 ! Moorthi
  do k=1,npz
@@ -6705,7 +6720,7 @@ subroutine set_delp_and_tracers(BC_side,npz,nwat)
  enddo
  enddo
 !
- else source ! This else block is for all sources other than FV3GFS GAUSSIAN NEMSIO FILE
+ else source ! This else block is for all sources other than FV3GFS GAUSSIAN NEMSIO/NETCDF and GRIB2 FILE
 !
 ! 20160928: Adjust the mixing ratios consistently...
  do k=1,npz

tools/external_ic.F90

@@ -196,8 +196,7 @@ module external_ic_mod
  real, parameter:: zvir = rvgas/rdgas - 1.
  real(kind=R_GRID), parameter :: cnst_0p20=0.20d0
  real, parameter :: deg2rad = pi/180.
- character (len = 80),public :: source ! This tells what the input source was for the data
- character(len=27), parameter :: source_fv3gfs = 'FV3GFS GAUSSIAN NEMSIO FILE'
+ logical :: data_source_fv3gfs
  public get_external_ic, get_cubed_sphere_terrain
  public remap_scalar, remap_dwinds
 
@@ -572,9 +571,9 @@ subroutine get_nggps_ic (Atm, fv_domain, dt_atmos )
  ! *DH 20200922
 
 !
- call get_data_source(source,Atm%flagstruct%regional)
- if (trim(source) == source_fv3gfs) then
- call mpp_error(NOTE, "READING FROM REGRIDDED FV3GFS NEMSIO FILE")
+ call get_data_source(data_source_fv3gfs,Atm%flagstruct%regional)
+ if ( data_source_fv3gfs ) then
+ call mpp_error(NOTE, "READING FROM REGRIDDED FV3GFS NEMSIO/NETCDF/GRIB2 FILE")
  endif
 !
 !--- read in ak and bk from the gfs control file using fms_io read_data ---
@@ -842,7 +841,7 @@ subroutine get_nggps_ic (Atm, fv_domain, dt_atmos )
  snowwat = get_tracer_index(MODEL_ATMOS, 'snowwat')
  graupel = get_tracer_index(MODEL_ATMOS, 'graupel')
  ntclamt = get_tracer_index(MODEL_ATMOS, 'cld_amt')
- if (trim(source) == source_fv3gfs) then
+ if ( data_source_fv3gfs ) then
  do k=1,npz
  do j=js,je
  do i=is,ie
@@ -861,7 +860,7 @@ subroutine get_nggps_ic (Atm, fv_domain, dt_atmos )
  enddo
  enddo
  enddo
- else
+ else ! not NEMSIO/NETCDF/GRIB2
 !--- Add cloud condensate from GFS to total MASS
 ! 20160928: Adjust the mixing ratios consistently...
  do k=1,npz
@@ -3231,7 +3230,7 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in)
  endif
 
 !$OMP parallel do default(none) &
-!$OMP shared(sphum,liq_wat,rainwat,ice_wat,snowwat,graupel,liq_aero,ice_aero,source,  &
+!$OMP shared(sphum,liq_wat,rainwat,ice_wat,snowwat,graupel,liq_aero,ice_aero,data_source_fv3gfs, &
 !$OMP cld_amt,ncnst,npz,is,ie,js,je,km,k2,ak0,bk0,psc,t_in,zh,omga,qa,Atm,z500) &
 !$OMP private(l,m,pst,pn,gz,pe0,pn0,pe1,pn1,dp2,qp,qn1,gz_fv)
  do 5000 j=js,je
@@ -3375,7 +3374,7 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in)
 !----------------------------------------------------
 ! Compute true temperature using hydrostatic balance
 !----------------------------------------------------
- if (trim(source) /= source_fv3gfs .or. .not. present(t_in)) then
+ if ( .not.data_source_fv3gfs .or. .not. present(t_in) ) then
  do k=1,npz
 #ifdef MULTI_GASES
  Atm%pt(i,j,k) = (gz_fv(k)-gz_fv(k+1))/( rdgas*(pn1(i,k+1)-pn1(i,k))*virq(Atm%q(i,j,k,:)) )
@@ -3410,7 +3409,7 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in)
 ! seperate cloud water and cloud ice from Jan-Huey Chen's HiRAM code
 ! only use for NCEP IC and GFDL microphy
 !-----------------------------------------------------------------------
- if (trim(source) /= source_fv3gfs) then
+ if ( .not. data_source_fv3gfs ) then
  if ((Atm%flagstruct%nwat .eq. 3 .or. Atm%flagstruct%nwat .eq. 6) .and. &
  (Atm%flagstruct%ncep_ic .or. Atm%flagstruct%nggps_ic)) then
  do k=1,npz
@@ -3459,7 +3458,7 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in)
  enddo
  enddo
  endif
- endif ! data source /= FV3GFS GAUSSIAN NEMSIO FILE
+ endif ! data source /= FV3GFS GAUSSIAN NEMSIO/NETCDF and GRIB2 FILE
 
 ! For GFS spectral input, omega in pa/sec is stored as w in the input data so actual w(m/s) is calculated
 ! For GFS nemsio input, omega is 0, so best not to use for input since boundary data will not exist for w
@@ -3474,7 +3473,7 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in)
  enddo
  enddo
  call mappm(km, pe0, qp, npz, pe1, qn1, is,ie, -1, 4, Atm%ptop)
- if (trim(source) == source_fv3gfs) then
+ if ( data_source_fv3gfs ) then
  do k=1,npz
  do i=is,ie
  atm%w(i,j,k) = qn1(i,k)