Context Navigation

source: tags/ORCHIDEE_4_1/ORCHIDEE/src_driver/forcingdaily_tools.f90 @ 7852

Last change on this file since 7852 was 6135, checked in by jan.polcher, 5 years ago
Adds the option for daily forcings (ISIMIP) to use relative humidity instead of specific air humidity. If both variables are available there is an option to allow to select the preferred method.
File size: 28.1 KB

Line
1	! ==============================================================================================================================\n
2	! MODULE forcingdaily_tools : The general idea of this module is to re-generate a diurnal cycle of forcing variables based
3	! on the daily mean values and Tairmin and Tairmax. The approach is to generate a temporal sub-domain
4	! of size "szsbd" which will correspond to "nbdays" days. Thus "nbdays" daily means will be used
5	! to regenerate the diurnal cycles. Doing more than one day allows to use higher order interpolations
6	! in case it is needed and avoid discontinuities. This process is performed by forcingdaily_gensubd when
7	! ever we come to a new day. The the subroutine forcingdaily_getvalues will extract from the "nbdays" of
8	! reconstructed diurnal cycle the values ORCHIDEE needs.
9	! For most variables we have a specific subroutine to re-generate the diurnal cycle with some specific
10	! parameters which allow to adjust the process.
11	!
12	! forcingdaily_gensubd :
13	! forcingdaily_getvalues :
14	!
15	! CONTACT : jan.polcher@lmd.jussieu.fr
16	!
17	! LICENCE : IPSL (2016)
18	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
19	!
20	!>\BRIEF
21	!!
22	!! RECENT CHANGE(S): None
23	!!
24	!! REFERENCE(S) : None
25	!!
26	!_ ================================================================================================================================
27	!!
28	MODULE forcingdaily_tools
29	!
30	USE defprec
31	USE netcdf
32	!
33	USE ioipsl
34	USE constantes
35	USE solar
36	USE qsat_moisture
37	!
38	USE mod_orchidee_para
39	!
40	IMPLICIT NONE
41	!
42	PRIVATE
43	PUBLIC :: forcingdaily_gensubd, forcingdaily_getvalues
44	PUBLIC :: choice_qair_interpol, qair_interpol
45	!
46	! This PARAMETER essentially manages the memory usage of the module as it
47	! determines how much of the forcing will be uploaded from the netCDF file into
48	! memory.
49	!
50	INTEGER(i_std), SAVE :: current_day = -1
51	REAL(r_std), PARAMETER :: dusk_angle = 0.01
52	INTEGER(i_std), PARAMETER :: nbdays = 3
53	INTEGER(i_std), PARAMETER :: spreadprec = 7200 !! Time over which the precipitation should be distributed (in sec.)
54	REAL(r_std), PARAMETER :: convprec_temp = 20.0 !! Temperature above which all precipitation is supposed to be convective.
55	!! i.e. rainfall occurs only over spreadprec. Below rainfall will last longer.
56	INTEGER(i_std), PARAMETER :: tmaxshift = 10800 !! How long after the solar noon should Tairmax occus ? Time in seconds.
57	INTEGER(i_std), SAVE :: seed = 7865439
58	!
59	INTEGER(i_std), SAVE :: szsubd
60	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: time_subd
61	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: tair_subd, qair_subd
62	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: hurs_subd
63	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: ztq_subd, zuv_subd
64	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: rainf_subd, snowf_subd
65	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: solarang_subd, swdown_subd, lwdown_subd
66	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: u_subd, v_subd, ps_subd
67	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: sinangles
68	INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: idusk, irise, inoon
69	!
70	LOGICAL, SAVE :: choice_qair_interpol = .TRUE.
71	LOGICAL, SAVE :: qair_interpol = .TRUE.
72	CONTAINS
73	!!
74	!! =============================================================================================================================
75	!! SUBROUTINE: forcingdaily_getval
76	!!
77	!!
78	!>\BRIEF Extracts the forcing values needed by ORCHIDEE from teh re-generated diurnal cycles.
79	!!
80	!! DESCRIPTION: As we are in a re-geneted diurnal cycle case, not a lot of precaution is taken to identify the value to be
81	!! used for the next integration interval of ORCHIDEE. Simply the values located closest to the middle of the
82	!! integration interval. This could be improved based on what is coded in forcing_tools.f90.
83	!!
84	!! \n
85	!_ ==============================================================================================================================
86	!
87	SUBROUTINE forcingdaily_getvalues(time_int, dt, ztq, zuv, tair, qair, rainf, snowf, &
88	& swdown, lwdown, solarang, u, v, ps)
89	!
90	! ARGUMENTS
91	!
92	REAL(r_std), INTENT(in) :: time_int(2) !! The time interval over which the forcing is needed.
93	REAL(r_std), INTENT(in) :: dt !! timestep, i.e. distance in seconds between time_int(1) and time_int(2)
94	REAL(r_std), INTENT(out) :: ztq(:), zuv(:)
95	REAL(r_std), INTENT(out) :: tair(:), qair(:), rainf(:), snowf(:)
96	REAL(r_std), INTENT(out) :: swdown(:), lwdown(:), solarang(:)
97	REAL(r_std), INTENT(out) :: u(:), v(:), ps(:)
98	!
99	! LOCAL
100	!
101	INTEGER(i_std) :: imin(1), i, nbpt
102	REAL(r_std) :: tloc
103	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: qair_sat, hurs
104	!
105	tloc = (time_int(1)+time_int(2))/2.0
106	imin = MINLOC(ABS((tloc - (time_subd(:)+current_day))))
107	!
108	nbpt = SIZE(qair)
109	ALLOCATE(qair_sat(nbpt))
110	ALLOCATE(hurs(nbpt))
111	!
112	!
113	tair(:) = tair_subd(:,imin(1))
114	hurs(:) = hurs_subd(:,imin(1))
115	rainf(:) = rainf_subd(:,imin(1))
116	snowf(:) = snowf_subd(:,imin(1))
117	swdown(:) = swdown_subd(:,imin(1))
118	lwdown(:) = lwdown_subd(:,imin(1))
119	solarang(:) = solarang_subd(:,imin(1))
120	ztq(:) = ztq_subd(:,imin(1))
121	zuv(:) = zuv_subd(:,imin(1))
122	u(:) = u_subd(:,imin(1))
123	v(:) = v_subd(:,imin(1))
124	ps(:) = ps_subd(:,imin(1))
125	!
126	! Different options exist in order to get air humidity :
127	! - A linear interpolation of the specific humidity (not ideal !!).
128	! - A linear interpolation of the relative humidity which is then applied to the saturated
129	! value at air temperature (This should be preferred).
130	!
131	IF ( qair_interpol ) THEN
132	qair(:) = qair_subd(:,imin(1))
133	ELSE
134	CALL qsatcalc (nbpt, tair, ps/100.0, qair_sat)
135	IF ( MAXVAL(hurs) > un ) THEN
136	! hurs is certainly in %
137	qair(:) = qair_sat(:)*hurs(:)/100.0
138	ELSE
139	! Here hurs is a ratio.
140	qair(:) = qair_sat(:)*hurs(:)
141	ENDIF
142	ENDIF
143	END SUBROUTINE forcingdaily_getvalues
144	!!
145	!! =============================================================================================================================
146	!! SUBROUTINE: forcingdaily_gensubd
147	!!
148	!!
149	!>\BRIEF generates the sub-diurnal cycle for a number of days around the current time step of the ORCHIDEE simulation.
150	!!
151	!! DESCRIPTION: This routine only works when we start a new day so that not too much work is done. At each new day first the
152	!! the diurnal evolution of the solar angle is computed and then from there all the rest is derived. For swdown
153	!! it is a trivial process but for the other variables more complex procedures are used. For those variables
154	!! nothing could be invented the daily mean value is places et the center of the day and a linear interpolation is
155	!! used.
156	!!
157	!! \n
158	!_ ==============================================================================================================================
159	!
160	SUBROUTINE forcingdaily_gensubd(time_int, dt, iim, jjm, lon, lat, gindex_proc, &
161	& szdom, szslab, time_slab, ztq, zuv, tair, tairmin, tairmax, &
162	& qair, hurs, rainf, snowf, swdown, lwdown, u, v, ps)
163
164	!
165	! ARGUMENTS
166	!
167	REAL(r_std), INTENT(in) :: time_int(2) !! The time interval over which the forcing is needed.
168	REAL(r_std), INTENT(in) :: dt !! timestep, i.e. distance in seconds between time_int(1) and time_int(2)
169	INTEGER(i_std), INTENT(in) :: szdom, szslab
170	INTEGER(i_std), INTENT(in) :: iim, jjm ! Size of 2D domain
171	REAL(r_std), INTENT(in) :: lon(iim,jjm), lat(iim,jjm) ! Longitude and latitude
172	INTEGER(i_std), INTENT(in) :: gindex_proc(szdom)
173	REAL(r_std), INTENT(in) :: time_slab(szslab)
174	REAL(r_std), INTENT(in) :: ztq(szdom,szslab), zuv(szdom,szslab)
175	REAL(r_std), INTENT(in) :: tair(szdom,szslab), tairmin(szdom,szslab), tairmax(szdom,szslab)
176	REAL(r_std), INTENT(in) :: qair(szdom,szslab), hurs(szdom,szslab)
177	REAL(r_std), INTENT(in) :: rainf(szdom,szslab), snowf(szdom,szslab)
178	REAL(r_std), INTENT(in) :: swdown(szdom,szslab), lwdown(szdom,szslab)
179	REAL(r_std), INTENT(in) :: u(szdom,szslab), v(szdom,szslab), ps(szdom,szslab)
180	!
181	! LOCAL
182	!
183	REAL(r_std) :: tloc
184	INTEGER(i_std) :: it, i, ist, imin(1), imax(1), tmin(1), iday
185	INTEGER(i_std) :: stpday, half_subd
186	REAL(r_std) :: julian
187	!
188	! Get the options chosen by the user
189	!
190	!Config Key = DAILY_QAIR_INTERPOL
191	!Config Desc = Decide if qair from the daily forcing should be interpolated
192	!Config Def = false
193	!Config If =
194	!Config Help = The daily forcing can generate sub-diurnal qair either from the
195	!Config daily mean qair through an interpolation, or generate it from the
196	!Config relative humidity and the reconstructed tair. The second option is
197	!Config preferred so the default value for DAILY_QIAR_INTERPOL=false.
198	!Config Units = [-]
199	!
200	!
201	IF ( choice_qair_interpol ) THEN
202	qair_interpol = .TRUE.
203	CALL getin_p('DAILY_QAIR_INTERPOL', qair_interpol)
204	ELSE
205	WRITE(,) "qair_interpol = ", qair_interpol
206	CALL ipslerr(2, 'forcingdaily_gensubd', 'The forcing file has constrained our choices and the user', &
207	& 'choice is not being used.', '')
208	ENDIF
209	!
210	! Set date to middle of requested interval.
211	!
212	tloc = (time_int(1)+time_int(2))/2.0
213	!
214	IF (INT(tloc) .NE. current_day) THEN
215	!
216	! Save the date on which the 3 days are centered
217	!
218	current_day = INT(tloc)
219	!
220	! The sub-diurnal cycle needs to be generated for the 3 days around the current time step.
221	! Allocate memory needed.
222	!
223	IF ( .NOT. ALLOCATED(time_subd) ) THEN
224	!
225	szsubd=INT(nbdays*one_day/dt)
226	CALL random_seed()
227	!
228	ALLOCATE(time_subd(szsubd))
229	ALLOCATE(tair_subd(szdom,szsubd))
230	ALLOCATE(qair_subd(szdom,szsubd))
231	ALLOCATE(hurs_subd(szdom,szsubd))
232	ALLOCATE(rainf_subd(szdom,szsubd))
233	ALLOCATE(snowf_subd(szdom,szsubd))
234	ALLOCATE(swdown_subd(szdom,szsubd))
235	ALLOCATE(lwdown_subd(szdom,szsubd))
236	ALLOCATE(solarang_subd(szdom,szsubd))
237	ALLOCATE(ztq_subd(szdom,szsubd))
238	ALLOCATE(zuv_subd(szdom,szsubd))
239	ALLOCATE(u_subd(szdom,szsubd))
240	ALLOCATE(v_subd(szdom,szsubd))
241	ALLOCATE(ps_subd(szdom,szsubd))
242	ALLOCATE(sinangles(iim,jjm,szsubd))
243	ALLOCATE(idusk(szdom,nbdays))
244	ALLOCATE(irise(szdom,nbdays))
245	ALLOCATE(inoon(szdom,nbdays))
246	ENDIF
247	!
248	! Number of sub-diurnal time steps per day.
249	! half_subd : number of days on either side of current date to be added.
250	stpday = NINT(REAL(szsubd/nbdays))
251	half_subd=INT(nbdays/2.0)
252	!
253	! Generate time axis for the sub-diurnal domain, over the number of days we have (nbdays)
254	!
255	DO it=1,szsubd
256	! Time will be relative to the day at which we are being called.
257	time_subd(it) = ((-half_subdone_day)+(it-0.5)dt)/one_day
258	ENDDO
259	!
260	CALL forcingdaily_solar(tloc, iim, jjm, lon, lat, gindex_proc, &
261	& szdom, szslab, time_slab, swdown, &
262	& szsubd, swdown_subd, sinangles, irise, inoon, idusk)
263	!
264	! Temperature
265	!
266	CALL forcingdaily_tair(tloc, dt, szdom, nbdays, szslab, time_slab, tairmin, tairmax, &
267	& irise, inoon, szsubd, tair_subd)
268	!
269	! LWdown
270	!
271	CALL forcingdaily_lwdown(tloc, dt, szdom, szslab, time_slab, &
272	& lwdown, szsubd, stpday, lwdown_subd)
273	!
274	! Precipitation
275	!
276	CALL forcingdaily_precip(tloc, dt, iim, jjm, lon, lat, gindex_proc, &
277	& szdom, szslab, time_slab, rainf, snowf, tair, &
278	& szsubd, rainf_subd, snowf_subd)
279	!
280	! Qair, U, V, PS
281	!
282	qair_subd(:,:) = undef_sechiba
283	hurs_subd(:,:) = undef_sechiba
284	ztq_subd(:,:) = undef_sechiba
285	zuv_subd(:,:) = undef_sechiba
286	u_subd(:,:) = undef_sechiba
287	v_subd(:,:) = undef_sechiba
288	ps_subd(:,:) = undef_sechiba
289	DO i=1,szdom
290	DO iday=1,nbdays
291	ist = (iday-1)*stpday+stpday/2
292	julian = INT(time_subd(ist) + current_day)+0.5
293	imin = MINLOC(ABS((julian - time_slab(1:szslab))))
294	qair_subd(i,ist) = qair(i,imin(1))
295	hurs_subd(i,ist) = hurs(i,imin(1))
296	ztq_subd(i,ist) = ztq(i,imin(1))
297	zuv_subd(i,ist) = zuv(i,imin(1))
298	u_subd(i,ist) = u(i,imin(1))
299	v_subd(i,ist) = v(i,imin(1))
300	ps_subd(i,ist) = ps(i,imin(1))
301	ENDDO
302	ENDDO
303	CALL forcingdaily_linint(szdom, szsubd, nbdays, qair_subd)
304	CALL forcingdaily_linint(szdom, szsubd, nbdays, hurs_subd)
305	CALL forcingdaily_linint(szdom, szsubd, nbdays, ztq_subd)
306	CALL forcingdaily_linint(szdom, szsubd, nbdays, zuv_subd)
307	CALL forcingdaily_linint(szdom, szsubd, nbdays, u_subd)
308	CALL forcingdaily_linint(szdom, szsubd, nbdays, v_subd)
309	CALL forcingdaily_linint(szdom, szsubd, nbdays, ps_subd)
310	ELSE
311	! Nothing to do as the sub-diurnal cycle has already been generated
312	ENDIF
313	END SUBROUTINE forcingdaily_gensubd
314	!!
315	!! =============================================================================================================================
316	!! SUBROUTINE: forcingdaily_linint
317	!!
318	!>\BRIEF Does a linear interpolation of the variables already placed in the variable.
319	!!
320	!! DESCRIPTION: It will replace the undefined values (x >= undef_sechiba) by the interpolated value. Before entering this
321	!! subroutine all time series need to be set to undef_sechiba with only the point between the interpolation should
322	!! be applied filled with actual values.
323	!!
324	!! \n
325	!_ ==============================================================================================================================
326	!
327	SUBROUTINE forcingdaily_linint(xsz, tsz, nbdays, x)
328	!!
329	!! Arguments
330	!!
331	INTEGER(i_std), INTENT(in) :: xsz, tsz, nbdays
332	REAL(r_std), INTENT(inout) :: x(xsz,tsz)
333	!!
334	!! Local
335	!!
336	INTEGER(i_std) :: ii, it, i
337	INTEGER(i_std) :: find, nind, sind
338	REAL(r_std) :: fval, nval, del
339	!!
340	!
341	DO ii=1,xsz
342	!
343	find = -1
344	fval = undef_sechiba
345	nind = 1
346	nval = undef_sechiba
347	!
348	DO WHILE (nind < tsz)
349	!
350	! Look for the next defined value
351	!
352	DO WHILE (nval >= undef_sechiba .AND. nind < tsz)
353	IF (x(ii,nind) < undef_sechiba) THEN
354	nval=x(ii,nind)
355	ELSE
356	nind=nind+1
357	ENDIF
358	IF (nind >= tsz) THEN
359	nval = undef_sechiba
360	nind = tsz
361	ENDIF
362	ENDDO
363	!
364	! Do the filling or interpolation between find and nind
365	!
366	IF ( find < 0 ) THEN
367	DO i=1,nind
368	x(ii,i) = nval
369	ENDDO
370	ELSE IF (nind == tsz) THEN
371	DO i=find,tsz
372	x(ii,i) = fval
373	ENDDO
374	ELSE
375	del = (nval-fval)/(nind-find)
376	DO i=find,nind
377	x(ii,i) = fval+del*(i-find)
378	ENDDO
379	ENDIF
380	!
381	! Move information to first index
382	!
383	find = nind
384	fval = nval
385	nval = undef_sechiba
386	nind = nind+1
387	ENDDO
388	ENDDO
389	END SUBROUTINE forcingdaily_linint
390	!!
391	!! =============================================================================================================================
392	!! SUBROUTINE: forcingdaily_solar
393	!!
394	!!
395	!>\BRIEF Computes the diurnal cycle of the solarangle and incident solar radiation.
396	!!
397	!! DESCRIPTION: This is very close to what is done in forcing_tools.f90 for the SWdown interpolation. The added code here
398	!! is to compute the time indicis of sun rise, noon and sun set. These will be important point in the
399	!! re-generation of the diurnal cycles.
400	!!
401	!! \n
402	!_ ==============================================================================================================================
403	!
404	SUBROUTINE forcingdaily_solar(tloc, iim, jjm, lon, lat, gindex_proc, &
405	& szdom, szslab, time_slab, swdown_loc, &
406	& szsubd, swdown_subd, sinangles, irise, inoon, idusk)
407	!
408	! Arguments
409	!
410	REAL(r_std), INTENT(in) :: tloc
411	INTEGER(i_std), INTENT(in) :: szdom, szslab
412	INTEGER(i_std), INTENT(in) :: iim, jjm ! Size of 2D domain
413	REAL(r_std), INTENT(in) :: lon(iim,jjm), lat(iim,jjm) ! Longitude and latitude
414	INTEGER(i_std), INTENT(in) :: gindex_proc(szdom)
415	REAL(r_std), INTENT(in) :: time_slab(szslab)
416	REAL(r_std), INTENT(in) :: swdown_loc(szdom,szslab)
417	INTEGER(i_std), INTENT(in) :: szsubd
418	REAL(r_std), INTENT(out) :: swdown_subd(szdom,szsubd)
419	REAL(r_std), INTENT(out) :: sinangles(iim,jjm,szsubd)
420	INTEGER(i_std), INTENT(out) :: idusk(szdom,nbdays)
421	INTEGER(i_std), INTENT(out) :: irise(szdom,nbdays)
422	INTEGER(i_std), INTENT(out) :: inoon(szdom,nbdays)
423	!
424	! Local
425	!
426	REAL(r_std) :: cval, lval
427	INTEGER(i_std) :: it, ii, jj, i, ist, imin(1), imax(1), tmin(1), iday
428	INTEGER(i_std) :: year, month, day, hours, minutes
429	INTEGER(i_std) :: stpday, half_subd
430	REAL(r_std) :: sec, julian
431	REAL(r_std), SAVE :: solaryearstart
432	REAL(r_std) :: sinang(iim,jjm), mean_sinang(iim,jjm,nbdays)
433	INTEGER(i_std) :: nbval(nbdays)
434	LOGICAL, SAVE, ALLOCATABLE, DIMENSION(:) :: mask
435	!
436	! Allocate memory
437	!
438	IF ( .NOT. ALLOCATED(mask) ) THEN
439	ALLOCATE(mask(szsubd))
440	ENDIF
441	!
442	! Get some basic dates and dimensions
443	!
444	stpday = NINT(REAL(szsubd/nbdays))
445	half_subd=INT(nbdays/2.0)
446	CALL ju2ymds (tloc, year, month, day, sec)
447	CALL ymds2ju (year, 1, 1, 0.0, solaryearstart)
448	!
449	mean_sinang(:,:,:) = 0.0
450	nbval(:) = 0
451	!
452	! Compute all solar angles
453	!
454	DO it=1,szsubd
455	julian=time_subd(it) + current_day
456	iday = INT(time_subd(it)+(nbdays-half_subd))
457	CALL solarang (julian, solaryearstart, iim, jjm, lon, lat, sinang)
458	DO ii=1,iim
459	DO jj=1,jjm
460	IF ( sinang(ii,jj) > zero .AND. sinang(ii,jj) < dusk_angle ) THEN
461	sinang(ii,jj) = dusk_angle
462	ENDIF
463	mean_sinang(ii,jj,iday) = mean_sinang(ii,jj,iday)+sinang(ii,jj)
464	ENDDO
465	ENDDO
466	nbval(iday) = nbval(iday)+1
467	!
468	sinangles(:,:,it) = sinang(:,:)
469	!
470	ENDDO
471	DO it=1,nbdays
472	mean_sinang(:,:,it) = mean_sinang(:,:,it)/nbval(it)
473	ENDDO
474	!
475	tmin = MINLOC(ABS((tloc-current_day)-time_subd(1:szsubd)))
476	imin = MINLOC(ABS((tloc - time_slab(1:szslab))))
477	sinang(:,:) = sinangles(:,:,tmin(1))
478	!
479	! Set default values for irise, inoon and idusk so that in polar
480	! regions we are not left with these indicis outside of range.
481	! These values only matter when the sun never rises (sinangles < dusk_angle).
482	!
483	DO iday=1,nbdays
484	irise(:,iday) = (iday-1)*stpday + INT(stpday/3)
485	inoon(:,iday) = (iday-1)*stpday + INT(stpday/2)
486	idusk(:,iday) = (iday-1)*stpday + INT(stpday/2) + INT(stpday/3)
487	ENDDO
488	!
489	DO it=1,szsubd
490	DO i=1,szdom
491	!
492	iday = INT(time_subd(it)+(nbdays-half_subd))
493	! Put Julian date to mid-day of current day
494	julian = INT(time_subd(it) + current_day)+0.5
495	imin = MINLOC(ABS((julian - time_slab(1:szslab))))
496	!
497	jj = ((gindex_proc(i)-1)/iim)+1
498	ii = (gindex_proc(i)-(jj-1)*iim)
499	!
500	IF ( mean_sinang(ii,jj,iday) > zero ) THEN
501	swdown_subd(i,it) = swdown_loc(i,imin(1))*sinangles(ii,jj,it)/mean_sinang(ii,jj,iday)
502	ELSE
503	swdown_subd(i,it) = zero
504	ENDIF
505	!
506	!
507	solarang_subd(i,it) = sinangles(ii,jj,it)
508	!
509	lval = sinangles(ii,jj,MAX(it-1,1))
510	cval = sinangles(ii,jj,it)
511	IF ( lval .LE. dusk_angle .AND. cval .GT. dusk_angle ) THEN
512	irise(i,iday) = it
513	ENDIF
514	IF ( lval .GT. dusk_angle .AND. cval .LE. dusk_angle ) THEN
515	idusk(i,iday) = it
516	ENDIF
517	ENDDO
518	ENDDO
519	!
520	! Position the solar noon in each day
521	!
522	DO i=1,szdom
523
524	jj = ((gindex_proc(i)-1)/iim)+1
525	ii = (gindex_proc(i)-(jj-1)*iim)
526
527	DO it=1,nbdays
528	mask=.FALSE.
529	DO ist=1,stpday
530	mask((it-1)*stpday+ist) = .TRUE.
531	ENDDO
532	imax = MAXLOC(swdown_subd(i,:), mask)
533	IF ( sinangles(ii,jj,imax(1)) > dusk_angle ) THEN
534	inoon(i,it) = imax(1)
535	ENDIF
536	ENDDO
537	ENDDO
538	!
539	END SUBROUTINE forcingdaily_solar
540	!! =============================================================================================================================
541	!! SUBROUTINE: forcingdaily_lwdown
542	!!
543	!!
544	!>\BRIEF Re-generates the LWdown variable.
545	!!
546	!! DESCRIPTION: For the moment LW down is maintained to the mean value over the full diurnal cycle. An idea would be to
547	!! to introduce a diurnal cycle related to the lower atmospheric temperature (Tair).
548	!!
549	!! \n
550	!_ ==============================================================================================================================
551	!
552	!! Interpolation of LWdown
553	!!
554	SUBROUTINE forcingdaily_lwdown(tloc, dt, szdom, szslab, time_slab, &
555	& lwdown, szsubd, stpday, lwdown_subd)
556	!
557	! Arguments
558	!
559	REAL(r_std), INTENT(in) :: tloc, dt
560	INTEGER(i_std), INTENT(in) :: szdom, szslab, stpday, szsubd
561	REAL(r_std), INTENT(in) :: time_slab(szslab)
562	REAL(r_std), INTENT(in) :: lwdown(szdom,szslab)
563	REAL(r_std), INTENT(out) :: lwdown_subd(szdom,szsubd)
564	!
565	! Local
566	!
567	INTEGER(i_std) :: i, iday, it, ist, imin(1)
568	REAL(r_std) :: julian
569	!
570	lwdown_subd(:,:) = undef_sechiba
571	DO i=1,szdom
572	DO iday=1,nbdays
573	DO it=1,stpday
574	ist = (iday-1)*stpday+it
575	julian = INT(time_subd(ist) + current_day)+0.5
576	imin = MINLOC(ABS((julian - time_slab(1:szslab))))
577	lwdown_subd(i,ist) = lwdown(i,imin(1))
578	ENDDO
579	ENDDO
580	ENDDO
581	END SUBROUTINE forcingdaily_lwdown
582	!! =============================================================================================================================
583	!! SUBROUTINE: forcingdaily_tair
584	!!
585	!!
586	!>\BRIEF Re-generates the Tair variable.
587	!!
588	!! DESCRIPTION:
589	!!
590	!! \n
591	!_ ==============================================================================================================================
592	!
593	!! Interpolation of Tair
594	!!
595	SUBROUTINE forcingdaily_tair(tloc, dt, szdom, nbdays, szslab, time_slab, tairmin, tairmax, &
596	& irise, inoon, szsubd, tair_subd)
597	!
598	! Arguments
599	!
600	REAL(r_std), INTENT(in) :: tloc, dt
601	INTEGER(i_std), INTENT(in) :: szdom, szslab, nbdays, szsubd
602	REAL(r_std), INTENT(in) :: time_slab(szslab)
603	REAL(r_std), INTENT(in) :: tairmin(szdom,szslab), tairmax(szdom,szslab)
604	INTEGER(i_std), INTENT(in) :: irise(szdom,nbdays), inoon(szdom,nbdays)
605	REAL(r_std), INTENT(out) :: tair_subd(szdom,szsubd)
606	!
607	! Local
608	!
609	INTEGER(i_std) :: i, it, imin(1), in, ir
610	INTEGER(i_std) :: stpday, nbshift
611	REAL(r_std) :: julian
612	!
613	stpday = NINT(REAL(szsubd/nbdays))
614	nbshift = NINT(REAL(tmaxshift/dt))
615	!
616	tair_subd(:,:) = undef_sechiba
617	DO i=1,szdom
618	DO it=1,nbdays
619	!
620	! Tairmax position
621	!
622	in = (it-1)*stpday + MOD(inoon(i,it)+nbshift-1, stpday)+1
623	!! in = inoon(i,it)+nbshift
624	julian = INT(time_subd(in) + current_day)+0.5
625	imin = MINLOC(ABS((julian - time_slab(1:szslab))))
626	tair_subd(i,in) = tairmax(i,imin(1))
627	!
628	! Tairmin position
629	!
630	ir = irise(i,it)
631	julian = INT(time_subd(ir) + current_day)+0.5
632	imin = MINLOC(ABS((julian - time_slab(1:szslab))))
633	tair_subd(i,ir) = tairmin(i,imin(1))
634	ENDDO
635	ENDDO
636	!
637	CALL forcingdaily_linint(szdom, szsubd, nbdays, tair_subd)
638	!
639	END SUBROUTINE forcingdaily_tair
640	!!
641	!! =============================================================================================================================
642	!! SUBROUTINE: forcingdaily_precip
643	!!
644	!>\BRIEF Distributes the rainfall of a day on a period spreadprec in a random place within the diurnal cycle.
645	!!
646	!! DESCRIPTION: Rainfall is distributed randomly within the day over the spreadprec period. A wide room for improvements here.
647	!! We need to think about the geographical variations of spreadprec and the most lieky time of day for rainfall.
648	!!
649	!! \n
650	!_ ==============================================================================================================================
651	!
652	SUBROUTINE forcingdaily_precip(tloc, dt, iim, jjm, lon, lat, gindex_proc, &
653	& szdom, szslab, time_slab, rainf, snowf, tair, &
654	& szsubd, rainf_subd, snowf_subd)
655	!
656	! Arguments
657	!
658	REAL(r_std), INTENT(in) :: tloc, dt
659	INTEGER(i_std), INTENT(in) :: szdom, szslab
660	INTEGER(i_std), INTENT(in) :: iim, jjm ! Size of 2D domain
661	REAL(r_std), INTENT(in) :: lon(iim,jjm), lat(iim,jjm) ! Longitude and latitude
662	INTEGER(i_std), INTENT(in) :: gindex_proc(szdom)
663	REAL(r_std), INTENT(in) :: time_slab(szslab)
664	REAL(r_std), INTENT(in) :: rainf(szdom,szslab), snowf(szdom,szslab)
665	REAL(r_std), INTENT(in) :: tair(szdom,szslab)
666	INTEGER(i_std), INTENT(in) :: szsubd
667	REAL(r_std), INTENT(out) :: rainf_subd(szdom,szsubd), snowf_subd(szdom,szsubd)
668	!
669	! Local
670	!
671	INTEGER(i_std) :: i, ist, iday, it, ip
672	INTEGER(i_std) :: imin(1)
673	INTEGER(i_std) :: stpday, nbstep
674	REAL(r_std) :: julian, rr, rainlength
675	!
676	stpday = NINT(REAL(szsubd/nbdays))
677	!
678	rainf_subd(:,:) = zero
679	snowf_subd(:,:) = zero
680	!
681	DO i=1,szdom
682	DO iday=1,nbdays
683	!
684	ist = (iday-1)*stpday+INT(stpday/2)
685	julian = INT(time_subd(ist) + current_day)+0.5
686	imin = MINLOC(ABS((julian - time_slab(1:szslab))))
687	!
688	! Air temperature decides how long the rainfall will be !
689	! When it is cold the precipitation will last 3 times longer.
690	!
691	rainlength=spreadprec + MIN(MAX(convprec_temp - (tair(i,imin(1)) - tp_00), 0.0), &
692	convprec_temp)/convprec_temp3.0spreadprec
693	nbstep = NINT(REAL(rainlength/dt))
694	!
695	IF (rainf(i,imin(1)) > zero .OR. snowf(i,imin(1)) > zero) THEN
696	CALL random_number(rr)
697	it = INT(rr*stpday)+1
698	DO ist=1,nbstep
699	ip = MOD(it-1+(ist-1),stpday)+1
700	rainf_subd(i,(iday-1)stpday+ip) = rainf(i,imin(1))one_day/rainlength*dt
701	snowf_subd(i,(iday-1)stpday+ip) = snowf(i,imin(1))one_day/rainlength*dt
702	ENDDO
703	ENDIF
704	ENDDO
705	ENDDO
706	!
707	END SUBROUTINE forcingdaily_precip
708	!
709	END MODULE forcingdaily_tools

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

Download in other formats: