;+ ; ; =================== ; d2m_to_q2m_erai.pro ; =================== ; ; .. function:: d2m_to_q2m_erai(yyyymmddb,yyyymmdde) ; ; DESCRIPTION ; =========== ; ; Conversion of d2m to q2m on OAFLUX grid ; ; :file:`${PROJECT_OD}/erai_d2m_{yyyymmdd}_{yyyymmdd}_oafluxgrid.nc` ; containing ; ++ ; has been produced by ; :func:`interp_erai_dewt`. ; ; :file:`${PROJECT_OD}/erai_t2m_{yyyymmdd}_{yyyymmdd}_oafluxgrid.nc` ; containing ; ++ ; has been produced by ; :func:`interp_erai_t2m`. ; ; :file:`${PROJECT_OD}/erai_msl_{yyyymmdd}_{yyyymmdd}_oafluxgrid.nc` ; containing ; ++ ; has been ; produced by ; :func:`interp_erai_msl`. ; ; q2m on OAFLUX grid ; is written in ; :file:`${PROJECT_OD}/erai_q2m_{yyyymmdd}_{yyyymmdd}_oafluxgrid.nc` ; if this file not already exists. ; ; This file ; :file:`${PROJECT_OD}/erai_q2m_{yyyymmdd}_{yyyymmdd}_oafluxgrid.nc` ; will be used by ; :func:`q2m_correction_ncdf`. ; ; .. only:: man ; ; Figure is visible on PDF and HTML documents only. ; ; .. only:: html or latex ; ; .. graphviz:: ; ; digraph d2m_to_q2m_erai { ; ; file_d2m [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/erai_d2m_{yyyymmdd}_{yyyymmdd}_oafluxgrid.nc"]; ; file_t2m [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/erai_t2m_{yyyymmdd}_{yyyymmdd}_oafluxgrid.nc"]; ; file_msl [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/erai_msl_{yyyymmdd}_{yyyymmdd}_oafluxgrid.nc"]; ; file_q2m [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/erai_q2m_{yyyymmdd}_{yyyymmdd}_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 ; ======== ; ; :ref:`mooring_corrections` ; ; [Gill:AP:1982]_ ; ; [AOMIP]_ ; ; Used by :ref:`tropflux.sh` ; ; :ref:`project_profile.sh` ; ; :func:`interp_erai_dewt` ; :func:`interp_erai_msl` ; :func:`interp_erai_t2m` ; ; :func:`q2m_correction_ncdf` ; ; :func:`rh_to_spechum` ; :func:`initncdf ` ; :func:`read_ncdf ` ; :func:`julday ` ; :func:`ncdf_quickwrite ` ; ; EXAMPLES ; ======== ; ; .. code-block:: idl ; ; yyyymmddb = 20000101L ; yyyymmdde = 20001231L ; result = d2m_to_q2m_erai(yyyymmddb, yyyymmdde) ; print, result ; ; TODO ; ==== ; ; fix time pb: ; ; .. parsed-literal:: ; ; % the time axis has no date before date 2: 20000302 ; d2m first date 29811221 ; d2m last date 29851224 ; ; compare time between d2m and t2m and msl after reading ; ; pb on cratos idl7 even if output file is written: ; ; .. parsed-literal:: ; ; % Program caused arithmetic error: Floating overflow ; % Program caused arithmetic error: Floating illegal operand ; ; related to occurrence 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 20120320 ; ; * try to add compile_opt seems to be incompatible with ncdf_quickwrite ; * pro -> function with yyyymmddb and yyyymmdde parameters ; * hard coded d1 and d2 replaced by yyyymmddb and yyyymmdde parameters ; * start to homogenize time handling regarding :func:`interp_erai_t2m` ; ; - 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 ; ;- function d2m_to_q2m_erai $ , yyyymmddb $ , yyyymmdde ; compile_opt idl2, strictarrsubs, logical_predicate ; @cm_4cal @cm_4data @cm_4mesh @cm_4data @cm_project ; ; Return to caller if errors ;+++ ON_ERROR, 2 ; result = -1 ; usage = 'result = d2m_to_q2m_erai(yyyymmddb, yyyymmdde)' nparam = N_PARAMS() IF (nparam NE 2) THEN BEGIN ras = report(['Incorrect number of arguments.' $ + '!C' $ + 'Usage : ' + usage]) return, result ENDIF ; ; test if ${PROJECT_OD} defined CASE project_od_env OF '' : BEGIN msg = 'eee : ${PROJECT_OD} is not defined' ras = report(msg) return, result 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) return, result 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) return, result ENDIF ; ; build d2m filename filename_d2m='erai_d2m_' + string(yyyymmddb,format='(I8.8)') + '_' + string(yyyymmdde,format='(I8.8)') + '_oafluxgrid.nc' ; ; check if this file exists msg='iii : looking for ' + filename_d2m ras = report(msg) 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) return, result ENDIF ; ; build t2m filename filename_t2m='erai_t2m_' + string(yyyymmddb,format='(I8.8)') + '_' + string(yyyymmdde,format='(I8.8)') + '_oafluxgrid.nc' ; ; check if this file exists msg='iii : looking for ' + filename_t2m ras = report(msg) 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) return, result ENDIF ; ; build msl filename filename_msl='erai_msl_' + string(yyyymmddb,format='(I8.8)') + '_' + string(yyyymmdde,format='(I8.8)') + '_oafluxgrid.nc' ; ; check if this file exists msg='iii : looking for ' + filename_msl ras = report(msg) 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) return, result ENDIF ; ; build output filename filename_out = 'erai_q2m_' + string(yyyymmddb,format='(I8.8)') + '_' + string(yyyymmdde,format='(I8.8)') + '_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) return, result ENDIF ; ; reading erai dew point temperature initncdf, fullfilename_d2m td=read_ncdf("d2m", yyyymmddb, yyyymmdde, file=fullfilename_d2m,/nostr) timein=24.d*(time-julday(1,1,1957,0,0,0)) jpt=n_elements(timein) da=jul2date(time[0]) cda0=string(da,format='(i8.8)') da=jul2date(time[jpt-1]) cda1=string(da,format='(i8.8)') print, 'd2m first date ', cda0 print, 'd2m last date ' , cda1 ; ; reading erai air temperature initncdf, fullfilename_t2m t=read_ncdf("t2m", yyyymmddb, yyymmdde, file=fullfilename_t2m,/nostr) timein=24.d*(time-julday(1,1,1957,0,0,0)) jpt=n_elements(timein) da=jul2date(time[0]) cda0=string(da,format='(i8.8)') da=jul2date(time[jpt-1]) cda1=string(da,format='(i8.8)') print, 't2m in d2m_to_q2m_erai first date ', cda0 print, 't2m in d2m_to_q2m_erai last date ' , cda1 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", yyyymmddb, yyyymmdde, file=fullfilename_msl,/nostr) timein=24.d*(time-julday(1,1,1957,0,0,0)) jpt=n_elements(timein) da=jul2date(time[0]) cda0=string(da,format='(i8.8)') da=jul2date(time[jpt-1]) cda1=string(da,format='(i8.8)') print, 'msl in d2m_to_q2m_erai first date ', cda0 print, 'msl in d2m_to_q2m_erai last date ' , cda1 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 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 = 'q2m[longitude,latitude,*time]=q2m:q2m_attr; ' $ + 'longitude[]=xlon:lon_attr; ' $ + 'latitude[]=ylat:lat_attr; ' $ + 'time[]=timein:time_attr ' $ + ' @ globattr' ; @ncdf_quickwrite ; result = 0 return, result ; end