source: trunk/LMDZ.GENERIC/libf/phystd/recombin_corrk_mod.F90 @ 2582

Last change on this file since 2582 was 2582, checked in by emillour, 3 years ago

Generic GCM:
Fixed bug in tpindex (for low temperatures, between first and second
reference temperatures, temperature was wrongly set to tref(1)).
The input temperature was also allowed to be modified by the routine,
which is probably not a good idea and no longer the case.
Took this opportunity to turn tpindex into a module.
GM+EM
File size: 45.1 KB
Line 
1MODULE 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(2*l) qvar(2*l+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      REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS)    :: ktwospec
397      REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS)    :: ktwospec_s
398      REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS)    :: wtwospec_s
399      REAL*8, DIMENSION(L_NGAUSS*L_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      REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS)    :: ktwospec
618      REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS)    :: ktwospec_s
619      REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS)    :: wtwospec_s
620      REAL*8, DIMENSION(L_NGAUSS*L_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      REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS)    :: ktwospec
839      REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS)    :: ktwospec_s
840      REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS)    :: wtwospec_s
841      REAL*8, DIMENSION(L_NGAUSS*L_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.