Context Navigation

sfluxi.F @ 2613

Last change on this file since 2613 was 2283, checked in by jvatant, 5 years ago

Set the temperature boundaries and step for Planck function integration as input in callphys.def, for more flexibility.

+ User can now set them by tplanckmin, tplanckmax and dtplanck
+ Default values are a wide range 30-1500K
+ Add a sanity check in callcorrk instead of leaving out-of-bounds planckir indexes.

--JVO

File size: 7.0 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	! AB : variables for interpolation
37	REAL*8 C1
38	REAL*8 C2
39	REAL*8 P1
40
41	!======================================================================C
42
43	NLEVRAD = L_NLEVRAD
44
45	! ZERO THE NET FLUXES
46	NFLUXTOPI = 0.0D0
47
48	DO NW=1,L_NSPECTI
49	NFLUXTOPI_nu(NW) = 0.0D0
50	DO L=1,L_NLAYRAD
51	FLUXUPI_nu(L,NW) = 0.0D0
52	fup_tmp(nw)=0.0D0
53	fdn_tmp(nw)=0.0D0
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
63	! WE NOW ENTER A MAJOR LOOP OVER SPECTRAL INTERVALS IN THE INFRARED
64	! 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)-NTstart+1
70	NTT = int(TTOP *NTfac)-NTstart+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
81	DO NW=1,L_NSPECTI
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)
86	ENDDO
87
88	PLANCKSUM = PLANCKREF / (PLANCKSUM * Pi)
89	!JL12
90
91	DO 501 NW=1,L_NSPECTI
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)
100
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),
127	* TAUCUMI(1,NW,NG),
128	* WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP,
129	* BSURF,FTOPUP,FMUPI,FMDI)
130
131	! NOW CALCULATE THE CUMULATIVE IR NET FLUX
132	NFLUXTOPI = NFLUXTOPI+FTOPUPDWNI(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),
170	* TAUCUMI(1,NW,NG),
171	* WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP,
172	* BSURF,FTOPUP,FMUPI,FMDI)
173
174	! NOW CALCULATE THE CUMULATIVE IR NET FLUX
175	NFLUXTOPI = NFLUXTOPI+FTOPUPDWNI(NW)FZERO
176
177	! and same thing by spectral band... (RW)
178	NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW)
179	* +FTOPUPDWNI(NW)FZERO
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
191	501 CONTINUE !End Spectral Interval LOOP
192	! * END OF MAJOR SPECTRAL INTERVAL LOOP IN THE INFRARED**
193
194	RETURN
195	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format