Context Navigation

weather.f90 @ 7262

Last change on this file since 7262 was 4646, checked in by josefine.ghattas, 7 years ago

Ticket #185 :
Created new module containing time variables and subroutines to calculate them. Variables in the time module are public and can be used elswhere in the model. They should not be modified outside the time module.

Note that now each time step have a time interval. xx_start or xx_end values for the date can be used. Previously the date values correspondend to the end of the interval. See description in module time for more information.

These changes do not make any difference in results for couled mode with LMDZ or dim2_driver offline driver.

These modifications also corrects the problem with do_slow related to orchideedriver, ticket #327

File size: 122.2 KB

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

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

Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/weather.f90 @ 7262

Download in other formats: