Context Navigation

callcorrk.F90 @ 848

Last change on this file since 848 was 787, checked in by aslmd, 12 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
RUN RCM1D : perfect match in normal mode.
RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

Additional arguments to some subroutines (ngrid and nq)
F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this

Property svn:executable set to *

File size: 29.0 KB

Line
1	subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf, &
2	albedo,emis,mu0,pplev,pplay,pt, &
3	tsurf,fract,dist_star,aerosol,muvar, &
4	dtlw,dtsw,fluxsurf_lw, &
5	fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn, &
6	OLR_nu,OSR_nu, &
7	reffrad,nueffrad,tau_col,cloudfrac,totcloudfrac, &
8	clearsky,firstcall,lastcall)
9
10	use radinc_h
11	use radcommon_h
12	use watercommon_h
13	use datafile_mod, only: datadir
14	use ioipsl_getincom
15	use gases_h
16	use radii_mod, only : su_aer_radii,co2_reffrad,h2o_reffrad,dust_reffrad,h2so4_reffrad
17	use aerosol_mod
18
19	USE tracer_h
20
21	implicit none
22
23	!==================================================================
24	!
25	! Purpose
26	! -------
27	! Solve the radiative transfer using the correlated-k method for
28	! the gaseous absorption and the Toon et al. (1989) method for
29	! scatttering due to aerosols.
30	!
31	! Authors
32	! -------
33	! Emmanuel 01/2001, Forget 09/2001
34	! Robin Wordsworth (2009)
35	!
36	!==================================================================
37
38	#include "dimphys.h"
39	#include "comcstfi.h"
40	#include "callkeys.h"
41
42	!-----------------------------------------------------------------------
43	! Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
44	! Layer #1 is the layer near the ground.
45	! Layer #nlayermx is the layer at the top.
46
47	! INPUT
48	INTEGER icount
49	INTEGER ngrid,nlayer
50	REAL aerosol(ngrid,nlayermx,naerkind) ! aerosol tau (kg/kg)
51	REAL albedo(ngrid) ! SW albedo
52	REAL emis(ngrid) ! LW emissivity
53	REAL pplay(ngrid,nlayermx) ! pres. level in GCM mid of layer
54	REAL pplev(ngrid,nlayermx+1) ! pres. level at GCM layer boundaries
55
56	REAL pt(ngrid,nlayermx) ! air temperature (K)
57	REAL tsurf(ngrid) ! surface temperature (K)
58	REAL dist_star,mu0(ngrid) ! distance star-planet (AU)
59	REAL fract(ngrid) ! fraction of day
60
61	! Globally varying aerosol optical properties on GCM grid
62	! Not needed everywhere so not in radcommon_h
63	REAL :: QVISsQREF3d(ngrid,nlayermx,L_NSPECTV,naerkind)
64	REAL :: omegaVIS3d(ngrid,nlayermx,L_NSPECTV,naerkind)
65	REAL :: gVIS3d(ngrid,nlayermx,L_NSPECTV,naerkind)
66
67	REAL :: QIRsQREF3d(ngrid,nlayermx,L_NSPECTI,naerkind)
68	REAL :: omegaIR3d(ngrid,nlayermx,L_NSPECTI,naerkind)
69	REAL :: gIR3d(ngrid,nlayermx,L_NSPECTI,naerkind)
70
71	! REAL :: omegaREFvis3d(ngrid,nlayermx,naerkind)
72	! REAL :: omegaREFir3d(ngrid,nlayermx,naerkind) ! not sure of the point of these...
73
74	!!! this is a local instance of a variable saved in physiq.F and being an argument
75	REAL, INTENT(INOUT) :: reffrad(ngrid,nlayer,naerkind)
76	REAL, INTENT(INOUT) :: nueffrad(ngrid,nlayermx,naerkind)
77
78	!! OUTPUT
79	REAL dtsw(ngrid,nlayermx) ! heating rate (K/s) due to SW
80	REAL dtlw(ngrid,nlayermx) ! heating rate (K/s) due to LW
81	REAL fluxsurf_lw(ngrid) ! incident LW flux to surf (W/m2)
82	REAL fluxtop_lw(ngrid) ! outgoing LW flux to space (W/m2)
83	REAL fluxsurf_sw(ngrid) ! incident SW flux to surf (W/m2)
84	REAL fluxabs_sw(ngrid) ! SW flux absorbed by planet (W/m2)
85	REAL fluxtop_dn(ngrid) ! incident top of atmosphere SW flux (W/m2)
86	REAL OLR_nu(ngrid,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
87	REAL OSR_nu(ngrid,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
88
89	!-----------------------------------------------------------------------
90	! Declaration of the variables required by correlated-k subroutines
91	! Numbered from top to bottom unlike in the GCM!
92
93	REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
94	REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
95
96	! Optical values for the optci/cv subroutines
97	REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
98	REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
99	REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
100	REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
101	REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
102	REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
103	REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
104	REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
105	REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
106	REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
107
108	REAL*8 tauaero(L_LEVELS+1,naerkind)
109	REAL*8 nfluxtopv,nfluxtopi,nfluxtop
110	real*8 nfluxoutv_nu(L_NSPECTV) ! outgoing band-resolved VI flux at TOA (W/m2)
111	real*8 nfluxtopi_nu(L_NSPECTI) ! net band-resolved IR flux at TOA (W/m2)
112	real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! for 1D diagnostic
113	REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
114	REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
115	REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
116	REAL*8 albi,albv,acosz
117
118	INTEGER ig,l,k,nw,iaer,irad
119
120	real fluxtoplanet
121	real szangle
122	logical global1d
123	save szangle,global1d
124	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
125	real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
126
127	real*8 qvar(L_LEVELS) ! mixing ratio of variable component (mol/mol)
128	REAL pq(ngrid,nlayer,nq)
129	REAL qsurf(ngrid,nq) ! tracer on surface (e.g. kg.m-2)
130	integer nq
131
132	! Local aerosol optical properties for each column on RADIATIVE grid
133	real*8 QXVAER(L_LEVELS+1,L_NSPECTV,naerkind)
134	real*8 QSVAER(L_LEVELS+1,L_NSPECTV,naerkind)
135	real*8 GVAER(L_LEVELS+1,L_NSPECTV,naerkind)
136	real*8 QXIAER(L_LEVELS+1,L_NSPECTI,naerkind)
137	real*8 QSIAER(L_LEVELS+1,L_NSPECTI,naerkind)
138	real*8 GIAER(L_LEVELS+1,L_NSPECTI,naerkind)
139
140	save qxvaer, qsvaer, gvaer
141	save qxiaer, qsiaer, giaer
142
143	!REAL :: QREFvis3d(ngrid,nlayermx,naerkind)
144	!REAL :: QREFir3d(ngrid,nlayermx,naerkind)
145	!save QREFvis3d, QREFir3d
146	real, dimension(:,:,:), save, allocatable :: QREFvis3d
147	real, dimension(:,:,:), save, allocatable :: QREFir3d
148
149	REAL tau_col(ngrid) ! diagnostic from aeropacity
150
151	! Misc.
152	logical firstcall, lastcall, nantest
153	real*8 tempv(L_NSPECTV)
154	real*8 tempi(L_NSPECTI)
155	real*8 temp,temp1,temp2,pweight
156	character(len=10) :: tmp1
157	character(len=10) :: tmp2
158
159	! for fixed water vapour profiles
160	integer i_var
161	real RH
162	real*8 pq_temp(nlayer)
163	real ptemp, Ttemp, qsat
164
165	! real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!
166
167	!real ptime, pday
168	logical OLRz
169	real*8 NFLUXGNDV_nu(L_NSPECTV)
170
171	! for H2O cloud fraction in aeropacity
172	real cloudfrac(ngrid,nlayermx)
173	real totcloudfrac(ngrid)
174	logical clearsky
175
176	! for weird cloud test
177	real pqtest(ngrid,nlayer,nq)
178
179	real maxrad, minrad
180
181	real CBRT
182	external CBRT
183
184	! included by RW for runaway greenhouse 1D study
185	real muvar(ngrid,nlayermx+1)
186	real vtmp(nlayermx)
187	REAL*8 muvarrad(L_LEVELS)
188
189
190	!===============================================================
191	! Initialization on first call
192
193	qxvaer(:,:,:)=0.0
194	qsvaer(:,:,:)=0.0
195	gvaer(:,:,:) =0.0
196
197	qxiaer(:,:,:)=0.0
198	qsiaer(:,:,:)=0.0
199	giaer(:,:,:) =0.0
200
201	if(firstcall) then
202
203	ALLOCATE(QREFvis3d(ngrid,nlayermx,naerkind))
204	ALLOCATE(QREFir3d(ngrid,nlayermx,naerkind))
205
206	call system('rm -f surf_vals_long.out')
207
208	!--------------------------------------------------
209	! Effective radius and variance of the aerosols
210	if(naerkind.gt.4)then
211	print*,'Code not general enough to deal with naerkind > 4 yet.'
212	call abort
213	endif
214	call su_aer_radii(ngrid,reffrad,nueffrad)
215
216
217
218	!--------------------------------------------------
219	! set up correlated k
220	print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
221	call getin("corrkdir",corrkdir)
222	print*, "corrkdir = ",corrkdir
223	write( tmp1, '(i3)' ) L_NSPECTI
224	write( tmp2, '(i3)' ) L_NSPECTV
225	banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
226	banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
227
228	call sugas_corrk ! set up gaseous absorption properties
229	call setspi ! basic infrared properties
230	call setspv ! basic visible properties
231	call suaer_corrk ! set up aerosol optical properties
232
233	Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed
234
235	if((igcm_h2o_vap.eq.0) .and. varactive)then
236	print*,'varactive in callcorrk but no h2o_vap tracer.'
237	stop
238	endif
239
240	OLR_nu(:,:) = 0.
241	OSR_nu(:,:) = 0.
242
243	if (ngrid.eq.1) then
244	PRINT*, 'Simulate global averaged conditions ?'
245	global1d = .false. ! default value
246	call getin("global1d",global1d)
247	write(,) "global1d = ",global1d
248	! Test of incompatibility:
249	! if global1d is true, there should not be any diurnal cycle
250	if (global1d.and.diurnal) then
251	print*,'if global1d is true, diurnal must be set to false'
252	stop
253	endif
254
255	if (global1d) then
256	PRINT *,'Solar Zenith angle (deg.) ?'
257	PRINT *,'(assumed for averaged solar flux S/4)'
258	szangle=60.0 ! default value
259	call getin("szangle",szangle)
260	write(,) "szangle = ",szangle
261	endif
262	endif
263
264	firstcall=.false.
265
266	end if
267
268	!=======================================================================
269	! Initialization on every call
270
271	!--------------------------------------------------
272	! Effective radius and variance of the aerosols
273	do iaer=1,naerkind
274
275	if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! treat condensed co2 particles.
276	call co2_reffrad(ngrid,nq,pq,reffrad)
277	print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
278	print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
279	end if
280	if ((iaer.eq.iaero_h2o).and.water) then ! treat condensed water particles. to be generalized for other aerosols
281	call h2o_reffrad(ngrid,nq,pq,pt,reffrad,nueffrad)
282	print*,'Max. H2O cloud particle size = ',maxval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
283	print*,'Min. H2O cloud particle size = ',minval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
284	endif
285	if(iaer.eq.iaero_dust)then
286	call dust_reffrad(ngrid,reffrad)
287	print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
288	endif
289	if(iaer.eq.iaero_h2so4)then
290	call h2so4_reffrad(ngrid,reffrad)
291	print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
292	endif
293	end do !iaer=1,naerkind
294
295
296	! how much light we get
297	fluxtoplanet=0
298	do nw=1,L_NSPECTV
299	stel(nw)=stellarf(nw)/(dist_star**2)
300	fluxtoplanet=fluxtoplanet + stel(nw)
301	end do
302
303	call aeroptproperties(ngrid,nlayer,reffrad,nueffrad, &
304	QVISsQREF3d,omegaVIS3d,gVIS3d, &
305	QIRsQREF3d,omegaIR3d,gIR3d, &
306	QREFvis3d,QREFir3d) ! get 3D aerosol optical properties
307
308	call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol, &
309	reffrad,QREFvis3d,QREFir3d, &
310	tau_col,cloudfrac,totcloudfrac,clearsky) ! get aerosol optical depths
311
312	!-----------------------------------------------------------------------
313	! Starting Big Loop over every GCM column
314	do ig=1,ngrid
315
316	!=======================================================================
317	! Transformation of the GCM variables
318
319	!-----------------------------------------------------------------------
320	! Aerosol optical properties Qext, Qscat and g
321	! The transformation in the vertical is the same as for temperature
322
323	! shortwave
324	do iaer=1,naerkind
325	DO nw=1,L_NSPECTV
326	do l=1,nlayermx
327
328	temp1=QVISsQREF3d(ig,nlayermx+1-l,nw,iaer) &
329	*QREFvis3d(ig,nlayermx+1-l,iaer)
330
331	temp2=QVISsQREF3d(ig,max(nlayermx-l,1),nw,iaer) &
332	*QREFvis3d(ig,max(nlayermx-l,1),iaer)
333
334	qxvaer(2*l,nw,iaer) = temp1
335	qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
336
337	temp1=temp1*omegavis3d(ig,nlayermx+1-l,nw,iaer)
338	temp2=temp2*omegavis3d(ig,max(nlayermx-l,1),nw,iaer)
339
340	qsvaer(2*l,nw,iaer) = temp1
341	qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
342
343	temp1=gvis3d(ig,nlayermx+1-l,nw,iaer)
344	temp2=gvis3d(ig,max(nlayermx-l,1),nw,iaer)
345
346	gvaer(2*l,nw,iaer) = temp1
347	gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
348
349	end do
350
351	qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
352	qxvaer(2*nlayermx+1,nw,iaer)=0.
353
354	qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
355	qsvaer(2*nlayermx+1,nw,iaer)=0.
356
357	gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
358	gvaer(2*nlayermx+1,nw,iaer)=0.
359
360	end do
361
362	! longwave
363	DO nw=1,L_NSPECTI
364	do l=1,nlayermx
365
366	temp1=QIRsQREF3d(ig,nlayermx+1-l,nw,iaer) &
367	*QREFir3d(ig,nlayermx+1-l,iaer)
368
369	temp2=QIRsQREF3d(ig,max(nlayermx-l,1),nw,iaer) &
370	*QREFir3d(ig,max(nlayermx-l,1),iaer)
371
372	qxiaer(2*l,nw,iaer) = temp1
373	qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
374
375	temp1=temp1*omegair3d(ig,nlayermx+1-l,nw,iaer)
376	temp2=temp2*omegair3d(ig,max(nlayermx-l,1),nw,iaer)
377
378	qsiaer(2*l,nw,iaer) = temp1
379	qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
380
381	temp1=gir3d(ig,nlayermx+1-l,nw,iaer)
382	temp2=gir3d(ig,max(nlayermx-l,1),nw,iaer)
383
384	giaer(2*l,nw,iaer) = temp1
385	giaer(2*l+1,nw,iaer)=(temp1+temp2)/2
386
387	end do
388
389	qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
390	qxiaer(2*nlayermx+1,nw,iaer)=0.
391
392	qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
393	qsiaer(2*nlayermx+1,nw,iaer)=0.
394
395	giaer(1,nw,iaer)=giaer(2,nw,iaer)
396	giaer(2*nlayermx+1,nw,iaer)=0.
397
398	end do
399	end do
400
401	! test / correct for freaky s. s. albedo values
402	do iaer=1,naerkind
403	do k=1,L_LEVELS+1
404
405	do nw=1,L_NSPECTV
406	if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
407	print*,'Serious problems with qsvaer values'
408	print*,'in callcorrk'
409	call abort
410	endif
411	if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
412	qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
413	endif
414	end do
415
416	do nw=1,L_NSPECTI
417	if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
418	print*,'Serious problems with qsiaer values'
419	print*,'in callcorrk'
420	call abort
421	endif
422	if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then
423	qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer)
424	endif
425	end do
426
427	end do
428	end do
429
430	!-----------------------------------------------------------------------
431	! Aerosol optical depths
432
433	do iaer=1,naerkind ! a bug was here
434	do k=0,nlayer-1
435
436	pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/ &
437	(pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
438
439	temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
440
441	tauaero(2k+2,iaer)=max(temppweight,0.d0)
442	tauaero(2k+3,iaer)=max(temp-tauaero(2k+2,iaer),0.d0)
443	!
444	end do
445	! boundary conditions
446	tauaero(1,iaer) = tauaero(2,iaer)
447	tauaero(L_LEVELS+1,iaer) = tauaero(L_LEVELS,iaer)
448	!tauaero(1,iaer) = 0.
449	!tauaero(L_LEVELS+1,iaer) = 0.
450	end do
451
452	! Albedo and emissivity
453	albi=1-emis(ig) ! longwave
454	albv=albedo(ig) ! shortwave
455
456	if(noradsurf.and.(albv.gt.0.0))then
457	print*,'For open lower boundary in callcorrk must'
458	print*,'have surface albedo set to zero!'
459	call abort
460	endif
461
462	if ((ngrid.eq.1).and.(global1d)) then ! fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight
463	acosz = cos(pi*szangle/180.0)
464	print*,'acosz=',acosz,', szangle=',szangle
465	else
466	acosz=mu0(ig) ! cosine of sun incident angle : 3D simulations or local 1D simulations using latitude
467	endif
468
469	!-----------------------------------------------------------------------
470	! Water vapour (to be generalised for other gases eventually)
471
472	if(varactive)then
473
474	i_var=igcm_h2o_vap
475	do l=1,nlayer
476	qvar(2*l) = pq(ig,nlayer+1-l,i_var)
477	qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2
478	! Average approximation as for temperature...
479	end do
480	qvar(1)=qvar(2)
481
482	elseif(varfixed)then
483
484	do l=1,nlayermx ! here we will assign fixed water vapour profiles globally
485	RH = satval * ((pplay(ig,l)/pplev(ig,1) - 0.02) / 0.98)
486	if(RH.lt.0.0) RH=0.0
487
488	ptemp=pplay(ig,l)
489	Ttemp=pt(ig,l)
490	call watersat(Ttemp,ptemp,qsat)
491
492	!pq_temp(l) = qsat ! fully saturated everywhere
493	pq_temp(l) = RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
494	end do
495
496	do l=1,nlayer
497	qvar(2*l) = pq_temp(nlayer+1-l)
498	qvar(2*l+1) = (pq_temp(nlayer+1-l)+pq_temp(max(nlayer-l,1)))/2
499	end do
500	qvar(1)=qvar(2)
501
502	! Lowest layer of atmosphere
503	RH = satval * (1 - 0.02) / 0.98
504	if(RH.lt.0.0) RH=0.0
505
506	ptemp = pplev(ig,1)
507	Ttemp = tsurf(ig)
508	call watersat(Ttemp,ptemp,qsat)
509
510	!qvar(2*nlayermx+1)=qsat ! fully saturated everywhere
511	qvar(2nlayermx+1)= RH qsat ! ~realistic profile (e.g. 80% saturation at ground)
512	!qvar=0.005 ! completely constant profile (JL)
513
514	else
515	do k=1,L_LEVELS
516	qvar(k) = 1.0D-7
517	end do
518	end if
519
520	if(.not.kastprof)then
521	! IMPORTANT: Now convert from kg/kg to mol/mol
522	do k=1,L_LEVELS
523	qvar(k) = qvar(k)/(epsi+qvar(k)*(1.-epsi))
524	end do
525	end if
526
527	!-----------------------------------------------------------------------
528	! kcm mode only
529	if(kastprof)then
530
531	! initial values equivalent to mugaz
532	DO l=1,nlayer
533	muvarrad(2*l) = mugaz
534	muvarrad(2*l+1) = mugaz
535	END DO
536
537	!do k=1,L_LEVELS
538	! qvar(k) = 0.0
539	!end do
540	!print*,'ASSUMING qH2O=0 EVERYWHERE IN CALLCORRK!'
541	endif
542
543
544	if(kastprof.and.(ngasmx.gt.1))then
545
546	DO l=1,nlayer
547	muvarrad(2*l) = muvar(ig,nlayer+2-l)
548	muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + &
549	muvar(ig,max(nlayer+1-l,1)))/2
550	END DO
551
552	muvarrad(1) = muvarrad(2)
553	muvarrad(2*nlayermx+1)=muvar(ig,1)
554
555	print*,'Recalculating qvar with VARIABLE epsi for kastprof'
556	print*,'Assumes that the variable gas is H2O!!!'
557	print*,'Assumes that there is only one tracer'
558	!i_var=igcm_h2o_vap
559	i_var=1
560	if(nq.gt.1)then
561	print*,'Need 1 tracer only to run kcm1d.e'
562	stop
563	endif
564	do l=1,nlayer
565	vtmp(l)=pq(ig,l,i_var)*muvar(ig,l+1)/mH2O
566	end do
567
568	do l=1,nlayer
569	qvar(2*l) = vtmp(nlayer+1-l)
570	qvar(2*l+1) = ( vtmp(nlayer+1-l) + vtmp(max(nlayer-l,1)) )/2
571	end do
572	qvar(1)=qvar(2)
573
574	print*,'Warning: reducing qvar in callcorrk.'
575	print*,'Temperature profile no longer consistent ', &
576	'with saturated H2O.'
577	do k=1,L_LEVELS
578	qvar(k) = qvar(k)*satval
579	end do
580
581	endif
582
583	! Keep values inside limits for which we have radiative transfer coefficients
584	if(L_REFVAR.gt.1)then ! there was a bug here!
585	do k=1,L_LEVELS
586	if(qvar(k).lt.wrefvar(1))then
587	qvar(k)=wrefvar(1)+1.0e-8
588	elseif(qvar(k).gt.wrefvar(L_REFVAR))then
589	qvar(k)=wrefvar(L_REFVAR)-1.0e-8
590	endif
591	end do
592	endif
593
594	!-----------------------------------------------------------------------
595	! Pressure and temperature
596
597	DO l=1,nlayer
598	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
599	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
600	tlevrad(2*l) = pt(ig,nlayer+1-l)
601	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
602	END DO
603
604	plevrad(1) = 0.
605	plevrad(2) = max(pgasmin,0.0001*plevrad(3))
606
607	tlevrad(1) = tlevrad(2)
608	tlevrad(2*nlayermx+1)=tsurf(ig)
609
610	tmid(1) = tlevrad(2)
611	tmid(2) = tlevrad(2)
612	pmid(1) = plevrad(2)
613	pmid(2) = plevrad(2)
614
615	DO l=1,L_NLAYRAD-1
616	tmid(2l+1) = tlevrad(2l+1)
617	tmid(2l+2) = tlevrad(2l+1)
618	pmid(2l+1) = plevrad(2l+1)
619	pmid(2l+2) = plevrad(2l+1)
620	END DO
621	pmid(L_LEVELS) = plevrad(L_LEVELS)
622	tmid(L_LEVELS) = tlevrad(L_LEVELS)
623
624	! test for out-of-bounds pressure
625	if(plevrad(3).lt.pgasmin)then
626	print*,'Minimum pressure is outside the radiative'
627	print*,'transfer kmatrix bounds, exiting.'
628	call abort
629	elseif(plevrad(L_LEVELS).gt.pgasmax)then
630	print*,'Maximum pressure is outside the radiative'
631	print*,'transfer kmatrix bounds, exiting.'
632	call abort
633	endif
634
635	! test for out-of-bounds temperature
636	do k=1,L_LEVELS
637	if(tlevrad(k).lt.tgasmin)then
638	print*,'Minimum temperature is outside the radiative'
639	print*,'transfer kmatrix bounds, exiting.'
640	!print*,'WARNING, OVERRIDING FOR TEST'
641	call abort
642	elseif(tlevrad(k).gt.tgasmax)then
643	print*,'Maximum temperature is outside the radiative'
644	print*,'transfer kmatrix bounds, exiting.'
645	print*,'level,grid,T',k,ig,tlevrad(k)
646	print*,'WARNING, OVERRIDING FOR TEST'
647	!call abort
648	tlevrad(k)=tgasmax
649	endif
650	enddo
651
652	!=======================================================================
653	! Calling the main radiative transfer subroutines
654
655
656	!-----------------------------------------------------------------------
657	! Shortwave
658
659	if(fract(ig) .ge. 1.0e-4) then ! only during daylight!
660
661	fluxtoplanet=0.
662
663	if((ngrid.eq.1).and.(global1d))then
664	do nw=1,L_NSPECTV
665	stel_fract(nw)= stel(nw) * 0.25 / acosz
666	fluxtoplanet=fluxtoplanet + stel_fract(nw)
667	! globally averaged = divide by 4
668	! but we correct for solar zenith angle
669	end do
670	else
671	do nw=1,L_NSPECTV
672	stel_fract(nw)= stel(nw) * fract(ig)
673	fluxtoplanet=fluxtoplanet + stel_fract(nw)
674	end do
675	endif
676
677	call optcv(dtauv,tauv,taucumv,plevrad, &
678	qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero, &
679	tmid,pmid,taugsurf,qvar,muvarrad)
680
681	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, &
682	acosz,stel_fract,gweight, &
683	nfluxtopv,nfluxoutv_nu,nfluxgndv_nu, &
684	fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
685
686	else ! during the night, fluxes = 0
687	nfluxtopv = 0.0
688	nfluxoutv_nu(:) = 0.0
689	nfluxgndv_nu(:) = 0.0
690	do l=1,L_NLAYRAD
691	fmnetv(l)=0.0
692	fluxupv(l)=0.0
693	fluxdnv(l)=0.0
694	end do
695	end if
696
697	!-----------------------------------------------------------------------
698	! Longwave
699
700	call optci(plevrad,tlevrad,dtaui,taucumi, &
701	qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid, &
702	taugsurfi,qvar,muvarrad)
703
704	call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, &
705	wnoi,dwni,cosbi,wbari,gweight,nfluxtopi,nfluxtopi_nu, &
706	fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
707
708	!-----------------------------------------------------------------------
709	! Transformation of the correlated-k code outputs
710	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
711
712	! Flux incident at the top of the atmosphere
713	fluxtop_dn(ig)=fluxdnv(1)
714
715	fluxtop_lw(ig) = real(nfluxtopi)
716	fluxabs_sw(ig) = real(-nfluxtopv)
717	fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
718	fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
719
720	if(fluxtop_dn(ig).lt.0.0)then
721	print*,'Achtung! fluxtop_dn has lost the plot!'
722	print*,'fluxtop_dn=',fluxtop_dn(ig)
723	print*,'acosz=',acosz
724	print*,'aerosol=',aerosol(ig,:,:)
725	print*,'temp= ',pt(ig,:)
726	print*,'pplay= ',pplay(ig,:)
727	call abort
728	endif
729
730	! Spectral output, for exoplanet observational comparison
731	if(specOLR)then
732	do nw=1,L_NSPECTI
733	OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
734	end do
735	do nw=1,L_NSPECTV
736	!GSR_nu(ig,nw)=nfluxgndv_nu(nw)
737	OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
738	end do
739	endif
740
741	! Finally, the heating rates
742
743	DO l=2,L_NLAYRAD
744	dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) &
745	g/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
746	dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) &
747	g/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
748	END DO
749
750	! These are values at top of atmosphere
751	dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) &
752	g/(cppscalep*(plevrad(3)-plevrad(1)))
753	dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) &
754	g/(cppscalep*(plevrad(3)-plevrad(1)))
755
756	! ---------------------------------------------------------------
757	end do ! end of big loop over every GCM column (ig = 1:ngrid)
758
759
760	!-----------------------------------------------------------------------
761	! Additional diagnostics
762
763	! IR spectral output, for exoplanet observational comparison
764
765
766	if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12
767
768	print*,'Saving scalar quantities in surf_vals.out...'
769	print*,'psurf = ', pplev(1,1),' Pa'
770	open(116,file='surf_vals.out')
771	write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), &
772	real(-nfluxtopv),real(nfluxtopi)
773	close(116)
774
775	! I am useful, please don`t remove me!
776	! if(specOLR)then
777	! open(117,file='OLRnu.out')
778	! do nw=1,L_NSPECTI
779	! write(117,*) OLR_nu(1,nw)
780	! enddo
781	! close(117)
782	!
783	! open(127,file='OSRnu.out')
784	! do nw=1,L_NSPECTV
785	! write(127,*) OSR_nu(1,nw)
786	! enddo
787	! close(127)
788	! endif
789
790	! OLR vs altitude: do it as a .txt file
791	OLRz=.false.
792	if(OLRz)then
793	print*,'saving IR vertical flux for OLRz...'
794	open(118,file='OLRz_plevs.out')
795	open(119,file='OLRz.out')
796	do l=1,L_NLAYRAD
797	write(118,) plevrad(2l)
798	do nw=1,L_NSPECTI
799	write(119,*) fluxupi_nu(l,nw)
800	enddo
801	enddo
802	close(118)
803	close(119)
804	endif
805
806	endif
807
808	! see physiq.F for explanations about CLFvarying. This is temporary.
809	if (lastcall .and. .not.CLFvarying) then
810	IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
811	IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
812	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
813	IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
814	IF( ALLOCATED( wrefvar ) ) DEALLOCATE( wrefvar )
815	IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
816	endif
817
818
819	end subroutine callcorrk

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

Download in other formats:

Original Format