source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/phys/sfluxi.f @ 224

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