source: trunk/LMDZ.MARS/libf/phymars/conc_mod.F90

Last change on this file was 3185, checked in by csegonne, 11 months ago

MARS PCM:

  • conc_mod.F : Added new subroutines init_r_cp_mu and update_r_cp_mu_ak to replace initialization of rnew, cpnew, mmean and akknew as constants and their update if callthermos or photochem (same update as in concentrations.F) that were done in physic_mod.F.
  • rnew, cpnew, mmean and akknew are now protected meaning they cannot be modified anywhere in the model other than in conc_mod.F .
  • concentrations.F has been removed.
  • Knowing where rnew, cpnew, mmean come from is now more explicit. These modifications respond to the need expressed in ticket #90.
File size: 10.1 KB
RevLine 
[1047]1module conc_mod
2
3implicit none
4
[3185]5 real,save,allocatable,protected :: mmean(:,:)  ! mean molecular mass of the atmosphere
6 real,save,allocatable,protected :: Akknew(:,:) ! thermal conductivity coefficient
7 real,save,allocatable,protected :: cpnew(:,:)  ! specific heat
8 real,save,allocatable,protected :: rnew(:,:)   ! specific gas constant
[2578]9
10!$OMP THREADPRIVATE(mmean,Akknew,cpnew,rnew)
[1047]11 
12contains
13
[3185]14 subroutine ini_conc_mod(ngrid,nlayer)
[1047]15 
[3185]16 implicit none
17 integer,intent(in) :: ngrid ! number of atmospheric columns
18 integer,intent(in) :: nlayer ! number of atmospheric levels
[1047]19 
[3185]20  allocate(mmean(ngrid,nlayer))
21  allocate(Akknew(ngrid,nlayer))
22  allocate(cpnew(ngrid,nlayer))
23  allocate(rnew(ngrid,nlayer))
[1047]24   
[3185]25 end subroutine ini_conc_mod
[1047]26
[1770]27
[3185]28 subroutine end_conc_mod
[1770]29
[3185]30 implicit none
[1770]31
[3185]32  if (allocated(mmean)) deallocate(mmean)
33  if (allocated(Akknew)) deallocate(Akknew)
34  if (allocated(cpnew)) deallocate(cpnew)
35  if (allocated(rnew)) deallocate(rnew)
[1770]36
[3185]37 end subroutine end_conc_mod
[1770]38
[3185]39 subroutine init_r_cp_mu(ngrid,nlayer)
40
41 USE comcstfi_h, only: r, cpp, mugaz
42 implicit none
43! Input 
44 integer,intent(in) :: ngrid ! number of atmospheric columns
45 integer,intent(in) :: nlayer ! number of atmospheric layers
46
47! local variables
48 integer l, ig
49
50! Initialize R, Cp and mu as constant
51!        do l=1,nlayer
52!          do ig=1,ngrid
53!            rnew(ig,l)=r
54!            cpnew(ig,l)=cpp
55 !           mmean(ig,l)=mugaz
56!          enddo
57 !       enddo
58         rnew(:,:)=r
59         cpnew(:,:)=cpp
60         mmean(:,:)=mugaz
61     return           
62 end subroutine init_r_cp_mu
63
64
65 subroutine update_r_cp_mu_ak(ngrid,nlayer,nq,pplay,pt,pdt,pq,pdq,ptimestep)
66
67 use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,             &
68&                      igcm_o2, igcm_o3, igcm_h, igcm_h2,               &
69&                      igcm_oh, igcm_ho2, igcm_n2, igcm_ar,             &
70&                      igcm_h2o_vap, igcm_n, igcm_no, igcm_no2, &
71&                      igcm_n2d, igcm_co2plus, igcm_oplus,              &
72&                      igcm_o2plus, igcm_coplus, igcm_cplus,            &
73&                      igcm_nplus, igcm_noplus, igcm_n2plus,            &
74&                      igcm_hplus, igcm_hco2plus, mmol,         &
75&                      igcm_he, igcm_elec
76 implicit none
77 !=======================================================================
78! CALCULATION OF MEAN MOLECULAR MASS, Cp, Akk and R
79!
80! mmean(ngrid,nlayer)   amu
81! cpnew(ngrid,nlayer)   J/kg/K
82! rnew(ngrid,nlayer)    J/kg/K
83! akknew(ngrid,nlayer)  coefficient of thermal conduction
84!
85! version: April 2012 - Franck Lefevre
86!=======================================================================
87
88!     declarations
89 
90    include "callkeys.h"
91
92!     input/output
93
94    integer,intent(in) :: ngrid ! number of atmospheric columns
95    integer,intent(in) :: nlayer ! number of atmospheric layers
96    integer,intent(in) :: nq ! number of tracers
97    real,intent(in) :: pplay(ngrid,nlayer)
98    real,intent(in) :: pt(ngrid,nlayer)
99    real,intent(in) :: pdt(ngrid,nlayer)
100    real,intent(in) :: pq(ngrid,nlayer,nq)
101    real,intent(in) :: pdq(ngrid,nlayer,nq)
102    real,intent(in) :: ptimestep
103
104!     local variables
105
106    integer       :: i, l, ig, iq
107    integer, save :: nbq
108    integer,allocatable,save :: niq(:)
109    real          :: ni(nq), ntot
110    real          :: zq(ngrid, nlayer, nq)
111    real          :: zt(ngrid, nlayer)
112    real,allocatable,save    :: aki(:)
113    real,allocatable,save    :: cpi(:)
114
115    logical, save :: firstcall = .true.
116
117!$OMP THREADPRIVATE(nbq,niq,aki,cpi,firstcall)
118
119    if (firstcall) then
120
121      ! allocate local saved arrays:
122       allocate(aki(nq))
123       allocate(cpi(nq))
124       allocate(niq(nq))
125!        find index of chemical tracers to use
126!        initialize thermal conductivity and specific heat coefficients
127!        !? values are estimated
128
129       nbq = 0 ! to count number of tracers used in this subroutine
130
131       if (igcm_co2 /= 0) then
132          nbq = nbq + 1
133          niq(nbq) = igcm_co2
134          aki(nbq) = 3.072e-4
135          cpi(nbq) = 0.834e3
136       end if
137       if (igcm_co /= 0) then
138          nbq = nbq + 1
139          niq(nbq) = igcm_co
140          aki(nbq) = 4.87e-4
141          cpi(nbq) = 1.034e3
142       end if
143       if (igcm_o /= 0) then
144          nbq = nbq + 1
145          niq(nbq) = igcm_o
146          aki(nbq) = 7.59e-4
147          cpi(nbq) = 1.3e3
148       end if
149       if (igcm_o1d /= 0) then
150          nbq = nbq + 1
151          niq(nbq) = igcm_o1d
152          aki(nbq) = 7.59e-4  !?
153          cpi(nbq) = 1.3e3    !?
154       end if
155       if (igcm_o2 /= 0) then
156          nbq = nbq + 1
157          niq(nbq) = igcm_o2
158          aki(nbq) = 5.68e-4
159          cpi(nbq) = 0.9194e3
160       end if
161       if (igcm_o3 /= 0) then
162          nbq = nbq + 1
163          niq(nbq) = igcm_o3
164          aki(nbq) = 3.00e-4  !?
165          cpi(nbq) = 0.800e3  !?
166       end if
167       if (igcm_h /= 0) then
168          nbq = nbq + 1
169          niq(nbq) = igcm_h
170          aki(nbq) = 0.0
171          cpi(nbq) = 20.780e3
172       end if
173       if (igcm_h2 /= 0) then
174          nbq = nbq + 1
175          niq(nbq) = igcm_h2
176          aki(nbq) = 36.314e-4
177          cpi(nbq) = 14.266e3
178       end if
179       if (igcm_oh /= 0) then
180          nbq = nbq + 1
181          niq(nbq) = igcm_oh
182          aki(nbq)  = 7.00e-4 !?
183          cpi(nbq)  = 1.045e3
184       end if
185       if (igcm_ho2 /= 0) then
186          nbq = nbq + 1
187          niq(nbq) = igcm_ho2
188          aki(nbq) = 0.0
189          cpi(nbq) = 1.065e3  !?
190       end if
191       if (igcm_n2 /= 0) then
192          nbq = nbq + 1
193          niq(nbq) = igcm_n2
194          aki(nbq) = 5.6e-4
195          cpi(nbq) = 1.034e3
196       end if
197       if (igcm_ar /= 0) then
198          nbq = nbq + 1
199          niq(nbq) = igcm_ar
200          aki(nbq) = 3.4e-4
201          cpi(nbq) = 5.2e2
202         ! aki/cpi from Leslie A. Guildner,
203         ! JOURNAL OF RESEARCH of the National Bureau of Standards-
204         ! A. Physics and Chemistry
205         ! Vol. 79A, No.2, March-April 1975
206       end if
207       if (igcm_h2o_vap /= 0) then
208          nbq = nbq + 1
209          niq(nbq) = igcm_h2o_vap
210          aki(nbq) = 0.0
211          cpi(nbq) = 1.870e3
212       end if
213       if (igcm_n /= 0) then
214          nbq = nbq + 1
215          niq(nbq) = igcm_n
216          aki(nbq) = 0.0
217          cpi(nbq) = 0.0
218       endif
219       if(igcm_no /= 0) then
220          nbq = nbq + 1
221          niq(nbq) = igcm_no
222          aki(nbq) = 0.0
223          cpi(nbq) = 0.0
224       endif
225       if(igcm_no2 /= 0) then
226          nbq = nbq + 1
227          niq(nbq) = igcm_no2
228          aki(nbq) = 0.0
229          cpi(nbq) = 0.0
230       endif
231       if(igcm_n2d /= 0) then
232          nbq = nbq + 1
233          niq(nbq) = igcm_n2d
234          aki(nbq) = 0.0
235          cpi(nbq) = 0.0
236       endif
237       if(igcm_he /= 0) then
238          nbq = nbq + 1
239          niq(nbq) = igcm_he
240          aki(nbq) = 30.e-4
241          cpi(nbq) = 5.2e3
242         ! aki/cpi from Leslie A. Guildner,
243         ! JOURNAL OF RESEARCH of the National Bureau of Standards-
244         ! A. Physics and Chemistry
245         ! Vol. 79A, No.2, March-April 1975
246       endif
247       if(igcm_co2plus /= 0) then
248          nbq = nbq + 1
249          niq(nbq) = igcm_co2plus
250          aki(nbq) = 0.0
251          cpi(nbq) = 0.0
252       endif
253       if(igcm_oplus /= 0) then
254          nbq = nbq + 1
255          niq(nbq) = igcm_oplus
256          aki(nbq) = 0.0
257          cpi(nbq) = 0.0
258       endif
259       if(igcm_o2plus /= 0) then
260          nbq = nbq + 1
261          niq(nbq) = igcm_o2plus
262          aki(nbq) = 0.0
263          cpi(nbq) = 0.0
264       endif
265       if(igcm_coplus /= 0) then
266          nbq = nbq + 1
267          niq(nbq) = igcm_coplus
268          aki(nbq) = 0.0
269          cpi(nbq) = 0.0
270       endif
271       if(igcm_cplus /= 0) then
272          nbq = nbq + 1
273          niq(nbq) = igcm_cplus
274          aki(nbq) = 0.0
275          cpi(nbq) = 0.0
276       endif
277       if(igcm_nplus /= 0) then
278          nbq = nbq + 1
279          niq(nbq) = igcm_nplus
280          aki(nbq) = 0.0
281          cpi(nbq) = 0.0
282       endif
283       if(igcm_noplus /= 0) then
284          nbq = nbq + 1
285          niq(nbq) = igcm_noplus
286          aki(nbq) = 0.0
287          cpi(nbq) = 0.0
288       endif
289       if(igcm_n2plus /= 0) then
290          nbq = nbq + 1
291          niq(nbq) = igcm_n2plus
292          aki(nbq) = 0.0
293          cpi(nbq) = 0.0
294       endif
295       if(igcm_hplus /= 0) then
296          nbq = nbq + 1
297          niq(nbq) = igcm_hplus
298          aki(nbq) = 0.0
299          cpi(nbq) = 0.0
300       endif
301       if(igcm_hco2plus /= 0) then
302          nbq = nbq + 1
303          niq(nbq) = igcm_hco2plus
304          aki(nbq) = 0.0
305          cpi(nbq) = 0.0
306       endif
307       if(igcm_elec /= 0) then
308          nbq = nbq + 1
309          niq(nbq) = igcm_elec
310          aki(nbq) = 0.0
311          cpi(nbq) = 0.0
312       endif
313
314      ! tell the world about it:
315       write(*,*) "concentrations: firstcall, nbq=",nbq
316       write(*,*) "  niq(1:nbq)=",niq(1:nbq)
317       write(*,*) "  aki(1:nbq)=",aki(1:nbq)
318       write(*,*) "  cpi(1:nbq)=",cpi(1:nbq)
319
320      firstcall = .false.
321
322    end if ! if (firstcall)
323
324!     update temperature
325
326    do l = 1,nlayer
327       do ig = 1,ngrid
328          zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
329       end do
330    end do
331
332!     update tracers
333
334    do l = 1,nlayer
335       do ig = 1,ngrid
336          do i = 1,nbq
337             iq = niq(i)
338             zq(ig,l,iq) = max(1.e-30, pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep)
339         end do
340      end do
341    end do
342
343!     mmean : mean molecular mass
344!     rnew  : specific gas constant
345
346    mmean(:,:)  = 0.
347
348    do l = 1,nlayer
349       do ig = 1,ngrid
350          do i = 1,nbq
351             iq = niq(i)
352             mmean(ig,l) = mmean(ig,l) + zq(ig,l,iq)/mmol(iq)
353          end do
354          mmean(ig,l) = 1./mmean(ig,l)
355          rnew(ig,l) = 8.314/mmean(ig,l)*1.e3     ! J/kg/K           
356       end do
357    end do
358
359!     cpnew  : specific heat
360!     akknew : thermal conductivity coefficient
361
362    cpnew(:,:)  = 0.
363    akknew(:,:) = 0.
364
365    do l = 1,nlayer
366       do ig = 1,ngrid
367          ntot = pplay(ig,l)/(1.381e-23*zt(ig,l))*1.e-6  ! in #/cm3
368          do i = 1,nbq
369             iq = niq(i)
370             ni(iq) = ntot*zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
371             cpnew(ig,l) = cpnew(ig,l) + ni(iq)*cpi(i)
372             akknew(ig,l) = akknew(ig,l) + ni(iq)*aki(i)
373          end do
374          cpnew(ig,l) = cpnew(ig,l)/ntot
375          akknew(ig,l) = akknew(ig,l)/ntot
376       end do
377    end do
378
379    return
380 end subroutine update_r_cp_mu_ak
381
[1047]382end module conc_mod
Note: See TracBrowser for help on using the repository browser.