Context Navigation

radiative_sw.F90 @ 4246

Last change on this file since 4246 was 4245, checked in by dubos, 6 years ago
simple_physics : turn zc, zd, capcal, fluxgrd into local temporaries
File size: 8.1 KB

Line
1	MODULE radiative_sw
2
3	#include "use_logging.h"
4
5	IMPLICIT NONE
6	SAVE
7
8	PRIVATE
9
10	PUBLIC :: sw
11
12	CONTAINS
13
14	PURE SUBROUTINE monGATHER(ngrid, n,index, a,b)
15	INTEGER, INTENT(IN) :: ngrid, n, index(n)
16	REAL, INTENT(IN) :: a(ngrid)
17	REAL, INTENT(OUT) :: b(n)
18	INTEGER :: i
19	IF(n<ngrid) THEN
20	DO i=1,n
21	b(i)=a(index(i))
22	END DO
23	ELSE
24	b(:)=a(:)
25	END IF
26	END SUBROUTINE monGATHER
27
28	PURE subroutine monscatter(ngrid, n,index, b,a)
29	INTEGER, INTENT(IN) :: ngrid, n,index(n)
30	REAL, INTENT(IN) :: b(n)
31	REAL, INTENT(OUT) :: a(ngrid)
32	INTEGER :: i
33	IF(n<ngrid) THEN
34	a(:)=0.
35	DO i=1,n
36	a(index(i))=b(i)
37	END DO
38	ELSE
39	a(:)=b(:)
40	END IF
41	end subroutine monscatter
42
43	SUBROUTINE sw(ngrid,nlayer,ldiurn, coefvis,albedo, &
44	& plevel,ps_rad,pmu,pfract,psolarf0, &
45	& fsrfvis,dtsw, lverbose, lwrite)
46	USE phys_const, ONLY : cpp, g
47	USE writefield_mod, ONLY : writefield
48
49	!=======================================================================
50	!
51	! Rayonnement solaire en atmosphere non diffusante avec un
52	! coefficient d absorption gris.
53	!
54	!=======================================================================
55
56	INTEGER, INTENT(IN) :: ngrid, nlayer
57	LOGICAL, INTENT(IN) :: ldiurn, lverbose, lwrite
58	REAL, INTENT(IN) :: &
59	psolarf0, & ! solar constant
60	ps_rad, coefvis, & ! coefvis = attenuation at p=ps_rad
61	albedo(ngrid), & ! albedo
62	pmu(ngrid), & ! cosine zenithal angle
63	pfract(ngrid), & ! day fraction
64	plevel(ngrid,nlayer+1) ! pressure at interfaces
65	REAL, INTENT(OUT) :: &
66	fsrfvis(ngrid), & ! net surface flux
67	dtsw(ngrid,nlayer) ! temperature tendency
68
69	REAL :: buf1(ngrid), buf2(ngrid, nlayer+1) ! buffers for output
70	! fluxes are non-zero only on those points where the sun shines (mu0>0)
71	! We compute only on those ncount points and gather them to vectorize
72	INTEGER :: ncount, index(ngrid)
73	! In the work arrays below, ngrid should be ncount but ncount is not known yet
74	REAL :: zalb(ngrid), & ! albedo
75	& zmu(ngrid), & ! cosine zenithal angle
76	& zfract(ngrid), & ! day fraction
77	& flux_in(ngrid), & ! incoming solar flux
78	& flux_down(ngrid, nlayer+1), & ! downward flux
79	& flux_up(ngrid, nlayer+1), & ! upward flux
80	& zplev(ngrid,nlayer+1), & ! pressure at interfaces
81	& zflux(ngrid), & ! net surface flux
82	& zdtsw(ngrid,nlayer), & ! temperature tendency
83	& zu(ngrid,nlayer+1)
84	INTEGER :: ig,l,nlevel,igout
85	REAL :: tau0
86
87	nlevel=nlayer+1
88
89	!-----------------------------------------------------------------------
90	! Definitions des tableaux locaux pour les points ensoleilles:
91	! ------------------------------------------------------------
92
93	IF (ldiurn) THEN
94	ncount=0
95	DO ig=1,ngrid
96	index(ig)=0
97	ENDDO
98	DO ig=1,ngrid
99	IF(pfract(ig).GT.1.e-6) THEN
100	ncount=ncount+1
101	index(ncount)=ig
102	ENDIF
103	ENDDO
104	ELSE
105	ncount=ngrid
106	ENDIF
107
108	igout=ncount/2+1
109	CALL monGATHER(ngrid,ncount,index, pfract,zfract)
110	CALL monGATHER(ngrid,ncount,index, pmu, zmu)
111	CALL monGATHER(ngrid,ncount,index, albedo,zalb)
112	DO l=1,nlevel
113	CALL monGATHER(ngrid,ncount,index, plevel(:,l),zplev(:,l))
114	ENDDO
115
116	!-----------------------------------------------------------------------
117	! calcul des profondeurs optiques integres depuis p=0:
118	! ----------------------------------------------------
119
120	! calcul de la partie homogene de l opacite
121	tau0=-.5*log(coefvis)/ps_rad
122	DO l=1,nlayer+1
123	DO ig=1,ncount
124	zu(ig,l)=tau0*zplev(ig,l)
125	ENDDO
126	ENDDO
127
128	!-----------------------------------------------------------------------
129	! 2. calcul de la transmission depuis le sommet de l atmosphere:
130	! -----------------------------------------------------------
131
132	DO ig=1,ncount
133	flux_in(ig) = psolarf0zfract(ig)zmu(ig)
134	ENDDO
135
136	DO l=1,nlevel
137	DO ig=1,ncount
138	flux_down(ig,l) = flux_in(ig)*exp(-zu(ig,l)/zmu(ig))
139	ENDDO
140	ENDDO
141
142	IF (lverbose) THEN
143	WRITELOG(,) 'Diagnostique des transmission dans le spectre solaire'
144	WRITELOG(,) 'zfract, zmu, zalb'
145	WRITELOG(,) zfract(igout), zmu(igout), zalb(igout)
146	WRITELOG(,) 'Pression, quantite d abs, transmission'
147	DO l=1,nlayer+1
148	WRITELOG(,) zplev(igout,l),zu(igout,l),flux_down(igout,l)/flux_in(igout)
149	ENDDO
150	LOG_INFO('rad_sw')
151	ENDIF
152
153	!-----------------------------------------------------------------------
154	! 4. calcul du flux solaire arrivant sur le sol:
155	! ----------------------------------------------
156
157	DO ig=1,ncount
158	zflux(ig) = (1.-zalb(ig))*flux_down(ig,1) ! absorbed (net)
159	flux_up(ig,1) = zalb(ig)*flux_down(ig,1) ! reflected (up)
160	ENDDO
161	IF (lverbose) THEN
162	WRITELOG(,) 'Diagnostique des taux de chauffage solaires:'
163	WRITELOG(,) ' 2 flux solaire net incident sur le sol'
164	WRITELOG(,) zflux(igout)
165	LOG_INFO('rad_sw')
166	ENDIF
167
168	!-----------------------------------------------------------------------
169	! 5.calcul des transmissions depuis le sol, cas diffus:
170	! ------------------------------------------------------
171
172	DO l=1,nlevel
173	DO ig=1,ncount
174	flux_up(ig,l)=flux_up(ig,1)exp(-(zu(ig,1)-zu(ig,l))1.66)
175	ENDDO
176	ENDDO
177
178	IF (lverbose) THEN
179	WRITELOG(,) 'Diagnostique des taux de chauffage solaires'
180	WRITELOG(,) ' 3 transmission avec les sol'
181	WRITELOG(,) 'niveau transmission'
182	DO l=1,nlevel
183	WRITELOG(,) l, flux_up(igout,l)/flux_up(igout,1)
184	ENDDO
185	LOG_INFO('rad_sw')
186	ENDIF
187
188	!-----------------------------------------------------------------------
189	! 10. sorties eventuelles:
190	! ------------------------
191
192	IF(lwrite) THEN
193	CALL monscatter(ngrid,ncount,index, flux_in,buf1)
194	CALL writefield('swtop','SW down TOA','W/m2',buf1)
195
196	DO l=1,nlevel
197	CALL monscatter(ngrid,ncount,index, flux_down(:,l),buf2(:,l))
198	ENDDO
199	CALL writefield('swflux_down','Downward SW flux','W/m2',buf2)
200
201	DO l=1,nlevel
202	CALL monscatter(ngrid,ncount,index, flux_up(:,l),buf2(:,l))
203	ENDDO
204	CALL writefield('swflux_up','Upward SW flux','W/m2',buf2)
205
206	END IF
207
208	!-----------------------------------------------------------------------
209	! 3. taux de chauffage, ray. solaire direct:
210	! ------------------------------------------
211
212	DO l=1,nlayer
213	DO ig=1,ncount
214	! m.cp.dT = dflux/dz
215	! m = -(dp/dz)/g
216	zdtsw(ig,l)=(g/cpp) &
217	& * (flux_down(ig,l+1)-flux_down(ig,l)) &
218	& / (zplev(ig,l)-zplev(ig,l+1))
219	ENDDO
220	ENDDO
221
222	IF (lverbose) THEN
223	WRITELOG(,) 'Diagnostique des taux de chauffage solaires:'
224	WRITELOG(,) ' 1 taux de chauffage lie au ray. solaire direct'
225	DO l=1,nlayer
226	WRITELOG(,) zdtsw(igout,l)
227	ENDDO
228	LOG_INFO('rad_sw')
229	ENDIF
230
231	!-----------------------------------------------------------------------
232	! 6.ajout a l echauffement de la contribution du ray. sol. reflechit:
233	! -------------------------------------------------------------------
234
235	DO l=1,nlayer
236	DO ig=1,ncount
237	zdtsw(ig,l)=zdtsw(ig,l)+ &
238	(g/cpp)*(flux_up(ig,l)-flux_up(ig,l+1)) &
239	/(zplev(ig,l)-zplev(ig,l+1))
240	ENDDO
241	ENDDO
242
243	IF (lverbose) THEN
244	WRITELOG(,) 'Diagnostique des taux de chauffage solaires:'
245	WRITELOG(,) ' 3 taux de chauffage total'
246	DO l=1,nlayer
247	WRITELOG(,) zdtsw(igout,l)
248	ENDDO
249	LOG_INFO('rad_sw')
250	ENDIF
251
252	CALL monscatter(ngrid,ncount,index, zflux,fsrfvis)
253	DO l=1,nlayer
254	CALL monscatter(ngrid,ncount,index, zdtsw(:,l),dtsw(:,l))
255	ENDDO
256
257	END SUBROUTINE sw
258
259	END MODULE radiative_sw

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

Download in other formats:

Original Format