! $Id$ ! ! ! module rtm_mod use nwp_data_mod use rtm_constants use geo_mod use rtm_tools_mod use planck_mod use sfc_data_mod, only: sfc_main_type use rtm_types_mod use science_tools_mod, only : locate use sat_mod real, parameter :: Pi = 3.14159 real, parameter :: DTOR = PI / 180. real, parameter, private:: Co2_Ratio = 380.0 !in ppmv contains ! ! ! subroutine compute_rtm ( nwp , geo , sfc, satzen , sat , rtm ) implicit none type(nwp_main_type ) :: nwp type(geo_type) , target :: geo type(sfc_main_type) :: sfc real , dimension(:,:) , intent ( in) :: satzen type (rtm_main_type) , target :: rtm type ( sat_type ) :: sat integer :: i_rtm integer :: j_rtm integer :: i_sat integer :: j_sat integer :: min_idx_x integer :: max_idx_x integer :: min_idx_y integer :: max_idx_y real, allocatable :: rad_sfc_bb(:,:,:) integer :: num_sat_pxls integer :: Idx_Ang real :: Satzen_Mid_bin character ( len = 128) :: Ancil_Data_Dir integer :: kbin integer :: chan_idx real , dimension(2,2):: t_midlev , rad_midlev_dum real , dimension(2,2):: t_lev , rad_lev_dum real :: rad_midlev , rad_lev real, pointer ::trans_prof(:), t_prof_pnt(:), rad_prof(:) , rad_bb_prof(:) integer :: k integer :: idx_ang_min,idx_ang_max real :: sfc_rad integer, pointer :: i_nwp_x , i_nwp_y, i_angl integer :: i_sfc , i_tropo type (rtm_ngrid_chn_prof_type) , pointer :: ngrid_p type (rtm_sgrid_chn_data_type) , pointer :: sgrid_p ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! - create mapping of pfaast numbers to sensor channels call map_pfaast_chn ! - loop over all rtm min_idx_x = minval ( geo % idx_nwp_x) max_idx_x = maxval ( geo % idx_nwp_x) min_idx_y = minval ( geo % idx_nwp_y) max_idx_y = maxval ( geo % idx_nwp_y) call rtm%allocate_ngrid( nwp % n_lon, nwp % n_lat) rtm_x_loop: do i_rtm = min_idx_x, max_idx_x num_sat_pxls = count ( geo % idx_nwp_x == i_rtm ) if ( num_sat_pxls == 0) cycle rtm_x_loop rtm_y_loop: do j_rtm = min_idx_y, max_idx_y num_sat_pxls = count ( geo % idx_nwp_y == j_rtm .and. geo % idx_nwp_x == i_rtm ) if ( num_sat_pxls == 0) cycle rtm_y_loop ! allocate the atmospheric profiles call rtm % ngrid (i_rtm , j_rtm ) % allocate( nlevels_rtm) ! - transform to RTM profiles ! t_nwp, z_nwp, ozmr_nwp, rh_nwp !TODO make abetter input and output strutucre for this routine! call convert_profiles_NWP_RTM ( & & nwp % p_std & & , nwp % sfc_level(i_rtm, j_rtm) & & , nwp % psfc (i_rtm , j_rtm ) & & , nwp % wvmrsfc (i_rtm, j_rtm) & & , nwp % tmpair (i_rtm, j_rtm) & & , nwp % t_prof (i_rtm, j_rtm, :) & & , nwp % ozone_prof (i_rtm, j_rtm, :) & & , nwp % wvmr_prof (i_rtm, j_rtm, :) & & , nwp % z_prof (i_rtm, j_rtm,:) & & , rtm % ngrid (i_rtm, j_rtm) % z_prof & & , rtm % ngrid (i_rtm, j_rtm) % t_prof & & , rtm % ngrid (i_rtm, j_rtm) % wvmr_prof & & , rtm % ngrid (i_rtm, j_rtm) % ozmr_prof & ) ! allocate and fill profiles in rtm structure idx_ang_min = floor(50*minval(cosd(satzen) & , geo % idx_nwp_y == j_rtm .and. & geo % idx_nwp_x == i_rtm )) idx_ang_max = ceiling(50* maxval(cosd(satzen) & , geo % idx_nwp_y == j_rtm .and. & geo % idx_nwp_x == i_rtm)) rtm_ang_loop: do Idx_ang = idx_ang_min , idx_ang_max ! - check if this was already computed ( is this possible?) if ( rtm % ngrid (i_rtm, j_rtm) % angl(idx_ang) % is_set) cycle rtm_ang_loop rtm % ngrid (i_rtm, j_rtm) % angl(idx_ang) % is_set = .true. ! call rtm % ngrid (i_rtm, j_rtm) % angl(idx_ang) % allocate() !--- determine the zenith angle for the RTM Satzen_Mid_Bin = acosd((Idx_Ang-1)*Rtm_vza_binsize + Rtm_vza_binsize/2.0) Ancil_Data_Dir = '/DATA/Ancil_Data/clavrx_ancil_data/' rtm_chn_loop: do Chan_Idx = 20 , 36 if ( sat % chan_on (chan_idx) .eqv. .false. ) cycle kbin = pfaast_chn_idx(chan_Idx) if (kbin == 0 ) cycle rtm_chn_loop call rtm % ngrid (i_rtm, j_rtm) % angl (Idx_ang) % chn (chan_idx) % allocate( nlevels_rtm) t_prof_pnt => rtm % ngrid (i_rtm, j_rtm) % t_prof trans_prof => rtm % ngrid (i_rtm, j_rtm) % angl (Idx_ang) % chn (chan_idx) % trans_prof rad_prof => rtm % ngrid (i_rtm, j_rtm) % angl (Idx_ang) % chn (chan_idx) % rad_prof rad_bb_prof => rtm % ngrid (i_rtm, j_rtm) % angl (Idx_ang) % chn (chan_idx) % rad_bb_prof ! pfaast call tran_viirsm(Ancil_Data_Dir & , rtm % ngrid (i_rtm, j_rtm) % T_Prof & , rtm % ngrid (i_rtm, j_rtm) % Wvmr_Prof & , rtm % ngrid (i_rtm, j_rtm) % Ozmr_Prof & , Satzen_Mid_Bin & , Co2_Ratio & , kbin & , trans_prof ) level_loop: do k = 2 , nLevels_rtm !_---TODOTODOTDO PLANCK must be also for skalars ! and has to be checked if the values are right ! this is the reason for the ._dum arrays ! acronym lev stands for representative value for this level ( the mean between k and k+1) t_midlev(1,1) = 0.5 * ( t_prof_pnt(k-1) + t_prof_pnt (k) ) call planck_rad ( t_midlev , 'VIIRS' , chan_idx , rad_midlev_dum ) ! - see planck comment rad_midlev = rad_midlev_dum(1,1) rad_prof(k) = rad_prof(k-1) + (trans_prof(k-1) - trans_prof(k)) * rad_midlev t_lev(1,1) = t_prof_pnt (k) call planck_rad ( t_lev , 'VIIRS' , chan_idx , rad_lev_dum ) rad_lev = rad_lev_dum(1,1) rad_bb_prof(k) = rad_prof(k) + trans_prof(k) * rad_lev end do level_loop trans_prof => null() t_prof_pnt => null() rad_prof => null() rad_bb_prof => null() end do rtm_chn_loop end do rtm_ang_loop end do rtm_y_loop end do rtm_x_loop ! set key levels call set_key_levels ( rtm , nwp ) ! ! - pixel stuff ! - pixels has specific ! surface idx ! ! 1. allocate call rtm % allocate_sgrid( geo % n_x, geo % n_y) ! - compute black body BT from surface allocate ( rad_sfc_bb (20:36, geo % n_x, geo % n_y )) do chan_idx = 21, 36 call planck_rad ( nwp % sgrid % t_sfc , 'VIIRS' , chan_idx, rad_sfc_bb(chan_idx,:,:) ) end do rtm % sgrid % idx_angl = floor( 50 * cosd(satzen )) ! - loop over sat grid do i_sat = 1 , geo % n_x do j_sat = 1, geo % n_y i_nwp_x => geo%idx_nwp_x(i_sat,j_sat) i_nwp_y => geo%idx_nwp_y(i_sat,j_sat) i_angl => rtm % sgrid(i_sat,j_sat) % idx_angl ! - determine sfc level in profile ! initial set to the coarse nwp grid sfc i_sfc = rtm % ngrid (i_nwp_x, i_nwp_y) % level_sfc call locate ( rtm % ngrid (i_nwp_x, i_nwp_y) % z_prof ,nlevels_rtm & , (real(sfc % elevation % data (i_sat,j_sat) ) )/1000.0 , i_sfc ) i_tropo = rtm % ngrid (i_nwp_x, i_nwp_y) % level_tropo call rtm % sgrid (i_sat,j_sat) % allocate() do chan_idx = 20, 36 ngrid_p => rtm % ngrid (i_nwp_x, i_nwp_y) % angl (I_angl) %chn(Chan_idx) sgrid_p => rtm % sgrid(i_sat,j_sat) % chn(Chan_idx) if ( sat % chan_on (chan_idx) .eqv. .false. ) cycle if ( .not. ngrid_p % is_set ) cycle sgrid_p % rad_atm = ngrid_p % Rad_Prof (i_sfc) sgrid_p % trans_atm = ngrid_p % Trans_Prof (i_sfc) !-to do weights for exact sfc_rad = sfc % emis ( chan_idx) % data (i_sat,j_sat) * rad_sfc_bb(chan_idx,i_sat , j_sat ) sgrid_p % rad_atm_sfc = & sgrid_p % rad_atm + sgrid_p % trans_atm * sfc_rad rtm % ngrid (i_nwp_x, i_nwp_y) % angl (I_angl) % chn (chan_idx) % is_set = .true. if ( sat % data (chan_idx) % rad (i_sat,j_sat) > 0. .and. sfc % emis ( chan_idx) % data (i_sat,j_sat) > 0 ) then sgrid_p % emiss_tropo = emissivity ( sat % data (chan_idx) % rad (i_sat,j_sat) & & , sgrid_p % rad_atm_sfc & & , ngrid_p % Rad_bb_Prof (i_tropo)) ! - compute emissivities at tropopause print*, 'emiss tropo: ', sgrid_p % emiss_tropo if ( sgrid_p % emiss_tropo .lt. -1 ) then print*,chan_idx, i_sat,j_sat,shape (nwp % lon) print*,'sfc temp: ',nwp % sgrid(i_sat,j_sat) % t_sfc print*,' sfc_bb: ',rad_sfc_bb(chan_idx,i_sat , j_sat ) print*,'surface radiance:',sfc_rad print*, 'surface emis:',sfc % emis ( chan_idx) % data (i_sat,j_sat) print*,'input1: ',sat % data (chan_idx) % rad (i_sat,j_sat) print*,'input2: ',sgrid_p % rad_atm_sfc print*,'input3: ', ngrid_p % Rad_bb_Prof (i_tropo) print*,'a',sgrid_p % rad_atm print*,'trans: ',sgrid_p % trans_atm stop end if end if end do end do end do i_nwp_x => null() i_nwp_y => null() i_angl => null() sgrid_p => null() ngrid_p => null() deallocate ( rad_sfc_bb) end subroutine ! ! subroutine map_pfaast_chn ! for viirs pfaast_chn_idx = 0 pfaast_chn_idx (20) = 1 pfaast_chn_idx (22) = 2 pfaast_chn_idx (29) = 3 pfaast_chn_idx (31) = 4 pfaast_chn_idx (32) = 5 end subroutine map_pfaast_chn end module rtm_mod