;******************************************************************************************* ; Z_MPI_ESM1_2_LR_historical.ncl ;******************************************************************************************* ; ; Script to compute geopotential height in hybrid coordinates (hybrid model levels), via ; function "cz2ccm" which is the exact routine used in the former Cray CESM Processor. ; ;******************************************************************************************* load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" begin DataDir = "/media/marques/d54b17d1-d931-4891-87a4-a80844295958/data/CMIP6/MPI_ESM1.2_LR/historical/" DataDirOut = DataDir ;do A=1979,2014 ; Loops over years do A=1980,2014 ; Loops over years do m=1,12 ; Loops over months ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% INPUT FILES TO READ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ============================================= TEMPERATURE ===================================================== Ta = addfile(DataDir+"ta_6hrLev_MPI-ESM1-2-LR_historical_r1i1p1f1_gn_"+sprinti("%0.4i",A)+sprinti("%0.2i",m)+".nc", "r") ; ========================================== SPECIFIC HUMIDITY =================================================== Ha = addfile(DataDir+"hus_6hrLev_MPI-ESM1-2-LR_historical_r1i1p1f1_gn_"+sprinti("%0.4i",A)+sprinti("%0.2i",m)+".nc", "r") ; ========================================= SURFACE GEOPOTENTIAL ================================================ Zs = addfile(DataDir+"orog_fx_MPI-ESM1-2-LR_historical_r1i1p1f1_gn.nc", "r") ; ========================================== SURFACE PRESSURE =================================================== PS = addfile(DataDir+"ta_6hrLev_MPI-ESM1-2-LR_historical_r1i1p1f1_gn_"+sprinti("%0.4i",A)+sprinti("%0.2i",m)+".nc", "r") ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; --------------------------------------- Constants ---------------------------------------------------- ;p0 = 101320 ; Reference pressure in Pa. p0 = 100000 ; Reference pressure in Pa. ; The normal acceleration of gravity: ; g = 9.80616 - This value has been chosen as the most representative of the acceleration of gravity ; at sea-level at latitude 45º in accordance with the opinion of the International ; Association of Geodesy (1950). This is the value generally used by physicists. ; ; g = 9.80665 - Value adopted by WMO (1966) which agrees with that recommended by the International ; Committee of Weights and Measures (1901, 1948). g = 9.80665 Rv = 461.51 ; (WMO, 1966) Rd = 287.05 ; (WMO, 1966) ; ------------------------------------------------------------------------------------------------------ ; ############################################# LOADING VARIABLES ############################################### ; =====================> These are needed to compute geopotential height <========================= ; NOTE: The latitudes in phis, ta, hus and ps have been reversed south-to-north: phis = Zs->orog*g ; surface geopotential (m^2/sec^2). ta = Ta->ta(:,::-1,:,:) ; air temperature (K) ordered top-to-bottom. hus = Ha->hus(:,::-1,:,:) ; Specific humidity (kg/kg) ordered top-to-bottom. ps = PS->ps(:,:,:) ; surface pressures (Pa). ; /////// Assignment of named dimensions /////// phis!0 = "lat" phis!1 = "lon" ta!0 = "time" ta!1 = "lev" ta!2 = "lat" ta!3 = "lon" hus!0 = "time" hus!1 = "lev" hus!2 = "lat" hus!3 = "lon" ps!0 = "time" ps!1 = "lat" ps!2 = "lon" ;///////////////////////////////////////////////// ; =================> These are needed to compute vertical pressure velocity <====================== ap = Ta->ap ; hybrid A coefficients at layer midpoints ordered bottom-to-top. ; These are expected to be normalized to one (1.0). If not, divide by p0. b = Ta->b ; hybrid B coefficients at layer midpoints ordered bottom-to-top. ap_bnds = Ta->ap_bnds ; hybrid A coefficients at layer interfaces ordered bottom-to-top. ; If not normalized to one (1.0), divide by p0. b_bnds = Ta->b_bnds ; hybrid B coefficients at layer interfaces ordered bottom-to-top. ; ================================================><================================================ ; ---------------------- Calculate virtual temperature ---------------------------- Tv = ta*(1.+(Rv/Rd-1)*hus) ; Ordered top-to-bottom ; Assign named dimensions Tv!0 = "time" Tv!1 = "lev" Tv!2 = "lat" Tv!3 = "lon" delete(ta) delete(hus) ; --------------------------------------------------------------------------------- ; ==========================================><====================================================== ; ========================> Remaining variable to include in output file <========================== ap = Ta->ap ; hybrid A coefficients at layer midpoints. ap_bnds = Ta->ap_bnds ; hybrid A coefficients at layer interfaces. b = Ta->b ; hybrid B coefficients at layer midpoints. b_bnds = Ta->b_bnds ; hybrid B coefficients at layer interfaces. lon = Ta->lon ; Longitudes. lon_bnds = Ta->lon_bnds lat = Ta->lat ; Latitudes. lat_bnds = Ta->lat_bnds lev = Ta->lev ; Hybrid sigma pressure coordinate. lev_bnds = Ta->lev_bnds time = Ta->time(:) ; time ; ==========================================><====================================================== ; Auxiliary manipulations lixo = dimsizes(ap_bnds) hyai = new(lixo(0)+1,typeof(ap_bnds)) hyai(0:lixo(0)-1) = ap_bnds(:,0) hyai(lixo(0)) = ap_bnds(lixo(0)-1,1) hybi = hyai hybi(0:lixo(0)-1) = b_bnds(:,0) hybi(lixo(0)) = b_bnds(lixo(0)-1,1) ; ################################################################################################################### ; »»»»»»»»»»» Finally compute geopotential height «««««««««« ;;;; ap, b, hyai and hybi MUST be ordered bottom-to-top ;;;;;;;;;;;;;; zg = cz2ccm(ps,phis,Tv,p0,ap/p0,b,hyai/p0,hybi) ; Output zg is ordered top-to-bottom zg = zg(:,::-1,:,:) ; Ordering bottom-to-top ; Attach some attributes zg@standard_name = "geopotential_height" zg@long_name = "Geopotential Height" zg@units = "m" zg@comment = "Computed offline via NCL function 'cz2ccm'" ;printVarSummary(zg) p0@long_name = "vertical coordinate formula term: reference pressure" ; p0@units = "Pa" ; ; »»»»»»»»»»»»»»»»»»»»»»»»»»»»«««««««««««««««««««««««««««««« delete(ps) delete(Tv) delete(phis) ; Get dimension sizes: ntim = dimsizes(time) nlev = dimsizes(lev) nlat = dimsizes(lat) nlon = dimsizes(lon) ;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ; Writting output to netCDF file with associated meta data: ;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;------ The output file name in native grid: Zout = "zg_6hrLev_MPI-ESM1-2-LR_historical_r1i1p1f1_gn_"+sprinti("%0.4i",A)+sprinti("%0.2i",m)+".nc" ; system("/bin/rm -f " + DataDirOut+Zout) ; remove if exists system("/bin/rm -f " + Zout) ; remove if exists ; fZout = addfile (DataDirOut+Zout, "c") ; open output file fZout = addfile (Zout, "c") ; open output file ;================================================================================================== ; explicitly declare file definition mode. Improve efficiency. ;================================================================================================== setfileoption(fZout,"DefineMode",True) ;================================================================================================== ; create global attributes of the file: ;================================================================================================== ;------ Assign file attributes: fAtt = True fAtt@institution = "Max Planck Institute for Meteorology, Hamburg 20146, Germany" fAtt@institution_id = "MPI-M" fAtt@experiment_id = "historical" fAtt@activity_id = "CMIP" fAtt@source = "MPI-ESM1.2-LR (2017): aerosol: none, prescribed MACv2-SP, atmos: ECHAM6.3 (spectral T63; 192 x 96 longitude/latitude; 47 levels; top level 0.01 hPa), atmosChem: none, land: JSBACH3.20, landIce: none/prescribed, ocean: MPIOM1.63 (bipolar GR1.5, approximately 1.5deg; 256 x 220 longitude/latitude; 40 levels; top grid cell 0-12 m), ocnBgchem: HAMOCC6, seaIce: unnamed (thermodynamic (Semtner zero-layer) dynamic (Hibler 79) sea ice model)" fAtt@source_id = "MPI-ESM1-2-LR" fAtt@source_type = "AOGCM" fAtt@branch_time_in_child = "0" fAtt@branch_time_in_parent = "0" fAtt@branch_method = "standard" fAtt@contact = "cmip6-mpi-esm@dkrz.de" fAtt@history = "2019-09-04T14:33:15Z ; CMOR rewrote data to be consistent with CMIP6, CF-1.7 CMIP-6.2 and CF standards." fAtt@references = "MPI-ESM: Mauritsen, T. et al. (2019), Developments in the MPI‐M Earth System Model version 1.2 (MPI‐ESM1.2) and Its Response to Increasing CO2, J. Adv. Model. Earth Syst.,11, 998-1038, doi:10.1029/2018MS001400, Mueller, W.A. et al. (2018): A high‐resolution version of the Max Planck Institute Earth System Model MPI‐ESM1.2‐HR. J. Adv. Model. EarthSyst.,10,1383–1413, doi:10.1029/2017MS001217" fAtt@experiment = "all-forcing simulation of the recent past" fAtt@Conventions = "CF-1.7 CMIP-6.2" fAtt@project_id = "CMIP6" fAtt@table_id = "Table 6hrLev" fAtt@title = "MPI-ESM1-2-LR output prepared for CMIP6" fAtt@creation_date = "2019-09-04T14:33:15Z" fAtt@data_specs_version = "01.00.30" fAtt@external_variables = "areacella" fAtt@forcing_index = "1" fAtt@frequency = "6hrPt" fAtt@further_info_url = "https://furtherinfo.es-doc.org/CMIP6.MPI-M.MPI-ESM1-2-LR.historical.none.r1i1p1f1" fAtt@grid = "gn" fAtt@grid_label = "gn" fAtt@initialization_index = "1" fAtt@mip_era = "CMIP6" fAtt@nominal_resolution = "250 km" fAtt@physics_index = "1" fAtt@product = "model-output" fAtt@parent_experiment_id = "piControl" fAtt@parent_activity_id = "CMIP" fAtt@parent_mip_era = "CMIP6" fAtt@parent_source_id = "MPI-ESM1-2-LR" fAtt@parent_time_units = "days since 1850-1-1 00:00:00" fAtt@parent_variant_label = "r1i1p1f1" fAtt@realm = "atmos" fAtt@realization_index = "1" fAtt@sub_experiment = "none" fAtt@sub_experiment_id = "none" fAtt@table_info = "Creation Date:(09 May 2019) MD5:e6ef8ececc8f338646ebfb3aeed36bfc" fAtt@variable_id = "zg" fAtt@variant_label = "r1i1p1f1" fAtt@license = "CMIP6 model data produced by MPI-M is licensed under a Creative Commons Attribution ShareAlike 4.0 International License (https://creativecommons.org/licenses). Consult https://pcmdi.llnl.gov/CMIP6/TermsOfUse for terms of use governing CMIP6 output, including citation requirements and proper acknowledgment. Further information about this data, including some limitations, can be found via the further_info_url (recorded as a global attribute in this file) and. The data producers and data providers make no warranty, either express or implied, including, but not limited to, warranties of merchantability and fitness for a particular purpose. All liabilities arising from the supply of the information (including any liability arising in negligence) are excluded to the fullest extent permitted by law." fAtt@cmor_version = "3.5.0" fAtt@tracking_id = "hdl:21.14100/c12b3600-3326-4e58-b258-7f0de40616c9" ;------ Copy file attributes: fileattdef( fZout, fAtt ) ;================================================================================================== ; predefine the coordinate variables and their dimensionality. ; To get an UNLIMITED record dimension, the dimensionality is set to -1 and the unlimited array to True. ;================================================================================================== dimNames = (/"time", "lev", "lat", "lon", "bnds"/) dimSizes = (/ -1 , nlev , nlat, nlon , 2 /) dimUnlim = (/ True , False, False, False, False /) filedimdef(fZout,dimNames,dimSizes,dimUnlim) ;================================================================================================== ; Defining the variable's dimension name(s). ;================================================================================================== filevardef(fZout, "time" ,typeof(time) ,"time") filevardef(fZout, "lev" ,typeof(lev) ,"lev") filevardef(fZout, "lev_bnds" ,typeof(lev_bnds) ,(/"lev","bnds"/)) filevardef(fZout, "p0" ,typeof(p0) ,"ncl_scalar") filevardef(fZout, "ap" ,typeof(ap) ,"lev") filevardef(fZout, "b" ,typeof(b) ,"lev") filevardef(fZout, "ap_bnds" ,typeof(ap_bnds) ,(/"lev","bnds"/)) filevardef(fZout, "b_bnds" ,typeof(b_bnds) ,(/"lev","bnds"/)) filevardef(fZout, "lat" ,typeof(lat) ,"lat") filevardef(fZout, "lat_bnds" ,typeof(lat_bnds) ,(/"lat","bnds"/)) filevardef(fZout, "lon" ,typeof(lon) ,"lon") filevardef(fZout, "lon_bnds" ,typeof(lon_bnds) ,(/"lon","bnds"/)) filevardef(fZout, "zg" ,typeof(zg) ,(/"time","lev","lat","lon"/)) ;================================================================================================== ; Copy attributes associated with each variable to the file ; All attributes associated with each variable will be copied. ;================================================================================================== filevarattdef(fZout,"time",time) ; copy time attributes filevarattdef(fZout,"lev",lev) ; copy lev attributes filevarattdef(fZout,"lev_bnds",lev_bnds) ; copy lev_bnds attributes filevarattdef(fZout,"p0",p0) ; copy p0 attributes filevarattdef(fZout,"ap",ap) ; copy ap attributes filevarattdef(fZout,"b",b) ; copy b attributes filevarattdef(fZout,"ap_bnds",ap_bnds) ; copy ap_bnds attributes filevarattdef(fZout,"b_bnds",b_bnds) ; copy b_bnds attributes filevarattdef(fZout,"lat",lat) ; copy lat attributes filevarattdef(fZout,"lat_bnds",lat_bnds) ; copy lat_bnds attributes filevarattdef(fZout,"lon",lon) ; copy lon attributes filevarattdef(fZout,"lon_bnds",lon_bnds) ; copy lon_bnds attributes filevarattdef(fZout,"zg",zg) ; copy zg attributes ;=================================================================== ; explicitly exit file definition mode. **NOT REQUIRED** ;=================================================================== setfileoption(fZout,"DefineMode",False) ;=================================================================== ; output only the data values since the dimensionality and such have ; been predefined. The "(/", "/)" syntax tells NCL to only output the ; data values to the predefined locations on the file. ;==================================================================== fZout->time = (/time/) fZout->lev = (/lev/) fZout->lev_bnds = (/lev_bnds/) fZout->p0 = (/p0/) fZout->ap = (/ap/) fZout->b = (/b/) fZout->ap_bnds = (/ap_bnds/) fZout->b_bnds = (/b_bnds/) fZout->lat = (/lat/) fZout->lat_bnds = (/lat_bnds/) fZout->lon = (/lon/) fZout->lon_bnds = (/lon_bnds/) fZout->zg = (/zg/) ; »»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»X«««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««« ; Deleting variables with a time dimension, because it may change in the next iteration delete(time) delete(zg) end do ; End looping over months end do ; End looping over years end