;+ ; ; .. _d2m_to_q2m_erai.pro: ; ; =================== ; d2m_to_q2m_erai.pro ; =================== ; ; Conversion of d2m to q2m on OAFLUX grid ; ; :file:`${PROJECT_OD}/erai_d2m_19890101_20091231_oafluxgrid.nc` have been ; produced by :ref:`interp_erai_dewt_1989_2009.pro`. ; ; :file:`${PROJECT_OD}/erai_t2m_19890101_20091231_oafluxgrid.nc` have been ; produced by :ref:`interp_erai_t2m_1989_2009.pro`. ; ; :file:`${PROJECT_OD}/erai_msl_19890101_20091231_oafluxgrid.nc` have been ; produced by :ref:`interp_erai_msl_1989_2009.pro`. ; ; q2m on OAFLUX grid is written in ; :file:`${PROJECT_OD}/erai_q2m_19890101_20091231_oafluxgrid.nc` ; if this file not already exists. ; ; This file will be used by :ref:`q2m_correction_ncdf.pro`. ; ; .. graphviz:: ; ; digraph d2m_to_q2m_erai { ; graph [ ; rankdir="LR", ; ] ; ; file_d2m [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/erai_d2m_19890101_20091231_oafluxgrid.nc"]; ; file_t2m [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/erai_t2m_19890101_20091231_oafluxgrid.nc"]; ; file_msl [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/erai_msl_19890101_20091231_oafluxgrid.nc"]; ; file_q2m [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/erai_q2m_19890101_20091231_oafluxgrid.nc"]; ; ; d2m_to_q2m_erai [shape=box, ; fontname=Courier, ; color=blue, ; URL="http://forge.ipsl.jussieu.fr/tropflux/browser/trunk/src/d2m_to_q2m_erai.pro", ; label="${PROJECT}/src/d2m_to_q2m_erai.pro"]; ; ; {file_d2m file_t2m file_msl} -> {d2m_to_q2m_erai} -> {file_q2m} ; } ; ; SEE ALSO ; ======== ; ; [Gill:AP:1982]_ ; ; [AOMIP]_ ; ; :ref:`project_profile.sh` ; ; :ref:`mooring_corrections` ; ; :ref:`interp_erai_dewt_1989_2009.pro` ; :ref:`interp_erai_msl_1989_2009.pro` ; :ref:`interp_erai_t2m_1989_2009.pro` ; ; :ref:`q2m_correction_ncdf.pro` ; ; :func:`rh_to_spechum` ; :func:`initncdf ` ; :func:`read_ncdf ` ; :func:`julday ` ; :func:`ncdf_quickwrite ` ; ; EXAMPLES ; ======== ; ; :: ; ; IDL> d2m_to_q2m_erai ; ; TODO ; ==== ; ; pb on cratos idl7 even if output file is written:: ; ; % Program caused arithmetic error: Floating overflow ; % Program caused arithmetic error: Floating illegal operand ; ; related to occurence of NaNf in ncdump output ? do no yet ; ; pk also see these error messages ... ; ; check for NaN value in P return by reading Netcdf msl file ; solution is may be to apply *mskout to q2m like in interp* tools ; not yet done to keep close as possible to pk processing ; ; coding rules ; ; EVOLUTIONS ; ========== ; ; $Id$ ; ; $URL$ ; ; - fplod 20110811T114932Z aedon.locean-ipsl.upmc.fr (Darwin) ; ; * add es0 value ; * add reference ; * replace q2m_erai by erai_q2m in output filename ; ; - fplod 20110809T154249Z aedon.locean-ipsl.upmc.fr (Darwin) ; ; * typo ; ; - fplod 20110808T120234Z cratos (Linux) ; ; * usage of ${PROJECT_OD} ; * complete description ; * remove return statement ; ; - fplod 20101215T162933Z aedon.locean-ipsl.upmc.fr (Darwin) ; ; * graph correction ; ; - fplod 20101215T093141Z aedon.locean-ipsl.upmc.fr (Darwin) ; ; * add graph in header ; ; - fplod 20101214T093615Z aedon.locean-ipsl.upmc.fr (Darwin) ; ; * minimal header ; ; - pbk 2008 ; ; * creation ; ;- pro d2m_to_q2m_erai ; @cm_4cal @cm_4data @cm_4mesh @cm_4data @cm_project ; ; test if ${PROJECT_OD} defined CASE project_od_env OF '' : BEGIN msg = 'eee : ${PROJECT_OD} is not defined' ras = report(msg) STOP END ELSE: BEGIN msg = 'iii : ${PROJECT_OD} is ' + project_od_env ras = report(msg) END ENDCASE ; ; check if output data will be possible iodirout = isadirectory(project_od_env) ; ; existence and protection for reading IF (FILE_TEST(iodirout, /DIRECTORY, /EXECUTABLE, /READ) EQ 0) THEN BEGIN msg = 'eee : the directory' + iodirout + ' is not accessible.' ras = report(msg) STOP ENDIF ; ; existence and protection for writing IF (FILE_TEST(iodirout, /DIRECTORY, /WRITE) EQ 0) THEN BEGIN msg = 'eee : the directory' + iodirout + ' was not found.' ras = report(msg) STOP ENDIF ; d1=19881001 & d2=20101231 ; ; build d2m filename filename_d2m='erai_d2m_19890101_20091231_oafluxgrid.nc' ; ; check if this file exists fullfilename_d2m = isafile(iodirout + filename_d2m, NEW=0, /MUST_EXIST) IF fullfilename_d2m[0] EQ '' THEN BEGIN msg = 'eee : the file ' + fullfilename_d2m + ' was not found.' ras = report(msg) STOP ENDIF ; ; build t2m filename filename_t2m='erai_t2m_19890101_20091231_oafluxgrid.nc' ; ; check if this file exists fullfilename_t2m = isafile(iodirout + filename_t2m, NEW=0, /MUST_EXIST) IF fullfilename_t2m[0] EQ '' THEN BEGIN msg = 'eee : the file ' + fullfilename_t2m + ' was not found.' ras = report(msg) STOP ENDIF ; ; build msl filename filename_msl='erai_msl_19890101_20091231_oafluxgrid.nc' ; ; check if this file exists fullfilename_msl = isafile(iodirout + filename_msl, NEW=0, /MUST_EXIST) IF fullfilename_msl[0] EQ '' THEN BEGIN msg = 'eee : the file ' + fullfilename_msl + ' was not found.' ras = report(msg) STOP ENDIF ; ; build output filename filename_out = 'erai_q2m_19890101_20091231_oafluxgrid.nc' fullfilename_out = iodirout + filename_out ; in order to avoid unexpected overwritten IF (FILE_TEST(fullfilename_out) EQ 1) THEN BEGIN msg = 'eee : the file ' + fullfilename_out + ' already exists.' ras = report(msg) STOP ENDIF ; ; reading erai dew point temperature initncdf, fullfilename_d2m td=read_ncdf("d2m", d1,d2, file=fullfilename_d2m,/nostr) ; reading erai air temperature initncdf, fullfilename_t2m t=read_ncdf("t2m", d1,d2, file=fullfilename_t2m,/nostr) tt=time &jpt=n_elements(time) help, t, td ;since the t2m and d2m are already in K, no need for conversion es0=6.11 e=es0*exp((25.e5/461.5)*((1/273.15)-(1/td))) es=es0*exp((25.e5/461.5)*((1/273.15)-(1/t))) rh=(e/es)*100 ;; conversion of relative humidity to sp. humidity ;; Ref.1 - Gill, Appendix 4, ;; Ref.2 - http://efdl.cims.nyu.edu/project_aomip/forcing_data/atmosphere/humidity.html t=t-273.15 initncdf, fullfilename_msl P=read_ncdf("msl", d1, d2, file=fullfilename_msl,/nostr) help, P ; ++ for debug ftp error ; ++print,P(0:1,0,0) ; ++stop q2m=rh_to_spechum(rh,t,P) ;writing field ncfile='!' + fullfilename_out time=timegen(7670, units='days', start=julday(1,1,1989)) & jpt=n_elements(time) cda0=string(jul2date(time(0)),format='(i8.8)') cda1=string(jul2date(time(jpt-1)),format='(i8.8)') tt=time-julday(1,1,1950,00,00,00) xlon=reform(glamt(*,0) ) & ylat=reform(gphit(0,*)) valmask=1.e20 lon_attr={units:'degrees_east',long_name:'Longitude'} lat_attr={units:'degrees_north',long_name:'Latitude'} q2m_attr={units:'kg/kg',missing_value:valmask,long_name:'Surface specific humidity at 2m ',short_name:'q2m',axis:'TYX'} time_attr={units:'days since 1950-01-01 00:00:00',long_name:'Time axis',time_origin:' 1950-JAN-01 00:00:00'} globattr={source:'Surface specific humidity at 2m calculated from dew point temperature',timerange:cda0+' - '+cda1} ncfields = i'q2m[longitude,latitude,time]=q2m:q2m_attr; ' $ + 'longitude[]=xlon:lon_attr; ' $ + 'latitude[]=ylat:lat_attr; ' $ + 'tt[*time]=tt:time_attr ' $ + ' @ globattr' @ncdf_quickwrite end