Context Navigation

source: TOOLS/MOZAIC/src/MOZAIC/mosaic.f90 @ 3329

Last change on this file since 3329 was 3329, checked in by omamce, 7 years ago
OM: essai propset
Property svn:keywords set to `Date Revision HeadURL Author`
File size: 27.2 KB

Line
1	! -- Mode: f90 --
2	PROGRAM mosaic
3	! -------------------------------------------------------------------------------
4	!
5	!> @brief Generates interpolation weights and adresses for Oasis Coupler, Mozaic interpolation
6	!
7	!>
8	!> Oceanic model : OPA (C Grid)
9	!> Atmospheric Model : LMD5, LMDz or spectral model with Gaussian grid
10	!>
11	!> SVN :
12	!> $Author$
13	!> $Date$
14	!> $Revision$
15	!> $Id$
16	!>
17	!> $Log$
18	!>
19	!>
20	!> -------------------------------------------------------------------------------
21	!>
22	!> http://www.ipsl.jussieu.fr/~omamce/IPSLCM4/MOZAIC
23	!>
24	!> -------------------------------------------------------------------------------
25	!>
26	!> Warning, to install, configure, run, use any of Olivier Marti's
27	!> software or to read the associated documentation you'll need at least
28	!> one (1) brain in a reasonably working order. Lack of this implement
29	!> will void any warranties (either express or implied).
30	!> O. Marti assumes no responsability for errors, omissions,
31	!> data loss, or any other consequences caused directly or indirectly by
32	!> the usage of his software by incorrectly or partially configured
33	!> personal.
34	!>
35	!> 1) Read mask, grid and areas files, at OASIS format
36	!> Convention for mask are those of Oasis : 0 for sea, 1 for land
37	!>
38	!> 2) Generates interpolation weights
39	!>
40	!> 3) Write them at Oasis format
41	!>
42	!> 4) Write some diagnostic fields in NetCDF format
43	!>
44	!>
45	!> @author Olivier MARTI, june 1996
46	!
47	!> Modification : Olivier MARTI, may 1997 : LMDZ case, from L. Li
48	!> Eric GUILYARDI, may 1999 : Gaussian Grid + ORCA
49	!> Olivier MARTI and Jacques BELLIER, april 2000 : use of polygon matching to compute weights
50	!>
51	! --------------------------------------------------------------------------------
52	! LSCE, CEA-Saclay
53	! --------------------------------------------------------------------------------
54	!
55	!
56	USE modeles
57	USE mod_inipar
58	USE mod_locio
59	USE mod_neuf
60	USE bords
61	USE mod_lbc
62	USE mod_wripoly
63	USE mod_polyw_a
64	USE mod_polyw_o
65	USE formula
66	USE mod_lecdata
67	USE mod_prih
68	USE mod_norma
69	USE mod_rectifie
70	USE math
71	USE mod_inter
72	USE mod_wri_wei
73	USE fliocom
74	USE getincom
75	USE errioipsl
76	USE mod_polom
77	USE mod_proj
78	USE fliocom
79	!
80	IMPLICIT NONE
81	!
82	REAL (kind=rl) :: zzz
83	REAL (kind=rd) :: zs, zs1, zs2
84	INTEGER (kind=il) :: ja, jo, jai, jaj, joi, joj, jno, jna, kn ! Dummy indexes
85	INTEGER (kind=il), DIMENSION (2) :: jind
86	CHARACTER (len = 80) :: cfname4, cfname8
87	!
88	INTEGER, DIMENSION (:,:), ALLOCATABLE :: mt1, mt2, mt3
89	!
90	! NetCDF
91	CHARACTER (len = 50) :: cldiag_o2a ! NetCDF file name for fields on atmospheric grid (diagnostics)
92	CHARACTER (len = 50) :: cldiag_a2o ! NetCDF file name for fields on oceanic grid (diagnostics)
93	CHARACTER (len = 50) :: clw_o2a, clw_a2o, clw_a2o_mct, clw_o2a_mct
94	CHARACTER (len = 25) :: mo_name, o2a_name, a2o_name
95	CHARACTER (len = 1) :: c_i, c_r
96	INTEGER :: il_mask
97	CHARACTER (len=10) :: cmodel
98	!
99	REAL (kind=rd), DIMENSION (:), ALLOCATABLE :: xa, ya, xo, yo
100	INTEGER (kind=il) :: ierr
101	INTEGER (kind=il), DIMENSION (:,:), ALLOCATABLE :: itmp
102	CHARACTER (len=80) :: clarg
103	INTEGER :: narg
104	!>
105	!! Lecture de la ligne de comande
106	narg = iargc()
107	IF ( narg > 0 ) THEN
108	CALL getarg ( 1, clarg)
109	SELECT CASE (TRIM(clarg))
110	CASE Default
111	IF ( INDEX ( TRIM(clarg), '.def' ) /= 0 ) THEN
112	WRITE (nout,*) 'Lecture des parametres dans ', TRIM(clarg)
113	CALL getin_name ( TRIM (clarg) )
114	ELSE
115	WRITE (nout,*) 'Lecture des parametres dans run.def'
116	ENDIF
117	END SELECT
118	END IF
119	!
120	CALL inipar
121	CALL alloc_modeles
122	!
123	!> @par Initialisation
124	!
125	ngrd = 11 ; nsrf = 12 ; nmsk = 13; nchk = 16 ; ndeb = 9 ; nbug1 = 10
126	nwei4o2a = 21 ; nwei8o2a = 22 ; nwei4a2o = 23 ; nwei8a2o = 24
127	! Diagnostics files
128	!
129	nout = 6
130	IF (nout /= 6 .AND. nout /= 0) THEN
131	OPEN (unit = nout, file = 'mosaic' // TRIM(c_suffix) // '.out', position = 'REWIND', action = 'WRITE', status = 'replace')
132	END IF
133	!>
134
135	WRITE (nout,*) 'm2ai ', MINVAL (m2ai), MAXVAL (m2ai)
136	WRITE (nout,*) 'm2aj ', MINVAL (m2aj), MAXVAL (m2aj)
137	!
138	ndebug = 2
139	!
140	IF (lwra2o) THEN
141	WRITE (unit = nout, fmt = *) ' Compute weights and adresses for interpolation atm -> oce '
142	ELSE
143	WRITE (unit = nout, fmt = *) ' Nothing for atm -> oce '
144	ENDIF
145	IF (lwro2a) THEN
146	WRITE (unit = nout, fmt = *) ' Compute weights and adresses for interpolation oce -> atm '
147	ELSE
148	WRITE (unit = nout, fmt = *) ' Nothing for oce -> atm '
149	ENDIF
150	!
151	WRITE (nout,*) 'epsfrac : ', epsfrac
152	!
153	ka2o = 0 ; wa2o = 0.0_rl ; ko2a = 0 ; wo2a = 0.0_rl
154	!
155	!> Read coordinates of all models
156	!
157	IF (l_dryrun .AND. lev_dry >= 1 ) THEN
158	CALL lecdata (lecweights = .TRUE. , leco2amask =.FALSE.)
159	ELSE
160	CALL lecdata (lecweights = .FALSE., leco2amask =.FALSE.)
161	ENDIF
162	!
163	!> Set periodicity
164	CALL lbc (xolont, noperio, 'T', jpoi, jpoj)
165	CALL lbc (xolonu, noperio, 'U', jpoi, jpoj)
166	CALL lbc (xolonv, noperio, 'V', jpoi, jpoj)
167	CALL lbc (xolonf, noperio, 'F', jpoi, jpoj)
168	CALL lbc (xolatt, noperio, 'T', jpoi, jpoj)
169	CALL lbc (xolatu, noperio, 'U', jpoi, jpoj)
170	CALL lbc (xolatv, noperio, 'V', jpoi, jpoj)
171	CALL lbc (xolatf, noperio, 'F', jpoi, jpoj)
172	CALL lbc (xosrft, noperio, 'T', jpoi, jpoj)
173	CALL lbc (iomskt, noperio, 'T', jpoi, jpoj)
174	!
175	!-$$ CALL lbc (xalatt, naperio, 'T', jpai, jpaj)
176	!-$$ CALL lbc (xalatu, naperio, 'U', jpai, jpaj)
177	!-$$ CALL lbc (xalatv, naperio, 'V', jpai, jpaj)
178	!-$$ CALL lbc (xasrft, naperio, 'T', jpai, jpaj)
179	!-$$ CALL lbc (iamskt, naperio, 'T', jpai, jpaj)
180	!>
181	CALL rectifie
182	!>
183	CALL bord_a
184	!
185	CALL bord_o
186	!
187	CALL bord_oa
188	!
189	!> Surfaces par les polygones
190	ALLOCATE (xa (jpae), STAT=ierr) ; CALL chk_allo (ierr, 'xa in mosaic')
191	ALLOCATE (ya (jpae), STAT=ierr) ; CALL chk_allo (ierr, 'ya in mosaic')
192	!
193	IF (.NOT. l_dryrun .AND. l_recalc_a) THEN
194	WRITE (nout,*) 'Recalcul les surfaces atm'
195	DO ja = 1, jpan
196	DO kn = 1, jpae
197	CALL proj (xa_ed (ja,kn), ya_ed (ja,kn), REAL(xa_ed(ja,jpae),rd), REAL(ya_ed(ja,jpae),rd), xa(kn), ya(kn), 2, ierr)
198	ENDDO
199	CALL polom (xa, ya, xa, ya, zs1, zs2, zs)
200	IF (MAX (zs1, zs2) .GT. MIN (zs1, zs2) + zeps .OR. ABS(zs-1.0_rl) > zeps) THEN
201	WRITE (nout,*) 'pb surface atm : ', ja, mai(ja), maj(ja), zs1, zs2, zs
202	ELSE
203	xasrft_pol (ja) = zs1
204	END IF
205	END DO
206	END IF
207	!
208	WRITE (nout,*) 'Perio sur atm'
209	CALL lbc (xasrft_pol, naperio, 'T', jpai, jpaj)
210	!
211	WRITE (nout, *) 'Points T Atmosphere'
212	WRITE (nout, *) ' Longitude Atm '
213	WRITE (nout, '(10F7.1)') (xalont(m1a(jai,jpaj/2)),jai=1,jpai) ! xalont (jpan/2+1:jpan/2+jpai)
214	WRITE (nout, *) ' Latitude Atm '
215	WRITE (nout, '(10F7.1)') (xalatt(m1a(jpai/2,jaj)),jaj=1,jpaj) ! xalatt (1:jpan:jpai)
216	WRITE (nout, *) 'surface Atm'
217	CALL prihre(RESHAPE(xasrft,(/jpai,jpaj/)), 0.0_rl, nout)!-$$ WRITE (nout, *) 'Points east Atm'
218
219	WRITE (nout, *) 'Mask Atm'
220	DO jaj = jpaj, 1, -1
221	WRITE (nout, '(144I1)' ) (iamskt (m1a (jai, jaj)), jai = 1, jpai)
222	END DO
223	WRITE (nout, *) 'Mask Atm perio'
224	DO jaj = jpaj, 1, -1
225	WRITE (nout, '(144I1)' ) (iamskp (m1a (jai, jaj)), jai = 1, jpai)
226	END DO
227	WRITE (nout, *) 'Mask Oce'
228	DO joj = jpoj, 1, -1
229	WRITE (nout, '(1I5, " ", 182I1)') joj, (iomskt (m1o (joi, joj)), joi = 1, jpoi )
230	END DO
231	WRITE (nout, '(1I5, 18I10)') joj, ( joi , joi = 1, jpoi, 10)
232	WRITE (nout, *) 'Mask Oce Perio'
233	DO joj = jpoj, 1, -1
234	WRITE (nout, '(1I5, " ", 182I1)') joj, (iomskp (m1o (joi, joj)), joi = 1, jpoi )
235	END DO
236	WRITE (nout, '(1I5, 18I10)') joj, ( joi , joi = 1, jpoi, 10)
237
238	!
239	ALLOCATE (xo (jpoe), STAT=ierr) ; CALL chk_allo (ierr, 'xo in mosaic')
240	ALLOCATE (yo (jpoe), STAT=ierr) ; CALL chk_allo (ierr, 'yo in mosaic')
241	!
242	IF (.NOT. l_dryrun .AND. l_recalc_o) THEN
243	WRITE (nout,*) 'Recalcul des surfaces oce'
244	ReOce: DO jo = 1, jpon
245	IF ( iomskp(jo) == 0 .AND. iomskt(jo) == 0 ) THEN
246	! WRITE (nout,*) jo, xo_ed (jo,jpoe), yo_ed (jo,jpoe)
247	DO kn = 1, jpoe
248	CALL proj (xo_ed (jo,kn), yo_ed (jo,kn), REAL(xo_ed(jo,jpoe),rd), REAL(yo_ed(jo,jpoe),rd), xo(kn), yo(kn), 2, ierr)
249	ENDDO
250	CALL polom (xo, yo, xo, yo, zs1, zs2, zs)
251	IF (MAX (zs1, zs2) .GT. MIN (zs1, zs2) + zeps ) THEN
252	WRITE (nout,*) 'pb 1 surface oce : ', jo, iomskp(jo), iomskt(jo), m2oi(jo), m2oj(jo), &
253	& xolont(jo), xolatt(jo), zs1, zs2, ABS(zs1-zs2), ABS(zs1-zs2)/MAX(zs1,zs2),zs, zs-1.0_rl
254	END IF
255	IF (ABS(zs-1.0_rl) > zeps) THEN
256	WRITE (nout,*) 'pb 2 surface oce : ', jo, iomskp(jo), iomskt(jo), m2oi(jo), m2oj(jo), &
257	& xolont(jo), xolatt(jo), zs1, zs2, ABS(zs1-zs2), ABS(zs1-zs2)/MAX(zs1,zs2),zs, zs-1.0_rl
258	END IF
259	xosrft_pol (jo) = zs1
260	END IF
261	END DO ReOce
262	END IF
263	!-$$ xosrft_pol = xosrft
264
265	CALL lbc (xosrft_pol, noperio, 'T', jpoi, jpoj)
266	!
267	WRITE (nout, *) 'Points T Ocean'
268	WRITE (nout, *) ' Longitudes Oce '
269	WRITE (nout, '(10F7.1)') xo_ed (1: jpoi, jpoe)
270	DO jno = 1, jpoe-1
271	WRITE (nout, *) 'Longitudes points Oce ', jno
272	WRITE (nout, '(10F7.1)') xo_ed (1: jpoi, jno)
273	WRITE (nout, *) 'Latitudes points Oce ', jno
274	WRITE (nout, '(10F7.1)') yo_ed (1: jpon: jpoj, jno)
275	ENDDO
276	!
277	DO jna = 1, jpae-1
278	WRITE (nout, *) 'Longitudes points Atm ', jna
279	WRITE (nout, '(10F7.1)') xa_ed (1: jpai, jna)
280	WRITE (nout, *) 'Latitudes points Atm ', jna
281	WRITE (nout, '(10F7.1)') ya_ed (1: jpan: jpaj, jna)
282	ENDDO
283	!
284	WRITE (nout,) 'total surface oce ', SUM(xosrftREAL(1-iomskp,rl))
285	WRITE (nout,) 'total surface atm ', SUM(xasrftREAL(1-iamskp,rl))
286	!
287	IF (l_dryrun) THEN
288	SELECT CASE (lev_dry)
289	CASE (0)
290	WRITE (nout,*) "Remplissage bidon des champs pour tester les I/O"
291	!-$$ jma2o = jpa2o ; jmo2a = jpo2a
292	nvo = jma2o ; nva = jmo2a
293	CALL RANDOM_SEED
294	CALL RANDOM_NUMBER (wa2o)
295	CALL RANDOM_NUMBER (wo2a)
296
297	WRITE (nout,*) 'wa2o : ', MINVAL (wa2o), MAXVAL (wa2o)
298	WRITE (nout,*) 'wo2a : ', MINVAL (wo2a), MAXVAL (wo2a)
299
300	ka2o (:,:) = NINT ( wa2o(:,:) * REAL(jpan,rl) )
301	ko2a (:,:) = NINT ( wo2a(:,:) * REAL(jpon,rl) )
302
303	WRITE (nout,*) 'ka2o : ', MINVAL (ka2o), MAXVAL (ka2o)
304	WRITE (nout,*) 'ko2a : ', MINVAL (ko2a), MAXVAL (ko2a)
305
306	wa2o(:,:) = ( 0.5 + wa2o(:,:) ) / REAL ( jma2o, KIND=rl)
307	wo2a(:,:) = ( 0.5 + wo2a(:,:) ) / REAL ( jmo2a, KIND=rl)
308
309	WRITE (nout,*) 'wa2o : ', MINVAL (wa2o), MAXVAL (wa2o)
310	WRITE (nout,*) 'wo2a : ', MINVAL (wo2a), MAXVAL (wo2a)
311
312	!! Masquage
313	DO ja = 1, jpan
314	IF ( iamskt (ja) == 1 ) THEN
315	wo2a (:, ja) = 0.0_rl ; ko2a (:, ja) = 0_il
316	END IF
317	END DO
318	DO jo = 1, jpon
319	IF ( iomskt (jo) == 1 ) THEN
320	wa2o (:, jo) = 0.0_rl ; ka2o (:, jo) = 0_il
321	END IF
322	END DO
323
324	DO ja = 1, jpan
325	DO jna = 1, jpo2a
326	jo = ko2a (jna, ja)
327	IF (jo > 0) THEN
328	IF (iomskt (jo) == 1) THEN
329	wo2a (jna, ja) = 0.0_rl ; ko2a (jna, ja) = 0_il
330	END IF
331	END IF
332	END DO
333	END DO
334
335	DO jo = 1, jpon
336	DO jno = 1, jpa2o
337	ja = ka2o (jno, jo)
338	IF (ja > 0) THEN
339	IF (iamskt (ja) == 1) THEN
340	wa2o (jno, jo) = 0.0_rl ; ka2o (jno, jo) = 0_il
341	END IF
342	END IF
343	END DO
344	END DO
345	!! Normalisation
346	DO jo = 1, jpon
347	wa2o(:,jo) = wa2o (:,jo) / SUM ( wa2o (:,jo))
348	END DO
349	DO ja = 1, jpan
350	wo2a(:,ja) = wo2a (:,ja) / SUM ( wo2a (:,ja))
351	END DO
352
353	CALL gather ('o2a')
354	!CALL inter_o2a (REAL (1-iomskt, KIND = rl), o2amask)
355	!CALL rem_small ('o2a')
356	CALL gather ('a2o')
357
358	WRITE (nout,*) 'ka2o : ', MINVAL (ka2o), MAXVAL (ka2o)
359	WRITE (nout,*) 'ko2a : ', MINVAL (ko2a), MAXVAL (ko2a)
360
361	WRITE (nout,*) 'wa2o : ', MINVAL (wa2o), MAXVAL (wa2o)
362	WRITE (nout,*) 'wo2a : ', MINVAL (wo2a), MAXVAL (wo2a)
363
364	CASE (1)
365	WRITE (nout,*) "Utilisation des poids lus"
366	END SELECT
367	ELSE
368	WRITE (nout,*) 'Appel de polyw_a'
369	CALL polyw_a (lo2a = .TRUE., la2o = .TRUE.)
370	WRITE (nout,*) 'Appel de gather(o2a) apres polyw_a'
371	CALL gather ('o2a')
372	WRITE (nout,*) 'Appel de gather(a2o) apres polyw_a'
373	CALL gather ('a2o')
374	ENDIF
375	!
376
377	!> Checking weights
378	WRITE (unit = nout, fmt = '("Neighbors o2a : ", 3I12 )') MINVAL (nva, nva /=0), MAXVAL (nva), jmo2a
379	WRITE (unit = nout, fmt = '("Index o2a : ", 2I12 )') MINVAL (ko2a, ko2a /= 0), MAXVAL (ko2a)
380	jind = MINLOC(wo2a, wo2a /= 0.0_rl) ; ja = jind(2) ; jna = jind(1); jai = m2ai(ja) ; jaj = m2aj(ja)
381	WRITE (unit = nout, fmt = '("Weights o2a mini : ", 2I8, 2I4, 1E13.4, 3I8)') &
382	& ja, jna, jai, jaj, wo2a(jna,ja), ko2a(jna,ja), m2oi(ko2a(jna,ja)), m2oj(ko2a(jna,ja))
383	jind = MAXLOC(wo2a, wo2a /= 0.0_rl) ; ja = jind(2) ; jna = jind(1); jai = m2ai(ja) ; jaj = m2aj(ja)
384	WRITE (unit = nout, fmt = '("Weights o2a maxi : ", 2I8, 2I4, 1E13.4, 3I8)') &
385	& ja, jna, jai, jaj, wo2a(jna,ja), ko2a(jna,ja), m2oi(ko2a(jna,ja)), m2oj(ko2a(jna,ja))
386	!> Ocean mask interpolated to atm
387	CALL inter_o2a (REAL (1-iomskt, KIND = rl), o2amask)
388	!-$$ WRITE (nout,*) 'Ocean mask interpolated to atm before normalisation'
389	!-$$ CALL prihre(RESHAPE(o2amask,(/jpai,jpaj/)),10.0_rl,INT(nout,kind=i_4))
390	!
391	!> Classes
392	WRITE (nout, *) 'Low sea fraction before normalization'
393	DO kn = 1, 16
394	zzz = 10.0_rl**(-kn)
395	WRITE (nout, *) 'Classe : ', kn, COUNT (o2amask .LE. zzz .AND. o2amask .GT. zzz/10.0_rl)
396	END DO
397	WRITE (nout, *) 'High sea fraction'
398	DO kn = 1, 16
399	zzz = 10.0_rl**(-kn)
400	WRITE (nout, *) 'Classe : ', kn, COUNT (1.0_rl-o2amask .LE. zzz .AND. 1.0_rl-o2amask .GT. zzz/10.0_rl)
401	END DO
402	!> Remove small fractions
403	CALL rem_small ('o2a')
404	!> Recalcule le masque. A faire avant la normalisation
405	CALL inter_o2a (REAL (1-iomskt, KIND = rl), o2amask)
406	! Petit bug ...
407	!IF ( catyp == 'lmdz' ) o2amask (jpan-jpai+1:jpan) = 1.0_rl
408	!
409	!-$$ WRITE (nout, *) 'ATM : Sea fraction < 0'
410	!-$$ DO ja = 1, jpan
411	!-$$ IF (o2amask (ja) < 0.0_rl) THEN
412	!-$$ WRITE (nout,*) ja, m2ai (ja), m2aj (ja), o2amask (ja)
413	!-$$ ENDIF
414	!-$$ ENDDO
415	!-$$ WRITE (nout, *) 'ATM : Sea fraction > 1'
416	!-$$ DO ja = 1, jpan
417	!-$$ IF (o2amask (ja) > 1.0_rl+eps) THEN
418	!-$$ WRITE (nout,*) ja, m2ai (ja), m2aj (ja), o2amask (ja)
419	!-$$ ENDIF
420	!-$$ ENDDO
421
422	o2amask_i_int (:) = 1
423	o2amask_i_ext (:) = 1
424	WHERE ( o2amask (:) > 0.0_rl ) o2amask_i_int = 0
425	WHERE ( o2amask (:) > 1.0_rl-epsfrac ) o2amask_i_ext = 0
426
427	!
428	WRITE (nout,*) 'Appel de gather a2o avant verif'
429	CALL gather ('a2o')
430	!
431	!> Checking weights
432	WRITE (unit = nout, fmt = '("Neighbors a2o : ", 3I9 )') MINVAL (nvo, nvo /=0), MAXVAL (nvo), jma2o
433	WRITE (unit = nout, fmt = '("Index a2o : ", 2I9 )') MINVAL (ka2o, ka2o /= 0), MAXVAL (ka2o)
434	jind = MINLOC(wa2o, wa2o /= 0.0_rl) ; jo = jind(2) ; jno = jind(1); joi = m2oi(jo) ; joj = m2oj(jo)
435	WRITE (unit = nout, fmt = '("Weights a2o mini : ", 2I8, 2I4, 1E13.4, 3I6)') &
436	& jo, jno, joi, joj, wa2o(jno,jo), ka2o(jno,jo), m2ai(ka2o(jno,jo)), m2aj(ka2o(jno,jo))
437	jind = MAXLOC(wa2o, wa2o /= 0.0_rl) ; jo = jind(2) ; jno = jind(1); joi = m2oi(jo) ; joj = m2oj(jo)
438	WRITE (unit = nout, fmt = '("Weights a2o maxi : ", 2I8, 2I4, 1E13.4, 3I6)') &
439	& jo, jno, joi, joj, wa2o(jno,jo), ka2o(jno,jo), m2ai(ka2o(jno,jo)), m2aj(ka2o(jno,jo))
440	!> Atm mask interpolated to oce
441	CALL inter_a2o (REAL (1-iamskt, KIND = rl), a2omask)
442	WRITE (nout,*) 'Atm mask interpolated to oce before normalization'
443	!$$$ CALL prihre(RESHAPE(a2omask,(/jpoi,jpoj/)),10.0_rl,INT(nout,kind=i_4))
444	!> Remove small fractions
445	CALL rem_small ('a2o')
446	CALL inter_a2o (REAL (1-iamskt, KIND = rl), a2omask)
447	!
448	!-$$ WRITE (nout, *) 'OCE : Sea fraction < 0'
449	!-$$ DO jo = 1, jpon
450	!-$$ IF (a2omask(jo) < 0.0) THEN
451	!-$$ WRITE (nout,*) jo, m2oi(jo), m2oj(jo), o2amask(jo)
452	!-$$ ENDIF
453	!-$$ ENDDO
454	!-$$ WRITE (nout, *) 'OCE : Sea fraction > 1'
455	!-$$ DO jo = 1, jpon
456	!-$$ IF (a2omask(jo) > 1.0+eps) THEN
457	!-$$ WRITE (nout,*) jo, m2oi(jo), m2oj(jo), a2omask(jo)
458	!-$$ ENDIF
459	!-$$ ENDDO
460	! Normalise
461	IF (TRIM (catyp) == 'lmdz' ) THEN
462	WRITE (unit=nout, FMT=*) 'Calling normalization '
463	CALL norma (knorma2o = norma2o, knormo2a = normo2a)
464	END IF
465	!
466	WRITE (unit=nout,fmt='("Maxval de sum(wo2a)-1 ", 1E15.4)') MAXVAL (SUM (wo2a(:,:),DIM=1)) -1.0_rl
467	!
468	WRITE (unit = nout, fmt = '("Surface atmos : ", 1PE15.4)') &
469	DOT_PRODUCT (xasrft, REAL(1-iamskp,kind=rl))
470	WRITE (unit = nout, fmt = '("Surface ocean : ", 1PE15.4)') &
471	DOT_PRODUCT (xosrft, REAL(1-iomskp,kind=rl))
472	!
473	WRITE (unit = nout, fmt = '("Oce->Atm : ", 1I9, " non-null adresses, ", 1I12, " non-null weights ")') &
474	& COUNT (wo2a > 0.0_rl), COUNT (ko2a > 0)
475	WRITE (unit = nout, fmt = '("Atm->Oce : ", 1I9, " non-null adresses, ", 1I12, " non-null weights ")') &
476	& COUNT (wa2o > 0.0_rl), COUNT (ka2o > 0)
477	!
478	xsurfa = DOT_PRODUCT (xasrft*REAL(1-iamskp,kind=rl), 1.0_rl - o2amask)
479	xsurfo = DOT_PRODUCT (xosrft, REAL((1-iomskt)*(1-iomskp), kind=rl))
480	WRITE (nout, fmt= '("Ocean surface on ocean and atmosphere : ", 2(1PE20.6))') xsurfo, xsurfa
481	!
482	!> Write Weights for OASIS
483	!
484	mo_name = 'mozaic'
485	o2a_name = "o2a" ; a2o_name = "a2o"
486	WRITE (unit = c_r, fmt = '(1i1)') rk_out
487	IF (lwro2a) THEN
488	IF (l_wei_i4) THEN
489	WRITE (unit = c_i, fmt = '(1i1)') i_4
490	cfname4 = TRIM(mo_name) // ".wo2a" // TRIM(c_suffix) // ".i" // TRIM(c_i) // ".r" // TRIM(c_r)
491	WRITE (nout, *) 'Poids i4 o -> a: ', TRIM(cfname4)
492	OPEN (unit = nwei4o2a, file = TRIM(cfname4), form = 'unformatted', status = 'replace', action = 'write')
493	ENDIF
494	!
495	IF (l_wei_i8) THEN
496	WRITE (unit = c_i, fmt = '(1i1)') i_8
497	cfname8 = TRIM(mo_name) // ".wo2a" // TRIM(c_suffix) // ".i" // TRIM(c_i) // ".r" // TRIM(c_r)
498	WRITE (nout, *) 'Poids i8 : ', TRIM(cfname8)
499	OPEN (unit = nwei8o2a, file = TRIM(cfname8), form = 'unformatted', status = 'replace', action = 'write')
500	ENDIF
501	!
502	cldiag_o2a = TRIM(o2a_name) // '.diag'
503	clw_o2a = TRIM(mo_name) // '.wo2a'
504	clw_o2a_mct = "rmp_"//TRIM(comod_t)//"_to_"//TRIM(camod_t)//"_MOSAIC"
505	CALL wri_weights_o2a (cldiag_o2a, clw_o2a, clw_o2a_mct, "1", l_fulldiag=.FALSE., &
506	& lo_src_grid_frac=.FALSE., lo_dst_grid_frac=.FALSE.)
507	CLOSE (unit = nwei8o2a) ; CLOSE (unit = nwei4o2a)
508
509	ENDIF
510
511	IF (lwra2o) THEN
512	!$$$ WRITE (nout,*) 'Ocean mask interpolated to atm, before writing'
513	!$$$ CALL prihre(RESHAPE(o2amask,(/jpai,jpaj/)),10.0_rl,INT(nout,kind=i_4))
514	IF (l_wei_i4) THEN
515	WRITE (unit = c_i, fmt = '(1i1)') i_4
516	cfname4 = TRIM(mo_name) // ".wa2o" // TRIM(c_suffix) // ".i" // TRIM(c_i) // ".r" // TRIM(c_r)
517	WRITE (nout, *) 'Poids i4 a -> o: ', TRIM(cfname4)
518	OPEN (unit = nwei4a2o, file = TRIM(cfname4), form = 'unformatted', status = 'replace', action = 'write')
519	END IF
520
521	IF (l_wei_i8) THEN
522	WRITE (unit = c_i, fmt = '(1i1)') i_8
523	cfname8 = TRIM(mo_name) // ".wa2o" // TRIM(c_suffix) // ".i" // TRIM(c_i) // ".r" // TRIM(c_r)
524	WRITE (nout, *) 'Poids i8 : ', TRIM(cfname8)
525	OPEN (unit = nwei8a2o, file = TRIM(cfname8), form = 'unformatted', status = 'replace', action = 'write')
526	ENDIF
527
528	cldiag_a2o = TRIM(a2o_name) // '.diag'
529	clw_a2o = TRIM(mo_name) // '.wa2o'
530	clw_a2o_mct = "rmp_"//TRIM(camod_t)//"_to_"//TRIM(comod_t)//"_MOSAIC"
531	CALL wri_weights_a2o (cldiag_a2o, clw_a2o, clw_a2o_mct, "2", l_fulldiag=.FALSE., &
532	& lo_src_grid_frac=.FALSE., lo_dst_grid_frac=.FALSE.)
533
534	IF (l_wei_i4) CLOSE (unit = nwei8a2o)
535	IF (l_wei_i8) CLOSE (unit = nwei4a2o)
536
537	ENDIF
538
539	!$$$ IF (lwro2a) THEN
540	!$$$ !> Compute o2a weights for interpolating fluxes
541	!$$$ !
542	!$$$ !
543	!$$$ !> Normalize weights
544	!$$$ WRITE (nout, *) 'Calling normalization '
545	!$$$ CALL norma (knorma2o = 0_il, knormo2a = 2_il)
546	!$$$ !
547	!$$$ WRITE (unit = cfname4, fmt = '("mozaic.wo2a.flx.i", 1I1, ".r", 1I1)') i_4, rk_out
548	!$$$ OPEN (unit = nwei4o2a, file = TRIM(cfname4), form = 'unformatted', status = 'replace', action = 'write')
549	!$$$ WRITE (nout, *) 'Poids i4 o -> a: ', TRIM(cfname4)
550	!$$$
551	!$$$ WRITE (unit = cfname8, fmt = '("mozaic.wo2a.flx.i", 1I1, ".r", 1I1)') i_8, rk_out
552	!$$$ OPEN (unit = nwei8o2a, file = TRIM(cfname8), form = 'unformatted', status = 'replace', action = 'write')
553	!$$$ WRITE (nout, *) 'Poids i8 : ', TRIM(cfname8)
554	!$$$
555	!$$$ cldiag_o2a = 'o2a_flx_diag' ; clw_o2a = 'mozaic.wo2a.flx'
556	!$$$
557	!$$$ nwei4o2a = 56 ; nwei8o2a = 57
558	!$$$ CALL wri_weights_o2a (cldiag_o2a, clw_o2a, "6")
559	!$$$
560	!$$$ CLOSE (unit = nwei8o2a)
561	!$$$ CLOSE (unit = nwei4o2a)
562	!$$$
563	!$$$ !
564	!$$$ ENDIF
565	!
566	!
567	CALL flioclo
568	WRITE (unit = nout, fmt = *) 'That''s all folks ! '
569	!
570	STOP
571	CONTAINS
572	!
573	SUBROUTINE perio_w (clname)
574	CHARACTER (len=*), INTENT (in) :: clname
575	!
576	SELECT CASE (TRIM(clname))
577	CASE ('a2o')
578	!> Periodicity ocean
579	jma2o = MAXVAL (nvo)
580	DO jno = 1, jma2o
581	CALL lbc (wa2o (jno,:) , noperio, 'T', jpoi, jpoj)
582	CALL lbc (ka2o (jno,:) , noperio, 'T', jpoi, jpoj)
583	END DO
584	CASE ('o2a')
585	!> Periodicity atm
586	jmo2a = MAXVAL (nva)
587	DO jna = 1, jmo2a
588	CALL lbc (wo2a (jna,:) , naperio, 'T', jpai, jpaj)
589	CALL lbc (ko2a (jna,:) , naperio, 'T', jpai, jpaj)
590	END DO
591	!
592	END SELECT
593	END SUBROUTINE perio_w
594	!
595	SUBROUTINE gather (clname)
596	!
597	CHARACTER (len=*), INTENT (in) :: clname
598	INTEGER (kind=il), DIMENSION (:), ALLOCATABLE :: ktemp
599	REAL (kind=rl), DIMENSION (:), ALLOCATABLE :: wtemp
600	INTEGER (kind=il) :: ierr, jn, jon1, jon2, ja1, ja2
601	!
602
603	!
604	SELECT CASE (TRIM(clname))
605	CASE ('o2a')
606	!
607	CALL perio_w (clname)
608	!
609	ALLOCATE (ktemp (jpo2a), STAT=ierr) ; CALL chk_allo (ierr, 'ktemp (jpo2a)', crout='mosaic, gather o2a')
610	ALLOCATE (wtemp (jpo2a), STAT=ierr) ; CALL chk_allo (ierr, 'wtemp (jpo2a)')
611	DO ja = 1, jpan
612	ktemp (:) = ko2a (:,ja) ; wtemp (:) = wo2a (:,ja)
613	ko2a (:,ja) = 0_il ; wo2a (:,ja) = 0.0_rl
614	kn = 1
615	DO jna = 1, jmo2a
616	IF ( ktemp (jna) /= 0_il ) THEN
617	ko2a (kn, ja) = ktemp (jna)
618	wo2a (kn, ja) = wtemp (jna)
619	kn = kn + 1
620	END IF
621	END DO
622	END DO
623	nva = COUNT (ko2a /= 0, DIM = 1)
624	IF ( MAXVAL (nva) > jmo2a ) THEN
625	WRITE (nout, FMT=*) 'Attention: augmentation de jmo2a dans gather o2a'
626	STOP
627	ELSE
628	jmo2a = MAXVAL (nva)
629	END IF
630	WRITE (unit = nout, fmt = *) 'Nombre de voisins o->a : ', jmo2a
631	DEALLOCATE (ktemp) ; DEALLOCATE(wtemp)
632	!
633	CASE ('a2o')
634	IF ( naperio == -1 ) THEN
635	! Pointage vers le point pole
636	DO jo = 1, jpon
637	DO jn = 1, jma2o
638	IF ( m2aj (ka2o (jn, jo)) == 1 ) ka2o (jn, jo) = m1a ( 1, 1 )
639	IF ( m2aj (ka2o (jn, jo)) == jpaj ) ka2o (jn, jo) = m1a ( 1, jpaj)
640	END DO
641	END DO
642	END IF
643	!
644	! Si deux adresses identiques, simplification et somme des poids
645	DO jo = 1, jpon
646	DO jon1 = 1, jma2o
647	ja1 = ka2o (jon1, jo)
648	IF ( ja1 == 0 ) CYCLE
649	DO jon2 = 1, jpa2o
650	IF ( jon1 == jon2 ) CYCLE
651	ja2 = ka2o (jon2, jo)
652	IF ( ja1 == ja2 ) THEN
653	wa2o (jon1, jo) = wa2o (jon1, jo) + wa2o (jon2, jo)
654	wa2o (jon2, jo) = 0.0_rl
655	ka2o (jon2, jo) = 0
656	END IF
657	END DO
658	END DO
659	END DO
660	!
661	CALL perio_w (clname)
662	!
663	ALLOCATE (ktemp (jpa2o), STAT=ierr) ; CALL chk_allo (ierr, 'ktemp (jpa2o)', crout='mosaic, gather a2o')
664	ALLOCATE (wtemp (jpa2o), STAT=ierr) ; CALL chk_allo (ierr, 'wtemp (jpa2o)')
665	DO jo = 1, jpon
666	ktemp = ka2o(:,jo) ; wtemp = wa2o(:,jo)
667	ka2o (:,jo) = 0_il ; wa2o (:,jo) = 0.0_rl
668	kn = 1_il
669	DO jno = 1, jma2o
670	IF ( ktemp (jno) /= 0_il ) THEN
671	ka2o (kn, jo) = ktemp (jno)
672	wa2o (kn, jo) = wtemp (jno)
673	kn = kn + 1_il
674	END IF
675	END DO
676	END DO
677	nvo = COUNT (ka2o /= 0, DIM = 1)
678	IF ( MAXVAL (nvo) > jma2o) THEN
679	WRITE (nout,*) 'Attention: augmentation de nvo dans gather'
680	STOP
681	ELSE
682	jma2o = MAXVAL (nvo)
683	ENDIF
684	WRITE (unit = nout, fmt = *) 'Nombre de voisins a->o : ', jma2o
685	DEALLOCATE (ktemp) ; DEALLOCATE (wtemp)
686	!
687	END SELECT
688
689	!
690	END SUBROUTINE gather
691	!
692	!-$$ SUBROUTINE pole_a2o
693	!-$$ !
694	!-$$ IMPLICIT NONE
695	!-$$ !
696	!-$$ DO jn = 1, jma2o
697	!-$$ DO jo =
698	!-$$
699	!
700	SUBROUTINE rem_small (clname)
701	!
702	CHARACTER (len=*), INTENT (in) :: clname
703	REAL (kind=rl) :: zzs
704	INTEGER :: nn_count
705	!
706	nn_count = 0
707	SELECT CASE (TRIM(clname))
708	CASE ( 'a2o' )
709	!> a -> o
710	DO jo = 1, jpon
711	zzs = SUM (wa2o (:,jo))
712	DO jno = 1, jma2o
713	ja = ka2o (jno, jo)
714	IF (ja > 0) THEN
715	IF (wa2o (jno, jo) < epsfrac) THEN
716	wa2o (jno, jo) = 0.0_rl
717	ka2o (jno, jo) = 0_il
718	nn_count = nn_count + 1
719	ENDIF
720	ENDIF
721	END DO
722	END DO
723	WRITE (nout,*) 'rem small a2o : ', nn_count
724	!
725	CASE ( 'o2a' )
726	DO ja = 1, jpan
727	zzs = SUM (wo2a (:, ja))
728	DO jna = 1, jmo2a
729	jo = ko2a (jna, ja)
730	IF (jo > 0) THEN
731	IF (wo2a (jna, ja) < epsfrac) THEN
732	wo2a (jna, ja) = 0.0_rl
733	ko2a (jna, ja) = 0_il
734	nn_count = nn_count + 1
735	END IF
736	END IF
737	END DO
738	END DO
739	WRITE (nout,*) 'rem small o2a : ', nn_count
740	!
741	END SELECT
742	!
743	WRITE (nout,*) 'Appel de gather (', TRIM(clname),') dans rem_small'
744	CALL gather (clname)
745	!
746	END SUBROUTINE rem_small
747	!
748	END PROGRAM mosaic
749	!

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

Download in other formats: