source: trunk/libf/phytitan/optci.F @ 110

      SUBROUTINE OPTCI(ykim,qaer,nmicro,IPRINT)
      use dimphy
#include "dimensions.h"
#include "microtab.h"
#include "numchimrad.h"
#include "clesphys.h"

c   Arguments:
c   ---------
      REAL    ykim(klon,klev,nqtot)
      real    qaer(klon,nlayer,nqtot)
      integer nmicro
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/ RADCLD(NLAYER), XNCLD(NLAYER)
     & , RCLDI(NSPECI), XICLDI(NSPECI), RCLDV(NSPECV), XICLDV(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 /TAUD/   TAUHID(ngrid,NLAYER,NSPECI)
     &               ,TAUGID(ngrid,NLAYER,NSPECI)
     &               ,TAUHVD(ngrid,NLAYER,NSPECV)
     &               ,TAUGVD(ngrid,NLAYER,NSPECV)


      COMMON /OPTICI/ DTAUI(ngrid,NLAYER,NSPECI)
     &               ,TAUI (ngrid,NLEVEL,NSPECI)
     &               ,WBARI(ngrid,NLAYER,NSPECI)
     &               ,COSBI(ngrid,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 /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
      save iopti,iwarning,seulmtunpt
      data iopti,iwarning,seulmtunpt/0,0,0/

      real   zqaer_1pt(NLAYER,nrad)
      real   TAUHID_1pt(NLAYER,NSPECI)
      real   TAUGID_1pt(NLAYER,NSPECI)
      real   TAUHI_1pt(NSPECI),TAUCI_1pt(NSPECI)
      real   TAUGI_1pt(NSPECI)
      real   DTAUI_1pt(NLAYER,NSPECI),TAUI_1pt(NLEVEL,NSPECI)
      real   WBARI_1pt(NLAYER,NSPECI)
      real   COSBI_1pt(NLAYER,NSPECI)
      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)
       if ((nmicro.ne.nrad).and.(microfi.eq.1)) then
          print*,"nmicro.ne.nrad",nmicro,nrad
          print*,"PROBLEME pour zqaer_1pt dans optci !!"
          stop
       endif

      
      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 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.eq.1) then
           do iq=1,nrad
              do j=1,NLAYER
                 zqaer_1pt(j,iq)=qaer(ig,j,iq)
              enddo
           enddo
        else
         if (ig.eq.1)  then
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
        
        if (seulmtunpt.eq.0) then
           call optci_1pt(zqaer_1pt,iopti,
     .            COSBI_1pt,DTAUI_1pt,TAUHI_1pt,TAUHID_1pt,TAUCI_1pt,
     .            TAUGI_1pt,TAUGID_1pt,WBARI_1pt,TAUI_1pt,IPRINT)
           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,:,:)  = COSBI_1pt(:,:) 
        WBARI(ig,:,:)  = WBARI_1pt(:,:) 
        DTAUI(ig,:,:)  = DTAUI_1pt(:,:) 
        TAUHI(ig,:)    = TAUHI_1pt(:) 
        TAUCI(ig,:)    = TAUCI_1pt(:) 
        TAUGI(ig,:)    = TAUGI_1pt(:) 
        TAUI(ig,:,:)   = TAUI_1pt(:,:) 
        TAUHID(ig,:,:) = TAUHID_1pt(:,:) 
        TAUGID(ig,:,:) = TAUGID_1pt(:,:) 

c************************************************************************
c************************************************************************
  79  CONTINUE   ! FIN BOUCLE GRILLE HORIZONTALE
c************************************************************************
c************************************************************************

      print*, 'FIN OPTCI'

      RETURN
      END