diff --git a/CMakeLists.txt b/CMakeLists.txt index 4d5c8eae4..793b360b3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -66,6 +66,20 @@ foreach(typedef_module ${TYPEDEFS}) list(APPEND MODULES_F90 ${CMAKE_CURRENT_BINARY_DIR}/${typedef_module}) endforeach() +# This is a hack because I don't know how to get the pre_build to add these +# modules to the install list --PAR +if (PROJECT STREQUAL "CCPP-NEPTUNE") + # added for CCPP 5: ozne_def.mod, h2o_def.mod, module_radiation_surface.mod + list(APPEND MODULES_F90 ${CMAKE_CURRENT_BINARY_DIR}/namelist_soilveg.mod + ${CMAKE_CURRENT_BINARY_DIR}/set_soilveg_mod.mod + ${CMAKE_CURRENT_BINARY_DIR}/physcons.mod + ${CMAKE_CURRENT_BINARY_DIR}/funcphys.mod + ${CMAKE_CURRENT_BINARY_DIR}/ozne_def.mod + ${CMAKE_CURRENT_BINARY_DIR}/h2o_def.mod + ${CMAKE_CURRENT_BINARY_DIR}/module_radiation_surface.mod + ) +endif (PROJECT STREQUAL "CCPP-NEPTUNE") + #------------------------------------------------------------------------------ # Set the sources: physics schemes set(SCHEMES $ENV{CCPP_SCHEMES}) @@ -263,7 +277,7 @@ target_include_directories(ccpp_physics PUBLIC target_link_libraries(ccpp_physics PUBLIC w3nco::w3nco_d NetCDF::NetCDF_Fortran) -if (PROJECT STREQUAL "CCPP-FV3") +if (PROJECT STREQUAL "CCPP-FV3" OR PROJECT STREQUAL "CCPP-NEPTUNE") # Define where to install the library install(TARGETS ccpp_physics EXPORT ccpp_physics-targets @@ -279,4 +293,4 @@ if (PROJECT STREQUAL "CCPP-FV3") # Define where to install the C headers and Fortran modules #install(FILES ${HEADERS_C} DESTINATION include) install(FILES ${MODULES_F90} DESTINATION include) -endif (PROJECT STREQUAL "CCPP-FV3") +endif (PROJECT STREQUAL "CCPP-FV3" OR PROJECT STREQUAL "CCPP-NEPTUNE") diff --git a/physics/GFS_MP_generic.F90 b/physics/GFS_MP_generic.F90 index 6a8d3bfcb..e3d096f81 100644 --- a/physics/GFS_MP_generic.F90 +++ b/physics/GFS_MP_generic.F90 @@ -96,6 +96,7 @@ subroutine GFS_MP_generic_post_run( errmsg, errflg) ! use machine, only: kind_phys + use mod_calpreciptype implicit none diff --git a/physics/GFS_rad_time_vary.fv3.F90 b/physics/GFS_rad_time_vary.fv3.F90 index 8dd070b12..ac2ae7c20 100644 --- a/physics/GFS_rad_time_vary.fv3.F90 +++ b/physics/GFS_rad_time_vary.fv3.F90 @@ -15,7 +15,8 @@ module GFS_rad_time_vary !> \section arg_table_GFS_rad_time_vary_timestep_init Argument Table !! \htmlinclude GFS_rad_time_vary_timestep_init.html !! - subroutine GFS_rad_time_vary_timestep_init ( & + subroutine GFS_rad_time_vary_timestep_init (nthrds, blksz, lrseeds, & + rseeds, & lslwr, lsswr, isubc_lw, isubc_sw, icsdsw, icsdlw, cnx, cny, isc, jsc, & imap, jmap, sec, kdt, imp_physics, imp_physics_zhao_carr, ps_2delt, & ps_1delt, t_2delt, t_1delt, qv_2delt, qv_1delt, t, qv, ps, errmsg, errflg) @@ -28,6 +29,10 @@ subroutine GFS_rad_time_vary_timestep_init ( implicit none ! Interface variables + integer, intent(in) :: nthrds + integer, intent(in) :: blksz(:) + logical, intent(in) :: lrseeds + integer, intent(in) :: rseeds(:,:) integer, intent(in) :: isubc_lw, isubc_sw, cnx, cny, isc, jsc, kdt integer, intent(in) :: imp_physics, imp_physics_zhao_carr logical, intent(in) :: lslwr, lsswr @@ -46,7 +51,7 @@ subroutine GFS_rad_time_vary_timestep_init ( ! Local variables type (random_stat) :: stat - integer :: ix, j, i, nblks, ipseed + integer :: ii, ix, nb, j, i, nblks, ipseed integer :: numrdm(cnx*cny*2) ! Initialize CCPP error handling variables @@ -55,24 +60,50 @@ subroutine GFS_rad_time_vary_timestep_init ( if (lsswr .or. lslwr) then - !--- call to GFS_radupdate_timestep_init is now in GFS_rrtmg_setup_timestep_init + nblks = size(blksz) + + !--- call to GFS_radupdate_run is now in GFS_rrtmg_setup_run + +!$OMP parallel num_threads(nthrds) default(none) & +!$OMP private (nb,ix,ii, i,j) & +!$OMP shared (lrseeds,isubc_lw,isubc_sw,ipsdlim,ipsd0,ipseed) & +!$OMP shared (cnx,cny,sec,numrdm,stat,nblks,isc,jsc) & +!$OMP shared (blksz,icsdsw,icsdlw,jmap,imap) !--- set up random seed index in a reproducible way for entire cubed-sphere face (lat-lon grid) if ((isubc_lw==2) .or. (isubc_sw==2)) then - ipseed = mod(nint(con_100*sqrt(sec)), ipsdlim) + 1 + ipsd0 - call random_setseed (ipseed, stat) - call random_index (ipsdlim, numrdm, stat) - - do ix=1,size(jmap) - j = jmap(ix) - i = imap(ix) - !--- for testing purposes, replace numrdm with '100' - icsdsw(ix) = numrdm(i+isc-1 + (j+jsc-2)*cnx) - icsdlw(ix) = numrdm(i+isc-1 + (j+jsc-2)*cnx + cnx*cny) - enddo - +!NRL If random seeds supplied by NEPTUNE + if(lrseeds) then + ii = 1 + do nb=1,nblks + do ix=1,blksz(nb) + icsdsw(ix) = rseeds(ix,1) + icsdlw(ix) = rseeds(ix,2) + end do + enddo + else +!$OMP single + ipseed = mod(nint(con_100*sqrt(sec)), ipsdlim) + 1 + ipsd0 + call random_setseed (ipseed, stat) + call random_index (ipsdlim, numrdm, stat) +!$OMP end single + +!$OMP do schedule (dynamic,1) + do nb=1,nblks + do ix=1,blksz(nb) + j = jmap(ix) + i = imap(ix) + !--- for testing purposes, replace numrdm with '100' + icsdsw(ix) = numrdm(i+isc-1 + (j+jsc-2)*cnx) + icsdlw(ix) = numrdm(i+isc-1 + (j+jsc-2)*cnx + cnx*cny) + enddo + enddo +!$OMP end do + end if endif ! isubc_lw and isubc_sw +!$OMP end parallel + if (imp_physics == imp_physics_zhao_carr) then if (kdt == 1) then t_2delt = t diff --git a/physics/GFS_rad_time_vary.fv3.meta b/physics/GFS_rad_time_vary.fv3.meta index 8a4938667..77ad8ebd8 100644 --- a/physics/GFS_rad_time_vary.fv3.meta +++ b/physics/GFS_rad_time_vary.fv3.meta @@ -7,6 +7,34 @@ [ccpp-arg-table] name = GFS_rad_time_vary_timestep_init type = scheme +[nthrds] + standard_name = omp_threads + long_name = number of OpenMP threads available for physics schemes + units = count + dimensions = () + type = integer + intent = in +[blksz] + standard_name = ccpp_block_sizes + long_name = for explicit data blocking: block sizes of all blocks + units = count + dimensions = (ccpp_block_count) + type = integer + intent = in +[lrseeds] + standard_name = flag_for_model_specific_random_seeds + long_name = flag for model specific random seeds + units = flag + dimensions = () + type = logical + intent = in +[rseeds] + standard_name = model_specific_random_seeds + long_name = model specific random seeds + units = none + dimensions = (horizontal_dimension,number_of_random_streams) + type = integer + intent = in [lslwr] standard_name = flag_for_calling_longwave_radiation long_name = logical flags for lw radiation calls @@ -207,4 +235,3 @@ dimensions = () type = integer intent = out - diff --git a/physics/GFS_rrtmg_post.meta b/physics/GFS_rrtmg_post.meta index 8677d9900..e67c76de4 100644 --- a/physics/GFS_rrtmg_post.meta +++ b/physics/GFS_rrtmg_post.meta @@ -38,7 +38,7 @@ [ltp] standard_name = extra_top_layer long_name = extra top layers - units = none + units = count dimensions = () type = integer intent = in diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index dbea66985..01cbc666f 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -16,7 +16,7 @@ end subroutine GFS_rrtmg_pre_init !! ! Attention - the output arguments lm, im, lmk, lmp must not be set ! in the CCPP version - they are defined in the interstitial_create routine - subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & + subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp, n_var_lndp, & imfdeepcnv, imfdeepcnv_gf, me, ncnd, ntrac, num_p3d, npdf3d, ncnvcld3d,& ntqv, ntcw,ntiw, ntlnc, ntinc, ntrw, ntsw, ntgl, ntwa, ntoz, & ntclamt, nleffr, nieffr, nseffr, lndp_type, kdt, imp_physics, & @@ -37,11 +37,10 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & clouds9, cldsa, cldfra, faersw1, faersw2, faersw3, faerlw1, faerlw2, & faerlw3, alpha, errmsg, errflg) - use machine, only: kind_phys + use machine, only: kind_phys, r8=>kind_dbl_prec use physparam - - use radcons, only: itsfc,ltp, lextop, qmin, & + use radcons, only: itsfc, qmin, & qme5, qme6, epsq, prsmin use funcphys, only: fpvs @@ -79,7 +78,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & implicit none - integer, intent(in) :: im, levs, lm, lmk, lmp, n_var_lndp, & + integer, intent(in) :: im, levs, lm, lmk, lmp, ltp, n_var_lndp, & imfdeepcnv, & imfdeepcnv_gf, me, ncnd, ntrac, & num_p3d, npdf3d, ncnvcld3d, ntqv, & @@ -98,7 +97,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & character(len=3), dimension(:), intent(in) :: lndp_var_list - logical, intent(in) :: lsswr, lslwr, ltaerosol, lgfdlmprad, & + logical, intent(in) :: lextop, lsswr, lslwr, ltaerosol, lgfdlmprad, & uni_cld, effr_in, do_mynnedmf, & lmfshal, lmfdeep2, pert_clds @@ -295,6 +294,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & plyr(i,k1) = prsl(i,k2) * 0.01 ! pa to mb (hpa) tlyr(i,k1) = tgrs(i,k2) prslk1(i,k1) = prslk(i,k2) + rho(i,k1) = prsl(i,k2)/(con_rd*tlyr(i,k1)) + orho(i,k1) = 1.0/rho(i,k1) !> - Compute relative humidity. es = min( prsl(i,k2), fpvs( tgrs(i,k2) ) ) ! fpvs and prsl in pa @@ -351,6 +352,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & plyr(i,lyb) = 0.5 * plvl(i,lla) tlyr(i,lyb) = tlyr(i,lya) prslk1(i,lyb) = (plyr(i,lyb)*0.001) ** rocp ! plyr in hPa + rho(i,lyb) = plyr(i,lyb) *100.0/(con_rd*tlyr(i,lyb)) + orho(i,lyb) = 1.0/rho(i,lyb) rhly(i,lyb) = rhly(i,lya) qstl(i,lyb) = qstl(i,lya) enddo diff --git a/physics/GFS_rrtmg_pre.meta b/physics/GFS_rrtmg_pre.meta index 09ed62f7c..d567afafc 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -44,6 +44,18 @@ dimensions = () type = integer intent = in +[lextop] + standard_name = flag_for_extra_top_layer_for_radiation_calculations + long_name = flag for extra top layer for radiation calculations + units = flag + dimensions = () + type = logical +[ltp] + standard_name = extra_top_layer + long_name = extra top layer + units = count + dimensions = () + type = integer [n_var_lndp] standard_name = number_of_perturbed_land_surface_variables long_name = number of land surface variables perturbed diff --git a/physics/GFS_rrtmg_setup.F90 b/physics/GFS_rrtmg_setup.F90 index 0e2d87feb..6df82672b 100644 --- a/physics/GFS_rrtmg_setup.F90 +++ b/physics/GFS_rrtmg_setup.F90 @@ -11,8 +11,6 @@ module GFS_rrtmg_setup & iswcliq, & & kind_phys - use radcons, only: ltp, lextop - implicit none public GFS_rrtmg_setup_init, GFS_rrtmg_setup_timestep_init, GFS_rrtmg_setup_finalize @@ -46,7 +44,7 @@ subroutine GFS_rrtmg_setup_init ( & si, levr, ictm, isol, ico2, iaer, ntcw, & num_p3d, npdf3d, ntoz, iovr, isubc_sw, isubc_lw, & icliq_sw, crick_proof, ccnorm, & - imp_physics, & + imp_physics, ltp, & norad_precip, idate, iflip, & do_RRTMGP, me, errmsg, errflg) ! ================= subprogram documentation block ================ ! @@ -164,6 +162,7 @@ subroutine GFS_rrtmg_setup_init ( & logical, intent(in) :: crick_proof logical, intent(in) :: ccnorm integer, intent(in) :: imp_physics + integer, intent(in) :: ltp logical, intent(in) :: norad_precip integer, intent(in) :: idate(:) integer, intent(in) :: iflip @@ -241,7 +240,7 @@ subroutine GFS_rrtmg_setup_init ( & call radinit & ! --- inputs: - & ( si, levr, imp_physics, me ) + & ( si, levr, imp_physics, ltp, me ) ! --- outputs: ! ( none ) @@ -322,7 +321,7 @@ end subroutine GFS_rrtmg_setup_finalize ! Private functions - subroutine radinit( si, NLAY, imp_physics, me ) + subroutine radinit( si, NLAY, imp_physics, ltp, me ) !................................... ! --- inputs: @@ -435,7 +434,7 @@ subroutine radinit( si, NLAY, imp_physics, me ) implicit none ! --- inputs: - integer, intent(in) :: NLAY, me, imp_physics + integer, intent(in) :: NLAY, me, imp_physics, ltp real (kind=kind_phys), intent(in) :: si(:) @@ -465,7 +464,7 @@ subroutine radinit( si, NLAY, imp_physics, me ) print *,' IVFLIP=',ivflip,' IOVR=',iovrRad, & & ' ISUBCSW=',isubcsw,' ISUBCLW=',isubclw print *,' LCRICK=',lcrick,' LCNORM=',lcnorm,' LNOPREC=',lnoprec - print *,' LTP =',ltp,', add extra top layer =',lextop + print *,' LTP =',ltp,', add extra top layer =', ltp == 1 if ( ictmflg==0 .or. ictmflg==-2 ) then print *,' Data usage is limited by initial condition!' diff --git a/physics/GFS_rrtmg_setup.meta b/physics/GFS_rrtmg_setup.meta index ecd849c48..3a93a023f 100644 --- a/physics/GFS_rrtmg_setup.meta +++ b/physics/GFS_rrtmg_setup.meta @@ -128,6 +128,13 @@ dimensions = () type = integer intent = in +[ltp] + standard_name = extra_top_layer + long_name = extra top layer for radiation calculations + units = count + dimensions = () + type = integer + intent = in [norad_precip] standard_name = flag_for_turning_off_precipitation_radiative_effect long_name = radiation precip flag for Ferrier/Moorthi diff --git a/physics/calpreciptype.f90 b/physics/calpreciptype.f90 index dcc8ed49b..3eb9eb2d1 100644 --- a/physics/calpreciptype.f90 +++ b/physics/calpreciptype.f90 @@ -1,3 +1,5 @@ + module mod_calpreciptype + contains !>\file calpreciptype.f90 !! This file contains the subroutines that calculates dominant precipitation type. @@ -1377,3 +1379,4 @@ subroutine calwxt_dominant(nalg,rain,freezr,sleet,snow, & return end !! @} + end module mod_calpreciptype diff --git a/physics/cires_orowam2017.f b/physics/cires_orowam2017.f index c20f98f42..c235cda87 100644 --- a/physics/cires_orowam2017.f +++ b/physics/cires_orowam2017.f @@ -1,3 +1,5 @@ + module mod_cires_orowam2017 + contains subroutine oro_wam_2017(im, levs,npt,ipt, kref,kdt,me,master, & dtp,dxres, taub, u1, v1, t1, xn, yn, bn2, rho, prsi, prsL, & del, sigma, hprime, gamma, theta, @@ -384,3 +386,4 @@ subroutine ugwpv0_tofd1d(levs, sigflt, elvmax, zsurf, enddo ! end subroutine ugwpv0_tofd1d + end module mod_cires_orowam2017 diff --git a/physics/cires_ugwp.F90 b/physics/cires_ugwp.F90 index 6efce96f5..f88389ba2 100644 --- a/physics/cires_ugwp.F90 +++ b/physics/cires_ugwp.F90 @@ -19,6 +19,10 @@ module cires_ugwp use gwdps, only: gwdps_run + use mod_cires_ugwp_triggers + + use mod_ugwp_driver_v0 + implicit none private diff --git a/physics/cires_ugwp_triggers.F90 b/physics/cires_ugwp_triggers.F90 index 4a8b97590..966dba575 100644 --- a/physics/cires_ugwp_triggers.F90 +++ b/physics/cires_ugwp_triggers.F90 @@ -1,3 +1,5 @@ + module mod_cires_ugwp_triggers + contains ! subroutine slat_geos5_tamp_v0(im, tau_amp, xlatdeg, tau_gw) !================= @@ -97,3 +99,4 @@ subroutine init_nazdir_v0(naz, xaz, yaz) yaz(4) =-1.0 !S endif end subroutine init_nazdir_v0 + end module mod_cires_ugwp_triggers diff --git a/physics/cires_ugwpv1_sporo.F90 b/physics/cires_ugwpv1_sporo.F90 index c840b49d8..639c9dcb0 100644 --- a/physics/cires_ugwpv1_sporo.F90 +++ b/physics/cires_ugwpv1_sporo.F90 @@ -1,4 +1,5 @@ - + module mod_cires_ugwpv1_sporo + contains subroutine oro_spectral_solver(im, levs,npt,ipt, kref,kdt,me,master, & dtp,dxres, taub, u1, v1, t1, xn, yn, bn2, rho, prsi, prsL, & del, sigma, hprime, gamma, theta, & @@ -348,4 +349,4 @@ subroutine oro_meanflow(nz, nzi, u1, v1, t1, pint, pmid, & dzi(k) = dzi(k-1) end subroutine oro_meanflow - + end module mod_cires_ugwpv1_sporo diff --git a/physics/cu_gf_deep.F90 b/physics/cu_gf_deep.F90 index a07523342..76adc3450 100644 --- a/physics/cu_gf_deep.F90 +++ b/physics/cu_gf_deep.F90 @@ -2623,10 +2623,10 @@ subroutine cup_env(z,qes,he,hes,t,q,p,z1, & p,t,q real(kind=kind_phys), dimension (its:ite,kts:kte) & ,intent (out ) :: & - he,hes,qes + hes,qes real(kind=kind_phys), dimension (its:ite,kts:kte) & ,intent (inout) :: & - z + he,z real(kind=kind_phys), dimension (its:ite) & ,intent (in ) :: & psur,z1 diff --git a/physics/gcycle.F90 b/physics/gcycle.F90 index 5f4f959c6..894d558ff 100644 --- a/physics/gcycle.F90 +++ b/physics/gcycle.F90 @@ -209,6 +209,7 @@ subroutine gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit, ! enddo ! +#if 0 #ifndef INTERNAL_FILE_NML inquire (file=trim(Model%fn_nml),exist=exists) if (.not. exists) then @@ -218,6 +219,7 @@ subroutine gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit, open (unit=Model%nlunit, file=trim(Model%fn_nml), action='READ', status='OLD', iostat=ios) rewind (Model%nlunit) endif +#endif #endif CALL SFCCYCLE (9998, npts, max(lsoil,lsoil_lsm), sig1t, fhcyc, & idate(4), idate(2), idate(3), idate(1), & @@ -232,9 +234,11 @@ subroutine gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit, nlunit, size(input_nml_file), input_nml_file, & min_ice, ialb, isot, ivegsrc, & trim(tile_num_ch), i_indx, j_indx) +#if 0 #ifndef INTERNAL_FILE_NML close (Model%nlunit) #endif +#endif ! if ( nsst > 0 ) then tref = TSFFCS diff --git a/physics/h2ophys.f b/physics/h2ophys.f index e21417e80..ce3792d6a 100644 --- a/physics/h2ophys.f +++ b/physics/h2ophys.f @@ -122,7 +122,7 @@ subroutine h2ophys_run(im, levs, kh2o, dt, h2o, ph2o, prsl, & enddo endif do i=1,im - if (prsl(i,l) < prsmax) then + if (prsl(i,l) < prsmax .and. pltc(i,2).ne.0.) then ! TODO JM 20211104 h2oib(i) = h2o(i,l) ! no filling tem = 1.0 / pltc(i,2) ! 1/teff h2o(i,l) = (h2oib(i) + (pltc(i,1)+pltc(i,3)*tem)*dt) diff --git a/physics/iounitdef.f b/physics/iounitdef.f index 61c711bb1..efffda365 100644 --- a/physics/iounitdef.f +++ b/physics/iounitdef.f @@ -57,7 +57,7 @@ module module_iounitdef ! integer, parameter :: NISIGI2 = 12 integer, parameter :: NISFCI = 14 integer, parameter :: NICO2TR = 15 - integer, parameter :: NICO2CN = 102 + integer, parameter :: NICO2CN = 113 ! CCE (Cray) forbids 100-102 20211112 JM integer, parameter :: NIMTNVR = 24 integer, parameter :: NIDTBTH = 27 integer, parameter :: NIO3PRD = 28 @@ -66,9 +66,9 @@ module module_iounitdef ! integer, parameter :: NICLTUN = 43 integer, parameter :: NIO3CLM = 48 integer, parameter :: NIMICPH = 1 - integer, parameter :: NISFCYC = 101 - integer, parameter :: NIAERCM = 102 - integer, parameter :: NIRADSF = 102 + integer, parameter :: NISFCYC = 111 ! CCE (Cray) forbids 100-102 20210701 JM + integer, parameter :: NIAERCM = 112 ! CCE (Cray) forbids 100-102 20210701 JM + integer, parameter :: NIRADSF = 112 ! CCE (Cray) forbids 100-102 20210701 JM ! --- ... output units diff --git a/physics/mfpbl.f b/physics/mfpbl.f index 2df84945b..286ceb097 100644 --- a/physics/mfpbl.f +++ b/physics/mfpbl.f @@ -1,3 +1,5 @@ + module mod_mfpbl + contains !> \file mfpbl.f !! This file contains the subroutine that calculates the updraft properties and mass flux for use in the Hybrid EDMF PBL scheme. @@ -396,3 +398,4 @@ subroutine mfpbl(im,ix,km,ntrac,delt,cnvflg, & return end !> @} + end module mod_mfpbl diff --git a/physics/mfpblt.f b/physics/mfpblt.f index bd0baf558..0dcef2646 100644 --- a/physics/mfpblt.f +++ b/physics/mfpblt.f @@ -1,3 +1,5 @@ + module mod_mfpblt + contains !>\file mfpblt.f !! This file contains the subroutine that calculates mass flux and !! updraft parcel properties for thermals driven by surface heating @@ -452,3 +454,4 @@ subroutine mfpblt(im,ix,km,kmpbl,ntcw,ntrac1,delt, & return end !> @} + end module mod_mfpblt diff --git a/physics/mfpbltq.f b/physics/mfpbltq.f index b906052cd..d506da00a 100644 --- a/physics/mfpbltq.f +++ b/physics/mfpbltq.f @@ -1,3 +1,5 @@ + module mod_mfpbltq + contains !>\file mfpbltq.f !! This file contains the subroutine that calculates mass flux and !! updraft parcel properties for thermals driven by surface heating @@ -468,3 +470,4 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, return end !> @} + end module mod_mfpbltq diff --git a/physics/mfscu.f b/physics/mfscu.f index 9128c7c10..c258b0c40 100644 --- a/physics/mfscu.f +++ b/physics/mfscu.f @@ -1,3 +1,5 @@ + module mod_mfscu + contains !>\file mfscu.f !! This file contains the mass flux and downdraft parcel preperties !! parameterization for stratocumulus-top-driven turbulence. @@ -554,3 +556,4 @@ subroutine mfscu(im,ix,km,kmscu,ntcw,ntrac1,delt, & return end !> @} + end module mod_mfscu diff --git a/physics/mfscuq.f b/physics/mfscuq.f index 3390c3e58..4c556de82 100644 --- a/physics/mfscuq.f +++ b/physics/mfscuq.f @@ -1,3 +1,5 @@ + module mod_mfscuq + contains !>\file mfscuq.f !! This file contains the mass flux and downdraft parcel preperties !! parameterization for stratocumulus-top-driven turbulence (updated version). @@ -548,3 +550,4 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, return end !> @} + end module mod_mfscuq diff --git a/physics/module_bl_mynn.F90 b/physics/module_bl_mynn.F90 index d691de909..05891fdb6 100644 --- a/physics/module_bl_mynn.F90 +++ b/physics/module_bl_mynn.F90 @@ -950,7 +950,7 @@ SUBROUTINE mym_length ( & alp2 = 0.65 alp3 = 3.0 alp4 = 20. - alp5 = 0.4 + alp5 = 1.2 ! Impose limits on the height integration for elt and the transition layer depth zi2=MAX(zi,minzi) @@ -1045,7 +1045,7 @@ SUBROUTINE mym_length ( & alp2 = 0.30 + 0.3*MIN(MAX((dx - 3000.)/10000., 0.0), 1.0) alp3 = 2.0 alp4 = 20. !10. - alp5 = alp2 !like alp2, but for free atmosphere + alp5 = 2.0*alp2 !like alp2, but for free atmosphere alp6 = 50.0 !used for MF mixing length ! Impose limits on the height integration for elt and the transition layer depth @@ -1056,13 +1056,15 @@ SUBROUTINE mym_length ( & h2=h1*0.5 ! 1/4 transition layer depth qtke(kts)=MAX(0.5*qke(kts),0.01) !tke at full sigma levels + thetaw(kts)=theta(kts) !theta at full-sigma levels qkw(kts) = SQRT(MAX(qke(kts),1.0e-10)) DO k = kts+1,kte afk = dz(k)/( dz(k)+dz(k-1) ) abk = 1.0 -afk qkw(k) = SQRT(MAX(qke(k)*abk+qke(k-1)*afk,1.0e-3)) - qtke(k) = 0.5*qkw(k) ! qkw -> TKE + qtke(k) = 0.5*qkw(k)**2 ! qkw -> TKE + thetaw(k)= theta(k)*abk + theta(k-1)*afk END DO elt = 1.0e-5 @@ -1089,6 +1091,9 @@ SUBROUTINE mym_length ( & el(kts) = 0.0 zwk1 = zw(kts+1) + ! COMPUTE BouLac mixing length + CALL boulac_length(kts,kte,zw,dz,qtke,thetaw,elBLmin,elBLavg) + DO k = kts+1,kte zwk = zw(k) !full-sigma levels cldavg = 0.5*(cldfra_bl1D(k-1)+cldfra_bl1D(k)) @@ -1102,7 +1107,9 @@ SUBROUTINE mym_length ( & & alp6*edmf_a1(k)*edmf_w1(k)) / bv & & *( 1.0 + alp3*SQRT( vsc/( bv*elt ) ) ) elb = MIN(alp5*qkw(k)/bv, zwk) - elf = elb/(1. + (elb/600.)) !bound free-atmos mixing length to < 600 m. + !elf = elb/(1. + (elb/600.)) !bound free-atmos mixing length to < 600 m. + elf = max(600.0,dz(k)) + elf = elf*tanh(elb/elf) !bound free-atmos mixing length to < 600 m. !IF (zwk > zi .AND. elf > 400.) THEN ! ! COMPUTE BouLac mixing length ! !CALL boulac_length0(k,kts,kte,zw,dz,qtke,thetaw,elBLmin0,elBLavg0) @@ -1124,10 +1131,10 @@ SUBROUTINE mym_length ( & tau_cloud = MIN(MAX(0.5*zi/((gtr*zi*MAX(flt,1.0e-4))**onethird),50.),150.) !minimize influence of surface heat flux on tau far away from the PBLH. wt=.5*TANH((zwk - (zi2+h1))/h2) + .5 - tau_cloud = tau_cloud*(1.-wt) + 50.*wt + tau_cloud = tau_cloud*(1.-wt) + 300.*wt elb = MIN(tau_cloud*SQRT(MIN(qtke(k),30.)), zwk) - elf = elb + elf = MIN(MAX(elb,dz(k)),zwk) elb_mf = elb END IF @@ -1148,7 +1155,8 @@ SUBROUTINE mym_length ( & ! "el_unstab" = blended els-elt el_unstab = els/(1. + (els1/elt)) el(k) = MIN(el_unstab, elb_mf) - el(k) = el(k)*(1.-wt) + elf*wt + !el(k) = el(k)*(1.-wt) + elf*wt + el(k) = el(k)*(1.-wt) + alp5*elBLmin(k)*wt ! include scale-awareness. For now, use simple asymptotic kz -> 12 m. el_les= MIN(els/(1. + (els1/12.)), elb_mf) @@ -1476,6 +1484,8 @@ SUBROUTINE boulac_length(kts,kte,zw,dz,qtke,theta,lb1,lb2) !lb2(iz) = 0.5*(dlu(iz)+dld(iz)) !average !Apply soft limit (only impacts very large lb; lb=100 by 5%, lb=500 by 20%). + !lb1(iz) = lb1(iz)*tanh(lb1(iz)/Lmax) + !lb2(iz) = lb2(iz)*tanh(lb2(iz)/Lmax) lb1(iz) = lb1(iz)/(1. + (lb1(iz)/Lmax)) lb2(iz) = lb2(iz)/(1. + (lb2(iz)/Lmax)) @@ -2511,7 +2521,7 @@ SUBROUTINE mym_condensation (kts,kte, & !CLOUD FRACTION. rr2 = 1/SQRT(2) = 0.707 cldfra_bl1D(k) = 0.5*( 1.0+erf( q1(k)*rr2 ) ) - eq1 = rrp*EXP( -0.5*q1k*q1k ) + eq1 = rrp*EXP( -0.5*q1k*q1k ) ! TODO q1k is used before defined, 20211112 JM qll = MAX( cldfra_bl1D(k)*q1k + eq1, 0.0 ) !ESTIMATED LIQUID WATER CONTENT (UNNORMALIZED) ql(k) = alp(k)*sgm(k)*qll diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index dbf554508..4d2cbf5fc 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -1318,7 +1318,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & m = RSHIFT(ABS(rand_perturb_on),1) if (MOD(m,2) .ne. 0) rand2 = rand_pert(i,1,j)*2. m = RSHIFT(ABS(rand_perturb_on),2) - if (MOD(m,2) .ne. 0) rand3 = 0.25*(rand_pert(i,1,j)+ABS(min_rand)) + if (MOD(m,2) .ne. 0) rand3 = 0.25*(rand_pert(i,1,j)+ABS(min_rand)) !TODO min_rand used before set 2021112 JM m = RSHIFT(ABS(rand_perturb_on),3) endif !+---+-----------------------------------------------------------------+ @@ -1505,69 +1505,69 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & jmax_qc = j kmax_qc = k qc_max = qc1d(k) - elseif (qc1d(k) .lt. 0.0) then - write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qc ', qc1d(k), & - ' at i,j,k=', i,j,k + !elseif (qc1d(k) .lt. 0.0) then + ! write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qc ', qc1d(k), & + ! ' at i,j,k=', i,j,k endif if (qr1d(k) .gt. qr_max) then imax_qr = i jmax_qr = j kmax_qr = k qr_max = qr1d(k) - elseif (qr1d(k) .lt. 0.0) then - write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qr ', qr1d(k), & - ' at i,j,k=', i,j,k + !elseif (qr1d(k) .lt. 0.0) then + ! write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qr ', qr1d(k), & + ! ' at i,j,k=', i,j,k endif if (nr1d(k) .gt. nr_max) then imax_nr = i jmax_nr = j kmax_nr = k nr_max = nr1d(k) - elseif (nr1d(k) .lt. 0.0) then - write(*,'(a,e16.7,a,3i8)') 'WARNING, negative nr ', nr1d(k), & - ' at i,j,k=', i,j,k + !elseif (nr1d(k) .lt. 0.0) then + ! write(*,'(a,e16.7,a,3i8)') 'WARNING, negative nr ', nr1d(k), & + ! ' at i,j,k=', i,j,k endif if (qs1d(k) .gt. qs_max) then imax_qs = i jmax_qs = j kmax_qs = k qs_max = qs1d(k) - elseif (qs1d(k) .lt. 0.0) then - write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qs ', qs1d(k), & - ' at i,j,k=', i,j,k + !elseif (qs1d(k) .lt. 0.0) then + ! write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qs ', qs1d(k), & + ! ' at i,j,k=', i,j,k endif if (qi1d(k) .gt. qi_max) then imax_qi = i jmax_qi = j kmax_qi = k qi_max = qi1d(k) - elseif (qi1d(k) .lt. 0.0) then - write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qi ', qi1d(k), & - ' at i,j,k=', i,j,k + !elseif (qi1d(k) .lt. 0.0) then + ! write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qi ', qi1d(k), & + ! ' at i,j,k=', i,j,k endif if (qg1d(k) .gt. qg_max) then imax_qg = i jmax_qg = j kmax_qg = k qg_max = qg1d(k) - elseif (qg1d(k) .lt. 0.0) then - write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qg ', qg1d(k), & - ' at i,j,k=', i,j,k + !elseif (qg1d(k) .lt. 0.0) then + ! write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qg ', qg1d(k), & + ! ' at i,j,k=', i,j,k endif if (ni1d(k) .gt. ni_max) then imax_ni = i jmax_ni = j kmax_ni = k ni_max = ni1d(k) - elseif (ni1d(k) .lt. 0.0) then - write(*,'(a,e16.7,a,3i8)') 'WARNING, negative ni ', ni1d(k), & - ' at i,j,k=', i,j,k + !elseif (ni1d(k) .lt. 0.0) then + ! write(*,'(a,e16.7,a,3i8)') 'WARNING, negative ni ', ni1d(k), & + ! ' at i,j,k=', i,j,k endif if (qv1d(k) .lt. 0.0) then - write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qv ', qv1d(k), & - ' at i,j,k=', i,j,k + ! write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qv ', qv1d(k), & + ! ' at i,j,k=', i,j,k if (k.lt.kte-2 .and. k.gt.kts+1) then - write(*,*) ' below and above are: ', qv(i,k-1,j), qv(i,k+1,j) + ! write(*,*) ' below and above are: ', qv(i,k-1,j), qv(i,k+1,j) qv(i,k,j) = MAX(1.E-7, 0.5*(qv(i,k-1,j) + qv(i,k+1,j))) else qv(i,k,j) = 1.E-7 diff --git a/physics/moninedmf.f b/physics/moninedmf.f index 19e055da4..d55057e1b 100644 --- a/physics/moninedmf.f +++ b/physics/moninedmf.f @@ -5,7 +5,10 @@ !> This module contains the CCPP-compliant hybrid eddy-diffusivity mass-flux !! scheme. module hedmf - + + use mod_mfpbl + use mod_tridi + contains !> \section arg_table_hedmf_init Argument Table diff --git a/physics/physcons.F90 b/physics/physcons.F90 index 41d37491a..dcfb3a26c 100644 --- a/physics/physcons.F90 +++ b/physics/physcons.F90 @@ -35,7 +35,7 @@ !! constants for GCM models. module physcons ! - use machine, only: kind_phys, kind_dyn + use machine ! implicit none ! diff --git a/physics/radiation_astronomy.f b/physics/radiation_astronomy.f index f1651ca84..5244cce94 100644 --- a/physics/radiation_astronomy.f +++ b/physics/radiation_astronomy.f @@ -898,7 +898,9 @@ subroutine coszmn & do i = 1, IM coszdg(i) = coszen(i) * rstp - if (istsun(i) > 0) coszen(i) = coszen(i) / istsun(i) + if (istsun(i) > 0 .and. coszen(i).ne.0) then + coszen(i) = coszen(i) / istsun(i) + endif enddo ! return diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 2f34041c2..10e121acc 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -1616,6 +1616,13 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & ktcon1(i) = k flg(i) = .false. endif +!NRL MNM: Limit overshooting not to be deeper than the actual cloud + tem = zi(i,ktcon(i))-zi(i,kbcon(i)) + tem1 = zi(i,ktcon1(i))-zi(i,ktcon(i)) + if(tem1.ge.tem) then + ktcon1(i) = k + flg(i) = .false. + endif endif endif enddo @@ -2950,6 +2957,9 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp t1(i,k) = t1(i,k) + tem2 * dellat q1(i,k) = q1(i,k) + tem2 * dellaq(i,k) +! NRL MNM: Limit q1 + val = epsilon(1.0_kind_phys) + q1(i,k) = max(q1(i,k), val ) ! tem = tem2 / rcs(i) ! u1(i,k) = u1(i,k) + dellau(i,k) * tem ! v1(i,k) = v1(i,k) + dellav(i,k) * tem diff --git a/physics/sascnvn.F b/physics/sascnvn.F index f9d91578b..8dc118ad6 100644 --- a/physics/sascnvn.F +++ b/physics/sascnvn.F @@ -1008,6 +1008,13 @@ subroutine sascnvn_run( ktcon1(i) = k flg(i) = .false. endif +!NRL MNM: Limit overshooting not to be deeper than the actual cloud + tem = zi(i,ktcon(i))-zi(i,kbcon(i)) + tem1 = zi(i,ktcon1(i))-zi(i,ktcon(i)) + if(tem1.ge.tem) then + ktcon1(i) = k + flg(i) = .false. + endif endif endif enddo @@ -1874,6 +1881,9 @@ subroutine sascnvn_run( dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2 q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2 +! NRL MNM: Limit q1 + val = epsilon(1.0_kind_phys) + q1(i,k) = max(q1(i,k), val ) ! tem = 1./rcs(i) ! u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem ! v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem diff --git a/physics/sfc_drv.f b/physics/sfc_drv.f index 817897fe7..3b4909594 100644 --- a/physics/sfc_drv.f +++ b/physics/sfc_drv.f @@ -7,6 +7,7 @@ module lsm_noah use machine, only: kind_phys use set_soilveg_mod, only: set_soilveg use namelist_soilveg + use mod_sflx implicit none diff --git a/physics/sflx.f b/physics/sflx.f index 61fe015cc..d0c6dfffd 100644 --- a/physics/sflx.f +++ b/physics/sflx.f @@ -1,3 +1,5 @@ + module mod_sflx + contains !>\file sflx.f !! This file is the entity of GFS Noah LSM Model(Version 2.7). @@ -906,7 +908,14 @@ subroutine gfssflx &! --- input eta = etp endif - beta = eta / etp +! beta = eta / etp +! guard against div zero (from WRF, JM 20211104) + IF (ETP == 0.0) THEN + BETA = 0.0 + ELSE + BETA = ETA/ETP + ENDIF + !> - Convert the sign of soil heat flux so that: !! - ssoil>0: warm the surface (night time) @@ -5801,3 +5810,4 @@ end subroutine wdfcnd end subroutine gfssflx !! @} !----------------------------------- + end module mod_sflx diff --git a/physics/tridi.f b/physics/tridi.f index 0103b388f..fefbc0b7e 100644 --- a/physics/tridi.f +++ b/physics/tridi.f @@ -1,3 +1,5 @@ + module mod_tridi + contains !>\file tridi.f !! These subroutines are originally internal subroutines in moninedmf.f @@ -220,3 +222,4 @@ subroutine tridit(l,n,nt,cl,cm,cu,rt,au,at) return end subroutine tridit !> @} + end module mod_tridi diff --git a/physics/ugwp_driver_v0.F b/physics/ugwp_driver_v0.F index 844acf722..2090951e6 100644 --- a/physics/ugwp_driver_v0.F +++ b/physics/ugwp_driver_v0.F @@ -1,3 +1,6 @@ + module mod_ugwp_driver_v0 + use mod_cires_orowam2017 + contains !>\file ugwp_driver_v0.F ! @@ -1336,7 +1339,7 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime, else kzw2 = zero endif - if ( kzw2 > zero ) then + if ( kzw2 > zero .and. cdf2 > zero ) then ! JM, should always be case but ... 20211020 v_kzw = sqrt(kzw2) ! !linsatdis: kzw2, kzw3, kdsat, c2f2, cdf2, cdf1 @@ -1484,4 +1487,4 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime, return end subroutine fv3_ugwp_solv2_v0 - + end module mod_ugwp_driver_v0