diff --git a/CMakeLists.txt b/CMakeLists.txt index e5c756d15..5252802ed 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,9 +7,9 @@ endif (NOT PROJECT) #------------------------------------------------------------------------------ cmake_minimum_required(VERSION 3.0) -project(ccppphys - VERSION 4.0.0 - LANGUAGES C CXX Fortran) +project(ccpp_physics + VERSION 5.0.0 + LANGUAGES Fortran) # Use rpaths on MacOSX set(CMAKE_MACOSX_RPATH 1) @@ -31,25 +31,9 @@ set(AUTHORS "Grant Firl" "Dom Heinzeller" "Man Zhang" "Laurie Carson") if (OPENMP) include(detect_openmp) detect_openmp() - #set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") - #set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") - #set (CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} ${OpenMP_Fortran_FLAGS}") - message(STATUS "Enable OpenMP support for C/C++/Fortran compiler") + message(STATUS "Enable OpenMP support") else (OPENMP) - message (STATUS "Disable OpenMP support for C/C++/Fortran compiler") -endif() - -#------------------------------------------------------------------------------ -# The Fortran compiler/linker flag inserted by cmake to create shared libraries -# with the Intel compiler is deprecated (-i_dynamic), correct here. -# CMAKE_Fortran_COMPILER_ID = {"Intel", "PGI", "GNU", "Clang", "MSVC", ...} -if ("${CMAKE_Fortran_COMPILER_ID}" STREQUAL "Intel") - string(REPLACE "-i_dynamic" "-shared-intel" - CMAKE_SHARED_LIBRARY_CREATE_Fortran_FLAGS - "${CMAKE_SHARED_LIBRARY_CREATE_Fortran_FLAGS}") - string(REPLACE "-i_dynamic" "-shared-intel" - CMAKE_SHARED_LIBRARY_LINK_Fortran_FLAGS - "${CMAKE_SHARED_LIBRARY_LINK_Fortran_FLAGS}") + message (STATUS "Disable OpenMP support") endif() #------------------------------------------------------------------------------ @@ -321,22 +305,8 @@ else() message (FATAL_ERROR "This program has only been compiled with gfortran, pgf90 and ifort. If another compiler is needed, the appropriate flags must be added in ${GFS_PHYS_SRC}/CMakeLists.txt") endif() -# The auto-generated caps can contain calls to physics schemes in -# which some of the arguments (pointers, arrays) are not associated/allocated. -# This is on purpose to avoid allocating fields that are not used inside the -# scheme if, for example, certain conditions are not met. To avoid -# Fortran runtime errors, it is necessary to remove checks for pointers -# that are not associated and for array bounds from the caps ONLY. For the -# physics schemes, these checks can and should remain enabled. Overwriting -# the pointer check flags explicitly works for Intel and GNU, but not for PGI. -if (${CMAKE_Fortran_COMPILER_ID} STREQUAL "GNU") - set_property(SOURCE ${CAPS} APPEND_STRING PROPERTY COMPILE_FLAGS " -fcheck=no-pointer,no-bounds ") -elseif (${CMAKE_Fortran_COMPILER_ID} STREQUAL "Intel") - set_property(SOURCE ${CAPS} APPEND_STRING PROPERTY COMPILE_FLAGS " -check nopointers,nobounds ") -endif () - #------------------------------------------------------------------------------ -add_library(ccppphys STATIC ${SCHEMES} ${SCHEMES_SFX_OPT} ${CAPS}) +add_library(ccpp_physics STATIC ${SCHEMES} ${SCHEMES_SFX_OPT} ${CAPS}) # Generate list of Fortran modules from defined sources foreach(source_f90 ${CAPS}) get_filename_component(tmp_source_f90 ${source_f90} NAME) @@ -345,20 +315,25 @@ foreach(source_f90 ${CAPS}) list(APPEND MODULES_F90 ${CMAKE_CURRENT_BINARY_DIR}/${module_f90}) endforeach() -set_target_properties(ccppphys PROPERTIES VERSION ${PROJECT_VERSION} +set_target_properties(ccpp_physics PROPERTIES VERSION ${PROJECT_VERSION} SOVERSION ${PROJECT_VERSION_MAJOR}) +target_include_directories(ccpp_physics PUBLIC + $) + +target_link_libraries(ccpp_physics PUBLIC w3nco::w3nco_d NetCDF::NetCDF_Fortran) + if (PROJECT STREQUAL "CCPP-FV3") # Define where to install the library - install(TARGETS ccppphys - EXPORT ccppphys-targets + install(TARGETS ccpp_physics + EXPORT ccpp_physics-targets ARCHIVE DESTINATION lib LIBRARY DESTINATION lib RUNTIME DESTINATION lib ) # Export our configuration - install(EXPORT ccppphys-targets - FILE ccppphys-config.cmake + install(EXPORT ccpp_physics-targets + FILE ccpp_physics-config.cmake DESTINATION lib/cmake ) # Define where to install the C headers and Fortran modules diff --git a/CODEOWNERS b/CODEOWNERS index 0d5230f89..d7c3658fd 100644 --- a/CODEOWNERS +++ b/CODEOWNERS @@ -3,7 +3,7 @@ # These owners will be the default owners for everything in the repo. #* @defunkt -* @climbfuji @llpcarson @grantfirl @JulieSchramm +* @climbfuji @llpcarson @grantfirl # Order is important. The last matching pattern has the most precedence. # So if a pull request only touches javascript files, only these owners diff --git a/physics/GFS_DCNV_generic.meta b/physics/GFS_DCNV_generic.meta index f9ea812b5..c64e1fadb 100644 --- a/physics/GFS_DCNV_generic.meta +++ b/physics/GFS_DCNV_generic.meta @@ -267,7 +267,7 @@ standard_name = convective_transportable_tracers long_name = array to contain cloud water and other convective trans. tracers units = kg kg-1 - dimensions = (horizontal_loop_extent,vertical_dimension,number_of_tracers_for_convective_transport) + dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_tracers_for_convective_transport) type = real kind = kind_phys intent = in @@ -477,7 +477,7 @@ standard_name = cumulative_change_of_state_variables long_name = diagnostic tendencies for state variables units = various - dimensions = (horizontal_loop_extent,vertical_dimension,number_of_cumulative_change_processes) + dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_cumulative_change_processes) type = real kind = kind_phys intent = inout @@ -534,7 +534,7 @@ standard_name = tracer_concentration_of_new_state long_name = tracer concentration updated by physics units = kg kg-1 - dimensions = (horizontal_loop_extent,vertical_dimension,number_of_tracers) + dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_tracers) type = real kind = kind_phys intent = in @@ -543,7 +543,7 @@ standard_name = tracer_concentration_save long_name = tracer concentration before entering a physics scheme units = kg kg-1 - dimensions = (horizontal_loop_extent,vertical_dimension,number_of_tracers) + dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_tracers) type = real kind = kind_phys intent = in @@ -778,7 +778,7 @@ standard_name = convective_transportable_tracers long_name = array to contain cloud water and other convective trans. tracers units = kg kg-1 - dimensions = (horizontal_loop_extent,vertical_dimension,number_of_tracers_for_convective_transport) + dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_tracers_for_convective_transport) type = real kind = kind_phys intent = in diff --git a/physics/GFS_MP_generic.meta b/physics/GFS_MP_generic.meta index 30de79201..a87cfe578 100644 --- a/physics/GFS_MP_generic.meta +++ b/physics/GFS_MP_generic.meta @@ -823,7 +823,7 @@ standard_name = cumulative_change_of_state_variables long_name = diagnostic tendencies for state variables units = various - dimensions = (horizontal_loop_extent,vertical_dimension,number_of_cumulative_change_processes) + dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_cumulative_change_processes) type = real kind = kind_phys intent = inout @@ -832,7 +832,7 @@ standard_name = cumulative_change_of_state_variables_outer_index long_name = index of state-variable and process in last dimension of diagnostic tendencies array AKA cumulative_change_index units = index - dimensions = (number_of_tracers_plus_one_hundred,number_of_causes) + dimensions = (number_of_tracers_plus_one_hundred,number_of_cumulative_change_processes) type = integer intent = in optional = F diff --git a/physics/GFS_PBL_generic.meta b/physics/GFS_PBL_generic.meta index a3eb00cd0..84614d8d0 100644 --- a/physics/GFS_PBL_generic.meta +++ b/physics/GFS_PBL_generic.meta @@ -836,7 +836,7 @@ standard_name = cumulative_change_of_state_variables long_name = diagnostic tendencies for state variables units = various - dimensions = (horizontal_loop_extent,vertical_dimension,number_of_cumulative_change_processes) + dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_cumulative_change_processes) type = real kind = kind_phys intent = inout diff --git a/physics/GFS_SCNV_generic.meta b/physics/GFS_SCNV_generic.meta index 5bed0b907..adf82ed92 100644 --- a/physics/GFS_SCNV_generic.meta +++ b/physics/GFS_SCNV_generic.meta @@ -259,7 +259,7 @@ standard_name = convective_transportable_tracers long_name = array to contain cloud water and other convective trans. tracers units = kg kg-1 - dimensions = (horizontal_loop_extent,vertical_dimension,number_of_tracers_for_convective_transport) + dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_tracers_for_convective_transport) type = real kind = kind_phys intent = in diff --git a/physics/GFS_debug.F90 b/physics/GFS_debug.F90 index 43f48dfed..deb88458b 100644 --- a/physics/GFS_debug.F90 +++ b/physics/GFS_debug.F90 @@ -455,55 +455,58 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, do impi=0,mpisize-1 do iomp=0,ompsize-1 if (mpirank==impi .and. omprank==iomp) then - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Model%kdt' , Model%kdt) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Model%kdt' , Model%kdt) ! Sfcprop - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%slmsk' , Sfcprop%slmsk) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%oceanfrac', Sfcprop%oceanfrac) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%landfrac' , Sfcprop%landfrac) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%lakefrac' , Sfcprop%lakefrac) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%tsfc' , Sfcprop%tsfc) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%tsfco' , Sfcprop%tsfco) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%tsfcl' , Sfcprop%tsfcl) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%tisfc' , Sfcprop%tisfc) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%snowd' , Sfcprop%snowd) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%zorl' , Sfcprop%zorl) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%zorlw' , Sfcprop%zorlw) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%zorll' , Sfcprop%zorll) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%zorli' , Sfcprop%zorli) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%zorlwav' , Sfcprop%zorlwav) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%fice' , Sfcprop%fice) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%hprime' , Sfcprop%hprime) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%sncovr' , Sfcprop%sncovr) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%snoalb' , Sfcprop%snoalb) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%alvsf' , Sfcprop%alvsf) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%alnsf' , Sfcprop%alnsf) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%alvwf' , Sfcprop%alvwf) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%alnwf' , Sfcprop%alnwf) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%facsf' , Sfcprop%facsf) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%facwf' , Sfcprop%facwf) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%slope' , Sfcprop%slope) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%shdmin' , Sfcprop%shdmin) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%shdmax' , Sfcprop%shdmax) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%tg3' , Sfcprop%tg3) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%vfrac' , Sfcprop%vfrac) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%vtype' , Sfcprop%vtype) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%stype' , Sfcprop%stype) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%uustar' , Sfcprop%uustar) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%oro' , Sfcprop%oro) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%oro_uf' , Sfcprop%oro_uf) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%hice' , Sfcprop%hice) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%weasd' , Sfcprop%weasd) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%canopy' , Sfcprop%canopy) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%ffmm' , Sfcprop%ffmm) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%ffhh' , Sfcprop%ffhh) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%f10m' , Sfcprop%f10m) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%tprcp' , Sfcprop%tprcp) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%srflag' , Sfcprop%srflag) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%slc' , Sfcprop%slc) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%smc' , Sfcprop%smc) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%stc' , Sfcprop%stc) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%t2m' , Sfcprop%t2m) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%q2m' , Sfcprop%q2m) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%slmsk' , Sfcprop%slmsk) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%oceanfrac' , Sfcprop%oceanfrac) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%landfrac' , Sfcprop%landfrac) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%lakefrac' , Sfcprop%lakefrac) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%tsfc' , Sfcprop%tsfc) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%tsfco' , Sfcprop%tsfco) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%tsfcl' , Sfcprop%tsfcl) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%tisfc' , Sfcprop%tisfc) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%snowd' , Sfcprop%snowd) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%zorl' , Sfcprop%zorl) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%zorlw' , Sfcprop%zorlw) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%zorll' , Sfcprop%zorll) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%zorli' , Sfcprop%zorli) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%zorlwav' , Sfcprop%zorlwav) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%fice' , Sfcprop%fice) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%hprime' , Sfcprop%hprime) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%sncovr' , Sfcprop%sncovr) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%snoalb' , Sfcprop%snoalb) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%alvsf' , Sfcprop%alvsf) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%alnsf' , Sfcprop%alnsf) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%alvwf' , Sfcprop%alvwf) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%alnwf' , Sfcprop%alnwf) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%facsf' , Sfcprop%facsf) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%facwf' , Sfcprop%facwf) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%slope' , Sfcprop%slope) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%slope_save', Sfcprop%slope_save) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%shdmin' , Sfcprop%shdmin) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%shdmax' , Sfcprop%shdmax) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%tg3' , Sfcprop%tg3) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%vfrac' , Sfcprop%vfrac) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%vtype' , Sfcprop%vtype) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%vtype_save', Sfcprop%vtype_save) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%stype' , Sfcprop%stype) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%stype_save', Sfcprop%stype_save) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%uustar' , Sfcprop%uustar) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%oro' , Sfcprop%oro) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%oro_uf' , Sfcprop%oro_uf) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%hice' , Sfcprop%hice) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%weasd' , Sfcprop%weasd) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%canopy' , Sfcprop%canopy) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%ffmm' , Sfcprop%ffmm) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%ffhh' , Sfcprop%ffhh) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%f10m' , Sfcprop%f10m) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%tprcp' , Sfcprop%tprcp) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%srflag' , Sfcprop%srflag) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%slc' , Sfcprop%slc) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%smc' , Sfcprop%smc) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%stc' , Sfcprop%stc) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%t2m' , Sfcprop%t2m) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%q2m' , Sfcprop%q2m) if (Model%nstf_name(1)>0) then call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%tref ', Sfcprop%tref) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%z_c ', Sfcprop%z_c) @@ -544,18 +547,18 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%emis_lnd', Sfcprop%emis_lnd) ! NoahMP and RUC if (Model%lsm == Model%lsm_ruc .or. Model%lsm == Model%lsm_noahmp) then - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%albdvis_lnd', Sfcprop%albdvis_lnd) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%albdnir_lnd', Sfcprop%albdnir_lnd) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%albivis_lnd', Sfcprop%albivis_lnd) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%albinir_lnd', Sfcprop%albinir_lnd) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%albdirvis_lnd', Sfcprop%albdirvis_lnd) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%albdirnir_lnd', Sfcprop%albdirnir_lnd) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%albdifvis_lnd', Sfcprop%albdifvis_lnd) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%albdifnir_lnd', Sfcprop%albdifnir_lnd) end if ! RUC only if (Model%lsm == Model%lsm_ruc) then call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%emis_ice', Sfcprop%emis_ice) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%albdvis_ice', Sfcprop%albdvis_ice) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%albdnir_ice', Sfcprop%albdnir_ice) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%albivis_ice', Sfcprop%albivis_ice) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%albinir_ice', Sfcprop%albinir_ice) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%albdirvis_ice', Sfcprop%albdirvis_ice) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%albdifvis_ice', Sfcprop%albdifvis_ice) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%albdirnir_ice', Sfcprop%albdirnir_ice) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%albdifnir_ice', Sfcprop%albdifnir_ice) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%sfalb_lnd', Sfcprop%sfalb_lnd) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%sfalb_ice', Sfcprop%sfalb_ice) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Sfcprop%sfalb_lnd_bck', Sfcprop%sfalb_lnd_bck) @@ -691,10 +694,10 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, else call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%wet1 ', Diag%wet1) end if - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%skebu_wts ', Diag%skebu_wts) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%skebv_wts ', Diag%skebv_wts) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%sppt_wts ', Diag%sppt_wts) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%shum_wts ', Diag%shum_wts) + !call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%skebu_wts ', Diag%skebu_wts) + !call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%skebv_wts ', Diag%skebv_wts) + !call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%sppt_wts ', Diag%sppt_wts) + !call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%shum_wts ', Diag%shum_wts) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%zmtnblck ', Diag%zmtnblck) if (Model%ldiag3d) then !do itracer=2,Model%ntracp100 @@ -1199,6 +1202,7 @@ subroutine GFS_interstitialtoscreen_run (Model, Statein, Stateout, Sfcprop, Coup call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%evap_water ', Interstitial%evap_water ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%evbs ', Interstitial%evbs ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%evcw ', Interstitial%evcw ) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%ext_diag_thompson_reset', Interstitial%ext_diag_thompson_reset) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%faerlw ', Interstitial%faerlw ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%faersw ', Interstitial%faersw ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%ffhh_ice ', Interstitial%ffhh_ice ) @@ -1262,6 +1266,7 @@ subroutine GFS_interstitialtoscreen_run (Model, Statein, Stateout, Sfcprop, Coup call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%kpbl ', Interstitial%kpbl ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%kt ', Interstitial%kt ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%ktop ', Interstitial%ktop ) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%max_hourly_reset ', Interstitial%max_hourly_reset ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%mbota ', Interstitial%mbota ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%mtopa ', Interstitial%mtopa ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%nday ', Interstitial%nday ) @@ -1286,7 +1291,6 @@ subroutine GFS_interstitialtoscreen_run (Model, Statein, Stateout, Sfcprop, Coup call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%rb_ice ', Interstitial%rb_ice ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%rb_land ', Interstitial%rb_land ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%rb_water ', Interstitial%rb_water ) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%reset ', Interstitial%reset ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%rhc ', Interstitial%rhc ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%runoff ', Interstitial%runoff ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%save_q ', Interstitial%save_q ) @@ -1309,20 +1313,17 @@ subroutine GFS_interstitialtoscreen_run (Model, Statein, Stateout, Sfcprop, Coup call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%sigmaf ', Interstitial%sigmaf ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%sigmafrac ', Interstitial%sigmafrac ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%sigmatot ', Interstitial%sigmatot ) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%slopetype ', Interstitial%slopetype ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%snowc ', Interstitial%snowc ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%snowd_ice ', Interstitial%snowd_ice ) ! call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%snowd_land ', Interstitial%snowd_land ) ! call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%snowd_water ', Interstitial%snowd_water ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%snohf ', Interstitial%snohf ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%snowmt ', Interstitial%snowmt ) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%soiltype ', Interstitial%soiltype ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%stress ', Interstitial%stress ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%stress_ice ', Interstitial%stress_ice ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%stress_land ', Interstitial%stress_land ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%stress_water ', Interstitial%stress_water ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%theta ', Interstitial%theta ) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%tice ', Interstitial%tice ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%tlvl ', Interstitial%tlvl ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%tlyr ', Interstitial%tlyr ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%tprcp_ice ', Interstitial%tprcp_ice ) @@ -1332,7 +1333,6 @@ subroutine GFS_interstitialtoscreen_run (Model, Statein, Stateout, Sfcprop, Coup call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%tseal ', Interstitial%tseal ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%tsfa ', Interstitial%tsfa ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%tsfc_ice ', Interstitial%tsfc_ice ) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%tsfc_land ', Interstitial%tsfc_land ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%tsfc_water ', Interstitial%tsfc_water ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%tsfg ', Interstitial%tsfg ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%tsurf_ice ', Interstitial%tsurf_ice ) @@ -1344,7 +1344,6 @@ subroutine GFS_interstitialtoscreen_run (Model, Statein, Stateout, Sfcprop, Coup call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%uustar_water ', Interstitial%uustar_water ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%vdftra ', Interstitial%vdftra ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%vegf1d ', Interstitial%vegf1d ) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%vegtype ', Interstitial%vegtype ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%wcbmax ', Interstitial%wcbmax ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%weasd_ice ', Interstitial%weasd_ice ) ! call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%weasd_land ', Interstitial%weasd_land ) @@ -1551,8 +1550,7 @@ end subroutine GFS_checkland_finalize !! subroutine GFS_checkland_run (me, master, blkno, im, kdt, iter, flag_iter, flag_guess, & flag_init, flag_restart, frac_grid, isot, ivegsrc, stype, vtype, slope, & - soiltyp, vegtype, slopetyp, dry, icy, wet, lake, ocean, & - oceanfrac, landfrac, lakefrac, slmsk, islmsk, & + dry, icy, wet, lake, ocean, oceanfrac, landfrac, lakefrac, slmsk, islmsk, & zorl, zorlw, zorll, zorli, fice, errmsg, errflg ) use machine, only: kind_phys @@ -1566,34 +1564,31 @@ subroutine GFS_checkland_run (me, master, blkno, im, kdt, iter, flag_iter, flag_ integer, intent(in ) :: im integer, intent(in ) :: kdt integer, intent(in ) :: iter - logical, intent(in ) :: flag_iter(im) - logical, intent(in ) :: flag_guess(im) + logical, intent(in ) :: flag_iter(:) + logical, intent(in ) :: flag_guess(:) logical, intent(in ) :: flag_init logical, intent(in ) :: flag_restart logical, intent(in ) :: frac_grid integer, intent(in ) :: isot integer, intent(in ) :: ivegsrc - real(kind_phys), intent(in ) :: stype(im) - real(kind_phys), intent(in ) :: vtype(im) - real(kind_phys), intent(in ) :: slope(im) - integer, intent(in ) :: soiltyp(im) - integer, intent(in ) :: vegtype(im) - integer, intent(in ) :: slopetyp(im) - logical, intent(in ) :: dry(im) - logical, intent(in ) :: icy(im) - logical, intent(in ) :: wet(im) - logical, intent(in ) :: lake(im) - logical, intent(in ) :: ocean(im) - real(kind_phys), intent(in ) :: oceanfrac(im) - real(kind_phys), intent(in ) :: landfrac(im) - real(kind_phys), intent(in ) :: lakefrac(im) - real(kind_phys), intent(in ) :: slmsk(im) - integer, intent(in ) :: islmsk(im) - real(kind_phys), intent(in ) :: zorl(im) - real(kind_phys), intent(in ) :: zorlw(im) - real(kind_phys), intent(in ) :: zorll(im) - real(kind_phys), intent(in ) :: zorli(im) - real(kind_phys), intent(in ) :: fice(im) + integer, intent(in ) :: stype(:) + integer, intent(in ) :: vtype(:) + integer, intent(in ) :: slope(:) + logical, intent(in ) :: dry(:) + logical, intent(in ) :: icy(:) + logical, intent(in ) :: wet(:) + logical, intent(in ) :: lake(:) + logical, intent(in ) :: ocean(:) + real(kind_phys), intent(in ) :: oceanfrac(:) + real(kind_phys), intent(in ) :: landfrac(:) + real(kind_phys), intent(in ) :: lakefrac(:) + real(kind_phys), intent(in ) :: slmsk(:) + integer, intent(in ) :: islmsk(:) + real(kind_phys), intent(in ) :: zorl(:) + real(kind_phys), intent(in ) :: zorlw(:) + real(kind_phys), intent(in ) :: zorll(:) + real(kind_phys), intent(in ) :: zorli(:) + real(kind_phys), intent(in ) :: fice(:) character(len=*), intent( out) :: errmsg integer, intent( out) :: errflg @@ -1623,9 +1618,6 @@ subroutine GFS_checkland_run (me, master, blkno, im, kdt, iter, flag_iter, flag_ write(0,'(a,2i5,1x,e16.7)')'YYY: i, blk, stype(i) :', i, blkno, stype(i) write(0,'(a,2i5,1x,e16.7)')'YYY: i, blk, vtype(i) :', i, blkno, vtype(i) write(0,'(a,2i5,1x,e16.7)')'YYY: i, blk, slope(i) :', i, blkno, slope(i) - write(0,'(a,2i5,1x,i5)') 'YYY: i, blk, soiltyp(i) :', i, blkno, soiltyp(i) - write(0,'(a,2i5,1x,i5)') 'YYY: i, blk, vegtype(i) :', i, blkno, vegtype(i) - write(0,'(a,2i5,1x,i5)') 'YYY: i, blk, slopetyp(i) :', i, blkno, slopetyp(i) write(0,'(a,2i5,1x,1x,l)') 'YYY: i, blk, dry(i) :', i, blkno, dry(i) write(0,'(a,2i5,1x,1x,l)') 'YYY: i, blk, icy(i) :', i, blkno, icy(i) write(0,'(a,2i5,1x,1x,l)') 'YYY: i, blk, wet(i) :', i, blkno, wet(i) diff --git a/physics/GFS_debug.meta b/physics/GFS_debug.meta index 4615e163a..fb77772eb 100644 --- a/physics/GFS_debug.meta +++ b/physics/GFS_debug.meta @@ -604,51 +604,24 @@ intent = in optional = F [stype] - standard_name = soil_type_classification_real - long_name = soil type for lsm - units = index - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = in - optional = F -[vtype] - standard_name = vegetation_type_classification_real - long_name = vegetation type for lsm - units = index - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = in - optional = F -[slope] - standard_name = surface_slope_classification_real - long_name = sfc slope type for lsm - units = index - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = in - optional = F -[soiltyp] standard_name = soil_type_classification - long_name = soil type at each grid cell + long_name = soil type for lsm units = index dimensions = (horizontal_loop_extent) type = integer intent = in optional = F -[vegtype] +[vtype] standard_name = vegetation_type_classification - long_name = vegetation type at each grid cell + long_name = vegetation type for lsm units = index dimensions = (horizontal_loop_extent) type = integer intent = in optional = F -[slopetyp] +[slope] standard_name = surface_slope_classification - long_name = surface slope type at each grid cell + long_name = sfc slope type for lsm units = index dimensions = (horizontal_loop_extent) type = integer diff --git a/physics/GFS_phys_time_vary.fv3.F90 b/physics/GFS_phys_time_vary.fv3.F90 index 20c6c68c3..a8ecc1a5e 100644 --- a/physics/GFS_phys_time_vary.fv3.F90 +++ b/physics/GFS_phys_time_vary.fv3.F90 @@ -105,9 +105,9 @@ subroutine GFS_phys_time_vary_init ( integer, intent(in) :: isot, ivegsrc, nlunit real(kind_phys), intent(inout) :: sncovr(:), sncovr_ice(:) - integer, intent(in) :: lsm, lsm_noahmp, lsm_ruc + integer, intent(in) :: lsm, lsm_noahmp, lsm_ruc, vtype(:) real(kind_phys), intent(in) :: min_seaice, fice(:) - real(kind_phys), intent(in) :: landfrac(:), vtype(:) + real(kind_phys), intent(in) :: landfrac(:) real(kind_phys), intent(inout) :: weasd(:) ! NoahMP - only allocated when NoahMP is used @@ -165,7 +165,7 @@ subroutine GFS_phys_time_vary_init ( real(kind_phys), intent(in) :: snowd(:) real(kind_phys), intent(in) :: canopy(:) real(kind_phys), intent(in) :: tg3(:) - real(kind_phys), intent(in) :: stype(:) + integer, intent(in) :: stype(:) real(kind_phys), intent(in) :: con_t0c integer, intent(in) :: nthrds @@ -765,9 +765,10 @@ subroutine GFS_phys_time_vary_timestep_init ( tslb(:,:), tiice(:,:), tg3(:), tref(:), & tsfc(:), tsfco(:), tisfc(:), hice(:), fice(:), & facsf(:), facwf(:), alvsf(:), alvwf(:), alnsf(:), alnwf(:), & - zorli(:), zorll(:), zorlo(:), weasd(:), slope(:), snoalb(:), & - canopy(:), vfrac(:), vtype(:), stype(:), shdmin(:), shdmax(:), & + zorli(:), zorll(:), zorlo(:), weasd(:), snoalb(:), & + canopy(:), vfrac(:), shdmin(:), shdmax(:), & snowd(:), cv(:), cvb(:), cvt(:), oro(:), oro_uf(:), slmsk(:) + integer, intent(inout) :: vtype(:), stype(:), slope(:) character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg diff --git a/physics/GFS_phys_time_vary.fv3.meta b/physics/GFS_phys_time_vary.fv3.meta index 3fb2473bd..979200a85 100644 --- a/physics/GFS_phys_time_vary.fv3.meta +++ b/physics/GFS_phys_time_vary.fv3.meta @@ -444,12 +444,11 @@ intent = in optional = F [vtype] - standard_name = vegetation_type_classification_real + standard_name = vegetation_type_classification long_name = vegetation type for lsm units = index dimensions = (horizontal_dimension) - type = real - kind = kind_phys + type = integer intent = in optional = F [weasd] @@ -963,12 +962,11 @@ intent = in optional = F [stype] - standard_name = soil_type_classification_real + standard_name = soil_type_classification long_name = soil type for lsm units = index dimensions = (horizontal_dimension) - type = real - kind = kind_phys + type = integer intent = in optional = F [con_t0c] @@ -1863,12 +1861,11 @@ intent = inout optional = F [slope] - standard_name = surface_slope_classification_real + standard_name = surface_slope_classification long_name = sfc slope type for lsm units = index dimensions = (horizontal_dimension) - type = real - kind = kind_phys + type = integer intent = inout optional = F [snoalb] @@ -1899,21 +1896,19 @@ intent = inout optional = F [vtype] - standard_name = vegetation_type_classification_real + standard_name = vegetation_type_classification long_name = vegetation type for lsm units = index dimensions = (horizontal_dimension) - type = real - kind = kind_phys + type = integer intent = inout optional = F [stype] - standard_name = soil_type_classification_real + standard_name = soil_type_classification long_name = soil type for lsm units = index dimensions = (horizontal_dimension) - type = real - kind = kind_phys + type = integer intent = inout optional = F [shdmin] diff --git a/physics/GFS_phys_time_vary.scm.F90 b/physics/GFS_phys_time_vary.scm.F90 index e0f380276..b06e46cdc 100644 --- a/physics/GFS_phys_time_vary.scm.F90 +++ b/physics/GFS_phys_time_vary.scm.F90 @@ -101,7 +101,8 @@ subroutine GFS_phys_time_vary_init ( real(kind_phys), intent(inout) :: sncovr(:), sncovr_ice(:) integer, intent(in) :: lsm, lsm_noahmp, lsm_ruc real(kind_phys), intent(in) :: min_seaice, fice(:) - real(kind_phys), intent(in) :: landfrac(:), vtype(:) + real(kind_phys), intent(in) :: landfrac(:) + integer, intent(in) :: vtype(:) real(kind_phys), intent(inout) :: weasd(:) ! NoahMP - only allocated when NoahMP is used @@ -159,7 +160,7 @@ subroutine GFS_phys_time_vary_init ( real(kind_phys), intent(in) :: snowd(:) real(kind_phys), intent(in) :: canopy(:) real(kind_phys), intent(in) :: tg3(:) - real(kind_phys), intent(in) :: stype(:) + integer, intent(in) :: stype(:) real(kind_phys), intent(in) :: con_t0c integer, intent(in) :: nthrds diff --git a/physics/GFS_phys_time_vary.scm.meta b/physics/GFS_phys_time_vary.scm.meta index 1dd03d6b7..a075e8d82 100644 --- a/physics/GFS_phys_time_vary.scm.meta +++ b/physics/GFS_phys_time_vary.scm.meta @@ -444,12 +444,11 @@ intent = in optional = F [vtype] - standard_name = vegetation_type_classification_real + standard_name = vegetation_type_classification long_name = vegetation type for lsm units = index dimensions = (horizontal_dimension) - type = real - kind = kind_phys + type = integer intent = in optional = F [weasd] @@ -963,12 +962,11 @@ intent = in optional = F [stype] - standard_name = soil_type_classification_real + standard_name = soil_type_classification long_name = soil type for lsm units = index dimensions = (horizontal_dimension) - type = real - kind = kind_phys + type = integer intent = in optional = F [con_t0c] diff --git a/physics/GFS_radiation_surface.F90 b/physics/GFS_radiation_surface.F90 index 11703c23c..f4b4b1e62 100644 --- a/physics/GFS_radiation_surface.F90 +++ b/physics/GFS_radiation_surface.F90 @@ -57,7 +57,7 @@ end subroutine GFS_radiation_surface_init !! subroutine GFS_radiation_surface_run ( & im, frac_grid, lslwr, lsswr, lsm, lsm_noahmp, lsm_ruc, & - vtype, xlat, xlon, slmsk, lndp_type, n_var_lndp, sfc_alb_pert, & + xlat, xlon, slmsk, lndp_type, n_var_lndp, sfc_alb_pert, & lndp_var_list, lndp_prt_list, landfrac, snowd, sncovr, & sncovr_ice, fice, zorl, hprime, tsfg, tsfa, tisfc, coszen, & min_seaice, min_lakeice, lakefrac, & @@ -78,9 +78,9 @@ subroutine GFS_radiation_surface_run ( & integer, intent(in) :: lsm, lsm_noahmp, lsm_ruc, lndp_type, n_var_lndp real(kind=kind_phys), intent(in) :: min_seaice, min_lakeice - real(kind=kind_phys), dimension(:), intent(in) :: xlat, xlon, vtype, slmsk, & + real(kind=kind_phys), dimension(:), intent(in) :: xlat, xlon, slmsk, & sfc_alb_pert, lndp_prt_list, & - landfrac, lakefrac, & + landfrac, lakefrac, & snowd, sncovr, & sncovr_ice, fice, zorl, & hprime, tsfg, tsfa, tisfc, & @@ -101,8 +101,7 @@ subroutine GFS_radiation_surface_run ( & ! Local variables integer :: i real(kind=kind_phys) :: lndp_alb - real(kind=kind_phys) :: cimin - real(kind=kind_phys), dimension(im) :: fracl, fraci, fraco + real(kind=kind_phys), dimension(im) :: cimin, fracl, fraci, fraco logical, dimension(im) :: icy ! Initialize CCPP error handling variables @@ -114,9 +113,9 @@ subroutine GFS_radiation_surface_run ( & do i=1,im if (lakefrac(i) > f_zero) then - cimin = min_lakeice + cimin(i) = min_lakeice else - cimin = min_seaice + cimin(i) = min_seaice endif enddo @@ -131,7 +130,7 @@ subroutine GFS_radiation_surface_run ( & else fracl(i) = f_zero fraco(i) = f_one - if(fice(i) < cimin) then + if(fice(i) < cimin(i)) then fraci(i) = f_zero icy(i) = .false. else @@ -145,7 +144,7 @@ subroutine GFS_radiation_surface_run ( & do i=1,im fracl(i) = landfrac(i) fraco(i) = max(f_zero, f_one - fracl(i)) - if(fice(i) < cimin) then + if(fice(i) < cimin(i)) then fraci(i) = f_zero icy(i) = .false. else @@ -159,7 +158,7 @@ subroutine GFS_radiation_surface_run ( & if (lslwr) then !> - Call module_radiation_surface::setemis(),to set up surface !! emissivity for LW radiation. - call setemis (lsm, lsm_noahmp, lsm_ruc, vtype, & + call setemis (lsm, lsm_noahmp, lsm_ruc, & frac_grid, xlon, xlat, slmsk, & ! frac_grid, min_seaice, xlon, xlat, slmsk, & snowd, sncovr, sncovr_ice, zorl, tsfg, tsfa, & diff --git a/physics/GFS_radiation_surface.meta b/physics/GFS_radiation_surface.meta index f021cfe4d..6b8fb1e18 100644 --- a/physics/GFS_radiation_surface.meta +++ b/physics/GFS_radiation_surface.meta @@ -118,15 +118,6 @@ type = integer intent = in optional = F -[vtype] - standard_name = vegetation_type_classification_real - long_name = vegetation type for lsm - units = index - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = in - optional = F [xlat] standard_name = latitude long_name = latitude diff --git a/physics/GFS_rrtmg_pre.meta b/physics/GFS_rrtmg_pre.meta index 76bed26d3..5aba39d98 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -612,7 +612,7 @@ standard_name = surface_stochastic_weights_from_coupled_process long_name = weights for stochastic surface physics perturbation units = none - dimensions = (horizontal_loop_extent,number_of_surface_perturbations) + dimensions = (horizontal_loop_extent,number_of_perturbed_land_surface_variables) type = real kind = kind_phys intent = in diff --git a/physics/GFS_rrtmgp_pre.meta b/physics/GFS_rrtmgp_pre.meta index ddda82bb6..e33663748 100644 --- a/physics/GFS_rrtmgp_pre.meta +++ b/physics/GFS_rrtmgp_pre.meta @@ -2,7 +2,7 @@ name = GFS_rrtmgp_pre type = scheme dependencies = funcphys.f90,iounitdef.f,machine.F,module_bfmicrophysics.f,physcons.F90,physparam.f,radcons.f90,radiation_aerosols.f - dependencies = radiation_astronomy.f,radiation_clouds.f,module_mp_thompson.F90,radiation_gases.f,radiation_surface.f,radiation_tools.F90,rrtmg_lw_cloud_optics.F90 + dependencies = radiation_astronomy.f,radiation_clouds.f,radiation_gases.f,radiation_tools.F90,rrtmg_lw_cloud_optics.F90 ######################################################################## [ccpp-arg-table] diff --git a/physics/GFS_stochastics.F90 b/physics/GFS_stochastics.F90 index 0f1124955..30ca67cd5 100644 --- a/physics/GFS_stochastics.F90 +++ b/physics/GFS_stochastics.F90 @@ -7,7 +7,43 @@ module GFS_stochastics contains - subroutine GFS_stochastics_init () +!> \section arg_table_GFS_stochastics_init Argument Table +!! \htmlinclude GFS_stochastics_init.html +!! +!>\section gfs_stochy_general GFS_stochastics_init General Algorithm +!! This is the GFS stochastic physics initialization. +!! -# define vertical tapering for CA global + subroutine GFS_stochastics_init (si,vfact_ca,km,do_ca,ca_global, errmsg, errflg) + + use machine, only: kind_phys + + implicit none + real(kind_phys), dimension(:), intent(in) :: si + real(kind_phys), dimension(:), intent(inout) :: vfact_ca + integer, intent(in) :: km + logical, intent(in) :: do_ca + logical, intent(in) :: ca_global + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + integer :: k,nz + + errmsg = '' + errflg = 0 + if (do_ca .and. ca_global) then + nz=min(km,size(vfact_ca)) + vfact_ca(:)=0.0 + do k=1,nz + if (si(k) .lt. 0.1 .and. si(k) .gt. 0.025) then + vfact_ca(k) = (si(k)-0.025)/(0.1-0.025) + else if (si(k) .lt. 0.025) then + vfact_ca(k) = 0.0 + else + vfact_ca(k) = 1.0 + endif + enddo + vfact_ca(2)=vfact_ca(3)*0.5 + vfact_ca(1)=0.0 + endif end subroutine GFS_stochastics_init subroutine GFS_stochastics_finalize() @@ -26,10 +62,9 @@ end subroutine GFS_stochastics_finalize !! -# interpolates coefficients for prognostic ozone calculation !! -# performs surface data cycling via the GFS gcycle routine subroutine GFS_stochastics_run (im, km, kdt, delt, do_sppt, pert_mp, use_zmtnblck, & - do_shum ,do_skeb, do_ca,ca_global,ca1,si,vfact_ca, & + do_shum ,do_skeb, do_ca,ca_global,ca1,vfact_ca, & zmtnblck, sppt_wts, skebu_wts, skebv_wts, shum_wts,& - sppt_wts_inv, skebu_wts_inv, skebv_wts_inv, & - shum_wts_inv, diss_est, ugrs, vgrs, tgrs, qgrs_wv, & + diss_est, ugrs, vgrs, tgrs, qgrs_wv, & qgrs_cw, qgrs_rw, qgrs_sw, qgrs_iw, qgrs_gl, & gu0, gv0, gt0, gq0_wv, dtdtnp, & gq0_cw, gq0_rw, gq0_sw, gq0_iw, gq0_gl, & @@ -62,11 +97,6 @@ subroutine GFS_stochastics_run (im, km, kdt, delt, do_sppt, pert_mp, use_zmtnblc real(kind_phys), dimension(:,:), intent(in) :: skebv_wts ! shum_wts only allocated if do_shum == .true. real(kind_phys), dimension(:,:), intent(in) :: shum_wts - ! inverse/flipped weights are always allocated - real(kind_phys), dimension(:,:), intent(inout) :: sppt_wts_inv - real(kind_phys), dimension(:,:), intent(inout) :: skebu_wts_inv - real(kind_phys), dimension(:,:), intent(inout) :: skebv_wts_inv - real(kind_phys), dimension(:,:), intent(inout) :: shum_wts_inv real(kind_phys), dimension(:,:), intent(in) :: diss_est real(kind_phys), dimension(:,:), intent(in) :: ugrs real(kind_phys), dimension(:,:), intent(in) :: vgrs @@ -106,8 +136,7 @@ subroutine GFS_stochastics_run (im, km, kdt, delt, do_sppt, pert_mp, use_zmtnblc ! drain_cpl, dsnow_cpl only allocated if cplflx == .true. or cplchm == .true. real(kind_phys), dimension(:), intent(in) :: drain_cpl real(kind_phys), dimension(:), intent(in) :: dsnow_cpl - real(kind_phys), dimension(:), intent(in) :: si - real(kind_phys), dimension(:), intent(inout) :: vfact_ca + real(kind_phys), dimension(:), intent(in) :: vfact_ca real(kind_phys), dimension(:), intent(in) :: ca1 character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -144,7 +173,6 @@ subroutine GFS_stochastics_run (im, km, kdt, delt, do_sppt, pert_mp, use_zmtnblc if (use_zmtnblck)then sppt_wts(i,k)=(sppt_wts(i,k)-1)*sppt_vwt+1.0 endif - sppt_wts_inv(i,k)=sppt_wts(i,k) upert = (gu0(i,k) - ugrs(i,k)) * sppt_wts(i,k) vpert = (gv0(i,k) - vgrs(i,k)) * sppt_wts(i,k) @@ -225,19 +253,8 @@ subroutine GFS_stochastics_run (im, km, kdt, delt, do_sppt, pert_mp, use_zmtnblc if (do_ca .and. ca_global) then - if(kdt == 1)then - do k=1,km - if (si(k) .lt. 0.1 .and. si(k) .gt. 0.025) then - vfact_ca(k) = (si(k)-0.025)/(0.1-0.025) - else if (si(k) .lt. 0.025) then - vfact_ca(k) = 0.0 - else - vfact_ca(k) = 1.0 - endif - enddo - vfact_ca(2)=vfact_ca(3)*0.5 - vfact_ca(1)=0.0 - endif + !if(kdt == 1)then + !endif do k = 1,km do i = 1,im @@ -340,7 +357,6 @@ subroutine GFS_stochastics_run (im, km, kdt, delt, do_sppt, pert_mp, use_zmtnblc if (do_shum) then do k=1,km gq0_wv(:,k) = gq0_wv(:,k)*(1.0 + shum_wts(:,k)) - shum_wts_inv(:,k) = shum_wts(:,k) end do endif @@ -348,8 +364,6 @@ subroutine GFS_stochastics_run (im, km, kdt, delt, do_sppt, pert_mp, use_zmtnblc do k=1,km gu0(:,k) = gu0(:,k)+skebu_wts(:,k)*(diss_est(:,k)) gv0(:,k) = gv0(:,k)+skebv_wts(:,k)*(diss_est(:,k)) - skebu_wts_inv(:,k) = skebu_wts(:,k) - skebv_wts_inv(:,k) = skebv_wts(:,k) enddo endif diff --git a/physics/GFS_stochastics.meta b/physics/GFS_stochastics.meta index dd7bda9c4..b5574d569 100644 --- a/physics/GFS_stochastics.meta +++ b/physics/GFS_stochastics.meta @@ -3,6 +3,68 @@ type = scheme dependencies = machine.F +[ccpp-arg-table] + name = GFS_stochastics_init + type = scheme +[km] + standard_name = vertical_dimension_for_radiation + long_name = number of vertical levels for radiation calculations + units = count + dimensions = () + type = integer + intent = in + optional = F +[do_ca] + standard_name = flag_for_cellular_automata + long_name = cellular automata main switch + units = flag + dimensions = () + type = logical + intent = in + optional = F +[ca_global] + standard_name = flag_for_global_cellular_automata + long_name = switch for global ca + units = flag + dimensions = () + type = logical + intent = in + optional = F +[si] + standard_name = sigma_pressure_hybrid_vertical_coordinate + long_name = vertical sigma coordinate for radiation initialization + units = none + dimensions = (vertical_interface_dimension_for_radiation) + type = real + kind = kind_phys + intent = in + optional = F +[vfact_ca] + standard_name = cellular_automata_vertical_weight + long_name = vertical weight for ca + units = frac + dimensions = (vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out + optional = F +[errflg] + standard_name = ccpp_error_flag + long_name = error flag for error handling in CCPP + units = flag + dimensions = () + type = integer + intent = out + optional = F ######################################################################## [ccpp-arg-table] name = GFS_stochastics_run @@ -105,15 +167,6 @@ kind = kind_phys intent = in optional = F -[si] - standard_name = sigma_pressure_hybrid_vertical_coordinate - long_name = vertical sigma coordinate for radiation initialization - units = none - dimensions = (vertical_interface_dimension_for_radiation) - type = real - kind = kind_phys - intent = in - optional = F [vfact_ca] standard_name = cellular_automata_vertical_weight long_name = vertical weight for ca @@ -121,7 +174,7 @@ dimensions = (vertical_layer_dimension) type = real kind = kind_phys - intent = inout + intent = in optional = F [zmtnblck] standard_name = level_of_dividing_streamline @@ -168,42 +221,6 @@ kind = kind_phys intent = in optional = F -[sppt_wts_inv] - standard_name = weights_for_stochastic_sppt_perturbation_flipped - long_name = weights for stochastic sppt perturbation, flipped - units = none - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = inout - optional = F -[skebu_wts_inv] - standard_name = weights_for_stochastic_skeb_perturbation_of_x_wind_flipped - long_name = weights for stochastic skeb perturbation of x wind, flipped - units = none - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = inout - optional = F -[skebv_wts_inv] - standard_name = weights_for_stochastic_skeb_perturbation_of_y_wind_flipped - long_name = weights for stochastic skeb perturbation of y wind, flipped - units = none - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = inout - optional = F -[shum_wts_inv] - standard_name = weights_for_stochastic_shum_perturbation_flipped - long_name = weights for stochastic shum perturbation, flipped - units = none - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = inout - optional = F [diss_est] standard_name = dissipation_estimate_of_air_temperature_at_model_layers long_name = dissipation estimate model layer mean temperature diff --git a/physics/GFS_suite_interstitial.meta b/physics/GFS_suite_interstitial.meta index 2510a5344..3735b887c 100644 --- a/physics/GFS_suite_interstitial.meta +++ b/physics/GFS_suite_interstitial.meta @@ -1613,7 +1613,7 @@ standard_name = liquid_cloud_number_concentration_save long_name = liquid cloud number concentration before entering a physics scheme units = kg kg-1 - dimensions = (horizontal_loop_extent,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = inout @@ -1622,7 +1622,7 @@ standard_name = ice_cloud_number_concentration_save long_name = ice cloud number concentration before entering a physics scheme units = kg kg-1 - dimensions = (horizontal_loop_extent,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = inout @@ -1902,7 +1902,7 @@ standard_name = liquid_cloud_number_concentration_save long_name = liquid cloud number concentration before entering a physics scheme units = kg kg-1 - dimensions = (horizontal_loop_extent,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = in @@ -1911,7 +1911,7 @@ standard_name = ice_cloud_number_concentration_save long_name = ice cloud number concentration before entering a physics scheme units = kg kg-1 - dimensions = (horizontal_loop_extent,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = in @@ -2009,7 +2009,7 @@ standard_name = cumulative_change_of_state_variables long_name = diagnostic tendencies for state variables units = various - dimensions = (horizontal_loop_extent,vertical_dimension,number_of_cumulative_change_processes) + dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_cumulative_change_processes) type = real kind = kind_phys intent = inout diff --git a/physics/GFS_surface_composites.F90 b/physics/GFS_surface_composites.F90 index 966393d03..2f0680541 100644 --- a/physics/GFS_surface_composites.F90 +++ b/physics/GFS_surface_composites.F90 @@ -160,7 +160,7 @@ subroutine GFS_surface_composites_pre_run (im, flag_init, flag_restart, lkm, lsm ! when cplflx is .true. (e.g., for the S2S application). ! Whereas, for the HAFS FV3ATM-HYCOM coupling, cplice is set as ! .false.. In the future HAFS FV3ATM-MOM6 coupling, the cplflx - ! could be .true., while cplice being .false.. + ! could be .true., while cplice being .false.. if (cplice .and. cplflx) then islmsk_cice(i) = 4 flag_cice(i) = .true. @@ -229,7 +229,7 @@ subroutine GFS_surface_composites_pre_run (im, flag_init, flag_restart, lkm, lsm uustar_lnd(i) = uustar(i) weasd_lnd(i) = weasd(i) tsurf_lnd(i) = tsfcl(i) - if (iemsflg == 2 .and. .not. flag_init) then + if (iemsflg == 2 .and. (.not.flag_init .or. flag_restart)) then !-- use land emissivity from the LSM semis_lnd(i) = emis_lnd(i) else @@ -455,7 +455,7 @@ subroutine GFS_surface_composites_post_run ( fh2, cmm, chh, gflx, ep1d, weasd, snowd, tprcp, evap, hflx, qss, tsfc, tsfco, tsfcl, tisfc real(kind=kind_phys), dimension(:), intent(inout) :: hice, cice - real(kind=kind_phys), dimension(:), intent(inout) :: sigmaf, zvfun, hflxq, hffac + real(kind=kind_phys), dimension(:), intent(inout) :: sigmaf, zvfun, hflxq, hffac real(kind=kind_phys), intent(in ) :: h0facu, h0facs real(kind=kind_phys), intent(in ) :: min_seaice real(kind=kind_phys), intent(in ) :: rd, rvrdm1 @@ -613,7 +613,7 @@ subroutine GFS_surface_composites_post_run ( endif ! call stability(z1(i), zvfun(i), gdx, tv1, thv1, wind(i), & ! inputs - z0max, ztmax, tvs, grav, thsfc_loc, & ! inputs + z0max, ztmax, tvs, grav, thsfc_loc, & ! inputs rb(i), ffmm(i), ffhh(i), fm10(i), fh2(i), cd(i), cdq(i), & ! outputs stress(i), uustar(i)) endif ! Checking to see if point is one or multiple surface types diff --git a/physics/GFS_surface_generic.F90 b/physics/GFS_surface_generic.F90 index 1ec7ff784..07d9c1770 100644 --- a/physics/GFS_surface_generic.F90 +++ b/physics/GFS_surface_generic.F90 @@ -16,7 +16,41 @@ module GFS_surface_generic_pre contains - subroutine GFS_surface_generic_pre_init () +!> \section arg_table_GFS_surface_generic_pre_init Argument Table +!! \htmlinclude GFS_surface_generic_pre_init.html +!! + subroutine GFS_surface_generic_pre_init (nthreads, im, slmsk, isot, ivegsrc, stype, vtype, slope, & + vtype_save, stype_save, slope_save, errmsg, errflg) + + implicit none + + ! Interface variables + integer, intent(in) :: nthreads, im, isot, ivegsrc + real(kind_phys), dimension(:), intent(in) :: slmsk + integer, dimension(:), intent(inout) :: vtype, stype, slope + integer, dimension(:), intent(out) :: vtype_save, stype_save, slope_save + + ! CCPP error handling + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Local variables + integer, dimension(1:im) :: islmsk + integer :: i + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + islmsk = nint(slmsk) + + ! Save current values of vegetation, soil and slope type + vtype_save(:) = vtype(:) + stype_save(:) = stype(:) + slope_save(:) = slope(:) + + call update_vegetation_soil_slope_type(nthreads, im, isot, ivegsrc, islmsk, vtype, stype, slope) + end subroutine GFS_surface_generic_pre_init subroutine GFS_surface_generic_pre_finalize() @@ -25,26 +59,27 @@ end subroutine GFS_surface_generic_pre_finalize !> \section arg_table_GFS_surface_generic_pre_run Argument Table !! \htmlinclude GFS_surface_generic_pre_run.html !! - subroutine GFS_surface_generic_pre_run (im, levs, vfrac, islmsk, isot, ivegsrc, stype, vtype, slope, & - prsik_1, prslk_1, tsfc, phil, con_g, & - sigmaf, soiltyp, vegtype, slopetyp, work3, zlvl, & + subroutine GFS_surface_generic_pre_run (nthreads, im, levs, vfrac, islmsk, isot, ivegsrc, stype, vtype, slope, & + prsik_1, prslk_1, tsfc, phil, con_g, sigmaf, work3, zlvl, & drain_cpl, dsnow_cpl, rain_cpl, snow_cpl, lndp_type, n_var_lndp, sfc_wts, & lndp_var_list, lndp_prt_list, & - z01d, zt1d, bexp1d, xlai1d, vegf1d, lndp_vgf, sfc_wts_inv, & + z01d, zt1d, bexp1d, xlai1d, vegf1d, lndp_vgf, & cplflx, flag_cice, islmsk_cice, slimskin_cpl, & - wind, u1, v1, cnvwind, smcwlt2, smcref2, errmsg, errflg) + wind, u1, v1, cnvwind, smcwlt2, smcref2, vtype_save, stype_save, slope_save, & + errmsg, errflg) use surface_perturbation, only: cdfnor implicit none ! Interface variables - integer, intent(in) :: im, levs, isot, ivegsrc + integer, intent(in) :: nthreads, im, levs, isot, ivegsrc integer, dimension(:), intent(in) :: islmsk - integer, dimension(:), intent(inout) :: soiltyp, vegtype, slopetyp real(kind=kind_phys), intent(in) :: con_g - real(kind=kind_phys), dimension(:), intent(in) :: vfrac, stype, vtype, slope, prsik_1, prslk_1 + real(kind=kind_phys), dimension(:), intent(in) :: vfrac, prsik_1, prslk_1 + integer, dimension(:), intent(inout) :: vtype, stype, slope + integer, dimension(:), intent(out) :: vtype_save(:), stype_save(:), slope_save(:) real(kind=kind_phys), dimension(:), intent(inout) :: tsfc real(kind=kind_phys), dimension(:,:), intent(in) :: phil @@ -66,7 +101,6 @@ subroutine GFS_surface_generic_pre_run (im, levs, vfrac, islmsk, isot, ivegsrc, real(kind=kind_phys), dimension(:), intent(out) :: xlai1d real(kind=kind_phys), dimension(:), intent(out) :: vegf1d real(kind=kind_phys), intent(out) :: lndp_vgf - real(kind=kind_phys), dimension(:,:), intent(inout) :: sfc_wts_inv logical, intent(in) :: cplflx real(kind=kind_phys), dimension(:), intent(in) :: slimskin_cpl @@ -99,9 +133,6 @@ subroutine GFS_surface_generic_pre_run (im, levs, vfrac, islmsk, isot, ivegsrc, ! Turn vegetation fraction pattern into percentile pattern lndp_vgf=-999. - if (lndp_type>0) then - sfc_wts_inv(:,:)=sfc_wts(:,:) - endif if (lndp_type==1) then do k =1,n_var_lndp select case(lndp_var_list(k)) @@ -126,32 +157,16 @@ subroutine GFS_surface_generic_pre_run (im, levs, vfrac, islmsk, isot, ivegsrc, ! End of stochastic physics / surface perturbation + ! Save current values of vegetation, soil and slope type + vtype_save(:) = vtype(:) + stype_save(:) = stype(:) + slope_save(:) = slope(:) + + call update_vegetation_soil_slope_type(nthreads, im, isot, ivegsrc, islmsk, vtype, stype, slope) + do i=1,im sigmaf(i) = max(vfrac(i), 0.01_kind_phys) islmsk_cice(i) = islmsk(i) - if (islmsk(i) == 2) then - if (isot == 1) then - soiltyp(i) = 16 - else - soiltyp(i) = 9 - endif - if (ivegsrc == 0 .or. ivegsrc == 4) then - vegtype(i) = 24 - elseif (ivegsrc == 1) then - vegtype(i) = 15 - elseif (ivegsrc == 2) then - vegtype(i) = 13 - elseif (ivegsrc == 3 .or. ivegsrc == 5) then - vegtype(i) = 15 - endif - slopetyp(i) = 9 - else - soiltyp(i) = int( stype(i)+0.5_kind_phys ) - vegtype(i) = int( vtype(i)+0.5_kind_phys ) - slopetyp(i) = int( slope(i)+0.5_kind_phys ) !! clu: slope -> slopetyp - if (vegtype(i) < 1) vegtype(i) = 17 - if (slopetyp(i) < 1) slopetyp(i) = 1 - endif work3(i) = prsik_1(i) / prslk_1(i) @@ -177,6 +192,42 @@ subroutine GFS_surface_generic_pre_run (im, levs, vfrac, islmsk, isot, ivegsrc, end subroutine GFS_surface_generic_pre_run + subroutine update_vegetation_soil_slope_type(nthreads, im, isot, ivegsrc, islmsk, vtype, stype, slope) + + implicit none + + integer, intent(in) :: nthreads, im, isot, ivegsrc, islmsk(:) + integer, intent(inout) :: vtype(:), stype(:), slope(:) + integer :: i + +!$OMP parallel do num_threads(nthreads) default(none) private(i) & +!$OMP shared(im, isot, ivegsrc, islmsk, vtype, stype, slope) + do i=1,im + if (islmsk(i) == 2) then + if (isot == 1) then + stype(i) = 16 + else + stype(i) = 9 + endif + if (ivegsrc == 0 .or. ivegsrc == 4) then + vtype(i) = 24 + elseif (ivegsrc == 1) then + vtype(i) = 15 + elseif (ivegsrc == 2) then + vtype(i) = 13 + elseif (ivegsrc == 3 .or. ivegsrc == 5) then + vtype(i) = 15 + endif + slope(i) = 9 + else + if (vtype(i) < 1) vtype(i) = 17 + if (slope(i) < 1) slope(i) = 1 + endif + enddo +!$OMP end parallel do + + end subroutine update_vegetation_soil_slope_type + end module GFS_surface_generic_pre @@ -194,7 +245,27 @@ module GFS_surface_generic_post contains - subroutine GFS_surface_generic_post_init () +!> \section arg_table_GFS_surface_generic_post_init Argument Table +!! \htmlinclude GFS_surface_generic_post_init.html +!! + subroutine GFS_surface_generic_post_init (vtype, stype, slope, vtype_save, stype_save, slope_save, errmsg, errflg) + + integer, dimension(:), intent(in) :: vtype_save, stype_save, slope_save + integer, dimension(:), intent(out) :: vtype, stype, slope + + ! CCPP error handling + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + ! Restore vegetation, soil and slope type + vtype(:) = vtype_save(:) + stype(:) = stype_save(:) + slope(:) = slope_save(:) + end subroutine GFS_surface_generic_post_init subroutine GFS_surface_generic_post_finalize() @@ -211,7 +282,8 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplchm, cplwav, lssav, dry, dvisdfi_cpl, dnirbm_cpl, dnirdf_cpl, dvisbm_cpl, dvisdf_cpl, nlwsfci_cpl, nlwsfc_cpl, t2mi_cpl, q2mi_cpl, u10mi_cpl, & v10mi_cpl, tsfci_cpl, psurfi_cpl, nnirbmi_cpl, nnirdfi_cpl, nvisbmi_cpl, nvisdfi_cpl, nswsfci_cpl, nswsfc_cpl, nnirbm_cpl, & nnirdf_cpl, nvisbm_cpl, nvisdf_cpl, gflux, evbsa, evcwa, transa, sbsnoa, snowca, snohfa, ep, & - runoff, srunoff, runof, drain, lheatstrg, h0facu, h0facs, zvfun, hflx, evap, hflxq, hffac, errmsg, errflg) + runoff, srunoff, runof, drain, lheatstrg, h0facu, h0facs, zvfun, hflx, evap, hflxq, hffac, & + isot, ivegsrc, islmsk, vtype, stype, slope, vtype_save, stype_save, slope_save, errmsg, errflg) implicit none @@ -241,6 +313,9 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplchm, cplwav, lssav, dry, real(kind=kind_phys), dimension(:), intent(out) :: hflxq real(kind=kind_phys), dimension(:), intent(out) :: hffac + integer, intent(in) :: isot, ivegsrc, islmsk(:), vtype_save(:), stype_save(:), slope_save(:) + integer, intent(out) :: vtype(:), stype(:), slope(:) + ! CCPP error handling variables character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -379,6 +454,11 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplchm, cplwav, lssav, dry, enddo endif + ! Restore vegetation, soil and slope type + vtype(:) = vtype_save(:) + stype(:) = stype_save(:) + slope(:) = slope_save(:) + end subroutine GFS_surface_generic_post_run end module GFS_surface_generic_post diff --git a/physics/GFS_surface_generic.meta b/physics/GFS_surface_generic.meta index 54756c1b4..69f38177d 100644 --- a/physics/GFS_surface_generic.meta +++ b/physics/GFS_surface_generic.meta @@ -3,10 +3,129 @@ type = scheme dependencies = machine.F,surface_perturbation.F90 +######################################################################## +[ccpp-arg-table] + name = GFS_surface_generic_pre_init + type = scheme +[nthreads] + standard_name = number_of_openmp_threads + long_name = number of OpenMP threads available for physics schemes + units = count + dimensions = () + type = integer + intent = in + optional = F +[im] + standard_name = horizontal_dimension + long_name = horizontal dimension + units = count + dimensions = () + type = integer + intent = in + optional = F +[slmsk] + standard_name = area_type + long_name = landmask: sea/land/ice=0/1/2 + units = flag + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[isot] + standard_name = control_for_soil_type_dataset + long_name = soil type dataset choice + units = index + dimensions = () + type = integer + intent = in + optional = F +[ivegsrc] + standard_name = control_for_vegetation_dataset + long_name = land use dataset choice + units = index + dimensions = () + type = integer + intent = in + optional = F +[stype] + standard_name = soil_type_classification + long_name = soil type for lsm + units = index + dimensions = (horizontal_dimension) + type = integer + intent = inout + optional = F +[vtype] + standard_name = vegetation_type_classification + long_name = vegetation type for lsm + units = index + dimensions = (horizontal_dimension) + type = integer + intent = inout + optional = F +[slope] + standard_name = surface_slope_classification + long_name = sfc slope type for lsm + units = index + dimensions = (horizontal_dimension) + type = integer + intent = inout + optional = F +[stype_save] + standard_name = soil_type_classification_save + long_name = soil type for lsm save + units = index + dimensions = (horizontal_dimension) + type = integer + intent = out + optional = F +[vtype_save] + standard_name = vegetation_type_classification_save + long_name = vegetation type for lsm save + units = index + dimensions = (horizontal_dimension) + type = integer + intent = out + optional = F +[slope_save] + standard_name = surface_slope_classification_save + long_name = sfc slope type for lsm save + units = index + dimensions = (horizontal_dimension) + type = integer + intent = out + optional = F +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out + optional = F +[errflg] + standard_name = ccpp_error_flag + long_name = error flag for error handling in CCPP + units = flag + dimensions = () + type = integer + intent = out + optional = F + ######################################################################## [ccpp-arg-table] name = GFS_surface_generic_pre_run type = scheme +[nthreads] + standard_name = number_of_openmp_threads + long_name = number of OpenMP threads available for physics schemes + units = count + dimensions = () + type = integer + intent = in + optional = F [im] standard_name = horizontal_loop_extent long_name = horizontal loop extent @@ -57,31 +176,52 @@ intent = in optional = F [stype] - standard_name = soil_type_classification_real + standard_name = soil_type_classification long_name = soil type for lsm units = index dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = in + type = integer + intent = inout optional = F [vtype] - standard_name = vegetation_type_classification_real + standard_name = vegetation_type_classification long_name = vegetation type for lsm units = index dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = in + type = integer + intent = inout optional = F [slope] - standard_name = surface_slope_classification_real + standard_name = surface_slope_classification long_name = sfc slope type for lsm units = index dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = in + type = integer + intent = inout + optional = F +[vtype_save] + standard_name = vegetation_type_classification_save + long_name = vegetation type for lsm save + units = index + dimensions = (horizontal_loop_extent) + type = integer + intent = out + optional = F +[stype_save] + standard_name = soil_type_classification_save + long_name = soil type for lsm save + units = index + dimensions = (horizontal_loop_extent) + type = integer + intent = out + optional = F +[slope_save] + standard_name = surface_slope_classification_save + long_name = sfc slope type for lsm save + units = index + dimensions = (horizontal_loop_extent) + type = integer + intent = out optional = F [prsik_1] standard_name = surface_dimensionless_exner_function @@ -137,30 +277,6 @@ kind = kind_phys intent = inout optional = F -[soiltyp] - standard_name = soil_type_classification - long_name = soil type at each grid cell - units = index - dimensions = (horizontal_loop_extent) - type = integer - intent = inout - optional = F -[vegtype] - standard_name = vegetation_type_classification - long_name = vegetation type at each grid cell - units = index - dimensions = (horizontal_loop_extent) - type = integer - intent = inout - optional = F -[slopetyp] - standard_name = surface_slope_classification - long_name = surface slope type at each grid cell - units = index - dimensions = (horizontal_loop_extent) - type = integer - intent = inout - optional = F [work3] standard_name = ratio_of_exner_function_between_midlayer_and_interface_at_lowest_model_layer long_name = Exner function ratio bt midlayer and interface at 1st layer @@ -312,15 +428,6 @@ kind = kind_phys intent = out optional = F -[sfc_wts_inv] - standard_name = weights_for_stochastic_surface_physics_perturbation_flipped - long_name = weights for stochastic surface physics perturbation, flipped - units = none - dimensions = (horizontal_loop_extent,number_of_perturbed_land_surface_variables) - type = real - kind = kind_phys - intent = inout - optional = F [cplflx] standard_name = flag_for_surface_flux_coupling long_name = flag controlling cplflx collection (default off) @@ -432,6 +539,76 @@ type = scheme dependencies = machine.F,surface_perturbation.F90 +######################################################################## +[ccpp-arg-table] + name = GFS_surface_generic_post_init + type = scheme +[vtype] + standard_name = vegetation_type_classification + long_name = vegetation type for lsm + units = index + dimensions = (horizontal_dimension) + type = integer + intent = out + optional = F +[stype] + standard_name = soil_type_classification + long_name = soil type for lsm + units = index + dimensions = (horizontal_dimension) + type = integer + intent = out + optional = F +[slope] + standard_name = surface_slope_classification + long_name = sfc slope type for lsm + units = index + dimensions = (horizontal_dimension) + type = integer + intent = out + optional = F +[vtype_save] + standard_name = vegetation_type_classification_save + long_name = vegetation type for lsm save + units = index + dimensions = (horizontal_dimension) + type = integer + intent = in + optional = F +[stype_save] + standard_name = soil_type_classification_save + long_name = soil type for lsm save + units = index + dimensions = (horizontal_dimension) + type = integer + intent = in + optional = F +[slope_save] + standard_name = surface_slope_classification_save + long_name = sfc slope type for lsm save + units = index + dimensions = (horizontal_dimension) + type = integer + intent = in + optional = F +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out + optional = F +[errflg] + standard_name = ccpp_error_flag + long_name = error flag for error handling in CCPP + units = flag + dimensions = () + type = integer + intent = out + optional = F + ######################################################################## [ccpp-arg-table] name = GFS_surface_generic_post_run @@ -1300,6 +1477,78 @@ kind = kind_phys intent = out optional = F +[isot] + standard_name = control_for_soil_type_dataset + long_name = soil type dataset choice + units = index + dimensions = () + type = integer + intent = in + optional = F +[ivegsrc] + standard_name = control_for_vegetation_dataset + long_name = land use dataset choice + units = index + dimensions = () + type = integer + intent = in + optional = F +[islmsk] + standard_name = sea_land_ice_mask + long_name = landmask: sea/land/ice=0/1/2 + units = flag + dimensions = (horizontal_loop_extent) + type = integer + intent = in + optional = F +[vtype] + standard_name = vegetation_type_classification + long_name = vegetation type for lsm + units = index + dimensions = (horizontal_loop_extent) + type = integer + intent = out + optional = F +[stype] + standard_name = soil_type_classification + long_name = soil type for lsm + units = index + dimensions = (horizontal_loop_extent) + type = integer + intent = out + optional = F +[slope] + standard_name = surface_slope_classification + long_name = sfc slope type for lsm + units = index + dimensions = (horizontal_loop_extent) + type = integer + intent = out + optional = F +[vtype_save] + standard_name = vegetation_type_classification_save + long_name = vegetation type for lsm save + units = index + dimensions = (horizontal_loop_extent) + type = integer + intent = in + optional = F +[stype_save] + standard_name = soil_type_classification_save + long_name = soil type for lsm save + units = index + dimensions = (horizontal_loop_extent) + type = integer + intent = in + optional = F +[slope_save] + standard_name = surface_slope_classification_save + long_name = sfc slope type for lsm save + units = index + dimensions = (horizontal_loop_extent) + type = integer + intent = in + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/cires_ugwpv1_oro.F90 b/physics/cires_ugwpv1_oro.F90 index 904731b16..247112bf1 100644 --- a/physics/cires_ugwpv1_oro.F90 +++ b/physics/cires_ugwpv1_oro.F90 @@ -908,8 +908,8 @@ subroutine orogw_v1 (im, km, imx, me, master, dtp, kdt, do_tofd, & dudt_obl(j,k) = -dbim * u1(j,k) dvdt_obl(j,k) = -dbim * v1(j,k) - pdvdt(j,k) = dudt_obl(j,k) +pdvdt(j,k) - pdudt(j,k) = dvdt_obl(j,k) +pdudt(j,k) + pdudt(j,k) = dudt_obl(j,k) +pdudt(j,k) + pdvdt(j,k) = dvdt_obl(j,k) +pdvdt(j,k) du_oblcol(j) = du_oblcol(j) + dudt_obl(j,k)* del(j,k) dv_oblcol(j) = dv_oblcol(j) + dvdt_obl(j,k)* del(j,k) dusfc(j) = dusfc(j) + du_oblcol(j) diff --git a/physics/drag_suite.F90 b/physics/drag_suite.F90 index 9b110d689..e997122f7 100644 --- a/physics/drag_suite.F90 +++ b/physics/drag_suite.F90 @@ -407,7 +407,7 @@ subroutine drag_suite_run( & ! Large-scale GWD + blocking real(kind=kind_phys), parameter :: dxmin_ls = 3000., & & dxmax_ls = 13000. ! min,max range of tapering (m) - real(kind=kind_phys) :: ss_taper, ls_taper ! small- and large-scale tapering factors (-) + real(kind=kind_phys), dimension(im) :: ss_taper, ls_taper ! small- and large-scale tapering factors (-) ! ! Variables for limiting topographic standard deviation (var) real(kind=kind_phys), parameter :: varmax_ss = 50., & @@ -488,7 +488,8 @@ subroutine drag_suite_run( & real(kind=kind_phys) :: cd real(kind=kind_phys) :: zblk,tautem real(kind=kind_phys) :: pe,ke - real(kind=kind_phys) :: delx,dely,dxy4(4),dxy4p(4) + real(kind=kind_phys) :: delx,dely + real(kind=kind_phys) :: dxy4(im,4),dxy4p(im,4) real(kind=kind_phys) :: dxy(im),dxyp(im) real(kind=kind_phys) :: ol4p(4),olp(im),od(im) real(kind=kind_phys) :: taufb(im,km+1) @@ -568,52 +569,46 @@ subroutine drag_suite_run( & RDXZB(i) = 0.0 enddo -!temporary use of large-scale data: -! do i=1,im -! varss(i)=var(i) -! oc1ss(i)=oc1(i) -! do j=1,4 -! oa4ss(i,j)=oa4(i,j) -! ol4ss(i,j)=ol4(i,j) -! enddo -! enddo -! !--- calculate scale-aware tapering factors -!NOTE: if dx(1) is not representative of most/all dx, this needs to change... -if ( dx(1) .ge. dxmax_ls ) then - ls_taper = 1. -else - if ( dx(1) .le. dxmin_ls) then - ls_taper = 0. +do i=1,im + if ( dx(i) .ge. dxmax_ls ) then + ls_taper(i) = 1. else - ls_taper = 0.5 * ( SIN(pi*(dx(1)-0.5*(dxmax_ls+dxmin_ls))/ & - (dxmax_ls-dxmin_ls)) + 1. ) - end if -end if -! if (me==master) print *,"in Drag Suite, dx(1:2):",dx(1),dx(2) -if ( dx(1) .ge. dxmax_ss ) then - ss_taper = 1. -else - if ( dx(1) .le. dxmin_ss) then - ss_taper = 0. + if ( dx(i) .le. dxmin_ls) then + ls_taper(i) = 0. + else + ls_taper(i) = 0.5 * ( SIN(pi*(dx(i)-0.5*(dxmax_ls+dxmin_ls))/ & + (dxmax_ls-dxmin_ls)) + 1. ) + endif + endif +enddo + +do i=1,im + if ( dx(i) .ge. dxmax_ss ) then + ss_taper(i) = 1. else - ss_taper = dxmax_ss * (1. - dxmin_ss/dx(1))/(dxmax_ss-dxmin_ss) - end if -end if -! if (me==master) print *,"in Drag Suite, ss_taper:",ss_taper + if ( dx(i) .le. dxmin_ss) then + ss_taper(i) = 0. + else + ss_taper(i) = dxmax_ss * (1. - dxmin_ss/dx(i))/(dxmax_ss-dxmin_ss) + endif + endif +enddo !--- calculate length of grid for flow-blocking drag ! - delx = dx(1) - dely = dx(1) - dxy4(1) = delx - dxy4(2) = dely - dxy4(3) = sqrt(delx*delx + dely*dely) - dxy4(4) = dxy4(3) - dxy4p(1) = dxy4(2) - dxy4p(2) = dxy4(1) - dxy4p(3) = dxy4(4) - dxy4p(4) = dxy4(3) +do i=1,im + delx = dx(i) + dely = dx(i) + dxy4(i,1) = delx + dxy4(i,2) = dely + dxy4(i,3) = sqrt(delx*delx + dely*dely) + dxy4(i,4) = dxy4(i,3) + dxy4p(i,1) = dxy4(i,2) + dxy4p(i,2) = dxy4(i,1) + dxy4p(i,3) = dxy4(i,4) + dxy4p(i,4) = dxy4(i,3) +enddo ! !-----initialize arrays ! @@ -794,148 +789,132 @@ subroutine drag_suite_run( & ! !----- compute length of grid in the along(dxy) and cross(dxyp) wind directions ! - dxy(i) = dxy4(MOD(nwd-1,4)+1) - dxyp(i) = dxy4p(MOD(nwd-1,4)+1) + dxy(i) = dxy4(i,MOD(nwd-1,4)+1) + dxyp(i) = dxy4p(i,MOD(nwd-1,4)+1) enddo ! ! END INITIALIZATION; BEGIN GWD CALCULATIONS: ! IF ( (do_gsl_drag_ls_bl).and. & - ((gwd_opt_ls .EQ. 1).or.(gwd_opt_bl .EQ. 1)).and. & - (ls_taper .GT. 1.E-02) ) THEN !==== + ((gwd_opt_ls .EQ. 1).or.(gwd_opt_bl .EQ. 1)) ) then + + do i=its,im + + if ( ls_taper(i).GT.1.E-02 ) then + ! !--- saving richardson number in usqj for migwdi ! - do k = kts,km-1 - do i = its,im - ti = 2.0 / (t1(i,k)+t1(i,k+1)) - rdz = 1./(zl(i,k+1) - zl(i,k)) - tem1 = u1(i,k) - u1(i,k+1) - tem2 = v1(i,k) - v1(i,k+1) - dw2 = rcl*(tem1*tem1 + tem2*tem2) - shr2 = max(dw2,dw2min) * rdz * rdz - bvf2 = g*(g/cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti - usqj(i,k) = max(bvf2/shr2,rimin) - bnv2(i,k) = 2.0*g*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k)) - bnv2(i,k) = max( bnv2(i,k), bnv2min ) - enddo - enddo + do k = kts,km-1 + ti = 2.0 / (t1(i,k)+t1(i,k+1)) + rdz = 1./(zl(i,k+1) - zl(i,k)) + tem1 = u1(i,k) - u1(i,k+1) + tem2 = v1(i,k) - v1(i,k+1) + dw2 = rcl*(tem1*tem1 + tem2*tem2) + shr2 = max(dw2,dw2min) * rdz * rdz + bvf2 = g*(g/cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti + usqj(i,k) = max(bvf2/shr2,rimin) + bnv2(i,k) = 2.0*g*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k)) + bnv2(i,k) = max( bnv2(i,k), bnv2min ) + enddo ! !----compute the "low level" or 1/3 wind magnitude (m/s) ! - do i = its,im - ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0) - rulow(i) = 1./ulow(i) - enddo + ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0) + rulow(i) = 1./ulow(i) ! - do k = kts,km-1 - do i = its,im - velco(i,k) = (0.5*rcs) * ((u1(i,k)+u1(i,k+1)) * ubar(i) & - + (v1(i,k)+v1(i,k+1)) * vbar(i)) - velco(i,k) = velco(i,k) * rulow(i) - if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0.)) then - velco(i,k) = veleps - endif - enddo - enddo + do k = kts,km-1 + velco(i,k) = (0.5*rcs) * ((u1(i,k)+u1(i,k+1)) * ubar(i) & + + (v1(i,k)+v1(i,k+1)) * vbar(i)) + velco(i,k) = velco(i,k) * rulow(i) + if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0.)) then + velco(i,k) = veleps + endif + enddo ! ! no drag when critical level in the base layer ! - do i = its,im - ldrag(i) = velco(i,1).le.0. - enddo + ldrag(i) = velco(i,1).le.0. ! ! no drag when velco.lt.0 ! - do k = kpblmin,kpblmax - do i = its,im - if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0. - enddo - enddo + do k = kpblmin,kpblmax + if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0. + enddo ! ! no drag when bnv2.lt.0 ! - do k = kts,kpblmax - do i = its,im - if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. bnv2(i,k).lt.0. - enddo - enddo + do k = kts,kpblmax + if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. bnv2(i,k).lt.0. + enddo ! !-----the low level weighted average ri is stored in usqj(1,1; im) !-----the low level weighted average n**2 is stored in bnv2(1,1; im) !---- this is called bnvl2 in phys_gwd_alpert_sub not bnv2 !---- rdelks (del(k)/delks) vert ave factor so we can * instead of / ! - do i = its,im - wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i) - bnv2(i,1) = wtkbj * bnv2(i,1) - usqj(i,1) = wtkbj * usqj(i,1) - enddo + wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i) + bnv2(i,1) = wtkbj * bnv2(i,1) + usqj(i,1) = wtkbj * usqj(i,1) ! - do k = kpblmin,kpblmax - do i = its,im - if (k .lt. kbl(i)) then - rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i) - bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks - usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks - endif - enddo - enddo + do k = kpblmin,kpblmax + if (k .lt. kbl(i)) then + rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i) + bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks + usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks + endif + enddo ! - do i = its,im - ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0 - ldrag(i) = ldrag(i) .or. ulow(i).eq.1.0 - ldrag(i) = ldrag(i) .or. var(i) .le. 0.0 - enddo + ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0 + ldrag(i) = ldrag(i) .or. ulow(i).eq.1.0 + ldrag(i) = ldrag(i) .or. var(i) .le. 0.0 ! ! set all ri low level values to the low level value ! - do k = kpblmin,kpblmax - do i = its,im - if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1) - enddo - enddo + do k = kpblmin,kpblmax + if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1) + enddo ! - do i = its,im - if (.not.ldrag(i)) then - bnv(i) = sqrt( bnv2(i,1) ) - fr(i) = bnv(i) * rulow(i) * 2. * var(i) * od(i) - fr(i) = min(fr(i),frmax) - xn(i) = ubar(i) * rulow(i) - yn(i) = vbar(i) * rulow(i) - endif - enddo + if (.not.ldrag(i)) then + bnv(i) = sqrt( bnv2(i,1) ) + fr(i) = bnv(i) * rulow(i) * 2. * var(i) * od(i) + fr(i) = min(fr(i),frmax) + xn(i) = ubar(i) * rulow(i) + yn(i) = vbar(i) * rulow(i) + endif ! ! compute the base level stress and store it in taub ! calculate enhancement factor, number of mountains & aspect ! ratio const. use simplified relationship between standard ! deviation & critical hgt - do i = its,im - if (.not. ldrag(i)) then - efact = (oa(i) + 2.) ** (ce*fr(i)/frc) - efact = min( max(efact,efmin), efmax ) + if (.not. ldrag(i)) then + efact = (oa(i) + 2.) ** (ce*fr(i)/frc) + efact = min( max(efact,efmin), efmax ) !!!!!!! cleff (effective grid length) is highly tunable parameter !!!!!!! the bigger (smaller) value produce weaker (stronger) wave drag -!WRF cleff = sqrt(dxy(i)**2. + dxyp(i)**2.) -!WRF cleff = 3. * max(dx(i),cleff) - coefm(i) = (1. + ol(i)) ** (oa(i)+1.) -!WRF xlinv(i) = coefm(i) / cleff - xlinv(i) = coefm(i) * cleff - tem = fr(i) * fr(i) * oc1(i) - gfobnv = gmax * tem / ((tem + cg)*bnv(i)) - if ( gwd_opt_ls .NE. 0 ) then - taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i) & - * ulow(i) * gfobnv * efact - else ! We've gotten what we need for the blocking scheme - taub(i) = 0.0 - end if - else - taub(i) = 0.0 - xn(i) = 0.0 - yn(i) = 0.0 - endif - enddo +!WRF cleff = sqrt(dxy(i)**2. + dxyp(i)**2.) +!WRF cleff = 3. * max(dx(i),cleff) + coefm(i) = (1. + ol(i)) ** (oa(i)+1.) +!WRF xlinv(i) = coefm(i) / cleff + xlinv(i) = coefm(i) * cleff + tem = fr(i) * fr(i) * oc1(i) + gfobnv = gmax * tem / ((tem + cg)*bnv(i)) + if ( gwd_opt_ls .NE. 0 ) then + taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i) & + * ulow(i) * gfobnv * efact + else ! We've gotten what we need for the blocking scheme + taub(i) = 0.0 + end if + else + taub(i) = 0.0 + xn(i) = 0.0 + yn(i) = 0.0 + endif + + endif ! (ls_taper(i).GT.1.E-02) + + enddo ! do i=its,im ENDIF ! (do_gsl_drag_ls_bl).and.((gwd_opt_ls .EQ. 1).or.(gwd_opt_bl .EQ. 1)) @@ -949,412 +928,409 @@ subroutine drag_suite_run( & utendwave=0. vtendwave=0. ! - IF ( (do_gsl_drag_ss).and.(ss_taper.GT.1.E-02) ) THEN - ! if (me==master) print *,"in Drag Suite: Running small-scale gravity wave drag" -! -! declaring potential temperature -! - do k = kts,km - do i = its,im - thx(i,k) = t1(i,k)/prslk(i,k) - enddo - enddo -! - do k = kts,km - do i = its,im - tvcon = (1.+fv*q1(i,k)) - thvx(i,k) = thx(i,k)*tvcon - enddo - enddo - - do i=its,im - hpbl2 = hpbl(i)+10. - kpbl2 = kpbl(i) - !kvar = MIN(kpbl, k-level of var) - kvar = 1 - do k=kts+1,MAX(kpbl(i),kts+1) -! IF (zl(i,k)>2.*var(i) .or. zl(i,k)>2*varmax) then - IF (zl(i,k)>300.) then - kpbl2 = k - IF (k == kpbl(i)) then - hpbl2 = hpbl(i)+10. - ELSE - hpbl2 = zl(i,k)+10. - ENDIF - exit - ENDIF - enddo - if((xland(i)-1.5).le.0. .and. 2.*varss(i).le.hpbl(i))then - if(br1(i).gt.0. .and. thvx(i,kpbl2)-thvx(i,kts) > 0.)then - cleff_ss = sqrt(dxy(i)**2 + dxyp(i)**2) ! WRF -! cleff_ss = 3. * max(dx(i),cleff_ss) -! cleff_ss = 10. * max(dxmax_ss,cleff_ss) - cleff_ss = 0.1 * max(dxmax_ss,cleff_ss) ! WRF -! cleff_ss = 0.1 * 12000. - coefm_ss(i) = (1. + olss(i)) ** (oass(i)+1.) - xlinv(i) = coefm_ss(i) / cleff_ss - !govrth(i)=g/(0.5*(thvx(i,kpbl(i))+thvx(i,kts))) - govrth(i)=g/(0.5*(thvx(i,kpbl2)+thvx(i,kts))) - !XNBV=sqrt(govrth(i)*(thvx(i,kpbl(i))-thvx(i,kts))/hpbl(i)) - XNBV=sqrt(govrth(i)*(thvx(i,kpbl2)-thvx(i,kts))/hpbl2) -! - !if(abs(XNBV/u1(i,kpbl(i))).gt.xlinv(i))then - if(abs(XNBV/u1(i,kpbl2)).gt.xlinv(i))then - !tauwavex0=0.5*XNBV*xlinv(i)*(2*MIN(varss(i),75.))**2*ro(i,kts)*u1(i,kpbl(i)) - !tauwavex0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*u1(i,kpbl2) - !tauwavex0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*u1(i,3) - var_temp = MIN(varss(i),varmax_ss) + & - MAX(0.,beta_ss*(varss(i)-varmax_ss)) - ! Note: This is a semi-implicit treatment of the time differencing - var_temp2 = 0.5*XNBV*xlinv(i)*(2.*var_temp)**2*ro(i,kvar) ! this is greater than zero - tauwavex0=-var_temp2*u1(i,kvar)/(1.+var_temp2*deltim) - tauwavex0=tauwavex0*ss_taper - else - tauwavex0=0. - endif +IF ( do_gsl_drag_ss ) THEN + + do i=its,im + + if ( ss_taper(i).GT.1.E-02 ) then + ! + ! calculating potential temperature + ! + do k = kts,km + thx(i,k) = t1(i,k)/prslk(i,k) + enddo + ! + do k = kts,km + tvcon = (1.+fv*q1(i,k)) + thvx(i,k) = thx(i,k)*tvcon + enddo + + hpbl2 = hpbl(i)+10. + kpbl2 = kpbl(i) + !kvar = MIN(kpbl, k-level of var) + kvar = 1 + do k=kts+1,MAX(kpbl(i),kts+1) +! IF (zl(i,k)>2.*var(i) .or. zl(i,k)>2*varmax) then + IF (zl(i,k)>300.) then + kpbl2 = k + IF (k == kpbl(i)) then + hpbl2 = hpbl(i)+10. + ELSE + hpbl2 = zl(i,k)+10. + ENDIF + exit + ENDIF + enddo + if((xland(i)-1.5).le.0. .and. 2.*varss(i).le.hpbl(i))then + if(br1(i).gt.0. .and. thvx(i,kpbl2)-thvx(i,kts) > 0.)then + cleff_ss = sqrt(dxy(i)**2 + dxyp(i)**2) ! WRF +! cleff_ss = 3. * max(dx(i),cleff_ss) +! cleff_ss = 10. * max(dxmax_ss,cleff_ss) + cleff_ss = 0.1 * max(dxmax_ss,cleff_ss) ! WRF +! cleff_ss = 0.1 * 12000. + coefm_ss(i) = (1. + olss(i)) ** (oass(i)+1.) + xlinv(i) = coefm_ss(i) / cleff_ss + !govrth(i)=g/(0.5*(thvx(i,kpbl(i))+thvx(i,kts))) + govrth(i)=g/(0.5*(thvx(i,kpbl2)+thvx(i,kts))) + !XNBV=sqrt(govrth(i)*(thvx(i,kpbl(i))-thvx(i,kts))/hpbl(i)) + XNBV=sqrt(govrth(i)*(thvx(i,kpbl2)-thvx(i,kts))/hpbl2) +! + !if(abs(XNBV/u1(i,kpbl(i))).gt.xlinv(i))then + if(abs(XNBV/u1(i,kpbl2)).gt.xlinv(i))then + !tauwavex0=0.5*XNBV*xlinv(i)*(2*MIN(varss(i),75.))**2*ro(i,kts)*u1(i,kpbl(i)) + !tauwavex0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*u1(i,kpbl2) + !tauwavex0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*u1(i,3) + var_temp = MIN(varss(i),varmax_ss) + & + MAX(0.,beta_ss*(varss(i)-varmax_ss)) + ! Note: This is a semi-implicit treatment of the time differencing + var_temp2 = 0.5*XNBV*xlinv(i)*(2.*var_temp)**2*ro(i,kvar) ! this is greater than zero + tauwavex0=-var_temp2*u1(i,kvar)/(1.+var_temp2*deltim) + tauwavex0=tauwavex0*ss_taper(i) + else + tauwavex0=0. + endif ! - !if(abs(XNBV/v1(i,kpbl(i))).gt.xlinv(i))then - if(abs(XNBV/v1(i,kpbl2)).gt.xlinv(i))then - !tauwavey0=0.5*XNBV*xlinv(i)*(2*MIN(varss(i),75.))**2*ro(i,kts)*v1(i,kpbl(i)) - !tauwavey0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*v1(i,kpbl2) - !tauwavey0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*v1(i,3) - var_temp = MIN(varss(i),varmax_ss) + & - MAX(0.,beta_ss*(varss(i)-varmax_ss)) - ! Note: This is a semi-implicit treatment of the time differencing - tauwavey0=-var_temp2*v1(i,kvar)/(1.+var_temp2*deltim) - tauwavey0=tauwavey0*ss_taper - else - tauwavey0=0. - endif + !if(abs(XNBV/v1(i,kpbl(i))).gt.xlinv(i))then + if(abs(XNBV/v1(i,kpbl2)).gt.xlinv(i))then + !tauwavey0=0.5*XNBV*xlinv(i)*(2*MIN(varss(i),75.))**2*ro(i,kts)*v1(i,kpbl(i)) + !tauwavey0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*v1(i,kpbl2) + !tauwavey0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*v1(i,3) + var_temp = MIN(varss(i),varmax_ss) + & + MAX(0.,beta_ss*(varss(i)-varmax_ss)) + ! Note: This is a semi-implicit treatment of the time differencing + var_temp2 = 0.5*XNBV*xlinv(i)*(2.*var_temp)**2*ro(i,kvar) ! this is greater than zero + tauwavey0=-var_temp2*v1(i,kvar)/(1.+var_temp2*deltim) + tauwavey0=tauwavey0*ss_taper(i) + else + tauwavey0=0. + endif - do k=kts,kpbl(i) !MIN(kpbl2+1,km-1) + do k=kts,kpbl(i) !MIN(kpbl2+1,km-1) !original - !utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl(i)),0.)/hpbl(i) - !vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl(i)),0.)/hpbl(i) + !utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl(i)),0.)/hpbl(i) + !vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl(i)),0.)/hpbl(i) !new - utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2 - vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2 + utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2 + vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2 !mod-to be used in HRRRv3/RAPv4 - !utendwave(i,k)=-1.*tauwavex0 * max((1.-zl(i,k)/hpbl2),0.)**2 - !vtendwave(i,k)=-1.*tauwavey0 * max((1.-zl(i,k)/hpbl2),0.)**2 + !utendwave(i,k)=-1.*tauwavex0 * max((1.-zl(i,k)/hpbl2),0.)**2 + !vtendwave(i,k)=-1.*tauwavey0 * max((1.-zl(i,k)/hpbl2),0.)**2 + enddo + endif + endif + + do k = kts,km + dudt(i,k) = dudt(i,k) + utendwave(i,k) + dvdt(i,k) = dvdt(i,k) + vtendwave(i,k) + dusfc(i) = dusfc(i) + utendwave(i,k) * del(i,k) + dvsfc(i) = dvsfc(i) + vtendwave(i,k) * del(i,k) + enddo + if(udtend>0) then + dtend(i,kts:km,udtend) = dtend(i,kts:km,udtend) + utendwave(i,kts:km)*deltim + endif + if(vdtend>0) then + dtend(i,kts:km,vdtend) = dtend(i,kts:km,vdtend) + vtendwave(i,kts:km)*deltim + endif + if ( (gwd_opt == 33).or.(gwd_opt == 22) ) then + do k = kts,km + dusfc_ss(i) = dusfc_ss(i) + utendwave(i,k) * del(i,k) + dvsfc_ss(i) = dvsfc_ss(i) + vtendwave(i,k) * del(i,k) + dtaux2d_ss(i,k) = utendwave(i,k) + dtauy2d_ss(i,k) = vtendwave(i,k) enddo - endif - endif - enddo ! end i loop + endif - do k = kts,km - do i = its,im - dudt(i,k) = dudt(i,k) + utendwave(i,k) - dvdt(i,k) = dvdt(i,k) + vtendwave(i,k) - dusfc(i) = dusfc(i) + utendwave(i,k) * del(i,k) - dvsfc(i) = dvsfc(i) + vtendwave(i,k) * del(i,k) - enddo - enddo - if(udtend>0) then - dtend(its:im,kts:km,udtend) = dtend(its:im,kts:km,udtend) + utendwave(its:im,kts:km)*deltim - endif - if(vdtend>0) then - dtend(its:im,kts:km,vdtend) = dtend(its:im,kts:km,vdtend) + vtendwave(its:im,kts:km)*deltim - endif - if ( (gwd_opt == 33).or.(gwd_opt == 22) ) then - do k = kts,km - do i = its,im - dusfc_ss(i) = dusfc_ss(i) + utendwave(i,k) * del(i,k) - dvsfc_ss(i) = dvsfc_ss(i) + vtendwave(i,k) * del(i,k) - dtaux2d_ss(i,k) = utendwave(i,k) - dtauy2d_ss(i,k) = vtendwave(i,k) - enddo - enddo - endif + endif ! if (ss_taper(i).GT.1.E-02) + + enddo ! i=its,im ENDIF ! if (do_gsl_drag_ss) !================================================================ ! Topographic Form Drag from Beljaars et al. (2004, QJRMS, equ. 16): !================================================================ -IF ( (do_gsl_drag_tofd).and.(ss_taper.GT.1.E-02) ) THEN - ! if (me==master) print *,"in Drag Suite: Running form drag" - - utendform=0. - vtendform=0. - - DO i=its,im - IF ((xland(i)-1.5) .le. 0.) then - !(IH*kflt**n1)**-1 = (0.00102*0.00035**-1.9)**-1 = 0.00026615161 - var_temp = MIN(varss(i),varmax_fd) + & - MAX(0.,beta_fd*(varss(i)-varmax_fd)) - var_temp = MIN(var_temp, 250.) - a1=0.00026615161*var_temp**2 -! a1=0.00026615161*MIN(varss(i),varmax)**2 -! a1=0.00026615161*(0.5*varss(i))**2 - ! k1**(n1-n2) = 0.003**(-1.9 - -2.8) = 0.003**0.9 = 0.005363 - a2=a1*0.005363 - ! Revise e-folding height based on PBL height and topographic std. dev. -- M. Toy 3/12/2018 - H_efold = max(2*varss(i),hpbl(i)) - H_efold = min(H_efold,1500.) - DO k=kts,km - wsp=SQRT(u1(i,k)**2 + v1(i,k)**2) - ! alpha*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759 - var_temp = 0.0759*EXP(-(zl(i,k)/H_efold)**1.5)*a2* & - zl(i,k)**(-1.2)*ss_taper ! this is greater than zero - ! Note: This is a semi-implicit treatment of the time differencing - ! per Beljaars et al. (2004, QJRMS) - utendform(i,k) = - var_temp*wsp*u1(i,k)/(1. + var_temp*deltim*wsp) - vtendform(i,k) = - var_temp*wsp*v1(i,k)/(1. + var_temp*deltim*wsp) - !IF(zl(i,k) > 4000.) exit - ENDDO - ENDIF - ENDDO +IF ( do_gsl_drag_tofd ) THEN - do k = kts,km - do i = its,im - dudt(i,k) = dudt(i,k) + utendform(i,k) - dvdt(i,k) = dvdt(i,k) + vtendform(i,k) - dusfc(i) = dusfc(i) + utendform(i,k) * del(i,k) - dvsfc(i) = dvsfc(i) + vtendform(i,k) * del(i,k) - enddo - enddo - if(udtend>0) then - dtend(its:im,kts:km,udtend) = dtend(its:im,kts:km,udtend) + utendform(its:im,kts:km)*deltim - endif - if(vdtend>0) then - dtend(its:im,kts:km,vdtend) = dtend(its:im,kts:km,vdtend) + vtendform(its:im,kts:km)*deltim - endif - if ( (gwd_opt == 33).or.(gwd_opt == 22) ) then - do k = kts,km - do i = its,im - dtaux2d_fd(i,k) = utendform(i,k) - dtauy2d_fd(i,k) = vtendform(i,k) - dusfc_fd(i) = dusfc_fd(i) + utendform(i,k) * del(i,k) - dvsfc_fd(i) = dvsfc_fd(i) + vtendform(i,k) * del(i,k) - enddo - enddo - endif + do i=its,im + + if ( ss_taper(i).GT.1.E-02 ) then + + utendform=0. + vtendform=0. + + IF ((xland(i)-1.5) .le. 0.) then + !(IH*kflt**n1)**-1 = (0.00102*0.00035**-1.9)**-1 = 0.00026615161 + var_temp = MIN(varss(i),varmax_fd) + & + MAX(0.,beta_fd*(varss(i)-varmax_fd)) + var_temp = MIN(var_temp, 250.) + a1=0.00026615161*var_temp**2 +! a1=0.00026615161*MIN(varss(i),varmax)**2 +! a1=0.00026615161*(0.5*varss(i))**2 + ! k1**(n1-n2) = 0.003**(-1.9 - -2.8) = 0.003**0.9 = 0.005363 + a2=a1*0.005363 + ! Revise e-folding height based on PBL height and topographic std. dev. -- M. Toy 3/12/2018 + H_efold = max(2*varss(i),hpbl(i)) + H_efold = min(H_efold,1500.) + DO k=kts,km + wsp=SQRT(u1(i,k)**2 + v1(i,k)**2) + ! alpha*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759 + var_temp = 0.0759*EXP(-(zl(i,k)/H_efold)**1.5)*a2* & + zl(i,k)**(-1.2)*ss_taper(i) ! this is greater than zero + ! Note: This is a semi-implicit treatment of the time differencing + ! per Beljaars et al. (2004, QJRMS) + utendform(i,k) = - var_temp*wsp*u1(i,k)/(1. + var_temp*deltim*wsp) + vtendform(i,k) = - var_temp*wsp*v1(i,k)/(1. + var_temp*deltim*wsp) + !IF(zl(i,k) > 4000.) exit + ENDDO + ENDIF + + do k = kts,km + dudt(i,k) = dudt(i,k) + utendform(i,k) + dvdt(i,k) = dvdt(i,k) + vtendform(i,k) + dusfc(i) = dusfc(i) + utendform(i,k) * del(i,k) + dvsfc(i) = dvsfc(i) + vtendform(i,k) * del(i,k) + enddo + if(udtend>0) then + dtend(i,kts:km,udtend) = dtend(i,kts:km,udtend) + utendform(i,kts:km)*deltim + endif + if(vdtend>0) then + dtend(i,kts:km,vdtend) = dtend(i,kts:km,vdtend) + vtendform(i,kts:km)*deltim + endif + if ( (gwd_opt == 33).or.(gwd_opt == 22) ) then + do k = kts,km + dtaux2d_fd(i,k) = utendform(i,k) + dtauy2d_fd(i,k) = vtendform(i,k) + dusfc_fd(i) = dusfc_fd(i) + utendform(i,k) * del(i,k) + dvsfc_fd(i) = dvsfc_fd(i) + vtendform(i,k) * del(i,k) + enddo + endif + + endif ! if (ss_taper(i).GT.1.E-02) + + enddo ! i=its,im ENDIF ! if (do_gsl_drag_tofd) !======================================================= ! More for the large-scale gwd component -IF ( (do_gsl_drag_ls_bl).and. & - (gwd_opt_ls .EQ. 1).and.(ls_taper.GT.1.E-02) ) THEN - ! if (me==master) print *,"in Drag Suite: Running large-scale gravity wave drag" +IF ( (do_gsl_drag_ls_bl).and.(gwd_opt_ls .EQ. 1) ) THEN + + do i=its,im + + if ( ls_taper(i).GT.1.E-02 ) then + ! ! now compute vertical structure of the stress. - do k = kts,kpblmax - do i = its,im - if (k .le. kbl(i)) taup(i,k) = taub(i) - enddo - enddo + do k = kts,kpblmax + if (k .le. kbl(i)) taup(i,k) = taub(i) + enddo ! - do k = kpblmin, km-1 ! vertical level k loop! - kp1 = k + 1 - do i = its,im + do k = kpblmin, km-1 ! vertical level k loop! + kp1 = k + 1 ! ! unstablelayer if ri < ric ! unstable layer if upper air vel comp along surf vel <=0 (crit lay) ! at (u-c)=0. crit layer exists and bit vector should be set (.le.) ! - if (k .ge. kbl(i)) then - icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric) & - .or. (velco(i,k) .le. 0.0) - brvf(i) = max(bnv2(i,k),bnv2min) ! brunt-vaisala frequency squared - brvf(i) = sqrt(brvf(i)) ! brunt-vaisala frequency - endif - enddo + if (k .ge. kbl(i)) then + icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric) & + .or. (velco(i,k) .le. 0.0) + brvf(i) = max(bnv2(i,k),bnv2min) ! brunt-vaisala frequency squared + brvf(i) = sqrt(brvf(i)) ! brunt-vaisala frequency + endif ! - do i = its,im - if (k .ge. kbl(i) .and. (.not. ldrag(i))) then - if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0 ) then - temv = 1.0 / velco(i,k) - tem1 = coefm(i)/dxy(i)*(ro(i,kp1)+ro(i,k))*brvf(i)*velco(i,k)*0.5 - hd = sqrt(taup(i,k) / tem1) - fro = brvf(i) * hd * temv + if (k .ge. kbl(i) .and. (.not. ldrag(i))) then + if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0 ) then + temv = 1.0 / velco(i,k) + tem1 = coefm(i)/dxy(i)*(ro(i,kp1)+ro(i,k))*brvf(i)* & + velco(i,k)*0.5 + hd = sqrt(taup(i,k) / tem1) + fro = brvf(i) * hd * temv ! ! rim is the minimum-richardson number by shutts (1985) - tem2 = sqrt(usqj(i,k)) - tem = 1. + tem2 * fro - rim = usqj(i,k) * (1.-fro) / (tem * tem) + tem2 = sqrt(usqj(i,k)) + tem = 1. + tem2 * fro + rim = usqj(i,k) * (1.-fro) / (tem * tem) ! ! check stability to employ the 'saturation hypothesis' ! of lindzen (1981) except at tropospheric downstream regions ! - if (rim .le. ric) then ! saturation hypothesis! - if ((oa(i) .le. 0.).or.(kp1 .ge. kpblmin )) then - temc = 2.0 + 1.0 / tem2 - hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf(i) - taup(i,kp1) = tem1 * hd * hd - endif - else ! no wavebreaking! - taup(i,kp1) = taup(i,k) + if (rim .le. ric) then ! saturation hypothesis! + if ((oa(i) .le. 0.).or.(kp1 .ge. kpblmin )) then + temc = 2.0 + 1.0 / tem2 + hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf(i) + taup(i,kp1) = tem1 * hd * hd + endif + else ! no wavebreaking! + taup(i,kp1) = taup(i,k) + endif + endif endif - endif - endif - enddo - enddo -! - if(lcap.lt.km) then - do klcap = lcapp1,km - do i = its,im - taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap) enddo - enddo - endif +! + if(lcap.lt.km) then + do klcap = lcapp1,km + taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap) + enddo + endif + + endif ! if ( ls_taper(i).GT.1.E-02 ) + + enddo ! do i=its,im -ENDIF !END LARGE-SCALE TAU CALCULATION +ENDIF ! (do_gsl_drag_ls_bl).and.(gwd_opt_ls .EQ. 1) !=============================================================== !COMPUTE BLOCKING COMPONENT !=============================================================== -IF ( (do_gsl_drag_ls_bl) .and. & - (gwd_opt_bl .EQ. 1) .and. (ls_taper .GT. 1.E-02) ) THEN - ! if (me==master) print *,"in Drag Suite: Running blocking drag" +IF ( (do_gsl_drag_ls_bl) .and. (gwd_opt_bl .EQ. 1) ) THEN - do i = its,im - if(.not.ldrag(i)) then + do i=its,im + + if ( ls_taper(i).GT.1.E-02 ) then + + if (.not.ldrag(i)) then ! !------- determine the height of flow-blocking layer ! - kblk = 0 - pe = 0.0 - do k = km, kpblmin, -1 - if(kblk.eq.0 .and. k.le.komax(i)) then - pe = pe + bnv2(i,k)*(zl(i,komax(i))-zl(i,k))*del(i,k)/g/ro(i,k) - ke = 0.5*((rcs*u1(i,k))**2.+(rcs*v1(i,k))**2.) + kblk = 0 + pe = 0.0 + do k = km, kpblmin, -1 + if(kblk.eq.0 .and. k.le.komax(i)) then + pe = pe + bnv2(i,k)*(zl(i,komax(i))-zl(i,k))* & + del(i,k)/g/ro(i,k) + ke = 0.5*((rcs*u1(i,k))**2.+(rcs*v1(i,k))**2.) ! !---------- apply flow-blocking drag when pe >= ke ! - if(pe.ge.ke) then - kblk = k - kblk = min(kblk,kbl(i)) - zblk = zl(i,kblk)-zl(i,kts) - RDXZB(i) = real(k,kind=kind_phys) - endif - endif - enddo - if(kblk.ne.0) then + if(pe.ge.ke) then + kblk = k + kblk = min(kblk,kbl(i)) + zblk = zl(i,kblk)-zl(i,kts) + RDXZB(i) = real(k,kind=kind_phys) + endif + endif + enddo + if(kblk.ne.0) then ! !--------- compute flow-blocking stress ! - cd = max(2.0-1.0/od(i),0.0) - taufb(i,kts) = 0.5 * roll(i) * coefm(i) / max(dxmax_ls,dxy(i))**2 * cd * dxyp(i) & - * olp(i) * zblk * ulow(i)**2 - tautem = taufb(i,kts)/float(kblk-kts) - do k = kts+1, kblk - taufb(i,k) = taufb(i,k-1) - tautem - enddo + cd = max(2.0-1.0/od(i),0.0) + taufb(i,kts) = 0.5 * roll(i) * coefm(i) / & + max(dxmax_ls,dxy(i))**2 * cd * dxyp(i) * & + olp(i) * zblk * ulow(i)**2 + tautem = taufb(i,kts)/float(kblk-kts) + do k = kts+1, kblk + taufb(i,k) = taufb(i,k-1) - tautem + enddo ! !----------sum orographic GW stress and flow-blocking stress ! - ! taup(i,:) = taup(i,:) + taufb(i,:) ! Keep taup and taufb separate for now - endif - endif - enddo + ! taup(i,:) = taup(i,:) + taufb(i,:) ! Keep taup and taufb separate for now + endif + + endif ! if (.not.ldrag(i)) + + endif ! if ( ls_taper(i).GT.1.E-02 ) + + enddo ! do i=its,im -ENDIF ! end blocking drag +ENDIF ! IF ( (do_gsl_drag_ls_bl) .and. (gwd_opt_bl .EQ. 1) ) !=========================================================== IF ( (do_gsl_drag_ls_bl) .and. & - (gwd_opt_ls .EQ. 1 .OR. gwd_opt_bl .EQ. 1) .and. (ls_taper .GT. 1.E-02) ) THEN + (gwd_opt_ls .EQ. 1 .OR. gwd_opt_bl .EQ. 1) ) THEN + + do i=its,im + + if ( ls_taper(i) .GT. 1.E-02 ) then + ! ! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy ! - do k = kts,km - do i = its,im - taud_ls(i,k) = 1. * (taup(i,k+1) - taup(i,k)) * csg / del(i,k) - taud_bl(i,k) = 1. * (taufb(i,k+1) - taufb(i,k)) * csg / del(i,k) - enddo - enddo + do k = kts,km + taud_ls(i,k) = 1. * (taup(i,k+1) - taup(i,k)) * csg / del(i,k) + taud_bl(i,k) = 1. * (taufb(i,k+1) - taufb(i,k)) * csg / del(i,k) + enddo ! ! limit de-acceleration (momentum deposition ) at top to 1/2 value ! the idea is some stuff must go out the 'top' - do klcap = lcap,km - do i = its,im - taud_ls(i,klcap) = taud_ls(i,klcap) * factop - taud_bl(i,klcap) = taud_bl(i,klcap) * factop - enddo - enddo + do klcap = lcap,km + taud_ls(i,klcap) = taud_ls(i,klcap) * factop + taud_bl(i,klcap) = taud_bl(i,klcap) * factop + enddo ! ! if the gravity wave drag would force a critical line ! in the lower ksmm1 layers during the next deltim timestep, ! then only apply drag until that critical line is reached. ! - do k = kts,kpblmax-1 - do i = its,im - if (k .le. kbl(i)) then - if((taud_ls(i,k)+taud_bl(i,k)).ne.0.) & - dtfac(i) = min(dtfac(i),abs(velco(i,k) & - /(deltim*rcs*(taud_ls(i,k)+taud_bl(i,k))))) - endif - enddo - enddo + do k = kts,kpblmax-1 + if (k .le. kbl(i)) then + if ((taud_ls(i,k)+taud_bl(i,k)).ne.0.) & + dtfac(i) = min(dtfac(i),abs(velco(i,k) & + /(deltim*rcs*(taud_ls(i,k)+taud_bl(i,k))))) + endif + enddo ! - do k = kts,km - do i = its,im - taud_ls(i,k) = taud_ls(i,k) * dtfac(i) * ls_taper *(1.-rstoch(i)) - taud_bl(i,k) = taud_bl(i,k) * dtfac(i) * ls_taper *(1.-rstoch(i)) - - dtaux = taud_ls(i,k) * xn(i) - dtauy = taud_ls(i,k) * yn(i) - dtauxb = taud_bl(i,k) * xn(i) - dtauyb = taud_bl(i,k) * yn(i) - - !add blocking and large-scale contributions to tendencies - dudt(i,k) = dtaux + dtauxb + dudt(i,k) - dvdt(i,k) = dtauy + dtauyb + dvdt(i,k) - - if ( gsd_diss_ht_opt .EQ. 1 ) then - ! Calculate dissipation heating - ! Initial kinetic energy (at t0-dt) - eng0 = 0.5*( (rcs*u1(i,k))**2. + (rcs*v1(i,k))**2. ) - ! Kinetic energy after wave-breaking/flow-blocking - eng1 = 0.5*( (rcs*(u1(i,k)+(dtaux+dtauxb)*deltim))**2 + & - (rcs*(v1(i,k)+(dtauy+dtauyb)*deltim))**2 ) - ! Modify theta tendency - dtdt(i,k) = dtdt(i,k) + max((eng0-eng1),0.0)/cp/deltim - end if - - dusfc(i) = dusfc(i) + taud_ls(i,k)*xn(i)*del(i,k) + taud_bl(i,k)*xn(i)*del(i,k) - dvsfc(i) = dvsfc(i) + taud_ls(i,k)*yn(i)*del(i,k) + taud_bl(i,k)*yn(i)*del(i,k) - enddo - if(udtend>0) then - dtend(its:im,k,udtend) = dtend(its:im,k,udtend) + (taud_ls(its:im,k) * xn(its:im) + & - taud_bl(its:im,k) * xn(its:im)) * deltim - endif - if(vdtend>0) then - dtend(its:im,k,vdtend) = dtend(its:im,k,vdtend) + (taud_ls(its:im,k) * yn(its:im) + & - taud_bl(its:im,k) * yn(its:im)) * deltim - endif - if(gsd_diss_ht_opt .EQ. 1 .and. Tdtend>0) then - do i=its,im - ! Calculate dissipation heating - ! Initial kinetic energy (at t0-dt) - eng0 = 0.5*( (rcs*u1(i,k))**2. + (rcs*v1(i,k))**2. ) - ! Kinetic energy after wave-breaking/flow-blocking - eng1 = 0.5*( (rcs*(u1(i,k)+(dtaux+dtauxb)*deltim))**2 + & - (rcs*(v1(i,k)+(dtauy+dtauyb)*deltim))**2 ) - ! Modify theta tendency - dtend(i,k,Tdtend) = dtend(i,k,Tdtend) + max((eng0-eng1),0.0)/cp + do k = kts,km + taud_ls(i,k) = taud_ls(i,k)*dtfac(i)* ls_taper(i) *(1.-rstoch(i)) + taud_bl(i,k) = taud_bl(i,k)*dtfac(i)* ls_taper(i) *(1.-rstoch(i)) + + dtaux = taud_ls(i,k) * xn(i) + dtauy = taud_ls(i,k) * yn(i) + dtauxb = taud_bl(i,k) * xn(i) + dtauyb = taud_bl(i,k) * yn(i) + + !add blocking and large-scale contributions to tendencies + dudt(i,k) = dtaux + dtauxb + dudt(i,k) + dvdt(i,k) = dtauy + dtauyb + dvdt(i,k) + + if ( gsd_diss_ht_opt .EQ. 1 ) then + ! Calculate dissipation heating + ! Initial kinetic energy (at t0-dt) + eng0 = 0.5*( (rcs*u1(i,k))**2. + (rcs*v1(i,k))**2. ) + ! Kinetic energy after wave-breaking/flow-blocking + eng1 = 0.5*( (rcs*(u1(i,k)+(dtaux+dtauxb)*deltim))**2 + & + (rcs*(v1(i,k)+(dtauy+dtauyb)*deltim))**2 ) + ! Modify theta tendency + dtdt(i,k) = dtdt(i,k) + max((eng0-eng1),0.0)/cp/deltim + if ( Tdtend>0 ) then + dtend(i,k,Tdtend) = dtend(i,k,Tdtend) + max((eng0-eng1),0.0)/cp + endif + endif + + dusfc(i) = dusfc(i) + taud_ls(i,k)*xn(i)*del(i,k) + & + taud_bl(i,k)*xn(i)*del(i,k) + dvsfc(i) = dvsfc(i) + taud_ls(i,k)*yn(i)*del(i,k) + & + taud_bl(i,k)*yn(i)*del(i,k) + if(udtend>0) then + dtend(i,k,udtend) = dtend(i,k,udtend) + (taud_ls(i,k) * & + xn(i) + taud_bl(i,k) * xn(i)) * deltim + endif + if(vdtend>0) then + dtend(i,k,vdtend) = dtend(i,k,vdtend) + (taud_ls(i,k) * & + yn(i) + taud_bl(i,k) * yn(i)) * deltim + endif + enddo - endif - enddo - ! Finalize dusfc and dvsfc diagnostics - do i = its,im - dusfc(i) = (-1./g*rcs) * dusfc(i) - dvsfc(i) = (-1./g*rcs) * dvsfc(i) - enddo + ! Finalize dusfc and dvsfc diagnostics + dusfc(i) = (-1./g*rcs) * dusfc(i) + dvsfc(i) = (-1./g*rcs) * dvsfc(i) - if ( (gwd_opt == 33).or.(gwd_opt == 22) ) then - do k = kts,km - do i = its,im - dtaux2d_ls(i,k) = taud_ls(i,k) * xn(i) - dtauy2d_ls(i,k) = taud_ls(i,k) * yn(i) - dtaux2d_bl(i,k) = taud_bl(i,k) * xn(i) - dtauy2d_bl(i,k) = taud_bl(i,k) * yn(i) - dusfc_ls(i) = dusfc_ls(i) + dtaux2d_ls(i,k) * del(i,k) - dvsfc_ls(i) = dvsfc_ls(i) + dtauy2d_ls(i,k) * del(i,k) - dusfc_bl(i) = dusfc_bl(i) + dtaux2d_bl(i,k) * del(i,k) - dvsfc_bl(i) = dvsfc_bl(i) + dtauy2d_bl(i,k) * del(i,k) - enddo - enddo - endif + if ( (gwd_opt == 33).or.(gwd_opt == 22) ) then + do k = kts,km + dtaux2d_ls(i,k) = taud_ls(i,k) * xn(i) + dtauy2d_ls(i,k) = taud_ls(i,k) * yn(i) + dtaux2d_bl(i,k) = taud_bl(i,k) * xn(i) + dtauy2d_bl(i,k) = taud_bl(i,k) * yn(i) + dusfc_ls(i) = dusfc_ls(i) + dtaux2d_ls(i,k) * del(i,k) + dvsfc_ls(i) = dvsfc_ls(i) + dtauy2d_ls(i,k) * del(i,k) + dusfc_bl(i) = dusfc_bl(i) + dtaux2d_bl(i,k) * del(i,k) + dvsfc_bl(i) = dvsfc_bl(i) + dtauy2d_bl(i,k) * del(i,k) + enddo + endif + + endif ! if ( ls_taper(i) .GT. 1.E-02 ) + + enddo ! do i=its,im ENDIF ! (do_gsl_drag_ls_bl).and.(gwd_opt_ls.EQ.1 .OR. gwd_opt_bl.EQ.1) diff --git a/physics/drag_suite.meta b/physics/drag_suite.meta index 0cdcc682f..a36fc1e82 100644 --- a/physics/drag_suite.meta +++ b/physics/drag_suite.meta @@ -649,7 +649,7 @@ standard_name = cumulative_change_of_state_variables long_name = diagnostic tendencies for state variables units = various - dimensions = (horizontal_loop_extent,vertical_dimension,number_of_cumulative_change_processes) + dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_cumulative_change_processes) type = real kind = kind_phys active = (flag_for_diagnostics_3D) diff --git a/physics/gcycle.F90 b/physics/gcycle.F90 index 608147c4b..5f4f959c6 100644 --- a/physics/gcycle.F90 +++ b/physics/gcycle.F90 @@ -60,12 +60,9 @@ subroutine gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit, zorll(:), & zorlo(:), & weasd(:), & - slope(:), & snoalb(:), & canopy(:), & vfrac(:), & - vtype(:), & - stype(:), & shdmin(:), & shdmax(:), & snowd(:), & @@ -75,6 +72,9 @@ subroutine gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit, oro(:), & oro_uf(:), & slmsk(:) + integer, intent(inout) :: vtype(:), & + stype(:), & + slope(:) integer, intent(in) :: imap(:), jmap(:) ! @@ -84,6 +84,9 @@ subroutine gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit, real(kind=kind_io8) :: & slmskl (nx*ny), & slmskw (nx*ny), & + slpfcs (nx*ny), & + vegfcs (nx*ny), & + sltfcs (nx*ny), & TSFFCS (nx*ny), & ZORFCS (nx*ny), & AISFCS (nx*ny), & @@ -121,6 +124,10 @@ subroutine gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit, else TSFFCS = tsfco end if +! integer to real/double precision + slpfcs = real(slope) + vegfcs = real(vtype) + sltfcs = real(stype) ! if (frac_grid) then do ix=1,npts @@ -217,10 +224,10 @@ subroutine gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit, phour, xlat_d, xlon_d, slmskl, slmskw, & oro, oro_uf, use_ufo, nst_anl, & hice, fice, tisfc, snowd, slcfc1, & - shdmin, shdmax, slope, snoalb, tsffcs, & + shdmin, shdmax, slpfcs, snoalb, tsffcs, & weasd, zorfcs, albfc1, tg3, canopy, & smcfc1, stcfc1, slmsk, aisfcs, & - vfrac, vtype, stype, alffc1, cv, & + vfrac, vegfcs, sltfcs, alffc1, cv, & cvb, cvt, me, nthrds, & nlunit, size(input_nml_file), input_nml_file, & min_ice, ialb, isot, ivegsrc, & @@ -235,6 +242,11 @@ subroutine gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit, tsfc = TSFFCS tsfco = TSFFCS endif +! +! real/double precision to integer + slope = int(slpfcs) + vtype = int(vegfcs) + stype = int(sltfcs) ! do ix=1,npts zorll(ix) = ZORFCS(ix) diff --git a/physics/module_SGSCloud_RadPre.F90 b/physics/module_SGSCloud_RadPre.F90 index 20136ba00..73b72d10a 100644 --- a/physics/module_SGSCloud_RadPre.F90 +++ b/physics/module_SGSCloud_RadPre.F90 @@ -61,7 +61,7 @@ subroutine sgscloud_radpre_run( & ! should be moved to inside the mynn: use machine , only : kind_phys use module_radiation_clouds, only : gethml - use radcons, only: qmin ! Minimum vlaues for varius calculations + use radcons, only: qmin ! Minimum values for various calculations use funcphys, only: fpvs ! Function ot compute sat. vapor pressure over liq. !------------------------------------------------------------------- implicit none diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index c91c28dc4..4e376c6eb 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -439,7 +439,7 @@ MODULE module_mp_thompson !>\section gen_thompson_init thompson_init General Algorithm !> @{ SUBROUTINE thompson_init(is_aerosol_aware_in, & - merra2_aerosol_aware_in, & + merra2_aerosol_aware_in, & mpicomm, mpirank, mpiroot, & threads, errmsg, errflg) @@ -457,11 +457,11 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, & real :: stime, etime LOGICAL, PARAMETER :: precomputed_tables = .FALSE. -! Set module variable is_aerosol_aware +! Set module variable is_aerosol_aware/merra2_aerosol_aware is_aerosol_aware = is_aerosol_aware_in merra2_aerosol_aware = merra2_aerosol_aware_in if (mpirank==mpiroot) then - if (is_aerosol_aware) then + if (is_aerosol_aware) then write (0,'(a)') 'Using aerosol-aware version of Thompson microphysics' else if(merra2_aerosol_aware) then write (0,'(a)') 'Using merra2 aerosol-aware version of Thompson microphysics' @@ -1160,6 +1160,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & stop end if end if + if (is_aerosol_aware .and. (.not.present(nc) .or. & .not.present(nwfa) .or. & .not.present(nifa) .or. & @@ -1178,11 +1179,25 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & else stop end if - else if (.not.is_aerosol_aware .and. (present(nwfa) .or. & - present(nifa) .or. & - present(nwfa2d) .or. & - present(nifa2d) )) then - write(*,*) 'WARNING, nc/nwfa/nifa/nwfa2d/nifa2d present but is_aerosol_aware is FALSE' + else if (merra2_aerosol_aware .and. (.not.present(nc) .or. & + .not.present(nwfa) .or. & + .not.present(nifa) )) then + if (present(errmsg)) then + write(errmsg, '(*(a))') 'Logic error in mp_gt_driver: provide nc, nwfa, and nifa', & + ' for merra2 aerosol-aware version of Thompson microphysics' + else + write(*, '(*(a))') 'Logic error in mp_gt_driver: provide nc, nwfa, and nifa', & + ' for merra2 aerosol-aware version of Thompson microphysics' + end if + if (present(errflg)) then + errflg = 1 + return + else + stop + end if + else if (.not.is_aerosol_aware .and. .not.merra2_aerosol_aware .and. & + (present(nwfa) .or. present(nifa) .or. present(nwfa2d) .or. present(nifa2d))) then + write(*,*) 'WARNING, nc/nwfa/nifa/nwfa2d/nifa2d present but is_aerosol_aware/merra2_aerosol_aware are FALSE' end if end if test_only_once diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index cb4e21e52..48e02f36a 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -108,6 +108,7 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & errflg = 0 if (is_initialized) return + ! Consistency checks if (imp_physics/=imp_physics_thompson) then write(errmsg,'(*(a))') "Logic error: namelist choice of microphysics is different from Thompson MP" @@ -123,12 +124,17 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & end if end if + if (is_aerosol_aware .and. merra2_aerosol_aware) then + write(errmsg,'(*(a))') "Logic error: Only one Thompson aerosol option can be true, either is_aerosol_aware or merra2_aerosol_aware)" + errflg = 1 + return + end if + ! Call Thompson init - call thompson_init(is_aerosol_aware_in=is_aerosol_aware, & - merra2_aerosol_aware_in=merra2_aerosol_aware, & - mpicomm=mpicomm, & - mpirank=mpirank, mpiroot=mpiroot, threads=threads, & - errmsg=errmsg, errflg=errflg) + call thompson_init(is_aerosol_aware_in=is_aerosol_aware, & + merra2_aerosol_aware_in=merra2_aerosol_aware, & + mpicomm=mpicomm, mpirank=mpirank, mpiroot=mpiroot, & + threads=threads, errmsg=errmsg, errflg=errflg) if (errflg /= 0) return ! For restart runs, the init is done here @@ -148,6 +154,10 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & where(qs<0) qs = 0.0 where(qg<0) qg = 0.0 + if (merra2_aerosol_aware) then + call get_niwfa(aerfld, nifa, nwfa, ncol, nlev) + end if + !> - Convert specific humidity to water vapor mixing ratio. !> - Also, hydrometeor variables are mass or number mixing ratio !> - either kg of species per kg of dry air, or per kg of (dry + vapor). @@ -262,7 +272,6 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & ! Ensure we have 1st guess cloud droplet number where mass non-zero but no number. where(qc .LE. 0.0) nc=0.0 where(qc .GT. 0 .and. nc .LE. 0.0) nc = make_DropletNumber(qc*rho, nwfa*rho) * orho - where(qc .EQ. 0.0 .and. nc .GT. 0.0) nc = 0.0 ! Ensure non-negative aerosol number concentrations. where(nwfa .LE. 0.0) nwfa = 1.1E6 @@ -270,12 +279,12 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & ! Copy to local array for calculating cloud effective radii below nc_local = nc - else if(merra2_aerosol_aware) then - call get_niwfa(aerfld, nifa, nwfa, ncol, nlev) + + else if (merra2_aerosol_aware) then + ! Ensure we have 1st guess cloud droplet number where mass non-zero but no number. where(qc .LE. 0.0) nc=0.0 where(qc .GT. 0 .and. nc .LE. 0.0) nc = make_DropletNumber(qc*rho, nwfa*rho) * orho - where(qc .EQ. 0.0 .and. nc .GT. 0.0) nc = 0.0 else @@ -314,11 +323,11 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & end if if (convert_dry_rho) then - !qc = qc/(1.0_kind_phys+qv) - !qr = qr/(1.0_kind_phys+qv) - !qi = qi/(1.0_kind_phys+qv) - !qs = qs/(1.0_kind_phys+qv) - !qg = qg/(1.0_kind_phys+qv) + qc = qc/(1.0_kind_phys+qv) + qr = qr/(1.0_kind_phys+qv) + qi = qi/(1.0_kind_phys+qv) + qs = qs/(1.0_kind_phys+qv) + qg = qg/(1.0_kind_phys+qv) ni = ni/(1.0_kind_phys+qv) nr = nr/(1.0_kind_phys+qv) @@ -344,8 +353,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & con_eps, convert_dry_rho, & spechum, qc, qr, qi, qs, qg, ni, nr, & is_aerosol_aware, merra2_aerosol_aware, & - nc, nwfa, nifa, & - nwfa2d, nifa2d, & + nc, nwfa, nifa, nwfa2d, nifa2d, & tgrs, prsl, phii, omega, dt_inner, & dtp, first_time_step, istep, nsteps, & prcp, rain, graupel, ice, snow, sr, & @@ -513,7 +521,6 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & return end if - ! Set reduced time step if subcycling is used if (nsteps>1) then dtstep = dtp/real(nsteps, kind=kind_phys) @@ -531,15 +538,27 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & present(nifa) .and. & present(nwfa2d) .and. & present(nifa2d) )) then - write(errmsg,fmt='(*(a))') 'Logic error in mp_thompson_run:', & + write(errmsg,fmt='(*(a))') 'Logic error in mp_thompson_run:', & ' aerosol-aware microphysics require all of the', & ' following optional arguments:', & ' nc, nwfa, nifa, nwfa2d, nifa2d' errflg = 1 return + else if (is_aerosol_aware .and. .not. (present(nc) .and. & + present(nwfa) .and. & + present(nifa) )) then + write(errmsg,fmt='(*(a))') 'Logic error in mp_thompson_run:', & + ' merra2 aerosol-aware microphysics require the', & + ' following optional arguments: nc, nwfa, nifa' + errflg = 1 + return end if end if + if (merra2_aerosol_aware) then + call get_niwfa(aerfld, nifa, nwfa, ncol, nlev) + end if + !> - Convert specific humidity to water vapor mixing ratio. !> - Also, hydrometeor variables are mass or number mixing ratio !> - either kg of species per kg of dry air, or per kg of (dry + vapor). @@ -562,8 +581,6 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & nc = nc/(1.0_kind_phys-spechum) nwfa = nwfa/(1.0_kind_phys-spechum) nifa = nifa/(1.0_kind_phys-spechum) - else if (merra2_aerosol_aware) then - call get_niwfa(aerfld, nifa, nwfa, ncol, nlev) end if end if ! *DH @@ -747,8 +764,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte, & reset_dBZ=reset_dBZ, istep=istep, nsteps=nsteps, & - first_time_step=first_time_step, & - errmsg=errmsg, errflg=errflg, & + first_time_step=first_time_step, errmsg=errmsg, errflg=errflg, & ! Extended diagnostics ext_diag=ext_diag, & ! vts1=vts1, txri=txri, txrc=txrc, & @@ -827,8 +843,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte, & reset_dBZ=reset_dBZ, istep=istep, nsteps=nsteps, & - first_time_step=first_time_step, & - errmsg=errmsg, errflg=errflg, & + first_time_step=first_time_step, errmsg=errmsg, errflg=errflg, & ! Extended diagnostics ext_diag=ext_diag, & ! vts1=vts1, txri=txri, txrc=txrc, & @@ -966,9 +981,8 @@ subroutine get_niwfa(aerfld, nifa, nwfa, ncol, nlev) nifa=(aerfld(:,:,1)/4.0737762+aerfld(:,:,2)/30.459203+aerfld(:,:,3)/153.45048+ & aerfld(:,:,4)/1011.5142+ aerfld(:,:,5)/5683.3501)*1.e15 nwfa=((aerfld(:,:,6)/0.0045435214+aerfld(:,:,7)/0.2907854+aerfld(:,:,8)/12.91224+ & - aerfld(:,:,9)/206.2216+ aerfld(:,:,10)/4326.23)*1.+aerfld(:,:,11)/0.3053104*5+ & - +aerfld(:,:,15)/0.3232698*1)*1.e15 + aerfld(:,:,9)/206.2216+ aerfld(:,:,10)/4326.23)*1.+aerfld(:,:,11)/0.3053104*5+ & + aerfld(:,:,15)/0.3232698*1)*1.e15 end subroutine get_niwfa - end module mp_thompson diff --git a/physics/radiation_surface.f b/physics/radiation_surface.f index 750c54dd6..fed29526c 100644 --- a/physics/radiation_surface.f +++ b/physics/radiation_surface.f @@ -726,7 +726,7 @@ end subroutine setalb !! @{ !----------------------------------- subroutine setemis & - & ( lsm,lsm_noahmp,lsm_ruc,vtype,frac_grid, & ! --- inputs: + & ( lsm,lsm_noahmp,lsm_ruc,frac_grid, & ! --- inputs: & xlon,xlat,slmsk,snowf,sncovr,sncovr_ice, & ! & min_seaice,xlon,xlat,slmsk,snowf,sncovr,sncovr_ice, & & zorlf,tsknf,tairf,hprif, & @@ -784,7 +784,6 @@ subroutine setemis & integer, intent(in) :: IMAX integer, intent(in) :: lsm, lsm_noahmp, lsm_ruc logical, intent(in) :: frac_grid - real (kind=kind_phys), dimension(:), intent(in) :: vtype ! real (kind=kind_phys), intent(in) :: min_seaice real (kind=kind_phys), dimension(:), intent(in) :: & diff --git a/physics/rrtmg_lw_pre.F90 b/physics/rrtmg_lw_pre.F90 index 3ace48c0b..6da0e3100 100644 --- a/physics/rrtmg_lw_pre.F90 +++ b/physics/rrtmg_lw_pre.F90 @@ -1,6 +1,4 @@ !>\file rrtmg_lw_pre.f90 -!! This file contains a call to module_radiation_surface::setemis() to -!! setup surface emissivity for LW radiation. module rrtmg_lw_pre contains diff --git a/physics/rrtmg_sw_pre.F90 b/physics/rrtmg_sw_pre.F90 index cc329f180..1c7d3d76b 100644 --- a/physics/rrtmg_sw_pre.F90 +++ b/physics/rrtmg_sw_pre.F90 @@ -1,6 +1,4 @@ !>\file rrtmg_sw_pre.f90 -!! This file contains a subroutine to module_radiation_surface::setalb() to -!! setup surface albedo for SW radiation. module rrtmg_sw_pre contains diff --git a/physics/rrtmg_sw_pre.meta b/physics/rrtmg_sw_pre.meta index 9edf59e73..abf63a447 100644 --- a/physics/rrtmg_sw_pre.meta +++ b/physics/rrtmg_sw_pre.meta @@ -1,7 +1,7 @@ [ccpp-table-properties] name = rrtmg_sw_pre type = scheme - dependencies = iounitdef.f,machine.F,radiation_surface.f + dependencies = iounitdef.f,machine.F ######################################################################## [ccpp-arg-table] diff --git a/physics/rrtmgp_lw_pre.F90 b/physics/rrtmgp_lw_pre.F90 index 99318c1b8..d33a4e52c 100644 --- a/physics/rrtmgp_lw_pre.F90 +++ b/physics/rrtmgp_lw_pre.F90 @@ -1,8 +1,6 @@ module rrtmgp_lw_pre use machine, only: & kind_phys ! Working type - use module_radiation_surface, only: & - setemis ! Routine to compute surface-emissivity use mo_gas_optics_rrtmgp, only: & ty_gas_optics_rrtmgp use rrtmgp_lw_gas_optics, only: lw_gas_props diff --git a/physics/rrtmgp_lw_pre.meta b/physics/rrtmgp_lw_pre.meta index 071f1944a..3918f85e4 100644 --- a/physics/rrtmgp_lw_pre.meta +++ b/physics/rrtmgp_lw_pre.meta @@ -1,7 +1,7 @@ [ccpp-table-properties] name = rrtmgp_lw_pre type = scheme - dependencies = iounitdef.f,machine.F,physparam.f,radiation_surface.f + dependencies = iounitdef.f,machine.F,physparam.f ######################################################################## [ccpp-arg-table] diff --git a/physics/scm_sfc_flux_spec.F90 b/physics/scm_sfc_flux_spec.F90 index 398bff5d6..e4f425eb2 100644 --- a/physics/scm_sfc_flux_spec.F90 +++ b/physics/scm_sfc_flux_spec.F90 @@ -53,7 +53,7 @@ end subroutine scm_sfc_flux_spec_finalize !! -# Calculate the surface drag coefficient for heat and moisture. !! -# Calculate the u and v wind at 10m. subroutine scm_sfc_flux_spec_run (u1, v1, z1, t1, q1, p1, roughness_length, spec_sh_flux, spec_lh_flux, & - exner_inverse, T_surf, cp, grav, hvap, rd, fvirt, vonKarman, sh_flux, lh_flux, sh_flux_chs, lh_flux_chs, u_star, sfc_stress, cm, ch, & + exner_inverse, T_surf, cp, grav, hvap, rd, fvirt, vonKarman, sh_flux, lh_flux, sh_flux_chs, u_star, sfc_stress, cm, ch, & fm, fh, rb, u10m, v10m, wind1, qss, t2m, q2m, errmsg, errflg) use machine, only: kind_phys @@ -63,7 +63,7 @@ subroutine scm_sfc_flux_spec_run (u1, v1, z1, t1, q1, p1, roughness_length, spec real(kind=kind_phys), intent(in) :: cp, grav, hvap, rd, fvirt, vonKarman real(kind=kind_phys), intent(out) :: sh_flux(:), lh_flux(:), u_star(:), sfc_stress(:), & cm(:), ch(:), fm(:), fh(:), rb(:), u10m(:), v10m(:), wind1(:), qss(:), t2m(:), q2m(:), & - sh_flux_chs(:), lh_flux_chs(:) + sh_flux_chs(:) character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -83,8 +83,7 @@ subroutine scm_sfc_flux_spec_run (u1, v1, z1, t1, q1, p1, roughness_length, spec sh_flux(i) = spec_sh_flux(i) lh_flux(i) = spec_lh_flux(i) sh_flux_chs(i) = sh_flux(i) - lh_flux_chs(i) = lh_flux(i) - + roughness_length_m = 0.01*roughness_length(i) wind1(i) = sqrt(u1(i)*u1(i) + v1(i)*v1(i)) diff --git a/physics/scm_sfc_flux_spec.meta b/physics/scm_sfc_flux_spec.meta index 2899bbdcc..d7b29a4c0 100644 --- a/physics/scm_sfc_flux_spec.meta +++ b/physics/scm_sfc_flux_spec.meta @@ -182,8 +182,8 @@ intent = in optional = F [vonKarman] - standard_name = vonKarman_constant - long_name = vonKarman constant + standard_name = von_karman_constant + long_name = Von Karman constant units = none dimensions = () type = real @@ -209,23 +209,14 @@ intent = out optional = F [sh_flux_chs] - standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward sensible heat flux reduced by surface roughness + standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness_and_vegetation + long_name = kinematic surface upward sensible heat flux reduced by surface roughness and vegetation units = K m s-1 dimensions = (horizontal_loop_extent) type = real kind = kind_phys intent = out optional = F -[lh_flux_chs] - standard_name = surface_upward_specific_humidity_flux_reduced_by_surface_roughness - long_name = kinematic surface upward latent heat flux reduced by surface roughness - units = kg kg-1 m s-1 - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out - optional = F [u_star] standard_name = surface_friction_velocity long_name = boundary layer parameter diff --git a/physics/sfc_drv_ruc.F90 b/physics/sfc_drv_ruc.F90 index f313f2fba..e426424a8 100644 --- a/physics/sfc_drv_ruc.F90 +++ b/physics/sfc_drv_ruc.F90 @@ -62,8 +62,8 @@ subroutine lsm_ruc_init (me, master, isot, ivegsrc, nlunit, & real (kind=kind_phys), dimension(:), intent(in) :: slmsk - real (kind=kind_phys), dimension(:), intent(in) :: stype - real (kind=kind_phys), dimension(:), intent(in) :: vtype + integer, dimension(:), intent(in) :: stype + integer, dimension(:), intent(in) :: vtype real (kind=kind_phys), dimension(:), intent(in) :: q1 real (kind=kind_phys), dimension(:), intent(in) :: prsl1 real (kind=kind_phys), dimension(:), intent(in) :: tsfc_lnd @@ -112,12 +112,11 @@ subroutine lsm_ruc_init (me, master, isot, ivegsrc, nlunit, & real (kind=kind_phys) :: q0, qs1 integer :: ipr, i, k logical :: debug_print - integer, dimension(im) :: soiltyp, vegtype ! Initialize CCPP error handling variables errmsg = '' errflg = 0 - + ! Consistency checks if (lsm/=lsm_ruc) then write(errmsg,'(*(a))') 'Logic error: namelist choice of ', & @@ -164,31 +163,10 @@ subroutine lsm_ruc_init (me, master, isot, ivegsrc, nlunit, & pores (:) = maxsmc (:) resid (:) = drysmc (:) - soiltyp(:) = 0 - vegtype(:) = 0 - do i = 1, im ! i - horizontal loop - if (slmsk(i) == 2.) then - !-- ice - if (isot == 1) then - soiltyp(i) = 16 - else - soiltyp(i) = 9 - endif - if (ivegsrc == 1) then - vegtype(i) = 15 - elseif(ivegsrc == 2) then - vegtype(i) = 13 - endif - else - !-- land or water - soiltyp(i) = int( stype(i)+0.5 ) - vegtype(i) = int( vtype(i)+0.5 ) - if (soiltyp(i) < 1) soiltyp(i) = 14 - if (vegtype(i) < 1) vegtype(i) = 17 - endif + !-- initialize background emissivity - semisbase(i) = lemitbl(vegtype(i)) ! no snow effect + semisbase(i) = lemitbl(vtype(i)) ! no snow effect if (.not.flag_restart) then !-- land @@ -225,7 +203,7 @@ subroutine lsm_ruc_init (me, master, isot, ivegsrc, nlunit, & call rucinit (flag_restart, im, lsoil_ruc, lsoil, nlev, & ! in me, master, lsm_ruc, lsm, slmsk, & ! in - soiltyp, vegtype, & ! in + stype, vtype, & ! in tsfc_lnd, tsfc_wat, tg3, & ! in zs, dzs, smc, slc, stc, & ! in sh2o, smfrkeep, tslb, smois, & ! out @@ -280,8 +258,8 @@ end subroutine lsm_ruc_finalize ! ps - real, surface pressure (pa) im ! ! t1 - real, surface layer mean temperature (k) im ! ! q1 - real, surface layer mean specific humidity im ! -! soiltyp - integer, soil type (integer index) im ! -! vegtype - integer, vegetation type (integer index) im ! +! stype - integer, soil type (integer index) im ! +! vtype - integer, vegetation type (integer index) im ! ! sigmaf - real, areal fractional cover of green vegetation im ! ! dlwflx - real, total sky sfc downward lw flux ( w/m**2 ) im ! ! dswflx - real, total sky sfc downward sw flux ( w/m**2 ) im ! @@ -343,7 +321,7 @@ subroutine lsm_ruc_run & ! inputs & ( iter, me, master, delt, kdt, im, nlev, lsm_ruc, lsm, & & imp_physics, imp_physics_gfdl, imp_physics_thompson, & & do_mynnsfclay, lsoil_ruc, lsoil, rdlai, zs, & - & t1, q1, qc, soiltyp, vegtype, sigmaf, laixy, & + & t1, q1, qc, stype, vtype, sigmaf, laixy, & & dlwflx, dswsfc, tg3, coszen, land, icy, lake, & & rainnc, rainc, ice, snow, graupel, & & prsl1, zf, wind, shdmin, shdmax, & @@ -420,7 +398,8 @@ subroutine lsm_ruc_run & ! inputs logical, intent(in) :: rdlai ! --- in/out: - integer, dimension(:), intent(inout) :: soiltyp, vegtype + integer, dimension(:), intent(inout) :: stype + integer, dimension(:), intent(in) :: vtype real (kind=kind_phys), dimension(:), intent(in) :: zs real (kind=kind_phys), dimension(:), intent(in) :: srflag real (kind=kind_phys), dimension(:), intent(inout) :: & @@ -538,7 +517,7 @@ subroutine lsm_ruc_run & ! inputs ! local integer :: ims,ime, its,ite, jms,jme, jts,jte, kms,kme, kts,kte integer :: l, k, i, j, fractional_seaice, ilst - real (kind=kind_phys) :: dm, cimin + real (kind=kind_phys) :: dm, cimin(im) logical :: flag(im), flag_ice(im), flag_ice_uncoupled(im) logical :: rdlai2d, myj, frpcpn logical :: debug_print @@ -558,11 +537,11 @@ subroutine lsm_ruc_run & ! inputs if (icy(i) .and. .not. flag_cice(i)) then ! - uncoupled ice model if (oceanfrac(i) > zero) then - cimin = min_seaice + cimin(i) = min_seaice else - cimin = min_lakeice + cimin(i) = min_lakeice endif - if (fice(i) >= cimin) then + if (fice(i) >= cimin(i)) then ! - ice fraction is above the threshold for ice flag_ice(i) = .true. endif @@ -591,8 +570,8 @@ subroutine lsm_ruc_run & ! inputs if(debug_print) then write (0,*)'RUC LSM run' - write (0,*)'soiltyp=',ipr,soiltyp(ipr) - write (0,*)'vegtype=',ipr,vegtype(ipr) + write (0,*)'stype=',ipr,stype(ipr) + write (0,*)'vtype=',ipr,vtype(ipr) write (0,*)'kdt, iter =',kdt,iter write (0,*)'flag_init =',flag_init write (0,*)'flag_restart =',flag_restart @@ -655,8 +634,8 @@ subroutine lsm_ruc_run & ! inputs smcwlt2 (i) = 0. else !land - smcref2 (i) = REFSMC(soiltyp(i)) - smcwlt2 (i) = WLTSMC(soiltyp(i)) + smcref2 (i) = REFSMC(stype(i)) + smcwlt2 (i) = WLTSMC(stype(i)) endif enddo @@ -858,18 +837,18 @@ subroutine lsm_ruc_run & ! inputs tbot(i,j) = tg3(i) !> - 3. canopy/soil characteristics (s): -!!\n \a vegtyp - vegetation type (integer index) -> vtype -!!\n \a soiltyp - soil type (integer index) -> stype -!!\n \a sfcems - surface emmisivity -> sfcemis -!!\n \a sfalb_lnd_bck - backround snow-free surface albedo (fraction) -> albbck_lnd -!!\n \a snoalb - upper bound on maximum albedo over deep snow -> snoalb1d_lnd +!!\n \a vtype - vegetation type (integer index) +!!\n \a stype - soil type (integer index) +!!\n \a sfcems - surface emmisivity -> sfcemis +!!\n \a sfalb_lnd_bck - backround snow-free surface albedo (fraction) -> albbck_lnd +!!\n \a snoalb - upper bound on maximum albedo over deep snow -> snoalb1d_lnd if(ivegsrc == 1) then ! IGBP - MODIS vtype_wat(i,j) = 17 ! 17 - water (oceans and lakes) in MODIS stype_wat(i,j) = 14 xland_wat(i,j) = 2. ! xland = 2 for water - vtype_lnd(i,j) = vegtype(i) - stype_lnd(i,j) = soiltyp(i) + vtype_lnd(i,j) = vtype(i) + stype_lnd(i,j) = stype(i) vtype_ice(i,j) = 15 ! MODIS if(isot == 0) then stype_ice(i,j) = 9 ! ZOBLER @@ -936,7 +915,7 @@ subroutine lsm_ruc_run & ! inputs if(coszen(i) > 0. .and. weasd_lnd(i) < 1.e-4) then !-- solar zenith angle dependence when no snow - ilst=istwe(vegtype(i)) ! 1 or 2 + ilst=istwe(vtype(i)) ! 1 or 2 dm = (1.+2.*d(ilst))/(1.+2.*d(ilst)*coszen(i)) albbcksol(i) = sfalb_lnd_bck(i)*dm else @@ -1400,10 +1379,10 @@ subroutine lsm_ruc_run & ! inputs fwat = 1.0 - fice(i) ! Check if ice fraction is below the minimum value: 15% in GFS ! physics. - if (fice(i) < cimin) then ! cimin - minimal ice fraction + if (fice(i) < cimin(i)) then ! cimin - minimal ice fraction write (0,*)'warning: ice fraction is low:', fice(i) - fice(i) = cimin - fwat = 1.0 - cimin + fice(i) = cimin(i) + fwat = 1.0 - cimin(i) write (0,*)'fix ice fraction: reset it to:', fice(i), tskin_wat(i) endif @@ -1484,7 +1463,7 @@ end subroutine lsm_ruc_run !! This subroutine contains RUC LSM initialization. subroutine rucinit (restart, im, lsoil_ruc, lsoil, nlev, & ! in me, master, lsm_ruc, lsm, slmsk, & ! in - soiltyp, vegtype, & ! in + stype, vtype, & ! in tskin_lnd, tskin_wat, tg3, & ! !in zs, dzs, smc, slc, stc, & ! in sh2o, smfrkeep, tslb, smois, & ! out @@ -1507,8 +1486,7 @@ subroutine rucinit (restart, im, lsoil_ruc, lsoil, nlev, & ! in real (kind=kind_phys), dimension(im,lsoil), intent(in ) :: stc ! Noah real (kind=kind_phys), dimension(im,lsoil), intent(in ) :: slc ! Noah - integer, dimension(im), intent(inout) :: soiltyp - integer, dimension(im), intent(inout) :: vegtype + integer, dimension(im), intent(in) :: stype, vtype real (kind=kind_phys), dimension(im), intent(inout) :: wetness real (kind=kind_phys), dimension(im,lsoil_ruc), intent(inout) :: smois! ruc real (kind=kind_phys), dimension(im,lsoil_ruc), intent(inout) :: tslb ! ruc @@ -1639,13 +1617,13 @@ subroutine rucinit (restart, im, lsoil_ruc, lsoil, nlev, & ! in endif if(debug_print) then - write (0,*)'smc(ipr,:) ==', ipr, smc(ipr,:) - write (0,*)'stc(ipr,:) ==', ipr, stc(ipr,:) - write (0,*)'tskin_lnd(:)=',tskin_lnd(:) - write (0,*)'tskin_wat(:)=',tskin_wat(:) - write (0,*)'vegtype(ipr) ==', ipr, vegtype(ipr) - write (0,*)'soiltyp(ipr) ==', ipr, soiltyp(ipr) - write (0,*)'its,ite,jts,jte ',its,ite,jts,jte + write (0,*)'smc(ipr,:) =', ipr, smc(ipr,:) + write (0,*)'stc(ipr,:) =', ipr, stc(ipr,:) + write (0,*)'tskin_lnd(ipr) =', tskin_lnd(ipr) + write (0,*)'tskin_wat(ipr) =', tskin_wat(ipr) + write (0,*)'vtype(ipr) =', ipr, vtype(ipr) + write (0,*)'stype(ipr) =', ipr, stype(ipr) + write (0,*)'its,ite,jts,jte =',its,ite,jts,jte endif @@ -1654,8 +1632,8 @@ subroutine rucinit (restart, im, lsoil_ruc, lsoil, nlev, & ! in sst(i,j) = tskin_wat(i) tbot(i,j) = tg3(i) - ivgtyp(i,j) = vegtype(i) - isltyp(i,j) = soiltyp(i) + ivgtyp(i,j) = vtype(i) + isltyp(i,j) = stype(i) if (slmsk(i) == 0.) then !-- water tsk(i,j) = tskin_wat(i) @@ -1679,8 +1657,8 @@ subroutine rucinit (restart, im, lsoil_ruc, lsoil, nlev, & ! in !--- initialize smcwlt2 and smcref2 with Noah values if(slmsk(i) == 1.) then - smcref2 (i) = REFSMCnoah(soiltyp(i)) - smcwlt2 (i) = WLTSMCnoah(soiltyp(i)) + smcref2 (i) = REFSMCnoah(stype(i)) + smcwlt2 (i) = WLTSMCnoah(stype(i)) else smcref2 (i) = 1. smcwlt2 (i) = 0. diff --git a/physics/sfc_drv_ruc.meta b/physics/sfc_drv_ruc.meta index bdb058343..c793b5b9a 100644 --- a/physics/sfc_drv_ruc.meta +++ b/physics/sfc_drv_ruc.meta @@ -147,21 +147,19 @@ intent = in optional = F [stype] - standard_name = soil_type_classification_real + standard_name = soil_type_classification long_name = soil type for lsm units = index dimensions = (horizontal_dimension) - type = real - kind = kind_phys + type = integer intent = in optional = F [vtype] - standard_name = vegetation_type_classification_real + standard_name = vegetation_type_classification long_name = vegetation type for lsm units = index dimensions = (horizontal_dimension) - type = real - kind = kind_phys + type = integer intent = in optional = F [q1] @@ -758,21 +756,21 @@ kind = kind_phys intent = in optional = F -[soiltyp] +[stype] standard_name = soil_type_classification long_name = soil type at each grid cell units = index dimensions = (horizontal_loop_extent) type = integer - intent = inout + intent = in optional = F -[vegtype] +[vtype] standard_name = vegetation_type_classification long_name = vegetation type at each grid cell units = index dimensions = (horizontal_loop_extent) type = integer - intent = inout + intent = in optional = F [sigmaf] standard_name = vegetation_area_fraction diff --git a/physics/sfc_noah_wrfv4_interstitial.F90 b/physics/sfc_noah_wrfv4_interstitial.F90 index 7b37de568..ee66b0092 100644 --- a/physics/sfc_noah_wrfv4_interstitial.F90 +++ b/physics/sfc_noah_wrfv4_interstitial.F90 @@ -150,7 +150,8 @@ subroutine sfc_noah_wrfv4_pre_run (im, nsoil, ialb, isice, land, & logical, dimension(:), intent(in) :: flag_guess, flag_iter, land real(kind=kind_phys), dimension(:), intent(in) :: sfcprs, tprcp, sfctmp, q1, prslki, wind, cm, ch, snwdph - real(kind=kind_phys), dimension(:), intent(in) :: weasd, tsfc, vtype + real(kind=kind_phys), dimension(:), intent(in) :: weasd, tsfc + integer , dimension(:), intent(in) :: vtype real(kind=kind_phys), dimension(:,:), intent(in) :: smc, stc, slc logical, dimension(:), intent(inout) :: flag_lsm, flag_lsm_glacier @@ -179,7 +180,7 @@ subroutine sfc_noah_wrfv4_pre_run (im, nsoil, ialb, isice, land, & !from module_sf_noahdrv.F/lsminit if (.not. restart .and. first_time_step .and. ialb == 0) then do i = 1, im - snoalb(i) = maxalb(int(0.5 + vtype(i)))*0.01 + snoalb(i) = maxalb(vtype(i))*0.01 end do end if diff --git a/physics/sfc_noah_wrfv4_interstitial.meta b/physics/sfc_noah_wrfv4_interstitial.meta index bff028fdc..e2d98e15a 100644 --- a/physics/sfc_noah_wrfv4_interstitial.meta +++ b/physics/sfc_noah_wrfv4_interstitial.meta @@ -349,12 +349,11 @@ intent = in optional = F [vtype] - standard_name = vegetation_type_classification_real + standard_name = vegetation_type_classification long_name = vegetation type for lsm units = index dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys + type = integer intent = in optional = F [smc] diff --git a/physics/sfcsub.F b/physics/sfcsub.F index b0aefb858..a9c219052 100644 --- a/physics/sfcsub.F +++ b/physics/sfcsub.F @@ -2019,12 +2019,14 @@ subroutine sfccycle(lugb,len,lsoil,sig1t,deltsfc & ! *,' tsffcs=',tsffcs(iprnt),' slianl=',slianl(iprnt) do i=1,len + if (nint(slmskl(i)) /= 1) then if (sicanl(i) >= min_ice(i)) then slianl(i) = 2.0_kind_io8 else slianl(i) = zero sicanl(i) = zero endif + endif enddo if (fh-deltsfc > -0.001 ) then diff --git a/physics/ugwpv1_gsldrag.F90 b/physics/ugwpv1_gsldrag.F90 index 104fc8e3f..71193ed88 100644 --- a/physics/ugwpv1_gsldrag.F90 +++ b/physics/ugwpv1_gsldrag.F90 @@ -581,9 +581,9 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ntrac, lonr, dtp, fhzero,kd ! endif ! endif - else + endif ! -! not gsldrag oro-scheme for example "do_ugwp_v1_orog_only" +! not gsldrag large-scale oro-scheme for example "do_ugwp_v1_orog_only" ! if ( do_ugwp_v1_orog_only ) then @@ -641,7 +641,6 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ntrac, lonr, dtp, fhzero,kd dtend(:,:,idtend) = dtend(:,:,idtend) + Pdtdt*dtp endif endif - ENDIF ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Begin non-stationary GW schemes