=============
20160104:

1. Updated Higgins precipitation data through 20151231.

2. December report.

3. Made maps to compare MERRA2 and Princeton forcing. Programs created:
rain_monthly.gs
tair_monthly.gs
m2_m0001s_02_rain.ctl
m2_m0001s_02_tair.ctl
p0007s_66_rain.ctl
p0007s_66_tair.ctl
compare_rain_maxFPAR_princeton_merra2.m
compare_tair_maxFPAR_princeton_merra2.m


=============
20160105:

1. Continued from 3 on 20160104. Programs created:
SWdown_monthly.gs
LWdown_monthly.gs
m2_m0001s_02_SWdown.ctl
m2_m0001s_02_LWdown.ctl
p0007s_66_SWdown.ctl
p0007s_66_LWdown.ctl
compare_SWdown_maxFPAR_princeton_merra2.m
compare_LWdown_maxFPAR_princeton_merra2.m
compare_rain_maxFPAR_princeton_merra2_lead1.m

Checked all the programs in 3 on 20160104 and in 1 today to make sure that they are correct.

2. Met with Randy to talk about the differences in the MERRA2 and Princeton forcing, and my work plans for 2016 (see below).

(1) (Short term) FPAR project: using the combination of MERRA2 + Princeton forcing in the offline runs 
(2) (Short term) Work with Ben and Sarith to investigate the GCM-CN restart issue
(3) (Long term) Upgrade the CN model from CLM4 to CLM4.5. This should be my priority in the long term plans. 
(4) (Long term) Learn CVS
(5) (Long term) Work on the paper with Jim Collatz

3. GMAO carbon meeting.

4. Processed GMAO GEOS5 seasonal forecast for January.

=============
20160106:

1. About the GCM-CN issue. Reviewed the emails from/to Sarith, Greg and Ben and read GEOS_CatchGridComp.F90. Talked to Sarith. He said the catchment model is supposed to remain the same but now there is non-zero diff. He has to figure out why this happens and make sure catchment is the same first before we can move on to next step.

2. Worked on downscaling Princeton forcing from 3-hourly to hourly. 

/discover/nobackup/projects/lis/MET_FORCING/PRINCETON/2001> ncdump -h prcp_3hourly_2001-2001.nc 
netcdf prcp_3hourly_2001-2001 {
dimensions:
	longitude = 360 ;
	latitude = 180 ;
	z = 1 ;
	time = UNLIMITED ; // (2920 currently)
variables:
	float longitude(longitude) ;
		longitude:units = "degrees_east" ;
		longitude:point_spacing = "even" ;
		longitude:modulo = " " ;
	float latitude(latitude) ;
		latitude:units = "degrees_north" ;
		latitude:point_spacing = "even" ;
	float z(z) ;
		z:units = "level" ;
		z:positive = "up" ;
	float time(time) ;
		time:units = "hours since 2001-01-01 00:00:00" ;
		time:time_origin = "01-JAN-2001:00:00:00" ;
	float prcp(time, z, latitude, longitude) ;
		prcp:source = "Reanalysis daily precipitation corrected for spurious rain day frequencies, downscaled in space to 1.0deg and time to 3-hourly, scaled to match the CRU TS3.0 monthly dataset" ;
		prcp:name = "prcp" ;
		prcp:title = "3-hourly bias corrected precipitation" ;
		prcp:date = "01/01/01" ;
		prcp:time = "00:00" ;
		prcp:long_name = "Precipitation" ;
		prcp:units = "kg m-2 s-1" ;
		prcp:missing_value = 2.e+20f ;
		prcp:_FillValue = 2.e+20f ;
		prcp:valid_min = 0.f ;
		prcp:valid_max = 0.005674095f ;

// global attributes:
		:history = "Thu Mar 27 15:14:48 EDT 2014: created by JS using convert2alma.sh" ;
		:title = "Princeton University Hydroclimatology Group Bias Corrected 59-yr (1948-2006) Meteorological Forcing Dataset" ;
		:institution = "Princeton University" ;
		:contact = "Justin Sheffield (justin@princeton.edu)" ;
		:source = "Forcings are a hybrid of NCEP/NCAR reanalysis and observations" ;
		:comment = "This is an update to the original version reported in Sheffield et al., J. Climate (2006). Updates include: i) extension to nearer real-time; ii) improved sampling procedure for correction of rain day statistics; iii) use of latest versions of CRU, SRB and TRMM products; iv) improved consistency between specific and relative humidity and air temperature. See Sheffield et al., J. Climate (2006) for details of the observations used and the bias correction and downscaling methodology." ;
}

Created program:
princeton_3hourly_to_hourly.m

3. Checked out the bma_H50_catchcn tag:

cvs co -r bma_H50_catchcn GEOSagcm

The /discover/nobackup/fzeng/bma_H50_catchcn/src/GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts/check_land_restarts.pro can check the restart files. 
However, I can't get it work. It keeps giving an error message "Variable in undefined: NCDF_ISNCDF".

4. My /discover/nobackup/fzeng/Heracles-4_0-CN/src/GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 allocated tlai twice (in L3983 and L4636). Deleted the first one (the one in L3983) to avoid future error. 

Compiled /discover/nobackup/fzeng/Heracles-4_0-CN/src. It looks to me that it compiled successfully.
Will set up a test run tomorrow.


=============
20160107:

1. Met with Sarith.

Sarith used the executable and gcm_setup from Heracles-5_0 (can choose to run catchcn when running gcm_setup) and the restart files from my h4cn-0000 run to set up a test run. It runs! So the problem should be in the catchcn_internal_rst he created. He will look into this.

The top part of my h4cn-0000/gcm_run.j is too old. It now should be:

#######################################################################
#                     Batch Parameters for Run Job
#######################################################################

#PBS -l walltime=12:00:00
#SBATCH --ntasks=28
#SBATCH --output=/discover/nobackup/fzeng/h4cn-0000/slurm-%j.out
#PBS -N h4cn-0000_RUN
#SBATCH --constraint=sp3
#SBATCH -A g0620

#######################################################################
#                  System Environment Variables
#######################################################################


In my h4cn-0000/CAP.rc:

BEG_DATE:     18910301 000000
END_DATE:     20150701 210000
JOB_SGMT:     00000090 000000
NUM_SGMT:     5

It tells us that, 
(1) GCM writes out the rst files every 90 days (JOB_SGMT);
(2) GCM processes everything (archive, plot, etc.) every 90 days/sgmt * 5 sgmt = 450 days. 

2. Reversed the time stamp of /discover/nobackup/fzeng/h4cn-0001/GEOSgcm.x.old to it's original time stamp.

The /discover/nobackup/fzeng/h4cn-0001/GEOSgcm.x.old is the same as h4cn-0000/GEOSgcm.x:
/discover/nobackup/fzeng/h4cn-0001> diff GEOSgcm.x.old ../h4cn-0000/GEOSgcm.x
But it has a different time stamp now:
-rwxr-xr-x 1 fzeng g0620   90252069 2015-12-10 17:14 GEOSgcm.x.old  (h4cn-0001)
-rwxr-xr-x 1 fzeng g0620   90252069 2015-08-10 16:45 GEOSgcm.x      (h4cn-0000)
Because when Ben copied it back in 20151210 he didn't preserved the time stamp. 
Copied h4cn-0000/GEOSgcm.x back and name it GEOSgcm.x.old so that it has it's original time stamp:
/discover/nobackup/fzeng/h4cn-0001> cp -p ../h4cn-0000/GEOSgcm.x GEOSgcm.x.old
cp: overwrite `GEOSgcm.x.old'? y

3. Continued from 4 on 20160106. Set up a test run h4cn-0001_cont to run GCM from 20150701 to 20150801.

(1) Ran gcm_setup to create paths and the files (gcm_run.j etc.) below. This should be done after compiling Heracles-4_0-CN.

/discover/nobackup/fzeng/Heracles-4_0-CN/src/Applications/GEOSgcm_App> gcm_setup

Enter the Experiment ID:
h4cn-0001-cont
Enter a 1-line Experiment Description:
Continue h4cn-0001 from 20150701 to 20150801 for testing purpose
Enter an Experiment Source Tag for History (Default: Heracles-4_0):

Do you wish to CLONE an old experiment? (Default: NO or FALSE)

Enter the Atmospheric Horizontal Resolution code:
-----------------------------------------------------------
     Lat/Lon                     Cubed-Sphere
-----------------------------------------------------------
   b --  2  deg                c48  --  2   deg 
   c --  1  deg                c90  --  1   deg 
   d -- 1/2 deg                c180 -- 1/2  deg (56-km) 
   e -- 1/4 deg (35-km)        c360 -- 1/4  deg (28-km)  
                               c720 -- 1/8  deg (14-km) 
                               c1440 - 1/16 deg ( 7-km) 
 
b
Enter the Atmospheric Model Vertical Resolution: LM (Default: 72)

Do you wish to run the COUPLED Ocean/Sea-Ice Model? (Default: NO or FALSE)

Enter the Data_Ocean Horizontal Resolution code: o1 (1  -deg,  360x180  Reynolds) Default
                                                 o2 (1/4-deg, 1440x720  MERRA-2)
                                                 o3 (1/8-deg, 2880x1440 OSTIA)

Do you wish to run GOCART? (Default: NO or FALSE)

Enter the tag or directory (/filename) of the HISTORY.AGCM.rc.tmpl to use
(To use HISTORY.AGCM.rc.tmpl from current build, Type:  Current         )
-------------------------------------------------------------------------
Hit ENTER to use Default Tag/Location: (Current)

 
Enter Desired Location for the HOME Directory (to contain scripts and RC files)
Hit ENTER to use Default Location:
----------------------------------
Default:  /home/fzeng/geos5/h4cn-0001-cont


Enter Desired Location for the EXPERIMENT Directory (to contain model output and restart files)
Hit ENTER to use Default Location:
----------------------------------
Default:  /discover/nobackup/fzeng/h4cn-0001-cont

Enter Location for Build directory containing:  src/ Linux/ etc...
Hit ENTER to use Default Location:
----------------------------------
Default:  /discover/nobackup/fzeng/Heracles-4_0-CN

 
Current GROUPS: g0620 gmaoint m2val
Enter your GROUP ID for Current EXP: (Default: g0620)
-----------------------------------

sending incremental file list
GEOSgcm.x

sent 90263182 bytes  received 31 bytes  180526426.00 bytes/sec
total size is 90252069  speedup is 1.00
 
Creating gcm_run.j for Experiment: h4cn-0001-cont 
Creating gcm_post.j for Experiment: h4cn-0001-cont 
Creating gcm_archive.j for Experiment: h4cn-0001-cont 
Creating gcm_regress.j for Experiment: h4cn-0001-cont 
Creating gcm_convert.j for Experiment: h4cn-0001-cont 
Creating gcm_plot.tmpl for Experiment: h4cn-0001-cont 
Creating gcm_quickplot.csh for Experiment: h4cn-0001-cont 
Creating gcm_moveplot.j for Experiment: h4cn-0001-cont 
Creating gcm_stats.j for Experiment: h4cn-0001-cont 
Creating gcm_forecast.tmpl for Experiment: h4cn-0001-cont 
Creating gcm_forecast.setup for Experiment: h4cn-0001-cont 
Creating CAP.rc.tmpl for Experiment: h4cn-0001-cont 
Creating AGCM.rc.tmpl for Experiment: h4cn-0001-cont 
Creating HISTORY.rc.tmpl for Experiment: h4cn-0001-cont 
Creating ExtData.rc.tmpl for Experiment: h4cn-0001-cont 
 
Done!
-----
 
Build Directory: /discover/nobackup/fzeng/Heracles-4_0-CN
----------------
 
 
The following executable has been placed in your Experiment Directory:
----------------------------------------------------------------------
/discover/nobackup/fzeng/Heracles-4_0-CN/Linux/bin/GEOSgcm.x
 
 
You must now copy your Initial Conditions into: 
----------------------------------------------- 
/discover/nobackup/fzeng/h4cn-0001-cont


This U.S. Government resource is for authorized users only.  By accessing
this system you are consenting to complete monitoring with no expectation
of privacy.  Unauthorized access or use may subject you to disciplinary
action and criminal prosecution.

(2) Copied restart files from h4cn-0001:

/discover/nobackup/fzeng/h4cn-0001> cp -p *_rst ../h4cn-0001-cont/.

(3) Move AGCM.rc  CAP.rc*  ExtData.rc*  gcm_run.j*  HISTORY.rc in ~/geos5/h4cn-0001-cont to h4cn-0001-cont in nobackup and create a symbolic link in ~/geos5 to save space in home directory:

/home/fzeng/geos5/h4cn-0001-cont> mv * /discover/nobackup/fzeng/h4cn-0001-cont
/home/fzeng/geos5/h4cn-0001-cont> cd ..
/home/fzeng/geos5> rmdir h4cn-0001-cont
/home/fzeng/geos5> ln -s /discover/nobackup/fzeng/h4cn-0001-cont h4cn-0001-cont

(4)

/discover/nobackup/fzeng/h4cn-0001-cont> mv AGCM.rc AGCM.rc.orig
/discover/nobackup/fzeng/h4cn-0001-cont> mv CAP.rc CAP.rc.orig
/discover/nobackup/fzeng/h4cn-0001-cont> mv gcm_run.j gcm_run.j.orig
/discover/nobackup/fzeng/h4cn-0001-cont> mv HISTORY.rc HISTORY.rc.orig

/discover/nobackup/fzeng/h4cn-0001> cp -p AGCM.rc CAP.rc gcm_run.j HISTORY.rc ../h4cn-0001-cont/.
/discover/nobackup/fzeng/h4cn-0001> cp -p cap_restart ../h4cn-0001-cont/.

/discover/nobackup/fzeng/h4cn-0001-cont> nedit AGCM.rc CAP.rc gcm_run.j HISTORY.rc cap_restart &

(4a) In CAP.rc:
changed "END_DATE:     20150701 210000" to "END_DATE:     20150801 210000".

(4b) In gcm_run.j:
changed 

"#PBS -l walltime=12:00:00
#SBATCH --ntasks=112 --ntasks-per-node=28
#SBATCH --output=/discover/nobackup/fzeng/h4cn-0001/slurm-%j.out
#PBS -N h4cn-0001_RUN
#SBATCH --constraint=hasw
#SBATCH -A g0620"

to 

"#PBS -l walltime=12:00:00
#SBATCH --ntasks=28
#SBATCH --output=/discover/nobackup/fzeng/h4cn-0001-cont-cont/slurm-%j.out
#PBS -N h4cn-0001-cont-cont_RUN
#SBATCH --constraint=sp3
#SBATCH -A g0620"

as suggested by Sarith (see my note 1 today above).

Do a global change of "h4cn-0001" to "h4cn-0001-cont".

(4c) In HISTORY.rc:

changed

"EXPID:  h4cn-0001
EXPDSC: h4cn-0000_with_fire_factor_pt1
EXPSRC: Heracles-4_0" 

to 

"EXPID:  h4cn-0001-cont
EXPDSC: Continue_h4cn-0001_from_20150701_to_20150801_for_testing_purpose
EXPSRC: Heracles-4_0" 

(same as the top part of HISTORY.rc.orig).

(4d) Now cat_restart is: 20150701 210000

4. Continued from 2 on 20160106. Worked on downscaling Princeton forcing from 3-hourly to hourly.
princeton_3hourly_to_hourly.m
downscale_1to3_1D.m 
Both in ~fzeng/Catchment/matlab. 

Need to improve downscale_1to3_1D.m. It's currently too slow!

5. Met with Sarith. 

Looked into my h4cn-0000 catch_internal_rst and Sarith's catchcn_internal_rst. Sarith fixed two things in his catchcn_internal_rst:
(1) NDEP = NDEP*1.e-9 so that the unit of ndep is g N ...
(2) albedo = max(albedo, 1.e-6) to avoid 0 in albedo.

I will use Sarith's catch_internal_rst (not including CN rst) and append_CN.f90 to make catchcn_internal_rst and see if it works.

=============
20160108:

Met with Sarith the whole day on the catchcn_internal_rst issue.

=============
20160111:

Continued to work on the catchcn_internal_rst issue with Sarith. Sarith will use the CN_restart from my SMAP_EASEv2_M09/e0004s to create catchcn_internal_rst.

=============
20160112:

1. Continued to work on the catchcn_internal_rst issue with Sarith.

2. Continued to work on downscaling Princeton forcing from 3-hourly to hourly.
Improved function downscale_1to3_1D to save run time. 

Asked Qing about how to write out the netCDF file. She shared with me one of her scripts.
~qliu/projects/matlab_code/precip_correction/CPCU_MLAND/SMAP/write_CPCUcorr_G5_nc4.m 

=============
20160113:

1. Continued to work on downscaling Princeton forcing from 3-hourly to hourly.

In matlab:

>> ncdisp('/discover/nobackup/projects/lis/MET_FORCING/PRINCETON/2012/prcp_3hourly_2012-2012.nc')
Source:
           /discover/nobackup/projects/lis/MET_FORCING/PRINCETON/2012/prcp_3hourly_2012-2012.nc
Format:
           classic
Global Attributes:
           history     = 'Thu Mar 27 16:17:01 EDT 2014: created by JS using convert2alma.sh'
           title       = 'Princeton University Hydroclimatology Group Bias Corrected 59-yr (1948-2006) Meteorological Forcing Dataset'
           institution = 'Princeton University'
           contact     = 'Justin Sheffield (justin@princeton.edu)'
           source      = 'Forcings are a hybrid of NCEP/NCAR reanalysis and observations'
           comment     = 'This is an update to the original version reported in Sheffield et al., J. Climate (2006). Updates include: i) extension to nearer real-time; ii) improved sampling procedure for correction of rain day statistics; iii) use of latest versions of CRU, SRB and TRMM products; iv) improved consistency between specific and relative humidity and air temperature. See Sheffield et al., J. Climate (2006) for details of the observations used and the bias correction and downscaling methodology.'
Dimensions:
           longitude = 360
           latitude  = 180
           z         = 1
           time      = 2928  (UNLIMITED)
Variables:
    longitude
           Size:       360x1
           Dimensions: longitude
           Datatype:   single
           Attributes:
                       units         = 'degrees_east'
                       point_spacing = 'even'
                       modulo        = ' '
    latitude 
           Size:       180x1
           Dimensions: latitude
           Datatype:   single
           Attributes:
                       units         = 'degrees_north'
                       point_spacing = 'even'
    z        
           Size:       1x1
           Dimensions: z
           Datatype:   single
           Attributes:
                       units    = 'level'
                       positive = 'up'
    time     
           Size:       2928x1
           Dimensions: time
           Datatype:   single
           Attributes:
                       units       = 'hours since 2012-01-01 00:00:00'
                       time_origin = '01-JAN-2012:00:00:00'
    prcp     
           Size:       360x180x1x2928
           Dimensions: longitude,latitude,z,time
           Datatype:   single
           Attributes:
                       source        = 'Reanalysis daily precipitation corrected for spurious rain day frequencies, downscaled in space to 1.0deg and time to 3-hourly, scaled to match the CRU TS3.0 monthly dataset'
                       name          = 'prcp'
                       title         = '3-hourly bias corrected precipitation'
                       date          = '01/01/12'
                       time          = '00:00'
                       long_name     = 'Precipitation'
                       units         = 'kg m-2 s-1'
                       missing_value = 2.000000040081755e+20
                       _FillValue    = 2.000000040081755e+20
                       valid_min     = 0
                       valid_max     = 0.006042


=============
20160114:

1. Continued to work on downscaling Princeton forcing from 3-hourly to hourly.

Double checked the global attributes and variable attributes to see if other modifications are needed. 

In princeton_3hourly_to_hourly.m (~fzeng/Catchment/matlab/fpar/):
Since ChunkSize is [] for all 5 variables
    data_chunksize = finfo.Variables(i).ChunkSize;
    netcdf.defVarChunking(fout_id, varid(i), 'CHUNKED', data_chunksize); 
    
    % ChunkSize is [] for all 5 variables, causing error message "The length of the chunking parameter (0) does not match the variable's number of dimensions."
    netcdf.defVarChunking(fout_id, varid(i), 'CHUNKED', dims_len(i));
    
    
=============
20160119:

1. Continued to work on downscaling Princeton forcing from 3-hourly to hourly.

Qing's script is for netCDF4 and the Princeton forcing is in netCDF3. Such difference caused some error messages. Resolved most of them. 

2. Melanie asked for these data to test out a new biomass burning algorithm (see her email today).

Temperature: Tair (2m air temperature?)
RH: we have specific humidity and surface pressure, need water vapor to convert to relative humidity
Precipitation: RainfC, Rainf, Snowf (separate)
LAI 

Checked my notes on 20150727: the data I sent her was from p0007s_37. Retrieved p0007s_37 from archive.  
  
    
=============
20160120:

1. Extracted the forcing data for Melanie. Scripts/files created:

extract_forcing.gs
read_forcing.f90
p0007s_37_forcing.gdat
p0007s_37_forcing.ctl  

Variables extracted:

Tair       2m air temperature (K) [forcing]
Qair      2m specific humidity (kg/kg) [forcing]
Psurf     surface pressure (Pa) [forcing]
RainfC   convective rainfall (mm/d) [forcing]
Rainf     largescale rainfall (mm/d) [forcing]
Snowf   snowfall (mm/d) [forcing]
tsurf      surface temperature (K)
lai          leaf-area index

See my email to her today.

2. Continued to work on downscaling Princeton forcing from 3-hourly to hourly. 

Wrote check_hourly_princeton.m to check the data and global/variable attributes in the hourly netCDF files created by princeton_3hourly_to_hourly.m.


=============
20160121:

1. Finished check_hourly_princeton.m and ran it on the 2004 precipitation data. The difference between the hourly data generated and the original 3hourly data is generally <0.1% or 4 mm/yr. 

Modified princeton_3hourly_to_hourly.m so that it also works for air temperature (tas). 
(1) Checked global attributes.
(2) Checked variable attributes.

Produced 2000 to 2012 hourly Princeton precipitation data (~29 GB) in netCDF files using princeton_3hourly_to_hourly.m. Files are currently located in /discover/nobackup/fzeng/princeton_hourly. 

2. Modified clsm_ensdrv_force_routines.F90 to use MERRA2 forcing and Princeton precipitation.

Tried parallel build: No need to ssh to discover-sp3.
/discover/nobackup/fzeng/LDASsa_m3-15_2-CN> setenv ESMADIR /discover/nobackup/fzeng/LDASsa_m3-15_2-CN
/discover/nobackup/fzeng/LDASsa_m3-15_2-CN> source $ESMADIR/src/g5_modules
/discover/nobackup/fzeng/LDASsa_m3-15_2-CN/src> ./parallel_build.csh

> q (to check the progress, just as when we check the run progress)

The log files are 
/discover/nobackup/fzeng/LDASsa_m3-15_2-CN/src/BUILD_LOG_DIR/LOG and
/discover/nobackup/fzeng/LDASsa_m3-15_2-CN/src/parallel_build.o*******.
The latter includes everything in the former. Can just look at the former since it already includes everything about the compiling process that we want to know.   

If the build was not successful, just do this again:
/discover/nobackup/fzeng/LDASsa_m3-15_2-CN/src> ./parallel_build.csh 

Build successful. Executable is 
/discover/nobackup/fzeng/LDASsa_m3-15_2-CN/Linux/bin/LDASsa_assim_mpi.x

=============
20160122: 

1. Copy the executable to /discover/nobackup/fzeng/LDASsa_m3-15_2-CN/exec and do a test run interactively. 

/discover/nobackup/fzeng/LDASsa_m3-15_2-CN/Linux/bin> cp -p LDASsa_assim_mpi.x ../../exec/LDASsa_mpi_73q_e0007s_old_merra2_Princeton_precip.x

Make sure the CN restart and land restart match:
/discover/nobackup/fzeng/Catchment/merra2/m0002s> ls -l CN_restart 
-rw-r--r-- 1 fzeng g0620 1533248288 2015-10-27 05:11 CN_restart
/discover/nobackup/fzeng/Catchment/merra2/m0002s> ls -l RUN/rs/ens0000/Y2001/M01/
total 40576
-rw-r--r-- 1 fzeng g0620 41519520 2015-10-27 05:11 m0002.ens0000.catch_ldas_rst.20010101_0000z.bin
    
/discover/nobackup/fzeng/Catchment/merra2/m0002s> nedit lenkf.j &
changed STOPYEAR from 2015 to 2013 because Princeton forcing is up to 2012. 

Interactive run not successful:

/discover/nobackup/fzeng/Catchment/merra2/m0002s> interative.py -A sp3 -n 96 -a g0620 -X -t 1:00:00
interative.py: Command not found.
/discover/nobackup/fzeng/Catchment/merra2/m0002s> ~fzeng/bin/interative.py -A sp3 -n 96 -a g0620 -X -t 1:00:00

CORRECT>~fzeng/bin/interactive.py -A sp3 -n 96 -a g0620 -X -t 1:00:00 (y|n|e|a)? yes
Running: /usr/bin/ssh -XYqt -p2255 dali12 salloc --ntasks=96 --constraint=sp3 --time=1:00:00 --account=g0620 --mail-type=BEGIN
Traceback (most recent call last):
  File "/home/fzeng/bin/interactive.py", line 163, in <module>
    main()
  File "/home/fzeng/bin/interactive.py", line 129, in main
    sp.check_call(cmd.split())
  File "/usr/local/other/SLES11/SIVO-PyD/1.10.0/lib/python2.7/subprocess.py", line 504, in check_call
    raise CalledProcessError(retcode, cmd)
subprocess.CalledProcessError: Command '['/usr/bin/ssh', '-XYqt', '-p2255', 'dali12', 'salloc', '--ntasks=96', '--constraint=sp3', '--time=1:00:00', '--account=g0620', '--mail-type=BEGIN']' returned non-zero exit status 255

/discover/nobackup/fzeng/Catchment/merra2/m0002s> ./lenkf.j
./lenkf.j: Permission denied. 

2. Went back to the usual way but run crashed:

/discover/nobackup/fzeng/Catchment/merra2/m0002s> sbatch lenkf.j

 LDAS ERROR (3000) from read_ens_upd_inputs: executable was built for assimilati
 on but update_type=0
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD 
with errorcode 3000.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------

3. Compiled in the usual way and ran in the usual way.

Got this error message:
 opening/discover/nobackup/fzeng/princeton_hourly//dswrf_hourly_2001-2001.nc
 LDAS ERROR (3000) from get_Princeton_netcdf: error opening netcdf file
 
This is because so far I have only downscaled precipitation. Need to downscale the other 6 variables as well before I can run the model with the current code. 

4. Modified clsm_ensdrv_force_routines.F90 to just read in one Princeton forcing variable at each run. Compiled and did a test run. Worked. Submitted the full run.

5. Ran grid_restore to process the output and checked the output in GrADS.

~/Catchment > grid_restore_merra2 0002 s
So far I've got output up to Nov 2002. 

In GrADS, compared merra2/m0002s to p0005s_63, looked at jul2001 and jul2002:
RainfC: same, both are 0 everywhere. This is good, but
Rainf and Snowf are different between the two runs. Difference in Rainf is large. Difference in Snowf is small. 
Rainf+Snowf are different between the two. Therefore, the difference couldn't be due to the difference in Tair. Something must be wrong!!

In GrADS, compared m0002s to m0001s both in merra2/, looked at jul2001 and jul2002:
Tair, Psurf, LWdown, SWdown and Wind: same. Good.
Qair: same except in a few grid cells where the difference is below 4e-05. Is this ok?

Differences in Rainf and Snowf between m0002s and m0001s are smaller than those between m0002s and p0005s.  

=============
20160127: 

1. Investigated why the Rainf and Snowf in merra2/m0002s are different from those in p0005s.

Questions:

In call get_GEOS (L318-407) in clsm_ensdrv_force_routines.F90, there are met_force_obs_tile_old and met_force_obs_tile_new. How are the hourly MERRA2 forcing used in the model run exactly? Is that forcing t and t+1 (or t-1 and t) are interpolated to obtain the meteorology between t and t+1? Which routine does this? The Princeton precip was given to met_force_obs_tile_old? What about met_force_obs_tile_new? Is it still using MERRA2 forcing?

When using MERRA2 (and MERRA as well), there are met_force_obs_tile_old and met_force_obs_tile_new. met_force_obs_tile_old has valid data while met_force_obs_tile_new is nodata_generic. (see L310-372 in the original clsm_ensdrv_force_routines.F90.m3-15_2).

When using Princeton forcing, the valid data is given to met_force_obs_tile_new. (see L282-292 in the original clsm_ensdrv_force_routines.F90.m3-15_2). What's in met_force_obs_tile_old in this case?

2. Another problem. When the model run reached the end year 2013, it asked for /discover/nobackup/fzeng/princeton_hourly//prcp_hourly_2013-2013.nc which does not exist. Therefore, the run failed at the end of the run in 2013. I should have set a symbolic link to make prcp_hourly_2013-2013.nc point to prcp_hourly_2001-2001.nc. Done.

3. Made plots showing Jul 2015 and Dec 2015 GPP, NPP, NEE, FPAR and latent heat of the SMAP run, the control run and the difference between the tow. These are for Randy's presentation at the SMAP meeting next week. 


=============
20160128: 

1. Continued 3 on 20160127.

2. Organized emails. Mail box almost full.

3. Continued to investigate why the Rainf and Snowf in merra2/m0002s are different from those in p0005s. Talked to Qing and she suggested/asked these:

(1) Does met_force_obs_tile_new in Princeton have the same structure as that in get_GEOS which reads MERRA2?
I think so. met_force_obs_tile_new is declared at the beginning of subroutine get_forcing in clsm_ensdrv_force_routines.F90 as:

    type(met_force_type), dimension(N_catd), intent(out) :: &
         met_force_obs_tile_new
It's then used for both subroutines get_Princeton_netcdf and get_GEOS.

And in driver_types.F90, met_force_type is declared as:

  type :: met_force_type
     real :: Tair                 ! air temperature at RefH                 [K]
     real :: Qair                 ! specific humidity at RefH               [kg/kg]
     real :: Psurf                ! surface pressure                        [Pa]
     real :: Rainf_C              ! convective rainfall                     [kg/m2/s]
     real :: Rainf                ! total rainfall                          [kg/m2/s]
     real :: Snowf                ! total snowfall                          [kg/m2/s]
     real :: LWdown               ! downward longwave radiation             [W/m2]
     real :: SWdown               ! downward shortwave radiation            [W/m2]
     real :: SWnet                ! downward net shortwave radiation        [W/m2]
     real :: PARdrct              ! Photosynth. Active Radiation (direct)   [W/m2]
     real :: PARdffs              ! Photosynth. Active Radiation (diffuse)  [W/m2]
     real :: Wind                 ! wind speed at RefH                      [m/s]
     real :: RefH                 ! reference height for Tair, Qair, Wind   [m]
  end type met_force_type
        
(2) Try reading Princeton precip inside get_GEOS()
(3) Print out the precip of the 1st time step of some tile to verify

(4) Could the difference in spatial resolution of Princeton and MERRA2 have caused the problem with the current code that I modified?

Princeton: 180x360 (1x1deg)
MERRA2: 361x576 (0.5x0.625deg)
But after the data is read in, it becomes tile basis. So the difference in spatial resolution should not be an issue, right?

(5) p0005s and merra2/m0002s use different tile coordinate systems. Could this have caused the difference in precip output?

4. To see if the difference in tile coordinate system has caused the difference in precip output ((5) above):

Set up a run p0005st ("t" means "test") same as p0005s but use the tile coordinate system in merra2/m0002s:
/discover/nobackup/fzeng/Catchment/princeton> mkdir -p p0005st/RUN/rs/ens0000/Y2001/M01
/discover/nobackup/fzeng/Catchment/princeton/p0005s_62> cp -p CN_restart lenkf.j ../p0005st/.
/discover/nobackup/fzeng/Catchment/princeton/p0005s_62> cp -p RUN/rs/ens0000/Y2013/M01/p0005.ens0000.catch_ldas_rst.20130101_0000z.bin ../p0005st/RUN/rs/ens0000/Y2001/p0005.ens0000.catch_ldas_rst.20010101_0000z.bin

When modifying lenkf.j, realized that the restart files should be in the same tile coordinate system as in merra2/m0002s. Therefore, copied the restart files from m0008s_31 as what's done for merra2/m0002s (see notes on 20151224):

/discover/nobackup/fzeng/Catchment/merra/m0008s_31> cp -p CN_restart /discover/nobackup/fzeng/Catchment/princeton/p0005st/.
/discover/nobackup/fzeng/Catchment/merra/m0008s_31> cp -p RUN/rs/ens0000/Y2015/M01/m0008.ens0000.catch_ldas_rst.20150101_0000z.bin /discover/nobackup/fzeng/Catchment/princeton/p0005st/RUN/rs/ens0000/Y2001/M01/p0005.ens0000.catch_ldas_rst.20010101_0000z.bin

In lenkf.j: 
(a) Changed all 0005s to 0005st
(b) Changed all the command lines related to tile systems to reflect the "new" tile system
(c) Created a new namelist /home/fzeng/Catchment/LDASsa_m3-15/p0000s/LDASsa_inputs_driver_p0005st.nml to reflect the "new" tile system

Looked at the (Rainf+Snowf) for Jan 2001 in GrADS:
Difference between p0005st and p0005s: within -0.2 to 0.2 mm/d everywhere;
Difference between p0005st and merra2/m0002s: within -0.2 and 0.2 mm/d in most places, but in some places the difference is up to 1.2 mm/d. 

Therefore, the difference in tile coordinate system only caused a very small difference in the precip output and is not the major reason.


=============
20160129: 

1. Continued to investigate why the Rainf and Snowf in merra2/m0002s are different from those in p0005s.

Princeton forcing does not have all the 13 variables in met_force_type defined in driver_types.F90. How does the model deal with those variables not available in Princeton forcing?  

Tair: Princeton     
Qair: Princeton     
Psurf: Princeton    
Rainf_C: set as 0 in subroutine get_Princeton_netcdf
Rainf: Princeton    
Snowf: Princeton    
LWdown: Princeton   
SWdown: Princeton   
SWnet: set as nodata_generic at the beginning of subroutine get_forcing    
PARdrct: set as nodata_generic at the beginning of subroutine get_forcing  
PARdffs: set as nodata_generic at the beginning of subroutine get_forcing  
Wind: Princeton     
RefH: set as DEFAULT_REFH at the beginning of subroutine get_forcing     

Questions:
(1) Is it ok that met_force_obs_tile_new from the subroutine get_Princeton_netcdf that I modified has only three variables?
Read subroutine get_Princeton_netcdf. Should be fine. 
(2) How does it decide whether the precip goes to Rainf or Snowf in the subroutine get_Princeton_netcdf that I modified?
After "call get_GEOS" to read in MERRA2 data, all 13 variables in met_force_obs_tile_new have valid data. So met_force_obs_tile_new%Tair (i.e. met_force_new%Tair inside subroutine get_Princeton_netcdf) is from MERRA2. Therefore,  if MERRA2 Tair is below freezing point, precip goes to Snowf. Otherwise to Rainf. 

Found a stipid mistake in my code under subroutine get_Princeton_netcdf that I modified:

    do k=1, N_catd
       
       met_force_new(k)%Rainf_C = 0.
       
       if ( met_force_new(k)%Tair < Tzero ) then
          met_force_new(k)%Rainf = 0.
!          met_force_new(k)%Snowf = force_array(k,3)  ! Comment this out to just use Princeton precipitation, F Zeng, Jan 2016
          met_force_new(k)%Snowf = force_array(k,1)
       else
!          met_force_new(k)%Rainf = force_array(k,3)  ! Comment this out to just use Princeton precipitation, F Zeng, Jan 2016
          met_force_new(k)%Snowf = force_array(k,1)
          met_force_new(k)%Snowf = 0.
       endif
       
    enddo

The above block should be:

    do k=1, N_catd
       
       met_force_new(k)%Rainf_C = 0.
       
       if ( met_force_new(k)%Tair < Tzero ) then
          met_force_new(k)%Rainf = 0.
!          met_force_new(k)%Snowf = force_array(k,3)  ! Comment this out to just use Princeton precipitation, F Zeng, Jan 2016
          met_force_new(k)%Snowf = force_array(k,1)
       else
!          met_force_new(k)%Rainf = force_array(k,3)  ! Comment this out to just use Princeton precipitation, F Zeng, Jan 2016
          met_force_new(k)%Rainf = force_array(k,1)
          met_force_new(k)%Snowf = 0.
       endif
       
    enddo   
    
Fixed this bug. Compiled. Set up a debug run: same as m0001s_01, got the restart files from m0008s_31.

/discover/nobackup/fzeng/Catchment/merra2 > mv m0002s m0002s_wrong    
/discover/nobackup/fzeng/Catchment/merra2 > mkdir -p m0002s/RUN/rs/ens0000/Y2001/M01
/discover/nobackup/fzeng/Catchment/merra/m0008s_31 > cp -p CN_restart ../../merra2/m0002s/.
/discover/nobackup/fzeng/Catchment/merra/m0008s_31 > cp -p RUN/rs/ens0000/Y2015/M01/m0008.ens0000.catch_ldas_rst.20150101_0000z.bin ../../merra2/m0002s/RUN/rs/ens0000/Y2001/M01/m0002.ens0000.catch_ldas_rst.20010101_0000z.bin
/discover/nobackup/fzeng/Catchment/merra2/m0002s > cp -p ../m0002s_wrong/lenkf.j .

Make sure that CN restart matches land restart:
dali11:/discover/nobackup/fzeng/Catchment/merra2/m0002s > ls -l CN_restart 
-rw-r--r-- 1 fzeng g0620 1533248288 2015-10-27 05:11 CN_restart
dali11:/discover/nobackup/fzeng/Catchment/merra2/m0002s > ls -l RUN/rs/ens0000/Y2001/M01
total 40576
-rw-r--r-- 1 fzeng g0620 41519520 2015-10-27 05:11 m0002.ens0000.catch_ldas_rst.20010101_0000z.bin

Submitted the job. It finished Jan-May 2001 within 1 hour. 
Looked at the Jan and May 2001 data in GrADS:
Jan 2001: 
m0002s vs. p0005s (or p0005st): difference in Rainf+Snowf within 0.01 mm/d in most places;
m0002s vs. m0001s: 0 diff for Tair, Psurf, SWdown and Wind; within 1e-5 (kg/kg) almost everywhere for Qair; within 0.1 W/m2 in most places for LWdown, but up to ~2 W/m2 in the high latitudes and in Tibetan plateau. 
May 2001: 
m0002s vs. p0005s: difference in Rainf+Snowf within 0.01 mm/d in most places;
m0002s vs. m0001s: 0 diff for Tair, Psurf, SWdown and Wind; within 1e-5 (kg/kg) almost everywhere for Qair; within 0.1 W/m2 in most places for LWdown.     
    
Sarith helped me to look at the differences above and the code that I modified. We found the problem is in this block of code that I modified in clsm_ensdrv_force_routines.F90:

       ! Replace precip with Princeton precip for sensitivity analysis, fzeng, Jan 2016
       !--------------------------------------------------------------------------------
       
       if (logit) write (logunit,*) 'replace MERRA2 precip with Princeton precip'
       
       met_path_hourly_Princeton = '/discover/nobackup/fzeng/princeton_hourly/'
       
       call get_Princeton_netcdf(  date_time_tmp, met_path_hourly_Princeton, N_catd, tile_coord, &
            met_force_obs_tile_new, nodata_forcing_Princeton)
       
       ! check for nodata values and unphysical values
       ! (check here, not outside "if" block, because of GEOSgcm case)
       
       call check_forcing_nodata( N_catd, tile_coord, nodata_forcing_Princeton, &
            met_force_obs_tile_new )
       
       call repair_forcing( N_catd, met_force_obs_tile_new, &
            echo=.true., tile_coord=tile_coord )

While for MERRA2, call repair_forcing is different:

       ! call repair_forcing with switch "unlimited_Qair=.true." for GEOS5 forcing
       ! (default is to limit Qair so that it does not exceed Qair_sat)
       ! reichle+qliu,  8 Oct 2008
       ! 
       ! likewise for "unlimited_LWdown=.true."
       ! reichle, 11 Feb 2009       
       
       call repair_forcing( N_catd, met_force_obs_tile_new, &
            echo=.true., tile_coord=tile_coord,             &
            unlimited_Qair=.true., unlimited_LWdown=.true. )

That's why there are differences in Qair and LWdown between m0002s and m0001s.

Need to read subroutine repair_forcing in clsm_ensdrv_force_routines.F90 (L4562-).   
    
Fixed the bug. Compiled. Set up a debug run:

/discover/nobackup/fzeng/Catchment/merra2 > mv m0002s m0002s_wrong
/discover/nobackup/fzeng/Catchment/merra2 > mkdir -p m0002s/RUN/rs/ens0000/Y2001/M01
/discover/nobackup/fzeng/Catchment/merra/m0008s_31 > cp -p CN_restart ../../merra2/m0002s/.
/discover/nobackup/fzeng/Catchment/merra/m0008s_31 > cp -p RUN/rs/ens0000/Y2015/M01/m0008.ens0000.catch_ldas_rst.20150101_0000z.bin ../../merra2/m0002s/RUN/rs/ens0000/Y2001/M01/m0002.ens0000.catch_ldas_rst.20010101_0000z.bin   
/discover/nobackup/fzeng/Catchment/merra2/m0002s > cp -p ../m0002s_wrong/lenkf.j .
/discover/nobackup/fzeng/Catchment/merra2/m0002s > sbatch lenkf.j

It finished Jan-Feb 2001 in half an hour. 
Looked at the data in GrADS, for both Jan and Feb 2001.
m0002s vs. p0005s (or p0005st): difference in Rainf+Snowf within 0.01 mm/d in most places;
m0002s vs. m0001s: 0 diff for all other variables (Tair, Psurf, SWdown, LWdown, Qair and Wind)
Great! This is correct now.

Cancel the debug run and set up a full run:
/discover/nobackup/fzeng/Catchment/merra2 > mv m0002s m0002s_wrong    
/discover/nobackup/fzeng/Catchment/merra2 > mkdir -p m0002s/RUN/rs/ens0000/Y2001/M01
/discover/nobackup/fzeng/Catchment/merra/m0008s_31 > cp -p CN_restart ../../merra2/m0002s/.
/discover/nobackup/fzeng/Catchment/merra/m0008s_31 > cp -p RUN/rs/ens0000/Y2015/M01/m0008.ens0000.catch_ldas_rst.20150101_0000z.bin ../../merra2/m0002s/RUN/rs/ens0000/Y2001/M01/m0002.ens0000.catch_ldas_rst.20010101_0000z.bin
/discover/nobackup/fzeng/Catchment/merra2/m0002s > cp -p ../m0002s_wrong/lenkf.j .

2. e0004s and e0005s FPAR vs. GIMMS FPAR:

Randy asked me to correlate spatially the e0004s FPAR and the e0005s FPAR against GIMMS FPAR for Africa. First I need to get the GIMMS NDVI data from Jorge, and then I need to convert the NDVI to FPAR. Talked to Jorge. He told me that Al has obtained 8km NDVI data up to Sep 2015 and has aggregated it to 0.5deg NDVI. 

Obtained from Al 0.5deg NDVI for Jan-Sep 2015.
Obtained from Jim the updated script that converts NDVI to FPAR.     
    
    
    
    
    
    
