Arm diags v4

arm_diags/arm_driver.py

@@ -4,6 +4,7 @@
 import numpy
 import cdutil
 import genutil
+import shutil
 import cdms2
 import MV2
 import glob
@@ -20,13 +21,18 @@
 from src.pdf_daily import pdf_daily_data, pdf_daily_plot # .src
 from src.convection_onset_driver import convection_onset # .src
 from src.aerosol_activation import aerosol_activation_density_plot # .src
-#from src.convection_onset_driver_todd import convection_onset
-from src.create_htmls import annual_cycle_zt_html,diurnal_cycle_zt_html,diurnal_cycle_html,seasonal_mean_table_html,annual_cycle_html,annual_cycle_aci_html,pdf_daily_html,convection_onset_html,aerosol_activation_html,diags_main_html
+from src.twolegged_metric import twolegged_metric_plot # .src
+from src.diurnal_cycle_LAcoupling import diurnal_cycle_LAcoupling_plot
 
+#from src.convection_onset_driver_todd import convection_onset
+from src.create_htmls import annual_cycle_zt_html,diurnal_cycle_zt_html,diurnal_cycle_html,seasonal_mean_table_html,annual_cycle_html,annual_cycle_aci_html,pdf_daily_html,convection_onset_html,aerosol_activation_html,twolegged_metric_html,diurnal_cycle_LAcoupling_html,diags_main_html
 def make_parameters(basic_parameter):
     #f_data = open('examples/diags_set3.json').read()
     #f_data = open('diags_all_multisites_for_cmip5.json').read()
     f_data = open('diags_all_multisites_for_cmip6.json').read()
+    #f_data = open('diags_set10_cmip6.json').read()
+    #f_data = open('diags_all_multisites_for_LAcoupling.json').read()
+    #f_data = open('diags_all_multisites_check.json').read()
     json_file = json.loads(f_data)
 
     parameters = [] 
@@ -48,9 +54,10 @@ def make_parameters(basic_parameter):
 #basic_parameter = parser.get_parameters()
 parameters = make_parameters(basic_parameter)
 
-
 case_id =  basic_parameter.case_id
 output_path = basic_parameter.output_path
+armdiags_path = basic_parameter.armdiags_path
+print('output_path: ',output_path)
 
 # Generate new case folder given case_id:
 if not os.path.exists(os.path.join(output_path)):
@@ -59,6 +66,13 @@ def make_parameters(basic_parameter):
     os.makedirs(os.path.join(output_path,'figures'))
     os.makedirs(os.path.join(output_path,'metrics'))
 
+# Copy the logo figures to the newly created html folder
+src = os.listdir(armdiags_path+'arm_diags/misc/')
+dst = output_path+'/html/'
+for ifile in range(len(src)):
+    src1 = armdiags_path+'arm_diags/misc/'+src[ifile]
+    shutil.copy(src1, dst)
+
 # Loop through diagnostic sets prespecified from diags_sets.json
 html_count = 0
 for parameter in parameters:
@@ -138,7 +152,18 @@ def make_parameters(basic_parameter):
         html_count = html_count + 1
         #except:
             #pass
-
+
+    if diags_set == 'set10_twolegged_metric':
+        twolegged_metric_plot(parameter)
+        twolegged_metric_html(parameter)
+        html_count = html_count + 1     
+
+    if diags_set == 'set11_diurnal_cycle_LAcoupling':
+        diurnal_cycle_LAcoupling_plot(parameter)
+        diurnal_cycle_LAcoupling_html(parameter)
+        html_count = html_count + 1
+
+
 #
 if html_count >= 1:
 # Creat the main html page hosting all sets of diagnostics

arm_diags/basicparameter.py

@@ -5,28 +5,39 @@
 #---------------------------------------------------------------------------------------------------------------------------
 #===========================================================================================================================
 # User defined case id
-case_id = 'V3_TestCMIP6'
+#case_id = 'output_cheng_20230705_mdtfv3_cesm'
+case_id = 'output_cheng_20240905_armdiags_v4'
+
+# User defined the ARM-Diags package path
+armdiags_path = '/Users/tao4/Documents/ARM_Infrastructure/ARM_DIAG/arm-gcm-diagnostics/'
 
 #--------------------------------------------------------------------------
 # Testing model dataset (User defined model)
 # Specify name of the test model to find the files
 test_data_set = 'testmodel'
 #specify the data starting/ending years in the testmodel file
 #default is 1979 - 2006 as in the CMIP file
-test_start_year = 1979 
+
+#NCAR: 2013, 2014
+#GFDL: 1980, 2000
+test_start_year = 1979
 test_end_year   = 2006
 
 #--------------------------------------------------------------------------
 # Set input path, where the model, observational and cmip data are located.
-base_path = '/DATA/ARM-Diag/arm-gcm-diagnostics/arm_diags/'
-
+#base_path = '/DATA/ARM-Diag/arm-gcm-diagnostics/arm_diags/'
+#base_path = '/Users/tao4/Documents/ARM_Infrastructure/ARM_DIAG/arm_diags_data_v3.1_07052023_mdtfv3_cesm/'
+#base_path = '/Users/tao4/Documents/ARM_Infrastructure/ARM_DIAG/arm_diags_data_v3.1_05192023/'
+#base_path = '/Users/tao4/Documents/ARM_Infrastructure/ARM_DIAG/arm_diags_data_v3.1_06292023_mdtfv3_gfdl/'
+base_path = '/Users/tao4/Documents/ARM_Infrastructure/ARM_DIAG/arm_diags_data_v3.1_06122023/'
 test_data_path = base_path+'testmodel'
 obs_path = base_path+'observation'
 cmip_path =  base_path+'cmip6'
 
 #--------------------------------------------------------------------------
 # Set output path, where the results will be saved
-output_path = '/DATA/ARM-Diag/Results/'+case_id
+#output_path = '/DATA/ARM-Diag/Results/'+case_id
+output_path = '/Users/tao4/Documents/ARM_Infrastructure/ARM_DIAG/'+case_id
 
 arm_filename = True
 #===========================================================================================================================

arm_diags/diags_all_multisites_for_cmip6.json

@@ -283,7 +283,22 @@
       "diags_set": "set9_aerosol_activation",
       "variables": ["cpc","ccn"],
       "sites": ["enac1"]
+                },
+          {
+      "diags_set": "set10_twolegged_metric",
+      "variables": ["Atmospheric Coupling Legs","Terrestrial Cooupling Legs"],
+      "sites": ["sgpc1"],
+      "season": ["DJF","MAM","JJA","SON"]
+                },
+          {
+      "diags_set": "set11_diurnal_cycle_LAcoupling",
+      "variables": ["SH","LH","T_srf","RH_srf","LCL","pbl"],
+      "varnames": ["Surface Sensible Heat Flux","Surface Latent Heat Flux","Surface Air Temperature","Surface Relative Humidity","Lifting Condensation Level","Planatary Boundary Layer"],
+      "units": ["W/m\u00b2","W/m\u00b2","deg C","%","meter","meter"],
+      "sites": ["sgpc1"],
+      "season": ["ANN","DJF","MAM","JJA","SON"]
                 }
+
           ]
 }
 

arm_diags/src/annual_cycle_zt.py

@@ -344,8 +344,10 @@ def annual_cycle_zt_plot(parameter):
     if test_findex == 1: 
         ct_up=np.nanmax([cl_ob,cl_p])
         tmpct=cl_p[:,::-1]-cl_ob[:,::-1]
-        ct_lo_diff=np.int(np.nanmin(tmpct)-1)
-        ct_up_diff=np.int(np.nanmax(tmpct)+1)
+
+        ct_lo_diff=int(np.nanmin(tmpct)-1)
+        ct_up_diff=int(np.nanmax(tmpct)+1)
+
     rlevel=np.arange(ct_lo,ct_up+1,0.5) #original
     drlevel=np.arange(ct_lo_diff,ct_up_diff,0.5)   #difference
     #---------------------------------------------- 

arm_diags/src/convection_onset_driver.py

@@ -86,10 +86,10 @@ def convection_onset(parameter):
             else:
                 filename = glob.glob(os.path.join(obs_path,site[:3]+'armdiags1hr' + site[3:5].upper()+'*.nc'))[0]
 
-            print(filename)
+            print('filename: ',filename)
             f_in=cdms2.open(filename)
             var=f_in(va)
