source: trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F @ 1157

Last change on this file since 1157 was 961, checked in by jleconte, 12 years ago

15/05/2013 == JL

  • correction in radiative scheme to enforce double precision
  • corrects calculation of ISR
File size: 6.1 KB
Line 
1      SUBROUTINE SFLUXI(PLEV,TLEV,DTAUI,TAUCUMI,UBARI,RSFI,WNOI,DWNI,
2     *                  COSBI,WBARI,GWEIGHT,NFLUXTOPI,NFLUXTOPI_nu,
3     *                  FMNETI,fluxupi,fluxdni,fluxupi_nu,
4     *                  FZEROI,TAUGSURF)
5
6      use radinc_h
7      use radcommon_h, only: planckir, tlimit,sigma
8
9      implicit none
10
11#include "comcstfi.h"
12
13      integer NLEVRAD, L, NW, NG, NTS, NTT
14
15      real*8 TLEV(L_LEVELS), PLEV(L_LEVELS)
16      real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
17      real*8 FMNETI(L_NLAYRAD)
18      real*8 WNOI(L_NSPECTI), DWNI(L_NSPECTI)
19      real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
20      real*8 FMUPI(L_NLEVRAD), FMDI(L_NLEVRAD)
21      real*8 COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
22      real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
23      real*8 GWEIGHT(L_NGAUSS), NFLUXTOPI
24      real*8 NFLUXTOPI_nu(L_NSPECTI)
25      real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI)
26      real*8 FTOPUP
27
28      real*8 UBARI, RSFI, TSURF, BSURF, TTOP, BTOP, TAUTOP
29      real*8 PLANCK, PLTOP
30      real*8 fluxupi(L_NLAYRAD), fluxdni(L_NLAYRAD)
31      real*8 FZEROI(L_NSPECTI)
32      real*8 taugsurf(L_NSPECTI,L_NGAUSS-1), fzero
33
34      real*8 fup_tmp(L_NSPECTI),fdn_tmp(L_NSPECTI)
35      real*8 PLANCKSUM,PLANCKREF
36
37
38C======================================================================C
39 
40      NLEVRAD = L_NLEVRAD
41 
42
43C     ZERO THE NET FLUXES
44   
45      NFLUXTOPI = 0.0D0
46
47      DO NW=1,L_NSPECTI
48        NFLUXTOPI_nu(NW) = 0.0D0
49        DO L=1,L_NLAYRAD
50           FLUXUPI_nu(L,NW) = 0.0D0
51
52              fup_tmp(nw)=0.0D0
53              fdn_tmp(nw)=0.0D0
54
55        END DO
56      END DO
57
58      DO L=1,L_NLAYRAD
59        FMNETI(L)  = 0.0D0
60        FLUXUPI(L) = 0.0D0
61        FLUXDNI(L) = 0.0D0
62      END DO
63 
64C     WE NOW ENTER A MAJOR LOOP OVER SPECTRAL INTERVALS IN THE INFRARED
65C     TO CALCULATE THE NET FLUX IN EACH SPECTRAL INTERVAL
66
67      TTOP  = TLEV(2)  ! JL12 why not (1) ???
68      TSURF = TLEV(L_LEVELS)
69
70      NTS   = int(TSURF*NTfac)-NTstar+1
71      NTT   = int(TTOP *NTfac)-NTstar+1
72
73!JL12 corrects the surface planck function so that its integral is equal to sigma Tsurf^4
74!JL12   this ensure that no flux is lost due to:
75!JL12          -truncation of the planck function at high/low wavenumber
76!JL12          -numerical error during first spectral integration
77!JL12          -discrepancy between Tsurf and NTS/NTfac
78      PLANCKSUM=0.d0
79      PLANCKREF=TSURF*TSURF
80      PLANCKREF=sigma*PLANCKREF*PLANCKREF
81      DO NW=1,L_NSPECTI
82         PLANCKSUM=PLANCKSUM+PLANCKIR(NW,NTS)*DWNI(NW)
83      ENDDO
84      PLANCKSUM=PLANCKREF/(PLANCKSUM*Pi)
85!JL12
86
87      DO 501 NW=1,L_NSPECTI
88
89C       SURFACE EMISSIONS - INDEPENDENT OF GAUSS POINTS
90        BSURF = (1.-RSFI)*PLANCKIR(NW,NTS)*PLANCKSUM !JL12 plancksum see above
91        PLTOP = PLANCKIR(NW,NTT)
92
93C  If FZEROI(NW) = 1, then the k-coefficients are zero - skip to the
94C  special Gauss point at the end.
95 
96        FZERO = FZEROI(NW)
97        IF(FZERO.ge.0.99) goto 40
98 
99        DO NG=1,L_NGAUSS-1
100         
101          if(TAUGSURF(NW,NG).lt. TLIMIT) then
102            fzero = fzero + (1.0D0-FZEROI(NW))*GWEIGHT(NG)
103            goto 30
104          end if
105
106C         SET UP THE UPPER AND LOWER BOUNDARY CONDITIONS ON THE IR
107C         CALCULATE THE DOWNWELLING RADIATION AT THE TOP OF THE MODEL
108C         OR THE TOP LAYER WILL COOL TO SPACE UNPHYSICALLY
109 
110!          TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
111          TAUTOP = TAUCUMI(2,NW,NG)
112          BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP
113 
114C         WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
115C         CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
116C         WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
117         
118          CALL GFLUXI(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG),
119     *                TAUCUMI(1,NW,NG),
120     *                WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP,
121     *                BSURF,FTOPUP,FMUPI,FMDI)
122
123
124
125C         NOW CALCULATE THE CUMULATIVE IR NET FLUX
126
127          NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*GWEIGHT(NG)*
128     *                           (1.0D0-FZEROI(NW))
129
130c         and same thing by spectral band... (RDW)
131          NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW)
132     *      +FTOPUP*DWNI(NW)*GWEIGHT(NG)*(1.0D0-FZEROI(NW))
133
134
135          DO L=1,L_NLEVRAD-1
136
137C           CORRECT FOR THE WAVENUMBER INTERVALS
138
139            FMNETI(L)  = FMNETI(L)+(FMUPI(L)-FMDI(L))*DWNI(NW)*
140     *                              GWEIGHT(NG)*(1.0D0-FZEROI(NW))
141            FLUXUPI(L) = FLUXUPI(L) + FMUPI(L)*DWNI(NW)*GWEIGHT(NG)*
142     *                                (1.0D0-FZEROI(NW))
143            FLUXDNI(L) = FLUXDNI(L) + FMDI(L)*DWNI(NW)*GWEIGHT(NG)*
144     *                                (1.0D0-FZEROI(NW))
145
146c         and same thing by spectral band... (RW)
147            FLUXUPI_nu(L,NW) = FLUXUPI_nu(L,NW) +
148     *                FMUPI(L)*DWNI(NW)*GWEIGHT(NG)*(1.0D0-FZEROI(NW))
149
150          END DO
151
152   30     CONTINUE
153
154       END DO       !End NGAUSS LOOP
155
156   40  CONTINUE
157
158C      SPECIAL 17th Gauss point
159
160       NG     = L_NGAUSS
161
162!       TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
163       TAUTOP = TAUCUMI(2,NW,NG)
164       BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP
165
166C      WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
167C      CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
168C      WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
169
170
171       CALL GFLUXI(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG),
172     *                TAUCUMI(1,NW,NG),
173     *                WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP,
174     *                BSURF,FTOPUP,FMUPI,FMDI)
175 
176C      NOW CALCULATE THE CUMULATIVE IR NET FLUX
177
178       NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*FZERO
179
180c         and same thing by spectral band... (RW)
181          NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW)
182     *      +FTOPUP*DWNI(NW)*FZERO
183
184       DO L=1,L_NLEVRAD-1
185
186C        CORRECT FOR THE WAVENUMBER INTERVALS
187
188         FMNETI(L)  = FMNETI(L)+(FMUPI(L)-FMDI(L))*DWNI(NW)*FZERO
189         FLUXUPI(L) = FLUXUPI(L) + FMUPI(L)*DWNI(NW)*FZERO
190         FLUXDNI(L) = FLUXDNI(L) + FMDI(L)*DWNI(NW)*FZERO
191
192c         and same thing by spectral band... (RW)
193         FLUXUPI_nu(L,NW) = FLUXUPI_nu(L,NW) + FMUPI(L)*DWNI(NW)*FZERO
194
195       END DO
196
197  501 CONTINUE      !End Spectral Interval LOOP
198
199C *** END OF MAJOR SPECTRAL INTERVAL LOOP IN THE INFRARED****
200
201      RETURN
202      END
Note: See TracBrowser for help on using the repository browser.