diff --git a/CMakeLists.txt b/CMakeLists.txt index 483d8defd..d9bfe482e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,7 +7,7 @@ file(STRINGS "VERSION" pVersion LIMIT_COUNT 1) project( upp VERSION ${pVersion} - LANGUAGES Fortran) + LANGUAGES Fortran C CXX) list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake") @@ -15,6 +15,8 @@ list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake") option(OPENMP "use OpenMP threading" ON) option(BUILD_POSTEXEC "Build NCEPpost executable" ON) option(BUILD_WITH_WRFIO "Build NCEPpost with WRF-IO library" OFF) +option(BUILD_WITH_IFI "Build NCEPpost with In-Flight Icing (IFI) library if present" OFF) +option(REQUIRE_IFI "Abort if libIFI is not found ; enables BUILD_WITH_IFI=ON" OFF) option(BUILD_WITH_GTG "Build NCEPpost with NCAR/GTG" OFF) option(ENABLE_DOCS "Enable generation of doxygen-based documentation." OFF) @@ -52,6 +54,18 @@ if(BUILD_WITH_GTG) find_package(ip REQUIRED) endif() +if(REQUIRE_IFI) + set(BUILD_WITH_IFI ON) +endif() + +if(BUILD_WITH_IFI) + if(REQUIRE_IFI) + find_package(IFI REQUIRED) + else() + find_package(IFI) + endif() +endif() + if(BUILD_POSTEXEC) find_package(nemsio REQUIRED) find_package(sfcio REQUIRED) @@ -61,6 +75,22 @@ if(BUILD_POSTEXEC) if(BUILD_WITH_WRFIO) find_package(wrf_io REQUIRED) endif() + if(IFI_FOUND) + if(CMAKE_Fortran_COMPILER_ID MATCHES "^Intel$") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -cxxlib") + elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^GNU$") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -lstdc++") + elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU|AppleClang)$") + # This one is a wild guess. I haven't tried linking a Fortran + # executable with C++ libraries in clang. + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -stdlib=libc++ -lc++abi") + endif() + endif() +endif() + +if(IFI_FOUND) + message(NOTICE "Enabling USE_IFI in the post.") + add_definitions(-DUSE_IFI=1) endif() add_subdirectory(sorc) diff --git a/modulefiles/hera.lua b/modulefiles/hera.lua index a0e6574a7..adc9b001c 100644 --- a/modulefiles/hera.lua +++ b/modulefiles/hera.lua @@ -51,8 +51,15 @@ load(pathJoin("sfcio", sfcio_ver)) wrf_io_ver=os.getenv("wrf_io_ver") or "1.2.0" load(pathJoin("wrf_io", wrf_io_ver)) +ifi_ver=os.getenv("ifi_ver") or "20220901-intel-2022.1.2" +prepend_path("MODULEPATH", "/scratch2/BMC/ifi/modulefiles") +try_load(pathJoin("ifi",ifi_ver)) + setenv("CC","mpiicc") setenv("CXX","mpiicpc") setenv("FC","mpiifort") +prepend_path("MODULEPATH", "/scratch2/BMC/ifi/modulefiles") +try_load("ifi/20220901-intel-2022.1.2") + whatis("Description: post build environment") diff --git a/modulefiles/jet b/modulefiles/jet index e789ff698..39695fcc0 100644 --- a/modulefiles/jet +++ b/modulefiles/jet @@ -32,3 +32,6 @@ module load sigio/2.3.2 module load sp/2.3.3 module load w3emc/2.9.2 module load wrf_io/1.1.1 + +module use /lfs4/BMC/ifi/modulefiles +module try-load ifi/20220906-intel-18.0.5.274 diff --git a/parm/fv3lam_post_avblflds.xml b/parm/fv3lam_post_avblflds.xml index 4f9b7b4ab..cbb28ecf9 100644 --- a/parm/fv3lam_post_avblflds.xml +++ b/parm/fv3lam_post_avblflds.xml @@ -5678,6 +5678,41 @@ 4.0 + + 1003 + ICE_PROB_IFI_FLIGHT_LEVEL + ICPRB + NCEP + spec_alt_above_mean_sea_lvl + 4.0 + + 1 + 5000. 10000. 15000. 20000. 25000. 30000. 35000. 40000. 45000. 50000. 55000. 60000. 65000. 70000. 75000. 80000. 85000. 90000. 95000. 100000. 105000. 110000. 115000. 120000. 125000. 130000. 135000. 140000. 145000. 150000. 155000. 160000. 165000. 170000. 175000. 180000. 185000. 190000. 195000. 200000. 205000. 210000. 215000. 220000. 225000. 230000. 235000. 240000. 245000. 250000. 255000. 260000. 265000. 270000. 275000. 280000. 285000. 290000. 295000. 300000. + + + + 1004 + SLD_IFI_FLIGHT_LEVEL + SIPD + NCEP + spec_alt_above_mean_sea_lvl + 4.0 + + 1 + 5000. 10000. 15000. 20000. 25000. 30000. 35000. 40000. 45000. 50000. 55000. 60000. 65000. 70000. 75000. 80000. 85000. 90000. 95000. 100000. 105000. 110000. 115000. 120000. 125000. 130000. 135000. 140000. 145000. 150000. 155000. 160000. 165000. 170000. 175000. 180000. 185000. 190000. 195000. 200000. 205000. 210000. 215000. 220000. 225000. 230000. 235000. 240000. 245000. 250000. 255000. 260000. 265000. 270000. 275000. 280000. 285000. 290000. 295000. 300000. + + + + 1005 + ICE_SEV_CAT_IFI_FLIGHT_LEVEL + ICSEV + NCEP + spec_alt_above_mean_sea_lvl + 4.0 + + 1 + 5000. 10000. 15000. 20000. 25000. 30000. 35000. 40000. 45000. 50000. 55000. 60000. 65000. 70000. 75000. 80000. 85000. 90000. 95000. 100000. 105000. 110000. 115000. 120000. 125000. 130000. 135000. 140000. 145000. 150000. 155000. 160000. 165000. 170000. 175000. 180000. 185000. 190000. 195000. 200000. 205000. 210000. 215000. 220000. 225000. 230000. 235000. 240000. 245000. 250000. 255000. 260000. 265000. 270000. 275000. 280000. 285000. 290000. 295000. 300000. + diff --git a/parm/fv3lam_rrfs.xml b/parm/fv3lam_rrfs.xml index 68db26703..79ca309d3 100755 --- a/parm/fv3lam_rrfs.xml +++ b/parm/fv3lam_rrfs.xml @@ -3539,4 +3539,49 @@ + + + IFIFIP + 4 + ncep_emc + v2003 + local_tab_yes1 + fcst + oper + fcst + fcst + hour + nws_ncep + hrrr + complex_packing_spatial_diff + 2nd_ord_sptdiff + fltng_pnt + lossless + + + + + ICE_PROB_IFI_FLIGHT_LEVEL + 4.0 + + 1 + 5000. 10000. 15000. 20000. 25000. 30000. 35000. 40000. 45000. 50000. 55000. 60000. 65000. 70000. 75000. 80000. 85000. 90000. 95000. 100000. 105000. 110000. 115000. 120000. 125000. 130000. 135000. 140000. 145000. 150000. 155000. 160000. 165000. 170000. 175000. 180000. 185000. 190000. 195000. 200000. 205000. 210000. 215000. 220000. 225000. 230000. 235000. 240000. 245000. 250000. 255000. 260000. 265000. 270000. 275000. 280000. 285000. 290000. 295000. 300000. + + + + SLD_IFI_FLIGHT_LEVEL + 4.0 + + 1 + 5000. 10000. 15000. 20000. 25000. 30000. 35000. 40000. 45000. 50000. 55000. 60000. 65000. 70000. 75000. 80000. 85000. 90000. 95000. 100000. 105000. 110000. 115000. 120000. 125000. 130000. 135000. 140000. 145000. 150000. 155000. 160000. 165000. 170000. 175000. 180000. 185000. 190000. 195000. 200000. 205000. 210000. 215000. 220000. 225000. 230000. 235000. 240000. 245000. 250000. 255000. 260000. 265000. 270000. 275000. 280000. 285000. 290000. 295000. 300000. + + + + ICE_SEV_CAT_IFI_FLIGHT_LEVEL + 4.0 + + 1 + 5000. 10000. 15000. 20000. 25000. 30000. 35000. 40000. 45000. 50000. 55000. 60000. 65000. 70000. 75000. 80000. 85000. 90000. 95000. 100000. 105000. 110000. 115000. 120000. 125000. 130000. 135000. 140000. 145000. 150000. 155000. 160000. 165000. 170000. 175000. 180000. 185000. 190000. 195000. 200000. 205000. 210000. 215000. 220000. 225000. 230000. 235000. 240000. 245000. 250000. 255000. 260000. 265000. 270000. 275000. 280000. 285000. 290000. 295000. 300000. + + diff --git a/parm/post_avblflds.xml b/parm/post_avblflds.xml index 89065f857..fcd01f2d4 100755 --- a/parm/post_avblflds.xml +++ b/parm/post_avblflds.xml @@ -8024,5 +8024,41 @@ surface 3.0 + + + 1003 + ICE_PROB_IFI_FLIGHT_LEVEL + ICPRB + NCEP + spec_alt_above_mean_sea_lvl + 4.0 + + 1 + 5000. 10000. 15000. 20000. 25000. 30000. 35000. 40000. 45000. 50000. 55000. 60000. 65000. 70000. 75000. 80000. 85000. 90000. 95000. 100000. 105000. 110000. 115000. 120000. 125000. 130000. 135000. 140000. 145000. 150000. 155000. 160000. 165000. 170000. 175000. 180000. 185000. 190000. 195000. 200000. 205000. 210000. 215000. 220000. 225000. 230000. 235000. 240000. 245000. 250000. 255000. 260000. 265000. 270000. 275000. 280000. 285000. 290000. 295000. 300000. + + + + 1004 + SLD_IFI_FLIGHT_LEVEL + SIPD + NCEP + spec_alt_above_mean_sea_lvl + 4.0 + + 1 + 5000. 10000. 15000. 20000. 25000. 30000. 35000. 40000. 45000. 50000. 55000. 60000. 65000. 70000. 75000. 80000. 85000. 90000. 95000. 100000. 105000. 110000. 115000. 120000. 125000. 130000. 135000. 140000. 145000. 150000. 155000. 160000. 165000. 170000. 175000. 180000. 185000. 190000. 195000. 200000. 205000. 210000. 215000. 220000. 225000. 230000. 235000. 240000. 245000. 250000. 255000. 260000. 265000. 270000. 275000. 280000. 285000. 290000. 295000. 300000. + + + + 1005 + ICE_SEV_CAT_IFI_FLIGHT_LEVEL + ICSEV + NCEP + spec_alt_above_mean_sea_lvl + 4.0 + + 1 + 5000. 10000. 15000. 20000. 25000. 30000. 35000. 40000. 45000. 50000. 55000. 60000. 65000. 70000. 75000. 80000. 85000. 90000. 95000. 100000. 105000. 110000. 115000. 120000. 125000. 130000. 135000. 140000. 145000. 150000. 155000. 160000. 165000. 170000. 175000. 180000. 185000. 190000. 195000. 200000. 205000. 210000. 215000. 220000. 225000. 230000. 235000. 240000. 245000. 250000. 255000. 260000. 265000. 270000. 275000. 280000. 285000. 290000. 295000. 300000. + diff --git a/parm/postxconfig-NT-fv3lam_rrfs.txt b/parm/postxconfig-NT-fv3lam_rrfs.txt index 7c5b29e33..185cc746f 100644 --- a/parm/postxconfig-NT-fv3lam_rrfs.txt +++ b/parm/postxconfig-NT-fv3lam_rrfs.txt @@ -1,4 +1,5 @@ -2 +3 +3 239 276 PRSLEV @@ -19088,3 +19089,130 @@ top_of_atmos ? ? ? +IFIFIP +4 +ncep_emc +v2003 +local_tab_yes1 +fcst +oper +fcst +fcst +hour +nws_ncep +hrrr +complex_packing_spatial_diff +2nd_ord_sptdiff +fltng_pnt +lossless +1003 +ICE_PROB_IFI_FLIGHT_LEVEL +? +1 +tmpl4_0 +ICPRB +NCEP +? +spec_alt_above_mean_sea_lvl +1 +1 +60 +5000. 10000. 15000. 20000. 25000. 30000. 35000. 40000. 45000. 50000. 55000. 60000. 65000. 70000. 75000. 80000. 85000. 90000. 95000. 100000. 105000. 110000. 115000. 120000. 125000. 130000. 135000. 140000. 145000. 150000. 155000. 160000. 165000. 170000. 175000. 180000. 185000. 190000. 195000. 200000. 205000. 210000. 215000. 220000. 225000. 230000. 235000. 240000. 245000. 250000. 255000. 260000. 265000. 270000. 275000. 280000. 285000. 290000. 295000. 300000. +? +0 +? +0 +? +? +? +0 +0.0 +0 +0.0 +? +0 +0.0 +0 +0.0 +1 +4.0 +0 +0 +0 +? +? +? +1004 +SLD_IFI_FLIGHT_LEVEL +? +1 +tmpl4_0 +SIPD +NCEP +? +spec_alt_above_mean_sea_lvl +1 +1 +60 +5000. 10000. 15000. 20000. 25000. 30000. 35000. 40000. 45000. 50000. 55000. 60000. 65000. 70000. 75000. 80000. 85000. 90000. 95000. 100000. 105000. 110000. 115000. 120000. 125000. 130000. 135000. 140000. 145000. 150000. 155000. 160000. 165000. 170000. 175000. 180000. 185000. 190000. 195000. 200000. 205000. 210000. 215000. 220000. 225000. 230000. 235000. 240000. 245000. 250000. 255000. 260000. 265000. 270000. 275000. 280000. 285000. 290000. 295000. 300000. +? +0 +? +0 +? +? +? +0 +0.0 +0 +0.0 +? +0 +0.0 +0 +0.0 +1 +4.0 +0 +0 +0 +? +? +? +1005 +ICE_SEV_CAT_IFI_FLIGHT_LEVEL +? +1 +tmpl4_0 +ICSEV +NCEP +? +spec_alt_above_mean_sea_lvl +1 +1 +60 +5000. 10000. 15000. 20000. 25000. 30000. 35000. 40000. 45000. 50000. 55000. 60000. 65000. 70000. 75000. 80000. 85000. 90000. 95000. 100000. 105000. 110000. 115000. 120000. 125000. 130000. 135000. 140000. 145000. 150000. 155000. 160000. 165000. 170000. 175000. 180000. 185000. 190000. 195000. 200000. 205000. 210000. 215000. 220000. 225000. 230000. 235000. 240000. 245000. 250000. 255000. 260000. 265000. 270000. 275000. 280000. 285000. 290000. 295000. 300000. +? +0 +? +0 +? +? +? +0 +0.0 +0 +0.0 +? +0 +0.0 +0 +0.0 +1 +4.0 +0 +0 +0 +? +? +? diff --git a/sorc/ncep_post.fd/ALLOCATE_ALL.f b/sorc/ncep_post.fd/ALLOCATE_ALL.f index a93741328..7773bdc5b 100644 --- a/sorc/ncep_post.fd/ALLOCATE_ALL.f +++ b/sorc/ncep_post.fd/ALLOCATE_ALL.f @@ -30,6 +30,7 @@ !! SUBROUTINE ALLOCATE_ALL() ! + use upp_ifi_mod, only: set_ifi_dims use vrbls4d use vrbls3d use vrbls2d @@ -76,6 +77,8 @@ SUBROUTINE ALLOCATE_ALL() allocate(tcucn(ista_2l:iend_2u,jsta_2l:jend_2u,lm)) allocate(EL_PBL(ista_2l:iend_2u,jsta_2l:jend_2u,lm)) + call set_ifi_dims() ! set ifi_nflight and ifi_flight_levels + !Initialization !$omp parallel do private(i,j,l) do l=1,lm @@ -115,7 +118,7 @@ SUBROUTINE ALLOCATE_ALL() exch_h(i,j,l)=spval train(i,j,l)=spval tcucn(i,j,l)=spval - EL_PBL(i,j,l)=spval + EL_PBL(i,j,l)=spval enddo enddo enddo @@ -335,6 +338,9 @@ SUBROUTINE ALLOCATE_ALL() ! ! FROM VRBLS2D ! + allocate(CAPE(ista_2l:iend_2u,jsta_2l:jend_2u)) + allocate(CIN(ista_2l:iend_2u,jsta_2l:jend_2u)) + allocate(IFI_APCP(ista_2l:iend_2u,jsta_2l:jend_2u)) ! SRD allocate(wspd10max(ista_2l:iend_2u,jsta_2l:jend_2u)) allocate(w_up_max(ista_2l:iend_2u,jsta_2l:jend_2u)) @@ -355,6 +361,9 @@ SUBROUTINE ALLOCATE_ALL() !$omp parallel do private(i,j) do j=jsta_2l,jend_2u do i=ista_2l,iend_2u + CAPE(i,j)=spval + CIN(i,j)=spval + IFI_APCP(i,j)=spval wspd10max(i,j)=spval w_up_max(i,j)=spval w_dn_max(i,j)=spval diff --git a/sorc/ncep_post.fd/CMakeLists.txt b/sorc/ncep_post.fd/CMakeLists.txt index 3bd218cf3..40721c758 100644 --- a/sorc/ncep_post.fd/CMakeLists.txt +++ b/sorc/ncep_post.fd/CMakeLists.txt @@ -46,6 +46,7 @@ list(APPEND LIB_SRC ETCALC.f ETAMP_Q2F.f EXCH.f + EXCH_c_float.f FDLVL.f FGAMMA.f FILL_PSETFLD.f @@ -62,6 +63,7 @@ list(APPEND LIB_SRC grib2_module.f GRIDSPEC.f ICAOHEIGHT.f + IFI.F kinds_mod.F LFMFLD.f LFMFLD_GFS.f @@ -169,7 +171,7 @@ if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$") set(CMAKE_Fortran_FLAGS "-g -traceback -fp-model source -free -convert big_endian") set(CMAKE_Fortran_FLAGS_RELEASE "-O3") - set(CMAKE_Fortran_FLAGS_DEBUG "-O0") + set(CMAKE_Fortran_FLAGS_DEBUG "-O0 -check bounds -check pointers -check shape -check stack -check uninit") elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU|Clang|AppleClang)$") set(CMAKE_Fortran_FLAGS "-g -fbacktrace -ffree-form -ffree-line-length-none -fconvert=big-endian") @@ -203,6 +205,11 @@ target_link_libraries(${LIBNAME} PUBLIC MPI::MPI_Fortran NetCDF::NetCDF_Fortran) +if(IFI_FOUND) + target_link_libraries(${LIBNAME} PUBLIC IFI) + target_link_libraries(${LIBNAME} PUBLIC tdrp) +endif() + if(BUILD_WITH_GTG) target_link_libraries(${LIBNAME} PUBLIC ip::ip_4) @@ -221,6 +228,10 @@ if(BUILD_POSTEXEC) sp::sp_4 sfcio::sfcio sigio::sigio) + if(IFI_FOUND) + target_link_libraries(${EXENAME} PRIVATE IFI) + target_link_libraries(${EXENAME} PRIVATE tdrp) + endif() if(wrf_io_FOUND) target_link_libraries(${EXENAME} PRIVATE wrf_io::wrf_io) diff --git a/sorc/ncep_post.fd/CTLBLK.f b/sorc/ncep_post.fd/CTLBLK.f index 5ca6a0f60..9d7442e3a 100644 --- a/sorc/ncep_post.fd/CTLBLK.f +++ b/sorc/ncep_post.fd/CTLBLK.f @@ -83,7 +83,8 @@ module CTLBLK_mod real(kind=8) :: ETAFLD2_tim=0.,ETA2P_tim=0.,SURFCE2_tim=0., & CLDRAD_tim=0.,MISCLN_tim=0.,FIXED_tim=0., & MDL2SIGMA_tim=0.,READxml_tim=0.,MDL2AGL_tim=0., & - MDL2STD_tim=0.,MDL2THANDPV_tim=0.,CALRAD_WCLOUD_tim=0.!comm tim_info + MDL2STD_tim=0.,MDL2THANDPV_tim=0., & + CALRAD_WCLOUD_tim=0.,RUN_IFI_TIM=0. !comm tim_info ! real(kind=8) :: time_output=0., time_e2out=0. !comm jjt ! @@ -114,5 +115,12 @@ module CTLBLK_mod DATA SIGBND / 0.985,0.955,0.925,0.895,0.865,0.835 / DATA PETABND / 15.,45.,75.,105.,135.,165./ ! +! Precipitation bucket time + real :: ITPREC=-1 +! +! Flight levels in feet, provided by libIFI + integer :: ifi_nflight = 0 + real, allocatable :: ifi_flight_levels(:) ! units are FEET +! !----------------------------------------------------------------------- end module CTLBLK_mod diff --git a/sorc/ncep_post.fd/DEALLOCATE.f b/sorc/ncep_post.fd/DEALLOCATE.f index ac9d5c648..ba5b57944 100644 --- a/sorc/ncep_post.fd/DEALLOCATE.f +++ b/sorc/ncep_post.fd/DEALLOCATE.f @@ -53,6 +53,7 @@ SUBROUTINE DE_ALLOCATE ! deallocate(rainw(im,jsta_2l:jend_2u,lm)) deallocate(q2) deallocate(omga) + deallocate(dpres) deallocate(T_ADJ) deallocate(ttnd) deallocate(rswtt) @@ -354,6 +355,7 @@ SUBROUTINE DE_ALLOCATE deallocate(avgetrans) deallocate(avgesnow) deallocate(avgpotevp) + deallocate(aod550) deallocate(ti) deallocate(du_aod550) deallocate(ss_aod550) @@ -500,7 +502,6 @@ SUBROUTINE DE_ALLOCATE deallocate(ssdp) deallocate(sswt) deallocate(sssv) - deallocate(dpres) deallocate(rhomid) ! vrbls2d deallocate(dusmass) diff --git a/sorc/ncep_post.fd/EXCH_c_float.f b/sorc/ncep_post.fd/EXCH_c_float.f new file mode 100644 index 000000000..8404f80f9 --- /dev/null +++ b/sorc/ncep_post.fd/EXCH_c_float.f @@ -0,0 +1,409 @@ +!> @file +!> @brief Subroutines that exchange one halo row. +!> +!> These routines are to exchange one halo row. +!> +!> @param[in] A Array to have halos exchanged. +!> @param[out] A Array with halos exchanged. +!> +!> ### Program history log: +!> Date | Programmer | Comments +!> -----------|---------------------|---------- +!> 2000-01-06 | Jim Tuccillo | Initial +!> 2021-06-01 | George Vandenberghe | 2D decomposition +!> +!> @note The 1st line is an inlined compiler directive that turns off -qcheck +!> during compilation, even if it's specified as a compiler option in the +!> makefile (Tuccillo, personal communication; Ferrier, Feb '02). +!> +!> @author Jim Tuccillo IBM @date 2000-01-06 + SUBROUTINE EXCH_c_float(A) + + use ctlblk_mod, only: num_procs, jend, iup, jsta, idn, mpi_comm_comp, im,& + icoords,ibcoords,bufs,ibufs,me,numx, & + jsta_2l, jend_2u,ileft,iright,ista_2l,iend_2u,ista,iend,jm,modelname +!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + use iso_c_binding, only: c_sizeof, c_float + use mpi + implicit none +! + real(kind=c_float),intent(inout) :: a ( ista_2l:iend_2u,jsta_2l:jend_2u ) + real(kind=c_float), allocatable :: coll(:), colr(:) + integer, allocatable :: icoll(:), icolr(:) + integer status(MPI_STATUS_SIZE) + integer ierr, jstam1, jendp1,j + integer size,ubound,lbound + integer msglenl, msglenr + integer i,ii,jj, ibl,ibu,jbl,jbu,icc,jcc + integer iwest,ieast + integer ifirst + integer mpi_kind + + logical, parameter :: checkcoords = .false. + + data ifirst/0/ + allocate(coll(jm)) + allocate(colr(jm)) + allocate(icolr(jm)) + allocate(icoll(jm)) + ibl=max(ista-1,1) + ibu=min(im,iend+1) + jbu=min(jm,jend+1) + jbl=max(jsta-1,1) +! + +! write(0,*) 'mype=',me,'num_procs=',num_procs,'im=',im,'jsta_2l=', & +! jsta_2l,'jend_2u=',jend_2u,'jend=',jend,'iup=',iup,'jsta=', & +! jsta,'idn=',idn + if ( num_procs <= 1 ) return +! +! for global model apply cyclic boundary condition + + IF(MODELNAME == 'GFS') then + if(ifirst .le. 0 .and. me .eq. 0) print *,' CYCLIC BC APPLIED' + if(ileft .eq. MPI_PROC_NULL) iwest=1 ! get eastern bc from western boundary of full domain + if(iright .eq. MPI_PROC_NULL) ieast=1 ! get western bc from eastern boundary of full domain + if(ileft .eq. MPI_PROC_NULL) ileft=me+(numx-1) + if(iright .eq. MPI_PROC_NULL) iright=(me-numx) +1 + endif + + jstam1 = max(jsta_2l,jsta-1) ! Moorthi + +! send last row to iup's first row+ and receive first row- from idn's last row + + call mpi_sendrecv(a(ista,jend),iend-ista+1,MPI_REAL4,iup,1, & + & a(ista,jstam1),iend-ista+1,MPI_REAL4,idn,1, & + & MPI_COMM_COMP,status,ierr) + + if ( ierr /= 0 ) then + print *, ' problem with first sendrecv in exch, ierr = ',ierr + stop 6661 + endif + + if (checkcoords) then + if(ifirst .le. 0) then !IFIRST ONLY + call mpi_sendrecv(ibcoords(ista,jend),iend-ista+1,MPI_INTEGER,iup,1, & + & ibcoords(ista,jstam1),iend-ista+1,MPI_INTEGER,idn,1, & + & MPI_COMM_COMP,status,ierr) + if ( ierr /= 0 ) then + print *, ' problem with second sendrecv in exch, ierr = ',ierr + stop 7661 + endif + do i=ista,iend + ii=ibcoords(i,jstam1)/10000 + jj=ibcoords(i,jstam1)-(ii*10000) + if(ii .ne. i .or. jj .ne. jstam1 ) print *,' GWVX JEXCH CHECK FAIL ',ii,jj,ibcoords(i,jstam1),i + end do + endif !IFIRST + endif !checkcoords + +! build the I columns to send and receive + + msglenl=jend-jsta+1 + msglenr=jend-jsta+1 + if(iright .lt. 0) msglenr=1 + if(ileft .lt. 0) msglenl=1 + + do j=jsta,jend + coll(j)=a(ista,j) + end do + + call mpi_barrier(mpi_comm_comp,ierr) + +! send first col to ileft last col+ and receive last col+ from ileft first col + + call mpi_sendrecv(coll(jsta),msglenl ,MPI_REAL4,ileft,1, & + & colr(jsta),msglenr ,MPI_REAL4,iright,1, & + & MPI_COMM_COMP,status,ierr) + + if ( ierr /= 0 ) then + print *, ' problem with third sendrecv in exch, ierr = ',ierr + stop 6662 + endif + + if(ifirst .le. 0) then ! IFIRST ONLY + call mpi_sendrecv(icoll(jsta),msglenl ,MPI_INTEGER,ileft,1, & + & icolr(jsta),msglenr ,MPI_INTEGER,iright,1, & + & MPI_COMM_COMP,status,ierr) + if ( ierr /= 0 ) then + print *, ' problem with fourth sendrecv in exch, ierr = ',ierr + stop 7662 + endif + endif !IFIRST + + if(iright .ge. 0) then + do j=jsta,jend + a(iend+1,j)=colr(j) + if(checkcoords) then + if(ifirst .le. 0) then !IFIRST ONLY + ibcoords(iend+1,j)=icolr(j) + ii=ibcoords(iend+1,j)/10000 + jj=ibcoords( iend+1,j)-(ii*10000) + if( j .ne. jj .or. ii .ne. iend+1 .and. ii .ne. im .and. ii .ne. 1) & + write(0,921) j,iend+1,ii,jj,ibcoords(iend+1,j),'IEXCH COORD FAIL j,iend+1,ii,jj,ibcoord ' + endif !IFIRST + endif !checkcoords + end do + endif ! for iright + + 921 format(5i10,a50) + +! print *,'mype=',me,'in EXCH, after first mpi_sendrecv' + + if ( ierr /= 0 ) then + print *, ' problem with fifth sendrecv in exch, ierr = ',ierr + stop 6663 + end if + jendp1 = min(jend+1,jend_2u) ! Moorthi + +!GWV. change from full im row exchange to iend-ista+1 subrow exchange, + + do j=jsta,jend + colr(j)=a(iend,j) + end do + +! send first row to idown's last row+ and receive last row+ from iup's first row + + call mpi_sendrecv(a(ista,jsta),iend-ista+1,MPI_REAL4,idn,1, & + & a(ista,jendp1),iend-ista+1,MPI_REAL4,iup,1, & + & MPI_COMM_COMP,status,ierr) + if ( ierr /= 0 ) then + print *, ' problem with sixth sendrecv in exch, ierr = ',ierr + stop 6664 + endif + + if (checkcoords) then + if (ifirst .le. 0) then + call mpi_sendrecv(ibcoords(ista,jsta),iend-ista+1,MPI_INTEGER,idn,1, & + & ibcoords(ista,jendp1),iend-ista+1,MPI_INTEGER,iup,1, & + & MPI_COMM_COMP,status,ierr) + if ( ierr /= 0 ) then + print *, ' problem with seventh sendrecv in exch, ierr = ',ierr + stop 7664 + endif + endif ! IFIRST + endif ! checkcoords + +! send last col to iright first col- and receive first col- from ileft last col + + call mpi_sendrecv(colr(jsta),msglenr ,MPI_REAL4,iright,1 , & + & coll(jsta),msglenl ,MPI_REAL4,ileft ,1, & + & MPI_COMM_COMP,status,ierr) + + if ( ierr /= 0 ) then + print *, ' problem with eighth sendrecv in exch, ierr = ',ierr + stop 6665 + endif + + if (ifirst .le. 0) then + call mpi_sendrecv(icolr(jsta),msglenr ,MPI_integer,iright,1 , & + & icoll(jsta),msglenl ,MPI_integer,ileft ,1, & + & MPI_COMM_COMP,status,ierr) + if ( ierr /= 0 ) then + print *, ' problem with ninth sendrecv in exch, ierr = ',ierr + stop 7665 + endif + endif !IFIRST + + if(ileft .ge. 0) then + do j=jsta,jend + a(ista-1,j)=coll(j) + if(checkcoords) then + if(ifirst .le. 0) then + ibcoords(ista-1,j)=icoll(j) + ii=ibcoords(ista-1,j)/10000 + jj=ibcoords( ista-1,j)-(ii*10000) + if( j .ne. jj .or. ii .ne. ista-1 .and. ii .ne. im .and. ii .ne. 1) & + write(0,921) j,ista-1,ii,jj,ibcoords(ista-1,j),'EXCH COORD FAIL j,ista-1,ii,jj,ibcoord ' + endif !IFIRST + endif !checkcoords + end do + endif + +! interior check + + if(checkcoords) then + if(ifirst .le. 0) then + do j=jsta,jend + do i=ista,iend + ii=ibcoords(i,j)/10000 + jj=ibcoords( i,j)-(ii*10000) + if(ii .ne. i .or. jj .ne. j) write(0,151) 'INFAILED IJ ',i,j,ibcoords(i,j),ibl,jbl,ibu,jbu + end do + end do + endif !IFIRST + endif !checkcoords + +!! corner points. After the exchanges above, corner points are replicated in +! neighbour halos so we can get them from the neighbors rather than +! calculating more corner neighbor numbers +! A(ista-1,jsta-1) is in the ileft a(iend,jsta-1) location +! A(ista-1,jend+1) is in the ileft a(iend,jend+1) location +! A(iend+1,jsta-1) is in the iright a(ista,jsta-1) location +! A(iend+1,jend+1) is in the iright a(ista,jend+1) location +!GWVx ibl=max(ista-1,1) +!GWVx ibu=min(im,iend+1) + + ibl=max(ista-1,1) + ibu=min(im,iend+1) + if(modelname == 'GFS') then + ibl=max(ista-1,0) + ibu=min(im+1,iend+1) + endif + + jbu=min(jm,jend+1) + jbl=max(jsta-1,1) + + call mpi_sendrecv(a(iend,jbl ),1, MPI_REAL4,iright,1 , & + & a(ibl ,jbl ),1, MPI_REAL4,ileft ,1, & + & MPI_COMM_COMP,status,ierr) + if ( ierr /= 0 ) then + print *, ' problem with tenth sendrecv in exch, ierr = ',ierr + stop 6771 + endif + + call mpi_sendrecv(a(iend,jbu ),1, MPI_REAL4,iright,1 , & + & a(ibl ,jbu ),1, MPI_REAL4,ileft ,1, & + & MPI_COMM_COMP,status,ierr) + + if ( ierr /= 0 ) then + print *, ' problem with eleventh sendrecv in exch, ierr = ',ierr + stop 6772 + endif + + call mpi_sendrecv(a(ista,jbl ),1, MPI_REAL4,ileft ,1, & + & a(ibu ,jbl ),1, MPI_REAL4,iright,1, & + & MPI_COMM_COMP,status,ierr) + + if ( ierr /= 0 ) then + print *, ' problem with twelft sendrecv in exch, ierr = ',ierr + stop 6773 + endif + + call mpi_sendrecv(a(ista,jbu ),1, MPI_REAL4,ileft ,1 , & + & a(ibu ,jbu ),1, MPI_REAL4,iright,1, & + & MPI_COMM_COMP,status,ierr) + + if ( ierr /= 0 ) then + print *, ' problem with thirteenth sendrecv in exch, ierr = ',ierr + stop 6774 + endif + + 139 format(a20,5(i10,i6,i6,'<>')) + + if(checkcoords) then + if(ifirst .le. 0) then + call mpi_sendrecv(ibcoords(iend,jbl ),1 ,MPI_INTEGER,iright,1 , & + & ibcoords(ibl ,jbl ),1 ,MPI_INTEGER,ileft ,1, & + & MPI_COMM_COMP,status,ierr) + + call mpi_sendrecv(ibcoords(iend,jbu ),1 ,MPI_INTEGER,iright,1, & + & ibcoords(ibl ,jbu ),1 ,MPI_INTEGER,ileft ,1, & + & MPI_COMM_COMP,status,ierr) + call mpi_sendrecv(ibcoords(ista,jbl ),1 ,MPI_INTEGER,ileft ,1, & + & ibcoords(ibu ,jbl ),1 ,MPI_INTEGER,iright,1, & + & MPI_COMM_COMP,status,ierr) + call mpi_sendrecv(ibcoords(ista,jbu ),1 ,MPI_INTEGER,ileft ,1 , & + & ibcoords(ibu ,jbu ),1 ,MPI_INTEGER,iright,1, & + MPI_COMM_COMP,status,ierr) + +! corner check for coordnates + + icc=ibl + jcc=jbl + ii=ibcoords(icc,jcc)/10000 + jj=ibcoords(icc,jcc)-(ii*10000) + + if(ii .ne. icc .and. icc .ne. 0) write(0,151) ' CORNER FAILI ilb ll ',icc,jcc,ibcoords(icc,jcc),ii,jj + if( jj .ne. jcc) write(0,151) ' CORNER FAILJ ilb ll ',icc,jcc,ibcoords(icc,jcc),ii,jj + + icc=ibu + jcc=jbl + ii=ibcoords(icc,jcc)/10000 + jj=ibcoords(icc,jcc)-(ii*10000) + if(ii .ne. icc .and. icc .ne. im+1 ) write(0,151) ' CORNER FAILI ilb ul ',icc,jcc,ibcoords(icc,jcc),ii,jj + if( jj .ne. jcc ) write(0,151) ' CORNER FAILJ ilb ul ',icc,jcc,ibcoords(icc,jcc),ii,jj + + icc=ibu + jcc=jbu + ii=ibcoords(icc,jcc)/10000 + jj=ibcoords(icc,jcc)-(ii*10000) + if(ii .ne. icc .and. icc .ne. im+1) write(0,151) ' CORNER FAILI ilb uu ',icc,jcc,ibcoords(icc,jcc),ii,jj + if( jj .ne. jcc ) write(0,151) ' CORNER FAILJ ilb ul ',icc,jcc,ibcoords(icc,jcc),ii,jj + + icc=ibl + jcc=jbu + ii=ibcoords(icc,jcc)/10000. + jj=ibcoords(icc,jcc)-(ii*10000) + if(ii .ne. icc .and. icc .ne. 0 ) write(0,151) ' CORNER FAILI ilb lu ',icc,jcc,ibcoords(icc,jcc),ii,jj + if( jj .ne. jcc ) write(0,151) ' CORNER FAILJ ilb ul ',icc,jcc,ibcoords(icc,jcc),ii,jj + +! if(ileft .ge. 0) then +!119 format(' GWX LEFT EXCHANGE ileft,me,ibcoords(ista-1,jend+1),ibcoords(ista-1,jend-1),ista-1,jend-1,jend+1', & +! 10i10) +! endif + +! if(iright .ge. 0) then +!! write(0,129) iright,me,ibcoords(ista+1,jend+1),ibcoords(ista+1,jend-1),ista-1,jend-1,jend+1 !GWVX +!129 format(' GWX RIGHT EXCHANGE iright,me,ibcoords(ista+1,jend+1),ibcoords(ista-1,jend+1),ista-1,jend-1,jend+1', & +! 10i10) +! endif + +! interior check + + do j=jsta,jend + do i=ista,iend + ii=ibcoords(i,j)/10000 + jj=ibcoords( i,j)-(ii*10000) + if(ii .ne. i .or. jj .ne. j) write(0,151) 'GWVX FAILED IJ ',i,j,ibcoords(i,j),ibl,jbl,ibu,jbu + end do + end do + + 151 format(a70,10i10) + +! bounds check +! first check top and bottom halo rows + + j=jbu + do i=ista,iend + ii=ibcoords(i,j)/10000 + jj=ibcoords( i,j)-(ii*10000) + if(ii .ne. i .or. jj .ne. j) write(0,151) 'GWVX FAILEDI JBU IJ ',i,j,ibcoords(i,j),ibl,jbl,ibu,jbu + end do + + j=jbl + do i=ista,iend + ii=ibcoords(i,j)/10000 + jj=ibcoords( i,j)-(ii*10000) + if(ii .ne. i .or. jj .ne. j) write(0,151) 'GWVX FAILEDI JBL IJ ',i,j,ibcoords(i,j),ibl,jbl,ibu,jbu + end do + +! second and last, check left and right halo columns + + i=ibl + do j=jsta,jend + ii=ibcoords(i,j)/10000 + jj=ibcoords( i,j)-(ii*10000) + if(ii .ne. i .and. ii .ne. im .or. jj .ne. j) write(0,151) 'GWVX FAILED IBL IJ ',ii,i,j,ibcoords(i,j),ibl,jbl,ibu,jbu + end do + + i=ibu + do j=jsta,jend + ii=ibcoords(i,j)/10000 + jj=ibcoords( i,j)-(ii*10000) + if(ii .ne. i .and. ii .ne. 1 .or. jj .ne. j) write(0,151) 'GWVX FAILED IBU ii i j ibcoords ibl,jbl,ibu,jbu',ii,i,j,ibcoords(i,j),ibl,jbl,ibu,jbu + end do + + if(me .eq. 0) write(0,*) ' IFIRST CHECK' + + endif ! IFIRST + endif !checkcoords + +! end halo checks + if ( ierr /= 0 ) then + print *, ' problem with second sendrecv in exch, ierr = ',ierr + stop + end if + call mpi_barrier(mpi_comm_comp,ierr) + ifirst=ifirst+1 + end diff --git a/sorc/ncep_post.fd/IFI.F b/sorc/ncep_post.fd/IFI.F new file mode 100644 index 000000000..c5346add7 --- /dev/null +++ b/sorc/ncep_post.fd/IFI.F @@ -0,0 +1,1424 @@ +#ifndef USE_IFI + +! IFI stubs + +module upp_ifi_mod + use iso_c_binding, only: c_float + implicit none + + private + public run_ifi, set_ifi_dims, ifi_real_t, write_ifi_debug_files, enable_ifi_smoother + + integer, parameter :: ifi_real_t = c_float + logical :: write_ifi_debug_files = .false. + logical :: enable_ifi_smoother = .false. + +contains !============================================================== + + subroutine set_ifi_dims() + use CTLBLK_mod, only: ifi_nflight, ifi_flight_levels + implicit none + integer :: i + + ! Bogus fill value for flight levels to prevent a crash + ifi_nflight = 60 + allocate(ifi_flight_levels(ifi_nflight)) + do i=1,ifi_nflight + ifi_flight_levels(i) = 500*i + enddo + end subroutine set_ifi_dims + + ! -------------------------------------------------------------------- + + subroutine run_ifi() + ! Fill any requested IFI fields with missing data. + call send_missing_data(1003) ! ICE_PROB = missing + call send_missing_data(1004) ! SLD = missing + call send_missing_data(1005) ! ICE_SEV_CAT = missing + end subroutine run_ifi + + ! -------------------------------------------------------------------- + + subroutine send_missing_data(ient) + use CTLBLK_mod, only: ifi_nflight, ifi_flight_levels + use ctlblk_mod, only: spval, ista, iend, jsta, jend, lm, im, cfld, datapd, fld_info, ifi_flight_levels, jsta_2l, jend_2u + use rqstfld_mod, only: iget, iavblfld, lvlsxml, lvls + implicit none + + integer, intent(in) :: ient + + ! Locals + integer :: i,j,k + + logical, save :: wrote_message = .false. ! guarded by an OMP CRITICAL below + + if(size(IGET)0) then + if(.not.wrote_message) then + !$OMP CRITICAL + if(.not.wrote_message) then + write(0,'(A)') 'This post cannot produce IFI icing products because it was not compiled with libIFI.' + wrote_message = .true. + endif + !$OMP END CRITICAL + endif + + cfld = cfld+1 + fld_info(cfld)%ifld = IAVBLFLD(IGET(ient)) + fld_info(cfld)%lvl = k + + !$OMP PARALLEL DO PRIVATE(i,j) COLLAPSE(2) + do j=jsta,jend + do i=ista,iend + datapd(i-ista+1,j-jsta+1,cfld) = spval + enddo + enddo + endif + enddo + end subroutine send_missing_data +end module upp_ifi_mod + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +#else + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +! Actual IFI code. +module upp_ifi_mod + use ifi_type_mod + use ifi_mod + use iso_c_binding, only: c_bool, c_int64_t, c_ptr, c_double + + implicit none + + private + public run_ifi, set_ifi_dims, ifi_real_t + + real, parameter :: feet2meters = 0.3048 + + interface + subroutine abort() bind(C) + end subroutine abort + end interface + + type(IFIConfig) :: ifi_config + logical :: have_read_ifi_config = .false. + logical, public :: write_ifi_debug_files = .false. + logical, public :: enable_ifi_smoother = .false. + + ! Communication-related variables for gathers and + ! scatters. Gathering is done in the i direction first, all to one + ! rank in each row. Then that rank in each row gathers to the + ! destination rank. Scatter follows the same process, in reverse. + + integer :: ifi_mpi_real_kind=-999 + integer :: row_comm = -999 ! communicator for i-direction messages + integer :: col_comm = -999 ! communicator for j-direction messages + integer :: row_comm_rank=-999, col_comm_rank=-999 ! Rank in the i & j comm communicators + integer :: row_comm_size=-999, col_comm_size=-999 ! Size of the i & j comm communicators + integer :: grid_rank ! rank within mpi_comm_comp + + integer, allocatable :: ista_grid(:) ! ista for each rank in mpi_comm_comp in i direction + integer, allocatable :: jsta_grid(:) ! jsta for each rank in mpi_comm_comp in j direction + + real, allocatable :: rearrange(:,:) ! for rearranging data on local processor + integer, allocatable :: row_comm_count(:) ! scatterv&gatherv receive counts for i-direction collectives + integer, allocatable :: row_comm_displ(:) ! scatterv&gatherv displacements for i-direction collectives + integer, allocatable :: row_ista(:) ! i starts for rearranging data (icomm_size) + integer, allocatable :: row_iend(:) ! i ends for rearranging data (icomm_size) + real, allocatable :: rearrange_row1d(:) ! for rearranging data for an entire row (roots of icomm communicators) + + real, allocatable :: rearrange_row2d(:,:) ! for rearranging data for an entire row (roots of icomm communicators) + integer, allocatable :: col_comm_count(:) ! scatterv&gatherv receive counts for j-direction collectives + integer, allocatable :: col_comm_displ(:) ! scatterv&gatherv displacements for j-direction collectives + integer, allocatable :: col_jsta(:) ! j starts for rearranging data (jcomm_size) + integer, allocatable :: col_jend(:) ! j ends for rearranging data (jcomm_size) + + integer, allocatable :: isize_grid(:) ! iend-ista+1 for each gridpoint + +contains !============================================================== + + subroutine make_communicators + use ctlblk_mod, only: jsta, jend, ista, iend, mpi_comm_comp, im, jm + use mpi + use iso_c_binding, only: c_sizeof + implicit none + + integer, allocatable :: key(:) ! ensure ranks are in a consistent order + integer :: size, ierr, local_count, accum, irank, i + real(kind=ifi_real_t) :: testreal + + if(allocated(jsta_grid)) then + return ! Already made communicators + endif + + ! Use the right MPI type for the ifi real kind + if(c_sizeof(testreal)==4) then + ifi_mpi_real_kind=MPI_REAL4 + else + ifi_mpi_real_kind=MPI_REAL8 + endif + + ! Get the correct size of the communicator before we try to split it. + call MPI_Comm_size(mpi_comm_comp,size,ierr) + call MPI_Comm_rank(mpi_comm_comp,grid_rank,ierr) + +! 929 format('Rank ',I0,' ista=',I0,' jsta=',I0) +! print 929,grid_rank,ista,jsta + + ! Get the start locations on every rank. We need this for the key, + ! and for determining root ranks. + allocate(ista_grid(size)) + allocate(jsta_grid(size)) + ista_grid=-1 + jsta_grid=-1 + call MPI_Allgather(ista,1,MPI_INTEGER,ista_grid,1,MPI_INTEGER,mpi_comm_comp,ierr) + call MPI_Allgather(jsta,1,MPI_INTEGER,jsta_grid,1,MPI_INTEGER,mpi_comm_comp,ierr) + + ! Ensure ranks are arranged in a consistent order. + allocate(key(size)) + do i=1,size + if(ista_grid(i)==-1) then + write(0,*) 'invalid ista_grid ',ista_grid(i) + call abort + endif + if(jsta_grid(i)==-1) then + write(0,*) 'invalid jsta_grid ',jsta_grid(i) + call abort + endif + key(i) = ista_grid(i) + im*(jsta_grid(i)-1) + enddo + + ! Create a communicator for row gather/scatter (all ranks that share a jsta) + if(all(jsta_grid==jsta_grid(1))) then + row_comm = mpi_comm_comp + else + call mpi_comm_split(mpi_comm_comp,jsta_grid(grid_rank+1),key(grid_rank+1),row_comm,ierr) + endif + if(row_comm==MPI_COMM_NULL .or. row_comm==0) then + write(0,*) 'MPI_Comm_split gave MPI_COMM_NULL for row_comm' + call abort + endif + call mpi_comm_rank(row_comm,row_comm_rank,ierr) + call mpi_comm_size(row_comm,row_comm_size,ierr) + + ! Create a communicator for column gather/scatter (all ranks that share an ista) + if(all(ista_grid==ista_grid(1))) then + col_comm = mpi_comm_comp + else + call mpi_comm_split(mpi_comm_comp,ista_grid(grid_rank+1),key(grid_rank+1),col_comm,ierr) + endif + if(col_comm==MPI_COMM_NULL .or. col_comm==0) then + write(0,*) 'MPI_Comm_split gave MPI_COMM_NULL for col_comm' + call abort + endif + call mpi_comm_rank(col_comm,col_comm_rank,ierr) + call mpi_comm_size(col_comm,col_comm_size,ierr) + + ! Done with the key. + deallocate(key) + + ! Allocate the arrays used to rearrange data. + allocate(rearrange(ista:iend,jsta:jend)) + allocate(rearrange_row1d(im*(jend-jsta+1))) + allocate(rearrange_row2d(1:im,jsta:jend)) + + ! Calculate information for row MPI_Gatherv and MPI_Scatterv calls. + allocate(row_ista(row_comm_size)) + call MPI_Allgather(ista,1,MPI_INTEGER,row_ista,1,MPI_INTEGER,row_comm,ierr) + allocate(row_iend(row_comm_size)) + call MPI_Allgather(iend,1,MPI_INTEGER,row_iend,1,MPI_INTEGER,row_comm,ierr) + allocate(row_comm_count(row_comm_size)) + allocate(row_comm_displ(row_comm_size)) + call get_count_and_displ(row_comm,row_comm_size,(iend-ista+1)*(jend-jsta+1), & + row_comm_count,row_comm_displ) + + ! Calculate information for j-direction MPI_Gatherv and MPI_Scatterv calls. + allocate(col_jsta(col_comm_size)) + call MPI_Allgather(jsta,1,MPI_INTEGER,col_jsta,1,MPI_INTEGER,col_comm,ierr) + allocate(col_jend(col_comm_size)) + call MPI_Allgather(jend,1,MPI_INTEGER,col_jend,1,MPI_INTEGER,col_comm,ierr) + allocate(col_comm_count(col_comm_size)) + allocate(col_comm_displ(col_comm_size)) + call get_count_and_displ(col_comm,col_comm_size,im*(jend-jsta+1), & + col_comm_count,col_comm_displ) + + end subroutine make_communicators + + ! -------------------------------------------------------------------- + + subroutine get_count_and_displ(comm,comm_size,here,counts,displacements) + use mpi + implicit none + integer, intent(in) :: comm,comm_size,here + integer, intent(inout) :: counts(comm_size),displacements(comm_size) + integer :: accum, ierr, i + + ! Get counts from all ranks + call MPI_Allgather(here,1,MPI_INTEGER,counts,1,MPI_INTEGER,comm,ierr) + + ! Displacements are a cumulative sum starting at 0 for rank 0 + accum=0 + do i=1,comm_size + displacements(i) = accum + accum = accum+counts(i) + end do + end subroutine get_count_and_displ + + ! -------------------------------------------------------------------- + + subroutine find_comm_roots(global_rank,row_root,col_root) + ! Find the roots of the row and column communicators for a gather + ! or scatter with the specified rank as the location of the global + ! array. + implicit none + + integer, intent(in) :: global_rank + integer, intent(inout) :: row_root,col_root + integer :: r, destination_ista, destination_jsta + + ! print *,'find roots for rank ',global_rank + + destination_ista=ista_grid(global_rank+1) + destination_jsta=jsta_grid(global_rank+1) + + row_root=-1 + do r=1,row_comm_size + if(row_ista(r)==destination_ista) then + row_root=r-1 + exit + endif + enddo + + col_root=-1 + do r=1,col_comm_size + if(col_jsta(r)==destination_jsta) then + col_root=r-1 + exit + endif + enddo + +! 201 format('root is row_root=',I0,' col_root=',I0,' for rank ',I0) +! print 201,row_root,col_root,global_rank + + if(col_root==-1) then + write(0,'(A,I0)') 'ABORT: Could not find column root rank for rank ',global_rank + call abort + endif + if(row_root==-1) then + write(0,'(A,I0)') 'ABORT: Could not find row root rank for rank ',global_rank + call abort + endif + end subroutine find_comm_roots + + ! -------------------------------------------------------------------- + + subroutine global_gather(local,global,global_rank,ista_local,iend_local,jsta_local,jend_local) + ! Gather from local arrays on all ranks to global array on specified rank + use mpi + use ctlblk_mod, only: jsta, jend, ista, iend, im, jm + use iso_c_binding, only: c_sizeof + implicit none + + real(kind=ifi_real_t), intent(in) :: local(ista_local:iend_local,jsta_local:jend_local) + real(kind=ifi_real_t), intent(out) :: global(im,jm) + integer, intent(in) :: global_rank,ista_local,iend_local,jsta_local,jend_local + + integer :: i,j,r,inindex,row_root,col_root, destination_ista, destination_jsta,ni_rank,nj_rank,idxstart, ierr + + if(col_comm==0 .or. col_comm==MPI_COMM_NULL) then + write(0,*) 'somehow, col_comm became invalid ',col_comm + call abort() + endif + + if(row_comm==0 .or. row_comm==MPI_COMM_NULL) then + write(0,*) 'somehow, row_comm became invalid ',row_comm + call abort() + endif + + ! Find out who is the root of the i and j direction communications. + call find_comm_roots(global_rank,row_root,col_root) + + ! Store the data in a contiguous local array + !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(i,j) + do j=jsta,jend + do i=ista,iend + rearrange(i,j) = local(i,j) + enddo + enddo + + ! Gather across all ranks in row + call MPI_Gatherv(rearrange,(jend-jsta+1)*(iend-ista+1),ifi_mpi_real_kind, & + rearrange_row1d,row_comm_count,row_comm_displ,ifi_mpi_real_kind,row_root,row_comm,ierr) + + if(row_comm_rank==row_root) then + ! Rearrange from i-j-rank dimensions to i-j dimensions. + idxstart=1 + do r=1,row_comm_size + + ni_rank=row_iend(r)-row_ista(r)+1 + nj_rank=jend-jsta+1 + ! print *,'r,ista,iend=',r,row_ista(r),row_iend(r) + + !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(inindex) + do j=0,nj_rank-1 + do i=0,ni_rank-1 + inindex = idxstart + i + ni_rank*j + rearrange_row2d(i+row_ista(r),j+jsta) = rearrange_row1d(inindex) + enddo + enddo + + idxstart = idxstart + ni_rank*nj_rank + enddo + + ! Gather along columns + ! print *,'global gatherv to rank ',grid_rank + call MPI_Gatherv(rearrange_row2d,im*(jend-jsta+1),ifi_mpi_real_kind, & + global,col_comm_count,col_comm_displ,ifi_mpi_real_kind,col_root,col_comm,ierr) + endif + + ! print *,'global_gather success' + end subroutine global_gather + + ! -------------------------------------------------------------------- + + subroutine global_scatter(local,global,global_rank,ista_local,iend_local,jsta_local,jend_local) + ! Scatter from global array at specified rank to local arrays on all ranks. + ! NOTE: Does NOT update halo regions! + use mpi + use ctlblk_mod, only: jsta, jend, ista, iend, im, jm + use iso_c_binding, only: c_sizeof + implicit none + + real(kind=ifi_real_t), intent(out) :: local(ista_local:iend_local,jsta_local:jend_local) + real(kind=ifi_real_t), intent(in) :: global(im,jm) + integer, intent(in) :: global_rank,ista_local,iend_local,jsta_local,jend_local + + integer :: i,j,r,outindex,idxstart,ni_rank,nj_rank,col_root,ierr,row_root + + if(col_comm==0 .or. col_comm==MPI_COMM_NULL) then + write(0,*) 'somehow, col_comm became invalid ',col_comm + call abort() + endif + + if(row_comm==0 .or. row_comm==MPI_COMM_NULL) then + write(0,*) 'somehow, row_comm became invalid ',row_comm + call abort() + endif + + ! Find out who is the root of the i and j direction communications. + call find_comm_roots(global_rank,row_root,col_root) + + if(row_comm_rank==row_root) then + ! Distribute in J direction. + ! print *,'global scatterv from rank ',grid_rank + call MPI_Scatterv(global, col_comm_count, col_comm_displ, ifi_mpi_real_kind, & + rearrange_row2d ,im*(jend-jsta+1), ifi_mpi_real_kind, col_root, col_comm, ierr) + ! Rearrange from i-j dimensions to i-j-rank dimensions + idxstart=1 + do r=1,row_comm_size + + ni_rank=row_iend(r)-row_ista(r)+1 + nj_rank=jend-jsta+1 + + !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(outindex) + do j=0,nj_rank-1 + do i=0,ni_rank-1 + outindex = idxstart + i + ni_rank*j + rearrange_row1d(outindex) = rearrange_row2d(i+row_ista(r),j+jsta) + enddo + enddo + + idxstart = idxstart + ni_rank*nj_rank + enddo + endif + + ! Distribute across row + call MPI_Scatterv(rearrange_row1d,row_comm_count,row_comm_displ,ifi_mpi_real_kind, & + rearrange,(iend-ista+1)*(jend-jsta+1),ifi_mpi_real_kind,row_root,row_comm,ierr) + + ! Copy back to contiguous local array. + ! NOTE: Does NOT update the halo! + !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(i,j) + do j=jsta,jend + do i=ista,iend + local(i,j) = rearrange(i,j) + enddo + enddo + + ! print *,'global_scatter success' + + end subroutine global_scatter + + ! -------------------------------------------------------------------- + + subroutine ifi_check(status,error_message) + use mpi, only: MPI_Abort, MPI_COMM_WORLD + implicit none + integer(c_int64_t), intent(in) :: status + character(*), intent(in) :: error_message + integer :: ierr + + ! Exit program if status is non-zero + if(status/=0) then + write(0,'("IFI Failed: ",A)') trim(error_message) + call MPI_Abort(MPI_COMM_WORLD,1,ierr) + endif + end subroutine ifi_check + + ! -------------------------------------------------------------------- + + subroutine read_ifi_config() + implicit none + ! Only read the config if we haven't already: + if(.not.have_read_ifi_config) then + call ifi_check(ifi_config%init(),'read IFI config') + have_read_ifi_config = .true. + endif + end subroutine read_ifi_config + + ! -------------------------------------------------------------------- + + subroutine ifi_minmax_callback(minv,maxv) bind(C) + use ctlblk_mod, only: mpi_comm_comp, me + use mpi + implicit none + real(kind=ifi_real_t) :: minv,maxv + real(kind=8) :: sendbuf(1),recvbuf(1) + integer :: ierr + + sendbuf(1)=minv + call MPI_Allreduce(sendbuf,recvbuf,1,MPI_REAL8,MPI_MIN,mpi_comm_comp,ierr) + minv=recvbuf(1) + + sendbuf(1)=maxv + call MPI_Allreduce(sendbuf,recvbuf,1,MPI_REAL8,MPI_MAX,mpi_comm_comp,ierr) + maxv=recvbuf(1) + end subroutine ifi_minmax_callback + + ! -------------------------------------------------------------------- + + subroutine ifi_halo_callback(d,ims1,ime1,jms1,jme1, & + ids1,ide1,jds1,jde1, & + ips1,ipe1,jps1,jpe1) BIND(C) + use ctlblk_mod, only: jsta,jend,im,jsta_2l,jend_2u,ista,iend,ista_2l,iend_2u + implicit none + real(kind=ifi_real_t) :: d(*) + integer(kind=c_int64_t), value :: & + ims1,ime1,jms1,jme1, & + ids1,ide1,jds1,jde1, & + ips1,ipe1,jps1,jpe1 + integer(kind=c_int64_t) :: & + ims,ime,jms,jme, & + ids,ide,jds,jde, & + ips,ipe,jps,jpe + + ! Convert from C++ to Fortran indexing: + ims=ims1+1 ; ime=ime1+1 ; jms=jms1+1 ; jme=jme1+1 + ids=ids1+1 ; ide=ide1+1 ; jds=jds1+1 ; jde=jde1+1 + ips=ips1+1 ; ipe=ipe1+1 ; jps=jps1+1 ; jpe=jpe1+1 + + if( jms/=jsta_2l .or. jme/=jend_2u .or. jsta/=jps .or. jend/=jpe & + .or. ims/=ista_2l .or. ime/=iend_2u .or. ista/=ips .or. iend/=ipe ) then +94 format(' ',A,' = ',I0,', ',I0) +95 format(' ',A,' = ',I0,', ',I0,', ',I0) + print '(A)','Warning: IFI halo callback j bounds do not match UPP j bounds' + print 94,'jps,jsta',jps,jsta + print 94,'jpe,jend',jpe,jend + print 94,'jms,jsta_2l',jms,jsta_2l + print 94,'jme,jend_2u',jme,jend_2u + print 94,'ips,ista',ips,ista + print 94,'ipe,iend',ipe,iend + print 94,'ims,ista_2l',ims,ista_2l + print 94,'ime,iend_2u',ime,iend_2u + !call ifi_check(-1,'Internal error: IFI halo callback i/j bounds do not match UPP i/j bounds') + end if + call EXCH_c_float(d) + end subroutine ifi_halo_callback + + ! -------------------------------------------------------------------- + + subroutine set_ifi_dims() + use CTLBLK_mod, only: ifi_nflight, ifi_flight_levels + implicit none + + ! Warning: do not deallocate config_flight_levels_feet. + ! The IFI C++ library manages that array. + integer(kind=c_int64_t), pointer :: config_flight_levels_feet(:) + integer :: i + + ! Make sure the config was read in + call read_ifi_config() + + ! Get the flight levels + call ifi_check(ifi_config%get_flight_levels_feet(config_flight_levels_feet), & + 'cannot get flight levels in feet from IFI config') + + ! Convert from integer to real: + ifi_nflight = size(config_flight_levels_feet) + allocate(ifi_flight_levels(ifi_nflight)) + ifi_flight_levels = config_flight_levels_feet + +! print '(A)','IFI flight levels:' +! 38 format(' ifi_flight_level[',I0,'] = ',F15.5) +! do i=1,ifi_nflight +! print 38,i,ifi_flight_levels(i) +! enddo + + end subroutine set_ifi_dims + + ! -------------------------------------------------------------------- + + subroutine send_data(vars,name,ient) + use ctlblk_mod, only: spval, jsta, jend, lm, cfld, datapd, fld_info, ifi_flight_levels, jsta_2l, jend_2u, ista, iend, ista_2l, iend_2u + use rqstfld_mod, only: iget, iavblfld, lvlsxml, lvls + implicit none + + class(IFIData) :: vars + integer, intent(in) :: ient + character(*), intent(in) :: name + + ! Locals + integer(c_int64_t) :: & + ims,ime,jms,jme,kms,kme, ids,ide,jds,jde,kds,kde, ips,ipe,jps,jpe,kps,kpe + real(kind=ifi_real_t), pointer :: data(:) + logical(c_bool) :: missing_value_is_set + real(kind=ifi_real_t) :: missing_value + integer :: i,j,k,nj_local,jpad,ilen,jlen,kstartm1,jstartm1,iloc,ndata,ji,ni_local,ipad + + if(.not. IGET(ient)>0) then + return + endif + + ! WARNING: do not deallocate the data pointer. It is managed by the IFI C++ library. + data => vars%get_data(trim(name),missing_value_is_set,missing_value, & + ims,ime,jms,jme,kms,kme, ids,ide,jds,jde,kds,kde, ips,ipe,jps,jpe,kps,kpe) + if(.not.missing_value_is_set) then + missing_value = MISSING + endif + + ! Get dimensions and do some sanity checks + nj_local = jpe-jps+1 + ni_local = ipe-ips+1 + jpad = jps-jms + ipad = ips-ims + if(nj_local/=jend-jsta+1 .or. jpad/=jsta-jsta_2l .or. ni_local/=iend-ista+1 .or. ipad/=ista-ista_2l) then +94 format(' ',A,' = ',I0,', ',I0) +83 format('Warning: ',A,': IFI output j bounds do not match UPP j bounds') + print 83,trim(name) + print 94,'jps,jsta',jps,jsta + print 94,'jpe,jend',jpe,jend + print 94,'jms,jsta_2l',jms,jsta_2l + print 94,'jme,jend_2u',jme,jend_2u + print 94,'jpad,jsta-jsta_2l',jpad,jsta-jsta_2l + print 94,'nj_local,jend-jsta+1',nj_local,jend-jsta+1 + print 94,'ips,ista',ips,ista + print 94,'ipe,iend',ipe,iend + print 94,'ims,ista_2l',ims,ista_2l + print 94,'ime,iend_2u',ime,iend_2u + print 94,'ipad,ista-ista_2l',ipad,ista-ista_2l + print 94,'ni_local,iend-ista+1',ni_local,iend-ista+1 + !call ifi_check(-1,'Internal error: IFI output j bounds do not match UPP j bounds') + end if + ilen=ime-ims+1 + jlen=jme-jms+1 + ndata=ilen*jlen*(kme-kms+1) + + ! Go level-by-level writing grib2 if requested + do k=kps,kpe + kstartm1=ilen*jlen*(k-kms) + if(LVLS(k,IGET(ient))>0) then + cfld = cfld+1 + fld_info(cfld)%ifld = IAVBLFLD(IGET(ient)) + fld_info(cfld)%lvl = k ! ifi_flight_levels(k)*feet2meters + !$OMP PARALLEL DO PRIVATE(i,j,iloc,jstartm1) COLLAPSE(2) + do j=jps,jpe + do i=ips,ipe + jstartm1 = kstartm1 + (j-jms)*ilen + iloc = jstartm1+(i-ims)+1 + if(iloc<1 .or. iloc>size(data)) then + call ifi_check(-1,'Internal error: out of bounds in send_data.') + endif + if(data(iloc)==missing_value) then + datapd(i-ips+1,j-jps+1,cfld) = spval + else + datapd(i-ips+1,j-jps+1,cfld) = data(iloc) + endif + enddo + enddo + endif + enddo + end subroutine send_data + + ! -------------------------------------------------------------------- + + subroutine gather_ifi(local_buf_c,global_buf_c,receiving_rank) bind(C) + use iso_c_binding, only: c_float, c_double,c_int64_t,c_f_pointer + use mpi + use ctlblk_mod, only: ME, MPI_COMM_COMP, JSTA, JEND, IM, JM, ICNT, IDSP, ISTA, IEND + implicit none + + integer(kind=c_int64_t), value :: receiving_rank + type(c_ptr), value :: global_buf_c, local_buf_c + real(kind=ifi_real_t), pointer :: global_buf(:,:), local_buf(:,:) + + integer :: type, iret + + call c_f_pointer(local_buf_c,local_buf,(/ iend-ista+1, jend-jsta+1 /)) + if(me==receiving_rank) then + call c_f_pointer(global_buf_c,global_buf,(/ im, jm /)) + else + call c_f_pointer(global_buf_c,global_buf,(/ 1,1 /)) + endif + + call gather_for_write(local_buf,global_buf,receiving_rank) + + end subroutine gather_ifi + + ! -------------------------------------------------------------------- + + subroutine gather_for_write(local_buf,global_buf,receiving_rank_c) bind(C) + use iso_c_binding, only: c_float, c_double,c_int64_t,c_f_pointer + use mpi + use ctlblk_mod, only: ME, MPI_COMM_COMP, JSTA, JEND, IM, JM, ISTA, ISTA, IEND + implicit none + + integer(kind=c_int64_t), value :: receiving_rank_c + real(kind=ifi_real_t) :: global_buf(:,:), local_buf(:,:) + integer :: receiving_rank + receiving_rank = receiving_rank_c + call global_gather(local_buf,global_buf,receiving_rank,ista,iend,jsta,jend) + end subroutine gather_for_write + + ! -------------------------------------------------------------------- + + subroutine scatter_ifi(local_buf_c,global_buf_c,sending_rank_c)bind(C) + use iso_c_binding, only: c_float, c_double, c_int64_t, c_f_pointer + use mpi + use ctlblk_mod, only: JSTA, JEND, IM, JM, ISTA, IEND, ME + implicit none + + integer(kind=c_int64_t), value :: sending_rank_c + type(c_ptr), value :: global_buf_c, local_buf_c + real(kind=ifi_real_t), pointer :: global_buf(:,:), local_buf(:,:) + integer :: sending_rank + + sending_rank=sending_rank_c + + call c_f_pointer(local_buf_c,local_buf,(/ iend-ista+1, jend-jsta+1 /)) + if(me==sending_rank) then + call c_f_pointer(global_buf_c,global_buf,(/ im, jm /)) + else + call c_f_pointer(global_buf_c,global_buf,(/ 1,1 /)) + endif + + call global_scatter(local_buf,global_buf,sending_rank,ista,iend,jsta,jend) + end subroutine scatter_ifi + + ! -------------------------------------------------------------------- + + subroutine run_ifi() + use ctlblk_mod, only: spval, lm, lp1, im, jsta_2l,jend_2u, ITPREC, IFHR, IFMIN, grib, & + jm, jsta,jend, me, num_procs, mpi_comm_comp, ista, iend, ista_2l, iend_2u, dtq2 + use vrbls3d, only: & + zmid, & ! = IFI "HGT" + zint, & ! zint(:,:,LM+1) = IFI "HGT_surface" + QQI, & ! = IFI "CIMIXR" + QQG, & ! = IFI "GRLE" + QQW, & ! = IFI "CLMR" + pmid, & ! = IFI "PRES" + QQR, & ! "RWMR" + q, & ! "SPFH" + t, & ! "TMP" + QQS, & ! "SNMR" + OMGA ! "VVEL" + use vrbls2d, only : & + IFI_APCP, & ! IFI "APCP_surface" over ITPREC bucket time + CAPE, CIN, & + AVGPREC_CONT ! Backup plan for points where IFI_APCP is invalid or 0 + use rqstfld_mod, only: IGET + use iso_c_binding, only: c_bool, c_int64_t + + implicit none + + INTERFACE ! implemented in EXCH_c_float.f + SUBROUTINE EXCH_c_float(A) BIND(C) + use ifi_type_mod, only: ifi_real_t + implicit none + real(kind=ifi_real_t) :: a(*) + END SUBROUTINE EXCH_c_float + END INTERFACE + + type(IFIData) :: hybr_vars, pres_vars, derived_vars, fip_algo_vars, flight_vars, cat_vars + type(IFIAlgo) :: algo + character(88) :: outfile + real :: to_hourly + real(c_double) :: fcst_lead_sec + integer :: i, j + + ! FIXME: Once the post is i-decomposed, replace the appropriate 1..im in this routine. + + if(grib/='grib2' .or. (IGET(1003)<=0 .and. IGET(1004)<=0 .and. IGET(1005)<=0)) then + if(me==0) then + print '(A)','IFI fields were not requested; skipping IFI' + endif + return ! nothing to do + endif + + if(me==0) then + print '(A)','Running libIFI to get icing products.' + call ifi_print_copyright() + + if(write_ifi_debug_files) then + print '(A)','Will output IFI debug files.' + else + print '(A)','Will NOT output IFI debug files.' + endif + endif + + call make_communicators() + + ! Read config and initialize input data structures: + + call read_ifi_config() + call ifi_check(hybr_vars%init(int(1,c_int64_t),int(im,c_int64_t),int(1,c_int64_t),int(jm,c_int64_t),& + int(1,c_int64_t),int(lm,c_int64_t),int( ista,c_int64_t),int(iend,c_int64_t),int(jsta,c_int64_t),& + int(jend,c_int64_t),int(1,c_int64_t),int(lm,c_int64_t)),& + 'could not initialize IFI input data structures') + + ! Copy 2D vars to IFI internal structures: + + call ifi_check(hybr_vars%add_ij_var('topography',int(ista_2l,c_int64_t),int(iend_2u,c_int64_t),& + int(jsta_2l,c_int64_t),int(jend_2u,c_int64_t), & + zint(:,:,LP1)),'could not send topography to IFI') + + to_hourly=3600*1000.0/dtq2 +!$OMP PARALLEL DO COLLAPSE(2) + do j=jsta,jend + do i=ista,iend + if(ifi_apcp(i,j)>9e9 .or. ifi_apcp(i,j)<-9e9 .or. ifi_apcp(i,j)==0) then + ifi_apcp(i,j) = avgprec_cont(i,j)*to_hourly + else if(ITPREC>1e-5) then + ifi_apcp(i,j) = ifi_apcp(i,j)/itprec + endif + enddo + enddo + + call ifi_check(hybr_vars%add_ij_var('APCP_surface',int(ista_2l,c_int64_t),int(iend_2u,c_int64_t),& + int(jsta_2l,c_int64_t),int(jend_2u,c_int64_t), & + IFI_APCP),'could not send IFI_APCP to IFI') + call ifi_check(hybr_vars%add_ij_var('CAPE_surface',int(ista_2l,c_int64_t),int(iend_2u,c_int64_t), & + int(jsta_2l,c_int64_t),int(jend_2u,c_int64_t),CAPE), & + 'could not send CAPE to IFI') + call ifi_check(hybr_vars%add_ij_var('CIN_surface',int(ista_2l,c_int64_t),int(iend_2u,c_int64_t), & + int(jsta_2l,c_int64_t),int(jend_2u,c_int64_t),CIN), & + 'could not send CIN to IFI') + + ! Copy 3D vars to IFI internal structures, inverting K dimension + + call add_ijk_var('HGT','zmid',zmid) + call add_ijk_var('CIMIXR','QQI',QQI) + call add_ijk_var('CLMR','QQW',QQW) + call add_ijk_var('GRLE','QQG',QQG) + call add_ijk_var('PRES','pmid',pmid) + call add_ijk_var('RWMR','QQR',QQR) + call add_ijk_var('SPFH','Q',Q) + call add_ijk_var('TMP','T',T) + call add_ijk_var('SNMR','QQS',QQS) + call add_ijk_var('VVEL','OMGA',OMGA) + +308 format(A,'_',I0,'.nc') + + fcst_lead_sec=IFHR*3600.0+IFMIN*60.0 + + if(write_ifi_debug_files) then + write(outfile,308) 'hybr_vars',me + call write_fip_output(hybr_vars,fcst_lead_sec,trim(outfile),.true.,'z0','z1') + endif + + ! Initialize the IFI algorithm + + call ifi_check(algo%init(ifi_config,fcst_lead_sec,hybr_vars,ME,NUM_PROCS,MPI_COMM_COMP), & + 'could not initialize IFI algorithm') + call ifi_check(algo%set_halo_callback(ifi_halo_callback), & + 'send halo callback function to IFI algorithm') + call ifi_check(algo%set_minmax_callback(ifi_minmax_callback), & + 'send minmax callback function to IFI algorithm') + call ifi_check(algo%set_gather_callback(gather_ifi), & + 'send gather callback function to IFI algorithm') + call ifi_check(algo%set_scatter_callback(scatter_ifi), & + 'send scatter callback function to IFI algorithm') + + ! Run the IFI algorithm + + call ifi_check(algo%calc_exner(),'calc_exner() failed') + call ifi_check(algo%hybrid_to_pressure(),'hybrid_to_pressure() failed') + + call ifi_check(algo%get_pres_vars(pres_vars),'get_pres_vars()') + if(write_ifi_debug_files) then + write(outfile,308) 'pres_vars',me + call write_fip_output(pres_vars,fcst_lead_sec,trim(outfile),.true.,'z1','z0') + endif + + call ifi_check(algo%discard_hybrid_level_vars(),'discard_hybrid_level_vars() failed') + + call ifi_check(algo%derive_fields(),'derive_fields() failed') + + call ifi_check(algo%get_derived_vars(derived_vars),'get_derived_vars()') + if(write_ifi_debug_files) then + write(outfile,308) 'derived_vars',me + call write_fip_output(derived_vars,fcst_lead_sec,trim(outfile),.false.,'z1','z0') + endif + + call ifi_check(algo%run_fip_algo(),'run_fip_algo() failed') + + call ifi_check(algo%get_fip_algo_vars(pres_vars),'get_fip_algo_vars()') + if(write_ifi_debug_files) then + write(outfile,308) 'fip_algo_vars',me + call write_fip_output(pres_vars,fcst_lead_sec,trim(outfile),.false.,'z1','z0') + endif + + call ifi_check(algo%discard_pressure_level_vars(),'discard_pressure_level_vars() failed') + call ifi_check(algo%discard_derived_vars(),'discard_derived_vars() failed') + call ifi_check(algo%pressure_to_flight(),'pressure_to_flight() failed') + + call ifi_check(algo%get_flight_vars(flight_vars),'get_flight_vars()') + if(write_ifi_debug_files) then + write(outfile,308) 'flight_vars',me + call write_fip_output(flight_vars,fcst_lead_sec,trim(outfile),.false.,'z1','z0') + endif + + call ifi_check(algo%discard_fip_algo_vars(),'discard_fip_algo_vars() failed') + call ifi_check(algo%make_icing_category(),'make_icing_category() failed') + + ! Get the final output fields: + + call ifi_check(algo%get_cat_vars(cat_vars),'get_cat_vars()') + + if(write_ifi_debug_files) then + write(outfile,308) 'cat_vars',me + call write_fip_output(cat_vars,fcst_lead_sec,trim(outfile),.false.,'z1','z0') + endif + + call send_data(cat_vars,'ICE_PROB',1003) + call send_data(cat_vars,'SLD',1004) + call send_data(cat_vars,'ICE_SEV_CAT',1005) + + ! When this subroutine ends, a Fortran-2003-compliant compiler + ! will free all memory IFI uses, by calling the destructors (final + ! routines) for cat_vars, hybr_vars, algo, and config. + + contains + + subroutine add_ijk_var(ifi_name,upp_name,upp_var) + implicit none + character(len=*), intent(in) :: ifi_name,upp_name + real, intent(in) :: upp_var(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U,LM) + real(kind=ifi_real_t) :: ifi_var(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U,LM) + integer i,j,k + +!$OMP PARALLEL DO COLLAPSE(2) + do k=1,lm + do j=jsta_2l,jend_2u + do i=ista_2l,iend_2u + ifi_var(i,j,k) = upp_var(i,j,lm-k+1) + enddo + enddo + enddo + + call ifi_check(hybr_vars%add_ijk_var(ifi_name,int(ista_2l,c_int64_t),int(iend_2u,c_int64_t),& + int(jsta_2l,c_int64_t),int(jend_2u,c_int64_t),& + int(1,c_int64_t),int(lm,c_int64_t),ifi_var), & + 'could not send '//ifi_name//' ('//upp_name//') to IFI') + end subroutine add_ijk_var + + end subroutine run_ifi + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +! FIXME: DELETE THIS TEST CODE + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + subroutine find_range(var,count,missing_value_is_set,missing_value,min_not_miss,max_not_miss,all_missing) + USE ieee_arithmetic + use mpi + use ctlblk_mod, only: mpi_comm_comp + implicit none + real(kind=ifi_real_t), intent(in) :: var(*) + real(kind=ifi_real_t), intent(in) :: missing_value + real(kind=ifi_real_t), intent(out) :: min_not_miss,max_not_miss + integer, intent(in) :: count + logical, intent(out) :: all_missing + logical(c_bool), intent(in) :: missing_value_is_set + + real(kind=ifi_real_t) :: minv,maxv,epsilon, global_minv,global_maxv + integer :: i,type,iret + real(kind=ifi_real_t), parameter :: zero = 0 + + print *,'find range begin' + + if(missing_value_is_set) then + epsilon = abs(missing_value)*1e-4 + if(.not. epsilon+1>epsilon) then + epsilon=1 + endif + endif + + ! Initialize to out-of-bounds values so they'll stay that way if no valid values are found: + minv = 1e30 + maxv = -1e30 + + ! Find the min and max values in the array: + if(missing_value_is_set) then + !$OMP PARALLEL DO REDUCTION(min:minv) REDUCTION(max:maxv) + do i=1,count + if(.not. abs(var(i)-missing_value)>epsilon .and. var(i)<9e9 .and. var(i)>-9e9) then + minv=min(minv,var(i)) + maxv=max(maxv,var(i)) + endif + end do + else + !$OMP PARALLEL DO REDUCTION(min:minv) REDUCTION(max:maxv) + do i=1,count + if(var(i)+1>var(i)) then + if(var(i)<9e9 .and. var(i)>-9e9) then + minv=min(minv,var(i)) + maxv=max(maxv,var(i)) + endif + endif + end do + endif + + if(ifi_real_t==c_double) then + type=MPI_REAL8 + else + type=MPI_REAL4 + endif + + ! call MPI_Allreduce(minv,global_minv,1,type,MPI_MIN,mpi_comm_comp,iret) + ! call MPI_Allreduce(maxv,global_maxv,1,type,MPI_MAX,mpi_comm_comp,iret) + global_minv=minv + global_maxv=maxv + + ! If min or max are inf, -inf, or NaN, assume all values are missing. + all_missing = global_minv<-9e9 .or. global_minv>9e9 .or. global_maxv<-9e9 .or. global_maxv>9e9 + + print *,'range min,max = ',global_minv,global_maxv + + if(all_missing) then + min_not_miss=-1 + max_not_miss=1 + else + min_not_miss=global_minv + max_not_miss=global_maxv + endif + + print *,'find range end' + end subroutine find_range + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + subroutine nc_check(code,file,message) + use netcdf + implicit none + integer, intent(in) :: code + character(*), intent(in) :: file, message + integer :: i + if(code/=0) then + write(0,93) trim(file),trim(message),trim(nf90_strerror(code)),code + call abort + endif + 93 format(A,': ',A,': ',A,'(',I0,')') + end subroutine nc_check + + subroutine write_fip_output(ifi_data,fcst_lead_sec,output_file,rename,z2dname,z3dname) + use ifi_mod + use iso_c_binding + use netcdf + use ctlblk_mod, only: me + implicit none + integer, parameter :: maxname=80 + integer, parameter :: maxvars=999 + logical, intent(in) :: rename + double precision, intent(in) :: fcst_lead_sec + character(len=*), intent(in) :: z2dname,z3dname + type var_info + character(len=maxname) :: varname, outname + integer :: varid, ndims, dims(4), dimids(4), count + real(ifi_real_t), pointer :: data(:) + integer :: ims,ime,jms,jme,kms,kme + integer :: ids,ide,jds,jde,kds,kde + integer :: ips,ipe,jps,jpe,kps,kpe + logical :: should_dealloc + real(kind=ifi_real_t) :: missing_value + real(kind=ifi_real_t) :: min_not_miss,max_not_miss + logical(c_bool) :: missing_value_is_set + logical :: all_missing + end type var_info + character(len=*), intent(in) :: output_file + + type(IFIData) :: ifi_data + integer(c_int64_t) :: ids,ide, jds,jde, kds,kde + integer(c_int64_t) :: ips,ipe, jps,jpe, kps,kpe + integer(c_int64_t) :: ims,ime, jms,jme, kms,kme + integer :: count + + character(len=200) :: varname,nextname + + integer :: dimids(4), ncid, ivar, id, nvars, nx0, ny0, nz0, ntime, dims(4) + type(var_info),target :: var_data(maxvars) + + if(me==0) then + write(0,'(A,A)') trim(output_file),': writing IFI debug data' + endif + + call ifi_check(ifi_data%get_dims(ids,ide, jds,jde, kds,kde, ips,ipe, jps,jpe, kps,kpe),"get_dims") + + nx0=ide-ids+1 + ny0=jde-jds+1 + nz0=kde-kds+1 + ntime=1 + + dims = (/nx0,ny0,nz0,ntime/) + + varname=' ' ! special value indicating "start at the beginning" + + ivar=0 + do while(ivar read_var(var_data(ivar),trim(varname),var_data(ivar)%should_dealloc) + if(me==0) then + call set_dims(var_data(ivar),dims,dimids) + var_data(ivar)%varid = def_var(var_data(ivar),rename) + endif + enddo + + if(me/=0) then + ! Ranks that are not writing are done now + return + endif + + nvars=ivar + + if(nvars<1) then + write(0,"(A,A)") output_file,': no data to write!' + return + endif + + call nc_check(nf90_enddef(ncid),output_file,'nf90_enddef') + + write(0,*) 'before write loop ',me + do ivar=1,nvars +24 format(' put var ',A) + print 24,trim(var_data(ivar)%varname) + + call find_range(var_data(ivar)%data,var_data(ivar)%count, & + var_data(ivar)%missing_value_is_set,var_data(ivar)%missing_value,& + var_data(ivar)%min_not_miss,var_data(ivar)%max_not_miss, & + var_data(ivar)%all_missing) + call write_var(var_data(ivar)) + if(var_data(ivar)%should_dealloc) then + deallocate(var_data(ivar)%data) + endif + nullify(var_data(ivar)%data) + enddo + + call nc_check(nf90_close(ncid),output_file,"nf90_close") + contains + + subroutine set_dims(var,dims,dimids) + implicit none + type(var_info), intent(inout), target :: var + integer, intent(in) :: dims(:), dimids(:) + character(len=:), pointer :: varname + + varname=>var%varname(1:len_trim(var%varname)) + + if(me/=0) then + ! Ranks that are not writing are done now + return + endif + + ! These must also be in IFITest.cc IFITest::write_netcdf + + if(varname=='x' .or. varname=='x0') then + var%dims=(/nx0,1,1,1/) + var%dimids=(/dimids(1),-1,-1,-1/) + var%ndims=1 + else if(varname=='y' .or. varname=='y0') then + var%dims=(/ny0,1,1,1/) + var%dimids=(/dimids(2),-1,-1,-1/) + var%ndims=1 + else if(varname=='z' .or. varname=='z0' .or. varname=='z1' .or. & + varname=='pressure_levels' .or. varname=='exner_levels') then + var%dims=(/nz0,1,1,1/) + var%dimids=(/dimids(3),-1,-1,-1/) + var%ndims=1 + else if(varname=='latitude' .or. varname=='longitude') then + var%dims=(/nx0,ny0,1,1/) + var%dimids=(/dimids(1),dimids(2),-1,-1/) + var%ndims=2 + else if(varname=='time') then + var%dims=(/ntime,1,1,1/) + var%dimids=(/dimids(4),-1,-1,-1/) + var%ndims=1 + else if(kme<=kms) then + var%dims=(/nx0,ny0,ntime,1/) + var%dimids=(/dimids(1),dimids(2),dimids(4),-1/) + var%ndims=3 + else + var%dims=(/nx0,ny0,nz0,ntime/) + var%dimids=dimids + var%ndims=4 + endif + + var%count=var%dims(1)*var%dims(2)*var%dims(3)*var%dims(4) + end subroutine set_dims + + subroutine write_var(var) + implicit none + type(var_info), intent(inout) :: var + + integer :: ones(var%ndims), dims(var%ndims) + real(ifi_real_t) :: put(var%count) + integer :: i,j,k,n,ilen,jlen,m,imlen,jmlen + + ones = 1 + dims = var%dims(1:var%ndims) + + if(me/=0) then + ! Ranks that are not writing are done now + return + endif + + call nc_check(nf90_put_var(ncid=ncid,varid=var%varid,values=var%data, & + start=ones,count=dims), & + output_file,"nf90_put_var "//trim(var%outname)) + + call nc_check(nf90_put_att(ncid,var%varid,"min_value",var%min_not_miss), & + output_file,"nf90_put_att "//trim(var%outname)//" min_value") + + call nc_check(nf90_put_att(ncid,var%varid,"max_value",var%max_not_miss), & + output_file,"nf90_put_att "//trim(var%outname)//" max_value") + + end subroutine write_var + + integer function def_var(var,rename) + use iso_c_binding, only: c_float + implicit none + logical :: rename + type(var_info), intent(inout) :: var + + integer :: varid, xtype, dimids(var%ndims) + character(len=100) :: outname + + if(rename) then + var%outname=var%varname + select case(trim(var%varname)) + case('CIMIXR') + var%outname = 'ICMR' + case('SPFH') + var%outname = 'MIXR' + case('GRLE') + var%outname = 'GRMR' + case('CAPE_surface') + var%outname = 'CAPE' + case('CIN_surface') + var%outname = 'CIN' + case('CLMR') + var%outname = 'CLWMR' + case('APCP_surface') + var%outname = 'APCP1Hr' + end select + endif + + var%ims=ims ; var%ime=ime ; var%jms=jms ; var%jme=jme ; var%kms=kms ; var%kme=kme + var%ids=ids ; var%ide=ide ; var%jds=jds ; var%jde=jde ; var%kds=kds ; var%kde=kde + var%ips=ips ; var%ipe=ipe ; var%jps=jps ; var%jpe=jpe ; var%kps=kps ; var%kpe=kpe + + dimids = var%dimids(1:var%ndims) + + if(ifi_real_t==c_float) then + xtype = NF90_FLOAT + else + xtype = NF90_DOUBLE + endif + + call nc_check(nf90_def_var(ncid,trim(var%outname),xtype,dimids,def_var), & + output_file,"nf90_def_var "//trim(var%outname)) + + call nc_check(nf90_put_att(ncid,def_var,"min_value",var%min_not_miss), & + output_file,"nf90_put_att "//trim(var%outname)//" min_value") + + call nc_check(nf90_put_att(ncid,def_var,"max_value",var%max_not_miss), & + output_file,"nf90_put_att "//trim(var%outname)//" max_value") + + if(var%missing_value_is_set) then + call nc_check(nf90_put_att(ncid,def_var,"_FillValue",-9999.0), & + output_file,"nf90_put_att "//trim(var%outname)//" _FillValue") + end if + + end function def_var + + function read_var(var,varname,should_dealloc) + implicit none + type(var_info), intent(inout) :: var + real(kind=ifi_real_t), pointer :: read_var(:) + real(kind=ifi_real_t), pointer :: local_data_1D(:) + character(len=*), intent(in) :: varname + logical, intent(out) :: should_dealloc + + real(kind=ifi_real_t), allocatable :: local_data(:,:),global_data(:,:) + real(kind=ifi_real_t), pointer :: global_data_1D(:) + + integer :: nxny_local,nxny_global,nz,count,i,j,k + integer :: local_ilen,local_jlen,local_klen,local_index + integer :: global_ilen,global_jlen,global_klen,global_index + type(c_ptr) :: global_cptr,local_cptr + + local_data_1D => ifi_data%get_data(trim(varname),var%missing_value_is_set,var%missing_value, & + ims,ime,jms,jme,kms,kme, ids,ide,jds,jde,kds,kde, ips,ipe,jps,jpe,kps,kpe) + + global_ilen=ide-ids+1 + global_jlen=jde-jds+1 + global_klen=kde-kds+1 + local_ilen=ipe-ips+1 + local_jlen=jpe-jps+1 + local_klen=kpe-kps+1 + + if(.not.associated(local_data_1D)) then + write(0,38) trim(varname) +38 format("IFI did not produce ",A," variable.") + call abort + end if + + count=global_ilen*global_jlen*global_klen + if(count<=0) then + write(0,39) trim(varname),count +39 format("IFI variable ",A," had no data (size=",I0,")") + call abort + endif + + if(ids==ide .or. jds==jde) then + ! This is a 1D variable so we write it as is from rank 0 + read_var => local_data_1D + should_dealloc = .false. + return + endif + +40 format('var ',A,' has bad ',A,'. Got: ',I0,' but expected: ',I0) + + if(global_ilen /= nx0) then + write(0,40) trim(varname),'global_ilen',global_ilen,nx0 + call abort + endif + + if(global_jlen /= ny0) then + write(0,40) trim(varname),'global_jlen',global_jlen,ny0 + call abort + endif + + if(global_klen/=1 .and. global_klen /= nz0) then + write(0,40) trim(varname),'global_klen',global_klen,nz0 + call abort + endif + + + allocate(local_data(local_ilen,local_jlen)) + if(me==0) then + allocate(global_data(global_ilen,global_jlen)) + allocate(global_data_1D(global_ilen*global_jlen*global_klen)) + else + allocate(global_data(1,1)) + allocate(global_data_1D(1)) + endif + + do k=kds,kde + do j=jps,jpe + do i=ips,ipe + local_data(i-ips+1,j-jps+1) = local_data_1D( 1 + (i-ims) + ((j-jms) + (k-kms)*(jme-jms+1))*(ime-ims+1) ) + enddo + enddo + call gather_for_write(local_data,global_data,0) + if(me==0) then + do j=1,global_jlen + do i=1,global_ilen + global_data_1D(1 + (i-1) + ((j-1) + (k-kds)*global_jlen)*global_ilen) = & + global_data(i,j) + enddo + enddo + endif + enddo + + ! Ranks that are not writing data are done. + if(me/=0) then + ! This rank does not have the global data, and will not use the data anyway, + ! but it wants an array, so we'll send it the local data that is managed + ! internally by libIFI. + read_var => local_data_1D + should_dealloc = .false. + deallocate(global_data_1D) + else + read_var => global_data_1D + should_dealloc = .true. + endif + + deallocate(global_data) + deallocate(local_data) + end function read_var + end subroutine write_fip_output + + +end module upp_ifi_mod + + +#endif diff --git a/sorc/ncep_post.fd/INITPOST.F b/sorc/ncep_post.fd/INITPOST.F index 5e74564ab..5bde320f4 100644 --- a/sorc/ncep_post.fd/INITPOST.F +++ b/sorc/ncep_post.fd/INITPOST.F @@ -373,15 +373,15 @@ SUBROUTINE INITPOST do j = jsta_2l, jend_2u do i = 1, im if((PMID(I,J,ll-1) - PMID(I,J,ll))>=0.) then - write(*,*) 'non-monotonic PMID, i,j,ll ', i,j,ll - write(*,*) 'PMID: ll-1,ll,ll+1', PMID(I,J,LL-1) & - ,PMID(I,J,LL), PMID(I,J,LL+1) +! write(*,*) 'non-monotonic PMID, i,j,ll ', i,j,ll +! write(*,*) 'PMID: ll-1,ll,ll+1', PMID(I,J,LL-1) & +! ,PMID(I,J,LL), PMID(I,J,LL+1) PMID(I,J,LL)=0.5*(PMID(I,J,LL+1)+PMID(I,J,LL-1)) - write(*,*) 'after adjustment-i,j,ll ', i,j,ll - write(*,*) 'PMID: ll-1,ll,ll+1', PMID(I,J,LL-1) & - ,PMID(I,J,LL), PMID(I,J,LL+1) +! write(*,*) 'after adjustment-i,j,ll ', i,j,ll +! write(*,*) 'PMID: ll-1,ll,ll+1', PMID(I,J,LL-1) & +! ,PMID(I,J,LL), PMID(I,J,LL+1) endif end do end do @@ -400,16 +400,16 @@ SUBROUTINE INITPOST do j = jsta_2l, jend_2u do i = 1, im if((PMID(I,J,ll-1) - PMID(I,J,ll))>=0.) then - write(*,*) 'non-monotonic PMID, i,j,ll ', i,j,ll - write(*,*) 'PMID: ll-2,ll-1,ll', PMID(I,J,LL-2) & - ,PMID(I,J,LL-1), PMID(I,J,LL) + ! write(*,*) 'non-monotonic PMID, i,j,ll ', i,j,ll + ! write(*,*) 'PMID: ll-2,ll-1,ll', PMID(I,J,LL-2) & + ! ,PMID(I,J,LL-1), PMID(I,J,LL) PMID(I,J,LL)=PMID(I,J,ll-1) + & fact*(PMID(I,J,LL-1)-PMID(I,J,LL-2)) - write(*,*) 'after adjustment-i,j,ll ', i,j,ll - write(*,*) 'PMID: ll-2,ll-1,ll', PMID(I,J,LL-2) & - ,PMID(I,J,LL-1), PMID(I,J,LL) + ! write(*,*) 'after adjustment-i,j,ll ', i,j,ll + ! write(*,*) 'PMID: ll-2,ll-1,ll', PMID(I,J,LL-2) & + ! ,PMID(I,J,LL-1), PMID(I,J,LL) endif end do end do @@ -846,10 +846,10 @@ SUBROUTINE INITPOST DO J=Jsta,jend DO I=1,IM if((PINT(I,J,lm) - DUMMY(I,J))>=0.) then - write(*,*) 'non-monotonic PINT, i,j,lm ', i,j,lm - write(*,*) 'PINT: lm,lm+1, PMID: lm', PINT(I,J,LM),DUMMY(I,J), PMID(I,J,LM) + ! write(*,*) 'non-monotonic PINT, i,j,lm ', i,j,lm + ! write(*,*) 'PINT: lm,lm+1, PMID: lm', PINT(I,J,LM),DUMMY(I,J), PMID(I,J,LM) DUMMY(I,J)=PMID(I,J,LM)*1.001 - write(*,*) 'after adjustment-i,j,lm+1 ', i,j,lm+1,DUMMY(I,J) + ! write(*,*) 'after adjustment-i,j,lm+1 ', i,j,lm+1,DUMMY(I,J) endif PINT(I,J,LM+1)=DUMMY(I,J) ALPINT(I,J,LM+1)=ALOG(PINT(I,J,LM+1)) @@ -944,10 +944,10 @@ SUBROUTINE INITPOST ! KRF - check surface pressure for monotonic correctness if((PINT(I,J,lm) - PINT(I,J,LM+1))>=0. ) then - write(*,*) 'non-monotonic PINT, i,j,lm ', i,j,lm - write(*,*) 'PINT: lm,lm+1, PMID: lm', PINT(I,J,LM), PINT(I,J,LM+1), PMID(I,J,LM) + ! write(*,*) 'non-monotonic PINT, i,j,lm ', i,j,lm + ! write(*,*) 'PINT: lm,lm+1, PMID: lm', PINT(I,J,LM), PINT(I,J,LM+1), PMID(I,J,LM) PINT(I,J,LM+1) = PINT(I,J,LM)*1.001 - write(*,*) 'after adjustment-i,j,lm+1 PINT ', i,j,lm+1, PINT(I,J,LM+1) + ! write(*,*) 'after adjustment-i,j,lm+1 PINT ', i,j,lm+1, PINT(I,J,LM+1) endif ALPINT(I,J,LM+1)=ALOG(PINT(I,J,LM+1)) ENDDO diff --git a/sorc/ncep_post.fd/MDL2P.f b/sorc/ncep_post.fd/MDL2P.f index 8a023af0e..003f429b2 100644 --- a/sorc/ncep_post.fd/MDL2P.f +++ b/sorc/ncep_post.fd/MDL2P.f @@ -30,6 +30,7 @@ !> 2021-03-11 | B Cui | Change local arrays to dimension (im,jsta:jend) !> 2021-04-01 | J Meng | Computation on defined points only !> 2021-07-07 | J MENG | 2D DECOMPOSITION +!> 2022-08-03 | W Meng | Modify total cloud fraction(331) !> !> @author T Black W/NP2 @date 1999-09-23 SUBROUTINE MDL2P(iostatusD3D) @@ -718,7 +719,7 @@ SUBROUTINE MDL2P(iostatusD3D) FRIME(I,J) = 1. RAD(I,J) = 0. O3SL(I,J) = O3(I,J,LLMH) - CFRSL(I,J) = 0. + IF(CFR(I,J,1) SMALL) & - GRID1(I,J) = CFRSL(I,J)*H100 + IF(abs(CFRSL(I,J)-SPVAL) > SMALL) THEN + CFRSL(I,J) = MIN(MAX(0.0,CFRSL(I,J)),1.0) + GRID1(I,J) = CFRSL(I,J)*H100 + ENDIF ENDDO ENDDO if(grib == 'grib2')then diff --git a/sorc/ncep_post.fd/MDL2SIGMA.f b/sorc/ncep_post.fd/MDL2SIGMA.f index e2fd97a36..6dcb2654e 100644 --- a/sorc/ncep_post.fd/MDL2SIGMA.f +++ b/sorc/ncep_post.fd/MDL2SIGMA.f @@ -21,6 +21,7 @@ !! 20-03-25 J MENG - remove grib1 !! 21-03-11 B Cui - change local arrays to dimension (im,jsta:jend) !! 21-10-14 J MENG - 2D DECOMPOSITION +!! 2022-09-01 S Trahan - fixed bugs where extreme atmospheric conditions can cause out-of-bounds access !! !! USAGE: CALL MDL2P !! INPUT ARGUMENT LIST: @@ -277,8 +278,9 @@ SUBROUTINE MDL2SIGMA AHF =D00 FAC =D00 DONEFSL1=.TRUE. - ELSEIF(T(I,J,NL1XF(I,J))0 .or. IGET(1004)>0 .or. IGET(1005)>0 allocate(USHR1(ista_2l:iend_2u,jsta_2l:jend_2u),VSHR1(ista_2l:iend_2u,jsta_2l:jend_2u), & USHR6(ista_2l:iend_2u,jsta_2l:jend_2u),VSHR6(ista_2l:iend_2u,jsta_2l:jend_2u)) allocate(UST(ista_2l:iend_2u,jsta_2l:jend_2u),VST(ista_2l:iend_2u,jsta_2l:jend_2u), & HELI(ista_2l:iend_2u,jsta_2l:jend_2u,2),FSHR(ista_2l:iend_2u,jsta_2l:jend_2u)) + + print *,'start MISCLN' + ! ! HELICITY AND STORM MOTION. iget1 = IGET(162) @@ -860,7 +865,7 @@ SUBROUTINE MISCLN CALL FDLVL(ITYPEFDLVL,T7D,Q7D,U7D,V6D,P7D,ICINGFD,AERFD) ! - DO 10 IFD = 1,NFD + loop_10: DO IFD = 1,NFD ! ! FD LEVEL TEMPERATURE. iget1 = IGET(059) @@ -1328,7 +1333,7 @@ SUBROUTINE MISCLN ENDIF ENDIF - 10 CONTINUE + END DO loop_10 DEALLOCATE(T7D,Q7D,U7D,V6D,P7D,ICINGFD,AERFD) ENDIF @@ -1771,7 +1776,6 @@ SUBROUTINE MISCLN QBND(ista:iend,jsta:jend,NBND), UBND(ista:iend,jsta:jend,NBND), & VBND(ista:iend,jsta:jend,NBND), RHBND(ista:iend,jsta:jend,NBND), & WBND(ista:iend,jsta:jend,NBND)) - ! ! ***BLOCK 5: BOUNDARY LAYER FIELDS. ! @@ -1783,17 +1787,14 @@ SUBROUTINE MISCLN (IGET(090)>0).OR.(IGET(075)>0).OR. & (IGET(109)>0).OR.(IGET(110)>0).OR. & (IGET(031)>0).OR.(IGET(032)>0).OR. & - (IGET(573)>0).OR. & + (IGET(573)>0).OR. NEED_IFI .OR. & (IGET(107)>0).OR.(IGET(091)>0).OR. & (IGET(092)>0).OR.(IGET(093)>0).OR. & (IGET(094)>0).OR.(IGET(095)>0).OR. & (IGET(096)>0).OR.(IGET(097)>0).OR. & (IGET(098)>0).OR.(IGET(221)>0) ) THEN ! - allocate(OMGBND(ista:iend,jsta:jend,NBND),PWTBND(ista:iend,jsta:jend,NBND), & - QCNVBND(ista:iend,jsta:jend,NBND),LVLBND(ista:iend,jsta:jend,NBND), & - LB2(ista:iend,jsta:jend)) - + call allocate_cape_arrays ! COMPUTE ETA BOUNDARY LAYER FIELDS. CALL BNDLYR(PBND,TBND,QBND,RHBND,UBND,VBND, & WBND,OMGBND,PWTBND,QCNVBND,LVLBND) @@ -1807,7 +1808,7 @@ SUBROUTINE MISCLN ! ! LOOP OVER NBND BOUNDARY LAYERS. - DO 20 LBND = 1,NBND + boundary_layer_loop: DO LBND = 1,NBND ! ! BOUNDARY LAYER PRESSURE. IF (IGET(067)>0) THEN @@ -2118,7 +2119,7 @@ SUBROUTINE MISCLN ENDIF ! ! END OF ETA BOUNDARY LAYER LOOP. - 20 CONTINUE + END DO boundary_layer_loop deallocate(OMGBND,PWTBND,QCNVBND) ! ! BEST LIFTED INDEX FROM BOUNDARY LAYER FIELDS. @@ -2197,8 +2198,9 @@ SUBROUTINE MISCLN ! LVLSXML(1,IGET(566)),'LVLSXML(1,IGET(567)=', & ! LVLSXML(1,IGET(567)),'field1=',field1,'field2=',field2 ! - IF(FIELD1.OR.FIELD2)THEN + IF(FIELD1.OR.FIELD2.OR.NEED_IFI)THEN ITYPE = 2 + call allocate_cape_arrays ! !$omp parallel do private(i,j) DO J=JSTA,JEND @@ -2208,7 +2210,7 @@ SUBROUTINE MISCLN ENDDO ENDDO ! - DO 80 LBND = 1,NBND + loop_80: DO LBND = 1,NBND CALL CALTHTE(PBND(ista,jsta,LBND),TBND(ista,jsta,LBND), & QBND(ista,jsta,LBND),EGRID1) !$omp parallel do private(i,j) @@ -2223,14 +2225,14 @@ SUBROUTINE MISCLN ENDIF ENDDO ENDDO - 80 CONTINUE + ENDDO loop_80 ! DPBND = 0. CALL CALCAPE(ITYPE,DPBND,P1D,T1D,Q1D,LB2,EGRID1, & EGRID2,EGRID3,EGRID4,EGRID5) + ! - IF (IGET(566)>0) THEN -! dong add missing value for cape + IF(IGET(566)>0 .or. NEED_IFI) THEN GRID1=spval !$omp parallel do private(i,j) DO J=JSTA,JEND @@ -2238,7 +2240,30 @@ SUBROUTINE MISCLN IF(T1D(I,J) < spval) GRID1(I,J) = EGRID1(I,J) ENDDO ENDDO + CALL BOUND(GRID1,D00,H99999) +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + CAPE(I,J) = GRID1(I,J) + ENDDO + ENDDO + ENDIF + + IF(IGET(567)>0 .or. NEED_IFI) THEN + GRID1=spval +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IF(T1D(I,J) < spval) THEN + GRID1(I,J) = - EGRID2(I,J) + ENDIF + ENDDO + ENDDO +! CALL BOUND(GRID1,D00,H99999) + ENDIF + + IF (IGET(566)>0) THEN if(grib=='grib2') then cfld=cfld+1 fld_info(cfld)%ifld=IAVBLFLD(IGET(566)) @@ -2248,14 +2273,14 @@ SUBROUTINE MISCLN jj = jsta+j-1 do i=1,iend-ista+1 ii = ista+i-1 - datapd(i,j,cfld) = GRID1(ii,jj) + datapd(i,j,cfld) = CAPE(ii,jj) enddo enddo endif ENDIF ! - IF (IGET(567) > 0) THEN -! dong add missing value for cape + IF (IGET(567) > 0 .or. NEED_IFI) THEN +! dong add missing value for CIN GRID1=spval !$omp parallel do private(i,j) DO J=JSTA,JEND @@ -2265,14 +2290,18 @@ SUBROUTINE MISCLN ENDDO ! CALL BOUND(GRID1,D00,H99999) + print *,"STORE CIN" ! !$omp parallel do private(i,j) DO J=JSTA,JEND DO I=ISTA,IEND IF(T1D(I,J) < spval) GRID1(I,J) = - GRID1(I,J) + CIN(I,J) = GRID1(I,J) ENDDO ENDDO -! + ENDIF + + IF(IGET(567) > 0) THEN if(grib=='grib2') then cfld=cfld+1 fld_info(cfld)%ifld=IAVBLFLD(IGET(567)) @@ -2282,7 +2311,7 @@ SUBROUTINE MISCLN jj = jsta+j-1 do i=1,iend-ista+1 ii = ista+i-1 - datapd(i,j,cfld) = GRID1(ii,jj) + datapd(i,j,cfld) = CIN(ii,jj) enddo enddo endif @@ -3365,8 +3394,9 @@ SUBROUTINE MISCLN FIELD2=.TRUE. ENDIF ! - IF(FIELD1.OR.FIELD2)THEN + IF(FIELD1.OR.FIELD2.OR.NEED_IFI)THEN ITYPE = 1 + call allocate_cape_arrays ! !$omp parallel do private(i,j) DO J=JSTA,JEND @@ -3379,15 +3409,15 @@ SUBROUTINE MISCLN DPBND = 300.E2 CALL CALCAPE(ITYPE,DPBND,P1D,T1D,Q1D,LB2,EGRID1, & EGRID2,EGRID3,EGRID4,EGRID5) - IF (IGET(584)>0) THEN + IF (IGET(584)>0 .or. NEED_IFI) THEN ! dong add missing value to cin GRID1 = spval !$omp parallel do private(i,j) DO J=JSTA,JEND DO I=ISTA,IEND IF(T1D(I,J) < spval) THEN - GRID1(I,J) = EGRID1(I,J) - IF (SUBMODELNAME == 'RTMA') MUCAPE(I,J)=GRID1(I,J) + GRID1(I,J) = EGRID1(I,J) + IF (SUBMODELNAME == 'RTMA') MUCAPE(I,J)=GRID1(I,J) ENDIF ENDDO ENDDO @@ -3395,7 +3425,13 @@ SUBROUTINE MISCLN ! IF (SUBMODELNAME == 'RTMA') THEN ! CALL BOUND(MUCAPE,D00,H99999) ! ENDIF - if(grib=='grib2') then +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + CAPE(I,J) = GRID1(I,J) + ENDDO + ENDDO + if(IGET(584)>0 .and. grib=='grib2') then cfld=cfld+1 fld_info(cfld)%ifld=IAVBLFLD(IGET(584)) fld_info(cfld)%lvl=LVLSXML(1,IGET(584)) @@ -3411,7 +3447,7 @@ SUBROUTINE MISCLN ENDIF - IF (IGET(585)>0) THEN + IF (IGET(585)>0 .or. NEED_IFI) THEN ! dong add missing value to cin GRID1 = spval !$omp parallel do private(i,j) @@ -3432,7 +3468,15 @@ SUBROUTINE MISCLN ENDIF ENDDO ENDDO - if(grib=='grib2') then + +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + CIN(I,J) = GRID1(I,J) + ENDDO + ENDDO + + if(IGET(585)>0 .and. grib=='grib2') then cfld=cfld+1 fld_info(cfld)%ifld=IAVBLFLD(IGET(585)) fld_info(cfld)%lvl=LVLSXML(1,IGET(585)) @@ -3707,6 +3751,8 @@ SUBROUTINE MISCLN IF(FIELD1.OR.FIELD2)THEN ITYPE = 2 + call allocate_cape_arrays + ! !$omp parallel do private(i,j) DO J=JSTA,JEND @@ -4681,4 +4727,14 @@ SUBROUTINE MISCLN ! END OF ROUTINE. ! RETURN - END + CONTAINS + + subroutine allocate_cape_arrays + if(.not.allocated(OMGBND)) allocate(OMGBND(ista:iend,jsta:jend,NBND)) + if(.not.allocated(PWTBND)) allocate(PWTBND(ista:iend,jsta:jend,NBND)) + if(.not.allocated(QCNVBND)) allocate(QCNVBND(ista:iend,jsta:jend,NBND)) + if(.not.allocated(LVLBND)) allocate(LVLBND(ista:iend,jsta:jend,NBND)) + if(.not.allocated(LB2)) allocate(LB2(ista:iend,jsta:jend)) + end subroutine allocate_cape_arrays + + END SUBROUTINE MISCLN diff --git a/sorc/ncep_post.fd/PROCESS.f b/sorc/ncep_post.fd/PROCESS.f index 64c9c35de..a969dac3d 100644 --- a/sorc/ncep_post.fd/PROCESS.f +++ b/sorc/ncep_post.fd/PROCESS.f @@ -30,10 +30,10 @@ SUBROUTINE PROCESS(kth,kpv,th,pv,iostatusD3D) !---------------------------------------------------------------------------- ! use mpi, only: mpi_wtime - + use upp_ifi_mod, only: run_ifi use CTLBLK_mod, only: cfld, etafld2_tim, eta2p_tim, mdl2sigma_tim, surfce2_tim,& mdl2agl_tim, mdl2std_tim, mdl2thandpv_tim, calrad_wcloud_tim,& - cldrad_tim, miscln_tim, fixed_tim, ntlfld, me + cldrad_tim, miscln_tim, fixed_tim, ntlfld, me, run_ifi_tim !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! @@ -128,6 +128,11 @@ SUBROUTINE PROCESS(kth,kpv,th,pv,iostatusD3D) if(me==0) write(0,*) "PROCESS CALRAD_WCLOUD done" CALRAD_WCLOUD_tim = CALRAD_WCLOUD_tim +(mpi_wtime() - btim) ! +! IN-FLIGHT ICING PRODUCTS + btim = mpi_wtime() + CALL RUN_IFI + RUN_IFI_tim = RUN_IFI_tim +(mpi_wtime()-btim) +! ! END OF ROUTINE. ! NTLFLD=cfld diff --git a/sorc/ncep_post.fd/SET_LVLSXML.f b/sorc/ncep_post.fd/SET_LVLSXML.f index ed51e12d1..c29c97488 100644 --- a/sorc/ncep_post.fd/SET_LVLSXML.f +++ b/sorc/ncep_post.fd/SET_LVLSXML.f @@ -47,7 +47,7 @@ subroutine SET_LVLSXML(param,ifld,irec,kpv,pv,kth,th) ! use xml_perl_data, only: param_t use ctlblk_mod, only: lsm, spl, nsoil, isf_surface_physics, nfd, htfd, & - petabnd, nbnd + petabnd, nbnd, ifi_nflight, ifi_flight_levels use soil, only: SLDPTH,SLLEVEL use rqstfld_mod,only : mxlvl,LVLS,LVLSXML implicit none @@ -65,7 +65,7 @@ subroutine SET_LVLSXML(param,ifld,irec,kpv,pv,kth,th) real,parameter :: small2=1 integer,parameter :: LSIG1=22,LSIG2=5 integer i,j,l,nlevel,scalef,lvlcape,lvlcin - logical READTHK,logrec + logical READTHK,logrec,found REAL :: SIGO2(LSIG2+1),ASIGO2(LSIG2),DSIGO2(LSIG2) REAL :: SIGO1(LSIG1+1),ASIGO1(LSIG1),DSIGO1(LSIG1) ! @@ -187,23 +187,48 @@ subroutine SET_LVLSXML(param,ifld,irec,kpv,pv,kth,th) endif ! if(trim(param%fixed_sfc1_type)=='spec_alt_above_mean_sea_lvl') then - if(index(param%shortname,"GTG_ON_SPEC_ALT_ABOVE_MEAN_SEA_LVL")<=0) then - do j=1, nlevel + if(index(param%shortname,"SPECIFIC_IFI_FLIGHT_LEVEL")>0) then + do j=1, nlevel + found=.false. + iloop411: do i=1, ifi_nflight + if(nint(param%level(j)/10)==nint(ifi_flight_levels(i)/10) )then + LVLS(i,ifld)=1 + LVLSXML(i,ifld)=j + !print *,'SPECIFIC_IFI_FLIGHT_LEVEL ',j,' is ',param%level(j) + irec=irec+1 + found=.true. + exit iloop411 + endif + enddo iloop411 + if(.not.found) then + write(0,*) 'ERROR: No such IFI flight level: ',param%level(j)/10 + LVLS(i,ifld)=0 + endif + enddo + else if(index(param%shortname,"IFI_FLIGHT_LEVEL")>0) then + do j=1, ifi_nflight + LVLS(j,ifld)=1 + LVLSXML(j,ifld)=j + !print *,'IFI_FLIGHT_LEVEL ',j,' is ',param%level(j) + irec=irec+1 + enddo + elseif(index(param%shortname,"GTG_ON_SPEC_ALT_ABOVE_MEAN_SEA_LVL")<=0) then + do j=1, nlevel iloop4: do i=1, NFD - if(nint(param%level(j))==nint(HTFD(i)) )then - if(HTFD(i)>300.) then - LVLS(i,ifld)=1 - else - LVLS(i,ifld)=2 + if(nint(param%level(j))==nint(HTFD(i)) )then + if(HTFD(i)>300.) then + LVLS(i,ifld)=1 + else + LVLS(i,ifld)=2 + endif + LVLSXML(i,ifld)=j + irec=irec+1 + exit iloop4 endif - LVLSXML(i,ifld)=j - irec=irec+1 - exit iloop4 - endif - enddo iloop4 - enddo - return - endif + enddo iloop4 + enddo + return + endif ! Allow inputs of FD levels from control file. For GTG (EDPARM CATEDR MWTURB) ! SET LVLS to 1 do j=1, nlevel @@ -277,15 +302,16 @@ subroutine SET_LVLSXML(param,ifld,irec,kpv,pv,kth,th) endif enddo iloop41 enddo - return + return endif do j=1, nlevel - LVLS(j,ifld)=1 - LVLSXML(j,ifld)=j - irec=irec+1 + LVLS(j,ifld)=1 + LVLSXML(j,ifld)=j + irec=irec+1 enddo return endif +! !for hpc tmp at sigma lvl if(trim(param%shortname)=='TMP_ON_SIGMA_LVL_HPC') then IF(READTHK)THEN ! EITHER READ DSG THICKNESS diff --git a/sorc/ncep_post.fd/SURFCE.f b/sorc/ncep_post.fd/SURFCE.f index cb6260612..d62d7583a 100644 --- a/sorc/ncep_post.fd/SURFCE.f +++ b/sorc/ncep_post.fd/SURFCE.f @@ -91,7 +91,7 @@ SUBROUTINE SURFCE AVGCPRATE_CONT,sst,pcp_bucket1,rainnc_bucket1, & snow_bucket1, rainc_bucket1, graup_bucket1, & shdmin, shdmax, lai, ch10,cd10,landfrac,paha,pahi, & - tecan,tetran,tedir,twa + tecan,tetran,tedir,twa,IFI_APCP use soil, only: stc, sllevel, sldpth, smc, sh2o use masks, only: lmh, sm, sice, htm, gdlat, gdlon use physcons_post,only: CON_EPS, CON_EPSM1 @@ -3689,17 +3689,21 @@ SUBROUTINE SURFCE ! PRECIPITATION BUCKETS - accumulated between output times ! 'BUCKET TOTAL PRECIP ' - IF (IGET(434)>0.) THEN + IF (IGET(434)>0. .or. IGET(1003)>0 .or. IGET(1004)>0 .or. IGET(1005)>0) THEN !$omp parallel do private(i,j) DO J=JSTA,JEND DO I=ISTA,IEND IF (IFHR == 0) THEN - GRID1(I,J) = 0.0 + IFI_APCP(I,J) = 0.0 ELSE - GRID1(I,J) = PCP_BUCKET(I,J) + IFI_APCP(I,J) = PCP_BUCKET(I,J) ENDIF ENDDO ENDDO + ! Note: IFI.F may replace IFI_APCP with other values where it is spval or 0 + ENDIF + + IF (IGET(434)>0.) THEN ID(1:25) = 0 ITPREC = NINT(TPREC) !mp @@ -3722,7 +3726,7 @@ SUBROUTINE SURFCE IF(IFMIN >= 1)ID(18)=IFHR*60+IFMIN-IFINCR ENDIF IF (ID(18)<0) ID(18) = 0 - if(grib=='grib2') then + if(grib=='grib2' .and. IGET(434)>0) then cfld=cfld+1 fld_info(cfld)%ifld=IAVBLFLD(IGET(434)) if(ITPREC>0) then @@ -3744,7 +3748,7 @@ SUBROUTINE SURFCE jj = jsta+j-1 do i=1,iend-ista+1 ii = ista+i-1 - datapd(i,j,cfld) = GRID1(ii,jj) + datapd(i,j,cfld) = IFI_APCP(i,jj) enddo enddo endif diff --git a/sorc/ncep_post.fd/UPP_PHYSICS.f b/sorc/ncep_post.fd/UPP_PHYSICS.f index cc609aabf..7367b3d15 100644 --- a/sorc/ncep_post.fd/UPP_PHYSICS.f +++ b/sorc/ncep_post.fd/UPP_PHYSICS.f @@ -998,6 +998,8 @@ END SUBROUTINE CALCAPE !> 2015-??-?? | S Moorthi | Optimization and threading !> 2021-09-03 | J Meng | Modified to add 0-3km CAPE/CINS, LFC, effective helicity, downdraft CAPE, dendritic growth layer depth, ESP !> 2021-09-01 | E Colon | Equivalent level height index for RTMA +!> 2022-08-27 | S Trahan | Fixed bug in CALCAPE2 where extreme atmospheric conditions cause an out-of-bounds access +!> 2022-09-01 | S Trahan | Fixed another bug in CALCAPE2 where extreme atmospheric conditions cause an out-of-bounds access !> !> @author Russ Treadon W/NP2 @date 1993-02-10 SUBROUTINE CALCAPE2(ITYPE,DPBND,P1D,T1D,Q1D,L1D, & @@ -1319,6 +1321,9 @@ SUBROUTINE CALCAPE2(ITYPE,DPBND,P1D,T1D,Q1D,L1D, & THUNDER(I,J) = .FALSE. ENDIF ENDIF + + ! Limit LCL to prevent out-of-bounds accesses later + LCL(I,J) = max(min(LCL(I,J),LM-1),1) ENDDO ENDDO !----------------------------------------------------------------------- @@ -1398,6 +1403,13 @@ SUBROUTINE CALCAPE2(ITYPE,DPBND,P1D,T1D,Q1D,L1D, & ENDIF ENDDO ENDDO + + LEND=MIN(LEND,LM-1) + +! +!Ensure later calculations do not access LM+1 +! + LEND=MIN(LEND,LM-1) ! !reverse L order from bottom up for ESRH calculation ! diff --git a/sorc/ncep_post.fd/VRBLS2D_mod.f b/sorc/ncep_post.fd/VRBLS2D_mod.f index 569f34ea5..e17ec8009 100644 --- a/sorc/ncep_post.fd/VRBLS2D_mod.f +++ b/sorc/ncep_post.fd/VRBLS2D_mod.f @@ -29,6 +29,8 @@ module vrbls2d ,CLDEFI(:,:),ALBASE(:,:),SI(:,:),LSPA(:,:) & ,RSWINC(:,:),VIS(:,:),PD(:,:),MXSNAL(:,:),MIXHT(:,:) & ,SNONC(:,:),EPSR(:,:),RSWTOA(:,:),TEQL(:,:) & +! Variables saved for input to IFI + ,IFI_APCP(:,:),CAPE(:,:),CIN(:,:) & ! HWRF additions ,MDLTAUX(:,:),MDLTAUY(:,:),CD10(:,:),CH10(:,:) & ,ACSWUPT(:,:),SWDNT(:,:),ACSWDNT(:,:) & diff --git a/sorc/ncep_post.fd/VRBLS3D_mod.f b/sorc/ncep_post.fd/VRBLS3D_mod.f index f27428f3a..bf98549a7 100644 --- a/sorc/ncep_post.fd/VRBLS3D_mod.f +++ b/sorc/ncep_post.fd/VRBLS3D_mod.f @@ -80,7 +80,6 @@ module vrbls3d ,icing_gfip(:,:,:),icing_gfis(:,:,:) & ! Add NCAR GTG turbulence ,catedr(:,:,:),mwt(:,:,:),gtg(:,:,:) & - ! AQF ,ozcon(:,:,:),pmtf(:,:,:) diff --git a/sorc/ncep_post.fd/WRFPOST.f b/sorc/ncep_post.fd/WRFPOST.f index c670150f9..a5633e421 100644 --- a/sorc/ncep_post.fd/WRFPOST.f +++ b/sorc/ncep_post.fd/WRFPOST.f @@ -115,8 +115,10 @@ PROGRAM WRFPOST lsm, fld_info, etafld2_tim, eta2p_tim, mdl2sigma_tim, cldrad_tim, miscln_tim, & mdl2agl_tim, mdl2std_tim, mdl2thandpv_tim, calrad_wcloud_tim, & fixed_tim, time_output, imin, surfce2_tim, komax, ivegsrc, d3d_on, gocart_on,rdaod, & - readxml_tim, spval, fullmodelname, submodelname, hyb_sigp, filenameflat, aqfcmaq_on,numx + readxml_tim, spval, fullmodelname, submodelname, hyb_sigp, filenameflat, aqfcmaq_on, & + numx, run_ifi_tim use grib2_module, only: gribit2,num_pset,nrecout,first_grbtbl,grib_info_finalize + use upp_ifi_mod, only: write_ifi_debug_files, enable_ifi_smoother !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! @@ -141,7 +143,8 @@ PROGRAM WRFPOST integer :: kpo,kth,kpv real,dimension(komax) :: po,th,pv namelist/nampgb/kpo,po,kth,th,kpv,pv,fileNameAER,d3d_on,gocart_on,popascal & - ,hyb_sigp,rdaod,aqfcmaq_on,vtimeunits,numx + ,hyb_sigp,rdaod,aqfcmaq_on,vtimeunits,numx & + ,write_ifi_debug_files,enable_ifi_smoother integer :: itag_ierr namelist/model_inputs/fileName,IOFORM,grib,DateStr,MODELNAME,SUBMODELNAME & ,fileNameFlux,fileNameFlat @@ -270,6 +273,7 @@ PROGRAM WRFPOST fileNameFlat='postxconfig-NT.txt' read(5,nampgb,iostat=iret,end=119) 119 continue + if (me==0) print*,'in itag, write_ifi_debug_files=', write_ifi_debug_files if (me==0) print*,'in itag, mod(num_procs,numx)=', mod(num_procs,numx) if(mod(num_procs,numx)/=0) then if (me==0) then @@ -721,7 +725,7 @@ PROGRAM WRFPOST CALL PROCESS(kth,kpv,th(1:kth),pv(1:kpv),iostatusD3D) IF(ME == 0) WRITE(6,*)'WRFPOST: PREPARE TO PROCESS NEXT GRID' ! -! write(0,*)'enter gribit2 before mpi_barrier' +! write(0,*)'enter mpi_barrier before gribit2' call mpi_barrier(mpi_comm_comp,ierr) ! if(me==0)call w3tage('bf grb2 ') @@ -764,6 +768,7 @@ PROGRAM WRFPOST print*, 'FIXED_tim = ',FIXED_tim print*, 'MDL2THANDPV_tim = ',MDL2THANDPV_tim print*, 'CALRAD_WCLOUD_tim = ',CALRAD_WCLOUD_tim + print*, 'RUN_IFI_tim = ',RUN_IFI_tim print*, 'Total time = ',(mpi_wtime() - bbtim) print*, 'Time for OUTPUT = ',time_output print*, 'Time for READxml = ',READxml_tim diff --git a/sorc/ncep_post.fd/grib2_module.f b/sorc/ncep_post.fd/grib2_module.f index bb24fc660..fadd07ec6 100644 --- a/sorc/ncep_post.fd/grib2_module.f +++ b/sorc/ncep_post.fd/grib2_module.f @@ -225,6 +225,7 @@ subroutine gribit2(post_fname) logical, parameter :: debugprint = .false. ! character(1), dimension(:), allocatable :: cgrib + real :: level_unit_conversion ! ! !---------------- code starts here -------------------------- @@ -324,10 +325,14 @@ subroutine gribit2(post_fname) ' category ',icatg, & ' parameter ',iparm, & ' for var ',trim(pset%param(nprm)%pname) - + if(index(pset%param(nprm)%shortname,'IFI_FLIGHT_LEVEL')>0) then + level_unit_conversion=0.3048 ! convert feet->meters + else + level_unit_conversion=1 + endif call gengrb2msg(idisc,icatg, iparm,nprm,nlvl,fldlvl1,fldlvl2, & fld_info(i)%ntrange,fld_info(i)%tinvstat,datafld(:,i), & - cgrib,clength) + cgrib,clength,level_unit_conversion) ! print *,'finished gengrb2msg field=',i,'ntlfld=',ntlfld,'clength=',clength call wryte(lunout, clength, cgrib) else @@ -419,8 +424,14 @@ subroutine gribit2(post_fname) ! !--- generate grib2 message --- ! + if(index(pset%param(nprm)%shortname,'IFI_FLIGHT_LEVEL')>0) then + level_unit_conversion=0.3048 ! convert feet->meters + else + level_unit_conversion=1 + endif call gengrb2msg(idisc,icatg, iparm,nprm,nlvl,fldlvl1,fldlvl2,ntrange, & - leng_time_range_stat,datafld(:,i),cgrib(cstart),clength) + leng_time_range_stat,datafld(:,i),cgrib(cstart),clength, & + level_unit_conversion) cstart=cstart+clength ! else @@ -479,7 +490,7 @@ end subroutine gribit2 !---------------------------------------------------------------------------------------- ! subroutine gengrb2msg(idisc,icatg, iparm,nprm,nlvl,fldlvl1,fldlvl2,ntrange,tinvstat, & - datafld1,cgrib,lengrib) + datafld1,cgrib,lengrib,level_unit_conversion) ! !---------------------------------------------------------------------------------------- ! @@ -498,6 +509,7 @@ subroutine gengrb2msg(idisc,icatg, iparm,nprm,nlvl,fldlvl1,fldlvl2,ntrange,tinvs real,dimension(:),intent(in) :: datafld1 character(1),intent(inout) :: cgrib(max_bytes) integer, intent(inout) :: lengrib + real, intent(in) :: level_unit_conversion ! integer, parameter :: igdsmaxlen=200 ! @@ -734,6 +746,14 @@ subroutine gengrb2msg(idisc,icatg, iparm,nprm,nlvl,fldlvl1,fldlvl2,ntrange,tinvs scale_fct_fixed_sfc2=0 endif + if(abs(level_unit_conversion-1)>1e-4) then + print *,'apply level unit conversion ',level_unit_conversion + print *,'scaled_val_fixed_sfc1 was ',scaled_val_fixed_sfc1 + scaled_val_fixed_sfc1=nint(scaled_val_fixed_sfc1*real(level_unit_conversion,kind=8)) + scaled_val_fixed_sfc2=nint(scaled_val_fixed_sfc2*real(level_unit_conversion,kind=8)) + print *,'scaled_val_fixed_sfc1 now ',scaled_val_fixed_sfc1 + endif + ihr_start = ifhr-tinvstat if(modelname=='RAPR'.and.vtimeunits=='FMIN') then ifhrorig = ifhr diff --git a/tests/compile_upp.sh b/tests/compile_upp.sh index 2c20f660c..90d6133c1 100755 --- a/tests/compile_upp.sh +++ b/tests/compile_upp.sh @@ -7,11 +7,13 @@ set -eu usage() { echo - echo "Usage: $0 [-p] [-g] [-w] [-v] [-c] -h" + echo "Usage: $0 [-p] [-g] [-w] [-v] [-c] [-i] -h" echo echo " -p installation prefix DEFAULT: ../install" echo " -g build with GTG(users with gtg repos. access only) DEFAULT: OFF" + echo " -i build with libIFI (users with ifi access only) DEFAULT: OFF" echo " -w build without WRF-IO DEFAULT: ON" + echo " -i build with IFI DEFAULT: OFF" echo " -v build with cmake verbose DEFAULT: NO" echo " -c Compiler to use for build DEFAULT: intel" echo " -h display this message and quit" @@ -20,11 +22,13 @@ usage() { } prefix="../install" +ifi_opt=" -DBUILD_WITH_IFI=OFF" gtg_opt=" -DBUILD_WITH_GTG=OFF" wrfio_opt=" -DBUILD_WITH_WRFIO=ON" +ifi_opt=" -DBUILD_WITH_IFI=OFF" compiler="intel" verbose_opt="" -while getopts ":p:gwc:vh" opt; do +while getopts ":p:gwc:vhi" opt; do case $opt in p) prefix=$OPTARG @@ -32,9 +36,15 @@ while getopts ":p:gwc:vh" opt; do g) gtg_opt=" -DBUILD_WITH_GTG=ON" ;; + i) + ifi_opt=" -DREQUIRE_IFI=ON" + ;; w) wrfio_opt=" -DBUILD_WITH_WRFIO=OFF" ;; + i) + ifi_opt=" -DREQUIRE_IFI=ON" + ;; c) compiler=$OPTARG ;; @@ -46,7 +56,7 @@ while getopts ":p:gwc:vh" opt; do ;; esac done -cmake_opts=" -DCMAKE_INSTALL_PREFIX=$prefix"${wrfio_opt}${gtg_opt} +cmake_opts=" -DCMAKE_INSTALL_PREFIX=$prefix"${wrfio_opt}${gtg_opt}${ifi_opt} source ./detect_machine.sh if [[ $(uname -s) == Darwin ]]; then