Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

recombin_corrk_mod.F90 @ 3558

Last change on this file since 3558 was 3184, checked in by afalco, 12 months ago
Pluto PCM: Added LMDZ.PLUTO, a copy of the generic model, cleaned from some unnecessary modules (water, ...) AF
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	USE tpindex_mod, only: tpindex
234
235	IMPLICIT NONE
236
237	! Inputs
238	INTEGER, INTENT(IN) :: igrid ! lon-lat grid point
239	INTEGER, INTENT(IN) :: nlayer ! Number of atmospheric layers.
240
241	REAL*8, DIMENSION(:,:), INTENT(IN) :: pq ! Tracers vertical profiles (kg/kg)
242	REAL*8, DIMENSION(nlayer), INTENT(IN) :: pplay ! Atmospheric pressure (Pa)
243	REAL*8, DIMENSION(nlayer), INTENT(IN) :: pt ! Atmospheric temperature (K)
244	REAL*8, DIMENSION(L_LEVELS), INTENT(IN) :: qvar ! Mixing ratio of variable component (mol/mol)
245	REAL*8, DIMENSION(L_LEVELS), INTENT(IN) :: tmid ! Temperature of layers, mid levels (K)
246	REAL*8, DIMENSION(L_LEVELS), INTENT(IN) :: pmid ! Pressure of layers, mid levels (mBar)
247
248	! NB : qvar is on L_LEVELS but it has been processed in callcork compared to the one in pq
249	! so we imperatively need to take this one. Note that there is no interpolation, so
250	! pq(nlayer+1-l,ivar) is broadcast in qvar(2l) qvar(2l+1)
251
252	! Local variables
253	INTEGER :: ng
254	INTEGER :: ig,l,k,nw,iq,ip,ilay,it,ix,iw
255
256	LOGICAL :: found
257
258	REAL*8 :: fact, tmin, tmax, qmin, qmax
259	REAL*8 :: LCOEF(4), WRATIO
260
261	LOGICAL,DIMENSION(:,:,:),ALLOCATABLE,save :: useptx ! Mask on which t-p-x ref grid point will be used
262
263	REAL*8, DIMENSION(:,:), ALLOCATABLE,save :: pqr ! Tracers abundances at ref pressures used for onthefly recombining (mol/mol).
264
265	LOGICAL,SAVE :: firstcall=.true.
266	!$OMP THREADPRIVATE(firstcall)
267
268	! At firstcall we precombine all what needs to be done only once (pre-mixed,forced profiles..), if needed.
269	IF (firstcall) THEN
270
271	IF(.NOT. ALLOCATED(useptx)) ALLOCATE(useptx(L_NTREF,L_PINT,L_REFVAR))
272	useptx(:,:,:) = .false.
273
274	IF(use_premix .or. (.not.use_premix .and. nrecomb_qotf/=nrecomb_tot)) THEN ! we skip this if all species are on-the-fly
275	all_otf=.false.
276	IF(.NOT. ALLOCATED(pqr)) ALLOCATE(pqr(nrecomb_tot,L_PINT))
277	! Default value for premix and for fixed species for whom vmr has been taken
278	! into account while computing high-resolution spectra
279	pqr(:,:) = 1.0
280
281	! TODO : Missing implementation here for the tracers where we read an input profile !!
282	do iq=1,nrecomb_tot
283	if (is_recomb_qset(rcb2tot_idx(iq))==0 .and. is_recomb_qotf(rcb2tot_idx(iq))==0) then
284	print*, 'Recombining tracer ', noms(rcb2tot_idx(iq)),' requires an input profile, this is not implemented yet !!'
285	call abort_physic('call_recombin','Missing implementation',1)
286	! Read pqr(:,iq)
287	endif
288	enddo
289
290	! Recombine for all T-P-Q as we do it only once for all.
291	call recombin_corrk_ini(pqr)
292	ELSE
293	all_otf=.true.
294	ENDIF
295
296	firstcall=.false.
297	IF (nrecomb_qotf==0) corrk_recombin = .false.
298	IF(ALLOCATED(pqr)) DEALLOCATE(pqr)
299	IF(.NOT. ALLOCATED(pqr)) ALLOCATE(pqr(nrecomb_qotf,L_PINT))
300
301	ENDIF ! firstcall
302
303	! NB : To have decent CPU time recombining is not done on all gridpoints and wavelenghts but we
304	! calculate a gasi/v_recomb variable on the reference corrk-k T,P,X grid (only for T,P,X values
305	! who match the atmospheric conditions) which is then processed as a standard pre-mix in
306	! optci/v routines, but updated every time tracers on the ref P grid have varied > 1%.
307
308	! Extract tracers for species which are recombined on-the-fly
309	do ip=1,L_PINT
310
311	ilay=0
312	found = .false.
313	do l=1,nlayer
314	if (pplay(l) .le. 10.0**(pfgasref(ip)+2.0)) then ! pfgasref=log(p[mbar])
315	found=.true.
316	ilay=l-1
317	exit
318	endif
319	enddo
320
321	if (.not. found) then ! set pq to top value
322	do iq=1,nrecomb_qotf
323	pqr(iq,ip) = pq(nlayer,otf2tot_idx(iq))*mugaz/mmol(otf2tot_idx(iq)) ! mol/mol
324	enddo
325	else
326	if (ilay==0) then ! set pq to bottom value
327	do iq=1,nrecomb_qotf
328	pqr(iq,ip) = pq(1,otf2tot_idx(iq))*mugaz/mmol(otf2tot_idx(iq)) ! mol/mol
329	enddo
330	else ! standard : interp pq between layers
331	fact = (10.0**(pfgasref(ip)+2.0) - pplay(ilay+1)) / (pplay(ilay) - pplay(ilay+1)) ! pfgasref=log(p[mbar])
332	do iq=1,nrecomb_qotf
333	pqr(iq,ip) = pq(ilay,otf2tot_idx(iq))*fact pq(ilay+1,otf2tot_idx(iq))**(1.0-fact)
334	pqr(iq,ip) = pqr(iq,ip)*mugaz/mmol(otf2tot_idx(iq)) ! mol/mol
335	enddo
336	endif ! if ilay==nlayer
337	endif ! if not found
338
339	enddo ! ip=1,L_PINT
340
341	! The following useptx is a trick to call recombine only for the reference T-P-X
342	! reference grid points that are useful given the temperature range (and variable specie amount) at a given altitude
343	! (cf optci/optcv routines where we interpolate corrk calling tpindex)
344	! It saves a looot of time - JVO 18
345
346	do K=2,L_LEVELS
347	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR,LCOEF,it,ip,ix,WRATIO)
348	useptx(it:it+1,ip:ip+1,ix:ix+1) = .true.
349	end do
350
351	if (.not.all_otf) then
352	call recombin_corrk_mix(igrid,pqr,useptx)
353	else
354	if (nrecomb_qotf.gt.1) then
355	call recombin_corrk_mix_allotf(igrid,pqr,useptx)
356	else
357	do ix=1,L_REFVAR
358	do ip=1,L_PINT
359	do it=1,L_NTREF
360	if (.not. useptx(it,ip,ix)) cycle
361	gasi_recomb(it,ip,ix,:,:) = pqr(1,ip)*gasi_otf(:,1,:,it,ip)
362	gasv_recomb(it,ip,ix,:,:) = pqr(1,ip)*gasv_otf(:,1,:,it,ip)
363	useptx(it,ip,ix) = .false.
364	enddo
365	enddo
366	enddo
367	endif
368	endif
369
370	END SUBROUTINE call_recombin
371
372	SUBROUTINE recombin_corrk_ini(pqr)
373
374	USE radinc_h
375	USE radcommon_h, only: gweight, tlimit, gasi, gasv
376	USE tracer_h, only: is_recomb_qotf
377	USE sort_mod, only: isort
378
379	!-----------------------------------------------------------------------
380	! Declarations:
381	! -------------
382
383	IMPLICIT NONE
384
385	! Arguments :
386	! -----------
387	REAL*8, DIMENSION(nrecomb_qotf,L_PINT), INTENT(IN) :: pqr ! otf species mixing ration
388
389	! Local variables :
390	! -----------------
391	INTEGER :: it, ip, ix, iw, ig, jg, ind, ibin, ispec
392
393	REAL*8, DIMENSION(L_NGAUSS) :: tmpk ! otf correlated-k by mixing ratio
394	REAL*8, DIMENSION(L_NGAUSS) :: krecomb
395
396	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: ktwospec
397	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: ktwospec_s
398	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: wtwospec_s
399	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: wtwospec_cum
400
401	REAL*8 :: wsplit
402
403	do ix=1,L_REFVAR
404	do ip=1,L_PINT
405	do it=1,L_NTREF
406
407	! -------------------
408	! I. INFRARED
409	! -------------------
410
411	DO iw=1,L_NSPECTI
412	DO ig=1,L_NGAUSS ! init correlated-k with premix
413	! utiliser directement gasi_recomb au lieu d'un variable intermediere krecomb ?
414	! Peut ok si L_NGAUSS première dimension de gasi_recomb
415	krecomb(ig) = gasi(it,ip,ix,iw,ig)
416	ENDDO
417	DO ispec=1,nrecomb_tot ! Loop on additional species
418
419	IF(is_recomb_qotf(rcb2tot_idx(ispec))==1) CYCLE
420	! takes all to recomb, the otf ones are skipped
421	DO ig=1,L_NGAUSS
422	tmpk(ig) = pqr(ispec,ip)*gasi(it,ip,L_REFVAR+ispec,iw,ig)
423	ENDDO
424
425	! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band
426	IF ( tmpk(L_NGAUSS-1) .LE. tlimit ) CYCLE
427	IF ( ALL( tmpk(1:L_NGAUSS-1) .LE. krecomb(1:L_NGAUSS-1)*0.1 ) ) CYCLE
428	IF ( ALL( krecomb(1:L_NGAUSS-1) .LE. tmpk(1:L_NGAUSS-1)*0.1 ) ) THEN
429	krecomb(1:L_NGAUSS-1) = tmpk(1:L_NGAUSS-1)
430	CYCLE
431	ENDIF
432
433	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
434	! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~
435	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
436	DO ig=1,L_NGAUSS
437	DO jg=1, L_NGAUSS
438	ind = jg+(ig-1)*L_NGAUSS
439	ktwospec(ind) = krecomb(ig)+tmpk(jg)
440	ENDDO
441	ENDDO
442
443	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
444	! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~
445	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
446
447	! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting
448	! NB : quite small array, quicker to permut with 2nd array than in place
449	DO ind=1,L_NGAUSS*L_NGAUSS
450	ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall
451	ENDDO
452
453	CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Insertion sort quicker because pre-sorted
454
455	! Sort w according to permutations of k.
456	! NB : quite small array, quicker to permut with 2nd array than in place
457	DO ind=1,L_NGAUSS*L_NGAUSS
458	wtwospec_s(ind) = wtwospec(permut_idx(ind))
459	ENDDO
460
461	wtwospec_cum(1) = wtwospec_s(1)
462	DO ind=2,L_NGAUSS*L_NGAUSS
463	wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind)
464	ENDDO
465
466	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
467	! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~
468	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
469	ibin=1
470
471	krecomb(:)=0.0
472
473	DO ig=1, L_NGAUSS-1
474
475	DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while
476	IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN
477	wsplit = w_cum(ig) - wtwospec_cum(ind-1)
478	krecomb(ig) = krecomb(ig) &
479	+ sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) &
480	+ wsplit*ktwospec_s(ind)
481	krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind)
482	ibin=ind+1
483	EXIT
484	ENDIF
485	ENDDO
486
487	krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) )
488
489	ENDDO
490
491	krecomb(1:L_NGAUSS-1) = krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0
492
493	ENDDO ! ispec=1,nrecomb_qotf
494	gasi(it,ip,ix,iw,:) = krecomb(:)
495	ENDDO ! iw=1,L_NSPECTI
496
497	! -------------------
498	! II. VISIBLE
499	! -------------------
500
501	DO iw=1,L_NSPECTV
502	DO ig=1,L_NGAUSS ! init correlated-k with premix
503	krecomb(ig) = gasv(it,ip,ix,iw,ig) ! there is a prerecombined cocktail
504	ENDDO
505	DO ispec=1,nrecomb_tot ! Loop on additional species
506
507	IF(is_recomb_qotf(rcb2tot_idx(ispec))==1) CYCLE
508	! takes all to recomb, the otf ones are skipped
509	DO ig=1,L_NGAUSS
510	tmpk(ig) = pqr(ispec,ip)*gasv(it,ip,L_REFVAR+ispec,iw,ig)
511	ENDDO
512
513	! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band
514	IF ( tmpk(L_NGAUSS-1) .LE. tlimit ) CYCLE
515	IF ( ALL( tmpk(1:L_NGAUSS-1) .LE. krecomb(1:L_NGAUSS-1)*0.1 ) ) CYCLE
516	IF ( ALL( krecomb(1:L_NGAUSS-1) .LE. tmpk(1:L_NGAUSS-1)*0.1 ) ) THEN
517	krecomb(1:L_NGAUSS-1) = tmpk(1:L_NGAUSS-1)
518	CYCLE
519	ENDIF
520
521	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
522	! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~
523	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
524	DO ig=1,L_NGAUSS
525	DO jg=1, L_NGAUSS
526	ind = jg+(ig-1)*L_NGAUSS
527	ktwospec(ind) = krecomb(ig)+tmpk(jg)
528	ENDDO
529	ENDDO
530
531	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
532	! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~
533	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
534
535	! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting
536	! NB : quite small array, quicker to permut with 2nd array than in place
537	DO ind=1,L_NGAUSS*L_NGAUSS
538	ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall
539	ENDDO
540
541	CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Insertion sort quicker because pre-sorted
542
543	! Sort w according to permutations of k.
544	! NB : quite small array, quicker to permut with 2nd array than in place
545	DO ind=1,L_NGAUSS*L_NGAUSS
546	wtwospec_s(ind) = wtwospec(permut_idx(ind))
547	ENDDO
548
549	wtwospec_cum(1) = wtwospec_s(1)
550	DO ind=2,L_NGAUSS*L_NGAUSS
551	wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind)
552	ENDDO
553
554	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
555	! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~
556	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
557	ibin=1
558
559	krecomb(:)=0.0
560
561	DO ig=1, L_NGAUSS-1
562
563	DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while
564	IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN
565	wsplit = w_cum(ig) - wtwospec_cum(ind-1)
566	krecomb(ig) = krecomb(ig) &
567	+ sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) &
568	+ wsplit*ktwospec_s(ind)
569	krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind)
570	ibin=ind+1
571	EXIT
572	ENDIF
573	ENDDO
574
575	krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) )
576
577	ENDDO
578
579	krecomb(1:L_NGAUSS-1) = krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0
580
581	ENDDO ! ispec=1,nrecomb_qotf
582	gasv(it,ip,ix,iw,:) = krecomb(:)
583	ENDDO ! iw=1,L_NSPECTV
584
585	enddo ! ix=1,L_NTREF
586	enddo ! it=1,L_PINT
587	enddo ! ip=1,L_REFVAR
588
589	END SUBROUTINE recombin_corrk_ini
590
591	SUBROUTINE recombin_corrk_mix(igrid,pqr,useptx)
592
593	USE radinc_h
594	USE radcommon_h, only: gweight, tlimit, gasi, gasv
595	USE sort_mod, only: isort
596	use comsaison_h, only: fract
597
598	!-----------------------------------------------------------------------
599	! Declarations:
600	! -------------
601
602	IMPLICIT NONE
603
604	! Arguments :
605	! -----------
606	INTEGER, INTENT(IN) :: igrid ! lon-lat grid point
607	REAL*8, DIMENSION(nrecomb_qotf,L_PINT), INTENT(IN) :: pqr ! otf species mixing ration
608	LOGICAL, INTENT(INOUT) :: useptx(L_NTREF,L_PINT,L_REFVAR) ! Mask on which t-p-x ref grid point will be used
609
610	! Local variables :
611	! -----------------
612	INTEGER :: it, ip, ix, iw, ig, jg, ind, ibin, ispec
613
614	REAL*8, DIMENSION(L_NGAUSS) :: tmpk ! otf correlated-k by mixing ratio
615	REAL*8, DIMENSION(L_NGAUSS) :: krecomb
616
617	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: ktwospec
618	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: ktwospec_s
619	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: wtwospec_s
620	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: wtwospec_cum
621
622	REAL*8 :: wsplit
623
624	do ix=1,L_REFVAR
625	do ip=1,L_PINT
626	do it=1,L_NTREF
627	if (.not. useptx(it,ip,ix)) cycle
628
629	! -------------------
630	! I. INFRARED
631	! -------------------
632
633	DO iw=1,L_NSPECTI
634	DO ig=1,L_NGAUSS ! init correlated-k with premix
635	! utiliser directement gasi_recomb au lieu d'un variable intermediere krecomb ?
636	! Peut ok si L_NGAUSS première dimension de gasi_recomb
637	krecomb(ig) = gasi(it,ip,ix,iw,ig)
638	ENDDO
639	DO ispec=1,nrecomb_qotf ! Loop on additional species
640
641	DO ig=1,L_NGAUSS
642	tmpk(ig) = pqr(ispec,ip)*gasi_otf(ig,ispec,iw,it,ip)
643	ENDDO
644
645	! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band
646	IF ( tmpk(L_NGAUSS-1) .LE. tlimit ) CYCLE
647	IF ( ALL( tmpk(1:L_NGAUSS-1) .LE. krecomb(1:L_NGAUSS-1)*0.1 ) ) CYCLE
648	IF ( ALL( krecomb(1:L_NGAUSS-1) .LE. tmpk(1:L_NGAUSS-1)*0.1 ) ) THEN
649	krecomb(1:L_NGAUSS-1) = tmpk(1:L_NGAUSS-1)
650	CYCLE
651	ENDIF
652
653	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
654	! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~
655	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
656	DO ig=1,L_NGAUSS
657	DO jg=1, L_NGAUSS
658	ind = jg+(ig-1)*L_NGAUSS
659	ktwospec(ind) = krecomb(ig)+tmpk(jg)
660	ENDDO
661	ENDDO
662
663	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
664	! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~
665	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
666
667	! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting
668	! NB : quite small array, quicker to permut with 2nd array than in place
669	DO ind=1,L_NGAUSS*L_NGAUSS
670	ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall
671	ENDDO
672
673	CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Insertion sort quicker because pre-sorted
674
675	! Sort w according to permutations of k.
676	! NB : quite small array, quicker to permut with 2nd array than in place
677	DO ind=1,L_NGAUSS*L_NGAUSS
678	wtwospec_s(ind) = wtwospec(permut_idx(ind))
679	ENDDO
680
681	wtwospec_cum(1) = wtwospec_s(1)
682	DO ind=2,L_NGAUSS*L_NGAUSS
683	wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind)
684	ENDDO
685
686	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
687	! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~
688	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
689	ibin=1
690
691	krecomb(:)=0.0
692
693	DO ig=1, L_NGAUSS-1
694
695	DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while
696	IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN
697	wsplit = w_cum(ig) - wtwospec_cum(ind-1)
698	krecomb(ig) = krecomb(ig) &
699	+ sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) &
700	+ wsplit*ktwospec_s(ind)
701	krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind)
702	ibin=ind+1
703	EXIT
704	ENDIF
705	ENDDO
706
707	krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) )
708
709	ENDDO
710
711	krecomb(1:L_NGAUSS-1) = krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0
712
713	ENDDO ! ispec=1,nrecomb_qotf
714	gasi_recomb(it,ip,ix,iw,:) = krecomb(:)
715	ENDDO ! iw=1,L_NSPECTI
716
717	! -------------------
718	! II. VISIBLE
719	! -------------------
720
721	if(fract(igrid) .lt. 1.0e-4) cycle ! Only during daylight.
722
723	DO iw=1,L_NSPECTV
724	DO ig=1,L_NGAUSS ! init correlated-k with premix
725	krecomb(ig) = gasv(it,ip,ix,iw,ig) ! gasv_loc order correctly for running time?
726	ENDDO
727	DO ispec=1,nrecomb_qotf ! Loop on additional species
728
729	DO ig=1,L_NGAUSS
730	tmpk(ig) = pqr(ispec,ip)*gasv_otf(ig,ispec,iw,it,ip)
731	ENDDO
732
733	! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band
734	IF ( tmpk(L_NGAUSS-1) .LE. tlimit ) CYCLE
735	IF ( ALL( tmpk(1:L_NGAUSS-1) .LE. krecomb(1:L_NGAUSS-1)*0.1 ) ) CYCLE
736	IF ( ALL( krecomb(1:L_NGAUSS-1) .LE. tmpk(1:L_NGAUSS-1)*0.1 ) ) THEN
737	krecomb(1:L_NGAUSS-1) = tmpk(1:L_NGAUSS-1)
738	CYCLE
739	ENDIF
740
741	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
742	! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~
743	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
744	DO ig=1,L_NGAUSS
745	DO jg=1, L_NGAUSS
746	ind = jg+(ig-1)*L_NGAUSS
747	ktwospec(ind) = krecomb(ig)+tmpk(jg)
748	ENDDO
749	ENDDO
750
751	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
752	! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~
753	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
754
755	! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting
756	! NB : quite small array, quicker to permut with 2nd array than in place
757	DO ind=1,L_NGAUSS*L_NGAUSS
758	ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall
759	ENDDO
760
761	CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Insertion sort quicker because pre-sorted
762
763	! Sort w according to permutations of k.
764	! NB : quite small array, quicker to permut with 2nd array than in place
765	DO ind=1,L_NGAUSS*L_NGAUSS
766	wtwospec_s(ind) = wtwospec(permut_idx(ind))
767	ENDDO
768
769	wtwospec_cum(1) = wtwospec_s(1)
770	DO ind=2,L_NGAUSS*L_NGAUSS
771	wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind)
772	ENDDO
773
774	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
775	! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~
776	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
777	ibin=1
778
779	krecomb(:)=0.0
780
781	DO ig=1, L_NGAUSS-1
782
783	DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while
784	IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN
785	wsplit = w_cum(ig) - wtwospec_cum(ind-1)
786	krecomb(ig) = krecomb(ig) &
787	+ sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) &
788	+ wsplit*ktwospec_s(ind)
789	krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind)
790	ibin=ind+1
791	EXIT
792	ENDIF
793	ENDDO
794
795	krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) )
796
797	ENDDO
798
799	krecomb(1:L_NGAUSS-1) = krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0
800
801	ENDDO ! ispec=1,nrecomb_qotf
802	gasv_recomb(it,ip,ix,iw,:) = krecomb(:)
803	ENDDO ! iw=1,L_NSPECTV
804
805	useptx(it,ip,ix) = .false.
806	enddo ! ix=1,L_NTREF
807	enddo ! it=1,L_PINT
808	enddo ! ip=1,L_REFVAR
809
810	END SUBROUTINE recombin_corrk_mix
811
812	SUBROUTINE recombin_corrk_mix_allotf(igrid,pqr,useptx)
813
814	USE radinc_h
815	USE radcommon_h, only: gweight, tlimit, gasi, gasv
816	USE sort_mod, only: isort
817	use comsaison_h, only: fract
818
819	!-----------------------------------------------------------------------
820	! Declarations:
821	! -------------
822
823	IMPLICIT NONE
824
825	! Arguments :
826	! -----------
827	INTEGER, INTENT(IN) :: igrid ! lon-lat grid point
828	REAL*8, DIMENSION(nrecomb_qotf,L_PINT), INTENT(IN) :: pqr ! otf species mixing ration
829	LOGICAL, INTENT(INOUT) :: useptx(L_NTREF,L_PINT,L_REFVAR) ! Mask on which t-p-x ref grid point will be used
830
831	! Local variables :
832	! -----------------
833	INTEGER :: it, ip, ix, iw, ig, jg, ind, ibin, ispec
834
835	REAL*8, DIMENSION(L_NGAUSS) :: tmpk ! otf correlated-k by mixing ratio
836	REAL*8, DIMENSION(L_NGAUSS) :: krecomb
837
838	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: ktwospec
839	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: ktwospec_s
840	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: wtwospec_s
841	REAL8, DIMENSION(L_NGAUSSL_NGAUSS) :: wtwospec_cum
842
843	REAL*8 :: wsplit
844
845	do ix=1,L_REFVAR
846	do ip=1,L_PINT
847	do it=1,L_NTREF
848	if (.not. useptx(it,ip,ix)) cycle
849
850	! -------------------
851	! I. INFRARED
852	! -------------------
853
854	DO iw=1,L_NSPECTI
855	DO ig=1,L_NGAUSS ! init correlated-k with first gas
856	krecomb(ig) = pqr(1,ip)*gasi_otf(ig,1,iw,it,ip)
857	ENDDO
858	DO ispec=2,nrecomb_qotf ! Loop on additional species
859
860	DO ig=1,L_NGAUSS
861	tmpk(ig) = pqr(ispec,ip)*gasi_otf(ig,ispec,iw,it,ip)
862	ENDDO
863
864	! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band
865	IF ( tmpk(L_NGAUSS-1) .LE. tlimit ) CYCLE
866	IF ( ALL( tmpk(1:L_NGAUSS-1) .LE. krecomb(1:L_NGAUSS-1)*0.01 ) ) CYCLE
867	IF ( ALL( krecomb(1:L_NGAUSS-1) .LE. tmpk(1:L_NGAUSS-1)*0.01 ) ) THEN
868	krecomb(1:L_NGAUSS-1) = tmpk(1:L_NGAUSS-1)
869	CYCLE
870	ENDIF
871
872	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
873	! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~
874	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
875	DO ig=1,L_NGAUSS
876	DO jg=1, L_NGAUSS
877	ind = jg+(ig-1)*L_NGAUSS
878	ktwospec(ind) = krecomb(ig)+tmpk(jg)
879	ENDDO
880	ENDDO
881
882	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
883	! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~
884	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
885
886	! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting
887	! NB : quite small array, quicker to permut with 2nd array than in place
888	DO ind=1,L_NGAUSS*L_NGAUSS
889	ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall
890	ENDDO
891
892	CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Insertion sort quicker because pre-sorted
893
894	! Sort w according to permutations of k.
895	! NB : quite small array, quicker to permut with 2nd array than in place
896	DO ind=1,L_NGAUSS*L_NGAUSS
897	wtwospec_s(ind) = wtwospec(permut_idx(ind))
898	ENDDO
899
900	wtwospec_cum(1) = wtwospec_s(1)
901	DO ind=2,L_NGAUSS*L_NGAUSS
902	wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind)
903	ENDDO
904
905	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
906	! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~
907	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
908	ibin=1
909
910	krecomb(:)=0.0
911
912	DO ig=1, L_NGAUSS-1
913
914	DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while
915	IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN
916	wsplit = w_cum(ig) - wtwospec_cum(ind-1)
917	krecomb(ig) = krecomb(ig) &
918	+ sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) &
919	+ wsplit*ktwospec_s(ind)
920	krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind)
921	ibin=ind+1
922	EXIT
923	ENDIF
924	ENDDO
925
926	krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) )
927
928	ENDDO
929
930	krecomb(1:L_NGAUSS-1) = krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0
931
932	ENDDO ! ispec=1,nrecomb_qotf
933	gasi_recomb(it,ip,ix,iw,:) = krecomb(:)
934	ENDDO ! iw=1,L_NSPECTI
935
936	! -------------------
937	! II. VISIBLE
938	! -------------------
939
940	if(fract(igrid) .lt. 1.0e-4) cycle ! Only during daylight.
941
942	DO iw=1,L_NSPECTV
943	DO ig=1,L_NGAUSS ! init correlated-k with first gas
944	krecomb(ig) = pqr(1,ip)*gasv_otf(ig,1,iw,it,ip)
945	ENDDO
946	DO ispec=2,nrecomb_qotf ! Loop on additional species
947
948	DO ig=1,L_NGAUSS
949	tmpk(ig) = pqr(ispec,ip)*gasv_otf(ig,ispec,iw,it,ip)
950	ENDDO
951
952	! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band
953	IF ( tmpk(L_NGAUSS-1) .LE. tlimit ) CYCLE
954	IF ( ALL( tmpk(1:L_NGAUSS-1) .LE. krecomb(1:L_NGAUSS-1)*0.01 ) ) CYCLE
955	IF ( ALL( krecomb(1:L_NGAUSS-1) .LE. tmpk(1:L_NGAUSS-1)*0.01 ) ) THEN
956	krecomb(1:L_NGAUSS-1) = tmpk(1:L_NGAUSS-1)
957	CYCLE
958	ENDIF
959
960	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
961	! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~
962	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
963	DO ig=1,L_NGAUSS
964	DO jg=1, L_NGAUSS
965	ind = jg+(ig-1)*L_NGAUSS
966	ktwospec(ind) = krecomb(ig)+tmpk(jg)
967	ENDDO
968	ENDDO
969
970	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
971	! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~
972	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
973
974	! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting
975	! NB : quite small array, quicker to permut with 2nd array than in place
976	DO ind=1,L_NGAUSS*L_NGAUSS
977	ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall
978	ENDDO
979
980	CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Insertion sort quicker because pre-sorted
981
982	! Sort w according to permutations of k.
983	! NB : quite small array, quicker to permut with 2nd array than in place
984	DO ind=1,L_NGAUSS*L_NGAUSS
985	wtwospec_s(ind) = wtwospec(permut_idx(ind))
986	ENDDO
987
988	wtwospec_cum(1) = wtwospec_s(1)
989	DO ind=2,L_NGAUSS*L_NGAUSS
990	wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind)
991	ENDDO
992
993	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
994	! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~
995	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
996	ibin=1
997
998	krecomb(:)=0.0
999
1000	DO ig=1, L_NGAUSS-1
1001
1002	DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while
1003	IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN
1004	wsplit = w_cum(ig) - wtwospec_cum(ind-1)
1005	krecomb(ig) = krecomb(ig) &
1006	+ sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) &
1007	+ wsplit*ktwospec_s(ind)
1008	krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind)
1009	ibin=ind+1
1010	EXIT
1011	ENDIF
1012	ENDDO
1013
1014	krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) )
1015
1016	ENDDO
1017
1018	krecomb(1:L_NGAUSS-1) = krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0
1019
1020	ENDDO ! ispec=1,nrecomb_qotf
1021	gasv_recomb(it,ip,ix,iw,:) = krecomb(:)
1022	ENDDO ! iw=1,L_NSPECTV
1023
1024	useptx(it,ip,ix) = .false.
1025	enddo ! ix=1,L_NTREF
1026	enddo ! it=1,L_PINT
1027	enddo ! ip=1,L_REFVAR
1028
1029	END SUBROUTINE recombin_corrk_mix_allotf
1030
1031	END MODULE recombin_corrk_mod
1032

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

Download in other formats:

Original Format