SUBROUTINE OPTCV(nmicro,IPRINT) #include "dimensions.h" #include "dimphy.h" #include "microtab.h" #include "clesphys.h" c Argument: c --------- integer nmicro c --------- PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1) PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25) COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL) COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL) & ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER) COMMON /VISGAS/SOLARF(NSPECV),NTERM(NSPECV),PEXPON(NSPECV), & ATERM(4,NSPECV),BTERM(4,NSPECV) COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER) & , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV) COMMON /CLOUD/ RADCLD(NLAYER), XNCLD(NLAYER) & , RCLDI(NSPECI), XICLDI(NSPECI) & , RCLDV(NSPECV), XICLDV(NSPECV) COMMON /TAUS/ TAUHI(klon,NSPECI), TAUCI(klon,NSPECI) & ,TAUGI(klon,NSPECI), TAURV(klon,NSPECV) & ,TAUHV(klon,NSPECV) ,TAUCV(klon,NSPECV) & ,TAUGV(klon,NSPECV) COMMON /TAUD/ TAUHID(klon,NLAYER,NSPECI) & ,TAUGID(klon,NLAYER,NSPECI) & ,TAUHVD(klon,NLAYER,NSPECV) & ,TAUGVD(klon,NLAYER,NSPECV) COMMON /OPTICV/ DTAUV(klon,NLAYER,NSPECV,4) & ,TAUV(klon,NLEVEL,NSPECV,4) & ,WBARV(klon,NLAYER,NSPECV,4) & ,COSBV(klon,NLAYER,NSPECV,4) COMMON /SPECTV/ BWNV(NSPC1V),WNOV(NSPECV) & ,DWNV(NSPECV),WLNV(NSPECV) COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON COMMON /CONST/ RGAS,RHOP,PI,SIGMA COMMON /traceurs/qaer(klon,nlayer,nqmx) COMMON /part/ v(nqmx),r(nqmx),vrat,dr(nqmx),dv(nqmx) REAL xv1(klev,NSPECV) REAL xv2(klev,NSPECV) REAL xv3(klev,NSPECV) REAL QF1(nqmx,NSPECV),QF2(nqmx,NSPECV) REAL QF3(nqmx,NSPECV),QF4(nqmx,NSPECV) REAL QM1(nqmx,NSPECV),QM2(nqmx,NSPECV) REAL QM3(nqmx,NSPECV),QM4(nqmx,NSPECV) save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4 integer ioptv,iwarning integer ig_ save ioptv,iwarning data ioptv,iwarning/0,0/ C* C* C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISIBLE C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VIS C LAYER: WBAR, DTAU, COSBAR C LEVEL: TAU C sum=0. PRINT*,'OPTCV' print*,'ATTENTION, TAU UNIFORME DANS OPTCV' c do nng=2,klon c do i=1,klev c do j=1,nqmx c sum=sum+qaer(nng,i,j)*r(j)**3.*1.3333*3.1415*1000. c enddo c enddo c enddo c print*,sum/(klon-1),'SOMME COLONNE/OPTCV' c open (unit=1,file='xsetupv') c do j=1,nspecv c read(1,*) a c do i=1,klev c read(1,*) xv1(i,j),xv2(i,j),xv3(i,j) c enddo c enddo c close(1) DO 130 K=1,NSPECV C LETS USE THE OPTICAL CONSTANTS FOR THOLIN CALL THOLIN(WLNV(K),TNR,TNI) REALV(K)=TNR XIMGV(K)=TNI*FHVIS C BUT WE NOW USE THE GEOMETRIC ALBEDO FITTED RESULTS C XIMGV(K)=FITEDT(WLNV(K)) C XIMGV(K)=FITEDN(WLNV(K)) C THE CLOUD IS CLEAR IN THE VISIBLE RCLDV(K)=1.27 XICLDV(K)=1.E-7 130 CONTINUE C c******* DEBUT DES BOUCLES ************************ c PRINT*, 'AEROSOLS EN VISIBLE' DO 101 ig=1,klon !c! BOUCLE SUR GRILLE HORIZONTALE DO 100 K=1,NSPECV !b! BOUCLE SUR LAMBDA C ZERO THE COLUMN OPTICAL DEPTHS OF EACH TYPE C ??FLAG? THE OPTICAL DEPTH OF THE TOP OF THE MODEL C MAY NOT BE ZERO. TAURV(ig,K)=0. TAUHV(ig,K)=0. ! INTEGRATED TAU.......INITIALIZATION. TAUCV(ig,K)=0. ! Rayleigh, Haze, Cloud, Gas TAUGV(ig,K)=0. ! sca, abs, abs , abs DO 100 J=1,NLAYER !a! BOUCLE SUR L"ALTITUDE C #1: HAZE c--------------------------- c CALL THE MIE CODE TO GIVE THE AEROSOL PROPERTIES c USE XFRAC FOR FRACTAL AEROSOLS PROPERTIES AT LAMBDA < 2. um c /\ c / \ c / \ c / _O \ c / |/ \ c / / \ \ c / |\ \/\ \ c / || / \ \ c ---------------- c | WARNING | c | SLOW DOWN | c ---------------- c*********** EN TRAVAUX *************************** TAEROS=0. TAEROSCAT=0. CBAR=0. DO inq=1,nmicro !BOUCLE SUR LES nmicro TAILLE D"AEROSOLS IF (R(inq).lt.RF(inq)) THEN ! aerosols spheriques if(ioptv.eq.0.and.J.eq.1) then c CALL XMIE(R(inq)*1.e6,REALV(K),XIMGV(K), c & QEXT,QSCT,QABS,QBAR,WNOV(K)) CALL CMIE(1.E-2/WNOV(K),REALV(K),XIMGV(K),R(inq), & QEXT,QSCT,QABS,QBAR) c print*,'inq=',inq,' QM1=',QM1(inq,K),' QEXT=',QEXT QM1(inq,K)=QEXT QM2(inq,K)=QSCT QM3(inq,K)=QABS QM4(inq,K)=QBAR endif if (microfi.eq.1) then ig_=ig else ig_=12 endif TAEROS=QM1(inq,K)*qaer(ig_,NLAYER+1-J,inq)*1.e-4+TAEROS TAEROSCAT=QM2(inq,K)*qaer(ig_,NLAYER+1-J,inq)*1.e-4+TAEROSCAT CBAR=CBAR+QM4(inq,K)*QM2(inq,K)*qaer(ig_,NLAYER+1-J,inq)*1.e-4 ELSE ! aerosols fractals XMONO=(R(inq)/RF(inq))**3. XRULE=1. if(XMONO.gt.16384./1.5) then XRULE=(XMONO/16384.) XMONO=16384. endif if(ioptv.eq.0.and.J.eq.1) then CALL OPTFRAC(XMONO,10000./WNOV(K) & ,QEXT,QSCT,QABS,QBAR) c CALL CFFFV11(1.e-2/WNOV(K),REALV(K),XIMGV(K),RF(inq),2. c & ,XMONO,QSCT,QEXT,QABS,QBAR) QF1(inq,K)=QEXT*XRULE QF2(inq,K)=QSCT*XRULE QF3(inq,K)=QABS*XRULE QF4(inq,K)=QBAR c print*,'inq=',inq,' QF1=',QF1(inq,K),' QEXT=',QEXT,' XRULE=',XRULE endif if (microfi.eq.1) then ig_=ig else ig_=12 endif TAEROS=QF1(inq,K)*qaer(ig_,NLAYER+1-J,inq)+TAEROS TAEROSCAT=QF2(inq,K)*qaer(ig_,NLAYER+1-J,inq)+TAEROSCAT CBAR=CBAR+QF4(inq,K)*QF2(inq,K)*qaer(ig_,NLAYER+1-J,inq) ENDIF ENDDO ! nmicro CBAR=CBAR/TAEROSCAT DELTAZ=Z(J)-Z(J+1) c -------------------------------------------------------------------- c profil brume Pascal: fit T (sauf tropopause) et albedo c ------------------- if( cutoff.eq.1) then IF(PRESS(J).gt.9.e-3) THEN TAEROS=TAEROSM1*DELTAZ/DELTAZM1*0.85 TAEROSCAT=TAEROSCATM1*DELTAZ/DELTAZM1*0.85 c TAEROS=0. c TAEROSCAT=0. ENDIF IF(PRESS(J).gt.1.e-1) THEN TAEROS=TAEROSM1*DELTAZ/DELTAZM1*1.15 TAEROSCAT=TAEROSCATM1*DELTAZ/DELTAZM1*1.15 c TAEROS=0. c TAEROSCAT=0. ENDIF endif !cutoff=1 c profil brume pour fit T (y compris tropopause), mais ne fit plus albedo... c ----------------------- if( cutoff.eq.2) then IF(PRESS(J).gt.1.e-1) THEN TAEROS=0. TAEROSCAT=0. ENDIF endif !cutoff=2 c -------------------------------------------------------------------- TAEROSM1=TAEROS TAEROSCATM1=TAEROSCAT DELTAZM1=DELTAZ IF (TAEROSCAT.le.0.) CBAR=0. c if(ig.eq.12) then c if(j.eq.1) print*,'NEWK',wlnv(k) c print*,j,TAEROS,xv1(j,k),' ', TAEROSCAT/TAEROS, c & xv2(j,k)/xv1(j,k),' ',CBAR,xv3(j,k) c print*,' ' c endif c print*,'HERE; MCKAY AEROSOLS VIS' c TAEROSCAT=xv2(j,k) c TAEROS=xv1(j,k) c CBAR=xv3(j,k) c if (ig.eq.1) then c if (k.eq.NSPECV/2) then c print*,'@VI',K,J,TAEROS,TAEROSCAT,CBAR c print*,'@ ',K,J,QF1(1,K),QF2(1,K),qaer(12,NLAYER+1-J,1) c print*,'@ ',K,J,QF1(3,K),QF2(3,K),qaer(12,NLAYER+1-J,3) c print*,'@ ',K,J,QF1(5,K),QF2(5,K),qaer(12,NLAYER+1-J,5) c print*,'@ ',K,J,QF1(7,K),QF2(7,K),qaer(12,NLAYER+1-J,7) c print*,'@ ',K,J,QF1(9,K),QF2(9,K),qaer(12,NLAYER+1-J,9) c print* c endif c endif c*********** EN TRAVAUX *************************** C #2: RAYLEIGH c------------------------------- C RAYLEIGH SCATTERING STRAIGHT FROM HANSEN AND TRAVIS...SEE NOTES C RATIOED BY THE LAYER COLUMN NUMBER TO THE TOTAL C COLUMN NUMBER ON EARTH. CM-2 C THIS IS THE SCATTERING BY THE ATMOSPHERE TAURAY=(COLDEN(J)*28.9/(XMU(J)*1013.25))* &(.008569/WLNV(K)**4)*(1.+.0113/WLNV(K)**2+.00013/WLNV(K)**4) c PRINT*,WLNV(K) c COLX=0. c COLP=0. c COLT=0. c DO IU=1,NLAYER c COLP=COLDEN(IU)*1.e+1*1.35+COLP c TAURAY=(COLDEN(IU)*28.9/(XMU(IU)*1013.25))* c & (.008569/WLNV(K)**4)*(1.+.0113/WLNV(K)**2 c & +.00013/WLNV(K)**4) c COLT=COLT+TAURAY c COLX=COLDEN(IU)*1.e+1/(1.E5*28./22.4E3)*1.e-1*0.0933e-1+COLX c | c | c g/cm2->kg/m2 | m2/kg c Print*,IU, tauray, c & COLDEN(IU)*1.e+1/(1.E5*28./22.4E3)*1.e-1*0.543e-1 c ENDDO c PRINT*,COLP,' PRESSURE AT GROUND;' c PRINT*,COLX,' TAU_GAS AT GROUND;' c print*,colt,colx,' COLT, COLX' c STOP c DZ=Z(J)-Z(J+1) c PRINT*, Z(J),WLNV(K), c &(28.9/(XMU(J)*1013.25))*(.008569/WLNV(K)**4)* c &(1.+.0113/WLNV(K)**2+.00013/WLNV(K)**4) c & ,COLDEN(J)/DZ/100000., c &(28.9/(XMU(J)*1013.25))*(.008569/WLNV(K)**4)* c &(1.+.0113/WLNV(K)**2+.00013/WLNV(K)**4) c & *COLDEN(J)/DZ/100000. C #3: CLOUD c---------------------------- C NEXT COMPUTE TAU CLOUD TAUCLD=0.0 c XNCLD(J)=0. IF ( XNCLD(J) .GT. 0. .and .taufac.gt.0.) THEN CALL XMIE(RADCLD(J),RCLDV(K),XICLDV(K), & QEXTC,QSCTC,QABSC,CBARC,WNOV(K)) TAUCLD=QEXTC*XNCLD(J) ENDIF C TAURV(ig,K)=TAURV(ig,K)+TAURAY TAUGVD(ig,J,K)=TAURV(ig,K) TAUHV(ig,K)=TAUHV(ig,K)+TAEROS ! INTEGRATED Quant. TAUHVD(ig,J,K)=TAUHV(ig,K) TAUCV(ig,K)=TAUCV(ig,K)+TAUCLD C #4: TAUGAS C---------------------------- C LOOP OVER THE NTERMS C THIS IS THE ABSORPTION BY THE ATMOSPHERE (METHANE) DO 909 NT=1,NTERM(K) TAUGAS=COLDEN(J)*GAS1(J)*BTERM(NT,K)* & ( (PRESS(J+1) + PRESS(J))*.5 )**PEXPON(K) C COMPUTE THE AVERAGE COSBAR AND WBAR C&& c CBAR=MIN(1.0,1.05*CBAR) ! THE HAZE FORWARD SCATTERING 5%(WHY?) COSBV(ig,J,K,NT)=(CBAR*TAEROSCAT + CBARC*TAUCLD) & /(TAEROSCAT+TAUCLD+TAURAY) !CBAR_RAY=0. c print*,'CBV',J,K,NT,CBAR,TAEROSCAT,CBARC,TAUCLD DTAUV(ig,J,K,NT)=TAUGAS+TAEROS+TAURAY+TAUCLD !TOTAL TAU_EXT TAUGV(ig,K)=TAUGV(ig,K)+TAUGAS*ATERM(NT,K) !TAU_ABS_METH INTEG. C WE LET W RAYLEIGH BE .999 OR W=1 WHEN ONLY RAYLEIGH=PROBLEM FOR TRID c WE HAVE ASSUMED ABOVE THAT COSBAR FOR RAYLEIGH IS ZERO. c if (ig.eq.1) then c if (k.eq.NSPECV/2) then c print*,'@VI',K,J,DTAUV(ig,J,K,1),TAUGAS,TAEROS,TAUCLD c endif c endif c***************** ECHANGE c WBARV(J,K,NT)=(QSCT*XNUMB(J)+TAURAY*0.9999999 + QSCTC*XNCLD(J) ) c**************** WBARV(ig,J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 + QSCTC*XNCLD(J) ) c WBARV(ig,J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 ) & /(TAUGAS+TAEROS+TAURAY+TAUCLD) c**************** IF((TAEROS+TAUCLD+TAURAY+TAUCLD).le.0.) WBARV(ig,J,K,NT)=0. IF((TAEROS+TAUCLD+TAURAY).le.0.) COSBV(ig,J,K,NT)=0. c print*,'WBV',J,K,NT,TAEROSCAT,TAURAY,QSCTC*XNCLD(J) c print*,'WBV',J,K,NT,TAEROS,TAUGAS,TAURAY,TAUCLD c print*,Z(j),J,K,NT,TAUV(1,j,K,NT),WBARV(1,j,K,NT),COSBV(1,j,K,NT) 909 CONTINUE TAUGVD(ig,J,K)=TAUGVD(ig,J,K)+TAUGV(ig,K) 100 CONTINUE ioptv=1 101 CONTINUE c HERE END OF THE LOOPS ******* c****************************** DO 102 ig=1,klon C TOTAL EXTINCTION OPTICAL DEPTHS DO 119 K=1,NSPECV C LOOP OVER NTERMS DO 119 NT=1,NTERM(K) TAUV(ig,1,K,NT)=0.0 DO 119 J=1,NLAYER TAUV(ig,J+1,K,NT)=TAUV(ig,J,K,NT)+DTAUV(ig,J,K,NT) 119 CONTINUE c print*,'SETUP' c do i=1,NSPECV c print*,WLNV(i) c do j=1,NLAYER+1 c print*,Z(j),TAUV(1,j,i,1),WBARV(1,j,i,1),COSBV(1,j,i,1) c enddo c enddo c c IF (IPRINT .GT. 1) THEN c NT=1 c IF (2 .GT. 1) THEN c WRITE (6,120) c 120 FORMAT(///' OPTICAL CONSTANTS IN THE VISIBLE (@EQUATOR) ') c WRITE(6,*) 'latitude:',ig c DO 200 K=1,NSPECV c WRITE (6,190) c WRITE (6,210)K,WLNV(K),WNOV(K),BWNV(K) c & ,BWNV(K)+DWNV(K),DWNV(K) c WRITE (6,230)REALV(K),XIMGV(K) c DO 195 J=1,NLAYER,NLAYER C RECALCULATE FOR PRINT OUT ONLY, ONLY FIRST NTERM AT ig=12 (EQUATOR) c WRITE (6,220)XNUMB(J), WBARV(ig,J,K,NT),COSBV(ig,J,K,NT) c & ,DTAUV(ig,J,K,NT),TAUV(ig,J,K,NT) c 195 CONTINUE c WRITE (6,240) TAUV(IG,NLEVEL,K,NT) c 200 CONTINUE c END IF 102 CONTINUE 210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3) 190 FORMAT(1X//' SNUM MICRONS WAVENU INTERVAL DELTA-WN') 230 FORMAT(1X,'NREAL(LAYER)= ',1PE10.3,' NIMG(LAYER)= ',E10.3/ &' #AEROSOLS WBAR COSBAR DTAU TAU' & ,9X,'RAY GAS AEROSOL') 220 FORMAT(8(1X,F9.3)) 240 FORMAT(41X,F9.3) print*,"TAUV(1400,:,10,2)=",TAUV(1400,:,10,2) print*,"DTAUV(1400,:,10,2)=",DTAUV(1400,:,10,2) c ioptv=1 PRINT*, 'FIN OPTCV' stop RETURN END