source: trunk/LMDZ.VENUS/libf/phyvenus/compo_hedin_mod2.F90 @ 3900

  module compo_hedin83_mod2
  !!05/2013 Laura Salmi
  !!03/2014 revision Gabriella Gilli
  !!Calcul des vmr pour CO2, CO, O, N et N2 en s'appuyant sur les tables donnees dans l'article Hedin (1983)
  
  implicit none      
  contains

        subroutine compo_hedin83_init2
        use YOMCST_mod
        implicit none
!#include "YOMCST.h"
#include "hedin.h"

        REAL :: alpha
        REAL :: mu(musize),z(zsize)
        REAL :: T(musize,zsize)
        REAL :: ntot(musize, zsize)
        REAL :: nco2(musize, zsize),nn2(musize, zsize)
        REAL :: nco(musize, zsize), no(musize, zsize)
        REAL :: nn(musize, zsize),nhe(musize, zsize)
        REAL :: sig_co2(zsize), sig_o(zsize), sig_n2(zsize),sig_co(zsize)
        integer :: i,j
        REAL factor_ox3p


! Initialisation des sza et de l'atlitude
        mu(1)=1.
        do i=2,musize-1
              mu(i)=mu(i-1)-(2./(1.*musize-1.))
        enddo
        mu(musize)=-1.

        z(1)=100.
        do i=2,zsize
              z(i)=z(i-1)+5.
        enddo

!Lecture des variables (tables de Hedin (1983))

       open( 11, file='HIGHATM/noon.txt')
      read (11, *)
      do i=1,zsize
         read (11,*) T(1,i), nco2(1,i),no(1,i), nco(1,i), nn2(1,i),nn(1,i), nhe(1,i)
      enddo
      close (11)
       open( 12, file='HIGHATM/midnight.txt')
      read (12, *)
      do i=1,zsize
         read (12,*) T(musize,i), nco2(musize,i),no(musize,i), nco(musize,i), nn2(musize,i), nn(musize,i), nhe(musize,i)
      enddo
      close (12)

!Dependance en SZA

        do i=2,musize-1
          alpha=0.5*tanh(4.*(mu(i)-0.5))+0.5
          do j=1,zsize 
                T(i,j)=alpha*T(1,j)+(1.-alpha)*T(musize,j)
                nco2(i,j)=alpha*nco2(1,j)+(1.-alpha)*nco2(musize,j)
                no(i,j)=alpha*no(1,j)+(1.-alpha)*no(musize,j)
                nco(i,j)=alpha*nco(1,j)+(1.-alpha)*nco(musize,j)
                nn2(i,j)=alpha*nn2(1,j)+(1.-alpha)*nn2(musize,j)
                nn(i,j)=alpha*nn(1,j)+(1.-alpha)*nn(musize,j)
                nhe(i,j)=alpha*nhe(1,j)+(1.-alpha)*nhe(musize,j)

          enddo
         enddo

!! Test: effect of varying atomic oxygen abundances

         factor_ox3p = 2.
           DO i=1,musize
             DO j=1,zsize
               no(i,j)    =  factor_ox3p* no(i,j)
             ENDDO
           ENDDO






!! Conversion en volume mixture ratio

           do i=1,musize
                  do j=1,zsize 
                  ntot(i,j)=nco2(i,j)+nn2(i,j)+nn(i,j)+nco(i,j)+no(i,j)+nhe(i,j)+2e-3*nco2(i,j) !for o2 and no si on change la proportion changer aussi dans euvheat
                  pres_hedin(i,j)=ntot(i,j)*1e6*RKBOL*T(i,j)
                 co2_hedin(i,j)=nco2(i,j)/ntot(i,j)
                 co_hedin(i,j)=nco(i,j)/ntot(i,j)
                 n2_hedin(i,j)=nn2(i,j)/ntot(i,j)
                 n_hedin(i,j)=nn(i,j)/ntot(i,j)
                 o_hedin(i,j)=no(i,j)/ntot(i,j)
                 mu_hedin(i)=mu(i)

          enddo
         enddo

