Ignore:
Timestamp:
Oct 28, 2011, 5:00:48 PM (13 years ago)
Author:
aslmd
Message:

LMDZ.MARS:

28/10/11 == FL + AS

ADDED THE NEW FRAMEWORK FOR PHOTOCHEMISTRY
This is not the last version. A new rewritten version of calchim.F [using LAPACK] is yet to be committed.
--> A new version for photochemistry routines : *_B.F no longer exists, new routines in aeronomars
D 333 libf/aeronomars/photochemist_B.F
D 333 libf/aeronomars/init_chimie_B.F
A 0 libf/aeronomars/read_phototable.F
M 333 libf/aeronomars/calchim.F
A 0 libf/aeronomars/photochemistry.F
M 333 libf/aeronomars/chimiedata.h
--> Changed calls to calchim and watercloud [surfdust and surfice needed] in physiq.F
--> Left commented code for outputs in physiq.F [search for FL]
--> Added settings which works for 35 levels in inidissip.F according to FL. Commented for the moment.
--> Checked compilation and run, looks fine. Note that 'jmars.20101220' is needed.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/aeronomars/calchim.F

    r38 r334  
    11      subroutine calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
    2      $                   zzlay,zday,pq,pdq,rice,
    3      &                   dqchim,dqschim,dqcloud,dqscloud)
     2     $                   zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud,
     3     $                   dqscloud,tauref,co2ice,
     4     $                   pu,pdu,pv,pdv,surfdust,surfice)
    45c
    56      implicit none
    … …  
    1920c    update 03/05/2005 cosmetic changes (Franck Lefevre)
    2021c    update sept. 2008 identify tracers by their names (Ehouarn Millour)
     22c    update 17/03/2011 synchronize with latest version of chemistry (Franck Lefevre)
    2123c
    2224c   Arguments:
    … …  
    3032c    pt(ngridmx,nlayermx)       Temperature (K)
    3133c    pdt(ngridmx,nlayermx)      Temperature tendency (K)
     34c    pu(ngridmx,nlayermx)       u component of the wind (ms-1)
     35c    pdu(ngridmx,nlayermx)      u component tendency (K)
     36c    pv(ngridmx,nlayermx)       v component of the wind (ms-1)
     37c    pdv(ngridmx,nlayermx)      v component tendency (K)
    3238c    dist_sol                   distance of the sun (AU)
    3339c    mu0(ngridmx)               cos of solar zenith angle (=1 when sun at zenith)
    3440c    pq(ngridmx,nlayermx,nqmx)  Advected fields, ie chemical species here
    3541c    pdq(ngridmx,nlayermx,nqmx) Previous tendencies on pq
    36 c    rice(ngridmx,nlayermx)     Estimated ice crystal radius (m)
     42c    tauref(ngridmx)            Optical depth at 7 hPa
     43c    co2ice(ngridmx)            co2 ice surface layer (kg.m-2)
     44c    surfdust(ngridmx,nlayermx) dust surface area (micron2/cm3)
     45c    surfice(ngridmx,nlayermx)  ice surface area (micron2/cm3)
    3746c
    3847c  Output:
    … …  
    6473      real    zzlay(ngridmx,nlayermx)    ! pressure at the middle of the layers
    6574      real    pplev(ngridmx,nlayermx+1)  ! intermediate pressure levels
     75      real    zzlev(ngridmx,nlayermx+1)  ! altitude at layer boundaries
    6676      real    pt(ngridmx,nlayermx)       ! temperature
    6777      real    pdt(ngridmx,nlayermx)      ! temperature tendency
     78      real    pu(ngridmx,nlayermx)       ! u component of the wind (m.s-1)
     79      real    pdu(ngridmx,nlayermx)      ! u component tendency
     80      real    pv(ngridmx,nlayermx)       ! v component of the wind (m.s-1)
     81      real    pdv(ngridmx,nlayermx)      ! v component tendency
    6882      real    dist_sol                   ! distance of the sun (AU)
    69       real    mu0(ngridmx)       ! cos of solar zenith angle (=1 when sun at zenith)
     83      real    mu0(ngridmx)               ! cos of solar zenith angle (=1 when sun at zenith)
    7084      real    pq(ngridmx,nlayermx,nqmx)  ! tracers mass mixing ratio
    7185      real    pdq(ngridmx,nlayermx,nqmx) ! previous tendencies
    7286      real    zday                       ! date (time since Ls=0, in martian days)
    73       real    rice(ngridmx,nlayermx)     ! Estimated ice crystal radius (m)
     87      real    tauref(ngridmx)            ! optical depth at 7 hPa
     88      real    co2ice(ngridmx)            ! co2 ice surface layer (kg.m-2)
     89      real    surfdust(ngridmx,nlayermx) ! dust surface area (micron2/cm3)
     90      real    surfice(ngridmx,nlayermx)  !  ice surface area (micron2/cm3)
    7491
    7592c   outputs:
    … …  
    87104      integer  ig,l,i,iq
    88105      integer  foundswitch, lswitch
     106      integer  ig_vl1
     107      real latvl1, lonvl1
    89108      real zq(ngridmx,nlayermx,nqmx) ! pq+pdq*ptimestep before chemistry
    90109                                     ! new mole fraction after
    91       real colden(ngridmx,nqmx)      ! Column densities (cm-2)
    92110      real zt(ngridmx,nlayermx)      ! temperature
     111      real zu(ngridmx,nlayermx)      ! u component of the wind
     112      real zv(ngridmx,nlayermx)      ! v component of the wind
     113      real taucol                    ! optical depth at 7 hPa
     114
     115      logical depos                  ! switch for dry deposition
     116      parameter (depos=.false.)
    93117c
    94118c  for each column of atmosphere:
    … …  
    100124      real zycol(nlayermx,nqmx)   !  Composition (mole fractions)
    101125      real szacol                 !  Solar zenith angle
     126      real surfice1d(nlayermx)    !  Ice surface area (cm2/cm3)
     127      real surfdust1d(nlayermx)   !  Dust surface area (cm2/cm3)
    102128      real jo3(nlayermx)          !  Photodissociation rate O3->O1D (s-1)
    103129c
    104130c for output:
    105131c
    106       real zdens3d(ngridmx,nlayermx) !  Density  (cm-3)
    107132      real jo3_3d(ngridmx,nlayermx)  !  Photodissociation rate O3->O1D (s-1)
    108       real surfice(ngridmx,nlayermx) !  Surface of ice particules (um2/cm3)
    109133
    110134      logical output              ! to issue calls to writediagfi and stats
    … …  
    112136
    113137      logical,save :: firstcall=.true.
    114       integer,save :: nbq  ! number of tracers used in the chemistry
     138      integer,save :: nbq       ! number of tracers used in the chemistry
    115139      integer,save :: niq(nqmx) ! array storing the indexes of the tracers
    116140
    117141! index of tracers:
     142
    118143      integer,save :: i_co2=0
    119144      integer,save :: i_co=0
    … …  
    127152      integer,save :: i_ho2=0
    128153      integer,save :: i_h2o2=0
     154      integer,save :: i_ch4=0
    129155      integer,save :: i_n2=0
    130156      integer,save :: i_ar=0
    131157      integer,save :: i_ice=0 ! water ice
    132158      integer,save :: i_h2o=0 ! water vapour
    133 
    134 c
    135 c  scheme A: 1 ; scheme B: 2
    136 c
    137       integer,parameter :: scheme=2
    138159c
    139160c=======================================================================
    … …  
    144165c
    145166         if (photochem) then
    146             print*,'calchim: INIT CHEMISTRY'
    147             if (scheme  .eq.  1) then
    148                print*,'calchim: Scheme A : A METTRE A JOUR !!'
    149                stop
    150 c              call init_chimie_A
    151             else
    152                print*,'calchim: Scheme B'
    153                call init_chimie_B
    154             end if
     167            print*,'calchim: Read photolysis lookup table'
     168            call read_phototable
    155169         end if
    156170
    … …  
    245259           niq(nbq)=i_h2o2
    246260         endif
     261         i_ch4=igcm_ch4
     262         if (i_ch4.eq.0) then
     263           write(*,*) "calchim: Error; no CH4 tracer !!!"
     264           stop
     265         else
     266           nbq=nbq+1
     267           niq(nbq)=i_ch4
     268         endif
    247269         i_n2=igcm_n2
    248270         if (i_n2.eq.0) then
    … …  
    278300         endif
    279301
    280          write(*,*) 'calchim: found nbq=',nbq,' tracers'
    281          write(*,*) '         i_co2=',i_co2
    282          write(*,*) '         i_co=',i_co
    283          write(*,*) '         i_o=',i_o
    284          write(*,*) '         i_o1d=',i_o1d
    285          write(*,*) '         i_o2=',i_o2
    286          write(*,*) '         i_o3=',i_o3
    287          write(*,*) '         i_h=',i_h
    288          write(*,*) '         i_h2=',i_h2
    289          write(*,*) '         i_oh=',i_oh
    290          write(*,*) '         i_ho2=',i_ho2
    291          write(*,*) '         i_h2o2=',i_h2o2
    292          write(*,*) '         i_n2=',i_n2
    293          write(*,*) '         i_ar=',i_ar
    294          write(*,*) '         i_ice=',i_ice
    295          write(*,*) '         i_h2o=',i_h2o
    296 !         write(*,*) '         niq(:)=',niq
    297 !         write(*,*) '     nqchem_min,nqmx=',nqchem_min,nqmx
     302         write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
     303         write(*,*) '               i_co2  = ',i_co2
     304         write(*,*) '               i_co   = ',i_co
     305         write(*,*) '               i_o    = ',i_o
     306         write(*,*) '               i_o1d  = ',i_o1d
     307         write(*,*) '               i_o2   = ',i_o2
     308         write(*,*) '               i_o3   = ',i_o3
     309         write(*,*) '               i_h    = ',i_h
     310         write(*,*) '               i_h2   = ',i_h2
     311         write(*,*) '               i_oh   = ',i_oh
     312         write(*,*) '               i_ho2  = ',i_ho2
     313         write(*,*) '               i_h2o2 = ',i_h2o2
     314         write(*,*) '               i_ch4  = ',i_ch4
     315         write(*,*) '               i_n2   = ',i_n2
     316         write(*,*) '               i_ar   = ',i_ar
     317         write(*,*) '               i_ice  = ',i_ice
     318         write(*,*) '               i_h2o  = ',i_h2o
    298319         
    299320         firstcall = .false.
    … …  
    302323! Initialize output tendencies to zero (to handle case of tracers which
    303324! are not used in the chemistry (e.g. dust))
    304       dqchim(:,:,:)=0
    305       dqschim(:,:)=0
    306 
    307 
     325
     326      dqchim(:,:,:) = 0
     327      dqschim(:,:)  = 0
     328c
     329!     latvl1= 22.27
     330!     lonvl1= -47.94
     331!     ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.)  -2 )*iim +
     332!    $        int(1.5+(lonvl1+180)*iim/360.)
    308333c
    309334c=======================================================================
    … …  
    317342         foundswitch = 0
    318343         do l = 1,nlayermx
    319             zt(ig,l)=pt(ig,l)+pdt(ig,l)*ptimestep
    320             do i=1,nbq
    321               iq=niq(i) ! get tracer index
     344            do i = 1,nbq
     345              iq = niq(i) ! get tracer index
    322346              zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
    323               zycol(l,iq) = zq(ig,l,iq) * mmean(ig,l)/mmol(iq)
    324             enddo
     347              zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
     348            end do
     349            zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
     350            zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
     351            zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
     352c
    325353            zpress(l) = pplay(ig,l)/100.
    326354            ztemp(l)  = zt(ig,l)
    … …  
    328356            zlocal(l) = zzlay(ig,l)/1000.
    329357c
    330 c     search for switch index between regions
     358c           surfdust1d and surfice1d in cm2/cm3
     359c
     360            surfdust1d(l) = surfdust(ig,l)*1.e-8
     361            surfice1d(l)  = surfice(ig,l)*1.e-8
     362c
     363c           search for switch index between regions
    331364c
    332365            if (photochem .and. thermochem) then
    … …  
    340373            end if
    341374            if (.not.  thermochem) then
    342                lswitch = min(33,nlayermx+1)
     375               lswitch = min(50,nlayermx+1)    ! FL (original value: 33)
    343376            end if
    344377c
    345 c     ice surface area  in microns^2/cm^3
    346 c
    347 c     = 4 pi r^2 * [ zq * mugaz/NA / (rhoice*4/3 pi r^3) ] *zdens
    348 c     = 3/r * [ zq * mugaz/NA / rhoice ] *zdens
    349 c     with r in microns, rhoice = 0.92e-12 g microns^-3 and zdens in cm^-3
    350 c
    351             if (water) then
    352                zycol(l,i_ice) = (3.e-6/rice(ig,l))*zq(ig,l,i_ice)
    353      $                           *(mugaz/6.022e23)*zdens(l)/0.92e-12
    354 c              write(*,*) "rice=",rice(ig,l)," m / zdens=",zdens(l),
    355 c    $        " cm-3 / icesurf=",zycol(l,nqmx-1)," microns^2/cm^3"
    356                surfice(ig,l) = zycol(l,i_ice)
    357             end if
    358 c
    359378         end do ! of do l=1,nlayermx
    360379c
    361380         szacol = acos(mu0(ig))*180./pi
     381         taucol = tauref(ig)
    362382c
    363383c=======================================================================
    … …  
    365385c=======================================================================
    366386c
    367         if (photochem) then
    368            if (scheme .eq. 1) then
    369               print*,'Scheme A : A METTRE A JOUR !!'
    370 c             call photochemist_A(zycol,szacol,ptimestep,
    371 c    $                            zpress,ztemp,zdens,dist_sol)
    372            else
    373               call photochemist_B(lswitch,zycol,szacol,ptimestep,
    374      $                            zpress,ztemp,zdens,dist_sol,jo3)
    375            end if
    376         end if
    377 
    378         if (thermochem) then
    379            call chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress,
    380      $                      zlocal,szacol,ptimestep,zday)
    381         end if
     387c        chemistry in lower atmosphere
     388c
     389         if (photochem) then
     390            call photochemistry(lswitch,zycol,szacol,ptimestep,
     391     $                          zpress,ztemp,zdens,dist_sol,
     392     $                          surfdust1d,surfice1d,jo3,taucol)
     393         end if
     394c
     395c        chemistry in upper atmosphere
     396c
     397         if (thermochem) then
     398            call chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress,
     399     $                       zlocal,szacol,ptimestep,zday)
     400         end if
     401c
     402c        dry deposition
     403c
     404         if (depos) then
     405            call deposition(ig, ig_vl1, pplay, pplev, zzlay, zzlev,
     406     $                      zu, zv, zt, zycol, ptimestep, co2ice)
     407         end if
    382408c
    383409c=======================================================================
    … …  
    389415         if (water) then
    390416            do l = 1,nlayermx
    391 !               dqchim(ig,l,nqmx-1) = 0.
    392417               dqchim(ig,l,i_ice) = 0.
    393418            end do
    … …  
    425450     $                                 - zq(ig,l,iq))/ptimestep
    426451                   else if (iq.eq.i_o) then
    427 !                      i_o = iq
    428452                      dqchim(ig,l,iq) = 0.
    429453                   end if
    … …  
    440464          end do ! of do l = 1,nlayermx
    441465c
    442 c     dust: This is now taken care of as a first step at beginning of routine
    443 c
    444 !          if (nqchem_min .gt. 1) then
    445 !             do iq = 1,nqchem_min-1
    446 !                do l = 1,nlayermx
    447 !                   dqchim(ig,l,iq) = 0.
    448 !                end do
    449 !                dqschim(ig,iq) = 0.
    450 !             end do
    451 !          end if
    452 c
    453466c     condensation of h2o2
    454467c
    … …  
    456469     $                 ztemp,zycol,dqcloud,dqscloud)
    457470c
    458 c     for outputs
    459 c
    460           do i=1,nbq
    461             iq=niq(i) ! get tracer index
    462              colden(ig,iq) = 0.
    463              do l = 1,nlayermx
    464 c
    465 c     column density converted in cm-2
    466 c     pplev en pa, mugaz en g.mol-1 et g en m.s-2
    467 c     not for ice
    468 c
    469                 if (iq.ne.i_h2o2) then
    470                   colden(ig,iq) = colden(ig,iq) + zycol(l,iq)
    471      $                         *6.022e22*(pplev(ig,l)-pplev(ig,l+1))
    472      $                         /(mmean(ig,l)*g)
    473                 else   ! for H2O2, remove condensation from zycol
    474                   colden(ig,iq) = colden(ig,iq) + (zycol(l,iq) +
    475      $               dqcloud(ig,l,iq)*ptimestep*mmean(ig,l)/mmol(iq))
    476      $                         *6.022e22*(pplev(ig,l)-pplev(ig,l+1))
    477      $                         /(mmean(ig,l)*g)
    478                 end if
    479 c
    480 c     local densities, for outputs (put in zq)
    481 c     not for ice
    482 c
    483                 zq(ig,l,iq) = zycol(l,iq)*zdens(l)
    484 c                        for H2O2, remove condensation from zycol
    485                 if (iq.eq.i_h2o2) then
    486                    zq(ig,l,iq) = zdens(l)*(zycol(l,iq) +
    487      $               dqcloud(ig,l,iq)*ptimestep*mmean(ig,l)/mmol(iq))
    488                 end if
    489              end do
    490           end do
    491 c
    492471c     density and j(o3->o1d), for outputs
    493472c
    494           zdens3d(ig,1) = zdens(1)
    495           jo3_3d(ig,1) = jo3(1)
    496           do l = 2,nlayermx
    497              zdens3d(ig,l) = zdens(l)
     473          do l = 1,nlayermx
    498474             jo3_3d(ig,l) = jo3(l)
    499475          end do
    … …  
    513489
    514490      if (output) then
    515 
    516491         if (ngridmx .gt. 1) then
    517 c           call writediagfi(ngridmx,'dens','atm dens.','cm-3',3,zdens3d(1,1))
    518 c           call writediagfi(ngridmx,'jo3','j o3->o1d','s-1',3,jo3_3d(1,1))
    519 c           call writediagfi(ngridmx,'sice','ice surf.','um2/cm3',3,surfice(1,1))
    520             do i=1,nbq
    521               iq=niq(i) ! get tracer index
    522                if (iq.ne.i_ice) then
    523                  write(str20(1:20),'(a20)') noms(iq)
    524                  call writediagfi(ngridmx,'n_'//trim(str20),'density',
    525      $                             'cm-3',3,zq(1,1,iq))
    526 c                 call writediagfi(ngridmx,'dqch_'//str5,'density','cm-3',3,dqchim(1,1,iq))
    527 c                 if (noms(iq) .eq. "h2o2" .or. noms(iq) .eq. "h2o") then
    528 c                    call writediagfi(ngridmx,'cl_'//str5,'density','cm-3',3,dqcloud(1,1,iq))
    529 c                 end if
    530                  call writediagfi(ngridmx,'c_'//trim(str20),
    531      $                            'col. dens.','cm-2',2,colden(1,iq))
    532                end if
    533             end do
    534 c
    535             if (callstats) then
    536 c
    537 c              convert to mole.cm-2 for the column densities
    538 c
    539                do i=1,nbq
    540                  iq=niq(i) ! get tracer index
    541                   do ig = 1,ngridmx
    542                      colden(ig,iq) = colden(ig,iq)/6.022e23
    543                   end do   
    544                end do
    545 c
    546 c              call wstats(ngridmx,"jo3","jo3->o1d","s-1",3,jo3_3d)
    547 c
    548                do i=1,nbq
    549                  iq=niq(i) ! get tracer index
    550                   if (iq.ne.i_ice) then
    551                      write(str20(1:20),'(a20)') noms(iq)
    552                      call wstats(ngridmx,"n_"//trim(str20),"density",
    553      &                           "cm-3",3,zq(1,1,iq))
    554                      call wstats(ngridmx,"c_"//trim(str20),"col. dens.",
    555      &                           "mol cm-2",2,colden(1,iq))
    556                   end if
    557                end do ! of i=1,nbq
    558             end if ! of if (callstats)
     492            call writediagfi(ngridmx,'jo3','j o3->o1d',
     493     $                       's-1',3,jo3_3d(1,1))
     494            call wstats(ngridmx,'jo3','j o3->o1d',
     495     $                       's-1',3,jo3_3d(1,1))
    559496         end if ! of if (ngridmx.gt.1)
    560 c
    561497      endif ! of if (output)
    562498c
Note: See TracChangeset for help on using the changeset viewer.