Context Navigation

lwflux.F @ 3027

Last change on this file since 3027 was 3004, checked in by emillour, 2 years ago
Mars PCM: Some code cleanup. Remove obsolete "comg1d.h" and "writeg1d.F" (were used to specifically output for Grads in 1D). Also turned lwi and lwflux into modules while at it. EM
File size: 11.3 KB

Rev	Line
[3004]	1	module lwflux_mod
	2
	3	implicit none
	4
	5	contains
	6
	7	subroutine lwflux (ig0,kdlon,kflev,dp
[38]	8	. ,bsurf,btop,blev,blay,dbsublay
	9	. ,tlay, tlev, dt0 ! pour sortie dans g2d uniquement
	10	. ,emis
	11	. , tautotal,omegtotal,gtotal
	12	. ,coolrate,fluxground,fluxtop
	13	. ,netrad)
	14
	15	c----------------------------------------------------------------------
	16	c LWFLUX computes the fluxes
	17	c----------------------------------------------------------------------
	18
[1047]	19	use dimradmars_mod, only: ndlo2, nir, ndlon, nuco2, nflev
	20	use yomlw_h, only: nlaylte, xi, xi_ground, gcp
[38]	21	implicit none
	22
[1917]	23	include "callkeys.h"
[38]	24
	25	c----------------------------------------------------------------------
	26	c 0.1 arguments
	27	c ---------
	28	c inputs:
	29	c -------
[1917]	30	integer,intent(in) :: ig0
	31	integer,intent(in) :: kdlon ! part of ngrid
	32	integer,intent(in) :: kflev ! part of nlayer
[38]	33
[1917]	34	real,intent(in) :: dp (ndlo2,kflev) ! layer pressure thickness (Pa)
[38]	35
[1917]	36	real,intent(in) :: bsurf (ndlo2,nir) ! surface spectral planck function
	37	real,intent(in) :: blev (ndlo2,nir,kflev+1) ! level spectral planck function
	38	real,intent(in) :: blay (ndlo2,nir,kflev) ! layer spectral planck function
	39	real,intent(in) :: btop (ndlo2,nir) ! top spectral planck function
	40	real,intent(in) :: dbsublay (ndlo2,nir,2*kflev) ! layer gradient spectral planck
	41	! function in sub layers
[38]	42
[1917]	43	real,intent(in) :: dt0 (ndlo2) ! surface temperature discontinuity
	44	real,intent(in) :: tlay (ndlo2,kflev) ! layer temperature
	45	real,intent(in) :: tlev (ndlo2,kflev+1) ! level temperature
[38]	46
[1917]	47	real,intent(in) :: emis (ndlo2) ! surface emissivity
[38]	48
[1917]	49	real,intent(in) :: tautotal(ndlo2,kflev,nir) ! \ Total single scattering
	50	real,intent(in) :: omegtotal(ndlo2,kflev,nir) ! > properties (Addition of the
	51	real,intent(in) :: gtotal(ndlo2,kflev,nir) ! / NAERKIND aerosols prop.)
[38]	52
	53
	54	c outputs:
	55	c --------
[1917]	56	real,intent(out) :: coolrate(ndlo2,kflev) ! radiative cooling rate (K/s)
	57	real,intent(out) :: netrad (ndlo2,kflev) ! radiative budget (W/m2)
	58	real,intent(out) :: fluxground(ndlo2) ! downward flux on the ground
	59	! for surface radiative budget
	60	real,intent(out) :: fluxtop(ndlo2) ! upward flux on the top of atm ("OLR")
[38]	61
	62
	63	c----------------------------------------------------------------------
	64	c 0.2 local arrays
	65	c ------------
	66
[1917]	67	integer ja,jl,j,i,ig1d,ig,l
[38]	68	real ksidb (ndlon,nuco2+1,0:nflev+1,0:nflev+1) ! net exchange rate (W/m2)
	69
	70	real dpsgcp (0:nflev+1,0:nflev+1) ! dp/(g.cp)
	71	real temp (0:nflev+1,0:nflev+1)
	72
	73	real fluxdiff(ndlon,2,nflev+1) ! diffusion flux: upward(1) downward(2)
	74
	75	real*4 reel4
	76
	77	c To compute IR flux in the atmosphere (For diagnostic only !!)
	78	logical computeflux
[1047]	79	real coefd(kdlon,nuco2,nflev+1,nflev+1)
	80	real coefu(kdlon,nuco2,0:nflev,nflev+1)
	81	real flw_up(kdlon,nflev+1), flw_dn(kdlon,nflev+1) ! fluxes (W/m2)
[38]	82
	83
