Context Navigation

mzcud.F @ 482

Last change on this file since 482 was 414, checked in by aslmd, 13 years ago

LMDZ.MARS : NEW NLTE MODEL FROM GRANADA AMIGOS

23/11/11 == FGG + MALV

New parameterization of the NLTE 15 micron cooling. The old parameterization is kept as an option, including or not variable atomic oxygen concentration. A new flag is introduced in callphys.def, nltemodel, to select which parameterization wants to be used (new one, old one with variable [O], or old one with fixed [O], see below). Includes many new subroutines and commons in phymars. Some existing routines are also modified:

-physiq.F. Call to the new subroutine NLTE_leedat in first call. Call to nltecool modified to allow for variable atomic oxygen. Depending on the value of nltemodel, the new subroutine NLTEdlvr09_TCOOL is called instead of nltecool.

-inifis.F. Reading of nltemodel is added.

-callkeys.h Declaration of nltemodel is added.

The following lines should be added to callphys.def (ideally after setting callnlte):

# NLTE 15um scheme to use.
# 0-> Old scheme, static oxygen
# 1-> Old scheme, dynamic oxygen
# 2-> New scheme
nltemodel = 2

A new directory, NLTEDAT, has to be included in datagcm.

Improvements into NLTE NIR heating parameterization to take into account variability with O/CO2 ratio and SZA. A new subroutine, NIR_leedat.F, and a new common, NIRdata.h, are included. A new flag, nircorr, is added in callphys.def, to include or not these corrections. The following files are modified:

-nirco2abs.F: nq and pq are added in the arguments. The corrections factors are interpolated to the GCM grid and included in the clculation. A new subroutine, interpnir, is added at the end of the file.

-physiq.F: Call to NIR_leedat added at first call. Modified call to nirco2abs

-inifis: Reading new flag nircorr.

-callkeys.h: Declaration of nircorr.

The following lines have to be added to callphys.def (ideally after callnirco2):

# NIR NLTE correction for variable SZA and O/CO2?
# matters only if callnirco2=T
# 0-> no correction
# 1-> include correction
nircorr=1

A new data file, NIRcorrection_feb2011.dat, has to be included in datagcm.

Small changes to the molecular diffusion scheme to fix the number of species considered, to avoid problems when compiling with more than 15 tracers (for example, when CH4 is included). Modified subroutines: aeronomars/moldiff.F and aeronomars/moldiffcoeff.F

File size: 17.1 KB

