source: trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/inichim_newstart.F90 @ 3546

Last change on this file since 3546 was 3184, checked in by afalco, 11 months ago

Pluto PCM:
Added LMDZ.PLUTO, a copy of the generic model,
cleaned from some unnecessary modules (water, ...)
AF

  • Property svn:executable set to *
File size: 7.9 KB
Line 
1      subroutine inichim_newstart(ngrid, nq, pq, qsurf, ps, &
2                                  flagh2o, flagthermo)
3
4      use tracer_h,          only: noms, mmol
5      use datafile_mod,      only: datadir
6      use comvert_mod,       only: aps,bps
7      use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat, nbp_lev
8
9      implicit none
10
11!=======================================================================
12!
13!  Purpose:
14!  --------
15!
16!  Initialization of the chemistry for use in the building of a new start file
17!  used by program newstart.F
18!  also used by program testphys1d.F
19!
20!  Authors:
21!  -------
22!  Initial version 11/2002 by Sebastien Lebonnois
23!  Revised 07/2003 by Monica Angelats-i-Coll to use input files
24!  Modified 10/2008 Identify tracers by their names Ehouarn Millour
25!  Modified 11/2011 Addition of methane Franck Lefevre
26!  Rewritten 04/2012 Franck Lefevre
27!  Rewritten 03/2021 Yassin Jaziri (Use of #Moderntrac-v1 to init thanks traceur.def)
28!
29!  Arguments:
30!  ----------
31!
32!  pq(nbp_lon+1,nbp_lat,nbp_lev,nq)   Advected fields, ie chemical species here
33!  qsurf(ngrid,nq)                    Amount of tracer on the surface (kg/m2)
34!  ps(nbp_lon+1,nbp_lat)              Surface pressure (Pa)
35!  flagh2o                            flag for initialisation of h2o (1: yes / 0: no)
36!  flagthermo                         flag for initialisation of thermosphere only (1: yes / 0: no)
37!
38!=======================================================================
39
40
41! inputs :
42
43      integer,intent(in) :: ngrid                  ! number of atmospheric columns in the physics
44      integer,intent(in) :: nq                     ! number of tracers
45      real   ,intent(in) :: ps(nbp_lon+1,nbp_lat)  ! surface pressure in the gcm (Pa)
46      integer,intent(in) :: flagh2o                ! flag for h2o initialisation
47      integer,intent(in) :: flagthermo             ! flag for thermosphere initialisation only
48
49! outputs :
50
51      real,intent(out)   :: pq(nbp_lon+1,nbp_lat,nbp_lev,nq) ! advected fields, ie chemical species
52      real,intent(out)   :: qsurf(ngrid,nq)                  ! surface values (kg/m2) of tracers
53
54! local :
55
56      real, allocatable  :: pf(:)                  ! pressure in vmr profile files set in traceur.def
57      real, allocatable  :: qf(:)                  ! vmr      in vmr profile files set in traceur.def
58
59      real    :: pgcm                              ! pressure at each layer in the gcm (Pa)
60      real    :: mmean(nbp_lev)                    ! mean molecular mass (g)
61      real    :: pqx(nbp_lon+1,nbp_lat,nbp_lev,nq) ! tracers (vmr)
62      real    :: qx(nq)                            ! constant vmr set in traceur.def
63      integer :: ilon, ilat, iq, ilay, iline, nlines, ierr
64
65      CHARACTER(len=100) :: qxf(nq)                ! vmr profile files set in traceur.def
66      CHARACTER(len=100) :: fil                    ! path files
67      character(len=500) :: tracline               ! to store a line of text
68
69      logical :: foundback = .false.
70
71! 1. initialization
72
73      pq(:,:,:,:)  = 0.
74      qsurf(:,:)   = 0.
75      pqx(:,:,:,:) = 0.
76      qx(:)        = 0.
77      qxf(:)       = 'None'
78      nlines       = 0
79
80! 2. load in traceur.def chemistry data for initialization:
81
82      ! Skip nq
83      open(90,file='traceur.def',status='old',form='formatted',iostat=ierr)
84      if (ierr.eq.0) then
85        READ(90,'(A)') tracline
86        IF (trim(tracline).eq.'#ModernTrac-v1') THEN ! Test modern traceur.def
87           DO
88             READ(90,'(A)',iostat=ierr) tracline
89             IF (ierr.eq.0) THEN
90               IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
91                 EXIT
92               ENDIF
93             ELSE ! If pb, or if reached EOF without having found number of tracer
94                write(*,*) "calchim: error reading line of tracers"
95                write(*,*) "   (first lines of traceur.def) "
96                stop
97             ENDIF
98           ENDDO
99        ENDIF ! if modern or standard traceur.def
100      else
101         write(*,*) "calchim: error opening traceur.def in inichim_newstart"
102         stop
103      endif
104
105      ! Get data of tracers
106      do iq=1,nq
107         read(90,'(A)') tracline
108         write(*,*)"inichim_newstart: iq=",iq,"noms(iq)=",trim(noms(iq))
109         if (index(tracline,'qx='    ) /= 0) then
110            read(tracline(index(tracline,'qx=')+len('qx='):),*) qx(iq)
111            write(*,*) ' Parameter value (traceur.def) : qx=', qx(iq)
112         else
113            write(*,*) ' Parameter value (default)     : qx=', qx(iq)
114         end if
115         if (index(tracline,'qxf='    ) /= 0) then
116            read(tracline(index(tracline,'qxf=')+len('qxf='):),*) qxf(iq)
117            write(*,*) ' Parameter value (traceur.def) : qxf=', qxf(iq)
118         else
119            write(*,*) ' Parameter value (default)     : qxf=', qxf(iq)
120         end if
121      end do
122
123      close(90)
124
125! 3. initialization of tracers
126
127! 3.1 vertical interpolation
128
129      do iq=1,nq
130         if (qx(iq) /= 0.) then
131            pqx(:,:,:,iq) = qx(iq)
132         else if (qxf(iq) /= 'None') then
133            ! Opening file
134            fil = trim(datadir)//'/chemical_profiles/'//qxf(iq)
135            print*, 'chemical pofile '//trim(noms(iq))//': ', fil
136            open(UNIT=90,FILE=fil,STATUS='old',iostat=ierr)
137            if (ierr.eq.0) then
138               read(90,*) ! read one header line
139               do         ! get number of lines
140                   read(90,*,iostat=ierr)
141                   if (ierr<0) exit
142                   nlines = nlines + 1
143               end do
144               ! allocate reading variable
145               allocate(pf(nlines))
146               allocate(qf(nlines))
147               ! read file
148               rewind(90) ! restart from the beggining of the file
149               read(90,*) ! read one header line
150               do iline=1,nlines
151                  read(90,*) pf(iline), qf(iline) ! pf [Pa], qf [vmr]
152               end do
153               ! interp in gcm grid
154               do ilon = 1,nbp_lon+1
155                  do ilat = 1,nbp_lat
156                     do ilay=1,nbp_lev
157                        pgcm = aps(ilay) + bps(ilay)*ps(ilon,ilat)  ! gcm pressure
158                        call intrplf(log(pgcm),pqx(ilon,ilat,ilay,iq),log(pf),qf,nlines)
159                     end do
160                  end do
161               end do
162               ! deallocate for next tracer
163               deallocate(pf)
164               deallocate(qf)
165            else
166               write(*,*) 'inichim_newstart: error opening ', fil
167               stop
168            endif
169            close(90)
170         end if
171      end do
172
173! 3.2 background gas
174
175      do iq=1,nq
176         if (qx(iq)==1.) then
177            pqx(:,:,:,iq) = 0.
178            do ilon = 1,nbp_lon+1
179               do ilat = 1,nbp_lat
180                  do ilay=1,nbp_lev
181                     pqx(ilon,ilat,ilay,iq) = 1-sum(pqx(ilon,ilat,ilay,:))
182                     if (pqx(ilon,ilat,ilay,iq)<=0.) then
183                        write(*,*) 'inichim_newstart: vmr tot > 1 not possible'
184                        stop
185                     end if
186                  end do
187               end do
188            end do
189            foundback = .true.
190            exit ! you found the background gas you can skip others
191         end if
192      end do
193      if (.not.foundback) then
194         write(*,*) 'inichim_newstart: you need to set a background gas'
195         write(*,*) '            by qx=1. in traceur.def'
196         stop
197      end if
198
199! 3.3 convert vmr to mmr
200
201      do ilon = 1,nbp_lon+1
202         do ilat = 1,nbp_lat
203            mmean(:) = 0.
204            do ilay=1,nbp_lev
205               do iq=1,nq
206                  mmean(ilay) = mmean(ilay) + pqx(ilon,ilat,ilay,iq)*mmol(iq)
207               end do
208               do iq=1,nq
209                  pq(ilon,ilat,ilay,iq) = pqx(ilon,ilat,ilay,iq)*mmol(iq)/mmean(ilay)
210               end do
211            end do
212         end do
213      end do
214
215! 4. Hard coding
216! Do whatever you want here to specify pq and qsurf
217! Or use #ModernTrac-v1 and add another option section 2.
218
219      end
Note: See TracBrowser for help on using the repository browser.