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/tools/fv_diagnostics.F90 b/tools/fv_diagnostics.F90 index 83d0b58a7..58dc2ef4d 100644 --- a/tools/fv_diagnostics.F90 +++ b/tools/fv_diagnostics.F90 @@ -708,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) @@ -1577,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) @@ -1732,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 @@ -4237,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 @@ -4920,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 @@ -4931,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 @@ -4942,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. @@ -4953,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 @@ -4966,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 @@ -4986,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_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