From 59d82e9502d8e8c4ac7eab39950fe666b97404c3 Mon Sep 17 00:00:00 2001 From: Soren Rasmussen Date: Fri, 19 Jul 2024 09:31:41 -0600 Subject: [PATCH 1/9] Installing OpenMPI, fixing ENV and LABEL warnings, splitting some commands into multiple RUN statements for quicker debugging by making the failing step easier to find Installing ccmake in cmake-curses-gui fixes some CMake issues that were occurring, unsure why. Changing venv install location, old one wasn't working. Python packages needed to be updated. Switching to parallel builds. Make sure python package f90nml can be found in paths --- docker/Dockerfile | 69 ++++++++++++++++++++++------------------------- 1 file changed, 32 insertions(+), 37 deletions(-) diff --git a/docker/Dockerfile b/docker/Dockerfile index 7d0586090..85ad4d856 100644 --- a/docker/Dockerfile +++ b/docker/Dockerfile @@ -1,5 +1,5 @@ FROM debian:12 -MAINTAINER Michael Kavulich +LABEL maintainer="Michael Kavulich " # Set up base OS environment @@ -7,21 +7,25 @@ RUN apt-get -y update # Get "essential" tools and libraries RUN apt-get -y install build-essential \ - && apt-get -y install cmake curl git file gfortran-12 ksh m4 python3 tcsh time wget vim \ + && apt-get -y install cmake cmake-curses-gui curl git file gfortran-12 ksh m4 python3 tcsh time wget vim emacs-nox \ && apt-get -y install libnetcdf-pnetcdf-19 libnetcdff7 libnetcdf-dev libnetcdff-dev libxml2 \ - && apt-get -y install python3-pip python3.11-venv + && apt-get -y install python3-pip python3.11-venv python3-netcdf4 \ + && apt-get -y install openmpi-bin libopenmpi-dev +RUN ln -s /usr/bin/python3 /usr/bin/python -MAINTAINER Grant Firl or Michael Kavulich +# Set up python needed packages, preferred Docker method is apt-get but +# f90nml can't be installed for debian that way +RUN pip install f90nml --break-system-packages #Compiler environment variables -ENV CC /usr/bin/gcc -ENV FC /usr/bin/gfortran -ENV CXX /usr/bin/g++ -ENV F77 /usr/bin/gfortran -ENV F90 /usr/bin/gfortran +ENV CC=/usr/bin/gcc +ENV FC=/usr/bin/gfortran +ENV CXX=/usr/bin/g++ +ENV F77=/usr/bin/gfortran +ENV F90=/usr/bin/gfortran # Other necessary environment variables -ENV LD_LIBRARY_PATH /usr/lib/ +ENV LD_LIBRARY_PATH=/usr/lib/ # Set up unpriviledged user account, set up user home space and make sure user has permissions on all stuff in /comsoftware RUN groupadd comusers -g 9999 \ @@ -32,57 +36,54 @@ RUN groupadd comusers -g 9999 \ && chown -R comuser:comusers /comsoftware \ && chmod -R 6755 /comsoftware -# Link version-specific aliases (python3 will be created later with virtual environment) -RUN ln -s ~comuser/.venv/bin/python3 /usr/local/bin/python -RUN ln -s /usr/bin/gfortran-12 /usr/bin/gfortran - # all root steps completed above, now continue below as regular userID comuser USER comuser WORKDIR /home # Build NCEP libraries we need for SCM -ENV NCEPLIBS_DIR /comsoftware/nceplibs +ENV NCEPLIBS_DIR=/comsoftware/nceplibs RUN mkdir -p $NCEPLIBS_DIR/src && cd $NCEPLIBS_DIR/src \ && git clone -b v2.4.1 --recursive https://github.com/NOAA-EMC/NCEPLIBS-bacio \ && mkdir NCEPLIBS-bacio/build && cd NCEPLIBS-bacio/build \ && cmake -DCMAKE_INSTALL_PREFIX=$NCEPLIBS_DIR .. \ - && make VERBOSE=1 \ - && make install + && make VERBOSE=1 -j \ + && make install RUN cd $NCEPLIBS_DIR/src \ && git clone -b v2.3.3 --recursive https://github.com/NOAA-EMC/NCEPLIBS-sp \ && mkdir NCEPLIBS-sp/build && cd NCEPLIBS-sp/build \ && cmake -DCMAKE_INSTALL_PREFIX=$NCEPLIBS_DIR .. \ - && make VERBOSE=1 \ + && make VERBOSE=1 -j \ && make install RUN cd $NCEPLIBS_DIR/src \ && git clone -b v2.11.0 --recursive https://github.com/NOAA-EMC/NCEPLIBS-w3emc \ && mkdir NCEPLIBS-w3emc/build && cd NCEPLIBS-w3emc/build \ && cmake -DCMAKE_INSTALL_PREFIX=$NCEPLIBS_DIR .. \ - && make VERBOSE=1 \ + && make VERBOSE=1 -j \ && make install -ENV bacio_ROOT /comsoftware/nceplibs -ENV sp_ROOT /comsoftware/nceplibs -ENV w3emc_ROOT /comsoftware/nceplibs +ENV bacio_ROOT=/comsoftware/nceplibs +ENV sp_ROOT=/comsoftware/nceplibs +ENV w3emc_ROOT=/comsoftware/nceplibs -# Obtain CCPP SCM source code and static data, build code +# Obtain CCPP SCM source code, build code, and download static data RUN cd /comsoftware \ && git clone --recursive -b main https://github.com/NCAR/ccpp-scm \ - && cd /comsoftware/ccpp-scm/ \ - && ./contrib/get_all_static_data.sh \ - && ./contrib/get_thompson_tables.sh \ - && cd /comsoftware/ccpp-scm/scm \ - && mkdir bin \ - && cd bin \ + && mkdir /comsoftware/ccpp-scm/scm/bin + +RUN cd /comsoftware/ccpp-scm/scm/bin \ && cmake ../src \ - && make -j4 + && make -j + +RUN cd /comsoftware/ccpp-scm/ \ + && ./contrib/get_all_static_data.sh \ + && ./contrib/get_thompson_tables.sh # The analysis scripts have options for using LaTeX when making figure labels. -# If you would like to install LaTeK, uncomment the section below. +# If you would like to install LaTeK, uncomment the section below. # Note: This will increase the image size by 1 GB. #USER root #RUN yum -y update \ @@ -96,9 +97,3 @@ ENV SCM_ROOT=/comsoftware/ccpp-scm/ # For interactive use, vim mouse settings are infuriating RUN echo "set mouse=" > ~/.vimrc - -# Set up python virtual environment and install needed packages -ENV VIRTUAL_ENV=~/.venv -RUN python3 -m venv $VIRTUAL_ENV -ENV PATH="$VIRTUAL_ENV/bin:$PATH" -RUN pip3 install f90nml==1.4.4 netcdf4==1.6.5 From d2a41d3664cc0094a816f92e0d9cb96415686cc2 Mon Sep 17 00:00:00 2001 From: Soren Rasmussen Date: Wed, 24 Jul 2024 18:13:44 -0600 Subject: [PATCH 2/9] Changing Fortran stop statements to error stop statements. Adding github actions flag to run_scm.py so it will stop on first error --- .github/workflows/ci_test_docker.yml | 2 +- scm/src/CCPP_typedefs.F90 | 19 +- scm/src/GFS_typedefs.F90 | 134 +++++----- scm/src/run_scm.py | 136 +++++----- scm/src/scm.F90 | 86 +++---- scm/src/scm_forcing.F90 | 168 ++++++------- scm/src/scm_input.F90 | 357 +++++++++++++-------------- scm/src/scm_setup.F90 | 62 ++--- scm/src/scm_time_integration.F90 | 27 +- scm/src/scm_utils.F90 | 28 +-- scm/src/scm_vgrid.F90 | 46 ++-- 11 files changed, 534 insertions(+), 531 deletions(-) diff --git a/.github/workflows/ci_test_docker.yml b/.github/workflows/ci_test_docker.yml index e6508a00d..51991a10d 100644 --- a/.github/workflows/ci_test_docker.yml +++ b/.github/workflows/ci_test_docker.yml @@ -26,4 +26,4 @@ jobs: run: | mkdir $HOME/output chmod a+rw $HOME/output - docker run --rm -v $HOME/output:/home ${{ env.TEST_TAG }} ./run_scm.py -f ../../test/rt_test_cases.py --runtime_mult 0.1 -d + docker run --rm -v $HOME/output:/home ${{ env.TEST_TAG }} ./run_scm.py -f ../../test/rt_test_cases.py --runtime_mult 0.1 -d -a diff --git a/scm/src/CCPP_typedefs.F90 b/scm/src/CCPP_typedefs.F90 index 0b49ed010..4cfaeacd7 100644 --- a/scm/src/CCPP_typedefs.F90 +++ b/scm/src/CCPP_typedefs.F90 @@ -355,8 +355,8 @@ module CCPP_typedefs real (kind=kind_phys), pointer :: qs_lay(:,:) => null() !< real (kind=kind_phys), pointer :: q_lay(:,:) => null() !< real (kind=kind_phys), pointer :: deltaZ(:,:) => null() !< - real (kind=kind_phys), pointer :: deltaZc(:,:) => null() !< - real (kind=kind_phys), pointer :: deltaP(:,:) => null() !< + real (kind=kind_phys), pointer :: deltaZc(:,:) => null() !< + real (kind=kind_phys), pointer :: deltaP(:,:) => null() !< real (kind=kind_phys), pointer :: cloud_overlap_param(:,:) => null() !< Cloud overlap parameter real (kind=kind_phys), pointer :: cnv_cloud_overlap_param(:,:) => null() !< Convective cloud overlap parameter real (kind=kind_phys), pointer :: precip_overlap_param(:,:) => null() !< Precipitation overlap parameter @@ -373,12 +373,12 @@ module CCPP_typedefs real (kind=kind_phys), pointer :: cld_rwp(:,:) => null() !< Cloud rain water path real (kind=kind_phys), pointer :: cld_rerain(:,:) => null() !< Cloud rain effective radius real (kind=kind_phys), pointer :: precip_frac(:,:) => null() !< Precipitation fraction - real (kind=kind_phys), pointer :: cld_cnv_frac(:,:) => null() !< SGS convective cloud fraction + real (kind=kind_phys), pointer :: cld_cnv_frac(:,:) => null() !< SGS convective cloud fraction real (kind=kind_phys), pointer :: cld_cnv_lwp(:,:) => null() !< SGS convective cloud liquid water path real (kind=kind_phys), pointer :: cld_cnv_reliq(:,:) => null() !< SGS convective cloud liquid effective radius real (kind=kind_phys), pointer :: cld_cnv_iwp(:,:) => null() !< SGS convective cloud ice water path real (kind=kind_phys), pointer :: cld_cnv_reice(:,:) => null() !< SGS convective cloud ice effecive radius - real (kind=kind_phys), pointer :: cld_pbl_lwp(:,:) => null() !< SGS PBL cloud liquid water path + real (kind=kind_phys), pointer :: cld_pbl_lwp(:,:) => null() !< SGS PBL cloud liquid water path real (kind=kind_phys), pointer :: cld_pbl_reliq(:,:) => null() !< SGS PBL cloud liquid effective radius real (kind=kind_phys), pointer :: cld_pbl_iwp(:,:) => null() !< SGS PBL cloud ice water path real (kind=kind_phys), pointer :: cld_pbl_reice(:,:) => null() !< SGS PBL cloud ice effecive radius @@ -390,7 +390,7 @@ module CCPP_typedefs real (kind=kind_phys), pointer :: fluxswDOWN_allsky(:,:) => null() !< RRTMGP downward shortwave all-sky flux profile real (kind=kind_phys), pointer :: fluxswUP_clrsky(:,:) => null() !< RRTMGP upward shortwave clr-sky flux profile real (kind=kind_phys), pointer :: fluxswDOWN_clrsky(:,:) => null() !< RRTMGP downward shortwave clr-sky flux profile - real (kind=kind_phys), pointer :: sfc_emiss_byband(:,:) => null() !< + real (kind=kind_phys), pointer :: sfc_emiss_byband(:,:) => null() !< real (kind=kind_phys), pointer :: sec_diff_byband(:,:) => null() !< real (kind=kind_phys), pointer :: sfc_alb_nir_dir(:,:) => null() !< real (kind=kind_phys), pointer :: sfc_alb_nir_dif(:,:) => null() !< @@ -888,7 +888,7 @@ subroutine gfs_interstitial_setup_tracers(Interstitial, Model) Interstitial%nvdiff = Interstitial%nvdiff + 1 ENDIF if (Model%me == Model%master) write(0,*) 'nssl_settings2: nvdiff,ntrac = ', Interstitial%nvdiff, Model%ntrac - + elseif (Model%imp_physics == Model%imp_physics_wsm6) then Interstitial%nvdiff = Model%ntrac -3 if (Model%satmedmf) Interstitial%nvdiff = Interstitial%nvdiff + 1 @@ -965,8 +965,7 @@ subroutine gfs_interstitial_setup_tracers(Interstitial, Model) Interstitial%nvdiff = 9 endif else - write(0,*) "Selected microphysics scheme is not supported when coupling with chemistry" - stop + error stop "Selected microphysics scheme is not supported when coupling with chemistry" endif if (Interstitial%trans_aero) Interstitial%nvdiff = Interstitial%nvdiff + Model%ntchm if (Model%ntke > 0) Interstitial%nvdiff = Interstitial%nvdiff + 1 ! adding tke to the list @@ -1419,11 +1418,11 @@ subroutine gfs_interstitial_phys_reset (Interstitial, Model) ! Use same logic in UFS to reset Thompson extended diagnostics Interstitial%ext_diag_thompson_reset = Interstitial%max_hourly_reset ! - ! Frequency flag for computing the full radar reflectivity (water coated ice) + ! Frequency flag for computing the full radar reflectivity (water coated ice) if (Model%nsfullradar_diag<0) then Interstitial%fullradar_diag = .true. else - Interstitial%fullradar_diag = (Model%kdt == 1 .or. mod(Model%kdt, nint(Model%nsfullradar_diag/Model%dtp)) == 0) + Interstitial%fullradar_diag = (Model%kdt == 1 .or. mod(Model%kdt, nint(Model%nsfullradar_diag/Model%dtp)) == 0) end if ! diff --git a/scm/src/GFS_typedefs.F90 b/scm/src/GFS_typedefs.F90 index d8fe5f446..35f325f03 100644 --- a/scm/src/GFS_typedefs.F90 +++ b/scm/src/GFS_typedefs.F90 @@ -2,7 +2,7 @@ module GFS_typedefs use mpi_f08 use machine, only: kind_phys, kind_dbl_prec, kind_sngl_prec - + use module_radsw_parameters, only: topfsw_type, sfcfsw_type use module_radlw_parameters, only: topflw_type, sfcflw_type use h2o_def, only: levh2o, h2o_coeff @@ -541,14 +541,14 @@ module GFS_typedefs real (kind=kind_phys), pointer :: evap_lnd(:) => null() !< sfc latent heat flux over land, converted to evaporative flux real (kind=kind_phys), pointer :: hflx_lnd(:) => null() !< sfc sensible heat flux over land real (kind=kind_phys), pointer :: ep_lnd(:) => null() !< sfc up pot latent heat flux over land - real (kind=kind_phys), pointer :: t2mmp_lnd(:) => null() !< 2 meter temperature over land + real (kind=kind_phys), pointer :: t2mmp_lnd(:) => null() !< 2 meter temperature over land real (kind=kind_phys), pointer :: q2mp_lnd(:) => null() !< 2 meter spec humidity over land real (kind=kind_phys), pointer :: gflux_lnd(:) => null() !< soil heat flux over land real (kind=kind_phys), pointer :: runoff_lnd(:) => null() !< surface runoff over land real (kind=kind_phys), pointer :: drain_lnd(:) => null() !< subsurface runoff over land real (kind=kind_phys), pointer :: cmm_lnd(:) => null() !< surface drag wind speed for momentum - real (kind=kind_phys), pointer :: chh_lnd(:) => null() !< surface drag mass flux for heat and moisture - real (kind=kind_phys), pointer :: zvfun_lnd(:) => null() !< function of surface roughness length and green vegetation fraction + real (kind=kind_phys), pointer :: chh_lnd(:) => null() !< surface drag mass flux for heat and moisture + real (kind=kind_phys), pointer :: zvfun_lnd(:) => null() !< function of surface roughness length and green vegetation fraction !--- outgoing accumulated quantities real (kind=kind_phys), pointer :: rain_cpl (:) => null() !< total rain precipitation @@ -754,7 +754,7 @@ module GFS_typedefs logical :: cplaqm !< default no cplaqm collection logical :: cplchm !< default no cplchm collection logical :: cpllnd !< default no cpllnd collection - logical :: cpllnd2atm !< default no lnd->atm coupling + logical :: cpllnd2atm !< default no lnd->atm coupling logical :: rrfs_sd !< default no rrfs_sd collection logical :: use_cice_alb !< default .false. - i.e. don't use albedo imported from the ice model logical :: cpl_imp_mrg !< default no merge import with internal forcings @@ -1225,8 +1225,8 @@ module GFS_typedefs real(kind=kind_phys) :: rbcr !< Critical Richardson Number in the PBL scheme real(kind=kind_phys) :: betascu !< Tuning parameter for prog. closure shallow clouds - real(kind=kind_phys) :: betamcu !< Tuning parameter for prog. closure midlevel clouds - real(kind=kind_phys) :: betadcu !< Tuning parameter for prog. closure deep clouds + real(kind=kind_phys) :: betamcu !< Tuning parameter for prog. closure midlevel clouds + real(kind=kind_phys) :: betadcu !< Tuning parameter for prog. closure deep clouds !--- MYNN parameters/switches logical :: do_mynnedmf @@ -3075,7 +3075,7 @@ subroutine coupling_create (Coupling, IM, Model) Coupling%slmsk_cpl = clear_val !< pointer to sfcprop%slmsk endif - ! -- Coupling options to retrive land fluxes from external land component + ! -- Coupling options to retrive land fluxes from external land component if (Model%cpllnd .and. Model%cpllnd2atm) then allocate (Coupling%sncovr1_lnd (IM)) allocate (Coupling%qsurf_lnd (IM)) @@ -4003,7 +4003,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- coupling parameters cplflx, cplice, cplocn2atm, cplwav, cplwav2atm, cplaqm, & cplchm, cpllnd, cpllnd2atm, cpl_imp_mrg, cpl_imp_dbg, & - rrfs_sd, use_cice_alb, & + rrfs_sd, use_cice_alb, & #ifdef IDEA_PHYS lsidea, weimer_model, f107_kp_size, f107_kp_interval, & f107_kp_skip_size, f107_kp_data_size, f107_kp_read_in_start, & @@ -4192,7 +4192,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & inquire (file=trim(fn_nml), exist=exists) if (.not. exists) then write(6,*) 'GFS_namelist_read:: namelist file: ',trim(fn_nml),' does not exist' - stop + error stop else open (unit=nlunit, file=fn_nml, action='READ', status='OLD', iostat=ios) endif @@ -4228,7 +4228,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%qdiag3d = qdiag3d if (qdiag3d .and. .not. ldiag3d) then write(0,*) 'Logic error in GFS_typedefs.F90: qdiag3d requires ldiag3d' - stop + error stop endif Model%flag_for_gwd_generic_tend = .true. Model%flag_for_pbl_generic_tend = .true. @@ -4255,12 +4255,12 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & else write(0,*) 'CCPP suite simulator turned on, but error encountered loading data.' write(0,*) errmsg - stop + error stop endif if(.not. qdiag3d .and. .not. ldiag3d) then write(0,*) 'CCPP suite simulator turned on, but qdiag3d and/or ldiag3d are not set to .true.' write(0,*) errmsg - stop + error stop endif endif @@ -4340,11 +4340,11 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & ! if (naux2d>naux2dmax) then write(0,*) "Error, number of requested auxiliary 2d arrays exceeds the maximum defined in GFS_typedefs.F90" - stop + error stop endif if (naux3d>naux3dmax) then write(0,*) "Error, number of requested auxiliary 3d arrays exceeds the maximum defined in GFS_typedefs.F90" - stop + error stop endif Model%naux2d = naux2d Model%naux3d = naux3d @@ -4358,7 +4358,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & end if if (any(aux2d_time_avg) .or. any(aux3d_time_avg)) then write(0,*) "Error, the SCM has not implemented time averaging of diagnostics in GFS_typedefs.F90" - stop + error stop end if Model%fhcyc = fhcyc @@ -4397,7 +4397,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & ! Model%cplflx == .true. and Model%cplice == .false. (HAFS FV3ATM-HYCOM) if (Model%cplice .and. .not. Model%cplflx) then print *,' Logic error: Model%cplflx==.false. and Model%cplice==.true. is currently not supported - shutting down' - stop + error stop endif Model%cplocn2atm = cplocn2atm Model%cplwav = cplwav @@ -4451,7 +4451,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%lsidea = lsidea if (Model%lsidea) then print *,' LSIDEA is active but needs to be reworked for FV3 - shutting down' - stop + error stop endif #ifdef IDEA_PHYS !--- integrated dynamics through earth's atmosphere @@ -4490,7 +4490,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%levr = levs else if (levr > levs) then write(0,*) "Logic error, number of radiation levels (levr) cannot exceed number of model levels (levs)" - stop + error stop else Model%levr = levr endif @@ -4498,11 +4498,11 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & if (isubc_sw < 0 .or. isubc_sw > 2) then write(0,'(a,i0)') 'ERROR: shortwave cloud-sampling (isubc_sw) scheme selected not valid: ',isubc_sw - stop + error stop endif if (isubc_lw < 0 .or. isubc_lw > 2) then write(0,'(a,i0)') 'ERROR: longwave cloud-sampling (isubc_lw) scheme selected not valid: ',isubc_lw - stop + error stop endif @@ -4510,7 +4510,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & (iovr .ne. Model%iovr_max) .and. (iovr .ne. Model%iovr_dcorr) .and. & (iovr .ne. Model%iovr_exp) .and. (iovr .ne. Model%iovr_exprand)) then write(0,'(a,i0)') 'ERROR: cloud-overlap (iovr) scheme selected not valid: ',iovr - stop + error stop endif if ((isubc_sw == 0 .or. isubc_lw == 0) .and. iovr > 2 ) then @@ -4615,16 +4615,16 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & ! RRTMGP incompatible with levr /= levs if (Model%levr /= Model%levs) then write(0,*) "Logic error, RRTMGP only works with levr = levs" - stop + error stop end if ! RRTMGP LW scattering calculation not supported w/ RRTMG cloud-optics if (Model%doGP_lwscat .and. Model%doG_cldoptics) then write(0,*) "Logic error, RRTMGP Longwave cloud-scattering not supported with RRTMG cloud-optics." - stop + error stop end if if (Model%doGP_sgs_mynn .and. .not. do_mynnedmf) then write(0,*) "Logic error, RRTMGP flag doGP_sgs_mynn only works with do_mynnedmf=.true." - stop + error stop endif if (Model%doGP_sgs_cnv .or. Model%doGP_sgs_mynn) then write(0,*) "RRTMGP explicit cloud scheme being used." @@ -4635,7 +4635,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & if (Model%doGP_cldoptics_PADE .and. Model%doGP_cldoptics_LUT) then write(0,*) "Logic error, Both RRTMGP cloud-optics options cannot be selected. " - stop + error stop end if if (.not. Model%doGP_cldoptics_PADE .and. .not. Model%doGP_cldoptics_LUT .and. .not. Model%doG_cldoptics) then write(0,*) "Logic error, No option for cloud-optics scheme provided. Using RRTMG cloud-optics" @@ -4644,7 +4644,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & if (Model%rrtmgp_nGptsSW .lt. 0 .or. Model%rrtmgp_nGptsLW .lt. 0 .or. & Model%rrtmgp_nBandsSW .lt. 0 .or. Model%rrtmgp_nBandsLW .lt. 0) then write(0,*) "Logic error, RRTMGP spectral dimensions (bands/gpts) need to be provided." - stop + error stop endif else if (Model%use_LW_jacobian) then @@ -4660,7 +4660,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & if (.not.lwhtr .or. .not.swhtr) then write(0,*) "Logic error, the CCPP version of RRTMG lwrad/swrad require the output" // & " of the lw/sw heating rates to be turned on (namelist options lwhtr and swhtr)" - stop + error stop end if !--- microphysical switch @@ -4723,7 +4723,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%mraerosol = mraerosol if (Model%ltaerosol .and. Model%mraerosol) then write(0,*) 'Logic error: Only one Thompson aerosol option can be true, either ltaerosol or mraerosol)' - stop + error stop end if Model%lradar = lradar Model%nsfullradar_diag = nsfullradar_diag @@ -4754,7 +4754,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%rdlai = rdlai if (Model%rdlai .and. .not. Model%lsm == Model%lsm_ruc) then write(0,*) 'Logic error: rdlai = .true. only works with RUC LSM' - stop + error stop end if ! Set surface layers for CCPP physics @@ -4772,7 +4772,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & if (Model%lsm==Model%lsm_noah .or. Model%lsm==Model%lsm_noahmp) then if (Model%lsoil_lsm/=4) then write(0,*) 'Error in GFS_typedefs.F90, number of soil layers must be 4 for Noah/NoahMP' - stop + error stop end if Model%zs = (/-0.1_kind_phys, -0.4_kind_phys, -1.0_kind_phys, -2.0_kind_phys/) Model%dzs = (/ 0.1_kind_phys, 0.3_kind_phys, 0.6_kind_phys, 1.0_kind_phys/) @@ -4785,7 +4785,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & if (Model%lsm==Model%lsm_ruc) then if (Model%lsoil_lsm/=9) then write(0,*) 'Error in GFS_typedefs.F90, number of soil layers must be 9 for RUC' - stop + error stop end if end if @@ -4795,12 +4795,12 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & if (Model%lsm==Model%lsm_noah .or. Model%lsm==Model%lsm_noahmp) then if (kice/=2) then write(0,*) 'Error in GFS_typedefs.F90, number of ice model layers must be 2 for Noah/NoahMP/Noah_WRFv4' - stop + error stop end if elseif (Model%lsm==Model%lsm_ruc) then if (kice/=9) then write(0,*) 'Error in GFS_typedefs.F90, number of ice model layers must be 9 for RUC' - stop + error stop end if end if @@ -4813,7 +4813,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & if (Model%lsm==Model%lsm_noahmp) then if (lsnow_lsm/=3) then write(0,*) 'Logic error: NoahMP expects the maximum number of snow layers to be exactly 3 (see sfc_noahmp_drv.f)' - stop + error stop else Model%lsnow_lsm = lsnow_lsm ! Set lower bound for LSM model, runs from negative (above surface) to surface (zero) @@ -4843,7 +4843,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !see GFS_MP_generic_post.F90; exticeden is only compatible with GFDL, !Thompson, or NSSL MP print *,' Using exticeden = T is only valid when using GFDL, Thompson, or NSSL microphysics.' - stop + error stop end if ! GFDL surface layer options Model%lcurr_sf = lcurr_sf @@ -4917,19 +4917,19 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !HWRF physics suite if (hwrf_samfdeep .and. imfdeepcnv/=2) then write(*,*) 'Logic error: hwrf_samfdeep requires imfdeepcnv=2' - stop + error stop end if if (hwrf_samfshal .and. imfshalcnv/=2) then write(*,*) 'Logic error: hwrf_samfshal requires imfshalcnv=2' - stop + error stop end if Model%hwrf_samfdeep = hwrf_samfdeep Model%hwrf_samfshal = hwrf_samfshal - !--prognostic closure - moisture coupling + !--prognostic closure - moisture coupling if ((progsigma .and. imfdeepcnv/=2) .and. (progsigma .and. imfdeepcnv/=5)) then write(*,*) 'Logic error: progsigma requires imfdeepcnv=2 or 5' - stop + error stop end if Model%progsigma = progsigma Model%betascu = betascu @@ -4938,7 +4938,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & if (oz_phys .and. oz_phys_2015) then write(*,*) 'Logic error: can only use one ozone physics option (oz_phys or oz_phys_2015), not both. Exiting.' - stop + error stop end if Model%oz_phys = oz_phys Model%oz_phys_2015 = oz_phys_2015 @@ -5556,7 +5556,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & write(0,*) 'NSSL micro: error! CCN is OFF (nssl_ccn_on = F) but ntccn > 1.' write(0,*) 'Should either remove ccn_nc from field_table or set nssl_ccn_on = .true.' ENDIF - stop + error stop ENDIF Model%ntccn = -99 Model%ntccna = -99 @@ -5565,7 +5565,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & write(*,*) 'NSSL micro: error! CCN is ON but ntccn < 1. Must have ccn_nc in field_table if nssl_ccn_on=T' write(0,*) 'NSSL micro: error! CCN is ON but ntccn < 1. Must have ccn_nc in field_table if nssl_ccn_on=T' ENDIF - stop + error stop ELSE if (Model%me == Model%master) then write(*,*) 'NSSL micro: CCN is ON' @@ -5601,7 +5601,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & write(0,*) 'missing needed tracers for NSSL hail! nthl > 1 but either volume or number is not in field_table' write(0,*) 'nthv, nthnc = ', Model%nthv, Model%nthnc ENDIF - stop + error stop ENDIF ENDIF @@ -5620,21 +5620,21 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & write(*,*) 'NSSL micro: CCN is ON' ENDIF ENDIF - + ! add checks for nssl_3moment - IF ( ( Model%nssl_3moment ) ) THEN + IF ( ( Model%nssl_3moment ) ) THEN IF ( Model%ntrz < 1 ) THEN write(*,*) 'NSSL micro: 3-moment is ON, but rain_ref tracer is missing' - stop + error stop ENDIF IF ( Model%ntgz < 1 ) THEN write(*,*) 'NSSL micro: 3-moment is ON, but graupel_ref tracer is missing' - stop + error stop ENDIF IF ( nssl_hail_on ) THEN IF ( Model%nthz < 1 ) THEN write(*,*) 'NSSL micro: 3-moment is ON, but hail_ref tracer is missing' - stop + error stop ENDIF ENDIF ENDIF @@ -5646,7 +5646,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%ntcw < 1 .or. Model%ntlnc < 1 & ) THEN if (Model%me == Model%master) write(0,*) 'missing needed tracers for NSSL!' - stop + error stop ENDIF @@ -5762,7 +5762,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & if (do_shoc) then if (Model%imp_physics == Model%imp_physics_thompson) then print *,'SHOC is not currently compatible with Thompson MP -- shutting down' - stop + error stop endif Model%nshoc_3d = 3 Model%nshoc_2d = 0 @@ -5781,7 +5781,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & if (Model%do_mynnedmf) then if (Model%do_shoc .or. Model%hybedmf .or. Model%satmedmf) then print *,' Logic error: MYNN EDMF cannot be run with SHOC, HEDMF or SATMEDMF' - stop + error stop end if ! Model%shal_cnv = .false. ! Model%imfshalcnv = -1 @@ -5827,14 +5827,14 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & elseif (Model%lsm == 0) then print *,' OSU no longer supported - job aborted' - stop + error stop elseif (Model%lsm == Model%lsm_noahmp) then if (Model%ivegsrc /= 1) then print *,'Vegetation type must be IGBP if Noah MP is used' - stop + error stop elseif (Model%isot /= 1) then print *,'Soil type must be STATSGO if Noah MP is used' - stop + error stop endif print *, 'New Noah MP Land Surface Model will be used' print *, 'The Physics options are' @@ -5863,7 +5863,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & print *,' add_fire_heat_flux = ',add_fire_heat_flux else print *,' Unsupported LSM type - job aborted - lsm=',Model%lsm - stop + error stop endif ! if (Model%lsm == Model%lsm_noahmp .and. Model%iopt_snf == 4) then @@ -5908,7 +5908,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & if ( (Model%imfdeepcnv == Model%imfdeepcnv_ntiedtke .or. Model%imfshalcnv == Model%imfshalcnv_ntiedtke) .and. & .not. (Model%imfdeepcnv == Model%imfdeepcnv_ntiedtke .and. Model%imfshalcnv == Model%imfshalcnv_ntiedtke) ) then write(0,*) "Logic error: if NTDK deep convection is used, must also use NTDK shallow convection (and vice versa)" - stop + error stop end if if (.not. Model%cscnv) then @@ -6048,7 +6048,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%npsdelt = 2 if (nwat /= 2) then print *,' Zhao-Carr MP requires nwat to be set to 2 - job aborted' - stop + error stop end if if (Model%me == Model%master) print *,' Using Zhao/Carr/Sundqvist Microphysics' @@ -6069,7 +6069,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%nseffr = 3 if (nwat /= 4) then print *,' Ferrier-Aligo MP requires nwat to be set to 4 - job aborted' - stop + error stop end if if (Model%me == Model%master) print *,' Using Ferrier-Aligo MP scheme', & ' microphysics', & @@ -6077,7 +6077,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & elseif (Model%imp_physics == Model%imp_physics_wsm6) then !WSM6 microphysics print *,' Error, WSM6 no longer supported - job aborted' - stop + error stop !Model%npdf3d = 0 !Model%num_p3d = 3 !Model%num_p2d = 1 @@ -6101,7 +6101,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & ENDIF if ( nwat /= 6+i ) then print *,' NSSL MP requires nwat to be set to ', 6+i,' - job aborted, nssl_hail_on = ',nssl_hail_on - stop + error stop end if Model%nleffr = 1 Model%nieffr = 2 @@ -6132,11 +6132,11 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%nseffr = 3 if (nwat /= 6) then print *,' Thompson MP requires nwat to be set to 6 - job aborted' - stop + error stop end if if (.not. Model%effr_in) then print *,' Thompson MP requires effr_in to be set to .true. - job aborted' - stop + error stop end if if (Model%me == Model%master) print *,' Using Thompson double moment microphysics', & ' ltaerosol = ',Model%ltaerosol, & @@ -6168,7 +6168,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & endif if (nwat /= 6 .and. Model%fprcp >= 2) then print *,' Morrison-Gettelman MP requires nwat to be set to 6 - job aborted' - stop + error stop end if if (Model%me == Model%master) & print *,' Using Morrison-Gettelman double moment microphysics', & @@ -6208,14 +6208,14 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%shcnvcw = .false. if (nwat /= 6) then print *,' GFDL MP requires nwat to be set to 6 - job aborted' - stop + error stop end if if (Model%me == Model%master) print *,' avg_max_length=',Model%avg_max_length if (Model%me == Model%master) print *,' Using GFDL Cloud Microphysics' else if (Model%me == Model%master) print *,'Wrong imp_physics value. Job abort.' - stop + error stop endif if(Model%ras .or. Model%cscnv) Model%cnvcld = .false. @@ -6349,7 +6349,7 @@ subroutine control_initialize_radar_tten(Model, radar_tten_limits) write(0,'(A,F12.4,A)') 'radar_tten_limits(1) = ',radar_tten_limits(1),' <-- lower limit' write(0,'(A,F12.4,A)') 'radar_tten_limits(2) = ',radar_tten_limits(2),' <-- upper limit' write(0,*) "If you do not want me to apply the prescribed tendencies, just say so! Remove fh_dfi_radar from your namelist." - stop + error stop endif else !o! Rejoice !o! Radar_tten_limits had lower and upper bounds. diff --git a/scm/src/run_scm.py b/scm/src/run_scm.py index b499fafac..35222f168 100755 --- a/scm/src/run_scm.py +++ b/scm/src/run_scm.py @@ -120,6 +120,7 @@ parser.add_argument('--n_itt_out', help='period of instantaneous output (number of timesteps)', required=False, type=int) parser.add_argument('--n_itt_diag', help='period of diagnostic output (number of timesteps)', required=False, type=int) parser.add_argument('-dt', '--timestep', help='timestep (s)', required=False, type=float) +parser.add_argument('-a', '--actions', help='if running from Github Actions this will select correct configuration', required=False, action='store_true') parser.add_argument('-v', '--verbose', help='set logging level to debug and write log to file', action='count', default=0) parser.add_argument('-f', '--file', help='name of file where SCM runs are defined') parser.add_argument('--mpi_command', help='command used to invoke the executable via MPI (including options)', required=False) @@ -190,13 +191,14 @@ def parse_arguments(): bin_dir = args.bin_dir timestep = args.timestep mpi_command = args.mpi_command - + github_actions = args.actions + if not sdf: sdf = DEFAULT_SUITE - + return (file, case, sdf, namelist, tracers, gdb, runtime, runtime_mult, docker, \ verbose, levels, npz_type, vert_coord_file, case_data_dir, n_itt_out, \ - n_itt_diag, run_dir, bin_dir, timestep, mpi_command) + n_itt_diag, run_dir, bin_dir, timestep, mpi_command, github_actions) def find_gdb(): """Detect gdb, abort if not found""" @@ -217,28 +219,28 @@ def __init__(self, case, suite, runtime, runtime_mult, levels, npz_type, vert_co """Initialize experiment. This routine does most of the work, including setting and checking the experiment configuration (namelist).""" - + self._case = case self._suite_obj = suite self._suite = suite._name self._name = case + '_' + self._suite - + self._physics_namelist = suite.namelist - + #check to see that the physics namelists exists in the right dir if not os.path.isfile(os.path.join(SCM_ROOT, PHYSICS_NAMELIST_DIR, self._physics_namelist)): message = 'The physics namelist {0} was not found'.format(os.path.join(SCM_ROOT, PHYSICS_NAMELIST_DIR, self._physics_namelist)) logging.critical(message) raise Exception(message) - + self._tracers = suite.tracers - + #check to see that the tracers exists in the right dir if not os.path.isfile(os.path.join(SCM_ROOT, TRACERS_DIR, self._tracers)): message = 'The tracer configuration {0} was not found'.format(os.path.join(TRACERS_DIR, self._tracers)) logging.critical(message) raise Exception(message) - + #check to see if the case namelists exists in the right dir self._namelist = os.path.join(SCM_ROOT, CASE_NAMELIST_DIR, self._case + '.nml') if not os.path.isfile(self._namelist): @@ -254,14 +256,14 @@ def __init__(self, case, suite, runtime, runtime_mult, levels, npz_type, vert_co self._runtime = None message = 'Namelist runtime adjustment {0} IS NOT applied'.format(self._runtime) logging.debug(message) - + if runtime_mult: self._runtime_mult = runtime_mult message = 'Existing case namelist runtime multiplied by {0}'.format(self._runtime_mult) logging.debug(message) else: self._runtime_mult = None - + if levels: self._levels = levels message = 'The number of vertical levels is set to {0}'.format(self._levels) @@ -270,7 +272,7 @@ def __init__(self, case, suite, runtime, runtime_mult, levels, npz_type, vert_co self._levels = None message = 'The number of vertical levels contained in the case configuration file is used if present, otherwise the default value in scm_input.F90 is used.' logging.debug(message) - + if npz_type: self._npz_type = npz_type message = 'The npz_type of vertical levels is set to {0}'.format(self._npz_type) @@ -294,7 +296,7 @@ def __init__(self, case, suite, runtime, runtime_mult, levels, npz_type, vert_co self._vert_coord_file = None message = 'The npz_type contained in the case configuration file is used if present, otherwise the default value in scm_input.F90 is used.' logging.debug(message) - + if case_data_dir: self._case_data_dir = case_data_dir else: @@ -304,12 +306,12 @@ def __init__(self, case, suite, runtime, runtime_mult, levels, npz_type, vert_co message = 'The case data directory {0} was not found'.format(self._case_data_dir) logging.critical(message) raise Exception(message) - + if n_itt_out: self._n_itt_out = n_itt_out else: self._n_itt_out = DEFAULT_OUTPUT_PERIOD - + if n_itt_diag: self._n_itt_diag = n_itt_diag else: @@ -319,7 +321,7 @@ def __init__(self, case, suite, runtime, runtime_mult, levels, npz_type, vert_co self._timestep = timestep else: self._timestep = suite.timestep - + @property def name(self): """Get the name of the experiment.""" @@ -339,27 +341,27 @@ def namelist(self): def namelist(self, value): """Set the case namelist of the experiment.""" self._namelist = value - + @property def case(self): """Get the case of the experiment.""" return self._case - + @name.setter def case(self, value): """Set the case of the experiment.""" self._case = value - + @property def suite(self): """Get the suite of the experiment.""" return self._suite - + @suite.setter def suite(self, value): """Set the suite of the experiment.""" self._suite = value - + @property def physics_namelist(self): """Get the physics namelist of the experiment.""" @@ -369,7 +371,7 @@ def physics_namelist(self): def physics_namelist(self, value): """Set the physics namelist of the experiment.""" self._physics_namelist = value - + @property def tracers(self): """Get the tracer file for the experiment.""" @@ -379,7 +381,7 @@ def tracers(self): def tracers(self, value): """Set the tracer file for the experiment.""" self._tracers = value - + @property def levels(self): """Get the number of vertical levels for the experiment.""" @@ -389,7 +391,7 @@ def levels(self): def levels(self, value): """Set the number of vertical levels for the experiment.""" self._levels = value - + @property def npz_type(self): """Get the vertical level type for the experiment.""" @@ -399,7 +401,7 @@ def npz_type(self): def npz_type(self, value): """Set the vertical level type for the experiment.""" self._npz_type = value - + @property def vert_coord_file(self): """Get the file containing vertical levels for the experiment.""" @@ -409,7 +411,7 @@ def vert_coord_file(self): def vert_coord_file(self, value): """Set the file containing vertical levels for the experiment.""" self._vert_coord_file = value - + @property def case_data_dir(self): """Get the case data directory for the experiment.""" @@ -419,7 +421,7 @@ def case_data_dir(self): def case_data_dir(self, value): """Set the case data directory for the experiment.""" self._case_data_dir = value - + @property def n_itt_out(self): """Get the output period (in timesteps) for the experiment.""" @@ -429,7 +431,7 @@ def n_itt_out(self): def n_itt_out(self, value): """Set the output period (in timesteps) for the experiment.""" self._n_itt_out = value - + @property def n_itt_diag(self): """Get the diagnostic period (in timesteps) for the experiment.""" @@ -439,10 +441,10 @@ def n_itt_diag(self): def n_itt_diag(self, value): """Set the diagnostic period (in timesteps) for the experiment.""" self._n_itt_diag = value - + def setup_rundir(self): """Set up run directory for this experiment.""" - + # Parse case configuration namelist and extract # - output directory # - surface_flux_spec @@ -475,7 +477,7 @@ def setup_rundir(self): case_nml['case_config']['n_itt_diag'] = self._n_itt_diag if self._timestep: case_nml['case_config']['dt'] = self._timestep - # look for the output_dir variable in the case configuration namelist and use it if it does; + # look for the output_dir variable in the case configuration namelist and use it if it does; # if it doesn't exist, create a default output directory name (from the case and suite names) and create a namelist patch try: output_dir = case_nml['case_config']['output_dir'] @@ -488,7 +490,7 @@ def setup_rundir(self): output_dir = 'output_' + self._case + '_' + self._suite + '_' + os.path.splitext(self._physics_namelist)[0] output_dir_patch_nml = {'case_config':{'output_dir':output_dir}} custom_output_dir = False - + #if using the DEPHY format, need to also check the case data file for the surfaceForcing global attribute for 'Flux' or 'surfaceFlux', which denotes prescribed surface fluxes try: input_type = case_nml['case_config']['input_type'] @@ -506,34 +508,34 @@ def setup_rundir(self): surface_flux_spec = case_nml['case_config']['sfc_flux_spec'] except KeyError: surface_flux_spec = False - + # If surface fluxes are specified for this case, use the SDF modified to use them if surface_flux_spec: logging.info('Specified surface fluxes are used for case {0}. Switching to SDF {1} from {2}'.format(self._case,'suite_' + self._suite + '_ps' + '.xml','suite_' + self._suite + '.xml')) self._suite = self._suite + '_ps' - + # Create physics_config namelist for experiment configuration file physics_config = {"physics_suite":self._suite, "physics_nml":self._physics_namelist,} physics_config_dict = {"physics_config":physics_config} physics_config_nml = f90nml.namelist.Namelist(physics_config_dict) - + # Create STANDARD_EXPERIMENT_NAMELIST in the run directory with the case configuration and physics configuration namelists logging.info('Creating experiment configuration namelist {0} in the run directory from {1} using {2} and {3} '.format(STANDARD_EXPERIMENT_NAMELIST,self._namelist,self._suite,self._physics_namelist)) - + with open(os.path.join(SCM_RUN, STANDARD_EXPERIMENT_NAMELIST), "w+") as nml_file: case_nml.write(nml_file) - + with open(os.path.join(SCM_RUN, STANDARD_EXPERIMENT_NAMELIST), "a") as nml_file: physics_config_nml.write(nml_file) - + # if using the default output dir name created in this script, patch the experiment namelist with the new output_dir variable if(not custom_output_dir): # GJF TODO: this implementation is clunky; newer versions of f90nml can handle this better, but this works with v0.19 so no need to require newer version f90nml.patch(os.path.join(SCM_RUN, STANDARD_EXPERIMENT_NAMELIST), output_dir_patch_nml, 'temp.nml') cmd = "mv {0} {1}".format('temp.nml', os.path.join(SCM_RUN, STANDARD_EXPERIMENT_NAMELIST)) execute(cmd) - + # Link physics namelist to run directory with its original name logging.debug('Linking physics namelist {0} to run directory'.format(self._physics_namelist)) if os.path.isfile(os.path.join(SCM_RUN, self._physics_namelist)): @@ -550,7 +552,7 @@ def setup_rundir(self): logging.info('Found optional prescribed surface physics namelist {0}; linking it to run directory'.format(opt_ps_nml_filename)) cmd = "ln -sf {0} {1}".format(opt_ps_nml_filename, os.path.join(SCM_RUN, self._physics_namelist)) execute(cmd) - + # Link tracer configuration to run directory with standard name logging.debug('Linking tracer configuration {0} to run directory'.format(self._tracers)) if os.path.isfile(os.path.join(SCM_RUN, self._tracers)): @@ -561,7 +563,7 @@ def setup_rundir(self): raise Exception(message) cmd = "ln -sf {0} {1}".format(os.path.join(SCM_ROOT, TRACERS_DIR, self._tracers), os.path.join(SCM_RUN, TRACERS_LINK)) execute(cmd) - + # Link case data file to run directory with original name try: input_type = case_nml['case_config']['input_type'] @@ -580,7 +582,7 @@ def setup_rundir(self): raise Exception(message) cmd = "ln -sf {0} {1}".format(os.path.join(SCM_ROOT, self._case_data_dir, case_data_netcdf_file), os.path.join(SCM_RUN, case_data_netcdf_file)) execute(cmd) - + # Link vertical coordinate file to run directory with its original name if (self._npz_type == 'input'): logging.debug('Linking vertical coordinate file {0} to run directory'.format(self._vert_coord_file)) @@ -592,7 +594,7 @@ def setup_rundir(self): raise Exception(message) cmd = "ln -sf {0} {1}".format(os.path.join(SCM_ROOT, VERT_COORD_DATA_DIR, self._vert_coord_file), os.path.join(SCM_RUN, self._vert_coord_file)) execute(cmd) - + # Link physics SDF to run directory physics_suite = 'suite_' + self._suite + '.xml' logging.debug('Linking physics suite {0} to run directory'.format(physics_suite)) @@ -604,7 +606,7 @@ def setup_rundir(self): raise Exception(message) cmd = "ln -sf {0} {1}".format(os.path.join(SCM_ROOT, PHYSICS_SUITE_DIR, physics_suite), os.path.join(SCM_RUN, physics_suite)) execute(cmd) - + # Link physics data needed for schemes to run directory logging.debug('Linking physics input data from {0} into run directory'.format(os.path.join(SCM_ROOT, PHYSICS_DATA_DIR))) for entry in os.listdir(os.path.join(SCM_ROOT, PHYSICS_DATA_DIR)): @@ -613,7 +615,7 @@ def setup_rundir(self): logging.debug('Linking file {0}'.format(entry)) cmd = 'ln -sf {0} {1}'.format(os.path.join(SCM_ROOT, PHYSICS_DATA_DIR, entry), os.path.join(SCM_RUN, entry)) execute(cmd) - + # Link reference profile data to run directory logging.debug('Linking reference profile data from {0} into run directory'.format(os.path.join(SCM_ROOT, REFERENCE_PROFILE_DIR))) for entry in REFERENCE_PROFILE_FILE_LIST: @@ -622,7 +624,7 @@ def setup_rundir(self): logging.debug('Linking file {0}'.format(entry)) cmd = 'ln -sf {0} {1}'.format(os.path.join(SCM_ROOT, REFERENCE_PROFILE_DIR, entry), os.path.join(SCM_RUN, entry)) execute(cmd) - + # Parse physics namelist and extract # - oz_phys # - oz_phys_2015 @@ -643,7 +645,7 @@ def setup_rundir(self): message = 'Logic error, both oz_phys and oz_phys_2015 are set to true in the physics namelist' logging.critical(message) raise Exception(message) - + # Link input data for oz_phys or oz_phys_2015 if os.path.exists(os.path.join(SCM_RUN, OZ_PHYS_LINK)): os.remove(os.path.join(SCM_RUN, OZ_PHYS_LINK)) @@ -655,13 +657,13 @@ def setup_rundir(self): logging.debug('Linking input data for oz_phys_2015') cmd = 'ln -sf {0} {1}'.format(os.path.join(SCM_RUN, OZ_PHYS_2015_TARGET), os.path.join(SCM_RUN, OZ_PHYS_LINK)) execute(cmd) - + # Look for do_ugwp_v1 try: do_ugwp_v1 = nml['gfs_physics_nml']['do_ugwp_v1'] except KeyError: do_ugwp_v1 = DEFAULT_DO_UGWP_V1 - + # Link the tau file if do_ugwp_v1 if do_ugwp_v1: if os.path.exists(os.path.join(SCM_RUN, TAU_LINK)): @@ -669,7 +671,7 @@ def setup_rundir(self): logging.debug('Linking input data for UGWP_v1') cmd = 'ln -sf {0} {1}'.format(os.path.join(SCM_RUN, TAU_TARGET), os.path.join(SCM_RUN, TAU_LINK)) execute(cmd) - + # Link scripts needed to run SCM analysis logging.debug('Linking analysis scripts from {0} into run directory'.format(os.path.join(SCM_ROOT, SCM_ANALYSIS_SCRIPT_DIR))) analysis_script_files = ['scm_analysis.py','configspec.ini'] @@ -679,7 +681,7 @@ def setup_rundir(self): logging.debug('Linking file {0}'.format(entry)) cmd = 'ln -sf {0} {1}'.format(os.path.join(SCM_ROOT, SCM_ANALYSIS_SCRIPT_DIR, entry), os.path.join(SCM_RUN, entry)) execute(cmd) - + # Link plot configuration files needed to run SCM analysis logging.debug('Linking plot configuration files from {0} into run directory'.format(os.path.join(SCM_ROOT, SCM_ANALYSIS_CONFIG_DIR))) for entry in os.listdir(os.path.join(SCM_ROOT, SCM_ANALYSIS_CONFIG_DIR)): @@ -688,18 +690,18 @@ def setup_rundir(self): logging.debug('Linking file {0}'.format(entry)) cmd = 'ln -sf {0} {1}'.format(os.path.join(SCM_ROOT, SCM_ANALYSIS_CONFIG_DIR, entry), os.path.join(SCM_RUN, entry)) execute(cmd) - + # Create output directory (delete existing directory) logging.debug('Creating output directory {0} in run directory'.format(output_dir)) if os.path.isdir(os.path.join(SCM_RUN, output_dir)): shutil.rmtree(os.path.join(SCM_RUN, output_dir)) os.makedirs(os.path.join(SCM_RUN, output_dir)) - + # Write experiment configuration file to output directory logging.debug('Writing experiment configuration {0}.nml to output directory'.format(self._name)) cmd = 'cp {0} {1}'.format(os.path.join(SCM_RUN, STANDARD_EXPERIMENT_NAMELIST), os.path.join(SCM_RUN, output_dir,self._name + '.nml')) execute(cmd) - + # Move executable to run dir if COPY_EXECUTABLE: logging.debug('Copying executable to run directory') @@ -709,7 +711,7 @@ def setup_rundir(self): logging.debug('Linking executable to run directory') cmd = 'ln -sf {0} {1}'.format(os.path.join(SCM_ROOT, SCM_BIN, EXECUTABLE_NAME), os.path.join(SCM_RUN, EXECUTABLE_NAME)) execute(cmd) - + #Inform user of timestep and output intervals if self._timestep: logging.info('Using {0}s as the timestep with an instantaneous output period of {1}s and a diagnostic output period of {2}s'.format( @@ -717,7 +719,7 @@ def setup_rundir(self): else: logging.info('Using the default timestep in src/scm_input.F90 with an instantaneous output period of {0}*dt and a diagnostic output period of {1}*dt'.format( case_nml['case_config']['n_itt_out'],case_nml['case_config']['n_itt_diag'])) - + return os.path.join(SCM_RUN, output_dir) def launch_executable(use_gdb, gdb, mpi_command, ignore_error = False): @@ -753,7 +755,7 @@ def launch_executable(use_gdb, gdb, mpi_command, ignore_error = False): else: time_elapsed = int(minutes)*60 + float(seconds) return (status, time_elapsed) - + def copy_outdir(exp_dir): """Copy output directory to /home for this experiment.""" dir_name = os.path.basename(exp_dir) @@ -765,8 +767,8 @@ def copy_outdir(exp_dir): def main(): (file, case, sdf, namelist, tracers, use_gdb, runtime, runtime_mult, docker, \ verbose, levels, npz_type, vert_coord_file, case_data_dir, n_itt_out, \ - n_itt_diag, run_dir, bin_dir, timestep, mpi_command) = parse_arguments() - + n_itt_diag, run_dir, bin_dir, timestep, mpi_command, github_actions) = parse_arguments() + setup_logging(verbose) global SCM_ROOT @@ -782,7 +784,7 @@ def main(): SCM_BIN = bin_dir else: SCM_BIN = os.path.join(SCM_ROOT, DEFAULT_BIN_DIR) - + global SCM_RUN if run_dir: SCM_RUN = run_dir @@ -790,10 +792,10 @@ def main(): SCM_RUN = os.path.join(SCM_ROOT, DEFAULT_RUN_DIR) if not os.path.isdir(SCM_RUN): os.makedirs(SCM_RUN) - + global EXECUTABLE EXECUTABLE = os.path.join(SCM_RUN, EXECUTABLE_NAME) - + # Debugger if use_gdb: gdb = find_gdb() @@ -844,9 +846,9 @@ def main(): logging.info('Using provided namelist for suite={0}'.format(run["suite"])) except: logging.info('Using default namelist for suite={0}'.format(run["suite"])) - + # If tracer file provided, use. Otherwise use default tracer file for this suite - try: + try: active_suite.tracers = run["tracers"] logging.info('Using provided tracer file for suite={0}'.format(run["suite"])) except: @@ -859,7 +861,7 @@ def main(): # If namelist and tracer file provided, use. Otherwise error. if ("namelist" in run) and ("tracers" in run): irun = irun + 1 - if timestep: + if timestep: active_suite = suite(run["suite"], run["tracers"], run["namelist"], timestep, -1, False) else: active_suite = suite(run["suite"], run["tracers"], run["namelist"], -1, -1, False) @@ -882,6 +884,10 @@ def main(): l_ignore_error = MULTIRUN_IGNORE_ERROR else: l_ignore_error = False + # need to correctly fail if running Github Actions + if github_actions: + l_ignore_error = False + (status, time_elapsed) = launch_executable(use_gdb, gdb, mpi_command, ignore_error = l_ignore_error) # if status == 0: diff --git a/scm/src/scm.F90 b/scm/src/scm.F90 index ff28593f0..b6ac80d4a 100644 --- a/scm/src/scm.F90 +++ b/scm/src/scm.F90 @@ -45,12 +45,12 @@ subroutine scm_main_sub() call MPI_INIT(ierr) if (ierr/=0) then write(*,*) 'An error occurred in MPI_INIT: ', ierr - stop 1 + error stop end if - fcst_mpi_comm = MPI_COMM_WORLD + fcst_mpi_comm = MPI_COMM_WORLD call get_config_nml(scm_state) - + select case(scm_state%input_type) case(0) call get_case_init(scm_state, scm_input_instance) @@ -58,9 +58,9 @@ subroutine scm_main_sub() call get_case_init_DEPHY(scm_state, scm_input_instance) case default write(*,*) 'An unrecognized specification of the input_type namelist variable is being used. Exiting...' - stop + error stop end select - + call get_reference_profile(scm_state, scm_reference) call get_FV3_vgrid(scm_input_instance, scm_state) @@ -75,15 +75,15 @@ subroutine scm_main_sub() scm_state%itt = 1 in_spinup = (scm_state%do_spinup .and. scm_state%itt <= scm_state%spinup_timesteps) - + if (in_spinup) then call set_spinup_nudging(scm_state) end if - + call interpolate_forcing(scm_input_instance, scm_state, in_spinup) call physics%create(scm_state%n_cols) - + !physics initialization section !set the array index of the time level of the state variables that the cdata @@ -97,13 +97,13 @@ subroutine scm_main_sub() case default cdata_time_index = 2 end select - + !open a logfile if (physics%Init_parm%me == physics%Init_parm%master .and. physics%Init_parm%logunit>=0) then write (logfile_name, '(A7,I0.5,A4)') 'logfile.out' open(unit=physics%Init_parm%logunit, file=trim(scm_state%output_dir)//'/'//logfile_name, action='write', status='replace') end if - + physics%Init_parm%fcst_mpi_comm = fcst_mpi_comm physics%Init_parm%levs = scm_state%n_levels physics%Init_parm%bdat(1) = scm_state%init_year @@ -126,7 +126,7 @@ subroutine scm_main_sub() physics%Init_parm%hydrostatic = .true. physics%Init_parm%restart = .false. physics%Init_parm%nwat = scm_state%nwat - + ! Allocate and initialize DDTs call GFS_suite_setup(physics%Model, physics%Statein, physics%Stateout, & physics%Sfcprop, physics%Coupling, physics%Grid, & @@ -134,7 +134,7 @@ subroutine scm_main_sub() physics%Diag, physics%Interstitial, 1, 1, & physics%Init_parm, scm_state%n_cols, scm_state%lon, & scm_state%lat, scm_state%area) - + !override radiation frequency if (scm_state%force_rad_T > 0) then !turn off radiation since it is already accounted for in the forcing @@ -142,22 +142,22 @@ subroutine scm_main_sub() physics%Model%nsswr = -1 physics%Model%nslwr = -1 end if - + !override fhzero in physics namelist if n_itt_diag is set in the case namelist if (scm_state%n_itt_diag >= 1) then physics%Model%nszero = scm_state%n_itt_diag physics%Model%fhzero = scm_state%n_itt_diag*scm_state%dt/3600.0 end if - + !check for problematic diagnostic and radiation periods if (mod(physics%Model%nszero,scm_state%n_itt_out) /= 0) then write(*,*) "***ERROR***: The diagnostic output period must be a multiple of the output period." write(*,*) "From ", adjustl(trim(scm_state%physics_nml)), ", fhzero = ",physics%Model%fhzero write(*,*) "implying a diagnostic output period of ", physics%Model%nszero*scm_state%dt, "seconds." write(*,*) "The given output period in the case configuration namelist is ", scm_state%output_period,"seconds." - STOP + error stop end if - + if (mod(physics%Model%nsswr,scm_state%n_itt_out) /= 0) then write(*,*) "***WARNING***: The shortwave radiation calling period is different than the output period." write(*,*) "From ", adjustl(trim(scm_state%physics_nml)), ", fhswr = ",physics%Model%fhswr @@ -165,7 +165,7 @@ subroutine scm_main_sub() write(*,*) "This will cause the effective output period of variables that are only given values during shortwave calls to be ",& lcm(scm_state%n_itt_out,physics%Model%nsswr)*scm_state%dt," seconds." end if - + if (mod(physics%Model%nslwr,scm_state%n_itt_out) /= 0) then write(*,*) "***WARNING***: The longwave radiation calling period is different than the output period." write(*,*) "From ", adjustl(trim(scm_state%physics_nml)), ", fhlwr = ",physics%Model%fhlwr @@ -173,14 +173,14 @@ subroutine scm_main_sub() write(*,*) "This will cause the effective output period of variables that are only given values during longwave calls to be ",& lcm(scm_state%n_itt_out,physics%Model%nslwr)*scm_state%dt," seconds." end if - + cdata%blk_no = 1 cdata%thrd_no = 1 cdata%thrd_cnt = 1 call physics%associate(scm_state) call physics%set(scm_input_instance, scm_state) - + ! When asked to calculate 3-dim. tendencies, set Stateout variables to ! Statein variables here in order to capture the first call to dycore if (physics%Model%ldiag3d) then @@ -189,7 +189,7 @@ subroutine scm_main_sub() physics%Stateout%gt0 = physics%Statein%tgrs physics%Stateout%gq0 = physics%Statein%qgrs endif - + !initialize the column's physics write(0,'(a,i0,a)') "Calling ccpp_physics_init with suite '" // trim(trim(adjustl(scm_state%physics_suite_name))) // "'" @@ -197,9 +197,9 @@ subroutine scm_main_sub() write(0,'(a,i0,a,i0)') "Called ccpp_physics_init with suite '" // trim(trim(adjustl(scm_state%physics_suite_name))) // "', ierr=", ierr if (ierr/=0) then write(*,'(a,i0,a)') 'An error occurred in ccpp_physics_init: ' // trim(cdata%errmsg) // '. Exiting...' - stop 1 + error stop end if - + physics%Model%first_time_step = .true. call output_init(scm_state, physics) @@ -212,7 +212,7 @@ subroutine scm_main_sub() scm_state%dt_now = scm_state%dt scm_state%model_time = scm_state%dt_now end if - + call interpolate_forcing(scm_input_instance, scm_state, in_spinup) if (.not. scm_state%model_ics) call calc_pres_exner_geopotential(1, scm_state) @@ -259,19 +259,19 @@ subroutine scm_main_sub() physics%Diag%dtend(i,:,idtend) = physics%Diag%dtend(i,:,idtend) & + (physics%Statein%ugrs(i,:) - physics%Stateout%gu0(i,:)) endif - + idtend = physics%Model%dtidx(physics%Model%index_of_y_wind,physics%Model%index_of_process_non_physics) if(idtend>=1) then physics%Diag%dtend(i,:,idtend) = physics%Diag%dtend(i,:,idtend) & + (physics%Statein%vgrs(i,:) - physics%Stateout%gv0(i,:)) endif - + idtend = physics%Model%dtidx(physics%Model%index_of_temperature,physics%Model%index_of_process_non_physics) if(idtend>=1) then physics%Diag%dtend(i,:,idtend) = physics%Diag%dtend(i,:,idtend) & + (physics%Statein%tgrs(i,:) - physics%Stateout%gt0(i,:)) endif - + if (physics%Model%qdiag3d) then do itrac=1,physics%Model%ntrac idtend = physics%Model%dtidx(itrac+100,physics%Model%index_of_process_non_physics) @@ -283,13 +283,13 @@ subroutine scm_main_sub() endif endif end do - + call ccpp_physics_timestep_init(cdata, suite_name=trim(adjustl(scm_state%physics_suite_name)), ierr=ierr) if (ierr/=0) then write(*,'(a,i0,a)') 'An error occurred in ccpp_physics_timestep_init: ' // trim(cdata%errmsg) // '. Exiting...' - stop 1 + error stop end if - + !--- determine if radiation diagnostics buckets need to be cleared if (nint(physics%Model%fhzero*3600) >= nint(max(physics%Model%fhswr,physics%Model%fhlwr))) then if (mod(physics%Model%kdt,physics%Model%nszero) == 1 .or. physics%Model%nszero == 1) then @@ -301,22 +301,22 @@ subroutine scm_main_sub() call physics%Diag%rad_zero (physics%Model) endif endif - + !--- determine if physics diagnostics buckets need to be cleared if (mod(physics%Model%kdt,physics%Model%nszero) == 1 .or. physics%Model%nszero == 1) then call physics%Diag%phys_zero (physics%Model) endif - + call ccpp_physics_run(cdata, suite_name=trim(adjustl(scm_state%physics_suite_name)), ierr=ierr) if (ierr/=0) then write(*,'(a,i0,a)') 'An error occurred in ccpp_physics_run: ' // trim(cdata%errmsg) // '. Exiting...' - stop 1 + error stop end if - + call ccpp_physics_timestep_finalize(cdata, suite_name=trim(adjustl(scm_state%physics_suite_name)), ierr=ierr) if (ierr/=0) then write(*,'(a,i0,a)') 'An error occurred in ccpp_physics_timestep_finalize: ' // trim(cdata%errmsg) // '. Exiting...' - stop 1 + error stop end if !the filter routine (called after the following leapfrog time step) expects time level 2 in temp_tracer to be the updated, unfiltered state after the previous time step @@ -356,9 +356,9 @@ subroutine scm_main_sub() !prepare for time loop scm_state%n_timesteps = ceiling(scm_state%runtime/scm_state%dt) + scm_state%spinup_timesteps - + scm_state%dt_now = scm_state%dt - + !if (.not. in_spinup) then physics%Model%first_time_step = .false. !end if @@ -371,7 +371,7 @@ subroutine scm_main_sub() in_spinup = .false. scm_state%itt = i - scm_state%spinup_timesteps end if - + !> - Calculate the elapsed model time. if (.not. in_spinup) then scm_state%model_time = scm_state%itt*scm_state%dt @@ -392,7 +392,7 @@ subroutine scm_main_sub() end if call interpolate_forcing(scm_input_instance, scm_state, in_spinup) - + call calc_pres_exner_geopotential(1, scm_state) !zero out diagnostics output on EVERY time step - breaks diagnostics averaged over many timesteps @@ -410,7 +410,7 @@ subroutine scm_main_sub() scm_state%state_tracer(:,:,scm_state%cloud_water_index,1) = scm_state%state_tracer(:,:,scm_state%cloud_water_index,2) scm_state%state_tracer(:,:,scm_state%ozone_index,1) = scm_state%state_tracer(:,:,scm_state%ozone_index,2) end if - + write(*,*) "itt = ",scm_state%itt write(*,*) "model time (s) = ",scm_state%model_time if (scm_state%lsm_ics .or. scm_state%model_ics) then @@ -418,7 +418,7 @@ subroutine scm_main_sub() write(*,*) "sensible heat flux (W m-2): ",physics%Interstitial%dtsfc1(1) write(*,*) "latent heat flux (W m-2): ",physics%Interstitial%dqsfc1(1) end if - + if (.not. in_spinup) then call output_append(scm_state, physics) end if @@ -428,13 +428,13 @@ subroutine scm_main_sub() if (ierr/=0) then write(*,'(a,i0,a)') 'An error occurred in ccpp_physics_finalize: ' // trim(cdata%errmsg) // '. Exiting...' - stop 1 + error stop end if - + call MPI_FINALIZE(ierr) if (ierr/=0) then write(*,*) 'An error occurred in MPI_FINALIZE: ', ierr - stop 1 + error stop end if end subroutine scm_main_sub diff --git a/scm/src/scm_forcing.F90 b/scm/src/scm_forcing.F90 index ad7a47b4f..db239aa8b 100644 --- a/scm/src/scm_forcing.F90 +++ b/scm/src/scm_forcing.F90 @@ -39,10 +39,10 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) tot_advec_theta_bracket(2,scm_state%n_levels), tot_advec_thetal_bracket(2,scm_state%n_levels), tot_advec_qv_bracket(2,scm_state%n_levels), & tot_advec_u_bracket(2,scm_state%n_levels), tot_advec_v_bracket(2,scm_state%n_levels) !< forcing terms that "bracket" around the model time real(kind=dp) :: rho - + !> \section interpolate_forcing_alg Algorithm !! @{ - + if (in_spinup) then do i=1, scm_state%n_cols scm_state%pres_surf(i) = scm_input%input_pres_surf(1) @@ -56,15 +56,15 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) end do return end if - + !> - Check for the case where the elapsed model time extends beyond the supplied forcing. if(scm_state%model_time >= scm_input%input_time(scm_input%input_ntimes)) then !> - If so, hold the forcing terms constant at the last supplied values. The forcing still needs to be interpolated to the grid. write(*,*) "The model_time has exceeded the specifed period of forcing. Forcing will now be held constant at the last & specified values." - + if(scm_state%input_type == 0) then - + !> - For all forcing terms, call interpolate_to_grid_centers from \ref utils for each variable. This subroutine returns the last vertical index calculated in case forcing terms above the case input needs to be specified. do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres, & @@ -107,7 +107,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres, & scm_input%input_v_advec_qt(scm_input%input_ntimes,:), scm_state%pres_l(i,:), scm_state%n_levels, & v_advec_qt_bracket(1,:), top_index, 3) - + !> - If the input forcing file does not reach to the model domain top, fill in values above the input forcing file domain with those from the top level. if (top_index < scm_state%n_levels .AND. top_index.GT.0) then w_ls_bracket(1,top_index+1:scm_state%n_levels) = 0.0!w_ls_bracket(1,top_index) @@ -125,7 +125,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) v_advec_thil_bracket(1,top_index+1:scm_state%n_levels) = v_advec_thil_bracket(1,top_index) v_advec_qt_bracket(1,top_index+1:scm_state%n_levels) = v_advec_qt_bracket(1,top_index) end if - + !> - For this case, no time interpolation is necessary; just set the forcing terms to the vertically-interpolated values. scm_state%w_ls(i,:) = w_ls_bracket(1,:) scm_state%omega(i,:) = omega_bracket(1,:) @@ -141,7 +141,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) scm_state%h_advec_qt(i,:) = h_advec_qt_bracket(1,:) scm_state%v_advec_thil(i,:) = v_advec_thil_bracket(1,:) scm_state%v_advec_qt(i,:) = v_advec_qt_bracket(1,:) - + !> - Set the surface parameters to the last available data. scm_state%pres_surf(i) = scm_input%input_pres_surf(scm_input%input_ntimes) scm_state%T_surf(i) = scm_input%input_T_surf(scm_input%input_ntimes) @@ -172,7 +172,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) scm_state%w_ls(i,:) = w_ls_bracket(1,:) end do end if - + if (scm_state%force_geo) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(scm_input%input_ntimes,:), & @@ -187,7 +187,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) scm_state%v_g(i,:) = v_g_bracket(1,:) end do end if - + if (scm_state%force_adv_T == 1) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(scm_input%input_ntimes,:), & @@ -219,7 +219,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) scm_state%tot_advec_thetal(i,:) = tot_advec_thetal_bracket(1,:) end do end if - + if (scm_state%force_adv_qv) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(scm_input%input_ntimes,:), & @@ -231,7 +231,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) scm_state%tot_advec_qv(i,:) = tot_advec_qv_bracket(1,:) end do end if - + if (scm_state%force_adv_u) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(scm_input%input_ntimes,:), scm_input%input_tot_advec_u(scm_input%input_ntimes,:), & @@ -242,7 +242,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) scm_state%tot_advec_u(i,:) = tot_advec_u_bracket(1,:) end do end if - + if (scm_state%force_adv_v) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(scm_input%input_ntimes,:), scm_input%input_tot_advec_v(scm_input%input_ntimes,:), & @@ -253,7 +253,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) scm_state%tot_advec_v(i,:) = tot_advec_v_bracket(1,:) end do end if - + if (scm_state%force_nudging_t == 1) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(scm_input%input_ntimes,:), scm_input%input_T_nudge(scm_input%input_ntimes,:), & @@ -264,7 +264,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) scm_state%T_nudge(i,:) = T_nudge_bracket(1,:) call find_vertical_index_pressure(scm_input%input_pres_forcing(scm_input%input_ntimes,scm_input%input_k_T_nudge(scm_input%input_ntimes)), scm_state%pres_l(i,:), scm_state%force_nudging_T_k(i)) end do - + else if (scm_state%force_nudging_T == 2 .or. scm_state%force_nudging_T == 3) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(scm_input%input_ntimes,:), scm_input%input_thil_nudge(scm_input%input_ntimes,:), & @@ -276,7 +276,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) call find_vertical_index_pressure(scm_input%input_pres_forcing(scm_input%input_ntimes,scm_input%input_k_thil_nudge(scm_input%input_ntimes)), scm_state%pres_l(i,:), scm_state%force_nudging_T_k(i)) end do end if - + if (scm_state%force_nudging_qv) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(scm_input%input_ntimes,:), scm_input%input_qt_nudge(scm_input%input_ntimes,:), & @@ -288,7 +288,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) call find_vertical_index_pressure(scm_input%input_pres_forcing(scm_input%input_ntimes,scm_input%input_k_qt_nudge(scm_input%input_ntimes)), scm_state%pres_l(i,:), scm_state%force_nudging_qv_k(i)) end do end if - + if (scm_state%force_nudging_u) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(scm_input%input_ntimes,:), scm_input%input_u_nudge(scm_input%input_ntimes,:), & @@ -300,7 +300,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) call find_vertical_index_pressure(scm_input%input_pres_forcing(scm_input%input_ntimes,scm_input%input_k_u_nudge(scm_input%input_ntimes)), scm_state%pres_l(i,:), scm_state%force_nudging_u_k(i)) end do end if - + if (scm_state%force_nudging_v) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(scm_input%input_ntimes,:), scm_input%input_v_nudge(scm_input%input_ntimes,:), & @@ -312,7 +312,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) call find_vertical_index_pressure(scm_input%input_pres_forcing(scm_input%input_ntimes,scm_input%input_k_v_nudge(scm_input%input_ntimes)), scm_state%pres_l(i,:), scm_state%force_nudging_v_k(i)) end do end if - + if (scm_state%force_rad_T == 1 .or. scm_state%force_rad_T == 2 .or. scm_state%force_rad_T == 3) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(scm_input%input_ntimes,:), scm_input%input_dT_dt_rad(scm_input%input_ntimes,:), & @@ -323,14 +323,14 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) scm_state%dT_dt_rad(i,:) = dT_dt_rad_bracket(1,:) end do end if - + if (scm_state%surface_thermo_control == 0 .or. scm_state%surface_thermo_control == 1 .or. scm_state%surface_thermo_control == 2) then !skin temperature is needed if surface fluxes are specified (for calculating bulk Richardson number in the specified surface flux scheme) and for simple ocean scheme do i=1, scm_state%n_cols scm_state%T_surf(i) = scm_input%input_T_surf(scm_input%input_ntimes) end do end if - + if (scm_state%surface_thermo_control == 0) then do i=1, scm_state%n_cols scm_state%sh_flux(i) = scm_input%input_sh_flux_sfc_kin(scm_input%input_ntimes) @@ -344,11 +344,11 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) scm_state%lh_flux(i) = (1.0/(con_hvap*rho))*scm_input%input_lh_flux_sfc(scm_input%input_ntimes) end do end if - + do i=1, scm_state%n_cols scm_state%pres_surf(i) = scm_input%input_pres_surf(scm_input%input_ntimes) end do - + end if else !> - When the model elapsed time is within the time-frame specified by the input forcing, the forcing must be interpolated in time and space. @@ -362,9 +362,9 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) exit end if end do - + if(scm_state%input_type == 0) then - + do i=1, scm_state%n_cols !> - For all forcing terms, call interpolate_to_grid_centers from \ref utils for each variable for each time level that "bracket" around !> the current model time. This subroutine returns the last vertical index calculated in case forcing terms above the case input needs @@ -499,7 +499,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) scm_state%omega(i,:) = (1.0 - lifrac)*omega_bracket(1,:) + lifrac*omega_bracket(2,:) end do end if - + if (scm_state%force_w) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(low_t_index,:), scm_input%input_w_ls(low_t_index,:), & @@ -513,7 +513,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) scm_state%w_ls(i,:) = (1.0 - lifrac)*w_ls_bracket(1,:) + lifrac*w_ls_bracket(2,:) end do end if - + if (scm_state%force_geo) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(low_t_index,:), scm_input%input_u_g(low_t_index,:), & @@ -534,7 +534,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) scm_state%v_g(i,:) = (1.0 - lifrac)*v_g_bracket(1,:) + lifrac*v_g_bracket(2,:) end do end if - + if (scm_state%force_adv_T == 1) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(low_t_index,:), scm_input%input_tot_advec_t(low_t_index,:), & @@ -572,7 +572,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) scm_state%tot_advec_thetal(i,:) = (1.0 - lifrac)*tot_advec_thetal_bracket(1,:) + lifrac*tot_advec_thetal_bracket(2,:) end do end if - + if (scm_state%force_adv_qv) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(low_t_index,:), scm_input%input_tot_advec_qv(low_t_index,:), & @@ -586,7 +586,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) scm_state%tot_advec_qv(i,:) = (1.0 - lifrac)*tot_advec_qv_bracket(1,:) + lifrac*tot_advec_qv_bracket(2,:) end do end if - + if (scm_state%force_adv_u) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(low_t_index,:), scm_input%input_tot_advec_u(low_t_index,:), & @@ -600,7 +600,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) scm_state%tot_advec_u(i,:) = (1.0 - lifrac)*tot_advec_u_bracket(1,:) + lifrac*tot_advec_u_bracket(2,:) end do end if - + if (scm_state%force_adv_v) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(low_t_index,:), scm_input%input_tot_advec_v(low_t_index,:), & @@ -614,7 +614,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) scm_state%tot_advec_v(i,:) = (1.0 - lifrac)*tot_advec_v_bracket(1,:) + lifrac*tot_advec_v_bracket(2,:) end do end if - + if (scm_state%force_nudging_t == 1) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(low_t_index,:), scm_input%input_T_nudge(low_t_index,:), & @@ -642,7 +642,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) call find_vertical_index_pressure(scm_input%input_pres_forcing(low_t_index,scm_input%input_k_thil_nudge(low_t_index)), scm_state%pres_l(i,:), scm_state%force_nudging_T_k(i)) end do end if - + if (scm_state%force_nudging_qv) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(low_t_index,:), scm_input%input_qt_nudge(low_t_index,:), & @@ -657,7 +657,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) call find_vertical_index_pressure(scm_input%input_pres_forcing(low_t_index,scm_input%input_k_qt_nudge(low_t_index)), scm_state%pres_l(i,:), scm_state%force_nudging_qv_k(i)) end do end if - + if (scm_state%force_nudging_u) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(low_t_index,:), scm_input%input_u_nudge(low_t_index,:), & @@ -672,7 +672,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) call find_vertical_index_pressure(scm_input%input_pres_forcing(low_t_index,scm_input%input_k_u_nudge(low_t_index)), scm_state%pres_l(i,:), scm_state%force_nudging_u_k(i)) end do end if - + if (scm_state%force_nudging_v) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(low_t_index,:), scm_input%input_v_nudge(low_t_index,:), & @@ -687,7 +687,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) call find_vertical_index_pressure(scm_input%input_pres_forcing(low_t_index,scm_input%input_k_v_nudge(low_t_index)), scm_state%pres_l(i,:), scm_state%force_nudging_v_k(i)) end do end if - + if (scm_state%force_rad_T == 1 .or. scm_state%force_rad_T == 2 .or. scm_state%force_rad_T == 3) then do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres_forcing(low_t_index,:), scm_input%input_dT_dt_rad(low_t_index,:), & @@ -701,14 +701,14 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) scm_state%dT_dt_rad(i,:) = (1.0 - lifrac)*dT_dt_rad_bracket(1,:) + lifrac*dT_dt_rad_bracket(2,:) end do end if - + if (scm_state%surface_thermo_control == 0 .or. scm_state%surface_thermo_control == 1 .or. scm_state%surface_thermo_control == 2) then !skin temperature is needed if surface fluxes are specified (for calculating bulk Richardson number in the specified surface flux scheme) and for simple ocean scheme do i=1, scm_state%n_cols scm_state%T_surf(i) = (1.0 - lifrac)*scm_input%input_T_surf(low_t_index) + lifrac*scm_input%input_T_surf(low_t_index+1) end do end if - + if (scm_state%surface_thermo_control == 0) then do i=1, scm_state%n_cols scm_state%sh_flux(i) = (1.0 - lifrac)*scm_input%input_sh_flux_sfc_kin(low_t_index) + & @@ -726,7 +726,7 @@ subroutine interpolate_forcing(scm_input, scm_state, in_spinup) lifrac*scm_input%input_lh_flux_sfc(low_t_index+1)) end do end if - + do i=1, scm_state%n_cols !Interpolate the surface parameters in time. scm_state%pres_surf(i) = (1.0 - lifrac)*scm_input%input_pres_surf(low_t_index) + & @@ -792,7 +792,7 @@ subroutine apply_forcing_leapfrog(scm_state) select case(scm_state%mom_forcing_type) case (1) write(*,*) 'momentum forcing type = 1 is not implemented. Pick 2 or 3. Stopping...' - stop + error stop case (2) !> - Calculate change in state momentum variables due to vertical advection (subsidence). @@ -935,12 +935,12 @@ subroutine apply_forcing_forward_Euler(scm_state, in_spinup) integer :: i,k real(kind=dp) :: f_coriolis, grav_inv, g_over_cp, omega_plus, omega_minus, dth_dp_plus, dth_dp_minus, & dqv_dp_plus, dqv_dp_minus, spinup_relax_time - + !> \section apply_leapfrog_forcing_alg Algorithm !! @{ - + spinup_relax_time = scm_state%dt - + grav_inv = 1.0/con_g g_over_cp = con_g/con_cp @@ -971,7 +971,7 @@ subroutine apply_forcing_forward_Euler(scm_state, in_spinup) zi(i,scm_state%n_levels+1) = scm_state%geopotential_i(i,scm_state%n_levels+1)*grav_inv end do !end if - + if (in_spinup) then do i=1, scm_state%n_cols do k=1, scm_state%n_levels @@ -991,10 +991,10 @@ subroutine apply_forcing_forward_Euler(scm_state, in_spinup) select case(scm_state%mom_forcing_type) case (1) write(*,*) 'momentum forcing type = 1 is not implemented. Pick 2 or 3. Stopping...' - stop + error stop case (2) !> - Calculate change in state momentum variables due to vertical advection (subsidence). - + !> - Calculate tendencies due to vertical advection using same discretization as in previous GFS SCM implmentation (staggered central difference) !! \f[ !! \frac{\partial x}{\partial t}|_{vert. advection} = \frac{w_{k+1}\left(x_{k+1} - x_{k}\right) + w_k\left(x_k - x_{k-1}\right)}{-2\left(z_{k+1}-z_{k}\right)} @@ -1014,7 +1014,7 @@ subroutine apply_forcing_forward_Euler(scm_state, in_spinup) scm_state%v_force_tend(i,1) = -w_ls_i(i,2)*(old_v(i,2) - old_v(i,1))/(zi(i,2)-zi(i,1)) scm_state%v_force_tend(i,scm_state%n_levels) = -w_ls_i(i,scm_state%n_levels)*& (old_v(i,scm_state%n_levels) - old_v(i,scm_state%n_levels-1))/(zi(i,scm_state%n_levels+1)-zi(i,scm_state%n_levels)) - + !> - Add forcing due to geostrophic wind !> - Calculate Coriolis parameter. f_coriolis = 2.0*con_omega*sin(scm_state%lat(i)) @@ -1037,7 +1037,7 @@ subroutine apply_forcing_forward_Euler(scm_state, in_spinup) scm_state%u_force_tend = 0.0 scm_state%v_force_tend = 0.0 end select - + select case (scm_state%thermo_forcing_type) case (1) do i=1, scm_state%n_cols @@ -1061,7 +1061,7 @@ subroutine apply_forcing_forward_Euler(scm_state, in_spinup) scm_state%qv_force_tend(i,k) = -omega_plus*dqv_dp_minus - omega_minus*dqv_dp_plus scm_state%T_force_tend(i,k) = scm_state%exner_l(i,k)*(-omega_plus*dth_dp_minus - omega_minus*dth_dp_plus) end do - + !> - Add forcing due to prescribed radiation and horizontal advection do k=1, scm_state%n_levels scm_state%T_force_tend(i,k) = scm_state%T_force_tend(i,k) + scm_state%dT_dt_rad(i,k) + & @@ -1077,7 +1077,7 @@ subroutine apply_forcing_forward_Euler(scm_state, in_spinup) scm_state%T_force_tend(i,k) = (scm_state%T_nudge(i,k) - old_T(i,k))/scm_state%relax_time scm_state%qv_force_tend(i,k) = (scm_state%qt_nudge(i,k) - old_qv(i,k))/scm_state%relax_time end do - + do k=2, scm_state%n_levels-1 !upstream scheme (for boundaries, assume vertical derivatives are 0 => no vertical advection) omega_plus = MAX(scm_state%omega(i,k), 0.0) @@ -1120,28 +1120,28 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) type(scm_state_type), intent(inout) :: scm_state logical, intent(in) :: in_spinup - + real(kind=dp) :: old_u(scm_state%n_cols, scm_state%n_levels), old_v(scm_state%n_cols, scm_state%n_levels), & old_T(scm_state%n_cols, scm_state%n_levels), old_qv(scm_state%n_cols, scm_state%n_levels) real(kind=dp) :: theta(scm_state%n_cols, scm_state%n_levels) - + real(kind=dp) :: spinup_relax_time, omega_asc, omega_des, w_asc, w_des, gradient_asc, gradient_des, rho, adiabatic_exp_comp_term, & f_coriolis - + integer :: i,k - - logical :: use_theta !formulations using potential temperature don't need adiabatic expansion/compression term (simpler), + + logical :: use_theta !formulations using potential temperature don't need adiabatic expansion/compression term (simpler), !but reqires conversion to/from since absolute temperature is state variable - + use_theta = .false. spinup_relax_time = scm_state%dt - + !Save old state variables old_u = scm_state%state_u(:,:,1) old_v = scm_state%state_v(:,:,1) old_T = scm_state%state_T(:,:,1) old_qv = scm_state%state_tracer(:,:,scm_state%water_vapor_index,1) - + if (use_theta) then theta = old_T/scm_state%exner_l(:,:) end if @@ -1151,7 +1151,7 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) scm_state%v_force_tend = 0.0 scm_state%T_force_tend = 0.0 scm_state%qv_force_tend = 0.0 - + if (in_spinup) then do i=1, scm_state%n_cols do k=1, scm_state%n_levels @@ -1194,7 +1194,7 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) end do end if !use_theta end if !force_sub_for_T - + if (scm_state%force_sub_for_qv) then do i=1, scm_state%n_cols do k=2, scm_state%n_levels-1 @@ -1206,7 +1206,7 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) end do end do end if - + if (scm_state%force_sub_for_u) then do i=1, scm_state%n_cols do k=2, scm_state%n_levels-1 @@ -1218,7 +1218,7 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) end do end do end if - + if (scm_state%force_sub_for_v) then do i=1, scm_state%n_cols do k=2, scm_state%n_levels-1 @@ -1230,7 +1230,7 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) end do end do end if - + else if (scm_state%force_w) then if (scm_state%force_sub_for_T) then if (use_theta) then @@ -1258,7 +1258,7 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) end do end if !use_theta end if !force_sub_for_T - + if (scm_state%force_sub_for_qv) then do i=1, scm_state%n_cols do k=2, scm_state%n_levels-1 @@ -1271,7 +1271,7 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) end do end do end if - + if (scm_state%force_sub_for_u) then do i=1, scm_state%n_cols do k=2, scm_state%n_levels-1 @@ -1284,7 +1284,7 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) end do end do end if - + if (scm_state%force_sub_for_v) then do i=1, scm_state%n_cols do k=2, scm_state%n_levels-1 @@ -1297,12 +1297,12 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) end do end do end if - + end if !force_omega or force_w - + if (scm_state%force_geo) then !Add forcing due to geostrophic wind - + !Calculate Coriolis parameter. do i=1, scm_state%n_cols f_coriolis = 2.0*con_omega*sin(scm_state%lat(i)) @@ -1313,7 +1313,7 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) end do end do end if !force_geo - + if (scm_state%force_adv_T == 1) then !advection term is in terms of absolute temperature do i=1, scm_state%n_cols @@ -1336,7 +1336,7 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) end do end do end if - + if (scm_state%force_adv_qv) then do i=1, scm_state%n_cols do k=1, scm_state%n_levels @@ -1344,7 +1344,7 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) end do end do end if - + if (scm_state%force_adv_u) then do i=1, scm_state%n_cols do k=1, scm_state%n_levels @@ -1352,7 +1352,7 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) end do end do end if - + if (scm_state%force_adv_v) then do i=1, scm_state%n_cols do k=1, scm_state%n_levels @@ -1360,7 +1360,7 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) end do end do end if - + if (scm_state%force_nudging_t == 1) then do i=1, scm_state%n_cols do k=scm_state%force_nudging_T_k(i), scm_state%n_levels @@ -1375,7 +1375,7 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) end do end do end if - + if (scm_state%force_nudging_qv) then do i=1, scm_state%n_cols do k=scm_state%force_nudging_qv_k(i), scm_state%n_levels @@ -1383,7 +1383,7 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) end do end do end if - + if (scm_state%force_nudging_u) then do i=1, scm_state%n_cols do k=scm_state%force_nudging_u_k(i), scm_state%n_levels @@ -1391,7 +1391,7 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) end do end do end if - + if (scm_state%force_nudging_v) then do i=1, scm_state%n_cols do k=scm_state%force_nudging_v_k(i), scm_state%n_levels @@ -1399,7 +1399,7 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) end do end do end if - + if (scm_state%force_rad_T == 1 .or. scm_state%force_rad_T == 2 .or. scm_state%force_rad_T == 3) then do i=1, scm_state%n_cols do k=1, scm_state%n_levels @@ -1407,9 +1407,9 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) end do end do end if - + end if !in_spinup - + do i=1, scm_state%n_cols do k=1, scm_state%n_levels !> - Update the state variables using the forward Euler scheme: @@ -1424,23 +1424,23 @@ subroutine apply_forcing_DEPHY(scm_state, in_spinup) scm_state%dt*(scm_state%qv_force_tend(i,k)) end do end do - + end subroutine apply_forcing_DEPHY subroutine set_spinup_nudging(scm_state) use scm_type_defs, only: scm_state_type type(scm_state_type), intent(inout) :: scm_state - + integer :: i - + do i=1, scm_state%n_cols scm_state%u_nudge(i,:) = scm_state%state_u(i,:,1) scm_state%v_nudge(i,:) = scm_state%state_v(i,:,1) scm_state%T_nudge(i,:) = scm_state%state_T(i,:,1) scm_state%qt_nudge(i,:) = scm_state%state_tracer(i,:,scm_state%water_vapor_index,1) end do - + end subroutine set_spinup_nudging !> @} diff --git a/scm/src/scm_input.F90 b/scm/src/scm_input.F90 index cc3effcd5..abaeedefc 100644 --- a/scm/src/scm_input.F90 +++ b/scm/src/scm_input.F90 @@ -67,7 +67,7 @@ subroutine get_config_nml(scm_state) character(len=character_length) :: physics_suite !< name of the physics suite name (currently only GFS_operational supported) character(len=character_length) :: physics_nml - + character(len=character_length), allocatable, dimension(:) :: tracer_names integer, allocatable, dimension(:) :: tracer_types @@ -80,7 +80,7 @@ subroutine get_config_nml(scm_state) lsm_ics, do_spinup, C_RES, spinup_timesteps, mom_forcing_type, relax_time, sfc_type, sfc_flux_spec, & sfc_roughness_length_cm, reference_profile_choice, year, month, day, hour, min, & column_area, input_type - + NAMELIST /physics_config/ physics_suite, physics_nml !> \section get_config_alg Algorithm @@ -119,12 +119,12 @@ subroutine get_config_nml(scm_state) hour = 3 min = 0 input_type = 0 - + open(unit=10, file=experiment_namelist, status='old', action='read', iostat=ioerror) if(ioerror /= 0) then write(*,'(a,i0)') 'There was an error opening the file ' // experiment_namelist // & '; error code = ', ioerror - STOP + error stop "error opening namelist" else read(10, NML=case_config, iostat=ioerror) end if @@ -132,7 +132,7 @@ subroutine get_config_nml(scm_state) if(ioerror /= 0) then write(*,'(a,i0)') 'There was an error reading the namelist case_config in the file '& // experiment_namelist // '; error code = ',ioerror - STOP + error stop "error opening namelist" end if !The current implementation of GFS physics does not support more than one column, since radiation sub schemes use @@ -140,10 +140,9 @@ subroutine get_config_nml(scm_state) !the code crashes in GFS_initialize_scm_run and later in radiation_gases.f, because it tries to allocate module !variables that are already allocated. For now, throw an error and abort. if (n_columns>1) then - write(*,'(a)') 'The current implementation does not allow to run more than one column at a time.' - STOP + error stop "The current implementation does not allow to run more than one column at a time." end if - + !read in the physics suite and namelist read(10, NML=physics_config, iostat=ioerror) close(10) @@ -156,9 +155,9 @@ subroutine get_config_nml(scm_state) case default n_time_levels = 2 end select - + call get_tracers(tracer_names, tracer_types) - + call scm_state%create(n_columns, n_levels, n_soil, n_snow, n_time_levels, tracer_names, tracer_types) scm_state%experiment_name = experiment_name @@ -199,7 +198,7 @@ subroutine get_config_nml(scm_state) scm_state%reference_profile_choice = reference_profile_choice scm_state%relax_time = relax_time scm_state%input_type = input_type - + deallocate(tracer_names) !> @} end subroutine get_config_nml @@ -212,7 +211,7 @@ subroutine get_case_init(scm_state, scm_input) use NetCDF_read, only: NetCDF_read_var, check, missing_value type(scm_state_type), intent(in) :: scm_state type(scm_input_type), target, intent(inout) :: scm_input - + integer :: input_nlev !< number of levels in the input file integer :: input_nsoil !< number of soil levels in the input file integer :: input_nsnow !< number of snow levels in the input file @@ -222,9 +221,9 @@ subroutine get_case_init(scm_state, scm_input) real(kind=dp) :: input_lat !< column latitude (deg) real(kind=dp) :: input_lon !< column longitude (deg) real(kind=dp) :: input_area !< surface area [m^2] - + integer :: input_vegsrc !< vegetation source - + real(kind=dp) :: input_slmsk !< sea land ice mask [0,1,2] real(kind=dp) :: input_tsfco !< input sea surface temperature OR surface skin temperature over land OR surface skin temperature over ice (depending on slmsk) (K) real(kind=dp) :: input_weasd !< water equivalent accumulated snow depth (mm) @@ -275,14 +274,14 @@ subroutine get_case_init(scm_state, scm_input) real(kind=dp) :: input_albdifvis_ice !< real(kind=dp) :: input_albdifnir_ice !< real(kind=dp) :: input_zorlwav !< surface roughness length from wave model (cm) - + real(kind=dp), allocatable :: input_stc(:) !< soil temperature (K) - real(kind=dp), allocatable :: input_smc(:) !< total soil moisture content (fraction) + real(kind=dp), allocatable :: input_smc(:) !< total soil moisture content (fraction) real(kind=dp), allocatable :: input_slc(:) !< liquid soil moisture content (fraction) real(kind=dp), allocatable :: input_tiice(:) !< sea ice internal temperature (K) real(kind=dp) :: input_stddev !< standard deviation of subgrid orography (m) - real(kind=dp) :: input_convexity !< convexity of subgrid orography + real(kind=dp) :: input_convexity !< convexity of subgrid orography real(kind=dp) :: input_ol1 !< fraction of grid box with subgrid orography higher than critical height 1 real(kind=dp) :: input_ol2 !< fraction of grid box with subgrid orography higher than critical height 2 real(kind=dp) :: input_ol3 !< fraction of grid box with subgrid orography higher than critical height 3 @@ -330,13 +329,13 @@ subroutine get_case_init(scm_state, scm_input) real(kind=dp) :: input_deeprechxy !< recharge to or from the water table when deep (m) real(kind=dp) :: input_rechxy !< recharge to or from the water table when shallow (m) real(kind=dp) :: input_snowxy !< number of snow layers - + real(kind=dp), allocatable :: input_snicexy(:) !< snow layer ice (mm) real(kind=dp), allocatable :: input_snliqxy(:) !< snow layer liquid (mm) real(kind=dp), allocatable :: input_tsnoxy(:) !< snow temperature (K) real(kind=dp), allocatable :: input_smoiseq(:) !< equilibrium soil water content (m3 m-3) real(kind=dp), allocatable :: input_zsnsoxy(:) !< layer bottom depth from snow surface (m) - + real(kind=dp) :: input_tref !< sea surface reference temperature for NSST (K) real(kind=dp) :: input_z_c !< sub-layer cooling thickness for NSST (m) real(kind=dp) :: input_c_0 !< coefficient 1 to calculate d(Tz)/d(Ts) for NSST @@ -355,7 +354,7 @@ subroutine get_case_init(scm_state, scm_input) real(kind=dp) :: input_ifd !< index to start DTM run for NSST real(kind=dp) :: input_dt_cool !< sub-layer cooling amount for NSST (K) real(kind=dp) :: input_qrain !< sensible heat due to rainfall for NSST (W) - + real(kind=dp) :: input_wetness !< normalized soil wetness for RUC LSM real(kind=dp) :: input_clw_surf_land !< cloud condensed water mixing ratio at surface over land for RUC LSM (kg kg-1) real(kind=dp) :: input_clw_surf_ice !< cloud condensed water mixing ratio at surface over ice for RUC LSM (kg kg-1) @@ -371,13 +370,13 @@ subroutine get_case_init(scm_state, scm_input) real(kind=dp) :: input_sfalb_ice real(kind=dp) :: input_emis_ice real(kind=dp) :: input_lai !< leaf area index for RUC LSM - + real(kind=dp), allocatable :: input_tslb(:) !< soil temperature for RUC LSM (K) real(kind=dp), allocatable :: input_smois(:) !< volume fraction of soil moisture for RUC LSM (frac) real(kind=dp), allocatable :: input_sh2o(:) !< volume fraction of unfrozen soil moisture for RUC LSM (frac) real(kind=dp), allocatable :: input_smfr(:) !< volume fraction of frozen soil moisture for RUC LSM (frac) real(kind=dp), allocatable :: input_flfr(:) !< flag for frozen soil physics - + ! dimension variables !real(kind=dp), allocatable :: input_pres_i(:) !< interface pressures !real(kind=dp), allocatable :: input_pres_l(:) !< layer pressures @@ -394,7 +393,7 @@ subroutine get_case_init(scm_state, scm_input) real(kind=dp), allocatable :: input_v(:) !< north-south horizontal wind profile (m s^-1) real(kind=dp), allocatable :: input_tke(:) !< TKE profile (m^2 s^-2) real(kind=dp), allocatable :: input_ozone(:) !< ozone profile (kg kg^-1) - + real(kind=dp), allocatable :: input_pres_surf(:) !< time-series of surface pressure (Pa) real(kind=dp), allocatable :: input_T_surf(:) !< time-series of surface temperature (K) @@ -415,7 +414,7 @@ subroutine get_case_init(scm_state, scm_input) real(kind=dp), allocatable :: input_sh_flux_sfc(:) !< time-series of surface sensible heat flux (K m s^-1) real(kind=dp), allocatable :: input_lh_flux_sfc(:) !< time-series of surface latent heat flux (kg kg^-1 m s^-1) - + CHARACTER(LEN=nf90_max_name) :: tmpName integer :: ncid, varID, grp_ncid, allocate_status,ierr real(kind=dp) :: nc_missing_value @@ -425,21 +424,21 @@ subroutine get_case_init(scm_state, scm_input) !> - Open the case input file found in the processed_case_input dir corresponding to the experiment name. call check(NF90_OPEN(trim(adjustl(scm_state%case_name))//'.nc',nf90_nowrite,ncid),"nf90_open()") - + !> - Read in missing value from file (replace module variable if present) ierr = NF90_GET_ATT(ncid, NF90_GLOBAL, 'missing_value', nc_missing_value) if(ierr == NF90_NOERR) then missing_value = nc_missing_value end if - + !> - Get the dimensions (global group). - + !required dimensions call check(NF90_INQ_DIMID(ncid,"levels",varID),"nf90_inq_dimid(levels)") call check(NF90_INQUIRE_DIMENSION(ncid, varID, tmpName, input_nlev),"nf90_inq_dim(levels)") call check(NF90_INQ_DIMID(ncid,"time",varID),"inq_dimid(time)") call check(NF90_INQUIRE_DIMENSION(ncid, varID, tmpName, input_ntimes),"nf90_inq_dim(time)") - + !possible dimensions (if using model ICs) ierr = NF90_INQ_DIMID(ncid,"nsoil",varID) if(ierr /= NF90_NOERR) then @@ -458,8 +457,8 @@ subroutine get_case_init(scm_state, scm_input) input_nice = missing_ice_layers else call check(NF90_INQUIRE_DIMENSION(ncid, varID, tmpName, input_nice),"nf90_inq_dim(nice)") - end if - + end if + !> - Allocate the dimension variables. allocate(input_pres(input_nlev),input_time(input_ntimes), stat=allocate_status) @@ -475,15 +474,15 @@ subroutine get_case_init(scm_state, scm_input) !> - Allocate the initial profiles (required). One of thetail or temp is required. allocate(input_thetail(input_nlev), input_temp(input_nlev), input_qt(input_nlev), input_ql(input_nlev), input_qi(input_nlev), & input_u(input_nlev), input_v(input_nlev), input_tke(input_nlev), input_ozone(input_nlev), stat=allocate_status) - + !> - Read in the initial profiles. The variable names in all input files are expected to be identical. - + !Either thetail or T must be present call NetCDF_read_var(grp_ncid, "thetail", .False., input_thetail) call NetCDF_read_var(grp_ncid, "temp", .False., input_temp) if (maxval(input_thetail) < 0 .and. maxval(input_temp) < 0) then write(*,*) "One of thetail or temp variables must be present in ",trim(adjustl(scm_state%case_name))//'.nc',". Stopping..." - STOP + error stop "One of thetail or temp variables" end if call NetCDF_read_var(grp_ncid, "qt", .True., input_qt ) call NetCDF_read_var(grp_ncid, "ql", .True., input_ql ) @@ -492,7 +491,7 @@ subroutine get_case_init(scm_state, scm_input) call NetCDF_read_var(grp_ncid, "v", .True., input_v ) call NetCDF_read_var(grp_ncid, "tke", .True., input_tke ) call NetCDF_read_var(grp_ncid, "ozone", .True., input_ozone) - + !possible initial profiles !needed for Noah LSM and others (when running with model ICs) allocate(input_stc(input_nsoil), input_smc(input_nsoil), input_slc(input_nsoil), & @@ -500,7 +499,7 @@ subroutine get_case_init(scm_state, scm_input) call NetCDF_read_var(grp_ncid, "stc", .False., input_stc) call NetCDF_read_var(grp_ncid, "smc", .False., input_smc) call NetCDF_read_var(grp_ncid, "slc", .False., input_slc) - + !needed for NoahMP LSM (when running with model ICs) allocate(input_snicexy(input_nsnow), input_snliqxy(input_nsnow), input_tsnoxy(input_nsnow), & input_smoiseq(input_nsoil), input_zsnsoxy(input_nsnow + input_nsoil)) @@ -509,11 +508,11 @@ subroutine get_case_init(scm_state, scm_input) call NetCDF_read_var(grp_ncid, "tsnoxy", .False., input_tsnoxy ) call NetCDF_read_var(grp_ncid, "smoiseq", .False., input_smoiseq) call NetCDF_read_var(grp_ncid, "zsnsoxy", .False., input_zsnsoxy) - + !needed for fractional grid (when running with model ICs) allocate(input_tiice(input_nice)) call NetCDF_read_var(grp_ncid, "tiice", .False., input_tiice) - + !needed for RUC LSM (when running with model ICs) allocate(input_tslb(input_nsoil), input_smois(input_nsoil), input_sh2o(input_nsoil), & input_smfr(input_nsoil), input_flfr(input_nsoil)) @@ -522,16 +521,16 @@ subroutine get_case_init(scm_state, scm_input) call NetCDF_read_var(grp_ncid, "sh2o", .False., input_sh2o ) call NetCDF_read_var(grp_ncid, "smfr", .False., input_smfr ) call NetCDF_read_var(grp_ncid, "flfr", .False., input_flfr ) - + !> - Find group ncid for scalar group. call check(NF90_INQ_GRP_NCID(ncid,"scalars",grp_ncid),"nf90_inq_grp_ncid(scalars)") - + !required call NetCDF_read_var(grp_ncid, "lat", .True., input_lat) call NetCDF_read_var(grp_ncid, "lon", .True., input_lon) !time data in file ignored? call NetCDF_read_var(grp_ncid, "area", .False., input_area) - + !possible scalars !Noah LSM parameters (when running with model ICs) call NetCDF_read_var(grp_ncid, "vegsrc", .False., input_vegsrc ) @@ -585,7 +584,7 @@ subroutine get_case_init(scm_state, scm_input) call NetCDF_read_var(grp_ncid, "albdifvis_ice", .False., input_albdifvis_ice) call NetCDF_read_var(grp_ncid, "albdifnir_ice", .False., input_albdifnir_ice) call NetCDF_read_var(grp_ncid, "zorlwav", .False., input_zorlwav) - + !orographic parameters call NetCDF_read_var(grp_ncid, "stddev", .False., input_stddev) call NetCDF_read_var(grp_ncid, "convexity", .False., input_convexity) @@ -606,7 +605,7 @@ subroutine get_case_init(scm_state, scm_input) call NetCDF_read_var(grp_ncid, "landfrac", .False., input_landfrac) call NetCDF_read_var(grp_ncid, "lakefrac", .False., input_lakefrac) call NetCDF_read_var(grp_ncid, "lakedepth", .False., input_lakedepth) - + !NoahMP parameters call NetCDF_read_var(grp_ncid, "tvxy", .False., input_tvxy) call NetCDF_read_var(grp_ncid, "tgxy", .False., input_tgxy) @@ -637,7 +636,7 @@ subroutine get_case_init(scm_state, scm_input) call NetCDF_read_var(grp_ncid, "deeprechxy",.False., input_deeprechxy) call NetCDF_read_var(grp_ncid, "rechxy", .False., input_rechxy) call NetCDF_read_var(grp_ncid, "snowxy", .False., input_snowxy) - + !NSST variables call NetCDF_read_var(grp_ncid, "tref", .False., input_tref) call NetCDF_read_var(grp_ncid, "z_c", .False., input_z_c) @@ -657,7 +656,7 @@ subroutine get_case_init(scm_state, scm_input) call NetCDF_read_var(grp_ncid, "ifd", .False., input_ifd) call NetCDF_read_var(grp_ncid, "dt_cool", .False., input_dt_cool) call NetCDF_read_var(grp_ncid, "qrain", .False., input_qrain) - + !RUC LSM variables call NetCDF_read_var(grp_ncid, "wetness", .False., input_wetness) call NetCDF_read_var(grp_ncid, "clw_surf_land", .False., input_clw_surf_land) @@ -673,7 +672,7 @@ subroutine get_case_init(scm_state, scm_input) call NetCDF_read_var(grp_ncid, "sfalb_lnd_bck", .False., input_sfalb_lnd_bck) call NetCDF_read_var(grp_ncid, "emis_ice", .False., input_emis_ice) call NetCDF_read_var(grp_ncid, "lai", .False., input_lai) - + !> - Read in the forcing data. !> - Find group ncid for forcing group. @@ -696,7 +695,7 @@ subroutine get_case_init(scm_state, scm_input) call NetCDF_read_var(grp_ncid, "T_surf", .True., input_T_surf) call NetCDF_read_var(grp_ncid, "sh_flux_sfc", .False., input_sh_flux_sfc) call NetCDF_read_var(grp_ncid, "lh_flux_sfc", .False., input_lh_flux_sfc) - + call NetCDF_read_var(grp_ncid, "w_ls", .True., input_w_ls) call NetCDF_read_var(grp_ncid, "omega", .True., input_omega) call NetCDF_read_var(grp_ncid, "u_g", .True., input_u_g) @@ -715,7 +714,7 @@ subroutine get_case_init(scm_state, scm_input) call check(NF90_CLOSE(NCID=ncid),"nf90_close()") call scm_input%create(input_ntimes, input_nlev, input_nsoil, input_nsnow, input_nice) - + ! GJF already done in scm_input%create routine !scm_input%input_nlev = input_nlev !scm_input%input_ntimes = input_ntimes @@ -751,26 +750,26 @@ subroutine get_case_init(scm_state, scm_input) scm_input%input_T_nudge = input_T_nudge scm_input%input_thil_nudge = input_thil_nudge scm_input%input_qt_nudge = input_qt_nudge - - scm_input%input_stc = input_stc - scm_input%input_smc = input_smc - scm_input%input_slc = input_slc - + + scm_input%input_stc = input_stc + scm_input%input_smc = input_smc + scm_input%input_slc = input_slc + scm_input%input_snicexy = input_snicexy scm_input%input_snliqxy = input_snliqxy scm_input%input_tsnoxy = input_tsnoxy scm_input%input_smoiseq = input_smoiseq scm_input%input_zsnsoxy = input_zsnsoxy - + scm_input%input_tiice = input_tiice scm_input%input_tslb = input_tslb scm_input%input_smois = input_smois scm_input%input_sh2o = input_sh2o scm_input%input_smfr = input_smfr scm_input%input_flfr = input_flfr - + scm_input%input_vegsrc = input_vegsrc - + scm_input%input_slmsk = input_slmsk scm_input%input_canopy = input_canopy scm_input%input_hice = input_hice @@ -831,7 +830,7 @@ subroutine get_case_init(scm_state, scm_input) scm_input%input_albdifvis_ice = input_albdifvis_ice scm_input%input_albdifnir_ice = input_albdifnir_ice scm_input%input_zorlwav = input_zorlwav - + scm_input%input_stddev = input_stddev scm_input%input_convexity= input_convexity scm_input%input_oa1 = input_oa1 @@ -851,7 +850,7 @@ subroutine get_case_init(scm_state, scm_input) scm_input%input_landfrac = input_landfrac scm_input%input_lakefrac = input_lakefrac scm_input%input_lakedepth= input_lakedepth - + scm_input%input_tvxy = input_tvxy scm_input%input_tgxy = input_tgxy scm_input%input_tahxy = input_tahxy @@ -881,7 +880,7 @@ subroutine get_case_init(scm_state, scm_input) scm_input%input_deeprechxy = input_deeprechxy scm_input%input_rechxy = input_rechxy scm_input%input_snowxy = input_snowxy - + scm_input%input_tref = input_tref scm_input%input_z_c = input_z_c scm_input%input_c_0 = input_c_0 @@ -900,9 +899,9 @@ subroutine get_case_init(scm_state, scm_input) scm_input%input_ifd = input_ifd scm_input%input_dt_cool = input_dt_cool scm_input%input_qrain = input_qrain - + scm_input%input_area = input_area - + scm_input%input_wetness = input_wetness scm_input%input_clw_surf_land = input_clw_surf_land scm_input%input_clw_surf_ice = input_clw_surf_ice @@ -917,29 +916,29 @@ subroutine get_case_init(scm_state, scm_input) scm_input%input_sfalb_lnd_bck = input_sfalb_lnd_bck scm_input%input_emis_ice = input_emis_ice scm_input%input_lai = input_lai - + !> @} end subroutine get_case_init subroutine get_case_init_DEPHY(scm_state, scm_input) !corresponds to the DEPHY-SCM specs, version 1 - + use scm_type_defs, only : scm_state_type, scm_input_type use NetCDF_read, only: NetCDF_read_var, NetCDF_read_att, NetCDF_conditionally_read_var, check, missing_value, missing_value_int use scm_physical_constants, only: con_hvap, con_hfus, con_cp, con_rocp, con_rd use scm_utils, only: find_vertical_index_pressure, find_vertical_index_height - + type(scm_state_type), intent(inout) :: scm_state type(scm_input_type), target, intent(inout) :: scm_input - + ! dimension variables real(kind=dp), allocatable :: input_t0(:) !< input initialization times (seconds since global attribute "startdate") real(kind=dp), allocatable :: input_time(:) !< input forcing times (seconds since the beginning of the case) real(kind=dp), allocatable :: input_lev(:) !< corresponds to either pressure or height (depending on attribute) - why is this needed when both pressure and height also provided in ICs? - + !non-standard dimensions (may or may not exist in the file) real(kind=dp), allocatable :: input_soil(:) !< soil depth - + ! global attributes character(len=19) :: char_startDate, char_endDate !format YYYY-MM-DD HH:MM:SS integer :: init_year, init_month, init_day, init_hour, init_min, init_sec, end_year, end_month, end_day, end_hour, end_min, end_sec @@ -951,7 +950,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) real(kind=sp) :: p_nudging_temp, p_nudging_theta, p_nudging_thetal, p_nudging_qv, p_nudging_qt, p_nudging_rv, p_nudging_rt, p_nudging_u, p_nudging_v character(len=5) :: input_surfaceType character(len=12) :: input_surfaceForcingWind='',input_surfaceForcingMoist='',input_surfaceForcingLSM='',input_surfaceForcingTemp='' - + ! initial variables (IC = Initial Condition) real(kind=dp), allocatable :: input_lat(:) !< column latitude (deg) real(kind=dp), allocatable :: input_lon(:) !< column longitude (deg) @@ -974,11 +973,11 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) real(kind=sp), allocatable :: input_ri(:,:) !< IC ice water mixing ratio profile (kg kg^-1) real(kind=sp), allocatable :: input_rh(:,:) !< IC relative humidity profile (%) real(kind=sp), allocatable :: input_tke(:,:) !< IC TKE profile (m^2 s^-2) - + ! Model ICs (extension of DEPHY format) real(kind=dp), allocatable :: input_ozone(:,:) !< ozone profile (kg kg^-1) real(kind=dp), allocatable :: input_stc(:,:) !< soil temperature (K) - real(kind=dp), allocatable :: input_smc(:,:) !< total soil moisture content (fraction) + real(kind=dp), allocatable :: input_smc(:,:) !< total soil moisture content (fraction) real(kind=dp), allocatable :: input_slc(:,:) !< liquid soil moisture content (fraction) real(kind=dp), allocatable :: input_snicexy(:,:) !< snow layer ice (mm) real(kind=dp), allocatable :: input_snliqxy(:,:) !< snow layer liquid (mm) @@ -991,7 +990,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) real(kind=dp), allocatable :: input_sh2o(:,:) !< volume fraction of unfrozen soil moisture for RUC LSM (frac) real(kind=dp), allocatable :: input_smfr(:,:) !< volume fraction of frozen soil moisture for RUC LSM (frac) real(kind=dp), allocatable :: input_flfr(:,:) !< flag for frozen soil physics - + real(kind=dp), allocatable :: input_area(:) !< surface area [m^2] real(kind=dp), allocatable :: input_tsfco(:) !< input sea surface temperature OR surface skin temperature over land OR surface skin temperature over ice (depending on slmsk) (K) integer , allocatable :: input_vegsrc(:) !< vegetation source @@ -1030,9 +1029,9 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) real(kind=dp), allocatable :: input_zorll(:) !< surface roughness length over land (cm) real(kind=dp), allocatable :: input_zorli(:) !< surface roughness length over ice (cm) real(kind=dp), allocatable :: input_zorlw(:) !< surface roughness length from wave model (cm) - + real(kind=dp), allocatable :: input_stddev(:) !< standard deviation of subgrid orography (m) - real(kind=dp), allocatable :: input_convexity(:) !< convexity of subgrid orography + real(kind=dp), allocatable :: input_convexity(:) !< convexity of subgrid orography real(kind=dp), allocatable :: input_ol1(:) !< fraction of grid box with subgrid orography higher than critical height 1 real(kind=dp), allocatable :: input_ol2(:) !< fraction of grid box with subgrid orography higher than critical height 2 real(kind=dp), allocatable :: input_ol3(:) !< fraction of grid box with subgrid orography higher than critical height 3 @@ -1050,7 +1049,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) real(kind=dp), allocatable :: input_landfrac(:) !< fraction of horizontal grid area occupied by land real(kind=dp), allocatable :: input_lakefrac(:) !< fraction of horizontal grid area occupied by lake real(kind=dp), allocatable :: input_lakedepth(:) !< lake depth (m) - + real(kind=dp), allocatable :: input_tvxy(:) !< vegetation temperature (K) real(kind=dp), allocatable :: input_tgxy(:) !< ground temperature for Noahmp (K) real(kind=dp), allocatable :: input_tahxy(:) !< canopy air temperature (K) @@ -1080,7 +1079,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) real(kind=dp), allocatable :: input_deeprechxy(:) !< recharge to or from the water table when deep (m) real(kind=dp), allocatable :: input_rechxy(:) !< recharge to or from the water table when shallow (m) real(kind=dp), allocatable :: input_snowxy(:) !< number of snow layers - + real(kind=dp), allocatable :: input_tref(:) !< sea surface reference temperature for NSST (K) real(kind=dp), allocatable :: input_z_c(:) !< sub-layer cooling thickness for NSST (m) real(kind=dp), allocatable :: input_c_0(:) !< coefficient 1 to calculate d(Tz)/d(Ts) for NSST @@ -1099,23 +1098,23 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) real(kind=dp), allocatable :: input_ifd(:) !< index to start DTM run for NSST real(kind=dp), allocatable :: input_dt_cool(:) !< sub-layer cooling amount for NSST (K) real(kind=dp), allocatable :: input_qrain(:) !< sensible heat due to rainfall for NSST (W) - + real(kind=dp), allocatable :: input_wetness(:) !< normalized soil wetness for RUC LSM real(kind=dp), allocatable :: input_lai(:) !< leaf area index for RUC LSM - real(kind=dp), allocatable :: input_clw_surf_land(:) !< cloud condensed water mixing ratio at surface over land for RUC LSM + real(kind=dp), allocatable :: input_clw_surf_land(:) !< cloud condensed water mixing ratio at surface over land for RUC LSM real(kind=dp), allocatable :: input_clw_surf_ice(:) !< cloud condensed water mixing ratio at surface over ice for RUC LSM - real(kind=dp), allocatable :: input_qwv_surf_land(:) !< water vapor mixing ratio at surface over land for RUC LSM - real(kind=dp), allocatable :: input_qwv_surf_ice(:) !< water vapor mixing ratio at surface over ice for RUC LSM - real(kind=dp), allocatable :: input_tsnow_land(:) !< snow temperature at the bottom of the first snow layer over land for RUC LSM - real(kind=dp), allocatable :: input_tsnow_ice(:) !< snow temperature at the bottom of the first snow layer over ice for RUC LSM - real(kind=dp), allocatable :: input_snowfallac_land(:) !< run-total snow accumulation on the ground over land for RUC LSM - real(kind=dp), allocatable :: input_snowfallac_ice(:) !< run-total snow accumulation on the ground over ice for RUC LSM + real(kind=dp), allocatable :: input_qwv_surf_land(:) !< water vapor mixing ratio at surface over land for RUC LSM + real(kind=dp), allocatable :: input_qwv_surf_ice(:) !< water vapor mixing ratio at surface over ice for RUC LSM + real(kind=dp), allocatable :: input_tsnow_land(:) !< snow temperature at the bottom of the first snow layer over land for RUC LSM + real(kind=dp), allocatable :: input_tsnow_ice(:) !< snow temperature at the bottom of the first snow layer over ice for RUC LSM + real(kind=dp), allocatable :: input_snowfallac_land(:) !< run-total snow accumulation on the ground over land for RUC LSM + real(kind=dp), allocatable :: input_snowfallac_ice(:) !< run-total snow accumulation on the ground over ice for RUC LSM real(kind=dp), allocatable :: input_sncovr_ice(:) !< real(kind=dp), allocatable :: input_sfalb_lnd(:) !< real(kind=dp), allocatable :: input_sfalb_lnd_bck(:) !< real(kind=dp), allocatable :: input_sfalb_ice(:) !< real(kind=dp), allocatable :: input_emis_ice(:) !< - + ! forcing variables real(kind=sp), allocatable :: input_force_pres_surf(:) !< forcing surface pressure (Pa) real(kind=sp), allocatable :: input_force_height(:,:) !< forcing height levels (m) @@ -1154,7 +1153,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) real(kind=sp), allocatable :: input_force_rt_nudging(:,:) real(kind=sp), allocatable :: input_force_u_nudging(:,:) real(kind=sp), allocatable :: input_force_v_nudging(:,:) - + integer :: ncid, varID, allocate_status, ierr, i, k integer :: active_lon, active_lat, active_init_time CHARACTER(LEN=nf90_max_name) :: tmpName @@ -1164,14 +1163,14 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) integer :: jdat(1:8), idat(1:8) !(yr, mon, day, t-zone, hr, min, sec, mil-sec) integer :: input_n_init_times, input_n_forcing_times, input_n_lev, input_n_snow, input_n_ice, input_n_soil - + missing_value_eps = missing_value + 0.01 - + !> - Open the case input file found in the processed_case_input dir corresponding to the experiment name. call check(NF90_OPEN(trim(adjustl(scm_state%case_name))//'_SCM_driver.nc',nf90_nowrite,ncid),"nf90_open()") - + !> - Get the dimensions. - + !required dimensions call check(NF90_INQ_DIMID(ncid,"t0",varID),"nf90_inq_dimid(t0)") call check(NF90_INQUIRE_DIMENSION(ncid, varID, tmpName, input_n_init_times),"nf90_inq_dim(t0)") @@ -1201,8 +1200,8 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) input_n_ice = missing_ice_layers else call check(NF90_INQUIRE_DIMENSION(ncid, varID, tmpName, input_n_ice),"nf90_inq_dim(nice)") - end if - + end if + !> - Allocate the dimension variables. allocate(input_t0 (input_n_init_times), & input_time (input_n_forcing_times), & @@ -1217,23 +1216,23 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) !> - Read in global attributes call NetCDF_read_att(ncid, NF90_GLOBAL, 'start_date', .True., char_startDate) - + read(char_startDate(1:4),'(i4)') init_year read(char_startDate(6:7),'(i2)') init_month read(char_startDate(9:10),'(i2)') init_day read(char_startDate(12:13),'(i2)') init_hour read(char_startDate(15:16),'(i2)') init_min read(char_startDate(18:19),'(i2)') init_sec - + call NetCDF_read_att(ncid, NF90_GLOBAL, 'end_date', .True., char_endDate) - + read(char_endDate(1:4),'(i4)') end_year read(char_endDate(6:7),'(i2)') end_month read(char_endDate(9:10),'(i2)') end_day read(char_endDate(12:13),'(i2)') end_hour read(char_endDate(15:16),'(i2)') end_min read(char_endDate(18:19),'(i2)') end_sec - + !compare init time to what was in case config file? replace? call NetCDF_read_att(ncid, NF90_GLOBAL, 'adv_ua', .False., adv_u) call NetCDF_read_att(ncid, NF90_GLOBAL, 'adv_va', .False., adv_v) @@ -1325,8 +1324,8 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) input_smfr (input_n_soil, input_n_init_times), & input_flfr (input_n_soil, input_n_init_times), & stat=allocate_status) - - + + !variables without vertical extent allocate(input_area ( input_n_init_times), & input_tsfco ( input_n_init_times), & @@ -1388,7 +1387,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) input_lakedepth ( input_n_init_times), & stat=allocate_status) allocate(input_tvxy ( input_n_init_times), & - input_tgxy ( input_n_init_times), & + input_tgxy ( input_n_init_times), & input_tahxy ( input_n_init_times), & input_canicexy ( input_n_init_times), & input_canliqxy ( input_n_init_times), & @@ -1414,8 +1413,8 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) input_fastcpxy ( input_n_init_times), & input_smcwtdxy ( input_n_init_times), & input_deeprechxy( input_n_init_times), & - input_rechxy ( input_n_init_times), & - input_snowxy ( input_n_init_times), & + input_rechxy ( input_n_init_times), & + input_snowxy ( input_n_init_times), & stat=allocate_status) allocate(input_tref ( input_n_init_times), & input_z_c ( input_n_init_times), & @@ -1460,12 +1459,12 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) call NetCDF_read_var(ncid, "ps", .True., input_pres_surf) call NetCDF_read_var(ncid, "ua", .True., input_u) call NetCDF_read_var(ncid, "va", .True., input_v) - + !one of the following should be present, but not all, hence they are not requried call NetCDF_read_var(ncid, "ta", .False., input_temp) call NetCDF_read_var(ncid, "theta", .False., input_theta) call NetCDF_read_var(ncid, "thetal", .False., input_thetal) - + !one or more of the following should be present, but not all, hence they are not requried call NetCDF_read_var(ncid, "qv", .False., input_qv) call NetCDF_read_var(ncid, "qt", .False., input_qt) @@ -1476,13 +1475,13 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) call NetCDF_read_var(ncid, "rl", .False., input_rl) call NetCDF_read_var(ncid, "ri", .False., input_ri) call NetCDF_read_var(ncid, "hur", .False., input_rh) - + call NetCDF_read_var(ncid, "tke", .True., input_tke) - + if (trim(input_surfaceForcingLSM) == "lsm") then call NetCDF_read_var(ncid, "o3", .True., input_ozone) call NetCDF_read_var(ncid, "area", .True., input_area) - + !orographic parameters call NetCDF_read_var(ncid, "stddev", .True., input_stddev) call NetCDF_read_var(ncid, "convexity", .True., input_convexity) @@ -1503,7 +1502,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) call NetCDF_read_var(ncid, "landfrac", .True., input_landfrac) call NetCDF_read_var(ncid, "lakefrac", .True., input_lakefrac) call NetCDF_read_var(ncid, "lakedepth", .True., input_lakedepth) - + !NSST variables call NetCDF_read_var(ncid, "tref", .True., input_tref) call NetCDF_read_var(ncid, "z_c", .True., input_z_c) @@ -1524,9 +1523,9 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) call NetCDF_read_var(ncid, "dt_cool", .True., input_dt_cool) call NetCDF_read_var(ncid, "qrain", .True., input_qrain) end if - + !> - Allocate the forcing variables. - + !allocate all, but conditionally read forcing variables given global atts; set unused forcing variables to missing allocate(input_lat (input_n_forcing_times), & @@ -1572,11 +1571,11 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) stat=allocate_status) call NetCDF_read_var(ncid, "lat", .True., input_lat) - call NetCDF_read_var(ncid, "lon", .True., input_lon) + call NetCDF_read_var(ncid, "lon", .True., input_lon) call NetCDF_read_var(ncid, "ps_forc", .True., input_force_pres_surf) call NetCDF_read_var(ncid, "zh_forc", .True., input_force_height) call NetCDF_read_var(ncid, "pa_forc", .True., input_force_pres) - + !conditionally read forcing vars (or set to missing); if the global attribute is set to expect a variable and it doesn't exist, stop the model call NetCDF_conditionally_read_var(adv_u, "adv_ua", "tnua_adv", trim(adjustl(scm_state%case_name))//'.nc', ncid, input_force_u_adv) call NetCDF_conditionally_read_var(adv_v, "adv_va", "tnva_adv", trim(adjustl(scm_state%case_name))//'.nc', ncid, input_force_v_adv) @@ -1592,12 +1591,12 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) call NetCDF_conditionally_read_var(rad_thetal, "rad_thetal", "tnthetal_rad", trim(adjustl(scm_state%case_name))//'.nc', ncid, input_force_thetal_rad) !need to also handle the case when rad_[temp,theta,thetal]_char = 'adv' (make sure [temp,theta,thetal]_adv is not missing) !need to also turn off radiation when radiation is being forced (put in a warning that this is not supported for now?) - + call NetCDF_conditionally_read_var(forc_w, "forc_w", "wa", trim(adjustl(scm_state%case_name))//'.nc', ncid, input_force_w) call NetCDF_conditionally_read_var(forc_omega, "forc_wap", "wap", trim(adjustl(scm_state%case_name))//'.nc', ncid, input_force_omega) call NetCDF_conditionally_read_var(forc_geo, "forc_geo", "ug", trim(adjustl(scm_state%case_name))//'.nc', ncid, input_force_u_g) call NetCDF_conditionally_read_var(forc_geo, "forc_geo", "vg", trim(adjustl(scm_state%case_name))//'.nc', ncid, input_force_v_g) - + call NetCDF_conditionally_read_var(nudging_u, "nudging_u", "ua_nud", trim(adjustl(scm_state%case_name))//'.nc', ncid, input_force_u_nudging) call NetCDF_conditionally_read_var(nudging_v, "nudging_v", "va_nud", trim(adjustl(scm_state%case_name))//'.nc', ncid, input_force_v_nudging) call NetCDF_conditionally_read_var(nudging_temp, "nudging_temp", "ta_nud", trim(adjustl(scm_state%case_name))//'.nc', ncid, input_force_temp_nudging) @@ -1660,7 +1659,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) call NetCDF_read_var(ncid, "sh2o", .True., input_sh2o ) call NetCDF_read_var(ncid, "smfr", .True., input_smfr ) call NetCDF_read_var(ncid, "flfr", .True., input_flfr ) - + call NetCDF_read_var(ncid, "vegsrc", .True., input_vegsrc ) call NetCDF_read_var(ncid, "vegtyp", .True., input_vegtyp ) call NetCDF_read_var(ncid, "soiltyp", .True., input_soiltyp ) @@ -1698,7 +1697,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) call NetCDF_read_var(ncid, "zorll", .True., input_zorll) call NetCDF_read_var(ncid, "zorli", .True., input_zorli) call NetCDF_read_var(ncid, "zorlw", .True., input_zorlw) - + !NoahMP parameters call NetCDF_read_var(ncid, "tvxy", .False., input_tvxy) call NetCDF_read_var(ncid, "tgxy", .False., input_tgxy) @@ -1745,18 +1744,18 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) call NetCDF_read_var(ncid, "emis_ice", .False., input_emis_ice) call NetCDF_read_var(ncid, "lai", .False., input_lai) end if - + call check(NF90_CLOSE(NCID=ncid),"nf90_close()") - + call scm_input%create(input_n_forcing_times, input_n_lev, input_n_soil, input_n_snow, input_n_ice) - + !fill the scm_input DDT - + !There may need to be logic to control which of the lon, lat, and init_times to use in the future, but for now, just take the first active_lon = 1 active_lat = 1 active_init_time = 1 - + rinc(1:5) = 0 idat = 0 jdat = 0 @@ -1774,7 +1773,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) jdat(7) = end_sec call w3difdat(jdat,idat,4,rinc) elapsed_sec = rinc(4) - + !the following variables replace what is in the case configuration file scm_state%init_year = init_year scm_state%init_month = init_month @@ -1782,14 +1781,14 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_state%init_hour = init_hour scm_state%init_min = init_min scm_state%runtime = elapsed_sec - + scm_input%input_time = input_time scm_input%input_pres_surf(1) = input_pres_surf(active_init_time) !perhaps input_pres_surf should only be equal to input_force_pres_surf? scm_input%input_pres = input_pres(:,active_init_time) scm_input%input_u = input_u(:,active_init_time) scm_input%input_v = input_v(:,active_init_time) scm_input%input_tke = input_tke(:,active_init_time) - + !if mixing ratios are present, and not specific humidities, convert from mixing ratio to specific humidities if ((maxval(input_qv(:,active_init_time)) < 0 .and. & maxval(input_qt(:,active_init_time)) < 0) .and. & @@ -1820,7 +1819,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) end do end if end if - + !make sure that one of qv or qt (and rv or rt due to above conversion) is present (add support for rh later) if (maxval(input_qv(:,active_init_time)) >= 0) then if (maxval(input_qt(:,active_init_time)) >= 0) then @@ -1840,7 +1839,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_qi(k) = max(0.0, scm_input%input_qt(k) - scm_input%input_qv(k) - scm_input%input_ql(k)) end do end if !qi test - else + else if (maxval(input_qi(:,active_init_time)) >= 0) then !qv, qt, qi, but no ql scm_input%input_qv = input_qv(:,active_init_time) scm_input%input_qt = input_qt(:,active_init_time) @@ -1876,7 +1875,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) do k=1, input_n_lev scm_input%input_qt(k) = max(0.0, scm_input%input_qv(k) + scm_input%input_ql(k)) end do - scm_input%input_qi = 0.0 + scm_input%input_qi = 0.0 end if ! qi test else if (maxval(input_qi(:,active_init_time)) >= 0) then !qv, no qt, no ql, qi @@ -1933,9 +1932,9 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) else !no qv or qt write(*,*) 'When reading '//trim(adjustl(scm_state%case_name))//'.nc, all of the supported moisture variables (qv, qt, rv, rt) were missing. Stopping...' - stop + error stop "Aall of the supported moisture variables (qv, qt, rv, rt) were missing" end if - + !make sure that at least one of temp, theta, thetal is present; !the priority for use is temp, thetal, theta if (maxval(input_temp(:,active_init_time)) > 0) then @@ -1963,13 +1962,13 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_temp = missing_value else write(*,*) 'When reading '//trim(adjustl(scm_state%case_name))//'.nc, all of the supported temperature variables (temp, theta, thetal) were missing. Stopping...' - stop + error stop "All of the supported temperature variables (temp, theta, thetal) were missing" end if if (trim(input_surfaceForcingLSM) == "lsm") then scm_input%input_ozone = input_ozone(:,active_init_time) scm_input%input_area = input_area(active_init_time) - + scm_input%input_stddev = input_stddev(active_init_time) scm_input%input_convexity= input_convexity(active_init_time) scm_input%input_oa1 = input_oa1(active_init_time) @@ -1989,7 +1988,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_landfrac = input_landfrac(active_init_time) scm_input%input_lakefrac = input_lakefrac(active_init_time) scm_input%input_lakedepth= input_lakedepth(active_init_time) - + scm_input%input_tref = input_tref(active_init_time) scm_input%input_z_c = input_z_c(active_init_time) scm_input%input_c_0 = input_c_0(active_init_time) @@ -2014,24 +2013,24 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) end if scm_input%input_lat = input_lat(active_lat) scm_input%input_lon = input_lon(active_lon) - + scm_input%input_pres_surf = input_force_pres_surf(:) - + do i=1, input_n_forcing_times scm_input%input_pres_forcing(i,:) = input_force_pres(:,i) end do - + if (input_SurfaceType == 'ocean') then scm_state%sfc_type = 0.0 else if (input_SurfaceType == 'land') then scm_state%sfc_type = 1.0 end if !no sea ice type? - + if (input_surfaceForcingTemp == 'ts') then if (maxval(input_force_ts) < 0) then write(*,*) 'The global attribute surfaceForcing in '//trim(adjustl(scm_state%case_name))//'.nc indicates that the variable ts should be present, but it is missing. Stopping ...' - stop + error stop "The global attribute surfaceForcing indicates that the variable ts should be present, but it is missing" else !overwrite sfc_flux_spec scm_state%sfc_flux_spec = .false. @@ -2042,7 +2041,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) !overwrite sfc_flux_spec scm_state%sfc_flux_spec = .true. scm_state%surface_thermo_control = 0 - + if (maxval(input_force_ts) < 0) then !since no surface temperature is given, assume that the surface temperature is equivalent to the static, surface-adjacent temperature in the initial profile if (maxval(scm_input%input_temp) > 0) then @@ -2059,11 +2058,11 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) else scm_input%input_T_surf = input_force_ts(:) end if - + !kinematic surface fluxes are specified (but may need to be converted) if (maxval(input_force_wpthetap(:)) < missing_value_eps) then write(*,*) 'The global attribute surfaceForcing in '//trim(adjustl(scm_state%case_name))//'.nc indicates that the variable wpthetap should be present, but it is missing. Stopping ...' - stop + error stop "The global attribute surfaceForcing indicates that the variable wpthetap should be present, but it is missing." else !convert from theta to T do i=1, input_n_forcing_times @@ -2071,7 +2070,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_sh_flux_sfc_kin(i) = exner*input_force_wpthetap(i) end do end if - + !if mixing ratios are present, and not specific humidities, convert from mixing ratio to specific humidities if ((maxval(input_force_wpqvp(:)) < missing_value_eps .and. & maxval(input_force_wpqtp(:)) < missing_value_eps) .and. & @@ -2090,10 +2089,10 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) end do end if end if - + if (maxval(input_force_wpqvp(:)) < missing_value_eps .and. maxval(input_force_wpqtp(:)) < missing_value_eps) then write(*,*) 'The global attribute surfaceForcing in '//trim(adjustl(scm_state%case_name))//'.nc indicates that the variable wpqvp, wpqtp, wprvp, or wprtp should be present, but all are missing. Stopping ...' - stop + error stop "The global attribute surfaceForcing indicates that the variable wpqvp, wpqtp, wprvp, or wprtp should be present, but all are missing." else if (maxval(input_force_wpqvp(:)) > missing_value_eps) then !use wpqvp if available scm_input%input_lh_flux_sfc_kin = input_force_wpqvp(:) @@ -2106,7 +2105,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) !overwrite sfc_flux_spec scm_state%sfc_flux_spec = .true. scm_state%surface_thermo_control = 1 - + if (maxval(input_force_ts) < 0) then !since no surface temperature is given, assume that the surface temperature is equivalent to the static, surface-adjacent temperature in the initial profile if (maxval(scm_input%input_temp) > 0) then @@ -2123,40 +2122,40 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) else scm_input%input_T_surf = input_force_ts(:) end if - - + + if (maxval(input_force_sfc_sens_flx(:)) < missing_value_eps) then write(*,*) 'The global attribute surfaceForcing in '//trim(adjustl(scm_state%case_name))//'.nc indicates that the variable sfc_sens_flx should be present, but it is missing. Stopping ...' - stop + error stop "The global attribute surfaceForcing in indicates that the variable sfc_sens_flx should be present, but it is missing." else scm_input%input_sh_flux_sfc = input_force_sfc_sens_flx(:) end if - + if (maxval(input_force_sfc_lat_flx(:)) < missing_value_eps) then write(*,*) 'The global attribute surfaceForcing in '//trim(adjustl(scm_state%case_name))//'.nc indicates that the variable sfc_lat_flx should be present, but it is missing. Stopping ...' - stop + error stop "The global attribute surfaceForcing indicates that the variable sfc_lat_flx should be present, but it is missing." else scm_input%input_lh_flux_sfc = input_force_sfc_lat_flx(:) end if else if (trim(input_surfaceForcingLSM) == 'lsm') then !these were considered required variables above, so they should not need to be checked for missing scm_input%input_stc = input_stc(:,active_init_time) - scm_input%input_smc = input_smc(:,active_init_time) - scm_input%input_slc = input_slc(:,active_init_time) - + scm_input%input_smc = input_smc(:,active_init_time) + scm_input%input_slc = input_slc(:,active_init_time) + scm_input%input_snicexy = input_snicexy(:,active_init_time) scm_input%input_snliqxy = input_snliqxy(:,active_init_time) scm_input%input_tsnoxy = input_tsnoxy(:,active_init_time) scm_input%input_smoiseq = input_smoiseq(:,active_init_time) scm_input%input_zsnsoxy = input_zsnsoxy(:,active_init_time) - + scm_input%input_tiice = input_tiice(:,active_init_time) scm_input%input_tslb = input_tslb(:,active_init_time) scm_input%input_smois = input_smois(:,active_init_time) scm_input%input_sh2o = input_sh2o(:,active_init_time) scm_input%input_smfr = input_smfr(:,active_init_time) scm_input%input_flfr = input_flfr(:,active_init_time) - + scm_input%input_vegsrc = input_vegsrc(active_init_time) scm_input%input_vegtyp = REAL(input_vegtyp(active_init_time), kind=dp) scm_input%input_soiltyp = REAL(input_soiltyp(active_init_time), kind=dp) @@ -2194,7 +2193,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_zorll = input_zorll(active_init_time) scm_input%input_zorli = input_zorli(active_init_time) scm_input%input_zorlw = input_zorlw(active_init_time) - + scm_input%input_tvxy = input_tvxy(active_init_time) scm_input%input_tgxy = input_tgxy(active_init_time) scm_input%input_tahxy = input_tahxy(active_init_time) @@ -2239,7 +2238,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_sfalb_lnd_bck = input_sfalb_lnd_bck(active_init_time) scm_input%input_emis_ice = input_emis_ice(active_init_time) end if - + if (input_surfaceForcingWind == 'z0') then scm_state%surface_momentum_control = 0 scm_state%sfc_roughness_length_cm = input_z0*100.0 !convert from m to cm @@ -2247,9 +2246,9 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) !not supported scm_state%surface_momentum_control = 1 write(*,*) 'The global attribute surfaceForcingWind in '//trim(adjustl(scm_state%case_name))//'.nc indicates that surface wind is controlled by a specified time-series of ustar. This is currently not supported. Stopping ...' - stop + error stop "The global attribute surfaceForcingWind indicates that surface wind is controlled by a specified time-series of ustar. This is currently not supported." end if - + if (forc_omega > 0) then do i=1, input_n_forcing_times scm_input%input_omega(i,:) = input_force_omega(:,i) @@ -2281,7 +2280,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) end do scm_state%force_geo = .true. end if - + if (adv_temp > 0) then do i=1, input_n_forcing_times scm_input%input_tot_advec_T(i,:) = input_force_temp_adv(:,i) @@ -2298,7 +2297,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) end do scm_state%force_adv_T = 3 end if - + if (adv_qv > 0) then do i=1, input_n_forcing_times scm_input%input_tot_advec_qv(i,:) = input_force_qv_adv(:,i) @@ -2330,26 +2329,26 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) end if scm_state%force_adv_qv = .true. end if - + if (adv_u > 0) then do i=1, input_n_forcing_times scm_input%input_tot_advec_u(i,:) = input_force_u_adv(:,i) end do scm_state%force_adv_u = .true. end if - + if (adv_v > 0) then do i=1, input_n_forcing_times scm_input%input_tot_advec_v(i,:) = input_force_v_adv(:,i) end do scm_state%force_adv_v = .true. end if - + if (char_rad_temp == 'adv' .or. char_rad_theta == 'adv' .or. char_rad_thetal == 'adv') then scm_state%force_rad_T = 4 if (scm_state%force_adv_T == 0) then write(*,*) 'The global attribute rad_temp, rad_theta, or rad_thetal in '//trim(adjustl(scm_state%case_name))//'.nc indicates that radiative forcing is included in the advection term, but there is no advection term. Stopping ...' - stop + error stop "The global attribute rad_temp, rad_theta, or rad_thetal indicates that radiative forcing is included in the advection term, but there is no advection term." end if else if (rad_temp > 0) then do i=1, input_n_forcing_times @@ -2375,7 +2374,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) else scm_state%force_rad_T = 0 end if - + if (nudging_temp > 0) then do i=1, input_n_forcing_times scm_input%input_T_nudge(i,:) = input_force_temp_nudging(:,i) @@ -2390,7 +2389,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_k_T_nudge(i) = input_n_lev end if end do - else if (z_nudging_temp > 0) then + else if (z_nudging_temp > 0) then do i=1, input_n_forcing_times call find_vertical_index_height(z_nudging_temp, input_force_height(:,i), scm_input%input_k_T_nudge(i)) if (scm_input%input_k_T_nudge(i) < 0) then @@ -2416,7 +2415,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_k_thil_nudge(i) = input_n_lev end if end do - else if (z_nudging_theta > 0) then + else if (z_nudging_theta > 0) then do i=1, input_n_forcing_times call find_vertical_index_height(z_nudging_theta, input_force_height(:,i), scm_input%input_k_thil_nudge(i)) if (scm_input%input_k_thil_nudge(i) < 0) then @@ -2442,7 +2441,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_k_thil_nudge(i) = input_n_lev end if end do - else if (z_nudging_thetal > 0) then + else if (z_nudging_thetal > 0) then do i=1, input_n_forcing_times call find_vertical_index_height(z_nudging_thetal, input_force_height(:,i), scm_input%input_k_thil_nudge(i)) if (scm_input%input_k_thil_nudge(i) < 0) then @@ -2454,7 +2453,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_k_thil_nudge = 1 end if end if - + if (nudging_qv > 0) then do i=1, input_n_forcing_times scm_input%input_qt_nudge(i,:) = input_force_qv_nudging(:,i) @@ -2469,7 +2468,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_k_qt_nudge(i) = input_n_lev end if end do - else if (z_nudging_qv > 0) then + else if (z_nudging_qv > 0) then do i=1, input_n_forcing_times call find_vertical_index_height(z_nudging_qv, input_force_height(:,i), scm_input%input_k_qt_nudge(i)) if (scm_input%input_k_qt_nudge(i) < 0) then @@ -2494,7 +2493,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_k_qt_nudge(i) = input_n_lev end if end do - else if (z_nudging_qt > 0) then + else if (z_nudging_qt > 0) then do i=1, input_n_forcing_times call find_vertical_index_height(z_nudging_qt, input_force_height(:,i), scm_input%input_k_qt_nudge(i)) if (scm_input%input_k_qt_nudge(i) < 0) then @@ -2522,7 +2521,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_k_qt_nudge(i) = input_n_lev end if end do - else if (z_nudging_rv > 0) then + else if (z_nudging_rv > 0) then do i=1, input_n_forcing_times call find_vertical_index_height(z_nudging_rv, input_force_height(:,i), scm_input%input_k_qt_nudge(i)) if (scm_input%input_k_qt_nudge(i) < 0) then @@ -2550,7 +2549,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_k_qt_nudge(i) = input_n_lev end if end do - else if (z_nudging_rt > 0) then + else if (z_nudging_rt > 0) then do i=1, input_n_forcing_times call find_vertical_index_height(z_nudging_rt, input_force_height(:,i), scm_input%input_k_qt_nudge(i)) if (scm_input%input_k_qt_nudge(i) < 0) then @@ -2562,7 +2561,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_k_qt_nudge = 1 end if end if - + if (nudging_u > 0) then do i=1, input_n_forcing_times scm_input%input_u_nudge(i,:) = input_force_u_nudging(:,i) @@ -2577,8 +2576,8 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_k_u_nudge(i) = input_n_lev end if end do - else if (z_nudging_u > 0) then - do i=1, input_n_forcing_times + else if (z_nudging_u > 0) then + do i=1, input_n_forcing_times call find_vertical_index_height(z_nudging_u, input_force_height(:,i), scm_input%input_k_u_nudge(i)) if (scm_input%input_k_u_nudge(i) < 0) then !if the vertical index is not found (when it is less than 0), set the nudging index to the top of the input profile so that nudging is turned off @@ -2589,7 +2588,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_k_u_nudge = 1 end if end if - + if (nudging_v > 0) then do i=1, input_n_forcing_times scm_input%input_v_nudge(i,:) = input_force_v_nudging(:,i) @@ -2604,7 +2603,7 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_k_v_nudge(i) = input_n_lev end if end do - else if (z_nudging_v > 0) then + else if (z_nudging_v > 0) then do i=1, input_n_forcing_times call find_vertical_index_height(z_nudging_v, input_force_height(:,i), scm_input%input_k_v_nudge(i)) if (scm_input%input_k_v_nudge(i) < 0) then @@ -2616,14 +2615,14 @@ subroutine get_case_init_DEPHY(scm_state, scm_input) scm_input%input_k_v_nudge = 1 end if end if - + end subroutine get_case_init_DEPHY !> Subroutine to get reference profile to use above the case data (temporarily hard-coded profile) subroutine get_reference_profile(scm_state, scm_reference) use scm_type_defs, only : scm_state_type, scm_reference_type use NetCDF_read, only: check - + type(scm_state_type), target, intent(in) :: scm_state type(scm_reference_type), target, intent(inout) :: scm_reference @@ -2646,7 +2645,7 @@ subroutine get_reference_profile(scm_state, scm_reference) if(ioerror /= 0) then write(*,*) 'There was an error opening the file McCprofiles.dat in the processed_case_input directory. & Error code = ',ioerror - stop + error stop "There was an error opening the file McCprofiles.dat in the processed_case_input directory." endif ! find number of records @@ -2756,7 +2755,7 @@ subroutine get_tracers(tracer_names, tracer_types) else write(*,'(a,i0)') 'There was an error opening the file ' // FILE_NAME // & '; error code = ', rc - stop + error stop "Error opening tracers file" end if close (fu) diff --git a/scm/src/scm_setup.F90 b/scm/src/scm_setup.F90 index 5f4423e29..811844dc0 100644 --- a/scm/src/scm_setup.F90 +++ b/scm/src/scm_setup.F90 @@ -31,34 +31,34 @@ subroutine set_state(scm_input, scm_reference, scm_state) real(kind=dp), dimension(scm_input%input_nlev) :: input_qv, input_T real(kind=dp), parameter :: p0 = 1.0E5 real(kind=dp) :: deg_to_rad_const - + deg_to_rad_const = con_pi/180.0 !> \section set_state_alg Algorithm !! @{ - + !> - Set the longitude and latitude and convert from degrees to radians do i=1, scm_state%n_cols scm_state%lon(i) = scm_input%input_lon*deg_to_rad_const scm_state%lat(i) = scm_input%input_lat*deg_to_rad_const end do - - + + !> - \todo When patching in a reference sounding, need to handle the case when the reference sounding is too short; patch_in_ref !! checks for the case, but as of now, it just extrapolates where it needs to and returns an error code; error should be handled !! here or passed up to the main program. - + if (.NOT. scm_state%model_ics) then ! not a model - + if (scm_state%input_type == 0) then !> - Calculate water vapor from total water, suspended liquid water, and suspended ice. input_qv = scm_input%input_qt - scm_input%input_ql - scm_input%input_qi else input_qv = scm_input%input_qv end if - + !> - For each column, interpolate the water vapor to the model grid. do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres, input_qv, scm_state%pres_l(i,:), & @@ -70,7 +70,7 @@ subroutine set_state(scm_input, scm_reference, scm_state) scm_state%state_tracer(i,:,scm_state%water_vapor_index,1), grid_error) end if end do - + if (scm_state%input_type == 0) then !> - Calculate the input absolute temperature from input pressure, theta_il, ql, and qi. input_T = (scm_input%input_pres/p0)**con_rocp*(scm_input%input_thetail + (con_hvap/con_cp)*scm_input%input_ql + & @@ -83,7 +83,7 @@ subroutine set_state(scm_input, scm_reference, scm_state) (con_hfus/con_cp)*scm_input%input_qi) end if end if - + !> - For each column, interpolate the temperature to the model grid. do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres, input_T, scm_state%pres_l(i,:), & @@ -94,7 +94,7 @@ subroutine set_state(scm_input, scm_reference, scm_state) scm_reference%ref_T, scm_state%pres_l(i,:), scm_state%n_levels, scm_state%state_T(i,:,1), grid_error) end if end do - + !> - For each column, interpolate the u-wind to the model grid. do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres, scm_input%input_u, scm_state%pres_l(i,:), & @@ -106,7 +106,7 @@ subroutine set_state(scm_input, scm_reference, scm_state) end do end if end do - + !> - For each column, interpolate the v-wind to the model grid. do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres, scm_input%input_v, scm_state%pres_l(i,:), & @@ -118,7 +118,7 @@ subroutine set_state(scm_input, scm_reference, scm_state) end do end if end do - + !> - For each column, interpolate the ozone to the model grid. if(scm_state%ozone_index > 0) then do i=1, scm_state%n_cols @@ -132,7 +132,7 @@ subroutine set_state(scm_input, scm_reference, scm_state) end if end do end if - + !> - For each column, interpolate the cloud liquid water to the model grid. do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres, scm_input%input_ql, scm_state%pres_l(i,:), & @@ -142,7 +142,7 @@ subroutine set_state(scm_input, scm_reference, scm_state) scm_state%state_tracer(i,last_index_init+1:,scm_state%cloud_water_index,1) = 0.0 end if end do - + !> - For each column, interpolate the cloud ice water to the model grid. do i=1, scm_state%n_cols call interpolate_to_grid_centers(scm_input%input_nlev, scm_input%input_pres, scm_input%input_qi, scm_state%pres_l(i,:), & @@ -160,14 +160,14 @@ subroutine set_state(scm_input, scm_reference, scm_state) scm_state%state_T(i,:,1) = scm_input%input_temp(:) scm_state%state_tracer(i,:,scm_state%water_vapor_index,1)=scm_input%input_qt scm_state%state_tracer(i,:,scm_state%ozone_index,1)=scm_input%input_ozone - scm_state%area(i) = scm_input%input_area - + scm_state%area(i) = scm_input%input_area + if (scm_input%input_pres_i(1).GT. 0.0) then ! pressure are read in, overwrite values scm_state%pres_i(i,:)=scm_input%input_pres_i scm_state%pres_l(i,:)=scm_input%input_pres_l endif enddo - + endif !> @} end subroutine set_state @@ -291,8 +291,8 @@ subroutine GFS_suite_setup (Model, Statein, Stateout, Sfcprop, GFS_radtend_type, GFS_diag_type use CCPP_typedefs, only: GFS_interstitial_type use physcons, only: pi => con_pi - - + + !use cldwat2m_micro, only: ini_micro !use aer_cloud, only: aer_cloud_init !use module_ras, only: ras_init @@ -312,13 +312,13 @@ subroutine GFS_suite_setup (Model, Statein, Stateout, Sfcprop, type(GFS_init_type), intent(in) :: Init_parm integer, intent(in) :: ntasks, nthreads, n_cols - + real(kind=dp), dimension(n_cols), intent(in) :: lon, lat, area - + real(kind=dp), parameter :: rad2deg = 180.0_dp/pi - + integer :: i - + !--- set control properties (including namelist read) call Model%init (Init_parm%nlunit, Init_parm%fn_nml, & Init_parm%me, Init_parm%master, & @@ -336,7 +336,7 @@ subroutine GFS_suite_setup (Model, Statein, Stateout, Sfcprop, Init_parm%fcst_mpi_comm, ntasks, nthreads) !--- initialize DDTs - + call Statein%create(n_cols, Model) call Stateout%create(n_cols, Model) call Sfcprop%create(n_cols, Model) @@ -349,10 +349,10 @@ subroutine GFS_suite_setup (Model, Statein, Stateout, Sfcprop, call Diag%create(n_cols, Model) !--- internal representation of interstitials for CCPP physics call Interstitial%create(n_cols, Model) - + !--- populate the grid components !call GFS_grid_populate (Grid(i), Init_parm%xlon, Init_parm%xlat, Init_parm%area) - + do i=1, n_cols Grid%xlon(i) = lon(i) Grid%xlat(i) = lat(i) @@ -363,7 +363,7 @@ subroutine GFS_suite_setup (Model, Statein, Stateout, Sfcprop, Grid%area(i) = area(i) Grid%dx(i) = sqrt(area(i)) end do - + !--- initialize Morrison-Gettleman microphysics !if (Model%ncld == 2) then @@ -377,17 +377,17 @@ subroutine GFS_suite_setup (Model, Statein, Stateout, Sfcprop, !--- lsidea initialization if (Model%lsidea) then print *,' LSIDEA is active but needs to be reworked for FV3 - shutting down' - stop + error stop !--- NEED TO get the logic from the old phys/gloopb.f initialization area endif - + if(Model%do_ca)then print *,'Cellular automata cannot be used when CCPP is turned on until' print *,'the stochastic physics pattern generation code has been pulled' print *,'out of the FV3 repository and updated with the CCPP version.' - stop + error stop endif - + !--- sncovr may not exist in ICs from chgres. !--- FV3GFS handles this as part of the IC ingest !--- this not is placed here to alert users to the need to study diff --git a/scm/src/scm_time_integration.F90 b/scm/src/scm_time_integration.F90 index 4e5743c14..895d70435 100644 --- a/scm/src/scm_time_integration.F90 +++ b/scm/src/scm_time_integration.F90 @@ -11,7 +11,7 @@ module scm_time_integration only: ccpp_physics_timestep_init, & ccpp_physics_run, & ccpp_physics_timestep_finalize - + implicit none @@ -78,10 +78,9 @@ subroutine do_time_step(scm_state, physics, cdata, in_spinup) if (scm_state%input_type == 0) then call apply_forcing_leapfrog(scm_state) else - write(*,*) 'The application of forcing terms from the DEPHY file format has not been implemented for the leapfrog time scheme.' - stop + error stop 'The application of forcing terms from the DEPHY file format has not been implemented for the leapfrog time scheme.' end if - + case default if (scm_state%input_type == 0) then call apply_forcing_forward_Euler(scm_state, in_spinup) @@ -108,19 +107,19 @@ subroutine do_time_step(scm_state, physics, cdata, in_spinup) physics%Diag%dtend(i,:,idtend) = physics%Diag%dtend(i,:,idtend) & + (physics%Statein%ugrs(i,:) - physics%Stateout%gu0(i,:)) endif - + idtend = physics%Model%dtidx(physics%Model%index_of_y_wind,physics%Model%index_of_process_non_physics) if(idtend>=1) then physics%Diag%dtend(i,:,idtend) = physics%Diag%dtend(i,:,idtend) & + (physics%Statein%vgrs(i,:) - physics%Stateout%gv0(i,:)) endif - + idtend = physics%Model%dtidx(physics%Model%index_of_temperature,physics%Model%index_of_process_non_physics) if(idtend>=1) then physics%Diag%dtend(i,:,idtend) = physics%Diag%dtend(i,:,idtend) & + (physics%Statein%tgrs(i,:) - physics%Stateout%gt0(i,:)) endif - + if (physics%Model%qdiag3d) then do itrac=1,physics%Model%ntrac idtend = physics%Model%dtidx(itrac+100,physics%Model%index_of_process_non_physics) @@ -132,13 +131,13 @@ subroutine do_time_step(scm_state, physics, cdata, in_spinup) endif endif enddo - + call ccpp_physics_timestep_init(cdata, suite_name=trim(adjustl(scm_state%physics_suite_name)), ierr=ierr) if (ierr/=0) then write(*,'(a,i0,a)') 'An error occurred in ccpp_physics_timestep_init: ' // trim(cdata%errmsg) // '. Exiting...' - stop + error stop trim(cdata%errmsg) end if - + !--- determine if radiation diagnostics buckets need to be cleared if (nint(physics%Model%fhzero*3600) >= nint(max(physics%Model%fhswr,physics%Model%fhlwr))) then if (mod(physics%Model%kdt,physics%Model%nszero) == 1 .or. physics%Model%nszero == 1) then @@ -150,22 +149,22 @@ subroutine do_time_step(scm_state, physics, cdata, in_spinup) call physics%Diag%rad_zero (physics%Model) endif endif - + !--- determine if physics diagnostics buckets need to be cleared if (mod(physics%Model%kdt,physics%Model%nszero) == 1 .or. physics%Model%nszero == 1) then call physics%Diag%phys_zero (physics%Model) endif - + call ccpp_physics_run(cdata, suite_name=trim(adjustl(scm_state%physics_suite_name)), ierr=ierr) if (ierr/=0) then write(*,'(a,i0,a)') 'An error occurred in ccpp_physics_run: ' // trim(cdata%errmsg) // '. Exiting...' - stop + error stop trim(cdata%errmsg) end if call ccpp_physics_timestep_finalize(cdata, suite_name=trim(adjustl(scm_state%physics_suite_name)), ierr=ierr) if (ierr/=0) then write(*,'(a,i0,a)') 'An error occurred in ccpp_physics_timestep_finalize: ' // trim(cdata%errmsg) // '. Exiting...' - stop + error stop trim(cdata%errmsg) end if !if no physics call, need to transfer state_variables(:,:,1) to state_variables (:,:,2) diff --git a/scm/src/scm_utils.F90 b/scm/src/scm_utils.F90 index ce2183afd..dee65c4d8 100644 --- a/scm/src/scm_utils.F90 +++ b/scm/src/scm_utils.F90 @@ -543,7 +543,7 @@ subroutine check(status, var_name) if(status /= NF90_NOERR) then print *, trim(nf90_strerror(status)),' ',var_name - stop "stopped" + error stop "stopped" end if end subroutine check @@ -560,7 +560,7 @@ subroutine NetCDF_read_att_int(ncid, var_id, att_name, req, att_data) ierr = NF90_INQUIRE_ATTRIBUTE(ncid, var_id, att_name) if (ierr /= NF90_NOERR) then write(*,*) 'There was an error reading the required '//adjustl(trim(att_name))//' attribute. Stopping...' - stop + error stop "There was an error reading the required attribute" else call check(NF90_GET_ATT(ncid, var_id, att_name, att_data),att_name) end if @@ -588,7 +588,7 @@ subroutine NetCDF_read_att_sp(ncid, var_id, att_name, req, att_data) ierr = NF90_INQUIRE_ATTRIBUTE(ncid, var_id, att_name) if (ierr /= NF90_NOERR) then write(*,*) 'There was an error reading the required '//adjustl(trim(att_name))//' attribute. Stopping...' - stop + error stop "There was an error reading the required attribute" else call check(NF90_GET_ATT(ncid, var_id, att_name, att_data),att_name) end if @@ -616,7 +616,7 @@ subroutine NetCDF_read_att_char(ncid, var_id, att_name, req, att_data) ierr = NF90_INQUIRE_ATTRIBUTE(ncid, var_id, att_name) if (ierr /= NF90_NOERR) then write(*,*) 'There was an error reading the required '//adjustl(trim(att_name))//' attribute. Stopping...' - stop + error stop "There was an error reading the required attribute" else call check(NF90_GET_ATT(ncid, var_id, att_name, att_data),att_name) end if @@ -645,7 +645,7 @@ subroutine NetCDF_read_att_char_or_int(ncid, var_id, att_name, req, att_data, ch ierr = NF90_INQUIRE_ATTRIBUTE(ncid, var_id, att_name, xtype = type) if (ierr /= NF90_NOERR) then write(*,*) 'There was an error reading the required '//adjustl(trim(att_name))//' attribute. Stopping...' - stop + error stop "There was an error reading the required attribute" else if (type == NF90_CHAR) then call check(NF90_GET_ATT(ncid, var_id, att_name, char_att_data),att_name) @@ -691,7 +691,7 @@ subroutine NetCDF_conditionally_read_var_1d_sp(var_ctl, var_att, var_name, filen call NetCDF_read_var(ncid, var_name, .False., var_data) if (maxval(var_data) < missing_value_eps) then write(*,*) 'The global attribute '//var_att//' in '//filename//' indicates that the variable '//var_name//' should be present, but it is missing. Stopping ...' - stop + error stop "Missing variable" end if else var_data = missing_value @@ -710,7 +710,7 @@ subroutine NetCDF_conditionally_read_var_2d_sp(var_ctl, var_att, var_name, filen call NetCDF_read_var(ncid, var_name, .False., var_data) if (maxval(var_data) < missing_value_eps) then write(*,*) 'The global attribute '//var_att//' in '//filename//' indicates that the variable '//var_name//' should be present, but it is missing. Stopping ...' - stop + error stop "Missing variable" end if else var_data = missing_value @@ -729,7 +729,7 @@ subroutine NetCDF_conditionally_read_var_3d_sp(var_ctl, var_att, var_name, filen call NetCDF_read_var(ncid, var_name, .False., var_data) if (maxval(var_data) < missing_value_eps) then write(*,*) 'The global attribute '//var_att//' in '//filename//' indicates that the variable '//var_name//' should be present, but it is missing. Stopping ...' - stop + error stop "Missing variable" end if else var_data = missing_value @@ -748,7 +748,7 @@ subroutine NetCDF_conditionally_read_var_4d_sp(var_ctl, var_att, var_name, filen call NetCDF_read_var(ncid, var_name, .False., var_data) if (maxval(var_data) < missing_value_eps) then write(*,*) 'The global attribute '//var_att//' in '//filename//' indicates that the variable '//var_name//' should be present, but it is missing. Stopping ...' - stop + error stop "Missing variable" end if else var_data = missing_value @@ -789,7 +789,7 @@ subroutine NetCDF_def_var(ncid, var_name, var_type, desc, unit, varid, dims) call CHECK(NF90_PUT_ATT(NCID=ncid,VARID=varid,NAME="_FillValue",VALUES=missing_value_int),var_name) else write(0,'(a,i0,a)') "The variable '" // var_name // "' is defined as a type other than NF90_FLOAT or NF90_INT. Stopping..." - STOP + error stop "Variable defined as a type other than NF90_FLOAT or NF90_INT." end if end subroutine NetCDF_def_var @@ -933,7 +933,7 @@ subroutine conditionally_set_int_var_0d(input, set_var, input_name, req, missing else if (req) then write(0,'(a,i0,a)') "The variable '" // input_name // "' in the case data file had missing data, but it is required for the given physics configuration. Stopping..." - STOP + error stop "Variable in the case data file had missing data, but it is required for the given physics configuration." end if end if @@ -953,7 +953,7 @@ subroutine conditionally_set_int_var_1d(input, set_var, input_name, req, missing else if (req) then write(0,'(a,i0,a)') "The variable '" // input_name // "' in the case data file had missing data, but it is required for the given physics configuration. Stopping..." - STOP + error stop "The variable in the case data file had missing data, but it is required for the given physics configuration" end if end if @@ -973,7 +973,7 @@ subroutine conditionally_set_real_dp_var_0d(input, set_var, input_name, req, mis else if (req) then write(0,'(a,i0,a)') "The variable '" // input_name // "' in the case data file had missing data, but it is required for the given physics configuration. Stopping..." - STOP + error stop "The variable in the case data file had missing data, but it is required for the given physics configuration" end if end if @@ -993,7 +993,7 @@ subroutine conditionally_set_real_dp_var_1d(input, set_var, input_name, req, mis else if (req) then write(0,'(a,i0,a)') "The variable '" // input_name // "' in the case data file had missing data, but it is required for the given physics configuration. Stopping..." - STOP + error stop "The variable in the case data file had missing data, but it is required for the given physics configuration" end if end if diff --git a/scm/src/scm_vgrid.F90 b/scm/src/scm_vgrid.F90 index 12ac24fc0..e492824fe 100644 --- a/scm/src/scm_vgrid.F90 +++ b/scm/src/scm_vgrid.F90 @@ -36,16 +36,16 @@ subroutine get_FV3_vgrid(scm_input, scm_state) ! local real(kind=dp) :: pres_sfc_inv, p_ref, ak_tmp, bk_tmp - + real:: p0=1000.E2 real:: pc=200.E2 - + real :: pint = 100.E2 real :: stretch_fac = 1.03 integer :: auto_routine = 0 - + integer k, last_index, mid_index, ierr, dummy, n_levels_file - + character(len=80) :: line character(len=16) :: file_format integer :: nx, ny @@ -55,13 +55,13 @@ subroutine get_FV3_vgrid(scm_input, scm_state) ! single and double precision runs real(kind_scm_sp), parameter :: zero_sp = 0.0 real(kind_scm_sp), allocatable :: pres_l_row_sp(:) - + #include "fv_eta.h" - + km = scm_state%n_levels ptop = 1. ! Definition: press(i,j,k) = ak(k) + bk(k) * ps(i,j) - + if (trim(scm_state%npz_type) == 'superC' .or. trim(scm_state%npz_type) == 'superK') then auto_routine = 1 select case (km) @@ -99,33 +99,33 @@ subroutine get_FV3_vgrid(scm_input, scm_state) auto_routine = 2 end select else if (trim(scm_state%npz_type) == 'input') then - + !> - Open the appropriate file. open(unit=1, file=scm_state%vert_coord_file, status='old', action='read', iostat=ierr) if(ierr /= 0) then write(*,*) 'There was an error opening the file ', scm_state%vert_coord_file, ' in the run directory. & Error code = ',ierr - stop + error stop endif - + !> The file being read in must have the following format: !! include a single line description: number of coefficients, number of layers !! ak/bk pairs, with each pair occupying a single line !! the pairs must be ordered from surface to TOA !! the pairs define the interfaces of the grid to create levels-1 layer - - !> - The first line contains the number of coefficients and number of levels + + !> - The first line contains the number of coefficients and number of levels read(1,*) dummy, n_levels_file if (n_levels_file /= scm_state%n_levels) then write(*,*) 'There is a mismatch in the number of levels expected and the number of coefficients supplied in the file ',scm_state%vert_coord_file - stop + error stop end if - !> - Read in the coefficient data. + !> - Read in the coefficient data. do k=1, km+1 read(1,*) scm_state%a_k(k), scm_state%b_k(k) end do close(1) - + ! flip scm_state%a_k, scm_state%b_k in vertical (a_k and b_k are expected to be TOA-to-surface at this point) mid_index = (km+1)/2 last_index = km+1 @@ -138,7 +138,7 @@ subroutine get_FV3_vgrid(scm_input, scm_state) scm_state%b_k(last_index) = bk_tmp last_index = last_index - 1 end do - + else select case (km) case (5,10) ! does this work???? @@ -217,7 +217,7 @@ subroutine get_FV3_vgrid(scm_input, scm_state) pint = 100.E2 stretch_fac = 1.035 auto_routine = 1 - case (47) + case (47) if (trim(scm_state%npz_type) == 'lowtop') then ptop = 100. stretch_fac = 1.035 @@ -504,7 +504,7 @@ subroutine get_FV3_vgrid(scm_input, scm_state) endif end select endif ! superC/superK - + select case (auto_routine) case (1) call var_hi(km, scm_state%a_k, scm_state%b_k, ptop, ks, pint, stretch_fac) @@ -519,7 +519,7 @@ subroutine get_FV3_vgrid(scm_input, scm_state) case (6) call var_gfs(km, scm_state%a_k, scm_state%b_k, ptop, ks, pint, stretch_fac) end select - + call check_eta_levels (scm_state%a_k, scm_state%b_k) if (verbose) then @@ -528,7 +528,7 @@ subroutine get_FV3_vgrid(scm_input, scm_state) write(*,'(I4, F13.5, F13.5, F11.2)') k, scm_state%a_k(k), scm_state%b_k(k), 1000.E2*scm_state%b_k(k) + scm_state%a_k(k) enddo endif - + !> - Calculate interface pressures, sigma, and exner function. ! flip scm_state%a_k, scm_state%b_k in vertical @@ -905,7 +905,7 @@ subroutine var_hi2(km, ak, bk, ptop, ks, pint, s_rate) end subroutine var_hi2 - + subroutine var_les(km, ak, bk, ptop, ks, pint, s_rate) implicit none integer, intent(in):: km @@ -1065,7 +1065,7 @@ subroutine var_les(km, ak, bk, ptop, ks, pint, s_rate) 800 format(1x,5(1x,f9.4)) end subroutine var_les - + subroutine mount_waves(km, ak, bk, ptop, ks, pint) integer, intent(in):: km real, intent(out):: ak(km+1), bk(km+1) @@ -1586,7 +1586,7 @@ subroutine check_eta_levels(ak, bk) enddo endif write(*,*) 'FV3 check_eta_levels: ak/bk pairs do not provide a monotonic vertical coordinate' - stop + error stop endif end subroutine check_eta_levels From ae787929b148f255fe25e83935c19b96c5246020 Mon Sep 17 00:00:00 2001 From: Soren Rasmussen Date: Thu, 25 Jul 2024 11:48:26 -0600 Subject: [PATCH 3/9] Docker will checkout PR if being run during PR CI --- .github/workflows/ci_test_docker.yml | 4 +++- docker/Dockerfile | 23 +++++++++++++++++------ scm/src/run_scm.py | 12 ++++++------ 3 files changed, 26 insertions(+), 13 deletions(-) diff --git a/.github/workflows/ci_test_docker.yml b/.github/workflows/ci_test_docker.yml index 51991a10d..730a5381c 100644 --- a/.github/workflows/ci_test_docker.yml +++ b/.github/workflows/ci_test_docker.yml @@ -4,6 +4,7 @@ on: [pull_request,workflow_dispatch] env: TEST_TAG: dtcenter/ccpp-scm:test + PR_NUMBER: ${{ github.event.number }} jobs: docker: @@ -22,8 +23,9 @@ jobs: file: docker/Dockerfile load: true tags: ${{ env.TEST_TAG }} + build-args: PR_NUMBER=${{ github.event.number }} - name: Test run: | mkdir $HOME/output chmod a+rw $HOME/output - docker run --rm -v $HOME/output:/home ${{ env.TEST_TAG }} ./run_scm.py -f ../../test/rt_test_cases.py --runtime_mult 0.1 -d -a + docker run --rm -v $HOME/output:/home ${{ env.TEST_TAG }} ./run_scm.py -f ../../test/rt_test_cases.py --runtime_mult 0.1 -d --stop_on_error diff --git a/docker/Dockerfile b/docker/Dockerfile index 85ad4d856..f269fff5f 100644 --- a/docker/Dockerfile +++ b/docker/Dockerfile @@ -1,8 +1,10 @@ FROM debian:12 LABEL maintainer="Michael Kavulich " -# Set up base OS environment +# arguments that can be passed in +ARG PR_NUMBER +# Set up base OS environment RUN apt-get -y update # Get "essential" tools and libraries @@ -70,11 +72,20 @@ ENV sp_ROOT=/comsoftware/nceplibs ENV w3emc_ROOT=/comsoftware/nceplibs # Obtain CCPP SCM source code, build code, and download static data -RUN cd /comsoftware \ - && git clone --recursive -b main https://github.com/NCAR/ccpp-scm \ - && mkdir /comsoftware/ccpp-scm/scm/bin - -RUN cd /comsoftware/ccpp-scm/scm/bin \ +RUN if [ -z "$PR_NUMBER" ]; then \ + cd /comsoftware \ + && git clone --recursive -b main https://github.com/NCAR/ccpp-scm; \ + else \ + cd /comsoftware \ + && git clone https://github.com/NCAR/ccpp-scm \ + && cd ccpp-scm \ + && git fetch origin pull/${PR_NUMBER}/head:test_pr \ + && git checkout test_pr \ + && git submodule update --init --recursive; \ + fi + +RUN mkdir /comsoftware/ccpp-scm/scm/bin \ + && cd /comsoftware/ccpp-scm/scm/bin \ && cmake ../src \ && make -j diff --git a/scm/src/run_scm.py b/scm/src/run_scm.py index 35222f168..40c658843 100755 --- a/scm/src/run_scm.py +++ b/scm/src/run_scm.py @@ -34,6 +34,7 @@ # Default command string to run MPI apps (number of processes should be 1 since SCM is not set up to use more than 1 yet) DEFAULT_MPI_COMMAND = 'mpirun -np 1' + # Copy executable to run directory if true (otherwise it will be linked) COPY_EXECUTABLE = False @@ -120,7 +121,7 @@ parser.add_argument('--n_itt_out', help='period of instantaneous output (number of timesteps)', required=False, type=int) parser.add_argument('--n_itt_diag', help='period of diagnostic output (number of timesteps)', required=False, type=int) parser.add_argument('-dt', '--timestep', help='timestep (s)', required=False, type=float) -parser.add_argument('-a', '--actions', help='if running from Github Actions this will select correct configuration', required=False, action='store_true') +parser.add_argument('--stop_on_error', help='when running multiple SCM runs, stop on first error', required=False, action='store_true') parser.add_argument('-v', '--verbose', help='set logging level to debug and write log to file', action='count', default=0) parser.add_argument('-f', '--file', help='name of file where SCM runs are defined') parser.add_argument('--mpi_command', help='command used to invoke the executable via MPI (including options)', required=False) @@ -191,14 +192,14 @@ def parse_arguments(): bin_dir = args.bin_dir timestep = args.timestep mpi_command = args.mpi_command - github_actions = args.actions + stop_on_error = args.stop_on_error if not sdf: sdf = DEFAULT_SUITE return (file, case, sdf, namelist, tracers, gdb, runtime, runtime_mult, docker, \ verbose, levels, npz_type, vert_coord_file, case_data_dir, n_itt_out, \ - n_itt_diag, run_dir, bin_dir, timestep, mpi_command, github_actions) + n_itt_diag, run_dir, bin_dir, timestep, mpi_command, stop_on_error) def find_gdb(): """Detect gdb, abort if not found""" @@ -767,7 +768,7 @@ def copy_outdir(exp_dir): def main(): (file, case, sdf, namelist, tracers, use_gdb, runtime, runtime_mult, docker, \ verbose, levels, npz_type, vert_coord_file, case_data_dir, n_itt_out, \ - n_itt_diag, run_dir, bin_dir, timestep, mpi_command, github_actions) = parse_arguments() + n_itt_diag, run_dir, bin_dir, timestep, mpi_command, stop_on_error) = parse_arguments() setup_logging(verbose) @@ -884,8 +885,7 @@ def main(): l_ignore_error = MULTIRUN_IGNORE_ERROR else: l_ignore_error = False - # need to correctly fail if running Github Actions - if github_actions: + if stop_on_error: l_ignore_error = False (status, time_elapsed) = launch_executable(use_gdb, gdb, mpi_command, ignore_error = l_ignore_error) From f935b13a09d19f784520232ee7a364c00ad9e72a Mon Sep 17 00:00:00 2001 From: Soren Rasmussen Date: Fri, 26 Jul 2024 12:55:26 -0600 Subject: [PATCH 4/9] Print summary of failed cases and exit it with an error if a failed cases exists Continue docker run so report of all failed tests gets produced. --- .github/workflows/ci_test_docker.yml | 2 +- scm/src/run_scm.py | 54 +++++++++++++++++++++++----- 2 files changed, 47 insertions(+), 9 deletions(-) diff --git a/.github/workflows/ci_test_docker.yml b/.github/workflows/ci_test_docker.yml index 730a5381c..b1dd5e747 100644 --- a/.github/workflows/ci_test_docker.yml +++ b/.github/workflows/ci_test_docker.yml @@ -28,4 +28,4 @@ jobs: run: | mkdir $HOME/output chmod a+rw $HOME/output - docker run --rm -v $HOME/output:/home ${{ env.TEST_TAG }} ./run_scm.py -f ../../test/rt_test_cases.py --runtime_mult 0.1 -d --stop_on_error + docker run --rm -v $HOME/output:/home ${{ env.TEST_TAG }} ./run_scm.py -f ../../test/rt_test_cases.py --runtime_mult 0.1 -d diff --git a/scm/src/run_scm.py b/scm/src/run_scm.py index 40c658843..c0febaa48 100755 --- a/scm/src/run_scm.py +++ b/scm/src/run_scm.py @@ -3,6 +3,7 @@ import argparse import f90nml import logging +import numpy as np import os import re import shutil @@ -156,17 +157,18 @@ def execute(cmd, ignore_error = False): stderr = subprocess.PIPE, shell = True) (stdout, stderr) = p.communicate() status = p.returncode + message = 'Execution of "{0}" returned with exit code {1}\n'.format(cmd, status) + message += ' stdout: "{0}"\n'.format(stdout.decode(encoding='ascii', errors='ignore').rstrip('\n')) + message += ' stderr: "{0}"'.format(stderr.decode(encoding='ascii', errors='ignore').rstrip('\n')) + if status == 0: - message = 'Execution of "{0}" returned with exit code {1}\n'.format(cmd, status) - message += ' stdout: "{0}"\n'.format(stdout.decode(encoding='ascii', errors='ignore').rstrip('\n')) - message += ' stderr: "{0}"'.format(stderr.decode(encoding='ascii', errors='ignore').rstrip('\n')) logging.debug(message) elif not ignore_error: - message = 'Execution of command "{0}" failed, exit code {1}\n'.format(cmd, status) - message += ' stdout: "{0}"\n'.format(stdout.decode(encoding='ascii', errors='ignore').rstrip('\n')) - message += ' stderr: "{0}"'.format(stderr.decode(encoding='ascii', errors='ignore').rstrip('\n')) logging.critical(message) raise Exception('Execution of command "{0}" failed, exit code {1}\n'.format(cmd, status)) + else: + print("SHOULD PRINT ERROR MESSAGE: status and ignore error ==", status, ignore_error) + logging.error(message) return (status, stdout.decode(encoding='ascii', errors='ignore').rstrip('\n'), stderr.decode(encoding='ascii', errors='ignore').rstrip('\n')) def parse_arguments(): @@ -765,6 +767,30 @@ def copy_outdir(exp_dir): shutil.rmtree(home_output_dir) shutil.copytree(exp_dir, home_output_dir) +def print_error_report(error_logs, total_count): + case_l = len(max(error_logs[:,0], key=len)) + suite_l = len(max(error_logs[:,1], key=len)) + namelist_l = len(max(error_logs[:,2], key=len)) + status_l = len(max(error_logs[:,3], key=len)) + # error_log contains header, subtracting 1 from error + error_count = error_logs.shape[0] - 1 + passing_count = total_count - error_count + header_printed = False + column_width = (case_l + suite_l + namelist_l + status_l + 13) + + # print formatted asummary + print("Failure Summary:") + print("-" * column_width) + for error_log in error_logs: + case_s, suite, namelist, status = error_log + print(f"| {case_s:<{case_l}} | {suite:<{suite_l}} | {namelist:<{namelist_l}} | {status:<{status_l}} |") + if not header_printed: + print("-" * column_width) + header_printed = True + print("-" * column_width) + print(f"[{error_count}/{total_count}] failed cases, [{passing_count}/{total_count}] passing cases") + + def main(): (file, case, sdf, namelist, tracers, use_gdb, runtime, runtime_mult, docker, \ verbose, levels, npz_type, vert_coord_file, case_data_dir, n_itt_out, \ @@ -822,6 +848,8 @@ def main(): if (tracers != None): run_list[0]["tracers"] = tracers # Loop through all input "run dictionaires" + error_logs = [["Failed Case", "Suite", "Namelist", "Status"]] + failed_case = False irun = 0 for run in run_list: @@ -894,13 +922,23 @@ def main(): logging.info('Process "(case={0}, suite={1}, namelist={2}" completed successfully'. \ format(run["case"], run["suite"], active_suite.namelist)) else: - logging.warning('Process "(case={0}, suite={1}, namelist={2}" exited with code {3}'. \ - format( run["case"], run["suite"], active_suite.namelist, status)) + failed_case = True + error_str = 'Process "(case={0}, suite={1}, namelist={2}" exited with code {3}'. \ + format( run["case"], run["suite"], active_suite.namelist, status) + logging.warning(error_str) + error_logs = np.append(error_logs, + [[run["case"], run["suite"], active_suite.namelist, status]], + axis=0) # if time_elapsed: logging.info(' Elapsed time: {0}s'.format(time_elapsed)) if docker: copy_outdir(exp_dir) + if (failed_case): + print_error_report(error_logs, len(run_list)) + sys.exit(1) + + if __name__ == '__main__': main() From b394f933c222d9b26f798b4e06dbe1b816d46e6d Mon Sep 17 00:00:00 2001 From: Soren Rasmussen Date: Fri, 26 Jul 2024 15:36:44 -0600 Subject: [PATCH 5/9] Printing formatted string for each experiment to making searching CI output easy --- scm/src/run_scm.py | 56 ++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 47 insertions(+), 9 deletions(-) diff --git a/scm/src/run_scm.py b/scm/src/run_scm.py index c0febaa48..87d2020c5 100755 --- a/scm/src/run_scm.py +++ b/scm/src/run_scm.py @@ -167,7 +167,6 @@ def execute(cmd, ignore_error = False): logging.critical(message) raise Exception('Execution of command "{0}" failed, exit code {1}\n'.format(cmd, status)) else: - print("SHOULD PRINT ERROR MESSAGE: status and ignore error ==", status, ignore_error) logging.error(message) return (status, stdout.decode(encoding='ascii', errors='ignore').rstrip('\n'), stderr.decode(encoding='ascii', errors='ignore').rstrip('\n')) @@ -767,11 +766,19 @@ def copy_outdir(exp_dir): shutil.rmtree(home_output_dir) shutil.copytree(exp_dir, home_output_dir) -def print_error_report(error_logs, total_count): - case_l = len(max(error_logs[:,0], key=len)) - suite_l = len(max(error_logs[:,1], key=len)) - namelist_l = len(max(error_logs[:,2], key=len)) - status_l = len(max(error_logs[:,3], key=len)) + +def print_report_line(case_s, suite, namelist, max_str_lens): + case_l = max_str_lens.case + suite_l = max_str_lens.suite + namelist_l = max_str_lens.namelist + logging.info(f"| {case_s:<{case_l}} | {suite:<{suite_l}} | {namelist:<{namelist_l}} |") + + +def print_error_report(error_logs, total_count, max_str_lens): + case_l = max_str_lens.case + suite_l = max_str_lens.suite + namelist_l = max_str_lens.namelist + status_l = max_str_lens.status # error_log contains header, subtracting 1 from error error_count = error_logs.shape[0] - 1 passing_count = total_count - error_count @@ -791,6 +798,32 @@ def print_error_report(error_logs, total_count): print(f"[{error_count}/{total_count}] failed cases, [{passing_count}/{total_count}] passing cases") +class MaxStrLengths: + def __init__(self, max_case_len, max_suite_len, + max_namelist_len, max_status_len): + self.case = max_case_len + self.suite = max_suite_len + self.namelist = max_namelist_len + self.status = max_status_len + + +def find_max_str_lengths(run_list): + max_case_len = 0 + max_suite_len = 0 + max_status_len = len('status') + + # Loop through the list of dictionaries to find the longest lengths + for item in run_list: + max_case_len = max(max_case_len, len(item['case'])) + max_suite_len = max(max_suite_len, len(item['suite'])) + + # add 6, e.g. diff between len('SCM_HRRR_gf') and len('input_HRRR_gf.nml') + max_namelist_len = max_suite_len + 6 + max_str_lens = MaxStrLengths(max_case_len, max_suite_len, + max_namelist_len, max_status_len) + return max_str_lens + + def main(): (file, case, sdf, namelist, tracers, use_gdb, runtime, runtime_mult, docker, \ verbose, levels, npz_type, vert_coord_file, case_data_dir, n_itt_out, \ @@ -847,15 +880,19 @@ def main(): if (namelist != None): run_list[0]["namelist"] = namelist if (tracers != None): run_list[0]["tracers"] = tracers - # Loop through all input "run dictionaires" + + # setup variables error_logs = [["Failed Case", "Suite", "Namelist", "Status"]] + max_str_lens = find_max_str_lengths(run_list) failed_case = False irun = 0 + + # Loop through all input "run dictionaires" for run in run_list: # # Is this a "supported" SCM configuration? - # (e.g Do we have defualt namelist and tracer files for this suite?) + # (e.g Do we have default namelist and tracer files for this suite?) # If supported, copy default configuration, modify below if necessary. # active_suite = None @@ -902,6 +939,7 @@ def main(): # # Run the SCM case # + print_report_line(run["case"], run["suite"], active_suite.namelist, max_str_lens) logging.info('Executing process {0} of {1}: case={2}, suite={3}, namelist={4}'.format( irun, len(run_list), run["case"], run["suite"], active_suite.namelist)) # @@ -936,7 +974,7 @@ def main(): copy_outdir(exp_dir) if (failed_case): - print_error_report(error_logs, len(run_list)) + print_error_report(error_logs, len(run_list), max_str_lens) sys.exit(1) From 1a54e8de500a2616ba8e3955fafbd24cd70f2a3d Mon Sep 17 00:00:00 2001 From: Soren Rasmussen Date: Fri, 2 Aug 2024 11:36:06 -0600 Subject: [PATCH 6/9] CCPP Physics branch: bugfix/scm-rts --- ccpp/physics | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ccpp/physics b/ccpp/physics index 62f7656eb..e2f4a6368 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 62f7656ebdaa39d989413e5d285ed6d6723eb4de +Subproject commit e2f4a63682c679660f0fd334b34e7c5c8c4042d4 From 4288c2a0f657a60e86c99390124da615b6ebe683 Mon Sep 17 00:00:00 2001 From: Soren Rasmussen Date: Fri, 9 Aug 2024 13:02:57 -0600 Subject: [PATCH 7/9] Print passing summary as well as the failing summary --- scm/src/run_scm.py | 32 +++++++++++++++++++++++--------- 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/scm/src/run_scm.py b/scm/src/run_scm.py index 87d2020c5..520e1fb93 100755 --- a/scm/src/run_scm.py +++ b/scm/src/run_scm.py @@ -774,28 +774,36 @@ def print_report_line(case_s, suite, namelist, max_str_lens): logging.info(f"| {case_s:<{case_l}} | {suite:<{suite_l}} | {namelist:<{namelist_l}} |") -def print_error_report(error_logs, total_count, max_str_lens): +def print_report(logs, total_count, max_str_lens, + passing=False, failing=False): case_l = max_str_lens.case suite_l = max_str_lens.suite namelist_l = max_str_lens.namelist status_l = max_str_lens.status # error_log contains header, subtracting 1 from error - error_count = error_logs.shape[0] - 1 + error_count = logs.shape[0] - 1 passing_count = total_count - error_count header_printed = False column_width = (case_l + suite_l + namelist_l + status_l + 13) - # print formatted asummary - print("Failure Summary:") + # print formatted summary + print("") + if (passing): + print("Passing Summary:") + if (failing): + print("Failure Summary:") print("-" * column_width) - for error_log in error_logs: - case_s, suite, namelist, status = error_log + for log in logs: + case_s, suite, namelist, status = log print(f"| {case_s:<{case_l}} | {suite:<{suite_l}} | {namelist:<{namelist_l}} | {status:<{status_l}} |") if not header_printed: print("-" * column_width) header_printed = True print("-" * column_width) - print(f"[{error_count}/{total_count}] failed cases, [{passing_count}/{total_count}] passing cases") + if (passing): + print(f"[{passing_count}/{total_count}] passing cases") + if (failing): + print(f"[{error_count}/{total_count}] failed cases") class MaxStrLengths: @@ -827,7 +835,8 @@ def find_max_str_lengths(run_list): def main(): (file, case, sdf, namelist, tracers, use_gdb, runtime, runtime_mult, docker, \ verbose, levels, npz_type, vert_coord_file, case_data_dir, n_itt_out, \ - n_itt_diag, run_dir, bin_dir, timestep, mpi_command, stop_on_error) = parse_arguments() + n_itt_diag, run_dir, bin_dir, timestep, mpi_command, stop_on_error \ + ) = parse_arguments() setup_logging(verbose) @@ -883,6 +892,7 @@ def main(): # setup variables error_logs = [["Failed Case", "Suite", "Namelist", "Status"]] + pass_logs = [["Passing Case", "Suite", "Namelist", "Status"]] max_str_lens = find_max_str_lengths(run_list) failed_case = False irun = 0 @@ -959,6 +969,9 @@ def main(): if status == 0: logging.info('Process "(case={0}, suite={1}, namelist={2}" completed successfully'. \ format(run["case"], run["suite"], active_suite.namelist)) + pass_logs = np.append(pass_logs, + [[run["case"], run["suite"], active_suite.namelist, status]], + axis=0) else: failed_case = True error_str = 'Process "(case={0}, suite={1}, namelist={2}" exited with code {3}'. \ @@ -973,8 +986,9 @@ def main(): if docker: copy_outdir(exp_dir) + print_report(pass_logs, len(run_list), max_str_lens, passing=True) if (failed_case): - print_error_report(error_logs, len(run_list), max_str_lens) + print_report(error_logs, len(run_list), max_str_lens, failing=True) sys.exit(1) From c2447bb6403bd7d49ec9322e7314addcac2a4ad9 Mon Sep 17 00:00:00 2001 From: Soren Rasmussen Date: Fri, 9 Aug 2024 13:06:33 -0600 Subject: [PATCH 8/9] Adding merge pull request #1081 --- ccpp/physics | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ccpp/physics b/ccpp/physics index e2f4a6368..7506fd110 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit e2f4a63682c679660f0fd334b34e7c5c8c4042d4 +Subproject commit 7506fd110dbeefca83ecc6ecd7f50c0b4989694a From 06ebb1761a5dbab4ccb6ae2658bc44c3cd0d80a2 Mon Sep 17 00:00:00 2001 From: Soren Rasmussen Date: Fri, 9 Aug 2024 13:50:50 -0600 Subject: [PATCH 9/9] fixing passing/failing counts in CI summary --- scm/src/run_scm.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/scm/src/run_scm.py b/scm/src/run_scm.py index 520e1fb93..93f09435f 100755 --- a/scm/src/run_scm.py +++ b/scm/src/run_scm.py @@ -780,9 +780,6 @@ def print_report(logs, total_count, max_str_lens, suite_l = max_str_lens.suite namelist_l = max_str_lens.namelist status_l = max_str_lens.status - # error_log contains header, subtracting 1 from error - error_count = logs.shape[0] - 1 - passing_count = total_count - error_count header_printed = False column_width = (case_l + suite_l + namelist_l + status_l + 13) @@ -799,10 +796,14 @@ def print_report(logs, total_count, max_str_lens, if not header_printed: print("-" * column_width) header_printed = True + print("-" * column_width) if (passing): + # error_log contains header, subtracting 1 from error + passing_count = logs.shape[0] - 1 print(f"[{passing_count}/{total_count}] passing cases") if (failing): + error_count = logs.shape[0] - 1 print(f"[{error_count}/{total_count}] failed cases")