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

Last change on this file since 220 was 135, checked in by aslmd, 14 years ago

CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET
File size: 6.9 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
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 BSURFtest ! by RW for test
35
36      real*8 fup_tmp(L_NSPECTI),fdn_tmp(L_NSPECTI)
37
38
39C======================================================================C
40 
41      NLEVRAD = L_NLEVRAD
42 
43
44C     ZERO THE NET FLUXES
45   
46      NFLUXTOPI = 0.0
47      DO NW=1,L_NSPECTI
48        NFLUXTOPI_nu(NW) = 0.0
49        DO L=1,L_NLAYRAD
50           FLUXUPI_nu(L,NW) = 0.0
51
52              fup_tmp(nw)=0.0
53              fdn_tmp(nw)=0.0
54
55        END DO
56      END DO
57
58      DO L=1,L_NLAYRAD
59        FMNETI(L)  = 0.0
60        FLUXUPI(L) = 0.0
61        FLUXDNI(L) = 0.0
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)
68      TSURF = TLEV(L_LEVELS)
69
70      NTS   = TSURF*10.0D0-499
71      NTT   = TTOP *10.0D0-499
72
73      BSURFtest=0.0
74
75      DO 501 NW=1,L_NSPECTI
76
77C       SURFACE EMISSIONS - INDEPENDENT OF GAUSS POINTS
78        BSURF = (1.-RSFI)*PLANCKIR(NW,NTS) ! interpolate for accuracy??
79        PLTOP = PLANCKIR(NW,NTT)
80
81        !BSURFtest=BSURFtest+BSURF*DWNI(NW)
82        !if(NW.eq.L_NSPECTI)then
83        !     print*,'eps*sigma*T^4',5.67e-8*(1.-RSFI)*TSURF**4
84        !     print*,'BSURF in sfluxi=',pi*BSURFtest
85        !endif
86
87
88C  If FZEROI(NW) = 1, then the k-coefficients are zero - skip to the
89C  special Gauss point at the end.
90 
91        FZERO = FZEROI(NW)
92        IF(FZERO.ge.0.99) goto 40
93 
94        DO NG=1,L_NGAUSS-1
95         
96          if(TAUGSURF(NW,NG).lt. TLIMIT) then
97            fzero = fzero + (1.0-FZEROI(NW))*GWEIGHT(NG)
98            goto 30
99          end if
100
101C         SET UP THE UPPER AND LOWER BOUNDARY CONDITIONS ON THE IR
102C         CALCULATE THE DOWNWELLING RADIATION AT THE TOP OF THE MODEL
103C         OR THE TOP LAYER WILL COOL TO SPACE UNPHYSICALLY
104 
105          TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
106          BTOP   = (1.0-EXP(-TAUTOP/UBARI))*PLTOP
107 
108C         WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
109C         CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
110C         WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
111         
112          CALL GFLUXI(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG),
113     *                TAUCUMI(1,NW,NG),
114     *                WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP,
115     *                BSURF,FTOPUP,FMUPI,FMDI)
116
117
118C         NOW CALCULATE THE CUMULATIVE IR NET FLUX
119
120          NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*GWEIGHT(NG)*
121     *                           (1.0-FZEROI(NW))
122
123c         and same thing by spectral band... (RDW)
124          NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW)
125     *      +FTOPUP*DWNI(NW)*GWEIGHT(NG)*(1.0-FZEROI(NW))
126
127
128          DO L=1,L_NLEVRAD-1
129
130C           CORRECT FOR THE WAVENUMBER INTERVALS
131
132            FMNETI(L)  = FMNETI(L)+(FMUPI(L)-FMDI(L))*DWNI(NW)*
133     *                              GWEIGHT(NG)*(1.0-FZEROI(NW))
134            FLUXUPI(L) = FLUXUPI(L) + FMUPI(L)*DWNI(NW)*GWEIGHT(NG)*
135     *                                (1.0-FZEROI(NW))
136            FLUXDNI(L) = FLUXDNI(L) + FMDI(L)*DWNI(NW)*GWEIGHT(NG)*
137     *                                (1.0-FZEROI(NW))
138
139
140c         and same thing by spectral band... (RW)
141            FLUXUPI_nu(L,NW) = FLUXUPI_nu(L,NW) +
142     *                FMUPI(L)*DWNI(NW)*GWEIGHT(NG)*(1.0-FZEROI(NW))
143
144
145          END DO
146
147         !fup_tmp(nw)=fup_tmp(nw)+FMUPI(L_NLEVRAD-1)*DWNI(NW)*GWEIGHT(NG)
148         !fdn_tmp(nw)=fdn_tmp(nw)+FMDI(L_NLEVRAD-1)*DWNI(NW)*GWEIGHT(NG)
149         !fup_tmp(nw)=fup_tmp(nw)+FMUPI(1)*DWNI(NW)*GWEIGHT(NG)
150         !fdn_tmp(nw)=fdn_tmp(nw)+FMDI(1)*DWNI(NW)*GWEIGHT(NG)
151
152   30     CONTINUE
153
154       END DO       !End NGAUSS LOOP
155
156           ! print*,'-----------------------------------'
157            !print*,'FMDI(',nw,')=',FMDI(L_NLEVRAD-1)
158            !print*,'FMUPI(',nw,')=',FMUPI(L_NLEVRAD-1)
159            !print*,'DWNI(',nw,')=',DWNI(nw)
160           ! print*,'-----------------------------------'
161
162   40  CONTINUE
163
164C      SPECIAL 17th Gauss point
165
166    !   print*,'fzero for ng=17',fzero
167
168
169       NG     = L_NGAUSS
170
171       TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
172       BTOP   = (1.0-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       !     print*,'-----------------------------------'
206       !     print*,'nw=',nw
207       !     print*,'ng=',ng
208       !     print*,'FMDI=',FMDI(L_NLEVRAD-1)
209       !     print*,'FMUPI=',FMUPI(L_NLEVRAD-1)
210       !     print*,'-----------------------------------'
211
212  501 CONTINUE      !End Spectral Interval LOOP
213
214C *** END OF MAJOR SPECTRAL INTERVAL LOOP IN THE INFRARED****
215
216            !print*,'-----------------------------------'
217            !print*,'gweight=',gweight
218            !print*,'FLUXDNI=',FLUXDNI(L_NLEVRAD-1)
219            !print*,'FLUXUPI=',FLUXUPI(L_NLEVRAD-1)
220            !print*,'-----------------------------------'
221
222            !do nw=1,L_NSPECTI
223            !   print*,'fup_tmp(',nw,')=',fup_tmp(nw)
224            !   print*,'fdn_tmp(',nw,')=',fdn_tmp(nw)
225            !enddo
226
227
228      RETURN
229      END
Note: See TracBrowser for help on using the repository browser.