SUBROUTINE OPTCI(ykim,qaer,nmicro,IPRINT) use dimphy use infotrac_phy, only: nqtot use common_mod, only:rmcbar,xfbar,ncount,TauHID,TauCID,TauGID USE TGMDAT_MOD, ONLY: RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC, & RCLOUD,FARGON #include "dimensions.h" #include "microtab.h" #include "numchimrad.h" #include "clesphys.h" c Arguments: c --------- REAL ykim(klon,klev,nqtot) real qaer(klon,klev,nqtot) integer nmicro,IPRINT c --------- c ASTUCE POUR EVITER klon... EN ATTENDANT MIEUX INTEGER ngrid PARAMETER (ngrid=(jjm-1)*iim+2) ! = klon 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/ & RCLDI(NSPECI), XICLDI(NSPECI) & , RCLDV(NSPECV), XICLDV(NSPECV) & , RCLDI2(NSPECI), XICLDI2(NSPECI) & , RCLDV2(NSPECV), XICLDV2(NSPECV) COMMON /TAUS/ TAUHI(ngrid,NSPECI),TAUCI(ngrid,NSPECI), & TAUGI(ngrid,NSPECI),TAURV(ngrid,NSPECV), & TAUHV(ngrid,NSPECV),TAUCV(ngrid,NSPECV), & TAUGV(ngrid,NSPECV) COMMON /OPTICI/ DTAUI(ngrid,NLAYER,NSPECI) & ,TAUI (ngrid,NLEVEL,NSPECI) & ,WBARI(ngrid,NLAYER,NSPECI) & ,COSBI(ngrid,NLAYER,NSPECI) & ,DTAUIP(ngrid,NLAYER,NSPECI) & ,TAUIP(ngrid,NLEVEL,NSPECI) & ,WBARIP(ngrid,NLAYER,NSPECI) & ,COSBIP(ngrid,NLAYER,NSPECI) COMMON /SPECTI/ BWNI(NSPC1I), WNOI(NSPECI), & DWNI(NSPECI), WLNI(NSPECI) REAL DTAUP(ngrid,NLAYER,NSPECI),DTAUPP(ngrid,NLAYER,NSPECI) COMMON /IRTAUS/ DTAUP,DTAUPP COMMON /part/v,rayon,vrat,dr,dv DIMENSION PROD(NLEVEL) * nrad dans microtab.h real v(nrad),rayon(nrad),vrat,dr(nrad),dv(nrad) real xv1(klev,nspeci),xv2(klev,nspeci) real xv3(klev,nspeci) REAL QF1(nrad,NSPECI),QF2(nrad,NSPECI) REAL QF3(nrad,NSPECI),QF4(nrad,NSPECI) REAL QM1(nrad,NSPECI),QM2(nrad,NSPECI) REAL QM3(nrad,NSPECI),QM4(nrad,NSPECI) real emu REAL TAEROSM1(NSPECI),TAEROSCATM1(NSPECI),DELTAZM1(NSPECI) save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4 integer iopti,iwarning ! iopti: premier appel, une seule boucle sur les l.d'o. integer ig,seulmtunpt,iout save iopti,iwarning,seulmtunpt data iopti,iwarning,seulmtunpt/0,0,0/ real zqaer_1pt(NLAYER,2*nrad) #include "optci_1pt.h" character*100 dummy real dummy2,dummy3 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 Diagnostic eventuellement: c if (nmicro.gt.0) then c sum=0. c do nng=2,klon c do i=1,klev c do j=1,nmicro c print*,'j,rj',j,rayon(j) c print*,'paer',qaer(nng,i,j) c sum=sum+qaer(nng,i,j)*rayon(j)**3.*1.3333*3.1415*1000. c enddo c enddo c enddo c print*,sum/(klon-1),'SOMME COLONNE/OPTCI' c endif c do inq=1,nrad c print*,inq,rayon(inq),vrat,qaer(12,25,inq) c enddo C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c INITIALISATIONS UNE SEULE FOIS C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ if (iopti.eq.0) then c verif pour taille zqaer_1pt, sachant que si microfi=0 et nqtot=1, c il faut quand meme qu on lise la look-up table de dim nrad=10 c et si microfi=1, on doit avoir nmicro=nrad (dans microtab.h) c c Nouvelle verif pour nuages : c La condition ci-dessus n'est plus realisable ! c nmicro comprend maintenant aussi des glaces c Donc on teste juste que nmicro soit > 2*nrad (ou nrad si on ne fait pas de nuages) if (microfi.ge.1) then if ((clouds.eq.1).and.(nmicro.lt.2*nrad)) then print*,"OPTCI :" print*,"clouds = 1 MAIS nmicro < 2*nrad" print*,"Probleme pour zqaer_1pt dans optci." stop endif if ((clouds.eq.0).and.(nmicro.lt.nrad)) then print*,"OPTCI :" print*,"nmicro < nrad" print*,"Probleme pour zqaer_1pt dans optci." stop endif endif DO 420 K=1,NSPECI C LETS USE THE THOLIN OPTICAL CONSTANTS FOR THE HAZE. c CALL THOLIN(WLNI(K),TNR,TNI) CALL THOLIN_CVD(WLNI(K),TNR,TNI) REALI(K)=TNR XIMGI(K)=TNI*FHIR C SET UP THE OPTICAL CONSTANTS FOR THE CLOUD CALL LIQCH4(WLNI(K),TNR,TNI) RCLDI(K)=TNR XICLDI(K)=TNI CALL LIQC2H6(WLNI(K),TNR,TNI) RCLDI2(K)=TNR XICLDI2(K)=TNI 420 CONTINUE c DEBUG c print*,"wnoi=",WNOI 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) endif ! fin initialisations premier appel c************************************************************************ c************************************************************************ DO 79 ig=1,klon ! BOUCLE SUR GRILLE HORIZONTALE c print*,'ig NEW optci',ig c************************************************************************ c************************************************************************ 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 if ((.not.ylellouch).and.(ig.eq.klon/2)) then c print*,' LAYER C2H2 C2H6 HCN masmix ratios' c do j=1,nlayer c print*,j,C2H2(j),C2H6(j),HCN(j) c enddo c endif if (microfi.ge.1) then do iq=1,2*nrad c si on ne fait pas de nuages on ne copie que les nrad premieres valeurs. if (clouds.eq.0.and.iq.gt.nrad) then zqaer_1pt(:,iq)=0. else do j=1,NLAYER zqaer_1pt(j,iq)=qaer(ig,j,iq) enddo endif enddo else if (ig.eq.1) then zqaer_1pt = 0. c initialisation zqaer_1pt a partir d une look-up table (uniforme en ig) c boucle sur nrad=10 (dans microtab.h) open(10,file="qaer_eq_1d.dat") do iq=1,15 read(10,'(A100)') dummy enddo do j=NLAYER,1,-1 read(10,*) dummy2,dummy3,(zqaer_1pt(j,iq),iq=1,nrad) enddo close(10) c ici, les tableaux definissant la structure des aerosols sont c remplis: rf,df(nq),rayon(nq,)v(nq)...... call rdf() endif endif c if ((ig.eq.klon/2).or.(microfi.eq.0)) then c print*,"Q01=",zqaer_1pt(:,1) c print*,"Q05=",zqaer_1pt(:,5) c print*,"Q10=",zqaer_1pt(:,10) c stop c endif iout=0 c if ((microfi.eq.0).or.(ig.eq.(klon/2+16))) iout=1 if (seulmtunpt.eq.0) then call optci_1pt3(zqaer_1pt,rmcbar(ig,:),xfbar(ig,:,:), & iopti,iout) iopti = 1 endif c Pas de microphysique, ni de composition variable: un seul passage c dans optci_1pt. if ((microfi.eq.0).and.(ylellouch)) then seulmtunpt = 1 endif COSBI(ig,:,:) = MAX(MIN(COSBI_1pt(:,:),0.999999),1e-6) WBARI(ig,:,:) = MAX(MIN(WBARI_1pt(:,:),0.999999),1e-6) DTAUI(ig,:,:) = DTAUI_1pt(:,:) TAUI(ig,:,:) = TAUI_1pt(:,:) COSBIP(ig,:,:) = MAX(MIN(COSBIP_1pt(:,:),0.999999),1e-6) WBARIP(ig,:,:) = MAX(MIN(WBARIP_1pt(:,:),0.999999),1e-6) DTAUIP(ig,:,:) = DTAUIP_1pt(:,:) TAUIP(ig,:,:) = TAUIP_1pt(:,:) TAUHI(ig,:) = TAUHI_1pt(:) TAUCI(ig,:) = TAUCI_1pt(:) TAUGI(ig,:) = TAUGI_1pt(:) TauHID(ig,:,:) = TAUHID_1pt(:,:) TauCID(ig,:,:) = TAUCID_1pt(:,:) TauGID(ig,:,:) = TAUGID_1pt(:,:) c DEBUG c if(ig.eq.(ngrid/2+16)) then c print*,ig,'/',KLON,':' c print*,'TauHID_1',TAUHID(ig,1,:) c print*,'TauGID_1',TAUGID(ig,1,:) c print*,'TauHID_50',TAUHID(ig,50,:) c print*,'TauGID_50',TAUGID(ig,50,:) c print*,"DTAUI_1=",DTAUI(ig,1,:) c print*,"DTAUI_50=",DTAUI(ig,50,:) c print*,'cosBI_1',COSBI(ig,1,:) c print*,'cosBI_50',COSBI(ig,50,:) c print*,'WBARI_1',WBARI(ig,1,:) c print*,'WBARI_50',WBARI(ig,50,:) c stop c endif c************************************************************************ c************************************************************************ 79 CONTINUE ! FIN BOUCLE GRILLE HORIZONTALE c************************************************************************ c************************************************************************ C THIS ROUTINE HAS ALREADY SET THE DTAUI(J,K) VALUES BUT MUST BE PASSED DO 225 IG=1,klon DO 220 J=1,NLAYER DO 230 K=1,NSPECI DTAUP(IG,J,K)=DTAUI(IG,J,K) DTAUPP(IG,J,K)=DTAUIP(IG,J,K) 230 CONTINUE 220 CONTINUE 225 CONTINUE print*, 'FIN OPTCI' RETURN END