Changeset 655 for trunk/LMDZ.MARS/libf/aeronomars

Timestamp:

May 11, 2012, 2:23:57 PM (13 years ago)

Author:

flefevre

Message:

• aeronomars/inichim_newstart.F : initialization of chemistry now handles

nitrogen species and ions (FGG, FL)

• aeronomars/conc.h : cleaning of obsolete variables (FL)
• aeronomars/chimiedata.h : cleaning of obsolete variables (FL)

Location:

trunk/LMDZ.MARS/libf/aeronomars

Files:

• chimiedata.h (modified) (4 diffs)
• conc.h (modified) (1 diff)
• inichim_newstart.F90 (modified) (12 diffs)

trunk/LMDZ.MARS/libf/aeronomars/chimiedata.h

@@ -1 @@
-! ------------------------------------
-! ---  Data for chemical routines  ---
-! ------------------------------------
-!
+! -----------------------------------------
+! ---  Data for photolysis lookup table ---
+! -----------------------------------------
+
       integer nd, nozo, nr, nsza, ntemp, ntau
-!
+
       parameter (nd    = 29)
       parameter (nozo  = 7)
@@ -11 +11 @@
       parameter (ntemp = 4)
       parameter (ntau  = 8)
-!
+
       real kb
       parameter (kb = 1.3806e-23)
-!
+
       common/chimiedata/jphot,colairtab,table_ozo
-!
+
       real jphot(ntemp,nsza,0:200,nozo,ntau,nd)
       real colairtab(0:200)
@@ -22 +22 @@
       real table_ozo(nozo)
       real tautab(ntau)
-!
+
       data szatab/0.,  5., 10., 15., 20., 25.,                          &
      &            30., 35., 40., 45., 50., 55.,                         &
@@ -28 +28 @@
      &            84., 86., 88., 90., 91., 92.,                         &
      &            93., 94., 95./
-!
+
       data tautab/0., 0.2, 0.4, 0.6, 0.8, 1., 2., 4./
-!
-      integer       ncomp
-      parameter    (ncomp = 15)
-      character*10  nomchem(ncomp)
-      data          nomchem/"co2","co","o","o1d","o2","o3","h","h2",    &
-     &                      "oh","ho2","h2o2", "ch4", "n2", "ar", "h2o"/
-!
-      real mmolchem(ncomp)
-      data mmolchem/44.,28.,16.,16.,32.,48., 1.,                        &
-     &               2.,17.,33.,34.,16.,28.,40.,18./

trunk/LMDZ.MARS/libf/aeronomars/conc.h

@@ -4 +4 @@
         
 !**********************************************************************
-
-      integer ncomptot
-
-      parameter (ncomptot=ncomp+1)
 
       common/conc/mmean,Akknew,cpnew,rnew

trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90

@@ -58 +58 @@
 
       integer, parameter         :: nalt = 252    ! number of altitudes in the initialization files
-      integer, parameter         :: nspe = 14     ! number of species in the initialization files
-      integer, dimension(nspe)   :: niq           ! local index of species in initialization files
+      integer                    :: nspe          ! number of species in the initialization files
+      integer, allocatable       :: niq(:)        ! local index of species in initialization files
       real, dimension(nalt)      :: tinit, zzfile ! temperature in initialization files
       real, dimension(nalt)      :: pinit         ! pressure in initialization files
       real, dimension(nalt)      :: densinit      ! total number density in initialization files
-      real, dimension(nalt,nspe) :: vmrinit       ! mixing ratios in initialization files
-      real, dimension(nspe)      :: vmrint        ! mixing ratio interpolated onto the gcm vertical grid
+      real, allocatable          :: vmrinit(:,:)  ! mixing ratios in initialization files
+      real, allocatable          :: vmrint(:)     ! mixing ratio interpolated onto the gcm vertical grid
       real                       :: vmr
 
       character(len=20)          :: txt           ! to store some text
+      logical                    :: flagnitro     ! checks if N species present
 
 ! 1. identify tracers by their names: (and set corresponding values of mmol)
@@ -107 +108 @@
       igcm_n2plus      = 0
       igcm_hplus       = 0
+      igcm_hco2plus    = 0
       igcm_elec        = 0
 
@@ -342 +344 @@
            nbqchem = nbqchem + 1
         end if
+        if (noms(iq) == "hco2plus") then
+           igcm_hco2plus = iq
+           mmol(igcm_hco2plus) = 45.
+           count = count + 1
+           nbqchem = nbqchem + 1
+        end if
         if (noms(iq) == "elec") then
            igcm_elec = iq
