Context Navigation

weather.f90

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

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

Context Navigation

source: tags/ORCHIDEE_1_9_5/ORCHIDEE_OL/weather.f90

Download in other formats: