Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

cotes.f90_20170616 @ 3928

Last change on this file since 3928 was 3326, checked in by omamce, 7 years ago
O.M. : Utility to generate interpolatio weights for OASIS-MCT
File size: 64.2 KB

Line
1	! -- Mode: f90 --
2	PROGRAM cote
3	USE declare
4	USE modeles
5	USE mod_lecdata
6	USE formula
7	USE math
8	USE mod_rectifie
9	USE mod_lbc
10	USE mod_norma
11	USE bords
12	USE mod_inter
13	USE mod_prih
14	USE mod_wri_wei
15	USE mod_inipar
16	USE fliocom
17	USE getincom
18	USE errioipsl
19	!!!
20	IMPLICIT NONE
21	!!!
22	!!! A partir d'un fichier de poids, complete avec les rivieres et
23	!!! les poids pour le run off
24	!! On met a zero les poids sur les points non cotiers, et on renormalise
25	!!
26	INTEGER (kind=il), PARAMETER :: jpr = 52 ! Nombre de rivieres
27	REAL (kind=rl) :: zsuma, zsumo
28	REAL (kind=rl), DIMENSION (:), ALLOCATABLE :: za
29	REAL (kind=rl), DIMENSION (:), ALLOCATABLE :: zo
30	!!
31	CHARACTER (LEN = 1) :: clriv
32	CHARACTER (LEN = 30) :: clname
33	INTEGER (kind=il) :: ja, jo, jo2, jo3, jn, jno, joi, joj, jai, jaj, jr, kn ! indices de boucle
34	INTEGER (kind=il) :: jai_p, jai_m, jaj_p, jaj_m, n1
35	INTEGER (kind=il) :: joi_p, joi_m, joj_p, joj_m
36	REAL (kind=rl) :: z_d_1, z_d_2, zzmin
37	INTEGER (kind=il) :: nriv = 32, nlmd = 35
38	INTEGER (kind=il), DIMENSION (:, :), ALLOCATABLE :: jx ! Point trouve sous une maille atm
39	INTEGER (kind=il) :: kom, ko
40	LOGICAL, DIMENSION (:), ALLOCATABLE :: lacot, laland, laoce, lanoroute
41	LOGICAL, DIMENSION (:), ALLOCATABLE :: locot, looce
42	INTEGER (kind=il), DIMENSION (:), ALLOCATABLE :: ktemp
43	REAL(kind=rl), DIMENSION (:), ALLOCATABLE :: wtemp, d_cote
44	INTEGER (kind=il), DIMENSION (:,:), ALLOCATABLE :: ktemp2
45	REAL(kind=rl), DIMENSION (:,:), ALLOCATABLE :: wtemp2
46	INTEGER (kind=il), DIMENSION (:), ALLOCATABLE :: iamsk
47	INTEGER (kind=il), DIMENSION (:), ALLOCATABLE :: iomsk
48	INTEGER (kind=il) :: jovoid, jot, inum
49	CHARACTER (len=80) :: cfname4, cfname8
50	LOGICAL :: ltrouve, ltrouve2
51	INTEGER (kind=il), DIMENSION (:), ALLOCATABLE :: itarget
52	REAL (kind=rl), DIMENSION (:), ALLOCATABLE :: wtarget
53	INTEGER (kind=il), DIMENSION (:,:), ALLOCATABLE :: i2a
54	INTEGER (kind=il) :: no, ntrouve
55	INTEGER (kind=il), DIMENSION (:), ALLOCATABLE :: jo_n
56	LOGICAL :: l_fulldiag = .TRUE.
57	CHARACTER (len=180) :: c_date, c_time, c_zone, c_tmp
58	CHARACTER (len=80) :: clarg
59	INTEGER :: narg
60	REAL(kind=rl), DIMENSION (:), ALLOCATABLE :: tmask_atl, tmask_noclose, tmask_nomed, tmask_med
61	REAL(kind=rl), DIMENSION (:, :), ALLOCATABLE :: z2o
62	CHARACTER (len=40) :: cl_bassin
63	INTEGER (kind=il) :: il_ncid
64	!!
65	INTEGER (kind=il), DIMENSION (2) :: jind
66	INTEGER (kind=il) :: ierr
67	!!
68	!! Recherche voisins proches
69	INTEGER (kind=il) :: jp_nb
70	INTEGER (kind=il), DIMENSION (:), ALLOCATABLE :: jo_tab
71	INTEGER (kind=il) :: jo_nb
72	!!
73	INTEGER :: nid_it
74	!!
75	CHARACTER (len = 50) :: cldiag_a2o, clw_a2o, clw_a2o_cdf, clw_a2o_mct
76	CHARACTER (len = 25) :: mo_name, a2o_name
77	CHARACTER (len = 1) :: c_i, c_r
78
79	!! Lecture de la ligne de comande
80	narg = iargc()
81	IF ( narg > 0 ) THEN
82	CALL getarg ( 1, clarg)
83	SELECT CASE (TRIM(clarg))
84	CASE Default
85	IF ( INDEX ( TRIM(clarg), '.def' ) /= 0 ) THEN
86	WRITE (nout,*) 'Lecture des parametres dans ', TRIM(clarg)
87	CALL getin_name ( TRIM (clarg) )
88	ELSE
89	WRITE (nout,*) 'Lecture des parametres dans run.def'
90	ENDIF
91	END SELECT
92	END IF
93	!!
94	CALL inipar
95	CALL alloc_modeles
96
97	!!
98	ALLOCATE (za (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'za', lreset=.TRUE., crout='cotes')
99	ALLOCATE (zo (jpon) , STAT=ierr) ; CALL chk_allo (ierr, 'zo')
100	ALLOCATE (jx (jpo2a,2) , STAT=ierr) ; CALL chk_allo (ierr, 'jx')
101	ALLOCATE (lacot (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'lacot') ; lacot = .FALSE.
102	ALLOCATE (laland (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'laland') ; laland = .FALSE.
103	ALLOCATE (laoce (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'laoce') ; laoce = .FALSE.
104	!
105	ALLOCATE (locot (jpon) , STAT=ierr) ; CALL chk_allo (ierr, 'locot') ; locot = .FALSE.
106	ALLOCATE (looce (jpon) , STAT=ierr) ; CALL chk_allo (ierr, 'looce') ; looce = .FALSE.
107	!
108	ALLOCATE (ktemp (jpa2o) , STAT=ierr) ; CALL chk_allo (ierr, 'ktemp')
109	ALLOCATE (wtemp (jpa2o) , STAT=ierr) ; CALL chk_allo (ierr, 'wtemp')
110	ALLOCATE (iamsk (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'iamsk')
111	ALLOCATE (iomsk (jpon) , STAT=ierr) ; CALL chk_allo (ierr, 'iomsk')
112	ALLOCATE (d_cote (jpon) , STAT=ierr) ; CALL chk_allo (ierr, 'd_cote')
113	!
114	ALLOCATE (itarget (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'itarget')
115	ALLOCATE (wtarget (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'wtarget')
116	ALLOCATE (i2a (jpai, jpaj), STAT=ierr) ; CALL chk_allo (ierr, 'i2a')
117	!
118	ALLOCATE (tmask_atl (jpon), STAT=ierr) ; CALL chk_allo (ierr, 'tmask_atl')
119	ALLOCATE (tmask_noclose (jpon), STAT=ierr) ; CALL chk_allo (ierr, 'tmask_noclose')
120	ALLOCATE (tmask_nomed (jpon), STAT=ierr) ; CALL chk_allo (ierr, 'tmasl_nomed')
121	ALLOCATE (tmask_med (jpon), STAT=ierr) ; CALL chk_allo (ierr, 'tmask_med')
122	ALLOCATE (z2o (jpoi, jpoj), STAT=ierr) ; CALL chk_allo (ierr, 'z2o')
123	!
124	jp_nb = jpon
125	ALLOCATE (jo_tab (jpon)) ; CALL chk_allo (ierr, 'jo_tab')
126	!
127	IF (l_limit_iosize) l_fulldiag = .FALSE.
128	!!
129	!! Initialisation
130	ngrd = 11 ; nsrf = 12 ; nmsk = 13 ; nchk = 14 ; ndeb = 9 ; nbug1 = 10
131	nwei4o2a = 81 ; nwei8o2a = 82 ; nwei4a2o = 83 ; nwei8a2o = 84
132	!
133	nout = 6
134	IF (nout /= 6 .AND. nout /= 0 ) &
135	& OPEN (unit = nout, file = 'cote' // TRIM(c_suffix) // '.out', action = 'WRITE', status = 'unknown', position = 'REWIND' )
136	!!
137	!! Diagnostics files
138	!!
139	IF (nout /= 6 ) &
140	& OPEN (unit = nout, file = 'cotes' // TRIM(c_suffix) // '.out', position = 'REWIND', action = 'WRITE', status = 'replace' )
141	!!
142	!! Read coordinates of all models, and weights
143	!!
144	CALL lecdata (lecweights = .TRUE., leco2amask = .TRUE. )
145
146	!! Lecture masques de bassin
147	IF ( TRIM (cotyp) == 'orca4' .OR. &
148	& TRIM (cotyp) == 'orca2' .OR. TRIM (cotyp) == 'orca2.0' .OR. TRIM (cotyp) == 'orca2.1' .OR. TRIM (cotyp) == 'orca2.3' .OR. &
149	& TRIM (cotyp) == 'orca1' &
150	& ) THEN
151	cl_bassin = TRIM (cotyp)
152	!-$$ IF (c_period /= 'none' ) cl_bassin = TRIM (cl_bassin) // TRIM (c_period)
153	cl_bassin = TRIM (cl_bassin) // '.nc'
154	WRITE (nout,*) 'Reading bassin masks in : ', TRIM (cl_bassin)
155	CALL flioopfd (TRIM (cl_bassin), il_ncid )
156	CALL fliogetv (il_ncid, 'mask_atl' , z2o )
157	tmask_atl = RESHAPE (z2o, (/ jpon /) )
158	CALL fliogetv (il_ncid, 'mask_noclose', z2o )
159	tmask_noclose = RESHAPE (z2o, (/ jpon /) )
160	CALL fliogetv (il_ncid, 'mask_nomed' , z2o )
161	tmask_nomed = RESHAPE (z2o, (/ jpon /) )
162	CALL flioclo (il_ncid)
163	ELSE
164	tmask_atl = REAL (1-iomskt, rl)
165	tmask_noclose = REAL (1-iomskt, rl)
166	tmask_nomed = REAL (1-iomskt, rl)
167	END IF
168	!
169	WRITE (unit=nout, fmt=*) 'iomskt'
170	CALL prihin (RESHAPE(iomskt ,(/jpoi,jpoj/)), 1_il, nout)
171	WRITE (unit=nout, fmt=*) 'tmask_atl'
172	CALL prihin (RESHAPE(NINT(tmask_atl) ,(/jpoi,jpoj/)), 1_il, nout)
173	WRITE (unit=nout, fmt=*) 'tmask_med'
174	!CALL prihin (RESHAPE(NINT(tmask_med) ,(/jpoi,jpoj/)), 1_il, nout)
175	WRITE (unit=nout, fmt=*) 'tmask_noclose'
176	!CALL prihin (RESHAPE(NINT(tmask_noclose),(/jpoi,jpoj/)), 1_il, nout)
177	!
178	! Set correct periodicity to grids
179	CALL lbc (xolont, noperio, 'T', jpoi, jpoj)
180	CALL lbc (xolonu, noperio, 'U', jpoi, jpoj)
181	CALL lbc (xolonv, noperio, 'V', jpoi, jpoj)
182	CALL lbc (xolonf, noperio, 'F', jpoi, jpoj)
183	CALL lbc (xolatt, noperio, 'T', jpoi, jpoj)
184	CALL lbc (xolatu, noperio, 'U', jpoi, jpoj)
185	CALL lbc (xolatv, noperio, 'V', jpoi, jpoj)
186	CALL lbc (xolatf, noperio, 'F', jpoi, jpoj)
187	CALL lbc (xosrft, noperio, 'T', jpoi, jpoj)
188	CALL lbc (iomskt, noperio, 'T', jpoi, jpoj)
189	!
190	CALL lbc (xalont, naperio, 'T', jpai, jpaj)
191	CALL lbc (xalonu, naperio, 'U', jpai, jpaj)
192	CALL lbc (xalonv, naperio, 'V', jpai, jpaj)
193	CALL lbc (xalatt, naperio, 'T', jpai, jpaj)
194	CALL lbc (xalatu, naperio, 'U', jpai, jpaj)
195	CALL lbc (xalatv, naperio, 'V', jpai, jpaj)
196	CALL lbc (xasrft, naperio, 'T', jpai, jpaj)
197	CALL lbc (iamskt, naperio, 'T', jpai, jpaj)
198	!!
199	CALL rectifie
200	!!
201	CALL fliocrfd ('itarget' // TRIM(c_suffix), (/'x', 'y'/), (/jpai, jpaj/), nid_it, mode='REPLACE')
202	CALL flioputa (nid_it, "?", "title" , "itarget" )
203	CALL flioputa (nid_it, '?', 'Comment', TRIM(c_comment) )
204	CALL DATE_AND_TIME (c_date, c_time, c_zone )
205	CALL flioputa (nid_it, "?", "history" , "Created: "//c_date(1:4)//"-"//c_date(5:6)//"-"// &
206	& c_date(7:8)//" "//c_time(1:2)//"h"//c_time(3:4)//" GMT"//TRIM(c_zone) )
207	CALL flioputa (nid_it, "?", "method" , "MOSAIC" )
208	CALL flioputa (nid_it, "?", "source_grid" , "curvilinear" )
209	CALL flioputa (nid_it, "?", "dest_grid" , "curvilinear" )
210	CALL flioputa (nid_it, "?", "Institution" , "IPSL" )
211	CALL flioputa (nid_it, "?", "Model" , "IPSL CM6" )
212	CALL GET_ENVIRONMENT_VARIABLE ( NAME="HOSTNAME", VALUE=c_tmp, TRIM_NAME=.TRUE., STATUS=ierr)
213	IF ( ierr == 0 ) THEN
214	CALL flioputa (nid_it, "?", "HOSTNAME" , TRIM(c_tmp) )
215	ELSE
216	WRITE (nout,*) 'Environment variable not found : $HOSTNAME'
217	END IF
218	CALL GET_ENVIRONMENT_VARIABLE ( NAME="LOGNAME" , VALUE=c_tmp, TRIM_NAME=.TRUE., STATUS=ierr)
219	IF ( ierr == 0 ) THEN
220	CALL flioputa (nid_it, "?", "LOGNAME" , TRIM(c_tmp) )
221	ELSE
222	WRITE (nout,*) 'Environment variable not found : $LOGNAME'
223	END IF
224
225	CALL fliodefv (nid_it, 'lon', (/1, 2/), units="degree_north", v_t=flio_r4 )
226	CALL fliodefv (nid_it, 'lat', (/1, 2/), units="degree_east", v_t=flio_r4 )
227	CALL fliodefv (nid_it, 'o2amask', (/1,2/), v_t=flio_r4)
228	CALL fliodefv (nid_it, 'laland', (/1,2/), v_t=flio_i4)
229	CALL fliodefv (nid_it, 'lacot' , (/1,2/), v_t=flio_i4)
230	CALL fliodefv (nid_it, 'laoce' , (/1,2/), v_t=flio_i4)
231	CALL fliodefv (nid_it, 'iamskp', (/1,2/), v_t=flio_i4)
232
233	CALL flioputv (nid_it, 'lon' , RESHAPE ( xalont, (/jpai, jpaj/)) )
234	CALL flioputv (nid_it, 'lat' , RESHAPE ( xalatt, (/jpai, jpaj/)) )
235	CALL flioputv (nid_it, 'o2amask', RESHAPE (o2amask, (/jpai, jpaj/)) )
236	CALL flioputv (nid_it, 'iamskp' , RESHAPE (iamskp , (/jpai, jpaj/)) )
237	!!
238	!! Compute edges of grid boxes
239	!!
240	CALL bord_a
241	CALL bord_o
242	CALL bord_oa
243	!!
244	mo_name = "mozaic"
245	a2o_name = "a2o"
246	!!
247	WRITE (nout,*) 'o2amask (>0, <=0, >1) ', COUNT (o2amask > 0.0_rl), COUNT (o2amask <= 0.0_rl) &
248	& , COUNT (o2amask >= 1.0_rl)
249	i2a = 0_il
250	WHERE (RESHAPE(o2amask,(/jpai,jpaj/)) .GT. 0.0_rl )
251	i2a = 1_il
252	END WHERE
253	WRITE (nout, *) 'i2a (o2amask > 0)'
254	CALL prihin (i2a , 1_il, nout)
255	!!
256	!! Verif poids avant reduction
257	!!
258	WRITE (unit=nout,fmt=*) 'Verif poids globale'
259	za = 0.0_rl
260	WRITE (nout, *) 'Comptage xasrft ' , COUNT ( xasrft < 1.0E-10_rl )
261	WRITE (nout, *) 'Comptage o2amask ' , COUNT ( o2amask < 0.5_rl )
262	WHERE ( xasrft * o2amask /= 0.0_rl )
263	za = 1.0_rl / (xasrft * o2amask)
264	END WHERE
265	WRITE (nout, *) 'Somme za ', SUM ( za)
266	WRITE (nout, *) 'Somme xasrft ', SUM ( xasrft)
267	zsuma = DOT_PRODUCT (za, xasrft * o2amask)
268	!!
269	WRITE (nout, *) 'jma2o for coastal runoff : ', jma2o
270	!! Atm mask interpolated to oce
271	za = REAL (1_il-iamskt, KIND = rl )
272	CALL inter_a2o (za, a2omask )
273	!! 1 on Atm interpolated to oce
274	za = REAL (1_il , KIND = rl )
275	CALL inter_a2o (za, a2ofull )
276	!!
277	CALL verif ('Apres lecture')
278	!!
279	CALL bilan ("Poids ", c_case="complet", c_int_atm="Local", c_int_oce="Local")
280	CALL bilan ("Poids ", c_case="ocean" , c_int_atm="Local", c_int_oce="Local")
281	!!
282	!! Remise des poids en integre
283	DO jn = 1, jpa2o
284	DO jo = 1, jpon
285	ja = ka2o (jn, jo)
286	IF (ja > 0) THEN
287	IF ( o2amask (ja) .GT. 0.0_rl) &
288	& wa2o (jn, jo) = wa2o (jn, jo) * xosrft (jo) / ( xasrft (ja) * o2amask (ja))
289	END IF
290	END DO
291	END DO
292	!
293	CALL bilan ("Integre ", c_case="complet", c_int_atm="Integral", c_int_oce="Integral")
294	CALL bilan ("Integre ", c_case="ocean" , c_int_atm="Integral", c_int_oce="Integral")
295	CALL bilan ("Integre ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral")
296	!!
297	E_DRYRUN: IF (l_dryrun) THEN
298	WRITE (nout,*) "Remplissage bidon des champs pour tester les I/O"
299	jma2o = jma2or
300	nvo = jma2or
301	WRITE (nout,*) 'jma2o pour test : ', jma2o
302	CALL RANDOM_SEED
303	CALL RANDOM_NUMBER (wa2o) ; ka2o = NINT (wa2o*REAL(jpan,rl) )
304
305	DO jo = 1, jpon
306	IF ( iomskt (jo) == 1 ) THEN
307	wa2o (:, jo) = 0.0_rl
308	ka2o (:, jo) = 0_il
309	END IF
310	END DO
311
312	DO jo = 1, jpon
313	DO jno = 1, jpa2o
314	ja = ka2o (jno, jo)
315	IF (ja > 0) THEN
316	IF (iamskt (ja) == 1) THEN
317	wa2o (jno, jo) = 0.0_rl ; ka2o (jno, jo) = 0_il
318	END IF
319	END IF
320	END DO
321	END DO
322	ELSE
323	L_lcoast: IF (lcoast) THEN
324	!! Dataset 4 : idem 2, mais en interpolant seulement des points cotes atmosphere vers
325	!! les points cotes ocean.
326	WRITE (nout,*) 'Demarrage L_coast'
327	!!
328	WRITE (nout,*) 'Calculs points terre/oce/cote sur grille atmosphere'
329	WHERE (o2amask .LT. epsfrac ) laland = .TRUE.
330	WHERE (o2amask .GE. epsfrac .AND. o2amask .LE. (1.0_rl-epsfrac) ) lacot = .TRUE.
331	WHERE ( o2amask .GT. (1.0_rl-epsfrac) ) laoce = .TRUE.
332
333	WRITE (UNIT=nout, FMT=*) 'Points cotiers ', COUNT(lacot)
334
335	WRITE (nout,*) 'Ecriture '
336	i2a = 0
337	WHERE (RESHAPE(laland,(/jpai,jpaj/))) i2a=1
338	CALL flioputv (nid_it, "laland", i2a)
339	i2a = 0
340	WHERE (RESHAPE(laoce ,(/jpai,jpaj/))) i2a=1
341	CALL flioputv (nid_it, "laoce", i2a)
342	i2a = 0
343	WHERE (RESHAPE(lacot ,(/jpai,jpaj/))) i2a=1
344	CALL flioputv (nid_it, "lacot", i2a)
345
346	WRITE (nout,*) 'Determination des points cotes ocean '
347	ALLOCATE (ktemp2(jpoi,jpoj))
348	ktemp2 = 0_il
349	locot = .FALSE.
350	looce = .FALSE.
351	IF ( lcoast ) THEN
352	DO jo = 1, jpon
353	joi = moi (jo) ; joj = moj (jo) ! Indices 2D
354	joi_p = joi+1 ; joi_m = joi-1
355	joj_p = MIN(joj+1, jpoj) ; joj_m = MAX (joj-1, 1)
356	IF ( iomskt (jo) == 0_il .AND. iomskp (jo) == 0_il ) THEN
357	IF ( iomskt (m1o (joi_p, joj )) == 1 &
358	.OR. iomskt (m1o (joi_m, joj )) == 1 &
359	.OR. iomskt (m1o (joi , joj_p)) == 1 &
360	.OR. iomskt (m1o (joi , joj_m)) == 1 &
361	.OR. iomskt (m1o (joi_p, joj_p)) == 1 &
362	.OR. iomskt (m1o (joi_p, joj_m)) == 1 &
363	.OR. iomskt (m1o (joi_m, joj_p)) == 1 &
364	.OR. iomskt (m1o (joi_m, joj_m)) == 1 ) THEN
365	locot (jo) = .TRUE.
366	ktemp2 (joi, joj) = 1_il
367	ELSE
368	looce (jo) = .TRUE.
369	END IF
370	END IF
371	END DO
372	ELSE
373	locot (:) = ( iomskt == 0 )
374	ENDIF
375	CALL lbc (locot, noperio, 'T', jpoi, jpoj)
376	CALL lbc (looce, noperio, 'T', jpoi, jpoj)
377	WRITE (nout,*) "iomskt"
378	CALL prihin (RESHAPE(iomskt,(/jpoi,jpoj/)), 1, nout)
379	WRITE (nout,*) "Locot"
380	CALL prihin (ktemp2, 1, nout)
381	DEALLOCATE (ktemp2)
382	WRITE (nout,*) 'Verif coherence (tout doit etre nul)'
383	WRITE (nout,*) COUNT ( laland .AND. lacot)
384	WRITE (nout,*) COUNT ( lacot .AND. laoce)
385	WRITE (nout,*) COUNT ( laoce .AND. laland)
386	WRITE (nout,*) COUNT ( .NOT. ( laland .OR. lacot .OR. laoce) )
387
388	!-$$ iamsk = 0 ; WHERE (lacot) iamsk = 1
389	!-$$ WRITE(nout,*) 'Fraction lacot'
390	!-$$ CALL prihin(RESHAPE(iamsk,(/jpai,jpaj/)),1_il,nout)
391	!-$$
392	!-$$ iamsk = 0 ; WHERE (laland) iamsk = 1
393	!-$$ WRITE(nout,*) 'Points avec terre laland'
394	!-$$ CALL prihin(RESHAPE(iamsk,(/jpai,jpaj/)),1_il,nout)
395	!-$$
396	!-$$ iamsk = 0 ; WHERE (laoce) iamsk = 1
397	!-$$ WRITE(nout,*) 'Fraction laoce'
398	!-$$ CALL prihin(RESHAPE(iamsk,(/jpai,jpaj/)),1_il,nout)
399	!-$$
400	!-$$ WRITE(nout,*) 'Periodicite atm'
401	!-$$ CALL prihin(RESHAPE(iamskp,(/jpai,jpaj/)),1_il,nout)
402	!!
403	WRITE (nout,*) "iamskt"
404	CALL prihin (RESHAPE(iamskt,(/jpai,jpaj/)), 1, nout)
405	WRITE (nout,*) "Lacot"
406	ALLOCATE (ktemp2 (jpai, jpaj))
407	ktemp2 = 0_il ; WHERE (RESHAPE(lacot, (/jpai, jpaj/))) ktemp2 = 1_il
408	CALL prihin (ktemp2, 1, nout)
409	!!
410	WRITE (unit = nout, fmt = *) 'Nombre de points cotes atmosphere : ', COUNT (lacot)
411	WRITE (unit = nout, fmt = *) 'Nombre de points cotes ocean : ', COUNT (locot)
412	!!
413	CALL bilan ("Avec cote ", c_case="complet", c_int_atm="Integral", c_int_oce="Integral")
414	CALL bilan ("Avec cote ", c_case="ocean" , c_int_atm="Integral", c_int_oce="Integral")
415	CALL bilan ("Avec cote ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral")
416	!
417	! Mise a zero des poids sur les points non cotiers
418	DO jo = 1, jpon
419	IF (.NOT. locot (jo) ) THEN
420	wa2o (:, jo) = 0.0e0_rl
421	ka2o (:, jo) = 0_il
422	ELSE
423	DO jn = 1, jpa2o
424	ja = ka2o (jn, jo)
425	IF (ja > 0 ) THEN
426	IF (.NOT. lacot (ja) ) THEN
427	wa2o (jn, jo) = 0.0e0_rl
428	ka2o (jn, jo) = 0_il
429	END IF
430	END IF
431	END DO
432	END IF
433	END DO
434	CALL verif ('Apres remise a zero')
435	!!
436	!! Elimination des points redondant par periodicite
437	DO jo = 1, jpon
438	IF (iomskp (jo) .EQ. 1_il ) THEN
439	wa2o (:,jo) = 0.0e0_rl
440	ka2o (:,jo) = 0_il
441	END IF
442	END DO
443	!!
444	!! Controle
445	!!
446	CALL verif ('Apres elimination de periodicite')
447	!!
448	CALL bilan ("Apres remise a zero ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral")
449	!!
450	!! Renormalize
451	! Cas 1
452	RenormA : DO ja = 1, jpan
453	IF (.NOT. lacot (ja) ) CYCLE RenormA
454	z_d_1 = 0.0e0_rl
455	kom = 0
456	RenormO : DO jo = 1, jpon
457	IF (.NOT. locot (jo) ) CYCLE RenormO
458	DO jn = 1, jpa2o
459	IF ( ka2o (jn,jo) == ja ) THEN
460	z_d_1 = z_d_1 + wa2o (jn,jo)
461	kom = kom + 1
462	jx (kom,:) = (/ jn, jo /)
463	!!$ WRITE (unit = nout, fmt = '(6I6,16F18.6) ') ja, jo, jn, kom, jx (kom, :), &
464	!!$ wa2o (jn,jo), wa2o (jx (kom, 1), jx (kom, 2)), z_d_1, &
465	!!$ xosrft (jo), xasrft (ja)
466	END IF
467	END DO
468	ENDDO RenormO
469	IF (kom /= 0 ) THEN
470	IF (z_d_1 == 0.0_rl ) THEN
471	WRITE (nout, *) 'z_d_1 nul ', ja, xasrft (ja), lacot (ja), iamskt (ja), &
472	& wa2o (:,ja), ka2o (:,ja)
473	STOP
474	ELSE
475	DO ko = 1, kom
476	wa2o (jx (ko,1), jx (ko,2)) = wa2o (jx (ko,1), jx (ko,2)) / z_d_1
477	END DO
478	ENDIF
479	END IF
480	ENDDO RenormA
481	!!
482	!! Calcul des distances à la cote des points oceans
483	!-$$ WRITE (nout,*) 'Distance cotes'
484	!-$$ d_cote = rpi
485	!-$$ Cal_dist_1: DO jo = 1, jpon
486	!-$$ IF ( iomskt (jo) == 1 .OR. locot (jo) ) THEN
487	!-$$ d_cote (jo) = 0.0E0_rl
488	!-$$ ELSE
489	!-$$ Cal_dist_2:DO jo2 = 1, jpon
490	!-$$ IF ( locot (jo2) ) THEN
491	!-$$ d_cote (jo) = MIN (d_cote (jo), REAL (geodist(xolont(jo), xolatt(jo), xolont(jo2), xolatt(jo2)), KIND=rl))
492	!-$$ END IF
493	!-$$ END DO Cal_dist_2
494	!-$$ END IF
495	!-$$ END DO Cal_dist_1
496	!-$$ d_cote = d_cote * r_earth
497	!CALL prihre (RESHAPE(d_cote,(/jpoi,jpoj/)), 1.0E-5_rl, nout)
498	WRITE (nout,*) 'Distance maxi a la cote (ocean) : ', MAXVAL (d_cote)
499	CALL fliocrfd ("dist_cote" // TRIM(c_suffix), (/ "x", "y"/), (/jpoi, jpoj/), n1, MODE="REPLACE")
500	CALL flioputa (n1, '?', 'Comment', TRIM(c_comment) )
501	CALL fliodefv (n1, "dist_cote", (/1, 2/))
502	CALL flioputv (n1, "dist_cote", RESHAPE (d_cote, (/jpoi, jpoj/) ))
503	CALL flioclo (n1)
504	!
505	!! Controle
506	!!
507	CALL verif ('Premiere normalisation')
508	CALL gather_wei ('Premiere etape')
509	CALL verif ('Premier gather')
510	!!
511	CALL bilan ("Renorm ", c_case="complet", c_int_atm="Integral", c_int_oce="Integral")
512	CALL bilan ("Renorm ", c_case="ocean" , c_int_atm="Integral", c_int_oce="Integral")
513	CALL bilan ("Renorm ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral")
514	!!
515	!! Traitement des points atmospheres non atteint
516	!!
517	total: IF (ltotal) THEN
518	WRITE (nout,*) 'Va chercher tout les points terre (ltotal=TRUE)'
519	SeekAtm_1_1: DO ja = 1, jpan
520	!-$$ WRITE(nout,*) 'Point atm : ', ja, m2ai(ja), m2aj(ja), itarget(ja), iamskp(ja), laland(ja), lacot(ja), laoce(ja)
521	IF (itarget (ja) .GE. 1_il) CYCLE SeekAtm_1_1
522	IF (iamskp (ja) .EQ. 1_il) CYCLE SeekAtm_1_1
523	IF (.NOT. laland (ja) ) CYCLE SeekAtm_1_1
524	zzmin = REAL (cartdist(0.0_rl, 90.0_rl, 0.0_rl, -90.0_rl), KIND=rl)
525	ltrouve = .FALSE. ; jot = 0_il
526	!-$$ WRITE(nout,*) 'Point atm : ', ja, m2ai(ja), m2aj(ja)
527	SeekOce_1_1: DO jo = 1, jpon ! Cherche point ocean le plus proche
528	IF (.NOT. locot (jo) ) CYCLE SeekOce_1_1
529	IF (iomskp(jo) .EQ. 1_il) CYCLE SeekOce_1_1
530	z_d_1 = REAL (cartdist(xalont(ja), xalatt(ja), xolont(jo), xolatt(jo)), KIND=rl)
531	IF (z_d_1 .LT. zzmin) THEN
532	zzmin = z_d_1
533	ltrouve = .TRUE. ; jot = jo
534	END IF
535	ENDDO SeekOce_1_1
536	Trouve_1_1: IF (ltrouve) THEN
537	! Calcul du poids
538	z_d_1 = 1.0_rl
539	! Rangement des poids
540	IF (nvo (jot) < jpa2o ) THEN ! Reste un poids libre
541	nvo (jot) = nvo (jot) + 1
542	wa2o (nvo (jot), jot) = z_d_1
543	ka2o (nvo (jot), jot) = ja
544	!-$$ WRITE (nout,*) ' Vers point oce ', jot, m2oi(jot), m2oj(jot), nvo(jot)
545	ELSE ! Trouve12
546	WRITE(nout,*) ' Plus de place ', ja, jot, nvo (jot), z_d_1
547	STOP
548	ENDIF
549	ELSE ! Trouve
550	!! Message d'erreur
551	WRITE (unit = nout, fmt = '("Pas de voisin proche pour ", 1I5,2I4)' ) ja, m2ai(ja), m2aj(ja)
552	END IF Trouve_1_1
553	END DO SeekAtm_1_1
554	!!
555	CALL verif ('ltotal avant gather')
556	CALL bilan ("Ltotal avant gather ", c_case="terre" , c_int_atm="Integral", c_int_oce="Integral")
557	CALL bilan ("Ltotal avant gather ", c_case="terre100", c_int_atm="Integral", c_int_oce="Integral")
558	CALL bilan ("Ltotal avant gather ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral")
559	!!
560	CALL gather_wei ('apres total')
561	CALL verif ( 'Cas total', l_non_a=.TRUE.)
562	!!
563	WRITE (unit = nout, fmt = * ) 'Nombre de voisins apres completion : ', jma2o
564	!!
565	CALL bilan ("Fin ltotal ", c_case="terre" , c_int_atm="Integral", c_int_oce="Integral")
566	CALL bilan ("Fin ltotal ", c_case="terre100", c_int_atm="Integral", c_int_oce="Integral")
567	CALL bilan ("Fin ltotal ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral")
568	!!
569	END IF Total
570	!!
571	Total_dist: IF (ltotal_dist) THEN
572	!! On limite la recherche à une certaine distance de la cote.
573	WRITE (nout,*) 'Cas ltotal_dist'
574	SeekAtm_2_1: DO ja = 1, jpan
575	IF (itarget (ja) .GE. 1_il) CYCLE SeekAtm_2_1
576	IF (iamskp (ja) .EQ. 1_il) CYCLE SeekAtm_2_1
577	IF (.NOT. laland (ja) ) CYCLE SeekAtm_2_1
578	zzmin = dist_max
579	ltrouve = .FALSE. ; jot = 0_il
580	!$$$ WRITE(nout,*) 'Point atm ', ja, m2ai(ja), m2aj(ja)
581	SeekOce_2_1: DO jo = 1, jpon ! Cherche point ocean le plus proche
582	IF (.NOT. locot (jo) ) CYCLE SeekOce_2_1
583	IF (iomskp(jo) .EQ. 1_il) CYCLE SeekOce_2_1
584	z_d_1 = r_earth * REAL (geodist(xalont(ja), xalatt(ja), xolont(jo), xolatt(jo)), KIND=rl)
585	IF (z_d_1 .LT. zzmin) THEN
586	zzmin = z_d_1
587	ltrouve = .TRUE. ; jot = jo
588	END IF
589	ENDDO SeekOce_2_1
590	Trouve_2_1: IF (ltrouve) THEN
591	! Calcul du poids
592	z_d_1 = 1.0_rl
593	IF (nvo (jot) < jpa2o ) THEN ! Reste un poids libre
594	nvo (jot) = nvo (jot) + 1
595	wa2o (nvo (jot), jot) = z_d_1
596	ka2o (nvo (jot), jot) = ja
597	!$$$ WRITE (nout,*) ' Vers point oce ', jot, m2oi(jot), m2oj(jot), nvo(jot)
598	ELSE
599	WRITE(nout,*) ' Plus de place ', ja, jot, nvo (jot), z_d_1
600	STOP
601	ENDIF
602	END IF Trouve_2_1
603	END DO SeekAtm_2_1
604	!!
605	CALL verif ('Ltotal_dist avant gather')
606	CALL bilan ("Ltotal_dist avant gather ", c_case="terre" , c_int_atm="Integral", c_int_oce="Integral")
607	CALL bilan ("Ltotal_dist avant gather ", c_case="terre100", c_int_atm="Integral", c_int_oce="Integral")
608	CALL bilan ("Ltotal_dist avant gather ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral")
609	!!
610	CALL gather_wei ('apres total_dist')
611	CALL verif ( 'Cas total_dist')
612	!!
613	WRITE (unit = nout, fmt = * ) 'Nombre de voisins apres completion : ', jma2o
614	!!
615	CALL bilan ("Fin ltotal_dist ", c_case="terre" , c_int_atm="Integral", c_int_oce="Integral")
616	CALL bilan ("Fin ltotal_dist ", c_case="terre100", c_int_atm="Integral", c_int_oce="Integral")
617	CALL bilan ("Fin ltotal_dist ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral")
618	!!
619	END IF Total_dist
620	!!
621	Total_dist_2: IF (ltotal_dist_2) THEN
622	!! On limite la recherche à une certaine distance de la cote.
623	!! On etale sur les point oceans cotiers proches
624	WRITE (nout,*) '-------------------'
625	WRITE (nout,*) 'Cas ltotal_dist_2'
626	WRITE (nout,*) 'Distance recherche : ', dist_max/1.0E3_rl, 'km'
627	WRITE (nout,*) 'Distance etalement ocean : ', dist_max_voisin/1.0E3_rl, 'km'
628	!!
629	SeekAtm_3_1: DO ja = 1, jpan
630	IF (itarget (ja) .GE. 1_il) CYCLE SeekAtm_3_1
631	IF (iamskp (ja) .EQ. 1_il) CYCLE SeekAtm_3_1
632	IF (.NOT. laland (ja) ) CYCLE SeekAtm_3_1
633	zzmin = dist_max
634	ltrouve = .FALSE. ; jot = 0_il ; ntrouve = 0
635	!WRITE(nout,*) 'Point atm : ', ja, m2ai(ja), m2aj(ja)
636	SeekOce_3_1: DO jo = 1, jpon ! Cherche point ocean le plus proche
637	IF (.NOT. locot (jo) ) CYCLE SeekOce_3_1
638	IF (iomskp (jo) .EQ. 1_il) CYCLE SeekOce_3_1
639	z_d_1 = r_earth * REAL (geodist (xalont (ja), xalatt (ja), xolont (jo), xolatt (jo)), KIND=rl)
640	IF (z_d_1 .LT. zzmin) THEN
641	zzmin = z_d_1 ; ltrouve = .TRUE. ; jot = jo ; ntrouve = ntrouve+1
642	END IF
643	ENDDO SeekOce_3_1
644	WRITE(nout,*) 'Point atm ', ja, m2ai(ja), m2aj(ja), jot, ntrouve
645	Trouve_3_1: IF (ltrouve) THEN
646	! On cherche les points proches du point ocean trouve
647	z_d_1 = 0.0_rl ; jo_tab = 0 ; jo_nb = 0
648	SeekOce_3_2: DO jo2 = 1, jpon
649	IF (.NOT. locot (jo2) ) CYCLE SeekOce_3_2
650	IF (iomskp (jo2) .EQ. 1_il) CYCLE SeekOce_3_2
651	IF ( r_earth * REAL (geodist (xolont (jot), xolatt (jot), xolont (jo2), xolatt (jo2)), KIND=rl) &
652	& < dist_max_voisin) THEN
653	IF (jo_nb >= jp_nb) THEN
654	WRITE (nout,*) 'Trop de point proches ', ja, jot, jo_tab
655	STOP
656	ELSE
657	jo_nb = jo_nb + 1
658	z_d_1 = z_d_1 + xosrft (jo2)
659	jo_tab (jo_nb) = jo2
660	END IF
661	END IF
662	! Calcul du poids
663	END DO SeekOce_3_2
664	IF (jo_nb == 0) THEN
665	WRITE (nout,*) 'Pas de voisin au voisin ! ', ja, jot, m2oi(jot), m2oj(jot), xolont(jot), xolatt(jot)
666	WRITE (nout,*) 'Point atm ', ja, m2ai(ja), m2aj(ja)
667	WRITE (nout,) r_earth REAL (geodist (xolont (jot), xolatt (jot), xolont (jot), xolatt (jot)), KIND=rl)
668	STOP
669	ELSE
670	!-$$ WRITE(nout,*) 'Point atm ', ja, m2ai(ja), m2aj(ja), jot, jo_nb, jo_tab (1:jo_nb)
671	DO jn = 1, jo_nb
672	jo3 = jo_tab (jn)
673	IF (nvo (jo3) < jpa2o ) THEN ! Reste un poids libre
674	nvo (jo3) = nvo (jo3) + 1
675	wa2o (nvo (jo3), jo3) = xosrft (jo3) / z_d_1
676	ka2o (nvo (jo3), jo3) = ja
677	!$$$ WRITE (nout,*) ' Vers point oce ', jot, m2oi(jot), m2oj(jot), nvo(jot)
678	ELSE
679	WRITE(nout,*) ' Plus de place ', ja, jo3, nvo (jo3), z_d_1
680	STOP
681	ENDIF
682	END DO
683	END IF
684	!-$$ ELSE
685	!-$$ WRITE (nout,*) 'Pas de voisins proche trouve '
686	!-$$ WRITE (nout,*) 'Point atm ', ja, m2ai(ja), m2aj(ja)
687	!-$$ WRITE (nout,*) dist_max, zzmin, z_d_1
688	!-$$ !STOP
689	END IF Trouve_3_1
690	END DO SeekAtm_3_1
691	!!
692	CALL verif ('Ltotal_dist_2 avant gather')
693	CALL bilan ("Ltotal_dist_2 avant gather ", c_case="terre" , c_int_atm="Integral", c_int_oce="Integral")
694	CALL bilan ("Ltotal_dist_2 avant gather ", c_case="terre100", c_int_atm="Integral", c_int_oce="Integral")
695	CALL bilan ("Ltotal_dist_2 avant gather ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral")
696	!!
697	CALL gather_wei ('apres total_dist_2')
698	CALL verif ( 'Cas total_dist_2')
699	!!
700	WRITE (unit = nout, fmt = * ) 'Nombre de voisins apres completion : ', jma2o
701	!!
702	CALL bilan ("Fin ltotal_dist_2 ", c_case="terre" , c_int_atm="Integral", c_int_oce="Integral")
703	CALL bilan ("Fin ltotal_dist_2 ", c_case="terre100", c_int_atm="Integral", c_int_oce="Integral")
704	CALL bilan ("Fin ltotal_dist_2 ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral")
705	!!
706	END IF Total_dist_2
707	!!
708	Total_dist_3: IF (ltotal_dist_3) THEN
709	!! On limite la recherche à une certaine distance de la cote.
710	!! On etale sur les point oceans proches, cote ou pas
711	WRITE (nout,*) '-------------------'
712	WRITE (nout,*) 'Cas ltotal_dist_3'
713	WRITE (nout,*) 'Distance recherche : ', dist_max/1.0E3_rl, 'km'
714	WRITE (nout,*) 'Distance etalement ocean : ', dist_max_voisin/1.0E3_rl, 'km'
715	!!
716	SeekAtm_4_1: DO ja = 1, jpan
717	IF (itarget (ja) .GE. 1_il) CYCLE SeekAtm_4_1
718	IF (iamskp (ja) .EQ. 1_il) CYCLE SeekAtm_4_1
719	IF (.NOT. laland (ja) ) CYCLE SeekAtm_4_1
720	zzmin = dist_max
721	ltrouve = .FALSE. ; jot = 0_il
722	!-$$ WRITE(nout,*) 'Point atm : ', ja, m2ai(ja), m2aj(ja)
723	SeekOce_4_1: DO jo = 1, jpon ! Cherche point ocean le plus proche
724	IF (iomskp (jo) .EQ. 1_il .OR. iomskt (jo) .EQ. 1_il ) CYCLE SeekOce_4_1
725	z_d_1 = r_earth * REAL (geodist (xalont (ja), xalatt (ja), xolont (jo), xolatt (jo)), KIND=rl)
726	IF (z_d_1 .LT. zzmin) THEN
727	zzmin = z_d_1 ; ltrouve = .TRUE. ; jot = jo
728	END IF
729	ENDDO SeekOce_4_1
730	!-$$ WRITE(nout,*) 'Point atm ', ja, m2ai(ja), m2aj(ja), jot
731	Trouve_4_1: IF (ltrouve) THEN
732	! On cherche les points proches du point ocean trouve, le long de la cote et vers le large
733	z_d_1 = 0.0_rl ; jo_tab = 0 ; jo_nb = 0
734	SeekOce_4_2: DO jo2 = 1, jpon
735	IF (iomskp (jo2) .EQ. 1_il .OR. iomskt (jo2) .EQ. 1_il) CYCLE SeekOce_4_2
736	IF ( .NOT. l_large .AND. .NOT. locot (jo2)) CYCLE SeekOce_4_2
737	z_d_2 = r_earth * REAL (geodist (xolont (jot), xolatt (jot), xolont (jo2), xolatt (jo2)), KIND=rl)
738	IF ( ( locot (jo2) .AND. z_d_2 < dist_max_cote ) &
739	.OR. ( .NOT. locot (jo2) .AND. z_d_2 < dist_max_large ) ) THEN
740	IF (jo_nb > jp_nb) THEN
741	WRITE (nout,*) 'Trop de point proches ', ja, jot, jo_tab
742	STOP
743	ELSE
744	jo_nb = jo_nb + 1
745	z_d_1 = z_d_1 + xosrft (jo2)
746	jo_tab (jo_nb) = jo2
747	END IF
748	END IF
749	! Calcul du poids
750	END DO SeekOce_4_2
751	IF (jo_nb == 0) THEN
752	WRITE (nout,*) 'Pas de voisin au voisin ! ', ja, jot, m2oi(jot), m2oj(jot), xolont(jot), xolatt(jot)
753	WRITE (nout,*) 'Point atm ', ja, m2ai(ja), m2aj(ja)
754	WRITE (nout,) r_earth REAL (geodist (xolont (jot), xolatt (jot), xolont (jot), xolatt (jot)), KIND=rl)
755	STOP
756	ELSE
757	!-$$ WRITE(nout,*) 'Point atm ', ja, m2ai(ja), m2aj(ja), jot, jo_nb, jo_tab (1:jo_nb)
758	DO jn = 1, jo_nb
759	jo3 = jo_tab (jn)
760	IF (nvo (jo3) < jpa2o ) THEN ! Reste un poids libre
761	nvo (jo3) = nvo (jo3) + 1
762	wa2o (nvo (jo3), jo3) = xosrft (jo3) / z_d_1
763	ka2o (nvo (jo3), jo3) = ja
764	!$$$ WRITE (nout,*) ' Vers point oce ', jot, m2oi(jot), m2oj(jot), nvo(jot)
765	ELSE
766	WRITE(nout,*) ' Plus de place ', ja, jo3, nvo (jo3), z_d_1
767	STOP
768	ENDIF
769	END DO
770	END IF
771	ELSE
772	WRITE (nout,*) 'Pas de voisins proche trouve '
773	WRITE (nout,*) 'Point atm ', ja, m2ai(ja), m2aj(ja)
774	WRITE (nout,*) dist_max, zzmin, z_d_1
775	STOP
776	END IF Trouve_4_1
777	END DO SeekAtm_4_1
778	!!
779	CALL verif ('Ltotal_dist_3 avant gather')
780	CALL bilan ("Ltotal_dist_3 avant gather ", c_case="terre" , c_int_atm="Integral", c_int_oce="Integral")
781	CALL bilan ("Ltotal_dist_3 avant gather ", c_case="terre100", c_int_atm="Integral", c_int_oce="Integral")
782	CALL bilan ("Ltotal_dist_3 avant gather ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral")
783	!!
784	CALL gather_wei ('apres total_dist_3')
785	CALL verif ( 'Cas total_dist_3')
786	!!
787	WRITE (unit = nout, fmt = * ) 'Nombre de voisins apres completion : ', jma2o
788	!!
789	CALL bilan ("Fin ltotal_dist_3 ", c_case="terre" , c_int_atm="Integral", c_int_oce="Integral")
790	CALL bilan ("Fin ltotal_dist_3 ", c_case="terre100", c_int_atm="Integral", c_int_oce="Integral")
791	CALL bilan ("Fin ltotal_dist_3 ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral")
792	!!
793	END IF Total_dist_3
794	!!
795	Essai: IF (lessai) THEN
796	WRITE (nout,*) '--- Methode essai ----'
797	no = COUNT ( laland .AND. (itarget == 0_il))
798	WRITE (nout,*) 'Nombre de points atm non vises : ', no
799	WRITE (nout,*) 'jma2o, jpa2o, jma2o+no ', jma2o, jpa2o, jma2o+no
800	IF ( jma2o + no .GT. jpa2o ) THEN
801	WRITE (nout, *) 'Pas assez de voisins pour la methode essai'
802	STOP
803	END IF
804	!! Liste des points oceans
805	no = COUNT ( looce .AND. iomskp == 0_il)
806	WRITE (nout,*) 'Nombre de points oce cibles : ', no
807	ALLOCATE ( jo_n (no), STAT=ierr) ; CALL chk_allo (ierr, 'jo_n')
808	!!
809	no = 0 ; z_d_1 = 0.0_rl
810	DO jo = 1, jpon
811	IF ( .NOT. looce (jo) ) CYCLE
812	IF ( iomskp (jo) == 1 ) CYCLE
813	no = no + 1
814	jo_n (no) = jo
815	z_d_1 = z_d_1 + xosrft (jo)
816	END DO
817	no = COUNT ( looce .AND. iomskp == 0_il)
818	WRITE (nout,*) 'Surface : ', z_d_1
819	SeekAtm_5_1: DO ja = 1, jpan
820	IF (itarget (ja) .GE. 1_il) CYCLE SeekAtm_5_1
821	IF (iamskp (ja) .EQ. 1_il) CYCLE SeekAtm_5_1
822	IF (.NOT. laland (ja) ) CYCLE SeekAtm_5_1
823	DO kn = 1, no
824	jo = jo_n (kn)
825	IF (nvo (jo) < jpa2o ) THEN ! Reste un poids libre
826	nvo (jo) = nvo (jo) + 1
827	wa2o (nvo (jo), jo) = xosrft (jo) / z_d_1
828	ka2o (nvo (jo), jo) = ja
829	ELSE
830	WRITE(nout,*) ' Plus de place ', ja, jot, nvo (jot), z_d_1
831	STOP
832	ENDIF
833	END DO
834	END DO SeekAtm_5_1
835	!!
836	CALL gather_wei ('apres essai')
837	!!
838	CALL verif ( 'Cas essai')
839	CALL bilan ("Lessai avant gather ", c_case="terre" , c_int_atm="Integral", c_int_oce="Integral")
840	CALL bilan ("Lessai avant gather ", c_case="terre100", c_int_atm="Integral", c_int_oce="Integral")
841	CALL bilan ("Lessai avant gather ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral")
842	!!
843	WRITE (unit = nout, fmt = * ) 'Nombre de voisins apres completion : ', jma2o
844	END IF Essai
845	!!
846	!! Traitement des points atmosphere en bord de cote
847	!!
848	!-$$ IfNear: IF (lnear) THEN
849	!-$$ WRITE (nout,*) 'Extension de 1 point a l''interieur, vers le point ocean le plus proche (lnear=TRUE)'
850	!-$$ SeekAtm1: DO ja = 1, jpan
851	!-$$ IF (.NOT. laland_bord(ja) ) CYCLE SeekAtm1
852	!-$$ IF (iamskp(ja) .EQ. 1_il ) CYCLE SeekAtm1
853	!-$$ IF (.NOT. laland(ja) ) CYCLE SeekAtm1
854	!-$$ zzmin = REAL (geodist(0.0_rl, 90.0_rl, 0.0_rl, -90.0_rl), KIND=rl)
855	!-$$ ltrouve = .FALSE. ; jot = 0_il
856	!-$$ SeekOce1: DO jo = 1, jpon
857	!-$$ IF (.NOT. locot(jo) ) CYCLE SeekOce1
858	!-$$ IF (iomskp(jo) .EQ. 1_il ) CYCLE SeekOce1
859	!-$$ z_d_1 = REAL (geodist(xalont(ja), xalatt(ja), xolont(jo), xolatt(jo)), KIND=rl)
860	!-$$ IF (z_d_1 .LT. zzmin) THEN
861	!-$$ zzmin = z_d_1
862	!-$$ ltrouve = .TRUE. ; jot = jo
863	!-$$ END IF
864	!-$$ ENDDO SeekOce1
865	!-$$ Trouve11: IF (ltrouve) THEN
866	!-$$ WRITE (nout,*) 'Point oce proche ', jot, m2oi(jot), m2oj(jot)
867	!-$$ !! Calcul du poids
868	!-$$ z_d_1 = xasrft(ja)/xosrft(jot)
869	!-$$ !! Rangement des poids
870	!-$$ ! On cherche si ce point atm est deja vise par ce point ocean (normalement non)
871	!-$$ ltrouve2 = .FALSE.
872	!-$$ L11: DO jn = 1, jpa2o
873	!-$$ IF (ka2o(jn,jot) == ja ) THEN
874	!-$$ wa2o(jn,jot) = wa2o(jn,jot) + z_d_1
875	!-$$ WRITE (nout,*) 'Ajout au point oce ', jo, m2oi(jo), m2oj(jo)
876	!-$$ ltrouve2 = .TRUE.
877	!-$$ EXIT L11
878	!-$$ ENDIF
879	!-$$ ENDDO L11
880	!-$$ Trouve12: IF (.NOT. ltrouve2 ) THEN
881	!-$$ IF (nvo(jot) < jpa2o ) THEN ! Reste un poids libre
882	!-$$ nvo(jot) = nvo(jot) + 1
883	!-$$ wa2o(nvo(jot),jot) = z_d_1
884	!-$$ ka2o(nvo(jot),jot) = ja
885	!-$$ WRITE (nout,*) 'Vers point oce ', jot, m2oi(jot), m2oj(jot)
886	!-$$ ELSE ! Trouve2
887	!-$$ WRITE(nout,*) 'Plus de place'
888	!-$$ STOP
889	!-$$ ENDIF
890	!-$$ END IF Trouve12
891	!-$$ ELSE ! Trouve
892	!-$$ !! Message d'erreur
893	!-$$ WRITE (unit = nout, fmt = '("Pas de voisin proche pour ", 1I5,2I4)' ) ja, m2ai(ja), m2aj(ja)
894	!-$$ END IF Trouve11
895	!-$$ END DO SeekAtm1
896	!-$$ !!
897	!-$$ CALL verif ( 'Cas Near' )
898	!-$$ !!
899	!-$$ WRITE (unit = nout, fmt = * ) 'Nombre de voisins apres completion : ', jma2o
900	!-$$ !!
901	!-$$ END IF IfNear
902	!!
903	!-$$ IfNei: IF (lnei) THEN
904	!-$$ WRITE(nout,*) 'Route le run-off vers les autres voisins mouilles (lnei=TRUE)'
905	!-$$ i2a = 0_il
906	!-$$ WHERE (RESHAPE(o2amask,(/jpai,jpaj/)) .GT. 0.0_rl ) i2a = 1_il
907	!-$$ SeekAtm2: DO ja = 1_il, jpan
908	!-$$ IF (.NOT. laland_bord(ja) ) CYCLE SeekAtm2
909	!-$$ IF (iamskp(ja) .EQ. 1_il ) CYCLE SeekAtm2
910	!-$$ IF (.NOT. laland(ja) ) CYCLE SeekAtm2
911	!-$$ !! Nombre de voisins
912	!-$$ jai = m2ai(ja) ; jaj = m2aj(ja) ! Indices 2D
913	!-$$ jai_p = mai(jai+1_il) ; jai_m = mai(jai-1_il)
914	!-$$ jaj_p = maj(jaj+1_il) ; jaj_m = maj(jaj-1_il)
915	!-$$ inum = &
916	!-$$ & i2a (jai_p, jaj ) &
917	!-$$ & + i2a (jai_m, jaj ) &
918	!-$$ & + i2a (jai , jaj_p) &
919	!-$$ & + i2a (jai , jaj_m) &
920	!-$$ & + i2a (jai_p, jaj_p) &
921	!-$$ & + i2a (jai_p, jaj_m) &
922	!-$$ & + i2a (jai_m, jaj_p) &
923	!-$$ & + i2a (jai_m, jaj_m)
924	!-$$ IF (inum == 0) WRITE(nout,*) "erreur ", ja
925	!-$$ !!
926	!-$$ IF (i2a(jai_p, jaj ) == 1_il ) CALL complet_run (m1a(jai_p, jaj ), m1a(jai,jaj), inum)
927	!-$$ IF (i2a(jai_m, jaj ) == 1_il ) CALL complet_run (m1a(jai_m, jaj ), m1a(jai,jaj), inum)
928	!-$$ IF (i2a(jai , jaj_p) == 1_il ) CALL complet_run (m1a(jai , jaj_p), m1a(jai,jaj), inum)
929	!-$$ IF (i2a(jai , jaj_m) == 1_il ) CALL complet_run (m1a(jai , jaj_m), m1a(jai,jaj), inum)
930	!-$$ IF (i2a(jai_p, jaj_p) == 1_il ) CALL complet_run (m1a(jai_p, jaj_p), m1a(jai,jaj), inum)
931	!-$$ IF (i2a(jai_p, jaj_m) == 1_il ) CALL complet_run (m1a(jai_p, jaj_m), m1a(jai,jaj), inum)
932	!-$$ IF (i2a(jai_m, jaj_p) == 1_il ) CALL complet_run (m1a(jai_m, jaj_p), m1a(jai,jaj), inum)
933	!-$$ IF (i2a(jai_m, jaj_m) == 1_il ) CALL complet_run (m1a(jai_m, jaj_m), m1a(jai,jaj), inum)
934	!-$$ END DO SeekAtm2
935	!-$$ !!
936	!-$$ CALL verif ('Cas Nei')
937	!-$$ !!
938	!-$$ WRITE (unit = nout, fmt = * ) 'Nombre de voisins apres completion : ', jma2o
939	!-$$ !!
940	!-$$ END IF IfNei
941	!!
942	CALL int_wei
943	!!
944	!! Controle
945	!!
946	CALL verif (' Normalisation')
947	!!
948	CALL bilan ("Apres lcoast ", c_case="terre" , c_int_atm="Integral", c_int_oce="Local")
949	CALL bilan ("Apres lcoast ", c_case="terre100", c_int_atm="Integral", c_int_oce="Local")
950	CALL bilan ("Apres lcoast ", c_case="cote" , c_int_atm="Integral", c_int_oce="Local")
951	!!
952	WRITE (unit = nout, fmt = * ) 'Nombre de voisins : ', jma2o
953	!!
954	CALL gather_wei ('apres normalisation')
955	CALL verif ('apres gather' )
956
957	!! Atm mask interpolated to oce
958	za = REAL (1_il-iamskt, KIND = rl )
959	CALL inter_a2o (za, a2omask )
960	!! 1 on Atm interpolated to oce
961	za = REAL (1_il , KIND = rl )
962	CALL inter_a2o (za, a2ofull )
963
964	!!
965	!! Checking weights
966	WRITE (unit = nout, fmt = '("Neighbors a2o : ", 3I9 )' ) &
967	& MINVAL (nvo, nvo /=0), MAXVAL (nvo), jma2o
968	WRITE (unit = nout, fmt = '("Index a2o : ", 2I9 )' ) &
969	& MINVAL (ka2o, ka2o /= 0), MAXVAL (ka2o)
970	jind = MINLOC(wa2o, wa2o /= 0.0_rl) ; jo = jind(2) ; jn = jind(1)
971	joi = m2oi(jo) ; joj = m2oj(jo)
972	WRITE (unit = nout, fmt = '("Weights a2o mini : ", 4I6, 1E13.4, 3I6)' ) &
973	& jo, jn, joi, joj, wa2o(jn,jo), ka2o(jn,jo), m2ai(ka2o(jn,jo)), m2aj(ka2o(jn,jo))
974	jind = MAXLOC(wa2o, wa2o /= 0.0_rl) ; jo = jind(2) ; jn = jind(1)
975	joi = m2oi(jo) ; joj = m2oj(jo)
976	WRITE (unit = nout, fmt = '("Weights a2o maxi : ", 4I6, 1E13.4, 3I6)' ) &
977	& jo, jn, joi, joj, wa2o(jn,jo), ka2o(jn,jo), m2ai(ka2o(jn,jo)), m2aj(ka2o(jn,jo))
978	!
979	jind(1:1) = INT(MAXLOC (nvo),il) ; jo = jind(1)
980	WRITE (unit = nout, fmt = *) jo, m2oi(jo), m2oj(jo), ka2o(1:nvo(jo),jo), &
981	& (m2ai(ka2o(jn,jo)), m2aj(ka2o(jn,jo)), jn = 1, nvo(jo))
982	!!
983	clweight = "WEIGHTS4" ; cladress = "ADRESSE4"
984	WRITE (unit = nout, fmt = *) cladress, " ", clweight, " ", jpa2o, jma2o
985	!!
986	!! Controle
987	!!
988	WRITE (nout, *) 'Verif avant ecriture'
989	WRITE (nout, *) 'wa2o ', MINVAL (wa2o), MAXVAL (wa2o), COUNT (wa2o > 0.0_rl)
990	WRITE (nout, *) 'ka2o ', MINVAL (ka2o), MAXVAL (ka2o), COUNT (ka2o > 0 )
991	WRITE (nout, *) 'xasrft ', MINVAL (xasrft), MAXVAL (xasrft)
992	WRITE (nout, *) 'o2amask ', MINVAL (o2amask), MAXVAL (o2amask), SUM (o2amask)
993	WRITE (nout, *) 'xosrft ', MINVAL (xosrft), MAXVAL (xosrft)
994	WRITE (nout, *) 'iomskt ', MINVAL (iomskt), MAXVAL (iomskt), COUNT (iomskt == 0)
995	WRITE (nout, *) 'iomskp ', MINVAL (iomskp), MAXVAL (iomskp), COUNT (iomskp == 0)
996	WRITE (nout, *) 'iamskt ', MINVAL (iamskt), MAXVAL (iamskt), COUNT (iamskt == 0)
997	WRITE (nout, *) 'locot ', COUNT (locot)
998	WRITE (nout, *) 'lacot ', COUNT (lacot)
999	!!
1000	ENDIF L_LCOAST
1001	!!
1002	END IF E_DRYRUN
1003	!!
1004	!!
1005	WRITE (nout, *) ' '
1006	WRITE (nout, *) ' '
1007	WRITE (nout, *) 'Eventuels traitement specifiques'
1008	WRITE (nout, *) ' '
1009	WRITE (nout, *) ' '
1010
1011
1012	!! Extension ??
1013
1014	!! Traitements specifiques, tres specifiques ... !!
1015	IF ( TRIM(c_suffix) == "_runoffNord09" ) THEN
1016	WRITE (nout, *) 'Traitement specifique _runoffNord09'
1017	DO jo = 1, jpon
1018	IF ( tmask_atl(jo) == 1 .AND. xolatt(jo) > 40.0_rl ) wa2o(:,jo) = wa2o(:,jo) * 0.9
1019	END DO
1020	END IF
1021	!
1022	IF ( TRIM(c_suffix) == "_runoffNord07" ) THEN
1023	WRITE (nout, *) 'Traitement specifique _runoffNord07'
1024	DO jo = 1, jpon
1025	IF ( tmask_atl(jo) == 1 .AND. xolatt(jo) > 40.0_rl ) wa2o(:,jo) = wa2o(:,jo) * 0.7
1026	END DO
1027	END IF
1028
1029	IF ( TRIM(c_suffix) == "_runoffNord00" ) THEN
1030	WRITE (nout, *) 'Traitement specifique _runoffNord00'
1031	DO jo = 1, jpon
1032	IF ( tmask_atl(jo) == 1 .AND. xolatt(jo) > 40.0_rl ) wa2o(:,jo) = wa2o(:,jo) * 0.0
1033	END DO
1034	END IF
1035	!!
1036	!!
1037	!! Verif nombre de voisins
1038	!!
1039	IF ( jma2or /= jma2o ) THEN
1040	WRITE (nout, *) 'Erreur nombre de voisins dans run.def : '
1041	WRITE (nout, *) ' jma2o calcule / lu : ', jma2o, jma2or
1042	!STOP
1043	ENDIF
1044	!!
1045	!! Ecriture des poids
1046	!!
1047	WRITE (unit = c_r, fmt = '(1i1)' ) rk_out
1048	!
1049	IF ( l_wei_i4) THEN
1050	WRITE (unit = c_i, fmt = '(1i1)' ) i_4
1051	cfname4 = TRIM(mo_name) // ".wa2o.runoff.i" // TRIM(c_i) // ".r" // TRIM(c_r) // TRIM(c_suffix)
1052	OPEN (unit = nwei4a2o, file = TRIM(cfname4), form = 'unformatted', status = 'replace',&
1053	& action = 'write')
1054	WRITE (nout, *) 'Poids i_4 a -> o: ', TRIM(cfname4)
1055	END IF
1056	!
1057	IF (l_wei_i8) THEN
1058	WRITE (unit = c_i, fmt = '(1i1)' ) i_8
1059	cfname8 = TRIM(mo_name) // ".wa2o.runoff.i" // TRIM(c_i) // ".r" // TRIM(c_r) // TRIM(c_suffix)
1060	OPEN (unit = nwei8a2o, file = TRIM(cfname8), form = 'unformatted', status = 'replace',&
1061	& action = 'write')
1062	WRITE (nout, *) 'Poids i_8 a -> o : ', TRIM(cfname8)
1063	END IF
1064	!
1065	cldiag_a2o = TRIM(a2o_name) // '.runcoa.diag'
1066	clw_a2o = TRIM(mo_name) // '.wa2o.runcoa'
1067	clw_a2o_mct= TRIM(mo_name) // '.wa2o_mct.runcoa'
1068	!
1069	CALL wri_weights_a2o (cldiag_a2o, clw_a2o, clw_a2o_mct, "4", l_fulldiag, lo_src_grid_frac=.FALSE., lo_dst_grid_frac=.FALSE., &
1070	& co_omsk=TRIM(cotes_omsk), co_amsk=TRIM(cotes_amsk) )
1071	!
1072	IF (lriv) THEN
1073	!! Dataset 3 : cas particulier du run off des rivieres
1074	!!
1075	ka2o = 0 ; wa2o = 0.0e0_rl ; nvo = 0 ; nva = 0 ; jma2o = 1
1076	OPEN (unit = nriv, file = '/home/clim/ocean/CPL/DATA/riviere.dat', status = 'old', action = 'read', &
1077	& form = 'formatted', position = 'rewind' )
1078	jovoid = 1
1079	DO jr = 1, jpr
1080	READ (unit = nriv, fmt = '(1A1,2I3,1I5,1A30)' ) clriv, joi, joj, nlmd, clname
1081	clname = ADJUSTL (TRIM (ADJUSTL (clname)))
1082	jai = m2ai (nlmd) ; jaj = m2aj (nlmd)
1083	jaj = jpaj - jaj + 1
1084	ja = m1a(jai,jaj)
1085	IF (joi == 99 .AND. joj == 99 ) THEN
1086	!! On ramene sur la ligne j=1 de l'ocean pour faire un bilan plus tard
1087	jovoid = jovoid + 1
1088	WRITE (unit = nout, fmt = &
1089	& '(1A, " (", 1A30, ") ", 1I3, " : ", 4I5, " : ", 2F6.1 )' ) &
1090	& clriv, clname, jr, nlmd, ja, jai, jaj, xalont (ja), xalatt (ja)
1091	joi = jovoid ; joj = 1
1092	jo = m1o (joi, joj)
1093	nva (ja) = nva (ja) + 1
1094	nvo (jo) = nvo (jo) + 1
1095	ka2o (nvo (jo), jo) = ja
1096	wa2o (nvo (jo), jo) = 1.0e0_rl
1097	ELSE
1098	jo = m1o (joi, joj)
1099	nva (ja) = nva (ja) + 1
1100	nvo (jo) = nvo (jo) + 1
1101	ka2o (nvo (jo), jo) = ja
1102	IF (lint_atm .AND. lint_oce ) THEN
1103	!! Atm envoie : kg.s-1, Oce reçoit : kg.s-1
1104	wa2o(nvo(jo),jo) = 1.0e0_rl
1105	ELSE
1106	!! Atm envoie kg/m ^2/s. Oce reçoit kg/m^2/s
1107	wa2o (nvo (jo), jo) = xasrft (ja) / xosrft (jo)
1108	END IF
1109	WRITE (unit = nout, fmt = &
1110	& '(1A, " (", 1A30, ") ", 1I3, " : ", 4I5, " -> ", 2I4, 1I2, " : ", 1E12.4, " ", 2F6.1, " -> ", 2F6.1)' ) &
1111	& clriv, clname, jr, nlmd, ja, jai, jaj, joi, joj, nvo(jo), wa2o (nvo(jo),jo), xalont(ja), xalatt(ja), &
1112	& clo_lon(xolont(jo),xalont(ja)), xolatt (jo)
1113	END IF
1114	END DO
1115	CLOSE (unit = nriv)
1116	jma2o = MAXVAL (nvo)
1117	WRITE (unit = nout, fmt = '("Voisins atm : ", 2I9 )' ) &
1118	& MINVAL(nva), MAXVAL(nva)
1119	WRITE (unit = nout, fmt = '("Voisins ocean : ", 2I9 )' ) &
1120	& MINVAL(nvo), MAXVAL(nvo)
1121	WRITE (unit = nout, fmt = '("Point bidons : ", 2I9 )' ) jovoid
1122	!!
1123	END IF
1124
1125	!!
1126	!! Verif nombre de voisins
1127	!!
1128	IF ( jma2or /= jma2o ) THEN
1129	WRITE (nout, *) 'Erreur nombre de voisins dans run.def : '
1130	WRITE (nout, *) ' jma2o calcule / lu : ', jma2o, jma2or
1131	!STOP
1132	ENDIF
1133	!!
1134	!! Ecriture des poids
1135	!!
1136	!! Ouverture des fichiers poids
1137	!!
1138	WRITE (unit = c_r, fmt = '(1i1)' ) rk_out
1139	!
1140	IF (l_wei_i4) THEN
1141	WRITE (unit = c_i, fmt = '(1i1)' ) i_4
1142	cfname4 = TRIM(mo_name) // ".wa2o.rivflu" // TRIM(c_suffix) // ".i" // TRIM(c_i) // ".r" // TRIM(c_r)
1143	OPEN (unit = nwei4a2o, file = TRIM(cfname4), form = 'unformatted', status = 'replace',&
1144	& action = 'write')
1145	WRITE (nout, *) 'Poids i4 a -> o: ', TRIM(cfname4)
1146	END IF
1147	!
1148	IF (l_wei_i4) THEN
1149	WRITE (unit = c_i, fmt = '(1i1)' ) i_8
1150	cfname8 = TRIM (mo_name) // ".wa2o.rivflu" // TRIM(c_suffix) // ".i" // TRIM(c_i) // ".r" // TRIM(c_r)
1151	OPEN (unit = nwei8a2o, file = TRIM(cfname8), form = 'unformatted', status = 'replace',&
1152	& action = 'write')
1153	WRITE (nout, *) 'Poids i8 a -> o : ', TRIM(cfname8)
1154	END IF
1155
1156	!
1157	cldiag_a2o = TRIM (a2o_name) // '.rivflu.diag'
1158	clw_a2o = TRIM (mo_name) // '.wa2o.rivflu'
1159	clw_a2o_mct = "rmp_"//TRIM(camod_t)//"_to_"//TRIM(comod_t)//"_MOSAIC_rivflu"
1160	!
1161	CALL wri_weights_a2o (cldiag_a2o, clw_a2o, clw_a2o_mct, "3", l_fulldiag, lo_src_grid_frac=.FALSE., lo_dst_grid_frac=.FALSE., &
1162	& co_omsk='noperio', co_amsk='full' )
1163
1164	!!
1165	IF (l_wei_i4) CLOSE (unit = nwei8a2o)
1166	IF (l_wei_i8) CLOSE (unit = nwei4a2o)
1167	!!
1168	WRITE (nout,*) 'Appel flioclo'
1169	CALL flioclo
1170	!!
1171	STOP
1172	!!
1173	CONTAINS
1174	!!
1175	SUBROUTINE diag_itarget (clmess, lwrite)
1176	CHARACTER (len=*), INTENT (in) :: clmess
1177	LOGICAL, INTENT (in) :: lwrite
1178	INTEGER, SAVE :: nappel = 0
1179	CHARACTER (len=2) :: cc
1180	!!
1181	nappel = nappel + 1
1182	!!
1183	itarget = 0_il ; wtarget = 0.0_rl
1184	DO jo = 1, jpon
1185	DO jn = 1, jpa2o
1186	ja = ka2o (jn, jo)
1187	IF (ja /= 0 ) THEN
1188	itarget (ja) = itarget (ja) + 1_il
1189	wtarget (ja) = wtarget (ja) + wa2o (jn, jo)
1190	ENDIF
1191	END DO
1192	END DO
1193	!!
1194	!CALL lbc (itarget, naperio, 'T', jpai, jpaj)
1195	!CALL lbc (wtarget, naperio, 'T', jpai, jpaj)
1196	!!
1197	IF (lwrite) THEN
1198	!$$$ WRITE(nout,*) 'itarget, etape '// clmess
1199	!$$$ CALL prihin (RESHAPE(itarget,(/jpai,jpaj/)),1_il,nout)
1200	WRITE (unit = cc, fmt = '(1I2.2)' ) nappel
1201	CALL fliodefv (nid_it, 'itarget' // cc, (/1, 2/), v_t=flio_i4)
1202	CALL fliodefv (nid_it, 'ntarget' // cc, (/1, 2/), v_t=flio_i4)
1203	CALL fliodefv (nid_it, 'wtarget' // cc, (/1, 2/), v_t=flio_r4)
1204	CALL fliodefv (nid_it, 'non_target' // cc, (/1, 2/), v_t=flio_i4)
1205
1206
1207	CALL flioputa (nid_it, 'itarget' // cc, 'long_name ', clmess // ' itarget' // cc)
1208	CALL flioputv (nid_it, 'itarget' // cc, MIN (1, RESHAPE (itarget, (/jpai, jpaj/))) )
1209
1210	CALL flioputa (nid_it, 'ntarget' // cc, 'long_name ', clmess // ' ntarget' // cc)
1211	CALL flioputv (nid_it, 'ntarget' // cc, RESHAPE (itarget, (/jpai, jpaj/)) )
1212
1213	CALL flioputa (nid_it, 'wtarget' // cc, 'long_name ', clmess // ' wtarget' // cc)
1214	CALL flioputv (nid_it, 'wtarget' // cc, RESHAPE (wtarget, (/jpai, jpaj/)) )
1215
1216	wtarget = 0
1217	WHERE (itarget .EQ. 0 .AND. (laland .OR. lacot)) wtarget = 1.0_rl
1218	CALL flioputa (nid_it, 'non_target' // cc, 'long_name ', clmess // ' non_target' // cc)
1219	CALL flioputv (nid_it, 'non_target' //cc, RESHAPE (NINT (wtarget), (/jpai, jpaj/)) )
1220	END IF
1221	!!
1222	END SUBROUTINE diag_itarget
1223	!!
1224	SUBROUTINE int_wei
1225	!!
1226	IMPLICIT NONE
1227	WRITE(nout,*) 'Debut normalisation'
1228	!!
1229	IF (.NOT. lint_atm) THEN
1230	WRITE(nout,*) 'Multiplication par les surfaces atm'
1231	DO jo = 1, jpon
1232	DO jn = 1, jma2o
1233	ja = ka2o (jn, jo)
1234	IF (ja > 0 ) wa2o (jn,jo) = wa2o (jn,jo) * xasrft (ja)
1235	END DO
1236	ENDDO
1237	ENDIF
1238	!!
1239	IF (.NOT. lint_oce) THEN
1240	WRITE (nout,*) 'Division par les surfaces oce'
1241	DO jo = 1, jpon
1242	DO jn = 1, jma2o
1243	wa2o (jn,jo) = wa2o (jn,jo) / xosrft (jo)
1244	END DO
1245	ENDDO
1246	ENDIF
1247	END SUBROUTINE int_wei
1248	!!
1249	SUBROUTINE complet_run (ka1, ka2, kn)
1250	INTEGER (kind=il), INTENT(in) :: ka1, ka2 ! Index of atm boxes
1251	INTEGER (kind=il), INTENT(in) :: kn ! Number of neighbors
1252	INTEGER (kind=il) :: jt
1253	!!
1254	DO jo = 1, jpon
1255	IF (iomskt(jo) == 1 ) CYCLE
1256	IF (iomskp(jo) .EQ. 1_il ) CYCLE
1257	DO jn = 1, jpa2o
1258	IF (ka2o(jn,jo) == ka1 ) THEN
1259	nvo(jo) = nvo(jo)+1
1260	IF (nvo(jo) > jpa2o ) THEN
1261	WRITE(nout,*) 'Pas assez de voisins pour la completion '
1262	WRITE(nout,fmt='(1I7,2I4,2F10.1)' ) ka1, m2ai(ka1), m2aj(ka1), &
1263	& xalont(ka1), xalatt(ka1)
1264	WRITE(nout,fmt='(1I7,2I4,2F10.1)' ) ka2, m2ai(ka2), m2aj(ka2), &
1265	& xalont(ka2), xalatt(ka2)
1266	WRITE(nout,fmt='(1I7,2I4,2F10.1,2I8)') jo, m2oi(jo), m2oj(jo), &
1267	& xolont(jo), xolatt(jo), jn, nvo(jo)
1268	DO jt = 1, jpa2o
1269	WRITE(nout,fmt=*) jt, ka2o(jt,jo), xalont(ka2o(jt,jo)), xalatt(ka2o(jt,jo)), &
1270	& lacot(ka2o(jt,jo))
1271	ENDDO
1272	END IF
1273	ka2o (nvo (jo), jo) = ka2
1274	wa2o (nvo (jo), jo) = wa2o(jn,jo) * xasrft(ka2) / xasrft(ka1) / REAL(kn,KIND=rl)
1275	END IF
1276	END DO
1277	ENDDO
1278	END SUBROUTINE complet_run
1279	!!
1280	SUBROUTINE gather_wei (clmess)
1281	CHARACTER (len=*), INTENT (in) :: clmess
1282	INTEGER (kind=il) :: jn1, jn2, jo2
1283	!! Eliminate redundancy
1284	DO jo = 1, jpon
1285	DO jn1 = 1, jpa2o-1
1286	IF (ka2o (jn1,jo) == 0 ) CYCLE
1287	DO jn2 = jn1+1, jpa2o
1288	IF (ka2o (jn1, jo) == ka2o (jn2, jo) ) THEN
1289	WRITE(nout,*) 'Redondance ', jo, jn1, jn2, ka2o(jn1,jo), ka2o (jn2,jo), wa2o(jn1,jo), wa2o(jn2,jo)
1290	wa2o (jn1, jo) = wa2o (jn1, jo) + wa2o (jn2, jo)
1291	ka2o (jn2, jo) = 0
1292	END IF
1293	END DO
1294	END DO
1295	END DO
1296	!!
1297	!! Seek for strangeness
1298	DO jo = 1, jpon
1299	DO jn = 1, jpa2o
1300	ja = ka2o (jn, jo)
1301	IF (ja == 0) CYCLE
1302	IF (wa2o (jn, jo) == 0.0E0 ) THEN
1303	WRITE (nout, *) 'Null weight : ', jo, jn, ja, locot (jo), lacot (ja)
1304	END IF
1305	END DO
1306	END DO
1307	!!
1308	!! Periodicity
1309	!DO jn = 1, jpa2o
1310	! CALL lbc (wa2o (jn,:), noperio, 'T', jpoi, jpoj)
1311	! CALL lbc (ka2o (jn,:), noperio, 'T', jpoi, jpoj)
1312	!END DO
1313	CALL diag_itarget (' gather/perio - ' // clmess, lwrite=.TRUE.)
1314	!!
1315	!! Gather weights
1316	DO jo = 1, jpon
1317	ktemp = ka2o (:,jo) ; wtemp = wa2o (:,jo)
1318	ka2o (:,jo) = 0_il ; wa2o (:,jo) = 0.0_rl
1319	kn = 1_il
1320	DO jo2 = 1, jpa2o
1321	IF (ktemp (jo2) /= 0_il ) THEN
1322	ka2o (kn, jo) = ktemp (jo2)
1323	wa2o (kn, jo) = wtemp (jo2)
1324	kn = kn + 1_il
1325	END IF
1326	END DO
1327	END DO
1328	nvo = COUNT (ka2o /= 0_il, DIM = 1)
1329	jma2o = MAXVAL (nvo)
1330	!
1331	!!
1332	!! Periodicity
1333	!DO jn = 1, jpa2o
1334	! CALL lbc (wa2o(jn,:), noperio, 'T', jpoi, jpoj)
1335	! CALL lbc (ka2o(jn,:), noperio, 'T', jpoi, jpoj)
1336	!END DO
1337	!!
1338	nvo = COUNT (ka2o /= 0_il, DIM = 1)
1339	jma2o = MAXVAL (nvo)
1340	!!
1341	END SUBROUTINE gather_wei
1342	!!
1343	SUBROUTINE verif (clmess, l_non_a)
1344	!!
1345	CHARACTER (len=*), INTENT (in) :: clmess
1346	INTEGER :: ja
1347	LOGICAL, INTENT (in), OPTIONAL :: l_non_a
1348	!!
1349	nvo = COUNT (ka2o /= 0_il, DIM = 1)
1350	jma2o = MAXVAL (nvo)
1351	CALL diag_itarget (clmess, lwrite=.TRUE.)
1352	!!
1353	WRITE (nout, *) ' --- Verifs -- : ', TRIM (clmess)
1354	WRITE (nout, *) 'wa2o ', MINVAL (wa2o), MAXVAL (wa2o), COUNT (wa2o > 0.0_rl)
1355	WRITE (nout, *) 'ka2o ', MINVAL (ka2o), MAXVAL (ka2o), COUNT (ka2o > 0 )
1356	WRITE (nout, *) 'xasrft ', MINVAL (xasrft), MAXVAL (xasrft)
1357	WRITE (nout, *) 'o2amask ', MINVAL (o2amask), MAXVAL (o2amask), SUM (o2amask)
1358	WRITE (nout, *) 'xosrft ', MINVAL (xosrft), MAXVAL (xosrft)
1359	WRITE (nout, *) 'iomskt ', MINVAL (iomskt), MAXVAL (iomskt), COUNT (iomskt == 0)
1360	WRITE (nout, *) 'iomskp ', MINVAL (iomskp), MAXVAL (iomskp), COUNT (iomskp == 0)
1361	WRITE (nout, *) 'iamskt ', MINVAL (iamskt), MAXVAL (iamskt), COUNT (iamskt == 0)
1362	WRITE (nout, *) 'locot ', COUNT (locot)
1363	WRITE (nout, *) 'lacot ', COUNT (lacot)
1364	WRITE (nout, *) 'Points non attribues : ', COUNT (itarget .EQ. 0 .AND. (laland .OR. lacot))
1365
1366	IF (PRESENT (l_non_a)) THEN
1367	DO ja = 1, jpan
1368	IF (itarget (ja) == 0 .AND. iamskp(ja) == 0 .AND. (laland (ja) .OR. lacot(ja) )) THEN
1369	!-$$ WRITE (UNIT=nout,FMT='("Point atm non attribues : ", 3I7, 2F7.1, 1I4)' ) ja, m2ai(ja), m2aj (ja), xalont(ja), xalatt(ja), itarget(ja)
1370	WRITE (nout,*) ja, m2ai(ja), m2aj (ja), xalont(ja), xalatt(ja), itarget(ja), laland(ja), lacot(ja), laoce(ja)
1371	ENDIF
1372	END DO
1373	END IF
1374	END SUBROUTINE verif
1375	!!
1376	!!
1377	SUBROUTINE noroute
1378	!! Determine les points atmospheres non routes
1379	lanoroute = .TRUE.
1380
1381	DO jn = 1, jpa2o
1382	DO jo = 1, jpon
1383	ja = ka2o ( jo, jn)
1384	IF ( ja /= 0 .AND. wa2o (jo,jn) > 0.0_rl ) THEN
1385	lanoroute = .FALSE.
1386	END IF
1387	END DO
1388	END DO
1389
1390	END SUBROUTINE noroute
1391	!!
1392	SUBROUTINE bilan (clmess, c_case, c_int_atm, c_int_oce)
1393	!!
1394	CHARACTER (len=*), INTENT (in) :: clmess
1395	CHARACTER (len=*), INTENT (in) :: c_case ! Type de champ d'entree
1396	CHARACTER (len=*), INTENT (in) :: c_int_atm ! Type d'integrale sur l'atm
1397	CHARACTER (len=*), INTENT (in) :: c_int_oce ! Type d'integrale sur l'ocean
1398	!!
1399	WRITE (nout,*) "Bilan :" // TRIM (clmess) // ", cas :", TRIM (c_case), ", c_int_atm :", TRIM (c_int_atm),&
1400	& ", c_int_oce :", TRIM(c_int_oce), "."
1401	!!
1402	za = 0.0_rl
1403	!!
1404	SELECT CASE (TRIM (c_case))
1405	CASE ('complet')
1406	za = 1.0_rl * REAL (1-iamskp, KIND=rl)
1407	SELECT CASE (TRIM (c_int_atm))
1408	CASE ('Integral')
1409	za = za * xasrft
1410	zsuma = SUM (za)
1411	CASE ('Local')
1412	zsuma = DOT_PRODUCT (za, xasrft)
1413	CASE default
1414	WRITE (nout,*) 'Erreur type de bilan atm' ; STOP
1415	END SELECT
1416	!!
1417	CASE ('terre')
1418	WHERE (laland .OR. lacot)
1419	za = 1.0_rl * REAL (1-iamskp, KIND=rl)
1420	END WHERE
1421	SELECT CASE (TRIM (c_int_atm))
1422	CASE ('Integral')
1423	za = za * xasrft
1424	zsuma = SUM (za)
1425	CASE ('Local')
1426	zsuma = DOT_PRODUCT (za, xasrft*(1.0-o2amask))
1427	CASE default
1428	WRITE (nout,*) 'Erreur type de bilan atm' ; STOP
1429	END SELECT
1430	!!
1431	CASE ('terre100')
1432	WHERE (laland)
1433	za = 1.0_rl * REAL (1-iamskp, KIND=rl)
1434	END WHERE
1435	SELECT CASE (TRIM (c_int_atm))
1436	CASE ('Integral')
1437	za = za * xasrft
1438	zsuma = SUM (za)
1439	CASE ('Local')
1440	zsuma = DOT_PRODUCT (za, xasrft)
1441	CASE default
1442	WRITE (nout,*) 'Erreur type de bilan atm' ; STOP
1443	END SELECT
1444	!!
1445	CASE ('ocean')
1446	WHERE (lacot)
1447	za = 1.0_rl * REAL (1-iamskp, KIND=rl)
1448	END WHERE
1449	SELECT CASE (TRIM (c_int_atm))
1450	CASE ('Integral')
1451	za = za * xasrft
1452	zsuma = SUM (za)
1453	CASE ('Local')
1454	zsuma = DOT_PRODUCT (za, xasrft*o2amask)
1455	CASE default
1456	WRITE (nout,*) TRIM(clmess), ' : Erreur type de bilan atm' ; STOP
1457	END SELECT
1458	!!
1459	CASE ('cote')
1460	WHERE (lacot)
1461	za = 1.0_rl * REAL (1-iamskp, KIND=rl)
1462	END WHERE
1463	SELECT CASE (TRIM (c_int_atm))
1464	CASE ('Integral')
1465	za = za * xasrft
1466	zsuma = SUM (za)
1467	CASE ('Local')
1468	zsuma = DOT_PRODUCT (za, xasrft*o2amask)
1469	CASE default
1470	WRITE (nout,*) 'Erreur type de bilan atm' ; STOP
1471	END SELECT
1472	!!
1473	CASE default
1474	WRITE (nout,*)'Erreur choix de cas' ; STOP
1475	STOP
1476	END SELECT
1477	!!
1478	CALL inter_a2o (za, zo)
1479	!!
1480	SELECT CASE (TRIM (c_int_oce))
1481	CASE ('Integral')
1482	zsumo = DOT_PRODUCT (zo, REAL(1-iomskp, KIND=rl))
1483	CASE ('Local')
1484	zsumo = DOT_PRODUCT (zo*xosrft,REAL(1-iomskp, KIND=rl))
1485	CASE default
1486	WRITE (nout,*) TRIM (clmess), ' : Erreur type de bilan ' ; STOP
1487	END SELECT
1488	!!
1489	WRITE (nout, fmt='(" bilan atm: ", 1E15.6, ", bilan oce: ", 1E15.6)' ) zsuma, zsumo
1490	!!
1491	RETURN
1492	!!
1493	END SUBROUTINE bilan
1494	!!
1495	END PROGRAM cote
1496

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

Download in other formats:

Original Format