[1917]	84	ksidb(:,:,:,:)=0
[38]	85
	86	c----------------------------------------------------------------------
	87	c 1.1 exchanges (layer i <--> all layers up to i)
	88	c -------------------------------------------
	89
	90	do i = 1,nlaylte
	91	do j = i+1,nlaylte
	92	do ja = 1,nuco2
	93	do jl = 1,kdlon
	94
	95	ksidb(jl,ja,i,j) = xi(ig0+jl,ja,i,j)
	96	. * (blay(jl,ja,j)-blay(jl,ja,i))
	97	c ksidb reciprocity
	98	c -----------------
	99	ksidb(jl,ja,j,i) = -ksidb(jl,ja,i,j)
	100
	101	enddo
	102	enddo
	103	enddo
	104	enddo
	105
	106	c----------------------------------------------------------------------
	107	c 1.2 exchanges (ground <--> all layers)
	108	c ----------------------------------
	109
	110	do i = 1,nlaylte
	111	do ja = 1,nuco2
	112	do jl = 1,kdlon
	113
	114	ksidb(jl,ja,i,0) = xi(ig0+jl,ja,0,i)
	115	. * (bsurf(jl,ja)-blay(jl,ja,i))
	116	c ksidb reciprocity
	117	c -----------------
	118	ksidb(jl,ja,0,i) = -ksidb(jl,ja,i,0)
	119
	120	enddo
	121	enddo
	122	enddo
	123
	124	c--------------------------------------------------------
	125	c Here we add the neighbour contribution
	126	c for exchanges between ground and first layer
	127	c--------------------------------------------------------
	128
	129	do ja = 1,nuco2
	130	do jl = 1,kdlon
	131
	132	ksidb(jl,ja,1,0) = ksidb(jl,ja,1,0)
	133	. - xi_ground(ig0+jl,ja)
	134	. * (blev(jl,ja,1)-blay(jl,ja,1))
	135
	136	cc ksidb reciprocity
	137	cc -----------------
	138	ksidb(jl,ja,0,1) = - ksidb(jl,ja,1,0)
	139
	140	enddo
	141	enddo
	142
	143	c----------------------------------------------------------------------
	144	c 1.3 exchanges (layer i <--> space)
	145	c ------------------------------
	146
	147	do i = 1,nlaylte
	148	do ja = 1,nuco2
	149	do jl = 1,kdlon
	150
	151	ksidb(jl,ja,i,nlaylte+1) = xi(ig0+jl,ja,i,nlaylte+1)
	152	. * (-blay(jl,ja,i))
	153	c ksidb reciprocity
	154	c -----------------
	155	ksidb(jl,ja,nlaylte+1,i) = - ksidb(jl,ja,i,nlaylte+1)
	156
	157	enddo
	158	enddo
	159	enddo
	160
	161	c----------------------------------------------------------------------
	162	c 1.4 exchanges (ground <--> space)
	163	c -----------------------------
	164
	165	do ja = 1,nuco2
	166	do jl = 1,kdlon
	167
	168	ksidb(jl,ja,0,nlaylte+1) = xi(ig0+jl,ja,0,nlaylte+1)
	169	. * (-bsurf(jl,ja))
	170
	171	c ksidb reciprocity
	172	c -----------------
	173	ksidb(jl,ja,nlaylte+1,0) = - ksidb(jl,ja,0,nlaylte+1)
	174
	175	enddo
	176	enddo
	177
	178	c----------------------------------------------------------------------
	179	c 2.0 sum of band 1 and 2 of co2 contribution
	180	c ---------------------------------------
	181
	182	do i = 0,nlaylte+1
	183	do j = 0,nlaylte+1
	184	do jl = 1,kdlon
	185
	186	ksidb(jl,3,i,j)= ksidb(jl,1,i,j) + ksidb(jl,2,i,j)
	187
	188	enddo
	189	enddo
	190	enddo
	191
	192	c----------------------------------------------------------------------
	193	c 3.0 Diffusion
	194	c ---------
	195
	196	i = nlaylte+1
	197	do jl = 1,kdlon
	198	fluxdiff(jl,1,i) = 0.
	199	fluxdiff(jl,2,i) = 0.
	200	enddo
	201
	202	call lwdiff (kdlon,kflev
	203	. ,bsurf,btop,dbsublay
	204	. ,tautotal,omegtotal,gtotal
	205	. ,emis,fluxdiff)
	206
	207	c----------------------------------------------------------------------
	208	c 4.0 Radiative Budget for each layer i
	209	c ---------------------------------
	210
	211	do i = 1,nlaylte
	212	do jl = 1,kdlon
	213	netrad(jl,i) = 0.
	214	enddo
	215	enddo
	216
	217	do i = 1,nlaylte
	218	do j = 0,nlaylte+1
	219	do jl = 1,kdlon
	220
	221	netrad(jl,i) = netrad(jl,i) + ksidb(jl,3,i,j)
	222
	223	enddo
	224	enddo
	225	enddo
	226	c diffusion contribution
	227	c ----------------------
	228	do i = 1,nlaylte
	229	do jl = 1,kdlon
	230
	231	netrad(jl,i) = netrad(jl,i)
	232	. - fluxdiff(jl,1,i+1) - fluxdiff(jl,2,i+1)
	233	. + fluxdiff(jl,1,i) + fluxdiff(jl,2,i)
	234
	235	enddo
	236	enddo
	237
	238	c----------------------------------------------------------------------
	239	c 4.0 cooling rate for each layer i
	240	c -----------------------------
	241
	242	do i = 1,nlaylte
	243	do jl = 1,kdlon
	244
	245	coolrate(jl,i) = gcp * netrad(jl,i) / dp(jl,i)
	246
	247	enddo
	248	enddo
	249
	250	c----------------------------------------------------------------------
	251	c 5.0 downward flux (all layers --> ground): "fluxground"
	252	c ---------------------------------------------------
	253
	254	do jl = 1,kdlon
	255	fluxground(jl) = 0.
	256	enddo
	257
	258	do i = 1,nlaylte
	259	do ja = 1,nuco2
	260	do jl = 1,kdlon
	261
	262	fluxground(jl) = fluxground(jl)
	263	. + xi(ig0+jl,ja,0,i) * (blay(jl,ja,i))
	264
	265	enddo
	266	enddo
	267	enddo
	268
	269	do jl = 1,kdlon
	270	fluxground(jl) = fluxground(jl) - fluxdiff(jl,2,1)
	271	enddo
	272
	273	c----------------------------------------------------------------------
	274	c 6.0 outgoing flux (all layers --> space): "fluxtop"
	275	c ---------------------------------------------------
	276
	277	do jl = 1,kdlon
	278	fluxtop(jl) = 0.
	279	enddo
	280
	281	do i = 0,nlaylte
	282	do jl = 1,kdlon
	283	fluxtop(jl) = fluxtop(jl)- ksidb(jl,3,i,nlaylte+1)
	284	enddo
	285	enddo
	286
	287	do jl = 1,kdlon
	288	fluxtop(jl) = fluxtop(jl) + fluxdiff(jl,1,nlaylte+1)
	289	enddo
	290
	291	c----------------------------------------------------------------------
	292	c 6.5 ONLY FOR DIAGNOSTIC : Compute IR flux in the atmosphere
	293	c -------------------
	294	c The broadband fluxes (W.m-2) at every level from surface level (l=1)
	295	c up the top of the upper layer (here: l=nlaylte+1) are:
	296	c upward : flw_up(ig1d,l) ; downward : flw_dn(ig1d,j)
	297	c
	298	computeflux = .false.
	299
	300	IF (computeflux) THEN ! not used by the GCM only for diagnostic !
	301	c upward flux
	302	c ~~~~~~~~~~~
	303	do i = 0,nlaylte
	304	do j = 1,nlaylte+1
	305	do ja = 1,nuco2
	306	do jl = 1,kdlon
	307	coefu(jl,ja,i,j) =0.
	308	do l=j,nlaylte+1
	309	coefu(jl,ja,i,j)=coefu(jl,ja,i,j)+xi(ig0+jl,ja,l,i)
	310	end do
	311
	312	enddo
	313	enddo
	314	enddo
	315	enddo
	316	do j = 1,nlaylte+1
	317	do jl = 1,kdlon
	318	flw_up(jl,j) = 0.
	319	do ja = 1,nuco2
	320	flw_up(jl,j)=flw_up(jl,j)+bsurf(jl,ja)*coefu(jl,ja,0,j)
	321	do i=1,j-1
	322	flw_up(jl,j)=flw_up(jl,j)+blay(jl,ja,i)*coefu(jl,ja,i,j)
	323	end do
	324	end do
	325	flw_up(jl,j)=flw_up(jl,j) + fluxdiff(jl,1,j)
	326	end do
	327	end do
	328
	329	c downward flux
	330	c ~~~~~~~~~~~~~
	331	do i = 1,nlaylte+1
	332	do j = 1,nlaylte+1
	333	do ja = 1,nuco2
	334	do jl = 1,kdlon
	335	coefd(jl,ja,i,j) =0.
	336	do l=0,j-1
	337	coefd(jl,ja,i,j)=coefd(jl,ja,i,j)+xi(ig0+jl,ja,l,i)
	338	end do
	339	enddo
	340	enddo
	341	enddo
	342	enddo
	343	do j = 1,nlaylte+1
	344	do jl = 1,kdlon
	345	flw_dn(jl,j) = 0.
	346	do ja = 1,nuco2
	347	do i=j,nlaylte
	348	flw_dn(jl,j)=flw_dn(jl,j)+blay(jl,ja,i)*coefd(jl,ja,i,j)
	349	end do
	350	end do
	351	flw_dn(jl,j)=flw_dn(jl,j) - fluxdiff(jl,2,j)
	352	end do
	353	end do
[3004]	354	END IF ! of IF (computeflux)
[38]	355
[3004]	356	end subroutine lwflux
	357
	358	end module lwflux_mod

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

Download in other formats:

Original Format