Context Navigation

callradite.F @ 1292

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

Line
1	SUBROUTINE callradite(icount,ngrid,nlayer,nq,zday,ls,pq,albedo,
2	$ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
3	$ dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
4	& tauref,tau,aerosol,tauscaling,taucloudtes,rdust,rice,
5	& nuice,co2ice)
6
7	use dimradmars_mod, only: ndomainsz, nflev, nsun, nir
8	use dimradmars_mod, only: naerkind, name_iaer,
9	& iaer_dust_conrath,iaer_dust_doubleq,
10	& iaer_dust_submicron,iaer_h2o_ice
11	use yomlw_h, only: gcp, nlaylte
12	USE comcstfi_h
13	IMPLICIT NONE
14	c=======================================================================
15	c subject:
16	c --------
17	c Subroutine designed to call the main canonic
18	c radiative transfer subroutine "lwmain" et "swmain"
19	c to compute radiative heating and cooling rate and
20	c radiative fluxes to the surface.
21	c
22	c These calculations are only valid on the part of the atmosphere
23	c where Local Thermal Equilibrium (NLTE) is verified. In practice
24	c The calculations are only performed for the first "nlaylte"
25	c parameters (nlaylte is calculated by subroutine "nlthermeq"
26	c and stored in module "yomlw_h").
27	c
28	c The purpose of this subroutine is to:
29	c 1) Make some initial calculation at first call
30	c 2) Split the calculation in several sub-grid
31	c ("sub-domain") to save memory and
32	c be able run on a workstation at high resolution
33	c The sub-grid size is defined in dimradmars_mod
34	c 3) Compute the 3D scattering parameters depending on the
35	c size distribution of the different tracers (added by JBM)
36	c 4) call "lwmain" and "swmain"
37	c
38	c
39	c authors:
40	c ------
41	c Francois Forget / Christophe Hourdin / J.-B. Madeleine (2009)
42	c
43	c
44	c 3D scattering scheme user's guide (J.-B. Madeleine)
45	c ---------------------------------
46	c
47	c This routine has been modified to take into account 3D, time
48	c dependent scattering properties of the aerosols.
49	c---- The look-up tables that contain the scattering parameters
50	c of a given tracer, for different sizes, are read by SUAER.F90.
51	c The names of the corresponding ASCII files have to be set in
52	c this subroutine (file_id variable), and files must be in the
53	c directory specified in datafile.h. Please make sure that the
54	c ASCII files are correctly written, and that the range
55	c of particle sizes is consistent with what you would expect.
56	c---- SUAER.F90 is in charge of reading the ASCII files and averaging
57	c the scattering parameters in each GCM channel, using the three last
58	c equations of Forget et al. 1998 (GRL 25, No.7, p.1105-1108).
59	c---- These look-up tables, loaded during the firstcall, are then
60	c constantly used by the subroutine "aeroptproperties.F" to compute,
61	c online, the 3D scattering parameters, based on the size distribution
62	c (reffrad and nueffrad) of the different tracers, in each grid box.
63	c These 3D size distributions are loaded by the "updatereffrad.F"
64	c subroutine. A log-normal distribution is then assumed in
65	c "aeroptproperties.F", along with a Gauss-Legendre integration.
66	c---- The optical depth at the visible reference wavelength (set in
67	c SUAER.F90, after the file_id variable) is then computed by
68	c the subroutine "aeropacity.F", by using the size and spatial
69	c distribution of the corresponding tracer. This connection has to
70	c be implemented in "aeropacity.F" when adding a new tracer. To do so,
71	c one can use equation 2 of Forget et al. 1998 (Icarus 131, p.302-316).
72	c---- The resulting variables "aerosol", "QVISsQREF3d", "omegaVIS3d" and
73	c "gVIS3d" (same in the infrared) are finally used by lwmain.F and
74	c swmain.F to solve the radiative transfer equation.
75	c
76	c changes:
77	c -------
78	c
79	c > SRL 7/2000
80	c
81	c This version has been modified to only calculate radiative tendencies
82	c over layers 1..NFLEV (set in dimradmars_mod). Returns zero for higher
83	c layers, if any.
84	c In other routines, nlayer -> nflev.
85	c Routines affected: lwflux, lwi, lwmain, lwxb, lwxd, lwxn.
86	c
87	c > J.-B. Madeleine 10W12
88	c This version uses the variable's splitting, which can be usefull
89	c when performing very high resolution simulation like LES.
90	c
91	c ----------
92	c Here, solar band#1 is spectral interval between "long1vis" and "long2vis"
93	c set in dimradmars_mod
94	c Here, solar band#2 is spectral interval between "long2vis" and "long3vis"
95	c set in dimradmars_mod
96	c
97	c input:
98	c -----
99	c icount counter of call to subroutine physic by gcm
100	c ngrid number of gridpoint of horizontal grid
101	c nlayer Number of layer
102	c nq Number of tracer
103	c ls Solar longitude (Ls) , radian
104	c zday Date (time since Ls=0, in martian days)
105	c pq(ngrid,nlayer,nq) Advected fields
106	c
107	c albedo (ngrid,2) hemispheric surface albedo
108	c albedo (i,1) : mean albedo for solar band#1
109	c (see below)
110	c albedo (i,2) : mean albedo for solar band#2
111	c (see below)
112	c emis Thermal IR surface emissivity (no unit)
113	c mu0(ngrid) cos of solar zenith angle
114	c (=1 when sun at zenith)
115	c pplay(ngrid,nlayer) pressure (Pa) in the middle of each layer
116	c pplev(ngrid,nlayer+1) pressure (Pa) at boundaries of each layer
117	c pt(ngrid,nlayer) atmospheric temperature in each layer (K)
118	c tsurf(ngrid) surface temperature (K)
119	c fract(ngrid) day fraction of the time interval
120	c =1 during the full day ; =0 during the night
121	c declin latitude of subsolar point
122	c dist_sol sun-Mars distance (AU)
123	c igout coordinate of analysed point for debugging
124	c reffrad(ngrid,nlayer,naerkind) Aerosol effective radius
125	c nueffrad(ngrid,nlayer,naerkind) Aerosol effective variance
126
127	c
128	c output:
129	c -------
130	c dtlw (ngrid,nlayer) longwave (IR) heating rate (K/s)
131	c dtsw(ngrid,nlayer) shortwave (Solar) heating rate (K/s)
132	c fluxsurf_lw(ngrid) surface downward flux tota LW (thermal IR) (W.m-2)
133	c fluxsurf_sw(ngrid,1) surface downward flux SW for solar band#1 (W.m-2)
134	c fluxsurf_sw(ngrid,2) surface downward flux SW for solar band#2 (W.m-2)
135	c
136	c fluxtop_lw(ngrid) outgoing upward flux tota LW (thermal IR) (W.m-2)
137	c fluxtop_sw(ngrid,1) outgoing upward flux SW for solar band#1 (W.m-2)
138	c fluxtop_sw(ngrid,2) outgoing upward flux SW for solar band#2 (W.m-2)
139
140	c tauref Prescribed mean column optical depth at 610 Pa
141	c tau Column total visible dust optical depth at each point
142	c aerosol(ngrid,nlayer,naerkind) aerosol extinction optical depth
143	c at reference wavelength "longrefvis" set
144	c in dimradmars_h , in each layer, for one of
145	c the "naerkind" kind of aerosol optical
146	c properties.
147
148	c=======================================================================
149	c
150	c Declarations :
151	c -------------
152	c
153	#include "callkeys.h"
154
155	c-----------------------------------------------------------------------
156	c Input/Output
157	c ------------
158	INTEGER,INTENT(IN) :: icount
159	INTEGER,INTENT(IN) :: ngrid,nlayer,nq
160	INTEGER,INTENT(IN) :: igout
161
162	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
163	REAL,INTENT(IN) :: tauscaling(ngrid) ! Conversion factor for
164	! qdust and Ndust
165	REAL,INTENT(IN) :: albedo(ngrid,2),emis(ngrid)
166	REAL,INTENT(IN) :: ls,zday
167
168	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
169	REAL,INTENT(IN) :: pt(ngrid,nlayer)
170	REAL,INTENT(IN) :: tsurf(ngrid)
171	REAL,INTENT(IN) :: dist_sol,mu0(ngrid),fract(ngrid)
172	REAL,INTENT(OUT) :: dtlw(ngrid,nlayer),dtsw(ngrid,nlayer)
173	REAL,INTENT(OUT) :: fluxsurf_lw(ngrid), fluxtop_lw(ngrid)
174	REAL,INTENT(OUT) :: fluxsurf_sw(ngrid,2), fluxtop_sw(ngrid,2)
175
176	REAL,INTENT(OUT) :: tauref(ngrid), tau(ngrid,naerkind)
177	REAL,INTENT(OUT) :: taucloudtes(ngrid)! Cloud opacity at infrared
178	! reference wavelength using
179	! Qabs instead of Qext
180	! (direct comparison with TES)
181	REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind)
182	REAL,INTENT(OUT) :: rdust(ngrid,nlayer) ! Dust geometric mean radius (m)
183	REAL,INTENT(OUT) :: rice(ngrid,nlayer) ! Ice geometric mean radius (m)
184	REAL,INTENT(OUT) :: nuice(ngrid,nlayer) ! Estimated effective variance
185	REAL,INTENT(IN) :: co2ice(ngrid) ! co2 ice surface layer (kg.m-2)
186
187	c
188	c Local variables :
189	c -----------------
190
191	INTEGER j,l,ig,n,ich
192	INTEGER aer_count,iaer
193	INTEGER jd,ig0,nd
194
195	real cste_mars ! solar constant on Mars (Wm-2)
196	REAL ptlev(ngrid,nlayer+1)
197
198	INTEGER,SAVE :: ndomain
199
200	c Thermal IR net radiative budget (W m-2)
201
202	real znetrad(ndomainsz,nflev)
203
204	real zfluxd_sw(ndomainsz,nflev+1,2)
205	real zfluxu_sw(ndomainsz,nflev+1,2)
206
207	REAL zplev(ndomainsz,nflev+1)
208	REAL zztlev(ndomainsz,nflev+1)
209	REAL zplay(ndomainsz,nflev)
210	REAL zt(ndomainsz,nflev)
211	REAL zaerosol(ndomainsz,nflev,naerkind)
212	REAL zalbedo(ndomainsz,2)
213	REAL zdp(ndomainsz,nflev)
214	REAL zdt0(ndomainsz)
215
216	REAL zzdtlw(ndomainsz,nflev)
217	REAL zzdtsw(ndomainsz,nflev)
218	REAL zzflux(ndomainsz,6)
219	real zrmuz
220
221	REAL :: zQVISsQREF3d(ndomainsz,nflev,nsun,naerkind)
222	REAL :: zomegaVIS3d(ndomainsz,nflev,nsun,naerkind)
223	REAL :: zgVIS3d(ndomainsz,nflev,nsun,naerkind)
224
225	REAL :: zQIRsQREF3d(ndomainsz,nflev,nir,naerkind)
226	REAL :: zomegaIR3d(ndomainsz,nflev,nir,naerkind)
227	REAL :: zgIR3d(ndomainsz,nflev,nir,naerkind)
228
229	c Aerosol size distribution
230	REAL :: reffrad(ngrid,nlayer,naerkind)
231	REAL :: nueffrad(ngrid,nlayer,naerkind)
232	c Aerosol optical properties
233	REAL :: QVISsQREF3d(ngrid,nlayer,nsun,naerkind)
234	REAL :: omegaVIS3d(ngrid,nlayer,nsun,naerkind)
235	REAL :: gVIS3d(ngrid,nlayer,nsun,naerkind)
236
237	REAL :: QIRsQREF3d(ngrid,nlayer,nir,naerkind)
238	REAL :: omegaIR3d(ngrid,nlayer,nir,naerkind)
239	REAL :: gIR3d(ngrid,nlayer,nir,naerkind)
240
241	REAL :: QREFvis3d(ngrid,nlayer,naerkind)
242	REAL :: QREFir3d(ngrid,nlayer,naerkind)
243
244	REAL :: omegaREFvis3d(ngrid,nlayer,naerkind)
245	REAL :: omegaREFir3d(ngrid,nlayer,naerkind)
246
247	c local saved variables
248	c ---------------------
249
250	real,save,allocatable :: pview(:)
251
252	real zco2 ! volume fraction of CO2 in Mars atmosphere
253	DATA zco2/0.95/
254	SAVE zco2
255
256	LOGICAL firstcall
257	DATA firstcall/.true./
258	SAVE firstcall
259
260	c----------------------------------------------------------------------
261
262	c Initialisation
263	c --------------
264
265	IF (firstcall) THEN
266	! compute ndomain and allocate local saved arrays
267	ndomain= (ngrid-1) / ndomainsz + 1
268	allocate(pview(ngrid))
269
270	c Assign a number to the different scatterers
271	c -------------------------------------------
272
273	iaer_dust_conrath=0
274	iaer_dust_doubleq=0
275	iaer_dust_submicron=0
276	iaer_h2o_ice=0
277
278	aer_count=0
279	if (.NOT.active) then
280	do iaer=1,naerkind
281	if (name_iaer(iaer).eq."dust_conrath") then
282	iaer_dust_conrath = iaer
283	aer_count = aer_count + 1
284	endif
285	enddo
286	endif
287	if (doubleq.AND.active) then
288	do iaer=1,naerkind
289	if (name_iaer(iaer).eq."dust_doubleq") then
290	iaer_dust_doubleq = iaer
291	aer_count = aer_count + 1
292	endif
293	enddo
294	endif
295	if (submicron.AND.active) then
296	do iaer=1,naerkind
297	if (name_iaer(iaer).eq."dust_submicron") then
298	iaer_dust_submicron = iaer
299	aer_count = aer_count + 1
300	endif
301	enddo
302	endif
303	if (water.AND.activice) then
304	do iaer=1,naerkind
305	if (name_iaer(iaer).eq."h2o_ice") then
306	iaer_h2o_ice = iaer
307	aer_count = aer_count + 1
308	endif
309	enddo
310	endif
311
312	c Check that we identified all tracers:
313	if (aer_count.ne.naerkind) then
314	write(,) "callradite: found only ",aer_count," scatterers"
315	write(,) " expected ",naerkind
316	write(,) "please make sure that the number of"
317	write(,) "scatterers in scatterers.h, the names"
318	write(,) "in callradite.F, and the flags in"
319	write(,) "callphys.def are all consistent!"
320	do iaer=1,naerkind
321	write(,)' ',iaer,' ',trim(name_iaer(iaer))
322	enddo
323	stop
324	else
325	write(,) "callradite: found all scatterers, namely:"
326	do iaer=1,naerkind
327	write(,)' ',iaer,' ',trim(name_iaer(iaer))
328	enddo
329	endif
330	c -------------------------------------------
331
332	DO ig=1,ngrid
333	pview(ig)=1.66 ! cosecant of viewing angle
334	ENDDO
335	gcp = g/cpp
336
337	c Loading the optical properties in external look-up tables:
338	CALL SUAER
339	! CALL SULW ! this step is now done in ini_yomlw_h
340
341	write(,) 'Splitting radiative calculations: ',
342	$ ' ngrid,ndomainsz,ndomain',
343	$ ngrid,ndomainsz,ndomain
344	if (ngrid .EQ. 1) then
345	if (ndomainsz .NE. 1) then
346	print*
347	print*,'ATTENTION !!!'
348	print*,'pour tourner en 1D, '
349	print*,'fixer ndomainsz=1 dans phymars/dimradmars_h'
350	print*
351	call exit(1)
352	endif
353	endif
354	firstcall=.false.
355	END IF
356
357	c Computing aerosol optical properties and opacity
358	c ------------------------------------------------
359
360	c Updating aerosol size distributions:
361	CALL updatereffrad(ngrid,nlayer,
362	& rdust,rice,nuice,
363	& reffrad,nueffrad,
364	& pq,tauscaling,tau,pplay)
365
366	c Computing 3D scattering parameters:
367	CALL aeroptproperties(ngrid,nlayer,reffrad,nueffrad,
368	& QVISsQREF3d,omegaVIS3d,gVIS3d,
369	& QIRsQREF3d,omegaIR3d,gIR3d,
370	& QREFvis3d,QREFir3d,
371	& omegaREFvis3d,omegaREFir3d)
372
373	c Computing aerosol optical depth in each layer:
374	CALL aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
375	& pq,tauscaling,tauref,tau,taucloudtes,aerosol,reffrad,
376	& nueffrad,QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d)
377
378	c Starting loop on sub-domain
379	c ----------------------------
380
381	DO jd=1,ndomain
382	ig0=(jd-1)*ndomainsz
383	if (jd.eq.ndomain) then
384	nd=ngrid-ig0
385	else
386	nd=ndomainsz
387	endif
388
389	c Spliting input variable in sub-domain input variables
390	c ---------------------------------------------------
391
392	do l=1,nlaylte
393	do ig = 1,nd
394	do iaer = 1, naerkind
395	do ich = 1, nsun
396	zQVISsQREF3d(ig,l,ich,iaer) =
397	& QVISsQREF3d(ig0+ig,l,ich,iaer)
398	zomegaVIS3d(ig,l,ich,iaer) =
399	& omegaVIS3d(ig0+ig,l,ich,iaer)
400	zgVIS3d(ig,l,ich,iaer) =
401	& gVIS3d(ig0+ig,l,ich,iaer)
402	enddo
403	do ich = 1, nir
404	zQIRsQREF3d(ig,l,ich,iaer) =
405	& QIRsQREF3d(ig0+ig,l,ich,iaer)
406	zomegaIR3d(ig,l,ich,iaer) =
407	& omegaIR3d(ig0+ig,l,ich,iaer)
408	zgIR3d(ig,l,ich,iaer) =
409	& gIR3d(ig0+ig,l,ich,iaer)
410	enddo
411	enddo
412	enddo
413	enddo
414
415	do l=1,nlaylte+1
416	do ig = 1,nd
417	zplev(ig,l) = pplev(ig0+ig,l)
418	enddo
419	enddo
420
421	do l=1,nlaylte
422	do ig = 1,nd
423	zplay(ig,l) = pplay(ig0+ig,l)
424	zt(ig,l) = pt(ig0+ig,l)
425	c Thickness of each layer (Pa) :
426	zdp(ig,l)= pplev(ig0+ig,l) - pplev(ig0+ig,l+1)
427	enddo
428	enddo
429
430	do n=1,naerkind
431	do l=1,nlaylte
432	do ig=1,nd
433	zaerosol(ig,l,n) = aerosol(ig0+ig,l,n)
434	enddo
435	enddo
436	enddo
437
438	do j=1,2
439	do ig = 1,nd
440	zalbedo(ig,j) = albedo(ig0+ig,j)
441	enddo
442	enddo
443
444	c Intermediate levels: (computing tlev)
445	c ---------------------------------------
446	c Extrapolation for the air temperature above the surface
447	DO ig=1,nd
448	zztlev(ig,1)=zt(ig,1)+
449	s (zplev(ig,1)-zplay(ig,1))*
450	s (zt(ig,1)-zt(ig,2))/(zplay(ig,1)-zplay(ig,2))
451
452	zdt0(ig) = tsurf(ig0+ig) - zztlev(ig,1)
453	ENDDO
454
455	DO l=2,nlaylte
456	DO ig=1, nd
457	zztlev(ig,l)=0.5*(zt(ig,l-1)+zt(ig,l))
458	ENDDO
459	ENDDO
460
461	DO ig=1, nd
462	zztlev(ig,nlaylte+1)=zt(ig,nlaylte)
463	ENDDO
464
465
466	c Longwave ("lw") radiative transfer (= thermal infrared)
467	c -------------------------------------------------------
468	call lwmain (ig0,icount,nd,nflev
469	. ,zdp,zdt0,emis(ig0+1),zplev,zztlev,zt
470	. ,zaerosol,zzdtlw
471	. ,fluxsurf_lw(ig0+1),fluxtop_lw(ig0+1)
472	. ,znetrad
473	& ,zQIRsQREF3d,zomegaIR3d,zgIR3d
474	& ,co2ice(ig0+1))
475
476	c Shortwave ("sw") radiative transfer (= solar radiation)
477	c -------------------------------------------------------
478	c Mars solar constant (W m-2)
479	c 1370 W.m-2 is the solar constant at 1 AU.
480	cste_mars=1370./(dist_sol*dist_sol)
481
482	call swmain ( nd, nflev,
483	S cste_mars, zalbedo,
484	S mu0(ig0+1), zdp, zplev, zaerosol, fract(ig0+1),
485	S zzdtsw, zfluxd_sw, zfluxu_sw,
486	& zQVISsQREF3d,zomegaVIS3d,zgVIS3d)
487
488	c ------------------------------------------------------------
489	c Un-spliting output variable from sub-domain input variables
490	c ------------------------------------------------------------
491
492	do l=1,nlaylte
493	do ig = 1,nd
494	dtlw(ig0+ig,l) = zzdtlw(ig,l)
495	dtsw(ig0+ig,l) = zzdtsw(ig,l)
496	enddo
497	enddo
498
499	do l=1,nlaylte+1
500	do ig = 1,nd
501	ptlev(ig0+ig,l) = zztlev(ig,l)
502	enddo
503	enddo
504
505	do ig = 1,nd
506	fluxsurf_sw(ig0+ig,1) = zfluxd_sw(ig,1,1)
507	fluxsurf_sw(ig0+ig,2) = zfluxd_sw(ig,1,2)
508	fluxtop_sw(ig0+ig,1) = zfluxu_sw(ig,nlaylte+1,1)
509	fluxtop_sw(ig0+ig,2) = zfluxu_sw(ig,nlaylte+1,2)
510	enddo
511
512	ENDDO ! (boucle jd=1, ndomain)
513
514	c Zero tendencies for any remaining layers between nlaylte and nlayer
515	if (nlayer.gt.nlaylte) then
516	do l = nlaylte+1, nlayer
517	do ig = 1, ngrid
518	dtlw(ig, l) = 0.
519	dtsw(ig, l) = 0.
520	enddo
521	enddo
522	endif
523
524	c Output for debugging if lwrite=T
525	c --------------------------------
526	c Write all nlayer layers, even though only nlaylte layers may have
527	c non-zero tendencies.
528
529	IF(lwrite) THEN
530	PRINT*,'Diagnotique for the radiation'
531	PRINT*,'albedo, emissiv, mu0,fract,fluxsurf_lw,fluxsurf_sw'
532	PRINT*,albedo(igout,1),emis(igout),mu0(igout),
533	s fract(igout), fluxsurf_lw(igout),
534	$ fluxsurf_sw(igout,1)+fluxsurf_sw(igout,2)
535	PRINT*,'Tlay Tlev Play Plev dT/dt SW dT/dt LW (K/s)'
536	PRINT*,'daysec',daysec
537	DO l=1,nlayer
538	PRINT*,pt(igout,l),ptlev(igout,l),
539	s pplay(igout,l),pplev(igout,l),
540	s dtsw(igout,l),dtlw(igout,l)
541	ENDDO
542	ENDIF
543
544
545	return
546	end

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

Download in other formats:

Original Format