Context Navigation

recombin_corrk_mod.F90 @ 2550

Last change on this file since 2550 was 2543, checked in by aslmd, 4 years ago

Generic GCM:

Adding k-coefficients mixing on the fly
Working with MordernTrac?

JVO + YJ

File size: 45.1 KB

Line
1	MODULE recombin_corrk_mod
2
3	!
4	! Author : Jan Vatant d'Ollone (2018-2020)
5	!
6	! This module contains the following subroutines :
7	! - ini_recombin : From modern traceur.def options check if we will use recombining and for which species. Called by initracer.
8	! - su_recombin : Initialise tables. Called by sugas_corrk
9	! - call_recombin : Test profile of species and decide whether to call recombin_corrk. Called by callcork
10	! - recombin_corrk : The core algorithm properly recombining corrk tables. Called by callrecombin_corrk.
11	!
12
13	! TODO : Think about the case where nqtot phys /= nqtot_dyn !!
14	! TODO : Add the possibility to read an input profile for a specie !!
15	! Also think about the hybrid case where we want force with a latitudinally variable input profile (likely need to tweak on-the-fly)
16
17	IMPLICIT NONE
18
19	LOGICAL, SAVE :: corrk_recombin=.false. ! Key boolean, is there any specie to recombin.
20	LOGICAL, SAVE :: use_premix=.true. ! Indicates if we recombin on top of a existing premix set of corrk.
21	!$OMP THREADPRIVATE(corrk_recombin,use_premix)
22
23	INTEGER, SAVE :: nrecomb_tot ! # (total) of compounds to recombine
24	! -> Linked to key is_recomb in tracer_h
25
26	INTEGER, SAVE :: nrecomb_qset ! # of compounds to recombine having preset abundances (ie spectra computed with true vmr and not vmr=1)
27	!-> Linked to key is_recomb_qset in tracer_h
28
29	INTEGER, SAVE :: nrecomb_qotf ! # of compounds to recombine on-the-fly (ie using value of pq, not input profile)
30	!-> Linked to key is_recomb_qotf in tracer_h
31	!$OMP THREADPRIVATE(nrecomb_tot,nrecomb_qset,nrecomb_qotf)
32
33	! Note : The variables above are not read in callphys.def but automatically defined in inirecombin called in initracer after processing traceur.def
34
35	LOGICAL, SAVE :: all_otf ! True if all species are recombined on-the-fly, no premix, no preset ..
36	!$OMP THREADPRIVATE(all_otf)
37
38	! Following arrays are allocated in su_recombin (excepted otf2tot_idx, in ini_recombin) and deallocated in callcork lastcall
39	REAL*8, save, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi_recomb, gasv_recomb
40	REAL*8, save, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi_otf, gasv_otf
41	REAL*8, save, DIMENSION(:), ALLOCATABLE :: w_cum
42	REAL*8, save, DIMENSION(:), ALLOCATABLE :: wtwospec
43
44	INTEGER, save, DIMENSION(:), ALLOCATABLE :: otf2tot_idx
45	INTEGER, save, DIMENSION(:), ALLOCATABLE :: rcb2tot_idx
46	INTEGER, save, DIMENSION(:), ALLOCATABLE :: otf2rcb_idx
47
48	INTEGER, save, DIMENSION(:), ALLOCATABLE :: permut_idx
49	!$OMP THREADPRIVATE(gasi_recomb,gasv_recomb,w_cum,otf2tot_idx,rcb2tot_idx,otf2rcb_idx,permut_idx)
50
51	CONTAINS
52
53
54	SUBROUTINE ini_recombin
55
56	USE tracer_h
57
58	IMPLICIT NONE
59
60	INTEGER :: nspcrad ! # of is_rad species (tempooray here, should be a radcommon variabke)
61	INTEGER :: iq, icount
62
63	! Default values
64	use_premix = .true.
65	corrk_recombin = .false.
66	nrecomb_tot = 0
67	nrecomb_qset = 0
68	nrecomb_qotf = 0
69
70	nspcrad = 0
71
72	do iq=1,nqtot
73
74	! Counter used to check if all rad. species are recombin then no premix
75	if(is_rad(iq)==1) nspcrad = nspcrad + 1
76
77	! Sanity checks
78	if (is_recomb(iq)==1 .and. is_rad(iq)==0) then
79	write(,) 'initracer : Error for tracer iq=',iq
80	write(,) 'is_recomb=1 but is_rad=0, this is not possible !'
81	call abort_physic('initrac','Conflicting traceur.def options',1)
82	endif
83	if (is_recomb_qset(iq)==1 .and. is_recomb(iq)==0) then
84	write(,) 'initracer : Error for tracer iq=',iq
85	write(,) 'is_recomb_qset=1 but is_recomb=0, this is not possible !'
86	call abort_physic('initrac','Conflicting traceur.def options',1)
87	endif
88	if (is_recomb_qotf(iq)==1 .and. is_recomb_qset(iq)==1) then
89	write(,) 'initracer : Error for tracer iq=',iq
90	write(,) 'is_recomb_qset=1 and is_recomb_qotf=1, this is not possible !'
91	call abort_physic('initrac','Conflicting traceur.def options',1)
92	endif
93
94	if(is_recomb(iq)==1) then
95	corrk_recombin = .true. ! activating on first one found would be enough actually but I'm lazy
96
97	nrecomb_tot = nrecomb_tot + 1
98
99	if(is_recomb_qset(iq)==1) nrecomb_qset = nrecomb_qset + 1
100	if(is_recomb_qotf(iq)==1) nrecomb_qotf = nrecomb_qotf + 1
101
102	write(,)
103	write(,) 'ini_recombin: found specie : Name = ',trim(noms(iq)), &
104	' ; Predefined vmr=', is_recomb_qset(iq), &
105	' ; On-the-fly=', is_recomb_qotf(iq)
106
107	endif
108
109	enddo
110
111	! Init. correspondance array of indexes between subset of tracers
112	IF(.NOT. ALLOCATED(otf2tot_idx)) ALLOCATE(otf2tot_idx(nrecomb_qotf))
113	icount=0
114	do iq=1,nqtot
115	if(is_recomb_qotf(iq)==1) then
116	icount=icount+1
117	otf2tot_idx(icount) = iq
118	endif
119	enddo
120
121	IF(.NOT. ALLOCATED(rcb2tot_idx)) ALLOCATE(rcb2tot_idx(nrecomb_tot))
122	icount=0
123	do iq=1,nqtot
124	if(is_recomb(iq)==1) then
125	icount=icount+1
126	rcb2tot_idx(icount) = iq
127	endif
128	enddo
129
130	IF(.NOT. ALLOCATED(otf2rcb_idx)) ALLOCATE(otf2rcb_idx(nrecomb_qotf))
131	icount=0
132	do iq=1,nrecomb_tot
133	if(is_recomb_qotf(rcb2tot_idx(iq))==1) then
134	icount=icount+1
135	otf2rcb_idx(icount) = iq
136	endif
137	enddo
138
139	! Use a premix set on top ?
140	if (nspcrad == nrecomb_tot .and. nspcrad /= 0) use_premix = .false. ! In this case all rad. species are recombined
141
142	! Summary
143	write(,)
144	write(,) 'ini_recombin: Total species found for corrk recombination =', nrecomb_tot
145
146	if (corrk_recombin) then
147	if (use_premix) then
148	write(,) 'ini_recombin: .. Total radiative species matching total species for recombination...'
149	write(,) 'ini_recombin: .. Any pre-mixed set of opacities will be ignored.'
150	else
151	write(,) 'ini_recombin: .. Found less species to recombine than total radiative species..'
152	write(,) 'ini_recombin: .. Recombination will occur ontop of premix set of opacities'
153	endif
154	else
155	write(,) 'ini_recombin: .. No species found for recombination, I will use premix set only.'
156	endif
157	write(,)
158
159	END SUBROUTINE ini_recombin
160
161
162
163	SUBROUTINE su_recombin
164	USE radinc_h
165	USE radcommon_h, only: gweight, gasi, gasv
166
167	IMPLICIT NONE
168
169	INTEGER :: i, ig, jg, ind, iw, it, ip
170
171	! Allocations
172	IF(.NOT. ALLOCATED(permut_idx)) ALLOCATE(permut_idx(L_NGAUSS*L_NGAUSS))
173	IF(.NOT. ALLOCATED(w_cum)) ALLOCATE(w_cum(L_NGAUSS))
174	IF(.NOT. ALLOCATED(gasi_recomb)) ALLOCATE(gasi_recomb(L_NTREF,L_PINT,L_REFVAR,L_NSPECTI,L_NGAUSS))
175	IF(.NOT. ALLOCATED(gasv_recomb)) ALLOCATE(gasv_recomb(L_NTREF,L_PINT,L_REFVAR,L_NSPECTV,L_NGAUSS))
176	IF(.NOT. ALLOCATED(gasi_otf)) ALLOCATE(gasi_otf(L_NGAUSS,nrecomb_qotf,L_NSPECTI,L_NTREF,L_PINT))
177	IF(.NOT. ALLOCATED(gasv_otf)) ALLOCATE(gasv_otf(L_NGAUSS,nrecomb_qotf,L_NSPECTI,L_NTREF,L_PINT))
178	IF(.NOT. ALLOCATED(wtwospec)) ALLOCATE(wtwospec(L_NGAUSS*L_NGAUSS))
179
180	! Init. for recombin_corrk firstcall
181	permut_idx = (/(i, i=1,L_NGAUSS*L_NGAUSS)/)
182
183	w_cum(1)= gweight(1)
184	DO i=2,L_NGAUSS
185	w_cum(i) = w_cum(i-1)+gweight(i)
186	ENDDO
187
188	! init wtwospec once for all
189	DO ig=1,L_NGAUSS
190	DO jg=1, L_NGAUSS
191	ind = jg+(ig-1)*L_NGAUSS
192	wtwospec(ind) = gweight(ig)*gweight(jg)
193	ENDDO
194	ENDDO
195
196	! init otf correlated-k array
197	do ip=1,L_PINT
198	do it=1,L_NTREF
199	do iw=1,L_NSPECTI
200	do i=1,nrecomb_qotf
201	do ig=1,L_NGAUSS
202	gasi_otf(ig,i,iw,it,ip) = gasi(it,ip,L_REFVAR+otf2rcb_idx(i),iw,ig) ! choose only idx corresponding to otf in gasi
203	enddo
204	enddo
205	enddo
206	enddo
207	enddo
208	do ip=1,L_PINT
209	do it=1,L_NTREF
210	do iw=1,L_NSPECTV
211	do i=1,nrecomb_qotf
212	do ig=1,L_NGAUSS
213	gasv_otf(ig,i,iw,it,ip) = gasv(it,ip,L_REFVAR+otf2rcb_idx(i),iw,ig) ! choose only idx corresponding to otf in gasv
214	enddo
215	enddo
216	enddo
217	enddo
218	enddo
219
220	gasi_recomb(:,:,:,:,:) = gasi(:,:,1:L_REFVAR,:,:) ! non-zero init (=kappa_ir)
221	gasv_recomb(:,:,:,:,:) = gasv(:,:,1:L_REFVAR,:,:) ! non-zero init (=kappa_vi)
222
223	END SUBROUTINE su_recombin
224
225
226
227	SUBROUTINE call_recombin(igrid,nlayer,pq,pplay,pt,qvar,tmid,pmid)
228
229	USE comcstfi_mod, only: mugaz
230	USE radinc_h
231	USE radcommon_h
232	USE tracer_h, only: noms, mmol, is_recomb_qotf, is_recomb_qset
233
234	IMPLICIT NONE
235
236	! Inputs
237	INTEGER, INTENT(IN) :: igrid ! lon-lat grid point
238	INTEGER, INTENT(IN) :: nlayer ! Number of atmospheric layers.
239
240	REAL*8, DIMENSION(:,:), INTENT(IN) :: pq ! Tracers vertical profiles (kg/kg)
241	REAL*8, DIMENSION(nlayer), INTENT(IN) :: pplay ! Atmospheric pressure (Pa)
242	REAL*8, DIMENSION(nlayer), INTENT(IN) :: pt ! Atmospheric temperature (K)
243	REAL*8, DIMENSION(L_LEVELS), INTENT(IN) :: qvar ! Mixing ratio of variable component (mol/mol)
244	REAL*8, DIMENSION(L_LEVELS), INTENT(IN) :: tmid ! Temperature of layers, mid levels (K)
245	REAL*8, DIMENSION(L_LEVELS), INTENT(IN) :: pmid ! Pressure of layers, mid levels (mBar)
246
247	! NB : qvar is on L_LEVELS but it has been processed in callcork compared to the one in pq
248	! so we imperatively need to take this one. Note that there is no interpolation, so
249	! pq(nlayer+1-l,ivar) is broadcast in qvar(2l) qvar(2l+1)
250
251	! Local variables
252	INTEGER :: ng
253	INTEGER :: ig,l,k,nw,iq,ip,ilay,it,ix,iw
254
255	LOGICAL :: found
256
257	REAL*8 :: fact, tmin, tmax, qmin, qmax
258	REAL*8 :: LCOEF(4), WRATIO
259
260	LOGICAL,DIMENSION(:,:,:),ALLOCATABLE,save :: useptx ! Mask on which t-p-x ref grid point will be used
261
262	REAL*8, DIMENSION(:,:), ALLOCATABLE,save :: pqr ! Tracers abundances at ref pressures used for onthefly recombining (mol/mol).
263
264	LOGICAL,SAVE :: firstcall=.true.
265	!$OMP THREADPRIVATE(firstcall)
266
267	! At firstcall we precombine all what needs to be done only once (pre-mixed,forced profiles..), if needed.
268	IF (firstcall) THEN
269
270	IF(.NOT. ALLOCATED(useptx)) ALLOCATE(useptx(L_NTREF,L_PINT,L_REFVAR))
271	useptx(:,:,:) = .false.
272
273	IF(use_premix .or. (.not.use_premix .and. nrecomb_qotf/=nrecomb_tot)) THEN ! we skip this if all species are on-the-fly
274	all_otf=.false.
275	IF(.NOT. ALLOCATED(pqr)) ALLOCATE(pqr(nrecomb_tot,L_PINT))
276	! Default value for premix and for fixed species for whom vmr has been taken
277	! into account while computing high-resolution spectra
278	pqr(:,:) = 1.0
279
280	! TODO : Missing implementation here for the tracers where we read an input profile !!
281	do iq=1,nrecomb_tot
282	if (is_recomb_qset(rcb2tot_idx(iq))==0 .and. is_recomb_qotf(rcb2tot_idx(iq))==0) then
283	print*, 'Recombining tracer ', noms(rcb2tot_idx(iq)),' requires an input profile, this is not implemented yet !!'
284	call abort_physic('call_recombin','Missing implementation',1)
285	! Read pqr(:,iq)
286	endif
287	enddo
288
289	! Recombine for all T-P-Q as we do it only once for all.
290	call recombin_corrk_ini(pqr)
291	ELSE
292	all_otf=.true.
293	ENDIF
294
295	firstcall=.false.
296	IF (nrecomb_qotf==0) corrk_recombin = .false.
297	IF(ALLOCATED(pqr)) DEALLOCATE(pqr)
298	IF(.NOT. ALLOCATED(pqr)) ALLOCATE(pqr(nrecomb_qotf,L_PINT))
299
300	ENDIF ! firstcall
301
302	! NB : To have decent CPU time recombining is not done on all gridpoints and wavelenghts but we
303	! calculate a gasi/v_recomb variable on the reference corrk-k T,P,X grid (only for T,P,X values
304	! who match the atmospheric conditions) which is then processed as a standard pre-mix in
305	! optci/v routines, but updated every time tracers on the ref P grid have varied > 1%.
306
307	! Extract tracers for species which are recombined on-the-fly
308	do ip=1,L_PINT
309
310	ilay=0
311	found = .false.
312	do l=1,nlayer
313	if (pplay(l) .le. 10.0**(pfgasref(ip)+2.0)) then ! pfgasref=log(p[mbar])
314	found=.true.
315	ilay=l-1
316	exit
317	endif
318	enddo
319
320	if (.not. found) then ! set pq to top value
321	do iq=1,nrecomb_qotf
322	pqr(iq,ip) = pq(nlayer,otf2tot_idx(iq))*mugaz/mmol(otf2tot_idx(iq)) ! mol/mol
323	enddo
324	else
325	if (ilay==0) then ! set pq to bottom value
326	do iq=1,nrecomb_qotf
327	pqr(iq,ip) = pq(1,otf2tot_idx(iq))*mugaz/mmol(otf2tot_idx(iq)) ! mol/mol
328	enddo
329	else ! standard : interp pq between layers
330	fact = (10.0**(pfgasref(ip)+2.0) - pplay(ilay+1)) / (pplay(ilay) - pplay(ilay+1)) ! pfgasref=log(p[mbar])
331	do iq=1,nrecomb_qotf
332	pqr(iq,ip) = pq(ilay,otf2tot_idx(iq))*fact pq(ilay+1,otf2tot_idx(iq))**(1.0-fact)
333	pqr(iq,ip) = pqr(iq,ip)*mugaz/mmol(otf2tot_idx(iq)) ! mol/mol
334	enddo
335	endif ! if ilay==nlayer
336	endif ! if not found
337
338	enddo ! ip=1,L_PINT
339
340	! The following useptx is a trick to call recombine only for the reference T-P-X
341	! reference grid points that are useful given the temperature range (and variable specie amount) at a given altitude
342	! (cf optci/optcv routines where we interpolate corrk calling tpindex)
343	! It saves a looot of time - JVO 18
344
345	do K=2,L_LEVELS
346	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR,LCOEF,it,ip,ix,WRATIO)
347	useptx(it:it+1,ip:ip+1,ix:ix+1) = .true.
348	end do
349
350	if (.not.all_otf) then
351	call recombin_corrk_mix(igrid,pqr,useptx)
352	else
353	if (nrecomb_qotf.gt.1) then
354	call recombin_corrk_mix_allotf(igrid,pqr,useptx)
355	else
356	do ix=1,L_REFVAR
357	do ip=1,L_PINT
358	do it=1,L_NTREF
359	if (.not. useptx(it,ip,ix)) cycle
360	gasi_recomb(it,ip,ix,:,:) = pqr(1,ip)*gasi_otf(:,1,:,it,ip)
361	gasv_recomb(it,ip,ix,:,:) = pqr(1,ip)*gasv_otf(:,1,:,it,ip)
362	useptx(it,ip,ix) = .false.
363	enddo
364	enddo
365	enddo
366	endif
367	endif
368
369	END SUBROUTINE call_recombin
370
371	SUBROUTINE recombin_corrk_ini(pqr)
372
373	USE radinc_h
374	USE radcommon_h, only: gweight, tlimit, gasi, gasv
375	USE tracer_h, only: is_recomb_qotf
376	USE sort_mod, only: isort
377
378	!-----------------------------------------------------------------------
379	! Declarations:
380	! -------------
381
382	IMPLICIT NONE
383
384	! Arguments :
385	! -----------
386	REAL*8, DIMENSION(nrecomb_qotf,L_PINT), INTENT(IN) :: pqr ! otf species mixing ration
387
388	! Local variables :
389	! -----------------
390	INTEGER :: it, ip, ix, iw, ig, jg, ind, ibin, ispec
391
392	REAL*8, DIMENSION(L_NGAUSS) :: tmpk ! otf correlated-k by mixing ratio
393	REAL*8, DIMENSION(L_NGAUSS) :: krecomb
394
395	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: ktwospec
396	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: ktwospec_s
397	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: wtwospec_s
398	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: wtwospec_cum
399
400	REAL*8 :: wsplit
401
402	do ix=1,L_REFVAR
403	do ip=1,L_PINT
404	do it=1,L_NTREF
405
406	! -------------------
407	! I. INFRARED
408	! -------------------
409
410	DO iw=1,L_NSPECTI
411	DO ig=1,L_NGAUSS ! init correlated-k with premix
412	! utiliser directement gasi_recomb au lieu d'un variable intermediere krecomb ?
413	! Peut ok si L_NGAUSS première dimension de gasi_recomb
414	krecomb(ig) = gasi(it,ip,ix,iw,ig)
415	ENDDO
416	DO ispec=1,nrecomb_tot ! Loop on additional species
417
418	IF(is_recomb_qotf(rcb2tot_idx(ispec))==1) CYCLE
419	! takes all to recomb, the otf ones are skipped
420	DO ig=1,L_NGAUSS
421	tmpk(ig) = pqr(ispec,ip)*gasi(it,ip,L_REFVAR+ispec,iw,ig)
422	ENDDO
423
424	! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band
425	IF ( tmpk(L_NGAUSS-1) .LE. tlimit ) CYCLE
426	IF ( ALL( tmpk(1:L_NGAUSS-1) .LE. krecomb(1:L_NGAUSS-1)*0.1 ) ) CYCLE
427	IF ( ALL( krecomb(1:L_NGAUSS-1) .LE. tmpk(1:L_NGAUSS-1)*0.1 ) ) THEN
428	krecomb(1:L_NGAUSS-1) = tmpk(1:L_NGAUSS-1)
429	CYCLE
430	ENDIF
431
432	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
433	! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~
434	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
435	DO ig=1,L_NGAUSS
436	DO jg=1, L_NGAUSS
437	ind = jg+(ig-1)*L_NGAUSS
438	ktwospec(ind) = krecomb(ig)+tmpk(jg)
439	ENDDO
440	ENDDO
441
442	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
443	! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~
444	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
445
446	! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting
447	! NB : quite small array, quicker to permut with 2nd array than in place
448	DO ind=1,L_NGAUSS*L_NGAUSS
449	ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall
450	ENDDO
451
452	CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Insertion sort quicker because pre-sorted
453
454	! Sort w according to permutations of k.
455	! NB : quite small array, quicker to permut with 2nd array than in place
456	DO ind=1,L_NGAUSS*L_NGAUSS
457	wtwospec_s(ind) = wtwospec(permut_idx(ind))
458	ENDDO
459
460	wtwospec_cum(1) = wtwospec_s(1)
461	DO ind=2,L_NGAUSS*L_NGAUSS
462	wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind)
463	ENDDO
464
465	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
466	! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~
467	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
468	ibin=1
469
470	krecomb(:)=0.0
471
472	DO ig=1, L_NGAUSS-1
473
474	DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while
475	IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN
476	wsplit = w_cum(ig) - wtwospec_cum(ind-1)
477	krecomb(ig) = krecomb(ig) &
478	+ sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) &
479	+ wsplit*ktwospec_s(ind)
480	krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind)
481	ibin=ind+1
482	EXIT
483	ENDIF
484	ENDDO
485
486	krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) )
487
488	ENDDO
489
490	krecomb(1:L_NGAUSS-1) = krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0
491
492	ENDDO ! ispec=1,nrecomb_qotf
493	gasi(it,ip,ix,iw,:) = krecomb(:)
494	ENDDO ! iw=1,L_NSPECTI
495
496	! -------------------
497	! II. VISIBLE
498	! -------------------
499
500	DO iw=1,L_NSPECTV
501	DO ig=1,L_NGAUSS ! init correlated-k with premix
502	krecomb(ig) = gasv(it,ip,ix,iw,ig) ! there is a prerecombined cocktail
503	ENDDO
504	DO ispec=1,nrecomb_tot ! Loop on additional species
505
506	IF(is_recomb_qotf(rcb2tot_idx(ispec))==1) CYCLE
507	! takes all to recomb, the otf ones are skipped
508	DO ig=1,L_NGAUSS
509	tmpk(ig) = pqr(ispec,ip)*gasv(it,ip,L_REFVAR+ispec,iw,ig)
510	ENDDO
511
512	! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band
513	IF ( tmpk(L_NGAUSS-1) .LE. tlimit ) CYCLE
514	IF ( ALL( tmpk(1:L_NGAUSS-1) .LE. krecomb(1:L_NGAUSS-1)*0.1 ) ) CYCLE
515	IF ( ALL( krecomb(1:L_NGAUSS-1) .LE. tmpk(1:L_NGAUSS-1)*0.1 ) ) THEN
516	krecomb(1:L_NGAUSS-1) = tmpk(1:L_NGAUSS-1)
517	CYCLE
518	ENDIF
519
520	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
521	! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~
522	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
523	DO ig=1,L_NGAUSS
524	DO jg=1, L_NGAUSS
525	ind = jg+(ig-1)*L_NGAUSS
526	ktwospec(ind) = krecomb(ig)+tmpk(jg)
527	ENDDO
528	ENDDO
529
530	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
531	! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~
532	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
533
534	! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting
535	! NB : quite small array, quicker to permut with 2nd array than in place
536	DO ind=1,L_NGAUSS*L_NGAUSS
537	ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall
538	ENDDO
539
540	CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Insertion sort quicker because pre-sorted
541
542	! Sort w according to permutations of k.
543	! NB : quite small array, quicker to permut with 2nd array than in place
544	DO ind=1,L_NGAUSS*L_NGAUSS
545	wtwospec_s(ind) = wtwospec(permut_idx(ind))
546	ENDDO
547
548	wtwospec_cum(1) = wtwospec_s(1)
549	DO ind=2,L_NGAUSS*L_NGAUSS
550	wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind)
551	ENDDO
552
553	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
554	! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~
555	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
556	ibin=1
557
558	krecomb(:)=0.0
559
560	DO ig=1, L_NGAUSS-1
561
562	DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while
563	IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN
564	wsplit = w_cum(ig) - wtwospec_cum(ind-1)
565	krecomb(ig) = krecomb(ig) &
566	+ sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) &
567	+ wsplit*ktwospec_s(ind)
568	krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind)
569	ibin=ind+1
570	EXIT
571	ENDIF
572	ENDDO
573
574	krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) )
575
576	ENDDO
577
578	krecomb(1:L_NGAUSS-1) = krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0
579
580	ENDDO ! ispec=1,nrecomb_qotf
581	gasv(it,ip,ix,iw,:) = krecomb(:)
582	ENDDO ! iw=1,L_NSPECTV
583
584	enddo ! ix=1,L_NTREF
585	enddo ! it=1,L_PINT
586	enddo ! ip=1,L_REFVAR
587
588	END SUBROUTINE recombin_corrk_ini
589
590	SUBROUTINE recombin_corrk_mix(igrid,pqr,useptx)
591
592	USE radinc_h
593	USE radcommon_h, only: gweight, tlimit, gasi, gasv
594	USE sort_mod, only: isort
595	use comsaison_h, only: fract
596
597	!-----------------------------------------------------------------------
598	! Declarations:
599	! -------------
600
601	IMPLICIT NONE
602
603	! Arguments :
604	! -----------
605	INTEGER, INTENT(IN) :: igrid ! lon-lat grid point
606	REAL*8, DIMENSION(nrecomb_qotf,L_PINT), INTENT(IN) :: pqr ! otf species mixing ration
607	LOGICAL, INTENT(INOUT) :: useptx(L_NTREF,L_PINT,L_REFVAR) ! Mask on which t-p-x ref grid point will be used
608
609	! Local variables :
610	! -----------------
611	INTEGER :: it, ip, ix, iw, ig, jg, ind, ibin, ispec
612
613	REAL*8, DIMENSION(L_NGAUSS) :: tmpk ! otf correlated-k by mixing ratio
614	REAL*8, DIMENSION(L_NGAUSS) :: krecomb
615
616	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: ktwospec
617	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: ktwospec_s
618	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: wtwospec_s
619	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: wtwospec_cum
620
621	REAL*8 :: wsplit
622
623	do ix=1,L_REFVAR
624	do ip=1,L_PINT
625	do it=1,L_NTREF
626	if (.not. useptx(it,ip,ix)) cycle
627
628	! -------------------
629	! I. INFRARED
630	! -------------------
631
632	DO iw=1,L_NSPECTI
633	DO ig=1,L_NGAUSS ! init correlated-k with premix
634	! utiliser directement gasi_recomb au lieu d'un variable intermediere krecomb ?
635	! Peut ok si L_NGAUSS première dimension de gasi_recomb
636	krecomb(ig) = gasi(it,ip,ix,iw,ig)
637	ENDDO
638	DO ispec=1,nrecomb_qotf ! Loop on additional species
639
640	DO ig=1,L_NGAUSS
641	tmpk(ig) = pqr(ispec,ip)*gasi_otf(ig,ispec,iw,it,ip)
642	ENDDO
643
644	! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band
645	IF ( tmpk(L_NGAUSS-1) .LE. tlimit ) CYCLE
646	IF ( ALL( tmpk(1:L_NGAUSS-1) .LE. krecomb(1:L_NGAUSS-1)*0.1 ) ) CYCLE
647	IF ( ALL( krecomb(1:L_NGAUSS-1) .LE. tmpk(1:L_NGAUSS-1)*0.1 ) ) THEN
648	krecomb(1:L_NGAUSS-1) = tmpk(1:L_NGAUSS-1)
649	CYCLE
650	ENDIF
651
652	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
653	! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~
654	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
655	DO ig=1,L_NGAUSS
656	DO jg=1, L_NGAUSS
657	ind = jg+(ig-1)*L_NGAUSS
658	ktwospec(ind) = krecomb(ig)+tmpk(jg)
659	ENDDO
660	ENDDO
661
662	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
663	! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~
664	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
665
666	! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting
667	! NB : quite small array, quicker to permut with 2nd array than in place
668	DO ind=1,L_NGAUSS*L_NGAUSS
669	ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall
670	ENDDO
671
672	CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Insertion sort quicker because pre-sorted
673
674	! Sort w according to permutations of k.
675	! NB : quite small array, quicker to permut with 2nd array than in place
676	DO ind=1,L_NGAUSS*L_NGAUSS
677	wtwospec_s(ind) = wtwospec(permut_idx(ind))
678	ENDDO
679
680	wtwospec_cum(1) = wtwospec_s(1)
681	DO ind=2,L_NGAUSS*L_NGAUSS
682	wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind)
683	ENDDO
684
685	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
686	! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~
687	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
688	ibin=1
689
690	krecomb(:)=0.0
691
692	DO ig=1, L_NGAUSS-1
693
694	DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while
695	IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN
696	wsplit = w_cum(ig) - wtwospec_cum(ind-1)
697	krecomb(ig) = krecomb(ig) &
698	+ sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) &
699	+ wsplit*ktwospec_s(ind)
700	krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind)
701	ibin=ind+1
702	EXIT
703	ENDIF
704	ENDDO
705
706	krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) )
707
708	ENDDO
709
710	krecomb(1:L_NGAUSS-1) = krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0
711
712	ENDDO ! ispec=1,nrecomb_qotf
713	gasi_recomb(it,ip,ix,iw,:) = krecomb(:)
714	ENDDO ! iw=1,L_NSPECTI
715
716	! -------------------
717	! II. VISIBLE
718	! -------------------
719
720	if(fract(igrid) .lt. 1.0e-4) cycle ! Only during daylight.
721
722	DO iw=1,L_NSPECTV
723	DO ig=1,L_NGAUSS ! init correlated-k with premix
724	krecomb(ig) = gasv(it,ip,ix,iw,ig) ! gasv_loc order correctly for running time?
725	ENDDO
726	DO ispec=1,nrecomb_qotf ! Loop on additional species
727
728	DO ig=1,L_NGAUSS
729	tmpk(ig) = pqr(ispec,ip)*gasv_otf(ig,ispec,iw,it,ip)
730	ENDDO
731
732	! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band
733	IF ( tmpk(L_NGAUSS-1) .LE. tlimit ) CYCLE
734	IF ( ALL( tmpk(1:L_NGAUSS-1) .LE. krecomb(1:L_NGAUSS-1)*0.1 ) ) CYCLE
735	IF ( ALL( krecomb(1:L_NGAUSS-1) .LE. tmpk(1:L_NGAUSS-1)*0.1 ) ) THEN
736	krecomb(1:L_NGAUSS-1) = tmpk(1:L_NGAUSS-1)
737	CYCLE
738	ENDIF
739
740	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
741	! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~
742	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
743	DO ig=1,L_NGAUSS
744	DO jg=1, L_NGAUSS
745	ind = jg+(ig-1)*L_NGAUSS
746	ktwospec(ind) = krecomb(ig)+tmpk(jg)
747	ENDDO
748	ENDDO
749
750	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
751	! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~
752	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
753
754	! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting
755	! NB : quite small array, quicker to permut with 2nd array than in place
756	DO ind=1,L_NGAUSS*L_NGAUSS
757	ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall
758	ENDDO
759
760	CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Insertion sort quicker because pre-sorted
761
762	! Sort w according to permutations of k.
763	! NB : quite small array, quicker to permut with 2nd array than in place
764	DO ind=1,L_NGAUSS*L_NGAUSS
765	wtwospec_s(ind) = wtwospec(permut_idx(ind))
766	ENDDO
767
768	wtwospec_cum(1) = wtwospec_s(1)
769	DO ind=2,L_NGAUSS*L_NGAUSS
770	wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind)
771	ENDDO
772
773	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
774	! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~
775	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
776	ibin=1
777
778	krecomb(:)=0.0
779
780	DO ig=1, L_NGAUSS-1
781
782	DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while
783	IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN
784	wsplit = w_cum(ig) - wtwospec_cum(ind-1)
785	krecomb(ig) = krecomb(ig) &
786	+ sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) &
787	+ wsplit*ktwospec_s(ind)
788	krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind)
789	ibin=ind+1
790	EXIT
791	ENDIF
792	ENDDO
793
794	krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) )
795
796	ENDDO
797
798	krecomb(1:L_NGAUSS-1) = krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0
799
800	ENDDO ! ispec=1,nrecomb_qotf
801	gasv_recomb(it,ip,ix,iw,:) = krecomb(:)
802	ENDDO ! iw=1,L_NSPECTV
803
804	useptx(it,ip,ix) = .false.
805	enddo ! ix=1,L_NTREF
806	enddo ! it=1,L_PINT
807	enddo ! ip=1,L_REFVAR
808
809	END SUBROUTINE recombin_corrk_mix
810
811	SUBROUTINE recombin_corrk_mix_allotf(igrid,pqr,useptx)
812
813	USE radinc_h
814	USE radcommon_h, only: gweight, tlimit, gasi, gasv
815	USE sort_mod, only: isort
816	use comsaison_h, only: fract
817
818	!-----------------------------------------------------------------------
819	! Declarations:
820	! -------------
821
822	IMPLICIT NONE
823
824	! Arguments :
825	! -----------
826	INTEGER, INTENT(IN) :: igrid ! lon-lat grid point
827	REAL*8, DIMENSION(nrecomb_qotf,L_PINT), INTENT(IN) :: pqr ! otf species mixing ration
828	LOGICAL, INTENT(INOUT) :: useptx(L_NTREF,L_PINT,L_REFVAR) ! Mask on which t-p-x ref grid point will be used
829
830	! Local variables :
831	! -----------------
832	INTEGER :: it, ip, ix, iw, ig, jg, ind, ibin, ispec
833
834	REAL*8, DIMENSION(L_NGAUSS) :: tmpk ! otf correlated-k by mixing ratio
835	REAL*8, DIMENSION(L_NGAUSS) :: krecomb
836
837	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: ktwospec
838	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: ktwospec_s
839	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: wtwospec_s
840	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: wtwospec_cum
841
842	REAL*8 :: wsplit
843
844	do ix=1,L_REFVAR
845	do ip=1,L_PINT
846	do it=1,L_NTREF
847	if (.not. useptx(it,ip,ix)) cycle
848
849	! -------------------
850	! I. INFRARED
851	! -------------------
852
853	DO iw=1,L_NSPECTI
854	DO ig=1,L_NGAUSS ! init correlated-k with first gas
855	krecomb(ig) = pqr(1,ip)*gasi_otf(ig,1,iw,it,ip)
856	ENDDO
857	DO ispec=2,nrecomb_qotf ! Loop on additional species
858
859	DO ig=1,L_NGAUSS
860	tmpk(ig) = pqr(ispec,ip)*gasi_otf(ig,ispec,iw,it,ip)
861	ENDDO
862
863	! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band
864	IF ( tmpk(L_NGAUSS-1) .LE. tlimit ) CYCLE
865	IF ( ALL( tmpk(1:L_NGAUSS-1) .LE. krecomb(1:L_NGAUSS-1)*0.01 ) ) CYCLE
866	IF ( ALL( krecomb(1:L_NGAUSS-1) .LE. tmpk(1:L_NGAUSS-1)*0.01 ) ) THEN
867	krecomb(1:L_NGAUSS-1) = tmpk(1:L_NGAUSS-1)
868	CYCLE
869	ENDIF
870
871	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
872	! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~
873	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
874	DO ig=1,L_NGAUSS
875	DO jg=1, L_NGAUSS
876	ind = jg+(ig-1)*L_NGAUSS
877	ktwospec(ind) = krecomb(ig)+tmpk(jg)
878	ENDDO
879	ENDDO
880
881	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
882	! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~
883	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
884
885	! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting
886	! NB : quite small array, quicker to permut with 2nd array than in place
887	DO ind=1,L_NGAUSS*L_NGAUSS
888	ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall
889	ENDDO
890
891	CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Insertion sort quicker because pre-sorted
892
893	! Sort w according to permutations of k.
894	! NB : quite small array, quicker to permut with 2nd array than in place
895	DO ind=1,L_NGAUSS*L_NGAUSS
896	wtwospec_s(ind) = wtwospec(permut_idx(ind))
897	ENDDO
898
899	wtwospec_cum(1) = wtwospec_s(1)
900	DO ind=2,L_NGAUSS*L_NGAUSS
901	wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind)
902	ENDDO
903
904	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
905	! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~
906	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
907	ibin=1
908
909	krecomb(:)=0.0
910
911	DO ig=1, L_NGAUSS-1
912
913	DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while
914	IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN
915	wsplit = w_cum(ig) - wtwospec_cum(ind-1)
916	krecomb(ig) = krecomb(ig) &
917	+ sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) &
918	+ wsplit*ktwospec_s(ind)
919	krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind)
920	ibin=ind+1
921	EXIT
922	ENDIF
923	ENDDO
924
925	krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) )
926
927	ENDDO
928
929	krecomb(1:L_NGAUSS-1) = krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0
930
931	ENDDO ! ispec=1,nrecomb_qotf
932	gasi_recomb(it,ip,ix,iw,:) = krecomb(:)
933	ENDDO ! iw=1,L_NSPECTI
934
935	! -------------------
936	! II. VISIBLE
937	! -------------------
938
939	if(fract(igrid) .lt. 1.0e-4) cycle ! Only during daylight.
940
941	DO iw=1,L_NSPECTV
942	DO ig=1,L_NGAUSS ! init correlated-k with first gas
943	krecomb(ig) = pqr(1,ip)*gasv_otf(ig,1,iw,it,ip)
944	ENDDO
945	DO ispec=2,nrecomb_qotf ! Loop on additional species
946
947	DO ig=1,L_NGAUSS
948	tmpk(ig) = pqr(ispec,ip)*gasv_otf(ig,ispec,iw,it,ip)
949	ENDDO
950
951	! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band
952	IF ( tmpk(L_NGAUSS-1) .LE. tlimit ) CYCLE
953	IF ( ALL( tmpk(1:L_NGAUSS-1) .LE. krecomb(1:L_NGAUSS-1)*0.01 ) ) CYCLE
954	IF ( ALL( krecomb(1:L_NGAUSS-1) .LE. tmpk(1:L_NGAUSS-1)*0.01 ) ) THEN
955	krecomb(1:L_NGAUSS-1) = tmpk(1:L_NGAUSS-1)
956	CYCLE
957	ENDIF
958
959	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
960	! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~
961	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
962	DO ig=1,L_NGAUSS
963	DO jg=1, L_NGAUSS
964	ind = jg+(ig-1)*L_NGAUSS
965	ktwospec(ind) = krecomb(ig)+tmpk(jg)
966	ENDDO
967	ENDDO
968
969	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
970	! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~
971	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
972
973	! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting
974	! NB : quite small array, quicker to permut with 2nd array than in place
975	DO ind=1,L_NGAUSS*L_NGAUSS
976	ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall
977	ENDDO
978
979	CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Insertion sort quicker because pre-sorted
980
981	! Sort w according to permutations of k.
982	! NB : quite small array, quicker to permut with 2nd array than in place
983	DO ind=1,L_NGAUSS*L_NGAUSS
984	wtwospec_s(ind) = wtwospec(permut_idx(ind))
985	ENDDO
986
987	wtwospec_cum(1) = wtwospec_s(1)
988	DO ind=2,L_NGAUSS*L_NGAUSS
989	wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind)
990	ENDDO
991
992	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
993	! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~
994	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
995	ibin=1
996
997	krecomb(:)=0.0
998
999	DO ig=1, L_NGAUSS-1
1000
1001	DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while
1002	IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN
1003	wsplit = w_cum(ig) - wtwospec_cum(ind-1)
1004	krecomb(ig) = krecomb(ig) &
1005	+ sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) &
1006	+ wsplit*ktwospec_s(ind)
1007	krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind)
1008	ibin=ind+1
1009	EXIT
1010	ENDIF
1011	ENDDO
1012
1013	krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) )
1014
1015	ENDDO
1016
1017	krecomb(1:L_NGAUSS-1) = krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0
1018
1019	ENDDO ! ispec=1,nrecomb_qotf
1020	gasv_recomb(it,ip,ix,iw,:) = krecomb(:)
1021	ENDDO ! iw=1,L_NSPECTV
1022
1023	useptx(it,ip,ix) = .false.
1024	enddo ! ix=1,L_NTREF
1025	enddo ! it=1,L_PINT
1026	enddo ! ip=1,L_REFVAR
1027
1028	END SUBROUTINE recombin_corrk_mix_allotf
1029
1030	END MODULE recombin_corrk_mod
1031

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

Download in other formats:

Original Format