diff --git a/.github/PULL_REQUEST_TEMPLATE/pull_request_template.md b/.github/pull_request_template.md similarity index 76% rename from .github/PULL_REQUEST_TEMPLATE/pull_request_template.md rename to .github/pull_request_template.md index 5968702e5..9f33913cd 100644 --- a/.github/PULL_REQUEST_TEMPLATE/pull_request_template.md +++ b/.github/pull_request_template.md @@ -1,22 +1,21 @@ -## Description of Changes: +## Description of Changes: One or more paragraphs describing the problem, solution, and required changes. -## Tests Conducted: -Explicitly state what tests were run on these changes, or if any are still pending (for README or other text-only changes, just put "None required". Make note of the compilers used, the platform/machine, and other relevant details as necessary. For more complicated changes, or those resulting in scientific changes, please be explicit! -**OR** Add any links to tests conducted. For example, "See ufs-community/ufs-weather-model/pull/" +## Tests Conducted: +Explicitly state what tests were run on these changes, or if any are still pending (for README or other text-only changes, just put "None required". Make note of the compilers used, the platform/machine, and other relevant details as necessary. For more complicated changes, or those resulting in scientific changes, please be explicit! +**OR** Add any links to tests conducted. For example, "See ufs-community/ufs-weather-model#" ## Dependencies: Add any links to parent PRs (e.g. SCM and/or UFS PRs) or submodules (e.g. rte-rrtmgp). For example: -- NCAR/ccpp-framework/pull/ -- NOAA-EMC/fv3atm/pull/ -- ufs-community/ufs-weather-model/pull/ +- NCAR/ccpp-framework# +- NOAA-EMC/fv3atm# +- ufs-community/ufs-weather-model/# ## Documentation: Does this PR add new capabilities that need to be documented or require modifications to the existing documentation? If so, brief supporting material can be provided here. Contact the CODEOWNERS if your PR requires extensive updates to the documentation. See https://github.com/NCAR/ccpp-doc for Technical Documentation or https://dtcenter.org/community-code/common-community-physics-package-ccpp/documentation for the latest Scientific Documentation. -## Issue (optional): -If this PR is resolving or referencing one or more issues, in this repository or elewhere, list them here. For example, "Fixes issue mentioned in #123" or "Related to bug in https://github.com/NCAR/other_repository/pull/63" +## Issue (optional): +If this PR is resolving or referencing one or more issues, in this repository or elewhere, list them here. For example, "Fixes issue mentioned in #123" or "Related to bug in NCAR/other_repository#123" -## Contributors (optional): +## Contributors (optional): If others have contributed to this work aside from the PR author, list them here - diff --git a/.github/workflows/basic_checks.yml b/.github/workflows/basic_checks.yml index 4e40790b5..7730b423d 100644 --- a/.github/workflows/basic_checks.yml +++ b/.github/workflows/basic_checks.yml @@ -4,6 +4,7 @@ on: [push, pull_request] jobs: build: + if: github.repository == 'NCAR/ccpp-physics' || github.repository == 'ufs-community/ccpp-physics' runs-on: macos-latest diff --git a/.github/workflows/ci_fv3_ccpp_prebuild.yml b/.github/workflows/ci_fv3_ccpp_prebuild.yml index a5c2f8092..5d87f7d50 100644 --- a/.github/workflows/ci_fv3_ccpp_prebuild.yml +++ b/.github/workflows/ci_fv3_ccpp_prebuild.yml @@ -4,9 +4,10 @@ on: [push, pull_request] jobs: ccpp-prebuild-FV3: + if: github.repository == 'NCAR/ccpp-physics' || github.repository == 'ufs-community/ccpp-physics' # The type of runner that the job will run on - runs-on: ubuntu-20.04 + runs-on: ubuntu-latest steps: - name: Checkout current ccpp-physics code @@ -38,11 +39,11 @@ jobs: git remote add remote_local $GIT_REMOTE_URL git fetch remote_local $GIT_REMOTE_BRANCH git checkout remote_local/$GIT_REMOTE_BRANCH - - - name: Set up Python 3.8.5 + + - name: Set up Python 3.10.13 uses: actions/setup-python@v3 with: - python-version: 3.8.5 + python-version: 3.10.13 - name: Add conda to system path run: | @@ -53,4 +54,4 @@ jobs: run: | cd /home/runner/work/ccpp-physics/ccpp-physics/fv3atm/ccpp/ mkdir -p /home/runner/work/ccpp-physics/ccpp-physics/fv3atm/bin/ccpp/physics/physics/ - ./framework/scripts/ccpp_prebuild.py --config config/ccpp_prebuild_config.py \ No newline at end of file + ./framework/scripts/ccpp_prebuild.py --config config/ccpp_prebuild_config.py diff --git a/.github/workflows/ci_scm_ccpp_prebuild.yml b/.github/workflows/ci_scm_ccpp_prebuild.yml index 7c9d2300a..d50dafcef 100644 --- a/.github/workflows/ci_scm_ccpp_prebuild.yml +++ b/.github/workflows/ci_scm_ccpp_prebuild.yml @@ -4,9 +4,10 @@ on: [push, pull_request] jobs: ccpp-prebuild-SCM: + if: github.repository == 'NCAR/ccpp-physics' || github.repository == 'ufs-community/ccpp-physics' # The type of runner that the job will run on - runs-on: ubuntu-20.04 + runs-on: ubuntu-latest steps: @@ -42,11 +43,11 @@ jobs: git remote add remote_local $GIT_REMOTE_URL git fetch remote_local $GIT_REMOTE_BRANCH git checkout remote_local/$GIT_REMOTE_BRANCH - - - name: Set up Python 3.8.5 + + - name: Set up Python 3.10.13 uses: actions/setup-python@v3 with: - python-version: 3.8.5 + python-version: 3.10.13 - name: Add conda to system path run: | @@ -58,4 +59,4 @@ jobs: cd /home/runner/work/ccpp-physics/ccpp-physics/ccpp-scm/ git status mkdir -p /home/runner/work/ccpp-physics/ccpp-physics/ccpp-scm/scm/bin/ccpp/physics/physics/ - ./ccpp/framework/scripts/ccpp_prebuild.py --config ccpp/config/ccpp_prebuild_config.py \ No newline at end of file + ./ccpp/framework/scripts/ccpp_prebuild.py --config ccpp/config/ccpp_prebuild_config.py diff --git a/CMakeLists.txt b/CMakeLists.txt index db7fae815..3d2a2971e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -196,7 +196,12 @@ target_include_directories(ccpp_physics PUBLIC target_link_libraries(ccpp_physics PRIVATE MPI::MPI_Fortran) target_link_libraries(ccpp_physics PUBLIC w3emc::w3emc_d sp::sp_d - NetCDF::NetCDF_Fortran) + NetCDF::NetCDF_Fortran + ) +#add FMS for FV3 only +if(FV3) + target_link_libraries(ccpp_physics PUBLIC fms) +endif() # Define where to install the library install(TARGETS ccpp_physics diff --git a/CODEOWNERS b/CODEOWNERS index 257be0187..953c77858 100644 --- a/CODEOWNERS +++ b/CODEOWNERS @@ -36,7 +36,6 @@ physics/GWD/ugwpv1_gsldrag_post.* @md physics/GWD/unified_ugwp* @mdtoyNOAA @grantfirl @rhaesung @Qingfu-Liu @dustinswales physics/MP/Ferrier_Aligo/module_MP_FER_HIRES.* @ericaligo-NOAA @grantfirl @rhaesung @Qingfu-Liu @dustinswales physics/MP/Ferrier_Aligo/mp_fer_hires.* @ericaligo-NOAA @grantfirl @rhaesung @Qingfu-Liu @dustinswales -physics/MP/GFDL/GFDL_parse_tracers.F90 @grantfirl @rhaesung @Qingfu-Liu @dustinswales physics/MP/GFDL/gfdl_cloud_microphys.* @RuiyuSun @grantfirl @rhaesung @Qingfu-Liu @dustinswales physics/MP/GFDL/module_gfdl_cloud_microphys.* @RuiyuSun @grantfirl @rhaesung @Qingfu-Liu @dustinswales physics/MP/GFDL/fv_sat_adj.* @RuiyuSun @grantfirl @rhaesung @Qingfu-Liu @dustinswales diff --git a/physics/CONV/C3/cu_c3_driver.F90 b/physics/CONV/C3/cu_c3_driver.F90 index 0ea8be075..c7e1c1f8c 100644 --- a/physics/CONV/C3/cu_c3_driver.F90 +++ b/physics/CONV/C3/cu_c3_driver.F90 @@ -657,7 +657,9 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,& enddo do i = its,itf if(mconv(i).lt.0.)mconv(i)=0. - if((dx(i)<6500.).and.do_mynnedmf.and.(maxMF(i).gt.0.))ierr(i)=555 + if(do_mynnedmf) then + if((dx(i)<6500.).and.(maxMF(i).gt.0.))ierr(i)=555 + endif enddo !$acc end kernels if (dx(its)<6500.) then diff --git a/physics/CONV/Grell_Freitas/cu_gf_deep.F90 b/physics/CONV/Grell_Freitas/cu_gf_deep.F90 index 74da8d2de..2fd963d57 100644 --- a/physics/CONV/Grell_Freitas/cu_gf_deep.F90 +++ b/physics/CONV/Grell_Freitas/cu_gf_deep.F90 @@ -5496,7 +5496,7 @@ subroutine get_partition_liq_ice(ierr,tn,po_cup, p_liq_ice,melting_layer ,itf,ktf,its,ite, kts,kte, cumulus ) implicit none character *(*), intent (in) :: cumulus - integer ,intent (in ) :: itf,ktf, its,ite, kts,kte + integer ,intent (in ) :: itf,ktf, its,ite, kts,kte real(kind=kind_phys), intent (in ), dimension(its:ite,kts:kte) :: tn,po_cup real(kind=kind_phys), intent (inout), dimension(its:ite,kts:kte) :: p_liq_ice,melting_layer !$acc declare copyin(tn,po_cup) copy(p_liq_ice,melting_layer) diff --git a/physics/GWD/cires_tauamf_data.F90 b/physics/GWD/cires_tauamf_data.F90 index 364c79409..323cea9a8 100644 --- a/physics/GWD/cires_tauamf_data.F90 +++ b/physics/GWD/cires_tauamf_data.F90 @@ -36,7 +36,7 @@ subroutine read_tau_amf(me, master, errmsg, errflg) if(iernc.ne.0) then write(errmsg,'(*(a))') "read_tau_amf: cannot open file_limb_tab data-file ", & trim(ugwp_taufile) - print *, 'cannot open ugwp-v1 tau-file=',trim(ugwp_taufile) + print *, 'cannot open ugwp-v1 tau-file=',trim(ugwp_taufile) errflg = 1 return else @@ -51,26 +51,26 @@ subroutine read_tau_amf(me, master, errmsg, errflg) status = nf90_inquire_dimension(ncid, DimID, len =ntau_d2t ) if (me == master) print *, ntau_d1y, ntau_d2t, ' dimd of tau_ngw ugwp-v1 ' - if (ntau_d2t .le. 0 .or. ntau_d1y .le. 0) then - print *, 'ugwp-v1 tau-file=', trim(ugwp_taufile) - print *, ' ugwp-v1: ', 'ntau_d2t=',ntau_d2t, 'ntau_d2t=',ntau_d1y - stop - endif + if (ntau_d2t .le. 0 .or. ntau_d1y .le. 0) then + print *, 'ugwp-v1 tau-file=', trim(ugwp_taufile) + print *, ' ugwp-v1: ', 'ntau_d2t=',ntau_d2t, 'ntau_d2t=',ntau_d1y + stop + endif if (.not.allocated(ugwp_taulat)) allocate (ugwp_taulat(ntau_d1y )) if (.not.allocated(days_limb)) allocate (days_limb(ntau_d2t)) - if (.not.allocated(tau_limb)) allocate (tau_limb(ntau_d1y, ntau_d2t )) + if (.not.allocated(tau_limb)) allocate (tau_limb(ntau_d1y, ntau_d2t )) - iernc=nf90_inq_varid( ncid, 'DAYS', vid ) + iernc=nf90_inq_varid( ncid, 'DAYS', vid ) iernc= nf90_get_var( ncid, vid, days_limb) - iernc=nf90_inq_varid( ncid, 'LATS', vid ) + iernc=nf90_inq_varid( ncid, 'LATS', vid ) iernc= nf90_get_var( ncid, vid, ugwp_taulat) - iernc=nf90_inq_varid( ncid, 'ABSMF', vid ) + iernc=nf90_inq_varid( ncid, 'ABSMF', vid ) iernc= nf90_get_var( ncid, vid, tau_limb) - iernc=nf90_close(ncid) + iernc=nf90_close(ncid) - endif + endif end subroutine read_tau_amf @@ -102,22 +102,22 @@ subroutine cires_indx_ugwp (npts, me, master, dlat,j1_tau,j2_tau, w1_j1tau, w2_j j2_tau(j) = min(j2_tau(j),ntau_d1y) - j1_tau(j) = max(j2_tau(j)-1,1) + j1_tau(j) = max(j2_tau(j)-1,1) if (j1_tau(j) /= j2_tau(j) ) then w2_j2tau(j) = (dlat(j) - ugwp_taulat(j1_tau(j))) & - / (ugwp_taulat(j2_tau(j))-ugwp_taulat(j1_tau(j))) + / (ugwp_taulat(j2_tau(j))-ugwp_taulat(j1_tau(j))) else w2_j2tau(j) = 1.0 endif - w1_j1tau(j) = 1.0 - w2_j2tau(j) + w1_j1tau(j) = 1.0 - w2_j2tau(j) enddo return end subroutine cires_indx_ugwp !> subroutine tau_amf_interp(me, master, im, idate, fhour, j1_tau,j2_tau, ddy_j1, ddy_j2, tau_ddd) - use machine, only: kind_phys + use machine, only: kind_phys implicit none !input @@ -141,30 +141,30 @@ subroutine tau_amf_interp(me, master, im, idate, fhour, j1_tau,j2_tau, ddy_j1, d it1 = 2 do iday=1, ntau_d2t - if (fddd .lt. days_limb(iday) ) then - it2 = iday - exit - endif - enddo + if (fddd .lt. days_limb(iday) ) then + it2 = iday + exit + endif + enddo - it2 = min(it2,ntau_d2t) - it1 = max(it2-1,1) - if (it2 > ntau_d2t ) then - print *, ' Error in time-interpolation for tau_amf_interp ' - print *, ' it1, it2, ntau_d2t ', it1, it2, ntau_d2t - print *, ' Error in time-interpolation see cires_tauamf_data.F90 ' - stop - endif + it2 = min(it2,ntau_d2t) + it1 = max(it2-1,1) + if (it2 > ntau_d2t ) then + print *, ' Error in time-interpolation for tau_amf_interp ' + print *, ' it1, it2, ntau_d2t ', it1, it2, ntau_d2t + print *, ' Error in time-interpolation see cires_tauamf_data.F90 ' + stop + endif - w2 = (fddd-days_limb(it1))/(days_limb(it2)-days_limb(it1)) - w1 = 1.0-w2 + w2 = (fddd-days_limb(it1))/(days_limb(it2)-days_limb(it1)) + w1 = 1.0-w2 - do i=1, im - j1 = j1_tau(i) - j2 = j2_tau(i) - tx1 = tau_limb(j1, it1)*ddy_j1(i)+tau_limb(j2, it1)*ddy_j2(i) - tx2 = tau_limb(j1, it2)*ddy_j1(i)+tau_limb(j2, it2)*ddy_j2(i) - tau_ddd(i) = tx1*w1 + w2*tx2 + do i=1, im + j1 = j1_tau(i) + j2 = j2_tau(i) + tx1 = tau_limb(j1, it1)*ddy_j1(i)+tau_limb(j2, it1)*ddy_j2(i) + tx2 = tau_limb(j1, it2)*ddy_j1(i)+tau_limb(j2, it2)*ddy_j2(i) + tau_ddd(i) = tx1*w1 + w2*tx2 enddo end subroutine tau_amf_interp @@ -172,7 +172,7 @@ end subroutine tau_amf_interp !> subroutine gfs_idate_calendar(idate, fhour, ddd, fddd) - use machine, only: kind_phys + use machine, only: kind_phys implicit none ! input integer, intent(in) :: idate(4) diff --git a/physics/GWD/cires_ugwpv1_oro.F90 b/physics/GWD/cires_ugwpv1_oro.F90 index 423a21348..12f9d0d13 100644 --- a/physics/GWD/cires_ugwpv1_oro.F90 +++ b/physics/GWD/cires_ugwpv1_oro.F90 @@ -222,7 +222,8 @@ subroutine orogw_v1 (im, km, imx, me, master, dtp, kdt, do_tofd, & dusfc(i) = 0.0 dvsfc(i) = 0.0 ipt(i) = 0 - enddo + enddo + zlwb(:) = 0.0 ! ---- for lm and gwd calculation points !cires_ugwp_initialize.F90: real, parameter :: hpmax=2400.0, hpmin=25.0 diff --git a/physics/GWD/cires_ugwpv1_triggers.F90 b/physics/GWD/cires_ugwpv1_triggers.F90 index 009d91775..29d66f823 100644 --- a/physics/GWD/cires_ugwpv1_triggers.F90 +++ b/physics/GWD/cires_ugwpv1_triggers.F90 @@ -298,7 +298,7 @@ subroutine get_spectra_tau_okw(nw, im, levs, trig_okw, xlatd, sinlat, coslat, t if (dmax >= tlim_okw) kex = kex+1 do k=klow+1, ktop dtot = abs(trig_okw(i,k)) - if (dtot >= tlim_fgf ) kex = kex+1 + if (dtot >= tlim_okw ) kex = kex+1 if ( dtot > dmax) then klev(i) = k dmax = dtot diff --git a/physics/GWD/ecmwf_ngw.F90 b/physics/GWD/ecmwf_ngw.F90 new file mode 100644 index 000000000..09d5b80cd --- /dev/null +++ b/physics/GWD/ecmwf_ngw.F90 @@ -0,0 +1,956 @@ +!>\file ecmwf_ngw.F90 +!! + +module ecmwf_ngw + + +contains +!------------------------------------------------------------------------------------------ +! April 2025 Adding ECMWF Non-stationary gravity wave scheme option by Bo Yang +!------------------------------------------------------------------------------------------ +! different orientation at vertical +! 1 is the highest level for ECMWF, 1 is the lowest level for GFS +! Vertical levels are reversed after entering the subroutine ecmwf_ngw_emc and reversed back before exiting +! This non-stationary GWD module was obtained by Fanglin Yang from ECMWF, with permission for operational use at NCEP. We would +! like to thank Andy Brown, Michael Sleigh, Peter Bechtold, and Nils Wedi at ECMWF for their support in porting this code to the +! UFS. +!! Original Fortran Code by J. SCINOCCIA +! Rewritten in IFS format by A. ORR E.C.M.W.F. August 2008 +! PURPOSE +! ------- + +! THIS ROUTINE COMPUTES NON-OROGRAPHIC GRAVITY WAVE DRAG +! AFTER SCINOCCA (2003) AND Mc LANDRESS AND SCINOCCIA (JAS 2005) +! HYDROSTATIC NON-ROTATIONAL SIMPLIFIED VERSION OF THE +! WARNER AND MCINTYRE (1996) NON-OROGRAPHIC GRAVITY WAVE PARAMETERIZATION +! CONSTANTS HAVE BEEN OPTIMIZED FOLLOWING M. ERN ET AL. (ATMOS. CHEM. PHYS. 2006) + +! REFERENCE: Orr, A., P. Bechtold, J. Scinoccia, M. Ern, M. Janiskova, 2010: +! Improved middle atmosphere climate and analysis in the ECMWF forecasting system +! through a non-orographic gravity wave parametrization. J. Climate., 23, 5905-5926. + +! LAUNCH SPECTRUM - GENERALIZED DESAUBIES +! INCLUDES A CRITICAL-LEVEL CORRECTION THAT PREVENTS THE +! MOMEMTUM DEPOSITION IN EACH AZIMUTH FROM DRIVING THE FLOW TO SPEEDS FASTER +! THAN THE PHASE SPEED OF THE WAVES, I.E. WHEN WAVES BREAK THEY DRAG THE MEAN +! FLOW TOWARDS THEIR PHASE SPEED - NOT PAST IT. +!!! different orientation for vertical +!!! 1 is the highest level for ECMWF, 1 is the lowest level for GFS +!--------------------------------------------------- +! subroutine cires_ugwpv1_ngw_solv2(mpi_id, master, im, levs, kdt, dtp, & +! tau_ngw, tm , um, vm, qm, prsl, prsi, zmet, zmeti, prslk, & +! xlatd, sinlat, coslat, & +! pdudt, pdvdt, pdtdt, dked, zngw) + +! (C) Copyright 1989- ECMWF. +! +! This software is licensed under the terms of the Apache Licence Version 2.0 +! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +! +! In applying this licence, ECMWF does not waive the privileges and immunities +! granted to it by virtue of its status as an intergovernmental organisation +! nor does it submit to any jurisdiction. + + subroutine ecmwf_ngw_emc(mpi_id, master, KLON, KLEV, kdt, PTSTEP, DX, & + tau_ngw, PTM11, PUM11, PVM11, qm1, PAPM11, PAPHM11, PGEO11, zmeti1, prslk1, & + xlatd, sinlat, coslat, & + PTENU, PTENV, pdtdt, dked, zngw) + + +! + use machine, only : kind_phys + + + use cires_ugwpv1_module,only : psrc => knob_ugwp_palaunch + + use cires_ugwpv1_module,only : maxdudt, maxdtdt, max_eps, dked_min, dked_max +!! maxdudt=250.e-5; maxdtdt=15.e-2; dked_min=0.01; dked_max=250.0 +!! max_eps=max_kdis*4.e-4=450*4.e-4 + + use ugwp_common , only : rgrav, grav, cpd, rd, rv, rcpdl, grav2cpd, & + omega2, rcpd, rcpd2, pi, pi2, fv, & + rad_to_deg, deg_to_rad, & + rdi, gor, grcp, gocp, & + bnv2min, bnv2max, dw2min, velmin, gr2, & + hpscale, rhp, rh4, grav2, rgrav2, mkzmin, mkz2min + +!!! grav=con_g; rgrav=1/grav; cpd=con_cp; rd=con_rd; rv=con_rv; rcpdl=cpd*rgrav=cpd/g +!!!grav2cpd= grav*grcp=grav*grav*rcpd= g**2/cpd omega=con_omega; omega2=2*omega1 +!! rcpd2=0.5*rcpd=1/(2*cpd) +!! pi=con_pi; pi2=2*pi +!! fv=con_fvirt; rad_to_deg=180.0/pi; deg_to_rad=pi2/180.0 +!! rdi= 1.0/rd; gor=grav/rd; grcp= grav*rcpd=g/Cpd; gocp=grcp=grav*rcpd= g/cpd +!! bnv2min=(pi2/1800.)*(pi2/1800.); bnv2max=(pi2/30.)*(pi2/30.) +!! dw2min=1.0, velmin=sqrt(dw2min) +!! gr2=grav*gor=grav*grav/rd=g**2/rd +!! hpscale=7000; rhp = 1./hpscale = 1/7000; rh4=rhp2*rhp2=(0.5*rhp)**2=1/4*1/7000=1/28000 +!! rgrav2=rgrav*rgrav=1/(g**2); grav2=grav+grav=2*g; mkzmin=pi2/80.0e-3=2*pi/80.0e3 +!! mkz2min= mkzmin*mkzmin=(2*pi/80.0E3)**2 + + + use ugwp_wmsdis_init, only : v_kxw, rv_kxw, v_kxw2, tamp_mpa, tau_min, ucrit, & + gw_eff, & +! zms, & + zci4, zci3, zci2, & + rimin, sc2, sc2u, ric +! v_kxw=kxw=pi2/lhmet=pi2/200e3; rv_kxw=200e3/pi2 +! v_kxw2=v_kxw*v_kxw=(pi2/200e3)**2; tamp_mpa=knob_ugwp_tauamp amplitude for GEOS-5/MERRA-2 +! tau_min=min of GW MF 0.25mPa +! ucrit=cdmin=2e-2/mkzmax=2e-2/((2pie)/500) = 10/(2pie) +! gw_eff=effac=1.0 +! zci4=(zms*zci(inc))**4; zci2=(zms*zci(inc))**2; zci3(inc)=(zms*zci(inc))**3 +! rimin=-10.0; sc2=lturb*lturb=30m*30m; sc2u=ulturb*ulturb=150*150 + + + implicit none + +!in +!work + + +! integer, intent(in) :: KLAUNCH ! index for launch level + integer, intent(in) :: KLEV ! vertical level + integer, intent(in) :: KLON ! horiz tiles + integer, intent(in) :: mpi_id, master, kdt + + real(kind=kind_phys) ,intent(in) :: PTSTEP ! model time step + real(kind=kind_phys) ,intent(in) :: DX(KLON) ! model grid size + + real(kind=kind_phys) ,intent(in) :: tau_ngw(KLON) + + real(kind=kind_phys) ,intent(in) :: PVM11(KLON,KLEV) ! meridional wind + real(kind=kind_phys) ,intent(in) :: PUM11(KLON,KLEV) ! zonal wind + real(kind=kind_phys) ,intent(in) :: qm1(KLON,KLEV) ! spec. humidity + real(kind=kind_phys) ,intent(in) :: PTM11(KLON,KLEV) ! kinetic temperature + real(kind=kind_phys) ,intent(in) :: PAPM11(KLON,KLEV) ! mid-layer pressure + real(kind=kind_phys) ,intent(in) :: PAPHM11(KLON,KLEV+1) ! interface pressure + real(kind=kind_phys) ,intent(in) :: PGEO11(KLON,KLEV) ! full model level geopotential in meters + + real(kind=kind_phys) :: PVM1(KLON,KLEV) ! meridional wind + real(kind=kind_phys) :: PUM1(KLON,KLEV) ! zonal wind + real(kind=kind_phys) :: qm(KLON,KLEV) ! spec. humidity + real(kind=kind_phys) :: PTM1(KLON,KLEV) ! kinetic temperature + real(kind=kind_phys) :: PAPM1(KLON,KLEV) ! mid-layer pressure + real(kind=kind_phys) :: PAPHM1(KLON,KLEV+1) ! interface pressure + real(kind=kind_phys) :: PGEO1(KLON,KLEV) ! full model level geopotential in meters + +! real(kind=kind_phys) :: PGAW(KLON) !normalised gaussian quadrature weight/nb of longitude pts + ! local sub-area == 4*RPI*RA**2 * PGAW +! real(kind=kind_phys) ,intent(in) :: PPRECIP(KLON) ! total surface precipitation + + + + + real(kind=kind_phys) ,intent(in) :: prslk1(KLON,KLEV) ! mid-layer exner function + real(kind=kind_phys) :: prslk(KLON,KLEV) ! mid-layer exner function +! real(kind=kind_phys) ,intent(in) :: zmet(KLON,KLEV) ! meters phil =philg/grav ! use PGEO1 instead + + real(kind=kind_phys) ,intent(in) :: zmeti1(KLON,KLEV+1) ! interface geopi/meters + real(kind=kind_phys) :: zmeti(KLON,KLEV+1) ! interface geopi/meters + real(kind=kind_phys) ,intent(in) :: xlatd(KLON) ! xlat_d in degrees + real(kind=kind_phys) ,intent(in) :: sinlat(KLON) + real(kind=kind_phys) ,intent(in) :: coslat(KLON) +! +! out-gw effects +! + real(kind=kind_phys) ,intent(out) :: PTENU(KLON,KLEV) ! zonal momentum tendency + real(kind=kind_phys) ,intent(out) :: PTENV(KLON,KLEV) ! meridional momentum tendency + real(kind=kind_phys) ,intent(out) :: pdtdt(KLON,KLEV) ! gw-heating (u*ax+v*ay)/cp and cooling + + + real(kind=kind_phys) :: PFLUXU(KLON,KLEV+1) ! zonal momentum tendency + real(kind=kind_phys) :: PFLUXV(KLON,KLEV+1) ! meridional momentum tendency + + real(kind=kind_phys) ,intent(out) :: dked(KLON,KLEV) ! gw-eddy diffusion + + real(kind=kind_phys) ,intent(out) :: zngw(KLON) ! launch height +! +!work + INTEGER, PARAMETER :: IAZIDIM=4 !number of azimuths + INTEGER, PARAMETER :: INCDIM=20 !number of discretized c spectral elements in launch spectrum + + REAL(kind=kind_phys), PARAMETER :: RA=6370.e3 !half-model level zonal velocity + REAL(kind=kind_phys), PARAMETER :: GPTWO=2.0 ! 2p in equation + + REAL(kind=kind_phys) :: ZUHM1(KLON,KLEV) !half-model level zonal velocity + REAL(kind=kind_phys) :: ZVHM1(KLON,KLEV) !half-model level meridional velocity + + REAL(kind=kind_phys) :: ZBVFHM1(KLON,KLEV) !half-model level Brunt-Vaisalla frequency + REAL(kind=kind_phys) :: ZRHOHM1(KLON,KLEV) !half-model level density + REAL(kind=kind_phys) :: ZX(INCDIM) !coordinate transformation + REAL(kind=kind_phys) :: ZCI(INCDIM) !phase speed element + REAL(kind=kind_phys) :: ZDCI(INCDIM) + REAL(kind=kind_phys) :: ZUI(KLON,KLEV,IAZIDIM) !intrinsic velocity + REAL(kind=kind_phys) :: ZUL(KLON,IAZIDIM) !velocity in azimuthal direction at launch level + REAL(kind=kind_phys) :: ZBVFL(KLON) !buoyancy at launch level + REAL(kind=kind_phys) :: ZCOSANG(IAZIDIM) !cos of azimuth angle + REAL(kind=kind_phys) :: ZSINANG(IAZIDIM) !sin of azimuth angle + REAL(kind=kind_phys) :: ZFCT(KLON,KLEV) + REAL(kind=kind_phys) :: ZFNORM(KLON) !normalisation factor (A) + REAL(kind=kind_phys) :: ZCI_MIN(KLON,IAZIDIM) + REAL(kind=kind_phys) :: ZTHM1(KLON,KLEV) !temperature on half-model levels + REAL(kind=kind_phys) :: ZFLUX(KLON,INCDIM,IAZIDIM) !momentum flux at each vertical level and azimuth + REAL(kind=kind_phys) :: ZPU(KLON,KLEV,IAZIDIM) !momentum flux + REAL(kind=kind_phys) :: ZDFL(KLON,KLEV,IAZIDIM) + REAL(kind=kind_phys) :: ZACT(KLON,INCDIM,IAZIDIM) !if =1 then critical level encountered + REAL(kind=kind_phys) :: ZACC(KLON,INCDIM,IAZIDIM) + REAL(kind=kind_phys) :: ZCRT(KLON,KLEV,IAZIDIM) + + INTEGER :: ILAUNCH !model level from which GW spectrum is launched + INTEGER :: INC, JK, JL, IAZI + + + REAL(KIND=kind_phys) :: ZRADTODEG, ZGELATDEG + REAL(KIND=kind_phys) :: ZCIMIN, ZCIMAX + REAL(KIND=kind_phys) :: ZGAM, ZPEXP, ZXMAX, ZXMIN, ZXRAN + REAL(KIND=kind_phys) :: ZDX + + REAL(KIND=kind_phys) :: ZX1, ZX2, ZDXA, ZDXB, ZDXS + REAL(KIND=kind_phys) :: ZANG, ZAZ_FCT, ZNORM, ZANG1, ZTX + REAL(KIND=kind_phys) :: ZU, ZCIN, ZCPEAK + REAL(KIND=kind_phys) :: ZCIN4, ZBVFL4, ZCIN2, ZBVFL2, ZCIN3, ZBVFL3, ZCINC + REAL(KIND=kind_phys) :: ZATMP, ZFLUXS, ZDEP, ZFLUXSQ, ZULM, ZDFT, ZE1, ZE2 + REAL(KIND=kind_phys) :: ZMS_L,ZMS, Z0P5, Z0P0, Z50S + REAL(KIND=kind_phys) :: ZGAUSS(KLON), ZFLUXLAUN(KLON), ZCNGL(KLON) + REAL(KIND=kind_phys) :: ZCONS1,ZCONS2,ZDELP,ZRGPTS + +!!! try to assign values for the following + + REAL(KIND=kind_phys) :: GSSEC + REAL(KIND=kind_phys) :: ZGAUSSB,ZFLUXGLOB + +! REAL(KIND=kind_phys) :: PGAW(KLON) don't need this value, set ZDX directly + + +! REAL(KIND=kind_phys) :: PGELAT(KLON) !!! use xlatd from gfs instead + REAL(KIND=kind_phys) :: GCOEFF, GGAUSSA !!! GCOEFF link to precip +! REAL(KIND=kind_phys) :: GGAUSSB !!! GGAUSSB->ZGAUSS + !!!! GGAUSSA + REAL(KIND=kind_phys) :: GCSTAR + + + LOGICAL :: LGACALC, LGSATL, LOZPR + + integer :: NGAUSS, NSLOPE + + + NSLOPE=1 + LGACALC=.false. + LGSATL=.false. + LOZPR=.true. +! NGAUSS=4 + NGAUSS=1 +! GGAUSSA=20._kind_phys +! GGAUSSA=10._kind_phys + GGAUSSA=5._kind_phys +! GGAUSSB=1.0_kind_phys +! ZGAUSSB=0.25_kind_phys +! ZGAUSSB=0.3_kind_phys +! ZGAUSSB=0.35_kind_phys + ZGAUSSB=0.38_kind_phys +! ZGAUSSB=-0.25_kind_phys +! ZGAUSSB=0.5_kind_phys +! ZGAUSSB=0.3_kind_phys + GCSTAR=1.0_kind_phys + + + +! GSSEC=(pi2/1800.)*(pi2/1800.) + GSSEC=1.e-24 + GCOEFF=1.0_kind_phys !!do not know the value, but never use it when NGAUSS is not equal 1 + + + +!! LOGIC :: LGINDL !!LGINDL=.true. using standard atm values to calculate!! comment out +!! REAL(KIND=kind_phys) :: STPHI(KM),STTEM(KM), STPREH(KM) +!! REAL(KIND=kind_phys) :: ZTHSTD,ZRHOSTD,ZBVFSTD + +! Set parameters which are a function of launch height +!!! need to set the parameter ILAUCH, ZFLUXGLOB,ZGAUSSB,ZMS_L directly +!ILAUNCH=NLAUNCHL(KLAUNCH) +!ZFLUXGLOB=GFLUXLAUNL(KLAUNCH) +!ZGAUSSB=GGAUSSB(KLAUNCH) +!ZMS_L=GMSTAR_L(KLAUNCH) + +!*INPUT PARAMETERS +!* ---------------- + ZRADTODEG=57.29577951_kind_phys + ZMS_L=2.e3_kind_phys +! ZFLUXGLOB=3.75e-3_kind_phys +! ZFLUXGLOB=3.7e-3_kind_phys +! ZFLUXGLOB=3.55e-3_kind_phys +! ZFLUXGLOB=3.65e-3_kind_phys !standard value +! ZFLUXGLOB=3.62e-3_kind_phys !standard value + ZFLUXGLOB=3.60e-3_kind_phys + +! ZFLUXGLOB=3.2e-3_kind_phys +! ZFLUXGLOB=3.0e-3_kind_phys +! ZFLUXGLOB=5.0e-3_kind_phys +! ZFLUXGLOB=4.0e-3_kind_phys +! ZFLUXGLOB=0.0_kind_phys + + ZMS=2*pi/ZMS_L + +!* INITIALIZE FIELDS TO ZERO +!* ------------------------- + + PTENU(:,:)=0.0_kind_phys + PTENV(:,:)=0.0_kind_phys + pdtdt(:,:)=0.0_kind_phys + dked(:,:)=0.0_kind_phys + + DO JK=1,KLEV+1 + DO JL=1,KLON + PFLUXU(JL,JK)=0.0_kind_phys + PFLUXV(JL,JK)=0.0_kind_phys + ENDDO + ENDDO + + + DO IAZI=1,IAZIDIM + DO JK=1,KLEV + DO JL=1,KLON + ZPU(JL,JK,IAZI)=0.0_kind_phys + ZCRT(JL,JK,IAZI)=0.0_kind_phys + ZDFL(JL,JK,IAZI)=0.0_kind_phys + ENDDO + ENDDO + ENDDO + +!!!!!!!!!!!!!!!!!!!!!!!! +! redefine ilaunch + + DO JK=1, KLEV + if (PAPM11(KLON,JK) .LT. psrc) exit + ENDDO + ILAUNCH = max(JK-1,3) + + DO JL=1,KLON + zngw(JL) = PGEO11(JL, ILAUNCH) + ENDDO + + ILAUNCH = KLEV + 1 - ILAUNCH + +!!!!!!!!!!!!!!!!!!!!!!!!!!! + + +!* reverse vertical coordinate to ECMWF + DO JL=1,KLON + PTM1(JL,:)=transfer(PTM11(JL,KLEV:1:-1),PTM11(JL,:)) + PUM1(JL,:)=transfer(PUM11(JL,KLEV:1:-1),PUM11(JL,:)) + PVM1(JL,:)=transfer(PVM11(JL,KLEV:1:-1),PVM11(JL,:)) + qm(JL,:)=transfer(qm1(JL,KLEV:1:-1),qm1(JL,:)) + PAPM1(JL,:)=transfer(PAPM11(JL,KLEV:1:-1),PAPM11(JL,:)) + PGEO1(JL,:)=transfer(PGEO11(JL,KLEV:1:-1),PGEO11(JL,:)) + prslk(JL,:)=transfer(prslk1(JL,KLEV:1:-1),prslk1(JL,:)) + + dked(JL,:)=transfer(dked(JL,KLEV:1:-1),dked(JL,:)) + PTENU(JL,:)=transfer(PTENU(JL,KLEV:1:-1),PTENU(JL,:)) + PTENV(JL,:)=transfer(PTENV(JL,KLEV:1:-1),PTENV(JL,:)) + pdtdt(JL,:)=transfer(pdtdt(JL,KLEV:1:-1),pdtdt(JL,:)) + + + PAPHM1(JL,:)=transfer(PAPHM11(JL,KLEV+1:1:-1),PAPHM11(JL,:)) + zmeti(JL,:)=transfer(zmeti1(JL,KLEV+1:1:-1),zmeti1(JL,:)) + ENDDO +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + + +!* INITIALIZE PARAMETERS FOR COORDINATE TRANSFORM +!* ---------------------------------------------- + +! ZCIMIN,ZCIMAX - min,max intrinsic launch-level phase speed (c-U_o) (m/s) +! ZGAM - half=width of coordinate stretch + + ZCIMIN=0.50_kind_phys + ZCIMAX=100.0_kind_phys + ZGAM=0.25_kind_phys + + ZPEXP=GPTWO/2.0_kind_phys + +! set initial min ci in each column and azimuth (used for critical levels) + + DO IAZI=1,IAZIDIM + DO JL=1,KLON + ZCI_MIN(JL,IAZI)=ZCIMIN + ENDDO + ENDDO + + +!* DEFINE HALF MODEL LEVEL WINDS AND TEMPERATURE +!* ----------------------------------- + + DO JK=2,KLEV + DO JL=1,KLON + ZTHM1(JL,JK) =0.5_kind_phys*(PTM1(JL,JK-1)+PTM1(JL,JK)) + ZUHM1(JL,JK) =0.5_kind_phys*(PUM1(JL,JK-1)+PUM1(JL,JK)) + ZVHM1(JL,JK) =0.5_kind_phys*(PVM1(JL,JK-1)+PVM1(JL,JK)) + ENDDO + ENDDO + + JK=1 + DO JL=1,KLON + ZTHM1(JL,JK)=PTM1(JL,JK) + ZUHM1(JL,JK)=PUM1(JL,JK) + ZVHM1(JL,JK)=PVM1(JL,JK) + ENDDO + + +!* DEFINE STATIC STABILITY AND AIR DENSITY ON HALF MODEL LEVELS +!* ------------------------------------------------------------ + +! ZCONS1=1.0_kind_phys/RD +! ZCONS2=RG**2/RCPD + + ZCONS1=1.0_kind_phys/rd + ZCONS2=grav2cpd + + + DO JK=KLEV,2,-1 + DO JL=1,KLON +! ZDELP=PAPM1(JL,JK)-PAPM1(JL,JK-1) + ZDELP=(PGEO1(JL,JK)-PGEO1(JL,JK-1))*grav !! times grav since import from zmet + ZRHOHM1(JL,JK)=PAPHM1(JL,JK)*ZCONS1/ZTHM1(JL,JK) + ZBVFHM1(JL,JK)=ZCONS2/ZTHM1(JL,JK)*& + & (1.0_kind_phys+cpd*(PTM1(JL,JK)-PTM1(JL,JK-1))/ZDELP) +! & (1.0_kind_phys-RCPD*ZRHOHM1(JL,JK)*(PTM1(JL,JK)-PTM1(JL,JK-1))/ZDELP) + ZBVFHM1(JL,JK)=MAX(ZBVFHM1(JL,JK),GSSEC) + ZBVFHM1(JL,JK)=SQRT(ZBVFHM1(JL,JK)) + ENDDO + ENDDO + +!* SET UP AZIMUTH DIRECTIONS AND SOME TRIG FACTORS +!* ----------------------------------------------- + + ZANG=2*pi/IAZIDIM + ZAZ_FCT=1.0_kind_phys + + +! get normalization factor to ensure that the same amount of momentum +! flux is directed (n,s,e,w) no mater how many azimuths are selected. +! note, however, the code below assumes a symmetric distribution of +! of azimuthal directions (ie 4,8,16,32,...) + + ZNORM=0.0_kind_phys + DO IAZI=1,IAZIDIM + ZANG1=(IAZI-1)*ZANG + ZCOSANG(IAZI)=COS(ZANG1) + ZSINANG(IAZI)=SIN(ZANG1) + ZNORM=ZNORM+ABS(ZCOSANG(IAZI)) + ENDDO + ZAZ_FCT=2._kind_phys*ZAZ_FCT/ZNORM + + +!* DEFINE COORDINATE TRANSFORM +!* ----------------------------------------------- + +! note that this is expresed in terms of the intrinsic phase speed +! at launch ci=c-u_o so that the transformation is identical at every +! launch site. +! See Eq. 28-30 of Scinocca 2003. + + ZXMAX=1.0_kind_phys/ZCIMIN + ZXMIN=1.0_kind_phys/ZCIMAX + + ZXRAN=ZXMAX-ZXMIN + ZDX=ZXRAN/REAL(INCDIM-1) + IF(LGACALC) ZGAM=(ZXMAX-ZXMIN)/LOG(ZXMAX/ZXMIN) +!! LGACALC=.false. ZGAM =0.25 ZX1=0.0007 +!! LGACALC=.true. ZGAM =0.37559 ZX1=0.01 + + ZX1=ZXRAN/(EXP(ZXRAN/ZGAM)-1.0_kind_phys) + ZX2=ZXMIN-ZX1 + + + + + DO INC=1,INCDIM + ZTX=REAL(INC-1)*ZDX+ZXMIN + ZX(INC)=ZX1*EXP((ZTX-ZXMIN)/ZGAM)+ZX2 !Eq. 29 of Scinocca 2003 + ZCI(INC)=1.0_kind_phys/ZX(INC) !Eq. 28 of Scinocca 2003 + ZDCI(INC)=ZCI(INC)**2*(ZX1/ZGAM)*EXP((ZTX-ZXMIN)/ZGAM)*ZDX !Eq. 30 of Scinocca 2003 + ENDDO + + +!* DEFINE INTRINSIC VELOCITY (RELATIVE TO LAUNCH LEVEL VELOCITY) U(Z)-U(Zo), AND COEFFICINETS +!* ------------------------------------------------------------------------------------------ + + DO IAZI=1,IAZIDIM + DO JL=1,KLON + ZUL(JL,IAZI)=ZCOSANG(IAZI)*ZUHM1(JL,ILAUNCH)& + & +ZSINANG(IAZI)*ZVHM1(JL,ILAUNCH) + ENDDO + ENDDO + DO JL=1,KLON + ZBVFL(JL)=ZBVFHM1(JL,ILAUNCH) + ENDDO + + DO JK=2,ILAUNCH + DO IAZI=1,IAZIDIM + DO JL=1,KLON + ZU=ZCOSANG(IAZI)*ZUHM1(JL,JK)+ZSINANG(IAZI)*ZVHM1(JL,JK) + ZUI(JL,JK,IAZI)=ZU-ZUL(JL,IAZI) + ENDDO + ENDDO + ENDDO + +!* DEFINE RHO(Zo)/N(Zo) +!* ------------------- + DO JK=2,ILAUNCH + DO JL=1,KLON + ZFCT(JL,JK)=ZRHOHM1(JL,JK)/ZBVFHM1(JL,JK) + ENDDO + ENDDO + +! Optionally set ZFCT at launch level using standard atmos values, to ensure saturation is +! independent of location +! IF (LGINDL) THEN +! ZCONS1=1.0_kind_phys/rd +! ZCONS2=grav2cpd +! ZDELP=STPHI(ILAUNCH)-STPHI(ILAUNCH-1) !!probably need to time grav depending on input +! THSTD=0.5_kind_phys*(STTEM(ILAUNCH-1)+STTEM(ILAUNCH)) +! ZRHOSTD=STPREH(ILAUNCH-1)*ZCONS1/ZTHSTD +! ZBVFSTD=ZCONS2/ZTHSTD*(1.0_kind_phys+rcpd* +! & (STTEM(ILAUNCH)-STTEM(ILAUNCH-1))/ZDELP) +! +! ZBVFSTD=MAX(ZBVFSTD,GSSEC) +! ZBVFSTD=SQRT(ZBVFSTD) +! DO JL=1,KLON +! ZFCT(JL,ILAUNCH)=ZRHOSTD/ZBVFSTD +! ENDDO +! ENDIF + +!* SET LAUNCH MOMENTUM FLUX SPECTRAL DENSITY +!* ----------------------------------------- + +! Eq. (25) of Scinocca 2003 (not including the 'A' component), and with U-Uo=0 +! do this for only one azimuth since it is identical to all azimuths, and it will be renormalized +! Initial spectrum fully saturated if LGSATL + + IF(NSLOPE==1) THEN +! s=1 case + DO INC=1,INCDIM + ZCIN=ZCI(INC) + ZCIN4=(ZMS*ZCIN)**4 + DO JL=1,KLON + ZBVFL4=ZBVFL(JL)**4 + IF(LGSATL) THEN + ZFLUX(JL,INC,1)=ZFCT(JL,ILAUNCH)*ZBVFL4*MIN(ZCIN/ZCIN4,& + & ZCIN/ZBVFL4) + ELSE + ZFLUX(JL,INC,1)=ZFCT(JL,ILAUNCH)*ZBVFL4*ZCIN/(ZBVFL4+ZCIN4) + ENDIF + ZACT(JL,INC,1)=1.0_kind_phys + ENDDO + ENDDO + + ELSEIF(NSLOPE==2) THEN +! s=2 case + DO INC=1,INCDIM + ZCIN=ZCI(INC) + ZCIN4=(ZMS*ZCIN)**4 + DO JL=1,KLON + ZBVFL4=ZBVFL(JL)**4 + ZCPEAK=ZBVFL(JL)/ZMS + IF(LGSATL) THEN + ZFLUX(JL,INC,1)=ZFCT(JL,ILAUNCH)*ZBVFL4*MIN& + & (ZCPEAK/ZCIN4,ZCIN/ZBVFL4) + ELSE + ZFLUX(JL,INC,1)=ZFCT(JL,ILAUNCH)*ZBVFL4*ZCIN*ZCPEAK/(ZBVFL4& + & *ZCPEAK+ZCIN4*ZCIN) + ENDIF + ZACT(JL,INC,1)=1.0_kind_phys + ENDDO + ENDDO + + ELSEIF(NSLOPE==-1) THEN +! s=-1 case + DO INC=1,INCDIM + ZCIN=ZCI(INC) + ZCIN2=(ZMS*ZCIN)**2 + DO JL=1,KLON + ZBVFL2=ZBVFL(JL)**2 + ZFLUX(JL,INC,1)=ZFCT(JL,ILAUNCH)*ZBVFL2*ZCIN/(ZBVFL2+ZCIN2) + ZACT(JL,INC,1)=1.0_kind_phys + ENDDO + ENDDO + + ELSEIF(NSLOPE==0) THEN +! s=0 case + DO INC=1,INCDIM + ZCIN=ZCI(INC) + ZCIN3=(ZMS*ZCIN)**3 + DO JL=1,KLON + ZBVFL3=ZBVFL(JL)**3 + ZFLUX(JL,INC,1)=ZFCT(JL,ILAUNCH)*ZBVFL3*ZCIN/(ZBVFL3+ZCIN3) + ZACT(JL,INC,1)=1.0_kind_phys + ZACC(JL,INC,1)=1.0_kind_phys + ENDDO + ENDDO + + ENDIF + +!* NORMALIZE LAUNCH MOMENTUM FLUX +!* ------------------------------ + +! (rho x F^H = rho_o x F_p^total) + +! integrate (ZFLUX x dX) + DO INC=1,INCDIM + ZCINC=ZDCI(INC) + DO JL=1,KLON + ZPU(JL,ILAUNCH,1)=ZPU(JL,ILAUNCH,1)+ZFLUX(JL,INC,1)*ZCINC + ENDDO + ENDDO + + +!* NORMALIZE GFLUXLAUN TO INCLUDE SENSITIVITY TO PRECIPITATION +!* ----------------------------------------------------------- + +! Also other options to alter tropical values + +! A=ZFNORM in Scinocca 2003. A is independent of height. + ZDXA=1.0_kind_phys/29.E3_kind_phys + ZDXB=1.0_kind_phys/3.5E3_kind_phys + DO JL=1,KLON +!! ZDX=MAX(1.E2_kind_phys,2*RA*SQRT(pi*PGAW(JL))) !grid resolution (m) +!! ZDX=50.E3 ! c192 for c192 usage +!!! ZDX=25.E3 !C384 +!!! ZDX=13.E3 !C768 + + + !Scaling factor for launch flux depending on grid resolution + ! smooth reduction below 30 km + ZDXS=1.0_kind_phys-MIN(1.0_kind_phys,ATAN((MAX(1.0_kind_phys& + & /DX(JL),ZDXA)-ZDXA)/(ZDXB-ZDXA))) + ZFLUXLAUN(JL)=ZFLUXGLOB*ZDXS + ZFNORM(JL)=ZFLUXLAUN(JL)/ZPU(JL,ILAUNCH,1) + ENDDO + +! If LOZPR=TRUE then vary EPLAUNCH over tropics + IF (LOZPR) THEN + IF (NGAUSS==1) THEN + Z50S=-50.0_kind_phys + + DO JL=1,KLON + +! ZFLUXLAUN(JL)=ZFLUXLAUN(JL)*(1.0_kind_phys+MIN& +! & (0.5_kind_phys,GCOEFF*PPRECIP(JL))) !precip +! ZFNORM(JL)=ZFLUXLAUN(JL)/ZPU(JL,ILAUNCH,1) + + ZGELATDEG=xlatd(JL)-Z50S + ZGAUSS(JL)=ZGAUSSB*EXP((-ZGELATDEG*ZGELATDEG)& + & /(2*GGAUSSA*GGAUSSA)) + +! ZGELATDEG=xlatd(JL) +! ZGAUSS(JL)=-0.1_kind_phys*EXP((-ZGELATDEG*ZGELATDEG)& +! & /(2*GGAUSSA*GGAUSSA))+ZGAUSS(JL) + + +! ZGELATDEG=xlatd(JL) +! ZGAUSS(JL)=-0.05_kind_phys*EXP((-ZGELATDEG*ZGELATDEG)& +! & /(2*GGAUSSA*GGAUSSA))+ZGAUSS(JL) + + ZGELATDEG=xlatd(JL) + ZGAUSS(JL)=0.08_kind_phys*EXP((-ZGELATDEG*ZGELATDEG)& + & /(2*GGAUSSA*GGAUSSA))+ZGAUSS(JL) + + ZGELATDEG=xlatd(JL)-50.0_kind_phys + ZGAUSS(JL)= 0.0022_kind_phys*EXP((-ZGELATDEG*ZGELATDEG)& + & /(2*10.*10.))+ZGAUSS(JL) + + + ZFLUXLAUN(JL)=(1.0_kind_phys+ZGAUSS(JL))*ZFLUXLAUN(JL) + + + + ZFNORM(JL)=ZFLUXLAUN(JL)/ZPU(JL,ILAUNCH,1) + + + ENDDO + ELSEIF (NGAUSS==2) THEN + + DO JL=1,KLON +! ZGELATDEG=PGELAT(JL)*ZRADTODEG + ZGELATDEG=xlatd(JL) + ZGAUSS(JL)=ZGAUSSB*EXP((-ZGELATDEG*ZGELATDEG)& + & /(2*GGAUSSA*GGAUSSA)) + ZFLUXLAUN(JL)=(1.0_kind_phys+ZGAUSS(JL))*ZFLUXLAUN(JL) + ZFNORM(JL)=ZFLUXLAUN(JL)/ZPU(JL,ILAUNCH,1) + ENDDO + + + + ELSEIF (NGAUSS==4) THEN + +! Set latitudinal dependence to optimize stratospheric winds for 36r1 + Z50S=-50.0_kind_phys + DO JL=1,KLON +! ZGELATDEG=PGELAT(JL)*ZRADTODEG-Z50S + ZGELATDEG=xlatd(JL)-Z50S + ZGAUSS(JL)=ZGAUSSB*EXP((-ZGELATDEG*ZGELATDEG)/(2*GGAUSSA*GGAUSSA)) + ZFLUXLAUN(JL)=(1.0_kind_phys+ZGAUSS(JL))*ZFLUXLAUN(JL) + ZFNORM(JL)=ZFLUXLAUN(JL)/ZPU(JL,ILAUNCH,1) + ENDDO + + ENDIF + ENDIF + + DO IAZI=1,IAZIDIM + DO JL=1,KLON + ZPU(JL,ILAUNCH,IAZI)=ZFLUXLAUN(JL) + ENDDO + ENDDO + +!* ADJUST CONSTANT ZFCT +!* -------------------- + DO JK=2,ILAUNCH + DO JL=1,KLON + ZFCT(JL,JK)=ZFNORM(JL)*ZFCT(JL,JK) + ENDDO + ENDDO + +!* RENORMALIZE EACH SPECTRAL ELEMENT IN FIRST AZIMUTH +!* -------------------------------------------------- + DO INC=1,INCDIM + DO JL=1,KLON + ZFLUX(JL,INC,1)=ZFNORM(JL)*ZFLUX(JL,INC,1) + ENDDO + ENDDO + + + +!* COPY ZFLUX INTO ALL OTHER AZIMUTHS +!* -------------------------------- + +! ZACT=1 then no critical level +! ZACT=0 then critical level + + DO IAZI=2,IAZIDIM + DO INC=1,INCDIM + DO JL=1,KLON + ZFLUX(JL,INC,IAZI)=ZFLUX(JL,INC,1) + ZACT(JL,INC,IAZI)=1.0_kind_phys + ZACC(JL,INC,IAZI)=1.0_kind_phys + ENDDO + ENDDO + ENDDO + +! ----------------------------------------------------------------------------- + +!* BEGIN MAIN LOOP OVER LEVELS +!* --------------------------- + +!* begin IAZIDIM do-loop +!* -------------------- + + DO IAZI=1,IAZIDIM + +!* begin JK do-loop +!* ---------------- + + DO JK=ILAUNCH-1,2,-1 + +!* first do critical levels +!* ------------------------ + + DO JL=1,KLON + ZCI_MIN(JL,IAZI)=MAX(ZCI_MIN(JL,IAZI),ZUI(JL,JK,IAZI)) + ENDDO + +!* set ZACT to zero if critical level encountered +!* ---------------------------------------------- + + Z0P5=0.5_kind_phys + DO INC=1,INCDIM + ZCIN=ZCI(INC) + DO JL=1,KLON + ZATMP=Z0P5+SIGN(Z0P5,ZCIN-ZCI_MIN(JL,IAZI)) + ZACC(JL,INC,IAZI)=ZACT(JL,INC,IAZI)-ZATMP + ZACT(JL,INC,IAZI)=ZATMP + ENDDO + ENDDO + +!* integrate to get critical-level contribution to mom deposition on this level, i.e. ZACC=1 +!* ---------------------------------------------------------------------------------------- + + DO INC=1,INCDIM + ZCINC=ZDCI(INC) + DO JL=1,KLON + ZDFL(JL,JK,IAZI)=ZDFL(JL,JK,IAZI)+& + & ZACC(JL,INC,IAZI)*ZFLUX(JL,INC,IAZI)*ZCINC + ENDDO + ENDDO + +!* get weighted average of phase speed in layer + + DO JL=1,KLON + IF(ZDFL(JL,JK,IAZI)>0.0_kind_phys) THEN + ZATMP=ZCRT(JL,JK,IAZI) + DO INC=1,INCDIM + ZATMP=ZATMP+ZCI(INC)*& + & ZACC(JL,INC,IAZI)*ZFLUX(JL,INC,IAZI)*ZDCI(INC) + ENDDO + ZCRT(JL,JK,IAZI)=ZATMP/ZDFL(JL,JK,IAZI) + ELSE + ZCRT(JL,JK,IAZI)=ZCRT(JL,JK+1,IAZI) + ENDIF + ENDDO + +!* do saturation (Eq. (26) and (27) of Scinocca 2003) +!* ------------------------------------------------- + + IF(GPTWO==3.0_kind_phys) THEN + DO INC=1,INCDIM + ZCIN=ZCI(INC) + ZCINC=1.0_kind_phys/ZCIN + DO JL=1,KLON + ZE1=ZCIN-ZUI(JL,JK,IAZI) + ZE2=GCSTAR*ZFCT(JL,JK)*ZE1 + ZFLUXSQ=ZE2*ZE2*ZE1*ZCINC + ! ZFLUXSQ=ZE2*ZE2*ZE1/ZCIN + ZDEP=ZACT(JL,INC,IAZI)*(ZFLUX(JL,INC,IAZI)**2-ZFLUXSQ) + IF(ZDEP>0.0_kind_phys) THEN + ZFLUX(JL,INC,IAZI)=SQRT(ZFLUXSQ) + ENDIF + ENDDO + ENDDO + ELSEIF(GPTWO==2.0_kind_phys) THEN + DO INC=1,INCDIM + ZCIN=ZCI(INC) + ZCINC=1.0_kind_phys/ZCIN + DO JL=1,KLON + ZFLUXS=GCSTAR*ZFCT(JL,JK)*& + & (ZCIN-ZUI(JL,JK,IAZI))**2*ZCINC + ! ZFLUXS=GCSTAR*ZFCT(JL,JK)*(ZCIN-ZUI(JL,JK,IAZI))**2/ZCIN + ZDEP=ZACT(JL,INC,IAZI)*(ZFLUX(JL,INC,IAZI)-ZFLUXS) + IF(ZDEP>0.0_kind_phys) THEN + ZFLUX(JL,INC,IAZI)=ZFLUXS + ENDIF + ENDDO + ENDDO + ENDIF + +!* integrate spectrum +!* ------------------ + + DO INC=1,INCDIM + ZCINC=ZDCI(INC) + DO JL=1,KLON + ZPU(JL,JK,IAZI)=ZPU(JL,JK,IAZI)+& + & ZACT(JL,INC,IAZI)*ZFLUX(JL,INC,IAZI)*ZCINC + ENDDO + ENDDO + +!* end JK do-loop +!* -------------- + + ENDDO +!* end IAZIDIM do-loop +!* --------------- + + ENDDO + +! ----------------------------------------------------------------------------- + +!* MAKE CORRECTION FOR CRITICAL-LEVEL MOMENTUM DEPOSITION +!* ------------------------------------------------------ + + Z0P0=0._kind_phys +! ZRGPTS=1.0_kind_phys/(RG*PTSTEP) + ZRGPTS=1.0_kind_phys/(grav*PTSTEP) + DO IAZI=1,IAZIDIM + DO JL=1,KLON + ZCNGL(JL)=0.0_kind_phys + ENDDO + DO JK=2,ILAUNCH + DO JL=1,KLON + ZULM=ZCOSANG(IAZI)*PUM1(JL,JK)+ZSINANG(IAZI)*& + & PVM1(JL,JK)-ZUL(JL,IAZI) + ZDFL(JL,JK-1,IAZI)=ZDFL(JL,JK-1,IAZI)+ZCNGL(JL) + ZDFT=MIN(ZDFL(JL,JK-1,IAZI),2.0_kind_phys*(PAPM1(JL,JK-1)-& + & PAPM1(JL,JK))*(ZCRT(JL,JK-1,IAZI)-ZULM)*ZRGPTS) + + ZDFT=MAX(ZDFT,Z0P0) + ZCNGL(JL)=(ZDFL(JL,JK-1,IAZI)-ZDFT) + ZPU(JL,JK,IAZI)=ZPU(JL,JK,IAZI)-ZCNGL(JL) + ENDDO + ENDDO + ENDDO + + +!* SUM CONTRIBUTION FOR TOTAL ZONAL AND MERIDIONAL FLUX +!* --------------------------------------------------- + + DO IAZI=1,IAZIDIM + DO JK=ILAUNCH,2,-1 + DO JL=1,KLON + PFLUXU(JL,JK)=PFLUXU(JL,JK)+ZPU(JL,JK,IAZI)*ZAZ_FCT*ZCOSANG(IAZI) + PFLUXV(JL,JK)=PFLUXV(JL,JK)+ZPU(JL,JK,IAZI)*ZAZ_FCT*ZSINANG(IAZI) + ENDDO + ENDDO + ENDDO + + +!* UPDATE U AND V TENDENCIES +!* ---------------------------- + +! ZCONS1=1.0_kind_phys/RCPD + ZCONS1=rcpd + DO JK=1,ILAUNCH + DO JL=1, KLON +! ZDELP= RG/(PAPHM1(JL,JK+1)-PAPHM1(JL,JK)) + ZDELP= grav/(PAPHM1(JL,JK+1)-PAPHM1(JL,JK)) + ZE1=(PFLUXU(JL,JK+1)-PFLUXU(JL,JK))*ZDELP + ZE2=(PFLUXV(JL,JK+1)-PFLUXV(JL,JK))*ZDELP + + if (abs(ZE1) >= maxdudt ) then + ZE1 = sign(maxdudt, ZE1) + endif + if (abs(ZE2) >= maxdudt ) then + ZE2 = sign(maxdudt, ZE2) + endif + + + PTENU(JL,JK)=ZE1 + PTENV(JL,JK)=ZE2 +! add the tendency of dT/dt + ZE2=-(PUM1(JL,JK)*PTENU(JL,JK)+PVM1(JL,JK)*PTENV(JL,JK))/cpd + if (abs(ZE2) >= max_eps) pdtdt(JL,JK) = sign(max_eps, ZE2) + +! end of the tendency of dT/dt + ENDDO + ENDDO + + +!* reverse vertical coordinate back to GFS for outbound variables only + DO JL=1,KLON + dked(JL,:)=transfer(dked(JL,KLEV:1:-1),dked(JL,:)) + PTENU(JL,:)=transfer(PTENU(JL,KLEV:1:-1),PTENU(JL,:)) + PTENV(JL,:)=transfer(PTENV(JL,KLEV:1:-1),PTENV(JL,:)) + pdtdt(JL,:)=transfer(pdtdt(JL,KLEV:1:-1),pdtdt(JL,:)) + ENDDO + + +! DO JL=1,KLON +! dked(JL,:)=0.0 +! PTENU(JL,:)=0.0 +! PTENV(JL,:)=0.0 +! pdtdt(JL,:)=0.0 +! ENDDO + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + + return + end subroutine ecmwf_ngw_emc + + +end module ecmwf_ngw + + + + diff --git a/physics/GWD/ugwp_driver_v0.F b/physics/GWD/ugwp_driver_v0.F index 0fc4a2274..845323e43 100644 --- a/physics/GWD/ugwp_driver_v0.F +++ b/physics/GWD/ugwp_driver_v0.F @@ -348,7 +348,9 @@ SUBROUTINE GWDPS_V0(IM, km, imx, do_tofd, BVF2 = grav2 * RDZ * (VTK(I,K+1)-VTK(I,K)) & / (VTK(I,K+1)+VTK(I,K)) bnv2(i,k+1) = max( BVF2, bnv2min ) - RI_N(I,K+1) = Bnv2(i,k)/SHR2 ! Richardson number consistent with BNV2 + ! https://github.com/NCAR/ccpp-physics/issues/1103 + !RI_N(I,K+1) = Bnv2(i,k)/SHR2 ! Richardson number consistent with BNV2 + RI_N(I,K+1) = Bnv2(i,max(k,2))/SHR2 ! Richardson number consistent with BNV2 ! ! add here computation for Ktur and OGW-dissipation fro VE-GFS ! @@ -814,7 +816,7 @@ SUBROUTINE GWDPS_V0(IM, km, imx, do_tofd, IF( do_tofd ) then axtms(:,:) = 0.0 ; aytms(:,:) = 0.0 - + DO I = 1,npt J = ipt(i) zpbl =rgrav*phil( j, kpbl(j) ) diff --git a/physics/GWD/ugwpv1_gsldrag.F90 b/physics/GWD/ugwpv1_gsldrag.F90 index 094588a35..459e29362 100644 --- a/physics/GWD/ugwpv1_gsldrag.F90 +++ b/physics/GWD/ugwpv1_gsldrag.F90 @@ -43,6 +43,9 @@ module ugwpv1_gsldrag use cires_ugwpv1_module, only: cires_ugwpv1_init, ngwflux_update, calendar_ugwp use cires_ugwpv1_module, only: knob_ugwp_version, cires_ugwp_dealloc, tamp_mpa use cires_ugwpv1_solv2, only: cires_ugwpv1_ngw_solv2 +! use cires_ugwpv1_solv2, only: cires_ugwpv1_ngw_solv2, ecmwf_ngw + use ecmwf_ngw, only: ecmwf_ngw_emc + use cires_ugwpv1_oro, only: orogw_v1 use drag_suite, only: drag_suite_run, drag_suite_psl @@ -70,7 +73,8 @@ subroutine ugwpv1_gsldrag_init ( & con_pi, con_rerth, con_p0, & con_g, con_omega, con_cp, con_rd, con_rv,con_fvirt, & do_ugwp,do_ugwp_v0, do_ugwp_v0_orog_only, do_gsl_drag_ls_bl, & - do_gsl_drag_ss, do_gsl_drag_tofd, do_ugwp_v1, & + do_gsl_drag_ss, do_gsl_drag_tofd, do_ngw_ec, do_ugwp_v1, & +!! do_gsl_drag_ss, do_gsl_drag_tofd, do_ugwp_v1, & do_ugwp_v1_orog_only, do_ugwp_v1_w_gsldrag, errmsg, errflg) use ugwp_common @@ -94,9 +98,10 @@ subroutine ugwpv1_gsldrag_init ( & real(kind=kind_phys), intent (in) :: con_g, con_cp, con_rd, con_rv, con_omega, con_fvirt logical, intent (in) :: do_ugwp - logical, intent (in) :: do_ugwp_v0, do_ugwp_v0_orog_only, & - do_gsl_drag_ls_bl, do_gsl_drag_ss, & - do_gsl_drag_tofd, do_ugwp_v1, & + logical, intent (in) :: do_ugwp_v0, do_ugwp_v0_orog_only, & + do_gsl_drag_ls_bl, do_gsl_drag_ss, & +!! do_gsl_drag_tofd, do_ugwp_v1, & + do_gsl_drag_tofd, do_ugwp_v1, do_ngw_ec, & do_ugwp_v1_orog_only,do_ugwp_v1_w_gsldrag character(len=*), intent (in) :: fn_nml2 @@ -302,7 +307,7 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp, fhzero, kdt, ldiag3d, lssav, flag_for_gwd_generic_tend, do_gsl_drag_ls_bl, & do_gsl_drag_ss, do_gsl_drag_tofd, & do_gwd_opt_psl, psl_gwd_dx_factor, & - do_ugwp_v1, do_ugwp_v1_orog_only, & + do_ngw_ec, do_ugwp_v1, do_ugwp_v1_orog_only, & do_ugwp_v1_w_gsldrag, gwd_opt, do_tofd, ldiag_ugwp, ugwp_seq_update, & cdmbgwd, alpha_fd, jdat, nmtvr, hprime, oc, theta, sigma, gamma, & elvmax, clx, oa4, varss,oc1ss,oa4ss,ol4ss, dx, xlat, xlat_d, sinlat, coslat, & @@ -356,7 +361,7 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp, ! flags for choosing combination of GW drag schemes to run logical, intent (in) :: do_gsl_drag_ls_bl, do_gsl_drag_ss, do_gsl_drag_tofd - logical, intent (in) :: do_ugwp_v1, do_ugwp_v1_orog_only, do_tofd + logical, intent (in) :: do_ugwp_v1, do_ngw_ec, do_ugwp_v1_orog_only, do_tofd logical, intent (in) :: ldiag_ugwp, ugwp_seq_update logical, intent (in) :: do_ugwp_v1_w_gsldrag ! combination of ORO and NGW schemes @@ -698,10 +703,21 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp, call ngwflux_update(me, master, im, levs, kdt, ddd_ugwp,curdate, & tau_amf, xlat_d, sinlat,coslat, rain, tau_ngw) - call cires_ugwpv1_ngw_solv2(me, master, im, levs, kdt, dtp, & + if (do_ngw_ec) then + + call ecmwf_ngw_emc(me, master, im, levs, kdt, dtp, dx, & + tau_ngw, tgrs, ugrs, vgrs, q1, prsl, prsi, & + zmet, zmeti,prslk, xlat_d, sinlat, coslat, & + dudt_ngw, dvdt_ngw, dtdt_ngw, kdis_ngw, zngw) + else + + call cires_ugwpv1_ngw_solv2(me, master, im, levs, kdt, dtp, & tau_ngw, tgrs, ugrs, vgrs, q1, prsl, prsi, & zmet, zmeti,prslk, xlat_d, sinlat, coslat, & dudt_ngw, dvdt_ngw, dtdt_ngw, kdis_ngw, zngw) + + endif + ! ! => con_g, con_cp, con_rd, con_rv, con_omega, con_pi, con_fvirt ! diff --git a/physics/GWD/ugwpv1_gsldrag.meta b/physics/GWD/ugwpv1_gsldrag.meta index 376a562bb..24d8b0688 100644 --- a/physics/GWD/ugwpv1_gsldrag.meta +++ b/physics/GWD/ugwpv1_gsldrag.meta @@ -2,7 +2,7 @@ name = ugwpv1_gsldrag type = scheme dependencies = ../hooks/machine.F,drag_suite.F90 - dependencies = cires_ugwpv1_module.F90,cires_ugwpv1_triggers.F90,cires_ugwpv1_initialize.F90,cires_ugwpv1_solv2.F90 + dependencies = cires_ugwpv1_module.F90,cires_ugwpv1_triggers.F90,cires_ugwpv1_initialize.F90,cires_ugwpv1_solv2.F90,ecmwf_ngw.F90 dependencies = cires_ugwpv1_sporo.F90,cires_ugwpv1_oro.F90 ######################################################################## [ccpp-arg-table] @@ -218,6 +218,13 @@ dimensions = () type = logical intent = in +[do_ngw_ec] + standard_name = flag_for_ngw_ec + long_name = flag to activate ecmwf ngwd + units = flag + dimensions = () + type = logical + intent = in [do_ugwp_v1] standard_name = flag_for_ugwp_version_1 long_name = flag to activate ver 1 CIRES UGWP @@ -417,6 +424,13 @@ type = real kind = kind_phys intent = in +[do_ngw_ec] + standard_name = flag_for_ngw_ec + long_name = flag to activate ecmwf ngwd + units = flag + dimensions = () + type = logical + intent = in [do_ugwp_v1] standard_name = flag_for_ugwp_version_1 long_name = flag to activate ver 1 CIRES UGWP diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_3.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_3.meta index c79182679..c627b2202 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_3.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_3.meta @@ -286,7 +286,7 @@ dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys - intent = out + intent = inout optional = True [imp_physics] standard_name = control_for_microphysics_scheme diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/module_ccpp_suite_simulator.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/module_ccpp_suite_simulator.F90 index 45d3dd4e0..56d1d0666 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/module_ccpp_suite_simulator.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/module_ccpp_suite_simulator.F90 @@ -121,7 +121,7 @@ function linterp_1D(this, var_name, year, month, day, hour, min, sec) result(err this%tend1d%q = this%tend2d%q(:,1) endif end select - + err_message = "" end function linterp_1D !> Type-bound procedure to compute tendency profile for time-of-day. @@ -153,6 +153,7 @@ function linterp_2D(this, var_name, lon, lat, year, month, day, hour, min, sec) case("q") this%tend1d%q = w1*this%tend3d%q(iNearest,:,ti(1)) + w2*this%tend3d%q(iNearest,:,tf(1)) end select + err_message = "" end function linterp_2D !> Type-bound procedure to find nearest location. diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/sfcsub.F b/physics/Interstitials/UFS_SCM_NEPTUNE/sfcsub.F index 78afbac42..d4dadfe0d 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/sfcsub.F +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/sfcsub.F @@ -3093,8 +3093,11 @@ subroutine subst(data,imax,jmax,dlon,dlat,ijordr) enddo enddo else - do i=1,imax - data(imax-i+1,jj) = work(i,j) + do j=1,jmax + jj = jmax - j + 1 + do i=1,imax + data(imax-i+1,jj) = work(i,j) + enddo enddo endif else diff --git a/physics/MP/GFDL/GFDL_parse_tracers.F90 b/physics/MP/GFDL/GFDL_parse_tracers.F90 deleted file mode 100644 index 670c292ee..000000000 --- a/physics/MP/GFDL/GFDL_parse_tracers.F90 +++ /dev/null @@ -1,45 +0,0 @@ -!>\file GFDL_parse_tracers.F90 -!! - -!> This module contains code to parse tracers in GFDL MP -module parse_tracers - - integer, parameter :: NO_TRACER = -99 - - public get_tracer_index, NO_TRACER - -CONTAINS - - function get_tracer_index (tracer_names, name, me, master, debug) - - character(len=32), intent(in) :: tracer_names(:) - character(len=*), intent(in) :: name - integer, intent(in) :: me - integer, intent(in) :: master - logical, intent(in) :: debug - !--- local variables - integer :: get_tracer_index - integer :: i - - get_tracer_index = NO_TRACER - - do i=1, size(tracer_names) - if (trim(name) == trim(tracer_names(i))) then - get_tracer_index = i - exit - endif - enddo - - if (debug .and. (me == master)) then - if (get_tracer_index == NO_TRACER) then - print *,' PE ',me,' tracer with name '//trim(name)//' not found' - else - print *,' PE ',me,' tracer FOUND:',trim(name) - endif - endif - - return - - end function get_tracer_index - -end module parse_tracers diff --git a/physics/MP/GFDL/module_gfdl_cloud_microphys.F90 b/physics/MP/GFDL/module_gfdl_cloud_microphys.F90 index 72f3211b5..09e3c4b31 100644 --- a/physics/MP/GFDL/module_gfdl_cloud_microphys.F90 +++ b/physics/MP/GFDL/module_gfdl_cloud_microphys.F90 @@ -57,7 +57,7 @@ module gfdl_cloud_microphys_mod logical :: module_is_initialized = .false. logical :: qsmith_tables_initialized = .false. - character (len = 17) :: mod_name = 'gfdl_cloud_microphys' + character (len = 20) :: mod_name = 'gfdl_cloud_microphys' real, parameter :: n0r = 8.0e6, n0s = 3.0e6, n0g = 4.0e6 real, parameter :: rhos = 0.1e3, rhog = 0.4e3 diff --git a/physics/MP/NSSL/module_mp_nssl_2mom.F90 b/physics/MP/NSSL/module_mp_nssl_2mom.F90 index 963d2531c..1afbc13cf 100644 --- a/physics/MP/NSSL/module_mp_nssl_2mom.F90 +++ b/physics/MP/NSSL/module_mp_nssl_2mom.F90 @@ -1275,7 +1275,7 @@ END SUBROUTINE nssl_2mom_init_const !! NSSL MP setup routine (sets local options and array indices) SUBROUTINE nssl_2mom_init( & & ims,ime, jms,jme, kms,kme, nssl_params, ipctmp, mixphase,ihvol,idoniconlytmp, & - & namelist_filename, & + & namelist_filename, internal_nml, & & nssl_graupelfallfac, & & nssl_hailfallfac, & & nssl_ehw0, & @@ -1319,7 +1319,8 @@ SUBROUTINE nssl_2mom_init( & integer, intent(in),optional :: infileunit - integer,parameter::strsize=512 + integer,parameter::strsize=64 + character(len=*), intent(in), optional :: internal_nml(:) character(len=strsize), intent(in), optional :: namelist_filename character(len=strsize) :: namelist_inputfile @@ -1474,22 +1475,24 @@ SUBROUTINE nssl_2mom_init( & irainbreak = 0 ENDIF ENDIF - - - - + +#ifdef INTERNAL_FILE_NML + read (internal_nml, nml = nssl_mp_params, iostat=istat) +#else + namelist_inputfile = 'namelist.input' ! default for WRF/cm1 - IF ( present( namelist_filename ) ) THEN + IF ( present( namelist_filename ) ) THEN namelist_inputfile = namelist_filename ELSE namelist_inputfile = 'input.nml' ENDIF - - IF ( .true. ) THEN ! set to true to enable internal namelist read - open(15,file=namelist_inputfile,status='old',form='formatted',action='read') - rewind(15) - read(15,NML=nssl_mp_params,iostat=istat) - close(15) + + open(15,file=trim(namelist_inputfile),status='old',form='formatted',action='read') + rewind(15) + read(15,NML=nssl_mp_params,iostat=istat) + close(15) +#endif + IF ( .true. ) THEN ! set to true to enable internal namelist read IF ( present ( myrank ) .and. present ( mpiroot ) ) THEN IF ( myrank == mpiroot ) THEN IF ( istat /= 0 ) THEN @@ -1501,10 +1504,9 @@ SUBROUTINE nssl_2mom_init( & open(15,file='nssl_mp_params.out',status='unknown',form='formatted') write(15,NML=nssl_mp_params) close(15) - ENDIF - ENDIF - ENDIF - + ENDIF + ENDIF + ENDIF IF ( iufccn > 0 ) THEN ! make sure to use option that uses UF ccn irenuc = 7 diff --git a/physics/MP/NSSL/mp_nssl.F90 b/physics/MP/NSSL/mp_nssl.F90 index db2edfb86..74e6c780f 100644 --- a/physics/MP/NSSL/mp_nssl.F90 +++ b/physics/MP/NSSL/mp_nssl.F90 @@ -28,8 +28,8 @@ module mp_nssl !! \htmlinclude mp_nssl_init.html !! subroutine mp_nssl_init(ncol, nlev, errflg, errmsg, threads, restart, & - mpirank, mpiroot,mpicomm, & - qc, qr, qi, qs, qh, & + fn_nml, input_nml_file, mpirank, mpiroot, & + mpicomm, qc, qr, qi, qs, qh, & ccw, crw, cci, csw, chw, vh, & con_g, con_rd, con_cp, con_rv, & con_t0c, con_cliq, con_csol, con_eps, & @@ -52,6 +52,8 @@ subroutine mp_nssl_init(ncol, nlev, errflg, errmsg, threads, restart, & integer, intent( out) :: errflg integer, intent(in) :: threads logical, intent(in) :: restart + character(len=*), intent(in) :: fn_nml + character(len=*), intent(in) :: input_nml_file(:) real(kind_phys), intent(in) :: con_g, con_rd, con_cp, con_rv, & con_t0c, con_cliq, con_csol, con_eps @@ -156,7 +158,9 @@ subroutine mp_nssl_init(ncol, nlev, errflg, errmsg, threads, restart, & ! write(0,*) 'call nssl_2mom_init' CALL nssl_2mom_init(ims,ime, jms,jme, kms,kme,nssl_params,ipctmp=ipc,mixphase=0, & - ihvol=ihailv,nssl_ehw0=nssl_ehw0,nssl_ehlw0=nssl_ehlw0,errmsg=errmsg, & + namelist_filename=fn_nml,internal_nml=input_nml_file, & + ihvol=ihailv,nssl_ehw0=nssl_ehw0, & + nssl_ehlw0=nssl_ehlw0,errmsg=errmsg, & nssl_alphar=nssl_alphar, & nssl_alphah=nssl_alphah, & nssl_alphahl=nssl_alphahl, & diff --git a/physics/MP/NSSL/mp_nssl.meta b/physics/MP/NSSL/mp_nssl.meta index e818e8fe5..93a5aa65b 100644 --- a/physics/MP/NSSL/mp_nssl.meta +++ b/physics/MP/NSSL/mp_nssl.meta @@ -49,6 +49,22 @@ dimensions = () type = logical intent = in +[fn_nml] + standard_name = filename_of_namelist + long_name = namelist filename + units = none + dimensions = () + type = character + kind = len=* + intent = in +[input_nml_file] + standard_name = filename_of_internal_namelist + long_name = character string to store full namelist contents + units = none + dimensions = (number_of_lines_in_internal_namelist) + type = character + kind = len=* + intent = in [mpirank] standard_name = mpi_rank long_name = current MPI-rank diff --git a/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 b/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 index 9e3d14cd7..509c50c45 100644 --- a/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 +++ b/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 @@ -2628,7 +2628,7 @@ SUBROUTINE mym_turbulence ( & & sh, sm, & & El, & & Dfm, Dfh, Dfq, Tcd, Qcd, Pdk, Pdt, Pdq, Pdc, & - & qWT1D,qSHEAR1D,qBUOY1D,qDISS1D, & + & qWT1D,qSHEAR1D,qBUOY1D,qDISS1D, & & tke_budget, & & Psig_bl,Psig_shcu,cldfra_bl1D, & & bl_mynn_mixlength, & diff --git a/physics/Radiation/mersenne_twister.f b/physics/Radiation/mersenne_twister.f index bdd0da080..f9f9668f1 100644 --- a/physics/Radiation/mersenne_twister.f +++ b/physics/Radiation/mersenne_twister.f @@ -182,7 +182,7 @@ module mersenne_twister integer,parameter:: tmaskc=-272236544 !< tempering parameter integer,parameter:: mag01(0:1)=(/0,mata/) integer,parameter:: iseed=4357 - integer,parameter:: nrest=n+6 + integer,parameter:: nrest=n+4 ! Defined types type random_stat !< Generator state private diff --git a/physics/Radiation/radiation_aerosols.f b/physics/Radiation/radiation_aerosols.f index 983d592c1..cd5dd49f1 100644 --- a/physics/Radiation/radiation_aerosols.f +++ b/physics/Radiation/radiation_aerosols.f @@ -2901,7 +2901,7 @@ subroutine aer_property & ! --- map grid in longitude direction, lon from 0 to 355 deg resolution ! print *,' Seeking lon index for point i =',i - i3 = i1 + i3 = 1 lab_do_IMXAE : do while ( i3 <= IMXAE ) tmp1 = dltg * (i3 - 1) dtmp = alon(i) - tmp1 @@ -2942,7 +2942,7 @@ subroutine aer_property & ! --- map grid in latitude direction, lat from 90n to 90s in 5 deg resolution ! print *,' Seeking lat index for point i =',i - j3 = j1 + j3 = 1 lab_do_JMXAE : do while ( j3 <= JMXAE ) tmp2 = 90.0 - dltg * (j3 - 1) dtmp = tmp2 - alat(i) diff --git a/physics/SFC_Models/Lake/CLM/clm_lake.f90 b/physics/SFC_Models/Lake/CLM/clm_lake.f90 index 6dd973c8d..f02ed9249 100644 --- a/physics/SFC_Models/Lake/CLM/clm_lake.f90 +++ b/physics/SFC_Models/Lake/CLM/clm_lake.f90 @@ -607,7 +607,7 @@ SUBROUTINE clm_lake_run( & enddo do k = -nlevsnow+1,nlevsoil t_soisno(c,k) = t_soisno3d(i,k) - h2osoi_ice(c,k) = h2osoi_ice3d(i,k) + h2osoi_ice(c,k) = h2osoi_ice3d(i,k) h2osoi_liq(c,k) = h2osoi_liq3d(i,k) h2osoi_vol(c,k) = h2osoi_vol3d(i,k) z(c,k) = z3d(i,k) @@ -678,20 +678,20 @@ SUBROUTINE clm_lake_run( & savedtke12d(i) = savedtke1(c) snowdp2d(i) = snowdp(c) h2osno2d(i) = h2osno(c) - snl2d(i) = snl(c) + snl2d(i) = snl(c) t_grnd2d(i) = t_grnd(c) do k = 1,nlevlake t_lake3d(i,k) = t_lake(c,k) - lake_icefrac3d(i,k) = lake_icefrac(c,k) + lake_icefrac3d(i,k) = lake_icefrac(c,k) enddo - do k = -nlevsnow+1,nlevsoil - z3d(i,k) = z(c,k) - dz3d(i,k) = dz(c,k) - t_soisno3d(i,k) = t_soisno(c,k) - h2osoi_liq3d(i,k) = h2osoi_liq(c,k) - h2osoi_ice3d(i,k) = h2osoi_ice(c,k) + do k = -nlevsnow+1,nlevsoil + z3d(i,k) = z(c,k) + dz3d(i,k) = dz(c,k) + t_soisno3d(i,k) = t_soisno(c,k) + h2osoi_liq3d(i,k) = h2osoi_liq(c,k) + h2osoi_ice3d(i,k) = h2osoi_ice(c,k) h2osoi_vol3d(i,k) = h2osoi_vol(c,k) - enddo + enddo do k = -nlevsnow+0,nlevsoil zi3d(i,k) = zi(c,k) enddo @@ -2305,7 +2305,7 @@ SUBROUTINE ShalLakeTemperature(t_grnd,h2osno,sabg,dz,dz_lake,z,zi, & ! ! unlike eflx_gnet if(abs(errsoi(c)) > .001_kind_lake) then ! 1.e-5_kind_lake) then WRITE( message,* )'Primary soil energy conservation error in shlake & - column during Tridiagonal Solution,', 'error (W/m^2):', c, errsoi(c) + &column during Tridiagonal Solution,', 'error (W/m^2):', c, errsoi(c) errmsg=trim(message) errflg=1 return @@ -5626,7 +5626,7 @@ SUBROUTINE lakeini(kdt, ISLTYP, gt0, snowd, ! initial t_soisno3d ! in snow if(snowdp2d(i) > 0.) then - do k = snl2d(i)+1, 0 + do k = nint(snl2d(i))+1, 0 t_soisno3d(i,k) =min(tfrz,tsfc(i)) enddo endif diff --git a/physics/SFC_Models/Land/Noah/set_soilveg.f b/physics/SFC_Models/Land/Noah/set_soilveg.f index 35f4ace37..8f9c4e782 100644 --- a/physics/SFC_Models/Land/Noah/set_soilveg.f +++ b/physics/SFC_Models/Land/Noah/set_soilveg.f @@ -52,35 +52,35 @@ subroutine set_soilveg(me,isot,ivet,nlunit,errmsg,errflg) !using umd veg table slope_data =(/0.1, 0.6, 1.0, 0.35, 0.55, 0.8, - & 0.63, 0.0, 0.0, 0.0, 0.0, 0.0, - & 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, - & 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, - & 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0/) + & 0.63, 0.0, 0.0, 0.0, 0.0, 0.0, + & 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, + & 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, + & 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0/) rsmtbl =(/300.0, 175.0, 175.0, 300.0, 300.0, 70.0, - & 20.0, 225.0, 225.0, 225.0, 400.0, 20.0, - & 150.0, 0.0, 0.0, 0.0, 0.0, 0.0, - & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) + & 20.0, 225.0, 225.0, 225.0, 400.0, 20.0, + & 150.0, 0.0, 0.0, 0.0, 0.0, 0.0, + & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) c----------------------------- rgltbl =(/30.0, 30.0, 30.0, 30.0, 30.0, 65.0, - & 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, - & 100.0, 0.0, 0.0, 0.0, 0.0, 0.0, - & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) + & 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, + & 100.0, 0.0, 0.0, 0.0, 0.0, 0.0, + & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) hstbl =(/41.69, 54.53, 51.93, 47.35, 47.35, 54.53, - & 36.35, 42.00, 42.00, 42.00, 42.00, 36.35, - & 42.00, 0.00, 0.00, 0.00, 0.00, 0.00, - & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, - & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/) + & 36.35, 42.00, 42.00, 42.00, 42.00, 36.35, + & 42.00, 0.00, 0.00, 0.00, 0.00, 0.00, + & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, + & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/) ! changed for version 2.6 on june 2nd 2003 ! data snupx /0.080, 0.080, 0.080, 0.080, 0.080, 0.080, ! & 0.040, 0.040, 0.040, 0.040, 0.025, 0.040, ! & 0.025, 0.000, 0.000, 0.000, 0.000, 0.000, snupx =(/0.040, 0.040, 0.040, 0.040, 0.040, 0.040, - * 0.020, 0.020, 0.020, 0.020, 0.013, 0.020, - * 0.013, 0.000, 0.000, 0.000, 0.000, 0.000, - & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, - & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/) + * 0.020, 0.020, 0.020, 0.020, 0.013, 0.020, + * 0.013, 0.000, 0.000, 0.000, 0.000, 0.000, + & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, + & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/) bare =11 diff --git a/physics/SFC_Models/Land/Noah/sflx.f b/physics/SFC_Models/Land/Noah/sflx.f index 5f0c6c747..c6822c0eb 100644 --- a/physics/SFC_Models/Land/Noah/sflx.f +++ b/physics/SFC_Models/Land/Noah/sflx.f @@ -2662,7 +2662,7 @@ subroutine snopac ! t1 = tfreez * sncovr**snoexp + t12 * (1.0 - sncovr**snoexp) t1 = tfreez * max(0.01,sncovr**snoexp) + & - & t12 * (1.0 - max(0.01,sncovr**snoexp)) + & t12 * (1.0 - max(0.01,sncovr**snoexp)) beta = 1.0 ssoil = df1 * (t1 - stc(1)) / dtot diff --git a/physics/SFC_Models/Land/Noahmp/module_sf_noahmp_glacier.F90 b/physics/SFC_Models/Land/Noahmp/module_sf_noahmp_glacier.F90 index 5b7aa584d..912516cf5 100644 --- a/physics/SFC_Models/Land/Noahmp/module_sf_noahmp_glacier.F90 +++ b/physics/SFC_Models/Land/Noahmp/module_sf_noahmp_glacier.F90 @@ -1,4 +1,6 @@ +#ifndef CCPP #define CCPP +#endif !> \file module_sf_noahmp_glacier.F90 !! This file contains the NoahMP Glacier scheme. @@ -326,7 +328,7 @@ subroutine noahmp_glacier (& isnow ,snowh ,sneqv ,snice ,snliq ,stc , & !inout dzsnso ,sh2o ,sice ,ponding ,zsnso ,fsh , & !inout runsrf ,runsub ,qsnow ,ponding1 ,ponding2 ,qsnbot , & !out - fpice ,esnow) !out + fpice ,esnow) !out if(opt_gla == 2) then edir = qvap - qdew @@ -638,7 +640,7 @@ subroutine energy_glacier (nsnow ,nsoil ,isnow ,dt ,qsnow ,rhoair call tsnosoi_glacier (nsoil ,nsnow ,isnow ,dt ,tbot , & !in ssoil ,snowh ,zbot ,zsnso ,df , & !in - hcpct , & !in + hcpct , & !in stc ) !inout ! adjusting snow surface temperature @@ -1340,11 +1342,11 @@ subroutine glacier_flux (nsoil ,nsnow ,emg ,isnow ,df ,dzsnso end if csh = rhoair*cpair/rahb - if(snowh > 0.0 .or. opt_gla == 1) then + if(snowh > 0.0 .or. opt_gla == 1) then cev = rhoair*cpair/gamma/(rsurf+rawb) - else - cev = 0.0 ! don't allow any sublimation of glacier in opt_gla=2 - end if + else + cev = 0.0 ! don't allow any sublimation of glacier in opt_gla=2 + end if ! surface fluxes and dtg @@ -1730,7 +1732,7 @@ end subroutine sfcdif1_glacier !>\ingroup NoahMP_LSM subroutine tsnosoi_glacier (nsoil ,nsnow ,isnow ,dt ,tbot , & !in ssoil ,snowh ,zbot ,zsnso ,df , & !in - hcpct , & !in + hcpct , & !in stc ) !inout ! -------------------------------------------------------------------------------------------------- !> compute snow (up to 3l) and soil (4l) temperature. note that snow temperatures @@ -2222,11 +2224,11 @@ subroutine phasechange_glacier (nsnow ,nsoil ,isnow ,dt ,fact , & if (heatr(1) > 0.) then xm(1) = heatr(1)*dt/hfus hm(1) = heatr(1) - imelt(1) = 1 + imelt(1) = 1 else xm(1) = 0. hm(1) = 0. - imelt(1) = 0 + imelt(1) = 0 endif qmelt = max(0.,(temp1-sneqv))/dt xmf = hfus*qmelt @@ -2273,21 +2275,21 @@ subroutine phasechange_glacier (nsnow ,nsoil ,isnow ,dt ,fact , & if (any(stc(1:4) > tfrz) .and. any(stc(1:4) < tfrz)) then do j = 1,nsoil if ( stc(j) > tfrz ) then - heatr(j) = (stc(j)-tfrz)/fact(j) + heatr(j) = (stc(j)-tfrz)/fact(j) do k = 1,nsoil - if (j .ne. k .and. stc(k) < tfrz .and. heatr(j) > 0.1) then - heatr(k) = (stc(k)-tfrz)/fact(k) - if (abs(heatr(k)) > heatr(j)) then ! layer absorbs all - heatr(k) = heatr(k) + heatr(j) - stc(k) = tfrz + heatr(k)*fact(k) - heatr(j) = 0.0 + if (j .ne. k .and. stc(k) < tfrz .and. heatr(j) > 0.1) then + heatr(k) = (stc(k)-tfrz)/fact(k) + if (abs(heatr(k)) > heatr(j)) then ! layer absorbs all + heatr(k) = heatr(k) + heatr(j) + stc(k) = tfrz + heatr(k)*fact(k) + heatr(j) = 0.0 else - heatr(j) = heatr(j) + heatr(k) - heatr(k) = 0.0 - stc(k) = tfrz + heatr(j) = heatr(j) + heatr(k) + heatr(k) = 0.0 + stc(k) = tfrz end if - end if - end do + end if + end do stc(j) = tfrz + heatr(j)*fact(j) end if end do @@ -2298,21 +2300,21 @@ subroutine phasechange_glacier (nsnow ,nsoil ,isnow ,dt ,fact , & if (any(stc(1:4) > tfrz) .and. any(stc(1:4) < tfrz)) then do j = 1,nsoil if ( stc(j) < tfrz ) then - heatr(j) = (stc(j)-tfrz)/fact(j) + heatr(j) = (stc(j)-tfrz)/fact(j) do k = 1,nsoil - if (j .ne. k .and. stc(k) > tfrz .and. heatr(j) < -0.1) then - heatr(k) = (stc(k)-tfrz)/fact(k) - if (heatr(k) > abs(heatr(j))) then ! layer absorbs all - heatr(k) = heatr(k) + heatr(j) - stc(k) = tfrz + heatr(k)*fact(k) - heatr(j) = 0.0 + if (j .ne. k .and. stc(k) > tfrz .and. heatr(j) < -0.1) then + heatr(k) = (stc(k)-tfrz)/fact(k) + if (heatr(k) > abs(heatr(j))) then ! layer absorbs all + heatr(k) = heatr(k) + heatr(j) + stc(k) = tfrz + heatr(k)*fact(k) + heatr(j) = 0.0 else - heatr(j) = heatr(j) + heatr(k) - heatr(k) = 0.0 - stc(k) = tfrz + heatr(j) = heatr(j) + heatr(k) + heatr(k) = 0.0 + stc(k) = tfrz end if - end if - end do + end if + end do stc(j) = tfrz + heatr(j)*fact(j) end if end do @@ -2323,25 +2325,25 @@ subroutine phasechange_glacier (nsnow ,nsoil ,isnow ,dt ,fact , & if (any(stc(1:4) > tfrz) .and. any(mice(1:4) > 0.)) then do j = 1,nsoil if ( stc(j) > tfrz ) then - heatr(j) = (stc(j)-tfrz)/fact(j) + heatr(j) = (stc(j)-tfrz)/fact(j) xm(j) = heatr(j)*dt/hfus do k = 1,nsoil - if (j .ne. k .and. mice(k) > 0. .and. xm(j) > 0.1) then - if (mice(k) > xm(j)) then ! layer absorbs all - mice(k) = mice(k) - xm(j) - xmf = xmf + hfus * xm(j)/dt - stc(k) = tfrz - xm(j) = 0.0 + if (j .ne. k .and. mice(k) > 0. .and. xm(j) > 0.1) then + if (mice(k) > xm(j)) then ! layer absorbs all + mice(k) = mice(k) - xm(j) + xmf = xmf + hfus * xm(j)/dt + stc(k) = tfrz + xm(j) = 0.0 else - xm(j) = xm(j) - mice(k) - xmf = xmf + hfus * mice(k)/dt - mice(k) = 0.0 - stc(k) = tfrz + xm(j) = xm(j) - mice(k) + xmf = xmf + hfus * mice(k)/dt + mice(k) = 0.0 + stc(k) = tfrz end if mliq(k) = max(0.,wmass0(k)-mice(k)) - end if - end do - heatr(j) = xm(j)*hfus/dt + end if + end do + heatr(j) = xm(j)*hfus/dt stc(j) = tfrz + heatr(j)*fact(j) end if end do @@ -2352,25 +2354,25 @@ subroutine phasechange_glacier (nsnow ,nsoil ,isnow ,dt ,fact , & if (any(stc(1:4) < tfrz) .and. any(mliq(1:4) > 0.)) then do j = 1,nsoil if ( stc(j) < tfrz ) then - heatr(j) = (stc(j)-tfrz)/fact(j) + heatr(j) = (stc(j)-tfrz)/fact(j) xm(j) = heatr(j)*dt/hfus do k = 1,nsoil - if (j .ne. k .and. mliq(k) > 0. .and. xm(j) < -0.1) then - if (mliq(k) > abs(xm(j))) then ! layer absorbs all - mice(k) = mice(k) - xm(j) - xmf = xmf + hfus * xm(j)/dt - stc(k) = tfrz - xm(j) = 0.0 + if (j .ne. k .and. mliq(k) > 0. .and. xm(j) < -0.1) then + if (mliq(k) > abs(xm(j))) then ! layer absorbs all + mice(k) = mice(k) - xm(j) + xmf = xmf + hfus * xm(j)/dt + stc(k) = tfrz + xm(j) = 0.0 else - xm(j) = xm(j) + mliq(k) - xmf = xmf - hfus * mliq(k)/dt - mice(k) = wmass0(k) - stc(k) = tfrz + xm(j) = xm(j) + mliq(k) + xmf = xmf - hfus * mliq(k)/dt + mice(k) = wmass0(k) + stc(k) = tfrz end if mliq(k) = max(0.,wmass0(k)-mice(k)) - end if - end do - heatr(j) = xm(j)*hfus/dt + end if + end do + heatr(j) = xm(j)*hfus/dt stc(j) = tfrz + heatr(j)*fact(j) end if end do @@ -2402,7 +2404,7 @@ subroutine water_glacier (nsnow ,nsoil ,imelt ,dt ,prcp ,sfctmp , isnow ,snowh ,sneqv ,snice ,snliq ,stc , & !inout dzsnso ,sh2o ,sice ,ponding ,zsnso ,fsh , & !inout runsrf ,runsub ,qsnow ,ponding1 ,ponding2 ,qsnbot , & !out - fpice ,esnow) !out + fpice ,esnow) !out ! ---------------------------------------------------------------------- ! code history: ! initial code: guo-yue niu, oct. 2007 @@ -2573,7 +2575,7 @@ subroutine snowwater_glacier (nsnow ,nsoil ,imelt ,dt ,sfctmp , & !in ficeold ,zsoil , & !in isnow ,snowh ,sneqv ,snice ,snliq , & !inout sh2o ,sice ,stc ,dzsnso ,zsnso , & !inout - fsh , & !inout + fsh , & !inout qsnbot ,snoflow ,ponding1 ,ponding2) !out ! ---------------------------------------------------------------------- implicit none diff --git a/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 b/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 index 3f27636ff..13decd0be 100644 --- a/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 +++ b/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 @@ -1,4 +1,6 @@ +#ifndef CCPP #define CCPP +#endif !> \file module_sf_noahmplsm.F90 !! This file contains the NoahMP land surface model. @@ -424,7 +426,7 @@ subroutine noahmp_sflx (parameters, & sfctmp , sfcprs , psfc , uu , vv , q2, garea1 , & ! in : forcing qc , soldn , lwdn,thsfc_loc, prslkix,prsik1x,prslk1x,& ! in : forcing pblhx , iz0tlnd , itime ,psi_opt ,& - prcpconv, prcpnonc, prcpshcv, prcpsnow, prcpgrpl, prcphail, & ! in : forcing + prcpconv, prcpnonc, prcpshcv, prcpsnow, prcpgrpl, prcphail, & ! in : forcing tbot , co2air , o2air , foln , ficeold , zlvl , & ! in : forcing ep_1 , ep_2 , epsm1 , cp , & ! in : constants albold , sneqvo , & ! in/out : @@ -436,7 +438,7 @@ subroutine noahmp_sflx (parameters, & cm , ch , tauss , & ! in/out : grain , gdd , pgs , & ! in/out smcwtd ,deeprech , rech , ustarx , & ! in/out : - z0wrf , z0hwrf , ts , & ! out : + z0wrf , z0hwrf , ts , & ! out : fsa , fsr , fira , fsh , ssoil , fcev , & ! out : fgev , fctr , ecan , etran , edir , trad , & ! out : tgb , tgv , t2mv , t2mb , q2v , q2b , & ! out : @@ -445,9 +447,9 @@ subroutine noahmp_sflx (parameters, & qsnbot , ponding , ponding1, ponding2, rssun , rssha , & ! out : albd , albi , albsnd , albsni , & ! out : bgap , wgap , chv , chb , emissi , & ! out : - shg , shc , shb , evg , evb , ghv , & ! out : - ghb , irg , irc , irb , tr , evc , & ! out : - chleaf , chuc , chv2 , chb2 , fpice , pahv , & + shg , shc , shb , evg , evb , ghv , & ! out : + ghb , irg , irc , irb , tr , evc , & ! out : + chleaf , chuc , chv2 , chb2 , fpice , pahv , & pahg , pahb , pah , esnow , canhs , laisun , & laisha , rb , qsfcveg , qsfcbare & #ifdef CCPP @@ -819,7 +821,7 @@ subroutine noahmp_sflx (parameters, & canliq ,canice ,tv ,sfctmp ,tg , & !in qintr ,qdripr ,qthror ,qints ,qdrips ,qthros , & !out pahv ,pahg ,pahb ,qrain ,qsnow ,snowhin, & !out - fwet ,cmc ) !out + fwet ,cmc ) !out ! compute energy budget (momentum & energy fluxes and phase changes) @@ -833,7 +835,7 @@ subroutine noahmp_sflx (parameters, & qsnow ,dzsnso ,lat ,canliq ,canice ,iloc, jloc , & !in thsfc_loc, prslkix,prsik1x,prslk1x,garea1, & !in pblhx ,iz0tlnd, itime ,psi_opt, ep_1, ep_2, epsm1,cp, & - z0wrf ,z0hwrf , & !out + z0wrf ,z0hwrf , & !out imelt ,snicev ,snliqv ,epore ,t2m ,fsno , & !out sav ,sag ,qmelt ,fsa ,fsr ,taux , & !out tauy ,fira ,fsh ,fcev ,fgev ,fctr , & !out @@ -854,7 +856,7 @@ subroutine noahmp_sflx (parameters, & fsrg ,rssun ,rssha ,albd ,albi ,albsnd,albsni, bgap ,wgap, tgv,tgb,& q1 ,q2v ,q2b ,q2e ,chv ,chb , & !out emissi ,pah ,canhs, & - shg,shc,shb,evg,evb,ghv,ghb,irg,irc,irb,tr,evc,chleaf,chuc,chv2,chb2 ) !out + shg,shc,shb,evg,evb,ghv,ghb,irg,irc,irb,tr,evc,chleaf,chuc,chv2,chb2 ) !out qsfcveg = eah*ep_2/(sfcprs + epsm1*eah) qsfcbare = qsfc @@ -877,7 +879,7 @@ subroutine noahmp_sflx (parameters, & esai ,sfctmp ,qvap ,qdew ,zsoil ,btrani , & !in ficeold,ponding,tg ,ist ,fveg ,iloc,jloc , smceq , & !in bdfall ,fp ,rain ,snow , & !in mb/an: v3.7 - qsnow ,qrain ,snowhin,latheav,latheag,frozen_canopy,frozen_ground, & !in mb + qsnow ,qrain ,snowhin,latheav,latheag,frozen_canopy,frozen_ground, & !in mb isnow ,canliq ,canice ,tv ,snowh ,sneqv , & !inout snice ,snliq ,stc ,zsnso ,sh2o ,smc , & !inout sice ,zwt ,wa ,wt ,dzsnso ,wslake , & !inout @@ -911,9 +913,9 @@ subroutine noahmp_sflx (parameters, & if (opt_crop == 1 .and. crop_active) then call carbon_crop (parameters,nsnow ,nsoil ,vegtyp ,dt ,zsoil ,julian , & !in dzsnso ,stc ,smc ,tv ,psn ,foln ,btran , & !in - soldn ,t2m , & !in + soldn ,t2m , & !in lfmass ,rtmass ,stmass ,wood ,stblcp ,fastcp ,grain , & !inout - lai ,sai ,gdd , & !inout + lai ,sai ,gdd , & !inout gpp ,npp ,nee ,autors ,heters ,totsc ,totlb, pgs ) !out end if @@ -964,7 +966,7 @@ subroutine atm (parameters,ep_2,epsm1,sfcprs ,sfctmp ,q2 , & prcpconv,prcpnonc ,prcpshcv,prcpsnow,prcpgrpl,prcphail , & soldn ,cosz ,thair ,qair , & eair ,rhoair ,qprecc ,qprecl ,solad , solai , & - swdown ,bdfall ,rain ,snow ,fp , fpice ,prcp ) + swdown ,bdfall ,rain ,snow ,fp , fpice ,prcp ) ! -------------------------------------------------------------------------------------------------- ! re-process atmospheric forcing ! ---------------------------------------------------------------------- @@ -1037,7 +1039,7 @@ subroutine atm (parameters,ep_2,epsm1,sfcprs ,sfctmp ,q2 , & if(opt_snf == 4) then qprecc = prcpconv + prcpshcv - qprecl = prcpnonc + qprecl = prcpnonc else qprecc = 0.10 * prcp ! should be from the atmospheric model qprecl = 0.90 * prcp ! should be from the atmospheric model @@ -1090,13 +1092,13 @@ subroutine atm (parameters,ep_2,epsm1,sfcprs ,sfctmp ,q2 , & if(opt_snf == 4 .or. opt_snf == 5) then prcp_frozen = prcpsnow + prcpgrpl + prcphail if(prcpnonc > 0. .and. prcp_frozen > 0.) then - fpice = min(1.0,prcp_frozen/prcpnonc) - fpice = max(0.0,fpice) + fpice = min(1.0,prcp_frozen/prcpnonc) + fpice = max(0.0,fpice) if(opt_snf==4) bdfall = bdfall*(prcpsnow/prcp_frozen) + rho_grpl*(prcpgrpl/prcp_frozen) + & rho_hail*(prcphail/prcp_frozen) if(opt_snf==5) bdfall = parameters%prcpiceden - else - fpice = 0.0 + else + fpice = 0.0 endif endif @@ -1233,8 +1235,8 @@ subroutine precip_heat (parameters,iloc ,jloc ,vegtyp ,dt ,uu ,vv bdfall ,rain ,snow ,fp , & !in canliq ,canice ,tv ,sfctmp ,tg , & !in qintr ,qdripr ,qthror ,qints ,qdrips ,qthros , & !out - pahv ,pahg ,pahb ,qrain ,qsnow ,snowhin, & !out - fwet ,cmc ) !out + pahv ,pahg ,pahb ,qrain ,qsnow ,snowhin, & !out + fwet ,cmc ) !out ! ------------------------ code history ------------------------------ ! michael barlage: oct 2013 - split canwater to calculate precip movement for @@ -1336,10 +1338,10 @@ subroutine precip_heat (parameters,iloc ,jloc ,vegtyp ,dt ,uu ,vv qintr = 0. qdripr = 0. qthror = rain - if(canliq > 0.) then ! for case of canopy getting buried - qdripr = qdripr + canliq/dt - canliq = 0.0 - end if + if(canliq > 0.) then ! for case of canopy getting buried + qdripr = qdripr + canliq/dt + canliq = 0.0 + end if end if ! heat transported by liquid water @@ -1363,7 +1365,7 @@ subroutine precip_heat (parameters,iloc ,jloc ,vegtyp ,dt ,uu ,vv ft = max(0.0,(tv - 270.15) / 1.87e5) fv = sqrt(uu*uu + vv*vv) / 1.56e5 ! mb: changed below to reflect the rain assumption that all precip gets intercepted - icedrip = max(0.,canice) * (fv+ft) !mb: removed /dt + icedrip = max(0.,canice) * (fv+ft) !mb: removed /dt qdrips = (fveg * snow - qints) + icedrip qthros = (1.0-fveg) * snow canice= max(0.,canice + (qints - icedrip)*dt) @@ -1371,10 +1373,10 @@ subroutine precip_heat (parameters,iloc ,jloc ,vegtyp ,dt ,uu ,vv qints = 0. qdrips = 0. qthros = snow - if(canice > 0.) then ! for case of canopy getting buried - qdrips = qdrips + canice/dt - canice = 0.0 - end if + if(canice > 0.) then ! for case of canopy getting buried + qdrips = qdrips + canice/dt + canice = 0.0 + end if endif ! print*, "precip_heat canopy through:",3600.0*(fveg * snow - qints) ! print*, "precip_heat canopy drip:",3600.0*max(0.,canice) * (fv+ft) @@ -1404,13 +1406,13 @@ subroutine precip_heat (parameters,iloc ,jloc ,vegtyp ,dt ,uu ,vv if (fveg > 0.0 .and. fveg < 1.0) then pahg = pahg / fveg ! these will be multiplied by fraction later - pahb = pahb / (1.0-fveg) + pahb = pahb / (1.0-fveg) elseif (fveg <= 0.0) then pahb = pahg + pahb ! for case of canopy getting buried pahg = 0.0 - pahv = 0.0 + pahv = 0.0 elseif (fveg >= 1.0) then - pahb = 0.0 + pahb = 0.0 end if pahv = max(pahv,-20.0) ! put some artificial limits here for stability @@ -1677,7 +1679,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in qsnow ,dzsnso ,lat ,canliq ,canice ,iloc , jloc, & !in thsfc_loc, prslkix,prsik1x,prslk1x,garea1, & !in pblhx , iz0tlnd, itime,psi_opt,ep_1, ep_2, epsm1, cp, & - z0wrf ,z0hwrf , & !out + z0wrf ,z0hwrf , & !out imelt ,snicev ,snliqv ,epore ,t2m ,fsno , & !out sav ,sag ,qmelt ,fsa ,fsr ,taux , & !out tauy ,fira ,fsh ,fcev ,fgev ,fctr , & !out @@ -1697,7 +1699,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in t2mv ,t2mb ,fsrv , & fsrg ,rssun ,rssha ,albd ,albi,albsnd ,albsni,bgap ,wgap,tgv,tgb,& q1 ,q2v ,q2b ,q2e ,chv ,chb, emissi,pah,canhs,& - shg,shc,shb,evg,evb,ghv,ghb,irg,irc,irb,tr,evc,chleaf,chuc,chv2,chb2 ) !out + shg,shc,shb,evg,evb,ghv,ghb,irg,irc,irb,tr,evc,chleaf,chuc,chv2,chb2 ) !out !jref:end ! -------------------------------------------------------------------------------------------------- @@ -2211,19 +2213,19 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in if (tv .gt. tfrz) then ! barlage: add distinction between ground and latheav = hvap ! vegetation in v3.6 - frozen_canopy = .false. + frozen_canopy = .false. else latheav = hsub - frozen_canopy = .true. + frozen_canopy = .true. end if gammav = cpair*sfcprs/(ep_2*latheav) if (tg .gt. tfrz) then latheag = hvap - frozen_ground = .false. + frozen_ground = .false. else latheag = hsub - frozen_ground = .true. + frozen_ground = .true. end if gammag = cpair*sfcprs/(ep_2*latheag) @@ -2334,7 +2336,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in ssoil = fveg * ghv + (1.0 - fveg) * ghb fcev = evc fctr = tr - pah = fveg * pahg + (1.0 - fveg) * pahb + pahv + pah = fveg * pahg + (1.0 - fveg) * pahb + pahv tg = fveg * tgv + (1.0 - fveg) * tgb t2m = fveg * t2mv + (1.0 - fveg) * t2mb ts = fveg * tah + (1.0 - fveg) * tgb @@ -2364,7 +2366,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in t2m = t2mb fcev = 0. fctr = 0. - pah = pahb + pah = pahb ts = tg cm = cmb ch = chb @@ -3551,7 +3553,7 @@ subroutine twostream (parameters,ib ,ic ,vegtyp ,cosz ,vai , & ! kopen = 1.0 else if(opt_rad == 1) then - denfveg = -log(max(1.0-fveg,0.01))/(pai*parameters%rc**2) + denfveg = -log(max(1.0-fveg,0.01))/(pai*parameters%rc**2) hd = parameters%hvt - parameters%hvb bb = 0.5 * hd thetap = atan(bb/parameters%rc * tan(acos(max(0.01,cosz))) ) @@ -4263,11 +4265,11 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & shc = fveg*rhoair*cpair*cvh * ( tv-tah) evc = fveg*rhoair*cpair*cew * (estv-eah) / gammav ! barlage: change to v in v3.6 tr = fveg*rhoair*cpair*ctw * (estv-eah) / gammav - if (tv > tfrz) then + if (tv > tfrz) then evc = min(canliq*latheav/dt,evc) ! barlage: add if block for canice in v3.6 - else + else evc = min(canice*latheav/dt,evc) - end if + end if ! canopy heat capacity hcv = fveg*(parameters%cbiom*vaie*cwat + canliq*cwat/denh2o + canice*cice/denice) !j/m2/k @@ -7034,7 +7036,7 @@ subroutine water (parameters,vegtyp ,nsnow ,nsoil ,imelt ,dt ,uu , & esai ,sfctmp ,qvap ,qdew ,zsoil ,btrani , & !in ficeold,ponding,tg ,ist ,fveg ,iloc ,jloc ,smceq , & !in bdfall ,fp ,rain ,snow, & !in mb/an: v3.7 - qsnow ,qrain ,snowhin,latheav,latheag,frozen_canopy,frozen_ground, & !in mb + qsnow ,qrain ,snowhin,latheav,latheag,frozen_canopy,frozen_ground, & !in mb isnow ,canliq ,canice ,tv ,snowh ,sneqv , & !inout snice ,snliq ,stc ,zsnso ,sh2o ,smc , & !inout sice ,zwt ,wa ,wt ,dzsnso ,wslake , & !inout @@ -7663,19 +7665,19 @@ subroutine combine (parameters,nsnow ,nsoil ,iloc ,jloc , & !in snice(j-1) = snice(j-1) + snice(j) dzsnso(j-1) = dzsnso(j-1) + dzsnso(j) else - if(snice(j) >= 0.) then + if(snice(j) >= 0.) then ponding1 = snliq(j) ! isnow will get set to zero below; ponding1 will get sneqv = snice(j) ! added to ponding from phasechange ponding should be snowh = dzsnso(j) ! zero here because it was calculated for thin snow - else ! snice over-sublimated earlier - ponding1 = snliq(j) + snice(j) - if(ponding1 < 0.) then ! if snice and snliq sublimates remove from soil - sice(1) = max(0.0,sice(1)+ponding1/(dzsnso(1)*1000.)) + else ! snice over-sublimated earlier + ponding1 = snliq(j) + snice(j) + if(ponding1 < 0.) then ! if snice and snliq sublimates remove from soil + sice(1) = max(0.0,sice(1)+ponding1/(dzsnso(1)*1000.)) ponding1 = 0.0 - end if + end if sneqv = 0.0 snowh = 0.0 - end if + end if snliq(j) = 0.0 snice(j) = 0.0 dzsnso(j) = 0.0 @@ -9762,7 +9764,7 @@ subroutine carbon_crop (parameters,nsnow ,nsoil ,vegtyp ,dt ,zsoil ,julia dzsnso ,stc ,smc ,tv ,psn ,foln ,btran , & !in soldn ,t2m , & !in lfmass ,rtmass ,stmass ,wood ,stblcp ,fastcp ,grain , & !inout - xlai ,xsai ,gdd , & !inout + xlai ,xsai ,gdd , & !inout gpp ,npp ,nee ,autors ,heters ,totsc ,totlb, pgs ) !out ! ------------------------------------------------------------------------------------------ ! initial crop version created by xing liu @@ -10441,7 +10443,7 @@ end subroutine psn_crop !! subroutine noahmp_options(idveg , iopt_crs , iopt_btr , iopt_run , iopt_sfc , iopt_frz , & iopt_inf, iopt_rad , iopt_alb , iopt_snf , iopt_tbot, iopt_stc , & - iopt_rsf, iopt_soil, iopt_pedo, iopt_crop, iopt_trs , iopt_diag, & + iopt_rsf, iopt_soil, iopt_pedo, iopt_crop, iopt_trs , iopt_diag, & iopt_z0m ) implicit none diff --git a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 index a36b7a60a..baa0171de 100644 --- a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 +++ b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 @@ -1,4 +1,6 @@ +#ifndef CCPP #define CCPP +#endif !> \file noahmpdrv.F90 !! This file contains the NoahMP land surface scheme driver. @@ -1152,7 +1154,7 @@ subroutine noahmpdrv_run & call noahmp_options(idveg ,iopt_crs, iopt_btr , iopt_run, iopt_sfc, & iopt_frz, iopt_inf , iopt_rad, iopt_alb, & iopt_snf, iopt_tbot, iopt_stc, iopt_rsf, & - iopt_soil,iopt_pedo, iopt_crop,iopt_trs, & + iopt_soil,iopt_pedo, iopt_crop,iopt_trs, & iopt_diag,iopt_z0m) if ( vegetation_category == isice_table ) then @@ -1177,7 +1179,7 @@ subroutine noahmpdrv_run & temperature_soil_bot ,forcing_height ,snow_ice_frac_old ,zsoil , & thsfc_loc ,prslkix ,prsik1x ,prslk1x , & air_pressure_surface ,pblhx ,iz0tlnd ,itime , & - vegetation_frac ,area_grid ,psi_opt , & + vegetation_frac ,area_grid ,psi_opt , & con_fvirt ,con_eps ,con_epsm1 ,con_cp , & snowfall ,snow_water_equiv_old ,snow_albedo_old , & cm_noahmp ,ch_noahmp ,snow_levels ,snow_water_equiv , & @@ -1659,7 +1661,7 @@ subroutine transfer_mp_parameters (vegtype,soiltype,slopetype, & parameters%den = den_table(vegtype) !tree density (no. of trunks per m2) parameters%rc = rc_table(vegtype) !tree crown radius (m) parameters%mfsno = mfsno_table(vegtype) !snowmelt m parameter () - parameters%scffac = scffac_table(vegtype) !snow cover factor + parameters%scffac = scffac_table(vegtype) !snow cover factor parameters%cbiom = cbiom_table(vegtype) !canopy biomass heat capacity parameter (m) parameters%saim = saim_table(vegtype,:) !monthly stem area index, one-sided parameters%laim = laim_table(vegtype,:) !monthly leaf area index, one-sided diff --git a/physics/SFC_Models/Land/RUC/set_soilveg_ruc.F90 b/physics/SFC_Models/Land/RUC/set_soilveg_ruc.F90 index 012c81323..8e2f3f54a 100644 --- a/physics/SFC_Models/Land/RUC/set_soilveg_ruc.F90 +++ b/physics/SFC_Models/Land/RUC/set_soilveg_ruc.F90 @@ -45,36 +45,36 @@ subroutine set_soilveg_ruc(me,isot,ivet,nlunit,errmsg,errflg) if(ivet.eq.2) then ! Using umd veg classification slope_data =(/0.1, 0.6, 1.0, 0.35, 0.55, 0.8, & - & 0.63, 0.0, 0.0, 0.0, 0.0, 0.0, & - & 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, & - & 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, & - & 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0/) + & 0.63, 0.0, 0.0, 0.0, 0.0, 0.0, & + & 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, & + & 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, & + & 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0/) ! ---------------------------------------------------------------------- ! vegetation class-related arrays ! ---------------------------------------------------------------------- rstbl =(/300.0, 175.0, 175.0, 300.0, 300.0, 70.0, & - & 20.0, 225.0, 225.0, 225.0, 400.0, 20.0, & - & 150.0, 0.0, 0.0, 0.0, 0.0, 0.0, & - & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & - & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) + & 20.0, 225.0, 225.0, 225.0, 400.0, 20.0, & + & 150.0, 0.0, 0.0, 0.0, 0.0, 0.0, & + & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & + & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) rgltbl =(/30.0, 30.0, 30.0, 30.0, 30.0, 65.0, & - & 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, & - & 100.0, 0.0, 0.0, 0.0, 0.0, 0.0, & - & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & - & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) + & 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, & + & 100.0, 0.0, 0.0, 0.0, 0.0, 0.0, & + & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & + & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) hstbl =(/41.69, 54.53, 51.93, 47.35, 47.35, 54.53, & - & 36.35, 42.00, 42.00, 42.00, 42.00, 36.35, & - & 42.00, 0.00, 0.00, 0.00, 0.00, 0.00, & - & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, & - & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/) + & 36.35, 42.00, 42.00, 42.00, 42.00, 36.35, & + & 42.00, 0.00, 0.00, 0.00, 0.00, 0.00, & + & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, & + & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/) snuptbl =(/0.040, 0.040, 0.040, 0.040, 0.040, 0.040, & & 0.020, 0.020, 0.020, 0.020, 0.013, 0.020, & & 0.013, 0.000, 0.000, 0.000, 0.000, 0.000, & - & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, & - & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/) + & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, & + & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/) bare =11 diff --git a/physics/SFC_Models/SeaIce/CICE/sfc_sice.f b/physics/SFC_Models/SeaIce/CICE/sfc_sice.f index a5904d67c..e5f2deae9 100644 --- a/physics/SFC_Models/SeaIce/CICE/sfc_sice.f +++ b/physics/SFC_Models/SeaIce/CICE/sfc_sice.f @@ -137,7 +137,8 @@ subroutine sfc_sice_run & ! ! - Define constant parameters integer, parameter :: kmi = 2 !< 2-layer of ice - real(kind=kind_phys), parameter :: zero = 0.0_kind_phys, one = 1.0_kind_phys + real(kind=kind_phys), parameter :: zero = 0.0_kind_phys + real(kind=kind_phys), parameter :: one = 1.0_kind_phys real(kind=kind_phys), parameter :: himax = 8.0_kind_phys !< maximum ice thickness allowed real(kind=kind_phys), parameter :: himin = 0.1_kind_phys !< minimum ice thickness required real(kind=kind_phys), parameter :: hsmax = 2.0_kind_phys !< maximum snow depth allowed @@ -541,8 +542,8 @@ subroutine ice3lay real (kind=kind_phys), parameter :: dili = di*li real (kind=kind_phys), parameter :: dsli = ds*li real (kind=kind_phys), parameter :: ki4 = ki*4.0_kind_phys - real (kind=kind_phys), parameter :: zero = 0.0_kind_phys, one = 1.0_kind_phys - + real (kind=kind_phys), parameter :: zero = 0.0_kind_phys + real (kind=kind_phys), parameter :: one = 1.0_kind_phys ! --- inputs: integer, intent(in) :: im, kmi, ipr logical :: lprnt diff --git a/physics/docs/ccpp_doxyfile b/physics/docs/ccpp_doxyfile index 7e40361ae..4f62402e6 100644 --- a/physics/docs/ccpp_doxyfile +++ b/physics/docs/ccpp_doxyfile @@ -1066,7 +1066,6 @@ EXCLUDE = ../Radiation/RRTMGP/rte-rrtmgp \ ../MP/Ferrier_Aligo \ ../MP/Zhao_Carr \ ../PBL/MYJ \ - ../MP/GFDL/GFDL_parse_tracers.F90 \ ../PBL/HEDMF \ ../PBL/SHOC \ ../PBL/saYSU \