source: trunk/mesoscale/LMDZ.MARS/libf_gcm/phymars/dustopacity.F @ 57

      SUBROUTINE dustopacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,pq,
     $    tauref,tau,aerosol)
                                                   
       IMPLICIT NONE
c=======================================================================
c   subject:
c   --------
c   Computing aerosol optical depth (dust opacity) 
c   In each layers
c
c   author: F.Forget 
c   ------
c   update F. Montmessin (water ice scheme) 
c      and S. Lebonnois (12/06/2003) compatibility dust/ice/chemistry
c   
c   input:
c   ----- 
c   ngrid             Number of gridpoint of horizontal grid
c   nlayer            Number of layer
c   nq                Number of tracer
c   ls                Solar longitude (Ls) , radian
c   pplay,pplev       pressure (Pa) in the middle and boundary of each layer
c   pq                Dust mixing ratio (used if tracer =T and active=T).
c
c   output:
c   -------
c   tauref       Prescribed mean column optical depth at 700 Pa 
c   tau          Column total visible dust optical depth at each point
c   aerosol      aerosol(ig,l,1) is the dust optical
c                depth in layer l, grid point ig

c
c=======================================================================
#include "dimensions.h"
#include "dimphys.h"
#include "callkeys.h"
#include "comcstfi.h"
#include "comgeomfi.h"
#include "dimradmars.h"
#include "yomaer.h"
#include "tracer.h"
#include "planete.h"

c-----------------------------------------------------------------------
c
c    Declarations :
c    --------------
c
c    Input/Output
c    ------------
      INTEGER ngrid,nlayer,nq

      REAL ls,zday,expfactor    
      REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
      REAL pq(ngrid,nlayer,nq)
      REAL tauref(ngrid), tau(ngrid,naerkind)
      REAL aerosol(ngrid,nlayer,naerkind)
c
c    Local variables :
c    -----------------
      INTEGER l,ig,iq
      real topdust(ngridmx)
      real zlsconst, zp
      real taueq,tauS,tauN
      real r0,reff,coefsize
c
c   local saved variables
c   ---------------------

      REAL topdust0(ngridmx) 
      SAVE topdust0

      LOGICAL firstcall
      DATA firstcall/.true./
      SAVE firstcall


c----------------------------------------------------------------------

c     Initialisation
c     --------------

      IF (firstcall) THEN

c        altitude of the top of the aerosol layer (km) at Ls=2.76rad:
c        in the Viking year scenario
         DO ig=1,ngrid
             topdust0(ig)=60. -22.*SIN(lati(ig))**2
         END DO
         firstcall=.false.
      END IF

c  -------------------------------------------------------------
c     1) Prescribed dust  (if tracer=F or active=F)
c  -------------------------------------------------------------
      IF ((.not.tracer) .or. (.not.active)) THEN 

c       Vertical column optical depth at 700.Pa 
c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        IF(iaervar.eq.1) THEN
           do ig=1, ngridmx
            tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste as read in starfi
          end do
        ELSE IF (iaervar.eq.2) THEN   ! << "Viking" Scenario>>

          tauref(1) = 0.7+.3*cos(ls+80.*pi/180.) ! like seen by VL1
          do ig=2,ngrid
            tauref(ig) = tauref(1)
          end do

        ELSE IF (iaervar.eq.3) THEN  ! << "MGS" scenario >>

           taueq= 0.2 +(0.5-0.2) *(cos(0.5*(ls-4.363)))**14
           tauS= 0.1 +(0.5-0.1)  *(cos(0.5*(ls-4.363)))**14
           tauN = 0.1
c          if (peri_day.eq.150) then
c            tauS=0.1
c            tauN=0.1 +(0.5-0.1)  *(cos(0.5*(ls+pi-4.363)))**14
c            taueq= 0.2 +(0.5-0.2) *(cos(0.5*(ls+pi-4.363)))**14
c           endif
             
           do ig=1,ngrid/2  ! Northern hemisphere
             tauref(ig)= tauN +
     &       (taueq-tauN)*0.5*(1+tanh((45-lati(ig)*180./pi)*6/60))
           end do
           do ig=ngrid/2+1, ngridmx  ! Southern hemisphere
             tauref(ig)= tauS +
     &       (taueq-tauS)*0.5*(1+tanh((45+lati(ig)*180./pi)*6/60))
           end do

        ELSE IF (iaervar.eq.4) THEN  ! << "TES scenario >>
           call readtesassim(ngrid,nlayer,zday,pplev,tauref)

        ELSE IF (iaervar.eq.5) THEN   ! << Escalier Scenario>>
c         tauref(1) = 0.2
c         if ((ls.ge.210.*pi/180.).and.(ls.le.330.*pi/180.))
c    &                              tauref(1) = 2.5
          tauref(1) = 2.5
          if ((ls.ge.30.*pi/180.).and.(ls.le.150.*pi/180.))
     &                              tauref(1) = .2

          do ig=2,ngrid
            tauref(ig) = tauref(1)
          end do
        ELSE
          stop 'problem with iaervar in dustopacity.F'
        ENDIF

c       Altitude of the top of the dust layer
c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        zlsconst=SIN(ls-2.76)
        if (iddist.eq.1) then
          do ig=1,ngrid
             topdust(ig)=topdustref         ! constant dust layer top
          end do

        else if (iddist.eq.2) then          ! "Viking" scenario
          do ig=1,ngrid
            topdust(ig)=topdust0(ig)+18.*zlsconst
          end do

        else if(iddist.eq.3) then         !"MGS" scenario
          do ig=1,ngrid
            topdust(ig)=60.+18.*zlsconst
     &                -(32+18*zlsconst)*sin(lati(ig))**4
     &                 - 8*zlsconst*(sin(lati(ig)))**5
          end do
        endif


c       Optical depth in each layer :
c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        if(iddist.ge.1) then
         expfactor=0.
          DO l=1,nlayer
             DO ig=1,ngrid
             if(pplay(ig,l).gt.700.
     $                        /(988.**(topdust(ig)/70.))) then
                zp=(700./pplay(ig,l))**(70./topdust(ig))
                 expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
              else    
               expfactor=1.e-3
              endif       
                aerosol(ig,l,1)= tauref(ig)/700. *
     s           (pplev(ig,l)-pplev(ig,l+1))
     &           *expfactor
c     s           *max( exp(.007*(1.-max(zp,1.))) , 1.E-3 )
             ENDDO
          ENDDO
c     changement dans le calcul de la distribution verticale          
c     dans le cas des scenarios de poussieres assimiles
c        if (iaervar.eq.4) THEN  ! TES
c        call zerophys(ngrid*naerkind,tau)
c
c        do l=1,nlayer
c           do ig=1,ngrid
c                tau(ig,1)=tau(ig,1)+ aerosol(ig,l,1)
c           end do
c        end do
c        do l=1,nlayer
c           do ig=1,ngrid
c               aerosol(ig,l,1)=aerosol(ig,l,1)*tauref(ig)/tau(ig,1)
c     $     *(pplev(ig,1)/700)
c           end do
c        end do
c        endif
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc       
        else if(iddist.eq.0) then   
c         old dust vertical distribution function (pollack90)
          DO l=1,nlayer
             DO ig=1,ngrid
                zp=700./pplay(ig,l)
                aerosol(ig,l,1)= tauref(ig)/700. *
     s           (pplev(ig,l)-pplev(ig,l+1))
     s           *max( exp(.03*(1.-max(zp,1.))) , 1.E-3 )
             ENDDO
          ENDDO
        end if

c  ---------------------------------------------------------------------
c     2) Transported radiatively active dust  (if tracer=T and active=T)
c ----------------------------------------------------------------------
      ELSE  IF ((tracer) .and. (active)) THEN 
c     The dust opacity is computed from q

c     a) "doubleq" technique (transport of mass and number mixing ratio)
c        ~~~~~~~~~~~~~~~~~~~
       if(doubleq) then 

         call zerophys(ngrid*nlayer*naerkind,aerosol)

c        Computing effective radius :
         do  l=1,nlayer
           do ig=1, ngrid
              r0=
     &        (r3n_q*pq(ig,l,1)/max(pq(ig,l,2),0.01))**(1./3.)
              r0=min(max(r0,1.e-10),500.e-6)
              reff= ref_r0 * r0
cc           If  reff is small, the transported dust mean Qext
c            is reduced from the reference dust Qext by a factor "coefsize"

             coefsize=min(max(2.52e6*reff-0.043 ,0.)   ,1.)

cc              It is added 1.e-8 to pq to avoid low

                aerosol(ig,l,1)=aerosol(ig,l,1)+ 1.E-8 +
     &         ( 0.75*Qext(1)*coefsize/(rho_dust*reff))
     &          * (pq(ig,l,1))*
     &          (pplev(ig,l)-pplev(ig,l+1))/g
           end do
         end do
         call zerophys(ngrid,tauref)

c     b) Size bin technique (each aerosol can contribute to opacity))
c        ~~~~~~~~~~~~~~~~~~
       else
c        The dust opacity is computed from q
         call zerophys(ngrid*nlayer*naerkind,aerosol)
         do iq=1,dustbin
           do l=1,nlayer
              do ig=1,ngrid
cc               qextrhor(iq) is  (3/4)*Qext/(rho*reff)
cc               It is added 1.e-8 to pq to avoid low
                 aerosol(ig,l,1)=aerosol(ig,l,1)+
     &           qextrhor(iq)* (pq(ig,l,iq) + 1.e-8)*
     &           (pplev(ig,l)-pplev(ig,l+1))/g
              end do
           end do
         end do
         call zerophys(ngrid,tauref)
       end if  ! (doubleq)
      END IF   ! (dust scenario)

c --------------------------------------------------------------------------
c Column integrated visible optical depth in each point (used for diagnostic)
c --------------------------------------------------------------------------
      call zerophys(ngrid*naerkind,tau)
      do l=1,nlayer
         do ig=1,ngrid
               tau(ig,1)=tau(ig,1)+ aerosol(ig,l,1)
         end do
      end do

      return
      end