Context Navigation

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

Last change on this file since 382 was 258, checked in by didier.solyga, 13 years ago
Externalized version of ORCHIDEE_OL files merged with the trunk
File size: 108.5 KB

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

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

Download in other formats: