Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

PHY_Atm_S0_RUN.f90 @ 3152

Last change on this file since 3152 was 2089, checked in by Laurent Fairhead, 10 years ago

Inclusion de la physique de MAR

Integration of MAR physics

File size: 17.4 KB

Rev	Line
[2089]	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