Source code for exo_k.two_stream.two_stream_crisp

# -*- coding: utf-8 -*-
"""
Created in Jan 2021

@author: jeremy leconte
"""
import numpy as np
import numba

jit=False

[docs] @numba.jit(nopython=True, fastmath=True, cache=True) def solve_2stream_nu_xsec(source_nu, dtau_nu, omega0_nu, g_asym_nu, flux_top_dw_nu, mu0=2./3., alb_surf=0., verbose=False, flux_at_level=False): """Deals with the spectral axis """ NLEV, NW = source_nu.shape flux_up=np.zeros((NLEV, NW)) flux_dw=np.zeros((NLEV, NW)) flux_net=np.zeros((NLEV, NW)) kernel=np.zeros((NLEV, NLEV, NW)) for iW in range(NW): flux_up[:,iW], flux_dw[:,iW], flux_net[:,iW], kernel[:,:,iW] = \ solve_2stream(source_nu[:,iW], dtau_nu[:,iW], omega0_nu[:,iW], g_asym_nu[:,iW], mu0=mu0, flux_top_dw=flux_top_dw_nu[iW], alb_surf=alb_surf, verbose=verbose) return flux_up, flux_dw, flux_up-flux_dw, kernel
[docs] @numba.jit(nopython=True, fastmath=True, cache=True) def solve_2stream_nu_corrk(source_nu, dtau_nu, omega0_nu, g_asym_nu, flux_top_dw_nu, mu0=2./3., alb_surf=0., verbose=False, flux_at_level=False): """Deals with the spectral axis """ NLEV, NW = source_nu.shape NG = dtau_nu.shape[-1] flux_up=np.zeros((NLEV, NW, NG)) flux_dw=np.zeros((NLEV, NW, NG)) flux_net=np.zeros((NLEV, NW, NG)) kernel=np.zeros((NLEV, NLEV, NW, NG)) for iW in range(NW): for iG in range(NG): flux_up[:,iW, iG], flux_dw[:,iW,iG], flux_net[:,iW,iG], kernel[:,:,iW,iG] = \ solve_2stream(source_nu[:,iW], dtau_nu[:,iW,iG], omega0_nu[:,iW,iG], g_asym_nu[:,iW,iG], mu0=mu0, flux_top_dw=flux_top_dw_nu[iW], alb_surf=alb_surf, verbose=verbose) return flux_up, flux_dw, flux_up-flux_dw, kernel
[docs] @numba.jit(nopython=True, fastmath=True, cache=True) def solve_2stream(source, dtau, omega0, g_asym, mu0=2./3., flux_top_dw=0., alb_surf=0., verbose=False): """ Inherited from ExoRem (Blain et al. 2020): https://gitlab.obspm.fr/dblain/exorem/-/blob/master/src/fortran/exorem/radiative_transfer.f90 PURPOSE: THIS subroutine COMPUTES THE UPWARD, DOWNWARD AND NET THERMAL FLUX IN AN INHOMOGENEOUS ABSORBING, SCATTERING ATMOSPHERE. THE TWO_STREAM APPROXIMATION (QUADRATURE WITH COS OF AVERAGE ANGLE = mu0) IS USED TO FIND THE DIFFUSE REFLECTIVITY AND TRANSMISSIVITY AND THE TOTAL UPWARD AND DOWNWARD FLUXES FOR EACH OF THE NLAY HOMOGENEOUS LAYERS. THE ADDING METHOD IS then USED TO COMBINE THESE LAYERS. IF ANY LAYER IS THICKER THAN DTAU = *EMAX1*, IT IS ASSUMED TO BE SEMI-INFINITE. LAYERS THICKER THAN DTAU = *EMAX2*, ARE TREATED AS INFINITE LAYERS. NOTE: TO ACCOUNT FOR A DIFFUSE FLUX AT THE TOP OF THE ATMOS- PHERE, THE USER MUST SET flux_dw(1) EQUAL TO THAT VALUE. FOR NO DOWNWARD DIFFUSE FLUX AT THE TOP OF THE ATMOSPHERE, THE USER MUST INITIALIZE flux_dw(1) TO ZERO IN THE CALLING PROGRAM. Parameters ---------- source: array, np.ndarray PLANCK function AT EACH LEVEL FROM TOP OF THE ATMOSPHERE (L=0) TO THE SURFACE (L=NLAY) (NLAY+1 VALUES) JL: actually seems to be pi times the planck function as defined in the rest of the code. dtau: array, np.ndarray ARRAY OF NORMAL INCIDENCE OPTICAL DEPTHS IN EACH HOMOGENEOUS MODEL LAYER. (NLAY VALUES) omega0: array, np.ndarray ARRAY OF SINGLE SCATTERING ALBEDOS FOR EACH HOMO- GENEOUS MODEL LAYER. (NLAY VALUES) g_asym: array, np.ndarray ARRAY OF ASSYMETRY parameterS FOR EACH HOMOGENEOUS MODEL LAYER. (NLAY VALUES) Returns ------- flux_up: array, np.ndarray UPWARD FLUX AT NLAY+1 LAYER BOUNDARIES. (flux_up(L) REFERS TO THE UPWARD FLUX AT THE TOP OF LAYER L) flux_dw: array, np.ndarray DOWNWARD FLUX AT NLAY+1 LAYER BOUNDARIES. (flux_dw(L) REFERS TO THE DOWNWARD FLUX AT THE BOTTOM OF LAYER L-1) flux_up-flux_dw: array, np.ndarray Net flux at the same levels Other Parameters ---------------- mu0: float Cos of quadrature angle (Original mu0 was 2/3.) flux_top_dw: float Diffuse downward flux at upper interface. """ # se atmosphere, only : n_levels # # implicit none # # # ********** COMMON BLOCKS USED IN DELTA-EDDINGTON ROUTINES. # # integer, intent(in) :: nlay # doubleprecision, intent(in) :: dtau(n_levels - 1), # g_asym(n_levels - 1), source(n_levels) # # doubleprecision, intent(inout) :: omega0(n_levels - 1) # # doubleprecision, intent(out) :: flux_up(n_levels), flux_dw(n_levels), # DKERNEL(n_levels, n_levels) # # doubleprecision :: source2(n_levels), OMP(n_levels - 1), dtauP(n_levels - 1), # RU(n_levels), & # DD(n_levels), DFLUX(n_levels - 1), G1(n_levels - 1), G2(n_levels - 1), & # DU(n_levels - 1), UFL(n_levels), DFL(n_levels), UFLUX(n_levels), # EKT(n_levels - 1), & # SK(n_levels - 1), RL(n_levels), temperatures_layers(n_levels), RS(n_levels) # # integer :: j, l, np2, nlev # doubleprecision :: skt, prec, emax1, emax2, alb, mu0, denom, e2ktm, emis, # flux_top_dw PREC = 1.e-10 EMAX1 = 8. EMAX2 = 24. #EMAX1, EMAX2: optical depth at which layers are treated as semi infinite # and infinite (resp). NLAY = dtau.size NLEV = NLAY + 1 NP2 = NLAY + 2 dtauP=np.zeros_like(dtau) OMP=np.zeros_like(dtau) G1=np.zeros_like(dtau) G2=np.zeros_like(dtau) SK=np.zeros_like(dtau) DU=np.zeros_like(dtau) EKT=np.zeros_like(dtau) DFLUX=np.zeros_like(dtau) RL=np.zeros_like(source) RS=np.zeros_like(source) RU=np.zeros_like(source) DD=np.zeros_like(source) DFL=np.zeros_like(source) UFLUX=np.zeros_like(source) UFL=np.zeros_like(source) source2=np.zeros_like(source) temperatures_layers=np.zeros_like(source) flux_up=np.zeros_like(source) flux_dw=np.zeros_like(source) DKERNEL=np.zeros((NLEV,NLEV)) # #**** SCALE THE OPTICAL DEPTHS, SINGLE SCATTERING ALBEDOS AND THE # SCATTERING ASSYMETRY FACTORS FOR USE IN THE TWO-STREAM # APPROXIMATION. INITIALIZE OTHER QUANTITIES. USE WISCOMBE'S # TRICK TO SUBTRACT A SMALL VALUE FROM THE SINGLE SCATTERING # FOR THE CASE OF A CONSERVATIVE ATMOSPHERE # for L in range(NLAY): if (1. - omega0[L]) < PREC: omega0[L] = 1. - PREC dtauP[L] = (1. - omega0[L] * g_asym[L]) * dtau[L] OMP[L] = (1. - g_asym[L]) * omega0[L] / (1. - omega0[L] * g_asym[L]) G1[L] = (1. - 0.5 * OMP[L]) / mu0 G2[L] = 0.5 * OMP[L] / mu0 SK[L] = np.sqrt(G1[L] * G1[L] - G2[L] * G2[L]) # #**** DIFFUSE TRANSMITTANCE AND REFLECTANCE OF EACH HOMOGENEOUS LAYER # # THE EQUATIONS FOR RL AND temperatures_layers WERE DERIVED BY SOLVING THE HOMOGENEOUS # PART OF THE EQUATION OF TRANSFER (IE. NO SOURCE TERM) # for L in range(NLAY): SKT = SK[L] * dtauP[L] #print(L) if SKT > EMAX2 : # # INFINITE LAYERS # #print('EMAX2', L, SKT, G1[L], SK[L], (G1[L] + SK[L])) RL[L] = G2[L] / (G1[L] + SK[L]) temperatures_layers[L] = 0.e0 continue EKT[L] = np.exp(SKT) if SKT > EMAX1 : # # SEMI-INFINITE LAYERS # #print('EMAX1', L, SKT, (G1[L] + SK[L]),(EKT[L] * (G1[L] + SK[L]))) RL[L] = G2[L] / (G1[L] + SK[L]) temperatures_layers[L] = 2.e0 * SK[L] / (EKT[L] * (G1[L] + SK[L])) continue E2KTM = EKT[L] * EKT[L] - 1. DENOM = G1[L] * E2KTM + SK[L] * (E2KTM + 2.e0) RL[L] = G2[L] * E2KTM / DENOM temperatures_layers[L] = 2.e0 * SK[L] * EKT[L] / DENOM # #**** SET THE "REFLECTIVITY", "TRANSMISSIVITY" AND "EMISSIVITY" OF # THE "SURFACE" ASSUMING SEMI-INFINITE LAYER WITH SAME OMP # AS LAYER NLAY. FOR "EMISSIVITY", ASSUME SAME dB/dTau AS LAYER NLAY #ALB = G2[NLAY-1] / (G1[NLAY-1] + SK[NLAY-1]) #RL[NLEV-1] = (1. - ALB) #EMIS = 1. - ALB + (1. + ALB) * mu0 * \ # (1. - source[NLEV - 2] / source[NLEV-1]) / dtauP[NLAY-1] ##JL21 below, I put a new boundary condition at the bottom to mimick a dark surface. ALB=alb_surf RL[NLEV-1]=1. EMIS=1.-ALB temperatures_layers[NLEV-1] = 0.e0 # #**** USE ADDING METHOD TO FIND THE REFLECTANCE AND TRANSMITTANCE # OF COMBINED LAYERS. ADD DOWNWARD FROM THE TOP AND UPWARD # FROM THE BOTTOM AT THE SAME TIME. # RS[0] = RL[0] RU[NLEV-1] = ALB for L in range(NLAY): DD[L] = 1. / (1. - RS[L] * RL[L + 1]) RS[L + 1] = RL[L + 1] + temperatures_layers[L + 1] * \ temperatures_layers[L + 1] * RS[L] * DD[L] DU[NLEV - L - 2] = 1. / (1. - RL[NLEV - L - 2] * RU[NP2 - L - 2]) RU[NLEV - L - 2] = RL[NLEV - L - 2] + temperatures_layers[NLEV - L - 2] * \ temperatures_layers[NLEV - L - 2] * RU[NP2 - L - 2] * DU[NLEV - L - 2] # #**** COMPUTE THE UPWARD AND DOWNWARD FLUX FOR EACH HOMOGENEOUS LAYER # #**** LOOP OVER J FOR KERNEL for J in range(NLEV + 1): DFL[NLEV-1] = 0.e0 if J <= NLEV-1: flux_dw[0] = 0.e0 for L in range(NLEV): DKERNEL[J, L] = 0.e0 source2[L] = 0.e0 source2[J] = 1. else: flux_dw[0] = flux_top_dw for L in range(NLEV): source2[L] = source[L] UFL[NLEV-1] = EMIS * source2[NLEV-1] for L in range(NLAY): UFL[L], DFL[L]=_DEDIR1(source2[L], source2[L + 1], omega0[L], g_asym[L], dtau[L], mu0) # #**** USE ADDING METHOD TO FIND UPWARD AND DOWNWARD FLUXES # FOR COMBINED LAYERS. START AT TOP # DFLUX[0] = DFL[0] + temperatures_layers[0] * flux_dw[0] for L in range(NLAY-1): DFLUX[L + 1] = temperatures_layers[L + 1] * \ (RS[L] * (RL[L + 1] * DFLUX[L] + UFL[L + 1]) * DD[L] + DFLUX[L]) + DFL[L + 1] flux_dw[NLEV-1] = (DFLUX[NLAY-1] + RS[NLAY-1] * EMIS * source2[NLEV-1]) / \ (1. - RS[NLAY-1] * ALB) # #**** USE ADDING METHOD TO FIND UPWARD AND DOWNWARD FLUXES # FOR COMBINED LAYERS. START AT BOTTOM. # UFLUX[NLEV-1] = UFL[NLEV-1] for L in range(NLAY): UFLUX[NLEV - L - 2] = temperatures_layers[NLEV - L - 2] * (UFLUX[NP2 - L - 2] + \ RU[NP2 - L - 2] * DFL[NLEV - L - 2]) * DU[NLEV - L - 2] + \ UFL[NLEV - L - 2] # #**** FIND THE TOTAL UPWARD AND DOWNWARD FLUXES AT INTERFACES # BETWEEN INHOMOGENEOUS LAYERS. # flux_up[0] = UFLUX[0] + RU[0] * flux_dw[0] flux_up[NLEV-1] = ALB * flux_dw[NLEV-1] + EMIS * source2[NLEV-1] for L in range(NLAY-1): flux_up[NLEV - L - 2] = (UFLUX[NLEV - L - 2] + RU[NLEV - L - 2] * \ DFLUX[NLAY - L - 2]) / (1. - RU[NLEV - L - 2] * \ RS[NLAY - L - 2]) for L in range(NLAY-1): flux_dw[L + 1] = DFLUX[L] + RS[L] * flux_up[L + 1] # if J <= NLEV-1: for L in range(NLEV): DKERNEL[J, L] = flux_up[L] - flux_dw[L] return flux_up, flux_dw, flux_up-flux_dw, DKERNEL
[docs] @numba.jit(nopython=True, fastmath=True, cache=True) def _DEDIR1(BATM1, BATM2, AIR, GIR, TAU, mu0): """ CCCCCCCCCCCCCCCCCCCCCCCCCCC D E D I R 1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CC C PURPOSE : CC C CC C THIS subroutine USES THE TWO-STREAM ROUTINE (QUADRATURE WITH CC C COS OF AVERAGE ANGLE = mu0) TO FIND THE UPWARD AND DOWNWARD CC C THERMAL FLUXES EMITTED FROM A HOMOGENEOUS LAYER. CC C CC C AUTHORS: DAVID PAIGE AND DAVID CRISP CC C CC C INPUT: CC C CC C BTATM1 IS ATMOSPHERIC EMISSION at level L (ARBITRARY UNITS) CC C BTATM2 IS ATMOSPHERIC EMISSION at level L+1 (ARBITRARY UNITS) CC C AIR IS IR SINGLE SCATTERING ALBEDO CC C GIR IS IR ASYMMETRY parameter CC C TAU IS OPTICAL DEPTH (EXTINCTION = ABSORPTION + SCATTERING) CC C CC C OUTPUT : CC C CC C FUPTOP IS UPWARD FLUX AT TOP OF ATM (SAME UNITS AS BTATM1 or 2) CC C FDNBOT IS DOWNWARD FLUX AT SURFACE CC C FUPBOT IS UPWARD FLUX AT SURFACE CC C CC CCCCCCCCCCCCCCCCCCCCCCCCCCC D E D I R 1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCC """ # doubleprecision, intent(in) :: BATM1, BATM2, air, gir, tau # doubleprecision, intent(out) :: fuptop, fdnbot # # doubleprecision :: mu0, c1, c2, cap, dair, db, dtauIR, emkt, \ # epkt, omttp, opttp, v1, v2, v3, v4, v5, v6 DAIR = AIR * (1. - GIR) / (1. - GIR * AIR) CAP = np.sqrt(1. - DAIR) / mu0 OPTTP = np.sqrt(1. - 0.5 * DAIR + np.sqrt(1. - DAIR)) OMTTP = np.sqrt(1. - 0.5 * DAIR - np.sqrt(1. - DAIR)) dtauIR = (1. - GIR * AIR) * TAU db = (BATM2 - BATM1) * mu0 / dtauIR if CAP * dtauIR < 24.: # #**** THE LAYER THICKNESS IS FINITE # EPKT = np.exp(+CAP * dtauIR) EMKT = 1. / EPKT # #**** SET TOP B.C. # V1 = OPTTP V2 = OMTTP V3 = -(BATM1 - db) # #**** SET LOWER B.C. # V4 = EMKT * OMTTP V5 = EPKT * OPTTP V6 = -(BATM2 + db) # #**** SOLVE SYSTEM OF EQUATIONS # C1 = (V3 * V5 - V6 * V2) / (V1 * V5 - V4 * V2) C2 = (V1 * V6 - V3 * V4) / (V1 * V5 - V4 * V2) FUPTOP = C1 * OMTTP + C2 * OPTTP + BATM1 + db FDNBOT = C1 * EMKT * OPTTP + C2 * EPKT * OMTTP + BATM2 - db else: # #**** ASSUME LAYER IS SEMI-INFINITE. # FUPTOP = BATM1 * (1. - OMTTP / OPTTP) + db * (1. + OMTTP / OPTTP) FDNBOT = BATM2 * (1. - OMTTP / OPTTP) - db * (1. + OMTTP / OPTTP) return FUPTOP, FDNBOT