Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

PHY_Atm_S0_RUN.f90 @ 5192

Last change on this file since 5192 was 2160, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand
File size: 17.4 KB

Line
1	subroutine PHY_Atm_S0_RUN
2
3	!------------------------------------------------------------------------------+
4	! Mon 17-Jun-2013 MAR \|
5	! MAR PHY_Atm_S0_RUN \|
6	! subroutine PHY_Atm_S0_RUN performs Insolation Computation \|
7	! \|
8	! version 3.p.4.1 created by H. Gallee, Thu 25-Apr-2013 \|
9	! Last Modification by H. Gallee, Mon 17-Jun-2013 \|
10	! \|
11	!------------------------------------------------------------------------------+
12	! \|
13	! REFER.: Ch. Tricot, personal communication \|
14	! ^^^^^^^ M.F. Loutre, personal communication and thesis (1993) \|
15	! \|
16	! INPUT : Mon_TU, Day_TU : Month and Day of the Year \|
17	! ^^^^^^^ HourTU, MinuTU, Sec_TU: Hour, Minute, and Second \|
18	! lon__r(kcolp) : Latitude [radians] \|
19	! lat__h(kcolp) : Longitude [hours] \|
20	! \|
21	! OUTPUT: rsunS0 : Insolation normal to Atmosphere Top (W/m2) \|
22	! ^^^^^^^ csz0S0 : Cosinus of the Zenithal Distance \|
23	! \|
24	!------------------------------------------------------------------------------+
25
26
27	! Global Variables
28	! ================
29
30	use Mod_Real
31	use Mod_PHY____dat
32	use Mod_PHY____grd
33	use Mod_PHY____kkl
34	use Mod_PHY_S0_ctr
35	use Mod_PHY_S0_dat
36	use Mod_PHY_S0_grd
37	use Mod_PHY_S0_kkl
38
39
40
41	! LOCAL VARIABLES
42	! ===============
43
44	use Mod_Atm_S0_RUN
45
46
47	IMPLICIT NONE
48
49
50	integer :: i, j !
51	integer :: ikl !
52	integer :: nj !
53	integer :: lhc !
54
55	real(kind=real8) :: dlamm ! mean long. sun for ma-ja
56	real(kind=real8) :: anm !
57	real(kind=real8) :: ranm !
58	real(kind=real8) :: ranv ! anomalie vraie [radian]
59	real(kind=real8) :: anv ! anomalie vraie [degree]
60	real(kind=real8) :: tls ! longitude vraie [degree]
61	real(kind=real8) :: rlam ! longitude vraie [radian]
62	real(kind=real8) :: sd ! sinus of solar declination angle (delta) [-]
63	real(kind=real8) :: cd ! cosin of solar declination angle (delta) [-]
64	real(kind=real8) :: deltar ! Solar Declination (angle sun at equator) [radian]
65	real(kind=real8) :: delta ! Solar Declination (angle sun at equator) [degree]
66
67	real(kind=real8) :: ddt !
68	real(kind=real8) :: arg !
69	real(kind=real8) :: et !
70	real(kind=real8) :: c1 !
71	real(kind=real8) :: c2 !
72	real(kind=real8) :: c3 !
73	real(kind=real8) :: s1 !
74	real(kind=real8) :: s2 !
75	real(kind=real8) :: s3 !
76
77	real(kind=real8) :: argc !
78	real(kind=real8) :: ahc !
79
80	real(kind=real8) :: LenDay ! Day Length [h]
81	real(kind=real8) :: SunRis ! Time of Sun Rise [h]
82	real(kind=real8) :: SunSet ! Time of Sun Set [h]
83	real(kind=real8) :: timh !
84	real(kind=real8) :: ahor ! Time Angle [hour]
85	real(kind=real8) :: ahorr ! Time Angle [radian]
86	real(kind=real8) :: chor !
87	real(kind=real8) :: zenitr ! RUN
88
89	real(kind=real8) :: omesun ! Sun Azimuth [radian] RUN
90
91	real(kind=real8) :: comes !
92	real(kind=real8) :: somes !
93	real(kind=real8) :: omeswr !
94
95	real(kind=real8) :: r_azim ! azimuth , fine resolution ALL
96
97	integer :: k_azim ! ALL
98	integer :: knazim ! RUN
99
100
101
102
103	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
104	! !
105	! ALLOCATION
106	! ==========
107
108	IF (it_RUN.EQ.1 .OR. FlagDALLOC) THEN !
109	allocate ( cszMIN(kcolp) ) ! MIN cos(zenith.Dist), if Mountain MASK RUN
110	allocate ( sin_sz(kcolp) ) ! sin(Sun Zenith.Dist) RUN
111	allocate ( sin_ho(kcolp) ) ! sin(Sun Hour Angle) RUN
112	ENDIF
113	! !
114	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
115
116
117
118
119	! Insolation at the Top of the Atmosphere (TIME PARAMETERS)
120	! ===============================================================
121
122	! Solar declination : delta
123	! -------------------------
124
125	nj = Day_TU+njYear(Mon_TU) ! Nb of Days since Begin of Year [day]
126
127	dlamm = xlam + (nj-80) * step ! mean long. sun for ma-ja
128	anm = dlamm - xl
129	ranm = anm * Dg2Rad
130	ranv = ranm +(2.0ecc-xe3/4.0)sin(ranm) &! anomalie vraie [radian]
131	& +5.0/4.0eccecc sin(2.0ranm) &!
132	& +13.0/12.0 xe3 sin(3.0*ranm) !
133	anv = ranv / Dg2Rad ! anomalie vraie [degree]
134	tls = anv + xl ! longitude vraie [degree]
135	rlam = tls * Dg2Rad ! longitude vraie [radian]
136
137	sd = so * sin(rlam) ! sinus of solar declination angle (delta) [-]
138	cd = sqrt(un_1 - sd * sd) ! cosin of solar declination angle (delta) [-]
139	! sinus delta = sin(obl)*sin(lambda)
140	! with lambda = real longitude
141	!(Phd.thesis Marie-France Loutre, ASTR-UCL, Belgium, 1993)
142	deltar = atan( sd / cd) !
143	delta = deltar/Dg2Rad ! Solar Declination (degrees, angle sun at equator)
144
145
146	! Eccentricity Effect
147	! -------------------
148
149	dST_UA =(1.0-eccecc)/(1.0+ecccos(ranv))
150	ddt = 1.0/dST_UA ! 1 / normalized earth's sun distance
151
152
153	! Insolation normal to the atmosphere (W/m2)
154	! ------------------------------------------
155
156	rsunS0 = ddt ddt 1368.
157
158
159	! Time Equation (Should maybe be modified in case other than present
160	! ------------- conditions are used, minor impact)
161
162	arg = om*nj
163	c1 = cos( arg)
164	c2 = cos(2.d0*arg)
165	c3 = cos(3.d0*arg)
166	s1 = sin( arg)
167	s2 = sin(2.d0*arg)
168	s3 = sin(3.d0*arg)
169
170	et=0.0072d0c1 -0.0528d0c2 -0.0012d0*c3 &! difference between true solar and mean solar hour angles
171	& -0.1229d0s1 -0.1565d0s2 -0.0041d0*s3 ! (connected to the earth orbital rotation speed)
172
173
174
175
176	! Insolation at the Top of the Troposphere (Auxiliary Variables)
177	! ==============================================================
178
179
180	! Day Length, Time Sunrise and Sunset at Sounding Grid Point (i_x0,j_y0)
181	! ----------------------------------------------------------------------
182
183	IF (LDayS0) THEN
184	i = i_x0
185	j = j_y0
186
187	argc = -tan(lat__r(ikl0))*tan(deltar)
188	if (abs(argc).gt. 1.d0) then
189	ahc = 0.d0
190	if (argc .gt. 1.d0) then
191	lhc = -1
192	LenDay = 0.d0 ! Polar Night
193	else
194	lhc = 1
195	LenDay = 24.d0 ! Midnight Sun
196	end if
197	SunRis = 0.d0
198	SunSet = 0.d0
199	else
200	ahc = acos(argc)
201	lhc = 0
202
203	if(ahc.lt.0.d0) ahc = -ahc
204	ahc = ahc * 12.d0 / piNmbr
205	LenDay = ahc * 2.d0
206	SunRis = 12.d0 - ahc - et
207	SunSet = SunRis + LenDay
208	end if
209
210	END IF
211
212	! Time Angle
213	! ----------
214
215	timh = HourTU + MinuTU / 60.d0 ! time [hour]
216
217	DO ikl=1,kcolp
218	ahor = timh + lon__h(ikl) - 12.d0 - et ! time angle [hour]
219	ahorr = ahor * piNmbr / 12.d0 ! time angle [radian]
220	chor = cos(ahorr)
221
222
223
224
225	! Solar Zenithal Distance zenitr (radians) and
226	! Insolation (W/m2) at the Atmosphere Top ===
227	! =======================================
228
229	csz0S0(ikl) = sinLat(ikl) *sd &
230	& + cosLat(ikl) cd chor
231	! Martin modif to avoid cos(sza)=0 for LMDZ:
232	! csz0SV(ikl) = max(1E-6,csz0S0(ikl))
233	csz0S0(ikl) = max(csz0S0(ikl) ,zer0)
234	csz_S0(ikl) = csz0S0(ikl)
235
236	! Martin control
237	! PRINT*,'csz0S0(',ikl,')=',csz0S0(ikl)
238	! Martin control
239
240	ENDDO
241
242
243
244	! Sun Azimuth
245	! =============
246
247	IF (FaceS0) THEN
248	DO ikl = 1,kcolp
249
250
251	! Slope Impact
252	! ------------
253
254	zenitr = acos(csz0S0(ikl)) ! Solar Zenithal Distance [radians]
255	sin_sz(ikl) = sin(zenitr)
256	sin_ho(ikl) = sin(ahorr)
257	sin_sz(ikl) = max(eps6 ,sin_sz(ikl))
258
259	comes =(sd-sinLat(ikl) *csz0S0(ikl)) &! Cosine of Sun Azimuth
260	& /( cosLat(ikl) *sin_sz(ikl)) !
261	somes =(cd*sin_ho(ikl)) /sin_sz(ikl) ! Sine of Sun Azimuth
262
263	IF (abs(comes).gt.zer0) THEN
264	omesun = atan(somes/comes)
265	IF (comes .lt.zer0) &
266	& omesun = omesun + piNmbr
267
268	IF (omesun.gt. piNmbr) &
269	& omesun = -2.0d0 * piNmbr + omesun
270	IF (omesun.lt. -piNmbr) &
271	& omesun = 2.0d0 * piNmbr + omesun
272
273	ELSE
274	IF (somes .gt.zer0) THEN
275	omesun = 0.5d0 * piNmbr
276	ELSE
277	omesun = 1.5d0 * piNmbr
278	END IF
279	END IF
280
281	IF (ikl.eq.ikl0) &
282	& omeswr = omesun / Dg2Rad
283	omesun = -2.0d0 * piNmbr + omesun + DirAxX * Dg2Rad ! Sun Azimuth (in MAR Reference Frame)
284	! (positive counterclockwise)
285
286	! Minimum Zenithal Distance
287	! ~~~~~~~~~~~~~~~~~~~~~~~~~
288	cszMIN(ikl) = 0.0d00
289	IF (MMskS0) THEN
290	r_azim = omesun / d_azim
291	k_azim = int(r_azim)
292
293	! Mountain S0 MASKING
294	! ~~~~~~~~~~~~~~~~~~~
295	IF(k_azim.le. 0) THEN
296	r_azim = r_azim + n_azim
297	k_azim = k_azim + n_azim
298	END IF
299	knazim = k_azim + 1
300	IF(knazim.gt.n_azim) knazim = knazim - n_azim
301
302	cszMIN(ikl) = cszkS0(ikl,k_azim)+(r_azim - k_azim) &
303	& *(cszkS0(ikl,knazim)-cszkS0(ikl,k_azim))
304	END IF ! (MMskS0)
305
306	! Cosine of Solar Normal Angle
307	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
308	csz_S0(ikl) = slopS0(ikl) *csz0S0(ikl) &
309	& + sin_sz(ikl) * slopS0(ikl)*sqrt(un_1 -slopS0(ikl)) &
310	& * cos(omesun-omenS0(ikl))
311	csz_S0(ikl) = max(zer0 ,csz_S0(ikl))
312	IF (csz0S0(ikl).le.cszMIN(ikl)) &
313	& csz_S0(ikl) = 0.0d00
314	END DO
315	END IF ! (FaceS0)
316
317
318
319
320	! Output Insolation Characteristics
321	! =================================
322
323	IF (MinuTU .eq. 0 .and. Sec_TU .eq. 0) THEN
324
325	ahor = timh + lon__h(ikl0) - 12.d0 - et
326	zenitr= acos(csz0S0(ikl0)) / Dg2Rad
327	! #AZ anormr= acos(csz_S0(ikl0)) / Dg2Rad
328	! #AZ omenwr= DirAxX - omenS0(ikl0) / Dg2Rad
329	! #AZ if (omenwr.lt. 0.) omenwr = omenwr + 360.d0
330	! #AZ if (omenwr.gt.360.) omenwr = omenwr - 360.d0
331	! #AZ omeswr = 360.d0 - omeswr
332	! #AZ if (omeswr.lt. 0.) omeswr = omeswr + 360.d0
333	! #AZ if (omeswr.gt.360.) omeswr = omeswr - 360.d0
334
335	write(4,1) lat__r(ikl0)/Dg2Rad,lon__h(ikl0) &
336	& ,Day_TU,Mon_TU,HourTU,MinuTU,Sec_TU
337	1 format(/,' lat.=',f6.1,3x,'long.=',f7.1,4x,'date :',i3,'-',i2, &
338	& ' / ',i2,' h.UT',i3,' min.',i3,' sec.')
339	write(4,2) i_x0,j_y0,lat__r(ikl0)/Dg2Rad,lon__h(ikl0)
340	2 format(' Sounding at (',i3,i3,') / (',f6.2,'dg,',f6.2,'ho)')
341	write(4,3) rsunS0*csz_S0(ikl0) ,ahor,zenitr &
342	! #AZ& ,omeswr,omenwr , anormr &
343	& ,delta
344	3 format(' Insolation [W/m2] = ',f7.2,' Hor.Angle = ',f7.2, &
345	& ' Zenith.Angle = ',f7.2 &
346	! #AZ& ,/,' ', 7x ,' Sol.Azim. = ',f7.2 &
347	! #AZ& ,/,' ', 7x ,' Nrm.Azim. = ',f7.2 &
348	! #AZ& ,' Normal Angle = ',f7.2 &
349	& ,/,' Solar Declination = ',f7.2)
350
351	IF (LDayS0) THEN
352	if (lhc.eq.-1) &
353	& write(4,4) SunRis,LenDay,SunSet
354	4 format(' Sun Rise Time [h] = ',f7.2,' Day Leng. = ',f7.2, &
355	& ' Sun Set Time = ',f7.2,' -- POLAR NIGHT --')
356	if (lhc.eq. 0) &
357	& write(4,5) SunRis,LenDay,SunSet
358	5 format(' Sun Rise Time [h] = ',f7.2,' Day Leng. = ',f7.2, &
359	& ' Sun Set Time = ',f7.2,' -- SOLAR TIME --')
360	if (lhc.eq. 1) &
361	& write(4,6) SunRis,LenDay,SunSet
362	6 format(' Sun Rise Time [h] = ',f7.2,' Day Leng. = ',f7.2, &
363	& ' Sun Set Time = ',f7.2,' -- MIDNIGHT SUN --')
364	END IF ! (LDayS0)
365
366	END IF !
367
368
369
370
371	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
372	! !
373	! DE-ALLOCATION
374	! =============
375
376	IF (FlagDALLOC) THEN !
377	deallocate ( cszMIN ) ! Sum of cos(zenith.Dist) on Azimut Intervals RUN
378	deallocate ( sin_sz ) ! sin(Sun zenith.Dist) RUN
379	deallocate ( sin_ho ) ! sin(Sun Hour Angle) RUN
380	ENDIF
381	! !
382	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
383
384
385
386
387	return
388	end subroutine PHY_Atm_S0_RUN

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

Download in other formats:

Original Format