subroutine coxmunk_new & (NSTOKES, NPARS, PARS,& !I XJ, SXJ, XI, SXI,& !I CKPHI_REF, SKPHI_REF,& !I R1, Ls_R1,& !I nopr,iout,& !I mre,mim,& !I wmic,wavenumber,shadow87,comparen,exvalsign,term4,term5,& factor62,dmodmt,factor62mult,& m11,m12,m21,m22,& m11div,m12div,m21div,m22div,& theta86,r62mtm11,r62mtm12,ratiomtm21m11) ! ********************* ! Purpose of coxmunk_new is to calculate ! the reflection matrix ratios from scratch ! using Appendix A of MT (see getmtmatrix.f90 called ! in coxmunk_new.f90 ! The ratio of m12/m11 is ratiomtm21m11 ! Other values from gixxcoxmunk_vfunction_plus can be compared ! to terms calculated from breon.f90 ! This helps to see why breon.f90 and guerin.f90 sun glint values ! differs from gisscoxmunk_vfundion_plus.f90 ! ********************* ! MT refers to Mishchenko and Travis, "Satellite retrieval ! of aerosol properties over the ocean using polarization as ! well as intensity sunlight" ! JGR, vol 102, No. D14, pages 16989-17013, July 27, 1997 ! ********************* ! Subroutine arguments INTEGER NSTOKES, NPARS INTEGER nopr,iout DOUBLE PRECISION PARS(NPARS) DOUBLE PRECISION XI, SXI, XJ, SXJ DOUBLE PRECISION CKPHI_REF, SKPHI_REF DOUBLE PRECISION R1(NSTOKES) DOUBLE PRECISION Ls_R1(NSTOKES,NPARS) ! Critical exponent taken out DOUBLE PRECISION CRITEXP, PIE PARAMETER (CRITEXP = 88.D0) PARAMETER (PIE = 180.d0*1.7453292519943d-2) ! Local variables INTEGER I DOUBLE PRECISION CKPHI_INC, SKPHI_INC DOUBLE PRECISION VI1, VI2, VI3, VR1, VR2, VR3 DOUBLE PRECISION unit1, unit2, unit3, fact1, factor DOUBLE PRECISION XI1, CN1, CN2, CXI2, C2, C1, CRPER, CRPAR DOUBLE PRECISION TI1, TI2, TI3, TR1, TR2, TR3 DOUBLE PRECISION PI1, PII2, PI3, PR1, PR2, PR3 DOUBLE PRECISION PIKR, PRKI, TIKR, TRKI DOUBLE PRECISION E1, E2, E3, E4, SIGMA2 DOUBLE PRECISION CF11, CF12, CF21, CF22, RDZ2, RDZ4 DOUBLE PRECISION VP1, VP2, VP3, DMOD, DEX, DCOEFF DOUBLE PRECISION AF, AF11, AF12, AF21, AF22 DOUBLE PRECISION CTTPT, CTPPP DOUBLE PRECISION S1, S2, S3, XXI, XXJ, T1, T2, DCOT DOUBLE PRECISION SHADOWI, SHADOWR, SHADOW, DCOEFF_0, ARGUMENT DOUBLE PRECISION D_S1, D_S2, D_T1, D_T2, DERFAC DOUBLE PRECISION D_SHADOWI, D_SHADOWR, D_SHADOW DOUBLE PRECISION D_DCOEFF_0, D_DCOEFF, D_DEX, D_ARGUMENT DOUBLE PRECISION L_CXI2, L_CRPER, L_CRPAR DOUBLE PRECISION L_CF11, L_CF12, L_CF21, L_CF22 DOUBLE PRECISION L_AF11, L_AF12, L_AF21, L_AF22 ! Added to the routine double precision wmic,wavenumber double precision mre,mim double precision m11,m12,m21,m22 double precision m11div,m12div,m21div,m22div double precision ratiom21m11 double precision exvalsign,term4,term5,factor62 double precision dmodmt,factor62mult double precision r62mtm11,r62mtm12,ratiomtm21m11 double precision shadow87,theta86 ! ********************* ! Transcription of the RMATR subroutine from Mishchenko/Travis code. ! CALCULATION OF THE STOKES REFLECTION MATRIX FOR ILLUMINATION FROM ! ABOVE FOR A STATISTICALLY ROUGH SURFACE SEPARATING TWO HALF-SPACES ! WITH REFRACTIVE INDICES OF THE UPPER AND LOWER HALF-SPACES EQUAL TO ! CN1 AND CN2, RESPECTIVELY. ! SIGMA2 = s**2 = MEAN SQUARE SURFACE SLOPE (EQ. (18) IN THE JGR PAPER) ! XI = ABS(COSINE OF THE INCIDENT ZENITH ANGLE) ! XJ = ABS(COSINE OF THE REFLECTION ZENITH ANGLE). ! SXI and SXJ are the respective SINES (input) ! XPHI_REF = REFLECTION AZIMUTH ANGLE ! R1 = REFLECTION MATRIX - (1,1) and (2,1) elements only !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! Remark on Use of shadow effect ! ------------------------------ ! Shadow effect is controlled by the third parameter. That is, if ! PARS(3) not equal to 0 then shadow effect will be included. ! --- NPARS should always be 3 for this Kernel. !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! For real case, incident azimuth taken to be zero CKPHI_INC = 1.d0 SKPHI_INC = 0.d0 ! Slope square is PARS(1) SIGMA2 = PARS(1) ! ******** if (nopr .eq. 1) then write(iout,10) nstokes,npars, XJ,SXJ,XI,SXI, CKPHI_INC,SKPHI_INC, & CKPHI_REF,SKPHI_REF, SIGMA2 10 format(/,2x,'coxmunk_new: nstokes,npars ',2x,2(1x,i2),/, & 2x,'coxmunk_new: XJ,SKJ,XI,SKI are ',/,& 2x,'coxmunk_new: abs cos of reflection zenith',/,& 2x,'coxmunk_new: sine of reflection zenith',/,& 2x,'coxmunk_new: abs cos of reflection zenith',/,& 2x,'coxmunk_new: sine of reflection zenith',/,& 2x,'coxmunk_new: XJ,SKJ,XI,SKI ',4(1x,f12.5),/, & 2x,'coxmunk_new: CKPHI_INC, SKPHI_INC are cos, sine ',/,& 2x,'of incident azimuth',/,& 2x,'coxmunk_new: CKPHI_INC, SKPHI_INC ',2(1x,f12.5),/, & 2x,'coxmunk_new: CKPHI_INC, SKPHI_REF are cos, sine ',/,& 2x,'of reflected azimuth',/,& 2x,'coxmunk_new: CKPHI_REF, SKPHI_REF ',2(1x,f12.5),/, & 2x,'coxmunk_new: SIGMA2 ',1(1x,f12.5)) endif ! ******** ! Check for limiting cases IF (DABS(XI-1.D0) .LT. 1.d-9) XI = 0.999999999999d0 IF (DABS(XJ-1.D0) .LT. 1.d-9) XJ = 0.999999999999d0 ! help variables (coordinate transformations) ! vector VI is (VT1, VT2, VT3) ! note thaat VT3 is negative, component below x-y plane) VI1 = SXI * CKPHI_INC VI2 = SXI * SKPHI_INC VI3 = -XI ! vector VR is (VR1, VR2, VR3) ! note thaat VR3 is positive, component above x-y plane) VR1 = SXJ * CKPHI_REF VR2 = SXJ * SKPHI_REF VR3 = XJ ! LOCAL SURFACE NORMAL FOR SPECULAR REFLECTION (normalized to 1) ! does not seem to make sense UNIT1 = VI1-VR1 UNIT2 = VI2-VR2 UNIT3 = VI3-VR3 FACT1 = UNIT1*UNIT1 + UNIT2*UNIT2 + UNIT3*UNIT3 FACTOR = DSQRT(1.d0/FACT1) ! ******** ! Added to code for disgnostics UNIT1n = UNIT1*FACTOR UNIT2n = UNIT2*FACTOR UNIT3n = UNIT3*FACTOR ! See MT paper ! You calculate reflection matrix ratio of terms ! from scratch using the equations in MT call getmtmatrix(nopr,iout,& wmic,wavenumber,mre,mim,& XI,XJ,& SIGMA2,& VI1,VI2,VI3,& VR1,VR2,VR3,& m11,m12,m21,m22,& m11div,m12div,m21div,m22div,& theta86,r62mtm11,r62mtm12,ratiomtm21m11,& dmodmt,factor62,factor62mult,shadow87) term5=dmodmt if (nopr .eq. 1) then write(iout,fmt=15) m11,m12,m21,m22,& r62mtm11,r62mtm12,ratiomtm21m11,& factor62,dmodmt,factor62mult,shadow87,theta86 15 format(/,2x,'coxmunk_new: m11,m12,m21,m22 ',/,4(1x,f12.5),/,& 2x,'coxmunk_new: r62mtm11,r62mtm12,ratiomtm21m11 ',/,3(1x,f12.5),/,& 2x,'coxmunk_new: factor62,dmodmt,factor62mult ',/,3(1x,f12.5),/,& 2x,'coxmunk_new: shadow87 ',1(1x,f12.5),/,& 2x,'coxmunk_new: theta86 ',1(1x,f12.5)) endif ! ******** ! Back to the gisscoxmunk_vfunction_plus.f90 routine if (nopr .eq. 1) then write(iout,fmt=20) VI1,VI2,Vi3,VR1,VR2,VR3,UNIT1,UNIT2,UNIT3,FACT1,FACTOR 20 format(/,2x,'coxmunk_new: VI1,VI2,VI3 ',2x,3(1x,f12.5),/, & 2x,'coxmunk_new: VR1,VR2,VR3 ',3(1x,f12.5),/, & 2x,'coxmunk_new: UNIT1,UNIT2,UNIT3 ',3(1x,f12.5),/, & 2x,'coxmunk_new: FACT1 ',1(1x,f12.5),/, & 2x,'coxmunk_new: FACTOR ',1(1x,f12.5)) write(iout,fmt=22) UNIT1n,UNIT2n,UNIT3n 22 format(/,2x,'coxmunk_new: UNIT1n,UNIT2n,UNIT3n ',/,3(1x,f12.5)) endif ! ******** ! FRESNEL REFLECTION COEFFICIENTS, assume only real for now ! --------------------------------------------------------- CN1 = 1.d0 CN2 = PARS(2) XI1 = FACTOR*(UNIT1*VI1+UNIT2*VI2+UNIT3*VI3) CXI2 = 1.d0 - (1.d0-XI1*XI1)*CN1*CN1/(CN2*CN2) ! CN2 L_CXI2 = 2.0d0*(1.d0-CXI2)/CN2 CXI2 = DSQRT(CXI2) ! CN2 L_CXI2 = 0.5d0/CXI2*L_CXI2 C1 = CN1*XI1 C2 = CN2*CXI2 ! CN2 CRPER = (C1-C2)/(C1+C2) ! CN2 L_CRPER = -2.d0*C1*(CXI2 + CN2 * L_CXI2)/(C1+C2)**2 C1 = CN2*XI1 ! CN2 C2 = CN1*CXI2 ! CN2 CRPAR = (C1-C2)/(C1+C2) ! CN2 L_CRPAR = ((C1+C2)*(XI1-CN1*L_CXI2) - (C1-C2)*(XI1+CN1*L_CXI2))/(C1+C2)**2 !CN2 ! ******** if (nopr .eq. 1) then write(iout,30) CN1,CN2,XI1,CXI2,CRPER,CRPAR 30 format(/,2x,'coxmunk_new: CN1,CN2 ',2x,2(1x,f12.5),/, & 2x,'coxmunk_new: XI1,CXI2 ',2(1x,f12.5),/, & 2x,'coxmunk_new: CRPER ',1(1x,f12.5),/, & 2x,'coxmunk_new: CRPAR ',1(1x,f12.5)) endif ! ******** ! CALCULATION OF THE AMPLITUDE SCATTERING MATRIX ! ---------------------------------------------- TI1 = - XI * CKPHI_INC TI2 = - XI * SKPHI_INC TI3 = - SXI TR1 = + XJ * CKPHI_REF TR2 = + XJ * SKPHI_REF TR3 = -SXJ PI1 = - SKPHI_INC PII2 = + CKPHI_INC PI3 = 0.d0 PR1 = - SKPHI_REF PR2 = + CKPHI_REF PR3 = 0.d0 PIKR = PI1*VR1 + PII2*VR2 + PI3*VR3 PRKI = PR1*VI1 + PR2*VI2 + PR3*VI3 TIKR = TI1*VR1 + TI2*VR2 + TI3*VR3 TRKI = TR1*VI1 + TR2*VI2 + TR3*VI3 E1 = PIKR*PRKI E2 = TIKR*TRKI E3 = TIKR*PRKI E4 = PIKR*TRKI CF11 = E1*CRPER+E2*CRPAR ! CN2 CF12 = -E3*CRPER+E4*CRPAR ! CN2 CF21 = -E4*CRPER+E3*CRPAR ! CN2 CF22 = E2*CRPER+E1*CRPAR ! CN2 ! derivatives wrt ri L_CF11 = E1*L_CRPER + E2*L_CRPAR L_CF12 = -E3*L_CRPER + E4*L_CRPAR L_CF21 = -E4*L_CRPER + E3*L_CRPAR L_CF22 = E2*L_CRPER + E1*L_CRPAR ! ******** if (nopr .eq. 1) then write(iout,40) CF11,CF12,CF21,CF22 40 format(/,2x,'coxmunk_new: CF11,CF12,CF21,CF22 ',2x,4(1x,f12.5)) endif ! ******** ! CALCULATION OF THE STOKES REFLECTION MATRIX ! ------------------------------------------- ! vector VP is the cross product of VI and VR VP1 = VI2*VR3-VI3*VR2 VP2 = VI3*VR1-VI1*VR3 VP3 = VI1*VR2-VI2*VR1 DMOD = VP1*VP1+VP2*VP2+VP3*VP3 DMOD = DMOD*DMOD ! if DMOD = 0, that is | n x n_0 | ^ 4 = 0 in M-T formula) ! Then we need to set the ratio CF11 / DMOD IF (DABS(DMOD) .LT. 1.d-8) THEN CF11 = CRPAR ! CN2 CF22 = CRPER ! CN2 L_CF11 = L_CRPAR L_CF22 = L_CRPER DMOD = 1.d0 ENDIF ! For diagnostics ! Calculated in getmtmatrix.f90 above ! Can compare m11div with Hafermanmatrix.f90 m11 ! m11div=m11/DMOD ! m12div=m12/DMOD ! m21div=m21/DMOD ! m22div=m22/DMOD RDZ2 = UNIT3*UNIT3 RDZ4 = RDZ2*RDZ2 DCOEFF_0 = 1.d0/(8.D0*XI*XJ*DMOD*RDZ4*SIGMA2) ARGUMENT = (UNIT1*UNIT1 + UNIT2*UNIT2) / (2.d0*SIGMA2*RDZ2) IF ( ARGUMENT .GT. CRITEXP ) THEN DEX = 0.d0 ELSE DEX = DEXP(-ARGUMENT) ENDIF ! For output disgnostics exvalsign=DEX/SIGMA2 term4=(RDZ4)/(FACT1*FACT1) ! term5=DMOD DCOEFF = DCOEFF_0 * FACT1 * FACT1 * DEX ! Can compare to breon.f90 compareb comparen= DMOD * DCOEFF_0 * DEX * RDZ4 ! Derivative of DCOEFF w.r.t. SIGMA^2 D_DCOEFF = 0.D0 IF ( ARGUMENT .LT. CRITEXP ) THEN D_ARGUMENT = - ARGUMENT / SIGMA2 D_DCOEFF_0 = - DCOEFF_0 / SIGMA2 D_DEX = - D_ARGUMENT * DEX D_DCOEFF = FACT1*FACT1 * ( DCOEFF_0 * D_DEX + D_DCOEFF_0 * DEX ) ENDIF AF = 0.5d0 * DCOEFF ! ******** if (nopr .eq. 1) then write(iout,50) VP1,VP2,VP3,DMOD,RDZ2,RDZ4,DCOEFF_0,ARGUMENT,DEX,DCOEFF,AF 50 format(/,2x,'coxmunk_new: VP1,VP2,VP3 ',2x,3(1x,f12.5),/, & 2x,'coxmunk_new: DMOD ',1(1x,f12.5),/, & 2x,'coxmunk_new: RDZ2 ',1(1x,f12.5),/, & 2x,'coxmunk_new: RDZ4 ',1(1x,f12.5),/, & 2x,'coxmunk_new: DCOEFF_0 ',1(1x,f12.5),/, & 2x,'coxmunk_new: ARGUMENT ',1(1x,f12.5),/, & 2x,'coxmunk_new: DEX ',1(1x,f12.5),/, & 2x,'coxmunk_new: DCOEFF ',1(1x,f12.5),/, & 2x,'coxmunk_new: AF ',1(1x,f12.5)) write(iout,52) comparen 52 format(/,2x,'coxmunk_new: comparen and comparebreon should be similar',/,& 2x,'coxmunk_new: comparen ',2x,1(1x,f12.5)) write(iout,55) m11,m12,m21,m22,m11div,m12div,m21div,m22div 55 format(/,2x,'coxmunk_new: m11,m12,m21,m22 ',/,2x,4(1x,f12.5),/, & 2x,'coxmunk_new: m11div is m11/DMOD etc ',/,& 2x,'coxmunk_new: m11div,m12div,m21div,m22div ',/,2x,4(1x,f12.5)) endif ! ******** AF11 = DABS(CF11) ! CN2 AF12 = DABS(CF12) ! CN2 AF21 = DABS(CF21) ! CN2 AF22 = DABS(CF22) ! CN2 AF11 = AF11*AF11 ! CN2 AF12 = AF12*AF12 ! CN2 AF21 = AF21*AF21 ! CN2 AF22 = AF22*AF22 ! CN2 !derivs wrt ri L_AF11 = 2.0d0 * CF11 * L_CF11 L_AF12 = 2.0d0 * CF12 * L_CF12 L_AF21 = 2.0d0 * CF21 * L_CF21 L_AF22 = 2.0d0 * CF22 * L_CF22 term1= AF11+AF12+AF21+AF22 term2= AF11-AF22+AF12-AF21 ratioterm2term1=term2/term1 R1(1)=(AF11+AF12+AF21+AF22)*AF ! CN2 R1(2)=(AF11-AF22+AF12-AF21)*AF ! CN2 IF (NSTOKES .EQ. 3) THEN CTTPT = CF11*CF21 ! CN2 CTPPP = CF12*CF22 ! CN2 ! R1(NSTOKES) = (-CTTPT-CTPPP)*DCOEFF ! CN2 R1(NSTOKES) = (CTTPT+CTPPP)*DCOEFF ! Change sign of U since this code uses the opposite sign convention as we do ! in 2OS and LIDORT ENDIF ! ****************** if (nopr .eq. 1) then write(iout,fmt=60) AF11,AF12,AF21,AF22 60 format(/,2x,'coxmunk_new: AF11,AF12,AF21,AF22',/,& 2x,4(1x,f12.5)) write(iout,fmt=62) term1,term2,ratioterm2term1 62 format(/,2x,'coxmunk_new: term1=AF11+AF12+AF21+AF22',/,& 2x,'coxmunk_new: term2=AF11-AF22+AF12-AF21',/,& 2x,'coxmunk_new: ratioterm2term1=term2/term1',/,& 2x,'coxmunk_new: term1,term2,ratioterm2term1 ',3(1x,f12.5)) write(iout,64) R1(1),R1(2),R1(3) 64 format(/,2x,'coxmunk_new: R1(1),R1(2),R1(3) ',2x,3(1x,f12.5)) endif ! ****************** ! derivs wrt ri Ls_R1(1,2) = (L_AF11+L_AF12+L_AF21+L_AF22)*AF Ls_R1(2,2) = (L_AF11-L_AF22+L_AF12-L_AF21)*AF IF (NSTOKES .EQ. 3) & ! Ls_R1(3,2) = -(L_CF11*CF21+CF11*L_CF21+& ! L_CF12*CF22+CF12*L_CF22)*AF Ls_R1(3,2) = (L_CF11*CF21+CF11*L_CF21+& ! Change sign of U since this code uses the opposite sign convention as we L_CF12*CF22+CF12*L_CF22)*AF ! do in 2OS and LIDORT ! Derivative before shadow effect DERFAC = 0.d0 IF (DABS(DCOEFF) .GT. 1.d-8) THEN DERFAC = D_DCOEFF/DCOEFF ENDIF Ls_R1(:,1) = R1(:)*DERFAC ! Ls_R1(1,1) = (AF11+AF12+AF21+AF22)*0.5d0*D_DCOEFF ! Ls_R1(2,1) = (AF11-AF22+AF12-AF21)*0.5d0*D_DCOEFF ! IF (NSTOKES .EQ. 3) THEN ! Ls_R1(NSTOKES,1) = (-CTTPT-CTPPP)*D_DCOEFF ! ENDIF ! No Shadow code if not flagged ! Shadow parameter changed to PARS(4), 8/17/2010, V. Natraj ! You set this to 0.0d0 for printout SHADOW = 0.d0 IF (DABS(PARS(4)) .GE. 1.d-8) THEN ! Shadow code (includes derivative if flagged) S1 = DSQRT(2.d0*SIGMA2/PIE) S3 = 1.d0/(DSQRT(2.d0*SIGMA2)) S2 = S3*S3 D_S1 = 1.d0 / PIE / S1 D_S2 = - S2 / SIGMA2 SHADOWI = 0.d0 D_SHADOWI = 0.d0 IF (DABS(XI-1.d0) .GT. 1.d-8) THEN XXI = XI*XI DCOT = XI/DSQRT(1.d0-XXI) T1 = DEXP(-DCOT*DCOT*S2) T2 = DERFC(DCOT*S3) SHADOWI = 0.5d0*(S1*T1/DCOT-T2) D_T1 = - T1 * DCOT * DCOT * D_S2 D_T2 = 2.d0 * S2 * DCOT * T1 / PIE / S1 D_SHADOWI = 0.5d0 * ( D_S1*T1/DCOT + S1*D_T1/DCOT - D_T2 ) if (nopr .eq. 1) then write(iout,66) S1,T1,DCOT,T2,SHADOWI 66 format(/,2x,'coxmunk_new: S1,T1,DCOT,T2,SHADOWI ',2x,5(1x,f12.5)) endif ENDIF SHADOWR = 0.d0 D_SHADOWR = 0.d0 IF (DABS(XJ-1.d0) .GT. 1.d-8) THEN XXJ = XJ*XJ DCOT = XJ/DSQRT(1.d0-XXJ) T1 = DEXP(-DCOT*DCOT*S2) T2 = DERFC(DCOT*S3) SHADOWR = 0.5d0*(S1*T1/DCOT-T2) D_T1 = - T1 * DCOT * DCOT * D_S2 D_T2 = 2.d0 * S2 * DCOT * T1 / PIE / S1 D_SHADOWR = 0.5d0 * ( D_S1*T1/DCOT + S1*D_T1/DCOT - D_T2 ) if (nopr .eq. 1) then write(iout,68) S1,T1,DCOT,T2,SHADOWR 68 format(/,2x,'coxmunk_new: S1,T1,DCOT,T2,SHADOWR ',2x,5(1x,f12.5)) endif ENDIF SHADOW = 1.d0/(1.d0+SHADOWI+SHADOWR) if (nopr .eq. 1) then write(iout,70) SHADOWI,SHADOWR,SHADOW 70 format(/,2x,'coxmunk_new: SHADOWI,SHADOWRd,SHADOW ',2x,3(1x,f12.5)) endif R1(:) = R1(:) * SHADOW D_SHADOW = - SHADOW * SHADOW * ( D_SHADOWI + D_SHADOWR ) DO I = 1, NSTOKES Ls_R1(I,1) = Ls_R1(I,1)*SHADOW+R1(I)*D_SHADOW/SHADOW ENDDO Ls_R1(:,2) = Ls_R1(:,2) * SHADOW ! deriv wrt ri just propagates ENDIF ! 8/17/2010 Lambertian component added, V. Natraj ! PARS(3) = ALBEDO R1(1) = R1(1) + PARS(3) ! 8/17/2010 Lambertian albedo derivative added, V. Natraj ! Ls_R1(:,3) = ALBEDO DERIVATIVE Ls_R1(1,3) = 1.d0 ! Finish ! ****************** if (nopr .eq. 1) then write(iout,72) R1(1),PARS(3) 72 format(/,2x,'coxmunk_new: R1(1),PARS(3) ',2x,2(1x,f12.5)) write(iout,74) PARS(4),SHADOW 74 format(/,2x,'coxmunk_new: PARS(4),SHADOW ',2x,2(1x,f12.5)) endif ! ****************** RETURN end subroutine coxmunk_new