Context Navigation

lwflux.F @ 1655

Last change on this file since 1655 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

Remove all reference/use of nlayermx and dimphys.h
Made use of automatic arrays whenever arrays are needed with dimension nlayer
Remove lots of obsolete reference to dimensions.h
Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

Deleted param.h and put contents into module param_v4_h
Adapted testphys1d, newstart, etc...
Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

Optimized computations in paramfoto_compact (twice less dlog10 calculations)
Checked consistency before/after modification in debug mode
Checked performance is not impacted (same as before)

File size: 15.3 KB

Line
1	subroutine lwflux (ig0,kdlon,kflev,dp
2	. ,bsurf,btop,blev,blay,dbsublay
3	. ,tlay, tlev, dt0 ! pour sortie dans g2d uniquement
4	. ,emis
5	. , tautotal,omegtotal,gtotal
6	. ,coolrate,fluxground,fluxtop
7	. ,netrad)
8
9	c----------------------------------------------------------------------
10	c LWFLUX computes the fluxes
11	c----------------------------------------------------------------------
12
13	use dimradmars_mod, only: ndlo2, nir, ndlon, nuco2, nflev
14	use yomlw_h, only: nlaylte, xi, xi_ground, gcp
15	implicit none
16
17	#include "callkeys.h"
18	#include "comg1d.h"
19
20	c----------------------------------------------------------------------
21	c 0.1 arguments
22	c ---------
23	c inputs:
24	c -------
25	integer ig0
26	integer kdlon ! part of ngrid
27	integer kflev ! part of nlayer
28
29	real dp (ndlo2,kflev) ! layer pressure thickness (Pa)
30
31	real bsurf (ndlo2,nir) ! surface spectral planck function
32	real blev (ndlo2,nir,kflev+1) ! level spectral planck function
33	real blay (ndlo2,nir,kflev) ! layer spectral planck function
34	real btop (ndlo2,nir) ! top spectral planck function
35	real dbsublay (ndlo2,nir,2*kflev) ! layer gradient spectral planck
36	! function in sub layers
37
38	real dt0 (ndlo2) ! surface temperature discontinuity
39	real tlay (ndlo2,kflev) ! layer temperature
40	real tlev (ndlo2,kflev+1) ! level temperature
41
42	real emis (ndlo2) ! surface emissivity
43
44	real tautotal(ndlo2,kflev,nir) ! \ Total single scattering
45	real omegtotal(ndlo2,kflev,nir) ! > properties (Addition of the
46	real gtotal(ndlo2,kflev,nir) ! / NAERKIND aerosols prop.)
47
48
49	c outputs:
50	c --------
51	real coolrate(ndlo2,kflev) ! radiative cooling rate (K/s)
52	real netrad (ndlo2,kflev) ! radiative budget (W/m2)
53	real fluxground(ndlo2) ! downward flux on the ground
54	! for surface radiative budget
55	real fluxtop(ndlo2) ! upward flux on the top of atm ("OLR")
56
57
58	c----------------------------------------------------------------------
59	c 0.2 local arrays
60	c ------------
61
62	integer ja,jl,j,i,ig1d,ig,l,ndim
63	! parameter(ndim = ndlon(nuco2+1)(nflev+2)*(nflev+2))
64	real ksidb (ndlon,nuco2+1,0:nflev+1,0:nflev+1) ! net exchange rate (W/m2)
65
66	real dpsgcp (0:nflev+1,0:nflev+1) ! dp/(g.cp)
67	real temp (0:nflev+1,0:nflev+1)
68
69	real fluxdiff(ndlon,2,nflev+1) ! diffusion flux: upward(1) downward(2)
70
71	real*4 reel4
72
73	c To compute IR flux in the atmosphere (For diagnostic only !!)
74	logical computeflux
75	real coefd(kdlon,nuco2,nflev+1,nflev+1)
76	real coefu(kdlon,nuco2,0:nflev,nflev+1)
77	real flw_up(kdlon,nflev+1), flw_dn(kdlon,nflev+1) ! fluxes (W/m2)
78
79
80	ndim = ndlon(nuco2+1)(nflev+2)*(nflev+2)
81	call zerophys(ndim, ksidb)
82
83	c----------------------------------------------------------------------
84	c 1.1 exchanges (layer i <--> all layers up to i)
85	c -------------------------------------------
86
87	do i = 1,nlaylte
88	do j = i+1,nlaylte
89	do ja = 1,nuco2
90	do jl = 1,kdlon
91
92	ksidb(jl,ja,i,j) = xi(ig0+jl,ja,i,j)
93	. * (blay(jl,ja,j)-blay(jl,ja,i))
94	c ksidb reciprocity
95	c -----------------
96	ksidb(jl,ja,j,i) = -ksidb(jl,ja,i,j)
97
98	enddo
99	enddo
100	enddo
101	enddo
102
103	c----------------------------------------------------------------------
104	c 1.2 exchanges (ground <--> all layers)
105	c ----------------------------------
106
107	do i = 1,nlaylte
108	do ja = 1,nuco2
109	do jl = 1,kdlon
110
111	ksidb(jl,ja,i,0) = xi(ig0+jl,ja,0,i)
112	. * (bsurf(jl,ja)-blay(jl,ja,i))
113	c ksidb reciprocity
114	c -----------------
115	ksidb(jl,ja,0,i) = -ksidb(jl,ja,i,0)
116
117	enddo
118	enddo
119	enddo
120
121	c--------------------------------------------------------
122	c Here we add the neighbour contribution
123	c for exchanges between ground and first layer
124	c--------------------------------------------------------
125
126	do ja = 1,nuco2
127	do jl = 1,kdlon
128
129	ksidb(jl,ja,1,0) = ksidb(jl,ja,1,0)
130	. - xi_ground(ig0+jl,ja)
131	. * (blev(jl,ja,1)-blay(jl,ja,1))
132
133	cc ksidb reciprocity
134	cc -----------------
135	ksidb(jl,ja,0,1) = - ksidb(jl,ja,1,0)
136
137	enddo
138	enddo
139
140	c----------------------------------------------------------------------
141	c 1.3 exchanges (layer i <--> space)
142	c ------------------------------
143
144	do i = 1,nlaylte
145	do ja = 1,nuco2
146	do jl = 1,kdlon
147
148	ksidb(jl,ja,i,nlaylte+1) = xi(ig0+jl,ja,i,nlaylte+1)
149	. * (-blay(jl,ja,i))
150	c ksidb reciprocity
151	c -----------------
152	ksidb(jl,ja,nlaylte+1,i) = - ksidb(jl,ja,i,nlaylte+1)
153
154	enddo
155	enddo
156	enddo
157
158	c----------------------------------------------------------------------
159	c 1.4 exchanges (ground <--> space)
160	c -----------------------------
161
162	do ja = 1,nuco2
163	do jl = 1,kdlon
164
165	ksidb(jl,ja,0,nlaylte+1) = xi(ig0+jl,ja,0,nlaylte+1)
166	. * (-bsurf(jl,ja))
167
168	c ksidb reciprocity
169	c -----------------
170	ksidb(jl,ja,nlaylte+1,0) = - ksidb(jl,ja,0,nlaylte+1)
171
172	enddo
173	enddo
174
175	c----------------------------------------------------------------------
176	c 2.0 sum of band 1 and 2 of co2 contribution
177	c ---------------------------------------
178
179	do i = 0,nlaylte+1
180	do j = 0,nlaylte+1
181	do jl = 1,kdlon
182
183	ksidb(jl,3,i,j)= ksidb(jl,1,i,j) + ksidb(jl,2,i,j)
184
185	enddo
186	enddo
187	enddo
188
189	c----------------------------------------------------------------------
190	c 3.0 Diffusion
191	c ---------
192
193	i = nlaylte+1
194	do jl = 1,kdlon
195	fluxdiff(jl,1,i) = 0.
196	fluxdiff(jl,2,i) = 0.
197	enddo
198
199	call lwdiff (kdlon,kflev
200	. ,bsurf,btop,dbsublay
201	. ,tautotal,omegtotal,gtotal
202	. ,emis,fluxdiff)
203
204	c----------------------------------------------------------------------
205	c 4.0 Radiative Budget for each layer i
206	c ---------------------------------
207
208	do i = 1,nlaylte
209	do jl = 1,kdlon
210	netrad(jl,i) = 0.
211	enddo
212	enddo
213
214	do i = 1,nlaylte
215	do j = 0,nlaylte+1
216	do jl = 1,kdlon
217
218	netrad(jl,i) = netrad(jl,i) + ksidb(jl,3,i,j)
219
220	enddo
221	enddo
222	enddo
223	c diffusion contribution
224	c ----------------------
225	do i = 1,nlaylte
226	do jl = 1,kdlon
227
228	netrad(jl,i) = netrad(jl,i)
229	. - fluxdiff(jl,1,i+1) - fluxdiff(jl,2,i+1)
230	. + fluxdiff(jl,1,i) + fluxdiff(jl,2,i)
231
232	enddo
233	enddo
234
235	c----------------------------------------------------------------------
236	c 4.0 cooling rate for each layer i
237	c -----------------------------
238
239	do i = 1,nlaylte
240	do jl = 1,kdlon
241
242	coolrate(jl,i) = gcp * netrad(jl,i) / dp(jl,i)
243
244	enddo
245	enddo
246
247	c----------------------------------------------------------------------
248	c 5.0 downward flux (all layers --> ground): "fluxground"
249	c ---------------------------------------------------
250
251	do jl = 1,kdlon
252	fluxground(jl) = 0.
253	enddo
254
255	do i = 1,nlaylte
256	do ja = 1,nuco2
257	do jl = 1,kdlon
258
259	fluxground(jl) = fluxground(jl)
260	. + xi(ig0+jl,ja,0,i) * (blay(jl,ja,i))
261
262	enddo
263	enddo
264	enddo
265
266	do jl = 1,kdlon
267	fluxground(jl) = fluxground(jl) - fluxdiff(jl,2,1)
268	enddo
269
270	c----------------------------------------------------------------------
271	c 6.0 outgoing flux (all layers --> space): "fluxtop"
272	c ---------------------------------------------------
273
274	do jl = 1,kdlon
275	fluxtop(jl) = 0.
276	enddo
277
278	do i = 0,nlaylte
279	do jl = 1,kdlon
280	fluxtop(jl) = fluxtop(jl)- ksidb(jl,3,i,nlaylte+1)
281	enddo
282	enddo
283
284	do jl = 1,kdlon
285	fluxtop(jl) = fluxtop(jl) + fluxdiff(jl,1,nlaylte+1)
286	enddo
287
288	c----------------------------------------------------------------------
289	c 6.5 ONLY FOR DIAGNOSTIC : Compute IR flux in the atmosphere
290	c -------------------
291	c The broadband fluxes (W.m-2) at every level from surface level (l=1)
292	c up the top of the upper layer (here: l=nlaylte+1) are:
293	c upward : flw_up(ig1d,l) ; downward : flw_dn(ig1d,j)
294	c
295	computeflux = .false.
296
297	IF (computeflux) THEN ! not used by the GCM only for diagnostic !
298	c upward flux
299	c ~~~~~~~~~~~
300	do i = 0,nlaylte
301	do j = 1,nlaylte+1
302	do ja = 1,nuco2
303	do jl = 1,kdlon
304	coefu(jl,ja,i,j) =0.
305	do l=j,nlaylte+1
306	coefu(jl,ja,i,j)=coefu(jl,ja,i,j)+xi(ig0+jl,ja,l,i)
307	end do
308
309	enddo
310	enddo
311	enddo
312	enddo
313	do j = 1,nlaylte+1
314	do jl = 1,kdlon
315	flw_up(jl,j) = 0.
316	do ja = 1,nuco2
317	flw_up(jl,j)=flw_up(jl,j)+bsurf(jl,ja)*coefu(jl,ja,0,j)
318	do i=1,j-1
319	flw_up(jl,j)=flw_up(jl,j)+blay(jl,ja,i)*coefu(jl,ja,i,j)
320	end do
321	end do
322	flw_up(jl,j)=flw_up(jl,j) + fluxdiff(jl,1,j)
323	end do
324	end do
325
326	c downward flux
327	c ~~~~~~~~~~~~~
328	do i = 1,nlaylte+1
329	do j = 1,nlaylte+1
330	do ja = 1,nuco2
331	do jl = 1,kdlon
332	coefd(jl,ja,i,j) =0.
333	do l=0,j-1
334	coefd(jl,ja,i,j)=coefd(jl,ja,i,j)+xi(ig0+jl,ja,l,i)
335	end do
336	enddo
337	enddo
338	enddo
339	enddo
340	do j = 1,nlaylte+1
341	do jl = 1,kdlon
342	flw_dn(jl,j) = 0.
343	do ja = 1,nuco2
344	do i=j,nlaylte
345	flw_dn(jl,j)=flw_dn(jl,j)+blay(jl,ja,i)*coefd(jl,ja,i,j)
346	end do
347	end do
348	flw_dn(jl,j)=flw_dn(jl,j) - fluxdiff(jl,2,j)
349	end do
350	end do
351	END IF
352
353	c----------------------------------------------------------------------
354	c 7.0 outputs Grads 2D
355	c ----------------
356
357	c ig1d: point de la grille physique ou on veut faire la sortie
358	c ig0+1: point du decoupage de la grille physique
359
360	c#ifdef undim
361	if (callg2d) then
362
363	ig1d = kdlon/2 + 1
364	c ig1d = kdlon
365
366	if ((ig0+1).LE.ig1d .and. ig1d.LE.(ig0+kdlon)
367	. .OR. kdlon.EQ.1 ) then
368
369	ig = ig1d-ig0
370	print*, 'Sortie g2d: ig1d, ig, ig0', ig1d, ig, ig0
371
372	c--------------------------------------------
373	c Ouverture de g2d.dat
374	c--------------------------------------------
375	if (g2d_premier) then
376	open (47,file='g2d.dat'
377	clmd & ,form='unformatted',access='direct',recl=4)
378	& ,form='unformatted',access='direct',recl=1
379	& ,status='unknown')
380	g2d_irec=0
381	g2d_appel=0
382	g2d_premier=.false.
383	endif
384	g2d_appel = g2d_appel+1
385
386	c--------------------------------------------
387	c Sortie g2d des xi proches + distants
388	c--------------------------------------------
389	cl if (nflev .NE. 500) then
390	do ja = 1,nuco2
391	do j = 0,nlaylte+1
392	do i = 0,nlaylte+1
393	g2d_irec=g2d_irec+1
394	reel4 = xi(ig1d,ja,i,j)
395	write(47,rec=g2d_irec) reel4
396	enddo
397	enddo
398	enddo
399	cl endif
400
401	c------------------------------------------------------
402	c Writeg2d des ksidb
403	c------------------------------------------------------
404	do ja = 1,nuco2
405	c ja=1
406	do j = 0,nlaylte+1
407	do i = 0,nlaylte+1
408	g2d_irec=g2d_irec+1
409	reel4 = ksidb(ig,ja,i,j)
410	write(47,rec=g2d_irec) reel4
411	enddo
412	enddo
413	enddo
414
415	do j = 0,nlaylte+1
416	do i = 0,nlaylte+1
417	g2d_irec=g2d_irec+1
418	reel4 = ksidb(ig,3,i,j)
419	write(47,rec=g2d_irec) reel4
420	enddo
421	enddo
422
423	c------------------------------------------------------
424	c Writeg2d dpsgcp
425	c------------------------------------------------------
426
427	do j = 1 , nlaylte
428	do i = 0 , nlaylte+1
429	dpsgcp(i,j) = dp(ig,j) / gcp
430	enddo
431	enddo
432
433	do i = 0 , nlaylte+1
434	c dpsgcp(i,0) = 0.0002 ! (rapport ~ entre 1000 et 10000 pour le sol)
435	dpsgcp(i,0) = 1. ! (pour regler l'echelle des sorties)
436	dpsgcp(i,nlaylte+1) = 0.
437	enddo
438
439	c print*
440	c print*,'gcp: ',gcp
441	c print*
442	c do i = 0 , nlaylte+1
443	c print*,i,'dp: ',dp(ig,i)
444	c enddo
445	c print*
446	c do i = 0 , nlaylte+1
447	c print*,i,'dpsgcp: ',dpsgcp(i,1)
448	c enddo
449
450	do j = 0,nlaylte+1
451	do i = 0,nlaylte+1
452	g2d_irec=g2d_irec+1
453	reel4 = dpsgcp(i,j)
454	write(47,rec=g2d_irec) reel4
455	enddo
456	enddo
457
458	c------------------------------------------------------
459	c Writeg2d temperature
460	c------------------------------------------------------
461
462	do j = 1 , nlaylte
463	do i = 0 , nlaylte+1
464	temp(i,j) = tlay(ig,j)
465	enddo
466	enddo
467
468	do i = 0 , nlaylte+1
469	temp(i,0) = tlev(ig,1)+dt0(ig) ! temperature surface
470	temp(i,nlaylte+1) = 0. ! temperature espace (=0)
471	enddo
472
473	do j = 0,nlaylte+1
474	do i = 0,nlaylte+1
475	g2d_irec=g2d_irec+1
476	reel4 = temp(i,j)
477	write(47,rec=g2d_irec) reel4
478	enddo
479	enddo
480
481	write(76,*) 'ig1d, ig, ig0', ig1d, ig, ig0
482	write(76,*) 'nlaylte', nlaylte
483	write(76,*) 'nflev', nflev
484	write(76,*) 'kdlon', kdlon
485	write(76,*) 'ndlo2', ndlo2
486	write(76,*) 'ndlon', ndlon
487	do ja=1,4
488	write(76,*) 'bsurf', ja, bsurf(ig,ja)
489	write(76,*) 'btop', ja, btop(ig,ja)
490
491	do j=1,nlaylte+1
492	write(76,*) 'blev', ja, j, blev(ig,ja,j)
493	enddo
494
495	do j=1,nlaylte
496	write(76,*) 'blay', ja, j, blay(ig,ja,j)
497	enddo
498
499	do j=1,2*nlaylte
500	write(76,*) 'dbsublay', ja, j, dbsublay(ig,ja,j)
501	enddo
502	enddo
503
504	endif
505	c************************************************************************
506	c#endif
507	endif ! callg2d
508
509	return
510	end

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

Download in other formats:

Original Format