module m_interface_crtm !$$$ module documentation block ! . . . . ! module: interface_crtm ! prgmmr: R.Yang date: 2008-03 ! ! abstract: compute brightness temperature using CRTM with profiles of Nature Run ! at the real satellite observing systems. Observing systems include ! hirs2_n14,hirs3_n16,hirs3_n17,amsua_n15-n18,amsub_n15-n17, Aqua_airs, ! and Aqua_amsua. ! ! program history log: ! 2008-03-15 Yang - initiate the code ! 2008-09-16 Errico - modify the setting of surface type (for MW computation) ! 2008-09-17 Yang - modify to include msu use crtm_module implicit none private public :: setupf_names public :: InitGetbt_crtm integer, save :: nf_zsfc, nf_u10m, nf_v10m integer, save :: nf_ismk, nf_almk integer, save :: nf_temp, nf_sphu, nf_ozon contains ! X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X ! subroutine setupf_names (nfields2d,nfields3d,nfields,f_names) implicit none integer :: nfields2d, nfields3d,nfields character(len=4) :: f_names(nfields) call find_name (nfields2d,f_names(1:nfields2d),'zsfc', nf_zsfc,.true.) call find_name (nfields2d,f_names(1:nfields2d),'u10m', nf_u10m,.true.) call find_name (nfields2d,f_names(1:nfields2d),'v10m', nf_v10m,.true.) call find_name (nfields2d,f_names(1:nfields2d),'ismk', nf_ismk,.false.) call find_name (nfields2d,f_names(1:nfields2d),'almk', nf_almk,.false.) call find_name (nfields3d,f_names(1+nfields2d:nfields),'temp', & nf_temp,.true.) call find_name (nfields3d,f_names(1+nfields2d:nfields),'sphu', & nf_sphu,.true.) call find_name (nfields3d,f_names(1+nfields2d:nfields),'ozon', & nf_ozon,.true.) ! end subroutine setupf_names ! X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X ! subroutine InitGetbt_crtm (prof_ndim4,nlevs,nfields2d,nfields3d, & prof_info,prof_plevs,prof_2d,prof_3d,lcldtop,pcldtop, & linit,d_type,subtype,btchan,nchan,channelinfo,soza, & sfc_type) use m_kinds ! ------------ ! Module usage ! ------------ USE CRTM_SpcCoeff USE SpcCoeff_Define ! -- Utility modules USE Type_Kinds USE Message_Handler ! -- CRTM module USE CRTM_Module USE CRTM_Atmosphere_Define ! -- Modules to read in Atmosphere and Surface data USE CRTM_Atmosphere_Binary_IO USE CRTM_Surface_Binary_IO USE CRTM_Surface_Define USE CRTM_SensorInfo USE CRTM_ChannelInfo USE CRTM_ChannelInfo_Define USE CRTM_surface_ir_emissivity implicit none integer,intent(in) :: prof_ndim4, nlevs, nfields2d, nfields3d integer,intent(in) :: lcldtop real(rkind1),intent(in) :: pcldtop real(rkind1),intent(in) :: prof_info(prof_ndim4) real(rkind1),intent(in) :: prof_2d(nfields2d) real(rkind1),intent(in) :: prof_3d(nlevs,nfields3d) real(rkind1),intent(in) :: prof_plevs(nlevs) integer,intent(in) :: nchan real(rkind3),intent(inout) :: btchan(nchan) real(rkind1) :: soza ! solar zenith angle TYPE(CRTM_ChannelInfo_type ),intent(inout) :: ChannelInfo character(len=*),intent (in) :: d_type logical, intent (inout) :: linit integer,intent(out) :: sfc_type ! flag (1=ocean, 2=ice, 3=ocean) integer,intent(in) :: subtype ! ID # for either satellite or instrument type ! local variable real(rkind1) :: lza,panglr,rato real(rkind1) :: step,start,pi_fac,deg2rad,rad2deg CHARACTER( * ), PARAMETER :: PROGRAM_NAME = 'Test_setup_crtm' character(len=*), parameter :: routine_name='interface_crtm' ! --------- ! Variables ! --------- INTEGER :: Error_Status, L, i, j, k, cloud_layer_idx, i_prof INTEGER :: Allocate_Status CHARACTER( 256 ) :: File_Prefix CHARACTER( 256 ) :: SpcCoeff_File CHARACTER( 256 ) :: TauCoeff_File CHARACTER( 256 ) :: AerosolCoeff_File CHARACTER( 256 ) :: CloudCoeff_File CHARACTER( 256 ) :: EmisCoeff_File TYPE( CRTM_Atmosphere_type ) :: Atmosphere TYPE( CRTM_Surface_type ) :: Surface TYPE( CRTM_GeometryInfo_type ) :: GeometryInfo TYPE( CRTM_RTSolution_type ), ALLOCATABLE, DIMENSION( : ) :: RTSolution TYPE( CRTM_Options_type ) :: Options INTEGER, PARAMETER :: UNSET = 0 INTEGER :: CloudType, AerosolType, stype REAL( fp_kind ) :: a, alat, alon, height, scover, u, v real :: rlat,rlon ! local variable !============================================================== !----- ONLY READ COEFFICIENT FILES ONCE FOR EACH INSTRUMENT !============================================================== if (linit) then ! d_type=airs and subtype=420 is the Aqua satellite !there are six d_type if (d_type == 'AIRS_' .and. subtype ==420) then File_Prefix='airs281SUBSET_aqua' elseif (d_type == 'AIRS_' .and. subtype==570) then File_Prefix='amsua_aqua' elseif (d_type == 'HIRS2' ) then File_Prefix='hirs2_n14' elseif (d_type == 'HIRS3') then if (subtype ==207) File_Prefix='hirs3_n16' if (subtype ==208) File_Prefix='hirs3_n17' elseif (d_type == 'AMSUA') then !for AMSU A&B the subtype is the satellite ID if (subtype ==206) then File_Prefix='amsua_n15' elseif (subtype ==207) then File_Prefix='amsua_n16' elseif (subtype ==208) then File_Prefix='amsua_n17' elseif (subtype ==209) then File_Prefix='amsua_n18' else print *,' ' print ('(a,a5,a,i5,a)'),'Option d_type=',d_type, & ' with subtype=',subtype,' not permitted' stop endif elseif (d_type == 'AMSUB') then if (subtype ==206) then File_Prefix='amsub_n15' elseif (subtype ==207) then File_Prefix='amsub_n16' elseif (subtype ==208) then File_Prefix='amsub_n17' else print *,' ' print ('(a,a5,a,i5,a)'),'Option d_type=',d_type, & ' with subtype=',subtype,' not permitted' stop endif elseif (d_type == '_MSU_' ) then if (subtype ==205 ) File_Prefix='msu_n14' else print *,' ' print ('(a,a5,a)'),'Option d_type=',d_type,' not permitted' stop endif File_Prefix = ADJUSTL( File_Prefix ) ! -------------------- ! Create the filenames ! -------------------- SpcCoeff_File = TRIM( File_Prefix )//'.SpcCoeff.bin.Big_Endian' TauCoeff_File = TRIM( File_Prefix )//'.TauCoeff.bin.Big_Endian' AerosolCoeff_File = 'AerosolCoeff.bin.Big_Endian' CloudCoeff_File = 'CloudCoeff.bin.Big_Endian' EmisCoeff_File = 'EmisCoeff.bin.Big_Endian' !#----------------------------------------------------------------------------# !# -- INITIALISE THE CRTM -- # !#----------------------------------------------------------------------------# WRITE( *, '( /5x, "Initializing the CRTM..." )' ) write (6,*) 'SpcCoeff_File=', SpcCoeff_File write (6,*) 'TauCoeff_File=', TauCoeff_File Error_Status = CRTM_Init(ChannelInfo, & SpcCoeff_File = SpcCoeff_File, & TauCoeff_File = TauCoeff_File, & AerosolCoeff_File = AerosolCoeff_File, & EmisCoeff_File = EmisCoeff_File, & CloudCoeff_File = CloudCoeff_File ) IF ( Error_Status /= SUCCESS ) THEN CALL Display_Message( PROGRAM_NAME, & 'Error initializing CRTM', & Error_Status) STOP END IF ! print out channel information WRITE( *, '( /5x, "Number of channels indexed: ", i5 )' ) ChannelInfo%n_Channels WRITE( *, '( /2x, "Channel Sensor NCEP WMO WMO Channel", & &/2x, " Index Descriptor Sensor ID Satellite ID Sensor ID Number", & &/2x, "------------------------------------------------------------------------------" )' ) DO l = 1, ChannelInfo%n_Channels WRITE( *, '( 2x, 2x, i4, 2x, ">", a, "<", 5x, i3, 11x, i3, 11x, i3, 7x, i4 )' ) & ChannelInfo%Channel_Index( l ), & ChannelInfo%Sensor_Descriptor( l ), & ChannelInfo%NCEP_Sensor_ID( l ), & ChannelInfo%WMO_Satellite_ID( l ), & ChannelInfo%WMO_Sensor_ID( l ), & ChannelInfo%Sensor_Channel( l ) END DO !set linit=false linit=.false. endif !end of if (init) pi_fac = acos(-one) deg2rad = pi_fac/180.0 rad2deg = one/deg2rad if (nchan .ne. ChannelInfo%n_Channels) then stop 900 endif !#----------------------------------------------------------------------------# !# -- ALLOCATE THE OUTPUT RTSolution ARRAY TO THE NUMBER OF CHANNELS -- # !#----------------------------------------------------------------------------# Error_Status = CRTM_Allocate_Options( ChannelInfo%n_Channels, Options ) Atmosphere%n_Clouds = 0 Atmosphere%n_Absorbers = 2 Atmosphere%Climatology = 3 Atmosphere%n_Aerosols = 0 Atmosphere%n_Layers = lcldtop Error_Status = CRTM_Allocate_Atmosphere( Atmosphere%n_Layers, & Atmosphere%n_Absorbers, & Atmosphere%n_Clouds, & Atmosphere%n_Aerosols, & Atmosphere) IF ( Error_Status /= SUCCESS ) THEN CALL Display_Message(routine_name, & 'Error allocating atmosphere', & Error_Status) STOP END IF ! allocate optical depth array ALLOCATE( RTSolution( ChannelInfo%n_Channels ), STAT = Allocate_Status ) Error_Status = CRTM_Allocate_RTSolution( Atmosphere%n_Layers, RTSolution) IF ( Error_Status /= SUCCESS ) THEN CALL Display_Message(routine_name, & 'Error allocating RT and Options structure array', & Error_Status) STOP END IF Atmosphere%Absorber_ID(1) = 1 Atmosphere%Absorber_ID(2) = 3 Atmosphere%Absorber_Units(1) = 3 Atmosphere%Absorber_Units(2) = 1 alat = prof_info(2) ! units degree alon = prof_info(3) ! units degree if (alon .ge. 360.0) alon = alon - 360.0 if (alon .lt. 0.00) alon = alon + 360.0 height = prof_2d(nf_zsfc) ! in meter !-------------------------------------------------------------------------- ! compute local zenith angle and other angles in a same way as in geos-5 GSI. ! Note on an angle description, variable names used by this code and by GSI setuprad.f90: ! (setuprad.f90 reads in the data sets created by read_airs.f90 and read_bufrtovs.f90) ! ! (1) local zenith angle, or satellite zenith angle, ! Variable name here: lza, variable name in setuprad.f90: data_s(5) ! (2) look angle (scan angle) ! Variable name here: panglr, variable name in setuprad.f90: data_s(7) ! (3) solar zenith angle: ! Variable name here: SOZA, variable name in setuprad.f90: data_s(9) ! all angles are in degree when entering the CRTM !----------------------------------------------------------- if ( d_type(1:4) == 'HIRS' ) then step = 1.80 start = -49.5 rato=1.1363987 else if (d_type == '_MSU_' ) then step = 9.474 start = -47.37 rato=1.1363987 else if (d_type == 'AMSUA' ) then step = 3.+ 1.0/3.0 start = -48. - 1.0/3.0 else if ( d_type == 'AMSUB' ) then step = 1.1 start = -48.95 else if (d_type (1:4) == 'AIRS') then if (subtype ==420) then ! this is for airs step = 1.1 start = -48.9 elseif (subtype == 570) then ! this is for aqua_amsua step = 3. + 1./3. start = -48. - 1./3. endif endif ! data from input: ! prof_info(10)=SAZA ! prof_info(11)=fovs ! for all instruments considered, the common predictor parameters: panglr = (start+(prof_info(11)-1.0)*step)*deg2rad ! in radiance if (d_type(1:4) == 'HIRS' .or. d_type == '_MSU_') then ! separate hirs2 (msu) from hirs3 if (d_type == 'HIRS2' .or. d_type == '_MSU_') then lza = asin(rato*sin(panglr)) lza = lza*rad2deg ! in degree endif if (d_type == 'HIRS3') then lza = prof_info(10) ! in degree endif ! for all hirs if (d_type(1:4) == 'HIRS') then if( prof_info(11) <= 28.) lza=-lza endif elseif (d_type == 'AMSUA' .or. d_type == 'AMSUB' ) then lza=prof_info(10) ! in degree if((d_type == 'AMSUA' .and. prof_info(11) <= 15.) .or. & (d_type == 'AMSUB' .and. prof_info(11) <= 45.) ) & lza=-lza elseif (d_type == 'AIRS_' ) then ! include both airs and Aqua_amusa lza = prof_info(10) ! in degree if (prof_info(11) <=45.0 ) lza=-lza end if ! change lza back to the degree !-------------------------------------------------- ! for all five instruments: ! geometryinfo%sensor_zenith_angle = lza geometryinfo%source_zenith_angle= soza geometryinfo%sensor_scan_angle = panglr*rad2deg ! scan angle ! azimuth angle seems not in use GeometryInfo%Sensor_Azimuth_Angle = 0.0 ! GeometryInfo%Source_Azimuth_Angle = 0.0 ! if(abs(geometryinfo%sensor_zenith_angle) > one ) then geometryinfo%distance_ratio = & abs( sin(geometryinfo%sensor_scan_angle*deg2rad)/ & sin(geometryinfo%sensor_zenith_angle*deg2rad) ) endif u=prof_2d(nf_u10m) v=prof_2d(nf_v10m) ! ! copy input profile (above cloud only) and rescale ! change p from Pa to hPa ! change q from kg/kg to g/kg ! change ozone units ! top pressure is fixed as 0.005 mb Atmosphere%Level_Pressure(0) = 0.005 do i=1,lcldtop Atmosphere%Pressure(i) = prof_plevs(i)*0.01 Atmosphere%Level_Temperature(i) = prof_3d(i,nf_temp) Atmosphere%Temperature(i) = Atmosphere%Level_Temperature(i) Atmosphere%Absorber(i,1) = 1000.*max(prof_3d(i,nf_sphu), 1.e-7) Atmosphere%Absorber(i,2) = 604229.*max(prof_3d(i,nf_ozon), 1.e-10) enddo ! ! set pressure at interfaces (edges) do i=1,lcldtop-1 Atmosphere%Level_Pressure(i) = 0.5*(Atmosphere%Pressure(i+1) & + Atmosphere%Pressure(i)) enddo Atmosphere%Level_Pressure(lcldtop)=pcldtop*0.01 ! ! ! Begin setting of surface type information ! Surface%Land_Coverage = 0.0 Surface%Water_Coverage = 0.0 Surface%Snow_Coverage = 0.0 Surface%Ice_Coverage = 0.0 ! if ((nf_ismk == 0) .or. (nf_almk == 0)) then print *,' ' print *,' Some 2-d fields missing: ismk, almk =',nf_ismk,nf_almk print *,' These need to be provided from the nature run for use', & ' in the CRTM for MW calcs.' stop endif if (prof_2d(nf_ismk) > 0.2) then ! treat as ice point Surface%Ice_Coverage = 1.0 Surface%Ice_Type = 1 Surface%Ice_Temperature = & Atmosphere%Level_Temperature(Atmosphere%n_Layers) Surface%Ice_Thickness = 0.1 Surface%Ice_Density = 0.9 Surface%Ice_Roughness = 0.1 sfc_type=2 elseif (prof_2d(nf_almk) > 0.2) then ! treat as grassland Surface%Land_Coverage = 1.0 Surface%Land_Type = 11 Surface%Land_Temperature = & Atmosphere%Level_Temperature(Atmosphere%n_Layers) Surface%Soil_Moisture_Content = 0.1 Surface%Canopy_Water_Content = 0.05 Surface%Vegetation_Fraction = 0.5 Surface%Soil_Temperature = & Atmosphere%Level_Temperature(Atmosphere%n_Layers) sfc_type=3 else ! treat as sea water Surface%Water_Coverage = 1.0 Surface%Water_Type = 1 Surface%Water_Temperature = & Atmosphere%Level_Temperature(Atmosphere%n_Layers) Surface%Wind_Speed = sqrt(u*u+v*v) Surface%Wind_Direction = ACOS(u/Surface%Wind_Speed) * 180.0/3.14159 Surface%Salinity = 33.0 sfc_type=1 endif ! End setting of surface type information ! ! Call CRTM Error_Status = CRTM_Forward( Atmosphere, & Surface, & GeometryInfo, & ChannelInfo, & RTSolution) ! Options = Options ) if ( Error_Status /= SUCCESS ) then call Display_Message(routine_name, & 'Error in CRTM Forward Model',Error_Status) stop endif ! ! copy output into obs_values(1:n_channels) ! Runhua: check to see whther this is the full set of channels in the obs ! report. Otherwise we must make sure we put the info in the correct place. do i = 1, ChannelInfo%n_Channels l = ChannelInfo%Channel_Index(i) btchan(l)= RTSolution(l)%Brightness_Temperature enddo ! destroy the RTSolution and atmosphere ! Error_Status = CRTM_Destroy_RTSolution( RTSolution ) deallocate ( RTSolution, STAT = Allocate_Status ) if ( Error_Status /= SUCCESS ) then call Display_Message(PROGRAM_NAME, & 'Error destroying RTSolution structure', & Error_Status ) stop endif Error_Status = CRTM_Destroy_Atmosphere( Atmosphere ) if ( Error_Status /= SUCCESS ) then call Display_Message(PROGRAM_NAME, & 'Error destroying Atmosphere structure', & Error_Status ) stop endif end subroutine InitGetbt_crtm ! ! X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X ! subroutine find_name (nfields,f_names,f_name,nf,lstop) ! ! Find the index (nf) corresponding to the first occurance of the name ! of the field (f_name) within a chacter string (f_names). ! The name is a single character. ! ! input logical :: lstop ! stop if name not found in list integer :: nfields character(len=*) :: f_names(nfields) character(len=*) :: f_name ! ! output integer :: nf ! ! local integer :: n ! nf=0 do n=1,nfields if ( (nf == 0) .and. (f_names(n) == f_name) ) then nf=n endif enddo ! if (lstop .and. (nf==0)) then print *,' ' print *,' requested field name not found:' print *,' requested name = ',f_name print *,' possible names = ' print ('(10(1x,a4))'),f_names(1:nfields) stop endif ! end subroutine find_name end module m_interface_crtm