source: trunk/LMDZ.TITAN/libf/phytitan/recombin_corrk.F90 @ 2241

Last change on this file since 2241 was 2050, checked in by jvatant, 7 years ago

Major radiative transfer contribution : Add the 'corrk_recombin' option that allows to use
correlated-k for single species instead of pre-mix and enables more flexiblity for variable species.
-> Algorithm inspired from Lacis and Oinas 1991 and Amundsen et al 2016

+ Added 'recombin_corrk.F90' - Important improvements in 'sugas_corrk.90' and 'callcork.F90'

  • Must have the desired variable species as tracers -> TBD : Enable to force composition even if no tracers
  • 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 grid (only for T,P 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%.


READ CAREFULY :

  • In case of 'corrk_recombin', the variable L_NREFVAR doesn't have the same meaning as before and doesn't stand for the different mixing ratios but the different species.
  • Input corr-k should be found in corrkdir within 'corrk_gcm_IR/VI_XXX.dat' and can contain a 'fixed' specie ( compulsory if you include self-broadening ) that MUST have been created with correct mixing ratios, or a variable specie for which mixing ratio MUST have been set to 1 ( no self-broadening then, assume it's a trace sepecie ) -> You can't neither have CIA of variable species included upstream in the corr-k
File size: 8.3 KB
Line 
1SUBROUTINE recombin_corrk(q,ip,it)
2
3  !     ===================================================================
4  !     Purpose
5  !     -------
6  !     Recombine correlated-k of individual species in a unique corr-k.
7  !
8  !       -> See Lacis and Oinas (1991), Amundsen et al (2016).
9  !
10  !     One important hypothesis is that we assume random overlap
11  !     Steps are - 1. Recombining ( Convolution spectra of two species )
12  !                 2. Resorting the ngauss*ngauss obtained data
13  !                 3. Rebin down to a smaller amount of bin ( here ngauss )
14  !                 4. Start again if another specie have to be included
15  !
16  !     In each band computation is done only for species with enough optical depth.
17  !     
18  !     Authors
19  !     -------
20  !     J. Vatant d'Ollone (2018)
21  !     ====================================================================
22
23  USE radinc_h
24  USE radcommon_h, only: gasi, gasv, gasi_recomb, gasv_recomb, &
25       gweight, pqrold, permut_idx, tlimit, w_cum
26  USE sort_mod, only: qsort, isort
27
28  !-----------------------------------------------------------------------
29  !     Declarations:
30  !     -------------
31
32  IMPLICIT NONE
33
34  !  Arguments :
35  !  -----------
36  REAL, DIMENSION(L_REFVAR), INTENT(IN)  :: q
37
38  INTEGER,                   INTENT(IN)  :: ip
39  INTEGER,                   INTENT(IN)  :: it
40
41  !  Local variables :
42  !  -----------------
43  INTEGER :: iw, ispec, ig, jg, ind, ibin
44
45  REAL, DIMENSION(L_NGAUSS)             :: krecomb
46
47  REAL, DIMENSION(L_NGAUSS*L_NGAUSS)    :: ktwospec
48  REAL, DIMENSION(L_NGAUSS*L_NGAUSS)    :: ktwospec_s
49  REAL, DIMENSION(L_NGAUSS*L_NGAUSS)    :: wtwospec
50  REAL, DIMENSION(L_NGAUSS*L_NGAUSS)    :: wtwospec_s
51  REAL, DIMENSION(L_NGAUSS*L_NGAUSS)    :: wtwospec_cum 
52
53  REAL, DIMENSION(L_REFVAR,L_NGAUSS)    :: tmpk
54
55  REAL :: wsplit
56
57
58  ! -------------------
59  ! I. INFRARED
60  ! -------------------
61
62  DO iw=1,L_NSPECTI
63
64     krecomb(:) = q(1)*gasi(it,ip,1,iw,:) ! init for > 1 loop and works also if only one active specie
65
66     tmpk(:,:) = gasi(it,ip,:,iw,:)
67
68     IF ( L_REFVAR .GT. 1 ) THEN
69
70        DO ispec=2,L_REFVAR
71
72           ! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band
73           IF ( q(ispec)*tmpk(ispec,L_NGAUSS-1) .LE. tlimit ) CYCLE
74           !              IF ( tmpk(ispec,L_NGAUSS-1) .LE. tlimit ) CYCLE
75
76           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
77           ! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~
78           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
79           DO ig=1,L_NGAUSS
80              DO jg=1, L_NGAUSS
81                 ind = jg+(ig-1)*L_NGAUSS
82                 ktwospec(ind) = krecomb(ig)+q(ispec)*tmpk(ispec,jg)
83                 wtwospec(ind) = gweight(ig)*gweight(jg)
84              ENDDO
85           ENDDO
86
87           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
88           ! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~
89           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
90           
91           ! Pre-sort from last step ( we have always a similar regular pattern ) to gain time for sorting
92           ! NB : quite small array, quicker to permut with 2nd array than in place
93           DO ind=1,L_NGAUSS*L_NGAUSS
94              ktwospec_s(ind) = ktwospec(permut_idx(ind)) ! NB : won't do anything at firstcall
95           ENDDO
96
97           CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx)  ! Insertion sort quicker because pre-sorted
98           !CALL qsort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Quicksort slower for pre-sorted
99
100           ! Sort w according to permutations of k.
101           ! NB : quite small array, quicker to permut with 2nd array than in place
102           DO ind=1,L_NGAUSS*L_NGAUSS
103              wtwospec_s(ind) = wtwospec(permut_idx(ind))
104           ENDDO
105
106           wtwospec_cum(1) = wtwospec_s(1)
107           DO ind=2,L_NGAUSS*L_NGAUSS
108              wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind)
109           ENDDO
110
111           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
112           ! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~
113           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
114           ibin=1
115
116           krecomb(:)=0.0
117
118           DO ig=1, L_NGAUSS-1
119
120              DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while   
121                 IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN
122                    wsplit =  w_cum(ig) - wtwospec_cum(ind-1)
123                    krecomb(ig)   = krecomb(ig)                                            &
124                         + sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) &
125                         + wsplit*ktwospec_s(ind)
126                    krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind)
127                    ibin=ind+1
128                    EXIT
129                 ENDIF
130              ENDDO
131
132              krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) )
133
134           ENDDO
135
136           krecomb(1:L_NGAUSS-1) =  krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0
137
138        ENDDO ! ispec=2,L_REFVAR
139
140     ENDIF ! if L_REFVAR .GT. 1
141
142     gasi_recomb(it,ip,iw,:) = krecomb(:)
143
144  ENDDO ! iw=1,L_NSPECTI
145
146
147  ! ----------------
148  ! II. VISIBLE
149  ! ----------------
150
151  DO iw=1,L_NSPECTV
152
153     krecomb(:) = q(1)*gasv(it,ip,1,iw,:) ! init for > 1 loop and works also if only one active specie
154     
155     tmpk(:,:) = gasv(it,ip,:,iw,:)
156
157     IF ( L_REFVAR .GT. 1 ) THEN
158
159        DO ispec=2,L_REFVAR
160
161           ! Save ( a lot of ) CPU time, we don't add the specie if negligible absorption in the band
162           IF ( q(ispec)*tmpk(ispec,L_NGAUSS-1) .LE. tlimit ) CYCLE
163           !              IF ( tmpk(ispec,L_NGAUSS-1) .LE. tlimit ) CYCLE
164
165           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
166           ! 1. Recombining ~~~~~~~~~~~~~~~~~~~~~~~~~
167           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
168           DO ig=1,L_NGAUSS
169              DO jg=1, L_NGAUSS
170                 ind = jg+(ig-1)*L_NGAUSS
171                 ktwospec(ind) = krecomb(ig)+q(ispec)*tmpk(ispec,jg)
172                 wtwospec(ind) = gweight(ig)*gweight(jg)
173              ENDDO
174           ENDDO
175
176           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
177           ! 2. Resorting ~~~~~~~~~~~~~~~~~~~~~~~
178           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
179
180           ! Pre-sort from last call ( we have always a similar regular pattern ) to gain time in sorting
181           ! NB : quite small array, quicker to permut with 2nd array than in place
182           DO ind=1,L_NGAUSS*L_NGAUSS
183              ktwospec_s(ind) = ktwospec(permut_idx(ind))
184           ENDDO
185
186           CALL isort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx)  ! Insertion sort quicker because pre-sorted
187           !CALL qsort(ktwospec_s,L_NGAUSS*L_NGAUSS,permut_idx) ! Quicksort slower for pre-sorted
188
189           ! Sort w according to permutations of k.
190           ! NB : quite small array, quicker to permut with copy than in place
191           DO ind=1,L_NGAUSS*L_NGAUSS
192              wtwospec_s(ind) = wtwospec(permut_idx(ind))
193           ENDDO
194
195           wtwospec_cum(1) = wtwospec_s(1)
196           DO ind=2,L_NGAUSS*L_NGAUSS
197              wtwospec_cum(ind)= wtwospec_cum(ind-1)+wtwospec_s(ind)
198           ENDDO
199
200           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
201           ! 3. Rebinning on smaller amount of Gauss points ~~~~~~~~~~~
202           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
203           ibin=1
204
205           krecomb(:)=0.0
206
207           DO ig=1, L_NGAUSS-1
208
209              DO ind=ibin,L_NGAUSS*L_NGAUSS ! rather than a while   
210                 IF ( wtwospec_cum(ind) .GT. w_cum(ig) ) THEN
211                    wsplit =  w_cum(ig) - wtwospec_cum(ind-1)
212                    krecomb(ig)   = krecomb(ig)                                             &
213                         + sum ( wtwospec_s(ibin:ind-1)*ktwospec_s(ibin:ind-1) ) &
214                         + wsplit*ktwospec_s(ind)
215                    krecomb(ig+1) = (wtwospec_s(ind)-wsplit)*ktwospec_s(ind)
216                    ibin=ind+1
217                    EXIT
218                 ENDIF
219              ENDDO
220
221              krecomb(L_NGAUSS) = krecomb(L_NGAUSS) + sum ( wtwospec_s(ibin:)*ktwospec_s(ibin:) )
222
223           ENDDO
224
225           krecomb(1:L_NGAUSS-1) =  krecomb(1:L_NGAUSS-1) / gweight(1:L_NGAUSS-1) ! gw(L_NGAUSS)=0
226
227        ENDDO ! ispec=2,L_REFVAR
228
229     ENDIF ! if L_REFVAR .GT. 1
230
231     gasv_recomb(it,ip,iw,:) = krecomb(:)
232
233  ENDDO
234
235  ! Update saved mixing ratios
236  pqrold(ip,:) = q(:)
237
238END SUBROUTINE recombin_corrk
239
Note: See TracBrowser for help on using the repository browser.