Changeset 2056 for trunk


Ignore:
Timestamp:
Jan 7, 2019, 12:26:09 PM (6 years ago)
Author:
aboissinot
Message:

Planck step function is replaced by a piecewise linear function in gfluxi.F and sfluxi.F

Location:
trunk/LMDZ.GENERIC
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/README

    r2055 r2056  
    14161416== 21/12/2018 == MT
    14171417- Officially add a CECILL licence to the LMD Generic GCM. The licence is the same than in the LMDz Earth GCM.
     1418
     1419== 07/01/2018 == AB
     1420- Planck step function is replaced by a piecewise linear function in gfluxi.F and sfluxi.F
     1421in the computation of B1, B0, PLTOP, PLANCKSUM and BSURF.

  • trunk/LMDZ.GENERIC/libf/phystd/gfluxi.F

    r1988 r2056  
    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 
    120 C       CPM1 AND CMM1 ARE THE CPLUS AND CMINUS TERMS EVALUATED
    121 C       AT THE TOP OF THE LAYER, THAT IS ZERO OPTICAL DEPTH
    122 
    123         CPM1(L) = B0(L)+B1(L)*TERM
    124         CMM1(L) = B0(L)-B1(L)*TERM
    125 
    126 C       CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
    127 C       BOTTOM OF THE LAYER.  THAT IS AT DTAU OPTICAL DEPTH.
    128 C       JL18 put CP and CM after the calculation of CPM1 and CMM1 to avoid unecessary calculations.
    129 
    130         CP(L) = CPM1(L) +B1(L)*DTAU(L)
    131         CM(L) = CMM1(L) +B1(L)*DTAU(L)
    132       END DO
    133  
    134 C     NOW CALCULATE THE EXPONENTIAL TERMS NEEDED
    135 C     FOR THE TRIDIAGONAL ROTATED LAYERED METHOD
    136 C     WARNING IF DTAU(J) IS GREATER THAN ABOUT 35 (VAX)
    137 C     WE CLIP IT TO AVOID OVERFLOW.
    138 
     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     
    139156      DO L=1,L_NLAYRAD
    140 
    141 C       CLIP THE EXPONENTIAL HERE.
    142 
    143         EP    = EXP( MIN((LAMDA(L)*DTAU(L)),TAUMAX))
     157        EP    = EXP( MIN((LAMDA(L)*DTAU(L)),TAUMAX)) ! CLIPPED EXPONENTIAL
    144158        EM    = 1.0D0/EP
    145159        E1(L) = EP+GAMA(L)*EM
    … …  
    148162        E4(L) = GAMA(L)*EP-EM
    149163      END DO
    150  
    151 c     B81=BTOP  ! RENAME BEFORE CALLING DSOLVER - used to be to set
    152 c     B82=BSURF ! them to real*8 - but now everything is real*8
    153 c     R81=RSF   ! so this may not be necessary
    154 
    155 C     Double precision tridiagonal solver
    156 
     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     
    157171      CALL DSOLVER(NLAYER,GAMA,CP,CM,CPM1,CMM1,E1,E2,E3,E4,BTOP,
    158172     *             BSURF,RSF,XK1,XK2)
    159  
    160 
    161 C     NOW WE CALCULATE THE FLUXES AT THE MIDPOINTS OF THE LAYERS.
    162  
     173     
     174! NOW WE CALCULATE THE FLUXES AT THE MIDPOINTS OF THE LAYERS.
     175     
    163176      DO L=1,L_NLAYRAD-1
    164         DTAUK = TAUCUM(2*L+1)-TAUCUM(2*L)
    165         EP    = EXP(MIN(LAMDA(L)*DTAUK,TAUMAX)) ! CLIPPED EXPONENTIAL
    166         EM    = 1.0D0/EP
    167         TERM  = UBARI/(1.D0-W0(L)*COSBAR(L))
    168 
    169 C      CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
    170 C      BOTTOM OF THE LAYER.  THAT IS AT DTAU  OPTICAL DEPTH
    171 
    172         CPMID    = B0(L)+B1(L)*DTAUK +B1(L)*TERM
    173         CMMID    = B0(L)+B1(L)*DTAUK -B1(L)*TERM
    174         FMIDP(L) = XK1(L)*EP + GAMA(L)*XK2(L)*EM + CPMID
    175         FMIDM(L) = XK1(L)*EP*GAMA(L) + XK2(L)*EM + CMMID
    176 
    177 C      FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT
    178 
    179         FMIDP(L) = FMIDP(L)*PI
    180         FMIDM(L) = FMIDM(L)*PI
    181       END DO
    182  
    183 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
    184197
    185198      L    = L_NLAYRAD
    … …  
    189202      TERM = UBARI/(1.D0-W0(L)*COSBAR(L))
    190203
    191 C    CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
    192 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
    193206
    194207      CPMID    = B0(L)+B1(L)*DTAU(L) +B1(L)*TERM
    … …  
    197210      FMIDM(L) = XK1(L)*EP*GAMA(L) + XK2(L)*EM + CMMID
    198211 
    199 C    FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT
    200 
     212! FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT
     213     
    201214      FMIDP(L) = FMIDP(L)*PI
    202215      FMIDM(L) = FMIDM(L)*PI
    203 
    204 C    FLUX AT THE PTOP LEVEL
    205 
     216     
     217! FLUX AT THE PTOP LEVEL
     218     
    206219      EP   = 1.0D0
    207220      EM   = 1.0D0
    208221      TERM = UBARI/(1.0D0-W0(1)*COSBAR(1))
    209 
    210 C    CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
    211 C    BOTTOM OF THE LAYER.  THAT IS AT DTAU  OPTICAL DEPTH
    212 
     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     
    213226      CPMID  = B0(1)+B1(1)*TERM
    214227      CMMID  = B0(1)-B1(1)*TERM
    215 
     228     
    216229      FLUXUP = XK1(1)*EP + GAMA(1)*XK2(1)*EM + CPMID
    217230      FLUXDN = XK1(1)*EP*GAMA(1) + XK2(1)*EM + CMMID
    218  
    219 C    FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT
    220 
     231      
     232! FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT
     233     
    221234      FTOPUP = (FLUXUP-FLUXDN)*PI
    222 
    223 
     235     
     236     
    224237      RETURN
    225238      END

  • trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F

    r1781 r2056  
    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.