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

Location:
trunk/LMDZ.TITAN/libf/phytitan
Files:
4 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

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

    r1420 r2095  
    4545!!      PARAMETER (NLP=101) ! MUST BE LARGER THAN NLEVEL
    4646
    47       REAL*8 EM, EP
     47      REAL*8 EM, EP, EXPTRM
    4848      REAL*8 W0(L_NLAYRAD), COSBAR(L_NLAYRAD), DTAU(L_NLAYRAD)
    4949      REAL*8 TAU(L_NLEVRAD), WDEL(L_NLAYRAD), CDEL(L_NLAYRAD)
    … …  
    5454      REAL*8 G3(L_NLAYRAD), GAMA(L_NLAYRAD),CP(L_NLAYRAD),CM(L_NLAYRAD)
    5555      REAL*8 CPM1(L_NLAYRAD),CMM1(L_NLAYRAD), E1(L_NLAYRAD)
    56       REAL*8 E2(L_NLAYRAD),E3(L_NLAYRAD),E4(L_NLAYRAD),EXPTRM(L_NLAYRAD)
     56      REAL*8 E2(L_NLAYRAD),E3(L_NLAYRAD),E4(L_NLAYRAD)
    5757      REAL*8 FLUXUP, FLUXDN
    5858      REAL*8 FACTOR, TAUCUMIN(L_LEVELS), TAUCUM(L_LEVELS)
    … …  
    184184
    185185      DO L=1,L_NLAYRAD
    186         EXPTRM(L) = MIN(TAUMAX,LAMDA(L)*DTAU(L))  ! CLIPPED EXPONENTIAL
    187         EP = EXP(EXPTRM(L))
     186        EXPTRM = MIN(TAUMAX,LAMDA(L)*DTAU(L))  ! CLIPPED EXPONENTIAL
     187        EP = EXP(EXPTRM)
    188188
    189189        EM        = 1.0/EP
    … …  
    200200 
    201201      DO L=1,L_NLAYRAD-1
    202         EXPTRM(L) = MIN(TAUMAX,LAMDA(L)*(TAUCUM(2*L+1)-TAUCUM(2*L)))
    203  
    204         EP = EXP(EXPTRM(L))
     202        EXPTRM = MIN(TAUMAX,LAMDA(L)*(TAUCUM(2*L+1)-TAUCUM(2*L)))
     203 
     204        EP = EXP(EXPTRM)
    205205
    206206        EM    = 1.0/EP
    … …  
    240240C     FLUX AT THE Ptop layer
    241241
    242       EP    = 1.0
    243       EM    = 1.0
     242!      EP    = 1.0
     243!      EM    = 1.0
     244C JL18 correction to account for the fact that the radiative top is not at zero optical depth.
     245      EXPTRM = MIN(TAUMAX,LAMDA(L)*(TAUCUM(2)))
     246      EP = EXP(EXPTRM)
     247      EM    = 1.0/EP
    244248      G4    = 1.0-G3(1)
    245249      DENOM = LAMDA(1)**2 - 1./UBAR0**2
    … …  
    260264C     AT THE MIDDLE OF THE LAYER.
    261265
    262       CPMID  = AP
    263       CMMID  = AM
     266C      CPMID  = AP
     267C      CMMID  = AM
     268C JL18 correction to account for the fact that the radiative top is not at zero optical depth.
     269      TAUMID   = TAUCUM(2)
     270      CPMID = AP*EXP(-TAUMID/UBAR0)
     271      CMMID = AM*EXP(-TAUMID/UBAR0)
    264272
    265273      FLUXUP = XK1(1)*EP + GAMA(1)*XK2(1)*EM + CPMID
    … …  
    268276C     ADD THE DIRECT FLUX TO THE DOWNWELLING TERM
    269277
    270       fluxdn = fluxdn+UBAR0*F0PI*EXP(-MIN(TAUCUM(1),TAUMAX)/UBAR0)
     278!      fluxdn = fluxdn+UBAR0*F0PI*EXP(-MIN(TAUCUM(1),TAUMAX)/UBAR0)
     279!JL18 the line above assumed that the top of the radiative model was P=0
     280!   it seems to be better for the IR to use the middle of the last physical layer as the radiative top.
     281!   so we correct the downwelling flux below for the calculation of the heating rate
     282      fluxdn = fluxdn+UBAR0*F0PI*EXP(-TAUCUM(2)/UBAR0)
    271283
    272284C     This is for the "special" bottom layer, where we take
    … …  
    274286
    275287      L     = L_NLAYRAD
    276       EXPTRM(L) = MIN(TAUMAX,LAMDA(L)*(TAUCUM(L_LEVELS)-
     288      EXPTRM = MIN(TAUMAX,LAMDA(L)*(TAUCUM(L_LEVELS)-
    277289     *                                 TAUCUM(L_LEVELS-1)))
    278290
    279       EP    = EXP(EXPTRM(L))
     291      EP    = EXP(EXPTRM)
    280292      EM    = 1.0/EP
    281293      G4    = 1.0-G3(L)

  • trunk/LMDZ.TITAN/libf/phytitan/optcv.F90

    r2050 r2095  
    154154  ! Rayleigh scattering
    155155  do NW=1,L_NSPECTV
    156      do K=2,L_LEVELS
     156     TRAY(1:4,NW)   = 1d-30
     157     do K=5,L_LEVELS
    157158        TRAY(K,NW)   = TAURAY(NW) * DPR(K)
    158159     end do                    ! levels
    … …  
    189190          if (seashaze) dhaze_T(k,nw) = dhaze_T(k,nw)*seashazefact(k)
    190191        ENDIF
     192       
     193        !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
     194        !   but visible does not handle very well diffusion in first layer.
     195        !   The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise)
     196        !   in the 4 first semilayers in optcv, but not optci.
     197        !   This solves random variations of the sw heating at the model top.
     198        if (k<5)  dhaze_T(K,:) = 0.0
    191199         
    192200        DRAYAER = TRAY(K,NW)
    … …  
    311319
    312320  ! Haze scattering
     321            !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
     322            !   but not in the visible
     323            !   The dhaze_s is thus set to 0 in the 4 first semilayers in optcv, but not optci.
     324            !   This solves random variations of the sw heating at the model top.
    313325  DO NW=1,L_NSPECTV
    314     DO K=2,L_LEVELS
     326      DHAZES_T(1:4,NW) = 0.d0
     327    DO K=5,L_LEVELS
    315328      DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW) ! effect of scattering albedo on haze
    316329    ENDDO
    … …  
    365378  DO NG=1,L_NGAUSS       ! full gauss loop
    366379     DO NW=1,L_NSPECTV       
    367         TAUV(1,NW,NG)=0.0D0
    368         DO L=1,L_NLAYRAD
    369            TAUV(L+1,NW,NG)=TAUV(L,NW,NG)+DTAUV(L,NW,NG)
    370         END DO
    371 
    372380        TAUCUMV(1,NW,NG)=0.0D0
    373381        DO K=2,L_LEVELS
    374382           TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
    375383        END DO
     384
     385        DO L=1,L_NLAYRAD
     386           TAUV(L,NW,NG)=TAUCUMV(2*L,NW,NG)
     387        END DO
     388        TAUV(L,NW,NG)=TAUCUMV(2*L_NLAYRAD+1,NW,NG)
    376389     END DO           
    377390  END DO                 ! end full gauss loop

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

    r1781 r2095  
    33     *                  FMNETI,fluxupi,fluxdni,fluxupi_nu,
    44     *                  FZEROI,TAUGSURF)
    5 
     5     
    66      use radinc_h
    77      use radcommon_h, only: planckir, tlimit,sigma, gweight
    88      use comcstfi_mod, only: pi
    9 
     9     
    1010      implicit none
    11 
     11     
    1212      integer NLEVRAD, L, NW, NG, NTS, NTT
    13 
     13     
    1414      real*8 TLEV(L_LEVELS), PLEV(L_LEVELS)
    1515      real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
    … …  
    2424      real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI)
    2525      real*8 FTOPUP
    26 
     26     
    2727      real*8 UBARI, RSFI, TSURF, BSURF, TTOP, BTOP, TAUTOP
    2828      real*8 PLANCK, PLTOP
    … …  
    3030      real*8 FZEROI(L_NSPECTI)
    3131      real*8 taugsurf(L_NSPECTI,L_NGAUSS-1), fzero
    32 
     32     
    3333      real*8 fup_tmp(L_NSPECTI),fdn_tmp(L_NSPECTI)
    3434      real*8 PLANCKSUM,PLANCKREF
    35 
    36 
    37 C======================================================================C
    38  
     35     
     36! AB : variables for interpolation
     37      REAL*8 C1
     38      REAL*8 C2
     39      REAL*8 P1
     40     
     41!======================================================================C
     42     
    3943      NLEVRAD = L_NLEVRAD
    40  
    41 
    42 C     ZERO THE NET FLUXES
    43    
     44     
     45! ZERO THE NET FLUXES
    4446      NFLUXTOPI = 0.0D0
    45 
     47     
    4648      DO NW=1,L_NSPECTI
    4749        NFLUXTOPI_nu(NW) = 0.0D0
    4850        DO L=1,L_NLAYRAD
    4951           FLUXUPI_nu(L,NW) = 0.0D0
    50 
    51               fup_tmp(nw)=0.0D0
    52               fdn_tmp(nw)=0.0D0
    53 
     52           fup_tmp(nw)=0.0D0
     53           fdn_tmp(nw)=0.0D0
    5454        END DO
    5555      END DO
    56 
     56     
    5757      DO L=1,L_NLAYRAD
    5858        FMNETI(L)  = 0.0D0
    … …  
    6060        FLUXDNI(L) = 0.0D0
    6161      END DO
    62  
    63 C    WE NOW ENTER A MAJOR LOOP OVER SPECTRAL INTERVALS IN THE INFRARED
    64 C    TO CALCULATE THE NET FLUX IN EACH SPECTRAL INTERVAL
    65 
     62      
     63! WE NOW ENTER A MAJOR LOOP OVER SPECTRAL INTERVALS IN THE INFRARED
     64! TO CALCULATE THE NET FLUX IN EACH SPECTRAL INTERVAL
     65     
    6666      TTOP  = TLEV(2)  ! JL12 why not (1) ???
    6767      TSURF = TLEV(L_LEVELS)
    … …  
    7171
    7272!JL12 corrects the surface planck function so that its integral is equal to sigma Tsurf^4
    73 !JL12   this ensure that no flux is lost due to:
     73!JL12 this ensure that no flux is lost due to:
    7474!JL12          -truncation of the planck function at high/low wavenumber
    7575!JL12          -numerical error during first spectral integration
    7676!JL12          -discrepancy between Tsurf and NTS/NTfac
    77       PLANCKSUM=0.d0
    78       PLANCKREF=TSURF*TSURF
    79       PLANCKREF=sigma*PLANCKREF*PLANCKREF
     77      PLANCKSUM = 0.d0
     78      PLANCKREF = TSURF * TSURF
     79      PLANCKREF = sigma * PLANCKREF * PLANCKREF
     80     
    8081      DO NW=1,L_NSPECTI
    81          PLANCKSUM=PLANCKSUM+PLANCKIR(NW,NTS)*DWNI(NW)
     82! AB : PLANCKIR(NW,NTS) is replaced by P1, the linear interpolation result for a temperature TSURF
     83         C1 = TSURF * NTfac - int(TSURF * NTfac)
     84         P1 = (1.0D0 - C1) * PLANCKIR(NW,NTS) + C1 * PLANCKIR(NW,NTS+1)
     85         PLANCKSUM = PLANCKSUM + P1 * DWNI(NW)
    8286      ENDDO
    83       PLANCKSUM=PLANCKREF/(PLANCKSUM*Pi)
     87     
     88      PLANCKSUM = PLANCKREF / (PLANCKSUM * Pi)
    8489!JL12
    85 
     90     
    8691      DO 501 NW=1,L_NSPECTI
    87 
    88 C       SURFACE EMISSIONS - INDEPENDENT OF GAUSS POINTS
    89         BSURF = (1.-RSFI)*PLANCKIR(NW,NTS)*PLANCKSUM !JL12 plancksum see above
    90         PLTOP = PLANCKIR(NW,NTT)
    91 
    92 C  If FZEROI(NW) = 1, then the k-coefficients are zero - skip to the
    93 C  special Gauss point at the end.
    94  
    95         FZERO = FZEROI(NW)
    96         IF(FZERO.ge.0.99) goto 40
    97  
    98         DO NG=1,L_NGAUSS-1
     92! SURFACE EMISSIONS - INDEPENDENT OF GAUSS POINTS
     93! AB : PLANCKIR(NW,NTS) is replaced by P1, the linear interpolation result for a temperature TSURF
     94! AB : idem for PLANCKIR(NW,NTT) and PLTOP
     95         C1 = TSURF * NTfac - int(TSURF * NTfac)
     96         C2 = TTOP  * NTfac - int(TTOP  * NTfac)
     97         P1 = (1.0D0 - C1) * PLANCKIR(NW,NTS) + C1 * PLANCKIR(NW,NTS+1)
     98         BSURF = (1. - RSFI) * P1 * PLANCKSUM
     99         PLTOP = (1.0D0 - C2) * PLANCKIR(NW,NTT) + C2*PLANCKIR(NW,NTT+1)
    99100         
    100           if(TAUGSURF(NW,NG).lt. TLIMIT) then
    101             fzero = fzero + (1.0D0-FZEROI(NW))*GWEIGHT(NG)
    102             goto 30
    103           end if
    104 
    105 C         SET UP THE UPPER AND LOWER BOUNDARY CONDITIONS ON THE IR
    106 C         CALCULATE THE DOWNWELLING RADIATION AT THE TOP OF THE MODEL
    107 C         OR THE TOP LAYER WILL COOL TO SPACE UNPHYSICALLY
    108  
    109 !          TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
    110           TAUTOP = TAUCUMI(2,NW,NG)
    111           BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP
    112  
    113 C         WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
    114 C         CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
    115 C         WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
    116          
    117           CALL GFLUXI(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG),
     101! If FZEROI(NW) = 1, then the k-coefficients are zero - skip to the
     102! special Gauss point at the end.
     103         FZERO = FZEROI(NW)
     104         
     105         IF(FZERO.ge.0.99) goto 40
     106         
     107         DO NG=1,L_NGAUSS-1
     108           
     109            if(TAUGSURF(NW,NG).lt. TLIMIT) then
     110               fzero = fzero + (1.0D0-FZEROI(NW))*GWEIGHT(NG)
     111               goto 30
     112            end if
     113           
     114! SET UP THE UPPER AND LOWER BOUNDARY CONDITIONS ON THE IR
     115! CALCULATE THE DOWNWELLING RADIATION AT THE TOP OF THE MODEL
     116! OR THE TOP LAYER WILL COOL TO SPACE UNPHYSICALLY
     117           
     118!            TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
     119            TAUTOP = TAUCUMI(2,NW,NG)
     120            BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP
     121           
     122! WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
     123! CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
     124! WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
     125           
     126            CALL GFLUXI(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG),
    118127     *                TAUCUMI(1,NW,NG),
    119128     *                WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP,
    120129     *                BSURF,FTOPUP,FMUPI,FMDI)
    121 
    122 
    123 
    124 C         NOW CALCULATE THE CUMULATIVE IR NET FLUX
    125 
    126           NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*GWEIGHT(NG)*
    127      *                           (1.0D0-FZEROI(NW))
    128 
    129 c         and same thing by spectral band... (RDW)
    130           NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW)
    131      *      +FTOPUP*DWNI(NW)*GWEIGHT(NG)*(1.0D0-FZEROI(NW))
    132 
    133 
    134           DO L=1,L_NLEVRAD-1
    135 
    136 C           CORRECT FOR THE WAVENUMBER INTERVALS
    137 
    138             FMNETI(L)  = FMNETI(L)+(FMUPI(L)-FMDI(L))*DWNI(NW)*
    139      *                              GWEIGHT(NG)*(1.0D0-FZEROI(NW))
    140             FLUXUPI(L) = FLUXUPI(L) + FMUPI(L)*DWNI(NW)*GWEIGHT(NG)*
    141      *                                (1.0D0-FZEROI(NW))
    142             FLUXDNI(L) = FLUXDNI(L) + FMDI(L)*DWNI(NW)*GWEIGHT(NG)*
    143      *                                (1.0D0-FZEROI(NW))
    144 
    145 c         and same thing by spectral band... (RW)
    146             FLUXUPI_nu(L,NW) = FLUXUPI_nu(L,NW) +
    147      *                FMUPI(L)*DWNI(NW)*GWEIGHT(NG)*(1.0D0-FZEROI(NW))
    148 
    149           END DO
    150 
    151    30     CONTINUE
    152 
    153        END DO       !End NGAUSS LOOP
    154 
    155    40  CONTINUE
    156 
    157 C      SPECIAL 17th Gauss point
    158 
    159        NG     = L_NGAUSS
    160 
    161 !       TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
    162        TAUTOP = TAUCUMI(2,NW,NG)
    163        BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP
    164 
    165 C      WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
    166 C      CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
    167 C      WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
    168 
    169 
    170        CALL GFLUXI(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG),
     130         
     131! NOW CALCULATE THE CUMULATIVE IR NET FLUX
     132            NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*GWEIGHT(NG)
     133     *                * (1.0D0-FZEROI(NW))
     134           
     135! and same thing by spectral band... (RDW)
     136            NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW) + FTOPUP * DWNI(NW)
     137     *                       * GWEIGHT(NG) * (1.0D0-FZEROI(NW))
     138           
     139            DO L=1,L_NLEVRAD-1
     140!           CORRECT FOR THE WAVENUMBER INTERVALS
     141               FMNETI(L)  = FMNETI(L) + (FMUPI(L)-FMDI(L)) * DWNI(NW)
     142     *                    * GWEIGHT(NG)*(1.0D0-FZEROI(NW))
     143               FLUXUPI(L) = FLUXUPI(L) + FMUPI(L)*DWNI(NW)*GWEIGHT(NG)
     144     *                    * (1.0D0-FZEROI(NW))
     145               FLUXDNI(L) = FLUXDNI(L) + FMDI(L)*DWNI(NW)*GWEIGHT(NG)
     146     *                    * (1.0D0-FZEROI(NW))
     147!         and same thing by spectral band... (RW)
     148               FLUXUPI_nu(L,NW) = FLUXUPI_nu(L,NW) + FMUPI(L)*DWNI(NW)
     149     *                          * GWEIGHT(NG) * (1.0D0 - FZEROI(NW))
     150            END DO
     151           
     152   30       CONTINUE
     153         
     154         END DO       !End NGAUSS LOOP
     155         
     156   40    CONTINUE
     157         
     158! SPECIAL 17th Gauss point
     159         NG     = L_NGAUSS
     160         
     161!         TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
     162         TAUTOP = TAUCUMI(2,NW,NG)
     163         BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP
     164         
     165! WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
     166! CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
     167! WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
     168         
     169         CALL GFLUXI(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG),
    171170     *                TAUCUMI(1,NW,NG),
    172171     *                WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP,
    173172     *                BSURF,FTOPUP,FMUPI,FMDI)
    174  
    175 C      NOW CALCULATE THE CUMULATIVE IR NET FLUX
    176 
    177        NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*FZERO
    178 
    179 c         and same thing by spectral band... (RW)
    180           NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW)
     173         
     174! NOW CALCULATE THE CUMULATIVE IR NET FLUX
     175         NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*FZERO
     176         
     177!         and same thing by spectral band... (RW)
     178         NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW)
    181179     *      +FTOPUP*DWNI(NW)*FZERO
    182 
    183        DO L=1,L_NLEVRAD-1
    184 
    185 C        CORRECT FOR THE WAVENUMBER INTERVALS
    186 
    187          FMNETI(L)  = FMNETI(L)+(FMUPI(L)-FMDI(L))*DWNI(NW)*FZERO
    188          FLUXUPI(L) = FLUXUPI(L) + FMUPI(L)*DWNI(NW)*FZERO
    189          FLUXDNI(L) = FLUXDNI(L) + FMDI(L)*DWNI(NW)*FZERO
    190 
    191 c         and same thing by spectral band... (RW)
    192          FLUXUPI_nu(L,NW) = FLUXUPI_nu(L,NW) + FMUPI(L)*DWNI(NW)*FZERO
    193 
    194        END DO
    195 
     180         
     181         DO L=1,L_NLEVRAD-1
     182! CORRECT FOR THE WAVENUMBER INTERVALS
     183            FMNETI(L)  = FMNETI(L)+(FMUPI(L)-FMDI(L))*DWNI(NW)*FZERO
     184            FLUXUPI(L) = FLUXUPI(L) + FMUPI(L)*DWNI(NW)*FZERO
     185            FLUXDNI(L) = FLUXDNI(L) + FMDI(L)*DWNI(NW)*FZERO
     186! and same thing by spectral band... (RW)
     187            FLUXUPI_nu(L,NW) = FLUXUPI_nu(L,NW)
     188     *                       + FMUPI(L) * DWNI(NW) * FZERO
     189         END DO
     190         
    196191  501 CONTINUE      !End Spectral Interval LOOP
    197 
    198 C *** END OF MAJOR SPECTRAL INTERVAL LOOP IN THE INFRARED****
    199 
     192! *** END OF MAJOR SPECTRAL INTERVAL LOOP IN THE INFRARED****
     193     
    200194      RETURN
    201195      END
Note: See TracChangeset for help on using the changeset viewer.