diff --git a/README.md b/README.md index 4aaef30f1..567020bca 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,6 @@ # GFDL_atmos_cubed_sphere -The source contained herein reflects the 201912 release of the Finite Volume Cubed-Sphere Dynamical Core (FV3) from GFDL - -The GFDL Microphysics is also available via this repository. +The source contained herein reflects the 201912 release of the Finite Volume Cubed-Sphere Dynamical Core (FV3) for use in the current GFDL models (AM4/CM4/ESM4/SPEAR). # Where to find information @@ -12,7 +10,6 @@ for more information. # Proper usage attribution Cite Putman and Lin (2007) and Harris and Lin (2013) when describing a model using the FV3 dynamical core. -Cite Chen et al (2013) and Zhou et al (2019) when using the GFDL Microphysics. # What files are what @@ -27,6 +24,13 @@ The top level directory structure groups source code and input files as follow: | ```tools/``` | contains source code of tools used within the core | | ```GFDL_tools/``` | contains source code of tools specific to GFDL models | +# Dependencies + +The source code in this repository requires other NOAA-GFDL projects in +order to compile. As stated above, this branch is for use with the current +NOAA-GFDL/AM4 and NOAA-GFDL/CM4 projects and is included as a submodule within +those projects. + # Disclaimer The United States Department of Commerce (DOC) GitHub project code is provided diff --git a/RELEASE.md b/RELEASE.md index 85f7df54d..e43b13b9f 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -1,31 +1,31 @@ -# RELEASE NOTES for FV3: Summary +# RELEASE NOTES for GFDL FV3: Summary -FV3-201912-public --- 10 January 2020 -Lucas Harris, GFDL +2020.02 --- 22 April 2020 -This version has been tested against the current SHiELD (formerly fvGFS) physics -and with FMS release candidate 2020.02 from https://github.com/NOAA-GFDL/FMS +This version has been tested within current GFDL Models (AM4+, CM4+, ESM4+, SPEAR, etc.) and requires the 2020.02 release of the [FMS infrastructure](https://github.com/NOAA-GFDL/FMS). -Includes all of the features of the GFDL Release to EMC, as well as: +Includes all of the features from the [201912 Public Release](https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/releases/tag/201912_public_release) which include: -- Updated 2017 GFDL Microphysics (from S-J Lin and L Zhou included in GFSv15) -- Updates for GFSv15 ICs (from T Black/J Abeles, EMC) -- Updates to support new nesting capabilities in FMS (from Z Liang) +- Updates to support new nesting capabilities in FMS - Re-written grid nesting code for efficiency and parallelization - Re-organized fv_eta for improved vertical level selection -- 2018 Stand-alone regional capabilities (from T Black/J Abeles, EMC) +- 2018 Stand-alone regional capabilities (from EMC) - Refactored model front-end (fv_control, fv_restart) - Support for point soundings +- full non-hydrostatic capability is now included as a runtime option - And other updates -# Interface changes +# Directory structure changes -drivers: renamed 'fvGFS' directory to SHiELD +***drivers/*** (important for those moving from the GFDL internal project) + - renamed ***fvGFS*** to ***SHiELD*** + - renamed ***coupled*** to ***GFDL*** -atmosphere.F90: 'mytile' is renamed 'mygrid' +***model_nh_null/*** + - has been removed -The non-functional gfdl_cloud_microphys.F90 has been removed and replaced with the 2017 public release given to EMC. Also added a proper initialization routine, that includes the use of INTERNAL_FILE_NML and thereby requires the input_nml_file argument. If you do not define the compiler flag INTERNAL_FILE_NML then you can delete this argument. +Update your build system as appropriate -The namelist nggps_diag_nml has been eliminated. 'fdiag' is no longer handled by the dynamical core, and should be handled by the physics driver. +# Documentation For a complete technical description see the NOAA Technical Memorandum OAR GFDL: https://repository.library.noaa.gov/view/noaa/23432 diff --git a/driver/GFDL/atmosphere.F90 b/driver/GFDL/atmosphere.F90 index 75fe0678f..3c43158b0 100644 --- a/driver/GFDL/atmosphere.F90 +++ b/driver/GFDL/atmosphere.F90 @@ -1264,7 +1264,7 @@ subroutine fv_compute_p_z (npz, phis, pe, peln, delp, delz, pt, q_sph, & endif if (do_uni_zfull) then do k=1,npz - z_full(:,:,k)=0.5*(z_half(:,:,k)+z_half(:,:,k+1)) + z_full(:,:,k)=0.5*(z_half(:,:,k)+z_half(:,:,k+1)) enddo endif end subroutine fv_compute_p_z diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index fddeaf635..6be84f66c 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -60,6 +60,9 @@ module fv_arrays_mod id_dbz, id_maxdbz, id_basedbz, id_dbz4km, id_dbztop, id_dbz_m10C, & id_ctz, id_w1km, id_wmaxup, id_wmaxdn, id_cape, id_cin +! Selected theta-level fields from 3D variables: + integer :: id_pv350K, id_pv550K + ! Selected p-level fields from 3D variables: integer :: id_vort200, id_vort500, id_w500, id_w700 integer :: id_vort850, id_w850, id_x850, id_srh25, & diff --git a/model/fv_sg.F90 b/model/fv_sg.F90 index d283d746d..460d5c989 100644 --- a/model/fv_sg.F90 +++ b/model/fv_sg.F90 @@ -26,7 +26,8 @@ module fv_sg_mod use constants_mod, only: rdgas, rvgas, cp_air, cp_vapor, hlv, hlf, kappa, grav use tracer_manager_mod, only: get_tracer_index use field_manager_mod, only: MODEL_ATMOS - use gfdl_cloud_microphys_mod, only: wqs1, wqs2, wqsat2_moist +! use gfdl_cloud_microphys_mod, only: wqs1, wqs2, wqsat2_moist + use lin_cld_microphys_mod, only: wqs2, wqsat2_moist use fv_mp_mod, only: mp_reduce_min, is_master implicit none diff --git a/model/sw_core.F90 b/model/sw_core.F90 index 99f079ad6..d5df6e8d7 100644 --- a/model/sw_core.F90 +++ b/model/sw_core.F90 @@ -979,8 +979,10 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & ! endif call fv_tp_2d(pt, crx_adv,cry_adv, npx, npy, hord_tm, gx, gy, & xfx_adv,yfx_adv, gridstruct, bd, ra_x, ra_y, & - mfx=fx, mfy=fy, mass=delp, nord=nord_v, damp_c=damp_v) -! mfx=fx, mfy=fy, mass=delp, nord=nord_t, damp_c=damp_t) +! This is a difference from the publicly released version in github master +! and is provided to ensure reproducibility with AM4/CM4 models from GFDL +! mfx=fx, mfy=fy, mass=delp, nord=nord_v, damp_c=damp_v) + mfx=fx, mfy=fy, mass=delp, nord=nord_t, damp_c=damp_t) #endif if ( inline_q ) then diff --git a/tools/fv_diagnostics.F90 b/tools/fv_diagnostics.F90 index bac8fa440..58dc2ef4d 100644 --- a/tools/fv_diagnostics.F90 +++ b/tools/fv_diagnostics.F90 @@ -46,7 +46,8 @@ module fv_diagnostics_mod use sat_vapor_pres_mod, only: compute_qs, lookup_es use fv_arrays_mod, only: max_step - use gfdl_cloud_microphys_mod, only: wqs1, qsmith_init +! use gfdl_cloud_microphys_mod, only: wqs1, qsmith_init + use lin_cld_microphys_mod, only: wqs1, qsmith_init use column_diagnostics_mod, only: column_diagnostics_init, & initialize_diagnostic_columns, & @@ -707,6 +708,10 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) !-------------------- idiag%id_pv = register_diag_field ( trim(field), 'pv', axes(1:3), Time, & 'potential vorticity', '1/s', missing_value=missing_value ) + idiag%id_pv350K = register_diag_field ( trim(field), 'pv350K', axes(1:2), Time, & + '350-K potential vorticity; needs x350 scaling', '(K m**2) / (kg s)', missing_value=missing_value) + idiag%id_pv550K = register_diag_field ( trim(field), 'pv550K', axes(1:2), Time, & + '550-K potential vorticity; needs x550 scaling', '(K m**2) / (kg s)', missing_value=missing_value) ! ------------------- ! Vertical flux correlation terms (good for averages) @@ -1576,8 +1581,8 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) endif if ( idiag%id_vort200>0 .or. idiag%id_vort500>0 .or. idiag%id_vort850>0 .or. idiag%id_vorts>0 & - .or. idiag%id_vort>0 .or. idiag%id_pv>0 .or. idiag%id_rh>0 .or. idiag%id_x850>0 .or. & - idiag%id_uh03>0 .or. idiag%id_uh25>0) then + .or. idiag%id_vort>0 .or. idiag%id_pv>0 .or. idiag%id_pv350K>0 .or. idiag%id_pv550K>0 & + .or. idiag%id_rh>0 .or. idiag%id_x850>0 .or. idiag%id_uh03>0 .or. idiag%id_uh25>0) then call get_vorticity(isc, iec, jsc, jec, isd, ied, jsd, jed, npz, Atm(n)%u, Atm(n)%v, wk, & Atm(n)%gridstruct%dx, Atm(n)%gridstruct%dy, Atm(n)%gridstruct%rarea) @@ -1731,11 +1736,36 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) endif - if ( idiag%id_pv > 0 ) then -! Note: this is expensive computation. + if ( idiag%id_pv > 0 .or. idiag%id_pv350K >0 .or. idiag%id_pv550K >0) then + if (allocated(a3)) deallocate(a3) + allocate ( a3(isc:iec,jsc:jec,npz+1) ) + ! Modified pv_entropy to get potential temperature at layer interfaces (last variable) + ! The values are needed for interpolate_z + ! Note: this is expensive computation. call pv_entropy(isc, iec, jsc, jec, ngc, npz, wk, & - Atm(n)%gridstruct%f0, Atm(n)%pt, Atm(n)%pkz, Atm(n)%delp, grav) - used = send_data ( idiag%id_pv, wk, Time ) + Atm(n)%gridstruct%f0, Atm(n)%pt, Atm(n)%pkz, Atm(n)%delp, grav,a3) + if ( idiag%id_pv > 0) then + used = send_data ( idiag%id_pv, wk, Time ) + endif + if( idiag%id_pv350K > 0 .or. idiag%id_pv550K >0 ) then + !"pot temp" from pv_entropy is only semi-finished; needs p0^kappa (pk0) + do k=1,npz+1 + do j=jsc,jec + do i=isc,iec + a3(i,j,k) = a3(i,j,k) * pk0 + enddo + enddo + enddo + !wk as input, which stores pv from pv_entropy; + !use z interpolation, both potential temp and z increase monotonically with height + !interpolate_vertical will apply log operation to 350, and also assumes a different vertical order + call interpolate_z(isc, iec, jsc, jec, npz, 350., a3, wk, a2) + used = send_data( idiag%id_pv350K, a2, Time) + !interpolate_vertical will apply log operation to 350, and also assumes a different vertical order + call interpolate_z(isc, iec, jsc, jec, npz, 550., a3, wk, a2) + used = send_data( idiag%id_pv550K, a2, Time) + endif + deallocate ( a3 ) if (prt_minmax) call prt_maxmin('PV', wk, isc, iec, jsc, jec, 0, 1, 1.) endif @@ -4236,6 +4266,7 @@ subroutine get_height_given_pressure(is, ie, js, je, km, wz, kd, id, log_p, peln go to 1000 endif enddo + a2(i,j,n) =missing_value 1000 continue enddo enddo @@ -4919,7 +4950,7 @@ end subroutine updraft_helicity - subroutine pv_entropy(is, ie, js, je, ng, km, vort, f_d, pt, pkz, delp, grav) + subroutine pv_entropy(is, ie, js, je, ng, km, vort, f_d, pt, pkz, delp, grav, te) ! !INPUT PARAMETERS: integer, intent(in):: is, ie, js, je, ng, km @@ -4930,7 +4961,9 @@ subroutine pv_entropy(is, ie, js, je, ng, km, vort, f_d, pt, pkz, delp, grav) real, intent(in):: f_d(is-ng:ie+ng,js-ng:je+ng) ! vort is relative vorticity as input. Becomes PV on output - real, intent(inout):: vort(is:ie,js:je,km) + real, intent(inout):: vort(is:ie,js:je,km) +! output potential temperature at the interface so it can be used for diagnostics + real, intent(out):: te(is:ie,js:je,km+1) ! !DESCRIPTION: ! EPV = 1/r * (vort+f_d) * d(S)/dz; where S is a conservative scalar @@ -4941,7 +4974,7 @@ subroutine pv_entropy(is, ie, js, je, ng, km, vort, f_d, pt, pkz, delp, grav) ! z-surface is not that different from the hybrid sigma-p coordinate. ! See page 39, Pedlosky 1979: Geophysical Fluid Dynamics ! -! The follwoing simplified form is strictly correct only if vort is computed on +! The following simplified form is strictly correct only if vort is computed on ! constant z surfaces. In addition hydrostatic approximation is made. ! EPV = - GRAV * (vort+f_d) / del(p) * del(pt) / pt ! where del() is the vertical difference operator. @@ -4952,7 +4985,7 @@ subroutine pv_entropy(is, ie, js, je, ng, km, vort, f_d, pt, pkz, delp, grav) !--------------------------------------------------------------------- !BOC real w3d(is:ie,js:je,km) - real te(is:ie,js:je,km+1), t2(is:ie,km), delp2(is:ie,km) + real t2(is:ie,km), delp2(is:ie,km) real te2(is:ie,km+1) integer i, j, k @@ -4965,8 +4998,8 @@ subroutine pv_entropy(is, ie, js, je, ng, km, vort, f_d, pt, pkz, delp, grav) enddo #else ! Compute PT at layer edges. -!$OMP parallel do default(none) shared(is,ie,js,je,km,pt,pkz,w3d,delp,te2,te) & -!$OMP private(t2, delp2) +!$OMP parallel do default(none) shared(is,ie,js,je,km,pt,pkz,w3d,delp,te) & +!$OMP private(t2, te2, delp2) do j=js,je do k=1,km do i=is,ie @@ -4985,7 +5018,8 @@ subroutine pv_entropy(is, ie, js, je, ng, km, vort, f_d, pt, pkz, delp, grav) enddo enddo -!$OMP parallel do default(none) shared(is,ie,js,je,km,vort,f_d,te,w3d,delp,grav) +!$OMP parallel do default(none) shared(is,ie,js,je,km,vort,te,w3d,delp,grav) & +!$OMP private(f_d) do k=1,km do j=js,je do i=is,ie diff --git a/tools/fv_mp_mod.F90 b/tools/fv_mp_mod.F90 index ced71e3f9..fd9881f68 100644 --- a/tools/fv_mp_mod.F90 +++ b/tools/fv_mp_mod.F90 @@ -28,7 +28,7 @@ module fv_mp_mod ! !USES: use fms_mod, only : fms_init, fms_end use mpp_mod, only : FATAL, MPP_DEBUG, NOTE, MPP_CLOCK_SYNC,MPP_CLOCK_DETAILED, WARNING - use mpp_mod, only : mpp_pe, mpp_npes, mpp_node, mpp_root_pe, mpp_error, mpp_set_warn_level + use mpp_mod, only : mpp_pe, mpp_npes, mpp_root_pe, mpp_error, mpp_set_warn_level use mpp_mod, only : mpp_declare_pelist, mpp_set_current_pelist, mpp_sync use mpp_mod, only : mpp_clock_begin, mpp_clock_end, mpp_clock_id use mpp_mod, only : mpp_chksum, stdout, stderr, mpp_broadcast diff --git a/tools/fv_nudge.F90 b/tools/fv_nudge.F90 index b1a5c1e88..057358136 100644 --- a/tools/fv_nudge.F90 +++ b/tools/fv_nudge.F90 @@ -61,6 +61,7 @@ module fv_nwp_nudge_mod public fv_nwp_nudge, fv_nwp_nudge_init, fv_nwp_nudge_end, breed_slp_inline, T_is_Tv public do_adiabatic_init + public nwp_nudge_int integer im ! Data x-dimension integer jm ! Data y-dimension integer km ! Data z-dimension @@ -112,6 +113,8 @@ module fv_nwp_nudge_mod real :: p_wvp = 100.E2 ! cutoff level for specific humidity nudging integer :: kord_data = 8 + integer :: nwp_nudge_int = 21600 ! 6 hours + real :: mask_fac = 0.25 ! [0,1] 0: no mask; 1: full strength logical :: T_is_Tv = .false. @@ -131,8 +134,11 @@ module fv_nwp_nudge_mod logical :: nudge_virt = .true. logical :: nudge_hght = .true. logical :: time_varying = .true. + logical :: time_varying_nwp = .false. logical :: print_end_breed = .true. logical :: print_end_nudge = .true. + logical :: using_merra2 = .false. ! Flag to allow avoidance of multiplicative factor if using MERRA2 data. + logical :: climate_nudging = .false. ! Flag to allow for climate nudging. ! Nudging time-scales (seconds): note, however, the effective time-scale is 2X smaller (stronger) due @@ -216,7 +222,8 @@ module fv_nwp_nudge_mod kbot_t, kbot_q, p_wvp, time_varying, time_interval, use_pt_inc, pt_lim, & tau_vt_rad, r_lo, r_hi, use_high_top, add_bg_wind, conserve_mom, conserve_hgt, & min_nobs, min_mslp, nudged_time, r_fac, r_min, r_inc, ibtrack, track_file_name, file_names, & - input_fname_list, analysis_file_first, analysis_file_last, P_relax, P_norelax !h1g, add 3 namelist variables, 2012-20-22 + input_fname_list, analysis_file_first, analysis_file_last, P_relax, P_norelax, & + nwp_nudge_int, time_varying_nwp, using_merra2, climate_nudging contains @@ -241,7 +248,8 @@ subroutine fv_nwp_nudge ( Time, dt, npx, npy, npz, ps_dt, u_dt, v_dt, t_dt, q_dt real, intent(inout), dimension(isd:ied,jsd:jed):: ps ! Accumulated tendencies real, intent(inout), dimension(isd:ied,jsd:jed,npz):: u_dt, v_dt - real, intent(out), dimension(is:ie,js:je,npz):: t_dt, q_dt + real, intent(out), dimension(is:ie,js:je,npz):: t_dt + real, intent(inout), dimension(is:ie,js:je,npz,1):: q_dt real, intent(out), dimension(is:ie,js:je):: ps_dt, ts type(fv_grid_type), intent(INOUT), target :: gridstruct @@ -262,6 +270,7 @@ subroutine fv_nwp_nudge ( Time, dt, npx, npy, npz, ps_dt, u_dt, v_dt, t_dt, q_dt integer :: seconds, days integer :: i,j,k, iq, kht real :: factor, rms, bias, co + real :: factor_nwp real :: rdt, press(npz), profile(npz), prof_t(npz), prof_q(npz), du, dv logical used @@ -339,22 +348,43 @@ subroutine fv_nwp_nudge ( Time, dt, npx, npy, npz, ps_dt, u_dt, v_dt, t_dt, q_dt ! Thermodynamics: prof_t(:) = 1. + if (climate_nudging) then +!$OMP parallel do default(none) shared(npz,press,P_norelax,prof_t) + do k=1,npz + if ( press(k) < 10.E2 ) then + prof_t(k) = max(0.01, press(k)/10.E2) + endif + ! above P_norelax, no nudging. + if( press(k) < P_norelax ) prof_t(k) = 0.0 + enddo + else !$OMP parallel do default(none) shared(npz,press,prof_t) - do k=1,npz - if ( press(k) < 10.E2 ) then + do k=1,npz + if ( press(k) < 10.E2 ) then prof_t(k) = max(0.01, press(k)/10.E2) - endif - enddo + endif + enddo + endif prof_t(1) = 0. ! Water vapor: prof_q(:) = 1. + if ( climate_nudging) then +!$OMP parallel do default(none) shared(npz,press,P_norelax,prof_q) + do k=1,npz + if ( press(k) < 200.E2 ) then + prof_q(k) = max(0., press(k)/200.E2) + endif + if( press(k) < P_norelax ) prof_q(k) = 0.0 + enddo + else !$OMP parallel do default(none) shared(npz,press,prof_q) - do k=1,npz - if ( press(k) < 300.E2 ) then + do k=1,npz + if ( press(k) < 300.E2 ) then prof_q(k) = max(0., press(k)/300.E2) - endif - enddo + endif + enddo + endif prof_q(1) = 0. ! Height @@ -382,6 +412,16 @@ subroutine fv_nwp_nudge ( Time, dt, npx, npy, npz, ps_dt, u_dt, v_dt, t_dt, q_dt factor = 1. endif + if ( time_varying_nwp ) then + if (mod(seconds, nwp_nudge_int) == 0) then + factor_nwp = 1.0 + else + factor_nwp = 0.0 + endif + else + factor_nwp = 1.0 + endif + if ( do_adiabatic_init ) factor = 2.*factor allocate (ps_obs(is:ie,js:je) ) @@ -395,7 +435,7 @@ subroutine fv_nwp_nudge ( Time, dt, npx, npy, npz, ps_dt, u_dt, v_dt, t_dt, q_dt call get_obs(Time, dt, zvir, ak, bk, ps, ts, ps_obs, delp, pt, nwat, q, u_obs, v_obs, t_obs, q_obs, & - phis, ua, va, u_dt, v_dt, npx, npy, npz, factor, mask, ptop, bd, gridstruct, domain) + phis, ua, va, u_dt, v_dt, npx, npy, npz, factor, factor_nwp, mask, ptop, bd, gridstruct, domain) ! *t_obs* is virtual temperature if ( no_obs ) then @@ -486,34 +526,48 @@ subroutine fv_nwp_nudge ( Time, dt, npx, npy, npz, ps_dt, u_dt, v_dt, t_dt, q_dt if ( nf_uv>0 ) call del2_uv(du_obs, dv_obs, del2_cd, npz, nf_uv, bd, npx, npy, gridstruct, domain) !$OMP parallel do default(none) shared(kstart,kbot_winds,npz,is,ie,js,je,du_obs,dv_obs, & -!$OMP mask,ps_fac,u_dt,v_dt,ua,va,dt) +!$OMP mask,ps_fac,u_dt,v_dt,ua,va,dt,climate_nudging) do k=kstart, npz - kbot_winds - if ( k==npz ) then - do j=js,je + if ( climate_nudging) then + do j=js,je do i=is,ie - du_obs(i,j,k) = du_obs(i,j,k) * mask(i,j) * ps_fac(i,j) - dv_obs(i,j,k) = dv_obs(i,j,k) * mask(i,j) * ps_fac(i,j) +! Apply TC mask + du_obs(i,j,k) = du_obs(i,j,k) * mask(i,j) + dv_obs(i,j,k) = dv_obs(i,j,k) * mask(i,j) + u_dt(i,j,k) = u_dt(i,j,k) + du_obs(i,j,k) + v_dt(i,j,k) = v_dt(i,j,k) + dv_obs(i,j,k) + ua(i,j,k) = ua(i,j,k) + du_obs(i,j,k)*dt + va(i,j,k) = va(i,j,k) + dv_obs(i,j,k)*dt + enddo + enddo + else + if ( k==npz ) then + do j=js,je + do i=is,ie + du_obs(i,j,k) = du_obs(i,j,k) * mask(i,j) * ps_fac(i,j) + dv_obs(i,j,k) = dv_obs(i,j,k) * mask(i,j) * ps_fac(i,j) ! - u_dt(i,j,k) = u_dt(i,j,k) + du_obs(i,j,k) - v_dt(i,j,k) = v_dt(i,j,k) + dv_obs(i,j,k) - ua(i,j,k) = ua(i,j,k) + du_obs(i,j,k)*dt - va(i,j,k) = va(i,j,k) + dv_obs(i,j,k)*dt + u_dt(i,j,k) = u_dt(i,j,k) + du_obs(i,j,k) + v_dt(i,j,k) = v_dt(i,j,k) + dv_obs(i,j,k) + ua(i,j,k) = ua(i,j,k) + du_obs(i,j,k)*dt + va(i,j,k) = va(i,j,k) + dv_obs(i,j,k)*dt + enddo enddo - enddo - else - do j=js,je - do i=is,ie + else + do j=js,je + do i=is,ie ! Apply TC mask - du_obs(i,j,k) = du_obs(i,j,k) * mask(i,j) - dv_obs(i,j,k) = dv_obs(i,j,k) * mask(i,j) + du_obs(i,j,k) = du_obs(i,j,k) * mask(i,j) + dv_obs(i,j,k) = dv_obs(i,j,k) * mask(i,j) ! - u_dt(i,j,k) = u_dt(i,j,k) + du_obs(i,j,k) - v_dt(i,j,k) = v_dt(i,j,k) + dv_obs(i,j,k) - ua(i,j,k) = ua(i,j,k) + du_obs(i,j,k)*dt - va(i,j,k) = va(i,j,k) + dv_obs(i,j,k)*dt + u_dt(i,j,k) = u_dt(i,j,k) + du_obs(i,j,k) + v_dt(i,j,k) = v_dt(i,j,k) + dv_obs(i,j,k) + ua(i,j,k) = ua(i,j,k) + du_obs(i,j,k)*dt + va(i,j,k) = va(i,j,k) + dv_obs(i,j,k)*dt + enddo enddo - enddo - endif + endif + endif ! climate_nudging enddo endif @@ -529,28 +583,36 @@ subroutine fv_nwp_nudge ( Time, dt, npx, npy, npz, ps_dt, u_dt, v_dt, t_dt, q_dt if ( nudge_virt ) then rdt = 1./(tau_virt/factor + dt) !$OMP parallel do default(none) shared(is,ie,js,je,npz,kstart,kht,t_dt,prof_t,t_obs,zvir, & -!$OMP q,pt,rdt,ps_fac) +!$OMP q,pt,rdt,ps_fac,climate_nudging) do k=kstart, kht - if ( k==npz ) then - do j=js,je + if ( climate_nudging) then + do j=js,je do i=is,ie - t_dt(i,j,k) = prof_t(k)*(t_obs(i,j,k)/(1.+zvir*q(i,j,k,1))-pt(i,j,k))*rdt*ps_fac(i,j) + t_dt(i,j,k) = prof_t(k)*(t_obs(i,j,k)-pt(i,j,k))*rdt enddo - enddo - else - do j=js,je - do i=is,ie - t_dt(i,j,k) = prof_t(k)*(t_obs(i,j,k)/(1.+zvir*q(i,j,k,1))-pt(i,j,k))*rdt + enddo + else + if ( k==npz ) then + do j=js,je + do i=is,ie + t_dt(i,j,k) = prof_t(k)*(t_obs(i,j,k)/(1.+zvir*q(i,j,k,1))-pt(i,j,k))*rdt*ps_fac(i,j) + enddo enddo - enddo - endif + else + do j=js,je + do i=is,ie + t_dt(i,j,k) = prof_t(k)*(t_obs(i,j,k)/(1.+zvir*q(i,j,k,1))-pt(i,j,k))*rdt + enddo + enddo + endif + endif enddo endif if ( nudge_hght .and. kht p_wvp ) then do iq=2,nwat @@ -618,8 +685,13 @@ subroutine fv_nwp_nudge ( Time, dt, npx, npy, npz, ps_dt, u_dt, v_dt, t_dt, q_dt do j=js,je do i=is,ie delp(i,j,k) = delp(i,j,k)*(1.-q(i,j,k,1)) - q_dt(i,j,k) = prof_q(k)*(max(q_min,q_obs(i,j,k))-q(i,j,k,1))*rdt*mask(i,j) - q(i,j,k,1) = q(i,j,k,1) + q_dt(i,j,k)*dt + if ( climate_nudging ) then + q_dt(i,j,k,1) = q_dt(i,j,k,1)+prof_q(k)*(max(q_min,q_obs(i,j,k))-q(i,j,k,1))*rdt + q(i,j,k,1) = q(i,j,k,1) + prof_q(k)*(max(q_min,q_obs(i,j,k))-q(i,j,k,1))*rdt*dt + else + q_dt(i,j,k,1) = prof_q(k)*(max(q_min,q_obs(i,j,k))-q(i,j,k,1))*rdt*mask(i,j) + q(i,j,k,1) = q(i,j,k,1) + q_dt(i,j,k,1)*dt + endif delp(i,j,k) = delp(i,j,k)/(1.-q(i,j,k,1)) enddo enddo @@ -681,9 +753,9 @@ subroutine fv_nwp_nudge ( Time, dt, npx, npy, npz, ps_dt, u_dt, v_dt, t_dt, q_dt end subroutine fv_nwp_nudge - subroutine ps_nudging(dt, factor, npz, ak, bk, ps_obs, mask, tm, ps, phis, delp, ua, va, pt, nwat, q, bd, npx, npy, gridstruct, domain) + subroutine ps_nudging(dt, factor, factor_nwp, npz, ak, bk, ps_obs, mask, tm, ps, phis, delp, ua, va, pt, nwat, q, bd, npx, npy, gridstruct, domain) ! Input - real, intent(in):: dt, factor + real, intent(in):: dt, factor, factor_nwp integer, intent(in):: npz, nwat, npx, npy real, intent(in), dimension(npz+1):: ak, bk type(fv_grid_bounds_type), intent(IN) :: bd @@ -791,7 +863,7 @@ subroutine ps_nudging(dt, factor, npz, ak, bk, ps_obs, mask, tm, ps, phis, delp, enddo enddo - rdt = dt / (tau_ps/factor + dt) + rdt = factor_nwp*dt / (tau_ps/factor + dt) do k=1,npz dbk = rdt*(bk(k+1) - bk(k)) do j=js,je @@ -957,11 +1029,11 @@ end subroutine compute_slp subroutine get_obs(Time, dt, zvir, ak, bk, ps, ts, ps_obs, delp, pt, nwat, q, u_obs, v_obs, t_obs, q_obs, & - phis, ua, va, u_dt, v_dt, npx, npy, npz, factor, mask, ptop, bd, gridstruct, domain) + phis, ua, va, u_dt, v_dt, npx, npy, npz, factor, factor_nwp, mask, ptop, bd, gridstruct, domain) type(time_type), intent(in):: Time integer, intent(in):: npz, nwat, npx, npy real, intent(in):: zvir, ptop - real, intent(in):: dt, factor + real, intent(in):: dt, factor, factor_nwp real, intent(in), dimension(npz+1):: ak, bk type(fv_grid_bounds_type), intent(IN) :: bd real, intent(in), dimension(isd:ied,jsd:jed):: phis @@ -1069,7 +1141,7 @@ subroutine get_obs(Time, dt, zvir, ak, bk, ps, ts, ps_obs, delp, pt, nwat, q, u_ allocate ( vv(isd:ied,jsd:jed,npz) ) uu = ua vv = va - call ps_nudging(dt, factor, npz, ak, bk, ps_obs, mask, tm, ps, phis, delp, uu, vv, pt, nwat, q, bd, npx, npy, gridstruct, domain) + call ps_nudging(dt, factor, factor_nwp, npz, ak, bk, ps_obs, mask, tm, ps, phis, delp, uu, vv, pt, nwat, q, bd, npx, npy, gridstruct, domain) do k=1,npz do j=js,je do i=is,ie @@ -1203,7 +1275,7 @@ subroutine fv_nwp_nudge_init(time, axes, npz, zvir, ak, bk, ts, phis, gridstruct if( trim(fname_tmp) .ne. "" ) then ! escape any empty record if ( trim(fname_tmp) == trim(analysis_file_last) ) then nt = nt + 1 - file_names(nt) = 'INPUT/'//trim(fname_tmp) + file_names(nt) = trim(fname_tmp) if(master .and. nudge_debug) write(*,*) 'From NCEP file list, last file: ', nt, file_names(nt) nt = 0 goto 101 ! read last analysis data and then close file @@ -1211,7 +1283,7 @@ subroutine fv_nwp_nudge_init(time, axes, npz, zvir, ak, bk, ts, phis, gridstruct if ( trim(analysis_file_first) == "" ) then nt = nt + 1 - file_names(nt) = 'INPUT/'//trim(fname_tmp) + file_names(nt) = trim(fname_tmp) if(master .and. nudge_debug) then if( nt .eq. 1 ) then write(*,*) 'From NCEP file list, first file: ', nt, file_names(nt),trim(analysis_file_first) @@ -1222,7 +1294,7 @@ subroutine fv_nwp_nudge_init(time, axes, npz, zvir, ak, bk, ts, phis, gridstruct else if ( trim(fname_tmp) == trim(analysis_file_first) .or. nt > 0 ) then nt = nt + 1 - file_names(nt) = 'INPUT/'//trim(fname_tmp) + file_names(nt) = trim(fname_tmp) if(master .and. nudge_debug) then if( nt .eq. 1 ) then write(*,*) 'From NCEP file list, first file: ', nt, file_names(nt),trim(analysis_file_first) @@ -1294,7 +1366,10 @@ subroutine fv_nwp_nudge_init(time, axes, npz, zvir, ak, bk, ts, phis, gridstruct call close_ncfile( ncid ) ! Note: definition of NCEP hybrid is p(k) = a(k)*1.E5 + b(k)*ps - ak0(:) = ak0(:) * 1.E5 + if ( .not. using_merra2) then + ! This is not needed for MERRA2 data + ak0(:) = ak0(:) * 1.E5 + endif ! Limiter to prevent NAN at top during remapping if ( bk0(1) < 1.E-9 ) ak0(1) = max(1.e-9, ak0(1)) @@ -1375,6 +1450,7 @@ subroutine get_ncep_analysis ( ps, u, v, t, q, zvir, ts, nfile, fname, bd ) if(master) write(*,*) 'Reading NCEP anlysis file:', fname endif + if ( climate_nudging ) read_ts =.false. if ( read_ts ) then ! read skin temperature; could be used for SST allocate ( wk1(im,jm) ) @@ -1572,6 +1648,7 @@ subroutine get_ncep_analysis ( ps, u, v, t, q, zvir, ts, nfile, fname, bd ) enddo enddo + if ( .not. climate_nudging) then if ( .not. T_is_Tv ) then do k=1,km do j=js,je @@ -1584,6 +1661,7 @@ subroutine get_ncep_analysis ( ps, u, v, t, q, zvir, ts, nfile, fname, bd ) enddo enddo endif + endif ! endif diff --git a/tools/test_cases.F90 b/tools/test_cases.F90 index 480247ae0..9311c5f18 100644 --- a/tools/test_cases.F90 +++ b/tools/test_cases.F90 @@ -7162,7 +7162,8 @@ end subroutine superK_u subroutine SuperCell_Sounding(km, ps, pk1, tp, qp) - use gfdl_cloud_microphys_mod, only: wqsat_moist, qsmith_init, qs_blend +! use gfdl_cloud_microphys_mod, only: wqsat_moist, qsmith_init, qs_blend + use lin_cld_microphys_mod, only: wqsat_moist, qsmith_init, qs_blend ! Morris Weisman & J. Klemp 2002 sounding ! Output sounding on pressure levels: integer, intent(in):: km