Context Navigation

mzcf.F @ 414

Last change on this file since 414 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: 16.4 KB

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

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

Download in other formats:

Original Format