MODULE recombin_corrk_mod ! ! Author : Jan Vatant d'Ollone (2018-2020) ! ! This module contains the following subroutines : ! - ini_recombin : From modern traceur.def options check if we will use recombining and for which species. Called by initracer. ! - su_recombin : Initialise tables. Called by sugas_corrk ! - call_recombin : Test profile of species and decide whether to call recombin_corrk. Called by callcork ! - recombin_corrk : The core algorithm properly recombining corrk tables. Called by callrecombin_corrk. ! ! TODO : Think about the case where nqtot phys /= nqtot_dyn !! ! TODO : Add the possibility to read an input profile for a specie !! ! Also think about the hybrid case where we want force with a latitudinally variable input profile (likely need to tweak on-the-fly) IMPLICIT NONE LOGICAL, SAVE :: corrk_recombin=.false. ! Key boolean, is there any specie to recombin. LOGICAL, SAVE :: use_premix=.true. ! Indicates if we recombin on top of a existing premix set of corrk. !$OMP THREADPRIVATE(corrk_recombin,use_premix) INTEGER, SAVE :: nrecomb_tot ! # (total) of compounds to recombine ! -> Linked to key is_recomb in tracer_h INTEGER, SAVE :: nrecomb_qset ! # of compounds to recombine having preset abundances (ie spectra computed with true vmr and not vmr=1) !-> Linked to key is_recomb_qset in tracer_h INTEGER, SAVE :: nrecomb_qotf ! # of compounds to recombine on-the-fly (ie using value of pq, not input profile) !-> Linked to key is_recomb_qotf in tracer_h !$OMP THREADPRIVATE(nrecomb_tot,nrecomb_qset,nrecomb_qotf) ! Note : The variables above are not read in callphys.def but automatically defined in inirecombin called in initracer after processing traceur.def LOGICAL, SAVE :: all_otf ! True if all species are recombined on-the-fly, no premix, no preset .. !$OMP THREADPRIVATE(all_otf) ! Following arrays are allocated in su_recombin (excepted otf2tot_idx, in ini_recombin) and deallocated in callcork lastcall REAL*8, save, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi_recomb, gasv_recomb REAL*8, save, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi_otf, gasv_otf REAL*8, save, DIMENSION(:), ALLOCATABLE :: w_cum REAL*8, save, DIMENSION(:), ALLOCATABLE :: wtwospec INTEGER, save, DIMENSION(:), ALLOCATABLE :: otf2tot_idx INTEGER, save, DIMENSION(:), ALLOCATABLE :: rcb2tot_idx INTEGER, save, DIMENSION(:), ALLOCATABLE :: otf2rcb_idx INTEGER, save, DIMENSION(:), ALLOCATABLE :: permut_idx !$OMP THREADPRIVATE(gasi_recomb,gasv_recomb,w_cum,otf2tot_idx,rcb2tot_idx,otf2rcb_idx,permut_idx) CONTAINS SUBROUTINE ini_recombin USE tracer_h IMPLICIT NONE INTEGER :: nspcrad ! # of is_rad species (tempooray here, should be a radcommon variabke) INTEGER :: iq, icount ! Default values use_premix = .true. corrk_recombin = .false. nrecomb_tot = 0 nrecomb_qset = 0 nrecomb_qotf = 0 nspcrad = 0 do iq=1,nqtot ! Counter used to check if all rad. species are recombin then no premix if(is_rad(iq)==1) nspcrad = nspcrad + 1 ! Sanity checks if (is_recomb(iq)==1 .and. is_rad(iq)==0) then write(*,*) 'initracer : Error for tracer iq=',iq write(*,*) 'is_recomb=1 but is_rad=0, this is not possible !' call abort_physic('initrac','Conflicting traceur.def options',1) endif if (is_recomb_qset(iq)==1 .and. is_recomb(iq)==0) then write(*,*) 'initracer : Error for tracer iq=',iq write(*,*) 'is_recomb_qset=1 but is_recomb=0, this is not possible !' call abort_physic('initrac','Conflicting traceur.def options',1) endif if (is_recomb_qotf(iq)==1 .and. is_recomb_qset(iq)==1) then write(*,*) 'initracer : Error for tracer iq=',iq write(*,*) 'is_recomb_qset=1 and is_recomb_qotf=1, this is not possible !' call abort_physic('initrac','Conflicting traceur.def options',1) endif if(is_recomb(iq)==1) then corrk_recombin = .true. ! activating on first one found would be enough actually but I'm lazy nrecomb_tot = nrecomb_tot + 1 if(is_recomb_qset(iq)==1) nrecomb_qset = nrecomb_qset + 1 if(is_recomb_qotf(iq)==1) nrecomb_qotf = nrecomb_qotf + 1 write(*,*) write(*,*) 'ini_recombin: found specie : Name = ',trim(noms(iq)), & ' ; Predefined vmr=', is_recomb_qset(iq), & ' ; On-the-fly=', is_recomb_qotf(iq) endif enddo ! Init. correspondance array of indexes between subset of tracers IF(.NOT. ALLOCATED(otf2tot_idx)) ALLOCATE(otf2tot_idx(nrecomb_qotf)) icount=0 do iq=1,nqtot if(is_recomb_qotf(iq)==1) then icount=icount+1 otf2tot_idx(icount) = iq endif enddo IF(.NOT. ALLOCATED(rcb2tot_idx)) ALLOCATE(rcb2tot_idx(nrecomb_tot)) icount=0 do iq=1,nqtot if(is_recomb(iq)==1) then icount=icount+1 rcb2tot_idx(icount) = iq endif enddo IF(.NOT. ALLOCATED(otf2rcb_idx)) ALLOCATE(otf2rcb_idx(nrecomb_qotf)) icount=0 do iq=1,nrecomb_tot if(is_recomb_qotf(rcb2tot_idx(iq))==1) then icount=icount+1 otf2rcb_idx(icount) = iq endif enddo ! Use a premix set on top ? if (nspcrad == nrecomb_tot .and. nspcrad /= 0) use_premix = .false. ! In this case all rad. species are recombined ! Summary write(*,*) write(*,*) 'ini_recombin: Total species found for corrk recombination =', nrecomb_tot if (corrk_recombin) then if (use_premix) then write(*,*) 'ini_recombin: .. Total radiative species matching total species for recombination...' write(*,*) 'ini_recombin: .. Any pre-mixed set of opacities will be ignored.' else write(*,*) 'ini_recombin: .. Found less species to recombine than total radiative species..' write(*,*) 'ini_recombin: .. Recombination will occur ontop of premix set of opacities' endif else write(*,*) 'ini_recombin: .. No species found for recombination, I will use premix set only.' endif write(*,*) END SUBROUTINE ini_recombin SUBROUTINE su_recombin USE radinc_h USE radcommon_h, only: gweight, gasi, gasv IMPLICIT NONE INTEGER :: i, ig, jg, ind, iw, it, ip ! Allocations IF(.NOT. ALLOCATED(permut_idx)) ALLOCATE(permut_idx(L_NGAUSS*L_NGAUSS)) IF(.NOT. ALLOCATED(w_cum)) ALLOCATE(w_cum(L_NGAUSS)) IF(.NOT. ALLOCATED(gasi_recomb)) ALLOCATE(gasi_recomb(L_NTREF,L_PINT,L_REFVAR,L_NSPECTI,L_NGAUSS)) IF(.NOT. ALLOCATED(gasv_recomb)) ALLOCATE(gasv_recomb(L_NTREF,L_PINT,L_REFVAR,L_NSPECTV,L_NGAUSS)) IF(.NOT. ALLOCATED(gasi_otf)) ALLOCATE(gasi_otf(L_NGAUSS,nrecomb_qotf,L_NSPECTI,L_NTREF,L_PINT)) IF(.NOT. ALLOCATED(gasv_otf)) ALLOCATE(gasv_otf(L_NGAUSS,nrecomb_qotf,L_NSPECTI,L_NTREF,L_PINT)) IF(.NOT. ALLOCATED(wtwospec)) ALLOCATE(wtwospec(L_NGAUSS*L_NGAUSS)) ! Init. for recombin_corrk firstcall permut_idx = (/(i, i=1,L_NGAUSS*L_NGAUSS)/) w_cum(1)= gweight(1) DO i=2,L_NGAUSS w_cum(i) = w_cum(i-1)+gweight(i) ENDDO ! init wtwospec once for all DO ig=1,L_NGAUSS DO jg=1, L_NGAUSS ind = jg+(ig-1)*L_NGAUSS wtwospec(ind) = gweight(ig)*gweight(jg) ENDDO ENDDO ! init otf correlated-k array do ip=1,L_PINT do it=1,L_NTREF do iw=1,L_NSPECTI do i=1,nrecomb_qotf do ig=1,L_NGAUSS 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 enddo enddo enddo enddo enddo do ip=1,L_PINT do it=1,L_NTREF do iw=1,L_NSPECTV do i=1,nrecomb_qotf do ig=1,L_NGAUSS 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 enddo enddo enddo enddo enddo gasi_recomb(:,:,:,:,:) = gasi(:,:,1:L_REFVAR,:,:) ! non-zero init (=kappa_ir) gasv_recomb(:,:,:,:,:) = gasv(:,:,1:L_REFVAR,:,:) ! non-zero init (=kappa_vi) END SUBROUTINE su_recombin SUBROUTINE call_recombin(igrid,nlayer,pq,pplay,pt,qvar,tmid,pmid) USE comcstfi_mod, only: mugaz USE radinc_h USE radcommon_h USE tracer_h, only: noms, mmol, is_recomb_qotf, is_recomb_qset USE tpindex_mod, only: tpindex IMPLICIT NONE ! Inputs INTEGER, INTENT(IN) :: igrid ! lon-lat grid point INTEGER, INTENT(IN) :: nlayer ! Number of atmospheric layers. REAL*8, DIMENSION(:,:), INTENT(IN) :: pq ! Tracers vertical profiles (kg/kg) REAL*8, DIMENSION(nlayer), INTENT(IN) :: pplay ! Atmospheric pressure (Pa) REAL*8, DIMENSION(nlayer), INTENT(IN) :: pt ! Atmospheric temperature (K) REAL*8, DIMENSION(L_LEVELS), INTENT(IN) :: qvar ! Mixing ratio of variable component (mol/mol) REAL*8, DIMENSION(L_LEVELS), INTENT(IN) :: tmid ! Temperature of layers, mid levels (K) REAL*8, DIMENSION(L_LEVELS), INTENT(IN) :: pmid ! Pressure of layers, mid levels (mBar) ! NB : qvar is on L_LEVELS but it has been processed in callcork compared to the one in pq ! so we imperatively need to take this one. Note that there is no interpolation, so ! pq(nlayer+1-l,ivar) is broadcast in qvar(2*l) qvar(2*l+1) ! Local variables INTEGER :: ng INTEGER :: ig,l,k,nw,iq,ip,ilay,it,ix,iw LOGICAL :: found REAL*8 :: fact, tmin, tmax, qmin, qmax REAL*8 :: LCOEF(4), WRATIO LOGICAL,DIMENSION(:,:,:),ALLOCATABLE,save :: useptx ! Mask on which t-p-x ref grid point will be used REAL*8, DIMENSION(:,:), ALLOCATABLE,save :: pqr ! Tracers abundances at ref pressures used for onthefly recombining (mol/mol). LOGICAL,SAVE :: firstcall=.true. !$OMP THREADPRIVATE(firstcall) ! At firstcall we precombine all what needs to be done only once (pre-mixed,forced profiles..), if needed. IF (firstcall) THEN IF(.NOT. ALLOCATED(useptx)) ALLOCATE(useptx(L_NTREF,L_PINT,L_REFVAR)) useptx(:,:,:) = .false. IF(use_premix .or. (.not.use_premix .and. nrecomb_qotf/=nrecomb_tot)) THEN ! we skip this if all species are on-the-fly all_otf=.false. IF(.NOT. ALLOCATED(pqr)) ALLOCATE(pqr(nrecomb_tot,L_PINT)) ! Default value for premix and for fixed species for whom vmr has been taken ! into account while computing high-resolution spectra pqr(:,:) = 1.0 ! TODO : Missing implementation here for the tracers where we read an input profile !! do iq=1,nrecomb_tot if (is_recomb_qset(rcb2tot_idx(iq))==0 .and. is_recomb_qotf(rcb2tot_idx(iq))==0) then print*, 'Recombining tracer ', noms(rcb2tot_idx(iq)),' requires an input profile, this is not implemented yet !!' call abort_physic('call_recombin','Missing implementation',1) ! Read pqr(:,iq) endif enddo ! Recombine for all T-P-Q as we do it only once for all. call recombin_corrk_ini(pqr) ELSE all_otf=.true. ENDIF firstcall=.false. IF (nrecomb_qotf==0) corrk_recombin = .false. IF(ALLOCATED(pqr)) DEALLOCATE(pqr) IF(.NOT. ALLOCATED(pqr)) ALLOCATE(pqr(nrecomb_qotf,L_PINT)) ENDIF ! firstcall ! NB : To have decent CPU time recombining is not done on all gridpoints and wavelenghts but we ! calculate a gasi/v_recomb variable on the reference corrk-k T,P,X grid (only for T,P,X values ! who match the atmospheric conditions) which is then processed as a standard pre-mix in ! optci/v routines, but updated every time tracers on the ref P grid have varied > 1%. ! Extract tracers for species which are recombined on-the-fly do ip=1,L_PINT ilay=0 found = .false. do l=1,nlayer if (pplay(l) .le. 10.0**(pfgasref(ip)+2.0)) then ! pfgasref=log(p[mbar]) found=.true. ilay=l-1 exit endif enddo if (.not. found) then ! set pq to top value do iq=1,nrecomb_qotf pqr(iq,ip) = pq(nlayer,otf2tot_idx(iq))*mugaz/mmol(otf2tot_idx(iq)) ! mol/mol enddo else if (ilay==0) then ! set pq to bottom value do iq=1,nrecomb_qotf pqr(iq,ip) = pq(1,otf2tot_idx(iq))*mugaz/mmol(otf2tot_idx(iq)) ! mol/mol enddo else ! standard : interp pq between layers fact = (10.0**(pfgasref(ip)+2.0) - pplay(ilay+1)) / (pplay(ilay) - pplay(ilay+1)) ! pfgasref=log(p[mbar]) do iq=1,nrecomb_qotf pqr(iq,ip) = pq(ilay,otf2tot_idx(iq))**fact * pq(ilay+1,otf2tot_idx(iq))**(1.0-fact) pqr(iq,ip) = pqr(iq,ip)*mugaz/mmol(otf2tot_idx(iq)) ! mol/mol enddo endif ! if ilay==nlayer endif ! if not found enddo ! ip=1,L_PINT ! The following useptx is a trick to call recombine only for the reference T-P-X ! reference grid points that are useful given the temperature range (and variable specie amount) at a given altitude ! (cf optci/optcv routines where we interpolate corrk calling tpindex) ! It saves a looot of time - JVO 18 do K=2,L_LEVELS call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR,LCOEF,it,ip,ix,WRATIO) useptx(it:it+1,ip:ip+1,ix:ix+1) = .true. end do if (.not.all_otf) then call recombin_corrk_mix(igrid,pqr,useptx) else if (nrecomb_qotf.gt.1) then call recombin_corrk_mix_allotf(igrid,pqr,useptx) else do ix=1,L_REFVAR do ip=1,L_PINT do it=1,L_NTREF if (.not. useptx(it,ip,ix)) cycle gasi_recomb(it,ip,ix,:,:) = pqr(1,ip)*gasi_otf(:,1,:,it,ip) gasv_recomb(it,ip,ix,:,:) = pqr(1,ip)*gasv_otf(:,1,:,it,ip) useptx(it,ip,ix) = .false. enddo enddo enddo endif endif END SUBROUTINE call_recombin SUBROUTINE recombin_corrk_ini(pqr) USE radinc_h USE radcommon_h, only: gweight, tlimit, gasi, gasv USE tracer_h, only: is_recomb_qotf USE sort_mod, only: isort !----------------------------------------------------------------------- ! Declarations: ! ------------- IMPLICIT NONE ! Arguments : ! ----------- REAL*8, DIMENSION(nrecomb_qotf,L_PINT), INTENT(IN) :: pqr ! otf species mixing ration ! Local variables : ! ----------------- INTEGER :: it, ip, ix, iw, ig, jg, ind, ibin, ispec REAL*8, DIMENSION(L_NGAUSS) :: tmpk ! otf correlated-k by mixing ratio REAL*8, DIMENSION(L_NGAUSS) :: krecomb REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS) :: ktwospec REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS) :: ktwospec_s REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS) :: wtwospec_s REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS) :: wtwospec_cum REAL*8 :: wsplit do ix=1,L_REFVAR do ip=1,L_PINT do it=1,L_NTREF ! ------------------- ! I. INFRARED ! ------------------- DO iw=1,L_NSPECTI DO ig=1,L_NGAUSS ! init correlated-k with premix ! utiliser directement gasi_recomb au lieu d'un variable intermediere krecomb ? ! Peut ok si L_NGAUSS première dimension de gasi_recomb krecomb(ig) = gasi(it,ip,ix,iw,ig) ENDDO DO ispec=1,nrecomb_tot ! Loop on additional species IF(is_recomb_qotf(rcb2tot_idx(ispec))==1) CYCLE ! takes all to recomb, the otf ones are skipped DO ig=1,L_NGAUSS tmpk(ig) = pqr(ispec,ip)*gasi(it,ip,L_REFVAR+ispec,iw,ig) ENDDO ! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band IF ( tmpk(L_NGAUSS-1) .LE. tlimit ) CYCLE IF ( ALL( tmpk(1:L_NGAUSS-1) .LE. krecomb(1:L_NGAUSS-1)*0.1 ) ) CYCLE IF ( ALL( krecomb(1:L_NGAUSS-1) .LE. tmpk(1:L_NGAUSS-1)*0.1 ) ) THEN krecomb(1:L_NGAUSS-1) = tmpk(1:L_NGAUSS-1) CYCLE ENDIF ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO ig=1,L_NGAUSS DO jg=1, L_NGAUSS ind = jg+(ig-1)*L_NGAUSS ktwospec(ind) = krecomb(ig)+tmpk(jg) ENDDO ENDDO ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting ! NB : quite small array, quicker to permut with 2nd array than in place DO ind=1,L_NGAUSS*L_NGAUSS ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall ENDDO CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Insertion sort quicker because pre-sorted ! Sort w according to permutations of k. ! NB : quite small array, quicker to permut with 2nd array than in place DO ind=1,L_NGAUSS*L_NGAUSS wtwospec_s(ind) = wtwospec(permut_idx(ind)) ENDDO wtwospec_cum(1) = wtwospec_s(1) DO ind=2,L_NGAUSS*L_NGAUSS wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind) ENDDO ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ibin=1 krecomb(:)=0.0 DO ig=1, L_NGAUSS-1 DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN wsplit = w_cum(ig) - wtwospec_cum(ind-1) krecomb(ig) = krecomb(ig) & + sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) & + wsplit*ktwospec_s(ind) krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind) ibin=ind+1 EXIT ENDIF ENDDO krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) ) ENDDO krecomb(1:L_NGAUSS-1) = krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0 ENDDO ! ispec=1,nrecomb_qotf gasi(it,ip,ix,iw,:) = krecomb(:) ENDDO ! iw=1,L_NSPECTI ! ------------------- ! II. VISIBLE ! ------------------- DO iw=1,L_NSPECTV DO ig=1,L_NGAUSS ! init correlated-k with premix krecomb(ig) = gasv(it,ip,ix,iw,ig) ! there is a prerecombined cocktail ENDDO DO ispec=1,nrecomb_tot ! Loop on additional species IF(is_recomb_qotf(rcb2tot_idx(ispec))==1) CYCLE ! takes all to recomb, the otf ones are skipped DO ig=1,L_NGAUSS tmpk(ig) = pqr(ispec,ip)*gasv(it,ip,L_REFVAR+ispec,iw,ig) ENDDO ! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band IF ( tmpk(L_NGAUSS-1) .LE. tlimit ) CYCLE IF ( ALL( tmpk(1:L_NGAUSS-1) .LE. krecomb(1:L_NGAUSS-1)*0.1 ) ) CYCLE IF ( ALL( krecomb(1:L_NGAUSS-1) .LE. tmpk(1:L_NGAUSS-1)*0.1 ) ) THEN krecomb(1:L_NGAUSS-1) = tmpk(1:L_NGAUSS-1) CYCLE ENDIF ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO ig=1,L_NGAUSS DO jg=1, L_NGAUSS ind = jg+(ig-1)*L_NGAUSS ktwospec(ind) = krecomb(ig)+tmpk(jg) ENDDO ENDDO ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting ! NB : quite small array, quicker to permut with 2nd array than in place DO ind=1,L_NGAUSS*L_NGAUSS ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall ENDDO CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Insertion sort quicker because pre-sorted ! Sort w according to permutations of k. ! NB : quite small array, quicker to permut with 2nd array than in place DO ind=1,L_NGAUSS*L_NGAUSS wtwospec_s(ind) = wtwospec(permut_idx(ind)) ENDDO wtwospec_cum(1) = wtwospec_s(1) DO ind=2,L_NGAUSS*L_NGAUSS wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind) ENDDO ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ibin=1 krecomb(:)=0.0 DO ig=1, L_NGAUSS-1 DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN wsplit = w_cum(ig) - wtwospec_cum(ind-1) krecomb(ig) = krecomb(ig) & + sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) & + wsplit*ktwospec_s(ind) krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind) ibin=ind+1 EXIT ENDIF ENDDO krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) ) ENDDO krecomb(1:L_NGAUSS-1) = krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0 ENDDO ! ispec=1,nrecomb_qotf gasv(it,ip,ix,iw,:) = krecomb(:) ENDDO ! iw=1,L_NSPECTV enddo ! ix=1,L_NTREF enddo ! it=1,L_PINT enddo ! ip=1,L_REFVAR END SUBROUTINE recombin_corrk_ini SUBROUTINE recombin_corrk_mix(igrid,pqr,useptx) USE radinc_h USE radcommon_h, only: gweight, tlimit, gasi, gasv USE sort_mod, only: isort use comsaison_h, only: fract !----------------------------------------------------------------------- ! Declarations: ! ------------- IMPLICIT NONE ! Arguments : ! ----------- INTEGER, INTENT(IN) :: igrid ! lon-lat grid point REAL*8, DIMENSION(nrecomb_qotf,L_PINT), INTENT(IN) :: pqr ! otf species mixing ration LOGICAL, INTENT(INOUT) :: useptx(L_NTREF,L_PINT,L_REFVAR) ! Mask on which t-p-x ref grid point will be used ! Local variables : ! ----------------- INTEGER :: it, ip, ix, iw, ig, jg, ind, ibin, ispec REAL*8, DIMENSION(L_NGAUSS) :: tmpk ! otf correlated-k by mixing ratio REAL*8, DIMENSION(L_NGAUSS) :: krecomb REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS) :: ktwospec REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS) :: ktwospec_s REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS) :: wtwospec_s REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS) :: wtwospec_cum REAL*8 :: wsplit do ix=1,L_REFVAR do ip=1,L_PINT do it=1,L_NTREF if (.not. useptx(it,ip,ix)) cycle ! ------------------- ! I. INFRARED ! ------------------- DO iw=1,L_NSPECTI DO ig=1,L_NGAUSS ! init correlated-k with premix ! utiliser directement gasi_recomb au lieu d'un variable intermediere krecomb ? ! Peut ok si L_NGAUSS première dimension de gasi_recomb krecomb(ig) = gasi(it,ip,ix,iw,ig) ENDDO DO ispec=1,nrecomb_qotf ! Loop on additional species DO ig=1,L_NGAUSS tmpk(ig) = pqr(ispec,ip)*gasi_otf(ig,ispec,iw,it,ip) ENDDO ! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band IF ( tmpk(L_NGAUSS-1) .LE. tlimit ) CYCLE IF ( ALL( tmpk(1:L_NGAUSS-1) .LE. krecomb(1:L_NGAUSS-1)*0.1 ) ) CYCLE IF ( ALL( krecomb(1:L_NGAUSS-1) .LE. tmpk(1:L_NGAUSS-1)*0.1 ) ) THEN krecomb(1:L_NGAUSS-1) = tmpk(1:L_NGAUSS-1) CYCLE ENDIF ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO ig=1,L_NGAUSS DO jg=1, L_NGAUSS ind = jg+(ig-1)*L_NGAUSS ktwospec(ind) = krecomb(ig)+tmpk(jg) ENDDO ENDDO ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting ! NB : quite small array, quicker to permut with 2nd array than in place DO ind=1,L_NGAUSS*L_NGAUSS ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall ENDDO CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Insertion sort quicker because pre-sorted ! Sort w according to permutations of k. ! NB : quite small array, quicker to permut with 2nd array than in place DO ind=1,L_NGAUSS*L_NGAUSS wtwospec_s(ind) = wtwospec(permut_idx(ind)) ENDDO wtwospec_cum(1) = wtwospec_s(1) DO ind=2,L_NGAUSS*L_NGAUSS wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind) ENDDO ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ibin=1 krecomb(:)=0.0 DO ig=1, L_NGAUSS-1 DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN wsplit = w_cum(ig) - wtwospec_cum(ind-1) krecomb(ig) = krecomb(ig) & + sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) & + wsplit*ktwospec_s(ind) krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind) ibin=ind+1 EXIT ENDIF ENDDO krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) ) ENDDO krecomb(1:L_NGAUSS-1) = krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0 ENDDO ! ispec=1,nrecomb_qotf gasi_recomb(it,ip,ix,iw,:) = krecomb(:) ENDDO ! iw=1,L_NSPECTI ! ------------------- ! II. VISIBLE ! ------------------- if(fract(igrid) .lt. 1.0e-4) cycle ! Only during daylight. DO iw=1,L_NSPECTV DO ig=1,L_NGAUSS ! init correlated-k with premix krecomb(ig) = gasv(it,ip,ix,iw,ig) ! gasv_loc order correctly for running time? ENDDO DO ispec=1,nrecomb_qotf ! Loop on additional species DO ig=1,L_NGAUSS tmpk(ig) = pqr(ispec,ip)*gasv_otf(ig,ispec,iw,it,ip) ENDDO ! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band IF ( tmpk(L_NGAUSS-1) .LE. tlimit ) CYCLE IF ( ALL( tmpk(1:L_NGAUSS-1) .LE. krecomb(1:L_NGAUSS-1)*0.1 ) ) CYCLE IF ( ALL( krecomb(1:L_NGAUSS-1) .LE. tmpk(1:L_NGAUSS-1)*0.1 ) ) THEN krecomb(1:L_NGAUSS-1) = tmpk(1:L_NGAUSS-1) CYCLE ENDIF ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO ig=1,L_NGAUSS DO jg=1, L_NGAUSS ind = jg+(ig-1)*L_NGAUSS ktwospec(ind) = krecomb(ig)+tmpk(jg) ENDDO ENDDO ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting ! NB : quite small array, quicker to permut with 2nd array than in place DO ind=1,L_NGAUSS*L_NGAUSS ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall ENDDO CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Insertion sort quicker because pre-sorted ! Sort w according to permutations of k. ! NB : quite small array, quicker to permut with 2nd array than in place DO ind=1,L_NGAUSS*L_NGAUSS wtwospec_s(ind) = wtwospec(permut_idx(ind)) ENDDO wtwospec_cum(1) = wtwospec_s(1) DO ind=2,L_NGAUSS*L_NGAUSS wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind) ENDDO ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ibin=1 krecomb(:)=0.0 DO ig=1, L_NGAUSS-1 DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN wsplit = w_cum(ig) - wtwospec_cum(ind-1) krecomb(ig) = krecomb(ig) & + sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) & + wsplit*ktwospec_s(ind) krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind) ibin=ind+1 EXIT ENDIF ENDDO krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) ) ENDDO krecomb(1:L_NGAUSS-1) = krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0 ENDDO ! ispec=1,nrecomb_qotf gasv_recomb(it,ip,ix,iw,:) = krecomb(:) ENDDO ! iw=1,L_NSPECTV useptx(it,ip,ix) = .false. enddo ! ix=1,L_NTREF enddo ! it=1,L_PINT enddo ! ip=1,L_REFVAR END SUBROUTINE recombin_corrk_mix SUBROUTINE recombin_corrk_mix_allotf(igrid,pqr,useptx) USE radinc_h USE radcommon_h, only: gweight, tlimit, gasi, gasv USE sort_mod, only: isort use comsaison_h, only: fract !----------------------------------------------------------------------- ! Declarations: ! ------------- IMPLICIT NONE ! Arguments : ! ----------- INTEGER, INTENT(IN) :: igrid ! lon-lat grid point REAL*8, DIMENSION(nrecomb_qotf,L_PINT), INTENT(IN) :: pqr ! otf species mixing ration LOGICAL, INTENT(INOUT) :: useptx(L_NTREF,L_PINT,L_REFVAR) ! Mask on which t-p-x ref grid point will be used ! Local variables : ! ----------------- INTEGER :: it, ip, ix, iw, ig, jg, ind, ibin, ispec REAL*8, DIMENSION(L_NGAUSS) :: tmpk ! otf correlated-k by mixing ratio REAL*8, DIMENSION(L_NGAUSS) :: krecomb REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS) :: ktwospec REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS) :: ktwospec_s REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS) :: wtwospec_s REAL*8, DIMENSION(L_NGAUSS*L_NGAUSS) :: wtwospec_cum REAL*8 :: wsplit do ix=1,L_REFVAR do ip=1,L_PINT do it=1,L_NTREF if (.not. useptx(it,ip,ix)) cycle ! ------------------- ! I. INFRARED ! ------------------- DO iw=1,L_NSPECTI DO ig=1,L_NGAUSS ! init correlated-k with first gas krecomb(ig) = pqr(1,ip)*gasi_otf(ig,1,iw,it,ip) ENDDO DO ispec=2,nrecomb_qotf ! Loop on additional species DO ig=1,L_NGAUSS tmpk(ig) = pqr(ispec,ip)*gasi_otf(ig,ispec,iw,it,ip) ENDDO ! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band IF ( tmpk(L_NGAUSS-1) .LE. tlimit ) CYCLE IF ( ALL( tmpk(1:L_NGAUSS-1) .LE. krecomb(1:L_NGAUSS-1)*0.01 ) ) CYCLE IF ( ALL( krecomb(1:L_NGAUSS-1) .LE. tmpk(1:L_NGAUSS-1)*0.01 ) ) THEN krecomb(1:L_NGAUSS-1) = tmpk(1:L_NGAUSS-1) CYCLE ENDIF ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO ig=1,L_NGAUSS DO jg=1, L_NGAUSS ind = jg+(ig-1)*L_NGAUSS ktwospec(ind) = krecomb(ig)+tmpk(jg) ENDDO ENDDO ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting ! NB : quite small array, quicker to permut with 2nd array than in place DO ind=1,L_NGAUSS*L_NGAUSS ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall ENDDO CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Insertion sort quicker because pre-sorted ! Sort w according to permutations of k. ! NB : quite small array, quicker to permut with 2nd array than in place DO ind=1,L_NGAUSS*L_NGAUSS wtwospec_s(ind) = wtwospec(permut_idx(ind)) ENDDO wtwospec_cum(1) = wtwospec_s(1) DO ind=2,L_NGAUSS*L_NGAUSS wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind) ENDDO ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ibin=1 krecomb(:)=0.0 DO ig=1, L_NGAUSS-1 DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN wsplit = w_cum(ig) - wtwospec_cum(ind-1) krecomb(ig) = krecomb(ig) & + sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) & + wsplit*ktwospec_s(ind) krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind) ibin=ind+1 EXIT ENDIF ENDDO krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) ) ENDDO krecomb(1:L_NGAUSS-1) = krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0 ENDDO ! ispec=1,nrecomb_qotf gasi_recomb(it,ip,ix,iw,:) = krecomb(:) ENDDO ! iw=1,L_NSPECTI ! ------------------- ! II. VISIBLE ! ------------------- if(fract(igrid) .lt. 1.0e-4) cycle ! Only during daylight. DO iw=1,L_NSPECTV DO ig=1,L_NGAUSS ! init correlated-k with first gas krecomb(ig) = pqr(1,ip)*gasv_otf(ig,1,iw,it,ip) ENDDO DO ispec=2,nrecomb_qotf ! Loop on additional species DO ig=1,L_NGAUSS tmpk(ig) = pqr(ispec,ip)*gasv_otf(ig,ispec,iw,it,ip) ENDDO ! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band IF ( tmpk(L_NGAUSS-1) .LE. tlimit ) CYCLE IF ( ALL( tmpk(1:L_NGAUSS-1) .LE. krecomb(1:L_NGAUSS-1)*0.01 ) ) CYCLE IF ( ALL( krecomb(1:L_NGAUSS-1) .LE. tmpk(1:L_NGAUSS-1)*0.01 ) ) THEN krecomb(1:L_NGAUSS-1) = tmpk(1:L_NGAUSS-1) CYCLE ENDIF ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO ig=1,L_NGAUSS DO jg=1, L_NGAUSS ind = jg+(ig-1)*L_NGAUSS ktwospec(ind) = krecomb(ig)+tmpk(jg) ENDDO ENDDO ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting ! NB : quite small array, quicker to permut with 2nd array than in place DO ind=1,L_NGAUSS*L_NGAUSS ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall ENDDO CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Insertion sort quicker because pre-sorted ! Sort w according to permutations of k. ! NB : quite small array, quicker to permut with 2nd array than in place DO ind=1,L_NGAUSS*L_NGAUSS wtwospec_s(ind) = wtwospec(permut_idx(ind)) ENDDO wtwospec_cum(1) = wtwospec_s(1) DO ind=2,L_NGAUSS*L_NGAUSS wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind) ENDDO ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ibin=1 krecomb(:)=0.0 DO ig=1, L_NGAUSS-1 DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN wsplit = w_cum(ig) - wtwospec_cum(ind-1) krecomb(ig) = krecomb(ig) & + sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) & + wsplit*ktwospec_s(ind) krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind) ibin=ind+1 EXIT ENDIF ENDDO krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) ) ENDDO krecomb(1:L_NGAUSS-1) = krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0 ENDDO ! ispec=1,nrecomb_qotf gasv_recomb(it,ip,ix,iw,:) = krecomb(:) ENDDO ! iw=1,L_NSPECTV useptx(it,ip,ix) = .false. enddo ! ix=1,L_NTREF enddo ! it=1,L_PINT enddo ! ip=1,L_REFVAR END SUBROUTINE recombin_corrk_mix_allotf END MODULE recombin_corrk_mod