Context Navigation

source: TOOLS/MOZAIC/src/MOZAIC/cotes_etal.f90 @ 3326

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

Line
1	! -- Mode: f90 --
2	PROGRAM cotes_etal
3	USE declare
4	USE modeles
5	USE mod_lecdata
6	USE formula
7	USE math
8	USE mod_cotes_a
9	USE mod_rectifie
10	USE mod_lbc
11	USE mod_norma
12	USE bords
13	USE mod_inter
14	USE mod_prih
15	USE mod_wri_wei
16	USE mod_inipar
17	USE fliocom
18	USE getincom
19	USE errioipsl
20	!!!
21	IMPLICIT NONE
22	!!!
23	!!! A partir d'un fichier de poids, complete avec les rivieres et
24	!!! les poids pour le run off
25	!! On met a zero les poids sur les points non cotiers, et on renormalise
26	!!
27	INTEGER (kind=il), PARAMETER :: jpr = 52 ! Nombre de rivieres
28	REAL (kind=rl) :: zsuma, zsumo
29	REAL (kind=rl), DIMENSION (:), ALLOCATABLE :: za
30	REAL (kind=rl), DIMENSION (:), ALLOCATABLE :: zo
31	!!
32	CHARACTER (LEN = 1) :: clriv
33	CHARACTER (LEN = 30) :: clname
34	INTEGER (kind=il) :: ja, jo, jo2, jo3, jn, jno, joi, joj, jai, jaj, jr, kn ! indices de boucle
35	INTEGER (kind=il) :: jai_p, jai_m, jaj_p, jaj_m, n1
36	INTEGER (kind=il) :: joi_p, joi_m, joj_p, joj_m
37	REAL (kind=rl) :: z_d_1, z_d_2, zzmin
38	INTEGER (kind=il) :: nriv = 32, nlmd = 35
39	INTEGER (kind=il), DIMENSION (:, :), ALLOCATABLE :: jx ! Point trouve sous une maille atm
40	INTEGER (kind=il) :: kom, ko, niter, nc
41	LOGICAL, DIMENSION (:), ALLOCATABLE :: la_temp, lo_temp
42	REAL (kind=rl), DIMENSION (:), ALLOCATABLE :: wtemp
43	INTEGER (kind=il), DIMENSION (:), ALLOCATABLE :: ktemp
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, idest
52	REAL (kind=rl), DIMENSION (:), ALLOCATABLE :: wtarget, wdest
53	INTEGER (kind=il), DIMENSION (:,:), ALLOCATABLE :: i2a, i2o
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	INTEGER :: nid_it, nid_id
69	!!
70	CHARACTER (len = 50) :: cldiag_a2o, clw_a2o, clw_a2o_cdf, clw_a2o_mct
71	CHARACTER (len = 25) :: mo_name, a2o_name
72	CHARACTER (len = 1) :: c_i, c_r
73
74	!! Lecture de la ligne de comande
75	narg = iargc()
76	IF ( narg > 0 ) THEN
77	CALL getarg ( 1, clarg)
78	SELECT CASE (TRIM(clarg))
79	CASE Default
80	IF ( INDEX ( TRIM(clarg), '.def' ) /= 0 ) THEN
81	WRITE (nout,*) 'Lecture des parametres dans ', TRIM(clarg)
82	CALL getin_name ( TRIM (clarg) )
83	ELSE
84	WRITE (nout,*) 'Lecture des parametres dans run.def'
85	ENDIF
86	END SELECT
87	END IF
88	!!
89	CALL inipar
90	CALL alloc_modeles
91
92	!!
93	ALLOCATE (za (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'za', lreset=.TRUE., crout='cotes')
94	ALLOCATE (zo (jpon) , STAT=ierr) ; CALL chk_allo (ierr, 'zo')
95	ALLOCATE (jx (jpa2o,2) , STAT=ierr) ; CALL chk_allo (ierr, 'jx')
96	!
97	ALLOCATE (lacot (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'lacot') ; lacot = .FALSE.
98	ALLOCATE (laland (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'laland') ; laland = .FALSE.
99	ALLOCATE (laoce (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'laoce') ; laoce = .FALSE.
100	!
101	ALLOCATE (locot (jpon) , STAT=ierr) ; CALL chk_allo (ierr, 'locot') ; locot = .FALSE.
102	ALLOCATE (looce (jpon) , STAT=ierr) ; CALL chk_allo (ierr, 'looce') ; looce = .FALSE.
103	ALLOCATE (lo_nearcoast(jpon) , STAT=ierr) ; CALL chk_allo (ierr, 'lonearcoast') ; looce = .FALSE.
104	ALLOCATE (la_nearcoast(jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'lanearcoast') ; looce = .FALSE.
105	!
106	ALLOCATE (ktemp (jpa2o) , STAT=ierr) ; CALL chk_allo (ierr, 'ktemp')
107	ALLOCATE (wtemp (jpa2o) , STAT=ierr) ; CALL chk_allo (ierr, 'wtemp')
108	ALLOCATE (dist_coast_o (jpon) , STAT=ierr) ; CALL chk_allo (ierr, 'dist_coast_o')
109	ALLOCATE (dist_coast_a (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'dist_coast_a')
110	!
111	ALLOCATE (itarget (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'itarget')
112	ALLOCATE (wtarget (jpan) , STAT=ierr) ; CALL chk_allo (ierr, 'wtarget')
113	ALLOCATE (i2a (jpai, jpaj), STAT=ierr) ; CALL chk_allo (ierr, 'i2a')
114	!
115	ALLOCATE (idest (jpon) , STAT=ierr) ; CALL chk_allo (ierr, 'idest')
116	ALLOCATE (wdest (jpon) , STAT=ierr) ; CALL chk_allo (ierr, 'wdest')
117	ALLOCATE (i2o (jpoi, jpoj), STAT=ierr) ; CALL chk_allo (ierr, 'i2o')
118	!
119	ALLOCATE (tmask_atl (jpon), STAT=ierr) ; CALL chk_allo (ierr, 'tmask_atl')
120	ALLOCATE (tmask_noclose (jpon), STAT=ierr) ; CALL chk_allo (ierr, 'tmask_noclose')
121	ALLOCATE (tmask_nomed (jpon), STAT=ierr) ; CALL chk_allo (ierr, 'tmasl_nomed')
122	ALLOCATE (tmask_med (jpon), STAT=ierr) ; CALL chk_allo (ierr, 'tmask_med')
123	ALLOCATE (z2o (jpoi, jpoj), STAT=ierr) ; CALL chk_allo (ierr, 'z2o')
124	!
125	IF (l_limit_iosize) l_fulldiag = .FALSE.
126	!!
127	!! Initialisation
128	ngrd = 11 ; nsrf = 12 ; nmsk = 13 ; nchk = 14 ; ndeb = 9 ; nbug1 = 10
129	nwei4o2a = 81 ; nwei8o2a = 82 ; nwei4a2o = 83 ; nwei8a2o = 84
130	!
131	nout = 6
132	IF (nout /= 6 .AND. nout /= 0 ) &
133	& OPEN (unit = nout, file = 'cote' // TRIM(c_suffix) // '.out', action = 'WRITE', status = 'unknown', position = 'REWIND' )
134	!!
135	!! Diagnostics files
136	!!
137	IF (nout /= 6 ) &
138	& OPEN (unit = nout, file = 'cotes' // TRIM(c_suffix) // '.out', position = 'REWIND', action = 'WRITE', status = 'replace' )
139	!!
140	!! Read coordinates of all models, and weights
141	!!
142	CALL lecdata (lecweights = .FALSE., leco2amask = .TRUE. )
143
144	!! Lecture masques de bassin
145	IF ( TRIM (cotyp) == 'orca4' .OR. &
146	& TRIM (cotyp) == 'orca2' .OR. TRIM (cotyp) == 'orca2.0' .OR. TRIM (cotyp) == 'orca2.1' .OR. TRIM (cotyp) == 'orca2.3' .OR. &
147	& TRIM (cotyp) == 'orca1' &
148	& ) THEN
149	!cl_bassin = TRIM (cotyp)
150	cl_bassin = TRIM(c_basins)
151	!-$$ IF (c_period /= 'none' ) cl_bassin = TRIM (cl_bassin) // TRIM (c_period)
152	WRITE (nout,*) 'Reading bassin masks in : ', TRIM (cl_bassin)
153	CALL flioopfd (TRIM (cl_bassin), il_ncid )
154	CALL fliogetv (il_ncid, 'mask_atl' , z2o )
155	tmask_atl = RESHAPE (z2o, (/ jpon /) )
156	CALL fliogetv (il_ncid, 'mask_noclose', z2o )
157	tmask_noclose = RESHAPE (z2o, (/ jpon /) )
158	CALL fliogetv (il_ncid, 'mask_nomed' , z2o )
159	tmask_nomed = RESHAPE (z2o, (/ jpon /) )
160	CALL flioclo (il_ncid)
161	ELSE
162	tmask_atl = REAL (1-iomskt, rl)
163	tmask_noclose = REAL (1-iomskt, rl)
164	tmask_nomed = REAL (1-iomskt, rl)
165	END IF
166	!
167	WRITE (unit=nout, fmt=*) 'iomskt'
168	CALL prihin (RESHAPE(iomskt ,(/jpoi,jpoj/)), 1_il, nout)
169	WRITE (unit=nout, fmt=*) 'tmask_atl'
170	CALL prihin (RESHAPE(NINT(tmask_atl) ,(/jpoi,jpoj/)), 1_il, nout)
171	WRITE (unit=nout, fmt=*) 'tmask_med'
172	!CALL prihin (RESHAPE(NINT(tmask_med) ,(/jpoi,jpoj/)), 1_il, nout)
173	WRITE (unit=nout, fmt=*) 'tmask_noclose'
174	!CALL prihin (RESHAPE(NINT(tmask_noclose),(/jpoi,jpoj/)), 1_il, nout)
175	!
176	WRITE (nout,*) "iamskt"
177	iamskt = 0
178	WHERE (o2amask .LT. epsfrac) iamskt = 1
179	CALL prihin (RESHAPE(iamskt,(/jpai,jpaj/)), 1, nout)
180
181	! Set correct periodicity to grids
182	CALL lbc (xolont, noperio, 'T', jpoi, jpoj)
183	CALL lbc (xolonu, noperio, 'U', jpoi, jpoj)
184	CALL lbc (xolonv, noperio, 'V', jpoi, jpoj)
185	CALL lbc (xolonf, noperio, 'F', jpoi, jpoj)
186	CALL lbc (xolatt, noperio, 'T', jpoi, jpoj)
187	CALL lbc (xolatu, noperio, 'U', jpoi, jpoj)
188	CALL lbc (xolatv, noperio, 'V', jpoi, jpoj)
189	CALL lbc (xolatf, noperio, 'F', jpoi, jpoj)
190	CALL lbc (xosrft, noperio, 'T', jpoi, jpoj)
191	CALL lbc (iomskt, noperio, 'T', jpoi, jpoj)
192	!
193	CALL lbc (xalont, naperio, 'T', jpai, jpaj)
194	CALL lbc (xalonu, naperio, 'U', jpai, jpaj)
195	CALL lbc (xalonv, naperio, 'V', jpai, jpaj)
196	CALL lbc (xalatt, naperio, 'T', jpai, jpaj)
197	CALL lbc (xalatu, naperio, 'U', jpai, jpaj)
198	CALL lbc (xalatv, naperio, 'V', jpai, jpaj)
199	CALL lbc (xasrft, naperio, 'T', jpai, jpaj)
200	CALL lbc (iamskt, naperio, 'T', jpai, jpaj)
201	!!
202	CALL rectifie
203	!!
204	CALL fliocrfd ('itarget' // TRIM(c_suffix), (/'x', 'y'/), (/jpai, jpaj/), nid_it, mode='REPLACE')
205	CALL flioputa (nid_it, "?", "title" , "itarget" )
206	CALL flioputa (nid_it, '?', 'Comment', TRIM(c_comment) )
207	CALL DATE_AND_TIME (c_date, c_time, c_zone )
208	CALL flioputa (nid_it, "?", "history" , "Created: "//c_date(1:4)//"-"//c_date(5:6)//"-"// &
209	& c_date(7:8)//" "//c_time(1:2)//"h"//c_time(3:4)//" GMT"//TRIM(c_zone) )
210	CALL flioputa (nid_it, "?", "method" , "MOSAIC" )
211	CALL flioputa (nid_it, "?", "source_grid" , "curvilinear" )
212	CALL flioputa (nid_it, "?", "dest_grid" , "curvilinear" )
213	CALL flioputa (nid_it, "?", "Institution" , "IPSL" )
214	CALL flioputa (nid_it, "?", "Model" , "IPSL CM6" )
215	CALL GET_ENVIRONMENT_VARIABLE ( NAME="HOSTNAME", VALUE=c_tmp, TRIM_NAME=.TRUE., STATUS=ierr)
216	IF ( ierr == 0 ) THEN
217	CALL flioputa (nid_it, "?", "HOSTNAME" , TRIM(c_tmp) )
218	ELSE
219	WRITE (nout,*) 'Environment variable not found : $HOSTNAME'
220	END IF
221	CALL GET_ENVIRONMENT_VARIABLE ( NAME="LOGNAME" , VALUE=c_tmp, TRIM_NAME=.TRUE., STATUS=ierr)
222	IF ( ierr == 0 ) THEN
223	CALL flioputa (nid_it, "?", "LOGNAME" , TRIM(c_tmp) )
224	ELSE
225	WRITE (nout,*) 'Environment variable not found : $LOGNAME'
226	END IF
227
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, 'la_nearcoast' , (/1,2/), v_t=flio_i4)
232	CALL fliodefv (nid_it, 'dist_coast_a' , (/1,2/), v_t=flio_r4)
233
234	CALL fliodefv (nid_it, 'lon', (/1, 2/), units="degree_north", v_t=flio_r4 )
235	CALL fliodefv (nid_it, 'lat', (/1, 2/), units="degree_east", v_t=flio_r4 )
236	CALL fliodefv (nid_it, 'o2amask', (/1,2/), v_t=flio_r4)
237	CALL fliodefv (nid_it, 'iamskt', (/1,2/), v_t=flio_i4)
238	CALL fliodefv (nid_it, 'iamskp', (/1,2/), v_t=flio_i4)
239
240	CALL flioputv (nid_it, 'lon' , RESHAPE ( xalont, (/jpai, jpaj/)) )
241	CALL flioputv (nid_it, 'lat' , RESHAPE ( xalatt, (/jpai, jpaj/)) )
242	CALL flioputv (nid_it, 'o2amask', RESHAPE (o2amask, (/jpai, jpaj/)) )
243	CALL flioputv (nid_it, 'iamskt' , RESHAPE (iamskt , (/jpai, jpaj/)) )
244	CALL flioputv (nid_it, 'iamskp' , RESHAPE (iamskp , (/jpai, jpaj/)) )
245	!!
246
247	!!
248	CALL fliocrfd ('idest' // TRIM(c_suffix), (/'x', 'y'/), (/jpoi, jpoj/), nid_id, mode='REPLACE')
249	CALL flioputa (nid_id, "?", "title" , "idest" )
250	CALL flioputa (nid_id, '?', 'Comment', TRIM(c_comment) )
251	CALL DATE_AND_TIME (c_date, c_time, c_zone )
252	CALL flioputa (nid_id, "?", "history" , "Created: "//c_date(1:4)//"-"//c_date(5:6)//"-"// &
253	& c_date(7:8)//" "//c_time(1:2)//"h"//c_time(3:4)//" GMT"//TRIM(c_zone) )
254	CALL flioputa (nid_id, "?", "method" , "MOSAIC" )
255	CALL flioputa (nid_id, "?", "source_grid" , "curvilinear" )
256	CALL flioputa (nid_id, "?", "dest_grid" , "curvilinear" )
257	CALL flioputa (nid_id, "?", "Institution" , "IPSL" )
258	CALL flioputa (nid_id, "?", "Model" , "IPSL CM6" )
259	CALL GET_ENVIRONMENT_VARIABLE ( NAME="HOSTNAME", VALUE=c_tmp, TRIM_NAME=.TRUE., STATUS=ierr)
260	IF ( ierr == 0 ) THEN
261	CALL flioputa (nid_id, "?", "HOSTNAME" , TRIM(c_tmp) )
262	ELSE
263	WRITE (nout,*) 'Environment variable not found : $HOSTNAME'
264	END IF
265	CALL GET_ENVIRONMENT_VARIABLE ( NAME="LOGNAME" , VALUE=c_tmp, TRIM_NAME=.TRUE., STATUS=ierr)
266	IF ( ierr == 0 ) THEN
267	CALL flioputa (nid_id, "?", "LOGNAME" , TRIM(c_tmp) )
268	ELSE
269	WRITE (nout,*) 'Environment variable not found : $LOGNAME'
270	END IF
271
272	CALL fliodefv (nid_id, 'locot' , (/1,2/), v_t=flio_i4)
273	CALL fliodefv (nid_id, 'looce' , (/1,2/), v_t=flio_i4)
274	CALL fliodefv (nid_id, 'lo_nearcoast' , (/1,2/), v_t=flio_i4)
275	CALL fliodefv (nid_id, 'dist_coast_o' , (/1,2/), v_t=flio_r4)
276
277	CALL fliodefv (nid_id, 'lon', (/1, 2/), units="degree_north", v_t=flio_r4 )
278	CALL fliodefv (nid_id, 'lat', (/1, 2/), units="degree_east", v_t=flio_r4 )
279	CALL fliodefv (nid_id, 'iomskt', (/1,2/), v_t=flio_i4)
280	CALL fliodefv (nid_id, 'iomskp', (/1,2/), v_t=flio_i4)
281
282	CALL flioputv (nid_id, 'lon' , RESHAPE ( xolont, (/jpoi, jpoj/)) )
283	CALL flioputv (nid_id, 'lat' , RESHAPE ( xolatt, (/jpoi, jpoj/)) )
284	CALL flioputv (nid_id, 'iomskt' , RESHAPE (iamskt , (/jpoi, jpoj/)) )
285	CALL flioputv (nid_id, 'iomskp' , RESHAPE (iamskp , (/jpoi, jpoj/)) )
286
287	!!
288	!! Compute edges of grid boxes
289	!!
290	CALL bord_a
291	CALL bord_o
292	CALL bord_oa
293	!!
294	mo_name = "mozaic"
295	a2o_name = "a2o"
296	!!
297	WRITE (nout,*) 'Calculs points terre/oce/cote sur grille atmosphere'
298	WHERE (o2amask .LT. epsfrac ) laland = .TRUE.
299	WHERE (o2amask .GE. epsfrac .AND. o2amask .LE. (1.0_rl-epsfrac) ) lacot = .TRUE.
300	WHERE ( o2amask .GT. (1.0_rl-epsfrac) ) laoce = .TRUE.
301	!!
302	WRITE (UNIT=nout, FMT=*) 'Points cotiers ', COUNT(lacot)
303	WRITE (nout,*) 'lacot '
304	CALL prihbo (RESHAPE(lacot, (/jpai, jpaj/)), nout, l_norsud=.FALSE.)
305
306	WRITE (nout,*) 'Ecriture '
307	i2a = 0
308	WHERE (RESHAPE(laland,(/jpai,jpaj/))) i2a=1
309	CALL flioputv (nid_it, "laland", i2a)
310	i2a = 0
311	WHERE (RESHAPE(laoce ,(/jpai,jpaj/))) i2a=1
312	CALL flioputv (nid_it, "laoce", i2a)
313	i2a = 0
314	WHERE (RESHAPE(lacot ,(/jpai,jpaj/))) i2a=1
315	CALL flioputv (nid_it, "lacot", i2a)
316	!!
317	WRITE (nout,*) 'Determination des points cotes ocean '
318	locot = .FALSE.
319	looce = .FALSE.
320	DO jo = 1, jpon
321	IF ( iomskt (jo) == 0_il .AND. iomskp (jo) == 0_il ) THEN
322	joi = moi (jo) ; joj = moj (jo) ! Indices 2D
323	joi_p = MIN(joi+1, jpoi) ; joi_m = MAX (joi-1, 1)
324	joj_p = MIN(joj+1, jpoj) ; joj_m = MAX (joj-1, 1)
325	IF ( iomskt (m1o (joi_p, joj )) == 1 &
326	.OR. iomskt (m1o (joi_m, joj )) == 1 &
327	.OR. iomskt (m1o (joi , joj_p)) == 1 &
328	.OR. iomskt (m1o (joi , joj_m)) == 1 &
329	.OR. iomskt (m1o (joi_p, joj_p)) == 1 &
330	.OR. iomskt (m1o (joi_p, joj_m)) == 1 &
331	.OR. iomskt (m1o (joi_m, joj_p)) == 1 &
332	.OR. iomskt (m1o (joi_m, joj_m)) == 1 ) THEN
333	locot (jo) = .TRUE.
334	ELSE
335	looce (jo) = .TRUE.
336	END IF
337	END IF
338	END DO
339	CALL lbc (locot, noperio, 'T', jpoi, jpoj)
340	CALL lbc (looce, noperio, 'T', jpoi, jpoj)
341
342	i2o = 0
343	WHERE (RESHAPE(locot,(/jpoi,jpoj/))) i2o=1
344	CALL flioputv (nid_id, "locot", i2o)
345	i2o = 0
346	WHERE (RESHAPE(looce,(/jpoi,jpoj/))) i2o=1
347	CALL flioputv (nid_id, "looce", i2o)
348
349	WRITE (nout,*) 'locot '
350	CALL prihbo (RESHAPE(locot, (/jpoi, jpoj/)), nout, l_norsud=.FALSE.)
351
352	!!
353	!! Determination d'une bande cotiere atm
354	la_nearcoast = lacot
355	ALLOCATE (la_temp(jpan), STAT=ierr) ; CALL chk_allo (ierr, 'la_temp')
356	DO niter = 1, 5
357	la_temp = la_nearcoast
358	DO ja = 1, jpan
359	IF ( iamskt(ja) == 1_il .AND. iamskp (ja) == 0_il ) THEN
360	jai = mai (ja) ; jaj = maj (ja) ! Indices 2D
361	jai_p = MIN(jai+1, jpai) ; jai_m = MAX (jai-1, 1)
362	jaj_p = MIN(jaj+1, jpaj) ; jaj_m = MAX (jaj-1, 1)
363	IF ( la_temp (m1a (jai_p, jaj )) &
364	.OR. la_temp (m1a (jai_m, jaj )) &
365	.OR. la_temp (m1a (jai , jaj_p)) &
366	.OR. la_temp (m1a (jai , jaj_m)) &
367	.OR. la_temp (m1a (jai_p, jaj_p)) &
368	.OR. la_temp (m1a (jai_p, jaj_m)) &
369	.OR. la_temp (m1a (jai_m, jaj_p)) &
370	.OR. la_temp (m1a (jai_m, jaj_m)) ) THEN
371	la_nearcoast (ja) = .TRUE.
372	END IF
373	END IF
374	END DO
375	CALL lbc (la_nearcoast, naperio, 'T', jpai, jpaj)
376	END DO
377
378	WRITE (nout,*) 'la_nearcoast '
379	CALL prihbo (RESHAPE(la_nearcoast, (/jpai, jpaj/)), nout, l_norsud=.TRUE.)
380	i2a = 0
381	WHERE (RESHAPE(la_nearcoast ,(/jpai,jpaj/))) i2a=1
382	CALL flioputv (nid_it, "la_nearcoast", i2a)
383
384	!!
385	!! Determination d'une bande cotiere ocean
386	lo_nearcoast = locot
387	ALLOCATE (lo_temp(jpon), STAT=ierr) ; CALL chk_allo (ierr, 'lo_temp')
388	DO niter = 1, 8
389	lo_temp = lo_nearcoast
390	DO jo = 1, jpon
391	IF ( iomskt(jo) == 0 .AND. iomskp (jo) == 0_il ) THEN
392	joi = moi (jo) ; joj = moj (jo) ! Indices 2D
393	joi_p = MIN(joi+1, jpoi) ; joi_m = MAX(joi-1, 1)
394	joj_p = MIN(joj+1, jpoj) ; joj_m = MAX(joj-1, 1)
395	IF ( lo_temp (m1o (joi_p, joj )) &
396	.OR. lo_temp (m1o (joi_m, joj )) &
397	.OR. lo_temp (m1o (joi , joj_p)) &
398	.OR. lo_temp (m1o (joi , joj_m)) &
399	.OR. lo_temp (m1o (joi_p, joj_p)) &
400	.OR. lo_temp (m1o (joi_p, joj_m)) &
401	.OR. lo_temp (m1o (joi_m, joj_p)) &
402	.OR. lo_temp (m1o (joi_m, joj_m)) ) THEN
403	lo_nearcoast (jo) = .TRUE.
404	END IF
405	END IF
406	END DO
407	CALL lbc (lo_nearcoast, noperio, 'T', jpoi, jpoj)
408	END DO
409
410	WRITE (nout,*) 'lo_nearcoast '
411	CALL prihbo (RESHAPE(lo_nearcoast, (/jpoi, jpoj/)), nout, l_norsud=.FALSE.)
412	i2o = 0
413	WHERE (RESHAPE(lo_nearcoast ,(/jpoi,jpoj/))) i2o=1
414	CALL flioputv (nid_id, "lo_nearcoast", i2o)
415
416	!!
417	WRITE (nout,*) 'Distance a la cote - points atmosphere'
418	DO ja = 1, jpan
419	dist_coast_a (ja) = rpi * ra
420	IF ( MOD(ja, 1000) .EQ. 0 ) WRITE (nout,*) ja, '/', jpan
421	IF ( iamskt (ja) .EQ. 0 ) THEN
422	dist_coast_a (ja) = 0.0E0_rl
423	CYCLE
424	END IF
425	IF ( .NOT. la_nearcoast (ja) ) THEN
426	dist_coast_a (ja) = 2.0_rl * dist_max_atm
427	CYCLE
428	END IF
429	DO jo = 1, jpon
430	IF ( locot(jo) ) THEN
431	dist_coast_a (ja) = MIN ( dist_coast_a (ja), &
432	& ra * geodist ( xalont(ja), xalatt(ja), xolont(jo), xolatt(jo) ) )
433	END IF
434	END DO
435	END DO
436	CALL lbc ( dist_coast_a, naperio, 'T', jpai, jpaj )
437	CALL prihre (RESHAPE(dist_coast_a, (/jpai,jpaj/)), 1.0E-4, nout, l_norsud=.TRUE.)
438
439	CALL flioputv (nid_it, "dist_coast_a", RESHAPE(dist_coast_a, (/jpai,jpaj/)))
440	!!
441	WRITE (nout,*) 'Distance a la cote - point ocean'
442	DO jo = 1, jpon
443	IF ( MOD(jo, 5000) .EQ. 0 ) WRITE (nout,*) jo, '/', jpon
444	dist_coast_o (jo) = rpi * ra
445	IF ( iomskt (jo) .EQ. 1 ) THEN
446	dist_coast_o (jo) = 0.0E0_rl
447	CYCLE
448	END IF
449	IF ( .NOT. lo_nearcoast (jo) ) THEN
450	dist_coast_o (jo) = 2.0_rl * dist_max_oce
451	CYCLE
452	END IF
453
454	DO jo2 = 1, jpon
455	IF ( locot(jo2) ) THEN
456	dist_coast_o (jo) = MIN ( dist_coast_o (jo), &
457	& ra * geodist ( xolont(jo), xolatt(jo), xolont(jo2), xolatt(jo2) ) )
458	END IF
459	END DO
460	END DO
461	CALL lbc ( dist_coast_a, naperio, 'T', jpoi, jpoj )
462	CALL prihre (RESHAPE(dist_coast_o, (/jpoi,jpoj/)), 1.0E-4, nout)
463
464	CALL flioputv (nid_id, "dist_coast_o", RESHAPE(dist_coast_o, (/jpoi,jpoj/)))
465
466	WRITE (nout,*) 'Verif coherence (tout doit etre nul)'
467	WRITE (nout,*) COUNT ( laland .AND. lacot)
468	WRITE (nout,*) COUNT ( lacot .AND. laoce)
469	WRITE (nout,*) COUNT ( laoce .AND. laland)
470	WRITE (nout,*) COUNT ( .NOT. ( laland .OR. lacot .OR. laoce) )
471
472
473
474	!!
475	WRITE (unit = nout, fmt = *) 'Nombre de points cotes atmosphere : ', COUNT (lacot)
476	WRITE (unit = nout, fmt = *) 'Nombre de points cotes ocean : ', COUNT (locot)
477	!!
478	CALL bilan ("Avec cote ", c_case="complet", c_int_atm="Integral", c_int_oce="Integral")
479	CALL bilan ("Avec cote ", c_case="ocean" , c_int_atm="Integral", c_int_oce="Integral")
480	CALL bilan ("Avec cote ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral")
481
482
483	E_DRYRUN: IF (l_dryrun) THEN
484	WRITE (nout,*) "Remplissage bidon des champs pour tester les I/O"
485	jma2o = jma2or
486	nvo = jma2or
487	WRITE (nout,*) 'jma2o pour test : ', jma2o
488	CALL RANDOM_SEED
489	CALL RANDOM_NUMBER (wa2o) ; ka2o = NINT (wa2o*REAL(jpan,rl) )
490
491	DO jo = 1, jpon
492	IF ( iomskt (jo) == 1 ) THEN
493	wa2o (:, jo) = 0.0_rl
494	ka2o (:, jo) = 0_il
495	END IF
496	END DO
497
498	DO jo = 1, jpon
499	DO jno = 1, jpa2o
500	ja = ka2o (jno, jo)
501	IF (ja > 0) THEN
502	IF (iamskt (ja) == 1) THEN
503	wa2o (jno, jo) = 0.0_rl ; ka2o (jno, jo) = 0_il
504	END IF
505	END IF
506	END DO
507	END DO
508	ELSE
509
510	CALL cotes_a
511
512	WRITE (nout, *) 'Renormalization'
513	!-$$ RenormA : DO ja = 1, jpan
514	!-$$ z_d_1 = 0.0e0_rl
515	!-$$ kom = 0
516	!-$$ RenormO : DO jo = 1, jpon
517	!-$$ DO jn = 1, jpa2o
518	!-$$ IF ( ka2o (jn,jo) == ja ) THEN
519	!-$$ z_d_1 = z_d_1 + wa2o (jn,jo)
520	!-$$ kom = kom + 1
521	!-$$ jx (kom,:) = (/ jn, jo /)
522	!-$$ !!$ WRITE (unit = nout, fmt = '(6I6,16F18.6) ') ja, jo, jn, kom, jx (kom, :), &
523	!-$$ !!$ wa2o (jn,jo), wa2o (jx (kom, 1), jx (kom, 2)), z_d_1, &
524	!-$$ !!$ xosrft (jo), xasrft (ja)
525	!-$$ END IF
526	!-$$ END DO
527	!-$$ ENDDO RenormO
528	!-$$ IF (kom /= 0 ) THEN
529	!-$$ IF ( z_d_1 == 0.0_rl ) THEN
530	!-$$ WRITE (nout, *) 'z_d_1 nul ', ja, xasrft (ja), lacot (ja), iamskt (ja), &
531	!-$$ & wa2o (:,ja), ka2o (:,ja)
532	!-$$ STOP
533	!-$$ ELSE
534	!-$$ DO ko = 1, kom
535	!-$$ wa2o (jx (ko,1), jx (ko,2)) = wa2o (jx (ko,1), jx (ko,2)) / z_d_1
536	!-$$ END DO
537	!-$$ ENDIF
538	!-$$ END IF
539	!-$$ ENDDO RenormA
540	wasum = 0.0_rl
541	DO jo = 1, jpon
542	IF ( nvo(jo) /= 0) THEN
543	DO jn = 1, nvo(jo)
544	ja = ka2o (jn,jo)
545	wasum (ja) = wasum (ja) + wa2o (jn,jo)
546	END DO
547	END IF
548	END DO
549	!!
550	DO jo = 1, jpon
551	IF ( nvo(jo) /= 0 ) THEN
552	DO jn = 1, nvo(jo)
553	ja = ka2o (jn,jo)
554	IF ( wasum (ja) /= 0.0_rl ) wa2o (jn,jo) = wa2o (jn,jo) / wasum (ja)
555	END DO
556	END IF
557	END DO
558
559	!!
560	CALL int_wei
561	!!
562	!! Controle
563	!!
564	CALL verif (' Normalisation')
565	!!
566	CALL bilan ("Final ", c_case="terre" , c_int_atm="Integral", c_int_oce="Local")
567	CALL bilan ("Final ", c_case="terre100", c_int_atm="Integral", c_int_oce="Local")
568	CALL bilan ("Final ", c_case="cote" , c_int_atm="Integral", c_int_oce="Local")
569	!!
570	WRITE (unit = nout, fmt = * ) 'Nombre de voisins : ', jma2o
571	!!
572	CALL gather_wei ('apres normalisation')
573	CALL verif ('apres gather' )
574
575	!! Atm mask interpolated to oce
576	za = REAL (1_il-iamskt, KIND = rl )
577	CALL inter_a2o (za, a2omask )
578	!! 1 on Atm interpolated to oce
579	za = REAL (1_il , KIND = rl )
580	CALL inter_a2o (za, a2ofull )
581
582	!!
583	!! Checking weights
584	WRITE (unit = nout, fmt = '("Neighbors a2o : ", 3I9 )' ) &
585	& MINVAL (nvo, nvo /=0), MAXVAL (nvo), jma2o
586	WRITE (unit = nout, fmt = '("Index a2o : ", 2I9 )' ) &
587	& MINVAL (ka2o, ka2o /= 0), MAXVAL (ka2o)
588	jind = MINLOC(wa2o, wa2o /= 0.0_rl) ; jo = jind(2) ; jn = jind(1)
589	joi = m2oi(jo) ; joj = m2oj(jo)
590	WRITE (unit = nout, fmt = '("Weights a2o mini : ", 4I6, 1E13.4, 3I6)' ) &
591	& jo, jn, joi, joj, wa2o(jn,jo), ka2o(jn,jo), m2ai(ka2o(jn,jo)), m2aj(ka2o(jn,jo))
592	jind = MAXLOC(wa2o, wa2o /= 0.0_rl) ; jo = jind(2) ; jn = jind(1)
593	joi = m2oi(jo) ; joj = m2oj(jo)
594	WRITE (unit = nout, fmt = '("Weights a2o maxi : ", 4I6, 1E13.4, 3I6)' ) &
595	& jo, jn, joi, joj, wa2o(jn,jo), ka2o(jn,jo), m2ai(ka2o(jn,jo)), m2aj(ka2o(jn,jo))
596	!
597	jind(1:1) = INT(MAXLOC (nvo),il) ; jo = jind(1)
598	WRITE (unit = nout, fmt = *) jo, m2oi(jo), m2oj(jo), ka2o(1:nvo(jo),jo), &
599	& (m2ai(ka2o(jn,jo)), m2aj(ka2o(jn,jo)), jn = 1, nvo(jo))
600	!!
601	clweight = "WEIGHTS4" ; cladress = "ADRESSE4"
602	WRITE (unit = nout, fmt = *) cladress, " ", clweight, " ", jpa2o, jma2o
603	!!
604	!! Controle
605	!!
606	WRITE (nout, *) 'Verif avant ecriture'
607	WRITE (nout, *) 'wa2o ', MINVAL (wa2o), MAXVAL (wa2o), COUNT (wa2o > 0.0_rl)
608	WRITE (nout, *) 'ka2o ', MINVAL (ka2o), MAXVAL (ka2o), COUNT (ka2o > 0 )
609	WRITE (nout, *) 'xasrft ', MINVAL (xasrft), MAXVAL (xasrft)
610	WRITE (nout, *) 'o2amask ', MINVAL (o2amask), MAXVAL (o2amask), SUM (o2amask)
611	WRITE (nout, *) 'xosrft ', MINVAL (xosrft), MAXVAL (xosrft)
612	WRITE (nout, *) 'iomskt ', MINVAL (iomskt), MAXVAL (iomskt), COUNT (iomskt == 0)
613	WRITE (nout, *) 'iomskp ', MINVAL (iomskp), MAXVAL (iomskp), COUNT (iomskp == 0)
614	WRITE (nout, *) 'iamskt ', MINVAL (iamskt), MAXVAL (iamskt), COUNT (iamskt == 0)
615	WRITE (nout, *) 'locot ', COUNT (locot)
616	!!
617	!!
618	END IF E_DRYRUN
619	!!
620	!!
621	WRITE (nout, *) ' '
622	WRITE (nout, *) ' '
623	WRITE (nout, *) 'Eventuels traitement specifiques'
624	WRITE (nout, *) ' '
625	WRITE (nout, *) ' '
626
627
628	!! Extension ??
629
630	!! Traitements specifiques, tres specifiques ... !!
631	IF ( TRIM(c_suffix) == "_runoffNord09" ) THEN
632	WRITE (nout, *) 'Traitement specifique _runoffNord09'
633	DO jo = 1, jpon
634	IF ( tmask_atl(jo) == 1 .AND. xolatt(jo) > 40.0_rl ) wa2o(:,jo) = wa2o(:,jo) * 0.9
635	END DO
636	END IF
637	!
638	IF ( TRIM(c_suffix) == "_runoffNord07" ) THEN
639	WRITE (nout, *) 'Traitement specifique _runoffNord07'
640	DO jo = 1, jpon
641	IF ( tmask_atl(jo) == 1 .AND. xolatt(jo) > 40.0_rl ) wa2o(:,jo) = wa2o(:,jo) * 0.7
642	END DO
643	END IF
644
645	IF ( TRIM(c_suffix) == "_runoffNord00" ) THEN
646	WRITE (nout, *) 'Traitement specifique _runoffNord00'
647	DO jo = 1, jpon
648	IF ( tmask_atl(jo) == 1 .AND. xolatt(jo) > 40.0_rl ) wa2o(:,jo) = wa2o(:,jo) * 0.0
649	END DO
650	END IF
651	!!
652	!!
653	!! Verif nombre de voisins
654	!!
655	IF ( jma2or /= jma2o ) THEN
656	WRITE (nout, *) 'Erreur nombre de voisins dans run.def : '
657	WRITE (nout, *) ' jma2o calcule / lu : ', jma2o, jma2or
658	!STOP
659	ENDIF
660	!!
661	!-$$ !! Ecriture des poids
662	!-$$ !!
663	!-$$ WRITE (unit = c_r, fmt = '(1i1)' ) rk_out
664	!-$$ !
665	!-$$ IF ( l_wei_i4) THEN
666	!-$$ WRITE (unit = c_i, fmt = '(1i1)' ) i_4
667	!-$$ cfname4 = TRIM(mo_name) // ".wa2o.runoff.i" // TRIM(c_i) // ".r" // TRIM(c_r) // TRIM(c_suffix)
668	!-$$ OPEN (unit = nwei4a2o, file = TRIM(cfname4), form = 'unformatted', status = 'replace',&
669	!-$$ & action = 'write')
670	!-$$ WRITE (nout, *) 'Poids i_4 a -> o: ', TRIM(cfname4)
671	!-$$ END IF
672	!-$$ !
673	!-$$ IF (l_wei_i8) THEN
674	!-$$ WRITE (unit = c_i, fmt = '(1i1)' ) i_8
675	!-$$ cfname8 = TRIM(mo_name) // ".wa2o.runoff.i" // TRIM(c_i) // ".r" // TRIM(c_r) // TRIM(c_suffix)
676	!-$$ OPEN (unit = nwei8a2o, file = TRIM(cfname8), form = 'unformatted', status = 'replace',&
677	!-$$ & action = 'write')
678	!-$$ WRITE (nout, *) 'Poids i_8 a -> o : ', TRIM(cfname8)
679	!-$$ END IF
680	!-$$ !
681	!-$$ cldiag_a2o = TRIM(a2o_name) // '.runcoa.diag'
682	!-$$ clw_a2o = TRIM(mo_name) // '.wa2o.runcoa'
683	!-$$ clw_a2o_mct= TRIM(mo_name) // '.wa2o_mct.runcoa'
684	!-$$ !
685	!-$$ CALL wri_weights_a2o (cldiag_a2o, clw_a2o, clw_a2o_mct, "4", l_fulldiag, lo_src_grid_frac=.FALSE., lo_dst_grid_frac=.FALSE., &
686	!-$$ & co_omsk=TRIM(cotes_omsk), co_amsk=TRIM(cotes_amsk) )
687	!-$$ !
688
689
690	!!
691	!! Ecriture des poids
692	!!
693	!! Ouverture des fichiers poids
694	!!
695	WRITE (unit = c_r, fmt = '(1i1)' ) rk_out
696	!
697	IF (l_wei_i4) THEN
698	WRITE (unit = c_i, fmt = '(1i1)' ) i_4
699	cfname4 = TRIM(mo_name) // ".wa2o.rivflu" // TRIM(c_suffix) // ".i" // TRIM(c_i) // ".r" // TRIM(c_r)
700	OPEN (unit = nwei4a2o, file = TRIM(cfname4), form = 'unformatted', status = 'replace',&
701	& action = 'write')
702	WRITE (nout, *) 'Poids i4 a -> o: ', TRIM(cfname4)
703	END IF
704	!
705	IF (l_wei_i4) THEN
706	WRITE (unit = c_i, fmt = '(1i1)' ) i_8
707	cfname8 = TRIM (mo_name) // ".wa2o.rivflu" // TRIM(c_suffix) // ".i" // TRIM(c_i) // ".r" // TRIM(c_r)
708	OPEN (unit = nwei8a2o, file = TRIM(cfname8), form = 'unformatted', status = 'replace',&
709	& action = 'write')
710	WRITE (nout, *) 'Poids i8 a -> o : ', TRIM(cfname8)
711	END IF
712
713	!
714	cldiag_a2o = TRIM (a2o_name) // '.rivflu.diag'
715	clw_a2o = TRIM (mo_name) // '.wa2o.rivflu'
716	clw_a2o_mct = "rmp_"//TRIM(camod_t)//"_to_"//TRIM(comod_t)//"_MOSAIC_rivflu"
717	!
718	CALL wri_weights_a2o (cldiag_a2o, clw_a2o, clw_a2o_mct, "3", l_fulldiag, lo_src_grid_frac=.FALSE., lo_dst_grid_frac=.FALSE., &
719	& co_omsk='noperio', co_amsk='full' )
720
721	!!
722	IF (l_wei_i4) THEN
723	WRITE (nout,*) 'Fermeture poids i4'
724	CLOSE (unit = nwei8a2o)
725	END IF
726	IF (l_wei_i8) THEN
727	WRITE (nout,*) 'Fermeture poids i8 '
728	CLOSE (unit = nwei4a2o)
729	END IF
730	!!
731	WRITE (nout,*) 'Appel flioclo'
732	CALL flioclo
733	!!
734	STOP
735	!!
736	CONTAINS
737	!!
738	SUBROUTINE diag_itarget (clmess, lwrite)
739	CHARACTER (len=*), INTENT (in) :: clmess
740	LOGICAL, INTENT (in) :: lwrite
741	INTEGER, SAVE :: nappel = 0
742	CHARACTER (len=2) :: cc
743	!!
744	nappel = nappel + 1
745	!!
746	itarget = 0_il ; wtarget = 0.0_rl
747	DO jo = 1, jpon
748	DO jn = 1, jpa2o
749	ja = ka2o (jn, jo)
750	IF (ja /= 0 ) THEN
751	itarget (ja) = itarget (ja) + 1_il
752	wtarget (ja) = wtarget (ja) + wa2o (jn, jo)
753	ENDIF
754	END DO
755	END DO
756	!!
757	!CALL lbc (itarget, naperio, 'T', jpai, jpaj)
758	!CALL lbc (wtarget, naperio, 'T', jpai, jpaj)
759	!!
760	IF (lwrite) THEN
761	!$$$ WRITE(nout,*) 'itarget, etape '// clmess
762	!$$$ CALL prihin (RESHAPE(itarget,(/jpai,jpaj/)),1_il,nout)
763	WRITE (unit = cc, fmt = '(1I2.2)' ) nappel
764	CALL fliodefv (nid_it, 'itarget' // cc, (/1, 2/), v_t=flio_i4)
765	CALL fliodefv (nid_it, 'ntarget' // cc, (/1, 2/), v_t=flio_i4)
766	CALL fliodefv (nid_it, 'wtarget' // cc, (/1, 2/), v_t=flio_r4)
767	CALL fliodefv (nid_it, 'non_target' // cc, (/1, 2/), v_t=flio_i4)
768
769
770	CALL flioputa (nid_it, 'itarget' // cc, 'long_name ', clmess // ' itarget' // cc)
771	CALL flioputv (nid_it, 'itarget' // cc, MIN (1, RESHAPE (itarget, (/jpai, jpaj/))) )
772
773	CALL flioputa (nid_it, 'ntarget' // cc, 'long_name ', clmess // ' ntarget' // cc)
774	CALL flioputv (nid_it, 'ntarget' // cc, RESHAPE (itarget, (/jpai, jpaj/)) )
775
776	CALL flioputa (nid_it, 'wtarget' // cc, 'long_name ', clmess // ' wtarget' // cc)
777	CALL flioputv (nid_it, 'wtarget' // cc, RESHAPE (wtarget, (/jpai, jpaj/)) )
778
779	wtarget = 0
780	WHERE (itarget .EQ. 0 ) wtarget = 1.0_rl
781	CALL flioputa (nid_it, 'non_target' // cc, 'long_name ', clmess // ' non_target' // cc)
782	CALL flioputv (nid_it, 'non_target' //cc, RESHAPE (NINT (wtarget), (/jpai, jpaj/)) )
783	END IF
784	!!
785	END SUBROUTINE diag_itarget
786	!!
787	SUBROUTINE diag_idest (clmess, lwrite)
788	CHARACTER (len=*), INTENT (in) :: clmess
789	LOGICAL, INTENT (in) :: lwrite
790	INTEGER, SAVE :: nappel = 0
791	CHARACTER (len=2) :: cc
792	!!
793	nappel = nappel + 1
794	!!
795	idest = 0_il ; wdest = 0.0_rl
796	DO jo = 1, jpon
797	DO jn = 1, jpa2o
798	ja = ka2o (jn, jo)
799	IF (ja /= 0 ) THEN
800	idest (jo) = idest (jo) + 1_il
801	wdest (jo) = wdest (jo) + wa2o (jn, jo)
802	ENDIF
803	END DO
804	END DO
805	!!
806	!CALL lbc (itarget, naperio, 'T', jpai, jpaj)
807	!CALL lbc (wtarget, naperio, 'T', jpai, jpaj)
808	!!
809	IF (lwrite) THEN
810	!$$$ WRITE(nout,*) 'idest, etape '// clmess
811	!$$$ CALL prihin (RESHAPE(idest,(/jpoi,jpoj/)),1_il,nout)
812	WRITE (unit = cc, fmt = '(1I2.2)' ) nappel
813	CALL fliodefv (nid_id, 'idest' // cc, (/1, 2/), v_t=flio_i4)
814	CALL fliodefv (nid_id, 'ndest' // cc, (/1, 2/), v_t=flio_i4)
815	CALL fliodefv (nid_id, 'wdest' // cc, (/1, 2/), v_t=flio_r4)
816	CALL fliodefv (nid_id, 'non_dest' // cc, (/1, 2/), v_t=flio_i4)
817
818
819	CALL flioputa (nid_id, 'idest' // cc, 'long_name ', clmess // ' idest' // cc)
820	CALL flioputv (nid_id, 'idest' // cc, MIN (1, RESHAPE (idest, (/jpoi, jpoj/))) )
821
822	CALL flioputa (nid_id, 'ndest' // cc, 'long_name ', clmess // ' ndest' // cc)
823	CALL flioputv (nid_id, 'ndest' // cc, RESHAPE (idest, (/jpoi, jpoj/)) )
824
825	CALL flioputa (nid_id, 'wdest' // cc, 'long_name ', clmess // ' wdest' // cc)
826	CALL flioputv (nid_id, 'wdest' // cc, RESHAPE (wdest, (/jpoi, jpoj/)) )
827
828	wdest = 0
829	WHERE (idest .EQ. 0 ) wdest = 1.0_rl
830	CALL flioputa (nid_id, 'non_dest' // cc, 'long_name ', clmess // ' non_dest' // cc)
831	CALL flioputv (nid_id, 'non_dest' //cc, RESHAPE (NINT (wdest), (/jpoi, jpoj/)) )
832	END IF
833	!!
834	END SUBROUTINE diag_idest
835	!!
836	!!
837	SUBROUTINE int_wei
838	!!
839	IMPLICIT NONE
840	WRITE(nout,*) 'Debut normalisation'
841	!!
842	IF (.NOT. lint_atm) THEN
843	WRITE(nout,*) 'Multiplication par les surfaces atm'
844	DO jo = 1, jpon
845	DO jn = 1, jma2o
846	ja = ka2o (jn, jo)
847	IF (ja > 0 ) wa2o (jn,jo) = wa2o (jn,jo) * xasrft (ja)
848	END DO
849	ENDDO
850	ENDIF
851	!!
852	IF (.NOT. lint_oce) THEN
853	WRITE (nout,*) 'Division par les surfaces oce'
854	DO jo = 1, jpon
855	DO jn = 1, jma2o
856	wa2o (jn,jo) = wa2o (jn,jo) / xosrft (jo)
857	END DO
858	ENDDO
859	ENDIF
860	END SUBROUTINE int_wei
861	!!
862	SUBROUTINE complet_run (ka1, ka2, kn)
863	INTEGER (kind=il), INTENT(in) :: ka1, ka2 ! Index of atm boxes
864	INTEGER (kind=il), INTENT(in) :: kn ! Number of neighbors
865	INTEGER (kind=il) :: jt
866	!!
867	DO jo = 1, jpon
868	IF (iomskt(jo) == 1 ) CYCLE
869	IF (iomskp(jo) .EQ. 1_il ) CYCLE
870	DO jn = 1, jpa2o
871	IF (ka2o(jn,jo) == ka1 ) THEN
872	nvo(jo) = nvo(jo)+1
873	IF (nvo(jo) > jpa2o ) THEN
874	WRITE(nout,*) 'Pas assez de voisins pour la completion '
875	WRITE(nout,fmt='(1I7,2I4,2F10.1)' ) ka1, m2ai(ka1), m2aj(ka1), &
876	& xalont(ka1), xalatt(ka1)
877	WRITE(nout,fmt='(1I7,2I4,2F10.1)' ) ka2, m2ai(ka2), m2aj(ka2), &
878	& xalont(ka2), xalatt(ka2)
879	WRITE(nout,fmt='(1I7,2I4,2F10.1,2I8)') jo, m2oi(jo), m2oj(jo), &
880	& xolont(jo), xolatt(jo), jn, nvo(jo)
881	DO jt = 1, jpa2o
882	WRITE(nout,fmt=*) jt, ka2o(jt,jo), xalont(ka2o(jt,jo)), xalatt(ka2o(jt,jo))
883	ENDDO
884	END IF
885	ka2o (nvo (jo), jo) = ka2
886	wa2o (nvo (jo), jo) = wa2o(jn,jo) * xasrft(ka2) / xasrft(ka1) / REAL(kn,KIND=rl)
887	END IF
888	END DO
889	ENDDO
890	END SUBROUTINE complet_run
891	!!
892	SUBROUTINE gather_wei (clmess)
893	CHARACTER (len=*), INTENT (in) :: clmess
894	INTEGER (kind=il) :: jn1, jn2, jo2
895	!! Eliminate redundancy
896	DO jo = 1, jpon
897	DO jn1 = 1, jpa2o-1
898	IF (ka2o (jn1,jo) == 0 ) CYCLE
899	DO jn2 = jn1+1, jpa2o
900	IF (ka2o (jn1, jo) == ka2o (jn2, jo) ) THEN
901	WRITE(nout,*) 'Redondance ', jo, jn1, jn2, ka2o(jn1,jo), ka2o (jn2,jo), wa2o(jn1,jo), wa2o(jn2,jo)
902	wa2o (jn1, jo) = wa2o (jn1, jo) + wa2o (jn2, jo)
903	ka2o (jn2, jo) = 0
904	END IF
905	END DO
906	END DO
907	END DO
908	!!
909	!! Seek for strangeness
910	DO jo = 1, jpon
911	DO jn = 1, jpa2o
912	ja = ka2o (jn, jo)
913	IF (ja == 0) CYCLE
914	IF (wa2o (jn, jo) == 0.0E0 ) THEN
915	WRITE (nout, *) 'Null weight : ', jo, jn, ja, locot (jo)
916	END IF
917	END DO
918	END DO
919	!!
920	!! Periodicity
921	!DO jn = 1, jpa2o
922	! CALL lbc (wa2o (jn,:), noperio, 'T', jpoi, jpoj)
923	! CALL lbc (ka2o (jn,:), noperio, 'T', jpoi, jpoj)
924	!END DO
925	CALL diag_itarget (' gather/perio - ' // clmess, lwrite=.TRUE.)
926	CALL diag_idest (' gather/perio - ' // clmess, lwrite=.TRUE.)
927	!!
928	!! Gather weights
929	DO jo = 1, jpon
930	ktemp = ka2o (:,jo) ; wtemp = wa2o (:,jo)
931	ka2o (:,jo) = 0_il ; wa2o (:,jo) = 0.0_rl
932	kn = 1_il
933	DO jo2 = 1, jpa2o
934	IF (ktemp (jo2) /= 0_il ) THEN
935	ka2o (kn, jo) = ktemp (jo2)
936	wa2o (kn, jo) = wtemp (jo2)
937	kn = kn + 1_il
938	END IF
939	END DO
940	END DO
941	nvo = COUNT (ka2o /= 0_il, DIM = 1)
942	jma2o = MAXVAL (nvo)
943	!
944	!!
945	!! Periodicity
946	!DO jn = 1, jpa2o
947	! CALL lbc (wa2o(jn,:), noperio, 'T', jpoi, jpoj)
948	! CALL lbc (ka2o(jn,:), noperio, 'T', jpoi, jpoj)
949	!END DO
950	!!
951	nvo = COUNT (ka2o /= 0_il, DIM = 1)
952	jma2o = MAXVAL (nvo)
953	!!
954	END SUBROUTINE gather_wei
955	!!
956	SUBROUTINE verif (clmess, l_non_a)
957	!!
958	CHARACTER (len=*), INTENT (in) :: clmess
959	INTEGER :: ja
960	LOGICAL, INTENT (in), OPTIONAL :: l_non_a
961	!!
962	nvo = COUNT (ka2o /= 0_il, DIM = 1)
963	jma2o = MAXVAL (nvo)
964	CALL diag_itarget (clmess, lwrite=.TRUE.)
965	CALL diag_idest (clmess, lwrite=.TRUE.)
966	!!
967	WRITE (nout, *) ' --- Verifs -- : ', TRIM (clmess)
968	WRITE (nout, *) 'wa2o ', MINVAL (wa2o), MAXVAL (wa2o), COUNT (wa2o > 0.0_rl)
969	WRITE (nout, *) 'ka2o ', MINVAL (ka2o), MAXVAL (ka2o), COUNT (ka2o > 0 )
970	WRITE (nout, *) 'xasrft ', MINVAL (xasrft), MAXVAL (xasrft)
971	WRITE (nout, *) 'o2amask ', MINVAL (o2amask), MAXVAL (o2amask), SUM (o2amask)
972	WRITE (nout, *) 'xosrft ', MINVAL (xosrft), MAXVAL (xosrft)
973	WRITE (nout, *) 'iomskt ', MINVAL (iomskt), MAXVAL (iomskt), COUNT (iomskt == 0)
974	WRITE (nout, *) 'iomskp ', MINVAL (iomskp), MAXVAL (iomskp), COUNT (iomskp == 0)
975	WRITE (nout, *) 'iamskt ', MINVAL (iamskt), MAXVAL (iamskt), COUNT (iamskt == 0)
976	WRITE (nout, *) 'locot ', COUNT (locot)
977	WRITE (nout, *) 'Points non attribues : ', COUNT (itarget .EQ. 0 )
978
979	IF (PRESENT (l_non_a)) THEN
980	DO ja = 1, jpan
981	IF (itarget (ja) == 0 .AND. iamskp(ja) == 0 .AND. (laland (ja) .OR. lacot(ja) )) THEN
982	!-$$ WRITE (UNIT=nout,FMT='("Point atm non attribues : ", 3I7, 2F7.1, 1I4)' ) ja, m2ai(ja), m2aj (ja), xalont(ja), xalatt(ja), itarget(ja)
983	WRITE (nout,*) ja, m2ai(ja), m2aj (ja), xalont(ja), xalatt(ja), itarget(ja), laland(ja), lacot(ja), laoce(ja)
984	ENDIF
985	END DO
986	END IF
987	END SUBROUTINE verif
988	!!
989	!!
990	SUBROUTINE noroute
991	!! Determine les points atmospheres non routes
992	lanoroute = .TRUE.
993
994	DO jn = 1, jpa2o
995	DO jo = 1, jpon
996	ja = ka2o ( jo, jn)
997	IF ( ja /= 0 .AND. wa2o (jo,jn) > 0.0_rl ) THEN
998	lanoroute = .FALSE.
999	END IF
1000	END DO
1001	END DO
1002
1003	END SUBROUTINE noroute
1004	!!
1005	SUBROUTINE bilan (clmess, c_case, c_int_atm, c_int_oce)
1006	!!
1007	CHARACTER (len=*), INTENT (in) :: clmess
1008	CHARACTER (len=*), INTENT (in) :: c_case ! Type de champ d'entree
1009	CHARACTER (len=*), INTENT (in) :: c_int_atm ! Type d'integrale sur l'atm
1010	CHARACTER (len=*), INTENT (in) :: c_int_oce ! Type d'integrale sur l'ocean
1011	!!
1012	WRITE (nout,*) "Bilan :" // TRIM (clmess) // ", cas :", TRIM (c_case), ", c_int_atm :", TRIM (c_int_atm),&
1013	& ", c_int_oce :", TRIM(c_int_oce), "."
1014	!!
1015	za = 0.0_rl
1016	!!
1017	SELECT CASE (TRIM (c_case))
1018	CASE ('complet')
1019	za = 1.0_rl * REAL (1-iamskp, KIND=rl)
1020	SELECT CASE (TRIM (c_int_atm))
1021	CASE ('Integral')
1022	za = za * xasrft
1023	zsuma = SUM (za)
1024	CASE ('Local')
1025	zsuma = DOT_PRODUCT (za, xasrft)
1026	CASE default
1027	WRITE (nout,*) 'Erreur type de bilan atm' ; STOP
1028	END SELECT
1029	!!
1030	CASE ('terre')
1031	WHERE (laland .OR. lacot)
1032	za = 1.0_rl * REAL (1-iamskp, KIND=rl)
1033	END WHERE
1034	SELECT CASE (TRIM (c_int_atm))
1035	CASE ('Integral')
1036	za = za * xasrft
1037	zsuma = SUM (za)
1038	CASE ('Local')
1039	zsuma = DOT_PRODUCT (za, xasrft*(1.0-o2amask))
1040	CASE default
1041	WRITE (nout,*) 'Erreur type de bilan atm' ; STOP
1042	END SELECT
1043	!!
1044	CASE ('terre100')
1045	WHERE (laland)
1046	za = 1.0_rl * REAL (1-iamskp, KIND=rl)
1047	END WHERE
1048	SELECT CASE (TRIM (c_int_atm))
1049	CASE ('Integral')
1050	za = za * xasrft
1051	zsuma = SUM (za)
1052	CASE ('Local')
1053	zsuma = DOT_PRODUCT (za, xasrft)
1054	CASE default
1055	WRITE (nout,*) 'Erreur type de bilan atm' ; STOP
1056	END SELECT
1057	!!
1058	CASE ('ocean')
1059	WHERE (lacot)
1060	za = 1.0_rl * REAL (1-iamskp, KIND=rl)
1061	END WHERE
1062	SELECT CASE (TRIM (c_int_atm))
1063	CASE ('Integral')
1064	za = za * xasrft
1065	zsuma = SUM (za)
1066	CASE ('Local')
1067	zsuma = DOT_PRODUCT (za, xasrft*o2amask)
1068	CASE default
1069	WRITE (nout,*) TRIM(clmess), ' : Erreur type de bilan atm' ; STOP
1070	END SELECT
1071	!!
1072	CASE ('cote')
1073	WHERE (lacot)
1074	za = 1.0_rl * REAL (1-iamskp, KIND=rl)
1075	END WHERE
1076	SELECT CASE (TRIM (c_int_atm))
1077	CASE ('Integral')
1078	za = za * xasrft
1079	zsuma = SUM (za)
1080	CASE ('Local')
1081	zsuma = DOT_PRODUCT (za, xasrft*o2amask)
1082	CASE default
1083	WRITE (nout,*) 'Erreur type de bilan atm' ; STOP
1084	END SELECT
1085	!!
1086	CASE default
1087	WRITE (nout,*)'Erreur choix de cas' ; STOP
1088	STOP
1089	END SELECT
1090	!!
1091	CALL inter_a2o (za, zo)
1092	!!
1093	SELECT CASE (TRIM (c_int_oce))
1094	CASE ('Integral')
1095	zsumo = DOT_PRODUCT (zo, REAL(1-iomskp, KIND=rl))
1096	CASE ('Local')
1097	zsumo = DOT_PRODUCT (zo*xosrft,REAL(1-iomskp, KIND=rl))
1098	CASE default
1099	WRITE (nout,*) TRIM (clmess), ' : Erreur type de bilan ' ; STOP
1100	END SELECT
1101	!!
1102	WRITE (nout, fmt='(" bilan atm: ", 1E15.6, ", bilan oce: ", 1E15.6)' ) zsuma, zsumo
1103	!!
1104	RETURN
1105	!!
1106	END SUBROUTINE bilan
1107	!!
1108	END PROGRAM cotes_etal
1109

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

Download in other formats: