Changeset 3083 for trunk/LMDZ.TITAN/libf/phytitan/optcv.F90
- Timestamp:
- Oct 12, 2023, 10:30:22 AM (14 months ago)
- File:
-
- 1 edited
-
trunk/LMDZ.TITAN/libf/phytitan/optcv.F90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.TITAN/libf/phytitan/optcv.F90
r2366 r3083 1 SUBROUTINE OPTCV(PQMO,NLAY,PLEV,TMID,PMID, & 2 DTAUV,TAUV,TAUCUMV,WBARV,COSBV,TAURAY,TAUGSURF,SEASHAZEFACT) 1 SUBROUTINE OPTCV(PQMO,NLAY,ZLEV,PLEV,TMID,PMID, & 2 DTAUV,TAUV,TAUCUMV,WBARV,COSBV,TAURAY,TAUGSURF,SEASHAZEFACT,& 3 POPT_HAZE,POPT_CLOUDS,CDCOLUMN) 3 4 4 5 use radinc_h … … 7 8 use gases_h 8 9 use datafile_mod, only: haze_opt_file 9 use comcstfi_mod, only: r 10 use callkeys_mod, only: continuum,graybody,callgasvis,corrk_recombin, & 11 callclouds,callmufi,seashaze,uncoupl_optic_haze 12 use tracer_h, only: nmicro,nice 10 use comcstfi_mod, only: pi,r 11 use callkeys_mod, only: continuum,graybody,callgasvis,corrk_recombin, & 12 callclouds,callmufi,seashaze,uncoupl_optic_haze,& 13 opt4clouds,FHVIS, Fssa 14 use tracer_h, only: nmicro,nice,ices_indx 13 15 14 16 implicit none … … 24 26 ! Adapted from the NASA Ames code by R. Wordsworth (2009) 25 27 ! Clean and adaptation to Titan by J. Vatant d'Ollone (2016-17) 28 ! Clean and correction to Titan by B. de Batz de Trenquelléon (2022) 29 ! New optics added for Titan's clouds by B. de Batz de Trenquelléon (2023) 30 ! The whole routine needs to be redone... 26 31 ! 27 32 !================================================================== … … 47 52 REAL*8, INTENT(IN) :: PQMO(nlay,nmicro) ! Tracers for microphysics optics (X/m2). 48 53 INTEGER, INTENT(IN) :: NLAY ! Number of pressure layers (for pqmo) 54 REAL*8, INTENT(IN) :: ZLEV(NLAY+1) 49 55 REAL*8, INTENT(IN) :: PLEV(L_LEVELS) 50 56 REAL*8, INTENT(IN) :: TMID(L_LEVELS), PMID(L_LEVELS) 51 57 REAL*8, INTENT(IN) :: TAURAY(L_NSPECTV) 52 58 REAL*8, INTENT(IN) :: SEASHAZEFACT(L_LEVELS) 59 INTEGER, INTENT(IN) :: CDCOLUMN 53 60 54 61 REAL*8, INTENT(OUT) :: DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS) … … 58 65 REAL*8, INTENT(OUT) :: WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS) 59 66 REAL*8, INTENT(OUT) :: TAUGSURF(L_NSPECTV,L_NGAUSS-1) 67 REAL*8, INTENT(OUT) :: POPT_HAZE(L_LEVELS,L_NSPECTI,3) 68 REAL*8, INTENT(OUT) :: POPT_CLOUDS(L_LEVELS,L_NSPECTI,3) 60 69 ! ========================================================== 61 70 … … 64 73 ! Titan customisation 65 74 ! J. Vatant d'Ollone (2016) 66 real*8 DHAZE_T(L_LEVELS,L_NSPECT I)67 real*8 DHAZES_T(L_LEVELS,L_NSPECT I)68 real*8 SSA_T(L_LEVELS,L_NSPECT I)69 real*8 ASF_T(L_LEVELS,L_NSPECT I)75 real*8 DHAZE_T(L_LEVELS,L_NSPECTV) 76 real*8 DHAZES_T(L_LEVELS,L_NSPECTV) 77 real*8 SSA_T(L_LEVELS,L_NSPECTV) 78 real*8 ASF_T(L_LEVELS,L_NSPECTV) 70 79 ! ========================== 71 80 … … 105 114 real*8 dumwvl 106 115 107 real*8 m3as,m3af 108 real*8 dtauaer_s,dtauaer_f 116 ! Variables for new optics 117 integer iq, iw, FTYPE, CTYPE 118 real*8 m0as,m0af,m0ccn,m3as,m3af,m3ccn,m3ices,m3cld 119 real*8 dtauaer_s,dtauaer_f,dtau_ccn,dtau_cld 109 120 real*8,save :: rhoaer_s(L_NSPECTV),ssa_s(L_NSPECTV),asf_s(L_NSPECTV) 110 121 real*8,save :: rhoaer_f(L_NSPECTV),ssa_f(L_NSPECTV),asf_f(L_NSPECTV) 111 !$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,asf_s,asf_f) 122 real*8,save :: ssa_ccn(L_NSPECTI),asf_ccn(L_NSPECTI) 123 real*8,save :: ssa_cld(L_NSPECTV),asf_cld(L_NSPECTV) 124 !$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,ssa_cld,asf_s,asf_f,asf_cld) 112 125 113 126 logical,save :: firstcall=.true. … … 141 154 ENDIF 142 155 143 !======================================================================= 144 ! Determine the total gas opacity throughout the column, for each 145 ! spectral interval, NW, and each Gauss point, NG. 146 ! Calculate the continuum opacities, i.e., those that do not depend on 147 ! NG, the Gauss index. 148 149 taugsurf(:,:) = 0.0 150 dpr(:) = 0.0 151 lkcoef(:,:) = 0.0 152 153 do K=2,L_LEVELS 154 155 ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed) 156 157 DPR(k) = PLEV(K)-PLEV(K-1) 158 159 ! if we have continuum opacities, we need dz 160 156 !======================================================================= 157 ! Determine the total gas opacity throughout the column, for each 158 ! spectral interval, NW, and each Gauss point, NG. 159 ! Calculate the continuum opacities, i.e., those that do not depend on 160 ! NG, the Gauss index. 161 162 taugsurf(:,:) = 0.0 163 dpr(:) = 0.0 164 lkcoef(:,:) = 0.0 165 166 do K=2,L_LEVELS 167 ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed) 168 DPR(k) = PLEV(K)-PLEV(K-1) 169 170 ! if we have continuum opacities, we need dz 161 171 dz(k) = dpr(k)*R*TMID(K)/(gzlat_ig(ilay)*PMID(K)) 162 172 U(k) = Cmk(ilay)*DPR(k) ! only Cmk line in optcv.F 163 173 164 call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K)) 165 166 do LK=1,4 167 LKCOEF(K,LK) = LCOEF(LK) 168 end do 169 end do ! levels 170 171 ! Rayleigh scattering 172 do NW=1,L_NSPECTV 173 TRAY(1:4,NW) = 1.d-30 174 do K=5,L_LEVELS 175 TRAY(K,NW) = TAURAY(NW) * DPR(K) 176 end do ! levels 177 end do 178 179 ! we ignore K=1... 180 do K=2,L_LEVELS 181 182 ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed) 183 184 do NW=1,L_NSPECTV 185 186 IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN 187 m3as = pqmo(ilay,2) / 2.0 188 m3af = pqmo(ilay,4) / 2.0 189 190 IF ( ilay .lt. 18 ) THEN 191 m3as = pqmo(18,2) / 2.0 * dz(k) / dz(18) 192 m3af = pqmo(18,4) / 2.0 * dz(k) / dz(18) 193 ENDIF 194 195 dtauaer_s = m3as*rhoaer_s(nw) 196 dtauaer_f = m3af*rhoaer_f(nw) 197 DHAZE_T(k,nw) = dtauaer_s + dtauaer_f 198 199 IF ( dtauaer_s + dtauaer_f .GT. 1.D-30 ) THEN 200 SSA_T(k,nw) = ( dtauaer_s*ssa_s(nw) + dtauaer_f*ssa_f(nw) ) / ( dtauaer_s+dtauaer_f ) 201 ASF_T(k,nw) = ( dtauaer_s*ssa_s(nw)*asf_s(nw) + dtauaer_f*ssa_f(nw)*asf_f(nw) ) & 202 / ( ssa_s(nw)*dtauaer_s + ssa_f(nw)*dtauaer_f ) 203 ELSE 204 DHAZE_T(k,nw) = 0.D0 205 SSA_T(k,nw) = 1.0 206 ASF_T(k,nw) = 1.0 207 ENDIF 208 209 IF (callclouds.and.firstcall) & 210 WRITE(*,*) 'WARNING: In optcv, optical properties & 211 &calculations are not implemented yet' 212 ELSE 213 ! Call fixed vertical haze profile of extinction - same for all columns 214 call disr_haze(dz(k),plev(k),wnov(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw)) 215 if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k) 216 ENDIF 217 218 !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR 219 ! but visible does not handle very well diffusion in first layer. 220 ! The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise) 221 ! in the 4 first semilayers in optcv, but not optci. 222 ! This solves random variations of the sw heating at the model top. 223 if (k<5) DHAZE_T(K,:) = 0.0 174 call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K)) 175 176 do LK=1,4 177 LKCOEF(K,LK) = LCOEF(LK) 178 end do 179 end do ! L_LEVELS 180 181 ! Rayleigh scattering 182 do NW=1,L_NSPECTV 183 TRAY(1:4,NW) = 1.d-30 184 do K=5,L_LEVELS 185 TRAY(K,NW) = TAURAY(NW) * DPR(K) 186 end do ! L_LEVELS 187 end do 188 189 do NW=1,L_NSPECTV 190 ! We ignore K=1... 191 do K=2,L_LEVELS 192 ! int. arithmetic => gives the gcm layer index (reversed) 193 ilay = L_NLAYRAD+1 - k/2 224 194 225 DRAYAER = TRAY(K,NW) 226 ! DRAYAER is Tau RAYleigh scattering, plus AERosol opacity 227 DRAYAER = DRAYAER + DHAZE_T(K,NW) ! Titan's aerosol 228 229 DCONT = 0.0 ! continuum absorption 230 231 if(continuum.and.(.not.graybody).and.callgasvis)then 232 ! include continua if necessary 233 wn_cont = dble(wnov(nw)) 234 T_cont = dble(TMID(k)) 235 do igas=1,ngasmx 236 237 p_cont = dble(PMID(k)*scalep*gfrac(igas,ilay)) 238 239 dtemp=0.0 240 if(igas.eq.igas_N2)then 241 242 interm = indv(nw,igas,igas) 243 ! call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm) 244 indv(nw,igas,igas) = interm 245 ! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible 246 247 elseif(igas.eq.igas_H2)then 248 249 ! first do self-induced absorption 250 interm = indv(nw,igas,igas) 251 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm) 252 indv(nw,igas,igas) = interm 253 254 ! then cross-interactions with other gases 255 do jgas=1,ngasmx 256 p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay)) 257 dtempc = 0.0 258 if(jgas.eq.igas_N2)then 259 interm = indv(nw,igas,jgas) 260 call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) 261 indv(nw,igas,jgas) = interm 262 ! should be irrelevant in the visible 263 endif 264 dtemp = dtemp + dtempc 265 enddo 266 267 elseif(igas.eq.igas_CH4)then 268 269 ! first do self-induced absorption 270 interm = indv(nw,igas,igas) 271 call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm) 272 indv(nw,igas,igas) = interm 273 274 ! then cross-interactions with other gases 275 do jgas=1,ngasmx 276 p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay)) 277 dtempc = 0.0 278 if(jgas.eq.igas_N2)then 279 interm = indv(nw,igas,jgas) 280 call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) 281 indv(nw,igas,jgas) = interm 282 endif 283 dtemp = dtemp + dtempc 284 enddo 285 286 endif 287 288 DCONT = DCONT + dtemp 289 290 enddo 291 292 DCONT = DCONT*dz(k) 293 294 endif 295 296 do ng=1,L_NGAUSS-1 297 298 ! Now compute TAUGAS 299 300 ! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically 301 ! the execution time of optci/v -> ~ factor 2 on the whole radiative 302 ! transfer on the tested simulations ! 303 304 if (corrk_recombin) then 305 tmpk = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NW,NG) 306 else 307 tmpk = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG) 308 endif 309 310 KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG) 311 KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG) 312 KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG) 313 KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG) 314 315 ! Interpolate the gaseous k-coefficients to the requested T,P values 316 317 ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) + & 318 LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4) 319 320 321 TAUGAS = U(k)*ANS 322 323 TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT 324 DTAUKV(K,nw,ng) = TAUGAS & 325 + DRAYAER & ! DRAYAER includes all scattering contributions 326 + DCONT ! For parameterized continuum aborption 327 328 end do 329 330 ! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS), 331 ! which holds continuum opacity only 332 333 NG = L_NGAUSS 334 DTAUKV(K,nw,ng) = DRAYAER + DCONT ! Scattering + parameterized continuum absorption, including Titan's haze 335 336 end do 337 end do 338 195 ! Optics coupled with the microphysics : 196 IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN 197 198 !========================================================================================== 199 ! Old optics (must have callclouds = .false.): 200 !========================================================================================== 201 IF (.NOT. opt4clouds) THEN 202 m3as = pqmo(ilay,2) / 2.0 203 m3af = pqmo(ilay,4) / 2.0 204 ! Cut-off 205 IF (ilay .lt. 19) THEN 206 m3as = pqmo(19,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19)) 207 m3af = pqmo(19,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19)) 208 ENDIF 209 210 dtauaer_s = m3as*rhoaer_s(nw) 211 dtauaer_f = m3af*rhoaer_f(nw) 212 213 !========================================================================================== 214 ! New optics : 215 !========================================================================================== 216 ELSE 217 iw = (L_NSPECTV + 1) - NW ! Visible first and return 218 !----------------------------- 219 ! HAZE (Spherical + Fractal) : 220 !----------------------------- 221 FTYPE = 1 222 223 ! Spherical aerosols : 224 !--------------------- 225 CTYPE = 5 226 m0as = pqmo(ilay,1) / 2.0 227 m3as = pqmo(ilay,2) / 2.0 228 ! If not callclouds : must have a cut-off 229 IF (.NOT. callclouds) THEN 230 IF (ilay .lt. 19) THEN 231 m0as = pqmo(19,1) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19)) 232 m3as = pqmo(19,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19)) 233 ENDIF 234 ENDIF 235 call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0as,m3as,iw,dtauaer_s,ssa_s(nw),asf_s(nw)) 236 237 ! Fractal aerosols : 238 !------------------- 239 CTYPE = FTYPE 240 m0af = pqmo(ilay,3) / 2.0 241 m3af = pqmo(ilay,4) / 2.0 242 ! If not callclouds : must have a cut-off 243 IF (.NOT. callclouds) THEN 244 IF (ilay .lt. 19) THEN 245 m0af = pqmo(19,3) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19)) 246 m3af = pqmo(19,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19)) 247 ENDIF 248 ENDIF 249 call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0af,m3af,iw,dtauaer_f,ssa_f(nw),asf_f(nw)) 250 ENDIF 251 252 ! Total of Haze opacity (dtau), SSA (w) and ASF (COS) : 253 DHAZE_T(k,nw) = dtauaer_s + dtauaer_f 254 IF (dtauaer_s + dtauaer_f .GT. 1.D-30) THEN 255 SSA_T(k,nw) = ( dtauaer_s*ssa_s(nw) + dtauaer_f*ssa_f(nw) ) / ( dtauaer_s+dtauaer_f ) 256 ASF_T(k,nw) = ( dtauaer_s*ssa_s(nw)*asf_s(nw) + dtauaer_f*ssa_f(nw)*asf_f(nw) ) & 257 / ( ssa_s(nw)*dtauaer_s + ssa_f(nw)*dtauaer_f ) 258 ELSE 259 DHAZE_T(k,nw) = 0.D0 260 SSA_T(k,nw) = 1.0 261 ASF_T(k,nw) = 1.0 262 ENDIF 263 ! Diagnostics for the haze : 264 POPT_HAZE(k,nw,1) = DHAZE_T(k,nw) ! dtau 265 POPT_HAZE(k,nw,2) = SSA_T(k,nw) ! wbar 266 POPT_HAZE(k,nw,3) = ASF_T(k,nw) ! gbar 267 268 !--------------------- 269 ! CLOUDS (Spherical) : 270 !--------------------- 271 IF (callclouds) THEN 272 CTYPE = 0 273 m0ccn = pqmo(ilay,5) / 2.0 274 m3ccn = pqmo(ilay,6) / 2.0 275 m3ices = 0.0d0 276 m3cld = m3ccn 277 278 ! Clear / Dark column method : 279 !----------------------------- 280 281 ! CCN's SSA : 282 call get_haze_and_cloud_opacity(FTYPE,FTYPE,m0ccn,m3ccn,iw,dtau_ccn,ssa_ccn(nw),asf_ccn(nw)) 283 284 ! Clear column (CCN, C2H2, C2H6, HCN) : 285 IF (CDCOLUMN == 0) THEN 286 DO iq = 2, nice 287 m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0) 288 m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0) 289 ENDDO 290 call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw)) 291 292 ! Dark column (CCN, CH4, C2H2, C2H6, HCN) : 293 ELSEIF (CDCOLUMN == 1) THEN 294 DO iq = 1, nice 295 m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0) 296 m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0) 297 ENDDO 298 call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw)) 299 300 ELSE 301 WRITE(*,*) 'WARNING OPTCV.F90 : CDCOLUMN MUST BE 0 OR 1' 302 WRITE(*,*) 'WE USE DARK COLUMN ...' 303 DO iq = 1, nice 304 m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0) 305 m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0) 306 ENDDO 307 call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw)) 308 ENDIF 309 310 ! For small dropplets, opacity of nucleus dominates 311 !ssa_cld(nw) = (ssa_ccn(nw)*m3ccn + ssa_cld(nw)*m3ices) / (m3ccn + m3ices) 312 ssa_cld(nw) = (ssa_ccn(nw)*m3ccn + ssa_cld(nw)*m3cld) / (m3ccn + m3cld) 313 ssa_cld(nw) = ssa_cld(nw) * Fssa 314 315 ! Total of Haze + Clouds opacity (dtau), SSA (w) and ASF (COS) : 316 DHAZE_T(k,nw) = dtauaer_s + dtauaer_f + dtau_cld 317 IF (DHAZE_T(k,nw) .GT. 1.D-30) THEN 318 SSA_T(k,nw) = ( dtauaer_s*ssa_s(nw) + dtauaer_f*ssa_f(nw) + dtau_cld*ssa_cld(nw) ) / ( dtauaer_s+dtauaer_f+dtau_cld ) 319 ASF_T(k,nw) = ( dtauaer_s*ssa_s(nw)*asf_s(nw) + dtauaer_f*ssa_f(nw)*asf_f(nw) + dtau_cld*ssa_cld(nw)*asf_cld(nw) ) & 320 / ( ssa_s(nw)*dtauaer_s + ssa_f(nw)*dtauaer_f + ssa_cld(nw)*dtau_cld ) 321 ELSE 322 DHAZE_T(k,nw) = 0.D0 323 SSA_T(k,nw) = 1.0 324 ASF_T(k,nw) = 1.0 325 ENDIF 326 327 ! Tuning of optical properties (now useless...) : 328 DHAZE_T(k,nw) = DHAZE_T(k,nw) * (FHVIS * (1-SSA_T(k,nw)) + SSA_T(k,nw)) 329 SSA_T(k,nw) = SSA_T(k,nw) / (FHVIS * (1-SSA_T(k,nw)) + SSA_T(k,nw)) 330 331 ! Diagnostics for clouds : 332 POPT_CLOUDS(k,nw,1) = DHAZE_T(k,nw) ! dtau 333 POPT_CLOUDS(k,nw,2) = SSA_T(k,nw) ! wbar 334 POPT_CLOUDS(k,nw,3) = ASF_T(k,nw) ! gbar 335 336 ELSE 337 ! Diagnostics for clouds : 338 POPT_CLOUDS(k,nw,1) = 0.D0 ! dtau 339 POPT_CLOUDS(k,nw,2) = 1.0 ! wbar 340 POPT_CLOUDS(k,nw,3) = 1.0 ! gbar 341 ENDIF 342 343 ! Optics and microphysics no coupled : 344 ELSE 345 ! Call fixed vertical haze profile of extinction - same for all columns 346 call disr_haze(dz(k),plev(k),wnov(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw)) 347 if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k) 348 ! Diagnostics for the haze : 349 POPT_HAZE(k,nw,1) = DHAZE_T(k,nw) ! dtau 350 POPT_HAZE(k,nw,2) = SSA_T(k,nw) ! wbar 351 POPT_HAZE(k,nw,3) = ASF_T(k,nw) ! gbar 352 ! Diagnostics for clouds : 353 POPT_CLOUDS(k,nw,1) = 0.D0 ! dtau 354 POPT_CLOUDS(k,nw,2) = 1.0 ! wbar 355 POPT_CLOUDS(k,nw,3) = 1.0 ! gbar 356 ENDIF ! ENDIF callmufi 357 358 !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR 359 ! but visible does not handle very well diffusion in first layer. 360 ! The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise) 361 ! in the 4 first semilayers in optcv, but not optci. 362 ! This solves random variations of the sw heating at the model top. 363 if (k<5) DHAZE_T(K,:) = 0.0 364 365 DRAYAER = TRAY(K,NW) 366 ! DRAYAER is Tau RAYleigh scattering, plus AERosol opacity 367 DRAYAER = DRAYAER + DHAZE_T(K,NW) ! Titan's aerosol 368 369 DCONT = 0.0 ! continuum absorption 370 371 if(continuum.and.(.not.graybody).and.callgasvis)then 372 ! include continua if necessary 373 wn_cont = dble(wnov(nw)) 374 T_cont = dble(TMID(k)) 375 do igas=1,ngasmx 376 377 p_cont = dble(PMID(k)*scalep*gfrac(igas,ilay)) 378 379 dtemp=0.0 380 if(igas.eq.igas_N2)then 381 382 interm = indv(nw,igas,igas) 383 ! call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm) 384 indv(nw,igas,igas) = interm 385 ! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible 386 387 elseif(igas.eq.igas_H2)then 388 389 ! first do self-induced absorption 390 interm = indv(nw,igas,igas) 391 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm) 392 indv(nw,igas,igas) = interm 393 394 ! then cross-interactions with other gases 395 do jgas=1,ngasmx 396 p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay)) 397 dtempc = 0.0 398 if(jgas.eq.igas_N2)then 399 interm = indv(nw,igas,jgas) 400 call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) 401 indv(nw,igas,jgas) = interm 402 ! should be irrelevant in the visible 403 endif 404 dtemp = dtemp + dtempc 405 enddo 406 407 elseif(igas.eq.igas_CH4)then 408 409 ! first do self-induced absorption 410 interm = indv(nw,igas,igas) 411 call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm) 412 indv(nw,igas,igas) = interm 413 414 ! then cross-interactions with other gases 415 do jgas=1,ngasmx 416 p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay)) 417 dtempc = 0.0 418 if(jgas.eq.igas_N2)then 419 interm = indv(nw,igas,jgas) 420 call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) 421 indv(nw,igas,jgas) = interm 422 endif 423 dtemp = dtemp + dtempc 424 enddo 425 426 endif 427 428 DCONT = DCONT + dtemp 429 430 enddo 431 432 DCONT = DCONT*dz(k) 433 434 endif 435 436 do ng=1,L_NGAUSS-1 437 438 ! Now compute TAUGAS 439 440 ! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically 441 ! the execution time of optci/v -> ~ factor 2 on the whole radiative 442 ! transfer on the tested simulations ! 443 444 if (corrk_recombin) then 445 tmpk = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NW,NG) 446 else 447 tmpk = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG) 448 endif 449 450 KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG) 451 KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG) 452 KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG) 453 KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG) 454 455 ! Interpolate the gaseous k-coefficients to the requested T,P values 456 457 ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) + & 458 LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4) 459 460 461 TAUGAS = U(k)*ANS 462 463 TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT 464 DTAUKV(K,nw,ng) = TAUGAS & 465 + DRAYAER & ! DRAYAER includes all scattering contributions 466 + DCONT ! For parameterized continuum aborption 467 end do 468 469 ! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS), 470 ! which holds continuum opacity only 471 472 NG = L_NGAUSS 473 DTAUKV(K,nw,ng) = DRAYAER + DCONT ! Scattering + parameterized continuum absorption, including Titan's haze 474 475 end do ! k = L_LEVELS 476 end do ! nw = L_NSPECTI 339 477 340 478 !======================================================================= … … 355 493 ENDDO 356 494 357 495 ! NW spectral loop 358 496 DO NW=1,L_NSPECTV 359 DO L=1,L_NLAYRAD-1 360 K = 2*L+1 361 atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW) + ASF_T(K+1,NW)*DHAZES_T(K+1,NW) 362 btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW) 363 ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) ! JVO 2017 : does this 0.999 is really meaningful ? 364 btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW) 365 COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW) 366 END DO ! L vertical loop 497 ! L vertical loop 498 DO L=1,L_NLAYRAD-1 499 K = 2*L+1 500 atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW) + ASF_T(K+1,NW)*DHAZES_T(K+1,NW) 501 btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW) 502 ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) ! JVO 2017 : does this 0.999 is really meaningful ? 503 btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW) 504 COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW) 505 END DO 367 506 368 ! Last level 369 L = L_NLAYRAD 370 K = 2*L+1 371 atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW) 372 btemp(L,NW) = DHAZES_T(K,NW) 373 ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW) ! JVO 2017 : does this 0.999 is really meaningful ? 374 btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) 375 COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW) 376 377 378 END DO ! NW spectral loop 379 507 ! Last level 508 L = L_NLAYRAD 509 K = 2*L+1 510 atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW) 511 btemp(L,NW) = DHAZES_T(K,NW) 512 ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW) ! JVO 2017 : does this 0.999 is really meaningful ? 513 btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) 514 COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW) 515 END DO 516 517 ! NG Gauss loop 380 518 DO NG=1,L_NGAUSS 381 DO NW=1,L_NSPECTV 382 DO L=1,L_NLAYRAD-1 383 384 K = 2*L+1 385 DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG) 386 WBARV(L,nw,ng) = ctemp(L,NW) / DTAUV(L,nw,ng) 387 388 END DO ! L vertical loop 389 390 ! Last level 391 392 L = L_NLAYRAD 393 K = 2*L+1 394 DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) 395 396 WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG) 397 398 END DO ! NW spectral loop 399 END DO ! NG Gauss loop 519 ! NW spectral loop 520 DO NW=1,L_NSPECTV 521 ! L vertical loop 522 DO L=1,L_NLAYRAD-1 523 K = 2*L+1 524 DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG) 525 WBARV(L,nw,ng) = ctemp(L,NW) / DTAUV(L,nw,ng) 526 END DO 527 528 ! Last level 529 L = L_NLAYRAD 530 K = 2*L+1 531 DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) 532 WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG) 533 END DO 534 END DO 400 535 401 536 ! Total extinction optical depths 402 403 DO NG=1,L_NGAUSS ! full gauss loop 404 DO NW=1,L_NSPECTV 405 TAUCUMV(1,NW,NG)=0.0D0 406 DO K=2,L_LEVELS 407 TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG) 408 END DO 409 410 DO L=1,L_NLAYRAD 411 TAUV(L,NW,NG)=TAUCUMV(2*L,NW,NG) 412 END DO 413 TAUV(L,NW,NG)=TAUCUMV(2*L_NLAYRAD+1,NW,NG) 414 END DO 415 END DO ! end full gauss loop 537 !DO NG=1,L_NGAUSS ! full gauss loop 538 ! DO NW=1,L_NSPECTV 539 ! TAUCUMV(1,NW,NG)=0.0D0 540 ! DO K=2,L_LEVELS 541 ! TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG) 542 ! END DO 543 ! DO L=1,L_NLAYRAD 544 ! TAUV(L,NW,NG)=TAUCUMV(2*L,NW,NG) 545 ! END DO 546 ! TAUV(L,NW,NG)=TAUCUMV(2*L_NLAYRAD+1,NW,NG) 547 ! END DO 548 !END DO ! end full gauss loop 549 550 TAUCUMV(:,:,:) = DTAUKV(:,:,:) 551 DO L=1,L_NLAYRAD 552 TAUV(L,:,:)=TAUCUMV(2*L,:,:) 553 END DO 554 TAUV(L,:,:)=TAUCUMV(2*L_NLAYRAD+1,:,:) 416 555 417 556 if(firstcall) firstcall = .false.
Note: See TracChangeset
for help on using the changeset viewer.
