Context Navigation

source: branches/publications/ORCHIDEE_CAN_r2290/src_driver/weather.f90

Last change on this file was 1586, checked in by matthew.mcgrath, 11 years ago
DEV: Compiles with NAG but does not link yet
File size: 122.1 KB

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

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

Download in other formats: