Context Navigation

source: branches/ORCHIDEE_EXT/ORCHIDEE_OL/weather.f90 @ 65

Last change on this file since 65 was 65, checked in by didier.solyga, 14 years ago
Import first version of ORCHIDEE_EXT
File size: 108.5 KB

Line
1	MODULE weather
2	!- IPSL (2006)
3	!- This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
4	!---------------------------------------------------------------------
5	USE netcdf
6	!-
7	USE defprec
8	USE ioipsl
9	USE constantes
10	USE parallel
11	USE grid, ONLY : year,month,day,sec
12	!-
13	IMPLICIT NONE
14	!-
15	PRIVATE
16	PUBLIC weathgen_main, weathgen_domain_size, weathgen_init, weathgen_read_file
17	!
18	! Only for root proc
19	INTEGER, SAVE :: iim_file, jjm_file, llm_file, ttm_file
20	INTEGER,DIMENSION(:,:),SAVE,ALLOCATABLE :: ncorr
21	INTEGER,DIMENSION(:,:,:),SAVE,ALLOCATABLE :: icorr,jcorr
22	INTEGER,SAVE :: i_cut, n_agg
23
24	! climatological wet days + anomaly (days/month)
25	REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xinwet
26	! climatological precipition + anomaly (mm/day)
27	REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xinprec
28	! climatological temp + anomaly (C)
29	REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xint
30	! climatological relative humidity + anomaly (%)
31	REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xinq
32	! climatological wind speed + anomaly (m s-1)
33	REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xinwind
34	! climatological cloudiness + anomaly(%)
35	REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xincld
36	! climatological temp range + anomaly(C)
37	REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xintrng
38	! topography (m)
39	REAL,DIMENSION(:),SAVE,ALLOCATABLE :: xintopo
40	! latitudes of land points
41	REAL,DIMENSION(:),SAVE,ALLOCATABLE :: lat_land
42	!-
43	! daily values
44	!-
45	REAL,SAVE :: julian_last
46	! flag for wet day / dry day
47	INTEGER,DIMENSION(:),SAVE,ALLOCATABLE :: iwet
48	!-
49	! today's values (m0 means "minus 0")
50	!-
51	! surface pressure (Pa)
52	REAL,DIMENSION(:),SAVE,ALLOCATABLE :: psurfm0
53	! cloud fraction
54	REAL,DIMENSION(:),SAVE,ALLOCATABLE :: cloudm0
55	! maximum daily temperature (K)
56	REAL,DIMENSION(:),SAVE,ALLOCATABLE :: tmaxm0
57	! minimum daily temperature (K)
58	REAL,DIMENSION(:),SAVE,ALLOCATABLE :: tminm0
59	! daily average specific humidity (kg_h2o/kg_air)
60	REAL,DIMENSION(:),SAVE,ALLOCATABLE :: qdm0
61	! daily average wind speed (m/sec)
62	REAL,DIMENSION(:),SAVE,ALLOCATABLE :: udm0
63	! daily precitation (mm/day)
64	REAL,DIMENSION(:),SAVE,ALLOCATABLE :: precipm0
65	!-
66	! yesterday's values (m1 means "minus 1")
67	!-
68	! surface pressure (Pa)
69	REAL,DIMENSION(:),SAVE,ALLOCATABLE :: psurfm1
70	! cloud fraction
71	REAL,DIMENSION(:),SAVE,ALLOCATABLE :: cloudm1
72	! maximum daily temperature (K)
73	REAL,DIMENSION(:),SAVE,ALLOCATABLE :: tmaxm1
74	! minimum daily temperature (K)
75	REAL,DIMENSION(:),SAVE,ALLOCATABLE :: tminm1
76	! daily average specific humidity (kg_h2o/kg_air)
77	REAL,DIMENSION(:),SAVE,ALLOCATABLE :: qdm1
78	! daily average wind speed (m/sec)
79	REAL,DIMENSION(:),SAVE,ALLOCATABLE :: udm1
80	! daily precitation (mm/day)
81	REAL,DIMENSION(:),SAVE,ALLOCATABLE :: precipm1
82	!-
83	! other
84	!-
85	! statistical (0) or prescribed (1) daily values
86	INTEGER,SAVE :: ipprec
87	! respect monthly precipitation
88	LOGICAL,SAVE :: precip_exact
89	INTEGER,DIMENSION(31,12),SAVE :: jour_precip
90	! max size of random seed
91	INTEGER,PARAMETER :: seedsize_max = 300
92	LOGICAL,SAVE :: dump_weather
93	CHARACTER(LEN=20),SAVE :: dump_weather_file
94	LOGICAL,SAVE :: gathered
95	INTEGER,SAVE :: dump_id
96	!
97	! Absolute zero
98	REAL,PARAMETER :: zero_t=273.16
99	!-
100	! Parametres orbitaux:
101	!-
102	! Eccentricity
103	REAL,SAVE :: ecc
104	! Longitude of perihelie
105	REAL,SAVE :: perh
106	! obliquity
107	REAL,SAVE :: xob
108	!-
109	INTEGER,PARAMETER :: nmon = 12
110	!
111	CHARACTER(LEN=3),DIMENSION(12) :: cal = &
112	& (/ 'JAN','FEB','MAR','APR','MAY','JUN', &
113	& 'JUL','AUG','SEP','OCT','NOV','DEC' /)
114	INTEGER,DIMENSION(12),SAVE :: ndaypm = &
115	& (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
116	INTEGER,SAVE :: soldownid, rainfid, snowfid, lwradid, &
117	& tairid, qairid, psolid, uid, vid, &
118	& time_id, timestp_id
119	!-
120	! Parameters for NETCDF :
121	!-
122	INTEGER,SAVE :: n_rtp = nf90_real4
123	!-
124	! Flag for dynamic allocations :
125	INTEGER :: ALLOC_ERR
126	!-
127	! Calendar type
128	CHARACTER(LEN=20),SAVE :: calendar_str
129	!-
130	! Land points index
131	INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: kindex_w
132	INTEGER, SAVE :: nbindex_w
133	!-
134	! Plot of projection file => grid, number of column on terminal
135	INTEGER, PARAMETER :: termcol = 100
136	!80
137	CONTAINS
138	!-
139	!===
140	!-
141	SUBROUTINE daily &
142	& (npoi, imonth, iday, cloud, tmax, tmin, precip, qd, ud, psurf)
143	!---------------------------------------------------------------------
144	! overview
145	!
146	! this routine generates daily weather conditions from monthly-mean
147	! climatic parameters
148	!
149	! specifically, this routine generates daily values of
150	!
151	! - daily total precipitation
152	! - daily maximum temperature
153	! - daily minimum temperature
154	! - daily average cloud cover
155	! - daily average relative humidity
156	! - daily average wind speed
157	!
158	! in order to generate daily weather conditions, the model uses
159	! a series of 'weather generator' approaches,
160	! which generate random combinations of weather conditions
161	! based upon the climatological conditions in general,
162	! this weather generator is based upon the so-called Richardson
163	! weather generator
164	!
165	! appropriate references include:
166	!
167	! Geng, S., F.W.T. Penning de Vries, and L. Supit, 1985: A simple
168	! method for generating rainfall data, Agricultural and Forest
169	! Meteorology, 36, 363-376.
170	!
171	! Richardson, C. W. and Wright, D. A., 1984: WGEN: A model for
172	! generating daily weather variables: U. S. Department of
173	! Agriculture, Agricultural Research Service.
174	!
175	! Richardson, C., 1981: Stochastic simulation of daily
176	! precipitation, temperature, and solar radiation. Water Resources
177	! Research 17, 182-190.
178	!---------------------------------------------------------------------
179	!-
180	! in & out: global variables
181	!-
182	! wet day / dry day flag
183	! INTEGER,INTENT(INOUT):: iwet(npoi)
184	!-
185	! input
186	!-
187	! total number of land points
188	INTEGER,INTENT(IN) :: npoi
189	INTEGER,INTENT(IN) :: imonth, iday
190	!-
191	! output
192	!-
193	! surface pressure (Pa)
194	REAL,INTENT(OUT) :: psurf(npoi)
195	! cloud fraction
196	REAL,INTENT(OUT) :: cloud(npoi)
197	! maximum daily temperature (K)
198	REAL,INTENT(OUT) :: tmax(npoi)
199	! maximum daily temperature (K)
200	REAL,INTENT(OUT) :: tmin(npoi)
201	! daily average specific humidity (kg_h2o/kg_air)
202	REAL,INTENT(OUT) :: qd(npoi)
203	! daily average wind speed (m/sec)
204	REAL,INTENT(OUT) :: ud(npoi)
205	! daily precitation (mm/day)
206	REAL,INTENT(OUT) :: precip(npoi)
207	!-
208	! local
209	!-
210	! daily average temperature (K)
211	REAL :: td(npoi)
212	REAL,PARAMETER :: rair = 287.
213	REAL,PARAMETER :: grav = 9.81
214	!!$ REAL,PARAMETER :: pi = 3.1415927
215	!-
216	! weather generator 'memory' matrix
217	!-
218	REAL,allocatable,save,dimension(:,:) :: xstore
219	!-
220	REAL :: ee(3), r(3), rr(npoi,3)
221	REAL :: alpha(npoi), rndnum, pwd, pww
222	REAL :: beta(npoi)
223	REAL :: pwet(npoi)
224	REAL :: rwork
225	REAL :: omcloud, omqd, omtmax
226	REAL :: cloudm, cloudw, cloudd
227	REAL :: cloude(npoi), clouds(npoi)
228	REAL :: tmaxd, tmaxw, tmaxm
229	REAL :: tminm
230	REAL :: tmins(npoi), tmaxs(npoi)
231	REAL :: tmine(npoi), tmaxe(npoi)
232	REAL :: qdm(npoi),qdd(npoi),qde(npoi),qdw(npoi),qdup(npoi),qdlow(npoi)
233	INTEGER :: i,j,k
234	REAL :: amn,b1,b2,b3,eud,rn,rn1,rn2,rn3,rn4,s1,s2,s12,x1,y, z(npoi)
235	REAL :: aa(npoi),ab(npoi),tr1(npoi), tr2(npoi)
236	REAL :: tdm,trngm,tdum
237	REAL :: qsattd(npoi)
238	INTEGER :: it1w, it2w
239	REAL :: dt
240	REAL :: rainpwd(npoi)
241	INTEGER :: not_ok(npoi)
242	INTEGER :: count_not_ok,count_not_ok_g
243	LOGICAL,SAVE :: firstcall = .TRUE.
244	INTEGER,save :: npoi0
245	REAL :: xx,e
246	!-
247	! define autocorrelation matrices for Richardson generator
248	!
249	! note that this matrix should be based upon a statistical
250	! analysis of regional weather patterns
251	!
252	! for global simulations, we use 'nominal' values
253	!-
254	REAL, DIMENSION(3,3) :: a,b
255	! Warnings
256	LOGICAL :: Warning_aa_ab(npoi), Warning_iwet(npoi)
257
258	!---------------------------------------------------------------------
259	!-
260	! initial setup for daily climate calculations
261	!-
262	a(1,:) = (/ 0.600, 0.500, 0.005 /)
263	a(2,:) = (/ 0.010, 0.250, 0.005 /)
264	a(3,:) = (/ 0.020, 0.125, 0.250 /)
265	!-
266	b(1,:) = (/ 0.500, 0.250, -0.250 /)
267	b(2,:) = (/ 0.000, 0.500, 0.250 /)
268	b(3,:) = (/ 0.000, 0.000, 0.500 /)
269	!-
270	e = EXP(1.)
271	! GK240100
272	IF (firstcall) THEN
273	firstcall = .FALSE.
274	ALLOC_ERR=-1
275	ALLOCATE(xstore(npoi,3), STAT=ALLOC_ERR)
276	IF (ALLOC_ERR/=0) THEN
277	WRITE(numout,*) "ERROR IN ALLOCATION of xstore : ",ALLOC_ERR
278	STOP
279	ENDIF
280	xstore(:,:) = 0.
281	npoi0 = npoi
282	ELSE IF (npoi /= npoi0) THEN
283	WRITE(numout,*) 'Domain size old, new: ',npoi0,npoi
284	STOP 'WG Daily: Problem: Domain has changed since last call'
285	ENDIF
286	!-
287	! define working variables
288	!-
289	rwork = (grav/rair/0.0065)
290	!-
291	! 'omega' parameters used to calculate differences in expected
292	! climatic parameters on wet and dry days
293	!
294	! following logic of weather generator used in the EPIC crop model
295	!
296	! omcloud -- cloud cover
297	! omqd -- humidity
298	! omtmax -- maximum temperature
299	!-
300	omcloud = 0.90 ! originally 0.90
301	omqd = 0.50 ! originally 0.50
302	omtmax = 0.75 ! originally 0.75
303	!-
304	! calculate weighting factors used in interpolating climatological
305	! monthly-mean input values to daily-mean values
306	!-
307	! this is a simple linear interpolation technique that takes into
308	! account the length of each month
309	!-
310	IF (ipprec == 0) THEN
311	IF (REAL(iday) < REAL(ndaypm(imonth)+1)/2.0) then
312	it1w = imonth-1
313	it2w = imonth
314	dt = (REAL(iday)-0.5)/ndaypm(imonth)+0.5
315	ELSE
316	it1w = imonth
317	it2w = imonth+1
318	dt = (REAL(iday)-0.5)/ndaypm(imonth)-0.5
319	ENDIF
320	if (it1w < 1) it1w = 12
321	if (it2w > 12) it2w = 1
322	ELSE
323	dt = -1.
324	it1w = -1
325	it2w = -1
326	ENDIF
327	!-
328	IF (ipprec == 0) THEN
329	!---
330	!-- use weather generator to create daily statistics
331	!---
332	! (1) determine if today will rain or not
333	!---
334	!-- calculate monthly-average probability of rainy day
335	!---
336	DO i=1,npoi
337	pwet(i) = xinwet(i,imonth)/ndaypm(imonth)
338	ENDDO
339	!---
340	IF (.NOT.precip_exact) THEN
341	!-----
342	!---- (1.1) following Geng et al.
343	!-----
344	IF (is_root_prc) THEN
345	CALL random_number (rndnum)
346	ENDIF
347	CALL bcast(rndnum)
348	!-----
349	DO i=1,npoi
350	IF (xinprec(i,imonth) > 1.e-6) THEN
351	!---------
352	!-------- implement simple first-order Markov-chain precipitation
353	!-------- generator logic based on Geng et al. (1986),
354	!-------- Richardson and Wright (1984), and Richardson (1981)
355	!---------
356	!-------- basically, this allows for the probability of today being
357	!-------- a wet day (a day with measureable precipitation)
358	!-------- to be a function of what yesterday was (wet or dry)
359	!---------
360	!-------- the logic here is that it is more likely that a wet day
361	!-------- will follow another wet day -- allowing
362	!-------- for 'storm events' to persist
363	!---------
364	!-------- estimate the probability of a wet day after a dry day
365	!---------
366	pwd = 0.75*pwet(i)
367	!---------
368	!-------- estimate the probability of a wet day after a wet day
369	!---------
370	pww = 0.25+pwd
371	!---------
372	!-------- decide if today is a wet day or a dry day
373	!-------- using a random number
374	!---------
375	!-------- call random_number(rndnum) ! done before
376	!---------
377	IF (iwet(i) == 0) then
378	IF (rndnum <= pwd) iwet(i) = 1
379	ELSE
380	IF (rndnum > pww) iwet(i) = 0
381	ENDIF
382	ELSE
383	iwet(i) = 0
384	ENDIF
385	ENDDO
386	ELSE
387	!-----
388	!---- (1.2) preserving the monthly number of precip days
389	!---- and monthly precip
390	!-----
391	DO i=1,npoi
392	IF (ABS(xinwet(i,imonth)) < 32.) THEN
393	IF (xinprec(i,imonth) > 1.e-6) THEN
394	IF ( jour_precip(iday,imonth) &
395	& <= NINT(MAX(1.,xinwet(i,imonth))) ) THEN
396	iwet(i) = 1
397	ELSE
398	iwet(i) = 0
399	ENDIF
400	ELSE
401	iwet(i) = 0
402	ENDIF
403	ENDIF
404	ENDDO
405	ENDIF
406	!---
407	!-- (2) determine today's precipitation amount
408	!---
409	IF (.not.precip_exact) THEN
410	Warning_aa_ab(:)=.FALSE.
411	Warning_iwet(:)=.FALSE.
412	!-----
413	!---- (2.1) following Geng et al.
414	!-----
415	aa(:) = zero
416	ab(:) = zero
417	tr2(:)= zero
418	tr1(:)= zero
419	beta(:) = un
420	DO i=1,npoi
421	!-------
422	!------ initialize daily precipitation to zero
423	!-------
424	precip(i) = 0.0
425	!-------
426	IF (xinprec(i,imonth) > 1.e-6) THEN
427	!---------
428	!-------- if it is going to rain today
429	!---------
430	IF (iwet(i) == 1) THEN
431	!-----------
432	!---------- calculate average rainfall per wet day
433	!-----------
434	rainpwd(i) = xinprec(i,imonth) &
435	& *ndaypm(imonth)/MAX(0.1,xinwet(i,imonth))
436	!-----------
437	!---------- randomly select a daily rainfall amount
438	!---------- from a probability density function of rainfall
439	!----------
440	!---------- method i --
441	!-----------
442	!---------- use the following technique from Geng et al. and Richardson
443	!---------- to distribute rainfall probabilities
444	!-----------
445	!---------- pick a random rainfall amount
446	!---------- from a two-parameter gamma function distribution function
447	!-----------
448	!---------- estimate two parameters for gamma function
449	!---------- (following Geng et al.)
450	!-----------
451	beta(i) = MAX(1.0,-2.16+1.83*rainpwd(i))
452	alpha(i) = rainpwd(i)/beta(i)
453	!-----------
454	!---------- determine daily precipitation amount
455	!---------- from gamma distribution function
456	!---------- (following WGEN code of Richardson and Wright (1984))
457	!-----------
458	IF (ABS(1.-alpha(i)) < 1.e-6) THEN
459	alpha(i) = 1.e-6*(alpha(i)/ABS(alpha(i)))
460	ENDIF
461	aa(i) = 1.0/alpha(i)
462	ab(i) = 1.0/(1.0-alpha(i))
463	!-----------
464	IF ( (ABS(aa(i)) < 1.e-6) .OR. (ABS(ab(i)) < 1.e-6) ) THEN
465	Warning_aa_ab(:)=.TRUE.
466	ENDIF
467	tr1(i) = exp(-18.42/aa(i))
468	tr2(i) = exp(-18.42/ab(i))
469	ENDIF
470	ELSE
471	IF (iwet(i) == 1) THEN
472	Warning_iwet(i)=.TRUE.
473	ENDIF
474	ENDIF
475	ENDDO
476
477	DO i=1,npoi
478	IF ( Warning_aa_ab(i) ) THEN
479	WRITE(numout,*) ' ATTENTION, aa ou ab:'
480	WRITE(numout,*) ' aa, ab = ',aa(i),ab(i)
481	WRITE(numout,*) ' alpha, rainpwd, beta =', &
482	& alpha(i),rainpwd(i),beta(i)
483	ENDIF
484	IF ( Warning_iwet(i) ) THEN
485	WRITE(numout,*) ' ATTENTION, iwet = 1 alors que xinprec = 0)'
486	WRITE(numout,*) ' xinprec, iwet = ',xinprec(i,imonth),iwet(i)
487	ENDIF
488	ENDDO
489	!-----
490	where ( iwet(:) == 1 )
491	not_ok(:) = 1
492	elsewhere
493	not_ok(:) = 0
494	endwhere
495	!-----
496	count_not_ok_g=0
497	count_not_ok=SUM(not_ok(:))
498	CALL reduce_sum(count_not_ok,count_not_ok_g)
499	CALL bcast(count_not_ok_g)
500	!-
501	z(:)=0.0
502	DO WHILE (count_not_ok_g > 0)
503	IF (is_root_prc) THEN
504	CALL random_number (rn1)
505	CALL random_number (rn2)
506	ENDIF
507	CALL bcast(rn1)
508	CALL bcast(rn2)
509
510	DO i=1,npoi
511	IF ((iwet(i) == 1).AND.(not_ok(i) == 1)) then
512	IF ( (rn1-tr1(i)) <= 0. ) THEN
513	s1 = 0.0
514	ELSE
515	s1 = rn1**aa(i)
516	ENDIF
517	!-----------
518	IF ((rn2-tr2(i)) <= 0.) THEN
519	s2 = 0.0
520	ELSE
521	s2 = rn2**ab(i)
522	ENDIF
523	!-----------
524	s12 = s1+s2
525	IF ((s12-1.0) <= 0.) THEN
526	z(i) = s1/s12
527	not_ok(i) = 0
528	ENDIF
529	ENDIF
530	ENDDO
531
532	count_not_ok_g=0
533	count_not_ok=SUM(not_ok(:))
534	CALL reduce_sum(count_not_ok,count_not_ok_g)
535	CALL bcast(count_not_ok_g)
536	ENDDO
537	!-----
538	IF (is_root_prc) THEN
539	CALL random_number (rn3)
540	ENDIF
541	CALL bcast(rn3)
542	! WRITE(,) mpi_rank,"rn3=",rn3
543	!-----
544	DO i=1,npoi
545	IF (iwet(i) == 1) then
546	precip(i) = -z(i)log(rn3)beta(i)
547	ENDIF
548	ENDDO
549	!-----
550	!---- method ii --
551	!-----
552	!---- here we use a one-parameter Weibull distribution function
553	!---- following the analysis of Selker and Haith (1990)
554	!-----
555	!---- Selker, J.S. and D.A. Haith, 1990: Development and testing
556	!---- of single- parameter precipitation distributions,
557	!---- Water Resources Research, 11, 2733-2740.
558	!-----
559	!---- this technique seems to have a significant advantage over other
560	!---- means of generating rainfall distribution functions
561	!-----
562	!---- by calibrating the Weibull function to U.S. precipitation records,
563	!---- Selker and Haith were able to establish the following relationship
564	!-----
565	!---- the cumulative probability of rainfall intensity x is given as:
566	!-----
567	!---- f(x) = 1.0-exp(-(1.191 x / rainpwd)**0.75)
568	!-----
569	!---- where x : rainfall intensity
570	!---- rainpwd : rainfall per wet day
571	!-----
572	!---- using transformation method, take uniform deviate and convert
573	!---- it to a random number weighted by the following Weibull function
574	!-----
575	!---- call random_number(rndnum)
576	!-----
577	!---- precip(i) = rainpwd / 1.191(-log(1.0-rndnum))*1.333333
578	!-----
579	!---- bound daily precipitation to "REAListic" range
580	!-----
581	DO i=1,npoi
582	IF (iwet(i) == 1) THEN
583	!---------
584	!-------- lower end is determined by definition of a 'wet day'
585	!-------- (at least 0.25 mm of total precipitation)
586	!---------
587	!-------- upper end is to prevent ibis from blowing up
588	!---------
589	precip(i) = MAX(precip(i), 0.25) ! min = 0.25 mm/day
590	precip(i) = MIN(precip(i),150.00) ! max = 150.00 mm/day
591	ENDIF
592	ENDDO
593	ELSE
594	!-----
595	!---- (2.2) preserving the monthly number of precip days
596	!---- and monthly precip
597	!-----
598	DO i=1,npoi
599	!-------
600	!------ Correction Nathalie. C'est abs(xinwet) < 32 qu'il faut tester
601	!------ et non pas abs(xinprec(i,imonth)) < 32.
602	!-------
603	!------ IF ( (iwet(i) == 1).and.(abs(xinprec(i,imonth)) < 32.) ) THEN
604	IF ( (iwet(i) == 1).and.(abs(xinwet(i,imonth)) < 32.) ) THEN
605	precip(i) = xinprec(i,imonth)*REAL(ndaypm(imonth)) &
606	& /REAL(NINT(MAX(1.,xinwet(i,imonth))))
607	ELSE
608	precip(i) = 0.
609	ENDIF
610	ENDDO
611	ENDIF
612	!---
613	!-- (3) estimate expected minimum and maximum temperatures
614	!---
615	DO i=1,npoi
616	!-----
617	!---- first determine the expected maximum and minimum temperatures
618	!---- (from climatological means) for this day of the year
619	!-----
620	!---- mean daily mean temperature (K)
621	tdm = xint(i,it1w)+dt*(xint(i,it2w)-xint(i,it1w))+zero_t
622	!---- mean daily temperature range (K)
623	trngm = xintrng(i,it1w)+dt*(xintrng(i,it2w)-xintrng(i,it1w))
624	!---- mean minimum and maximum temperatures
625	tmaxm = tdm+0.56*trngm
626	tminm = tdm-0.44*trngm
627	!-----
628	!---- modify maximum temperatures for wet and dry days
629	!-----
630	IF (pwet(i) /= 0.0) THEN
631	tmaxd = tmaxm+pwet(i)omtmaxtrngm
632	tmaxw = tmaxd- omtmax*trngm
633	ELSE
634	tmaxd = tmaxm
635	tmaxw = tmaxm
636	ENDIF
637	!-----
638	!---- set the 'expected' maximum and minimum temperatures for today
639	!-----
640	!---- note that the expected minimum temperatures are the same for
641	!---- both wet and dry days
642	!-----
643	if (iwet(i) == 0) tmaxe(i) = tmaxd
644	if (iwet(i) == 1) tmaxe(i) = tmaxw
645	!-----
646	tmine(i) = tminm
647	!-----
648	!---- estimate variability in minimum and maximum temperatures
649	!-----
650	!---- tmaxs : standard deviation in maximum temperature (K)
651	!---- tmins : standard deviation in minimum temperature (K)
652	!-----
653	!---- Regression is based on analysis of 2-m air temperature data
654	!---- from the NCEP/NCAR reanalysis (1958-1997) for 294 land points
655	!---- over central North America
656	!---- (24N-52N, 130W-60W, 0.5-degree resolution):
657	!---- Daily max and min temperatures were calculated for each
658	!---- land point from daily mean temperature and temperature range.
659	!---- Anomalies were calculated
660	!---- by subtracting similar max and min temperatures calculated from
661	!---- monthly mean temperature and range (interpolated to daily).
662	!---- The 40 years of anomalies were then binned by month
663	!---- and the standard deviation calculated for each month.
664	!---- The 294 x 12 standard deviations were then regressed
665	!---- against the 3528 long-term monthly mean temperatures.
666	!-----
667	!---- note: the values are bound to be greater than 1.0 K
668	!---- (at the very least they must be bound
669	!---- so they don't go below zero)
670	!-----
671	tmaxs(i) = MAX(1.0,-0.0713*(tdm-zero_t)+4.89)
672	tmins(i) = MAX(1.0,-0.1280*(tdm-zero_t)+5.73)
673	ENDDO
674	!---
675	!-- (4) estimate expected cloud cover
676	!---
677	!---
678	!-- the formulation of dry and wet cloud cover has been
679	!-- derived from the weather generator used in the epic crop model
680	!---
681	DO i=1,npoi
682	!-----
683	!---- cloudm : mean cloud cover for today
684	!---- cloudd : dry day cloud cover
685	!---- cloudw : dry day cloud cover
686	!---- cloude : expected cloud cover today
687	!-----
688	!---- mean cloud cover (%)
689	!-----
690	cloudm = xincld(i,it1w)+dt*(xincld(i,it2w)-xincld(i,it1w))
691	!-----
692	!---- convert from percent to fraction
693	!-----
694	cloudm = cloudm/100.0
695	!-----
696	!---- adjust cloud cover depending on dry day / rainy day
697	!---- following logic of the EPIC weather generator code
698	!-----
699	IF (pwet(i) /= 0.0) THEN
700	cloudd = (cloudm-pwet(i)omcloud)/(1.0-pwet(i)omcloud)
701	cloudd = MIN(1.0,MAX(0.0,cloudd))
702	cloudw = (cloudm-(1.0-pwet(i))*cloudd)/pwet(i)
703	ELSE
704	cloudd = cloudm
705	cloudw = cloudm
706	ENDIF
707	IF (iwet(i) == 0) cloude(i) = cloudd
708	IF (iwet(i) == 1) cloude(i) = cloudw
709	!-----
710	!---- estimate variability in cloud cover for wet and dry days
711	!---- following numbers proposed by Richardson
712	!-----
713	!---- clouds : standard deviation of cloud fraction
714	!-----
715	IF (iwet(i) == 0) clouds(i) = 0.24*cloude(i)
716	IF (iwet(i) == 1) clouds(i) = 0.48*cloude(i)
717	ENDDO
718	!---
719	! (5) determine today's temperatures and cloud cover using
720	! first-order serial autoregressive technique
721	!---
722	!-- use the Richardson (1981) weather generator approach to simulate the
723	!-- daily values of minimum / maximum temperature and cloud cover
724	!---
725	!-- following the implementation of the Richardson WGEN weather
726	!-- generator used in the EPIC crop model
727	!---
728	!-- this approach uses a multivariate generator, which assumes that the
729	!-- perturbation of minimum / maximum temperature and cloud cover are
730	!-- normally distributed and that the serial correlation of each
731	!-- variable may be described by a first-order autoregressive model
732	!---
733	!-- generate standard deviates for weather generator
734	!---
735	DO j=1,3
736	ee(j) = 9999.
737	DO WHILE (ee(j) > 2.5)
738	IF (is_root_prc) THEN
739	CALL random_number (rn1)
740	CALL random_number (rn2)
741	ENDIF
742	CALL bcast(rn1)
743	CALL bcast(rn2)
744	ee(j) = SQRT(-2.0LOG(rn1))COS(6.283185*rn2)
745	ENDDO
746	ENDDO
747	!---
748	!-- zero out vectors
749	!---
750	r(1:3) = 0.0
751	rr(1:npoi,1:3) = 0.0
752	!---
753	!-- update working vectors
754	!---
755	do j=1,3
756	do k=1,3
757	r(j) = r(j)+b(j,k)*ee(j)
758	enddo
759	enddo
760	!---
761	do j=1,3
762	do k=1,3
763	do i=1,npoi
764	rr(i,j) = rr(i,j)+a(j,k)*xstore(i,k)
765	enddo
766	enddo
767	enddo
768	!---
769	!-- solve for x() perturbation vector and save current vector
770	!-- into the xim1() storage vector (saved for each point)
771	!---
772	do j=1,3
773	do i=1,npoi
774	xstore(i,j) = r(j)+rr(i,j)
775	enddo
776	enddo
777	!---
778	!-- determine today's minimum and maximum temperature
779	!--
780	do i=1,npoi
781	tmax(i) = tmaxe(i)+tmaxs(i)*xstore(i,1)
782	tmin(i) = tmine(i)+tmins(i)*xstore(i,2)
783	!-----
784	!---- if tmin > tmax, then switch the two around
785	!-----
786	if (tmin(i) > tmax(i)) then
787	tdum = tmax(i)
788	tmax(i) = tmin(i)
789	tmin(i) = tdum
790	ENDIF
791	!---- daily average temperature
792	td(i) = 0.44tmax(i)+0.56tmin(i)
793	!---- determine today's cloud cover
794	cloud(i) = cloude(i)+clouds(i)*xstore(i,3)
795	!---- constrain cloud cover to be between 0 and 100%
796	cloud(i) = MAX(0.0,MIN(1.0,cloud(i)))
797	enddo
798	!---
799	!-- (6) estimate today's surface atmospheric pressure
800	!---
801	do i=1,npoi
802	!-----
803	!---- simply a function of the daily average temperature
804	!---- and topographic height -- nothing fancy here
805	!-----
806	psurf(i) = 101325.0(td(i)/(td(i)+0.0065xintopo(i)))**rwork
807	enddo
808	!---
809	!-- (7) estimate today's relative humidity
810	!---
811	IF (is_root_prc) THEN
812	CALL random_number (rn)
813	ENDIF
814	CALL bcast(rn)
815	!---
816	CALL weathgen_qsat (npoi,td,psurf,qsattd)
817	!---
818	do i=1,npoi
819	!-----
820	!---- the formulation of dry and wet relative humidities has been
821	!---- derived from the weather generator used in the epic crop model
822	!-----
823	!---- qdm : mean relative humidity
824	!---- qdd : dry day relative humidity
825	!---- qdw : rainy day relative humidity
826	!---- qde : expected relative humidity (based on wet/dry decision)
827	!-----
828	!---- mean relative humidity (%)
829	qdm(i) = xinq(i,it1w)+dt*(xinq(i,it2w)-xinq(i,it1w))
830	!---- convert from percent to fraction
831	qdm(i) = qdm(i)/100.0
832	!-----
833	!---- adjust humidity depending on dry day / rainy day
834	!---- following logic of the EPIC weather generator code
835	!-----
836	if (pwet(i) /= 0.0) then
837	qdd(i) = (qdm(i)-pwet(i)omqd)/(1.0-pwet(i)omqd)
838	if (qdd(i) < 0.2) then
839	qdd(i) = 0.2
840	if (qdd(i) > qdm(i)) qdm(i) = qdd(i)
841	ENDIF
842	qdd(i) = MIN(1.0,qdd(i))
843	qdw(i) = (qdm(i)-(1.0-pwet(i))*qdd(i))/pwet(i)
844	ELSE
845	qdd(i) = qdm(i)
846	qdw(i) = qdm(i)
847	ENDIF
848	!-----
849	if (iwet(i) == 0) qde(i) = qdd(i)
850	if (iwet(i) == 1) qde(i) = qdw(i)
851	!-----
852	!---- estimate lower and upper bounds of humidity distribution function
853	!---- following logic of the EPIC weather generator code
854	!-----
855	xx = exp(qde(i))
856	qdup(i) = qde(i)+(1.0-qde(i))*xx/e
857	qdlow(i) = qde(i)*(1.0-1./xx)
858	!-----
859	!---- randomlly select humidity from triangular distribution function
860	!---- following logic of the EPIC weather generator code
861	!-----
862	!---- call random_number(rn) ! GK done before
863	!-----
864	y = 2.0/(qdup(i)-qdlow(i))
865	!-----
866	b3 = qde(i)-qdlow(i)
867	b2 = qdup(i)-qde(i)
868	b1 = rn/y
869	!-----
870	x1 = y*b3/2.0
871	!-----
872	if (rn > x1) then
873	qd(i) = qdup(i)-sqrt (b2b2-2.0b2(b1-0.5b3))
874	ELSE
875	qd(i) = qdlow(i)+sqrt (2.0b1b3)
876	ENDIF
877	!-----
878	!---- adjust daily humidity to conserve monthly mean values
879	!-----
880	!---- note that this adjustment sometimes gives rise to humidity
881	!---- values greater than 1.0 -- which is corrected below
882	!-----
883	amn = (qdup(i)+qde(i)+qdlow(i))/3.0
884	qd(i) = qd(i)*qde(i)/amn
885	!-----
886	!---- constrain daily average relative humidity
887	!-----
888	qd(i) = MAX(0.30,MIN(qd(i),0.99))
889	!-----
890	!---- convert from relative humidity to specific humidity at
891	!---- daily mean temperature
892	!-----
893	qd(i) = qd(i)*qsattd(i)
894	enddo
895	!---
896	!-- (8) estimate today's daily average wind speed
897	!---
898	IF (is_root_prc) THEN
899	CALL random_number (rn4)
900	ENDIF
901	CALL bcast(rn4)
902
903	do i=1,npoi
904	!-----
905	!---- first estimate the expected daily average wind speed (from monthly
906	!---- means)
907	!-----
908	eud = xinwind(i,it1w)+dt*(xinwind(i,it2w)-xinwind(i,it1w))
909	!-----
910	!---- following logic of the EPIC weather generator
911	!---- select random wind speed following this equation
912	!-----
913	!---- call random_number(rn4)
914	!-----
915	ud(i) = 1.13989eud(-log(rn4))**0.30
916	!---- constrain daily wind speeds to be between 2.5 and 10.0 m/sec
917	ud(i) = MAX(2.5,MIN(ud(i),10.0))
918	ENDDO
919	ELSE
920	!---
921	!-- use REAL daily climate data
922	!---
923	DO i=1,npoi
924	!-----
925	!---- use basic daily climate data, converting units
926	!-----
927	!---- daily total precipitation
928	precip(i) = xinprec(i,imonth)
929	!---- daily average temperatures
930	td(i) = xint(i,imonth)+zero_t
931	trngm = MIN(44.0,xintrng(i,imonth))
932	!-----
933	tmax(i) = td(i)+0.56*trngm
934	tmin(i) = td(i)-0.44*trngm
935	!---- daily average cloud cover
936	cloud(i) = xincld(i,imonth)*0.01
937	!---- daily average specific humidity
938	qd(i) = xinq(i,imonth)
939	!---- daily average wind speed
940	ud(i) = xinwind(i,imonth)
941	!-----
942	!---- compute surface atmospheric pressure
943	!-----
944	psurf(i) = 101325.0(td(i)/(td(i)+0.0065xintopo(i)))**rwork
945	ENDDO
946	ENDIF
947	!-------------------
948	END SUBROUTINE daily
949	!-
950	!===
951	!-
952	SUBROUTINE diurnal &
953	& (npoi, nband, time, jday, plens, startp, endp, latitude, &
954	& cloud, tmax, tmin, precip, qd, ud, psurf, &
955	& fira, solad, solai, ua, ta, qa, raina, snowa, rh)
956	!---------------------------------------------------------------------
957	IMPLICIT NONE
958	!-
959	! input
960	!-
961	! number of grid points
962	INTEGER,INTENT(IN) :: npoi
963	! number of visible bands
964	INTEGER,INTENT(IN) :: nband
965	REAL,INTENT(IN) :: time
966	INTEGER, INTENT(IN) :: jday
967	REAL,INTENT(IN) :: plens,startp,endp
968	! latitude in degrees
969	REAL,INTENT(IN) :: latitude(npoi)
970	! cloud fraction [0,1]
971	REAL,INTENT(IN) :: cloud(npoi)
972	! maximum daily temperature (K)
973	REAL,INTENT(IN) :: tmax(npoi)
974	! maximum daily temperature (K)
975	REAL,INTENT(IN) :: tmin(npoi)
976	! daily precitation (mm/day)
977	REAL,INTENT(IN) :: precip(npoi)
978	! daily average specific humidity (kg_h2o/kg_air)
979	REAL,INTENT(IN) :: qd(npoi)
980	! daily average wind speed (m/sec)
981	REAL,INTENT(IN) :: ud(npoi)
982	! surface pressure (Pa)
983	REAL,INTENT(IN) :: psurf(npoi)
984	!-
985	! output
986	!-
987	! incoming ir flux (W m-2)
988	REAL,INTENT(OUT) :: fira(npoi)
989	! direct downward solar flux (W m-2)
990	REAL,INTENT(OUT) :: solad(npoi,nband)
991	! diffuse downward solar flux (W m-2)
992	REAL,INTENT(OUT) :: solai(npoi,nband)
993	! wind speed (m s-1)
994	REAL,INTENT(OUT) :: ua(npoi)
995	! air temperature (K)
996	REAL,INTENT(OUT) :: ta(npoi)
997	! specific humidity (kg_h2o/kg_air)
998	REAL,INTENT(OUT) :: qa(npoi)
999	! rainfall rate (mm/day)
1000	REAL,INTENT(OUT) :: raina(npoi)
1001	! snowfall rate (mm/day)
1002	REAL,INTENT(OUT) :: snowa(npoi)
1003	! relative humidity(%)
1004	REAL,INTENT(OUT) :: rh(npoi)
1005	!-
1006	! local
1007	!-
1008	REAL,PARAMETER :: stef = 5.67051E-8
1009	REAL,PARAMETER :: pi = 3.1415927
1010	REAL,SAVE :: step
1011	REAL,PARAMETER :: pir = pi/180.
1012	REAL :: xl,so,xllp,xee,xse
1013	REAL :: xlam,dlamm,anm,ranm,ranv,anv,tls,rlam
1014	REAL :: sd,cd,deltar,delta,Dis_ST,ddt
1015	!-
1016	REAL :: coszen(npoi) ! cosine of solar zenith angle
1017	REAL :: rtime
1018	REAL :: orbit,angle,xdecl,xlat
1019	REAL :: sw,frac,gamma,qmin,qmax,qsa,emb,ea,ec,dtair,dtcloud,rn
1020	REAL :: trans(npoi), fdiffuse(npoi), qsatta(npoi), qsattmin(npoi)
1021	INTEGER :: i,ib
1022	INTEGER,SAVE :: npoi0
1023	LOGICAL,SAVE :: firstcall = .TRUE.
1024	!---------------------------------------------------------------------
1025	! GK240100
1026	IF (firstcall) THEN
1027	IF ( TRIM(calendar_str) .EQ. 'gregorian' ) THEN
1028	step = 1.0
1029	ELSE
1030	step = one_year/365.2425
1031	ENDIF
1032	firstcall = .FALSE.
1033	npoi0 = npoi
1034	ELSE IF (npoi /= npoi0) THEN
1035	WRITE(numout,*) 'Domain size old, new: ',npoi0,npoi
1036	STOP 'WG Diurnal: Problem: Domain has changed since last call'
1037	ENDIF
1038	!-
1039	! calendar and orbital calculations
1040	!-
1041	! calculate time in hours
1042	rtime = time/3600.0
1043	!-
1044	! calculate the earth's orbital angle (around the sun) in radians
1045	orbit = 2.0piREAL(jday)/365.2425
1046	!-
1047	! calculate the solar hour angle in radians
1048	angle = 2.0pi(rtime-12.0)/24.0
1049	!-
1050	! calculate the current solar declination angle
1051	! ref: global physical climatology, hartmann, appendix a
1052	!
1053	! xdecl = 0.006918 &
1054	! -0.399912*cos(orbit) &
1055	! +0.070257*sin(orbit) &
1056	! -0.006758cos(2.0orbit) &
1057	! +0.000907sin(2.0orbit) &
1058	! -0.002697cos(3.0orbit) &
1059	! +0.001480sin(3.0orbit)
1060	!
1061	! calculate the effective solar constant,
1062	! including effects of eccentricity
1063	! ref: global physical climatology, hartmann, appendix a
1064	!
1065	! sw = 1370.*( 1.000110 &
1066	! +0.034221*cos(orbit) &
1067	! +0.001280*sin(orbit) &
1068	! +0.000719cos(2.0orbit) &
1069	! +0.000077sin(2.0orbit))
1070	!
1071	! correction Marie-France Loutre
1072	!
1073	! orbital parameters
1074	!
1075	! ecc = 0.016724
1076	! perh = 102.04
1077	! xob = 23.446
1078	!-
1079	xl = perh+180.0
1080	! so : sinus of obliquity
1081	so = sin(xob*pir)
1082	!-
1083	xllp = xl*pir
1084	xee = ecc*ecc
1085	xse = sqrt(1.0d0-xee)
1086	! xlam : true long. sun for mean long. = 0
1087	xlam = (ecc/2.0+eccxee/8.0d0)(1.0+xse)*sin(xllp)-xee/4.0 &
1088	& (0.5+xse)sin(2.0xllp)+eccxee/8.0*(1.0/3.0+xse) &
1089	& sin(3.0xllp)
1090	xlam =2.0d0*xlam/pir
1091	! dlamm : mean long. sun for ma-ja
1092	dlamm =xlam+(jday-79)*step
1093	anm = dlamm-xl
1094	ranm = anm*pir
1095	xee = xee*ecc
1096	! ranv : anomalie vraie (radian)
1097	ranv = ranm+(2.0ecc-xee/4.0)sin(ranm)+5.0/4.0eccecc &
1098	& sin(2.0ranm)+13.0/12.0xeesin(3.0*ranm)
1099	! anv : anomalie vraie (degrees)
1100	anv = ranv/pir
1101	! tls : longitude vraie (degrees)
1102	tls = anv+xl
1103	! rlam : longitude vraie (radian)
1104	rlam = tls*pir
1105	! sd and cd: cosinus and sinus of solar declination angle (delta)
1106	! sinus delta = sin (obl)*sin(lambda) with lambda = real longitude
1107	! (Phd. thesis of Marie-France Loutre, ASTR-UCL, Belgium, 1993)
1108	sd = so*sin(rlam)
1109	cd = sqrt(1.0d0-sd*sd)
1110	! delta : Solar Declination (degrees, angle sun at equator)
1111	deltar = atan(sd/cd)
1112	delta = deltar/pir
1113	!-
1114	! Eccentricity Effect
1115	!-
1116	Dis_ST =(1.0-eccecc)/(1.0+ecccos(ranv))
1117	! ddt : 1 / normalized earth's sun distance
1118	ddt = 1.0/Dis_ST
1119	!-
1120	! Insolation normal to the atmosphere (W/m2)
1121	!-
1122	sw = ddt ddt 1370.d0
1123	!-
1124	xdecl = deltar
1125	!-
1126	! solar calculations
1127	!-
1128	do i=1,npoi
1129	!---
1130	!-- calculate the latitude in radians
1131	!---
1132	xlat = latitude(i)*pi/180.0
1133	!---
1134	!-- calculate the cosine of the solar zenith angle
1135	!---
1136	coszen(i) = MAX(0.0, (sin(xlat)*sin(xdecl) &
1137	& + cos(xlat)cos(xdecl)cos(angle)))
1138	!---
1139	!-- calculate the solar transmission through the atmosphere
1140	!-- using simple linear function of tranmission and cloud cover
1141	!---
1142	!-- note that the 'cloud cover' data is typically obtained from
1143	!-- sunshine hours -- not direct cloud observations
1144	!---
1145	!-- where, cloud cover = 1 - sunshine fraction
1146	!---
1147	!-- different authors present different values for the slope and
1148	!-- intercept terms of this equation
1149	!---
1150	!-- Friend, A: Parameterization of a global daily weather generator
1151	!-- for terrestrial ecosystem and biogeochemical modelling,
1152	!-- Ecological Modelling
1153	!---
1154	!-- Spitters et al., 1986: Separating the diffuse and direct component
1155	!-- of global radiation and its implications for modeling canopy
1156	!-- photosynthesis, Part I: Components of incoming radiation,
1157	!-- Agricultural and Forest Meteorology, 38, 217-229.
1158	!---
1159	!-- A. Friend : trans = 0.251+0.509*(1.0-cloud(i))
1160	!-- Spitters et al. : trans = 0.200+0.560*(1.0-cloud(i))
1161	!---
1162	!-- we are using the values from A. Friend
1163	!---
1164	trans(i) = 0.251+0.509*(1.0-cloud(i))
1165	!---
1166	!-- calculate the fraction of indirect (diffuse) solar radiation
1167	!-- based upon the cloud cover
1168	!---
1169	!-- note that these relationships typically are measured for either
1170	!-- monthly or daily timescales, and may not be exactly appropriate
1171	!-- for hourly calculations -- however, in ibis, cloud cover is fixed
1172	!-- through the entire day so it may not make much difference
1173	!---
1174	!-- method i --
1175	!---
1176	!-- we use a simple empirical relationships
1177	!-- from Nikolov and Zeller (1992)
1178	!---
1179	!-- Nikolov, N. and K.F. Zeller, 1992: A solar radiation algorithm
1180	!-- for ecosystem dynamics models, Ecological Modelling, 61, 149-168.
1181	!---
1182	fdiffuse(i) = 1.0045+trans(i)*( 0.0435+trans(i) &
1183	& (-3.5227+trans(i)2.6313))
1184	!---
1185	IF (trans(i) > 0.75) fdiffuse(i) = 0.166
1186	!---
1187	!-- method ii --
1188	!---
1189	!-- another method was suggested by Spitters et al. (1986) based on
1190	!-- long-term data from the Netherlands
1191	!--
1192	!-- Spitters et al., 1986: Separating the diffuse and direct component
1193	!-- of global radiation and its implications for modeling canopy
1194	!-- photosynthesis, Part I: Components of incoming radiation,
1195	!-- Agricultural and Forest Meteorology, 38, 217-229.
1196	!---
1197	!-- if ((trans == 0.00).and.(trans < 0.07)) then
1198	!-- fdiffuse = 1.0
1199	!-- else if ((trans >= 0.07).and.(trans < 0.35)) then
1200	!-- fdiffuse = 1.0-2.3(trans-0.07)*2
1201	!-- else if ((trans >= 0.35).and.(trans < 0.75)) then
1202	!-- fdiffuse = 1.33-1.46*trans
1203	!-- ELSE
1204	!-- fdiffuse = 0.23
1205	!-- endif
1206	!---
1207	ENDDO
1208	!-
1209	! do for each waveband
1210	!-
1211	DO ib=1,nband
1212	!---
1213	!-- calculate the fraction in each waveband
1214	!---
1215	!-- GK010200
1216	IF (nband == 2) then
1217	frac = 0.46+0.08*REAL(ib-1)
1218	ELSE
1219	frac = 1./REAL(nband)
1220	ENDIF
1221	!---
1222	DO i=1,npoi
1223	!-----
1224	!---- calculate the direct and indirect solar radiation
1225	!-----
1226	solad(i,ib) = swcoszen(i)fractrans(i)(1.-fdiffuse(i))
1227	solai(i,ib) = swcoszen(i)fractrans(i)fdiffuse(i)
1228	ENDDO
1229	ENDDO
1230	!-
1231	! temperature calculations
1232	!-
1233	!---
1234	!-- assign hourly temperatures using tmax and tmin
1235	!-- following Environmental Biophysics, by Campbell and Norman, p.23
1236	!---
1237	!-- this function fits a fourier series to the diurnal temperature cycle
1238	!-- note that the maximum temperature occurs at 2:00 pm local solar time
1239	!---
1240	!-- note that the daily mean value of gamma is 0.44,
1241	!-- so td = 0.44tmax+0.56tmin, instead of
1242	!-- td = 0.50tmax+0.50tmin
1243	!---
1244	gamma = 0.44-0.46SIN( pi/12.0rtime+0.9) &
1245	+0.11SIN(2.0pi/12.0*rtime+0.9)
1246	DO i=1,npoi
1247	ta(i) = tmax(i)gamma+tmin(i)(1.0-gamma)
1248	ENDDO
1249	!-
1250	! humidity calculations
1251	!-
1252	CALL weathgen_qsat (npoi,tmin,psurf,qsattmin)
1253	CALL weathgen_qsat (npoi,ta,psurf,qsatta)
1254	!-
1255	DO i=1,npoi
1256	!---
1257	!-- adjust specific humidity against daily minimum temperatures
1258	!---
1259	!-- To do this, qa is written as an approximate sine function
1260	!-- (same as ta) to preserve the daily mean specific humidity,
1261	!-- while also preventing rh from exceeding 99% at night
1262	!---
1263	!-- Note that the daily mean RH is not preserved, and therefore the
1264	!-- output RH will be slightly different from what was read in.
1265	!---
1266	!-- first adjust so that maximum RH cannot exceed 99% at night
1267	!---
1268	qmin = MIN(qd(i),0.99*qsattmin(i))
1269	qmax = (qd(i)-0.56*qmin)/0.44
1270	!---
1271	!-- if needed, adjust again to 99% at other times of the day (in which
1272	!-- case the daily mean specific humidity is also not preserved)
1273	!---
1274	qsa = 0.99*qsatta(i)
1275	!---
1276	!-- calculate the hourly specific humidity, using the above adjustments
1277	!---
1278	qa(i) = MIN(qsa,qmaxgamma+qmin(1.0-gamma))
1279	!---
1280	!-- calculate the hourly relative humidity
1281	!--
1282	rh(i) = 100.*qa(i)/qsatta(i)
1283	ENDDO
1284	!-
1285	! wind speed calculations
1286	!-
1287	IF (is_root_prc) THEN
1288	CALL random_number (rn)
1289	ENDIF
1290	CALL bcast(rn)
1291	!-
1292	DO i=1,npoi
1293	!---
1294	!-- following logic of the EPIC weather generator
1295	!-- select random wind speed following this equation
1296	!---
1297	!-- call random_number(rn) ! done before
1298	!---
1299	ua(i) = 1.13989ud(i)(-log(rn))**0.30
1300	!---
1301	!-- fix wind speeds to always be above 2.5 m/sec and below 10.0 m/sec
1302	!---
1303	ua(i) = MAX(2.5,MIN(10.0,ua(i)))
1304	ENDDO
1305	!-
1306	! ir flux calculations
1307	!-
1308	DO i=1,npoi
1309	!---
1310	!-- clear-sky emissivity as a function of water vapor pressure
1311	!-- and atmospheric temperature
1312	!---
1313	!-- calculate the ir emissivity of the clear sky
1314	!-- using equation from idso (1981) water resources res., 17, 295-304
1315	!---
1316	emb = 0.01(psurf(i)qa(i)/(0.622+qa(i)))
1317	ea = 0.70+5.95e-5embEXP(1500.0/ta(i))
1318	!---
1319	!-- assign the ir emissivity of clouds
1320	!-- (assume they are ~black in the ir)
1321	!---
1322	ec = 0.950
1323	!---
1324	!-- assign the temperature difference of emissions (air+cloud) from
1325	!-- the surface air temperature
1326	!---
1327	dtair = 0.0
1328	dtcloud = 0.0
1329	!---
1330	!-- total downward ir is equal to the sum of:
1331	!---
1332	!-- (1) clear sky contribution to downward ir radiation flux
1333	!-- (2) cloud contribution to downward ir radiation flux
1334	!---
1335	fira(i) = (1.-cloud(i))eastef(ta(i)-dtair)*4 &
1336	& +cloud(i)ecstef(ta(i)-dtcloud)*4
1337	ENDDO
1338	!-
1339	! snow and rain calculations
1340	!-
1341	DO i=1,npoi
1342	!---
1343	!-- reset snow and rain to zero
1344	!---
1345	snowa(i) = 0.0
1346	raina(i) = 0.0
1347	!---
1348	!-- if precipitation event then calculate
1349	!---
1350	IF (time >= startp .and. time < endp) then
1351	!-----
1352	!---- for rain / snow partitioning, make it all rain if
1353	!---- ta > 2.5 C and all snow if ta <= 2.5 C
1354	!-----
1355	!---- reference:
1356	!-----
1357	!---- Auer, A. H., 1974: The rain versus snow threshold temperatures,
1358	!---- Weatherwise, 27, 67.
1359	!-----
1360	IF (ta(i)-273.15 > 2.5) then
1361	raina(i) = precip(i)/plens
1362	ELSE
1363	snowa(i) = precip(i)/plens
1364	ENDIF
1365	ENDIF
1366	ENDDO
1367	!---------------------
1368	END SUBROUTINE diurnal
1369	!-
1370	!===
1371	!-
1372	SUBROUTINE weathgen_main &
1373	& (itau, istp, filename, force_id, iim, jjm, nband, &
1374	& rest_id, lrstread, lrstwrite, &
1375	& limit_west, limit_east, limit_north, limit_south, &
1376	& zonal_res, merid_res, lon, lat, date0, dt_force, &
1377	& kindex, nbindex, &
1378	& swdown, rainf, snowf, tair, u, v, qair, pb, lwdown)
1379	!---------------------------------------------------------------------
1380	IMPLICIT NONE
1381	!-
1382	INTEGER,INTENT(IN) :: itau,istp
1383	CHARACTER(LEN=*),INTENT(IN) :: filename
1384	INTEGER,INTENT(IN) :: force_id
1385	INTEGER,INTENT(IN) :: iim,jjm
1386	INTEGER,INTENT(IN) :: nband
1387	INTEGER,INTENT(IN) :: rest_id
1388	LOGICAL,INTENT(IN) :: lrstread, lrstwrite
1389	REAL,INTENT(IN) :: limit_west,limit_east
1390	REAL,INTENT(IN) :: limit_north,limit_south
1391	REAL,INTENT(IN) :: zonal_res,merid_res
1392	REAL,DIMENSION(iim,jjm),INTENT(IN) :: lon,lat
1393	REAL,INTENT(IN) :: date0,dt_force
1394	!-
1395	INTEGER,DIMENSION(iim*jjm),INTENT(INOUT) :: kindex
1396	INTEGER,INTENT(INOUT) :: nbindex
1397	!-
1398	REAL,DIMENSION(iim,jjm),INTENT(OUT) :: &
1399	& tair,pb,qair,swdown,rainf,snowf, u,v,lwdown
1400	REAL, ALLOCATABLE, DIMENSION(:,:) :: &
1401	& tair_g,pb_g,qair_g,swdown_g,rainf_g,snowf_g, u_g,v_g,lwdown_g
1402	REAL,DIMENSION(nbindex) :: &
1403	& tairl,pbl,qairl,swdownl,rainfl,snowfl, ul,vl,lwdownl
1404	INTEGER :: i,j,ij
1405	!---------------------------------------------------------------------
1406	IF (lrstread) THEN
1407	CALL weathgen_begin ( &
1408	& dt_force,itau, date0, &
1409	& rest_id,iim,jjm, &
1410	& lon,lat,tair,pb,qair,kindex,nbindex)
1411	RETURN
1412	ELSE
1413	CALL weathgen_get &
1414	& (itau, date0, dt_force, nbindex, nband, lat_land, &
1415	& swdownl, rainfl, snowfl, tairl, ul, vl, qairl, pbl, lwdownl)
1416
1417	tair(:,:)=val_exp
1418	qair(:,:)=val_exp
1419	pb(:,:)=val_exp
1420	!!$ tair(:,:)=280.
1421	!!$ qair(:,:)=1.E-03
1422	!!$ pb(:,:)=101325
1423	DO ij=1,nbindex
1424	j = (((kindex(ij)-1)/iim) + 1)
1425	i = (kindex(ij) - (j-1)*iim)
1426	!
1427	swdown(i,j)=swdownl(ij)
1428	rainf(i,j)=rainfl(ij)
1429	snowf(i,j)=snowfl(ij)
1430	tair(i,j)=tairl(ij)
1431	u(i,j)=ul(ij)
1432	v(i,j)=vl(ij)
1433	qair(i,j)=qairl(ij)
1434	pb(i,j)=pbl(ij)
1435	lwdown(i,j)=lwdownl(ij)
1436	ENDDO
1437	!---
1438	IF (dump_weather) THEN
1439	ALLOCATE(tair_g(iim_g,jjm_g))
1440	ALLOCATE(pb_g(iim_g,jjm_g))
1441	ALLOCATE(qair_g(iim_g,jjm_g))
1442	ALLOCATE(swdown_g(iim_g,jjm_g))
1443	ALLOCATE(rainf_g(iim_g,jjm_g))
1444	ALLOCATE(snowf_g(iim_g,jjm_g))
1445	ALLOCATE(u_g(iim_g,jjm_g))
1446	ALLOCATE(v_g(iim_g,jjm_g))
1447	ALLOCATE(lwdown_g(iim_g,jjm_g))
1448
1449	CALL gather2D(tair, tair_g)
1450	CALL gather2D(pb, pb_g)
1451	CALL gather2D(qair, qair_g)
1452	CALL gather2D(swdown, swdown_g)
1453	CALL gather2D(rainf, rainf_g)
1454	CALL gather2D(snowf, snowf_g)
1455	CALL gather2D(u, u_g)
1456	CALL gather2D(v, v_g)
1457	CALL gather2D(lwdown, lwdown_g)
1458	IF (is_root_prc) THEN
1459	CALL weathgen_dump &
1460	& (itau, dt_force, iim_g, jjm_g, nbp_glo, index_g, lrstwrite, &
1461	& swdown_g, rainf_g, snowf_g, tair_g, u_g, v_g, qair_g, pb_g, lwdown_g)
1462	ENDIF
1463	ENDIF
1464	!---
1465	IF (lrstwrite) THEN
1466	CALL weathgen_restwrite (rest_id,istp,iim,jjm,nbindex,kindex)
1467	ENDIF
1468	ENDIF
1469	!---------------------------
1470	END SUBROUTINE weathgen_main
1471	!-
1472	!===
1473	!-
1474	SUBROUTINE weathgen_init &
1475	& (filename,dt_force,force_id,iim,jjm, &
1476	& zonal_res,merid_res,lon,lat,kindex,nbindex,iind,jind)
1477	!---------------------------------------------------------------------
1478	IMPLICIT NONE
1479	!-
1480	CHARACTER(LEN=*),INTENT(IN) :: filename
1481	REAL,INTENT(IN) :: dt_force
1482	INTEGER,INTENT(IN) :: force_id
1483	INTEGER,INTENT(IN) :: iim, jjm
1484	REAL,INTENT(IN) :: zonal_res,merid_res
1485	REAL,DIMENSION(iim,jjm),INTENT(IN) :: lon,lat
1486	!-
1487	INTEGER,DIMENSION(iim*jjm),INTENT(OUT) :: kindex
1488	INTEGER,INTENT(OUT) :: nbindex
1489	INTEGER,DIMENSION(iim),INTENT(OUT) :: iind
1490	INTEGER,DIMENSION(jjm),INTENT(OUT) :: jind
1491	!-
1492	REAL,PARAMETER :: fcrit = .5
1493	REAL,DIMENSION(:),ALLOCATABLE :: lon_file, lon_temp
1494	REAL,DIMENSION(:,:),ALLOCATABLE :: nav_lon, nav_lat
1495	REAL,DIMENSION(:),ALLOCATABLE :: lat_file, lat_temp
1496	REAL,DIMENSION(:,:),ALLOCATABLE :: lsm_file
1497	REAL :: xextent_file, yextent_file, xres_file, yres_file
1498	INTEGER :: i, j, plotstep
1499	REAL,DIMENSION(iim,jjm) :: mask
1500	CHARACTER(LEN=1),DIMENSION(0:1) :: maskchar
1501	CHARACTER(LEN=30) :: var_name
1502	REAL :: x_cut
1503
1504	REAL,DIMENSION(iim) :: tmp_lon
1505	REAL,DIMENSION(jjm) :: tmp_lat
1506	!---------------------------------------------------------------------
1507	!-
1508	! 0. messages, initialisations
1509	!-
1510	WRITE(numout,*) &
1511	& 'weathgen_init: Minimum land fraction on original grid =',fcrit
1512	CALL ioget_calendar(calendar_str)
1513	!-
1514	! 1. on lit les longitudes et latitudes de la grille de depart
1515	! ainsi que le masque terre-ocean
1516	!-
1517	CALL flinclo(force_id)
1518	CALL flininfo (filename,iim_file,jjm_file,llm_file,ttm_file,force_id)
1519	!-
1520	ALLOC_ERR=-1
1521	ALLOCATE(nav_lon(iim_file,jjm_file), STAT=ALLOC_ERR)
1522	IF (ALLOC_ERR/=0) THEN
1523	WRITE(numout,*) "ERROR IN ALLOCATION of nav_lon : ",ALLOC_ERR
1524	STOP
1525	ENDIF
1526
1527	ALLOC_ERR=-1
1528	ALLOCATE(nav_lat(iim_file,jjm_file), STAT=ALLOC_ERR)
1529	IF (ALLOC_ERR/=0) THEN
1530	WRITE(numout,*) "ERROR IN ALLOCATION of nav_lat : ",ALLOC_ERR
1531	STOP
1532	ENDIF
1533	!-
1534	var_name='nav_lon'
1535	CALL flinget (force_id,var_name,iim_file,jjm_file,0,0,1,1,nav_lon)
1536	var_name='nav_lat'
1537	CALL flinget (force_id,var_name,iim_file,jjm_file,0,0,1,1,nav_lat)
1538	!-
1539
1540	ALLOC_ERR=-1
1541	ALLOCATE(lon_file(iim_file), STAT=ALLOC_ERR)
1542	IF (ALLOC_ERR/=0) THEN
1543	WRITE(numout,*) "ERROR IN ALLOCATION of lon_file : ",ALLOC_ERR
1544	STOP
1545	ENDIF
1546
1547	ALLOC_ERR=-1
1548	ALLOCATE(lat_file(jjm_file), STAT=ALLOC_ERR)
1549	IF (ALLOC_ERR/=0) THEN
1550	WRITE(numout,*) "ERROR IN ALLOCATION of lat_file : ",ALLOC_ERR
1551	STOP
1552	ENDIF
1553	!-
1554	DO i=1,iim_file
1555	lon_file(i) = nav_lon(i,1)
1556	ENDDO
1557	DO j=1,jjm_file
1558	lat_file(j) = nav_lat(1,j)
1559	ENDDO
1560	!-
1561
1562	ALLOC_ERR=-1
1563	ALLOCATE(lsm_file(iim_file,jjm_file), STAT=ALLOC_ERR)
1564	IF (ALLOC_ERR/=0) THEN
1565	WRITE(numout,*) "ERROR IN ALLOCATION of lsm_file : ",ALLOC_ERR
1566	STOP
1567	ENDIF
1568	!-
1569	var_name='lsmera'
1570	CALL flinget (force_id,var_name,iim_file,jjm_file,0,0,1,1,lsm_file)
1571	!-
1572	! 2. La resolution du modele ne doit pas etre superieure
1573	! a celle de la grille de depart
1574	!-
1575	xextent_file = lon_file(iim_file)-lon_file(1)
1576	yextent_file = lat_file(1)-lat_file(jjm_file)
1577	xres_file = xextent_file/REAL(iim_file-1)
1578	yres_file = yextent_file/REAL(jjm_file-1)
1579	!-
1580	IF (xres_file > zonal_res) THEN
1581	WRITE(numout,*) 'Zonal resolution of model grid too fine.'
1582	WRITE(numout,*) 'Resolution of original data (deg): ', xres_file
1583	STOP
1584	ENDIF
1585	!-
1586	IF (yres_file > merid_res) THEN
1587	WRITE(numout,*) 'Meridional resolution of model grid too fine.'
1588	WRITE(numout,*) 'Resolution of original data (deg): ', yres_file
1589	STOP
1590	ENDIF
1591	!-
1592	! 3. On verifie la coherence des coordonnees de depart et d'arrivee.
1593	! Sinon, il faut deplacer une partie du monde (rien que ca).
1594	!-
1595	i_cut = 0
1596	!-
1597
1598	ALLOC_ERR=-1
1599	ALLOCATE(lon_temp(iim_file), STAT=ALLOC_ERR)
1600	IF (ALLOC_ERR/=0) THEN
1601	WRITE(numout,*) "ERROR IN ALLOCATION of lon_temp : ",ALLOC_ERR
1602	STOP
1603	ENDIF
1604
1605	ALLOC_ERR=-1
1606	ALLOCATE(lat_temp(jjm_file), STAT=ALLOC_ERR)
1607	IF (ALLOC_ERR/=0) THEN
1608	WRITE(numout,*) "ERROR IN ALLOCATION of lat_temp : ",ALLOC_ERR
1609	STOP
1610	ENDIF
1611	!-
1612	IF (lon(iim,1) > lon_file(iim_file)) THEN
1613	!-- A l'Est de la region d'interet, il n'y a pas de donnees
1614	!-- le bout a l'Ouest de la region d'interet est deplace de 360 deg
1615	!-- vers l'Est
1616	x_cut = lon(1,1)
1617	DO i=1,iim_file
1618	IF (lon_file(i) < x_cut) i_cut = i
1619	ENDDO
1620	IF ((i_cut < 1).OR.(i_cut >= iim_file)) THEN
1621	STOP 'Cannot find longitude for domain shift'
1622	ENDIF
1623	!---
1624	lon_temp(1:iim_file-i_cut-1) = lon_file(i_cut:iim_file)
1625	lon_temp(iim_file-i_cut:iim_file) = lon_file(1:i_cut+1)+360.
1626	lon_file(:) = lon_temp(:)
1627	ELSEIF (lon(1,1) < lon_file(1)) THEN
1628	!-- A l'Ouest de la region d'interet, il n'y a pas de donnees
1629	!-- le bout a l'Est de la region d'interet est deplace de 360 deg
1630	!-- vers l'Ouest
1631	x_cut = lon(iim,1)
1632	DO i=1,iim_file
1633	IF (lon_file(i) < x_cut) i_cut = i
1634	ENDDO
1635	IF ( ( i_cut < 1 ) .OR. ( i_cut >= iim_file ) ) THEN
1636	STOP 'Cannot find longitude for domain shift'
1637	ENDIF
1638	!---
1639	lon_temp(1:iim_file-i_cut-1) = lon_file(i_cut:iim_file) -360.
1640	lon_temp(iim_file-i_cut:iim_file) = lon_file(1:i_cut+1)
1641	lon_file(:) = lon_temp(:)
1642	ENDIF
1643	!-
1644	DEALLOCATE (lon_temp)
1645	DEALLOCATE (lat_temp)
1646	!-
1647	IF ( (lon(1,1) < lon_file(1)) &
1648	& .OR.( (lon(iim,1) > lon_file(iim_file)) &
1649	& .AND.(lon(iim,1) > lon_file(1)+360.) ) ) THEN
1650	WRITE(numout,*) lon(:,1)
1651	WRITE(numout,*)
1652	WRITE(numout,*) lon_file(:)
1653	STOP 'weathgen_init: cannot find coherent x-coordinates'
1654	ENDIF
1655	!-
1656	IF (i_cut /= 0) THEN
1657	CALL shift_field (iim_file,jjm_file,i_cut,lsm_file)
1658	ENDIF
1659	!-
1660	! 4. Verification
1661	!-
1662	WRITE(numout,*)
1663	WRITE(numout,*) 'Input File: (Shifted) Global Land-Sea Mask'
1664	WRITE(numout,*)
1665	maskchar(0) = '-'
1666	maskchar(1) = 'X'
1667	plotstep = INT(REAL(iim_file-1)/termcol)+1
1668	DO j=1,jjm_file,plotstep
1669	DO i=1,iim_file,plotstep
1670	WRITE(numout,'(a1,$)') maskchar(NINT(lsm_file(i,j)))
1671	ENDDO
1672	WRITE(numout,*)
1673	ENDDO
1674	WRITE(numout,*)
1675	!-
1676	! 5. Prepare interpolation from fine grid land points to model grid
1677	!-
1678	! 5.1 how many grid points of the original grid fall into one grid
1679	! box of the model grid?
1680	!-
1681	n_agg = NINT((zonal_res/xres_file*merid_res/yres_file )+1.)
1682	!-
1683	! au diable l'avarice !
1684	!-
1685	n_agg = n_agg*2
1686	!-
1687	! 5.2 Allocate arrays where we store information about which
1688	! (and how many) points on the original grid fall
1689	! into which box on the model grid
1690	!-
1691
1692	ALLOC_ERR=-1
1693	ALLOCATE(ncorr(iim,jjm), STAT=ALLOC_ERR)
1694	IF (ALLOC_ERR/=0) THEN
1695	WRITE(numout,*) "ERROR IN ALLOCATION of ncorr : ",ALLOC_ERR
1696	STOP
1697	ENDIF
1698
1699	ALLOC_ERR=-1
1700	ALLOCATE(icorr(iim,jjm,n_agg), STAT=ALLOC_ERR)
1701	IF (ALLOC_ERR/=0) THEN
1702	WRITE(numout,*) "ERROR IN ALLOCATION of icorr : ",ALLOC_ERR
1703	STOP
1704	ENDIF
1705
1706	ALLOC_ERR=-1
1707	ALLOCATE(jcorr(iim,jjm,n_agg), STAT=ALLOC_ERR)
1708	IF (ALLOC_ERR/=0) THEN
1709	WRITE(numout,*) "ERROR IN ALLOCATION of jcorr : ",ALLOC_ERR
1710	STOP
1711	ENDIF
1712	!-
1713	! 6. Land-Ocean Mask on model grid
1714	!-
1715	tmp_lon = lon(:,1)
1716	tmp_lat = lat(1,:)
1717
1718	CALL mask_c_o &
1719	& (iim_file, jjm_file, lon_file, lat_file, lsm_file, fcrit, &
1720	& iim, jjm, zonal_res, merid_res, n_agg, tmp_lon, tmp_lat, &
1721	! & iim, jjm, zonal_res, merid_res, n_agg, lon(:,1), lat(1,:), &
1722	& mask, ncorr, icorr, jcorr)
1723	!-
1724	WRITE(numout,*) 'Land-Sea Mask on Model Grid'
1725	WRITE(numout,*)
1726	plotstep = INT(REAL(iim-1)/termcol)+1
1727	DO j=1,jjm,plotstep
1728	DO i=1,iim,plotstep
1729	WRITE(numout,'(a1,$)') maskchar(NINT(mask(i,j)))
1730	ENDDO
1731	WRITE(numout,*)
1732	ENDDO
1733	WRITE(numout,*)
1734	!-
1735	! 7. kindex table.
1736	!-
1737	nbindex = 0
1738	DO j=1,jjm
1739	DO i=1,iim
1740	IF (NINT(mask(i,j)) == 1) THEN
1741	nbindex = nbindex+1
1742	kindex(nbindex) = (j-1)*iim+i
1743	ENDIF
1744	ENDDO
1745	ENDDO
1746	nbindex_w = nbindex
1747	ALLOC_ERR=-1
1748	ALLOCATE(kindex_w(nbindex_w), STAT=ALLOC_ERR)
1749	IF (ALLOC_ERR/=0) THEN
1750	WRITE(numout,*) "ERROR IN ALLOCATION of kindex_w : ",ALLOC_ERR
1751	STOP
1752	ENDIF
1753	kindex_w(:)=kindex(1:nbindex)
1754	!-
1755	IF ( nbindex == 0 ) THEN
1756	WRITE(numout,*) 'Couillon! On est au plein milieu de l''ocean.'
1757	STOP 'Ou est-ce un bug?'
1758	ELSE
1759	WRITE(numout,*) 'Number of continental points: ',nbindex
1760	ENDIF
1761
1762	!-
1763	! 10. clean up
1764	!-
1765	DEALLOCATE (lon_file)
1766	DEALLOCATE (lat_file)
1767	DEALLOCATE (lsm_file)
1768
1769	END SUBROUTINE weathgen_init
1770
1771	SUBROUTINE weathgen_read_file &
1772	& (force_id,iim,jjm)
1773	!---------------------------------------------------------------------
1774	IMPLICIT NONE
1775	!-
1776	INTEGER,INTENT(IN) :: force_id
1777	INTEGER,INTENT(IN) :: iim, jjm
1778	!-
1779	REAL,PARAMETER :: fcrit = .5
1780
1781	CHARACTER(LEN=30) :: var_name
1782
1783	INTEGER,DIMENSION(iim*jjm) :: kindex
1784	INTEGER :: nbindex
1785
1786	REAL,DIMENSION(iim*jjm) :: xchamp
1787	REAL,DIMENSION(iim*jjm,nmon) :: xchampm
1788
1789	kindex(:)=0
1790
1791	#ifdef CPP_PARA
1792	nbindex=nbp_loc
1793	CALL scatter(index_g,kindex)
1794	kindex(1:nbindex)=kindex(1:nbindex)-(jj_begin-1)*iim_g
1795	#else
1796	! Copy saved land points index
1797	nbindex = nbindex_w
1798	kindex(1:nbindex_w) = kindex_w(:)
1799	#endif
1800
1801	!-
1802
1803	ALLOC_ERR=-1
1804	ALLOCATE(lat_land(nbindex), STAT=ALLOC_ERR)
1805	IF (ALLOC_ERR/=0) THEN
1806	WRITE(numout,*) "ERROR IN ALLOCATION of lat_land : ",ALLOC_ERR
1807	STOP
1808	ENDIF
1809
1810	!-
1811	! 8 topography
1812	!-
1813
1814	ALLOC_ERR=-1
1815	ALLOCATE(xintopo(nbindex), STAT=ALLOC_ERR)
1816	IF (ALLOC_ERR/=0) THEN
1817	WRITE(numout,*) "ERROR IN ALLOCATION of xintopo : ",ALLOC_ERR
1818	STOP
1819	ENDIF
1820	!-
1821	var_name='altitude'
1822	CALL weather_read (force_id,var_name,iim_file,jjm_file,1,i_cut, &
1823	iim,jjm,n_agg,ncorr,icorr,jcorr,xchamp)
1824	xintopo(:)=xchamp(kindex(1:nbindex))
1825	!-
1826	! 9. Allocate and read the monthly fields
1827	!-
1828
1829	ALLOC_ERR=-1
1830	ALLOCATE(xinwet(nbindex,nmon), STAT=ALLOC_ERR)
1831	IF (ALLOC_ERR/=0) THEN
1832	WRITE(numout,*) "ERROR IN ALLOCATION of xinwet : ",ALLOC_ERR
1833	STOP
1834	ENDIF
1835	var_name='prs'
1836	CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, &
1837	iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm)
1838	xinwet(:,:)=xchampm(kindex(1:nbindex),:)
1839	!-
1840
1841	ALLOC_ERR=-1
1842	ALLOCATE(xinprec(nbindex,nmon), STAT=ALLOC_ERR)
1843	IF (ALLOC_ERR/=0) THEN
1844	WRITE(numout,*) "ERROR IN ALLOCATION of xinprec : ",ALLOC_ERR
1845	STOP
1846	ENDIF
1847	var_name='prm'
1848	CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, &
1849	iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm)
1850	xinprec(:,:)=xchampm(kindex(1:nbindex),:)
1851	!-
1852
1853	ALLOC_ERR=-1
1854	ALLOCATE(xint(nbindex,nmon), STAT=ALLOC_ERR)
1855	IF (ALLOC_ERR/=0) THEN
1856	WRITE(numout,*) "ERROR IN ALLOCATION of xint : ",ALLOC_ERR
1857	STOP
1858	ENDIF
1859	var_name='t2m'
1860	CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, &
1861	iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm)
1862	xint(:,:)=xchampm(kindex(1:nbindex),:)
1863	!-
1864
1865	ALLOC_ERR=-1
1866	ALLOCATE(xinq(nbindex,nmon), STAT=ALLOC_ERR)
1867	IF (ALLOC_ERR/=0) THEN
1868	WRITE(numout,*) "ERROR IN ALLOCATION of xinq : ",ALLOC_ERR
1869	STOP
1870	ENDIF
1871	var_name='r2m'
1872	CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, &
1873	iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm)
1874	xinq(:,:)=xchampm(kindex(1:nbindex),:)
1875	!-
1876
1877	ALLOC_ERR=-1
1878	ALLOCATE(xinwind(nbindex,nmon), STAT=ALLOC_ERR)
1879	IF (ALLOC_ERR/=0) THEN
1880	WRITE(numout,*) "ERROR IN ALLOCATION of xinwind : ",ALLOC_ERR
1881	STOP
1882	ENDIF
1883	var_name='uv10m'
1884	CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, &
1885	iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm)
1886	xinwind(:,:)=xchampm(kindex(1:nbindex),:)
1887	!-
1888
1889	ALLOC_ERR=-1
1890	ALLOCATE(xincld(nbindex,nmon), STAT=ALLOC_ERR)
1891	IF (ALLOC_ERR/=0) THEN
1892	WRITE(numout,*) "ERROR IN ALLOCATION of xincld : ",ALLOC_ERR
1893	STOP
1894	ENDIF
1895	var_name='tc'
1896	CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, &
1897	iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm)
1898	xincld(:,:)=xchampm(kindex(1:nbindex),:)
1899	!-
1900
1901	ALLOC_ERR=-1
1902	ALLOCATE(xintrng(nbindex,nmon), STAT=ALLOC_ERR)
1903	IF (ALLOC_ERR/=0) THEN
1904	WRITE(numout,*) "ERROR IN ALLOCATION of xintrng : ",ALLOC_ERR
1905	STOP
1906	ENDIF
1907	var_name='t2a'
1908	CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, &
1909	iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm)
1910	xintrng(:,:)=xchampm(kindex(1:nbindex),:)
1911	!-
1912	! 10. clean up
1913	!-
1914	IF (is_root_prc) THEN
1915	DEALLOCATE (ncorr)
1916	DEALLOCATE (icorr)
1917	DEALLOCATE (jcorr)
1918	ENDIF
1919
1920	!-
1921	! 12. Allocate space for daily mean values
1922	!-
1923
1924	ALLOC_ERR=-1
1925	ALLOCATE(iwet(nbindex), STAT=ALLOC_ERR)
1926	IF (ALLOC_ERR/=0) THEN
1927	WRITE(numout,*) "ERROR IN ALLOCATION of iwet : ",ALLOC_ERR
1928	STOP
1929	ENDIF
1930	!-
1931
1932	ALLOC_ERR=-1
1933	ALLOCATE(psurfm0(nbindex), STAT=ALLOC_ERR)
1934	IF (ALLOC_ERR/=0) THEN
1935	WRITE(numout,*) "ERROR IN ALLOCATION of psurfm0 : ",ALLOC_ERR
1936	STOP
1937	ENDIF
1938
1939	ALLOC_ERR=-1
1940	ALLOCATE(cloudm0(nbindex), STAT=ALLOC_ERR)
1941	IF (ALLOC_ERR/=0) THEN
1942	WRITE(numout,*) "ERROR IN ALLOCATION of cloudm0 : ",ALLOC_ERR
1943	STOP
1944	ENDIF
1945
1946	ALLOC_ERR=-1
1947	ALLOCATE(tmaxm0(nbindex), STAT=ALLOC_ERR)
1948	IF (ALLOC_ERR/=0) THEN
1949	WRITE(numout,*) "ERROR IN ALLOCATION of tmaxm0 : ",ALLOC_ERR
1950	STOP
1951	ENDIF
1952
1953	ALLOC_ERR=-1
1954	ALLOCATE(tminm0(nbindex), STAT=ALLOC_ERR)
1955	IF (ALLOC_ERR/=0) THEN
1956	WRITE(numout,*) "ERROR IN ALLOCATION of tminm0 : ",ALLOC_ERR
1957	STOP
1958	ENDIF
1959
1960	ALLOC_ERR=-1
1961	ALLOCATE(qdm0(nbindex), STAT=ALLOC_ERR)
1962	IF (ALLOC_ERR/=0) THEN
1963	WRITE(numout,*) "ERROR IN ALLOCATION of qdm0 : ",ALLOC_ERR
1964	STOP
1965	ENDIF
1966
1967	ALLOC_ERR=-1
1968	ALLOCATE(udm0(nbindex), STAT=ALLOC_ERR)
1969	IF (ALLOC_ERR/=0) THEN
1970	WRITE(numout,*) "ERROR IN ALLOCATION of udm0 : ",ALLOC_ERR
1971	STOP
1972	ENDIF
1973
1974	ALLOC_ERR=-1
1975	ALLOCATE(precipm0(nbindex), STAT=ALLOC_ERR)
1976	IF (ALLOC_ERR/=0) THEN
1977	WRITE(numout,*) "ERROR IN ALLOCATION of precipm0 : ",ALLOC_ERR
1978	STOP
1979	ENDIF
1980	!-
1981
1982	ALLOC_ERR=-1
1983	ALLOCATE(psurfm1(nbindex), STAT=ALLOC_ERR)
1984	IF (ALLOC_ERR/=0) THEN
1985	WRITE(numout,*) "ERROR IN ALLOCATION of psurfm1 : ",ALLOC_ERR
1986	STOP
1987	ENDIF
1988
1989	ALLOC_ERR=-1
1990	ALLOCATE(cloudm1(nbindex), STAT=ALLOC_ERR)
1991	IF (ALLOC_ERR/=0) THEN
1992	WRITE(numout,*) "ERROR IN ALLOCATION of cloudm1 : ",ALLOC_ERR
1993	STOP
1994	ENDIF
1995
1996	ALLOC_ERR=-1
1997	ALLOCATE(tmaxm1(nbindex), STAT=ALLOC_ERR)
1998	IF (ALLOC_ERR/=0) THEN
1999	WRITE(numout,*) "ERROR IN ALLOCATION of tmaxm1 : ",ALLOC_ERR
2000	STOP
2001	ENDIF
2002
2003	ALLOC_ERR=-1
2004	ALLOCATE(tminm1(nbindex), STAT=ALLOC_ERR)
2005	IF (ALLOC_ERR/=0) THEN
2006	WRITE(numout,*) "ERROR IN ALLOCATION of tminm1 : ",ALLOC_ERR
2007	STOP
2008	ENDIF
2009
2010	ALLOC_ERR=-1
2011	ALLOCATE(qdm1(nbindex), STAT=ALLOC_ERR)
2012	IF (ALLOC_ERR/=0) THEN
2013	WRITE(numout,*) "ERROR IN ALLOCATION of qdm1 : ",ALLOC_ERR
2014	STOP
2015	ENDIF
2016
2017	ALLOC_ERR=-1
2018	ALLOCATE(udm1(nbindex), STAT=ALLOC_ERR)
2019	IF (ALLOC_ERR/=0) THEN
2020	WRITE(numout,*) "ERROR IN ALLOCATION of udm1 : ",ALLOC_ERR
2021	STOP
2022	ENDIF
2023
2024	ALLOC_ERR=-1
2025	ALLOCATE(precipm1(nbindex), STAT=ALLOC_ERR)
2026	IF (ALLOC_ERR/=0) THEN
2027	WRITE(numout,*) "ERROR IN ALLOCATION of precipm1 : ",ALLOC_ERR
2028	STOP
2029	ENDIF
2030	END SUBROUTINE weathgen_read_file
2031
2032	SUBROUTINE weathgen_begin ( &
2033	& dt_force,itau, date0, &
2034	& rest_id,iim,jjm, &
2035	& lon,lat,tair,pb,qair,kindex,nbindex)
2036	!---------------------------------------------------------------------
2037	IMPLICIT NONE
2038
2039	!-
2040	REAL,INTENT(IN) :: dt_force, date0
2041	INTEGER,INTENT(IN) :: itau
2042	INTEGER,INTENT(IN) :: rest_id
2043	INTEGER,INTENT(IN) :: iim, jjm
2044	REAL,DIMENSION(iim,jjm),INTENT(IN) :: lon,lat
2045	INTEGER,DIMENSION(iim*jjm),INTENT(OUT) :: kindex
2046	INTEGER,INTENT(OUT) :: nbindex
2047	!-
2048	REAL,DIMENSION(iim,jjm),INTENT(OUT) :: tair,pb,qair
2049	INTEGER :: i, j, ij
2050	REAL :: val_exp
2051	REAL,DIMENSION(iim*jjm) :: xchamp
2052	REAL,DIMENSION(iim_g*jjm_g) :: xchamp_g
2053	CHARACTER(LEN=30) :: var_name
2054	REAL,DIMENSION(1) :: jullasttab
2055	REAL,DIMENSION(seedsize_max) :: seed_in_file
2056	INTEGER,DIMENSION(:), ALLOCATABLE :: seed_in_proc
2057	INTEGER :: seedsize, iret
2058	INTEGER :: nlonid, nlatid, nlonid1, nlatid1, tdimid1, varid
2059	INTEGER :: ndim, dims(4)
2060	CHARACTER(LEN=30) :: assoc
2061	CHARACTER(LEN=70) :: str70
2062	CHARACTER(LEN=80) :: stamp
2063	INTEGER :: yy_b, mm_b, dd_b, hh, mm
2064	REAL :: ss
2065	CHARACTER(LEN=10) :: today, att
2066	INTEGER :: nlandid1, nlandid, nlevid, nlevid1
2067	REAL :: lev_max, lev_min
2068	REAL :: height_lev1
2069	INTEGER :: imois
2070	REAL :: xx, td
2071
2072	kindex(:)=0
2073
2074	#ifdef CPP_PARA
2075	nbindex=nbp_loc
2076	CALL scatter(index_g,kindex)
2077	kindex(1:nbindex)=kindex(1:nbindex)-(jj_begin-1)*iim_g
2078	#else
2079	! Copy saved land points index
2080	nbindex = nbindex_w
2081	kindex(1:nbindex_w) = kindex_w(:)
2082	#endif
2083	!-
2084	! 13. Prescribed or statistical values?
2085	!-
2086	!Config Key = IPPREC
2087	!Config Desc = Use prescribed values
2088	!Config If = ALLOW_WEATHERGEN
2089	!Config Def = 0
2090	!Config Help = If this is set to 1, the weather generator
2091	!Config uses the monthly mean values for daily means.
2092	!Config If it is set to 0, the weather generator
2093	!Config uses statistical relationships to derive daily
2094	!Config values from monthly means.
2095	ipprec = 0
2096	CALL getin_p ('IPPREC',ipprec)
2097	WRITE(numout,*) 'IPPREC: ',ipprec
2098	!-
2099	! 14. Do we want exact monthly precipitations even with ipprec=0 ?
2100	!-
2101	!Config Key = WEATHGEN_PRECIP_EXACT
2102	!Config Desc = Exact monthly precipitation
2103	!Config If = ALLOW_WEATHERGEN
2104	!Config Def = n
2105	!Config Help = If this is set to y, the weather generator
2106	!Config will generate pseudo-random precipitations
2107	!Config whose monthly mean is exactly the prescribed one.
2108	!Config In this case, the daily precipitation (for rainy
2109	!Config days) is constant (that is, some days have 0 precip,
2110	!Config the other days have precip=Precip_month/n_precip,
2111	!Config where n_precip is the prescribed number of rainy days
2112	!Config per month).
2113	precip_exact = .FALSE.
2114	CALL getin_p ('WEATHGEN_PRECIP_EXACT',precip_exact)
2115	WRITE(numout,*) 'PRECIP_EXACT: ',precip_exact
2116	!-
2117	IF (precip_exact) THEN
2118	!---
2119	!-- preparer un tableau utilise pour determiner s'il pleut ou pas
2120	!-- (en fct. du nombre de jours de pluie par mois).
2121	!---
2122	IF (is_root_prc) THEN
2123	DO imois=1,12
2124	CALL permute (ndaypm(imois),jour_precip(:,imois))
2125	ENDDO
2126	ENDIF
2127	CALL bcast(jour_precip)
2128	ENDIF
2129	!-
2130	! Read Orbital Parameters
2131	!-
2132	!Config Key = ECCENTRICITY
2133	!Config Desc = Use prescribed values
2134	!Config If = ALLOW_WEATHERGEN
2135	!Config Def = 0.016724
2136	ecc = 0.016724
2137	CALL getin_p ('ECCENTRICITY',ecc)
2138	WRITE(numout,*) 'ECCENTRICITY: ',ecc
2139	!
2140	!Config Key = PERIHELIE
2141	!Config Desc = Use prescribed values
2142	!Config If = ALLOW_WEATHERGEN
2143	!Config Def = 102.04
2144	perh = 102.04
2145	CALL getin_p ('PERIHELIE',perh)
2146	WRITE(numout,*) 'PERIHELIE: ',perh
2147	!
2148	!Config Key = OBLIQUITY
2149	!Config Desc = Use prescribed values
2150	!Config If = ALLOW_WEATHERGEN
2151	!Config Def = 23.446
2152	xob = 23.446
2153	CALL getin_p ('OBLIQUITY',xob)
2154	WRITE(numout,*) 'OBLIQUITY: ',xob
2155	!-
2156	! 15. Read restart file
2157	!-
2158	CALL ioget_expval (val_exp)
2159	!-
2160	var_name= 'julian'
2161	IF (is_root_prc) THEN
2162	CALL restget (rest_id,var_name,1,1,1,itau,.TRUE.,jullasttab)
2163	IF (jullasttab(1) == val_exp) THEN
2164	jullasttab(1) = itau2date(itau-1, date0, dt_force)
2165	ENDIF
2166	ENDIF
2167	CALL bcast(jullasttab)
2168	julian_last = jullasttab(1)
2169	!-
2170	var_name= 'seed'
2171	IF (is_root_prc) &
2172	CALL restget (rest_id,var_name,seedsize_max, &
2173	& 1,1,itau,.TRUE.,seed_in_file)
2174	CALL bcast(seed_in_file)
2175	IF (ALL(seed_in_file(:) == val_exp)) THEN
2176	!---
2177	!-- there is no need to reinitialize the random number generator as
2178	!-- this does not seem to be a restart
2179	!---
2180	CONTINUE
2181	ELSE
2182	!---
2183	!-- reinitialize the random number generator
2184	!---
2185	IF (is_root_prc) &
2186	CALL RANDOM_SEED( SIZE = seedsize )
2187	CALL bcast(seedsize)
2188	IF (seedsize > seedsize_max) THEN
2189	STOP 'weathgen_begin: increase seedsize_max'
2190	ENDIF
2191
2192	ALLOC_ERR=-1
2193	ALLOCATE(seed_in_proc(seedsize), STAT=ALLOC_ERR)
2194	IF (ALLOC_ERR/=0) THEN
2195	WRITE(numout,*) "ERROR IN ALLOCATION of seed_in_proc : ",ALLOC_ERR
2196	STOP
2197	ENDIF
2198	seed_in_proc(1:seedsize) = NINT( seed_in_file(1:seedsize) )
2199	CALL RANDOM_SEED (PUT = seed_in_proc)
2200	DEALLOCATE( seed_in_proc )
2201	ENDIF
2202	!-
2203	var_name= 'iwet'
2204	IF (is_root_prc) THEN
2205	CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
2206	IF (ALL(xchamp_g(:) == val_exp)) THEN
2207	xchamp_g(:) = 0.
2208	ENDIF
2209	ENDIF
2210	CALL scatter2D(xchamp_g,xchamp)
2211	iwet(:) = NINT(xchamp(kindex(1:nbindex)))
2212	!-
2213	var_name= 'psurfm0'
2214	IF (is_root_prc) THEN
2215	CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
2216	IF (ALL(xchamp_g(:) == val_exp)) THEN
2217	xchamp_g(:) = 100000.
2218	ENDIF
2219	ENDIF
2220	CALL scatter2D(xchamp_g,xchamp)
2221	psurfm0(:) = xchamp(kindex(1:nbindex))
2222	!-
2223	var_name= 'cloudm0'
2224	IF (is_root_prc) THEN
2225	CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
2226	IF (ALL(xchamp_g(:) == val_exp)) THEN
2227	xchamp_g(:) = .5
2228	ENDIF
2229	ENDIF
2230	CALL scatter2D(xchamp_g,xchamp)
2231	cloudm0(:) = xchamp(kindex(1:nbindex))
2232	!-
2233	var_name= 'tmaxm0'
2234	IF (is_root_prc) THEN
2235	CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
2236	IF (ALL(xchamp_g(:) == val_exp)) THEN
2237	xchamp_g(:) = 285.
2238	ENDIF
2239	ENDIF
2240	CALL scatter2D(xchamp_g,xchamp)
2241	tmaxm0(:) = xchamp(kindex(1:nbindex))
2242	!-
2243	var_name= 'tminm0'
2244	IF (is_root_prc) THEN
2245	CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
2246	IF (ALL(xchamp_g(:) == val_exp)) THEN
2247	xchamp_g(:) = 275.
2248	ENDIF
2249	ENDIF
2250	CALL scatter2D(xchamp_g,xchamp)
2251	tminm0(:) = xchamp(kindex(1:nbindex))
2252	!-
2253	var_name= 'qdm0'
2254	IF (is_root_prc) THEN
2255	CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
2256	IF (ALL(xchamp_g(:) == val_exp)) THEN
2257	xchamp_g(:) = 1.E-03
2258	ENDIF
2259	ENDIF
2260	CALL scatter2D(xchamp_g,xchamp)
2261	qdm0(:) = xchamp(kindex(1:nbindex))
2262	!-
2263	var_name= 'udm0'
2264	IF (is_root_prc) THEN
2265	CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
2266	IF (ALL(xchamp_g(:) == val_exp)) THEN
2267	xchamp_g(:) = 2.
2268	ENDIF
2269	ENDIF
2270	CALL scatter2D(xchamp_g,xchamp)
2271	udm0(:) = xchamp(kindex(1:nbindex))
2272	!-
2273	var_name= 'precipm0'
2274	IF (is_root_prc) THEN
2275	CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
2276	IF (ALL(xchamp_g(:) == val_exp)) THEN
2277	xchamp_g(:) = 1.
2278	ENDIF
2279	ENDIF
2280	CALL scatter2D(xchamp_g,xchamp)
2281	precipm0(:) = xchamp(kindex(1:nbindex))
2282	!-
2283	var_name= 'psurfm1'
2284	IF (is_root_prc) THEN
2285	CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
2286	IF (ALL(xchamp_g(:) == val_exp)) THEN
2287	xchamp_g(:) = 100000.
2288	ENDIF
2289	ENDIF
2290	CALL scatter2D(xchamp_g,xchamp)
2291	psurfm1(:) = xchamp(kindex(1:nbindex))
2292	!-
2293	var_name= 'cloudm1'
2294	IF (is_root_prc) THEN
2295	CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
2296	IF (ALL(xchamp_g(:) == val_exp)) THEN
2297	xchamp_g(:) = .5
2298	ENDIF
2299	ENDIF
2300	CALL scatter2D(xchamp_g,xchamp)
2301	cloudm1(:) = xchamp(kindex(1:nbindex))
2302	!-
2303	var_name= 'tmaxm1'
2304	IF (is_root_prc) THEN
2305	CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
2306	IF (ALL(xchamp_g(:) == val_exp)) THEN
2307	xchamp_g(:) = 285.
2308	ENDIF
2309	ENDIF
2310	CALL scatter2D(xchamp_g,xchamp)
2311	tmaxm1(:) = xchamp(kindex(1:nbindex))
2312	!-
2313	var_name= 'tminm1'
2314	IF (is_root_prc) THEN
2315	CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
2316	IF (ALL(xchamp_g(:) == val_exp)) THEN
2317	xchamp_g(:) = 275.
2318	ENDIF
2319	ENDIF
2320	CALL scatter2D(xchamp_g,xchamp)
2321	tminm1(:) = xchamp(kindex(1:nbindex))
2322	!-
2323	var_name= 'qdm1'
2324	IF (is_root_prc) THEN
2325	CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
2326	IF (ALL(xchamp_g(:) == val_exp)) THEN
2327	xchamp_g(:) = 1.E-03
2328	ENDIF
2329	ENDIF
2330	CALL scatter2D(xchamp_g,xchamp)
2331	qdm1(:) = xchamp(kindex(1:nbindex))
2332	!-
2333	var_name= 'udm1'
2334	IF (is_root_prc) THEN
2335	CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
2336	IF (ALL(xchamp_g(:) == val_exp)) THEN
2337	xchamp_g(:) = 2.
2338	ENDIF
2339	ENDIF
2340	CALL scatter2D(xchamp_g,xchamp)
2341	udm1(:) = xchamp(kindex(1:nbindex))
2342	!-
2343	var_name= 'precipm1'
2344	IF (is_root_prc) THEN
2345	CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
2346	IF (ALL(xchamp_g(:) == val_exp)) THEN
2347	xchamp_g(:) = 1.
2348	ENDIF
2349	ENDIF
2350	CALL scatter2D(xchamp_g,xchamp)
2351	precipm1(:) = xchamp(kindex(1:nbindex))
2352	!-
2353	! 16. We still need instantaneous tair, qair, and the surface pressure
2354	! We take daily mean values read from the restart file
2355	!-
2356	!!$ tair(:,:)=280.
2357	!!$ qair(:,:)=1.E-03
2358	!!$ pb(:,:)=101325
2359	tair(:,:)=val_exp
2360	qair(:,:)=val_exp
2361	pb(:,:)=val_exp
2362	xx = 9.81/287./0.0065
2363	DO ij=1,nbindex
2364	j = ((kindex(ij)-1)/iim) + 1
2365	i = kindex(ij) - (j-1)*iim
2366
2367	lat_land(ij) = lat(i,j)
2368
2369	td = (tmaxm0(ij)+tminm0(ij))/2.
2370	tair(i,j) = td
2371	qair(i,j) = qdm1(ij)
2372	pb(i,j) = 101325.(td/(td+0.0065xintopo(ij)))**xx
2373	ENDDO
2374	!-
2375	! 17. We can write a forcing file for Orchidee
2376	! from this weather Generator run.
2377	!-
2378	!Config Key = DUMP_WEATHER
2379	!Config Desc = Write weather from generator into a forcing file
2380	!Config Def = n
2381	!Config Help = This flag makes the weather generator dump its
2382	! generated weather into a forcing file which can
2383	! then be used to get the same forcing on different
2384	! machines. This only works correctly if there is
2385	! a restart file (otherwise the forcing at the first
2386	! time step is slightly wrong).
2387	dump_weather = .FALSE.
2388	CALL getin_p ('DUMP_WEATHER',dump_weather)
2389	!-
2390	IF (dump_weather .AND. is_root_prc) THEN
2391	!---
2392	!-- Initialize the file
2393	!---
2394	!Config Key = DUMP_WEATHER_FILE
2395	!Config Desc = Name of the file that contains
2396	! the weather from generator
2397	!Config Def = 'weather_dump.nc'
2398	!Config If = DUMP_WEATHER
2399	!Config Help =
2400	dump_weather_file = 'weather_dump.nc'
2401	CALL getin ('DUMP_WEATHER_FILE',dump_weather_file)
2402	!---
2403	!Config Key = DUMP_WEATHER_GATHERED
2404	!Config Desc = Dump weather data on gathered grid
2405	!Config Def = y
2406	!Config If = DUMP_WEATHER
2407	!Config Help = If 'y', the weather data are gathered
2408	! for all land points.
2409	gathered = .TRUE.
2410	CALL getin ('DUMP_WEATHER_GATHERED',gathered)
2411	!---
2412	iret = NF90_CREATE (TRIM(dump_weather_file),NF90_CLOBBER,dump_id)
2413	!---
2414	!-- Dimensions
2415	!---
2416	iret = NF90_DEF_DIM (dump_id,'x',iim_g,nlonid1)
2417	iret = NF90_DEF_DIM (dump_id,'y',jjm_g,nlatid1)
2418	iret = NF90_DEF_DIM (dump_id,'z', 1,nlevid1)
2419	!---
2420	IF (gathered) THEN
2421	iret = NF90_DEF_DIM (dump_id,'land',nbp_glo,nlandid1)
2422	ENDIF
2423	iret = NF90_DEF_DIM (dump_id,'tstep',NF90_UNLIMITED,tdimid1)
2424	!---
2425	!-- Coordinate variables
2426	!---
2427	dims(1:2) = (/ nlonid1, nlatid1 /)
2428	!---
2429	iret = NF90_DEF_VAR (dump_id,'nav_lon',n_rtp,dims(1:2),nlonid)
2430	iret = NF90_PUT_ATT (dump_id,nlonid,'units',"degrees_east")
2431	iret = NF90_PUT_ATT (dump_id,nlonid,'valid_min',MINVAL(lon_g))
2432	iret = NF90_PUT_ATT (dump_id,nlonid,'valid_max',MAXVAL(lon_g))
2433	iret = NF90_PUT_ATT (dump_id,nlonid,'long_name',"Longitude")
2434	!---
2435	iret = NF90_DEF_VAR (dump_id,'nav_lat',n_rtp,dims(1:2),nlatid)
2436	iret = NF90_PUT_ATT (dump_id,nlatid,'units',"degrees_north")
2437	iret = NF90_PUT_ATT (dump_id,nlatid,'valid_min',MINVAL(lat_g))
2438	iret = NF90_PUT_ATT (dump_id,nlatid,'valid_max',MAXVAL(lat_g))
2439	iret = NF90_PUT_ATT (dump_id,nlatid,'long_name',"Latitude")
2440	!---
2441	height_lev1 = 10.
2442	CALL getin ('HEIGHT_LEV1',height_lev1)
2443	lev_min = height_lev1
2444	lev_max = height_lev1
2445	!---
2446	iret = NF90_DEF_VAR (dump_id,'level',n_rtp,(/ nlevid1 /),nlevid)
2447	iret = NF90_PUT_ATT (dump_id,nlevid,'units',"m")
2448	iret = NF90_PUT_ATT (dump_id,nlevid,'valid_min',lev_min)
2449	iret = NF90_PUT_ATT (dump_id,nlevid,'valid_max',lev_max)
2450	iret = NF90_PUT_ATT (dump_id,nlevid,'long_name',"Vertical levels")
2451	!---
2452	IF (gathered) THEN
2453	iret = NF90_DEF_VAR (dump_id,'land',NF90_INT,(/ nlandid1 /),nlandid)
2454	iret = NF90_PUT_ATT (dump_id,nlandid,'compress',"y x")
2455	ENDIF
2456	!---
2457	!-- Store the time axes.
2458	!---
2459	iret = NF90_DEF_VAR (dump_id,'time',n_rtp,tdimid1,time_id)
2460
2461	yy_b=0
2462	mm_b=1
2463	dd_b=1
2464	hh=00
2465	mm=00
2466	ss=0.
2467
2468	WRITE (str70,7000) yy_b, mm_b, dd_b, hh, mm, INT(ss)
2469	iret = NF90_PUT_ATT (dump_id,time_id,'units',TRIM(str70))
2470	iret = NF90_PUT_ATT (dump_id,time_id,'calendar',TRIM(calendar_str))
2471	iret = NF90_PUT_ATT (dump_id,time_id,'title','Time')
2472	iret = NF90_PUT_ATT (dump_id,time_id,'long_name','Time axis')
2473	WRITE(str70,7001) yy_b, cal(mm_b), dd_b, hh, mm, INT(ss)
2474	iret = NF90_PUT_ATT (dump_id,time_id,'time_origin',TRIM(str70))
2475	!---
2476	!-- Time steps
2477	!---
2478	iret = NF90_DEF_VAR (dump_id,'timestp',NF90_INT,tdimid1,timestp_id)
2479	WRITE(str70,7002) yy_b, mm_b, dd_b, hh, mm, INT(ss)
2480	iret = NF90_PUT_ATT (dump_id,timestp_id,'units',TRIM(str70))
2481	iret = NF90_PUT_ATT (dump_id,timestp_id,'title','Time steps')
2482	iret = NF90_PUT_ATT (dump_id,timestp_id,'tstep_sec',dt_force)
2483	iret = NF90_PUT_ATT &
2484	& (dump_id,timestp_id,'long_name','Time step axis')
2485	WRITE(str70,7001) yy_b, cal(mm_b), dd_b, hh, mm, INT(ss)
2486	iret = NF90_PUT_ATT (dump_id,timestp_id,'time_origin',TRIM(str70))
2487	!---
2488	7000 FORMAT('seconds since ',I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2)
2489	7001 FORMAT(' ',I4.4,'-',A3,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2)
2490	7002 FORMAT('timesteps since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2)
2491	!---
2492	!-- The variables in the file are declared and their attributes
2493	!-- written.
2494	!---
2495	IF (gathered) THEN
2496	ndim = 2
2497	dims(1:2) = (/ nlandid1, tdimid1 /)
2498	assoc = 'time (nav_lat nav_lon)'
2499	ELSE
2500	ndim = 3
2501	dims(1:3) = (/ nlonid1, nlatid1, tdimid1 /)
2502	assoc = 'time nav_lat nav_lon'
2503	ENDIF
2504	!---
2505	iret = NF90_DEF_VAR (dump_id,'SWdown',n_rtp,dims(1:ndim),varid)
2506	iret = NF90_PUT_ATT (dump_id,varid,'axis','TYX')
2507	iret = NF90_PUT_ATT (dump_id,varid,'units','W/m^2')
2508	iret = NF90_PUT_ATT (dump_id,varid,'long_name', &
2509	& 'Surface incident shortwave radiation')
2510	iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc))
2511	iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef_sechiba)
2512	soldownid = varid
2513	!---
2514	iret = NF90_DEF_VAR (dump_id,'Rainf',n_rtp,dims(1:ndim),varid)
2515	iret = NF90_PUT_ATT (dump_id,varid,'axis','TYX')
2516	iret = NF90_PUT_ATT (dump_id,varid,'units','kg/m^2s')
2517	iret = NF90_PUT_ATT (dump_id,varid,'long_name', &
2518	& 'Rainfall rate')
2519	iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc))
2520	iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef_sechiba)
2521	rainfid = varid
2522	!---
2523	iret = NF90_DEF_VAR (dump_id,'Snowf',n_rtp,dims(1:ndim),varid)
2524	iret = NF90_PUT_ATT (dump_id,varid,'axis','TYX')
2525	iret = NF90_PUT_ATT (dump_id,varid,'units','kg/m^2s')
2526	iret = NF90_PUT_ATT (dump_id,varid,'long_name', &
2527	& 'Snowfall rate')
2528	iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc))
2529	iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef_sechiba)
2530	snowfid = varid
2531	!---
2532	iret = NF90_DEF_VAR (dump_id,'LWdown',n_rtp,dims(1:ndim),varid)
2533	iret = NF90_PUT_ATT (dump_id,varid,'axis','TYX')
2534	iret = NF90_PUT_ATT (dump_id,varid,'units','W/m^2')
2535	iret = NF90_PUT_ATT (dump_id,varid,'long_name', &
2536	& 'Surface incident longwave radiation')
2537	iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc))
2538	iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef_sechiba)
2539	lwradid = varid
2540	!---
2541	iret = NF90_DEF_VAR (dump_id,'PSurf',n_rtp,dims(1:ndim),varid)
2542	iret = NF90_PUT_ATT (dump_id,varid,'axis','TYX')
2543	iret = NF90_PUT_ATT (dump_id,varid,'units','Pa')
2544	iret = NF90_PUT_ATT (dump_id,varid,'long_name', &
2545	& 'Surface pressure')
2546	iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc))
2547	iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef_sechiba)
2548	psolid = varid
2549	!---
2550	!-- 3D Variables to be written
2551	!---
2552	IF (gathered) THEN
2553	ndim = 3
2554	dims(1:3) = (/ nlandid1, nlevid1, tdimid1 /)
2555	assoc = 'time level (nav_lat nav_lon)'
2556	ELSE
2557	ndim = 4
2558	dims(1:4) = (/ nlonid1, nlatid1, nlevid1, tdimid1 /)
2559	assoc = 'time level nav_lat nav_lon'
2560	ENDIF
2561	!---
2562	iret = NF90_DEF_VAR (dump_id,'Tair',n_rtp,dims(1:ndim),varid)
2563	iret = NF90_PUT_ATT (dump_id,varid,'axis','TZYX')
2564	iret = NF90_PUT_ATT (dump_id,varid,'units','K')
2565	iret = NF90_PUT_ATT (dump_id,varid,'long_name', &
2566	& 'Near surface air temperature')
2567	iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc))
2568	iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef_sechiba)
2569	tairid = varid
2570	!---
2571	iret = NF90_DEF_VAR (dump_id,'Qair',n_rtp,dims(1:ndim),varid)
2572	iret = NF90_PUT_ATT (dump_id,varid,'axis','TZYX')
2573	iret = NF90_PUT_ATT (dump_id,varid,'units','kg/kg')
2574	iret = NF90_PUT_ATT (dump_id,varid,'long_name', &
2575	& 'Near surface specific humidity')
2576	iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc))
2577	iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef_sechiba)
2578	qairid = varid
2579	!---
2580	iret = NF90_DEF_VAR (dump_id,'Wind_N',n_rtp,dims(1:ndim),varid)
2581	iret = NF90_PUT_ATT (dump_id,varid,'axis','TZYX')
2582	iret = NF90_PUT_ATT (dump_id,varid,'units','m/s')
2583	iret = NF90_PUT_ATT (dump_id,varid,'long_name', &
2584	& 'Near surface northward wind component')
2585	iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc))
2586	iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef_sechiba)
2587	uid = varid
2588	!---
2589	iret = NF90_DEF_VAR (dump_id,'Wind_E',n_rtp,dims(1:ndim),varid)
2590	iret = NF90_PUT_ATT (dump_id,varid,'axis','TZYX')
2591	iret = NF90_PUT_ATT (dump_id,varid,'units','m/s')
2592	iret = NF90_PUT_ATT (dump_id,varid,'long_name', &
2593	& 'Near surface eastward wind component')
2594	iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc))
2595	iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef_sechiba)
2596	vid = varid
2597	!---
2598	!-- Global attributes
2599	!---
2600	CALL DATE_AND_TIME (today, att)
2601	stamp = "WG, date: "//TRIM(today)//" at "//TRIM(att)
2602	iret = NF90_PUT_ATT (dump_id,NF90_GLOBAL,'Conventions',"GDT 1.2")
2603	iret = NF90_PUT_ATT (dump_id,NF90_GLOBAL,'file_name', &
2604	& TRIM(dump_weather_file))
2605	iret = NF90_PUT_ATT (dump_id,NF90_GLOBAL,'production',TRIM(stamp))
2606	!---
2607	!-- Finish the definition phase
2608	!---
2609	iret = NF90_ENDDEF (dump_id)
2610	!---
2611	!-- Write coordinates
2612	!---
2613	iret = NF90_PUT_VAR (dump_id,nlonid,lon_g)
2614	IF (iret /= NF90_NOERR) THEN
2615	WRITE(numout,*) iret
2616	CALL ipslerr (3,'weathgen_begin', &
2617	& 'Could not put variable nav_lon in the file : ', &
2618	& TRIM(dump_weather_file),'(Solution ?)')
2619	ENDIF
2620	iret = NF90_PUT_VAR (dump_id,nlatid,lat_g)
2621	IF (iret /= NF90_NOERR) THEN
2622	WRITE(numout,*) iret
2623	CALL ipslerr (3,'weathgen_begin', &
2624	& 'Could not put variable nav_lat in the file : ', &
2625	& TRIM(dump_weather_file),'(Solution ?)')
2626	ENDIF
2627	iret = NF90_PUT_VAR (dump_id,nlevid,height_lev1)
2628	IF (iret /= NF90_NOERR) THEN
2629	WRITE(numout,*) iret
2630	CALL ipslerr (3,'weathgen_begin', &
2631	& 'Could not put variable level in the file : ', &
2632	& TRIM(dump_weather_file),'(Solution ?)')
2633	ENDIF
2634	!---
2635	IF (gathered) THEN
2636	iret = NF90_PUT_VAR (dump_id,nlandid,index_g(1:nbp_glo))
2637	IF (iret /= NF90_NOERR) THEN
2638	WRITE(numout,*) iret
2639	CALL ipslerr (3,'weathgen_begin', &
2640	& 'Could not put variable land in the file : ', &
2641	& TRIM(dump_weather_file),'(Solution ?)')
2642	ENDIF
2643	ENDIF
2644	!---
2645	ENDIF ! dump_weather
2646	!-----------------------------
2647	END SUBROUTINE weathgen_begin
2648	!-
2649	!===
2650	!-
2651	SUBROUTINE weathgen_get &
2652	& (itau, date0, dt_force, nbindex, nband, lat, &
2653	& swdown, raina, snowa, tair, u, v, qair, psurf, lwdown)
2654	!---------------------------------------------------------------------
2655	IMPLICIT NONE
2656	! number of time step
2657	INTEGER,INTENT(IN) :: itau
2658	! date when itau was 0
2659	REAL,INTENT(IN) :: date0
2660	! time step (s)
2661	REAL,INTENT(IN) :: dt_force
2662	! number of land points
2663	INTEGER,INTENT(IN) :: nbindex
2664	! number of visible bands
2665	INTEGER,INTENT(IN) :: nband
2666	! latitude (deg)
2667	REAL,DIMENSION(nbindex),INTENT(IN) :: lat
2668	!-
2669	REAL,DIMENSION(nbindex,nband),INTENT(OUT) :: swdown
2670	REAL,DIMENSION(nbindex),INTENT(OUT) :: raina, snowa
2671	REAL,DIMENSION(nbindex),INTENT(OUT) :: tair
2672	REAL,DIMENSION(nbindex),INTENT(OUT) :: u,v
2673	REAL,DIMENSION(nbindex),INTENT(OUT) :: qair
2674	REAL,DIMENSION(nbindex),INTENT(OUT) :: psurf
2675	REAL,DIMENSION(nbindex),INTENT(OUT) :: lwdown
2676	!-
2677	REAL,DIMENSION(nbindex) :: cloud, tmax, tmin, precipd, qd, ud
2678	REAL,DIMENSION(nbindex) :: rh
2679	REAL,DIMENSION(nbindex,nband) :: solai, solad
2680	REAL :: julian, jur
2681	REAL :: x
2682	INTEGER :: yy, mm, dd
2683	REAL :: ss, plens, time
2684	!---------------------------------------------------------------------
2685	!-
2686	! 1. get a reduced julian day
2687	!-
2688	julian = itau2date(itau-1, date0, dt_force)
2689	!SZ, test: solar noon at 12 o'clock!
2690	! julian = itau2date(itau, date0, dt_force)
2691	CALL ju2ymds (julian, yy, mm, dd, ss)
2692	CALL ymds2ju (yy,1,1,0.0, jur)
2693	julian = julian-jur
2694	CALL ju2ymds (julian, yy, mm, dd, ss)
2695	!-
2696	! 2. daily values
2697	!-
2698	IF (INT(julian_last) /= INT(julian)) THEN
2699	!-- today's values become yesterday's values
2700	cloudm1(:) = cloudm0(:)
2701	tmaxm1(:) = tmaxm0(:)
2702	tminm1(:) = tminm0(:)
2703	precipm1(:) = precipm0(:)
2704	qdm1(:) = qdm0(:)
2705	udm1(:) = udm0(:)
2706	psurfm1(:) = psurfm0(:)
2707	!-- we have to get new daily values
2708	!!$ WRITE(,) mpi_rank, "weathgen_get : date ",yy, mm, dd, ss
2709	!!$ WRITE(,) mpi_rank, "weathgen_get : grid date ",year, month, day, sec
2710	CALL daily (nbindex, mm, dd, cloudm0, tmaxm0, tminm0, &
2711	& precipm0, qdm0, udm0, psurfm0)
2712	ENDIF
2713	!-
2714	! 3. interpolate daily values
2715	! (otherwise we get ugly temperature jumps at midnight)
2716	!-
2717	x = (julian-INT(julian))
2718	!-
2719	cloud(:) = (1.-x)cloudm1(:)+xcloudm0(:)
2720	tmax(:) = (1.-x)tmaxm1(:)+xtmaxm0(:)
2721	tmin(:) = (1.-x)tminm1(:)+xtminm0(:)
2722	precipd(:) = (1.-x)precipm1(:)+xprecipm0(:)
2723	qd(:) = (1.-x)qdm1(:)+xqdm0(:)
2724	ud(:) = (1.-x)udm1(:)+xudm0(:)
2725	psurf(:) = (1.-x)psurfm1(:)+xpsurfm0(:)
2726	!-
2727	! 4. read instantaneous values
2728	!-
2729	plens = one_day/dt_force
2730	time = (julian-REAL(INT(julian)))*one_day
2731	!-
2732	CALL diurnal (nbindex, nband, time, NINT(julian), plens, 0., one_day, &
2733	& lat, cloud, tmax, tmin, precipd, qd, ud, psurf, &
2734	& lwdown, solad, solai, u, tair, qair, raina, snowa, rh)
2735	!-
2736	raina(:) = raina(:)/dt_force
2737	snowa(:) = snowa(:)/dt_force
2738	!-
2739	swdown(:,:) = solad(:,:)+solai(:,:)
2740	!-
2741	v(:) = 0.
2742	!-
2743	! 5. Store date
2744	!-
2745	julian_last = julian
2746	!--------------------------
2747	END SUBROUTINE weathgen_get
2748	!-
2749	!===
2750	!-
2751	SUBROUTINE weathgen_restwrite (rest_id,itau,iim,jjm,nbindex,kindex)
2752	!---------------------------------------------------------------------
2753	IMPLICIT NONE
2754	!-
2755	INTEGER,INTENT(IN) :: rest_id,itau,iim,jjm,nbindex
2756	INTEGER,DIMENSION(iim*jjm),INTENT(IN) :: kindex
2757	!-
2758	CHARACTER(LEN=30) :: var_name
2759	INTEGER :: i,j,ij
2760	REAL,DIMENSION(1) :: jullasttab
2761	REAL,DIMENSION(seedsize_max) :: seed_in_file
2762	INTEGER,DIMENSION(:),ALLOCATABLE :: seed_in_proc
2763	INTEGER :: seedsize
2764	REAL :: rndnum
2765	REAL,DIMENSION(iim*jjm) :: xchamp
2766	REAL,DIMENSION(iim_g*jjm_g) :: xchamp_g
2767	!---------------------------------------------------------------------
2768	var_name= 'julian'
2769	jullasttab(:) = julian_last
2770	IF (is_root_prc) CALL restput (rest_id, var_name, 1, 1, 1, itau, jullasttab)
2771	!-
2772	IF (is_root_prc) THEN
2773	CALL RANDOM_SEED( SIZE = seedsize )
2774	IF (seedsize > seedsize_max) THEN
2775	STOP 'weathgen_restwrite: increase seedsize_max'
2776	ENDIF
2777	ENDIF
2778	CALL bcast(seedsize)
2779
2780	IF (is_root_prc) THEN
2781	ALLOC_ERR=-1
2782	ALLOCATE(seed_in_proc(seedsize), STAT=ALLOC_ERR)
2783	IF (ALLOC_ERR/=0) THEN
2784	WRITE(numout,*) "ERROR IN ALLOCATION of seed_in_proc : ",ALLOC_ERR
2785	STOP
2786	ENDIF
2787	!-
2788	CALL RANDOM_SEED (GET = seed_in_proc)
2789	!-
2790	seed_in_file(1:seedsize) = REAL(seed_in_proc(1:seedsize))
2791	!-
2792	! fill in the seed up to seedsize_max
2793	! (useful in the case we restart on
2794	! a machine with a longer seed vector:
2795	! we do not want a degenerated seed)
2796	!-
2797	DO i=seedsize+1,seedsize_max
2798	CALL RANDOM_NUMBER( rndnum )
2799	seed_in_file(i) = 100000.*rndnum
2800	ENDDO
2801	ENDIF
2802	CALL bcast (seed_in_file)
2803	!-
2804	IF (is_root_prc) THEN
2805	DEALLOCATE( seed_in_proc )
2806	!-
2807	var_name= 'seed'
2808	CALL restput (rest_id,var_name,seedsize_max,1,1,itau,seed_in_file)
2809	ENDIF
2810	!-
2811
2812	xchamp(:) = val_exp
2813
2814	!!$ DO j=1,jjm
2815	!!$ DO i=1,iim
2816	!!$ ij = (j-1)*iim+i
2817	!!$ xchamp(i,j) = REAL(iwet(ij))
2818	!!$ ENDDO
2819	!!$ ENDDO
2820	DO ij=1,nbindex
2821	xchamp(kindex(ij)) = REAL(iwet(ij))
2822	ENDDO
2823	var_name= 'iwet'
2824	CALL gather2D(xchamp,xchamp_g)
2825	IF (is_root_prc) THEN
2826	CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
2827	ENDIF
2828	!-
2829	DO ij=1,nbindex
2830	xchamp(kindex(ij)) = psurfm0(ij)
2831	ENDDO
2832	var_name= 'psurfm0'
2833	CALL gather2D(xchamp,xchamp_g)
2834	IF (is_root_prc) THEN
2835	CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
2836	ENDIF
2837	!-
2838	DO ij=1,nbindex
2839	xchamp(kindex(ij)) = cloudm0(ij)
2840	ENDDO
2841	var_name= 'cloudm0'
2842	CALL gather2D(xchamp,xchamp_g)
2843	IF (is_root_prc) THEN
2844	CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
2845	ENDIF
2846	!-
2847	DO ij=1,nbindex
2848	xchamp(kindex(ij)) = tmaxm0(ij)
2849	ENDDO
2850	var_name= 'tmaxm0'
2851	CALL gather2D(xchamp,xchamp_g)
2852	IF (is_root_prc) THEN
2853	CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
2854	ENDIF
2855	!-
2856	DO ij=1,nbindex
2857	xchamp(kindex(ij)) = tminm0(ij)
2858	ENDDO
2859	var_name= 'tminm0'
2860	CALL gather2D(xchamp,xchamp_g)
2861	IF (is_root_prc) THEN
2862	CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
2863	ENDIF
2864	!-
2865	DO ij=1,nbindex
2866	xchamp(kindex(ij)) = qdm0(ij)
2867	ENDDO
2868	var_name= 'qdm0'
2869	CALL gather2D(xchamp,xchamp_g)
2870	IF (is_root_prc) THEN
2871	CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
2872	ENDIF
2873	!-
2874	DO ij=1,nbindex
2875	xchamp(kindex(ij)) = udm0(ij)
2876	ENDDO
2877	var_name= 'udm0'
2878	CALL gather2D(xchamp,xchamp_g)
2879	IF (is_root_prc) THEN
2880	CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
2881	ENDIF
2882	!-
2883	DO ij=1,nbindex
2884	xchamp(kindex(ij)) = precipm0(ij)
2885	ENDDO
2886	var_name= 'precipm0'
2887	CALL gather2D(xchamp,xchamp_g)
2888	IF (is_root_prc) THEN
2889	CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
2890	ENDIF
2891	!-
2892	DO ij=1,nbindex
2893	xchamp(kindex(ij)) = psurfm1(ij)
2894	ENDDO
2895	var_name= 'psurfm1'
2896	CALL gather2D(xchamp,xchamp_g)
2897	IF (is_root_prc) THEN
2898	CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
2899	ENDIF
2900	!-
2901	DO ij=1,nbindex
2902	xchamp(kindex(ij)) = cloudm1(ij)
2903	ENDDO
2904	var_name= 'cloudm1'
2905	CALL gather2D(xchamp,xchamp_g)
2906	IF (is_root_prc) THEN
2907	CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
2908	ENDIF
2909	!-
2910	DO ij=1,nbindex
2911	xchamp(kindex(ij)) = tmaxm1(ij)
2912	ENDDO
2913	var_name= 'tmaxm1'
2914	CALL gather2D(xchamp,xchamp_g)
2915	IF (is_root_prc) THEN
2916	CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
2917	ENDIF
2918	!-
2919	DO ij=1,nbindex
2920	xchamp(kindex(ij)) = tminm1(ij)
2921	ENDDO
2922	var_name= 'tminm1'
2923	CALL gather2D(xchamp,xchamp_g)
2924	IF (is_root_prc) THEN
2925	CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
2926	ENDIF
2927	!-
2928	DO ij=1,nbindex
2929	xchamp(kindex(ij)) = qdm1(ij)
2930	ENDDO
2931	var_name= 'qdm1'
2932	CALL gather2D(xchamp,xchamp_g)
2933	IF (is_root_prc) THEN
2934	CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
2935	ENDIF
2936	!-
2937	DO ij=1,nbindex
2938	xchamp(kindex(ij)) = udm1(ij)
2939	ENDDO
2940	var_name= 'udm1'
2941	CALL gather2D(xchamp,xchamp_g)
2942	IF (is_root_prc) THEN
2943	CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
2944	ENDIF
2945	!-
2946	DO ij=1,nbindex
2947	xchamp(kindex(ij)) = precipm1(ij)
2948	ENDDO
2949	var_name= 'precipm1'
2950	CALL gather2D(xchamp,xchamp_g)
2951	IF (is_root_prc) THEN
2952	CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
2953	ENDIF
2954	!--------------------------------
2955	END SUBROUTINE weathgen_restwrite
2956	!-
2957	!===
2958	!-
2959	SUBROUTINE weather_read &
2960	& (force_id,nomvar,iim_file,jjm_file,n3,i_cut, &
2961	& iim,jjm,n_agg,ncorr,icorr,jcorr,champout)
2962	!---------------------------------------------------------------------
2963	IMPLICIT NONE
2964	!-
2965	INTEGER,INTENT(IN) :: force_id
2966	CHARACTER(LEN=*),INTENT(IN) :: nomvar
2967	INTEGER,INTENT(IN) :: iim_file,jjm_file
2968	INTEGER,INTENT(IN) :: n3
2969	INTEGER,INTENT(IN) :: i_cut
2970	INTEGER,INTENT(IN) :: iim,jjm
2971	INTEGER,INTENT(IN) :: n_agg
2972	INTEGER,DIMENSION(:,:),INTENT(IN) :: ncorr
2973	INTEGER,DIMENSION(:,:,:),INTENT(IN) :: icorr,jcorr
2974	!-
2975	REAL,DIMENSION(iim*jjm,n3),INTENT(OUT) :: champout
2976	!-
2977	REAL,DIMENSION(iim_file,jjm_file,n3) :: champ_file
2978	REAL,ALLOCATABLE,DIMENSION(:,:) :: champout_g
2979	INTEGER :: i,j,ij,l,m
2980	!---------------------------------------------------------------------
2981	WRITE(numout,*) 'Lecture ',TRIM(nomvar)
2982	!-
2983	IF (is_root_prc) THEN
2984	ALLOCATE(champout_g(iim_g*jjm_g,n3))
2985	IF ( n3 == 1 ) THEN
2986	CALL flinget (force_id,nomvar(1:LEN_TRIM(nomvar)), &
2987	& iim_file, jjm_file, 0, 0, 1, 1, champ_file)
2988	ELSE
2989	DO l=1,n3
2990	CALL flinget &
2991	& (force_id,nomvar(1:LEN_TRIM(nomvar)), &
2992	& iim_file, jjm_file, 0, n3, l, l, champ_file(:,:,l))
2993	ENDDO
2994	ENDIF
2995	! shift if necessary
2996	IF (i_cut /= 0) THEN
2997	DO l=1,n3
2998	CALL shift_field (iim_file,jjm_file,i_cut,champ_file(:,:,l))
2999	ENDDO
3000	ENDIF
3001	! interpolate onto the model grid
3002	DO l=1,n3
3003	DO j=1,jjm_g
3004	DO i=1,iim_g
3005	ij = i+(j-1)*iim_g
3006	champout_g(ij,l) = 0.
3007	DO m=1,ncorr(i,j)
3008	champout_g(ij,l) = champout_g(ij,l) &
3009	& +champ_file(icorr(i,j,m),jcorr(i,j,m),l)
3010	ENDDO
3011	champout_g(ij,l) = champout_g(ij,l)/REAL(ncorr(i,j))
3012	ENDDO
3013	ENDDO
3014	ENDDO
3015	!!$ DO l=1,n3
3016	!!$ DO j=1,jjm_g
3017	!!$ WRITE(numout,) j,(/ ( champout_g((j-1)iim_g+i,l), i=1,iim_g ) /)
3018	!!$ ENDDO
3019	!!$ ENDDO
3020	ELSE
3021	ALLOCATE(champout_g(1,1))
3022	ENDIF
3023	!!$ CALL scatter2D(champout_g,champout)
3024	#ifndef CPP_PARA
3025	champout(:,:)=champout_g(:,:)
3026	#else
3027	CALL scatter2D_rgen(champout_g,champout,n3)
3028	#endif
3029
3030	!!$ DO l=1,n3
3031	!!$ DO j=1,jjm
3032	!!$ WRITE(numout,) j,(/ ( champout((j-1)iim_g+i,l), i=1,iim_g ) /)
3033	!!$ ENDDO
3034	!!$ ENDDO
3035	!----------------------------
3036	END SUBROUTINE weather_read
3037	!-
3038	!===
3039	!-
3040	SUBROUTINE weathgen_dump &
3041	& (itau, dt_force, iim, jjm, nbindex, kindex, lrstwrite, &
3042	& swdown, rainf, snowf, tair, u, v, qair, pb, lwdown )
3043	!---------------------------------------------------------------------
3044	IMPLICIT NONE
3045	!-
3046	INTEGER,INTENT(IN) :: itau
3047	REAL,INTENT(IN) :: dt_force
3048	INTEGER,INTENT(IN) :: iim,jjm
3049	INTEGER,INTENT(IN) :: nbindex
3050	INTEGER,DIMENSION(iim*jjm),INTENT(IN) :: kindex
3051	LOGICAL,INTENT(IN) :: lrstwrite
3052	REAL,DIMENSION(iim*jjm),INTENT(IN) :: &
3053	& swdown, rainf, snowf, tair, u, v, qair, pb, lwdown
3054	!-
3055	INTEGER :: iret,ndim
3056	INTEGER,DIMENSION(4) :: corner,edges
3057	REAL,DIMENSION(iim*jjm) :: var_gather
3058	!---------------------------------------------------------------------
3059	!-
3060	! time dimension
3061	!-
3062	iret = NF90_PUT_VAR (dump_id,timestp_id,(/ REAL(itau) /), &
3063	& start=(/ itau /),count=(/ 1 /))
3064	iret = NF90_PUT_VAR (dump_id,time_id,(/ REAL(itau)*dt_force /), &
3065	& start=(/ itau /),count=(/ 1 /))
3066	!-
3067	! 2D variables: pas de dimension verticale
3068	!-
3069	IF (gathered) THEN
3070	ndim = 2
3071	corner(1:2) = (/ 1, itau /)
3072	edges(1:2) = (/ nbindex, 1 /)
3073	ELSE
3074	ndim = 3
3075	corner(1:3) = (/ 1, 1, itau /)
3076	edges(1:3) = (/ iim, jjm, 1 /)
3077	ENDIF
3078	!-
3079	CALL gather_weather (iim*jjm,nbindex,kindex,swdown, var_gather)
3080	iret = NF90_PUT_VAR (dump_id,soldownid, var_gather, &
3081	& start=corner(1:ndim), count=edges(1:ndim))
3082	CALL gather_weather (iim*jjm,nbindex,kindex,rainf, var_gather)
3083	iret = NF90_PUT_VAR (dump_id,rainfid, var_gather, &
3084	& start=corner(1:ndim), count=edges(1:ndim))
3085	CALL gather_weather (iim*jjm,nbindex,kindex,snowf, var_gather)
3086	iret = NF90_PUT_VAR (dump_id,snowfid, var_gather, &
3087	& start=corner(1:ndim), count=edges(1:ndim))
3088	CALL gather_weather (iim*jjm,nbindex,kindex,pb, var_gather)
3089	iret = NF90_PUT_VAR (dump_id,psolid, var_gather, &
3090	& start=corner(1:ndim), count=edges(1:ndim))
3091	CALL gather_weather (iim*jjm,nbindex,kindex,lwdown, var_gather)
3092	iret = NF90_PUT_VAR (dump_id,lwradid, var_gather, &
3093	& start=corner(1:ndim), count=edges(1:ndim))
3094	!-
3095	! 3D variables
3096	!-
3097	IF (gathered) THEN
3098	ndim = 3
3099	corner(1:3) = (/ 1, 1, itau /)
3100	edges(1:3) = (/ nbindex, 1, 1 /)
3101	ELSE
3102	ndim = 4
3103	corner(1:4) = (/ 1, 1, 1, itau /)
3104	edges(1:4) = (/ iim, jjm, 1, 1 /)
3105	ENDIF
3106	!-
3107	CALL gather_weather (iim*jjm,nbindex,kindex,u, var_gather)
3108	iret = NF90_PUT_VAR (dump_id,uid, var_gather, &
3109	& start=corner(1:ndim), count=edges(1:ndim))
3110	CALL gather_weather (iim*jjm,nbindex,kindex,v, var_gather)
3111	iret = NF90_PUT_VAR (dump_id,vid, var_gather, &
3112	& start=corner(1:ndim), count=edges(1:ndim))
3113	CALL gather_weather (iim*jjm,nbindex,kindex,tair, var_gather)
3114	iret = NF90_PUT_VAR (dump_id,tairid, var_gather, &
3115	& start=corner(1:ndim), count=edges(1:ndim))
3116	CALL gather_weather (iim*jjm,nbindex,kindex,qair, var_gather)
3117	iret = NF90_PUT_VAR (dump_id,qairid, var_gather, &
3118	& start=corner(1:ndim), count=edges(1:ndim))
3119	!-
3120	IF (lrstwrite) THEN
3121	iret = NF90_CLOSE (dump_id)
3122	ENDIF
3123	!---------------------------
3124	END SUBROUTINE weathgen_dump
3125	!-
3126	!===
3127	!-
3128	SUBROUTINE gather_weather (iimjjm, nbindex, kindex, var_in, var_out)
3129	!---------------------------------------------------------------------
3130	IMPLICIT NONE
3131	!-
3132	INTEGER,INTENT(IN) :: iimjjm,nbindex
3133	INTEGER,DIMENSION(iimjjm),INTENT(IN) :: kindex
3134	REAL,DIMENSION(iimjjm),INTENT(IN) :: var_in
3135	!-
3136	REAL,DIMENSION(iimjjm),INTENT(OUT) :: var_out
3137	!-
3138	INTEGER :: i
3139	LOGICAL,SAVE :: firstcall = .TRUE.
3140	INTEGER,SAVE :: nb_outside
3141	INTEGER,ALLOCATABLE,SAVE,DIMENSION(:) :: outside
3142	!---------------------------------------------------------------------
3143	IF (firstcall) THEN
3144	!---
3145	!-- determine which points are not in the computational domain and
3146	!-- create a mask for these points
3147	!---
3148	firstcall = .FALSE.
3149
3150	ALLOC_ERR=-1
3151	ALLOCATE(outside(iimjjm), STAT=ALLOC_ERR)
3152	IF (ALLOC_ERR/=0) THEN
3153	WRITE(numout,*) "ERROR IN ALLOCATION of outside : ",ALLOC_ERR
3154	STOP
3155	ENDIF
3156	outside(:) = 0.
3157	nb_outside = 0
3158	DO i=1,iimjjm
3159	IF ( ALL( kindex(:) /= i ) ) THEN
3160	nb_outside = nb_outside+1
3161	outside(nb_outside) = i
3162	ENDIF
3163	ENDDO
3164	ENDIF
3165	!-
3166	IF ( gathered ) THEN
3167	DO i=1,nbindex
3168	var_out(i) = var_in(kindex(i))
3169	ENDDO
3170	ELSE
3171	var_out(:) = var_in(:)
3172	DO i=1,nb_outside
3173	var_out(outside(i)) = undef_sechiba
3174	ENDDO
3175	ENDIF
3176	!--------------------
3177	END SUBROUTINE gather_weather
3178	!-
3179	!===
3180	!-
3181	SUBROUTINE shift_field (im,jm,i_cut,champ)
3182	!---------------------------------------------------------------------
3183	INTEGER,INTENT(IN) :: im,jm,i_cut
3184	REAL,DIMENSION(im,jm),INTENT(INOUT) :: champ
3185	!-
3186	REAL,DIMENSION(im,jm) :: champ_temp
3187	!---------------------------------------------------------------------
3188	IF ( (i_cut >= 1).AND.(i_cut <= im) ) THEN
3189	champ_temp(1:im-i_cut-1,:) = champ(i_cut:im,:)
3190	champ_temp(im-i_cut:im,:) = champ(1:i_cut+1,:)
3191	champ(:,:) = champ_temp(:,:)
3192	ENDIF
3193	!-------------------------
3194	END SUBROUTINE shift_field
3195	!-
3196	!===
3197	!-
3198	SUBROUTINE weathgen_domain_size &
3199	& (limit_west,limit_east,limit_north,limit_south, &
3200	& zonal_res,merid_res,iim,jjm)
3201	!---------------------------------------------------------------------
3202	IMPLICIT NONE
3203	!-
3204	REAL,INTENT(INOUT) :: limit_west,limit_east,limit_north,limit_south
3205	REAL,INTENT(IN) :: zonal_res,merid_res
3206	INTEGER,INTENT(OUT) :: iim,jjm
3207	!---------------------------------------------------------------------
3208	IF (limit_west > limit_east) limit_east = limit_east+360.
3209	!-
3210	IF ( (limit_west >= limit_east) &
3211	& .OR.(limit_east > 360.) &
3212	& .OR.(limit_west < -180.) &
3213	& .OR.(limit_east-limit_west > 360.) ) THEN
3214	WRITE(numout,*) 'PROBLEME longitudes.'
3215	WRITE(numout,*) 'Limites Ouest, Est: ',limit_west,limit_east
3216	STOP
3217	ENDIF
3218	!-
3219	IF ( (limit_south < -90.) &
3220	& .OR.(limit_north > 90.) &
3221	& .OR.(limit_south >= limit_north ) ) THEN
3222	WRITE(numout,*) 'PROBLEME latitudes.'
3223	WRITE(numout,*) 'Limites Nord, Sud: ',limit_north,limit_south
3224	STOP
3225	ENDIF
3226	!-
3227	IF ( (zonal_res <= 0. ) &
3228	& .OR.(zonal_res > limit_east-limit_west) ) THEN
3229	WRITE(numout,*) 'PROBLEME resolution zonale.'
3230	WRITE(numout,*) 'Limites Ouest, Est, Resolution: ', &
3231	& limit_west,limit_east,zonal_res
3232	STOP
3233	ENDIF
3234	!-
3235	IF ( (merid_res <= 0.) &
3236	& .OR.(merid_res > limit_north-limit_south) ) THEN
3237	WRITE(numout,*) 'PROBLEME resolution meridionale.'
3238	WRITE(numout,*) 'Limites Nord, Sud, Resolution: ', &
3239	& limit_north,limit_south,merid_res
3240	STOP
3241	ENDIF
3242	!-
3243	iim = NINT(MAX((limit_east-limit_west)/zonal_res,1.))
3244	jjm = NINT(MAX((limit_north-limit_south)/merid_res,1.))
3245	!-
3246	WRITE(numout,*) 'Domain size: iim, jjm = ', iim, jjm
3247	!----------------------------------
3248	END SUBROUTINE weathgen_domain_size
3249	!-
3250	!===
3251	!-
3252	FUNCTION tsatl (t) RESULT (tsat)
3253	!---------------------------------------------------------------------
3254	! statement functions tsatl,tsati are used below so that lowe's
3255	! polyomial for liquid is used if t gt 273.16, or for ice if
3256	! t lt 273.16. also impose range of validity for lowe's polys.
3257	!---------------------------------------------------------------------
3258	REAL,INTENT(IN) :: t
3259	REAL :: tsat
3260	!---------------------------------------------------------------------
3261	tsat = MIN(100.,MAX(t-zero_t,0.))
3262	!-----------------
3263	END FUNCTION tsatl
3264	!-
3265	!===
3266	!-
3267	FUNCTION tsati (t) RESULT (tsat)
3268	!---------------------------------------------------------------------
3269	! statement functions tsatl,tsati are used below so that lowe's
3270	! polyomial for liquid is used if t gt 273.16, or for ice if
3271	! t lt 273.16. also impose range of validity for lowe's polys.
3272	!---------------------------------------------------------------------
3273	REAL,INTENT(IN) :: t
3274	REAL :: tsat
3275	!---------------------------------------------------------------------
3276	tsat = MAX(-60.,MIN(t-zero_t,0.))
3277	!-----------------
3278	END FUNCTION tsati
3279	!-
3280	!===
3281	!-
3282	FUNCTION esat (t) RESULT (esatout)
3283	!---------------------------------------------------------------------
3284	! statement function esat is svp in n/m**2, with t in deg k.
3285	! (100 * lowe's poly since 1 mb = 100 n/m**2.)
3286	!---------------------------------------------------------------------
3287	REAL,INTENT(IN) :: t
3288	REAL :: esatout
3289	REAL :: x
3290	!-
3291	! polynomials for svp(t), d(svp)/dt over water and ice are from
3292	! lowe(1977),jam,16,101-103.
3293	!-
3294	REAL,PARAMETER :: &
3295	& asat0 = 6.1078000, asat1 = 4.4365185e-1, asat2 = 1.4289458e-2, &
3296	& asat3 = 2.6506485e-4, asat4 = 3.0312404e-6, asat5 = 2.0340809e-8, &
3297	& asat6 = 6.1368209e-11, &
3298	& bsat0 = 6.1091780, bsat1 = 5.0346990e-1, bsat2 = 1.8860134e-2, &
3299	& bsat3 = 4.1762237e-4, bsat4 = 5.8247203e-6, bsat5 = 4.8388032e-8, &
3300	& bsat6 = 1.8388269e-10
3301	!---------------------------------------------------------------------
3302	IF (t >= zero_t) THEN
3303	x = asat0
3304	ELSE
3305	x = bsat0
3306	ENDIF
3307	!-
3308	esatout = 100.* &
3309	( x &
3310	+tsatl(t)(asat1+tsatl(t)(asat2+tsatl(t)*(asat3 &
3311	+tsatl(t)(asat4+tsatl(t)(asat5+tsatl(t)* asat6))))) &
3312	+tsati(t)(bsat1+tsati(t)(bsat2+tsati(t)*(bsat3 &
3313	+tsati(t)(bsat4+tsati(t)(bsat5+tsati(t)* bsat6))))) )
3314	!----------------
3315	END FUNCTION esat
3316	!-
3317	!===
3318	!-
3319	FUNCTION qsat (e,p) RESULT (qsatout)
3320	!---------------------------------------------------------------------
3321	! statement function qsat is saturation specific humidity,
3322	! with svp e and ambient pressure p in n/m**2. impose an upper
3323	! limit of 1 to avoid spurious values for very high svp
3324	! and/or small p
3325	!---------------------------------------------------------------------
3326	REAL, INTENT(IN) :: e,p
3327	REAL :: qsatout
3328	!---------------------------------------------------------------------
3329	qsatout = 0.622e/MAX(p-(1.0-0.622)e,0.622*e)
3330	!----------------
3331	END FUNCTION qsat
3332	!-
3333	!===
3334	!-
3335	SUBROUTINE weathgen_qsat (npoi,t,p,qsat)
3336	!---------------------------------------------------------------------
3337	! vectorized version of functions esat and qsat.
3338	! statement function esat is svp in n/m**2, with t in deg k.
3339	! (100 * lowe's poly since 1 mb = 100 n/m**2.)
3340	!---------------------------------------------------------------------
3341	INTEGER,INTENT(IN) :: npoi
3342	REAL,DIMENSION(npoi),INTENT(IN) :: t,p
3343	REAL,DIMENSION(npoi),INTENT(OUT):: qsat
3344	!-
3345	REAL,DIMENSION(npoi) :: x, tl, ti, e
3346	!-
3347	! polynomials for svp(t), d(svp)/dt over water and ice
3348	! are from lowe(1977),jam,16,101-103.
3349	!-
3350	REAL,PARAMETER :: &
3351	& asat0 = 6.1078000, asat1 = 4.4365185e-1, asat2 = 1.4289458e-2, &
3352	& asat3 = 2.6506485e-4, asat4 = 3.0312404e-6, asat5 = 2.0340809e-8, &
3353	& asat6 = 6.1368209e-11, &
3354	& bsat0 = 6.1091780, bsat1 = 5.0346990e-1, bsat2 = 1.8860134e-2, &
3355	& bsat3 = 4.1762237e-4, bsat4 = 5.8247203e-6, bsat5 = 4.8388032e-8, &
3356	& bsat6 = 1.8388269e-10
3357	!---------------------------------------------------------------------
3358	WHERE (t(:) > zero_t)
3359	x(:) = asat0
3360	ELSEWHERE
3361	x(:) = bsat0
3362	ENDWHERE
3363	!-
3364	tl(:) = MIN(100.,MAX(t(:)-zero_t,0.))
3365	ti(:) = MAX(-60.,MIN(t(:)-zero_t,0.))
3366	!-
3367	e(:) = 100.* &
3368	( x(:) &
3369	+tl(:)(asat1+tl(:)(asat2+tl(:)*(asat3 &
3370	+tl(:)(asat4+tl(:)(asat5+tl(:)* asat6))))) &
3371	+ti(:)(bsat1+ti(:)(bsat2+ti(:)*(bsat3 &
3372	+ti(:)(bsat4+ti(:)(bsat5+ti(:)* bsat6))))) )
3373	!-
3374	qsat(:) = 0.622e(:)/MAX(p(:)-(1.0-0.622)e(:),0.622*e(:))
3375	!---------------------------
3376	END SUBROUTINE weathgen_qsat
3377	!-
3378	!===
3379	!-
3380	SUBROUTINE mask_c_o &
3381	& (imdep, jmdep, xdata, ydata, mask_in, fcrit, &
3382	& imar, jmar, zonal_res, merid_res, n_agg, x, y, mask, &
3383	& ncorr, icorr, jcorr)
3384	!---------------------------------------------------------------------
3385	! z.x.li (le 01/04/1994) :
3386	! A partir du champ de masque, on fabrique
3387	! un champ indicateur (masque) terre/ocean
3388	! terre:1; ocean:0
3389	!---------------------------------------------------------------------
3390	INTEGER :: imdep,jmdep
3391	REAL :: xdata(imdep),ydata(jmdep)
3392	REAL :: mask_in(imdep,jmdep)
3393	REAL :: fcrit
3394	INTEGER :: imar,jmar
3395	REAL :: zonal_res,merid_res
3396	INTEGER :: n_agg
3397	REAL :: x(imar),y(jmar)
3398	REAL, INTENT(OUT) :: mask(imar,jmar)
3399	INTEGER :: ncorr(imar,jmar)
3400	INTEGER,DIMENSION(imar,jmar,n_agg) :: icorr,jcorr
3401	!-
3402	INTEGER i, j, ii, jj
3403	REAL a(imar),b(imar),c(jmar),d(jmar)
3404	INTEGER num_tot(imar,jmar), num_oce(imar,jmar)
3405	REAL,ALLOCATABLE :: distans(:)
3406	INTEGER ij_proche(1),i_proche,j_proche
3407	!-
3408	INTEGER,DIMENSION(imar,jmar) :: ncorr_oce , ncorr_land
3409	INTEGER,DIMENSION(imar,jmar,n_agg) :: &
3410	& icorr_oce, jcorr_oce , icorr_land, jcorr_land
3411	!-
3412	INTEGER imdepp
3413	REAL,ALLOCATABLE :: xdatap(:)
3414	REAL,ALLOCATABLE :: mask_inp(:,:)
3415	LOGICAL :: extend
3416	!---------------------------------------------------------------------
3417	ncorr(:,:) = 0
3418	icorr(:,:,:) = -1; jcorr(:,:,:) = -1
3419	ncorr_land(:,:) = 0
3420	icorr_land(:,:,:) = -1; jcorr_land(:,:,:) = -1
3421	ncorr_oce(:,:) = 0
3422	icorr_oce(:,:,:) = -1; jcorr_oce(:,:,:) = -1
3423	! do we have to extend the domain (-x...-x+360)?
3424	IF ( xdata(1)+360.-xdata(imdep) > 0.01 ) THEN
3425	extend = .TRUE.
3426	imdepp = imdep+1
3427	ELSE
3428	extend = .FALSE.
3429	imdepp = imdep
3430	ENDIF
3431	!-
3432
3433	ALLOC_ERR=-1
3434	ALLOCATE(xdatap(imdepp), STAT=ALLOC_ERR)
3435	IF (ALLOC_ERR/=0) THEN
3436	WRITE(numout,*) "ERROR IN ALLOCATION of xdatap : ",ALLOC_ERR
3437	STOP
3438	ENDIF
3439
3440	ALLOC_ERR=-1
3441	ALLOCATE(mask_inp(imdepp,jmdep), STAT=ALLOC_ERR)
3442	IF (ALLOC_ERR/=0) THEN
3443	WRITE(numout,*) "ERROR IN ALLOCATION of mask_inp : ",ALLOC_ERR
3444	STOP
3445	ENDIF
3446	!-
3447	xdatap(1:imdep) = xdata(1:imdep)
3448	mask_inp(1:imdep,:) = mask_in(1:imdep,:)
3449	!-
3450	IF (extend) THEN
3451	xdatap(imdepp) = xdatap(1)+360.
3452	mask_inp(imdepp,:) = mask_inp(1,:)
3453	ENDIF
3454	!-
3455
3456	ALLOC_ERR=-1
3457	ALLOCATE(distans(imdepp*jmdep), STAT=ALLOC_ERR)
3458	IF (ALLOC_ERR/=0) THEN
3459	WRITE(numout,*) "ERROR IN ALLOCATION of distans : ",ALLOC_ERR
3460	STOP
3461	ENDIF
3462	! Definition des limites des boites de la grille d'arrivee.
3463	IF (imar > 1) THEN
3464	a(1) = x(1)-(x(2)-x(1))/2.0
3465	b(1) = (x(1)+x(2))/2.0
3466	DO i=2,imar-1
3467	a(i) = b(i-1)
3468	b(i) = (x(i)+x(i+1))/2.0
3469	ENDDO
3470	a(imar) = b(imar-1)
3471	b(imar) = x(imar)+(x(imar)-x(imar-1))/2.0
3472	ELSE
3473	a(1) = x(1)-zonal_res/2.
3474	b(1) = x(1)+zonal_res/2.
3475	ENDIF
3476	!-
3477	IF (jmar > 1) THEN
3478	c(1) = y(1)-(y(2)-y(1))/2.0
3479	d(1) = (y(1)+y(2))/2.0
3480	DO j=2,jmar-1
3481	c(j) = d(j-1)
3482	d(j) = (y(j)+y(j+1))/2.0
3483	ENDDO
3484	c(jmar) = d(jmar-1)
3485	d(jmar) = y(jmar)+(y(jmar)-y(jmar-1))/2.0
3486	ELSE
3487	c(1) = y(1)-merid_res/2.
3488	d(1) = y(1)+merid_res/2.
3489	ENDIF
3490	!-
3491	num_oce(1:imar,1:jmar) = 0
3492	num_tot(1:imar,1:jmar) = 0
3493	!-
3494	! ..... Modif P. Le Van ( 23/08/95 ) ....
3495	!-
3496	DO ii=1,imar
3497	DO jj=1,jmar
3498	DO i=1,imdepp
3499	IF ( ( (xdatap(i)-a(ii) >= 1.e-5) &
3500	& .AND.(xdatap(i)-b(ii) <= 1.e-5) ) &
3501	& .OR.( (xdatap(i)-a(ii) <= 1.e-5) &
3502	& .AND.(xdatap(i)-b(ii) >= 1.e-5) ) ) THEN
3503	DO j=1,jmdep
3504	IF ( ( (ydata(j)-c(jj) >= 1.e-5) &
3505	& .AND.(ydata(j)-d(jj) <= 1.e-5) ) &
3506	& .OR.( (ydata(j)-c(jj) <= 1.e-5) &
3507	& .AND.(ydata(j)-d(jj) >= 1.e-5) ) ) THEN
3508	num_tot(ii,jj) = num_tot(ii,jj)+1
3509	IF (mask_inp(i,j) < 0.5) THEN
3510	num_oce(ii,jj) = num_oce(ii,jj)+1
3511	!-------------- on a trouve un point oceanique. On le memorise.
3512	ncorr_oce(ii,jj) = ncorr_oce(ii,jj)+1
3513	IF ((i == imdepp).AND.extend) THEN
3514	icorr_oce(ii,jj,ncorr_oce(ii,jj)) = 1
3515	ELSE
3516	icorr_oce(ii,jj,ncorr_oce(ii,jj)) = i
3517	ENDIF
3518	jcorr_oce(ii,jj,ncorr_oce(ii,jj)) = j
3519	ELSE
3520	!-------------- on a trouve un point continental. On le memorise.
3521	ncorr_land(ii,jj) = ncorr_land(ii,jj)+1
3522	IF ((i == imdepp).AND.extend) THEN
3523	icorr_land(ii,jj,ncorr_land(ii,jj)) = 1
3524	ELSE
3525	icorr_land(ii,jj,ncorr_land(ii,jj)) = i
3526	ENDIF
3527	jcorr_land(ii,jj,ncorr_land(ii,jj)) = j
3528	ENDIF
3529	ENDIF
3530	ENDDO
3531	ENDIF
3532	ENDDO
3533	ENDDO
3534	ENDDO
3535	!-
3536	DO i=1,imar
3537	DO j=1,jmar
3538	IF (num_tot(i,j) > 0) THEN
3539	IF ( ( (num_oce(i,j) == 0) &
3540	& .AND.(num_tot(i,j) > 0) ) &
3541	& .OR.( (num_oce(i,j) > 0) &
3542	& .AND.( REAL(num_oce(i,j)) &
3543	& <= REAL(num_tot(i,j))*(1.-fcrit) ) ) ) THEN
3544	mask(i,j) = 1.
3545	ncorr(i,j) = ncorr_land(i,j)
3546	icorr(i,j,:) = icorr_land(i,j,:)
3547	jcorr(i,j,:) = jcorr_land(i,j,:)
3548	ELSE
3549	mask(i,j) = 0.
3550	ncorr(i,j) = ncorr_oce(i,j)
3551	icorr(i,j,:) = icorr_oce(i,j,:)
3552	jcorr(i,j,:) = jcorr_oce(i,j,:)
3553	ENDIF
3554	ELSE
3555	CALL dist_sphe(x(i),y(j),xdatap,ydata,imdepp,jmdep,distans)
3556	ij_proche(:) = MINLOC(distans)
3557	j_proche = (ij_proche(1)-1)/imdepp+1
3558	i_proche = ij_proche(1)-(j_proche-1)*imdepp
3559	mask(i,j) = mask_inp(i_proche,j_proche)
3560	IF ( (i_proche == imdepp).AND.extend) i_proche=1
3561	ncorr(i,j) = 1
3562	icorr(i,j,1) = i_proche
3563	jcorr(i,j,1) = j_proche
3564	ENDIF
3565	ENDDO
3566	ENDDO
3567	!----------------------
3568	END SUBROUTINE mask_c_o
3569	!-
3570	!===
3571	!-
3572	SUBROUTINE dist_sphe (rf_lon,rf_lat,rlon,rlat,im,jm,distance)
3573	!---------------------------------------------------------------------
3574	! Auteur: Laurent Li (le 30/12/1996)
3575	!
3576	! Ce programme calcule la distance minimale (selon le grand cercle)
3577	! entre deux points sur la terre
3578	!---------------------------------------------------------------------
3579	!-
3580	! Input:
3581	!-
3582	INTEGER im, jm ! dimensions
3583	REAL rf_lon ! longitude du point de reference (degres)
3584	REAL rf_lat ! latitude du point de reference (degres)
3585	REAL rlon(im), rlat(jm) ! longitude et latitude des points
3586	!-
3587	! Output:
3588	!-
3589	REAL distance(im,jm) ! distances en metre
3590	!-
3591	REAL rlon1, rlat1
3592	REAL rlon2, rlat2
3593	REAL dist
3594	REAL pa, pb, p, pi
3595	INTEGER i,j
3596	!-
3597	REAL radius
3598	PARAMETER (radius=6371229.)
3599	!---------------------------------------------------------------------
3600	pi = 4.0*ATAN(1.0)
3601	!-
3602	DO j=1,jm
3603	DO i=1,im
3604	rlon1=rf_lon
3605	rlat1=rf_lat
3606	rlon2=rlon(i)
3607	rlat2=rlat(j)
3608	pa = pi/2.0-rlat1*pi/180.0 ! dist. entre pole n et point a
3609	pb = pi/2.0-rlat2*pi/180.0 ! dist. entre pole n et point b
3610	!-----
3611	p = (rlon1-rlon2)*pi/180.0 ! angle entre a et b (leurs meridiens)
3612	!-----
3613	dist = ACOS( COS(pa)COS(pb)+SIN(pa)SIN(pb)*COS(p))
3614	dist = radius*dist
3615	distance(i,j) = dist
3616	ENDDO
3617	ENDDO
3618	!-----------------------
3619	END SUBROUTINE dist_sphe
3620	!-
3621	!===
3622	!-
3623	SUBROUTINE permute (n,ordre)
3624	!---------------------------------------------------------------------
3625	INTEGER,INTENT(IN) :: n
3626	INTEGER,DIMENSION(n),INTENT(OUT) :: ordre
3627	!-
3628	INTEGER,DIMENSION(n) :: restant
3629	INTEGER :: ipique, i, n_rest
3630	REAL :: rndnum
3631	!---------------------------------------------------------------------
3632	n_rest = n
3633	restant(:) = (/ (i, i=1,n) /)
3634	!-
3635	DO i=1,n
3636	CALL random_number (rndnum)
3637	ipique = INT(rndnum*n_rest)+1
3638	ordre(i) = restant(ipique)
3639	restant(ipique:n_rest-1) = restant(ipique+1:n_rest)
3640	n_rest = n_rest-1
3641	ENDDO
3642	!---------------------
3643	END SUBROUTINE permute
3644	!-
3645	!===
3646	!-----------------
3647	END MODULE weather

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

Download in other formats: