MODULE module_diagnostics !------------------------------------------------------------------------------ ! Compute wind components, dew point, relative humidity and mixing ratio ! when the corresponding input variables are present. ! ! D. GILL, April 1998 ! F. VANDENBERGHE, March 2001 !------------------------------------------------------------------------------ CONTAINS !------------------------------------------------------------------------------ SUBROUTINE derived_quantities (max_number_of_obs , obs , number_of_obs) !------------------------------------------------------------------------------ USE module_type USE module_func USE module_mm5 USE module_map USE module_intp IMPLICIT NONE INTEGER, INTENT (in) :: max_number_of_obs TYPE (report), DIMENSION (max_number_of_obs) :: obs INTEGER , INTENT (in) :: number_of_obs TYPE (measurement), POINTER :: current INTEGER :: i, fm REAL :: xj, yi CHARACTER (LEN= 5) :: bogus_type INCLUDE 'missing.inc' !------------------------------------------------------------------- WRITE (0,'(A)') & '------------------------------------------------------------------------------' WRITE (UNIT = 0, FMT = '(A,/)') 'DIAGNOSTICS: U, V, RH, QV etc...:' ! 1. DERIVED MODEL DEPENDENT VARIABLES ! ===================================== loop_all_obs: & DO i = 1, number_of_obs READ (obs(i) % info % platform (4:6), '(I3)') fm if (fm == 135) bogus_type = obs (i) % info % platform (8:12) valid_obs1:& IF (.NOT. obs(i)%info%discard) THEN ! 1.2 Loop over upper levels ! ---------------------- valid_obs2:& IF (ASSOCIATED (obs(i)%surface)) THEN current => obs(i)%surface all_levels:& DO WHILE (ASSOCIATED (current)) if ( fm == 281 ) then ! Set missing_r to rh (originally it is max_likehood_est for Quikscat): current%meas%rh%data = missing_r current%meas%rh%qc = missing else ! With fm = 116 (gpsref) and fm = 118 (gpseph), the field "dew_point" was ! used to store gpsref, and ssmi_retrieval fm=125, no "diagnostics_moist" is allowed. if (fm /= 116 .and. fm /= 118.and.fm /=125) then if (fm ==135.and.bogus_type /='BOGUS') then ! nothing to do if TCBOG... else CALL diagnostics_moist (current%meas) CALL diagnostics_wind (current%meas, & obs(i)%location%longitude) endif endif endif current => current%next ENDDO all_levels ENDIF valid_obs2 ENDIF valid_obs1 ENDDO loop_all_obs END SUBROUTINE derived_quantities !------------------------------------------------------------------------------ SUBROUTINE diagnostics_moist (new) !------------------------------------------------------------------------------ ! This routine computes the derived diagnostic moisture fields Td, Rh, Qv USE module_type USE module_func IMPLICIT NONE TYPE ( meas_data ) :: new REAL :: t , td, p REAL :: es, qs INCLUDE 'constants.inc' INCLUDE 'missing.inc' !------------------------------------------------------------------------------ ! Given a valid temperature and dew point, the relative humidity is ! computed. The relative humidity QC flag is arbitrarily set to the ! temperature flag value. ! 1. DISCARD NEGATIVE OR NULL TEMPERATURE AND DEW POINT ! ====================================================== IF ((eps_equal (new%temperature%data, 0., 1.)) .OR. & (new%temperature%data .LT. 0.)) THEN new%temperature%data = missing_r new%temperature%qc = zero_t_td ENDIF ! Discard negative or null dew point IF ((eps_equal (new%dew_point%data, 0., 1.)) .OR. & (new%dew_point%data .LT. 0.)) THEN new%dew_point%data = missing_r new%dew_point%qc = zero_t_td ENDIF ! 2. COMPUTE DIAGNOSTICS ! ======================== ! Compute RH if missing when both temperature and dew point are not missing IF ((.NOT. eps_equal (new%temperature%data, missing_r , 1.)) .AND. & (.NOT. eps_equal (new%dew_point%data, missing_r , 1.)) .AND. & ( eps_equal (new%rh%data, missing_r , 1.))) THEN t = new%temperature%data td = new%dew_point%data new%rh%data = 100. * exp ( L_over_Rv * (1./T - 1./Td) ) new%rh%qc = new%temperature%qc ! Compute Qv if missing when pressure is not missing ! IF ((.NOT. eps_equal (new%pressure%data, missing_r , 1.)) .AND. & IF ((new%pressure%qc >= 0) .AND. & ( eps_equal (new%qv%data, missing_r , 1.))) THEN p = new%pressure%data / 100 ! P in mb es = 6.112 * EXP (17.67*(t-273.15) & /(t-273.15+243.5)) qs = 0.622 * es /(p-es) new%qv%data = MAX (0.01*new%rh%data*qs,1.E-06); new%qv%qc = new%temperature%qc ELSE new%qv%data = missing_r new%qv%qc = zero_t_td ENDIF ELSE IF (eps_equal (new%rh%data, missing_r , 1.)) THEN new%rh%data = missing_r new%rh%qc = zero_t_td ENDIF new%qv%data = missing_r new%qv%qc = zero_t_td END IF ! 3. CHECK ON PERMISSIBLE VALUES ! ============================== ! Rh must be between and 100 IF (.NOT. eps_equal (new%rh%data, missing_r , 1.)) THEN new%rh%data = MAX (new%rh%data, 1.E-06 ) new%rh%data = MIN (new%rh%data, 100.) ENDIF ! Qv must be positive IF (.NOT. eps_equal (new%qv%data, missing_r , 1.)) THEN new%qv%data = MAX (new%qv%data, 1.E-06) ENDIF ! Td must be lower than T IF (.NOT. eps_equal (new%dew_point%data, missing_r , 1.)) THEN new%dew_point%data = MIN (new%dew_point%data, new%temperature%data) ENDIF END SUBROUTINE diagnostics_moist SUBROUTINE diagnostics_wind (new, longitude) !------------------------------------------------------------------------------ ! This routine performs some check on wind speed and direction, ! and compute the wind components on the MM5 grid. USE module_type USE module_func USE module_map IMPLICIT NONE TYPE ( meas_data ) :: new REAL :: longitude LOGICAL :: bad_wind INCLUDE 'constants.inc' INCLUDE 'missing.inc' !------------------------------------------------------------------------------ bad_wind = .false. ! The wind components are computed in two steps. First we need the ! meteorological u and v from the speed and direction. Second, those ! values of u and v are rotated to the model grid. IF ((eps_equal (new%speed%data, missing_r , 1.)) .OR. & (eps_equal (new%direction%data, missing_r , 1.))) THEN new%u%data = missing_r new%u%qc = missing new%v%data = missing_r new%v%qc = missing bad_wind = .true. ENDIF IF (eps_equal (new%speed%data, 0., 0.1)) THEN new%speed%qc = zero_spd new%u%qc = zero_spd new%v%qc = zero_spd new%u%data = missing_r new%v%data = missing_r bad_wind = .true. ENDIF IF (new%speed%data .LT. 0.) THEN new%speed%qc = negative_spd new%u%qc = negative_spd new%v%qc = negative_spd new%u%data = missing_r new%v%data = missing_r bad_wind = .true. ENDIF IF ((new%direction%data .LE. 1.e-3) .OR. & (new%direction%data .GT. 360.0)) THEN ! ! Y.-R. Guo 04/22/2004: Rizvi said WMO don't allow to report ! a zero of wind direction, here modified to be a value of 1.e-3. new%direction%qc = wrong_direction new%u%qc = wrong_direction new%v%qc = wrong_direction new%u%data = missing_r new%v%data = missing_r bad_wind = .true. ENDIF ! Y.-R. Guo 04/27/2004: to avoid assimilation of a bad wind ! with a meaningless value, set them to be missing. if (bad_wind) then new%direction%data = missing_r new%direction%qc = missing new%speed %data = missing_r new%speed %qc = missing return endif CALL FFDDUV (new%speed%data, new%direction%data,& new%u%data, new%v%data, longitude, 1) END SUBROUTINE diagnostics_wind SUBROUTINE FFDDUV (F,D,U,V,YLON,ID) !------------------------------------------------------------------------------! ! When ID = 1 ! Convert wind speed (F in m/s) and direction (D in degree 0-360) into MM5 ! wind (U-V in m/s) components ! ! When ID = -1 ! Convert MM5 wind (U-V in m/s) components into wind speed (F in m/s) and ! direction (D in degree 0-360) ! ! Need MM5 map projection parameters from module modulke_map.F ! ! IPROJ: Projection type ! PHIC: Central latitude ! XLONC: Central longitude ! XN: Cone projection ! CONV: 180/Pi ! !------------------------------------------------------------------------------! USE MODULE_MAP IMPLICIT NONE REAL, INTENT (inout) :: F,D REAL, INTENT (inout) :: U, V REAL, INTENT (in) :: YLON INTEGER, INTENT (in) :: ID REAL :: AEARTH, UEARTH, VEARTH REAL :: XLONRT, ANG INCLUDE 'constants.inc' !------------------------------------------------------------------------------! SELECT CASE (ID) CASE (1); ! CONVERT WIND MODULE/DIRECTION INTO U/V WIND COMPONENTS ON EARTH, ! THEN CONVERT U/V WIND COMPONENTS ON EARTH INTO LAMBERT CONFORMAL OR ! POLAR STEREOGRAPHIC PROJECTION U/V WIND COMPONENTS. ! ! PROJECTIONS CHANGE REQUIRES ONLY A CHANGE OF THE CONE CONSTANT, XN ! EQUATIONS REMAIN THE SAME. AEARTH = D/CONV UEARTH = -F*SIN(AEARTH) VEARTH = -F*COS(AEARTH) ! FOR CONVERSION TO GRID COORDINATES, ! SEE PROGRAM DATAMAP, SUBR VECT, AND ! ANTHES' METEO. 597 NOTES, EQUA. 2.23, 2.25, 2.28. XLONRT = XLONC-YLON IF(XLONRT .GT. 180.) XLONRT=XLONRT-360. IF(XLONRT .LT.-180.) XLONRT=XLONRT+360. ANG=XLONRT*XN/CONV ! FOR MERCATOR PROJECTION, THE WINDS ARE AS IN EARTH COORDINATES IF(IPROJ.EQ.3) ANG=0. IF(PHIC.LT.0.0) ANG=-ANG U = VEARTH*SIN(ANG) + UEARTH*COS(ANG) V = VEARTH*COS(ANG) - UEARTH*SIN(ANG) ! CONVERT LAMBERT CONFORMAL OR POLAR STEREOGRAPHIC PROJECTION U/V ! WIND COMPONENTS INTO U/V WIND COMPONENTS ON EART ! THEN CONVERT U/V WIND COMPONENTS ON EARTH INTO WIND MODULE/DIRECTION ! ! PROJECTIONS CHANGE REQUIRES ONLY A CHANGE OF THE CONE CONSTANT, XN CASE (-1); XLONRT = XLONC-YLON IF(XLONRT .GT. 180.) XLONRT=XLONRT-360. IF(XLONRT .LT.-180.) XLONRT=XLONRT+360. ANG=XLONRT*XN/CONV ! ! FOR MERCATOR PROJECTION, THE WINDS ARE AS IN EARTH COORDINATES ! IF(IPROJ .EQ. 3) ANG = 0. IF(PHIC .LT. 0.) ANG = -ANG UEARTH = U*COS(ANG) - V*SIN(ANG) VEARTH = U*SIN(ANG) + V*COS(ANG) F = SQRT(UEARTH*UEARTH + VEARTH*VEARTH) IF (F .EQ. 0.0) THEN D = 0. RETURN ENDIF IF (VEARTH .EQ. 0.) THEN IF (UEARTH .GT. 0.) D = 270. IF (UEARTH .LT. 0.) D = 90. ELSE AEARTH = ATAN (UEARTH/VEARTH)*CONV IF (UEARTH .LE. 0.0 .AND. VEARTH .LE. 0.0 ) D = AEARTH IF (UEARTH .LE. 0.0 .AND. VEARTH .GE. 0.0 ) D = AEARTH + 180.0 IF (UEARTH .GE. 0.0 .AND. VEARTH .GE. 0.0 ) D = AEARTH + 180.0 IF (UEARTH .GE. 0.0 .AND. VEARTH .LE. 0.0 ) D = AEARTH + 360.0 ENDIF CASE DEFAULT WRITE (0,'(/,A,I2,/)') ' UNKNOWN OPTION ',ID STOP ' in ffdduv.F' END SELECT END SUBROUTINE FFDDUV !------------------------------------------------------------------------------! END MODULE module_diagnostics