;+ ; NAME: ; Create_GridFile ; ; PURPOSE: ; This procedure generates a netcdf file with grid information for ; regular grids or gaussian grids corresponding to a spectral ; truncation of T21, T31, T42, T63, T85, T106, T126, T170, T213, or ; TL319. The netcdf file contains the three variable lon, lat, and ; lat_weight. Lat_weight is the area weighting factor for latitudes ; (integral = 1). ; ; CATEGORY: ; Geography ; ; AUTHOR: ; Martin Schultz ; Max-Planck-Institute for Meteorology ; Bundesstr. 55, 20146 Hamburg, Germany ; martin.schultz@dkrz.de ; ; MODIFICATION HISTORY: ; mgs, 24 Aug 2000: VERSION 1.00 ; mgs, 11 Feb 2008: - replaced gauss weights with area ;- ; ;########################################################################### ; ; LICENSE ; ; This software is OSI Certified Open Source Software. ; OSI Certified is a certification mark of the Open Source Initiative. ; ; Copyright © 2000 Martin Schultz ; ; This software is provided "as-is", without any express or ; implied warranty. In no event will the authors be held liable ; for any damages arising from the use of this software. ; ; Permission is granted to anyone to use this software for any ; purpose, including commercial applications, and to alter it and ; redistribute it freely, subject to the following restrictions: ; ; 1. The origin of this software must not be misrepresented; you must ; not claim you wrote the original software. If you use this software ; in a product, an acknowledgment in the product documentation ; would be appreciated, but is not required. ; ; 2. Altered source versions must be plainly marked as such, and must ; not be misrepresented as being the original software. ; ; 3. This notice may not be removed or altered from any source distribution. ; ; For more information on Open Source Software, visit the Open Source ; web site: http://www.opensource.org. ; ;########################################################################### pro gauaw, nlat, lats=lats, weights=weights ;; This routine computes the gaussian latitude abscissas and ;; weights. ;; ;; Translated to IDL from ECHAM4 code written by ;; M. Jarraud, ECMWF, Dec 1982, ;; L. Kornblueh, U. Schulzweida, and A. Rhodin, MPIfM, 1999 ;; Martin Schultz, 19 JUL 2000 and bug fix 24 Aug 2000 ;; INPUT: nlat = The number of latitude coordinates that shall be ;; returned ;; OUTPUT: lats = vector of latitudes (~ grid box centers) ;; weights = vector of gaussian weights ;; Comments: ;; Gaussian latitudes are defined as latitudes where the Legendre ;; functions are zero for a given truncation. They do not ;; correspond directly to the center coordinates of a grid box. To ;; compute the grid box edges, use the gaussian weights: ;; sinlat = fltarr(nlat+1) ;; sinlat[0] = -1. ; South Pole ;; for i=1L,nlats do sinlat[i] = sinlat[i-1] + weights[i-1] ;; latedges = asin(sinlat)*180./!PI ;; The gaussian weights are also used for area integration. The ;; total of weights yields 2, thus: ;; totflux = total(weights*fluxes) * 0.5 ; !!!!!!!!!!!!!! ;; Perform calculation in double precision eps = 1.D-30 ; machine precision ITEMAX = 20 ; maximum number of iterations ;; result arrays need one extra element because FORTRAN starts at ;; index 1, and code wasn't "rewritten". lats = dblarr(nlat+1) weights = dblarr(nlat+1) ins2 = fix( nlat/2. + nlat MOD 2. ) ;; Find first approximation of the roots of the ;; Legendre polynomial of degree nlat FOR jgl = 1L, ins2 DO BEGIN z = double(4*jgl-1)*!dpi/double(4.*long(nlat)+2) lats(jgl) = COS(z+1./(TAN(z)*double(8.*long(nlat)^2))) ENDFOR ;; Computes roots and weights ;; Perform the Newton loop ;; Find 0 of Legendre polynomial with Newton loop FOR jgl = 1L, ins2 DO BEGIN za = lats[jgl] FOR iter = 1L, itemax+1 DO BEGIN zk = 0.0D ;; Newton iteration step zkm2 = 1.0D zkm1 = za zx = za FOR jn = 2L, nlat DO BEGIN zk = (double(2*jn-1)*zx*zkm1-double(jn-1)*zkm2)/double(jn) zkm2 = zkm1 zkm1 = zk ENDFOR zkm1 = zkm2 zldn = (double(nlat)*(zkm1-zx*zk))/(1.-zx*zx) zmod = -zk/zldn zxn = zx+zmod zan = zxn ;; computes weight zkm2 = 1.0 zkm1 = zxn zx = zxn FOR jn = 2L,nlat DO BEGIN zk = (double(2*jn-1)*zx*zkm1-double(jn-1)*zkm2)/double(jn) zkm2 = zkm1 zkm1 = zk ENDFOR zkm1 = zkm2 zw = (1.0-zx*zx)/(double(long(nlat)*long(nlat))*zkm1*zkm1) za = zan IF (ABS(zmod) LE eps) THEN GOTO,finished ;; Check for math error IF Check_Math(/PRINT) NE 0 THEN BEGIN help,iter,zmod,zldn stop ENDIF ENDFOR finished: lats[jgl] = zan weights[jgl] = zw * 2. ENDFOR FOR jgl = 1L, fix(nlat/2.) DO BEGIN isym = nlat-jgl+1 lats(isym) = -lats(jgl) weights(isym) = weights(jgl) ENDFOR ;; since IDL starts at index 0 but code was written in FORTRAN ;; with start index 1, omit the first value. lats = lats[1:*] weights = weights[1:*] end pro Get_Hybrid_Levels,hyai,hybi,hyam,hybm ;; read hybrid level indices and compute hyai, hybi, hyam, hybm ;; for use with Mozart ;; reference pressure p0 = 1.0e5 ;; file with ECMWF 31 level indices file = 'hybrid_indices_ecmwf.txt' readdata,file,rind,/noheader,cols=2,skp1=4 help,rind oria = reform(rind[0,*]) orib = reform(rind[1,*]) ;; ECMWF hybrid levels form pressure as p = A + B*ps ;; Mozart wants p = A*p0 + B*ps hyai = oria/p0 hybi = orib ;; !!! to prevent potential division by zero, "fix" highest level hyai[0] = 1.e-6 ;; compute mid-level indices tmp = 0.5*( hyai + shift(hyai,1) ) hyam = tmp[1:*] tmp = 0.5*( hybi + shift(hybi,1) ) hybm = tmp[1:*] end pro Get_Horizontal_Grid,resolution,lon,lat,gw ;; Get number of longitudes and latitudes for spectral truncation print,'checking for resolution ',resolution CASE resolution OF 21 : BEGIN nlon = 64 nlat = 32 END 31 : BEGIN nlon = 96 nlat = 48 END 42 : BEGIN nlon = 128 nlat = 64 END 63 : BEGIN nlon = 192 nlat = 96 END 85 : BEGIN nlon = 256 nlat = 128 END 106 : BEGIN nlon = 320 nlat = 160 END 126 : BEGIN nlon = 384 nlat = 192 END 170 : BEGIN nlon = 512 nlat = 256 END 213 : BEGIN nlon = 640 nlat = 320 END 319 : BEGIN nlon = 640 nlat = 320 END ;; continues with 511, 639, ... ELSE : BEGIN read,nlon,nlat, prompt='Please enter nlon, nlat: ' END ENDCASE ;; Compute longitudes dlon = 360./float(nlon) ;; lon = (lindgen(nlon)+0.5)*dlon ;; ### error! ECHAM grid starts with center value 0. lon = (lindgen(nlon))*dlon ;; Compute gaussian weights and latitudes ;; ATTENTION: integral(gw) = 2 ! gauaw, nlat, lats=sinlat, weights=gw ;; Convert sin(lat) to lat lat = asin(sinlat) * 180./!PI end pro Create_GridFile,resolution,regular=regular,gaussgrid=gaussgrid, $ horizontal_only=horonly ;; For now restrict to horizontal grid info horonly = 1 ;; Test grid type (default = gauss) gridtype = 'gauss' IF Keyword_Set(regular) THEN gridtype = 'regular' ;; Check resolution and compose name ;; For gauss grid, resolution must be spectral truncation number IF gridtype EQ 'gauss' THEN BEGIN IF N_Elements(resolution) NE 1 THEN BEGIN Message, 'For gaussian grid, resolution must be spectral truncation number.', $ /Continue RETURN ENDIF ENDIF ELSE BEGIN ;; For regular grid, resolution must be two-element vector with ;; delta-lon and delta-lat or four-element vector with ;; additional fields for lon0 and lat0 IF N_Elements(resolution) NE 2 AND N_Elements(resolution) NE 4 THEN BEGIN Message,'For regular grid resolution must be 2- or 4-element vector!', $ /Continue RETURN ENDIF ;; Compute lon0 and lat0 if not given IF N_Elements(resolution) EQ 2 THEN BEGIN lon0 = resolution[0]/2. lat0 = -90. + resolution[1]/2. resolution = [ resolution, lon0, lat0 ] ENDIF ENDELSE ;; ------------------------------------------- ;; get vertical grid information ;; ------------------------------------------- IF Keyword_Set(horonly) EQ 0 THEN BEGIN Get_Hybrid_Levels,hyai,hybi,hyam,hybm lev = 1000.*(hyam+hybm) ilev = 1000.*(hyai+hybi) ENDIF ;; ------------------------------------------- ;; get horizontal grid information ;; ------------------------------------------- IF gridtype EQ 'gauss' THEN BEGIN Get_Horizontal_Grid,resolution,lon,lat,latw ;; Note that latitude weights sum up to two latw = latw/2. ENDIF ELSE BEGIN ;; Get number of longitude and latitude values ;; *** assuming no offset in lon or lat nlon = long(360./resolution[0]) nlat = long(180./resolution[1]) lon = lindgen(nlon)*resolution[0] + resolution[2] lat = lindgen(nlat)*resolution[1] + resolution[3] ;; Crop arrays wx = Where(lon GE 0. AND lon LT 360., cnt) IF cnt GT 0 THEN lon = lon[wx] wy = Where(lat GT -90. AND lat LT 90., cnt) IF cnt GT 0 THEN lat = lat[wy] ;; Compute latitude interfaces lati = lat-resolution[1]/2. lati[0] = -90. lati = [ lati, 90. ] ;; Compute latitude weights latw = sin(lati * !PI/180.) latw = ( latw - Shift(latw,1) )[1:*] lattot = total(latw) latw = latw/lattot ENDELSE ;; Compute area nlon = n_elements(lon) REarth = 6371.00e3 ;; geometric mean radius of Earth in m area = ( latw ## (fltarr(nlon)+1.) ) /float(nlon) * 4.d0 * !DPI * REarth * REarth if gridtype eq 'gauss' then begin resname = 'T' + StrTrim(fix(resolution),2) endif else begin resname = 'r' + strtrim(nlon,2) + 'x' + strtrim(nlat,2) endelse ;; ------------------------------------------- ; create netCDF file ;; ------------------------------------------- ; overwite if it already exists id = NCDF_CREATE('gb'+resname+'.nc',/CLOBBER) ; define dimensions dlon = NCDF_DIMDEF(id,'lon',n_elements(lon)) dlat = NCDF_DIMDEF(id,'lat',n_elements(lat)) IF Keyword_Set(horonly) EQ 0 THEN BEGIN dlev = NCDF_DIMDEF(id,'lev',n_elements(hyam)) dilev = NCDF_DIMDEF(id,'ilev',n_elements(hyai)) ENDIF ; define variables IF Keyword_Set(horonly) EQ 0 THEN BEGIN vp0 = NCDF_VARDEF(id, 'P0', /FLOAT) NCDF_ATTPUT,id,vp0,'long_name','reference pressure',/CHAR NCDF_ATTPUT,id,vp0,'units','Pa',/CHAR ; vntrm = NCDF_VARDEF(id, 'ntrm', /LONG) ; NCDF_ATTPUT,id,vntrm,'long_name','spectral truncation parameter M',/CHAR ; vntrn = NCDF_VARDEF(id, 'ntrn', /LONG) ; NCDF_ATTPUT,id,vntrn,'long_name','spectral truncation parameter N',/CHAR ; vntrk = NCDF_VARDEF(id, 'ntrk', /LONG) ; NCDF_ATTPUT,id,vntrk,'long_name','spectral truncation parameter K',/CHAR vhyai = NCDF_VARDEF(id, 'hyai',[dilev], /FLOAT) NCDF_ATTPUT,id,vhyai,'long_name','hybrid A coefficient at layer interfaces',/CHAR vhybi = NCDF_VARDEF(id, 'hybi',[dilev], /FLOAT) NCDF_ATTPUT,id,vhybi,'long_name','hybrid B coefficient at layer interfaces',/CHAR vhyam = NCDF_VARDEF(id, 'hyam',[dlev], /FLOAT) NCDF_ATTPUT,id,vhyam,'long_name','hybrid A coefficient at layer midpoints',/CHAR vhybm = NCDF_VARDEF(id, 'hybm',[dlev], /FLOAT) NCDF_ATTPUT,id,vhybm,'long_name','hybrid B coefficient at layer midpoints',/CHAR vlev = NCDF_VARDEF(id, 'lev', [dlev], /FLOAT) NCDF_ATTPUT,id,vlev,'long_name','hybrid level at layer midpoints (1000*(A+B))',/CHAR NCDF_ATTPUT,id,vlev,'units','hybrid_sigma_pressure',/CHAR NCDF_ATTPUT,id,vlev,'positive','down',/CHAR NCDF_ATTPUT,id,vlev,'A_var','hyam',/CHAR NCDF_ATTPUT,id,vlev,'B_var','hybm',/CHAR NCDF_ATTPUT,id,vlev,'P0_var','P0',/CHAR NCDF_ATTPUT,id,vlev,'PS_var','PS',/CHAR NCDF_ATTPUT,id,vlev,'edges','ilev',/CHAR vilev = NCDF_VARDEF(id, 'ilev', [dilev], /FLOAT) NCDF_ATTPUT,id,vilev,'long_name','hybrid level at layer interfaces (1000*(A+B))',/CHAR NCDF_ATTPUT,id,vilev,'units','hybrid_sigma_pressure',/CHAR NCDF_ATTPUT,id,vilev,'positive','down',/CHAR NCDF_ATTPUT,id,vilev,'A_var','hyam',/CHAR NCDF_ATTPUT,id,vilev,'B_var','hybm',/CHAR NCDF_ATTPUT,id,vilev,'P0_var','P0',/CHAR NCDF_ATTPUT,id,vilev,'PS_var','PS',/CHAR ENDIF vlon = NCDF_VARDEF(id, 'lon', [dlon], /FLOAT) NCDF_ATTPUT,id,vlon,'long_name','longitude',/CHAR NCDF_ATTPUT,id,vlon,'units','degrees_east',/CHAR vlat = NCDF_VARDEF(id, 'lat', [dlat], /FLOAT) NCDF_ATTPUT,id,vlat,'long_name','latitude',/CHAR NCDF_ATTPUT,id,vlat,'units','degrees_north',/CHAR ; vgw = NCDF_VARDEF(id, 'lat_weight', [dlat], /FLOAT) ; NCDF_ATTPUT,id,vgw,'long_name','latitude weights',/CHAR varea = NCDF_VARDEF(id, 'grid_area', [dlon, dlat], /FLOAT) print,'varea = ',varea NCDF_ATTPUT,id,varea,'long_name','surface area of grid box',/CHAR print,'*1' NCDF_ATTPUT,id,varea,'units','m^2',/CHAR print,'*2' ; define global attributes NCDF_ATTPUT,id,/GLOBAL,'title', 'Definition for '+resname+' '+gridtype+' grid',/CHAR print,'*3' NCDF_ATTPUT,id,/GLOBAL,'author', 'Martin Schultz',/CHAR print,'*4' NCDF_ATTPUT,id,/GLOBAL,'affiliation', 'ICG-2, Research Centre Juelich',/CHAR print,'*5' NCDF_ATTPUT,id,/GLOBAL,'creation_date', strdate(),/CHAR print,'*6' ; change to DATA mode and put data into file NCDF_CONTROL,id,/ENDEF IF Keyword_Set(horonly) EQ 0 THEN BEGIN NCDF_VARPUT,id,vp0,1.0E5 ; NCDF_VARPUT,id,vntrm,long(resolution) ; NCDF_VARPUT,id,vntrn,long(resolution) ; NCDF_VARPUT,id,vntrk,long(resolution) NCDF_VARPUT,id,vhyai,hyai NCDF_VARPUT,id,vhybi,hybi NCDF_VARPUT,id,vhyam,hyam NCDF_VARPUT,id,vhybm,hybm NCDF_VARPUT,id,vlev,lev NCDF_VARPUT,id,vilev,ilev ENDIF NCDF_VARPUT,id,vlon,lon NCDF_VARPUT,id,vlat,lat ; NCDF_VARPUT,id,vgw,latw print,id,varea help,area NCDF_VARPUT,id,varea,area ; close file NCDF_CLOSE,id return end pro create_retro_templates ;; Gaussian grids: triangular spectral truncation as resolution parameter Create_GridFile,319, /gaussgrid, /horizontal_only Create_GridFile,106, /gaussgrid, /horizontal_only Create_GridFile,85, /gaussgrid, /horizontal_only Create_GridFile,63, /gaussgrid, /horizontal_only Create_GridFile,42, /gaussgrid, /horizontal_only Create_GridFile,31, /gaussgrid, /horizontal_only Create_GridFile,21, /gaussgrid, /horizontal_only ;; regular grids: delta-lon, delta-lat, lon0, lat0 as reolution parameter Create_GridFile, [5., 4., 2.5, -88.], /regular, /horizontal_only Create_GridFile, [3., 2., 1.5, -89.], /regular, /horizontal_only Create_GridFile, [2.5, 2., 1.25, -89.], /regular, /horizontal_only Create_GridFile, [2., 2., 1., -89.], /regular, /horizontal_only Create_GridFile, [1., 1., 0.5, -89.5], /regular, /horizontal_only Create_GridFile, [0.5, 0.5, 0.25, -89.75], /regular, /horizontal_only end