!           print*, pres_hedin(1,:)
!           print*, ' '
!           print*, T(musize,:)  
!           print*, ' '
!            print*, T(1,:)
!               print*, nn2(10,:)
!               print*, ' '
!               print*,nco2(10,:)
!               print*, ' '
!               print*,nco(10,:)
!               print*, ' '
!               print*,no(10,:)

!               print*, " "
!               print*, ntot(1,:)
!               print*, co2_hedin(1,:)
!               print*, mu_hedin(:)
!               print*, " "
!               print*, n2_hedin(1,:)
!               print*, " "
!               print*, o_hedin(1,:)
!               print*, " "
!               print*, co_hedin(1,:)
!               print*, ' '
!           stop

       end subroutine compo_hedin83_init2



SUBROUTINE compo_hedin83_mod(pression,mu0,co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)


!!Interpolation des profils de Hedin (1983) sur la grille du modele

        use dimphy
        implicit none
        include 'hedin.h'
        REAL, intent(in) :: pression(klon,klev), mu0(klon)
        integer :: i,k,iz,jmu,jz,z_ok, mu_ok
        REAL :: factp(klev)
        REAL, intent(out) :: co2vmr_gcm(klon,klev),covmr_gcm(klon,klev)
        REAL, intent(out) ::ovmr_gcm(klon,klev), n2vmr_gcm(klon,klev), nvmr_gcm(klon,klev)
        REAL :: ang0, ang_hedin
!        print*, pres_hedin(1,:)
!       print*, ' '
!       print*, pres_hedin(10,:)
!       stop
       do i=1,klon
         ang0=acos(mu0(i))*180./3.1415
         do jmu=1,musize
            ang_hedin=acos(mu_hedin(jmu))*180./3.1415
         if (ang_hedin.le.ang0 ) then   
!          if (mu_hedin(jmu) .le. mu0(i) ) then  
              mu_ok = jmu          
         endif
         enddo 
!       print*, ang_hedin, ang0, mu_ok
!       STOP

        do k=1,klev
         z_ok=2         
         do jz=2,zsize-1
           if (pres_hedin(mu_ok,jz).ge.pression(i,k)) then
                 z_ok = jz+1
           endif     
          enddo

        if (pres_hedin(mu_ok,zsize).ge.pression(i,k)) then
! avoid the extrapolation... Problem when p < 1.E-6 Pa...
         factp(k) = 1.
        else
         factp(k) = (log10(pression(i,k))-log10(pres_hedin(mu_ok,z_ok-1)))/(log10(pres_hedin(mu_ok,z_ok))-log10(pres_hedin(mu_ok,z_ok-1)))
        endif

        ovmr_gcm(i,k) = (10.**(log10(o_hedin(mu_ok,z_ok))*factp(k)+log10(o_hedin(mu_ok,z_ok-1))*(1.-factp(k))))
        n2vmr_gcm(i,k) = 10.**(log10(n2_hedin(mu_ok,z_ok))*factp(k)+log10(n2_hedin(mu_ok,z_ok-1))*(1.-factp(k)))
        nvmr_gcm(i,k) = 10.**(log10(n_hedin(mu_ok,z_ok))*factp(k)+log10(n_hedin(mu_ok,z_ok-1))*(1.-factp(k)))
        covmr_gcm(i,k) = 10.**(log10(co_hedin(mu_ok,z_ok))*factp(k)+log10(co_hedin(mu_ok,z_ok-1))*(1.-factp(k)))
        co2vmr_gcm(i,k) =1.- (ovmr_gcm(i,k)+n2vmr_gcm(i,k)+covmr_gcm(i,k))
!        co2vmr_gcm(i,k) =1.- (ovmr_gcm(i,k)+n2vmr_gcm(i,k)+covmr_gcm(i,k)+nvmr_gcm(i,k))

        enddo   
       enddo
!      stop
        

END SUBROUTINE compo_hedin83_mod


  end module compo_hedin83_mod2