@@ -373 +381 @@
          end do
       end if
+
+! 1.7 check if nitrogen species are present:
+
+      if(igcm_no == 0) then
+         !check that no N species is in traceur.def
+         if(igcm_n /= 0 .or. igcm_no2 /= 0 .or. igcm_n2d /= 0) then
+            write(*,*)'inichim_newstart error:'
+            write(*,*)'N, NO2 and/or N2D are in traceur.def, but not NO'
+            write(*,*)'stop'
+            stop
+         endif
+         flagnitro = .false.
+         nspe = 14
+      else
+         !check that all N species are in traceur.def
+         if(igcm_n == 0 .or. igcm_no2 == 0 .or. igcm_n2d == 0) then
+            write(*,*)'inichim_newstart error:'
+            write(*,*)'if NO is in traceur.def, N, NO2 and N2D must also be'
+            write(*,*)'stop'
+            stop
+         endif
+         flagnitro = .true.
+         nspe = 18
+      endif
+
+! 1.8 allocate arrays
+
+      allocate(niq(nspe))
+      allocate(vmrinit(nalt,nspe))
+      allocate(vmrint(nspe))
 
 ! 2. load in chemistry data for initialization:
@@ -395 +433 @@
 !    6: o1d
 !    7: o3
+!
+! order of nitrogen species in initialization file:
+!
+!    1: n
+!    2: no
+!    3: no2
+!    4: n2d
 
 ! major species:
@@ -416 +461 @@
       niq(14) = igcm_o3
 
+! nitrogen species:
+
+      if (flagnitro) then
+         niq(15) = igcm_n
+         niq(16) = igcm_no
+         niq(17) = igcm_no2
+         niq(18) = igcm_n2d         
+      end if
+
 ! 2.1 open initialization files
 