Line
1	c***********************************************************************
2
3	subroutine mzcud( tauinf,tau, c,cup,cdw,vc,taugr,
4	@ ib,isot,icfout,itableout )
5
6	c old times mlp first version of mzcf
7	c a.k.murphy method to avoid extrapolation in the curtis matrix
8	c feb-89 malv AKM method to avoid extrapolation in C.M.
9	c 25-sept-96 cristina dejar las matrices en doble precision
10	c jan 98 malv version para mz1d
11	c oct 01 malv update version for fluxes for hb and fb
12	c jul 2011 malv+fgg Adapted to LMD-MGCM
13	c***********************************************************************
14
15	implicit none
16
17	include 'comcstfi.h'
18	include 'nltedefs.h'
19	include 'nlte_atm.h'
20	include 'nlte_data.h'
21	include 'tcr_15um.h'
22	include 'nlte_curtis.h'
23
24	c arguments
25	real*8 c(nl,nl), cup(nl,nl), cdw(nl,nl) ! o
26	real*8 vc(nl), taugr(nl) ! o
27	real*8 tau(nl,nl) ! i
28	real*8 tauinf(nl) ! i
29	integer ib ! i
30	integer isot ! i
31	integer icfout, itableout ! i
32
33	c external
34	external bandid
35	character*2 bandid
36
37	c local variables
38	integer i, in, ir, iw, itblout
39	real*8 cfup(nl,nl), cfdw(nl,nl)
40	real*8 a(nl,nl), cf(nl,nl)
41	character isotcode2, bcode2
42
43	c formats
44	101 format(i1)
45	202 format(i2)
46	180 format(a80)
47	181 format(a80)
48	c***********************************************************************
49
50	if (isot.eq.1) isotcode = '26'
51	if (isot.eq.2) isotcode = '28'
52	if (isot.eq.3) isotcode = '36'
53	if (isot.eq.4) isotcode = '27'
54	if (isot.eq.5) isotcode = 'co'
55	bcode = bandid( ib )
56
57	! write (,) ' '
58
59	do in=1,nl
60
61	do ir=1,nl
62
63	cf(in,ir) = 0.0d0
64	cfup(in,ir) = 0.0d0
65	cfdw(in,ir) = 0.0d0
66	c(in,ir) = 0.0d0
67	cup(in,ir) = 0.0d0
68	cdw(in,ir) = 0.0d0
69	a(in,ir) = 0.0d0
70
71	end do
72
73	vc(in) = 0.0d0
74	taugr(in) = 0.0d0
75
76	end do
77
78
79	c the next lines are a reduced and equivalent way of calculating
80	c the c(in,ir) elements for n=2,nl1 and r=1,nl
81
82
83	c do in=2,nl1
84	c do ir=1,nl
85	c if(ir.eq.1)then
86	c c(in,ir)=tau(in-1,1)-tau(in+1,1)
87	c elseif(ir.eq.nl)then
88	c c(in,ir)=tau(in+1,nl1)-tauinf(in+1)-tau(in-1,nl1)+tauinf(in-1)
89	c else
90	c c(in,ir)=tau(in+1,ir-1)-tau(in+1,ir)-tau(in-1,ir-1)+tau(in-1,ir)
91	c end if
92	c c(in,ir)=c(in,ir)pideltanu(ib)/(2.deltaz1.0e5)
93	c end do
94	c end do
95	c go to 1000
96
97	c calculation of the matrix cfup(nl,nl)
98
99	cfup(1,1) = 1.d0 - tau(1,1)
100
101	do in=2,nl
102	do ir=1,in
103
104	if (ir.eq.1) then
105	cfup(in,ir) = tau(in,ir) - tau(in,1)
106	elseif (ir.eq.in) then
107	cfup(in,ir) = 1.d0 - tau(in,ir-1)
108	else
109	cfup(in,ir) = tau(in,ir) - tau(in,ir-1)
110	end if
111
112	end do
113	end do
114
115	! contribution to upwards fluxes from bb at bottom :
116	do in=1,nl
117	taugr(in) = tau(in,1)
118	enddo
119
120	c calculation of the matrix cfdw(nl,nl)
121
122	cfdw(nl,nl) = 1.d0 - tauinf(nl)
123
124	do in=1,nl-1
125	do ir=in,nl
126
127	if (ir.eq.in) then
128	cfdw(in,ir) = 1.d0 - tau(in,ir)
129	elseif (ir.eq.nl) then
130	cfdw(in,ir) = tau(in,ir-1) - tauinf(in)
131	else
132	cfdw(in,ir) = tau(in,ir-1) - tau(in,ir)
133	end if
134
135	end do
136	end do
137
138
139	c calculation of the matrix cf(nl,nl)
140
141	do in=1,nl
142	do ir=1,nl
143
144	if (ir.eq.1) then
145	! version con l_bb(tg) = l_bb(t(1))=j(1) (see also vc below)
146	! cf(in,ir) = tau(in,ir)
147	! version con l_bb(tg) =/= l_bb(t(1))=j(1) (see also vc below)
148	cf(in,ir) = tau(in,ir) - tau(in,1)
149	elseif (ir.eq.nl) then
150	cf(in,ir) = tauinf(in) - tau(in,ir-1)
151	else
152	cf(in,ir) = tau(in,ir) - tau(in,ir-1)
153	end if
154
155	end do
156	end do
157
158
159	c definition of the a(nl,nl) matrix
160
161	do in=2,nl-1
162	do ir=1,nl
163	if (ir.eq.in+1) a(in,ir) = -1.d0
164	if (ir.eq.in-1) a(in,ir) = +1.d0
165	a(in,ir) = a(in,ir) / ( 2.d0deltaz1.d5 )
166	end do
167	end do
168	! this is not needed anymore in the akm scheme
169	! a(1,1) = +3.d0
170	! a(1,2) = -4.d0
171	! a(1,3) = +1.d0
172	! a(nl,nl) = -3.d0
173	! a(nl,nl1) = +4.d0
174	! a(nl,nl2) = -1.d0
175
176	c calculation of the final curtis matrix ("reduced" by murphy's method)
177
178	if (isot.ne.5) then
179	do in=1,nl
180	do ir=1,nl
181	cf(in,ir) = cf(in,ir) * pi*deltanu(isot,ib)
182	cfup(in,ir) = cfup(in,ir) * pi*deltanu(isot,ib)
183	cfdw(in,ir) = cfdw(in,ir) * pi*deltanu(isot,ib)
184	end do
185	taugr(in) = taugr(in) * pi*deltanu(isot,ib)
186	end do
187	else
188	do in=1,nl
189	do ir=1,nl
190	cf(in,ir) = cf(in,ir) * pi*deltanuco
191	enddo
192	taugr(in) = taugr(in) * pi*deltanuco
193	enddo
194	endif
195
196	do in=2,nl-1
197
198	do ir=1,nl
199
200	do i=1,nl
201	! only c contains the matrix a. matrixes cup,cdw dont because
202	! these two will be used for flux calculations, not
203	! only for flux divergencies
204
205	c(in,ir) = c(in,ir) + a(in,i) * cf(i,ir)
206	! from this matrix we will extract (see below) the
207	! nl2 x nl2 "core" for the "reduced" final curtis matrix.
208
209	end do
210	cup(in,ir) = cfup(in,ir)
211	cdw(in,ir) = cfdw(in,ir)
212
213	end do
214	! version con l_bb(tg) = l_bb(t(1))=j(1) (see cf above)
215	!vc(in) = c(in,1)
216	! version con l_bb(tg) =/= l_bb(t(1))=j(1) (see cf above)
217	if (isot.ne.5) then
218	vc(in) = pideltanu(isot,ib)/( 2.d0deltaz1.d5 )
219	@ ( tau(in-1,1) - tau(in+1,1) )
220	else
221	vc(in) = pideltanuco/( 2.d0deltaz1.d5 )
222	@ ( tau(in-1,1) - tau(in+1,1) )
223	endif
224
225	end do
226
227	5 continue
228
229	! write (,) 'mztf/1/ c(2,*) =', (c(2,i), i=1,nl)
230
231	! call elimin_dibuja(c,nl,itableout)
232
233	c ventana del smoothing de c es nw=3 y de vc es 5 (puesto en lisa):
234	c subroutine elimin_mz4(c,vc,ilayer,nl,nan,iw, nw)
235
236	iw = nan
237	if (isot.eq.4) iw = 5 ! eliminates values < 1.d-19
238	if (itableout.eq.30) then
239	itblout = 0
240	else
241	itblout = itableout
242	endif
243	call elimin_mz1d (c,vc,0,iw,itblout,nw)
244
245	! upper boundary condition
246	! j'(nl) = j'(nl1) ==> j(nl) = 2j(nl1) - j(nl2) ==>
247	do in=2,nl-1
248	c(in,nl-2) = c(in,nl-2) - c(in,nl)
249	c(in,nl-1) = c(in,nl-1) + 2.d0*c(in,nl)
250	cup(in,nl-2) = cup(in,nl-2) - cup(in,nl)
251	cup(in,nl-1) = cup(in,nl-1) + 2.d0*cup(in,nl)
252	cdw(in,nl-2) = cdw(in,nl-2) - cdw(in,nl)
253	cdw(in,nl-1) = cdw(in,nl-1) + 2.d0*cdw(in,nl)
254	end do
255	! j(nl) = j(nl1) ==>
256	! do in=2,nl1
257	! c(in,nl1) = c(in,nl1) + c(in,nl)
258	! end do
259
260	! 1000 continue
261
262
263	if (icfout.eq.1) then
264
265	! if (ib.eq.1 .or. ib.eq.12 .or. ib.eq.16 .or. ib.eq.18) then
266	! codmatrx = codmatrx_fb
267	! else
268	! codmatrx = codmatrx_hot
269	! end if
270	! if (ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
271	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
272	! ibcode2 = '0'//ibcode1
273	! else
274	! write ( ibcode2, 202) ib
275	! endif
276
277	! open ( 1, access='sequential', form='unformatted', file=
278	! @ dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat')
279	! open ( 2, access='sequential', form='unformatted', file=
280	! @ dircurtis//'cflup'//isotcode//dn//ibcode2//codmatrx//'.dat')
281	! open ( 3, access='sequential', form='unformatted', file=
282	! @ dircurtis//'cfldw'//isotcode//dn//ibcode2//codmatrx//'.dat')
283	! open ( 4, access='sequential', form='unformatted', file=
284	! @ dircurtis//'cflgr'//isotcode//dn//ibcode2//codmatrx//'.dat')
285
286	! write(1) dummy
287	! write(1) ' format: (vc(n),(ch(n,r),r=2,nl-1),n=2,nl-1)'
288	! do in=2,nl-1
289	! write(1) vc(in), (c(in,ir) , ir=2,nl-1 )
290	!! write (,) in, vc(in)
291	! end do
292
293	! write(2) dummy
294	! write(2)' format: (cfup(n,r),r=1,nl), n=1,nl)'
295	! do in=1,nl
296	! write(2) ( cup(in,ir) , ir=1,nl )
297	! end do
298
299	! write(3) dummy
300	! write(3)' format: (cfdw(n,r),r=1,nl), n=1,nl)'
301	! do in=1,nl
302	! write(3) (cdw(in,ir) , ir=1,nl )
303	! end do
304
305	! write(4) dummy
306	! write(4)' format: (taugr(n), n=1,nl)'
307	! do in=1,nl
308	! write(4) (taugr(in), ir=1,nl )
309	! end do
310	! !write (,) ' Last value in file: ', taugr(nl)
311
312	! write (*,'(1x,30hcurtis matrix written out in: ,a,a,a,a)' )
313	! @ dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat',
314	! @ dircurtis//'cflup'//isotcode//dn//ibcode2//codmatrx//'.dat',
315	! @ dircurtis//'cfldw'//isotcode//dn//ibcode2//codmatrx//'.dat',
316	! @ dircurtis//'cflgr'//isotcode//dn//ibcode2//codmatrx//'.dat'
317
318	! close (1)
319	! close (2)
320	! close (3)
321	! close (4)
322
323	else
324
325	! write (,) ' no curtis matrix output file ', char(10)
326
327	end if
328
329	if (itableout.eq.30) then ! Force output of C.M. in ascii file
330
331	! if (ib.eq.1 .or. ib.eq.12 .or. ib.eq.16 .or. ib.eq.18) then
332	! codmatrx = codmatrx_fb
333	! else
334	! codmatrx = codmatrx_hot
335	! end if
336	! if (ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
337	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
338	! ibcode2 = '0'//ibcode1
339	! else
340	! write ( ibcode2, 202) ib
341	! endif
342
343	! open (10, file=
344	! & dircurtis//'table'//isotcode//dn//ibcode2//codmatrx//'.dat')
345	! write(10,*) nl, ' = number of layers '
346	! write(10,*) ' format: (vc(n),(ch(n,r),r=2,nl-1),n=2,nl-1)'
347	! do in=2,nl-1
348	! write(10,*) vc(in), (c(in,ir) , ir=2,nl-1 )
349	! enddo
350	! close (10)
351	endif
352
353	c end
354	return
355	end

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format