Context Navigation

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

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

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

Download in other formats: