import netCDF4 import sys import os import xarray as xr import numpy as np from netCDF4 import Dataset as NetCDFFile from netCDF4 import num2date,date2num import datetime import pandas as pd import pylab as plt import matplotlib.cm as cm from matplotlib import colors from datetime import datetime def shr_orb_params(iyear_AD,log_print): #------------------------------------------------------------------------------- # # Calculate earths orbital parameters using Dave Threshers formula which # came from Berger, Andre. 1978 "A Simple Algorithm to Compute Long-Term # Variations of Daily Insolation". Contribution 18, Institute of Astronomy # and Geophysics, Universite Catholique de Louvain, Louvain-la-Neuve, Belgium #------------------------------Code history------------------------------------- #constant used by E3SM SHR_ORB_ECCEN_MIN = np.float64(0.0) # min value for eccen SHR_ORB_ECCEN_MAX = np.float64(0.1) # max value for eccen SHR_ORB_OBLIQ_MIN = np.float64(-90.0) # min value for obliq SHR_ORB_OBLIQ_MAX = np.float64(90.0) # max value for obliq SHR_ORB_MVELP_MIN = np.float64(0.0) # min value for mvelp SHR_ORB_MVELP_MAX = np.float64(360.0) # max value for mvelp SHR_ORB_UNDEF_REAL = np.float64(1.e36) # max value for mvelp SHR_ORB_UNDEF_INT = np.int64(2000000000) subname = '(shr_orb_params)' pi = np.pi # pi parameter poblen = 47 # of elements in series wrt obliquity pecclen = 19 # of elements in series wrt eccentricity pmvelen = 78 # of elements in series wrt vernal equinox psecdeg = 1.0/3600.0 # arc sec to deg conversion degrad = np.pi/180. # degree to radians raddeg = 180./np.pi # radians to degree # Cosine series data for computation of obliquity: amplitude (arc seconds), # rate (arc seconds/year), phase (degrees). # size : poblen obamp = [ -2462.2214466, -857.3232075, -629.3231835, -414.2804924, -311.7632587, 308.9408604, -162.5533601, -116.1077911, 101.1189923, -67.6856209, 24.9079067, 22.5811241, -21.1648355, -15.6549876, 15.3936813, 14.6660938, -11.7273029, 10.2742696, 6.4914588, 5.8539148, -5.4872205, -5.4290191, 5.1609570, 5.0786314, -4.0735782, 3.7227167, 3.3971932, -2.8347004, -2.6550721, -2.5717867, -2.4712188, 2.4625410, 2.2464112, -2.0755511, -1.9713669, -1.8813061, -1.8468785, 1.8186742, 1.7601888, -1.5428851, 1.4738838, -1.4593669, 1.4192259, -1.1818980, 1.1756474, -1.1316126, 1.0896928 ] obrate = [31.609974, 32.620504, 24.172203, 31.983787, 44.828336, 30.973257, 43.668246, 32.246691, 30.599444, 42.681324, 43.836462, 47.439436, 63.219948, 64.230478, 1.010530, 7.437771, 55.782177, 0.373813, 13.218362, 62.583231, 63.593761, 76.438310, 45.815258, 8.448301, 56.792707, 49.747842, 12.058272, 75.278220, 65.241008, 64.604291, 1.647247, 7.811584, 12.207832, 63.856665, 56.155990, 77.448840, 6.801054, 62.209418, 20.656133, 48.344406, 55.145460, 69.000539, 11.071350, 74.291298, 11.047742, 0.636717, 12.844549 ] obphas = [251.9025, 280.8325, 128.3057, 292.7252, 15.3747, 263.7951, 308.4258, 240.0099, 222.9725, 268.7809, 316.7998, 319.6024, 143.8050, 172.7351, 28.9300, 123.5968, 20.2082, 40.8226, 123.4722, 155.6977, 184.6277, 267.2772, 55.0196, 152.5268, 49.1382, 204.6609, 56.5233, 200.3284, 201.6651, 213.5577, 17.0374, 164.4194, 94.5422, 131.9124, 61.0309, 296.2073, 135.4894, 114.8750, 247.0691, 256.6114, 32.1008, 143.6804, 16.8784, 160.6835, 27.5932, 348.1074, 82.6496 ] #Cosine/sine series data for computation of eccentricity and fixed vernal #equinox longitude of perihelion (fvelp): amplitude, #rate (arc seconds/year), phase (degrees). #size: pecclen ecamp = [ 0.01860798, 0.01627522, -0.01300660, 0.00988829, -0.00336700, 0.00333077, -0.00235400, 0.00140015, 0.00100700, 0.00085700, 0.00064990, 0.00059900, 0.00037800, -0.00033700, 0.00027600, 0.00018200, -0.00017400, -0.00012400, 0.00001250] ecrate = [4.2072050, 7.3460910, 17.8572630, 17.2205460, 16.8467330, 5.1990790, 18.2310760, 26.2167580, 6.3591690, 16.2100160, 3.0651810, 16.5838290, 18.4939800, 6.1909530, 18.8677930, 17.4255670, 6.1860010, 18.4174410, 0.6678630] ecphas = [28.620089, 193.788772, 308.307024, 320.199637, 279.376984, 87.195000, 349.129677, 128.443387, 154.143880, 291.269597, 114.860583, 332.092251, 296.414411, 145.769910, 337.237063, 152.092288, 126.839891, 210.667199, 72.108838 ] # Sine series data for computation of moving vernal equinox longitude of # perihelion: amplitude (arc seconds), rate (arc sec/year), phase (degrees). # size: pmvelen mvamp = [ 7391.0225890, 2555.1526947, 2022.7629188, -1973.6517951, 1240.2321818, 953.8679112, -931.7537108, 872.3795383, 606.3544732, -496.0274038, 456.9608039, 346.9462320, -305.8412902, 249.6173246, -199.1027200, 191.0560889, -175.2936572, 165.9068833, 161.1285917, 139.7878093, -133.5228399, 117.0673811, 104.6907281, 95.3227476, 86.7824524, 86.0857729, 70.5893698, -69.9719343, -62.5817473, 61.5450059, -57.9364011, 57.1899832, -57.0236109, -54.2119253, 53.2834147, 52.1223575, -49.0059908, -48.3118757, -45.4191685, -42.2357920, -34.7971099, 34.4623613, -33.8356643, 33.6689362, -31.2521586, -30.8798701, 28.4640769, -27.1960802, 27.0860736, -26.3437456, 24.7253740, 24.6732126, 24.4272733, 24.0127327, 21.7150294, -21.5375347, 18.1148363, -16.9603104, -16.1765215, 15.5567653, 15.4846529, 15.2150632, 14.5047426, -14.3873316, 13.1351419, 12.8776311, 11.9867234, 11.9385578, 11.7030822, 11.6018181, -11.2617293, -10.4664199, 10.4333970, -10.2377466, 10.1934446, -10.1280191, 10.0289441, -10.0034259] mvrate = [31.609974, 32.620504, 24.172203, 0.636717, 31.983787, 3.138886, 30.973257, 44.828336, 0.991874, 0.373813, 43.668246, 32.246691, 30.599444, 2.147012, 10.511172, 42.681324, 13.650058, 0.986922, 9.874455, 13.013341, 0.262904, 0.004952, 1.142024, 63.219948, 0.205021, 2.151964, 64.230478, 43.836462, 47.439436, 1.384343, 7.437771, 18.829299, 9.500642, 0.431696, 1.160090, 55.782177, 12.639528, 1.155138, 0.168216, 1.647247, 10.884985, 5.610937, 12.658184, 1.010530, 1.983748, 14.023871, 0.560178, 1.273434, 12.021467, 62.583231, 63.593761, 76.438310, 4.280910, 13.218362, 17.818769, 8.359495, 56.792707, 8.448301, 1.978796, 8.863925, 0.186365, 8.996212, 6.771027, 45.815258, 12.002811, 75.278220, 65.241008, 18.870667, 22.009553, 64.604291, 11.498094, 0.578834, 9.237738, 49.747842, 2.147012, 1.196895, 2.133898, 0.173168 ] mvphas = [251.9025, 280.8325, 128.3057, 348.1074, 292.7252, 165.1686, 263.7951, 15.3747, 58.5749, 40.8226, 308.4258, 240.0099, 222.9725, 106.5937, 114.5182, 268.7809, 279.6869, 39.6448, 126.4108, 291.5795, 307.2848, 18.9300, 273.7596, 143.8050, 191.8927, 125.5237, 172.7351, 316.7998, 319.6024, 69.7526, 123.5968, 217.6432, 85.5882, 156.2147, 66.9489, 20.2082, 250.7568, 48.0188, 8.3739, 17.0374, 155.3409, 94.1709, 221.1120, 28.9300, 117.1498, 320.5095, 262.3602, 336.2148, 233.0046, 155.6977, 184.6277, 267.2772, 78.9281, 123.4722, 188.7132, 180.1364, 49.1382, 152.5268, 98.2198, 97.4808, 221.5376, 168.2438, 161.1199, 55.0196, 262.6495, 200.3284, 201.6651, 294.6547, 99.8233, 213.5577, 154.1631, 232.7153, 138.3034, 204.6609, 106.5938, 250.4676, 332.3345, 27.3039 ] #---------------------------------------------------------------------------- # radinp and algorithms below will need a degree to radian conversion factor if log_print: print('Calculate characteristics of the orbit:') # Check for flag to use input orbit parameters if iyear_AD == SHR_ORB_UNDEF_INT: # Check input obliq, eccen, and mvelp to ensure reasonable if obliq == SHR_ORB_UNDEF_REAL: print('{} Have to specify orbital parameters:'.format(subname)) print('Either set: iyear_AD, OR [obliq, eccen, and mvelp]:') print('iyear_AD is the year to simulate orbit for (ie. 1950): ') print('obliq, eccen, mvelp specify the orbit directly:') print('The AMIP II settings (for a 1995 orbit) are: ') print(' obliq = 23.4441') print(' eccen = 0.016715') print(' mvelp = 102.7') exit('{} ERROR: unreasonable obliq'.format(subname)) elif log_print: print('Use input orbital parameters: ') if obliq < SHR_ORB_OBLIQ_MIN or obliq > SHR_ORB_OBLIQ_MAX: print('Input obliquity unreasonable: {}'.format(obliq)) exit('{} ERROR: unreasonable obliq'.format(subname)) if eccen < SHR_ORB_ECCEN_MIN or eccen > SHR_ORB_ECCEN_MAX: print('Input eccentricity unreasonable: {}'.format(eccen)) exit('{} ERROR: unreasonable eccen'.format(subname)) if mvelp < SHR_ORB_MVELP_MIN or mvelp > SHR_ORB_MVELP_MAX: print('Input mvelp unreasonable: {}'.format(eccen)) exit('{} ERROR: unreasonable mvelp'.format(subname)) eccen2 = eccen*eccen eccen3 = eccen2*eccen #Otherwise calculate based on years before present else: if log_print: print('Calculate orbit for year: {}'.format(iyear_AD)) yb4_1950AD = 1950.0 - np.float64(iyear_AD) if np.abs(yb4_1950AD) > 1000000.0: print('orbit only valid for years+-1000000') print('Relative to 1950 AD') print('# of years before 1950: {}'.format(yb4_1950AD)) print('Year to simulate was : {}'.format(iyear_AD)) exit('{} ERROR: unreasonable year'.format(subname)) # The following calculates the earths obliquity, orbital eccentricity # (and various powers of it) and vernal equinox mean longitude of # perihelion for years in the past (future = negative of years past), # using constants (see parameter section) given in the program of: years = - yb4_1950AD #In the summations below, cosine or sine arguments, which end up in # degrees, must be converted to radians via multiplication by degrad. obsum = 0.0 for i in range(poblen): obsum = obsum + obamp[i]*psecdeg*np.cos((obrate[i]*psecdeg*years + obphas[i])*degrad) obliq = 23.320556 + obsum # Summation of cosine and sine series for computation of eccentricity # (eccen; e in Berger 1978) and fixed vernal equinox longitude of # perihelion (fvelp; pi in Berger 1978), which is used for computation # of moving vernal equinox longitude of perihelion. Convert the rates, # which are in arc seconds, into degrees via multiplication by psecdeg. cossum = 0.0 for i in range(pecclen): cossum = cossum+ecamp[i]*np.cos((ecrate[i]*psecdeg*years+ecphas[i])*degrad) sinsum = 0.0 for i in range(pecclen): sinsum = sinsum+ecamp[i]*np.sin((ecrate[i]*psecdeg*years+ecphas[i])*degrad) # Use summations to calculate eccentricity eccen2 = cossum * cossum + sinsum * sinsum eccen = np.sqrt(eccen2) eccen3 = eccen2 * eccen # A series of cases for fvelp, which is in radians. if np.abs(cossum) <= 1.0E-8: if sinsum == 0.0: fvelp = 0.0 elif sinsum < 0.0: fvelp = 1.5 * pi elif sinsum > 0.0: fvelp = 0.5 * pi elif cossum < 0.0: fvelp = np.arctan(sinsum/cossum) + pi else: # cossum > 1.0e-8 if sinsum < 0.0: fvelp = np.arctan(sinsum/cossum) + 2.0 * pi else: fvelp = np.arctan(sinsum/cossum) # Summation of sin series for computation of moving vernal equinox long # of perihelion (mvelp; omega bar in Berger 1978) in degrees. For mvelp, # first term is fvelp in degrees; second term is Berger 1978 psi bar # times years and in degrees; third term is Berger 1978 zeta; fourth # term is series summation in degrees. Convert the amplitudes and rates, # which are in arc seconds, into degrees via multiplication by psecdeg. # Series summation plus second and third terms constitute Berger 1978 # psi, which is the general precession. mvsum = 0.0 for i in range(pmvelen): mvsum = mvsum + mvamp[i]*psecdeg*np.sin((mvrate[i]*psecdeg*years \ + mvphas[i])*degrad) mvelp = fvelp/degrad + 50.439273*psecdeg*years + 3.392506 + mvsum # Cases to make sure mvelp is between 0 and 360. if mvelp < 0.0: mvelp = mvelp + 360.0 if mvelp >= 360.0: mvelp = mvelp - 360.0 # Orbit needs the obliquity in radians obliqr = obliq * degrad #180 degrees must be added to mvelp since observations are made from the # earth and the sun is considered (wrongly for the algorithm) to go around # the earth. For a more graphic explanation see Appendix B in: # # A. Berger, M. Loutre and C. Tricot. 1993. Insolation and Earth Orbital # Periods. J. of Geophysical Research 98:10,341-10,362. # # Additionally, orbit will need this value in radians. So mvelp becomes # mvelpp (mvelp plus pi) mvelpp = (mvelp + 180.0)*degrad #Set up an argument used several times in lambm0 calculation ahead. beta = np.sqrt(1.0 - eccen2) # The mean longitude at the vernal equinox (lambda m nought in Berger # 1978; in radians) is calculated from the following formula given in # Berger 1978. At the vernal equinox the true longitude (lambda in Berger # 1978) is 0. lambm0 = 2.0*((0.5*eccen + 0.125*eccen3)*(1.0 + beta)*np.sin(mvelpp) \ - 0.250*eccen2*(0.5 + beta)*np.sin(2.0*mvelpp) \ + 0.125*eccen3*(1.0/3.0 + beta)*np.sin(3.0*mvelpp)) if log_print: print('------ Computed Orbital Parameters ------') print('Eccentricity = {}'.format(eccen)) print('Obliquity (deg) = {}'.format(obliq)) print('Obliquity (rad) = {}'.format(obliqr)) print('Long of perh(deg) = {}'.format(mvelp)) print('Long of perh(rad) = {}'.format(mvelpp)) print('Long at v.e.(rad) = {}'.format(lambm0)) print('-----------------------------------------') return eccen,obliq,mvelp,obliqr,lambm0,mvelpp def shr_orb_decl(calday,eccen,mvelpp,lambm0,obliqr): #------------------------------------------------------------------------------- # Compute earth/orbit parameters using formula suggested by Duane Thresher. dayspy = 365.0 # days per year ve = 80.5 # Calday of vernal equinox pi = np.pi # pi parameter degrad = np.pi/180. # degree to radians raddeg = 180./np.pi # radians to degree # Compute eccentricity factor and solar declination using # day value where a round day (such as 213.0) refers to 0z at # Greenwich longitude. # Use formulas from Berger, Andre 1978: Long-Term Variations of Daily # Insolation and Quaternary Climatic Changes. J. of the Atmo. Sci. # 35:2362-2367. lambm = lambm0 + (calday - ve)*2.0*pi/dayspy lmm = lambm - mvelpp # The earths true longitude, in radians, is then found from # the formula in Berger 1978: sinl = np.sin(lmm) lamb = lambm + eccen*(2.0*sinl + eccen*(1.25*np.sin(2.0*lmm) \ + eccen*((13.0/12.0)*np.sin(3.0*lmm) - 0.25*sinl))) # Using the obliquity, eccentricity, moving vernal equinox longitude of # perihelion (plus), and earths true longitude, the declination (delta) # and the normalized earth/sun distance (rho in Berger 1978; actually inverse # rho will be used), and thus the eccentricity factor (eccf), can be # calculated from formulas given in Berger 1978. invrho = (1.0 + eccen*np.cos(lamb - mvelpp)) / (1.0 - eccen*eccen) # Set solar declination and eccentricity factor delta = np.arcsin(np.sin(obliqr)*np.sin(lamb)) eccf = invrho*invrho return delta, eccf def shr_orb_print( iyear_AD, eccen, obliq, mvelp ): #------------------------------------------------------------------------------- # Print out the information on the Earths input orbital characteristics #------------------------------------------------------------------------------- if iyear_AD != SHR_ORB_UNDEF_INT: if iyear_AD > 0: print('Orbital parameters calculated for year: AD '.format(iyear_AD)) else: print('Orbital parameters calculated for year: BC '.format(iyear_AD)) elif obliq != SHR_ORB_UNDEF_REAL: print('Orbital parameters: ') print('Obliquity (degree): {}'.format(obliq)) print('Eccentricity (unitless): {}'.format(eccen)) print('Long. of moving Perhelion (deg): {}'.format(mvelp)) else: print('Orbit parameters not set!') def shr_orb_avg_cosz(jday, lat, lon, declin, dt_avg): SHR_ORB_UNDEF_REAL = np.float64(1.e36) # max value for mvelp SHR_ORB_UNDEF_INT = np.int64(2000000000) pi = np.pi piover2 = np.float64(pi/2.0) twopi = np.float64(pi*2.0) #----------------------------------------------------------------------- #Compute Half-day Length #adjust latitude so that its tangent will be defined if lat == piover2: xdel = lat - np.float64(1.0e-5) elif lat == -piover2: xdel = lat + np.float64(1.0e-5) else: xdel = lat # adjust declination so that its tangent will be defined if declin == piover2: phi = declin - np.float64(1.0e-5) elif declin == -piover2: phi = declin + np.float64(1.0e-5) else: phi = declin # define the cosine of the half-day length # adjust for cases of all daylight or all night cos_h = - np.tan(xdel) * np.tan(phi) if cos_h <= -1.0: h = pi elif cos_h >= 1.0: h = 0.0 else: h = np.arccos(cos_h) #----------------------------------------------------------------------- # Define Local Time t and t + dt # adjust t to be between -pi and pi t1 = (jday - np.int(jday)) * twopi + lon - pi if t1 >= pi: t1 = t1 - twopi elif t1 < -pi: t1 = t1 + twopi dt = dt_avg / 86400.0 * twopi t2 = t1 + dt #----------------------------------------------------------------------- # Compute Cosine Solar Zenith angle # define terms needed in the cosine zenith angle equation aa = np.sin(lat) * np.sin(declin) bb = np.cos(lat) * np.cos(declin) # define the hour angle # force it to be between -h and h # consider the situation when the night period is too short if t2 >= pi and t1 <= pi and (pi - h) <= dt: tt2 = h tt1 = np.minimum(np.maximum(t1, -h) , h) tt4 = np.minimum(np.maximum(t2, twopi - h), twopi + h) tt3 = twopi - h elif t2 >= -pi and t1 <= -pi and (pi - h) <= dt: tt2 = - twopi + h tt1 = np.minimum(np.maximum(t1, -twopi - h), -twopi + h) tt4 = np.minimum(np.maximum(t2, -h) , h) tt3 = -h else: if t2 > pi: tt2 = np.minimum(np.maximum(t2 - twopi, -h), h) elif t2 < - pi: tt2 = np.minimum(np.maximum(t2 + twopi, -h), h) else: tt2 = np.minimum(np.maximum(t2 , -h), h) if t1 > pi: tt1 = np.minimum(np.maximum(t1 - twopi, -h), h) elif t1 < - pi: tt1 = np.minimum(np.maximum(t1 + twopi, -h), h) else: tt1 = np.minimum(np.maximum(t1 , -h), h) tt4 = 0.0 tt3 = 0.0 #! perform a time integration to obtain cosz if desired # output is valid over the period from t to t + dt if tt2 > tt1 or tt4 > tt3: shr_orb_avg_cosz = (aa * (tt2 - tt1) + bb * (np.sin(tt2) - np.sin(tt1))) / dt \ + (aa * (tt4 - tt3) + bb * (np.sin(tt4) - np.sin(tt3))) / dt else: shr_orb_avg_cosz = 0.0 return def main(exp,year,smdh,emdh,tunit,datadir,dataref,outpath): #----------------------------------------------------------------------------- # This script contains a post-processing calculation of the cosine # of the solar zenith angle following the same method as in # shr_orb_mod.F90 for E3SM. Assumes 365.0 days/year. #----------------------------------------------------------------------------- var_out = "COSSZA" #constant values pi = np.pi deg2rad=np.pi/180. # degree to radians rad2deg=180./np.pi # radians to degree #setup followed E3SM use_dt_avg = False l_uniform_angle = False if l_uniform_angle: uniform_angle = custom_value # value to apply uniform angle dt_avg = 0.0 # variable to calculate time average (0.0 for false ) if dt_avg != 0.0: use_dt_avg = True print('use_dt_avg:',use_dt_avg) print('uniform_angle',l_uniform_angle) #This variable overrides the behavior of shr_orb_cosz() when >=0 #this is be set by calling set_constant_zenith_angle_deg() constant_zenith_angle_deg = -1 # constant, uniform zneith angle [degrees] #select the formula used for solar zenith angle calculation l_atm_formula = True #False #select the formula used for solar zenith angle calculation l_icepack_formula = False #sanity check l_sanity_check = False #print info during calculation log_print = False #note: we use this simulation to extract the time information #note: to keep the consistenty time format used for other quantity fref = "{}_3hourly_{:04d}{}-{:04d}{}.nc".format(exp,year,smdh,year,emdh) fr = xr.open_dataset(os.path.join(dataref,fref),decode_times=True,decode_timedelta=True) lat_val = fr['lat'] lon_val = fr['lon'] tim_val = fr['time'].astype('datetime64[ns]') #.to_pandas() #.astype('datetime64[ns]') #to_pandas() print('latitude range:', lat_val[:].min(),lat_val[:].max()) print('longitude range:', lon_val[:].min(),lon_val[:].max()) print('time range:', tim_val[0],tim_val[len(tim_val)-1]) ax1 = len(tim_val) # itim2 - itim1 + 1 ax2 = len(lat_val) # ilat2 - ilat1 + 1 ax3 = len(lon_val) # ilon2 - ilon1 + 1 fr.close() #----------------------------------------------------------------------- # ... Calculate cosine of zenith angle # then cast back to angle (radians) #----------------------------------------------------------------------- rlat = lat_val.copy() * deg2rad #(to radians) rlon = lon_val.copy() * deg2rad #(to radians) print("min/max longitude (radians):",rlon.min().values,rlon.max().values) print("min/max longitude (degree):",lon_val.min().values,lon_val.max().values) print("constant_zenith_angle_deg",np.cos( constant_zenith_angle_deg * deg2rad )) xdcln = np.zeros(shape=(ax1, ax2)) xeccf = np.zeros(shape=(ax1, ax2)) xcday = np.zeros(shape=(ax1, ax2)) xclat = np.zeros(shape=(ax1, ax2)) xclon = np.zeros(shape=(ax1, ax2)) data_np = np.zeros(shape=(ax1, ax2)) eccen,obliq,mvelp,obliqr,lambm0,mvelpp = shr_orb_params(year,log_print) for it in range(len(tim_val)): tcurr = str(tim_val.to_pandas().iloc[it]) doyr = datetime.strptime(tcurr,'%Y-%m-%d %H:%M:%S').timetuple().tm_yday clhr = datetime.strptime(tcurr,'%Y-%m-%d %H:%M:%S').timetuple().tm_hour calday = doyr + np.float64(clhr) / 24.0 #print(tcurr,calday) # calculate solar decline angle declin, eccf = shr_orb_decl (calday,eccen,mvelpp,lambm0,obliqr) #print("current time:",tcurr) #print("calendar day:",calday) #print("calendar day(diff):",calday-np.floor(calday)) #print("decline angle:",declin*rad2deg) #eccen,obliq,mvelp,obliqr,lambm0,mvelpp) #print("eccentricity:",eccf) xcday[it,:] = calday xdcln[it,:] = declin xeccf[it,:] = eccf xclat[it,:] = rlat[:] xclon[it,:] = rlon[:] del(tcurr,doyr,clhr,calday,declin,eccf) print("min/max decline angle:",xdcln.min()*rad2deg,xdcln.max()*rad2deg) print("min/max eccentricity:",xeccf.min(),xeccf.max()) #start to calculate solar zenith angle: if constant_zenith_angle_deg >= 0: data_np[:,:] = np.cos( constant_zenith_angle_deg * deg2rad ) if l_uniform_angle: data_np[:,:] = np.cos(uniform_angle * deg2rad) if use_dt_avg: for it in range(len(tim_val)): for ic in range(len(lat_val)): data_np[it,ic] = shr_orb_avg_cosz(xcday[it,ic], rlat[ic], rlon[ic],xdcln[it,ic], dt_avg) elif l_atm_formula: #formula in orbit.F90 xfact = (xcday-np.floor(xcday))*2.0*pi + xclon data_np = np.sin(xclat)*np.sin(xdcln) \ - np.cos(xclat)*np.cos(xdcln) * np.cos(xfact) del(xfact) elif l_icepack_formula: #formula in ice_orbital.F90 secday = 86400.0 p5 = 0.5 sec = xcday * secday xfact = (sec / secday - p5)*2.0*pi + xclon data_np = np.sin(xclat)*np.sin(xdcln) \ + np.cos(xclat)*np.cos(xdcln)*np.cos(xfact) del(secday,p5,sec,xfact) print("min/max cossza:",data_np.min(),data_np.max()) exit() eccen,obliq,mvelp,obliqr,lambm0,mvelpp = shr_orb_params(year,log_print) for it in range(len(tim_val)): tcurr = str(tim_val.to_pandas().iloc[it]) doyr = datetime.strptime(tcurr,'%Y-%m-%d %H:%M:%S').timetuple().tm_yday clhr = datetime.strptime(tcurr,'%Y-%m-%d %H:%M:%S').timetuple().tm_hour calday = doyr + np.float64(clhr) / 24.0 #print(tcurr,calday) # calculate solar decline angle declin, eccf = shr_orb_decl (calday,eccen,mvelpp,lambm0,obliqr) print("current time:",tcurr) print("calendar day:",calday) print("calendar day(diff):",calday-np.floor(calday)) print("decline angle:",declin*rad2deg) #eccen,obliq,mvelp,obliqr,lambm0,mvelpp) print("eccentricity:",eccf) #start to calculate solar zenith angle: if constant_zenith_angle_deg >= 0: data_np[it,:] = np.cos( constant_zenith_angle_deg * deg2rad ) if l_uniform_angle: data_np[it,:] = np.cos(uniform_angle * deg2rad) if use_dt_avg: for ic in range(len(lat_val)): data_np[it,ic] = shr_orb_avg_cosz(calday, rlat[ic], rlon[ic], declin, dt_avg) else: if l_atm_formula: #formula in orbit.F90 #data_np[it,:] = np.sin(rlat[:])*np.sin(declin) \ # - np.cos(rlat[:])*np.cos(declin) \ # * np.cos((calday-np.floor(calday))*2.0*pi + rlon[:]) data_np[it,:] = np.sin(rlat[:])*np.sin(declin) \ + np.cos(rlat[:])*np.cos(declin) \ * np.cos(calday*2.0*pi + rlon[:]) print("min/max cosine solar zenith angle:",data_np[it,:].min(),data_np[it,:].max()) elif l_icepack_formula: #formula in ice_orbital.F90 secday = 86400.0 p5 = 0.5 sec = calday * secday data_np[it,:] = np.sin(rlat[:])*np.sin(declin) \ + np.cos(rlat[:])*np.cos(declin) \ * np.cos((sec/secday - p5)*2.0*pi + rlon[:]) else: exit("need to speify an option for solar zenith angle calculation") #sanity check or output data if l_sanity_check: var_ref = "LINOZ_SZA" fname = "{}_Scalar_monthly_{}.nc".format(exp,year) df = xr.open_dataset(os.path.join(datadir,fname),decode_times=True,decode_timedelta=True) data = df[var_ref][:,:] ref_time = df['time'].astype('datetime64[ns]') str_time = df['time_bnds'][:,0].astype('datetime64[ns]') end_time = df['time_bnds'][:,1].astype('datetime64[ns]') df.close() #convert diagnosed cossza to solar zenith angle data_np = np.arccos(data_np) * rad2deg #- 180.0 # create data for plot vtim1 = ref_time.values.astype('datetime64[ms]').astype('O') vtim2 = tim_val.values.astype('datetime64[ms]').astype('O') xvals = np.arange(len(lat_val)) da = xr.DataArray(data_np, dims=("time", "ncol"),coords={"time": tim_val}) dam = da.groupby("time.month").mean("time") print(data.data.min(),data.data.max()) print(dam.data.min(),dam.data.max()) print(da.data.min(),da.data.max()) #start to plot data clevs = np.arange(0,180,11) fig, axs = plt.subplots(3,1, figsize=(10,15)) #fig.suptitle('{}'.format(var_out)) cntr1 = axs[0].contourf(xvals,vtim1, data.data, clevs, cmap=plt.cm.jet) axs[0].set_title('Monthly Mean (Model Output)', loc='left') axs[0].set_xlabel("model columns (#)") cntr2 = axs[1].contourf(xvals,vtim1, dam.data, clevs, cmap=plt.cm.jet) axs[1].set_title('3hourly to Monthly Mean (offline calculation)', loc='left') axs[1].set_xlabel("model columns (#)") cntr3 = axs[2].contourf(xvals,vtim2, da.data, clevs, cmap=plt.cm.jet) axs[2].set_title('3hourly (offline calculation)', loc='left') axs[2].set_xlabel("model columns (#)") plt.subplots_adjust(bottom=0.2, right = 0.8, hspace=0.5) #, right=0.9, top=0.9) cbar_ax = fig.add_axes([0.82, 0.19, 0.02, 0.65]) #(left, bottom, width, height) cbar = fig.colorbar(cntr3, cax=cbar_ax, shrink=0.6, orientation='vertical', pad=0.1, aspect=10, extendrect=True) cbar.set_label('Solar Zenith Angle (degree)') #fig.colorbar(cntr1[0], ax=axs, orientation='horizontal', fraction=.1) figname = "eam_{}_3hourly_{:04d}.png".format("solar_zenith_angle",year) if l_icepack_formula: figname = "icepack_{}_3hourly_{:04d}.png".format("solar_zenith_angle",year) plt.savefig(os.path.join("./diag_figure",figname)) plt.close() del(vtim1,vtim2,xvals,da,dam,clevs,fig,axs,cntr1,cntr2,cntr3,cbar_ax,cbar) del(df,fname,data) else: outname = "{}_3hourly_{:04d}.npy".format(var_out,year) fout = os.path.join(outpath,outname) if os.path.exists(fout): os.remove(fout) np.save(fout, data_np, allow_pickle=True,fix_imports=True) del(outname,fout,data_np) del(ax1,ax2,ax3,fr,lat_val,lon_val,tim_val) return if __name__ == "__main__": topdir = "/global/cfs/cdirs/e3sm/www/zhan391/darpa_temporary_data_share" outdir = "/global/cfs/cdirs/e3sm/www/zhan391/SEA_CROGS" syear = 2007 eyear = 2017 tunit = "hours since 2007-01-01 00:00:00" exp = "E3SMv2_NDGUVTQ_SRF1_tau6" datadir = os.path.join(topdir,"New_Training","nc_data",exp) dataref = os.path.join(topdir,"SE_PG2","before_nudging") smdh = "01010000" emdh = "12312100" #output path outpath = os.path.join(outdir,"New_Training","E3SMv2_Scalar") if not os.path.exists(outpath): os.makedirs(outpath) for year in range(syear,eyear+1): main(exp,year,smdh,emdh,tunit,datadir,dataref,outpath)