PROGRAM ETL_RTTOV91 ! ! NAME: etl_rttov91.F90 ! ! PURPOSE: Uses rttov91 to perform forward radiative transfer of etl-formatted profiles ! ! INPUTS: (from command line) ! ! Platform name (ie. noaa,ssmi) ! satellite id (noaaxx = xx, ssmi = fxx) ! instrument number (determined by rttov_id routine) ! model type (Forward model only (0), full gradient (1), full radiance (2) ! profile id (integer number identifying profile number) ! surface type (0=land,1=sea,2=ice/snow) ! Zenith angle (local zenith angle in degrees) ! Azimuth angle (local azimuthal angle in degrees) ! ! NOTES: ! Uses ETL-formatted profile data as input. Each sounding is in individual file and read using ! etl_sounding.f90. ! ! CREATION HISTORY: ! - Taken from tstrad.F90 from version rttov8_5 ! - Darren Jackson CIRES/ESRL July 2008 ! !--------------------------------------------------------------------------------- ! Notes from tstrad.F90 code ! ! ! This software was developed within the context of ! the EUMETSAT Satellite Application Facility on ! Numerical Weather Prediction (NWP SAF), under the ! Cooperation Agreement dated 25 November 1998, between ! EUMETSAT and the Met Office, UK, by one or more partners ! within the NWP SAF. The partners in the NWP SAF are ! the Met Office, ECMWF, KNMI and MeteoFrance. ! ! Copyright 2002, EUMETSAT, All Rights Reserved. ! ! ************************************************************* ! ! TEST PROGRAM FOR RTTOV SUITE. ! RTTOV VERSION 9 ! ! Description: This program is the test harness for RTTOV-8. There ! are 3 options: ! option = 0 to test forward model only ! option = 1 to test the full model ie TL/AD/K ! ! To run this program you must have the following files ! either resident in the same directory or set up as a ! symbolic link: ! refprof.dat -- reference profile ! prof.dat -- input profile ! input.dat -- file to select channels and surface emis ! rtcoef_platform_id_sensor.dat -- coefficient file to match ! the sensor you request in the input dialogue ! There are unix scripts available to set up the files above and ! run this program (e.g. tstrad_full.scr) ! The output is generated in a file called print.dat. ! This output can be compared with reference output generated ! by the code developers and included with the export package. ! ! Current Code Owner: SAF NWP ! ! History: ! Version Date Comment ! ------- ---- ------- ! 25/01/2002 Initial version (R. Saunders) ! 01/05/2002 Updated for NOAA-17 (R. Saunders) ! 01/12/2002 New F90 code with structures (P Brunel A Smith) ! 02/01/2003 Comments added (R Saunders) ! 10/12/2003 Updated for polarimetric changes (S. English/R.Saunders) ! 01/04/2004 Updated for chan setup routines (R.Saunders) ! 01/06/2005 Marco Matricardi (ECMWF): ! -- IASI capability added. ! -- Variable trace gases added for IASI and AIRS. ! -- Solar radiation added for IASI and AIRS ! -- Altitude dependent local zenith angle added. ! -- Linear in tau approximation for RTequation added. ! 22/01/2007 Removed polarisation indexing. R Saunders ! 14/03/2007 Cleaned code + added more options (R Saunders) ! 11/10/2007 Make parallel version. Extended interface for ! rttov_alloc_prof. P.Marguinaud ! 18/10/2007 Add extra outputs for Phil Watts ! ! Code Description: ! Language: Fortran 90. ! Software Standards: "European Standards for Writing and ! Documenting Exchangeable Fortran 90 Code". ! Use rttov_const, only : & & nplatforms ,& & ninst ,& & pi ,& & errorstatus_fatal ,& & errorstatus_warning ,& & errorstatus_success ,& & platform_name , & & sensor_id_mw ,& & inst_name, & & ncldtyp , & & naer_cl, & & q_mixratio_to_ppmv USE YOMHOOK ,ONLY : LHOOK, DR_HOOK Use rttov_types, only : & & rttov_coef , & & predictors_type , & & profile_type , & & transmission_Type , & & radiance_type , & & geometry_type , & & rttov_coef_scatt_ir , & & rttov_optpar_ir Use mod_tstrad ! Use parkind1, Only : jpim ,jprb Use etl_sounding Implicit None #include "rttov_errorreport.interface" #include "rttov_setup.interface" #include "rttov_errorhandling.interface" #ifdef _RTTOV_TSTRAD_PARALLEL #include "rttov_parallel_direct.interface" #define rttov_direct rttov_parallel_direct #else #include "rttov_direct.interface" #endif #include "rttov_alloc_prof.interface" #include "rttov_alloc_rad.interface" #include "rttov_dealloc_coef.interface" #include "tstrad_tl.interface" #include "tstrad_ad.interface" #include "tstrad_k.interface" ! type( rttov_coef ), allocatable :: coef(:) ! coefficients type(profile_type), allocatable :: profiles(:) type(transmission_type) :: transmission type(radiance_type) :: radiance type(predictors_type) :: predictors type(rttov_coef_scatt_ir),allocatable:: coef_scatt_ir(:) type(rttov_optpar_ir),allocatable :: optp(:) ! Logical :: addinterp ! switch for the interpolator Integer(Kind=jpim), Allocatable :: instrument(:,:) ! instrument id Integer(Kind=jpim), Allocatable :: nchan(:,:) ! number of channels per instrument and profile Integer(Kind=jpim), Allocatable :: nchan1(:) ! number of channels per instrument (first profile) Integer(Kind=jpim), Allocatable :: nchannels(:) ! number of channels per instrument Integer(Kind=jpim), Allocatable :: nchan_byprof(:) ! number of channels per profile (specified instrument) Integer(Kind=jpim), Allocatable :: ifull(:) ! full test (with TL,AD,K) per instrument Integer(Kind=jpim), Allocatable :: nprof(:) ! number of profiles per instrument Integer(Kind=jpim), Allocatable :: nsurf(:) ! surface id number per instrument Integer(Kind=jpim), Allocatable :: nwater(:) ! water id number per instrument Real(Kind=jprb), Allocatable :: surfem(:,:) ! surface input emissivity per channel , instrument Integer(Kind=jpim), Allocatable :: ichan(:,:) ! channel list per instrument Integer(Kind=jpim), Allocatable :: channels(:) ! channel list per instrument*profiles Integer(Kind=jpim), Allocatable :: lprofiles (:) Real(Kind=jprb), Allocatable :: emissivity (:) Real(Kind=jprb), Allocatable :: fresnrefl (:) Real(Kind=jprb), Allocatable :: input_emissivity (:) logical, Allocatable :: calcemis (:) Integer(Kind=jpim) :: nprofiles Integer(Kind=jpim) :: iref Integer(Kind=jpim) :: isun Integer(Kind=jpim) :: iaer Integer(Kind=jpim) :: icld Integer(Kind=jpim) :: iclddg Integer(Kind=jpim) :: iclshp Integer(Kind=jpim) :: icla Integer(Kind=jpim) :: icoa Integer(Kind=jpim) :: asw Integer(Kind=jpim) :: coef_errorstatus ! read coeffs error return code Integer(Kind=jpim), Allocatable :: rttov_errorstatus(:) ! rttov error return code Integer(Kind=jpim), Allocatable :: setup_errorstatus(:) ! setup return code ! min and max satellite id for each platform Integer(Kind=jpim), dimension(nplatforms) :: max_satid Integer(Kind=jpim), dimension(nplatforms) :: min_satid ! min and max channel numbers for each instrument integer(Kind=jpim), dimension(0:ninst-1) :: max_channel_old integer(Kind=jpim), dimension(0:ninst-1) :: max_channel_new integer(Kind=jpim), dimension(0:ninst-1) :: max_channel integer(Kind=jpim), parameter :: mxchn = 9000 ! max number of channels per instruments allowed in one run ! printing arrays real(Kind=jprb), Allocatable :: pr_radcld(:) real(Kind=jprb), Allocatable :: pr_trans(:) real(Kind=jprb), Allocatable :: pr_emis(:) real(Kind=jprb), Allocatable :: pr_trans_lev(:,:) real(Kind=jprb), Allocatable :: pr_upclr(:) real(Kind=jprb), Allocatable :: pr_ovcst(:,:) real(Kind=jprb), Allocatable :: pr_up(:,:) real(Kind=jprb), Allocatable :: pr_down(:,:) real(Kind=jprb), Allocatable :: pr_surf(:,:) integer(Kind=jpim), dimension(1:mxchn) :: pr_pol ! Arrays from profile data REAL(Kind=jprb), ALLOCATABLE, DIMENSION ( : ), TARGET :: p REAL(Kind=jprb), ALLOCATABLE, DIMENSION ( : ), TARGET :: t REAL(Kind=jprb), ALLOCATABLE, DIMENSION ( : ), TARGET :: q REAL(Kind=jprb), ALLOCATABLE, DIMENSION ( : ), TARGET :: o3 data min_satid / 1, 8, 1, 8, 5, 2, 1, 1, 1, 1, 1, 1, 3, 2, 1, 1, 1, 1, 0, 1 / data max_satid /18,16, 7,12, 5, 3, 1, 2, 3, 2, 1, 2, 4, 2, 1, 1, 1, 1, 0, 1/ data max_channel_old / 20, 4, 3, 15, 5, 3, 7, 8, 8, 9,& & 24, 2378, 4, 16, 3, 5, 8461,14, 4,22,& & 2, 8, 4, 18, 4, 3, 3,1000, 40, 22, & & 16, 3000, 0, 0, 6, 9/ data max_channel_new / 20, 4, 3, 15, 5, 3, 4, 8, 8, 9,& & 21, 2378, 4, 16, 3, 5, 8461,14, 4,22,& & 2, 8, 4, 18, 4, 3, 3,1000, 40, 22, & & 16, 3000, 0, 0, 6, 9/ Character (len=3) :: cref Character (len=3) :: caer Character (len=3) :: ccld Character (len=3) :: csun Character (len=80) :: errMessage Character (len=6) :: NameOfRoutine = 'frt_tigr91' Character (len=3) :: coeff_version = 'old' ! Integer(Kind=jpim) :: Err_Unit ! Logical error unit (<0 for default) Integer(Kind=jpim) :: verbosity_level ! (<0 for default) Integer(Kind=jpim) :: nrttovid ! maximum number of instruments Integer(Kind=jpim) :: no_id ! instrument loop index Integer(Kind=jpim) :: nlevels, nformat Integer(Kind=jpim) :: ios integer(Kind=jpim) :: i,pol_id,ich2 integer(Kind=jpim) :: ichannels Integer(Kind=jpim) :: j Integer(Kind=jpim) :: jjm, ira, jj,icl integer(Kind=jpim) :: jch, jpol integer(Kind=jpim) :: jn, joff1, joff2, joff3 Integer(Kind=jpim) :: nprint Integer(Kind=jpim) :: ncomp Integer(Kind=jpim) :: np, ilev Integer(Kind=jpim) :: n Integer(Kind=jpim) :: nch ! intermediate variable Integer(Kind=jpim) :: ich ! intermediate variable Integer(Kind=jpim) :: ii ! intermediate variable Integer(Kind=jpim) :: errorstatus Integer(Kind=jpim) :: io_status Real(Kind=jprb) :: s Real(Kind=jprb) :: del_q,del_t Real(Kind=jprb) :: zenang Real(Kind=jprb) :: azang Real(Kind=jprb) :: lat Real(Kind=jprb) :: zerht Real(Kind=jprb) :: sunzang Real(Kind=jprb) :: sunazang logical :: refrac logical :: solrad logical :: laerosl logical :: lclouds logical :: lsun logical :: all_channels Integer(Kind=jpim) :: iua Integer(Kind=jpim) :: ioout Integer(Kind=jpim) :: iue Integer(Kind=jpim) :: interp ! local Integer(Kind=jpim), Parameter :: jpnav = 31 ! no. of profile variables Integer(Kind=jpim), Parameter :: jpnsav = 6 ! no. of surface air variables Integer(Kind=jpim), Parameter :: jpnssv = 6 ! no. of skin variables Integer(Kind=jpim), Parameter :: jpncv = 2 ! no. of cloud variables Integer(Kind=jpim), Parameter :: sscvar = jpnsav+jpnssv+jpncv ! no of surface,skin,cloud vars ! Unit numbers for input/output DATA IUA/1/,IOOUT/2/,IUE/56/ Integer(Kind=jpim) :: alloc_status(60) REAL(KIND=JPRB) :: ZHOOK_HANDLE !- End of header -------------------------------------------------------- IF (LHOOK) CALL DR_HOOK('TSTRAD',0,ZHOOK_HANDLE) errorstatus = 0 alloc_status(:) = 0 all_channels = .false. sunzang = 0._jprb sunazang = 0._jprb ! !Initialise error management with default value for ! the error unit number and ! Fatal error message output Err_unit = -1 verbosity_level = 1 ! All error message output verbosity_level = 3 call rttov_errorhandling(Err_unit, verbosity_level, print_checkinput_warnings=.false.) io_status = 0 errMessage = '' ! Beginning of Routine. ! --------------------- OPEN(IOOUT,file='print.dat',status='unknown',form='formatted') ! For more than one satellite ! comment out the next line and uncomment the following two. NRTTOVID = 1 ! PRINT *, 'How many satellites do you want?' ! READ *, NRTTOVID allocate (coef(nrttovid),stat= alloc_status(1)) allocate (coef_scatt_ir(nrttovid),stat= alloc_status(2)) allocate (optp(nrttovid),stat= alloc_status(3)) allocate (instrument(3,nrttovid),stat= alloc_status(4)) allocate (ifull(nrttovid),stat= alloc_status(5)) allocate (nprof(nrttovid),stat= alloc_status(6)) allocate (nsurf(nrttovid),stat= alloc_status(7)) allocate (nwater(nrttovid),stat= alloc_status(8)) allocate (nchannels(nrttovid),stat= alloc_status(9)) allocate (nchan1(nrttovid),stat= alloc_status(10)) !maximum number of channels allowed for one instrument is mxchn allocate (surfem(mxchn,nrttovid),stat= alloc_status(11)) allocate (ichan (mxchn,nrttovid),stat= alloc_status(12)) If( any(alloc_status /= 0) ) then errorstatus = errorstatus_fatal Write( errMessage, '( "mem allocation error")' ) Call Rttov_ErrorReport (errorstatus, errMessage, NameOfRoutine) Stop End If surfem(:,:) = 0.0_JPRB ichan(:,:) = 0 DO NO_ID = 1, NRTTOVID write(*,*) 'Which satellite platform do you want?' WRITE(*,'(4(2x,i3,2x,a8))') (i,platform_name(i), i = 1, nplatforms) READ *, Instrument(1,no_id) IF ( Instrument(1,no_id) <= 0 .OR. & & Instrument(1,no_id) > nplatforms) STOP 'Platform number not allowed' WRITE(*,*) 'Which satellite id do you want for this platform?' WRITE(*,*) 'Noaaxx = xx GOESyy = yy' READ *, instrument(2,no_id) if( instrument(2,no_id) < min_satid(Instrument(1,no_id)) .or. & & instrument(2,no_id) > max_satid(Instrument(1,no_id)) ) & & STOP 'Satellite id not allowed' WRITE(*,*) 'Which instrument type do you want for this satellite?' write(*, '(4(2x,i3,2x,a8))') (i, inst_name(i), i = 0, ninst-1) READ *, instrument(3,no_id) IF ( instrument(3,no_id) < 0 .OR. & & instrument(3,no_id) > ninst-1)& & STOP 'instrument number not allowed' WRITE(*,*) ' Forward model only (0) or full gradient test (1)',& & ' or full radiance output (2)?' READ *, IFULL(no_id) IFULL(no_id)=0 PRINT *, ' Surface type (0=land, 1=sea, 2=ice/snow)? ' READ *, NSURF(no_id) PRINT *, ' Water type (0=fresh water, 1=ocean water)' READ *, NWATER(no_id) ! PRINT *, ' Number of profiles to test code? ' NPROF(no_id) = 1 nprofiles = NPROF(no_id) ! !..SET UP CHANNEL NUMBERS ! ! .. DEFAULT MAXIMUMS if (coeff_version == 'old') max_channel(:)=max_channel_old(:) if (coeff_version == 'new') max_channel(:)=max_channel_new(:) allocate (nchan(nprof(no_id),nrttovid), stat= alloc_status(1)) nchan(1:nprof(no_id),no_id) = max_channel(instrument(3,no_id)) ! ! Note that channels are the same for all instruments ! and all profiles because the filename is the same OPEN (IUE,FILE='input.dat',status='old') READ(IUE,*) NCH = 0 DO ICH = 1 , nchan(1,no_id) READ(IUE,*,iostat=ios)I,II,S if(ios /= 0 ) then write (*,*) ' TOO FEW CHANNELS IN INPUT FILE ' write (*,*) ' nchan(1,no_id),no_id ',nchan(1,no_id),no_id stop endif IF(II.GT.0)THEN NCH = NCH + 1 ICHAN(nch,no_id) = I SURFEM(nch,no_id) = s ENDIF ENDDO ! CLOSE(IUE) ! nchan(1,no_id) is now the real number of channels selected do j = 1 , nprof(no_id) nchan(j,no_id) = MIN(max_channel(instrument(3,no_id)),NCH) enddo write(6,*)' Number of channels selected = ',nchan(1,no_id) ! compute channels*profiles nchannels(no_id) = 0 write(6,*)'nprofiles = ',nprofiles Do j = 1 , nprofiles nchannels(no_id) = nchannels(no_id) + nchan (j,no_id) write(6,*) 'j,nchannels(no_id)=',j,nchannels(no_id) End Do nchan1(no_id) = nchan(1,no_id) ! END DO ! write(6,'(/)') ! write(6,*)'Do you want to include aerosols?' ! write(6,*)'0=No , 1=Yes' ! read(5,*)iaer iaer=0 if(iaer==1)then write(6,'(a)')'Do you want to use the default optical parameters?' write(6,*)'0=No , 1=Yes' read(5,*)icoa if(icoa==1)then WRITE(6,'(a)')'Do you want to use a climatological profile for aerosols' WRITE(6,'(a)')' 0=No' WRITE(6,'(a)')' 1=Continental clean' WRITE(6,'(a)')' 2=Continental average' WRITE(6,'(a)')' 3=Continental polluted' WRITE(6,'(a)')' 4=Urban' WRITE(6,'(a)')' 5=Desert' WRITE(6,'(a)')' 6=Maritime clean' WRITE(6,'(a)')' 7=Maritime polluted' WRITE(6,'(a)')' 8=Maritime tropical' WRITE(6,'(a)')' 9=Arctic' WRITE(6,'(a)')'10=Antarctic' read(5,*)icla endif endif ! write(6,'(/)') ! write(6,*)'Do you want to include clouds?' ! write(6,*)'0=No , 1=Yes' ! read(5,*)icld icld=0 if(iaer==0)then laerosl=.False. else if(iaer==1)then laerosl=.True. endif if(icld==1)then lclouds=.True. write(6,'(/)') write(6,*)'Choose the scheme to convert IWC to effective generalized diameter,Dg' write(6,*)'1=Ou and Liuou scheme' write(6,*)'2=Wyser et a. scheme ' write(6,*)'3=Boudala et al.scheme' write(6,*)'4=McFarquhar et al. scheme' read(5,*)iclddg write(6,'(/)') write(6,*)'Choose the shape of ice crystals' write(6,*)'1=Hexagonal ice crystals' write(6,*)'2=Ice aggregates' read(5,*)iclshp else if (icld==0)then lclouds=.False. endif !--------------------------------------------------------- ! Beginning of rttov_setup test !--------------------------------------------------------- alloc_status = 0 allocate ( setup_errorstatus(nrttovid),stat= alloc_status(1)) If( any(alloc_status /= 0) ) then errorstatus = errorstatus_fatal Write( errMessage, '( "mem allocation error for errorsetup")' ) Call Rttov_ErrorReport (errorstatus, errMessage, NameOfRoutine) Stop End If If (all_channels)Then Call rttov_setup (& & setup_errorstatus, &! out & Err_unit, &! in & verbosity_level, &! in & nrttovid, &! in & laerosl, &! in & lclouds, &! in & coef, &! out & coef_scatt_ir, &! out & optp, & & instrument) ! in Else Call rttov_setup (& & setup_errorstatus, &! out & Err_unit, &! in & verbosity_level, &! in & nrttovid, &! in & laerosl, &! in & lclouds, &! in & coef, &! out & coef_scatt_ir, &! out & optp, & & instrument, &! in & ichan ) ! in Optional Endif if(any(setup_errorstatus(:) /= errorstatus_success ) ) then write ( ioout, * ) 'rttov_setup fatal error' stop endif deallocate( setup_errorstatus ,stat=alloc_status(1)) If( any(alloc_status /= 0) ) then errorstatus = errorstatus_fatal Write( errMessage, '( "mem deallocation error for setup_errorstatus")' ) Call Rttov_ErrorReport (errorstatus, errMessage, NameOfRoutine) Stop End If DO no_id = 1, NRTTOVID if( any(coef(no_id)%ff_val_chn( : ) /= 1 )) then WRITE(*,*) ' some requested channels have bad validity parameter' do i = 1, nchan1(no_id) write(*,*) i, coef(no_id)%ff_val_chn(i) end do endif End Do !--------------------------------------------------------- ! End of rttov_setup test !--------------------------------------------------------- ! DO no_id = 1, NRTTOVID !----------------------------------- ! Memory allocation for RTTOV_Direct !----------------------------------- allocate( rttov_errorstatus(nprof(no_id)),stat= alloc_status(3)) ! Allocate profiles asw = 1 ! allocate ! On USER levs allocate( profiles(nprof(no_id)),stat= alloc_status(1)) ! Turn off option for for reading this information ! OPEN (IUA,FILE='plevs.dat',status='old') ! read(IUA,*) ! read(IUA,*) interp, & ! logical switch for the interpolator ! & nlevels ! number of USER pressure levels for this run ! close(IUA) interp = 0 nlevels = 43_jpim if(interp == 0) addinterp = .false. if(interp == 1) addinterp = .true. call rttov_alloc_prof ( & & errorstatus, & & nprof(no_id), & & profiles, & & nlevels, & & coef_scatt_ir(no_id), & & asw, & & addclouds = lclouds, & & addaerosl = laerosl, & & init=.true. ) profiles(1) % nlevels = nlevels Do j = 1 , nprof(no_id) profiles(j) % nlevels = profiles(1) % nlevels Enddo alloc_status = 0_jpim ! number of channels per RTTOV call is only nchannels allocate( lprofiles ( nchannels(no_id) ) ,stat= alloc_status(9)) allocate( channels ( nchannels(no_id) ) ,stat= alloc_status(10)) allocate( nchan_byprof ( nprof(no_id) ) ,stat= alloc_status(11)) allocate( emissivity ( nchannels(no_id) ) ,stat= alloc_status(12)) allocate( fresnrefl ( nchannels(no_id) ) ,stat= alloc_status(13)) allocate( input_emissivity ( nchannels(no_id) ) ,stat= alloc_status(14)) allocate( calcemis ( nchannels(no_id) ) ,stat= alloc_status(15)) ! allocate transmittance arrays with number of channels allocate( transmission % tau_layers (profiles(1) % nlevels,nchannels(no_id) ) ,stat= alloc_status(11)) allocate( transmission % tau_total (nchannels(no_id) ) ,stat= alloc_status(12)) If( Any(alloc_status /= 0) ) Then errorstatus = errorstatus_fatal Write( errMessage, '( "allocation of transmission")' ) Call Rttov_ErrorReport (errorstatus_fatal, errMessage, NameOfRoutine) IF (LHOOK) CALL DR_HOOK('RTTOV_DIRECT',1,ZHOOK_HANDLE) Stop End If transmission % tau_layers = 0._jprb transmission % tau_total = 0._jprb ! allocate radiance results arrays with number of channels asw = 1 ! allocate call rttov_alloc_rad (errorstatus,nchannels(no_id),radiance,profiles(1)%nlevels,asw) ! allocate the cloudy radiances with full size even if not used ! Save input values of emissivities for all calculations. Allocate(pr_radcld(nchannels(no_id)) ,stat= alloc_status(33)) Allocate(pr_trans(nchannels(no_id)) ,stat= alloc_status(34)) Allocate(pr_emis(nchannels(no_id)) ,stat= alloc_status(35)) Allocate(pr_trans_lev(profiles(1) % nlevels,nchannels(no_id)) ,stat= alloc_status(36)) Allocate(pr_upclr(nchannels(no_id)) ,stat= alloc_status(37)) Allocate(pr_ovcst(profiles(1) % nlevels,nchannels(no_id)) ,stat= alloc_status(40)) Allocate(pr_up(profiles(1) % nlevels,nchannels(no_id)) ,stat= alloc_status(41)) Allocate(pr_down(profiles(1) % nlevels,nchannels(no_id)) ,stat= alloc_status(42)) Allocate(pr_surf(profiles(1) % nlevels,nchannels(no_id)) ,stat= alloc_status(43)) If( any(alloc_status /= 0) ) then errorstatus = errorstatus_fatal Write( errMessage, '( "mem allocation error")' ) Call Rttov_ErrorReport (errorstatus, errMessage, NameOfRoutine) Stop End If WRITE(6,*)'Zenith angle (degrees)?' READ(5,*)ZENANG WRITE(6,*)'Azimuth angle (degrees)?' READ(5,*)AZANG ! WRITE(6,*)'Do you want to include solar radiation?' ! WRITE(6,*)'0=No , 1=Yes' ! READ(5,*)ISUN isun=0 if(isun==1)then WRITE(6,*)'Solar Zenith angle (degrees)?' READ(5,*)SUNZANG WRITE(6,*)'Solar Azimuth angle (degrees)?' READ(5,*)SUNAZANG endif ! WRITE(6,*)'Latitude of sounding (degrees)?' ! READ(5,*)LAT lat=0._jprb ! WRITE(6,*)'Elevation of sounding (km)?' ! READ(5,*)ZERHT zerht=0._jprb ! WRITE(6,*)'Do you want to include atmospheric refraction?' ! WRITE(6,*)'0=No , 1=Yes' ! READ(5,*)IREF iref=0 if(iref==0)then cref='NO' refrac=.False. else if(iref==1)then cref='YES' refrac=.True. endif if(sunzang<=87._JPRB)then solrad=.True. else solrad=.False. endif if(isun==1)then lsun=.true. if(sunzang<=87._JPRB)then csun='YES' solrad=.True. else csun='NO' solrad=.False. endif else csun='NO' solrad=.False. endif if(iaer==1)then caer='YES' profiles(1)%aer_data=.true. if(icla==0)then WRITE(6,*)' Aerosols components =',coef_scatt_ir(no_id) % fmv_aer_comp OPEN (IUA,FILE='prof_aerosl.dat',status='old') DO IRA=1,profiles(1)% nlevels read(iua,*)(profiles(1)%aerosols(I,IRA),I=1,coef_scatt_ir(no_id) % fmv_aer_comp) ! write(*,'(11f11.3)')(profiles(1)%aerosols(I,IRA),I=1,coef_scatt_ir(no_id) % fmv_aer_comp) ENDDO close(iua) else if(icla>=1.and.icla<=10)then WRITE(6,*)' Aerosols types =',naer_cl OPEN (IUA,FILE='prof_aerosl_cl.dat',status='old') DO ICL=1,naer_cl DO IRA=1,profiles(1)% nlevels if(icl/=icla)then READ(iua,*) else READ(iua,*)(profiles(1)%aerosols(I,IRA),I=1,coef_scatt_ir(no_id) % fmv_aer_comp) ! write(*,'(11f11.3)')(profiles(1)%aerosols(I,IRA),I=1,coef_scatt_ir(no_id) % fmv_aer_comp) endif ENDDO ENDDO close(iua) endif else caer='NO' profiles(1)%aer_data=.false. if( laerosl ) & & profiles(1)%aerosols(:,:)=0._jprb endif if(icld==1)then ccld='YES' profiles(1)%idg=iclddg profiles(1)%ish=iclshp profiles(1)%cld_data=.true. WRITE(6,*)' Cloud types =',coef_scatt_ir(no_id) % fmv_wcl_comp+1 OPEN (IUA,FILE='prof_cloud.dat',status='old') DO IRA=1,profiles(1)% nlevels read(iua,*)(profiles(1)%cloud(I,IRA),I=1,coef_scatt_ir(no_id) % fmv_wcl_comp+1) ! write(*,'(11f11.3)')(profiles(1)%cloud(I,IRA),I=1,coef_scatt_ir(no_id) % fmv_wcl_comp+1) ENDDO write(*,'(/)') LOOP1: DO ira=1,profiles(1)% nlevels DO I=1,coef_scatt_ir(no_id) % fmv_wcl_comp+1 DO J=I,coef_scatt_ir(no_id) % fmv_wcl_comp+1 IF(I/=J)THEN IF(profiles(1)%cloud(I,ira)/= 0.AND.profiles(1)%cloud(J,ira)/= 0.)THEN WRITE(*,'(A65,I4)')'Error:two or more different cloud typ& & es are not allowed in layer',ira io_status=2 CYCLE LOOP1 ENDIF ENDIF ENDDO ENDDO ENDDO LOOP1 if(io_status/=0)then Call Rttov_ErrorReport (io_status, errMessage, NameOfRoutine) errorstatus = errorstatus_fatal IF (LHOOK) CALL DR_HOOK('TSTRAD',1,ZHOOK_HANDLE) Stop endif close(iua) OPEN (IUA,FILE='prof_cfrac.dat',status='old') DO IRA=1,profiles(1)% nlevels read(iua,*)(profiles(1)%cfrac(I,IRA),I=1,coef_scatt_ir(no_id) % fmv_wcl_comp+1) ! write(*,'(11f11.3)')(profiles(1)%cfrac(I,IRA),I=1,coef_scatt_ir(no_id) % fmv_wcl_comp+1) ENDDO LOOP2: DO ira=1,profiles(1)% nlevels DO I=1,coef_scatt_ir(no_id) % fmv_wcl_comp+1 DO J=I,coef_scatt_ir(no_id) % fmv_wcl_comp+1 IF(I/=J)THEN IF(profiles(1)%cfrac(I,ira)/= 0.AND.profiles(1)%cfrac(J,ira)/= 0.)THEN WRITE(*,'(A71,I4)')'Error:fractional coverage for two or & & more cloud types is given in layer',& & ira io_status=2 CYCLE LOOP2 ENDIF ENDIF ENDDO ENDDO ENDDO LOOP2 if(io_status/=0)then Call Rttov_ErrorReport (io_status, errMessage, NameOfRoutine) errorstatus = errorstatus_fatal IF (LHOOK) CALL DR_HOOK('TSTRAD',1,ZHOOK_HANDLE) Stop endif LOOP3: DO ira=1,profiles(1)% nlevels DO I=1,coef_scatt_ir(no_id) % fmv_wcl_comp+1 IF(profiles(1)%cfrac(I,ira)== 0.AND.profiles(1)%cloud(I,ira)/= 0.)THEN WRITE(*,'(A54,I4)')'Error:fractional coverage not given & & for cloud in layer',ira io_status=2 CYCLE LOOP3 ENDIF ENDDO ENDDO LOOP3 if(io_status/=0)then Call Rttov_ErrorReport (io_status, errMessage, NameOfRoutine) errorstatus = errorstatus_fatal IF (LHOOK) CALL DR_HOOK('TSTRAD',1,ZHOOK_HANDLE) Stop endif LOOP4: DO ira=1,profiles(1)% nlevels DO I=1,coef_scatt_ir(no_id) % fmv_wcl_comp+1 IF(profiles(1)%cfrac(I,ira)/= 0.AND.profiles(1)%cloud(I,ira)== 0.)THEN WRITE(*,'(A55,I4)')'Error:liquid water content not given & & for cloud in layer',ira io_status=2 CYCLE LOOP4 ENDIF ENDDO ENDDO LOOP4 if(io_status/=0)then Call Rttov_ErrorReport (io_status, errMessage, NameOfRoutine) errorstatus = errorstatus_fatal IF (LHOOK) CALL DR_HOOK('TSTRAD',1,ZHOOK_HANDLE) Stop endif close(iua) else ccld='NO' profiles(1)%idg =0._jprb profiles(1)%ish =0._jprb profiles(1)%cld_data =.false. endif WRITE(6,*) WRITE(6,'('' number of levels input = '',i4)') profiles(1)% nlevels if(addinterp) then WRITE(6,*) 'call interpolator - Yes' else WRITE(6,*) 'call interpolator - No' endif ! Read pressure levels from separate file ! OPEN (IUA,FILE='plevs.dat',status='old') ! read(IUA,*) ! aleady read ! nformat=10 ! how many fields in one line ! JJM = 0 ! DO IRA = 1 , 1+(profiles(1)% nlevels-1)/nformat ! JJ = 1+JJM ! JJM = MIN(JJ+nformat-1,profiles(1)% nlevels) ! READ(IUA,*) (profiles(1) % p(J),J=JJ,JJM) ! enddo ! CLOSE(IUA) ! ! Read ETL formatted profile data ! ALLOCATE( p(nlevels), & t(nlevels), & o3(nlevels), & q(nlevels) ) CALL read_etl_sounding('prof.dat', & nlevels, & p, & t, & q, & o3=o3 ) ! ! Change units of specific humidity to ppmv, temperature to Kelvin ! q(:) = q(:) / 1000._JPRB * q_mixratio_to_ppmv t(:) = t(:) + 273.15_JPRB ! ! Place TIGR profile information into profile structure ! profiles(1) % t => t profiles(1) % o3 => o3 profiles(1) % q => q profiles(1) % p => p profiles(1) % ozone_data = .true. profiles(1) % co2_data = .false. profiles(1) % clw_data = .false. profiles(1) % ch4_data = .false. profiles(1) % n2o_data = .false. profiles(1) % ch4_data = .false. del_t = profiles(1) % t(43) - profiles(1) % t(42) del_q = profiles(1) % q(43) - profiles(1) % q(42) profiles(1) % s2m % t = t(43) - del_t * 0.76087 profiles(1) % s2m % q = q(43) - del_q * 0.76087 profiles(1) % s2m % p = 1007.3 profiles(1) % s2m % u = 8. profiles(1) % s2m % v = 8. profiles(1) % skin % t = t(43) profiles(1) % skin % fastem(1) = 3.0 profiles(1) % skin % fastem(2) = 5.0 profiles(1) % skin % fastem(3) = 15.0 profiles(1) % skin % fastem(4) = 0.1 profiles(1) % skin % fastem(5) = 0.3 profiles(1) % ctp = 500. profiles(1) % cfraction = 0. !------------------------------------------------------------------------------------------------ ! This is used for profile data supplied by RTTOV authors ! ! OPEN (IUA,FILE='prof.dat',status='old') ! Read profile ONE and fill other profiles with profile one ! Then read rest of input profile one (NB gases in HITRAN order) ! nformat=10 ! how many fields in one line ! t ------------------------------------------------------------ ! JJM = 0 ! DO IRA = 1 , 1+(profiles(1)% nlevels-1)/nformat ! JJ = 1+JJM ! JJM = MIN(JJ+nformat-1,profiles(1)% nlevels) ! READ(IUA,*) (profiles(1) % t(J),J=JJ,JJM) ! end do ! q ------------------------------------------------------------ ! JJM = 0 ! DO IRA = 1 , 1+(profiles(1)% nlevels-1)/nformat ! JJ = 1+JJM ! JJM = MIN(JJ+nformat-1,profiles(1)% nlevels) ! READ(IUA,*) (profiles(1) % q(J),J=JJ,JJM) ! end do ! co2 ---------------------------------------------------------- ! if(coef(no_id)%nco2 > 0)then ! JJM = 0 ! DO IRA = 1 , 1+(profiles(1)% nlevels-1)/nformat ! JJ = 1+JJM ! JJM = MIN(JJ+nformat-1,profiles(1)% nlevels) ! READ(IUA,*) (profiles(1) % co2(J),J=JJ,JJM) ! end do ! profiles(1) % co2_data = .true. ! else ! profiles(1) % co2_data = .false. ! endif ! o3 ----------------------------------------------------------- ! JJM = 0 ! DO IRA = 1 , 1+(profiles(1)% nlevels-1)/nformat ! JJ = 1+JJM ! JJM = MIN(JJ+nformat-1,profiles(1)% nlevels) ! READ(IUA,*) (profiles(1) % o3(J),J=JJ,JJM) ! end do ! if (profiles(1) % o3(1) > 0._jprb .and. coef(no_id)%nozone > 0 )Then ! profiles(1) % ozone_data = .true. ! else ! profiles(1) % ozone_data = .false. ! endif ! n2o ----------------------------------------------------------- ! if(coef(no_id)%nn2o > 0)then ! JJM = 0 ! DO IRA = 1 , 1+(profiles(1)% nlevels-1)/nformat ! JJ = 1+JJM ! JJM = MIN(JJ+nformat-1,profiles(1)% nlevels) ! READ(IUA,*) (profiles(1) % n2o(J),J=JJ,JJM) ! end do ! profiles(1) % n2o_data = .true. ! else ! profiles(1) % n2o_data = .false. ! endif ! ! co ------------------------------------------------------------ ! if(coef(no_id)%nco > 0)then ! JJM = 0 ! DO IRA = 1 , 1+(profiles(1)% nlevels-1)/nformat ! JJ = 1+JJM ! JJM = MIN(JJ+nformat-1,profiles(1)% nlevels) ! READ(IUA,*) (profiles(1) % co(J),J=JJ,JJM) ! end do ! profiles(1) % co_data = .true. ! else ! profiles(1) % co_data = .false. ! endif ! ch4 ------------------------------------------------------------ ! if(coef(no_id)%nch4 > 0)then ! JJM = 0 ! DO IRA = 1 , 1+(profiles(1)% nlevels-1)/nformat ! JJ = 1+JJM ! JJM = MIN(JJ+nformat-1,profiles(1)% nlevels) ! READ(IUA,*) (profiles(1) % ch4(J),J=JJ,JJM) ! end do ! profiles(1) % ch4_data = .true. ! else ! profiles(1) % ch4_data = .false. ! endif ! clw ---------------------------------------------------------- ! JJM = 0 ! DO IRA = 1 , 1+(profiles(1)% nlevels-1)/nformat ! JJ = 1+JJM ! JJM = MIN(JJ+nformat-1,profiles(1)% nlevels) ! READ(IUA,*) (profiles(1) % clw(J),J=JJ,JJM) ! end do ! check value of first level ! profiles(1) % clw_data = profiles(1) % clw(1) >= 0.0_JPRB ! other -------------------------------------------------------- ! READ(IUA,*) profiles(1) % s2m % t ,& ! & profiles(1) % s2m % q ,& ! & profiles(1) % s2m % p ,& ! & profiles(1) % s2m % u ,& ! & profiles(1) % s2m % v, & ! & profiles(1) % s2m % wfetc ! profiles(1) % s2m % o = 0._JPRB ! ! READ(IUA,*) profiles(1) % skin % t ,& ! & profiles(1) % skin % fastem ! ! READ(IUA,*) profiles(1) % ctp,& ! & profiles(1) % cfraction ! ! CLOSE(IUA) ! ! WRITE(IOOUT,*)' INPUT PROFILE' do jj = 1,nlevels write(ioout,445) profiles(1) % p(jj), p(jj), profiles(1) % t(jj),profiles(1) % q(jj), & profiles(1) % o3(jj) enddo write(ioout,*) ' ' WRITE(IOOUT,445) profiles(1) % s2m % p ,& & profiles(1) % s2m % t ,& & profiles(1) % s2m % q ,& & profiles(1) % s2m % u ,& & profiles(1) % s2m % v WRITE(IOOUT,445) profiles(1) % skin % t ,& & profiles(1) % skin % fastem WRITE(IOOUT,445) profiles(1) % ctp,& & profiles(1) % cfraction WRITE(IOOUT,*)' ' ! Output information for original input file from RTTOV authors ! WRITE(IOOUT,*)' INPUT PROFILE' ! WRITE(IOOUT,444) (profiles(1) % t(JJ) ,JJ=1,profiles(1)% nlevels) ! WRITE(IOOUT,444) (profiles(1) % q(JJ) ,JJ=1,profiles(1)% nlevels) ! WRITE(IOOUT,444) (profiles(1) % o3(JJ) ,JJ=1,profiles(1)% nlevels) ! if(coef(no_id)%nco2> 0)then ! WRITE(IOOUT,444) (profiles(1) % co2(JJ) ,JJ=1,profiles(1)% nlevels) ! else ! profiles(1) % co2(:) =0.0 ! endif ! if(coef(no_id)%nn2o> 0)then ! WRITE(IOOUT,444) (profiles(1) % n2o(JJ) ,JJ=1,profiles(1)% nlevels) ! else ! profiles(1) % n2o(:) =0.0 ! endif ! if(coef(no_id)%nco> 0)then ! WRITE(IOOUT,444) (profiles(1) % co(JJ) ,JJ=1,profiles(1)% nlevels) ! else ! profiles(1) % co(:) =0.0 ! endif ! if(coef(no_id)%nch4> 0)then ! WRITE(IOOUT,444) (profiles(1) % ch4(JJ) ,JJ=1,profiles(1)% nlevels) ! else ! profiles(1) % ch4(:) =0.0 ! endif ! WRITE(IOOUT,444) (profiles(1) % clw(JJ) ,JJ=1,profiles(1)% nlevels) ! WRITE(IOOUT,444) profiles(1) % s2m % t ,& ! & profiles(1) % s2m % q ,& ! & profiles(1) % s2m % p ,& ! & profiles(1) % s2m % u ,& ! & profiles(1) % s2m % v ,& ! & profiles(1) % s2m % wfetc ! WRITE(IOOUT,444) profiles(1) % skin % t ,& ! & profiles(1) % skin % fastem ! WRITE(IOOUT,444) profiles(1) % ctp,& ! & profiles(1) % cfraction ! WRITE(IOOUT,*)' ' ! Convert lnq to q in ppmv for profile ! ! profiles(1) % q(:) = (exp(profiles(1) % q(:)) / 1000._JPRB) * q_mixratio_to_ppmv ! profiles(1) % s2m % q = (exp(profiles(1) % s2m % q) / 1000._JPRB) * q_mixratio_to_ppmv ! Keep Ozone in ppmv ! viewing geometry profiles(1) % zenangle = ZENANG profiles(1) % azangle = AZANG profiles(1) % latitude = LAT profiles(1) % elevation = ZERHT profiles(1) % sunzenangle = SUNZANG profiles(1) % sunazangle = SUNAZANG profiles(1) % addsolar = solrad profiles(1) % addrefrac = refrac ! surface type profiles(1) % skin % surftype = nsurf(no_id) profiles(1) % skin % watertype = nwater(no_id) !.. Fill profile arrays with the 1 profile NPROF times DO J = 1 , NPROF(no_id) profiles(j) % p(:) = profiles(1) % p(:) profiles(j) % t(:) = profiles(1) % t(:) profiles(j) % q(:) = profiles(1) % q(:) profiles(j) % o3(:) = profiles(1) % o3(:) profiles(j) % co2(:) = profiles(1) % co2(:) profiles(j) % n2o(:) = profiles(1) % n2o(:) profiles(j) % co(:) = profiles(1) % co(:) profiles(j) % ch4(:) = profiles(1) % ch4(:) if( laerosl ) then profiles(j) % aerosols(:,:)=profiles(1) % aerosols(:,:) endif if( lclouds ) then profiles(j) % cloud(:,:)=profiles(1) % cloud(:,:) profiles(j) % cfrac(:,:)=profiles(1) % cfrac(:,:) endif profiles(j) % clw(:) = profiles(1) % clw(:) profiles(j) % idg = profiles(1) % idg profiles(j) % ish = profiles(1) % ish profiles(j) % s2m = profiles(1) % s2m profiles(j) % skin = profiles(1) % skin profiles(j) % ctp = profiles(1) % ctp profiles(j) % cfraction = profiles(1) % cfraction profiles(j) % ozone_data = profiles(1) % ozone_data profiles(j) % co2_data = profiles(1) % co2_data profiles(j) % n2o_data = profiles(1) % n2o_data profiles(j) % co_data = profiles(1) % co_data profiles(j) % ch4_data = profiles(1) % ch4_data profiles(j) % aer_data = profiles(1) % aer_data profiles(j) % cld_data = profiles(1) % cld_data profiles(j) % clw_data = profiles(1) % clw_data profiles(j) % zenangle = profiles(1) % zenangle profiles(j) % azangle = profiles(1) % azangle profiles(j) % sunzenangle = profiles(1) % sunzenangle profiles(j) % sunazangle = profiles(1) % sunazangle profiles(j) % latitude = profiles(1) % latitude profiles(j) % elevation = profiles(1) % elevation profiles(j) % addsolar = profiles(1) % addsolar profiles(j) % addrefrac = profiles(1) % addrefrac end do ! Build the list of channels/profiles indices emissivity(:) = 0.0_JPRB channels(:) = 0_jpim lprofiles(:) = 0_jpim nch = 0_jpim Do j = 1 , nprof(no_id) nchan_byprof(j) = nchan(j,no_id) DO jch = 1,nchan1(no_id) nch = nch +1_jpim lprofiles ( nch ) = j if (all_channels)then channels( nch ) = ichan(jch,no_id) else channels( nch ) = jch endif emissivity( nch ) = surfem(jch,no_id) End Do End Do ! !! write(6,*)'channels=',(channels(ich2),ich2=1,nchannels(no_id)) ! save input values of emissivities for all calculations ! calculate emissivity where the input emissivity value is less than 0.01 input_emissivity(:) = emissivity(:) calcemis(:) = emissivity(:) < 0.01_JPRB WRITE(IOOUT,*)' NUMBER OF PROFILES PROCESSED=',NPROF(no_id) WRITE(IOOUT,*)' ' ! WRITE(IOOUT,*)'CHANNELS PROCESSED:' WRITE(IOOUT,111) (ichan(J,no_id), J = 1,NCHAN1(no_id)) WRITE (IOOUT,*)' ' WRITE(IOOUT,*)'INPUT SURFACE EMISSIVITIES '& & ,'SAT =', instrument(2,no_id) WRITE(IOOUT,444) (emissivity(J),J=1,NCHAN1(no_id)) WRITE(IOOUT,*)' ' ! PERFORM RADIATIVE TRANSFER CALCULATIONS call rttov_direct( & & rttov_errorstatus, &! out & nprof(no_id), &! in & nchannels(no_id),&! in & channels, &! in & lprofiles, &! in & addinterp, &! in & profiles, &! in & coef(no_id), &! in & coef_scatt_ir(no_id) , & & optp(no_id) , & & lsun, &! in & laerosl, &! in & lclouds, &! in & calcemis, &! in & emissivity, &! inout & transmission, &! out & radiance ) ! inout If ( any( rttov_errorstatus(:) == errorstatus_warning ) ) Then Do j = 1, nprof(no_id) If ( rttov_errorstatus(j) == errorstatus_warning ) Then write ( ioout, * ) 'rttov warning for profile',j End If End Do End If If ( any( rttov_errorstatus(:) == errorstatus_fatal ) ) Then Do j = 1, nprof(no_id) If ( rttov_errorstatus(j) == errorstatus_fatal ) Then write ( ioout, * ) 'rttov error for profile',j End If End Do Stop End If ! transfer data to printing arrays pr_radcld(:) = 0.0_JPRB pr_trans(:) = 0.0_JPRB pr_emis(:) = 0.0_JPRB pr_trans_lev(:,:) = 0.0_JPRB pr_upclr(:) = 0.0_JPRB pr_ovcst(:,:) = 0.0_JPRB pr_up(:,:) = 0.0_JPRB pr_down(:,:) = 0.0_JPRB pr_surf(:,:) = 0.0_JPRB ! do j = 1 , nchannels(no_id) pr_radcld(j) = radiance % cloudy(j) pr_trans(j) = Transmission % tau_total(J) pr_emis(j) = emissivity(j) pr_upclr(j) = radiance % upclear(J) do ilev = 1 , nlevels pr_trans_lev(ilev,j) = Transmission % tau_layers(ilev,J) pr_ovcst(ilev,j) = radiance % overcast(ILEV,J) pr_up(ilev,j) = radiance % up(ILEV,J) pr_down(ilev,j) = radiance % down(ILEV,J) pr_surf(ilev,j) = radiance % surf(ILEV,J) enddo enddo ! OUTPUT RESULTS ! NPRINT = 1+ INT((nchan_byprof(1)-1)/10) DO JN=1,NPROF(no_id) WRITE(IOOUT,*)' -----------------' WRITE(IOOUT,*)' Profile number ',JN, 'Instrument ',& & instrument(3,no_id) WRITE(IOOUT,*)' -----------------' WRITE(IOOUT,*)' ' ! JOFF=NCHAN(no_id)*(JN-1) JOFF1=nchannels(no_id)/nprof(no_id)*(JN-1) JOFF2=nchannels(no_id)/nprof(no_id)*(JN-1) JOFF3=nchannels(no_id)/nprof(no_id)*(JN-1) if(isun==1)then WRITE(IOOUT,777) profiles(jn)%zenangle,profiles(jn)%azangle, & & csun,profiles(jn)%sunzenangle,profiles(jn)%sunazangle,profiles(jn)%skin%surftype,& & profiles(jn)%skin%watertype,profiles(jn)%latitude,profiles(jn)%elevation,cref,caer,ccld,& & instrument(2,no_id) else WRITE(IOOUT,776)profiles(jn)%zenangle,profiles(jn)%azangle,csun,profiles(jn)%skin%surftype,& & profiles(jn)%skin%watertype,profiles(jn)%latitude,profiles(jn)%elevation,cref,caer,ccld,& & instrument(2,no_id) endif WRITE(IOOUT,222) (radiance % bt(J+JOFF1),J=1,nchan1(no_id)) WRITE(IOOUT,*)' ' WRITE(IOOUT,*)'CALCULATED RADIANCES: SAT =', instrument(2,no_id) WRITE(IOOUT,222) (radiance % total(J+JOFF1),J=1,nchan1(no_id)) WRITE(IOOUT,*)' ' WRITE(IOOUT,*)'CALCULATED OVERCAST RADIANCES: SAT =', instrument(2,no_id) WRITE(IOOUT,222) (pr_radcld(J+JOFF2),J=1,nchan1(no_id)) WRITE (IOOUT,*)' ' WRITE(IOOUT,*)'CALCULATED SURFACE TO SPACE TRANSMITTANCE: S'& & ,'AT =',instrument(2,no_id) WRITE(IOOUT,4444) (pr_trans(J+JOFF2),J=1,nchan1(no_id)) WRITE (IOOUT,*)' ' WRITE(IOOUT,*)'CALCULATED SURFACE EMISSIVITIES '& & ,'SAT =',instrument(2,no_id) WRITE(IOOUT,444) (pr_emis(J+JOFF2),J=1,nchan1(no_id)) ! IF(JN.EQ.1 .AND. nchan(1,no_id) .LE. 20)THEN DO NP = 1 , NPRINT WRITE (IOOUT,*)' ' WRITE (IOOUT,*)'Level to space transmittances for channels' ! WRITE(IOOUT,1115)(pr_pol(j),& WRITE(IOOUT,1115) (ICHAN(J,no_id),& & J = 1+(NP-1)*10,MIN(10+(NP-1)*10,nchan1(no_id))) DO ILEV = 1 , profiles(1) % nlevels WRITE(IOOUT,4445)ILEV,(pr_trans_lev(ilev,J+JOFF2),& & J = 1+(NP-1)*10,MIN(10+(NP-1)*10,nchan1(no_id))) end do WRITE(IOOUT,1115) (ICHAN(J,no_id),& & J = 1+(NP-1)*10,MIN(10+(NP-1)*10,nchan1(no_id))) end do ENDIF ! ! Print overcast radiance upwelling arrays IF(JN==1 .AND. IFULL(no_id)==2 .AND. nchan(1,no_id)<=20)THEN DO NP = 1 , NPRINT WRITE (IOOUT,*)' ' WRITE (IOOUT,*)'Level to space upwelling overcast radiances for channels' WRITE(IOOUT,1115) (ICHAN(J,no_id),& & J = 1+(NP-1)*10,MIN(10+(NP-1)*10,nchan1(no_id))) DO ILEV = 1 , profiles(1) % nlevels WRITE(IOOUT,4446)ILEV,(pr_ovcst(ILEV,J+JOFF2),& & J = 1+(NP-1)*10,MIN(10+(NP-1)*10,nchan1(no_id))) end do WRITE(IOOUT,1115) (ICHAN(J,no_id),& & J = 1+(NP-1)*10,MIN(10+(NP-1)*10,nchan1(no_id))) end do ENDIF IF(JN==1 .AND. IFULL(no_id)==2 .AND. nchan(1,no_id)<=20)THEN DO NP = 1 , NPRINT WRITE (IOOUT,*)' ' WRITE (IOOUT,*)'Above cloud upwelling radiances for channels' WRITE(IOOUT,1115) (ICHAN(J,no_id),& & J = 1+(NP-1)*10,MIN(10+(NP-1)*10,nchan1(no_id))) DO ILEV = 1 , NLEVELS WRITE(IOOUT,4446)ILEV,(pr_up(ILEV,J+JOFF2),& & J = 1+(NP-1)*10,MIN(10+(NP-1)*10,nchan1(no_id))) end do WRITE(IOOUT,1115) (ICHAN(J,no_id),& & J = 1+(NP-1)*10,MIN(10+(NP-1)*10,nchan1(no_id))) end do ! DO NP = 1 , NPRINT WRITE (IOOUT,*)' ' WRITE (IOOUT,*)'Downwelling radiances for channels' WRITE(IOOUT,1115) (ICHAN(J,no_id),& & J = 1+(NP-1)*10,MIN(10+(NP-1)*10,nchan1(no_id))) DO ILEV = 1 , NLEVELS WRITE(IOOUT,4446)ILEV,(pr_down(ILEV,J+JOFF2),& & J = 1+(NP-1)*10,MIN(10+(NP-1)*10,nchan1(no_id))) end do WRITE(IOOUT,1115) (ICHAN(J,no_id),& & J = 1+(NP-1)*10,MIN(10+(NP-1)*10,nchan1(no_id))) end do DO NP = 1 , NPRINT WRITE (IOOUT,*)' ' WRITE (IOOUT,*)'Radiances from cloud for channels' WRITE(IOOUT,1115) (ICHAN(J,no_id),& & J = 1+(NP-1)*10,MIN(10+(NP-1)*10,nchan1(no_id))) DO ILEV = 1 , NLEVELS WRITE(IOOUT,4446)ILEV,(pr_surf(ILEV,J+JOFF2),& & J = 1+(NP-1)*10,MIN(10+(NP-1)*10,nchan1(no_id))) end do WRITE(IOOUT,1115) (ICHAN(J,no_id),& & J = 1+(NP-1)*10,MIN(10+(NP-1)*10,nchan1(no_id))) end do Endif End do ! End of profile loop WRITE(*,*) ' FORWARD MODEL FINISHED' ! deallocate transmittances Deallocate( transmission % tau_total ,stat= alloc_status(7)) Deallocate( transmission % tau_layers ,stat= alloc_status(8)) ! IF (IFULL(no_id).GE.1)THEN ! allocate arrays used for TL/AD/K tests allocate(xkbav(nlevels,jpnav,nchannels(no_id)),stat= alloc_status(1)) allocate(xkradovu(nlevels,nchannels(no_id)),stat= alloc_status(2)) allocate(xkbsav(sscvar,nchannels(no_id)),stat= alloc_status(6)) allocate(xkbem(nchannels(no_id)),stat= alloc_status(7)) If( any(alloc_status /= 0) ) then errorstatus = errorstatus_fatal Write( errMessage, '( "mem deallocation error")' ) Call Rttov_ErrorReport (errorstatus, errMessage, NameOfRoutine) Stop End If xkbav(:,:,:)=0._JPRB ! !----------------------------------------------------------- ! Test tangent linear !----------------------------------------------------------- write(*,*) 'Tangent linear'! call TSTRAD_TL( & & nprof(no_id), &! in & nchannels(no_id),&! in & channels, &! in & lprofiles, &! in & addinterp, &! in & profiles, &! in & coef(no_id), &! in & coef_scatt_ir(no_id),& & optp(no_id), & & lsun, & & lclouds, &! in & laerosl, & & calcemis, &! in & input_emissivity) ! in write(*,*) 'Adjoint'! call TSTRAD_AD( & & nprof(no_id), &! in & nchannels(no_id),&! in & channels, &! in & lprofiles, &! in & addinterp, &! in & profiles, &! in & coef(no_id), &! in & coef_scatt_ir(no_id),& & optp(no_id), & & lsun, & & lclouds, &! in & laerosl, & & calcemis, &! in & input_emissivity) ! in write(*,*) 'K' call TSTRAD_K( & & nprof(no_id), &! in & nchannels(no_id),&! in & channels, &! in & lprofiles, &! in & addinterp, &! in & profiles, &! in & coef(no_id), &! in & coef_scatt_ir(no_id),& & optp(no_id), & & lsun, & & lclouds, &! in & laerosl, & & calcemis, &! in & input_emissivity) ! in ! Deallocate test arrays deallocate(xkbav ,stat= alloc_status(1)) deallocate(xkradovu,stat= alloc_status(2)) deallocate(xkbsav ,stat= alloc_status(6)) deallocate(xkbem ,stat= alloc_status(7)) If( any(alloc_status /= 0) ) then errorstatus = errorstatus_fatal Write( errMessage, '( "mem deallocation error for xkb arrays")' ) Call Rttov_ErrorReport (errorstatus, errMessage, NameOfRoutine) Stop End If ENDIF ! End of TL/AD/K Block do j = 1, nprof(no_id) ! deallocate model profiles atmospheric arrays deallocate( profiles(j) % p ,stat=alloc_status(1)) deallocate( profiles(j) % t ,stat=alloc_status(2)) deallocate( profiles(j) % q ,stat=alloc_status(3)) deallocate( profiles(j) % o3 ,stat=alloc_status(4)) deallocate( profiles(j) % clw ,stat=alloc_status(5)) If( any(alloc_status /= 0) ) then errorstatus = errorstatus_fatal Write( errMessage, '( "mem deallocation error")' ) Call Rttov_ErrorReport (errorstatus, errMessage, NameOfRoutine) Stop End If end do deallocate( profiles,stat=alloc_status(1)) ! number of channels per RTTOV call is only nchannels deallocate( channels ,stat=alloc_status(1)) deallocate( nchan_byprof ,stat=alloc_status(2)) deallocate( lprofiles ,stat=alloc_status(3)) deallocate( emissivity ,stat=alloc_status(4)) deallocate( fresnrefl ,stat=alloc_status(5)) deallocate( calcemis ,stat=alloc_status(6)) deallocate( input_emissivity ,stat= alloc_status(14)) If( any(alloc_status /= 0) ) then errorstatus = errorstatus_fatal Write( errMessage, '( "mem deallocation error for channels etc")' ) Call Rttov_ErrorReport (errorstatus, errMessage, NameOfRoutine) Stop End If deallocate(pr_radcld ,stat= alloc_status(31)) deallocate(pr_trans ,stat= alloc_status(32)) deallocate(pr_emis ,stat= alloc_status(33)) deallocate(pr_trans_lev ,stat= alloc_status(34)) deallocate(pr_upclr ,stat= alloc_status(37)) deallocate(pr_ovcst ,stat= alloc_status(40)) If( any(alloc_status /= 0) ) then errorstatus = errorstatus_fatal Write( errMessage, '( "mem deallocation error for printing arrays")' ) Call Rttov_ErrorReport (errorstatus, errMessage, NameOfRoutine) Stop End If asw = 0 ! deallocate radiance arrays call rttov_alloc_rad (errorstatus,nchan1(no_id),radiance,profiles(1) % nlevels,asw) If(errorstatus /= errorstatus_success) Then Write( errMessage, '( "deallocation error for radiances")' ) Call Rttov_ErrorReport (errorstatus, errMessage, NameOfRoutine) Endif asw = 0 ! deallocate profile arrays ! Memory error when using this deallocation routine ! ! call rttov_alloc_prof (errorstatus,nprof(no_id),profiles,nlevels,coef_scatt_ir(no_id),asw, & ! & addclouds = lclouds, addaerosl = laerosl ) ! ! deallocate( profiles,stat=alloc_status(1)) If( any(alloc_status /= 0) ) then errorstatus = errorstatus_fatal Write( errMessage, '( "mem deallocation error for profiles")' ) Call Rttov_ErrorReport (errorstatus, errMessage, NameOfRoutine) Stop End If ENDDO ! End of sensor id loop ! ! Deallocate coeffs ! Deallocate (nchan1,stat= alloc_status(1)) Do no_id = 1, nrttovid Call rttov_dealloc_coef (errorstatus, coef(no_id),coef_scatt_ir(no_id),optp(no_id)) If(errorstatus /= errorstatus_success) Then Write( errMessage, '( "deallocation error for coeffs")' ) Call Rttov_ErrorReport (errorstatus, errMessage, NameOfRoutine) Endif End Do IF (LHOOK) CALL DR_HOOK('TSTRAD',1,ZHOOK_HANDLE) 111 FORMAT(1X,10I8) 1115 FORMAT(3X,10I8) 2222 FORMAT(1X,10(1x,F8.6)) 222 FORMAT(1X,10F8.2) 333 FORMAT(1X,I3,20I5) 3333 FORMAT(1X,I3,2I5) 444 FORMAT(1X,10F8.3) 445 FORMAT(1X,5F12.4) 4444 FORMAT(1X,10F8.4) !4445 FORMAT(1X,I2,10F8.4) 4445 FORMAT(1X,I2,10F13.9) 4446 FORMAT(1X,I2,10F8.3) 555 FORMAT(1X,10E8.2) 776 FORMAT( & & ' ZENITH ANGLE =',F7.2,/ & & ' AZIMUTH ANGLE =',F7.2,/& & ' SOLAR RADIATION =',A7,/& & ' SURFACE TYPE =',I7,/& & ' WATER TYPE =',I7,/ & & ' LATITUDE =',F7.2,/& & ' ELEVATION =',F7.2/& & ' REFRACTION =',A7,/& & ' AEROSOLS =',A7,/& & ' CLOUDS =',A7,//& & ' CALCULATED BRIGHTNESS TEMPERATURES: SAT =',I2) 777 FORMAT( & & ' ZENITH ANGLE =',F7.2,/ & & ' AZIMUTH ANGLE =',F7.2,/& & ' SOLAR RADIATION =',A7,/& & ' SOLAR ZENITH ANGLE =',F7.2,/& & ' SOLAR AZIMUTH ANGLE=',F7.2,/ & & ' SURFACE TYPE =',I7,/& & ' WATER TYPE =',I7,/ & & ' LATITUDE =',F7.2,/& & ' ELEVATION =',F7.2/& & ' REFRACTION =',A7,/& & ' AEROSOLS =',A7,/& & ' CLOUDS =',A7,//,& &'CALCULATED BRIGHTNESS TEMPERATURES: SAT =',I2 ) END PROGRAM ETL_RTTOV91