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

Last change on this file since 1834 was 1781, checked in by jvatant, 7 years ago

Useless argument Gweight in rad. tr. routines ( present in radcommon.h )
-JVO

File size: 6.1 KB
Line 
1      SUBROUTINE SFLUXI(PLEV,TLEV,DTAUI,TAUCUMI,UBARI,RSFI,WNOI,DWNI,
2     *                  COSBI,WBARI,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, gweight
8      use comcstfi_mod, only: pi
9
10      implicit none
11
12      integer NLEVRAD, L, NW, NG, NTS, NTT
13
14      real*8 TLEV(L_LEVELS), PLEV(L_LEVELS)
15      real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
16      real*8 FMNETI(L_NLAYRAD)
17      real*8 WNOI(L_NSPECTI), DWNI(L_NSPECTI)
18      real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
19      real*8 FMUPI(L_NLEVRAD), FMDI(L_NLEVRAD)
20      real*8 COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
21      real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
22      real*8 NFLUXTOPI
23      real*8 NFLUXTOPI_nu(L_NSPECTI)
24      real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI)
25      real*8 FTOPUP
26
27      real*8 UBARI, RSFI, TSURF, BSURF, TTOP, BTOP, TAUTOP
28      real*8 PLANCK, PLTOP
29      real*8 fluxupi(L_NLAYRAD), fluxdni(L_NLAYRAD)
30      real*8 FZEROI(L_NSPECTI)
31      real*8 taugsurf(L_NSPECTI,L_NGAUSS-1), fzero
32
33      real*8 fup_tmp(L_NSPECTI),fdn_tmp(L_NSPECTI)
34      real*8 PLANCKSUM,PLANCKREF
35
36
37C======================================================================C
38 
39      NLEVRAD = L_NLEVRAD
40 
41
42C     ZERO THE NET FLUXES
43   
44      NFLUXTOPI = 0.0D0
45
46      DO NW=1,L_NSPECTI
47        NFLUXTOPI_nu(NW) = 0.0D0
48        DO L=1,L_NLAYRAD
49           FLUXUPI_nu(L,NW) = 0.0D0
50
51              fup_tmp(nw)=0.0D0
52              fdn_tmp(nw)=0.0D0
53
54        END DO
55      END DO
56
57      DO L=1,L_NLAYRAD
58        FMNETI(L)  = 0.0D0
59        FLUXUPI(L) = 0.0D0
60        FLUXDNI(L) = 0.0D0
61      END DO
62 
63C     WE NOW ENTER A MAJOR LOOP OVER SPECTRAL INTERVALS IN THE INFRARED
64C     TO CALCULATE THE NET FLUX IN EACH SPECTRAL INTERVAL
65
66      TTOP  = TLEV(2)  ! JL12 why not (1) ???
67      TSURF = TLEV(L_LEVELS)
68
69      NTS   = int(TSURF*NTfac)-NTstar+1
70      NTT   = int(TTOP *NTfac)-NTstar+1
71
72!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:
74!JL12          -truncation of the planck function at high/low wavenumber
75!JL12          -numerical error during first spectral integration
76!JL12          -discrepancy between Tsurf and NTS/NTfac
77      PLANCKSUM=0.d0
78      PLANCKREF=TSURF*TSURF
79      PLANCKREF=sigma*PLANCKREF*PLANCKREF
80      DO NW=1,L_NSPECTI
81         PLANCKSUM=PLANCKSUM+PLANCKIR(NW,NTS)*DWNI(NW)
82      ENDDO
83      PLANCKSUM=PLANCKREF/(PLANCKSUM*Pi)
84!JL12
85
86      DO 501 NW=1,L_NSPECTI
87
88C       SURFACE EMISSIONS - INDEPENDENT OF GAUSS POINTS
89        BSURF = (1.-RSFI)*PLANCKIR(NW,NTS)*PLANCKSUM !JL12 plancksum see above
90        PLTOP = PLANCKIR(NW,NTT)
91
92C  If FZEROI(NW) = 1, then the k-coefficients are zero - skip to the
93C  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
99         
100          if(TAUGSURF(NW,NG).lt. TLIMIT) then
101            fzero = fzero + (1.0D0-FZEROI(NW))*GWEIGHT(NG)
102            goto 30
103          end if
104
105C         SET UP THE UPPER AND LOWER BOUNDARY CONDITIONS ON THE IR
106C         CALCULATE THE DOWNWELLING RADIATION AT THE TOP OF THE MODEL
107C         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 
113C         WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
114C         CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
115C         WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
116         
117          CALL GFLUXI(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG),
118     *                TAUCUMI(1,NW,NG),
119     *                WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP,
120     *                BSURF,FTOPUP,FMUPI,FMDI)
121
122
123
124C         NOW CALCULATE THE CUMULATIVE IR NET FLUX
125
126          NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*GWEIGHT(NG)*
127     *                           (1.0D0-FZEROI(NW))
128
129c         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
136C           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
145c         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
157C      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
165C      WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
166C      CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
167C      WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
168
169
170       CALL GFLUXI(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG),
171     *                TAUCUMI(1,NW,NG),
172     *                WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP,
173     *                BSURF,FTOPUP,FMUPI,FMDI)
174 
175C      NOW CALCULATE THE CUMULATIVE IR NET FLUX
176
177       NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*FZERO
178
179c         and same thing by spectral band... (RW)
180          NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW)
181     *      +FTOPUP*DWNI(NW)*FZERO
182
183       DO L=1,L_NLEVRAD-1
184
185C        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
191c         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
196  501 CONTINUE      !End Spectral Interval LOOP
197
198C *** END OF MAJOR SPECTRAL INTERVAL LOOP IN THE INFRARED****
199
200      RETURN
201      END
Note: See TracBrowser for help on using the repository browser.