SUBROUTINE OPTCI(ykim,nmicro,IPRINT) #include "dimensions.h" #include "dimphy.h" #include "microtab.h" #include "numchimrad.h" #include "clesphys.h" c Arguments: c --------- REAL ykim(klon,klev,nqmx) 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 /STRATO/ C2H2(NLAYER),C2H6(NLAYER) COMMON /STRAT2/ HCN(NLAYER) 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 /OPTICI/ DTAUI(klon,NLAYER,NSPECI) & ,TAUI (klon,NLEVEL,NSPECI) & ,WBARI(klon,NLAYER,NSPECI) & ,COSBI(klon,NLAYER,NSPECI) COMMON /SPECTI/ BWNI(NSPC1I), WNOI(NSPECI), & DWNI(NSPECI), WLNI(NSPECI) 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 COMMON /part/v,r,vrat,dr,dv DIMENSION PROD(NLEVEL) real qaer(klon,nlayer,nqmx) real v(nqmx),r(nqmx),vrat,dr(nqmx),dv(nqmx) real xv1(klev,nspeci),xv2(klev,nspeci) real xv3(klev,nspeci) REAL QF1(nqmx,NSPECI),QF2(nqmx,NSPECI) REAL QF3(nqmx,NSPECI),QF4(nqmx,NSPECI) REAL QM1(nqmx,NSPECI),QM2(nqmx,NSPECI) REAL QM3(nqmx,NSPECI),QM4(nqmx,NSPECI) real emu REAL TAEROSM1(NSPECI),TAEROSCATM1(NSPECI),DELTAZM1(NSPECI) save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4 integer iopti,iwarning integer ig_ save iopti,iwarning data iopti,iwarning/0,0/ C THE PRESSURE INDUCED TRANSITIONS ARE FROM REGIS C THE LAST SEVENTEEN INTERVALS ARE THE BANDS FROM GNF. C C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE INFRARED C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE IR C LAYER: WBAR, DTAU, COSBAR C LEVEL: TAU C print*,'START OPTCI' c sum=0. c do nng=2,klon c do i=1,klev c do j=1,nqmx c print*,'j,rj',j,r(j) c print*,'paer',qaer(nng,i,j) 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/OPTCI' c do inq=1,nqmx c print*,inq,r(inq),vrat,qaer(12,25,inq) c enddo DO 420 K=1,NSPECI C LETS USE THE THOLIN OPTICAL CONSTANTS FOR THE HAZE. CALL THOLIN(WLNI(K),TNR,TNI) REALI(K)=TNR XIMGI(K)=TNI*FHIR C SET UP THE OPTICAL CONSTANTS FOR THE CLOUD RCLDI(K)=1.27 XICLDI(K)=REFLIQ(WNOI(K)) 420 CONTINUE C C C ZERO ALL OPTICAL DEPTHS. C ??FLAG? FOR SOME APPLCIATIONS THE TOP OPACITY MAY NOT VANISH c open (unit=1,file='xsetupi') c do i=1,klev c read(1,*) a c do j=1,nspeci c read(1,*) xv1(i,j),xv2(i,j),xv3(i,j) c enddo c enddo c close(1) c************************************************************************ c************************************************************************ DO 79 ig=1,klon ! BOUCLE SUR GRILLE HORIZONTALE c print*,'ig NEW optci',ig c************************************************************************ c************************************************************************ DO 80 K=1,NSPECI TAUHI(ig,K)=0. TAUCI(ig,K)=0. TAUGI(ig,K)=0. 80 CONTINUE if (.not.ylellouch) then XN2(1) = ykim(ig,1,iradn2) CH4(1) = ykim(ig,1,iradch4) H2(1) = ykim(ig,1,iradh2) do j=2,nlayer XN2(j) = (ykim(ig,j,iradn2)+ykim(ig,j-1,iradn2))/2. CH4(j) = (ykim(ig,j,iradch4)+ykim(ig,j-1,iradch4))/2. H2(j) = (ykim(ig,j,iradh2)+ykim(ig,j-1,iradh2))/2. enddo XN2(nlevel) = ykim(ig,nlayer,iradn2) CH4(nlevel) = ykim(ig,nlayer,iradch4) H2(nlevel) = ykim(ig,nlayer,iradh2) do j=1,nlayer emu = ( xmu(j) + xmu(j+1) )/2. C2H2(j) = ykim(ig,j,iradc2h2) * 26./emu C2H6(j) = ykim(ig,j,iradc2h6) * 30./emu HCN(j) = ykim(ig,j,iradhcn ) * 27./emu enddo endif c************************************************************************ DO 100 J=1,NLAYER ! BOUCLE SUR L'ALTITUDE c************************************************************************ c print*,'ig,k,j ',ig,k,j C SET UP THE COEFFICIENT TO REDUCE MASS PATH TO STP ...SEE NOTES C T0 =273.15 PO=1.01325 BAR TBAR=0.5*(TEMP(J)+TEMP(J+1)) PBAR=SQRT(PRESS(J)*PRESS(J+1)) BMU=0.5*(XMU(J+1)+XMU(J)) COEF1=RGAS*273.15**2*.5E5* (PRESS(J+1)**2 - PRESS(J)**2) & /(1.01325**2 *EFFG(Z(J))*TBAR*BMU) IF (IPRINT .GT. 9) WRITE(6,21) J,EFFG(Z(J)),TBAR,BMU,COEF1 21 FORMAT(' J, EFFG, TBAR, BMU, COEF1,: ',I3,1P6E10.3) c------------------------------------------------------------------------ DO 101 K=1,NSPECI ! BOUCLE SUR LES L.D'O c------------------------------------------------------------------------ C #1: HAZE C--------------------- C FIRST COMPUTE TAU AEROSOL c 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 (WNOI(K).lt.wco) THEN ! lamda > 56 um if (iopti.eq.0) then c CALL XMIE(R(inq)*1.e6,REALI(K),XIMGI(K), c & QEXT,QSCT,QABS,QBAR,WNOI(K)) CALL CMIE(1.E-2/WNOI(K),REALI(K),XIMGI(K),R(inq), & QEXT,QSCT,QABS,QBAR) QM1(inq,K)=QEXT QM2(inq,K)=QSCT QM3(inq,K)=QABS QM4(inq,K)=QBAR endif ! end iopti 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 ! 0.2 < lambda < 56 um if(R(inq).lt.RF(inq)) THEN if (iopti.eq.0) then CALL XMIE(R(inq)*1.e6,REALI(K),XIMGI(K), & QEXT,QSCT,QABS,QBAR,WNOI(K)) 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 XMONO=(r(inq)/RF(inq))**3. XRULE=1. if(XMONO.gt.16384./1.5) then XRULE=(XMONO/16384.) XMONO=16384. endif if (iopti.eq.0) then CALL OPTFRAC(XMONO,10000./WNOI(K) & ,QEXT,QSCT,QABS,QBAR) c CALL CFFFV11(1.e-2/WNOI(K),REALI(K),XIMGI(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 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 IF(TAEROS.LT.1.e-10) TAEROS=1.e-10 ENDIF ENDDO ! FIN DE LA BOUCLE SUR 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(K)*DELTAZ/DELTAZM1(K)*0.85 TAEROSCAT=TAEROSCATM1(K)*DELTAZ/DELTAZM1(K)*0.85 c TAEROS=0. c TAEROSCAT=0. ENDIF IF(PRESS(J).gt.1.e-1) THEN TAEROS=TAEROSM1(K)*DELTAZ/DELTAZM1(K)*1.15 TAEROSCAT=TAEROSCATM1(K)*DELTAZ/DELTAZM1(K)*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(K)=TAEROS TAEROSCATM1(K)=TAEROSCAT DELTAZM1(K)=DELTAZ IF(TAEROSCAT.le.0.) CBAR=0. c print*,'HERE, MCKAY AEROSOLS IR' c TAEROS=xv1(j,k) c TAEROSCAT=xv2(j,k) c CBAR=xv3(j,k) c if (ig.eq.1) then c if (k.eq.NSPECV/2) then c print*,'@IR',K,J,TAEROS,TAEROSCAT,CBAR c stop'Pour faire des comparaisons' c endif c endif c*********** EN TRAVAUX *************************** C #2: CLOUD C------------------ C NEXT COMPUTE TAU CLOUD TAUCLD=0.0 IF ( XNCLD(J) .GT. 0..and .taufac.gt.0.) THEN print*,'Appel a xmie avec radcld=',radcld(j) CALL XMIE(RADCLD(J),RCLDI(K),XICLDI(K), & QEXTC,QSCTC,QABSC,CBARC,WNOI(K)) TAUCLD=QEXTC*XNCLD(J) ENDIF C #3: GAZ C------------------ C NOW COMPUTE TAUGAS DUE TO THE PIA TERM ONLY FOR LAMDA LT 940 TAUGAS=0.0 IF (WNOI(K) .LT. 940. ) THEN c if(ig.eq.1.and.k.eq.nspecv/2) print*,'avant PIA' CALL PIA(K,TBAR,PNN,PCC,PCN,PHN) c if(ig.eq.1.and.k.eq.nspecv/2) print*,'apres PIA' C HERE IS WHERE WE COULD SCALE THE PIA COEFFICEINTS TO FIT DATA C BASED ON REGIS' NOTES. ---TGM HAS THIS ADJUST IN IT AS DEFAULT PCN=PCN*MIN(1.75 , AMAX1(1.0,WNOI(K)/200.)) C***REPLACE ABOVE WITH: PCN=PCN*1.25*MIN(1.75 , AMAX1(1.0,WNOI(K)/200.)) C 1.25 FACTOR (NOT FROM DATA) SUGGESTED BY TOON et al. (1988) TAUGAS=COEF1* & (XN2(J)*XN2(J)*PNN + CH4(J)*CH4(J)*PCC & + XN2(J)*CH4(J)*PCN + XN2(J)*H2(J)*PHN) IF (J .EQ. NLAYER .AND. IPRINT .GT. 9) & WRITE (6,22) WNOI(K),TAUGAS,XN2(J),CH4(J),H2(J), & TBAR, PNN,PCC,PCN, PHN, & XN2(J)*XN2(J)*PNN , CH4(J)*CH4(J)*PCC , & XN2(J)*CH4(J)*PCN , XN2(J)*H2(J)*PHN 22 FORMAT(1X,1P8E10.2) ENDIF IF (K .GT. 28) THEN KGAS=K-28 C ??FLAG? HERE MUST BE WATCHED CAREFULLY U=COLDEN(J)*6.02204E23/BMU if(ig.eq.1.and.k.eq.nspecv/2) print*,'Avant GAS2' if((ylellouch).or.(.not.hcnrad)) then CALL GAS2_NOHCN(J, KGAS,TBAR,PBAR,U,TAU2) else CALL GAS2(J, KGAS,TBAR,PBAR,U,TAU2) endif if(ig.eq.1.and.k.eq.nspecv/2) print*,'Apres GAS2' TAUGAS=TAUGAS+TAU2 ENDIF C IF (TAEROS + TAUCLD .GT. 0.) THEN COSBI(ig,J,K)=(CBAR*TAEROSCAT + CBARC*TAUCLD ) & /(TAEROSCAT + TAUCLD) ELSE COSBI(ig,J,K)=0.0 ENDIF DTAUI(ig,J,K)=TAUGAS+TAEROS+TAUCLD TAUHI(ig,K)=TAUHI(ig,K) + TAEROS TAUHID(ig,J,K)=TAUHI(ig,K) TAUGI(ig,K)=TAUGI(ig,K) + TAUGAS TAUGID(ig,J,K)=TAUGI(ig,K) TAUCI(ig,K)=TAUCI(ig,K) + TAUCLD c if (ig.eq.1) then c if (k.eq.NSPECI/2.or.k.eq.1.or.k.eq.NSPECI) then c print*,'@IR',K,J,DTAUI(ig,J,K),TAUGAS,TAEROS,TAUCLD c endif c endif C ??FLAG? SERIOUS PROBLEM WITH THE CODE HERE! TLIMIT=1.E-16 IF (DTAUI(ig,J,K) .GT. TLIMIT) THEN c***************** ECHANGE c WBARI(J,K)=(QSCT*XNUMB(J) + QSCTC*XNCLD(J)) c**************** CFC WBARI(ig,J,K)=(TAEROSCAT + QSCTC*XNCLD(J)) c**************** WBARI(ig,J,K)=TAEROSCAT & /DTAUI(ig,J,K) if(iwarning.eq.0) s print*,'WARNING!!! dans optci xncld non initialise' iwarning=1 ELSE WBARI(ig,J,K)=0.0 c PRINT*,'gaz ',taugas,'aerosols',taeros,'nuages',taucld c WRITE (6,999) J,K,DTAUI(ig,J,K) 999 FORMAT (' WARNING TAU MINIMUM AT J,K,DTAU:',2I3,1PE10.3) DTAUI(ig,J,K)=TLIMIT ENDIF c IF (IPRINT .GT. 9) c & WRITE(6,73)J,K,TAUGAS,TAEROS,QEXT,QSCT 73 FORMAT(2I3,1P8E10.3) c------------------------------------------------------------------------ 101 CONTINUE ! FIN BOUCLE L D'O c------------------------------------------------------------------------ iopti=1 c************************************************************************ 100 CONTINUE ! FIN BOUCLE ALTITUDE c************************************************************************ c************************************************************************ c************************************************************************ 79 CONTINUE ! FIN BOUCLE GRILLE HORIZONTALE c************************************************************************ c************************************************************************ C TOTAL EXTINCTION OPTICAL DEPTHS DO 78 ig=1,klon ! #1 DO 119 K=1,NSPECI TAUI(ig,1,K)=0.0 DO 119 J=1,NLAYER TAUI(ig,J+1,K)=TAUI(ig,J,K)+DTAUI(ig,J,K) 119 CONTINUE IF (IPRINT .GT. 2) THEN c WRITE (6,120) 120 FORMAT(///' OPTICAL CONSTANTS IN THE INFRARED') c DO 200 K=1,NSPECI ! #2 c WRITE (6,190) c WRITE (6,210)K,WLNI(K),WNOI(K),BWNI(K) c & ,BWNI(K)+DWNI(K),DWNI(K) c WRITE (6,230)REALI(K),XIMGI(K) c DO 195 J=1,NLAYER ! #3 c WRITE (6,220)XNUMB(J), WBARI(ig,J,K),COSBI(ig,J,K) c & , DTAUI(ig,J,K),TAUI(ig,J,K) c 195 CONTINUE c IF(ig.eq.12) WRITE (6,240) TAUI(ig,NLEVEL,K) c 200 CONTINUE END IF 78 CONTINUE print*,"TAUI(1400,:,10)=",TAUI(1400,:,10) print*,"DTAUI(1400,:,10)=",DTAUI(1400,:,10) 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') 220 FORMAT(5(1X,G9.3)) 240 FORMAT(41X,G9.3) c iopti=1 print*, 'FIN OPTCI' stop RETURN END