Context Navigation

source: tags/ORCHIDEE_1_9_5_1/ORCHIDEE_OL/weather.f90 @ 118

Last change on this file since 118 was 48, checked in by mmaipsl, 14 years ago

MM: Tests with lf95 compiler : correct f95 strict norm problems.

There is no change in numerical result after these commits.

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

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

Download in other formats: