From f02d958423432ba119a83211f8b3523476e27e52 Mon Sep 17 00:00:00 2001 From: Rahul Mahajan Date: Sun, 5 Mar 2023 00:48:55 -0500 Subject: [PATCH] remove deprecated code --- src/CMakeLists.txt | 3 - src/enkf_chgres_recenter.fd/.gitignore | 3 - src/enkf_chgres_recenter.fd/CMakeLists.txt | 27 - src/enkf_chgres_recenter.fd/driver.f90 | 65 -- src/enkf_chgres_recenter.fd/input_data.f90 | 383 --------- src/enkf_chgres_recenter.fd/interp.f90 | 552 ------------ src/enkf_chgres_recenter.fd/output_data.f90 | 396 --------- src/enkf_chgres_recenter.fd/setup.f90 | 53 -- src/enkf_chgres_recenter.fd/utils.f90 | 783 ------------------ src/fv3nc2nemsio.fd/0readme | 23 - src/fv3nc2nemsio.fd/CMakeLists.txt | 15 - src/fv3nc2nemsio.fd/constants.f90 | 314 ------- src/fv3nc2nemsio.fd/fv3_main.f90 | 215 ----- src/fv3nc2nemsio.fd/fv3_module.f90 | 372 --------- src/fv3nc2nemsio.fd/kinds.f90 | 107 --- src/gaussian_sfcanl.fd/CMakeLists.txt | 3 +- src/gaussian_sfcanl.fd/gaussian_sfcanl.f90 | 571 +------------ src/regrid_nemsio.fd/CMakeLists.txt | 25 - src/regrid_nemsio.fd/constants.f90 | 314 ------- src/regrid_nemsio.fd/fv3_interface.f90 | 779 ----------------- src/regrid_nemsio.fd/gfs_nems_interface.f90 | 595 ------------- .../interpolation_interface.f90 | 335 -------- src/regrid_nemsio.fd/kinds.f90 | 107 --- src/regrid_nemsio.fd/main.f90 | 92 -- src/regrid_nemsio.fd/mpi_interface.f90 | 89 -- src/regrid_nemsio.fd/namelist_def.f90 | 181 ---- src/regrid_nemsio.fd/netcdfio_interface.f90 | 592 ------------- src/regrid_nemsio.fd/physcons.f90 | 77 -- .../regrid_nemsio_interface.f90 | 50 -- src/regrid_nemsio.fd/variable_interface.f90 | 66 -- 30 files changed, 22 insertions(+), 7165 deletions(-) delete mode 100644 src/enkf_chgres_recenter.fd/.gitignore delete mode 100644 src/enkf_chgres_recenter.fd/CMakeLists.txt delete mode 100644 src/enkf_chgres_recenter.fd/driver.f90 delete mode 100644 src/enkf_chgres_recenter.fd/input_data.f90 delete mode 100644 src/enkf_chgres_recenter.fd/interp.f90 delete mode 100644 src/enkf_chgres_recenter.fd/output_data.f90 delete mode 100644 src/enkf_chgres_recenter.fd/setup.f90 delete mode 100644 src/enkf_chgres_recenter.fd/utils.f90 delete mode 100644 src/fv3nc2nemsio.fd/0readme delete mode 100644 src/fv3nc2nemsio.fd/CMakeLists.txt delete mode 100644 src/fv3nc2nemsio.fd/constants.f90 delete mode 100644 src/fv3nc2nemsio.fd/fv3_main.f90 delete mode 100644 src/fv3nc2nemsio.fd/fv3_module.f90 delete mode 100644 src/fv3nc2nemsio.fd/kinds.f90 delete mode 100644 src/regrid_nemsio.fd/CMakeLists.txt delete mode 100644 src/regrid_nemsio.fd/constants.f90 delete mode 100644 src/regrid_nemsio.fd/fv3_interface.f90 delete mode 100644 src/regrid_nemsio.fd/gfs_nems_interface.f90 delete mode 100644 src/regrid_nemsio.fd/interpolation_interface.f90 delete mode 100644 src/regrid_nemsio.fd/kinds.f90 delete mode 100644 src/regrid_nemsio.fd/main.f90 delete mode 100644 src/regrid_nemsio.fd/mpi_interface.f90 delete mode 100644 src/regrid_nemsio.fd/namelist_def.f90 delete mode 100644 src/regrid_nemsio.fd/netcdfio_interface.f90 delete mode 100644 src/regrid_nemsio.fd/physcons.f90 delete mode 100644 src/regrid_nemsio.fd/regrid_nemsio_interface.f90 delete mode 100644 src/regrid_nemsio.fd/variable_interface.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f0719c52..3c024ea0 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,7 +1,4 @@ -add_subdirectory(enkf_chgres_recenter.fd) add_subdirectory(enkf_chgres_recenter_nc.fd) -add_subdirectory(fv3nc2nemsio.fd) -add_subdirectory(regrid_nemsio.fd) add_subdirectory(gaussian_sfcanl.fd) add_subdirectory(gfs_bufr.fd) add_subdirectory(reg2grb2.fd) diff --git a/src/enkf_chgres_recenter.fd/.gitignore b/src/enkf_chgres_recenter.fd/.gitignore deleted file mode 100644 index 544aec4c..00000000 --- a/src/enkf_chgres_recenter.fd/.gitignore +++ /dev/null @@ -1,3 +0,0 @@ -*.exe -*.o -*.mod diff --git a/src/enkf_chgres_recenter.fd/CMakeLists.txt b/src/enkf_chgres_recenter.fd/CMakeLists.txt deleted file mode 100644 index 3bbdbbf0..00000000 --- a/src/enkf_chgres_recenter.fd/CMakeLists.txt +++ /dev/null @@ -1,27 +0,0 @@ -list(APPEND fortran_src -driver.f90 -input_data.f90 -interp.f90 -output_data.f90 -setup.f90 -utils.f90 -) - -if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$") - set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -r8 -i4 -fp-model precise") -elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$") - set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fdefault-real-8") -endif() - -set(exe_name enkf_chgres_recenter.x) -add_executable(${exe_name} ${fortran_src}) -target_link_libraries(${exe_name} PRIVATE nemsio::nemsio - bacio::bacio_4 - ip::ip_d - sp::sp_d - w3emc::w3emc_d) -if(OpenMP_Fortran_FOUND) - target_link_libraries(${exe_name} PRIVATE OpenMP::OpenMP_Fortran) -endif() - -install(TARGETS ${exe_name} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/src/enkf_chgres_recenter.fd/driver.f90 b/src/enkf_chgres_recenter.fd/driver.f90 deleted file mode 100644 index 02a138ae..00000000 --- a/src/enkf_chgres_recenter.fd/driver.f90 +++ /dev/null @@ -1,65 +0,0 @@ - program recenter - - use setup, only : program_setup - use interp, only : gaus_to_gaus, adjust_for_terrain - use input_data, only : read_input_data, & - read_vcoord_info - use output_data, only : set_output_grid, write_output_data - - implicit none - - call w3tagb('CHGRES_RECENTER',2018,0179,0055,'NP20') - - print*,"STARTING PROGRAM" - -!-------------------------------------------------------- -! Read configuration namelist. -!-------------------------------------------------------- - - call program_setup - -!-------------------------------------------------------- -! Read input grid data -!-------------------------------------------------------- - - call read_input_data - -!-------------------------------------------------------- -! Read vertical coordinate info -!-------------------------------------------------------- - - call read_vcoord_info - -!-------------------------------------------------------- -! Get output grid specs -!-------------------------------------------------------- - - call set_output_grid - -!-------------------------------------------------------- -! Interpolate data to output grid -!-------------------------------------------------------- - - call gaus_to_gaus - -!-------------------------------------------------------- -! Adjust output fields for differences between -! interpolated and external terrain. -!-------------------------------------------------------- - - call adjust_for_terrain - -!-------------------------------------------------------- -! Write output data to file. -!-------------------------------------------------------- - - call write_output_data - - print* - print*,"PROGRAM FINISHED NORMALLY!" - - call w3tage('CHGRES_RECENTER') - - stop - - end program recenter diff --git a/src/enkf_chgres_recenter.fd/input_data.f90 b/src/enkf_chgres_recenter.fd/input_data.f90 deleted file mode 100644 index 704aa58c..00000000 --- a/src/enkf_chgres_recenter.fd/input_data.f90 +++ /dev/null @@ -1,383 +0,0 @@ - module input_data - - use nemsio_module - use utils - use setup - - implicit none - - private - - integer, public :: idvc, idsl, idvm, nvcoord - integer, public :: ntrac, ncldt,icldamt - integer, public :: ij_input, kgds_input(200) - integer(nemsio_intkind), public :: i_input, j_input, lev - integer(nemsio_intkind), public :: idate(7) - - logical, public :: gfdl_mp - - real, allocatable, public :: vcoord(:,:) - real, allocatable, public :: clwmr_input(:,:) - real, allocatable, public :: dzdt_input(:,:) - real, allocatable, public :: grle_input(:,:) - real, allocatable, public :: cldamt_input(:,:) - real, allocatable, public :: hgt_input(:) - real, allocatable, public :: icmr_input(:,:) - real, allocatable, public :: o3mr_input(:,:) - real, allocatable, public :: rwmr_input(:,:) - real, allocatable, public :: sfcp_input(:) - real, allocatable, public :: snmr_input(:,:) - real, allocatable, public :: spfh_input(:,:) - real, allocatable, public :: tmp_input(:,:) - real, allocatable, public :: ugrd_input(:,:) - real, allocatable, public :: vgrd_input(:,:) - - public :: read_input_data - public :: read_vcoord_info - - contains - - subroutine read_input_data - -!------------------------------------------------------------------------------------- -! Read input grid data from a nemsio file. -!------------------------------------------------------------------------------------- - - implicit none - - character(len=20) :: vlevtyp, vname - character(len=50), allocatable :: recname(:) - - integer(nemsio_intkind) :: vlev, iret, idum, nrec - integer :: n - - real(nemsio_realkind), allocatable :: dummy(:) - - type(nemsio_gfile) :: gfile - - call nemsio_init(iret) - - print* - print*,"OPEN INPUT FILE: ",trim(input_file) - call nemsio_open(gfile, input_file, "read", iret=iret) - if (iret /= 0) then - print*,"FATAL ERROR OPENING FILE: ",trim(input_file) - print*,"IRET IS: ", iret - call errexit(2) - endif - - print*,"GET INPUT FILE HEADER" - call nemsio_getfilehead(gfile, iret=iret, nrec=nrec, idate=idate, & - dimx=i_input, dimy=j_input, dimz=lev) - if (iret /= 0) goto 67 - - print*,'DIMENSIONS OF DATA ARE: ', i_input, j_input, lev - print*,'DATE OF DATA IS: ', idate - - ij_input = i_input * j_input - - allocate(recname(nrec)) - - call nemsio_getfilehead(gfile, iret=iret, recname=recname) - if (iret /= 0) goto 67 - - gfdl_mp = .false. ! Zhao-Carr MP - do n = 1, nrec - if (trim(recname(n)) == "icmr") then - gfdl_mp = .true. ! GFDL MP - exit - endif - enddo - - icldamt = 0 - do n = 1, nrec - if (trim(recname(n)) == "cld_amt") then - icldamt = 1 ! 3D cloud amount present - exit - endif - enddo - - call nemsio_getfilehead(gfile, iret=iret, idvc=idum) - if (iret /= 0) goto 67 - idvc = idum - print*,'IDVC IS: ', idvc - - call nemsio_getfilehead(gfile, iret=iret, idsl=idum) - if (iret /= 0) goto 67 - idsl = idum - print*,'IDSL IS: ', idsl - - call nemsio_getfilehead(gfile, iret=iret, idvm=idum) - if (iret /= 0) goto 67 - idvm = idum - print*,'IDVM IS: ', idvm - - if (gfdl_mp) then - ntrac = 7 + icldamt - ncldt = 5 - else - ntrac = 3 - ncldt = 1 - endif - - allocate(dummy(ij_input)) - - ! figure out the sign of delz - print*,"READ DELZ FOR SIGN CHECK" - vlev = 1 - vlevtyp = "mid layer" - vname = "delz" - call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret) - if (iret /= 0) goto 67 - if ( sum(dummy) > 0 ) then - flipdelz = .false. - print*,"DELZ IS POSITIVE" - else - flipdelz = .true. - print*,"DELZ IS NEGATIVE" - end if - - print* - print*,"READ SURFACE PRESSURE" - vlev = 1 - vlevtyp = "sfc" - vname = "pres" - call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret) - if (iret /= 0) goto 67 - - allocate(sfcp_input(ij_input)) - sfcp_input = dummy - print*,'MAX/MIN SURFACE PRESSURE: ',maxval(sfcp_input), minval(sfcp_input) - - print* - print*,"READ SURFACE HEIGHT" - vlev = 1 - vlevtyp = "sfc" - vname = "hgt" - call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret) - if (iret /= 0) goto 67 - - allocate(hgt_input(ij_input)) - hgt_input = dummy - print*,'MAX/MIN SURFACE HEIGHT: ',maxval(hgt_input), minval(hgt_input) - - print* - print*,"READ U WIND" - vname = "ugrd" - vlevtyp = "mid layer" - allocate(ugrd_input(ij_input,lev)) - do vlev = 1, lev - call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret) - if (iret /= 0) goto 67 - ugrd_input(:,vlev) = dummy - print*,'MAX/MIN U WIND AT LEVEL ',vlev, "IS: ", maxval(ugrd_input(:,vlev)), minval(ugrd_input(:,vlev)) - enddo - - print* - print*,"READ V WIND" - vname = "vgrd" - vlevtyp = "mid layer" - allocate(vgrd_input(ij_input,lev)) - do vlev = 1, lev - call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret) - if (iret /= 0) goto 67 - vgrd_input(:,vlev) = dummy - print*,'MAX/MIN V WIND AT LEVEL ', vlev, "IS: ", maxval(vgrd_input(:,vlev)), minval(vgrd_input(:,vlev)) - enddo - - print* - print*,"READ TEMPERATURE" - vname = "tmp" - vlevtyp = "mid layer" - allocate(tmp_input(ij_input,lev)) - do vlev = 1, lev - call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret) - if (iret /= 0) goto 67 - tmp_input(:,vlev) = dummy(:) - print*,'MAX/MIN TEMPERATURE AT LEVEL ', vlev, 'IS: ', maxval(tmp_input(:,vlev)), minval(tmp_input(:,vlev)) - enddo - - print* - print*,"READ SPECIFIC HUMIDITY" - vname = "spfh" - vlevtyp = "mid layer" - allocate(spfh_input(ij_input,lev)) - do vlev = 1, lev - call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret) - if (iret /= 0) goto 67 - spfh_input(:,vlev) = dummy - print*,'MAX/MIN SPECIFIC HUMIDITY AT LEVEL ', vlev, 'IS: ', maxval(spfh_input(:,vlev)), minval(spfh_input(:,vlev)) - enddo - - print* - print*,"READ CLOUD LIQUID WATER" - vname = "clwmr" - vlevtyp = "mid layer" - allocate(clwmr_input(ij_input,lev)) - do vlev = 1, lev - call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret) - if (iret /= 0) goto 67 - clwmr_input(:,vlev) = dummy - print*,'MAX/MIN CLOUD LIQUID WATER AT LEVEL ', vlev, 'IS: ', maxval(clwmr_input(:,vlev)), minval(clwmr_input(:,vlev)) - enddo - - print* - print*,"READ OZONE" - vname = "o3mr" - vlevtyp = "mid layer" - allocate(o3mr_input(ij_input,lev)) - do vlev = 1, lev - call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret) - if (iret /= 0) goto 67 - o3mr_input(:,vlev) = dummy - print*,'MAX/MIN OZONE AT LEVEL ', vlev, 'IS: ', maxval(o3mr_input(:,vlev)), minval(o3mr_input(:,vlev)) - enddo - - print* - print*,"READ DZDT" - vname = "dzdt" - vlevtyp = "mid layer" - allocate(dzdt_input(ij_input,lev)) - do vlev = 1, lev - call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret) - if (iret /= 0) goto 67 - dzdt_input(:,vlev) = dummy - print*,'MAX/MIN DZDT AT LEVEL ', vlev, 'IS: ', maxval(dzdt_input(:,vlev)), minval(dzdt_input(:,vlev)) - enddo - - if (gfdl_mp) then - - print* - print*,"READ RWMR" - vname = "rwmr" - vlevtyp = "mid layer" - allocate(rwmr_input(ij_input,lev)) - do vlev = 1, lev - call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret) - if (iret /= 0) goto 67 - rwmr_input(:,vlev) = dummy - print*,'MAX/MIN RWMR AT LEVEL ', vlev, 'IS: ', maxval(rwmr_input(:,vlev)), minval(rwmr_input(:,vlev)) - enddo - - print* - print*,"READ ICMR" - vname = "icmr" - vlevtyp = "mid layer" - allocate(icmr_input(ij_input,lev)) - do vlev = 1, lev - call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret) - if (iret /= 0) goto 67 - icmr_input(:,vlev) = dummy - print*,'MAX/MIN ICMR AT LEVEL ', vlev, 'IS: ', maxval(icmr_input(:,vlev)), minval(icmr_input(:,vlev)) - enddo - - print* - print*,"READ SNMR" - vname = "snmr" - vlevtyp = "mid layer" - allocate(snmr_input(ij_input,lev)) - do vlev = 1, lev - call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret) - if (iret /= 0) goto 67 - snmr_input(:,vlev) = dummy - print*,'MAX/MIN SNMR AT LEVEL ', vlev, 'IS: ', maxval(snmr_input(:,vlev)), minval(snmr_input(:,vlev)) - enddo - - print* - print*,"READ GRLE" - vname = "grle" - vlevtyp = "mid layer" - allocate(grle_input(ij_input,lev)) - do vlev = 1, lev - call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret) - if (iret /= 0) goto 67 - grle_input(:,vlev) = dummy - print*,'MAX/MIN GRLE AT LEVEL ', vlev, 'IS: ', maxval(grle_input(:,vlev)), minval(grle_input(:,vlev)) - enddo - - if (icldamt == 1) then - print* - print*,"READ CLD_AMT" - vname = "cld_amt" - vlevtyp = "mid layer" - allocate(cldamt_input(ij_input,lev)) - do vlev = 1, lev - write(6,*) 'read ',vname,' on ',vlev - call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret) - if (iret /= 0) goto 67 - cldamt_input(:,vlev) = dummy - print*,'MAX/MIN CLD_AMT AT LEVEL ', vlev, 'IS: ', maxval(cldamt_input(:,vlev)), minval(cldamt_input(:,vlev)) - enddo - endif - - endif - - deallocate(dummy) - - print*,"CLOSE FILE" - call nemsio_close(gfile, iret=iret) - - call nemsio_finalize() - -!--------------------------------------------------------------------------------------- -! Set the grib 1 grid description array need by the NCEP IPOLATES library. -!--------------------------------------------------------------------------------------- - - call calc_kgds(i_input, j_input, kgds_input) - - return - - 67 continue - - print*,"FATAL ERROR READING FILE: ", trim(input_file) - print*,"IRET IS: ", iret - call errexit(3) - - end subroutine read_input_data - - subroutine read_vcoord_info - -!--------------------------------------------------------------------------------- -! Read vertical coordinate information. -!--------------------------------------------------------------------------------- - - 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 - -!--------------------------------------------------------------------------------- -! The last value in the file is not used for the fv3 core. Only read the first -! (lev + 1) values. -!--------------------------------------------------------------------------------- - - 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 - - print* - do k = 1, (lev+1) - print*,'VCOORD FOR LEV ', k, 'IS: ', vcoord(k,:) - enddo - - close(14) - - end subroutine read_vcoord_info - - end module input_data diff --git a/src/enkf_chgres_recenter.fd/interp.f90 b/src/enkf_chgres_recenter.fd/interp.f90 deleted file mode 100644 index bb2afedb..00000000 --- a/src/enkf_chgres_recenter.fd/interp.f90 +++ /dev/null @@ -1,552 +0,0 @@ - module interp - - use nemsio_module - - implicit none - - private - - real, allocatable :: sfcp_b4_adj_output(:) - real, allocatable :: clwmr_b4_adj_output(:,:) - real, allocatable :: dzdt_b4_adj_output(:,:) - real, allocatable :: grle_b4_adj_output(:,:) - real, allocatable :: cldamt_b4_adj_output(:,:) - real, allocatable :: icmr_b4_adj_output(:,:) - real, allocatable :: o3mr_b4_adj_output(:,:) - real, allocatable :: rwmr_b4_adj_output(:,:) - real, allocatable :: snmr_b4_adj_output(:,:) - real, allocatable :: spfh_b4_adj_output(:,:) - real, allocatable :: tmp_b4_adj_output(:,:) - real, allocatable :: ugrd_b4_adj_output(:,:) - real, allocatable :: vgrd_b4_adj_output(:,:) - - public :: adjust_for_terrain - public :: gaus_to_gaus - - contains - - subroutine adjust_for_terrain - -!--------------------------------------------------------------------------------- -! Adjust fields based on differences between the interpolated and external -! terrain. -!--------------------------------------------------------------------------------- - - use input_data - use output_data - use utils - use setup - - implicit none - - integer :: k - - real, allocatable :: pres_b4_adj_output(:,:) - real, allocatable :: pres_output(:,:) - real, allocatable :: q_b4_adj_output(:,:,:), q_output(:,:,:) - -!--------------------------------------------------------------------------------- -! First, compute the mid-layer pressure using the interpolated surface pressure. -!--------------------------------------------------------------------------------- - - allocate(pres_b4_adj_output(ij_output,lev)) - pres_b4_adj_output = 0.0 - - print* - print*,"COMPUTE MID-LAYER PRESSURE FROM INTERPOLATED SURFACE PRESSURE." - call newpr1(ij_output, lev, idvc, idsl, nvcoord, vcoord, & - sfcp_b4_adj_output, pres_b4_adj_output) - -!print*,'after newpr1, pres b4 adj: ', pres_b4_adj_output(ij_output/2,:) - -!--------------------------------------------------------------------------------- -! Adjust surface pressure based on differences between interpolated and -! grid terrain. -!--------------------------------------------------------------------------------- - - allocate(sfcp_output(ij_output)) - sfcp_output = 0.0 - - print*,"ADJUST SURFACE PRESSURE BASED ON TERRAIN DIFFERENCES" - call newps(hgt_output, sfcp_b4_adj_output, ij_output, & - lev, pres_b4_adj_output, tmp_b4_adj_output, & - spfh_b4_adj_output, hgt_external_output, sfcp_output) - -!print*,'after newps ',sfcp_b4_adj_output(ij_output/2),sfcp_output(ij_output/2) - - deallocate(sfcp_b4_adj_output) - -!--------------------------------------------------------------------------------- -! Recompute mid-layer pressure based on the adjusted surface pressure. -!--------------------------------------------------------------------------------- - - allocate(pres_output(ij_output, lev)) - pres_output = 0.0 - - allocate(dpres_output(ij_output, lev)) - dpres_output = 0.0 - - print*,"RECOMPUTE MID-LAYER PRESSURE." - call newpr1(ij_output, lev, idvc, idsl, nvcoord, vcoord, & - sfcp_output, pres_output, dpres_output) - -!do k = 1, lev -! print*,'after newpr1 ',pres_b4_adj_output(ij_output/2,k),pres_output(ij_output/2,k), dpres_output(ij_output/2,k) -!enddo - -!--------------------------------------------------------------------------------- -! Vertically interpolate from the pre-adjusted to the adjusted mid-layer -! pressures. -!--------------------------------------------------------------------------------- - - allocate(q_b4_adj_output(ij_output,lev,ntrac)) - q_b4_adj_output(:,:,1) = spfh_b4_adj_output(:,:) - q_b4_adj_output(:,:,2) = o3mr_b4_adj_output(:,:) - q_b4_adj_output(:,:,3) = clwmr_b4_adj_output(:,:) - if (gfdl_mp) then - q_b4_adj_output(:,:,4) = rwmr_b4_adj_output(:,:) - q_b4_adj_output(:,:,5) = icmr_b4_adj_output(:,:) - q_b4_adj_output(:,:,6) = snmr_b4_adj_output(:,:) - q_b4_adj_output(:,:,7) = grle_b4_adj_output(:,:) - if (icldamt == 1) q_b4_adj_output(:,:,8) = cldamt_b4_adj_output(:,:) - endif - - allocate(q_output(ij_output,lev,ntrac)) - q_output = 0.0 - - allocate(dzdt_output(ij_output,lev)) - dzdt_output = 0.0 - - allocate(ugrd_output(ij_output,lev)) - ugrd_output=0.0 - - allocate(vgrd_output(ij_output,lev)) - vgrd_output=0.0 - - allocate(tmp_output(ij_output,lev)) - tmp_output=0.0 - - print*,"VERTICALLY INTERPOLATE TO NEW PRESSURE LEVELS" - call vintg(ij_output, lev, lev, ntrac, pres_b4_adj_output, & - ugrd_b4_adj_output, vgrd_b4_adj_output, tmp_b4_adj_output, q_b4_adj_output, & - dzdt_b4_adj_output, pres_output, ugrd_output, vgrd_output, tmp_output, & - q_output, dzdt_output) - - deallocate (dzdt_b4_adj_output, q_b4_adj_output) - deallocate (pres_b4_adj_output, pres_output) - - allocate(spfh_output(ij_output,lev)) - spfh_output = q_output(:,:,1) - allocate(o3mr_output(ij_output,lev)) - o3mr_output = q_output(:,:,2) - allocate(clwmr_output(ij_output,lev)) - clwmr_output = q_output(:,:,3) - if (gfdl_mp) then - allocate(rwmr_output(ij_output,lev)) - rwmr_output = q_output(:,:,4) - allocate(icmr_output(ij_output,lev)) - icmr_output = q_output(:,:,5) - allocate(snmr_output(ij_output,lev)) - snmr_output = q_output(:,:,6) - allocate(grle_output(ij_output,lev)) - grle_output = q_output(:,:,7) - if (icldamt == 1) then - allocate(cldamt_output(ij_output,lev)) - cldamt_output = q_output(:,:,8) - endif - endif - - deallocate(q_output) - -!do k = 1, lev -!print*,'after vintg tmp ',tmp_b4_adj_output(ij_output/2,k),tmp_output(ij_output/2,k) -!enddo - - deallocate(tmp_b4_adj_output) - -!do k = 1, lev -!print*,'after vintg u ',ugrd_b4_adj_output(ij_output/2,k),ugrd_output(ij_output/2,k) -!enddo - - deallocate(ugrd_b4_adj_output) - -!do k = 1, lev -!print*,'after vintg v ',vgrd_b4_adj_output(ij_output/2,k),vgrd_output(ij_output/2,k) -!enddo - - deallocate(vgrd_b4_adj_output) - -!do k = 1, lev -!print*,'after vintg spfh ',spfh_b4_adj_output(ij_output/2,k),spfh_output(ij_output/2,k) -!enddo - - deallocate(spfh_b4_adj_output) - -!do k = 1, lev -!print*,'after vintg o3 ',o3mr_b4_adj_output(ij_output/2,k),o3mr_output(ij_output/2,k) -!enddo - - deallocate(o3mr_b4_adj_output) - -!do k = 1, lev -!print*,'after vintg clw ',clwmr_b4_adj_output(ij_output/2,k),clwmr_output(ij_output/2,k) -!enddo - - deallocate(clwmr_b4_adj_output) - - if (gfdl_mp) then - -! do k = 1, lev -! print*,'after vintg rw ',rwmr_b4_adj_output(ij_output/2,k),rwmr_output(ij_output/2,k) -! enddo - - deallocate(rwmr_b4_adj_output) - -! do k = 1, lev -! print*,'after vintg ic ',icmr_b4_adj_output(ij_output/2,k),icmr_output(ij_output/2,k) -! enddo - - deallocate(icmr_b4_adj_output) - -! do k = 1, lev -! print*,'after vintg sn ',snmr_b4_adj_output(ij_output/2,k),snmr_output(ij_output/2,k) -! enddo - - deallocate(snmr_b4_adj_output) - -! do k = 1, lev -! print*,'after vintg grle ',grle_b4_adj_output(ij_output/2,k),grle_output(ij_output/2,k) -! enddo - - deallocate(grle_b4_adj_output) - - if (icldamt == 1) then -! do k = 1, lev -! print*,'after vintg cld_amt ',cldamt_b4_adj_output(ij_output/2,k),cldamt_output(ij_output/2,k) -! enddo - - deallocate(cldamt_b4_adj_output) - endif - - - endif - - allocate(delz_output(ij_output, lev)) - delz_output = 0.0 - - call compute_delz(ij_output, lev, vcoord(:,1), vcoord(:,2), sfcp_output, hgt_output, & - tmp_output, spfh_output, delz_output, flipdelz) - - deallocate(hgt_output) - - end subroutine adjust_for_terrain - - subroutine gaus_to_gaus - -!---------------------------------------------------------------------------------- -! Interpolate data from the input to output grid using IPOLATES library. -!---------------------------------------------------------------------------------- - - use output_data - use input_data - use setup - - implicit none - - integer :: ip, ipopt(20) - integer :: num_fields - integer :: iret, numpts - integer, allocatable :: ibi(:), ibo(:) - - logical*1, allocatable :: bitmap_input(:,:), bitmap_output(:,:) - - real, allocatable :: data_input(:,:) - real, allocatable :: data_output(:,:), crot(:), srot(:) - - print* - print*,'INTERPOLATE DATA TO OUTPUT GRID' - - ip = 0 ! bilinear - ipopt = 0 - -!---------------------------------------------------------------------------------- -! Do 2-D fields first -!---------------------------------------------------------------------------------- - - num_fields = 1 - - allocate(ibi(num_fields)) - ibi = 0 ! no bitmap - allocate(ibo(num_fields)) - ibo = 0 ! no bitmap - - allocate(bitmap_input(ij_input,num_fields)) - bitmap_input = .true. - allocate(bitmap_output(ij_output,num_fields)) - bitmap_output = .true. - - allocate(rlat_output(ij_output)) - rlat_output = 0.0 - allocate(rlon_output(ij_output)) - rlon_output = 0.0 - -!---------------- -! Surface height -!---------------- - - allocate(data_input(ij_input,num_fields)) - data_input(:,num_fields) = hgt_input(:) - deallocate(hgt_input) - - allocate(data_output(ij_output,num_fields)) - data_output = 0 - - print*,"INTERPOLATE SURFACE HEIGHT" - call ipolates(ip, ipopt, kgds_input, kgds_output, ij_input, ij_output,& - num_fields, ibi, bitmap_input, data_input, & - numpts, rlat_output, rlon_output, ibo, bitmap_output, & - data_output, iret) - if (iret /= 0) goto 89 - - allocate(hgt_output(ij_output)) - hgt_output = data_output(:,num_fields) - -!------------------ -! surface pressure -!------------------ - - data_input(:,num_fields) = sfcp_input(:) - deallocate(sfcp_input) - - print*,"INTERPOLATE SURFACE PRESSURE" - call ipolates(ip, ipopt, kgds_input, kgds_output, ij_input, ij_output,& - num_fields, ibi, bitmap_input, data_input, & - numpts, rlat_output, rlon_output, ibo, bitmap_output, & - data_output, iret) - if (iret /= 0) goto 89 - - allocate(sfcp_b4_adj_output(ij_output)) - sfcp_b4_adj_output = data_output(:,num_fields) - - deallocate(ibi, ibo, bitmap_input, bitmap_output, data_input, data_output) - -!---------------------------------------------------------------------------------- -! 3d scalars -!---------------------------------------------------------------------------------- - - num_fields = lev - - allocate(ibi(num_fields)) - ibi = 0 ! no bitmap - allocate(ibo(num_fields)) - ibo = 0 ! no bitmap - - allocate(bitmap_input(ij_input,num_fields)) - bitmap_input = .true. - allocate(bitmap_output(ij_output,num_fields)) - bitmap_output = .true. - -!------------- -! Temperature -!------------- - - allocate(tmp_b4_adj_output(ij_output,num_fields)) - tmp_b4_adj_output = 0 - - print*,'INTERPOLATE TEMPERATURE' - call ipolates(ip, ipopt, kgds_input, kgds_output, ij_input, ij_output,& - num_fields, ibi, bitmap_input, tmp_input, & - numpts, rlat_output, rlon_output, ibo, bitmap_output, & - tmp_b4_adj_output, iret) - if (iret /= 0) goto 89 - - deallocate(tmp_input) - -!-------------------- -! Cloud liquid water -!-------------------- - - allocate(clwmr_b4_adj_output(ij_output,num_fields)) - clwmr_b4_adj_output = 0 - - print*,'INTERPOLATE CLOUD LIQUID WATER' - call ipolates(ip, ipopt, kgds_input, kgds_output, ij_input, ij_output,& - num_fields, ibi, bitmap_input, clwmr_input, & - numpts, rlat_output, rlon_output, ibo, bitmap_output, & - clwmr_b4_adj_output, iret) - if (iret /= 0) goto 89 - - deallocate(clwmr_input) - -!-------------------- -! Specific humidity -!-------------------- - - allocate(spfh_b4_adj_output(ij_output,num_fields)) - spfh_b4_adj_output = 0 - - print*,'INTERPOLATE SPECIFIC HUMIDITY' - call ipolates(ip, ipopt, kgds_input, kgds_output, ij_input, ij_output,& - num_fields, ibi, bitmap_input, spfh_input, & - numpts, rlat_output, rlon_output, ibo, bitmap_output, & - spfh_b4_adj_output, iret) - if (iret /= 0) goto 89 - - deallocate(spfh_input) - -!----------- -! Ozone -!----------- - - allocate(o3mr_b4_adj_output(ij_output,num_fields)) - o3mr_b4_adj_output = 0 - - print*,'INTERPOLATE OZONE' - call ipolates(ip, ipopt, kgds_input, kgds_output, ij_input, ij_output,& - num_fields, ibi, bitmap_input, o3mr_input, & - numpts, rlat_output, rlon_output, ibo, bitmap_output, & - o3mr_b4_adj_output, iret) - if (iret /= 0) goto 89 - - deallocate(o3mr_input) - -!----------- -! DZDT -!----------- - - allocate(dzdt_b4_adj_output(ij_output,num_fields)) - dzdt_b4_adj_output = 0 - - print*,'INTERPOLATE DZDT' - call ipolates(ip, ipopt, kgds_input, kgds_output, ij_input, ij_output,& - num_fields, ibi, bitmap_input, dzdt_input, & - numpts, rlat_output, rlon_output, ibo, bitmap_output, & - dzdt_b4_adj_output, iret) - if (iret /= 0) goto 89 - - deallocate(dzdt_input) - -!---------------------------------------------------------------------------------- -! Interpolate additional 3-d scalars for GFDL microphysics. -!---------------------------------------------------------------------------------- - - if (gfdl_mp) then - -!------------- -! Rain water -!------------- - - allocate(rwmr_b4_adj_output(ij_output,num_fields)) - rwmr_b4_adj_output = 0 - - print*,'INTERPOLATE RWMR' - call ipolates(ip, ipopt, kgds_input, kgds_output, ij_input, ij_output,& - num_fields, ibi, bitmap_input, rwmr_input, & - numpts, rlat_output, rlon_output, ibo, bitmap_output, & - rwmr_b4_adj_output, iret) - if (iret /= 0) goto 89 - - deallocate(rwmr_input) - -!------------- -! Snow water -!------------- - - allocate(snmr_b4_adj_output(ij_output,num_fields)) - snmr_b4_adj_output = 0 - - print*,'INTERPOLATE SNMR' - call ipolates(ip, ipopt, kgds_input, kgds_output, ij_input, ij_output,& - num_fields, ibi, bitmap_input, snmr_input, & - numpts, rlat_output, rlon_output, ibo, bitmap_output, & - snmr_b4_adj_output, iret) - if (iret /= 0) goto 89 - - deallocate(snmr_input) - -!------------- -! Ice water -!------------- - - allocate(icmr_b4_adj_output(ij_output,num_fields)) - icmr_b4_adj_output = 0 - - print*,'INTERPOLATE ICMR' - call ipolates(ip, ipopt, kgds_input, kgds_output, ij_input, ij_output,& - num_fields, ibi, bitmap_input, icmr_input, & - numpts, rlat_output, rlon_output, ibo, bitmap_output, & - icmr_b4_adj_output, iret) - if (iret /= 0) goto 89 - - deallocate(icmr_input) - -!------------- -! Graupel -!------------- - - allocate(grle_b4_adj_output(ij_output,num_fields)) - grle_b4_adj_output = 0 - - print*,'INTERPOLATE GRLE' - call ipolates(ip, ipopt, kgds_input, kgds_output, ij_input, ij_output,& - num_fields, ibi, bitmap_input, grle_input, & - numpts, rlat_output, rlon_output, ibo, bitmap_output, & - grle_b4_adj_output, iret) - if (iret /= 0) goto 89 - - deallocate(grle_input) - -!--------------------------- -! Cloud amount (if present) -!--------------------------- - - if (icldamt == 1) then - allocate(cldamt_b4_adj_output(ij_output,num_fields)) - cldamt_b4_adj_output = 0 - - print*,'INTERPOLATE CLD_AMT' - call ipolates(ip, ipopt, kgds_input, kgds_output, ij_input, ij_output,& - num_fields, ibi, bitmap_input, cldamt_input, & - numpts, rlat_output, rlon_output, ibo, bitmap_output, & - cldamt_b4_adj_output, iret) - if (iret /= 0) goto 89 - - deallocate(cldamt_input) - endif - - - endif - -!---------------------------------------------------------------------------------- -! 3d u/v winds -!---------------------------------------------------------------------------------- - - allocate(crot(ij_output), srot(ij_output)) - crot = 0. - srot = 0. - - allocate(ugrd_b4_adj_output(ij_output,num_fields)) - ugrd_b4_adj_output = 0 - allocate(vgrd_b4_adj_output(ij_output,num_fields)) - vgrd_b4_adj_output = 0 - - print*,'INTERPOLATE WINDS' - call ipolatev(ip, ipopt, kgds_input, kgds_output, ij_input, ij_output,& - num_fields, ibi, bitmap_input, ugrd_input, vgrd_input, & - numpts, rlat_output, rlon_output, crot, srot, ibo, bitmap_output, & - ugrd_b4_adj_output, vgrd_b4_adj_output, iret) - if (iret /= 0) goto 89 - - deallocate (ugrd_input, vgrd_input) - deallocate (crot, srot) - deallocate (ibi, ibo, bitmap_input, bitmap_output) - - return - - 89 continue - print*,"FATAL ERROR IN IPOLATES. IRET IS: ", iret - call errexit(23) - - end subroutine gaus_to_gaus - - end module interp diff --git a/src/enkf_chgres_recenter.fd/output_data.f90 b/src/enkf_chgres_recenter.fd/output_data.f90 deleted file mode 100644 index 36063d3a..00000000 --- a/src/enkf_chgres_recenter.fd/output_data.f90 +++ /dev/null @@ -1,396 +0,0 @@ - module output_data - - use nemsio_module - - implicit none - - private - - integer, public :: kgds_output(200) - -! data on the output grid. - real, allocatable, public :: hgt_output(:) ! interpolated from input grid - real, allocatable, public :: hgt_external_output(:) - real, allocatable, public :: sfcp_output(:) - real, allocatable, public :: tmp_output(:,:) - real, allocatable, public :: clwmr_output(:,:) - real, allocatable, public :: delz_output(:,:) - real, allocatable, public :: dpres_output(:,:) - real, allocatable, public :: dzdt_output(:,:) - real, allocatable, public :: o3mr_output(:,:) - real, allocatable, public :: spfh_output(:,:) - real, allocatable, public :: ugrd_output(:,:) - real, allocatable, public :: vgrd_output(:,:) - real, allocatable, public :: rwmr_output(:,:) - real, allocatable, public :: icmr_output(:,:) - real, allocatable, public :: snmr_output(:,:) - real, allocatable, public :: grle_output(:,:) - real, allocatable, public :: cldamt_output(:,:) - real, allocatable, public :: rlat_output(:) - real, allocatable, public :: rlon_output(:) - - public :: set_output_grid - public :: write_output_data - - character(len=50), allocatable :: recname(:) - character(len=50), allocatable :: reclevtyp(:) - - integer(nemsio_intkind) :: nrec - integer(nemsio_intkind), allocatable :: reclev(:) - - real(nemsio_realkind), allocatable :: vcoord_header(:,:,:) - real(nemsio_realkind), allocatable :: lat(:), lon(:) - - contains - - subroutine set_output_grid - -!------------------------------------------------------------------- -! Set grid specs on the output grid. -!------------------------------------------------------------------- - - use setup - use input_data - use utils - - implicit none - - character(len=20) :: vlevtyp, vname - - integer(nemsio_intkind) :: vlev - integer :: iret - - real(nemsio_realkind), allocatable :: dummy(:) - - type(nemsio_gfile) :: gfile - - print* - print*,"OUTPUT GRID I/J DIMENSIONS: ", i_output, j_output - -!------------------------------------------------------------------- -! Set the grib 1 grid description section, which is needed -! by the IPOLATES library. -!------------------------------------------------------------------- - - kgds_output = 0 - - call calc_kgds(i_output, j_output, kgds_output) - -!------------------------------------------------------------------- -! Read the terrain on the output grid. To ensure exact match, -! read it from an existing enkf nemsio restart file. -!------------------------------------------------------------------- - - call nemsio_init(iret) - - print* - print*,"OPEN OUTPUT GRID TERRAIN FILE: ", trim(terrain_file) - call nemsio_open(gfile, terrain_file, "read", iret=iret) - if (iret /= 0) then - print*,"FATAL ERROR OPENING FILE: ",trim(terrain_file) - print*,"IRET IS: ", iret - call errexit(50) - endif - - allocate(dummy(ij_output)) - allocate(hgt_external_output(ij_output)) - - print* - print*,"READ SURFACE HEIGHT" - vlev = 1 - vlevtyp = "sfc" - vname = "hgt" - call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret) - if (iret /= 0) then - print*,"FATAL ERROR READING FILE: ",trim(terrain_file) - print*,"IRET IS: ", iret - call errexit(51) - endif - - hgt_external_output = dummy - - deallocate(dummy) - - call nemsio_close(gfile, iret=iret) - - call nemsio_finalize() - - end subroutine set_output_grid - - subroutine write_output_data - -!------------------------------------------------------------------- -! Write output grid data to a nemsio file. -!------------------------------------------------------------------- - - use input_data - use setup - - implicit none - - character(len=5) :: gaction - - integer :: n, iret - - real(nemsio_realkind), allocatable :: dummy(:) - - type(nemsio_gfile) :: gfile - -!------------------------------------------------------------------- -! Set up some header info. -!------------------------------------------------------------------- - - call header_set - -!------------------------------------------------------------------- -! Open and write file. -!------------------------------------------------------------------- - - call nemsio_init(iret) - - gaction="write" - - print* - print*,'OPEN OUTPUT FILE: ',trim(output_file) - call nemsio_open(gfile, output_file, gaction, iret=iret, gdatatype="bin4", & - nmeta=8, modelname="FV3GFS", nrec=nrec, & - idate=idate, dimx=i_output, & - dimy=j_output, dimz=lev, ntrac=ntrac, & - ncldt=ncldt, idvc=idvc, idsl=idsl, idvm=idvm, & - idrt=4, recname=recname, reclevtyp=reclevtyp, & - reclev=reclev,vcoord=vcoord_header, & - lat=lat, lon=lon) - if (iret/=0) then - print*,"FATAL ERROR OPENING FILE. IRET IS: ", iret - call errexit(9) - endif - - deallocate(lon, lat, recname, reclevtyp, reclev, vcoord_header) - - allocate(dummy(i_output*j_output)) - - print*,"WRITE SURFACE HEIGHT" - dummy = hgt_external_output - call nemsio_writerecv(gfile, "hgt", "sfc", 1, dummy, iret=iret) - if (iret/=0) goto 88 - deallocate(hgt_external_output) - - print*,"WRITE SURFACE PRESSURE" - dummy = sfcp_output - call nemsio_writerecv(gfile, "pres", "sfc", 1, dummy, iret=iret) - if (iret/=0) goto 88 - deallocate(sfcp_output) - - print*,"WRITE TEMPERATURE" - do n = 1, lev - dummy = tmp_output(:,n) - call nemsio_writerecv(gfile, "tmp", "mid layer", n, dummy, iret=iret) - if (iret/=0) goto 88 - enddo - deallocate(tmp_output) - - print*,"WRITE CLOUD LIQUID WATER" - do n = 1, lev - dummy = clwmr_output(:,n) - call nemsio_writerecv(gfile, "clwmr", "mid layer", n, dummy, iret=iret) - if (iret/=0) goto 88 - enddo - deallocate(clwmr_output) - - print*,"WRITE SPECIFIC HUMIDITY" - do n = 1, lev - dummy = spfh_output(:,n) - call nemsio_writerecv(gfile, "spfh", "mid layer", n, dummy, iret=iret) - if (iret/=0) goto 88 - enddo - deallocate(spfh_output) - - print*,"WRITE OZONE" - do n = 1, lev - dummy = o3mr_output(:,n) - call nemsio_writerecv(gfile, "o3mr", "mid layer", n, dummy, iret=iret) - if (iret/=0) goto 88 - enddo - deallocate(o3mr_output) - - print*,"WRITE U-WINDS" - do n = 1, lev - dummy = ugrd_output(:,n) - call nemsio_writerecv(gfile, "ugrd", "mid layer", n, dummy, iret=iret) - if (iret/=0) goto 88 - enddo - deallocate(ugrd_output) - - print*,"WRITE V-WINDS" - do n = 1, lev - dummy = vgrd_output(:,n) - call nemsio_writerecv(gfile, "vgrd", "mid layer", n, dummy, iret=iret) - if (iret/=0) goto 88 - enddo - deallocate(vgrd_output) - - print*,"WRITE DZDT" - do n = 1, lev - dummy = dzdt_output(:,n) - call nemsio_writerecv(gfile, "dzdt", "mid layer", n, dummy, iret=iret) - if (iret/=0) goto 88 - enddo - deallocate(dzdt_output) - - print*,"WRITE DPRES" - do n = 1, lev - dummy = dpres_output(:,n) - call nemsio_writerecv(gfile, "dpres", "mid layer", n, dummy, iret=iret) - if (iret/=0) goto 88 - enddo - deallocate(dpres_output) - - print*,"WRITE DELZ" - do n = 1, lev - dummy = delz_output(:,n) - call nemsio_writerecv(gfile, "delz", "mid layer", n, dummy, iret=iret) - if (iret/=0) goto 88 - enddo - deallocate(delz_output) - - if (gfdl_mp) then - - print*,"WRITE RAIN WATER" - do n = 1, lev - dummy = rwmr_output(:,n) - call nemsio_writerecv(gfile, "rwmr", "mid layer", n, dummy, iret=iret) - if (iret/=0) goto 88 - enddo - deallocate(rwmr_output) - - print*,"WRITE SNOW WATER" - do n = 1, lev - dummy = snmr_output(:,n) - call nemsio_writerecv(gfile, "snmr", "mid layer", n, dummy, iret=iret) - if (iret/=0) goto 88 - enddo - deallocate(snmr_output) - - print*,"WRITE ICE WATER" - do n = 1, lev - dummy = icmr_output(:,n) - call nemsio_writerecv(gfile, "icmr", "mid layer", n, dummy, iret=iret) - if (iret/=0) goto 88 - enddo - deallocate(icmr_output) - - print*,"WRITE GRAUPEL" - do n = 1, lev - dummy = grle_output(:,n) - call nemsio_writerecv(gfile, "grle", "mid layer", n, dummy, iret=iret) - if (iret/=0) goto 88 - enddo - deallocate(grle_output) - - if (icldamt == 1) then - print*,"WRITE CLD_AMT" - do n = 1, lev - dummy = cldamt_output(:,n) - call nemsio_writerecv(gfile, "cld_amt", "mid layer", n, dummy, iret=iret) - if (iret/=0) goto 88 - enddo - deallocate(cldamt_output) - endif - - - endif - - deallocate(dummy) - - call nemsio_close(gfile, iret=iret) - - call nemsio_finalize() - - return - - 88 continue - print*,"FATAL ERROR WRITING FILE. IRET IS: ", iret - call errexit(10) - - end subroutine write_output_data - - subroutine header_set - -!------------------------------------------------------------------- -! Set header information for the output nemsio file. -!------------------------------------------------------------------- - - use input_data - use setup - - implicit none - - character(len=8) :: fields(9) - character(len=8) :: fields_gfdl_mp(5) - - integer :: count, l, n - -! Fields common to Zhao-Carr and GFDL microphysics - data fields /'ugrd', 'vgrd', 'dzdt', 'dpres', 'delz', & - 'tmp', 'spfh', 'clwmr', 'o3mr'/ - -! Fields for GFDL microphysics - data fields_gfdl_mp /'rwmr', 'icmr', 'snmr', 'grle', 'cld_amt'/ - - print* - print*,"SET HEADER INFO FOR OUTPUT FILE." - - if (gfdl_mp) then - nrec = ((13+icldamt) * lev) + 2 - else - nrec = (9 * lev) + 2 - endif - - allocate(recname(nrec)) - allocate(reclev(nrec)) - allocate(reclevtyp(nrec)) - - count = 0 - do n = 1, 9 - do l = 1, lev - count = count + 1 - recname(count) = fields(n) - reclev(count) = l - reclevtyp(count) = "mid layer" - enddo - enddo - - if (gfdl_mp) then - do n = 1, 4 + icldamt - do l = 1, lev - count = count + 1 - recname(count) = fields_gfdl_mp(n) - reclev(count) = l - reclevtyp(count) = "mid layer" - enddo - enddo - endif - - recname(nrec-1) = "pres" - reclev(nrec-1) = 1 - reclevtyp(nrec-1) = "sfc" - - recname(nrec) = "hgt" - reclev(nrec) = 1 - reclevtyp(nrec) = "sfc" - - allocate(vcoord_header(lev+1,3,2)) - vcoord_header = 0.0 - vcoord_header(:,1,1) = vcoord(:,1) - vcoord_header(:,2,1) = vcoord(:,2) - - allocate(lat(ij_output), lon(ij_output)) - - lat = rlat_output - lon = rlon_output - - deallocate(rlat_output, rlon_output) - - end subroutine header_set - - end module output_data diff --git a/src/enkf_chgres_recenter.fd/setup.f90 b/src/enkf_chgres_recenter.fd/setup.f90 deleted file mode 100644 index c2c2dc45..00000000 --- a/src/enkf_chgres_recenter.fd/setup.f90 +++ /dev/null @@ -1,53 +0,0 @@ - module setup - - use nemsio_module - - implicit none - - private - - character(len=300), public :: input_file - character(len=300), public :: output_file - character(len=300), public :: terrain_file - character(len=300), public :: vcoord_file - - integer(nemsio_intkind), public :: i_output - integer(nemsio_intkind), public :: j_output - integer , public :: ij_output - logical, public :: flipdelz - - public :: program_setup - - contains - - subroutine program_setup - - implicit none - - integer :: istat - - namelist /nam_setup/ i_output, j_output, input_file, output_file, & - terrain_file, vcoord_file - - print* - print*,"OPEN SETUP NAMELIST." - open(43, file="./fort.43", iostat=istat) - if (istat /= 0) then - print*,"FATAL ERROR OPENING NAMELIST FILE. ISTAT IS: ",istat - call errexit(30) - endif - - print*,"READ SETUP NAMELIST." - read(43, nml=nam_setup, iostat=istat) - if (istat /= 0) then - print*,"FATAL ERROR READING NAMELIST FILE. ISTAT IS: ",istat - call errexit(31) - endif - - ij_output = i_output * j_output - - close(43) - - end subroutine program_setup - - end module setup diff --git a/src/enkf_chgres_recenter.fd/utils.f90 b/src/enkf_chgres_recenter.fd/utils.f90 deleted file mode 100644 index e09c75b0..00000000 --- a/src/enkf_chgres_recenter.fd/utils.f90 +++ /dev/null @@ -1,783 +0,0 @@ - module utils - - private - - public :: calc_kgds - public :: newps - public :: newpr1 - public :: vintg - public :: compute_delz - - contains - - subroutine compute_delz(ijm, levp, ak_in, bk_in, ps, zs, t, sphum, delz, flipsign) - - implicit none - integer, intent(in):: levp, ijm - real, intent(in), dimension(levp+1):: ak_in, bk_in - real, intent(in), dimension(ijm):: ps, zs - real, intent(in), dimension(ijm,levp):: t - real, intent(in), dimension(ijm,levp):: sphum - real, intent(out), dimension(ijm,levp):: delz - logical, intent(in) :: flipsign -! Local: - real, dimension(ijm,levp+1):: zh - real, dimension(ijm,levp+1):: pe0, pn0 - real, dimension(levp+1) :: ak, bk - integer i,k - real, parameter :: GRAV = 9.80665 - real, parameter :: RDGAS = 287.05 - real, parameter :: RVGAS = 461.50 - real :: zvir - real:: grd - - print*,"COMPUTE LAYER THICKNESS." - - grd = grav/rdgas - zvir = rvgas/rdgas - 1. - ak = ak_in - bk = bk_in - ak(levp+1) = max(1.e-9, ak(levp+1)) - - do i=1, ijm - pe0(i,levp+1) = ak(levp+1) - pn0(i,levp+1) = log(pe0(i,levp+1)) - enddo - - do k=levp,1, -1 - do i=1,ijm - pe0(i,k) = ak(k) + bk(k)*ps(i) - pn0(i,k) = log(pe0(i,k)) - enddo - enddo - - do i = 1, ijm - zh(i,1) = zs(i) - enddo - - do k = 2, levp+1 - do i = 1, ijm - zh(i,k) = zh(i,k-1)+t(i,k-1)*(1.+zvir*sphum(i,k-1))* & - (pn0(i,k-1)-pn0(i,k))/grd - enddo - enddo - - do k = 1, levp - do i = 1, ijm - if (flipsign) then - delz(i,k) = zh(i,k) - zh(i,k+1) - else - delz(i,k) = zh(i,k+1) - zh(i,k) - end if - enddo - enddo - - end subroutine compute_delz - - subroutine calc_kgds(idim, jdim, kgds) - - use nemsio_module - - implicit none - - integer(nemsio_intkind), intent(in) :: idim, jdim - - integer, intent(out) :: kgds(200) - - kgds = 0 - kgds(1) = 4 ! OCT 6 - TYPE OF GRID (GAUSSIAN) - kgds(2) = idim ! OCT 7-8 - # PTS ON LATITUDE CIRCLE - kgds(3) = jdim ! OCT 9-10 - # PTS ON LONGITUDE CIRCLE - kgds(4) = 90000 ! OCT 11-13 - LAT OF ORIGIN - kgds(5) = 0 ! OCT 14-16 - LON OF ORIGIN - kgds(6) = 128 ! OCT 17 - RESOLUTION FLAG - kgds(7) = -90000 ! OCT 18-20 - LAT OF EXTREME POINT - kgds(8) = nint(-360000./idim) ! OCT 21-23 - LON OF EXTREME POINT - kgds(9) = nint((360.0 / float(idim))*1000.0) - ! OCT 24-25 - LONGITUDE DIRECTION INCR. - kgds(10) = jdim/2 ! OCT 26-27 - NUMBER OF CIRCLES POLE TO EQUATOR - kgds(12) = 255 ! OCT 29 - RESERVED - kgds(20) = 255 ! OCT 5 - NOT USED, SET TO 255 - - end subroutine calc_kgds - - SUBROUTINE NEWPS(ZS,PS,IM,KM,P,T,Q,ZSNEW,PSNEW) -!$$$ SUBPROGRAM DOCUMENTATION BLOCK -! -! SUBPROGRAM: NEWPS COMPUTE NEW SURFACE PRESSURE -! PRGMMR: IREDELL ORG: W/NMC23 DATE: 92-10-31 -! -! ABSTRACT: COMPUTES A NEW SURFACE PRESSURE GIVEN A NEW OROGRAPHY. -! THE NEW PRESSURE IS COMPUTED ASSUMING A HYDROSTATIC BALANCE -! AND A CONSTANT TEMPERATURE LAPSE RATE. BELOW GROUND, THE -! LAPSE RATE IS ASSUMED TO BE -6.5 K/KM. -! -! PROGRAM HISTORY LOG: -! 91-10-31 MARK IREDELL -! -! USAGE: CALL NEWPS(ZS,PS,IM,KM,P,T,Q,ZSNEW,PSNEW) -! INPUT ARGUMENT LIST: -! IM INTEGER NUMBER OF POINTS TO COMPUTE -! ZS REAL (IM) OLD OROGRAPHY (M) -! PS REAL (IM) OLD SURFACE PRESSURE (PA) -! KM INTEGER NUMBER OF LEVELS -! P REAL (IM,KM) PRESSURES (PA) -! T REAL (IM,KM) TEMPERATURES (K) -! Q REAL (IM,KM) SPECIFIC HUMIDITIES (KG/KG) -! ZSNEW REAL (IM) NEW OROGRAPHY (M) -! OUTPUT ARGUMENT LIST: -! PSNEW REAL (IM) NEW SURFACE PRESSURE (PA) -! -! ATTRIBUTES: -! LANGUAGE: FORTRAN -! -!C$$$ - REAL ZS(IM),PS(IM),P(IM,KM),T(IM,KM),Q(IM,KM) - REAL ZSNEW(IM),PSNEW(IM) - PARAMETER(BETA=-6.5E-3,EPSILON=1.E-9) - PARAMETER(G=9.80665,RD=287.05,RV=461.50) - PARAMETER(GOR=G/RD,FV=RV/RD-1.) - REAL ZU(IM) - FTV(AT,AQ)=AT*(1+FV*AQ) - FGAM(APU,ATVU,APD,ATVD)=-GOR*LOG(ATVD/ATVU)/LOG(APD/APU) - FZ0(AP,ATV,AZD,APD)=AZD+ATV/GOR*LOG(APD/AP) - FZ1(AP,ATV,AZD,APD,AGAM)=AZD-ATV/AGAM*((APD/AP)**(-AGAM/GOR)-1) - FP0(AZ,AZU,APU,ATVU)=APU*EXP(-GOR/ATVU*(AZ-AZU)) - FP1(AZ,AZU,APU,ATVU,AGAM)=APU*(1+AGAM/ATVU*(AZ-AZU))**(-GOR/AGAM) -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -! COMPUTE SURFACE PRESSURE BELOW THE ORIGINAL GROUND - LS=0 - K=1 - GAMMA=BETA - DO I=1,IM - PU=P(I,K) - TVU=FTV(T(I,K),Q(I,K)) - ZU(I)=FZ1(PU,TVU,ZS(I),PS(I),GAMMA) - IF(ZSNEW(I).LE.ZU(I)) THEN - PU=P(I,K) - TVU=FTV(T(I,K),Q(I,K)) - IF(ABS(GAMMA).GT.EPSILON) THEN - PSNEW(I)=FP1(ZSNEW(I),ZU(I),PU,TVU,GAMMA) - ELSE - PSNEW(I)=FP0(ZSNEW(I),ZU(I),PU,TVU) - ENDIF - ELSE - PSNEW(I)=0 - LS=LS+1 - ENDIF -! endif - ENDDO -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -! COMPUTE SURFACE PRESSURE ABOVE THE ORIGINAL GROUND - DO K=2,KM - IF(LS.GT.0) THEN - DO I=1,IM - IF(PSNEW(I).EQ.0) THEN - PU=P(I,K) - TVU=FTV(T(I,K),Q(I,K)) - PD=P(I,K-1) - TVD=FTV(T(I,K-1),Q(I,K-1)) - GAMMA=FGAM(PU,TVU,PD,TVD) - IF(ABS(GAMMA).GT.EPSILON) THEN - ZU(I)=FZ1(PU,TVU,ZU(I),PD,GAMMA) - ELSE - ZU(I)=FZ0(PU,TVU,ZU(I),PD) - ENDIF - IF(ZSNEW(I).LE.ZU(I)) THEN - IF(ABS(GAMMA).GT.EPSILON) THEN - PSNEW(I)=FP1(ZSNEW(I),ZU(I),PU,TVU,GAMMA) - ELSE - PSNEW(I)=FP0(ZSNEW(I),ZU(I),PU,TVU) - ENDIF - LS=LS-1 - ENDIF - ENDIF - ENDDO - ENDIF - ENDDO -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -! COMPUTE SURFACE PRESSURE OVER THE TOP - IF(LS.GT.0) THEN - K=KM - GAMMA=0 - DO I=1,IM - IF(PSNEW(I).EQ.0) THEN - PU=P(I,K) - TVU=FTV(T(I,K),Q(I,K)) - PSNEW(I)=FP0(ZSNEW(I),ZU(I),PU,TVU) - ENDIF - ENDDO - ENDIF - END SUBROUTINE NEWPS - - SUBROUTINE NEWPR1(IM,KM,IDVC,IDSL,NVCOORD,VCOORD, & - PS,PM,DP) -!$$$ SUBPROGRAM DOCUMENTATION BLOCK -! -! SUBPROGRAM: NEWPR1 COMPUTE MODEL PRESSURES -! PRGMMR: JUANG ORG: W/NMC23 DATE: 2005-04-11 -! PRGMMR: Fanglin Yang ORG: W/NMC23 DATE: 2006-11-28 -! PRGMMR: S. Moorthi ORG: NCEP/EMC DATE: 2006-12-12 -! PRGMMR: S. Moorthi ORG: NCEP/EMC DATE: 2007-01-02 -! -! ABSTRACT: COMPUTE MODEL PRESSURES. -! -! PROGRAM HISTORY LOG: -! 2005-04-11 HANN_MING HENRY JUANG hybrid sigma, sigma-p, and sigma- -! -! USAGE: CALL NEWPR1(IM,IX,KM,KMP,IDVC,IDSL,NVCOORD,VCOORD,PP,TP,QP,P -! INPUT ARGUMENT LIST: -! IM INTEGER NUMBER OF POINTS TO COMPUTE -! KM INTEGER NUMBER OF LEVELS -! IDVC INTEGER VERTICAL COORDINATE ID -! (1 FOR SIGMA AND 2 FOR HYBRID) -! IDSL INTEGER TYPE OF SIGMA STRUCTURE -! (1 FOR PHILLIPS OR 2 FOR MEAN) -! NVCOORD INTEGER NUMBER OF VERTICAL COORDINATES -! VCOORD REAL (KM+1,NVCOORD) VERTICAL COORDINATE VALUES -! FOR IDVC=1, NVCOORD=1: SIGMA INTERFACE -! FOR IDVC=2, NVCOORD=2: HYBRID INTERFACE A AND B -! FOR IDVC=3, NVCOORD=3: JUANG GENERAL HYBRID INTERFACE -! AK REAL (KM+1) HYBRID INTERFACE A -! BK REAL (KM+1) HYBRID INTERFACE B -! PS REAL (IX) SURFACE PRESSURE (PA) -! OUTPUT ARGUMENT LIST: -! PM REAL (IX,KM) MID-LAYER PRESSURE (PA) -! DP REAL (IX,KM) LAYER DELTA PRESSURE (PA) -! -! ATTRIBUTES: -! LANGUAGE: FORTRAN -! -!C$$$ - IMPLICIT NONE - - INTEGER, INTENT(IN) :: IM, KM, NVCOORD, IDVC, IDSL - - REAL, INTENT(IN) :: VCOORD(KM+1,NVCOORD) - REAL, INTENT(IN) :: PS(IM) - - REAL, INTENT(OUT) :: PM(IM,KM) - REAL, OPTIONAL, INTENT(OUT) :: DP(IM,KM) - - REAL, PARAMETER :: RD=287.05, RV=461.50, CP=1004.6, & - ROCP=RD/CP, ROCP1=ROCP+1, ROCPR=1/ROCP, & - FV=RV/RD-1. - - INTEGER :: I, K - - REAL :: AK(KM+1), BK(KM+1), PI(IM,KM+1) - - IF(IDVC.EQ.2) THEN - DO K=1,KM+1 - AK(K) = VCOORD(K,1) - BK(K) = VCOORD(K,2) - PI(1:IM,K) = AK(K) + BK(K)*PS(1:IM) - ENDDO - ELSE - print*,'routine only works for idvc 2' - stop - ENDIF - - IF(IDSL.EQ.2) THEN - DO K=1,KM - PM(1:IM,K) = (PI(1:IM,K)+PI(1:IM,K+1))/2 - ENDDO - ELSE - DO K=1,KM - PM(1:IM,K) = ((PI(1:IM,K)**ROCP1-PI(1:IM,K+1)**ROCP1)/ & - (ROCP1*(PI(1:IM,K)-PI(1:IM,K+1))))**ROCPR - ENDDO - ENDIF - - IF(PRESENT(DP))THEN - DO K=1,KM - DO I=1,IM - DP(I,K) = PI(I,K) - PI(I,K+1) - ENDDO - ENDDO - ENDIF - - END SUBROUTINE NEWPR1 - - SUBROUTINE TERP3(IM,IXZ1,IXQ1,IXZ2,IXQ2,NM,NXQ1,NXQ2, & - KM1,KXZ1,KXQ1,Z1,Q1,KM2,KXZ2,KXQ2,Z2,Q2,J2) -!$$$ SUBPROGRAM DOCUMENTATION BLOCK -! -! SUBPROGRAM: TERP3 CUBICALLY INTERPOLATE IN ONE DIMENSION -! PRGMMR: IREDELL ORG: W/NMC23 DATE: 98-05-01 -! -! ABSTRACT: INTERPOLATE FIELD(S) IN ONE DIMENSION ALONG THE COLUMN(S). -! THE INTERPOLATION IS CUBIC LAGRANGIAN WITH A MONOTONIC CONSTRAINT -! IN THE CENTER OF THE DOMAIN. IN THE OUTER INTERVALS IT IS LINEAR. -! OUTSIDE THE DOMAIN, FIELDS ARE HELD CONSTANT. -! -! PROGRAM HISTORY LOG: -! 98-05-01 MARK IREDELL -! 1999-01-04 IREDELL USE ESSL SEARCH -! -! USAGE: CALL TERP3(IM,IXZ1,IXQ1,IXZ2,IXQ2,NM,NXQ1,NXQ2, -! & KM1,KXZ1,KXQ1,Z1,Q1,KM2,KXZ2,KXQ2,Z2,Q2,J2) -! INPUT ARGUMENT LIST: -! IM INTEGER NUMBER OF COLUMNS -! IXZ1 INTEGER COLUMN SKIP NUMBER FOR Z1 -! IXQ1 INTEGER COLUMN SKIP NUMBER FOR Q1 -! IXZ2 INTEGER COLUMN SKIP NUMBER FOR Z2 -! IXQ2 INTEGER COLUMN SKIP NUMBER FOR Q2 -! NM INTEGER NUMBER OF FIELDS PER COLUMN -! NXQ1 INTEGER FIELD SKIP NUMBER FOR Q1 -! NXQ2 INTEGER FIELD SKIP NUMBER FOR Q2 -! KM1 INTEGER NUMBER OF INPUT POINTS -! KXZ1 INTEGER POINT SKIP NUMBER FOR Z1 -! KXQ1 INTEGER POINT SKIP NUMBER FOR Q1 -! Z1 REAL (1+(IM-1)*IXZ1+(KM1-1)*KXZ1) -! INPUT COORDINATE VALUES IN WHICH TO INTERPOLATE -! (Z1 MUST BE STRICTLY MONOTONIC IN EITHER DIRECTION) -! Q1 REAL (1+(IM-1)*IXQ1+(KM1-1)*KXQ1+(NM-1)*NXQ1) -! INPUT FIELDS TO INTERPOLATE -! KM2 INTEGER NUMBER OF OUTPUT POINTS -! KXZ2 INTEGER POINT SKIP NUMBER FOR Z2 -! KXQ2 INTEGER POINT SKIP NUMBER FOR Q2 -! Z2 REAL (1+(IM-1)*IXZ2+(KM2-1)*KXZ2) -! OUTPUT COORDINATE VALUES TO WHICH TO INTERPOLATE -! (Z2 NEED NOT BE MONOTONIC) -! -! OUTPUT ARGUMENT LIST: -! Q2 REAL (1+(IM-1)*IXQ2+(KM2-1)*KXQ2+(NM-1)*NXQ2) -! OUTPUT INTERPOLATED FIELDS -! J2 REAL (1+(IM-1)*IXQ2+(KM2-1)*KXQ2+(NM-1)*NXQ2) -! OUTPUT INTERPOLATED FIELDS CHANGE WRT Z2 -! -! SUBPROGRAMS CALLED: -! RSEARCH SEARCH FOR A SURROUNDING REAL INTERVAL -! -! ATTRIBUTES: -! LANGUAGE: FORTRAN -! -!C$$$ - IMPLICIT NONE - INTEGER IM,IXZ1,IXQ1,IXZ2,IXQ2,NM,NXQ1,NXQ2 - INTEGER KM1,KXZ1,KXQ1,KM2,KXZ2,KXQ2 - INTEGER I,K1,K2,N - REAL Z1(1+(IM-1)*IXZ1+(KM1-1)*KXZ1) - REAL Q1(1+(IM-1)*IXQ1+(KM1-1)*KXQ1+(NM-1)*NXQ1) - REAL Z2(1+(IM-1)*IXZ2+(KM2-1)*KXZ2) - REAL Q2(1+(IM-1)*IXQ2+(KM2-1)*KXQ2+(NM-1)*NXQ2) - REAL J2(1+(IM-1)*IXQ2+(KM2-1)*KXQ2+(NM-1)*NXQ2) - REAL FFA(IM),FFB(IM),FFC(IM),FFD(IM) - REAL GGA(IM),GGB(IM),GGC(IM),GGD(IM) - INTEGER K1S(IM,KM2) - REAL Z1A,Z1B,Z1C,Z1D,Q1A,Q1B,Q1C,Q1D,Z2S,Q2S,J2S -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -! FIND THE SURROUNDING INPUT INTERVAL FOR EACH OUTPUT POINT. - CALL RSEARCH(IM,KM1,IXZ1,KXZ1,Z1,KM2,IXZ2,KXZ2,Z2,1,IM,K1S) -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -! GENERALLY INTERPOLATE CUBICALLY WITH MONOTONIC CONSTRAINT -! FROM TWO NEAREST INPUT POINTS ON EITHER SIDE OF THE OUTPUT POINT, -! BUT WITHIN THE TWO EDGE INTERVALS INTERPOLATE LINEARLY. -! KEEP THE OUTPUT FIELDS CONSTANT OUTSIDE THE INPUT DOMAIN. - -!!$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(IM,IXZ1,IXQ1,IXZ2), & -!!$OMP& SHARED(IXQ2,NM,NXQ1,NXQ2,KM1,KXZ1,KXQ1,Z1,Q1,KM2,KXZ2), & -!!$OMP& SHARED(KXQ2,Z2,Q2,J2,K1S) - - DO K2=1,KM2 - DO I=1,IM - K1=K1S(I,K2) - IF(K1.EQ.1.OR.K1.EQ.KM1-1) THEN - Z2S=Z2(1+(I-1)*IXZ2+(K2-1)*KXZ2) - Z1A=Z1(1+(I-1)*IXZ1+(K1-1)*KXZ1) - Z1B=Z1(1+(I-1)*IXZ1+(K1+0)*KXZ1) - FFA(I)=(Z2S-Z1B)/(Z1A-Z1B) - FFB(I)=(Z2S-Z1A)/(Z1B-Z1A) - GGA(I)=1/(Z1A-Z1B) - GGB(I)=1/(Z1B-Z1A) - ELSEIF(K1.GT.1.AND.K1.LT.KM1-1) THEN - Z2S=Z2(1+(I-1)*IXZ2+(K2-1)*KXZ2) - Z1A=Z1(1+(I-1)*IXZ1+(K1-2)*KXZ1) - Z1B=Z1(1+(I-1)*IXZ1+(K1-1)*KXZ1) - Z1C=Z1(1+(I-1)*IXZ1+(K1+0)*KXZ1) - Z1D=Z1(1+(I-1)*IXZ1+(K1+1)*KXZ1) - FFA(I)=(Z2S-Z1B)/(Z1A-Z1B)* & - (Z2S-Z1C)/(Z1A-Z1C)* & - (Z2S-Z1D)/(Z1A-Z1D) - FFB(I)=(Z2S-Z1A)/(Z1B-Z1A)* & - (Z2S-Z1C)/(Z1B-Z1C)* & - (Z2S-Z1D)/(Z1B-Z1D) - FFC(I)=(Z2S-Z1A)/(Z1C-Z1A)* & - (Z2S-Z1B)/(Z1C-Z1B)* & - (Z2S-Z1D)/(Z1C-Z1D) - FFD(I)=(Z2S-Z1A)/(Z1D-Z1A)* & - (Z2S-Z1B)/(Z1D-Z1B)* & - (Z2S-Z1C)/(Z1D-Z1C) - GGA(I)= 1/(Z1A-Z1B)* & - (Z2S-Z1C)/(Z1A-Z1C)* & - (Z2S-Z1D)/(Z1A-Z1D)+ & - (Z2S-Z1B)/(Z1A-Z1B)* & - 1/(Z1A-Z1C)* & - (Z2S-Z1D)/(Z1A-Z1D)+ & - (Z2S-Z1B)/(Z1A-Z1B)* & - (Z2S-Z1C)/(Z1A-Z1C)* & - 1/(Z1A-Z1D) - GGB(I)= 1/(Z1B-Z1A)* & - (Z2S-Z1C)/(Z1B-Z1C)* & - (Z2S-Z1D)/(Z1B-Z1D)+ & - (Z2S-Z1A)/(Z1B-Z1A)* & - 1/(Z1B-Z1C)* & - (Z2S-Z1D)/(Z1B-Z1D)+ & - (Z2S-Z1A)/(Z1B-Z1A)* & - (Z2S-Z1C)/(Z1B-Z1C)* & - 1/(Z1B-Z1D) - GGC(I)= 1/(Z1C-Z1A)* & - (Z2S-Z1B)/(Z1C-Z1B)* & - (Z2S-Z1D)/(Z1C-Z1D)+ & - (Z2S-Z1A)/(Z1C-Z1A)* & - 1/(Z1C-Z1B)* & - (Z2S-Z1D)/(Z1C-Z1D)+ & - (Z2S-Z1A)/(Z1C-Z1A)* & - (Z2S-Z1B)/(Z1C-Z1B)* & - 1/(Z1C-Z1D) - GGD(I)= 1/(Z1D-Z1A)* & - (Z2S-Z1B)/(Z1D-Z1B)* & - (Z2S-Z1C)/(Z1D-Z1C)+ & - (Z2S-Z1A)/(Z1D-Z1A)* & - 1/(Z1D-Z1B)* & - (Z2S-Z1C)/(Z1D-Z1C)+ & - (Z2S-Z1A)/(Z1D-Z1A)* & - (Z2S-Z1B)/(Z1D-Z1B)* & - 1/(Z1D-Z1C) - ENDIF - ENDDO -! INTERPOLATE. - DO N=1,NM - DO I=1,IM - K1=K1S(I,K2) - IF(K1.EQ.0) THEN - Q2S=Q1(1+(I-1)*IXQ1+(N-1)*NXQ1) - J2S=0 - ELSEIF(K1.EQ.KM1) THEN - Q2S=Q1(1+(I-1)*IXQ1+(KM1-1)*KXQ1+(N-1)*NXQ1) - J2S=0 - ELSEIF(K1.EQ.1.OR.K1.EQ.KM1-1) THEN - Q1A=Q1(1+(I-1)*IXQ1+(K1-1)*KXQ1+(N-1)*NXQ1) - Q1B=Q1(1+(I-1)*IXQ1+(K1+0)*KXQ1+(N-1)*NXQ1) - Q2S=FFA(I)*Q1A+FFB(I)*Q1B - J2S=GGA(I)*Q1A+GGB(I)*Q1B - ELSE - Q1A=Q1(1+(I-1)*IXQ1+(K1-2)*KXQ1+(N-1)*NXQ1) - Q1B=Q1(1+(I-1)*IXQ1+(K1-1)*KXQ1+(N-1)*NXQ1) - Q1C=Q1(1+(I-1)*IXQ1+(K1+0)*KXQ1+(N-1)*NXQ1) - Q1D=Q1(1+(I-1)*IXQ1+(K1+1)*KXQ1+(N-1)*NXQ1) - Q2S=FFA(I)*Q1A+FFB(I)*Q1B+FFC(I)*Q1C+FFD(I)*Q1D - J2S=GGA(I)*Q1A+GGB(I)*Q1B+GGC(I)*Q1C+GGD(I)*Q1D - IF(Q2S.LT.MIN(Q1B,Q1C)) THEN - Q2S=MIN(Q1B,Q1C) - J2S=0 - ELSEIF(Q2S.GT.MAX(Q1B,Q1C)) THEN - Q2S=MAX(Q1B,Q1C) - J2S=0 - ENDIF - ENDIF - Q2(1+(I-1)*IXQ2+(K2-1)*KXQ2+(N-1)*NXQ2)=Q2S - J2(1+(I-1)*IXQ2+(K2-1)*KXQ2+(N-1)*NXQ2)=J2S - ENDDO - ENDDO - ENDDO -!!$OMP END PARALLEL DO - - END SUBROUTINE TERP3 - - SUBROUTINE RSEARCH(IM,KM1,IXZ1,KXZ1,Z1,KM2,IXZ2,KXZ2,Z2,IXL2,KXL2,& - L2) -!$$$ SUBPROGRAM DOCUMENTATION BLOCK -! -! SUBPROGRAM: RSEARCH SEARCH FOR A SURROUNDING REAL INTERVAL -! PRGMMR: IREDELL ORG: W/NMC23 DATE: 98-05-01 -! -! ABSTRACT: THIS SUBPROGRAM SEARCHES MONOTONIC SEQUENCES OF REAL NUMBERS -! FOR INTERVALS THAT SURROUND A GIVEN SEARCH SET OF REAL NUMBERS. -! THE SEQUENCES MAY BE MONOTONIC IN EITHER DIRECTION; THE REAL NUMBERS -! MAY BE SINGLE OR DOUBLE PRECISION; THE INPUT SEQUENCES AND SETS -! AND THE OUTPUT LOCATIONS MAY BE ARBITRARILY DIMENSIONED. -! -! PROGRAM HISTORY LOG: -! 1999-01-05 MARK IREDELL -! -! USAGE: CALL RSEARCH(IM,KM1,IXZ1,KXZ1,Z1,KM2,IXZ2,KXZ2,Z2,IXL2,KXL2, -! & L2) -! INPUT ARGUMENT LIST: -! IM INTEGER NUMBER OF SEQUENCES TO SEARCH -! KM1 INTEGER NUMBER OF POINTS IN EACH SEQUENCE -! IXZ1 INTEGER SEQUENCE SKIP NUMBER FOR Z1 -! KXZ1 INTEGER POINT SKIP NUMBER FOR Z1 -! Z1 REAL (1+(IM-1)*IXZ1+(KM1-1)*KXZ1) -! SEQUENCE VALUES TO SEARCH -! (Z1 MUST BE MONOTONIC IN EITHER DIRECTION) -! KM2 INTEGER NUMBER OF POINTS TO SEARCH FOR -! IN EACH RESPECTIVE SEQUENCE -! IXZ2 INTEGER SEQUENCE SKIP NUMBER FOR Z2 -! KXZ2 INTEGER POINT SKIP NUMBER FOR Z2 -! Z2 REAL (1+(IM-1)*IXZ2+(KM2-1)*KXZ2) -! SET OF VALUES TO SEARCH FOR -! (Z2 NEED NOT BE MONOTONIC) -! IXL2 INTEGER SEQUENCE SKIP NUMBER FOR L2 -! KXL2 INTEGER POINT SKIP NUMBER FOR L2 -! -! OUTPUT ARGUMENT LIST: -! L2 INTEGER (1+(IM-1)*IXL2+(KM2-1)*KXL2) -! INTERVAL LOCATIONS HAVING VALUES FROM 0 TO KM1 -! (Z2 WILL BE BETWEEN Z1(L2) AND Z1(L2+1)) -! -! SUBPROGRAMS CALLED: -! SBSRCH ESSL BINARY SEARCH -! DBSRCH ESSL BINARY SEARCH -! -! REMARKS: -! IF THE ARRAY Z1 IS DIMENSIONED (IM,KM1), THEN THE SKIP NUMBERS ARE -! IXZ1=1 AND KXZ1=IM; IF IT IS DIMENSIONED (KM1,IM), THEN THE SKIP -! NUMBERS ARE IXZ1=KM1 AND KXZ1=1; IF IT IS DIMENSIONED (IM,JM,KM1), -! THEN THE SKIP NUMBERS ARE IXZ1=1 AND KXZ1=IM*JM; ETCETERA. -! SIMILAR EXAMPLES APPLY TO THE SKIP NUMBERS FOR Z2 AND L2. -! -! RETURNED VALUES OF 0 OR KM1 INDICATE THAT THE GIVEN SEARCH VALUE -! IS OUTSIDE THE RANGE OF THE SEQUENCE. -! -! IF A SEARCH VALUE IS IDENTICAL TO ONE OF THE SEQUENCE VALUES -! THEN THE LOCATION RETURNED POINTS TO THE IDENTICAL VALUE. -! IF THE SEQUENCE IS NOT STRICTLY MONOTONIC AND A SEARCH VALUE IS -! IDENTICAL TO MORE THAN ONE OF THE SEQUENCE VALUES, THEN THE -! LOCATION RETURNED MAY POINT TO ANY OF THE IDENTICAL VALUES. -! -! TO BE EXACT, FOR EACH I FROM 1 TO IM AND FOR EACH K FROM 1 TO KM2, -! Z=Z2(1+(I-1)*IXZ2+(K-1)*KXZ2) IS THE SEARCH VALUE AND -! L=L2(1+(I-1)*IXL2+(K-1)*KXL2) IS THE LOCATION RETURNED. -! IF L=0, THEN Z IS LESS THAN THE START POINT Z1(1+(I-1)*IXZ1) -! FOR ASCENDING SEQUENCES (OR GREATER THAN FOR DESCENDING SEQUENCES). -! IF L=KM1, THEN Z IS GREATER THAN OR EQUAL TO THE END POINT -! Z1(1+(I-1)*IXZ1+(KM1-1)*KXZ1) FOR ASCENDING SEQUENCES -! (OR LESS THAN OR EQUAL TO FOR DESCENDING SEQUENCES). -! OTHERWISE Z IS BETWEEN THE VALUES Z1(1+(I-1)*IXZ1+(L-1)*KXZ1) AND -! Z1(1+(I-1)*IXZ1+(L-0)*KXZ1) AND MAY EQUAL THE FORMER. -! -! ATTRIBUTES: -! LANGUAGE: FORTRAN -! -!C$$$ -! IMPLICIT NONE -! INTEGER,INTENT(IN):: IM,KM1,IXZ1,KXZ1,KM2,IXZ2,KXZ2,IXL2,KXL2 -! REAL,INTENT(IN):: Z1(1+(IM-1)*IXZ1+(KM1-1)*KXZ1) -! REAL,INTENT(IN):: Z2(1+(IM-1)*IXZ2+(KM2-1)*KXZ2) -! INTEGER,INTENT(OUT):: L2(1+(IM-1)*IXL2+(KM2-1)*KXL2) -! INTEGER(4) INCX,N,INCY,M,INDX(KM2),RC(KM2),IOPT -! INTEGER I,K2 -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -! FIND THE SURROUNDING INPUT INTERVAL FOR EACH OUTPUT POINT. -! DO I=1,IM -! IF(Z1(1+(I-1)*IXZ1).LE.Z1(1+(I-1)*IXZ1+(KM1-1)*KXZ1)) THEN -! INPUT COORDINATE IS MONOTONICALLY ASCENDING. -! INCX=KXZ2 -! N=KM2 -! INCY=KXZ1 -! M=KM1 -! IOPT=1 -! IF(DIGITS(1.).LT.DIGITS(1._8)) THEN -! CALL SBSRCH(Z2(1+(I-1)*IXZ2),INCX,N, -! & Z1(1+(I-1)*IXZ1),INCY,M,INDX,RC,IOPT) -! ELSE -! CALL DBSRCH(Z2(1+(I-1)*IXZ2),INCX,N, -! & Z1(1+(I-1)*IXZ1),INCY,M,INDX,RC,IOPT) -! ENDIF -! DO K2=1,KM2 -! L2(1+(I-1)*IXL2+(K2-1)*KXL2)=INDX(K2)-RC(K2) -! ENDDO -! ELSE -! INPUT COORDINATE IS MONOTONICALLY DESCENDING. -! INCX=KXZ2 -! N=KM2 -! INCY=-KXZ1 -! M=KM1 -! IOPT=0 -! IF(DIGITS(1.).LT.DIGITS(1._8)) THEN -! CALL SBSRCH(Z2(1+(I-1)*IXZ2),INCX,N, -! & Z1(1+(I-1)*IXZ1),INCY,M,INDX,RC,IOPT) -! ELSE -! CALL DBSRCH(Z2(1+(I-1)*IXZ2),INCX,N, -! & Z1(1+(I-1)*IXZ1),INCY,M,INDX,RC,IOPT) -! ENDIF -! DO K2=1,KM2 -! L2(1+(I-1)*IXL2+(K2-1)*KXL2)=KM1+1-INDX(K2) -! ENDDO -! ENDIF -! ENDDO -! - IMPLICIT NONE - INTEGER,INTENT(IN):: IM,KM1,IXZ1,KXZ1,KM2,IXZ2,KXZ2,IXL2,KXL2 - REAL,INTENT(IN):: Z1(1+(IM-1)*IXZ1+(KM1-1)*KXZ1) - REAL,INTENT(IN):: Z2(1+(IM-1)*IXZ2+(KM2-1)*KXZ2) - INTEGER,INTENT(OUT):: L2(1+(IM-1)*IXL2+(KM2-1)*KXL2) - INTEGER I,K2,L - REAL Z -!C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!C FIND THE SURROUNDING INPUT INTERVAL FOR EACH OUTPUT POINT. - DO I=1,IM - IF(Z1(1+(I-1)*IXZ1).LE.Z1(1+(I-1)*IXZ1+(KM1-1)*KXZ1)) THEN -!C INPUT COORDINATE IS MONOTONICALLY ASCENDING. - DO K2=1,KM2 - Z=Z2(1+(I-1)*IXZ2+(K2-1)*KXZ2) - L=0 - DO - IF(Z.LT.Z1(1+(I-1)*IXZ1+L*KXZ1)) EXIT - L=L+1 - IF(L.EQ.KM1) EXIT - ENDDO - L2(1+(I-1)*IXL2+(K2-1)*KXL2)=L - ENDDO - ELSE -!C INPUT COORDINATE IS MONOTONICALLY DESCENDING. - DO K2=1,KM2 - Z=Z2(1+(I-1)*IXZ2+(K2-1)*KXZ2) - L=0 - DO - IF(Z.GT.Z1(1+(I-1)*IXZ1+L*KXZ1)) EXIT - L=L+1 - IF(L.EQ.KM1) EXIT - ENDDO - L2(1+(I-1)*IXL2+(K2-1)*KXL2)=L - ENDDO - ENDIF - ENDDO - - END SUBROUTINE RSEARCH - - SUBROUTINE VINTG(IM,KM1,KM2,NT,P1,U1,V1,T1,Q1,W1,P2, & - U2,V2,T2,Q2,W2) -!$$$ SUBPROGRAM DOCUMENTATION BLOCK -! -! SUBPROGRAM: VINTG VERTICALLY INTERPOLATE UPPER-AIR FIELDS -! PRGMMR: IREDELL ORG: W/NMC23 DATE: 92-10-31 -! -! ABSTRACT: VERTICALLY INTERPOLATE UPPER-AIR FIELDS. -! WIND, TEMPERATURE, HUMIDITY AND OTHER TRACERS ARE INTERPOLATED. -! THE INTERPOLATION IS CUBIC LAGRANGIAN IN LOG PRESSURE -! WITH A MONOTONIC CONSTRAINT IN THE CENTER OF THE DOMAIN. -! IN THE OUTER INTERVALS IT IS LINEAR IN LOG PRESSURE. -! OUTSIDE THE DOMAIN, FIELDS ARE GENERALLY HELD CONSTANT, -! EXCEPT FOR TEMPERATURE AND HUMIDITY BELOW THE INPUT DOMAIN, -! WHERE THE TEMPERATURE LAPSE RATE IS HELD FIXED AT -6.5 K/KM AND -! THE RELATIVE HUMIDITY IS HELD CONSTANT. -! -! PROGRAM HISTORY LOG: -! 91-10-31 MARK IREDELL -! -! USAGE: CALL VINTG(IM,KM1,KM2,NT,P1,U1,V1,T1,Q1,P2, -! & U2,V2,T2,Q2) -! INPUT ARGUMENT LIST: -! IM INTEGER NUMBER OF POINTS TO COMPUTE -! KM1 INTEGER NUMBER OF INPUT LEVELS -! KM2 INTEGER NUMBER OF OUTPUT LEVELS -! NT INTEGER NUMBER OF TRACERS -! P1 REAL (IM,KM1) INPUT PRESSURES -! ORDERED FROM BOTTOM TO TOP OF ATMOSPHERE -! U1 REAL (IM,KM1) INPUT ZONAL WIND -! V1 REAL (IM,KM1) INPUT MERIDIONAL WIND -! T1 REAL (IM,KM1) INPUT TEMPERATURE (K) -! Q1 REAL (IM,KM1,NT) INPUT TRACERS (HUMIDITY FIRST) -! P2 REAL (IM,KM2) OUTPUT PRESSURES -! OUTPUT ARGUMENT LIST: -! U2 REAL (IM,KM2) OUTPUT ZONAL WIND -! V2 REAL (IM,KM2) OUTPUT MERIDIONAL WIND -! T2 REAL (IM,KM2) OUTPUT TEMPERATURE (K) -! Q2 REAL (IM,KM2,NT) OUTPUT TRACERS (HUMIDITY FIRST) -! -! SUBPROGRAMS CALLED: -! TERP3 CUBICALLY INTERPOLATE IN ONE DIMENSION -! -! ATTRIBUTES: -! LANGUAGE: FORTRAN -! -!C$$$ - IMPLICIT NONE - - INTEGER, INTENT(IN) :: IM, KM1, KM2, NT - - REAL, INTENT(IN) :: P1(IM,KM1),U1(IM,KM1),V1(IM,KM1) - REAL, INTENT(IN) :: T1(IM,KM1),Q1(IM,KM1,NT) - REAL, INTENT(IN) :: W1(IM,KM1),P2(IM,KM2) - REAL, INTENT(OUT) :: U2(IM,KM2),V2(IM,KM2) - REAL, INTENT(OUT) :: T2(IM,KM2),Q2(IM,KM2,NT) - REAL, INTENT(OUT) :: W2(IM,KM2) - - REAL, PARAMETER :: DLTDZ=-6.5E-3*287.05/9.80665 - REAL, PARAMETER :: DLPVDRT=-2.5E6/461.50 - - INTEGER :: I, K, N - - REAL :: DZ - REAL,ALLOCATABLE :: Z1(:,:),Z2(:,:) - REAL,ALLOCATABLE :: C1(:,:,:),C2(:,:,:),J2(:,:,:) - - ALLOCATE (Z1(IM+1,KM1),Z2(IM+1,KM2)) - ALLOCATE (C1(IM+1,KM1,4+NT),C2(IM+1,KM2,4+NT),J2(IM+1,KM2,4+NT)) -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -! COMPUTE LOG PRESSURE INTERPOLATING COORDINATE -! AND COPY INPUT WIND, TEMPERATURE, HUMIDITY AND OTHER TRACERS -!$OMP PARALLEL DO PRIVATE(K,I) - DO K=1,KM1 - DO I=1,IM - Z1(I,K) = -LOG(P1(I,K)) - C1(I,K,1) = U1(I,K) - C1(I,K,2) = V1(I,K) - C1(I,K,3) = W1(I,K) - C1(I,K,4) = T1(I,K) - C1(I,K,5) = Q1(I,K,1) - ENDDO - ENDDO -!$OMP END PARALLEL DO - DO N=2,NT - DO K=1,KM1 - DO I=1,IM - C1(I,K,4+N) = Q1(I,K,N) - ENDDO - ENDDO - ENDDO -!$OMP PARALLEL DO PRIVATE(K,I) - DO K=1,KM2 - DO I=1,IM - Z2(I,K) = -LOG(P2(I,K)) - ENDDO - ENDDO -!$OMP END PARALLEL DO -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -! PERFORM LAGRANGIAN ONE-DIMENSIONAL INTERPOLATION -! THAT IS 4TH-ORDER IN INTERIOR, 2ND-ORDER IN OUTSIDE INTERVALS -! AND 1ST-ORDER FOR EXTRAPOLATION. - CALL TERP3(IM,1,1,1,1,4+NT,(IM+1)*KM1,(IM+1)*KM2, & - KM1,IM+1,IM+1,Z1,C1,KM2,IM+1,IM+1,Z2,C2,J2) -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -! COPY OUTPUT WIND, TEMPERATURE, HUMIDITY AND OTHER TRACERS -! EXCEPT BELOW THE INPUT DOMAIN, LET TEMPERATURE INCREASE WITH A FIXED -! LAPSE RATE AND LET THE RELATIVE HUMIDITY REMAIN CONSTANT. - DO K=1,KM2 - DO I=1,IM - U2(I,K)=C2(I,K,1) - V2(I,K)=C2(I,K,2) - W2(I,K)=C2(I,K,3) - DZ=Z2(I,K)-Z1(I,1) - IF(DZ.GE.0) THEN - T2(I,K)=C2(I,K,4) - Q2(I,K,1)=C2(I,K,5) - ELSE - T2(I,K)=T1(I,1)*EXP(DLTDZ*DZ) - Q2(I,K,1)=Q1(I,1,1)*EXP(DLPVDRT*(1/T2(I,K)-1/T1(I,1))-DZ) - ENDIF - ENDDO - ENDDO - DO N=2,NT - DO K=1,KM2 - DO I=1,IM - Q2(I,K,N)=C2(I,K,4+N) - ENDDO - ENDDO - ENDDO - DEALLOCATE (Z1,Z2,C1,C2,J2) - END SUBROUTINE VINTG - end module utils diff --git a/src/fv3nc2nemsio.fd/0readme b/src/fv3nc2nemsio.fd/0readme deleted file mode 100644 index 7be2fbcd..00000000 --- a/src/fv3nc2nemsio.fd/0readme +++ /dev/null @@ -1,23 +0,0 @@ -The first version of this program was provided by Jeff Whitaker and Philip Pegion from ESRL. -Fanglin Ynag has subsequently made a few revsions. - -10/20/2016, Fanglin Yang -Note that FV3 lat-lon grids are located at the center of each grid box, -start from south to north, and from east to west. -For example, for a 0.5-deg uniform grid, -nlon=720, nlat=360 -X(1,1)=[0.25E,89.75S] -X(nlon,nlat)=[359.75E,89.75N] - -write out nemsio, S->N is reversed to N->S to follow NCEP convention - -12/18/2016 Fanglin Yang -updated to handle output of any frequency and any accumulation bucket - - -01/10/2017 Fanglin Yang -updated to handle both hydrostatic and nonhydrostatic cases. They have different output numbers and variable names. - -10/07/2017 Fanglin Yang -In FV3 tic26 branch which includes the lastest Write Component, hgtsfc has been defined as [m] instead of [gpm]. -The scaling by 1/grav in fv3nc2nemsio.fd needs to be removed. diff --git a/src/fv3nc2nemsio.fd/CMakeLists.txt b/src/fv3nc2nemsio.fd/CMakeLists.txt deleted file mode 100644 index 3c6d6171..00000000 --- a/src/fv3nc2nemsio.fd/CMakeLists.txt +++ /dev/null @@ -1,15 +0,0 @@ -list(APPEND fortran_src -constants.f90 -fv3_main.f90 -fv3_module.f90 -kinds.f90 -) - -set(exe_name fv3nc2nemsio.x) -add_executable(${exe_name} ${fortran_src}) -target_link_libraries(${exe_name} PRIVATE NetCDF::NetCDF_Fortran - bacio::bacio_4 - w3emc::w3emc_d - nemsio::nemsio) - -install(TARGETS ${exe_name} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/src/fv3nc2nemsio.fd/constants.f90 b/src/fv3nc2nemsio.fd/constants.f90 deleted file mode 100644 index c0a066ee..00000000 --- a/src/fv3nc2nemsio.fd/constants.f90 +++ /dev/null @@ -1,314 +0,0 @@ -! this module was extracted from the GSI version operational -! at NCEP in Dec. 2007. -module constants -!$$$ module documentation block -! . . . . -! module: constants -! prgmmr: treadon org: np23 date: 2003-09-25 -! -! abstract: This module contains the definition of various constants -! used in the gsi code -! -! program history log: -! 2003-09-25 treadon - original code -! 2004-03-02 treadon - allow global and regional constants to differ -! 2004-06-16 treadon - update documentation -! 2004-10-28 treadon - replace parameter tiny=1.e-12 with tiny_r_kind -! and tiny_single -! 2004-11-16 treadon - add huge_single, huge_r_kind parameters -! 2005-01-27 cucurull - add ione -! 2005-08-24 derber - move cg_term to constants from qcmod -! 2006-03-07 treadon - add rd_over_cp_mass -! 2006-05-18 treadon - add huge_i_kind -! 2006-06-06 su - add var-qc wgtlim, change value to 0.25 (ECMWF) -! 2006-07-28 derber - add r1000 -! -! Subroutines Included: -! sub init_constants - compute derived constants, set regional/global constants -! -! Variable Definitions: -! see below -! -! attributes: -! language: f90 -! machine: ibm RS/6000 SP -! -!$$$ - use kinds, only: r_single,r_kind,i_kind - implicit none - -! Declare constants - integer(i_kind) izero,ione - real(r_kind) rearth,grav,omega,rd,rv,cp,cv,cvap,cliq - real(r_kind) csol,hvap,hfus,psat,t0c,ttp,jcal,cp_mass,cg_term - real(r_kind) fv,deg2rad,rad2deg,pi,tiny_r_kind,huge_r_kind,huge_i_kind - real(r_kind) ozcon,rozcon,tpwcon,rd_over_g,rd_over_cp,g_over_rd - real(r_kind) amsua_clw_d1,amsua_clw_d2,constoz,zero,one,two,four - real(r_kind) one_tenth,quarter,three,five,rd_over_cp_mass, gamma - real(r_kind) rearth_equator,stndrd_atmos_ps,r1000 - real(r_kind) semi_major_axis,semi_minor_axis,n_a,n_b - real(r_kind) eccentricity,grav_polar,grav_ratio - real(r_kind) grav_equator,earth_omega,grav_constant - real(r_kind) flattening,eccentricity_linear,somigliana - real(r_kind) dldt,dldti,hsub,psatk,tmix,xa,xai,xb,xbi - real(r_kind) eps,epsm1,omeps,wgtlim - real(r_kind) elocp,cpr,el2orc,cclimit,climit,epsq - real(r_kind) pcpeff0,pcpeff1,pcpeff2,pcpeff3,rcp,c0,delta - real(r_kind) h1000,factor1,factor2,rhcbot,rhctop,dx_max,dx_min,dx_inv - real(r_kind) h300,half,cmr,cws,ke2,row,rrow - real(r_single) zero_single,tiny_single,huge_single - real(r_single) rmw_mean_distance, roic_mean_distance - logical :: constants_initialized = .true. - - -! Define constants common to global and regional applications -! name value description units -! ---- ----- ----------- ----- - parameter(rearth_equator= 6.37813662e6_r_kind) ! equatorial earth radius (m) - parameter(omega = 7.2921e-5_r_kind) ! angular velocity of earth (1/s) - parameter(cp = 1.0046e+3_r_kind) ! specific heat of air @pressure (J/kg/K) - parameter(cvap = 1.8460e+3_r_kind) ! specific heat of h2o vapor (J/kg/K) - parameter(csol = 2.1060e+3_r_kind) ! specific heat of solid h2o (ice)(J/kg/K) - parameter(hvap = 2.5000e+6_r_kind) ! latent heat of h2o condensation (J/kg) - parameter(hfus = 3.3358e+5_r_kind) ! latent heat of h2o fusion (J/kg) - parameter(psat = 6.1078e+2_r_kind) ! pressure at h2o triple point (Pa) - parameter(t0c = 2.7315e+2_r_kind) ! temperature at zero celsius (K) - parameter(ttp = 2.7316e+2_r_kind) ! temperature at h2o triple point (K) - parameter(jcal = 4.1855e+0_r_kind) ! joules per calorie () - parameter(stndrd_atmos_ps = 1013.25e2_r_kind) ! 1976 US standard atmosphere ps (Pa) - -! Numeric constants - parameter(izero = 0) - parameter(ione = 1) - parameter(zero_single = 0.0_r_single) - parameter(zero = 0.0_r_kind) - parameter(one_tenth = 0.10_r_kind) - parameter(quarter= 0.25_r_kind) - parameter(one = 1.0_r_kind) - parameter(two = 2.0_r_kind) - parameter(three = 3.0_r_kind) - parameter(four = 4.0_r_kind) - parameter(five = 5.0_r_kind) - parameter(r1000 = 1000.0_r_kind) - -! Constants for gps refractivity - parameter(n_a=77.6_r_kind) !K/mb - parameter(n_b=3.73e+5_r_kind) !K^2/mb - -! Parameters below from WGS-84 model software inside GPS receivers. - parameter(semi_major_axis = 6378.1370e3_r_kind) ! (m) - parameter(semi_minor_axis = 6356.7523142e3_r_kind) ! (m) - parameter(grav_polar = 9.8321849378_r_kind) ! (m/s2) - parameter(grav_equator = 9.7803253359_r_kind) ! (m/s2) - parameter(earth_omega = 7.292115e-5_r_kind) ! (rad/s) - parameter(grav_constant = 3.986004418e14_r_kind) ! (m3/s2) - -! Derived geophysical constants - parameter(flattening = (semi_major_axis-semi_minor_axis)/semi_major_axis)!() - parameter(somigliana = & - (semi_minor_axis/semi_major_axis) * (grav_polar/grav_equator) - one)!() - parameter(grav_ratio = (earth_omega*earth_omega * & - semi_major_axis*semi_major_axis * semi_minor_axis) / grav_constant) !() - -! Derived thermodynamic constants - parameter ( dldti = cvap-csol ) - parameter ( hsub = hvap+hfus ) - parameter ( psatk = psat*0.001_r_kind ) - parameter ( tmix = ttp-20._r_kind ) - parameter ( elocp = hvap/cp ) - parameter ( rcp = one/cp ) - -! Constants used in GFS moist physics - parameter ( h300 = 300._r_kind ) - parameter ( half = 0.5_r_kind ) - parameter ( cclimit = 0.001_r_kind ) - parameter ( climit = 1.e-20_r_kind) - parameter ( epsq = 2.e-12_r_kind ) - parameter ( h1000 = 1000.0_r_kind) - parameter ( rhcbot=0.85_r_kind ) - parameter ( rhctop=0.85_r_kind ) - parameter ( dx_max=-8.8818363_r_kind ) - parameter ( dx_min=-5.2574954_r_kind ) - parameter ( dx_inv=one/(dx_max-dx_min) ) - parameter ( c0=0.002_r_kind ) - parameter ( delta=0.6077338_r_kind ) - parameter ( pcpeff0=1.591_r_kind ) - parameter ( pcpeff1=-0.639_r_kind ) - parameter ( pcpeff2=0.0953_r_kind ) - parameter ( pcpeff3=-0.00496_r_kind ) - parameter ( cmr = one/0.0003_r_kind ) - parameter ( cws = 0.025_r_kind ) - parameter ( ke2 = 0.00002_r_kind ) - parameter ( row = 1000._r_kind ) - parameter ( rrow = one/row ) - -! Constant used to process ozone - parameter ( constoz = 604229.0_r_kind) - -! Constants used in cloud liquid water correction for AMSU-A -! brightness temperatures - parameter ( amsua_clw_d1 = 0.754_r_kind ) - parameter ( amsua_clw_d2 = -2.265_r_kind ) - -! Constants used for variational qc - parameter ( wgtlim = 0.25_r_kind) ! Cutoff weight for concluding that obs has been - ! rejected by nonlinear qc. This limit is arbitrary - ! and DOES NOT affect nonlinear qc. It only affects - ! the printout which "counts" the number of obs that - ! "fail" nonlinear qc. Observations counted as failing - ! nonlinear qc are still assimilated. Their weight - ! relative to other observations is reduced. Changing - ! wgtlim does not alter the analysis, only - ! the nonlinear qc data "count" - -! Constants describing the Extended Best-Track Reanalysis [Demuth et -! al., 2008] tropical cyclone (TC) distance for regions relative to TC -! track position; units are in kilometers - - parameter (rmw_mean_distance = 64.5479412) - parameter (roic_mean_distance = 338.319656) - -contains - subroutine init_constants_derived -!$$$ subprogram documentation block -! . . . . -! subprogram: init_constants_derived set derived constants -! prgmmr: treadon org: np23 date: 2004-12-02 -! -! abstract: This routine sets derived constants -! -! program history log: -! 2004-12-02 treadon -! 2005-03-03 treadon - add implicit none -! -! input argument list: -! -! output argument list: -! -! attributes: -! language: f90 -! machine: ibm rs/6000 sp -! -!$$$ - implicit none - -! Trigonometric constants - pi = acos(-one) - deg2rad = pi/180.0_r_kind - rad2deg = one/deg2rad - cg_term = (sqrt(two*pi))/two ! constant for variational qc - tiny_r_kind = tiny(zero) - huge_r_kind = huge(zero) - tiny_single = tiny(zero_single) - huge_single = huge(zero_single) - huge_i_kind = huge(izero) - -! Geophysical parameters used in conversion of geopotential to -! geometric height - eccentricity_linear = sqrt(semi_major_axis**2 - semi_minor_axis**2) - eccentricity = eccentricity_linear / semi_major_axis - constants_initialized = .true. - - return - end subroutine init_constants_derived - - subroutine init_constants(regional) -!$$$ subprogram documentation block -! . . . . -! subprogram: init_constants set regional or global constants -! prgmmr: treadon org: np23 date: 2004-03-02 -! -! abstract: This routine sets constants specific to regional or global -! applications of the gsi -! -! program history log: -! 2004-03-02 treadon -! 2004-06-16 treadon, documentation -! 2004-10-28 treadon - use intrinsic TINY function to set value -! for smallest machine representable positive -! number -! 2004-12-03 treadon - move derived constants to init_constants_derived -! 2005-03-03 treadon - add implicit none -! -! input argument list: -! regional - if .true., set regional gsi constants; -! otherwise (.false.), use global constants -! -! output argument list: -! -! attributes: -! language: f90 -! machine: ibm rs/6000 sp -! -!$$$ - implicit none - logical regional - real(r_kind) reradius,g,r_d,r_v,cliq_wrf - - gamma = 0.0065 - -! Define regional constants here - if (regional) then - -! Name given to WRF constants - reradius = one/6370.e03_r_kind - g = 9.81_r_kind - r_d = 287.04_r_kind - r_v = 461.6_r_kind - cliq_wrf = 4190.0_r_kind - cp_mass = 1004.67_r_kind - -! Transfer WRF constants into unified GSI constants - rearth = one/reradius - grav = g - rd = r_d - rv = r_v - cv = cp-r_d - cliq = cliq_wrf - rd_over_cp_mass = rd / cp_mass - -! Define global constants here - else - rearth = 6.3712e+6_r_kind - grav = 9.80665e+0_r_kind - rd = 2.8705e+2_r_kind - rv = 4.6150e+2_r_kind - cv = 7.1760e+2_r_kind - cliq = 4.1855e+3_r_kind - cp_mass= zero - rd_over_cp_mass = zero - endif - - -! Now define derived constants which depend on constants -! which differ between global and regional applications. - -! Constants related to ozone assimilation - ozcon = grav*21.4e-9_r_kind - rozcon= one/ozcon - -! Constant used in vertical integral for precipitable water - tpwcon = 100.0_r_kind/grav - -! Derived atmospheric constants - fv = rv/rd-one ! used in virtual temperature equation - dldt = cvap-cliq - xa = -(dldt/rv) - xai = -(dldti/rv) - xb = xa+hvap/(rv*ttp) - xbi = xai+hsub/(rv*ttp) - eps = rd/rv - epsm1 = rd/rv-one - omeps = one-eps - factor1 = (cvap-cliq)/rv - factor2 = hvap/rv-factor1*t0c - cpr = cp*rd - el2orc = hvap*hvap/(rv*cp) - rd_over_g = rd/grav - rd_over_cp = rd/cp - g_over_rd = grav/rd - - return - end subroutine init_constants - -end module constants diff --git a/src/fv3nc2nemsio.fd/fv3_main.f90 b/src/fv3nc2nemsio.fd/fv3_main.f90 deleted file mode 100644 index 48c7440b..00000000 --- a/src/fv3nc2nemsio.fd/fv3_main.f90 +++ /dev/null @@ -1,215 +0,0 @@ -program fv3_main - use fv3_module - use netcdf - use nemsio_module - implicit none - - type(nemsio_gfile) :: gfile - type(nemsio_meta) :: meta_nemsio - integer,parameter :: nvar2d=48 - character(nemsio_charkind) :: name2d(nvar2d) - integer :: nvar3d - character(nemsio_charkind), allocatable :: name3din(:), name3dout(:) - character(nemsio_charkind) :: varname,levtype - character(len=300) :: inpath,outpath - character(len=100) :: infile2d,infile3d,outfile - character(len=10) :: analdate, cfhour - character(len=5) :: cfhr,cfhzh - character(len=2) :: nhcase - real , allocatable :: lons(:),lats(:),tmp2d(:,:), tmp2dx(:,:) - real*8,allocatable :: tmp1d(:),tmp1dx(:),fhours(:) - real*4 :: fhour - integer :: fhzh, nhcas - - integer :: ii,i,j,k,ncid2d,ncid3d,ifhr,nlevs,nlons,nlats,ntimes,nargs,iargc,YYYY,MM,DD,HH,stat,varid - - data name2d /'ALBDOsfc','CPRATsfc','PRATEsfc','DLWRFsfc','ULWRFsfc','DSWRFsfc','USWRFsfc','DSWRFtoa','USWRFtoa',& - 'ULWRFtoa','GFLUXsfc','HGTsfc','HPBLsfc',& - 'ICECsfc','SLMSKsfc','LHTFLsfc','SHTFLsfc','PRESsfc','PWATclm','SOILM','SOILW1','SOILW2','SOILW3','SOILW4','SPFH2m',& - 'SOILT1','SOILT2','SOILT3','SOILT4','TMP2m','TMPsfc','UGWDsfc','VGWDsfc','UFLXsfc','VFLXsfc','UGRD10m','VGRD10m',& - 'WEASDsfc','SNODsfc','ZORLsfc','VFRACsfc','F10Msfc','VTYPEsfc','STYPEsfc',& - 'TCDCclm', 'TCDChcl', 'TCDCmcl', 'TCDClcl'/ - - !===================================================================== - - ! read in from command line - nargs=iargc() - IF (nargs .NE. 10) THEN - print*,'usage fv3_interface analdate ifhr fhzh fhour inpath infile2d infile3d outpath,outfile,nhcase' - STOP 1 - ENDIF - call getarg(1,analdate) - call getarg(2,cfhr) - call getarg(3,cfhzh) - call getarg(4,cfhour) - call getarg(5,inpath) - call getarg(6,infile2d) - call getarg(7,infile3d) - call getarg(8,outpath) - call getarg(9,outfile) - call getarg(10,nhcase) -! print*,analdate,cfhr,cfhzh,cfhour,inpath,infile2d,infile3d,outpath,outfile,nhcase - - read(nhcase,'(i2.1)') nhcas - read(cfhr,'(i5.1)') ifhr - read(cfhzh,'(i5.1)') fhzh - read(cfhour,*) fhour - read(analdate(1:4),'(i4)') YYYY - read(analdate(5:6),'(i2)') MM - read(analdate(7:8),'(i2)') DD - read(analdate(9:10),'(i2)') HH - print*,"ifhr,fhzh,fhour,analdate ",ifhr,fhzh,fhour,analdate - - if (nhcas == 0 ) then !non-hydrostatic case - nvar3d=9 - allocate (name3din(nvar3d), name3dout(nvar3d)) - name3din=(/'ucomp ','vcomp ','temp ','sphum ','o3mr ','nhpres','w ','clwmr ','delp '/) - name3dout=(/'ugrd ','vgrd ','tmp ','spfh ','o3mr ','pres ','vvel ','clwmr','dpres'/) - else - nvar3d=8 - allocate (name3din(nvar3d), name3dout(nvar3d)) - name3din=(/'ucomp ','vcomp ','temp ','sphum ','o3mr ','hypres','clwmr ','delp '/) - name3dout=(/'ugrd ','vgrd ','tmp ','spfh ','o3mr ','pres ','clwmr','dpres'/) - endif - - ! open netcdf files - print*,'reading',trim(inpath)//'/'//trim(infile2d) - stat = nf90_open(trim(inpath)//'/'//trim(infile2d),NF90_NOWRITE, ncid2d) - if (stat .NE.0) print*,stat - print*,'reading',trim(inpath)//'/'//trim(infile3d) - stat = nf90_open(trim(inpath)//'/'//trim(infile3d),NF90_NOWRITE, ncid3d) - if (stat .NE.0) print*,stat - ! get dimesions - - stat = nf90_inq_dimid(ncid2d,'time',varid) - if (stat .NE.0) print*,stat,varid - if (stat .NE. 0) STOP 1 - stat = nf90_inquire_dimension(ncid2d,varid,len=ntimes) - if (stat .NE.0) print*,stat,ntimes - if (stat .NE. 0) STOP 1 - allocate(fhours(ntimes)) - stat = nf90_inq_varid(ncid2d,'time',varid) - if (stat .NE. 0) STOP 1 - stat = nf90_get_var(ncid2d,varid,fhours) - if (stat .NE.0) print*,stat,fhours - if (stat .NE. 0) STOP 1 - - stat = nf90_inq_dimid(ncid3d,'grid_xt',varid) - if (stat .NE.0) print*,stat,varid - if (stat .NE. 0) STOP 1 - stat = nf90_inquire_dimension(ncid3d,varid,len=nlons) - if (stat .NE.0) print*,stat,nlons - if (stat .NE. 0) STOP 1 - allocate(lons(nlons)) - allocate(tmp1d(nlons)) - stat = nf90_inq_varid(ncid3d,'grid_xt',varid) - if (stat .NE. 0) STOP 1 - stat = nf90_get_var(ncid3d,varid,tmp1d) - if (stat .NE.0) print*,stat - if (stat .NE. 0) STOP 1 - - lons=real(tmp1d,kind=4) - !print*,lons(1),lons(3072) - deallocate(tmp1d) - - stat = nf90_inq_dimid(ncid3d,'grid_yt',varid) - if (stat .NE.0) print*,stat - if (stat .NE. 0) STOP 1 - stat = nf90_inquire_dimension(ncid3d,varid,len=nlats) - if (stat .NE.0) print*,stat - if (stat .NE. 0) STOP 1 - allocate(lats(nlats)) - allocate(tmp1d(nlats)) - allocate(tmp1dx(nlats)) - stat = nf90_inq_varid(ncid3d,'grid_yt',varid) - stat = nf90_get_var(ncid3d,varid,tmp1dx,start=(/1/),count=(/nlats/)) - if (stat .NE.0) print*,stat - if (stat .NE. 0) STOP 1 - do j=1,nlats - tmp1d(j)=tmp1dx(nlats-j+1) - enddo - lats=real(tmp1d,kind=4) - print*,"lats_beg, lats_end",lats(1),lats(nlats) - deallocate(tmp1d, tmp1dx) - - stat = nf90_inq_dimid(ncid3d,'pfull',varid) - if (stat .NE.0) print*,stat - if (stat .NE. 0) STOP 1 - stat = nf90_inquire_dimension(ncid3d,varid,len=nlevs) - if (stat .NE.0) print*,stat - if (stat .NE. 0) STOP 1 - - call define_nemsio_meta(meta_nemsio,nlons,nlats,nlevs,nvar2d,nvar3d,lons,lats) - - allocate (tmp2d(nlons,nlats)) - allocate (tmp2dx(nlons,nlats)) - - meta_nemsio%idate(1)=YYYY - meta_nemsio%idate(2)=MM - meta_nemsio%idate(3)=DD - meta_nemsio%idate(4)=HH - - meta_nemsio%varrval(1)=float(fhzh) -! if (ifhr.EQ.0) then -! meta_nemsio%varrval(1)=0.0 -! else -! meta_nemsio%varrval(1)=(ifhr-1.0)*6.0 -! endif - - ! read in data - meta_nemsio%nfhour= fhours(ifhr) - meta_nemsio%fhour= fhours(ifhr) - print*,fhours(ifhr),ifhr,'calling netcdf read' -!--for ifhr=1, fhours=dt but fhour=00 if diag is determined by FHOUT - if (fhour .ne. fhours(ifhr) .and. ifhr.gt.1 )then - print*, 'requested ',fhour, ' not equal to fhours(ifhr) ', fhours(ifhr) - print*, 'abort ! ' - stop 1 - endif - - call nems_write_init(outpath,outfile,meta_nemsio,gfile) -! read in all of the 2d variables and write out - print*,'calling write',meta_nemsio%rlat_min,meta_nemsio%rlat_max - print*,'lats',minval(meta_nemsio%lat),maxval(meta_nemsio%lat) - print *,'loop over 2d variables' - DO i=1,nvar2d - print *,i,trim(name2d(i)) - call fv3_netcdf_read_2d(ncid2d,ifhr,meta_nemsio,name2d(i),tmp2dx) - do ii=1,nlons - do j=1,nlats - tmp2d(ii,j)=tmp2dx(ii,nlats-j+1) - enddo - enddo - call nems_write(gfile,meta_nemsio%recname(i),meta_nemsio%reclevtyp(i),meta_nemsio%reclev(i), & - nlons*nlats,tmp2d,stat) - ENDDO - levtype='mid layer' -! loop through 3d fields - print *,'loop over 3d variables' - DO i=1,nvar3d - print*,i,trim(name3din(i)) - DO k=1,nlevs -! print*,k - call fv3_netcdf_read_3d(ncid3d,ifhr,meta_nemsio,name3din(i),k,tmp2dx) - do ii=1,nlons - do j=1,nlats - tmp2d(ii,j)=tmp2dx(ii,nlats-j+1) - enddo - enddo - call nems_write(gfile,name3dout(i),levtype,nlevs-k+1,nlons*nlats,tmp2d(:,:),stat) - IF (stat .NE. 0) then - print*,'error writing ,named3dout(i)',stat - STOP 1 - ENDIF - ENDDO - ENDDO - - call nemsio_close(gfile,iret=stat) - stat = nf90_close(ncid2d) - stat = nf90_close(ncid3d) - - deallocate(tmp2dx,tmp2d) - deallocate(name3din,name3dout) - - stop -end program fv3_main diff --git a/src/fv3nc2nemsio.fd/fv3_module.f90 b/src/fv3nc2nemsio.fd/fv3_module.f90 deleted file mode 100644 index 8d161acf..00000000 --- a/src/fv3nc2nemsio.fd/fv3_module.f90 +++ /dev/null @@ -1,372 +0,0 @@ -module fv3_module - - - !======================================================================= - - ! Define associated modules and subroutines - - !----------------------------------------------------------------------- - use netcdf - use constants - use kinds - use nemsio_module - - type nemsio_meta - character(nemsio_charkind), dimension(:), allocatable :: recname - character(nemsio_charkind), dimension(:), allocatable :: reclevtyp - character(16), dimension(:), allocatable :: variname - character(16), dimension(:), allocatable :: varrname - character(16), dimension(:), allocatable :: varr8name - character(16), dimension(:), allocatable :: aryiname - character(16), dimension(:), allocatable :: aryr8name - character(nemsio_charkind8) :: gdatatype - character(nemsio_charkind8) :: modelname - real(nemsio_realkind) :: rlon_min - real(nemsio_realkind) :: rlon_max - real(nemsio_realkind) :: rlat_min - real(nemsio_realkind) :: rlat_max - real(nemsio_realkind), dimension(:), allocatable :: lon - real(nemsio_realkind), dimension(:), allocatable :: lat - real(nemsio_realkind), dimension(:), allocatable :: varrval - integer(nemsio_intkind), dimension(:,:), allocatable :: aryival - integer(nemsio_intkind), dimension(:), allocatable :: reclev - integer(nemsio_intkind), dimension(:), allocatable :: varival - integer(nemsio_intkind), dimension(:), allocatable :: aryilen - integer(nemsio_intkind), dimension(:), allocatable :: aryr8len - integer(nemsio_intkind) :: idate(7) - integer(nemsio_intkind) :: version - integer(nemsio_intkind) :: nreo_vc - integer(nemsio_intkind) :: nrec - integer(nemsio_intkind) :: nmeta - integer(nemsio_intkind) :: nmetavari - integer(nemsio_intkind) :: nmetaaryi - integer(nemsio_intkind) :: nmetavarr - integer(nemsio_intkind) :: nfhour - integer(nemsio_intkind) :: nfminute - integer(nemsio_intkind) :: nfsecondn - integer(nemsio_intkind) :: nfsecondd - integer(nemsio_intkind) :: dimx - integer(nemsio_intkind) :: dimy - integer(nemsio_intkind) :: dimz - integer(nemsio_intkind) :: nframe - integer(nemsio_intkind) :: nsoil - integer(nemsio_intkind) :: ntrac - integer(nemsio_intkind) :: ncldt - integer(nemsio_intkind) :: idvc - integer(nemsio_intkind) :: idsl - integer(nemsio_intkind) :: idvm - integer(nemsio_intkind) :: idrt - integer(nemsio_intkind) :: fhour - - end type nemsio_meta ! type nemsio_meta - contains -!----------------------------------------------------------------------- - subroutine fv3_netcdf_read_2d(ncid2d,ifhr,meta_nemsio,varname,data2d) - - implicit none - type(nemsio_meta) :: meta_nemsio - integer :: ncid2d - integer :: ifhr,varid,stat - real :: data2d(meta_nemsio%dimx,meta_nemsio%dimy) - character(nemsio_charkind) :: varname - - ! loop through 2d data - stat = nf90_inq_varid(ncid2d,trim(varname),varid) - !print*,stat,varid,trim(varname) - stat = nf90_get_var(ncid2d,varid,data2d,start=(/1,1,ifhr/),count=(/meta_nemsio%dimx,meta_nemsio%dimy,1/)) - IF (stat .NE. 0 ) THEN - print*,'error reading ',varname - STOP - ENDIF - -end subroutine fv3_netcdf_read_2d -!----------------------------------------------------------------------- - - subroutine fv3_netcdf_read_3d(ncid3d,ifhr,meta_nemsio,varname,k,data2d) - - implicit none - - type(nemsio_meta) :: meta_nemsio - integer :: ncid3d - integer :: k - integer :: ifhr,varid,stat - character(nemsio_charkind) :: varname - !real :: data3d(meta_nemsio%dimx,meta_nemsio%dimy,meta_nemsio%dimz) - real :: data2d(meta_nemsio%dimx,meta_nemsio%dimy) - - - stat = nf90_inq_varid(ncid3d,trim(varname),varid) - !print*,stat,varname,varid - !stat = nf90_get_var(ncid3d,varid,data3d,start=(/1,1,1,ifhr/),count=(/meta_nemsio%dimx,meta_nemsio%dimy,meta_nemsio%dimz,1/)) - stat = nf90_get_var(ncid3d,varid,data2d,start=(/1,1,k,ifhr/),count=(/meta_nemsio%dimx,meta_nemsio%dimy,1,1/)) - - IF (stat .NE. 0 ) THEN - print*,'error reading ',varname - STOP - ENDIF - -end subroutine fv3_netcdf_read_3d -!----------------------------------------------------------------------- - - subroutine define_nemsio_meta(meta_nemsio,nlons,nlats,nlevs,nvar2d,nvar3d,lons,lats) - implicit none - type(nemsio_meta) :: meta_nemsio - integer :: nlons,nlats,nlevs,i,j,k,nvar2d,nvar3d - integer*8 :: ct - real :: lons(nlons),lats(nlats) -! local - - meta_nemsio%idate(1:6) = 0 - meta_nemsio%idate(7) = 1 - meta_nemsio%modelname = 'GFS' - meta_nemsio%version = 198410 - meta_nemsio%nrec = nvar2d + nlevs*nvar3d - meta_nemsio%nmeta = 8 - meta_nemsio%nmetavari = 3 - meta_nemsio%nmetavarr = 1 - meta_nemsio%nmetaaryi = 1 - meta_nemsio%dimx = nlons - meta_nemsio%dimy = nlats - meta_nemsio%dimz = nlevs - meta_nemsio%rlon_min = minval(lons) - meta_nemsio%rlon_max = maxval(lons) - meta_nemsio%rlat_min = minval(lats) - meta_nemsio%rlat_max = maxval(lats) - meta_nemsio%nsoil = 4 - meta_nemsio%nframe = 0 - meta_nemsio%nfminute = 0 - meta_nemsio%nfsecondn = 0 - meta_nemsio%nfsecondd = 1 - meta_nemsio%ntrac = 3 - meta_nemsio%idrt = 0 - meta_nemsio%ncldt = 3 - meta_nemsio%idvc = 2 - - - allocate(meta_nemsio%recname(meta_nemsio%nrec)) - allocate(meta_nemsio%reclevtyp(meta_nemsio%nrec)) - allocate(meta_nemsio%reclev(meta_nemsio%nrec)) - allocate(meta_nemsio%variname(meta_nemsio%nmetavari)) - allocate(meta_nemsio%varival(meta_nemsio%nmetavari)) - allocate(meta_nemsio%aryiname(meta_nemsio%nmetavari)) - allocate(meta_nemsio%aryilen(meta_nemsio%nmetavari)) - allocate(meta_nemsio%varrname(meta_nemsio%nmetavarr)) - allocate(meta_nemsio%varrval(meta_nemsio%nmetavarr)) - allocate(meta_nemsio%lon(nlons*nlats)) - allocate(meta_nemsio%lat(nlons*nlats)) - - meta_nemsio%varrname(1)='zhour' - meta_nemsio%variname(1)='cu_physics' - meta_nemsio%varival(1)=4 - meta_nemsio%variname(2)='mp_physics' - meta_nemsio%varival(2)=1000 - meta_nemsio%variname(3)='IVEGSRC' - meta_nemsio%varival(3)=2 - ct=1 - DO j=1,nlats - DO i=1,nlons - meta_nemsio%lon(ct) = lons(i) - meta_nemsio%lat(ct) = lats(j) - ct=ct+1 - ENDDO - ENDDO - - meta_nemsio%aryilen(1) = nlats/2 - meta_nemsio%aryiname(1) = 'lpl' - meta_nemsio%reclev(:)=1 - meta_nemsio%recname(1) = 'albdo_ave' - meta_nemsio%reclevtyp(1) = 'sfc' - meta_nemsio%recname(2) = 'cprat_ave' - meta_nemsio%reclevtyp(2) = 'sfc' - meta_nemsio%recname(3) = 'prate_ave' - meta_nemsio%reclevtyp(3) = 'sfc' - meta_nemsio%recname(4) = 'dlwrf_ave' - meta_nemsio%reclevtyp(4) = 'sfc' - meta_nemsio%recname(5) = 'ulwrf_ave' - meta_nemsio%reclevtyp(5) = 'sfc' - meta_nemsio%recname(6) = 'dswrf_ave' - meta_nemsio%reclevtyp(6) = 'sfc' - meta_nemsio%recname(7) = 'uswrf_ave' - meta_nemsio%reclevtyp(7) = 'sfc' - meta_nemsio%recname(8) = 'dswrf_ave' - meta_nemsio%reclevtyp(8) = 'nom. top' - meta_nemsio%recname(9) = 'uswrf_ave' - meta_nemsio%reclevtyp(9) = 'nom. top' - meta_nemsio%recname(10) = 'ulwrf_ave' - meta_nemsio%reclevtyp(10) = 'nom. top' - meta_nemsio%recname(11) = 'gflux_ave' - meta_nemsio%reclevtyp(11) = 'sfc' - meta_nemsio%recname(12) = 'hgt' - meta_nemsio%reclevtyp(12) = 'sfc' - meta_nemsio%recname(13) = 'hpbl' - meta_nemsio%reclevtyp(13) = 'sfc' - meta_nemsio%recname(14) = 'icec' - meta_nemsio%reclevtyp(14) = 'sfc' - meta_nemsio%recname(15) = 'land' - meta_nemsio%reclevtyp(15) = 'sfc' - meta_nemsio%recname(16) = 'lhtfl_ave' - meta_nemsio%reclevtyp(16) = 'sfc' - meta_nemsio%recname(17) = 'shtfl_ave' - meta_nemsio%reclevtyp(17) = 'sfc' - meta_nemsio%recname(18) = 'pres' - meta_nemsio%reclevtyp(18) = 'sfc' - meta_nemsio%recname(19) = 'pwat' - meta_nemsio%reclevtyp(19) = 'atmos col' - meta_nemsio%recname(20) = 'soilm' - meta_nemsio%reclevtyp(20) = '0-200 cm down' - meta_nemsio%recname(21) = 'soilw' - meta_nemsio%reclevtyp(21) = '0-10 cm down' - meta_nemsio%recname(22) = 'soilw' - meta_nemsio%reclevtyp(22) = '10-40 cm down' - meta_nemsio%recname(23) = 'soilw' - meta_nemsio%reclevtyp(23) = '40-100 cm down' - meta_nemsio%recname(24) = 'soilw' - meta_nemsio%reclevtyp(24) = '100-200 cm down' - meta_nemsio%recname(25) = 'spfh' - meta_nemsio%reclevtyp(25) = '2 m above gnd' - meta_nemsio%recname(26) = 'tmp' - meta_nemsio%reclevtyp(26) = '0-10 cm down' - meta_nemsio%recname(27) = 'tmp' - meta_nemsio%reclevtyp(27) = '10-40 cm down' - meta_nemsio%recname(28) = 'tmp' - meta_nemsio%reclevtyp(28) = '40-100 cm down' - meta_nemsio%recname(29) = 'tmp' - meta_nemsio%reclevtyp(29) = '100-200 cm down' - meta_nemsio%recname(30) = 'tmp' - meta_nemsio%reclevtyp(30) = '2 m above gnd' - meta_nemsio%recname(31) = 'tmp' - meta_nemsio%reclevtyp(31) = 'sfc' - meta_nemsio%recname(32) = 'ugwd' - meta_nemsio%reclevtyp(32) = 'sfc' - meta_nemsio%recname(33) = 'vgwd' - meta_nemsio%reclevtyp(33) = 'sfc' - meta_nemsio%recname(34) = 'uflx_ave' - meta_nemsio%reclevtyp(34) = 'sfc' - meta_nemsio%recname(35) = 'vflx_ave' - meta_nemsio%reclevtyp(35) = 'sfc' - meta_nemsio%recname(36) = 'ugrd' - meta_nemsio%reclevtyp(36) = '10 m above gnd' - meta_nemsio%recname(37) = 'vgrd' - meta_nemsio%reclevtyp(37) = '10 m above gnd' - meta_nemsio%recname(38) = 'weasd' - meta_nemsio%reclevtyp(38) = 'sfc' - meta_nemsio%recname(39) = 'snod' - meta_nemsio%reclevtyp(39) = 'sfc' - meta_nemsio%recname(40) = 'zorl' - meta_nemsio%reclevtyp(40) = 'sfc' - meta_nemsio%recname(41) = 'vfrac' - meta_nemsio%reclevtyp(41) = 'sfc' - meta_nemsio%recname(42) = 'f10m' - meta_nemsio%reclevtyp(42) = 'sfc' - meta_nemsio%recname(43) = 'vtype' - meta_nemsio%reclevtyp(43) = 'sfc' - meta_nemsio%recname(44) = 'stype' - meta_nemsio%reclevtyp(44) = 'sfc' - meta_nemsio%recname(45) = 'tcdc_ave' - meta_nemsio%reclevtyp(45) = 'atmos col' - meta_nemsio%recname(46) = 'tcdc_ave' - meta_nemsio%reclevtyp(46) = 'high cld lay' - meta_nemsio%recname(47) = 'tcdc_ave' - meta_nemsio%reclevtyp(47) = 'mid cld lay' - meta_nemsio%recname(48) = 'tcdc_ave' - meta_nemsio%reclevtyp(48) = 'low cld lay' -! loop through 3d variables - DO k = 1, nlevs - meta_nemsio%recname(k+nvar2d) = 'ugrd' - meta_nemsio%reclevtyp(k+nvar2d) = 'mid layer' - meta_nemsio%reclev(k+nvar2d) = k - meta_nemsio%recname(k+nvar2d+nlevs) = 'vgrd' - meta_nemsio%reclevtyp(k+nvar2d+nlevs) = 'mid layer' - meta_nemsio%reclev(k+nvar2d+nlevs) = k - meta_nemsio%recname(k+nvar2d+nlevs*2) = 'tmp' - meta_nemsio%reclevtyp(k+nvar2d+nlevs*2) = 'mid layer' - meta_nemsio%reclev(k+nvar2d+nlevs*2) = k - meta_nemsio%recname(k+nvar2d+nlevs*3) = 'spfh' - meta_nemsio%reclevtyp(k+nvar2d+nlevs*3) = 'mid layer' - meta_nemsio%reclev(k+nvar2d+nlevs*3) = k - meta_nemsio%recname(k+nvar2d+nlevs*4) = 'o3mr' - meta_nemsio%reclevtyp(k+nvar2d+nlevs*4) = 'mid layer' - meta_nemsio%reclev(k+nvar2d+nlevs*4) = k - meta_nemsio%recname(k+nvar2d+nlevs*5) = 'pres' - meta_nemsio%reclevtyp(k+nvar2d+nlevs*5) = 'mid layer' - meta_nemsio%reclev(k+nvar2d+nlevs*5) = k - meta_nemsio%recname(k+nvar2d+nlevs*6) = 'clwmr' - meta_nemsio%reclevtyp(k+nvar2d+nlevs*6) = 'mid layer' - meta_nemsio%reclev(k+nvar2d+nlevs*6) = k - meta_nemsio%recname(k+nvar2d+nlevs*7) = 'dpres' - meta_nemsio%reclevtyp(k+nvar2d+nlevs*7) = 'mid layer' - meta_nemsio%reclev(k+nvar2d+nlevs*7) = k - if (nvar3d == 9) then - meta_nemsio%recname(k+nvar2d+nlevs*8) = 'vvel' - meta_nemsio%reclevtyp(k+nvar2d+nlevs*8) = 'mid layer' - meta_nemsio%reclev(k+nvar2d+nlevs*8) = k - endif - ENDDO - - end subroutine define_nemsio_meta - - subroutine nems_write_init(datapath,filename_base,meta_nemsio,gfile) - - - implicit none - - type(nemsio_meta) :: meta_nemsio - character(len=200) :: datapath - character(len=100) :: filename_base - character(len=400) :: filename - type(nemsio_gfile) :: gfile - integer :: nemsio_iret - integer :: i, j, k - - write(filename,500) trim(datapath)//'/'//trim(filename_base) -500 format(a,i3.3) - print*,trim(filename) - call nemsio_init(iret=nemsio_iret) - print*,'iret=',nemsio_iret - !gfile%gtype = 'NEMSIO' - meta_nemsio%gdatatype = 'bin4' - call nemsio_open(gfile,trim(filename),'write', & - & iret=nemsio_iret, & - & modelname=trim(meta_nemsio%modelname), & - & version=meta_nemsio%version,gdatatype=meta_nemsio%gdatatype, & - & dimx=meta_nemsio%dimx,dimy=meta_nemsio%dimy, & - & dimz=meta_nemsio%dimz,rlon_min=meta_nemsio%rlon_min, & - & rlon_max=meta_nemsio%rlon_max,rlat_min=meta_nemsio%rlat_min, & - & rlat_max=meta_nemsio%rlat_max, & - & lon=meta_nemsio%lon,lat=meta_nemsio%lat, & - & idate=meta_nemsio%idate,nrec=meta_nemsio%nrec, & - & nframe=meta_nemsio%nframe,idrt=meta_nemsio%idrt,ncldt= & - & meta_nemsio%ncldt,idvc=meta_nemsio%idvc, & - & nfhour=meta_nemsio%nfhour,nfminute=meta_nemsio%nfminute, & - & nfsecondn=meta_nemsio%nfsecondn,nmeta=meta_nemsio%nmeta, & - & nfsecondd=meta_nemsio%nfsecondd,extrameta=.true., & - & nmetaaryi=meta_nemsio%nmetaaryi,recname=meta_nemsio%recname, & - & nmetavari=meta_nemsio%nmetavari,variname=meta_nemsio%variname, & - & varival=meta_nemsio%varival,varrval=meta_nemsio%varrval, & - & nmetavarr=meta_nemsio%nmetavarr,varrname=meta_nemsio%varrname, & - & reclevtyp=meta_nemsio%reclevtyp, & - & reclev=meta_nemsio%reclev,aryiname=meta_nemsio%aryiname, & - & aryilen=meta_nemsio%aryilen) - print*,'iret=',nemsio_iret - end subroutine nems_write_init - - -!------------------------------------------------------ - subroutine nems_write(gfile,recname,reclevtyp,level,dimx,data2d,iret) - - implicit none - type(nemsio_gfile) :: gfile - integer :: iret,level,dimx - real :: data2d(dimx) - character(nemsio_charkind) :: recname, reclevtyp - - call nemsio_writerecv(gfile,recname,levtyp=reclevtyp,lev=level,data=data2d,iret=iret) - if (iret.NE.0) then - print*,'error writing',recname,level,iret - STOP - ENDIF - - end subroutine nems_write - - -end module fv3_module diff --git a/src/fv3nc2nemsio.fd/kinds.f90 b/src/fv3nc2nemsio.fd/kinds.f90 deleted file mode 100644 index b3378bfc..00000000 --- a/src/fv3nc2nemsio.fd/kinds.f90 +++ /dev/null @@ -1,107 +0,0 @@ -! this module was extracted from the GSI version operational -! at NCEP in Dec. 2007. -module kinds -!$$$ module documentation block -! . . . . -! module: kinds -! prgmmr: treadon org: np23 date: 2004-08-15 -! -! abstract: Module to hold specification kinds for variable declaration. -! This module is based on (copied from) Paul vanDelst's -! type_kinds module found in the community radiative transfer -! model -! -! module history log: -! 2004-08-15 treadon -! -! Subroutines Included: -! -! Functions Included: -! -! remarks: -! The numerical data types defined in this module are: -! i_byte - specification kind for byte (1-byte) integer variable -! i_short - specification kind for short (2-byte) integer variable -! i_long - specification kind for long (4-byte) integer variable -! i_llong - specification kind for double long (8-byte) integer variable -! r_single - specification kind for single precision (4-byte) real variable -! r_double - specification kind for double precision (8-byte) real variable -! r_quad - specification kind for quad precision (16-byte) real variable -! -! i_kind - generic specification kind for default integer -! r_kind - generic specification kind for default floating point -! -! -! attributes: -! language: f90 -! machine: ibm RS/6000 SP -! -!$$$ end documentation block - implicit none - private - -! Integer type definitions below - -! Integer types - integer, parameter, public :: i_byte = selected_int_kind(1) ! byte integer - integer, parameter, public :: i_short = selected_int_kind(4) ! short integer - integer, parameter, public :: i_long = selected_int_kind(8) ! long integer - integer, parameter, private :: llong_t = selected_int_kind(16) ! llong integer - integer, parameter, public :: i_llong = max( llong_t, i_long ) - -! Expected 8-bit byte sizes of the integer kinds - integer, parameter, public :: num_bytes_for_i_byte = 1 - integer, parameter, public :: num_bytes_for_i_short = 2 - integer, parameter, public :: num_bytes_for_i_long = 4 - integer, parameter, public :: num_bytes_for_i_llong = 8 - -! Define arrays for default definition - integer, parameter, private :: num_i_kinds = 4 - integer, parameter, dimension( num_i_kinds ), private :: integer_types = (/ & - i_byte, i_short, i_long, i_llong /) - integer, parameter, dimension( num_i_kinds ), private :: integer_byte_sizes = (/ & - num_bytes_for_i_byte, num_bytes_for_i_short, & - num_bytes_for_i_long, num_bytes_for_i_llong /) - -! Default values -! **** CHANGE THE FOLLOWING TO CHANGE THE DEFAULT INTEGER TYPE KIND *** - integer, parameter, private :: default_integer = 3 ! 1=byte, - ! 2=short, - ! 3=long, - ! 4=llong - integer, parameter, public :: i_kind = integer_types( default_integer ) - integer, parameter, public :: num_bytes_for_i_kind = & - integer_byte_sizes( default_integer ) - - -! Real definitions below - -! Real types - integer, parameter, public :: r_single = selected_real_kind(6) ! single precision - integer, parameter, public :: r_double = selected_real_kind(15) ! double precision - integer, parameter, private :: quad_t = selected_real_kind(20) ! quad precision - integer, parameter, public :: r_quad = max( quad_t, r_double ) - -! Expected 8-bit byte sizes of the real kinds - integer, parameter, public :: num_bytes_for_r_single = 4 - integer, parameter, public :: num_bytes_for_r_double = 8 - integer, parameter, public :: num_bytes_for_r_quad = 16 - -! Define arrays for default definition - integer, parameter, private :: num_r_kinds = 3 - integer, parameter, dimension( num_r_kinds ), private :: real_kinds = (/ & - r_single, r_double, r_quad /) - integer, parameter, dimension( num_r_kinds ), private :: real_byte_sizes = (/ & - num_bytes_for_r_single, num_bytes_for_r_double, & - num_bytes_for_r_quad /) - -! Default values -! **** CHANGE THE FOLLOWING TO CHANGE THE DEFAULT REAL TYPE KIND *** - integer, parameter, private :: default_real = 1 ! 1=single, - ! 2=double, - ! 3=quad - integer, parameter, public :: r_kind = real_kinds( default_real ) - integer, parameter, public :: num_bytes_for_r_kind = & - real_byte_sizes( default_real ) - -end module kinds diff --git a/src/gaussian_sfcanl.fd/CMakeLists.txt b/src/gaussian_sfcanl.fd/CMakeLists.txt index 6bfe8387..7e37fab0 100644 --- a/src/gaussian_sfcanl.fd/CMakeLists.txt +++ b/src/gaussian_sfcanl.fd/CMakeLists.txt @@ -13,7 +13,6 @@ add_executable(${exe_name} ${fortran_src}) target_link_libraries(${exe_name} PRIVATE NetCDF::NetCDF_Fortran bacio::bacio_4 sp::sp_4 - w3emc::w3emc_d - nemsio::nemsio) + w3emc::w3emc_d) install(TARGETS ${exe_name} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/src/gaussian_sfcanl.fd/gaussian_sfcanl.f90 b/src/gaussian_sfcanl.fd/gaussian_sfcanl.f90 index acce575c..97de504e 100644 --- a/src/gaussian_sfcanl.fd/gaussian_sfcanl.f90 +++ b/src/gaussian_sfcanl.fd/gaussian_sfcanl.f90 @@ -1,11 +1,9 @@ !------------------------------------------------------------------ -! +! ! Read in surface and nst data on the cubed-sphere grid, ! interpolate it to the gaussian grid, and output the result -! to a nemsio or netcdf file. To not process nst data, -! set flag 'donst' to 'no'. To process nst, set to 'yes'. -! To output gaussian file in netcdf, set netcdf_out=.true. -! Otherwise, nemsio format will be output. +! to a netcdf file. To not process nst data, +! set flag 'donst' to '.false.'. To process nst, set to '.true.'. ! ! Input files: ! ------------ @@ -18,37 +16,33 @@ ! ! Output files: ! ------------- -! sfc.gaussian.analysis.file surface data on gaussian grid - -! nemsio or netcdf. +! sfc.gaussian.analysis.file surface data on gaussian grid - netcdf. ! ! Namelist variables: ! ------------------- ! yy/mm/dd/hh year/month/day/hour of data. ! i/jgaus i/j dimension of gaussian grid. -! donst When 'no' do not process nst data. -! When 'yes' process nst data. -! netcdf_out When 'true', output gaussian file in -! netcdf. Otherwise output nemsio format. +! donst When '.true.' process nst data, default; .false. ! ! 2018-Jan-30 Gayno Initial version ! 2019-Oct-30 Gayno Option to output gaussian analysis file ! in netcdf. +! 2023-Mar-04 Mahajan Deprecate option to output gaussian analysis file +! in nemsio. ! !------------------------------------------------------------------ module io - use nemsio_module - implicit none - character(len=3) :: donst + logical :: donst integer, parameter :: num_tiles = 6 integer :: itile, jtile, igaus, jgaus - integer(nemsio_intkind) :: idate(8) + integer :: idate(8) type :: sfc_data ! surface variables @@ -129,19 +123,15 @@ program main integer :: yy, mm, dd, hh integer, allocatable :: col(:), row(:) - logical :: netcdf_out - real(kind=8), allocatable :: s(:) - namelist /setup/ yy, mm, dd, hh, igaus, jgaus, donst, netcdf_out + namelist /setup/ yy, mm, dd, hh, igaus, jgaus, donst call w3tagb('GAUSSIAN_SFCANL',2018,0179,0055,'NP20') print*,"- BEGIN EXECUTION" - netcdf_out = .true. - - donst = 'no' + donst = .false. print* print*,"- READ SETUP NAMELIST" @@ -242,7 +232,7 @@ program main allocate(gaussian_data%smc(igaus*jgaus,4)) allocate(gaussian_data%stc(igaus*jgaus,4)) - if (trim(donst) == "yes" .or. trim(donst) == "YES") then + if (donst) then allocate(gaussian_data%c0(igaus*jgaus)) ! nst allocate(gaussian_data%cd(igaus*jgaus)) allocate(gaussian_data%dconv(igaus*jgaus)) @@ -295,7 +285,7 @@ program main gaussian_data%hice(row(i)) = gaussian_data%hice(row(i)) + s(i)*tile_data%hice(col(i)) gaussian_data%snoalb(row(i)) = gaussian_data%snoalb(row(i)) + s(i)*tile_data%snoalb(col(i)) gaussian_data%srflag(row(i)) = gaussian_data%srflag(row(i)) + s(i)*tile_data%srflag(col(i)) - if (trim(donst) == "yes" .or. trim(donst) == "YES") then + if (donst) then gaussian_data%c0(row(i)) = gaussian_data%c0(row(i)) + s(i)*tile_data%c0(col(i)) gaussian_data%cd(row(i)) = gaussian_data%cd(row(i)) + s(i)*tile_data%cd(col(i)) gaussian_data%dconv(row(i)) = gaussian_data%dconv(row(i)) + s(i)*tile_data%dconv(col(i)) @@ -358,7 +348,7 @@ program main deallocate(tile_data%smc) deallocate(tile_data%stc) - if (trim(donst) == "yes" .or. trim(donst) == "YES") then + if (donst) then deallocate(tile_data%c0) deallocate(tile_data%cd) deallocate(tile_data%dconv) @@ -378,14 +368,10 @@ program main endif !------------------------------------------------------------------------------ -! Write gaussian data to either netcdf or nemsio file. +! Write gaussian data to netcdf file. !------------------------------------------------------------------------------ - if (netcdf_out) then - call write_sfc_data_netcdf - else - call write_sfc_data_nemsio - endif + call write_sfc_data_netcdf deallocate(gaussian_data%orog) deallocate(gaussian_data%t2m) @@ -423,7 +409,7 @@ program main deallocate(gaussian_data%smc) deallocate(gaussian_data%stc) - if (trim(donst) == "yes" .or. trim(donst) == "YES") then + if (donst) then deallocate(gaussian_data%c0) deallocate(gaussian_data%cd) deallocate(gaussian_data%dconv) @@ -802,7 +788,7 @@ subroutine write_sfc_data_netcdf ! Determine what variables to output (noah, or noah plus nst). !------------------------------------------------------------------------------------------- - if (trim(donst) == "yes" .or. trim(donst) == "YES") then + if (donst) then num_vars = num_noah + num_nst else num_vars = num_noah @@ -817,7 +803,7 @@ subroutine write_sfc_data_netcdf name(1:num_noah) = noah_name units(1:num_noah) = noah_units - if (trim(donst) == "yes" .or. trim(donst) == "YES") then + if (donst) then do n = 1, num_nst var(n+num_noah) = nst_var(n) name(n+num_noah) = nst_name(n) @@ -1056,521 +1042,6 @@ subroutine get_netcdf_var(var, dummy) end subroutine get_netcdf_var -!------------------------------------------------------------------------------------------- -! Write gaussian surface data to nemsio file. -!------------------------------------------------------------------------------------------- - - subroutine write_sfc_data_nemsio - - use nemsio_module - use io - - implicit none - - integer(nemsio_intkind), parameter :: nrec_all=60 - integer(nemsio_intkind), parameter :: nmetaaryi=1 - integer(nemsio_intkind), parameter :: nmetavari=4 - integer(nemsio_intkind), parameter :: nmetavarr=1 - integer(nemsio_intkind), parameter :: nmetavarc=2 - - character(nemsio_charkind) :: recname_all(nrec_all) - character(nemsio_charkind) :: reclevtyp_all(nrec_all) - character(nemsio_charkind) :: aryiname(nmetaaryi) - character(nemsio_charkind) :: variname(nmetavari) - character(nemsio_charkind) :: varrname(nmetavarr) - character(nemsio_charkind) :: varcname(nmetavarc) - character(nemsio_charkind) :: varcval(nmetavarc) - character(nemsio_charkind), allocatable :: recname(:) - character(nemsio_charkind), allocatable :: reclevtyp(:) - - integer(nemsio_intkind) :: iret, version, nrec - integer(nemsio_intkind) :: reclev_all(nrec_all) - integer(nemsio_intkind) :: aryival(jgaus,nmetaaryi) - integer(nemsio_intkind) :: aryilen(nmetaaryi) - integer(nemsio_intkind) :: varival(nmetavari) - integer :: i, k, n, nvcoord, levs_vcoord - integer(nemsio_intkind), allocatable :: reclev(:) - - real(nemsio_realkind), allocatable :: the_data(:) - real(nemsio_realkind) :: varrval(nmetavarr) - real(nemsio_realkind), allocatable :: lat(:), lon(:) - real(kind=4), allocatable :: dummy(:,:), slat(:), wlat(:) - real(nemsio_realkind), allocatable :: vcoord(:,:,:) - - type(nemsio_gfile) :: gfileo - - data recname_all /'alnsf', 'alnwf', 'alvsf', 'alvwf', & - 'cnwat', 'crain', 'f10m', 'facsf', & - 'facwf', 'ffhh', 'ffmm', 'fricv', & - 'icec', 'icetk', 'land', 'orog', & - 'snoalb', 'sfcr', 'shdmax', 'shdmin', & - 'soill', 'soill', 'soill', 'soill', & - 'sltyp', 'soilw', 'soilw', 'soilw', & - 'soilw', 'snod', 'sotyp', 'spfh', & - 'tmp', 'tmp', 'tmp', 'tmp', & - 'tg3', 'ti', 'tmp', 'tmp', & - 'tprcp', 'veg', 'vtype', 'weasd', & - 'c0', 'cd', 'dconv', 'dtcool', & - 'qrain', 'tref', & - 'w0', 'wd', 'xs', 'xt', & - 'xtts', 'xu', 'xv', 'xz', & - 'xzts', 'zc'/ - - data reclevtyp_all /'sfc', 'sfc', 'sfc', 'sfc', & - 'sfc', 'sfc', '10 m above gnd', 'sfc', & - 'sfc', 'sfc', 'sfc', 'sfc', & - 'sfc', 'sfc', 'sfc', 'sfc', & - 'sfc', 'sfc', 'sfc', 'sfc', & - '0-10 cm down', '10-40 cm down', '40-100 cm down', '100-200 cm down', & - 'sfc', '0-10 cm down', '10-40 cm down', '40-100 cm down', & - '100-200 cm down', 'sfc', 'sfc', '2 m above gnd', & - '0-10 cm down', '10-40 cm down', '40-100 cm down', '100-200 cm down', & - 'sfc', 'sfc', '2 m above gnd', 'sfc', & - 'sfc', 'sfc', 'sfc', 'sfc', & - 'sfc', 'sfc', 'sfc', 'sfc', & - 'sfc', 'sfc', 'sfc', & - 'sfc', 'sfc', 'sfc', 'sfc', & - 'sfc', 'sfc', 'sfc', 'sfc', & - 'sfc'/ - - data reclev_all /1, 1, 1, 1, 1, & - 1, 1, 1, 1, 1, 1, & - 1, 1, 1, 1, 1, 1, & - 1, 1, 1, 1, 1, 1, & - 1, 1, 1, 1, 1, 1, & - 1, 1, 1, 1, 1, 1, & - 1, 1, 1, 1, 1, 1, & - 1, 1, 1, 1, 1, 1, & - 1, 1, 1, 1, 1, 1, & - 1, 1, 1, 1, 1, 1, 1/ - - data aryiname /'lpl'/ - - data variname /'fhzero', 'ncld', 'nsoil', 'imp_physics'/ - - data varival /6, 5, 4, 11/ - - data varrname /'dtp'/ - - data varrval /225.0/ - - data varcname /"y-direction", "z-direction"/ - - data varcval /"north2south", "bottom2top"/ - - version = 200809 - - aryival = igaus ! reduced grid definition - aryilen = jgaus - - allocate(dummy(igaus,jgaus)) - do i = 1, igaus - dummy(i,:) = float(i-1) * 360.0 / float(igaus) - enddo - - allocate(lon(igaus*jgaus)) - lon = reshape (dummy, (/igaus*jgaus/) ) - -! Call 4-byte version of splib to match latitudes in history files. - - allocate(slat(jgaus)) - allocate(wlat(jgaus)) - call splat(4, jgaus, slat, wlat) - - do i = (jgaus/2+1), jgaus - dummy(:,i) = 90.0 - (acos(slat(i)) * 180.0 / (4.0*atan(1.0))) - enddo - - do i = 1, (jgaus/2) - dummy(:,i) = -(dummy(:,(jgaus-i+1))) - enddo - - deallocate(slat, wlat) - - allocate(lat(igaus*jgaus)) - lat = reshape (dummy, (/igaus*jgaus/) ) - - deallocate(dummy) - - print* - print*, "- OPEN VCOORD FILE." - open(14, file="vcoord.txt", form='formatted', iostat=iret) - if (iret /= 0) goto 43 - - print*, "- READ VCOORD FILE." - read(14, *, iostat=iret) nvcoord, levs_vcoord - if (iret /= 0) goto 43 - - allocate(vcoord(levs_vcoord,3,2)) - vcoord = 0.0 - read(14, *, iostat=iret) ((vcoord(n,k,1), k=1,nvcoord), n=1,levs_vcoord) - if (iret /= 0) goto 43 - - close (14) - - if (trim(donst) == "yes" .or. trim(donst) == "YES") then - nrec = nrec_all - allocate(recname(nrec)) - recname = recname_all - allocate(reclevtyp(nrec)) - reclevtyp = reclevtyp_all - allocate(reclev(nrec)) - reclev = reclev_all - else - nrec = 44 - allocate(recname(nrec)) - recname = recname_all(1:nrec) - allocate(reclevtyp(nrec)) - reclevtyp = reclevtyp_all(1:nrec) - allocate(reclev(nrec)) - reclev = reclev_all(1:nrec) - endif - - call nemsio_init(iret=iret) - - print* - print*,"- OPEN GAUSSIAN NEMSIO SURFACE FILE" - - call nemsio_open(gfileo, "sfc.gaussian.analysis.file", 'write', & - modelname="FV3GFS", gdatatype="bin4", version=version, & - nmeta=8, nrec=nrec, dimx=igaus, dimy=jgaus, dimz=(levs_vcoord-1), & - nframe=0, nsoil=4, ntrac=8, jcap=-9999, & - ncldt=5, idvc=-9999, idsl=-9999, idvm=-9999, & - idrt=4, lat=lat, lon=lon, vcoord=vcoord, & - nfhour=0, nfminute=0, nfsecondn=0, & - nfsecondd=1, nfday=0, idate=idate, & - recname=recname, reclevtyp=reclevtyp, & - reclev=reclev, extrameta=.true., & - nmetavari=nmetavari, variname=variname, varival=varival, & - nmetavarr=nmetavarr, varrname=varrname, varrval=varrval, & - nmetavarc=nmetavarc, varcname=varcname, varcval=varcval, & - nmetaaryi=nmetaaryi, aryiname=aryiname, & - aryival=aryival, aryilen=aryilen, iret=iret) - if (iret /= 0) goto 44 - - deallocate (lat, lon, vcoord, recname, reclevtyp, reclev) - - allocate(the_data(igaus*jgaus)) - - print*,"- WRITE GAUSSIAN NEMSIO SURFACE FILE" - - print*,"- WRITE ALNSF" - the_data = gaussian_data%alnsf - call nemsio_writerec(gfileo, 1, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE ALNWF" - the_data = gaussian_data%alnwf - call nemsio_writerec(gfileo, 2, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE ALVSF" - the_data = gaussian_data%alvsf - call nemsio_writerec(gfileo, 3, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE ALVWF" - the_data = gaussian_data%alvwf - call nemsio_writerec(gfileo, 4, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE CANOPY" - the_data = gaussian_data%canopy - call nemsio_writerec(gfileo, 5, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE CRAIN (SRFLAG)" - the_data = gaussian_data%srflag - call nemsio_writerec(gfileo, 6, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE F10M" - the_data = gaussian_data%f10m - call nemsio_writerec(gfileo, 7, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE FACSF" - the_data = gaussian_data%facsf - call nemsio_writerec(gfileo, 8, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE FACWF" - the_data = gaussian_data%facwf - call nemsio_writerec(gfileo, 9, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE FFHH" - the_data = gaussian_data%ffhh - call nemsio_writerec(gfileo, 10, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE FFMM" - the_data = gaussian_data%ffmm - call nemsio_writerec(gfileo, 11, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE UUSTAR" - the_data = gaussian_data%uustar - call nemsio_writerec(gfileo, 12, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE FICE" - the_data = gaussian_data%fice - call nemsio_writerec(gfileo, 13, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE HICE" - the_data = gaussian_data%hice - call nemsio_writerec(gfileo, 14, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE SLMSK" - the_data = gaussian_data%slmask - call nemsio_writerec(gfileo, 15, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE OROG" - the_data = gaussian_data%orog - call nemsio_writerec(gfileo, 16, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE SNOALB" - the_data = gaussian_data%snoalb - call nemsio_writerec(gfileo, 17, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE ZORL" - the_data = gaussian_data%zorl * 0.01 ! meters - call nemsio_writerec(gfileo, 18, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE SHDMAX" - the_data = gaussian_data%shdmax - call nemsio_writerec(gfileo, 19, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE SHDMIN" - the_data = gaussian_data%shdmin - call nemsio_writerec(gfileo, 20, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE SLC" - the_data = gaussian_data%slc(:,1) - call nemsio_writerec(gfileo, 21, the_data, iret=iret) - if (iret /= 0) goto 44 - - the_data = gaussian_data%slc(:,2) - call nemsio_writerec(gfileo, 22, the_data, iret=iret) - if (iret /= 0) goto 44 - - the_data = gaussian_data%slc(:,3) - call nemsio_writerec(gfileo, 23, the_data, iret=iret) - if (iret /= 0) goto 44 - - the_data = gaussian_data%slc(:,4) - call nemsio_writerec(gfileo, 24, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE SLOPE" - the_data = gaussian_data%slope - call nemsio_writerec(gfileo, 25, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE SMC" - the_data = gaussian_data%smc(:,1) - call nemsio_writerec(gfileo, 26, the_data, iret=iret) - if (iret /= 0) goto 44 - - the_data = gaussian_data%smc(:,2) - call nemsio_writerec(gfileo, 27, the_data, iret=iret) - if (iret /= 0) goto 44 - - the_data = gaussian_data%smc(:,3) - call nemsio_writerec(gfileo, 28, the_data, iret=iret) - if (iret /= 0) goto 44 - - the_data = gaussian_data%smc(:,4) - call nemsio_writerec(gfileo, 29, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE SNWDPH" - the_data = gaussian_data%snwdph * 0.001 ! meters - call nemsio_writerec(gfileo, 30, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE STYPE" - the_data = gaussian_data%stype - call nemsio_writerec(gfileo, 31, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE Q2M" - the_data = gaussian_data%q2m - call nemsio_writerec(gfileo, 32, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE STC" - the_data = gaussian_data%stc(:,1) - call nemsio_writerec(gfileo, 33, the_data, iret=iret) - if (iret /= 0) goto 44 - - the_data = gaussian_data%stc(:,2) - call nemsio_writerec(gfileo, 34, the_data, iret=iret) - if (iret /= 0) goto 44 - - the_data = gaussian_data%stc(:,3) - call nemsio_writerec(gfileo, 35, the_data, iret=iret) - if (iret /= 0) goto 44 - - the_data = gaussian_data%stc(:,4) - call nemsio_writerec(gfileo, 36, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE TG3" - the_data = gaussian_data%tg3 - call nemsio_writerec(gfileo, 37, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE TISFC" - the_data = gaussian_data%tisfc - call nemsio_writerec(gfileo, 38, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE T2M" - the_data = gaussian_data%t2m - call nemsio_writerec(gfileo, 39, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE TSEA" - the_data = gaussian_data%tsea - call nemsio_writerec(gfileo, 40, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE TPRCP" - the_data = gaussian_data%tprcp - call nemsio_writerec(gfileo, 41, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE VFRAC" - the_data = gaussian_data%vfrac * 100.0 ! whole percent - call nemsio_writerec(gfileo, 42, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE VTYPE" - the_data = gaussian_data%vtype - call nemsio_writerec(gfileo, 43, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE SHELEG" - the_data = gaussian_data%sheleg - call nemsio_writerec(gfileo, 44, the_data, iret=iret) - if (iret /= 0) goto 44 - - if (trim(donst) == "yes" .or. trim(donst) == "YES") then - - print*,"- WRITE C0" - the_data = gaussian_data%c0 - call nemsio_writerec(gfileo, 45, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE CD" - the_data = gaussian_data%cd - call nemsio_writerec(gfileo, 46, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE DCONV" - the_data = gaussian_data%dconv - call nemsio_writerec(gfileo, 47, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE DTCOOL" - the_data = gaussian_data%dtcool - call nemsio_writerec(gfileo, 48, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE QRAIN" - the_data = gaussian_data%qrain - call nemsio_writerec(gfileo, 49, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE TREF" - the_data = gaussian_data%tref - call nemsio_writerec(gfileo, 50, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE W0" - the_data = gaussian_data%w0 - call nemsio_writerec(gfileo, 51, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE WD" - the_data = gaussian_data%wd - call nemsio_writerec(gfileo, 52, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE XS" - the_data = gaussian_data%xs - call nemsio_writerec(gfileo, 53, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE XT" - the_data = gaussian_data%xt - call nemsio_writerec(gfileo, 54, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE XTTS" - the_data = gaussian_data%xtts - call nemsio_writerec(gfileo, 55, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE XU" - the_data = gaussian_data%xu - call nemsio_writerec(gfileo, 56, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE XV" - the_data = gaussian_data%xv - call nemsio_writerec(gfileo, 57, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE XZ" - the_data = gaussian_data%xz - call nemsio_writerec(gfileo, 58, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE XZTS" - the_data = gaussian_data%xzts - call nemsio_writerec(gfileo, 59, the_data, iret=iret) - if (iret /= 0) goto 44 - - print*,"- WRITE ZC" - the_data = gaussian_data%zc - call nemsio_writerec(gfileo, 60, the_data, iret=iret) - if (iret /= 0) goto 44 - - endif - - call nemsio_close(gfileo,iret=iret) - - call nemsio_finalize() - - deallocate(the_data) - - return - - 43 continue - print*,"- ** FATAL ERROR OPENING/READING VCOORD FILE." - print*,"- IRET IS: ", iret - call errexit(17) - stop - - 44 continue - print*,"- ** FATAL ERROR WRITING GAUSSIAN NEMSIO FILE." - print*,"- IRET IS: ", iret - call errexit(15) - stop - - end subroutine write_sfc_data_nemsio - !------------------------------------------------------------------------------------------- ! Read tile data. !------------------------------------------------------------------------------------------- @@ -1646,7 +1117,7 @@ subroutine read_data_anl allocate(tile_data%smc(ijtile*num_tiles,4)) allocate(tile_data%stc(ijtile*num_tiles,4)) ! nst - if (trim(donst) == "yes" .or. trim(donst) == "YES") then + if (donst) then allocate(tile_data%c0(ijtile*num_tiles)) allocate(tile_data%cd(ijtile*num_tiles)) allocate(tile_data%dconv(ijtile*num_tiles)) @@ -1900,7 +1371,7 @@ subroutine read_data_anl print*,'- SNOALB: ',maxval(dummy),minval(dummy) tile_data%snoalb(istart:iend) = reshape(dummy, (/ijtile/)) - if (trim(donst) == "yes" .or. trim(donst) == "YES") then + if (donst) then error=nf90_inq_varid(ncid, "c_0", id_var) call netcdf_err(error, 'READING c_0 ID' ) diff --git a/src/regrid_nemsio.fd/CMakeLists.txt b/src/regrid_nemsio.fd/CMakeLists.txt deleted file mode 100644 index dbd5138c..00000000 --- a/src/regrid_nemsio.fd/CMakeLists.txt +++ /dev/null @@ -1,25 +0,0 @@ -list(APPEND fortran_src -constants.f90 -fv3_interface.f90 -gfs_nems_interface.f90 -interpolation_interface.f90 -kinds.f90 -main.f90 -mpi_interface.f90 -namelist_def.f90 -netcdfio_interface.f90 -physcons.f90 -regrid_nemsio_interface.f90 -variable_interface.f90 -) - -set(exe_name regrid_nemsio.x) -add_executable(${exe_name} ${fortran_src}) -target_link_libraries(${exe_name} PRIVATE NetCDF::NetCDF_Fortran - MPI::MPI_Fortran - bacio::bacio_4 - sp::sp_d - w3emc::w3emc_d - nemsio::nemsio) - -install(TARGETS ${exe_name} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/src/regrid_nemsio.fd/constants.f90 b/src/regrid_nemsio.fd/constants.f90 deleted file mode 100644 index 8627358e..00000000 --- a/src/regrid_nemsio.fd/constants.f90 +++ /dev/null @@ -1,314 +0,0 @@ -! this module was extracted from the GSI version operational -! at NCEP in Dec. 2007. -module constants -!$$$ module documentation block -! . . . . -! module: constants -! prgmmr: treadon org: np23 date: 2003-09-25 -! -! abstract: This module contains the definition of various constants -! used in the gsi code -! -! program history log: -! 2003-09-25 treadon - original code -! 2004-03-02 treadon - allow global and regional constants to differ -! 2004-06-16 treadon - update documentation -! 2004-10-28 treadon - replace parameter tiny=1.e-12 with tiny_r_kind -! and tiny_single -! 2004-11-16 treadon - add huge_single, huge_r_kind parameters -! 2005-01-27 cucurull - add ione -! 2005-08-24 derber - move cg_term to constants from qcmod -! 2006-03-07 treadon - add rd_over_cp_mass -! 2006-05-18 treadon - add huge_i_kind -! 2006-06-06 su - add var-qc wgtlim, change value to 0.25 (ECMWF) -! 2006-07-28 derber - add r1000 -! -! Subroutines Included: -! sub init_constants - compute derived constants, set regional/global constants -! -! Variable Definitions: -! see below -! -! attributes: -! language: f90 -! machine: ibm RS/6000 SP -! -!$$$ - use kinds, only: r_single,r_kind,i_kind - implicit none - -! Declare constants - integer(i_kind) izero,ione - real(r_kind) rearth,grav,omega,rd,rv,cp,cv,cvap,cliq - real(r_kind) csol,hvap,hfus,psat,t0c,ttp,jcal,cp_mass,cg_term - real(r_kind) fv,deg2rad,rad2deg,pi,tiny_r_kind,huge_r_kind,huge_i_kind - real(r_kind) ozcon,rozcon,tpwcon,rd_over_g,rd_over_cp,g_over_rd - real(r_kind) amsua_clw_d1,amsua_clw_d2,constoz,zero,one,two,four - real(r_kind) one_tenth,quarter,three,five,rd_over_cp_mass - real(r_kind) rearth_equator,stndrd_atmos_ps,r1000,stndrd_atmos_lapse - real(r_kind) semi_major_axis,semi_minor_axis,n_a,n_b - real(r_kind) eccentricity,grav_polar,grav_ratio - real(r_kind) grav_equator,earth_omega,grav_constant - real(r_kind) flattening,eccentricity_linear,somigliana - real(r_kind) dldt,dldti,hsub,psatk,tmix,xa,xai,xb,xbi - real(r_kind) eps,epsm1,omeps,wgtlim - real(r_kind) elocp,cpr,el2orc,cclimit,climit,epsq - real(r_kind) pcpeff0,pcpeff1,pcpeff2,pcpeff3,rcp,c0,delta - real(r_kind) h1000,factor1,factor2,rhcbot,rhctop,dx_max,dx_min,dx_inv - real(r_kind) h300,half,cmr,cws,ke2,row,rrow - real(r_single) zero_single,tiny_single,huge_single - real(r_single) rmw_mean_distance, roic_mean_distance - logical :: constants_initialized = .true. - - -! Define constants common to global and regional applications -! name value description units -! ---- ----- ----------- ----- - parameter(rearth_equator= 6.37813662e6_r_kind) ! equatorial earth radius (m) - parameter(omega = 7.2921e-5_r_kind) ! angular velocity of earth (1/s) - parameter(cp = 1.0046e+3_r_kind) ! specific heat of air @pressure (J/kg/K) - parameter(cvap = 1.8460e+3_r_kind) ! specific heat of h2o vapor (J/kg/K) - parameter(csol = 2.1060e+3_r_kind) ! specific heat of solid h2o (ice)(J/kg/K) - parameter(hvap = 2.5000e+6_r_kind) ! latent heat of h2o condensation (J/kg) - parameter(hfus = 3.3358e+5_r_kind) ! latent heat of h2o fusion (J/kg) - parameter(psat = 6.1078e+2_r_kind) ! pressure at h2o triple point (Pa) - parameter(t0c = 2.7315e+2_r_kind) ! temperature at zero celsius (K) - parameter(ttp = 2.7316e+2_r_kind) ! temperature at h2o triple point (K) - parameter(jcal = 4.1855e+0_r_kind) ! joules per calorie () - parameter(stndrd_atmos_ps = 1013.25e2_r_kind) ! 1976 US standard atmosphere ps (Pa) - -! Numeric constants - parameter(izero = 0) - parameter(ione = 1) - parameter(zero_single = 0.0_r_single) - parameter(zero = 0.0_r_kind) - parameter(one_tenth = 0.10_r_kind) - parameter(quarter= 0.25_r_kind) - parameter(one = 1.0_r_kind) - parameter(two = 2.0_r_kind) - parameter(three = 3.0_r_kind) - parameter(four = 4.0_r_kind) - parameter(five = 5.0_r_kind) - parameter(r1000 = 1000.0_r_kind) - -! Constants for gps refractivity - parameter(n_a=77.6_r_kind) !K/mb - parameter(n_b=3.73e+5_r_kind) !K^2/mb - -! Parameters below from WGS-84 model software inside GPS receivers. - parameter(semi_major_axis = 6378.1370e3_r_kind) ! (m) - parameter(semi_minor_axis = 6356.7523142e3_r_kind) ! (m) - parameter(grav_polar = 9.8321849378_r_kind) ! (m/s2) - parameter(grav_equator = 9.7803253359_r_kind) ! (m/s2) - parameter(earth_omega = 7.292115e-5_r_kind) ! (rad/s) - parameter(grav_constant = 3.986004418e14_r_kind) ! (m3/s2) - -! Derived geophysical constants - parameter(flattening = (semi_major_axis-semi_minor_axis)/semi_major_axis)!() - parameter(somigliana = & - (semi_minor_axis/semi_major_axis) * (grav_polar/grav_equator) - one)!() - parameter(grav_ratio = (earth_omega*earth_omega * & - semi_major_axis*semi_major_axis * semi_minor_axis) / grav_constant) !() - -! Derived thermodynamic constants - parameter ( dldti = cvap-csol ) - parameter ( hsub = hvap+hfus ) - parameter ( psatk = psat*0.001_r_kind ) - parameter ( tmix = ttp-20._r_kind ) - parameter ( elocp = hvap/cp ) - parameter ( rcp = one/cp ) - -! Constants used in GFS moist physics - parameter ( h300 = 300._r_kind ) - parameter ( half = 0.5_r_kind ) - parameter ( cclimit = 0.001_r_kind ) - parameter ( climit = 1.e-20_r_kind) - parameter ( epsq = 2.e-12_r_kind ) - parameter ( h1000 = 1000.0_r_kind) - parameter ( rhcbot=0.85_r_kind ) - parameter ( rhctop=0.85_r_kind ) - parameter ( dx_max=-8.8818363_r_kind ) - parameter ( dx_min=-5.2574954_r_kind ) - parameter ( dx_inv=one/(dx_max-dx_min) ) - parameter ( c0=0.002_r_kind ) - parameter ( delta=0.6077338_r_kind ) - parameter ( pcpeff0=1.591_r_kind ) - parameter ( pcpeff1=-0.639_r_kind ) - parameter ( pcpeff2=0.0953_r_kind ) - parameter ( pcpeff3=-0.00496_r_kind ) - parameter ( cmr = one/0.0003_r_kind ) - parameter ( cws = 0.025_r_kind ) - parameter ( ke2 = 0.00002_r_kind ) - parameter ( row = 1000._r_kind ) - parameter ( rrow = one/row ) - -! Constant used to process ozone - parameter ( constoz = 604229.0_r_kind) - -! Constants used in cloud liquid water correction for AMSU-A -! brightness temperatures - parameter ( amsua_clw_d1 = 0.754_r_kind ) - parameter ( amsua_clw_d2 = -2.265_r_kind ) - -! Constants used for variational qc - parameter ( wgtlim = 0.25_r_kind) ! Cutoff weight for concluding that obs has been - ! rejected by nonlinear qc. This limit is arbitrary - ! and DOES NOT affect nonlinear qc. It only affects - ! the printout which "counts" the number of obs that - ! "fail" nonlinear qc. Observations counted as failing - ! nonlinear qc are still assimilated. Their weight - ! relative to other observations is reduced. Changing - ! wgtlim does not alter the analysis, only - ! the nonlinear qc data "count" - -! Constants describing the Extended Best-Track Reanalysis [Demuth et -! al., 2008] tropical cyclone (TC) distance for regions relative to TC -! track position; units are in kilometers - - parameter (rmw_mean_distance = 64.5479412) - parameter (roic_mean_distance = 338.319656) - -contains - subroutine init_constants_derived -!$$$ subprogram documentation block -! . . . . -! subprogram: init_constants_derived set derived constants -! prgmmr: treadon org: np23 date: 2004-12-02 -! -! abstract: This routine sets derived constants -! -! program history log: -! 2004-12-02 treadon -! 2005-03-03 treadon - add implicit none -! -! input argument list: -! -! output argument list: -! -! attributes: -! language: f90 -! machine: ibm rs/6000 sp -! -!$$$ - implicit none - -! Trigonometric constants - pi = acos(-one) - deg2rad = pi/180.0_r_kind - rad2deg = one/deg2rad - cg_term = (sqrt(two*pi))/two ! constant for variational qc - tiny_r_kind = tiny(zero) - huge_r_kind = huge(zero) - tiny_single = tiny(zero_single) - huge_single = huge(zero_single) - huge_i_kind = huge(izero) - -! Geophysical parameters used in conversion of geopotential to -! geometric height - eccentricity_linear = sqrt(semi_major_axis**2 - semi_minor_axis**2) - eccentricity = eccentricity_linear / semi_major_axis - constants_initialized = .true. - - return - end subroutine init_constants_derived - - subroutine init_constants(regional) -!$$$ subprogram documentation block -! . . . . -! subprogram: init_constants set regional or global constants -! prgmmr: treadon org: np23 date: 2004-03-02 -! -! abstract: This routine sets constants specific to regional or global -! applications of the gsi -! -! program history log: -! 2004-03-02 treadon -! 2004-06-16 treadon, documentation -! 2004-10-28 treadon - use intrinsic TINY function to set value -! for smallest machine representable positive -! number -! 2004-12-03 treadon - move derived constants to init_constants_derived -! 2005-03-03 treadon - add implicit none -! -! input argument list: -! regional - if .true., set regional gsi constants; -! otherwise (.false.), use global constants -! -! output argument list: -! -! attributes: -! language: f90 -! machine: ibm rs/6000 sp -! -!$$$ - implicit none - logical regional - real(r_kind) reradius,g,r_d,r_v,cliq_wrf - - stndrd_atmos_lapse = 0.0065 - -! Define regional constants here - if (regional) then - -! Name given to WRF constants - reradius = one/6370.e03_r_kind - g = 9.81_r_kind - r_d = 287.04_r_kind - r_v = 461.6_r_kind - cliq_wrf = 4190.0_r_kind - cp_mass = 1004.67_r_kind - -! Transfer WRF constants into unified GSI constants - rearth = one/reradius - grav = g - rd = r_d - rv = r_v - cv = cp-r_d - cliq = cliq_wrf - rd_over_cp_mass = rd / cp_mass - -! Define global constants here - else - rearth = 6.3712e+6_r_kind - grav = 9.80665e+0_r_kind - rd = 2.8705e+2_r_kind - rv = 4.6150e+2_r_kind - cv = 7.1760e+2_r_kind - cliq = 4.1855e+3_r_kind - cp_mass= zero - rd_over_cp_mass = zero - endif - - -! Now define derived constants which depend on constants -! which differ between global and regional applications. - -! Constants related to ozone assimilation - ozcon = grav*21.4e-9_r_kind - rozcon= one/ozcon - -! Constant used in vertical integral for precipitable water - tpwcon = 100.0_r_kind/grav - -! Derived atmospheric constants - fv = rv/rd-one ! used in virtual temperature equation - dldt = cvap-cliq - xa = -(dldt/rv) - xai = -(dldti/rv) - xb = xa+hvap/(rv*ttp) - xbi = xai+hsub/(rv*ttp) - eps = rd/rv - epsm1 = rd/rv-one - omeps = one-eps - factor1 = (cvap-cliq)/rv - factor2 = hvap/rv-factor1*t0c - cpr = cp*rd - el2orc = hvap*hvap/(rv*cp) - rd_over_g = rd/grav - rd_over_cp = rd/cp - g_over_rd = grav/rd - - return - end subroutine init_constants - -end module constants diff --git a/src/regrid_nemsio.fd/fv3_interface.f90 b/src/regrid_nemsio.fd/fv3_interface.f90 deleted file mode 100644 index bbe558e4..00000000 --- a/src/regrid_nemsio.fd/fv3_interface.f90 +++ /dev/null @@ -1,779 +0,0 @@ -module fv3_interface - - !======================================================================= - - ! Define associated modules and subroutines - - !----------------------------------------------------------------------- - - use constants - - !----------------------------------------------------------------------- - - use gfs_nems_interface - use interpolation_interface - use mpi_interface - use namelist_def - use netcdfio_interface - use variable_interface - use nemsio_module - - !----------------------------------------------------------------------- - - implicit none - - !----------------------------------------------------------------------- - - ! Define all data and structure types for routine; these variables - ! are variables required by the subroutines within this module - - type analysis_grid - character(len=500) :: filename - character(len=500) :: filename2d - integer :: nx - integer :: ny - integer :: nz - integer :: ntime - end type analysis_grid ! type analysis_grid - - ! Define global variables - - integer n2dvar,n3dvar,ntvars,nrecs,nvvars - real(nemsio_realkind), dimension(:,:,:,:), allocatable :: fv3_var_3d - real(nemsio_realkind), dimension(:,:,:), allocatable :: fv3_var_2d - - !----------------------------------------------------------------------- - - ! Define interfaces and attributes for module routines - - private - public :: fv3_regrid_nemsio - - !----------------------------------------------------------------------- - -contains - - !----------------------------------------------------------------------- - - subroutine fv3_regrid_nemsio() - - ! Define variables computed within routine - - implicit none - type(analysis_grid) :: anlygrd(ngrids) - type(varinfo), allocatable, dimension(:) :: var_info,var_info2d,var_info3d - type(gfs_grid) :: gfs_grid - type(gridvar) :: invar,invar2 - type(gridvar) :: outvar,outvar2 - type(nemsio_meta) :: meta_nemsio2d, meta_nemsio3d - - type(esmfgrid) :: grid_bilin - type(esmfgrid) :: grid_nn - - character(len=20) :: var_name - character(len=20) :: nems_name - character(len=20) :: nems_levtyp - character(len=20) :: itrptyp - logical :: itrp_bilinear - logical :: itrp_nrstnghbr - real(nemsio_realkind), dimension(:,:), allocatable :: workgrid - real(nemsio_realkind), dimension(:), allocatable :: pk - real(nemsio_realkind), dimension(:), allocatable :: bk - real, dimension(:), allocatable :: sendbuffer,recvbuffer - integer :: fhour - integer :: ncoords - integer nems_lev,ndims,istatus,ncol,levs_fix - logical clip - - ! Define counting variables - - integer :: i, j, k, l,nlev,k2,k3,nrec - - !===================================================================== - - ! Define local variables - - call init_constants_derived() - call gfs_grid_initialize(gfs_grid) - - ! Loop through local variables - - if(mpi_procid .eq. mpi_masternode) then - print *,'variable table' - print *,'--------------' - open(912,file=trim(variable_table),form='formatted') - ntvars=0; n2dvar=0; n3dvar=0 - nrecs = 0 - loop_read: do while (istatus == 0) - read(912,199,iostat=istatus) var_name,nems_name,nems_levtyp,nems_lev,itrptyp,clip,ndims - if( istatus /= 0 ) exit loop_read - nrecs = nrecs + 1 - if(var_name(1:1) .ne. "#") then - ntvars = ntvars + 1 - ntvars = ntvars + 1 - if (ndims == 2) then - n2dvar = n2dvar+1 - else if (ndims == 3) then - n3dvar = n3dvar+1 - else - print *,'ndims must be 2 or 3 in variable_table.txt' - call mpi_abort(mpi_comm_world,-91,mpi_ierror) - stop - endif - !print *,'ntvars,n2dvar,n3dvar',ntvars,n2dvar,n3dvar - !write(6,199) var_name, nems_name,nems_levtyp,nems_lev,itrptyp,clip,ndims - endif - enddo loop_read - close(912) - print *,'nrecs,ntvars,n2dvar,n3dvar',nrecs,ntvars,n2dvar,n3dvar - endif - call mpi_bcast(nrecs,1,mpi_integer,mpi_masternode,mpi_comm_world,mpi_ierror) - call mpi_bcast(n2dvar,1,mpi_integer,mpi_masternode,mpi_comm_world,mpi_ierror) - call mpi_bcast(n3dvar,1,mpi_integer,mpi_masternode,mpi_comm_world,mpi_ierror) - call mpi_bcast(ntvars,1,mpi_integer,mpi_masternode,mpi_comm_world,mpi_ierror) - if (ntvars == 0) then - print *,'empty variable_table.txt!' - call mpi_interface_terminate() - stop - endif - allocate(var_info(ntvars)) - open(912,file=trim(variable_table),form='formatted') - k = 0 - nvvars = 0 ! number of vector variables - do nrec = 1, nrecs - read(912,199,iostat=istatus) var_name,nems_name,nems_levtyp,nems_lev,itrptyp,clip,ndims - if (var_name(1:1) .ne. "#") then - k = k + 1 - var_info(k)%var_name = var_name - var_info(k)%nems_name = nems_name - var_info(k)%nems_levtyp = nems_levtyp - var_info(k)%nems_lev = nems_lev - var_info(k)%itrptyp = itrptyp - if (itrptyp.EQ.'vector') then - nvvars=nvvars+1 - endif - var_info(k)%clip = clip - var_info(k)%ndims = ndims - if(mpi_procid .eq. mpi_masternode) then - write(6,199) var_info(k)%var_name, var_info(k)%nems_name,var_info(k)%nems_levtyp, & - var_info(k)%nems_lev,var_info(k)%itrptyp,var_info(k)%clip,var_info(k)%ndims - endif - endif - end do ! do k = 1, ntvars - ! assume vectors are in pairs - nvvars=nvvars/2 - call mpi_bcast(nvvars,1,mpi_integer,mpi_masternode,mpi_comm_world,mpi_ierror) - close(912) -199 format(a20,1x,a20,1x,a20,1x,i1,1x,a20,1x,l1,1x,i1) - allocate(var_info3d(n3dvar+2)) - allocate(var_info2d(n2dvar)) - k2 = 0 - k3 = 0 - do k=1,ntvars - if (var_info(k)%ndims == 2) then - k2 = k2 + 1 - var_info2d(k2) = var_info(k) - endif - if (var_info(k)%ndims == 3 .or. & - trim(var_info(k)%nems_name) == 'pres' .or. & - trim(var_info(k)%nems_name) == 'orog') then - k3 = k3 + 1 - var_info3d(k3) = var_info(k) - ! orography called 'hgt' in 3d file, not 'orog' - if (trim(var_info(k)%nems_name) == 'orog') then - var_info3d(k3)%nems_name = 'hgt ' - endif - endif - enddo - - - do i = 1, ngrids - anlygrd(i)%filename = analysis_filename(i) - anlygrd(i)%filename2d = analysis_filename2d(i) - call fv3_regrid_initialize(anlygrd(i)) - end do ! do i = 1, ngrids - - ! Define local variables - - ncxdim = anlygrd(1)%nx - ncydim = anlygrd(1)%ny - if (n3dvar > 0) then - nczdim = anlygrd(1)%nz - else - nczdim = 0 - endif - nctdim = anlygrd(1)%ntime - ncoords = ncxdim*ncydim - invar%ncoords = ncoords*ngrids - invar2%ncoords = ncoords*ngrids - outvar%ncoords = gfs_grid%ncoords - outvar2%ncoords = gfs_grid%ncoords - call interpolation_initialize_gridvar(invar) - call interpolation_initialize_gridvar(invar2) - call interpolation_initialize_gridvar(outvar) - call interpolation_initialize_gridvar(outvar2) - meta_nemsio3d%modelname = 'GFS' - meta_nemsio3d%version = 200509 - meta_nemsio3d%nrec = 2 + nczdim*n3dvar - meta_nemsio3d%nmeta = 5 - meta_nemsio3d%nmetavari = 3 - meta_nemsio3d%nmetaaryi = 1 - meta_nemsio3d%dimx = gfs_grid%nlons - meta_nemsio3d%dimy = gfs_grid%nlats - meta_nemsio3d%dimz = nczdim - meta_nemsio3d%jcap = ntrunc - meta_nemsio3d%nsoil = 4 - meta_nemsio3d%nframe = 0 - meta_nemsio3d%ntrac = 3 - meta_nemsio3d%idrt = 4 - meta_nemsio3d%ncldt = 3 - meta_nemsio3d%idvc = 2 - meta_nemsio3d%idvm = 2 - meta_nemsio3d%idsl = 1 - meta_nemsio3d%idate(1:6) = 0 - meta_nemsio3d%idate(7) = 1 - read(forecast_timestamp(9:10),'(i2)') meta_nemsio3d%idate(4) - read(forecast_timestamp(7:8), '(i2)') meta_nemsio3d%idate(3) - read(forecast_timestamp(5:6), '(i2)') meta_nemsio3d%idate(2) - read(forecast_timestamp(1:4), '(i4)') meta_nemsio3d%idate(1) - meta_nemsio2d = meta_nemsio3d - meta_nemsio2d%nrec = n2dvar - call mpi_barrier(mpi_comm_world,mpi_ierror) - call gfs_nems_meta_initialization(meta_nemsio2d,var_info2d,gfs_grid) - call gfs_nems_meta_initialization(meta_nemsio3d,var_info3d,gfs_grid) - - ! Allocate memory for local variables - - if(.not. allocated(fv3_var_2d) .and. n2dvar > 0) & - & allocate(fv3_var_2d(ngrids,ncxdim,ncydim)) - if (mpi_nprocs /= nczdim+1) then - call mpi_barrier(mpi_comm_world, mpi_ierror) - if (mpi_procid .eq. mpi_masternode) then - print *,'number of mpi tasks must be equal to number of levels + 1' - print *,'mpi procs = ',mpi_nprocs,' levels = ',nczdim - endif - call mpi_interface_terminate() - stop - endif - !print *,'allocate fv3_var_3d',ngrids,ncxdim,ncydim,nczdim,mpi_procid - if(.not. allocated(fv3_var_3d) .and. n3dvar > 0) & - & allocate(fv3_var_3d(ngrids,ncxdim,ncydim,nczdim)) - !print *,'done allocating fv3_var_3d',ngrids,ncxdim,ncydim,nczdim,mpi_procid - - ! Check local variable and proceed accordingly - - call mpi_barrier(mpi_comm_world,mpi_ierror) - if(mpi_procid .eq. mpi_masternode) then - - ! Allocate memory for local variables - - if (n3dvar > 0) then - if(.not. allocated(pk)) allocate(pk(nczdim+1)) - if(.not. allocated(bk)) allocate(bk(nczdim+1)) - - ! Define local variables - - if (trim(gfs_hyblevs_filename) == 'NOT USED' ) then - call netcdfio_values_1d(anlygrd(1)%filename,'pk',pk) - call netcdfio_values_1d(anlygrd(1)%filename,'bk',bk) - else - open(913,file=trim(gfs_hyblevs_filename),form='formatted') - read(913,*) ncol, levs_fix - if (levs_fix /= (nczdim+1) ) then - call mpi_barrier(mpi_comm_world, mpi_ierror) - print *,'levs in ', trim(gfs_hyblevs_filename), ' not equal to',(nczdim+1) - call mpi_interface_terminate() - stop - endif - do k=nczdim+1,1,-1 - read(913,*) pk(k),bk(k) - enddo - close(913) - endif - if (minval(pk) < -1.e10 .or. minval(bk) < -1.e10) then - print *,'pk,bk not found in netcdf file..' - meta_nemsio3d%vcoord = -9999._nemsio_realkind - else - ! Loop through local variable - - do k = 1, nczdim + 1 - - ! Define local variables - - meta_nemsio3d%vcoord((nczdim + 1) - k + 1,1,1) = pk(k) - meta_nemsio3d%vcoord((nczdim + 1) - k + 1,2,1) = bk(k) - - end do ! do k = 1, nczdim + 1 - endif - endif - - end if ! if(mpi_procid .eq. mpi_masternode) - - ! initialize/read in interpolation weight - - grid_bilin%filename = esmf_bilinear_filename - call interpolation_initialize_esmf(grid_bilin) - - grid_nn%filename = esmf_neareststod_filename - call interpolation_initialize_esmf(grid_nn) - - do l = 1, nctdim - - ncrec = l ! time level to read from netcdf file - - ! Define local variables - - call fv3_grid_fhour(anlygrd(1),meta_nemsio2d%nfhour) - call fv3_grid_fhour(anlygrd(1),meta_nemsio3d%nfhour) - meta_nemsio3d%nfminute = int(0.0) - meta_nemsio3d%nfsecondn = int(0.0) - meta_nemsio3d%nfsecondd = int(1.0) - meta_nemsio3d%fhour = meta_nemsio3d%nfhour - meta_nemsio2d%nfminute = int(0.0) - meta_nemsio2d%nfsecondn = int(0.0) - meta_nemsio2d%nfsecondd = int(1.0) - meta_nemsio2d%fhour = meta_nemsio2d%nfhour - - ! initialize nemsio file. - if(mpi_procid .eq. mpi_masternode) then - call gfs_nems_initialize(meta_nemsio2d, meta_nemsio3d) - end if - - ! wait here. - call mpi_barrier(mpi_comm_world,mpi_ierror) - - ! Loop through local variables - k2=1 - do k = 1, ntvars - nvvars - - ! Define local variables - - itrp_bilinear = .false. - itrp_nrstnghbr = .false. - - ! Do 2D variables. - - if(var_info(k2)%ndims .eq. 2) then - - ! Check local variable and proceed accordingly - - if(mpi_procid .eq. mpi_masternode) then - - ! Check local variable and proceed accordingly - - call fv3_grid_read(anlygrd(1:ngrids), var_info(k2)%var_name,.true.,.false.) - - call interpolation_define_gridvar(invar,ncxdim,ncydim, ngrids,fv3_var_2d) - if (trim(var_info(k2)%nems_name) == 'pres') then - ! interpolate in exner(pressure) - invar%var = (invar%var/stndrd_atmos_ps)**(rd_over_g*stndrd_atmos_lapse) - end if - - if(var_info(k2)%itrptyp .eq. 'bilinear') then - call interpolation_esmf(invar,outvar,grid_bilin, .false.) - end if - - if(var_info(k2)%itrptyp .eq. 'nrstnghbr') then - call interpolation_esmf(invar,outvar,grid_nn, .true.) - end if - - if (trim(var_info(k2)%nems_name) == 'pres') then - outvar%var = stndrd_atmos_ps*(outvar%var**(g_over_rd/stndrd_atmos_lapse)) - end if - - if(var_info(k2)%itrptyp .eq. 'vector') then - ! read in u winds - call fv3_grid_read(anlygrd(1:ngrids), var_info(k2)%var_name,.true.,.false.) - call interpolation_define_gridvar(invar,ncxdim,ncydim,ngrids,fv3_var_2d) - ! read in v winds - call fv3_grid_read(anlygrd(1:ngrids), var_info(k2+1)%var_name,.true.,.false.) - call interpolation_define_gridvar(invar2,ncxdim,ncydim,ngrids,fv3_var_2d) - call interpolation_esmf_vect(invar,invar2,grid_bilin,outvar,outvar2) - end if - - ! Clip variable to zero if desired. - if(var_info(k2)%clip) call variable_clip(outvar%var) - - ! Write to NEMSIO file. - call gfs_nems_write('2d',real(outvar%var), & - var_info(k2)%nems_name,var_info(k2)%nems_levtyp,var_info(k2)%nems_lev) - if (trim(var_info(k2)%nems_name) == 'pres' .or. & - trim(var_info(k2)%nems_name) == 'orog' .or. & - trim(var_info(k2)%nems_name) == 'hgt') then - ! write surface height and surface pressure to 3d file. - ! (surface height called 'orog' in nemsio bin4, 'hgt' in - ! grib) - if (trim(var_info(k2)%nems_name) == 'orog') then - call gfs_nems_write('3d',real(outvar%var), & - 'hgt ',var_info(k2)%nems_levtyp,1) - else - call gfs_nems_write('3d',real(outvar%var), & - var_info(k2)%nems_name,var_info(k2)%nems_levtyp,1) - endif - endif - if(var_info(k2)%itrptyp .eq. 'vector') then ! write v winds - call gfs_nems_write('2d',real(outvar2%var), & - var_info(k2+1)%nems_name,var_info(k2+1)%nems_levtyp,var_info(k2+1)%nems_lev) - endif - end if ! if(mpi_procid .eq. mpi_masternode) - - ! Define local variables - call mpi_barrier(mpi_comm_world,mpi_ierror) - - end if ! if(var_info(k2)%ndims .eq. 2) - - ! Do 3D variables. - - if(var_info(k2)%ndims .eq. 3) then - - ! read 3d grid on master node, send to other tasks - if(mpi_procid .eq. mpi_masternode) then - call fv3_grid_read(anlygrd(1:ngrids), var_info(k2)%var_name,.false.,.true.) - do nlev=1,nczdim - call mpi_send(fv3_var_3d(1,1,1,nlev),ngrids*ncxdim*ncydim,mpi_real,& - nlev,1,mpi_comm_world,mpi_errorstatus,mpi_ierror) - enddo - if(trim(adjustl(var_info(k2)%itrptyp)) .eq. 'vector') then ! winds - call mpi_barrier(mpi_comm_world,mpi_ierror) - call fv3_grid_read(anlygrd(1:ngrids), var_info(k2+1)%var_name,.false.,.true.) - do nlev=1,nczdim - call mpi_send(fv3_var_3d(1,1,1,nlev),ngrids*ncxdim*ncydim,mpi_real,& - nlev,1,mpi_comm_world,mpi_errorstatus,mpi_ierror) - enddo - endif - else if (mpi_procid .le. nczdim) then - ! do interpolation, one level on each task. - call mpi_recv(fv3_var_3d(1,1,1,mpi_procid),ngrids*ncxdim*ncydim,mpi_real,& - 0,1,mpi_comm_world,mpi_errorstatus,mpi_ierror) - - call interpolation_define_gridvar(invar,ncxdim,ncydim, ngrids,fv3_var_3d(:,:,:,mpi_procid)) - - if(var_info(k2)%itrptyp .eq. 'bilinear') then - call interpolation_esmf(invar,outvar,grid_bilin, .false.) - end if ! if(var_info(k2)%itrptyp .eq. 'bilinear') - - if(var_info(k2)%itrptyp .eq. 'nrstnghbr') then - call interpolation_esmf(invar,outvar,grid_nn, .true.) - end if ! if(var_info(k2)%itrptyp .eq. 'nrstnghbr') - - if(trim(adjustl(var_info(k2)%itrptyp)) .eq. 'vector') then ! winds - call mpi_barrier(mpi_comm_world,mpi_ierror) - call mpi_recv(fv3_var_3d(1,1,1,mpi_procid),ngrids*ncxdim*ncydim,mpi_real,& - 0,1,mpi_comm_world,mpi_errorstatus,mpi_ierror) - call interpolation_define_gridvar(invar2,ncxdim,ncydim,ngrids,fv3_var_3d(:,:,:,mpi_procid)) - call interpolation_esmf_vect(invar,invar2,grid_bilin,outvar,outvar2) - endif - - if(var_info(k2)%clip) call variable_clip(outvar%var(:)) - - end if ! if(mpi_procid .ne. mpi_masternode .and. & - ! mpi_procid .le. nczdim) - - ! gather results back on root node to write out. - - if (mpi_procid == mpi_masternode) then - ! receive one level of interpolated data on root task. - if (.not. allocated(workgrid)) allocate(workgrid(gfs_grid%ncoords,nczdim)) - if (.not. allocated(recvbuffer)) allocate(recvbuffer(gfs_grid%ncoords)) - do nlev=1,nczdim - call mpi_recv(recvbuffer,gfs_grid%ncoords,mpi_real,& - nlev,1,mpi_comm_world,mpi_errorstatus,mpi_ierror) - workgrid(:,nlev) = recvbuffer - enddo - deallocate(recvbuffer) - else - ! send one level of interpolated data to root task. - if (.not. allocated(sendbuffer)) allocate(sendbuffer(gfs_grid%ncoords)) - sendbuffer(:) = outvar%var(:) - call mpi_send(sendbuffer,gfs_grid%ncoords,mpi_real,& - 0,1,mpi_comm_world,mpi_errorstatus,mpi_ierror) - endif - - ! Write to NEMSIO file. - - if(mpi_procid .eq. mpi_masternode) then - - ! Loop through local variable - - do j = 1, nczdim - - ! Define local variables - - call gfs_nems_write('3d',workgrid(:,nczdim - j + 1), & - & var_info(k2)%nems_name,var_info(k2)%nems_levtyp, & - & j) - - end do ! do j = 1, nczdim - - end if ! if(mpi_procid .eq. mpi_masternode) - - if(trim(adjustl(var_info(k2)%itrptyp)) .eq. 'vector') then ! winds - if (mpi_procid == mpi_masternode) then - ! receive one level of interpolated data on root task. - if (.not. allocated(workgrid)) allocate(workgrid(gfs_grid%ncoords,nczdim)) - if (.not. allocated(recvbuffer)) allocate(recvbuffer(gfs_grid%ncoords)) - do nlev=1,nczdim - call mpi_recv(recvbuffer,gfs_grid%ncoords,mpi_real,& - nlev,1,mpi_comm_world,mpi_errorstatus,mpi_ierror) - workgrid(:,nlev) = recvbuffer - enddo - deallocate(recvbuffer) - else - ! send one level of interpolated data to root task. - if (.not. allocated(sendbuffer)) allocate(sendbuffer(gfs_grid%ncoords)) - sendbuffer(:) = outvar2%var(:) - call mpi_send(sendbuffer,gfs_grid%ncoords,mpi_real,& - 0,1,mpi_comm_world,mpi_errorstatus,mpi_ierror) - endif - - ! Write to NEMSIO file. - - if(mpi_procid .eq. mpi_masternode) then - - do j = 1, nczdim - - call gfs_nems_write('3d',workgrid(:,nczdim - j + 1), & - & var_info(k2+1)%nems_name,var_info(k2+1)%nems_levtyp, & - & j) - end do ! do j = 1, nczdim - - end if ! if(mpi_procid .eq. mpi_masternode) - endif - - ! wait here - - call mpi_barrier(mpi_comm_world,mpi_ierror) - - end if ! if(var_info(k2)%ndims .eq. 3) - if(var_info(k2)%itrptyp .eq. 'vector') then ! skip v record here - k2=k2+1 - endif - k2=k2+1 - end do ! do k = 1, ntvars - - ! Wait here. - - call mpi_barrier(mpi_comm_world,mpi_ierror) - - ! Finalize and cleanup - - if(mpi_procid .eq. mpi_masternode) then - call gfs_nems_finalize() - end if - call mpi_barrier(mpi_comm_world,mpi_ierror) - if(allocated(workgrid)) deallocate(workgrid) - - end do ! do l = 1, nctdim - - -!===================================================================== - - end subroutine fv3_regrid_nemsio - - !======================================================================= - - ! fv3_regrid_initialize.f90: - - !----------------------------------------------------------------------- - - subroutine fv3_regrid_initialize(grid) - - ! Define variables passed to routine - - implicit none - type(analysis_grid) :: grid - - !===================================================================== - - ! Define local variables - - call netcdfio_dimension(grid%filename,'grid_xt',grid%nx) - call netcdfio_dimension(grid%filename,'grid_yt',grid%ny) - if (n3dvar > 0) then - call netcdfio_dimension(grid%filename,'pfull',grid%nz) - else - grid%nz = 0 - endif - call netcdfio_dimension(grid%filename,'time',grid%ntime) - - !===================================================================== - - end subroutine fv3_regrid_initialize - - !======================================================================= - - ! fv3_grid_read.f90: - - !----------------------------------------------------------------------- - - subroutine fv3_grid_read(anlygrd,varname,is_2d,is_3d) - - ! Define variables passed to subroutine - - type(analysis_grid) :: anlygrd(ngrids) - character(len=20) :: varname - logical :: is_2d - logical :: is_3d - - ! Define counting variables - - integer :: i, j, k - - !===================================================================== - - ! Loop through local variable - - do k = 1, ngrids - - ! Check local variable and proceed accordingly - - if(debug) write(6,500) ncrec, k - if(is_2d) then - - ! Define local variables - - ! orog and psfc are in 3d file. - if (trim(varname) == 'orog' .or. trim(varname) == 'psfc') then - call netcdfio_values_2d(anlygrd(k)%filename,varname, & - & fv3_var_2d(k,:,:)) - else - call netcdfio_values_2d(anlygrd(k)%filename2d,varname, & - & fv3_var_2d(k,:,:)) - endif - - end if ! if(is_2d) - - ! Check local variable and proceed accordingly - - if(is_3d) then - - ! Define local variables - - call netcdfio_values_3d(anlygrd(k)%filename,varname, & - & fv3_var_3d(k,:,:,:)) - - end if ! if(is_3d) - - end do ! do k = 1, ngrids - - !===================================================================== - - ! Define format statements - -500 format('FV3_GRID_READ: Time record = ', i6, '; Cubed sphere face = ', & - & i1,'.') - - !===================================================================== - - end subroutine fv3_grid_read - - !======================================================================= - - ! fv3_grid_fhour.f90: - - !----------------------------------------------------------------------- - - subroutine fv3_grid_fhour(grid,fhour) - - ! Define variables passed to routine - - implicit none - type(analysis_grid) :: grid - integer :: fhour - - ! Define variables computed within routine - - real(nemsio_realkind) :: workgrid(grid%ntime) - real(nemsio_realkind) :: start_jday - real(nemsio_realkind) :: fcst_jday - integer :: year - integer :: month - integer :: day - integer :: hour - integer :: minute - integer :: second, iw3jdn - character(len=80) timeunits - - !===================================================================== - - ! Define local variables - - read(forecast_timestamp(1:4), '(i4)') year - read(forecast_timestamp(5:6), '(i2)') month - read(forecast_timestamp(7:8), '(i2)') day - read(forecast_timestamp(9:10),'(i2)') hour - minute = 0; second = 0 - - ! Compute local variables - - ! 'flux day' (days since dec 31 1900) - !call date2wnday(start_jday,year,month,day) - ! same as above, but valid after 2099 - start_jday=real(iw3jdn(year,month,day)-iw3jdn(1900,12,31)) - start_jday = start_jday + real(hour)/24.0 + real(minute)/1440.0 + & - & real(second)/86400.0 - - ! Define local variables - - call netcdfio_values_1d(grid%filename,'time',workgrid) - call netcdfio_variable_attr(grid%filename,'time','units',timeunits) - - ! Compute local variables - - ! ncrec is a global variable in the netcdfio-interface module - if (timeunits(1:4) == 'days') then - fcst_jday = start_jday + workgrid(ncrec) - else if (timeunits(1:5) == 'hours') then - fcst_jday = start_jday + workgrid(ncrec)/24. - else if (timeunits(1:7) == 'seconds') then - fcst_jday = start_jday + workgrid(ncrec)/86400.0 - else - print *,'unrecognized time units',trim(timeunits) - call mpi_interface_terminate() - stop - endif - fhour = nint((86400*(fcst_jday - start_jday))/3600.0) - - !===================================================================== - - end subroutine fv3_grid_fhour - -! SUBROUTINE DATE2WNDAY(WDAY, IYR,MON,IDY) -! IMPLICIT NONE -! INTEGER IYR,MON,IDY -! REAL WDAY -!! -!!********** -!!* -!! 1) CONVERT DATE INTO 'FLUX DAY'. -!! -!! 2) THE 'FLUX DAY' IS THE NUMBER OF DAYS SINCE 001/1901 (WHICH IS -!! FLUX DAY 1.0). -!! FOR EXAMPLE: -!! A) IYR=1901,MON=1,IDY=1, REPRESENTS 0000Z HRS ON 01/01/1901 -!! SO WDAY WOULD BE 1.0. -!! A) IYR=1901,MON=1,IDY=2, REPRESENTS 0000Z HRS ON 02/01/1901 -!! SO WDAY WOULD BE 2.0. -!! YEAR MUST BE NO LESS THAN 1901.0, AND NO GREATER THAN 2099.0. -!! NOTE THAT YEAR 2000 IS A LEAP YEAR (BUT 1900 AND 2100 ARE NOT). -!! -!! 3) ALAN J. WALLCRAFT, NAVAL RESEARCH LABORATORY, JULY 2002. -!!* -!!********** -!! -! INTEGER NLEAP -! REAL WDAY1 -! REAL MONTH(13) -! DATA MONTH / 0, 31, 59, 90, 120, 151, 181, & -! 212, 243, 273, 304, 334, 365 / -!! FIND THE RIGHT YEAR. -! NLEAP = (IYR-1901)/4 -! WDAY = 365.0*(IYR-1901) + NLEAP + MONTH(MON) + IDY -! IF (MOD(IYR,4).EQ.0 .AND. MON.GT.2) THEN -! WDAY = WDAY + 1.0 -! ENDIF -! END SUBROUTINE DATE2WNDAY - - !======================================================================= - -end module fv3_interface diff --git a/src/regrid_nemsio.fd/gfs_nems_interface.f90 b/src/regrid_nemsio.fd/gfs_nems_interface.f90 deleted file mode 100644 index aa1305dc..00000000 --- a/src/regrid_nemsio.fd/gfs_nems_interface.f90 +++ /dev/null @@ -1,595 +0,0 @@ -module gfs_nems_interface - - !======================================================================= - - ! Define associated modules and subroutines - - !----------------------------------------------------------------------- - - use constants - use kinds - - !----------------------------------------------------------------------- - - use interpolation_interface - use mpi_interface - use namelist_def - use nemsio_module - use netcdfio_interface - use variable_interface - - !----------------------------------------------------------------------- - - implicit none - - !----------------------------------------------------------------------- - - ! Define all data and structure types for routine; these variables - ! are variables required by the subroutines within this module - - type gfs_grid - real(r_kind), dimension(:,:), allocatable :: rlon - real(r_kind), dimension(:,:), allocatable :: rlat - integer :: ncoords - integer :: nlons - integer :: nlats - integer :: nz - end type gfs_grid ! type gfs_grid - - type nemsio_meta - character(nemsio_charkind), dimension(:), allocatable :: recname - character(nemsio_charkind), dimension(:), allocatable :: reclevtyp - character(nemsio_charkind), dimension(:), allocatable :: variname - character(nemsio_charkind), dimension(:), allocatable :: varr8name - character(nemsio_charkind), dimension(:), allocatable :: aryiname - character(nemsio_charkind), dimension(:), allocatable :: aryr8name - character(nemsio_charkind8) :: gdatatype - character(nemsio_charkind8) :: modelname - real(nemsio_realkind), dimension(:,:,:), allocatable :: vcoord - real(nemsio_realkind), dimension(:), allocatable :: lon - real(nemsio_realkind), dimension(:), allocatable :: lat - integer(nemsio_intkind), dimension(:,:), allocatable :: aryival - integer(nemsio_intkind), dimension(:), allocatable :: reclev - integer(nemsio_intkind), dimension(:), allocatable :: varival - integer(nemsio_intkind), dimension(:), allocatable :: aryilen - integer(nemsio_intkind), dimension(:), allocatable :: aryr8len - integer(nemsio_intkind) :: idate(7) - integer(nemsio_intkind) :: version - integer(nemsio_intkind) :: nreo_vc - integer(nemsio_intkind) :: nrec - integer(nemsio_intkind) :: nmeta - integer(nemsio_intkind) :: nmetavari - integer(nemsio_intkind) :: nmetaaryi - integer(nemsio_intkind) :: nfhour - integer(nemsio_intkind) :: nfminute - integer(nemsio_intkind) :: nfsecondn - integer(nemsio_intkind) :: nfsecondd - integer(nemsio_intkind) :: jcap - integer(nemsio_intkind) :: dimx - integer(nemsio_intkind) :: dimy - integer(nemsio_intkind) :: dimz - integer(nemsio_intkind) :: nframe - integer(nemsio_intkind) :: nsoil - integer(nemsio_intkind) :: ntrac - integer(nemsio_intkind) :: ncldt - integer(nemsio_intkind) :: idvc - integer(nemsio_intkind) :: idsl - integer(nemsio_intkind) :: idvm - integer(nemsio_intkind) :: idrt - integer(nemsio_intkind) :: fhour - end type nemsio_meta ! type nemsio_meta - - !----------------------------------------------------------------------- - - ! Define global variables - - type(nemsio_gfile) :: gfile2d,gfile3d - integer :: nemsio_iret - - !----------------------------------------------------------------------- - - ! Define interfaces and attributes for module routines - - private - public :: gfs_grid_initialize - public :: gfs_grid_cleanup - public :: gfs_grid - public :: gfs_nems_meta_initialization - public :: gfs_nems_meta_cleanup - public :: gfs_nems_initialize - public :: gfs_nems_finalize - public :: gfs_nems_write - public :: nemsio_meta - -contains - - !======================================================================= - - ! gfs_nems_write.f90: - - !----------------------------------------------------------------------- - - subroutine gfs_nems_write(c2dor3d,nems_data,nems_varname,nems_levtyp,nems_lev) - - ! Define variables passed to routine - - character(nemsio_charkind) :: nems_varname - character(nemsio_charkind) :: nems_levtyp - real(nemsio_realkind) :: nems_data(:) - integer(nemsio_intkind) :: nems_lev - character(len=2) :: c2dor3d - - !===================================================================== - - ! Define local variables - - if (c2dor3d == '2d') then - call nemsio_writerecv(gfile2d,trim(adjustl(nems_varname)),levtyp= & - & trim(adjustl(nems_levtyp)),lev=nems_lev,data=nems_data, & - & iret=nemsio_iret) - else if (c2dor3d == '3d') then - call nemsio_writerecv(gfile3d,trim(adjustl(nems_varname)),levtyp= & - & trim(adjustl(nems_levtyp)),lev=nems_lev,data=nems_data, & - & iret=nemsio_iret) - else - nemsio_iret=-99 - endif - - ! Check local variable and proceed accordingly - - if(debug) write(6,500) c2dor3d,trim(adjustl(nems_varname)), nemsio_iret, & - & nems_lev, minval(nems_data), maxval(nems_data) - - !===================================================================== - - ! Define format statements - -500 format('GFS_NEMS_WRITE',a2,': NEMS I/O name = ', a, '; writerecv return ', & - & 'code = ', i5,'; level = ', i3, '; (min,max) = (', f13.5,f13.5, & - & ').') - if (nemsio_iret /= 0) then - print *,'nemsio_writerecv failed, stopping...' - call mpi_interface_terminate() - stop - endif - - !===================================================================== - - end subroutine gfs_nems_write - - !======================================================================= - - ! gfs_nems_meta_initialization.f90: - - !----------------------------------------------------------------------- - - subroutine gfs_nems_meta_initialization(meta_nemsio,var_info,grid) - - ! Define variables passed to routine - - type(nemsio_meta) :: meta_nemsio - type(varinfo) :: var_info(:) - type(gfs_grid) :: grid - - ! Define variables computed within routine - - integer :: offset - integer :: n2dvar - integer :: n3dvar - - ! Define counting variables - - integer :: i, j, k - - !===================================================================== - - ! Allocate memory for local variables - - if(.not. allocated(meta_nemsio%recname)) & - & allocate(meta_nemsio%recname(meta_nemsio%nrec)) - if(.not. allocated(meta_nemsio%reclevtyp)) & - & allocate(meta_nemsio%reclevtyp(meta_nemsio%nrec)) - if(.not. allocated(meta_nemsio%reclev)) & - & allocate(meta_nemsio%reclev(meta_nemsio%nrec)) - if(.not. allocated(meta_nemsio%variname)) & - & allocate(meta_nemsio%variname(meta_nemsio%nmetavari)) - if(.not. allocated(meta_nemsio%varival)) & - & allocate(meta_nemsio%varival(meta_nemsio%nmetavari)) - if(.not. allocated(meta_nemsio%aryiname)) & - & allocate(meta_nemsio%aryiname(meta_nemsio%nmetaaryi)) - if(.not. allocated(meta_nemsio%aryilen)) & - & allocate(meta_nemsio%aryilen(meta_nemsio%nmetaaryi)) - if(.not. allocated(meta_nemsio%vcoord)) & - & allocate(meta_nemsio%vcoord(meta_nemsio%dimz+1,3,2)) - if(.not. allocated(meta_nemsio%aryival)) & - & allocate(meta_nemsio%aryival(grid%nlats/2, & - & meta_nemsio%nmetaaryi)) - if(.not. allocated(meta_nemsio%lon)) & - & allocate(meta_nemsio%lon(grid%ncoords)) - if(.not. allocated(meta_nemsio%lat)) & - & allocate(meta_nemsio%lat(grid%ncoords)) - meta_nemsio%vcoord(:,:,:)=0.0 - ! Define local variables - - meta_nemsio%lon = & - & reshape(grid%rlon,(/grid%ncoords/)) - meta_nemsio%lat = & - & reshape(grid%rlat,(/grid%ncoords/)) - meta_nemsio%aryilen(1) = grid%nlats/2 - meta_nemsio%aryiname(1) = 'lpl' - meta_nemsio%aryival(1:grid%nlats/2,1) = grid%nlons - k = 0 - - ! Loop through local variable - offset = 0 - n3dvar = 0 - n2dvar = 0 - - - do i = 1, size(var_info) - - ! Check local variable and proceed accordingly - - if(var_info(i)%ndims .eq. 2) then - - ! Define local variables - - k = k + 1 - meta_nemsio%reclev(k) = var_info(i)%nems_lev - meta_nemsio%recname(k) = trim(adjustl(var_info(i)%nems_name)) - meta_nemsio%reclevtyp(k) = trim(adjustl(var_info(i)%nems_levtyp)) - n2dvar = k - - else if(var_info(i)%ndims .eq. 3) then - - ! Loop through local variable - - meta_nemsio%variname(1) = 'LEVS' - meta_nemsio%varival(1) = meta_nemsio%dimz - meta_nemsio%variname(2) = 'NVCOORD' - meta_nemsio%varival(2) = 2 - meta_nemsio%variname(3) = 'IVS' - meta_nemsio%varival(3) = 200509 - do k = 1, meta_nemsio%dimz - - ! Define local variables - - meta_nemsio%reclev(k+n2dvar+offset) = k - meta_nemsio%recname(k+n2dvar+offset) = & - & trim(adjustl(var_info(i)%nems_name)) - meta_nemsio%reclevtyp(k+n2dvar+offset) = & - & trim(adjustl(var_info(i)%nems_levtyp)) - - end do ! do k = 1, nczdim - - ! Define local variables - - n3dvar = n3dvar + 1 - offset = nczdim*n3dvar - - end if ! if(var_info(i)%ndims .eq. 3) - - end do ! do i = 1, size(var_info) - - !===================================================================== - - end subroutine gfs_nems_meta_initialization - - !======================================================================= - - ! gfs_nems_meta_cleanup.f90: - - !----------------------------------------------------------------------- - - subroutine gfs_nems_meta_cleanup(meta_nemsio2d,meta_nemsio3d) - - ! Define variables passed to routine - - type(nemsio_meta) :: meta_nemsio2d,meta_nemsio3d - - !===================================================================== - - ! Deallocate memory for local variables - - if(allocated(meta_nemsio2d%recname)) & - & deallocate(meta_nemsio2d%recname) - if(allocated(meta_nemsio2d%reclevtyp)) & - & deallocate(meta_nemsio2d%reclevtyp) - if(allocated(meta_nemsio2d%reclev)) & - & deallocate(meta_nemsio2d%reclev) - if(allocated(meta_nemsio2d%variname)) & - & deallocate(meta_nemsio2d%variname) - if(allocated(meta_nemsio2d%aryiname)) & - & deallocate(meta_nemsio2d%aryiname) - if(allocated(meta_nemsio2d%aryival)) & - & deallocate(meta_nemsio2d%aryival) - if(allocated(meta_nemsio2d%aryilen)) & - & deallocate(meta_nemsio2d%aryilen) - if(allocated(meta_nemsio2d%vcoord)) & - & deallocate(meta_nemsio2d%vcoord) - if(allocated(meta_nemsio2d%lon)) & - & deallocate(meta_nemsio2d%lon) - if(allocated(meta_nemsio2d%lat)) & - & deallocate(meta_nemsio2d%lat) - if(allocated(meta_nemsio3d%recname)) & - & deallocate(meta_nemsio3d%recname) - if(allocated(meta_nemsio3d%reclevtyp)) & - & deallocate(meta_nemsio3d%reclevtyp) - if(allocated(meta_nemsio3d%reclev)) & - & deallocate(meta_nemsio3d%reclev) - if(allocated(meta_nemsio3d%variname)) & - & deallocate(meta_nemsio3d%variname) - if(allocated(meta_nemsio3d%aryiname)) & - & deallocate(meta_nemsio3d%aryiname) - if(allocated(meta_nemsio3d%aryival)) & - & deallocate(meta_nemsio3d%aryival) - if(allocated(meta_nemsio3d%aryilen)) & - & deallocate(meta_nemsio3d%aryilen) - if(allocated(meta_nemsio3d%vcoord)) & - & deallocate(meta_nemsio3d%vcoord) - if(allocated(meta_nemsio3d%lon)) & - & deallocate(meta_nemsio3d%lon) - if(allocated(meta_nemsio3d%lat)) & - & deallocate(meta_nemsio3d%lat) - - !===================================================================== - - end subroutine gfs_nems_meta_cleanup - - !======================================================================= - - ! gfs_nems_initialize.f90: - - !----------------------------------------------------------------------- - - subroutine gfs_nems_initialize(meta_nemsio2d, meta_nemsio3d) - - ! Define variables passed to routine - - type(nemsio_meta) :: meta_nemsio2d,meta_nemsio3d - character(len=500) :: filename - character(len=7) :: suffix - - !===================================================================== - - ! Define local variables - - call nemsio_init(iret=nemsio_iret) - write(suffix,500) meta_nemsio2d%nfhour - filename = trim(adjustl(datapathout2d))//suffix - meta_nemsio2d%gdatatype = trim(adjustl(nemsio_opt2d)) - meta_nemsio3d%gdatatype = trim(adjustl(nemsio_opt3d)) - call nemsio_open(gfile2d,trim(adjustl(filename)),'write', & - & iret=nemsio_iret, & - & modelname=trim(adjustl(meta_nemsio2d%modelname)), & - & version=meta_nemsio2d%version, & - & gdatatype=meta_nemsio2d%gdatatype, & - & jcap=meta_nemsio2d%jcap, & - & dimx=meta_nemsio2d%dimx, & - & dimy=meta_nemsio2d%dimy, & - & dimz=meta_nemsio2d%dimz, & - & idate=meta_nemsio2d%idate, & - & nrec=meta_nemsio2d%nrec, & - & nframe=meta_nemsio2d%nframe, & - & idrt=meta_nemsio2d%idrt, & - & ncldt=meta_nemsio2d%ncldt, & - & idvc=meta_nemsio2d%idvc, & - & idvm=meta_nemsio2d%idvm, & - & idsl=meta_nemsio2d%idsl, & - & nfhour=meta_nemsio2d%fhour, & - & nfminute=meta_nemsio2d%nfminute, & - & nfsecondn=meta_nemsio2d%nfsecondn, & - & nfsecondd=meta_nemsio2d%nfsecondd, & - & extrameta=.true., & - & nmetaaryi=meta_nemsio2d%nmetaaryi, & - & recname=meta_nemsio2d%recname, & - & reclevtyp=meta_nemsio2d%reclevtyp, & - & reclev=meta_nemsio2d%reclev, & - & aryiname=meta_nemsio2d%aryiname, & - & aryilen=meta_nemsio2d%aryilen, & - & aryival=meta_nemsio2d%aryival, & - & vcoord=meta_nemsio2d%vcoord) - write(suffix,500) meta_nemsio3d%nfhour - filename = trim(adjustl(datapathout3d))//suffix - call nemsio_open(gfile3d,trim(adjustl(filename)),'write', & - & iret=nemsio_iret, & - & modelname=trim(adjustl(meta_nemsio3d%modelname)), & - & version=meta_nemsio3d%version, & - & gdatatype=meta_nemsio3d%gdatatype, & - & jcap=meta_nemsio3d%jcap, & - & dimx=meta_nemsio3d%dimx, & - & dimy=meta_nemsio3d%dimy, & - & dimz=meta_nemsio3d%dimz, & - & idate=meta_nemsio3d%idate, & - & nrec=meta_nemsio3d%nrec, & - & nframe=meta_nemsio3d%nframe, & - & idrt=meta_nemsio3d%idrt, & - & ncldt=meta_nemsio3d%ncldt, & - & idvc=meta_nemsio3d%idvc, & - & idvm=meta_nemsio3d%idvm, & - & idsl=meta_nemsio3d%idsl, & - & nfhour=meta_nemsio3d%fhour, & - & nfminute=meta_nemsio3d%nfminute, & - & nfsecondn=meta_nemsio3d%nfsecondn, & - & nfsecondd=meta_nemsio3d%nfsecondd, & - & extrameta=.true., & - & nmetaaryi=meta_nemsio3d%nmetaaryi, & - & recname=meta_nemsio3d%recname, & - & reclevtyp=meta_nemsio3d%reclevtyp, & - & reclev=meta_nemsio3d%reclev, & - & aryiname=meta_nemsio3d%aryiname, & - & aryilen=meta_nemsio3d%aryilen, & - & aryival=meta_nemsio3d%aryival, & - & variname=meta_nemsio3d%variname, & - & varival=meta_nemsio3d%varival, & - & nmetavari=meta_nemsio3d%nmetavari, & - & vcoord=meta_nemsio3d%vcoord) - - !===================================================================== - - ! Define format statements - -500 format('.fhr',i3.3) - - !===================================================================== - - end subroutine gfs_nems_initialize - - !======================================================================= - - ! gfs_nems_finalize.f90: - - !----------------------------------------------------------------------- - - subroutine gfs_nems_finalize() - - !===================================================================== - - ! Define local variables - - call nemsio_close(gfile2d,iret=nemsio_iret) - call nemsio_close(gfile3d,iret=nemsio_iret) - - !===================================================================== - - end subroutine gfs_nems_finalize - - !======================================================================= - - ! gfs_grid_initialize.f90: - - !----------------------------------------------------------------------- - - subroutine gfs_grid_initialize(grid) - - ! Define variables passed to routine - - type(gfs_grid) :: grid - - ! Define variables computed within routine - - real(r_kind), dimension(:), allocatable :: slat - real(r_kind), dimension(:), allocatable :: wlat - real(r_kind), dimension(:), allocatable :: workgrid - - ! Define counting variables - - integer :: i, j, k - - !===================================================================== - - ! Check local variable and proceed accordingly - - if(mpi_procid .eq. mpi_masternode) then - - ! Define local variables - - call init_constants_derived() - - ! Check local variable and proceed accordingly - - ! Define local variables - - grid%nlons = nlons - grid%nlats = nlats - - end if ! if(mpi_procid .eq. mpi_masternode) - - ! Define local variables - - call mpi_barrier(mpi_comm_world,mpi_ierror) - - ! Broadcast all necessary variables to compute nodes - - call mpi_bcast(grid%nlons,1,mpi_integer,mpi_masternode,mpi_comm_world, & - & mpi_ierror) - call mpi_bcast(grid%nlats,1,mpi_integer,mpi_masternode,mpi_comm_world, & - & mpi_ierror) - - ! Allocate memory for local variables - - if(.not. allocated(grid%rlon)) & - & allocate(grid%rlon(grid%nlons,grid%nlats)) - if(.not. allocated(grid%rlat)) & - & allocate(grid%rlat(grid%nlons,grid%nlats)) - - ! Check local variable and proceed accordingly - - if(mpi_procid .eq. mpi_masternode) then - - ! Allocate memory for local variables - - if(.not. allocated(slat)) allocate(slat(grid%nlats)) - if(.not. allocated(wlat)) allocate(wlat(grid%nlats)) - if(.not. allocated(workgrid)) allocate(workgrid(grid%nlats)) - - ! Compute local variables - - grid%ncoords = grid%nlons*grid%nlats - call splat(grid%nlats,slat,wlat) - workgrid = acos(slat) - pi/2.0 - - ! Loop through local variable - - do j = 1, grid%nlats - - ! Loop through local variable - - do i = 1, grid%nlons - - ! Compute local variables - - grid%rlon(i,j) = (i-1)*(360./grid%nlons)*deg2rad - grid%rlat(i,j) = workgrid(grid%nlats - j + 1) - - end do ! do i = 1, grid%nlons - - end do ! do j = 1, grid%nlats - - ! Deallocate memory for local variables - - if(allocated(slat)) deallocate(slat) - if(allocated(wlat)) deallocate(wlat) - if(allocated(workgrid)) deallocate(workgrid) - - endif ! if(mpi_procid .eq. mpi_masternode) - - ! Broadcast all necessary variables to compute nodes - - call mpi_bcast(grid%ncoords,1,mpi_integer,mpi_masternode, & - & mpi_comm_world,mpi_ierror) - call mpi_bcast(grid%rlon,grid%ncoords,mpi_real,mpi_masternode, & - & mpi_comm_world,mpi_ierror) - call mpi_bcast(grid%rlat,grid%ncoords,mpi_real,mpi_masternode, & - & mpi_comm_world,mpi_ierror) - - !===================================================================== - - end subroutine gfs_grid_initialize - - !======================================================================= - - ! gfs_grid_cleanup.f90: - - !----------------------------------------------------------------------- - - subroutine gfs_grid_cleanup(grid) - - ! Define variables passed to routine - - type(gfs_grid) :: grid - - !===================================================================== - - ! Deallocate memory for local variables - - if(allocated(grid%rlon)) deallocate(grid%rlon) - if(allocated(grid%rlat)) deallocate(grid%rlat) - - !===================================================================== - - end subroutine gfs_grid_cleanup - - !======================================================================= - -end module gfs_nems_interface diff --git a/src/regrid_nemsio.fd/interpolation_interface.f90 b/src/regrid_nemsio.fd/interpolation_interface.f90 deleted file mode 100644 index 775d1a7c..00000000 --- a/src/regrid_nemsio.fd/interpolation_interface.f90 +++ /dev/null @@ -1,335 +0,0 @@ -module interpolation_interface - - !======================================================================= - - ! Define associated modules and subroutines - - !----------------------------------------------------------------------- - - use constants - use kinds - - !----------------------------------------------------------------------- - - use namelist_def - use netcdf - use netcdfio_interface - use mpi_interface - - !----------------------------------------------------------------------- - - implicit none - - !----------------------------------------------------------------------- - - ! Define interfaces and attributes for module routines - - private - public :: interpolation_initialize_gridvar - public :: interpolation_initialize_esmf - public :: interpolation_define_gridvar - public :: interpolation_define_gridvar_out - public :: interpolation_esmf - public :: interpolation_esmf_vect - public :: gridvar - public :: esmfgrid - - !----------------------------------------------------------------------- - - ! Define all data and structure types for routine; these variables - ! are variables required by the subroutines within this module - - type esmfgrid - character(len=500) :: filename - real(r_double), dimension(:), allocatable :: s - integer, dimension(:), allocatable :: col - integer, dimension(:), allocatable :: row - real(r_double), dimension(:), allocatable :: inlats - real(r_double), dimension(:), allocatable :: inlons - real(r_double), dimension(:), allocatable :: outlats - real(r_double), dimension(:), allocatable :: outlons - integer :: n_s,n_a,n_b - end type esmfgrid ! type esmfgrid - - type gridvar - logical, dimension(:), allocatable :: check - real(r_double), dimension(:), allocatable :: var - integer :: ncoords - integer :: nx - integer :: ny - end type gridvar ! type gridvar - - ! Define global variables - - integer :: ncfileid - integer :: ncvarid - integer :: ncdimid - integer :: ncstatus - - !----------------------------------------------------------------------- - -contains - - !======================================================================= - - subroutine interpolation_define_gridvar(grid,xdim,ydim,ngrid,input) -! collapses the cubed grid into a 1-d array -! Define variables passed to routine - - use nemsio_module, only: nemsio_realkind - integer,intent(in) :: ngrid - integer,intent(in) :: xdim,ydim - type(gridvar),intent(inout) :: grid - real(nemsio_realkind),intent(in) :: input(ngrid,xdim,ydim) - -! locals - integer :: i,j,k,ncount - - ncount = 1 - do k = 1, ngrid - do j = 1, ydim - do i = 1, xdim - grid%var(ncount) = input(k,i,j) - ncount = ncount + 1 - end do - end do - end do - - - end subroutine interpolation_define_gridvar - -!======================================================================= - - - subroutine interpolation_define_gridvar_out(grid,xdim,ydim,output) -! make a 2-d array for output - ! Define variables passed to routine - - integer,intent(in) :: xdim,ydim - type(gridvar),intent(in) :: grid - real(r_double),intent(out) :: output(xdim,ydim) - -! locals - integer :: i,j,ncount - - ncount = 1 - do j = 1, ydim - do i = 1, xdim - output(j,i) = grid%var(ncount) - ncount = ncount + 1 - enddo - enddo - - end subroutine interpolation_define_gridvar_out - - !======================================================================= - - subroutine interpolation_initialize_gridvar(grid) - - ! Define variables passed to routine - - type(gridvar) :: grid - - allocate(grid%var(grid%ncoords)) - - end subroutine interpolation_initialize_gridvar - - -!======================================================================= - - subroutine interpolation_initialize_esmf(grid) - - ! Define variables passed to routine - - type(esmfgrid) :: grid - - !===================================================================== - - ! Define local variables - - ncstatus = nf90_open(path=trim(adjustl(grid%filename)),mode= & - & nf90_nowrite,ncid=ncfileid) - ncstatus = nf90_inq_dimid(ncfileid,'n_s',ncdimid) - ncstatus = nf90_inquire_dimension(ncfileid,ncdimid,len=grid%n_s) - ncstatus = nf90_inq_dimid(ncfileid,'n_a',ncdimid) - ncstatus = nf90_inquire_dimension(ncfileid,ncdimid,len=grid%n_a) - ncstatus = nf90_inq_dimid(ncfileid,'n_b',ncdimid) - ncstatus = nf90_inquire_dimension(ncfileid,ncdimid,len=grid%n_b) - - - ! Allocate memory for local variables - - allocate(grid%s(grid%n_s)) - allocate(grid%row(grid%n_s)) - allocate(grid%col(grid%n_s)) - - allocate(grid%inlats(grid%n_a)) - allocate(grid%inlons(grid%n_a)) - allocate(grid%outlats(grid%n_b)) - allocate(grid%outlons(grid%n_b)) - - ncstatus = nf90_inq_varid(ncfileid,'col',ncvarid) - ncstatus = nf90_get_var(ncfileid,ncvarid,grid%col) - ncstatus = nf90_inq_varid(ncfileid,'row',ncvarid) - ncstatus = nf90_get_var(ncfileid,ncvarid,grid%row) - ncstatus = nf90_inq_varid(ncfileid,'S',ncvarid) - ncstatus = nf90_get_var(ncfileid,ncvarid,grid%s) - ncstatus = nf90_inq_varid(ncfileid,'yc_a',ncvarid) - ncstatus = nf90_get_var(ncfileid,ncvarid,grid%inlats) - ncstatus = nf90_inq_varid(ncfileid,'xc_a',ncvarid) - ncstatus = nf90_get_var(ncfileid,ncvarid,grid%inlons) - where(grid%inlons .LT. 0.0) - grid%inlons=360+grid%inlons - endwhere - ncstatus = nf90_inq_varid(ncfileid,'yc_b',ncvarid) - ncstatus = nf90_get_var(ncfileid,ncvarid,grid%outlats) - ncstatus = nf90_inq_varid(ncfileid,'xc_b',ncvarid) - ncstatus = nf90_get_var(ncfileid,ncvarid,grid%outlons) - ncstatus = nf90_close(ncfileid) - - !===================================================================== - - end subroutine interpolation_initialize_esmf - - -!======================================================================= - - - subroutine interpolation_esmf(invar,outvar,grid,is_nrstnghbr) - - ! Define variables passed to routine - - type(gridvar) :: invar - type(gridvar) :: outvar - logical :: is_nrstnghbr - - type(esmfgrid) :: grid - - integer :: i, j, k, l - - outvar%var = dble(0.0) - - if(is_nrstnghbr) then - do i = 1, grid%n_s - outvar%var(grid%row(i)) = invar%var(grid%col(i)) - enddo - else - do i = 1, grid%n_s - outvar%var(grid%row(i)) = outvar%var(grid%row(i)) + grid%s(i)*invar%var(grid%col(i)) - end do - end if - - end subroutine interpolation_esmf -!===================================================================== - - subroutine interpolation_esmf_vect(invaru,invarv,grid,outvaru,outvarv) - - ! Define variables passed to routine - - type(gridvar) :: invaru,invarv - type(gridvar) :: outvaru,outvarv - type(esmfgrid) :: grid - - integer :: i, j, k, l - real(r_double) :: cxy,sxy,urot,vrot - - - outvaru%var = dble(0.0) - outvarv%var = dble(0.0) - - do i = 1, grid%n_s - CALL MOVECT(grid%inlats(grid%col(i)),grid%inlons(grid%col(i)),& - grid%outlats(grid%row(i)),grid%outlons(grid%row(i)),& - cxy,sxy) - urot=cxy*invaru%var(grid%col(i))-sxy*invarv%var(grid%col(i)) - vrot=sxy*invaru%var(grid%col(i))+cxy*invarv%var(grid%col(i)) - outvaru%var(grid%row(i)) = outvaru%var(grid%row(i)) + grid%s(i)*urot - outvarv%var(grid%row(i)) = outvarv%var(grid%row(i)) + grid%s(i)*vrot - - end do - - end subroutine interpolation_esmf_vect - -!===================================================================== - - SUBROUTINE MOVECT(FLAT,FLON,TLAT,TLON,CROT,SROT) -!$$$ SUBPROGRAM DOCUMENTATION BLOCK -! -! SUBPROGRAM: MOVECT MOVE A VECTOR ALONG A GREAT CIRCLE -! PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10 -! -! ABSTRACT: THIS SUBPROGRAM PROVIDES THE ROTATION PARAMETERS -! TO MOVE A VECTOR ALONG A GREAT CIRCLE FROM ONE -! POSITION TO ANOTHER WHILE CONSERVING ITS ORIENTATION -! WITH RESPECT TO THE GREAT CIRCLE. THESE ROTATION -! PARAMETERS ARE USEFUL FOR VECTOR INTERPOLATION. -! -! PROGRAM HISTORY LOG: -! 96-04-10 IREDELL -! 1999-04-08 IREDELL GENERALIZE PRECISION -! -! USAGE: CALL MOVECT(FLAT,FLON,TLAT,TLON,CROT,SROT) -! -! INPUT ARGUMENT LIST: -! FLAT - REAL LATITUDE IN DEGREES FROM WHICH TO MOVE THE VECTOR -! FLON - REAL LONGITUDE IN DEGREES FROM WHICH TO MOVE THE VECTOR -! TLAT - REAL LATITUDE IN DEGREES TO WHICH TO MOVE THE VECTOR -! TLON - REAL LONGITUDE IN DEGREES TO WHICH TO MOVE THE VECTOR -! -! OUTPUT ARGUMENT LIST: -! CROT - REAL CLOCKWISE VECTOR ROTATION COSINE -! SROT - REAL CLOCKWISE VECTOR ROTATION SINE -! (UTO=CROT*UFROM-SROT*VFROM; -! VTO=SROT*UFROM+CROT*VFROM) -! -! ATTRIBUTES: -! LANGUAGE: FORTRAN 90 -! -!$$$ - IMPLICIT NONE -! - INTEGER, PARAMETER :: KD=SELECTED_REAL_KIND(15,45) -! - REAL(KIND=r_double), INTENT(IN ) :: FLAT, FLON - REAL(KIND=r_double), INTENT(IN ) :: TLAT, TLON - REAL(KIND=r_double), INTENT( OUT) :: CROT, SROT -! - REAL(KIND=r_double), PARAMETER :: CRDLIM=0.9999999 - REAL(KIND=r_double), PARAMETER :: PI=3.14159265358979 - REAL(KIND=r_double), PARAMETER :: DPR=180./PI -! - REAL(KIND=r_double) :: CTLAT,STLAT,CFLAT,SFLAT - REAL(KIND=r_double) :: CDLON,SDLON,CRD - REAL(KIND=r_double) :: SRD2RN,STR,CTR,SFR,CFR -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -! COMPUTE COSINE OF THE RADIAL DISTANCE BETWEEN THE POINTS. - CTLAT=COS(TLAT/DPR) - STLAT=SIN(TLAT/DPR) -CFLAT=COS(FLAT/DPR) - SFLAT=SIN(FLAT/DPR) - CDLON=COS((FLON-TLON)/DPR) - SDLON=SIN((FLON-TLON)/DPR) - CRD=STLAT*SFLAT+CTLAT*CFLAT*CDLON -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -! COMPUTE ROTATIONS AT BOTH POINTS WITH RESPECT TO THE GREAT CIRCLE -! AND COMBINE THEM TO GIVE THE TOTAL VECTOR ROTATION PARAMETERS. - IF(ABS(CRD).LE.CRDLIM) THEN - SRD2RN=-1/(1-CRD**2) - STR=CFLAT*SDLON - CTR=CFLAT*STLAT*CDLON-SFLAT*CTLAT - SFR=CTLAT*SDLON - CFR=CTLAT*SFLAT*CDLON-STLAT*CFLAT - CROT=SRD2RN*(CTR*CFR-STR*SFR) - SROT=SRD2RN*(CTR*SFR+STR*CFR) -! USE A DIFFERENT APPROXIMATION FOR NEARLY COINCIDENT POINTS. -! MOVING VECTORS TO ANTIPODAL POINTS IS AMBIGUOUS ANYWAY. - ELSE - CROT=CDLON - SROT=SDLON*STLAT - ENDIF -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - END SUBROUTINE MOVECT - - !======================================================================= - -end module interpolation_interface diff --git a/src/regrid_nemsio.fd/kinds.f90 b/src/regrid_nemsio.fd/kinds.f90 deleted file mode 100644 index 11c93b98..00000000 --- a/src/regrid_nemsio.fd/kinds.f90 +++ /dev/null @@ -1,107 +0,0 @@ -! this module was extracted from the GSI version operational -! at NCEP in Dec. 2007. -module kinds -!$$$ module documentation block -! . . . . -! module: kinds -! prgmmr: treadon org: np23 date: 2004-08-15 -! -! abstract: Module to hold specification kinds for variable declaration. -! This module is based on (copied from) Paul vanDelst's -! type_kinds module found in the community radiative transfer -! model -! -! module history log: -! 2004-08-15 treadon -! -! Subroutines Included: -! -! Functions Included: -! -! remarks: -! The numerical data types defined in this module are: -! i_byte - specification kind for byte (1-byte) integer variable -! i_short - specification kind for short (2-byte) integer variable -! i_long - specification kind for long (4-byte) integer variable -! i_llong - specification kind for double long (8-byte) integer variable -! r_single - specification kind for single precision (4-byte) real variable -! r_double - specification kind for double precision (8-byte) real variable -! r_quad - specification kind for quad precision (16-byte) real variable -! -! i_kind - generic specification kind for default integer -! r_kind - generic specification kind for default floating point -! -! -! attributes: -! language: f90 -! machine: ibm RS/6000 SP -! -!$$$ end documentation block - implicit none - private - -! Integer type definitions below - -! Integer types - integer, parameter, public :: i_byte = selected_int_kind(1) ! byte integer - integer, parameter, public :: i_short = selected_int_kind(4) ! short integer - integer, parameter, public :: i_long = selected_int_kind(8) ! long integer - integer, parameter, private :: llong_t = selected_int_kind(16) ! llong integer - integer, parameter, public :: i_llong = max( llong_t, i_long ) - -! Expected 8-bit byte sizes of the integer kinds - integer, parameter, public :: num_bytes_for_i_byte = 1 - integer, parameter, public :: num_bytes_for_i_short = 2 - integer, parameter, public :: num_bytes_for_i_long = 4 - integer, parameter, public :: num_bytes_for_i_llong = 8 - -! Define arrays for default definition - integer, parameter, private :: num_i_kinds = 4 - integer, parameter, dimension( num_i_kinds ), private :: integer_types = (/ & - i_byte, i_short, i_long, i_llong /) - integer, parameter, dimension( num_i_kinds ), private :: integer_byte_sizes = (/ & - num_bytes_for_i_byte, num_bytes_for_i_short, & - num_bytes_for_i_long, num_bytes_for_i_llong /) - -! Default values -! **** CHANGE THE FOLLOWING TO CHANGE THE DEFAULT INTEGER TYPE KIND *** - integer, parameter, public :: default_integer = 3 ! 1=byte, - ! 2=short, - ! 3=long, - ! 4=llong - integer, parameter, public :: i_kind = integer_types( default_integer ) - integer, parameter, public :: num_bytes_for_i_kind = & - integer_byte_sizes( default_integer ) - - -! Real definitions below - -! Real types - integer, parameter, public :: r_single = selected_real_kind(6) ! single precision - integer, parameter, public :: r_double = selected_real_kind(15) ! double precision - integer, parameter, private :: quad_t = selected_real_kind(20) ! quad precision - integer, parameter, public :: r_quad = max( quad_t, r_double ) - -! Expected 8-bit byte sizes of the real kinds - integer, parameter, public :: num_bytes_for_r_single = 4 - integer, parameter, public :: num_bytes_for_r_double = 8 - integer, parameter, public :: num_bytes_for_r_quad = 16 - -! Define arrays for default definition - integer, parameter, private :: num_r_kinds = 3 - integer, parameter, dimension( num_r_kinds ), private :: real_kinds = (/ & - r_single, r_double, r_quad /) - integer, parameter, dimension( num_r_kinds ), private :: real_byte_sizes = (/ & - num_bytes_for_r_single, num_bytes_for_r_double, & - num_bytes_for_r_quad /) - -! Default values -! **** CHANGE THE FOLLOWING TO CHANGE THE DEFAULT REAL TYPE KIND *** - integer, parameter, public :: default_real = 1 ! 1=single, - ! 2=double, - ! 3=quad - integer, parameter, public :: r_kind = real_kinds( default_real ) - integer, parameter, public :: num_bytes_for_r_kind = & - real_byte_sizes( default_real ) - -end module kinds diff --git a/src/regrid_nemsio.fd/main.f90 b/src/regrid_nemsio.fd/main.f90 deleted file mode 100644 index f3dfe4ef..00000000 --- a/src/regrid_nemsio.fd/main.f90 +++ /dev/null @@ -1,92 +0,0 @@ -program regrid_nemsio_main - - !===================================================================== - - !$$$ PROGRAM DOCUMENTATION BLOCK - ! - ! ABSTRACT: - ! - ! This routine provides an interface between the National Oceanic - ! and Atmospheric Administration (NOAA) National Centers for - ! Environmental Prediction (NCEP) implemented NOAA Environmental - ! Modeling System (NEMS) input/output file format and the native - ! FV3 cubed sphere grid. - ! - ! NOTES: - ! - ! * Uses interpolation weights generated by - ! Earth-System Modeling Framework (ESMF) remapping utilities. - ! - ! PRGMMR: Winterbottom - ! ORG: ESRL/PSD1 - ! DATE: 2016-02-02 - ! - ! PROGRAM HISTORY LOG: - ! - ! 2016-02-02 Initial version. Henry R. Winterbottom - ! 2016-11-01 Modifed by Jeff Whitaker. - ! - !$$$ - - !===================================================================== - - ! Define associated modules and subroutines - - !--------------------------------------------------------------------- - - use kinds - - !--------------------------------------------------------------------- - - use mpi_interface - use fv3_interface - use namelist_def - use constants - - !--------------------------------------------------------------------- - - implicit none - - !===================================================================== - - ! Define variables computed within routine - - real(r_kind) :: exectime_start - real(r_kind) :: exectime_finish - - !===================================================================== - - ! Define local variables - - call mpi_interface_initialize() - call init_constants(.false.) - - if(mpi_procid .eq. mpi_masternode) then - - call cpu_time(exectime_start) - - end if - call mpi_barrier(mpi_comm_world,mpi_ierror) - - call namelistparams() - call fv3_regrid_nemsio() - - - if(mpi_procid .eq. mpi_masternode) then - - call cpu_time(exectime_finish) - write(6,500) exectime_finish - exectime_start - - end if ! if(mpi_procid .eq. mpi_masternode) - - call mpi_barrier(mpi_comm_world,mpi_ierror) - call mpi_interface_terminate() - - !===================================================================== - ! Define format statements - -500 format('MAIN: Execution time: ', f13.5, ' seconds.') - - !===================================================================== - -end program regrid_nemsio_main diff --git a/src/regrid_nemsio.fd/mpi_interface.f90 b/src/regrid_nemsio.fd/mpi_interface.f90 deleted file mode 100644 index 2e6c5c7a..00000000 --- a/src/regrid_nemsio.fd/mpi_interface.f90 +++ /dev/null @@ -1,89 +0,0 @@ -module mpi_interface - - !======================================================================= - - use kinds - - !----------------------------------------------------------------------- - - implicit none - - !----------------------------------------------------------------------- - - ! Define necessary include files - - include "mpif.h" - - !----------------------------------------------------------------------- - - ! Define global variables - - character :: mpi_nodename(mpi_max_processor_name) - character :: mpi_noderequest - logical :: abort_mpi - integer(kind=4), dimension(:), allocatable :: mpi_ranks - integer(kind=4) :: mpi_errorstatus(mpi_status_size) - integer(kind=4) :: mpi_masternode - integer(kind=4) :: mpi_slavenode - integer(kind=4) :: mpi_ierror - integer(kind=4) :: mpi_ierrorcode - integer(kind=4) :: mpi_procid - integer(kind=4) :: mpi_nprocs - integer(kind=4) :: mpi_node_source - integer(kind=4) :: mpi_node_destination - integer(kind=4) :: mpi_loopcount - integer(kind=4) :: mpi_request - integer(kind=4) :: mpi_group_user - integer(kind=4) :: mpi_group_nprocs - integer(kind=4) :: mpi_group_procid - integer(kind=4) :: mpi_group_begin - integer(kind=4) :: mpi_group_end - - !----------------------------------------------------------------------- - -contains - - !======================================================================= - - ! mpi_interface_initialize.f90: - - !----------------------------------------------------------------------- - - subroutine mpi_interface_initialize() - - !===================================================================== - - ! Define local variables - - call mpi_init(mpi_ierror) - call mpi_comm_rank(mpi_comm_world,mpi_procid,mpi_ierror) - call mpi_comm_size(mpi_comm_world,mpi_nprocs,mpi_ierror) - mpi_masternode = 0 - abort_mpi = .false. - - !===================================================================== - - end subroutine mpi_interface_initialize - - !======================================================================= - - ! mpi_interface_terminate.f90: - - !----------------------------------------------------------------------- - - subroutine mpi_interface_terminate() - - !===================================================================== - - ! Define local variables - - !call mpi_abort(mpi_comm_world,ierror_code,mpi_ierror) - call mpi_finalize(mpi_ierror) - - !===================================================================== - - end subroutine mpi_interface_terminate - - !======================================================================= - -end module mpi_interface diff --git a/src/regrid_nemsio.fd/namelist_def.f90 b/src/regrid_nemsio.fd/namelist_def.f90 deleted file mode 100644 index ff15a335..00000000 --- a/src/regrid_nemsio.fd/namelist_def.f90 +++ /dev/null @@ -1,181 +0,0 @@ -module namelist_def - - !======================================================================= - - ! Define associated modules and subroutines - - !----------------------------------------------------------------------- - - use kinds - - !----------------------------------------------------------------------- - - use mpi_interface - - !----------------------------------------------------------------------- - - implicit none - - !----------------------------------------------------------------------- - - ! Define global variables - - integer, parameter :: max_ngrids = 12 - character(len=500) :: analysis_filename(max_ngrids) = 'NOT USED' - character(len=500) :: analysis_filename2d(max_ngrids) = 'NOT USED' - character(len=500) :: gfs_hyblevs_filename = 'NOT USED' - character(len=500) :: esmf_neareststod_filename = 'NOT USED' - character(len=500) :: esmf_bilinear_filename = 'NOT USED' - character(len=500) :: variable_table = 'NOT USED' - character(len=500) :: datapathout2d = './' - character(len=500) :: datapathout3d = './' - character(len=19) :: forecast_timestamp = '0000-00-00_00:00:00' - character(len=4) :: nemsio_opt = 'bin4' - character(len=4) :: nemsio_opt2d = 'none' - character(len=4) :: nemsio_opt3d = 'none' - logical :: is_ugrid2sgrid = .false. - logical :: debug = .false. - integer :: nlons = 0 - integer :: nlats = 0 - integer :: ntrunc = 0 - integer :: ngrids = 0 - namelist /share/ debug, nlons,nlats,ntrunc,datapathout2d,datapathout3d, & - analysis_filename, forecast_timestamp, nemsio_opt, nemsio_opt2d, nemsio_opt3d, & - analysis_filename2d, variable_table - - namelist /interpio/ esmf_bilinear_filename, esmf_neareststod_filename, gfs_hyblevs_filename - - !--------------------------------------------------------------------- - -contains - - !===================================================================== - - ! namelistparams.f90: - - !--------------------------------------------------------------------- - - subroutine namelistparams() - - ! Define variables computed within routine - - logical :: is_it_there - integer :: unit_nml - - ! Define counting variables - - integer :: i, j, k - - !=================================================================== - - ! Define local variables - - unit_nml = 9 - is_it_there = .false. - inquire(file='regrid-nemsio.input',exist = is_it_there) - - ! Check local variable and proceed accordingly - - if(is_it_there) then - - ! Define local variables - - open(file = 'regrid-nemsio.input', & - unit = unit_nml , & - status = 'old' , & - form = 'formatted' , & - action = 'read' , & - access = 'sequential' ) - read(unit_nml,NML = share) - read(unit_nml,NML = interpio) - close(unit_nml) - if (nemsio_opt2d == 'none') nemsio_opt2d=nemsio_opt - if (nemsio_opt3d == 'none') nemsio_opt3d=nemsio_opt - - ! Loop through local variable - - do i = 1, max_ngrids - - ! Check local variable and proceed accordingly - - if(analysis_filename(i) .ne. 'NOT USED') then - - ! Define local variables - - ngrids = ngrids + 1 - - end if ! if(analysis_filename(i) .ne. 'NOT USED') - - end do ! do i = 1, max_ngrids - - else ! if(is_it_there) - - ! Define local variables - - if(mpi_procid .eq. mpi_masternode) write(6,500) - call mpi_barrier(mpi_comm_world,mpi_ierror) - call mpi_interface_terminate() - - end if ! if(.not. is_it_there) - - !=================================================================== - - ! Check local variable and proceed accordingly - - if(mpi_procid .eq. mpi_masternode) then - - ! Define local variables - - write(6,*) '&SHARE' - write(6,*) 'DEBUG = ', debug - write(6,*) 'ANALYSIS_FILENAME = ' - do k = 1, ngrids - write(6,*) trim(adjustl(analysis_filename(k))) - ! if analysis_filename2d not specified, set to analysis_filename - if (trim(analysis_filename2d(k)) == 'NOT USED') then - analysis_filename2d(k) = analysis_filename(k) - endif - end do ! do k = 1, ngrids - write(6,*) 'ANALYSIS_FILENAME2D = ' - do k = 1, ngrids - write(6,*) trim(adjustl(analysis_filename2d(k))) - end do ! do k = 1, ngrids - write(6,*) 'VARIABLE_TABLE = ', & - & trim(adjustl(variable_table)) - write(6,*) 'FORECAST_TIMESTAMP = ', forecast_timestamp - write(6,*) 'OUTPUT DATAPATH (2d) = ', & - & trim(adjustl(datapathout2d)) - write(6,*) 'OUTPUT DATAPATH (3d) = ', & - & trim(adjustl(datapathout3d)) - write(6,*) 'NEMSIO_OPT (2d) = ', nemsio_opt2d - write(6,*) 'NEMSIO_OPT (3d) = ', nemsio_opt3d - write(6,*) '/' - write(6,*) '&INTERPIO' - write(6,*) 'ESMF_BILINEAR_FILENAME = ', & - & trim(adjustl(esmf_bilinear_filename)) - write(6,*) 'ESMF_NEARESTSTOD_FILENAME = ', & - & trim(adjustl(esmf_neareststod_filename)) - write(6,*) 'GFS_HYBLEVS_FILENAME = ', & - & trim(adjustl(gfs_hyblevs_filename)) - write(6,*) '/' - - end if ! if(mpi_procid .eq. mpi_masternode) - - ! Define local variables - - call mpi_barrier(mpi_comm_world,mpi_ierror) - - !=================================================================== - - ! Define format statements - -500 format('NAMELISTPARAMS: regrid-nemsio.input not found in the', & - & ' current working directory. ABORTING!!!!') - - !=================================================================== - - end subroutine namelistparams - - !===================================================================== - -end module namelist_def diff --git a/src/regrid_nemsio.fd/netcdfio_interface.f90 b/src/regrid_nemsio.fd/netcdfio_interface.f90 deleted file mode 100644 index 473b765c..00000000 --- a/src/regrid_nemsio.fd/netcdfio_interface.f90 +++ /dev/null @@ -1,592 +0,0 @@ -module netcdfio_interface - - !======================================================================= - - ! Define associated modules and subroutines - - !----------------------------------------------------------------------- - - use kinds - - !----------------------------------------------------------------------- - - use namelist_def - use netcdf - use mpi_interface - - !----------------------------------------------------------------------- - - implicit none - - !----------------------------------------------------------------------- - - ! Define global variables - - logical :: ncstatic - integer :: ncrec - integer :: ncxdim - integer :: ncydim - integer :: nczdim - integer :: nctdim - integer :: ncfileid - integer :: ncvarid - integer :: ncdimid - integer :: ncstatus - - !----------------------------------------------------------------------- - - ! Define interfaces and attributes for module routines - - private - interface netcdfio_values_1d - module procedure netcdfio_values_1d_dblepr - module procedure netcdfio_values_1d_realpr - module procedure netcdfio_values_1d_intepr - end interface ! interface netcdfio_values_2d - interface netcdfio_values_2d - module procedure netcdfio_values_2d_dblepr - module procedure netcdfio_values_2d_realpr - module procedure netcdfio_values_2d_intepr - end interface ! interface netcdfio_values_2d - interface netcdfio_values_3d - module procedure netcdfio_values_3d_dblepr - module procedure netcdfio_values_3d_realpr - module procedure netcdfio_values_3d_intepr - end interface ! interface netcdfio_values_3d - interface netcdfio_global_attr - module procedure netcdfio_global_attr_char - end interface ! interface netcdfio_global_attr - interface netcdfio_variable_attr - module procedure netcdfio_variable_attr_char - end interface ! interface netcdfio_variable_attr - public :: netcdfio_values_1d - public :: netcdfio_values_2d - public :: netcdfio_values_3d - public :: netcdfio_dimension - public :: netcdfio_global_attr - public :: netcdfio_variable_attr - public :: ncrec - public :: ncxdim - public :: ncydim - public :: nczdim - public :: nctdim - public :: ncstatic - - !----------------------------------------------------------------------- - -contains - - !======================================================================= - - ! netcdfio_global_attr.f90: - - !----------------------------------------------------------------------- - - subroutine netcdfio_global_attr_char(filename,varname,varvalue) - - ! Define variables passed to routine - - character(len=500) :: filename - character(len=*) :: varname - character(len=*) :: varvalue - - !===================================================================== - - ! Define local variables - - ncstatus = nf90_open(path=trim(adjustl(filename)),mode=nf90_nowrite, & - & ncid=ncfileid) - ncstatus = nf90_get_att(ncfileid,nf90_global,trim(adjustl(varname)), & - & varvalue) - ncstatus = nf90_close(ncfileid) - - !===================================================================== - - end subroutine netcdfio_global_attr_char - - subroutine netcdfio_variable_attr_char(filename,varname,attribute,varvalue) - - implicit none - - !======================================================================= - - ! Define variables passed to subroutine - - character(len=500),intent(in) :: filename - character(len=*),intent(in) :: attribute - character(len=*),intent(in) :: varname - - ! Define variables returned by subroutine - - character(len=80),intent(out) :: varvalue - - ! Define variables for decoding netCDF data - - integer ncid, varid, ncstatus - - ncstatus = nf90_open(path=trim(adjustl(filename)),mode=nf90_nowrite,ncid=ncid) - ncstatus = nf90_inq_varid(ncid,trim(adjustl(varname)),varid) - ncstatus = nf90_get_att(ncid,varid,trim(adjustl(attribute)),varvalue) - ncstatus = nf90_close(ncfileid) - - !===================================================================== - - end subroutine netcdfio_variable_attr_char - - !======================================================================= - - ! netcdfio_values_1d_dblepr.f90: - - !----------------------------------------------------------------------- - - subroutine netcdfio_values_1d_dblepr(filename,varname,varvalue) - - ! Define variables passed to routine - - character(len=500) :: filename - character(len=*) :: varname - real(r_double) :: varvalue(:) - - !===================================================================== - - ! Define local variables - - ncstatus = nf90_open(path=trim(adjustl(filename)),mode=nf90_nowrite, & - & ncid=ncfileid) - ncstatus = nf90_inq_varid(ncfileid,trim(adjustl(varname)),ncvarid) - if (ncstatus /= 0) then - varvalue = -1.e30 - else - ncstatus = nf90_get_var(ncfileid,ncvarid,varvalue) - if (ncstatus .ne. 0) then - print *,'fv3 read failed for ',trim(adjustl(varname)) - call mpi_interface_terminate() - stop - endif - endif - ncstatus = nf90_close(ncfileid) - - !===================================================================== - - end subroutine netcdfio_values_1d_dblepr - - !======================================================================= - - ! netcdfio_values_2d_dblepr.f90: - - !----------------------------------------------------------------------- - - subroutine netcdfio_values_2d_dblepr(filename,varname,varvalue) - - ! Define variables passed to routine - - character(len=500) :: filename - character(len=*) :: varname - real(r_double), dimension(ncxdim,ncydim) :: varvalue - - ! Define variables computed within routine - - integer, dimension(3) :: start - integer, dimension(3) :: count - - !===================================================================== - - ! Define local variables - - ncstatus = nf90_open(path=trim(adjustl(filename)),mode=nf90_nowrite, & - & ncid=ncfileid) - ncstatus = nf90_inq_varid(ncfileid,trim(adjustl(varname)),ncvarid) - if (ncstatus .ne. 0) then - print *,'fv3 read failed for ',trim(adjustl(varname)) - call mpi_interface_terminate() - stop - endif - if(ncstatic) start = (/1,1,1/) - if(.not. ncstatic) start = (/1,1,ncrec/) - count = (/ncxdim,ncydim,1/) - ncstatus = nf90_get_var(ncfileid,ncvarid,varvalue,start,count) - if (ncstatus .ne. 0) then - print *,'fv3 read failed for ',trim(adjustl(varname)) - call mpi_interface_terminate() - stop - endif - if(debug) write(6,500) trim(adjustl(varname)), minval(varvalue), & - & maxval(varvalue) - ncstatus = nf90_close(ncfileid) - - !===================================================================== - - ! Define format statements - -500 format('NETCDFIO_VALUES_2D: Variable name = ', a, '; (min,max) = (', & - & f13.5,',',f13.5,').') - - !===================================================================== - - end subroutine netcdfio_values_2d_dblepr - - !======================================================================= - - ! netcdfio_values_3d_dblepr.f90: - - !----------------------------------------------------------------------- - - subroutine netcdfio_values_3d_dblepr(filename,varname,varvalue) - - ! Define variables passed to routine - - character(len=500) :: filename - character(len=*) :: varname - real(r_double), dimension(ncxdim,ncydim,nczdim) :: varvalue - - ! Define variables computed within routine - - integer, dimension(4) :: start - integer, dimension(4) :: count - - !===================================================================== - - ! Define local variables - - ncstatus = nf90_open(path=trim(adjustl(filename)),mode=nf90_nowrite, & - & ncid=ncfileid) - ncstatus = nf90_inq_varid(ncfileid,trim(adjustl(varname)),ncvarid) - if (ncstatus .ne. 0) then - print *,'fv3 read failed for ',trim(adjustl(varname)) - call mpi_interface_terminate() - stop - endif - if(ncstatic) start = (/1,1,1,1/) - if(.not. ncstatic) start = (/1,1,1,ncrec/) - count = (/ncxdim,ncydim,nczdim,1/) - ncstatus = nf90_get_var(ncfileid,ncvarid,varvalue,start,count) - if (ncstatus .ne. 0) then - print *,'fv3 read failed for ',trim(adjustl(varname)) - call mpi_interface_terminate() - stop - endif - if(debug) write(6,500) trim(adjustl(varname)), minval(varvalue), & - & maxval(varvalue) - ncstatus = nf90_close(ncfileid) - - !===================================================================== - - ! Define format statements - -500 format('NETCDFIO_VALUES_3D: Variable name = ', a, '; (min,max) = (', & - & f13.5,',',f13.5,').') - - !===================================================================== - - end subroutine netcdfio_values_3d_dblepr - - !======================================================================= - - ! netcdfio_values_1d_realpr.f90: - - !----------------------------------------------------------------------- - - subroutine netcdfio_values_1d_realpr(filename,varname,varvalue) - - ! Define variables passed to routine - - character(len=500) :: filename - character(len=*) :: varname - real(r_kind) :: varvalue(:) - - !===================================================================== - - ! Define local variables - - ncstatus = nf90_open(path=trim(adjustl(filename)),mode=nf90_nowrite, & - & ncid=ncfileid) - ncstatus = nf90_inq_varid(ncfileid,trim(adjustl(varname)),ncvarid) - if (ncstatus .ne. 0) then - print *,'fv3 read failed for ',trim(adjustl(varname)) - call mpi_interface_terminate() - stop - endif - if (ncstatus /= 0) then - varvalue = -1.e30 - else - ncstatus = nf90_get_var(ncfileid,ncvarid,varvalue) - if (ncstatus .ne. 0) then - print *,'fv3 read failed for ',trim(adjustl(varname)) - call mpi_interface_terminate() - stop - endif - endif - ncstatus = nf90_close(ncfileid) - - !===================================================================== - - end subroutine netcdfio_values_1d_realpr - - !======================================================================= - - ! netcdfio_values_2d_realpr.f90: - - !----------------------------------------------------------------------- - - subroutine netcdfio_values_2d_realpr(filename,varname,varvalue) - - ! Define variables passed to routine - - character(len=500) :: filename - character(len=*) :: varname - real(r_kind), dimension(ncxdim,ncydim) :: varvalue - - ! Define variables computed within routine - - integer, dimension(3) :: start - integer, dimension(3) :: count - - !===================================================================== - - ! Define local variables - - ncstatus = nf90_open(path=trim(adjustl(filename)),mode=nf90_nowrite, & - & ncid=ncfileid) - ncstatus = nf90_inq_varid(ncfileid,trim(adjustl(varname)),ncvarid) - if (ncstatus .ne. 0) then - print *,'fv3 read failed for ',trim(adjustl(varname)) - call mpi_interface_terminate() - stop - endif - if(ncstatic) start = (/1,1,1/) - if(.not. ncstatic) start = (/1,1,ncrec/) - count = (/ncxdim,ncydim,1/) - ncstatus = nf90_get_var(ncfileid,ncvarid,varvalue,start,count) - if (ncstatus .ne. 0) then - print *,'fv3 read failed for ',trim(adjustl(varname)) - call mpi_interface_terminate() - stop - endif - if(debug) write(6,500) trim(adjustl(varname)), minval(varvalue), & - & maxval(varvalue) - ncstatus = nf90_close(ncfileid) - - !===================================================================== - - ! Define format statements - -500 format('NETCDFIO_VALUES_2D: Variable name = ', a, '; (min,max) = (', & - & f13.5,',',f13.5,').') - - !===================================================================== - - end subroutine netcdfio_values_2d_realpr - - !======================================================================= - - ! netcdfio_values_3d_realpr.f90: - - !----------------------------------------------------------------------- - - subroutine netcdfio_values_3d_realpr(filename,varname,varvalue) - - ! Define variables passed to routine - - character(len=500) :: filename - character(len=*) :: varname - real(r_kind), dimension(ncxdim,ncydim,nczdim) :: varvalue - - ! Define variables computed within routine - - integer, dimension(4) :: start - integer, dimension(4) :: count - - !===================================================================== - - ! Define local variables - - ncstatus = nf90_open(path=trim(adjustl(filename)),mode=nf90_nowrite, & - & ncid=ncfileid) - ncstatus = nf90_inq_varid(ncfileid,trim(adjustl(varname)),ncvarid) - if (ncstatus .ne. 0) then - print *,'fv3 read failed for ',trim(adjustl(varname)) - call mpi_interface_terminate() - stop - endif - if(ncstatic) start = (/1,1,1,1/) - if(.not. ncstatic) start = (/1,1,1,ncrec/) - count = (/ncxdim,ncydim,nczdim,1/) - ncstatus = nf90_get_var(ncfileid,ncvarid,varvalue,start,count) - if (ncstatus .ne. 0) then - print *,'fv3 read failed for ',trim(adjustl(varname)) - call mpi_interface_terminate() - stop - endif - if(debug) write(6,500) trim(adjustl(varname)), minval(varvalue), & - & maxval(varvalue) - ncstatus = nf90_close(ncfileid) - - !===================================================================== - - ! Define format statements - -500 format('NETCDFIO_VALUES_3D: Variable name = ', a, '; (min,max) = (', & - & f13.5,',',f13.5,').') - - !===================================================================== - - end subroutine netcdfio_values_3d_realpr - - !======================================================================= - - ! netcdfio_values_1d_intepr.f90: - - !----------------------------------------------------------------------- - - subroutine netcdfio_values_1d_intepr(filename,varname,varvalue) - - ! Define variables passed to routine - - character(len=500) :: filename - character(len=*) :: varname - integer :: varvalue(:) - - !===================================================================== - - ! Define local variables - - ncstatus = nf90_open(path=trim(adjustl(filename)),mode=nf90_nowrite, & - & ncid=ncfileid) - ncstatus = nf90_inq_varid(ncfileid,trim(adjustl(varname)),ncvarid) - if (ncstatus /= 0) then - varvalue = -9999 - else - ncstatus = nf90_get_var(ncfileid,ncvarid,varvalue) - endif - ncstatus = nf90_close(ncfileid) - - !===================================================================== - - end subroutine netcdfio_values_1d_intepr - - !======================================================================= - - ! netcdfio_values_2d_intepr.f90: - - !----------------------------------------------------------------------- - - subroutine netcdfio_values_2d_intepr(filename,varname,varvalue) - - ! Define variables passed to routine - - character(len=500) :: filename - character(len=*) :: varname - integer, dimension(ncxdim,ncydim) :: varvalue - - ! Define variables computed within routine - - integer, dimension(3) :: start - integer, dimension(3) :: count - - !===================================================================== - - ! Define local variables - - ncstatus = nf90_open(path=trim(adjustl(filename)),mode=nf90_nowrite, & - & ncid=ncfileid) - ncstatus = nf90_inq_varid(ncfileid,trim(adjustl(varname)),ncvarid) - if(ncstatic) start = (/1,1,1/) - if(.not. ncstatic) start = (/1,1,ncrec/) - count = (/ncxdim,ncydim,1/) - ncstatus = nf90_get_var(ncfileid,ncvarid,varvalue,start,count) - if(debug) write(6,500) trim(adjustl(varname)), minval(varvalue), & - & maxval(varvalue) - ncstatus = nf90_close(ncfileid) - - !===================================================================== - - ! Define format statements - -500 format('NETCDFIO_VALUES_2D: Variable name = ', a, '; (min,max) = (', & - & f13.5,',',f13.5,').') - - !===================================================================== - - end subroutine netcdfio_values_2d_intepr - - !======================================================================= - - ! netcdfio_values_3d_intepr.f90: - - !----------------------------------------------------------------------- - - subroutine netcdfio_values_3d_intepr(filename,varname,varvalue) - - ! Define variables passed to routine - - character(len=500) :: filename - character(len=*) :: varname - integer, dimension(ncxdim,ncydim,nczdim) :: varvalue - - ! Define variables computed within routine - - integer, dimension(4) :: start - integer, dimension(4) :: count - - !===================================================================== - - ! Define local variables - - ncstatus = nf90_open(path=trim(adjustl(filename)),mode=nf90_nowrite, & - & ncid=ncfileid) - ncstatus = nf90_inq_varid(ncfileid,trim(adjustl(varname)),ncvarid) - if (ncstatus .ne. 0) then - print *,'fv3 read failed for ',trim(adjustl(varname)) - call mpi_interface_terminate() - stop - endif - if(ncstatic) start = (/1,1,1,1/) - if(.not. ncstatic) start = (/1,1,1,ncrec/) - count = (/ncxdim,ncydim,nczdim,1/) - ncstatus = nf90_get_var(ncfileid,ncvarid,varvalue,start,count) - if(debug) write(6,500) trim(adjustl(varname)), minval(varvalue), & - & maxval(varvalue) - ncstatus = nf90_close(ncfileid) - - !===================================================================== - - ! Define format statements - -500 format('NETCDFIO_VALUES_3D: Variable name = ', a, '; (min,max) = (', & - & f13.5,',',f13.5,').') - - !===================================================================== - - end subroutine netcdfio_values_3d_intepr - - !======================================================================= - - ! netcdfio_dimension.f90: - - !----------------------------------------------------------------------- - - subroutine netcdfio_dimension(filename,dimname,dimvalue) - - ! Define variables passed to routine - - character(len=500) :: filename - character(len=*) :: dimname - integer :: dimvalue - - !===================================================================== - - ! Define local variables - - ncstatus = nf90_open(path=trim(adjustl(filename)),mode=nf90_nowrite, & - & ncid=ncfileid) - ncstatus = nf90_inq_dimid(ncfileid,trim(adjustl(dimname)),ncdimid) - ncstatus = nf90_inquire_dimension(ncfileid,ncdimid,len=dimvalue) - ncstatus = nf90_close(ncfileid) - - !===================================================================== - - end subroutine netcdfio_dimension - - !======================================================================= - -end module netcdfio_interface diff --git a/src/regrid_nemsio.fd/physcons.f90 b/src/regrid_nemsio.fd/physcons.f90 deleted file mode 100644 index 4e69dca3..00000000 --- a/src/regrid_nemsio.fd/physcons.f90 +++ /dev/null @@ -1,77 +0,0 @@ -! this module contains some the most frequently used math and ! -! physics constatns for gcm models. ! -! ! -! references: ! -! as set in NMC handbook from Smithsonian tables. ! -! ! - module physcons -! - use kinds, only : r_kind -! - implicit none -! - public - -! --- ... Math constants - - real(r_kind),parameter:: con_pi =3.1415926535897931 ! pi - real(r_kind),parameter:: con_sqrt2 =1.414214e+0 ! square root of 2 - real(r_kind),parameter:: con_sqrt3 =1.732051e+0 ! square root of 3 - -! --- ... Geophysics/Astronomy constants - - real(r_kind),parameter:: con_rerth =6.3712e+6 ! radius of earth (m) - real(r_kind),parameter:: con_g =9.80665e+0 ! gravity (m/s2) - real(r_kind),parameter:: con_omega =7.2921e-5 ! ang vel of earth (1/s) - real(r_kind),parameter:: con_p0 =1.01325e5 ! std atms pressure (pa) - real(r_kind),parameter:: con_solr =1.3660e+3 ! solar constant (W/m2)-liu(2002) - -! --- ... Thermodynamics constants - - real(r_kind),parameter:: con_rgas =8.314472 ! molar gas constant (J/mol/K) - real(r_kind),parameter:: con_rd =2.8705e+2 ! gas constant air (J/kg/K) - real(r_kind),parameter:: con_rv =4.6150e+2 ! gas constant H2O (J/kg/K) - real(r_kind),parameter:: con_cp =1.0046e+3 ! spec heat air @p (J/kg/K) - real(r_kind),parameter:: con_cv =7.1760e+2 ! spec heat air @v (J/kg/K) - real(r_kind),parameter:: con_cvap =1.8460e+3 ! spec heat H2O gas (J/kg/K) - real(r_kind),parameter:: con_cliq =4.1855e+3 ! spec heat H2O liq (J/kg/K) - real(r_kind),parameter:: con_csol =2.1060e+3 ! spec heat H2O ice (J/kg/K) - real(r_kind),parameter:: con_hvap =2.5000e+6 ! lat heat H2O cond (J/kg) - real(r_kind),parameter:: con_hfus =3.3358e+5 ! lat heat H2O fusion (J/kg) - real(r_kind),parameter:: con_psat =6.1078e+2 ! pres at H2O 3pt (Pa) - real(r_kind),parameter:: con_t0c =2.7315e+2 ! temp at 0C (K) - real(r_kind),parameter:: con_ttp =2.7316e+2 ! temp at H2O 3pt (K) - real(r_kind),parameter:: con_tice =2.7120e+2 ! temp freezing sea (K) - real(r_kind),parameter:: con_jcal =4.1855E+0 ! joules per calorie () - real(r_kind),parameter:: con_rhw0 =1022.0 ! sea water reference density (kg/m^3) - real(r_kind),parameter:: con_epsq =1.0E-12 ! min q for computing precip type - -! Secondary constants - - real(r_kind),parameter:: con_rocp =con_rd/con_cp - real(r_kind),parameter:: con_cpor =con_cp/con_rd - real(r_kind),parameter:: con_rog =con_rd/con_g - real(r_kind),parameter:: con_fvirt =con_rv/con_rd-1. - real(r_kind),parameter:: con_eps =con_rd/con_rv - real(r_kind),parameter:: con_epsm1 =con_rd/con_rv-1. - real(r_kind),parameter:: con_dldt =con_cvap-con_cliq - real(r_kind),parameter:: con_xpona =-con_dldt/con_rv - real(r_kind),parameter:: con_xponb =-con_dldt/con_rv+con_hvap/(con_rv*con_ttp) - -! --- ... Other Physics/Chemistry constants (source: 2002 CODATA) - - real(r_kind),parameter:: con_c =2.99792458e+8 ! speed of light (m/s) - real(r_kind),parameter:: con_plnk =6.6260693e-34 ! planck constatn (J/s) - real(r_kind),parameter:: con_boltz =1.3806505e-23 ! boltzmann constant (J/K) - real(r_kind),parameter:: con_sbc =5.670400e-8 ! stefan-boltzmann (W/m2/K4) - real(r_kind),parameter:: con_avgd =6.0221415e23 ! avogadro constant (1/mol) - real(r_kind),parameter:: con_gasv =22413.996e-6 ! vol of ideal gas at 273.15k, 101.325kpa (m3/mol) - real(r_kind),parameter:: con_amd =28.9644 ! molecular wght of dry air (g/mol) - real(r_kind),parameter:: con_amw =18.0154 ! molecular wght of water vapor (g/mol) - real(r_kind),parameter:: con_amo3 =47.9982 ! molecular wght of o3 (g/mol) - real(r_kind),parameter:: con_amco2 =44.011 ! molecular wght of co2 (g/mol) - real(r_kind),parameter:: con_amo2 =31.9999 ! molecular wght of o2 (g/mol) - real(r_kind),parameter:: con_amch4 =16.043 ! molecular wght of ch4 (g/mol) - real(r_kind),parameter:: con_amn2o =44.013 ! molecular wght of n2o (g/mol) - -end module physcons diff --git a/src/regrid_nemsio.fd/regrid_nemsio_interface.f90 b/src/regrid_nemsio.fd/regrid_nemsio_interface.f90 deleted file mode 100644 index 9ab5597a..00000000 --- a/src/regrid_nemsio.fd/regrid_nemsio_interface.f90 +++ /dev/null @@ -1,50 +0,0 @@ -module regrid_nemsio_interface - - !======================================================================= - - ! Define associated modules and subroutines - - !----------------------------------------------------------------------- - - use constants - use kinds - - !----------------------------------------------------------------------- - - use fv3_interface - use gfs_nems_interface - use namelist_def - - !----------------------------------------------------------------------- - - implicit none - - !----------------------------------------------------------------------- - -contains - - !======================================================================= - - ! regrid_nemsio.f90: - - !----------------------------------------------------------------------- - - subroutine regrid_nemsio() - - !===================================================================== - - ! Define local variables - - call namelistparams() - - ! Check local variable and proceed accordingly - - call fv3_regrid_nemsio() - - !===================================================================== - - end subroutine regrid_nemsio - - !======================================================================= - -end module regrid_nemsio_interface diff --git a/src/regrid_nemsio.fd/variable_interface.f90 b/src/regrid_nemsio.fd/variable_interface.f90 deleted file mode 100644 index d0d56842..00000000 --- a/src/regrid_nemsio.fd/variable_interface.f90 +++ /dev/null @@ -1,66 +0,0 @@ -module variable_interface - - !======================================================================= - - ! Define associated modules and subroutines - - !----------------------------------------------------------------------- - - use kinds - use physcons, only: rgas => con_rd, cp => con_cp, grav => con_g, & - & rerth => con_rerth, rocp => con_rocp, & - & pi => con_pi, con_rog - - !----------------------------------------------------------------------- - - use mpi_interface - use namelist_def - - !----------------------------------------------------------------------- - - implicit none - - !----------------------------------------------------------------------- - - ! Define interfaces and attributes for module routines - - private - public :: varinfo - !public :: variable_lookup - public :: variable_clip - - !----------------------------------------------------------------------- - - ! Define all data and structure types for routine; these variables - ! are variables required by the subroutines within this module - - type varinfo - character(len=20) :: var_name - character(len=20) :: nems_name - character(len=20) :: nems_levtyp - integer :: nems_lev - character(len=20) :: itrptyp - logical :: clip - integer :: ndims - end type varinfo ! type varinfo - - !----------------------------------------------------------------------- - -contains - - !======================================================================= - - subroutine variable_clip(grid) - - - real(r_double) :: grid(:) - real(r_double) :: clip - - clip = tiny(grid(1)) - where(grid .le. dble(0.0)) grid = clip - - end subroutine variable_clip - - !======================================================================= - -end module variable_interface