Ignore:
Timestamp:
Feb 8, 2019, 5:45:36 PM (6 years ago)
Author:
jvatant
Message:

Update radiative transfer with some of the latest updates in the generic model
Cf r2056 by AB and r1987-1991 by JL
--JVO

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.TITAN/libf/phytitan/gfluxi.F

    r1420 r2095  
    11      SUBROUTINE GFLUXI(NLL,TLEV,NW,DW,DTAU,TAUCUM,W0,COSBAR,UBARI,
    22     *                  RSF,BTOP,BSURF,FTOPUP,FMIDP,FMIDM)
    3 
    4 C  THIS SUBROUTINE TAKES THE OPTICAL CONSTANTS AND BOUNDARY CONDITIONS
    5 C  FOR THE INFRARED FLUX AT ONE WAVELENGTH AND SOLVES FOR THE FLUXES AT
    6 C  THE LEVELS. THIS VERSION IS SET UP TO WORK WITH LAYER OPTICAL DEPTHS
    7 C  MEASURED FROM THE TOP OF EACH LAYER.  THE TOP OF EACH LAYER HAS 
    8 C  OPTICAL DEPTH ZERO.  IN THIS SUB LEVEL N IS ABOVE LAYER N. THAT IS LAYER N
    9 C  HAS LEVEL N ON TOP AND LEVEL N+1 ON BOTTOM. OPTICAL DEPTH INCREASES
    10 C  FROM TOP TO BOTTOM. SEE C.P. MCKAY, TGM NOTES.
    11 C THE TRI-DIAGONAL MATRIX SOLVER IS DSOLVER AND IS DOUBLE PRECISION SO MANY
    12 C VARIABLES ARE PASSED AS SINGLE THEN BECOME DOUBLE IN DSOLVER
    13 C
    14 C NLL            = NUMBER OF LEVELS (NLAYERS + 1) MUST BE LESS THAT NL (101)
    15 C TLEV(L_LEVELS) = ARRAY OF TEMPERATURES AT GCM LEVELS
    16 C WAVEN          = WAVELENGTH FOR THE COMPUTATION
    17 C DW             = WAVENUMBER INTERVAL
    18 C DTAU(NLAYER)   = ARRAY OPTICAL DEPTH OF THE LAYERS
    19 C W0(NLEVEL)     = SINGLE SCATTERING ALBEDO
    20 C COSBAR(NLEVEL) = ASYMMETRY FACTORS, 0=ISOTROPIC
    21 C UBARI          = AVERAGE ANGLE, MUST BE EQUAL TO 0.5 IN IR
    22 C RSF            = SURFACE REFLECTANCE
    23 C BTOP           = UPPER BOUNDARY CONDITION ON IR INTENSITY (NOT FLUX)
    24 C BSURF          = SURFACE EMISSION = (1-RSFI)*PLANCK, INTENSITY (NOT FLUX)
    25 C FP(NLEVEL)     = UPWARD FLUX AT LEVELS
    26 C FM(NLEVEL)     = DOWNWARD FLUX AT LEVELS
    27 C FMIDP(NLAYER)  = UPWARD FLUX AT LAYER MIDPOINTS
    28 C FMIDM(NLAYER)  = DOWNWARD FLUX AT LAYER MIDPOINTS
    29 C
    30 C----------------------------------------------------------------------C
    31 
     3     
    324      use radinc_h
    335      use radcommon_h, only: planckir
    346      use comcstfi_mod, only: pi
    35 
    36       implicit none
    37 
    38 
     7     
     8      IMPLICIT NONE
     9     
     10!-----------------------------------------------------------------------
     11!  THIS SUBROUTINE TAKES THE OPTICAL CONSTANTS AND BOUNDARY CONDITIONS
     12!  FOR THE INFRARED FLUX AT ONE WAVELENGTH AND SOLVES FOR THE FLUXES AT
     13!  THE LEVELS.  THIS VERSION IS SET UP TO WORK WITH LAYER OPTICAL DEPTHS
     14!  MEASURED FROM THE TOP OF EACH LAYER.  THE TOP OF EACH LAYER HAS 
     15!  OPTICAL DEPTH ZERO.  IN THIS SUB LEVEL N IS ABOVE LAYER N. THAT IS LAYER N
     16!  HAS LEVEL N ON TOP AND LEVEL N+1 ON BOTTOM. OPTICAL DEPTH INCREASES
     17!  FROM TOP TO BOTTOM.  SEE C.P. MCKAY, TGM NOTES.
     18!  THE TRI-DIAGONAL MATRIX SOLVER IS DSOLVER AND IS DOUBLE PRECISION SO MANY
     19!  VARIABLES ARE PASSED AS SINGLE THEN BECOME DOUBLE IN DSOLVER
     20!
     21! NLL            = NUMBER OF LEVELS (NLAYERS + 1) MUST BE LESS THAT NL (101)
     22! TLEV(L_LEVELS) = ARRAY OF TEMPERATURES AT GCM LEVELS
     23! WAVEN          = WAVELENGTH FOR THE COMPUTATION
     24! DW             = WAVENUMBER INTERVAL
     25! DTAU(NLAYER)   = ARRAY OPTICAL DEPTH OF THE LAYERS
     26! W0(NLEVEL)     = SINGLE SCATTERING ALBEDO
     27! COSBAR(NLEVEL) = ASYMMETRY FACTORS, 0=ISOTROPIC
     28! UBARI          = AVERAGE ANGLE, MUST BE EQUAL TO 0.5 IN IR
     29! RSF            = SURFACE REFLECTANCE
     30! BTOP           = UPPER BOUNDARY CONDITION ON IR INTENSITY (NOT FLUX)
     31! BSURF          = SURFACE EMISSION = (1-RSFI)*PLANCK, INTENSITY (NOT FLUX)
     32! FP(NLEVEL)     = UPWARD FLUX AT LEVELS
     33! FM(NLEVEL)     = DOWNWARD FLUX AT LEVELS
     34! FMIDP(NLAYER)  = UPWARD FLUX AT LAYER MIDPOINTS
     35! FMIDM(NLAYER)  = DOWNWARD FLUX AT LAYER MIDPOINTS
     36!-----------------------------------------------------------------------
     37     
    3938      INTEGER NLL, NLAYER, L, NW, NT, NT2
    4039      REAL*8  TERM, CPMID, CMMID
    … …  
    4645      REAL*8  WAVEN, DW, UBARI, RSF
    4746      REAL*8  BTOP, BSURF, FMIDP(L_NLAYRAD), FMIDM(L_NLAYRAD)
    48       REAL*8  B0(L_NLAYRAD),B1(L_NLAYRAD),ALPHA(L_NLAYRAD)
     47      REAL*8  B0(L_NLAYRAD)
     48      REAL*8  B1(L_NLAYRAD)
     49      REAL*8  ALPHA(L_NLAYRAD)
    4950      REAL*8  LAMDA(L_NLAYRAD),XK1(L_NLAYRAD),XK2(L_NLAYRAD)
    5051      REAL*8  GAMA(L_NLAYRAD),CP(L_NLAYRAD),CM(L_NLAYRAD)
    5152      REAL*8  CPM1(L_NLAYRAD),CMM1(L_NLAYRAD),E1(L_NLAYRAD)
    52       REAL*8  E2(L_NLAYRAD),E3(L_NLAYRAD),E4(L_NLAYRAD)
    53 
     53      REAL*8  E2(L_NLAYRAD)
     54      REAL*8  E3(L_NLAYRAD)
     55      REAL*8  E4(L_NLAYRAD)
    5456      REAL*8  FTOPUP, FLUXUP, FLUXDN
    55 
    56       real*8 :: TAUMAX = L_TAUMAX
    57 
    58 C======================================================================C
    59  
    60 C     WE GO WITH THE HEMISPHERIC CONSTANT APPROACH IN THE INFRARED
    61  
    62 
     57      REAL*8 :: TAUMAX = L_TAUMAX
     58
     59! AB : variables for interpolation
     60      REAL*8 C1
     61      REAL*8 C2
     62      REAL*8 P1
     63      REAL*8 P2
     64     
     65!=======================================================================
     66!     WE GO WITH THE HEMISPHERIC CONSTANT APPROACH IN THE INFRARED
     67     
    6368      NLAYER = L_NLAYRAD
    6469
    6570      DO L=1,L_NLAYRAD-1
    6671
    67 !----------------------------------------------------
     72!-----------------------------------------------------------------------
    6873! There is a problem when W0 = 1
    6974!         open(888,file='W0')
    … …  
    7277!           endif
    7378! Prevent this with an if statement:
    74         if (W0(L).eq.1.D0) then
    75            W0(L) = 0.99999D0
    76         endif
    77 !----------------------------------------------------
    78 
    79         ALPHA(L) = SQRT( (1.0D0-W0(L))/(1.0D0-W0(L)*COSBAR(L)) )
    80         LAMDA(L) = ALPHA(L)*(1.0D0-W0(L)*COSBAR(L))/UBARI
    81 
    82         NT    = int(TLEV(2*L)*NTfac)   - NTstar+1
    83         NT2   = int(TLEV(2*L+2)*NTfac) - NTstar+1
    84 
    85         B1(L) = (PLANCKIR(NW,NT2)-PLANCKIR(NW,NT))/DTAU(L)
    86         B0(L) = PLANCKIR(NW,NT)
    87       END DO
    88 
    89 C     Take care of special lower layer
    90 
     79!-----------------------------------------------------------------------
     80         if (W0(L).eq.1.D0) then
     81            W0(L) = 0.99999D0
     82         endif
     83         
     84         ALPHA(L) = SQRT( (1.0D0-W0(L))/(1.0D0-W0(L)*COSBAR(L)) )
     85         LAMDA(L) = ALPHA(L)*(1.0D0-W0(L)*COSBAR(L))/UBARI
     86         
     87         NT    = int(TLEV(2*L)*NTfac)   - NTstar+1
     88         NT2   = int(TLEV(2*L+2)*NTfac) - NTstar+1
     89         
     90! AB : PLANCKIR(NW,NT) is replaced by P1, the linear interpolation result for a temperature NT
     91! AB : idem for PLANCKIR(NW,NT2) and P2
     92         C1 = TLEV(2*L) * NTfac - int(TLEV(2*L) * NTfac)
     93         C2 = TLEV(2*L+2)*NTfac - int(TLEV(2*L+2)*NTfac)
     94         P1 = (1.0D0 - C1) * PLANCKIR(NW,NT) + C1 * PLANCKIR(NW,NT+1)
     95         P2 = (1.0D0 - C2) * PLANCKIR(NW,NT2) + C2 * PLANCKIR(NW,NT2+1)
     96         B1(L) = (P2 - P1) / DTAU(L)
     97         B0(L) = P1
     98      END DO
     99     
     100!     Take care of special lower layer
     101     
    91102      L        = L_NLAYRAD
    92103
    … …  
    94105          W0(L) = 0.99999D0
    95106      end if
    96 
     107     
    97108      ALPHA(L) = SQRT( (1.0D0-W0(L))/(1.0D0-W0(L)*COSBAR(L)) )
    98109      LAMDA(L) = ALPHA(L)*(1.0D0-W0(L)*COSBAR(L))/UBARI
    99 
     110     
    100111      ! Tsurf is used for 1st layer source function
    101112      ! -- same results for most thin atmospheres
    … …  
    108119      !! (or this one yields same results
    109120      !NT    = int( (TLEV(2*L)+TLEV(2*L+1))*0.5*NTfac ) - NTstar+1
    110 
     121     
    111122      NT2   = int(TLEV(2*L)*NTfac)   - NTstar+1
    112       B1(L) = (PLANCKIR(NW,NT)-PLANCKIR(NW,NT2))/DTAU(L)
    113       B0(L) = PLANCKIR(NW,NT2)
    114 
     123     
     124! AB : PLANCKIR(NW,NT) is replaced by P1, the linear interpolation result for a temperature NT
     125! AB : idem for PLANCKIR(NW,NT2) and P2
     126      C1 = TLEV(2*L+1)*NTfac - int(TLEV(2*L+1)*NTfac)
     127      C2 = TLEV(2*L) * NTfac - int(TLEV(2*L) * NTfac)
     128      P1 = (1.0D0 - C1) * PLANCKIR(NW,NT) + C1 * PLANCKIR(NW,NT+1)
     129      P2 = (1.0D0 - C2) * PLANCKIR(NW,NT2) + C2 * PLANCKIR(NW,NT2+1)
     130      B1(L) = (P1 - P2) / DTAU(L)
     131      B0(L) = P2
     132     
    115133      DO L=1,L_NLAYRAD
    116         GAMA(L) = (1.0D0-ALPHA(L))/(1.0D0+ALPHA(L))
    117         TERM    = UBARI/(1.0D0-W0(L)*COSBAR(L))
    118 
    119 C       CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
    120 C       BOTTOM OF THE LAYER.  THAT IS AT DTAU OPTICAL DEPTH
    121 
    122         CP(L) = B0(L)+B1(L)*DTAU(L) +B1(L)*TERM
    123         CM(L) = B0(L)+B1(L)*DTAU(L) -B1(L)*TERM
    124 
    125 C       CPM1 AND CMM1 ARE THE CPLUS AND CMINUS TERMS EVALUATED
    126 C       AT THE TOP OF THE LAYER, THAT IS ZERO OPTICAL DEPTH
    127 
    128         CPM1(L) = B0(L)+B1(L)*TERM
    129         CMM1(L) = B0(L)-B1(L)*TERM
    130       END DO
    131  
    132 C     NOW CALCULATE THE EXPONENTIAL TERMS NEEDED
    133 C     FOR THE TRIDIAGONAL ROTATED LAYERED METHOD
    134 C     WARNING IF DTAU(J) IS GREATER THAN ABOUT 35 (VAX)
    135 C     WE CLIP IT TO AVOID OVERFLOW.
    136 
     134         GAMA(L) = (1.0D0-ALPHA(L))/(1.0D0+ALPHA(L))
     135         TERM    = UBARI/(1.0D0-W0(L)*COSBAR(L))
     136         
     137! CPM1 AND CMM1 ARE THE CPLUS AND CMINUS TERMS EVALUATED
     138! AT THE TOP OF THE LAYER, THAT IS ZERO OPTICAL DEPTH
     139         
     140         CPM1(L) = B0(L)+B1(L)*TERM
     141         CMM1(L) = B0(L)-B1(L)*TERM
     142         
     143! CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
     144! BOTTOM OF THE LAYER.  THAT IS AT DTAU OPTICAL DEPTH.
     145! JL18 put CP and CM after the calculation of CPM1 and CMM1 to avoid unecessary calculations.
     146         
     147         CP(L) = CPM1(L) +B1(L)*DTAU(L)
     148         CM(L) = CMM1(L) +B1(L)*DTAU(L)
     149      END DO
     150     
     151! NOW CALCULATE THE EXPONENTIAL TERMS NEEDED
     152! FOR THE TRIDIAGONAL ROTATED LAYERED METHOD
     153! WARNING IF DTAU(J) IS GREATER THAN ABOUT 35 (VAX)
     154! WE CLIP IT TO AVOID OVERFLOW.
     155     
    137156      DO L=1,L_NLAYRAD
    138 
    139 C       CLIP THE EXPONENTIAL HERE.
    140 
    141         EP    = EXP( MIN((LAMDA(L)*DTAU(L)),TAUMAX))
     157        EP    = EXP( MIN((LAMDA(L)*DTAU(L)),TAUMAX)) ! CLIPPED EXPONENTIAL
    142158        EM    = 1.0D0/EP
    143159        E1(L) = EP+GAMA(L)*EM
    … …  
    146162        E4(L) = GAMA(L)*EP-EM
    147163      END DO
    148  
    149 c     B81=BTOP  ! RENAME BEFORE CALLING DSOLVER - used to be to set
    150 c     B82=BSURF ! them to real*8 - but now everything is real*8
    151 c     R81=RSF   ! so this may not be necessary
    152 
    153 C     Double precision tridiagonal solver
    154 
     164      
     165!      B81=BTOP  ! RENAME BEFORE CALLING DSOLVER - used to be to set
     166!      B82=BSURF ! them to real*8 - but now everything is real*8
     167!      R81=RSF   ! so this may not be necessary
     168
     169! DOUBLE PRECISION TRIDIAGONAL SOLVER
     170     
    155171      CALL DSOLVER(NLAYER,GAMA,CP,CM,CPM1,CMM1,E1,E2,E3,E4,BTOP,
    156172     *             BSURF,RSF,XK1,XK2)
    157  
    158 
    159 C     NOW WE CALCULATE THE FLUXES AT THE MIDPOINTS OF THE LAYERS.
    160  
     173     
     174! NOW WE CALCULATE THE FLUXES AT THE MIDPOINTS OF THE LAYERS.
     175     
    161176      DO L=1,L_NLAYRAD-1
    162         DTAUK = TAUCUM(2*L+1)-TAUCUM(2*L)
    163         EP    = EXP(MIN(LAMDA(L)*DTAUK,TAUMAX)) ! CLIPPED EXPONENTIAL
    164         EM    = 1.0D0/EP
    165         TERM  = UBARI/(1.D0-W0(L)*COSBAR(L))
    166 
    167 C      CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
    168 C      BOTTOM OF THE LAYER.  THAT IS AT DTAU  OPTICAL DEPTH
    169 
    170         CPMID    = B0(L)+B1(L)*DTAUK +B1(L)*TERM
    171         CMMID    = B0(L)+B1(L)*DTAUK -B1(L)*TERM
    172         FMIDP(L) = XK1(L)*EP + GAMA(L)*XK2(L)*EM + CPMID
    173         FMIDM(L) = XK1(L)*EP*GAMA(L) + XK2(L)*EM + CMMID
    174 
    175 C      FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT
    176 
    177         FMIDP(L) = FMIDP(L)*PI
    178         FMIDM(L) = FMIDM(L)*PI
    179       END DO
    180  
    181 C    And now, for the special bottom layer
     177         DTAUK = TAUCUM(2*L+1)-TAUCUM(2*L)
     178         EP    = EXP(MIN(LAMDA(L)*DTAUK,TAUMAX)) ! CLIPPED EXPONENTIAL
     179         EM    = 1.0D0/EP
     180         TERM  = UBARI/(1.D0-W0(L)*COSBAR(L))
     181         
     182! CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
     183! BOTTOM OF THE LAYER.  THAT IS AT DTAU  OPTICAL DEPTH
     184         
     185         CPMID    = B0(L)+B1(L)*DTAUK +B1(L)*TERM
     186         CMMID    = B0(L)+B1(L)*DTAUK -B1(L)*TERM
     187         FMIDP(L) = XK1(L)*EP + GAMA(L)*XK2(L)*EM + CPMID
     188         FMIDM(L) = XK1(L)*EP*GAMA(L) + XK2(L)*EM + CMMID
     189         
     190! FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT
     191         
     192         FMIDP(L) = FMIDP(L)*PI
     193         FMIDM(L) = FMIDM(L)*PI
     194      END DO
     195      
     196! And now, for the special bottom layer
    182197
    183198      L    = L_NLAYRAD
    … …  
    187202      TERM = UBARI/(1.D0-W0(L)*COSBAR(L))
    188203
    189 C    CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
    190 C    BOTTOM OF THE LAYER.  THAT IS AT DTAU  OPTICAL DEPTH
     204! CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
     205! BOTTOM OF THE LAYER.  THAT IS AT DTAU  OPTICAL DEPTH
    191206
    192207      CPMID    = B0(L)+B1(L)*DTAU(L) +B1(L)*TERM
    … …  
    195210      FMIDM(L) = XK1(L)*EP*GAMA(L) + XK2(L)*EM + CMMID
    196211 
    197 C    FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT
    198 
     212! FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT
     213     
    199214      FMIDP(L) = FMIDP(L)*PI
    200215      FMIDM(L) = FMIDM(L)*PI
    201 
    202 C    FLUX AT THE PTOP LEVEL
    203 
     216     
     217! FLUX AT THE PTOP LEVEL
     218     
    204219      EP   = 1.0D0
    205220      EM   = 1.0D0
    206221      TERM = UBARI/(1.0D0-W0(1)*COSBAR(1))
    207 
    208 C    CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
    209 C    BOTTOM OF THE LAYER.  THAT IS AT DTAU  OPTICAL DEPTH
    210 
     222     
     223! CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
     224! BOTTOM OF THE LAYER.  THAT IS AT DTAU  OPTICAL DEPTH
     225     
    211226      CPMID  = B0(1)+B1(1)*TERM
    212227      CMMID  = B0(1)-B1(1)*TERM
    213 
     228     
    214229      FLUXUP = XK1(1)*EP + GAMA(1)*XK2(1)*EM + CPMID
    215230      FLUXDN = XK1(1)*EP*GAMA(1) + XK2(1)*EM + CMMID
    216  
    217 C    FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT
    218 
     231      
     232! FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT
     233     
    219234      FTOPUP = (FLUXUP-FLUXDN)*PI
    220 
    221 
     235     
     236     
    222237      RETURN
    223238      END
Note: See TracChangeset for help on using the changeset viewer.