Context Navigation

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

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

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

Download in other formats: