Context Navigation

weather.f90 @ 7329

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

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

Context Navigation

source: branches/publications/ORCHIDEE_CAN_r3069/src_driver/weather.f90 @ 7329

Download in other formats: