MODULE MODULE_ERR_NCEP !!!--------------------------------------------------------------------------CCC !!! !!! MODULE MODULE_NCEP_ERR !!! ************************** !!! !!! PURPOSE: !!! ------- !!! OBSERVATIONAL ERROR DEFINITION !!! !! METHOD: !! ------ !! OBSERVATIONAL ERROR IS DEFINED AT 1000, 700, 500, 300, 100, 50mb !! FOR U, V, T, RH AND P !! VALUES AT ACTUAL OBS POSITION ARE OBTAINED BY INTERPOLATION BETWEEN !! THESE LEVELS (LOGARITHMIC FOR U ANF V) AND LINEAR FOR (T AND RH) !! ERRORS ARE EXTRAPOLATED BELOW 1000mb AND ABOVE 50mb !! !! INPUT: !! ----- !! OUTPUT: !! ------ !! COMMON: !! ------- !! EXTERNAL: NO !! -------- !! !! REFERENCES: !! ----------- !! PARRISH S. AND J. DERBER 1992: THE NATIONAL METEOROLOGICAL CENTERS !! SPECTRAL STATISTICAL ANALYSIS-INTERPOLATION !! SYSTEM. MONTH. WEA. REV., 120,1747-1763. !! !! !! MODIFICATIONS: !! ------------ !! ORIGINAL : 99-08 (F. VANDENBERGHE) !!----------------------------------------------------------------------------CC ! INTEGER JPERR PARAMETER (JPERR = 6) REAL ERR_K (0:JPERR+1) REAL ERR_U (0:JPERR+1), ERR_V (0:JPERR+1) REAL ERR_T (0:JPERR+1), ERR_RH (0:JPERR+1) REAL ERR_P (0:JPERR+1) REAL ERR_H (0:JPERR+1) ! PARAMETER (JPERR = 6) ! REAL ERR_K (1:JPERR) ! REAL ERR_U (1:JPERR), ERR_V (1:JPERR) ! REAL ERR_T (1:JPERR), ERR_RH (1:JPERR) ! REAL ERR_P (1:JPERR) ! !...NCEP ERRORS (U in m/s, V in m/s, T in K, H in %, P in Pa) !...PRESSURE IS OBTAINED FROM SURFACE VALUE LINEARLY DECREASING !...TO 10% OF ITS VALUE AT 10hPA !...RH HAS BEEN DIVIDED BY 2 ! DATA ERR_K (1:JPERR) & / 100000.,70000.,50000.,30000.,10000.,5000. / DATA ERR_U (1:JPERR) & / 1.4, 2.4, 2.8, 3.4, 2.5, 2.7 / DATA ERR_V (1:JPERR) & / 1.4, 2.4, 2.8, 3.4, 2.5, 2.7 / DATA ERR_T (1:JPERR) & / 1.8, 1.3, 1.3, 2.0, 3.1, 4.0 / DATA ERR_RH (1:JPERR) & / 10.0, 10.0, 10.0, 10.0, 10.0, 10.0 / DATA ERR_P (1:JPERR) & / 100.0, 72.7, 54.5, 36.4, 18.2, 13.6 / DATA ERR_H (1:JPERR) & / 7.0, 8.0, 8.6, 18.8, 39.4, 59.3 / ! REFRACTIVITY ERROR INTERPOLATION dN = NBOT*EXP[-Z/AA] ! (Adapted from Francois Vandenberghe by Y.-R. Guo 06/05/2004) ! REAL, PARAMETER :: PMAX = 1000., PMSL = 100000. ! TOP AND BOTTOM PRES PA ! REAL, PARAMETER :: NTOP = 3., NBOT = 10. ! TOP AND BOTTOM ERROR IN ! REFRACTIVITY IN N UNITS ! REAL :: AA real:: ntop , nbot real:: aa0 , aa90 , nbot0 ,nbot90 real, parameter :: pmax = 15000. , pmsl = 100000. real, parameter :: ntopp = 0.27 , nbotp0 = 3.0 , nbotp90=1.5 REAL :: AA real :: erh90 ,erh0,err0,err90,hh real,parameter :: ha = 12000. , hb = 5500. , hc = 2500., hd=0. real,parameter :: ea = 0.3 , eb = 1.3 , ec = 2.5 , ed = 1.5 CONTAINS SUBROUTINE obs_err_ncep (nobs_max, obs, number_of_obs) !------------------------------------------------------------------------------- USE module_type USE module_func USE module_intp IMPLICIT NONE INTEGER, INTENT (in) :: nobs_max TYPE (report), INTENT (inout), DIMENSION (nobs_max) :: obs TYPE (measurement ) , POINTER :: current INTEGER, INTENT (in) :: number_of_obs INTEGER :: loop_index, is_sound INTEGER :: nvalids, nsoundings INTEGER :: nsurfaces, nlevels REAL :: pres, temp, error_rh REAL :: t9, p9, rh9, qv9 REAL :: t, p, rh, qv REAL :: es9, qs9 REAL :: es, qs INCLUDE 'missing.inc' INCLUDE 'constants.inc' !------------------------------------------------------------------------------! WRITE (UNIT = 0, FMT = '(A)') & '------------------------------------------------------------------------------' WRITE ( UNIT = 0, FMT = '(A,/)') ' OBSERVATIONAL ERROR:' ! REFRACTIVITY ERROR INTERPOLATION dN = NBOT*EXP[-Z/AA] ! (Adapted from Francois Vandenberghe by Y.-R. Guo 06/05/2004) ! AA = (PMAX - PMSL) / LOG (NTOP/NBOT) ! 1. LIST RECORDS ! =============== nvalids = 0 nsoundings = 0 nsurfaces = 0 nlevels = 0 record: DO loop_index = 1, number_of_obs ! 1.1 CHECK IF RECORD VALID ! --------------------- record_valid: IF (obs(loop_index)%info%discard ) THEN CYCLE record ELSE record_valid ! 1.2 COUNT VALID RECORD ! ------------------ nvalids = nvalids + 1 ! 1.3 2D FIELDS ! --------- obs (loop_index) % ground % slp % error = 200. ! 2hPa ! Some PW or ZTD data are read with their errors, don't modify them IF (obs (loop_index) % ground % pw % error .LE. 0.) THEN IF (obs(loop_index)%info%platform(1:6) == 'FM-114') THEN obs (loop_index) % ground % pw % error = 0.005 ! .5 cm for GPSZTD ELSE obs (loop_index) % ground % pw % error = 0.2 ! 2. mm for GPSPW ENDIF if (eps_equal(obs(loop_index)%ground %pw %data,missing_r,1.)) THEN obs (loop_index) % ground % pw % qc = missing else obs (loop_index) % ground % pw % qc = 1 ! error assigned. endif ENDIF ! SSMI Tb ! obs (loop_index) % ground % tb19v % error = err_t (1) ! obs (loop_index) % ground % tb19h % error = err_t (1) ! obs (loop_index) % ground % tb22v % error = err_t (1) ! obs (loop_index) % ground % tb37v % error = err_t (1) ! obs (loop_index) % ground % tb37h % error = err_t (1) ! obs (loop_index) % ground % tb85v % error = err_t (1) ! obs (loop_index) % ground % tb85h % error = err_t (1) ! SSMI Tb varies from 2 to 5 K from channel 19 to 85 obs (loop_index) % ground % tb19v % error = 1.00 obs (loop_index) % ground % tb19h % error = 1.00 obs (loop_index) % ground % tb22v % error = 2.33 obs (loop_index) % ground % tb37v % error = 3.66 obs (loop_index) % ground % tb37h % error = 3.66 obs (loop_index) % ground % tb85v % error = 5.00 obs (loop_index) % ground % tb85h % error = 5.00 ! 1.4 INITIALISE UPPER LEVEL POINTER TO SURFACE LEVEL ! ----------------------------------------------- current => obs (loop_index) % surface ! 2. LOOP ON UPPER-AIR LEVELS (FIRST LEVEL IS SURFACE) ! ==================================================== is_sound = -1 upper_level: DO WHILE (ASSOCIATED (current)) ! 2.1 TURN ON THE SOUNDING FLAG AND COUNT THE NUMBER OF LEVEL ! ------------------------------------------------------- is_sound = is_sound + 1 nlevels = nlevels + 1 ! 2.2 CHECK IF OBS PRESSURE IS NOT MISSING ! ------------------------------------ IF ((eps_equal (current%meas%pressure%data, missing_r, 1.)) .OR. & (eps_equal (current%meas%pressure%data, 0., 1.))) THEN WRITE (0,'(A,A,1X,A)') 'Internal error obs ', & TRIM (obs (loop_index) % location % id), & TRIM (obs (loop_index) % location % name) WRITE (0,'(A,F12.3)') 'Pressure = ', current%meas%pressure%data STOP 'in obs_err_ncep.F' ENDIF ! 2.3 PRESSURE ! -------- pres = current%meas%pressure%data ! 3. VERTICAL INTERPOLATION OF NCEP OBSERVATIONAL ERROR ! ====================================================== ! 3.1 WIND DIRECTION FIXE TO 5 DEGREES ! -------------------------------- current % meas % direction % error = 5. ! 3.2 WIND SPEED ! ---------- ! Some wind speed data are read with their errors, don't modify them IF (current % meas % speed % error .LE. 0.) THEN current % meas % speed % error = intplin (pres, err_k (1:JPERR), & err_u (1:JPERR)) ENDIF ! 3.3 U-WIND COMPONENTS AS WIND SPEED ! ------------------------------- current % meas % u % error = intplin (pres, err_k (1:JPERR), & err_u (1:JPERR)) ! 3.4 V-WIND COMPONENT AS WIND SPEED ! ------------------------------ current % meas % v % error = intplin (pres, err_k (1:JPERR), & err_v (1:JPERR)) ! 3.5 PRESSURE ! -------- current % meas % pressure % error = intplin (pres, err_k (1:JPERR), & err_p (1:JPERR)) ! 3.6 HEIGHT ! ------- current % meas % height % error = intplog (pres, err_k (1:JPERR), & err_h (1:JPERR)) ! 3.7 TEMPERATURE ! ----------- current % meas % temperature % error = intplog (pres, err_k (1:JPERR), & err_t (1:JPERR)) ! 3.8 DEW POINT AS TEMPERATURE ! ------------------------ current % meas % dew_point % error = current % meas % temperature % error ! 3.9 RELATIVE HUMIDITY ! ----------------- ! current % meas % rh % error = ERROR_INTP (pres,'RH','LINEAR') current % meas % rh % error = intplog (pres, err_k (1:JPERR), & err_rh (1:JPERR)) ! 3.10 MIXING RATIO ERROR DERIVED FROM RELATIVE HUMIDITY, TEMPERATURE, PRESSURE ! ============================================================================= IF ((.NOT.eps_equal (current % meas % pressure % data,missing_r,1.)) .AND.& (.NOT.eps_equal (current % meas % temperature % data,missing_r,1.)).AND.& (.NOT.eps_equal (current % meas % rh % data, missing_r, 1.))) THEN ! Linearize the Qv calculation p9 = current % meas % pressure % data / 100 p = current % meas % pressure % error / 100 t9 = current % meas % temperature % data t = current % meas % temperature % error rh9 = current % meas % rh % data rh = current % meas % rh % error es9 = 6.112 * EXP (17.67*(t9-273.15) & /(t9-273.15+243.5)) es = 6.112 & * EXP ( 17.67*(t9-273.15) / (t9-273.15+243.5) ) & * (+17.67* t / (t9-273.15+243.5) & -17.67*(t9-273.15) / ((t9-273.15+243.5)*(t9-273.15+243.5))) qs9 = 0.622 * es9 / (p9-es9) qs = 0.622 * es / (p9-es9) & - 0.622 * es9 * (p -es ) / ((p9-es9)*(p9-es9)) qv9 = 0.01*rh9*qs9 qv = 0.01*rh *qs9 + 0.01*rh9 *qs ! Error should not be lower than 1g/kg current % meas % qv % error = MAX (qv, 0.001) ELSE current % meas % qv % error = missing_r ENDIF ! 4. GO TOP NEXT LEVEL ! ===================== current => current%next ENDDO upper_level ! 5. GO TO NEXT RECORD ! ====================== ! 5.1 INCREMENT THE SURFACE OR SOUNDING COUNTER ! ----------------------------------------- if (is_sound .gt. 0) then nsoundings = nsoundings + 1 else nsurfaces = nsurfaces + 1 endif ! 5.2 GO TO NEXT VALID RECORD ! ----------------------- ENDIF record_valid ! 5.3 GO TO NEXT RECORD ! ----------------- ENDDO record ! 6. PRINT DIAGNOSTIC ! ==================== WRITE (UNIT = 0 , FMT = '(A,I6,A,I6,A)' ) & "Number of processed stations: ",nvalids,& " = ",nlevels," levels." WRITE (UNIT = 0 , FMT = '(A,I6,A,I6,A)' ) & "Number of processed surface stations: ",nsurfaces,& " = ",nsurfaces," surface levels." WRITE (UNIT = 0 , FMT = '(A,I6,A,I6,A,/)' ) & "Number of processed upper-air stations: ",nsoundings,& " = ",nlevels-nsurfaces," upper-air levels." ! 7. END ! ======= RETURN END SUBROUTINE obs_err_ncep FUNCTION ERROR_INTP (PZ, CDVAR, CDTYPE) !!!--------------------------------------------------------------------------CCC !!! !!! FUNCTION ERROR_INTP !!! *********************** !!! !!! PURPOSE: !!! ------- !!! VERTICAL INTERPOLATION OF OBSERVATIONAL ERROR !!! !! METHOD: !! ------ !! LINEAR OR LOGARITHMIC VERTICAL INTERPOLATION !! OUT OF BOUND LOCATIONS ARE EXTRAPOLATED !! !! INPUT: !! ----- !! PZ: VERTICAL POSITION WHERE O INTERPOLATE !! CDTYPE TYPE OF VARIABLE (CURRENTLY 'U','V','T','RH','P') !! CDTYPE TYPE OF INTERPOLATION ('LOGAR' OR 'LINEAR') !! !! OUTPUT: !! ------ !! ERROR_INTP: VALUE OF ERROR AT VERTICAL LOCATION PZ !! !! COMMON: !! ------- !! SEE MODULE MODULE_NCEP_ERR !! !! EXTERNAL: NO !! -------- !! !! REFERENCES: !! ----------- !! PARRISH S. AND J. DERBER 1992: THE NATIONAL METEOROLOGICAL CENTERS !! SPECTRAL STATISTICAL ANALYSIS-INTERPOLATION !! SYSTEM. MONTH. WEA. REV., 120,1747-1763. !! MODIFICATIONS: !! ------------ !! ORIGINAL : 99-08 (F. VANDENBERGHE) !!----------------------------------------------------------------------------CC ! ! ARGUMENT ! ------- REAL PZ CHARACTER*(*) CDVAR, CDTYPE REAL ERROR_INTP INTEGER IK, JK REAL DZINF, DZSUP, DZDIF ! ! ERRORS ON MANDATORY LEVELS DEFINED IN MODULE MODULE_NCEP_ERR ! !------------------------------------------------------------------------------C ! ! ! 1. EXTRAPOLATE ERROR DEFINITION BELOW AND BEYOND LEVELS ! -------------------------------------------------------- ! ERR_K (0) = 200000. ! SET HIGHEST PRESSURE TO 20000 PA ERR_K (JPERR+1) = 1. ! SET LOWEST PRESSURE TO 1 Pa ! ERR_U (0) = ERR_U (1) ERR_U (JPERR+1) = ERR_U (JPERR) ERR_V (0) = ERR_V (1) ERR_V (JPERR+1) = ERR_V (JPERR) ERR_T (0) = ERR_T (1) ERR_T (JPERR+1) = ERR_T (JPERR) ERR_RH (0) = ERR_RH (1) ERR_RH (JPERR+1) = ERR_RH (JPERR) ERR_P (0) = ERR_P (1) ERR_P (JPERR+1) = ERR_P (JPERR) ERR_H (0) = ERR_H (1) ERR_H (JPERR+1) = ERR_H (JPERR) ! ! ! 2. FIND INDEX FOR INTERPOLATION (PRESSURE IS DECREASING WITH INDEX K) ! --------------------------------------------------------------------- ! Z ! 2000 1000 + 700 500 300 100 50 1 P in mb ! |-----------------|--------------------------------------------------> ! 0 K | K+1 3 4 5 6 7 K ! |<----->+<----->| ! DZINF DZSUP ! ! VAL = DZINF * ERR_NCEP (K+1) + DZSUP * ERR_NCEP (K) ! --------------------------------------------------------------------- ! IK = -1 ! DO JK = 0, JPERR ! IF ((ERR_K (JK) .GE. PZ) .AND. (PZ .GT. ERR_K (JK+1))) & THEN IK = JK EXIT ENDIF ! ENDDO ! IF ((IK .LT. 0) .OR. (IK .GT. JPERR+1)) THEN WRITE (0,'(/,A,A,A,F10.2,/)') ' ERROR OBS ',CDVAR, & ' POSITION PZ = ',PZ STOP ' in error_intp.F' ENDIF ! ! ! 3. WHEN SPECIFIED: LOGARITHMIC INTERPOLATION COEFICIENTS ! --------------------------------------------------------- ! IF ((CDTYPE (1:5) .EQ. 'LOGAR') .OR. & (CDTYPE (1:5) .EQ. 'logar')) THEN ! DZINF = -LOG (PZ / ERR_K (IK )) DZSUP = LOG (PZ / ERR_K (IK+1)) DZDIF = LOG (ERR_K (IK) / ERR_K (IK+1)) ! ELSE ! ! 4. OR DEFAULT IS LINEAR INTERPOLATION COEFICIENTS ! ---------------------------------------------- ! DZINF = -(PZ - ERR_K (IK )) DZSUP = PZ - ERR_K (IK+1) DZDIF = ERR_K (IK) - ERR_K (IK+1) ! ENDIF ! ! 5. INTERPOLATION ! ------------- ! DZINF = DZINF / DZDIF DZSUP = DZSUP / DZDIF ! IF (CDVAR (1:1) .EQ. 'U') THEN ! ERROR_INTP = DZSUP * ERR_U (IK) & + DZINF * ERR_U (IK+1) ! ELSE IF (CDVAR (1:1) .EQ. 'V') THEN ! ERROR_INTP = DZSUP * ERR_V (IK) & + DZINF * ERR_V (IK+1) ! ELSE IF (CDVAR (1:1) .EQ. 'T') THEN ! ERROR_INTP = DZSUP * ERR_T (IK) & + DZINF * ERR_T (IK+1) ! ELSE IF (CDVAR (1:2) .EQ. 'RH') THEN ! ERROR_INTP = DZSUP * ERR_RH (IK) & + DZINF * ERR_RH (IK+1) ! ELSE IF (CDVAR (1:1) .EQ. 'P') THEN ! ERROR_INTP = DZSUP * ERR_P (IK) & + DZINF * ERR_P (IK+1) ! ELSE IF (CDVAR (1:1) .EQ. 'H') THEN ! ERROR_INTP = DZSUP * ERR_H (IK) & + DZINF * ERR_H (IK+1) ! ELSE WRITE (0,'(/,A,A,/)') ' ERROR VARIABLE TYPE ',CDVAR STOP ' in error_intp.F' ENDIF ! ! ! 6. EXIT ! --------- ! RETURN END FUNCTION ERROR_INTP END MODULE MODULE_ERR_NCEP