diff --git a/CMakeLists.txt b/CMakeLists.txt index f1740d46e..051a4bdf0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -40,6 +40,7 @@ find_package(Jasper REQUIRED) find_package(PNG REQUIRED) find_package(ESMF MODULE REQUIRED) find_package(ZLIB REQUIRED) +find_package(WGRIB2 REQUIRED) if(NOT TARGET landsfcutil_4) find_package(landsfcutil REQUIRED) diff --git a/cmake b/cmake index 7f5b1b3f5..d6d834fa1 160000 --- a/cmake +++ b/cmake @@ -1 +1 @@ -Subproject commit 7f5b1b3f51ddaf2d33b0e658587b978887c524f7 +Subproject commit d6d834fa1c2023248308643cffb4bb08fd3a9856 diff --git a/modulefiles/chgres_cube.hera b/modulefiles/chgres_cube.hera index 920d1efff..8fbd5ba4c 100644 --- a/modulefiles/chgres_cube.hera +++ b/modulefiles/chgres_cube.hera @@ -19,6 +19,9 @@ module load sfcio//v1.1.0-intel-19.0.4.243 module load sigio/v2.1.0-intel-19.0.4.243-impi module load esmf/8.0.0-intel-19.0.4.243-impi +export WGRIB2API_INC="/apps/wgrib2/2.0.8/intel/18.0.3.222/lib" +export WGRIB2_LIB="/apps/wgrib2/2.0.8/intel/18.0.3.222/lib/libwgrib2.a" + export FCOMP=mpiifort export FFLAGS="-O3 -fp-model precise -g -traceback -r8 -i4 -qopenmp -convert big_endian -assume byterecl" # for debugging diff --git a/modulefiles/chgres_cube.jet b/modulefiles/chgres_cube.jet index 9ec1069dc..6a1f9340e 100644 --- a/modulefiles/chgres_cube.jet +++ b/modulefiles/chgres_cube.jet @@ -4,8 +4,8 @@ module load intel/18.0.5.274 module load impi/2018.4.274 -module load szip -module load hdf5 +module load szip/2.1 +module load hdf5/1.8.9 module load netcdf/4.2.1.1 module load w3nco/v2.0.6 @@ -22,3 +22,6 @@ export FCOMP=mpiifort export FFLAGS="-O3 -fp-model precise -g -traceback -r8 -i4 -qopenmp -convert big_endian -assume byterecl" #debug #export FFLAGS="-O0 -g -traceback -r8 -i4 -qopenmp -convert big_endian -check bounds -warn unused -assume byterecl" + +export WGRIB2API_INC="/mnt/lfs3/projects/hwrfv3/Jili.Dong/wgrib2-2.0.8/grib2/lib" +export WGRIB2_LIB="/mnt/lfs3/projects/hwrfv3/Jili.Dong/wgrib2-2.0.8/grib2/lib/libwgrib2.a" diff --git a/modulefiles/chgres_cube.linux.gnu b/modulefiles/chgres_cube.linux.gnu new file mode 100644 index 000000000..ba76044c2 --- /dev/null +++ b/modulefiles/chgres_cube.linux.gnu @@ -0,0 +1,24 @@ +############################################################# +## chgres_cube component - linux.gnu +############################################################# + +export IP_INCd=${NCEPLIBS}/ip/include_d +export NEMSIO_INC=${NCEPLIBS}/nemsio/include +export SFCIO_INC4=${NCEPLIBS}/sfcio/include_4 +export SIGIO_INC4=${NCEPLIBS}/sigio/include_4 + +export BACIO_LIB4=${NCEPLIBS}/bacio/lib/libbacio_v2.1.0_4.a +export IP_LIBd=${NCEPLIBS}/ip/lib/libip_v3.0.0_d.a +export NEMSIO_LIB=${NCEPLIBS}/nemsio/lib/libnemsio_v2.2.3.a +export SFCIO_LIB4=${NCEPLIBS}/sfcio/lib/libsfcio_v1.1.0_4.a +export SIGIO_LIB4=${NCEPLIBS}/sigio/lib/libsigio_v2.1.0_4.a +export SP_LIBd=${NCEPLIBS}/sp/lib/libsp_v2.0.2_d.a +export W3NCO_LIBd=${NCEPLIBS}/w3nco/lib/libw3nco_v2.0.6_d.a + +export WGRIB2API_INC=${WGRIB2_DIR}/include +export WGRIB2_LIB=${WGRIB2_DIR}/lib/libwgrib2.a + +export FCOMP=mpif90 +export FFLAGS="-O3 -g -fbacktrace -fdefault-real-8 -ffree-line-length-none -fopenmp -fconvert=big-endian" +# for debugging +#export FFLAGS="-O0 -g -fbacktrace -fdefault-real-8 -ffree-line-length-none -fopenmp -fconvert=big-endian" diff --git a/modulefiles/chgres_cube.linux.intel b/modulefiles/chgres_cube.linux.intel new file mode 100644 index 000000000..9af64a3b9 --- /dev/null +++ b/modulefiles/chgres_cube.linux.intel @@ -0,0 +1,24 @@ +############################################################# +## chgres_cube component - linux.intel +############################################################# + +export IP_INCd=${NCEPLIBS}/ip/include_d +export NEMSIO_INC=${NCEPLIBS}/nemsio/include +export SFCIO_INC4=${NCEPLIBS}/sfcio/include_4 +export SIGIO_INC4=${NCEPLIBS}/sigio/include_4 + +export BACIO_LIB4=${NCEPLIBS}/bacio/lib/libbacio_v2.1.0_4.a +export IP_LIBd=${NCEPLIBS}/ip/lib/libip_v3.0.0_d.a +export NEMSIO_LIB=${NCEPLIBS}/nemsio/lib/libnemsio_v2.2.3.a +export SFCIO_LIB4=${NCEPLIBS}/sfcio/lib/libsfcio_v1.1.0_4.a +export SIGIO_LIB4=${NCEPLIBS}/sigio/lib/libsigio_v2.1.0_4.a +export SP_LIBd=${NCEPLIBS}/sp/lib/libsp_v2.0.2_d.a +export W3NCO_LIBd=${NCEPLIBS}/w3nco/lib/libw3nco_v2.0.6_d.a + +export WGRIB2API_INC=${WGRIB2_DIR}/include +export WGRIB2_LIB=${WGRIB2_DIR}/lib/libwgrib2.a + +export FCOMP=mpif90 +export FFLAGS="-O3 -fp-model source -g -traceback -r8 -i4 -qopenmp -convert big_endian -assume byterecl" +# for debugging +#export FFLAGS="-O0 -g -traceback -r8 -i4 -qopenmp -convert big_endian -check bounds -warn unused -assume byterecl" diff --git a/modulefiles/chgres_cube.odin b/modulefiles/chgres_cube.odin new file mode 100644 index 000000000..ecbc18977 --- /dev/null +++ b/modulefiles/chgres_cube.odin @@ -0,0 +1,27 @@ +#%Module##################################################### +## chgres build module for Odin +############################################################# + +module use /oldscratch/ywang/external/modulefiles +module load esmf/8.0.0bs30 + +module load cray-netcdf-hdf5parallel +module load cray-parallel-netcdf +module load cray-hdf5-parallel +module load w3nco/v2.0.6 +module load nemsio/v2.2.2 +module load bacio/v2.0.2 +module load sp/v2.0.2 +module load sfcio/v1.0.0 +module load sigio/v2.0.1 + + +export FCOMP=ftn +export FFLAGS="-O3 -fp-model precise -g -traceback -r8 -i4 -qopenmp -convert big_endian -assume byterecl" +export WGRIB2API_LIB="/home/larissa.reames/tmp/wgrib2-2/grib2/lib/libwgrib2_api.a" +export WGRIB2API_INC="/home/larissa.reames/tmp/wgrib2-2/grib2/lib" +export WGRIB2_LIB="/home/larissa.reames/tmp/wgrib2-2/grib2/lib/libwgrib2.a" + + +# for debugging +#export FFLAGS="-O0 -g -traceback -r8 -i4 -qopenmp -convert big_endian -check bounds -warn unused -assume byterecl" diff --git a/modulefiles/chgres_cube.wcoss_cray b/modulefiles/chgres_cube.wcoss_cray index 95098edc5..1a9df23b0 100644 --- a/modulefiles/chgres_cube.wcoss_cray +++ b/modulefiles/chgres_cube.wcoss_cray @@ -19,6 +19,9 @@ module load sfcio-intel/1.0.0 # module load esmf/7.1.0r export ESMFMKFILE=/gpfs/hps3/emc/global/noscrub/George.Gayno/esmf/8_0_0_bs20/lib/esmf.mk +export WGRIB2API_INC=/gpfs/hps3/emc/meso/save/Dusan.Jovic/wgrib2/include +export WGRIB2_LIB=/gpfs/hps3/emc/meso/save/Dusan.Jovic/wgrib2/lib/libwgrib2.a + export FCOMP=ftn export FFLAGS="-O3 -fp-model precise -g -r8 -i4 -qopenmp -convert big_endian -assume byterecl" # for debugging diff --git a/modulefiles/chgres_cube.wcoss_dell_p3 b/modulefiles/chgres_cube.wcoss_dell_p3 index 80ccf6589..61297e0cd 100644 --- a/modulefiles/chgres_cube.wcoss_dell_p3 +++ b/modulefiles/chgres_cube.wcoss_dell_p3 @@ -14,6 +14,9 @@ module load bacio/2.0.2 module load sfcio/1.0.0 module load sigio/2.1.0 +export WGRIB2API_INC=/u/Wesley.Ebisuzaki/home/grib2.v2.0.8.intel/lib +export WGRIB2_LIB=/u/Wesley.Ebisuzaki/home/grib2.v2.0.8.intel/lib/libwgrib2.a + export FCOMP=mpif90 export FFLAGS="-O3 -fp-model precise -g -traceback -r8 -i4 -qopenmp -convert big_endian -assume byterecl" # for debugging diff --git a/modulefiles/fv3gfs/fre-nctools.wcoss_cray b/modulefiles/fv3gfs/fre-nctools.wcoss_cray old mode 100755 new mode 100644 diff --git a/modulefiles/fv3gfs/fre-nctools.wcoss_dell_p3 b/modulefiles/fv3gfs/fre-nctools.wcoss_dell_p3 old mode 100755 new mode 100644 diff --git a/modulefiles/modulefile.global_emcsfc_ice_blend.wcoss b/modulefiles/modulefile.global_emcsfc_ice_blend.wcoss old mode 100755 new mode 100644 diff --git a/modulefiles/modulefile.global_emcsfc_ice_blend.wcoss_cray b/modulefiles/modulefile.global_emcsfc_ice_blend.wcoss_cray old mode 100755 new mode 100644 diff --git a/modulefiles/modulefile.global_emcsfc_ice_blend.wcoss_cray_userlib b/modulefiles/modulefile.global_emcsfc_ice_blend.wcoss_cray_userlib old mode 100755 new mode 100644 diff --git a/modulefiles/modulefile.global_emcsfc_ice_blend.wcoss_dell_p3 b/modulefiles/modulefile.global_emcsfc_ice_blend.wcoss_dell_p3 old mode 100755 new mode 100644 diff --git a/modulefiles/modulefile.global_emcsfc_snow2mdl.wcoss b/modulefiles/modulefile.global_emcsfc_snow2mdl.wcoss old mode 100755 new mode 100644 diff --git a/modulefiles/modulefile.global_emcsfc_snow2mdl.wcoss_cray b/modulefiles/modulefile.global_emcsfc_snow2mdl.wcoss_cray old mode 100755 new mode 100644 diff --git a/modulefiles/modulefile.global_emcsfc_snow2mdl.wcoss_cray_userlib b/modulefiles/modulefile.global_emcsfc_snow2mdl.wcoss_cray_userlib old mode 100755 new mode 100644 diff --git a/modulefiles/modulefile.global_emcsfc_snow2mdl.wcoss_dell_p3 b/modulefiles/modulefile.global_emcsfc_snow2mdl.wcoss_dell_p3 old mode 100755 new mode 100644 diff --git a/parm/varmap_tables/GFSphys_var_map.txt b/parm/varmap_tables/GFSphys_var_map.txt new file mode 100644 index 000000000..838110174 --- /dev/null +++ b/parm/varmap_tables/GFSphys_var_map.txt @@ -0,0 +1,23 @@ +dzdt dzdt set_to_fill 0 D +sphum sphum set_to_fill 1E-7 T +liq_wat liq_wat set_to_fill 0 T +o3mr o3mr set_to_fill 1E-7 T +ice_wat ice_wat set_to_fill 0 T +rainwat rainwat set_to_fill 0 T +snowwat snowwat set_to_fill 0 T +graupel graupel set_to_fill 0 T +vtype vtype skip 0 S +sotype stype skip 0 S +vfrac vfrac skip 0 S +fricv uustar set_to_fill 0 S +sfcr zorl set_to_fill 0.01 S +tprcp tprcp set_to_fill 0.00 S +ffmm ffmm set_to_fill 0.00 S +f10m f10m set_to_fill 0.00 S +soilw smc stop 0 S +soill slc set_to_fill 0.0 S +soilt stc stop 0 S +cnwat cnwat set_to_fill 0.0 S +hice icetk set_to_fill 1.5 S +weasd weasd set_to_fill 0.0 S +snod snod set_to_fill 0.0 S diff --git a/parm/varmap_tables/README b/parm/varmap_tables/README new file mode 100644 index 000000000..360341777 --- /dev/null +++ b/parm/varmap_tables/README @@ -0,0 +1,24 @@ +####################################################### +# Description of varmap_tables # +####################################################### + +These files, each named for the phys_suite variable set in the chgres_cube namelist, +control how chgres_cube, when processing grib2 files, handles variables that might +be missing from the grib2 files. Since there are so many different version of grib2 +files, it's often uncertain what fields are available even if you know what source +model the data is coming from. Each file contains : + +Line 1: number of entries in the table +Column 1: Name the code searches for in the table. Do not change. +Column 2: Name the code will use to save the variable in the output file. Unimplemented. +Comumn 3: Behavior when the code can't find the variable in the input file. Options are: + "skip": Don't write to output file. + "set_to_fill": Set to user-specified field value (see column 4). + "stop": Force an exception and stop code execution. Use this if you absolutely + require a field to be present. +Column 4: If column 3 = "set_to_fill", then this value is used to fill in all points + in the input field. These values may be over-written by the code before + output depending on the variable (esp. for surface variables). Be careful + with these values for surface variables. If you set this value too low + (e.g., -100000), the code may run extremely slowly due to variable replacment + at "missing" points. diff --git a/sorc/build_chgres_cube.sh b/sorc/build_chgres_cube.sh new file mode 100755 index 000000000..c13704e1c --- /dev/null +++ b/sorc/build_chgres_cube.sh @@ -0,0 +1,41 @@ +#! /usr/bin/env bash +set -eux + +target=${target:-"NULL"} + +if [[ $target == "linux.gnu" || $target == "linux.intel" ]]; then + unset -f module +else + source ./machine-setup.sh > /dev/null 2>&1 +fi + +cwd=`pwd` + +USE_PREINST_LIBS=${USE_PREINST_LIBS:-"true"} +if [ $USE_PREINST_LIBS = true ]; then + export MOD_PATH + source ../modulefiles/chgres_cube.$target > /dev/null 2>&1 +else + export MOD_PATH=${cwd}/lib/modulefiles + if [ $target = wcoss_cray ]; then + source ../modulefiles/chgres_cube.${target}_userlib > /dev/null 2>&1 + else + source ../modulefiles/chgres_cube.$target > /dev/null 2>&1 + fi +fi + +# Check final exec folder exists +if [ ! -d "../exec" ]; then + mkdir ../exec +fi + +# +# --- Chgres part +# +cd chgres_cube.fd + +make clean +make +make install + +exit diff --git a/sorc/chgres_cube.fd/CMakeLists.txt b/sorc/chgres_cube.fd/CMakeLists.txt index 25338a86e..b9b3ab126 100644 --- a/sorc/chgres_cube.fd/CMakeLists.txt +++ b/sorc/chgres_cube.fd/CMakeLists.txt @@ -32,11 +32,13 @@ target_include_directories(${exe_name} PUBLIC target_include_directories(${exe_name} PRIVATE + ${WGRIB2_INCLUDES} ${NETCDF_INCLUDES} ${NETCDF_INCLUDES_F90} ${MPI_Fortran_INCLUDE_PATH}) target_link_libraries(${exe_name} + ${WGRIB2_LIBRARIES} ${NETCDF_LIBRARIES_F90} ${NETCDF_LIBRARIES} nemsio sfcio_4 ${MPI_Fortran_LIBRARIES} diff --git a/sorc/chgres_cube.fd/atmosphere.F90 b/sorc/chgres_cube.fd/atmosphere.F90 index 98881d95e..e3cebba4e 100644 --- a/sorc/chgres_cube.fd/atmosphere.F90 +++ b/sorc/chgres_cube.fd/atmosphere.F90 @@ -47,7 +47,7 @@ module atmosphere terrain_target_grid use program_setup, only : vcoord_file_target_grid, & - regional, & + regional, & tracers, num_tracers, & atm_weight_file @@ -218,6 +218,7 @@ subroutine atmosphere_driver(localpet) termorderflag=ESMF_TERMORDER_SRCSEQ, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldRegrid", rc) + enddo print*,"- CALL Field_Regrid FOR VERTICAL VELOCITY." @@ -386,7 +387,7 @@ subroutine atmosphere_driver(localpet) !----------------------------------------------------------------------------------- call convert_winds - + !----------------------------------------------------------------------------------- ! Write target data to file. !----------------------------------------------------------------------------------- @@ -492,7 +493,7 @@ subroutine create_atm_esmf_fields allocate(tracers_target_grid(num_tracers)) do n = 1, num_tracers - print*,"- CALL FieldCreate FOR TARGET GRID TRACERS ", trim(tracers(n)) + print*,"- CALL FieldCreate FOR TARGET GRID TRACERS ", trim(tracers(n)) tracers_target_grid(n) = ESMF_FieldCreate(target_grid, & typekind=ESMF_TYPEKIND_R8, & staggerloc=ESMF_STAGGERLOC_CENTER, & @@ -816,7 +817,7 @@ subroutine newpr1(localpet) call error_handler("IN FieldGet", rc) allocate(pi(clb(1):cub(1),clb(2):cub(2),1:levp1_target)) - + if(idvc.eq.2) then do k=1,levp1_target ak = vcoord_target(k,1) @@ -938,7 +939,7 @@ subroutine newps(localpet) farrayPtr=tptr, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGet", rc) - + ! Find specific humidity in the array of tracer fields. do ii = 1, num_tracers @@ -950,7 +951,7 @@ subroutine newps(localpet) farrayPtr=qptr, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGet", rc) - + print*,"- CALL FieldGet FOR SURFACE PRESSURE BEFORE ADJUSTMENT" call ESMF_FieldGet(ps_b4adj_target_grid, & farrayPtr=psptr, rc=rc) @@ -1046,6 +1047,7 @@ subroutine newps(localpet) ! Compute surface pressure over the top. !----------------------------------------------------------------------------------- + if(ls.gt.0) then k=cub(3) gamma=0 @@ -1106,9 +1108,9 @@ subroutine read_vcoord_info print* do k = 1, levp1_target - print*,'VCOORD FOR LEV ', k, 'IS: ', vcoord_target(k,:) + print*,'VCOORD FOR LEV ', k, 'IS: ', vcoord_target(k,:) enddo - + close(14) end subroutine read_vcoord_info @@ -1255,7 +1257,6 @@ SUBROUTINE VINTG CALL TERP3(IM,1,1,1,1,4+NT,(IM*KM1),(IM*KM2), & KM1,IM,IM,Z1,C1,KM2,IM,IM,Z2,C2) - ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! COPY OUTPUT WIND, TEMPERATURE, HUMIDITY AND OTHER TRACERS ! EXCEPT BELOW THE INPUT DOMAIN, LET TEMPERATURE INCREASE WITH A FIXED @@ -1410,7 +1411,7 @@ SUBROUTINE TERP3(IM,IXZ1,IXQ1,IXZ2,IXQ2,NM,NXQ1,NXQ2, & ! REAL(ESMF_KIND_R8) :: J2S ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -! FIND THE SURROUNDING INPUT INTERVAL FOR EACH OUTPUT POINT. +! FIND THE SURROUNDING INPUT INTERVAL FOR EACH OUTPUT POINT. CALL RSEARCH(IM,KM1,IXZ1,KXZ1,Z1,KM2,IXZ2,KXZ2,Z2,1,IM,K1S) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -1422,7 +1423,6 @@ SUBROUTINE TERP3(IM,IXZ1,IXQ1,IXZ2,IXQ2,NM,NXQ1,NXQ2, & !$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(IM,IXZ1,IXQ1,IXZ2), & !$OMP& SHARED(IXQ2,NM,NXQ1,NXQ2,KM1,KXZ1,KXQ1,Z1,Q1,KM2,KXZ2), & !$OMP& SHARED(KXQ2,Z2,Q2,K1S) - DO K2=1,KM2 DO I=1,IM K1=K1S(I,K2) @@ -1490,6 +1490,7 @@ SUBROUTINE TERP3(IM,IXZ1,IXQ1,IXZ2,IXQ2,NM,NXQ1,NXQ2, & ONE/(Z1D-Z1C) ENDIF ENDDO + ! INTERPOLATE. DO N=1,NM DO I=1,IM @@ -1608,23 +1609,24 @@ SUBROUTINE RSEARCH(IM,KM1,IXZ1,KXZ1,Z1,KM2,IXZ2,KXZ2,Z2,IXL2,KXL2,L2) REAL(ESMF_KIND_R8),INTENT(IN) :: Z1(1+(IM-1)*IXZ1+(KM1-1)*KXZ1) REAL(ESMF_KIND_R8),INTENT(IN) :: Z2(1+(IM-1)*IXZ2+(KM2-1)*KXZ2) - INTEGER :: I,K2,L + INTEGER :: I,K2,L REAL(ESMF_KIND_R8) :: Z + ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! FIND THE SURROUNDING INPUT INTERVAL FOR EACH OUTPUT POINT. DO I=1,IM IF (Z1(1+(I-1)*IXZ1).LE.Z1(1+(I-1)*IXZ1+(KM1-1)*KXZ1)) THEN ! INPUT COORDINATE IS MONOTONICALLY ASCENDING. - DO K2=1,KM2 - Z=Z2(1+(I-1)*IXZ2+(K2-1)*KXZ2) + DO K2=1,KM2 + Z=Z2(1+(I-1)*IXZ2+(K2-1)*KXZ2) L=0 DO IF(Z.LT.Z1(1+(I-1)*IXZ1+L*KXZ1)) EXIT L=L+1 IF(L.EQ.KM1) EXIT - ENDDO + ENDDO L2(1+(I-1)*IXL2+(K2-1)*KXL2)=L ENDDO ELSE @@ -1636,7 +1638,7 @@ SUBROUTINE RSEARCH(IM,KM1,IXZ1,KXZ1,Z1,KM2,IXZ2,KXZ2,Z2,IXL2,KXL2,L2) IF(Z.GT.Z1(1+(I-1)*IXZ1+L*KXZ1)) EXIT L=L+1 IF(L.EQ.KM1) EXIT - ENDDO + ENDDO L2(1+(I-1)*IXL2+(K2-1)*KXL2)=L ENDDO ENDIF @@ -1730,7 +1732,7 @@ subroutine compute_zh deallocate(pe0, pn0) end subroutine compute_zh - + subroutine cleanup_target_atm_b4adj_data implicit none diff --git a/sorc/chgres_cube.fd/chgres.F90 b/sorc/chgres_cube.fd/chgres.F90 index 8dedc1122..0d3fcecbd 100644 --- a/sorc/chgres_cube.fd/chgres.F90 +++ b/sorc/chgres_cube.fd/chgres.F90 @@ -14,8 +14,9 @@ program chgres use atmosphere, only : atmosphere_driver use program_setup, only : read_setup_namelist, & + read_varmap, & convert_atm, & - convert_sfc + convert_sfc use model_grid, only : define_target_grid, & define_input_grid, & @@ -60,14 +61,20 @@ program chgres !------------------------------------------------------------------------- call read_setup_namelist + +!------------------------------------------------------------------------- +! Read variable mapping file (used for grib2 input data only). +!------------------------------------------------------------------------- + + call read_varmap !------------------------------------------------------------------------- ! Create esmf grid objects for input and target grids. !------------------------------------------------------------------------- - call define_input_grid(localpet, npets) - call define_target_grid(localpet, npets) + + call define_input_grid(localpet, npets) !------------------------------------------------------------------------- ! Convert atmospheric fields diff --git a/sorc/chgres_cube.fd/grib2_util.F90 b/sorc/chgres_cube.fd/grib2_util.F90 new file mode 100644 index 000000000..64a1cdb71 --- /dev/null +++ b/sorc/chgres_cube.fd/grib2_util.F90 @@ -0,0 +1,77 @@ +module grib2_util + +!-------------------------------------------------------------------------- +! Module: grib2_util +! +! Abstract: Utilities for use when reading grib2 data. +! +!-------------------------------------------------------------------------- + +use esmf + +use model_grid, only : i_input, j_input + +implicit none + +contains + + subroutine rh2spfh(rh_sphum,p,t) + + implicit none + real,parameter :: alpha=-9.477E-4 , & !K^-1, + Tnot=273.15, & !K + Lnot=2.5008E6, & !JKg^-1 + Rv=461.51, & !JKg^-1K^-1 + esnot=611.21 !Pa + + real(esmf_kind_r4), intent(inout), dimension(i_input,j_input) ::rh_sphum + real(esmf_kind_r8), intent(in) :: p, t(i_input,j_input) + + real, dimension(i_input,j_input) :: es, e, rh + + print*,"- CONVERT RH TO SPFH AT LEVEL ", p + + rh = rh_sphum + !print *, 'T = ', T, ' RH = ', RH, ' P = ', P + es = esnot * exp( Lnot/Rv * ((t-Tnot)/(t*tnot) + alpha * LOG(t/Tnot) - alpha * (t-Tnot)/ t)) + !print *, 'es = ', es + e = rh * es / 100.0 + !print *, 'e = ', e + rh_sphum = 0.622 * e / p + !print *, 'q = ', sphum + + !if (P .eq. 100000.0) THEN + ! print *, 'T = ', T, ' RH = ', RH, ' P = ', P, ' es = ', es, ' e = ', e, ' q = ', sphum + !end if + +end subroutine RH2SPFH + +subroutine convert_omega(omega,p,t,q,clb,cub) + + implicit none + real(esmf_kind_r8), pointer :: omega(:,:,:), p(:,:,:), t(:,:,:), q(:,:,:),omtmp,ptmp + + integer :: clb(3), cub(3), i ,j, k + + real, parameter :: Rd = 287.15_esmf_kind_r8, & !JKg^-1K^-1 + Rv=461.51_esmf_kind_r8, & !JKg^-1K^-1 + g = 9.81_esmf_kind_r8 ! ms^-2 + + real(esmf_kind_r8) :: tv, w + + do k = clb(3),cub(3) + do j = clb(2),cub(2) + do i = clb(1),cub(1) + tv = t(i,j,k)*(1+Rd/Rv*q(i,j,k)) + omtmp=>omega(i,j,k) + ptmp=>p(i,j,k) + + w = -1 * omtmp * Rd * tv / (ptmp * g) + omega(i,j,k)=w + enddo + enddo + enddo + +end subroutine convert_omega + + end module grib2_util diff --git a/sorc/chgres_cube.fd/input_data.F90 b/sorc/chgres_cube.fd/input_data.F90 index c813e0a6a..5795fd6f2 100644 --- a/sorc/chgres_cube.fd/input_data.F90 +++ b/sorc/chgres_cube.fd/input_data.F90 @@ -31,20 +31,23 @@ module input_data nst_files_input_grid, & sfc_files_input_grid, & atm_files_input_grid, & + grib2_file_input_grid, & atm_core_files_input_grid, & atm_tracer_files_input_grid, & convert_nst, & orog_dir_input_grid, & orog_files_input_grid, & tracers_input, num_tracers, & - input_type + input_type, tracers, & + get_var_cond, read_from_input use model_grid, only : input_grid, & i_input, j_input, & ip1_input, jp1_input, & num_tiles_input_grid, & latitude_input_grid, & - longitude_input_grid + longitude_input_grid, & + inv_file implicit none @@ -96,8 +99,10 @@ module input_data type(esmf_field), public :: veg_type_input_grid ! vegetation type type(esmf_field), public :: z0_input_grid ! roughness length - integer, parameter, public :: lsoil_input=4 ! # of soil layers, + integer, public :: lsoil_input=4 ! # of soil layers, ! # hardwire for now + + character(len=50), private, allocatable :: slevs(:) ! Fields associated with the nst model. @@ -150,6 +155,8 @@ subroutine read_input_atm_data(localpet) call read_input_atm_gfs_gaussian_file(localpet) elseif (trim(input_type) == "gfs_spectral") then ! spectral gfs sigio format. call read_input_atm_gfs_spectral_file(localpet) + elseif (trim(input_type) == "grib2") then + call read_input_atm_grib2_file(localpet) endif end subroutine read_input_atm_data @@ -484,6 +491,8 @@ subroutine read_input_sfc_data(localpet) call read_input_sfc_gfs_gaussian_file(localpet) elseif (trim(input_type) == "gfs_spectral") then call read_input_sfc_gfs_sfcio_file(localpet) + elseif (trim(input_type) == "grib2") then + call read_input_sfc_grib2_file(localpet) endif end subroutine read_input_sfc_data @@ -2149,6 +2158,567 @@ subroutine read_input_atm_history_file(localpet) call ESMF_FieldDestroy(dpres_input_grid, rc=rc) end subroutine read_input_atm_history_file + +!--------------------------------------------------------------------------- +! Read input grid atmospheric grib2 files. +!--------------------------------------------------------------------------- + + subroutine read_input_atm_grib2_file(localpet) + + use wgrib2api + + use grib2_util, only : rh2spfh, convert_omega + + implicit none + + integer, intent(in) :: localpet + + integer, parameter :: ntrac_max=14 + + character(len=300) :: the_file + character(len=20) :: vlevtyp, vname, lvl_str,lvl_str_space, & + trac_names_grib_1(ntrac_max), & + trac_names_grib_2(ntrac_max), & + trac_names_vmap(ntrac_max), & + tracers_input_grib_1(num_tracers), & + tracers_input_grib_2(num_tracers), & + tmpstr, & + method, tracers_input_vmap(num_tracers), & + tracers_default(ntrac_max), vname2 + character (len=500) :: metadata + + integer :: i, j, k, n, lvl_str_space_len + integer :: rc, clb(3), cub(3) + integer :: vlev, iret,varnum + + integer :: len_str + logical :: lret + + logical :: conv_omega=.false., & + hasspfh=.true. + + real(esmf_kind_r8), allocatable :: rlevs(:) + real(esmf_kind_r4), allocatable :: dummy2d(:,:) + real(esmf_kind_r8), allocatable :: dummy3d(:,:,:), dummy2d_8(:,:) + real(esmf_kind_r8), pointer :: presptr(:,:,:), psptr(:,:),tptr(:,:,:), & + qptr(:,:,:), wptr(:,:,:), & + uptr(:,:,:), vptr(:,:,:) + real(esmf_kind_r4) :: value + real(esmf_kind_r8), parameter :: p0 = 100000.0 + + + tracers(:) = "NULL" + !trac_names_grib = (/":SPFH:",":CLWR:", "O3MR",":CICE:", ":RWMR:",":SNMR:",":GRLE:", & + ! ":TCDC:", ":NCCICE:",":SPNCR:", ":NCONCD:",":PMTF:",":PMTC:",":TKE:"/) + trac_names_grib_1 = (/":var0_2", ":var0_2", ":var0_2", ":var0_2", ":var0_2",":var0_2", \ + ":var0_2", ":var0_2", ":var0_2", ":var0_2", ":var0_2",":var0_2", \ + ":var0_2", ":var0_2"/) + trac_names_grib_2 = (/"_1_0: ", "_1_22: ", "_14_192:", "_1_23: ", "_1_24: ","_1_25: ", \ + "_1_32: ", "_6_1: ", "_6_29: ", "_1_100: ", "_6_28: ","_13_193:", \ + "_13_192:", "_2_2: "/) + trac_names_vmap = (/"sphum ", "liq_wat ", "o3mr ", "ice_wat ", & + "rainwat ", "snowwat ", "graupel ", "cld_amt ", "ice_nc ", & + "rain_nc ", "water_nc", "liq_aero", "ice_aero", & + "sgs_tke "/) + tracers_default = (/"sphum ", "liq_wat ", "o3mr ", "ice_wat ", & + "rainwat ", "snowwat ", "graupel ", "cld_amt ", "ice_nc ", & + "rain_nc ", "water_nc", "liq_aero", "ice_aero", & + "sgs_tke "/) + + the_file = trim(data_dir_input_grid) // "/" // trim(grib2_file_input_grid) + + print*,"- READ ATMOS DATA FROM GRIB2 FILE: ", trim(the_file) + print*,"- USE INVENTORY FILE ", inv_file + + print*,"- OPEN FILE." + inquire(file=the_file,exist=lret) + if (.not.lret) call error_handler("OPENING GRIB2 ATM FILE.", iret) + + print*,"- READ VERTICAL COORDINATE." + iret = grb2_inq(the_file,inv_file,":var_0_2","_0_0:"," hybrid level:") + + if (iret <= 0) then + lvl_str = "mb:" + lvl_str_space = " mb:" + lvl_str_space_len = 4 + iret = grb2_inq(the_file,inv_file,":UGRD:",lvl_str_space) + lev_input=iret + if (localpet == 0) print*,"- DATA IS ON ", lev_input, " ISOBARIC LEVELS." + else + call error_handler("HYBRID VERTICAL COORD DATA NOT SUPPORTED", -1) + endif + + allocate(slevs(lev_input)) + allocate(rlevs(lev_input)) + levp1_input = lev_input + 1 + +! Get the vertical levels, and search string by sequential reads + + do i = 1,lev_input + iret=grb2_inq(the_file,inv_file,':UGRD:',trim(lvl_str),sequential=i-1,desc=metadata) + if (iret.ne.1) call error_handler(" IN SEQUENTIAL FILE READ.", iret) + + j = index(metadata,':UGRD:') + len(':UGRD:') + k = index(metadata,trim(lvl_str_space)) + len(trim(lvl_str_space))-1 + + read(metadata(j:k),*) rlevs(i) + + slevs(i) = metadata(j-1:k) + rlevs(i) = rlevs(i) * 100.0 + if (localpet==0) print*, "- LEVEL = ", slevs(i) + enddo + +! Jili Dong add sort to re-order isobaric levels. + + call quicksort(rlevs,1,lev_input) + + do i = 1,lev_input + write(slevs(i),"(F20.10)") rlevs(i)/100.0 + len_str = len_trim(slevs(i)) + + do while (slevs(i)(len_str:len_str) .eq. '0') + slevs(i) = slevs(i)(:len_str-1) + len_str = len_str - 1 + end do + + if (slevs(i)(len_str:len_str) .eq. '.') then + slevs(i) = slevs(i)(:len_str-1) + len_str = len_str - 1 + end if + + slevs(i) = trim(slevs(i)) + + slevs(i) = ":"//trim(adjustl(slevs(i)))//" mb:" + if (localpet==0) print*, "- LEVEL AFTER SORT = ",slevs(i) + enddo + + if (localpet == 0) print*,"- FIND SPFH OR RH IN FILE" + iret = grb2_inq(the_file,inv_file,trac_names_grib_1(1),trac_names_grib_2(1),lvl_str_space) + + if (iret <= 0) then + iret = grb2_inq(the_file,inv_file, ':var0_2','_1_1:',lvl_str_space) + if (iret <= 0) call error_handler("READING ATMOSPHERIC WATER VAPOR VARIABLE.", iret) + hasspfh = .false. + trac_names_grib_2(1)='_1_1:' + if (localpet == 0) print*,"- FILE CONTAINS RH." + else + if (localpet == 0) print*,"- FILE CONTAINS SPFH." + endif + + print*,"- COUNT NUMBER OF TRACERS TO BE READ IN BASED ON PHYSICS SUITE TABLE" + do n = 1, num_tracers + + vname = tracers_input(n) + + i = maxloc(merge(1.,0.,trac_names_vmap == vname),dim=1) + + tracers_input_grib_1(n) = trac_names_grib_1(i) + tracers_input_grib_2(n) = trac_names_grib_2(i) + tracers_input_vmap(n)=trac_names_vmap(i) + tracers(n)=tracers_default(i) + + enddo + + if (localpet==0) print*, "- NUMBER OF TRACERS IN FILE = ", num_tracers + + if (localpet == 0) print*,"- CALL FieldCreate FOR INPUT GRID SURFACE PRESSURE." + ps_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + if (localpet == 0) print*,"- CALL FieldCreate FOR INPUT GRID TEMPERATURE." + temp_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/lev_input/), rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + allocate(tracers_input_grid(num_tracers)) + + do i = 1,num_tracers + if (localpet == 0) print*,"- CALL FieldCreate FOR INPUT GRID TRACER ", trim(tracers_input(i)) + + tracers_input_grid(i) = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/lev_input/), rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + enddo + + if (localpet == 0) print*,"- CALL FieldCreate FOR INPUT GRID U." + u_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/lev_input/), rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + if (localpet == 0) print*,"- CALL FieldCreate FOR INPUT GRID V." + v_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/lev_input/), rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + if (localpet == 0) print*,"- CALL FieldCreate FOR INPUT GRID DZDT." + dzdt_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/lev_input/), rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + if (localpet == 0) print*,"- CALL FieldCreate FOR INPUT GRID TERRAIN." + terrain_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + if (localpet == 0) then + allocate(dummy2d(i_input,j_input)) + allocate(dummy2d_8(i_input,j_input)) + allocate(dummy3d(i_input,j_input,lev_input)) + else + allocate(dummy2d(0,0)) + allocate(dummy2d_8(0,0)) + allocate(dummy3d(0,0,0)) + endif + +!----------------------------------------------------------------------- +! Fields in non-native files read in from top to bottom. We will +! flip indices later. This program expects bottom to top. +!----------------------------------------------------------------------- + + if (localpet == 0) then + print*,"- READ TEMPERATURE." + vname = ":TMP:" + do vlev = 1, lev_input + iret = grb2_inq(the_file,inv_file,vname,slevs(vlev),data2=dummy2d) + if (iret<=0) then + call error_handler("READING IN TEMPERATURE AT LEVEL "//trim(slevs(vlev)),iret) + endif + dummy3d(:,:,vlev) = real(dummy2d,esmf_kind_r8) + print*,'temp check after read ',vlev, dummy3d(1,1,vlev) + enddo + endif + + if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT GRID TEMPERATURE." + call ESMF_FieldScatter(temp_input_grid, dummy3d, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + do n = 1, num_tracers + + if (localpet == 0) print*,"- READ ", trim(tracers_input_vmap(n)) + vname = tracers_input_vmap(n) + call get_var_cond(vname,this_miss_var_method=method, this_miss_var_value=value, & + this_field_var_name=tmpstr,loc=varnum) + if (n==1 .and. .not. hasspfh) then + print*,"- CALL FieldGather TEMPERATURE." + call ESMF_FieldGather(temp_input_grid,dummy3d,rootPet=0, tile=1, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + endif + if (localpet == 0) then + vname = trim(tracers_input_grib_1(n)) + vname2 = trim(tracers_input_grib_2(n)) + + do vlev = 1, lev_input + iret = grb2_inq(the_file,inv_file,vname,slevs(vlev),vname2,data2=dummy2d) + + if (iret <= 0) then + call handle_grib_error(vname, slevs(vlev),method,value,varnum,iret,var=dummy2d) + if (iret==1) then ! missing_var_method == skip or no entry + if (trim(vname2)=="_1_0:" .or. trim(vname2) == "_1_1:" .or. & + trim(vname2) == ":14:192:") then + call error_handler("READING IN "//trim(vname)//" AT LEVEL "//trim(slevs(vlev))& + //". SET A FILL VALUE IN THE VARMAP TABLE IF THIS ERROR IS NOT DESIRABLE.",iret) + else + exit + endif + endif + endif + + if (n==1 .and. .not. hasspfh) then + call rh2spfh(dummy2d,rlevs(vlev),dummy3d(:,:,vlev)) + endif + + print*,'tracer ',vlev, maxval(dummy2d),minval(dummy2d) + dummy3d(:,:,vlev) = real(dummy2d,esmf_kind_r8) + enddo + endif + + if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT ", trim(tracers_input_vmap(n)) + call ESMF_FieldScatter(tracers_input_grid(n), dummy3d, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + enddo + + if (localpet==0) then + do vlev = 1, lev_input + + vname = ":var0_2" + vname2 = "_2_2:" + iret = grb2_inq(the_file,inv_file,vname,vname2,slevs(vlev),data2=dummy2d) + if (iret<=0) then + call error_handler("READING UWIND AT LEVEL "//trim(slevs(vlev)),iret) + endif + + print*, 'max, min U ', minval(dummy2d), maxval(dummy2d) + dummy3d(:,:,vlev) = real(dummy2d,esmf_kind_r8) + + enddo + endif + + if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT U-WIND." + call ESMF_FieldScatter(u_input_grid, dummy3d, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + if (localpet==0) then + do vlev = 1, lev_input + + vname = ":var0_2" + vname2 = "_2_3:" + iret = grb2_inq(the_file,inv_file,vname,vname2,slevs(vlev),data2=dummy2d) + if (iret<=0) then + call error_handler("READING VWIND AT LEVEL "//trim(slevs(vlev)),iret) + endif + + print*, 'max, min V ', minval(dummy2d), maxval(dummy2d) + dummy3d(:,:,vlev) = real(dummy2d,esmf_kind_r8) + + enddo + endif + + if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT V-WIND." + call ESMF_FieldScatter(v_input_grid, dummy3d, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) then + print*,"- READ SURFACE PRESSURE." + !vname = ":PRES:" + vname = ":var0_2" + vname2 = "_3_0:" + vlevtyp = ":surface:" + iret = grb2_inq(the_file,inv_file,vname,vname2,vlevtyp,data2=dummy2d) + if (iret <= 0) call error_handler("READING SURFACE PRESSURE RECORD.", iret) + dummy2d_8 = real(dummy2d,esmf_kind_r8) + endif + + if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT GRID SURFACE PRESSURE." + call ESMF_FieldScatter(ps_input_grid, dummy2d_8, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) then + print*,"- READ DZDT." + vname = "dzdt" + call get_var_cond(vname,this_miss_var_method=method, this_miss_var_value=value, & + loc=varnum) + !vname = ":DZDT:" + vname = ":var0_2" + vname2 = "_2_9:" + do vlev = 1, lev_input + iret = grb2_inq(the_file,inv_file,vname,vname2,slevs(vlev),data2=dummy2d) + if (iret <= 0 ) then + print*,"DZDT not available at level ", trim(slevs(vlev)), " so checking for VVEL" + !vname = ":VVEL:" + vname2 = "_2_8:" + iret = grb2_inq(the_file,inv_file,vname,vname2,slevs(vlev),data2=dummy2d) + + + if (iret <= 0) then + call handle_grib_error(vname, slevs(vlev),method,value,varnum,iret,var=dummy2d) + if (iret==1) then ! missing_var_method == skip + cycle + endif + else + conv_omega = .true. + endif + + endif + print*,'dzdt ',vlev, maxval(dummy2d),minval(dummy2d) + dummy3d(:,:,vlev) = dummy2d + enddo + endif + + if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT DZDT." + call ESMF_FieldScatter(dzdt_input_grid, dummy3d, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) then + print*,"- READ TERRAIN." + !vname = ":HGT:" + vname = ":var0_2" + vname2 = "_3_5:" + vlevtyp = ":surface:" + iret = grb2_inq(the_file,inv_file,vname,vname2,vlevtyp,data2=dummy2d) + if (iret <= 0) call error_handler("READING TERRAIN HEIGHT RECORD.", iret) + dummy2d_8 = real(dummy2d,esmf_kind_r8) + endif + + if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT GRID TERRAIN." + call ESMF_FieldScatter(terrain_input_grid, dummy2d_8, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) print*,"- CALL FieldCreate FOR INPUT GRID PRESSURE." + pres_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/lev_input/), rc=rc) + + deallocate(dummy2d, dummy3d, dummy2d_8) + +!--------------------------------------------------------------------------- +! Flip 'z' indices to all 3-d variables. Data is read in from model +! top to surface. This program expects surface to model top. +!--------------------------------------------------------------------------- + + if (localpet == 0) print*,"- CALL FieldGet FOR SURFACE PRESSURE." + nullify(psptr) + call ESMF_FieldGet(ps_input_grid, & + farrayPtr=psptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + nullify(presptr) + if (localpet == 0) print*,"- CALL FieldGet FOR 3-D PRESSURE." + call ESMF_FieldGet(pres_input_grid, & + computationalLBound=clb, & + computationalUBound=cub, & + farrayPtr=presptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + nullify(tptr) + if (localpet == 0) print*,"- CALL FieldGet TEMPERATURE." + call ESMF_FieldGet(temp_input_grid, & + farrayPtr=tptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + nullify(uptr) + if (localpet == 0) print*,"- CALL FieldGet FOR U" + call ESMF_FieldGet(u_input_grid, & + farrayPtr=uptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + nullify(vptr) + if (localpet == 0) print*,"- CALL FieldGet FOR V" + call ESMF_FieldGet(v_input_grid, & + farrayPtr=vptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + nullify(wptr) + if (localpet == 0) print*,"- CALL FieldGet FOR W" + call ESMF_FieldGet(dzdt_input_grid, & + farrayPtr=wptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + if (localpet == 0) print*,"- CALL FieldGet FOR TRACERS." + do n=1,num_tracers + nullify(qptr) + call ESMF_FieldGet(tracers_input_grid(n), & + farrayPtr=qptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + do i = clb(1),cub(1) + do j = clb(2),cub(2) + qptr(i,j,:) = qptr(i,j,lev_input:1:-1) + end do + end do + end do + + do i = clb(1),cub(1) + do j = clb(2),cub(2) + presptr(i,j,:) = rlevs(lev_input:1:-1) + tptr(i,j,:) = tptr(i,j,lev_input:1:-1) + uptr(i,j,:) = uptr(i,j,lev_input:1:-1) + vptr(i,j,:) = vptr(i,j,lev_input:1:-1) + wptr(i,j,:) = wptr(i,j,lev_input:1:-1) + end do + end do + + if (localpet == 0) then + print*,'psfc is ',clb(1),clb(2),psptr(clb(1),clb(2)) + print*,'pres is ',cub(1),cub(2),presptr(cub(1),cub(2),:) + + print*,'pres check 1',localpet,maxval(presptr(clb(1):cub(1),clb(2):cub(2),1)), & + minval(presptr(clb(1):cub(1),clb(2):cub(2),1)) + print*,'pres check lev',localpet,maxval(presptr(clb(1):cub(1),clb(2):cub(2), & + lev_input)),minval(presptr(clb(1):cub(1),clb(2):cub(2),lev_input)) + endif + +!--------------------------------------------------------------------------- +! Convert from 2-d to 3-d component winds. +!--------------------------------------------------------------------------- + + call convert_winds + +!--------------------------------------------------------------------------- +! Convert dpdt to dzdt if needed +!--------------------------------------------------------------------------- + + if (conv_omega) then + + if (localpet == 0) print*,"- CONVERT FROM OMEGA TO DZDT." + + nullify(tptr) + if (localpet == 0) print*,"- CALL FieldGet TEMPERATURE." + call ESMF_FieldGet(temp_input_grid, & + farrayPtr=tptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + nullify(qptr) + if (localpet == 0) print*,"- CALL FieldGet SPECIFIC HUMIDITY." + call ESMF_FieldGet(tracers_input_grid(1), & + computationalLBound=clb, & + computationalUBound=cub, & + farrayPtr=qptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + nullify(wptr) + if (localpet == 0) print*,"- CALL FieldGet DZDT." + call ESMF_FieldGet(dzdt_input_grid, & + computationalLBound=clb, & + computationalUBound=cub, & + farrayPtr=wptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + nullify(presptr) + call ESMF_FieldGet(pres_input_grid, & + farrayPtr=presptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + call convert_omega(wptr,presptr,tptr,qptr,clb,cub) + + endif + + end subroutine read_input_atm_grib2_file !--------------------------------------------------------------------------- ! Read input grid surface data from a spectral gfs gaussian sfcio file. @@ -3738,6 +4308,516 @@ subroutine read_input_sfc_history_file(localpet) end subroutine read_input_sfc_history_file + subroutine read_input_sfc_grib2_file(localpet) + + use wgrib2api + + implicit none + + integer, intent(in) :: localpet + + character(len=250) :: the_file + character(len=20) :: vname, vname_file,slev + + character(len=50) :: method + + integer :: rc, varnum, iret, i, j,k + integer, parameter :: icet_default = 265.0 + + logical :: exist + + real(esmf_kind_r4) :: value + + real(esmf_kind_r4), allocatable :: dummy2d(:,:),tsk_save(:,:),icec_save(:,:) + real(esmf_kind_r8), allocatable :: dummy2d_8(:,:) + real(esmf_kind_r8), allocatable :: dummy3d(:,:,:) + integer(esmf_kind_i4), allocatable :: slmsk_save(:,:) + + + the_file = trim(data_dir_input_grid) // "/" // trim(grib2_file_input_grid) + + print*,"- READ SFC DATA FROM GRIB2 FILE: ", trim(the_file) + inquire(file=the_file,exist=exist) + if (.not.exist) then + iret = 1 + call error_handler("OPENING GRIB2 FILE.", iret) + end if + + lsoil_input = grb2_inq(the_file, inv_file, ':TSOIL:',' below ground:') + print*, "- FILE HAS ", lsoil_input, " SOIL LEVELS" + if (lsoil_input <= 0) call error_handler("COUNTING SOIL LEVELS.", rc) + + if (localpet == 0) then + allocate(dummy2d(i_input,j_input)) + allocate(slmsk_save(i_input,j_input)) + allocate(tsk_save(i_input,j_input)) + allocate(icec_save(i_input,j_input)) + allocate(dummy2d_8(i_input,j_input)) + allocate(dummy3d(i_input,j_input,lsoil_input)) + else + allocate(dummy3d(0,0,0)) + allocate(dummy2d_8(0,0)) + allocate(dummy2d(0,0)) + + endif + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! These variables are always in grib files, or are required, so no need to check for them + ! in the varmap table. If they can't be found in the input file, then stop the program. + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + if (localpet == 0) then + print*,"- READ TERRAIN." + rc = grb2_inq(the_file, inv_file, ':HGT:',':surface:', data2=dummy2d) + if (rc /= 1) call error_handler("READING TERRAIN.", rc) + print*,'orog ',maxval(dummy2d),minval(dummy2d) + endif + + print*,"- CALL FieldScatter FOR INPUT TERRAIN." + call ESMF_FieldScatter(terrain_input_grid, real(dummy2d,esmf_kind_r8),rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + +if (localpet == 0) then + print*,"- READ SEAICE FRACTION." + rc = grb2_inq(the_file, inv_file, ':ICEC:',':surface:', data2=dummy2d) + if (rc /= 1) call error_handler("READING SEAICE FRACTION.", rc) + !dummy2d = dummy2d(i_input:1:-1,j_input:1:-1) + print*,'icec ',maxval(dummy2d),minval(dummy2d) + icec_save = dummy2d + endif + + print*,"- CALL FieldScatter FOR INPUT GRID SEAICE FRACTION." + call ESMF_FieldScatter(seaice_fract_input_grid,real(dummy2d,esmf_kind_r8),rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + +!---------------------------------------------------------------------------------- +! GFS v14 and v15.2 grib data has two land masks. LANDN is created by +! nearest neighbor interpolation. LAND is created by bilinear interpolation. +! LANDN matches the bitmap. So use it first. For other GFS versions, use LAND. +! Mask in grib file is '1' (land), '0' (not land). Add sea/lake ice category +! '2' based on ice concentration. +!---------------------------------------------------------------------------------- + + if (localpet == 0) then + print*,"- READ LANDSEA MASK." + rc = grb2_inq(the_file, inv_file, ':LANDN:',':surface:', data2=dummy2d) + + if (rc /= 1) then + rc = grb2_inq(the_file, inv_file, ':LAND:',':surface:', data2=dummy2d) + if (rc /= 1) call error_handler("READING LANDSEA MASK.", rc) + endif + + do j = 1, j_input + do i = 1, i_input + if(dummy2d(i,j) < 0.5_esmf_kind_r4) dummy2d(i,j)=0.0_esmf_kind_r4 + if(icec_save(i,j) > 0.15_esmf_kind_r4) then + !if (dummy2d(i,j) == 0.0_esmf_kind_r4) print*, "CONVERTING WATER TO SEA/LAKE ICE AT ", i, j + dummy2d(i,j) = 2.0_esmf_kind_r4 + endif + enddo + enddo + + slmsk_save = nint(dummy2d) + + deallocate(icec_save) + endif + + print*,"- CALL FieldScatter FOR INPUT LANDSEA MASK." + call ESMF_FieldScatter(landsea_mask_input_grid,real(dummy2d,esmf_kind_r8),rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) then + print*,"- READ SEAICE SKIN TEMPERATURE." + rc = grb2_inq(the_file, inv_file, ':TMP:',':surface:', data2=dummy2d) + if (rc /= 1) call error_handler("READING SEAICE SKIN TEMP.", rc) + print*,'ti ',maxval(dummy2d),minval(dummy2d) + endif + + print*,"- CALL FieldScatter FOR INPUT GRID SEAICE SKIN TEMPERATURE." + call ESMF_FieldScatter(seaice_skin_temp_input_grid,real(dummy2d,esmf_kind_r8),rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + +!---------------------------------------------------------------------------------- +! Read snow fields. Zero out at non-land points and undefined points (points +! removed using the bitmap). Program expects depth and liquid equivalent +! in mm. +!---------------------------------------------------------------------------------- + + if (localpet == 0) then + print*,"- READ SNOW LIQUID EQUIVALENT." + rc = grb2_inq(the_file, inv_file, ':WEASD:',':surface:',':anl:',data2=dummy2d) + if (rc /= 1) then + rc = grb2_inq(the_file, inv_file, ':WEASD:',':surface:','hour fcst:',data2=dummy2d) + if (rc /= 1) call error_handler("READING SNOW LIQUID EQUIVALENT.", rc) + endif + do j = 1, j_input + do i = 1, i_input + if(slmsk_save(i,j) == 0) dummy2d(i,j) = 0.0_esmf_kind_r4 + if(dummy2d(i,j) == grb2_UNDEFINED) dummy2d(i,j) = 0.0_esmf_kind_r4 + enddo + enddo +! print*,'weasd ',maxval(dummy2d),minval(dummy2d) + endif + + print*,"- CALL FieldScatter FOR INPUT GRID SNOW LIQUID EQUIVALENT." + call ESMF_FieldScatter(snow_liq_equiv_input_grid,real(dummy2d,esmf_kind_r8),rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) then + print*,"- READ SNOW DEPTH." + rc = grb2_inq(the_file, inv_file, ':SNOD:',':surface:', data2=dummy2d) + if (rc /= 1) call error_handler("READING SNOW DEPTH.", rc) + where(dummy2d == grb2_UNDEFINED) dummy2d = 0.0_esmf_kind_r4 + dummy2d = dummy2d*1000.0 ! Grib2 files have snow depth in (m), fv3 expects it in mm + where(slmsk_save == 0) dummy2d = 0.0_esmf_kind_r4 +! print*,'snod ',maxval(dummy2d),minval(dummy2d) + endif + + print*,"- CALL FieldScatter FOR INPUT GRID SNOW DEPTH." + call ESMF_FieldScatter(snow_depth_input_grid,real(dummy2d,esmf_kind_r8),rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) then + print*,"- READ T2M." + rc = grb2_inq(the_file, inv_file, ':TMP:',':2 m above ground:',data2=dummy2d) + if (rc <= 0) call error_handler("READING T2M.", rc) + + print*,'t2m ',maxval(dummy2d),minval(dummy2d) + endif + + print*,"- CALL FieldScatter FOR INPUT GRID T2M." + call ESMF_FieldScatter(t2m_input_grid,real(dummy2d,esmf_kind_r8), rootpet=0,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) then + print*,"- READ Q2M." + rc = grb2_inq(the_file, inv_file, ':SPFH:',':2 m above ground:',data2=dummy2d) + if (rc <=0) call error_handler("READING Q2M.", rc) + print*,'q2m ',maxval(dummy2d),minval(dummy2d) + endif + + print*,"- CALL FieldScatter FOR INPUT GRID Q2M." + call ESMF_FieldScatter(q2m_input_grid,real(dummy2d,esmf_kind_r8), rootpet=0,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + +if (localpet == 0) then + print*,"- READ SKIN TEMPERATURE." + rc = grb2_inq(the_file, inv_file, ':TMP:',':surface:', data2=dummy2d) + if (rc <= 0 ) call error_handler("READING SKIN TEMPERATURE.", rc) + tsk_save(:,:) = dummy2d + dummy2d_8 = real(dummy2d,esmf_kind_r8) + + print*,'tmp ',maxval(dummy2d),minval(dummy2d) + endif + + print*,"- CALL FieldScatter FOR INPUT GRID SKIN TEMPERATURE" + call ESMF_FieldScatter(skin_temp_input_grid,real(dummy2d,esmf_kind_r8),rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) dummy2d = 0.0 + + print*,"- CALL FieldScatter FOR INPUT GRID SRFLAG" + call ESMF_FieldScatter(srflag_input_grid,real(dummy2d,esmf_kind_r8), rootpet=0,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + +! Soil type is not available. Set to a large negative fill value. + + dummy2d_8 = -99999.0_esmf_kind_r8 + + print*,"- CALL FieldScatter FOR INPUT GRID SOIL TYPE." + call ESMF_FieldScatter(soil_type_input_grid,dummy2d_8, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Begin variables whose presence in grib2 files varies, but no climatological + ! data is + ! available, so we have to account for values in the varmap table + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + if (localpet == 0) then + print*,"- READ SEAICE DEPTH." + vname="hice" + slev=":surface:" + call get_var_cond(vname,this_miss_var_method=method,this_miss_var_value=value, & + loc=varnum) + vname=":ICETK:" + rc= grb2_inq(the_file, inv_file, vname,slev, data2=dummy2d) + if (rc <= 0) then + call handle_grib_error(vname, slev ,method,value,varnum,rc, var= dummy2d) + if (rc==1) then ! missing_var_method == skip or no entry in varmap table + print*, "WARNING: "//trim(vname)//" NOT AVAILABLE IN FILE. THIS FIELD WILL BE"//& + " REPLACED WITH CLIMO. SET A FILL "// & + "VALUE IN THE VARMAP TABLE IF THIS IS NOT DESIRABLE." + dummy2d(:,:) = 0.0_esmf_kind_r4 + endif + endif + dummy2d_8= real(dummy2d,esmf_kind_r8) + print*,'hice ',maxval(dummy2d),minval(dummy2d) + + endif + + print*,"- CALL FieldScatter FOR INPUT GRID SEAICE DEPTH." + call ESMF_FieldScatter(seaice_depth_input_grid,dummy2d_8, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) then + print*,"- READ TPRCP." + vname="tprcp" + slev=":surface:" + call get_var_cond(vname,this_miss_var_method=method,this_miss_var_value=value, & + loc=varnum) + vname=":TPRCP:" + rc= grb2_inq(the_file, inv_file, vname,slev, data2=dummy2d) + if (rc <= 0) then + call handle_grib_error(vname, slev ,method,value,varnum,rc, var= dummy2d) + if (rc==1) then ! missing_var_method == skip or no entry in varmap table + print*, "WARNING: "//trim(vname)//" NOT AVAILABLE IN FILE. THIS FIELD WILL NOT"//& + " BE WRITTEN TO THE INPUT FILE. SET A FILL "// & + "VALUE IN THE VARMAP TABLE IF THIS IS NOT DESIRABLE." + dummy2d(:,:) = 0.0_esmf_kind_r4 + endif + endif + dummy2d_8= real(dummy2d,esmf_kind_r8) + print*,'tprcp ',maxval(dummy2d),minval(dummy2d) + endif + + print*,"- CALL FieldScatter FOR INPUT GRID TPRCP." + call ESMF_FieldScatter(tprcp_input_grid,dummy2d_8, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) then + print*,"- READ FFMM." + vname="ffmm" + slev=":surface:" + call get_var_cond(vname,this_miss_var_method=method,this_miss_var_value=value, & + loc=varnum) + vname=":FFMM:" + rc= grb2_inq(the_file, inv_file, vname,slev, data2=dummy2d) + if (rc <= 0) then + call handle_grib_error(vname, slev ,method,value,varnum,rc, var= dummy2d) + if (rc==1) then ! missing_var_method == skip or no entry in varmap table + print*, "WARNING: "//trim(vname)//" NOT AVAILABLE IN FILE. THIS FIELD WILL NOT"//& + " BE WRITTEN TO THE INPUT FILE. SET A FILL "// & + "VALUE IN THE VARMAP TABLE IF THIS IS NOT DESIRABLE." + dummy2d(:,:) = 0.0_esmf_kind_r4 + endif + endif + dummy2d_8= real(dummy2d,esmf_kind_r8) + print*,'ffmm ',maxval(dummy2d),minval(dummy2d) + endif + + print*,"- CALL FieldScatter FOR INPUT GRID FFMM" + call ESMF_FieldScatter(ffmm_input_grid,dummy2d_8, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) then + print*,"- READ USTAR." + vname="fricv" + slev=":surface:" + call get_var_cond(vname,this_miss_var_method=method,this_miss_var_value=value, & + loc=varnum) + vname=":FRICV:" + rc= grb2_inq(the_file, inv_file, vname,slev, data2=dummy2d) + if (rc <= 0) then + call handle_grib_error(vname, slev ,method,value,varnum,rc, var= dummy2d) + if (rc==1) then ! missing_var_method == skip or no entry in varmap table + print*, "WARNING: "//trim(vname)//" NOT AVAILABLE IN FILE. THIS FIELD WILL "//& + "REPLACED WITH CLIMO. SET A FILL "// & + "VALUE IN THE VARMAP TABLE IF THIS IS NOT DESIRABLE." + dummy2d(:,:) = 0.0_esmf_kind_r4 + endif + endif + dummy2d_8= real(dummy2d,esmf_kind_r8) + print*,'fricv ',maxval(dummy2d),minval(dummy2d) + endif + + print*,"- CALL FieldScatter FOR INPUT GRID USTAR" + call ESMF_FieldScatter(ustar_input_grid,dummy2d_8, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) then + print*,"- READ F10M." + vname="f10m" + slev=":10 m above ground:" + call get_var_cond(vname,this_miss_var_method=method,this_miss_var_value=value, & + loc=varnum) + vname=":F10M:" + rc= grb2_inq(the_file, inv_file, vname,slev, data2=dummy2d) + if (rc <= 0) then + call handle_grib_error(vname, slev ,method,value,varnum,rc, var= dummy2d) + if (rc==1) then ! missing_var_method == skip or no entry in varmap table + print*, "WARNING: "//trim(vname)//" NOT AVAILABLE IN FILE. THIS FIELD WILL NOT"//& + " BE WRITTEN TO THE INPUT FILE. SET A FILL "// & + "VALUE IN THE VARMAP TABLE IF THIS IS NOT DESIRABLE." + dummy2d(:,:) = 0.0_esmf_kind_r4 + endif + endif + dummy2d_8= real(dummy2d,esmf_kind_r8) + print*,'f10m ',maxval(dummy2d),minval(dummy2d) + endif + + print*,"- CALL FieldScatter FOR INPUT GRID F10M." + call ESMF_FieldScatter(f10m_input_grid,dummy2d_8, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) then + print*,"- READ CANOPY MOISTURE CONTENT." + vname="cnwat" + slev=":surface:" + call get_var_cond(vname,this_miss_var_method=method,this_miss_var_value=value, & + loc=varnum) + vname=":CNWAT:" + rc= grb2_inq(the_file, inv_file, vname,slev, data2=dummy2d) + if (rc <= 0) then + call handle_grib_error(vname, slev ,method,value,varnum,rc, var= dummy2d) + if (rc==1) then ! missing_var_method == skip or no entry in varmap table + print*, "WARNING: "//trim(vname)//" NOT AVAILABLE IN FILE. THIS FIELD WILL"//& + " REPLACED WITH CLIMO. SET A FILL "// & + "VALUE IN THE VARMAP TABLE IF THIS IS NOT DESIRABLE." + dummy2d(:,:) = 0.0_esmf_kind_r4 + endif + endif + dummy2d_8= real(dummy2d,esmf_kind_r8) + print*,'cnwat ',maxval(dummy2d),minval(dummy2d) + endif + + print*,"- CALL FieldScatter FOR INPUT GRID CANOPY MOISTURE CONTENT." + call ESMF_FieldScatter(canopy_mc_input_grid,dummy2d_8, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) then + print*,"- READ Z0." + vname="sfcr" + slev=":surface:" + call get_var_cond(vname,this_miss_var_method=method,this_miss_var_value=value, & + loc=varnum) + vname=":SFCR:" + rc= grb2_inq(the_file, inv_file, vname,slev, data2=dummy2d) + if (rc <= 0) then + call handle_grib_error(vname, slev ,method,value,varnum,rc, var= dummy2d) + if (rc==1) then ! missing_var_method == skip or no entry in varmap table + print*, "WARNING: "//trim(vname)//" NOT AVAILABLE IN FILE. THIS FIELD WILL BE"//& + " REPLACED WITH CLIMO. SET A FILL "// & + "VALUE IN THE VARMAP TABLE IF THIS IS NOT DESIRABLE." + dummy2d(:,:) = 0.0_esmf_kind_r4 + endif + else + ! Grib files have z0 (m), but fv3 expects z0(cm) + dummy2d(:,:) = dummy2d(:,:)*10.0 + endif + dummy2d_8= real(dummy2d,esmf_kind_r8) + print*,'sfcr ',maxval(dummy2d),minval(dummy2d) + + endif + + print*,"- CALL FieldScatter FOR INPUT GRID Z0." + call ESMF_FieldScatter(z0_input_grid,dummy2d_8, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + + deallocate(dummy2d) + + if (localpet == 0) then + print*,"- READ LIQUID SOIL MOISTURE." + vname = "soill" + vname_file = ":SOILL:" + call read_grib_soil(the_file,inv_file,vname,vname_file,dummy3d,rc) !!! NEEDTO HANDLE + !!! SOIL LEVELS + print*,'soill ',maxval(dummy3d),minval(dummy3d) + endif + + print*,"- CALL FieldScatter FOR INPUT LIQUID SOIL MOISTURE." + call ESMF_FieldScatter(soilm_liq_input_grid, dummy3d, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) then + print*,"- READ TOTAL SOIL MOISTURE." + vname = "soilw" + !vname_file = "var2_2_1_7_0_192" !Some files don't recognize this as soilw,so use + vname_file = "var2_2_1_" ! the var number instead + call read_grib_soil(the_file,inv_file,vname,vname_file,dummy3d,rc) + print*,'soilm ',maxval(dummy3d),minval(dummy3d) + endif + +!----------------------------------------------------------------------- +! Vegetation type is not available. However, it is needed to identify +! permanent land ice points. At land ice, the total soil moisture +! is a flag value of '1'. Use this flag as a temporary solution. +!----------------------------------------------------------------------- + + if (localpet == 0) then + dummy2d_8(:,:) = 0.0_esmf_kind_r8 + do j = 1, j_input + do i = 1, i_input + if(slmsk_save(i,j) == 1_esmf_kind_i4 .and. dummy3d(i,j,1) > 0.99) & + dummy2d_8(i,j) = real(veg_type_landice_input,esmf_kind_r8) + enddo + enddo + endif + + print*,"- CALL FieldScatter FOR INPUT VEG TYPE." + call ESMF_FieldScatter(veg_type_input_grid, dummy2d_8, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + + print*,"- CALL FieldScatter FOR INPUT TOTAL SOIL MOISTURE." + call ESMF_FieldScatter(soilm_tot_input_grid, dummy3d, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + +!--------------------------------------------------------------------------------- +! At open water (slmsk==0), the soil temperature array is not used and set +! to the filler value of SST. At lake/sea ice points (slmsk=2), the soil +! temperature array holds ice column temperature. That field is not available +! in GFS grib data, so set to a default value. +!--------------------------------------------------------------------------------- + + if (localpet == 0) then + print*,"- READ SOIL TEMPERATURE." + vname = "soilt" + vname_file = ":TSOIL:" + call read_grib_soil(the_file,inv_file,vname,vname_file,dummy3d,rc) + do k=1,lsoil_input + do j = 1, j_input + do i = 1, i_input + if (slmsk_save(i,j) == 0_esmf_kind_i4 ) dummy3d(i,j,k) = tsk_save(i,j) + if (slmsk_save(i,j) == 2_esmf_kind_i4 ) dummy3d(i,j,k) = icet_default + enddo + enddo + enddo + print*,'soilt ',maxval(dummy3d),minval(dummy3d) + + deallocate(tsk_save, slmsk_save) + endif + + print*,"- CALL FieldScatter FOR INPUT SOIL TEMPERATURE." + call ESMF_FieldScatter(soil_temp_input_grid, dummy3d, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + + deallocate(dummy3d) + deallocate(dummy2d_8) + + end subroutine read_input_sfc_grib2_file + !--------------------------------------------------------------------------- ! Read nst data from tiled history or restart files. !--------------------------------------------------------------------------- @@ -4402,6 +5482,114 @@ subroutine convert_winds call ESMF_FieldDestroy(v_input_grid, rc=rc) end subroutine convert_winds + +subroutine handle_grib_error(vname,lev,method,value,varnum, iret,var,var8,var3d) + + use, intrinsic :: ieee_arithmetic + + implicit none + + real(esmf_kind_r4), intent(in) :: value + real(esmf_kind_r4), intent(inout), optional :: var(:,:) + real(esmf_kind_r8), intent(inout), optional :: var8(:,:) + real(esmf_kind_r8), intent(inout), optional :: var3d(:,:,:) + + character(len=20), intent(in) :: vname, lev, method + + integer, intent(in) :: varnum + integer, intent(inout) :: iret + + iret = 0 + if (varnum == 9999) then + print*, "WARNING: ", trim(vname), " NOT FOUND AT LEVEL ", lev, " IN EXTERNAL FILE ", & + "AND NO ENTRY EXISTS IN VARMAP TABLE. VARIABLE WILL NOT BE USED." + iret = 1 + + return + endif + + if (trim(method) == "skip" ) then + print*, "WARNING: SKIPPING ", trim(vname), " IN FILE" + read_from_input(varnum) = .false. + iret = 1 + elseif (trim(method) == "set_to_fill") then + print*, "WARNING: ,", trim(vname), " NOT AVILABLE AT LEVEL ", trim(lev), & + ". SETTING EQUAL TO FILL VALUE OF ", value + if(present(var)) var(:,:) = value + if(present(var8)) var8(:,:) = value + if(present(var3d)) var3d(:,:,:) = value + elseif (trim(method) == "set_to_NaN") then + print*, "WARNING: ,", trim(vname), " NOT AVILABLE AT LEVEL ", trim(lev), & + ". SETTING EQUAL TO NaNs" + if(present(var)) var(:,:) = ieee_value(var,IEEE_QUIET_NAN) + if(present(var8)) var8(:,:) = ieee_value(var8,IEEE_QUIET_NAN) + if(present(var3d)) var3d(:,:,:) = ieee_value(var3d,IEEE_QUIET_NAN) + elseif (trim(method) == "stop") then + call error_handler("READING "//trim(vname)// " at level "//lev//". TO MAKE THIS NON- & + FATAL, CHANGE STOP TO SKIP FOR THIS VARIABLE IN YOUR VARMAP & + FILE.", iret) + else + call error_handler("ERROR USING MISSING_VAR_METHOD. PLEASE SET VALUES IN" // & + " VARMAP TABLE TO ONE OF: set_to_fill, set_to_NaN,"// & + " , skip, or stop.") + endif + +end subroutine handle_grib_error + +subroutine read_grib_soil(the_file,inv_file,vname,vname_file,dummy3d,rc) + + use wgrib2api + implicit none + + character(len=*), intent(in) :: the_file, inv_file + character(len=20), intent(in) :: vname,vname_file + + integer, intent(out) :: rc + + real(esmf_kind_r8), intent(inout) :: dummy3d(:,:,:) + + real(esmf_kind_r4), allocatable :: dummy2d(:,:) + real(esmf_kind_r4) :: value + integer :: varnum,i + character(len=50) :: slevs(lsoil_input) + character(len=50) :: method + + allocate(dummy2d(i_input,j_input)) + + if(lsoil_input == 4) then + slevs = (/character(24)::':0-0.1 m below ground:', ':0.1-0.4 m below ground:', & + ':0.4-1 m below ground:', ':1-2 m below ground:'/) + else + rc = -1 + call error_handler("reading soil levels. File must have 4 soil levels.") + endif + + call get_var_cond(vname,this_miss_var_method=method,this_miss_var_value=value, & + loc=varnum) + do i = 1,lsoil_input + if (vname_file=="var2_2_1_") then + rc = grb2_inq(the_file,inv_file,vname_file,"_0_192:",slevs(i),data2=dummy2d) + else + rc = grb2_inq(the_file,inv_file,vname_file,slevs(i),data2=dummy2d) + endif + if (rc <= 0) then + call handle_grib_error(vname_file, slevs(i),method,value,varnum,rc,var=dummy2d) + if (rc==1 .and. trim(vname) /= "soill") then + ! missing_var_method == skip or no entry in varmap table + call error_handler("READING IN "//trim(vname)//". SET A FILL "// & + "VALUE IN THE VARMAP TABLE IF THIS ERROR IS NOT DESIRABLE.",rc) + elseif (rc==1) then + dummy3d(:,:,:) = 0.0_esmf_kind_r8 + exit + endif + endif + + dummy3d(:,:,i) = real(dummy2d,esmf_kind_r8) + end do + + deallocate(dummy2d) + + end subroutine read_grib_soil subroutine cleanup_input_atm_data @@ -4490,4 +5678,35 @@ subroutine cleanup_input_sfc_data end subroutine cleanup_input_sfc_data +! Jili Dong add sort subroutine +! quicksort.f -*-f90-*- +! Author: t-nissie +! License: GPLv3 +! Gist: https://gist.github.com/t-nissie/479f0f16966925fa29ea +!! +recursive subroutine quicksort(a, first, last) + implicit none + real*8 a(*), x, t + integer first, last + integer i, j + + x = a( (first+last) / 2 ) + i = first + j = last + do + do while (a(i) < x) + i=i+1 + end do + do while (x < a(j)) + j=j-1 + end do + if (i >= j) exit + t = a(i); a(i) = a(j); a(j) = t + i=i+1 + j=j-1 + end do + if (first < i-1) call quicksort(a, first, i-1) + if (j+1 < last) call quicksort(a, j+1, last) +end subroutine quicksort + end module input_data diff --git a/sorc/chgres_cube.fd/makefile b/sorc/chgres_cube.fd/makefile old mode 100755 new mode 100644 index d7aca1d9b..a247c85b7 --- a/sorc/chgres_cube.fd/makefile +++ b/sorc/chgres_cube.fd/makefile @@ -13,17 +13,21 @@ OBJS = chgres.o \ write_data.o \ search_util.o \ static_data.o \ - utils.o + utils.o \ + grib2_util.o $(CMD): $(OBJS) - $(FCOMP) $(FFLAGS) $(ESMF_F90COMPILEPATHS) -o $(CMD) $(OBJS) $(SP_LIBd) $(NEMSIO_LIB) $(BACIO_LIB4) $(W3NCO_LIBd) $(SFCIO_LIB4) $(SIGIO_LIB4) $(ESMF_F90LINKPATHS) $(ESMF_F90ESMFLINKRPATHS) $(ESMF_F90ESMFLINKLIBS) -g -traceback + $(FCOMP) $(FFLAGS) $(ESMF_F90COMPILEPATHS) -o $(CMD) $(OBJS) $(SP_LIBd) $(NEMSIO_LIB) $(BACIO_LIB4) $(W3NCO_LIBd) $(SFCIO_LIB4) $(SIGIO_LIB4) $(ESMF_F90LINKPATHS) $(ESMF_F90ESMFLINKRPATHS) $(ESMF_F90ESMFLINKLIBS) $(WGRIB2_LIB) model_grid.o: program_setup.o model_grid.F90 - $(FCOMP) $(FFLAGS) -I$(NEMSIO_INC) -I$(SFCIO_INC4) -I$(SIGIO_INC4) $(ESMF_F90COMPILEPATHS) -c model_grid.F90 + $(FCOMP) $(FFLAGS) -I$(NEMSIO_INC) -I$(SFCIO_INC4) -I$(SIGIO_INC4) -I$(WGRIB2API_INC) $(ESMF_F90COMPILEPATHS) -c model_grid.F90 utils.o: utils.f90 $(FCOMP) $(FFLAGS) -c $(ESMF_F90COMPILEPATHS) utils.f90 +grib2_util.o: model_grid.o program_setup.o grib2_util.F90 + $(FCOMP) $(FFLAGS) -c $(ESMF_F90COMPILEPATHS) grib2_util.F90 + program_setup.o: program_setup.f90 $(FCOMP) $(FFLAGS) $(ESMF_F90COMPILEPATHS) -c program_setup.f90 @@ -36,8 +40,8 @@ chgres.o: atmosphere.o model_grid.o program_setup.o surface.o chgres.F90 write_data.o: atmosphere.o model_grid.o program_setup.o surface.o static_data.o write_data.F90 $(FCOMP) $(FFLAGS) $(ESMF_F90COMPILEPATHS) -c write_data.F90 -input_data.o: program_setup.o model_grid.o input_data.F90 - $(FCOMP) $(FFLAGS) -I$(NEMSIO_INC) -I$(SIGIO_INC4) -I$(SFCIO_INC4) $(ESMF_F90COMPILEPATHS) -c input_data.F90 +input_data.o: program_setup.o model_grid.o grib2_util.o input_data.F90 + $(FCOMP) $(FFLAGS) -I$(NEMSIO_INC) -I$(SIGIO_INC4) -I$(SFCIO_INC4) -I$(WGRIB2API_INC) $(ESMF_F90COMPILEPATHS) -c input_data.F90 surface.o: search_util.o model_grid.o input_data.o program_setup.o static_data.o surface.F90 $(FCOMP) $(FFLAGS) $(ESMF_F90COMPILEPATHS) -c surface.F90 diff --git a/sorc/chgres_cube.fd/model_grid.F90 b/sorc/chgres_cube.fd/model_grid.F90 index 48b95f7f4..fa945cbc6 100644 --- a/sorc/chgres_cube.fd/model_grid.F90 +++ b/sorc/chgres_cube.fd/model_grid.F90 @@ -63,6 +63,7 @@ module model_grid private character(len=5), allocatable, public :: tiles_target_grid(:) + character(len=10), public :: inv_file = "chgres.inv" integer, parameter, public :: lsoil_target = 4 ! # soil layers integer, public :: i_input, j_input @@ -119,6 +120,8 @@ subroutine define_input_grid(localpet, npets) trim(input_type) == "gfs_gaussian" .or. & trim(input_type) == "gfs_spectral") then call define_input_grid_gaussian(localpet, npets) + elseif (trim(input_type) == "grib2") then + call define_input_grid_gfs_grib2(localpet,npets) else call define_input_grid_mosaic(localpet, npets) endif @@ -234,9 +237,7 @@ subroutine define_input_grid_gaussian(localpet, npets) coordSys=ESMF_COORDSYS_SPH_DEG, & regDecomp=(/1,npets/), & indexflag=ESMF_INDEX_GLOBAL, rc=rc) - -if(ESMF_logFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__,file=__FILE__)) & + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN GridCreate1PeriDim", rc) print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE." @@ -573,6 +574,192 @@ subroutine define_input_grid_mosaic(localpet, npets) end subroutine define_input_grid_mosaic +!-------------------------------------------------------------------------- +! Define grid object for GFS grib2 data. Only works for data on +! global lat/lon or gaussian grids. +!-------------------------------------------------------------------------- + + subroutine define_input_grid_gfs_grib2(localpet, npets) + + use wgrib2api + + use program_setup, only : data_dir_input_grid, & + grib2_file_input_grid + + implicit none + + integer, intent(in) :: localpet, npets + + character(len=250) :: the_file + + integer :: i, j, rc, clb(2), cub(2) + + real(esmf_kind_r8), allocatable :: latitude(:,:) + real(esmf_kind_r8), allocatable :: longitude(:,:) + real(esmf_kind_r4), allocatable :: lat4(:,:), lon4(:,:) + real(esmf_kind_r8), pointer :: lat_src_ptr(:,:) + real(esmf_kind_r8), pointer :: lon_src_ptr(:,:) + real(esmf_kind_r8), pointer :: lat_corner_src_ptr(:,:) + real(esmf_kind_r8), pointer :: lon_corner_src_ptr(:,:) + real(esmf_kind_r8) :: deltalon + + type(esmf_polekind_flag) :: polekindflag(2) + + print*,"- DEFINE INPUT GRID OBJECT FOR INPUT GRIB2 DATA." + + num_tiles_input_grid = 1 + + the_file = trim(data_dir_input_grid) // "/" // grib2_file_input_grid + print*,'- OPEN AND INVENTORY GRIB2 FILE: ',trim(the_file) + rc=grb2_mk_inv(the_file,inv_file) + if (rc /=0) call error_handler("OPENING GRIB2 FILE",rc) + + rc = grb2_inq(the_file,inv_file,':PRES:',':surface:',nx=i_input, ny=j_input, & + lat=lat4, lon=lon4) + if (rc /= 1) call error_handler("READING GRIB2 FILE", rc) + + ip1_input = i_input + 1 + jp1_input = j_input + 1 + + polekindflag(1:2) = ESMF_POLEKIND_MONOPOLE + + print*,"- CALL GridCreate1PeriDim FOR INPUT GRID." + input_grid = ESMF_GridCreate1PeriDim(minIndex=(/1,1/), & + maxIndex=(/i_input,j_input/), & + polekindflag=polekindflag, & + periodicDim=1, & + poleDim=2, & + coordSys=ESMF_COORDSYS_SPH_DEG, & + regDecomp=(/1,npets/), & + indexflag=ESMF_INDEX_GLOBAL, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridCreate1PeriDim", rc) + + print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE." + latitude_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + name="input_grid_latitude", rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE." + longitude_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + name="input_grid_longitude", rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + allocate(longitude(i_input,j_input)) + allocate(latitude(i_input,j_input)) + + do i = 1, i_input + longitude(i,:) = real(lon4(i,:),kind=esmf_kind_r8) + enddo + + do i = 1, j_input + latitude(:,i) = real(lat4(:,i),kind=esmf_kind_r8) + enddo + + deallocate(lat4, lon4) + + deltalon = abs(longitude(2,1)-longitude(1,1)) + if(localpet==0) print*, "deltalon = ", deltalon + + print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE." + call ESMF_FieldScatter(longitude_input_grid, longitude, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE." + call ESMF_FieldScatter(latitude_input_grid, latitude, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + print*,"- CALL GridAddCoord FOR INPUT GRID." + call ESMF_GridAddCoord(input_grid, & + staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridAddCoord", rc) + + print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD." + nullify(lon_src_ptr) + call ESMF_GridGetCoord(input_grid, & + staggerLoc=ESMF_STAGGERLOC_CENTER, & + coordDim=1, & + farrayPtr=lon_src_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord", rc) + + print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD." + nullify(lat_src_ptr) + call ESMF_GridGetCoord(input_grid, & + staggerLoc=ESMF_STAGGERLOC_CENTER, & + coordDim=2, & + computationalLBound=clb, & + computationalUBound=cub, & + farrayPtr=lat_src_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord", rc) + + do j = clb(2), cub(2) + do i = clb(1), cub(1) + lon_src_ptr(i,j) = longitude(i,j) + if (lon_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_src_ptr(i,j) = lon_src_ptr(i,j) - 360.0_esmf_kind_r8 + lat_src_ptr(i,j) = latitude(i,j) + enddo + enddo + + if(localpet==0) print*, "lon first = ", lon_src_ptr(1:10,1) + if(localpet==0) print*, "lat first = ", lat_src_ptr(1,1:10) + + print*,"- CALL GridAddCoord FOR INPUT GRID." + call ESMF_GridAddCoord(input_grid, & + staggerloc=ESMF_STAGGERLOC_CORNER, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridAddCoord", rc) + + print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD." + nullify(lon_corner_src_ptr) + call ESMF_GridGetCoord(input_grid, & + staggerLoc=ESMF_STAGGERLOC_CORNER, & + coordDim=1, & + farrayPtr=lon_corner_src_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord", rc) + + print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD." + nullify(lat_corner_src_ptr) + call ESMF_GridGetCoord(input_grid, & + staggerLoc=ESMF_STAGGERLOC_CORNER, & + coordDim=2, & + computationalLBound=clb, & + computationalUBound=cub, & + farrayPtr=lat_corner_src_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord", rc) + + do j = clb(2), cub(2) + do i = clb(1), cub(1) + lon_corner_src_ptr(i,j) = longitude(i,1) - (0.5_esmf_kind_r8*deltalon) + if (lon_corner_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_corner_src_ptr(i,j) = lon_corner_src_ptr(i,j) - 360.0_esmf_kind_r8 + if (j == 1) then + lat_corner_src_ptr(i,j) = -90.0_esmf_kind_r8 + cycle + endif + if (j == jp1_input) then + lat_corner_src_ptr(i,j) = +90.0_esmf_kind_r8 + cycle + endif + lat_corner_src_ptr(i,j) = 0.5_esmf_kind_r8 * (latitude(i,j-1)+ latitude(i,j)) + enddo + enddo + + deallocate(latitude,longitude) + + end subroutine define_input_grid_gfs_grib2 + subroutine define_target_grid(localpet, npets) use netcdf diff --git a/sorc/chgres_cube.fd/program_setup.f90 b/sorc/chgres_cube.fd/program_setup.f90 index 001be72d4..25ac2f33f 100644 --- a/sorc/chgres_cube.fd/program_setup.f90 +++ b/sorc/chgres_cube.fd/program_setup.f90 @@ -94,7 +94,8 @@ module program_setup implicit none private - + + character(len=500), public :: varmap_file = "NULL" character(len=500), public :: atm_files_input_grid(6) = "NULL" character(len=500), public :: atm_core_files_input_grid(7) = "NULL" character(len=500), public :: atm_tracer_files_input_grid(6) = "NULL" @@ -103,21 +104,30 @@ module program_setup character(len=500), public :: mosaic_file_input_grid = "NULL" character(len=500), public :: mosaic_file_target_grid = "NULL" character(len=500), public :: nst_files_input_grid = "NULL" + character(len=500), public :: grib2_file_input_grid = "NULL" character(len=500), public :: orog_dir_input_grid = "NULL" character(len=500), public :: orog_files_input_grid(6) = "NULL" character(len=500), public :: orog_dir_target_grid = "NULL" character(len=500), public :: orog_files_target_grid(6) = "NULL" character(len=500), public :: sfc_files_input_grid(6) = "NULL" character(len=500), public :: vcoord_file_target_grid = "NULL" - character(len=6), public :: cres_target_grid = " " + character(len=6), public :: cres_target_grid = "NULL" character(len=500), public :: atm_weight_file="NULL" character(len=20), public :: input_type="restart" + character(len=20), public :: phys_suite="GFS" !Default to gfs physics suite integer, parameter, public :: max_tracers=100 - integer, public :: num_tracers + integer, public :: num_tracers, num_tracers_input + + logical, allocatable, public :: read_from_input(:) + character(len=20), public :: tracers(max_tracers)="NULL" character(len=20), public :: tracers_input(max_tracers)="NULL" - + character(len=20), allocatable, public :: missing_var_methods(:) + character(len=20), allocatable, public :: chgres_var_names(:) + character(len=20), allocatable, public :: field_var_names(:) + + integer, public :: cycle_mon = -999 integer, public :: cycle_day = -999 integer, public :: cycle_hour = -999 @@ -129,24 +139,33 @@ module program_setup logical, public :: convert_nst = .false. logical, public :: convert_sfc = .false. + real, allocatable, public :: drysmc_input(:), drysmc_target(:) real, allocatable, public :: maxsmc_input(:), maxsmc_target(:) real, allocatable, public :: refsmc_input(:), refsmc_target(:) real, allocatable, public :: wltsmc_input(:), wltsmc_target(:) real, allocatable, public :: bb_target(:), satpsi_target(:) + real, allocatable, public :: missing_var_values(:) + public :: read_setup_namelist public :: calc_soil_params_driver + public :: read_varmap + public :: get_var_cond contains subroutine read_setup_namelist - + implicit none + + integer :: is, ie, ierr - namelist /config/ mosaic_file_target_grid, & + + namelist /config/ varmap_file, & + mosaic_file_target_grid, & fix_dir_target_grid, & orog_dir_target_grid, & orog_files_target_grid, & @@ -158,6 +177,7 @@ subroutine read_setup_namelist atm_files_input_grid, & atm_core_files_input_grid, & atm_tracer_files_input_grid, & + grib2_file_input_grid, & data_dir_input_grid, & vcoord_file_target_grid, & cycle_mon, cycle_day, & @@ -165,7 +185,9 @@ subroutine read_setup_namelist convert_nst, convert_sfc, & regional, input_type, & atm_weight_file, tracers, & - tracers_input, halo_bndy, halo_blend + tracers_input,phys_suite, & + halo_bndy, & + halo_blend print*,"- READ SETUP NAMELIST" @@ -174,7 +196,10 @@ subroutine read_setup_namelist read(41, nml=config, iostat=ierr) if (ierr /= 0) call error_handler("READING SETUP NAMELIST.", ierr) close (41) - + + call to_lower(input_type) + call to_upper(phys_suite) + orog_dir_target_grid = trim(orog_dir_target_grid) // '/' orog_dir_input_grid = trim(orog_dir_input_grid) // '/' @@ -217,6 +242,13 @@ subroutine read_setup_namelist num_tracers = num_tracers + 1 print*,"- WILL PROCESS TRACER ", trim(tracers(is)) enddo + + num_tracers_input = 0 + do is = 1, max_tracers + if (trim(tracers_input(is)) == "NULL") exit + num_tracers_input = num_tracers_input + 1 + print*,"- WILL PROCESS INPUT TRACER ", trim(tracers_input(is)) + enddo !------------------------------------------------------------------------- ! Ensure program recognizes the input data type. @@ -233,14 +265,118 @@ subroutine read_setup_namelist print*,'- INPUT DATA FROM SPECTRAL GFS GAUSSIAN NEMSIO FILE.' case ("gfs_spectral") print*,'- INPUT DATA FROM SPECTRAL GFS SIGIO/SFCIO FILE.' + case ("grib2") + print*,'- INPUT DATA FROM A GRIB2 FILE' case default call error_handler("UNRECOGNIZED INPUT DATA TYPE.", 1) end select + +!------------------------------------------------------------------------- +! Ensure proper file variable provided for grib2 input +!------------------------------------------------------------------------- - return - + if (trim(input_type) == "grib2") then + if (trim(grib2_file_input_grid) == "NULL" .or. trim(grib2_file_input_grid) == "") then + call error_handler("FOR GRIB2 DATA, PLEASE PROVIDE GRIB2_FILE_INPUT_GRID") + endif + endif + end subroutine read_setup_namelist +subroutine read_varmap + + implicit none + + integer :: istat, k, nvars + character(len=500) :: line + character(len=20),allocatable :: var_type(:) + + if (trim(input_type) == "grib2") then + + print*,"OPEN VARIABLE MAPPING FILE: ", trim(varmap_file) + open(14, file=trim(varmap_file), form='formatted', iostat=istat) + if (istat /= 0) then + call error_handler("OPENING VARIABLE MAPPING FILE", istat) + endif + + num_tracers = 0 + nvars = 0 + + !Loop over lines of file to count the number of variables + do + read(14, '(A)', iostat=istat) line !chgres_var_names_tmp(k)!, field_var_names(k) , & + ! missing_var_methods(k), missing_var_values(k), var_type(k) + if (istat/=0) exit + if ( trim(line) .eq. '' ) cycle + nvars = nvars+1 + enddo + + + allocate(chgres_var_names(nvars)) + allocate(field_var_names(nvars)) + allocate(missing_var_methods(nvars)) + allocate(missing_var_values(nvars)) + allocate(read_from_input(nvars)) + allocate(var_type(nvars)) + + read_from_input(:) = .true. + rewind(14) + do k = 1,nvars + read(14, *, iostat=istat) chgres_var_names(k), field_var_names(k) , & + missing_var_methods(k), missing_var_values(k), var_type(k) + if (istat /= 0) call error_handler("READING VARIABLE MAPPING FILE", istat) + if(trim(var_type(k))=='T') then + num_tracers = num_tracers + 1 + tracers_input(num_tracers)=chgres_var_names(k) + endif + enddo + close(14) + endif +end subroutine read_varmap + +! ---------------------------------------------------------------------------------------- +! Find conditions for handling missing variables from varmap arrays +! ---------------------------------------------------------------------------------------- + +subroutine get_var_cond(var_name,this_miss_var_method,this_miss_var_value, & + this_field_var_name, loc) + use esmf + + implicit none + character(len=20) :: var_name + + character(len=20), optional, intent(out) :: this_miss_var_method, & + this_field_var_name + real(esmf_kind_r4), optional, intent(out):: this_miss_var_value + + integer, optional, intent(out) :: loc + + integer :: i, tmp(size(missing_var_methods)) + + i=0 + + tmp(:)=0 + where(chgres_var_names==var_name) tmp=1 + + i = maxloc(merge(1.,0.,chgres_var_names == var_name),dim=1) !findloc(chgres_var_names,var_name) + print*, i + if (maxval(tmp).eq.0) then + print*, "WARNING: NO ENTRY FOR ", trim(var_name), " IN VARMAP TABLE. WILL SKIP " // & + "VARIABLE IF NOT FOUND IN EXTERNAL MODEL FILE" + + if(present(this_miss_var_method)) this_miss_var_method = "skip" + if(present(this_miss_var_value)) this_miss_var_value = -9999.9_esmf_kind_r4 + if(present(this_field_var_name)) this_field_var_name = "NULL" + if(present(loc)) loc = 9999 + else + if(present(this_miss_var_method)) this_miss_var_method = missing_var_methods(i) + if(present(this_miss_var_value)) this_miss_var_value = missing_var_values(i) + if(present(this_field_var_name)) this_field_var_name = field_var_names(i) + if(present(loc)) loc = i + endif + +end subroutine get_var_cond + subroutine calc_soil_params_driver(localpet) implicit none diff --git a/sorc/chgres_cube.fd/search_util.f90 b/sorc/chgres_cube.fd/search_util.f90 index 7f67286f9..0810b0b9b 100644 --- a/sorc/chgres_cube.fd/search_util.f90 +++ b/sorc/chgres_cube.fd/search_util.f90 @@ -133,7 +133,7 @@ subroutine search (field, mask, idim, jdim, tile, field_num, latitude) if (mask(ii,jj) == 1 .and. field_save(ii,jj) > -9999.0) then field(i,j) = field_save(ii,jj) -! write(6,100) tile,i,j,ii,jj,field(i,j) + write(6,100) tile,i,j,ii,jj,field(i,j) cycle I_LOOP endif diff --git a/sorc/chgres_cube.fd/surface.F90 b/sorc/chgres_cube.fd/surface.F90 index 7eecf12fa..02e64dd1b 100644 --- a/sorc/chgres_cube.fd/surface.F90 +++ b/sorc/chgres_cube.fd/surface.F90 @@ -292,9 +292,10 @@ subroutine interp(localpet) seamask_target_grid, & latitude_target_grid - use program_setup, only : convert_nst + use program_setup, only : convert_nst, input_type - use static_data, only : veg_type_target_grid + use static_data, only : veg_type_target_grid, & + soil_type_target_grid use search_util @@ -318,6 +319,7 @@ subroutine interp(localpet) integer(esmf_kind_i8), pointer :: seamask_target_ptr(:,:) real(esmf_kind_r8), allocatable :: data_one_tile(:,:) + real(esmf_kind_r8), allocatable :: data_one_tile2(:,:) real(esmf_kind_r8), allocatable :: data_one_tile_3d(:,:,:) real(esmf_kind_r8), allocatable :: latitude_one_tile(:,:) real(esmf_kind_r8), pointer :: canopy_mc_target_ptr(:,:) @@ -514,6 +516,18 @@ subroutine interp(localpet) ! Interpolate. !----------------------------------------------------------------------- + if (localpet == 0) then + allocate(data_one_tile(i_target,j_target)) + allocate(data_one_tile2(i_target,j_target)) + allocate(data_one_tile_3d(i_target,j_target,lsoil_target)) + allocate(mask_target_one_tile(i_target,j_target)) + else + allocate(data_one_tile(0,0)) + allocate(data_one_tile2(0,0)) + allocate(data_one_tile_3d(0,0,0)) + allocate(mask_target_one_tile(0,0)) + endif + method=ESMF_REGRIDMETHOD_CONSERVE isrctermprocessing = 1 @@ -541,16 +555,6 @@ subroutine interp(localpet) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldRegrid", rc) - if (localpet == 0) then - allocate(data_one_tile(i_target,j_target)) - allocate(data_one_tile_3d(i_target,j_target,lsoil_target)) - allocate(mask_target_one_tile(i_target,j_target)) - else - allocate(data_one_tile(0,0)) - allocate(data_one_tile_3d(0,0,0)) - allocate(mask_target_one_tile(0,0)) - endif - print*,"- CALL FieldGet FOR TARGET grid sea ice fraction." call ESMF_FieldGet(seaice_fract_target_grid, & farrayPtr=seaice_fract_target_ptr, rc=rc) @@ -600,6 +604,7 @@ subroutine interp(localpet) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGather", rc) + if (localpet == 0) then do j = 1, j_target do i = 1, i_target @@ -1969,7 +1974,7 @@ subroutine interp(localpet) do ij = l(1), u(1) call ij_to_i_j(unmapped_ptr(ij), i_target, j_target, i, j) - soilm_tot_target_ptr(i,j,:) = -9999.9 + soilm_tot_target_ptr(i,j,:) = -9999.9 soil_temp_target_ptr(i,j,:) = -9999.9 skin_temp_target_ptr(i,j) = -9999.9 terrain_from_input_ptr(i,j) = -9999.9 @@ -2028,8 +2033,23 @@ subroutine interp(localpet) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGather", rc) + print*,"- CALL FieldGather FOR SOIL TYPE TARGET GRID, TILE: ", tile + call ESMF_FieldGather(soil_type_target_grid, data_one_tile2, rootPet=0,tile=tile, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldGather", rc) + +!--------------------------------------------------------------------------------------- +! grib2 data does not have soil type. Set soil type interpolated from input +! grid to the target (model) grid soil type. This turns off the soil moisture +! rescaling. +!--------------------------------------------------------------------------------------- + if (localpet == 0) then - call search(data_one_tile, mask_target_one_tile, i_target, j_target, tile, 224) + if (trim(input_type) .ne. "grib2") then + call search(data_one_tile, mask_target_one_tile, i_target, j_target, tile, 224) + else + data_one_tile = data_one_tile2 + endif endif print*,"- CALL FieldScatter FOR SOIL TYPE FROM INPUT GRID, TILE: ", tile @@ -2082,7 +2102,7 @@ subroutine interp(localpet) deallocate(veg_type_target_one_tile) - deallocate(data_one_tile) + deallocate(data_one_tile, data_one_tile2) deallocate(data_one_tile_3d) deallocate(mask_target_one_tile) diff --git a/sorc/chgres_cube.fd/utils.f90 b/sorc/chgres_cube.fd/utils.f90 index a99b0b5d7..ddda0bb90 100644 --- a/sorc/chgres_cube.fd/utils.f90 +++ b/sorc/chgres_cube.fd/utils.f90 @@ -32,3 +32,45 @@ subroutine netcdf_err( err, string ) return end subroutine netcdf_err + + subroutine to_upper(strIn) +! Adapted from http://www.star.le.ac.uk/~cgp/fortran.html (25 May 2012) +! Original author: Clive Page + + implicit none + + character(len=*), intent(inout) :: strIn + character(len=len(strIn)) :: strOut + integer :: i,j + + do i = 1, len(strIn) + j = iachar(strIn(i:i)) + if (j>= iachar("a") .and. j<=iachar("z") ) then + strOut(i:i) = achar(iachar(strIn(i:i))-32) + else + strOut(i:i) = strIn(i:i) + end if + end do + strIn(:) = strOut(:) +end subroutine to_upper + +subroutine to_lower(strIn) +! Adapted from http://www.star.le.ac.uk/~cgp/fortran.html (25 May 2012) +! Original author: Clive Page + + implicit none + + character(len=*), intent(inout) :: strIn + character(len=len(strIn)) :: strOut + integer :: i,j + + do i = 1, len(strIn) + j = iachar(strIn(i:i)) + if (j>= iachar("A") .and. j<=iachar("Z") ) then + strOut(i:i) = achar(iachar(strIn(i:i))+32) + else + strOut(i:i) = strIn(i:i) + end if + end do + strIn(:) = strOut(:) +end subroutine to_lower diff --git a/sorc/machine-setup.sh b/sorc/machine-setup.sh new file mode 100644 index 000000000..48b2ee1bf --- /dev/null +++ b/sorc/machine-setup.sh @@ -0,0 +1,122 @@ +# Create a test function for sh vs. bash detection. The name is +# randomly generated to reduce the chances of name collision. +__ms_function_name="setup__test_function__$$" +eval "$__ms_function_name() { /bin/true ; }" + +# Determine which shell we are using +__ms_ksh_test=$( eval '__text="text" ; if [[ $__text =~ ^(t).* ]] ; then printf "%s" ${.sh.match[1]} ; fi' 2> /dev/null | cat ) +__ms_bash_test=$( eval 'if ( set | grep '$__ms_function_name' | grep -v name > /dev/null 2>&1 ) ; then echo t ; fi ' 2> /dev/null | cat ) + +if [[ ! -z "$__ms_ksh_test" ]] ; then + __ms_shell=ksh +elif [[ ! -z "$__ms_bash_test" ]] ; then + __ms_shell=bash +else + # Not bash or ksh, so assume sh. + __ms_shell=sh +fi + +target="" +USERNAME=`echo $LOGNAME | awk '{ print tolower($0)'}` + +if [[ -d /lfs3 ]] ; then + # We are on NOAA Jet + if ( ! eval module help > /dev/null 2>&1 ) ; then + echo load the module command 1>&2 + source /apps/lmod/lmod/init/$__ms_shell + fi + target=jet + module purge + module use /lfs3/projects/hfv3gfs/nwprod/NCEPLIBS/modulefiles +elif [[ -d /scratch1 ]] ; then + # We are on NOAA Hera + if ( ! eval module help > /dev/null 2>&1 ) ; then + echo load the module command 1>&2 + source /apps/lmod/lmod/init/$__ms_shell + fi + target=hera + module purge + MOD_PATH=/scratch2/NCEPDEV/nwprod/NCEPLIBS/modulefiles +elif [[ -d /gpfs/hps && -e /etc/SuSE-release ]] ; then + # We are on NOAA Luna or Surge + if ( ! eval module help > /dev/null 2>&1 ) ; then + echo load the module command 1>&2 + source /opt/modules/default/init/$__ms_shell + fi + target=wcoss_cray + + # Silence the "module purge" to avoid the expected error messages + # related to modules that load modules. + module purge > /dev/null 2>&1 + module use /usrx/local/prod/modulefiles + module use /gpfs/hps/nco/ops/nwprod/lib/modulefiles + module use /gpfs/hps/nco/ops/nwprod/modulefiles + module use /opt/cray/alt-modulefiles + module use /opt/cray/craype/default/alt-modulefiles + module use /opt/cray/ari/modulefiles + module use /opt/modulefiles + module purge > /dev/null 2>&1 + + # Workaround until module issues are fixed: + #unset _LMFILES_ + #unset LOADEDMODULES + echo y 2> /dev/null | module clear > /dev/null 2>&1 + + module use /usrx/local/prod/modulefiles + module use /gpfs/hps/nco/ops/nwprod/lib/modulefiles + module use /gpfs/hps/nco/ops/nwprod/modulefiles + module use /opt/cray/alt-modulefiles + module use /opt/cray/craype/default/alt-modulefiles + module use /opt/cray/ari/modulefiles + module use /opt/modulefiles + module load modules + +elif [[ -L /usrx && "$( readlink /usrx 2> /dev/null )" =~ dell ]] ; then + # We are on NOAA Venus or Mars + if ( ! eval module help > /dev/null 2>&1 ) ; then + echo load the module command 1>&2 + source /usrx/local/prod/lmod/lmod/init/$__ms_shell + fi + target=wcoss_dell_p3 + module purge + +elif [[ -d /dcom && -d /hwrf ]] ; then + # We are on NOAA Tide or Gyre + if ( ! eval module help > /dev/null 2>&1 ) ; then + echo load the module command 1>&2 + source /usrx/local/Modules/default/init/$__ms_shell + fi + target=wcoss + module purge +elif [[ -d /glade ]] ; then + # We are on NCAR Yellowstone + if ( ! eval module help > /dev/null 2>&1 ) ; then + echo load the module command 1>&2 + . /usr/share/Modules/init/$__ms_shell + fi + target=yellowstone + module purge +elif [[ -d /lustre && -d /ncrc ]] ; then + # We are on GAEA. + if ( ! eval module help > /dev/null 2>&1 ) ; then + # We cannot simply load the module command. The GAEA + # /etc/profile modifies a number of module-related variables + # before loading the module command. Without those variables, + # the module command fails. Hence we actually have to source + # /etc/profile here. + echo load the module command 1>&2 + source /etc/profile + fi + target=gaea + module purge +elif [[ "$(hostname)" =~ "odin" ]]; then + target="odin" +else + echo WARNING: UNKNOWN PLATFORM 1>&2 +fi + +unset __ms_shell +unset __ms_ksh_test +unset __ms_bash_test +unset $__ms_function_name +unset __ms_function_name diff --git a/sorc/sfc_climo_gen.fd/driver.F90 b/sorc/sfc_climo_gen.fd/driver.F90 new file mode 100644 index 000000000..5e5207190 --- /dev/null +++ b/sorc/sfc_climo_gen.fd/driver.F90 @@ -0,0 +1,171 @@ + program driver + +!------------------------------------------------------------------------- +! program documentation block +! +! Abstract: Reads static surface data on a global lat/lon grid, +! interpolates the data to the fv3 model grid, and outputs the +! result in netcdf format. +! +! Program execution is controlled by variables defined in the +! program configuration namelist (see module program_setup for +! details). +! +! Requires the following input files: +! 1) Model mosaic file (netcdf) +! 2) Model orography files (netcdf) +! 3) Model grid files (netcdf) +! 4) Source data file on global lat/lon grid (netcdf) +! +! Outputs surface data on the model tiles in netcdf format. +!------------------------------------------------------------------------- + + use model_grid + use source_grid + use program_setup + use utils + use esmf + + implicit none + + integer :: rc + integer :: localpet + integer :: npets + + type(esmf_vm) :: vm + type(esmf_regridmethod_flag) :: method + +!------------------------------------------------------------------------- +! Initialize MPI and ESMF environments. +!------------------------------------------------------------------------- + + call mpi_init(rc) + + print*,"- INITIALIZE ESMF" + call ESMF_Initialize(rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("INITIALIZING ESMF", rc) + + print*,"- CALL VMGetGlobal" + call ESMF_VMGetGlobal(vm, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN VMGetGlobal.", rc) + + print*,"- CALL VMGet" + call ESMF_VMGet(vm, localPet=localpet, petCount=npets, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN VMGet.", rc) + + print*,'- NPETS IS ',npets + print*,'- LOCAL PET ',localpet + +!------------------------------------------------------------------------- +! Program set up step. +!------------------------------------------------------------------------- + + call read_setup_namelist(localpet) + +!------------------------------------------------------------------------- +! Define fv3 model grid. +!------------------------------------------------------------------------- + + call define_model_grid(localpet, npets) + +!------------------------------------------------------------------------- +! Interpolate fields. Vegetation type must be done first as +! it defines which points are permanent land ice. +!------------------------------------------------------------------------- + + call define_source_grid(localpet, npets, input_vegetation_type_file) + method=ESMF_REGRIDMETHOD_NEAREST_STOD + call interp(localpet, method, input_vegetation_type_file) + call source_grid_cleanup + +! Snow free albedo + + if (trim(input_snowfree_albedo_file) /= "NULL") then + call define_source_grid(localpet, npets, input_snowfree_albedo_file) + method=ESMF_REGRIDMETHOD_BILINEAR + if (trim(snowfree_albedo_method)=="conserve") method=ESMF_REGRIDMETHOD_CONSERVE + call interp(localpet, method, input_snowfree_albedo_file) + call source_grid_cleanup + endif + +! Maximum snow albedo + + if (trim(input_maximum_snow_albedo_file) /= "NULL") then + call define_source_grid(localpet, npets, input_maximum_snow_albedo_file) + method=ESMF_REGRIDMETHOD_BILINEAR + if (trim(maximum_snow_albedo_method)=="conserve") method=ESMF_REGRIDMETHOD_CONSERVE + call interp(localpet, method, input_maximum_snow_albedo_file) + call source_grid_cleanup + endif + +! FACSF - fractional coverage for strong zenith angle +! dependant albedo. + + if (trim(input_facsf_file) /= "NULL") then + call define_source_grid(localpet, npets, input_facsf_file) + method=ESMF_REGRIDMETHOD_BILINEAR + call interp(localpet, method, input_facsf_file) + call source_grid_cleanup + endif + +! Soil substrate temperature + + if (trim(input_substrate_temperature_file) /= "NULL") then + call define_source_grid(localpet, npets, input_substrate_temperature_file) + method=ESMF_REGRIDMETHOD_BILINEAR + call interp(localpet, method, input_substrate_temperature_file) + call source_grid_cleanup + endif + +! Slope type + + if (trim(input_slope_type_file) /= "NULL") then + call define_source_grid(localpet, npets, input_slope_type_file) + method=ESMF_REGRIDMETHOD_NEAREST_STOD + call interp(localpet, method, input_slope_type_file) + call source_grid_cleanup + endif + +! Soil type + + if (trim(input_soil_type_file) /= "NULL") then + call define_source_grid(localpet, npets, input_soil_type_file) + method=ESMF_REGRIDMETHOD_NEAREST_STOD + call interp(localpet, method, input_soil_type_file) + call source_grid_cleanup + endif + +! Vegetation greenness + + if (trim(input_vegetation_greenness_file) /= "NULL") then + call define_source_grid(localpet, npets, input_vegetation_greenness_file) + method=ESMF_REGRIDMETHOD_BILINEAR + if (trim(vegetation_greenness_method)=="conserve") method=ESMF_REGRIDMETHOD_CONSERVE + call interp(localpet, method, input_vegetation_greenness_file) + call source_grid_cleanup + endif + +! Leaf Area Index + + if (trim(input_leaf_area_index_file) /= "NULL") then + call define_source_grid(localpet, npets, input_leaf_area_index_file) + method=ESMF_REGRIDMETHOD_BILINEAR + if (trim(leaf_area_index_method)=="conserve") method=ESMF_REGRIDMETHOD_CONSERVE + call interp(localpet, method, input_leaf_area_index_file) + call source_grid_cleanup + endif + + call model_grid_cleanup + + print*,"- CALL ESMF_finalize" + call ESMF_finalize(endflag=ESMF_END_KEEPMPI, rc=rc) + + call mpi_finalize(rc) + + print*,'- DONE.' + stop + + end program driver diff --git a/sorc/sfc_climo_gen.fd/interp.F90 b/sorc/sfc_climo_gen.fd/interp.F90 new file mode 100644 index 000000000..5ed80f14d --- /dev/null +++ b/sorc/sfc_climo_gen.fd/interp.F90 @@ -0,0 +1,354 @@ + subroutine interp(localpet, method, input_file) + +!----------------------------------------------------------------------- +! subroutine documentation block +! +! Subroutine: interp +! prgmmr: gayno org: w/np2 date: 2018 +! +! Abstract: Read the input source data and interpolate it to the +! model grid. +! +! Usage: call interp(localpet, method, input_file) +! +! input argument list: +! input_flle filename of input source data. +! localpet this mpi task +! method interpolation method.defined where mask=1 +! +!----------------------------------------------------------------------- + + use esmf + use netcdf + use model_grid + use source_grid + use utils + + implicit none + + include 'mpif.h' + + character(len=*), intent(in) :: input_file + + integer :: rc, localpet + integer :: i, j, ij, tile, n, ncid, status + integer :: l(1), u(1), t + integer :: clb_mdl(2), cub_mdl(2) + integer :: varid, record + integer :: tile_num, pt_loc_this_tile + integer :: isrctermprocessing + + integer(esmf_kind_i4), allocatable :: mask_mdl_one_tile(:,:) + integer(esmf_kind_i4), pointer :: unmapped_ptr(:) + + real(esmf_kind_r4), pointer :: data_mdl_ptr(:,:) + real(esmf_kind_r4), pointer :: data_src_ptr(:,:) + real(esmf_kind_r4), allocatable :: data_src_global(:,:) + real(esmf_kind_r4), allocatable :: data_mdl_one_tile(:,:) + real(esmf_kind_r4), allocatable :: vegt_mdl_one_tile(:,:) + + type(esmf_regridmethod_flag),intent(in) :: method + type(esmf_field) :: data_field_src + type(esmf_routehandle) :: regrid_data + type(esmf_polemethod_flag) :: pole + + print*,"- CALL FieldCreate FOR SOURCE GRID DATA." + data_field_src = ESMF_FieldCreate(grid_src, & + typekind=ESMF_TYPEKIND_R4, & + indexflag=ESMF_INDEX_GLOBAL, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + name="source data", & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + print*,"- CALL FieldGet FOR SOURCE GRID DATA." + nullify(data_src_ptr) + call ESMF_FieldGet(data_field_src, & + farrayPtr=data_src_ptr, & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + print*,"- CALL FieldGet FOR MODEL GRID DATA." + nullify(data_mdl_ptr) + call ESMF_FieldGet(data_field_mdl, & + farrayPtr=data_mdl_ptr, & + computationalLBound=clb_mdl, & + computationalUBound=cub_mdl, & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + if (localpet == 0) then + allocate(data_src_global(i_src,j_src)) + else + allocate(data_src_global(0,0)) + endif + + print*,'- OPEN SOURCE FILE ', trim(input_file) + status = nf90_open(input_file, nf90_nowrite, ncid) + call netcdf_err(status, "IN ROUTINE INTERP OPENING SOURCE FILE") + + if (localpet == 0) then + allocate(data_mdl_one_tile(i_mdl,j_mdl)) + allocate(mask_mdl_one_tile(i_mdl,j_mdl)) + else + allocate(data_mdl_one_tile(0,0)) + allocate(mask_mdl_one_tile(0,0)) + endif + + record = 0 + + TIME_LOOP : do t = 1, num_time_recs ! loop over each time period + + FIELD_LOOP : do n = 1, num_fields ! loop over each surface field. + + record = record + 1 + + if (localpet == 0) then + status = nf90_inq_varid(ncid, field_names(n), varid) + call netcdf_err(status, "IN ROUTINE INTERP READING FIELD ID") + status = nf90_get_var(ncid, varid, data_src_global, start=(/1,1,t/), count=(/i_src,j_src,1/)) + call netcdf_err(status, "IN ROUTINE INTERP READING FIELD") + endif + + print*,"- CALL FieldScatter FOR SOURCE GRID DATA." + call ESMF_FieldScatter(data_field_src, data_src_global, rootPet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter.", rc) + + isrctermprocessing = 1 + + if (record == 1) then + + if (method == ESMF_REGRIDMETHOD_BILINEAR) then + pole = ESMF_POLEMETHOD_ALLAVG + else + pole = ESMF_POLEMETHOD_NONE + endif + + print*,"- CALL FieldRegridStore." + nullify(unmapped_ptr) + call ESMF_FieldRegridStore(data_field_src, & + data_field_mdl, & + srcmaskvalues=(/0/), & + dstmaskvalues=(/0/), & + polemethod=pole, & + unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, & + normtype=ESMF_NORMTYPE_FRACAREA, & + srctermprocessing=isrctermprocessing, & + routehandle=regrid_data, & + regridmethod=method, & + unmappedDstList=unmapped_ptr, & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldRegridStore.", rc) + + endif + + print*,"- CALL Field_Regrid." + call ESMF_FieldRegrid(data_field_src, & + data_field_mdl, & + routehandle=regrid_data, & + termorderflag=ESMF_TERMORDER_SRCSEQ, & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldRegrid.", rc) + +!----------------------------------------------------------------------- +! Unmapped points are stored in "unmapped_ptr". The pointer contains +! "ij" global indices as follows: +! +! tile 1: 1 thru (itile*jtile) +! tile n: (n-1)*(itile*jtile) thru n*(itile*jtile) +! +! This "ij" index is converted to the tile number and i/j index of the +! field object. This logic assumes the model grid object was +! created using "GLOBAL" indices. +! +! Unmapped data points are given the flag value of -9999.9, which +! will be replaced in routine "search". +!----------------------------------------------------------------------- + + l = lbound(unmapped_ptr) + u = ubound(unmapped_ptr) + do ij = l(1), u(1) + + tile_num = ((unmapped_ptr(ij)-1) / (i_mdl*j_mdl)) ! tile number minus 1 + pt_loc_this_tile = unmapped_ptr(ij) - (tile_num * i_mdl * j_mdl) + ! "ij" location of point within tile. + + j = (pt_loc_this_tile - 1) / i_mdl + 1 + i = mod(pt_loc_this_tile, i_mdl) + if (i==0) i = i_mdl + data_mdl_ptr(i,j) = -9999.9 + + enddo + +! These fields are adjusted at landice. + + select case (trim(field_names(n))) + case ('substrate_temperature','vegetation_greenness','leaf_area_index','slope_type','soil_type') + if (localpet == 0) then + allocate(vegt_mdl_one_tile(i_mdl,j_mdl)) + else + allocate(vegt_mdl_one_tile(0,0)) + endif + end select + + OUTPUT_LOOP : do tile = 1, num_tiles + + print*,"- CALL FieldGather FOR MODEL GRID DATA." + call ESMF_FieldGather(data_field_mdl, data_mdl_one_tile, rootPet=0, tile=tile, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather.", rc) + + print*,"- CALL FieldGather FOR MODEL GRID LAND MASK." + call ESMF_FieldGather(mask_field_mdl, mask_mdl_one_tile, rootPet=0, tile=tile, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather.", rc) + + select case (trim(field_names(n))) + case ('substrate_temperature','vegetation_greenness','leaf_area_index','slope_type','soil_type') + print*,"- CALL FieldGather FOR MODEL GRID VEG TYPE." + call ESMF_FieldGather(vegt_field_mdl, vegt_mdl_one_tile, rootPet=0, tile=tile, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather.", rc) + end select + + if (localpet == 0) then + print*,'- CALL SEARCH FOR TILE ',tile + call search (data_mdl_one_tile, mask_mdl_one_tile, i_mdl, j_mdl, tile, field_names(n)) + select case (field_names(n)) + case ('substrate_temperature','vegetation_greenness','leaf_area_index','slope_type','soil_type') + call adjust_for_landice (data_mdl_one_tile, vegt_mdl_one_tile, i_mdl, j_mdl, field_names(n)) + end select + where(mask_mdl_one_tile == 0) data_mdl_one_tile = missing + call output (data_mdl_one_tile, i_mdl, j_mdl, tile, record, t, n) + endif + + if (field_names(n) == 'vegetation_type') then + print*,"- CALL FieldScatter FOR MODEL GRID VEGETATION TYPE." + call ESMF_FieldScatter(vegt_field_mdl, data_mdl_one_tile, rootPet=0, tile=tile, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter.", rc) + endif + + enddo OUTPUT_LOOP + + if (allocated(vegt_mdl_one_tile)) deallocate(vegt_mdl_one_tile) + + enddo FIELD_LOOP + enddo TIME_LOOP + + status=nf90_close(ncid) + + deallocate(data_mdl_one_tile, mask_mdl_one_tile) + deallocate(data_src_global) + + print*,"- CALL FieldRegridRelease." + call ESMF_FieldRegridRelease(routehandle=regrid_data, rc=rc) + + end subroutine interp + + subroutine adjust_for_landice(field, vegt, idim, jdim, field_ch) + +!----------------------------------------------------------------------- +! subroutine documentation block +! +! Subroutine: adjust for landice +! prgmmr: gayno org: w/np2 date: 2018 +! +! Abstract: Ensure consistent fields at land ice points. +! Land ice is vegetation type 15 (variable landice). +! +! Usage: call adjust_for_landice(field, vegt, idim, jdim, field_ch) +! +! input argument list: +! field Model field before adjustments for +! land ice. +! field_ch Field name +! i/jdim i/j dimension of model tile. +! vegt Vegetation type on the model tile +! +! output argument list: +! field Model field after adjustments for +! land ice. +!----------------------------------------------------------------------- + + use esmf + + implicit none + + include 'mpif.h' + + character(len=*), intent(in) :: field_ch + + integer, intent(in) :: idim, jdim + + real(esmf_kind_i4), intent(in) :: vegt(idim,jdim) + real(esmf_kind_r4), intent(inout) :: field(idim,jdim) + + integer, parameter :: landice=15 + + integer :: i, j + + real :: landice_value + + select case (field_ch) + case ('substrate_temperature') ! soil substrate temp + landice_value = 273.15 + do j = 1, jdim + do i = 1, idim + if (nint(vegt(i,j)) == landice) then + field(i,j) = min(field(i,j), landice_value) + endif + enddo + enddo + case ('vegetation_greenness') ! vegetation greenness + landice_value = 0.01 ! 1.0% is bare ground + do j = 1, jdim + do i = 1, idim + if (nint(vegt(i,j)) == landice) then + field(i,j) = landice_value + endif + enddo + enddo + case ('leaf_area_index') ! leaf area index + landice_value = 0.0 ! bare ground + do j = 1, jdim + do i = 1, idim + if (nint(vegt(i,j)) == landice) then + field(i,j) = landice_value + endif + enddo + enddo + case ('slope_type') ! slope type + landice_value = 9.0 + do j = 1, jdim + do i = 1, idim + if (nint(vegt(i,j)) == landice) then + field(i,j) = landice_value + else + if (nint(field(i,j)) == nint(landice_value)) field(i,j) = 2.0 + endif + enddo + enddo + case ('soil_type') ! soil type + landice_value = 16.0 + do j = 1, jdim + do i = 1, idim + if (nint(vegt(i,j)) == landice) then + field(i,j) = landice_value + else + if (nint(field(i,j)) == nint(landice_value)) field(i,j) = 6.0 + endif + enddo + enddo + case default + print*,'- FATAL ERROR IN ROUTINE ADJUST_FOR_LANDICE. UNIDENTIFIED FIELD : ', field_ch + call mpi_abort(mpi_comm_world, 57) + end select + + end subroutine adjust_for_landice diff --git a/sorc/sfc_climo_gen.fd/model_grid.F90 b/sorc/sfc_climo_gen.fd/model_grid.F90 new file mode 100644 index 000000000..f281b64ea --- /dev/null +++ b/sorc/sfc_climo_gen.fd/model_grid.F90 @@ -0,0 +1,366 @@ + module model_grid + +!-------------------------------------------------------------------------- +! module documentation block +! +! Module: model grid +! pgrmmr: gayno org: w/np2 date: 2018 +! +! Abstract: Defines the model grid. +! +! Usage: use model_grid +! +! Public Subroutines: +! ------------------- +! define_model_grid Defines esmf grid object for the +! model grid. +! model_grid_cleanup Free up memory used in this module. +! +! Public variables: +! ----------------- +! +! Variables named with 'mdl' refer to the model grid. +! +! data_field_mdl ESMF field object that holds the +! data interpolated to model grid. +! grid_mdl ESMF grid object for the model grid. +! grid_tiles Array of model grid tile names. +! i/j_mdl i/j dimensions of model tile. +! mdl_field_mdl ESMF field object that holds the +! model land mask. +! missing Value assigned to undefined points +! (i.e., ocean points for a land +! field). +! num_tiles Total number of model grid tiles. +! vegt_field_mdl ESMF field object that holds the +! vegetation type on the model grid. +! +!----------------------------------------------------------------------- + + use esmf + + implicit none + + private + + character(len=5), allocatable, public :: grid_tiles(:) + + integer, public :: i_mdl, j_mdl, ij_mdl, num_tiles + + real(kind=4), public :: missing = -999. + + type(esmf_grid), public :: grid_mdl + type(esmf_field), public :: data_field_mdl, mask_field_mdl + type(esmf_field), public :: vegt_field_mdl + + public :: define_model_grid + public :: model_grid_cleanup + + contains + + subroutine define_model_grid(localpet, npets) + +!----------------------------------------------------------------------- +! subroutine documentation block +! +! Subroutine: define model grid +! prgmmr: gayno org: w/np2 date: 2018 +! +! Abstract: Define the model grid from the mosaic and orography +! files. Create the ESMF grid object for the model grid. +! +! Usage: define_model_grid(localpet, npets) +! +! input argument list: +! localpet this mpi task +! npets total number of mpi tasks +!----------------------------------------------------------------------- + + use esmf + use netcdf + use program_setup + use utils + + implicit none + + integer, intent(in) :: localpet, npets + + character(len=500) :: the_file + + integer :: error, id_dim, id_tiles, ncid + integer :: id_grid_tiles + integer :: extra, rc, tile + integer, allocatable :: decomptile(:,:) + + integer, allocatable :: mask_mdl_one_tile(:,:) + integer(esmf_kind_i4), pointer :: mask_field_mdl_ptr(:,:) + integer(esmf_kind_i4), pointer :: mask_mdl_ptr(:,:) + +!----------------------------------------------------------------------- +! Get the number of tiles from the mosaic file. +!----------------------------------------------------------------------- + + print*,'- OPEN MODEL GRID MOSAIC FILE: ',trim(mosaic_file_mdl) + error=nf90_open(trim(mosaic_file_mdl),nf90_nowrite,ncid) + call netcdf_err(error, "OPENING MODEL GRID MOSAIC FILE") + + print*,"- READ NUMBER OF TILES" + error=nf90_inq_dimid(ncid, 'ntiles', id_tiles) + call netcdf_err(error, "READING NTILES ID") + error=nf90_inquire_dimension(ncid,id_tiles,len=num_tiles) + call netcdf_err(error, "READING NTILES") + error=nf90_inq_varid(ncid, 'gridtiles', id_grid_tiles) + call netcdf_err(error, "READING GRIDTILES ID") + allocate(grid_tiles(num_tiles)) + grid_tiles="NULL" + print*,"- READ TILE NAMES" + error=nf90_get_var(ncid, id_grid_tiles, grid_tiles) + call netcdf_err(error, "READING GRIDTILES") + + error = nf90_close(ncid) + + print*,'- NUMBER OF TILES, MODEL GRID IS ', num_tiles + + if (mod(npets,num_tiles) /= 0) then + print*,'- FATAL ERROR: MUST RUN THIS PROGRAM WITH A TASK COUNT THAT' + print*,'- IS A MULTIPLE OF THE NUMBER OF TILES.' + call mpi_abort + endif + +!----------------------------------------------------------------------- +! Get the model grid specs and land mask from the orography files. +!----------------------------------------------------------------------- + + orog_dir_mdl = trim(orog_dir_mdl) // '/' + the_file = trim(orog_dir_mdl) // trim(orog_files_mdl(1)) + + print*,'- OPEN FIRST MODEL GRID OROGRAPHY FILE: ',trim(the_file) + error=nf90_open(trim(the_file),nf90_nowrite,ncid) + call netcdf_err(error, "OPENING MODEL GRID OROGRAPHY FILE") + print*,"- READ GRID DIMENSIONS" + error=nf90_inq_dimid(ncid, 'lon', id_dim) + call netcdf_err(error, "READING MODEL LON ID") + error=nf90_inquire_dimension(ncid,id_dim,len=i_mdl) + call netcdf_err(error, "READING MODEL LON") + error=nf90_inq_dimid(ncid, 'lat', id_dim) + call netcdf_err(error, "READING MODEL LAT ID") + error=nf90_inquire_dimension(ncid,id_dim,len=j_mdl) + call netcdf_err(error, "READING MODEL LAT") + error = nf90_close(ncid) + + print*,"- I/J DIMENSIONS OF THE MODEL GRID TILES ", i_mdl, j_mdl + + ij_mdl = i_mdl * j_mdl + +!----------------------------------------------------------------------- +! Create ESMF grid object for the model grid. +!----------------------------------------------------------------------- + + extra = npets / num_tiles + + allocate(decomptile(2,num_tiles)) + + do tile = 1, num_tiles + decomptile(:,tile)=(/1,extra/) + enddo + + print*,"- CALL GridCreateMosaic FOR MODEL GRID" + grid_mdl = ESMF_GridCreateMosaic(filename=trim(mosaic_file_mdl), & + regDecompPTile=decomptile, & + staggerLocList=(/ESMF_STAGGERLOC_CENTER, & + ESMF_STAGGERLOC_CORNER/), & + indexflag=ESMF_INDEX_GLOBAL, & + tileFilePath=trim(orog_dir_mdl), & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridCreateMosaic", rc) + + print*,"- CALL FieldCreate FOR DATA INTERPOLATED TO MODEL GRID." + data_field_mdl = ESMF_FieldCreate(grid_mdl, & + typekind=ESMF_TYPEKIND_R4, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + name="data on model grid", & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + print*,"- CALL FieldCreate FOR VEGETATION TYPE INTERPOLATED TO MODEL GRID." + vegt_field_mdl = ESMF_FieldCreate(grid_mdl, & + typekind=ESMF_TYPEKIND_R4, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + name="veg type on model grid", & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + +!----------------------------------------------------------------------- +! Set model land mask. +!----------------------------------------------------------------------- + + print*,"- CALL FieldCreate FOR MODEL GRID LANDMASK." + mask_field_mdl = ESMF_FieldCreate(grid_mdl, & + typekind=ESMF_TYPEKIND_I4, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + name="model grid land mask", & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + print*,"- CALL FieldGet FOR MODEL GRID LANDMASK." + nullify(mask_field_mdl_ptr) + call ESMF_FieldGet(mask_field_mdl, & + farrayPtr=mask_field_mdl_ptr, & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + if (localpet == 0) then + allocate(mask_mdl_one_tile(i_mdl,j_mdl)) + else + allocate(mask_mdl_one_tile(0,0)) + endif + + do tile = 1, num_tiles + if (localpet == 0) then + the_file = trim(orog_dir_mdl) // trim(orog_files_mdl(tile)) + call get_model_mask(trim(the_file), mask_mdl_one_tile, i_mdl, j_mdl) + endif + print*,"- CALL FieldScatter FOR MODEL GRID MASK. TILE IS: ", tile + call ESMF_FieldScatter(mask_field_mdl, mask_mdl_one_tile, rootpet=0, tile=tile, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + enddo + + deallocate(mask_mdl_one_tile) + + print*,"- CALL GridAddItem FOR MODEL GRID." + call ESMF_GridAddItem(grid_mdl, & + itemflag=ESMF_GRIDITEM_MASK, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridAddItem", rc) + + print*,"- CALL GridGetItem FOR MODEL GRID." + nullify(mask_mdl_ptr) + call ESMF_GridGetItem(grid_mdl, & + itemflag=ESMF_GRIDITEM_MASK, & + farrayPtr=mask_mdl_ptr, & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetItem", rc) + + mask_mdl_ptr = mask_field_mdl_ptr + + end subroutine define_model_grid + + subroutine get_model_mask(orog_file, mask, idim, jdim) + +!----------------------------------------------------------------------- +! subroutine documentation block +! +! Subroutine: get model mask +! prgmmr: gayno org: w/np2 date: 2018 +! +! Abstract: Read model land/sea mask from the orography file. +! +! Usage: call get_model_mask(orog_file, mask, idim, jdim) +! +! input argument list: +! orog_file the orography file +! i/jdim i/j dimension of the model tile +! +! output argument list: +! mask land/sea mask +! +!----------------------------------------------------------------------- + + use esmf + use netcdf + use utils + + implicit none + + character(len=*), intent(in) :: orog_file + + integer, intent(in) :: idim, jdim + integer(esmf_kind_i4), intent(out) :: mask(idim,jdim) + + integer :: error, lat, lon + integer :: ncid, id_dim, id_var + + real(kind=4), allocatable :: dummy(:,:) + + print*,"- READ MODEL LAND MASK FILE" + + print*,'- OPEN LAND MASK FILE: ', orog_file + error=nf90_open(orog_file,nf90_nowrite,ncid) + call netcdf_err(error, "OPENING MODEL LAND MASK FILE") + + print*,"- READ I-DIMENSION" + error=nf90_inq_dimid(ncid, 'lon', id_dim) + call netcdf_err(error, "READING LON ID") + error=nf90_inquire_dimension(ncid,id_dim,len=lon) + call netcdf_err(error, "READING LON") + + print*,"- READ J-DIMENSION" + error=nf90_inq_dimid(ncid, 'lat', id_dim) + call netcdf_err(error, "READING LAT ID") + error=nf90_inquire_dimension(ncid,id_dim,len=lat) + call netcdf_err(error, "READING LAT") + + print*,"- I/J DIMENSIONS: ", lon, lat + + if ((lon /= idim) .or. (lat /= jdim)) then + call error_handler("MISMATCH IN DIMENSIONS.") + endif + + allocate(dummy(idim,jdim)) + + print*,"- READ LAND MASK" + error=nf90_inq_varid(ncid, 'slmsk', id_var) + call netcdf_err(error, "READING SLMSK ID") + error=nf90_get_var(ncid, id_var, dummy) + call netcdf_err(error, "READING SLMSK") + + error = nf90_close(ncid) + + mask = nint(dummy) + + deallocate (dummy) + + end subroutine get_model_mask + + subroutine model_grid_cleanup + +!----------------------------------------------------------------------- +! subroutine documentation block +! +! Subroutine: model grid cleanup +! prgmmr: gayno org: w/np2 date: 2018 +! +! Abstract: Free up memory associated with this module. +! +! Usage: call model_grid_cleanup +! +!----------------------------------------------------------------------- + + implicit none + + integer :: rc + + print*,"- CALL GridDestroy FOR MODEL GRID." + call ESMF_GridDestroy(grid_mdl,rc=rc) + + print*,"- CALL FieldDestroy FOR MODEL GRID LAND MASK." + call ESMF_FieldDestroy(mask_field_mdl,rc=rc) + + print*,"- CALL FieldDestroy FOR MODEL GRID DATA FIELD." + call ESMF_FieldDestroy(data_field_mdl,rc=rc) + + print*,"- CALL FieldDestroy FOR MODEL GRID VEGETATION TYPE." + call ESMF_FieldDestroy(vegt_field_mdl,rc=rc) + + end subroutine model_grid_cleanup + + end module model_grid diff --git a/sorc/sfc_climo_gen.fd/source_grid.F90 b/sorc/sfc_climo_gen.fd/source_grid.F90 new file mode 100644 index 000000000..d73e850a1 --- /dev/null +++ b/sorc/sfc_climo_gen.fd/source_grid.F90 @@ -0,0 +1,434 @@ + module source_grid + +!-------------------------------------------------------------------------- +! module documentation block +! +! Module: source_grid +! pgrmmr: gayno org: w/np2 date: 2018 +! +! Abstract: Read grid specs, date information and land/sea mask for +! the source data that will be interpolated to the model grid. +! Also, sets up the ESMF grid object for the source grid. +! Source grid is assumed to be global lat/lon. +! +! Usage: use source_grid +! +! Public Subroutines: +! ------------------- +! define_source_grid Defines esmf grid object for source +! grid. Retrieves date and field +! information from source file. +! source_grid_cleanup Free up memory used in this module. +! +! Public variables: +! ----------------- +! +! day_of_rec Day of each time record with +! respect to Jan 1. +! field_names Names of fields to be processed. +! grid_src ESMF grid object for the source grid. +! i/j_src i/j dimensions of the source grid. +! num_fields Number of fields in the file. Some +! files have more than one (ex: +! the 4-component albedo). +! num_records Number of fields times time records. +! num_time_recs Number of time records. +! source Original source of the data. +! +!------------------------------------------------------------------ + + use esmf + use utils + + implicit none + + private + + character(len=50), allocatable, public :: field_names(:) + character(len=75), public :: source + + integer, public :: i_src, j_src, num_records + integer, public :: num_time_recs + integer, public :: num_fields + integer, allocatable, public :: day_of_rec(:) + + type(esmf_grid), public :: grid_src + + public :: define_source_grid + public :: source_grid_cleanup + + contains + + subroutine define_source_grid(localpet, npets, input_file) + +!----------------------------------------------------------------------- +! subroutine documentation block +! +! Subroutine: define source grid +! prgmmr: gayno org: w/np2 date: 2018 +! +! Abstract: Read date information from input source data file. +! Create esmf grid object for the source grid. +! +! Usage: call define_source_grid(localpet, npets, input_file) +! +! input argument list: +! localpet mpi task number +! npets total number mpi tasks +! input_file file containing the source grid data. +! +!----------------------------------------------------------------------- + + use netcdf + + implicit none + + include 'mpif.h' + + character(len=*), intent(in) :: input_file + + integer, intent(in) :: localpet, npets + + character(len=50) :: vname + + integer :: dimid, dims(1), ncid, status + integer :: count, num_vars, n + integer :: rc, varid, i, j, i_src_corner + integer(esmf_kind_i4), pointer :: mask_ptr(:,:) + integer :: clb(2), cub(2) + integer :: clb_corner(2), cub_corner(2) + + real(esmf_kind_r4), allocatable :: mask_global(:,:) + real(esmf_kind_r8), allocatable :: lat_global(:) + real(esmf_kind_r8), allocatable :: lon_global(:) + real(esmf_kind_r8), allocatable :: lat_corner_global(:) + real(esmf_kind_r8), allocatable :: lon_corner_global(:) + real(esmf_kind_r4), pointer :: mask_field_ptr(:,:) + real(esmf_kind_r8), pointer :: lat_ptr(:,:) + real(esmf_kind_r8), pointer :: lon_ptr(:,:) + real(esmf_kind_r8), pointer :: lat_corner_ptr(:,:) + real(esmf_kind_r8), pointer :: lon_corner_ptr(:,:) + real :: lon_extent + + type(esmf_field) :: mask_field + type(esmf_polekind_flag) :: polekindflag(2) + + print*,"- OPEN FILE: ", trim(input_file) + status = nf90_open(input_file, nf90_nowrite, ncid) + call netcdf_err(status, "OPENING INPUT SOURCE FILE") + + status = nf90_get_att(ncid, nf90_global, 'source', source) + if (status /= nf90_noerr) source="unknown" + + if(localpet == 0) print*,'- SOURCE OF DATA IS: ', trim(source) + + status = nf90_inq_dimid(ncid, 'idim', dimid) + call netcdf_err(status, "READING I DIMENSION ID OF SOURCE DATA") + status = nf90_inquire_dimension(ncid, dimid, len=i_src) + call netcdf_err(status, "READING I DIMENSION OF SOURCE DATA") + + status = nf90_inq_dimid(ncid, 'jdim', dimid) + call netcdf_err(status, "READING J DIMENSION ID OF SOURCE DATA") + status = nf90_inquire_dimension(ncid, dimid, len=j_src) + call netcdf_err(status, "READING J DIMENSION OF SOURCE DATA") + + if(localpet == 0) print*,'- I/J DIMENSION OF SOURCE DATA: ', i_src, j_src + + allocate(lat_global(j_src)) + status = nf90_inq_varid(ncid, 'lat', varid) + call netcdf_err(status, "READING SOURCE DATA LATITUDE ID") + status = nf90_get_var(ncid, varid, lat_global) + call netcdf_err(status, "READING SOURCE DATA LATITUDES") + + allocate(lon_global(i_src)) + status = nf90_inq_varid(ncid, 'lon', varid) + call netcdf_err(status, "READING SOURCE DATA LONGITUDE ID") + status = nf90_get_var(ncid, varid, lon_global) + call netcdf_err(status, "READING SOURCE DATA LONGITUDE") + + allocate(lat_corner_global(j_src+1)) + status = nf90_inq_varid(ncid, 'lat_corner', varid) + call netcdf_err(status, "READING SOURCE DATA CORNER LATITUDE ID") + status = nf90_get_var(ncid, varid, lat_corner_global) + call netcdf_err(status, "READING SOURCE DATA CORNER LATITUDE") + +!----------------------------------------------------------------------- +! Dimension of lon_corner depends on whether input data is periodic +! (global) or regional. +!----------------------------------------------------------------------- + + status = nf90_inq_varid(ncid, 'lon_corner', varid) + call netcdf_err(status, "READING SOURCE DATA CORNER LONGITUDE ID") + status = nf90_inquire_variable(ncid, varid, dimids=dims) + call netcdf_err(status, "READING SOURCE DATA CORNER LONGITUDE DIMIDS") + status = nf90_inquire_dimension(ncid, dims(1), len=i_src_corner) + call netcdf_err(status, "READING SOURCE DATA CORNER LONGITUDE DIMS") + allocate(lon_corner_global(i_src_corner)) + status = nf90_get_var(ncid, varid, lon_corner_global) + call netcdf_err(status, "READING SOURCE DATA CORNER LONGITUDE") + + status = nf90_inq_dimid(ncid, 'time', dimid) + call netcdf_err(status, "READING SOURCE DATA TIME ID") + status = nf90_inquire_dimension(ncid, dimid, len=num_time_recs) + call netcdf_err(status, "READING SOURCE DATA NUM TIME PERIODS") + + allocate(day_of_rec(num_time_recs)) + status = nf90_inq_varid(ncid, 'time', varid) + call netcdf_err(status, "READING SOURCE DATA TIME VARID") + status = nf90_get_var(ncid, varid, day_of_rec) + call netcdf_err(status, "READING SOURCE DATA DAY OF RECORD") + + print*,'- SOURCE DATA DAYS OF RECORD: ', day_of_rec + + status = nf90_inquire(ncid, nVariables=num_vars) + call netcdf_err(status, "READING NUMBER OF VARIABLES SOURCE DATA") + +!----------------------------------------------------------------------- +! Assumes files contain records for time, lat, lon, lat_corner, +! lon_corner. So number of fields processed will be the total +! number of variables minus 5. +!----------------------------------------------------------------------- + + num_fields = num_vars - 5 + num_records = num_vars * num_time_recs + + allocate(field_names(num_fields)) + + count = 0 + do n = 1, num_vars + + status = nf90_inquire_variable(ncid, n, name=vname) + call netcdf_err(status, "READING FIELD NAMES") + + if (trim(vname) == 'time') cycle + if (trim(vname) == 'lon') cycle + if (trim(vname) == 'lon_corner') cycle + if (trim(vname) == 'lat') cycle + if (trim(vname) == 'lat_corner') cycle + + count = count + 1 + field_names(count) = vname + + enddo + + if(localpet==0) print*,'- FIELDS TO BE PROCESSED: ', field_names + + if (localpet == 0) then + allocate(mask_global(i_src,j_src)) + status = nf90_inq_varid(ncid, field_names(1), varid) + call netcdf_err(status, "READING FIELD ID") + status = nf90_get_var(ncid, varid, mask_global) + call netcdf_err(status, "READING FIELD") + else + allocate(mask_global(0,0)) + endif + + status = nf90_close(ncid) + +!-------------------------------------------------------------------------- +! Create ESMF grid object for the source data grid. Check if +! data is periodic in the east/west direction. +! +! Note: When using regional data, there is always the chance of +! the target grid being located outside the input grid. ESMF +! support recommends passing back the dstField (esmf_typekind_i4) from +! all calls to FieldRegridStore and checking for the +! ESMF_REGRIDSTATUS_OUTSIDE flag. +!-------------------------------------------------------------------------- + + lon_extent = mod((lon_global(i_src)-lon_global(1))-1+3600,360.)+1.0 + + if (lon_extent < 350.0) then + + print*,"- CALL GridCreateNoPeriDim FOR SOURCE GRID" + grid_src = ESMF_GridCreateNoPeriDim(minIndex=(/1,1/), & + maxIndex=(/i_src,j_src/), & + coordSys=ESMF_COORDSYS_SPH_DEG, & + regDecomp=(/1,npets/), & + name="source_grid", & + indexflag=ESMF_INDEX_GLOBAL, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridCreateNoPeriDim.", rc) + + else + + polekindflag(1:2) = ESMF_POLEKIND_MONOPOLE + + print*,"- CALL GridCreate1PeriDim FOR SOURCE GRID" + grid_src = ESMF_GridCreate1PeriDim(minIndex=(/1,1/), & + maxIndex=(/i_src,j_src/), & + polekindflag=polekindflag, & + periodicDim=1, & + poleDim=2, & + coordSys=ESMF_COORDSYS_SPH_DEG, & + regDecomp=(/1,npets/), & + name="source_grid", & + indexflag=ESMF_INDEX_GLOBAL, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridCreate1PeriDim.", rc) + + endif + + print*,"- CALL FieldCreate FOR SOURCE GRID LANDMASK." + mask_field = ESMF_FieldCreate(grid_src, & + typekind=ESMF_TYPEKIND_R4, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + name="source grid land mask", & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate.", rc) + + print*,"- CALL FieldScatter FOR SOURCE GRID MASK." + call ESMF_FieldScatter(mask_field, mask_global, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter.", rc) + + print*,"- CALL GridAddItem FOR SOURCE GRID MASK." + call ESMF_GridAddItem(grid_src, & + itemflag=ESMF_GRIDITEM_MASK, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridAddItem.", rc) + + print*,"- CALL GridGetItem FOR SOURCE GRID MASK." + nullify(mask_ptr) + call ESMF_GridGetItem(grid_src, & + itemflag=ESMF_GRIDITEM_MASK, & + farrayPtr=mask_ptr, & + totalLBound=clb, & + totalUBound=cub, & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetItem", rc) + + print*,"- CALL FieldGet FOR SOURCE GRID LANDMASK." + nullify(mask_field_ptr) + call ESMF_FieldGet(mask_field, & + farrayPtr=mask_field_ptr, & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet.", rc) + + do j = clb(2), cub(2) + do i = clb(1), cub(1) + if (mask_field_ptr(i,j) < -1.0) then + mask_ptr(i,j) = 0 + else + mask_ptr(i,j) = 1 + endif + enddo + enddo + + deallocate(mask_global) + + print*,"- CALL FieldDestroy FOR SOURCE GRID LAND MASK." + call ESMF_FieldDestroy(mask_field,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldDestroy.", rc) + +! Set lat/lons of grid points + + print*,"- CALL GridAddCoord FOR SOURCE GRID CENTER LOCATION." + call ESMF_GridAddCoord(grid_src, & + staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridAddCoord.", rc) + + print*,"- CALL GridGetCoord FOR SOURCE GRID CENTER LONGITUDE." + nullify(lon_ptr) + call ESMF_GridGetCoord(grid_src, & + staggerLoc=ESMF_STAGGERLOC_CENTER, & + coordDim=1, & + farrayPtr=lon_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord.", rc) + + print*,"- CALL GridGetCoord FOR SOURCE GRID CENTER LATITUDE." + nullify(lat_ptr) + call ESMF_GridGetCoord(grid_src, & + staggerLoc=ESMF_STAGGERLOC_CENTER, & + coordDim=2, & + farrayPtr=lat_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord.", rc) + + do j = clb(2), cub(2) + lat_ptr(:,j) = lat_global(j) + enddo + + do i = clb(1), cub(1) + lon_ptr(i,:) = lon_global(i) + enddo + + print*,"- CALL GridAddCoord FOR SOURCE GRID CORNER LOCATION." + call ESMF_GridAddCoord(grid_src, & + staggerloc=ESMF_STAGGERLOC_CORNER, & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridAddCoord.", rc) + + print*,"- CALL GridGetCoord FOR SOURCE GRID CORNER LONGITUDE." + nullify(lon_corner_ptr) + call ESMF_GridGetCoord(grid_src, & + staggerLoc=ESMF_STAGGERLOC_CORNER, & + coordDim=1, & + farrayPtr=lon_corner_ptr, & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord.", rc) + + print*,"- CALL GridGetCoord FOR SOURCE GRID CORNER LATITUDE." + nullify(lat_corner_ptr) + call ESMF_GridGetCoord(grid_src, & + staggerLoc=ESMF_STAGGERLOC_CORNER, & + coordDim=2, & + computationalLBound=clb_corner, & + computationalUBound=cub_corner, & + farrayPtr=lat_corner_ptr, & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord.", rc) + + do j = clb_corner(2), cub_corner(2) + lat_corner_ptr(:,j) = lat_corner_global(j) + enddo + + do i = clb_corner(1), cub_corner(1) + lon_corner_ptr(i,:) = lon_corner_global(i) + enddo + + deallocate(lat_global) + deallocate(lon_global) + deallocate(lat_corner_global) + deallocate(lon_corner_global) + + end subroutine define_source_grid + + subroutine source_grid_cleanup + +!----------------------------------------------------------------------- +! subroutine documentation block +! +! Subroutine: source grid cleanup +! prgmmr: gayno org: w/np2 date: 2018 +! +! Abstract: Free up memory associated with this module. +! +! Usage: call source_grid_cleanup +! +!----------------------------------------------------------------------- + + implicit none + + integer :: rc + + print*,"- CALL GridDestroy FOR SOURCE GRID." + call ESMF_GridDestroy(grid_src,rc=rc) + + deallocate(day_of_rec) + deallocate(field_names) + + end subroutine source_grid_cleanup + + end module source_grid