diff --git a/Externals_CAM.cfg b/Externals_CAM.cfg
index 90fe4700d5..cd7c03e349 100644
--- a/Externals_CAM.cfg
+++ b/Externals_CAM.cfg
@@ -21,16 +21,16 @@ required = True
[clubb]
local_path = src/physics/clubb
-protocol = git
-repo_url = https://github.com/ESCOMP/CLUBB_CESM
-tag = clubb_release_b76a124_20200220_c20200320
+protocol = svn
+repo_url = https://github.com/larson-group/clubb_release/tags/
+tag = clubb_4ncar_20220311_f51de38/src/CLUBB_core
required = True
[silhs]
local_path = src/physics/silhs
-protocol = git
-repo_url = https://github.com/ESCOMP/SILHS_CESM
-tag = silhs_clubb_release_b76a124_20200220_c20200320
+protocol = svn
+repo_url = https://github.com/larson-group/clubb_release/tags/
+tag = clubb_4ncar_20220311_f51de38/src/SILHS
required = True
[pumas]
diff --git a/bld/build-namelist b/bld/build-namelist
index 6d260e601d..5023742114 100755
--- a/bld/build-namelist
+++ b/bld/build-namelist
@@ -3080,55 +3080,112 @@ if ($clubb_sgs =~ /$TRUE/io) {
}
add_default($nl, 'clubb_do_icesuper');
-
- add_default($nl, 'clubb_expldiff');
- add_default($nl, 'clubb_rainevap_turb');
+ add_default($nl, 'clubb_do_energyfix');
add_default($nl, 'clubb_cloudtop_cooling');
- add_default($nl, 'clubb_timestep');
+ add_default($nl, 'clubb_rainevap_turb');
add_default($nl, 'clubb_rnevap_effic');
- add_default($nl, 'clubb_beta');
- add_default($nl, 'clubb_c1');
- add_default($nl, 'clubb_c1b');
- add_default($nl, 'clubb_c11');
- add_default($nl, 'clubb_c11b');
- add_default($nl, 'clubb_c14');
- add_default($nl, 'clubb_C2rt');
- add_default($nl, 'clubb_C2thl');
- add_default($nl, 'clubb_C2rtthl');
- add_default($nl, 'clubb_C4');
- add_default($nl, 'clubb_c6rt');
- add_default($nl, 'clubb_c6rtb');
- add_default($nl, 'clubb_c6rtc');
- add_default($nl, 'clubb_c6thl');
- add_default($nl, 'clubb_c6thlb');
- add_default($nl, 'clubb_c6thlc');
+ add_default($nl, 'clubb_timestep');
+ add_default($nl, 'clubb_l_diag_Lscale_from_tau');
+
+ my $clubb_Lscale_from_tau = $nl->get_value('clubb_l_diag_Lscale_from_tau');
+
+ if($clubb_Lscale_from_tau =~ "true") {
+ add_default($nl, 'clubb_c1', 'val'=>1.0);
+ add_default($nl, 'clubb_c1b', 'val'=>1.0);
+ add_default($nl, 'clubb_C2rt', 'val'=>1.0);
+ add_default($nl, 'clubb_C2thl', 'val'=>1.0);
+ add_default($nl, 'clubb_C2rtthl', 'val'=>1.0);
+ add_default($nl, 'clubb_C4', 'val'=>5.2);
+ add_default($nl, 'clubb_C_uu_shr', 'val'=>0.1076484659222455);
+ add_default($nl, 'clubb_C_uu_buoy', 'val'=>0.3);
+ add_default($nl, 'clubb_c6rt', 'val'=>2.0);
+ add_default($nl, 'clubb_c6rtb', 'val'=>2.0);
+ add_default($nl, 'clubb_c6rtc', 'val'=>1.0);
+ add_default($nl, 'clubb_c6thl', 'val'=>2.0);
+ add_default($nl, 'clubb_c6thlb', 'val'=>2.0);
+ add_default($nl, 'clubb_c6thlc', 'val'=>1.0);
+ add_default($nl, 'clubb_C8', 'val'=>3.440377776099962);
+ add_default($nl, 'clubb_C8b', 'val'=>0.0);
+ add_default($nl, 'clubb_c11', 'val'=>0.31057411754034614);
+ add_default($nl, 'clubb_c11b', 'val'=>0.3250718127387944);
+ add_default($nl, 'clubb_c14', 'val'=>1.0);
+ add_default($nl, 'clubb_C_invrs_tau_bkgnd', 'val'=>3.727123755772682);
+ add_default($nl, 'clubb_C_invrs_tau_sfc', 'val'=>0.12743072568015346);
+ add_default($nl, 'clubb_C_invrs_tau_shear', 'val'=>0.12502726304767026);
+ add_default($nl, 'clubb_C_invrs_tau_N2', 'val'=>0.08122667220596895);
+ add_default($nl, 'clubb_C_invrs_tau_N2_wp2', 'val'=>0.1);
+ add_default($nl, 'clubb_C_invrs_tau_N2_xp2', 'val'=>0.05);
+ add_default($nl, 'clubb_C_invrs_tau_N2_wpxp', 'val'=>0.0);
+ add_default($nl, 'clubb_C_invrs_tau_N2_clear_wp3', 'val'=>1.0);
+ add_default($nl, 'clubb_gamma_coef', 'val'=>0.5492223674353673);
+ add_default($nl, 'clubb_gamma_coefb', 'val'=>0.2531868210746816);
+ add_default($nl, 'clubb_beta', 'val'=>2.27756371212011);
+ } else {
+ add_default($nl, 'clubb_c1');
+ add_default($nl, 'clubb_c1b');
+ add_default($nl, 'clubb_C2rt');
+ add_default($nl, 'clubb_C2thl');
+ add_default($nl, 'clubb_C2rtthl');
+ add_default($nl, 'clubb_C4');
+ add_default($nl, 'clubb_C_uu_shr');
+ add_default($nl, 'clubb_C_uu_buoy');
+ add_default($nl, 'clubb_c6rt');
+ add_default($nl, 'clubb_c6rtb');
+ add_default($nl, 'clubb_c6rtc');
+ add_default($nl, 'clubb_c6thl');
+ add_default($nl, 'clubb_c6thlb');
+ add_default($nl, 'clubb_c6thlc');
+ add_default($nl, 'clubb_C8');
+ add_default($nl, 'clubb_C8b');
+ add_default($nl, 'clubb_c11');
+ add_default($nl, 'clubb_c11b');
+ add_default($nl, 'clubb_c14');
+ add_default($nl, 'clubb_C_invrs_tau_bkgnd');
+ add_default($nl, 'clubb_C_invrs_tau_sfc');
+ add_default($nl, 'clubb_C_invrs_tau_shear');
+ add_default($nl, 'clubb_C_invrs_tau_N2');
+ add_default($nl, 'clubb_C_invrs_tau_N2_wp2');
+ add_default($nl, 'clubb_C_invrs_tau_N2_xp2');
+ add_default($nl, 'clubb_C_invrs_tau_N2_wpxp');
+ add_default($nl, 'clubb_C_invrs_tau_N2_clear_wp3');
+ add_default($nl, 'clubb_gamma_coef');
+ add_default($nl, 'clubb_gamma_coefb');
+ add_default($nl, 'clubb_beta');
+ }
+
add_default($nl, 'clubb_C7');
add_default($nl, 'clubb_C7b');
- add_default($nl, 'clubb_C8');
- add_default($nl, 'clubb_C8b');
+
+ add_default($nl, 'clubb_C_wp3_pr_turb');
+ add_default($nl, 'clubb_c_K1');
+ add_default($nl, 'clubb_c_K2');
+ add_default($nl, 'clubb_nu2');
+ add_default($nl, 'clubb_c_K8');
add_default($nl, 'clubb_c_K9');
add_default($nl, 'clubb_nu9');
add_default($nl, 'clubb_c_K10');
add_default($nl, 'clubb_c_K10h');
add_default($nl, 'clubb_do_liqsupersat');
- add_default($nl, 'clubb_gamma_coef');
- add_default($nl, 'clubb_gamma_coefb');
+ add_default($nl, 'clubb_wpxp_L_thresh');
+
add_default($nl, 'clubb_lambda0_stability_coef');
add_default($nl, 'clubb_lmin_coef');
add_default($nl, 'clubb_mult_coef');
add_default($nl, 'clubb_Skw_denom_coef');
add_default($nl, 'clubb_skw_max_mag');
- add_default($nl, 'clubb_up2_vp2_factor');
+ add_default($nl, 'clubb_up2_sfc_coef');
add_default($nl, 'clubb_C_wp2_splat');
- add_default($nl, 'clubb_wpxp_L_thresh');
add_default($nl, 'clubb_detliq_rad');
add_default($nl, 'clubb_detice_rad');
add_default($nl, 'clubb_detphase_lowtemp');
+ add_default($nl, 'clubb_ipdf_call_placement');
add_default($nl, 'clubb_l_brunt_vaisala_freq_moist');
add_default($nl, 'clubb_l_call_pdf_closure_twice');
add_default($nl, 'clubb_l_damp_wp3_Skw_squared');
+ add_default($nl, 'clubb_l_lmm_stepping');
+ add_default($nl, 'clubb_l_e3sm_config');
add_default($nl, 'clubb_l_lscale_plume_centered');
add_default($nl, 'clubb_l_min_wp2_from_corr_wx');
add_default($nl, 'clubb_l_min_xp2_from_corr_wx');
@@ -3141,11 +3198,18 @@ if ($clubb_sgs =~ /$TRUE/io) {
add_default($nl, 'clubb_l_use_C7_Richardson');
add_default($nl, 'clubb_l_use_C11_Richardson');
add_default($nl, 'clubb_l_use_cloud_cover');
- add_default($nl, 'clubb_l_use_ice_latent');
add_default($nl, 'clubb_l_use_thvm_in_bv_freq');
add_default($nl, 'clubb_l_vert_avg_closure');
- add_default($nl, 'clubb_l_diag_Lscale_from_tau');
add_default($nl, 'clubb_l_damp_wp2_using_em');
+ add_default($nl, 'clubb_l_godunov_upwind_wpxp_ta');
+ add_default($nl, 'clubb_l_godunov_upwind_xpyp_ta');
+ add_default($nl, 'clubb_l_use_shear_Richardson');
+ add_default($nl, 'clubb_l_use_tke_in_wp3_pr_turb_term');
+ add_default($nl, 'clubb_l_use_tke_in_wp2_wp3_K_dfsn');
+ add_default($nl, 'clubb_l_smooth_Heaviside_tau_wpxp');
+ add_default($nl, 'clubb_l_do_expldiff_rtm_thlm');
+
+ #CLUBB+MF options
add_default($nl, 'do_clubb_mf');
add_default($nl, 'do_clubb_mf_diag');
add_default($nl, 'clubb_mf_L0');
diff --git a/bld/namelist_files/namelist_defaults_cam.xml b/bld/namelist_files/namelist_defaults_cam.xml
index 8366f2d9b7..b5ce306c32 100644
--- a/bld/namelist_files/namelist_defaults_cam.xml
+++ b/bld/namelist_files/namelist_defaults_cam.xml
@@ -1823,9 +1823,6 @@
.false.
.false.
- .false.
- .true.
- .false.
.false.
1.0D0
@@ -1836,12 +1833,17 @@
150.0D0
75.0D0
+
+ .false.
+
1.0
1.0
1.0
1.0
1.3
5.2
+ 0.3
+ 0.3
4.0
6.0
1.0
@@ -1856,52 +1858,73 @@
0.35D0
2.2D0
1.6D0
+ 0.4
+ 0.75
+ 0.125
+ 5.0
+ 1.25
0.25
20.0
0.5
0.3
.false.
60.0
- 0.308
- 0.280
- 0.270
- 0.32
- 2.4
- 0.04
- 0.1
- 1.0D0
- 0.0
- 4.5
- 2.0
- 0.0
- 8.0D-6
- 25.0D-6
- 238.15D0
-
+ 1.0
+ 0.1
+ 0.02
+ 0.1
+ 0.2
+ 0.2
+ 0.0
+ 0.0
+ 0.308
+ 0.280
+ 0.270
+ 0.32
+ 2.4
+ 0.04
+ 0.1
+ 1.0D0
+ 0.0
+ 4.5
+ 2.0
+ 0.0
+ 8.0D-6
+ 25.0D-6
+ 238.15D0
+ 1
+
+ .true.
.false.
.true.
.false.
.false.
- .false.
- .false.
+ .true.
+ .true.
.false.
.false.
.true.
.false.
.false.
+ .false.
+ .false.
.true.
.false.
.false.
+ .false.
.true.
- .false.
.false.
.true.
- .false.
.false.
+ .false.
+ .false.
+ .false.
+ .false.
+ .false.
+ .false.
.true.
-
0.2
0.2
0.2
@@ -1911,7 +1934,7 @@
0.02
0.5
0.5
- 1.0
+ 0.5
0.25
20.0
2.0
@@ -1924,7 +1947,7 @@
1.5
4.0
10.0
- 4.0
+ 4.0
0.0
.true.
@@ -1943,7 +1966,7 @@
.true.
.false.
.false.
- .true.
+ .false.
.false.
@@ -2250,6 +2273,17 @@
0.0
0.0
+ .true.
+ .false.
+ .false.
+ .true.
+ .true.
+ .false.
+ .true.
+ .true.
+ .true.
+ .false.
+ .true.
NONE
diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml
index 8a89f7a76f..72ef68e5d5 100644
--- a/bld/namelist_files/namelist_definition.xml
+++ b/bld/namelist_files/namelist_definition.xml
@@ -3273,6 +3273,61 @@ Intercept of linear equation that calculates precribed in-cloud ice mixing ratio
Intercept of linear equation that calculates precribed in-cloud ice concentration ratio [N_i'^2] / [N_i]^2 [-]
+
+Enables importance sampling for SILHS subcolumns
+
+
+
+Enables calculation of Lscale_vert_avg, used to generate SILHS samples.
+
+
+
+Enables straight Monte Carlo sampling, this overrides l_lh_importance_sampling.
+
+
+
+Enables the "new" SILHS importance sampling scheme with prescribed probabilities. Requires l_lh_importance_sampling.
+
+
+
+Determine starting SILHS first sampling level (k_lh_start) based on maximum within-cloud rcm. If false, and if l_random_k_lh_start is also false, then the SILHS first sampling level is the column maximum of liquid cloud water.
+
+
+
+Determine starting SILHS first sampling level (k_lh_start) based on random choice. Overrides l_rcm_in_cloud_k_lh_start if true.
+
+
+
+Assumption of maximum vertical overlap when grid-box rcm exceeds cloud threshold.
+
+
+
+Produces "instantaneous" variance-covariance microphysical source terms, ignoring discretization effects.
+
+
+
+Limit SILHS sample point weights for stability.
+
+
+
+Prescribe variance fractions.
+
+
+
+Scale sample point weights to sum to num_samples (the "ratio estimate").
+
+
-
-Explicit diffusion on temperature and moisture when CLUBB is on
-Default: .false.
-
+
+Option for the placement of the call to CLUBB's PDF closure. The options include: ipdf_pre_advance_fields (1) calls the PDF closure before advancing prognostic fields. ipdf_post_advance_fields (2) calls after advancing prognostic fields, and ipdf_pre_post_advance_fields (3) calls both before and after advancing prognostic fields.
+Default: 1
+
+
Limiting value of C1 when skewness of w (vertical velocity) is small in
@@ -3583,6 +3639,16 @@ C4 coefficient in the wp2 return-to-isotropy term. A higher value of C4
tends wp2 more towards the value of subgrid TKE.
+
+Coefficient in the wp2 (variance of vertical velocity) pressure terms opposing shear production.
+
+
+
+Coefficient in the wp2 (variance of vertical velocity) pressure terms opposing buoyancy production.
+
+
Low Skewness in C7 Skw. Function
@@ -3605,11 +3671,45 @@ the damping of CLUBB's wp3 when skewness of w (vertical velocity) is large in
magnitude.
+
+Coefficient in the pressure-turbulence term of CLUBB's wp3 predictive equation.
+
+
+
+Coefficient of Kh_zm (diffusivity on momentum grid levels) in the wp2 (variance
+of vertical velocity) predictive equation.
+Default: 0.75
+
+
+
+Coefficient of Kh_zm (diffusivity on momentum grid levels) in the scalar
+variance predictive equations (e.g. rtp2, variance of total water).
+Default: 0.125
+
+
+
+Constant in the diffusivity term in the scalar variance predictive equations
+(e.g. rtp2, variance of total water).
+Default: 5.0
+
+
+
+Coefficient of Kh_zt (diffusivity on thermodynamic grid levels) in the wp3
+(third-order moment of vertical velocity) predictive equation.
+Default: 1.25
+
+
Coefficient of Kh_zm (diffusivity on momentum grid levels) in the up2 (variance
of the west-east wind component) and vp2 (variance of the south-north wind
component) predictive equations.
+Default: 0.25
CLUBB tunable parameter - Lscale threshold: damp C6 and C7 (units: m)
+Default: 60.0
+
+Coefficient of inverse tau term contributed by background constant value (units: none)
+Default: 1.0
+
+
+
+Coefficient of inverse tau term contributed by surface log law (units: none)
+Default: 0.1
+
+
+
+Coefficient of inverse tau term contributed by vertical wind shear (units: none)
+Default: 0.02
+
+
+
+Coefficient of inverse tau term contributed by Brunt Vaisala frequency (units: none)
+Default: 0.1
+
+
+
+Coefficient of inverse tau term contributed by Brunt Vaisala frequency but for wp3_wp2 (units: none)
+Default: 0.2
+
+
+
+Coefficient of inverse tau term contributed by Brunt Vaisala frequency but for xp2_wpxp (units: none)
+Default: 0.2
+
+
+
+Coefficient of inverse tau term contributed by Brunt Vaisala frequency but for xm_wpxp (units: none)
+Default: 0.0
+
+
+
+Coefficient of inverse tau term contributed by Brunt Vaisala frequency but for wp3 (units: none)
+Default: 0.0
+
+
-Low Skw.: gamma coef. Skw. Fnct.
+Low Skewness in gamma coefficient Skewness Function (units: none)
+Default: Changes depending on grid and physics options
-
Factor used in calculating the surface values of up2 (variance of the u wind
component) and vp2 (variance of the v wind component). Increasing
-clubb_up2_vp2_factor increases the values of up2 and vp2 at the surface.
+clubb_up2_sfc_coef increases the values of up2 and vp2 at the surface.
+Default: 2.0
Coefficient for gustiness near ground.
+Default: 0.0
+
+Flag to use shear in the calculation of Richardson number.
+Default: .false.
+
+
Flag to allow cloud fraction and mean cloud water at adjacent vertical grid
levels influence the amount of cloudiness and amount of cloud water in a
grid box.
-
-
-
-Include the effects of ice latent heating in turbulence terms
-Default: .false.
+Default: .true.
Flag to use mean theta-v in the calculation of Brunt-Vaisala frequency.
+Default: .false.
Flag that, when it is enabled, automatically enables CLUBB's
l_trapezoidal_rule_zt, l_trapezoidal_rule_zm, and l_call_pdf_closure_twice.
+Default: .true.
+
+
+
+Flag to apply Linear Multistep Method (LMM) stepping in CLUBB.
+Default: .false.
+
+
+
+Flag to run CLUBB with E3SM settings.
+Default: .true.
+
+Flag to use Total Kenetic Energy (TKE) in eddy diffusion for wp2 and wp3.
+Default: .false.
+
+
+
+Flag to use Total Kenetic Energy (TKE) formulation for wp3 pr_turb (turbulent
+production) term.
+Default: .false.
+
+
+
+Flag to use smooth Heaviside 'Peskin' in computation of invrs_tau.
+Default: .false.
+
+
+
+This flag determines whether we want to use an upwind differencing approximation
+rather than a centered differencing for turbulent advection terms. It affects
+wpxp only.
+Default: .false.
+
+
+
+This flag determines whether we want to use an upwind differencing approximation
+rather than a centered differencing for turbulent advection terms. It affects
+xpyp only.
+Default: .false.
+
+
Flag to use a dissipation formula of -(2/3)*em/tau_zm, as in Bougeault (1981),
in the wp2 (variance of vertical velocity) predictive equation.
+Default: .false.
+
+
+
+Explicit diffusion on temperature and moisture by CLUBB, in addition to CLUBB's
+normal prognostic equations for rtm and thlm.
+Default: .false.
+
+SCAM to calculate or read tendencies from a global ana/dycore
+Default: FALSE
+
+
+
+Use 1st order upwind for ana tendencies (instead of 2nd order space centered)
+Default: FALSE
+
+
+
+Use scam state as center column in stencil for ana adv tendencies
+Default: FALSE
+
+
+
+Use scam state as center column in stencil for ana adv tendencies
+Default: FALSE
+
+
+
+Use scam state as center column in stencil for ana adv tendencies
+Default: FALSE
+
+
+
+Use scam state as center column in stencil for ana adv tendencies
+Default: FALSE
+
+
+
+Force scam to use tendencies directly from dycore or ana (not recalculated)
+Default: FALSE
+
+
+
+Force scam to use omega directly from dycore or ana (not recalculated)
+Default: FALSE
+
+
+
+Interpolate ana fields to constant pressure surfaces
+Default: FALSE
+
+
+
+
+template for analysis forcing dataset.
+Default: set by build-namelist.
+
+
+
+templatefull path for analysis forcing dataset.
+Default: set by build-namelist.
+
+
+
+Force scam to compute large-scale forcing from renalysis or 3D model output
+Default: FALSE
+
+
0):
+ freqw=freqc*h0frq
+
+ freqs = str( freqw )
+ freqs = freqs.strip()+'S'
+ print('write interval='+freqs)
+
+ fl = sorted( glob.glob( dir +'/*cam.h0*') )
+ nf = len( fl )
+
+ #fl=fl[0:10]
+ ird=0
+ for f in fl:
+ print(f)
+ try:
+ a=xr.open_dataset( f )
+ print("Successfully opened"+f)
+ nx=a.lon.size
+ ny=a.lat.size
+ nl=a.lev.size
+ nli=a.ilev.size
+ nt=a.time.size
+
+ aa=a[fld] #.isel(time=0)
+ print(aa.dims)
+
+ if ('lon' in aa.dims) and ('lat' in aa.dims) and ('lev' not in aa.dims) and ('ilev' not in aa.dims):
+ print(fld+' is a surface var ' )
+ varType = 'surface'
+ if ('lon' in aa.dims) and ('lat' in aa.dims) and ('lev' in aa.dims):
+ print(fld+' is a profile var ' )
+ varType = 'profile'
+ if ('lon' in aa.dims) and ('lat' in aa.dims) and ('ilev' in aa.dims):
+ print(fld+' is an Iprofile var ' )
+ varType = 'iprofile'
+
+ if (ird == 0):
+ #Cook up time array
+ timeData=a.time.data
+ #interval=a.time[1]-a.time[0]
+ #intersec= ( interval.data.astype(int) / 10**9 ) # this is in NANOseconds
+ bigTime = xr.cftime_range( timeData[0] , periods=nt*nf, freq=freqs , calendar="noleap" )
+ if varType == 'surface':
+ dummy=np.zeros( nt*nf )
+ dummy=dummy.reshape( nt*nf ,1,1)
+ cu=xr.DataArray( dummy , coords=[bigTime,a.lat,a.lon], dims=['time','lat','lon'] , name=fld )
+ if varType == 'profile':
+ dummy=np.zeros( nt*nf*nl )
+ dummy=dummy.reshape( nt*nf,nl ,1,1)
+ cu=xr.DataArray( dummy , coords=[bigTime,a.lev,a.lat,a.lon], dims=['time','lev','lat','lon'], name=fld )
+ if varType == 'iprofile':
+ dummy=np.zeros( nt*nf*nli )
+ dummy=dummy.reshape( nt*nf,nli ,1,1)
+ cu=xr.DataArray( dummy , coords=[bigTime,a.ilev,a.lat,a.lon], dims=['time','ilev','lat','lon'], name=fld )
+
+ print('Prepped XARRAY for data with time')
+
+ if varType == 'surface':
+ #cu.values[ ird*nt :(ird+1)*nt-1, 0, 0] = aa[:,0,0] # why the heck is this wrong?????!!!!!
+ cu.values[ ird*nt :(ird+1)*nt , 0, 0] = aa[:,0,0]
+ if varType == 'profile' or varType == 'iprofile':
+ cu.values[ ird*nt :(ird+1)*nt , : , 0, 0] = aa[:,:,0,0]
+
+ ird=ird+1
+
+ except ValueError:
+ print('******** VALUE ERROR *************')
+ print('File \n'+f+'\n probably no good')
+ if (ird == 0):
+ exit('First dataset not valid')
+ else:
+ ird=ird+1
+
+ cu.to_netcdf('/project/amp/juliob/scam/'+xp+'_'+fld+'.nc')
+
+ return cu
+
+
+ def ncmerge(self):
+ import numpy as np
+ import xarray as xr
+ import glob
+ import txtutil as tx
+
+
+
+ xp=self.case
+ dir=self.archdir
+
+ fl = sorted( glob.glob( dir +'/*cam.h0*') )
+ nf = len( fl )
+
+ flt = fl[0:10]
+
+ #ds=xr.open_mfdataset( flt , concat_dim='time')
+ #print("merged nc files")
+ #ds.to_netcdf(path = dir +'/merged_h0.nc')
+ #exit()
+
+ ird=0
+ for f in flt:
+ print(f)
+ try:
+ a=xr.open_dataset( f )
+ print("Successfully opened"+f)
+
+ if (ird == 0):
+ b=a
+ elif (ird>0):
+ b.merge(a,join='inner' )
+ b.to_netcdf( 'testmerge.nc' )
+
+ ird=ird+1
+
+ except ValueError:
+ print('******** VALUE ERROR *************')
+ print('File \n'+f+'\n probably no good')
+ if (ird == 0):
+ exit('First dataset not valid')
+ else:
+ ird=ird+1
+
+ return
diff --git a/myPythonTools/scam_case.py b/myPythonTools/scam_case.py
new file mode 100644
index 0000000000..af5eb97ac6
--- /dev/null
+++ b/myPythonTools/scam_case.py
@@ -0,0 +1,342 @@
+class scam_case:
+ #import sys
+ #import getopt as go
+ #import os
+
+ def __init__(self):
+ import os
+
+ self.scmlat=36.6
+ self.scmlon=270.
+ self.nlev=58
+ self.name="base_case"
+ self.tag="case_tag"
+ self.startdate=[2010,4,1]
+ self.atm_ncpl=192
+ self.mfilt=1
+ self.nsteps=31*self.atm_ncpl
+ self.coupler="nuopc"
+ self.compiler="intel"
+ self.basecase="base_case"
+ self.isbasecase=True
+ self.NameByBuild=False
+ self.cime_output_root="dir"
+ self.ensemble_root="none"
+ self.project="P93300642" # Needed for Cheyenne
+
+ host=os.environ['HOST']
+
+ if ('izumi' in host):
+ self.machine="izumi"
+ elif ('cheyenne' in host):
+ self.machine="cheyenne"
+
+ def base_case(self):
+ import subprocess as sp
+ import os
+ import pickle
+
+ user=os.environ['USER']
+
+ y = self.startdate[0]
+ m = self.startdate[1]
+ d = self.startdate[2]
+ lat = self.scmlat
+ lon = self.scmlon
+ lev = self.nlev
+ tag = self.tag
+
+
+ case_day = str(d).zfill(2)
+ case_mon = str(m).zfill(2)
+ case_yr = str(y).zfill(4)
+
+ case_date = case_yr+case_mon+case_day
+ case_sdate = case_yr+"-"+case_mon+"-"+case_day
+
+ case_lat = str(abs(lat)).zfill(4)
+ if lat < 0:
+ case_lat = case_lat+'S'
+ elif lat > 0:
+ case_lat = case_lat+'N'
+ if lon < 0:
+ lon = lon+360.
+
+ case_lon = str(abs(lon)).zfill(5)+'E'
+ case_lev = 'L'+str(abs(lev))
+
+ lonstr = str(lon)
+ latstr = str(lat)
+ levstr = str(lev)
+
+ NstepStr = str( self.nsteps )
+ NcplStr = str( self.atm_ncpl )
+ CplStr = str( self.coupler )
+ CompilerStr = str( self.compiler )
+ MachStr = str( self.machine )
+ ProjStr = str( self.project )
+
+ if (self.NameByBuild == True):
+ case_tag = tag+'_'+MachStr+'_'+CompilerStr+'_'+CplStr
+ else:
+ case_tag = tag+'_'+case_lev+'_'+case_lon+'_'+case_lat+'_'+case_yr+'-'+case_mon+'-'+case_day
+
+ if ( self.coupler=="nuopc"):
+ COMPSET="FSCAM"
+ elif ( self.coupler=="mct"):
+ COMPSET="2000_CAM60%SCAM_CLM50%SP_CICE5%PRES_DOCN%DOM_SROF_SGLC_SWAV"
+
+
+ #This is the actual case name for the 'base' case
+ print(case_tag)
+ self.name=case_tag
+ self.basecase=case_tag
+ self.isbasecase=True
+
+ #-----------------------------------------------------
+ # This script assumes that you will create new cases in
+ # in a directory '../../cases' from
+ # ${CESMROOT}/cime/scripts
+ #-----------------------------------------------------
+ cmd1="mkdir -p ../../cases/"+case_tag
+
+ if (MachStr == 'cheyenne'):
+ cmd2="./create_newcase --case ../../cases/"+case_tag+ " --compset "+ COMPSET + " --res T42_T42 --driver " + CplStr + " --user-mods-dir ../../cime_config/usermods_dirs/scam_STUB --walltime 01:00:00 --mach " + MachStr + " --pecount 1 --compiler "+ CompilerStr + " --project "+ ProjStr + " --run-unsupported"
+ else:
+ cmd2="./create_newcase --case ../../cases/"+case_tag+ " --compset "+ COMPSET + " --res T42_T42 --driver " + CplStr + " --user-mods-dir ../../cime_config/usermods_dirs/scam_STUB --walltime 01:00:00 --mach " + MachStr + " --pecount 1 --compiler "+ CompilerStr + " --run-unsupported"
+
+ cd0 = 'cd ../../cases/'+case_tag +';'
+
+ #sp.run('date')
+ sp.run(cmd2,shell=True)
+
+ cmd = ( "./case.setup" )
+ sp.run(cd0 + cmd , shell=True )
+
+ #cmd = ( "cp ../../cime_config/usermods_dirs/scam_STUB/scripts/STUB_iop.nc ./")
+ #sp.run(cd0 + cmd , shell=True )
+
+ cmd = ( "cp ../../myPythonTools/STUB_iop.nc ./")
+ sp.run(cd0 + cmd , shell=True )
+ cmd = ( "cp ../../myPythonTools/user_nl_cice ./")
+ sp.run(cd0 + cmd , shell=True )
+
+ if (self.nlev==58) and (self.machine=='cheyenne'):
+ cmd = ( "cp ../../myPythonTools/user_nl_cam_L58_cheyenne ./user_nl_cam")
+ elif (self.nlev==93) and (self.machine=='cheyenne'):
+ cmd = ( "cp ../../myPythonTools/user_nl_cam_L93_cheyenne ./user_nl_cam")
+ else:
+ cmd = ( "cp ../../myPythonTools/user_nl_cam ./")
+
+ sp.run(cd0 + cmd , shell=True )
+
+ cmd = ( "./xmlchange DOUT_S_ROOT='/project/amp/"+user+"/scam/archive/"+case_tag+"'" + ";" +
+ "./xmlchange CAM_CONFIG_OPTS='-dyn eul -scam -phys cam_dev -nlev "+levstr+"'"+ ";" +
+ "./xmlchange ATM_NCPL="+NcplStr+";"+
+ "./xmlchange STOP_N="+NstepStr+";"+
+ "./xmlchange START_TOD=00000"+";"+
+ "./xmlchange STOP_OPTION=nsteps"+";"+
+ "./xmlchange PTS_LAT="+latstr+";"+
+ "./xmlchange PTS_LON="+lonstr+";"+
+ "./xmlchange RUN_STARTDATE="+case_sdate+";"
+ )
+ sp.run(cd0 + cmd , shell=True )
+
+
+ cmd = (
+ "ncap2 --overwrite -s bdate="+case_date+" STUB_iop.nc STUB_iop.nc"+";"+
+ "ncap2 --overwrite -s lat[lat]="+latstr+" STUB_iop.nc STUB_iop.nc"+";"+
+ "ncap2 --overwrite -s lon[lon]="+lonstr+" STUB_iop.nc STUB_iop.nc"+";"
+ )
+ sp.run(cd0 + cmd , shell=True )
+
+
+ print("Machine = "+self.machine)
+ print("Compiler = "+self.compiler)
+ print("created and setup case=")
+ print(" ../../cases/"+case_tag )
+ print("Should be ready to build and submit")
+ fname = '../../cases/'+case_tag+'/'+'env_build.xml'
+
+ # find CIME_OUTPUT_ROOT
+ fob=open(fname,"r")
+ linin = fob.readlines()
+ for line in linin:
+ if (line.find('entry id="CIME_OUTPUT_ROOT"') !=-1):
+ linspl1 = line.split('value=')
+ linspl2 = linspl1[1].split('"')
+ self.cime_output_root = linspl2[1]
+ fob.close()
+
+ # write 'self' to pickle file in casedir
+ fname = '../../cases/'+case_tag+'/'+'BaseCaseSelf.pkL'
+ with open( fname, 'wb') as fob:
+ pickle.dump( self, fob )
+ fob.close()
+
+
+ def spawn_case(self,basecase):
+ #---------------------------------------
+ # This function is still under development
+ #---------------------------------------
+ import subprocess as sp
+ import os
+ import pickle
+ import txtutil as tx
+
+ user=os.environ['USER']
+
+ fname = '../../cases/'+basecase+'/'+'BaseCaseSelf.pkL'
+ with open( fname, 'rb') as fob:
+ base=pickle.load( fob )
+ fob.close()
+
+ #---------------------------------
+ # These sould be inherited from
+ # base case or hardwried here
+ #---------------------------------
+ self.basecase = base.name
+ self.isbasecase = False
+ self.nlev = base.nlev
+ self.coupler = base.coupler
+ self.compiler = base.compiler
+ self.atm_ncpl = base.atm_ncpl
+ self.cime_output_root = base.cime_output_root
+ lev = base.nlev
+
+ #-----------------------------------
+ # These are specified in scam_drv.py
+ # or scam_ens.py
+ #-----------------------------------
+ y = self.startdate[0]
+ m = self.startdate[1]
+ d = self.startdate[2]
+ lat = self.scmlat
+ lon = self.scmlon
+ tag = self.tag
+
+
+ case_day = str(d).zfill(2)
+ case_mon = str(m).zfill(2)
+ case_yr = str(y).zfill(4)
+
+ case_date = case_yr+case_mon+case_day
+ case_sdate = case_yr+"-"+case_mon+"-"+case_day
+
+ case_lat = str(abs(lat)).zfill(4)
+ if lat < 0:
+ case_lat = case_lat+'S'
+ elif lat > 0:
+ case_lat = case_lat+'N'
+ if lon < 0:
+ lon = lon+360.
+
+ case_lon = str(abs(lon)).zfill(5)+'E'
+ case_lev = 'L'+str(abs(lev))
+
+ lonstr = str(lon)
+ latstr = str(lat)
+ levstr = str(lev)
+
+ case_tag = tag+'_'+case_lev+'_'+case_lon+'_'+case_lat+'_'+case_yr+'-'+case_mon+'-'+case_day
+ self.name = case_tag
+
+ ensemble_root = self.cime_output_root + '/' + self.basecase +'_ENS'
+ self.ensemble_root = ensemble_root
+
+
+ # --------------
+ # Clean before making new directory
+ #----------------
+ cmd = ("rm -rf "+ ensemble_root + "/"+case_tag)
+ sp.run( cmd , shell=True )
+
+ cmd = ("mkdir -p "+ ensemble_root + "/"+case_tag+"/bld")
+ sp.run( cmd , shell=True )
+
+ cmd = ("cp -r "+ self.cime_output_root + "/" + base.name +"/run" + " "
+ + ensemble_root + "/"+case_tag+"/run")
+ sp.run( cmd , shell=True )
+
+ cmd = ("cp "+ self.cime_output_root + "/" + base.name +"/bld/cesm.exe" + " "
+ + ensemble_root + "/"+case_tag+"/bld/")
+ sp.run( cmd , shell=True )
+
+ cmd = ( "cp ../../myPythonTools/STUB_iop.nc" + " "
+ + ensemble_root + "/"+case_tag+"/run/")
+ sp.run(cmd , shell=True )
+ cmd = ( "cp ../../myPythonTools/ens_run.sh" + " "
+ + ensemble_root + "/"+case_tag+"/run/")
+ sp.run(cmd , shell=True )
+
+ cd0 ="cd "+ ensemble_root + "/" + case_tag +"/run ;"
+
+ cmd = (
+ "ncap2 --overwrite -s bdate="+case_date+" STUB_iop.nc STUB_iop.nc"+";"+
+ "ncap2 --overwrite -s lat[lat]="+latstr+" STUB_iop.nc STUB_iop.nc"+";"+
+ "ncap2 --overwrite -s lon[lon]="+lonstr+" STUB_iop.nc STUB_iop.nc"+";"
+ )
+ sp.run(cd0 + cmd , shell=True )
+
+ fili= ensemble_root + "/" + case_tag +"/run/atm_in"
+ tx.nmled(fili,'iopfile','"STUB_iop.nc"')
+ # Set history to make one file per day
+ tx.nmled(fili,'mfilt',str(self.mfilt) )
+
+ if (base.coupler=='nuopc'):
+ fili= ensemble_root + "/" + case_tag +"/run/nuopc.runconfig"
+ tx.nmled(fili,'case_name',case_tag)
+ tx.nmled(fili,'start_ymd',case_date)
+ tx.nmled(fili,'scol_lat',latstr)
+ tx.nmled(fili,'scol_lon',lonstr)
+
+ if (base.coupler=='mct'):
+ fili= ensemble_root + "/" + case_tag +"/run/drv_in"
+ tx.nmled(fili,'case_name',case_tag)
+ tx.nmled(fili,'start_ymd',case_date)
+ tx.nmled(fili,'scmlat',latstr)
+ tx.nmled(fili,'scmlon',lonstr)
+
+ def ensemble_member_run(self):
+ import subprocess as sp
+ import os
+
+ cd0 ="cd "+ self.ensemble_root + "/" + self.name +"/run ;"
+ cmd ="/usr/local/torque/bin/qsub ens_run.sh"
+
+ sp.run(cd0 + cmd , shell=True )
+
+ def unpickle_base(self,basecase):
+ #---------------------------------------
+ # This function is still under development
+ #---------------------------------------
+ import pickle
+
+ fname = '../../cases/'+basecase+'/'+'BaseCaseSelf.pkL'
+ with open( fname, 'rb') as fob:
+ base=pickle.load( fob )
+ fob.close()
+
+ return base
+
+ def changeTag(self,newtag):
+ #---------------------------------------
+ # This function is still under development
+ #---------------------------------------
+
+ self.tag = newtag
+
+ def changeLon(self,newlon):
+ #---------------------------------------
+ # This function is still under development
+ #---------------------------------------
+
+ self.scmlon = newlon
+
+ def changeLat(self,newlat):
+ #---------------------------------------
+ # This function is still under development
+ #---------------------------------------
+
+ self.scmlat = newlat
+
diff --git a/myPythonTools/scam_drv.py b/myPythonTools/scam_drv.py
new file mode 100755
index 0000000000..904a9fd3d4
--- /dev/null
+++ b/myPythonTools/scam_drv.py
@@ -0,0 +1,69 @@
+#!/usr/bin/env python
+
+import getopt as go
+import scam_case as scm
+import sys
+import os
+
+argv=sys.argv
+
+case=scm.scam_case()
+user=os.environ['USER']
+spawncase=False
+
+print(case.scmlon)
+
+try:
+ opts, args = go.getopt( argv[1:], "i:j:y:m:d:t:l:x:n:q:c:S:M:N",
+ ["lon=","lat=","year=","month=","day=","tag=","nlev=","coupler=","nsteps="
+ ,"atm-ncpl=","compiler=","spawn=","machine=","NameByBuild="] )
+except:
+ print( "something is wrong")
+ exit()
+
+date=case.startdate
+print(date)
+
+for opt, arg in opts:
+ if opt in ("-i","--lon"):
+ case.scmlon = float(arg)
+ elif opt in ("-j","--lat"):
+ case.scmlat = float(arg)
+ elif opt in ("-y","--year"):
+ case.startdate[0] = int(arg)
+ elif opt in ("-m","--month"):
+ case.startdate[1] = int(arg)
+ elif opt in ("-d","--day"):
+ case.startdate[2] = int(arg)
+ elif opt in ("-l","--nlev"):
+ case.nlev = int(arg)
+ elif opt in ("-t","--tag"):
+ case.tag = arg
+ elif opt in ("-x","--coupler"):
+ case.coupler = arg
+ elif opt in ("-n","--nsteps"):
+ case.nsteps = int(arg)
+ elif opt in ("-q","--atm-ncpl"):
+ case.atm_ncpl = int(arg)
+ elif opt in ("-c","--compiler"):
+ case.compiler = arg
+ elif opt in ("-M","--machine"):
+ case.machine = arg
+ elif opt in ("-S","--spawn"):
+ basecase = arg
+ spawncase = True
+ elif opt in ("-N","--NameByBuild"):
+ case.NameByBuild=True
+
+
+date=case.startdate
+print(date)
+
+if (spawncase == False):
+ case.base_case()
+else:
+ case.spawn_case(basecase)
+
+#fname = '../../cases/'+case_tag+'/'+'CaseInst.pysave'
+#with open( fname, 'wb') as fob:
+# pickle.dump( case, fob )
diff --git a/myPythonTools/scam_ens.py b/myPythonTools/scam_ens.py
new file mode 100755
index 0000000000..be5502a8f1
--- /dev/null
+++ b/myPythonTools/scam_ens.py
@@ -0,0 +1,49 @@
+#!/usr/bin/env python
+
+import scam_case as scm
+import numpy as np
+
+#basecase = 'nCTOPb3_L58_080.0E_30.0N_2010-07-01'
+basecase = 'NCPL96_izumi_intel_nuopc'
+base=scm.scam_case()
+base=base.unpickle_base(basecase)
+
+print(base.__dict__)
+
+#exit()
+
+starting_ens=21
+lats = 32.+ np.arange(2)
+lons = 85.+ np.arange(10)
+
+Lons,Lats = np.meshgrid(lons,lats)
+
+
+dims=Lons.shape
+
+
+x=[]
+n=0
+for j in range( dims[1] ):
+ for i in range( dims[0]):
+ x.append( scm.scam_case() )
+ n=n+1
+
+Lons=Lons.reshape( dims[0]*dims[1] )
+Lats=Lats.reshape( dims[0]*dims[1] )
+
+N=n
+for n in range(N):
+ nens=n+starting_ens
+ ee = 'x_E'+str(nens).zfill(2)
+ x[n].changeTag( ee )
+ x[n].changeLon( Lons[n] )
+ x[n].changeLat( Lats[n] )
+ x[n].startdate = [2010,6,1] #base.startdate
+
+for n in range(N):
+ x[n].spawn_case(basecase)
+
+for n in range(N):
+ x[n].ensemble_member_run()
+
diff --git a/myPythonTools/testh0.py b/myPythonTools/testh0.py
new file mode 100644
index 0000000000..67394bf105
--- /dev/null
+++ b/myPythonTools/testh0.py
@@ -0,0 +1,20 @@
+import h0scam as h0
+
+#----------------------------------
+# Run from python prompt like this:
+# exec(open("./testh0.py").read())
+
+#dir='/scratch/cluster/juliob/nCTOPb3_L58_080.0E_30.0N_2010-07-01_ENS/x_E03_L58_081.0E_33.0N_2010-07-01/run/'
+#xp='x_E03_L58_081.0E_33.0N_2010-07-01'
+
+dir='/scratch/cluster/juliob/NCPL96_izumi_intel_nuopc_ENS/x_E01_L58_085.0E_32.0N_2010-04-01/run/'
+xp='x_E01_L58_085.0E_32.0N_2010-04-01'
+
+#dir='/scratch/cluster/juliob/nCTOPb2_L58_080.0E_30.0N_2010-07-01/run/'
+#xp='nCTOPb2_L58_080.0E_30.0N_2010-07-01'
+
+x = h0.h0scam(xp=xp,dir=dir)
+
+#x.ncmerge()
+
+cu=x.curtain('rtm_mfl')
diff --git a/myPythonTools/txtutil.py b/myPythonTools/txtutil.py
new file mode 100644
index 0000000000..415ade01cd
--- /dev/null
+++ b/myPythonTools/txtutil.py
@@ -0,0 +1,43 @@
+import subprocess as sp
+
+# which namelist parameter to change
+
+def nmled(fili,parm,valu):
+ #which file
+ #fili = "atm_in"
+ filo = fili+"_edit"
+
+
+ fin = open( fili ,"r")
+ fex = open( filo ,"w")
+ linin = fin.readlines()
+ for line in linin:
+ spl = line.split("=")
+ #if (line.find("zmconv_ke") !=-1):
+ if (spl[0].strip() == parm):
+ print(spl)
+ spl[1] = " "+ valu +" \n"
+ print(spl)
+ line="=".join(spl)
+
+ fex.write(line)
+
+ fin.close()
+ fex.close()
+
+ cmd = "mv "+filo+" "+fili
+ sp.run( cmd, shell=True)
+
+def nmlread(fili,parm):
+
+ valu=-99999
+ fin = open( fili ,"r")
+ linin = fin.readlines()
+ for line in linin:
+ spl = line.split("=")
+ if (spl[0].strip() == parm):
+ valu=spl[1]
+
+ fin.close()
+
+ return valu
diff --git a/myPythonTools/user_nl_cam b/myPythonTools/user_nl_cam
new file mode 100644
index 0000000000..eac76629f5
--- /dev/null
+++ b/myPythonTools/user_nl_cam
@@ -0,0 +1,93 @@
+!scmlon=$PTS_LON
+!scmlat=$PTS_LAT
+iopfile="$CASEROOT/STUB_iop.nc"
+ncdata="/project/amp/juliob/scam/inputdata/SCAM_IC_288x192_L58_48_BL10.nc"
+
+
+bnd_topo="/fs/cgd/csm/inputdata/atm/cam/topo/fv_0.9x1.25_nc3000_Nsw042_Nrs008_Co060_Fi001_ZR_sgh30_24km_GRNL_c170103.nc"
+
+mfilt=5760
+
+nhtfrq=1
+scm_use_obs_uv = .false.
+scm_relaxation = .false.
+scm_relax_fincl = 'T', 'bc_a1', 'bc_a4', 'dst_a1', 'dst_a2', 'dst_a3', 'ncl_a1', 'ncl_a2',
+ 'ncl_a3', 'num_a1', 'num_a2', 'num_a3',
+ 'num_a4', 'pom_a1', 'pom_a4', 'so4_a1', 'so4_a2', 'so4_a3', 'soa_a1', 'soa_a2'
+scm_relax_bot_p = 105000.
+scm_relax_top_p = 200.
+scm_relax_linear = .true.
+scm_relax_tau_bot_sec = 864000.
+scm_relax_tau_top_sec = 172800.
+
+use_scm_ana_frc = .true.
+scm_ana_frc_path = '/project/amp/juliob/ERAI/f09_omega/L58/'
+scm_ana_frc_file_template = '%y/ERAI_fv09_L58.cam2.i.%y-%m-%d-%s.nc'
+
+scm_ana_x_plevels = .true.
+scm_ana_direct_omega = .true.
+scm_ana_direct_ttend = .false.
+scm_ana_t_react = .true.
+scm_ana_q_react = .true.
+scm_ana_u_react = .true.
+scm_ana_v_react = .true.
+scm_ana_upwind = .false.
+
+fincl1 = 'Target_U','Target_V','Target_T','Target_Q',
+ 'Nudge_U','Nudge_V','Nudge_T','Nudge_Q',
+ 'OMEGA_ANA','ETAD_ANA','T_ANA','Q_ANA','U_ANA','V_ANA',
+ 'UTEND_PHYSTOT', 'VTEND_PHYSTOT', 'TTEN_PHYS',
+ 'UTEND_DCONV','UTEND_SHCONV','UTEND_MACROP','UTEND_VDIFF','UTEND_RAYLEIGH',
+ 'UTEND_GWDTOT','UTEND_QBORLX','UTEND_LUNART','UTEND_IONDRG','UTEND_NDG',
+ 'VTEND_DCONV','VTEND_SHCONV','VTEND_MACROP','VTEND_VDIFF','VTEND_RAYLEIGH',
+ 'VTEND_GWDTOT','VTEND_QBORLX','VTEND_LUNART','VTEND_IONDRG','VTEND_NDG',
+ 'KVH_CLUBB','TAUARDGBETAX','TAUARDGBETAY','TAU1RDGBETAX','TAU1RDGBETAY',
+ 'UBT1RDGBETA','TAU1RDGBETAM','RVMTEND_CLUBB'
+
+
+gw_drag_file = '/fs/cgd/inputdata/inputdata/atm/waccm/gw/newmfspectra40_dc25.nc'
+use_gw_convect_dp = .false.
+use_gw_convect_sh = .false.
+use_gw_front = .false.
+scm_use_ana_iop = .true.
+cld_macmic_num_steps=3
+do_clubb_mf = .false.
+do_clubb_mf_diag = .false.
+
+&nudging_nl
+Nudge_Model = .true.
+Nudge_Path = '/project/amp/juliob/ERAI/f09_omega/L58/'
+Nudge_File_Template = '%y/ERAI_fv09_L58.cam2.i.%y-%m-%d-%s.nc'
+Nudge_Force_Opt = 0
+Nudge_TimeScale_Opt = 0
+Nudge_Times_Per_Day = 4
+Model_Times_Per_Day = 192
+Nudge_Uprof = 2
+Nudge_Ucoef = 1.0
+Nudge_Vprof = 2
+Nudge_Vcoef = 1.0
+Nudge_Tprof = 2
+Nudge_Tcoef = 1.0
+Nudge_Qprof = 2
+Nudge_Qcoef = 0.0
+Nudge_PSprof = 0
+Nudge_PScoef = 0
+Nudge_Beg_Year = 2009
+Nudge_Beg_Month = 1
+Nudge_Beg_Day = 1
+Nudge_End_Year = 2018
+Nudge_End_Month = 4
+Nudge_End_Day = 5
+Nudge_Hwin_lat0 = 0.0
+Nudge_Hwin_lon0 = 180.
+Nudge_Hwin_latWidth = 9999.0
+Nudge_Hwin_lonWidth = 9999.0
+Nudge_Hwin_latDelta = 1.0
+Nudge_Hwin_lonDelta = 1.0
+Nudge_Hwin_Invert = .false.
+Nudge_Vwin_Hindex = 0.
+Nudge_Vwin_Lindex = 0.
+Nudge_Vwin_Hdelta = 0.001
+Nudge_Vwin_Ldelta = 0.001
+Nudge_Vwin_Invert = .true.
+/
diff --git a/myPythonTools/user_nl_cam_L58_cheyenne b/myPythonTools/user_nl_cam_L58_cheyenne
new file mode 100644
index 0000000000..195097a88a
--- /dev/null
+++ b/myPythonTools/user_nl_cam_L58_cheyenne
@@ -0,0 +1,92 @@
+!scmlon=$PTS_LON
+!scmlat=$PTS_LAT
+iopfile="$CASEROOT/STUB_iop.nc"
+ncdata="/glade/p/cgd/amp/juliob/scam/inputdata/SCAM_IC_288x192_L58_48_BL10.nc"
+
+bnd_topo='/glade/p/cesmdata/cseg/inputdata/atm/cam/topo/fv_0.9x1.25_nc3000_Nsw042_Nrs008_Co060_Fi001_ZR_sgh30_24km_GRNL_c170103.nc'
+
+mfilt=5760
+
+nhtfrq=1
+scm_use_obs_uv = .false.
+scm_relaxation = .false.
+scm_relax_fincl = 'T', 'bc_a1', 'bc_a4', 'dst_a1', 'dst_a2', 'dst_a3', 'ncl_a1', 'ncl_a2',
+ 'ncl_a3', 'num_a1', 'num_a2', 'num_a3',
+ 'num_a4', 'pom_a1', 'pom_a4', 'so4_a1', 'so4_a2', 'so4_a3', 'soa_a1', 'soa_a2'
+scm_relax_bot_p = 105000.
+scm_relax_top_p = 200.
+scm_relax_linear = .true.
+scm_relax_tau_bot_sec = 864000.
+scm_relax_tau_top_sec = 172800.
+
+use_scm_ana_frc = .true.
+scm_ana_frc_path = '/glade/p/cgd/amp/juliob/ERAI/f09_omega/L58/'
+scm_ana_frc_file_template = '%y/ERAI_fv09_L58.cam2.i.%y-%m-%d-%s.nc'
+
+scm_ana_x_plevels = .true.
+scm_ana_direct_omega = .true.
+scm_ana_direct_ttend = .false.
+scm_ana_t_react = .true.
+scm_ana_q_react = .true.
+scm_ana_u_react = .true.
+scm_ana_v_react = .true.
+scm_ana_upwind = .false.
+
+fincl1 = 'Target_U','Target_V','Target_T','Target_Q',
+ 'Nudge_U','Nudge_V','Nudge_T','Nudge_Q',
+ 'OMEGA_ANA','ETAD_ANA','T_ANA','Q_ANA','U_ANA','V_ANA',
+ 'UTEND_PHYSTOT', 'VTEND_PHYSTOT', 'TTEN_PHYS',
+ 'UTEND_DCONV','UTEND_SHCONV','UTEND_MACROP','UTEND_VDIFF','UTEND_RAYLEIGH',
+ 'UTEND_GWDTOT','UTEND_QBORLX','UTEND_LUNART','UTEND_IONDRG','UTEND_NDG',
+ 'VTEND_DCONV','VTEND_SHCONV','VTEND_MACROP','VTEND_VDIFF','VTEND_RAYLEIGH',
+ 'VTEND_GWDTOT','VTEND_QBORLX','VTEND_LUNART','VTEND_IONDRG','VTEND_NDG',
+ 'KVH_CLUBB','TAUARDGBETAX','TAUARDGBETAY','TAU1RDGBETAX','TAU1RDGBETAY',
+ 'UBT1RDGBETA','TAU1RDGBETAM','RVMTEND_CLUBB'
+
+
+gw_drag_file = '/glade/p/cesmdata/cseg/inputdata/atm/waccm/gw/newmfspectra40_dc25.nc'
+use_gw_convect_dp = .true.
+use_gw_convect_sh = .false.
+use_gw_front = .false.
+scm_use_ana_iop = .true.
+cld_macmic_num_steps=3
+do_clubb_mf = .false.
+do_clubb_mf_diag = .false.
+
+&nudging_nl
+Nudge_Model = .true.
+Nudge_Path = '/glade/p/cgd/amp/juliob/ERAI/f09_omega/L58/'
+Nudge_File_Template = '%y/ERAI_fv09_L58.cam2.i.%y-%m-%d-%s.nc'
+Nudge_Force_Opt = 0
+Nudge_TimeScale_Opt = 0
+Nudge_Times_Per_Day = 4
+Model_Times_Per_Day = 192
+Nudge_Uprof = 2
+Nudge_Ucoef = 1.0
+Nudge_Vprof = 2
+Nudge_Vcoef = 1.0
+Nudge_Tprof = 2
+Nudge_Tcoef = 1.0
+Nudge_Qprof = 2
+Nudge_Qcoef = 0.0
+Nudge_PSprof = 0
+Nudge_PScoef = 0
+Nudge_Beg_Year = 2009
+Nudge_Beg_Month = 1
+Nudge_Beg_Day = 1
+Nudge_End_Year = 2018
+Nudge_End_Month = 4
+Nudge_End_Day = 5
+Nudge_Hwin_lat0 = 0.0
+Nudge_Hwin_lon0 = 180.
+Nudge_Hwin_latWidth = 9999.0
+Nudge_Hwin_lonWidth = 9999.0
+Nudge_Hwin_latDelta = 1.0
+Nudge_Hwin_lonDelta = 1.0
+Nudge_Hwin_Invert = .false.
+Nudge_Vwin_Hindex = 0.
+Nudge_Vwin_Lindex = 0.
+Nudge_Vwin_Hdelta = 0.001
+Nudge_Vwin_Ldelta = 0.001
+Nudge_Vwin_Invert = .true.
+/
diff --git a/myPythonTools/user_nl_cam_L93_cheyenne b/myPythonTools/user_nl_cam_L93_cheyenne
new file mode 100644
index 0000000000..4b94b10b8c
--- /dev/null
+++ b/myPythonTools/user_nl_cam_L93_cheyenne
@@ -0,0 +1,92 @@
+!scmlon=$PTS_LON
+!scmlat=$PTS_LAT
+iopfile="$CASEROOT/STUB_iop.nc"
+ncdata="/glade/p/cgd/amp/juliob/scam/inputdata/SCAM_IC_288x192_L93.nc"
+
+bnd_topo='/glade/p/cesmdata/cseg/inputdata/atm/cam/topo/fv_0.9x1.25_nc3000_Nsw042_Nrs008_Co060_Fi001_ZR_sgh30_24km_GRNL_c170103.nc'
+
+mfilt=5760
+
+nhtfrq=1
+scm_use_obs_uv = .false.
+scm_relaxation = .false.
+scm_relax_fincl = 'T', 'bc_a1', 'bc_a4', 'dst_a1', 'dst_a2', 'dst_a3', 'ncl_a1', 'ncl_a2',
+ 'ncl_a3', 'num_a1', 'num_a2', 'num_a3',
+ 'num_a4', 'pom_a1', 'pom_a4', 'so4_a1', 'so4_a2', 'so4_a3', 'soa_a1', 'soa_a2'
+scm_relax_bot_p = 105000.
+scm_relax_top_p = 200.
+scm_relax_linear = .true.
+scm_relax_tau_bot_sec = 864000.
+scm_relax_tau_top_sec = 172800.
+
+use_scm_ana_frc = .true.
+scm_ana_frc_path = '/glade/p/cgd/amp/juliob/ERAI/f09_omega/L93/'
+scm_ana_frc_file_template = '%y/ERAI_fv09_L93.cam2.i.%y-%m-%d-%s.nc'
+
+scm_ana_x_plevels = .true.
+scm_ana_direct_omega = .true.
+scm_ana_direct_ttend = .false.
+scm_ana_t_react = .true.
+scm_ana_q_react = .true.
+scm_ana_u_react = .true.
+scm_ana_v_react = .true.
+scm_ana_upwind = .false.
+
+fincl1 = 'Target_U','Target_V','Target_T','Target_Q',
+ 'Nudge_U','Nudge_V','Nudge_T','Nudge_Q',
+ 'OMEGA_ANA','ETAD_ANA','T_ANA','Q_ANA','U_ANA','V_ANA',
+ 'UTEND_PHYSTOT', 'VTEND_PHYSTOT', 'TTEN_PHYS',
+ 'UTEND_DCONV','UTEND_SHCONV','UTEND_MACROP','UTEND_VDIFF','UTEND_RAYLEIGH',
+ 'UTEND_GWDTOT','UTEND_QBORLX','UTEND_LUNART','UTEND_IONDRG','UTEND_NDG',
+ 'VTEND_DCONV','VTEND_SHCONV','VTEND_MACROP','VTEND_VDIFF','VTEND_RAYLEIGH',
+ 'VTEND_GWDTOT','VTEND_QBORLX','VTEND_LUNART','VTEND_IONDRG','VTEND_NDG',
+ 'KVH_CLUBB','TAUARDGBETAX','TAUARDGBETAY','TAU1RDGBETAX','TAU1RDGBETAY',
+ 'UBT1RDGBETA','TAU1RDGBETAM','RVMTEND_CLUBB'
+
+
+gw_drag_file = '/glade/p/cesmdata/cseg/inputdata/atm/waccm/gw/newmfspectra40_dc25.nc'
+use_gw_convect_dp = .true.
+use_gw_convect_sh = .false.
+use_gw_front = .false.
+scm_use_ana_iop = .true.
+cld_macmic_num_steps=3
+do_clubb_mf = .false.
+do_clubb_mf_diag = .false.
+
+&nudging_nl
+Nudge_Model = .true.
+Nudge_Path = '/glade/p/cgd/amp/juliob/ERAI/f09_omega/L93/'
+Nudge_File_Template = '%y/ERAI_fv09_L93.cam2.i.%y-%m-%d-%s.nc'
+Nudge_Force_Opt = 0
+Nudge_TimeScale_Opt = 0
+Nudge_Times_Per_Day = 4
+Model_Times_Per_Day = 192
+Nudge_Uprof = 2
+Nudge_Ucoef = 1.0
+Nudge_Vprof = 2
+Nudge_Vcoef = 1.0
+Nudge_Tprof = 2
+Nudge_Tcoef = 1.0
+Nudge_Qprof = 2
+Nudge_Qcoef = 0.0
+Nudge_PSprof = 0
+Nudge_PScoef = 0
+Nudge_Beg_Year = 2009
+Nudge_Beg_Month = 1
+Nudge_Beg_Day = 1
+Nudge_End_Year = 2018
+Nudge_End_Month = 4
+Nudge_End_Day = 5
+Nudge_Hwin_lat0 = 0.0
+Nudge_Hwin_lon0 = 180.
+Nudge_Hwin_latWidth = 9999.0
+Nudge_Hwin_lonWidth = 9999.0
+Nudge_Hwin_latDelta = 1.0
+Nudge_Hwin_lonDelta = 1.0
+Nudge_Hwin_Invert = .false.
+Nudge_Vwin_Hindex = 0.
+Nudge_Vwin_Lindex = 0.
+Nudge_Vwin_Hdelta = 0.001
+Nudge_Vwin_Ldelta = 0.001
+Nudge_Vwin_Invert = .true.
+/
diff --git a/myPythonTools/user_nl_cice b/myPythonTools/user_nl_cice
new file mode 100644
index 0000000000..ce8b72f238
--- /dev/null
+++ b/myPythonTools/user_nl_cice
@@ -0,0 +1,8 @@
+!----------------------------------------------------------------------------------
+! Users should add all user specific namelist changes below in the form of
+! namelist_var = new_namelist_value
+! Note - that it does not matter what namelist group the namelist_var belongs to
+!----------------------------------------------------------------------------------
+
+ dumpfreq_n = 0
+ histfreq_n = 0, 0, 0, 0, 0
diff --git a/src/chemistry/mozart/chemistry.F90 b/src/chemistry/mozart/chemistry.F90
index 920d96a5e3..28127c5ce5 100644
--- a/src/chemistry/mozart/chemistry.F90
+++ b/src/chemistry/mozart/chemistry.F90
@@ -707,6 +707,7 @@ subroutine chem_init(phys_state, pbuf2d)
use fire_emissions, only : fire_emissions_init
use short_lived_species, only : short_lived_species_initic
use ocean_emis, only : ocean_emis_init
+ use scamMod, only : single_column
type(physics_buffer_desc), pointer :: pbuf2d(:,:)
type(physics_state), intent(in):: phys_state(begchunk:endchunk)
@@ -809,7 +810,11 @@ subroutine chem_init(phys_state, pbuf2d)
if ( 1.e-2_r8 >= ptop_ref .and. ptop_ref > 1.e-5_r8 ) then ! around waccm top, below top of waccmx
cnst_fixed_ubc(1) = .true.
else if ( 1.e1_r8 > ptop_ref .and. ptop_ref > 1.e-2_r8 ) then ! well above top of cam and below top of waccm
- call endrun('chem_init: do not know how to set water vapor upper boundary when model top is near mesopause')
+ if(.not.(single_column)) then
+ call endrun('chem_init: do not know how to set water vapor upper boundary when model top is near mesopause')
+ else
+ cnst_fixed_ubc(1) = .true.
+ endif
endif
if ( masterproc ) write(iulog,*) 'chem_init: addfld done'
diff --git a/src/chemistry/mozart/upper_bc.F90 b/src/chemistry/mozart/upper_bc.F90
index 71a4a65b0c..8b076623e7 100644
--- a/src/chemistry/mozart/upper_bc.F90
+++ b/src/chemistry/mozart/upper_bc.F90
@@ -157,11 +157,16 @@ subroutine ubc_init()
use mo_snoe, only: snoe_inti
use mo_msis_ubc, only: msis_ubc_inti
use constituents,only: cnst_get_ind
+ use scamMod,only: single_column
!---------------------------Local workspace-----------------------------
logical :: zonal_avg
!-----------------------------------------------------------------------
- apply_upper_bc = ptop_ref<1._r8 ! Pa
+ if (single_column) then
+ apply_upper_bc = .FALSE.
+ else
+ apply_upper_bc = ptop_ref<1._r8 ! Pa
+ endif
if (.not.apply_upper_bc) return
diff --git a/src/control/history_scam.F90 b/src/control/history_scam.F90
index 2c81ce1a78..3288bfc7ca 100644
--- a/src/control/history_scam.F90
+++ b/src/control/history_scam.F90
@@ -45,7 +45,10 @@ subroutine scm_intht()
call addfld ('UDIFF', (/ 'lev' /), 'A', 'K','difference from observed u wind', gridname='gauss_grid')
call addfld ('VDIFF', (/ 'lev' /), 'A', 'K','difference from observed v wind', gridname='gauss_grid')
- call addfld ('TOBS', (/ 'lev' /), 'A', 'K','observed temp')
+ call addfld ('TOBS', (/ 'lev' /), 'A', 'K','observed temp', gridname='gauss_grid')
+ call addfld ('UOBS', (/ 'lev' /), 'A', 'm/s','observed zonal wind', gridname='gauss_grid')
+ call addfld ('VOBS', (/ 'lev' /), 'A', 'm/s','observed meridional wind', gridname='gauss_grid')
+
call addfld ('QDIFF', (/ 'lev' /), 'A', 'kg/kg','difference from observed water', gridname='gauss_grid')
call addfld ('QOBS', (/ 'lev' /), 'A', 'kg/kg','observed water', gridname='physgrid')
@@ -100,6 +103,59 @@ subroutine scm_intht()
call addfld ('NLTEN_PHYS', (/ 'lev' /), 'I','#/kg/s', 'NL vertical advective forcing', gridname='gauss_grid' )
call addfld ('NITEN_PHYS', (/ 'lev' /), 'I','#/kg/s', 'NI vertical advective forcing', gridname='gauss_grid' )
+!++jtb
+ call addfld ('U_IOP', (/ 'lev' /), 'I', 'm/s', 'Zonal Wind from IOP ', gridname='gauss_grid' )
+ call addfld ('V_IOP', (/ 'lev' /), 'I', 'm/s', 'Mer. Wind from IOP ', gridname='gauss_grid' )
+ call addfld ('OMEGA_IOP', (/ 'lev' /), 'I', 'Pa/s', 'Vertical velocity (from IOP) ', gridname='gauss_grid' )
+ call addfld ('OMEGA_ANA', (/ 'lev' /), 'I', 'Pa/s', 'Vertical velocity (analysis) ', gridname='gauss_grid' )
+ call addfld ('ETAD_ANA', (/ 'lev' /), 'I', 'Pa/s', 'Eta_dot (analysis) ', gridname='gauss_grid' )
+ call addfld ('ZETA_ANA', (/ 'lev' /), 'I', '1/s', 'Rel. Vorticity (analysis) ', gridname='gauss_grid' )
+ call addfld ('T_ANA', (/ 'lev' /), 'I', 'K', 'Temperature (analysis) ', gridname='gauss_grid' )
+ call addfld ('Q_ANA', (/ 'lev' /), 'I', 'g/g', 'Spec. humidity (analysis) ', gridname='gauss_grid' )
+ call addfld ('U_ANA', (/ 'lev' /), 'I', 'm/s', 'Zonal wind (analysis) ', gridname='gauss_grid' )
+ call addfld ('V_ANA', (/ 'lev' /), 'I', 'm/s', 'Mer. Wind (analysis) ', gridname='gauss_grid' )
+ call addfld ('TV_ANA', (/ 'lev' /), 'I', 'K', 'Temperature (analysis) ', gridname='gauss_grid' )
+ call addfld ('TTEN_TOTDYN_ANA', (/ 'lev' /), 'I', 'K/s', 'Temperature tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('UTEN_TOTDYN_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('VTEN_TOTDYN_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Meridional wind tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('QTEN_TOTDYN_ANA', (/ 'lev' /), 'I', 'kg/kg/s', 'tracer tendency (analysis)', gridname='gauss_grid' )
+
+ call addfld ('UTEN_TOTDYN_ANAR', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('VTEN_TOTDYN_ANAR', (/ 'lev' /), 'I', 'm/s/s', 'Meridional wind tendency (analysis)', gridname='gauss_grid' )
+
+ call addfld ('UTEN_DYCORE_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('VTEN_DYCORE_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Meridional wind tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('TTEN_DYCORE_ANA', (/ 'lev' /), 'I', 'K/s', 'Temperature tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('OMEGA_DYCORE_ANA', (/ 'lev' /), 'I', 'Pa/s','Pressure tendency/velocity (analysis)', gridname='gauss_grid' )
+ call addfld ('OMEGA_RECALC_ANA', (/ 'lev' /), 'I', 'Pa/s','Pressure tendency/velocity (analysis)', gridname='gauss_grid' )
+
+ call addfld ('UTEN_PRG_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('UTEN_PHIG_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('UTEN_KEG_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('UTEN_VORT_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('UTEN_PFRC_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('UTEN_VADV_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('UTEN_HADV_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('UTEN_CORIOL', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' )
+
+
+ call addfld ('VTEN_VORT_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Meridional wind tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('VTEN_PFRC_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Meridional wind tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('VTEN_VADV_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Meridional wind tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('VTEN_HADV_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Meridional wind tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('VTEN_CORIOL', (/ 'lev' /), 'I', 'm/s/s', 'Meridional wind tendency (analysis)', gridname='gauss_grid' )
+
+ call addfld ('TTEN_VADV_ANA', (/ 'lev' /), 'I', 'K/s', 'Temperature tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('TTEN_HADV_ANA', (/ 'lev' /), 'I', 'K/s', 'Temperature tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('TTEN_COMP_ANA', (/ 'lev' /), 'I', 'K/s', 'Temperature tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('TTEN_COMP_IOP', (/ 'lev' /), 'I', 'K/s', 'Temperature tendency (analysis)', gridname='gauss_grid' )
+
+ call addfld ('QTEN_VADV_ANA', (/ 'lev' /), 'I', '1/s', 'Temperature tendency (analysis)', gridname='gauss_grid' )
+ call addfld ('QTEN_HADV_ANA', (/ 'lev' /), 'I', '1/s', 'Temperature tendency (analysis)', gridname='gauss_grid' )
+
+!--jtb
+
+
end subroutine scm_intht
!#######################################################################
diff --git a/src/control/scamMod.F90 b/src/control/scamMod.F90
index b18169b340..b5506184cb 100644
--- a/src/control/scamMod.F90
+++ b/src/control/scamMod.F90
@@ -76,6 +76,24 @@ module scamMod
character*(max_path_len), public :: lsmsurffile
character*(max_path_len), public :: lsminifile
+!++jtb
+logical, public :: use_scm_ana_frc = .false.
+character*(max_path_len), public :: scm_ana_frc_file_template
+character*(max_path_len), public :: scm_ana_frc_path
+
+logical, public :: scm_ana_x_plevels = .true.
+logical, public :: scm_ana_direct_omega = .false.
+logical, public :: scm_ana_direct_ttend = .false.
+logical, public :: scm_ana_t_react = .false.
+logical, public :: scm_ana_q_react = .false.
+logical, public :: scm_ana_u_react = .false.
+logical, public :: scm_ana_v_react = .false.
+logical, public :: scm_ana_upwind = .false.
+!+++ARH
+logical, public :: scm_use_ana_iop = .false.
+!---ARH
+!--jtb
+
! note that scm_zadv_q is set to slt to be consistent with CAM BFB testing
@@ -250,7 +268,13 @@ subroutine scam_readnl(nlfile,single_column_in,scmlat_in,scmlon_in)
scm_cambfb_mode,scm_crm_mode,scm_zadv_uv,scm_zadv_T,scm_zadv_q,&
scm_use_obs_T, scm_use_obs_uv, scm_use_obs_qv, &
scm_relax_linear, scm_relax_tau_top_sec, &
- scm_relax_tau_bot_sec, scm_force_latlon, scm_relax_fincl, scm_backfill_iop_w_init
+ scm_relax_tau_bot_sec, scm_force_latlon, scm_relax_fincl, scm_backfill_iop_w_init, &
+!+jtb
+ use_scm_ana_frc, scm_ana_frc_path, scm_ana_frc_file_template, &
+ scm_ana_x_plevels, scm_ana_direct_omega, &
+ scm_ana_t_react, scm_ana_q_react, scm_ana_u_react, scm_ana_v_react, &
+ scm_ana_upwind, scm_ana_direct_ttend, scm_use_ana_iop
+!--jtb
single_column=single_column_in
@@ -306,6 +330,9 @@ subroutine scam_readnl(nlfile,single_column_in,scmlat_in,scmlon_in)
use_camiop = .false.
endif
+write(*,*) "!!!!!!!!!! ScamMod !!!!!!!! "
+write(*,*) scm_force_latlon , scmlon, scmlat
+
! If we are not forcing the lat and lon from the namelist use the closest lat and lon that is found in the IOP file.
if (.not.scm_force_latlon) then
call shr_scam_GetCloseLatLon( ncid, scmlat, scmlon, ioplat, ioplon, latidx, lonidx )
@@ -316,7 +343,9 @@ subroutine scam_readnl(nlfile,single_column_in,scmlat_in,scmlon_in)
scmlat = ioplat
scmlon = ioplon
end if
-
+write(*,*) " after " , scmlon, scmlat
+
+
if (masterproc) then
write (iulog,*) 'Single Column Model Options: '
write (iulog,*) '============================='
diff --git a/src/dynamics/eul/dyn_comp.F90 b/src/dynamics/eul/dyn_comp.F90
index 442c9f3228..0ba8285207 100644
--- a/src/dynamics/eul/dyn_comp.F90
+++ b/src/dynamics/eul/dyn_comp.F90
@@ -842,8 +842,10 @@ subroutine process_inidat(fieldname, m_cnst, fh)
ret = pio_inq_varid(fh, cnst_name(m_cnst), varid)
ret = pio_get_att(fh, varid, 'units', trunits)
if (trunits(1:5) .ne. 'KG/KG' .and. trunits(1:5) .ne. 'kg/kg') then
- call endrun(sub//': ERROR: Units for tracer ' &
- //trim(cnst_name(m_cnst))//' must be in KG/KG')
+!+++ARH
+! call endrun(sub//': ERROR: Units for tracer ' &
+! //trim(cnst_name(m_cnst))//' must be in KG/KG')
+!---ARH
end if
else if (.not. analytic_ic_active()) then
diff --git a/src/dynamics/eul/get_ana_dynfrc_4scam.F90 b/src/dynamics/eul/get_ana_dynfrc_4scam.F90
new file mode 100644
index 0000000000..1fdd1c58e0
--- /dev/null
+++ b/src/dynamics/eul/get_ana_dynfrc_4scam.F90
@@ -0,0 +1,1650 @@
+module get_ana_dynfrc_4scam
+
+ use spmd_utils, only: masterproc
+ use cam_logfile, only: iulog
+ use shr_kind_mod, only: r8 => shr_kind_r8, i8 => shr_kind_i8, &
+ cs=>SHR_KIND_CS,cl=>SHR_KIND_CL
+ use shr_const_mod, only: rearth => shr_const_rearth , & ! =6.37122e6_R8 meters
+ pi => shr_const_pi , &
+ OOmega => shr_const_omega , &
+ rdair => shr_const_rdair , &
+ cpair => shr_const_cpdair
+
+ use scamMod, only: use_scm_ana_frc, &
+ scm_ana_frc_path, &
+ scm_ana_frc_file_template, &
+ scm_ana_x_plevels, &
+ scm_ana_direct_omega, &
+ scm_ana_t_react, &
+ scm_ana_q_react, &
+ scm_ana_u_react, &
+ scm_ana_v_react, &
+ scm_ana_upwind, &
+ scm_ana_direct_ttend
+
+
+
+ ! shr_const_mod is in ${CESMROOT}/cime/src/share/util/
+
+ implicit none
+ private
+ save
+
+ public get_ana_dynfrc_fv
+!
+! Private module data
+!
+
+ real(r8) , save , allocatable :: T_1(:,:,:) , U_1(:,:,:), V_1(:,:,:), Q_1(:,:,:),PS_1(:,:),PHIS_1(:,:)
+ real(r8) , save , allocatable :: T_2(:,:,:) , U_2(:,:,:), V_2(:,:,:), Q_2(:,:,:),PS_2(:,:),PHIS_2(:,:)
+ real(r8) , save , allocatable :: UTCORE_1(:,:,:) , UTCORE_2(:,:,:)
+ real(r8) , save , allocatable :: VTCORE_1(:,:,:) , VTCORE_2(:,:,:)
+ real(r8) , save , allocatable :: TTCORE_1(:,:,:) , TTCORE_2(:,:,:)
+ real(r8) , save , allocatable :: OGCORE_1(:,:,:) , OGCORE_2(:,:,:)
+ real(r8) , save , allocatable :: lat_ana(:),lon_ana(:),lev_ana(:)
+ integer , save :: nlev_ana, nlon_ana, nlat_ana
+
+ real(r8) , save , allocatable :: To_1(:,:,:) , Uo_1(:,:,:), Vo_1(:,:,:), Qo_1(:,:,:),PSo_1(:,:),PHISo_1(:,:)
+ real(r8) , save , allocatable :: To_2(:,:,:) , Uo_2(:,:,:), Vo_2(:,:,:), Qo_2(:,:,:),PSo_2(:,:),PHISo_2(:,:)
+ real(r8) , save , allocatable :: UTCOREo_1(:,:,:) , UTCOREo_2(:,:,:), UTCOREo_X(:,:,:)
+ real(r8) , save , allocatable :: VTCOREo_1(:,:,:) , VTCOREo_2(:,:,:), VTCOREo_X(:,:,:)
+ real(r8) , save , allocatable :: TTCOREo_1(:,:,:) , TTCOREo_2(:,:,:), TTCOREo_X(:,:,:)
+ real(r8) , save , allocatable :: OGCOREo_1(:,:,:) , OGCOREo_2(:,:,:), OGCOREo_X(:,:,:)
+
+
+
+ real(r8) , save , allocatable :: ETAD_X(:,:,:) , OMG_X(:,:,:)
+ real(r8) , save , allocatable :: ZETA_X(:)
+ real(r8) , save , allocatable :: KEh_X(:,:,:)
+ real(r8) , save , allocatable :: Tv_X(:,:,:)
+
+ real(r8) , save , allocatable :: pke_X(:,:,:),pko_X(:,:,:),phik_X(:,:,:),Thv_X(:,:,:)
+ real(r8) , save , allocatable :: ple_X(:,:,:) , plo_X(:,:,:), phi_X(:,:,:)
+
+ real(r8) , save , allocatable :: To_X(:,:,:) , Uo_X(:,:,:), Vo_X(:,:,:), Qo_X(:,:,:),PSo_X(:,:),PHISo_X(:,:)
+
+
+!=======================================================================
+contains
+!=======================================================================
+
+subroutine get_ana_dynfrc_fv ( scmlon, scmlat , &
+ omega_ana, etad_ana, zeta_ana, &
+ t_ana , tv_ana , &
+ q_ana , &
+ u_ana , &
+ v_ana , &
+ ps_ana , &
+ uten_hadv_ana , &
+ vten_hadv_ana , &
+ uten_pfrc_ana , &
+ vten_pfrc_ana , &
+ uten_vort_ana , &
+ vten_vort_ana , &
+ qten_hadv_ana , &
+ tten_hadv_ana , &
+ uten_vadv_ana , &
+ vten_vadv_ana , &
+ tten_vadv_ana , &
+ qten_vadv_ana , &
+ tten_comp_ana , &
+ uten_keg_ana , &
+ uten_phig_ana , &
+ uten_prg_ana , &
+ uten_dycore_ana , &
+ vten_dycore_ana , &
+ tten_dycore_ana , &
+ omega_dycore_ana , &
+ omega_recalc_ana , &
+ u_scm, v_scm, t_scm, q_scm, &
+ u_ana_diag, v_ana_diag, t_ana_diag, q_ana_diag )
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! US and VS are input (D-grid velocities)
+!--------------------------------------------
+! ub(i,j,L)= 0.5*(us(i-1,j,L) + us(i,j,L))
+! vb(i,j,L)= 0.5*(vs(i,j,L) + vs(i,j+1,L))
+!
+! uc(i,j,L)= 0.5*(ub(i,j,L) + ub(i,j-1,L))
+! vc(i,j,L)= 0.5*(vb(i,j-1,L) + vb(i+1,j-1,L))
+!---------------------------------------------
+! Grid arrangement in FV latlon h,i-files
+!---------------------------------------------
+! J=NY
+! ...
+!
+! ub,vb(I,J) us(I,J),vc(I,J+1)
+!
+!
+! vs(I,J),uc(I,J) ua,va,T,p(I,J) vs(I+1,J),uc(I+1,J)
+!
+!
+! vc(I,J)
+!
+!
+! ua,va,T,p(I,J-1)
+!
+! ...
+! J=1 ...
+!----------------------------------------------
+
+ use pmgrid, only : plev, plat, plevp, plon
+ use hycoef, only: hyai, hybi, ps0, hyam, hybm
+ use filenames, only: interpret_filename_spec
+ use time_manager, only: timemgr_time_ge,timemgr_time_inc,get_curr_date,get_step_size,is_first_step
+ use netcdf
+ use cam_abortutils, only: endrun
+ use ref_pres, only: pref_mid ! In Pascal
+
+ real(r8), intent(in) :: scmlon, scmlat
+ real(r8), intent(out) :: omega_ana( plev )
+ real(r8), intent(out) :: etad_ana(plev)
+ real(r8), intent(out) :: t_ana(plev) , tv_ana(plev)
+ real(r8), intent(out) :: zeta_ana(plev)
+ real(r8), intent(out) :: u_ana(plev)
+ real(r8), intent(out) :: v_ana(plev)
+ real(r8), intent(out) :: q_ana(plev)
+ real(r8), intent(out) :: ps_ana
+ real(r8), intent(out) :: uten_hadv_ana( plev )
+ real(r8), intent(out) :: vten_hadv_ana( plev )
+ real(r8), intent(out) :: uten_pfrc_ana( plev )
+ real(r8), intent(out) :: vten_pfrc_ana( plev )
+ real(r8), intent(out) :: qten_hadv_ana( plev )
+ real(r8), intent(out) :: tten_hadv_ana( plev )
+ real(r8), intent(out) :: qten_vadv_ana( plev )
+ real(r8), intent(out) :: tten_vadv_ana( plev )
+ real(r8), intent(out) :: uten_vadv_ana( plev )
+ real(r8), intent(out) :: vten_vadv_ana( plev )
+
+ real(r8), intent(out) :: tten_comp_ana( plev )
+
+ real(r8), intent(out) :: uten_keg_ana( plev )
+ real(r8), intent(out) :: uten_prg_ana( plev )
+ real(r8), intent(out) :: uten_phig_ana( plev )
+ real(r8), intent(out) :: uten_vort_ana( plev )
+ real(r8), intent(out) :: vten_vort_ana( plev )
+ real(r8), intent(out) :: uten_dycore_ana( plev )
+ real(r8), intent(out) :: vten_dycore_ana( plev )
+ real(r8), intent(out) :: tten_dycore_ana( plev )
+ real(r8), intent(out) :: omega_recalc_ana( plev )
+ real(r8), intent(out) :: omega_dycore_ana( plev )
+
+ real(r8), intent(in) :: u_scm(plev)
+ real(r8), intent(in) :: v_scm(plev)
+ real(r8), intent(in) :: t_scm(plev)
+ real(r8), intent(in) :: q_scm(plev)
+
+ real(r8), intent(out) :: u_ana_diag(plev)
+ real(r8), intent(out) :: v_ana_diag(plev)
+ real(r8), intent(out) :: t_ana_diag(plev)
+ real(r8), intent(out) :: q_ana_diag(plev)
+
+ integer, save :: iax, jax
+ integer, save :: Read_year2, Read_month2, Read_day2, Read_sec2, Read_YMD2
+ integer, save :: nlev_alc, nlon_alc, nlat_alc
+
+ !!logical , parameter :: l_vectinv = .FALSE.
+ !!real(r8) :: tv_ana(plev)
+ real(r8) :: rho_ana( plev ), plo_ana(plev)
+
+
+
+ real(r8) :: scmlonx
+
+ real(r8) :: ana_wgt1 , ana_wgt2 , dx0, dy, darea
+
+ integer :: nx, ny,i,j,k,L,LM, iav(1),jav(1),iac,jac
+
+ real(r8) , allocatable :: rlats(:),rlons(:)
+ real(r8) :: zeta(plev),absvo(plev)
+ ! Horz. gradient profiles (1=X, 2=Y)
+ real(r8) :: kehg_ana(plev,2),kehg_X(plev,2)
+ real(r8) :: phig_ana(plev,2),phig_X(plev,2)
+ real(r8) :: plog_ana(plev,2),plog_X(plev,2)
+ real(r8) :: teg_ana(plev,2), teg_X(plev,2)
+ real(r8) :: qg_ana(plev,2), qg_X(plev,2)
+ real(r8) :: ug_ana(plev,2), ug_X(plev,2)
+ real(r8) :: vg_ana(plev,2), vg_X(plev,2)
+ real(r8) :: lin_pfc_ana(plev,2) , lin_pfc_X(plev,2)
+
+ real(r8) :: omega_ana_x(plev)
+ real(r8) :: alpha_react(plev)
+
+ real(r8) :: lat_alc(3) , lon_alc(3)
+ real(r8) :: aalc(3,3,plev)
+
+
+ character(len=CL):: Ana_File_Template,Ana_file1,Ana_file2,Ana_Path
+
+
+ integer :: dyn_year,dyn_month,dyn_day,dyn_sec,year,month,day,sec
+ integer :: dyn_step,ymd1,ymd2,curr_sec,next_sec,curr_year,curr_month,curr_day,curr_ymd
+
+ integer :: analysis_step
+ integer :: ana_year1, ana_month1, ana_day1, ana_sec1
+ integer :: ana_year2, ana_month2, ana_day2, ana_sec2
+
+ logical :: l_Read_next_Ana, Alarm_Read_ana, Alarm_Bump_ana, initialize
+
+ write(iulog,*) " version 07 of get_ana_dynfrc_4scam ... "
+
+
+ Alarm_Read_Ana = .FALSE.
+ Alarm_Bump_Ana = .FALSE.
+
+ if ( scmlon < 0 ) then
+ scmlonx = scmlon + 360._r8
+ else
+ scmlonx = scmlon
+ end if
+
+ ! Default to 6 hour steps between ana
+ analysis_step = 6 * 3600
+
+
+ Ana_path = trim(scm_ana_frc_path)
+ Ana_File_Template = trim(Ana_path)//trim(scm_ana_frc_file_template)
+
+
+ call get_curr_date(Year,Month,Day,Sec)
+
+ curr_ymd = (Year*10000) + (Month*100) + Day
+ curr_sec = Sec
+
+ ana_sec1 = ( Sec / analysis_step ) * analysis_step
+ ana_day1 = Day
+ ana_month1 = Month
+ ana_year1 = Year
+
+ YMD1=(Ana_Year1*10000) + (Ana_Month1*100) + Ana_Day1
+
+
+ call timemgr_time_inc(YMD1,Ana_Sec1, &
+ YMD2,Ana_Sec2,Analysis_Step,0,0)
+
+ Ana_Year2 = YMD2 / 10000
+ Ana_Month2 = (YMD2 - Ana_Year2*10000)/100
+ Ana_Day2 = YMD2 - Ana_Year2*10000 - Ana_Month2*100
+
+ Ana_File1 = interpret_filename_spec(Ana_File_Template , &
+ yr_spec=Ana_Year1 , &
+ mon_spec=Ana_Month1, &
+ day_spec=Ana_Day1 , &
+ sec_spec=Ana_Sec1 )
+
+ Ana_File2 = interpret_filename_spec(Ana_File_Template , &
+ yr_spec=Ana_Year2 , &
+ mon_spec=Ana_Month2, &
+ day_spec=Ana_Day2 , &
+ sec_spec=Ana_Sec2 )
+
+
+ l_Read_next_Ana = .FALSE.
+ ! On first time step, read in 2 analysis files
+ if (is_first_step().and.masterproc) then
+ write(iulog,*) " It's now (First time step):" , curr_YMD, curr_sec
+ write(iulog,*) "Read Initial ana files "
+ write(iulog,*) Ana_file1
+ write(iulog,*) Ana_file2
+ Alarm_Read_Ana = .TRUE.
+ Alarm_Bump_Ana = .FALSE.
+ else
+ ! On subsequent steps test to see if "Curr" date is later or same as "Read".
+ ! If it is, then l_read_next_ana=.TRUE.
+ call timemgr_time_ge(Read_ymd2, Read_Sec2, curr_YMD, curr_Sec, l_Read_next_ana )
+ endif
+
+ if (l_Read_next_Ana) then
+ Alarm_Read_Ana = .TRUE.
+ Alarm_Bump_Ana = .TRUE.
+ endif
+
+ ! Aloocate space for analysis fields.
+ ! Read in both Initial Analysis files. Nothing to bump yet
+ if ( (Alarm_Read_Ana ) .AND. .NOT.(Alarm_Bump_Ana) ) then
+ initialize=.TRUE.
+ call read_netcdf_ana_fv_ini ( Ana_File1, nlon_ana, nlat_ana, nlev_ana ,iax, jax )
+
+ if ( plev /= nlev_ana) then
+ call endrun ("SCAM plev NE nlev_ana")
+ end if
+
+ ! Full global fields
+ allocate( lat_ana(nlat_ana) , lon_ana(nlon_ana), lev_ana(nlev_ana) )
+ allocate( U_1(nlon_ana, nlat_ana, nlev_ana), V_1(nlon_ana, nlat_ana, nlev_ana), T_1(nlon_ana, nlat_ana, nlev_ana), &
+ Q_1(nlon_ana, nlat_ana, nlev_ana), PS_1 (nlon_ana, nlat_ana ), PHIS_1 (nlon_ana, nlat_ana ) )
+ allocate( U_2(nlon_ana, nlat_ana, nlev_ana), V_2(nlon_ana, nlat_ana, nlev_ana), T_2(nlon_ana, nlat_ana, nlev_ana), &
+ Q_2(nlon_ana, nlat_ana, nlev_ana), PS_2 (nlon_ana, nlat_ana ), PHIS_2 (nlon_ana, nlat_ana ) )
+
+ allocate( UTCORE_1(nlon_ana, nlat_ana, nlev_ana), UTCORE_2(nlon_ana, nlat_ana, nlev_ana) )
+ allocate( VTCORE_1(nlon_ana, nlat_ana, nlev_ana), VTCORE_2(nlon_ana, nlat_ana, nlev_ana) )
+ allocate( TTCORE_1(nlon_ana, nlat_ana, nlev_ana), TTCORE_2(nlon_ana, nlat_ana, nlev_ana) )
+ allocate( OGCORE_1(nlon_ana, nlat_ana, nlev_ana), OGCORE_2(nlon_ana, nlat_ana, nlev_ana) )
+
+ ! SCM "patches"
+ nlon_alc=3
+ nlat_alc=3
+ nlev_alc=nlev_ana
+
+
+
+ ! Patches of full global fields
+ allocate( Uo_1(nlon_alc, nlat_alc, nlev_alc), Vo_1(nlon_alc, nlat_alc, nlev_alc), To_1(nlon_alc, nlat_alc, nlev_alc), &
+ Qo_1(nlon_alc, nlat_alc, nlev_alc), PSo_1 (nlon_alc, nlat_alc ), PHISo_1 (nlon_alc, nlat_alc ) )
+ allocate( Uo_2(nlon_alc, nlat_alc, nlev_alc), Vo_2(nlon_alc, nlat_alc, nlev_alc), To_2(nlon_alc, nlat_alc, nlev_alc), &
+ Qo_2(nlon_alc, nlat_alc, nlev_alc), PSo_2 (nlon_alc, nlat_alc ), PHISo_2 (nlon_alc, nlat_alc ) )
+
+ allocate( UTCOREo_1(nlon_alc, nlat_alc, nlev_alc), UTCOREo_2(nlon_alc, nlat_alc, nlev_alc), UTCOREo_X(nlon_alc, nlat_alc, nlev_alc) )
+ allocate( VTCOREo_1(nlon_alc, nlat_alc, nlev_alc), VTCOREo_2(nlon_alc, nlat_alc, nlev_alc), VTCOREo_X(nlon_alc, nlat_alc, nlev_alc) )
+ allocate( TTCOREo_1(nlon_alc, nlat_alc, nlev_alc), TTCOREo_2(nlon_alc, nlat_alc, nlev_alc), TTCOREo_X(nlon_alc, nlat_alc, nlev_alc) )
+ allocate( OGCOREo_1(nlon_alc, nlat_alc, nlev_alc), OGCOREo_2(nlon_alc, nlat_alc, nlev_alc), OGCOREo_X(nlon_alc, nlat_alc, nlev_alc) )
+
+ allocate( Uo_X(nlon_alc, nlat_alc, nlev_alc), Vo_X(nlon_alc, nlat_alc, nlev_alc), To_X(nlon_alc, nlat_alc, nlev_alc), &
+ Qo_X(nlon_alc, nlat_alc, nlev_alc), PSo_X (nlon_alc, nlat_alc ), PHISo_X (nlon_alc, nlat_alc ) )
+ allocate( ETAD_X(nlon_alc,nlat_alc,nlev_alc) )
+ allocate( OMG_X(nlon_alc,nlat_alc,nlev_alc) )
+ allocate( ple_X(nlon_alc, nlat_alc, nlev_alc+1), plo_X(nlon_alc, nlat_alc, nlev_alc), phi_X(nlon_alc, nlat_alc, nlev_alc+1) )
+ allocate( pke_X(nlon_alc, nlat_alc, nlev_alc+1), pko_X(nlon_alc, nlat_alc, nlev_alc), phik_X(nlon_alc, nlat_alc, nlev_alc+1) )
+ allocate( THv_X(nlon_alc, nlat_alc, nlev_alc ) )
+ allocate( zeta_X(nlev_alc) )
+ allocate( KEh_X(nlon_alc, nlat_alc, nlev_alc ) )
+ allocate( Tv_X(nlon_alc, nlat_alc, nlev_alc ) )
+
+ call read_netcdf_ana_fv ( Ana_File1, nlon_ana, nlat_ana, nlev_ana, &
+ U_1, V_1, &
+ T_1, Q_1, PS_1, PHIS_1, &
+ lon_ana, lat_ana, lev_ana &
+ , utcore_1, vtcore_1, ttcore_1, ogcore_1 &
+ )
+ write(*,*) " checks ... "
+ write(*,*) iax, jax
+
+ call read_netcdf_ana_fv ( Ana_File2, nlon_ana, nlat_ana, nlev_ana, &
+ U_2, V_2, &
+ T_2, Q_2, PS_2, PHIS_2, &
+ lon_ana, lat_ana, lev_ana &
+ , utcore_2, vtcore_2, ttcore_2, ogcore_2 &
+ )
+
+ ! Make patches
+ Uo_1 = U_1(iax-1:iax+1,jax-1:jax+1,:)
+ Vo_1 = V_1(iax-1:iax+1,jax-1:jax+1,:)
+ To_1 = T_1(iax-1:iax+1,jax-1:jax+1,:)
+ Qo_1 = Q_1(iax-1:iax+1,jax-1:jax+1,:)
+ PSo_1 = PS_1(iax-1:iax+1,jax-1:jax+1 )
+ PHISo_1 = PHIS_1(iax-1:iax+1,jax-1:jax+1 )
+ UTCOREo_1 = UTCORE_1(iax-1:iax+1,jax-1:jax+1,:)
+ VTCOREo_1 = VTCORE_1(iax-1:iax+1,jax-1:jax+1,:)
+ TTCOREo_1 = TTCORE_1(iax-1:iax+1,jax-1:jax+1,:)
+ OGCOREo_1 = OGCORE_1(iax-1:iax+1,jax-1:jax+1,:)
+
+ Uo_2 = U_2(iax-1:iax+1,jax-1:jax+1,:)
+ Vo_2 = V_2(iax-1:iax+1,jax-1:jax+1,:)
+ To_2 = T_2(iax-1:iax+1,jax-1:jax+1,:)
+ Qo_2 = Q_2(iax-1:iax+1,jax-1:jax+1,:)
+ PSo_2 = PS_2(iax-1:iax+1,jax-1:jax+1 )
+ PHISo_2 = PHIS_2(iax-1:iax+1,jax-1:jax+1 )
+ UTCOREo_2 = UTCORE_2(iax-1:iax+1,jax-1:jax+1,:)
+ VTCOREo_2 = VTCORE_2(iax-1:iax+1,jax-1:jax+1,:)
+ TTCOREo_2 = TTCORE_2(iax-1:iax+1,jax-1:jax+1,:)
+ OGCOREo_2 = OGCORE_2(iax-1:iax+1,jax-1:jax+1,:)
+
+
+ ! Mark Ana date as read
+ Read_year2 = Ana_year2
+ Read_month2 = Ana_month2
+ Read_day2 = Ana_day2
+ Read_sec2 = Ana_sec2
+ Read_YMD2 =(Ana_Year2*10000) + (Ana_Month2*100) + Ana_Day2
+
+ end if
+
+ ! Bump second analysis to first postion, and read in next analysis
+ if ( (Alarm_Read_Ana ) .AND. (Alarm_Bump_Ana) ) then
+
+ Uo_1 = Uo_2
+ Vo_1 = Vo_2
+ To_1 = To_2
+ Qo_1 = Qo_2
+ PSo_1 = PSo_2
+ PHISo_1 = PHISo_2
+ UTCOREo_1 = UTCOREo_2
+ VTCOREo_1 = VTCOREo_2
+ TTCOREo_1 = TTCOREo_2
+
+ call read_netcdf_ana_fv ( Ana_File2, nlon_ana, nlat_ana, nlev_ana, &
+ U_2, V_2, &
+ T_2, Q_2, PS_2, PHIS_2, &
+ lon_ana, lat_ana, lev_ana &
+ , utcore_2, vtcore_2, ttcore_2, ogcore_2 &
+ )
+
+
+ ! Make patches
+ Uo_2 = U_2(iax-1:iax+1,jax-1:jax+1,:)
+ Vo_2 = V_2(iax-1:iax+1,jax-1:jax+1,:)
+ To_2 = T_2(iax-1:iax+1,jax-1:jax+1,:)
+ Qo_2 = Q_2(iax-1:iax+1,jax-1:jax+1,:)
+ PSo_2 = PS_2(iax-1:iax+1,jax-1:jax+1 )
+ PHISo_2 = PHIS_2(iax-1:iax+1,jax-1:jax+1 )
+ UTCOREo_2 = UTCORE_2(iax-1:iax+1,jax-1:jax+1,:)
+ VTCOREo_2 = VTCORE_2(iax-1:iax+1,jax-1:jax+1,:)
+ TTCOREo_2 = TTCORE_2(iax-1:iax+1,jax-1:jax+1,:)
+ OGCOREo_2 = OGCORE_2(iax-1:iax+1,jax-1:jax+1,:)
+
+
+ ! Mark Ana date as read
+ Read_year2 = Ana_year2
+ Read_month2 = Ana_month2
+ Read_day2 = Ana_day2
+ Read_sec2 = Ana_sec2
+ Read_YMD2=(Ana_Year2*10000) + (Ana_Month2*100) + Ana_Day2
+ end if
+
+ Alarm_Read_Ana = .FALSE.
+ Alarm_Bump_Ana = .FALSE.
+
+
+
+
+#if 0
+ call dynfrc_timewgts( &
+ (/ Ana_Year1, Ana_Month1, Ana_day1, Ana_sec1 /) , &
+ (/ Ana_Year2, Ana_Month2, Ana_day2, Ana_sec2 /) , &
+ ana_wgt1 , ana_wgt2 )
+#else
+ ana_wgt1 = 0._r8 ! 0=all weight on t+1
+ ana_wgt2 = 1._r8 - ana_wgt1
+#endif
+ if (masterproc) write(iulog,*) " Ana forcing time wgts ",ana_wgt1,ana_wgt2
+
+ iac=2
+ jac=2
+
+
+
+ Uo_X = ana_wgt1 * Uo_1 + ana_wgt2 * Uo_2
+ Vo_X = ana_wgt1 * Vo_1 + ana_wgt2 * Vo_2
+ To_X = ana_wgt1 * To_1 + ana_wgt2 * To_2
+ Qo_X = ana_wgt1 * Qo_1 + ana_wgt2 * Qo_2
+ PSo_X = ana_wgt1 * PSo_1 + ana_wgt2 * PSo_2
+ PHISo_X = ana_wgt1 * PHISo_1 + ana_wgt2 * PHISo_2
+
+ UTCOREo_X = ana_wgt1 * UTCOREo_1 + ana_wgt2 * UTCOREo_2
+ VTCOREo_X = ana_wgt1 * VTCOREo_1 + ana_wgt2 * VTCOREo_2
+ TTCOREo_X = ana_wgt1 * TTCOREo_1 + ana_wgt2 * TTCOREo_2
+ OGCOREo_X = ana_wgt1 * OGCOREo_1 + ana_wgt2 * OGCOREo_2
+
+ lon_alc = lon_ana(iax-1:iax+1)
+ lat_alc = lat_ana(jax-1:jax+1)
+
+ if(masterproc) write(iulog,*) " SCM lon lat: ",scmlonx,scmlat
+ if(masterproc) write(iulog,*) " Closest Ana lon lat: ",lon_ana( iax ) , lat_ana( jax )
+
+
+ ! Save off analysis fields for diagnostics and
+ ! other purposes
+ T_ana_diag(:) = To_X( iac, jac, :)
+ Q_ana_diag(:) = Qo_X( iac, jac, :)
+ U_ana_diag(:) = Uo_X( iac, jac, :)
+ V_ana_diag(:) = Vo_X( iac, jac, :)
+
+ !================================================
+ ! Patch in SCM profiles here if wanted.
+ ! This acts as "dynamical nudging", since
+ ! horizontal advective tendencies will become
+ ! stronger if SCM state drifts away from re-ana.
+ ! Note, this will only be effective w/ upwind
+ ! scheme, since 2nd order cntrd skips over central
+ ! point in stencil.
+ !----
+ ! For stability it turns out may be good to scale
+ ! with pressure so that high-velocity strato winds
+ ! don't lead to CFL violations. So, as a bad, dirty,
+ ! dirty short term solution, weight "reaction" by
+ ! pref_mid. Clearly, better soln would be to
+ ! sub-step this part of the dynamics as is done
+ ! for the other "dycores".
+ !=================================================
+ ! Calculate "reaction coefficient"
+ !---------------------------------
+ alpha_react(:)=1.0_r8 !1._r8
+
+ ! Adjust central profiles in stencils
+ !------------------------------------
+ if (scm_ana_t_react) then
+ To_X( iac, jac, :) = alpha_react(:) * T_scm(:) &
+ + ( 1._r8-alpha_react(:) ) * To_X( iac, jac, :)
+ if(masterproc) write(iulog,*) " REACTING to SCM T-state ..... "
+ else
+ if(masterproc) write(iulog,*) " No reaction to SCM T-state ..... "
+ endif
+ if (scm_ana_q_react) then
+ Qo_X( iac, jac, :) = alpha_react(:) * Q_scm(:) &
+ + ( 1._r8-alpha_react(:) ) * Qo_X( iac, jac, :)
+ if(masterproc) write(iulog,*) " REACTING to SCM Q-state ..... "
+ else
+ if(masterproc) write(iulog,*) " No reaction to SCM Q-state ..... "
+ endif
+ if (scm_ana_u_react) then
+ Uo_X( iac, jac, :) = alpha_react(:) * U_scm(:) &
+ + ( 1._r8-alpha_react(:) ) * Uo_X( iac, jac, :)
+ if(masterproc) write(iulog,*) " REACTING to SCM U-state ..... "
+ else
+ if(masterproc) write(iulog,*) " No reaction to SCM U-state ..... "
+ endif
+ if (scm_ana_v_react) then
+ Vo_X( iac, jac, :) = alpha_react(:) * V_scm(:) &
+ + ( 1._r8-alpha_react(:) ) * Vo_X( iac, jac, :)
+ if(masterproc) write(iulog,*) " REACTING to SCM V-state ..... "
+ else
+ if(masterproc) write(iulog,*) " No reaction to SCM V-state ..... "
+ endif
+
+
+
+ !=========================================
+
+ call virtual_t( nlon_alc,nlat_alc,nlev_alc, &
+ To_X , Qo_X , Tv_X )
+
+ call makepr_fv( nlon_alc,nlat_alc,nlev_alc, &
+ tv_X , pso_X , phiso_X , &
+ plo_X, ple_X, phi_X )
+ call etadot_fv ( nlon_alc , nlat_alc , nlev_alc , lon_alc , lat_alc , &
+ uo_X , &
+ vo_X , &
+ plo_X, ple_X , etad_X , omg_X )
+ call zeta_fv( nlon_alc,nlat_alc,nlev_alc, &
+ lon_alc ,lat_alc , &
+ uo_X , vo_X , zeta_X )
+
+ call makepk_fv( nlon_alc,nlat_alc,nlev_alc, &
+ To_X , Qo_X , &
+ pso_X , phiso_X , &
+ pko_X, pke_X, phik_X, thv_X )
+
+ KEh_X = 0.5 * ( Uo_X**2 + Vo_X**2 )
+
+
+ if (scm_ana_x_plevels) then
+ call patch_eta_x_plv ( nlon_alc , nlat_alc , nlev_alc, &
+ iac, jac, uo_X , plo_X )
+ call patch_eta_x_plv ( nlon_alc , nlat_alc , nlev_alc, &
+ iac, jac, vo_X , plo_X )
+ call patch_eta_x_plv ( nlon_alc , nlat_alc , nlev_alc, &
+ iac, jac, to_X , plo_X )
+ call patch_eta_x_plv ( nlon_alc , nlat_alc , nlev_alc, &
+ iac, jac, qo_X , plo_X )
+ call patch_eta_x_plv ( nlon_alc , nlat_alc , nlev_alc, &
+ iac, jac, tv_X , plo_X )
+ !Retain p-frc calculation on eta???
+ !call patch_eta_x_plv ( nlon_alc , nlat_alc , nlev_alc+1, &
+ ! iac, jac, phi_X , ple_X )
+ if(masterproc) write(iulog,*) " calcs on PRESSURE levels "
+ else
+ if(masterproc) write(iulog,*) " calcs on ETA levels "
+ end if
+
+
+ zeta_ana = zeta_X
+ omega_recalc_ana = omg_X( iac,jac,:)
+ etad_ana = etad_X( iac,jac,:)
+ plo_ana = plo_X( iac,jac,:)
+ t_ana = To_X( iac,jac,:)
+ tv_ana = Tv_X( iac,jac,:)
+ q_ana = Qo_X( iac,jac,:)
+ ps_ana = PSo_X( iac,jac )
+
+ u_ana = Uo_X( iac,jac,:)
+ v_ana = Vo_X( iac,jac,:)
+
+ rho_ana = plo_ana / ( Rdair * tv_ana )
+
+ uten_dycore_ana = UTCOREo_X( iac,jac,:)
+ vten_dycore_ana = VTCOREo_X( iac,jac,:)
+ tten_dycore_ana = TTCOREo_X( iac,jac,:)
+ omega_dycore_ana = OGCOREo_X( iac,jac,:)
+
+
+ ! Horz. gradient calcs
+
+ kehg_X = grad_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, KEh_X )
+
+ ! T_x, T_y should be straight T (not virtual)
+ !!teg_X = grad_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, To_X )
+ teg_X = grad_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, Tv_X ) !test 05-31-21
+
+ qg_X = grad_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, Qo_X )
+
+ ug_X = grad_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, Uo_X )
+
+ vg_X = grad_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, Vo_X )
+
+ aalc = 0.5*( PHI_X( :, :, 2:nlev_alc+1) + PHI_X(: , : ,1:nlev_alc) )
+ !!aalc = PHI_X( :, :, 2:nlev_alc+1)
+ !!aalc = PHI_X(: , : ,1:nlev_alc)
+ phig_X = grad_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, aalc )
+
+ !Retain p-frc calculation on eta???
+ !if (scm_ana_x_plevels) then ! No horz. p-gradient in p-coords
+ ! plog_X(:,1:2) = 0._r8
+ !else
+ plog_X = grad_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, plo_X(:,:,1:nlev_alc) )
+ !plog_X = grad_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, ple_X(:,:,1:nlev_alc) )
+ !end if
+
+
+
+#if 1
+ lin_pfc_X = lin_pfc_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, ple_X, phi_X )
+#else
+ lin_pfc_X = lin_pfc_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, pke_X, phik_X )
+#endif
+
+ kehg_ana = kehg_X
+ plog_ana = plog_X
+ phig_ana = phig_X
+ teg_ana = teg_X
+ qg_ana = qg_X
+ ug_ana = ug_X
+ vg_ana = vg_X
+ lin_pfc_ana = lin_pfc_X
+
+ !put together pieces for u*grad(u) form of U and V adv tendencies
+
+ if ( scm_ana_upwind .OR. scm_ana_u_react ) then
+ uten_hadv_ana = upwind_hadv(nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, u_ana, v_ana, Uo_X )
+ else
+ uten_hadv_ana = -u_ana * ug_ana(:,1) - v_ana * ug_ana(:,2)
+ end if
+ if ( scm_ana_upwind .OR. scm_ana_v_react ) then
+ vten_hadv_ana = upwind_hadv(nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, u_ana, v_ana, Vo_X )
+ else
+ vten_hadv_ana = -u_ana * vg_ana(:,1) - v_ana * vg_ana(:,2)
+ end if
+
+ ! Coriolis terms
+ !======================================
+ absvo = 2._r8 * OOmega * sin( lat_ana(jax) * PI/180._r8 )
+ !Allow Coriolis to react to SCM winds
+ uten_vort_ana = absvo * v_ana
+ vten_vort_ana = -absvo * u_ana
+ ! Force Coriolis to ALWAYS be calc w/ analysis winds
+ !!uten_vort_ana = absvo * v_ana_diag
+ !!vten_vort_ana = -absvo * u_ana_diag
+ ! ----- Diags for VI form (0-out)
+ uten_keg_ana = 0._r8 ! fill with 0
+
+ !!if (scm_ana_x_plevels) then ! No horz. p-gradient in p-coords
+ if (.FALSE.) then ! No horz. p-gradient in p-coords
+ uten_pfrc_ana = - phig_ana(:,1)
+ vten_pfrc_ana = - phig_ana(:,2)
+ else
+#if 1
+ !put together pieces for Pressure and Phi gradient tencency terms
+ uten_pfrc_ana = -(1._r8/rho_ana) * plog_ana(:,1) - phig_ana(:,1)
+ vten_pfrc_ana = -(1._r8/rho_ana) * plog_ana(:,2) - phig_ana(:,2)
+#else
+ !Lin(1997) QJRMS pfrc tendency terms
+ uten_pfrc_ana = lin_pfc_ana(:,1)
+ vten_pfrc_ana = lin_pfc_ana(:,2)
+#endif
+ end if
+
+
+ if ( scm_ana_upwind .OR. scm_ana_t_react ) then
+ tten_hadv_ana = upwind_hadv(nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, u_ana, v_ana, Tv_X )
+ else
+ tten_hadv_ana = -u_ana * teg_ana(:,1) - v_ana * teg_ana(:,2) ! should be straight T (not virtual)
+ end if
+ if ( scm_ana_upwind .OR. scm_ana_q_react ) then
+ qten_hadv_ana = upwind_hadv(nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, u_ana, v_ana, Qo_X )
+ else
+ qten_hadv_ana = -u_ana * qg_ana(:,1) - v_ana * qg_ana(:,2)
+ end if
+
+ if (.not.(scm_ana_direct_omega)) then
+ omega_ana = omega_recalc_ana ! use reconstructed omega
+ if(masterproc) write(iulog,*) " Omega recalc from ana U,V etc."
+ else
+ omega_ana = omega_dycore_ana ! use direct omega from dycore/ana
+ if(masterproc) write(iulog,*) " Omega direct from ana"
+ end if
+
+
+ if (.not.(scm_ana_x_plevels)) then
+ !Tendencies due to vertical advection (etadot * D_eta ... )
+ uten_vadv_ana = vadv_fv( nlev_alc, etad_ana, u_ana )
+ vten_vadv_ana = vadv_fv( nlev_alc, etad_ana, v_ana )
+ tten_vadv_ana = vadv_fv( nlev_alc, etad_ana, tv_ana ) ! should be straight T (not virtual)
+ qten_vadv_ana = vadv_fv( nlev_alc, etad_ana, q_ana )
+ else
+ !Tendencies due to vertical advection (Omega * D_p ... )
+ uten_vadv_ana = vadv_fv_press( nlev_alc, omega_ana, plo_ana, u_ana )
+ vten_vadv_ana = vadv_fv_press( nlev_alc, omega_ana, plo_ana, v_ana )
+ tten_vadv_ana = vadv_fv_press( nlev_alc, omega_ana, plo_ana, t_ana ) ! should be straight T (not virtual)
+ qten_vadv_ana = vadv_fv_press( nlev_alc, omega_ana, plo_ana, q_ana )
+ end if
+
+ tten_comp_ana = (1./cpair)*( omega_ana / rho_ana )
+
+ !DIags for pressure/geop grad forces
+ uten_phig_ana = - phig_ana(:,1)
+ uten_prg_ana = - (1._r8/rho_ana) * plog_ana(:,1)
+
+ end subroutine get_ana_dynfrc_fv
+
+!-----------------------------------------------------
+! Stuff ... useful ojala
+!-----------------------------------------------------
+ !-------------------------
+ function vadv_fv( nlev, etad, aa ) result( tend )
+ use hycoef, only: hyai, hybi, ps0, hyam, hybm
+ integer, intent(in) :: nlev
+ real(r8), intent(in) :: etad(nlev) , aa(nlev)
+ real(r8) :: tend(nlev)
+ real(r8) :: eta(nlev)
+ integer :: L
+
+ eta = hybm+hyam
+
+ do L=2,nlev-1
+ tend(L) = etad(L)* ( aa(L+1) - aa(L-1) ) / ( eta(L+1) - eta(L-1) )
+ end do
+ L=1
+ tend(L) = etad(L)* ( aa(L+1) - aa(L) ) / ( eta(L+1) - eta(L) )
+ L=nlev
+ tend(L) = etad(L)* ( aa(L) - aa(L-1) ) / ( eta(L) - eta(L-1) )
+
+ tend = -1.*tend ! for RHS consistency
+
+ end function vadv_fv
+!---------------------------
+ !-------------------------
+ function vadv_fv_press( nlev, omega, plo, aa ) result( tend )
+ integer, intent(in) :: nlev
+ real(r8), intent(in) :: omega(nlev) , aa(nlev),plo(nlev)
+ real(r8) :: tend(nlev)
+ integer :: L
+
+ do L=2,nlev-1
+ tend(L) = omega(L)* ( aa(L+1) - aa(L-1) ) / ( plo(L+1) - plo(L-1) )
+ end do
+ L=1
+ tend(L) = omega(L)* ( aa(L+1) - aa(L) ) / ( plo(L+1) - plo(L) )
+ L=nlev
+ tend(L) = omega(L)* ( aa(L) - aa(L-1) ) / ( plo(L) - plo(L-1) )
+
+ tend = -1.*tend ! for RHS consistency
+
+ end function vadv_fv_press
+!---------------------------
+ function lin_pfc_fv( nlon,nlat,nlev,iax,jax,lons,lats, pre, phi ) result( pfc )
+ !use shr_kind_mod, only: r8 => shr_kind_r8
+ !use shr_const_mod, only: rearth => shr_const_rearth , & ! =6.37122e6_R8 meters
+ ! pi => shr_const_pi , &
+ ! omega => shr_const_omega
+
+ integer, intent(in) :: nlon,nlat,nlev,iax,jax
+ real(r8), intent(in) :: pre(nlon,nlat,nlev+1),phi(nlon,nlat,nlev+1)
+ real(r8), intent(in) :: lats(nlat),lons(nlon)
+ real(r8) :: pfc(nlev,2)
+ real(r8) :: pfxW(nlev) , pfxE(nlev)
+ real(r8) :: pfyS(nlev) , pfyN(nlev)
+ real(r8) :: rlats(nlat),rlons(nlon),dx,dy,ds
+ real(r8) :: pr1,pr2,pr3,pr4, ph1,ph2,ph3,ph4
+ integer :: L , igg
+
+ ! Begin
+ rlons(:) = lons(:) * PI/180._r8
+ rlats(:) = lats(:) * PI/180._r8
+
+ dx=( rlons(2)-rlons(1) ) * Rearth
+ dy=( rlats(2)-rlats(1) ) * Rearth
+
+ ds = MAX( dx*cos(rlats(jax)) , .1 )
+ igg = iax
+ do L=1,nlev
+ pr1 = pre(igg-1,jax,L+1)
+ pr2 = pre(igg ,jax,L+1)
+ pr3 = pre(igg ,jax,L )
+ pr4 = pre(igg-1,jax,L )
+ ph1 = phi(igg-1,jax,L+1)
+ ph2 = phi(igg ,jax,L+1)
+ ph3 = phi(igg ,jax,L )
+ ph4 = phi(igg-1,jax,L )
+ pfxW(L) = ( (pr2-pr4)*(ph1-ph3) + (pr1-pr3)*(ph4-ph2) ) /( ds * ( (pr2-pr4) + (pr1-pr3) ) )
+ end do
+ igg = iax +1
+ do L=1,nlev
+ pr1 = pre(igg-1,jax,L+1)
+ pr2 = pre(igg ,jax,L+1)
+ pr3 = pre(igg ,jax,L )
+ pr4 = pre(igg-1,jax,L )
+ ph1 = phi(igg-1,jax,L+1)
+ ph2 = phi(igg ,jax,L+1)
+ ph3 = phi(igg ,jax,L )
+ ph4 = phi(igg-1,jax,L )
+ pfxE(L) = ( (pr2-pr4)*(ph1-ph3) + (pr1-pr3)*(ph4-ph2) ) /( ds * ( (pr2-pr4) + (pr1-pr3) ) )
+ end do
+ ds = dy
+ igg = jax
+ do L=1,nlev
+ pr1 = pre(iax,igg-1,L+1)
+ pr2 = pre(iax,igg ,L+1)
+ pr3 = pre(iax,igg ,L )
+ pr4 = pre(iax,igg-1,L )
+ ph1 = phi(iax,igg-1,L+1)
+ ph2 = phi(iax,igg ,L+1)
+ ph3 = phi(iax,igg ,L )
+ ph4 = phi(iax,igg-1,L )
+ pfyS(L) = ( (pr2-pr4)*(ph1-ph3) + (pr1-pr3)*(ph4-ph2) ) /( ds * ( (pr2-pr4) + (pr1-pr3) ) )
+ end do
+ igg = jax +1
+ do L=1,nlev
+ pr1 = pre(iax,igg-1,L+1)
+ pr2 = pre(iax,igg ,L+1)
+ pr3 = pre(iax,igg ,L )
+ pr4 = pre(iax,igg-1,L )
+ ph1 = phi(iax,igg-1,L+1)
+ ph2 = phi(iax,igg ,L+1)
+ ph3 = phi(iax,igg ,L )
+ ph4 = phi(iax,igg-1,L )
+ pfyN(L) = ( (pr2-pr4)*(ph1-ph3) + (pr1-pr3)*(ph4-ph2) ) /( ds * ( (pr2-pr4) + (pr1-pr3) ) )
+ end do
+
+
+ do L=1,nlev
+ pfc(L,1) = 0.5*( pfxW(L) + pfxE(L) )
+ pfc(L,2) = 0.5*( pfyS(L) + pfyN(L) )
+ end do
+
+
+
+ end function lin_pfc_fv
+ !-------------------------
+ function grad_fv( nlon,nlat,nlev,iax,jax,lons,lats, aa ) result( ga )
+ !use shr_kind_mod, only: r8 => shr_kind_r8
+ !use shr_const_mod, only: rearth => shr_const_rearth , & ! =6.37122e6_R8 meters
+ ! pi => shr_const_pi , &
+ ! omega => shr_const_omega
+
+ integer, intent(in) :: nlon,nlat,nlev,iax,jax
+ real(r8), intent(in) :: aa(nlon,nlat,nlev)
+ real(r8), intent(in) :: lats(nlat),lons(nlon)
+ real(r8) :: ga(nlev,2)
+ real(r8) :: rlats(nlat),rlons(nlon),dx,dy
+ integer :: L
+
+ ! Begin
+ rlons(:) = lons(:) * PI/180._r8
+ rlats(:) = lats(:) * PI/180._r8
+
+ dx=( rlons(2)-rlons(1) ) * Rearth
+ dy=( rlats(2)-rlats(1) ) * Rearth
+
+ do L=1,nlev
+ ga(L,1) = (aa(iax+1,jax,L) - aa(iax-1,jax,L))/( 2.*dx*cos(rlats(jax)) + 0.1 )
+ ga(L,2) = (aa(iax,jax+1,L) - aa(iax,jax-1,L))/( 2.*dy )
+ end do
+
+
+
+ end function grad_fv
+ !-------------------------
+ function upwind_hadv( nlon,nlat,nlev,iax,jax,lons,lats,u,v, aa ) result( hadv_tend )
+ !use shr_kind_mod, only: r8 => shr_kind_r8
+ !use shr_const_mod, only: rearth => shr_const_rearth , & ! =6.37122e6_R8 meters
+ ! pi => shr_const_pi , &
+ ! omega => shr_const_omega
+
+ integer, intent(in) :: nlon,nlat,nlev,iax,jax
+ real(r8), intent(in) :: aa(nlon,nlat,nlev)
+ real(r8), intent(in) :: lats(nlat),lons(nlon),u(nlev),v(nlev)
+ real(r8) :: hadv_tend(nlev)
+ real(r8) :: rlats(nlat),rlons(nlon),dx,dy,xten(nlev),yten(nlev)
+ integer :: L
+
+ ! Begin
+ rlons(:) = lons(:) * PI/180._r8
+ rlats(:) = lats(:) * PI/180._r8
+
+ dx=( rlons(2)-rlons(1) ) * Rearth
+ dy=( rlats(2)-rlats(1) ) * Rearth
+
+ do L=1,nlev
+ if ( u(L) >= 0._r8 ) then
+ xten(L) = u(L) * ( aa(iax,jax,L) - aa(iax-1,jax,L))/( dx*cos(rlats(jax)) + 0.1 )
+ else
+ xten(L) = u(L) * ( aa(iax+1,jax,L) - aa(iax,jax,L))/( dx*cos(rlats(jax)) + 0.1 )
+ end if
+ end do
+ do L=1,nlev
+ if ( v(L) >= 0._r8 ) then
+ yten(L) = v(L) * ( aa(iax,jax,L) - aa(iax,jax-1,L))/( dy )
+ else
+ yten(L) = v(L) * ( aa(iax,jax+1,L) - aa(iax,jax,L))/( dy )
+ end if
+ end do
+
+ hadv_tend(:) = -1._r8 * ( xten(:) + yten(:) )
+
+
+ end function upwind_hadv
+!=========================================
+ subroutine makepk_fv( nlon,nlat,nlev, t, q, ps, phis, pko, pke, phi, th )
+ use hycoef, only: hyai, hybi, ps0, hyam, hybm
+ !!use shr_const_mod, only: rdair => shr_const_rdair, cpair => shr_const_cpdair,
+ integer, intent(in) :: nlon,nlat,nlev
+ real(r8), intent(in) :: t(nlon,nlat,nlev),q(nlon,nlat,nlev),ps(nlon,nlat),phis(nlon,nlat)
+ real(r8), intent(out) :: pko(nlon,nlat,nlev),th(nlon,nlat,nlev),pke(nlon,nlat,nlev+1), phi(nlon,nlat,nlev+1)
+ real(r8) :: ple(nlon,nlat,nlev+1),plo(nlon,nlat,nlev),rv(nlon,nlat,nlev)
+ real(r8) :: kappa, p00
+ integer :: L
+
+ do L=1,nlev+1
+ ple(:,:,L) = hyai(L)*ps0 + hybi(L)*ps(:,:)
+ end do
+ do L=1,nlev
+ plo(:,:,L) = hyam(L)*ps0 + hybm(L)*ps(:,:)
+ end do
+
+ kappa=rdair/cpair
+
+ pko = plo**kappa
+ pke = ple**kappa
+
+ p00 = 100000._r8
+ th = ( ( p00 / plo)**kappa ) * t
+
+ rv = 1._r8/(1._r8 - q) - 1._r8
+ th = th*(1._r8 + 0.61_r8 * rv )
+
+ phi(:,:,nlev+1) = phis(:,:)
+ do L=nlev,1,-1
+ phi(:,:,L) = phi(:,:,L+1) - ( CpAir * Th(:,:,L) ) * ( pke(:,:,L) - pke(:,:,L+1) ) / (p00**kappa )
+ end do
+
+
+ end subroutine makepk_fv
+
+!=============================================================================
+ subroutine makepr_fv( nlon,nlat,nlev, t, ps, phis, plo, ple, phi )
+ use hycoef, only: hyai, hybi, ps0, hyam, hybm
+ use shr_const_mod, only: rdair => shr_const_rdair
+ integer, intent(in) :: nlon,nlat,nlev
+ real(r8), intent(in) :: t(nlon,nlat,nlev),ps(nlon,nlat),phis(nlon,nlat)
+ real(r8), intent(out) :: plo(nlon,nlat,nlev), ple(nlon,nlat,nlev+1), phi(nlon,nlat,nlev+1)
+ real(r8) :: lnple(nlon,nlat,nlev+1)
+ integer :: L
+
+ do L=1,nlev+1
+ ple(:,:,L) = hyai(L)*ps0 + hybi(L)*ps(:,:)
+ end do
+ do L=1,nlev
+ plo(:,:,L) = hyam(L)*ps0 + hybm(L)*ps(:,:)
+ end do
+
+ lnple = log( ple )
+ phi(:,:,nlev+1) = phis(:,:)
+ do L=nlev,1,-1
+ phi(:,:,L) = phi(:,:,L+1) - (RdAir * T(:,:,L) ) * ( lnple(:,:,L) - lnple(:,:,L+1) )
+ !phi(:,:,L) = phi(:,:,L+1) - (RdAir * T(:,:,L) / plo(:,:,L) ) * ( ple(:,:,L) - ple(:,:,L+1) )
+ end do
+
+ end subroutine makepr_fv
+
+!=============================================================================
+ subroutine virtual_t( nlon,nlat,nlev, t, q, tv )
+ use hycoef, only: hyai, hybi, ps0, hyam, hybm
+ use shr_const_mod, only: rdair => shr_const_rdair
+ integer, intent(in) :: nlon,nlat,nlev
+ real(r8), intent(in) :: t(nlon,nlat,nlev),q(nlon,nlat,nlev)
+ real(r8), intent(out) :: tv(nlon,nlat,nlev)
+ real(r8) :: rv(nlon,nlat,nlev)
+ integer :: L
+
+
+ rv = 1._r8/(1._r8 - q) - 1._r8
+ tv = t*(1._r8 + 0.61_r8 * rv )
+
+
+ end subroutine virtual_t
+
+ !-------------------------
+ subroutine zeta_fv( nlon,nlat,nlev,lons,lats, u,v, zeta )
+ !use shr_kind_mod, only: r8 => shr_kind_r8
+ !use shr_const_mod, only: rearth => shr_const_rearth , & ! =6.37122e6_R8 meters
+ ! pi => shr_const_pi , &
+ ! omega => shr_const_omega
+
+ integer, intent(in) :: nlon,nlat,nlev
+ real(r8), intent(in) :: u(nlon,nlat,nlev),v(nlon,nlat,nlev)
+ real(r8), intent(out) :: zeta(nlev)
+ !real(r8), intent(in) :: u(iax-1:iax+1,jax-1:jax+1,nlev)
+ !real(r8), intent(in) :: v(iax-1:iax+1,jax-1:jax+1,nlev)
+ real(r8), intent(in) :: lats(nlat),lons(nlon)
+ real(r8) :: rlats(nlat),rlons(nlon)
+ real(r8) :: dy,dx0,dx,darea,voo,voo2
+
+ integer :: iap,jap,iam,jam,i,j,L,iax,jax
+
+ iax=2
+ jax=2
+ write(*,*) " we're in subr. zeta_fv Lon Lat: "
+ write(*,*) lons(iax),lats(jax)
+
+ rlons(:) = lons(:) * PI/180._r8
+ rlats(:) = lats(:) * PI/180._r8
+
+ dx0 = rearth* ( rlons(2)-rlons(1) )
+ dy = rearth* ( rlats(2)-rlats(1) )
+
+ darea = dy*dx0*cos( rlats(jax) )
+
+ write(*,*) dx0,dy,cos( rlats(jax) )
+
+ do L =1,nlev
+ zeta(L) = &
+ ( V(iax+1,jax, L) - V(iax-1,jax,L) ) / ( 2*dx0*cos( rlats(jax) ) ) &
+ - ( U(iax,jax+1, L) - U(iax,jax-1,L) ) / ( 2*dy )
+ end do
+
+ write(*,*) " vorticity est. ",zeta(nlev)
+
+ end subroutine zeta_fv
+!================================================================
+ subroutine etadot_fv ( nlon, nlat, nlev, lons, lats, u, v, plo, ple, etadot , omega )
+ use shr_kind_mod, only: r8 => shr_kind_r8, i8 => shr_kind_i8
+ use shr_const_mod, only: rearth => shr_const_rearth , & ! =6.37122e6_R8 meters
+ pi => shr_const_pi
+ use hycoef, only: hyai, hybi, ps0, hyam, hybm
+
+ integer, intent(in) :: nlon,nlat,nlev
+ real(r8), intent(in) :: lons(nlon),lats(nlat)
+ real(r8), intent(in) :: u(nlon,nlat,nlev) , v(nlon,nlat,nlev) , plo( nlon,nlat,nlev) , ple( nlon,nlat,nlev+1)
+ real(r8), intent(out) :: etadot( nlon,nlat,nlev) ,omega(nlon,nlat,nlev)
+ !real(r8), intent(in) :: uc(:,:,:) , vc(:,:,:) , ple(:,:,:)
+
+ ! Local variables
+ real(r8),allocatable :: div(:,:,:)
+ real(r8),allocatable :: mass(:,:,:), fuc(:,:,:),fvc(:,:,:)
+ real(r8) :: rlats(nlat), rlons(nlon), rcos1, eta(nlev+1) , dx,dy! radians
+ real(r8), allocatable :: etadot_t1(:,:), etadot_t2(:,:,:)
+ integer :: i,j,L,im1,jm1,ip1,jp1
+ real :: uc_ijL , vc_ijL
+
+ allocate ( div(nlon,nlat,nlev) )
+ allocate ( mass(nlon,nlat,nlev), fuc(nlon,nlat,nlev),fvc(nlon,nlat+1,nlev) )
+ allocate ( etadot_t1(nlon,nlat), etadot_t2(nlon,nlat,nlev) )
+
+ div = 0._r8
+ fuc = 0._r8
+ fvc = 0._r8
+ mass = 0._r8
+ etadot = 0._r8
+ etadot_t1 = 0._r8
+ etadot_t2 = 0._r8
+
+ rlons(:) = lons(:) * PI/180._r8
+ rlats(:) = lats(:) * PI/180._r8
+
+ do L=1,nlev+1
+ eta(L) = hyai(L) + hybi(L) ! 1._r8*L/(nlev+1)
+ end do
+ do L=1,nlev
+ mass(:,:,L) = ( ple(:,:,L+1)-ple(:,:,L) )/( eta(L+1)-eta(L) )
+ end do
+
+ ! calculate mass fluxes at gridbox edges, using upwind algorithm
+ do L=1,nlev
+ do j=1,nlat
+ do i=2,nlon
+ im1=i-1
+ !if ( i == 1) im1=nlon
+ uc_ijL = 0.5*( u(im1,j,L) + u(i,j,L) )
+ if ( uc_ijL < 0. ) fuc(i,j,L)= uc_ijL * mass(i,j,L)
+ if ( uc_ijL >= 0. ) fuc(i,j,L)= uc_ijL * mass(im1,j,L)
+ end do
+ end do
+ end do
+ ! Note: cos(lat) term incorporated into fluxes
+ do L=1,nlev
+ do j=2,nlat
+ do i=1,nlon
+ jm1=j-1
+ vc_ijL = 0.5 * ( v(i,jm1,L)+v(i,j,L) )
+ if ( vc_ijL < 0. ) fvc(i,j,L)= vc_ijL * mass(i,j,L) *cos( rlats(j) )
+ if ( vc_ijL >= 0. ) fvc(i,j,L)= vc_ijL * mass(i,jm1,L) *cos( rlats(jm1) )
+ end do
+ end do
+ end do
+
+
+ ! now calculate HORZ divergence of (FUC,FVC). Note coslat term already
+ ! incorporated in FVC.
+ do L=1,nlev
+ do j=1,nlat-1
+ do i=1,nlon-1
+ ip1=i+1
+ jp1=j+1
+ rcos1 = 1. /( Rearth*cos( rlats(j) ) )
+ div(i,j,L) = rcos1 * ( FUC(ip1,j,L)-FUC(i,j,L) ) / (rlons(ip1)-rlons(i) ) &
+ + rcos1 * ( FVC(i,jp1,L)-FVC(i,j,L) ) / (rlats(jp1)-rlats(j) )
+ end do
+ end do
+ end do
+
+
+ etadot_t1(:,:)=0._r8
+ etadot_t2(:,:,:)=0._r8
+ do L=1,nlev
+ etadot_t1(:,:) = etadot_t1(:,:) + div(:,:,L)*(eta(L+1)-eta(L))
+ end do
+ do L=2,nlev
+ etadot_t2(:,:,L) = etadot_t2(:,:,L-1) + div(:,:,L)*(eta(L+1)-eta(L))
+ end do
+ do L=1,nlev
+ etadot(:,:,L) = ( hybm(L)*etadot_t1(:,:) - etadot_t2(:,:,L) ) / mass(:,:,L)
+ end do
+
+ dx=( rlons(2)-rlons(1) ) * Rearth
+ dy=( rlats(2)-rlats(1) ) * Rearth
+ omega = 0._r8
+
+#if 1
+ do L=1,nlev
+ do j=2,nlat-1
+ do i=2,nlon-1
+ omega(i,j,L) = u(i,j,L) * (plo(i+1,j,L)-plo(i-1,j,L))/( 2.*dx*cos(rlats(j)) + 0.1 ) &
+ + v(i,j,L) * (plo(i,j+1,L)-plo(i,j-1,L))/( 2.*dy ) &
+ - etadot_t2(i,j,L)
+ end do
+ end do
+ end do
+#else
+ do L=1,nlev
+ do j=2,nlat-1
+ do i=2,nlon-1
+ omega(i,j,L) = etadot(i,j,L)*mass(i,j,L)
+ end do
+ end do
+ end do
+#endif
+
+
+end subroutine etadot_fv
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!! Reading netcdf files
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!================================================================
+ subroutine read_netcdf_ana_fv_ini( anal_file , nlon, nlat, nlev,lonidx,latidx )
+ !
+ ! READ_NETCDF_ANAL_INI:
+ ! Open the given analyses data file. Query dimesnisons.
+ ! Close.
+ !===============================================================
+ use cam_abortutils, only : endrun
+ use netcdf
+ use shr_kind_mod, only: r8 => shr_kind_r8, i8 => shr_kind_i8
+ use hycoef, only: hyai, hybi, ps0, hyam, hybm
+ use shr_const_mod, only: rdair => shr_const_rdair
+ use scammod, only: scmlon,scmlat
+ use shr_scam_mod, only: shr_scam_getCloseLatLon ! Standardized system subroutines
+
+ !-------------
+ character(len=*),intent(in):: anal_file
+
+ integer, intent(out) :: nlon,nlat,nlev,latidx,lonidx
+
+ ! Local values
+ !-------------
+ integer :: ncid,varid,istat
+ integer :: ilat,ilon,ilev
+ integer :: i,j,L
+
+ real(r8) :: closelon,closelat
+
+ logical :: l_have_us , l_have_vs
+
+ l_have_us = .FALSE.
+ l_have_vs = .FALSE.
+
+ ! masterporc does all of the work here
+ !-----------------------------------------
+ if(masterproc) then
+
+ ! Open the given file
+ !-----------------------
+ istat=nf90_open(trim(anal_file),NF90_NOWRITE,ncid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*)'NF90_OPEN: failed for file ',trim(anal_file)
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+
+ ! Read in Dimensions
+ !--------------------
+ istat=nf90_inq_dimid(ncid,'lon',varid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+ istat=nf90_inquire_dimension(ncid,varid,len=nlon)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+
+ istat=nf90_inq_dimid(ncid,'lat',varid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+ istat=nf90_inquire_dimension(ncid,varid,len=nlat)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+
+ istat=nf90_inq_dimid(ncid,'lev',varid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+ istat=nf90_inquire_dimension(ncid,varid,len=nlev)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+
+ call shr_scam_getCloseLatLon(ncid ,scmlat,scmlon,closelat,closelon,latidx,lonidx)
+
+ ! Close the analyses file and exit
+ !-----------------------
+ istat=nf90_close(ncid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_EUL')
+ endif
+
+ endif ! (masterproc) then
+
+
+ end subroutine read_netcdf_ana_fv_ini
+
+!================================================================
+ subroutine read_netcdf_ana_fv( anal_file , nlon, nlat, nlev, &
+ u, v, &
+ t, q, ps, phis, &
+ lons, lats, levs &
+ , utcore, vtcore, ttcore, ogcore &
+ )
+ !
+ ! READ_NETCDF_ANAL :
+ ! Open the given analyses data file, read in
+ ! U,V,T,Q, and PS values as well as Lons, Lats.
+ !===============================================================
+ use cam_abortutils, only : endrun
+ use netcdf
+ use shr_kind_mod, only: r8 => shr_kind_r8, i8 => shr_kind_i8
+ use hycoef, only: hyai, hybi, ps0, hyam, hybm
+ use shr_const_mod, only: rdair => shr_const_rdair
+ ! Arguments
+ !-------------
+ character(len=*),intent(in):: anal_file
+
+ integer, intent(in ) :: nlon,nlat,nlev
+ real(r8), intent(out) :: U(nlon,nlat,nlev), V(nlon,nlat,nlev)
+ real(r8), intent(out) :: T(nlon,nlat,nlev), Q(nlon,nlat,nlev)
+ real(r8), intent(out) :: PS(nlon,nlat), PHIS(nlon,nlat)
+ !real(r8), intent(out) :: PHI(nlon,nlat,nlev+1),PLE(nlon,nlat,nlev+1),PLO(nlon,nlat,nlev)
+ real(r8), intent(out) :: Lats(nlat),Lons(nlon),Levs(nlev)
+
+ real(r8), intent(out) :: UTCORE(nlon,nlat,nlev), VTCORE(nlon,nlat,nlev), TTCORE(nlon,nlat,nlev)
+ real(r8), intent(out) :: OGCORE(nlon,nlat,nlev)
+
+ ! Local values
+ !-------------
+ integer :: ncid,varid,istat
+ integer :: ilat,ilon,ilev
+ integer :: i,j,L
+
+ logical :: l_have_us , l_have_vs
+
+ l_have_us = .FALSE.
+ l_have_vs = .FALSE.
+
+ ! masterporc does all of the work here
+ !-----------------------------------------
+ if(masterproc) then
+
+ ! Open the given file
+ !-----------------------
+ istat=nf90_open(trim(anal_file),NF90_NOWRITE,ncid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*)'NF90_OPEN: failed for file ',trim(anal_file)
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ end if
+ end if
+
+
+
+ if(masterproc) then
+
+ istat=nf90_inq_varid(ncid,'lon',varid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+ istat=nf90_get_var(ncid,varid,Lons)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+
+ istat=nf90_inq_varid(ncid,'lat',varid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+ istat=nf90_get_var(ncid,varid,Lats)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+
+ istat=nf90_inq_varid(ncid,'lev',varid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+ istat=nf90_get_var(ncid,varid,Levs)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+ endif ! (masterproc) then
+
+
+ if(masterproc) then
+ ! Read in, transpose lat/lev indices,
+ ! and scatter data arrays
+ !----------------------------------
+ ! First block reads U
+ !----------------------------------
+ istat=nf90_inq_varid(ncid,'U',varid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+ istat=nf90_get_var(ncid,varid, U )
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+ endif ! (masterproc) then
+
+ if(masterproc) then
+ istat=nf90_inq_varid(ncid,'V',varid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+ istat=nf90_get_var(ncid,varid, V )
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+ endif ! (masterproc) then
+
+
+
+
+!!!!!!!!!!!!!!
+ if(masterproc) then
+ istat=nf90_inq_varid(ncid,'T',varid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+ istat=nf90_get_var(ncid,varid, T )
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+ endif ! (masterproc) then
+
+ if(masterproc) then
+ istat=nf90_inq_varid(ncid,'Q',varid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+ istat=nf90_get_var(ncid,varid, Q )
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+ endif ! (masterproc) then
+
+ if(masterproc) then
+ istat=nf90_inq_varid(ncid,'PS',varid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+ istat=nf90_get_var(ncid,varid,PS )
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+ endif ! (masterproc) then
+
+ if(masterproc) then
+ istat=nf90_inq_varid(ncid,'PHIS',varid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_SE')
+ endif
+ istat=nf90_get_var(ncid,varid,PHIS )
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ endif
+ endif ! (masterproc) then
+
+ if(masterproc) then
+ istat=nf90_inq_varid(ncid,'UTEND_CORE',varid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) "No UTEND_CORE on file: "
+ write(iulog,*) trim(anal_file)
+ utcore(:,:,:)=-9999._r8
+ else
+ istat=nf90_get_var(ncid,varid,utcore )
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ end if
+ end if
+ end if ! (masterproc) then
+
+ if(masterproc) then
+ istat=nf90_inq_varid(ncid,'VTEND_CORE',varid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) "No VTEND_CORE on file: "
+ write(iulog,*) trim(anal_file)
+ vtcore(:,:,:)=-9999._r8
+ else
+ istat=nf90_get_var(ncid,varid,vtcore )
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ end if
+ end if
+ end if ! (masterproc) then
+
+ if(masterproc) then
+ istat=nf90_inq_varid(ncid,'DTCORE',varid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) "No TTEND_CORE on file: "
+ write(iulog,*) trim(anal_file)
+ ttcore(:,:,:)=-9999._r8
+ else
+ istat=nf90_get_var(ncid,varid,ttcore )
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ end if
+ end if
+ end if ! (masterproc) then
+
+ if(masterproc) then
+ istat=nf90_inq_varid(ncid,'OMEGA',varid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) "No OMEGA (core) on file: "
+ write(iulog,*) trim(anal_file)
+ ogcore(:,:,:)=-9999._r8
+ else
+ istat=nf90_get_var(ncid,varid,ogcore )
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_FV')
+ end if
+ end if
+ end if ! (masterproc) then
+
+
+ if(masterproc) then
+ ! Close the analysis file
+ !-----------------------
+ istat=nf90_close(ncid)
+ if(istat.ne.NF90_NOERR) then
+ write(iulog,*) nf90_strerror(istat)
+ call endrun ('UPDATE_ANALYSES_EUL')
+ endif
+ end if
+ !------------
+#if 0
+! Block winds at 45 m/s for increased stability. Kluge if Jets nt important
+ where(U > 45._r8)
+ U = 45._r8
+ end where
+ where(V > 45._r8)
+ V = 45._r8
+ end where
+ where(U < -45._r8)
+ U = -45._r8
+ end where
+ where(V < -45._r8)
+ V = -45._r8
+ end where
+#endif
+
+
+ write(*,*) "In read_netcdf_anal "
+ write(*,*) "Reading: ",anal_file
+ write(*,*) "Lons ..."
+ write(*,*) "Shape: ",shape(Lons)
+ write(*,*) "MinMax: ",minval(Lons),maxval(Lons)
+ write(*,*) "US and VS are presnt on file: ",l_have_us, l_have_vs
+
+
+ return
+ end subroutine read_netcdf_ana_fv
+!================================================================
+!================================================================
+ subroutine dynfrc_timewgts ( &
+ ana_prev_date, ana_next_date, &
+ wgt1 , wgt2 )
+
+
+ use shr_kind_mod, only: r8 => shr_kind_r8, i8 => shr_kind_i8
+ use ESMF
+ use time_manager, only:timemgr_time_ge,timemgr_time_inc,get_curr_date,get_step_size
+
+ integer, intent(in) :: ana_prev_date(4), ana_next_date(4)
+ real(r8) , intent(out) :: wgt1,wgt2
+
+ type(ESMF_Time) :: Date1,Date2,Date0
+ type(ESMF_TimeInterval) :: DateDiff2,DateDiff0,DateDiff, AnaDiff
+ integer :: DeltaT0, DeltaT2 , YMD, Year,Month,Day,Sec, Ana_interval, rc
+
+ call get_curr_date(Year,Month,Day,Sec)
+ YMD=(Year*10000) + (Month*100) + Day
+
+ call ESMF_TimeSet(Date0,YY=Ana_prev_date(1), MM=Ana_prev_date(2) , &
+ DD= Ana_prev_date(3) , S= Ana_prev_date(4) )
+ call ESMF_TimeSet(Date1,YY=Year,MM=Month,DD=Day,S=Sec)
+
+ call ESMF_TimeSet(Date2,YY=Ana_next_date(1), MM=Ana_next_date(2) , &
+ DD= Ana_next_date(3) , S= Ana_next_date(4) )
+ AnaDiff =Date2-Date0
+ call ESMF_TimeIntervalGet(AnaDiff,S=Ana_interval ,rc=rc)
+
+ DateDiff2 =Date2-Date1
+ call ESMF_TimeIntervalGet(DateDiff2,S=DeltaT2,rc=rc)
+ DateDiff0 =Date1-Date0
+ call ESMF_TimeIntervalGet(DateDiff0,S=DeltaT0,rc=rc)
+
+ wgt1 = 1._r8 - ( 1._r8 * DeltaT0 ) / Ana_interval
+ wgt2 = 1._r8 - ( 1._r8 * DeltaT2 ) / Ana_interval
+
+end subroutine dynfrc_timewgts
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine patch_eta_x_plv ( nx , ny, nL,ix, jx, aa, plo )
+ integer, intent(in) :: nx,ny,nl,ix,jx
+ real(r8), intent(in) :: plo(nx,ny,nL)
+ real(r8), intent(inout) :: aa(nx,ny,nL)
+
+ real(r8) :: plx(nL),plq(nL),aax(nL),aaq(nL),aat(nx,ny,nL)
+ real(r8) :: dp,dpk,dpk1,wtk,wtk1
+ integer :: i,j,L,k
+
+
+ plx(:) = plo(ix,jx,:) ! target pressures
+
+ do j=1,ny
+ do i=1,nx
+ plq(:) = plo(i,j,:)
+ aaq(:) = aa(i,j,:)
+ !if (plq(1) <= MINVAL(plx) ) aax(1) = aaq(1)
+ !if (plq(nl) > MAXVAL(plx) ) aax(nl) = aaq(nl)
+ do L=1,nl
+ do k=2,nl
+ if ( ( plx(L) <= plq(k) ).AND.(plx(L) > plq(k-1) ) ) then
+ dp = plq(k)-plq(k-1)
+ dpk1 = plx(L)-plq(k-1)
+ dpk = plq(k)-plx(L)
+ wtk1 = 1._r8 - dpk1 / dp
+ wtk = 1._r8 - dpk / dp
+ aax(L) = wtk * aaq(k) + wtk1 * aaq(k-1)
+ end if
+ end do
+ if ( plx(L) <= plq(1) ) aax(L)=aaq(1)
+ if ( plx(L) > plq(NL) ) aax(L)=aaq(NL)
+ end do
+
+ aat(i,j,:)=aax(:)
+ end do
+ end do
+
+ aa=aat
+
+!write(*,*) " mod "
+!write(411) nx,ny,nL
+!write(411) plo,aa,aat
+!PAUSE
+
+
+ end subroutine patch_eta_x_plv
+
+
+end module get_ana_dynfrc_4scam
diff --git a/src/dynamics/eul/scmforecast.F90 b/src/dynamics/eul/scmforecast.F90
index f9c0cbc6a8..b52a3bd92c 100644
--- a/src/dynamics/eul/scmforecast.F90
+++ b/src/dynamics/eul/scmforecast.F90
@@ -1,3 +1,4 @@
+#define SCAMNUDGERUN
module scmforecast
! --------------------------------------------------------------------------- !
! !
@@ -9,7 +10,11 @@ module scmforecast
use spmd_utils, only: masterproc
use cam_logfile, only: iulog
use cam_control_mod, only: adiabatic
-
+!++jtb
+#ifdef SCAMNUDGERUN
+ use get_ana_dynfrc_4scam, only: get_ana_dynfrc_fv
+#endif
+!--jtb
implicit none
private
save
@@ -59,10 +64,20 @@ subroutine forecast( lat , nlon , ztodt , &
scm_relax_tau_sec,scm_relax_tau_top_sec,scm_relax_top_p, &
scm_relaxation,scm_use_obs_qv,scm_use_obs_t,scm_use_obs_uv,scm_zadv_q,scm_zadv_t, &
scm_zadv_uv,tdiff,tobs,uobs,use_3dfrc,use_camiop,vertdivq, &
- vertdivt,vertdivu,vertdivv,vobs,wfld,qinitobs,scm_relax_fincl
+ vertdivt,vertdivu,vertdivv,vobs,wfld,qinitobs,scm_relax_fincl, &
+!++jtb
+ scmlon,scmlat, &
+ scm_ana_direct_ttend, &
+ scm_use_ana_iop
+!--jtb
use time_manager, only : get_curr_calday, get_nstep, get_step_size, is_first_step
use cam_abortutils, only : endrun
use string_utils, only: to_upper
+!++jtb
+ use shr_const_mod, only: rearth => shr_const_rearth , & ! =6.37122e6_R8 meters
+ pi => shr_const_pi , &
+ OOmega => shr_const_omega
+!--jtb
implicit none
@@ -71,6 +86,7 @@ subroutine forecast( lat , nlon , ztodt , &
! ---------------------- !
character(len=*), parameter :: subname = "forecast"
+ real(r8),parameter :: hugebad=9.99e12_r8
! --------------------------------------------------- !
! x = t, u, v, q !
@@ -83,16 +99,16 @@ subroutine forecast( lat , nlon , ztodt , &
integer, intent(in) :: nlon
real(r8), intent(in) :: ztodt ! Twice time step unless nstep = 0 [ s ]
- real(r8), intent(in) :: ps(plon) ! Surface pressure [ Pa ]
- real(r8), intent(in) :: psm1(plon) ! Surface pressure [ Pa ]
- real(r8), intent(in) :: psm2(plon) ! Surface pressure [ Pa ]
+ real(r8), intent(inout) :: ps(plon) ! Surface pressure [ Pa ]
+ real(r8), intent(inout) :: psm1(plon) ! Surface pressure [ Pa ]
+ real(r8), intent(inout) :: psm2(plon) ! Surface pressure [ Pa ]
real(r8), intent(in) :: t3m1(plev) ! Temperature [ K ]
- real(r8), intent(in) :: t3m2(plev) ! Temperature [ K ]
+ real(r8), intent(inout) :: t3m2(plev) ! Temperature [ K ]
real(r8), intent(in) :: u3m1(plev) ! Zonal wind [ m/s ]
- real(r8), intent(in) :: u3m2(plev) ! Zonal wind [ m/s ]
+ real(r8), intent(inout) :: u3m2(plev) ! Zonal wind [ m/s ]
real(r8), intent(in) :: v3m1(plev) ! Meridional wind [ m/s ]
- real(r8), intent(in) :: v3m2(plev) ! Meridional wind [ m/s ]
+ real(r8), intent(inout) :: v3m2(plev) ! Meridional wind [ m/s ]
real(r8), intent(inout) :: q3m1(plev,pcnst) ! Tracers [ kg/kg, #/kg ]
real(r8), intent(inout) :: q3m2(plev,pcnst) ! Tracers [ kg/kg, #/kg ]
@@ -156,6 +172,7 @@ subroutine forecast( lat , nlon , ztodt , &
real(r8) vten_zadv(plev) ! Vertical advective forcing of v [ m/s/s ]
real(r8) qten_zadv(plev,pcnst) ! Vertical advective forcing of tracers [ #/kg/s, kg/kg/s ]
+
! --------------------------- !
! For 'scm_relaxation' switch !
! --------------------------- !
@@ -169,6 +186,88 @@ subroutine forecast( lat , nlon , ztodt , &
real(r8) rslope ! [optional] slope for linear relaxation profile
real(r8) rycept ! [optional] y-intercept for linear relaxtion profile
+
+!++jtb
+! ------------------------------------ !
+! Quantities derived from Analyses !
+! ------------------------------------ !
+!======================================!
+ real(r8) dynfrcp(plev) ! Scaling factor for ana-derived tends
+ logical l_vectinv
+ real(r8) omega_ana(plev) ! Vertical pressure velocity [ Pa/s ]
+ real(r8) etad_ana(plev) ! "Eta dot" velocity [ Pa/s ]
+ real(r8) T_ana(plev), Q_ana(plev) , Tv_ana(plev) !
+ real(r8) u_ana(plev), v_ana(plev) !
+ real(r8) zeta_ana(plev) !
+ real(r8) ps_ana
+ real(r8) T_ana_diag(plev), Q_ana_diag(plev) !
+ real(r8) u_ana_diag(plev), v_ana_diag(plev) !
+ ! ----------------------------------- !
+ ! vertical advective tendencies !
+ ! ----------------------------------- !
+ real(r8) tten_vadv_ana(plev) ! Vertical advective forcing of t [ K/s ]
+ real(r8) uten_vadv_ana(plev) ! Vertical advective forcing of u [ m/s/s ]
+ real(r8) vten_vadv_ana(plev) ! Vertical advective forcing of v [ m/s/s ]
+ real(r8) qten_vadv_ana(plev) ! Vertical advective forcing of tracers [ #/kg/s, kg/kg/s ]
+ ! ------------------------------------- !
+ ! Horizontal advective/other tendencies !
+ ! ------------------------------------- !
+ real(r8) uten_hadv_ana(plev) ! of u [ m/s/s ]
+ real(r8) vten_hadv_ana(plev) ! of v [ m/s/s ]
+ real(r8) uten_pfrc_ana(plev) ! of u [ m/s/s ]
+ real(r8) vten_pfrc_ana(plev) ! of v [ m/s/s ]
+ real(r8) uten_vort_ana(plev) ! of u [ m/s/s ]
+ real(r8) vten_vort_ana(plev) ! of v [ m/s/s ]
+ real(r8) tten_hadv_ana(plev) ! of t [ K/s ]
+ real(r8) qten_hadv_ana(plev) ! of tracers [ #/kg/s, kg/kg/s ]
+
+ !---------------------------------!
+ ! Adiabatic compression tendency !
+ !---------------------------------!
+ real(r8) tten_comp_ana(plev) ! of t [ K/s ]
+
+
+ real(r8) uten_keg_ana(plev) ! of u [ m/s/s ]
+ real(r8) uten_prg_ana(plev) ! of u [ m/s/s ]
+ real(r8) uten_phig_ana(plev) ! of u [ m/s/s ]
+ ! ------------------------------------------ !
+ ! Direct dycore or ana tendencies or quants !
+ ! Not recalculated. !
+ ! (not usually available, !
+ ! set=-9999 if missing ) !
+ ! ------------------------------------------ !
+ real(r8) tten_dycore_ana(plev) ! Total direct Ana forcing of t [ K/s ]
+ real(r8) vten_dycore_ana(plev) ! Total direct Ana forcing of v [ m/s/s ]
+ real(r8) uten_dycore_ana(plev) ! Total direct Ana forcing of u [ m/s/s ]
+ real(r8) omega_dycore_ana(plev) ! Omega direct from Ana/dycore (not recalc) [ Pa/s ]
+ ! ----------------------------------- !
+ ! total recalc. "dycore" tendencies !
+ ! ----------------------------------- !
+ real(r8) omega_recalc_ana(plev) ! Omega from Ana/dycore (recalculated) [ Pa/s ]
+ real(r8) tten_totdyn_ana(plev) ! Total Ana forcing of t [ K/s ]
+ real(r8) uten_totdyn_ana(plev) ! Total Ana forcing of u [ m/s/s ]
+ real(r8) vten_totdyn_ana(plev) ! Total Ana forcing of v [ m/s/s ]
+ real(r8) qten_totdyn_ana(plev) ! Total Ana forcing of tracers [ #/kg/s, kg/kg/s ]
+ real(r8) fcoriol,uten_coriol(plev),vten_coriol(plev)
+ real(r8) ufcstm2(plev),vfcstm2(plev)
+ real(r8) ufcor_0(plev),vfcor_0(plev)
+ real(r8) uten_totdyn_anax(plev) ! Total Ana forcing of u [ m/s/s ]
+ real(r8) vten_totdyn_anax(plev) ! Total Ana forcing of v [ m/s/s ]
+ real(r8) tfw0, tfw1, tfw2, tftotw,ztodtn,AA
+ integer nsubdyn,nt,nstep_curr
+
+!+++ARH
+ !logical use_ana_iop
+!---ARH
+ logical l_use_reconst_ttend ! use reconstructed T-tendency based on analysis
+ logical l_use_direct_ttend ! use T-tendency direct from dycore
+
+
+ l_use_reconst_ttend = .NOT.( scm_ana_direct_ttend )
+ l_use_direct_ttend = .NOT.( l_use_reconst_ttend )
+
+!--jtb
+
!+++ BPM check what we have:
if (masterproc .and. is_first_step()) write(iulog,*) 'SCAM FORECAST REPORT: ' , &
'have_divq ', have_divq , &
@@ -253,8 +352,122 @@ subroutine forecast( lat , nlon , ztodt , &
! = .false. : Use User-generated SCAM IOP file !
! ------------------------------------------------------- !
-
- if( use_camiop ) then
+#ifdef SCAMNUDGERUN
+ !!! use_ana_iop needs to get into namelist!! !!!!
+!+++ARH
+ !use_ana_iop=.TRUE.
+ !!use_ana_iop=.FALSE.
+!---ARH
+ l_vectinv =.FALSE.
+
+!+++ARH
+ !if (use_ana_iop) then
+ if (scm_use_ana_iop) then
+!---ARH
+ call get_ana_dynfrc_fv ( scmlon, scmlat , &
+ omega_ana, etad_ana, zeta_ana, &
+ t_ana , tv_ana , &
+ q_ana , &
+ u_ana , &
+ v_ana , &
+ ps_ana , &
+ uten_hadv_ana , &
+ vten_hadv_ana , &
+ uten_pfrc_ana , &
+ vten_pfrc_ana , &
+ uten_vort_ana , &
+ vten_vort_ana , &
+ qten_hadv_ana , &
+ tten_hadv_ana , &
+ uten_vadv_ana , &
+ vten_vadv_ana , &
+ tten_vadv_ana , &
+ qten_vadv_ana , &
+ tten_comp_ana , &
+ uten_keg_ana , &
+ uten_phig_ana , &
+ uten_prg_ana , &
+ uten_dycore_ana , &
+ vten_dycore_ana , &
+ tten_dycore_ana , &
+ omega_dycore_ana , &
+ omega_recalc_ana , &
+ u3m2, v3m2, t3m2, q3m2(:,1), &
+ u_ana_diag, v_ana_diag, t_ana_diag, q_ana_diag )
+ else
+ ! set these to a "bad" value
+ omega_ana = HugeBad
+ etad_ana = HugeBad
+ zeta_ana = HugeBad
+ t_ana = HugeBad
+ tv_ana = HugeBad
+ q_ana = HugeBad
+ u_ana = HugeBad
+ v_ana = HugeBad
+ t_ana_diag = HugeBad
+ q_ana_diag = HugeBad
+ u_ana_diag = HugeBad
+ v_ana_diag = HugeBad
+ ps_ana = HugeBad
+ uten_hadv_ana = HugeBad
+ vten_hadv_ana = HugeBad
+ uten_pfrc_ana = HugeBad
+ vten_pfrc_ana = HugeBad
+ uten_vort_ana = HugeBad
+ vten_vort_ana = HugeBad
+ qten_hadv_ana = HugeBad
+ tten_hadv_ana = HugeBad
+ uten_vadv_ana = HugeBad
+ vten_vadv_ana = HugeBad
+ tten_vadv_ana = HugeBad
+ qten_vadv_ana = HugeBad
+ tten_comp_ana = HugeBad
+ uten_keg_ana = HugeBad
+ uten_phig_ana = HugeBad
+ uten_prg_ana = HugeBad
+ uten_dycore_ana = HugeBad
+ vten_dycore_ana = HugeBad
+ tten_dycore_ana = HugeBad
+ omega_dycore_ana = HugeBad
+ omega_recalc_ana = HugeBad
+ endif
+
+ ! -------------------------------------------------------------- !
+ ! Re-Calculate midpoint pressure levels if PS_ANA is reasonable !
+ ! -------------------------------------------------------------- !
+ if (ps_ana < 500000._r8 ) then
+ psm1=ps_ana
+ call plevs0( nlon, plon, plev, psm1, pintm1, pmidm1, pdelm1 )
+ end if
+ if(l_vectinv) then
+ uten_totdyn_ana = uten_hadv_ana + uten_pfrc_ana + uten_vadv_ana
+ vten_totdyn_ana = vten_hadv_ana + vten_pfrc_ana + vten_vadv_ana
+ uten_totdyn_anax = uten_hadv_ana + uten_pfrc_ana + uten_vadv_ana
+ vten_totdyn_anax = vten_hadv_ana + vten_pfrc_ana + vten_vadv_ana
+ else
+ uten_totdyn_ana = uten_hadv_ana + uten_vort_ana + uten_pfrc_ana + uten_vadv_ana
+ vten_totdyn_ana = vten_hadv_ana + vten_vort_ana + vten_pfrc_ana + vten_vadv_ana
+ uten_totdyn_anax = uten_hadv_ana + uten_vort_ana + uten_pfrc_ana + uten_vadv_ana
+ vten_totdyn_anax = vten_hadv_ana + vten_vort_ana + vten_pfrc_ana + vten_vadv_ana
+ endif
+
+ tten_totdyn_ana = tten_hadv_ana + tten_vadv_ana + tten_comp_ana
+ qten_totdyn_ana = qten_hadv_ana + qten_vadv_ana
+#else
+!+++ARH
+ !use_ana_iop=.FALSE.
+!---ARH
+#endif
+
+!++jtb
+ ! Need 3rd option 'use_ana_iop'
+ ! - suboption: use {u,v,t,q}ten_vadv_ana OR recalculate with etad_ana
+ ! - what about other species in q?
+ ! - we might want to calculate fu,fv using evolving (local) u's and v's
+ ! to allow geostrophic adjustment.
+!--jtb
+
+if( use_camiop ) then
do k = 1, plev
tfcst(k) = t3m2(k) + ztodt * tten_phys(k) + ztodt * divt3d(k)
ufcst(k) = u3m2(k) + ztodt * uten_phys(k) + ztodt * divu3d(k)
@@ -269,8 +482,11 @@ subroutine forecast( lat , nlon , ztodt , &
enddo
enddo
- else
-
+else ! when use_camiop =.FALSE.
+!+++ARH
+ !if( .NOT.(use_ana_iop) ) then
+ if( .NOT.(scm_use_ana_iop) ) then
+!---ARH
! ---------------------------------------------------------------------------- !
! Compute 'omega'( wfldint ) at the interface from the value at the mid-point. !
! SCAM-IOP file must provide omega at the mid-point not at the interface. !
@@ -403,19 +619,197 @@ subroutine forecast( lat , nlon , ztodt , &
call endrun( subname//':: divq not on the dataset. Unable to forecast Humidity. Stopping')
end if
- do k = 1, plev
- tfcst(k) = t3m2(k) + ztodt * ( tten_phys(k) + divt(k) + tten_zadv(k) )
- ufcst(k) = u3m2(k) + ztodt * ( uten_phys(k) + divu(k) + uten_zadv(k) )
- vfcst(k) = v3m2(k) + ztodt * ( vten_phys(k) + divv(k) + vten_zadv(k) )
- do m = 1, pcnst
- qfcst(1,k,m) = q3m2(k,m) + ztodt * ( qten_phys(k,m) + divq(k,m) + qten_zadv(k,m) )
+
+ nstep_curr = get_nstep()
+
+ do k = 1, plev
+ tfcst(k) = t3m2(k) + ztodt * ( tten_phys(k) + divt(k) + tten_zadv(k) )
+ ufcst(k) = u3m2(k) + ztodt * ( uten_phys(k) + divu(k) + uten_zadv(k) )
+ vfcst(k) = v3m2(k) + ztodt * ( vten_phys(k) + divv(k) + vten_zadv(k) )
+ do m = 1, pcnst
+ qfcst(1,k,m) = q3m2(k,m) + ztodt * ( qten_phys(k,m) + divq(k,m) + qten_zadv(k,m) )
+ enddo
enddo
- enddo
+ else
+ !-------------------------------------
+ ! This is the use_ana_iop=.TRUE. block
+ !-------------------------------------
+
+ nstep_curr = get_nstep()
+
+ if (is_first_step()) then
+ u3m2 = u_ana
+ v3m2 = v_ana
+ t3m2 = t_ana
+ q3m2(:,1) = q_ana
+ psm2 = ps_ana
+ endif
+
+
+ ! -----------------------------------------------------
+ ! Applied tendencies are in two
+ ! categories: 1) physics (includes nudging);
+ ! and 2) dynamics. Dynamics tendencies are
+ ! grouped and then scaled by dynfrcp. This is
+ ! to allow removal of unreliable analysis driven
+ ! dynamics tendencies above some pressure,
+ ! typically <~ 10Pa.
+ !------------------------------------------------------
+ dynfrcp(:) = 1._r8
+ where( pmidm1 < 10._r8) ! changed from 10. Test.
+ dynfrcp = 0._r8
+ end where
+ !------------------------------------------------------
+ fcoriol = 2._r8 * OOmega * sin( scmlat * PI/180._r8 )
+ uten_coriol = fcoriol * v3m2
+ vten_coriol = -fcoriol * u3m2
+ nsubdyn = 1
+ vfcst = v3m2
+ ufcst = u3m2
+ ztodtn = ztodt/nsubdyn
+ do nt= 1, nsubdyn
+ do k = 1, plev
+ ufcst(k) = ufcst(k) + ztodtn * ( uten_phys(k) &
+ + dynfrcp(k) * &
+ ( uten_hadv_ana(k) + uten_vadv_ana(k) &
+ + uten_vort_ana(k) &
+ !! + fcoriol * vfcstm2(k) &
+ + uten_pfrc_ana(k) ) )
+ vfcst(k) = vfcst(k) + ztodtn * ( vten_phys(k) &
+ + dynfrcp(k) * &
+ ( vten_hadv_ana(k) + vten_vadv_ana(k) &
+ + vten_vort_ana(k) &
+ !! - fcoriol * ufcstm2(k) &
+ + vten_pfrc_ana(k) ) )
+ end do
+ ufcstm2 = ufcst
+ vfcstm2 = vfcst
+ end do
+
+
+ ufcor_0 = ufcst
+ vfcor_0 = vfcst
+
+#if 0
+ ! Implicit formulation of Coriolis terms
+ nsubdyn = 1
+ ztodtn = ztodt/nsubdyn
+ AA = 1._r8/(1._r8 + (ztodtn*fcoriol)**2 )
+ do nt= 1, nsubdyn
+ do k = 1, plev
+ ufcst(k) = dynfrcp(k) * AA * ( ufcstm2(k) + ztodtn*fcoriol*vfcstm2(k) ) &
+ + (1._r8 - dynfrcp(k) )*ufcst(k)
+ vfcst(k) = dynfrcp(k) * AA * ( vfcstm2(k) - ztodtn*fcoriol*ufcstm2(k) ) &
+ + (1._r8 - dynfrcp(k) )*vfcst(k)
+ end do
+ ufcstm2 = ufcst
+ vfcstm2 = vfcst
+ end do
+
+ uten_vort_ana = (ufcst - ufcor_0 )/ztodt
+ vten_vort_ana = (vfcst - vfcor_0 )/ztodt
+#endif
+
+ uten_totdyn_ana = uten_hadv_ana + uten_vort_ana + uten_pfrc_ana + uten_vadv_ana
+ vten_totdyn_ana = vten_hadv_ana + vten_vort_ana + vten_pfrc_ana + vten_vadv_ana
+
+#if 1
+ !----------------------------
+ ! Calculate "usual" T-tendencies from complete IOP-file anyway
+ !----------------------------
+ ! ---------------------------------------------------------------------------- !
+ ! Compute 'omega'( wfldint ) at the interface from the value at the mid-point. !
+ ! SCAM-IOP file must provide omega at the mid-point not at the interface. !
+ ! ---------------------------------------------------------------------------- !
+ wfldint(1) = 0._r8
+ do k = 2, plev
+ weight = ( pintm1(k) - pmidm1(k-1) ) / ( pmidm1(k) - pmidm1(k-1) )
+ wfldint(k) = ( 1._r8 - weight ) * wfld(k-1) + weight * wfld(k)
+ enddo
+ wfldint(plevp) = 0._r8
+ ! ------------------------------------------------------------ !
+ ! Compute Eulerian compression heating due to vertical motion. !
+ ! ------------------------------------------------------------ !
+ do k = 1, plev
+ tten_comp_EUL(k) = wfld(k) * t3m1(k) * rair / ( cpair * pmidm1(k) )
+ enddo
+ ! ---------------------------------------------------------------------------- !
+ ! Compute Centered Eulerian vertical advective tendencies for all 't, u, v, q' !
+ ! ---------------------------------------------------------------------------- !
+ do k = 2, plev - 1
+ fac = 1._r8 / ( 2.0_r8 * pdelm1(k) )
+ tten_zadv_EULc(k) = -fac * ( wfldint(k+1) * ( t3m1(k+1) - t3m1(k) ) + wfldint(k) * ( t3m1(k) - t3m1(k-1) ) )
+ end do
+ k = 1
+ fac = 1._r8 / ( 2.0_r8 * pdelm1(k) )
+ tten_zadv_EULc(k) = -fac * ( wfldint(k+1) * ( t3m1(k+1) - t3m1(k) ) )
+ k = plev
+ fac = 1._r8 / ( 2.0_r8 * pdelm1(k) )
+ tten_zadv_EULc(k) = -fac * ( wfldint(k) * ( t3m1(k) - t3m1(k-1) ) )
+ !----------------------------------------
+ ! Replace ERA-derived T-tendencies with
+ ! IOP-file derived T-tendencies
+ !----------------------------------------
+ !!tten_vadv_ana(:) = tten_zadv_EULc(:)
+ !!tten_comp_ana(:) = tten_comp_EUL(:)
+ !!tten_hadv_ana(:) = divt(:)
+ !-------------------
+ ! For output
+ !--------------------
+ tten_zadv(:) = tten_zadv_EULc(:)
+ !----------------------------
+ ! End of Calculate "usual" T-tendencies from complete IOP-file anyway
+ !----------------------------
+#endif
+
+
+
+ if (l_use_reconst_ttend) then
+ do k=1,plev
+ tfcst(k) = t3m2(k) + ztodt * ( tten_phys(k) &
+ + dynfrcp(k) * &
+ ( tten_hadv_ana(k) + tten_vadv_ana(k) &
+ + tten_comp_ana(k) ) )
+ end do
+ end if
+
+ if (l_use_direct_ttend) then
+ do k=1,plev
+ tfcst(k) = t3m2(k) + ztodt * ( tten_phys(k) &
+ + dynfrcp(k) * &
+ ( tten_dycore_ana(k) ) )
+ end do
+ end if
+
+ do k=1,plev
+ do m = 1, 1
+ qfcst(1,k,m) = q3m2(k,m) + ztodt * ( qten_phys(k,m) &
+ + dynfrcp(k) * &
+ ( qten_hadv_ana(k) + qten_vadv_ana(k) ) )
+ enddo
+ enddo
+
+ ps = ps_ana
+
+ write(*,*) " Nstep " ,nstep_curr
+ if (mod( nstep_curr,10)==0) then
+ !ufcst = 0.5*(ufcst+u3m1)
+ !vfcst = 0.5*(vfcst+v3m1)
+ endif
+
+ ! Zero-out NON ana_iop diagnostics
+ ! ????
+
+ end if ! END use_ana_iop IF block
+
+ ! This code is executed regardless of use_ana_iop value
! ------------------ !
! Diagnostic Outputs !
! ------------------ !
-
+ call outfld( 'TOBS' , tobs, plon, dummy_dyndecomp )
+ call outfld( 'UOBS' , uobs, plon, dummy_dyndecomp )
+ call outfld( 'VOBS' , vobs, plon, dummy_dyndecomp )
call outfld( 'TTEN_XYADV' , divt, plon, dummy_dyndecomp )
call outfld( 'UTEN_XYADV' , divu, plon, dummy_dyndecomp )
call outfld( 'VTEN_XYADV' , divv, plon, dummy_dyndecomp )
@@ -438,29 +832,35 @@ subroutine forecast( lat , nlon , ztodt , &
call outfld( 'UTEN_ZADV' , uten_zadv, plon, dummy_dyndecomp )
call outfld( 'VTEN_ZADV' , vten_zadv, plon, dummy_dyndecomp )
call outfld( 'QVTEN_ZADV' , qten_zadv(:,1), plon, dummy_dyndecomp )
- call outfld( 'TTEN_ZADV' , vertdivt, plon, dummy_dyndecomp )
- call outfld( 'QVTEN_ZADV' , vertdivq(:,1), plon, dummy_dyndecomp )
+ !call outfld( 'TTEN_ZADV' , vertdivt, plon, dummy_dyndecomp )
+ !call outfld( 'QVTEN_ZADV' , vertdivq(:,1), plon, dummy_dyndecomp )
- call outfld( 'TTEN_PHYS' , tten_phys, plon, dummy )
- call outfld( 'UTEN_PHYS' , uten_phys, plon, dummy )
- call outfld( 'VTEN_PHYS' , vten_phys, plon, dummy )
- call outfld( 'QVTEN_PHYS' , qten_phys(:,1), plon, dummy )
+ call outfld( 'TTEN_COMP_IOP', tten_comp_eul, plon, dummy_dyndecomp )
- endif
+ call outfld( 'TTEN_PHYS' , tten_phys, plon, dummy_dyndecomp )
+ call outfld( 'UTEN_PHYS' , uten_phys, plon, dummy_dyndecomp )
+ call outfld( 'VTEN_PHYS' , vten_phys, plon, dummy_dyndecomp )
+ call outfld( 'QVTEN_PHYS' , qten_phys(:,1), plon, dummy_dyndecomp )
+
+ end if ! END of use_camiop IF BLOCK
+!!!!#if 0
+!+++ARH
+ !if( .NOT.(use_ana_iop) ) then
+ if( .NOT.(scm_use_ana_iop) ) then
+!---ARH
! ---------------------------------------------------------------- !
! Used the SCAM-IOP-specified state instead of forecasted state !
! at each time step if specified by the switch. !
! If SCAM-IOP has 't,u,v,q' profile at a single initial time step. !
- ! ---------------------------------------------------------------- !
-
+ ! ---------------------------------------------------------------- !
if( scm_use_obs_T .and. have_t ) then
do k = 1, plev
tfcst(k) = tobs(k)
enddo
endif
- if( scm_use_obs_uv .and. have_u .and. have_v ) then
+ if( scm_use_obs_uv .and. have_u .and. have_v ) then
do k = 1, plev
ufcst(k) = uobs(k)
vfcst(k) = vobs(k)
@@ -540,7 +940,9 @@ subroutine forecast( lat , nlon , ztodt , &
call outfld( 'TRELAX' , relax_T , plon, dummy )
call outfld( 'QRELAX' , relax_q(1:plev,1) , plon, dummy )
call outfld( 'TAURELAX' , rtau , plon, dummy )
-
+!!!#endif
+ end if ! END of 2nd use_ana_iop BLOCK (exec for use_ana_iop=.F.)
+
! --------------------------------------------------------- !
! Assign the final forecasted state to the output variables !
! --------------------------------------------------------- !
@@ -548,15 +950,79 @@ subroutine forecast( lat , nlon , ztodt , &
t3(1:plev) = tfcst(1:plev)
u3(1:plev) = ufcst(1:plev)
v3(1:plev) = vfcst(1:plev)
- q3(1:plev,1:pcnst) = qfcst(1,1:plev,1:pcnst)
-
+
+!+++ARH
+ !if (use_ana_iop) then
+ if (scm_use_ana_iop) then
+!---ARH
+ q3(1:plev,1:1) = qfcst(1,1:plev,1:1)
+ else
+ q3(1:plev,1:pcnst) = qfcst(1,1:plev,1:pcnst)
+ endif
+
tdiff(1:plev) = t3(1:plev) - tobs(1:plev)
qdiff(1:plev) = q3(1:plev,1) - qobs(1:plev)
+
call outfld( 'QDIFF' , qdiff, plon, dummy_dyndecomp )
call outfld( 'TDIFF' , tdiff, plon, dummy_dyndecomp )
+
+#ifdef SCAMNUDGERUN
+ call outfld( 'OMEGA_IOP' , wfld, plon, dummy_dyndecomp )
+ call outfld( 'OMEGA_ANA' , omega_ana, plon, dummy_dyndecomp )
+ call outfld( 'ETAD_ANA' , etad_ana, plon, dummy_dyndecomp )
+ call outfld( 'ZETA_ANA' , zeta_ana, plon, dummy_dyndecomp )
+ call outfld( 'T_ANA' , T_ana_diag, plon, dummy_dyndecomp )
+ call outfld( 'Q_ANA' , Q_ana_diag, plon, dummy_dyndecomp )
+ call outfld( 'TV_ANA' , Tv_ana, plon, dummy_dyndecomp )
+ call outfld( 'U_ANA' , U_ana_diag, plon, dummy_dyndecomp )
+ call outfld( 'V_ANA' , V_ana_diag, plon, dummy_dyndecomp )
+
+ call outfld( 'UTEN_CORIOL' , uten_coriol, plon, dummy_dyndecomp )
+ call outfld( 'VTEN_CORIOL' , vten_coriol, plon, dummy_dyndecomp )
+
+ call outfld( 'UTEN_TOTDYN_ANA' , uten_totdyn_ana, plon, dummy_dyndecomp )
+ call outfld( 'VTEN_TOTDYN_ANA' , vten_totdyn_ana, plon, dummy_dyndecomp )
+ call outfld( 'TTEN_TOTDYN_ANA' , tten_totdyn_ana, plon, dummy_dyndecomp )
+ call outfld( 'QTEN_TOTDYN_ANA' , qten_totdyn_ana, plon, dummy_dyndecomp )
+
+ call outfld( 'UTEN_TOTDYN_ANAR' , uten_totdyn_anax, plon, dummy_dyndecomp )
+ call outfld( 'VTEN_TOTDYN_ANAR' , vten_totdyn_anax, plon, dummy_dyndecomp )
+
+ call outfld( 'UTEN_DYCORE_ANA' , uten_dycore_ana, plon, dummy_dyndecomp )
+ call outfld( 'VTEN_DYCORE_ANA' , vten_dycore_ana, plon, dummy_dyndecomp )
+ call outfld( 'TTEN_DYCORE_ANA' , tten_dycore_ana, plon, dummy_dyndecomp )
+ call outfld( 'OMEGA_DYCORE_ANA', omega_dycore_ana, plon, dummy_dyndecomp )
+ call outfld( 'OMEGA_RECALC_ANA', omega_recalc_ana, plon, dummy_dyndecomp )
+
+ call outfld( 'UTEN_HADV_ANA' , uten_hadv_ana, plon, dummy_dyndecomp )
+ call outfld( 'UTEN_VADV_ANA' , uten_vadv_ana, plon, dummy_dyndecomp )
+ call outfld( 'UTEN_VORT_ANA' , uten_vort_ana, plon, dummy_dyndecomp )
+ call outfld( 'UTEN_KEG_ANA' , uten_keg_ana, plon, dummy_dyndecomp )
+ call outfld( 'UTEN_PFRC_ANA' , uten_pfrc_ana, plon, dummy_dyndecomp )
+ call outfld( 'UTEN_PRG_ANA' , uten_prg_ana, plon, dummy_dyndecomp )
+ call outfld( 'UTEN_PHIG_ANA' , uten_phig_ana, plon, dummy_dyndecomp )
+
+ call outfld( 'VTEN_HADV_ANA' , vten_hadv_ana, plon, dummy_dyndecomp )
+ call outfld( 'VTEN_VADV_ANA' , vten_vadv_ana, plon, dummy_dyndecomp )
+ call outfld( 'VTEN_VORT_ANA' , vten_vort_ana, plon, dummy_dyndecomp )
+ call outfld( 'VTEN_PFRC_ANA' , vten_pfrc_ana, plon, dummy_dyndecomp )
+
+ call outfld( 'TTEN_HADV_ANA' , tten_hadv_ana, plon, dummy_dyndecomp )
+ call outfld( 'TTEN_VADV_ANA' , tten_vadv_ana, plon, dummy_dyndecomp )
+ call outfld( 'TTEN_COMP_ANA' , tten_comp_ana, plon, dummy_dyndecomp )
+
+ call outfld( 'QTEN_HADV_ANA' , qten_hadv_ana, plon, dummy_dyndecomp )
+ call outfld( 'QTEN_VADV_ANA' , qten_vadv_ana, plon, dummy_dyndecomp )
+
+ if (have_u) call outfld( 'U_IOP' , uobs, plon, dummy_dyndecomp )
+ if (have_u) call outfld( 'V_IOP' , vobs, plon, dummy_dyndecomp )
+
+#endif
return
end subroutine forecast
- end module scmforecast
+
+
+end module scmforecast
diff --git a/src/physics/cam/check_energy.F90 b/src/physics/cam/check_energy.F90
index 3d0232b356..b9290e7319 100644
--- a/src/physics/cam/check_energy.F90
+++ b/src/physics/cam/check_energy.F90
@@ -66,6 +66,9 @@ module check_energy
integer :: teout_idx = 0 ! teout index in physics buffer
integer :: dtcore_idx = 0 ! dtcore index in physics buffer
+!+++ARH
+ integer :: dqcore_idx = 0
+!---ARH
integer :: ducore_idx = 0 ! ducore index in physics buffer
integer :: dvcore_idx = 0 ! dvcore index in physics buffer
@@ -139,6 +142,9 @@ subroutine check_energy_register()
call pbuf_add_field('TEOUT', 'global',dtype_r8 , (/pcols,dyn_time_lvls/), teout_idx)
call pbuf_add_field('DTCORE','global',dtype_r8, (/pcols,pver,dyn_time_lvls/),dtcore_idx)
+!+++ARH
+ call pbuf_add_field('DQCORE','global',dtype_r8, (/pcols,pver,dyn_time_lvls/),dqcore_idx)
+!---ARH
call pbuf_add_field('DUCORE','global',dtype_r8, (/pcols,pver,dyn_time_lvls/),ducore_idx)
call pbuf_add_field('DVCORE','global',dtype_r8, (/pcols,pver,dyn_time_lvls/),dvcore_idx)
if(is_subcol_on()) then
@@ -199,12 +205,18 @@ subroutine check_energy_init()
call addfld('TEFIX', horiz_only, 'A', 'J/m2', 'Total energy after fixer')
call addfld('EFIX', horiz_only, 'A', 'W/m2', 'Effective sensible heat flux due to energy fixer')
call addfld('DTCORE', (/ 'lev' /), 'A', 'K/s' , 'T tendency due to dynamical core')
+!+++ARH
+ call addfld('DQCORE', (/ 'lev' /), 'A', 'kg/kg/s' , 'Q tendency due to dynamical core')
+!---ARH
if ( history_budget ) then
call add_default ('DTCORE', history_budget_histfile_num, ' ')
end if
if ( history_waccm ) then
call add_default ('DTCORE', 1, ' ')
+!+++ARH
+ call add_default ('DQCORE', history_budget_histfile_num, ' ')
+!---ARH
end if
end subroutine check_energy_init
diff --git a/src/physics/cam/clubb_intr.F90 b/src/physics/cam/clubb_intr.F90
index fe263a70c7..6a9d9d5d1b 100644
--- a/src/physics/cam/clubb_intr.F90
+++ b/src/physics/cam/clubb_intr.F90
@@ -1,3 +1,4 @@
+#undef CLUBBTOP_OFF
module clubb_intr
!----------------------------------------------------------------------------------------------------- !
@@ -29,12 +30,32 @@ module clubb_intr
use zm_conv_intr, only: zmconv_microp
#ifdef CLUBB_SGS
use clubb_api_module, only: pdf_parameter, implicit_coefs_terms
- use clubb_api_module, only: clubb_config_flags_type
+ use clubb_api_module, only: clubb_config_flags_type, grid, stats, nu_vertical_res_dep
+ use clubb_api_module, only: nparams
use clubb_mf, only: do_clubb_mf, do_clubb_mf_diag
use cloud_fraction, only: dp1, dp2
#endif
implicit none
+#ifdef CLUBB_SGS
+ ! Variables that contains all the statistics
+
+ type (stats), target, save :: stats_zt(pcols), & ! stats_zt grid
+ stats_zm(pcols), & ! stats_zm grid
+ stats_rad_zt(pcols), & ! stats_rad_zt grid
+ stats_rad_zm(pcols), & ! stats_rad_zm grid
+ stats_sfc(pcols) ! stats_sfc
+
+!$omp threadprivate(stats_zt, stats_zm, stats_rad_zt, stats_rad_zm, stats_sfc)
+
+ type(grid), target :: dummy_gr
+ type(nu_vertical_res_dep), private, save :: dummy_nu_vert_res_dep ! Vertical resolution dependent nu values
+ real(r8), private, save :: dummy_lmin
+
+!$omp threadprivate(dummy_gr)
+
+#endif
+
private
save
@@ -46,9 +67,10 @@ module clubb_intr
#ifdef CLUBB_SGS
! This utilizes CLUBB specific variables in its interface
stats_init_clubb, &
- init_clubb_config_flags, &
-#endif
+ stats_zt, stats_zm, stats_sfc, &
+ stats_rad_zt, stats_rad_zm, &
stats_end_timestep_clubb, &
+#endif
clubb_readnl, &
clubb_init_cnst, &
clubb_implements_cnst
@@ -63,6 +85,7 @@ module clubb_intr
#ifdef CLUBB_SGS
type(clubb_config_flags_type), public :: clubb_config_flags
+ real(r8), dimension(nparams), public :: clubb_params ! Adjustable CLUBB parameters (C1, C2 ...)
#endif
! ------------ !
@@ -73,7 +96,9 @@ module clubb_intr
grid_type = 3, & ! The 2 option specifies stretched thermodynamic levels
hydromet_dim = 0 ! The hydromet array in SAM-CLUBB is currently 0 elements
- real(r8), parameter, dimension(0) :: &
+ ! Even though sclr_dim is set to 0, the dimension here is set to 1 to prevent compiler errors
+ ! See github ticket larson-group/cam#133 for details
+ real(r8), parameter, dimension(1) :: &
sclr_tol = 1.e-8_r8 ! Total water in kg/kg
character(len=6) :: saturation_equation
@@ -99,6 +124,7 @@ module clubb_intr
rtpthlp_const = 0.01_r8 ! Constant to add to rtpthlp when moments are advected
real(r8), parameter :: unset_r8 = huge(1.0_r8)
+ integer, parameter :: unset_i = huge(1)
! Commonly used temperature for the melting temp of ice crystals [K]
real(r8), parameter :: meltpt_temp = 268.15_r8
@@ -125,10 +151,25 @@ module clubb_intr
real(r8) :: clubb_c11 = unset_r8
real(r8) :: clubb_c11b = unset_r8
real(r8) :: clubb_c14 = unset_r8
+ real(r8) :: clubb_C_wp3_pr_turb = unset_r8
+ real(r8) :: clubb_c_K1 = unset_r8
+ real(r8) :: clubb_c_K2 = unset_r8
+ real(r8) :: clubb_nu2 = unset_r8
+ real(r8) :: clubb_c_K8 = unset_r8
real(r8) :: clubb_c_K9 = unset_r8
real(r8) :: clubb_nu9 = unset_r8
real(r8) :: clubb_c_K10 = unset_r8
real(r8) :: clubb_c_K10h = unset_r8
+ real(r8) :: clubb_C_invrs_tau_bkgnd = unset_r8
+ real(r8) :: clubb_C_invrs_tau_sfc = unset_r8
+ real(r8) :: clubb_C_invrs_tau_shear = unset_r8
+ real(r8) :: clubb_C_invrs_tau_N2 = unset_r8
+ real(r8) :: clubb_C_invrs_tau_N2_wp2 = unset_r8
+ real(r8) :: clubb_C_invrs_tau_N2_xp2 = unset_r8
+ real(r8) :: clubb_C_invrs_tau_N2_wpxp = unset_r8
+ real(r8) :: clubb_C_invrs_tau_N2_clear_wp3 = unset_r8
+ real(r8) :: clubb_C_uu_shr = unset_r8
+ real(r8) :: clubb_C_uu_buoy = unset_r8
real(r8) :: clubb_gamma_coef = unset_r8
real(r8) :: clubb_gamma_coefb = unset_r8
real(r8) :: clubb_beta = unset_r8
@@ -137,31 +178,129 @@ module clubb_intr
real(r8) :: clubb_mult_coef = unset_r8
real(r8) :: clubb_Skw_denom_coef = unset_r8
real(r8) :: clubb_skw_max_mag = unset_r8
- real(r8) :: clubb_up2_vp2_factor = unset_r8
+ real(r8) :: clubb_up2_sfc_coef = unset_r8
real(r8) :: clubb_C_wp2_splat = unset_r8
real(r8) :: clubb_wpxp_L_thresh = unset_r8
real(r8) :: clubb_detliq_rad = unset_r8
real(r8) :: clubb_detice_rad = unset_r8
real(r8) :: clubb_detphase_lowtemp = unset_r8
- logical :: clubb_l_brunt_vaisala_freq_moist = .false.
- logical :: clubb_l_call_pdf_closure_twice = .false.
- logical :: clubb_l_damp_wp3_Skw_squared = .false.
- logical :: clubb_l_min_wp2_from_corr_wx = .false.
- logical :: clubb_l_min_xp2_from_corr_wx = .false.
- logical :: clubb_l_predict_upwp_vpwp = .false.
- logical :: clubb_l_rcm_supersat_adj = .false.
- logical :: clubb_l_stability_correct_tau_zm = .false.
- logical :: clubb_l_trapezoidal_rule_zt = .false.
- logical :: clubb_l_trapezoidal_rule_zm = .false.
- logical :: clubb_l_upwind_xpyp_ta = .false.
- logical :: clubb_l_use_C7_Richardson = .false.
- logical :: clubb_l_use_C11_Richardson = .false.
- logical :: clubb_l_use_cloud_cover = .false.
- logical :: clubb_l_use_thvm_in_bv_freq = .false.
- logical :: clubb_l_vert_avg_closure = .false.
- logical :: clubb_l_diag_Lscale_from_tau = .false.
- logical :: clubb_l_damp_wp2_using_em = .false.
+
+ integer :: &
+ clubb_iiPDF_type, & ! Selected option for the two-component normal
+ ! (double Gaussian) PDF type to use for the w, rt,
+ ! and theta-l (or w, chi, and eta) portion of
+ ! CLUBB's multivariate, two-component PDF.
+ clubb_ipdf_call_placement = unset_i ! Selected option for the placement of the call to
+ ! CLUBB's PDF.
+
+ logical :: &
+ clubb_l_use_precip_frac, & ! Flag to use precipitation fraction in KK microphysics. The
+ ! precipitation fraction is automatically set to 1 when this
+ ! flag is turned off.
+ clubb_l_predict_upwp_vpwp, & ! Flag to predict and along with and
+ ! alongside the advancement of