Context Navigation

cotes_etal.f90 @ 3363

Last change on this file since 3363 was 3363, checked in by omamce, 7 years ago
O.M. : add SVN keywords in files
File size: 42.2 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	CALL flioputa (il_ncid, "?", "SVN_Author" , TRIM (SVN_Author) )
228	CALL flioputa (il_ncid, "?", "SVN_Date" , TRIM (SVN_Date) )
229	CALL flioputa (il_ncid, "?", "SVN_Revision" , TRIM (SVN_Revision) )
230
231	CALL fliodefv (nid_it, 'laland' , (/1,2/), v_t=flio_i4)
232	CALL fliodefv (nid_it, 'lacot' , (/1,2/), v_t=flio_i4)
233	CALL fliodefv (nid_it, 'laoce' , (/1,2/), v_t=flio_i4)
234	CALL fliodefv (nid_it, 'la_nearcoast' , (/1,2/), v_t=flio_i4)
235	CALL fliodefv (nid_it, 'dist_coast_a' , (/1,2/), v_t=flio_r4)
236
237	CALL fliodefv (nid_it, 'lon', (/1, 2/), units="degree_north", v_t=flio_r4 )
238	CALL fliodefv (nid_it, 'lat', (/1, 2/), units="degree_east", v_t=flio_r4 )
239	CALL fliodefv (nid_it, 'o2amask', (/1,2/), v_t=flio_r4)
240	CALL fliodefv (nid_it, 'iamskt', (/1,2/), v_t=flio_i4)
241	CALL fliodefv (nid_it, 'iamskp', (/1,2/), v_t=flio_i4)
242
243	CALL flioputv (nid_it, 'lon' , RESHAPE ( xalont, (/jpai, jpaj/)) )
244	CALL flioputv (nid_it, 'lat' , RESHAPE ( xalatt, (/jpai, jpaj/)) )
245	CALL flioputv (nid_it, 'o2amask', RESHAPE (o2amask, (/jpai, jpaj/)) )
246	CALL flioputv (nid_it, 'iamskt' , RESHAPE (iamskt , (/jpai, jpaj/)) )
247	CALL flioputv (nid_it, 'iamskp' , RESHAPE (iamskp , (/jpai, jpaj/)) )
248	!!
249
250	!!
251	CALL fliocrfd ('idest' // TRIM(c_suffix), (/'x', 'y'/), (/jpoi, jpoj/), nid_id, mode='REPLACE')
252	CALL flioputa (nid_id, "?", "title" , "idest" )
253	CALL flioputa (nid_id, '?', 'Comment', TRIM(c_comment) )
254	CALL DATE_AND_TIME (c_date, c_time, c_zone )
255	CALL flioputa (nid_id, "?", "history" , "Created: "//c_date(1:4)//"-"//c_date(5:6)//"-"// &
256	& c_date(7:8)//" "//c_time(1:2)//"h"//c_time(3:4)//" GMT"//TRIM(c_zone) )
257	CALL flioputa (nid_id, "?", "method" , "MOSAIC" )
258	CALL flioputa (nid_id, "?", "source_grid" , "curvilinear" )
259	CALL flioputa (nid_id, "?", "dest_grid" , "curvilinear" )
260	CALL flioputa (nid_id, "?", "Institution" , "IPSL" )
261	CALL flioputa (nid_id, "?", "Model" , "IPSL CM6" )
262	CALL GET_ENVIRONMENT_VARIABLE ( NAME="HOSTNAME", VALUE=c_tmp, TRIM_NAME=.TRUE., STATUS=ierr)
263	IF ( ierr == 0 ) THEN
264	CALL flioputa (nid_id, "?", "HOSTNAME" , TRIM(c_tmp) )
265	ELSE
266	WRITE (nout,*) 'Environment variable not found : $HOSTNAME'
267	END IF
268	CALL GET_ENVIRONMENT_VARIABLE ( NAME="LOGNAME" , VALUE=c_tmp, TRIM_NAME=.TRUE., STATUS=ierr)
269	IF ( ierr == 0 ) THEN
270	CALL flioputa (nid_id, "?", "LOGNAME" , TRIM(c_tmp) )
271	ELSE
272	WRITE (nout,*) 'Environment variable not found : $LOGNAME'
273	END IF
274	CALL flioputa (il_ncid, "?", "SVN_Author" , TRIM (SVN_Author) )
275	CALL flioputa (il_ncid, "?", "SVN_Date" , TRIM (SVN_Date) )
276	CALL flioputa (il_ncid, "?", "SVN_Revision" , TRIM (SVN_Revision) )
277
278	CALL fliodefv (nid_id, 'locot' , (/1,2/), v_t=flio_i4)
279	CALL fliodefv (nid_id, 'looce' , (/1,2/), v_t=flio_i4)
280	CALL fliodefv (nid_id, 'lo_nearcoast' , (/1,2/), v_t=flio_i4)
281	CALL fliodefv (nid_id, 'dist_coast_o' , (/1,2/), v_t=flio_r4)
282
283	CALL fliodefv (nid_id, 'lon', (/1, 2/), units="degree_north", v_t=flio_r4 )
284	CALL fliodefv (nid_id, 'lat', (/1, 2/), units="degree_east", v_t=flio_r4 )
285	CALL fliodefv (nid_id, 'iomskt', (/1,2/), v_t=flio_i4)
286	CALL fliodefv (nid_id, 'iomskp', (/1,2/), v_t=flio_i4)
287
288	CALL flioputv (nid_id, 'lon' , RESHAPE ( xolont, (/jpoi, jpoj/)) )
289	CALL flioputv (nid_id, 'lat' , RESHAPE ( xolatt, (/jpoi, jpoj/)) )
290	CALL flioputv (nid_id, 'iomskt' , RESHAPE (iamskt , (/jpoi, jpoj/)) )
291	CALL flioputv (nid_id, 'iomskp' , RESHAPE (iamskp , (/jpoi, jpoj/)) )
292
293	!!
294	!! Compute edges of grid boxes
295	!!
296	CALL bord_a
297	CALL bord_o
298	CALL bord_oa
299	!!
300	mo_name = "mozaic"
301	a2o_name = "a2o"
302	!!
303	WRITE (nout,*) 'Calculs points terre/oce/cote sur grille atmosphere'
304	WHERE (o2amask .LT. epsfrac ) laland = .TRUE.
305	WHERE (o2amask .GE. epsfrac .AND. o2amask .LE. (1.0_rl-epsfrac) ) lacot = .TRUE.
306	WHERE ( o2amask .GT. (1.0_rl-epsfrac) ) laoce = .TRUE.
307	!!
308	WRITE (UNIT=nout, FMT=*) 'Points cotiers ', COUNT(lacot)
309	WRITE (nout,*) 'lacot '
310	CALL prihbo (RESHAPE(lacot, (/jpai, jpaj/)), nout, l_norsud=.FALSE.)
311
312	WRITE (nout,*) 'Ecriture '
313	i2a = 0
314	WHERE (RESHAPE(laland,(/jpai,jpaj/))) i2a=1
315	CALL flioputv (nid_it, "laland", i2a)
316	i2a = 0
317	WHERE (RESHAPE(laoce ,(/jpai,jpaj/))) i2a=1
318	CALL flioputv (nid_it, "laoce", i2a)
319	i2a = 0
320	WHERE (RESHAPE(lacot ,(/jpai,jpaj/))) i2a=1
321	CALL flioputv (nid_it, "lacot", i2a)
322	!!
323	WRITE (nout,*) 'Determination des points cotes ocean '
324	locot = .FALSE.
325	looce = .FALSE.
326	DO jo = 1, jpon
327	IF ( iomskt (jo) == 0_il .AND. iomskp (jo) == 0_il ) THEN
328	joi = moi (jo) ; joj = moj (jo) ! Indices 2D
329	joi_p = MIN(joi+1, jpoi) ; joi_m = MAX (joi-1, 1)
330	joj_p = MIN(joj+1, jpoj) ; joj_m = MAX (joj-1, 1)
331	IF ( iomskt (m1o (joi_p, joj )) == 1 &
332	.OR. iomskt (m1o (joi_m, joj )) == 1 &
333	.OR. iomskt (m1o (joi , joj_p)) == 1 &
334	.OR. iomskt (m1o (joi , joj_m)) == 1 &
335	.OR. iomskt (m1o (joi_p, joj_p)) == 1 &
336	.OR. iomskt (m1o (joi_p, joj_m)) == 1 &
337	.OR. iomskt (m1o (joi_m, joj_p)) == 1 &
338	.OR. iomskt (m1o (joi_m, joj_m)) == 1 ) THEN
339	locot (jo) = .TRUE.
340	ELSE
341	looce (jo) = .TRUE.
342	END IF
343	END IF
344	END DO
345	CALL lbc (locot, noperio, 'T', jpoi, jpoj)
346	CALL lbc (looce, noperio, 'T', jpoi, jpoj)
347
348	i2o = 0
349	WHERE (RESHAPE(locot,(/jpoi,jpoj/))) i2o=1
350	CALL flioputv (nid_id, "locot", i2o)
351	i2o = 0
352	WHERE (RESHAPE(looce,(/jpoi,jpoj/))) i2o=1
353	CALL flioputv (nid_id, "looce", i2o)
354
355	WRITE (nout,*) 'locot '
356	CALL prihbo (RESHAPE(locot, (/jpoi, jpoj/)), nout, l_norsud=.FALSE.)
357
358	!!
359	!! Determination d'une bande cotiere atm
360	la_nearcoast = lacot
361	ALLOCATE (la_temp(jpan), STAT=ierr) ; CALL chk_allo (ierr, 'la_temp')
362	DO niter = 1, 5
363	la_temp = la_nearcoast
364	DO ja = 1, jpan
365	IF ( iamskt(ja) == 1_il .AND. iamskp (ja) == 0_il ) THEN
366	jai = mai (ja) ; jaj = maj (ja) ! Indices 2D
367	jai_p = MIN(jai+1, jpai) ; jai_m = MAX (jai-1, 1)
368	jaj_p = MIN(jaj+1, jpaj) ; jaj_m = MAX (jaj-1, 1)
369	IF ( la_temp (m1a (jai_p, jaj )) &
370	.OR. la_temp (m1a (jai_m, jaj )) &
371	.OR. la_temp (m1a (jai , jaj_p)) &
372	.OR. la_temp (m1a (jai , jaj_m)) &
373	.OR. la_temp (m1a (jai_p, jaj_p)) &
374	.OR. la_temp (m1a (jai_p, jaj_m)) &
375	.OR. la_temp (m1a (jai_m, jaj_p)) &
376	.OR. la_temp (m1a (jai_m, jaj_m)) ) THEN
377	la_nearcoast (ja) = .TRUE.
378	END IF
379	END IF
380	END DO
381	CALL lbc (la_nearcoast, naperio, 'T', jpai, jpaj)
382	END DO
383
384	WRITE (nout,*) 'la_nearcoast '
385	CALL prihbo (RESHAPE(la_nearcoast, (/jpai, jpaj/)), nout, l_norsud=.TRUE.)
386	i2a = 0
387	WHERE (RESHAPE(la_nearcoast ,(/jpai,jpaj/))) i2a=1
388	CALL flioputv (nid_it, "la_nearcoast", i2a)
389
390	!!
391	!! Determination d'une bande cotiere ocean
392	lo_nearcoast = locot
393	ALLOCATE (lo_temp(jpon), STAT=ierr) ; CALL chk_allo (ierr, 'lo_temp')
394	DO niter = 1, 8
395	lo_temp = lo_nearcoast
396	DO jo = 1, jpon
397	IF ( iomskt(jo) == 0 .AND. iomskp (jo) == 0_il ) THEN
398	joi = moi (jo) ; joj = moj (jo) ! Indices 2D
399	joi_p = MIN(joi+1, jpoi) ; joi_m = MAX(joi-1, 1)
400	joj_p = MIN(joj+1, jpoj) ; joj_m = MAX(joj-1, 1)
401	IF ( lo_temp (m1o (joi_p, joj )) &
402	.OR. lo_temp (m1o (joi_m, joj )) &
403	.OR. lo_temp (m1o (joi , joj_p)) &
404	.OR. lo_temp (m1o (joi , joj_m)) &
405	.OR. lo_temp (m1o (joi_p, joj_p)) &
406	.OR. lo_temp (m1o (joi_p, joj_m)) &
407	.OR. lo_temp (m1o (joi_m, joj_p)) &
408	.OR. lo_temp (m1o (joi_m, joj_m)) ) THEN
409	lo_nearcoast (jo) = .TRUE.
410	END IF
411	END IF
412	END DO
413	CALL lbc (lo_nearcoast, noperio, 'T', jpoi, jpoj)
414	END DO
415
416	WRITE (nout,*) 'lo_nearcoast '
417	CALL prihbo (RESHAPE(lo_nearcoast, (/jpoi, jpoj/)), nout, l_norsud=.FALSE.)
418	i2o = 0
419	WHERE (RESHAPE(lo_nearcoast ,(/jpoi,jpoj/))) i2o=1
420	CALL flioputv (nid_id, "lo_nearcoast", i2o)
421
422	!!
423	WRITE (nout,*) 'Distance a la cote - points atmosphere'
424	DO ja = 1, jpan
425	dist_coast_a (ja) = rpi * ra
426	IF ( MOD(ja, 1000) .EQ. 0 ) WRITE (nout,*) ja, '/', jpan
427	IF ( iamskt (ja) .EQ. 0 ) THEN
428	dist_coast_a (ja) = 0.0E0_rl
429	CYCLE
430	END IF
431	IF ( .NOT. la_nearcoast (ja) ) THEN
432	dist_coast_a (ja) = 2.0_rl * dist_max_atm
433	CYCLE
434	END IF
435	DO jo = 1, jpon
436	IF ( locot(jo) ) THEN
437	dist_coast_a (ja) = MIN ( dist_coast_a (ja), &
438	& ra * geodist ( xalont(ja), xalatt(ja), xolont(jo), xolatt(jo) ) )
439	END IF
440	END DO
441	END DO
442	CALL lbc ( dist_coast_a, naperio, 'T', jpai, jpaj )
443	CALL prihre (RESHAPE(dist_coast_a, (/jpai,jpaj/)), 1.0E-4, nout, l_norsud=.TRUE.)
444
445	CALL flioputv (nid_it, "dist_coast_a", RESHAPE(dist_coast_a, (/jpai,jpaj/)))
446	!!
447	WRITE (nout,*) 'Distance a la cote - point ocean'
448	DO jo = 1, jpon
449	IF ( MOD(jo, 5000) .EQ. 0 ) WRITE (nout,*) jo, '/', jpon
450	dist_coast_o (jo) = rpi * ra
451	IF ( iomskt (jo) .EQ. 1 ) THEN
452	dist_coast_o (jo) = 0.0E0_rl
453	CYCLE
454	END IF
455	IF ( .NOT. lo_nearcoast (jo) ) THEN
456	dist_coast_o (jo) = 2.0_rl * dist_max_oce
457	CYCLE
458	END IF
459
460	DO jo2 = 1, jpon
461	IF ( locot(jo2) ) THEN
462	dist_coast_o (jo) = MIN ( dist_coast_o (jo), &
463	& ra * geodist ( xolont(jo), xolatt(jo), xolont(jo2), xolatt(jo2) ) )
464	END IF
465	END DO
466	END DO
467	CALL lbc ( dist_coast_a, naperio, 'T', jpoi, jpoj )
468	CALL prihre (RESHAPE(dist_coast_o, (/jpoi,jpoj/)), 1.0E-4, nout)
469
470	CALL flioputv (nid_id, "dist_coast_o", RESHAPE(dist_coast_o, (/jpoi,jpoj/)))
471
472	WRITE (nout,*) 'Verif coherence (tout doit etre nul)'
473	WRITE (nout,*) COUNT ( laland .AND. lacot)
474	WRITE (nout,*) COUNT ( lacot .AND. laoce)
475	WRITE (nout,*) COUNT ( laoce .AND. laland)
476	WRITE (nout,*) COUNT ( .NOT. ( laland .OR. lacot .OR. laoce) )
477
478
479
480	!!
481	WRITE (unit = nout, fmt = *) 'Nombre de points cotes atmosphere : ', COUNT (lacot)
482	WRITE (unit = nout, fmt = *) 'Nombre de points cotes ocean : ', COUNT (locot)
483	!!
484	CALL bilan ("Avec cote ", c_case="complet", c_int_atm="Integral", c_int_oce="Integral")
485	CALL bilan ("Avec cote ", c_case="ocean" , c_int_atm="Integral", c_int_oce="Integral")
486	CALL bilan ("Avec cote ", c_case="cote" , c_int_atm="Integral", c_int_oce="Integral")
487
488
489	E_DRYRUN: IF (l_dryrun) THEN
490	WRITE (nout,*) "Remplissage bidon des champs pour tester les I/O"
491	jma2o = jma2or
492	nvo = jma2or
493	WRITE (nout,*) 'jma2o pour test : ', jma2o
494	CALL RANDOM_SEED
495	CALL RANDOM_NUMBER (wa2o) ; ka2o = NINT (wa2o*REAL(jpan,rl) )
496
497	DO jo = 1, jpon
498	IF ( iomskt (jo) == 1 ) THEN
499	wa2o (:, jo) = 0.0_rl
500	ka2o (:, jo) = 0_il
501	END IF
502	END DO
503
504	DO jo = 1, jpon
505	DO jno = 1, jpa2o
506	ja = ka2o (jno, jo)
507	IF (ja > 0) THEN
508	IF (iamskt (ja) == 1) THEN
509	wa2o (jno, jo) = 0.0_rl ; ka2o (jno, jo) = 0_il
510	END IF
511	END IF
512	END DO
513	END DO
514	ELSE
515
516	CALL cotes_a
517
518	WRITE (nout, *) 'Renormalization'
519	!-$$ RenormA : DO ja = 1, jpan
520	!-$$ z_d_1 = 0.0e0_rl
521	!-$$ kom = 0
522	!-$$ RenormO : DO jo = 1, jpon
523	!-$$ DO jn = 1, jpa2o
524	!-$$ IF ( ka2o (jn,jo) == ja ) THEN
525	!-$$ z_d_1 = z_d_1 + wa2o (jn,jo)
526	!-$$ kom = kom + 1
527	!-$$ jx (kom,:) = (/ jn, jo /)
528	!-$$ !!$ WRITE (unit = nout, fmt = '(6I6,16F18.6) ') ja, jo, jn, kom, jx (kom, :), &
529	!-$$ !!$ wa2o (jn,jo), wa2o (jx (kom, 1), jx (kom, 2)), z_d_1, &
530	!-$$ !!$ xosrft (jo), xasrft (ja)
531	!-$$ END IF
532	!-$$ END DO
533	!-$$ ENDDO RenormO
534	!-$$ IF (kom /= 0 ) THEN
535	!-$$ IF ( z_d_1 == 0.0_rl ) THEN
536	!-$$ WRITE (nout, *) 'z_d_1 nul ', ja, xasrft (ja), lacot (ja), iamskt (ja), &
537	!-$$ & wa2o (:,ja), ka2o (:,ja)
538	!-$$ STOP
539	!-$$ ELSE
540	!-$$ DO ko = 1, kom
541	!-$$ wa2o (jx (ko,1), jx (ko,2)) = wa2o (jx (ko,1), jx (ko,2)) / z_d_1
542	!-$$ END DO
543	!-$$ ENDIF
544	!-$$ END IF
545	!-$$ ENDDO RenormA
546	wasum = 0.0_rl
547	DO jo = 1, jpon
548	IF ( nvo(jo) /= 0) THEN
549	DO jn = 1, nvo(jo)
550	ja = ka2o (jn,jo)
551	wasum (ja) = wasum (ja) + wa2o (jn,jo)
552	END DO
553	END IF
554	END DO
555	!!
556	DO jo = 1, jpon
557	IF ( nvo(jo) /= 0 ) THEN
558	DO jn = 1, nvo(jo)
559	ja = ka2o (jn,jo)
560	IF ( wasum (ja) /= 0.0_rl ) wa2o (jn,jo) = wa2o (jn,jo) / wasum (ja)
561	END DO
562	END IF
563	END DO
564
565	!!
566	CALL int_wei
567	!!
568	!! Controle
569	!!
570	CALL verif (' Normalisation')
571	!!
572	CALL bilan ("Final ", c_case="terre" , c_int_atm="Integral", c_int_oce="Local")
573	CALL bilan ("Final ", c_case="terre100", c_int_atm="Integral", c_int_oce="Local")
574	CALL bilan ("Final ", c_case="cote" , c_int_atm="Integral", c_int_oce="Local")
575	!!
576	WRITE (unit = nout, fmt = * ) 'Nombre de voisins : ', jma2o
577	!!
578	CALL gather_wei ('apres normalisation')
579	CALL verif ('apres gather' )
580
581	!! Atm mask interpolated to oce
582	za = REAL (1_il-iamskt, KIND = rl )
583	CALL inter_a2o (za, a2omask )
584	!! 1 on Atm interpolated to oce
585	za = REAL (1_il , KIND = rl )
586	CALL inter_a2o (za, a2ofull )
587
588	!!
589	!! Checking weights
590	WRITE (unit = nout, fmt = '("Neighbors a2o : ", 3I9 )' ) &
591	& MINVAL (nvo, nvo /=0), MAXVAL (nvo), jma2o
592	WRITE (unit = nout, fmt = '("Index a2o : ", 2I9 )' ) &
593	& MINVAL (ka2o, ka2o /= 0), MAXVAL (ka2o)
594	jind = MINLOC(wa2o, wa2o /= 0.0_rl) ; jo = jind(2) ; jn = jind(1)
595	joi = m2oi(jo) ; joj = m2oj(jo)
596	WRITE (unit = nout, fmt = '("Weights a2o mini : ", 4I6, 1E13.4, 3I6)' ) &
597	& jo, jn, joi, joj, wa2o(jn,jo), ka2o(jn,jo), m2ai(ka2o(jn,jo)), m2aj(ka2o(jn,jo))
598	jind = MAXLOC(wa2o, wa2o /= 0.0_rl) ; jo = jind(2) ; jn = jind(1)
599	joi = m2oi(jo) ; joj = m2oj(jo)
600	WRITE (unit = nout, fmt = '("Weights a2o maxi : ", 4I6, 1E13.4, 3I6)' ) &
601	& jo, jn, joi, joj, wa2o(jn,jo), ka2o(jn,jo), m2ai(ka2o(jn,jo)), m2aj(ka2o(jn,jo))
602	!
603	jind(1:1) = INT(MAXLOC (nvo),il) ; jo = jind(1)
604	WRITE (unit = nout, fmt = *) jo, m2oi(jo), m2oj(jo), ka2o(1:nvo(jo),jo), &
605	& (m2ai(ka2o(jn,jo)), m2aj(ka2o(jn,jo)), jn = 1, nvo(jo))
606	!!
607	clweight = "WEIGHTS4" ; cladress = "ADRESSE4"
608	WRITE (unit = nout, fmt = *) cladress, " ", clweight, " ", jpa2o, jma2o
609	!!
610	!! Controle
611	!!
612	WRITE (nout, *) 'Verif avant ecriture'
613	WRITE (nout, *) 'wa2o ', MINVAL (wa2o), MAXVAL (wa2o), COUNT (wa2o > 0.0_rl)
614	WRITE (nout, *) 'ka2o ', MINVAL (ka2o), MAXVAL (ka2o), COUNT (ka2o > 0 )
615	WRITE (nout, *) 'xasrft ', MINVAL (xasrft), MAXVAL (xasrft)
616	WRITE (nout, *) 'o2amask ', MINVAL (o2amask), MAXVAL (o2amask), SUM (o2amask)
617	WRITE (nout, *) 'xosrft ', MINVAL (xosrft), MAXVAL (xosrft)
618	WRITE (nout, *) 'iomskt ', MINVAL (iomskt), MAXVAL (iomskt), COUNT (iomskt == 0)
619	WRITE (nout, *) 'iomskp ', MINVAL (iomskp), MAXVAL (iomskp), COUNT (iomskp == 0)
620	WRITE (nout, *) 'iamskt ', MINVAL (iamskt), MAXVAL (iamskt), COUNT (iamskt == 0)
621	WRITE (nout, *) 'locot ', COUNT (locot)
622	!!
623	!!
624	END IF E_DRYRUN
625	!!
626	!!
627	WRITE (nout, *) ' '
628	WRITE (nout, *) ' '
629	WRITE (nout, *) 'Eventuels traitement specifiques'
630	WRITE (nout, *) ' '
631	WRITE (nout, *) ' '
632
633
634	!! Extension ??
635
636	!! Traitements specifiques, tres specifiques ... !!
637	IF ( TRIM(c_suffix) == "_runoffNord09" ) THEN
638	WRITE (nout, *) 'Traitement specifique _runoffNord09'
639	DO jo = 1, jpon
640	IF ( tmask_atl(jo) == 1 .AND. xolatt(jo) > 40.0_rl ) wa2o(:,jo) = wa2o(:,jo) * 0.9
641	END DO
642	END IF
643	!
644	IF ( TRIM(c_suffix) == "_runoffNord07" ) THEN
645	WRITE (nout, *) 'Traitement specifique _runoffNord07'
646	DO jo = 1, jpon
647	IF ( tmask_atl(jo) == 1 .AND. xolatt(jo) > 40.0_rl ) wa2o(:,jo) = wa2o(:,jo) * 0.7
648	END DO
649	END IF
650
651	IF ( TRIM(c_suffix) == "_runoffNord00" ) THEN
652	WRITE (nout, *) 'Traitement specifique _runoffNord00'
653	DO jo = 1, jpon
654	IF ( tmask_atl(jo) == 1 .AND. xolatt(jo) > 40.0_rl ) wa2o(:,jo) = wa2o(:,jo) * 0.0
655	END DO
656	END IF
657	!!
658	!!
659	!! Verif nombre de voisins
660	!!
661	IF ( jma2or /= jma2o ) THEN
662	WRITE (nout, *) 'Erreur nombre de voisins dans run.def : '
663	WRITE (nout, *) ' jma2o calcule / lu : ', jma2o, jma2or
664	!STOP
665	ENDIF
666	!!
667	!-$$ !! Ecriture des poids
668	!-$$ !!
669	!-$$ WRITE (unit = c_r, fmt = '(1i1)' ) rk_out
670	!-$$ !
671	!-$$ IF ( l_wei_i4) THEN
672	!-$$ WRITE (unit = c_i, fmt = '(1i1)' ) i_4
673	!-$$ cfname4 = TRIM(mo_name) // ".wa2o.runoff.i" // TRIM(c_i) // ".r" // TRIM(c_r) // TRIM(c_suffix)
674	!-$$ OPEN (unit = nwei4a2o, file = TRIM(cfname4), form = 'unformatted', status = 'replace',&
675	!-$$ & action = 'write')
676	!-$$ WRITE (nout, *) 'Poids i_4 a -> o: ', TRIM(cfname4)
677	!-$$ END IF
678	!-$$ !
679	!-$$ IF (l_wei_i8) THEN
680	!-$$ WRITE (unit = c_i, fmt = '(1i1)' ) i_8
681	!-$$ cfname8 = TRIM(mo_name) // ".wa2o.runoff.i" // TRIM(c_i) // ".r" // TRIM(c_r) // TRIM(c_suffix)
682	!-$$ OPEN (unit = nwei8a2o, file = TRIM(cfname8), form = 'unformatted', status = 'replace',&
683	!-$$ & action = 'write')
684	!-$$ WRITE (nout, *) 'Poids i_8 a -> o : ', TRIM(cfname8)
685	!-$$ END IF
686	!-$$ !
687	!-$$ cldiag_a2o = TRIM(a2o_name) // '.runcoa.diag'
688	!-$$ clw_a2o = TRIM(mo_name) // '.wa2o.runcoa'
689	!-$$ clw_a2o_mct= TRIM(mo_name) // '.wa2o_mct.runcoa'
690	!-$$ !
691	!-$$ CALL wri_weights_a2o (cldiag_a2o, clw_a2o, clw_a2o_mct, "4", l_fulldiag, lo_src_grid_frac=.FALSE., lo_dst_grid_frac=.FALSE., &
692	!-$$ & co_omsk=TRIM(cotes_omsk), co_amsk=TRIM(cotes_amsk) )
693	!-$$ !
694
695
696	!!
697	!! Ecriture des poids
698	!!
699	!! Ouverture des fichiers poids
700	!!
701	WRITE (unit = c_r, fmt = '(1i1)' ) rk_out
702	!
703	IF (l_wei_i4) THEN
704	WRITE (unit = c_i, fmt = '(1i1)' ) i_4
705	cfname4 = TRIM(mo_name) // ".wa2o.rivflu" // TRIM(c_suffix) // ".i" // TRIM(c_i) // ".r" // TRIM(c_r)
706	OPEN (unit = nwei4a2o, file = TRIM(cfname4), form = 'unformatted', status = 'replace',&
707	& action = 'write')
708	WRITE (nout, *) 'Poids i4 a -> o: ', TRIM(cfname4)
709	END IF
710	!
711	IF (l_wei_i4) THEN
712	WRITE (unit = c_i, fmt = '(1i1)' ) i_8
713	cfname8 = TRIM (mo_name) // ".wa2o.rivflu" // TRIM(c_suffix) // ".i" // TRIM(c_i) // ".r" // TRIM(c_r)
714	OPEN (unit = nwei8a2o, file = TRIM(cfname8), form = 'unformatted', status = 'replace',&
715	& action = 'write')
716	WRITE (nout, *) 'Poids i8 a -> o : ', TRIM(cfname8)
717	END IF
718
719	!
720	cldiag_a2o = TRIM (a2o_name) // '.rivflu.diag'
721	clw_a2o = TRIM (mo_name) // '.wa2o.rivflu'
722	clw_a2o_mct = "rmp_"//TRIM(camod_t)//"_to_"//TRIM(comod_t)//"_MOSAIC_rivflu"
723	!
724	CALL wri_weights_a2o (cldiag_a2o, clw_a2o, clw_a2o_mct, "3", l_fulldiag, lo_src_grid_frac=.FALSE., lo_dst_grid_frac=.FALSE., &
725	& co_omsk='noperio', co_amsk='full' )
726
727	!!
728	IF (l_wei_i4) THEN
729	WRITE (nout,*) 'Fermeture poids i4'
730	CLOSE (unit = nwei8a2o)
731	END IF
732	IF (l_wei_i8) THEN
733	WRITE (nout,*) 'Fermeture poids i8 '
734	CLOSE (unit = nwei4a2o)
735	END IF
736	!!
737	WRITE (nout,*) 'Appel flioclo'
738	CALL flioclo
739	!!
740	STOP
741	!!
742	CONTAINS
743	!!
744	SUBROUTINE diag_itarget (clmess, lwrite)
745	CHARACTER (len=*), INTENT (in) :: clmess
746	LOGICAL, INTENT (in) :: lwrite
747	INTEGER, SAVE :: nappel = 0
748	CHARACTER (len=2) :: cc
749	!!
750	nappel = nappel + 1
751	!!
752	itarget = 0_il ; wtarget = 0.0_rl
753	DO jo = 1, jpon
754	DO jn = 1, jpa2o
755	ja = ka2o (jn, jo)
756	IF (ja /= 0 ) THEN
757	itarget (ja) = itarget (ja) + 1_il
758	wtarget (ja) = wtarget (ja) + wa2o (jn, jo)
759	ENDIF
760	END DO
761	END DO
762	!!
763	!CALL lbc (itarget, naperio, 'T', jpai, jpaj)
764	!CALL lbc (wtarget, naperio, 'T', jpai, jpaj)
765	!!
766	IF (lwrite) THEN
767	!$$$ WRITE(nout,*) 'itarget, etape '// clmess
768	!$$$ CALL prihin (RESHAPE(itarget,(/jpai,jpaj/)),1_il,nout)
769	WRITE (unit = cc, fmt = '(1I2.2)' ) nappel
770	CALL fliodefv (nid_it, 'itarget' // cc, (/1, 2/), v_t=flio_i4)
771	CALL fliodefv (nid_it, 'ntarget' // cc, (/1, 2/), v_t=flio_i4)
772	CALL fliodefv (nid_it, 'wtarget' // cc, (/1, 2/), v_t=flio_r4)
773	CALL fliodefv (nid_it, 'non_target' // cc, (/1, 2/), v_t=flio_i4)
774
775
776	CALL flioputa (nid_it, 'itarget' // cc, 'long_name ', clmess // ' itarget' // cc)
777	CALL flioputv (nid_it, 'itarget' // cc, MIN (1, RESHAPE (itarget, (/jpai, jpaj/))) )
778
779	CALL flioputa (nid_it, 'ntarget' // cc, 'long_name ', clmess // ' ntarget' // cc)
780	CALL flioputv (nid_it, 'ntarget' // cc, RESHAPE (itarget, (/jpai, jpaj/)) )
781
782	CALL flioputa (nid_it, 'wtarget' // cc, 'long_name ', clmess // ' wtarget' // cc)
783	CALL flioputv (nid_it, 'wtarget' // cc, RESHAPE (wtarget, (/jpai, jpaj/)) )
784
785	wtarget = 0
786	WHERE (itarget .EQ. 0 ) wtarget = 1.0_rl
787	CALL flioputa (nid_it, 'non_target' // cc, 'long_name ', clmess // ' non_target' // cc)
788	CALL flioputv (nid_it, 'non_target' //cc, RESHAPE (NINT (wtarget), (/jpai, jpaj/)) )
789	END IF
790	!!
791	END SUBROUTINE diag_itarget
792	!!
793	SUBROUTINE diag_idest (clmess, lwrite)
794	CHARACTER (len=*), INTENT (in) :: clmess
795	LOGICAL, INTENT (in) :: lwrite
796	INTEGER, SAVE :: nappel = 0
797	CHARACTER (len=2) :: cc
798	!!
799	nappel = nappel + 1
800	!!
801	idest = 0_il ; wdest = 0.0_rl
802	DO jo = 1, jpon
803	DO jn = 1, jpa2o
804	ja = ka2o (jn, jo)
805	IF (ja /= 0 ) THEN
806	idest (jo) = idest (jo) + 1_il
807	wdest (jo) = wdest (jo) + wa2o (jn, jo)
808	ENDIF
809	END DO
810	END DO
811	!!
812	!CALL lbc (itarget, naperio, 'T', jpai, jpaj)
813	!CALL lbc (wtarget, naperio, 'T', jpai, jpaj)
814	!!
815	IF (lwrite) THEN
816	!$$$ WRITE(nout,*) 'idest, etape '// clmess
817	!$$$ CALL prihin (RESHAPE(idest,(/jpoi,jpoj/)),1_il,nout)
818	WRITE (unit = cc, fmt = '(1I2.2)' ) nappel
819	CALL fliodefv (nid_id, 'idest' // cc, (/1, 2/), v_t=flio_i4)
820	CALL fliodefv (nid_id, 'ndest' // cc, (/1, 2/), v_t=flio_i4)
821	CALL fliodefv (nid_id, 'wdest' // cc, (/1, 2/), v_t=flio_r4)
822	CALL fliodefv (nid_id, 'non_dest' // cc, (/1, 2/), v_t=flio_i4)
823
824
825	CALL flioputa (nid_id, 'idest' // cc, 'long_name ', clmess // ' idest' // cc)
826	CALL flioputv (nid_id, 'idest' // cc, MIN (1, RESHAPE (idest, (/jpoi, jpoj/))) )
827
828	CALL flioputa (nid_id, 'ndest' // cc, 'long_name ', clmess // ' ndest' // cc)
829	CALL flioputv (nid_id, 'ndest' // cc, RESHAPE (idest, (/jpoi, jpoj/)) )
830
831	CALL flioputa (nid_id, 'wdest' // cc, 'long_name ', clmess // ' wdest' // cc)
832	CALL flioputv (nid_id, 'wdest' // cc, RESHAPE (wdest, (/jpoi, jpoj/)) )
833
834	wdest = 0
835	WHERE (idest .EQ. 0 ) wdest = 1.0_rl
836	CALL flioputa (nid_id, 'non_dest' // cc, 'long_name ', clmess // ' non_dest' // cc)
837	CALL flioputv (nid_id, 'non_dest' //cc, RESHAPE (NINT (wdest), (/jpoi, jpoj/)) )
838	END IF
839	!!
840	END SUBROUTINE diag_idest
841	!!
842	!!
843	SUBROUTINE int_wei
844	!!
845	IMPLICIT NONE
846	WRITE(nout,*) 'Debut normalisation'
847	!!
848	IF (.NOT. lint_atm) THEN
849	WRITE(nout,*) 'Multiplication par les surfaces atm'
850	DO jo = 1, jpon
851	DO jn = 1, jma2o
852	ja = ka2o (jn, jo)
853	IF (ja > 0 ) wa2o (jn,jo) = wa2o (jn,jo) * xasrft (ja)
854	END DO
855	ENDDO
856	ENDIF
857	!!
858	IF (.NOT. lint_oce) THEN
859	WRITE (nout,*) 'Division par les surfaces oce'
860	DO jo = 1, jpon
861	DO jn = 1, jma2o
862	wa2o (jn,jo) = wa2o (jn,jo) / xosrft (jo)
863	END DO
864	ENDDO
865	ENDIF
866	END SUBROUTINE int_wei
867	!!
868	SUBROUTINE complet_run (ka1, ka2, kn)
869	INTEGER (kind=il), INTENT(in) :: ka1, ka2 ! Index of atm boxes
870	INTEGER (kind=il), INTENT(in) :: kn ! Number of neighbors
871	INTEGER (kind=il) :: jt
872	!!
873	DO jo = 1, jpon
874	IF (iomskt(jo) == 1 ) CYCLE
875	IF (iomskp(jo) .EQ. 1_il ) CYCLE
876	DO jn = 1, jpa2o
877	IF (ka2o(jn,jo) == ka1 ) THEN
878	nvo(jo) = nvo(jo)+1
879	IF (nvo(jo) > jpa2o ) THEN
880	WRITE(nout,*) 'Pas assez de voisins pour la completion '
881	WRITE(nout,fmt='(1I7,2I4,2F10.1)' ) ka1, m2ai(ka1), m2aj(ka1), &
882	& xalont(ka1), xalatt(ka1)
883	WRITE(nout,fmt='(1I7,2I4,2F10.1)' ) ka2, m2ai(ka2), m2aj(ka2), &
884	& xalont(ka2), xalatt(ka2)
885	WRITE(nout,fmt='(1I7,2I4,2F10.1,2I8)') jo, m2oi(jo), m2oj(jo), &
886	& xolont(jo), xolatt(jo), jn, nvo(jo)
887	DO jt = 1, jpa2o
888	WRITE(nout,fmt=*) jt, ka2o(jt,jo), xalont(ka2o(jt,jo)), xalatt(ka2o(jt,jo))
889	ENDDO
890	END IF
891	ka2o (nvo (jo), jo) = ka2
892	wa2o (nvo (jo), jo) = wa2o(jn,jo) * xasrft(ka2) / xasrft(ka1) / REAL(kn,KIND=rl)
893	END IF
894	END DO
895	ENDDO
896	END SUBROUTINE complet_run
897	!!
898	SUBROUTINE gather_wei (clmess)
899	CHARACTER (len=*), INTENT (in) :: clmess
900	INTEGER (kind=il) :: jn1, jn2, jo2
901	!! Eliminate redundancy
902	DO jo = 1, jpon
903	DO jn1 = 1, jpa2o-1
904	IF (ka2o (jn1,jo) == 0 ) CYCLE
905	DO jn2 = jn1+1, jpa2o
906	IF (ka2o (jn1, jo) == ka2o (jn2, jo) ) THEN
907	WRITE(nout,*) 'Redondance ', jo, jn1, jn2, ka2o(jn1,jo), ka2o (jn2,jo), wa2o(jn1,jo), wa2o(jn2,jo)
908	wa2o (jn1, jo) = wa2o (jn1, jo) + wa2o (jn2, jo)
909	ka2o (jn2, jo) = 0
910	END IF
911	END DO
912	END DO
913	END DO
914	!!
915	!! Seek for strangeness
916	DO jo = 1, jpon
917	DO jn = 1, jpa2o
918	ja = ka2o (jn, jo)
919	IF (ja == 0) CYCLE
920	IF (wa2o (jn, jo) == 0.0E0 ) THEN
921	WRITE (nout, *) 'Null weight : ', jo, jn, ja, locot (jo)
922	END IF
923	END DO
924	END DO
925	!!
926	!! Periodicity
927	!DO jn = 1, jpa2o
928	! CALL lbc (wa2o (jn,:), noperio, 'T', jpoi, jpoj)
929	! CALL lbc (ka2o (jn,:), noperio, 'T', jpoi, jpoj)
930	!END DO
931	CALL diag_itarget (' gather/perio - ' // clmess, lwrite=.TRUE.)
932	CALL diag_idest (' gather/perio - ' // clmess, lwrite=.TRUE.)
933	!!
934	!! Gather weights
935	DO jo = 1, jpon
936	ktemp = ka2o (:,jo) ; wtemp = wa2o (:,jo)
937	ka2o (:,jo) = 0_il ; wa2o (:,jo) = 0.0_rl
938	kn = 1_il
939	DO jo2 = 1, jpa2o
940	IF (ktemp (jo2) /= 0_il ) THEN
941	ka2o (kn, jo) = ktemp (jo2)
942	wa2o (kn, jo) = wtemp (jo2)
943	kn = kn + 1_il
944	END IF
945	END DO
946	END DO
947	nvo = COUNT (ka2o /= 0_il, DIM = 1)
948	jma2o = MAXVAL (nvo)
949	!
950	!!
951	!! Periodicity
952	!DO jn = 1, jpa2o
953	! CALL lbc (wa2o(jn,:), noperio, 'T', jpoi, jpoj)
954	! CALL lbc (ka2o(jn,:), noperio, 'T', jpoi, jpoj)
955	!END DO
956	!!
957	nvo = COUNT (ka2o /= 0_il, DIM = 1)
958	jma2o = MAXVAL (nvo)
959	!!
960	END SUBROUTINE gather_wei
961	!!
962	SUBROUTINE verif (clmess, l_non_a)
963	!!
964	CHARACTER (len=*), INTENT (in) :: clmess
965	INTEGER :: ja
966	LOGICAL, INTENT (in), OPTIONAL :: l_non_a
967	!!
968	nvo = COUNT (ka2o /= 0_il, DIM = 1)
969	jma2o = MAXVAL (nvo)
970	CALL diag_itarget (clmess, lwrite=.TRUE.)
971	CALL diag_idest (clmess, lwrite=.TRUE.)
972	!!
973	WRITE (nout, *) ' --- Verifs -- : ', TRIM (clmess)
974	WRITE (nout, *) 'wa2o ', MINVAL (wa2o), MAXVAL (wa2o), COUNT (wa2o > 0.0_rl)
975	WRITE (nout, *) 'ka2o ', MINVAL (ka2o), MAXVAL (ka2o), COUNT (ka2o > 0 )
976	WRITE (nout, *) 'xasrft ', MINVAL (xasrft), MAXVAL (xasrft)
977	WRITE (nout, *) 'o2amask ', MINVAL (o2amask), MAXVAL (o2amask), SUM (o2amask)
978	WRITE (nout, *) 'xosrft ', MINVAL (xosrft), MAXVAL (xosrft)
979	WRITE (nout, *) 'iomskt ', MINVAL (iomskt), MAXVAL (iomskt), COUNT (iomskt == 0)
980	WRITE (nout, *) 'iomskp ', MINVAL (iomskp), MAXVAL (iomskp), COUNT (iomskp == 0)
981	WRITE (nout, *) 'iamskt ', MINVAL (iamskt), MAXVAL (iamskt), COUNT (iamskt == 0)
982	WRITE (nout, *) 'locot ', COUNT (locot)
983	WRITE (nout, *) 'Points non attribues : ', COUNT (itarget .EQ. 0 )
984
985	IF (PRESENT (l_non_a)) THEN
986	DO ja = 1, jpan
987	IF (itarget (ja) == 0 .AND. iamskp(ja) == 0 .AND. (laland (ja) .OR. lacot(ja) )) THEN
988	!-$$ WRITE (UNIT=nout,FMT='("Point atm non attribues : ", 3I7, 2F7.1, 1I4)' ) ja, m2ai(ja), m2aj (ja), xalont(ja), xalatt(ja), itarget(ja)
989	WRITE (nout,*) ja, m2ai(ja), m2aj (ja), xalont(ja), xalatt(ja), itarget(ja), laland(ja), lacot(ja), laoce(ja)
990	ENDIF
991	END DO
992	END IF
993	END SUBROUTINE verif
994	!!
995	!!
996	SUBROUTINE noroute
997	!! Determine les points atmospheres non routes
998	lanoroute = .TRUE.
999
1000	DO jn = 1, jpa2o
1001	DO jo = 1, jpon
1002	ja = ka2o ( jo, jn)
1003	IF ( ja /= 0 .AND. wa2o (jo,jn) > 0.0_rl ) THEN
1004	lanoroute = .FALSE.
1005	END IF
1006	END DO
1007	END DO
1008
1009	END SUBROUTINE noroute
1010	!!
1011	SUBROUTINE bilan (clmess, c_case, c_int_atm, c_int_oce)
1012	!!
1013	CHARACTER (len=*), INTENT (in) :: clmess
1014	CHARACTER (len=*), INTENT (in) :: c_case ! Type de champ d'entree
1015	CHARACTER (len=*), INTENT (in) :: c_int_atm ! Type d'integrale sur l'atm
1016	CHARACTER (len=*), INTENT (in) :: c_int_oce ! Type d'integrale sur l'ocean
1017	!!
1018	WRITE (nout,*) "Bilan :" // TRIM (clmess) // ", cas :", TRIM (c_case), ", c_int_atm :", TRIM (c_int_atm),&
1019	& ", c_int_oce :", TRIM(c_int_oce), "."
1020	!!
1021	za = 0.0_rl
1022	!!
1023	SELECT CASE (TRIM (c_case))
1024	CASE ('complet')
1025	za = 1.0_rl * REAL (1-iamskp, KIND=rl)
1026	SELECT CASE (TRIM (c_int_atm))
1027	CASE ('Integral')
1028	za = za * xasrft
1029	zsuma = SUM (za)
1030	CASE ('Local')
1031	zsuma = DOT_PRODUCT (za, xasrft)
1032	CASE default
1033	WRITE (nout,*) 'Erreur type de bilan atm' ; STOP
1034	END SELECT
1035	!!
1036	CASE ('terre')
1037	WHERE (laland .OR. lacot)
1038	za = 1.0_rl * REAL (1-iamskp, KIND=rl)
1039	END WHERE
1040	SELECT CASE (TRIM (c_int_atm))
1041	CASE ('Integral')
1042	za = za * xasrft
1043	zsuma = SUM (za)
1044	CASE ('Local')
1045	zsuma = DOT_PRODUCT (za, xasrft*(1.0-o2amask))
1046	CASE default
1047	WRITE (nout,*) 'Erreur type de bilan atm' ; STOP
1048	END SELECT
1049	!!
1050	CASE ('terre100')
1051	WHERE (laland)
1052	za = 1.0_rl * REAL (1-iamskp, KIND=rl)
1053	END WHERE
1054	SELECT CASE (TRIM (c_int_atm))
1055	CASE ('Integral')
1056	za = za * xasrft
1057	zsuma = SUM (za)
1058	CASE ('Local')
1059	zsuma = DOT_PRODUCT (za, xasrft)
1060	CASE default
1061	WRITE (nout,*) 'Erreur type de bilan atm' ; STOP
1062	END SELECT
1063	!!
1064	CASE ('ocean')
1065	WHERE (lacot)
1066	za = 1.0_rl * REAL (1-iamskp, KIND=rl)
1067	END WHERE
1068	SELECT CASE (TRIM (c_int_atm))
1069	CASE ('Integral')
1070	za = za * xasrft
1071	zsuma = SUM (za)
1072	CASE ('Local')
1073	zsuma = DOT_PRODUCT (za, xasrft*o2amask)
1074	CASE default
1075	WRITE (nout,*) TRIM(clmess), ' : Erreur type de bilan atm' ; STOP
1076	END SELECT
1077	!!
1078	CASE ('cote')
1079	WHERE (lacot)
1080	za = 1.0_rl * REAL (1-iamskp, KIND=rl)
1081	END WHERE
1082	SELECT CASE (TRIM (c_int_atm))
1083	CASE ('Integral')
1084	za = za * xasrft
1085	zsuma = SUM (za)
1086	CASE ('Local')
1087	zsuma = DOT_PRODUCT (za, xasrft*o2amask)
1088	CASE default
1089	WRITE (nout,*) 'Erreur type de bilan atm' ; STOP
1090	END SELECT
1091	!!
1092	CASE default
1093	WRITE (nout,*)'Erreur choix de cas' ; STOP
1094	STOP
1095	END SELECT
1096	!!
1097	CALL inter_a2o (za, zo)
1098	!!
1099	SELECT CASE (TRIM (c_int_oce))
1100	CASE ('Integral')
1101	zsumo = DOT_PRODUCT (zo, REAL(1-iomskp, KIND=rl))
1102	CASE ('Local')
1103	zsumo = DOT_PRODUCT (zo*xosrft,REAL(1-iomskp, KIND=rl))
1104	CASE default
1105	WRITE (nout,*) TRIM (clmess), ' : Erreur type de bilan ' ; STOP
1106	END SELECT
1107	!!
1108	WRITE (nout, fmt='(" bilan atm: ", 1E15.6, ", bilan oce: ", 1E15.6)' ) zsuma, zsumo
1109	!!
1110	RETURN
1111	!!
1112	END SUBROUTINE bilan
1113	!!
1114	END PROGRAM cotes_etal
1115

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

Context Navigation

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

Download in other formats: