1! ======================================================================================================
2! This kernel represents a distillation of *part* of
3! the taumol03 calculation in the gas optics part of the PSRAD
4! atmospheric
5! radiation code.
6!
7! It is meant to show conceptually how one might "SIMD-ize" swaths of
8! the taumol03 code related to calculating the
9! taug term, so that the impact of the conditional expression on
10! specparm could be reduced and at least partial vectorization
11! across columns could be achieved.
12!
13! I consider it at this point to be "compiling pseudo-code".
14!
15! By this I mean that the code as written compiles under ifort, but has
16! not been tested
17! for correctness, nor I have written a driver routine for it. It does
18! not contain everything
19! that is going on in the taug parent taumol03 code, but I don't claim
20! to actually completely
21! understand the physical meaning of all or even most of the inputs
22! required to make it run.
23!
24! It has been written to vectorize, but apparently does not actually do
25! that
26! under the ifort V13 compiler with the -xHost -O3 level of
27! optimization, even with !dir$ assume_aligned directives.
28! I hypothesize that the compiler is baulking to do so for the indirect
29! addressed calls into the absa
30! look-up table, either that or 4 byte integers may be causing alignment
31! issues relative to 8 byte reals. Anyway,
32! it seems to complain about the key loop being too complex.
33! ======================================================================================================
34MODULE mo_taumol03
35    USE mo_kind, only:wp
36    USE mo_lrtm_setup, ONLY: nspa
37    USE mo_lrtm_setup, ONLY: nspb
38    USE rrlw_planck, ONLY: chi_mls
39    USE rrlw_kg03, ONLY: selfref
40    USE rrlw_kg03, ONLY: forref
41    USE rrlw_kg03, ONLY: ka_mn2o
42    USE rrlw_kg03, ONLY: absa
43    USE rrlw_kg03, ONLY: fracrefa
44    USE rrlw_kg03, ONLY: kb_mn2o
45    USE rrlw_kg03, ONLY: absb
46    USE rrlw_kg03, ONLY: fracrefb
47    IMPLICIT NONE
48    PRIVATE
49    PUBLIC taumol03_lwr,taumol03_upr
50    CONTAINS
51    SUBROUTINE taumol03_lwr(ncol, startCol, laytrop, nlayers,         &
52                rat_h2oco2, colco2, colh2o, coln2o, coldry, &
53                fac0, fac1, minorfrac,                      &
54                selffac,selffrac,forfac,forfrac,            &
55                jp, jt, ig, indself, &
56                indfor, indminor, &
57                taug, fracs)
58        IMPLICIT NONE
59
60        real(kind=wp), PARAMETER :: oneminus = 1.0_wp - 1.0e-06_wp
61        integer, intent(in) :: ncol                  ! number of simd columns
62        integer, intent(in) :: startCol              !starting index column
63        integer, intent(in) :: laytrop               ! number of layers forwer atmosphere kernel
64        integer, intent(in) :: nlayers               ! total number of layers
65        real(kind=wp), intent(in), dimension(startCol:ncol,0:1,nlayers) :: rat_h2oco2,fac0,fac1    ! not sure of ncol depend
66        real(kind=wp), intent(in), dimension(startCol:ncol,nlayers) :: colco2,colh2o,coln2o,coldry ! these appear to be gas concentrations
67
68        real(kind=wp), intent(in), dimension(startCol:ncol,nlayers) :: selffac,selffrac            ! not sure of ncol depend
69        real(kind=wp), intent(in), dimension(startCol:ncol,nlayers) :: forfac,forfrac              ! not sure of ncol depend
70        real(kind=wp), intent(in), dimension(startCol:ncol,nlayers) :: minorfrac                   ! not sure of ncol depend
71
72        ! Look up tables and related lookup indices
73        ! I assume all lookup indices depend on 3D position
74        ! =================================================
75
76        integer, intent(in) :: jp(startCol:ncol,nlayers)      ! I assume jp depends on ncol
77        integer, intent(in) :: jt(startCol:ncol,0:1,nlayers)  ! likewise for jt
78        integer, intent(in) :: ig                    ! ig indexes into lookup tables
79        integer, intent(in) :: indself(startCol:ncol,nlayers) ! self index array
80        integer, intent(in) :: indfor(startCol:ncol,nlayers)  ! for index array
81        integer, intent(in) :: indminor(startCol:ncol,nlayers) ! ka_mn2o index array
82        real(kind=wp), intent(out), dimension(startCol:ncol,nlayers) :: taug  ! kernel result
83        real(kind=wp), intent(out), dimension(startCol:ncol,nlayers) :: fracs ! kernel result
84
85        ! Local variable
86        ! ==============
87
88        integer :: lay   ! layer index
89        integer :: i     ! specparm types index
90        integer :: icol  ! column index
91        ! vector temporaries
92        ! ====================
93
94        integer, dimension(1:3,1:3) :: caseTypeOperations
95        integer, dimension(startCol:ncol)    :: caseType
96        real(kind=wp), dimension(startCol:ncol) :: p, p4, fs
97        real(kind=wp), dimension(startCol:ncol) :: fmn2o,fmn2omf
98        real(kind=wp), dimension(startCol:ncol) :: fpl
99        real(kind=wp), dimension(startCol:ncol) :: specmult, speccomb, specparm
100        real(kind=wp), dimension(startCol:ncol) :: specmult_mn2o, speccomb_mn2o,specparm_mn2o
101        real(kind=wp), dimension(startCol:ncol) :: specmult_planck, speccomb_planck,specparm_planck
102        real(kind=wp), dimension(startCol:ncol) :: n2om1,n2om2,absn2o,adjcoln2o,adjfac,chi_n2o,ratn2o
103        real(kind=wp), dimension(startCol:ncol,0:1) :: tau_major
104        real(kind=wp), dimension(startCol:ncol) :: taufor,tauself
105        integer, dimension(startCol:ncol) :: js, ind0, ind00, ind01, ind02, jmn2o, jpl
106
107        ! Register temporaries
108        ! ====================
109
110        real(kind=wp) :: p2,fk0,fk1,fk2
111        real(kind=wp) :: refrat_planck_a, refrat_planck_b
112        real(kind=wp) :: refrat_m_a, refrat_m_b
113        integer ::       rrpk_counter=0
114
115        !dir$ assume_aligned jp:64
116        !dir$ assume_aligned jt:64
117        !dir$ assume_aligned rat_h2oco2:64
118        !dir$ assume_aligned colco2:64
119        !dir$ assume_aligned colh2o:64
120        !dir$ assume_aligned fac0:64
121        !dir$ assume_aligned fac1:64
122        !dir$ assume_aligned taug:64
123
124        !dir$ assume_aligned p:64
125        !dir$ assume_aligned p4:64
126        !dir$ assume_aligned specmult:64
127        !dir$ assume_aligned speccomb:64
128        !dir$ assume_aligned specparm:64
129        !dir$ assume_aligned specmult_mn2o:64
130        !dir$ assume_aligned speccomb_mn2o:64
131        !dir$ assume_aligned specparm_mn2o:64
132        !dir$ assume_aligned specmult_planck:64
133        !dir$ assume_aligned speccomb_planck:64
134        !dir$ assume_aligned specparm_planck:64
135        !dir$ assume_aligned indself:64
136        !dir$ assume_aligned indfor:64
137        !dir$ assume_aligned indminor:64
138        !dir$ assume_aligned fs:64
139        !dir$ assume_aligned tau_major:64
140
141        !dir$ assume_aligned js:64
142        !dir$ assume_aligned ind0:64
143        !dir$ assume_aligned ind00:64
144        !dir$ assume_aligned ind01:64
145        !dir$ assume_aligned ind02:64
146
147        !dir$ assume_aligned caseTypeOperations:64
148        !dir$ assume_aligned caseType:64
149
150        ! Initialize Case type operations
151        !=================================
152
153        caseTypeOperations(1:3,1) = (/0, 1, 2/)
154        caseTypeOperations(1:3,2) = (/1, 0,-1/)
155        caseTypeOperations(1:3,3) = (/0, 1, 1/)
156
157        ! Minor gas mapping levels:
158        !     lower - n2o, p = 706.272 mbar, t = 278.94 k
159        !     upper - n2o, p = 95.58 mbar, t = 215.7 k
160
161        !  P = 212.725 mb
162        refrat_planck_a = chi_mls(1,9)/chi_mls(2,9)
163
164        !  P = 95.58 mb
165        refrat_planck_b = chi_mls(1,13)/chi_mls(2,13)
166
167        !  P = 706.270mb
168        refrat_m_a = chi_mls(1,3)/chi_mls(2,3)
169
170        !  P = 95.58 mb
171        refrat_m_b = chi_mls(1,13)/chi_mls(2,13)
172
173        ! Lower atmosphere loop
174        ! =====================
175
176        DO lay = 1,laytrop  ! loop over layers
177
178        ! Compute tau_major term
179        ! ======================
180
181            DO i=0,1         ! loop over specparm types
182
183                ! This loop should vectorize
184                ! =============================
185                !dir$ SIMD
186                DO icol=startCol,ncol  ! Vectorizes with dir - 14.0.2
187                    speccomb(icol) = colh2o(icol,lay) + rat_h2oco2(icol,i,lay)*colco2(icol,lay)
188                    specparm(icol) = colh2o(icol,lay)/speccomb(icol)
189                    IF (specparm(icol) .GE. oneminus) specparm(icol) = oneminus
190                    specmult(icol) = 8._wp*(specparm(icol))
191                    js(icol) = 1 + INT(specmult(icol))
192                    fs(icol) = MOD(specmult(icol),1.0_wp)
193                END DO
194
195                ! The only conditional loop
196                ! =========================
197
198                DO icol=startCol,ncol  ! Vectorizes as is 14.0.2
199                    IF (specparm(icol) .LT. 0.125_wp) THEN
200                        caseType(icol)=1
201                        p(icol)  = fs(icol) - 1.0_wp
202                        p2 = p(icol)*p(icol)
203                        p4(icol) = p2*p2
204                    ELSE IF (specparm(icol) .GT. 0.875_wp) THEN
205                        caseType(icol)=2
206                        p(icol)  = -fs(icol)
207                        p2 = p(icol)*p(icol)
208                        p4(icol) = p2*p2
209                    ELSE
210                        caseType(icol)=3
211                        ! SIMD way of writing fk0= 1-fs and fk1 = fs, fk2=zero
212                        ! ===========================================================
213
214                        p4(icol) = 1.0_wp - fs(icol)
215                        p(icol) = -p4(icol)          ! complicated way of getting fk2 = zero for ELSE case
216                    ENDIF
217                END DO
218
219                ! Vector/SIMD index loop calculation
220                ! ==================================
221
222                !dir$ SIMD
223                DO icol=startCol,ncol  ! Vectorizes with dir 14.0.2
224                    ind0(icol) = ((jp(icol,lay)-(1*MOD(i+1,2)))*5+(jt(icol,i,lay)-1))*nspa(3) +js(icol)
225                    ind00(icol) = ind0(icol) + caseTypeOperations(1,caseType(icol))
226                    ind01(icol) = ind0(icol) + caseTypeOperations(2,caseType(icol))
227                    ind02(icol) = ind0(icol) + caseTypeOperations(3,caseType(icol))
228                END DO
229
230                ! What we've been aiming for a nice flop intensive
231                ! SIMD/vectorizable loop!
232                ! 17 flops
233                !
234                ! Albeit at the cost of a couple extra flops for the fk2 term
235                ! and a repeated lookup table access for the fk2 term in the
236                ! the ELSE case when specparm or specparm1 is  (> .125 && < .875)
237                ! ===============================================================
238
239                !dir$ SIMD
240                DO icol=startCol,ncol ! Vectorizes with dir 14.0.2
241
242                    fk0 = p4(icol)
243                    fk1 = 1.0_wp - p(icol) - 2.0_wp*p4(icol)
244                    fk2 = p(icol) + p4(icol)
245                    tau_major(icol,i) = speccomb(icol) * ( &
246                    fac0(icol,i,lay)*(fk0*absa(ind00(icol),ig)  + &
247                    fk1*absa(ind01(icol),ig)  + &
248                    fk2*absa(ind02(icol),ig)) + &
249                    fac1(icol,i,lay)*(fk0*absa(ind00(icol)+9,ig)  + &
250                    fk1*absa(ind01(icol)+9,ig)  + &
251                    fk2*absa(ind02(icol)+9,ig)))
252                END DO
253
254            END DO ! end loop over specparm types for tau_major calculation
255
256            ! Compute taufor and tauself terms:
257            ! Note the use of 1D bilinear interpolation of selfref and forref
258            ! lookup table values
259            ! ===================================================================================
260
261            !dir$ SIMD
262            DO icol=startCol,ncol  ! Vectorizes with dir 14.0.2
263                tauself(icol) = selffac(icol,lay)*(selfref(indself(icol,lay),ig) +&
264                selffrac(icol,lay)*(selfref(indself(icol,lay)+1,ig)- selfref(indself(icol,lay),ig)))
265                taufor(icol) =  forfac(icol,lay)*( forref(indfor(icol,lay),ig) +&
266                forfrac(icol,lay)*( forref(indfor(icol,lay)+1,ig) -forref(indfor(icol,lay),ig)))
267            END DO
268
269            ! Compute absn2o term:
270            ! Note the use of 2D bilinear interpolation ka_mn2o lookup table
271            ! values
272            ! =====================================================================
273
274            !dir$ SIMD
275            DO icol=startCol,ncol !vectorizes with dir 14.0.2
276                speccomb_mn2o(icol) = colh2o(icol,lay) +refrat_m_a*colco2(icol,lay)
277                specparm_mn2o(icol) = colh2o(icol,lay)/speccomb_mn2o(icol)
278            END DO
279
280            do icol=startCol,ncol ! vectorizes as is 14.0.2
281                IF (specparm_mn2o(icol) .GE. oneminus) specparm_mn2o(icol) =oneminus
282            end do
283
284            !dir$ SIMD ! vectorizes with dir 14.0.2
285            DO icol=startCol,ncol
286                specmult_mn2o(icol) = 8.0_wp*specparm_mn2o(icol)
287                jmn2o(icol) = 1 + INT(specmult_mn2o(icol))
288                fmn2o(icol) = MOD(specmult_mn2o(icol),1.0_wp)
289                fmn2omf(icol) = minorfrac(icol,lay)*fmn2o(icol)
290            END DO
291
292            !
293            ! 2D bilinear interpolation
294            ! =========================
295
296            !dir$ SIMD
297            do icol=startCol,ncol ! vectorizes with dir 14.0.2
298                n2om1(icol)   = ka_mn2o(jmn2o(icol),indminor(icol,lay)  ,ig) + &
299                fmn2o(icol)*(ka_mn2o(jmn2o(icol)+1,indminor(icol,lay),ig) - &
300                ka_mn2o(jmn2o(icol),indminor(icol,lay)  ,ig))
301                n2om2(icol)   = ka_mn2o(jmn2o(icol),indminor(icol,lay)+1,ig) + &
302                fmn2o(icol)*(ka_mn2o(jmn2o(icol)+1,indminor(icol,lay)+1,ig)- &
303                ka_mn2o(jmn2o(icol),indminor(icol,lay)+1,ig))
304                absn2o(icol)  = n2om1(icol) + minorfrac(icol,lay)*(n2om2(icol) -n2om1(icol))
305            end do
306
307            !  In atmospheres where the amount of N2O is too great to be
308            !  considered
309            !  a minor species, adjust the column amount of N2O by an empirical
310            !  factor
311            !  to obtain the proper contribution.
312            ! ========================================================================
313
314            !dir$ SIMD
315            do icol=startCol,ncol  ! vectorized with dir 14.0.2
316                chi_n2o(icol) = coln2o(icol,lay)/coldry(icol,lay)
317                ratn2o(icol) = 1.e20*chi_n2o(icol)/chi_mls(4,jp(icol,lay)+1)
318            end do
319
320            do icol=startCol,ncol ! vectorizes as is 14.0.2
321                IF (ratn2o(icol) .GT. 1.5_wp) THEN
322                    adjfac(icol) = 0.5_wp+(ratn2o(icol)-0.5_wp)**0.65_wp
323                    adjcoln2o(icol) =adjfac(icol)*chi_mls(4,jp(icol,lay)+1)*coldry(icol,lay)*1.e-20_wp
324                ELSE
325                    adjcoln2o(icol) = coln2o(icol,lay)
326                ENDIF
327            end do
328
329            ! Compute taug, one of two terms returned by the lower atmosphere
330            ! kernel (the other is fracs)
331            ! This loop could be parallelized over specparm types (i) but might
332            ! produce
333            ! different results for different thread counts
334            ! ===========================================================================================
335
336            !dir$ SIMD  ! DOES NOT VECTORIZE even with SIMD dir 14.0.2
337            DO icol=startCol,ncol
338                taug(icol,lay) = tau_major(icol,0) + tau_major(icol,1) +tauself(icol) + taufor(icol) + adjcoln2o(icol)*absn2o(icol)
339            END DO
340
341            !dir$ SIMD ! vectorizes with dir 14.0.2
342            DO icol=startCol,ncol
343                speccomb_planck(icol) = colh2o(icol,lay)+refrat_planck_a*colco2(icol,lay)
344                specparm_planck(icol) = colh2o(icol,lay)/speccomb_planck(icol)
345            END DO
346
347            DO icol=startCol,ncol ! vectorizes as is 14.0.2
348                IF (specparm_planck(icol) .GE. oneminus) specparm_planck(icol)=oneminus
349            END DO
350
351            !dir$ SIMD
352            DO icol=startCol,ncol !vectorizes with dir 14.0.2
353                specmult_planck(icol) = 8.0_wp*specparm_planck(icol)
354                jpl(icol)= 1 + INT(specmult_planck(icol))
355                fpl(icol) = MOD(specmult_planck(icol),1.0_wp)
356                fracs(icol,lay) = fracrefa(ig,jpl(icol)) + fpl(icol)*(fracrefa(ig,jpl(icol)+1)-fracrefa(ig,jpl(icol)))
357            END DO
358            rrpk_counter=rrpk_counter+1
359        END DO ! end lower atmosphere loop
360    END SUBROUTINE taumol03_lwr
361
362
363    SUBROUTINE taumol03_upr(ncol, startCol, laytrop, nlayers,                     &
364                rat_h2oco2, colco2, colh2o, coln2o, coldry, &
365                fac0, fac1, minorfrac, &
366                forfac,forfrac,        &
367                jp, jt, ig,        &
368                indfor, indminor, &
369                taug, fracs)
370
371        use mo_kind, only : wp
372        USE mo_lrtm_setup, ONLY: nspa
373        USE mo_lrtm_setup, ONLY: nspb
374        USE rrlw_planck, ONLY: chi_mls
375        USE rrlw_kg03, ONLY: selfref
376        USE rrlw_kg03, ONLY: forref
377        USE rrlw_kg03, ONLY: ka_mn2o
378        USE rrlw_kg03, ONLY: absa
379        USE rrlw_kg03, ONLY: fracrefa
380        USE rrlw_kg03, ONLY: kb_mn2o
381        USE rrlw_kg03, ONLY: absb
382        USE rrlw_kg03, ONLY: fracrefb
383
384        IMPLICIT NONE
385
386        real(kind=wp), PARAMETER :: oneminus = 1.0_wp - 1.0e-06_wp
387
388        integer, intent(in) :: ncol                  ! number of simd columns
389        integer, intent(in) :: startCol              ! starting index for iterations in order to support parallelization across architectures
390        integer, intent(in) :: laytrop               ! number of layers for lower atmosphere kernel
391        integer, intent(in) :: nlayers               ! total number of layers
392        real(kind=wp), intent(in), dimension(startCol:ncol,0:1,nlayers) :: rat_h2oco2,fac0,fac1    ! not sure of ncol depend
393        real(kind=wp), intent(in), dimension(startCol:ncol,nlayers) :: colco2,colh2o,coln2o,coldry ! these appear to be gas concentrations
394        real(kind=wp), intent(in), dimension(startCol:ncol,nlayers) :: forfac,forfrac ! not sure of ncol depend
395        real(kind=wp), intent(in), dimension(startCol:ncol,nlayers) :: minorfrac      ! not sure of ncol depend
396
397        ! Look up tables and related lookup indices
398        ! I assume all lookup indices depend on 3D position
399        ! =================================================
400
401        integer, intent(in) :: jp(startCol:ncol,nlayers)      ! I assume jp depends on ncol
402        integer, intent(in) :: jt(startCol:ncol,0:1,nlayers)  ! likewise for jt
403        integer, intent(in) :: ig                    ! ig indexes into lookup tables
404        integer, intent(in) :: indfor(startCol:ncol,nlayers)  ! for index array
405        integer, intent(in) :: indminor(startCol:ncol,nlayers) ! ka_mn2o index array
406        real(kind=wp), intent(out), dimension(startCol:ncol,nlayers) :: taug  ! kernel result
407        real(kind=wp), intent(out), dimension(startCol:ncol,nlayers) :: fracs ! kernel result
408
409        ! Local variable
410        ! ==============
411
412        integer :: lay   ! layer index
413        integer :: i     ! specparm types index
414        integer :: icol  ! column index
415
416        ! vector temporaries
417        ! ====================
418
419        real(kind=wp), dimension(startCol:ncol) :: fs
420        real(kind=wp), dimension(startCol:ncol) :: fmn2o,fmn2omf
421        real(kind=wp), dimension(startCol:ncol) :: fpl
422        real(kind=wp), dimension(startCol:ncol) :: specmult, speccomb, specparm
423        real(kind=wp), dimension(startCol:ncol) :: specmult_mn2o, speccomb_mn2o, specparm_mn2o
424        real(kind=wp), dimension(startCol:ncol) :: specmult_planck, speccomb_planck,specparm_planck
425        real(kind=wp), dimension(startCol:ncol) :: n2om1,n2om2,absn2o,adjcoln2o,adjfac,chi_n2o,ratn2o
426        real(kind=wp), dimension(startCol:ncol,0:1) :: tau_major
427        real(kind=wp), dimension(startCol:ncol) :: taufor,tauself
428        integer, dimension(startCol:ncol) :: js, ind0, jmn2o, jpl
429
430        ! Register temporaries
431        ! ====================
432
433        real(kind=wp) :: p2,fk0,fk1,fk2
434        real(kind=wp) :: refrat_planck_a, refrat_planck_b
435        real(kind=wp) :: refrat_m_a, refrat_m_b
436
437        !dir$ assume_aligned jp:64
438        !dir$ assume_aligned jt:64
439        !dir$ assume_aligned rat_h2oco2:64
440        !dir$ assume_aligned colco2:64
441        !dir$ assume_aligned colh2o:64
442        !dir$ assume_aligned fac0:64
443        !dir$ assume_aligned fac1:64
444        !dir$ assume_aligned taug:64
445
446        !dir$ assume_aligned specmult:64
447        !dir$ assume_aligned speccomb:64
448        !dir$ assume_aligned specparm:64
449        !dir$ assume_aligned specmult_mn2o:64
450        !dir$ assume_aligned speccomb_mn2o:64
451        !dir$ assume_aligned specparm_mn2o:64
452        !dir$ assume_aligned specmult_planck:64
453        !dir$ assume_aligned speccomb_planck:64
454        !dir$ assume_aligned specparm_planck:64
455        !dir$ assume_aligned fs:64
456        !dir$ assume_aligned tau_major:64
457        !dir$ assume_aligned chi_n2o:64
458
459        !dir$ assume_aligned js:64
460        !dir$ assume_aligned ind0:64
461        !dir$ assume_aligned jpl:64
462        !dir$ assume_aligned fpl:64
463
464        !dir$ assume_aligned absn2o:64
465        !dir$ assume_aligned adjcoln2o:64
466        !dir$ assume_aligned adjfac:64
467        !dir$ assume_aligned ratn2o:64
468
469        ! Upper atmosphere loop
470        ! ========================
471        refrat_planck_b = chi_mls(1,13)/chi_mls(2,13)
472        refrat_m_b = chi_mls(1,13)/chi_mls(2,13)
473        DO lay = laytrop+1, nlayers
474
475            DO i=0,1         ! loop over specparm types
476
477            ! This loop should vectorize
478            ! =============================
479
480            !dir$ SIMD
481                DO icol=startCol,ncol ! Vectorizes with dir 14.0.2
482                    speccomb(icol) = colh2o(icol,lay) + rat_h2oco2(icol,i,lay)*colco2(icol,lay)
483                    specparm(icol) = colh2o(icol,lay)/speccomb(icol)
484                    IF (specparm(icol) .ge. oneminus) specparm(icol) = oneminus
485                    specmult(icol) = 4.0_wp*(specparm(icol))
486                    js(icol) = 1 + INT(specmult(icol))
487                    fs(icol) = MOD(specmult(icol),1.0_wp)
488                    ind0(icol) = ((jp(icol,lay)-13+i)*5+(jt(icol,i,lay)-1))*nspb(3) +js(icol)
489                END DO
490
491                !dir$ SIMD
492                DO icol=startCol,ncol ! Vectorizes with dir 14.0.2
493                    tau_major(icol,i) = speccomb(icol) * &
494                    ((1.0_wp - fs(icol))*fac0(icol,i,lay)*absb(ind0(icol)  ,ig) + &
495                    fs(icol) *fac0(icol,i,lay)*absb(ind0(icol)+1,ig) + &
496                    (1.0_wp - fs(icol))*fac1(icol,i,lay)*absb(ind0(icol)+5,ig) + &
497                    fs(icol) *fac1(icol,i,lay)*absb(ind0(icol)+6,ig))
498                END DO
499
500            END DO ! end loop over specparm types for tau_major calculation
501
502            ! Compute taufor terms
503            ! Note the use of 1D bilinear interpolation of selfref and forref lookup
504            ! table values
505            ! ===================================================================================
506            !dir$ SIMD
507            DO icol=startCol,ncol ! Vectorizes with dir 14.0.2
508                taufor(icol)  =  forfac(icol,lay)*( forref(indfor(icol,lay),ig) +  &
509                forfrac(icol,lay)*( forref(indfor(icol,lay)+1,ig) - forref(indfor(icol,lay),ig)))
510            END DO
511
512            ! Compute absn2o term:
513            ! Note the use of 2D bilinear interpolation ka_mn2o lookup table values
514            ! =====================================================================
515            !$DIR SIMD
516            DO icol=startCol,ncol ! Vectorizes with dir 14.0.2
517                speccomb_mn2o(icol) = colh2o(icol,lay) + refrat_m_b*colco2(icol,lay)
518                specparm_mn2o(icol) = colh2o(icol,lay)/speccomb_mn2o(icol)
519                IF (specparm_mn2o(icol) .GE. oneminus) specparm_mn2o(icol) = oneminus
520                specmult_mn2o(icol) = 4.0_wp*specparm_mn2o(icol)
521                jmn2o(icol) = 1 + INT(specmult_mn2o(icol))
522                fmn2o(icol) = MOD(specmult_mn2o(icol),1.0_wp)
523                fmn2omf(icol) = minorfrac(icol,lay)*fmn2o(icol)
524            END DO
525
526            !  In atmospheres where the amount of N2O is too great to be considered
527            !  a minor species, adjust the column amount of N2O by an empirical factor
528            !  to obtain the proper contribution.
529            ! ========================================================================
530
531            !dir$ SIMD
532            DO icol=startCol,ncol ! loop vectorized with directive 14.0.2
533                chi_n2o(icol) = coln2o(icol,lay)/coldry(icol,lay)
534                ratn2o(icol) = 1.e20*chi_n2o(icol)/chi_mls(4,jp(icol,lay)+1)
535            END DO
536
537            DO icol=startCol,ncol ! Loop vectorized as is 14.0.2
538                IF (ratn2o(icol) .GT. 1.5_wp) THEN
539                    adjfac(icol) = 0.5_wp+(ratn2o(icol)-0.5_wp)**0.65_wp
540                    adjcoln2o(icol) = adjfac(icol)*chi_mls(4,jp(icol,lay)+1)*coldry(icol,lay)*1.e-20_wp
541                ELSE
542                    adjcoln2o(icol) = coln2o(icol,lay)
543                ENDIF
544            END DO
545
546            !
547            ! 2D bilinear interpolation
548            ! =========================
549
550            !dir$ SIMD
551            DO icol=startCol,ncol ! loop vectorizes with directive 14.0.2
552                n2om1(icol)   = kb_mn2o(jmn2o(icol),indminor(icol,lay)  ,ig) + &
553                fmn2o(icol)*(kb_mn2o(jmn2o(icol)+1,indminor(icol,lay),ig) - &
554                kb_mn2o(jmn2o(icol)  ,indminor(icol,lay)  ,ig))
555                n2om2(icol)   = kb_mn2o(jmn2o(icol),indminor(icol,lay)+1,ig) + &
556                fmn2o(icol)*(kb_mn2o(jmn2o(icol)+1,indminor(icol,lay)+1,ig)- &
557                kb_mn2o(jmn2o(icol),indminor(icol,lay)+1,ig))
558                absn2o(icol)  = n2om1(icol) + minorfrac(icol,lay)*(n2om2(icol) -n2om1(icol))
559            END DO
560
561            !dir$ SIMD
562            DO icol=startCol,ncol ! loop vectorizes with directive 14.0.2
563                taug(icol,lay) =  tau_major(icol,0) + tau_major(icol,1) + taufor(icol) + adjcoln2o(icol)*absn2o(icol)
564            END DO
565
566            !dir$ SIMD
567            DO icol=startCol,ncol ! loop vectorizes with directive 14.0.2
568                speccomb_planck(icol) = colh2o(icol,lay)+refrat_planck_b*colco2(icol,lay)
569                specparm_planck(icol) = colh2o(icol,lay)/speccomb_planck(icol)
570            END DO
571
572            !dir$ SIMD
573            DO icol=startCol,ncol  ! loop vectorizes with directive 14.0.2
574                IF (specparm_planck(icol) .GE. oneminus) specparm_planck(icol)=oneminus
575                specmult_planck(icol) = 4.0_wp*specparm_planck(icol)
576                jpl(icol)= 1 + INT(specmult_planck(icol))
577                fpl(icol) = MOD(specmult_planck(icol),1.0_wp)
578                fracs(icol,lay) =  fracrefb(ig,jpl(icol)) + fpl(icol)*(fracrefb(ig,jpl(icol)+1)-fracrefb(ig,jpl(icol)))
579            END DO
580        END DO  ! nlayers loop
581
582    END SUBROUTINE taumol03_upr
583
584END MODULE mo_taumol03
585