Context Navigation

mzcud.F @ 469

Last change on this file since 469 was 414, checked in by aslmd, 14 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

Rev	Line
[414]	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