Thursday, August 2, 2012

Jacob's MATLAB functions

I've used MATLAB for most of my plotting needs for over a year, so most of my "final" plotting scripts are written in MATLAB.  During that time, I wrote a number of custom functions that get shared across all my scripts.  Unfortunately, they aren't included in the scripts themselves, so if you want to use my scripts you'll need to manually add them to your MATLAB path so matlab can see them.  If you don't, you'll see errors like "Undefined function or variable 'strrep_many'" when you try to run my scripts.

These functions are stored on Sox.  They can be found in 3 directories:
  • /Users/oberman/Documents/MATLAB/ColorUtils
  • /Users/oberman/Documents/MATLAB/MiscUtils
  • /Users/oberman/Documents/MATLAB/TimeUtils

If you already know how to add folders to the MATLAB path (and/or already have a folder where you keep custom functions), just take the scripts out of those 3 directories and put them somewhere on your path.  You'll be good to go.


If not, you'll need to either modify your startup.m file or create one if you don't have it.  startup.m is a MATLAB script that runs as soon as the program boots up.  It can execute any commands you like; here, we're going to use it to add my folders to your path so MATLAB can see the functions within them. 

The default location that MATLAB looks for a startup.m file is in the folder ~/Documents/MATLAB/, where ~ is your home directory.  If you don't have a /Documents/MATLAB/ folder, try starting MATLAB and then exiting.  It should exist after you've started MATLAB once.  Create a file in that folder called "startup.m" and in it put the following lines:
addpath('/Users/oberman/Documents/MATLAB/ColorUtils');
addpath('/Users/oberman/Documents/MATLAB/TimeUtils');
addpath('/Users/oberman/Documents/MATLAB/MiscUtils');
Now, next time you start MATLAB, my functions should work just as if they were built in.  You are of course welcome to add other lines to your startup.m (it uses the MATLAB scripting language).  You can also copy my scripts to another location and add that folder instead - this might be helpful if you want to modify or expand on something I did.

Friday, February 17, 2012

To ensure wet deposition --or at least the chance for it-- would be realistically modeled, we investigated different physics settings of WRF. Here are a couple notes about our findings:


(1) Be careful how you plot WRF precipitation. It is accumulated precipitation, written every hour (as currently used, "history_interval" in namelist.input = 60 minutes). Some models will accumulate precip. in "buckets" with the bucket being emptied out every {interval of history write}. WRF is not like that, it has an accumulated precip. bucket that grows with time. WRF accumulated precipitation for any particular time = the accumulated precipitation from the start of the run to that particular time. 

(2) There is a higher agreement between NARR precipitation and WRF precipitation when using the Kain-Fritsch cumulus parameterization (cu_physics = 1 in namelist.input). This is not the same cu. param. that Claus and Caitlin used (Grell, cu_physics = 5).


(3) The sensitivity to land surface model is very small, so continuing to use the Pleim-Xu LSM (sf_sfclay_physics = 7 in namelist.input), as Claus and Caitlin had used, is fine. 


(4) Changing the cumulus parameterization has very little influence on nudged temperatures because they are nudged. Changing the cu. param. is just for precipitation's sake.

Thursday, February 16, 2012

Configure MCIP3.6 to output full layer pressure (PRESF)

To compare CMAQ model output against OMI satellite retrievals using Jacob's processing scheme, it may be necessary to interpolate vertical layers between CMAQ and satellite vertical columns. MCIP v3.6 by default outputs pressure levels at mid-layer (one value for each layer), however pressure levels at the full layer height may be needed for some interpolation schemes.

Here are instructions from Tanya Otte (EPA MCIP developer) to modify MCIP to output PRESF, the full layer pressure values. This was tested in my copy of MCIP 3.6 on nitrate (/home/bickford/MCIP3.6), and worked on the first try.

*************************************************************************************

From Tanya Otte (Otte.Tanya@epamail.epa.gov):

Here are the mods that you need to make to MCIP to add PRESF to the output. As I explained in the reply via m3list, you will only be outputting the levels above the surface (such that the number of full levels output is really one fewer than exist, but matches the number of mid-layers in the model). The pressure at the lowest full-level (the surface) is in METCRO2D in variable PRSFC.

Caveat 1: I have not tested this; it's mostly just looking at what I believe the changes are that need to be made.

Caveat 2: A new version of MCIP will be released with CMAQv5.0 (any day now). You may want to wait for it. **MCIP v4.0 was released Feb 2012**

Caveat 3: If you are using the tipping bucket for precip in WRF, let me know, and I'll give you the mods. The current release of MCIP does not handle the tipping bucket; it's off by default, and it's fairly new.

Also, I think there's a bug in what I put into MCIPv4.0, so I'm working to fix that now. If you have the tipping bucket on, you'll need different code in other parts of MCIP (not related to PRESF).

** NOTE: When I tested this with MCIP 3.6, met was from WRF v3.0, so no tipping bucket**

In mcoutcom_mod.f90:

1. Increase parameter MC3INDEX by 1.

2. Define a REAL, POINTER for "pressf_c" and "pressf_b" (follow x3htf_c and x3htf_b).

3. Add "PRESF" to the array for MC3VNAME.

4. Add units (Pa) to the array for MC3UNITS.

5. Add a text description to the array MC3VDESC.

In alloc_ctm.f90

1. Associate the pointers for PRESSF_C and PRESSF_B under the

allocate statements for MC3 and MB3. Again, follow x3htf_c and x3htf_b.

In dealloc_ctm.f90

1. Nullify the pointers for PRESSF_C and PRESSF_B just before the deallocate statements for MC3 and MB3. Follow x3htf_c and x3htf_b.

In metcro.f90

1. Change the allocate size for DUMARAY0 so that the fourth dimension is changed from "4+" to "5+" in three places.

2. In the triple-loop that fills DUMARAY0, fill dumaray0(c,r,k,5)

with xpresf(c,r,k).

3. Change the filling of XWWIND into the 6th element of DUMARAY0.

4. Change the filling of XTKE into the 6+iwout element of DUMARAY0.

5. Add a call to collapx for xpresf; follow the entry for

xdensaf.

6. Add an entry to fill pressf_c from xpresf; follow x3htm_c. (Note that xpresf is indexed 0:NLAYS rather than 1:NLAYS+1.)

7. Four times: Add an entry to fill pressf_b from xpresf; follow

x3htm_b. (One time for each of the four boundaries.)

8. Refill xpresf with elements from DUMARAY0 using the fourth

dimension of 5.

9. Change the links on the fourth dimension of DUMARAY0 for xwwind and xtke from 5 to 6 and 5+iwout, respectively (toward bottom of code).