Context Navigation

mzcf.F @ 473

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

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