-            print('var_shape',va,var.shape)
+            print('var_shape: ',va,var.shape)
             if va == 'pr':
                 precip = var
                 precip[precip<-900] = np.nan

arm_diags/src/convection_onset_statistics.py

@@ -70,13 +70,14 @@ def convection_onset_statistics(precip_threshold,cwv_max,cwv_min,bin_width,cwv,p
     #cwv_max = 85 # default = 70 (in mm)
     #cwv_min = 28 # default = 28 (in mm)
     #bin_width = 1.5 # default = 1.5
+    #print('inputs: ',cwv_max,cwv_min,bin_width)
     number_of_bins = int(np.ceil((cwv_max-cwv_min)/bin_width))
-    #print('cwv_max',cwv_max)
-    #print(number_of_bins)
+    #print('number_of_bins: ',number_of_bins)
     bin_center = np.arange((cwv_min+(bin_width/2)), (cwv_max-(bin_width/2))+bin_width, bin_width)
     #print('bin_center',bin_center,'bin_width',bin_width)
     if len(bin_center)!=number_of_bins:
         bin_center = np.arange((cwv_min+(bin_width/2)), (cwv_max-(bin_width/2)), bin_width)
+    #print('bin_center range: ',np.min(bin_center),np.max(bin_center))
 
     # Define precip threshold
     #precip_threshold = 0.5 # default 0.5 (in mm/hr)
@@ -86,8 +87,14 @@ def convection_onset_statistics(precip_threshold,cwv_max,cwv_min,bin_width,cwv,p
     precip_binned = np.empty([number_of_bins,cwv.size]) * np.nan
     precip_counts = np.zeros([number_of_bins,cwv.size])
 
-    np.warnings.filterwarnings('ignore')
+    #============================================
+    # Cheng 09/05/2024
+    # Comment out the following line.
+    # Otherwise, the figures cannot be generated
+    #==========================================
+    #np.warnings.filterwarnings('ignore')
     # Bin the data by CWV value as specified above
+    print("Bin the data by CWV value as specified above")
     for i in range (0,number_of_bins):
         tmp1 = np.where(cwv > cwv_min+(i*bin_width))
         bin_index[i,tmp1] = 1
@@ -116,6 +123,7 @@ def convection_onset_statistics(precip_threshold,cwv_max,cwv_min,bin_width,cwv,p
             #    precip_counts[i,j] = 1
             if np.isnan(precip_binned[i,j]):
                 precip_counts[i,j] = np.nan
+    #print('Check')
 
     # Create binned arrays
     hist_cwv = np.empty([number_of_bins,1]) * np.nan
@@ -175,19 +183,19 @@ def convection_onset_statistics(precip_threshold,cwv_max,cwv_min,bin_width,cwv,p
     xulim = 5*np.ceil(np.max(np.round(bin_center+bin_width/2))/5)
     xllim = 5*np.floor(np.min(np.round(bin_center-bin_width/2))/5)
     ax1.set_xlim(xllim-10,xulim+15)
-    ax1.set_ylim(0,5)
+    #ax1.set_ylim(0,5)
     ax1.set_xticks(np.arange(np.ceil(xllim/10)*10-10,np.ceil(xulim/10)*10+15,15))
     #print(np.arange(np.ceil(xllim/10)*10-10,np.ceil(xulim/10)*10+15,15))
     #print(np.ceil(xllim/10)*10)
     #print(np.ceil(xulim/10)*10)
     ulim = np.nanmax(pr_binned_mean)
     #print('max precip',ulim)
-    #ax1.set_yticks(np.arange(0,np.ceil(ulim)))
-    ax1.set_yticks(np.arange(0,5))
+    ax1.set_yticks(np.arange(0,np.ceil(ulim)))
+    #ax1.set_yticks(np.arange(0,5))
     #ax1.set_xlim(25,72)
     #ax1.set_ylim(0,6)
     #ax1.set_xticks([30,40,50,60,70])
-    #ax1.set_yticks([0,1,2,3,4,5,6])
+    ax1.set_yticks([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15])
     ax1.tick_params(labelsize=axes_fontsize)
     ax1.tick_params(axis='x', pad=10)
     error = [errorbar_precip,errorbar_precip]