Context Navigation

callradite.F @ 1009

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

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

Download in other formats:

Original Format