MODULE NOAA_LoadAuxData ! Routine to read in the two formats of model data used at NOAA ! GRID and HDF USE GbcsConstants USE GbcsTypes USE GbcsDateTime USE MCIDAS_LoadImagery USE GbcsAtmPhysics USE GbcsPDFLoaders USE GbcsCovMatLoaders USE GbcsErrorHandler USE GbcsForecastModel USE GbcsSystemTools, ONLY: Get_New_File_Unit USE MCIDAS_LoadImagery USE NOAA_OzoneClimatology USE NOAA_GetAncilliaryData ! Netcdf #ifdef USE_NETCDF USE NetCDF #endif #ifdef USE_GRIB USE NOAA_ReadGrib #endif ! Units conversion USE Type_Kinds USE Units_Conversion, rh_to_mr_units=>rh_to_mr,& mr_to_pp_units=>mr_to_pp,& pp_to_md_units=>pp_to_md IMPLICIT NONE integer, save:: allocate_grid_header = 0 integer, parameter :: ngrid_header = 64 integer :: grid_header(ngrid_header) ! GRID file header character(len=4) :: cgrid_header(ngrid_header) CHARACTER(LEN=256) :: errorString PRIVATE PUBLIC :: Load_Aux_Data PUBLIC :: NOAA_Load_Cov_Mat ! PUBLIC :: TOA_PRESSURE ! PUBLIC :: TOA_TEMPERATURE ! PUBLIC :: TOA_RH #ifdef USE_HDF PUBLIC :: checkHDF #endif CONTAINS ! Load auxiliary data SUBROUTINE Load_Aux_Data( SAT , IMG , AUX , RTM, renorm_highres_sst ) #ifndef OLD_PDF_FILES USE GbcsPDFLookup #endif #ifdef USE_NETCDF USE BOM_LoadAuxData #endif TYPE(Satellite), INTENT(INOUT) :: SAT TYPE(Imagery), INTENT(INOUT) :: IMG TYPE(Aux_Data), INTENT(INOUT) :: AUX TYPE(RTM_Interface), INTENT(INOUT) :: RTM LOGICAL, OPTIONAL, INTENT(IN) :: renorm_highres_sst INTEGER :: nelems INTEGER :: nlines INTEGER :: nlevels INTEGER :: STAT INTEGER :: I_Chan INTEGER :: I_Map INTEGER :: RTM_Num INTEGER :: RTM_ID INTEGER :: ii INTEGER :: i,j LOGICAL :: Get_Aerosol LOGICAL :: renorm_hres_sst renorm_hres_sst = .FALSE. IF( PRESENT(renorm_highres_sst) )THEN renorm_hres_sst = renorm_highres_sst ENDIF if( AUX%Model%DataFrmt .eq. "MCIDAS" )then ! MCIDAS GRID files call MCIDAS_Load_Aux_Data( SAT, IMG, AUX, RTM ) else if( AUX%Model%DataFrmt .eq. "MCIDAS26" )then ! MCIDAS GRID files call MCIDAS_Load_Aux_Data( SAT, IMG, AUX, RTM, Level26=.TRUE. ) #ifdef USE_NETCDF else if( AUX%Model%DataFrmt .eq. 'BOM-NETCDF-GASP' )then call BOM_NetCDF_Load_Aux_Data( SAT, IMG, AUX, RTM ) #endif #ifdef USE_HDF else if( AUX%Model%DataFrmt .eq. "NOAA-HDF-ARCHIVE-OPE" )then ! Read archived HDF files but using same levels as operational code CALL NOAA_HDF_Load_Aux_Data( SAT, IMG, AUX, RTM, .TRUE. ) else if( AUX%Model%DataFrmt .eq. "NOAA-HDF-ARCHIVE" )then ! Read archived HDF files CALL NOAA_HDF_Load_Aux_Data( SAT, IMG, AUX, RTM, .FALSE. ) #endif #ifdef USE_GRIB else if( AUX%Model%DataFrmt .eq. "NOAA-GRIB-ECMWF" )then ! Read archived ECMWF 60 level data CALL NOAA_Grib_ECMWF_Load_Aux_Data( SAT, IMG, AUX, RTM ) #endif else WRITE(errorString,'('' ERROR: Cannot find '',a,'' as a valid format'')')& TRIM(AUX%Model%DataFrmt) CALL Gbcs_Critical( .TRUE., & errorString, & 'Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) endif nlevels = AUX%Model%Grid%NLevels nelems = AUX%Model%Grid%NElems nlines = AUX%Model%Grid%NLines ! Load Aerosol if required ! Note that this required NetCDF #ifdef USE_NETCDF ! Check rtm_id Get_Aerosol = .FALSE. DO RTM_Num = 1,RTM%N_RTMs ! Process location using all RTMs RTM_id = RTM%ID_Map( RTM_Num ) IF( Gbcs_CRTM_VIS_ID .eq. RTM_Id )THEN Get_Aerosol = .TRUE. ENDIF IF( 1 .eq. RTM%RTM_Set(RTM_id)%Use_Solar_Angle )THEN Get_Aerosol = .TRUE. ENDIF END DO ! If needed, load Aerosol into the model grid IF( Get_Aerosol )THEN IF ( Gbcs_Verbosity > 2 ) & WRITE(Gbcs_Log_File_Unit,*)& ' Load Forecast Model Aerosol from climatology' CALL Get_Aerosol_from_Clim( AUX%GbcsDataPath, IMG%Mean_Time, AUX%Model ) ELSE AUX%Model%NAerosol = 0 ENDIF #else AUX%Model%NAerosol = 0 #endif ! Covariance Matrix (B) IF ( Gbcs_Verbosity > 2 ) & WRITE(Gbcs_Log_File_Unit,*)' Load Forecast Model Covariance Matrix' STAT = NOAA_Load_Cov_Mat( AUX%GbcsDataPath , AUX%CovMatFile , AUX%CovMat ) IF ( Gbcs_Verbosity > 3 ) THEN WRITE(Gbcs_Log_File_Unit,*)' Matrix Size = ',AUX%CovMat%MSize,' x ',AUX%CovMat%MSize WRITE(Gbcs_Log_File_Unit,*)' Number of Bands = ',AUX%CovMat%NBands DO ii = 1,AUX%CovMat%NBands WRITE(Gbcs_Log_File_Unit,FMT='(6X,I3,2F7.2)')ii,AUX%CovMat%Limits(ii,:) END DO END IF ! RTM Tangent Linear Field ! PDFs #ifdef OLD_PDF_FILES ! Read in pre netCDF files IF ( Gbcs_Verbosity > 2 ) & WRITE(Gbcs_Log_File_Unit,*)& ' Load Cloudy Sky PDF (Old) Look Up Tables' IF ( Gbcs_Verbosity > 3 ) THEN WRITE(Gbcs_Log_File_Unit,*)' Spectral PDF Numbers = ',AUX%SN_PDF_Num_Str,' ',AUX%ST_PDF_Num_Str,' ',AUX%SD_PDF_Num_Str WRITE(Gbcs_Log_File_Unit,*)' Textural PDF Numbers = ',AUX%TN_PDF_Num_Str,' ',AUX%TT_PDF_Num_Str,' ',AUX%TD_PDF_Num_Str END IF STAT = Load_PDF_LUT( SAT , AUX , 0 , Gbcs_Nighttime , AUX%PDF_SN ) ! Spectral Nighttime IF ( STAT > 0 ) THEN CALL Gbcs_Critical( .TRUE., & 'Error loading PDF LUT', & 'Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF STAT = Load_PDF_LUT( SAT , AUX , 0 , Gbcs_Twilight , AUX%PDF_ST ) ! Spectral Twilight IF ( STAT > 0 ) THEN CALL Gbcs_Critical( .TRUE., & 'Error loading PDF LUT', & 'Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF STAT = Load_PDF_LUT( SAT , AUX , 0 , Gbcs_Daytime , AUX%PDF_SD ) ! Spectral Daytime IF ( STAT > 0 ) THEN CALL Gbcs_Critical( .TRUE., & 'Error loading PDF LUT', & 'Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF STAT = Load_PDF_LUT( SAT , AUX , 1 , Gbcs_Nighttime , AUX%PDF_TN ) ! Textural Nighttime IF ( STAT > 0 ) THEN CALL Gbcs_Critical( .TRUE., & 'Error loading PDF LUT', & 'Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF STAT = Load_PDF_LUT( SAT , AUX , 1 , Gbcs_Twilight , AUX%PDF_TT ) ! Textural Twilight IF ( STAT > 0 ) THEN CALL Gbcs_Critical( .TRUE., & 'Error loading PDF LUT', & 'Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF STAT = Load_PDF_LUT( SAT , AUX , 1 , Gbcs_Daytime , AUX%PDF_TD ) ! Textural Daytime IF ( STAT > 0 ) THEN CALL Gbcs_Critical( .TRUE., & 'Error loading PDF LUT', & 'Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF STAT = Load_CLR_PDF_LUT( SAT , AUX , 1 , Gbcs_Nighttime , AUX%PDF_TN_CLR ) ! Textural Clear Scene Nighttime IF ( STAT > 0 ) THEN CALL Gbcs_Warning( .TRUE., & 'Problem loading PDF LUT (Clear Scene Nighttime) - using theoretical dist', & 'Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ! Set to same set of filters as Cloudy textural AUX%PDF_TN_CLR%Spectral = .FALSE. AUX%PDF_TN_CLR%NChans = AUX%PDF_TN%NChans AUX%PDF_TN_CLR%Chan_ID = AUX%PDF_TN%Chan_ID ENDIF STAT = Load_CLR_PDF_LUT( SAT , AUX , 1 , Gbcs_Twilight , AUX%PDF_TT_CLR ) ! Textural Clear Scene Twilight IF ( STAT > 0 ) THEN CALL Gbcs_Warning( .TRUE., & 'Problem loading PDF LUT (Clear Scene Twilight) - using theoretical dist', & 'Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ! Set to same set of filters as Cloudy textural AUX%PDF_TT_CLR%Spectral = .FALSE. AUX%PDF_TT_CLR%NChans = AUX%PDF_TT%NChans AUX%PDF_TT_CLR%Chan_ID = AUX%PDF_TT%Chan_ID ENDIF STAT = Load_CLR_PDF_LUT( SAT , AUX , 1 , Gbcs_Daytime , AUX%PDF_TD_CLR ) ! Textural Clear Scene Daytime IF ( STAT > 0 ) THEN CALL Gbcs_Warning( .TRUE., & 'Problem loading PDF LUT (Clear Scene Daytime) - using theoretical dist', & 'Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ! Set to same set of filters as Cloudy textural AUX%PDF_TD_CLR%Spectral = .FALSE. AUX%PDF_TD_CLR%NChans = AUX%PDF_TD%NChans AUX%PDF_TD_CLR%Chan_ID = AUX%PDF_TD%Chan_ID ENDIF #else ! New NetCDF PDF files IF ( Gbcs_Verbosity > 2 ) & WRITE(Gbcs_Log_File_Unit,*)' Load Cloudy Sky PDF Look Up Tables' IF ( Gbcs_Verbosity > 3 ) THEN WRITE(Gbcs_Log_File_Unit,*)' Spectral PDF Numbers = ', & AUX%SN_PDF_Num_Str,' ',AUX%ST_PDF_Num_Str,' ',AUX%SD_PDF_Num_Str WRITE(Gbcs_Log_File_Unit,*)' Textural PDF Numbers = ', & AUX%TN_PDF_Num_Str,' ',AUX%TT_PDF_Num_Str,' ',AUX%TD_PDF_Num_Str END IF IF( Gbcs_GOES .eq. SAT%Platform_ID )THEN SELECT CASE( SAT%Platform_Num ) CASE(11) AUX%pdf_filename = TRIM(AUX%GbcsDataPath) // & 'GOES11_cloud_pdf.nc' CASE(12) AUX%pdf_filename = TRIM(AUX%GbcsDataPath) // & 'GOES12_cloud_pdf.nc' CASE(13) AUX%pdf_filename = TRIM(AUX%GbcsDataPath) // & 'GOES13_cloud_pdf.nc' CASE(14) AUX%pdf_filename = TRIM(AUX%GbcsDataPath) // & 'GOES14_cloud_pdf.nc' END SELECT ELSE IF( Gbcs_MTSAT .eq. SAT%Platform_ID )THEN SELECT CASE( SAT%Platform_Num ) CASE(1) AUX%pdf_filename = TRIM(AUX%GbcsDataPath) // & 'GOES11_cloud_pdf.nc' END SELECT ELSE IF( Gbcs_MSG .eq. SAT%Platform_ID )THEN SELECT CASE( SAT%Platform_Num ) CASE(1) AUX%pdf_filename = TRIM(AUX%GbcsDataPath) // & 'MSG1_cloud_pdf.nc' CASE(2) AUX%pdf_filename = TRIM(AUX%GbcsDataPath) // & 'MSG2_cloud_pdf.nc' CASE(3) AUX%pdf_filename = TRIM(AUX%GbcsDataPath) // & 'MSG3_cloud_pdf.nc' END SELECT ENDIF ! -- List of channels for D/T/N using Gbcs indicies ! WRITE(*,*) '---' ! WRITE(*,*) RTM%Gbcs_Map(AUX%Bayes_Spec_Chans(1:AUX%N_Bayes_Chans(1,1),1)) ! WRITE(*,*) RTM%Gbcs_Map(AUX%Bayes_Spec_Chans(1:AUX%N_Bayes_Chans(1,2),2)) ! WRITE(*,*) RTM%Gbcs_Map(AUX%Bayes_Spec_Chans(1:AUX%N_Bayes_Chans(1,3),3)) ! WRITE(*,*) '---' STAT = Load_PDF(AUX%pdf_filename, AUX%pdf_TIR_min, 'tir_110_120') IF ( STAT > 0 ) STOP STAT = Load_PDF(AUX%pdf_filename, AUX%pdf_TIR_max, 'tir_037_110_120') IF ( STAT > 0 ) STOP STAT = Load_PDF(AUX%pdf_filename, AUX%pdf_VIS, AUX%VIS_PDF_ID ) IF ( STAT > 0 ) STOP STAT = Load_PDF(AUX%pdf_filename, AUX%pdf_TXT, 'txt_110') IF ( STAT > 0 ) STOP ! ! -- Load clear-sky texture PDF ! AUX%pdf_filename = TRIM(AUX%GbcsDataPath) // 'AATSR_clear.nc' ! STAT = Load_PDF(AUX%pdf_filename, AUX%pdf_TXT_clr, 'txt_110') ! IF ( STAT > 0 ) STOP STAT = Load_CLR_PDF_LUT( SAT , AUX , 1 , Gbcs_Nighttime , AUX%PDF_TN_CLR ) ! Textural Clear Scene Nighttime IF ( STAT > 0 ) STOP STAT = Load_CLR_PDF_LUT( SAT , AUX , 1 , Gbcs_Twilight , AUX%PDF_TT_CLR ) ! Textural Clear Scene Twilight IF ( STAT > 0 ) STOP STAT = Load_CLR_PDF_LUT( SAT , AUX , 1 , Gbcs_Daytime , AUX%PDF_TD_CLR ) ! Textural Clear Scene Daytime IF ( STAT > 0 ) STOP #endif AUX%PDF_Loaded = .TRUE. ! Now reset ocean pixels in model using climatology CALL Reset_Model_Surface_Temp( AUX, IMG ) ! Load in SST climatology IF( AUX%Use_HighRes_SST )THEN CALL Load_SST_Climatology( IMG, AUX, AUX%SST_Clim, IMG%Start_Time ) ! Must always do for Casey & Cornillon IF( 'NODC Pathfinder climatology (Casey & Cornillon 1999)'.eq.& AUX%SST_Clim%Name )THEN ! Re-normalize climatology to MODEL surface temperature CALL Renormalize_HighRes_SST( AUX, AUX%SST_Clim%SST(:,:,1), & AUX%SST_Clim ) ELSE IF( renorm_hres_sst )THEN ! If wanbt to renorm other inputs ! Re-normalize climatology to MODEL surface temperature CALL Renormalize_HighRes_SST( AUX, AUX%SST_Clim%SST(:,:,1), & AUX%SST_Clim ) ENDIF ENDIF AUX%Aux_Data_Loaded = .TRUE. END SUBROUTINE Load_Aux_Data #ifdef USE_NETCDF ! Read in Aerosol data from a climatology ! Based on Prabhat Koner dataset and matlab code SUBROUTINE Get_Aerosol_from_Clim( GbcsDataDir, Image_Time, Model ) USE NETCDF USE GbcsProfileGenerator USE GbcsInterpolators USE GbcsDateTime CHARACTER(LEN=*), INTENT(IN) :: GbcsDataDir TYPE(DateTime), INTENT(IN) :: Image_Time TYPE(Forecast_Model), INTENT(INOUT) :: Model ! Local variables INTEGER :: I,J,K INTEGER :: stat INTEGER :: ncid INTEGER :: nlon, nlat INTEGER :: nlev INTEGER :: ntime INTEGER :: time_pos INTEGER :: LonDimID INTEGER :: LatDimID INTEGER :: LevDimID INTEGER :: TimeDimID INTEGER :: varID INTEGER :: Aerosol_Indx(8) REAL(GbcsDble) :: pressure_val REAL(GbcsDble) :: airdn REAL(GbcsDble) :: dust_r REAL(GbcsDble) :: dust_r_weight INTEGER, ALLOCATABLE :: int_lev(:) REAL, ALLOCATABLE :: int_fact(:) INTEGER :: int_srf REAL, ALLOCATABLE :: press(:) REAL, ALLOCATABLE :: prof(:) REAL :: prof_2m REAL(GbcsDble) :: P0 ! REAL(GbcsDble), ALLOCATABLE :: P0(:) REAL(GbcsDble), ALLOCATABLE :: lon(:) REAL(GbcsDble), ALLOCATABLE :: lat(:) REAL(GbcsDble), ALLOCATABLE :: time(:) REAL(GbcsDble), ALLOCATABLE :: hyam(:) REAL(GbcsDble), ALLOCATABLE :: hybm(:) REAL(GbcsDble), ALLOCATABLE :: surface_pressure(:,:,:) REAL(GbcsDble), ALLOCATABLE :: temp(:,:,:,:) REAL(GbcsDble), ALLOCATABLE :: cb1(:,:,:,:) REAL(GbcsDble), ALLOCATABLE :: cb2(:,:,:,:) REAL(GbcsDble), ALLOCATABLE :: dst01(:,:,:,:) REAL(GbcsDble), ALLOCATABLE :: dst02(:,:,:,:) REAL(GbcsDble), ALLOCATABLE :: dst03(:,:,:,:) REAL(GbcsDble), ALLOCATABLE :: dst04(:,:,:,:) REAL(GbcsDble), ALLOCATABLE :: oc1(:,:,:,:) REAL(GbcsDble), ALLOCATABLE :: oc2(:,:,:,:) REAL(GbcsDble), ALLOCATABLE :: so4(:,:,:,:) REAL(GbcsDble), ALLOCATABLE :: sslt01(:,:,:,:) REAL(GbcsDble), ALLOCATABLE :: sslt02(:,:,:,:) REAL(GbcsDble), ALLOCATABLE :: sslt03(:,:,:,:) REAL(GbcsDble), ALLOCATABLE :: sslt04(:,:,:,:) REAL(GbcsDble), ALLOCATABLE :: temp_cb(:,:,:) REAL(GbcsDble), ALLOCATABLE :: temp_dst(:,:,:) REAL(GbcsDble), ALLOCATABLE :: temp_oc(:,:,:) REAL(GbcsDble), ALLOCATABLE :: temp_so4(:,:,:) REAL(GbcsDble), ALLOCATABLE :: temp_sslt01(:,:,:) REAL(GbcsDble), ALLOCATABLE :: temp_sslt02(:,:,:) REAL(GbcsDble), ALLOCATABLE :: temp_sslt03(:,:,:) REAL(GbcsDble), ALLOCATABLE :: temp_sslt04(:,:,:) REAL(GbcsDble), ALLOCATABLE :: temp_sfc_pressure(:,:) REAL(GbcsDble), ALLOCATABLE :: temp_dustr(:,:) REAL, ALLOCATABLE :: temp_arr1(:,:) REAL, ALLOCATABLE :: temp_arr2(:,:) REAL, ALLOCATABLE :: temp_arr3(:,:) REAL, ALLOCATABLE :: temp_arr4(:,:) TYPE(DateTime) :: Start_Time TYPE(DateTime) :: Time_Val CHARACTER(LEN=512) :: Filename INTEGER, ALLOCATABLE :: surface_type(:,:) TYPE(Grid_Descriptor) :: Grid TYPE(Atm_Profile) :: Profile TYPE(ImagePixel) :: PX ! Read in Aerosol data Filename = TRIM(GbcsDataDir)//'/aero_1.9x2.5_L26_1990-1999.nc' stat = nf90_open(Filename,NF90_NOWRITE,ncid) IF( NF90_NOERR .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error opening NETCDF aerosol file',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Get dimension values stat = nf90_inq_dimid(ncid, "lon", LonDimID) IF(stat .ne. nf90_noerr)THEN CALL Gbcs_Critical( .TRUE., & 'Error getting longitude dim id',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inquire_dimension(ncid, LonDimID, len = nlon) IF(stat .ne. nf90_noerr)THEN CALL Gbcs_Critical( .TRUE., & 'Error getting longitude length',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inq_dimid(ncid, "lat", LatDimID) IF(stat .ne. nf90_noerr)THEN CALL Gbcs_Critical( .TRUE., & 'Error getting latitude dim id',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inquire_dimension(ncid, LatDimID, len = nlat) IF(stat .ne. nf90_noerr)THEN CALL Gbcs_Critical( .TRUE., & 'Error getting latitude length',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inq_dimid(ncid, "lev", LevDimID) IF(stat .ne. nf90_noerr)THEN CALL Gbcs_Critical( .TRUE., & 'Error getting level dim id',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inquire_dimension(ncid, LevDimID, len = nlev) IF(stat .ne. nf90_noerr)THEN CALL Gbcs_Critical( .TRUE., & 'Error getting nlevels length',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inq_dimid(ncid, "time", TimeDimID) IF(stat .ne. nf90_noerr)THEN CALL Gbcs_Critical( .TRUE., & 'Error getting time dim id',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inquire_dimension(ncid, TimeDimID, len = ntime) IF(stat .ne. nf90_noerr)THEN CALL Gbcs_Critical( .TRUE., & 'Error getting time length',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Get reference pressure stat = nf90_inq_varid(ncid,"P0",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for P0',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, P0) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting P0',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Get Lat/Lon/hyam,hyab IF( ALLOCATED(lon) )THEN DEALLOCATE(lon,lat,hyam,hybm) ENDIF ALLOCATE(lon(nlon),lat(nlat),hyam(nlev),hybm(nlev),STAT=stat) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error allocataing lon/lat/hyam/hybm arrays',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Longitude stat = nf90_inq_varid(ncid,"lon",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for lon',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, lon) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting lon',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Latitude stat = nf90_inq_varid(ncid,"lat",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for lat',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, lat) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting lat',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! hyam stat = nf90_inq_varid(ncid,"hyam",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for hyam',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, hyam) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting hyam',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! hybm stat = nf90_inq_varid(ncid,"hybm",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for hybm',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, hybm) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting hybm',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Now get time data and work out which elements we need for ! time IF( ALLOCATED(time) )THEN DEALLOCATE(time) ENDIF ALLOCATE(time(ntime),STAT=stat) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error allocataing time array',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inq_varid(ncid,"time",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for time',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, time) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting time',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Now convert time to get month - use data from same month Start_Time%Year = 1850 Start_Time%Month = 1 Start_Time%Day = 2 Start_Time%Hour = 0 Start_Time%Minute = 0 Start_Time%Seconds = 0 Start_Time%Sec1000 = 0 time_pos = -1 TimeLoop: DO I=1,nTime Time_Val = JD_to_Date(Time(I) + Date_to_JD(Start_Time)) IF( Time_Val%Month .eq. Image_Time%Month )THEN time_pos = I EXIT TimeLoop ENDIF END DO TimeLoop IF( -1 .eq. time_pos )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting position in array based on time',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Now get arrays ! First allocate everything we need IF( ALLOCATED(surface_pressure) )THEN DEALLOCATE(surface_pressure,temp,cb1,cb2,dst01,dst02,dst03,dst04,oc1,& oc2,so4,sslt01,sslt02,sslt03,sslt04) ENDIF ALLOCATE(surface_pressure(nlon,nlat,ntime),& temp(nlon,nlat,nlev,ntime),& cb1(nlon,nlat,nlev,ntime),& cb2(nlon,nlat,nlev,ntime),& dst01(nlon,nlat,nlev,ntime),& dst02(nlon,nlat,nlev,ntime),& dst03(nlon,nlat,nlev,ntime),& dst04(nlon,nlat,nlev,ntime),& oc1(nlon,nlat,nlev,ntime),& oc2(nlon,nlat,nlev,ntime),& so4(nlon,nlat,nlev,ntime),& sslt01(nlon,nlat,nlev,ntime),& sslt02(nlon,nlat,nlev,ntime),& sslt03(nlon,nlat,nlev,ntime),& sslt04(nlon,nlat,nlev,ntime),STAT=stat) IF(0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error allocating level arrays',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Actually get variables from NetCDF file ! Surface pressure stat = nf90_inq_varid(ncid,"PS",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for surface_pressure',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, surface_pressure) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting surface_pressure',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Temperature stat = nf90_inq_varid(ncid,"T",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for temperature',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, temp) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting temperature',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Aerosol data stat = nf90_inq_varid(ncid,"CB1",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for CB1',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, cb1) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting CB1',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inq_varid(ncid,"CB2",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for CB2',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, cb2) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting CB2',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inq_varid(ncid,"DST01",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for DST01',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, dst01) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting DST01',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inq_varid(ncid,"DST02",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for DST02',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, dst02) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting DST02',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inq_varid(ncid,"DST03",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for DST03',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, dst03) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting DST03',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inq_varid(ncid,"DST04",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for DST04',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, dst04) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting DST04',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inq_varid(ncid,"OC1",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for OC1',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, oc1) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting OC1',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inq_varid(ncid,"OC2",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for OC2',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, oc2) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting OC2',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inq_varid(ncid,"SO4",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for SO4',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, so4) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting SO4',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inq_varid(ncid,"SSLT01",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for SSLT01',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, sslt01) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting SSLT01',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inq_varid(ncid,"SSLT02",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for SSLT02',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, sslt02) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting SSLT02',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inq_varid(ncid,"SSLT03",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for SSLT03',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, sslt03) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting SSLT03',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_inq_varid(ncid,"SSLT04",varID) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting varID for SSLT04',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF stat = nf90_get_var(ncid, varID, sslt04) IF( 0 .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error getting SSLT04',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Close NETCDF file stat = nf90_close( ncid ) IF( NF90_NOERR .ne. stat )THEN CALL Gbcs_Critical( .TRUE., & 'Error closing NETCDF aerosol file',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! make Aerosol data - combine elements ! this part is taken from Prabhat Koners MATLAB code ! Get air density to convert aerosol data ! Then scale aerosol data DO I=1,nlev DO J=1,nlat DO K=1,nlon ! Note pressure is in Pa pressure_val = hyam(I)*P0 + hybm(I)*surface_pressure(K,J,time_pos) ! Units are kg/m^3 (28.97 is g/mol so 1000. converts to kg/mol) airdn = 28.97d0*pressure_val/& (8.314472d0*temp(K,J,I,time_pos)*1000.d0) ! Scale - only scale needed entries cb1(K,J,I,time_pos) = cb1(K,J,I,time_pos)*airdn cb2(K,J,I,time_pos) = cb2(K,J,I,time_pos)*airdn dst01(K,J,I,time_pos) = dst01(K,J,I,time_pos)*airdn dst02(K,J,I,time_pos) = dst02(K,J,I,time_pos)*airdn dst03(K,J,I,time_pos) = dst03(K,J,I,time_pos)*airdn dst04(K,J,I,time_pos) = dst04(K,J,I,time_pos)*airdn oc1(K,J,I,time_pos) = oc1(K,J,I,time_pos)*airdn oc2(K,J,I,time_pos) = oc2(K,J,I,time_pos)*airdn so4(K,J,I,time_pos) = so4(K,J,I,time_pos)*airdn sslt01(K,J,I,time_pos) = sslt01(K,J,I,time_pos)*airdn sslt02(K,J,I,time_pos) = sslt02(K,J,I,time_pos)*airdn sslt03(K,J,I,time_pos) = sslt03(K,J,I,time_pos)*airdn sslt04(K,J,I,time_pos) = sslt04(K,J,I,time_pos)*airdn END DO END DO END DO DEALLOCATE(temp) ! Bi-linear interpolation for spatial ! Setup GRID structure for Aerosol data ! Definition of grid Grid%Earth_Coord = .TRUE. Grid%RTM_Format = .FALSE. Grid%Full_Globe = .TRUE. Grid%NElems = nlon Grid%NLines = nlat Grid%NLevels = 1 Grid%Upwards = .FALSE. Grid%LLCrnr%lat = -90. Grid%LLCrnr%lon = 0. Grid%URCrnr%lat = 90. Grid%URCrnr%lon = 357.5 Grid%Incr%lat = 1.88542 Grid%Incr%lon = 2.5 Grid%I_LLCrnr%elem = 1 Grid%I_LLCrnr%line = 1 Grid%I_URCrnr%elem = nlon Grid%I_URCrnr%line = nlat ! Grid defined by positions at edge? Goes from 0,-90 Grid%DataLoc = 0 ! Surface Type - set to the same IF( ALLOCATED(surface_type) )THEN DEALLOCATE(surface_type) ENDIF ALLOCATE(surface_type(nlon,nlat),STAT=stat) IF( 0 .ne. stat )THEN CALL Gbcs_Critical(.TRUE.,'Cannot allocate surface_type',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF surface_type = 0 ! Interpolate onto model grid (levels) ! Allocate temp array for each type to hold interpolated values ! at Aerosol levels IF( ALLOCATED(temp_cb) )THEN DEALLOCATE(temp_cb,temp_dst,temp_oc,temp_so4,temp_sslt01,temp_sslt02,& temp_sslt03,temp_sslt04,temp_sfc_pressure,temp_dustr) ENDIF ALLOCATE(temp_cb(Model%Grid%NElems,Model%Grid%NLines,nlev),& temp_dst(Model%Grid%NElems,Model%Grid%NLines,nlev),& temp_oc(Model%Grid%NElems,Model%Grid%NLines,nlev),& temp_so4(Model%Grid%NElems,Model%Grid%NLines,nlev),& temp_sslt01(Model%Grid%NElems,Model%Grid%NLines,nlev),& temp_sslt02(Model%Grid%NElems,Model%Grid%NLines,nlev),& temp_sslt03(Model%Grid%NElems,Model%Grid%NLines,nlev),& temp_sslt04(Model%Grid%NElems,Model%Grid%NLines,nlev),& temp_sfc_pressure(Model%Grid%NElems,Model%Grid%NLines),& temp_dustr(Model%Grid%NElems,Model%Grid%NLines),& STAT=stat) IF( 0 .ne. stat )THEN CALL Gbcs_Critical(.TRUE.,'Cannot allocate temp arrays',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Do spatial interpolation ! Allocate memory so have correct array size to pass to InterpolateField ALLOCATE(temp_arr1(nlon,nlat),temp_arr2(nlon,nlat),temp_arr3(nlon,nlat),& temp_arr4(nlon,nlat),STAT=stat) ! This is done at all levels PX%Surface_Type = 0 DO I=1,Model%Grid%NLines PX%EarthCoord%lat = Model%Grid%LLCrnr%Lat + (I-1)*Grid%Incr%lat DO J=1,Model%Grid%NElems PX%EarthCoord%lon = Model%Grid%LLCrnr%Lon + (J-1)*Grid%Incr%lon CALL ModelIndices(PX,Grid,surface_type) temp_arr1 = surface_pressure(:,:,time_pos) temp_sfc_pressure(J,I) = InterpolateField(temp_arr1,PX%LocMP) dust_r = 0. dust_r_weight = 0. DO K=1,nlev temp_arr1 = cb1(:,:,K,time_pos) temp_arr2 = cb2(:,:,K,time_pos) temp_cb(J,I,K) = & InterpolateField(temp_arr1,PX%LocMP) + & InterpolateField(temp_arr2,PX%LocMP) temp_arr1 = dst01(:,:,K,time_pos) temp_arr2 = dst02(:,:,K,time_pos) temp_arr3 = dst03(:,:,K,time_pos) temp_arr4 = dst04(:,:,K,time_pos) temp_dst(J,I,K) = & InterpolateField(temp_arr1,PX%LocMP) + & InterpolateField(temp_arr2,PX%LocMP) + & InterpolateField(temp_arr3,PX%LocMP) + & InterpolateField(temp_arr4,PX%LocMP) ! Get dust radius dust_r = & 0.275*InterpolateField(temp_arr1,PX%LocMP)+& 0.875*InterpolateField(temp_arr2,PX%LocMP)+& 0.187*InterpolateField(temp_arr3,PX%LocMP)+& 3.75*InterpolateField(temp_arr4,PX%LocMP) dust_r_weight = & InterpolateField(temp_arr1,PX%LocMP)+& InterpolateField(temp_arr2,PX%LocMP)+& InterpolateField(temp_arr3,PX%LocMP)+& InterpolateField(temp_arr4,PX%LocMP) temp_arr1 = oc1(:,:,K,time_pos) temp_arr2 = oc2(:,:,K,time_pos) temp_oc(J,I,K) = & InterpolateField(temp_arr1,PX%LocMP) + & InterpolateField(temp_arr2,PX%LocMP) temp_arr1 = so4(:,:,K,time_pos) temp_so4(J,I,K) = InterpolateField(temp_arr1,PX%LocMP) temp_arr1 = sslt01(:,:,K,time_pos) temp_sslt01(J,I,K) = & InterpolateField(temp_arr1,PX%LocMP) temp_arr1 = sslt01(:,:,K,time_pos) temp_sslt02(J,I,K) = & InterpolateField(temp_arr1,PX%LocMP) temp_arr1 = sslt01(:,:,K,time_pos) temp_sslt03(J,I,K) = & InterpolateField(temp_arr1,PX%LocMP) temp_arr1 = sslt01(:,:,K,time_pos) temp_sslt04(J,I,K) = & InterpolateField(temp_arr1,PX%LocMP) END DO temp_dustr(J,I) = dust_r/dust_r_weight END DO END DO DEALLOCATE(temp_arr1,temp_arr2,temp_arr3,temp_arr4) DEALLOCATE(surface_type) DEALLOCATE(cb1,cb2,dst01,dst02,dst03,dst04,oc1,oc2,so4,sslt01,sslt02,& sslt03,sslt04,surface_pressure) ! Set model aerosol IDs Model%Aerosol_ID(1:8) = (/Gbcs_Aerosol_Dust,& Gbcs_Aerosol_SeaSalt_SAM,& Gbcs_Aerosol_SeaSalt_SSCM1,& Gbcs_Aerosol_SeaSalt_SSCM2,& Gbcs_Aerosol_SeaSalt_SSCM3,& Gbcs_Aerosol_Organic_Carbon,& Gbcs_Aerosol_Black_Carbon,& Gbcs_Aerosol_Sulfate/) ! Now have to interpolate to model levels ! Allocate memort for Aerosol data in model structure Aerosol_Indx = (/Gbcs_Map_Aerosol_Dust,& Gbcs_Map_Aerosol_SeaSalt_SAM,& Gbcs_Map_Aerosol_SeaSalt_SSCM1,& Gbcs_Map_Aerosol_SeaSalt_SSCM2,& Gbcs_Map_Aerosol_SeaSalt_SSCM3,& Gbcs_Map_Aerosol_Organic_Carbon,& Gbcs_Map_Aerosol_Black_Carbon,& Gbcs_Map_Aerosol_Sulfate/) DO I=1,8 IF( ALLOCATED(Model%Atm(I)%d) )THEN DEALLOCATE(Model%Atm(I)%d) ENDIF ALLOCATE(Model%Atm(I)%d(Model%Grid%NLevels,Model%Grid%NElems,& Model%Grid%NLines),STAT=stat) IF( 0 .ne. stat )THEN CALL Gbcs_Critical(.TRUE.,'Cannot allocate model aerosol arrays',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END DO ! Have to do this for all profiles individually since pressure levels ! are not fixed IF( ALLOCATED(Profile%P) )THEN DEALLOCATE(Profile%P) ENDIF ALLOCATE(Profile%P(nlev),STAT=stat) IF( 0 .ne. stat )THEN CALL Gbcs_Critical(.TRUE.,'Cannot allocate Profilt%P array',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF Profile%Nlevels = nlev ALLOCATE(press(Model%Grid%Nlevels),prof(nlev),int_lev(nlev),& int_fact(nlev),STAT=stat) IF( 0 .ne. stat )THEN CALL Gbcs_Critical(.TRUE.,'Cannot allocate prof array',& 'Get_Aerosol_from_Clim', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF DO I=1,Model%Grid%NLines DO J=1,Model%Grid%NElems ! Set Aerosol radius's Model%Srf(Gbcs_Map_Aerosol_Black_Carbon)%d(J,I) = 0.05 Model%Srf(Gbcs_Map_Aerosol_Dust)%d(J,I) = temp_dustr(J,I) Model%Srf(Gbcs_Map_Aerosol_Organic_Carbon)%d(J,I) = 0.15 Model%Srf(Gbcs_Map_Aerosol_Sulfate)%d(J,I) = 0.52 Model%Srf(Gbcs_Map_Aerosol_SeaSalt_SAM)%d(J,I) = 0.6 Model%Srf(Gbcs_Map_Aerosol_SeaSalt_SSCM1)%d(J,I) = 1.5 Model%Srf(Gbcs_Map_Aerosol_SeaSalt_SSCM2)%d(J,I) = 3.25 Model%Srf(Gbcs_Map_Aerosol_SeaSalt_SSCM3)%d(J,I) = 7.5 ! Make profile ! Pressure based on interplated surface pressure Profile%P = hyam*P0 + hybm*temp_sfc_pressure(J,I) ! Set 2m pressure to lowest pressure value ! 2m pressure is used by Get_Interpol_Weights Profile%P_2m = Profile%P(nlev) ! Get weights press = Model%Atm(Gbcs_Map_P_3d)%d(:,J,I) CALL Gen_Interpol_Weights( Profile, Model%Grid%NLevels, & press,int_lev,int_fact,int_srf) ! Now interpolate aerosol data prof = temp_cb(J,I,:) prof_2m = temp_cb(J,I,nlev) Model%Atm(Gbcs_Map_Aerosol_Black_Carbon)%d(:,J,I) = & Apply_Interpol_Weights(prof,prof_2m,& Model%Grid%NLevels,int_lev,int_fact,int_srf) prof = temp_dst(J,I,:) prof_2m = temp_dst(J,I,nlev) Model%Atm(Gbcs_Map_Aerosol_Dust)%d(:,J,I) = & Apply_Interpol_Weights(prof,prof_2m,& Model%Grid%NLevels,int_lev,int_fact,int_srf) prof = temp_oc(J,I,:) prof_2m = temp_oc(J,I,nlev) Model%Atm(Gbcs_Map_Aerosol_Organic_Carbon)%d(:,J,I) = & Apply_Interpol_Weights(prof,prof_2m,& Model%Grid%NLevels,int_lev,int_fact,int_srf) prof = temp_so4(J,I,:) prof_2m = temp_so4(J,I,nlev) Model%Atm(Gbcs_Map_Aerosol_Sulfate)%d(:,J,I) = & Apply_Interpol_Weights(prof,prof_2m,& Model%Grid%NLevels,int_lev,int_fact,int_srf) prof = temp_sslt01(J,I,:) prof_2m = temp_sslt01(J,I,nlev) Model%Atm(Gbcs_Map_Aerosol_SeaSalt_SAM)%d(:,J,I) = & Apply_Interpol_Weights(prof,prof_2m,Model%Grid%NLevels,& int_lev,int_fact,int_srf) prof = temp_sslt02(J,I,:) prof_2m = temp_sslt02(J,I,nlev) Model%Atm(Gbcs_Map_Aerosol_SeaSalt_SSCM1)%d(:,J,I) = & Apply_Interpol_Weights(prof,prof_2m,Model%Grid%NLevels,& int_lev,int_fact,int_srf) prof = temp_sslt03(J,I,:) prof_2m = temp_sslt03(J,I,nlev) Model%Atm(Gbcs_Map_Aerosol_SeaSalt_SSCM2)%d(:,J,I) = & Apply_Interpol_Weights(prof,prof_2m,Model%Grid%NLevels,& int_lev,int_fact,int_srf) prof = temp_sslt04(J,I,:) prof_2m = temp_sslt04(J,I,nlev) Model%Atm(Gbcs_Map_Aerosol_SeaSalt_SSCM3)%d(:,J,I) = & Apply_Interpol_Weights(prof,prof_2m,Model%Grid%NLevels,& int_lev,int_fact,int_srf) END DO END DO DEALLOCATE(press,prof,int_lev,int_fact) ! Register fields STAT = Register_3D_Field(Model, Gbcs_Map_Aerosol_Dust, & Gbcs_Units_kg_per_m3) STAT = Register_3D_Field(Model, Gbcs_Map_Aerosol_SeaSalt_SAM, & Gbcs_Units_kg_per_m3) STAT = Register_3D_Field(Model, Gbcs_Map_Aerosol_SeaSalt_SSCM1, & Gbcs_Units_kg_per_m3) STAT = Register_3D_Field(Model, Gbcs_Map_Aerosol_SeaSalt_SSCM2, & Gbcs_Units_kg_per_m3) STAT = Register_3D_Field(Model, Gbcs_Map_Aerosol_SeaSalt_SSCM3, & Gbcs_Units_kg_per_m3) STAT = Register_3D_Field(Model, Gbcs_Map_Aerosol_Organic_Carbon, & Gbcs_Units_kg_per_m3) STAT = Register_3D_Field(Model, Gbcs_Map_Aerosol_Black_Carbon, & Gbcs_Units_kg_per_m3) STAT = Register_3D_Field(Model, Gbcs_Map_Aerosol_Sulfate, & Gbcs_Units_kg_per_m3) STAT = Register_2D_Field(Model, Gbcs_Map_Aerosol_Dust_Rad, & Gbcs_Units_micron) STAT = Register_2D_Field(Model, Gbcs_Map_Aerosol_SeaSalt_SAM_Rad, & Gbcs_Units_micron) STAT = Register_2D_Field(Model, Gbcs_Map_Aerosol_SeaSalt_SSCM1_Rad, & Gbcs_Units_micron) STAT = Register_2D_Field(Model, Gbcs_Map_Aerosol_SeaSalt_SSCM2_Rad, & Gbcs_Units_micron) STAT = Register_2D_Field(Model, Gbcs_Map_Aerosol_SeaSalt_SSCM3_Rad, & Gbcs_Units_micron) STAT = Register_2D_Field(Model, Gbcs_Map_Aerosol_Organic_Carbon_Rad, & Gbcs_Units_micron) STAT = Register_2D_Field(Model, Gbcs_Map_Aerosol_Black_Carbon_Rad, & Gbcs_Units_micron) STAT = Register_2D_Field(Model, Gbcs_Map_Aerosol_Sulfate_Rad, & Gbcs_Units_micron) ! Deallocate DEALLOCATE(lon,lat,hyam,hybm,time,surface_pressure) DEALLOCATE(temp_cb,temp_dst,temp_oc,temp_so4,temp_sslt01,temp_sslt02,& temp_sslt03,temp_sslt04,temp_dustr,temp_sfc_pressure) END SUBROUTINE Get_Aerosol_from_Clim #endif ! Reset ocean pixels using high res sst and land mask ! Also set land pixels nearest ocean to ocean values if ocean onty ! retrievals are required SUBROUTINE Reset_Model_Surface_Temp( AUX, IMG ) USE MCIDAS_LoadImagery TYPE(Aux_Data), INTENT(INOUT) :: AUX TYPE(Imagery), INTENT(INOUT) :: IMG ! Local variables INTEGER :: I, J, K INTEGER :: ilat, ilon INTEGER :: index INTEGER :: Xpos, Ypos INTEGER :: STAT LOGICAL, ALLOCATABLE :: Reset_To_Ocean(:,:) ! Reset land near ocean - logical to tell you when it's been done ALLOCATE(Reset_To_Ocean(AUX%Model%Grid%Nelems,AUX%Model%Grid%Nlines),& STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical(.TRUE.,'Cannot allocate Reset_To_Ocean array',& 'Reset_Model_Surface_Type', 'NOAA/NOAA_LoadAuxData.f90') ENDIF Reset_to_Ocean = .FALSE. ! Landmask should be read in but if not (determined by individual image ! read routines ! make sure we don't re-do the image landmask if we are reading it in index = Get_SurfTemp_Index(AUX%Model, Gbcs_Surface_Sea) IF( -1 .eq. index )THEN call Gbcs_Critical(.TRUE.,& 'No AUX%Model sea temperature avail',& 'Reset_Model_Surface_Temp',& 'NOAA/NOAA_LoadAuxData.f90') ENDIF ! Read in each half/segment (lon based) to reduce memory load DO K=1,8 CALL Read_In_Half_Landmask( IMG%GbcsDataPath, K, 8 ) ! First reset model surface temperatures using high res sst DO I=1,AUX%Model%Grid%Nlines DO J=1,AUX%Model%Grid%Nelems ! If not in right half IF( -1 .eq. HR_LandMask(AUX%Model%Lon(J,I),AUX%Model%Lat(J,I),& use_half=.TRUE.) )THEN CYCLE ENDIF IF( ALLOCATED(AUX%SST_Clim%SST) )THEN IF( Gbcs_Surface_Sea .eq. & HR_Landmask(AUX%Model%Lon(J,I),AUX%Model%Lat(J,I)) .and. & Gbcs_Surface_Sea .ne. AUX%Model%Surface_type(J,I) )THEN ! Get high res SST location IF( 0 .eq. AUX%SST_Clim%Grid%DataLoc )THEN ilat = FLOOR( ( AUX%Model%Lat(J,I) - & AUX%SST_Clim%Grid%LLCrnr%lat )& / AUX%SST_Clim%Grid%Incr%lat ) + & AUX%SST_Clim%Grid%I_LLCrnr%line ELSE ilat = FLOOR( ( AUX%Model%Lat(J,I) - & AUX%SST_Clim%Grid%LLCrnr%lat - & 0.5*AUX%SST_Clim%Grid%Incr%lat )& / AUX%SST_Clim%Grid%Incr%lat ) + & AUX%SST_Clim%Grid%I_LLCrnr%line ENDIF IF ( ilat < 1 ) ilat = 1 IF ( ilat > AUX%SST_Clim%Grid%NLines ) ilat = & AUX%SST_Clim%Grid%NLines IF( 0 .eq. AUX%SST_Clim%Grid%DataLoc )THEN ilon = FLOOR( ( AUX%Model%Lon(J,I) - & AUX%SST_Clim%Grid%LLCrnr%lon )& / AUX%SST_Clim%Grid%Incr%lon ) + & AUX%SST_Clim%Grid%I_LLCrnr%elem ELSE ilon = FLOOR( ( AUX%Model%Lon(J,I) - & AUX%SST_Clim%Grid%LLCrnr%lon - & 0.5*AUX%SST_Clim%Grid%Incr%lon)& / AUX%SST_Clim%Grid%Incr%lon ) + & AUX%SST_Clim%Grid%I_LLCrnr%elem ENDIF IF ( ilon < 1 ) ilon = ilon + AUX%SST_Clim%Grid%NElems IF ( ilon > AUX%SST_Clim%Grid%NElems ) ilon = ilon - & AUX%SST_Clim%Grid%NElems ! If NODC Climate data then landmask looks to be ! slightly different so if over land use NCEP model data ! i.e. do nothing IF( 'NODC Pathfinder climatology (Casey & Cornillon 1999)' & .eq. AUX%SST_Clim%Name )THEN IF( 271.3 .le. AUX%SST_Clim%SST(ilon,ilat,1) )THEN AUX%Model%Srf(index)%d(J,I) = & AUX%SST_Clim%SST(ilon,ilat,1) ELSE ! Find nearest ocean within climatology ! -1 means no limit to search IF( Get_Nearest_Ocean(AUX,-1,ilon,ilat,XPos,YPos,& use_clim=.TRUE.) )THEN AUX%Model%Srf(index)%d(J,I) = & AUX%SST_Clim%SST(XPos,YPos,1) ENDIF ENDIF ELSE IF( 'NOAA RTG data file' .eq. AUX%SST_Clim%Name )THEN ! RTG - has data over whole globe AUX%Model%Srf(index)%d(J,I) = & AUX%SST_Clim%SST(ilon,ilat,1) ELSE IF( 'NOAA Reynolds OI v2' .eq. AUX%SST_Clim%Name )THEN ! Reynolds - has only over ocean IF( 271.3 .le. AUX%SST_Clim%SST(ilon,ilat,1) )THEN AUX%Model%Srf(index)%d(J,I) = & AUX%SST_Clim%SST(ilon,ilat,1) ELSE ! Find nearest ocean within climatology ! -1 means no limit to search IF( Get_Nearest_Ocean(AUX,-1,ilon,ilat,XPos,YPos,& use_clim=.TRUE.) )THEN AUX%Model%Srf(index)%d(J,I) = & AUX%SST_Clim%SST(XPos,YPos,1) ENDIF ENDIF ELSE CALL Gbcs_Critical(.TRUE., & 'No matching high res data to be applied', & 'Reset_Model_Surface_Type', & 'NOAA/NOAA_LoadAuxData.f90') ENDIF Reset_To_Ocean(J,I) = .TRUE. ENDIF ELSE IF( Gbcs_Surface_Sea .eq. & HR_Landmask(AUX%Model%Lon(J,I),AUX%Model%Lat(J,I)) .and. & Gbcs_Surface_Sea .ne. AUX%Model%Surface_type(J,I) )THEN ! Find nearest ocean within 3 pixels IF( Get_Nearest_Ocean(AUX,3,J,I,Xpos,Ypos) )THEN AUX%Model%Srf(index)%d(J,I) = & AUX%Model%Srf(index)%d(XPos,YPos) Reset_To_Ocean(J,I) = .TRUE. ENDIF ENDIF ENDIF END DO END DO END DO ! Set flags after first run DO I=1,AUX%Model%Grid%Nlines DO J=1,AUX%Model%Grid%Nelems IF( Reset_To_Ocean(J,I) )THEN AUX%Model%Surface_Type(J,I) = Gbcs_Surface_Sea ENDIF END DO END DO ! WRITE(101,*)AUX%Model%Grid%Nelems,AUX%Model%Grid%Nlines ! WRITE(101,*)AUX%Model%Surface_Type ! STOP ! This part sets all land near ocean to ocean Reset_To_Ocean = .FALSE. DO I=1,AUX%Model%Grid%Nlines DO J=1,AUX%Model%Grid%Nelems ! Reset only if next to ocean so maxstep = 1 IF( Get_Nearest_Ocean(AUX,1,J,I,Xpos,Ypos) )THEN AUX%Model%Srf(index)%d(J,I) = & AUX%Model%Srf(index)%d(XPos,YPos) Reset_To_Ocean(J,I) = .TRUE. ENDIF END DO END DO ! Set flags for final time DO I=1,AUX%Model%Grid%Nlines DO J=1,AUX%Model%Grid%Nelems IF( Reset_To_Ocean(J,I) )THEN AUX%Model%Surface_Type(J,I) = Gbcs_Surface_Sea ENDIF END DO END DO DEALLOCATE( Reset_To_Ocean ) END SUBROUTINE Reset_Model_Surface_Temp ! Function to return position of nearest ocean pixel on model to point LOGICAL FUNCTION Get_Nearest_Ocean( AUX, MaxStepIn, I, J, & XposOut, YposOut, use_clim )RESULT(Found) TYPE(Aux_Data), INTENT(IN) :: AUX INTEGER, INTENT(IN) :: MaxStepIn INTEGER, INTENT(IN) :: I, J INTEGER, INTENT(OUT) :: XposOut, YposOut LOGICAL, INTENT(IN), OPTIONAL :: use_clim ! Local Variables INTEGER :: Step INTEGER :: MaxStep LOGICAL :: use_climatology INTEGER :: nlines, nelems INTEGER :: XPos, YPos INTEGER :: XMin, XMax INTEGER :: YMin, YMax Found = .FALSE. IF( PRESENT(use_clim) )THEN use_climatology = use_clim ELSE use_climatology = .FALSE. ENDIF IF( use_climatology )THEN nelems = AUX%SST_Clim%Grid%Nelems nlines = AUX%SST_Clim%Grid%Nlines ELSE ! If we are on a sea pixel return IF( Gbcs_Surface_Sea .eq. & AUX%Model%Surface_Type(I,J) )THEN RETURN ENDIF nelems = AUX%Model%Grid%Nelems nlines = AUX%Model%Grid%Nlines ENDIF ! If no limit on search IF( -1 .eq. MaxStepIn )THEN MaxStep = MIN(nelems,nlines) ELSE MaxStep = MaxStepIn ENDIF ! Break up of loop means it's faster for large MaxStep's doLoop: DO Step=1,MaxStep ! Constant Y edges Ypos=J-Step IF( 1 .le. Ypos )THEN XMin=I-Step IF( 0 .ge. XMin)XMin = 1 XMax=I+Step IF( Nelems .lt. XMax)XMax = Nelems DO XPos=XMin,XMax IF( use_climatology )THEN IF( 271.3 .le. AUX%SST_Clim%SST(XPos,YPos,1) )THEN XPosOut = XPos YPosOut = YPos Found = .TRUE. EXIT DoLoop ENDIF ELSE IF( Gbcs_Surface_Sea .eq. & AUX%Model%Surface_Type(XPos,Ypos) )THEN XPosOut = XPos YPosOut = YPos Found = .TRUE. EXIT DoLoop ENDIF ENDIF END DO ENDIF Ypos=J+Step IF( Nlines .ge. Ypos )THEN XMin=I-Step IF( 0 .ge. XMin)XMin = 1 XMax=I+Step IF( Nelems .lt. XMax)XMax = Nelems DO XPos=XMin,XMax IF( use_climatology )THEN IF( 271.3 .le. AUX%SST_Clim%SST(XPos,YPos,1) )THEN XPosOut = XPos YPosOut = YPos Found = .TRUE. EXIT DoLoop ENDIF ELSE IF( Gbcs_Surface_Sea .eq. & AUX%Model%Surface_Type(XPos,Ypos) )THEN XPosOut = XPos YPosOut = YPos Found = .TRUE. EXIT DoLoop ENDIF ENDIF END DO ENDIF ! Constant X edges Xpos=I-Step IF( 1 .le. Xpos )THEN YMin=J-Step IF( 0 .ge. YMin)YMin = 1 YMax=J+Step IF( Nlines .lt. YMax)YMax = Nlines DO YPos=YMin,YMax IF( use_climatology )THEN IF( 271.3 .le. AUX%SST_Clim%SST(XPos,YPos,1) )THEN XPosOut = XPos YPosOut = YPos Found = .TRUE. EXIT DoLoop ENDIF ELSE IF( Gbcs_Surface_Sea .eq. & AUX%Model%Surface_Type(XPos,Ypos) )THEN XPosOut = XPos YPosOut = YPos Found = .TRUE. EXIT DoLoop ENDIF ENDIF END DO ENDIF Xpos=I+Step IF( Nelems .ge. Xpos )THEN YMin=J-Step IF( 0 .ge. YMin)YMin = 1 YMax=J+Step IF( Nlines .lt. YMax)YMax = Nlines DO YPos=YMin,YMax IF( use_climatology )THEN IF( 271.3 .le. AUX%SST_Clim%SST(XPos,YPos,1) )THEN XPosOut = XPos YPosOut = YPos Found = .TRUE. EXIT DoLoop ENDIF ELSE IF( Gbcs_Surface_Sea .eq. & AUX%Model%Surface_Type(XPos,Ypos) )THEN XPosOut = XPos YPosOut = YPos Found = .TRUE. EXIT DoLoop ENDIF ENDIF END DO ENDIF END DO doLoop END FUNCTION Get_Nearest_Ocean ! Read Met Office ATOVS Covariance Matrix ! Read in NOAA version of load covariance matrix INTEGER FUNCTION NOAA_Load_Cov_Mat( Cov_Mat_Path , Cov_Mat_File , Cov_Mat ) ! ! Description: ! ! Loads the covariance matrix including extra flag (if present) to ! deal with multaplicative covariances (not dealt with in standard ! GBCS code) ! ! Method: ! ! Owner: Manager of TRICS Project ! ! History: ! Version Date Comment ! ------- ---- ------- ! 0.0 14/12/2009 Creation JPDM ! ! ! Code Description: ! Language: Fortran 90 ! ! ! Modules used: USE GbcsTypes USE GbcsSystemTools, ONLY: Get_New_File_Unit USE GbcsErrorHandler, ONLY: Gbcs_Log_File_Unit IMPLICIT NONE ! Define Input/Output Variabiles CHARACTER(LEN=*), INTENT (IN ) :: Cov_Mat_Path ! Covariance matrix data path CHARACTER(LEN=*), INTENT (IN ) :: Cov_Mat_File ! Covariance matrix data file TYPE(Cov_Matrix), INTENT (INOUT) :: Cov_Mat ! Covariance matrix structure ! Define Local Variabiles INTEGER :: matrixsize, nbands, nfields INTEGER :: field,nlevels,itype INTEGER :: nvars INTEGER :: nband REAL :: lat1,lat2 INTEGER :: STAT INTEGER :: ii, jj CHARACTER :: Str LOGICAL :: File_Exists INTEGER :: File_Unit CHARACTER(LEN=256) :: File_Name CHARACTER(LEN=256) :: inLine CHARACTER(LEN=50), PARAMETER :: Routine_Name = 'Load_Cov_Mat' nvars = 0 File_Unit = Get_New_File_Unit() File_Name = TRIM(Cov_Mat_Path)//TRIM(Cov_Mat_File) INQUIRE( FILE=TRIM( File_Name ) , EXIST=File_Exists ) IF ( File_Exists ) THEN OPEN( UNIT=File_Unit , FILE=File_Name , FORM="FORMATTED" ) DO ii = 1 , 17 READ( UNIT=File_Unit , FMT=* ) Str END DO READ( UNIT=File_Unit , FMT=* ) matrixsize , nbands , nfields Cov_Mat%NBands = nbands Cov_Mat%MSize = matrixsize Cov_Mat%NFields = nfields ! Read in field, nlevels and type (1=constant, 2=mutliplicative) DO ii = 1 , nfields ! READ( UNIT=File_Unit , FMT=* ) field , nlevels READ( UNIT=File_Unit, FMT='(a)' )inLine READ( inLine, FMT=*, IOSTAT=STAT ) field, nlevels, itype IF( 0 .eq. STAT )THEN Cov_Mat%Field_ID( ii ) = field Cov_Mat%NLevels( ii ) = nlevels Cov_Mat%Type( ii ) = itype ELSE READ( inLine, FMT=*, IOSTAT=STAT ) field, nlevels IF( 0 .ne. STAT )THEN CALL Gbcs_Critical(.TRUE.,'Error reading line (COV MATRIX)',& 'NOAA_Load_Cov_Mat','NOAA/NOAA_LoadAuxData.f90') ENDIF Cov_Mat%Field_ID( ii ) = field Cov_Mat%NLevels( ii ) = nlevels Cov_Mat%Type( ii ) = 1 ! Just a constant value ENDIF nvars = nvars + nlevels END DO IF ( ALLOCATED( Cov_Mat%Data ) ) DEALLOCATE( Cov_Mat%Data , STAT=STAT ) ALLOCATE( Cov_Mat%Data( matrixsize , nvars , nbands ) , STAT=STAT ) DO jj = 1 , nbands READ( UNIT=File_Unit , FMT=* ) nband , lat1 , lat2 Cov_Mat%Limits( nband , 1 ) = lat1 Cov_Mat%Limits( nband , 2 ) = lat2 READ( UNIT=File_Unit , FMT=* ) Cov_Mat%Data( 1:matrixsize , 1:nvars , jj ) END DO CLOSE( File_Unit ) NOAA_Load_Cov_Mat = STAT ELSE WRITE(Gbcs_Log_File_Unit,*)' ERROR : File does not exist' WRITE(Gbcs_Log_File_Unit,*)' ',TRIM(File_Name) WRITE(Gbcs_Log_File_Unit,*)' SUBROUTINE Load_Cov_Matrix' WRITE(Gbcs_Log_File_Unit,*)' CODE - /GbcsPrg_DataReaders/GbcsCovMatLoaders.f90' STOP END IF RETURN END FUNCTION NOAA_Load_Cov_Mat !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Routines to read in a climatology (SST) #ifdef USE_NETCDF SUBROUTINE Read_Reynolds_SST(filename,nx,ny,SST) CHARACTER(LEN=*), INTENT(IN) :: filename INTEGER, INTENT(IN) :: nx INTEGER, INTENT(IN) :: ny REAL(GbcsReal), INTENT(INOUT) :: SST(nx,ny) INTEGER(GbcsInt2), ALLOCATABLE :: inSST(:,:,:,:) INTEGER :: I,J INTEGER :: ncid INTEGER :: dim_id INTEGER :: var_id INTEGER :: nlon INTEGER :: nlat INTEGER :: status ! Open file status = nf90_open( path=filename, mode=nf90_nowrite, ncid=ncid) if (status /= nf90_noerr) then call Gbcs_Critical(.TRUE.,'File open error','Read_Reynolds_SST',& 'NOAA/NOAA_LoadAuxData.f90') endif ! Get dimension ID for longitude status = nf90_inq_dimid(ncid, "lon", dim_id) if (status /= nf90_noerr) then call Gbcs_Critical(.TRUE.,'Getting lon dimension ID',& 'Read_Reynolds_SST','NOAA/NOAA_LoadAuxData.f90') endif ! Get size status = nf90_inquire_dimension(ncid, dim_id, len = nlon) if (status /= nf90_noerr) then call Gbcs_Critical(.TRUE.,'Getting nlon',& 'Read_Reynolds_SST','NOAA/NOAA_LoadAuxData.f90') endif ! Get dimension ID for latitude status = nf90_inq_dimid(ncid, "lat", dim_id) if (status /= nf90_noerr) then call Gbcs_Critical(.TRUE.,'Getting lat dimension ID',& 'Read_Reynolds_SST','NOAA/NOAA_LoadAuxData.f90') endif ! Get size status = nf90_inquire_dimension(ncid, dim_id, len = nlat) if (status /= nf90_noerr) then call Gbcs_Critical(.TRUE.,'Getting nlat',& 'Read_Reynolds_SST','NOAA/NOAA_LoadAuxData.f90') endif IF( nx .ne. nlon .or. ny .ne. nlat )THEN call Gbcs_Critical(.TRUE.,& 'Mismatch between input dims and netcdf dims',& 'Read_Reynolds_SST','NOAA/NOAA_LoadAuxData.f90') ENDIF ! Allocate 2 byte int array to read SST ALLOCATE(inSST(nlon,nlat,1,1),STAT=status) IF( 0 .ne. status )THEN call Gbcs_Critical(.TRUE.,& 'Cannot allocate inSST',& 'Read_Reynolds_SST','NOAA/NOAA_LoadAuxData.f90') ENDIF ! Get SST status = nf90_inq_varid(ncid, "sst", var_id) if (status /= nf90_noerr) then call Gbcs_Critical(.TRUE.,& 'Getting SST array (var_id)',& 'Read_Reynolds_SST','NOAA/NOAA_LoadAuxData.f90') endif status = nf90_get_var(ncid, var_id, inSST ) if (status /= nf90_noerr) then call Gbcs_Critical(.TRUE.,& 'Getting SST array (read)',& 'Read_Reynolds_SST','NOAA/NOAA_LoadAuxData.f90') endif ! Write SST to output array DO I=1,nlat DO J=1,nlon SST(J,I) = inSST(J,I,1,1)*0.01+273.15 END DO END DO ! Deallocate read array DEALLOCATE(inSST) ! Close file status = nf90_close(ncid=ncid) if (status /= nf90_noerr) then call Gbcs_Critical(.TRUE.,'Closing file',& 'Read_Reynolds_SST','NOAA/NOAA_LoadAuxData.f90') endif END SUBROUTINE Read_Reynolds_SST #endif ! Read in and renormalize SST Climatology from Pathfinder 9km climatology SUBROUTINE Load_SST_Climatology( IMG, AUX, SST_Clim, Time, use_other_anal ) TYPE(Imagery), INTENT(IN) :: IMG TYPE(Aux_Data), INTENT(IN) :: AUX TYPE(High_Res_SST_Struct), INTENT(OUT) :: SST_Clim TYPE(DateTime), INTENT(IN) :: Time LOGICAL, OPTIONAL :: use_other_anal INTEGER :: Unit INTEGER :: Ios CHARACTER(LEN=256) :: climate_filename LOGICAL :: exists CHARACTER(LEN=200) :: L2P_DIR CHARACTER(LEN=200) :: rtg_filename CHARACTER(LEN=200) :: reynolds_filename CHARACTER(LEN=200) :: command TYPE(DateTime) :: inTime TYPE(DateTime) :: New_Time INTEGER, PARAMETER :: nx_climate=4096 INTEGER, PARAMETER :: ny_climate=2048 INTEGER(GbcsInt1) :: climate_sst(nx_climate,ny_climate) INTEGER, PARAMETER :: nx_rtg=4320 INTEGER, PARAMETER :: ny_rtg=2160 REAL(GbcsReal) :: rtg_sst(nx_rtg,ny_rtg) LOGICAL :: Found_RTG INTEGER, PARAMETER :: nx_reynolds=1440 INTEGER, PARAMETER :: ny_reynolds=720 INTEGER :: I,J INTEGER :: II,JJ REAL(GbcsReal), ALLOCATABLE :: climate_data(:,:) INTEGER :: STAT INTEGER(GbcsInt2) :: climate_sst_value LOGICAL :: USE_RTG USE_RTG = .TRUE. IF( PRESENT(use_other_anal) )THEN USE_RTG = .not.use_other_anal ENDIF Found_RTG = .TRUE. ! Now try and read in latest RTG data/Analysis #ifdef USE_NETCDF ! Reynolds is stored in NetCDF format CALL GETENV('SST_REYNOLDS_DAILY_DIR',L2P_DIR) IF( ' ' .ne. L2P_DIR )THEN ! Read in REYNOLDS Daily SST ! make filename WRITE(rtg_filename,& '(a,''/'',i4.4,''/amsr-avhrr-v2.'',i4.4,i2.2,i2.2,''.nc.gz'')')& TRIM(L2P_DIR),IMG%Mean_Time%Year,IMG%Mean_Time%Year,& IMG%Mean_Time%Month,IMG%Mean_Time%Day INQUIRE(File=rtg_filename,EXIST=exists) IF( .not.exists )THEN CALL Gbcs_Critical(.TRUE.,'Cannot find Reynolds Daily SST',& 'Load_SST_Climatology','NOAA/NOAA_LoadAuxData.f90') ENDIF ! Get uncompressed version WRITE(reynolds_filename,'(''reynolds.'',a,''.nc'')')& TRIM(Get_Unique_String()) WRITE(command,'(''cp -f '',a,'' '',a,''.gz'')')& TRIM(rtg_filename),TRIM(reynolds_filename) CALL SYSTEM(command) WRITE(command,'(''gunzip -f '',a,''.gz'')')TRIM(reynolds_filename) CALL SYSTEM(command) IF( ALLOCATED(SST_Clim%SST) )THEN DEALLOCATE(SST_Clim%SST) ENDIF write(6,*) 'DEBUG------------------------------------------' write(6,*) ' System status before allocating SST_Clim%SST' call system(' free -m ') write(6,*) 'DEBUG------------------------------------------' ALLOCATE(SST_Clim%SST(nx_reynolds,ny_reynolds,1),STAT=STAT) IF( 0 /= STAT )THEN write(6,*) 'DEBUG------------------------------------------' write(6,*) ' Error return from allocate: ',STAT call system(' free -m ') write(6,*) 'DEBUG------------------------------------------' CALL Gbcs_Critical(.TRUE.,& 'Error allocating memory (SST_Clim%SST)',& 'Load_SST_Climatology','NOAA/NOAA_LoadAuxData') ENDIF write(6,*) 'DEBUG------------------------------------------' write(6,*) ' System status after allocating SST_Clim%SST' call system(' free -m ') write(6,*) 'DEBUG------------------------------------------' CALL Read_Reynolds_SST(TRIM(reynolds_filename),nx_reynolds,ny_reynolds,& SST_Clim%SST) WRITE(command,'(''rm -f '',a)')TRIM(reynolds_filename) CALL SYSTEM(command) ! Name SST_Clim%Name = 'NOAA Reynolds OI v2' ! Units SST_Clim%Units = Gbcs_Units_Kelvin ! Definition of grid SST_Clim%Grid%Earth_Coord = .TRUE. SST_Clim%Grid%RTM_Format = .FALSE. SST_Clim%Grid%Full_Globe = .TRUE. SST_Clim%Grid%NElems = nx_reynolds SST_Clim%Grid%NLines = ny_reynolds SST_Clim%Grid%NLevels = 1 SST_Clim%Grid%Upwards = .FALSE. SST_Clim%Grid%LLCrnr%lat = -89.875 SST_Clim%Grid%LLCrnr%lon = 0.125 SST_Clim%Grid%URCrnr%lat = 89.875 SST_Clim%Grid%URCrnr%lon = 359.875 SST_Clim%Grid%Incr%lat = 0.25 SST_Clim%Grid%Incr%lon = 0.25 SST_Clim%Grid%I_LLCrnr%elem = 1 SST_Clim%Grid%I_LLCrnr%line = 1 SST_Clim%Grid%I_URCrnr%elem = nx_reynolds SST_Clim%Grid%I_URCrnr%line = ny_reynolds SST_Clim%Grid%DataLoc = 1 write(Gbcs_Log_File_Unit,& '('' USING Reynolds AVHRR-AMSR OI SST'')') RETURN ENDIF #endif CALL GETENV('SST_ANAL',L2P_DIR) IF( ' ' .eq. L2P_DIR )THEN CALL Gbcs_Warning(.TRUE.,& 'SST_ANAL environment variable is not defined - using Pathfinder climatology',& 'Load_SST_Climatology','NOAA/NOAA_LoadAuxData') Found_RTG = .FALSE. ELSE IF( USE_RTG )THEN L2P_DIR = TRIM(L2P_DIR)//'/rtg/' ! Look for current date version - not always back 1 day CALL Time_Back_One_Day( Time, New_Time ) write(rtg_filename,& '(''rtg_sst_grb_hr_0.083.'',i4.4,i2.2,i2.2,''.dat'')')& New_Time%Year,New_Time%Month,New_Time%Day climate_filename = TRIM(L2P_DIR)//TRIM(rtg_filename) INQUIRE(File=climate_filename,EXIST=exists) IF( .not.exists )THEN inTime = New_Time CALL Gbcs_Warning(.TRUE.,& 'Curent data RTG file does not exist - looking back 1 day',& 'Load_SST_Climatology','NOAA/NOAA_LoadAuxData') ! Go back 1 day CALL Time_Back_One_Day( inTime, New_Time ) write(rtg_filename,& '(''rtg_sst_grb_hr_0.083.'',i4.4,i2.2,i2.2,''.dat'')')& New_Time%Year,New_Time%Month,New_Time%Day climate_filename = TRIM(L2P_DIR)//TRIM(rtg_filename) INQUIRE(File=climate_filename,EXIST=exists) IF( .not.exists )THEN CALL Gbcs_Warning(.TRUE.,& 'Cannot find RTG - not using high res. SST',& 'Load_SST_Climatology','NOAA/NOAA_LoadAuxData') Found_RTG = .FALSE. ENDIF ENDIF IF( Found_RTG )THEN Unit = Get_New_File_Unit() OPEN(Unit,FILE=climate_filename,STATUS='OLD',ACCESS='STREAM',& FORM='UNFORMATTED') READ(Unit,IOSTAT=IOS)rtg_sst IF( 0 /= IOS )THEN CALL Gbcs_Warning(.TRUE.,& 'Cannot read in RTG file - read error - using climatology',& 'Load_SST_Climatology','NOAA/NOAA_LoadAuxData') Found_RTG = .FALSE. ENDIF CLOSE(Unit) ! Map to required structure IF( Found_RTG )THEN SST_Clim%Name = 'NOAA RTG data file' SST_Clim%DataFrmt = 'L2Praw' SST_Clim%DataDir = TRIM(L2P_DIR) SST_Clim%DataFile = TRIM(rtg_filename) ! Copy data to array IF( ALLOCATED(SST_Clim%SST) )THEN DEALLOCATE(SST_Clim%SST) ENDIF ALLOCATE(SST_Clim%SST(nx_rtg,ny_rtg,1),STAT=STAT) IF( 0 /= STAT )THEN CALL Gbcs_Critical(.TRUE.,& 'Error allocating memory (SST_Clim%SST)',& 'Load_SST_Climatology','NOAA/NOAA_LoadAuxData') ENDIF SST_Clim%sst(:,:,1) = 0. SST_Clim%SST(:,:,1) = rtg_sst ! Units SST_Clim%Units = Gbcs_Units_Kelvin ! Definition of grid SST_Clim%Grid%Earth_Coord = .TRUE. SST_Clim%Grid%RTM_Format = .FALSE. SST_Clim%Grid%Full_Globe = .TRUE. SST_Clim%Grid%NElems = nx_rtg SST_Clim%Grid%NLines = ny_rtg SST_Clim%Grid%NLevels = 1 SST_Clim%Grid%Upwards = .FALSE. SST_Clim%Grid%LLCrnr%lat = 89.95833 SST_Clim%Grid%LLCrnr%lon = 0.04166 SST_Clim%Grid%URCrnr%lat = -89.95833 SST_Clim%Grid%URCrnr%lon = 359.95833 SST_Clim%Grid%Incr%lat = -0.08333 SST_Clim%Grid%Incr%lon = 0.08333 SST_Clim%Grid%I_LLCrnr%elem = 1 SST_Clim%Grid%I_LLCrnr%line = 1 SST_Clim%Grid%I_URCrnr%elem = nx_rtg SST_Clim%Grid%I_URCrnr%line = ny_rtg SST_Clim%Grid%DataLoc = 1 write(Gbcs_Log_File_Unit,& '('' USING High Resolution SST'')') write(Gbcs_Log_File_Unit,& '('' Read in most recent RTG climatology'')') ENDIF ENDIF ELSE CALL Gbcs_Critical(.TRUE.,& 'No other analysis code included apart from RTG',& 'Load_SST_Climatology','NOAA/NOAA_LoadAuxData') ENDIF ENDIF IF( .not.Found_RTG )THEN ! Find and read in relevant climatology - 5 day version from Pathfinder I = (Time%DayNum/5+1) IF( 73 .lt. I )I=73 write(rtg_filename,'(''clim'',i2.2,''e.sst'')')I climate_filename = TRIM(AUX%GbcsDataPath)//'/'//TRIM(rtg_filename) INQUIRE(File=climate_filename,EXIST=exists) IF( .not.exists )THEN CALL Gbcs_Warning(.TRUE.,& 'Climatology SST file does not exist - not using high res. SST',& 'Load_SST_Climatology','NOAA/NOAA_LoadAuxData') RETURN ENDIF write(Gbcs_Log_File_Unit,& '('' USING High Resolution SST'')') write(Gbcs_Log_File_Unit,& '('' Reading in climatology file : '',a)')& TRIM(climate_filename) ! Read in file Unit = Get_New_File_Unit() OPEN(Unit=Unit,FILE=climate_filename,STATUS='OLD',ACCESS='STREAM',& FORM='UNFORMATTED') READ(Unit,IOSTAT=IOS)climate_sst IF( 0 .ne. IOS )THEN CALL Gbcs_Warning(.TRUE.,& 'Failure in reading Climatology SST file - not using high res. SST',& 'Load_SST_Climatology','NOAA/NOAA_LoadAuxData') CLOSE(Unit) RETURN ENDIF CLOSE(Unit) ! Copy data structure SST_Clim%Name = 'NODC Pathfinder climatology (Casey & Cornillon 1999)' SST_Clim%DataFrmt = 'Binary' SST_Clim%DataDir = AUX%GbcsDataPath SST_Clim%DataFile = TRIM(rtg_filename) ! Copy data to array IF( ALLOCATED(SST_Clim%SST) )THEN DEALLOCATE(SST_Clim%SST) ENDIF write(6,*) 'DEBUG------------------------------------------' write(6,*) ' System status before allocat NoaaClimatology' call system(' free -m ') write(6,*) 'DEBUG------------------------------------------' ALLOCATE(SST_Clim%SST(nx_climate,ny_climate,1),& climate_data(nx_climate,ny_climate),STAT=STAT) IF( 0 /= STAT )THEN write(6,*) 'DEBUG------------------------------------------' write(6,*) ' System status after FAILing ',STAT call system(' free -m ') write(6,*) 'DEBUG------------------------------------------' CALL Gbcs_Critical(.TRUE.,& 'Error allocating memory climatology file (SST_Clim%SST)',& 'Load_SST_Climatology','NOAA/NOAA_LoadAuxData') ENDIF write(6,*) 'DEBUG------------------------------------------' write(6,*) ' System status after SUCCESS',STAT call system(' free -m ') write(6,*) 'DEBUG------------------------------------------' SST_Clim%SST(:,:,1) = 0. ! Note that this dataset is stored E+ve from 180 degrees so force ! switch to E+ve from 0 degrees ! Note also that the data is assumed unsigned int*1 so make sure ! we get the right number.... DO I=1,ny_climate DO J=1,nx_climate JJ = J+nx_climate/2 IF( JJ .gt. nx_climate )THEN JJ = J-nx_climate/2 ENDIF climate_sst_value = climate_sst(J,I) IF( 0 .gt. climate_sst_value )THEN climate_sst_value = climate_sst_value + 256 ENDIF IF( 4 .lt. climate_sst_value )THEN climate_data(JJ,I) = climate_sst_value*0.15+270.15 ELSE climate_data(JJ,I) = 0. ENDIF END DO END DO ! Units SST_Clim%Units = Gbcs_Units_Kelvin ! Definition of grid SST_Clim%Grid%Earth_Coord = .TRUE. SST_Clim%Grid%RTM_Format = .FALSE. SST_Clim%Grid%Full_Globe = .TRUE. SST_Clim%Grid%NElems = nx_climate SST_Clim%Grid%NLines = ny_climate SST_Clim%Grid%NLevels = 1 SST_Clim%Grid%Upwards = .FALSE. SST_Clim%Grid%LLCrnr%lat = 90.0 SST_Clim%Grid%LLCrnr%lon = 0.0 SST_Clim%Grid%URCrnr%lat = -89.9117 SST_Clim%Grid%URCrnr%lon = 359.912 SST_Clim%Grid%Incr%lat = -0.0883789 SST_Clim%Grid%Incr%lon = 0.0878906 SST_Clim%Grid%I_LLCrnr%elem = 1 SST_Clim%Grid%I_LLCrnr%line = 1 SST_Clim%Grid%I_URCrnr%elem = nx_climate SST_Clim%Grid%I_URCrnr%line = ny_climate SST_Clim%Grid%DataLoc = 0 ! ! Re-normalize climatology to MODEL surface temperature ! CALL Renormalize_HighRes_SST( AUX, climate_data, SST_Clim ) SST_Clim%SST(:,:,1) = climate_data DEALLOCATE(climate_data) ! write(Gbcs_Log_File_Unit,& ! '('' Finished renormalizing climatology to model'')') ENDIF END SUBROUTINE Load_SST_Climatology ! Renormalise climatological SST (5 day average) to model data SUBROUTINE Renormalize_HighRes_SST( AUX, climate_data, SST_Clim ) USE GbcsInterpolators TYPE(AUX_Data), INTENT(IN) :: AUX REAL(GbcsReal), INTENT(IN) :: climate_data(:,:) TYPE(High_Res_SST_Struct), INTENT(INOUT) :: SST_Clim INTEGER :: I,J INTEGER :: index REAL(GbcsReal) :: model_sst INTEGER :: nx,ny REAL(GbcsReal) :: av_sst INTEGER :: nxBox, nyBox INTEGER :: STAT REAL(GbcsReal), ALLOCATABLE :: store(:) INTEGER(GbcsInt4), ALLOCATABLE :: nstore(:) LOGICAL :: First TYPE(ImagePixel) :: PX nx = SST_Clim%Grid%NElems ny = SST_Clim%Grid%NLines ! Get 1/2 size of box for averaging climate data nxBox = NINT((AUX%Model%Grid%Incr%lon/SST_Clim%Grid%Incr%lon)/2.) nyBox = NINT((AUX%Model%Grid%Incr%lat/SST_Clim%Grid%Incr%lat)/2.) ! Allocate storage arrays ALLOCATE(store(2*nyBox+1),nstore(2*nyBox+1),STAT=STAT) IF( 0 /= STAT )THEN CALL Gbcs_Critical(.TRUE.,'Cannot allocate storage arrays',& 'Renormalized_HighRes_SST','NOAA/NOAA_LoadAuxData.f90') ENDIF store = 0. nstore = 0 ! Loop round image DO I=1,SST_Clim%Grid%NLines PX%EarthCoord%lat = SST_Clim%Grid%LLCrnr%lat + & (I-1)*SST_Clim%Grid%Incr%lat First = .TRUE. DO J=1,SST_Clim%Grid%NElems IF( 270.9 .le. climate_data(J,I) )THEN PX%EarthCoord%lon = SST_Clim%Grid%LLCrnr%lon + & (J-1)*SST_Clim%Grid%Incr%lon PX%Surface_Type = Gbcs_Surface_Sea ! Get interpolated value from model IF ( AUX%Model%Grid%Earth_Coord ) THEN CALL ModelIndices(PX,AUX%Model%Grid,AUX%Model%Surface_Type) ELSE CALL Gbcs_Critical(.TRUE.,& 'Can only deal with Earth Coord defined models',& 'Renormalize_HighRes_SST','NOAA/NOAA_LoadAuxData.f90') END IF index = Get_SurfTemp_Index(Aux%Model, PX%Surface_Type) IF (index /= -1 ) THEN model_sst = InterpolateField(Aux%Model%Srf(index)%d, PX%LocMP) ELSE CALL Gbcs_Critical(.TRUE.,'Cannot get model SST',& 'Renormalize_HighRes_SST','NOAA/NOAA_LoadAuxData.f90') ENDIF ! Now have model answer - get averaged climate answer av_sst = Get_Average_ClimSST_Slow( J, I, nx, ny, climate_data, & nxBox, nyBox ) ! av_sst = Get_Average_ClimSST( J, I, nx, ny, climate_data, & ! nxBox, nyBox, store, nstore, First ) ! Scale pixel at J,I by ratio ! If no climatology available use simple intepolation IF( 0 .lt. av_sst )THEN SST_Clim%SST(J,I,1) = (model_sst/av_sst)*climate_data(J,I) ELSE SST_Clim%SST(J,I,1) = model_sst ENDIF First = .FALSE. ENDIF END DO END DO ! Deallocate memory DEALLOCATE(store,nstore) END SUBROUTINE Renormalize_HighRes_SST ! Brute force method of getting average REAL(GbcsReal) FUNCTION Get_Average_ClimSST_Slow( I, J, nx, ny, SST, nxBox, & nyBox )RESULT(avSST) INTEGER, INTENT(IN) :: I INTEGER, INTENT(IN) :: J INTEGER, INTENT(IN) :: nx INTEGER, INTENT(IN) :: ny REAL(GbcsReal), INTENT(IN) :: SST(nx,ny) INTEGER, INTENT(IN) :: nxBox INTEGER, INTENT(IN) :: nyBox INTEGER :: II, JJ INTEGER :: Imin,Imax INTEGER :: Jmin,Jmax REAL(GbcsReal) :: total INTEGER :: ntotal IF( 0 .ge. I-nxBox )THEN Imin=1 Imax=I+nxBox ELSE IF( nx .lt. I+nxBox )THEN Imin=I-nxBox Imax=nx ELSE Imin = I-nxBox Imax = I+nxBox ENDIF IF( 0 .ge. J-nyBox )THEN Jmin=1 Jmax=J+nyBox ELSE IF( ny .lt. J+nyBox )THEN Jmin=J-nyBox Jmax=ny ELSE Jmin = J-nyBox Jmax = J+nyBox ENDIF total = 0. ntotal = 0 DO JJ=Jmin,Jmax DO II=Imin,Imax IF( 270.9 .le. SST(II,JJ) )THEN total = total + SST(II,JJ) ntotal = ntotal + 1 ENDIF END DO END DO IF( 0 .lt. ntotal )THEN avSST = total/ntotal ELSE ! If no data the return -1 which means use the interpolated value avSST = -1 ENDIF END FUNCTION Get_Average_ClimSST_Slow ! Meant to be smarter and faster way of getting average ! doesn't work yet, but working on it... REAL(GbcsReal) FUNCTION Get_Average_ClimSST( I, J, nx, ny, SST, nxBox, & nyBox, store, nstore, First )RESULT(avSST) INTEGER, INTENT(IN) :: I INTEGER, INTENT(IN) :: J INTEGER, INTENT(IN) :: nx INTEGER, INTENT(IN) :: ny REAL(GbcsReal), INTENT(IN) :: SST(nx,ny) INTEGER, INTENT(IN) :: nxBox INTEGER, INTENT(IN) :: nyBox REAL(GbcsReal), INTENT(INOUT) :: store(:) INTEGER(GbcsInt4), INTENT(INOUT) :: nstore(:) LOGICAL, INTENT(IN) :: First INTEGER :: II, JJ INTEGER :: Imin, Imax INTEGER :: Jmin, Jmax REAL(GbcsReal) :: total INTEGER :: ntotal INTEGER :: xpos,ypos IF( 0 .ge. I-nxBox )THEN Imin=1 Imax=I+nxBox ELSE IF( nx .lt. I+nxBox )THEN Imin=I-nxBox Imax=nx ELSE Imin = I-nxBox Imax = I+nxBox ENDIF IF( 0 .ge. J-nyBox )THEN Jmin=1 Jmax=J+nyBox ELSE IF( ny .lt. J+nyBox )THEN Jmin=J-nyBox Jmax=ny ELSE Jmin = J-nyBox Jmax = J+nyBox ENDIF ! print *,'===========================' ! print *,I,J ! print *,Imin,Imax,Jmin,Jmax ! print *,nxBox,nyBox ! print *,First ! print *,store,nstore total = 0. ntotal = 0 IF( First )THEN store = 0. nstore = 0 DO JJ=Jmin,Jmax ypos=(JJ-Jmin)+1 DO II=Imin,Imax IF( 270.9 .le. SST(II,JJ) )THEN store(yPos) = store(yPos) + SST(II,JJ) nstore(yPos) = nstore(ypos) + 1 ENDIF END DO total = total + store(ypos) ntotal = ntotal + nstore(ypos) END DO ELSE ! Only allow SST values >= 270.9K DO JJ=Jmin,Jmax ypos=(JJ-Jmin)+1 IF( 0 .ge. I-nxBox )THEN ! In LH Border IF( 270.9 .le. SST(Imax,JJ) )THEN store(yPos) = store(ypos) + SST(Imax,JJ) nstore(yPos) = nstore(yPos) + 1 ENDIF ELSE IF( nx .lt. I+nxBox )THEN !In RH Border IF( 270.9 .le. SST(Imin-1,JJ) )THEN store(yPos) = store(yPos) - SST(Imin-1,JJ) nstore(yPos) = nstore(yPos) - 1 if( 0 .gt. nstore(yPos) )then print *,'1:',Imin-1,JJ print *,SST(Imin-1,JJ),store(yPos),nstore(yPos) stop endif ENDIF ELSE IF( 1 .le. Imin-1 )THEN IF( 270.9 .le. SST(Imin-1,JJ) )THEN store(yPos) = store(yPos) - SST(Imin-1,JJ) nstore(yPos) = nstore(yPos) - 1 if( 0 .gt. nstore(yPos) )then print *,'2:',Imin-1,JJ print *,SST(Imin-1,JJ),store(yPos),nstore(yPos) stop endif ENDIF ENDIF IF( 270.9 .le. SST(Imax,JJ) )THEN store(yPos) = store(yPos) + SST(Imax,JJ) nstore(yPos) = nstore(yPos) + 1 ENDIF ENDIF total = total + store(ypos) ntotal = ntotal + nstore(ypos) END DO ENDIF IF( 0 .lt. ntotal )THEN avSST = total/ntotal ELSE call Gbcs_Critical(.TRUE.,& 'Getting average SST has no available entries',& 'Get_Average_ClimSST','NOAA/NOAA_LoadAuxData.f90') ENDIF END FUNCTION Get_Average_ClimSST ! Shift back time one day SUBROUTINE Time_Back_One_Day( inTime, outTime ) TYPE(DateTime), INTENT(IN) :: inTime TYPE(DateTime), INTENT(OUT) :: outTime INTEGER :: dayNumber ! Subtract 1 day from daynumber dayNumber = inTime%DayNum-1 IF( 0 .ge. dayNumber )THEN CALL date_to_day_number_mcidas(inTime%Year-1,12,31,dayNumber) outTime%Year = inTime%Year-1 outTime%Month = 12 outTime%Day = 31 outTime%DayNum = dayNumber outTime%Hour = inTime%Hour outTime%Minute = inTime%Minute outTime%Seconds = inTime%Seconds outTime%Sec1000 = inTime%Sec1000 ELSE CALL day_number_to_date_mcidas(inTime%Year,dayNumber,outTime%Month,& outTime%Day) outTime%Year = inTime%Year outTime%DayNum = dayNumber outTime%Hour = inTime%Hour outTime%Minute = inTime%Minute outTime%Seconds = inTime%Seconds outTime%Sec1000 = inTime%Sec1000 ENDIF END SUBROUTINE Time_Back_One_Day !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Routines to read in model files such as GFS (NCEP) ! This reads in the MCIDAS GRID file format SUBROUTINE MCIDAS_Load_Aux_Data( SAT , IMG , AUX , RTM, Level26 ) ! ! Description: ! ! Method: ! ! Owner: Manager of TRICS Project ! ! History: ! Version Date Comment ! ------- ---- ------- ! 0.0 05/09/2006 Creation JPDM ! ! ! Code Description: ! Language: Fortran 90 ! ! Based on code from : ! Institute of Atmospheric and Environmental Sciences ! The University of Edinburgh, The Crew Building, The King's Buildings ! West Mains Road, Edinburgh, UK EH9 3JN ! ! Written by Jonathan Mittaz ! NOAA/CICS UMD ! IMPLICIT NONE TYPE(Satellite), INTENT(INOUT) :: SAT TYPE(Imagery), INTENT(INOUT) :: IMG TYPE(Aux_Data), INTENT(INOUT) :: AUX TYPE(RTM_Interface), INTENT(INOUT) :: RTM LOGICAL, OPTIONAL, INTENT(IN) :: Level26 INTEGER :: ii, jj, kk INTEGER :: iline , ielem , ilevel INTEGER :: I_Map INTEGER :: I_Chan INTEGER :: NCoeff, NFuncs, NLen CHARACTER(LEN=256) :: FNAME, CBUFF INTEGER :: STAT INTEGER :: RTM_Num , RTM_ID INTEGER :: NElems,NLines,nlevels,findx INTEGER :: NElems_in, NLines_in INTEGER :: FieldID INTEGER(GbcsInt1), ALLOCATABLE, DIMENSION(:,:) :: surface_type_array ! equivalence(grid_header,cgrid_header) integer, ALLOCATABLE, DIMENSION(:,:) :: integer_in ! From Mcidas include file integer, parameter :: maxgridpt = 2000000 ! IGGET status integer :: status ! Top of atmisphere numbers integer :: iggmax integer :: IGGET character(LEN=4) :: cvar character(LEN=4) :: clu character(LEN=4) :: clev integer :: ngrids integer :: start logical :: file_exists integer :: Grid_File_Number integer :: ios logical :: use_level26 logical :: Found_Z1000 = .FALSE. real :: grid_scale real :: grid_level ! Pressure levels integer, parameter :: nwantedLevels_18 = 18 INTEGER, parameter, dimension(nwantedLevels_18) :: WantedLevels_18 = & (/150,200,250,300,350,400,450,500,550,600,650,& & 700,750,800,850,900,950,1000/) REAL(KIND=GbcsReal), parameter, dimension(19) :: Plevels_18 = & (/ 0.005, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0,& & 450.0, 500.0, 550.0, 600.0, 650.0, 700.0, 750.0, 800.0, 850.0,& & 900.0, 950.0, 1000.0 /) integer, parameter :: nwantedLevels_26 = 26 INTEGER, parameter, dimension(nwantedLevels_26) :: WantedLevels_26 = & (/10,20,30,50,70,100,150,200,250,300,350,400,450,500,550,600,650,& & 700,750,800,850,900,925,950,975,1000/) REAL(KIND=GbcsReal), parameter, dimension(27) :: Plevels_26 = & (/0.005, 10., 20. ,30. ,50. ,70. ,100. ,150. ,200. ,250. ,& & 300. ,350. ,400. ,450. ,500. ,550. ,600. ,650. ,& & 700. ,750. ,800. ,850. ,900. ,925. ,950. ,975. ,1000. /) ! Ozone mixing ratio (US Std Atm) REAL(KIND=GbcsReal), parameter, dimension(19) :: Olevels = & (/ 0.57, 0.4622, 0.2923, 0.1643, 0.09905, 0.06395, 0.05204, & & 0.04438,0.03972, 0.03720, 0.03470, 0.03360, 0.03319, 0.03276, & & 0.03222, 0.03075, 0.02928, 0.02810, 0.02691 /) integer :: nwantedLevels integer, allocatable :: wantedLevels(:) integer, allocatable :: foundLevels(:) integer :: scaleOffset logical :: found_950 use_level26 = .FALSE. IF( PRESENT(Level26) )THEN IF( Level26 )THEN use_level26 = .TRUE. ENDIF ENDIF IF( use_level26 )THEN nwantedLevels = nwantedLevels_26 ALLOCATE(wantedLevels(nwantedLevels),foundLevels(nwantedLevels),& STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical(.TRUE.,'Cannot allocate wantedLevels',& 'MCIDAS_Load_Aux_Data','NOAA_LoadAuxData.f90') ENDIF wantedLevels = wantedLevels_26 foundLevels = 0 ELSE nwantedLevels = nwantedLevels_18 ALLOCATE(wantedLevels(nwantedLevels),foundLevels(nwantedLevels),& STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical(.TRUE.,'Cannot allocate wantedLevels',& 'MCIDAS_Load_Aux_Data','NOAA_LoadAuxData.f90') ENDIF wantedLevels = wantedLevels_18 foundLevels = 0 ENDIF nlevels = AUX%Model%Grid%NLevels nelems = AUX%Model%Grid%NElems nlines = AUX%Model%Grid%NLines IF( nlevels .le. 1 .or. nelems .le. 1 .or. nlines .le. 1 )THEN WRITE(*,*)nlevels,nelems,nlines CALL Gbcs_Critical(.TRUE.,& 'Model (AUX) Nlevels/Nelems/Nlines not sensible',& 'MCIDAS_Load_Aux_Data','NOAA_LoadAuxData.f90') ENDIF IF( nlevels .ne. nwantedLevels )THEN WRITE(*,*)nlevels,nwantedLevels CALL Gbcs_Critical(.TRUE.,& 'Model (AUX) Nlevels mismatch between AUX$Model%Grid and wanterLevels',& 'MCIDAS_Load_Aux_Data','NOAA_LoadAuxData.f90') ENDIF nelems_in = nelems nlines_in = nlines ! Get SZA and landmask - commong to HDF and MCIDAS routines if( IMG%DataFrmt .eq. "MCIDAS" )then CALL MCIDAS_Get_Ancilliary_Data( AUX, IMG ) else CALL NotMCIDAS_Get_Ancilliary_Data( AUX, IMG, SAT ) endif ! Forecast Model 3D Fields !---------------------------- ! Get GRID file number from AUX Data file name read(AUX%Model%DataFile(5:8),*,IOSTAT=IOS)Grid_File_Number IF( 0 .NE. IOS )THEN WRITE(errorString,'('' ERROR: Cannot parse GRID number from '',a)')& TRIM(AUX%Model%DataFile) CALL Gbcs_Critical( .TRUE., & errorString, & 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Now load WRITE(*,'('' Loading AUX data from GRID'',i4.4)')Grid_File_Number IF ( Gbcs_Verbosity > 2 ) & WRITE(Gbcs_Log_File_Unit,*)' Load Forecast Model 3-D Fields' AUX%Model%N_3D_Model_Fields = 0 FieldID = 0 ! Allocate integer output from IGGET ALLOCATE( integer_in(NLines,NElems), STAT=STAT ) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating integer_in', & 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Pressure findx = Gbcs_Map_P_3D IF (ALLOCATED(AUX%Model%Atm(Gbcs_Map_P_3D)%d)) THEN DEALLOCATE(AUX%Model%Atm(Gbcs_Map_P_3D)%d, STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Deallocating AUX%Model%Atm', & 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF ALLOCATE(AUX%Model%Atm(Gbcs_Map_P_3D)%d(nlevels+1,NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating P 3D', & 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Load pressure level data from one of the 3-D fields. IF( nwantedLevels .eq. nwantedLevels_18 )THEN DO jj = 1,nlines DO ii = 1,nelems AUX%Model%Atm(findx)%d(:,ii,jj) = PLevels_18 END DO END DO ELSE DO jj = 1,nlines DO ii = 1,nelems AUX%Model%Atm(findx)%d(:,ii,jj) = PLevels_26 END DO END DO ENDIF ! FieldID = FieldID + 1 ! AUX%Model%N_3D_Model_Fields = AUX%Model%N_3D_Model_Fields + 1 ! AUX%Model%Field_3D_ID( AUX%Model%N_3D_Model_Fields ) = FieldID ! AUX%Model%Gbcs_Map_3D(AUX%Model%N_3D_Fields) = Gbcs_Map_P_3D ! AUX%Model%Field_3D_Units( AUX%Model%N_3D_Model_Fields ) = Gbcs_Units_hPa STAT = Register_3D_Field(AUX%Model, Gbcs_Map_P_3D, Gbcs_Units_hPa) ! Temperature ngrids = iggmax(Grid_File_Number,start) IF( 0 .eq. ngrids )THEN WRITE(*,'('' Grid file = '',a)')TRIM(AUX%Model%DataFile) CALL Gbcs_Critical(.TRUE.,'Cannot read GRID%%%% file - no grids found',& 'MCIDAS_Load_Aux_Data','NOAA/NOAA_LoadAuxData.f90') ENDIF ! Make sure we have a 950 mb level using scaling ScaleOffset = 0 found_950 = .FALSE. GetScale0: DO ii = 1,ngrids status = IGGET(Grid_File_Number,ii,maxgridpt,integer_in,& NLines_in,NElems_in,grid_header) if( 0 .ne. status )then exit GetScale0 endif grid_level = grid_header(10)*(10**(grid_header(11)-ScaleOffset)) IF( 950. .eq. grid_level )THEN found_950 = .TRUE. ENDIF END DO GetScale0 IF( .not. found_950 )THEN ScaleOffset = 1 GetScale1: DO ii = 1,ngrids status = IGGET(Grid_File_Number,ii,maxgridpt,integer_in,& NLines_in,NElems_in,grid_header) if( 0 .ne. status )then exit GetScale1 endif grid_level = grid_header(10)*(10**(grid_header(11)-ScaleOffset)) IF( 950. .eq. grid_level )THEN found_950 = .TRUE. ENDIF END DO GetScale1 ENDIF IF( .not. found_950 )THEN CALL Gbcs_Critical(.TRUE.,'Cannot determine scale for 950 MB level','MCIDAS_Load_Aux_Data',& 'NOAA_LoadAuxData.f90') ENDIF found_Z1000 = .FALSE. ReadLoop: DO ii = 1,ngrids ! read in grids status = IGGET(Grid_File_Number,ii,maxgridpt,integer_in,& NLines_in,NElems_in,grid_header) if( 0 .ne. status )then nlines_in = nlines nelems_in = nelems exit ReadLoop endif grid_scale = 10**grid_header(8) write(6,*) "Grid header(8)", grid_header(8) write(6,*) "Grid Scale", grid_scale cgrid_header = TRANSFER(grid_header,cgrid_header) cvar = cgrid_header(7) grid_level = grid_header(10)*(10**(grid_header(11)-ScaleOffset)) ! grid_level = grid_header(10) do jj=1,nwantedLevels if( 'T ' .eq. cvar(1:2) .and. & wantedLevels(jj) .eq. grid_level )then findx = Gbcs_Map_T_3D IF (.NOT.ALLOCATED(AUX%Model%Atm(findx)%d)) THEN ALLOCATE(AUX%Model%Atm(findx)%d(nlevels+1,NElems,NLines), & STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating T 3D', & 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF AUX%Model%Atm(findx)%d = 0. END IF AUX%Model%Atm(findx)%d(jj+1,:,:) = & TRANSPOSE( REAL(integer_in) ) /grid_scale if( 1 .eq. jj )then AUX%Model%Atm(findx)%d(1,:,:) = TOA_TEMPERATURE ! FieldID = FieldID + 1 ! AUX%Model%N_3D_Model_Fields = AUX%Model%N_3D_Model_Fields + 1 ! AUX%Model%Field_3D_ID( AUX%Model%N_3D_Model_Fields ) = FieldID ! AUX%Model%Gbcs_Map_3D(AUX%Model%N_3D_Fields) = Gbcs_Map_T_3D ! AUX%Model%Field_3D_Units( AUX%Model%N_3D_Model_Fields ) = & ! Gbcs_Units_Kelvin STAT = Register_3D_Field(AUX%Model, Gbcs_Map_T_3D, Gbcs_Units_Kelvin) endif foundlevels(jj) = 1 else if( 'RH' .eq. cvar(1:2) .and. & wantedLevels(jj) .eq. grid_level )then ! Water Vapour findx = Gbcs_Map_WV_3D IF (.NOT.ALLOCATED(AUX%Model%Atm(findx)%d)) THEN ALLOCATE(AUX%Model%Atm(findx)%d(nlevels+1,NElems,NLines), & STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating WV 3D', & 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF AUX%Model%Atm(findx)%d = 0. END IF AUX%Model%Atm(findx)%d(jj+1,:,:) = & TRANSPOSE( REAL(integer_in) ) /grid_scale if( 1 .eq. jj )then AUX%Model%Atm(findx)%d(1,:,:) = TOA_RH ! FieldID = FieldID + 1 ! AUX%Model%N_3D_Model_Fields = AUX%Model%N_3D_Model_Fields + 1 ! AUX%Model%Field_3D_ID( AUX%Model%N_3D_Model_Fields ) = FieldID ! AUX%Model%Gbcs_Map_3D(AUX%Model%N_3D_Fields) = Gbcs_Map_WV_3D ! AUX%Model%Field_3D_Units( AUX%Model%N_3D_Model_Fields ) = & ! Gbcs_Units_RH STAT = Register_3D_Field(AUX%Model, Gbcs_Map_WV_3D, Gbcs_Units_RH) endif else if('RH' .eq. cvar(1:2) .and. 1000. .eq. grid_level )then ! 1000mb height findx = Gbcs_Map_z1000 IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating z1000', & 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF AUX%Model%Srf(findx)%d = TRANSPOSE( REAL(integer_in) ) /grid_scale STAT = Register_2D_Field(AUX%Model, Gbcs_Map_z1000, Gbcs_Units_m) Found_Z1000 = .TRUE. endif end do end do ReadLoop ! Check to see if we've found all the levels IF( nwantedLevels .ne. SUM(foundLevels) )THEN CALL Gbcs_Critical(.TRUE.,'Missing levels in McIDAS GRID file',& 'MCIDAS_Load_Aux_Data','NOAA/NOAA_LoadAuxData.f90') ENDIF IF( .not. Found_Z1000 )THEN CALL Gbcs_Critical(.TRUE.,'Cannot find Z (1000MB) in GRID file',& 'MCIDAS_Load_Aux_Data','NOAA/NOAA_LoadAuxData.f90') ENDIF WRITE(*,'('' Read in '',i2,'' levels from GRID file'')')nwantedLevels ! findx = Gbcs_Map_T_3D ! print *, AUX%Model%Atm(findx)%d(:,50,50) ! findx = Gbcs_Map_WV_3D ! print *, AUX%Model%Atm(findx)%d(:,50,50) ! Ozone findx = Gbcs_Map_O3_3D IF (.NOT.ALLOCATED(AUX%Model%Atm(findx)%d)) THEN ALLOCATE(AUX%Model%Atm(findx)%d(nlevels+1,NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating O3 3D', & 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF ! Set NLevels + 1 due to TOA values AUX%Model%Grid%NLevels = nlevels+1 ! Need everything above to be setup so now get O3 CALL Get_Ozone_climatology( IMG%Mean_Time%Month, AUX ) ! DO ii = 1,nlevels+1 ! ! Get from above values ! AUX%Model%Atm(findx)%d(ii,:,:) = OLevels(ii) ! END DO ! FieldID = FieldID + 1 ! AUX%Model%N_3D_Model_Fields = AUX%Model%N_3D_Model_Fields + 1 ! AUX%Model%Field_3D_ID( AUX%Model%N_3D_Model_Fields ) = FieldID ! AUX%Model%Gbcs_Map_3D(AUX%Model%N_3D_Fields) = Gbcs_Map_O3_3D ! AUX%Model%Field_3D_Units( AUX%Model%N_3D_Model_Fields ) = Gbcs_Units_ppmv STAT = Register_3D_Field(AUX%Model, Gbcs_Map_O3_3D, Gbcs_Units_ppmv) ! Forecast Model 2D Fields !---------------------------- IF ( Gbcs_Verbosity > 2 ) & WRITE(Gbcs_Log_File_Unit,*)' Load Forecast Model 2-D Fields' AUX%Model%N_2D_Model_Fields = 0 FieldID = 0 ! Orography findx = Gbcs_Map_OROG IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating OROG', & 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF AUX%Model%Srf(findx)%d = 0.0 ! FieldID = FieldID + 1 ! AUX%Model%N_2D_Model_Fields = AUX%Model%N_2D_Model_Fields + 1 ! AUX%Model%Srf_Index(findx) = AUX%Model%N_2D_Model_Fields ! AUX%Model%Field_2D_ID( AUX%Model%N_2D_Model_Fields ) = FieldID ! AUX%Model%Field_2D_Map( AUX%Model%N_2D_Model_Fields ) = Gbcs_Map_OROG ! AUX%Model%Field_2D_Units( AUX%Model%N_2D_Model_Fields ) = Gbcs_Units_m STAT = Register_2D_Field(AUX%Model, Gbcs_Map_OROG, Gbcs_Units_m) ! Surface temperature ! Check difference between SST/Surface etc findx = Gbcs_Map_T_SKIN IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating T Skin', & 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF status = IGGET(Grid_File_Number+1,2,maxgridpt,integer_in,& NLines_in,NElems_in,grid_header) if( 0 .ne. status )then write(errorString,'(''IGGET (GRID1001): T Skin - Error status returned : '',i5)')status CALL Gbcs_Critical( .TRUE., & errorString,& 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) endif grid_scale = 10**grid_header(8) cgrid_header = TRANSFER(grid_header,cgrid_header) cvar = cgrid_header(7) clu = cgrid_header(12) clev = cgrid_header(10) write(6,*) "Surface Temperature" write(6,*) "cvar=",cvar," clu=",clu," clev=",clev write(6,*) "(1)",cgrid_header(1),grid_header(1) write(6,*) "(2)",cgrid_header(2),grid_header(2) write(6,*) "(3)",cgrid_header(3),grid_header(3) write(6,*) "(4)",cgrid_header(4),grid_header(4) write(6,*) "(5)",cgrid_header(5),grid_header(5) write(6,*) "(6)",cgrid_header(6),grid_header(6) write(6,*) "(7)",cgrid_header(7),grid_header(7) write(6,*) "(8)",cgrid_header(8),grid_header(8) write(6,*) "(9)",cgrid_header(9),grid_header(9) write(6,*) "(10)",cgrid_header(10),grid_header(10) write(6,*) "(11)",cgrid_header(11),grid_header(11) write(6,*) "(12)",cgrid_header(12),grid_header(12) write(6,*) "(13)",cgrid_header(13),grid_header(13) write(6,*) "(14)",cgrid_header(14),grid_header(14) write(6,*) "(15)",cgrid_header(15),grid_header(15) write(6,*) "(16)",cgrid_header(16),grid_header(16) write(6,*) "(17)",cgrid_header(17),grid_header(17) write(6,*) "(18)",cgrid_header(18),grid_header(18) write(6,*) "(19)",cgrid_header(19),grid_header(19) write(6,*) "(20)",cgrid_header(20),grid_header(20) write(6,*) "(21)",cgrid_header(21),grid_header(21) write(6,*) "(22)",cgrid_header(22),grid_header(22) write(6,*) "(23)",cgrid_header(23),grid_header(23) write(6,*) "(24)",cgrid_header(24),grid_header(24) write(6,*) "(25)",cgrid_header(25),grid_header(25) write(6,*) "(26)",cgrid_header(26),grid_header(26) write(6,*) "(27)",cgrid_header(27),grid_header(27) write(6,*) "(28)",cgrid_header(28),grid_header(28) write(6,*) "(29)",cgrid_header(29),grid_header(29) write(6,*) "(30)",cgrid_header(30) write(6,*) "(31)",cgrid_header(31) write(6,*) "(32)",cgrid_header(32) write(6,*) "(33)",cgrid_header(33) write(6,*) "(34)",cgrid_header(34) write(6,*) "(35)",cgrid_header(35) write(6,*) "(36)",cgrid_header(36) write(6,*) "(37)",cgrid_header(37) write(6,*) "(38)",cgrid_header(38) write(6,*) "(39)",cgrid_header(39) write(6,*) "(40)",cgrid_header(40) write(6,*) "(41)",cgrid_header(41) write(6,*) "(42)",cgrid_header(42) write(6,*) "(43)",cgrid_header(43) write(6,*) "(44)",cgrid_header(44) write(6,*) "(45)",cgrid_header(45) write(6,*) "(46)",cgrid_header(46) write(6,*) "(47)",cgrid_header(47) write(6,*) "(48)",cgrid_header(48) write(6,*) "(49)",cgrid_header(49) write(6,*) "(50)",cgrid_header(50) write(6,*) "(51)",cgrid_header(51) write(6,*) "(52)",cgrid_header(52) write(6,*) "(53)",cgrid_header(53) write(6,*) "(54)",cgrid_header(54) write(6,*) "(55)",cgrid_header(55) write(6,*) "(56)",cgrid_header(56) write(6,*) "(57)",cgrid_header(57) write(6,*) "(58)",cgrid_header(58) write(6,*) "(59)",cgrid_header(59) write(6,*) "(61)",cgrid_header(61) write(6,*) "(62)",cgrid_header(62) write(6,*) "(63)",cgrid_header(63) write(6,*) "(64)",cgrid_header(64) ! Necessary change for GRIB2 decoded GFS surface parameters ! SFC is defined as 1001='SFC' not equivalent number for 'SFC' ! if( 'SFC ' .ne. clev .and. 'T ' .ne. cvar )then if( 1001 .ne. grid_header(10) .and. 'T ' .ne. cvar )then CALL Gbcs_Critical( .TRUE., & 'IGGET (GRID1000): Cannot find T (SURF)',& 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) endif AUX%Model%Srf(findx)%d = TRANSPOSE( REAL(integer_in) ) /grid_scale ! FieldID = FieldID + 1 ! AUX%Model%N_2D_Model_Fields = AUX%Model%N_2D_Model_Fields + 1 ! AUX%Model%Srf_Index(findx) = AUX%Model%N_2D_Model_Fields ! AUX%Model%Field_2D_ID( AUX%Model%N_2D_Model_Fields ) = FieldID ! AUX%Model%Field_2D_Map( AUX%Model%N_2D_Model_Fields ) = Gbcs_Map_T_SKIN ! AUX%Model%Field_2D_Units( AUX%Model%N_2D_Model_Fields ) = Gbcs_Units_Kelvin STAT = Register_2D_Field(AUX%Model, Gbcs_Map_T_SKIN, Gbcs_Units_Kelvin) ! Remove anomolous Tskin values from ocean (min allowed value 271.35) DO jj=1,nlines DO ii=1,nelems IF( AUX%Model%Surface_Type(ii,jj) .eq. Gbcs_Surface_Sea .and. & 271.35 .gt. AUX%Model%Srf(findx)%d(ii,jj) )THEN AUX%Model%Srf(findx)%d(ii,jj) = 271.35 ENDIF END DO END DO ! Total Column Water Vapour findx = Gbcs_Map_TCWV IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating TCWV', & 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF status = IGGET(Grid_File_Number+1,1,maxgridpt,integer_in,& NLines_in,NElems_in,grid_header) if( 0 .ne. status )then write(errorString,'(''IGGET (GRID1001): TCWV - Error status returned : '',i5)')status CALL Gbcs_Critical( .TRUE., & errorString,& 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) endif grid_scale = 10**grid_header(8) write(6,*) "Grid 1001 header(8)", grid_header(8) write(6,*) "Grid 1001 Scale", grid_scale cgrid_header = TRANSFER(grid_header,cgrid_header) cvar = cgrid_header(7) clu = cgrid_header(12) clev = cgrid_header(10) write(6,*) "PWAT " write(6,*) "cvar=",cvar," clu=",clu," clev=",clev write(6,*) "(1)",cgrid_header(1),grid_header(1) write(6,*) "(2)",cgrid_header(2),grid_header(2) write(6,*) "(3)",cgrid_header(3),grid_header(3) write(6,*) "(4)",cgrid_header(4),grid_header(4) write(6,*) "(5)",cgrid_header(5),grid_header(5) write(6,*) "(6)",cgrid_header(6),grid_header(6) write(6,*) "(7)",cgrid_header(7),grid_header(7) write(6,*) "(8)",cgrid_header(8),grid_header(8) write(6,*) "(9)",cgrid_header(9),grid_header(9) write(6,*) "(10)",cgrid_header(10),grid_header(10) write(6,*) "(11)",cgrid_header(11),grid_header(11) write(6,*) "(12)",cgrid_header(12),grid_header(12) write(6,*) "(13)",cgrid_header(13),grid_header(13) write(6,*) "(14)",cgrid_header(14),grid_header(14) write(6,*) "(15)",cgrid_header(15),grid_header(15) write(6,*) "(16)",cgrid_header(16),grid_header(16) write(6,*) "(17)",cgrid_header(17),grid_header(17) write(6,*) "(18)",cgrid_header(18),grid_header(18) write(6,*) "(19)",cgrid_header(19),grid_header(19) write(6,*) "(20)",cgrid_header(20),grid_header(20) write(6,*) "(21)",cgrid_header(21),grid_header(21) write(6,*) "(22)",cgrid_header(22),grid_header(22) write(6,*) "(23)",cgrid_header(23),grid_header(23) write(6,*) "(24)",cgrid_header(24),grid_header(24) write(6,*) "(25)",cgrid_header(25),grid_header(25) write(6,*) "(26)",cgrid_header(26),grid_header(26) write(6,*) "(27)",cgrid_header(27),grid_header(27) write(6,*) "(28)",cgrid_header(28),grid_header(28) write(6,*) "(29)",cgrid_header(29),grid_header(29) if( 'EATM' .ne. clev .and. 'PWAT' .ne. cvar )then CALL Gbcs_Critical( .TRUE., & 'IGGET (GRID1000): Cannot find TCWV (SURF)',& 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) endif AUX%Model%Srf(findx)%d = TRANSPOSE( REAL(integer_in) ) /grid_scale ! FieldID = FieldID + 1 ! AUX%Model%N_2D_Model_Fields = AUX%Model%N_2D_Model_Fields + 1 ! AUX%Model%Srf_Index(findx) = AUX%Model%N_2D_Model_Fields ! AUX%Model%Field_2D_ID( AUX%Model%N_2D_Model_Fields ) = FieldID ! AUX%Model%Field_2D_Map( AUX%Model%N_2D_Model_Fields ) = Gbcs_Map_TCWV ! AUX%Model%Field_2D_Units( AUX%Model%N_2D_Model_Fields ) = Gbcs_Units_Kg_per_m2 STAT = Register_2D_Field(AUX%Model, Gbcs_Map_TCWV, Gbcs_Units_Kg_per_m2) ! Wind - U findx = Gbcs_Map_U_SURF IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating U Surf', & 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF status = IGGET(Grid_File_Number+1,4,maxgridpt,integer_in,& NLines_in,NElems_in,grid_header) if( 0 .ne. status )then write(errorString,'(''IGGET (GRID1001): Wind U - Error status returned : '',i5)')status CALL Gbcs_Critical( .TRUE., & errorString,& 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) endif grid_scale = 10**grid_header(8) write(6,*) "Grid 1001 header(8)", grid_header(8) write(6,*) "Grid 1001 Scale", grid_scale cgrid_header = TRANSFER(grid_header,cgrid_header) cvar = cgrid_header(7) clu = cgrid_header(12) clev = cgrid_header(10) write(6,*) "(1)",cgrid_header(1) write(6,*) "(2)",cgrid_header(2) write(6,*) "(3)",cgrid_header(3) write(6,*) "(4)",cgrid_header(4) write(6,*) "(5)",cgrid_header(5) write(6,*) "(6)",cgrid_header(6) write(6,*) "(7)",cgrid_header(7) write(6,*) "(8)",cgrid_header(8) write(6,*) "(9)",cgrid_header(9) write(6,*) "(10)",cgrid_header(10) write(6,*) "(11)",cgrid_header(11) write(6,*) "(12)",cgrid_header(12) write(6,*) "(13)",cgrid_header(13) write(6,*) "(14)",cgrid_header(14) write(6,*) "(15)",cgrid_header(15) write(6,*) "(16)",cgrid_header(16) write(6,*) "(17)",cgrid_header(17) write(6,*) "(18)",cgrid_header(18) write(6,*) "(19)",cgrid_header(19) write(6,*) "(20)",cgrid_header(20) write(6,*) "(21)",cgrid_header(21) write(6,*) "(22)",cgrid_header(22) write(6,*) "(23)",cgrid_header(23) write(6,*) "(24)",cgrid_header(24) write(6,*) "(25)",cgrid_header(25) write(6,*) "(26)",cgrid_header(26) write(6,*) "(27)",cgrid_header(27) write(6,*) "(28)",cgrid_header(28) write(6,*) "(29)",cgrid_header(29) write(6,*) "U " write(6,*) "cvar=",cvar," clu=",clu," clev=",clev if( 10 .ne. grid_header(10) .and. 'U ' .ne. cvar )then CALL Gbcs_Critical( .TRUE., & 'IGGET (GRID1001): Cannot find U (10M)',& 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) endif AUX%Model%Srf(findx)%d = TRANSPOSE( REAL(integer_in) ) /grid_scale ! print *,'U:',minval(integer_in)/100.0,maxval(integer_in)/100.0 ! FieldID = FieldID + 1 ! AUX%Model%N_2D_Model_Fields = AUX%Model%N_2D_Model_Fields + 1 ! AUX%Model%Srf_Index(findx) = AUX%Model%N_2D_Model_Fields ! AUX%Model%Field_2D_ID( AUX%Model%N_2D_Model_Fields ) = FieldID ! AUX%Model%Field_2D_Map( AUX%Model%N_2D_Model_Fields ) = Gbcs_Map_U_SURF ! AUX%Model%Field_2D_Units( AUX%Model%N_2D_Model_Fields ) = Gbcs_Units_m_per_s STAT = Register_2D_Field(AUX%Model, Gbcs_Map_U_SURF, Gbcs_Units_m_per_s) ! Wind - V findx = Gbcs_Map_V_SURF IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating V Surf', & 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF status = IGGET(Grid_File_Number+1,5,maxgridpt,integer_in,& NLines_in,NElems_in,grid_header) if( 0 .ne. status )then write(errorString,'(''IGGET (GRID1001): Wind V - Error status returned : '',i5)')status CALL Gbcs_Critical( .TRUE., & errorString,& 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) endif grid_scale = 10**grid_header(8) write(6,*) "Grid 1001 header(8)", grid_header(8) write(6,*) "Grid 1001 Scale", grid_scale cgrid_header = TRANSFER(grid_header,cgrid_header) cvar = cgrid_header(7) clu = cgrid_header(12) clev = cgrid_header(10) write(6,*) "V " write(6,*) "cvar=",cvar," clu=",clu," clev=",clev write(6,*) "(1)",cgrid_header(1) write(6,*) "(2)",cgrid_header(2) write(6,*) "(3)",cgrid_header(3) write(6,*) "(4)",cgrid_header(4) write(6,*) "(5)",cgrid_header(5) write(6,*) "(6)",cgrid_header(6) write(6,*) "(7)",cgrid_header(7) write(6,*) "(8)",cgrid_header(8) write(6,*) "(9)",cgrid_header(9) write(6,*) "(10)",cgrid_header(10) write(6,*) "(11)",cgrid_header(11) write(6,*) "(12)",cgrid_header(12) write(6,*) "(13)",cgrid_header(13) write(6,*) "(14)",cgrid_header(14) write(6,*) "(15)",cgrid_header(15) write(6,*) "(16)",cgrid_header(16) write(6,*) "(17)",cgrid_header(17) write(6,*) "(18)",cgrid_header(18) write(6,*) "(19)",cgrid_header(19) write(6,*) "(20)",cgrid_header(20) write(6,*) "(21)",cgrid_header(21) write(6,*) "(22)",cgrid_header(22) write(6,*) "(23)",cgrid_header(23) write(6,*) "(24)",cgrid_header(24) write(6,*) "(25)",cgrid_header(25) write(6,*) "(26)",cgrid_header(26) write(6,*) "(27)",cgrid_header(27) write(6,*) "(28)",cgrid_header(28) write(6,*) "(29)",cgrid_header(29) if( 10 .ne. grid_header(10) .and. 'V ' .ne. cvar )then CALL Gbcs_Critical( .TRUE., & 'IGGET (GRID1000): Cannot find V (10M)',& 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) endif AUX%Model%Srf(findx)%d = TRANSPOSE( REAL(integer_in) ) /grid_scale ! FieldID = FieldID + 1 ! AUX%Model%N_2D_Model_Fields = AUX%Model%N_2D_Model_Fields + 1 ! AUX%Model%Srf_Index(findx) = AUX%Model%N_2D_Model_Fields ! AUX%Model%Field_2D_ID( AUX%Model%N_2D_Model_Fields ) = FieldID ! AUX%Model%Field_2D_Map( AUX%Model%N_2D_Model_Fields ) = Gbcs_Map_V_SURF ! AUX%Model%Field_2D_Units( AUX%Model%N_2D_Model_Fields ) = Gbcs_Units_m_per_s STAT = Register_2D_Field(AUX%Model, Gbcs_Map_V_SURF, Gbcs_Units_m_per_s) ! print *,'V:',minval(integer_in)/100.0,maxval(integer_in)/100.0 ! Wind - 10m findx = Gbcs_Map_Wind_10M IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating 10m Wind', & 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF do jj=1,nlines do ii=1,nelems if( 25 .lt. AUX%Model%Srf(Gbcs_Map_U_SURF)%d(ii,jj) .or. & 25 .lt. AUX%Model%Srf(Gbcs_Map_V_SURF)%d(ii,jj) )then if( 25 .lt. AUX%Model%Srf(Gbcs_Map_U_SURF)%d(ii,jj) .and. & 25 .lt. AUX%Model%Srf(Gbcs_Map_V_SURF)%d(ii,jj) )then AUX%Model%Srf(findx)%d(ii,jj) = 8 else if( 25 .lt. AUX%Model%Srf(Gbcs_Map_U_SURF)%d(ii,jj) )then AUX%Model%Srf(findx)%d(ii,jj) = & AUX%Model%Srf(Gbcs_Map_V_SURF)%d(ii,jj) else if( 25 .lt. AUX%Model%Srf(Gbcs_Map_V_SURF)%d(ii,jj) )then AUX%Model%Srf(findx)%d(ii,jj) = & AUX%Model%Srf(Gbcs_Map_U_SURF)%d(ii,jj) endif else AUX%Model%Srf(findx)%d(ii,jj) = sqrt(& AUX%Model%Srf(Gbcs_Map_V_SURF)%d(ii,jj)**2 + & AUX%Model%Srf(Gbcs_Map_U_SURF)%d(ii,jj)**2) endif end do end do ! FieldID = FieldID + 1 ! AUX%Model%N_2D_Model_Fields = AUX%Model%N_2D_Model_Fields + 1 ! AUX%Model%Srf_Index(findx) = AUX%Model%N_2D_Model_Fields ! AUX%Model%Field_2D_ID( AUX%Model%N_2D_Model_Fields ) = FieldID ! AUX%Model%Field_2D_Map( AUX%Model%N_2D_Model_Fields ) = Gbcs_Map_Wind_10M ! AUX%Model%Field_2D_Units( AUX%Model%N_2D_Model_Fields ) = Gbcs_Units_m_per_s STAT = Register_2D_Field(AUX%Model, Gbcs_Map_Wind_10M, Gbcs_Units_m_per_s) ! Get surface pressure findx = Gbcs_Map_P_Surf IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating P Surf', & 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF ! This is taken from old operational code - but surface pressure exists in ! GRID file #ifndef OLDPRESSURE ! Read in from GRID file status = IGGET(Grid_File_Number+1,8,maxgridpt,integer_in,& NLines_in,NElems_in,grid_header) if( 0 .ne. status )then write(errorString,'(''IGGET (GRID1001): Surface Pressure - Error status returned : '',i5)')status CALL Gbcs_Critical( .TRUE., & errorString,& 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) endif grid_scale = 10**grid_header(8) write(6,*) "Grid 1001 header(8)", grid_header(8) write(6,*) "Grid 1001 Scale", grid_scale cgrid_header = TRANSFER(grid_header,cgrid_header) cvar = cgrid_header(7) clu = cgrid_header(12) clev = cgrid_header(10) write(6,*) "Surface Preassure " write(6,*) "cvar=",cvar," clu=",clu," clev=",clev write(6,*) "(1)",cgrid_header(1),grid_header(1) write(6,*) "(2)",cgrid_header(2),grid_header(2) write(6,*) "(3)",cgrid_header(3),grid_header(3) write(6,*) "(4)",cgrid_header(4),grid_header(4) write(6,*) "(5)",cgrid_header(5),grid_header(5) write(6,*) "(6)",cgrid_header(6),grid_header(6) write(6,*) "(7)",cgrid_header(7),grid_header(7) write(6,*) "(8)",cgrid_header(8),grid_header(8) write(6,*) "(9)",cgrid_header(9),grid_header(9) write(6,*) "(10)",cgrid_header(10),grid_header(10) write(6,*) "(11)",cgrid_header(11),grid_header(11) write(6,*) "(12)",cgrid_header(12),grid_header(12) write(6,*) "(13)",cgrid_header(13),grid_header(13) write(6,*) "(14)",cgrid_header(14),grid_header(14) write(6,*) "(15)",cgrid_header(15),grid_header(15) write(6,*) "(16)",cgrid_header(16),grid_header(16) write(6,*) "(17)",cgrid_header(17),grid_header(17) write(6,*) "(18)",cgrid_header(18),grid_header(18) write(6,*) "(19)",cgrid_header(19),grid_header(19) write(6,*) "(20)",cgrid_header(20),grid_header(20) write(6,*) "(21)",cgrid_header(21),grid_header(21) write(6,*) "(22)",cgrid_header(22),grid_header(22) write(6,*) "(23)",cgrid_header(23),grid_header(23) write(6,*) "(24)",cgrid_header(24),grid_header(24) write(6,*) "(25)",cgrid_header(25),grid_header(25) write(6,*) "(26)",cgrid_header(26),grid_header(26) write(6,*) "(27)",cgrid_header(27),grid_header(27) write(6,*) "(28)",cgrid_header(28),grid_header(28) write(6,*) "(29)",cgrid_header(29),grid_header(29) ! Necessary change for GRIB2 decoded GFS surface parameters ! SFC is defined as 1001='SFC' not equivalent number for 'SFC' ! if( 'SFC ' .ne. clev .and. 'Pres' .ne. cvar )then if( 1001 .ne. grid_header(10) .and. 'P ' .ne. cvar )then CALL Gbcs_Critical( .TRUE., & 'IGGET (GRID1000): Cannot find Surface Pressure',& 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) endif AUX%Model%Srf(findx)%d = TRANSPOSE( REAL(integer_in) ) /grid_scale #else ! Calculate surface pressure call surfpressure_GOES(nlines, nelems, & AUX%Model%Srf(Gbcs_Map_z1000)%d, & AUX%Model%Atm(Gbcs_Map_T_3D)%d(nlevels,:,:), & AUX%Model%Atm(Gbcs_Map_P_3D)%d(nlevels,:,:), & AUX%Model%Srf(Gbcs_Map_P_Surf)%d ) #endif ! FieldID = FieldID + 1 ! AUX%Model%N_2D_Model_Fields = AUX%Model%N_2D_Model_Fields + 1 ! AUX%Model%Srf_Index(findx) = AUX%Model%N_2D_Model_Fields ! AUX%Model%Field_2D_ID( AUX%Model%N_2D_Model_Fields ) = FieldID ! AUX%Model%Field_2D_Map( AUX%Model%N_2D_Model_Fields ) = Gbcs_Map_P_Surf ! AUX%Model%Field_2D_Units( AUX%Model%N_2D_Model_Fields ) = Gbcs_Units_hPa STAT = Register_2D_Field(AUX%Model, Gbcs_Map_P_Surf, Gbcs_Units_hPa) ! Total surface RH findx = Gbcs_Map_WV_2M IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating 2n WV', & 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF status = IGGET(Grid_File_Number+1,3,maxgridpt,integer_in,& NLines_in,NElems_in,grid_header) if( 0 .ne. status )then write(errorString,'(''IGGET (GRID1001): 2m WV - Error status returned : '',i5)')status CALL Gbcs_Critical( .TRUE., & errorString,& 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) endif grid_scale = 10**grid_header(8) write(6,*) "Grid 1001 header(8)", grid_header(8) write(6,*) "Grid 1001 Scale", grid_scale cgrid_header = TRANSFER(grid_header,cgrid_header) cvar = cgrid_header(7) clu = cgrid_header(12) clev = cgrid_header(10) write(6,*) "Surface RH " write(6,*) "cvar=",cvar," clu=",clu," clev=",clev write(6,*) "(1)",cgrid_header(1),grid_header(1) write(6,*) "(2)",cgrid_header(2),grid_header(2) write(6,*) "(3)",cgrid_header(3),grid_header(3) write(6,*) "(4)",cgrid_header(4),grid_header(4) write(6,*) "(5)",cgrid_header(5),grid_header(5) write(6,*) "(6)",cgrid_header(6),grid_header(6) write(6,*) "(7)",cgrid_header(7),grid_header(7) write(6,*) "(8)",cgrid_header(8),grid_header(8) write(6,*) "(9)",cgrid_header(9),grid_header(9) write(6,*) "(10)",cgrid_header(10),grid_header(10) write(6,*) "(11)",cgrid_header(11),grid_header(11) write(6,*) "(12)",cgrid_header(12),grid_header(12) write(6,*) "(13)",cgrid_header(13),grid_header(13) write(6,*) "(14)",cgrid_header(14),grid_header(14) write(6,*) "(15)",cgrid_header(15),grid_header(15) write(6,*) "(16)",cgrid_header(16),grid_header(16) write(6,*) "(17)",cgrid_header(17),grid_header(17) write(6,*) "(18)",cgrid_header(18),grid_header(18) write(6,*) "(19)",cgrid_header(19),grid_header(19) write(6,*) "(20)",cgrid_header(20),grid_header(20) write(6,*) "(21)",cgrid_header(21),grid_header(21) write(6,*) "(22)",cgrid_header(22),grid_header(22) write(6,*) "(23)",cgrid_header(23),grid_header(23) write(6,*) "(24)",cgrid_header(24),grid_header(24) write(6,*) "(25)",cgrid_header(25),grid_header(25) write(6,*) "(26)",cgrid_header(26),grid_header(26) write(6,*) "(27)",cgrid_header(27),grid_header(27) write(6,*) "(28)",cgrid_header(28),grid_header(28) write(6,*) "(29)",cgrid_header(29),grid_header(29) ! Necessary change for GRIB2 decoded GFS surface parameters ! SFC is defined as 1001='SFC' not equivalent number for 'SFC' ! if( 'SFC ' .ne. clev .and. 'RH ' .ne. cvar )then if( 1001 .ne. grid_header(10) .and. 'RH ' .ne. cvar )then CALL Gbcs_Critical( .TRUE., & 'IGGET (GRID1001): Cannot find T (SURF)',& 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) endif AUX%Model%Srf(findx)%d = TRANSPOSE( REAL(integer_in) ) /grid_scale ! print *,'RH:',minval(integer_in)/100.0,maxval(integer_in)/100.0 ! FieldID = FieldID + 1 ! AUX%Model%N_2D_Model_Fields = AUX%Model%N_2D_Model_Fields + 1 ! AUX%Model%Srf_Index(findx) = AUX%Model%N_2D_Model_Fields ! AUX%Model%Field_2D_ID( AUX%Model%N_2D_Model_Fields ) = FieldID ! AUX%Model%Field_2D_Map( AUX%Model%N_2D_Model_Fields ) = Gbcs_Map_WV_2M ! AUX%Model%Field_2D_Units( AUX%Model%N_2D_Model_Fields ) = Gbcs_Units_RH STAT = Register_2D_Field(AUX%Model, Gbcs_Map_WV_2M, Gbcs_Units_RH) #ifdef NEWCODE ! 2m temperature findx = Gbcs_Map_T_2M IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating 2n WV', & 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF status = IGGET(Grid_File_Number+1,9,maxgridpt,integer_in,& NLines_in,NElems_in,grid_header) if( 0 .ne. status )then write(errorString,'(''IGGET (GRID1001): 2m T - Error status returned : '',i5)')status CALL Gbcs_Critical( .TRUE., & errorString,& 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) endif grid_scale = 10**grid_header(8) write(6,*) "Grid 1001 header(8)", grid_header(8) write(6,*) "Grid 1001 Scale", grid_scale cgrid_header = TRANSFER(grid_header,cgrid_header) cvar = cgrid_header(7) clu = cgrid_header(12) clev = cgrid_header(10) write(6,*) "2M Temperature " write(6,*) "cvar=",cvar," clu=",clu," clev=",clev write(6,*) "(1)",cgrid_header(1),grid_header(1) write(6,*) "(2)",cgrid_header(2),grid_header(2) write(6,*) "(3)",cgrid_header(3),grid_header(3) write(6,*) "(4)",cgrid_header(4),grid_header(4) write(6,*) "(5)",cgrid_header(5),grid_header(5) write(6,*) "(6)",cgrid_header(6),grid_header(6) write(6,*) "(7)",cgrid_header(7),grid_header(7) write(6,*) "(8)",cgrid_header(8),grid_header(8) write(6,*) "(9)",cgrid_header(9),grid_header(9) write(6,*) "(10)",cgrid_header(10),grid_header(10) write(6,*) "(11)",cgrid_header(11),grid_header(11) write(6,*) "(12)",cgrid_header(12),grid_header(12) write(6,*) "(13)",cgrid_header(13),grid_header(13) write(6,*) "(14)",cgrid_header(14),grid_header(14) write(6,*) "(15)",cgrid_header(15),grid_header(15) write(6,*) "(16)",cgrid_header(16),grid_header(16) write(6,*) "(17)",cgrid_header(17),grid_header(17) write(6,*) "(18)",cgrid_header(18),grid_header(18) write(6,*) "(19)",cgrid_header(19),grid_header(19) write(6,*) "(20)",cgrid_header(20),grid_header(20) write(6,*) "(21)",cgrid_header(21),grid_header(21) write(6,*) "(22)",cgrid_header(22),grid_header(22) write(6,*) "(23)",cgrid_header(23),grid_header(23) write(6,*) "(24)",cgrid_header(24),grid_header(24) write(6,*) "(25)",cgrid_header(25),grid_header(25) write(6,*) "(26)",cgrid_header(26),grid_header(26) write(6,*) "(27)",cgrid_header(27),grid_header(27) write(6,*) "(28)",cgrid_header(28),grid_header(28) write(6,*) "(29)",cgrid_header(29),grid_header(29) if( 2 .ne. clev .and. 'T ' .ne. cvar )then CALL Gbcs_Critical( .TRUE., & 'IGGET (GRID1001): Cannot find T (2m)',& 'MCIDAS_Load_AUX_Data', & 'NOAA/NOAA_LoadAuxData.f90' ) endif AUX%Model%Srf(findx)%d = TRANSPOSE( REAL(integer_in) ) /grid_scale STAT = Register_2D_Field(AUX%Model, Gbcs_Map_T_2M, Gbcs_Units_Kelvin) #endif RETURN DEALLOCATE(wantedLevels,foundlevels) END SUBROUTINE MCIDAS_Load_Aux_Data SUBROUTINE surfpressure_GOES(nlines, nelems, z_lev, t_lev, p_lev, p_surf) ! --------------------------- ! FFM.1.1.1. Surface Pressure ! --------------------------- ! 03/09/03 ! SNM ! Function to calculate the surface pressure, p_surf (mB), ! given the temperature t_lev (K*100) and altitude, z_lev (m*100) ! at pressure p_lev (mB) ! ---------------------- Declare Variables ---------------------------- USE GbcsTypes ! INPUTS integer, intent(in) :: nlines, nelems ! Altitude (m) at pressure level, p_lev real(KIND=GbcsReal), intent(in) :: z_lev(nelems,nlines) ! Temperature (K) at pressure level, p_lev real(KIND=GbcsReal), intent(in) :: t_lev(nelems,nlines) ! Pressure (mB) at bottom level of profile real(KIND=GbcsReal), intent(in) :: p_lev(nelems,nlines) ! OUTPUTS ! Pressure (mB) at surface(z_surf= 0m) real (KIND=GbcsReal), intent(out) :: p_surf(nelems,nlines) ! LOCAL VARIABLES ! Gas constant for dry air (J/K/kg) real, parameter :: R = 287.04 !parameter (R = 286.9968933) !parameter (R = 287.04) ! Acceleration due to gravity (m/s^2) real, parameter :: grav = 9.80616 !parameter (grav = 9.80616) ! Altitude at surface (m) real, parameter :: z_surf = 0.0 real t_lev2 ! Temperature (K) at pressure level, p_lev real z_lev2 ! Altitude (m) at pressure level, p_lev ! ------------------------- Processing -------------------------------- ! Calculate the surface pressure p_surf = p_lev * EXP((-1.0) * (grav / (R * t_lev)) * (z_surf - z_lev)) RETURN END SUBROUTINE surfpressure_GOES #ifdef USE_HDF SUBROUTINE NOAA_HDF_Load_Aux_Data( SAT, IMG, AUX, RTM, OPER_VALS ) USE GbcsAtmPhysics USE GbcsPDFLoaders USE GbcsCovMatLoaders IMPLICIT NONE TYPE(Satellite), INTENT(INOUT) :: SAT TYPE(Imagery), INTENT(INOUT) :: IMG TYPE(Aux_Data), INTENT(INOUT) :: AUX TYPE(RTM_Interface), INTENT(INOUT) :: RTM LOGICAL, INTENT(IN) :: OPER_VALS ! Use operational levels INTEGER :: start_year INTEGER :: start_month INTEGER :: start_day INTEGER :: start_hour_anal INTEGER :: start_hour_forecast INTEGER :: end_year INTEGER :: end_month INTEGER :: end_day INTEGER :: end_hour_anal INTEGER :: end_hour_forecast CHARACTER(LEN=512) :: gfs_hdf_archive CHARACTER(LEN=512) :: file1 CHARACTER(LEN=512) :: file2 REAL, ALLOCATABLE :: landseamask(:,:) INTEGER :: day_num REAL :: hour INTEGER :: I, J ! Get SZA/Land mask for model - must have access to MCIDAS data files if( IMG%DataFrmt .eq. "MCIDAS" )then CALL MCIDAS_Get_Ancilliary_Data( AUX, IMG ) else CALL NotMCIDAS_Get_Ancilliary_Data( AUX, IMG, SAT ) endif ! Get name of HDF file - assumes Wen Meng archive is in GFS_HDF_ARCHIVE ! directory ! Get surrounding model files ! Start start_year = IMG%Mean_Time%Year start_month = IMG%Mean_Time%Month start_day = IMG%Mean_Time%Day start_hour_anal = (IMG%Mean_Time%Hour/6)*6 start_hour_forecast = ((IMG%Mean_Time%Hour - start_hour_anal)/3)*3 ! Next file in line if( 3 .eq. start_hour_forecast )then end_hour_anal = start_hour_anal+6 end_hour_forecast = 0 else end_hour_anal = start_hour_anal end_hour_forecast = 3 endif if( end_hour_anal .lt. 24 )then end_day = start_day end_month = start_month end_year = start_year else end_hour_anal = 0 end_hour_forecast = 0 call date_to_day_number_mcidas(start_year,start_month,start_day,day_num) day_num = day_num+1 if( (0 .eq. MOD(start_year,4) .and. 0 .ne. MOD(start_year,100)) .or. & (0 .eq. MOD(start_year,400)) )then if( day_num .gt. 366 )then end_year = start_year+1 end_month = 1 end_day = 1 else end_year = start_year call day_number_to_date_mcidas( end_year, day_num, end_month, end_day ) endif else if( day_num .gt. 365 )then end_year = start_year+1 end_month = 1 end_day = 1 else end_year = start_year call day_number_to_date_mcidas( end_year, day_num, end_month, end_day ) endif endif endif ! Get floating fraction of 3 hour window hour = IMG%Mean_Time%Hour + IMG%Mean_Time%Minute/60. + & IMG%Mean_Time%Seconds/3600. hour = hour - start_hour_anal - start_hour_forecast ! Now have dates and times of bracketing model files ! make filenames write(file1,& '(a,''/'',i4.4,''/'',i2.2,''/'',i4.4,i2.2,i2.2,''gdas1.t'',i2.2,''z.pgrbf'',i2.2,''.hdf'')')& TRIM(AUX%Model%DataDir),start_year,start_month,start_year,start_month,& start_day,start_hour_anal,start_hour_forecast write(file2,& '(a,''/'',i4.4,''/'',i2.2,''/'',i4.4,i2.2,i2.2,''gdas1.t'',i2.2,''z.pgrbf'',i2.2,''.hdf'')')& TRIM(AUX%Model%DataDir),end_year,end_month,end_year,end_month,& end_day,end_hour_anal,end_hour_forecast write(*,'('' Interpolating between : '')') write(*,'('' '',a)')TRIM(file1) write(*,'('' '',a)')TRIM(file2) IF( OPER_VALS )THEN CALL NOAA_HDF_Get_Aux_OPER( file1, file2, hour, AUX, IMG ) ELSE CALL NOAA_HDF_Get_Aux( file1, file2, hour, AUX, IMG, LandSeaMask ) ! Reset land sea mask given files (sets ice mask) DO I=1,AUX%Model%Grid%Nlines DO J=1,AUX%Model%Grid%NElems if( 1 .eq. LandSeaMask(J,I) )then AUX%Model%Surface_Type(J,I) = Gbcs_Surface_Land else if( 0 .eq. LandSeaMask(J,I) )then AUX%Model%Surface_Type(J,I) = Gbcs_Surface_Sea else if( 2 .eq. LandSeaMask(J,I) )then AUX%Model%Surface_Type(J,I) = Gbcs_Surface_Ice endif END DO END DO DEALLOCATE(LandSeaMask) ENDIF END SUBROUTINE NOAA_HDF_Load_Aux_Data ! Fill AUX arrays from HDF version of NCEP files - full data ! Make it generic SUBROUTINE NOAA_HDF_Get_AuX( file1, file2, hour, AUX, IMG, LandSeaMask ) ! INPUT CHARACTER(LEN=*), INTENT(IN) :: file1 CHARACTER(LEN=*), INTENT(IN) :: file2 REAL, INTENT(IN) :: hour TYPE(Aux_Data), INTENT(INOUT) :: AUX TYPE(Imagery), INTENT(IN) :: IMG REAL, ALLOCATABLE, INTENT(OUT) :: landseamask(:,:) logical :: file_there character(len=256) :: file1in character(len=256) :: file2in character(len=256) :: command INTEGER :: DateTime(8) CHARACTER(LEN=80) :: dateString INTEGER :: nelems1 INTEGER :: nlines1 INTEGER :: nlevels1 INTEGER :: nelems2 INTEGER :: nlines2 INTEGER :: nlevels2 REAL, ALLOCATABLE :: pressure1(:) REAL, ALLOCATABLE :: height1(:,:,:) REAL, ALLOCATABLE :: temperature1(:,:,:) REAL, ALLOCATABLE :: rh1(:,:,:) REAL, ALLOCATABLE :: ozone1(:,:,:) REAL, ALLOCATABLE :: sfc_t1(:,:) REAL, ALLOCATABLE :: sfc_rh1(:,:) REAL, ALLOCATABLE :: sfc_press1(:,:) REAL, ALLOCATABLE :: tcwv1(:,:) REAL, ALLOCATABLE :: wind_u1(:,:) REAL, ALLOCATABLE :: wind_v1(:,:) REAL, ALLOCATABLE :: height_1000mb_1(:,:) REAL, ALLOCATABLE :: two_m_temp1(:,:) REAL, ALLOCATABLE :: pressure2(:) REAL, ALLOCATABLE :: height2(:,:,:) REAL, ALLOCATABLE :: temperature2(:,:,:) REAL, ALLOCATABLE :: rh2(:,:,:) REAL, ALLOCATABLE :: ozone2(:,:,:) REAL, ALLOCATABLE :: sfc_t2(:,:) REAL, ALLOCATABLE :: sfc_rh2(:,:) REAL, ALLOCATABLE :: sfc_press2(:,:) REAL, ALLOCATABLE :: tcwv2(:,:) REAL, ALLOCATABLE :: wind_u2(:,:) REAL, ALLOCATABLE :: wind_v2(:,:) REAL, ALLOCATABLE :: height_1000mb_2(:,:) REAL, ALLOCATABLE :: two_m_temp2(:,:) REAL, ALLOCATABLE :: landseamask2(:,:) LOGICAL, ALLOCATABLE :: ozone_there1(:) LOGICAL, ALLOCATABLE :: ozone_there2(:) REAL :: min_Long1, max_Long1, long_step1 REAL :: min_Lat1, max_Lat1, lat_step1 REAL :: min_Long2, max_Long2, long_step2 REAL :: min_Lat2, max_Lat2, lat_step2 INTEGER :: findx INTEGER :: i,j INTEGER :: STAT ! Remove any previous HDF files call SYSTEM('rm -rf *.hdf') ! Check file for compression - if so copy and uncompress INQUIRE(FILE=file1,EXIST=file_there) if( .not. file_there )then CALL DATE_AND_TIME(VALUES=DateTime) WRITE(dateString,'(''model1_'',i4.4,''-'',i2.2,''-'',i2.2,'':'',i2.2,'':'',i2.2,'':'',i2.2,''.hdf'')')& DateTime(1),DateTime(2),DateTime(3),DateTime(5),DateTime(6),& DateTime(7) file1in = TRIM(file1)//'.gz' INQUIRE(FILE=file1in,EXIST=file_there) if( .not. file_there )then write(errorString,'(''ERROR: File '',a,'' not found'')')& TRIM(file1) CALL Gbcs_Critical( .TRUE., & errorString,& 'NOAA_HDF_Get_AuX', & 'NOAA/NOAA_LoadAuxData.f90' ) endif command = 'cp -f '//TRIM(file1in)//' '//TRIM(dateString)//'.gz' call SYSTEM(command) command = 'gunzip -f '//TRIM(dateString)//'.gz' call SYSTEM(command) file1in = TRIM(dateString) else file1in = file1 endif INQUIRE(FILE=file2,EXIST=file_there) if( .not. file_there )then CALL DATE_AND_TIME(VALUES=DateTime) WRITE(dateString,'(''model2_'',i4.4,''-'',i2.2,''-'',i2.2,'':'',i2.2,'':'',i2.2,'':'',i2.2,''.hdf'')')& DateTime(1),DateTime(2),DateTime(3),DateTime(5),DateTime(6),& DateTime(7) file2in = TRIM(file2)//'.gz' INQUIRE(FILE=file2in,EXIST=file_there) if( .not. file_there )then write(errorString,'(''ERROR: File '',a,'' not found'')')& TRIM(file2) CALL Gbcs_Critical( .TRUE., & errorString,& 'NOAA_HDF_Get_AuX', & 'NOAA/NOAA_LoadAuxData.f90' ) endif command = 'cp -f '//TRIM(file2in)//' '//TRIM(dateString)//'.gz' call SYSTEM(command) command = 'gunzip -f '//TRIM(dateString)//'.gz' call SYSTEM(command) file2in = TRIM(dateString) else file2in = file2 endif CALL Read_HDF_File( file1in, nelems1, nlines1, nlevels1, pressure1, & height1, temperature1, rh1, ozone1, sfc_t1, sfc_rh1, sfc_press1, & tcwv1, wind_u1, wind_v1, height_1000mb_1, min_long1, max_long1, & long_step1, min_lat1, max_lat1, lat_step1, ozone_there1, & landseamask, two_m_temp1 ) CALL Read_HDF_File( file2in, nelems2, nlines2, nlevels2, pressure2, & height2, temperature2, rh2, ozone2, sfc_t2, sfc_rh2, sfc_press2, & tcwv2, wind_u2, wind_v2, height_1000mb_2, min_long2, max_long2, & long_step2, min_lat2, max_lat2, lat_step2, ozone_there2, & landseamask2, two_m_temp2 ) ! Check grid IF( min_long1 .ne. min_long2 .or. & max_long1 .ne. max_long2 .or. & long_step1 .ne. long_step2 .or. & min_lat1 .ne. min_lat2 .or. & max_lat1 .ne. max_lat2 .or. & lat_step1 .ne. lat_step2 .or. & nelems1 .ne. nelems2 .or. & nlines1 .ne. nlines2 .or. & nlevels1 .ne. nlevels2 )THEN CALL Gbcs_Critical( .TRUE., & 'Long/Lat grid mismatch',& 'NOAA_HDF_Get_Aux', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF IF( AUX%Model%Grid%NElems .ne. nelems1 .or. & AUX%Model%Grid%NLines .ne. nlines1 .or. & AUX%Model%Grid%NLevels .ne. nlevels1 .or. & AUX%Model%Grid%incr%lon .ne. long_step1 .or. & AUX%Model%Grid%incr%lat .ne. lat_step1 .or. & AUX%Model%Grid%LLCrnr%lat .ne. min_lat1 .or. & AUX%Model%Grid%LLCrnr%lon .ne. min_long1 )THEN print *,AUX%Model%Grid%NElems,nelems1,& AUX%Model%Grid%NLines,nlines1,& AUX%Model%Grid%NLevels,nlevels1,& AUX%Model%Grid%incr%lon,long_step1,& AUX%Model%Grid%incr%lat,lat_step1,& AUX%Model%Grid%LLCrnr%lat,min_lat1,& AUX%Model%Grid%LLCrnr%lon,min_long1 CALL Gbcs_Critical( .TRUE., & 'HDF Grid does not match grid in .INF file',& 'NOAA_HDF_Get_Aux', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Do time interpolation pressure1 = ((3.-hour)/3.)*pressure1 + (hour/3.)*pressure2 height1 = ((3.-hour)/3.)*height1 + (hour/3.)*height2 temperature1 = ((3.-hour)/3.)*temperature1 + (hour/3.)*temperature2 rh1 = ((3.-hour)/3.)*rh1 + (hour/3.)*rh2 ozone1 = ((3.-hour)/3.)*ozone1 + (hour/3.)*ozone2 sfc_t1 = ((3.-hour)/3.)*sfc_t1 + (hour/3.)*sfc_t2 sfc_rh1 = ((3.-hour)/3.)*sfc_rh1 + (hour/3.)*sfc_rh2 sfc_press1 = ((3.-hour)/3.)*sfc_press1 + (hour/3.)*sfc_press2 tcwv1 = ((3.-hour)/3.)*tcwv1 + (hour/3.)*tcwv2 wind_u1 = ((3.-hour)/3.)*wind_u1 + (hour/3.)*wind_u2 wind_v1 = ((3.-hour)/3.)*wind_v1 + (hour/3.)*wind_v2 height_1000mb_1 = ((3.-hour)/3.)*height_1000mb_1 + (hour/3.)*height_1000mb_2 two_m_temp1 = ((3.-hour)/3.)*two_m_temp1 + (hour/3.)*two_m_temp2 DO I=1,nlines1 DO J=1,nelems1 ! If any ice mask as iced IF( 2. .eq. landseamask(J,I) .or. 2 .eq. landseamask2(J,I) )THEN landseamask(J,I) = 2 ENDIF END DO END DO DEALLOCATE(landseamask2) ! Forecast Model 3D Fields !---------------------------- IF ( Gbcs_Verbosity > 2 ) & WRITE(Gbcs_Log_File_Unit,*)& ' Load Forecast Model 3-D Fields (HDF file)' AUX%Model%N_3D_Model_Fields = 0 ! Pressure findx = Gbcs_Map_P_3D IF (ALLOCATED(AUX%Model%Atm(Gbcs_Map_P_3D)%d)) THEN DEALLOCATE(AUX%Model%Atm(Gbcs_Map_P_3D)%d, STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Deallocating P 3D', & 'NOAA_HDF_Get_Aux', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF ALLOCATE(AUX%Model%Atm(Gbcs_Map_P_3D)%d(nlevels1+1,NElems1,NLines1), & STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating P 3D', & 'NOAA_HDF_Get_Aux', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Load pressure level data from one of the 3-D fields. DO j = 1,nlines1 DO i = 1,nelems1 AUX%Model%Atm(findx)%d(2:nlevels1+1,i,j) = pressure1 END DO END DO AUX%Model%Atm(findx)%d(1,:,:) = TOA_PRESSURE STAT = Register_3D_Field(AUX%Model, Gbcs_Map_P_3D, Gbcs_Units_hPa) ! Temperature findx = Gbcs_Map_T_3D IF (.NOT.ALLOCATED(AUX%Model%Atm(findx)%d)) THEN ALLOCATE(AUX%Model%Atm(findx)%d(nlevels1+1,NElems1,NLines1), & STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating T 3D', & 'NOAA_HDF_Get_Aux', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF do i=2,nlevels1+1 AUX%Model%Atm(findx)%d(i,:,:) = temperature1(:,:,i-1) end do AUX%Model%Atm(findx)%d(1,:,:) = TOA_TEMPERATURE STAT = Register_3D_Field(AUX%Model, Gbcs_Map_T_3D, Gbcs_Units_Kelvin) ! RH findx = Gbcs_Map_WV_3D IF (.NOT.ALLOCATED(AUX%Model%Atm(findx)%d)) THEN ALLOCATE(AUX%Model%Atm(findx)%d(nlevels1+1,NElems1,NLines1), & STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating WV 3D', & 'NOAA_HDF_Get_Aux', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF do i=2,nlevels1+1 AUX%Model%Atm(findx)%d(i,:,:) = rh1(:,:,i-1) end do AUX%Model%Atm(findx)%d(1,:,:) = TOA_RH STAT = Register_3D_Field(AUX%Model, Gbcs_Map_WV_3D, Gbcs_Units_RH) ! Ozone findx = Gbcs_Map_O3_3D IF (.NOT.ALLOCATED(AUX%Model%Atm(findx)%d)) THEN ALLOCATE(AUX%Model%Atm(findx)%d(nlevels1+1,NElems1,NLines1), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating O3 3D', & 'NOAA_HDF_Get_Aux', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF ! Set NLevels + 1 due to TOA values AUX%Model%Grid%NLevels = nlevels1+1 ! Need everything above to be setup so now get O3 CALL Get_Ozone_climatology( IMG%Mean_Time%Month, AUX ) ! Now add in model values do i=2,nlevels1+1 IF( ozone_there1(i-1) )THEN AUX%Model%Atm(findx)%d(i,:,:) = ozone1(:,:,i-1) ENDIF end do STAT = Register_3D_Field(AUX%Model, Gbcs_Map_O3_3D, Gbcs_Units_ppmv) ! Forecast Model 2D Fields !---------------------------- IF ( Gbcs_Verbosity > 2 ) & WRITE(Gbcs_Log_File_Unit,*)' Load Forecast Model 2-D Fields (HDF file)' AUX%Model%N_2D_Model_Fields = 0 ! 1000mb height findx = Gbcs_Map_z1000 IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems1,NLines1), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating z1000', & 'NOAA_HDF_Get_Aux', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF AUX%Model%Srf(findx)%d = height_1000mb_1 STAT = Register_2D_Field(AUX%Model, Gbcs_Map_z1000, Gbcs_Units_m) ! Orography findx = Gbcs_Map_OROG IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems1,NLines1), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating OROG', & 'NOAA_HDF_Get_Aux', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF AUX%Model%Srf(findx)%d = 0.0 STAT = Register_2D_Field(AUX%Model, Gbcs_Map_OROG, Gbcs_Units_m) ! Surface temperature ! Check difference between SST/Surface etc findx = Gbcs_Map_T_SKIN IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems1,NLines1), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating T Skin', & 'NOAA_HDF_Get_Aux', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF AUX%Model%Srf(findx)%d = sfc_t1 STAT = Register_2D_Field(AUX%Model, Gbcs_Map_T_SKIN, Gbcs_Units_Kelvin) ! Total Column Water Vapour findx = Gbcs_Map_TCWV IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems1,NLines1), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating TCWV', & 'NOAA_HDF_Get_Aux', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF AUX%Model%Srf(findx)%d = tcwv1 STAT = Register_2D_Field(AUX%Model, Gbcs_Map_TCWV, Gbcs_Units_Kg_Per_m2) ! Wind - U findx = Gbcs_Map_U_SURF IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems1,NLines1), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating U Surf', & 'NOAA_HDF_Get_Aux', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF AUX%Model%Srf(findx)%d = wind_u1 STAT = Register_2D_Field(AUX%Model, Gbcs_Map_U_Surf, Gbcs_Units_m_per_s) ! Wind - V findx = Gbcs_Map_V_SURF IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems1,NLines1), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating V Surf', & 'NOAA_HDF_Get_Aux', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF AUX%Model%Srf(findx)%d = wind_v1 STAT = Register_2D_Field(AUX%Model, Gbcs_Map_V_Surf, Gbcs_Units_m_per_s) ! Wind - 10m findx = Gbcs_Map_Wind_10M IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems1,NLines1), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating 10m Wind', & 'NOAA_HDF_Get_Aux', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF do j=1,nlines1 do i=1,nelems1 if( 25 .lt. AUX%Model%Srf(Gbcs_Map_U_SURF)%d(i,j) .or. & 25 .lt. AUX%Model%Srf(Gbcs_Map_V_SURF)%d(i,j) )then if( 25 .lt. AUX%Model%Srf(Gbcs_Map_U_SURF)%d(i,j) .and. & 25 .lt. AUX%Model%Srf(Gbcs_Map_V_SURF)%d(i,j) )then AUX%Model%Srf(findx)%d(i,j) = 8 else if( 25 .lt. AUX%Model%Srf(Gbcs_Map_U_SURF)%d(i,j) )then AUX%Model%Srf(findx)%d(i,j) = & ABS(AUX%Model%Srf(Gbcs_Map_V_SURF)%d(i,j)) else if( 25 .lt. AUX%Model%Srf(Gbcs_Map_V_SURF)%d(i,j) )then AUX%Model%Srf(findx)%d(i,j) = & ABS(AUX%Model%Srf(Gbcs_Map_U_SURF)%d(i,j)) endif else AUX%Model%Srf(findx)%d(i,j) = sqrt(& AUX%Model%Srf(Gbcs_Map_V_SURF)%d(i,j)**2 + & AUX%Model%Srf(Gbcs_Map_U_SURF)%d(i,j)**2) endif end do end do STAT = Register_2D_Field(AUX%Model, Gbcs_Map_Wind_10M, Gbcs_Units_m_per_s) ! Get surface pressure findx = Gbcs_Map_P_Surf IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems1,NLines1), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating P Surf', & 'NOAA_HDF_Get_Aux', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF ! surface pressure AUX%Model%Srf(findx)%d = sfc_press1 STAT = Register_2D_Field(AUX%Model, Gbcs_Map_P_Surf, Gbcs_Units_hPa) ! Total surface RH findx = Gbcs_Map_WV_2M IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems1,NLines1), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating 2m WV', & 'NOAA_HDF_Get_Aux', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF AUX%Model%Srf(findx)%d = sfc_rh1 STAT = Register_2D_Field(AUX%Model, Gbcs_Map_WV_2M, Gbcs_Units_RH) ! 2m temperature - air temperature findx = Gbcs_Map_T_2M IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems1,NLines1), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating 2m WV', & 'NOAA_HDF_Get_Aux', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF AUX%Model%Srf(findx)%d = two_m_temp1 STAT = Register_2D_Field(AUX%Model, Gbcs_Map_T_2M, Gbcs_Units_kelvin) DEALLOCATE( pressure1, height1, temperature1, rh1, ozone1, sfc_t1, & sfc_rh1, sfc_press1, tcwv1, wind_u1, wind_v1, height_1000mb_1, & pressure2, height2, temperature2, rh2, ozone2, sfc_t2, sfc_rh2, & sfc_press2, tcwv2, wind_u2, wind_v2, height_1000mb_2, ozone_there1, & ozone_there2, two_m_temp1, two_m_temp2 ) END SUBROUTINE NOAA_HDF_Get_AuX SUBROUTINE Read_HDF_File( filename, nelems, nlines, nlevels, & pressure, height, temperature, rh, ozone, sfc_t, sfc_rh, sfc_press, & tcwv, wind_u, wind_v, height_1000mb, min_Long, max_Long, & long_step, min_Lat, max_Lat, lat_step, ozone_there, ocean_sfc, & two_m_temp ) CHARACTER(LEN=*), INTENT(IN) :: filename INTEGER, INTENT(OUT) :: nelems INTEGER, INTENT(OUT) :: nlines INTEGER, INTENT(OUT) :: nlevels REAL, ALLOCATABLE, INTENT(OUT) :: pressure(:) REAL, ALLOCATABLE, INTENT(OUT) :: height(:,:,:) REAL, ALLOCATABLE, INTENT(OUT) :: temperature(:,:,:) REAL, ALLOCATABLE, INTENT(OUT) :: rh(:,:,:) REAL, ALLOCATABLE, INTENT(OUT) :: ozone(:,:,:) REAL, ALLOCATABLE, INTENT(OUT) :: sfc_t(:,:) REAL, ALLOCATABLE, INTENT(OUT) :: sfc_rh(:,:) REAL, ALLOCATABLE, INTENT(OUT) :: sfc_press(:,:) REAL, ALLOCATABLE, INTENT(OUT) :: tcwv(:,:) REAL, ALLOCATABLE, INTENT(OUT) :: wind_u(:,:) REAL, ALLOCATABLE, INTENT(OUT) :: wind_v(:,:) REAL, ALLOCATABLE, INTENT(OUT) :: height_1000mb(:,:) REAL, INTENT(OUT) :: min_Long, max_Long, long_step REAL, INTENT(OUT) :: min_Lat, max_Lat, lat_step LOGICAL, ALLOCATABLE, INTENT(OUT) :: ozone_there(:) REAL, ALLOCATABLE, INTENT(OUT) :: ocean_sfc(:,:) REAL, ALLOCATABLE, INTENT(OUT) :: two_m_temp(:,:) ! Local variables INTEGER :: sd_id INTEGER :: sds_id INTEGER :: attr_index INTEGER :: status INTEGER :: ndatasets INTEGER :: nattributes INTEGER, parameter :: DFACC_RDONLY = 1 INTEGER :: I, J CHARACTER(LEN=256) :: sds_name INTEGER :: rank,dim_sizes(32),data_type,n_attrs LOGICAL :: set_latlon REAL :: min_Longitude, max_Longitude, longitude_step REAL :: min_Latitude, max_Latitude, latitude_step INTEGER :: nlines_in, nelems_in INTEGER :: level_type, variable_number REAL :: levels_value REAL :: temp LOGICAL :: ok INTEGER :: start(2), edges(2), stride(2) INTEGER :: Indx REAL(fp_kind), ALLOCATABLE :: pressure_temp(:) REAL(fp_kind), ALLOCATABLE :: rh_temp(:) REAL(fp_kind), ALLOCATABLE :: temperature_temp(:) REAL(fp_kind), ALLOCATABLE :: height_temp(:) REAL(fp_kind), ALLOCATABLE :: water_vapour_pressure(:) REAL(fp_kind), ALLOCATABLE :: mass_density(:) REAL(fp_kind), ALLOCATABLE :: mixing_ratio(:) REAL, ALLOCATABLE :: land(:,:) REAL, ALLOCATABLE :: ice(:,:) INTEGER, PARAMETER :: PROFILE_TYPE = 100 INTEGER, PARAMETER :: TEMPERATURE_TYPE = 11 INTEGER, PARAMETER :: HEIGHT_TYPE = 7 INTEGER, PARAMETER :: OZONE_TYPE = 154 INTEGER, PARAMETER :: RH_TYPE = 52 INTEGER, PARAMETER :: TWO_AND_TEN_M_TYPE = 105 INTEGER, PARAMETER :: WIND_U_TYPE = 33 INTEGER, PARAMETER :: WIND_V_TYPE = 34 INTEGER, PARAMETER :: SFC_TYPE = 1 INTEGER, PARAMETER :: PRESSURE_TYPE = 1 INTEGER, PARAMETER :: LAND_TYPE = 81 INTEGER, PARAMETER :: ICE_TYPE = 91 ! Tests to make sure we have everything LOGICAL :: temp_profile_logical = .FALSE. LOGICAL :: height_profile_logical = .FALSE. LOGICAL :: ozone_profile_logical = .FALSE. LOGICAL :: rh_profile_logical = .FALSE. LOGICAL :: wind_u_logical = .FALSE. LOGICAL :: wind_v_logical = .FALSE. LOGICAL :: temp_sfc_logical = .FALSE. LOGICAL :: rh_2m_logical = .FALSE. LOGICAL :: press_sfc_logical = .FALSE. LOGICAL :: land_sfc_logical = .FALSE. LOGICAL :: ice_sfc_logical = .FALSE. LOGICAL :: two_m_temp_logical = .FALSE. ! HDF routines INTEGER :: sfstart INTEGER :: sffinfo INTEGER :: sfselect INTEGER :: sfginfo INTEGER :: sffattr INTEGER :: sfrattr INTEGER :: sfendacc INTEGER :: sfrdata INTEGER :: sfend set_latlon = .FALSE. level_type = 0 variable_number = 0 ! get indexes of things we need ! Get number of levels, nelem, nlines sd_id = sfstart(filename, DFACC_RDONLY) CALL checkHDF(sd_id,"sfstart") ! Get number of entries status = sffinfo(sd_id,ndatasets,nattributes) CALL checkHDF(sd_id,"sffinfo") ! Loop round variables and find all TMP profile data ! Level_Type = 100 Variable_Number = 11 nlevels = 0 DO I=1,ndatasets ! select dataset - index from 0 sds_id = sfselect(sd_id,I-1) CALL checkHDF(sds_id,"sfselect") ! Get name/size etc. status = sfginfo(sds_id,sds_name,rank,dim_sizes,data_type,n_attrs) CALL checkHDF(status,"sfginfo") IF( 2 .ne. rank )THEN CALL Gbcs_Critical( .TRUE., & 'Rank != 2',& 'Read_HDF_File', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF nelems_in = dim_sizes(1) nlines_in = dim_sizes(2) ! Get attributes - get min/max lat/lon attr_index = sffattr(sds_id,"Min_Longitude") CALL checkHDF(attr_index,"sffattr") status = sfrattr(sds_id, attr_index, min_longitude) CALL checkHDF(status,"sfrattr") attr_index = sffattr(sds_id,"Max_Longitude") CALL checkHDF(attr_index,"sffattr") status = sfrattr(sds_id, attr_index, max_longitude) CALL checkHDF(status,"sfrattr") attr_index = sffattr(sds_id,"Longitude_Step") CALL checkHDF(attr_index,"sffattr") status = sfrattr(sds_id, attr_index, longitude_step) CALL checkHDF(status,"sfrattr") attr_index = sffattr(sds_id,"Min_Latitude") CALL checkHDF(attr_index,"sffattr") status = sfrattr(sds_id, attr_index, max_latitude) CALL checkHDF(status,"sfrattr") attr_index = sffattr(sds_id,"Max_Latitude") CALL checkHDF(attr_index,"sffattr") status = sfrattr(sds_id, attr_index, min_latitude) CALL checkHDF(status,"sfrattr") attr_index = sffattr(sds_id,"Latitude_Step") CALL checkHDF(attr_index,"sffattr") status = sfrattr(sds_id, attr_index, latitude_step) CALL checkHDF(status,"sfrattr") ! ! Fix 'bug' in HDF files as currently written ! IF( max_latitude .lt. min_latitude .and. latitude_step .gt. 0 )THEN ! latitude_step = -latitude_step ! ENDIF ! IF( max_latitude .gt. min_latitude .and. latitude_step .lt. 0 )THEN ! latitude_step = -latitude_step ! ENDIF ! Check grid size IF( set_latlon )THEN IF( min_long .ne. min_longitude .or. & max_long .ne. max_longitude .or. & long_step .ne. longitude_step .or. & nelems .ne. nelems_in .or. & min_lat .ne. min_latitude .or. & max_lat .ne. max_latitude .or. & lat_step .ne. latitude_step .or. & nlines .ne. nlines_in )THEN CALL Gbcs_Critical( .TRUE., & 'Mismatch in lat/lon grid definition',& 'Read_HDF_File', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ELSE set_latlon = .TRUE. min_long = min_longitude max_long = max_longitude long_step = longitude_step nelems = nelems_in min_lat = min_latitude max_lat = max_latitude lat_step = latitude_step nlines = nlines_in ENDIF ! Now find TMP levels attr_index = sffattr(sds_id,"Level_Type") CALL checkHDF(attr_index,"sffattr") status = sfrattr(sds_id, attr_index, level_type) CALL checkHDF(status,"sfrattr") attr_index = sffattr(sds_id,"Variable_Number") CALL checkHDF(attr_index,"sffattr") status = sfrattr(sds_id, attr_index, variable_number) CALL checkHDF(status,"sfrattr") IF( level_type .eq. PROFILE_TYPE .and. & variable_number .eq. TEMPERATURE_TYPE )THEN nlevels = nlevels + 1 ENDIF ! terminate access to current data set status = sfendacc(sds_id) CALL checkHDF(status,"sfendacc") END DO ! Allocate output arrays ALLOCATE(pressure(nlevels),& ozone_there(nlevels),& height(nelems,nlines,nlevels),& temperature(nelems,nlines,nlevels),& rh(nelems,nlines,nlevels),& ozone(nelems,nlines,nlevels),& sfc_t(nelems,nlines),& sfc_rh(nelems,nlines),& sfc_press(nelems,nlines),& two_m_temp(nelems,nlines),& tcwv(nelems,nlines),& wind_u(nelems,nlines),& wind_v(nelems,nlines),& land(nelems,nlines),& ice(nelems,nlines),& ocean_sfc(nelems,nlines),& height_1000mb(nelems,nlines),& pressure_temp(nlevels),& rh_temp(nlevels),& temperature_temp(nlevels),& height_temp(nlevels),& water_vapour_pressure(nlevels),& mass_density(nlevels),& mixing_ratio(nlevels),& STAT=status) IF( 0 .ne. status )THEN CALL Gbcs_Critical( .TRUE., & 'Allocation Error',& 'Read_HDF_File', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ozone = 0. ozone_there = .FALSE. ! Get pressure levels J = 0 DO I=1,ndatasets ! select dataset - index from 0 sds_id = sfselect(sd_id,I-1) CALL checkHDF(sds_id,"sfselect") ! Now find TMP levels attr_index = sffattr(sds_id,"Level_Type") CALL checkHDF(attr_index,"sffattr") status = sfrattr(sds_id, attr_index, level_type) CALL checkHDF(status,"sfrattr") attr_index = sffattr(sds_id,"Variable_Number") CALL checkHDF(attr_index,"sffattr") status = sfrattr(sds_id, attr_index, variable_number) CALL checkHDF(status,"sfrattr") IF( level_type .eq. PROFILE_TYPE .and. & variable_number .eq. TEMPERATURE_TYPE )THEN J = J + 1 attr_index = sffattr(sds_id,"Levels_Value") CALL checkHDF(attr_index,"sffattr") status = sfrattr(sds_id, attr_index, levels_value) pressure(J) = levels_value ENDIF ! terminate access to current data set status = sfendacc(sds_id) CALL checkHDF(status,"sfendacc") END DO IF( J .ne. nlevels )THEN CALL Gbcs_Critical( .TRUE., & 'J ne nlevels',& 'Read_HDF_File', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Sort pressure array OK = .FALSE. DO WHILE(.not.OK) OK = .TRUE. DO I=1,nlevels-1 IF( pressure(I+1) .lt. pressure(I) )THEN temp = pressure(I+1) pressure(I+1) = pressure(I) pressure(I) = temp OK = .FALSE. ENDIF END DO END DO temp_profile_logical = .FALSE. height_profile_logical = .FALSE. ozone_profile_logical = .FALSE. rh_profile_logical = .FALSE. wind_u_logical = .FALSE. wind_v_logical = .FALSE. temp_sfc_logical = .FALSE. rh_2m_logical = .FALSE. press_sfc_logical = .FALSE. land_sfc_logical = .FALSE. ice_sfc_logical = .FALSE. two_m_temp_logical = .FALSE. ! Now fill arrays start(1) = 0 start(2) = 0 edges(1) = nelems edges(2) = nlines stride(1) = 1 stride(2) = 1 DO I=1,ndatasets ! select dataset - index from 0 sds_id = sfselect(sd_id,I-1) CALL checkHDF(sds_id,"sfselect") attr_index = sffattr(sds_id,"Level_Type") CALL checkHDF(attr_index,"sffattr") status = sfrattr(sds_id, attr_index, level_type) CALL checkHDF(status,"sfrattr") attr_index = sffattr(sds_id,"Variable_Number") CALL checkHDF(attr_index,"sffattr") status = sfrattr(sds_id, attr_index, variable_number) CALL checkHDF(status,"sfrattr") IF( PROFILE_TYPE .eq. Level_Type .and. & Variable_Number .eq. TEMPERATURE_TYPE )THEN ! TMP attr_index = sffattr(sds_id,"Levels_Value") CALL checkHDF(attr_index,"sffattr") status = sfrattr(sds_id, attr_index, levels_value) Indx = Get_Pressure_Index_HDF( nlevels, pressure, levels_value ) ! Get values status = sfrdata(sds_id,start,stride,edges,temperature(:,:,Indx)) CALL checkHDF(status,"sfrdata") IF( 0 .eq. SUM(temperature(:,:,Indx)) )THEN CALL Gbcs_Critical(.TRUE.,& 'Temperature profile entry is all zero',& 'Read_HDF_File','NOAA/NOAA_LoadAuxData.f90') ENDIF temp_profile_logical = .TRUE. ELSE IF( PROFILE_TYPE .eq. Level_Type .and. & Variable_Number .eq. HEIGHT_TYPE )THEN ! Height attr_index = sffattr(sds_id,"Levels_Value") CALL checkHDF(attr_index,"sffattr") status = sfrattr(sds_id, attr_index, levels_value) Indx = Get_Pressure_Index_HDF( nlevels, pressure, levels_value ) ! Get values status = sfrdata(sds_id,start,stride,edges,height(:,:,Indx)) CALL checkHDF(status,"sfrdata") IF( 0 .eq. SUM(height(:,:,Indx)) )THEN CALL Gbcs_Critical(.TRUE.,& 'Height profile entry is all zero',& 'Read_HDF_File','NOAA/NOAA_LoadAuxData.f90') ENDIF height_profile_logical = .TRUE. ELSE IF( PROFILE_TYPE .eq. Level_Type .and. & Variable_Number .eq. OZONE_TYPE )THEN ! Ozone attr_index = sffattr(sds_id,"Levels_Value") CALL checkHDF(attr_index,"sffattr") status = sfrattr(sds_id, attr_index, levels_value) Indx = Get_Pressure_Index_HDF( nlevels, pressure, levels_value ) ! Get values status = sfrdata(sds_id,start,stride,edges,ozone(:,:,Indx)) ozone_there(Indx) = .TRUE. CALL checkHDF(status,"sfrdata") IF( 0 .eq. SUM(ozone(:,:,Indx)) )THEN CALL Gbcs_Critical(.TRUE.,& 'Ozone profile entry is all zero',& 'Read_HDF_File','NOAA/NOAA_LoadAuxData.f90') ENDIF ozone_profile_logical = .TRUE. ELSE IF( PROFILE_TYPE .eq. Level_Type .and. & Variable_Number .eq. RH_TYPE )THEN ! RH attr_index = sffattr(sds_id,"Levels_Value") CALL checkHDF(attr_index,"sffattr") status = sfrattr(sds_id, attr_index, levels_value) Indx = Get_Pressure_Index_HDF( nlevels, pressure, levels_value ) ! Get values status = sfrdata(sds_id,start,stride,edges,rh(:,:,Indx)) CALL checkHDF(status,"sfrdata") IF( 0 .eq. SUM(rh(:,:,Indx)) )THEN CALL Gbcs_Critical(.TRUE.,& 'RH profile entry is all zero',& 'Read_HDF_File','NOAA/NOAA_LoadAuxData.f90') ENDIF rh_profile_logical = .TRUE. ELSE IF( TWO_AND_TEN_M_TYPE .eq. Level_Type .and. & Variable_Number .eq. WIND_U_TYPE )THEN ! Wind U ! Get values status = sfrdata(sds_id,start,stride,edges,wind_u(:,:)) CALL checkHDF(status,"sfrdata") IF( 0 .eq. SUM(wind_u(:,:)) )THEN CALL Gbcs_Critical(.TRUE.,& 'Wind U entry is all zero',& 'Read_HDF_File','NOAA/NOAA_LoadAuxData.f90') ENDIF wind_u_logical = .TRUE. ELSE IF( TWO_AND_TEN_M_TYPE .eq. Level_Type .and. & Variable_Number .eq. WIND_V_TYPE )THEN ! Wind U ! Get values status = sfrdata(sds_id,start,stride,edges,wind_v(:,:)) CALL checkHDF(status,"sfrdata") IF( 0 .eq. SUM(wind_v(:,:)) )THEN CALL Gbcs_Critical(.TRUE.,& 'Wind V entry is all zero',& 'Read_HDF_File','NOAA/NOAA_LoadAuxData.f90') ENDIF wind_v_logical = .TRUE. ELSE IF( SFC_TYPE .eq. Level_Type .and. & Variable_Number .eq. TEMPERATURE_TYPE )THEN ! Sfc temperature ! Get values status = sfrdata(sds_id,start,stride,edges,sfc_t(:,:)) CALL checkHDF(status,"sfrdata") IF( 0 .eq. SUM(sfc_t(:,:)) )THEN CALL Gbcs_Critical(.TRUE.,& 'Surface Temperature entry is all zero',& 'Read_HDF_File','NOAA/NOAA_LoadAuxData.f90') ENDIF temp_sfc_logical = .TRUE. ELSE IF( TWO_AND_TEN_M_TYPE .eq. Level_Type .and. & Variable_Number .eq. RH_TYPE )THEN ! 2m RH ! Get values status = sfrdata(sds_id,start,stride,edges,sfc_rh(:,:)) CALL checkHDF(status,"sfrdata") IF( 0 .eq. SUM(sfc_rh(:,:)) )THEN CALL Gbcs_Critical(.TRUE.,& '2m RH entry is all zero',& 'Read_HDF_File','NOAA/NOAA_LoadAuxData.f90') ENDIF rh_2m_logical = .TRUE. ELSE IF( TWO_AND_TEN_M_TYPE .eq. Level_Type .and. & Variable_Number .eq. TEMPERATURE_TYPE )THEN ! 2m Temperature ! Get values status = sfrdata(sds_id,start,stride,edges,two_m_temp(:,:)) CALL checkHDF(status,"sfrdata") IF( 0 .eq. SUM(sfc_rh(:,:)) )THEN CALL Gbcs_Critical(.TRUE.,& '2m Temperature entry is all zero',& 'Read_HDF_File','NOAA/NOAA_LoadAuxData.f90') ENDIF two_m_temp_logical = .TRUE. ELSE IF( SFC_TYPE .eq. Level_Type .and. & Variable_Number .eq. PRESSURE_TYPE )THEN ! Sfc pressure ! Get values status = sfrdata(sds_id,start,stride,edges,sfc_press(:,:)) CALL checkHDF(status,"sfrdata") IF( 0 .eq. SUM(sfc_press(:,:)) )THEN CALL Gbcs_Critical(.TRUE.,& 'Surface Pressure entry is all zero',& 'Read_HDF_File','NOAA/NOAA_LoadAuxData.f90') ENDIF press_sfc_logical = .TRUE. ELSE IF( SFC_TYPE .eq. Level_Type .and. & Variable_Number .eq. LAND_TYPE )THEN ! Sfc pressure ! Get values status = sfrdata(sds_id,start,stride,edges,land(:,:)) CALL checkHDF(status,"sfrdata") IF( 0 .eq. SUM(land(:,:)) )THEN CALL Gbcs_Critical(.TRUE.,& 'Land Mask entry is all zero',& 'Read_HDF_File','NOAA/NOAA_LoadAuxData.f90') ENDIF land_sfc_logical = .TRUE. ELSE IF( SFC_TYPE .eq. Level_Type .and. & Variable_Number .eq. ICE_TYPE )THEN ! Sfc pressure ! Get values status = sfrdata(sds_id,start,stride,edges,ice(:,:)) CALL checkHDF(status,"sfrdata") IF( 0 .eq. SUM(ice(:,:)) )THEN CALL Gbcs_Critical(.TRUE.,& 'Ice Mask entry is all zero',& 'Read_HDF_File','NOAA/NOAA_LoadAuxData.f90') ENDIF ice_sfc_logical = .TRUE. ENDIF ! terminate access to current data set status = sfendacc(sds_id) CALL checkHDF(status,"sfendacc") END DO ! Close file status = sfend(sd_id) CALL checkHDF(status,"sfend") I = 0 IF( .not. temp_profile_logical )THEN I = 1 WRITE(*,'('' ERROR: Missing temperature profile data'')') ENDIF IF( .not. height_profile_logical )THEN I = 1 WRITE(*,'('' ERROR: Missing height profile data'')') ENDIF IF( .not. ozone_profile_logical )THEN I = 1 WRITE(*,'('' ERROR: Missing ozone profile data'')') ENDIF IF( .not. rh_profile_logical )THEN I = 1 WRITE(*,'('' ERROR: Missing RH profile data'')') ENDIF IF( .not. wind_u_logical )THEN I = 1 WRITE(*,'('' ERROR: Missing Wind U data'')') ENDIF IF( .not. wind_v_logical )THEN I = 1 WRITE(*,'('' ERROR: Missing Wind V data'')') ENDIF IF( .not. temp_sfc_logical )THEN I = 1 WRITE(*,'('' ERROR: Missing Surface Temperature data'')') ENDIF IF( .not. rh_2m_logical )THEN I = 1 WRITE(*,'('' ERROR: Missing 2m RH data'')') ENDIF IF( .not. two_m_temp_logical )THEN I = 1 WRITE(*,'('' ERROR: Missing 2m Temperature data'')') ENDIF IF( .not. press_sfc_logical )THEN I = 1 WRITE(*,'('' ERROR: Missing Surface Pressure data'')') ENDIF IF( .not. land_sfc_logical )THEN I = 1 WRITE(*,'('' ERROR: Missing Land Mask data'')') ENDIF IF( .not. ice_sfc_logical )THEN I = 1 WRITE(*,'('' ERROR: Missing Ice Mask data'')') ENDIF IF( 1 .eq. I )THEN call Gbcs_Critical(.TRUE.,'Missing data entries from HDF Model file',& 'Read_HDF_File','NOAA/NOAA_LoadAuxData.f90') ENDIF ! Derive ocean_sfc - 0 = sea, 1 = land 2= ice DO I=1,nlines DO J=1,nelems IF( 0. .eq. land(J,I) .and. 0. .eq. ice(J,I) )THEN ocean_sfc(J,I) = 0. ELSE IF( 0. .lt. land(J,I) )THEN ocean_sfc(J,I) = 1. ELSE ocean_sfc(J,I) = 2. ENDIF END DO END DO ! Deallocate land and ice fractions DEALLOCATE(land,ice) ! Now get derived info ! 1000mb height Indx = Get_Pressure_Index_HDF( nlevels, pressure, 1000. ) height_1000mb = height(:,:,Indx) ! Get TCWV pressure_temp = pressure DO I=1,nlines DO J=1,nelems rh_temp = rh(j,i,:) temperature_temp = temperature(j,i,:) height_temp = height(j,i,:) mixing_ratio = RH_to_MR_Units(pressure_temp,temperature_temp,rh_temp,& Message_Log='/dev/null') water_vapour_pressure = MR_to_PP_Units(pressure_temp,mixing_ratio,& Molecule_ID=1,Message_Log='/dev/null') ! Get mass density mass_density = PP_to_MD_Units(water_vapour_pressure,temperature_temp,& Molecule_ID=1,Message_log='/dev/null') ! Now integrate over height to get TCWV ! Get TCWV kg/m^3 tcwv(J,I) = integrate(nlevels,height_temp,mass_density)/1000. END DO END DO DEALLOCATE(pressure_temp,rh_temp,temperature_temp,height_temp,& mixing_ratio,water_vapour_pressure,mass_density) END SUBROUTINE Read_HDF_File real(fp_kind) function integrate( n, x, y ) integer, intent(in) :: n real(fp_kind), intent(in) :: x(n) real(fp_kind), intent(in) :: y(n) integer :: i real(fp_kind) :: sum sum = 0. do i=1,n-1 sum = sum - 0.5*(y(i+1)+y(i))*(x(i+1)-x(i)) end do integrate = sum return end function integrate INTEGER FUNCTION Get_Pressure_Index_HDF( nlevels, pressure, value )& RESULT(Indx) INTEGER, INTENT(IN) :: nlevels REAL, INTENT(IN) :: pressure(nlevels) REAL, INTENT(IN) :: value ! Local variables INTEGER :: I Indx = -1 Loop: DO I=1,nlevels IF( pressure(I) .eq. value )THEN Indx = I EXIT Loop END IF END DO Loop IF( -1 .eq. Indx )THEN print*,value,nlevels print*,pressure CALL Gbcs_Critical( .TRUE., & 'Index not found',& 'Get_Pressure_Index_HDF', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END FUNCTION Get_Pressure_Index_HDF ! Fill AUX arrays from operational version of HDF read SUBROUTINE NOAA_HDF_Get_Aux_OPER( file1, file2, hour, AUX, IMG ) ! INPUT CHARACTER(LEN=*), INTENT(IN) :: file1 CHARACTER(LEN=*), INTENT(IN) :: file2 REAL, INTENT(IN) :: hour TYPE(Aux_Data), INTENT(INOUT) :: AUX TYPE(Imagery), INTENT(IN) :: IMG integer, parameter :: onlevels = 18 integer, parameter :: onlines = 181 integer, parameter :: onelems = 360 ! Profiles of temperature on levels (K*100) integer:: lev_t(onlines,onelems,onlevels) ! Profiles of relative humidity on levels (%) integer:: lev_rh(onlines,onelems,onlevels) ! Altitude (m) at 1000.0 mB integer:: z1(onlines,onelems) ! Surface temperature (K) integer:: surf_t(onlines,onelems) ! Surface relative humidity (%) integer:: surf_rh(onlines,onelems) ! Surface windspeed, u (m / s) integer:: surf_u(onlines,onelems) ! Surface windspeed, v (m / s) integer:: surf_v(onlines,onelems) ! TCWV integer:: tcwv(onlines,onelems) integer :: grid(onlines,onelems) integer :: grid1(onlines,onelems) real :: data(onlines,onelems) real :: data1(onlines,onelems) integer :: ii, jj integer :: findx integer :: fieldid integer :: stat integer :: nlevels integer :: nelems integer :: nlines logical :: file_there character(len=256) :: file1in character(len=256) :: file2in character(len=256) :: command INTEGER :: DateTime(8) CHARACTER(LEN=80) :: dateString CHARACTER(LEN=6) :: hms integer, parameter :: oper_total_levels = onlevels+1 REAL(KIND=GbcsReal), parameter, dimension(oper_total_levels) :: & oper_Plevels = & (/ 0.005, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0,& & 450.0, 500.0, 550.0, 600.0, 650.0, 700.0, 750.0, 800.0, 850.0,& & 900.0, 950.0, 1000.0 /) ! Ozone mixing ratio (US Std Atm) ! Allocate integer output from IGGET REAL(KIND=GbcsReal), parameter, dimension(oper_total_levels) :: & oper_Olevels = & (/ 0.57, 0.4622, 0.2923, 0.1643, 0.09905, 0.06395, 0.05204, & & 0.04438,0.03972, 0.03720, 0.03470, 0.03360, 0.03319, 0.03276, & & 0.03222, 0.03075, 0.02928, 0.02810, 0.02691 /) REAL(KIND=GbcsReal), dimension(oper_total_levels) :: Plevels REAL(KIND=GbcsReal), dimension(oper_total_levels) :: Olevels ! extracted parameter list integer, parameter :: number_fields = 2*onlevels+6 character(len=20) :: oper_fieldlist(number_fields) DATA oper_fieldlist/'TMP 150 mb','TMP 200 mb','TMP 250 mb','TMP 300 mb', & 'TMP 350 mb','TMP 400 mb','TMP 450 mb','TMP 500 mb', & 'TMP 550 mb','TMP 600 mb','TMP 650 mb','TMP 700 mb', & 'TMP 750 mb','TMP 800 mb','TMP 850 mb','TMP 900 mb', & 'TMP 950 mb','TMP 1000 mb', & 'RH 150 mb','RH 200 mb','RH 250 mb','RH 300 mb', & 'RH 350 mb','RH 400 mb','RH 450 mb','RH 500 mb', & 'RH 550 mb','RH 600 mb','RH 650 mb','RH 700 mb', & 'RH 750 mb','RH 800 mb','RH 850 mb','RH 900 mb', & 'RH 950 mb','RH 1000 mb', & 'HGT 1000 mb', 'PWAT atmos col', 'TMP sfc','RH 2 m above gnd', & 'UGRD 10 m above gnd', 'VGRD 10 m above gnd'/ ! Remove any previous HDF files call SYSTEM('rm -rf *.hdf') ! Check file for compression - if so copy and uncompress INQUIRE(FILE=file1,EXIST=file_there) if( .not. file_there )then CALL DATE_AND_TIME(VALUES=DateTime) WRITE(dateString,'(''model1_'',i4.4,''-'',i2.2,''-'',i2.2,'':'',i2.2,'':'',i2.2,'':'',i2.2,''.hdf'')')& DateTime(1),DateTime(2),DateTime(3),DateTime(5),DateTime(6),& DateTime(7) file1in = TRIM(file1)//'.gz' INQUIRE(FILE=file1in,EXIST=file_there) if( .not. file_there )then write(errorString,'(''ERROR: File '',a,'' not found'')')& TRIM(file1) CALL Gbcs_Critical( .TRUE., & errorString,& 'NOAA_HDF_Get_AuX_OPER', & 'NOAA/NOAA_LoadAuxData.f90' ) endif command = 'cp -f '//TRIM(file1in)//' '//TRIM(dateString)//'.gz' call SYSTEM(command) command = 'gunzip -f '//TRIM(dateString)//'.gz' call SYSTEM(command) file1in = TRIM(dateString) else file1in = file1 endif INQUIRE(FILE=file2,EXIST=file_there) if( .not. file_there )then CALL DATE_AND_TIME(VALUES=DateTime) WRITE(dateString,'(''model2_'',i4.4,''-'',i2.2,''-'',i2.2,'':'',i2.2,'':'',i2.2,'':'',i2.2,''.hdf'')')& DateTime(1),DateTime(2),DateTime(3),DateTime(5),DateTime(6),& DateTime(7) file2in = TRIM(file2)//'.gz' INQUIRE(FILE=file2in,EXIST=file_there) if( .not. file_there )then write(errorString,'(''ERROR: File '',a,'' not found'')')& TRIM(file2) CALL Gbcs_Critical( .TRUE., & errorString,& 'NOAA_HDF_Get_AuX_OPER', & 'NOAA/NOAA_LoadAuxData.f90' ) endif command = 'cp -f '//TRIM(file2in)//' '//TRIM(dateString)//'.gz' call SYSTEM(command) command = 'gunzip -f '//TRIM(dateString)//'.gz' call SYSTEM(command) file2in = TRIM(dateString) else file2in = file2 endif ! Noe have files - can read then in if( onlevels .ne. AUX%Model%Grid%NLevels .or. & onelems .ne. AUX%Model%Grid%NElems .or. & onlines .ne. AUX%Model%Grid%NLines )then CALL Gbcs_Critical( .TRUE., & 'Mismatch between array bounds (nlems, nlines, nlevels )',& 'NOAA_HDF_Get_AuX_OPER', & 'NOAA/NOAA_LoadAuxData.f90' ) endif if( onlevels .ne. 18 )then CALL Gbcs_Critical( .TRUE., & 'Number of levels (total) must be = 18 (from HDF file)',& 'NOAA_HDF_Get_AuX_OPER', & 'NOAA/NOAA_LoadAuxData.f90' ) endif nelems = onelems nlines = onlines nlevels = onlevels call ncephdf_read1( file1in, file2in, nelems, nlines, nlevels, & lev_t, lev_rh, z1, surf_t, surf_rh, surf_u, surf_v, tcwv, hour, & grid, grid1, data, data1, oper_fieldlist) Plevels = oper_Plevels(:) Olevels = oper_Olevels(:) ! Forecast Model 3D Fields !---------------------------- IF ( Gbcs_Verbosity > 2 ) & WRITE(Gbcs_Log_File_Unit,*)& ' Load Forecast Model 3-D Fields (HDF file)' AUX%Model%N_3D_Model_Fields = 0 FieldID = 0 ! Pressure findx = Gbcs_Map_P_3D IF (ALLOCATED(AUX%Model%Atm(Gbcs_Map_P_3D)%d)) THEN DEALLOCATE(AUX%Model%Atm(Gbcs_Map_P_3D)%d, STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Deallocating P 3D', & 'NOAA_HDF_Get_Aux_OPER', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF ALLOCATE(AUX%Model%Atm(Gbcs_Map_P_3D)%d(nlevels+1,NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating P 3D', & 'NOAA_HDF_Get_Aux_OPER', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF ! Load pressure level data from one of the 3-D fields. DO jj = 1,nlines DO ii = 1,nelems AUX%Model%Atm(findx)%d(:,ii,jj) = Plevels END DO END DO ! FieldID = FieldID + 1 ! AUX%Model%N_3D_Model_Fields = AUX%Model%N_3D_Model_Fields + 1 ! AUX%Model%Field_3D_ID( AUX%Model%N_3D_Model_Fields ) = FieldID ! AUX%Model%Gbcs_Map_3D(AUX%Model%N_3D_Fields) = Gbcs_Map_P_3D ! AUX%Model%Field_3D_Units( AUX%Model%N_3D_Model_Fields ) = Gbcs_Units_hPa STAT = Register_3D_Field(AUX%Model, Gbcs_Map_P_3D, Gbcs_Units_hPa) ! Temperature findx = Gbcs_Map_T_3D IF (.NOT.ALLOCATED(AUX%Model%Atm(findx)%d)) THEN ALLOCATE(AUX%Model%Atm(findx)%d(nlevels+1,NElems,NLines), & STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating T 3D', & 'NOAA_HDF_Get_Aux_OPER', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF do ii=2,nlevels+1 AUX%Model%Atm(findx)%d(ii,:,:) = TRANSPOSE( REAL(lev_t(:,:,ii-1)) ) /100. end do AUX%Model%Atm(findx)%d(1,:,:) = TOA_TEMPERATURE ! FieldID = FieldID + 1 ! AUX%Model%N_3D_Model_Fields = AUX%Model%N_3D_Model_Fields + 1 ! AUX%Model%Field_3D_ID( AUX%Model%N_3D_Model_Fields ) = FieldID ! AUX%Model%Gbcs_Map_3D(AUX%Model%N_3D_Fields) = Gbcs_Map_T_3D ! AUX%Model%Field_3D_Units( AUX%Model%N_3D_Model_Fields ) = & ! Gbcs_Units_Kelvin STAT = Register_3D_Field(AUX%Model, Gbcs_Map_T_3D, Gbcs_Units_Kelvin) ! RH findx = Gbcs_Map_WV_3D IF (.NOT.ALLOCATED(AUX%Model%Atm(findx)%d)) THEN ALLOCATE(AUX%Model%Atm(findx)%d(nlevels+1,NElems,NLines), & STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating WV 3D', & 'NOAA_HDF_Get_Aux_OPER', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF do ii=2,nlevels+1 AUX%Model%Atm(findx)%d(ii,:,:) = & TRANSPOSE( REAL(lev_rh(:,:,ii-1)) ) /100. end do AUX%Model%Atm(findx)%d(1,:,:) = TOA_RH ! FieldID = FieldID + 1 ! AUX%Model%N_3D_Model_Fields = AUX%Model%N_3D_Model_Fields + 1 ! AUX%Model%Field_3D_ID( AUX%Model%N_3D_Model_Fields ) = FieldID ! AUX%Model%Gbcs_Map_3D(AUX%Model%N_3D_Fields) = Gbcs_Map_WV_3D ! AUX%Model%Field_3D_Units( AUX%Model%N_3D_Model_Fields ) = & ! Gbcs_Units_RH STAT = Register_3D_Field(AUX%Model, Gbcs_Map_WV_3D, Gbcs_Units_RH) ! Ozone findx = Gbcs_Map_O3_3D IF (.NOT.ALLOCATED(AUX%Model%Atm(findx)%d)) THEN ALLOCATE(AUX%Model%Atm(findx)%d(nlevels+1,NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating O3 3D', & 'NOAA_HDF_Get_Aux_OPER', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF ! Need everything above to be setup so now get O3 CALL Get_Ozone_climatology( IMG%Mean_Time%Month, AUX ) ! DO ii = 1,nlevels+1 ! ! Get from above values ! AUX%Model%Atm(findx)%d(ii,:,:) = OLevels(ii) ! END DO ! FieldID = FieldID + 1 ! AUX%Model%N_3D_Model_Fields = AUX%Model%N_3D_Model_Fields + 1 ! AUX%Model%Field_3D_ID( AUX%Model%N_3D_Model_Fields ) = FieldID ! AUX%Model%Gbcs_Map_3D(AUX%Model%N_3D_Fields) = Gbcs_Map_O3_3D ! AUX%Model%Field_3D_Units( AUX%Model%N_3D_Model_Fields ) = Gbcs_Units_ppmv STAT = Register_3D_Field(AUX%Model, Gbcs_Map_O3_3D, Gbcs_Units_ppmv) ! Set NLevels + 1 due to TOA values AUX%Model%Grid%NLevels = nlevels+1 ! Forecast Model 2D Fields !---------------------------- IF ( Gbcs_Verbosity > 2 ) & WRITE(Gbcs_Log_File_Unit,*)' Load Forecast Model 2-D Fields (HDF file)' AUX%Model%N_2D_Model_Fields = 0 FieldID = 0 ! 1000mb height findx = Gbcs_Map_z1000 IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating z1000', & 'NOAA_HDF_Get_Aux_OPER', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF AUX%Model%Srf(findx)%d = TRANSPOSE( REAL(z1) ) /100. ! FieldID = FieldID + 1 ! AUX%Model%N_2D_Model_Fields = AUX%Model%N_2D_Model_Fields + 1 ! AUX%Model%Srf_Index(findx) = AUX%Model%N_2D_Model_Fields ! AUX%Model%Field_2D_ID( AUX%Model%N_2D_Model_Fields ) = FieldID ! AUX%Model%Field_2D_Map( AUX%Model%N_2D_Model_Fields ) = Gbcs_Map_z1000 ! AUX%Model%Field_2D_Units( AUX%Model%N_2D_Model_Fields ) = Gbcs_Units_m STAT = Register_2D_Field(AUX%Model, Gbcs_Map_z1000, Gbcs_Units_m) ! Orography findx = Gbcs_Map_OROG IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating OROG', & 'NOAA_HDF_Get_Aux_OPER', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF AUX%Model%Srf(findx)%d = 0.0 ! FieldID = FieldID + 1 ! AUX%Model%N_2D_Model_Fields = AUX%Model%N_2D_Model_Fields + 1 ! AUX%Model%Srf_Index(findx) = AUX%Model%N_2D_Model_Fields ! AUX%Model%Field_2D_ID( AUX%Model%N_2D_Model_Fields ) = FieldID ! AUX%Model%Field_2D_Map( AUX%Model%N_2D_Model_Fields ) = Gbcs_Map_OROG ! AUX%Model%Field_2D_Units( AUX%Model%N_2D_Model_Fields ) = Gbcs_Units_m STAT = Register_2D_Field(AUX%Model, Gbcs_Map_OROG, Gbcs_Units_m) ! Surface temperature ! Check difference between SST/Surface etc findx = Gbcs_Map_T_SKIN IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating T Skin', & 'NOAA_HDF_Get_Aux_OPER', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF AUX%Model%Srf(findx)%d = TRANSPOSE( REAL(surf_t) ) /100. ! FieldID = FieldID + 1 ! AUX%Model%N_2D_Model_Fields = AUX%Model%N_2D_Model_Fields + 1 ! AUX%Model%Srf_Index(findx) = AUX%Model%N_2D_Model_Fields ! AUX%Model%Field_2D_ID( AUX%Model%N_2D_Model_Fields ) = FieldID ! AUX%Model%Field_2D_Map( AUX%Model%N_2D_Model_Fields ) = Gbcs_Map_T_SKIN ! AUX%Model%Field_2D_Units( AUX%Model%N_2D_Model_Fields ) = Gbcs_Units_Kelvin STAT = Register_2D_Field(AUX%Model, Gbcs_Map_T_SKIN, Gbcs_Units_Kelvin) ! Total Column Water Vapour findx = Gbcs_Map_TCWV IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating TCWV', & 'NOAA_HDF_Get_Aux_OPER', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF AUX%Model%Srf(findx)%d = TRANSPOSE( REAL(tcwv) ) /100. FieldID = FieldID + 1 ! AUX%Model%N_2D_Model_Fields = AUX%Model%N_2D_Model_Fields + 1 ! AUX%Model%Srf_Index(findx) = AUX%Model%N_2D_Model_Fields ! AUX%Model%Field_2D_ID( AUX%Model%N_2D_Model_Fields ) = FieldID ! AUX%Model%Field_2D_Map( AUX%Model%N_2D_Model_Fields ) = Gbcs_Map_TCWV ! AUX%Model%Field_2D_Units( AUX%Model%N_2D_Model_Fields ) = Gbcs_Units_Kg_per_m2 STAT = Register_2D_Field(AUX%Model, Gbcs_Map_TCWV, Gbcs_Units_Kg_Per_m2) ! Wind - U findx = Gbcs_Map_U_SURF IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating U Surf', & 'NOAA_HDF_Get_Aux_OPER', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF AUX%Model%Srf(findx)%d = TRANSPOSE( REAL(surf_u) ) /100. ! FieldID = FieldID + 1 ! AUX%Model%N_2D_Model_Fields = AUX%Model%N_2D_Model_Fields + 1 ! AUX%Model%Srf_Index(findx) = AUX%Model%N_2D_Model_Fields ! AUX%Model%Field_2D_ID( AUX%Model%N_2D_Model_Fields ) = FieldID ! AUX%Model%Field_2D_Map( AUX%Model%N_2D_Model_Fields ) = Gbcs_Map_U_SURF ! AUX%Model%Field_2D_Units( AUX%Model%N_2D_Model_Fields ) = Gbcs_Units_m_per_s STAT = Register_2D_Field(AUX%Model, Gbcs_Map_U_Surf, Gbcs_Units_m_per_s) ! Wind - V findx = Gbcs_Map_V_SURF IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating V Surf', & 'NOAA_HDF_Get_Aux_OPER', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF AUX%Model%Srf(findx)%d = TRANSPOSE( REAL(surf_v) ) /100. ! FieldID = FieldID + 1 ! AUX%Model%N_2D_Model_Fields = AUX%Model%N_2D_Model_Fields + 1 ! AUX%Model%Srf_Index(findx) = AUX%Model%N_2D_Model_Fields ! AUX%Model%Field_2D_ID( AUX%Model%N_2D_Model_Fields ) = FieldID ! AUX%Model%Field_2D_Map( AUX%Model%N_2D_Model_Fields ) = Gbcs_Map_V_SURF ! AUX%Model%Field_2D_Units( AUX%Model%N_2D_Model_Fields ) = Gbcs_Units_m_per_s STAT = Register_2D_Field(AUX%Model, Gbcs_Map_V_Surf, Gbcs_Units_m_per_s) ! Wind - 10m findx = Gbcs_Map_Wind_10M IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating 10m Wind', & 'NOAA_HDF_Get_Aux_OPER', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF do jj=1,nlines do ii=1,nelems if( 25 .lt. AUX%Model%Srf(Gbcs_Map_U_SURF)%d(ii,jj) .or. & 25 .lt. AUX%Model%Srf(Gbcs_Map_V_SURF)%d(ii,jj) )then if( 25 .lt. AUX%Model%Srf(Gbcs_Map_U_SURF)%d(ii,jj) .and. & 25 .lt. AUX%Model%Srf(Gbcs_Map_V_SURF)%d(ii,jj) )then AUX%Model%Srf(findx)%d(ii,jj) = 8 else if( 25 .lt. AUX%Model%Srf(Gbcs_Map_U_SURF)%d(ii,jj) )then AUX%Model%Srf(findx)%d(ii,jj) = & ABS(AUX%Model%Srf(Gbcs_Map_V_SURF)%d(ii,jj)) else if( 25 .lt. AUX%Model%Srf(Gbcs_Map_V_SURF)%d(ii,jj) )then AUX%Model%Srf(findx)%d(ii,jj) = & ABS(AUX%Model%Srf(Gbcs_Map_U_SURF)%d(ii,jj)) endif else AUX%Model%Srf(findx)%d(ii,jj) = sqrt(& AUX%Model%Srf(Gbcs_Map_V_SURF)%d(ii,jj)**2 + & AUX%Model%Srf(Gbcs_Map_U_SURF)%d(ii,jj)**2) endif end do end do ! FieldID = FieldID + 1 ! AUX%Model%N_2D_Model_Fields = AUX%Model%N_2D_Model_Fields + 1 ! AUX%Model%Srf_Index(findx) = AUX%Model%N_2D_Model_Fields ! AUX%Model%Field_2D_ID( AUX%Model%N_2D_Model_Fields ) = FieldID ! AUX%Model%Field_2D_Map( AUX%Model%N_2D_Model_Fields ) = Gbcs_Map_Wind_10M ! AUX%Model%Field_2D_Units( AUX%Model%N_2D_Model_Fields ) = Gbcs_Units_m_per_s STAT = Register_2D_Field(AUX%Model, Gbcs_Map_Wind_10M, Gbcs_Units_m_per_s) ! Get surface pressure findx = Gbcs_Map_P_Surf IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating P Surf', & 'NOAA_HDF_Get_Aux_OPER', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF ! Calculate surface pressure call surfpressure_GOES(nlines, nelems, & AUX%Model%Srf(Gbcs_Map_z1000)%d, & AUX%Model%Atm(Gbcs_Map_T_3D)%d(nlevels,:,:), & AUX%Model%Atm(Gbcs_Map_P_3D)%d(nlevels,:,:), & AUX%Model%Srf(Gbcs_Map_P_Surf)%d ) ! FieldID = FieldID + 1 ! AUX%Model%N_2D_Model_Fields = AUX%Model%N_2D_Model_Fields + 1 ! AUX%Model%Srf_Index(findx) = AUX%Model%N_2D_Model_Fields ! AUX%Model%Field_2D_ID( AUX%Model%N_2D_Model_Fields ) = FieldID ! AUX%Model%Field_2D_Map( AUX%Model%N_2D_Model_Fields ) = Gbcs_Map_P_Surf ! AUX%Model%Field_2D_Units( AUX%Model%N_2D_Model_Fields ) = Gbcs_Units_hPa STAT = Register_2D_Field(AUX%Model, Gbcs_Map_P_Surf, Gbcs_Units_hPa) ! Total surface RH findx = Gbcs_Map_WV_2M IF (.NOT.ALLOCATED(AUX%Model%Srf(findx)%d)) THEN ALLOCATE(AUX%Model%Srf(findx)%d(NElems,NLines), STAT=STAT) IF( 0 .ne. STAT )THEN CALL Gbcs_Critical( .TRUE., & 'Allocating 2m WV', & 'NOAA_HDF_Get_Aux_OPER', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END IF AUX%Model%Srf(findx)%d = TRANSPOSE( REAL(surf_rh) ) /100. ! FieldID = FieldID + 1 ! AUX%Model%N_2D_Model_Fields = AUX%Model%N_2D_Model_Fields + 1 ! AUX%Model%Srf_Index(findx) = AUX%Model%N_2D_Model_Fields ! AUX%Model%Field_2D_ID( AUX%Model%N_2D_Model_Fields ) = FieldID ! AUX%Model%Field_2D_Map( AUX%Model%N_2D_Model_Fields ) = Gbcs_Map_WV_2M ! AUX%Model%Field_2D_Units( AUX%Model%N_2D_Model_Fields ) = Gbcs_Units_RH STAT = Register_2D_Field(AUX%Model, Gbcs_Map_WV_2M, Gbcs_Units_RH) RETURN END SUBROUTINE NOAA_HDF_Get_Aux_OPER ! Code taken from Wen Meng (NOAA) SUBROUTINE ncephdf_read1(file1, file2, nelems, nlines, nlevels, & lev_t, lev_rh, z1, surf_t, surf_rh, surf_u, surf_v, tcwv, hour,& grid, grid1, data, data1, fieldlist) ! Function to read in NCEP data from binary files ! Requires the files... ! ---------------------- Declare Variables ---------------------------- ! INPUT CHARACTER(LEN=*), INTENT(IN) :: file1 CHARACTER(LEN=*), INTENT(IN) :: file2 INTEGER, INTENT(IN) :: nelems INTEGER, INTENT(IN) :: nlines INTEGER, INTENT(IN) :: nlevels ! OUTPUT integer, intent(out) :: lev_t(nlines,nelems,nlevels) ! Profiles of temperature on levels (K) integer, intent(out) :: lev_rh(nlines,nelems,nlevels) ! Profiles of relative humidity on levels (%) integer, intent(out) :: z1(nlines,nelems) ! Altitude (m) at 1000.0 mB integer, intent(out) :: surf_t(nlines,nelems) ! Surface temperature (K) integer, intent(out) :: surf_rh(nlines,nelems) ! Surface relative humidity (%) integer, intent(out) :: surf_u(nlines,nelems) ! Surface windspeed, u (m / s) integer, intent(out) :: surf_v(nlines,nelems) ! Surface windspeed, v (m / s) integer, intent(out) :: tcwv(nlines,nelems) ! TCWV real, intent(in) :: hour ! hour for process integer, intent(inout) :: grid(nlines,nelems) ! Temporary array for storing grids integer, intent(inout) :: grid1(nlines,nelems) ! Temporary array for storing grids real, intent(inout) :: data(nelems,nlines) real, intent(inout) :: data1(nelems,nlines) character(len=20), intent(in) :: fieldlist(42) integer k ! All outputs are on a 1x1 degree grid ! LOCAL VARIABLES integer grid_header(64) ! GRID file header integer status ! Status flag for IGGET integer nrows, ncols ! Number of rows and columns in GRID file integer ii, i,j ! Counter ! Variables about HDF files integer sfstart, sfn2index, sfselect integer sfrdata, sfendacc, sfend integer sd_id, sds_id, sds_index integer sd_id1, sds_id1, sds_index1 integer start(2), stride(2),edge(2) integer ed_end_status, sds_end_status, ret character(len=20) sds_name INTEGER, parameter :: DFACC_RDONLY = 1 ! ------------------------- Processing -------------------------------- ! Read in atmospheric profile data for temperature and relative humidity ! from the file GRID1000 ! print*, "open hour is ", hour sd_id=sfstart(file1, DFACC_RDONLY) call checkHDF(sd_id,"sdstart") start(1)=0 start(2)=0 stride(1)=1 stride(2)=1 edge(1)=360 edge(2)=181 ! read in temperature grids do ii = 1, nlevels sds_name=fieldlist(ii) sds_index=sfn2index(sd_id,sds_name) call checkHDF(sds_index,"sfn2index") sds_id = sfselect(sd_id,sds_index) call checkHDF(sds_id,"sfselect") ret=sfrdata(sds_id, start, stride, edge, data) call checkHDF(ret,"sfrdata") sds_end_status=sfendacc(sds_id) call checkHDF(ret,"sfendacc") do i = 1, nlines do j = 1, nelems grid(i,j)=Nint(data(j,nlines+1-i)*(3.-hour)/3.*100) enddo enddo lev_t(:,:,ii) = grid enddo ! read in relative humidity grids do ii = nlevels+1, 2*nlevels sds_name=fieldlist(ii) sds_index=sfn2index(sd_id,sds_name) call checkHDF(sds_index,"sfn2index") sds_id = sfselect(sd_id,sds_index) call checkHDF(sds_id,"sfselect") ret=sfrdata(sds_id, start, stride, edge, data) call checkHDF(ret,"sfrdata") sds_end_status=sfendacc(sds_id) call checkHDF(ret,"sfendacc") do i = 1, nlines do j = 1, nelems grid(i,j)=Nint(data(j,nlines+1-i)*(3.-hour)/3.*100) enddo enddo lev_rh(:,:,ii-nlevels) = grid enddo ! Read in the 1000mB altitude, z1 sds_name=fieldlist(2*nlevels+1) sds_index=sfn2index(sd_id,sds_name) call checkHDF(sds_index,"sfn2index") sds_id = sfselect(sd_id,sds_index) call checkHDF(sds_id,"sfselect") ret=sfrdata(sds_id, start, stride, edge, data) call checkHDF(ret,"sfrdata") sds_end_status=sfendacc(sds_id) call checkHDF(ret,"sfendacc") do i = 1, nlines do j = 1, nelems z1(i,j)=Nint(data(j,nlines+1-i)*(3.-hour)/3.*100) enddo enddo ! read in total column water vapour, tcwv sds_name=fieldlist(2*nlevels+2) sds_index=sfn2index(sd_id,sds_name) call checkHDF(sds_index,"sfn2index") sds_id = sfselect(sd_id,sds_index) call checkHDF(sds_id,"sfselect") ret=sfrdata(sds_id, start, stride, edge, data) call checkHDF(ret,"sfrdata") sds_end_status=sfendacc(sds_id) call checkHDF(ret,"sfendacc") do i = 1, nlines do j = 1, nelems tcwv(i,j)=Nint(data(j,nlines+1-i)*(3.-hour)/3.*100) enddo enddo ! read in surface temperature, surf_t sds_name=fieldlist(2*nlevels+3) sds_index=sfn2index(sd_id,sds_name) call checkHDF(sds_index,"sfn2index") sds_id = sfselect(sd_id,sds_index) call checkHDF(sds_id,"sfselect") ret=sfrdata(sds_id, start, stride, edge, data) call checkHDF(ret,"sfrdata") sds_end_status=sfendacc(sds_id) call checkHDF(ret,"sfendacc") do i = 1, nlines do j = 1, nelems surf_t(i,j)=Nint(data(j,nlines+1-i)*(3.-hour)/3.*100) enddo enddo ! read in surface RH, surf_rh sds_name=fieldlist(2*nlevels+4) sds_index=sfn2index(sd_id,sds_name) call checkHDF(sds_index,"sfn2index") sds_id = sfselect(sd_id,sds_index) call checkHDF(sds_id,"sfselect") ret=sfrdata(sds_id, start, stride, edge, data) call checkHDF(ret,"sfrdata") sds_end_status=sfendacc(sds_id) call checkHDF(ret,"sfendacc") do i = 1, nlines do j = 1, nelems surf_rh(i,j)=Nint(data(j,nlines+1-i)*(3.-hour)/3.*100) enddo enddo ! read in surface windspeed, surf_u sds_name=fieldlist(2*nlevels+5) sds_index=sfn2index(sd_id,sds_name) call checkHDF(sds_index,"sfn2index") sds_id = sfselect(sd_id,sds_index) call checkHDF(sds_id,"sfselect") ret=sfrdata(sds_id, start, stride, edge, data) call checkHDF(ret,"sfrdata") sds_end_status=sfendacc(sds_id) call checkHDF(ret,"sfendacc") do i = 1, nlines do j = 1, nelems surf_u(i,j)=Nint(data(j,nlines+1-i)*(3.-hour)/3.*100) enddo enddo ! read in surface windspeed, surf_v sds_name=fieldlist(2*nlevels+6) sds_index=sfn2index(sd_id,sds_name) call checkHDF(sds_index,"sfn2index") sds_id = sfselect(sd_id,sds_index) call checkHDF(sds_id,"sfselect") ret=sfrdata(sds_id, start, stride, edge, data) call checkHDF(ret,"sfrdata") sds_end_status=sfendacc(sds_id) call checkHDF(ret,"sfendacc") do i = 1, nlines do j = 1, nelems surf_v(i,j)=Nint(data(j,nlines+1-i)*(3.-hour)/3.*100) enddo enddo ed_end_status=sfend(sd_id) sd_id=sfstart(file2, DFACC_RDONLY) do ii = 1, nlevels sds_name=fieldlist(ii) sds_index=sfn2index(sd_id,sds_name) call checkHDF(sds_index,"sfn2index") sds_id = sfselect(sd_id,sds_index) call checkHDF(sds_id,"sfselect") ret=sfrdata(sds_id, start, stride, edge, data) call checkHDF(ret,"sfrdata") sds_end_status=sfendacc(sds_id) call checkHDF(ret,"sfendacc") do i = 1, nlines do j = 1, nelems grid(i,j)=Nint(data(j,nlines+1-i)*hour/3.*100) enddo enddo lev_t(:,:,ii) = grid+lev_t(:,:,ii) enddo ! read in relative humidity grids do ii = nlevels+1, 2*nlevels sds_name=fieldlist(ii) sds_index=sfn2index(sd_id,sds_name) call checkHDF(sds_index,"sfn2index") sds_id = sfselect(sd_id,sds_index) call checkHDF(sds_id,"sfselect") ret=sfrdata(sds_id, start, stride, edge, data) call checkHDF(ret,"sfrdata") sds_end_status=sfendacc(sds_id) call checkHDF(ret,"sfendacc") do i = 1, nlines do j = 1, nelems grid(i,j)=Nint(data(j,nlines+1-i)*hour/3.*100) enddo enddo lev_rh(:,:,ii-nlevels) = grid+lev_rh(:,:,ii-nlevels) enddo ! Read in the 1000mB altitude, z1 sds_name=fieldlist(2*nlevels+1) sds_index=sfn2index(sd_id,sds_name) call checkHDF(sds_index,"sfn2index") sds_id = sfselect(sd_id,sds_index) call checkHDF(sds_id,"sfselect") ret=sfrdata(sds_id, start, stride, edge, data) call checkHDF(ret,"sfrdata") sds_end_status=sfendacc(sds_id) call checkHDF(ret,"sfendacc") do i = 1, nlines do j = 1, nelems z1(i,j)=Nint(data(j,nlines+1-i)*hour/3.*100)+z1(i,j) enddo enddo ! read in total column water vapour, tcwv sds_name=fieldlist(2*nlevels+2) sds_index=sfn2index(sd_id,sds_name) call checkHDF(sds_index,"sfn2index") sds_id = sfselect(sd_id,sds_index) call checkHDF(sds_id,"sfselect") ret=sfrdata(sds_id, start, stride, edge, data) call checkHDF(ret,"sfrdata") sds_end_status=sfendacc(sds_id) call checkHDF(ret,"sfendacc") do i = 1, nlines do j = 1, nelems tcwv(i,j)=Nint(data(j,nlines+1-i)*hour/3.*100)+tcwv(i,j) enddo enddo ! read in surface temperature, surf_t sds_name=fieldlist(2*nlevels+3) sds_index=sfn2index(sd_id,sds_name) call checkHDF(sds_index,"sfn2index") sds_id = sfselect(sd_id,sds_index) call checkHDF(sds_id,"sfselect") ret=sfrdata(sds_id, start, stride, edge, data) call checkHDF(ret,"sfrdata") sds_end_status=sfendacc(sds_id) call checkHDF(ret,"sfendacc") do i = 1, nlines do j = 1, nelems surf_t(i,j)=Nint(data(j,nlines+1-i)*hour/3.*100)+surf_t(i,j) enddo enddo ! read in surface RH, surf_rh sds_name=fieldlist(2*nlevels+4) sds_index=sfn2index(sd_id,sds_name) call checkHDF(sds_index,"sfn2index") sds_id = sfselect(sd_id,sds_index) call checkHDF(sds_id,"sfselect") ret=sfrdata(sds_id, start, stride, edge, data) call checkHDF(ret,"sfrdata") sds_end_status=sfendacc(sds_id) call checkHDF(ret,"sfendacc") do i = 1, nlines do j = 1, nelems surf_rh(i,j)=Nint(data(j,nlines+1-i)*hour/3.*100)+surf_rh(i,j) enddo enddo ! read in surface windspeed, surf_u sds_name=fieldlist(2*nlevels+5) sds_index=sfn2index(sd_id,sds_name) call checkHDF(sds_index,"sfn2index") sds_id = sfselect(sd_id,sds_index) call checkHDF(sds_id,"sfselect") ret=sfrdata(sds_id, start, stride, edge, data) call checkHDF(ret,"sfrdata") sds_end_status=sfendacc(sds_id) call checkHDF(ret,"sfendacc") do i = 1, nlines do j = 1, nelems surf_u(i,j)=Nint(data(j,nlines+1-i)*hour/3.*100)+surf_u(i,j) enddo enddo ! read in surface windspeed, surf_v sds_name=fieldlist(2*nlevels+6) sds_index=sfn2index(sd_id,sds_name) call checkHDF(sds_index,"sfn2index") sds_id = sfselect(sd_id,sds_index) call checkHDF(sds_id,"sfselect") ret=sfrdata(sds_id, start, stride, edge, data) call checkHDF(ret,"sfrdata") sds_end_status=sfendacc(sds_id) call checkHDF(ret,"sfendacc") do i = 1, nlines do j = 1, nelems surf_v(i,j)=Nint(data(j,nlines+1-i)*hour/3.*100)+surf_v(i,j) enddo enddo ed_end_status=sfend(sd_id) call checkHDF(ret,"sfend") END SUBROUTINE ncephdf_read1 SUBROUTINE checkHDF( status, program_string ) INTEGER, INTENT(IN) :: status CHARACTER(LEN=*), INTENT(IN) :: program_string INTEGER :: ret integer :: heprnt CHARACTER(LEN=256) :: errorString IF( -1 .EQ. status )THEN ! We have an error write(errorString,'(''HDF Error Calling '',a)')TRIM(program_string) CALL Gbcs_Critical( .TRUE., & TRIM(errorString),& 'checkHDF', & 'NOAA/NOAA_LoadAuxData.f90' ) ENDIF END SUBROUTINE checkHDF #endif ! Load in data from a grib file - use wgrib/wgrib2 to decode ! Reads in two model GFS files to make up the 47 level data SUBROUTINE NOAA_GRIB_GFS_Load_Aux_Data( SAT, IMG, AUX, RTM) USE GbcsAtmPhysics USE GbcsPDFLoaders USE GbcsCovMatLoaders IMPLICIT NONE TYPE(Satellite), INTENT(INOUT) :: SAT TYPE(Imagery), INTENT(INOUT) :: IMG TYPE(Aux_Data), INTENT(INOUT) :: AUX TYPE(RTM_Interface), INTENT(INOUT) :: RTM INTEGER :: start_year INTEGER :: start_month INTEGER :: start_day INTEGER :: start_hour_anal INTEGER :: start_hour_forecast INTEGER :: end_year INTEGER :: end_month INTEGER :: end_day INTEGER :: end_hour_anal INTEGER :: end_hour_forecast CHARACTER(LEN=512) :: gfs_hdf_archive CHARACTER(LEN=512) :: file1 CHARACTER(LEN=512) :: file2 REAL, ALLOCATABLE :: landseamask(:,:) INTEGER :: day_num REAL :: hour INTEGER :: I, J ! Get SZA/Land mask for model - must have access to MCIDAS data files if( IMG%DataFrmt .eq. "MCIDAS" )then CALL MCIDAS_Get_Ancilliary_Data( AUX, IMG ) else CALL NotMCIDAS_Get_Ancilliary_Data( AUX, IMG, SAT ) endif ! Get name of input file ! Get surrounding model files ! Start start_year = IMG%Mean_Time%Year start_month = IMG%Mean_Time%Month start_day = IMG%Mean_Time%Day start_hour_anal = (IMG%Mean_Time%Hour/6)*6 start_hour_forecast = ((IMG%Mean_Time%Hour - start_hour_anal)/3)*3 ! Next file in line if( 3 .eq. start_hour_forecast )then end_hour_anal = start_hour_anal+6 end_hour_forecast = 0 else end_hour_anal = start_hour_anal end_hour_forecast = 3 endif if( end_hour_anal .lt. 24 )then end_day = start_day end_month = start_month end_year = start_year else end_hour_anal = 0 end_hour_forecast = 0 call date_to_day_number_mcidas(start_year,start_month,start_day,day_num) day_num = day_num+1 if( (0 .eq. MOD(start_year,4) .and. 0 .ne. MOD(start_year,100)) .or. & (0 .eq. MOD(start_year,400)) )then if( day_num .gt. 366 )then end_year = end_year+1 end_month = 1 end_day = 1 else end_year = start_year call day_number_to_date_mcidas( end_year, day_num, end_month, end_day ) endif else if( day_num .gt. 365 )then end_year = end_year+1 end_month = 1 end_day = 1 else end_year = start_year call day_number_to_date_mcidas( end_year, day_num, end_month, end_day ) endif endif endif ! Get floating fraction of 3 hour window hour = IMG%Mean_Time%Hour + IMG%Mean_Time%Minute/60. + & IMG%Mean_Time%Seconds/3600. hour = hour - start_hour_anal - start_hour_forecast ! Now have dates and times of bracketing model files ! make filenames write(file1,& '(a,''/'',i4.4,''/'',i2.2,''/'',i4.4,i2.2,i2.2,''gdas1.t'',i2.2,''z.pgrbf'',i2.2,''.hdf'')')& TRIM(AUX%Model%DataDir),start_year,start_month,start_year,start_month,& start_day,start_hour_anal,start_hour_forecast write(file2,& '(a,''/'',i4.4,''/'',i2.2,''/'',i4.4,i2.2,i2.2,''gdas1.t'',i2.2,''z.pgrbf'',i2.2,''.hdf'')')& TRIM(AUX%Model%DataDir),end_year,end_month,end_year,end_month,& end_day,end_hour_anal,end_hour_forecast write(*,'('' Interpolating between : '')') write(*,'('' '',a)')TRIM(file1) write(*,'('' '',a)')TRIM(file2) ! CALL NOAA_GRIB_Get_Aux( file1, file2, hour, AUX, IMG, LandSeaMask ) END SUBROUTINE NOAA_GRIB_GFS_Load_Aux_Data ! ! Routine to find relevant GFS data (forecast or analysis) ! which is closest in time to requred times (start/end) ! LOGICAL FUNCTION Find_GFS_Data( AUX, start_year,start_month,start_day,& start_hours,start_hour_forecast,end_year,end_month,& end_day,end_hours,end_hour_forecast,file1,file2 ) TYPE(Aux_Data), INTENT(IN) :: AUX INTEGER, INTENT(IN) :: start_year INTEGER, INTENT(IN) :: start_month INTEGER, INTENT(IN) :: start_day INTEGER, INTENT(IN) :: start_hours INTEGER, INTENT(IN) :: start_hour_forecast INTEGER, INTENT(IN) :: end_year INTEGER, INTENT(IN) :: end_month INTEGER, INTENT(IN) :: end_day INTEGER, INTENT(IN) :: end_hours INTEGER, INTENT(IN) :: end_hour_forecast CHARACTER(LEN=*), INTENT(OUT) :: file1 CHARACTER(LEN=*), INTENT(OUT) :: file2 ! Local variables INTEGER :: time, forecast CHARACTER(LEN=256) :: Filename LOGICAL :: Found_File ! Start time time = start_hours forecast = start_hour_forecast Found_file = .FALSE. DO WHILE( .not. Found_File ) ! Make filename WRITE(Filename,'(a,''/gfs.t'',i2.2,''z.pgrbf'',i2.2,''.grib2'')')& TRIM(AUX%Model%DataDir),time,forecast INQUIRE(FILE=Filename,EXIST=Found_File) IF( .not.Found_File )THEN ! Increment times forecast = forecast+6 time = time-6 ENDIF END DO file1 = Filename END FUNCTION Find_GFS_Data ! Routine to get executable directory so we can find where wgrib.exe and ! wgrib2.exe sit SUBROUTINE Get_Executable_Directory( exe_dir ) CHARACTER(LEN=*), INTENT(OUT) :: exe_dir ! Local variable INTEGER :: I INTEGER :: ios INTEGER :: pid INTEGER :: lun INTEGER :: start INTEGER :: end CHARACTER(LEN=256) :: line CHARACTER(LEN=256) :: command ! Works on Linux - use /proc/process_id/exe link to work out where executable ! is pid = GETPID() IF( 10 .gt. pid )THEN WRITE(command,'(''/bin/ls -l /proc/'',i1,''/exe > /tmp/exe_dir.link'')') ELSE IF( 10 .le. pid .and. 100 .gt. pid )THEN WRITE(command,'(''/bin/ls -l /proc/'',i2,''/exe > /tmp/exe_dir.link'')') ELSE IF( 100 .le. pid .and. 1000 .gt. pid )THEN WRITE(command,'(''/bin/ls -l /proc/'',i3,''/exe > /tmp/exe_dir.link'')') ELSE IF( 1000 .le. pid .and. 10000 .gt. pid )THEN WRITE(command,'(''/bin/ls -l /proc/'',i4,''/exe > /tmp/exe_dir.link'')') ELSE IF( 10000 .le. pid .and. 100000 .gt. pid )THEN WRITE(command,'(''/bin/ls -l /proc/'',i5,''/exe > /tmp/exe_dir.link'')') ELSE CALL Gbcs_Critical(.TRUE.,'PID out of range','Get_Executable_Direcory',& 'NOAA_LoadAuxData.f90') ENDIF CALL SYSTEM(command) lun = Get_New_File_Unit() OPEN(UNIT=lun,FILE='/tmp/exe_dir.link',STATUS='OLD',IOSTAT=ios) IF( 0 .ne. ios )THEN CALL Gbcs_Critical(.TRUE.,'Cannot open /tmp/exe_dir.link',& 'Get_Executable_Direcory','NOAA_LoadAuxData.f90') ENDIF READ(lun,'(a)',IOSTAT=ios)line IF( 0 .ne. ios )THEN CALL Gbcs_Critical(.TRUE.,'Cannot read /tmp/exe_dir.link',& 'Get_Executable_Direcory','NOAA_LoadAuxData.f90') ENDIF CLOSE(lun) ! Find extent of link name start = -1 Loop: DO I=LEN_TRIM(line),1,-1 IF( ' ' .eq. line(I:I) )THEN start = I+1 EXIT Loop ENDIF END DO Loop IF( -1 .ne. start )THEN CALL Gbcs_Critical(.TRUE.,'Cannot find start of executable string',& 'Get_Executable_Direcory','NOAA_LoadAuxData.f90') ENDIF end = -1 Loop2: DO I=LEN_TRIM(line),1,-1 IF( '/' .eq. line(I:I) )THEN end = I EXIT Loop2 ENDIF END DO Loop2 IF( -1 .ne. end )THEN CALL Gbcs_Critical(.TRUE.,'Cannot find end of executable dir',& 'Get_Executable_Direcory','NOAA_LoadAuxData.f90') ENDIF exe_dir = line(start:end) END SUBROUTINE Get_Executable_Directory END MODULE NOAA_LoadAuxData