========== 20170103: 1. Check out GEOS5 fresh because Sarith has updated it (including a tiny bug fix) since last time I checked it out on 20161129. (1) Check out: cd /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5 cvs co -r SM-Heracles-5_4_p3-M3_V24_C05 GEOSagcm mv GEOSagcm Heracles-5_4_p3-M3_V24_C05 xdiff Heracles-5_4_p3-M3_V24_C05/src/GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/compute_rc.F90 /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/src/Components/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/compute_rc.F90.clm4 & They are identical. Good. (2) Tag (so that I know which version of this GEOS5 tag I used): cd /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/src cvs tag FZ_ctrl_clm4 (3) Build the executable: cd /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/src source g5_modules ./parallel_build.csh 2. Saved a copy of the restart files created on 20161205. This can be used in CLM4.5 test runs. cd /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/ mkdir restarts cd test_clm4/restarts cp -p restarts.e20001230_21z.tar ../../restarts/. 3. Sarith said: If I need to repeat an GEOS5 experiment, I just need to copy the executable from the experiment directory back to Linux/bin and then run gcm_setup from Applications/GEOSgcm_App to set up another experiment. ========== 20170104: 1. Set up an experiment: (1) Run gcm_setup to set up an experiment: cd /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/src/Applications/GEOSgcm_App ./gcm_setup Enter the Experiment ID: clm4 Enter a 1-line Experiment Description: Test run of original clm4 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) c90 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) o2 Enter the choice of Land Surface Model: use 1 (Default: Catchment Model), 2 (CatchmentCN Model) 2 Do you wish to run the Runoff Routing Model? (Default: NO or FALSE) Do you wish to run GOCART with Actual or Climatological Aerosols? (Enter: A (Default) or C) Enter the GOCART Emission Files to use: MERRA2 (Default), PIESA, CMIP, NR, MERRA2-DD or OPS: 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: /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4 Enter Desired Location for the EXPERIMENT Directory (to contain model output and restart files) Hit ENTER to use Default Location: ---------------------------------- Default: /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4 Enter Location for Build directory containing: src/ Linux/ etc... Hit ENTER to use Default Location: ---------------------------------- Default: /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05 Current GROUPS: g0620 gmaoint m2val Enter your GROUP ID for Current EXP: (Default: g0620) ----------------------------------- sending incremental file list GEOSgcm.x sent 100153687 bytes received 31 bytes 66769145.33 bytes/sec total size is 100141366 speedup is 1.00 Creating gcm_run.j for Experiment: clm4 Creating gcm_post.j for Experiment: clm4 Creating gcm_archive.j for Experiment: clm4 Creating gcm_regress.j for Experiment: clm4 Creating gcm_convert.j for Experiment: clm4 Creating gcm_plot.tmpl for Experiment: clm4 Creating gcm_quickplot.csh for Experiment: clm4 Creating gcm_moveplot.j for Experiment: clm4 Creating gcm_stats.j for Experiment: clm4 Creating gcm_forecast.tmpl for Experiment: clm4 Creating gcm_forecast.setup for Experiment: clm4 Creating CAP.rc.tmpl for Experiment: clm4 Creating AGCM.rc.tmpl for Experiment: clm4 Creating HISTORY.rc.tmpl for Experiment: clm4 Creating fvcore_layout.rc for Experiment: clm4 Done! ----- Build Directory: /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05 ---------------- The following executable has been placed in your Experiment Directory: ---------------------------------------------------------------------- /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/Linux/bin/GEOSgcm.x You must now copy your Initial Conditions into: ----------------------------------------------- /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4 Copying Sandbox Source Code (including non-committed files) into /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4/src/GEOSagcm/src ... Enter CVS Tag to COMMIT this Experiment (Default: clm4__fzeng) Enter q or quit to QUIT or SKIP the CVS COMMIT: q 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. Could not chdir to home directory /home/fzeng: No such file or directory Non-Committed Sandbox Files: --------------------------------- Applications/GEOSgcm_App/.AGCM_VERSION Tarring Experiment Source Code into Single File ... (2) Copy and rename restart files: cd /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/ cp restarts/restarts.e20001230_21z.tar clm4/. cd clm4 tar -xf restarts.e20001230_21z.tar /bin/rm restarts.e20001230_21z.tar nedit move & move: ~~~~~~~~~~~~~~~~~~~~~~ mv test_clm4.agcm_import_rst.e20001230_21z.Heracles-5_4_p3.Heracles-4_3_MERRA-3_CF0090x6C_DE1440xPE0720 agcm_import_rst mv test_clm4.catchcn_internal_rst.e20001230_21z.Heracles-5_4_p3.Heracles-4_3_MERRA-3_CF0090x6C_DE1440xPE0720 catchcn_internal_rst mv test_clm4.fvcore_internal_rst.e20001230_21z.Heracles-5_4_p3.Heracles-4_3_MERRA-3_CF0090x6C_DE1440xPE0720 fvcore_internal_rst mv test_clm4.gocart_internal_rst.e20001230_21z.Heracles-5_4_p3.Heracles-4_3_MERRA-3_CF0090x6C_DE1440xPE0720 gocart_internal_rst mv test_clm4.lake_internal_rst.e20001230_21z.Heracles-5_4_p3.Heracles-4_3_MERRA-3_CF0090x6C_DE1440xPE0720 lake_internal_rst mv test_clm4.landice_internal_rst.e20001230_21z.Heracles-5_4_p3.Heracles-4_3_MERRA-3_CF0090x6C_DE1440xPE0720 landice_internal_rst mv test_clm4.moist_internal_rst.e20001230_21z.Heracles-5_4_p3.Heracles-4_3_MERRA-3_CF0090x6C_DE1440xPE0720 moist_internal_rst mv test_clm4.pchem_internal_rst.e20001230_21z.Heracles-5_4_p3.Heracles-4_3_MERRA-3_CF0090x6C_DE1440xPE0720 pchem_internal_rst mv test_clm4.saltwater_internal_rst.e20001230_21z.Heracles-5_4_p3.Heracles-4_3_MERRA-3_CF0090x6C_DE1440xPE0720 saltwater_internal_rst ~~~~~~~~~~~~~~~~~~~~~~ chmod 755 move move (3) Modify some files: Create cap_restart (which has only this line below) 20001230 210000 Edit CAP.rc END_DATE: 20010201 210000 (run 1 month for testing first) Edit HISTORY.rc In COLLECTIONS, comment out everything other than "geosgcm_surf" [Line 57] geosgcm_surf.frequency: 240000, [Lines 206-242] Uncomment geosgcm_surf collection for carbon variables (e.g., CNLAI, ICESOI) Edit gcm_run.j [line 7] #PBS -l walltime=06:00:00 [Line 161] setenv BCSDIR /discover/nobackup/ltakacs/bcs/Ganymed-4_0/Ganymed-4_0_MERRA-2 => setenv BCSDIR /discover/nobackup/smahanam/bcs/Heracles-4_3/Heracles-4_3_MERRA-3/ Comment out L1185 "qsub $FILE" in /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/src/GMAO_Shared/GEOS_Util/post/gcmpost.script to avoid submitting the post plotting job. sbatch gcm_run.j ========== 20170105: 1. The experiment stopped at "AGCM Date: 2001/01/01 Time: 20:45:00". Sarith told me to comment out CNLAI through ICESOI in HISTORY.rc and try again. /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests > mv clm4 clm4_old Repeat the steps in 1 on 20170104 to re-do the experiment. Again the experiment crashed at the same time step above. Sarith found the bug in GEOS_CatchCNGridComp.F90. Fixed it and re-compiled. L5236-5252 in GEOS_CatchCNGridComp.F90 is now: ! copy CN_restart vars to catchcn_internal_rst gkw: only do if stopping ! --------------------------------------------------------------------- if(NextTime == StopTime) then call CN_exit(ntiles,nveg,nzone,ityp,fveg,cncol,var_col,cnpft,var_pft) i = 1 do iv = 1,VAR_PFT do nv = 1,NUM_VEG do nz = 1, NUM_ZON do n = 1,ntiles ! to ensure unused array elements don't have crazy numbers in restart files. if(fveg (n,nv,nz) == 0.) cnpft (n,i) = 0. end do i = i + 1 end do end do end do endif cd /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/src source g5_modules ./parallel_build.csh Chose "realclean". Set up an experiment following the steps in 1 on 20170104. ========== 20170106: 1. The experiment finished successfully this time. Checked the output on GrADS. Looks good. However, Sarith said (see his email on 20170106): "Don't you think that the program should not have been inside that block at the end of the 2nd day, in the first place. I was expecting the code is supposed to visit that block (call cn_exit) at the end of the 32nd day (the length of job_sgmt in CAP.rc)". This makes sense too. More investigation is needed. Sarith told me to add a print statement, re-compiled and run. The statements printed out suggest that "StopTime" is not the EndTime of the job segment as we thought. Sarith said he will fix this next week. 2. Contacted Jinyun Tang about the bug in subroutine hybrid yesterday. He replied the same day with some suggestions: ~~~~~~~~~~~~~~~~~~~~~~ Fanwei I should say that I did not look into this issue, probably you can add a safety screening among these lines to ensure no negative ci is allowed. For instance, you can set the smallest number as ca*0.01, and maximum number as ca. Then directly call brent, which will fix the issue. Note, you should still check if f(ca*0.01)*f(ca)<0. Please let me know if it works. I will then forward your discovery to my colleagues in NCAR. ~~~~~~~~~~~~~~~~~~~~~~ Don't fully understand what he meant. Modified subroutine ci_func to make it stop if ci is <0.01*cair or >cair. Modified subroutine hybrid to make sure x falls in the bounds of 0.01*cair and cair. Compiled without BOPT=g. cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/exec /bin/cp -pr ../Linux clm4.5/. Re-do the global run starting from 20000101 and ran for one year. cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests mv m0001hdeg_clm4.5 m0001hdeg_clm4.5_old cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/exec/clm4.5/Linux/bin source /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/src/g5_modules ./ldsetup setup /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/ /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/run/D0.5_clm4.5.exe /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/run/D0.5.bat --runmodel --monthsperjob 12 --landmodel catchCN cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/input/restart/ ln -s /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts/0.5D/catchcn_internal_rst /bin/rm regridded_restarts cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/run nedit lenkf.0.j & (change job name to "clm4.5" change the line for restart_path to "-restart_path ../input/restart/ \") qsub lenkf.0.j It crashed in L1715 of compute_rc.F90 "dx = - f1 * (x1-x0)/(f1-f0)" saying "floating invalid". Added a print statement. Re-compiled without BOPT=g. cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/exec /bin/cp -pr ../Linux clm4.5/. Re-run. ssh -XY discover-sp3 interactive.py -A sp3 -n 96 -a g0620 -X --debug cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/run /bin/rm ../output/global/rc_out/* ./lenkf.0.k Crashed: date_time_new 20000101_010730z get_forcing(): assuming GEOS-5 forcing data set read SWGDN from ../input/met_forcing/MERRA2_land_forcing///MERRA2_200/diag/Y2000/M01/MERRA2_200.tavg1_2d_rad_Nx.20000101.nc4 read SWLAND from ../input/met_forcing/MERRA2_land_forcing///MERRA2_200/diag/Y2000/M01/MERRA2_200.tavg1_2d_lfo_Nx.20000101.nc4 tile id: 4542 x0: 0.3718660 x1: 0.3718660 f0: 6.993237 f1: 6.993237 hybrid: f1==f0 Why is x0 so high? Next: print more information and run only on this tile 4542. ========== 20170109: 1. /discover/nobackup/fzeng/bcs/Heracles-4_3/Heracles-4_3_MERRA-3/DE_00720x00360_PE_0720x0360/DE_00720x00360_PE_0720x0360.til: 100 280070 117.2500 -24.2500 595 132 1.000000000000 280070 0.1282E+03 633428100000 4542 so tile 4542: lon = 117.2500; lat = -24.2500. /discover/nobackup/smahanam/bcs/Heracles-4_3/Heracles-4_3_MERRA-3/DE_00720x00360_PE_0720x0360/clsm/CLM_veg_typs_fracs: 4542 280070 16 17 14 15 0.00 67.75 0.00 32.25 15 14 67.75 32.25 so tile 4542 has both c3 and c4 plants. Re-run on this tile only. ssh -XY discover-sp3 (interactive.py -A sp3 -n 12 -a g0620 -X --debug) cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/run /bin/rm ../output/global/rc_out/* ./lenkf.0.k (use only 1 node, so no need the "interactive.py -A sp3 -n 12 -a g0620 -X --debug" above) test #1: x0: 14.87464 f0: 21.49606 x1: 14.72589 f1: 21.34732 iteration starts iter: 1 x0: 14.87464 f0: 21.49606 x1: 14.72589 f1: 21.34732 dx: -21.34732 x: -6.621425 After do while, iter: 1 x0: 14.87464 f0: 21.49606 x1: 14.72589 f1: 21.34732 dx: -21.34732 x: 4.052234 After ci_func, iter: 1 x0: 14.72589 f0: 21.34732 x1: 4.052234 f1: 10.67366 dx: -21.34732 x: 4.052234 iter: 2 x0: 14.72589 f0: 21.34732 x1: 4.052234 f1: 10.67366 dx: -10.67365 x: -6.621421 After do while, iter: 2 x0: 14.72589 f0: 21.34732 x1: 4.052234 f1: 10.67366 dx: -10.67365 x: 0.4943488 After ci_func, iter: 2 x0: 4.052234 f0: 10.67366 x1: 0.4943488 f1: 7.115738 dx: -10.67365 x: 0.4943488 iter: 3 x0: 4.052234 f0: 10.67366 x1: 0.4943488 f1: 7.115738 dx: -7.115671 x: -6.621322 After do while, iter: 3 x0: 4.052234 f0: 10.67366 x1: 0.4943488 f1: 7.115738 dx: -7.115671 x: 0.3737442 After ci_func, iter: 3 x0: 0.4943488 f0: 7.115738 x1: 0.3737442 f1: 6.995117 dx: -7.115671 x: 0.3737442 iter: 4 x0: 0.4943488 f0: 7.115738 x1: 0.3737442 f1: 6.995117 dx: -6.994180 x: -6.620436 After do while, iter: 4 x0: 0.4943488 f0: 7.115738 x1: 0.3737442 f1: 6.995117 dx: -6.994180 x: 0.3718660 After ci_func, iter: 4 x0: 0.3737442 f0: 6.995117 x1: 0.3718660 f1: 6.993237 dx: -6.994180 x: 0.3718660 iter: 5 x0: 0.3737442 f0: 6.995117 x1: 0.3718660 f1: 6.993237 dx: -6.983928 x: -6.612062 After do while, iter: 5 x0: 0.3737442 f0: 6.995117 x1: 0.3718660 f1: 6.993237 dx: -6.983928 x: 0.3718660 After ci_func, iter: 5 x0: 0.3718660 f0: 6.993237 x1: 0.3718660 f1: 6.993237 dx: -6.983928 x: 0.3718660 forrtl: error (73): floating divide by zero looks like the ci value couldn't go below 0.3718660, making x0=x1=0.3718660 (and thus f0=f1) in the next iteration, causing floating divide by zero. So is cair 37.18660, because the lower bound of ci is 0.01*cair? test #2: (modified print statements for better understanding the issue) x0: 14.87464 f0: 21.49606 cair: 37.18660 x1: 14.72589 f1: 21.34732 iteration starts iter: 1 x0: 14.87464 f0: 21.49606 x1: 14.72589 f1: 21.34732 dx: -21.34732 x: -6.621425 After do while, n: 1 x: 4.052234 After ci_func, iter: 1 x0: 14.72589 f0: 21.34732 x1: 4.052234 f1: 10.67366 iter: 2 x0: 14.72589 f0: 21.34732 x1: 4.052234 f1: 10.67366 dx: -10.67365 x: -6.621421 After do while, n: 2 x: 0.4943488 After ci_func, iter: 2 x0: 4.052234 f0: 10.67366 x1: 0.4943488 f1: 7.115738 iter: 3 x0: 4.052234 f0: 10.67366 x1: 0.4943488 f1: 7.115738 dx: -7.115671 x: -6.621322 After do while, n: 58 x: 0.3737442 After ci_func, iter: 3 x0: 0.4943488 f0: 7.115738 x1: 0.3737442 f1: 6.995117 iter: 4 x0: 0.4943488 f0: 7.115738 x1: 0.3737442 f1: 6.995117 dx: -6.994180 x: -6.620436 After do while, n: 3723 x: 0.3718660 After ci_func, iter: 4 x0: 0.3737442 f0: 6.995117 x1: 0.3718660 f1: 6.993237 iter: 5 x0: 0.3737442 f0: 6.995117 x1: 0.3718660 f1: 6.993237 dx: -6.983928 x: -6.612062 After do while, n: 66954781 x: 0.3718660 After ci_func, iter: 5 x0: 0.3718660 f0: 6.993237 x1: 0.3718660 f1: 6.993237 forrtl: error (73): floating divide by zero As expected, cair=37.18660. The initial x0/cair = 14.87464/37.18660 = 0.4, so it's a c4 plant. In subroutine hybrid, there is this block of code: if(iter>itmax)then !in case of failing to converge within itmax iterations !stop at the minimum function !this happens because of some other issues besides the stomatal conductance calculation !and it happens usually in very dry places and more likely with c4 plants. call ci_func(tile_id, minx, f1, gb_mol, je, cair, oair, lmr_z, par_z, rh_can, gs_mol, & c3flag, vcmax_z, cp, kc, ko, qe, tpu_z, kp_z, theta_cj, forc_pbot, bbb, mbb, ac, aj, ap, ag, an) exit endif So for c4 plants, ci value has to be very close to 0 for the f value to be close 0 enough. However, here we make ci value always >= 0.01*cair, which is 0.3718660. When ci = 0.3718660, f value = 6.993237, still quite large. test #3: (modified code trying to fix the bug) x0: 14.87464 f0: 21.49606 cair: 37.18660 x1: 14.72589 f1: 21.34732 iteration starts iter: 1 x0: 14.87464 f0: 21.49606 x1: 14.72589 f1: 21.34732 dx: -21.34732 x: -6.621425 After do while, n: 1 x: 4.052234 After ci_func, x0: 14.72589 f0: 21.34732 x1: 4.052234 f1: 10.67366 iter: 2 x0: 14.72589 f0: 21.34732 x1: 4.052234 f1: 10.67366 dx: -10.67365 x: -6.621421 After do while, n: 2 x: 0.4943488 After ci_func, x0: 4.052234 f0: 10.67366 x1: 0.4943488 f1: 7.115738 iter: 3 x0: 4.052234 f0: 10.67366 x1: 0.4943488 f1: 7.115738 dx: -7.115671 x: -6.621322 After do while, n: 58 x: 0.3737442 After ci_func, x0: 0.4943488 f0: 7.115738 x1: 0.3737442 f1: 6.995117 iter: 4 x0: 0.4943488 f0: 7.115738 x1: 0.3737442 f1: 6.995117 dx: -6.994180 x: -6.620436 After do while, n: 3723 x: 0.3718660 After ci_func, x0: 0.3737442 f0: 6.995117 x1: 0.3718660 f1: 6.993237 iter: 5 x0: 0.3737442 f0: 6.995117 x1: 0.3718660 f1: 6.993237 dx: -6.983928 x: -6.612062 After do while, n: 66954781 x: 0.3718660 After ci_func, x0: 0.3718660 f0: 6.993237 x1: 0.3718660 f1: 6.993237 final x0: 0.3718660 minx: 0.3718660 f0: 6.993237 x0: 26.03062 f0: 1.573047 cair: 37.18660 x1: 25.77031 f1: 1.236067 iteration starts iter: 1 x0: 26.03062 f0: 1.573047 x1: 25.77031 f1: 1.236067 dx: -0.9548193 x: 24.81549 After do while, n: 0 x: 24.81549 After ci_func, x0: 25.77031 f0: 1.236067 x1: 24.81549 f1: -1.1626244E-02 x0: 14.87464 f0: 7.589281 cair: 37.18660 x1: 14.72589 f1: 7.440535 iteration starts iter: 1 x0: 14.87464 f0: 7.589281 x1: 14.72589 f1: 7.440535 dx: -7.440534 x: 7.285358 After do while, n: 0 x: 7.285358 After ci_func, x0: 14.72589 f0: 7.440535 x1: 7.285358 f1: 1.9073486E-06 x0: 26.03062 f0: -10.71537 cair: 37.18660 x1: 25.77031 f1: -10.98143 iteration starts iter: 1 x0: 26.03062 f0: -10.71537 x1: 25.77031 f1: -10.98143 dx: 10.74389 x: 36.51420 After do while, n: 0 x: 36.51420 After ci_func, x0: 25.77031 f0: -10.98143 x1: 36.51420 f1: -1.3003469E-02 iter: 2 x0: 25.77031 f0: -10.98143 x1: 36.51420 f1: -1.3003469E-02 dx: 1.2737262E-02 x: 36.52694 x0: 14.87464 f0: 21.49403 cair: 37.18660 x1: 14.72589 f1: 21.34529 iteration starts iter: 1 x0: 14.87464 f0: 21.49403 x1: 14.72589 f1: 21.34529 dx: -21.34529 x: -6.619395 After do while, n: 1 x: 4.053248 After ci_func, x0: 14.72589 f0: 21.34529 x1: 4.053248 f1: 10.67264 iter: 2 x0: 14.72589 f0: 21.34529 x1: 4.053248 f1: 10.67264 dx: -10.67264 x: -6.619391 After do while, n: 2 x: 0.4957018 After ci_func, x0: 4.053248 f0: 10.67264 x1: 0.4957018 f1: 7.115059 iter: 3 x0: 4.053248 f0: 10.67264 x1: 0.4957018 f1: 7.115059 dx: -7.114986 x: -6.619285 After do while, n: 57 x: 0.3730296 After ci_func, x0: 0.4957018 f0: 7.115059 x1: 0.3730296 f1: 6.992367 iter: 4 x0: 0.4957018 f0: 7.115059 x1: 0.3730296 f1: 6.992367 dx: -6.991230 x: -6.618201 After do while, n: 6007 x: 0.3718660 After ci_func, x0: 0.3730296 f0: 6.992367 x1: 0.3718660 f1: 6.991207 iter: 5 x0: 0.3730296 f0: 6.992367 x1: 0.3718660 f1: 6.991207 dx: -7.015283 x: -6.643417 After do while, n: 470787695 x: 0.3718660 After ci_func, x0: 0.3718660 f0: 6.991207 x1: 0.3718660 f1: 6.991207 final x0: 0.3718660 minx: 0.3718660 f0: 6.991207 Looks correct now. 2. Did a global test run to run. cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/run /bin/rm ../output/global/rc_out/* qsub lenkf.0.j final x0: 0.3954279 minx: 0.3954386 f0: 54.29114 final x0: 0.3848898 minx: 0.3848898 f0: 26.74969 final x0: 0.3877727 minx: 0.3877727 f0: 2.125164 final x0: 0.3778148 minx: 0.3778148 f0: 2.413769 final x0: 0.3811840 minx: 0.3811840 f0: 2.154686 final x0: 0.3742619 minx: 0.3742619 f0: 2.932575 final x0: 0.3815705 minx: 0.3815705 f0: 19.61778 final x0: 0.3874066 minx: 0.3874066 f0: 0.1562386 final x0: 0.3735659 minx: 0.3735659 f0: 18.91019 final x0: 0.3937417 minx: 0.3937417 f0: 44.02775 final x0: 0.3847470 minx: 0.3847470 f0: 23.08481 final x0: 0.3804815 minx: 0.3804815 f0: 16.60377 final x0: 0.3926388 minx: 0.3926388 f0: 24.43692 final x0: 0.3898864 minx: 0.3898864 f0: 0.5546684 final x0: 0.3859620 minx: 0.3859620 f0: 9.898438 final x0: 0.3815705 minx: 0.3815705 f0: 10.71646 final x0: 0.3859719 minx: 0.3859720 f0: 7.902290 final x0: 0.3868071 minx: 0.3868071 f0: 0.2300949 final x0: 0.3915845 minx: 0.3915845 f0: 1.193447 final x0: 0.3874078 minx: 0.3874078 f0: 3.724277 final x0: 0.3862179 minx: 0.3862179 f0: 3.905331 final x0: 0.3877814 minx: 0.3877814 f0: 18.58075 final x0: 0.3826896 minx: 0.3826896 f0: 18.09647 final x0: 0.3737405 minx: 0.3737405 f0: 33.73558 final x0: 0.3899175 minx: 0.3899175 f0: 6.809063 final x0: 0.3888870 minx: 0.3888870 f0: 0.7636833 final x0: 0.3892798 minx: 0.3892798 f0: 12.03304 final x0: 0.3903190 minx: 0.3903190 f0: 17.67704 final x0: 0.3897085 minx: 0.3897085 f0: 7.769131 final x0: 0.3852355 minx: 0.3852347 f0: 1.072891 final x0: 0.3922512 minx: 0.3922512 f0: 25.93940 final x0: 0.3842580 minx: 0.3842580 f0: 42.07596 final x0: 0.3907049 minx: 0.3907049 f0: 8.940388 final x0: 0.3866659 minx: 0.3866659 f0: 0.1971626 final x0: 0.3912155 minx: 0.3912155 f0: 11.84724 final x0: 0.3902541 minx: 0.3902541 f0: 9.567188 final x0: 0.3872351 minx: 0.3872351 f0: 4.769440 final x0: 0.3885926 minx: 0.3885926 f0: 8.081329 final x0: 0.3892255 minx: 0.3892255 f0: 3.405193 final x0: 0.3886082 minx: 0.3886082 f0: 3.289249 final x0: 0.3886082 minx: 0.3886082 f0: 6.220699 final x0: 0.3873333 minx: 0.3873333 f0: 2.309937 final x0: 0.3847470 minx: 0.3847470 f0: 32.95952 final x0: 0.3737405 minx: 0.3737405 f0: 33.13913 When ci=0.01*cair (the lower bound), f0 could still be as large as 54.29114 or even larger occasionally (in these cases ci and thus photosynthesise would be overestimated too much). fval = ci - cair + an(p,iv) * forc_pbot(g) * (1.4_r8*gs_mol+1.6_r8*gb_mol) / (gb_mol*gs_mol) 3. Tried using 0 as the lower bound for ci and did a global test run. cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/run /bin/rm ../output/global/rc_out/* /bin/rm ../output/global/rc_out/Y2000/M01/* qsub lenkf.0.j (1) In this case, brent is called because f0*f1<0. x0: 14.87464 f0: 21.49606 cair: 37.18660 x1: 14.72589 f1: 21.34732 iteration starts iter: 1 x0: 14.87464 f0: 21.49606 x1: 14.72589 f1: 21.34732 dx: -21.34732 x: -6.621425 After do while, n: 1 x: 4.052234 After ci_func, x0: 14.72589 f0: 21.34732 x1: 4.052234 f1: 10.67366 iter: 2 x0: 14.72589 f0: 21.34732 x1: 4.052234 f1: 10.67366 dx: -10.67365 x: -6.621421 After do while, n: 2 x: 0.4943488 After ci_func, x0: 4.052234 f0: 10.67366 x1: 0.4943488 f1: 7.115738 iter: 3 x0: 4.052234 f0: 10.67366 x1: 0.4943488 f1: 7.115738 dx: -7.115671 x: -6.621322 After do while, n: 14 x: 1.9970715E-02 After ci_func, x0: 0.4943488 f0: 7.115738 x1: 1.9970715E-02 f1: 6.640331 iter: 4 x0: 0.4943488 f0: 7.115738 x1: 1.9970715E-02 f1: 6.640331 dx: -6.625964 x: -6.605993 After do while, n: 331 x: 1.2993813E-05 After ci_func, x0: 1.9970715E-02 f0: 6.640331 x1: 1.2993813E-05 f1: -10.24713 (2) In this case, the ci value is increased because the initial value guessed is too low. x0: 26.03062 f0: -10.71537 cair: 37.18660 x1: 25.77031 f1: -10.98143 iteration starts iter: 1 x0: 26.03062 f0: -10.71537 x1: 25.77031 f1: -10.98143 dx: 10.74389 x: 36.51420 After do while, n: 0 x: 36.51420 After ci_func, x0: 25.77031 f0: -10.98143 x1: 36.51420 f1: -1.3003469E-02 iter: 2 x0: 25.77031 f0: -10.98143 x1: 36.51420 f1: -1.3003469E-02 dx: 1.2737262E-02 x: 36.52694 (3) In this case, the f value reaches 0 at a very small ci value (1.6989361E-08). x0: 14.87464 f0: 21.49403 cair: 37.18660 x1: 14.72589 f1: 21.34529 iteration starts iter: 1 x0: 14.87464 f0: 21.49403 x1: 14.72589 f1: 21.34529 dx: -21.34529 x: -6.619395 After do while, n: 1 x: 4.053248 After ci_func, x0: 14.72589 f0: 21.34529 x1: 4.053248 f1: 10.67264 iter: 2 x0: 14.72589 f0: 21.34529 x1: 4.053248 f1: 10.67264 dx: -10.67264 x: -6.619391 After do while, n: 2 x: 0.4957018 After ci_func, x0: 4.053248 f0: 10.67264 x1: 0.4957018 f1: 7.115059 iter: 3 x0: 4.053248 f0: 10.67264 x1: 0.4957018 f1: 7.115059 dx: -7.114986 x: -6.619285 After do while, n: 14 x: 2.1369368E-02 After ci_func, x0: 0.4957018 f0: 7.115059 x1: 2.1369368E-02 f1: 6.639767 iter: 4 x0: 0.4957018 f0: 7.115059 x1: 2.1369368E-02 f1: 6.639767 dx: -6.626359 x: -6.604989 After do while, n: 310 x: 6.2748790E-05 After ci_func, x0: 2.1369368E-02 f0: 6.639767 x1: 6.2748790E-05 f1: 5.191959 iter: 5 x0: 2.1369368E-02 f0: 6.639767 x1: 6.2748790E-05 f1: 5.191959 dx: -7.6407336E-02 x: -7.6344587E-02 After do while, n: 1217 x: 1.6989361E-08 After ci_func, x0: 6.2748790E-05 f0: 5.191959 x1: 1.6989361E-08 f1: 0.0000000E+00 4. Added print statements so that we know where the iteration ends and the final ci and f values. Next, compile and test run tomorrow. ========== 20170110: 1. Following 4 on 20170109 above, compiled and did a global test run. cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/run /bin/rm ../output/global/rc_out/* /bin/rm ../output/global/rc_out/Y2000/M01/* qsub lenkf.0.j (1) iteration ends after call brent; final f value is small enough. Good! x0: 26.03062 f0: 1.573047 cair: 37.18660 c3flag: T x1: 25.77031 f1: 1.236067 iteration starts iter: 1 x0: 26.03062 f0: 1.573047 x1: 25.77031 f1: 1.236067 dx: -0.9548193 x: 24.81549 After do while, n: 0 x: 24.81549 After ci_func, x0: 25.77031 f0: 1.236067 x1: 24.81549 f1: -1.1626244E-02 Before call brent, x0: 25.77031 f0: 1.236067 x1: 24.81549 f1: -1.1626244E-02 iteration ends after call brent, final x0: 24.81549 final fval: -1.1626244E-02 (2) iteration ends: abs(f1)<=eps1; final f value is small enough. Good! x0: 14.87464 f0: 7.589281 cair: 37.18660 c3flag: F x1: 14.72589 f1: 7.440535 iteration starts iter: 1 x0: 14.87464 f0: 7.589281 x1: 14.72589 f1: 7.440535 dx: -7.440534 x: 7.285358 After do while, n: 0 x: 7.285358 After ci_func, x0: 14.72589 f0: 7.440535 x1: 7.285358 f1: 1.9073486E-06 iteration ends: abs(f1)<=eps1, final x0: 7.285358 final fval: 1.9073486E-06 (3) iteration ends: abs(dx)cair(p) tile_id: 20094 ci out of bound: -6.4979553E-02 cair: 36.49923 an: 5.3104747E-04 forc_pbot: 93587.76 gs_mol: 2.174786 gb_mol: 3691266. bbb: 1.000000 tile_id: 20517 ci out of bound: -5.7456970E-02 cair: 36.05539 an: 4.9542205E-04 forc_pbot: 92449.72 gs_mol: 2.029268 gb_mol: 3676656. bbb: 1.000000 tile_id: 20063 ci out of bound: -5.7666779E-02 cair: 35.92472 an: 4.6601827E-04 forc_pbot: 92114.67 gs_mol: 1.908806 gb_mol: 3705335. bbb: 1.000000 tile_id: 20517 ci out of bound: -5.5194855E-02 cair: 36.05539 an: 4.9542205E-04 forc_pbot: 92449.72 gs_mol: 2.029395 gb_mol: 3676656. bbb: 1.000000 tile_id: 20063 ci out of bound: -5.5034637E-02 cair: 35.92472 an: 4.6601827E-04 forc_pbot: 92114.67 gs_mol: 1.908946 gb_mol: 3705335. bbb: 1.000000 tile_id: 20041 ci out of bound: -5.7182312E-02 cair: 37.51011 an: 3.9273634E-04 forc_pbot: 96179.77 gs_mol: 1.608774 gb_mol: 3903380. bbb: 1.000000 tile_id: 20041 ci out of bound: -5.3703308E-02 cair: 37.51011 an: 3.9273634E-04 forc_pbot: 96179.77 gs_mol: 1.608923 gb_mol: 3705335. bbb: 1.000000 tile_id: 20261 ci out of bound: -6.9377899E-02 cair: 37.23918 an: 3.0432193E-04 forc_pbot: 95485.07 gs_mol: 1.246179 gb_mol: 3853556. bbb: 1.000000 tile_id: 20261 ci out of bound: -6.6482544E-02 cair: 37.23918 an: 3.0432193E-04 forc_pbot: 95485.07 gs_mol: 1.246276 gb_mol: 3853556. bbb: 1.000000 tile_id: 19184 ci out of bound: -9.7640991E-02 cair: 38.19430 an: 3.3659654E-04 forc_pbot: 97934.09 gs_mol: 1.377388 gb_mol: 3868556. bbb: 1.000000 tile_id: 20244 ci out of bound: -7.1460724E-02 cair: 38.08064 an: 2.9918458E-04 forc_pbot: 97642.66 gs_mol: 1.225125 gb_mol: 3891990. bbb: 1.000000 tile_id: 19180 ci out of bound: -9.8712921E-02 cair: 37.94380 an: 3.4795460E-04 forc_pbot: 97291.79 gs_mol: 1.423802 gb_mol: 3843935. bbb: 1.000000 tile_id: 19385 ci out of bound: -5.8773041E-02 cair: 37.96838 an: 3.3059213E-04 forc_pbot: 97354.83 gs_mol: 1.354180 gb_mol: 3841635. bbb: 1.000000 tile_id: 19180 ci out of bound: -9.6485138E-02 cair: 37.94380 an: 3.4795460E-04 forc_pbot: 97291.79 gs_mol: 1.423886 gb_mol: 3843935. bbb: 1.000000 tile_id: 19385 ci out of bound: -5.6678772E-02 cair: 37.96838 an: 3.3059213E-04 forc_pbot: 97354.83 gs_mol: 1.354254 gb_mol: 3841635. bbb: 1.000000 tile_id: 19175 ci out of bound: -0.1103439 cair: 37.99253 an: 3.7936724E-04 forc_pbot: 97416.74 gs_mol: 1.551872 gb_mol: 3835980. bbb: 1.000000 tile_id: 19176 ci out of bound: -8.9740753E-02 cair: 37.90034 an: 3.7915097E-04 forc_pbot: 97180.37 gs_mol: 1.551817 gb_mol: 3834127. bbb: 1.000000 tile_id: 19383 ci out of bound: -0.1559334 cair: 37.98077 an: 3.3688440E-04 forc_pbot: 97386.59 gs_mol: 1.376439 gb_mol: 3838588. bbb: 1.000000 tile_id: 19175 ci out of bound: -0.1081848 cair: 37.99253 an: 3.7936724E-04 forc_pbot: 97416.74 gs_mol: 1.551960 gb_mol: 3835980. bbb: 1.000000 tile_id: 19176 ci out of bound: -8.7581635E-02 cair: 37.90034 an: 3.7915097E-04 forc_pbot: 97180.37 gs_mol: 1.551906 gb_mol: 3834127. bbb: 1.000000 tile_id: 19383 ci out of bound: -0.1538315 cair: 37.98077 an: 3.3688440E-04 forc_pbot: 97386.59 gs_mol: 1.376515 gb_mol: 3838588. bbb: 1.000000 tile_id: 19383 ci out of bound: -0.1828613 cair: 38.08456 an: 3.1981565E-04 forc_pbot: 97652.73 gs_mol: 1.305795 gb_mol: 3884449. bbb: 1.000000 tile_id: 4258 ci out of bound: -42.77575 cair: 23.93880 an: 1.1728124E-03 forc_pbot: 61381.54 gs_mol: 1.726498 gb_mol: 2592558. bbb: 1.000000 tile_id: 4258 ci out of bound: -42.76075 cair: 23.93880 an: 1.1728124E-03 forc_pbot: 61381.54 gs_mol: 1.726886 gb_mol: 2592558. bbb: 1.000000 tile_id: 4258 ci out of bound: -42.93273 cair: 23.93880 an: 1.1760489E-03 forc_pbot: 61381.54 gs_mol: 1.727198 gb_mol: 2592558. bbb: 1.000000 Crashed at L1473 in CNPhenologyMod.F90 "leafc_to_litter(p) = prev_leafc_to_litter(p) + t1*(leafc(p) - prev_leafc_to_litter(p)*offset_counter(p))". Error message is "floating overflow". (6) Jinyun just got back to me: "Fanwei Your were right. I too found the negative ci. However, this may come as the cost of the limited precision of the computer can do. The ci as returned from the brent solver is positive, when it is used to update an again, it will produce an that is not compatible with the stomata conductance that was calculated from the last iteration. One possible to further reduce such occurrence is to decrease the error tolerance in exiting from brent. However, the cost would be significant increase in computational cost. I also added one more update of an, after the call of brent, however, this resulted in breaking down of the ball-berry equation, which is again due to the limited precision of the computer can achieve. On way to overcome this is actually use the ci returned from hybrid when (gs_mol /= bbb), then you would get rid of the negative ci problem." In subroutine Photosynthesis, made this change: ! ci_z(p,iv) = cair(p) - an(p,iv) * forc_pbot(g) * (1.4_r8*gs_mol(p,iv)+1.6_r8*gb_mol(p)) / (gb_mol(p)*gs_mol(p,iv)) ! suggested by Jinyun Tang, Jan 2017: only update ci_z when gs_mol(p,iv) is very close to bbb if(abs(gs_mol(p,iv)-bbb(p))<1.e-14_r8) then ci_z(p,iv) = cair(p) - an(p,iv) * forc_pbot(g) * (1.4_r8*gs_mol(p,iv)+1.6_r8*gb_mol(p)) / (gb_mol(p)*gs_mol(p,iv)) end if And did a global test run and ran for one month. Now I don't see negative ci values any more. However, it still crashed at L1473 in CNPhenologyMod.F90 "leafc_to_litter(p) = prev_leafc_to_litter(p) + t1*(leafc(p) - prev_leafc_to_litter(p)*offset_counter(p))". Error message is "floating overflow". (7) Jinyun pointed out: "Also remember to update with the following if(abs(gs_mol(p,iv)-bbb(p))<1.e-14_r8) then ci_z(p,iv) = cair(p) - an(p,iv) * forc_pbot(g) * (1.4_r8*gs_mol(p,iv)+1.6_r8*gb_mol(p)) / (gb_mol(p)*gs_mol(p,iv)) else ci_z(p,iv)=ciold (!if I remember correctly). end if" Made this change and did a global test run and ran for one month. It still crashed at date_time_new 20000111_150730z at L1473 in CNPhenologyMod.F90 "leafc_to_litter(p) = prev_leafc_to_litter(p) + t1*(leafc(p) - prev_leafc_to_litter(p)*offset_counter(p))". Error message is "floating overflow". Before these above runs that crashed, I compiled with BOPT=g. Try compile without BOPT=g tomorrow. ========== 20170126: 1. Continued to investigate why the test runs crashed in CNPhenologyMod.F90. Compiled without BOPT=g. The global test run still crashed at date_time_new 20000111_150730z saying floating overflow, this time at L1474 in CNPhenologyMod.F90 "frootc_to_litter(p) = prev_frootc_to_litter(p) + t1*(frootc(p) - prev_frootc_to_litter(p)*offset_counter(p))" though. Tried running on these regions: (1) 15N-35N, 80E-100E: finished successfully. (2) 90S-90N, 80E-180E: finished successfully. (3) 90S-90N, 0-80E: finished successfully. (4) 90S-90N, 180W-90W: finished successfully. (5) 90S-90N, 90W-0: MINLON = -90.00000 , MAXLON = 0.0000000E+00, MINLAT = -90.00000 , MAXLAT = 90.00000 , Crashed. (6) 0-90N, 90W-0: MINLON = -90.00000 , MAXLON = 0.0000000E+00, MINLAT = 0.0000000E+00, MAXLAT = 90.00000 Crashed. (7) MINLON = -90.00000 , MAXLON = -45.00000 , MINLAT = 0.0000000E+00, MAXLAT = 90.00000 Crashed. (8) MINLON = -90.00000 , MAXLON = -45.00000 , MINLAT = 45.00000 , MAXLAT = 90.00000 , Finished successfully. (9) MINLON = -90.00000 , MAXLON = -60.00000 , MINLAT = 0.0000000E+00, MAXLAT = 45.00000 , Crashed. (10) MINLON = -90.00000 , MAXLON = -60.00000 , MINLAT = 20.00000 , MAXLAT = 45.00000 , Finished successfully. So one of the problematic regions is 0-20N, 90W-60W. More regional runs to further narrow down. (11) MINLON = -90.00000 , MAXLON = -75.00000 , MINLAT = 0.0000000E+00, MAXLAT = 10.00000 , Crashed. 2. Sarith suggested to change if(abs(gs_mol(p,iv)-bbb(p))<1.e-14_r8) then an(p,iv)=0._r8 ag(p,iv)=lmr_z(p,iv) end if to if(abs(gs_mol(p,iv)-bbb(p))<1.e-14_r8) then an(p,iv)=1.e-20 ! 0._r8 ag(p,iv)=max(lmr_z(p,iv),1.e-20) ! lmr_z(p,iv) end if and see if it would help. Still crashed. 3. Tried passing tile_id to subroutine CNOffsetLitterfall in CNPhenolobyMod but don't know how to do that, because the CLM structure is not on tile basis. Sarith modified some routines below to help with this. GEOSland_GridComp/Shared/update_model_paras.F90 CNPhenologyMod.F90 process_cn.F90 ========== 20170130: 1. Modified these subroutines to identify the problematic tile: CN_DriverMod.F90 CNEcosystemDynMod.F90 CNPhenologyMod.F90 Compiled with BOPT=g. Debug using Allinea: 90W-75W, 0-10N cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/run /bin/rm ../output/global/rc_out/* ddt /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/exec/clm4.5/Linux/bin/LDASsaCN_mpi.x & copy the block from work_path to the end in the lenkf.0.j file to auguments in the ddt window. click "run" click the play button on the upper left to make it run until it stops (1) Ran on 90W-75W, 0-10N: Identified one of the problematic tiles: tile 13966 /discover/nobackup/smahanam/bcs/Heracles-4_3/Heracles-4_3_MERRA-3/DE_00720x00360_PE_0720x0360/DE_00720x00360_PE_0720x0360.til: 100 245556 -75.2500 0.7500 210 182 1.000000000000 245556 0.1730E+03 549497000000 13966 /discover/nobackup/smahanam/bcs/Heracles-4_3/Heracles-4_3_MERRA-3/DE_00720x00360_PE_0720x0360/clsm/CLM_veg_typs_fracs: 13966 245556 4 4 18 19 94.02 0.00 0.00 5.98 5 16 94.02 5.98 (2) Ran on 75.5W-75W, 0.5N-1N (single tile 13966) It crashed in this line in subroutine CNOffsetLitterfall in CNPhenologyMod.F90: leafc_to_litter(p) = prev_leafc_to_litter(p) + t1*(leafc(p) - prev_leafc_to_litter(p)*offset_counter(p)) leafc is a small negative number ~-0.035; offset_counter = 5400; t1 = dt * 2.0 / (offset_counter(p) * offset_counter(p)) = 10800 * 2.0 / (5400*5400) so leafc_to_litter(p) = t1*leafc(p) + (1-t1*offset_counter(p))*prev_leafc_to_litter(p) = t1*leafc(p) + (1-4)*prev_leafc_to_litter(p) = t1*leafc(p) - 3*prev_leafc_to_litter(p) ~= - 3*prev_leafc_to_litter(p) and prev_leafc_to_litter(p) = leafc_to_litter(p) in subroutine CNOffsetLitterfall. Therefore, prev_leafc_to_litter keeps increasing with flipping signs in between negative and positive. At some point it reaches 1.6e+35, causing floating overflow! Read the CN routines to understand how leafc is updated. ========== 20170131: 1. Tried using the compute_rc.F90 with my fix earlier. Surprisingly, it also has the problem now. It was able to simulate the entire year of 2000 earlier and the output looked correct. Why it crashed now? One difference is, I added t10 here, while the previous run with my fix that finished year 2000 successfully used tm as t10 in compute_rc.F90. Would the crash here be due to the t10 I added? Restored catch_types.F90, clsm_ensdrv_init_routines.F90, process_cn.F90, and compute_rc.F90 with my fix to before t10 and used the restart file without t10. It finished successfully on this problematic tile 13966. Now, only change compute_rc.F90 to Jinyun's fix without t10: successful! Now include t10 in catch_types.F90, clsm_ensdrv_init_routines.F90, process_cn.F90 (init_t10 = .false.), and compute_rc.F90, use Jinyun's fix in compute_rc.F90, use /discover/nobackup/rreichle/l_data/LandRestarts_for_Regridding/CatchCN/catchcn_internal_clm45 as the restart file: crashed in L1156 process_cn.F90 "cn_progn(:)%T10 = ((n10-1)*cn_progn(:)%T10 + met_force(:)%Tair) / n10", floating overflow. cn_progn(:)%T10 is 9.97e+36!! Why? Did "ncdump -v T10 /discover/nobackup/rreichle/l_data/LandRestarts_for_Regridding/CatchCN/catchcn_internal_clm45" and found that there is no value for T10 in this file. The restart file I created after adding T10 /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts/0.5D_t10/catchcn_internal_rst has valid T10 values. Now change the restart file to /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts/0.5D_t10/catchcn_internal_rst. As seen earlier, it crashed in L1485 CNPhenologyMod.F90 "leafc_to_litter(p) = prev_leafc_to_litter(p) + t1*(leafc(p) - prev_leafc_to_litter(p)*offset_counter(p))". Allinea shows that, in the first time step, cn_progn(:)%T10 is 261.1K while met_force(:)%Tair is 299.7K. cn_progn(:)%T10 seems too low for this tile located in the tropics (75.25W, 0.75N). This could be why leafc is so low. met_force(:)%Tair is between 296K and 300K on 20000101 for this tile. I wonder why cn_progn(:)%T10 in the restart file I created would be so low (261.1K). Check if I did something wrong in creating the restart file tomorrow. ========== 20170201: 1. Checked process_cn.F90 to see what I did wrong in creating the restart file with t10. Couldn't find the problem by just reading the code. Run on Allinea to find out. First, set init_t10 = .true. in process_cn.F90 and compiled with BOPT=g. Then use the restart file without T10 to run on Allinea. Ran for one day 20000101 (as I did earlier to add t10 to the restart file). cn_progn%T10 values for this tile look correct, never falling below 297K or going above 300K, consistent with the met_force%Tair of this tile. I wonder why I got wrong cn_progn%T10 value for this tile when I did the global run earlier. Will do the same run globally and see if I will get wrong t10 number for this tile again. Ran globally and print out t10 for tile 13966 every 100 model time steps: if(mod(istep,100)==1 .and. tile_coord(n)%tile_id==13966) print *,'t10:',cn_progn(n)%T10 t10: 299.6906 t10: 297.4673 The cn_progn%T10 for this tile looks correct here. TILE_ID in the 0.5D restart file without t10 Sarith created IS NOT re-ordered. TILE_ID = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, ... Since I simply copied the restart file at the end of the global run ending 20000102_00z (see my notes on 20170120), TILE_ID in the 0.5D restart file with t10 I created IS re-ordered. TILE_ID = 35089, 35088, 34729, 34730, 34728, 35092, 35093, 35094, 35456, 35809, 35457, 35808, 35458, 35459, 36881, 37264, 37265, 37651, 37266, 37267, 36886, 37268, 37269, 36887, 37270, 36888, 36889, 36510, 36157, 36508, 36883, 36882, 36884, 36885, 36505, 36504, 36509, 36507, 36506, 36161, 35812, 36160, 35811, 36156, 36155, 35805, 35454, 35806, 35807, 35455, 35090, 35091, 36513, 36159, 35810, 36512, 36511, 36158, 35452, 35802, 36154, 36153, 35451, 35804, 36152, 35803, 35087, 35453, 43159, Therefore, wrong restart information is assigned to each tile when I use this restart file. 2. Modified Sarith's program /discover/nobackup/smahanam/reorder_LDASsa_rst.F90 to put the tile in order in /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts/0.5D_t10/catchcn_internal_rst. /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts/0.5D_t10 > mv catchcn_internal_rst catchcn_internal_rst_re-ordered cd /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/src/GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts cp /discover/nobackup/smahanam/reorder_LDASsa_rst.F90 . nedit reorder_LDASsa_rst.F90 & setenv ESMADIR /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/ source $ESMADIR/src/g5_modules gmake install reorder_LDASsa_rst (to run it) cp -p reorder_LDASsa_rst.F90 reorder_LDASsa_rst /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts/0.5D_t10/. Checked the TILE_ID in the newly created restart file /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts/0.5D_t10/catchcn_internal_rst (ncdump -v TILE_ID catchcn_internal_rst). They are in ascending order now. 3. Checked the t10 values in the newly created restart file /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts/0.5D_t10/catchcn_internal_rst using ~/Catchment/matlab/ read_catchcn_internal_rst.m. The t10 value for tile 13966 is 297.8205. This looks correct. 4. Restored init_t10 to .false. in process_cn.F90 and compiled with BOPT=g. Run on Allinea on tile 13966. Now cn_progn%T10 looks correct for the first few time steps checked. The run finished one month (200001) of simulation successfully. Wanted to check leafc in subroutine CNOffsetLitterfall in CNPhenolobyMod.F90. However, the run didn't stop at the breakpoint I set and just went through. Tried again and it didn't stop either, so was not able to check leafc values. 5. Do a global run and run one year (2000) to take a look at the output gpp and lai. ========== 20170202: 1. The run ended at date_time_new 20001127_210730z due to time limit (couldn't finish 12 months of simulation within 8 hours). Found that the output is daily - I forgot to change the namelist (from the test runs) to make it output monthly. That's why it's taking so long. Changed the namelist; removed the daily output; and remove the tile info in rc_out (/bin/rm ../output/global/rc_out/*). Submitted the job again. Crashed: date_time_new 20000113_060730z get_forcing(): assuming GEOS-5 forcing data set get_GEOS_gfio_openfile(): Unsuccessfully tried to open files: ../input/met_forcing/MERRA2_land_forcing///MERRA2_200.tavg1_2d_rad_Nx.20000113.nc4 ../input/met_forcing/MERRA2_land_forcing///MERRA2_200/diag/Y2000/M01/MERRA2_200.tavg1_2d_rad_Nx.20000113.nc4 LDAS ERROR (3000) from get_GEOS: error opening file forrtl: error (78): process killed (SIGTERM) LDASsaCN_mpi.x 0000000001928A5A cnsummarymod_mp_c 555 CNSummaryMod.F90 LDASsaCN_mpi.x 0000000001491CD9 cnecosystemdynmod 163 CNEcosystemDynMod.F90 LDASsaCN_mpi.x 0000000001045835 cn_drivermod_mp_c 555 CN_DriverMod.F90 LDASsaCN_mpi.x 0000000000B112F8 process_cn_mp_pro 1406 process_cn.F90 Why? Noticed that my SMAP M09 run is also terminated. date_time_new 20041219_210730z get_forcing(): assuming GEOS-5 forcing data set get_GEOS_gfio_openfile(): Unsuccessfully tried to open files: ../input/met_forcing/MERRA2_land_forcing///MERRA2_300.tavg1_2d_rad_Nx.20041219.nc4 ../input/met_forcing/MERRA2_land_forcing///MERRA2_300/diag/Y2004/M12/MERRA2_300.tavg1_2d_rad_Nx.20041219.nc4 LDAS ERROR (3000) from get_GEOS: error opening file So it seems that the MERRA-2 data is not there anymore? Talked to Qing. She checked the files but they are there and are "rw-r-r". Don't know what happened. Simply re-submitted the CLM4.5 test run again and see how it goes. It passed 20000113_060730z and is still running. 2. Test the modifications on GEOS-5. Don't include T10 yet because it requires T10 in the restart file. So far I have modified these offline routines. compute_rc.F90 process_cat.F90 in Greg's point driver, equivalent to process_cn.F90 in LDAS, and GEOS_CatchCNGridComp.F90 in GEOS5. unified_rc3f_matrix_calc.F90 in Greg's point driver, equivalent to catchmentCN.F90 in LDAS and GEOS5. pftvarcon.F90 clm_varpar.F90 CN_DriverMod.F90 in offline code (not this point driver): get_CN_LAI. catch_types.F90 clsm_ensdrv_init_routines.F90 GEOSland_GridComp/Shared/update_model_paras.F90 No change has been made to catchmentCN.F90 so far. Copied compute_rc.F90, pftvarcon.F90, clm_varpar.F90, CN_DriverMod.F90, and update_model_paras.F90 from offline sand box to GEOS-5 sand box. Modified compute_rc to use tm as t10 for now. Commented out all "tile_id" in compute_rc and the subroutines therein. L4813 in GEOS_CatchCNGridComp.F90 " cat_id = nint(tile_id) ! gkw: temporary for debugging" has cat_id. Maybe I don't need to remove tile_id from compute_rc. catch_types.F90 and clsm_ensdrv_init_routines.F90 were modified for adding T10, so don't change them yet. Also, it seems that they don't exist in GEOS-5 code. What's their equivalent routines? Will explore this later. Included the changes (not including those for adding t10) in process_cn.F90 to GEOS_CatchCNGridComp.F90. These changes in process_cn.F90 have not been included in GEOS_CatchCNGridComp.F90 yet: (line numbers in process_cn.F90) Done! (1) L351: integer ityp_cat1(N_cat), ityp_cat2(N_cat), ityp_tmp(N_cat) Done! (2) deallocate: tlaz, albdir, albdif Done! (3) This block below needs modifications. ! compute snow-free albedo for each PFT in each zone gkw: assume the snow albedo is not very important ! -------------------------------------------------- do nv = 1,num_veg ityp_tmp(:) = map_cat(ityz(:,nv)) ! fzeng: note that this is not exactly the same as calling sibalb_vis in the unified model because the ! "if(fveg(i)>1.e-4 .and. zth(i)>0.01)" branch in subroutine sibalb_vis is absent in the current subroutine sibalb. ! ----------------------------------------------------------------------------------------------------------------- call SIBALB(N_cat, ityp_tmp, elaz(:,nv), veg_param%grn, sunang, & cn_param%BGALBVR, cn_param%BGALBVF, cn_param%BGALBNR, cn_param%BGALBNF, & ! MODIS soil background albedo ALBVR, ALBNR, ALBVF, ALBNF, MODIS_SCALE=.TRUE. ) ! inst snow-free albedos ! fzeng: is fsnow the same as asnow? ! ---------------------------------- where(tlaz(:,nv) > 0.) fsnow(:) = 1. - elaz(:,nv)/tlaz(:,nv) fsnow(:) = min(max(fsnow(:),0.),1.) elsewhere fsnow(:) = 0. endwhere albdir(:,nv) = albvr(:)*(1.-fsnow(:)) + snovr(:)*fsnow(:) albdif(:,nv) = albvf(:)*(1.-fsnow(:)) + snovf(:)*fsnow(:) end do Done! (4) Put tile_id back to compute_rc, but not the subroutines therein yet. Will do that later when needed. Done! (5) Remove allocate and deallocate VEG_tmp from RUN1 and add them to RUN2. Ask Sarith what's the difference between RUN1 and RUN2? (will ask later) Need to add allocate and deallocate VEG_tmp to RUN1 as well? Yes. ========== 20170203: 1. The CLM4.5 test run saving monthly output still couldn't finish 12 months of simulation within 8 hours. It ended at "date_time_new 20001212_050730z". Will just work with the 11 months of output. Processed the output: /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/analysis > grid_restore_hdeg m0001hdeg_clm4.5 Made plots comparing FPAR, LAI and GPP between CLM4.5 and CLM4. Emailed results to Randy and Sarith to update them on the progress. 2. Worked on (1)-(5) in 2 on 20170202 above. 3. Build the executable: cd /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/src source g5_modules ./parallel_build.csh (chose realclean) 4. Set up an experiment: (1) Run gcm_setup to set up an experiment: cd /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/src/Applications/GEOSgcm_App ./gcm_setup Enter the Experiment ID: clm4.5 Enter a 1-line Experiment Description: Test run of clm4.5 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) c90 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) o2 Enter the choice of Land Surface Model: use 1 (Default: Catchment Model), 2 (CatchmentCN Model) 2 Do you wish to run the Runoff Routing Model? (Default: NO or FALSE) Do you wish to run GOCART with Actual or Climatological Aerosols? (Enter: A (Default) or C) Enter the GOCART Emission Files to use: MERRA2 (Default), PIESA, CMIP, NR, MERRA2-DD or OPS: 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: /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4.5 Enter Desired Location for the EXPERIMENT Directory (to contain model output and restart files) Hit ENTER to use Default Location: ---------------------------------- Default: /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4.5 Enter Location for Build directory containing: src/ Linux/ etc... Hit ENTER to use Default Location: ---------------------------------- Default: /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05 Current GROUPS: g0620 gmaoint m2val Enter your GROUP ID for Current EXP: (Default: g0620) ----------------------------------- sending incremental file list GEOSgcm.x sent 100189308 bytes received 31 bytes 200378678.00 bytes/sec total size is 100176983 speedup is 1.00 Creating gcm_run.j for Experiment: clm4.5 Creating gcm_post.j for Experiment: clm4.5 Creating gcm_archive.j for Experiment: clm4.5 Creating gcm_regress.j for Experiment: clm4.5 Creating gcm_convert.j for Experiment: clm4.5 Creating gcm_plot.tmpl for Experiment: clm4.5 Creating gcm_quickplot.csh for Experiment: clm4.5 Creating gcm_moveplot.j for Experiment: clm4.5 Creating gcm_stats.j for Experiment: clm4.5 Creating gcm_forecast.tmpl for Experiment: clm4.5 Creating gcm_forecast.setup for Experiment: clm4.5 Creating CAP.rc.tmpl for Experiment: clm4.5 Creating AGCM.rc.tmpl for Experiment: clm4.5 Creating HISTORY.rc.tmpl for Experiment: clm4.5 Creating fvcore_layout.rc for Experiment: clm4.5 Done! ----- Build Directory: /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05 ---------------- The following executable has been placed in your Experiment Directory: ---------------------------------------------------------------------- /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/Linux/bin/GEOSgcm.x You must now copy your Initial Conditions into: ----------------------------------------------- /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4.5 Copying Sandbox Source Code (including non-committed files) into /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4.5/src/GEOSagcm/src ... Enter CVS Tag to COMMIT this Experiment (Default: clm4.5__fzeng) Enter q or quit to QUIT or SKIP the CVS COMMIT: q 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. Could not chdir to home directory /home/fzeng: No such file or directory rcsmerge: warning: conflicts during merge cvs update: conflicts found in GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 Non-Committed Sandbox Files: --------------------------------- Applications/GEOSgcm_App/.AGCM_VERSION GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/CN_DriverMod.F90 file: /bin/cp: cannot stat `/discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4.5/src/GEOSagcm/src/file:': No such file or directory revision /bin/cp: cannot stat `/discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4.5/src/GEOSagcm/src/revision': No such file or directory revision /bin/cp: cannot stat `/discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4.5/src/GEOSagcm/src/revision': No such file or directory differences /bin/cp: cannot stat `/discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4.5/src/GEOSagcm/src/differences': No such file or directory GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/clm_varpar.F90 GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/compute_rc.F90 GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/pftvarcon.F90 GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts/GNUmakefile GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/update_model_paras.F90 file: /bin/cp: cannot stat `/discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4.5/src/GEOSagcm/src/file:': No such file or directory revision /bin/cp: cannot stat `/discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4.5/src/GEOSagcm/src/revision': No such file or directory revision /bin/cp: cannot stat `/discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4.5/src/GEOSagcm/src/revision': No such file or directory differences /bin/cp: cannot stat `/discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4.5/src/GEOSagcm/src/differences': No such file or directory GMAO_Shared/GEOS_Util/post/gcmpost.script Tarring Experiment Source Code into Single File ... (2) Copy and rename restart files: cd /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/ cp restarts/restarts.e20001230_21z.tar clm4.5/. cd clm4.5 tar -xf restarts.e20001230_21z.tar /bin/rm restarts.e20001230_21z.tar nedit move & move: ~~~~~~~~~~~~~~~~~~~~~~ mv test_clm4.agcm_import_rst.e20001230_21z.Heracles-5_4_p3.Heracles-4_3_MERRA-3_CF0090x6C_DE1440xPE0720 agcm_import_rst mv test_clm4.catchcn_internal_rst.e20001230_21z.Heracles-5_4_p3.Heracles-4_3_MERRA-3_CF0090x6C_DE1440xPE0720 catchcn_internal_rst mv test_clm4.fvcore_internal_rst.e20001230_21z.Heracles-5_4_p3.Heracles-4_3_MERRA-3_CF0090x6C_DE1440xPE0720 fvcore_internal_rst mv test_clm4.gocart_internal_rst.e20001230_21z.Heracles-5_4_p3.Heracles-4_3_MERRA-3_CF0090x6C_DE1440xPE0720 gocart_internal_rst mv test_clm4.lake_internal_rst.e20001230_21z.Heracles-5_4_p3.Heracles-4_3_MERRA-3_CF0090x6C_DE1440xPE0720 lake_internal_rst mv test_clm4.landice_internal_rst.e20001230_21z.Heracles-5_4_p3.Heracles-4_3_MERRA-3_CF0090x6C_DE1440xPE0720 landice_internal_rst mv test_clm4.moist_internal_rst.e20001230_21z.Heracles-5_4_p3.Heracles-4_3_MERRA-3_CF0090x6C_DE1440xPE0720 moist_internal_rst mv test_clm4.pchem_internal_rst.e20001230_21z.Heracles-5_4_p3.Heracles-4_3_MERRA-3_CF0090x6C_DE1440xPE0720 pchem_internal_rst mv test_clm4.saltwater_internal_rst.e20001230_21z.Heracles-5_4_p3.Heracles-4_3_MERRA-3_CF0090x6C_DE1440xPE0720 saltwater_internal_rst ~~~~~~~~~~~~~~~~~~~~~~ chmod 755 move move (3) Modify some files: Create cap_restart (which has only this line below) 20001230 210000 Edit CAP.rc END_DATE: 20010201 210000 (run 1 month for testing first) Edit HISTORY.rc In COLLECTIONS, comment out everything other than "geosgcm_surf" [Line 57] geosgcm_surf.frequency: 240000, [Lines 206-242] Uncomment geosgcm_surf collection for carbon variables (e.g., CNLAI, ICESOI) Edit gcm_run.j [line 7] #PBS -l walltime=06:00:00 [Line 161] setenv BCSDIR /discover/nobackup/ltakacs/bcs/Ganymed-4_0/Ganymed-4_0_MERRA-2 => setenv BCSDIR /discover/nobackup/smahanam/bcs/Heracles-4_3/Heracles-4_3_MERRA-3/ Comment out L1185 "qsub $FILE" in /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/src/GMAO_Shared/GEOS_Util/post/gcmpost.script to avoid submitting the post plotting job. sbatch gcm_run.j It finished successfully. ========== 20170209: 1. Procedures to check in (need to double check with Sarith when need to do this): see /discover/nobackup/fzeng/notes/cvs First check the difference between my sand box and the tag (the latter may be updated by others): LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/src > ~ltakacs/cvstools/cvscmp SM-LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3 Create a temporary file and record the files (including the paths) to be checked in. Then check in: (1) Create branch and decide what files to include in this branch? > cvs tag -b FZ_LandModeling file1 file2 ... fileN (filenames include path) (2) Update branch: > cvs upd -r FZ_LandModeling file1 file2 ... fileN (filenames include path) (3) Commit: cvs ci file1 file2 ... fileN (filenames include path) [NOTE: "ci" can be replaced by "commit"] (4) Check if files are committed: > tkcvs file1 (filename includes path) (5) Done and exit: > exit Files modified: (1) GEOScatchCN_GridComp: compute_rc.F90 process_cn.F90 in LDAS, and GEOS_CatchCNGridComp.F90 in GEOS5. pftvarcon.F90 clm_varpar.F90 CN_DriverMod.F90 in offline code (not this point driver): get_CN_LAI. (2) GEOSlana_GridComp: catch_types.F90 clsm_ensdrv_init_routines.F90 (3) GEOSland_GridComp/Shared/update_model_paras.F90 2. Sarith showed me how to do the above step by step: see 6 in see /discover/nobackup/fzeng/notes/cvs Sarith modifed GEOSland_GridComp/Shared/update_model_paras.F90 before checking it in. For my GEOS-5 sand box: Sarith updated a few routines that he had modified: cvs upd -r H54p3_v24C05_GOSWIM_zerodiff filenames These files include catchmentCN.F90, StieglitzSnow.F90, GOCART_GridComp.F90, AGCM.rc.tmp1, gcm_setup, MAPL_sun_uc.P90 (not F90), and a couple of more? ========== 20170210: 1. /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/src > ~ltakacs/cvstools/cvscmp H54p3_v24C05_GOSWIM_zerodiff The following files are locally MODIFIED: - GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/CN_DriverMod.F90 - GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 - GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/clm_varpar.F90 - GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/compute_rc.F90 - GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/pftvarcon.F90 - GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts/GNUmakefile - GMAO_Shared/GEOS_Util/post/gcmpost.script The following files need to be UPDATED for tag H54p3_v24C05_GOSWIM_zerodiff |-------------|------------------------------------------------------------ | Newest | |-------------| File Name | Here | Tag | |------|------|------------------------------------------------------------ | x | | GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/update_model_paras.F90 |------|------|------------------------------------------------------------ Compiled: /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/src > ./parallel_build.csh (chose realclean) Failed because there is a bug in GEOSland_GridComp/Shared/update_model_paras.F90: LocalTileID was saved twice while it should be saved only once Talked to Sarith. What he did after he fixed the bug is, in the offline code: (1) compiled "source g5_modules" and "gmake install" GEOSland_GridComp/Shared to make sure it compiles successfully (2) cvs status update_model_paras.F90 (Now update_model_paras.F90 is not updated in the repository yet) (3) module load other/tkcvs-8.2.3 (4) tkcvs update_model_paras.F90 (Now update_model_paras.F90 is not updated in the repository yet) (5) cvs ci update_model_paras.F90 (6) cvs tag -F FZ_CLM4toCLM45 update_model_paras.F90 Then in the GCM code: cvs upd -r FZ_CLM4toCLM45 update_model_paras.F90 ========== 20170213: 1. Compiled the GCM code. Finished successfully. 2. Followed steps in 4 on 20170203 to set up a test run clm4.5_20170213. Non-Committed Sandbox Files: --------------------------------- Applications/GEOSgcm_App/.AGCM_VERSION GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/CN_DriverMod.F90 GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/clm_varpar.F90 GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/compute_rc.F90 GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/pftvarcon.F90 GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts/GNUmakefile file: /bin/cp: cannot stat `/discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4.5_20170213/src/GEOSagcm/src/file:': No such file or directory revision /bin/cp: cannot stat `/discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4.5_20170213/src/GEOSagcm/src/revision': No such file or directory revision /bin/cp: cannot stat `/discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4.5_20170213/src/GEOSagcm/src/revision': No such file or directory differences /bin/cp: cannot stat `/discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/tests/clm4.5_20170213/src/GEOSagcm/src/differences': No such file or directory GMAO_Shared/GEOS_Util/post/gcmpost.script EXP: clm4.5_20170213 failed to be committed with TAG: clm4.5_20170213__fzeng ------------------------------------------------------- mv clm4.5 clm4.5_old rm -rf clm4_old It finished one month of simulation successfully. 3. Worked on subroutine Fractionation. ========== 20170214: 1. Continued to work on subroutine Fractionation. Where is downreg calculated? In /discover/nobackup/fzeng/clm4-to-clm4.5/clm4_5/main/histFldsMod.F90: call hist_addfld1d (fname='DOWNREG', units='proportion', & avgflag='A', long_name='fractional reduction in GPP due to N limitation', & ptr_pft=pepv%downreg, default='inactive') if ( use_c13 ) then call hist_addfld1d (fname='RC13_CANAIR', units='proportion', & avgflag='A', long_name='C13/C(12+13) for canopy air', & ptr_pft=pepv%rc13_canair) call hist_addfld1d (fname='RC13_PSNSUN', units='proportion', & avgflag='A', long_name='C13/C(12+13) for sunlit photosynthesis', & ptr_pft=pepv%rc13_psnsun) call hist_addfld1d (fname='RC13_PSNSHA', units='proportion', & avgflag='A', long_name='C13/C(12+13) for shaded photosynthesis', & ptr_pft=pepv%rc13_psnsha) endif subroutine hist_addfld1d is in /discover/nobackup/fzeng/clm4-to-clm4.5/clm4_5/main/histFileMod.F90. downreg is computed in CNAllocationMod.F90: excess_cflux(p) = availc(p) - plant_calloc(p) downreg(p) = excess_cflux(p)/gpp(p) Decided to add anything about carbon isotopes later. 2. Started to work on CN_DriverMod.F90. ========== 20170215: 1. Continued working on CN_DriverMod.F90: updated the data table in subroutine CN_init. 2. Started to work on clmtype.F90. Since not sure which variables should be removed at this time, just use the CLM4.5 clmtype.F90 as it is for now and will delete those not needed later. 3. Updated pftvarcon.F90. ========== 20170216: 1. Copied CLM4.5 clm_varctl.F90. Need to decide if this is really needed later. 2. Copied clmtypeInitMod.F90 and clm_varcon.F90 from CLM4.5. Will modify them later. 3. Modified clm_varpar.F90. ========== 20170217: 1. Compared all the CN routines between Catchment-CN and CCSM4_0 (the CLM4 version Greg implemented in Catchment-CN) to see how many modifications Greg made. This gives an idea on how many modifications I need to make to the CLM4.5 CN routines. 2. Read CNFireMod.F90. 3. Looked for the input data such as gdp etc. ========== 20170221: 1. The informataion below is from the global attributables of surfdata_1.9x2.5_mp24_simyr2000_c130419.nc and surfdata_0.125x0.125_mp24_simyr2000_c150114.nc in land01:/terra/fzeng/clm/surfdata_map, and CLM4.5 CNFireMod.F90. raw surfdata needed: (have to obtain from the CLM development people) ---------------------------------- agfirepkmon_raw_data_file_name = "mksrf_abm_0.5x0.5_AVHRR_simyr2000.c130201.nc" ; gdp_raw_data_file_name = "mksrf_gdp_0.5x0.5_AVHRR_simyr2000.c130228.nc" ; peatland_raw_data_file_name = "mksrf_peatf_0.5x0.5_AVHRR_simyr2000.c130228.nc" ; VOC_EF_raw_data_file_name = "mksrf_vocef_0.5x0.5_simyr2000.c110531.nc" ; (may not needed) new forcing data needed: (have to obtain from the CLM development people) ---------------------------------- human density data lightning frequency 2. perf_mod used in subroutine CNEcosystemDyn is located in: /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/utils/timing/perf_mod.F90 pftdynMod used in subroutine CNEcosystemDyn is located in main/pftdynMod.F90 surfrdMod used in subroutine CNEcosystemDyn is located in main/surfrdMod.F90 3. Replaced the CN routines in Catchment-CN with those in CLM4.5. From /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/src/Components/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp: mv CNAllocationMod.F90 clm4/CNAllocationMod.F90.clm4 mv CNAnnualUpdateMod.F90 clm4/CNAnnualUpdateMod.F90.clm4 mv CNBalanceCheckMod.F90 clm4/CNBalanceCheckMod.F90.clm4 mv CNCStateUpdate1Mod.F90 clm4/CNCStateUpdate1Mod.F90.clm4 mv CNCStateUpdate2Mod.F90 clm4/CNCStateUpdate2Mod.F90.clm4 mv CNCStateUpdate3Mod.F90 clm4/CNCStateUpdate3Mod.F90.clm4 mv CNDecompMod.F90 clm4/CNDecompMod.F90.clm4 mv CNFireMod.F90 clm4/CNFireMod.F90.clm4 mv CNGapMortalityMod.F90 clm4/CNGapMortalityMod.F90.clm4 mv CNGRespMod.F90 clm4/CNGRespMod.F90.clm4 mv CNHarvestMod.F90 clm4/CNHarvestMod.F90.clm4 mv CNiniTimeVar.F90 clm4/CNiniTimeVar.F90.clm4 mv CNMRespMod.F90 clm4/CNMRespMod.F90.clm4 mv CNNDynamicsMod.F90 clm4/CNNDynamicsMod.F90.clm4 mv CNNStateUpdate1Mod.F90 clm4/CNNStateUpdate1Mod.F90.clm4 mv CNNStateUpdate2Mod.F90 clm4/CNNStateUpdate2Mod.F90.clm4 mv CNNStateUpdate3Mod.F90 clm4/CNNStateUpdate3Mod.F90.clm4 mv CNPrecisionControlMod.F90 clm4/CNPrecisionControlMod.F90.clm4 mv CNSetValueMod.F90 clm4/CNSetValueMod.F90.clm4 mv CNSummaryMod.F90 clm4/CNSummaryMod.F90.clm4 mv CNVegStructUpdateMod.F90 clm4/CNVegStructUpdateMod.F90.clm4 mv CNWoodProductsMod.F90 clm4/CNWoodProductsMod.F90.clm4 mv nanMod.F90 clm4/nanMod.F90.clm4 mv clm_time_manager.F90 clm4/clm_time_manager.F90.clm4 cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/ch4Mod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/ch4RestMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/ch4varcon.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNAllocationMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNAnnualUpdateMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNBalanceCheckMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNC14DecayMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNCIsoFluxMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNCStateUpdate1Mod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNCStateUpdate2Mod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNCStateUpdate3Mod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNDecompCascadeMod_BGC.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNDecompCascadeMod_CENTURY.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNDecompMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNDVEcosystemDynIniMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNDVEstablishmentMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNDVLightMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNDVMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNEcosystemDynMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNFireMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNGapMortalityMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNGRespMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNMRespMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNNDynamicsMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNNitrifDenitrifMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNNStateUpdate1Mod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNNStateUpdate2Mod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNNStateUpdate3Mod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNPhenologyMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNPrecisionControlMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNrestMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNSetValueMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNSoilLittVertTranspMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNSummaryMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNVegStructUpdateMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNVerticalProfileMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNWoodProductsMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CropRestMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/DryDepVelocity.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/DUSTMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/initch4Mod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/MEGANFactorsMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/STATICEcosysDynMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/VOCEmissionMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/main/pftdynMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/main/surfrdMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/utils/timing/perf_mod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/util_share/clm_time_manager.F90 . 4. In CNNDynamicsMod.F90: gridcell =>col%gridcell (was gridcell => clm3%g%l%c%gridcell in CLM4) ndep_to_sminn => cnf%ndep_to_sminn (was ndep_to_sminn => clm3%g%l%c%cnf%ndep_to_sminn in CLM4) Need "qflx_surf(:) ! sub-surface runoff (mm H2O /s)" from the Catchment model for subroutine CNNLeaching. Also need "dz(:,:) !layer thickness (m)" Followed qflux_drain to add qflx_surf to CN_DriverMod.F90 (NOT finished yet!!) bflow in process_cn.F90 is passed to CN_DriverMod.F90 as qflx_drain. Why is bflow calculated in process_cn.F90 in this way below, instead of just using cat_diagF%bflow? bflow(n) = (1.-frice)*1000.* & cat_param(n)%cond*exp(-(cat_param(n)%bf3-ashift)-cat_param(n)%gnu*zbar)/cat_param(n)%gnu IF(cat_progn(n)%catdef >= cat_param(n)%cdcr1) bflow(n) = 0. bflow(n) = min(cat_param(n)%cond,bflow(n)) ========== 20170222: 1. Talk to Sarith about questions above! Calculation of bflow in process_cn.F90 for the CN model is from subroutine BASE in lsm_routines.F90. Calculation of surface runoff in process_cn.F90 for the CN model should follow subroutine SRUNOFF in catchmentCN.F90. It's more complicated than computing bflow because RUNSRF is an inout argument in subroutine SRUNOFF in catchmentCN. ar1, ar2 and ar4 have already been calculated in "call catch_calc_soil_moist" in process_cn.F90, so no need to worry about PARTITION which calculates ar1, ar2 and ar4. dz is already in CN_DriverMod.F90. 2. Added values for parameter fertnitro (used to compute fert nitrogen fertilizer rate (gN/m2/s)) in CN_DriverMod.F90. fert is used in subroutine CNNFert. Conversion from fertnitro to fert is done in CNPhenologyMod.F90. 3. subroutine CNSoyfix needs hui (=gdd since planting (gddplant)) but couldn't find where it's computed so far. ========== 20170223: 1. For calculation of runsrf (surface runoff) in process_cn.F90, followed Sarith's suggestion on using the output from "call CATCHCN" for now. By doing so, the first time step of every job segment will have 0 runsrf. Will improve this later. 2. Found gddplant which is passed onto hui used in subroutine CNSoyfix. gddplant is calculated in accFldsMod.F90. 3. Added values for these parameters to CN_DriverMod.F90: pftcon%lfemerg pftcon%grnfill pftcon%mxmat ========== 20170224: 1. cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/main/clm_initializeMod.F90 . 2. Added values for these parameters to CN_DriverMod.F90: pftcon%hybgdd 3. Where and how are gdd0, gdd8, gdd10, a5tmin, a10tmin, t_ref2m_min, and gddtsoi used in CNPhenologyClimate calculated? Check main/accFldsMod.F90 and /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/util_share/accumulMod.F90. May need these two routines. L191 in CNPheologyMod.F90: call decomp_rate_constants snow_depth is a field in a type now (see subroutine vernalization in CNPhenologyMod.F90). ========== 20170227: 1. Checked decomp_rate_constants in CNPheologyMod.F90. The constants are already set in the subroutine. 2. Calculation of gdd0, gdd8, gdd10, a5tmin, a10tmin, t_ref2m_min, t_ref2m_max, and gddtsoi is in these two modules: /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/main/accFldsMod.F90 /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/util_share/accumulMod.F90 Copied them. cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/main/accFldsMod.F90 . cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/util_share/accumulMod.F90 . 3. Subroutine CNHarvest is now in pftdynMod.F90. It needs variable harvest(n_cat), which is calculated in subroutine pftdyn_getharvest in pftdynMod.F90. See below. call ncd_io(ncid=ncid, varname= 'HARVEST_VH1', flag='read', data=arrayl, dim1name=grlnd, & nt=ntime, readvar=readvar) if (.not. readvar) call endrun( trim(subname)//' ERROR: HARVEST_VH1 not on pftdyn file' ) harvest(begg:endg) = arrayl(begg:endg) call ncd_io(ncid=ncid, varname= 'HARVEST_VH2', flag='read', data=arrayl, dim1name=grlnd, & nt=ntime, readvar=readvar) if (.not. readvar) call endrun( trim(subname)//' ERROR: HARVEST_VH2 not on pftdyn file' ) harvest(begg:endg) = harvest(begg:endg) + arrayl(begg:endg) call ncd_io(ncid=ncid, varname= 'HARVEST_SH1', flag='read', data=arrayl, dim1name=grlnd, & nt=ntime, readvar=readvar) if (.not. readvar) call endrun( trim(subname)//' ERROR: HARVEST_SH1 not on pftdyn file' ) harvest(begg:endg) = harvest(begg:endg) + arrayl(begg:endg) call ncd_io(ncid=ncid, varname= 'HARVEST_SH2', flag='read', data=arrayl, dim1name=grlnd, & nt=ntime, readvar=readvar) if (.not. readvar) call endrun( trim(subname)//' ERROR: HARVEST_SH2 not on pftdyn file' ) harvest(begg:endg) = harvest(begg:endg) + arrayl(begg:endg) call ncd_io(ncid=ncid, varname= 'HARVEST_SH3', flag='read', data=arrayl, dim1name=grlnd, & nt=ntime, readvar=readvar) if (.not. readvar) call endrun( trim(subname)//' ERROR: HARVEST_SH3 not on pftdyn file' ) harvest(begg:endg) = harvest(begg:endg) + arrayl(begg:endg) Where is "ncid" from? In subroutine pftdyn_init: ! Obtain file call getfil (fpftdyn, locfn, 0) call ncd_pio_openfile (ncid, locfn, 0) Where is subroutine getfil? models/atm/cam/src/control/ioFileMod.F90 37:subroutine getfil(fulpath, locfn, iflag, lexist) 131:end subroutine getfil models/lnd/clm/tools/shared/mkprocdata_map/src/fileutils.F90 108: subroutine getfil (fulpath, locfn, iflag) 172: end subroutine getfil models/lnd/clm/tools/clm4_5/mksurfdata_map/src/fileutils.F90 78: subroutine getfil (fulpath, locfn, iflag) 140: end subroutine getfil models/lnd/clm/src/util_share/fileutils.F90 79: subroutine getfil (fulpath, locfn, iflag) 141: end subroutine getfil models/rof/rtm/src/riverroute/RtmFileUtils.F90 59: subroutine getfil (fulpath, locfn, iflag) 111: end subroutine getfil Is the one in models/lnd/clm/tools/clm4_5/mksurfdata_map/src/fileutils.F90 being used here? Also in subroutine pftdyn_init: use clm_varctl , only : fpftdyn ========== 20170228: 1. About calling CNHarvest: In CNEcosystemDynMod.F90: if (fpftdyn /= ' ') then call CNHarvest(num_soilc, filter_soilc, num_soilp, filter_soilp) end if and in clm_varctl.F90: character(len=256), public :: fpftdyn = ' ' ! dynamic landuse dataset So CNHarvest is actually not called at the time. So no need to worry about this now. 2. CNFireMod.F90: prec10 and prec60 are computed by accFldsMod.F90 and accumulMod.F90. prec10_col and prec60_col are computed from prec10 and prec60 in CNFireMod.F90. Looked for lightning data. In CLM4.5 Tech Note, it says: Climatological 3-hourly lightning frequency at ~1.8o resolution is provided, which was calculated via bilinear interpolation from 1995-2011 NASA LIS/OTD grid product v2.2 (http://ghrc.msfc.nasa.gov) 2-hourly, 2.5o lightning frequency data. In future versions of the model, lightning data may be obtained directly from the atmosphere model. ========== 20170301: 1. Continued looking for lightning data. https://ghrc.nsstc.nasa.gov/hydro/ -> "Dataset List" on the left -> "Lightning Products": LIS 0.1 Degree Very High Resolution Gridded Lightning Annual Climatology (VHRAC) 38S-38N only LIS 0.1 Degree Very High Resolution Gridded Lightning Diurnal Climatology (VHRDC) 38S-38N only LIS 0.1 Degree Very High Resolution Gridded Lightning Full Climatology (VHRFC) 38S-38N only LIS 0.1 Degree Very High Resolution Gridded Lightning Monthly Climatology (VHRMC) 38S-38N only LIS 0.1 Degree Very High Resolution Gridded Lightning Seasonal Climatology (VHRSC) 38S-38N only LIS/OTD 0.5 Degree High Resolution Annual Climatology (HRAC) LIS/OTD 0.5 Degree High Resolution Full Climatology (HRFC) LIS/OTD 0.5 Degree High Resolution Monthly Climatology (HRMC) LIS/OTD 2.5 Degree Low Res Annual Diurnal Climatology (LRADC) LIS/OTD 2.5 Degree Low Resolution Annual Climatology (LRAC) LIS/OTD 2.5 Degree Low Resolution Annual Climatology Time Series LIS/OTD 2.5 Degree Low Resolution Diurnal Climatology (LRDC) LIS/OTD 2.5 Degree Low Resolution Full Climatology (LRFC) LIS/OTD 2.5 Degree Low Resolution Monthly Time Series (LRMTS) LIS/OTD 2.5 Degree Low Resolution Time Series (LRTS) On this webpage https://lightning.nsstc.nasa.gov/data/data_lis-otd-climatology.html it says: Users are strongly encouraged to choose the LIS/OTD 2.5 Degree Low Resolution Monthly Time Series (LRMTS) data set over the LIS/OTD 2.5 Degree Low Resolution Time Series (LRTS) data set. Both products have ~3 month smoothing, so the 'daily data' adds little useful information compared to that in the much smaller LRMTS files. So it looks like the finest temporal resolution we can get is monthly. 2. This routine is needed in CNFireMod.F90. Copied it to GEOScatchCN_GridComp. /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/csm_share/shr/shr_strdata_mod.F90 ========== 20170302: 1. Emailed support-ghrc@earthdata.nasa.gov for the 1995-2011 NASA LIS/OTD grid product v2.2 (http://ghrc.msfc.nasa.gov) 2-hourly 2.5 Degree lightning frequency data that NCAR used. 2. Read about QFED (Quick Fire Emission Dataset) and emailed Sarith. 3. btran2 (root zone soil wetness) needed in CNFireMod.F90 is new. It's computed in CanopyFluxesMod.F90. do j = 1,nlevgrnd do f = 1, fn p = filterp(f) c = pcolumn(p) l = plandunit(p) smp_node_lf = max(smpsc(ivt(p)), -sucsat(c,j)*(h2osoi_vol(c,j)/watsat(c,j))**(-bsw(c,j))) btran2(p) = btran2(p) +rootfr(p,j)*min((smp_node_lf - smpsc(ivt(p))) / (smpso(ivt(p)) - smpsc(ivt(p))), 1._r8) real(r8), pointer :: sucsat(:,:) ! minimum soil suction (mm) real(r8), pointer :: bsw(:,:) ! Clapp and Hornberger "b" real(r8), pointer :: watsat(:,:) ! volumetric soil water at saturation (porosity) real(r8), pointer :: h2osoi_vol(:,:) ! volumetric soil water (0<=h2osoi_vol<=watsat) [m3/m3] by F. Li and S. Levis real(r8), pointer :: rootfr(:,:) ! fraction of roots in each soil layer real(r8), pointer :: smpso(:) ! soil water potential at full stomatal opening (mm) real(r8), pointer :: smpsc(:) ! soil water potential at full stomatal closure (mm) sucsat => cps%sucsat bsw => cps%bsw watsat => cps%watsat h2osoi_vol => cws%h2osoi_vol rootfr => pps%rootfr smpsc and smpso are pft-dependent parameters. Values are in /discover/nobackup/fzeng/clm4-to-clm4.5/data/pftdata4.5/pft-physiology.c130503.nc j = 1 in our case. sucsat = cat_param%psis in Catchment watsat = cat_param%poros in Catchment bsw = cat_param%bee in Catchment h2osoi_vol = cat_diagS%rzmc in Catchment Need to check with Randy about this! Confirmed by Randy. For sucsat, in CLM4.5 main/iniTimeConst.F90: 131: real(r8), pointer :: sucsat(:,:) ! minimum soil suction (mm) (nlevgrnd) 348: sucsat => cps%sucsat 1179: sucsat(c,lev) = 10._r8 * ( 10._r8**(1.88_r8-0.0131_r8*sand) ) 1182: om_sucsat = min(10.3_r8 - 0.2_r8*(zsoi(lev)/zsapric), 10.1_r8) 1189: sucsat(c,lev) = (1._r8-om_frac)*sucsat(c,lev) + om_sucsat*om_frac Will ask Randy. For rootfr: in CLM4.5 main/iniTimeConst.F90: 14:! 3) rootfr only initialized for soil points 121: real(r8), pointer :: rootfr(:,:) ! fraction of roots in each soil layer 122: real(r8), pointer :: rootfr_road_perv(:,:) ! fraction of roots in each soil layer for urban pervious road 346: rootfr_road_perv => cps%rootfr_road_perv 398: rootfr => pps%rootfr 1237: ! pervious road layers -- same as above except also set rootfr_road_perv 1241: rootfr_road_perv(c,lev) = 0._r8 1244: rootfr_road_perv(c,lev) = 0.1_r8 ! uniform profile 1340: rootfr(p,lev) = 0._r8 1343: rootfr(p,lev) = .5_r8*( exp(-roota_par(ivt(p)) * zi(c,lev-1)) & 1348: rootfr(p,nlevsoi) = .5_r8*( exp(-roota_par(ivt(p)) * zi(c,nlevsoi-1)) & 1350: rootfr(p,nlevsoi+1:nlevgrnd) = 0.0_r8 1363:! rootfr(p,lev) = dz(c,lev) * 0.5_r8 * ((intercept+slope*zi(c,lev-1)) + (intercept+slope*zi(c,lev))) 1367:! rootfr(p,lev) = 0._r8 1375:! rootfr(p,lev) = dz(c,lev) * 0.5_r8 * ((intercept+slope*zi(c,lev-1)) + (intercept+slope*zi(c,lev))) 1379:! rootfr(p,lev) = 0._r8 1385: rootfr(p,1:nlevsoi) = 0._r8 in CLM4.5 biogeochem/CNMRespMod.F90: 72: real(r8), pointer :: rootfr(:,:) ! fraction of roots in each soil layer (nlevgrnd) 115: rootfr => pps%rootfr 190: ! rootfr(j) sums to 1.0 over all soil layers, and 196: froot_mr(p) = froot_mr(p) + frootn(p)*br*tcsoi(c,j)*rootfr(p,j) Does this mean that rootfr is 1.0 if there is only one soil layer in our case? Yes. In CN_DriverMod.F90: rootfr(p,1) = 1.0 ! gkw: affects maint resp. test sensitivity ========== 20170303: 1. Continued working on 3 on 20170302. ========== 20170306: 1. fsat is needed in CNFireMod.F90. It's computed in biogeophys/SoilHydrologyMod.F90 159: real(r8), pointer :: fsat(:) ! fractional area with water table at surface 229: fsat => cws%fsat 291: fsat(c) = A(c) !for output 293: fsat(c) = wtfact(c) * exp(-0.5_r8*fff(c)*zwt(c)) 296: ! use perched water table to determine fsat (if present) 299: fsat(c) = A(c) 301: fsat(c) = wtfact(c) * exp(-0.5_r8*fff(c)*zwt(c)) 305: fsat(c) = wtfact(c) * exp(-0.5_r8*fff(c)*zwt_perched(c))!*( frost_table(c) - zwt_perched(c))/4.0 312: fcov(c) = (1._r8 - fracice(c,1)) * fsat(c) + fracice(c,1) 315: fcov(c) = fsat(c) 327: qflx_surf(c) = fsat(c) * qflx_top_soil(c) 422: real(r8), pointer :: fsat(:) ! fractional area with water table at surface 516: fsat => cws%fsat 620: qinmax = (1._r8 - fsat(c)) * 10._r8**(-e_ice*top_icefrac)*(qflx_in_soil(c) - rsurf_vic) 622: qinmax=(1._r8 - fsat(c)) * minval(10._r8**(-e_ice*(icefrac(c,1:3)))*hksat(c,1:3)) Will ask Randy. 2. In CNFireMod: real(r8), pointer :: wf(:) ! soil water as frac. of whc for top 0.05 m real(r8), pointer :: wf2(:) ! soil water as frac. of whc for top 0.17 m wf is in CLM4, so no need to worry about it. wf2 is new. 3. Talked to Randy about lightning data, sucsat for btran2 calculation, fsat, wf and wf2 needed in CNFireMod. Added to CN_DriverMod: wf2(n) = rzm(nc,nz)/poros(nc) Added to process_cn: car1m and its calculation Added to CN_DriverMod: car1m and calculation of fsat ========== 20170307: 1. btran and btran2 in CLM4.5: what's the difference? In CanopyFluxesMod.F90: real(r8), parameter :: btran0 = 0.0_r8 ! initial value real(r8), pointer :: rootr(:,:) ! effective fraction of roots in each soil layer ! Initialize btran(p) = btran0 btran2(p) = btran0 do j = 1,nlevgrnd do f = 1, fn p = filterp(f) c = pcolumn(p) l = plandunit(p) ! Root resistance factors vol_ice = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice)) eff_porosity = watsat(c,j)-vol_ice vol_liq = min(eff_porosity, h2osoi_liq(c,j)/(dz(c,j)*denh2o)) if (vol_liq .le. 0._r8 .or. t_soisno(c,j) .le. tfrz-2._r8) then rootr(p,j) = 0._r8 else s_node = max(vol_liq/eff_porosity,0.01_r8) smp_node = max(smpsc(ivt(p)), -sucsat(c,j)*s_node**(-bsw(c,j))) rresis(p,j) = min( (eff_porosity/watsat(c,j))* & (smp_node - smpsc(ivt(p))) / (smpso(ivt(p)) - smpsc(ivt(p))), 1._r8) if (.not. (perchroot .or. perchroot_alt) ) then rootr(p,j) = rootfr(p,j)*rresis(p,j) else rootr(p,j) = rootfr_unf(p,j)*rresis(p,j) end if btran(p) = btran(p) + rootr(p,j) smp_node_lf = max(smpsc(ivt(p)), -sucsat(c,j)*(h2osoi_vol(c,j)/watsat(c,j))**(-bsw(c,j))) btran2(p) = btran2(p) +rootfr(p,j)*min((smp_node_lf - smpsc(ivt(p))) / (smpso(ivt(p)) - smpsc(ivt(p))), 1._r8) endif end do end do Emailed the info to Randy. 2. Downloaded LISOTD_LRMTS_V2.3.2014.hdf, LISOTD_HRMC_V2.3.2014.hdf, and LISOTD_HRAC_V2.3.2014.hdf from https://ghrc.nsstc.nasa.gov/pub/lis/climatology/LIS-OTD/. File info from the HDF Import Tool of MATLAB (by dragging each file into MATLAB command window): LISOTD_LRMTS_V2.3.2014.hdf: LIS/OTD 2.5 Degree Low Resolution Monthly Time Series (LRMTS) Name: LRMTS_COM_FR Dimensions: Name: Latitude Size: 72 Name: Longitude Size: 144 Name: Month since Jan 95 Size: 240 Precision: float valid_range: 0 0.27408 long_name: Combined Flash Rate Monthly Time Series units: Flashes/km^2/day _FillValue: 0 Period: 1995-2014 LISOTD_HRMC_V2.3.2014.hdf: LIS/OTD 0.5 Degree High Resolution Monthly Climatology (HRMC) Name: HRMC_COM_FR Dimensions: Name: Latitude Size: 360 Name: Longitude Size: 720 Name: Month Size: 12 Precision: float valid_range: 0 1.1855 long_name: Combined Flash Rate Monthly Climatology units: Flashes/km^2/day _FillValue: 0 Averaged over 1995-2013 LISOTD_HRAC_V2.3.2014.hdf: LIS/OTD 0.5 Degree High Resolution Annual Climatology (HRAC) Name: HRAC_COM_FR Dimensions: Name: Latitude Size: 360 Name: Longitude Size: 720 Name: Day of year Size: 365 Precision: float valid_range: 0 2.8115 long_name: Combined Flash Rate Annual Climatology units: Flashes/km^2/day _FillValue: 0 Averaged over 1995-2013 This webpage describes the lightning data sets in a better way for us to pick the one that is most suitable for us. https://ghrc.nsstc.nasa.gov/uso/ds_docs/lis_climatology/LISOTD_climatology_dataset.html The time period used in the climatology calculation can be found by clicking on the link under each data set on this webpage: https://lightning.nsstc.nasa.gov/data/data_lis-otd-climatology.html Emailed Randy and asked which one works best for us. ========== 20170308: 1. tsoi17 needed in CNFireMod is calculated in beogeophys/Hydrology2Mod: ! column-level real(r8), pointer :: tsoi17(:) ! soil temperature in top 17cm of soil (Kelvin) added by F. Li and S. Levis ! Determine ground temperature, ending water balance and volumetric soil water ! Calculate soil temperature and total water (liq+ice) in top 10cm of soil ! Calculate soil temperature and total water (liq+ice) in top 17cm of soil do fc = 1, num_nolakec c = filter_nolakec(fc) l = clandunit(c) if (ityplun(l) /= isturb) then t_soi_10cm(c) = 0._r8 tsoi17(c) = 0._r8 h2osoi_liqice_10cm(c) = 0._r8 end if end do do j = 1, nlevsoi do fc = 1, num_nolakec c = filter_nolakec(fc) l = clandunit(c) if (ityplun(l) /= isturb) then ! soil T at top 17 cm added by F. Li and S. Levis if (zi(c,j) <= 0.17_r8) then fracl = 1._r8 tsoi17(c) = tsoi17(c) + t_soisno(c,j)*dz(c,j)*fracl else if (zi(c,j) > 0.17_r8 .and. zi(c,j-1) .lt. 0.17_r8) then fracl = (0.17_r8 - zi(c,j-1))/dz(c,j) tsoi17(c) = tsoi17(c) + t_soisno(c,j)*dz(c,j)*fracl end if end if if (zi(c,j) <= 0.1_r8) then fracl = 1._r8 t_soi_10cm(c) = t_soi_10cm(c) + t_soisno(c,j)*dz(c,j)*fracl h2osoi_liqice_10cm(c) = h2osoi_liqice_10cm(c) + & (h2osoi_liq(c,j)+h2osoi_ice(c,j))* & fracl else if (zi(c,j) > 0.1_r8 .and. zi(c,j-1) .lt. 0.1_r8) then fracl = (0.1_r8 - zi(c,j-1))/dz(c,j) t_soi_10cm(c) = t_soi_10cm(c) + t_soisno(c,j)*dz(c,j)*fracl h2osoi_liqice_10cm(c) = h2osoi_liqice_10cm(c) + & (h2osoi_liq(c,j)+h2osoi_ice(c,j))* & fracl end if end if end if end do end do if (ityplun(l)==isturb) then t_grnd_u(c) = t_soisno(c,snl(c)+1) else t_soi_10cm(c) = t_soi_10cm(c)/0.1_r8 tsoi17(c) = tsoi17(c)/0.17_r8 ! F. Li and S. Levis end if ========== 20170309: 1. Use btran2 = btran, and tsoi17 = tgw for now. Need to improve these later!! 2. Subroutine CNFireFluxes: (1) wtcol: already in CN_DriverMod 282: pwtcol => clm3%g%l%c%p%wtcol 473: pwtcol(p) = 0. ! weight will be zero unless otherwise set 483: pwtcol(p) = bare 501: pwtcol(p) = fveg(nc,nv,nz) (2) latdeg and londeg: added calculation to CN_DriverMod (3) Added the following pft-dependent parameters and their values to CN_DriverMod cc_dstem,cc_leaf,cc_lstem,cc_other fm_leaf,fm_lstem,fm_other,fm_root,fm_lroot,fm_droot Allocated and initialized them in clmtypeInitMod.F90. Modified subroutine CNFireFluxes to use these parameters as well as the following six parameters: lf_flab,lf_fcel,lf_flig,fr_flab,fr_fcel,fr_flig 3. Subroutine CNVegStructUpdate: (1) pft-dependent parameters: laimx, ztopmx Added to CN_DriverMod Added to pftcon type in clmtype.F90 Allocated and initialized in clmtypeInitMod.F90. Modified subroutine CNVegStructUpdate to use these parameters. (2) dwood: Added to CN_DriverMod. Values are from CLM4.5 pftvarcon.F90. ========== 20170310: 1. Subroutine decomp_rate_constants in CNDecompCascadeMod_BGC.F90 and CNDecompCascadeMod_CENTURY.F90: (1) nlev_soildecomp_standard=1 ! 5 originally in CLM4.5, changed to 1 for now because we have only 1 hydrologically active soil layer (1m thick), fzeng, 10 Mar 2017 (2) sucsat is needed. Added to CN_DriverMod. (3) alt_indx (current depth of thaw) is needed. biogeophys/ActiveLayerMod.F90: do fc = 1,num_soilc c = filter_soilc(fc) ! calculate alt for a given timestep ! start from base of soil and search upwards for first thawed layer. ! note that this will put talik in with active layer ! a different way of doing this could be to keep track of how long a given layer has ben frozen for, and define ALT as the first layer that has been frozen for less than 2 years. if (t_soisno(c,nlevgrnd) > SHR_CONST_TKFRZ ) then alt(c) = zsoi(nlevgrnd) alt_indx(c) = nlevgrnd else k_frz=0 found_thawlayer = .false. do j=nlevgrnd-1,1,-1 if ( ( t_soisno(c,j) > SHR_CONST_TKFRZ ) .and. .not. found_thawlayer ) then k_frz=j found_thawlayer = .true. endif end do if ( k_frz .gt. 0 ) then ! define active layer as the depth at which the linearly interpolated temperature line intersects with zero z1 = zsoi(k_frz) z2 = zsoi(k_frz+1) t1 = t_soisno(c,k_frz) t2 = t_soisno(c,k_frz+1) alt(c) = z1 + (t1-SHR_CONST_TKFRZ)*(z2-z1)/(t1-t2) alt_indx(c) = k_frz else alt(c)=0._r8 alt_indx(c) = 0 endif endif Currently in our clm_varpar.F90: nlevgrnd = 1 ! number of ground layers (includes lower layers that are hydrologically inactive) 2. Subroutine decomp_vertprofiles in CNVerticalProfileMod.F90: use clm_varcon, only: zsoi, dzsoi, zisoi in clm_varcon.F90: real(r8), pointer :: zsoi(:) !soil z (layers) real(r8), pointer :: dzsoi(:) !soil dz (thickness) real(r8), pointer :: zisoi(:) !soil zi (interfaces) They are calculated in CLM4.5 main/iniTimeConst.F90. ========== 20170313: 1. Continued 2 on 20170310. In subroutine clm_varcon_init in clm_varcon.F90: allocate( zsoi(1:nlevgrnd) ) allocate( dzsoi(1:nlevgrnd) ) allocate( zisoi(0:nlevgrnd) ) allocate( dzsoi_decomp(1:nlevdecomp_full) ) In main/iniTimeConst.F90: ! -------------------------------------------------------------------- ! Define layer structure for soil, lakes, urban walls and roof ! Vertical profile of snow is not initialized here ! -------------------------------------------------------------------- ! Soil layers and interfaces (assumed same for all non-lake patches) ! "0" refers to soil surface and "nlevsoi" refers to the bottom of model soil if ( more_vertlayers )then ! replace standard exponential grid with a grid that starts out exponential, then has several evenly spaced layers, then finishes off exponential. ! this allows the upper soil to behave as standard, but then continues with higher resolution to a deeper depth, so that, for example, permafrost ! dynamics are not lost due to an inability to resolve temperature, moisture, and biogeochemical dynamics at the base of the active layer do j = 1, toplev_equalspace zsoi(j) = scalez*(exp(0.5_r8*(j-0.5_r8))-1._r8) !node depths enddo do j = toplev_equalspace+1,toplev_equalspace + nlev_equalspace zsoi(j) = zsoi(j-1) + thick_equal enddo do j = toplev_equalspace + nlev_equalspace +1, nlevgrnd zsoi(j) = scalez*(exp(0.5_r8*((j - nlev_equalspace)-0.5_r8))-1._r8) + nlev_equalspace * thick_equal enddo else do j = 1, nlevgrnd zsoi(j) = scalez*(exp(0.5_r8*(j-0.5_r8))-1._r8) !node depths enddo end if dzsoi(1) = 0.5_r8*(zsoi(1)+zsoi(2)) !thickness b/n two interfaces do j = 2,nlevgrnd-1 dzsoi(j)= 0.5_r8*(zsoi(j+1)-zsoi(j-1)) enddo dzsoi(nlevgrnd) = zsoi(nlevgrnd)-zsoi(nlevgrnd-1) zisoi(0) = 0._r8 do j = 1, nlevgrnd-1 zisoi(j) = 0.5_r8*(zsoi(j)+zsoi(j+1)) !interface depths enddo zisoi(nlevgrnd) = zsoi(nlevgrnd) + 0.5_r8*dzsoi(nlevgrnd) In CLM4.5 clm_varpar.F90: logical, public :: more_vertlayers = .false. ! true => run with more vertical soil layers so do j = 1, nlevgrnd zsoi(j) = scalez*(exp(0.5_r8*(j-0.5_r8))-1._r8) !node depths enddo In CLM4.5 biogeochem/CNSoilLittVertTranspMod.F90: Looks like zsoi of multiple soil layer is only used when VERTSOILC is defined (i.e. turned on). Temporary solution: (1) Turn off VERTSOILC and (2) set these in clm_varcon.F90: zsoi(1) = 0.5 dzsoi(1) = 1. zisoi(0) = 0. zisoi(1) = 1. dzsoi_decomp(1) = dzsoi(1) Ask Randy later about these!! 2. Subroutine nitrif_denitrif in CNNitrifDenitrifMod: (1) Added nlevdecomp and nlevdecomp_full to clm_varpar.F90 (followed CLM4.5 clm_varpar.F90). (2) Needs bd and watsat: real(r8), pointer :: bd(:,:) ! bulk density of dry soil material [kg/m3] bd => cps%bd do j = 1, nlevdecomp do fc = 1,num_soilc c = filter_soilc(fc) soil_bulkdensity(c,j) = bd(c,j) + h2osoi_liq(c,j)/dz(c,j) so bd(ncolumns,nlevdecomp). main/iniTimeConst.F90 126: real(r8), pointer :: watsat(:,:) ! volumetric soil water at saturation (porosity) (nlevgrnd) 142: real(r8), pointer :: bd(:,:) ! bulk density of dry soil material [kg/m^3] 342: watsat => cps%watsat 356: bd => cps%bd 1177: watsat(c,lev) = 0.489_r8 - 0.00126_r8*sand 1185: bd(c,lev) = (1._r8-watsat(c,lev))*2.7e3_r8 watsat = cat_param%poros in Catchment (see notes for 20170302 above). Added calculation of bd to CN_DriverMod. ========== 20170315: 1. Continued working on subroutine nitrif_denitrif in CNNitrifDenitrifMod: (3) Needs tmean_monthly_max_vr and tmean_monthly_vr: real(r8), pointer :: tmean_monthly_max_vr(:,:) ! maximumn monthly-mean soil temperature real(r8), pointer :: tmean_monthly_vr(:,:) ! monthly-mean soil temperature No. They are declared but not used at all in CNNitrifDenitrifMod. (4) Needs watfc and bsw: real(r8), pointer :: watfc(:,:) ! volumetric soil water at field capacity (nlevsoi) real(r8), pointer :: bsw(:,:) ! Clapp and Hornberger "b" (nlevgrnd) bsw => cps%bsw bsw = cat_param%bee in Catchment (see notes on 20170302) Which variable in Catchment is this "watfc (volumetric soil water at field capacity)"? Will ask Sarith. Sarith does not know about this. He said our model does not use field capacity. Need to ask Randy. Since we will turn on NITRIF_DENITRIF at the early testing, don't worry about this now. Will ask Randy later. Asked Gabrielle on 3/17/2017 and she replied on 3/18/2017: "Indeed, by "field capacity" we mean "volumetric water at field capacity". Both are the same." Sarith said he has field capacity data. He will process it. (5) Needs cellorg: real(r8), pointer :: cellorg(:,:) ! column 3D org (kg/m3 organic matter) (nlevgrnd) Looks like it's used in lake CH4 calculation. Turn off lake CH4 for now and worry about it later. 2. Subroutine CNAllocation in CNAllocationMod: (1) lgsf (long growing season factor [0-1]) is declared but not used. (2) Needs gddmaturity: real(r8), pointer :: gddmaturity(:)! gdd needed to harvest gddmaturity is calculated from hybgdd (which is pft-dependent) in subroutine CropPhenology in CNPhenologyMod. Added hybgdd to CN_DriverMod. (3) leafout: 235: real(r8), pointer :: leafout(:) ! =gdd from top soil layer temperature 448: leafout => pps%gddtsoi gddtsoi is calculated in accFldsMod.F90. (4) Added graincn, fleafcn, fstemcn and ffrootcn to CN_DriverMod: real(r8), pointer :: graincn(:) ! grain C:N (gC/gN) real(r8), pointer :: fleafcn(:) ! leaf c:n during organ fill real(r8), pointer :: fstemcn(:) ! stem c:n during organ fill real(r8), pointer :: ffrootcn(:) ! froot c:n during organ fill graincn => pftcon%graincn fleafcn => pftcon%fleafcn ffrootcn => pftcon%ffrootcn fstemcn => pftcon%fstemcn 3. Subroutine CNPhenologyClimate called from subroutine Phenology: (1) gdd020, gdd820, gdd1020: calculated within subroutine CNPhenologyClimate. real(r8), pointer :: gdd020(:) ! 20-yr mean of gdd0 (ddays) real(r8), pointer :: gdd820(:) ! 20-yr mean of gdd8 (ddays) real(r8), pointer :: gdd1020(:) ! 20-yr mean of gdd10 (ddays) ========== 20170316: 1. Added to CN_DriverMod: pftcon%fd_pft = fd_pfx ! Fire duration (hr) pftcon%fsr_pft = fsr_pfx ! Fire spread rate (m/s) pftcon%grperc = grperx ! Growth respiration factor (unitless) pftcon%grpnow = grpnox ! Growth respiration factor (unitless) pftcon%gddmin = gddmix ! Minimim growing degree days used in CNPhenology pftcon%crop = crox ! Binary crop PFT flag (0: not crop; 1: crop) pftcon%mnNHplantdate = mnNHplantdatx ! Minimum planting date for the Northern Hemipsphere (YYYYMMDD) ! for CNAllocation pftcon%aleaff = aleafx ! Leaf Allocation coefficient parameter used in CNAllocation pftcon%allconsl = allconslx ! Leaf Allocation coefficient parameter power used in CNAllocation pftcon%allconss = allconssx ! Stem Allocation coefficient parameter power used in CNAllocation pftcon%arootf = arootfx ! Root Allocation coefficient parameter used in CNAllocation pftcon%arooti = arootix ! Root Allocation coefficient parameter used in CNAllocation pftcon%astemf = astemx ! Stem Allocation coefficient parameter used in CNAllocation pftcon%bfact = bfacx ! Exponential factor used in CNAllocation for fraction allocated to leaf pftcon%declfact = declfacx ! Decline factor for gddmaturity used in CNAllocation pftcon%fleafi = fleafx ! Leaf Allocation coefficient parameter fraction used in CNAllocation pftcon%pconv = pconx ! Deadstem proportions to send to conversion flux (fraction, 0-1), pconv+pprod10+pprod100 must sum to 1.0 pftcon%pprod10 = pprod10x ! Deadstem proportions to send to 10 year product pool (fraction, 0-1) pftcon%pprod100 = pprod100x ! Deadstem proportions to send to 100 year product pool (fraction, 0-1) pftcon%baset = basex ! Base Temperature, parameter used in accFlds (degree C) pftcon%mxtmp = mxtmx ! Max Temperature, parameter used in accFlds (degree C) 2. Changed "use pft2colMod , only: p2c" to "use subgridAveMod, only: p2c" in these modules: CNFireMod.F90 CNNDynamicsMod.F90 CNAnnualUpdateMod.F90 CNDecompMod.F90 (actually p2c is not being used. removed) ch4Mod.F90 CNSummaryMod.F90 CNVerticalProfileMod.F90 CNAllocationMod.F90 ========== 20170317: 1. Change "use shr_kind_mod , only: r8 => shr_kind_r8" to "use MAPL_ConstantsMod, ONLY: r8 => MAPL_R4" in these modules: CNFireMod.F90 CNGapMortalityMod.F90 CNGRespMod.F90 CNMRespMod.F90 CNNDynamicsMod.F90 CNNitrifDenitrifMod.F90 CNEcosystemDynMod.F90 CNCStateUpdate1Mod.F90 CNCStateUpdate2Mod.F90 CNCStateUpdate3Mod.F90 CNDecompCascadeMod_BGC.F90 CNDecompCascadeMod_CENTURY.F90 CNAnnualUpdateMod.F90 CNBalanceCheckMod.F90 CNDecompMod.F90 CNPhenologyMod.F90 clm_time_manager.F90 clm_varcon.F90 clm_varctl.F90 clmtypeInitMod.F90 CropRestMod.F90 CNNStateUpdate1Mod.F90 CNNStateUpdate2Mod.F90 CNSummaryMod.F90 clm_atmlnd.F90 CNNStateUpdate3Mod.F90 CNSoilLittVertTranspMod.F90 CNVegStructUpdateMod.F90 CNVerticalProfileMod.F90 CNWoodProductsMod.F90 CNAllocationMod.F90 pftdynMod.F90 surfrdMod.F90 perf_mod.F90 accFldsMod.F90 accumulMod.F90 CNPrecisionControlMod.F90 CNSetValueMod.F90 MAPL_ConstantsMod is in /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/src/GMAO_Shared/MAPL_Base/MAPL_Constants.F90. Our CLM4 version of CN model does not use CNPrecisionControlMod.F90. shr_strdata_mod.F90 is used in CNFireMod.F90 to get human density and lightning forcing. Need to think about a different way to read in these forcing data sets, so won't use shr_strdata_mod.F90. In CNFireMod.F90: use shr_kind_mod , only: r8 => shr_kind_r8, CL => shr_kind_CL In subroutine hdm_init: character(len=CL) :: stream_fldFileName_popdens ! population density streams filename character(len=CL) :: popdensmapalgo = 'bilinear' ! mapping alogrithm for population density In subroutine lnfm_init: character(len=CL) :: stream_fldFileName_lightng ! lightning stream filename to read character(len=CL) :: lightngmapalgo = 'bilinear'! Mapping alogrithm What's shr_kind_CL? In /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/csm_share/shr/shr_kind_mod.F90: integer,parameter :: SHR_KIND_CL = 256 ! long char Don't worry about it now since we will use a different way to read in the data. CNrestMod.F90 is what CLM4.5 uses to read in the restart file. We don't need it. 2. Deleted "use abortutils , only: endrun" and replaced "call endrun()" with "stop" in these routines: CNFireMod.F90 CNNitrifDenitrifMod.F90 CNCStateUpdate1Mod.F90 CNCStateUpdate2Mod.F90 CNCStateUpdate3Mod.F90 CNDecompCascadeMod_CENTURY.F90 CNBalanceCheckMod.F90 CNPhenologyMod.F90 clmtypeInitMod.F90 CropRestMod.F90 CNSummaryMod.F90 clm_atmlnd.F90 CNSoilLittVertTranspMod.F90 CNVerticalProfileMod.F90 CNAllocationMod.F90 pftdynMod.F90 surfrdMod.F90 accFldsMod.F90 accumulMod.F90 CNPrecisionControlMod.F90 /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/atm/cam/src/utils/abortutils.F90 clm_time_manager.F90: use the one in our CLM4 version of CN model. Will modify later if needed. pftdynMod.F90, surfrdMod.F90: too many "call endrun", how to replace in a fast way? pftdynMod.F90: CNEcosystemDynMod.F90 115: use pftdynMod , only: CNHarvest clm_initializeMod.F90 285: use pftdynMod , only : pftdyn_init, pftdyn_interp 291: use pftdynMod , only : pftwt_init Looks like only these few subroutines in pftdynMod.F90 are used. Stripped out the rest subroutines and then replaced "call endrun" with "stop". surfrdMod.F90: Looks like only subroutine crop_prog is used. Subroutines surfrd_get_globmask, surfrd_get_grid, surfrd_get_topo, and surfrd_get_data are used in clm_initializeMod.F90 but we may get rid of clm_initializeMod.F90. These subroutines are for reading global land mask, grid/ladnfrac data, grid topography and surface dataset. We don't need these. Removed them. ! !PUBLIC DATA MEMBERS: logical, public :: crop_prog = .false. ! If prognostic crops is turned on Moved this to clm_varctl.F90. Removed surfrdMod.F90. ========== 20170320: 1. Went through these routines to see what changes Greg made. Use this as a guide to change these routines in the CLM4.5 version of our CN model. clm_time_manager.F90 clm_varcon.F90 clm_varpar.F90 CNAllocationMod.F90 CNAnnualUpdateMod.F90 CNBalanceCheckMod.F90 CNCStateUpdate1Mod.F90 CNCStateUpdate2Mod.F90 CNCStateUpdate3Mod.F90 CNDecompMod.F90 CNEcosystemDynMod.F90 CNFireMod.F90 CNGapMortalityMod.F90 CNGRespMod.F90 CNiniTimeVar.F90 CNMRespMod.F90 CNNDynamicsMod.F90 CNNStateUpdate1Mod.F90 CNNStateUpdate2Mod.F90 CNNStateUpdate3Mod.F90 CNPhenologyMod.F90.clm4 CNPrecisionControlMod.F90 CNSetValueMod.F90 CNSummaryMod.F90 CNVegStructUpdateMod.F90 CNWoodProductsMod.F90 pftvarcon.F90 clm_initializeMod.F90 clmtype.F90 clmtypeInitMod.F90 accFldsMod.F90 accumulMod.F90 (1) clm_time_manager.F90: Starting from the clm_time_manager.F90 in our CLM4 version of CN. Need to add these subroutines (from /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/util_share/clm_time_manager.F90): get_ref_date, get_curr_date, get_curr_time: for CNPhenologyMod.F90, CNFireMod.F90, CNDVMod.F90 and others get_calendar: for CNFireMod.F90 get_curr_calday, get_calday: for CNPhenologyMod.F90 is_restart, get_driver_start_ymd, get_start_date: for CropRestMod.F90 is_end_curr_day: for accFldsMod.F90 (1.1) Modified function get_nstep. (1.2) For get_rad_step_size, change "get_rad_step_size" to "get_rad_step_size => get_step_size" in CNVegStructUpdateMod.F90. (1.3) For subroutine get_curr_date: Added subroutine get_CurrDateTime to CN_DriverMod.F90; Added subroutine get_curr_date to clm_time_manager.F90. 2. Added num_pcropp and filter_pcropp (used in CNEcosystemDyn and maybe more) to CN_DriverMod. ========== 20170321: 1. Continued 1 on 20170320. (1) Modified clm_time_manager.F90: (1.4) Skip get_start_date for now. It's only called from subroutine checkDates (in CropRestMod), and checkDates is only called from subroutine CropRest which will be skipped at least at the early testings. (1.5) Skip get_driver_start_ymd for now, for the same reason as in 1.4. (1.6) Skip get_ref_date for now. It's only called from CNDVMod.F90. (1.7) Same as 1.6, skip get_curr_time for now. (1.8) Added function get_curr_calday. (1.9) Skip get_calday for now. It's only called from subroutine CropPhenologyInit. Since crop_prog = .false. in CLM4.5 surfrdMod.F90 (and now in clm_varctl.F90 in Catchment-CN), CropPhenologyInit is not called at all, so no need to worry about get_calday unless crop_prog is turned on. (1.10) Skip get_calendar. It's only used in subroutines hdm_init and lnfm_init in CNFireMod. We will read in hdm and lnfm in a different way, so these two subroutines will not be used. (1.11) Added is_end_curr_day and is_restart. (2) Modified clm_varcon.F90 (3) Modified clm_varpar.F90 (4) Modified CNAllocationMod.F90 (5) CNAnnualUpdateMod.F90: no change made (6) Modified CNBalanceCheckMod.F90 (7) CNCStateUpdate1Mod.F90, CNCStateUpdate2Mod.F90 and CNCStateUpdate3Mod.F90: no change made (8) CNDecompMod.F90: no change made /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/tools/clm4_5/mksurfdata_map/src/nanMod.F90 and /discover/nobackup/fzeng/clm_orig/ccsm4_0/models/lnd/clm/src/main/nanMod.F90 are identical. cp -p clm4/nanMod.F90.clm4 nanMod.F90 ========== 20170322: 1. Continued 1 on 20170321. (9) Modified CNEcosystemDynMod.F90. (10) CNFireMod.F90: no change (11) CNGapMortalityMod.F90: no change (12) CNGRespMod.F90: no change (13) Copied and modified CNiniTimeVar.F90 cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/main/CNiniTimeVar.F90 . (14) CNMRespMod.F90: no change (15) Modified CNNDynamicsMod.F90 (16) CNNStateUpdate1Mod.F90, CNNStateUpdate2Mod.F90, CNNStateUpdate3Mod.F90: no change (17) Modified CNPhenologyMod.F90 (18) CNPrecisionControlMod.F90: no change (19) Modified CNSetValueMod.F90 (20) Modified CNSummaryMod.F90. It seems that the subroutine CAdjust Greg added is not called in the CLM4 version of our CN model, so this subroutine is not added to this new CNSummaryMod.F90. (21) Modified CNVegStructUpdateMod.F90 (22) Modified CNWoodProductsMod.F90 (23) pftvarcon.F90: no change 2. Organized the routines. rm DryDepVelocity.F90 rm DUSTMod.F90 rm CNrestMod.F90 rm clm_initializeMod.F90 mkdir ch4-isotope-CNDV mv ch4Mod.F90 ch4RestMod.F90 ch4varcon.F90 CNC14DecayMod.F90 CNCIsoFluxMod.F90 CNDVEcosystemDynIniMod.F90 CNDVEstablishmentMod.F90 CNDVLightMod.F90 CNDVMod.F90 initch4Mod.F90 ch4-isotope-CNDV 3. Subroutines in STATICEcosysDynMod.F90 are only used in clm_initializeMod.F90. Since clm_initializeMod.F90 is not needed, deleted STATICEcosysDynMod.F90. rm STATICEcosysDynMod.F90 Subroutines in shr_strdata_mod.F90 are only used in CNFireMod.F90 to read in hdm and lnfm data. We will read these data sets in a different way, so deleted it. rm shr_strdata_mod.F90 ========== 20170323: 1. Keep CropRestMod.F90 for future. m_dbg_routines.F90 is from Catchment. Keep it. Subroutines in MEGANFactorsMod.F90 are only called from VOCEmissionMod.F90. Removed it. rm MEGANFactorsMod.F90 rm VOCEmissionMod.F90 2. Checked the CN routines if the ones listed in "use" are actually used and exist or not. (1) clm_time_manager.F90 (2) clm_varcon.F90 Need SHR_CONST_PDB. Found it in /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/csm_share/shr/shr_const_mod.F90. Added it to our shr_const_mod.F90. (3) clm_varpar.F90 ========== 20170324: 1. For this block in CNFireMod.F90: !assign local pointers to derived type members (grid-level) forc_rh => clm_a2l%forc_rh forc_wind => clm_a2l%forc_wind forc_t => clm_a2l%forc_t forc_rain => clm_a2l%forc_rain forc_snow => clm_a2l%forc_snow latdeg => grc%latdeg ??do we have relative humidity from Catchment?? In compute_rc.F90: real(r8) :: ceair ! vapor pressure of air, constrained (Pa) real(r8), intent(in) :: eair(nch) ! vapor pressure of canopy air (Pa) real(r8), intent(in) :: esat_tv(nch) ! saturation vapor pressure at t_veg (Pa) real ei(nch) ! vapor pressure inside leaf (sat vapor press at tc) (Pa) real ea(nch) ! vapor pressure of canopy air (Pa) real(r8), intent(in) :: rh_can ! canopy air realtive humidity ei(n) = MAPL_EQsat(tc(n),DQ=deldT(n)) ea(:) = pbot(:) * qa(:) / (0.622 + qa(:)) ! canopy air vapor pressure (Pa) eair=ei esat_tv=ea ceair = min( eair(p), esat_tv(p) ) rh_can = ceair / esat_tv(p) MAPL_EQsat is in src/GMAO_Shared/MAPL_Base/MAPL_SatVapor.F90. In clsm_ensdrv_force_routines.F90, saturated specific humidity is calculated as: Qair_sat = MAPL_EQsat(met_force(i)%Tair, met_force(i)%Psurf ) Since clm_a2l%forc_rh is only used in CNFireMod, maybe I can do the calculation in CNFireMod? No, because dtcn does not equal to dt. Changes made: (1) process_cn: Added L291-292: real :: Qair_sat ! saturated specific humidity (kg/kg) real, dimension(N_cat) :: Qair_relative ! relative humidity (%) Added L1377-1383: ! compute relative humidity (%) used in CNFireMod ! ----------------------------------------------- do n = 1,N_cat Qair_sat = MAPL_EQsat(met_force(n)%Tair, met_force(n)%Psurf ) Qair_relative(n) = met_force(n)%Qair / Qair_sat * 100. end do Qair_relative(:) = min(max(0., Qair_relative(:)), 100.) Added L1397-1400: cn_progn(:)%rhm = cn_progn(:)%rhm + Qair_relative(:) cn_progn(:)%windm = cn_progn(:)%windm + met_force(:)%Wind cn_progn(:)%rainfm = cn_progn(:)%rainfm + met_force(:)%Rainf cn_progn(:)%snowfm = cn_progn(:)%snowfm + met_force(:)%Snowf Added L1448-1451: cn_progn%rhm = cn_progn%rhm * scale cn_progn%windm = cn_progn%windm * scale cn_progn%rainfm = cn_progn%rainfm * scale cn_progn%snowfm = cn_progn%snowfm * scale Changed L1479-1480 to pass these fields to CN_Driver Added L1503-1506: cn_progn%rhm = 0. cn_progn%windm = 0. cn_progn%rainfm = 0. cn_progn%snowfm = 0. (2) CN_DriverMod: L47: modified arguments Added L68-71: real*4, dimension(nch), intent(in) :: rhm ! relative humidity (%) real*4, dimension(nch), intent(in) :: windm ! wind speed (m/s) real*4, dimension(nch), intent(in) :: rainfm ! rainfall (convective + largescale) (kg/m2/s) real*4, dimension(nch), intent(in) :: snowfm ! snowfall (kg/m2/s) Added L153-157: real, pointer :: forc_t(:) !air temperature (K) real, pointer :: forc_rh(:) !relative humidity (%) real, pointer :: forc_wind(:) !wind speed (m/s) real, pointer :: forc_rain(:) !rainfall (convective + largescale) (kg/m2/s) real, pointer :: forc_snow(:) !snowfall (kg/m2/s) Added L437-441: forc_rh => grc%forc_rh forc_wind => grc%forc_wind forc_t => grc%forc_t forc_rain => grc%forc_rain forc_snow => grc%forc_snow Added L620-626: ! forcing for CNFireMod ! --------------------- forc_rh(:) = rhm(:) forc_wind(:) = windm(:) forc_t(:) = tm(:) forc_rain(:) = rainfm(:) forc_snow(:) = snowfm(:) (3) CNFireMod: Changed this block !assign local pointers to derived type members (grid-level) forc_rh => clm_a2l%forc_rh forc_wind => clm_a2l%forc_wind forc_t => clm_a2l%forc_t forc_rain => clm_a2l%forc_rain forc_snow => clm_a2l%forc_snow latdeg => grc%latdeg to !assign local pointers to derived type members (grid-level) forc_rh => grc%forc_rh forc_wind => grc%forc_wind forc_t => grc%forc_t forc_rain => grc%forc_rain forc_snow => grc%forc_snow latdeg => grc%latdeg 2. Greg added this under "define the gridcell structure" in clmtype.F90: real, pointer :: forc_ndep(:) ! nitrogen deposition rate (gN/m2/s) Following what he did, added these to "define the gridcell structure" in clmtype.F90 (for the upgraded code): real(r8), pointer :: forc_ndep(:) ! nitrogen deposition rate (gN/m2/s) real(r8), pointer :: forc_rh(:) ! relative humidity (%) real(r8), pointer :: forc_wind(:) ! wind speed (m/s) real(r8), pointer :: forc_t(:) ! air temperature (K) real(r8), pointer :: forc_rain(:) ! rainfall (convective + largescale) (mm/s) real(r8), pointer :: forc_snow(:) ! snowfall (mm/s) ========== 20170327: 1. Checked the CN routines if the ones listed in "use" are actually used and exist or not. Deleted if not used. Added if not exist. (1) CNAllocationMod.F90 (2) CNAnnualUpdateMod.F90 (3) CNBalanceCheckMod.F90 (4) CNCStateUpdate1Mod.F90 (5) CNCStateUpdate2Mod.F90 (6) CNCStateUpdate3Mod.F90 (7) CNDecompCascadeMod_BGC.F90 (8) CNDecompCascadeMod_CENTURY.F90 (9) CNDecompMod.F90 (10) CN_DriverMod.F90 (11) CNEcosystemDynMod.F90 (12) CNFireMod.F90 (13) CNGapMortalityMod.F90 (14) CNGRespMod.F90 (15) CNiniTimeVar.F90 (16) CNMRespMod.F90 (17) CNNDynamicsMod.F90 (18) CNNStateUpdate1Mod.F90 (19) CNNStateUpdate2Mod.F90 (20) CNNStateUpdate3Mod.F90 (21) CNPhenologyMod.F90 (22) CNPrecisionControlMod.F90 (23) CNSetValueMod.F90 CNFireMod: commented out anything used in hdm_init, hdm_interp, lnfm_init and lnfm_interp commented out "call CNFireInit" in CNEcosystemDynMod.F90. How to avoid using this "use shr_infnan_mod , only: shr_infnan_isnan"? Asked Sarith for help. Sarith suggests to try (see his email on 3/28/2017): use, intrinsic :: ieee_arithmetic, only: & isnan => ieee_is_nan ========== 20170328: 1. Checked the CN routines if the ones listed in "use" are actually used and exist or not. Deleted if not used. Added if not exist. (24) CNSoilLittVertTranspMod.F90 Needs TridiagonalMod, copied it: cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeophys/TridiagonalMod.F90 . Modified TridiagonalMod.F90. Added this below (needed in TridiagonalMod.F90) to clm_varcon.F90: ! urban column types integer, parameter :: icol_roof = 61 integer, parameter :: icol_sunwall = 62 integer, parameter :: icol_shadewall = 63 integer, parameter :: icol_road_imperv = 64 integer, parameter :: icol_road_perv = 65 (25) CNSummaryMod.F90 (26) CNVegStructUpdateMod.F90 (27) CNVerticalProfileMod.F90 (28) CNWoodProductsMod.F90 (29) compute_rc.F90 (30) CropRestMod.F90 In CLM4.5 subroutine CropRest is called from main/restFileMod.F90. It's currrently not called in Catchment-CN, so don't worry about it yet. (31) nanMod.F90 (32) pftdynMod.F90 "use clm_varsur , only : pctspec", so copied clm_varsur.F90 and modified it (in case it's needed to get the code compile): cp -p /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/util_share/clm_varsur.F90 . pftdyn_init and pftdyn_interp are not called in Catchment-CN, so skipped them now. pftwt_init and pftwt_interp are for CNDV. Skipped them now. CNHarvest is actually not called at the time. So no need to worry about this now. (see notes on 20170228) Commented out "use spmdMod" and "use ncdio_pio" that we don't have now. (33) pftvarcon.F90 (34) shr_const_mod.F90 (35) subgridAveMod.F90 2. Anything to remove from these files below? (1) rm clm_atmlnd.F90 (2) accFldsMod.F90: don't remove anything for now. For subroutine updateAccFlds forc_t => clm_a2l%forc_t forc_rain => clm_a2l%forc_rain forc_snow => clm_a2l%forc_snow forc_solad => clm_a2l%forc_solad ! (heald 04/06) forc_solai => clm_a2l%forc_solai ! (heald 04/06) Need to change clm_a2l to grc. Wanted to add forc_solad and forc_solai to CN_DriverMod. real, pointer :: forc_solad(:,:) !direct beam radiation (numrad), (visible only) real, pointer :: forc_solai(:,:) !diffuse radiation (numrad), (visible only) However, it looks like Catchment doesn't have direct beam radiation and diffuse radiation. It has PARdir and PARdif, but they are not the same as solad and solai. Variables of which calculation needs forc_solad or forc_solai include: FSD240, FSD24, FSI24 and FSI240 (in accFldsMod.F90) These variable are for VOC emissions which we will skip at least for now. (Checked and made sure that these variables are not used in our CN routines.) Therefore, don't add them to CN_DriverMod yet and comment the lines with these two forcing varialbes and the four variables FSD240, FSD24, FSI24 and FSI240 in accFldsMod.F90. (3) accumulMod.F90 Do we need /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/util_share/decompMod.F90? The subroutines get_proc_bounds and get_proc_global are used in accFldsMod.F90, accumulMod.F90 and a few others. How to pass numg, numl, numc and nump from subroutine CN_Driver to subroutine init_accum_field in accumulMod.F90? Also, where are these cumulative fields computed in CLM4.5, before the CN routines? (see CLM4 vs CLM4.5.ppt on my local machine) Added "call initAccFlds()" to process_cn.F90. Added an input argument "nch - number of tiles" to subroutine initAccFlds in accFldsMod.F90 and subroutine init_accum_field in accumulMod.F90. 'TREFAV', t_ref2m: needed in CNPhenologyMod.F90 'TREFAV_U', t_ref2m_u: not used in our CN model, skip it in accFldsMod.F90 'TREFAV_R', t_ref2m_r: not used in our CN model, skip it in accFldsMod.F90 'T_VEG24', t_veg24: not used in our CN model, skip it in accFldsMod.F90 'T_VEG240', t_veg240: not used in our CN model, skip it in accFldsMod.F90 'FSUN24', fsun24: not used in our CN model, skip it in accFldsMod.F90 'FSUN240', fsun240: not used in our CN model, skip it in accFldsMod.F90 'LAIP', elai_p: not used in our CN model, skip it in accFldsMod.F90 Decided not to use their accumulMod.F90 and accFldsMod.F90, and just add the calculation to process_cn.F90 like what I did for t10. (4) clmtype.F90 (5) clmtypeInitMod.F90 Modified subroutine initClmtype to avoid using get_proc_bounds and get_proc_global. clm_varctl.F90 perf_mod.F90 ========== 20170329: 1. Continued working on 2 on 20170328. ========== 20170330: 1. Continued working on 2 on 20170328. ========== 20170331: 1. Added calculation of prec10, prec60, TREFAV and the related t_ref2m_min and t_ref2m_max to process_cn.F90 and CN_DriverMod. ========== 20170403: 1. The prognostic crop model: According to p347 in CLM4.5 tech note, if prognostic crop is on, the managed crop types, corn, soybean and temperate cereals will be included, in addition to the c3_crop which is unmanaged and in the default list of PFTs. To include these managed crop types, I need to do these: (1) Add the variables in crop restart file (a separate restart file in CLM4.5) to our catchcn_internal_rst (2) Add the computation of the accumulated variables needed by the crop model (3) Change the way of reading crop rest info Can I leave prog_crop off for now? Since Sarith has included these managed crop types in the new FPT maps, I might have to just turn on crop_prog. Went through CropRestMod.F90 and made a list of variables that need to be added to the restart file, and guessed their initial values for cold start from model code. Added calculation of TDM10 and TDM5 to process_cn.F90. 2. CNiniTimeVar.F90: initialize the variables in the restart file if a restart file does not exist. ========== 20170404: 1. Added calculation of GDD0, GDD8, GDD10, GDDPLANT and GDDTSOI to process_cn.F90 and CN_DriverMod. 2. Compared CNiniTimeVar.F90 between CLM4 and CLM4.5 to see what need to be added to the restart file!! ========== 20170406: 1. Summarized what to be added or removed from the current catchcn_internal_rst (see CLM4 vs CLM4.5.pptx in my local machine). ========== 20170407: 1. CLM4 also has clm_initializeMod.F90, but we didn't use it. How do we avoid using it now? Greg avoided using clm_initializeMod.F90 by calling subroutine initClmtype (in clmtypeInitMod.F90) from subroutine CN_init in CN_DriverMod. Do the same. Anything else? 2. Identified and removed those not needed in clmtypeInitMod.F90. 3. CNPhenologyMod.F90 also needs t10, so passed t10 to CN_Driver in process_cn.F90. ========== 20170410: 1. Updated the CLM types in CN_DriverMod.F90. ========== 20170411: 1. Updated the CLM types in CN_DriverMod.F90 and subgridAveMod.F90. ========== 20170412: 1. The CLM4.5 CN routines do not have cwdc, litr1c, litr2c, litr3c, soil1c, soil2c, soil3c, soil4c, litr1n, litr2n, litr3n, soil1n, soil2n, soil3n, or soil4n any more. See what replaced them! In clm_varpar.F90: #ifndef CENTURY_DECOMP ! parameters for decomposition cascade integer, parameter :: ndecomp_pools = 8 integer, parameter :: ndecomp_cascade_transitions = 9 integer, parameter :: i_met_lit = 1 integer, parameter :: i_cel_lit = 2 integer, parameter :: i_lig_lit = 3 integer, parameter :: i_cwd = 4 integer, parameter :: nsompools = 4 #else ! parameters for decomposition cascade integer, parameter :: ndecomp_pools = 7 integer, parameter :: ndecomp_cascade_transitions = 10 integer, parameter :: i_met_lit = 1 integer, parameter :: i_cel_lit = 2 integer, parameter :: i_lig_lit = 3 integer, parameter :: i_cwd = 4 integer, parameter :: nsompools = 3 #endif So if we choose NOT to use CENTURY_DECOMP, there will be 8 decomposition carbon pools, like what we have in our CLM4 version of Catchment-CN. This is confirmed by CLM4.5 Tech Note (p276): "CLM allows the user to define, at compile time, between 2 contrasting hypotheses of decomposition as embodied by two separate decomposition submodels: the CLM-CN pool structure used in CLM4.0, or a second pool structure, characterized by slower decomposition rates, based on the Century model (Parton et al. 1988). In addition, the user can choose, at compile time, whether to allow nlev to equal 1, as in CLM4.0, or to equal the number of soil levels used for the soil hydrology (default 10)." In CNCStateUpdate2Mod.F90: ! column level carbon fluxes from gap-phase mortality do j = 1,nlevdecomp ! column loop do fc = 1,num_soilc c = filter_soilc(fc) ! column gap mortality fluxes decomp_cpools_vr(c,j,i_met_lit) = decomp_cpools_vr(c,j,i_met_lit) + gap_mortality_c_to_litr_met_c(c,j) * dt decomp_cpools_vr(c,j,i_cel_lit) = decomp_cpools_vr(c,j,i_cel_lit) + gap_mortality_c_to_litr_cel_c(c,j) * dt decomp_cpools_vr(c,j,i_lig_lit) = decomp_cpools_vr(c,j,i_lig_lit) + gap_mortality_c_to_litr_lig_c(c,j) * dt decomp_cpools_vr(c,j,i_cwd) = decomp_cpools_vr(c,j,i_cwd) + gap_mortality_c_to_cwdc(c,j) * dt end do end do So decomp_cpools_vr (and therefore decomp_cpools which is the sum of decomp_cpools_vr of different soil layers, in our case only one soil layer so nlevdecomp = 1) is: decomp_cpools_vr(:,1,1): metabolic (equivalant to "Labile" in our CLM4 version of Catchment-CN), = litr1c decomp_cpools_vr(:,1,2): cellulose, = litr2c decomp_cpools_vr(:,1,3): lignin, = litr3c decomp_cpools_vr(:,1,4): cwd What about the order of the soil OM pools? In CNDecompCascadeMod_BGC.F90: i_soil1 = 5 i_soil2 = 6 i_soil3 = 7 i_soil4 = 8 Therefore: decomp_cpools_vr(:,1,5): soil1c decomp_cpools_vr(:,1,6): soil2c decomp_cpools_vr(:,1,7): soil3c decomp_cpools_vr(:,1,8): soil4c 2. Why are totlitn and totsomn not in the rst? Because they can be computed by summing up all the litter or SOM pools that are already in the rst. 3. Modified subroutines CN_init and CN_exit: listed the variables needed in the new offline restart file Referred to the files below: (1) /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CNrestMod.F90 (2) /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeochem/CropRestMod.F90 (3) /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeophys/BiogeophysRestMod.F90 (4) /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/util_share/accumulMod.F90 and /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/main/accFldsMod.F90 alt and alt_indx are not needed in the rst. altmax, altmax_lastyear, altmax_indx, and altmax_lastyear_indx are calculated in /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeophys/ActiveLayerMod.F90. ! on a set annual timestep, update annual maxima ! make this 1 January for NH columns, 1 July for SH columns call get_curr_date(year, mon, day, sec) dtime = get_step_size() if ( (mon .eq. 1) .and. (day .eq. 1) .and. ( sec / dtime .eq. 1) ) then do fc = 1,num_soilc c = filter_soilc(fc) g = cgridcell(c) if ( lat(g) .gt. 0. ) then altmax_lastyear(c) = altmax(c) altmax_lastyear_indx(c) = altmax_indx(c) altmax(c) = 0. altmax_indx(c) = 0 endif end do endif if ( (mon .eq. 7) .and. (day .eq. 1) .and. ( sec / dtime .eq. 1) ) then do fc = 1,num_soilc c = filter_soilc(fc) g = cgridcell(c) if ( lat(g) .le. 0. ) then altmax_lastyear(c) = altmax(c) altmax_lastyear_indx(c) = altmax_indx(c) altmax(c) = 0. altmax_indx(c) = 0 endif end do endif Added this calculation to CN_DriverMod. 4. Some constants: clm_varcon.F90 49: real(r8), public, parameter :: spval = 1.e36_r8 ! special value for real data 50: integer , public, parameter :: ispval = -9999 ! special value for int data ========== 20170413: 1. Continue working on 3 above. t10, prec10, prec60, trefav, t_ref2m_max_inst, t_ref2m_min_inst, t_ref2m_max, t_ref2m_min, tdm5, tdm10, gdd0, gdd8, and gdd10 are not in either CNrestMod.F90 or CropRestMod.F90. Are they unneeded in the rst file? Yes. They need to be in the rst file. They are in /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/util_share/accumulMod.F90 and /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/clm4_5/main/accFldsMod.F90. Also, istep not needed in the rst file? Avoid using it by doing a run longer than 60 days at the very beginning. Editing the variables in cncol and cnpft in CN_init and CN_exit: (1) column-level: Remove these: 30: ann_farea_burned 35: fire_prob 36: fireseasonl 39: me 40: mean_fire_prob Added these: altmax altmax_lastyear altmax_indx altmax_lastyear_indx (2) pft-level: Added these: gddplant gddtsoi peaklai idop aleaf aleafi astem astemi htmx hdidx vf cumvd croplive cropplant harvdate gdd1020 gdd820 gdd020 gddmaturity huileaf huigrain grainc grainc_storage grainc_xfer grainn grainn_storage grainn_xfer fert_counter fert grain_flag 2. Changed VAR_COL and VAR_PFT in clm_varpar.F90. 3. For the grid-level variables to be added to the rst file: modified catch_types.F90 and clsm_ensdrv_init_routines.F90 following what I did for T10. 4. Asked Sarith about how to set the #if directives when compiling the code. Can do this: gmake install FOPT='-DXXX' Example, for "#ifdef CN", if I want to use and compile CN, then "XXX" above is "CN" (i.e. the part after "#ifdef") and I should do this: gmake install FOPT='-DCN' For "#ifndef CENTURY_DECOMP" and if I don't want to use CENTURY_DECOMP, just don't include it and simply do "gmake install". Same thing for "#if (defined NOFIRE)" and if I want FIRE ON, just don't do anything about it and simply do "gmake install". A better way is to modify GNUmakefile_openmp to make this permanent, take "#ifdef CN" and if I want to use CN as an example: Change this line: USER_FFLAGS = $(BIG_ENDIAN) $(OPENMPFLAG) \ to: USER_FFLAGS = $(BIG_ENDIAN) $(OPENMPFLAG) $(D)CN \ Overall, this is not a good practice. Sarith suggests me to just delete "#ifdef CN" in the CN routines. 5. Also asked Sarith if the order of the files in GNUmakefile_openmp important. He said NO, which is good for me! ========== 20170414: 1. Continued working on catch_types.F90 in 3 on 20170413 above. Updated the description of variables in cncol and cnpft in catch_types.F90. 2. Before compiling the code, commented out "#ifdef CN" in all CN routines (see 4 on 20170413): ========== 20170417: 1. Continued working on 2 on 20170414. 2. GNUmakefile_openmp: Added: clm_varctl.F90 CNDecompCascadeMod_BGC.F90 CNSoilLittVertTranspMod.F90 CNVerticalProfileMod.F90 TridiagonalMod.F90 Commented out GEOS_CatchCNGridComp.F90 for now. 3. masterproc In CNFireMod and clm_varctl.F90, "masterproc" is in the subroutines that we are not going to use, so just leave it there. 4. Tried to compile, fixing bugs (1) For subroutine get_curr_date in clm_time_manager.F90: Can't do "use CN_DriverMod, only : curr_date_time" Can't do "use date_time_util, only: date_time_type" Can't do "use date_time_util, only: date_time_type" in GEOSsurface_GridComp/GEOSland_GridComp/Shared/update_model_paras.F90 See current solution in subroutine upd_curr_date_time in update_model_paras.F90. (2) Declared, allocated and initialized these in clmtype.F90 and clmtypeInitMod.F90: pftcon%declfact pftcon%bfact pftcon%aleaff pftcon%arootf pftcon%astemf pftcon%arooti pftcon%fleafi pftcon%allconsl pftcon%allconss pftcon%grperc ========== 20170418: 1. Continued working on 4 on 20170417 above. (1) Declared, allocated and initialized these in clmtype.F90 and clmtypeInitMod.F90: pftcon%fsr_pft pftcon%fd_pft pftcon%gddmin pftcon%hybgdd pftcon%lfemerg pftcon%grnfill pftcon%mxmat pftcon%minplanttemp pftcon%planttemp cps%psiwilt (2) CNPhenologyMod: Temporarily fixed the nyrs (number of crop spinup years) to make it compile. Added mnNHplantdate, mnSHplantdate, mxNHplantdate, and mxSHplantdate to CN_DriverMod, clmtype.F90 and clmtypeInitMod.F90. Added function get_calday to clm_time_manager.F90. (3) Commented out all "call t_startf" and "call t_stopf" in CNEcosystemDynMod.F90. They are performance timers. (4) Added CNSoilLittVertTranspMod.F90 to GNUmakefile_openmp. Removed CNHarvestMod.F90 from GNUmakefile_openmp. ========== 20170419: 1. Continued fixing the bugs found from compiling the offline code: (1) CN_DriverMod: Can't use date_time_util here. Call get_curr_date to get the date_time info. Commented out pftcon%crop, ftcon%pconv, pftcon%pprod10, and pftcon%pprod100, which are used. (2) Declared, allocated and initialized these in clmtype.F90 and clmtypeInitMod.F90: cps%psisat pftcon%mxtmp pftcon%baset (3) qe25 is needed to compute sun induced fluorescence (SIF) in compute_rc.F90 but it's no longer available in CLM4.5. What should I do? (3.1) For j (electron transport, umol co2/m**2/s): According to CLM4 Tech Note (p163-165): j in CLM4 version of Catchment-CN = 4.6 * APAR * quantum efficiency where 4.6: 4.6 umol photons per Joule APAR: absorbed photosynthetically active radiation (W m-2) quantum efficiency: umol CO2 per umol photons (not necessary at 25 degree C as in the model code) According to CLM4.5 Tech Note (p186-187): je in CLM4.5 version of Catchment-CN is electron transport rate (umol electrons/m**2/s) To calculate j (umol co2/m**2/s) from je (umol electrons/m**2/s), it seems that quantum efficiency (umol CO2 per umol photons) is needed. photon = electron? However, in CLM4.5 subroutine Photosynthesis: if (c3flag(p)) then qe(p) = 0._r8 theta_cj(p) = 0.98_r8 bbbopt(p) = 10000._r8 mbbopt(p) = 9._r8 else qe(p) = 0.05_r8 theta_cj(p) = 0.80_r8 bbbopt(p) = 40000._r8 mbbopt(p) = 4._r8 end if Can't just do j = je * qe !!!!DON'T know how to calculate j!!!! (3.2) For je (actual electron transport) in CLM4 version and je_sif in CLM4.5 version: According to CLM4 version of Catchment-CN: je = max(psn*(ci+2*cp)/max(ci-cp,1.e-8),0.), for c3 plants = psn, for c4 plants where ci: internal leaf CO2 partial pressure (Pa) cp: CO2 compensation point (Pa) Temporary solution: In CN_init, give qe25 the same numbers (0.06 for c3 plants and 0.04 for c4 plants) as the CLM4 version of Catchment-CN just to get the code compiling. Need to ask Jung-Eun Lee for the CLM4.5 version of fluorescence calculation!! Was able to compile both GEOScatchCN_GridComp and GEOSlana_GridComp! ========== 20170420: 1. Organized the slides about what to add to the rst file, and sent them to Sarith. 2. Created /home/fzeng/Catchment/matlab/clm4-to-clm4.5/m_map/plot_bcs_fz.m to plot the variables in /gpfsm/dnb02/smahanam/bcs/Heracles-4_3/Heracles-4_3_MERRA-3/DE_00720x00360_PE_0720x0360/clsm/CLM4.5_abm_peatf_gdp_hdm_fc and compared them to the original maps made using Panoply. Learned from Eunjee how to display continental boundaries. 3. Helped Eunjee with questions about the manuscript. ========== 20170421: 1. Modified the code to avoid having istep and restyear in the restart file. See istep in process_cn.F90 and nyrs in CNPhenologyMod.F90. Met with Sarith and asked him to take a look at the modifications. Asked Sarith about the unit of our field capacity data. He said it's m3/m3. 2. Talked with Eunjee and Brad about the NBP in the manuscript. ========== 20170426: 1. /gpfsm/dnb02/smahanam/bcs/Heracles-4_3/Heracles-4_3_MERRA-3/DE_00720x00360_PE_0720x0360/clsm/README: Please see below Table 2 for CLM (CLM4.5) and CLM-Carbon (CLM4.5-Carbon) land cover classification =================================================================================== Land Cover CLM (CLM4.5) CLM-Carbon Map (CLM4.5-Carbon) Class Class Legend ----------------------------------------------------------------------------------- Bare 1 - BARE Needleleaf evergreen temperate tree 2 1 NLEt Needleleaf evergreen boreal tree 3 2 NLEB Needleleaf deciduous boreal tree 4 3 NLDB Broadleaf evergreen tropical tree 5 4 BLET Broadleaf evergreen temperate tree 6 5 BLEt Broadleaf deciduous tropical tree 7 6 BLDT Broadleaf deciduous temperate tree 8 7 BLDt Broadleaf deciduous boreal tree 9 8 BLDB Broadleaf evergreen temperate shrub 10 9 BLEtS Broadleaf deciduous temperate shrub 11 10 BLDtS Broadleaf deciduous temperate shrub[moisture stress only] - 11 BLDtSm Broadleaf deciduous boreal shrub 12 12 BLDBS Arctic c3 grass 13 13 AC3G Cool c3 grass 14 14 CC3G Cool c3 grass [moisture stress only] - 15 CC3Gm Warm c4 grass 15 16 WC4G Warm c4 grass [moisture stress only] - 17 WC4Gm Crop 16 18 CROP (-) Crop [moisture stress only] - 19 CROPm(-) (C3_crop) (16) (18) C3CROP (C3_irrigated) (17) (19) C3IRR (Corn) (18) (20) CORN (Irrigated corn) (19) (21) ICORN (Spring temperate cereal) (20) (22) STCER (Irrigated spring temperate cereal) (21) (23) ISTCER (winter temperate cereal) (22) (24) WTCER (Irrigated winter temperate cereal) (23) (25) IWTCER (Soybean) (24) (26) SOYB (Irrigated Soybean) (25) (27) ISOYB Water 17 - ----------------------------------------------------------------------------------- Table 2: CLM and CLM-Carbon land cover classification description. CLM-4.5 and CLM-4.5-Carbon types are in brackets. These two Crop 16 18 CROP (-) Crop [moisture stress only] - 19 CROPm(-) are not in the CLM4.5-Carbon class. They are replaced by (C3_crop) (16) (18) C3CROP (C3_irrigated) (17) (19) C3IRR Confirmed with Sarith on 20170427. ========== 20170428: 1. Modified the data table in CN_init to match the the table above (on 20170426) because Sarith includes split PFT types in the veg_type_fracs file. Updated numpft to 27 (from 24) in clm_varpar.F90. Updated pftvarcon.F90. Not sure if we need to do the same for the split types. e.g. In compute_rc.F90: 503: use pftvarcon , only : nbrdlf_dcd_tmp_shrub 799: if (ivt == nbrdlf_dcd_tmp_shrub) btran(p) = min(1._r8, btran(p) * 3.33_r8) ! gkw: should we do this for the seasonal deciduous split type? 800:! if (ivt == nbrdlf_dcd_tmp_shrub2) btran(p) = min(1._r8, btran(p) * 3.33_r8) Will think about this later. ========== 20170502: 1. Followed Sarith's instruction to modify the code to read the lightning frequency data. Basically, followed what we do for grn and lai. (1) clsm_ensdrv_mpi.F90: Added to L603: "! real :: lnfm ! lightening frequency [Flashes/km^2/day]" Change L605 "N_real = 2" to "N_real = 3". (2) clsm_ensdrv_drv_routines.F90: subroutine interpolate_to_timestep (3) driver_types.F90 (4) clsm_ensdrv_vegalb_routines.F90: subroutine get_veg_and_alb_times, subroutine get_VEG, subroutine read_veg_or_alb_clim (5) cnlsm_ensdrv_main.F90 2. Search for the routine that calls interpolate_to_timestep and provide this info lnfm_time_old and lnfm_time_new. In ens_driver_routines.F90: Subroutine propagate_ensemble calls interpolate_to_timestep. grn_time_old etc. are input argument in subroutine propagate_ensemble. Subroutine propagate_ensemble is called from Applications/LDAS_App/cnlsm_ensdrv_main.F90 (L2033). ========== 20170503: 1. Continued working on 2 on 20170502. In cnlsm_ensdrv_main.F90: "call get_VEG" to get grn_time_new and grn_time_old etc. Talked to Sarith. He said I can just use grn_time_new and grn_time_old for lnfm_time_new and lnfm_time_old, respectively, because they are both monthly. So modified clsm_ensdrv_drv_routines.F90 to remove lnfm_time_new and lnfm_time_old. 2. Sarith said the case name has to be 3 characters in clsm_ensdrv_vegalb_routines.F90, so modified clsm_ensdrv_vegalb_routines.F90 because I used "LNFM" earlier. Follow what Sarith did for FPAR in /discover/nobackup/smahanam/LDAS/MERRA3_Phase2/LDASsa/src/. ========== 20170504: 1. Changed "lightening" to "lightning" in the four routines below that I modified in the past two days. clsm_ensdrv_mpi.F90 clsm_ensdrv_drv_routines.F90 driver_types.F90 clsm_ensdrv_vegalb_routines.F90 2. Modified the code to read the ABM, PEATF, GDP, HDM, and FC data. (1) catch_types.F90 (2) clsm_ensdrv_mpi.F90 (3) clsm_ensdrv_init_routines.F90: subroutine read_cn_param Below is Sarith's instruction: ! fzeng YOU read CLM4.5_abm_peatf_gdp_hdm_fc read (10,'(2I8, i3, f8.4, f8.2, f10.2, f8.4)' ) i, j, abm(n), peatf(n), & gdp(n), hdm(n), fc(n) 3. Generated the new restart file for the DE_00720x00360_PE_0720x0360.til tile system following Sarith's instruction: cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/src/Components/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts (1) Get the latest copy of mk_LDASsaRestarts.F90 that Sarith modified for CLM4.5 cp -p mk_LDASsaRestarts.F90 mk_LDASsaRestarts.F90.clm4 cvs upd -r 1.3 mk_LDASsaRestarts.F90 (2) Modified subroutine regrid_carbon_vars_clm45 in mk_LDASsaRestarts.F90 to change the initial value for T10 etc. (3) Compiled mk_LDASsaRestarts.F90: Copied mk_LDASsaRestarts.F90 to a compiled GEOS-5 code: cd /discover/nobackup/fzeng/geos5_code/Heracles-5_2/src/GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts cp -p /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/src/Components/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts/mk_LDASsaRestarts.F90 . Compiled: Add mk_LDASsaRestarts to the BINS list in GNUmakefile if it's not there yet. setenv ESMADIR /discover/nobackup/fzeng/geos5_code/Heracles-5_2/ source $ESMADIR/src/g5_modules gmake install (4) Run mk_LDASsaRestarts interactively following the procedures in mk_LDASsaRestarts.F90: cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts mkdir 0.5D ssh discover-sp3 interactive.py -A sp3 -n 96 -a g0620 -X --debug setenv ESMADIR /discover/nobackup/fzeng/geos5_code/Heracles-5_2/ source $ESMADIR/src/g5_modules limit stacksize unlimited cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts/0.5D (removed everything, if any, under this directory) ln -s $ESMADIR/Linux/bin mkdir -p OutData1/ OutData2/ setenv BCSDIR /discover/nobackup/smahanam/bcs/Heracles-4_3/Heracles-4_3_MERRA-3/DE_00720x00360_PE_0720x0360/ ln -s $BCSDIR/DE_00720x00360_PE_0720x0360.til OutData1/OutTileFile ln -s $BCSDIR/DE_00720x00360_PE_0720x0360.til OutData2/OutTileFile ln -s $BCSDIR/clsm OutData2/clsm mpirun -np 96 bin/mk_LDASsaRestarts 50 When done, scale: bin/Scale_CatchCN OutData1/catchcn_internal_rst OutData2/catchcn_internal_rst catchcn_internal_rst 50 SURFLAY: 50.00000 Total Tiles: 67704 Scaled Tiles: 14762 (21%) Performing Sanity Check ... CatDef Tiles: 0 (00%) SrfExc Tiles: 0 (00%) SrfExc Tiles: 0 (00%) bin/Scale_CatchCN OutData1/catchcn_internal_clm45 OutData2/catchcn_internal_clm45 catchcn_internal_clm45 50 SURFLAY: 50.00000 Total Tiles: 67704 Scaled Tiles: 14762 (21%) Performing Sanity Check ... CatDef Tiles: 0 (00%) SrfExc Tiles: 0 (00%) SrfExc Tiles: 0 (00%) (5) Saved a copy of the program and the executable: In /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/src/Components/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts/mk_LDASsaRestarts.F90: L1724: changed " if(MOD(n,10000) == 0) print *,'In Carbon', myid,n" to " if(MOD(n,10000) == 0) print *,'In Carbon CLM4.5', myid,n" cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts cp -p /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/src/Components/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts/mk_LDASsaRestarts.F90 . cp -p /discover/nobackup/fzeng/geos5_code/Heracles-5_2/src/GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts/mk_LDASsaRestarts . NOTE: the executable is older than the program because of the above change which does NOT affect the output restart files. 4. Update the variable names and order in cn_progn_type, cn_param_type, cn_diag_type in catch_types.F90, in clsm_ensdrv_init_routines.F90 and in CN_init (CN_DriverMod.F90). CN_init and CN_exit in CN_DriverMod.F90: updated clsm_ensdrv_init_routines.F90: updated catch_types.F90: updated Searched for T10, PREC10, PREC60 and the ones below and updated them. call MAPL_VarRead ( InNCIO,'TREFAV' , cn_progn(:)%TREFAV ) call MAPL_VarRead ( InNCIO,'t_ref2m_max_inst', cn_progn(:)%t_ref2m_max_inst ) call MAPL_VarRead ( InNCIO,'t_ref2m_min_inst', cn_progn(:)%t_ref2m_min_inst ) call MAPL_VarRead ( InNCIO,'t_ref2m_max', cn_progn(:)%t_ref2m_max ) call MAPL_VarRead ( InNCIO,'t_ref2m_min', cn_progn(:)%t_ref2m_min ) call MAPL_VarRead ( InNCIO,'TDM5' , cn_progn(:)%TDM5 ) call MAPL_VarRead ( InNCIO,'TDM10' , cn_progn(:)%TDM10 ) ========== 20170508: 1. Email from Sarith about mk_LDASsaRestarts.F90: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Hi Fanwei, I made the modifications that I mentioned on top of your modifications and checked in. cvs upd -r 1.3.2.1 mk_LDASsaRestarts.F90 I decided to keep it in a branch for a while until we are confident that this is bug free because the modifications were bit too many to blindly check in as version 1.4 to the trunk. Now this works differently (no need to get an interactive node etc) From the directory on discover where ever you work: (1) setenv ESMADIR /discover/nobackup...... (2) ln -s $ESMADIR/Linux/bin/ (3) mpirun -np 1 bin/mk_LDASsaRestarts -h will print the below message to show the two options that you can run: (1) to create an initial catch(cn)_internal_rst file : ------------------------------------------------------ mk_LDASsaRestarts -a SPONSORCODE -b BCSDIR -m MODEL -s SURFLAY(20/50) -t TILFIL E where MODEL : catch or catchcn (2) to reorder a LDASsa restart file to the order of the BCs : -------------------------------------------------------------- mk_LDASsaRestarts -b BCSDIR -d YYYYMMDD -e EXPNAME -l EXPDIR -m MODEL -s SURFL AY(20/50) -r Y -t TILFILE (4) so for example, if you want to produce restarts for 1/2 degree grid (with your sponsor code) mpirun -np 1 bin/mk_LDASsaRestarts -a s1583 -b /discover/nobackup/smahanam/bcs/Heracles-4_3/Heracles-4_3_MERRA-3/DE_00720x00360_PE_0720x0360/ -t DE_00720x00360_PE_0720x0360.til -m catchcn -s 50 (5) this will quickly produce a job script "mkLDASsa.j" and create OutDatat1/2 directories - you just sbatch mkLDASsa.j ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Fanwei's notes (after checking with Sarith): By doing (1) setenv ESMADIR /discover/nobackup......, (2) ln -s $ESMADIR/Linux/bin/ and (4) in Sarith's email above and submitting mkLDASsa.j, I will create a catchcn_internal_rst file for the 1/2 degree grid using SMAP_M09 restart. This catchcn_internal_rst is ready for use in an offline run. The below command is if I want to reorder an offline restart file to create a restart file for the GCM: mpirun -np 1 bin/mk_LDASsaRestarts -b /discover/nobackup/smahanam/bcs/Heracles-4_3/Heracles-4_3_MERRA-3/DE_00720x00360_PE_0720x0360/ -d YYYYMMDD -e EXPNAME -l EXPDIR -m catchcn -s 50 -r Y -t DE_00720x00360_PE_0720x0360.til YYYYMMDD the start year, month and day of the GCM experiment. If my HHMM is not 0000 (hour/min), I will need to make a small modification to the code. "-r Y" means reorder YES. Here is an example from Sarith - he reordered Qing's SMAP Nature v5: mpirun -np 1 bin/mk_LDASsaRestarts -b /discover/nobackup/projects/gmao/ssd/land/l_data/geos5/bcs/CLSM_params/mkCatchParam_SMAP_L4SM_v002//SMAP_EASEv2_M09/ -t SMAP_EASEv2_M09_3856x1624.til -r Y -m catch -e SMAP_Nature_v05 -l /gpfsm/dnb02/projects/p51/SMAP_Nature/SMAP_Nature_v05/output/SMAP_EASEv2_M09_GLOBAL/ -d 20170101 -s 50 2. Updated my /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/src/Components/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts/mk_LDASsaRestarts.F90: cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/src/Components/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts cvs upd -r 1.3.2.1 mk_LDASsaRestarts.F90 3. Modified the description of CNCOL and CNPFT, and the use/help instructions in mk_LDASsaRestarts.F90. 4. Checked in my changes to mk_LDASsaRestarts.F90 following Sarith's instruction (see his email today): Under /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/src/Components/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts: Checked the history/information of the mk_LDASsaRestarts.F90 I modified: (1) cvs status mk_LDASsaRestarts.F90 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. Could not chdir to home directory /home/fzeng: No such file or directory =================================================================== File: mk_LDASsaRestarts.F90 Status: Locally Modified Working revision: 1.3.2.1 Repository revision: 1.3.2.1 /cvsroot/esma/esma/src/Components/GEOScatch_GridComp/mk_restarts/mk_LDASsaRestarts.F90,v Commit Identifier: 10059106E066FD94CD4 Sticky Tag: 1.3.2.1 Sticky Date: (none) Sticky Options: (none) (2) module load other/tkcvs-8.2.3 (3) tkcvs mk_LDASsaRestarts.F90 & This file is under a branch called "b_mkLDAS". (4) cvs upd -r b_mkLDAS mk_LDASsaRestarts.F90 (5) cvs ci mk_LDASsaRestarts.F90 (6) cvs tag -F H54p3_v24C05_GOSWIM_zerodiff mk_LDASsaRestarts.F90 (7) cvs status mk_LDASsaRestarts.F90 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. Could not chdir to home directory /home/fzeng: No such file or directory =================================================================== File: mk_LDASsaRestarts.F90 Status: Up-to-date Working revision: 1.3.2.2 Repository revision: 1.3.2.2 /cvsroot/esma/esma/src/Components/GEOScatch_GridComp/mk_restarts/mk_LDASsaRestarts.F90,v Commit Identifier: 1005910A9E2553FABF4 Sticky Tag: b_mkLDAS (branch: 1.3.2) Sticky Date: (none) Sticky Options: (none) (1), (2), (3) and (7) are for checking the history/information of the mk_LDASsaRestarts.F90 in /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/src/Components/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts. (4) to (6) are the steps to check in the modifications I made. Sarith pointed out that I need to add "sbatch mkLDAS.j" to the USAGE/HELP and the print statement. Added it and checked in the change using (5) and (6) above (see below. NO need to do (4) again!): cvs ci mk_LDASsaRestarts.F90 cvs tag -F H54p3_v24C05_GOSWIM_zerodiff mk_LDASsaRestarts.F90 ###Add this to my CVS notes### ========== 20170509: 1. Continue to work on 4 on 20170504 above. ========== 20170516: 1. Sarith made further modifications to mk_LDASsaRestarts.F90. Updated my sandbox: /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/src/Components/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts > cvs upd -r 1.4 mk_LDASsaRestarts.F90 Copied the file to /discover/nobackup/fzeng/geos5_code/Heracles-5_2/src/GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts and compiled. Got error message: gmake[2]: *** No rule to make target `stieglitzsnow.mod', needed by `mk_LDASsaRestarts.o'. Stop. Tried using Sarith's GNUmakefile, still failed to compile. Tried copying the file to /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/src/GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts and using Sarith's GNUmakefile. Got error message: gmake[2]: *** No rule to make target `ldasrst2_catch_internal', needed by `bins'. Stop. ========== 20170517: 1. Checked out Icarus-1_0_p1 and compiled it. Will try compiling mk_LDASsaRestarts.F90 here. ========== 20170519: 1. Sarith said mk_LDASsaRestarts.F90 doesn't work in Icarus-1_0_p1. It doesn't include a lot of the development Sarith made. Sarith removed ldasrst2_catch_internal.F90 from GNUmakefile in my /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/src/GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts and now I am able to compile mk_restarts. The executable is /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/Linux/bin/mk_LDASsaRestarts. 2. Run mk_LDASsaRestarts to generate catchcn_internal_rst for 0.5D. cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts mv 0.5D 0.5D_old mkdir 0.5D cd 0.5D setenv ESMADIR /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/ ln -s $ESMADIR/Linux/bin bin mpirun -np 1 bin/mk_LDASsaRestarts -a g0620 -b /discover/nobackup/smahanam/bcs/Heracles-4_3/Heracles-4_3_MERRA-3/DE_00720x00360_PE_0720x0360/ -t DE_00720x00360_PE_0720x0360.til -m catchcn -s 50 sbatch mkLDASsa.j catchcn_internal_rst (for CLM4) and catchcn_internal_clm45 are created successfully. 3. Code added to read abm, gdp, etc from the rst file? Yes. In clsm_ensdrv_init_routines.F90. ========== 20170522: 1. Modified process_cn.F90 and CN_DriverMod.F90 to link the new input data (ABM, PEATF, GDP, HDM, FC) to the corresponding variables in the CN routines. 2. Checked the unit of field capacity in CN routines. It's "watfc(:,:) ! volumetric soil water at field capacity (nlevsoi)" (CNNitrifDenitrifMod.F90 and clmtype.F90). The field capacity we provide has a unit of m3/m3. They look the same. 3. Remove the diagnostics that no longer exist (e.g. fsel) and add the ones needed. Files need to be modified: process_cn.F90: catch_types.F90 (L598-627) clsm_ensdrv_out_routines.F90 (L1501-1535) clsm_ensdrv_glob_param.F90: change N_CatchCN_OutFields to the correct number Diagnostics removed: fsel ========== 20170523: 1. Tried compiling the upgraded offline code and fixed bugs. Was able to compile GEOScatchCN_GridComp and GEOSlana_GridComp, but failed in Applications/LDAS_App. It complained about the C isotope code. Talked to Sarith. Sarith suggests to comment out all the C isotope code. ========== 20170530: 1. Commented out all the C isotope code: CNEcosystemDynMod.F90 clmtypeInitMod.F90 CNSetValueMod.F90 CNPrecisionControlMod.F90 CNSoilLittVertTranspMod.F90 CNWoodProductsMod.F90 CNAllocationMod.F90 CNiniTimeVar.F90 Compile. GEOScatchCN_GridComp: passed GEOSlana_GridComp: passed Applications/LDAS_App: CNEcosystemDynMod.F90:(.text+0x667): undefined reference to `cnharvest_' Commented out "call CNHarvest" in CNEcosystemDynMod.F90. Compile again. Passed! Do test runs tomorrow!! ========== 20170531: 1. CLM4.5 version of Catchment-CN test runs. Domain: 82W-77W, 37N-42N /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/exec/clm4.5/Linux/bin: ./ldsetup setup /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/run/D0.5_clm4.5.exe /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/run/D0.5.bat --runmodel --monthsperjob 24 --landmodel catchCN /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5 > ls -l build/Linux/bin/LDASsaCN_mpi.x -rwxr-xr-x 1 fzeng g0620 70985576 2017-05-30 16:49 build/Linux/bin/LDASsaCN_mpi.x cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/input/restart ln -s /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts/0.5D/catchcn_internal_clm45 catchcn_internal_rst /bin/rm regridded_restarts ssh discover-sp3 interactive.py -A sp3 -n 96 -a g0620 -X --debug setenv ESMADIR /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/ source $ESMADIR/src/g5_modules cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/run ./lenkf.0.i date_time_new 20000101_013000z forrtl: severe (174): SIGSEGV, segmentation fault occurred Image PC Routine Line Source libirc.so 00002AAAAB58E961 Unknown Unknown Unknown libirc.so 00002AAAAB58D0B7 Unknown Unknown Unknown libmpi_usempif08. 00002AAAAC2E1362 Unknown Unknown Unknown libmpi_usempif08. 00002AAAAC2E11B6 Unknown Unknown Unknown libmpi_usempif08. 00002AAAAC2CFA0C Unknown Unknown Unknown libmpi_usempif08. 00002AAAAC2B0BA8 Unknown Unknown Unknown libc.so.6 00002AAAABA14910 Unknown Unknown Unknown LDASsaCN_mpi.x 0000000000761696 cn_drivermod_mp_c 645 CN_DriverMod.F90 L645 in CN_DriverMod.F90: fsat(n) = min(max(0.,car1m(nc)/wtzone(nc,nz)),1.) ========== 20170601: 1. Used Allinea to debug the test run above. See /discover/nobackup/fzeng/notes/debug_procedues for instructions. date_time_new 20000101_013000z forrtl: severe (408): fort: (2): Subscript #1 of the array FSAT has value 1 which is greater than the upper bound of -1 Image PC Routine Line Source libmpi_usempif08. 00002AAAAC2D0093 Unknown Unknown Unknown LDASsaCN_mpi.x 0000000001089AED cn_drivermod_mp_c 645 CN_DriverMod.F90 But L645 in CN_DriverMod.F90 is fsat(n) = min(max(0.,car1m(nc)/wtzone(nc,nz)),1.) The upper bound is 1, not -1. Why? I forgot to allocate and initialize cws%fsat in clmtypeInitMod.F90. Fixed this and tried. It passed that point but crashed here: date_time_new 20000101_013000z forrtl: error (65): floating invalid Image PC Routine Line Source libirc.so 00002AAAAB58E961 Unknown Unknown Unknown libirc.so 00002AAAAB58D0B7 Unknown Unknown Unknown libmpi_usempif08. 00002AAAAC2E1362 Unknown Unknown Unknown libmpi_usempif08. 00002AAAAC2E11B6 Unknown Unknown Unknown libmpi_usempif08. 00002AAAAC2CFA0C Unknown Unknown Unknown libmpi_usempif08. 00002AAAAC2B0E82 Unknown Unknown Unknown libc.so.6 00002AAAABA14910 Unknown Unknown Unknown LDASsaCN_mpi.x 0000000002D50C56 Unknown Unknown Unknown LDASsaCN_mpi.x 0000000001DCE99B cndecompcascademo 629 CNDecompCascadeMod_BGC.F90 L629 in CNDecompCascadeMod_BGC.F90: w_scalar(c,1) = w_scalar(c,1) + (log(minpsi/psi)/log(minpsi/maxpsi))*fr(c,j) maxpsi = sucsat(c,j) * (-9.8e-6_r8) maxpsi = 0.00243236008 (NOTE: positive) minpsi = -10 psi = -0.00243401062 sucsat: -500 to -150 psis is soil water potential. sucsat is minimum soil suction (mm). In CN_DriverMod.F90: sucsat(n,1) = psis(nc) * 1e3 ! minimum soil suction (mm), psis is in (m) ========== 20170606: 1. Tried changing "sucsat(n,1) = psis(nc) * 1e3" to "sucsat(n,1) = psis(nc) * 1e3 * (-1)" in CN_DriverMod. forrtl: severe (408): fort: (3): Subscript #1 of the array TILE_COORD has value 0 which is less than the lower bound of 1 Image PC Routine Line Source libmpi_usempif08. 00002AAAAC2D0093 Unknown Unknown Unknown LDASsaCN_mpi.x 0000000000AFD68C tile_coord_routin 1305 tile_coord.F90 LDASsaCN_mpi.x 0000000000763F97 clsm_ensdrv_init_ 944 clsm_ensdrv_init_routines.F90 LDASsaCN_mpi.x 00000000004A2F71 MAIN__ 828 cnlsm_ensdrv_main.F90 ========== 20170607: 1. More CLM4.5 Catchment-CN test runs and debugging. Created a new experiment by running ldsetup under LDAS_Apps and the above issue about tile_coord.F90 is resolved. But new error: date_time_new 20000101_013000z forrtl: severe (408): fort: (3): Subscript #3 of the array DECOMP_CPOOLS_VR has value 0 which is less than the lower bound of 1 Image PC Routine Line Source libmpi_usempif08. 00002AAAAC2D0093 Unknown Unknown Unknown LDASsaCN_mpi.x 00000000015BF4B1 cndecompmod_mp_cn 233 CNDecompMod.F90 L233 CNDecompMod.F90: if (decomp_cpools_vr(c,j,cascade_donor_pool(k)) > 0._r8 .and. & Allinea DDT shows that: cascade_donor_pool, dimension 1x9, is all 0. In CNDecompCascadeMod_BGC.F90, subroutine init_decompcascade: i_litr1 = i_met_lit i_litr2 = i_cel_lit i_litr3 = i_lig_lit i_soil1 = 5 i_soil2 = 6 i_soil3 = 7 i_soil4 = 8 i_l1s1 = 1 cascade_donor_pool(i_l1s1) = i_litr1 cascade_receiver_pool(i_l1s1) = i_soil1 pathfrac_decomp_cascade(begc:endc,1:nlevdecomp,i_l1s1) = 1.0_r8 i_l2s2 = 2 cascade_donor_pool(i_l2s2) = i_litr2 cascade_receiver_pool(i_l2s2) = i_soil2 pathfrac_decomp_cascade(begc:endc,1:nlevdecomp,i_l2s2) = 1.0_r8 i_l3s3 = 3 cascade_donor_pool(i_l3s3) = i_litr3 cascade_receiver_pool(i_l3s3) = i_soil3 pathfrac_decomp_cascade(begc:endc,1:nlevdecomp,i_l3s3) = 1.0_r8 i_s1s2 = 4 cascade_donor_pool(i_s1s2) = i_soil1 cascade_receiver_pool(i_s1s2) = i_soil2 pathfrac_decomp_cascade(begc:endc,1:nlevdecomp,i_s1s2) = 1.0_r8 i_s2s3 = 5 cascade_donor_pool(i_s2s3) = i_soil2 cascade_receiver_pool(i_s2s3) = i_soil3 pathfrac_decomp_cascade(begc:endc,1:nlevdecomp,i_s2s3) = 1.0_r8 i_s3s4 = 6 cascade_donor_pool(i_s3s4) = i_soil3 cascade_receiver_pool(i_s3s4) = i_soil4 pathfrac_decomp_cascade(begc:endc,1:nlevdecomp,i_s3s4) = 1.0_r8 i_s4atm = 7 cascade_donor_pool(i_s4atm) = i_soil4 cascade_receiver_pool(i_s4atm) = i_atm pathfrac_decomp_cascade(begc:endc,1:nlevdecomp,i_s4atm) = 1.0_r8 i_cwdl2 = 8 cascade_donor_pool(i_cwdl2) = i_cwd cascade_receiver_pool(i_cwdl2) = i_litr2 pathfrac_decomp_cascade(begc:endc,1:nlevdecomp,i_cwdl2) = cwd_fcel i_cwdl3 = 9 cascade_donor_pool(i_cwdl3) = i_cwd cascade_receiver_pool(i_cwdl3) = i_litr3 pathfrac_decomp_cascade(begc:endc,1:nlevdecomp,i_cwdl3) = cwd_flig Tried in clm_varpar.F90: !#ifndef CENTURY_DECOMP ! ! parameters for decomposition cascade ! integer, parameter :: ndecomp_pools = 8 ! integer, parameter :: ndecomp_cascade_transitions = 9 ! integer, parameter :: i_met_lit = 1 ! integer, parameter :: i_cel_lit = 2 ! integer, parameter :: i_lig_lit = 3 ! integer, parameter :: i_cwd = 4 ! integer, parameter :: nsompools = 4 !#else ! ! parameters for decomposition cascade ! integer, parameter :: ndecomp_pools = 7 ! integer, parameter :: ndecomp_cascade_transitions = 10 ! integer, parameter :: i_met_lit = 1 ! integer, parameter :: i_cel_lit = 2 ! integer, parameter :: i_lig_lit = 3 ! integer, parameter :: i_cwd = 4 ! integer, parameter :: nsompools = 3 !#endif ! repeat the first block above here in case the above directive is not compiled/used at all ! fzeng, 7 June 2017 ! parameters for decomposition cascade integer, parameter, public :: ndecomp_pools = 8 integer, parameter, public :: ndecomp_cascade_transitions = 9 integer, parameter, public :: i_met_lit = 1 integer, parameter, public :: i_cel_lit = 2 integer, parameter, public :: i_lig_lit = 3 integer, parameter, public :: i_cwd = 4 integer, parameter, public :: nsompools = 4 Still crashed giving the same error message. 2. Met with Randy about the sucsat issue above. ========== 20170608: 1. Met with Randy about the sucsat issue above. 2. Found that subroutine init_decompcascade is not called anywhere in Catchment-CN!!! In CLM4.5, it's called in main/iniTimeConst.F90. There is no iniTimeConst.F90 in Catchment-CN. Tried calling subroutine init_decompcascade in CN_DriverMod. This seems to solve the problem, but the test run crashed somewhere else. date_time_new 20000101_013000z forrtl: severe (408): fort: (2): Subscript #1 of the array DZSOI_DECOMP has value 1 which is greater than the upper bound of -1 Image PC Routine Line Source libmpi_usempif08. 00002AAAAC2D0093 Unknown Unknown Unknown LDASsaCN_mpi.x 0000000001D011EE cnverticalprofile 271 CNVerticalProfileMod.F90 LDASsaCN_mpi.x 0000000001A3EF41 cndecompmod_mp_cn 307 CNDecompMod.F90 L271 CNVerticalProfileMod.F90: ndep_prof_sum = ndep_prof_sum + ndep_prof(c,j) * dzsoi_decomp(j) In clm_varcon.F90 subroutine clm_varcon_init: allocate( dzsoi_decomp(1:nlevdecomp_full) ) dzsoi_decomp(1) = dzsoi(1) But, subroutine clm_varcon_init is not called (it's called in clm_initializeMod.F90 in CLM4.5, and there is no clm_initializeMod.F90 in Catchment-CN). clm_varpar_init is called in clm_initializeMod.F90 in CLM4.5 and it's called in subroutine CN_init in CN_DriverMod.F90. Therefore, also call clm_varcon_init in subroutine CN_init in CN_DriverMod.F90. Also moved "call init_decompcascade" from subroutine CN_Driver to CN_init. This seems to solve the problem, but the test run crashed with a new error message: date_time_new 20000101_013000z forrtl: severe (408): fort: (2): Subscript #1 of the array COLARR has value 1 which is greater than the upper bound of -1 Image PC Routine Line Source libmpi_usempif08. 00002AAAAC2D0093 Unknown Unknown Unknown LDASsaCN_mpi.x 00000000018127B7 subgridavemod_mp_ 245 subgridAveMod.F90 LDASsaCN_mpi.x 0000000001858E00 cnallocationmod_m 847 CNAllocationMod.F90 LDASsaCN_mpi.x 0000000001A3F019 cndecompmod_mp_cn 319 CNDecompMod.F90 L245 subgridAveMod.F90: colarr(c) = 0. There is little change to subgridAveMod.F90 in Catchment-CN from CLM4 version to CLM4.5. Why does it crash here? ========== 20170609: 1. Continued working on the issue about subgridAveMod.F90 above. Checked the difference in this file between original CLM4 and CLM4.5 code. Found that in CLM4.5, pactive is needed in subgridAveMod.F90, and a few other routines as well such as CNFireMod.F90. In CLM4.5, pactive is computed in reweightMod.F90, which needs filterMod, domainMod and decompMod which are not included in Catchment-CN. These should be added, but not sure if this is causing the issue. Found that in CLM4 version of Catchment-CN CNAllocationMod.F90 (L419-421): ! now use the p2c routine to get the column-averaged plant_ndemand allocate(col_plant_ndemand(lbc:ubc)) call p2c(num_soilc,filter_soilc,plant_ndemand,col_plant_ndemand) col_plant_ndemand is allocated before "call p2c". This is the same as the original CLM4 code. While in the original CLM4.5 CNAllocationMod.F90 (L831-832): ! now use the p2c routine to get the column-averaged plant_ndemand call p2c(num_soilc,filter_soilc,plant_ndemand,col_plant_ndemand) col_plant_ndemand is NOT allocated before "call p2c". Tried adding "allocate(col_plant_ndemand(lbc:ubc))" before "call p2c(num_soilc,filter_soilc,plant_ndemand,col_plant_ndemand)" in CLM4.5 version of Catchment-CN CNAllocationMod.F90. This solved the issue. 2. New issue: L323 in CNPhenologyMod.F90: tempavg_t2m(p) = tempavg_t2m(p) + t_ref2m(p) * (fracday/dayspyr) date_time_new 20000101_013000z dayspyr: 366 (correct) fracday: nan (floating invalid) Why is fracday nan? L180-181 in subroutine CNPhenologyInit in CNPhenologyMod.F90: dt = real( get_step_size(), r8 ) fracday = dt/secspday CLM4.5 version of Catchment-CN: subroutine CNPhenologyInit is called from subroutine CNEcosystemDynInit in CNEcosystemDynMod.F90, but subroutine CNEcosystemDynInit is not called anywhere. Original CLM4.5 code: subroutine CNPhenologyInit is called from subroutine CNEcosystemDynInit in CNEcosystemDynMod.F90, and subroutine CNEcosystemDynInit is called from clm_initializeMod.F90 (which is not included in Catchment-CN). subroutine CNEcosystemDynInit does not exist in the original CLM4 code. Tried adding "call CNEcosystemDynInit" to CN_init. date_time_new 20000101_000730z CN: dt_default < 0 Tried moving "call CNEcosystemDynInit" from CN_init to CN_Driver. The issue is resolved. 3. New issue: date_time_new 20000101_013000z forrtl: severe (408): fort: (2): Subscript #1 of the array FORC_T has value 1 which is greater than the upper bound of -1 Image PC Routine Line Source libmpi_usempif08. 00002AAAAC2D0093 Unknown Unknown Unknown LDASsaCN_mpi.x 0000000001A8D120 cnfiremod_mp_cnfi 568 CNFireMod.F90 LDASsaCN_mpi.x 00000000015ED7B2 cnecosystemdynmod 254 CNEcosystemDynMod.F90 LDASsaCN_mpi.x 00000000010AFE3A cn_drivermod_mp_c 819 CN_DriverMod.F90 LDASsaCN_mpi.x 0000000000B4F103 process_cn_mp_pro 1644 process_cn.F90 Added this block below to clmtypeInitMod.F90: allocate(grc%forc_rh (beg:end)) allocate(grc%forc_wind(beg:end)) allocate(grc%forc_t (beg:end)) allocate(grc%forc_rain(beg:end)) allocate(grc%forc_snow(beg:end)) 4. New issue: date_time_new 20000101_013000z floating point invalid operation in L636 CNFireMod.F90: fb = max(0.0_r8,min(1.0_r8,(fuelc(c)-lfuel)/(ufuel-lfuel))) fuelc(1) = nan fuelc(2:end) = 9.99999962e+35 L632-634: do j = 1, nlevdecomp fuelc(c) = fuelc(c)+decomp_cpools_vr(c,j,i_cwd) * dzsoi_decomp(j) end do j=2, nlevdecomp=1. How could this happen? L631: fuelc(c) = totlitc(c)+totvegc_col(c)-rootc_col(c)-fuelc_crop(c)*cropf_col(c) totvegc_col and rootc_col are nan. L359: call p2c(num_soilc, filter_soilc,totvegc, totvegc_col) Both totvegc and totvegc_col are nan. Checked the totvegc in the restart file using ~/matlab/plot_catchcn_internal_rst_vars.m and found it 0 everywhere. Why? ========== 20170612: 1. Used "ncdump -v" command to check the values in /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts/0.5D/catchcn_internal_clm45. ITY: good FVG: good NDEP: good HDM: 0 GDP: 0 PEATF: 0 CNCOL: 0 CNPFT: 0 Asked Sarith. He is working from home today. We will look into this tomorrow. ========== 20170616: 1. /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/src/Components/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts: cvs upd -r 1.4.2.1 mk_LDASsaRestarts.F90 NOTE: use the catchcn_internal_rst in OutData2 (not the final one) for now ========== 20170620: 1. Copied /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/src/Components/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts/mk_LDASsaRestarts.F90 to /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/src/GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/mk_restarts . The executable is /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/Linux/bin/mk_LDASsaRestarts. 2. Run mk_LDASsaRestarts to generate catchcn_internal_rst for 0.5D. cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts mv 0.5D 0.5D_wrong mkdir 0.5D cd 0.5D setenv ESMADIR /discover/nobackup/fzeng/clm4-to-clm4.5/GEOS5/Heracles-5_4_p3-M3_V24_C05/ ln -s $ESMADIR/Linux/bin bin mpirun -np 1 bin/mk_LDASsaRestarts -a g0620 -b /discover/nobackup/smahanam/bcs/Heracles-4_3/Heracles-4_3_MERRA-3/DE_00720x00360_PE_0720x0360/ -t DE_00720x00360_PE_0720x0360.til -m catchcn -s 50 sbatch mkLDASsa.j NOTE: use /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts/0.5D/OutData2/catchcn_internal_clm45 for now. ========== 20170621: 1. Continued working on CLM4.5 Catchment-CN test runs. cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/input/restart rm catchcn_internal_rst ln -s /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts/0.5D/OutData2/catchcn_internal_clm45 catchcn_internal_rst Followed notes in /discover/nobackup/fzeng/notes/debug_procedures: cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/run /bin/rm ../output/global/rc_out/* module load tool/allinea-tools-7.0.4 ddt /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/build/Linux/bin/LDASsaCN_mpi.x & (Note: the order is important. If do "module load tool/allinea-tools-6.1.0" before "source ../build/Linux/bin/g5_modules", the latter will unload the module.) Copy the part of lenkf.1.i from work_path to the end (no need to be in one line) to arguments. Copy this part of lenkf.1.i to Environment variables: ~~~~~~~~~~~~~~~~~~~~~~ setenv MKL_CBWR SSE4_2 # ensure zero-diff across archs setenv MV2_ON_DEMAND_THRESHOLD 8192 # MVAPICH2 setenv MYNAME `finger $USER | cut -d':' -f3 | head -1` ~~~~~~~~~~~~~~~~~~~~~~ click "Run" click the play button on the upper left to make it run until it stops. Stopped at L1235 in CN_DriverMod at the first time step: ccs%prod100c (n) = cncol(nc,nz, 7) Checked the restart file. prod100c is 0 everywhere, but it's supposed to be like this. totcolc, cannsum_npp look correct. Found that it's a breakpoint (I didn't set this though) and that's why it stopped. Clicked "continue" and the run crashed at L1877 in catchmentCN.F90: DTCWET = ( Q1*A21 - A11*Q2 ) / ( A12*A21 - A11*A22 ) A11 = 7.74059213e+09 A12 = -7.74059213e+09 A21 = -7.74059213e+09 A22 = 7.74059213e+09 So A12*A21 - A11*A22 = 0. Repeated the run (removed the breakpoint at L1235 in CN_DriverMod) and got the same error message. It's because I didn't set the restart path correctly!! Now it crashed at L636 in CNFireMod, saying floating invalid: fb = max(0.0_r8,min(1.0_r8,(fuelc(c)-lfuel)/(ufuel-lfuel))) In L631 in CNFireMod: fuelc(c) = totlitc(c)+totvegc_col(c)-rootc_col(c)-fuelc_crop(c)*cropf_col(c) totlitc looks correct: 180-260 gC/m2. totvegc_col: nan rootc_col: nan fuelc_crop: 0 cropf_col: 0 ========== 20170622: 1. Continued working on the issue above. totvegc_col => ccs%totvegc_col is from the restart file (see subroutine CN_init), but it's then computed in CNFireMod in original CLM4.5: call p2c(num_soilc, filter_soilc,totvegc, totvegc_col) totvegc => pcs%totvegc should be from the restart file, see below the original CLM4.5 code in CNrestMod: call ncd_io(varname='totvegc', data=pcs%totvegc, & dim1name=namep, ncid=ncid, flag=flag, readvar=readvar) However, we don't have pft-level totvegc in the restart file. That's why the totvegc_col became nan after call p2c in CNFireMod. In CLM4 CNrestMod, totvegc is also a pft-level variable in the restart file (same as in CLM4.5). Greg didn't include it in the restart file probably because it's not used anywhere before it's computed in subroutine CSummary (called from CNEcosystemDyn). In CLM4.5 CNiniTimeVar.F90: totvegc => pcs%totvegc ! calculate totvegc explicitly so that it is available for the isotope ! code on the first time step. totvegc(p) = leafc(p) + leafc_storage(p) + leafc_xfer(p) + frootc(p) + & frootc_storage(p) + frootc_xfer(p) + livestemc(p) + livestemc_storage(p) + & livestemc_xfer(p) + deadstemc(p) + deadstemc_storage(p) + deadstemc_xfer(p) + & livecrootc(p) + livecrootc_storage(p) + livecrootc_xfer(p) + deadcrootc(p) + & deadcrootc_storage(p) + deadcrootc_xfer(p) + gresp_storage(p) + & gresp_xfer(p) + cpool(p) if ( crop_prog )then totvegc(p) = totvegc(p) + grainc(p) + grainc_storage(p) + grainc_xfer(p) end if Need to go over CLM4.5 CNrestMod to see what else to be included in or excluded from our restart file? ========== 20170623: 1. Since pft-level totvegc can be computed from pft-level leafc etc (see above), added a block in subroutine CN_init to compute pft-level totvegc from leafc, leafc_storage etc. Now totvegc_col in CNFireMod looks correct. ccs%totvegc_col from the restart file: 27232.4336 24376.1758 19026.2969 ccs%totvegc_col computed by p2c in CNFireArea: 26568.6953 24030.6133 19063.4062 They are close to each other. However, rootc_col in CNFireMod is nan. L469 in CNFireMod: rootc_col(c) = rootc_col(c) + (frootc(p) + frootc_storage(p) + & frootc_xfer(p) + deadcrootc(p) + & deadcrootc_storage(p) + deadcrootc_xfer(p) + & livecrootc(p)+livecrootc_storage(p) + & livecrootc_xfer(p))*wtcol(p) Because not all 28 PFTs co-exist in each column, frootc(p) is nan if the corresponding PFT does not exist in the column, and therefore doing sum this way leads to nan in rootc_col. Added "if (wtcol(p) > 0.) then" before the block of rootc_col calculation above in CNFireMod. Now rootc_col in CNFireMod looks correct. However, fuelc is still nan. L635 in CNFireMod: fuelc(c) = fuelc(c)+decomp_cpools_vr(c,j,i_cwd) * dzsoi_decomp(j) decomp_cpools_vr is nan in the 1st and 3rd columns for cwdc, litr1c, litr2c and litr3c here. decomp_cpools_vr in CN_init looks correct. Why does it become nan in CNFireMod? In subroutine CNSoilLittVertTransp (called before CNFireArea is called from CNEcosystemDyn): conc_ptr => ccs%decomp_cpools_vr source => ccf%decomp_cpools_sourcesink !! for single level case, no transport; just update the fluxes calculated in the StateUpdate1 subroutines do l = 1, ndecomp_pools do j = 1,nlevdecomp do fc = 1, num_soilc c = filter_soilc (fc) conc_ptr(c,j,l) = conc_ptr(c,j,l) + source(c,j,l) trcr_tendency_ptr(c,j,l) = 0._r8 end do end do end do ========== 20170626: 1. Continued looking for the code that computes source above. subroutine CStateUpdate1 in CNCStateUpdate1Mod.F90: case ('bulk') pcisof => pcf pcisos => pcs ccisof => ccf ccisos => ccs decomp_cpools_vr => ccisos%decomp_cpools_vr decomp_cpools_sourcesink => ccisof%decomp_cpools_sourcesink ! plant to litter fluxes do j = 1,nlevdecomp ! column loop do fc = 1,num_soilc c = filter_soilc(fc) ! phenology and dynamic land cover fluxes decomp_cpools_sourcesink(c,j,i_met_lit) = ( phenology_c_to_litr_met_c(c,j) + dwt_frootc_to_litr_met_c(c,j) ) *dt decomp_cpools_sourcesink(c,j,i_cel_lit) = ( phenology_c_to_litr_cel_c(c,j) + dwt_frootc_to_litr_cel_c(c,j) ) *dt decomp_cpools_sourcesink(c,j,i_lig_lit) = ( phenology_c_to_litr_lig_c(c,j) + dwt_frootc_to_litr_lig_c(c,j) ) *dt decomp_cpools_sourcesink(c,j,i_cwd) = ( dwt_livecrootc_to_cwdc(c,j) + dwt_deadcrootc_to_cwdc(c,j) ) *dt end do end do ! litter and SOM HR fluxes do k = 1, ndecomp_cascade_transitions do j = 1,nlevdecomp ! column loop do fc = 1,num_soilc c = filter_soilc(fc) decomp_cpools_sourcesink(c,j,cascade_donor_pool(k)) = & decomp_cpools_sourcesink(c,j,cascade_donor_pool(k)) & - ( decomp_cascade_hr_vr(c,j,k) + decomp_cascade_ctransfer_vr(c,j,k)) *dt end do end do end do do k = 1, ndecomp_cascade_transitions if ( cascade_receiver_pool(k) .ne. 0 ) then ! skip terminal transitions do j = 1,nlevdecomp ! column loop do fc = 1,num_soilc c = filter_soilc(fc) decomp_cpools_sourcesink(c,j,cascade_receiver_pool(k)) = & decomp_cpools_sourcesink(c,j,cascade_receiver_pool(k)) & + decomp_cascade_ctransfer_vr(c,j,k)*dt end do end do end if end do decomp_cpools_sourcesink: only has valid values for the 2nd column, and is nan for the 1st and 3rd columns. phenology_c_to_litr_met_c, phenology_c_to_litr_cel_c and phenology_c_to_litr_lig_c: (1,1): nan (2,1): 0 (3,1): nan dwt_frootc_to_litr_met_c: (1,1): 0 (2,1): 0 (3,1): 0 So the problem is in phenology_c_to_litr_met_c, phenology_c_to_litr_cel_c and phenology_c_to_litr_lig_c. phenology_c_to_litr_met_c => ccisof%phenology_c_to_litr_met_c phenology_c_to_litr_cel_c => ccisof%phenology_c_to_litr_cel_c phenology_c_to_litr_lig_c => ccisof%phenology_c_to_litr_lig_c Subroutine CNLitterToColumn in CNPhenologyMod: phenology_c_to_litr_met_c => ccf%phenology_c_to_litr_met_c phenology_c_to_litr_cel_c => ccf%phenology_c_to_litr_cel_c phenology_c_to_litr_lig_c => ccf%phenology_c_to_litr_lig_c do j = 1, nlevdecomp do pi = 1,max_pft_per_col do fc = 1,num_soilc c = filter_soilc(fc) if ( pi <= npfts(c) ) then p = pfti(c) + pi - 1 if (pactive(p)) then ! leaf litter carbon fluxes phenology_c_to_litr_met_c(c,j) = phenology_c_to_litr_met_c(c,j) & + leafc_to_litter(p) * lf_flab(ivt(p)) * wtcol(p) * leaf_prof(p,j) phenology_c_to_litr_cel_c(c,j) = phenology_c_to_litr_cel_c(c,j) & + leafc_to_litter(p) * lf_fcel(ivt(p)) * wtcol(p) * leaf_prof(p,j) phenology_c_to_litr_lig_c(c,j) = phenology_c_to_litr_lig_c(c,j) & + leafc_to_litter(p) * lf_flig(ivt(p)) * wtcol(p) * leaf_prof(p,j) ! fine root litter carbon fluxes phenology_c_to_litr_met_c(c,j) = phenology_c_to_litr_met_c(c,j) & + frootc_to_litter(p) * fr_flab(ivt(p)) * wtcol(p) * froot_prof(p,j) phenology_c_to_litr_cel_c(c,j) = phenology_c_to_litr_cel_c(c,j) & + frootc_to_litter(p) * fr_fcel(ivt(p)) * wtcol(p) * froot_prof(p,j) phenology_c_to_litr_lig_c(c,j) = phenology_c_to_litr_lig_c(c,j) & + frootc_to_litter(p) * fr_flig(ivt(p)) * wtcol(p) * froot_prof(p,j) ! agroibis puts crop stem litter together with leaf litter ! so I've used the leaf lf_f* parameters instead of making ! new ones for now (slevis) ! also for simplicity I've put "food" into the litter pools if (ivt(p) >= npcropmin) then ! add livestemc to litter ! stem litter carbon fluxes phenology_c_to_litr_met_c(c,j) = phenology_c_to_litr_met_c(c,j) & + livestemc_to_litter(p) * lf_flab(ivt(p)) * wtcol(p) * leaf_prof(p,j) phenology_c_to_litr_cel_c(c,j) = phenology_c_to_litr_cel_c(c,j) & + livestemc_to_litter(p) * lf_fcel(ivt(p)) * wtcol(p) * leaf_prof(p,j) phenology_c_to_litr_lig_c(c,j) = phenology_c_to_litr_lig_c(c,j) & + livestemc_to_litter(p) * lf_flig(ivt(p)) * wtcol(p) * leaf_prof(p,j) ! grain litter carbon fluxes phenology_c_to_litr_met_c(c,j) = phenology_c_to_litr_met_c(c,j) & + grainc_to_food(p) * lf_flab(ivt(p)) * wtcol(p) * leaf_prof(p,j) phenology_c_to_litr_cel_c(c,j) = phenology_c_to_litr_cel_c(c,j) & + grainc_to_food(p) * lf_fcel(ivt(p)) * wtcol(p) * leaf_prof(p,j) phenology_c_to_litr_lig_c(c,j) = phenology_c_to_litr_lig_c(c,j) & + grainc_to_food(p) * lf_flig(ivt(p)) * wtcol(p) * leaf_prof(p,j) At the beginning (before going into the do loops above): phenology_c_to_litr_met_c, phenology_c_to_litr_cel_c and phenology_c_to_litr_lig_c: (1,1): 0 (2,1): 0 (3,1): 0 In the first step of the above do loop: pi = 5 (why?) pfti(1) = 1; pfti(2) = 29; pfti(3) = 57; p = 5 pactive: (1): True (2): True (3): True (4): False (5): True ... However, leafc_to_litter has valid value in 2, 8, 30, 36, 58 and 64, and is nan everywhere else. Vegetation types in this tile: 1 and 7. leafc_to_litter looks correct. pactive may be wrong. Check our code and can't find the calculation of pactive anywhere. So its values above may be randomly assigned and therefore incorrect. In CLM4.5, pactive is computed in main/reweightMod.F90, which needs main/filterMod, cesm1_2_2/models/lnd/clm/src/util_share/domainMod and cesm1_2_2/models/lnd/clm/src/util_share/decompMod. May not need decompMod/get_clump_bounds -- we provide begg, endg, begl, endl, begc, endc, begp, endp domainMod: ldomain is defined but seems not used (?) Tomorrow, try setting pactive in CN_Driver: as long as fveg>0, pactive = .true. ========== 20170627: 1. Set pactive to true in CN_Driver if fveg > 1.e-4. This seems to solve the problem above. Now the problem is in L303 in CNVegStructUpdateMod.F90: floating point invalid ol = min( max(snow_depth(c)-hbot(p), 0._r8), htop(p)-hbot(p)) hbot and htop have valid values. snow_depth is nan for all 3 columns. We have calculation of snowdp (gridcell averaged snow height) in CN_Driver, but no calculation of snow_depth (snow height of snow covered area) anywhere yet. Add calculation of snow_depth to CN_Driver. In CLM4, there is only snowdp. In CLM4.5, there are both snowdp and snow_depth. From code comparison, snowdp in CLM4 is the same as snow_depth in CLM4.5. snowdp in CLM4.5 seems not used anywhere. Use snow_depth = sndzn in CN_Driver for now. It's running! Talked to Randy about this and we decided to use the column averaged snow hight for snow_depth. In CN_Driver, changed snow_depth = sndzn to snow_depth = sndzn*asnow and deleted snowdp. ========== 20170628: 1. Checked below: (1) runsrf in the Catchment model has a unit of kg/m2/s, same as mm H2O/s used in CLM4.5 qflx_surf. See GEOSlana_GridComp/clsm_ensdrv_out_routines.F90. (2) In the N leaching code, both sub-surface runoff (qflx_drain = bflow) and surface runoff (qflx_surf = runsrf) are used. (3) sndz used in CN_Driver is the sum the thickness of 3 snow layers, see the code in process_cn.F90 below: do n = 1,N_snow cn_progn(:)%sndzm = cn_progn(:)%sndzm + sndz(n,:) end do 2. Went through the To Do list in /discover/nobackup/fzeng/clm4-to-clm4.5/notes/notes: (1) runsrf: Checked with Randy yesterday. It's ok that runsrf is 0 at the very first model time step of very job segment. Therefore, the current implementation is fine. (2) btran2x (used as btran2 in CNFireMod): currently it's the same as btran which is from the CLM4 version of Catchment-CN, since btran is already the same as btran2 in CNFireMod. (3) tsoi17: need to ask Randy. (4) Correct in CN_DriverMod? npp => pcf%npp? gpp => pcf%gpp? I think so. (5) organized the pointers in CN_DriverMod. (6) About zsoi, dzsoi and zisoi used in subroutine decomp_vertprofiles in CNVerticalProfileMod.F90: Temporary solution: (a) Turn off VERTSOILC and (b) set these in clm_varcon.F90: zsoi(1) = 0.5 dzsoi(1) = 1. zisoi(0) = 0. zisoi(1) = 1. dzsoi_decomp(1) = dzsoi(1) dzsoi, zisoi and dzsoi_decomp look correct (i.e. same as the CLM4.5) for just one soil layer with thickness of 1m. The value of zsoi 0.5m set here is different from that (7mm) if we set nlevgrnd=1 in CLM4.5. zsoi is only used if VERTSOILC is turned on in CNSoilLittVertTranspMod.F90. So don’t worry about it for now. (7) Accumulated fields in CLM4.5: main/restFileMod.F90, subroutine restFile_read: call accumulRest( ncid, flag='read' ) Subroutine accumulRest is located in /discover/nobackup/fzeng/clm_orig/cesm1_2_2/models/lnd/clm/src/util_share/accumulMod.F90. It reads/writes accumulation restart data. So yes the accumulation data is in the restart file in CLM4.5. ========== 20170629: 1. Continued going through the To Do list in /discover/nobackup/fzeng/clm4-to-clm4.5/notes/notes: (1) Continued working on accumulated fields in CLM4.5: In main/clm_initializeMod.F90, subroutine initialize2: L493: call restFile_read( fnamer ) ! fzeng: this calls accumulRest( ncid, flag='read' ) L549-554: ! Initialize clmtype variables that are obtained from accumulated fields. ! This routine is called in an initial run at nstep=0 ! This routine is also always called for a restart run and must ! therefore be called after the restart file is read in call initAccClmtype() In subroutine initAccClmtype (in main/accFldsMod.F90): ! Initialize 2m ref temperature max and min values if (nsrest == nsrStartup) then ! Why not restart&branch? These vars are not in clmr. do p = begp,endp t_ref2m_max(p) = spval t_ref2m_min(p) = spval t_ref2m_max_inst(p) = -spval t_ref2m_min_inst(p) = spval t_ref2m_max_u(p) = spval t_ref2m_max_r(p) = spval t_ref2m_min_u(p) = spval t_ref2m_min_r(p) = spval t_ref2m_max_inst_u(p) = -spval t_ref2m_max_inst_r(p) = -spval t_ref2m_min_inst_u(p) = spval t_ref2m_min_inst_r(p) = spval end do end if t_ref2m_max, t_ref2m_min, t_ref2m_max_inst, and t_ref2m_min_inst are not in the accumulRest file. So first the accumul rest file is read in, then t_ref2m_max, t_ref2m_min, t_ref2m_max_inst, and t_ref2m_min_inst are initialized to spval or -spval as above. Then in main/accFldsMod.F90 subroutine updateAccFlds: L477-501: ! Accumulate and extract TREFAV - hourly average 2m air temperature ! Used to compute maximum and minimum of hourly averaged 2m reference ! temperature over a day. Note that "spval" is returned by the call to ! accext if the time step does not correspond to the end of an ! accumulation interval. First, initialize the necessary values for ! an initial run at the first time step the accumulator is called call update_accum_field ('TREFAV', t_ref2m, nstep) call extract_accum_field ('TREFAV', rbufslp, nstep) end_cd = is_end_curr_day() do p = begp,endp if (rbufslp(p) /= spval) then t_ref2m_max_inst(p) = max(rbufslp(p), t_ref2m_max_inst(p)) t_ref2m_min_inst(p) = min(rbufslp(p), t_ref2m_min_inst(p)) endif if (end_cd) then t_ref2m_max(p) = t_ref2m_max_inst(p) t_ref2m_min(p) = t_ref2m_min_inst(p) t_ref2m_max_inst(p) = -spval t_ref2m_min_inst(p) = spval else if (secs == int(dtime)) then t_ref2m_max(p) = spval t_ref2m_min(p) = spval endif end do L677-693: ! Accumulate and extract TDM10 do p = begp,endp rbufslp(p) = min(t_ref2m_min(p),t_ref2m_min_inst(p)) !slevis: ok choice? if (rbufslp(p) > 1.e30_r8) rbufslp(p) = SHR_CONST_TKFRZ !and were 'min'& end do !'min_inst' not initialized? call update_accum_field ('TDM10', rbufslp, nstep) call extract_accum_field ('TDM10', a10tmin, nstep) ! Accumulate and extract TDM5 do p = begp,endp rbufslp(p) = min(t_ref2m_min(p),t_ref2m_min_inst(p)) !slevis: ok choice? if (rbufslp(p) > 1.e30_r8) rbufslp(p) = SHR_CONST_TKFRZ !and were 'min'& end do !'min_inst' not initialized? call update_accum_field ('TDM5', rbufslp, nstep) call extract_accum_field ('TDM5', a5tmin, nstep) so rbufslp will be spval (=1.e36) and then becomes SHR_CONST_TKFRZ due to rbufslp(p) > 1.e30_r8. Therefore, I need to: (a) initialize t_ref2m_max, t_ref2m_min, t_ref2m_max_inst, and t_ref2m_min_inst as CLM4.5 above; (b) remove them from the restart file. Modified process_cn.F90, catch_types.F90, clsm_ensdrv_init_routines.F90 Sarith modified mk_LDASsaRestarts.F90 and the dummy file for generating catchcn_internal_rst. (2) Checked the modified crop years in CNPhenologyMod. ========== 20170630: 1. Continued going through the To Do list in /discover/nobackup/fzeng/clm4-to-clm4.5/notes/notes: (1) In subroutine Photosynthesis, "use CNAllocationMod, only : CNAllocation_Carbon_only" was commented out. Put it back. (2) Double check this below in subroutine decomp_rate_constants in CNDecompCascadeMod_BGC.F90: nlev_soildecomp_standard=1 ! 5 originally in CLM4.5, changed to 1 for now because we have only 1 hydrologically active soil layer (1m thick), fzeng, 10 Mar 2017 (3) Why do we have grid/tile level, instead of column level, input data for sfmc to CN_Driver (see process_cn.F90), while column level wf (wf = sfmc/poros) is needed in CNFireMod? Can we change it to column level like rzmm? Will ask Randy this afternoon. (4) In process_cn.F90: fsnow is pft-level, while asnow is grid-level. (5) L828-853 in SurfaceAlbedoMod.F90 considered? No need, because in compute_rc.F90 we only do calculations on vegetated areas. (6) Checked nrad in compute_rc.F90: ok for nlevcan==1. Need modification if nlevcan>1. In SurfaceAlbedoMod.F90: ! Sun/shade big leaf code uses only one layer (nrad = ncan = 1), triggered by ! nlevcan = 1 Since our current model is big leaf model, use nlevcan = 1 for now and thus nrad = ncan = 1. Explore to improve this to multiple canopy layers later. Yes. In /terra/fzeng/clm/cesm1_2_2/models/lnd/clm/src/clm4_5/biogeophys/BiogeophysRestMod.F90, nrad is set to nlevcan (which is 1 in clm_varpar) if can't find nrad in the restart file. (7) Check the calculation of albgrd and albgri in process_cn.F90. (8) Is it correct to add "nlevcan" to clm_varpar (=1)? Yes. That's what's done in CLM4.5 itself. (9) Double checked: tlai_z in subroutine Photosynthesis is actually elai. In SurfaceAlbedoMod, if (nlevcan == 1) then tlai_z(p,1) = elai(p), and tlai_z => pps%tlai_z. Asked Randy about: calculation of tsoi17 zsoi, dzsoi, zisoi and dzsoi_decomp calculation of bd and watsat why grid-level sfmc is used in the CN model which needs column-level sfmc? ========== 20170703: 1. Continued going through the To Do list in /discover/nobackup/fzeng/clm4-to-clm4.5/notes/notes: (1) About why grid-level sfmc is used in the CN model which needs column-level sfmc. Randy's reply on 20170630: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Fanwei, I can see why Greg didn't to surface moisture in the same way as root zone moisture, when assigning values to the different zones. I think, though, that you can produce the answers you want from the following code, which was taken from Subroutine Partition. I annotated it a little with red comments. It's late; hopefully I'm not screwing this up! -- Randy SWSRF1(N)=1. ! = SM1 from the code in your powerpoint (i.e., start with this value) SWSRF2(N)=AMIN1(1., AMAX1(0.01, RZEQYI)) ! = SM2 SWSRF4(N)=AMIN1(1., AMAX1(0.01, WILT)) ! = SM4 !**** EXTRAPOLATION OF THE SURFACE WETNESSES ! 1st step: surface wetness in the unstressed fraction without considering ! the surface excess; we just assume an equilibrium profile from ! the middle of the root zone to the surface. SWSRF2(N)=((SWSRF2(N)**(-BEE(N))) - (.5/PSIS(N)))**(-1./BEE(N)) ! You are already using bee in that code; PSIS should be available too. SWSRF4(N)=((SWSRF4(N)**(-BEE(N))) - (.5/PSIS(N)))**(-1./BEE(N)) SWSRF2(N)=SWSRF2(N)+SRFEXC(N)/(surflay*poros(n)*(1.-ar1(n))+1.e-20) ! Surflay = 50 (kg/m2 of water); the same SRFEXC is used for both SWSRF2 and SWSRF4. SWSRF2(N)=AMIN1(1., AMAX1(1.E-5, SWSRF2(N))) swsrf4(n)=swsrf4(n)+srfexc(n)/(surflay*poros(n)*(1.-ar1(n))+1.e-20) SWSRF4(N)=AMIN1(1., AMAX1(1.E-5, SWSRF4(N))) ! The values [1, SWSRF2, SWSRF4] are then the surface layer equivalents to [SM1, SM2, SM4] ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Read subroutine PARTITION in GEOSland_GridComp/Shared/lsm_rouitnes.F90. WILT is WPWET. However, computation of RZEQYI varies for different situations, e.g. CATDEF .le. CDCR1 etc. See below. L533-573 lsm_rouitnes.F90: !**** SITUATION 1: CATDEF LE CDCR1 IF(CATDEF(N) .LE. CDCR1(N)) THEN RZEQY=RZEQX+WRZ RZEQYI=RZEQXI+WRZ WMNEW=WMIN+WRZ ARG3=AMAX1(-40., AMIN1(40., -AX*(1.-WMNEW))) EXPARG3=EXP(ARG3) AREA1=(1.+AX-AX*WMIN)*EXPARG1 AREA2=(1.+AX-AX*WMNEW)*EXPARG3 IF(WMNEW .GE. WILT) THEN AR1(N)=AR1(N)+AREA2-AREA1 AR2(N)=1.-AR1(N) AR4(N)=0. ENDIF IF(WMNEW .LT. WILT) THEN AREA3=(1.+AX-AX*WILT)*EXPARG2 AR1(N)=AR1(N)+AREA3-AREA1 AR2(N)=1.-AR1(N) FACTOR=(RZEQX+WRZ-WILT)/(RZEQW-WILT) AR1(N)=AR1(N)*FACTOR AR2(N)=AR2(N)*FACTOR AR4(N)=1.-FACTOR ENDIF ENDIF !**** SITUATION 2: CATDEF GT CDCR1 IF(CATDEF(N) .GT. CDCR1(N)) THEN FACTOR=(CDCR2(N)-CATDEF(N))/(CDCR2(N)-CDCR1(N)) RZEQY=WILT+(RZEQX-WILT)*FACTOR+WRZ RZEQYI=WILT+(RZEQXI-WILT)*FACTOR+WRZ IF(RZEQY .LT. WILT) THEN IF(RZEQY .LT. WILT-.001) THEN !rr WRITE(*,*) 'RZEXC WAY TOO LOW! N=',N,' RZEQY=',RZEQY !rr WRITE(*,*) 'SRFEXC=',SRFEXC(N),'RZEXC=',RZEXC(N), !rr & 'CATDEF=',CATDEF(N) ! ELSE ! WRITE(*,*) 'RZEXC TOO LOW N=',N ENDIF RZEQY=WILT RZEQYI=WILT ENDIF lsm_routines.F90: L437-681: SUBROUTINE PARTITION L1374-1597: subroutine catch_calc_soil_moist L1528: subroutine catch_calc_soil_moist calls partition process_cn.F90: L676: call catch_calc_soil_moist Can we add swsrf1, swsrf2 and swsrf4 to the argument list of subroutine catch_calc_soil_moist as optional output variables? Or any other idea? ========== 20170705: 1. Continued going through the To Do list in /discover/nobackup/fzeng/clm4-to-clm4.5/notes/notes: (1) Asked Sarith about ashift (set in subroutine BASE in lsm_routines.F90 and process_CN.F90) and the soil layer depths for soil temperature. The soil layer depths for soil temperature are in lsm_routines.F90 and they are public, so I can use them in CN_Driver. Sarith suggests to just remove ashift, and Randy thinks that's a good idea. Improved the calculation of tsoi17 in CN_DriverMod.F90. (2) About zsoi, dzsoi and zisoi used in subroutine decomp_vertprofiles in CNVerticalProfileMod.F90: Randy said let's set VERTSOILC asdide (i.e. turn it off) for now. So it's ok to set these in clm_varcon.F90: zsoi(1) = 0.5 dzsoi(1) = 1. zisoi(0) = 0. zisoi(1) = 1. dzsoi_decomp(1) = dzsoi(1) They are not used anyways if VERTSOILC is off. (3) Continued working on why grid-level sfmc is used in the CN model which needs column-level sfmc above. ========== 20170706: 1. Continued going through the To Do list in /discover/nobackup/fzeng/clm4-to-clm4.5/notes/notes: (1) For use of column level sfmc as % of WHC in the CN code, modified: subroutine catch_calc_soil_moist in lsm_routines.F90: added SWSRF1, SWSRF2 and SWSRF4 as optional output variables. process_cn.F90: added sfm and sfmm following rzm and rzmm catch_types.F90: added sfmm following rzmm clsm_ensdrv_init_routines.F90: added sfmm following rzmm CN_DriverMod.F90: added sfm and removed sfmc With Sarith's help, modified mk_LDASsaRestarts.F90: added sfmm and removed sfmcm (2) alphapsn is not used in subroutine Photosynthesis at all. Why is it there? Fine. It may be used for the isotope code in the future. (3) Should we allow btran to be modified in subroutine Photosynthesis? Will ask Randy. (4) Double checked the calculation of vcmaxcint in compute_rc.F90. (5) (Greg) find out where and how downreg is calculated: see CNAllocationMod.F90. (6) In compute_rc.F90: parabs(n) = parabs(n) + (pardir(n)*fabd + pardif(n)*fabi)*fveg(n,nv)*wl ! save absorbed PAR for FPAR calculation Should I remove "*wl"? Yes. CLM4.5 uses both leaf and stem. See CLM4.5 SurfaceRadiationMod.F90: ! Absorbed PAR profile through canopy ! If sun/shade big leaf code, nrad=1 and fluxes from SurfaceAlbedo ! are canopy integrated so that layer values equal big leaf values. if (ib == 1) then do iv = 1, nrad(p) parsun_z(p,iv) = forc_solad(g,ib)*fabd_sun_z(p,iv) + forc_solai(g,ib)*fabi_sun_z(p,iv) parsha_z(p,iv) = forc_solad(g,ib)*fabd_sha_z(p,iv) + forc_solai(g,ib)*fabi_sha_z(p,iv) end do end if (7) Double checked the radiation calculation in compute_rc.F90. ========== 20170707: 1. Continued going through the To Do list in /discover/nobackup/fzeng/clm4-to-clm4.5/notes/notes: (1) Added description to variables in compute_rc.F90 (2) Double checked CLM4.5 SurfaceRadiationMod.F90, SurfaceAlbedoMod.F90 and CanopyFluxesMod.F90, and didn't find anything needed to add to compute_rc.F90. (3) (Greg) Looked at CanopyFluxesMod.F90 on how CLM4.5 combines resistances as conductance ========== 20170710: 1. Continued going through the To Do list in /discover/nobackup/fzeng/clm4-to-clm4.5/notes/notes: (1) (Greg) Looked at CanopyFluxesMod.F90 on how CLM4.5 combines resistances as conductance: what we are doing now is correct. (2) Is it correct to interpolate lightening frequency to each time step following what we do for grn? Yes. This is the best to do now. (3) ashift: remove it from process_cn.F90, lsm_routines.F90 and GEOS_CatchCNGridComp.F90 (4) In subroutine Photosynthesis, should we allow btran to be modified within this routine? No. Make a note in the code about what btran in the output is exactly. 2. Compiled the code. 3. Created restart file for the 0.5D tile system. ========== 20170711: 1. Added print statements to process_cn.F90 to print out the values of accumulated fields for verification. Compiled the code. 2. Test run on single point: minlon = -77.5 ! min longitude maxlon = -77. ! max longitude minlat = 37. ! min latitude maxlat = 37.5 ! max latitude /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5 > ls -l build lrwxrwxrwx 1 fzeng g0620 92 2017-06-07 14:06 build -> /gpfsm/dnb31/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/exec/clm4.5/ /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5 > ls -l build/Linux/bin/LDASsaCN_mpi.x -rwxr-xr-x 1 fzeng g0620 75706050 2017-07-11 10:44 build/Linux/bin/LDASsaCN_mpi.x* /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/input/restart > ls -l /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts/0.5D/OutData2/catchcn_internal_clm45 -rw-r--r-- 1 fzeng g0620 406279770 2017-07-10 16:47 /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/restarts/0.5D/OutData2/catchcn_internal_clm45 (1) ssh -XY discover-sp3 interactive.py -A sp3 -n 96 -a g0620 -X --debug cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/run /bin/rm ../output/global/rc_out/* ./lenkf.0.i > log Crashed. forrtl: severe (174): SIGSEGV, segmentation fault occurred Image PC Routine Line Source libirc.so 00002AAAAB58E961 Unknown Unknown Unknown libirc.so 00002AAAAB58D0B7 Unknown Unknown Unknown libmpi_usempif08. 00002AAAAC2E1362 Unknown Unknown Unknown libmpi_usempif08. 00002AAAAC2E11B6 Unknown Unknown Unknown libmpi_usempif08. 00002AAAAC2CFA0C Unknown Unknown Unknown libmpi_usempif08. 00002AAAAC2B0BA8 Unknown Unknown Unknown libc.so.6 00002AAAABA14910 Unknown Unknown Unknown LDASsaCN_mpi.x 0000000001C53364 lsm_routines_mp_p 623 lsm_routines.F90 L623 SUBROUTINE PARTITION in lsm_routines.F90: SWSRF2(N)=((SWSRF2(N)**(-BEE(N))) - (.5/PSIS(N)))**(-1./BEE(N)) Checked the restart file. BEE and PSIS look correct there and no 0 is found. (2) Compiled the code with BOPT=g. "gmake clean" before compile. /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3 > /bin/cp -pr Linux/ exec/clm4.5_debug/. /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5 > ls -l build lrwxrwxrwx 1 fzeng g0620 104 2017-06-07 17:12 build -> /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/exec/clm4.5_debug/ Followed notes in /discover/nobackup/fzeng/notes/debug_procedures: cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/run setenv ESMADIR /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3 source $ESMADIR/src/g5_modules /bin/rm ../output/global/rc_out/* module load tool/allinea-tools-7.0.4 ddt /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/build/Linux/bin/LDASsaCN_mpi.x & (Note: the order is important. If do "module load tool/allinea-tools-6.1.0" before "source ../build/Linux/bin/g5_modules", the latter will unload the module.) Copy the part of lenkf.1.i from work_path to the end (no need to be in one line) to arguments. Copy this part of lenkf.1.i to Environment variables: ~~~~~~~~~~~~~~~~~~~~~~ setenv MKL_CBWR SSE4_2 # ensure zero-diff across archs setenv MV2_ON_DEMAND_THRESHOLD 8192 # MVAPICH2 setenv MYNAME `finger $USER | cut -d':' -f3 | head -1` ~~~~~~~~~~~~~~~~~~~~~~ click "Run" click the play button on the upper left to make it run until it stops. Memory error detected in lsm​_routines::partition ​(lsm​_routines.F90:612)​: null pointer dereference or unaligned memory access. (3) Repeated (2), made sure "gmake clean" before "gmake install BOPT=g". /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/run > ls -l ../build/Linux/bin/LDASsaCN_mpi.x -rwxr-xr-x 1 fzeng g0620 99506404 2017-07-11 12:39 ../build/Linux/bin/LDASsaCN_mpi.x* Gave the same error message. (3) Sarith suggests to make these changes. These changes resolved the above issue. in subroutine catch_calc_soil_moist (lsm_routines.F90) change Line 1447 to : real, dimension(NTILES), intent(out), optional :: swsrf1out, swsrf2out, swsrf4out change Line 1459 to : real, dimension(NTILES) :: srfmn, srfmx, swsrf1, swsrf2, swsrf4, rzi And around Line1542 after call partition line put these 3 lines if(present(swsrf1out)) swsrf1out = swsrf1 if(present(swsrf2out)) swsrf2out = swsrf2 if(present(swsrf4out)) swsrf4out = swsrf4 and change Line 1381 to swsrf1out, swsrf2out, swsrf4out) And in process_cn.F90 when you call catch_calc_moist swsrf1out = swsrf1, swsrf2out=swsrf2, swsrf4out= swsrf4) /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/run > ls -l ../build/Linux/bin/LDASsaCN_mpi.x -rwxr-xr-x 1 fzeng g0620 99514956 2017-07-11 13:16 ../build/Linux/bin/LDASsaCN_mpi.x* (4) The values of t_ref2m_max and t_ref2m_min don't look correct. Found that it's because "end_cd" is always true (T), which is wrong!! In process_cn.F90: end_cd = is_end_curr_day() Changed it to if (date_time%hour == 0 .and. date_time%min == 0 .and. date_time%sec == 0) then end_cd = .true. else end_cd = .false. endif 3. t_ref2m_max and t_ref2m_min are 1.e36 except for the last time step of each day. Double checked them in the CLM4.5 code. t_ref2m_max and t_ref2m_min are listed as variables in the restart file, see BiogeophysRestMod.F90 subroutine BiogeophysRest. Subroutine BiogeophysRest is called from subroutine restFile_read in restFileMod.F90, and subroutine restFile_read is called from subroutine initialize2 in clm_initializeMod.F90 (L493). However, in L554 subroutine initialize2 calls initAccClmtype, and there is this block in subroutine initAccClmtype: ! Initialize 2m ref temperature max and min values if (nsrest == nsrStartup) then ! Why not restart&branch? These vars are not in clmr. do p = begp,endp t_ref2m_max(p) = spval t_ref2m_min(p) = spval t_ref2m_max_inst(p) = -spval t_ref2m_min_inst(p) = spval t_ref2m_max_u(p) = spval t_ref2m_max_r(p) = spval t_ref2m_min_u(p) = spval t_ref2m_min_r(p) = spval t_ref2m_max_inst_u(p) = -spval t_ref2m_max_inst_r(p) = -spval t_ref2m_min_inst_u(p) = spval t_ref2m_min_inst_r(p) = spval end do end if nsrest and nsrStartup are set in clm_varctl.F90. integer, parameter, private :: iundef = -9999999 integer, public :: nsrest = iundef ! Type of run integer, public, parameter :: nsrStartup = 0 ! Startup from initial conditions Since nsrest /= nsrStartup, so the above block is not executed. So is there a bug in the following block of code in process_cn.F90? if(first) then ! see call initAccClmtype in CLM4.5 clm_initializeMod.F90 subroutine initialize2 allocate( t_ref2m_max(N_cat) ) allocate( t_ref2m_min(N_cat) ) allocate( t_ref2m_max_inst(N_cat) ) allocate( t_ref2m_min_inst(N_cat) ) t_ref2m_max(:) = spval t_ref2m_min(:) = spval t_ref2m_max_inst(:) = -spval t_ref2m_min_inst(:) = spval Probably not, because of this code below (in subroutine updateAccFlds in accFldsMod.F90, also in process_cn.F90): call update_accum_field ('TREFAV', t_ref2m, nstep) call extract_accum_field ('TREFAV', rbufslp, nstep) end_cd = is_end_curr_day() do p = begp,endp if (rbufslp(p) /= spval) then t_ref2m_max_inst(p) = max(rbufslp(p), t_ref2m_max_inst(p)) t_ref2m_min_inst(p) = min(rbufslp(p), t_ref2m_min_inst(p)) endif if (end_cd) then t_ref2m_max(p) = t_ref2m_max_inst(p) t_ref2m_min(p) = t_ref2m_min_inst(p) t_ref2m_max_inst(p) = -spval t_ref2m_min_inst(p) = spval else if (secs == int(dtime)) then t_ref2m_max(p) = spval t_ref2m_min(p) = spval endif end do Our Catchment-CN simulations alway start from the beginning of the day and end at the end of the day. If we write out t_ref2m_max, t_ref2m_min, t_ref2m_max_inst, and t_ref2m_min_inst to the restart file, because of the if (end_cd)" branch, the t_ref2m_max_inst and t_ref2m_min_inst will be -spval and spval, respectively, in the restart file. t_ref2m_max and t_ref2m_min will have valid values in the restart file, but as soon as the simulation starts, t_ref2m_max and t_ref2m_min will become spval at the first time step. the if(first) block above is doing the same thing. t_ref2m_max and t_ref2m_min are used in subroutine vernalization in CNPhenologyMod.F90. However, subroutine vernalization is only called when t_ref2m_min(p) < 1.e30_r8, i.e. it's only called at the end of each day when t_ref2m_max and t_ref2m_min have valid values, so no need to worry about t_ref2m_max = spval and t_ref2m_min = spval in the rest time of the day. NOTE: if Catchment-CN simulations do not end at the end of the day and start at the beginning of the day, t_ref2m_max, t_ref2m_min, t_ref2m_max_inst, and t_ref2m_min_inst have to be in the restart file. This is dangerous if the users do not know about this. Add t_ref2m_max, t_ref2m_min, t_ref2m_max_inst, and t_ref2m_min_inst back to the restart file. Files that need modificaitons: mk_LDASsaRestarts.F90: done process_cn.F90: done catch_types.F90: done clsm_ensdrv_init_routines.F90: done tomorrow: compile and test run ========== 20170712: 1. Compiled the code with BOPT=g. /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3 > /bin/cp -pr Linux/ exec/clm4.5_debug/. 2. Sarith helped change the dummy restart file. Ran mk_LDASsaRestarts.F90 to create the new restart file which includes T2MMAX, T2MMIN, T2MMAXI and T2MMINI. 3. Checked the time stemps of the executable and the restart file. Look correct. Did test runs. cd /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/run setenv ESMADIR /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3 source $ESMADIR/src/g5_modules /bin/rm ../output/global/rc_out/* module load tool/allinea-tools-7.0.4 ddt /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/build/Linux/bin/LDASsaCN_mpi.x & (Note: the order is important. If do "module load tool/allinea-tools-6.1.0" before "source ../build/Linux/bin/g5_modules", the latter will unload the module.) Copy the part of lenkf.1.i from work_path to the end (no need to be in one line) to arguments. Copy this part of lenkf.1.i to Environment variables: ~~~~~~~~~~~~~~~~~~~~~~ setenv MKL_CBWR SSE4_2 # ensure zero-diff across archs setenv MV2_ON_DEMAND_THRESHOLD 8192 # MVAPICH2 setenv MYNAME `finger $USER | cut -d':' -f3 | head -1` ~~~~~~~~~~~~~~~~~~~~~~ click "Run" click the play button on the upper left to make it run until it stops. (1) Reading restart file ../input/restart/catchcn_internal_rst *** glibc detected *** /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/build/Linux/bin/LDASsaCN_mpi.x: free(): invalid pointer: 0x0000000006565990 *** using catchcn_internal_rst file for all ensemble members! ======= Backtrace: ========= /lib64/libc.so.6(+0x790c8)[0x2aaaaba5b0c8] /lib64/libc.so.6(cfree+0x6c)[0x2aaaaba6010c] /usr/local/other/SLES11.3/openmpi/1.10.0/intel-15.0.3.187/lib/libmpi_usempif08.so.11(for_deallocate+0xed)[0x2aaaac2bcb8d] /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/build/Linux/bin/LDASsaCN_mpi.x[0x84ef9a] /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/build/Linux/bin/LDASsaCN_mpi.x[0x4c282e] /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/build/Linux/bin/LDASsaCN_mpi.x[0x49d88e] /lib64/libc.so.6(__libc_start_main+0xe6)[0x2aaaaba00c36] /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/m0001hdeg_clm4.5/build/Linux/bin/LDASsaCN_mpi.x[0x49d799] Is it the values -1e36 and 1e36 used for T2MMAX, T2MMIN, T2MMAXI and T2MMINI in the restart file are longer than 32-bit (Catchment-CN is always 32-bit)? (2) Tried using -500 and 500 instead of -1e36 and 1e36. Got the same error message. Checked the process_cn.F90, catch_types.F90 and clsm_ensdrv_init_routines.F90 but didn't see anything wrong with my modifications. (3) Recompiled. "gmake clean" and "gmake install BOPT=g" in (1) GEOSland_GridComp/Shared, (2) GEOSland_GridComp/GEOScatchCN_GridComp, (3) GEOSlana_GridComp, and then (4) Applications/LDAS_App. /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3 > /bin/cp -pr Linux/ exec/clm4.5_debug/. Made sure that the executable is correct (i.e. has the right time stamp). Use the restart file with 500 and -500 for T2MMAX, T2MMIN, T2MMAXI and T2MMINI. Works!! And it also works for the restart file with -1e36 and 1e36 for T2MMAX, T2MMIN, T2MMAXI and T2MMINI! So it seems that it's because I didn't do "gmake clean" before I "gmake install BOPT=g". 4. Automate the run for spin up. First recompiled the code by doing "gmake clean" and "gmake install" in (1) GEOSland_GridComp/Shared, (2) GEOSland_GridComp/GEOScatchCN_GridComp, (3) GEOSlana_GridComp, and then (4) Applications/LDAS_App. The option "BOPT=g" used earlier will slow down the simulation. ========== 20170713: 1. Decided to set init_gdd20 = .false. in CNPhenologyMod.F90 and init_accum = .false. in process_cn.F90. 2. Continued working on automating the spinup run. Proposed steps to automate spinup run: (1) modify the .bat, .exe and driver_*.nml file, and edit_lenkf.csh as well (2) run ldsetup: /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/LDASsa_m3-16_0_p2_CatchCatchCN_for_MERRA3/exec/clm4.5/Linux/bin copy edit_lenkf.csh and run it to add submit_next_batch to the end of the last lenkf.*.j (3) edit submit_next_batch (4) submit lenkf.0.j ========== 20170714: 1. Continued working on automating the spinup run. Find a way to put the restart file for the beginning of each cycle. ========== 20170717: 1. Continued working on automating the spinup run. Modified ldsetup so that there is no need to modify restart_path in lenkf.0.j, which can't be automated. Added these two lines if iseg==0: # the first job segment myOptExeInp['restart_path'] = '../input/restart/' right after this block: if iseg>0: # multiple job segments myOptExeInp['restart'] = '.true.' myOptExeInp['restart_id'] = self.rqdExeInp['exp_id'] myOptExeInp['restart_domain'] = self.rqdExeInp['exp_domain'] myOptExeInp['restart_path'] = os.path.relpath(self.wrkdir, self.rundir) ========== 20170718: 1. Continued working on automating the spinup run. Modified Sarith's edit_lenkf.csh. Worked on creating ~/Catchment/CLM4.5/submit_next_batch_point. ========== 20170719: 1. Continued working on automating the spinup run. Continued working on ~/Catchment/CLM4.5/submit_next_batch_point. Testing: /discover/nobackup/fzeng/clm4-to-clm4.5/LDAS/tests/clm4.5/clm4.5_point/run > qsub lenkf.0.j Running. Checked the accumulated field values printed out. Look correct to me. Plotted monthly time series of GPP, NPP, NEE, closs etc. Look correct. However, when this cycle is done, the command line "/home/fzeng/Catchment/CLM4.5/submit_next_batch_point clm4.5_point" in lenkf.3.j (the last job script) seemed to be skipped. Investigating ... Found that I forgot to "chmod 755 submit_next_batch_point". Fixed one bug in submit_next_batch_point. Now it seems that the problem is in "source $ESMADIR/src/g5_modules". Error message: ------------------------------------------ started time loop at 20170719, 165256.582 ended successfully at 20170719, 165256.778 ok NAME g5_modules - Script to handle BASEDIR and library module definitions DESCRIPTION This script provides a single location for storing the BASEDIR and library module name information. Previous to this script, this information was coded in multiple files. This script will set the BASEDIR environment variable to its proper value, add the BASEDIR lib directory to LD_LIBRARY_PATH (if necessary), and will load library modules when sourced. If the script is called with "basedir", "modules", "modinit", or "loadmodules", then it will echo the values to standard output without modifying the environment. The "modinit" and "loadmodules" options are primarily for use with the g5_modules_perl_wrapper script. SYNOPSIS source g5_modules or g5_modules