Context Navigation

source: tags/ORCHIDEE_1_9_6/ORCHIDEE_OL/weather.f90 @ 3850

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

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

Download in other formats: