SUBROUTINE SFIT4( ) ! SFIT4 VERSION 003.82 2011 ! SEE FILE COMMENTS.SAVED FOR OLDER COMMENTS AND HISTORY USE PARAMS USE VIBFCN USE RETVPARAM USE FRWDMDL USE BANDPARAM USE MOLCPARAM USE DATAFILES USE SOLAR USE OPT USE DIAGNOSTIC USE SYNSPEC USE LINEPARAM USE INITIALIZE USE CHANNEL USE READIN USE WRITEOUT USE BINPUT_4_0 USE RAYTRACE IMPLICIT NONE LOGICAL :: CONVERGE=.FALSE., DEBUG=.FALSE., DIVWARN=.FALSE., FILOPEN=.FALSE. CHARACTER, DIMENSION(5) :: CPNAM*14 INTEGER :: NFIT = 0, NLEV = 0, NVAR = 0, NAERR = 0 INTEGER :: NMASS, KZMAXLAY, ISMIX, INDXX, KVERT, ITER=0, NEGFLAG, NEGLOOP INTEGER :: I, IA, IB, IBAND, IND, IOUT, J, JROW, JCOL, JSCAN, K, KK, N INTEGER, DIMENSION(MAXSPE) :: NLAY REAL(DOUBLE) :: DELZ = 0.0D0 REAL(DOUBLE) :: RHO = 0.0D0 REAL(DOUBLE) :: YMAX = 0.0D0 REAL(DOUBLE) :: SERR = 0.0D0 REAL(DOUBLE) :: SIGB = 0.0D0 REAL(DOUBLE) :: SIGC = 0.0D0 REAL(DOUBLE) :: WAVE = 0.0D0 REAL(DOUBLE) :: SIGZ = 0.0D0 REAL(DOUBLE) :: AIRCOL = 0.0D0 REAL(DOUBLE), DIMENSION(MMAX) :: YHAT = 0.0D0 REAL(DOUBLE), DIMENSION(NMAX) :: XHAT = 0.0D0 REAL(DOUBLE), DIMENSION(12) :: FX = 0.0D0 REAL(DOUBLE), DIMENSION(MOLMAX,LAYMAX) :: VERSUM = 0.0D0, VOSUM = 0.0D0 REAL(DOUBLE), DIMENSION(MMAX) :: Y = 0.0D0 REAL(DOUBLE), DIMENSION(NMAX) :: PARM = 0.0D0 REAL(DOUBLE), DIMENSION(NMAX) :: XAPR = 0.0D0 REAL(DOUBLE), DIMENSION(NMAX) :: SPARM = 0.0D0 REAL(DOUBLE), DIMENSION(:), ALLOCATABLE :: SED REAL(DOUBLE), DIMENSION(:,:), ALLOCATABLE :: SAOUT ! ------------------------------------------------------------------------------ CALL FILEOPEN( 18, 2 ) CALL FILEOPEN( 20, 2 ) CALL FILEOPEN( 31, 2 ) CALL FILEOPEN( 73, 2 ) WRITE (31, *) VERSION ! --- READ IN RETRIEVAL LAYERING FROM STATION.LAYERS FILE CALL READLAYRS( NLEV ) KMAX = NLEV ! --- PRINT INVERSION PARAMETERS NOW THAT WE HAVE NLEV ! ISOTOPE CALL READCK1( NLEV, NEGFLAG ) ! --- PRINT OUT FM, RT PARAMETERS CALL READCK2( CPNAM ) ! --- CHECK WE ONLY FIT PHASE AND MODULATION FUNCTIONS TYPE = 2 IF( IEAP.NE.2 .AND. IRTEAP.EQ.1 ) GOTO 667 IF( IEPHS.NE.2 .AND. IRTEPHS.EQ.1 ) GOTO 668 ! --- CHECK THAT BACKGROUND AND SHIFT PARMS ARE IN RANGE IF( NBACK.LT.1 .OR. NBACK.GT.3 ) GOTO 664 IF( ISPARM.LT.0 .OR. ISPARM.GT.3 ) GOTO 665 ! --- PRINT OUT BAND PARAMETERS CALL READCK3( ) ! --- READ IN SPECTRA CALL SETUP1( ) ! --- PERFORM RAYTRACE FOR ALL SZA (NSPEC) CALL LBLATM( 0, NLEV ) IF( RAYTONLY ) STOP ': RAYTRACE ONLY FLAG SET.' IF( NLEV .LT. KMAX ) KMAX = NLEV KVERT = NSPEC +1 ! --- FIND MAXIMUM MASS PATH FOR EACH LAYER DO KK = 1, KMAX PMASMX(KK) = 0.D0 PMASMX(KK) = DMAX1(MAXVAL(CCC(:NSPEC,KK)),PMASMX(KK)) !print *, PMASMX(KK) END DO ! --- MORE SETUP GET HITRAN, CHECK MODULATION EFF & PHASE, SOLAR SPECTRUM... CALL SETUP2 ! --- FINISH SETUP CALCULATE CROSS SECTIONS, -1 means, crosssections for all levels. CALL SETUP3( .TRUE., -1 ) ! --- IF EMISSION, SNR IS THE NOISE ON THE MEASURMENT GIVEN IN (MW/(CM^2*SR*CM-1)) MP ! MUSTFIX ! COMMENTED OUT JWH JUN 2012 ! check initialize:filse ! DO I=1, NBAND ! IF (IEMISSION/=1.AND.(.NOT.IENORM(i)/=1)) SNR = 1/SNR ! ENDDO ! DO K = 1, NSPEC +1 ! print*, astang(k), appang(k), ispec(k), rearth(k), reflat(k), xvb(k) ! enddo ! do j=1, nband ! do k=1, NSCAN(j) ! print*, j, k, ISCAN(j,k), appang(ISCAN(j,k)) ! enddo ! enddo ! --- THIS LOOP DOES NOT DO MUCH BUT REMINDS OF WHERE TO SETUP VARIABLE LAYER SCHEMES ! --- KZTAN STILL USED DO K = 1, NSPEC NLAY(K) = KMAX KZTAN(K) = NLAY(K) ! # LAYERS IN THIS MASSPATH KZMAXLAY = KZTAN(K) END DO ! --- INITIALIZE PARAMETERS FOR RETRIEVAL NVAR = 0 ! --- BACKGROUND FITTING IF (NBACK == 1) THEN ELSE IF (NBACK == 2) THEN IF (NFITS > 0) THEN ! --- BACKGROUND SLOPE - NBACK=2 PNAME(:NFITS) = 'BckGrdSlp' PARM(:NFITS) = BCKSL SPARM(:NFITS) = SBCKSL ! --- BACKGROUND CURVATURE - NBACK=3 NVAR = NFITS ENDIF ELSE IF (NBACK == 3) THEN IF (NFITS > 0) THEN ! --- BACKGROUND SLOPE - NBACK=2 PNAME(NVAR+1:NFITS*2-1+NVAR:2) = 'BckGrdSlp' PARM(NVAR+1:NFITS*2-1+NVAR:2) = BCKSL SPARM(NVAR+1:NFITS*2-1+NVAR:2) = SBCKSL ! --- BACKGROUND CURVATURE - NBACK=3 PNAME(NVAR+2:NFITS*2+NVAR:2) = 'BckGrdCur' PARM(NVAR+2:NFITS*2+NVAR:2) = BCKCRV SPARM(NVAR+2:NFITS*2+NVAR:2) = SBCKCRV NVAR = NFITS*2 + NVAR ENDIF ENDIF ENDIF ENDIF NBKFIT = (NBACK - 1)*NFITS ! --- NO WAVENUMBER SHIFT IF (ISPARM /= 0) THEN ! --- SINGLE WAVELENGTH SHIFT FOR ALL BANDPASSES IF (ISPARM <= 1) THEN NVAR = NVAR + 1 PNAME(NVAR) = 'SWNumShft' PARM(NVAR) = WSHFT SPARM(NVAR) = SWSHFT NSHIFT = 1 ELSE ! --- INDEPENDENT WAVELENGTH SHIFT FOR EACH BANDPASS (ISPARM=2) ! --- OR INDEPENDENT WAVELENGTH SHIFT FOR EACH FITTING (ISPARM=3) IF (ISPARM == 2) N = NBAND IF (ISPARM == 3) N = NFITS IF (N > 0) THEN PNAME(NVAR+1:N+NVAR) = 'IWNumShft' PARM(NVAR+1:N+NVAR) = WSHFT SPARM(NVAR+1:N+NVAR) = SWSHFT NVAR = N + NVAR ENDIF NSHIFT = N ENDIF ENDIF ! --- FIT ZERO LEVELS IF IZERO=1 ! --- TOTAL NUMBER OF ZERO LEVEL FITS=NZERO NZERO = 0 DO I = 1, NBAND IF (IZERO(I) .NE. 1 ) CYCLE N = NSCAN(I) IF (N > 0) THEN PNAME(NVAR+1:N+NVAR) = 'ZeroLev' PARM(NVAR+1:N+NVAR) = ZSHIFT(I,1) SPARM(NVAR+1:N+NVAR) = SZERO(I) NVAR = N + NVAR NZERO = N + NZERO ENDIF END DO ! --- SOLAR LINES INCLUSION NSOLAR = 0 IF( IFCO )THEN DO I = 1, 5 IF ( .NOT. ICFIX(I) ) CYCLE NVAR = NVAR + 1 NSOLAR = NSOLAR + 1 IF (NSOLAR == 1) NSOLAR1 = NVAR SELECT CASE (I) CASE (1) PNAME(NVAR) = CPNAM(1) CASE (2) PNAME(NVAR) = CPNAM(2) CASE (3) PNAME(NVAR) = CPNAM(3) CASE (4) PNAME(NVAR) = CPNAM(4) CASE (5) PNAME(NVAR) = CPNAM(5) END SELECT PARM(NVAR) = 1.0d0 !COPAR +1.0d0 SPARM(NVAR) = SCOPAR END DO ENDIF ! --- EMPIRICAL APODIZATION IF (IRTEAP /= 0) THEN NEAPRT = NEAP IF (NEAPRT > 0) THEN EAPF0(:NEAPRT) = EAPF(:NEAPRT) PNAME(NVAR+1:NEAPRT+NVAR) = 'EmpApdFcn' PARM(NVAR+1:NEAPRT+NVAR) = EAPPAR SPARM(NVAR+1:NEAPRT+NVAR) = SEAPPAR NVAR = NEAPRT + NVAR ENDIF ENDIF ! --- EMPIRICAL PHASE FUNCTION IF (IRTEPHS /= 0) THEN NEPHSRT = NEPHS IF (NEPHSRT > 0) THEN EPHSF0(:NEPHSRT) = EPHSF(:NEPHSRT) PNAME(NVAR+1:NEPHSRT+NVAR) = 'EmpPhsFnc' PARM(NVAR+1:NEPHSRT+NVAR) = EPHSPAR SPARM(NVAR+1:NEPHSRT+NVAR) = SEPHSPAR NVAR = NEPHSRT + NVAR ENDIF ENDIF ! --- DIFFERENTIAL WAVENUMBER SHIFT FOR RETRIEVAL GASES NDIFF = 0 IF (IFDIFF .NE.0 .AND. NRET.GT.1) THEN PNAME(NVAR+1:NRET-1+NVAR) = 'DWNumShft' PARM(NVAR+1:NRET-1+NVAR) = WSHFT SPARM(NVAR+1:NRET-1+NVAR) = SWSHFT NDIFF = NRET - 1 NVAR = NRET - 1 + NVAR ENDIF ! --- SKIP IF NOT IFPHASE; FIT PHASE ERRORS IF IFPHASE=T ! --- TOTAL NUMBER OF PHASE ERROR FITS=NPHASE NPHASE = 0 IF( IFPHASE )THEN DO I = 1, NBAND N = NSCAN(I) IF (N > 0) THEN PNAME(NVAR+1:N+NVAR) = 'SPhsErr' PARM(NVAR+1:N+NVAR) = PHS SPARM(NVAR+1:N+NVAR) = SPHS NVAR = N + NVAR NPHASE = N + NPHASE ENDIF END DO ENDIF ! --- RETRIEVAL GAS MIXING RATIOS ! --- MIXING RATIO AT ISMIX +1 ISMIX = NVAR DO KK = 1, NRET NGIDX(KK,1,0) = NVAR + 1 IF( IFPRF(KK) )THEN ! --- RETRIEVING VERTICAL PROFILE N = NLAYERS PNAME(NVAR+1:N+NVAR) = NAME(IGAS(KK)) IF (ILOGRETRIEVAL(KK)/=0) THEN !MP ! WE WANT A LOGARITHMIC PROFILE PARM(NVAR+1:N+NVAR) = LOG(XORG(KK,:N)) ELSE PARM(NVAR+1:N+NVAR) = 1.D0 ENDIF SPARM(NVAR+1:N+NVAR) = SIG(:N,KK) NVAR = NVAR + N ELSE ! --- SCALING VERTICAL DISTRIBUTION NVAR = NVAR + 1 PNAME(NVAR) = NAME(IGAS(KK)) PARM(NVAR) = COLSF(KK) SPARM(NVAR) = SCOLSF(KK) ENDIF NGIDX(KK,2,0) = NVAR END DO ! --- TEMPERATURE RETRIEVAL IF( IFTEMP )THEN NTEMP = NLAYERS PNAME(NVAR+1:NVAR+NTEMP) = 'TEMPERAT' PARM(NVAR+1:NVAR+NTEMP) = 1.D0 SPARM(NVAR+1:NVAR+NTEMP) = TSIGMA(1:NTEMP) NTEMP1 = NVAR + 1 NVAR = NVAR + NTEMP ENDIF ! --- INSERT CHANNEL PARAMETERS INTO STATE VECTOR PARM() CALL INSERT_CHANNEL_PARMS (NVAR, PARM, PNAME, SPARM) ! --- TEST FOR OVERFLOWS IF (NVAR > NMAX) GO TO 556 NFIT = NATMOS IF (NFIT > MMAX) GO TO 557 ! --- SETUP ARRAYS FOR OPTIMAL ESTIMATION SUBROUTINE ! --- FILL MATRIX SED COVARIANCE OF MEASURED SPECTRUM - OFFDIAGNAL ELEMENTS ASSUMED TO BE ZERO ! --- ALLOCATE ALLOCATE( SE(NFIT), SED(NFIT), STAT=NAERR ) IF (NAERR /= 0) THEN WRITE (6, *) 'COULD NOT ALLOCATE SE ARRAY' WRITE (6, *) 'ERROR NUMBER = ', NAERR STOP 'SFIT ALLOCATION' ENDIF SE(:) = 0.0D0 SED(:) = 0.0D0 WRITE (*, *) ' INITIALIZING VARIANCE VECTOR...' CALL FILSE (SED, NFIT) WRITE (16, 3695) NFIT, NVAR WRITE (16, 3696) ICOUNT = 0 ALLOCATE( SA(NVAR,NVAR), SAINV(NVAR,NVAR), SHAT(NVAR,NVAR), STAT=NAERR ) IF (NAERR /= 0) THEN WRITE (6, *) 'COULD NOT ALLOCATE SA ARRAY' WRITE (6, *) 'ERROR NUMBER = ', NAERR STOP 'SFIT ALLOCATION' ENDIF SA(:,:) = 0.0D0 SHAT(:,:) = 0.0D0 SAINV(:,:) = 0.0D0 ! --- FILL DIAGONAL ELEMENTS OF SA DO I = 1, NVAR SA(I,I) = SPARM(I)*SPARM(I) ENDDO ! --- OFF DIAGONAL ELEMENTS OF SA MATRIX (A PRIORI COMPONENTS) INDXX = ISMIX DO KK = 1, NRET !print *, 'fill off diag ', kk, IFPRF(KK) N = 1 IF( IFPRF(KK) ) THEN N = NLAYERS SELECT CASE ( IFOFF(KK) ) CASE ( 1:3 ) ! --- FILL OFF DIAGONAL ELEMENTS OF SA DO I = 1, NLAYERS DO J = 1, NLAYERS IF (I == J) CYCLE JROW = I + INDXX JCOL = J + INDXX IF( ZBAR(I) < ZGMIN(KK) ) CYCLE IF( ZBAR(J) < ZGMIN(KK) ) CYCLE IF( ZBAR(I) > ZGMAX(KK) ) CYCLE IF( ZBAR(J) > ZGMAX(KK) ) CYCLE DELZ = ZBAR(I) - ZBAR(J) SELECT CASE ( IFOFF(KK) ) CASE (1) !gaussian RHO = (ALOGSQ*DELZ/ZWID(KK))**2 RHO = MIN( RHO, 90.0D0 ) SA(JROW,JCOL) = SPARM(JROW)*SPARM(JCOL)*EXP(-RHO) CASE (2) !exponential RHO = ABS(ALOGSQ*DELZ/ZWID(KK)) RHO = MIN( RHO, 90.0D0 ) !664 SA(JROW,JCOL) = SPARM(JROW)*SPARM(JCOL)*EXP(-RHO) CASE (3) STOP "IFOFF=3 not supported" END SELECT END DO END DO ! --- READ IN FULL COVARIANCE FROM FILE CASE ( 4 ) INQUIRE( UNIT=64, OPENED=FILOPEN ) IF ( .NOT. FILOPEN )CALL FILEOPEN( 64, 3 ) DO I = 1, NLAYERS READ( 64,* ) (SA( I+INDXX, J+INDXX), J=1, N) END DO CASE ( 0 ) PRINT *, ' PROFILE RETRIEVAL GAS: ', NAME(IGAS(KK)), ' NO OFF DIAGONAL VALUES SET.' END SELECT ENDIF INDXX = INDXX + N END DO CALL FILECLOSE( 64, 2 ) Y(:NFIT) = TOBS(:NFIT) XAPR = PARM ! --- WRITE OUT COVARIENCE MATRIX ! --- NOT! WRITE OUT IN REVERSE, REVERSE FORM IF (F_WRTSA) THEN CALL FILEOPEN( 61, 1 ) INDXX = ISMIX DO KK = 1, NRET N = 1 IF( IFPRF(KK) ) THEN N = NLAYERS IF ( .NOT. ALLOCATED( SAOUT )) ALLOCATE( SAOUT( N, N )) WRITE (61, '(3X,45A14)') (PNAME(I),I = INDXX+1, INDXX + N) DO I = 1, N DO J = 1, N JROW = I + INDXX JCOL = J + INDXX SELECT CASE ( TRIM( TFILE(61))) CASE ( 'Sa.norm' ) ! fractional form SAOUT(I,J) = SA(JROW,JCOL) CASE ( 'Sa.pcol' ) ! partial column [(molec/cm^2)^2] RHO = XORG(KK,I)*XORG(KK,J)*CCC(NMASS,I)*CCC(NMASS,J) SAOUT(I,J) = SA(JROW,JCOL) * RHO CASE ( 'Sa.vmr' ) ! vmr [ppm^2] RHO = XORG(KK,I)*XORG(KK,J)*1.0D12 SAOUT(I,J) = SA(JROW,JCOL) * RHO CASE DEFAULT SAOUT(I,J) = SA(JROW,JCOL) END SELECT END DO END DO DO I=1, N WRITE (61, '(45ES14.6)') (SAOUT(I,J),J=1,N) END DO DEALLOCATE( SAOUT, STAT=K ) ENDIF INDXX = INDXX + N END DO CALL FILECLOSE( 61, 1 ) ENDIF ! --- WRITE OUT FULL SA MATRIX - OR COMMENT OUT - now in K file CALL FILEOPEN( 63, 1 ) WRITE (63, '(3X,255A14)') (PNAME(I),I = 1, NVAR) DO I=1,NVAR WRITE (63, '(255ES14.6)') (SA(I,J),J=1,NVAR) END DO CALL FILECLOSE( 63, 1 ) ! --- SAVE NVAR KEY TO DETAIL WRITE(16,502) (PNAME(I),I=1,NVAR) WRITE(16,*)'' ! WRITE(16,888)(XAPR(I),I=1,NVAR) NEGLOOP = 0 SE = SED(:NFIT) CALL FILEOPEN( 89, 1 ) WRITE(89, *) NVAR WRITE(89, '(255(I14,12X))') (I,I = 1, NVAR) WRITE(89, '(255(12X,A14))') (PNAME(I),I = 1, NVAR) IF ( ALL_SPEC_OUT ) THEN CALL FILEOPEN( 91, 1 ) WRITE(91, *) NFIT WRITE(91, *) END IF ! --- CALL OPTIMAL ESTIMATION SUBROUTINE WRITE (*, *) ' BEGIN ITERATIVE RETRIEVAL LOOP...' !CALL OPT_3(Y, PARM, XHAT, YHAT, NFIT, NVAR, CONVERGE, ITRMAX, TOL, DEBUG, RETFLG, DIVWARN, ITER, SNR, ISMIX, NLEV ) CALL OPT_3(Y, PARM, XHAT, YHAT, NFIT, NVAR, CONVERGE, ITRMAX, TOL, DEBUG, RETFLG, DIVWARN, ITER, ISMIX, NLEV ) CALL FILECLOSE( 91, 1 ) CALL FILECLOSE( 89, 1 ) ! --- IF FM HAD ABNORMAL TERMINATION, THEN JUST LEAVE IF( Tflg )GOTO 700 GOTO 202 !- DOESNT WORK WITH SNR/BAND !! --- DECREASE SNR 2 TIMES TO TRY TO GET A FIT !! --- AT START IF NEGFLAG == -2 THEN END HERE !! --- IF 0 THEN ITERATE ON NEG VMR ! J = 1 ! DO KK=1,NLEV ! IF((X(J,KK).LT.0.0).AND.(IFPRF(J))) THEN ! NEGFLAG = NEGFLAG +1 ! NEGLOOP = NEGLOOP +1 ! PRINT *,'' ! PRINT *," Final vmr profile of target contains negative values." ! PRINT *,'' ! EXIT ! ENDIF ! ENDDO ! ! IF( NEGLOOP .EQ. 5 ) GOTO 202 ! IF( NEGFLAG .LE. 0 ) THEN ! GOTO 202 ! ELSE ! SED(:NFIT) = SED(:NFIT) / (0.36) ! SNR = SNR * 0.6 ! PRINT *," Reducing SNR...", NEGFLAG, NEGLOOP, NINT(SNR) ! NEGFLAG = 0 ! GOTO 200 ! ENDIF 202 IF (.NOT.RETFLG) THEN CALL FILEOPEN( 19, 2 ) IND = 0 DO I = 1, NBAND DO K = 1, NSPEC WRITE (19, *) 'Simulation from opt_fm' WRITE (19, *) WSTART(I), WSTOP(I), SPAC WRITE (19, *) NPRIM(I) DO J = 1, NPRIM(I) WRITE (19, *) YHAT(IND+J) END DO IND = IND + NPRIM(I) END DO END DO CALL FILECLOSE( 19, 1 ) ENDIF ! --- RENORMALIZATION OF FINAL SPECTRA IND = 0 print *, maxval(yhat(:nprim(1))) print *, maxval(TOBS(:nprim(1))) DO I = 1, NBAND IF (.NOT.IEMISSION/=1.OR.IENORM(I)/=1) THEN DO K = 1, NSPEC ! WRITE(*,*)(Y(IND+J),J=1,NPRIM(I)) ! WRITE(*,*)(YHAT(IND+J),J=1,NPRIM(I)) ! YMAX=0. YMAX = 0.D0 DO J = 1, NPRIM(I) YMAX = MAX(YMAX,YHAT(IND+J)) END DO WRITE(*,*)YMAX !PRINT *, I, IND, NPRIM(I) YHAT(IND+1:NPRIM(I)+IND) = YHAT(IND+1:NPRIM(I)+IND)/YMAX TOBS(IND+1:NPRIM(I)+IND) = TOBS(IND+1:NPRIM(I)+IND)/YMAX WRITE(*,*)YMAX IND = IND + NPRIM(I) END DO ENDIF END DO ! --- FINAL UPDATE OF UNCERTAINTIES OF MIXING RATIOS WHEN APPROPRIATE INDXX = ISMIX DO KK = 1, NRET IF( IFPRF(KK) )THEN N = NLAYERS DO K = 1, N INDXX = INDXX + 1 IF (ILOGRETRIEVAL(KK)/=0) THEN !MP X(KK,K) = EXP(XHAT(INDXX)) ELSE X(KK,K) = XHAT(INDXX)*XORG(KK,K) END IF !PRINT *, 'VMR',INDXX,SHAT(INDXX,INDXX) SIG(K,KK) = SQRT(ABS(SHAT(INDXX,INDXX)))*X(KK,K) END DO ELSE INDXX = INDXX + 1 X(KK,:KMAX) = XHAT(INDXX)*XORG(KK,:KMAX) !print *, 'vmr',indxx,shat(indxx,indxx) SERR = SQRT(ABS(SHAT(INDXX,INDXX))) WRITE (16, 350) NAME(IGAS(KK)), XHAT(INDXX), SERR ENDIF END DO ! --- PRINT OUT FINAL BACKGROUND PARAMETERS IF (NBACK > 1) THEN WRITE (16, 3006) DO I = 1, NFITS SELECT CASE (NBACK) CASE (2) J = I*(NBACK - 1) !print *, 'back',j,shat(j,j) SIGB = SQRT(ABS(SHAT(J,J))) WRITE (16, 3007) I, XHAT(J), SIGB CASE (3) J = I*(NBACK - 2) !print *, 'back',j,shat(j,j) SIGB = SQRT(ABS(SHAT(J,J))) SIGC = SQRT(ABS(SHAT(J+1,J+1))) WRITE (16, 3007) I, XHAT(J), SIGB, XHAT(J+1), SIGC END SELECT END DO ENDIF ! --- PRINT OUT FINAL WAVENUMBER SHIFTS ! --- CONVERT PARAMETER TO WAVENUMBERS !DO IBAND=1, NBAND IF( ISPARM > 0 ) THEN WRITE( 16, 3008) SELECT CASE (ISPARM) CASE (1) ! SINGLE SHIFT FOR ALL BANDS (FITS = BANDS * SPECS) WRITE( 16, 3011) J = NBKFIT+1 WSHFT = 0.5D0*(WAVE3(1)+WAVE4(1))*((WAVFAC(IBAND) + XHAT(J)) - 1.D0) SWSHFT = SQRT(ABS(SHAT(J,J))) SWSHFT = 0.5D0*(WAVE3(1)+WAVE4(1))*((WAVFAC(IBAND) + SWSHFT) - 1.D0) WRITE (16, 3009) WSHFT, SWSHFT CASE (2) ! INDEPENDENT SHIFT FOR EACH BAND WRITE( 16, 3012) DO I=1, NBAND J = NBKFIT+I WSHFT = 0.5D0*(WAVE3(I)+WAVE4(I))*((WAVFAC(IBAND) + XHAT(J)) - 1.D0) SWSHFT = SQRT(ABS(SHAT(J,J))) SWSHFT = 0.5D0*(WAVE3(1)+WAVE4(1))*((WAVFAC(IBAND) + SWSHFT) - 1.D0) WRITE (16, 3014) I, WSHFT, SWSHFT ENDDO CASE (3) ! INDEPENDENT SHIFT FOR EACH FIT (BAND * SPEC) WRITE( 16, 3013) DO IBAND=1, NBAND N = NSCAN(IBAND) DO I=1, N J = NBKFIT + IBAND + I - 1 WSHFT = 0.5D0*(WAVE3(I)+WAVE4(I))*((WAVFAC(IBAND) + XHAT(J)) - 1.D0) SWSHFT = SQRT(ABS(SHAT(J,J))) SWSHFT = 0.5D0*(WAVE3(1)+WAVE4(1))*((WAVFAC(IBAND) + SWSHFT) - 1.D0) WRITE (16, 3015) IBAND, I, WSHFT, SWSHFT ENDDO ENDDO END SELECT ENDIF !ENDDO ! --- PRINT OUT ZERO LEVEL OFFSETS - RETRIEVED OR NOT IF (NZERO /= 0) THEN K = NBKFIT + NSHIFT WRITE (16, 3001) DO I = 1, NBAND !IF (IZERO(I) == 0) CYCLE N = NSCAN(I) DO J = 1, N IF( IZERO(I) == 1 ) THEN !print *, 'zero',j,shat(j,j) K = K + 1 SIGZ = SQRT(ABS(SHAT(K,K))) WRITE (16, 3002) I, J, XHAT(K), SIGZ ELSE WRITE (16, 3010) I, J, ZSHIFT(I,J), " NOT FIT" ENDIF END DO END DO ENDIF ! --- PRINT OUT SOLAR CO PARAMETERS ! --- PRINTS IF THEY ARE USED REGARDLESS OF WHETHER THEY ARE FIT IF( IFCO )THEN CPARM(4) = CPARM(4) - 1.D0 CPARM(5) = CPARM(5) - 1.D0 WRITE (16, 3003) WRITE (16, 3004) (CPNAM(I),CPARM(I),I=1,5) ENDIF ! --- PRINT OUT FINAL CHANNEL PARMS CALL PRINT_CHANNEL_PARMS( 16 ) ! --- SUM UP COLUMNS DO I = 1, NRET VERSUM(I,1) = X(I,1)*CCC(KVERT,1) VOSUM(I,1) = XORG(I,1)*CCC(KVERT,1) IF( IFPRF(I) ) FX(I) = (SIG(1,I)*CCC(KVERT,1))**2 DO K = 2, KMAX VERSUM(I,K) = VERSUM(I,K-1) + X(I,K)*CCC(KVERT,K) VOSUM(I,K) = VOSUM(I,K-1) + XORG(I,K)*CCC(KVERT,K) IF( IFPRF(I) ) FX(I) = FX(I) + (SIG(K,I)*CCC(KVERT,K))**2 END DO IF( IFPRF(I) ) FX(I) = SQRT(ABS(FX(I))) IF( I == 1 ) SERR = 100.0D0*FX(I)/VERSUM(I,KMAX) END DO AIRCOL = 0.0D0 DO K = 1, KMAX AIRCOL = AIRCOL + CCC(KVERT,K) ENDDO ! --- PRINT OUT RETRIEVED PROFILES PRINT *,'' DO K = NRET, 1, -1 IF( .NOT. IFPRF(K) )CYCLE N = NLAYERS !PRINT 408, NAME(IGAS(K)), VERSUM(K,N), FX(K), 100.0D0*FX(K)/VERSUM(K,N) !WRITE(16,408) NAME(IGAS(K)), VERSUM(K,N), FX(K), 100.0D0*FX(K)/VERSUM(K,N) WRITE(16,406) !NAME(IGAS(K)) DO KK = 1, N WRITE (16, 407) Z(KK), ZBAR(KK), XORG(K,KK), X(K,KK), SIG(KK,K), VOSUM(K,KK), VERSUM(K,KK) END DO END DO ! --- PRFS.out IF( IFPRF(1) )THEN CALL FILEOPEN( 88, 2 ) WRITE (88, 409) DO KK = 1, N WRITE (88, 407) Z(KK), ZBAR(KK), XORG(1,KK), X(1,KK), SIG(KK,1), & XORG(1,KK)*CCC(KVERT,KK), X(1,KK)*CCC(KVERT,KK), & SQRT( SA(ISMIX+KK,ISMIX+KK) ) END DO ENDIF CALL FILECLOSE( 88, 1 ) ! --- WRITING OUT OBSERVED, CALCULATED, AND DIFFERENCES - pbpfile ! --- AT EACH POINT TO TAPE8 CALL FILEOPEN( 8, 2 ) WRITE (8, *) NFITS K = 0 DO IBAND = 1, NBAND N = NSCAN(IBAND) IF (N == 0) CYCLE DO J = 1, N K = K + 1 JSCAN = ISCAN(IBAND,J) ! WRITE (8, 37905) ISPEC(JSCAN), SPAC, NPRIM(IBAND), WAVE3(IBAND), & ! WAVE4(IBAND), Z(KZTAN(JSCAN)) WRITE(8,'(A80)') STITLE(K) WRITE(8, 37905) ISPEC(JSCAN), SPAC(IBAND), NPRIM(IBAND), WSTART(IBAND), & WSTOP(IBAND), Z(KZTAN(JSCAN)) IB = 0 INDXX = ISCNDX(1,IBAND,J) -1 3790 CONTINUE IA = IB + 1 IB = IA + 11 IB = MIN0(NPRIM(IBAND),IB) !print *, indxx,ia,ib FX(:IB-IA+1) = TOBS(IA+INDXX:IB+INDXX) - YHAT(IA+INDXX:IB+INDXX) WAVE = WSTART(IBAND) + (IA - 1)*SPAC(IBAND) WRITE (8, 3794) WAVE, (TOBS(IOUT+INDXX),IOUT=IA,IB) WRITE (8, 3795) (YHAT(IOUT+INDXX),IOUT=IA,IB) WRITE (8, 3795) (FX (IOUT-IA+1), IOUT=IA,IB) IF (IB < NPRIM(IBAND)) GO TO 3790 !INDXX = INDXX + NPRIM(IBAND) END DO END DO CALL FILECLOSE( 8, 1 ) ! --- FILE OF FINAL MIXING RATIOS AND COLUMNS TO TAPE18 - statevec WRITE (18, 504) NLEV, ITER, ITRMAX, IFTEMP, CONVERGE, DIVWARN WRITE (18, *) 'z (Mid-points, Find layer boundaries in station.layers)' ! WRITE(18,*)'z' WRITE (18, 508) (ZBAR(I),I=1,NLEV) WRITE (18, *) 'p' WRITE (18, 506) (PMB(I),I=1,NLEV) WRITE (18, *) 't' WRITE (18, 506) (TORG(I),I=1,NLEV) IF( IFTEMP )THEN WRITE (18, *) 'tr' WRITE (18, 506) (T(I),I=1,NLEV) ENDIF WRITE (18, *) WRITE (18, *) NRET DO I = 1, NRET WRITE (18, 507) 'A Priori', NAME(IGAS(I)) WRITE (18, 506) VOSUM(I,NLEV) WRITE (18, 506) (XORG(I,J),J=1,NLEV) WRITE (18, 507) 'Retrieved', NAME(IGAS(I)) WRITE (18, 506) VERSUM(I,NLEV) WRITE (18, 506) (X(I,J),J=1,NLEV) END DO WRITE (18, *) WRITE (18, *) ISMIX WRITE (18, 502) (PNAME(I),I=1,ISMIX) WRITE (18, 506) (XAPR(I),I=1,ISMIX) WRITE (18, 506) (XHAT(I),I=1,ISMIX) ! --- CALCULATE DEGREES OF FREEDOM FOR SIGNAL USING APOSTERIORI SOLUTION !print *, NVAR, 'NVAR' IF( IFPRF(1) ) CALL DOFS(NFIT,NVAR,ISMIX,NLEV) ! --- SUMMARY FILE WRITE (20, *) NRET WRITE (20, 502) (NAME(IGAS(I)),I=1,NRET) WRITE (20, *) (IFPRF(I),I=1,NRET) WRITE (20, 506) (VOSUM(I,NLEV),I=1,NRET) ! prior WRITE (20, *) NBAND WRITE (20, 505) (WAVE3(I),WAVE4(I),I=1,NBAND) WRITE (20, *) WRITE (20, 505) PMAX, FOVDIA, RMS, DOF WRITE (20, 506) (VERSUM(I,NLEV),I=1,NRET) ! retrieved WRITE (20, *) ITER, ITRMAX, CONVERGE, DIVWARN WRITE (20, 3886) TFILE(18) WRITE (20, 3886) TFILE(08) WRITE (20, 3886) TFILE(16) ! --- SHORT FILE WRITE (31, 460) (NAME(IGAS(I)),IFPRF(I),I=1,NRET) ! WRITE(31,461) ( VOSUM(I,NLEV), I=1, NRET ) WRITE (31, 461) (VERSUM(I,NLEV),I=1,NRET) !WRITE (31, 462) ITER, ITRMAX, RMS, NVAR, CONVERGE, DIVWARN, DOF, SERR, SNR, AIRCOL WRITE (31, 462) ITER, ITRMAX, RMS, NVAR, CONVERGE, DIVWARN, DOF, SNR, AIRCOL DO K = NRET, 1, -1 IF( IFPRF(K) )THEN WRITE (31, 464) NAME(IGAS(K)), RMS, DOF, SNR, VERSUM(K,NLEV), AIRCOL !WRITE (31, 464) NAME(IGAS(K)), RMS, DOF, SERR, SNR, VERSUM(K,NLEV), AIRCOL ENDIF END DO ! --- PRINT A SHORT SUMMARY TO THE CONSOLE PRINT *,'' PRINT 460, ( NAME(IGAS(I)), IFPRF(I), I=1, NRET ) PRINT 461, ( VOSUM(I,NLEV), I=1, NRET ) PRINT 461, ( VERSUM(I,NLEV), I=1, NRET ) !PRINT 462, ITER, ITRMAX, RMS, NVAR, CONVERGE, DIVWARN, DOF, SERR, SNR, AIRCOL PRINT 462, ITER, ITRMAX, RMS, NVAR, CONVERGE, DIVWARN, DOF, SNR, AIRCOL PRINT *,'' ! --- MEAN ALTITUDE, VMRS, AND RETRIEVED COLUMNS TO TAPE18 ! WRITE(18,570) ZBAR(1),(X(I,1),VERSUM(I,1),I=1,NRET) ! DO 471 K=2,KZMAXLAY ! WRITE(18,570) ZBAR(K),(X(I,K),VERSUM(I,K),I=1,NRET) ! 471 CONTINUE ! --- UNCOMMENT NEXT LINE TO ACTIVATE OUTPUT OF RETRIEVED MIX FILE ! CALL MIXOUT (KZMAXLAY, KVERT) 700 CONTINUE CALL FILECLOSE( 18, 1 ) CALL FILECLOSE( 20, 1 ) CALL FILECLOSE( 31, 1 ) CALL FILECLOSE( 73, 1 ) ! --- DEALLOCATE ARRAYS DEALLOCATE( SA ) DEALLOCATE( SHAT ) DEALLOCATE( SAINV ) DEALLOCATE( KS ) DEALLOCATE( SE ) DEALLOCATE( SED ) DEALLOCATE( KSK ) DEALLOCATE( KHAT ) DEALLOCATE( G_LM ) DEALLOCATE( Z ) DEALLOCATE( ZBAR ) CALL RELEASE_MEM_LP CALL RELEASE_MEM IF( IFCO )CALL SOLARFH ( 2 ) RETURN 556 CONTINUE WRITE (16, 223) WRITE (6, *) 'NVAR =', NVAR, ' TO BIG, NMAX = ', NMAX RETURN 557 CONTINUE WRITE (16, 224) WRITE (6, *) 'NFIT=', NFIT, 'TO BIG, MAX FIT PARMS =', MMAX RETURN 664 CONTINUE WRITE (16, 591) RETURN 665 CONTINUE WRITE (16, 592) RETURN 667 CONTINUE WRITE (16, 250) RETURN 668 CONTINUE WRITE (16, 251) RETURN 223 FORMAT(/,' ABORT--TOO MANY VARIABLES') 224 FORMAT(/,' ABORT--TOO MANY SPECTRAL DATA POINTS') 250 FORMAT(/,' Can only retrieve apodization parameters for polynomial'/,& ' (irteap must be 0 unless ieap=2)') 251 FORMAT(/,' Can only retrieve phase parameters for polynomial'/,& ' (irtephs must be 0 unless iephs=2)') 350 FORMAT(/,' MOLECULE = ',A7,' PROFILE SCALE FACTOR =',F7.4,' +/-',F7.4) 406 FORMAT(' RETRIEVED VERTICAL PROFILE:',/,& ' Z[km] ZBAR[km] APRIORI_VMR RETRIEVED_VMR SIGMA_VMR APRIORI_COL RETR& &IEVE_COL') 407 FORMAT(2(F7.2),6(ES13.4)) !408 FORMAT(/' GAS: ', A7, ' COLUMN: ', ES12.4, ' +/- ',ES12.4,1X,F6.3,'%') 409 FORMAT(' Z[km] ZBAR[km] APRIORI_VMR RETRIEVED_VMR SIGMA_VMR APRIORI_COL RETR& &IEVE_COL VARIANCE') 460 FORMAT(10(2X,A7,':',L2)) 461 FORMAT(10(1P,E12.4)) 462 FORMAT(/"Itr/Mx:",I2.2,"/",I2.2," %RMS=",F5.3," FitPrm=",I3,' Cvrg:',L1, & ' Dvrg:',L1," DOFS=",F5.3, " SNR=",F5.0, " AIRCOL=",ES12.4) !' Dvrg:',L1," DOFS=",F5.3, " %CERR=",F4.1, " SNR=",F5.0, " AIRCOL=",ES12.4) !464 FORMAT('pcol ', A7, 2(F6.3,1X), 2(F6.1,1X), 7ES12.4) 464 FORMAT('pcol ', A7, 2(F6.3,1X), (F6.1,1X), 7ES12.4) 502 FORMAT(5(1X,A14)) 504 FORMAT(3I5,3L5) 505 FORMAT(5(1X,F9.3)) 506 FORMAT(5(1X,1P,E14.4)) 507 FORMAT(A20) 508 FORMAT(5(1X,F14.4)) 591 FORMAT(/,' ABORT NBACK OUT OF RANGE (1-3 VALID)') 592 FORMAT('/ ABORT ISPARM OUT OF RANGE (0-3 VALID)') 3001 FORMAT(/,' RETRIEVED ZERO LEVEL OFFSETS',/,& 'BAND ISCAN ZERO SIGMA') 3002 FORMAT(I3,I9,2ES13.5) 3003 FORMAT(/,' RETRIEVED SOLAR SIMULATION PARAMETERS:') 3004 FORMAT(2(1X,A14,'=',1P,E13.4,1X)) ! 3005 FORMAT(/,' RETRIEVED WAVENUMBER SCALE FACTOR =',F12.9) 3006 FORMAT(/' RETRIEVED BACKGROUND PARAMETERS:'/,& ' BAND SLOPE SIGMA CURVATURE SIGMA') 3007 FORMAT(I3,4(1P,E12.4)) 3008 FORMAT(/,' RETRIEVED WAVENUMBER SHIFTS:') 3009 FORMAT(2ES13.5) 3010 FORMAT(I3,I9,ES13.5,A) 3011 FORMAT(' SHIFT SIGMA') 3012 FORMAT(' BAND SHIFT SIGMA') 3013 FORMAT(' BAND SCAN SHIFT SIGMA') 3014 FORMAT(I3, 3X, 2ES13.5) 3015 FORMAT(I5, I5, 3X, 2ES13.5) 3695 FORMAT(/,' NFIT =',I5,' NVAR =',I3) 3696 FORMAT(/,' PRINT OUT OF PARAMETERS FOR EACH ITERATION:',/) 37905 FORMAT(I12,1P,E25.15,I12,0P,F25.15,/,2F25.15) 3886 FORMAT(A) 3794 FORMAT(F10.4,12ES15.5) !mp 3795 FORMAT(10X,12ES15.5) !mp RETURN END SUBROUTINE SFIT4