-      open(210, iostat=ierr,file=trim(datafile)// &
-                            '/atmosfera_LMD_may.dat')
+      open(210, iostat=ierr,file=trim(datafile)//'/atmosfera_LMD_may.dat')
       if (ierr /= 0) then
          write(*,*)'Error : cannot open file atmosfera_LMD_may.dat '
-         write(*,*)'(in aeronomars/inichim.F)'
+         write(*,*)'(in aeronomars/inichim_newstart.F)'
          write(*,*)'It should be in :', trim(datafile),'/'
          write(*,*)'1) You can change this path in callphys.def with'
@@ -431 +484 @@
          stop
       end if
-      open(220, iostat=ierr,file=trim(datafile)// &
-                            '/atmosfera_LMD_min.dat')
+      open(220, iostat=ierr,file=trim(datafile)//'/atmosfera_LMD_min.dat')
       if (ierr /= 0) then
          write(*,*)'Error : cannot open file atmosfera_LMD_min.dat '
-         write(*,*)'(in aeronomars/inichim.F)'
+         write(*,*)'(in aeronomars/inichim_newstart.F)'
          write(*,*)'It should be in :', trim(datafile),'/'
          write(*,*)'1) You can change this path in callphys.def with'
@@ -444 +496 @@
          stop
       end if
+      if(flagnitro) then
+         open(230, iostat=ierr,file=trim(datafile)//'/atmosfera_LMD_nitr.dat')
+         if (ierr.ne.0) then
+            write(*,*)'Error : cannot open file atmosfera_LMD_nitr.dat '
+            write(*,*)'(in aeronomars/inichim_newstart.F)'
+            write(*,*)'It should be in :', datafile
+            write(*,*)'1) You can change this directory address in '
+            write(*,*)'   file phymars/datafile.h'
+            write(*,*)'2) If necessary atmosfera_LMD_nitr.dat (and others)'
+            write(*,*)'   can be obtained online on:'
+            write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
+            STOP
+         endif
+      endif   ! Of if(flagnitro)
 
 ! 2.2 read initialization files
@@ -466 +532 @@
       close(220)
 
+! nitrogen species
+
+      if (flagnitro) then
+         read(230,*)
+         do l = 1,nalt
+            read(230,*) dummy, (vmrinit(l,n), n = 15,18)
+         end do
+         close(230)
+      end if
+      
 ! 3. initialization of tracers
 
@@ -488 +564 @@
 !                             - with our without h2o 
 
-               if (flagthermo == 0 .or. (flagthermo == 1 .and. exp(pgcm) < 1.e-3)) then
+               if (flagthermo == 0 .or. (flagthermo == 1 .and. exp(pgcm) < 0.1)) then
                   do n = 1,nspe
                      iq = niq(n)
@@ -501 +577 @@
       end do
 
-      ! set surface values of chemistry tracers to zero
+! set surface values of chemistry tracers to zero
+
       if (flagthermo == 0) then
-        ! NB: no problem for "surface water vapour" tracer which is always 0
-        do n=1,nspe
-          iq=niq(n)
-          qsurf(1:ngridmx,iq)=0
-        enddo
-      endif
-
+         ! NB: no problem for "surface water vapour" tracer which is always 0
+         do n = 1,nspe
+            iq = niq(n)
+            qsurf(1:ngridmx,iq) = 0.
+         end do
+      end if
 
 ! 3.3 initialization of tracers not contained in the initialization files
@@ -525 +601 @@
          end do
          ! set surface value to zero
-         qsurf(1:ngridmx,igcm_ch4)=0
+         qsurf(1:ngridmx,igcm_ch4) = 0.
       end if
 
+! ions: 0
+
+      if (igcm_co2plus /= 0) then
+         !check that all required ions are in traceur.def
+         if (igcm_o2plus == 0 .or. igcm_oplus == 0 .or. igcm_coplus == 0          &
+              .or. igcm_cplus == 0 .or. igcm_nplus == 0 .or. igcm_noplus == 0    & 
+              .or. igcm_n2plus == 0 .or. igcm_hplus == 0 .or. igcm_hco2plus == 0 &
+              .or. igcm_elec == 0) then
+            write(*,*)'inichim_newstart error:'
+            write(*,*)'if co2plus is in traceur.def, all other ions must also be'
+            write(*,*)'o2plus, oplus, coplus, cplus, nplus, noplus, n2plus'
+            write(*,*)'hplus, hco2plus and elec'
+            write(*,*)'stop'
+            stop
+         end if
+
+         do i = 1,iip1
+            do j = 1,jjp1
+               do l = 1,llm
+                  ! all ions to 0     
+                  pq(i,j,l,igcm_co2plus)  = 0.
+                  pq(i,j,l,igcm_o2plus)   = 0.
+                  pq(i,j,l,igcm_oplus)    = 0.
+                  pq(i,j,l,igcm_coplus)   = 0.
+                  pq(i,j,l,igcm_cplus)    = 0.
+                  pq(i,j,l,igcm_nplus)    = 0.
+                  pq(i,j,l,igcm_noplus)   = 0.
+                  pq(i,j,l,igcm_n2plus)   = 0.
+                  pq(i,j,l,igcm_hplus)    = 0.
+                  pq(i,j,l,igcm_hco2plus) = 0.
+                  pq(i,j,l,igcm_elec)     = 0.
+               end do
+            end do
+         end do
+
+         ! surface value to 0
+
+         qsurf(1:ngridmx,igcm_co2plus)  = 0.
+         qsurf(1:ngridmx,igcm_o2plus)   = 0.
+         qsurf(1:ngridmx,igcm_oplus)    = 0.
+         qsurf(1:ngridmx,igcm_coplus)   = 0.
+         qsurf(1:ngridmx,igcm_cplus)    = 0.
+         qsurf(1:ngridmx,igcm_nplus)    = 0.
+         qsurf(1:ngridmx,igcm_noplus)   = 0.
+         qsurf(1:ngridmx,igcm_n2plus)   = 0.
+         qsurf(1:ngridmx,igcm_hplus)    = 0.
+         qsurf(1:ngridmx,igcm_hco2plus) = 0.
+         qsurf(1:ngridmx,igcm_elec)     = 0.
+
+      else
+
+         if (igcm_o2plus /= 0 .or. igcm_oplus /= 0 .or. igcm_coplus /= 0          &
+              .or. igcm_cplus /= 0 .or. igcm_nplus /= 0 .or. igcm_noplus /= 0    & 
+              .or. igcm_n2plus /= 0 .or. igcm_hplus /= 0 .or. igcm_hco2plus /= 0 &
+              .or. igcm_elec /= 0) then
+            write(*,*)'inichim_newstart error:'
+            write(*,*)'some ions are in traceur.def, but not co2plus'
+            write(*,*)'stop'
+            stop
+         end if
+      end if    ! of if(igcm_co2 /= 0)
+      
+      ! deallocations
+
+      deallocate(niq)
+      deallocate(vmrinit)
+      deallocate(vmrint)
+
       end