Context Navigation

forcingdaily_tools.f90

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

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

Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/forcingdaily_tools.f90

Download in other formats: